Main Content

本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

洛马普列塔地震分析

以下示例演示如何分析和以可视方式呈现地震数据。

加载地震数据

文件 quake.mat 包含圣克鲁斯山脉在 1989 年 10 月 17 日发生的洛马普列塔地震的 200Hz 数据。该数据由加州大学圣克鲁斯分校的 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)
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);

下面我们对全部三个变量执行积分来计算该速度。可以方便地创建函数,并使用 varfun 将该函数应用于 timetable 中的变量。在本例中,我们在文件末尾处包括了该函数,并将其命名为 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')

以可视方式呈现轨迹

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

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

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 表示网格中的每个坐标区,可用于确定要使用 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

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

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

另请参阅

| | | | | |

相关主题