基于局部一维模型与走时标定的区域三维速度模型构建技术研究

张慧民 靳 平 刘文学 李 欣 张锐迎

1）中国西安710024西北核技术研究所

2）中国西安710049西安交通大学

3）中国西安710043北方光电股份有限公司

引言

准确定位对于《全面禁止核试验条约》（comprehensive nuclear test ban treaty，简写为CTBT）的地震监测工作具有重要意义.实践中面临的主要挑战之一是稀疏台网条件下对小震级地震事件进行精确定位.目前地震定位中采用的主要是基于IASP91（Kennett，Engdahl，1991）或AK135（Kennett et al，1995）一维地球模型的走时表.虽然这些一维模型很好地反映了全球的平均性质，但实际地球是三维的，穿过三维介质的实际震相走时与由一维模型导出的理论走时之间往往存在偏差，尤其是对区域震相，在很多地区这种偏差可能很严重.在我国的新疆地区，地震事件在部分区域台站上的实际走时与理论走时之间就存在着较为明显的偏差.这种区域震相走时的偏差可能会影响到地震台网的定位精度.为解决这一问题，目前通常采用两种方法：① 利用爆炸或精确定位的大地震进行区域走时标定；②构建监测区域的地壳与上地幔三维模型，并据此得到高精度区域震相走时.相比之下，第二种方法更能反映问题的物理实质，并具有更强的实用性.

为得到能够提高区域地震定位精度的三维速度模型，以往主要通过两种方式进行模型构建：一种方式是直接反演法，即将模型区域划分为三维网格并反演其速度结构，或反演模型区域内部分观测点下方的一维模型后进行插值得到三维模型；另一种方式是利用已有的地震学研究成果，结合该地区的地质构造等资料，以震中与发震时间已知的地震事件（ground truth，简写为GT）为参考，对已有的介质模型进行选择、组合及必要的修正，从而得到准确反映区域走时特点的三维介质模型（Russian Federation／United States Calibration Working Group，2002；Steck et al，2004；Cormier，2001；Ryaboy et al，2001）.对禁核试地震监测工作而言，构建区域性地壳与上地幔三维速度模型的主要目的在于提高区域地震定位能力.出于此目的，Ryaboy等（2001）采用第二种方式构建了北美、北欧地区的三维介质模型，有效减小了区域地震定位误差.其作法的核心为：根据标定事件绘制区域走时偏差图（即从由标定事件获取的走时图上减去由IASP91模型得出的走时图），将其与现有模型的理论走时偏差图进行对比，找出两者的主要差异，然后结合附近的一维模型、有关地质构造信息等，确定如何修正现有模型以减小两走时偏差图的差异，并据此得到新的模型.如此不断修正现有模型以减小理论走时偏差图与标定事件走时偏差图之间的差异，最终即可获得满意的区域三维模型.与直接反演法不同，这一构建方法工作量很大且对人员的专业经验要求较高，实现难度较大.为此我们对其进行了改进，以类似直接反演的方法为构建过程提供数值参考，从而降低了此类构建过程的难度并减小了最终模型对人员主观因素的依赖.

1 方法

在已往的此类模型构建方法中，如何对已有模型进行修正，主要是由构建人员根据走时偏差图和自身的经验确定的.在实践中，这种做法需要构建人员具有非常丰富的经验，而且其主观因素会对模型的形成产生比较大的影响，这使得模型构建工作实施起来比较困难，且结果可能不够客观.为此，本文提出了借助于体波走时反演方法，为构建过程中如何修正模型提供较为客观的、半定量的参考意见，据此协助构建人员实现区域模型的构建，以降低构建难度并减少对构建人员个人经验的依赖程度.与一般直接反演方法不同的是，我们只是利用GT事件和走时反演方法，为构建人员修正现有模型提供参考数据，而并非直接进行介质速度结构反演.

本文的三维模型构建工作主要按以下步骤进行：① 根据已有的模型数据和地质资料，给出区域初始模型；② 计算现有模型走时偏差图，并以走时反演法对该模型进行计算，记录此方法对现有模型的修正量；③观察走时偏差图，找出主要的走时偏差异常，对比走时反演法对现有模型的修正数据，找出其中与主要走时偏差异常相关的部分，根据现有资料判断是否需要对此修正数据再进行必要的校正（例如，反演法所得修正量过大，会导致模型与实际介质不符等），之后以校正后的修正数据对现有模型进行调整，得出一个新模型；④反复进行上述②、③两步骤，直至走时偏差图上不再有明显的异常；⑤以走时反演法对已无明显走时偏差异常的新模型进行微调，得到最终的区域三维介质模型；⑥ 检验所得三维模型在区域地震定位中的作用.

