0 引言

地震定位的目的是确定震中位置、震源深度、发震时刻和震级等地震基本参数，是地震学的基本问题之一，对于开展地震学的各项研究具有重要意义，是奠基性工作。基于快速准确的地震定位的地震速报，对于震后防震减灾、抗震救险等工作至关重要。因此，对于地震定位方法的研究从未间断，研究重点主要在于如何提高地震定位精度，并提出新的地震定位方法。

中国早期地震定位工作始于北京鹫峰地震台，由李善邦先生于1930年率先开展，为单台地震定位，1953年采用多台站大规模观测数据确定震中（白玲等，2003）。早期地震定位以几何作图法为主（傅淑芳，1991），受观测仪器及硬件等设备制约，发展缓慢。随着科学技术的进步，计算机技术在地震定位中得到广泛发展和应用，计算机软件、硬件方面的升级，使快速计算成为可能，基于科学计算和计算机技术的智能化数值自动定位方法得到迅速发展。

随着人类生活范围的不断扩展，人工地震日益增多，特点是震级偏小。与震级较大的事件相比，微震事件由于释放能量小、信号微弱、接收的波形信噪比差等原因，通过宏观特征很难予以排除或识别，不同分析人员的识别结果会有差异。这些事件混淆于地震目录之中，影响地震记录分析的有效性和准确性。从大量记录信号中识别并剔除人工爆破，建立准确的天然地震事件目录，有助于强震预测及地震学研究。核武器因其毁灭性威力受国际社会关注，监测核试验成为地震学的一个重要任务。当前核试验正向小当量、隐蔽式的地下爆破形式发展，产生的震动也是1种微震事件。为实行全面禁止核试验条约（CTBT），强化中国对于全球核试验的监测、巩固国防，地震定位尤其是微震定位技术的研究，具有更高层次的战略地位和意义。

1 常规地震定位方法 1.1 线性定位法 1.1.1 绝对定位法

绝对定位方法是常规地震定位方法。Geiger（1912）以1905年3月4日印度地震为例，提出P波到时定位概率方法。Geiger定位法可做以下归纳：①给定一个震源位置和发震时刻的初始值；②在当前震源位置和发震时刻建立线性化的观测方程组；③计算残差均方根（RMS）极小化方向上的一个校正矢量；④采用该校正矢量或1个修正量，得到新的震源位置和发震时刻；⑤重复步骤②—④直至得到的解收敛，或满足某些迭代终止判据。

Geiger定位方法的实质是将非线性方程组线性化，并通过最小二乘法求解，该定位方法的思想对地震定位研究产生了重大影响，线性定位方法大多源于对Geiger经典方法的改进，统称为绝对定位方法（田玥等，2002）。如：①联合震源定位法（Joint Hypocenter Determination，简写为JHD），可由地震波走时反演单个地震的震源参数，或一批地震的震源参数及台站校正值，或同时反演震源参数及地球结构（白玲等，2003）；②Lee等（1975）连续给出中国地震台临时报告HYPO78—81系列程序，赵仲和（1983）参与了80、81版本程序的研制，并将HYPO81程序用于北京地震台网；③Backus和Gilbert提出新的反演理论后，Klein（1978）提出HYPOINVERSE算法；④Lienert等（1986）在此基础上进一步得到HYPOCENTER算法；⑤Nelson和Vidale（1990）改进了HYPOINVERSE算法，提出三维速度模型下的QUAKE3D方法。

吴明熙等（1990）和赵卫明等（1992）分别将Geiger经典方法用于禄劝地震和灵武地震序列定位。前者对禄劝地震的部分余震做出精确定位，通过28个余震震中分布成功得出余震发生与北东走向断层密切相关；后者精确定位出灵武地震M_L≥2.0余震，并确定主要断层分布。美国国家地震情报中心的EDR报告、中国地震局地球物理研究所出版的《中国地震台临时报告》和《中国地震年报》及各省地震编目工作，至今仍然采用Geiger经典定位法。

1.1.2 相对定位法

地震台网的分布、地震波到时读数的精确度、可以用来定位的震相、对地壳结构掌握的精细程度等因素均会影响地震绝对定位精度，为进一步提高地震定位精度，常采用相对定位法。

