Этот пример показывает, как анализировать и визуализировать данные о землетрясении.
Файл 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. Акселерометр записал основную амплитуду шока волны землетрясения. Переменные n, e, v относится к трем направленным компонентам, измеренным инструментом, который был выровнен параллельный отказу с его направлением N, указывающим в направлении Сакраменто. Данные не исправляются для ответа инструмента.
Создайте переменную, Time, содержа метки времени, выбранные на уровне 200 Гц с той же длиной как другие векторы. Представляйте правильные модули с функцией seconds и умножением, чтобы достигнуть Гц () выборка уровня. Это приводит к переменной 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×4 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;
Мы интересуемся областью времени, где амплитуда ударной взрывной волны начинает увеличиваться от почти нуля до максимальных уровней. Визуальный осмотр вышеупомянутого графика показывает, что временной интервал с 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×4 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×4 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×4 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