南北地震带南段地壳厚度重震联合最优化反演

陈石, 郑秋月, 徐伟民

中国地震局地球物理研究所， 北京 100081



重震联合反演； 地壳厚度； 多约束反演； 异常分离； 界面反演； 最优解

南北地震带南段地壳厚度重震联合最优化反演

陈石, 郑秋月*, 徐伟民

中国地震局地球物理研究所， 北京 100081

重力反演方法是研究地壳结构和物性界面起伏的有效地球物理手段之一.本文收集了南北地震带南段67个已有的固定台站接收函数反演的Moho面深度结果，并使用基于EGM2008重力异常模型计算的布格重力异常，验证了本文提出的重震联合密度界面反演方法的有效性.利用接收函数对台站下方Moho面深度估计作为先验约束，定义了一类评价函数，通过对重力反演算法中尺度因子，平移因子和稳定性因子的最优选择，最小化重力反演结果与接收函数模型之间的差异.结果表明，本文提出的方法，可以有效地同化不同地球物理方法获得的反演模型，且通过重震联合反演可以改进由于对空间分布不均匀的接收函数结果插值可能而引起的误差.本文还通过引入Crust1.0的Moho面深度为初值，同时考虑地壳密度的横向不均匀分布，通过模型之间的联合反演有效改善了地球物理反演模型间的不一致性问题.本文反演得到的最优化Moho面深度模型与已知67个台站位置接收函数模型之间的标准差约1.9 km，小于Crust1.0与接收函数结果模型之间标准差为3.73 km的统计结果.本文研究结果对于同化重震反演结果、精化地壳密度界面模型，都具有十分重要的参考意义.

1 引言

地球物理反演方法是了解地球内部结构的重要手段，Moho面起伏的确定是构建地壳模型的基础(地壳厚度等于Moho面深度与地形高程之和).采用布格重力异常和地震接收函数法，作为两种独立的地球物理手段，都可以用于研究地壳结构和确定Moho面深度特征(Langston, 1977;Guo et al., 2012; Meng et al., 2012).根据现有的研究结果表明，由于地球物理反演的多解性，用这两种方法得到的地壳厚度分布具有一定差异，这些差异性与反演的数据、方法都有着密切的关系(Li et al., 2014).一般认为，接收函数法对于台站下方的地壳速度变化界面或速度间断面较敏感，通常可以用于研究地壳结构(Langston, 1977).而布格重力异常也包含了各种壳内密度界面或不均匀体产生的异常信息，通过合适的异常分离方法，可以获取由Moho面起伏引起的异常最优频段，并反演地壳厚度模型.但是，接收函数法受地震台站分布位置的限制，反演结果的空间分布不均匀，而通过插值计算得到的地壳厚度分布在没有台站覆盖的位置，可靠性不强.相比之下，重力反演方法在数据上可以保证更好的空间覆盖性，近年来特别是基于卫星和陆地重力数据联合解算的EGM2008全球重力异常模型(Pavlis et al., 2008, 2012)，其精度和空间分辨率都有很高的进步(Zhang et al., 2009；付广裕等, 2013;)，且在地壳结构解释中得到广泛应用(陈石等, 2011a， 2011b， 2014，2015; Xu et al., 2014)，特别是对于一些地形起伏剧烈和地震台站分布稀疏地区，该模型也提供了较好的信息覆盖.但由于重力异常反映的是一种相对变化，且包含深浅部各种密度异常体的信息，因此，反演结果的精度有限，其绝对参考意义不强.

目前，接收函数方法已被广泛应用于研究台站下方的地壳结构，是获取岩石圈上地幔速度间断面深度的有效手段之一.接收函数可以被认为是除去震源、地震波传播路径和仪器响应等因素后的时间序列，包含了地震台站下方地壳和岩石圈上地幔范围内各个速度间断面所产生的转换波及多次反射波的信息(Xu and Zhao, 2008).从远震体波中提取的P波接收函数和S波接收函数可以广泛用于研究从地壳至地幔各种深度的壳内速度间断面和结构，并根据各种震相到时特点，形成优势互补.特别是随着流动宽频带地震仪的广泛应用，已经有很多研究结果可以作为约束性条件，为更精细地壳厚度模型的构建和修正提供参考(Xu et al., 2007；Chen et al., 2010; Wang et al., 2010; Li et al., 2014;).

但是由于接收函数计算本身也是一类典型的非线性反演问题，在反演过程中，受地壳结构的复杂性，滤波参数选择，一定时间段内地震事件的选取，H-k叠加过程中参数选择等多种因素的影响.对于同一个区域，不同来源观测数据采用相同的接收函数方法计算得到的地壳厚度差异通常也在0～4 km误差范围(Xu et al., 2007; Zheng et al., 2013；Li et al., 2014).

