地球物理反演方法是了解地球内部结构的重要手段，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或地壳厚度变化产生的重力异常成份.常用的重力异常分离方法，一般都是通过异常的频率特性而设计的，但是其实没有严格的物理定义场源深度和异常频率之间的关系.而这种分离异常的不准确性，是制约提高重力反演地壳厚度精度的主要障碍之一.本文的方法原理将进一步讨论这个问题，并提出解决方案.

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

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

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

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

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

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

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

图 1

Fig. 1

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

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

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

通过重力异常反演密度界面，在反演过程中，我们需要调整的参数主要有3个，可以用式(5)表示.3个参数组合的反演结果在第i个约束点上的反演值E_i^C和约束值E_i，通过式(6)来评判最优反演结果.第一个参数是需要给定界面密度差ρ；第二个需要对反演重力异常进行去“直流分量”处理，也就是等价于对异常进行一个布格板模型校正g_DC；第三个参数是需要指定一个平均界面起伏h₀，参与反演.这三个参数都影响最后密度界面反演的结果.

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

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

图 2

Fig. 2

图 2 重力反演密度界面的多解性参数 (a)不同密度差参数即尺度因子对结果的影响；(b)不同平均深度参数即平移因子对结果的影响.Fig. 2 The uncertainty parameters on the density interface inversed by gravity anomaly (a) Shows the effect of scale factor with difference density contrast; (b) shows the effect of shift factor with the difference averaged depth.

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

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

图 3

Fig. 3

图 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.

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

图 4

Fig. 4

图 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

Fig. 5

图 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′，可以给出波长约10 km以上的重力信号(Pavlis et al.，2012;Hirt et al.，2013)，采用该模型为基础的各种扩展重力模型也在不断发展(S and well et al.，1997; Smith et al.，1997; S and well et al.，2009; Hirt，2013;).在计算布格重力异常的过程中，首先使用Lambert坐标投影算法，将经纬度转换为大地坐标；然后，选择了FA2BOUG标准程序(Fullea et al.，2008)用于地形校正和相关重力改正，得到布格重力异常.另外，由于Moho面起伏引起的重力异常主要以低频分量体现，我们还采取高斯低通滤波器用于区域重力异常分离，图 6是我们计算得到的区域布格重力异常.异常在整个研究区内都呈现负值特征，最低异常高达-600 mGal，青藏高原东南缘的异常值最低.通过与其他来源的重力异常进行对比，其低频分量部分，都具有较好的一致性，我们认为该数据对于反演本文研究的Moho面深度问题，完全可以满足异常精度要求.

图 6

Fig. 6

图 6 基于EGM2008重力异常模型计算得到的研究区布格异常 布格重力异常计算采用FA2BOUG程序，异常分离采用50 km高斯滤波，网格间距重采样为20 km.Fig. 6 The Bouguer gravity anomaly derived by EGM2008 model Bouguer 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.

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

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

图 7

Fig. 7

图 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.

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

图 8

Fig. 8

图 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.

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

图 9

Fig. 9

图 9 多种方法反演的Moho面深度结果对比 (a)重震联合最优化反演方法与接收函数结果对比；(b) 重震联合最优化与Crust1.0模型结果对比；图中虚线表示偏离±5 km的区间范围.Fig. 9 Moho depth results analysis for the multiple approaches (a) shows the results for the joint approach of gravity-seismic and receiver function inversion; (b) shows the joint approach of gravity-seismic and Crust1.0 model. The dotted lines are the ±5 km interval range for result difference

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

表 1(Table 1)

表 1 反演结果对比分析 Table 1 Analysis on inversion 模型 最大值 最小值 平均值 均方差 C1-RF 6.20 -9.84 -1.03 3.73 GI-C1 7.56 -8.60 0.37 3.19 GI-RF 2.06 -4.87 -0.65 1.90 注：GI为重力反演结果，RF为接收函数结果，C1为Crust1.0模型结果，单位为km.

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

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

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

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

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

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

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

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

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

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

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