为实现上述模型构建步骤，需解决如下几个问题：

1）三维介质模型中的射线追踪.我们实现了二维格点化介质中计算初至波走时的初至波逐列递推算法（Schneider et al，1992），通过做剖面将其推广到三维介质的初至波走时计算中，并针对介质速度结构反演、地震定位等不同应用对其进行了相应的改进.实验表明，此算法计算速度快，且相对计算误差小于1‰，完全可以满足后续工作的需要.

2）理论走时偏差图计算.实践中直接逐点计算走时偏差图计算量较大.为此，我们结合数值插值设计了合理的格点计算顺序，使区域走时偏差图的计算量降低了约一个量级.

3）走时法介质三维速度模型反演.以射线追踪为基础，采用联合迭代重建方法进行.为提高射线追踪算法效率，本文设计了基于逐列递推算法的反方向射线追踪算法，即首先由逐列递推算法计算出地震射线在接收点的出射角；然后由接收点出发，以出射角的反方向作为入射角，计算出射线传播路径；之后，以此射线为参考，采用线性修正的方式修改反方向入射角，直至所得射线路径满足预设要求.

4）三维介质中的区域地震定位.本文采用Powell共轭方向算法进行.

2 数据

本文所用数据主要包括GT事件、部分台站下方的一维速度模型、相关地震地质资料等.

GT事件数据引自新疆维吾尔自治区地震局2006—2008年的地震公报.在其标称定位精度为1级（误差小于5km）的地震事件中，选取满足数据质量与台网评定质量均为1、定位台站数不少于8、台站最大张角小于90°、定位台站最小震中距小于80km等条件的286个地震事件作为GT事件.后面的工作中我们认为这些事件的震中参数是准确的.我们将挑选出的GT事件分为两组：其中一组包含146个事件（图1），含P波走时数据2 004条，S波走时数据1 883条，用于三维介质模型构建；另外一组的140个GT事件则作为震中参数已知的参考事件用于区域地震定位实验，以检验所构建三维模型用于定位时的效果.模型构建过程中使用的地震台站均为新疆地震台网台站.

图1 模型构建所用数据及台站分布图蓝色圆点表示GT事件，红色三角表示地震台站，黄色方块则表示该位置已有一维模型Fig.1 Data and stations used in the construction of 3D velocity model Blue dots denote GT events，red triangles represent stations，yellow squares stand for locations of 1Dpartial velocity model

已有的地震地质资料包括Crustal 2.0全球地壳模型数据，新疆地区莫霍面深度等高线图，以及新疆及周边地区介质结构的相关研究成果（胥颐等，2006；姜枚等，1999；李顺成等，2005；米宁等，2005；李海鸥等，2006；钱辉等，2006；高锐等，2001，2002；贺日政等，2001；赵俊猛等，2003，2008；冯梅，安美建，2007；李秋生等，2001；王有学等，2004；杨主恩等，2005；杨少敏等，2008；周仕勇，1999）等.

区域内的部分一维速度模型是通过接收函数与面波频散曲线联合反演方法获得的，其分布如图1所示.

3 模型构建

首先根据GT事件走时数据得到各个台站的实测走时偏差图，然后将新疆地区划分为1°×1°的网格，根据已有的部分一维速度模型生成初始三维介质模型.对于存在多个一维模型的网格，取其所有一维模型的均值；对于无一维模型与之对应的网格，则采用克里金空间插值方法获得.之后，计算现有模型走时偏差图，计算并记录体波走时反演法对现有模型的修正量；结合现有模型走时偏差图与实测走时偏差图的主要差异，找出相应的反投影法对现有模型进行修正的修正量，结合其它资料判断是否需要对这些修正数据再进行必要的校正；最后以校正后的修正数据对现有模型进行调整，得出一个新模型.反复进行上述步骤，直至所得三维模型的走时偏差图与实测走时偏差图之间不再有明显的不同.

图2给出了所构建三维模型在不同深度上的S波速度分布.

图2 所构建的三维速度模型切面图Fig.2 Profiles of modified 3Dvelocity model

4 定位实验