而重力反演方法主要是应用布格重力异常数据来反演Moho面深度.其中，除了考虑不同来源的重力异常模型观测误差可能引起的反演结果具有不一致性之外，更重要的是要考虑如何对布格重力异常进行分离，得到Moho或地壳厚度变化产生的重力异常成份.常用的重力异常分离方法，一般都是通过异常的频率特性而设计的，但是其实没有严格的物理定义场源深度和异常频率之间的关系.而这种分离异常的不准确性，是制约提高重力反演地壳厚度精度的主要障碍之一.本文的方法原理将进一步讨论这个问题，并提出解决方案.

接收函数法和重力反演方法对建立地壳厚度模型问题具有各自的优缺点.而本文的研究思路是通过重震联合反演方式，采用地震接收函数反演的多点地壳厚度结果作为已知条件，优化重力反演中不确定性参数的选取.主要对重力异常分离、反演尺度因子和平移因子等参数进行优化选择，将重力反演结果和已知多个约束点位置点信息的最小二乘解作为评价函数，通过利用两种方法的各自优点，减少地球物理反演多解性影响，反演得到最优化地壳厚度模型.

在研究区域位置上，我们选择了南北地震带南段地区为实验场，主要是因为该区地震活动性强，地震台站较多，且已经积累了较多的数据资料和研究成果，适合于对本文提出方法进行有效性验证.

2 方法原理

在地壳和地幔的分界面周围，由于介质物性存在较大差异，因此地震波速度会发生突变，这种在速度上的界面，同样也是一个密度差异性界面，通常称为Moho面.在重力反演中采用一个等效的密度界面来估算Moho面起伏引起的重力异常.而计算这个密度界面一般采用频率域快速正演算法(Parker,1973).该方法将重力位场积分核函数在频率域内通过Taylor公式展开，建立重力异常谱与不同阶次起伏密度界面谱之间的级数关系，如公式(1)所示.

(1)

其中，G为万有引力常数，Δρ=Δρ(x,y)为界面密度差(如果密度差为常数可以将其提出到括号外)，h0为起伏密度界面的平均深度，界面起伏表示为Δh=Δh(x,y)=h(x,y)-h0,ω为波数，F[*]为Fourier谱.对式中F[Δg]进行快速Fourier逆变换就可以得到空间域一个密度界面起伏的重力异常.一般由密度界面起伏计算重力异常称为正演问题，而由重力异常反演密度界面起伏称为反演问题(Oldenburg, 1974)，如果密度为常数，异常谱和界面谱位置交换表示为式(2)：

(2)

而反演公式(2)可以进一步表示成如(3)和(4)式所示的迭代形式，其中，第k步正演的异常谱表示为Fk[Δgc]，每一步计算得到(3)式的修正界面起伏，采用公式(4)迭代，直至算法收敛，完成反演：

(3)

Δhk+1=Δhk+δhk.

(4)

对于该方法的研究和应用很多(Gomez-Ortizetal., 2005;Shinetal., 2006)，通常正演计算较稳定，反演计算由于高频放大作用，通常需要引入低通滤波器和稳定性估计因子，辅助完成反演(Oldenburg, 1974).本文不具体讨论如何有效地完成反演等技术问题，而是讨论实际重力反演过程中，那些参数影响并如何影响反演结果.

2.1 重力异常的成份

重力异常表现为不同深度各种场源密度体信号的叠加形式.一般在反演过程中使用的布格重力异常可以看作地表以下偏离正常密度结构各个异常体效应的总和，如图1所示，其中这些场源几何形态、与正常结构之间的密度差、埋藏深度都不相同.在研究不同问题之前，有必要采用一定方法手段进行异常分离.由于重力异常的大小是以场源距离平方形式衰减，因此，一般认为距离近地表观测深度的密度异常体产生的异常具有高频特性，而深源位置产生的重力异常具有低频特性.所以，一般采用低通滤波器或者等效形式的方法进行重力异常分离.但是，各种研究都没有给出重力场源深度与异常频率之间的理论关系.

另外，由于重力异常校正中，采用的参考密度模型也可能与实际模型不符，主要是对于如图1所示的水平均匀密度层(Bouguer plate)，在地表产生的重力异常为一个常数，因此，不同观测得到的布格重力异常的绝对数值之间相差的常数项并没有意义，也可以理解为重力异常并不能分辨垂向水平层状异常体形式的密度结构.而在重力反演之前，一般也要对布格异常进行处理得到零平均形式的重力异常.

