洛马普列塔地震分析
此示例说明如何在时间表中存储带时间戳的地震数据,以及如何使用时间表函数来分析和可视化数据。
加载地震数据
示例文件 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 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)
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)
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)
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 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 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
另请参阅
timetable
| head
| summary
| varfun
| duration
| seconds
| timerange