地震序列分布比较集中，震源之间距离远小于它们到台站的距离，同时小于震源区速度不均匀尺度，每次地震到同一台站的路径几乎相同，此时可以采用相对定位法对地震进行精确定位。相对地震定位法可以用于研究与水库断层有关的浅源地震活动（Fitch et al，1974），还可广泛用于研究地震震源深度，例如中东地震震源深度（James and Fitch, 1979）。

Got等（1994）引入多重相对定位法，成功避免了主地震误差传递到待定地震中。每2次地震组成1个地震对，观测并计算走时方程，对于时空丛集的一群地震，运用波形互相关技术，对走时差读数进行校正，精度可达毫秒级（Nicholas and Mariano, 1992），定位误差达到10 m的数量级。马秀芳（1982）使用相对地震定位法，研究安徽霍山震群的震源位置；章光月等（1983）采用该方法对海城地震前震序列进行重新定位；周仕勇等（1999）引入参考台站改进此方法，用主地震定位法重新定位1997年新疆伽师强震群，定位走时残差降到传统方法的30%，提高了地震定位精度。另外，HDC、PEML等处理多事件的定位方法，在实际研究中得到应用。

目前应用较为广泛的相对地震定位方法是由Waldhauser和Ellsworth（2000）提出的双重残差定位法（Double Difference，简写为DD）。其基本原理如下：对于台站k，引入“事件对”（i，j）及双重残差

$ {\rm{d}}r_{ij}^k = r_i^k - r_j^k = \left({t_i^k - {t_{0i}} - T_i^k} \right) - \left({t_j^k - {t_{0j}} - T_j^k} \right) = \left({t_i^k - t_j^k} \right) - \left({{t_{0i}} + T_i^k - {t_{0j}} - T_j^k} \right) $

(1)

计算$ {\rm{d}}r_{ij}^k$可以用绝对到时，也可以用2个事件的到时差。对单个事件定位，有

$ {r_i} = {t_i} - t_0^* - {T_i}\left( {{h^*}} \right) = \delta {t_0} + \frac{{\partial {T_i}}}{{\partial {x_0}}}\delta {x_0} - \frac{{\partial {T_i}}}{{\partial {y_0}}}\delta {y_0} + \frac{{\partial {T_i}}}{{\partial {z_0}}}\delta {z_0} $

(2)

将式（2）分别用于事件i、j，并相减，得到

$ \delta {t_{0i}} + \frac{{\partial T_i^k}}{{\partial {x_0}}}\delta {x_{0i}} + \frac{{\partial T_i^k}}{{\partial {y_0}}}\delta {y_{0i}} + \frac{{\partial T_i^k}}{{\partial {z_0}}}\delta {z_{0i}} - \delta {t_{0j}} - \frac{{\partial T_j^k}}{{\partial {x_0}}}\delta {x_{0j}} - \frac{{\partial T_j^k}}{{\partial {y_0}}}\delta {y_{0j}} + \frac{{\partial T_j^k}}{{\partial {z_0}}}\delta {z_{0j}} = {\rm{d}}r_{ij}^k $

(3)

将式（3）用于所有台站和事件对，反演得到震源的绝对位置，即双重残差法（DDA），也称为“双差定位法”。

双差定位的优势在于不要求主事件，结果不依赖于初始震源位置，可对大空间范围的地震进行精确定位；可以利用互相关分析读取事件到时差，到时数据精确度得到提高；事件对的距离不受限制，适用性增强。双差定位抗干扰性强，如果使用多种相位到时差，定位效果更为显著，是目前较为重要的1种定位方法。

白玲等（2003）翻译了Waldhauser等（2003）的论文，并采用多重相对定位法和波形相关校正，对发生在顺义的地震进行重新定位（白玲等，2003）。杨智娴等（2003）用双差定位法对中国中西部地区地震进行精确定位。吕鹏等（2011）在双差定位中引入波形互相关，得到结合波形互相关的双差定位法，用于2008年汶川地震余震序列定位。黄浩等（2014）将此方法用于2011年盈江地震序列，提高了地震定位精度。黄媛（2008）还探讨了结合波形互相关技术的双差算法在地震定位中的应用及其前景。

1.2 非线性定位法