通过上面对重力异常组成成份的描述，可以对用于重力反演使用的布格异常有一定认识.首先，如果要使反演结果准确，对布格重力异常进行场分离是十分必要的，但是因为没有场源和异常特征之间物理关系，从而目前也没有统一的场源分离标准，因此，对布格重力异常的场源分离是第一个不确定性因素.其次，由于布格重力异常只是相对变化值有意义，也就是说对横向密度变化敏感，具有无法分辨垂向密度变化的局限性，所以，如果用于密度界面反演，会产生由于估计的界面起伏平均深度不准确而产生多解性，或以一个系统性误差形式体现.所以，我们需要进一步结合重力异常和反演算法，来讨论密度界面反演中的不确定性问题.

2.2 密度界面反演的不确定性

(5)

(6)

其中，对于给定密度差ρ不同，在反演结果中体现密度界面起伏程度差异，在这里称为“尺度因子”(Scalefactor)；而gDC值的选择，体现在界面反演结果中的平均深度不同，可称为“平移因子”(Shiftfactor)；而平均起伏h0的作用，与gDC等价；但是选择h0的数值大小要与重力异常网格间距相匹配，主要影响反演计算的稳定性，称为“稳定因子”(Stabilizedfactor).首先通过一组理论数据，模拟由于选择不同ρ和gDC对于反演结果的影响.

图2是采用相同的重力异常剖面数据，采取不同参数的反演结果.在图2中选取了10组不同参数模型，通过对比说明，上述尺度因子和平移因子参数选择不同对反演结果的影响.图2a和2b所示的反演结果参数不同，但反演过程中都满足收敛条件，参数变化范围如图例所示，结果除了可以说明重力反演结果的非唯一性特点，还可以说明不同模型参数选择对反演结果的规律性影响特征.在没有约束条件的情况下，尺度因子和平移因子的选择存在较大不确定性，而如果通过一些已知约束点信息，在合理的范围内优选两个因子参数，就可以提高不同地球物理资料反演结果的一致性程度，减小不同方法手段之间结果的差异，实现多场反演数据相互同化的目的.而这种“同化”的过程，可以相互利用各自独立地球物理场的优点，达到联合反演的目的.

图1 影响重力反演密度界面结果的因素示意图Fig.1 Sketch map on interface inversion using gravity anomaly and influence factors

图2 重力反演密度界面的多解性参数

2.3 重力反演中的高频放大因子

在频率域密度界面反演过程中，平均界面起伏h0的选择通常会影响反演的稳定性，同一个重力异常可能由于选择h0的不合适而影响收敛性.其主要由于eωh0项对异常谱产生高频放大影响.通常解决这种高频放大影响是通过施加一个合适的余弦低通滤波器来实现(Oldenburg, 1974)，但是有时尽管使用低通滤波器仍然可能由于h0选择不合适，使反演过程发散.在2.2节我们将h0的选择成为稳定因子，其中e指数项部分的最大值，可以由(7)式估算：

(7)

式中，Δxy=min(Δx,Δy)，Δx和Δy分别为重力异常网格数据在x和y方向上的采样间距.从式(7)中，我们可以看出所谓高频放到因子对反演稳定性的影响受到h0和Δxy之间的比值共同影响，如果h0远大于Δxy，这个高频放大作用将以指数形式增强，而如果选取与重力异常网格间距合适的平均深度h0会提高反演的收敛速度.

3 南北地震带南段重力异常和地壳厚度模型对比

南北地震带南段地区构造复杂、地形起伏大、强震活动剧烈，如图3所示.图中给出了1900年以来该区发生的6.0级以上地震分布，由于该区域受喜马拉雅东构造节的动力作用且作为青藏高原物质运移的通道之一，构造运动十分活跃.这使得该区的Moho面深度起伏很大，位于西北部的青藏高原地区Moho面深度可以超过60km，而在昆明东南部地区的Moho面深度小于30km.

本文首先基于Crust1.0给出的1°×1°分辨率的Moho面模型，提取了本文研究区范围内的Moho面深度数据(图4).在图4中也将本文收集的南北地震带南段范围内，使用远震接收函数方法获得的67个地震台站位置做了标注.图5给出了台站位置接收函数与Crust1.0两个数据源(Bassinetal., 2000;Laskeetal., 2013) 的Moho面深度对比(图5).图中用黑色虚线标注了±10km误差范围的区域.从这些点的对比结果看，Crust1.0的1°×1°空间尺度地壳厚度模型，是一种平均意义下的Moho面深度，相比接收函数的结果反演的地壳厚度结果存在较大的差异，最大差异约10km，标准差为3.73km.因此，如果直接用Crust1.0模型的地壳厚度去估算或填补接收函数台站位置没有覆盖区域的地壳厚度，也将会带来较大的系统误差.

