Загрузка данных о землетрясениях

Файл 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 функция и умножение для достижения Гц (${\mathit{s}}^{-1}$) частота дискретизации. Это приводит к duration переменная, которая полезна для представления истекшего времени.

Time = (1/200)*seconds(1:length(e))';
whos Time

  Name          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