非线性定位方法对初始位置依赖性弱，不需要求偏导数或逆矩阵，可以避免陷入局部极小点。但缺点也很明显，主要表现为运算量大，效率相对较低，使用频率与线性定位方法相比较低。

Powell法是直接搜索目标函数极小值的有效方法，是1种改进的共轭梯度法（Powell，1964）。此方法不要求计算导数，对初值要求低，一般只需要使用到时最早的台站位置作为初始值，在地震台网快速自动测报中有一定优势。网格搜索法适合求解多个台站和震相，来自多个事件的地震到时、方位和观测的任意数据组（Rodi and Toksoz, 2000）。Rodi等（2004）因此研究了1种新的网格多事件定位方法——GMEL。人工智能的崛起，带动了交叉学科的发展，基于模仿生物界遗传过程的遗传算法逐渐应用于非线性地震定位，万永革等（1995）和周民都等（1999）对遗传算法在确定震源位置中的应用进行了探讨。遗传算法解决了常规定位方法在搜索方面的漏洞，使定位过程中事件的搜索和匹配过程更加高效、稳定。除以上方法外，非线性定位方法还包括网格搜索法、单纯形法、模拟退火法、蒙特卡罗法等。

非线性定位方法大部分仍在研究分析阶段，只有Powell法和遗传算法在实际地震台网观测定位中得到应用。随着科学技术的发展，计算机性能逐渐提高，非线性定位方法运算量大、效率较低的缺点将被克服，可能会成为常规地震定位方法。

2 微震定位方法

传统地震监测方法主要通过震相识别和到时拾取来实现。微震事件因震级小、信号能量弱（甚至噪音量级），导致初至时间拾取困难，传统方法在微震识别上经常面临失败。因此，微震事件的识别和定位被认为是当前地震学的热点和难点之一。

目前微震定位技术主要有3种基本方法：P波定位法、地震波射线法、P波射线传播方向交汇点法。P波在微震定位中有重要作用，因为P波传播速度快，初至时间易于识别，所以一般情况下应采用P波定位（姜秀娣等，2005）。采用P波定位法进行微地震定位时，可以假设P波传播速度为已知量，采用均匀速度模型，将地震监测台站同时布设在至少4个以上不同地点（潘科等，1997；平键等，2010）。根据各个监测台站坐标和P波到达台站的时刻列出方程组，并通过最小二乘法求解，即可得到震源坐标和发震时刻。目前较为普遍的方法是利用监测台站地表阵列勘测微震的发射断层扫描技术。Kiselevitch等（1991）定义了1种通过将时移信号规范化为平均时间信号的方法，用到的时移是阵列中地震台站的到时差，是P波定位法的最早应用。

本文将重点介绍目前常用的2种微震定位方法：波形匹配滤波器法（Matched Filter Technique）和震源扫描叠加算法（Source-scanning algorithm, SSA）。

2.1 波形匹配滤波器法和波形互相关技术

波形匹配滤波器法以参考信号f(t)为模板，对时间序列u(t)做互相关，其中u(t)是连续的原始三分量数据，r(t)是包含P波、S波到时信息的加窗数据。对1个台阵数据接收器而言，每个接收器均要在3个分量上做互相关，建立1个多分量互相关记录剖面。对2个具有相似波形的事件，通过一定时间间隔的时移，使各事件离散的P波和S波相互重叠成1个单独的互相关峰值，然后利用自动增益方程（Automatic gain control，AGC）避免错误判断，它是原始数据经三角化滤波后的复振幅包络的卷积。对接收器各分量方向建立的多分量互相关记录剖面求和，可以得到叠加后的互相关记录，超过指定阈值的事件会被锁定。若有多个相关参考事件，使用相关系数最高的。

波形匹配滤波器法将事件波形信号叠加后的互相关图像与模板事件匹配，定位微震事件（Gibbons and Ringdal, 2006；Gibbons et al，2007；Shelly et al，2007；Peng and Zhao, 2009）。利用波形互相关技术检测具有相似波形的微弱信号，对于侦测低信噪比事件较有效，已经广泛应用于侦测低频地震震动（Shelly et al，2007）、余震（Peng and Zhao, 2009）和水力压裂引起的微震事件（Eisner et al，2008）。

