在地震事件发生后，快速准确地确定震中位置对于地震应急响应至关重要^[1]。目前，地震学研究主要采用高频全球卫星导航系统(HR-GNSS)和地震传感器^[2]，后者包含强震仪和测震仪，其中强震仪获取的强震动记录因其仪器动态范围大、灵敏度高、分辨率高、数据存储方便等特点，成为地震学领域研究的重要数据来源^[3]。

利用强震动记录计算震中位置的方法中，以传统Geiger法为代表的线性定位方法^[4]及格网搜索^[5]等非线性定位方法一直被广泛应用，然而Geiger方法依赖P波到时信息，可能受P波到时粗差的影响，从而降低震中位置计算的准确性^[6]。抗差估计方法可通过构造等价权函数的形式，基于迭代计算，抑制观测值粗差对参数估计结果的影响^[7]，在GNSS定位^[8]、水下声学定位^[9]等领域应用广泛且效果显著，但在震中定位方面应用较少。本文在Geiger方法的基础上，基于IGGⅢ等价权选择准则及抗差估计原理，提出一种基于选权迭代法的震中估计新方法，通过构造等价权函数和迭代计算，抵抗P波到时信息潜在的粗差影响，实现更为准确可靠的地震震中位置确定。

1 震中估计方法 1.1 P波到时提取方法

长短时平均方法是目前常用的P波到时提取方法之一^[4]。该方法通过计算信号的短时平均值STA和长时平均值LTA及二者之比STA/LTA来判断地震信号中P波震相的到达时刻。特征函数是反映信号到达时频率或幅值特征变化的时间序列函数，能够对异常信号进行放大，在利用STA/LTA方法识别P波震相时，通常需要引入特征函数。本文使用Allen^[10-11]提出的特征函数，该特征函数能够突出垂直向的记录特征：

$ \mathrm{CF}_i=X_i^2-\left(X_i-X_{i-1}\right)^2 $

(1)

式中，i为观测历元数，CF_i为地震信号X_i的特征函数，X_i为i时刻垂直向的加速度记录值。

确定特征函数后，利用其计算长短时信号的比值：

$ \frac{\mathrm{STA}_i}{\mathrm{LTA}_i}=\frac{\sum\limits_{j=i+1-N_S}^i \mathrm{CF}_j / N_S}{\sum\limits_{k=i+1-N_L}^i \mathrm{CF}_k / N_L} $

(2)

式中，N_S为短窗包含的记录点数，N_L为长窗包含的记录点数。当$\frac{\mathrm{STA}_i}{\mathrm{LTA}_i} $超过触发阈值时，说明有异常信号出现，即可判定为发生地震事件。

在P波到时拾取过程中，STA/LTA算法主要用于P波初至时刻的粗略提取，而赤池信息量准则(AIC)方法可用于P波初至时刻的精细提取。AIC准则的基本思想是：选定一定的时间窗口，求解窗口内信号的AIC曲线，曲线极小值点对应的时刻即为背景噪声与地震信号的最佳划分点，即P波震相到时点^[12]。对于地震信号而言，AIC函数的计算公式为：

$ \begin{gathered} \mathrm{AIC}(k)=k \log _{10}\{\operatorname{var}(X[1, k])\}+ \\ (N-k-1) \log _{10}\{\operatorname{var}(X[k+1, N])\} \end{gathered} $

(3)

式中，var表示方差计算函数，k为加速度记录的序列数，N为整个记录窗内包含的加速度记录点数，X[1, k]表示的是[1, k]区间的加速度值组成的数组。

1.2 传统Geiger震中估计方法

传统Geiger方法假定地震波传播速度具有区域一致性，将地震到时函数进行一阶泰勒级数展开，通过迭代计算使各台站到时残差平方和达到最小，进而计算出震中位置^[13]。各台站至震中的空间距离(震中距)计算公式及地震波到时函数模型为：

