Main Content

洛马普列塔地震分析

此示例说明如何在时间表中存储带时间戳的地震数据,以及如何使用时间表函数来分析和可视化数据。

加载地震数据

示例文件 quake.mat 包含圣克鲁斯山脉在 1989 年 10 月 17 日发生的洛马普列塔地震的 200 Hz 数据。该数据由加州大学圣克鲁斯分校的 Joel Yellin 通过 Charles F. Richter 地震实验室提供。

首先加载数据。

load quake 
whos e n v
  Name          Size            Bytes  Class     Attributes

  e         10001x1             80008  double              
  n         10001x1             80008  double              
  v         10001x1             80008  double              

工作区中有三个变量,包含由位于加州大学圣克鲁斯分校的自然科学大楼的一个加速度计记录的时间跟踪。加速度计记录了地震波的主震振幅。变量 nev 指代该工具测量的三个定向分量,工具已调整为与断层平行,其 N 方向指向萨克拉门托方向。数据未针对工具响应进行纠正。

创建一个变量 Time,其中包含以 200Hz 采样的时间戳,并且长度与其他向量相同。使用 seconds 函数和乘法表示正确的单位以实现该 Hz (s1) 采样率。这将生成一个适用于表示已用时间的 duration 变量。

Time = (1/200)*seconds(1:length(e))';
whos Time
  Name          Size            Bytes  Class       Attributes

  Time      10001x1             80010  duration              

组织时间表中的数据

可以将这些单独的变量放到 tabletimetable 中,以便于使用。timetable 提供了灵活丰富的时间戳数据处理功能。使用时间和三个加速度变量创建 timetable,并提供更有意义的变量名称。使用 head 函数来显示前八行。

varNames = {'EastWest', 'NorthSouth', 'Vertical'};
quakeData = timetable(Time, e, n, v, 'VariableNames', varNames);
head(quakeData)
      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')

Figure contains an axes object. The axes object with title East-West Acceleration contains an object of type line.

缩放数据

按照重力加速度缩放数据或者将表中的每个变量与常量相乘。由于这些变量都具有相同类型 (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

Figure contains an axes object. The axes object with title East-West Acceleration contains 3 objects of type line.

存储相关数据

使用此区间中的数据创建另一个 timetable。使用 timerange 选择感兴趣的行。

tr = timerange(seconds(8),seconds(15));
dataOfInterest = quakeData(tr,:);
head(dataOfInterest)
      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')

Figure contains 3 axes objects. Axes object 1 with title Acceleration, ylabel East-West contains an object of type line. Axes object 2 with ylabel North-South contains an object of type line. Axes object 3 with ylabel Vertical contains an object of type line.

计算摘要统计量

使用 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);

接下来对全部三个变量执行积分来计算该速度。可以方便地创建函数,并使用 varfun 将该函数应用于 timetable 中的变量。在此示例中,该函数包含在此示例的末尾并命名为 velFun

vel = varfun(@velFun,dataOfInterest);
head(vel)
      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)
      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')

Figure contains 2 axes objects. Axes object 1 with title Velocity contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical. Axes object 2 with title Position contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical.

以可视方式呈现轨迹

可使用分量值绘制二维或三维轨迹。此图显示了以可视方式呈现这些数据的不同方式。

首先使用二维投影。下面是第一张标注有几个时间值的图。

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)))

Figure contains an axes object. The axes object with xlabel North-South, ylabel Vertical contains 5 objects of type line, text.

使用 plotmatrix 以可视方式呈现一个网格,其中包含每个变量相对另一变量的散点图,对角线上则是每个变量自身的直方图。输出变量 Ax 表示网格中的每个坐标区,可用于确定要使用 xlabelylabel 标记的坐标区。

figure
[S,Ax] = plotmatrix(pos.Variables);

for ii = 1:length(varNames)
    xlabel(Ax(end,ii),varNames{ii})
    ylabel(Ax(ii,1),varNames{ii})
end

MATLAB figure

绘制轨迹的三维视图并每隔十个位置点绘制一个圆点。点之间的间距即表示速度。

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')

Figure contains an axes object. The axes object with title Position, xlabel North-South, ylabel East-West contains 2 objects of type line. One or more of the lines displays its values using only markers

支持函数

下面定义了这些函数。

function y = velFun(x)
    y = (1/200)*cumsum(x);
    y = y - mean(y);
end

另请参阅

| | | | | |

相关主题