考虑到时误差的地震定位算法及其在四川地区2001—2008年地震定位的应用

万永革 盛书中 程万正 张致伟 武 晔赵晓燕 卜玉菲 薛志芳 刘金丽

1)防灾科技学院，河北三河燕郊开发区 065201

2)四川省地震局，成都 610041

3)甘肃省地震局，兰州 730000

4)河北省地震局秦皇岛中心台，秦皇岛 066100

5)河北省曲周县实验中学，邯郸 057250

考虑到时误差的地震定位算法及其在四川地区2001—2008年地震定位的应用

万永革1)盛书中1)程万正2)张致伟2)武 晔1)赵晓燕1)卜玉菲3)薛志芳4)刘金丽5)

1)防灾科技学院，河北三河燕郊开发区 065201

2)四川省地震局，成都 610041

3)甘肃省地震局，兰州 730000

4)河北省地震局秦皇岛中心台，秦皇岛 066100

5)河北省曲周县实验中学，邯郸 057250

精确的地震位置对于地震活动性、地震层析成像和地壳应力场反演具有相当重要的意义，对于地震速报也具有重要的应用价值。将观测到时的不确定性、台站高程、地震震源深度进行约束的同时，根据反演理论给出了地震震源位置精确估计和误差估计的方法。该算法联合考虑Pg波、Sg波、Pn波和Sn波的到时进行反演，数据量的增加可以增强地震位置的准确性，并可同时应用于地方震和区域地震。采用模拟数据对该地震定位算法进行检验发现，该算法在观测数据的不确定性不等时明显优于其他方法。将该算法应用于四川地区2001—2008年间的地震定位，得到的地震位置更加符合地震的丛集性并集中于断裂带附近。这些结果为四川地区的地震活动性、断层构造以及地震层析成像研究打下了基础，并且为汶川地震之前的地震活动前兆研究也提供了有益帮助。

地震定位 不确定性 台站高程 四川地区

0 引言

地震定位是地震学中最基本的问题之一，对于研究地震活动构造、地球内部结构、震源的几何构造等地震学基本问题都有重要意义。另外，基于快速准确的地震定位的地震速报，对于震后的减灾、救灾工作也具有重要的实际意义。因此，地震学家一直在不断改进或提出新的定位方法。

地震定位算法有很多种(Lee et al.，1975;赵仲和，1983)。一些研究者采用S波与P波的到时差和虚波速度确定地震的位置，如Inglada法。该方法的优点在于不需要给定初始震中位置，但采用了与实际相差较大的均匀地壳模型(徐果明等，1982)。最常用的地震定位算法为Geiger法(Geiger，1912)。这种方法可采用较为复杂的地壳模型，甚至三维不均匀速度结构(Wu et al.，2008)，其好处是能够根据与实际相近的地壳模型反演震源位置，但必须给出与真实震源位置相差不大的初始震源位置和发震时刻。还有的研究者是采用遗传算法(万永革等，1995，1997)和单纯形法(赵珠等，1994)，这类方法的好处是不需要求理论走时对震源参数的偏导数，缺点是搜索较慢，且难以估计参数的不确定性。还有一些研究者采用相对定位法，如主事件定位法(Spence et al.，1980;周仕勇等，1999)、双差定位法(Waldhauser et al.，2000)。主事件定位法必须给出较为准确的主事件位置，并且其他事件与主事件的空间距离要比震源到台站的距离小很多。双差定位法可以得到事件之间的相对位置，其绝对位置的准确性也能得到保证，但该方法不能对零散地震进行定位。

牟磊育等(2006)研究了原始的Geiger法，给出了联合采用Inglada法和Geiger法对区域台网地震到时数据的精确定位，编写了MATLAB程序。虽然该方法是一种较为理想的区域台网定位法，但该程序没有利用Pn波到时，没有考虑各种震相到时拾取准确性的差别，并且不能给出定位参数的不确定范围。本文在他们工作的基础上增加了Pn波到时的利用，并且考虑了不同震相到时拾取的不同误差，利用反演理论给出了估计震源参数的标准差。将我们研制的方法应用于四川地区2001—2008年的地震定位，得到了较为精确的地震位置。

1 定位算法及其改进

1.1 考虑观测误差的Geiger法反演地震震源参数

4个参数描述了发震时间和震源位置，我们把这些参数称为模型，定义模型矢量

