Анализ землетрясения Лома-Приета
В этом примере показано, как анализировать и визуализировать данные о землетрясении.
Загрузка данных землетрясения
Файл quake.mat содержит данные на 200 Гц из землетрясения Лома-Приета 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. Акселерометр записал основную амплитуду шока волны землетрясения. Переменные nEV обратитесь к трем направленным компонентам, измеренным инструментом, который был выровнен параллельный отказу с его направлением N, указывающим в направлении Сакраменто. Данные не корректируются для ответа инструмента.
Создайте переменную, Time, содержа метки времени, произведенные на уровне 200 Гц с той же длиной как другие векторы. Представляйте правильные модули seconds функция и умножение, чтобы достигнуть Гц () частота дискретизации. Это приводит к a duration переменная, которая полезна для представления прошедшего времени.
Time = (1/200)*seconds(1:length(e))';
whos TimeName Size Bytes Class Attributes Time 10001x1 80010 duration
Организация данных в расписании
Отдельные переменные могут быть организованы в table или timetable для большего удобства. 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 с точечным индексированием. (Для получения дополнительной информации о точечном индексировании см. Доступ к данным в Таблице.) Мы выбрали амплитуду "East-West" и plot это как функция длительности.
plot(quakeData.Time,quakeData.EastWest)
title('East-West Acceleration')
Масштабирование данных
Масштабируйте данные гравитационным ускорением или умножьте каждую переменную в таблице константой. Поскольку переменные все имеют тот же тип (double), можно получить доступ ко всем переменным с помощью имени размерности, Variables. Обратите внимание на то, что quakeData.Variables обеспечивает прямой способ изменить численные значения в рамках расписания.
quakeData.Variables = 0.098*quakeData.Variables;
Выбор Subset of Data for Exploration
Мы интересуемся областью времени, где амплитуда ударной взрывной волны начинает увеличиваться от почти нуля до максимальных уровней. Визуальный осмотр вышеупомянутого графика показывает, что временной интервал с 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.This полезен для применения функций к каждой переменной в таблице или расписании. Функция, чтобы применяться передается 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
Смотрите также
timetable | head | summary | varfun | duration | seconds | timerange