山东地区地震波一维速度模型研究

李 铂，崔 鑫，叶庆东，于 澄

（1.中国科学技术大学，合肥 230026；2.山东省地震局，济南 250014；3.中国地震局地球物理研究所，北京 100080）

0 引言

地震波速度结构的反演是地震学的经典问题，精准的速度模型是许多研究的基础性资料，如地震定位、大尺度构造应力场研究、地球动力学等。速度结构研究主要可分为一维速度结构和三维速度结构研究，一维速度结构经典的研究方法主要有Herglotz－Wiechert方 法、Gutenberg 反 演 方 法、τ（p）法［1］等，主要依据大量地震资料的走时曲线计算速度结构，目前应用较少。地震波速度结构研究有2个主要的难点，一是地震波速度结构和地震定位的耦合问题，两者会互相影响其误差和准确性；二是地震数据的结构问题，即地震台站深度普遍位于近地表，难于深入地下，从而难以在深度上进行有效约束。为了解决地震波速度结构和地震定位的耦合性问题，在Geiger定位方法的基础上，1976年，Crosson［2］首次提出了震源参数和速度结构联合反演方法，将速度结构作为未知参数与震源参数同时反演，由此弥补选取速度模型引起的误差，得到了广泛的应用。Aki［3－4］等人将该方法推广到三维速度结构，Pavlis［5］，Spencer［6］，刘福田［7－8］等人相继提出了参数分离法，大大提高了运算效率。Kissling［9］等于1994年提出可以基于走时残差均方根最小作为目标函数进行一维速度结构和震源位置联合反演，定位精度更高。

基于现在全球数字地震观测技术的发展和观测台网的大量建设，积累了大量高质量的地震波走时数据和数字波形记录，使得三维速度结构的反演方法取得了长足的进步。目前层析成像方法在数据处理时分为2组：波速分析和地震偏移。波速分析的目的是反演地球内部介质的波速结构，其主要方法有：体波走时层析成像［10］、面波走时层析成像［11］及接收函数法［12］。

1 区域地质构造特征

山东省位于我国地势划分中的第三大阶梯中［13］，海拔高度除一小部分山地超过千米以外，大部分山地丘陵都在500m 左右，地势起伏较小。山东省地层在全国地层区划中，属华北地层大区的晋冀鲁豫地层区和秦祁昆地层区，前者可分为华北平原、鲁西和鲁东3个地层分区，后者为祁连－北秦岭地层分区（鲁东南地层分区）。鲁西北－鲁西南平原区属华北平原地层分区，鲁中南山地丘陵区属鲁西地层分区，鲁东丘陵区的西北部（鲁东北地区）属鲁东地层分区，鲁东丘陵区的东南部（鲁东南地区）属祁连－北秦岭地层分区（鲁东南地层分区）。山东省各断代地层发育比较齐全，自中太古代至新生代地层都有分布，地表出露以中、新生代地层为主，其次为古生代地层，元古宙地层分布局限，太古宙地层零星出露。华北平原地层分区以发育大量新生代地层区别于鲁西地层分区，鲁东地层分区和祁连－北秦岭地层分区以没有古生代地层区别于华北平原地层分区和鲁西地层分区，鲁东地层分区、祁连－北秦岭地层分区与鲁西地层分区（包括华北平原地层分区）三者的前寒武纪基底岩系组成特征明显不同。

山东省表层构造首先反映在地貌特征上，是由鲁中和半岛地区的低山丘陵及环绕的堆积平原、陆架海域［14］构成的现代地貌格局。沂沭断裂带纵贯山东中部，如“刀劈斧砍”将山东一分为二。苏鲁造山带则奠定了鲁东地区基底构造线的总体格局。齐河－广饶断裂和聊城－兰考断裂则是分划鲁西地块和华北坳陷平原的构造带。因此，山东的地质块体所反映的地表构造格局具有一坳（济阳坳陷）、两块（鲁西地块、胶北地块）、两带（沂沭断裂带、苏鲁造山带）及一域（黄、渤海陆架海域）6大构造块体格局。

图1 山东地区断层分布图

2 理论与算法

地震波的到时是台站坐标（s）、震源参数（h，包括发震时刻和地理坐标）和速度场（m）的非线性函数，