为检验所构建三维模型是否能够有效地减小区域震相走时残差并提高区域地震定位精度，我们选取了4种不同的介质模型以便进行对比：IASP91一维地球模型，Crustal 2.0模型，目标区域初始三维模型，修正后的三维介质模型.

重新定位后各事件震中相对于参考事件原震中位置的偏差如图3所示.图中各参考事件原震中位置均位于极坐标原点（未显示）.重新定位后的震中位置则如图3中小圆圈所示，即图中小圆圈的空间坐标代表各地震事件重新定位后的震中与参考事件震中的偏差.上述4种模型的平均定位偏差分别为5.0，5.8，5.1和3.8km，修正模型中的定位偏差明显小于其它模型.另外，从图3中不难发现，修正模型中定位偏差更加集中于原点且空间分布比较均匀，不像在Crustal 2.0模型中系统偏南，而IASP91模型中则系统偏东.多少有些让人意外的是，对于以局部一维模型为基础得到的初始三维模型，其定位偏差与一维的IASP91模型相仿，并未像预期的那样对减小定位偏差起到作用.

对于重新定位后各事件的震源深度偏差（相对于参考事件），对IASP91模型定位结果而言，几乎所有的定位深度都明显小于参考事件深度，而修正模型的定位深度与参考事件深度最为接近，且基本上不再有系统偏差，如图4所示.就重新定位后各事件的发震时间偏差而言，IASP91模型、Crustal 2.0模型与初始模型的时间偏差相当；而修正模型则小得多，与震源深度偏差类似.修正模型消除了其它几个模型中发震时间系统偏晚的现象.对定位实验中走时残差的统计表明，修正模型的走时残差分布明显优于其它几个模型.图5给出了初始模型与修正后模型的走时残差统计图.由图5可以看出，修正后模型的走时残差明显小于初始模型，它们更加对称地集中于零残差附近，消除了初始模型中震相走时系统偏小的情形.

图5 初始模型与修正后模型的走时残差统计图Fig.5 Residuals in initial 3Dvelocity model and modified 3Dvelocity model

为检验所构建模型在稀疏台网条件下对小震级事件的定位能力，我们进行了相应的模拟实验.所谓稀疏台网条件下的小震级事件定位，形式上基本可以理解为：记录到的地震事件的台站比较少，因此参与定位的台站就比较少，台站对事件的张角可能比较大，可能没有距震源很近的台站等.为此，我们依然采用前面实验中使用的参考事件，只是随机剔除其一部分记录台站，使其台站最大张角与最小震中距扩大.模拟实验表明，此时各模型的定位偏差均较前面实验有所扩大.相比而言，修正模型中扩大得最小，且其定位结果没有系统偏差，其走时残差的分布也较其它模型更为合理.

5 讨论与结论

综上所述，通过综合利用目标区域内的局部一维模型与相关资料，并结合走时标定方法，我们构建了新疆地区三维介质模型.用于地震定位时，该模型的走时残差明显小于初始模型，且所得结果与参考事件各参数的相似程度明显优于初始模型及全球平均模型.其中初始模型中对各参考事件的定位偏差平均为5.1km，而修正后模型的定位偏差平均为3.8km，即从统计意义上讲，修正后模型的定位偏差较初始模型可减少约20%.这说明本文的构建方法是行之有效的，确能达到预期目的——通过构建区域三维介质模型，减小区域震相走时残差，提高区域地震定位精度.

令人有些意外的是，初始三维模型中的定位偏差与一维IASP91模型相仿，并未像预期的那样对减小定位偏差起到明显作用.但结合定位深度偏差与发震时刻偏差来看，初始三维模型仍有明显改善.如图4所示，其定位深度偏差明显要比IASP91模型小很多，而且消除了系统偏差，对于发震时刻偏差也有类似结果.从图5可以看出，初始三维模型中的走时有明显的系统偏差，这应该是造成其震中位置偏差较大的主要原因.

另外，从原理上讲，三维模型地震定位只需将普通定位方法中的走时计算部分的一维模型更换为三维模型即可，但实践中直接计算三维模型中的走时十分耗时，实用性很差，应以查表法代替.

在模型构建的过程中，构建人员综合其它相关资料，增加了所构建模型的合理性.但这一过程带给构建人员的工作量相当大，尽管我们引入的对模型修正的数值参考已经大大降低了这一工作量.这也是这一方法仍有待继续改进的方面.