Grigoli和Cesca用Waveform Coherence Analysis对地震进行自动定位，基本原理如下。

（1）计算特征方程。特征方程通常对每个台站的能量变化、频率组成、偏振或目标信号的背景噪声较敏感（Lomax et al，2012）。这里需要2个特征方程，分别对P波到时、S波初至震相敏感。三分量地震道可被视为离散的时间序列。将EW分量表示为x(j)，NS分量表示为y (j)，垂直分量表示为z(j)，其中j = 1，j为整数，n_samples是序列的时间指标，同时代表了最后1个样本轨迹。根据Grigoli等（2013）提出的方法，定义P波特征方程CF^P为

$ C{F^{\rm{P}}}\left(j \right) = {z^2}\left(j \right) $

(4)

计算S波特征方程需先求水平分量曲线的解析式，即

$ X\left(j \right) = x\left(j \right) + iH\left\{ {x\left(j \right)} \right\} $

(5)

$ Y\left(j \right) = y\left(j \right) + iH\left\{ {y\left(j \right)} \right\} $

(6)

式中，H为希尔伯特变换，i² = -1。根据Vidale（1986）提出的理论，瞬时协方差矩阵Q(j)为

$ Q\left(j \right) = \left({\begin{array}{*{20}{c}} {X\left(j \right)\hat X\left(j \right)}&{X\left(j \right)\hat Y\left(j \right)}\\ {Y\left(j \right)\hat X\left(j \right)}&{Y\left(j \right)\hat Y\left(j \right)} \end{array}} \right) $

(7)

式中，^代表复共轭性。对每个样本j，协方差矩阵Q(j)有2个真实的正特征值λ₁和λ₂，且λ₁≥λ₂。在区域尺度上并不能保证入射的S波垂直，计算协方差矩阵应使用三分量（Rowe et al，2002），但为获取更好的结果依然仅使用基于水平分量的协方差矩阵。则S波特征方程为

$ C{F^{\rm{P}}}(j) = \lambda {(j)^2} + \varepsilon $

(8)

式中，ε是数值小的整数，可以克服计算STA/LTA时λ₁(j)趋于0的问题。最后，采用Allen（1978，1982）提出的算法，分别用CF^P和CF^S计算STA/LTA，采用递归方法，降低内存需求，在缺乏信号的情况下，结果比标准的STA/LTA更平滑。设n_s和n₁分别表示短时窗和长时窗样本数，则STA/LTA的递归算法可以写作

$ {\rm{STA}}(j){\rm{ }} = {K_{\rm{s}}}[CF(j)]{\rm{ }} + {\rm{ }}(1 - {K_{\rm{s}}}){\rm{STA}}(j - 1) $

(9)

$ {\rm{LTA}}(j) = {K_{\rm{l}}}[CF(j - {n_{\rm{s}}} - 1)]{\rm{ }} + {\rm{ }}(1 - {K_{\rm{l}}}){\rm{LTA}}(j - 1) $

(10)

式中，j在h = n_s + n_l和特征方程最后一个样本nsamples之间变化。衰减常数K_l和K_s分别设置为1/n_l和1/n_s（Withers，1999）。STA和LTA等式方程分别表示2个单极、低频带滤波器，在时间域中，滤波器常数分别为K_s和K_l（Baer and Kradolfer, 1987）。则

$ {W^{\rm{P}}}\left(j \right) = \frac{{{\rm{ST}}{{\rm{A}}^{\rm{P}}}\left(j \right)}}{{{\rm{LT}}{{\rm{A}}^{\rm{P}}}\left(j \right)}} $

(11)

$ {W^{\rm{S}}}\left(j \right) = \frac{{{\rm{ST}}{{\rm{A}}^{\rm{S}}}\left(j \right)}}{{{\rm{LT}}{{\rm{A}}^{\rm{S}}}\left(j \right)}} $

(12)

式中，W^P和W^S分别为P波和S波特征方程的STA/LTA，注意传播效应，将STA/LTA标准化，避免台站靠近震源而控制叠加过程，也平衡了P波和S波的权重。