一般而言，真正的震源参数和速度场都是未知的。仅有到时和台站坐标可以测量得到，我们不能直接解出（1）式。所以我们必须对未知参数做有根据的推测。使用初始速度模型，我们可以对初始震中到台站进行射线追踪，并计算理论到达时间（tcalc）.观测到时与理论到时的差，即为走时残差（tres），可以将其视为当前震中和速度模型与真实震中和速度模型的差（Δ）。为了计算震源和模型参数之间的偏差，我们需要知道观测走时对所有参数的相关性。对于（1）式进行泰勒展开，即可获得基于震中和速度模型的走时残差。

根据矩阵理论，震源参数和速度模型之间的耦合关系写为如下形式：

式中：t为走时残差向量；

H为走时对震源参数的偏导数矩阵；

h为震源参数的扰动量向量；

M为走时对速度模型的偏导数矩阵；

m为速度模型的扰动量向量；

e为走时误差向量，包括观测走时拾取误差以及由于台站坐标造成的理论走时误差、不适当的速度模型和震源参数的误差、方程线性化等的误差；

A为所有偏导数矩阵。A的结构为

d为震源参数和速度模型的扰动量，d＝（Δh1，Δh2，Δh3，Δh4，Δm1，…，Δmn）。

在定位地震时忽略公式（3）中Mm 的作用可能引入震源位置的系统误差［15－16］。与此类似，在等式中忽略Hh可能导致速度参数的误差［17－18］。因此，同时使用震源和模型参数并不一定能获得真正的震中和速度模型。

除非我们已经“猜”出了正确的震源坐标。地震数据和层析成像会更新震源和速度参数。Thurber［15］在1992 指出，联合反演震源－速度模型，要比轮换使用独立的震中和速度模型扰动可信度更高，即先用速度模型进行地震定位，再采用震源位置修正震源模型，反复迭代得到最终的震中和速度模型。Kissling［19］1988年在加利福尼亚长谷地区数据测试中展示了这种分离变量法在进行地震层析成像中的高效率性。

通常，公式（3）的解是个最小二乘解，该解使得误差（eTe）的加权组合最小，eTe即为整个算法的目标函数［20］。基于正态分布误差和扰动模型假设，最小二乘解能给出符合初始模型的可信度最好的解。通过搜索初始参考模型的邻域，初始模型的缺陷会导致三维层析成像的人工缺陷，特别是解不通过模型更新和反复迭代优化的时候。因为参数空间常常包含几千个未知数，简单的迭代过程可能不能避免局部最小，从而不能给出全局最小。而且A 的维数很大，分离变量法也可能会有困难。若有可能，计算模型、数据分辨率矩阵和模型协方差矩阵，能够对解了解更详细。否则震中参数和速度模型的权衡就很难。

由于山东地区地质构造比较复杂、地壳上的地幔速度结构区分明显［21］，所以如果采用简单的一维速度模型，实际速度模型差异对定位结果肯定会有较大的影响。利用Kissling方法计算得到一维速度模型，包括台站校正值，应用于地震定位中，很大程度上可消除初始模型与真实速度模型的差异，进而提高地震定位的精度。本文采用Kissling 计算方法，利用Kissling等人编写的程序Velest来计算山东地区P波的一维速度模型。

3 观测资料

2007年山东数字地震台网开始运行，截止目前，接入山东地震台网的地震台站为63个，其中宽频带地震仪39台，短周期地震仪24台，台站分布较均匀，平均台间距50km。

为了获得尽可能好的反演效果，我们通过以下标准进行地震事件的选择：

（1）震级大于ML2.0级；

（2）震中位于山东及其邻区，且震中位于陆上；

（3）地震定位精度为1。

同时，为了取得山东局部地质地区的地震波速度模型，我们针对3 个地区选取了部分地震事件。主要是：河南范县和山东鄄城交界地区，郯庐带两侧地区。河南范县和山东鄄城交界地区主要是以聊城－兰考断裂作为研究区域，郯庐带则以昌邑－大店断裂作为中心选取地震事件。

经过筛选 总计有45次地震事件，919条地震射线。震中分布图如图2所示，其中红色实心圆为地震震中，绿色三角形为地震台站。图中用颜色标示了山东及其邻区的高程信息，对照颜色与高程色块图，可以看出山东及其邻区大部分都是平原（紫色区域），中部有山脉分布（绿色区域）：

图2 震中与台站分布图

4 速度模型确定

Kissling方法基本计算过程为：在参考已有先验信息的基础上，选取不同的速度模型作为初始模型，同时反演震源参数和速度模型，包括台站校正值，再在初始模型的基础上对层厚和速度进行微调，经过反复计算得到走时均方根残差最小的一维速度模型作为我们最后求得的一维速度模型［22］。选取不同的速度模型反复计算旨在求得多值解空间中的全局最小速度模型，避免陷入局部最小。

