洛马普列塔地震分析
此示例说明如何在时间表中存储带时间戳的地震数据，以及如何使用时间表函数来分析和可视化数据。
加载地震数据
示例文件 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
工作区中有三个变量，包含由位于加州大学圣克鲁斯分校的自然科学大楼的一个加速度计记录的时间跟踪。加速度计记录了地震波的主震振幅。变量 n、e、v 指代该工具测量的三个定向分量，工具已调整为与断层平行，其 N 方向指向萨克拉门托方向。数据未针对工具响应进行纠正。
创建一个变量 Time，其中包含以 200Hz 采样的时间戳，并且长度与其他向量相同。使用 seconds 函数和乘法表示正确的单位以实现该 Hz () 采样率。这将生成一个适用于表示已用时间的 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)
quakeData=10001×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
0.045 sec 5 1 0
0.05 sec 5 0 0
0.055 sec 5 0 0
0.06 sec 5 0 0
0.065 sec 5 0 0
0.07 sec 5 0 0
0.075 sec 5 0 0
0.08 sec 5 0 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); t2 = seconds(15); xline(t1,LineWidth=2) xline(t2,LineWidth=2)
存储相关数据
使用此区间中的数据创建另一个 timetable。使用 timerange 选择感兴趣的行。
tr = timerange(t1,t2); quakeData8to15 = quakeData(tr,:)
quakeData8to15=1400×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
8.04 sec -6.37 2.646 -6.566
8.045 sec -4.508 2.744 -4.9
8.05 sec -2.646 2.842 -3.234
8.055 sec -0.784 2.842 -1.568
8.06 sec 1.078 2.548 0.098
8.065 sec 2.94 2.254 1.764
8.07 sec 4.802 1.96 3.43
8.075 sec 6.664 1.666 4.998
⋮
在单独的坐标区上以可视方式呈现三个加速度变量。要使用以 3×1 图块排列的三个图创建分块布局，请使用 tiledlayout 函数。
tiledlayout(3,1) % Tile 1 nexttile plot(quakeData8to15.Time,quakeData8to15.EastWest) ylabel("East-West") title("Acceleration") % tile 2 nexttile plot(quakeData8to15.Time,quakeData8to15.NorthSouth) ylabel("North-South") % Tile 3 nexttile plot(quakeData8to15.Time,quakeData8to15.Vertical) ylabel("Vertical")
计算摘要统计量
使用 summary 函数显示有关数据的统计信息。
summary(quakeData8to15)
quakeData8to15: 1400×3 timetable
Row Times:
Time: duration
Variables:
EastWest: double
NorthSouth: double
Vertical: double
Statistics for applicable variables and row times:
NumMissing Min Median Max Mean Std
Time 0 8 sec 11.498 sec 14.995 sec 11.498 sec 2.0214 sec
EastWest 0 -255.0940 -0.0980 244.5100 0.9338 66.1188
NorthSouth 0 -198.5480 1.0780 204.3300 -0.1028 42.7980
Vertical 0 -157.8780 0.9800 134.4560 -0.5254 35.8554
可以使用 varfun 计算其他数据统计信息。在对表或时间表中的每个变量应用函数时，这会非常有用。要应用的函数将以函数句柄的形式传递到 varfun。将 mean 函数应用于全部三个变量并以表格式输出结果，因为在计算了时间均值之后，时间会变得没有意义。
mn = varfun(@mean,quakeData8to15,OutputFormat="table")mn=1×3 table
mean_EastWest mean_NorthSouth mean_Vertical
_____________ _______________ _____________
0.9338 -0.10276 -0.52542
计算速度和位置
为了确定冲击波的传播速度，对加速度进行一次积分。使用沿时间变量的累积和来计算波前速度。
edot = (1/200)*cumsum(quakeData8to15.EastWest); edot = edot - mean(edot);
接下来对全部三个变量执行积分来计算该速度。可以方便地创建函数，并使用 varfun 将该函数应用于 timetable 中的变量。在此示例中，该函数包含在此示例的末尾并命名为 velFun。
vel = varfun(@velFun,quakeData8to15)
vel=1400×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
8.04 sec -0.79273 0.5454 1.5929
8.045 sec -0.81527 0.55912 1.5684
8.05 sec -0.8285 0.57333 1.5522
8.055 sec -0.83242 0.58754 1.5444
8.06 sec -0.82703 0.60028 1.5449
8.065 sec -0.81233 0.61155 1.5537
8.07 sec -0.78832 0.62135 1.5708
8.075 sec -0.755 0.62968 1.5958
⋮
将同一函数 velFun 应用于速度以确定位置。
pos = varfun(@velFun,vel)
pos=1400×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
8.04 sec 2.0923 -2.1592 -3.0128
8.045 sec 2.0883 -2.1564 -3.005
8.05 sec 2.0841 -2.1536 -2.9972
8.055 sec 2.08 -2.1506 -2.9895
8.06 sec 2.0758 -2.1476 -2.9818
8.065 sec 2.0718 -2.1446 -2.974
8.07 sec 2.0678 -2.1415 -2.9662
8.075 sec 2.064 -2.1383 -2.9582
⋮
请注意 varfun 所创建的时间表中的变量名称是如何包含所使用函数的名称的。此方法有助于跟踪已对原始数据执行了哪些运算。使用圆点表示法将变量名称重新调整回其原始值。
pos.Properties.VariableNames = varNames; vel.Properties.VariableNames = varNames;
分别绘制所关注时间段内速度和位置的 3 个分量。在 2×1 图块排列中绘制速度和位置。
tiledlayout(2,1) % Tile 1 nexttile plot(vel.Time,vel.Variables) legend(vel.Properties.VariableNames,Location="bestoutside") title("Velocity") % Tile 2 nexttile plot(pos.Time,pos.Variables) 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 表示网格中的每个坐标区，可用于确定要使用 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
绘制轨迹的三维视图并每隔十个位置点绘制一个竖线。垂直线之间的间距表示速度。
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 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
另请参阅
timetable | head | summary | varfun | duration | seconds | timerange