图3 研究区位置、地形、断裂分布和地震活动图中背景为地形起伏、红色实线为断裂带位置、黄色圆圈为1900年以来6.0级以上地震位置、蓝色三角为接收函数台站位置；附图中桔黄色实线为活动地块边界；WCB为西域地块，TPB为青藏地块，SCB为华南地块，方框为研究区位置.Fig.3 The topography, faults and seismicity of studying area The background is the topography and relief. The red solid lines are faults. The yellow cicles are the epicenters of M6.0+ events. The blue triangles show the seiemic stations for receiver function inversion. The inset of lower left includes the active block boundaries using orange solid lines. WCB means West China Block. TPB means Tibetan Plateau Block. SCB means South China Block.

图4 研究区Crust1.0模型的Moho面深度分布和接收函数台站位置图蓝色三角形表示接收函数台站位置，黑色实线为已知断裂构造位置.Fig.4 The Moho depth of Crust1.0 model and the locations of seismic station The blue triangles show the locations of seiemic station for receiver function inversion. The black solid lines are the known tectonic and faults.

图5 接收函数台站位置不同模型间的Moho面深度差异对比Fig.5 The Moho depths at the locations of seismic station derived by Crust1.0 and received function model

为了解决以上问题，我们采用重力反演方法联合地震学给出的模型结果，改进Moho面深度的反演问题.Moho面深度的反演，在重力反演算法中可以等同于一个密度起伏界面看待，在数据选择方面一般需要使用布格重力异常.在重力异常模型选择方面，我们使用EGM2008全球重力模型数据，该模型提供的球谐函数高达2160阶，提供的重力信号空间分辨率能力约为5′×5′，可以给出波长约10km以上的重力信号(Pavlisetal., 2012;Hirtetal., 2013)，采用该模型为基础的各种扩展重力模型也在不断发展(Sandwelletal., 1997;Smithetal., 1997;Sandwelletal., 2009;Hirt, 2013;).在计算布格重力异常的过程中，首先使用Lambert坐标投影算法，将经纬度转换为大地坐标；然后，选择了FA2BOUG标准程序(Fulleaetal., 2008)用于地形校正和相关重力改正，得到布格重力异常.另外，由于Moho面起伏引起的重力异常主要以低频分量体现，我们还采取高斯低通滤波器用于区域重力异常分离，图6是我们计算得到的区域布格重力异常.异常在整个研究区内都呈现负值特征，最低异常高达-600mGal，青藏高原东南缘的异常值最低.通过与其他来源的重力异常进行对比，其低频分量部分，都具有较好的一致性，我们认为该数据对于反演本文研究的Moho面深度问题，完全可以满足异常精度要求.

4 重震联合反演结果

应用图6给出的区域布格重力异常，我们对重力异常数据的处理和分析，应用本文在第2节所述的方法原理进行了反演.反演过程中重力异常网格单元间隔统一采样为20km，界面起伏平均深度设置为46km(参考67个固定地震台站接收函数地壳厚度平均值，网格单元的空间间隔与平均深度比值根据(7)式原理不应太大，否则不利于反演收敛).为了减小快速Fourier变换引起的边缘效应，我们对异常网格进行了最小曲率扩边处理，扩边后的反演网格尺度为256×256.

如果采用常密度反演，得到的反演结果不能优于实际接收函数与已知Crust1.0模型之间标准差(图5)，因此，我们认为如果不考虑地壳密度横向不均匀分布是不够的.特别对于本文选择的南北地震带南段地区，由于构造动力学作用十分复杂，地壳密度的横向不均匀性较强，通过常规的异常平滑与分离技术，不能完全实现由Moho面起伏引起的重力异常分离.特别要获得重震一致性较好的高精度反演结果，地壳横向密度的不均匀分布是要特殊考虑的因素.因此，我们先应用Crust1.0模型给出的Moho面起伏，并通过变换公式(1)，将密度界面起伏作为已知量，来反演地壳的横向密度变化(图7).

图7给出的结果反映的是地壳平均的横向密度变化，仅具有相对意义，结果表明，在青藏高原东南缘即图中西北部地区密度最低，昆明西北部的攀枝花地区存在一个地壳密度相对高值区，喜马拉雅东构造节和西北部的四川盆地位置都是相对地壳平均密度高值区，地壳密度的横向密度差异性可以超过0.2g·cm-3.图7给出的地壳横向密度差异特征，主要依赖于Crust1.0的Moho面深度结果，而该模型给出的Moho面深度的空间分辨率为1°×1°，因此，应用该结果正演则可以得到一个以低频异常为主的区域布格重力场，在此基础上与图6重力异常相减，就可以得到一个剩余扰动重力异常，将其作为重力反演与Crust1.0平均Moho面起伏模型之间的修正量.对其采用本文第2节的方法进行反演，经过多次迭代，并用公式(6)该选择最优结果，最终得到图8的Moho面深度结果.