$ D_i=\sqrt{\left(x_i-x_0\right)^2+\left(y_i-y_0\right)^2+\left(z_i-z_0\right)^2} $

(4)

$ t_i=\frac{D_i}{V_{\mathrm{P}}}+t_0 $

(5)

式中，下标i表示台站个数，D_i为震中距，(x_i, y_i, z_i)为台站的位置，(x₀, y₀, z₀)为地心地固坐标系下的震中初始位置，t_i为各台站触发的P波到时，t₀为发震时刻，V_P为地震波传播速度。

采用台站间差分方法消去发震时刻t₀^[13]：

$ \left\{\begin{array}{c} D_2-D_1-V_{\mathrm{P}}\left(t_2-t_1\right)=0 \\ \vdots \\ D_n-D_1-V_{\mathrm{P}}\left(t_n-t_1\right)=0 \end{array}\right. $

(6)

式中，n为触发的台站总数目。

将式(6)进行线性化，得到对应的误差方程：

$ \boldsymbol{V}=\boldsymbol{A} \mathrm{d} \boldsymbol{x}-\boldsymbol{l} $

(7)

式中，A为设计矩阵，l为观测向量，dx为改正数。具体可分别表示为：

$ \left\{\begin{aligned} \boldsymbol{l}= & {\left[\begin{array}{c} \left(D_2^0-D_1^0\right)-V_{\mathrm{P}, 0}\left(t_2-t_1\right) \\ \vdots \\ \left(D_n^0-D_1^0\right)-V_{\mathrm{P}, 0}\left(t_n-t_1\right) \end{array}\right], \mathrm{d} \boldsymbol{x}=\left[\mathrm{d} X, \mathrm{~d} Y, \mathrm{~d} Z, \mathrm{~d} V_{\mathrm{P}}\right]^{\mathrm{T}} } \\ \boldsymbol{A}= & {\left[\begin{array}{cccc} \left(\frac{x_2-x_0}{D_2^0}-\frac{x_1-x_0}{D_1^0}\right), & \left(\frac{y_2-y_0}{D_2^0}-\frac{y_1-y_0}{D_1^0}\right), & \left(\frac{z_2-z_0}{D_2^0}-\frac{z_1-z_0}{D_1^0}\right), & \left(t_2-t_1\right) \\ \vdots & \vdots & \vdots \\ \left(\frac{x_n-x_0}{D_n^0}-\frac{x_1-x_0}{D_1^0}\right), & \left(\frac{y_n-y_0}{D_n^0}-\frac{y_1-y_0}{D_1^0}\right), & \left(\frac{z_n-z_0}{D_n^0}-\frac{z_1-z_0}{D_1^0}\right), & \left(t_n-t_1\right) \end{array}\right] } \end{aligned}\right. $

(8)

式中，dX、dY、dZ为震中位置的改正量，dV_P为P波速度的改正量。

当触发地震的台站大于等于4个时，计算4个台站的重心位置作为震中初始位置(x₀, y₀, z₀)，综合各台站的坐标及P波震相到时，利用最小二乘原理求解改正数dx，其表达式为：

$ \mathrm{d} \boldsymbol{x}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{l} $

(9)

式中，P为各台站地震波提取的权矩阵，此处默认各台站提取的地震波到时是等精度的，即P = I_n－1，I为n－1维的单位矩阵。获得dx后，则可得到改正后的震中位置和P波速度；然后将改正后的震中位置由地心地固坐标系转换至大地坐标系，即可获得震中位置的大地经纬度及大地高。

1.3 基于选权迭代法的震中估计方法

在实际震中定位中，仪器噪声、背景噪声等多种因素会对强震动记录产生影响，从而导致STA/LTA+AIC拾取到的P波震相到时不可避免地存在粗差，进一步导致基于到时残差计算得到的震中结果偏差较大。抗差估计是指在存在粗差观测值的情况下，通过改变观测值的权重，使观测值在避免粗差干扰的同时更好地反映真实的测量值^[14]。因此，在传统Geiger方法的基础上，提出一种基于选权迭代法的震中估计新方法。

选权迭代法估计准则为：

$ \sum \boldsymbol{P}_i \boldsymbol{V}_i^2=\min $

(10)

式中，权函数P_i^(m+1)=f(V_i^(m), ⋯)，P_i^(m+1)为第(m+1)步迭代权阵的i分量，V_i^(m)为第m步迭代残差向量的i分量。在选权迭代法估计准则中，通常使用的权函数为IGGⅢ等价权选择准则，具体公式参考文献[15]。基于IGGⅢ抗差估计的震中估计方法的解算步骤为：

1) 首先利用STA/LTA+AIC方法对各强震动台站的加速度记录进行P波到时拾取，并将拾取到P波到时的台站标记为触发状态。

2) 当地震触发的台站大于等于4个时，计算4个台站的重心位置作为震中初始位置(x₀, y₀, z₀)；综合各台站的坐标、P波的理论传播速度V_{P, 0}及台站间差分的P波震相到时观测值，根据IGGⅢ等价权选择准则，计算每个台站的权重；利用每个台站的加权值计算观测值改正数和残差：

