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

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

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

  Name          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