在图8中，将每个地震台站位置的重力反演Moho面深度结果和接收函数结果之差进行了对比，按照差异性大小共分成三类，用蓝、绿和黄三种颜色表示，不同颜色之间图标空间分布并没有出现典型规律特征.进一步采用图9来分别对比几个模型之间的差异.从图9a结果可以看出，在67个地震台站位置的Moho面深度结果，通过公式(6)优选的最优重力反演结果与接收函数结果具有较好的线性相关性，且优于图9b的结果.图9a给出的重震联合反演结果，模型间的Moho面深度差异优于±5km，优于图5所示的Crust1.0结果与实际接收函数的差异性对比结果.

图9模型的详细对比统计结果见表1.重力反演结果与实际接收函数反演结果之间差异性最小，均方差约为1.9km，模型间差异性的最大、最小值也优于表1中的前两组结果.

因此，本文提出的重震联合最优化反演方法，可以提供一个区域尺度内优于已知接收函数插值和Crust1.0模型的精化Moho面深度模型.本文的联合反演结果，通过引入已知接收函数结果作为评价函数，可以对重力界面反演算法中的非确定性参数进行优选，实现两种地球物理手段之间的数据同化，给出的结果具有更好的一致性.

图6 基于EGM2008重力异常模型计算得到的研究区布格异常布格重力异常计算采用FA2BOUG程序，异常分离采用50 km高斯滤波，网格间距重采样为20 km.Fig.6 The Bouguer gravity anomaly derived by EGM2008 modelBouguer gravity anomaly is computed by FA2BOUG program. The Guass filter with 50 km is used for the anomaly sepeartion. The re-sampled distance is 20 km.

图7 基于Crust1.0模型反演的地壳横向密度差异图中黑色实线为在该区域内的断裂构造带分布，蓝色三角形为接收函数台站位置.Fig.7 The horzontial density variations on the basis of Crust1.0 model The black solid lines are the known tectonic and faults. The blue triangles show the locations of seiemic station for receiver function inversion.

图8 重震联合最优化Moho面深度反演结果彩色圆圈为接收函数台站位置，不同颜色表示联合优化方法与接收函数得到的Moho面深度结果之间差异.Fig.8 The Moho depth derived by the joint approach of gravity-seismic The color circles show the locations of seiemic station for receiver function inversion. The different colors mean the diffence of Moho depth between the joint approach and receiver functions result.

表1 反演结果对比分析Table 1 Analysis on inversion

5 结论与讨论

根据本文前述的重力反演密度界面起伏的方法原理，我们从原理上阐述了重力异常反演多解性问题产生的三个主要来源，即异常分离的不确定性，平移因子与尺度因子对反演结果形态的影响.重力反演与接收函数反演作为两种相互独立的地球物理观测手段，对于地壳厚度模型的估计应该具有一致性和差异性，如果能联合两种方法各自的优点，即重力异常的空间覆盖率高，接收函数对台站位置下方的估算准确，那么就可以优化反演模型，实现联合反

演，达到数据模型之间相互“同化”的目的.

本文的研究结果表明，这种重震联合方法采用接收函数结果作为空间多约束信息来评价重力反演结果合理性，可以优化重力反演过程中不确定性参数的选择，减少反演结果的非唯一性.最后从不同模型之间的结果对比表明，重震联合反演结果误差最小(图9)，该方法可以有效减小不同模型之间的差异.综上所述，本文得出以下研究结论：

1)重力反演密度界面的多解性主要体现在尺度因子与平移因子的选择问题上，不同参数组合都可以得到完全相同的与观测异常拟合一致的模型结果.通过多约束点方法构建重力反演模型评价函数，可以达到优选反演模型的目的.

2)重力异常分离并不能完全分离不同埋深的场源异常，在构造复杂地区，利用Crust1.0模型提供的Moho面深度模型，可以用于初步估计地壳内部物质密度的横向不均匀性分布程度.

3)重震联合反演模型相比其他两个模型之间的相互对比，结果显示重震联合方法反演得到的结果更好(表1)，这主要是由于联合模型具有更好的同化效应，这有利于消除各种地球物理反演手段之间的误差影响.

4)通过对接收函数、Crust1.0和重力最优化反演三种模型结果之间的对比表明，本文的联合反演方法可以优化模型质量，但是仍存在一定的相互差异，这主要与EGM2008重力模型的精度，浅部物质的不均匀分布和地壳的横向密度变化估算精度有限等相关.因此，仍然存在进一步改进反演结果的可能.