采用Geiger法是通过考虑目标位置的扰动把问题线性化

这里m0是猜想的最佳位置，可采用Inglada法求出粗略的震源发震时间和位置(徐果明等，1982)。m是离m0很小距离的新的位置。可以用泰勒级数展开式的第1项来近似表示在m处所预测的时间

为了使这些残差最小，我们试图用下列方式求解模型的修正量Δm

这里G 是偏导数Gij=∂ti/∂mj的矩阵，i=1，2，…，n，j=1，…，4。 为了得到位置的调整值 Δm，用标准的最小二乘法得到方程(4)的最佳拟合。下一步，我们使m0处于m0+Δm的位置，重复这一过程，直到定位收敛。假如初始猜想的位置离实际位置不是太远，迭代过程一般能相当快地收敛。

考虑观测误差时，式(6)可以写为

其中，δd为观测误差，对于Pg波，Sg波和Pn波可以有不同的误差，Cd为观测误差的协方差矩阵。利用最大似然解估计公式可得(Jackson et al.，1985;Shen et al.，1996，万永革等，2004，2008a，b)

解的协方差矩阵可表示为

该矩阵对角线元素的开方即为解的标准差。

平均残差定义为

其中m为数据个数，即方程的个数。

1.2 Pn，Sn波走时及其偏导数

设地球平层模型的地壳层的序列编号为l，最下一层为地幔层m，地壳层厚度为dl，速度为Vl，入射角为θl，地幔层的速度为Vm。震源位于j层中，dz为震源层中震源距该层顶部的距离。根据该地球模型，只有震源层的入射角满足时才会出现首波，在其他层中的入射角为

临界观测到首波的距离为

震中距为Δ的首波的走时为

Pn波走时的偏导数

1.3 台站高程的校正

在平原或高原地区，台站高程相差不大，对定位没有太大影响。但如果在地形起伏比较大的地区，台站高程往往相差较大，如果不考虑台站高程会导致定位结果误差较大。牟磊育等(2006)的程序因用于平原地区而没有采用台站高程校正。本研究在计算某一台站走时时，采用地壳模型第1层厚度增加台站高程值的方法来校正台站高程的影响，使得定位结果更为准确。

1.4 地震震源深度的约束

其他定位方法一般没有对深度进行约束。只是限定如果深度不满足一定的范围(如空中地震或深度特别大地震)，则强制给出某一深度值。这样处理的结果会导致震源参数不满足最小残差的要求。我们在每一次迭代过程中判断修正后的震源深度是否在合理的范围之内，如果不合理，只将参数增加到合理深度的边界，进行下一步迭代。下一步迭代还有可能使得震源深度向较为合理的方向进展。最终得到残差最小且深度在合理范围内的震源参数。

1.5 地震深度的精确求解

由于地震震相数据观测数据有限，并且地震台站大多位于同一水平面上，对深度的分辨率往往不高。为了使深度有更好的约束，我们分别在初始水平位置的基础上设置每一层的中点作为深度的初始深度进行迭代求解，选择残差最小的解为最优解。

2 数值实验

为了测试以上的定位方法，假设一个由9个台站组成的虚拟台网，虚拟台网之形状为边长300km的正方形(图1)。虚拟台网下的速度模型按照地球物理较为全面的观测结果得到的Crust5.1(Mooney et al.，1998)给定。具体参数见表 1。

图1 虚拟台站及虚拟地震位置Fig.1 Locations of the virtual stations and earthquake.台站的编号同表1

假设在虚拟台网内发生一个地震，发震时间 0时刻，位于横轴 75.071km，纵轴50km处，深度为19km(图1)。为简单起见，虚拟台网的高程全部设为0，根据虚拟台网内台站的位置坐标、速度模型、震中位置和发震时间，计算出虚拟台网内各台站的Pg、Sg和Pn的到时，分别随机产生以0.5，2.5，1.5幅度的随机干扰，得到模拟 Pg、Sg和 Pn的到时(表2)，然后采用我们修订的程序进行定位计算。初始位置X，Y，Z，t分别设定为50，30，18，0。当采用相同的到时拾取误差时，得到的地震震源 X，Y，Z，t分别为 76.0866，50.1371，18.1681，0.1544;采用与虚拟随机扰动相同的误差进行反演得到的 X，Y，Z，t分别为 75.6551km、50.1251km、17.6933km、0.1366s，得到的该参数的标准差为 1.4504km、1.5413km、4.2750km 和0.1958s。采用不同的到时拾取误差得到的结果更加接近虚拟设置的地震位置，可见引用相对到时标准差的重要性。

