收稿日期 2008-06-23收到初稿， 2008-12-04决定采用修改稿. 摘要 ：讨论了全局选择震源初始位置下的矩阵反演求取全局解的问题. 与流行的全局优化搜寻直接定位方法相比， 该反演算法是一种行之有效的定位方法， 具有定位精度高、 计算时间少以及对噪声数据不敏感等特点. 其突出的优点是在不增加计算难度和计算时间的前提下， 用矩阵反演的方法确保得到全局最小值解， 可适应于地震早期预警， 海啸早期预警， 以及大震速报等实际工作. 关键词 ： 射线追踪 全局直接节点搜寻 矩阵反演 全局解 到时残差空间 Fast and accurate earthquake location in 3-D complex medium

1. College of Geology Engineering and Geomatics, Chang’an University, Xi’an 710054, China;
2. Key Laboratory of Western China’s Mineral Resources and Geological Engineering, Ministry of Education, Xi’an 710054, China

Abstract : A new technique of matrix inversion with initial global sampling was developed to seek a global solution. Comparing with the global grid search method it is a practical selection for earthquake location procedure in the sense of solution accuracy, computational efficiency and non-sensitivity to noise data. The dominant feature of the new method is that a global solution can be found with matrix inversion without extra computational time and mathematic complexity required. It is a good choice to the earthquake early warning, tsunami early warning and rapid hazard assessment and emergency response after strong earthquake occurrence. Key words : shortest path ray tracing global direct grid search matrix inversion a global solution RMS error space 地震定位是地震学中经典而又基本的问题之一，以往的一维层状速度模型下的线性地震定位方法已不再能够满足研究诸如地震活动规律、 地球内部结构、 震源破裂几何尺度及构造等的精度要求. 随着宽频带数字地震仪大范围的布设以及多震相到时资料的积累，目前亟需解决复杂速度模型下三维地震的精确定位问题. 地震定位主要可分为两大系列，即直接定位法和反演定位法两种. 所谓直接定位法就是在速度模型已知的条件下，由地震观测到时采用全局搜寻的方法进行地震三要素的直接求解. 例如，Monte Carlo类搜寻及改进方法[如，直接网络结点搜寻法(direct grid search method，简写为DGS)( Nelson，Vidale，1990 )； 相邻搜寻算法(neighbourhood algorithm，简写为NA)( Sambridge，Kennett，2001 )；八象限分割重要采样法(oct-tree important sampling method，简写为OTIS)( Lomax，Curtis，2001 )]. 这种全局优化搜寻的直接定位法可以得到精确的全局解，但缺点是费时.

反演定位法就是在速度模型已知(或未知)的条件下，通过反演的方法求解残差泛函(实际观测到时与理论走时之差)在最小二乘意义下为最小值的解. 例如，多地震同时定位法( Douglas，1976 )，联合定位法( Crosson，1976 )，相对定位或主地震定位法( Spencer，Gubbins，1980 )，以及近年来的双差法( Waldhauser，Ellsworth，2000 )等. 国内学者在此基础上也做了不同程度的改进. 例如，联合定位法的改进( 刘福田，1984 )，昆明地震台网多地震事件定位问题的初步研究( 王椿镛等，1993 )，主地震定位法的修正( 周仕勇等，1999 )，双差定位法的应用研究( 杨智娴，陈运泰，2004 ；Yang et al ，2005)等. 反演方法的优点是省时，缺点是有时只能得到局部的极小值解，特别是在复杂速度模型下更是如此.

然而地震早期预警，海啸早期预警，以及大震速报等则要求定位算法不仅具有定位精度高，而且具有定位速度快的特点. 从某种意义上来说，这里速度是第一位的. 任何不必要的延误将直接导致对受伤群众应急救护的不及时，从而可能造成较大的人员伤亡; 二是要准确，较大的位置偏离，有可能影响震区政府应急救援措施的有效实施，从而造成不必要的社会秩序混乱，以及财产损失和人员伤亡. 因此这里的问题是能否综合全局直接优化搜寻算法全局解的特点，同时又能够吸纳反演方法速度快的优点，从而得到一种不仅定位精度高，而且速度快的算法. 换句话说，就是用矩阵反演算法如何去得到全局极小值解的问题. 这正是本文提出的全局初始值选择条件下的矩阵反演算法. 下面阐述这一反演算法的基本原理.

1 矩阵反演定位全局精确解算法原理