在地球物理资料不断积累的同时，联合不同地球物理反演资料优化反演模型是一项重要的科学研究任务.本文着重讨论了应用重力异常反演密度界面过程中遇到的问题和多解性产生因素，及如何从算法原理上改善反演过程的稳定性入手.选择南北地震带南段作为研究区域，主要是由于该地区已有研究结果丰富，特别是有较多的地震台站观测结果，这对于研究重震联合反演方法是非常有帮助的.本文方法也可以进一步推广到其他区域开展研究.特别是如果对范围更大的区域开展反演工作，应该进一步考虑应用变密度模型参与反演.

另外，通过本文提出的方法可以不断地改善不同空间尺度的地壳厚度模型估计，也可以有效地改进我们对均衡重力异常计算的方法，更好地分离Moho面异常与壳内不均匀体异常，这对于我们以后研究地壳内部三维密度结构，发现潜在的闭锁单元，反演壳内康拉德界面起伏特征及研究各种壳内变形界面形态等问题都具有重要的参考意义.

在后续研究中，将尝试参考其他地球物理资料，进一步对壳内浅部物质的不均匀分布和地壳的横向密度变化给出定量约束，并尝试联合Crust1.0以外的其他地壳模型参与“同化”反演过程和对比检验，为本文反演方法中的参数优化选择和检验多手段地球物理反演结果的合理性提供技术方案.

图9 多种方法反演的Moho面深度结果对比

致谢 本文研究过程中得到了长安大学王万银教授在频率域反演方法和程序编制等方面的指导；中国地震局地球物理研究所李永华研究员提供了接收函数反演结果，在此一并致谢.

Bassin C, Laske G, Masters G. 2000. The current limits of resolution for surface wave tomography in North America.EOSTransAGU, 81: F897.

Chen S, Wang Q S, Xu W M, et al. 2011a. Thermal isostasy of North China and its gravity isostasy and deep structure.ChineseJ.Geophys. (in Chinese), 54(11): 2864-2875, doi: 10.3969/j.issn.0001-5733.2011.11.016.

Chen S, Wang Q S, Zhu Y Q, et al. 2011b. Temporal and spatial features of isostasy anomaly using gravitational admittance model at eastern margin of Tibetan Plateau.ChineseJ.Geophys. (in Chinese), 54(1): 22-34, doi: 10.3969/j.issn.0001-5733.2011.01.004.

Chen S, Wang Q H, Wang Q S, et al. 2014. The 3D density structure and gravity change of LudianMS6.5 Yunnan epicenter and surrounding regions.ChineseJ.Geophys. (in Chinese), 57(9): 3080-3090, doi: 10.6038/cjg20140933.

Chen S, Wang Q S. 2015. Gravity anomalies and the distributions of inhomogeneous masses in the crust of Mongolia and its surrounding regions.ChineseJ.Geophys. (in Chinese), 58(1): 79-91, doi: 10.6038/cjg20150107.

Chen Y L, Niu F L, Liu R F, et al. 2010. Crustal structure beneath China from receiver function analysis.JournalofGeophysicalResearch:SolidEarth, 115(B3): B03307.

Fu G Y, Zhu Y Q, Gao S H, et al. 2013. Discrepancies between free air gravity anomalies from EGM2008 and the ones from dense gravity/GPS observations at west Sichuan Basin.ChineseJ.Geophys.(in Chinese), 56(11): 3761-3769, doi: 10.6038/cjg20131117.

Fullea J, Fernàndez M, Zeyen H. 2008. FA2BOUG-A FORTRAN 90 code to compute Bouguer gravity anomalies from gridded free-air anomalies: Application to the Atlantic-Mediterranean transition zone.Computers&Geosciences, 34(12): 1665-1681.

Gómez-Ortiz D, Agarwal B N P. 2005. 3DINVER. M: a MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by Parker-Oldenburg′s algorithm.Computers&Geosciences, 31(4): 513-520.

Guo L H, Meng X H, Shi L, et al. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent.ChineseJ.Geophys.(in Chinese), 55(12): 4078-4088, doi: 10.6038/j.issn.0001-5733.2012.12.020.

Hirt C. 2013. RTM gravity forward-modeling using topography/bathmetry data to improve high-degree global geopotential models in the coastal zone.MarineGeod., 36(2): 183-202.

Hirt C, Claessens S, Fecher T, et al. 2013. New ultra-high resolution picture of Earth′s gravity field.GeophysicalResearchLetters, 40(16): 4279-4283, doi: 10.1002/grl.50838.

Langston C A. 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves.BulletinoftheSeismologicalSocietyofAmerica, 67(3): 713-724.

Laske G, Masters. G, Ma Z, et al. 2013. Update on CRUST1.0—A 1-degree Global Model of Earth′s Crust. ∥EGU General Assembly 2013, Abstract EGU2013-2658, Vienna, Austria.