$ \left\{\begin{array}{l} \mathrm{d} \boldsymbol{x}=\left(\boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}} \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}} \boldsymbol{l} \\ \boldsymbol{V}=\boldsymbol{A} \mathrm{d} \boldsymbol{x}-\boldsymbol{l} \end{array}\right. $

(11)

3) 由V求解观测值的等价权矩阵$\overline{\boldsymbol{P}} $，应用抗差估计$\mathrm{d} \boldsymbol{x}=\left(\boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}} \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}} \boldsymbol{l} $进行迭代计算。设定阈值ε，当第k次与第k－1次迭代所得改正数之差的绝对值小于等于阈值ε时终止迭代，即

$ \left|\mathrm{d} \boldsymbol{x}^{(k)}-\mathrm{d} \boldsymbol{x}^{(k-1)}\right| \leqslant \boldsymbol{\varepsilon} $

(12)

4) 迭代结束时的改正数、残差计算如下：

$ \left\{\begin{array}{l} \mathrm{d} \boldsymbol{x}^{(k)}=\left(\boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}}^{(k-1)} \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \overline{\boldsymbol{P}}^{(k-1)} \boldsymbol{l} \\ \boldsymbol{V}^{(k)}=\boldsymbol{A} \mathbf{d}\boldsymbol{x}^{(k)}-\boldsymbol{l} \end{array}\right. $

(13)

5) 获得改正数dx后，即可得到震中位置和P波传播速度：

$ \left\{\begin{array}{l} x=x_0+\mathrm{d} X \\ y=y_0+\mathrm{d} Y \\ z=z_0+\mathrm{d} Z \\ V_{\mathrm{P}}=V_{p, 0}+\mathrm{d} V_{\mathrm{P}} \end{array}\right. $

(14)

6) 得到震中位置、地震波传播速度和各台站的震中距，再结合STA/LTA+AIC方法拾取的地震波P波到时信息，即可通过平均法计算主震的发震时刻^[5]：

$ t_0=\frac{1}{n} \sum\limits_{i=1}^n\left(t_i-\frac{D_i}{V_{\mathrm{P}}}\right) $

(15)

基于选权迭代法的震中估计方法流程如图 1所示。

2 震例分析

利用2021年希腊Damasi M_W6.3、2017年中国九寨沟M_W6.5、2021年日本Namie M_W7.1、2012年哥斯达黎加Nicoya M_W7.6、2008年中国汶川M_W7.9和2014年智利Iquique M_W8.2等6次地震事件计算各台站的P波到时、震中位置及发震时刻。以美国地质调查局(USGS)发布的震中位置及发震时刻为标准，评估方法的准确性。6次地震事件详细信息见表 1。