在Monte Carlo类全局搜寻方法基础上发展起来的全局网格优化搜寻直接定位算法，如有限差分算法下的直接网格节点搜寻定位方法(DGS)( Nelson，Vidale，1990 )和八象限分割重要采样法(OTIS)( Lomax，Curtis，2001 ).其基本原理是： 首先进行粗网格单元(或节点)的全局搜寻得到潜在震源节点(或单元)，然后以此节点(或单元)为中心进行局部模型的细化从而搜寻到时残差为最小的节点. 此节点所在的位置就是所求震源空间位置，而发震时间则由各地震台站上的观测到时减去相应的理论走时的平均值给出. 此定位搜寻结果可由到时均方根残差分布表示，也可由概率分布给出( Moser et al， 1992 ). 这类方法的优点是算法简单，但缺点是粗网格单元(或节点)全局搜寻时并不能确保锁定潜在震源节点(或单元). 因此为确保潜在震源节点(或单元)位于局部细化的小模型内，则需减小粗化网格单元(或节点)时的尺度或增加局部小模型的尺度(或范围)，从而增加搜寻所用时间. 鉴于此需另辟蹊径.

然而矩阵反演本质上具有的局域解一直是困扰着地震学乃至地球物理学者的问题，当然人们可用遗传算法、 模拟退火法等具有全局最佳估计的方法去寻找最佳解(全局解). 然而上述全局最佳估计方法在实际应用中存在着这样或那样的问题，如计算量过大、 学习样本太少等，都将增加解的不确定性. 这方面的讨论已超出本文所涉及的问题，这里不再重复. 我们这里的问题是: 在矩阵反演定位算法中是否可以寻求全局解? 答案是肯定的. 众 所周知，非线性反演问题的困难在于其解对初始值的依赖性，即不同的初始值反演后所得的局域解有可能各不相同. 这种非线性反演问题可简化为 图 1 所示，它具有一个全局最小值和多个局部最小值.如果初始值选取不当，则反演后很难收敛至全局最小值，除非采用某种扰动理论克服势垒(局部最大值)的屏蔽从而使初始值接近全局最小值. 然而我们并不知道解空间的细节，因此无法判定初始值是否在全局最小值附近.当然可采用最佳全局搜寻的方法来确定这一点，但这需要大量的尝试实验，这样问题又回到了Monte Carlo实验，在此不讨论. 从 图 1 可看出，只要初始值选取在全局最小值所在沟谷内任意一点，则任何反演算法都能轻松地收敛至全局最小值. 这样问题的关键在于选取位于全局最小值附近的初始值，同时要避免费时的类似Monte Carlo搜寻的方法. 设想采用一定的尺度(或标度)将解空间详细化分成 N 个子空间(或区域)，且这 N 个子空间(或区域)可确保其内涵盖了所有的局域(包括全局)最小值区域，则在这 N 个子空间(或区域)内所选取的 N 个不同的初始值，通过反演可确保其反演解之一为我们所要的全局最小值解. 这就是本文提出的通过矩阵反演寻求全局最小值解的基本原理(global minimum solution，简写为GMS).

然而要实现上述算法，用传统的正演方法(如： shooting and bending methods)( Julian，Gubbins，1977 ； Um，Thurber，1987 )则问题又回到了 N 次单一地震的反演定位问题，无疑大大增加了反演定位所用时间. 如果能将这 N 个不同的初始值定位问题看作是 N 次地震的同时定位问题，则所需时间与单一地震定位时间大体相同.

网络算法的地震射线追踪(modified shortest path method，简写为MSPM)( Bai et al， 2007 )和多地震同时定位方法(multiple event location，简写为MEL)( Bai，Greenhalgh，2006 )为我们解决上述关键问题提供了科学思路. 其具体作法是:

1)用一定的尺度(或规则)划分三维速度模型(实际上也就是三维解空间)得到 N 个子域解空间，并从中选取 N 个相应的待定震源位置初始值.

2)将这 N 个待定震源初始值用多地震同时反演方法(MEL)进行同时定位，此过程中使用相同的观测到时资料.

3)在这 N 个反演解中寻找到时残差最小的解，此解即为可能的全局极小值解.

4)为确保上述解是全局最佳解，在到时残差(如： RMS)空间中进一步改善定位精度.

2 方法有效性验证

下面举例说明这种算法的有效性. 为了模拟实际上较为复杂的速度分布情况，我们选取了 图 2 所示的三维速度模型. 这实际是一个三维变化的速度场，在水平面内由等距相间的高、 低速度异常区组成，并且在各异常区内速度也是变化的. 对低速区而言，由异常中心向外线性递增; 而对高速区来说，则是由异常中心向外线性递减. 模型底面背景速度值为5.75 km/s( 图 2a ). 在垂直方向上速度由模型顶面( z =0.0 km，4.0 km/s)随深度线性增加至模型底面( z =-14.0 km，5.75 km/s)( 图 2b ).

Fig.2 (a)A plan view of the velocity model in horizontal plane( z =-14.0 km). Dashed contours: velocity high; solid contours: velocity low; open circles: surface seismic array; black circles: cross source line projected to the plane; black star: the location of a vertical source line in horizontal plane;(b)velocity contour of vertical section.
Circles: earthquakes projected to the vertical section