3 2001—2008年的四川地区地震定位实例

采用我们的程序对四川地区2001—2008年的地震进行了研究，对震相进行核对，得到1,831个地震的位置(图2)。原来所有地震震相的平均残差减少为原来的0.35%。定位后的位置更加集中在断裂带上，并具有更强的丛集性，反映了小震分布的基本特征。重新定位后的小震丛集区主要为龙门山断裂带、鲜水河断裂带及其西南区域、马边-盐津断裂带附近及华蓥山断裂带北部和东南区域。很多研究者采用双差定位法来确定地震的丛集特性(Waldhauser et al.，2000)，其重定位是将地震组成地震簇而精确确定地震之前的相对位置，由于去除了不能组成地震簇的离散地震，因此能反映地震的丛集特性从而确定断层的形状(万永革等，2008a;王福昌等，2008)。而本研究没有去除离散地震，也明显表现出地震的丛集特性，说明了研究结果的客观性。赵珠等(1997)采用该地区早期较长时间的地震资料确定的地震位置也表现为地震丛集在龙门山断裂带附近的特性，与本文的研究结果有类似之处，表明了本文处理方法的正确性。

在深度上，地震目录列出的某些震源深度达60km左右(已达不易发生地震的上地幔)，重新定位后的震源深度主要集中在8km以上的上地壳(图3)，下地壳的地震较少。这与Takagi等(1997)的结果类似。另外，由于这期间大部分为小震，均发生在浅部的4～8km的范围内，深部缺少地震，似乎可以称之为深部“地震空区”，这也许是地震预测研究人员要注意的现象。

4 结论与讨论

与其他地震定位方法对比，我们的方法有如下优势:1)可以利用Pn波到时数据;2)可以对数据给予不同的观测误差，从而确定地震定位的不确定性;3)考虑了台站高程;4)可以约束震源深度的搜索范围;5)通过对各个层中分别设置震源来去除由于地壳结构各层速度差别较大而造成非最优解的情况。实验结果表明，在数据具有较大误差的情况下，我们的程序仍然可以得到较为精确的震源位置。将该算法用于2001年至汶川地震之前的四川地区地震记录，得到了残差更小的震源位置。

图2 定位的地震震源位置Fig.2 The hypocentral locations of earthquakes located.

将该算法用于2001年至汶川地震之前的四川地区1831个小震震源精确定位，结果更为精确地显示了地震的丛集性和在断裂带附近发生的特性，并且震源深度集中在浅部4～8km的范围内。这些结果为进一步研究汶川地震之前的地震活动特征及分析地震前兆现象打下了基础。

本文提出了考虑P波、S波和Pn波到时拾取误差来进行地震定位的方法。实际上地震波到时的读取和地震定位由不同的人员进行，仅根据地震观测报告无法得到确切的到时读取的误差。因此，本文采用直达P波到时的精度为0.5、S波到时拾取误差为1.5s、Pn波到时的误差为1.0s的原则来进行定位。更为严格的地震定位需要给定不同台站不同拾取震相不同到时的误差，这样才能更好地利用到时数据的拾取误差。因此，为了更好确定地震位置，建议在地震观测报告中增加震相拾取的误差参数。

图3 1831个地震的深度分布Fig.3 The depth distribution of the 1831 earthquakes.

我们编写的MATLAB程序可向感兴趣的读者提供。

致谢 审稿人对本文提出了建设性修改意见，特此致谢。

牟磊育，赵仲和，张伟，等.2006.用INGLADA与GEIGER方法实现近震精定位［J］.中国地震，22(3):294—302.

MOU Lei-yu，ZHAO Zhong-he，ZHANGWei，et al.2006.Using the INGLADA and GEIGER method to realize regional earthquake location accurately［J］.Earthquake Research in China，22(3):294—302(in Chinese).

万永革，李鸿吉.1995.遗传算法在确定震源位置中的应用［J］.地震地磁观测与研究，16(6):1—7.

WAN Yong-ge，LI Hong-ji.1995.The preliminary study on the seismic hypocenter location using genetic algorithms［J］.Seismological and Geomagnetic Observation and Research，16(6):1—7(in Chinese).

