3.11地震基于主成分分析法的震中方位角估算

基于主成分分析法的震中方位角估算分析马亮卢建旗李山有马强(中国地震局工程力学研究所，中国地震局强震动观测研究室，哈尔滨150080)摘要：地震定位的实现可以有多种方法，主要分为单台定位法、双台定位法、三台定位法和台网定位。单台定位法用于快速定位，优点是快速及时，能实现地震预警；缺点是精度小误差较大。本文使用主成分分析法对日本311地震多个台站测得的数据计算得出的方位角与实际方位角做比较误差分析，希望对实现快速定位有所裨益。关键词:地震预警，震中方位角，主成分分析法，地震定位中图分类号：P315.61文献标识码：A文章编号:引言:1988年,日本铁路技术研究所学者Nakamura在第九届世界地震工程大会上提出了UrEDAS单台地震预警系统，这个系统最初是为了让新干线列车遭遇地震时紧急制动不至于出轨。其中包含了地震的震中方位角确定方法，其原理是使用单个台站接收的最先到达的P波“EW分量到NS分量的振幅比值”导出的估计方位角。这个方法认为，在震源释放的能量的驱使下，地面上的质点M累积的运动轨迹和从震源→质点M的向量所在线段概率性重合，而且质点的累积位移向量与震源→质点向量一致，其的缺陷是它们的方向有可能相反，而且累计计算法得到的震中方位角离散型较大。1992年，吴兆熊、刘希强等的论文里提出了与“基于协方差矩阵的主成分分析法估算震中方位角”的方法有关的数学依据，但并未计算。2006年，刘建华、刘福田、婿颐在其论文里提出了“三分量地震记录的偏振分析”，其中的方法可以用于震中方位角估算。2008年，工程力学研究所学者马强在其博士论文里使用以上方法进行了大量的统计计算，显示了偏振分析法的良好价值。05-07年，李伟，李峰，唐子波等人的主持下完成了《兖州矿区数字遥测矿震台网建设及技术研究》的报告，其中应用了协方差主成分分析法。其主要目的是减轻地震造成的矿难引起的人员伤亡以及财产损失，并做了编程。2010年，山东省地震局学者周彦文、刘希强等再次在论文《基于单台P波记录的快速地震定位方法研究》中明确了协方差矩阵主成分分析法，具有很好的价值。学者发布的论文为震中方位角的估算提供了价值可观的论文与资料，在此基础上，本文对基于强震加速度三分量相关系数矩阵的主成分分析法做以下分析。相关系数矩阵的主成分分析：假如台站接收到的数据（可以使加速度、速度、位移）的垂直、南北、东西三分量存储在列向量（一维数组）D(t)、D2(t)、D3(t)中，令矩阵A=[D(t)，D2(t)，D3(t)]为三维数组。假如强震仪监测到的P波到时为Ka，而且在时间窗[Ka，Ka+F]内P波信息没有收到杂波污染，那么我们就可以对三维数组A1=[D(Ka,Ka+F)，D2(Ka,Ka+F)，D3(Ka,Ka+F)]计算而求得震中方位角，常取时间窗长F∈[0，0.6]秒。计算步骤：1.对三维数组A1做标准差标准化处理A2=zscore(A1)。并且令A2=(a2ij)3*3，（i，j=1，2，3）用zscore标准化的目的是：使得平均值为0，标准差为1，这样可以使不同量纲的数据放在同一个矩阵里。2.求数组A2的相关系数矩阵：R3=(rij)3*3rij（i，j=1，2，3）为原变量a2i与a2j的相关系数，其计算公式为：注意到rij=rji，R3是三阶对称方阵。3.计算R3的特征值与特征向量：用雅可比法（Jacobi）求出特征值，并使其按大小顺序排列λ1≥λ2≥λ3，再求出它们对应的特征向量e1、e2、e3，e1=(e11，e12，e13)，e2=(e21，e22，e23)，e3=(e31，e32，e33)，注意e1、e2、e3的排列顺序对应于λ，而不是从大到小。要求‖ei‖=1，即e2i1+e2i2+e2i3=1，按这样要求，特征向量就是唯一的。4.计算震中方位角azi：[2]azi=f{arctan(|e13|/|e12|)}，具体地当e110，e120和e130时,方位角azi=arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=π-arctan(|e13|/|e12|)12211()()()()nkiikjjkijnnkiikjjkkxxxxrxxxx当e110，e120和e130时，方位角azi=π+arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=2π-arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=π+arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=2π-arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=arctan(|e13|/|e12|)当e110，e120和e130时，方位角azi=π-arctan(|e13|/|e12|)震中方位角是经过台站的本初子午线与台站到震中连线之间的夹角，规定顺时针方向为正。算例：本文使用2011年日本东北青森、岩手、宫城、福岛等县的台站记录的3.11地震及其余震的加速度数据。由于主震和余震皆位于海岸线以东，所以实际方位角az＜180°。图1.本文使用数据对应的K-NET强震台站和震中在地图中的位置K-NET台站使用的强震仪频率为100Hz，在1s时窗内有100个采样点。程序采用的窗口长度是震中方位角取得稳定值时的加速度时段。在求震中方位角azi时，若取窗口长度太短，那么P波初至阶段其幅值小故信噪比低，那么程序算出的方位角azi就很不稳定而且会剧烈跳动；若窗口太长，则在P波初至时段过后受到后续反射折射波的干扰，那么程序计算的角度就会偏离真实方向，使误差增大。本文采用的时间窗口长度大比分满足F≤0.6秒，也就是说在地震预警中在线运行时至少有0.6秒的延误。选用的台站在距离震中350km的范围内，最近的台站距离3.11震中约100km，求取的结果如表1所示，误差的绝对值平均为28.8998度。主震的矩震级为9.0，震源深度为24.4km，因此各台站数据的震相都很明显，东京都到震中距离较大，其境内的台站取得的数据的噪声很大，所以没有选取。从主震获取数据的台站来看，东北地方有一些台站没有记录，而在它们在余震中却取得了地震数据。从表1中可以看出，宫城县（台站代码皆为MYG开头）的四个台站取得的数据误差较小，而青森县台站取得的数据误差较大，这说明距离震中越远，台站测得方位角误差越大。这条规律的实际意义告诉我们：要想取得更精确的方位角，台站的密度必须足够大。编码地震编码:台站名+11年03月11日14时46秒台站所在位置窗口长度F(s)真实震中方位角az(°)计算震中方位角azi(°)误差的绝对值s=|azi-az|备注1MYG0041103111446宫城县0.6112.9243112.62410.3002手动捡拾震相2MYG0061103111446宫城县0.6107.2108121.032713.82193MYG0131103111446宫城县0.695.5443105.27759.73324MYG0161103111446宫城县0.386.24291.17124.92925FKS0031103111446福岛县0.378.9539141.204862.25096FKS0161103111446福岛县0.2964.313360.15714.15627FKS0171103111446福岛县0.666.684589.500722.81628FKS0181103111446福岛县0.2269.544883.140313.59559FKS0191103111446福岛县0.0974.613391.362216.748910IWT0101103111446岩手县0.3120.7973144.117123.319811IWT0111103111446岩手县0.09127.4608177.403549.942712IWT0121103111446岩手县0.19131.6444189.35957.714613IWT0141103111446岩手县0.75135.4364201.693166.256714IWT0181103111446岩手县0.07139.5349185.841646.306715IWT0241103111446岩手县0.35180.3406150.107730.232916AOM0111103111446青森县0.3155.404138.2743117.129817AOM0121103111446青森县0.2155.6789197.556641.877718SIT0021103111446崎玉县0.653.722779.461225.738519IBR0091103111446茨城县0.1251.732742.80338.929420IBR0091103111515茨城县0.693.025293.49420.469余震21FKS0161103141002福岛县0.6131.3624128.66932.6931手动捡拾震相22GNM0131103191856群马县0.1869.094198.395329.3012余震23TCG0011103111515栃木県0.6130.916114.484216.4318余震误差的绝对值的平均值av=average(|azi-az|)28.8998表1.计算震中方位角的误差统计误差分析：为了更直观地表示主基于成分分析法的估算方位角离差，我们把表1的结果描述到图2中，散点的标签与台站的编码相对应。散点离对角线越近意味着对应台站的数据估算的方位角越准确。但上述结果显示，误差的绝对值的平均高达28度，这样的误差只能对震中的方位进行粗略估计，要想用于实际预警，还要结合场地的地质特性、断裂带特征等对计算方法加以修改。方位角偏差的产生原因我们现解释如下：[1]一.系统误差，可以通过一定的手段来消除1.仪器影响：包括增益、布设方位的影响。计算震中方位角时需要三分量记录保持一致，但是三分量记录不同频带本身的增益具有差异，而且由于三分量记录事件的频率成分并不完全一致，放大倍数也有差异；由于仪器的坐标和大地的坐标不能完全重合，仪器布设时会有一定的方位偏差。图2.由表1绘制的计算震中方位角误差分析示意图二．随机误差，使得方位角的偏差服从正态分布，如图3所示1.信噪比影响：在P波初至阶段，地震波幅值小，信噪比低，噪声会引起方位角的误差。2.震相捡拾误差的影响：大多数情况下P波捡拾到的震相都有足够的精度以满足预警的要求，个别记录由于软件选取的AIC计算时间窗过长导致触发点提前，反而不如STA/LTA法捡拾的结果准确，所以表1中对个别记录进行了手动捡拾。3.时间窗口影响：每条地震记录取得最佳方位角的时间窗起点和时窗口长度都不一样，我们尽可能地使它们一致为0.6秒。由于我们选用的是固定窗，窗口长度太短就包含不了信噪比较大的有用信息，只包含信噪比低的无价值信息；窗口长度太长则会引入尾波并且使延误时间增长，这样对于地震预警就失去了意义。由于震源和台站之间的介质不均匀，P波在传播过程中穿过不同介质的边界030609012015018021003060901201501802101234567891011121314151617181920212223az()azi()Epicenterazimuth面时，会产生折射、反射、衍射现象，加之介质的密度和弹性的差异，结果后续到达台站的P波(尾波)轨迹没有和震中-台站之间的连线吻合。初到台站的P波是沿着波速最快的介质传播的，它与震中-台站之间的连线吻合最好；4.加速度记录的零线误差的影响：加速度的零线是计算出的总体平均值，每一时段的零线可能不尽相同，尤其是东西、南北两分量的在时间窗内的零线稍有偏差就会导致方位角产生很大的偏差，这需要寻找更好的算法。5.地震方位影响：地震仪在布设时其坐标对应于南北和东西方向，若震中方位角趋近于90°的倍数