Li Y H, Gao M T, Wu Q J. 2014. Crustal thickness map of the Chinese mainland from teleseismic receiver functions.Tectonophysics, 611: 51-60.

Meng X H, Shi L, Guo L H, et al. 2012. Multi-scale analyses of transverse structures based on gravity anomalies in the northeastern margin of the Tibetan Plateau.ChineseJ.Geophys.(in Chinese), 55(12): 3933-3941, doi: 10.6038/j.issn.0001-5733.2012.12.006.

Oldenburg D W. 1974. The inversion and interpretation of gravtiy anomalies.Geophysics, 39(4): 526-536.

Parker R L. 1973. The rapid calculation of potential anomalies.Geophys.J.R.Astron.Soc., 31(4): 447-455.

Pavlis N K, Holmes S A, Kenyon S C, et al. 2008. An Earth Gravitational Model to Degree 2160: EGM2008. ∥Presented at the 2008 General Assembly of the European Geoscience Union, Vienna, Austria, 13-18.

Pavlis N K, Holmes S A, Kenyon S C, et al. 2012. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008).JournalofGeophysicalResearch, 117(B4): B04406.

Sandwell D T, Smith W H F. 1997. Marine gravity anomaly from Geosat and ERS 1 satellite altimetry.JournalofGeophysicalResearch:SolidEarth, 102(B5): 10039-10054.

Sandwell D T, Smith W H F. 2009. Global marine gravity from retracked Geosat and ERS-1 altimetry: Ridge segmentation versus spreading rate.JournalofGeophysicalResearch, 114(B4): B01411

Shin Y H, Choi K S, Xu H Z. 2006. Three-dimensional forward and inverse models for gravity fields based on the Fast Fourier Transform.Computers&Geosciences, 32(6): 727-738.

Smith W H F, Sandwell D T. 1997. Global sea floor topography from satellite altimetry and ship depth soundings.Science, 277(5334): 1956-1962.

Wang P, Wang L S, Mi N, et al. 2010. Crustal thickness and averageVp/Vsratio variations in southwest Yunnan, China, from teleseismic receiver functions.JournalofGeophysicalResearch:SolidEarth, 115(B11): B11308.

Xu L L, Rondenay S, van der Hilst R D. 2007. Structure of the crust beneath the southeastern Tibetan Plateau from teleseismic receiver functions.PhysicsoftheEarthandPlanetaryInteriors, 165(3-4): 176-193.

Xu Q, Zhao J M. 2008. A review of the receiver function method.ProgressinGeophysics(in Chinese), 23(6): 1709-1716.

Xu W M, Chen S, Shi L. 2014. Seismic Activity and Gravity Anomaly Characteristics of Yutian in Xinjiang and Surrounding Regions.EarthScience-JournalofChinaUniversityofGeosciences(in Chinese), 39(12): 1831-1841.

Zhang C Y, Guo C X, Chen J Y, et al. 2009. EGM 2008 and its application analysis in Chinese mainland.ActaGeodaeticaetCartographicaSinica(in Chinese), 38(4): 283-289.

Zheng Y, Ge C, Xie Z J, et al. 2013. Crustal and upper mantle structure and the deep seismogenic environment in the source regions of the Lushan earthquake and the Wenchuan earthquake.ScienceChina:EarthSciences, 56(7): 1158-1168.

附中文参考文献

陈石, 王谦身, 徐伟民等. 2011a. 华北地区热均衡、重力均衡与深部构造. 地球物理学报, 54(11): 2864-2875, doi: 10.3969/j.issn.0001-5733.2011.11.016.

陈石, 王谦身, 祝意青等. 2011b. 青藏高原东缘重力导纳模型均衡异常时空特征. 地球物理学报, 54(1): 22-34, doi: 10.3969/j.issn.0001-5733.2011.01.004.

陈石, 王青华, 王谦身等. 2014. 云南鲁甸MS6.5地震震源区和周边三维密度结构及重力场变化. 地球物理学报, 57(9): 3080-

3090, doi: 10.6038/cjg20140933.

陈石, 王谦身. 2015. 蒙古及周边地区重力异常和地壳不均匀体分布. 地球物理学报, 58(1): 79-91, doi: 10.6038/cjg20150107.

付广裕, 祝意青, 高尚华等. 2013. 川西地区实测自由空气重力异常与EGM2008模型结果的差异. 地球物理学报, 56(11): 3761-3769, doi: 10.6038/cjg20131117.

郭良辉, 孟小红, 石磊等. 2012. 优化滤波方法及其在中国大陆布格重力异常数据处理中的应用. 地球物理学报, 55(12): 4078-4088, doi: 10.6038/j.issn.0001-5733.2012.12.020