山东及其邻区相关的速度模型主要有：赵仲和［23］1983年给出 的MDBJ模型，陈立华［24］1990年给出的华北地区速度模型（下述简称“陈立华模型”），以及IASP91全球速度模型，如图3所示：

图3 3种不同的P波速度模型

由于陈立华模型相对其他模型较新，而且使用了更多的地震数据，其他模型则或者是全球模型，较为粗糙，又或者研究对象为首都圈地区，因而山东地区的地震数据结构较差。基于以上考虑，我们最终使用陈立华模型作为初始模型进行反演。由于山东的地震事件深度普遍较浅，我们选取的事件中，最深的地震位于地下10km，这自然影响了地震速度模型的反演，因为我们只能得到地震射线最低点以上的速度结构，对于更深的地方，对速度模型没有分辨。为此，我们需要求得地震射线的最低点，以便参考确定速度模型的深度下限［25］。经计算，地震射线的平均深度为17.1km，我们以此作为参考，取略深于这一深度的速度模型，开始进行迭代计算。断层层数为8层，最下面一层为均匀半无限空间，反演过程不改变层厚度，只进行层中波速的计算。计算到同时满足下列3项条件下迭代中止：①震源位置，台站校正值和速度值已没有很大的变化；②所有地震均方根残差相对于第一次重新定位后的结果明显地减小；③计算的一维速度模型和台站校正值有实际的地质意义且没有违反先验信息。据此，我们最终计算出一维P波速度模型［26］，如图4中的黑色加粗实线，该速度模型从地表到30km 深度处分为5层。从图4中可以发现，该模型与陈立华模型在浅层近似程度较高，但是陈立华模型在21～26km 处有低速层，本文模型无法反映该低速层的存在，这可能源于山东地区的特殊地球物理背景，也有可能是因为该深度范围的地震样本数较少，无法给出较好的速度约束。初始的地震到时数据的残差均方根为1.08s，我们给出的震中和速度模型使残差均方根降低至0.62s。

图4 本文得到的P波速度模型

为了研究山东各地区的速度模型，我们选取了3个地区分别计算其对应的速度模型，分别为：鲁西模型（震中位于河南范县及其邻区），郯庐带西侧模型，郯庐带东侧模型。经计算可得到图5所示的速度模型。

图5 鲁西、郯庐带P波速度模型对比

如图5 所示，在地壳浅层，郯庐带西侧速度稍低，郯庐带东侧则相对具有稍高的速度，这可能与区域的地质背景有关，如郯庐带西侧地区的沉积层覆盖较厚，因而速度较低，东侧地区则相对较薄，基岩出露的情况也较多，因而速度较高。

鲁西地区相对于郯庐带两侧地区，其浅层速度较低，而在10km 深度左右速度则较高。浅层的低速可能源于地表较厚的沉积层，而10km 深度的速度差异则反映了鲁西和鲁东地区花岗岩层岩石的差异。

图6 原始震中与修订震中图

我们再将原始的地震震中位置和与最终速度模型联合给出的震中位置进行做图比较，如图6所示，红色圆圈为原始地震震中，橙色圆圈为最终速度模型对应的震中位置，绿色三角形为地震台站。可以发现，在郯庐断裂带两侧区域，新的震中位置相对于原始震中位置距离断层更近。

5 讨论与结论

（1）通过Kissling提出的计算最小一维速度模型方法，利用山东地震台网的观测走时资料，确定了山东及其邻区的P 波一维速度模型以及台站校正值。

（2）采用一维速度模型，对选取的地震重新定位。定位后的地震走时残差有一定的降低，地震定位精度在经度、纬度、深度方向上有了提高。表明了一维速度模型要优于其他速度模型。同时走时残差降低的幅度有限，说明用一维速度结构近似真实的地壳结构仍然过于简化，应该使用更复杂更合理的速度模型进行联合定位，才能期望得到更高精度的结果。

（3）地震深度和震相的选择限制了速度结构的深度。由于山东地区地震事件普遍位于近地表，深度普遍小于10km，而我们选取的是直达P震相，其射线的最低点往往不超过20km，因而无法获取更深处的地震波速度信息。在这个意义上，PmP 与Pn震相具有重要意义，因为这2个震相携带了地壳深部、莫霍面以及上地幔顶部的相关信息。未来的研究工作中应考虑加入这2个震相的到时数据。

［1］ Geiger，L.（1912）.Probability method for the determination of earthquake epicenters from arrival time only.Bull.St.Louis.Univ.，8，60－71.