冯梅，安美建.2007.中国大陆中上地壳剪切波速结构［J］.地震学报，29（4）：337－－347.

高锐，肖序常，刘训，管烨，李秋生，卢德源，李朋武.2001.新疆地学断面深地震反射剖面揭示的西昆仑—塔里木结合带岩石圈细结构［J］.地球学报，22（6）：547－－552.

高锐，肖序常，高弘.2002.西昆仑—塔里木—天山岩石圈深地震探测综述［J］.地质通报，21（1）：11－－18.

贺日政，高锐，李秋生，管烨，李朋武.2001.新疆天山（独山子）—西昆仑（泉水沟）地学断面地震与重力联合反演地壳构造特征［J］.地球学报，22（6）：553－－558.

姜枚，许志琴，薛光琦，史大年.1999.青海茫崖—新疆若羌地震探测剖面及其深部构造的研究［J］.地质学报，73（2）：153－－161.

李海鸥，姜枚，王亚军，张立树，于更新.2006.新疆富蕴—库尔勒剖面接收函数方法获得的地壳上地幔结构成像［J］.地质学报，80（1）：135－－141.

李秋生，卢德源，高锐.2001.新疆地学断面（泉水沟—独山子）深地震探测成果综合研究［J］.地球学报，22（6）：534－－540.

李顺成，刘启元，陈九辉，郭飙，赖院根，王继.2005.横跨天山的宽频带流动地震台阵观测［J］.地球物理学进展，20（4）：955－－960.

米宁，王良书，李华，徐鸣洁，李成，张勇，陈运平，于大勇.2005.天山和塔里木盆地结合部地壳上地幔速度结构［J］.科学通报，50（4）：363－－368.

钱辉，许志琴，姜枚，宿和平.2006.西昆仑接收函数反演与构造解析［J］.中国地质，33（2）：309－－316.

王有学，韩果花，姜枚.2004.阿尔泰—阿尔金地学断面地壳结构［J］.地球物理学报，47（2）：46－－54.

胥颐，刘建华，刘福田.2006.天山—帕米尔结合带的地壳速度结构及地震活动研究［J］.地球物理学报，49（6）：1693－－1700.

杨少敏，李杰，王琪.2008.GPS研究天山现今变形与断层活动［J］.中国科学：D辑，38（7）：872－－880.

杨主恩，张先康，赵瑞武，周伟新.2005.天山中段的深浅构造［J］.地震地质，27（1）：11－－19.

赵俊猛，李植纯，马宗晋.2003.天山分段性的地球物理分析［J］.地学前沿，10（S1）：125－－131.

赵俊猛，程宏岗，裴顺平，刘宏兵，张健狮，刘宝丰.2008.塔里木盆地北缘的深部结构［J］.科学通报，53（8）：946－－955.

周仕勇.1999.主地震定位法分析以及1997年新疆伽师强震群高精度定位［J］.地震学报，21（3）：43－－51.

Cormier V F.2001.Construction of 3-D earth models for station specific path corrections by dynamic ray tracing［C］∥Proceedings of the 23rdSeismic Research Review.Tucson，Arizona，USA，NNSA，1：30－－36.

Kennett B L N，Engdahl E R.1991.Travel times for global earthquake location and phase identification［J］.Geophys J Int，105（3）：429－－465.

Kennett B L N，Engdahl E R，Buland R.1995.Constraints on seismic velocities in the Earth from travel times［J］.Geophys J Int，122（1）：108－－124.

Russian Federation／United States Calibration Working Group.2002.Refinement of regional seismic event location in northern Eurasia using 3-D crustal and upper mantle velocity model［C］∥Proceedings of the 24thSeismic Research Review.Tucson，Arizona，USA，NNSA，2：176－－185.

Ryaboy V，Baumgardt D R，Firbas P，Dainty A M.2001.Application of 3-D crustal and upper mantle velocity model of North America for location of regional seismic events［J］.Pure Appl Geophys，158（1）：79－－103.

Schneider W A，Ranzinger K A，Balch A H，Kruse C.1992.A dynamic programming approach to first arrival traveltime computation in media with arbitrarily distributed velocities［J］.Geophysics，57（1）：39－－50.

Steck L K，Rowe C A，Begnaud M L，Phillips W S，Gee V L，Velasco A A.2004.Advancing seismic event location through difference constraints and three-dimensional models［C］∥Proceedings of the 26thSeismic Research Review.Tucson，Arizona，USA，NNSA，3：346－－355.