（2）进行定位。将地震台网中所有地震台站的P波和S波特征方程的STA/LTA作为输入值，假设1个地震事件被N个三分量地震台站接收，首先定义1个包含全部孕震区的3D笛卡尔网格空间（Allen，1982）。对坐标为（n_x，n_y，n_z）的网格点，有

$ x\left(l \right) = {x_{{\rm{ref}}}} + l\delta x \left({l = 0, 1, \cdots, {n_x}} \right) $

(13)

$ y\left(m \right) = {y_{{\rm{ref}}}} + m\delta y \left({m = 0, 1, \cdots, {n_y}} \right) $

(14)

$ z\left(n \right) = {z_{{\rm{ref}}}} + n\delta z \left({n = 0, 1, \cdots, {n_z}} \right) $

(15)

式中，δ_x、δ_y和δ_z为x、y和z方向的网格间距（分别与东西向、南北向、垂直向一致），而x_ref、y_ref和z_ref表示参考点的笛卡尔坐标。每个网格点代表 1个潜在的震源定位，对于每个可能的试验震源位置，计算记录台网中N个台站的P波和S波初至的理论到时分别为[τ_k^P(x, y, z)]和[τ_k^S(x, y, z)]（k为台站号）。定义τ_min和τ_max为P波到时最小值和S波到时最大值，则

$ {\tau _{\min }}\left({l, m, n} \right) = \min \left({\left\{ {\tau _k^{\rm{P}}\left({l, m, n} \right)} \right\}_{k = 1}^N} \right) $

(16)

$ {\tau _{\max }}\left({l, m, n} \right) = \max \left({\left\{ {\tau _k^{\rm{S}}\left({l, m, n} \right)} \right\}_{k = 1}^N} \right) $

(17)

引入$ T_k^{\rm{P}}\left({l, m, n} \right)$、$ T_k^{\rm{S}}\left({l, m, n} \right)$，则

$ T_k^{\rm{P}}\left({l, m, n} \right) = \tau _k^{\rm{P}}\left({l, m, n} \right) - {\tau _{\min }}\left({l, m, n} \right) $

(18)

$ T_k^{\rm{S}}\left({l, m, n} \right) = \tau _k^{\rm{S}}\left({l, m, n} \right) - {\tau _{\max }}\left({l, m, n} \right) $

(19)

考虑到记录波形取样率δt，将上述方程离散化，为

$ \Delta T_k^{\rm{P}}\left({l, m, n} \right) = {\rm{round}}\left\{ {\frac{{T_k^{\rm{P}}\left({l, m, n} \right)}}{{\delta t}}} \right\} $

(20)

$ \Delta T_k^{\rm{S}}\left({l, m, n} \right) = {\rm{round}}\left\{ {\frac{{T_k^{\rm{S}}\left({l, m, n} \right)}}{{\delta t}}} \right\} $

(21)

利用等式（20）和（21），采用以下方程求得各网格点和时间样本的相关方程C^P和C^S，即

$ {C^{\rm{P}}}\left({l, m, n, j} \right) = \sum {_{k = 1}^NW_k^{\rm{P}}\left[ {j + \Delta T_k^{\rm{P}}\left({x, y, z} \right)} \right]} $

(22)

$ {C^{\rm{S}}}\left({l, m, n, j} \right) = \sum {_{k = 1}^NW_k^{\rm{S}}\left[ {j + \Delta T_k^{\rm{S}}\left({x, y, z} \right)} \right]} $

(23)

式中，W_k^P为第k个台站的P波特征方程规范化的STA/LTA，W_k^S类似。将Kao等（2004）的方法加以修改，选择1个以ΔT_k^P（或ΔT_k^S）为中心的时窗，将所有样本进行叠加，一致性矩阵最终可以定义为

$ C\left({l, m, n, j} \right) = \sqrt {\frac{{{C^{\rm{P}}}\left({l, m, n, j} \right){C^{\rm{S}}}\left({l, m, n, j} \right)}}{N}} $

(24)

由上述方程可以清晰看出C(l, m, n, j)是边界方程，理论极限是0（不一致）和1（与P波、S波初至震相完美一致）。取矩阵最大值

$ C\left({\hat l, \hat m, \hat n, \hat j} \right) = \max \left\{ {C\left({l, m, n, j} \right)} \right\} $