［2］ Crosson，R.S.（1976）.Crustal structure modeling of earthquake data，1，Simultaneous least squares estimation of hypocenter and velocity parameters.Geophys.Res.，81（17）：3036－3046.

［3］ Aki，K.，Lee，W.H.K.（1976）.Determination of three－dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes，part 1：A homogeneous initial model.Geophys.Res.，81（23）：4381－4399.

［4］ Aki，K.，et al.（1977）.Determination of the three－dimensional seismic structure of the lithosphere.Geophys.Res.，82（2）：277－296.

［5］ Pavlis，G.，Booker，J.R.（1980）.The mixed discrete－continuous inverse problem：Application of the simultaneous determination of earthquake hypocenters and velocity structure.Geophys.Res.，85（B9）：4801－4810.

［6］ Spencer，C.，Gubbins，D.（1980）.Travel－time inversion for simultaneous earthquake location and velocity structure determination in laterally varying media.Geophys.J.Roy.Astr.Soc.，63（1）：95－116.

［7］ 刘福田.震源位置和速度结构的联合反演（I）——理论和方法［J］.地球物理学报，1984，27（2）：167－175.

［8］ 李强，刘福田.一种横向不均匀介质中地震基本参数的测定方法［J］.中国地震，1991，7（3）：54－63.

［9］ Kissling E，et al（1994）.Initial reference models in local earthquake tomography.J.Geophys.Res.，99，19 635－19 646.

［10］ 雷建设.浅析中国地区体波走时层析成像研究进展［J］.西北地震学报，2000，22（4）：471－475.

［11］ 黄忠贤.中国东部海域岩石圈结构面波层析成像［J］.地球物理学报，2009，52（3）：653－662.

［12］ 黄海波.利用远震接收函数方法研究南海西沙群岛下方地壳结构［J］.地球物理学报，2011，54（11）：2788－2798.

［13］ 蒋科植，朱介寿，程先琼.利用面波方法研究欧亚大陆及周边海域岩石圈结构［J］.华北地震科学，2012，30（2）：1－7.

［14］ 王帅军，王夫运，张成科，等.北京及附近地区地壳速度结构与构造［J］.华北地震科学，2011，29（2）：6－12.

［15］ Thurber，C.H.（1992）.Hypocenter－velocity structure coupling in local earthquake tomography.Phys.EarthPlanet.Inter.，75：55－62.

［16］ Eberhart－Phillips，D.，and A.J.Michael（1993）.Three－dimensional velocity structure，seismicity，and fault structure in the Parkfield region，central California.Geophys.Res.，98，737－758.

［17］ Michael，A.J.（1988）.Effects of three－dimensional velocity structure on the seismicity of the 1984 Morgan Hill，California，aftershock sequence.Bull.Seismol.Soc.Am.，78：1199－1221.

［18］ Vander Hilst，R.D.，and W.Spakman（1989）.Importance of the reference model in linearized tomography and images of subduction below the Caribbean plate.Geophys.Res.Lett.，16，1093－1096.

［19］ Kissling，E.，et a1（1995）.Improved seismic velocity reference model from local earthquake data in Northwestern Italy.Terra.Nova.，7，528－534.

［20］ Spakman，W.，and G.Nolet（1988）.Imaging algorithms，accuracy and resolution in delay time tomography in Mathematical Geophysics.In Vlaar N.J.，et al.（Eds.），Mathematicalgeophysics：asurveyofrecentdevelopmentsinseismologyandgeodynamics（pp.155－187）.Dordrecht：Reidel Publ.Co..

［21］ 于利民，刁桂苓，李钦祖.由深源远震体波记录反演华北北部地壳上地幔速度结构［J］.华北地震科学，1983，（3）：11－20.

［22］ 傅淑芳，刘宝诚.地震学教程［M］.北京：地震出版社，1991：447－480.

［23］ 赵仲和.北京地区地震参数与速度结构的联合测定［J］.地球物理学报，1983，26（2）：131－139.

［24］ 陈立华，宋仲和.华北地区地壳上地幔P波速度结构［J］.地球物理学报，1990，33（5）：540－547.

［25］ 万永革，李鸿吉.京津唐张地区速度结构和震源位置联合反演的遗传算法［J］.地震地磁观测与研究，1996，（3）：57－66.

［26］ 王小龙，刘渊源，余国政，等.利用远震接收函数反演乌江彭水电站地震台下方地壳厚度［J］.华北地震科学，2011，29（2）：13－18.