万永革，李清河，李鸿吉，等.1997.用遗传算法确定三维横向不均匀介质中的近震震源位置［J］.西北地震学报，19(2):7—14.

WAN Yong-ge，LIQing-he，LIHong-ji，et al.1997.Hypocentral locations of the near earthquakes in 3-D lateral heterogeneousmedium determined by the genetic algorithms［J］.Northwestern Seismological Journal，19(2):7—14(in Chinese).

万永革，王敏，沈正康，等.2004.利用GPS和水准测量资料反演2001年昆仑山口西8.1级地震的同震滑动分布［J］.地震地质，26(3):393—404.

WAN Yong-ge，WANG Min，SHEN Zheng-kang，et al.2004.Co-seismic slip distribution of the 2001 west of Kunlun Mountain Pass earthquake inverted by GPSand leveling data［J］.Seismology and Geology，26(3):393—404(in Chinese).

万永革，沈正康，刁桂苓，等.2008a.利用小震分布和区域应力场确定大震断层面参数方法及其在唐山地震序列中的应用［J］.地球物理学报，51(3):793—804.

WAN Yong-ge，SHEN Zheng-kang，DIAO Gui-ling，et al.2008a.An algorithm of fault parameter determination using distribution of small earthquakes and parameters of regional stress field and its application to Tangshan earthquake sequence［J］.Chinese JGeophys，51(3):569—583(in Chinese).

万永革，沈正康，王敏，等.2008b.根据GPS和InSAR数据反演2001年昆仑山口西地震同震破裂分布［J］.地球物理学报，51(4):1074—1084.

WAN Yong-ge，SHEN Zheng-kang，WANG Min，et al.，2008b.Coseismic slip distribution of the 2001 West Kunlun Mountain Pass earthquake constrained using GPSand InSAR data［J］.Chinese JGeophysics，51(4):753—764.

万永革，沈正康，盛书中，等.2009.2008年汶川大地震对周围断层的影响［J］.地震学报，31(2):128—139.

WAN Yong-ge，SHEN Zheng-kang，SHENG Shu-zhong，etal.2009.The influence of2008Wenchuan earthquake on surrounding faults［J］.Acta Seismologica Sinica，31(2):128—139(in Chinese).

王福昌，万永革，胡顺田.2008.粒子群算法在主震断层面参数估计中的应用［J］.地震研究，31(2):149—154.

WANG Fu-chang，WAN Yong-ge，HU Shun-tian.2008.Application of particle swarm optimization to the estimation of mainshock fault plane parameters［J］.Journal of Seismological Research，31(2):149—154(in Chinese).

徐果明，周惠兰.1982.地震学原理［M］.北京:科学出版社.

XU Guo-ming，ZHOU Hui-lan.1982.Principle of Seismology［M］.Science Press，Beijing(in Chinese).

赵珠，丁志峰，易桂喜，等.1994.西藏地震定位:一种使用单纯形优化的非线性方法［J］.地震学报，16(2):212—219.

ZHAO Zhu，DING Zhi-feng，YIGui-xi，et al.1994.Location of Tibetan earthquakes:A nonlinear approach by a simplex optimized technique［J］.Acta Seismologica Sinica，16(2):212—219(in Chinese).

赵珠，范军，郑斯华，等.1997.龙门山断裂带地壳速度结构和震源位置的精确修定［J］.地震学报，19(6):615—622.

ZHAO Zhu，FAN Jun，ZHENG Si-hua，et al.1997.Crustal structure and accurate hypocenter determination along the Longmenshan Fault zone［J］.Acta Seismologica Sinica，19(6):615—622(in Chinese).

赵仲和.1983.多重模型地震定位程序及其在北京台网的应用［J］.地震学报，5(2):242—254.

ZHONG Zhong-he.1983.An earthquake location program withmultiple velocitymodel and its application in the Beijing seismic network［J］.Acta Seismologica Sinica，5(2):242—254(in Chinese).

周仕勇，许忠淮，韩京，等.1999.主地震定位方法分析以及1997年新疆伽师震群高精度定位［J］.地震学报，21(3):258—265.

ZHOU Shi-yong，XU Zhong-huai，HAN Jing，et al.1999.Analysis on the master eventmethod and precise location of 1997 Jiashi strong earthquake swarm in western China［J］.Acta Seismologica Sinica，21(3):258—265(in Chinese).