(25)

最终得到地震时间坐标为$\left({\hat x, \hat y, \hat z} \right) = \left({x\left(l \right), y\left({\hat m} \right), z\left({\hat n} \right)} \right) $，发震时间可以通过以下方程反演得到

$ \hat t = \hat j\delta t - {\tau _{\min }}\left({\hat l, \hat m, \hat n} \right) $

(26)

此方法用于意大利南部的Irpinia区域，取得较好结果。Enrico Caffagni和David W Eaton（2015）使用加拿大西部2个真实微震的数据库对MFA进行检验，并用自动增益函数（Automatic Gain Control Function）去除错误结果。这种基于模板、监测小震级事件的方法被广泛用于地震断层系统（Nadeau and Johnson, 1998；Igarashi et al，2003）和大型内陆板块研究（Schaff and Richard, 2004）。相似方法在微震活动背景中也有广泛应用，如：在Salah石油开采项目中，Goertz-Allman等（2004）用波形互相关方法检测并定位超过5 000多个微震事件；Steven J Gibbons和Frode Ringdal（2006）用基于台阵的波形互相关法检测并定位小震级事件，以Barentsburg煤矿为案例讲述该方法；Oprsal和Eisner（2014）利用规范化互相关（Normalized Cross-Correlation，NNC）作为判据识别天然地震和诱发地震。在中国，类似方法也在各个领域广泛应用，郑晨等（2015）利用波形互相关法识别并重新定位灌县—安县断裂带重复地震；张淼和温联星（2014）将多台站互相关波形图进行自动叠加，用Stacking Multiple-Station Autocorrelograms（SMAC）方法测量地震震源深度；Wen Lianxing等（2010）、Zhang Miao等（2013）用类似方法对朝鲜多次核试验进行重新定位。

2.2 震源扫描叠加算法（Source-Scanning Algorithm，SSA）

SSA方法基于潜在震相信号特征能量叠加，如：绝对振幅，能量包络或地震波短时平均和长时平均比（STA/LTA）等，是1种新的可以对震源分布进行成像的方法（Kao and Shan, 2004）。该方法仅利用相对振幅和到时等波形数据，就可以知道目的区域内是否有地震震源，而不需预先知道断层的产状和规模。其优势在于，不需要拾取到时，也不用计算理论地震图，只需通过系统扫描一定位置范围和时段，就可以获得整个震源序列分布。

Kao和Shan（2004，2007）详细介绍了震源扫描叠加算法，该方法基本原理如下。

SSA通过系统扫描整个模型空间寻找“亮”点来识别震源。空间中1个给定点（η）在特定时间（τ）的“亮度”定义为

$ br\left({\eta, \tau } \right) = \frac{1}{N}\sum {_{n = 1}^N} \left| {{u_n}\left({\tau + {t_{\eta n}}} \right)} \right| $

(27)

式中，u_n是台站n记录到的标准化地震记录振幅，t_ηn是某一震相（在地震动定位里，使用S震相）从点η到台站n的预测走时。

实际上，由于对速度模型了解不完全，在计算t_ηn时，每个台站的预测到时和观测到时可能有微小误差。因此，使用预测到时振幅修改式（26），选定1个时窗，考虑其周围点，得

$ br\left({\eta, \tau } \right) = \frac{1}{N}\sum {_{n = 1}^N} \left\{ {\frac{{\sum {_{m = - M}^M{W_m}\left| {{u_n}\left({\tau + {t_{\eta n}} + m\delta t} \right)} \right|} }}{{\sum {_{m = - M}^M{W_m}} }}} \right\} $

(28)

式中，2M是以预测到时为中心选定时窗内的样本数，δt为取样间隔，W_m为根据观测到时和预测到时差决定的权重因子。

将SSA应用于识别地震破裂时，原始方程（26）、（27）不再合适，主要有2个原因：①不像大多数具有相对较短振动时间的震动，地震震源结构在时间和空间上更加复杂，那些亮度不太高的点并不能完全忽略掉；②因构造产生的信号波形，如尾波、不连续面产生反射、折射波等，会混在震源信号中产生虚像。为克服这些问题，有必要对震源引入先验约束。对于第1个问题，可以对等式（27）添加1个台站校正项，这样源项的预测到时和观测值完美一致；对于第2个问题，可以通过仔细选取跟构造震相无关的波形成分来解决。综合这2项改动，亮度函数写成