2.1 P波到时提取

对6次地震事件采用STA/LTA+AIC方法提取各台站的P波到时。在进行P波到时处理之前，对原始的加速度记录进行零偏纠正，即将6次事件的全部加速度记录减去各地震前5 s的加速度均值。实际处理中，时间窗长N_S、N_L、触发阈值及AIC方法提取区间的选择会直接影响震相识别的准确性。本文采用STP/LTP+AIC方法的短窗设为1 s，长窗设为10 s，P波到时拾取的阈值设置为4，根据粗略拾取的P波到时T_p将AIC方法区间设置为(T_p-1.5 s，T_p+0.5 s)。以2021年希腊Damasi M_W6.3地震为例，以人工拾取结果为参考，给出各台站的STA/LTA+AIC方法拾取的P波到时结果。为判断P波正确识别的台站数，将识别准确率定义为：

$ \text { ratio }=\frac{M_{\text {Precise }}}{M_{\text {Trigger }}} $

(16)

式中，M_Precise为识别正确的台站数量，M_Trigger为已触发的台站数量，ratio为识别准确率。

在全部30个强震动台站的P波到时拾取结果中，共发生2次误报和1次漏报，ratio值为0.93。误报结果如图 2所示，其中MALA台站的噪声急剧且持续时间较长，而ARS1台站的噪声急剧但持续时间较短。由此可知，STA/LTA+AIC方法对这两类噪声的过滤效果较差，使噪声超阈值而引起误触发。去除误报和漏报的台站后，剩余27个强震动台站的P波到时拾取结果如表 2所示，STA/LTA+AIC方法拾取的P波到时与人工拾取结果相比最大偏差为1.17 s，最小偏差为0.02 s，平均绝对偏差为0.38 s，均方根误差(RMSE)为0.49 s。

2.2 震中估计

将STA/LTA+AIC方法提取的P波到时作为观测值，采用传统Geiger方法和基于选权迭代法的震中估计方法估计6次事件的震中位置。图 3分别给出使用传统Geiger方法和新方法计算得到的震中位置的经纬度及强震动台站分布。

6次事件中，2014年Iquique M_W8.2地震的震中偏差最大(图 3(f))，该事件强震动台站全都位于主震的同一侧；2017年中国九寨沟M_W6.5地震的台站(图 3(b))与2012年哥斯达黎加Nicoya M_W7.6地震的台站数(图 3(d))均为14个，但九寨沟M_W6.5地震的震中偏差要大于Nicoya M_W7.6地震，这是由于Nicoya M_W7.6地震的台站分布在海岸线东北侧的陆地区域，且较为均匀集中，而九寨沟M_W6.5地震虽然在震中附近区域有台站分布，但主震东北侧两个台站距离震中较远，台站分布并不均匀，不均匀的台站分布会导致震中估计结果出现较大偏差。2021年日本Namie M_W7.1地震的台站虽然分布于主震同一侧(图 3(c))，但Namie M_W7.1地震估计的震中位置较为准确，这是由于日本K-NET网络的强震动台站分布密集，所采用的66个触发台站为最小二乘方法提供了更多的多余观测；而2008年中国汶川M_W7.9地震的偏差最小(图 3(e))，这同样是由于汶川地震有50个台站纳入震中估计，且靠近震中的触发台站数目较多。因此，在降低粗差影响后，较多的触发台站数目会提高震中定位的精度。

