Анализ землетрясений в Лома-Приете
Этот пример показывает, как анализировать и визуализировать данные о землетрясении.
Загрузка данных о землетрясениях
Файл 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
В рабочей области есть три переменные, содержащие временные следы от акселерометра, расположенного в создании естественных наук в UC Santa Cruz. Акселерометр зафиксировал основную амплитуду удара волны землетрясения. Переменные n, e, v относятся к трем направленным компонентам, измеренным прибором, который был выровнен параллельно отказу с направлением N в направлении Сакраменто. Данные не скорректированы для отклика инструмента.
Создайте переменную, Time, содержащая временные метки, дискретизированные в 200Hz с той же длиной, что и другие векторы. Представление правильных модулей измерения с seconds функция и умножение для достижения Гц () частота дискретизации. Это приводит к duration переменная, которая полезна для представления истекшего времени.
Time = (1/200)*seconds(1:length(e))';
whos TimeName Size Bytes Class Attributes Time 10001x1 80010 duration
Организуйте данные в Timetable
Отдельные переменные могут быть организованы в 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 с помощью значения компонента. Ниже мы покажем различные способы визуализации этих данных.
Начнем с 2-мерных проекций. Вот первый с несколькими значениями времени, аннотированными.
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