$ br\left({\eta, \tau } \right) = \frac{1}{N}\sum {_{n = 1}^N} \left\{ {\frac{{\sum {_{m = - M}^M{W_m}\left| {{U_n}\left({\tau + {t_{\eta n}} + m\delta t + t_m^{{\rm{corr}}}} \right)} \right|} }}{{\sum {_{m = - M}^M{W_m}} }}} \right\} $

(29)

式中，U_n是台站n记录的P波规范波形包络，是对台站n的时间校正，如：用P波观测到时减去预测到时。

中国南北地震带部分地区震源扫描表明，用SSA定位地震的准确性较高，定位结果与实际波形较符合，说明SSA方法具有一定稳定性，可以用于复杂的实际地震定位，特别是在未拾取震相情况下，SSA方法更具优势。丁宁霞和张晓清（2013）利用震源扫描算法，计算2008年大柴旦6.3级地震震源破裂过程，得到理想结果。然而，较少单独使用震源扫描叠加算法进行地震定位，实际应用中常结合其他方法，辅助定位。

2.3 Match and Locate（M & ;L）和HypoRelocate法

基于波形互相关技术的波形匹配滤波器法和震源扫描叠加算法，在高精度定位中均具有一定优越性，是目前应用较为普遍的地震定位方法。然而，这2种方法各自存在局限性。当地震信号被噪声淹没或与其他信号混合后，可能会限制SSA对小震级事件的定位能力。此外，震源定位中对于走时的校正，过于依赖速度模型的精准性。另一方面，波形匹配滤波器法要求小震级事件与模板事件位置相关。使用单一方法进行地震定位的精度和可靠性有待改进，多种定位方法相结合，取长补短，成为一种解决办法。

张淼和温联星提出一种新的定位方法——Match and Locate（M & ;L），就是将波形匹配滤波器法和震源扫描叠加方法相结合，利用各自优点，对微震事件进行识别定位，进一步提高地震定位精度。与波形匹配滤波器法相似，均需大量样本事件，使用多台多分量连续波形，将样本事件和潜在小震级事件互相关图进行叠加；不同之处在于，M & ;L通过扫描模板事件周围的3D区域，确定模板事件和潜在事件之间的相对位置，进行相对走时校正，再进行叠加。将M & ;L应用于日本3次地震事件，得到较好的定位效果（Zhang and Wen, 2015）。

近期，张淼和温联星利用传统地震波绝对走时和走时残差，加入震源间叠加的尾波互相关信息作为约束，提出地震精定位新方法HypoRelocate（Sun et al，2016）。该方法通过模拟退火法，寻找以上3种约束总残差的最小值，实现对地震最佳位置和发震时刻的反演。在地震波走时差精确约束地震相对位置基础上，通过构建震源间叠加的地震尾波互相关波形，更加精确地约束震源间距；通过加入地震波绝对走时，削弱地震定位结果对初始地震目录位置的依赖。可以预见，多种定位方法相互约束、相互补充的联合定位法，将成为今后地震定位，特别是微震定位的一个重要发展方向，多种方法相结合的精定位方法将大量涌现。

3 微震定位方法在核试验监测领域的应用

核试验可以在大气层中进行，也可以在水中进行，而更为安全有效、隐蔽的则是在地下进行。《全面禁止核试验条约》设立以来，提出可以使用不同技术（地震波、水声、次声以及放射性核素）监测核试验。监测地下核试验是所有类型核试验探测中具挑战性的工作，而地震监测被认为是重要的有效监测手段。Waldhauser等（2004）就曾使用双差法进行核爆事件识别定位工作。对于在远距离上观测到的震级较大的地震事件，可以使用传统地震定位方法。然而，伴随着高质量区域地震数据的出现，很多人工震源引发的小地震事件（m_b等于4或更小）逐渐被检测到，混淆了人们对于目标事件的判断。当前，对当量较大的核爆炸已有较为可靠的监测与识别方法，但对于小当量核爆，由于地壳介质的复杂性及非弹性特征引起几何扩散和频散作用，导致地震波传播过程并不稳定，特别是介质对不同频率的波的吸收和衰减作用，致使小当量核爆炸识别更加困难。识别核爆炸并定位，尤其是小当量核爆炸，已成为爆破地震学的重要内容。在现代核试验背景下，各国进行核试验的方式更加隐蔽，当量越来越小，趋于微震化。因此，增强微震识别定位能力对核爆破研究显得尤为重要。