为定量描述震中估计的结果，以USGS发布的震中位置作为参考，将估计的震中位置与之作差，采用真值“0”作为震中偏差参考值，统计6次地震事件使用两种方法获取的震中偏差和RMSE值。如表 3所示，使用Geiger方法计算2021年Damasi M_W6.3地震的震中偏差为21.26 km、2017年中国九寨沟M_W6.5地震震中偏差为36.79 km、2021年日本Namie M_W7.1地震震中偏差为11.59 km、2012年哥斯达黎加Nicoya M_W7.6地震震中偏差为13.94 km、2008年汶川M_W7.9地震震中偏差为17.45 km、2014年Iquique M_W8.2地震震中偏差为33.88 km，6次地震事件的震中平均偏差为22.48 km，RMSE为19.44 km。使用基于选权迭代法的震中估计新方法计算2021年希腊Damasi M_W6.3地震的震中偏差为13.06 km、2017年中国九寨沟M_W6.5地震震中偏差为22.17 km、2021年日本Namie M_W7.1地震震中偏差为8.46 km、2012年Nicoya M_W7.6地震震中偏差为7.09 km、2008年汶川M_W7.9地震震中偏差为5.08 km、2014年Iquique M_W8.2地震震中偏差为19.73 km，6次事件的震中平均偏差为12.59 km，RMSE为11.00 km。结果表明，相较于Geiger方法，使用选权迭代法计算的震中结果精度提高了43%。

使用Geiger方法和选权迭代法获取震中位置、地震波传播速度和各台站的震中距后，根据式(15)可计算主震的发震时刻，进一步验证新方法的可行性。表 4给出6次地震事件的发震时刻精度统计，可以看出，与USGS公布的理论发震时刻相比，Geiger方法获取的发震时刻最大偏差为6 s，平均偏差为4.1 s，RMSE为4 s。选权迭代法获取的发震时刻最大偏差为4 s，平均偏差为2.3 s，RMSE为2 s。结果表明，基于选权迭代法的震中估计方法得到的发震时刻，无论是最大偏差还是平均偏差都优于传统的Geiger方法，发震时刻获取精度提升50%。

本文采用的基于选权迭代法的震中估计方法对于靠近震中台站分布密集且均匀的汶川M_W7.9地震和台站数量最多的Namie M_W7.1地震的应用效果良好，震中偏差均在9 km以内。但在台站分布较为稀疏不均的九寨沟M_W6.5地震和台站分布集中于主震一侧的Iquique M_W8.2地震的震中估计结果，与参考值仍存在19 km以上的偏差。这是由于采用的6个震例分布于全球不同区域，而不同场地、基墩及环境因素下各震例台站获取的到时精度略有不同，且不同震中距下各震例台站的P波到时观测值对震中位置估计的贡献率也略有不同，从而使得选权迭代法的适用性降低。此外，本文假设地震波在各个方向的传播速度相等，但实际的传播速度因地质结构的不同存在差异，所以震中定位包含由匀速传播假设带来的误差^[13]。综上，在实际应用中, 还需考虑传播速度的各向异性，采用距离震中较近的触发台站精确估计震中位置，对距离震中较远的台站进行剔除或降权处理。

3 结语

本文针对Geiger方法利用的P波到时信息可能存在粗差，导致地震震中估计准确性降低的问题，提出一种基于选权迭代法的震中估计新方法。利用6次地震事件(2021年希腊Damasi M_W6.3、2017年中国九寨沟M_W6.5、2021年日本Namie M_W7.1、2012年哥斯达黎加Nicoya M_W7.6、2008年中国汶川M_W7.9和2014年智利Iquique M_W8.2地震)的实测数据，评估震中位置计算精度。结果表明，使用Geiger方法计算6次地震的震中位置平均偏差为22.48 km、RMSE为19.44 km，获取的发震时刻平均偏差为4.1 s、RMSE为4 s；使用基于选权迭代法的震中估计新方法计算的震中位置平均偏差为12.59 km、RMSE为11.00 km，获取的发震时刻平均偏差为2.3 s、RMSE为2 s。相较于传统Geiger方法，基于选权迭代法的震中估计新方法的精度提高43%，发震时刻获取精度提高50%，可有效减弱P波震相到时识别粗差对震中定位结果的影响，提高震中定位结果的准确性。