孟小红, 石磊, 郭良辉等. 2012. 青藏高原东北缘重力异常多尺度横向构造分析. 地球物理学报, 55(12): 3933-3941, doi: 10.6038/j.issn.0001-5733.2012.12.006.

徐强, 赵俊猛. 2008. 接收函数方法的研究综述. 地球物理学进展, 23(6): 1709-1716.

徐伟民, 陈石, 石磊. 2014. 新疆于田及周边地区地震活动性与重力异常特征. 地球科学-中国地质大学学报, 39(12): 1831-1841.

章传银, 郭春喜, 陈俊勇. 2009. EGM 2008地球重力场模型在中国大陆适用性分析. 测绘学报, 38(4): 283-289.

郑勇, 葛粲, 谢祖军等. 2013. 芦山与汶川地震震区地壳上地幔结构及深部孕震环境. 中国科学: 地球科学, 43(6): 1027-1037.

(本文编辑 汪海英)

Joint optimal inversion of gravity and seismic data to estimate crustal thickness of the southern section of the north-south seismic belt

CHEN Shi, ZHENG Qiu-Yue*, XU Wei-Min

InstituteofGeophysics,ChinaEarthquakeAdministration,Beijing100081,China

The south section of the north-south seismic belt is located in the southeast of the eastern Himalayan syntaxis, which has a complicated crustal structure, crustal deformation history and fault system, where earthquakes frequently occur. We used the joint optimal inversion approach, which includes gravity and seismology, with the known observation dataset and model to study the Moho interface deformation in the crust of this area. The aim of our study is to minimize the uncertainties and model discrepancies in different geophysical approaches. The methodology and results would help understand the mechanics on deformation in deep crust. The refined model of Moho interface may help us to reveal earthquake generation conditions, which have important geological and geophysical significance.

In this paper, we propose a kind of the joint optimal inversion method using gravity and seismic data. We take advantage of the receiver function which accurately estimates the crustal thickness beneath seismic stations and define a new kind of evaluation function. We also find an effective way to control the inversion process by adjusting the scale and shifting factors. This new approach can effectively reduce the discrepancies between the gravity inversion and the receiver function. The results show that this approach can reduce the uncertainty in the results derived by the interpolated method due to the spatial non-uniform distribution of seismic stations.

On the basis of above data and methods, our results show that the regional Bouguer gravity anomalies along the south section of the north-south seismic belt can be estimated by means of the EGM2008 gravity anomaly model. The constrained crustal thickness derived from the seismic receiver function at the 67 seismic stations are compared with the Crust 1.0 model. The standard errors in the 67 crustal thickness values are about 3.73 km between the Crust 1.0 model and the receiver functions. We present a joint model based on gravity and seismic data. Its standard error is about 1.9 km between the gravity inversion and the receiver function model using our approach. This optimal model also introduces the Moho depth of Crust1.0 crustal model as the initial values to estimate the lateral density distribution.

The joint optimal inversion in this study can help us overcome the known non-uniqueness problem in geophysical inversion models and refine the crustal subsurface model derived by the geophysical methods. It is very difficult to build seismic stations with the uniform distribution in the region of the north-south seismic belt due to the limit of topographic relief. Our approach can provide a solution to optimally estimate the Moho interface. And the results based on our test at the south section of the north-south seismic belt are better than the existing models. It can serve as an alternative solution for refining the crustal density model. Its efficiency is better than the known interpolation method. The refined Moho interface model at the study area has the minimum discrepancies and will help us to better study deep crustal deformation and the seismogenic mechanism in this region.

Joint inversion of gravity and seismic; Crustal thickness; Anomaly separation; Subsurface inversion; Optimal model

10.6038/cjg20151105.

10.6038/cjg20151105

P315

2014-12-25，2015-09-28收修定稿

地震行业科研专项(201508006)，中国地震科学台阵探测项目(201308011)，国家自然科学基金青年项目(41104046)，中国地震局地球物理研究所基本科研业务专项(DQJB13B19)，地震科技星火计划项目(XH14014Y)联合资助．

陈石，男，1980年生，固体地球物理学博士，副研究员，主要从事重力位场方法和地球动力学方面的研究．E-mail:chenshi@cea-igp.ac.cn*通讯作者 郑秋月，女，1992年生，硕士研究生，主要从事重力位场反演方法和应用研究.E-mail:qiuyuezheng@gmail.com

陈石, 郑秋月, 徐伟民. 2015. 南北地震带南段地壳厚度重震联合最优化反演.地球物理学报，58(11):3941-3951,

Chen S, Zheng Q Y, Xu W M. 2015. Joint optimal inversion of gravity and seismic data to estimate crustal thickness of the southern section of the north-south seismic belt.ChineseJ.Geophys. (in Chinese),58(11):3941-3951,doi:10.6038/cjg20151105.