由于核试验的特殊性质，一个国家的核试验往往在固定区域进行，如：美国内华达、巴基斯坦的科-坎巴拉、印度的波克兰以及朝鲜丰溪里核试验场等，地域上存在极高重复性。这些特点与基于波形匹配技术的各种微震定位方法具有相当高的契合度：①在同一区域反复进行核试验，地质背景相同，可以很大程度地解决对速度模型依赖的问题；②事件之间匹配度较高，事件性质相同。因此，微震定位研究对于核试验监测具有重要意义。

2006年10月9日以来朝鲜相继开展几次核试验，张淼、温联星等利用微震识别方法，对几次朝鲜核爆事件进行重新定位（Wen and Long；2010；张淼等，2013；Zhang and Wen, 2015），结果认为：2009年朝鲜核爆事件的精确定位为（41.293 9°N，129.081 7°E），精度达到140 m；2013年核爆事件定位为（41.2908°N，129.076 3°E），精度为94 m，较2009年定位精度有所提升；针对2016年核试验，温联星研究组在其网站上给出准确定位结果，且精度得到进一步提升，达到50 m（http://seis.ustc.edu.cn/News/201601/t20160106_234889.html，http://seis.ustc.edu.cn/News/201609/t20160909_253322.html）。中国科学院地质与地球物理研究所赵连锋等(2013, 2014)对2009年朝鲜核爆事件的精定位结果为（41.292 8°N，129.084 9°E），与温联星研究组的结果接近。潘常州等（2014）也对朝鲜核爆事件做了大量研究。

高精度定位结果充分说明了微震定位方法在核试验监测工作中的重要作用。在小震级事件难以被检测、识别和定位的背景下，微震定位技术可能是目前核爆炸监测的重要手段之一。

4 结束语

地震定位作为地震研究的一项核心工作，具有重要意义。常规地震定位方法主要有线性定位和非线性定位法，其中线性定位法中又包括绝对定位法和相对定位法。Geiger经典定位法对地震定位影响较大，许多定位方法在此基础上发展而来。目前，双差地震定位法应用比较广泛。对于微震定位，基于波形互相关技术和震源扫描叠加算法，或单独使用或相互结合的定位方法较为常见，具有更高的定位精度。类似于M & L联合定位法的出现，将地震定位精度推向新的高度，将成为今后地震定位发展的主流方向之一。

传统地震定位方法的实质在于求目标函数极小值，各种方法对于目标函数的构造、处理以及求极小值的方法不同。影响定位精度的因素多种多样，客观因素有地震台站分布和密度、观测仪器精度、速度模型的可靠性等，主观因素有波形和震相拾取判别、阈值设定等。各种定位方法针对其中某个问题进行设计，各有优缺点。现在比较流行的定位方法多依赖波形匹配，对于波形记录比较完整、连续的区域可以较好地应用，而对于一些波形记录较少或难以筛选样本事件的区域，可能需要先进行相当时间的事件观测和记录，联合定位方法的出现，会使地震定位中存在的一些问题得到解决。微震监测识别定位技术逐渐显示其重要地位，但在中国起步较晚，在今后研究中将面临许多问题，如：微震数据采集系统的建立、对微震波形和产生机制的研究以及震源位置及发生时刻的精确判断等，微震定位要作为一种成熟的、常规的监测手段还有很长的路要走。

总体来讲，地震定位方法的发展日趋精细化、综合化：由传统的使用单一方法进行地震定位，发展为多种方法联合定位，提高了定位精度，减小了误差，消除了速度模型等因素对地震定位的影响。微震定位也成为了地震定位方法的热点之一。地震定位的功能性变的更加多样，微震定位在核试验监测方面将有更深入的应用。在今后地震学研究发展中，地震定位方法将更加多样。