Geiger L.1912.Probability method for the determination of earthquake epicenters from arrival time only ［J］.Bull St Louis Univ，8:60—71.

Jackson D D，Matsu'ura M.1985.A Bayesian approach to nonlinear inversion［J］.JGeophys Res，90(B1):581—591.

Lee W H K，Lahr JC.1975.A computer program for determining hypocenter，magnitude，and firstmotion pattern of local earthquakes［J］.U SGeol Surv Open-file Rep，75—311.

Mooney W D，Laske G，Masters G.1998.A global crustalmodel at5°x 5°［J］.JGeophys Re，103:727—747.

Shen Z K，Jackson D，Ge B.1996.Crustal deformation across and beyond the Los Angeles Basin from geodeticmeasurements［J］.JGeophys Res，101(B12):27957—27980.

SpenceW.1980.Relative epicenter determination using P-wave arrival-time differences［J］.Bull Seism Soc Amer，70(1):171—183.

Takagi A，Hasegawa A，Umino N.1977.Seismic activity in the Northeastern Janpan Arc［J］.JPhys Earth，25(Suppl):95—104.

Waldhauser F，Ellsworth W L.2000.A double-difference earthquake location algorithm:Method and application to the Northern Hayward Fault，California［J］.Bull Seism Soc Amer，90(6):1353—1368.

Wu Y M，Chang C H，Zhao L，et al.2008.A comprehensive relocation of earthquakes in Taiwan from 1991 to 2005［J］.Bull Seism Soc Amer，98(3):1471—1481.

EARTHQUAKE LOCATION METHODW ITH ARRIVAL TIME UNCERTAINTY CONSIDERED AND ITS APPLICATION TO LOCATION OF EARTHQUAKES FROM 2001 TO 2008 IN SICHUAN AREA

WAN Yong-ge1)SHENG Shu-zhong1)CHENGWan-zheng2)ZHANG Zhi-wei2)WU Ye1)ZHAO Xiao-yan1)BU Yu-fei3)XUE Zhi-fang4)LIU Jin-li5)

1)Institute of Disaster Prevention，Yanjiao，Hebei 065201，China

2)Sichuan Earthquake Administration，Chengdu 610041，China

3)Gansu Earthquake Administration，Lanzhou 730000，China

4)Qinhuangdao Center Station，Hebei Earthquake Administration，Qinhuangdao 066100，China

5)Quzhou Experimental High School，Hebei Province，Quzhou 057250，China

The precisely located earthquake catalogue is important to seismicity，seismic tomography and crustal stress inversion studies.It also has greatapplication value in rapid report of an earthquake that just occurred.By considering the arrival time uncertainty，and the constraints on station elevation and seismic depth，we propose a relatively accuratemethod to estimate hypocentral location and its uncertainty based on inversion theory.Ourmethod can combine the arrival times of Pg wave，Sg wave，Pn wave and Sn wave in hypocenter location，so it increases the location accuracy by involvingmore data;and it can be also used in local and regional earthquake location simultaneously.In order to test our location method，we located earthquakes by using the simulated data with different uncertainty of Pg，Sg，Pn，Sn arrivals.The result shows that the location determined by using ourmethod ismore accurate than that by using othermethod.We apply it to earthquakes occurring in the period from 2001 to 2008 in Sichuan area，and obtained amore clustered hypocentral distribution convergent to the fault zones.The result provides a solid foundation for studies of seismicity，geometry of the active faults and seismic tomography in Sichuan region.It is also helpful to study the seismicity precursors before theWenchuan earthquake.

earthquake location，uncertainty，station elevation，Sichuan area

P315.5

A

0253-4967(2012)01-0001-10

10.3969/j.issn.0253-4967.2012.01.001

2011-05-03收稿，2011-07-30改回。

国家自然科学基金(41074072、40874022)、中国地震局教师科研基金(20100101)、中央高校基本科研业务费专项(ZY20110101)和地震行业专项(200808053，201208009)共同资助。

万永革，男，1967年生，2001年在中国地震局地球物理研究所获得博士学位，2002—2006年在中国地震局地质研究所进行博士后合作研究，研究员，主要从事地震学和地球动力学方面的研究，电话:010-61597607;E-mail:wanyg217217@vip.sina.com;wanyongge8@gmail.com。