在模型顶部中心附近分布有11个地震记录仪( 图 2a )，36个待定震源沿两条构造线分布，其中有21个模拟穿越地震台网下的俯冲断层上的地震分布，而另外15个则模拟垂直断层上的地震分布( 图 2b ). 大多数模拟地震的震源位置位于地震台网之外，特别是垂直断层上的地震分布是远离地震台网的. 这种三维速度模型中地震与台网的极端分布实际上是检验各种定位方法有效性的例子. 为了模拟实际地震的到时记录，可产生36个正的随机数，其最大值为10.0 s. 将这一随机数序列按震源顺序加到相应的理论走时记录中，从而模拟生成36次地震到这11个地震记录仪的无噪声到时记录，而这36个正的随机数即为相应地震的发震时刻. 采用不同的长方体单元均匀地将三维速度模型分别划分成1，8，27，64，125，216，343个和512个子空间，并选取各个长方体单元的几何中心所在的位置作为地震反演定位时的震源初始位置，这样可分别得到1，8，27，64，125，216，343个和512个震源初始位置. 将这一系列震源初始位置用MEL法分别进行反演(初始发震时刻均置于0.0 s)，反演循环次数均为10次. 图 3 给出了28号震源(实际位置： x =5.0 km， y =5.0 km， z =-6.0 km和 t ₀ =5.011 s)一系列震源初始位置反演后的震中位置和各记录台站的均方根(RMS)误差在 x-z 垂直剖面的结果( y-z 垂直剖面形态与此类似，图略). 图 4 则为水平面上的结果.

从 图 3 和 图 4 中可看出，仅用一个随机选取的初始震源位置进行反演，其解的误差是比较大的( 图 3b 和 图 4b )，这也从侧面说明了反演解对初始震源位置的依赖性. 随着初始震源位置数目 N 的不断增加(当 N ≥27时)，反演解的不确定性得到了很好的约束(或控制). 均方根(RMS)误差空间的表现形式则为椭圆型( 图 3 )和直线型( 图 4 )，且具有中间小两头逐渐增大的到时均方根误差分布形态. 这种特点为我们锁定唯一的震源位置提供了科学依据. 另一方面也从本质上验证了这一方法的有效性，即通过这种全局性选取初始震源位置而进行矩阵反演的方法是能够确保所得到的解是全局意义下的最佳解. 实际上随着初始震源位置数目 N 的不断增加(当 N ≥27时)，从反演解中得到的最小到时残差解与最终由到时均方根误差分布锁定的震源空间位置基本相同.

3 定位精度对比分析

为了进一步验证该方法的有效性及解的全局性，我们选取了两种目前较为流行的全局 搜寻定位方法： 一种是有限差分算法下的DGS方法( Nelson，Vidale，1990 )；另一种则是OTIS方法( Lomax，Curtis，2001 ). 其原理在第一部分中已简述，这里不再重复. 这里要说明的是，在OTIS方法中其结果是以概率形式给出，为了与其它两种方法对比，我们将其结果也以到时均方根误差的形式给出. 图 5 给出了用全局反演定位(DGS)方法得到 的定位结果的一个例子(35号震源位置： x =5.0 km， y =5.0 km， z =-13.0 km和 t ₀ =0.414 s)，从 图 5 可看出由到时均方根误差空间分布可锁定唯一的震源位置. 这种方法的优点是简单直观，但缺点是粗网格单元全局搜寻时要采用适度的单元尺度，过小增加不必要的搜寻时间，过大则有可能漏掉潜在震源所在单元节点，这种情况在较为复杂的速度模型及台网布局不合理时容易发生. 另一方面为了确保所得的解是全局最佳解，局部小模型细化时要囊括尽量多的单元(一般是以该单元或节点为中心，向外扩展至少3个单元来建立局部细化的小模型，这样又增加了搜寻时间. 关于这个问题我们在随后讨论的计算运行时中还将涉及.

图 6 给出了用OTIS方法定位得到的一个例子(26号震源位置： x =5.0 km， y =5.0 km， z =-4.0 km和 t ₀ =2.508 s)，从到时均方根误差的空间分布上也很容易锁定待定震源空间位置. 这种方法的优点是每一次细化过程均是在最小到时均方根残差为最小的单元内进行，其搜寻的速度仅为Monte Carlo类全局搜寻方法的1%，但缺点依然是速度太慢. 而上述两种方法不像我们提出的通过全局采样的矩阵反演法那样直接给出震源空间位置和发震时刻，而是先给出震源空间位置，然后发震时刻的确定则是由各记录台站的观测到时减去理论走时的平均值确定.

表 1 给出了以上3种方法对上述36个模拟震源定位的平均误差，从 表 1 可以看出3种方法的定位精度应属同一量级. 这样的定位精度应该说是相当高的，可以满足精定位的要求. 下面我们讨论以上3种方法对噪声到时数据的敏感性，这也是人们关心的问题. 因为好的定位方法应对在容许范围内的具有一定误差的到时资料并不十分敏感，而定位精度的好坏应与使用速度模型的精确程度有关.