В этом примере показано, как анализировать и визуализировать данные о землетрясениях.
Файл quake.mat содержит 200Hz данные землетрясения Лома-Приета 17 октября 1989 года в горах Санта-Крус. Данные предоставлены Джоэлом Йеллином в Сейсмологической лаборатории Чарльза Ф. Рихтера, Калифорнийский университет, Санта-Крус.
Начните с загрузки данных.
load quake whos e n v
Name Size Bytes Class Attributes e 10001x1 80008 double n 10001x1 80008 double v 10001x1 80008 double
В рабочем пространстве есть три переменные, содержащие временные следы от акселерометра, расположенного в здании естественных наук в Калифорнийском университете в Санта-Крус. Акселерометр зафиксировал основную ударную амплитуду волны землетрясения. Переменные n, e, v относятся к трем направленным компонентам, измеренным прибором, который был выровнен параллельно разлому, при этом его N-направление было направлено в направлении Сакраменто. Данные не корректируются для ответа прибора.
Создайте переменную, Time, содержащий временные метки, дискретизированные в 200Hz с той же длиной, что и другие векторы. Представление правильных единиц измерения с помощью seconds функция и умножение для достижения частоты дискретизации Гц (s − 1). Это приводит к duration переменная, которая полезна для представления прошедшего времени.
Time = (1/200)*seconds(1:length(e))';
whos TimeName Size Bytes Class Attributes Time 10001x1 80010 duration
Отдельные переменные могут быть организованы в table или timetable для большего удобства. A timetable обеспечивает гибкость и функциональность для работы с данными с отметками времени. Создать timetable со временем и тремя переменными ускорения и предоставить более значимые имена переменных. Отображение первых восьми строк с помощью head функция.
varNames = {'EastWest', 'NorthSouth', 'Vertical'};
quakeData = timetable(Time, e, n, v, 'VariableNames', varNames);
head(quakeData)ans=8×3 timetable
Time EastWest NorthSouth Vertical
_________ ________ __________ ________
0.005 sec 5 3 0
0.01 sec 5 3 0
0.015 sec 5 2 0
0.02 sec 5 2 0
0.025 sec 5 2 0
0.03 sec 5 2 0
0.035 sec 5 1 0
0.04 sec 5 1 0
Изучение данных путем доступа к переменным в timetable с точечным подстрочным индексом. (Дополнительные сведения о подстрочном индексировании точек см. в разделе Доступ к данным в таблице.) Мы выбрали амплитуду «Восток-Запад» и plot это как функция длительности.
plot(quakeData.Time,quakeData.EastWest)
title('East-West Acceleration')
Масштабируйте данные на гравитационное ускорение или умножьте каждую переменную в таблице на константу. Поскольку все переменные одного типа (double), вы можете получить доступ ко всем переменным, используя имя измерения, Variables. Обратите внимание, что quakeData.Variables предоставляет прямой способ изменения числовых значений в расписании.
quakeData.Variables = 0.098*quakeData.Variables;
Нас интересует временная область, где амплитуда ударной волны начинает увеличиваться от почти нуля до максимальных уровней. Визуальный осмотр вышеуказанного графика показывает, что представляет интерес интервал времени от 8 до 15 секунд. Для лучшей визуализации мы рисуем черные линии в выбранных временных точках, чтобы привлечь внимание к этому интервалу. Все последующие расчеты будут включать этот интервал.
t1 = seconds(8)*[1;1]; t2 = seconds(15)*[1;1]; hold on plot([t1 t2],ylim,'k','LineWidth',2) hold off
Создать другое timetable с данными в этом интервале. Использовать timerange для выбора интересующих строк.
tr = timerange(seconds(8),seconds(15)); dataOfInterest = quakeData(tr,:); head(dataOfInterest)
ans=8×3 timetable
Time EastWest NorthSouth Vertical
_________ ________ __________ ________
8 sec -0.098 2.254 5.88
8.005 sec 0 2.254 3.332
8.01 sec -2.058 2.352 -0.392
8.015 sec -4.018 2.352 -4.116
8.02 sec -6.076 2.45 -7.742
8.025 sec -8.036 2.548 -11.466
8.03 sec -10.094 2.548 -9.8
8.035 sec -8.232 2.646 -8.134
Визуализируйте три переменные ускорения на трех отдельных осях.
figure subplot(3,1,1) plot(dataOfInterest.Time,dataOfInterest.EastWest) ylabel('East-West') title('Acceleration') subplot(3,1,2) plot(dataOfInterest.Time,dataOfInterest.NorthSouth) ylabel('North-South') subplot(3,1,3) plot(dataOfInterest.Time,dataOfInterest.Vertical) ylabel('Vertical')
Для отображения статистической информации о данных используется summary функция.
summary(dataOfInterest)
RowTimes:
Time: 1400x1 duration
Values:
Min 8 sec
Median 11.498 sec
Max 14.995 sec
TimeStep 0.005 sec
Variables:
EastWest: 1400x1 double
Values:
Min -255.09
Median -0.098
Max 244.51
NorthSouth: 1400x1 double
Values:
Min -198.55
Median 1.078
Max 204.33
Vertical: 1400x1 double
Values:
Min -157.88
Median 0.98
Max 134.46
Дополнительная статистическая информация о данных может быть рассчитана с помощью varfun.Это полезно для применения функций к каждой переменной в таблице или расписании. Применяемая функция передается varfun как дескриптор функции. Ниже мы применяем mean функция для всех трех переменных и вывод результата в формате таблицы, так как время не имеет значения после вычисления временного средства.
mn = varfun(@mean,dataOfInterest,'OutputFormat','table')
mn=1×3 table
mean_EastWest mean_NorthSouth mean_Vertical
_____________ _______________ _____________
0.9338 -0.10276 -0.52542
Чтобы определить скорость распространения ударной волны, мы интегрируем ускорения один раз. Мы используем кумулятивные суммы вдоль переменной времени, чтобы получить скорость волнового фронта.
edot = (1/200)*cumsum(dataOfInterest.EastWest); edot = edot - mean(edot);
Ниже выполняется интегрирование всех трех переменных для вычисления скорости. Удобно создать функцию и применить ее к переменным в timetable с varfun. В этом примере мы включили функцию в конец этого файла и назвали ее velFun.
vel = varfun(@velFun,dataOfInterest); head(vel)
ans=8×3 timetable
Time velFun_EastWest velFun_NorthSouth velFun_Vertical
_________ _______________ _________________ _______________
8 sec -0.56831 0.44642 1.8173
8.005 sec -0.56831 0.45769 1.834
8.01 sec -0.5786 0.46945 1.832
8.015 sec -0.59869 0.48121 1.8114
8.02 sec -0.62907 0.49346 1.7727
8.025 sec -0.66925 0.5062 1.7154
8.03 sec -0.71972 0.51894 1.6664
8.035 sec -0.76088 0.53217 1.6257
Применить ту же функцию velFun к скоростям для определения положения.
pos = varfun(@velFun,vel); head(pos)
ans=8×3 timetable
Time velFun_velFun_EastWest velFun_velFun_NorthSouth velFun_velFun_Vertical
_________ ______________________ ________________________ ______________________
8 sec 2.1189 -2.1793 -3.0821
8.005 sec 2.1161 -2.177 -3.0729
8.01 sec 2.1132 -2.1746 -3.0638
8.015 sec 2.1102 -2.1722 -3.0547
8.02 sec 2.107 -2.1698 -3.0458
8.025 sec 2.1037 -2.1672 -3.0373
8.03 sec 2.1001 -2.1646 -3.0289
8.035 sec 2.0963 -2.162 -3.0208
Обратите внимание, как имена переменных в расписании создаются varfun укажите имя используемой функции. Полезно отслеживать операции, выполненные с исходными данными. Откорректируйте имена переменных до их исходных значений с помощью точечной нотации.
pos.Properties.VariableNames = varNames;
Ниже мы построим график 3 компонентов скорости и положения для интересующего интервала времени.
figure subplot(2,1,1) plot(vel.Time,vel.Variables) legend(quakeData.Properties.VariableNames,'Location','NorthWest') title('Velocity') subplot(2,1,2) plot(vel.Time,pos.Variables) legend(quakeData.Properties.VariableNames,'Location','NorthWest') title('Position')
Траектории можно нарисовать в 2D или 3D с помощью значения компонента. Далее мы покажем различные способы визуализации этих данных.
Начните с двумерных проекций. Вот первый с несколькими значениями времени, аннотированными.
figure plot(pos.NorthSouth,pos.Vertical) xlabel('North-South') ylabel('Vertical') % Select locations and label nt = ceil((max(pos.Time) - min(pos.Time))/6); idx = find(fix(pos.Time/nt) == (pos.Time/nt))'; text(pos.NorthSouth(idx),pos.Vertical(idx),char(pos.Time(idx)))
Использовать plotmatrix для визуализации сетки графиков рассеяния всех переменных относительно друг друга и гистограмм каждой переменной на диагонали. Выходная переменная Ax, представляет каждую ось в сетке и может использоваться для определения осей для маркировки с помощью xlabel и ylabel.
figure [S,Ax] = plotmatrix(pos.Variables); for ii = 1:length(varNames) xlabel(Ax(end,ii),varNames{ii}) ylabel(Ax(ii,1),varNames{ii}) end
Постройте график 3-D вида траектории и постройте график точки в каждой десятой точке положения. Интервал между точками указывает на скорость.
step = 10; figure plot3(pos.NorthSouth,pos.EastWest,pos.Vertical,'r') hold on plot3(pos.NorthSouth(1:step:end),pos.EastWest(1:step:end),pos.Vertical(1:step:end),'.') hold off box on axis tight xlabel('North-South') ylabel('East-West') zlabel('Vertical') title('Position')
Функции определены ниже.
function y = velFun(x) y = (1/200)*cumsum(x); y = y - mean(y); end
duration | head | seconds | summary | timerange | timetable | varfun