CN106886047A - 一种接收函数和重力联合反演地壳厚度和波速比的方法 - Google Patents

Abstract

本发明属于地球动力学研究技术领域，具体涉及一种接收函数和重力反演地壳厚度和波速比的方法。该方法包括如下步骤：读入接收函数H‑κ叠加反演的H与κ及区域重力异常数据；选择待处理台站，设置窗口大小及网格间距；将窗口范围内的所有H、κ及重力数据网格化；求取密度参数及噪声参数；采用重力反演的似然估计法计算H‑κ似然谱；将接收函数H‑κ叠加图与重力H‑κ似然图归一化并点乘，得到新的H‑κ谱，从中拾取最优的H和κ作为联合反演结果；选择下一个待处理台站继续前述步骤。本发明通过重力似然图与接收函数叠加图相互加权实现最优H和κ的拾取，为解决复杂情况下接收函数反演的多解性问题提供有效途径。

Description

一种接收函数和重力联合反演地壳厚度和波速比的方法

技术领域

本发明属于地球动力学研究技术领域，具体涉及一种接收函数和重力反演地壳厚度和波速比的方法。

背景技术

地壳厚度和波速比等参数是研究地壳结构和物质组成的重要参数,对研究地球动力学过程具有重要意义。20世纪80年代以来，随着宽频带地震台网观测技术的不断完善和发展，远震体波接收函数方法的发展与应用使地壳上地幔速度结构的认识得到了空前的提升。

接收函数是由远震径向分量和垂直分量反褶积得到的时间序列，由于扣除了震源效应和传播路径的影响，因此更适合用来研究地壳和上地幔间断面的结构信息。接收函数的H-κ叠加法,用于估计地壳厚度与波速比参数，通过波速比与泊松比的关系式可换算地壳泊松比参数。该方法不用人工挑选震相走时,对间断面上的速度差敏感，简单快速，近年来已广泛用于研究全球多个地区的地壳厚度和波速比，是目前获取地壳厚度(H)和波速比(κ)参数的常用方法。

然而实际应用中，有些地区往往会由于台基不稳定、地震数据质量差、台站下方构造复杂等问题，造成接收函数的多次反射震相不清晰、拾取困难，从而导致无法或很难从H-κ叠加图中拾取出合理的H和κ参数值。如图1(a)中最大值所在是一个区域，无法确定出相对应的H和κ，图1(b)中置信区间有三个，同样无法确定出合理的H和κ值。其中，台站下方构造复杂包括较厚沉积层、壳内低速层和倾斜地层等地质情况。沉积层和基底之间存在明显的速度界面，这在接收函数中会出现介于初至P波和转换Ps波之间的一个波峰，从而对H-κ叠加结果产生较大影响，且沉积层越厚，这种影响就越大。

目前减少沉积层影响的方法有固定值域模拟搜索技术、添加沉积层速度结构等。壳内低速层引起的问题在接收函数的H-κ叠加中主要表现为出现双Ps波峰，这种情况下就会出现叠加结果的严重偏差。

重力方法是研究地质结构和区域地质构造的重要方法之一，水平分辨率高。布格重力异常通常包含了地壳内部物质成分密度不均匀引起的重力异常(即地壳重力异常)，也包括了莫霍面起伏引起的重力异常(即莫霍面重力异常)。地壳厚度参数可描述莫霍面起伏变化，波速比参数则可描述地壳内部物质成分并与地壳密度参数相关，也即它可描述地壳内部密度分布的不均匀性。因此,布格重力异常除了与地壳密度和厚度直接相关外，还与地壳波速比密切相关，通过布格重力异常反演理论上可以获取这些参数。但是仅仅依靠重力反演一种方法不能获得较为可信的地壳厚度与波速比。

总体来讲，在地壳结构复杂的情况下接收函数与重力反演任何一种方法单方面都不能获得较为可信的地壳厚度与波速比。而重力可以为接收函数H-κ叠加法提供约束信息，降低反演的不确定性或多解性。

发明内容

为了解决由于台基不稳定、地震数据质量差、台站下方构造复杂等问题造成接收函数的多次反射震相不清晰、拾取困难而无法拾取出合理地壳厚度H和波速比κ的技术问题。本发明提供一种接收函数和重力联合反演获得地壳厚度及波速比的方法；本发明通过利用布格重力异常反演莫霍面深度并转换成地壳厚度，在接收函数H-k叠加图中拾取地壳厚度相对应的波速比，以区域内所有台站的H和κ作为初值，计算密度参数与噪声参数，然后在H和κ的参数空间内扫面H和κ，并用该H和κ替换窗口中心台站的H和κ，计算相应的重力似然图，然后将该重力似然图与接收函数H-κ叠加图归一化并点乘即可得到联合反演的最优H和κ。

本发明是通过以下技术方案实现的：

一种准确获取地壳厚度和波速比的方法，所述方法用于降低接收函数H-κ叠加法的不确定性，提高获取地壳厚度H和波速比κ的精度和效率，所述方法基于接收函数和重力进行联合反演，所述方法首先利用布格重力异常反演获取地壳厚度H初值，进而获得波速比κ初值；然后，对选择的任意一个待反演地震台站为中心进行单一台站联合反演，得到单一台站的联合反演的H-κ图，从而获取该地震台站的最佳估计的地壳厚度H和波速比κ值；最后，对研究区域内所有待反演站台进行上述单一台站联合反演过程，以获得研究区域内所有待反演地震台站的地壳厚度H和波速比κ；所述联合反演的H-κ图中地壳厚度H和波速比κ拾取简单；所述方法能够为地壳增厚和地壳变形机制模式的研究提供高精度的地壳厚度和波速比数值。

进一步地，所述方法包括如下步骤：

利用布格重力异常反演研究区域内的莫霍面深度并转换成地壳厚度H，提取待反演台站位置所对应的H作为联合反演的地壳厚度H初值；

以所述地壳厚度H为引导，在接收函数H-κ叠加图中拾取与所有待反演地震台站的地壳厚度相对应的波速比κ，作为联合反演的波速比κ初值；

选择任意一个待反演的地震台站为中心，搜索一定窗口内所有地震台站的地壳厚度H、波速比κ和布格重力异常数据；

对上述选择的地震台站进行单一台站联合反演，使用单一台站联合反演中获得的最佳估计的H和κ值替换中心地震台站的原H和κ值；其中，所述单一台站联合反演基于重力反演的H-κ似然图与接收函数反演的H-κ叠加图；

然后选择下一个待联合反演的台站，重复所述单一台站联合反演的过程，直到研究区域内所有待反演地震台站的H和κ被替换为最佳估计的H和κ值，从而获得研究区域内所有待反演地震台站的地壳厚度H和波速比κ。

进一步地，所述单一台站联合反演的过程具体为：

(1)将以所述待反演地震台站为中心的一定窗口内的所有地震台站的地壳厚度H、波速比κ和布格重力异常数据，分别网格化成相同大小的网格数据；

(2)用步骤(1)获得的网格化后的地壳厚度H、波速比κ和布格重力异常，结合线性迭代算法计算出莫霍面上下密度差Δρ_Moho和地壳密度随波速比的变化率

(3)利用求得的所述莫霍面上下密度差Δρ_Moho和地壳密度随波速比的变化率计算出模型理论重力异常；

(4)将布格重力异常与所述模型理论重力异常作差，求得异常残差，采用似然估计法分别计算出异常残差的均值μ与方差σ²，再计算似然函数值；

(5)以接收函数H-κ叠加图相同的范围和步长选取H和κ值，重复步骤(3)和(4)，形成重力反演的H-κ似然图；

(6)将接收函数反演的H-κ叠加图与重力反演的H-κ似然图相乘，得到联合反演的H-κ图，从联合反演的H-κ图中拾取最佳估计的H和κ值。

进一步地，利用布格重力异常反演研究区域内的莫霍面深度并转换成地壳厚度H，具体为：

用频率域滤波法对研究区域内的布格重力异常作异常分离获得莫霍面重力异常，再应用密度界面反演方法获得研究区域内的莫霍面深度分布，并加上地形换算为地壳厚度分布，提取研究区内所有待反演台站的位置所对应的地壳厚度值作为联合反演的地壳厚度H初值。

进一步地，步骤(2)中，所述布格重力异常由莫霍面重力异常与地壳重力异常组成，所述莫霍面上下密度差Δρ_Moho与地壳密度随波速比的变化率的计算使用以下公式：

Δg_B＝Δg_Moho+Δg_Crust (1)

其中,Δg_B、Δg_Moho及Δg_Crust分别代表布格重力异常、莫霍面重力异常及地壳重力异常；Δρ_Moho为莫霍面上下密度差，为地壳密度随波速比的变化率，为研究区地壳厚度的平均值，为研究区波速比的平均值，F{·}和F^-1{·}分别表示傅里叶变换和反傅里叶变换，G为万有引力常数，f为波数，c为地壳厚度变化的校正因子；H为地壳厚度，κ为波速比。

进一步地，步骤(3)中计算出模型理论重力异常具体为：将步骤(2)中计算获得的密度参数Δρ_Moho与反代入步骤(2)的公式(2)、(3)中计算出Δg_Moho和Δg_Crust，计算所述模型理论重力异常Δg_Model的公式为：

Δg_Model＝Δg_Moho+Δg_Crust (4)

其中：Δg_Moho是莫霍面重力异常，Δg_Crust为地壳重力异常。

进一步地，步骤(4)中，假设异常残差服从高斯分布，因而采用下式计算μ与σ²：

其中：是待反演地震台站所在点上的重力异常残差，n代表计算窗口内的网格点数。

进一步地，步骤(5)中，在与接收函数H-κ叠加图相同的H和κ参数范围内，变动所述待反演地震台站的H和κ，当其似然值最大时，拟合出的模型重力异常与布格重力异常最接近，故此时的H和κ为最佳估计值，重力反演的H-κ似然图采用如下算法：

其中：L(μ,σ²)表示最大似然值；μ与σ²分别表示异常残差的均值和方差，是待反演地震台站所在点上的重力异常残差；n代表计算窗口内的网格点数。

进一步地，步骤(6)具体为：将重力反演的H-κ似然图与接收函数H-k叠加图均归一化，通过两矩阵点乘实现相互加权，得到联合反演的H-κ图，从联合反演的H-κ图中拾取最佳估计的H和κ值。

进一步地，不断更新并迭代研究区域内所有待反演地震台站的H和κ，以避免由于地壳厚度H和波速比κ初值误差太大而引起的联合反演结果不准确的情况。

本发明的有益技术效果：

(1)重力似然反演采用滑动窗口，每次反演只对窗口内的小面积进行，提高计算精度；

(2)考虑布格重力异常由莫霍面重力异常和地壳重力异常组成，简化计算算法；

(3)针对接收函数计算中多次波不明显而无法估计H、κ的问题，采用重力异常界面反演得到莫霍面深度，进一步换算出地壳厚度初值，再引导从接收函数H-κ叠加图中拾取出波速比参数，减少计算时间；

(4)布格重力异常反演地壳厚度和波速比过程中先用似然估计估计出噪声的均值和方差，然后在H和κ扫描过程中计算相应的最大似然值；

(5)区域内所有待反演地震台站的多次迭代计算有效提高了地壳厚度与波速比的反演精度；解决了在天然地震接收函数数据量少、数据质量差、台站下方存在较厚沉积层或者地壳内部存在低速层时，计算地壳厚度和波速比不准确的问题，能够为地壳增厚和地壳变形机制模式的研究提供相对可靠的地壳厚度和波速比数值，为检验地壳增厚和变形机制模式提供重要依据。

(6)本发明通过重力似然图与接收函数叠加图相互加权实现最优H和κ的拾取，为解决复杂情况下接收函数反演的多解性问题提供有效途径。

附图说明

图1：a图为接收函数数据少或者质量较差时的接收函数H-κ叠加图；b图为存在壳内低速层时的接收函数H-κ叠加图；

图2为重力反演的H-κ似然图；

图3分别为简单模型和复杂模型的接收函数与重力联合反演图；

图4为单一台站处理流程图；

图5为研究区域内所有待反演地震台站处理流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细描述。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。

实施例1

地壳厚度和波速比是研究地壳结构和物质组成的重要参数，对研究地球动力学机制有着重要的意义。在检验地壳增厚和地壳变形机制模式的研究中，需要以地壳厚度和波速比作为研究基础。以地壳厚度和波速比在青藏高原东北缘的研究中的应用为例，青藏高原东北缘是青藏高原向东北方向生长的最新的和正在变形的重要组成部分，广泛发育的第四系褶皱和逆冲断裂等表明该区正遭受着地壳缩短，并且伴随着垂直差异性隆升作用，地壳明显增厚。国内外学者通过研究利用地壳厚度和波速比对青藏高原东北缘地壳增厚和变形机制提出了很多模式，具有代表性的有：垂直连续变形说、上地壳增厚说和下地壳流说。其中，垂直连续变形说对应的耦合地壳增厚模型假设上地壳和下地壳的增厚速率一致，波速比和密度不会随地壳的增厚有明显变化。上地壳增厚说对应的上地壳增厚模型假设上地壳增厚速率明显大于下地壳，平均波速比和密度随地壳增厚而降低。下地壳流说对应的下地壳流模型假设地壳增厚主要是因为相邻板块侵入到下地壳而引起，地壳平均波速比和密度随地壳增厚而增加。

由以上青藏高原东北缘的例子可知，地壳厚度和波速比是地壳增厚和地壳变形机制模式的研究的基础。并且，地壳厚度和波速比在地壳增厚和地壳变形机制模式的研究中的应用并不以地区为限制，可以被广泛应用于不同地区的研究。

另外，地壳厚度和波速比还能够用于研究地壳组分，例如，有学者对华北克拉通中西部地区的地壳厚度与波速比进行研究，结果显示该地区山区平均波速比主要集中分布于1.70～1.77之间,暗示山区块体较东部盆地地区地壳组成更富长英质,而缺少铁镁质成分。

而目前使用最广泛的天然地震接收函数H-κ叠加法在实际应用中，往往会由于台基不稳定、地震数据质量差、台站下方构造复杂等问题，造成接收函数的多次反射震相不清晰、拾取困难，从而导致无法或很难从H-κ叠加图中拾取出合理的地壳厚度H和波速比κ参数值。错误或不准确的地壳厚度与波速比对地壳增厚和地壳变形机制模式带来了很大的阻碍，会造成不能够获得正确的地壳增厚和地壳变形机制模式。

综上所述，在地壳增厚和地壳变形机制模式的研究中，为了获得正确的地壳增厚和地壳变形机制模式，亟需一种方法能够准确测定地壳厚度和波速比。

本实施例提供一种一种接收函数和重力反演地壳厚度和波速比的方法，所述方法能够更加准确地获取地壳厚度和波速比，能够确保为地壳增厚和地壳变形机制模式的研究提供高精度的地壳厚度和波速比数值，以获得正确的地壳增厚和地壳变形机制模式。所述方法包括如下步骤：

利用布格重力异常反演研究区域内的莫霍面深度并转换成地壳厚度H，提取待反演台站位置所对应的H作为联合反演的地壳厚度H初值；具体为：用频率域滤波法对研究区域内的布格重力异常作异常分离获得莫霍面重力异常，再应用密度界面反演方法获得研究区域内的莫霍面深度分布，提取待反演台站位置的莫霍面深度并加上地形换算为地壳厚度分布，作为联合反演的地壳厚度H初值。

以所述地壳厚度H为引导，在接收函数H-κ叠加图中拾取与所有待反演地震台站的地壳厚度相对应的波速比κ，作为联合反演的波速比κ初值；

选择任意一个待反演的地震台站为中心，搜索一定窗口内所有地震台站的地壳厚度H、波速比κ和布格重力异常数据；

对上述选择的地震台站进行单一台站联合反演，使用单一台站联合反演中获得的最佳估计的H和κ值替换中心地震台站的原H和κ值；其中，所述单一台站联合反演基于重力反演的H-κ似然图与接收函数反演的H-κ叠加图；

然后选择下一个待联合反演的台站，重复所述单一台站联合反演的过程，直到研究区域内所有待反演地震台站的H和κ被替换为最佳估计的H和κ值，从而获得研究区域内所有待反演地震台站的地壳厚度H和波速比κ。

其中，所述单一台站联合反演的过程具体为：

(1)将以所述待反演地震台站为中心的一定窗口内的所有地震台站的地壳厚度H、波速比κ和布格重力异常数据，分别网格化成相同大小的网格数据；

窗口大小设置一般遵循如下规律：网格数大于11×11网格，网格间距大于1/4的台站平均间距且小于台站平均间距；

(2)用步骤(1)获得的网格化后的地壳厚度H、波速比κ和布格重力异常，结合线性迭代算法计算出莫霍面密度差Δρ_Moho和地壳密度随波速比变化率

所述地壳密度随波速比变化率的计算使用以下公式：

Δg_B＝Δg_Moho+Δg_Crust (1)

其中,Δg_B、Δg_Moho及Δg_Crust分别代表布格重力异常、莫霍面重力异常及地壳重力异常；Δρ_Moho为莫霍面上下密度差,为地壳密度随波速比的变化率,为研究区地壳厚度的平均值,为研究区波速比的平均值,F{·}和F^-1{·}分别表示傅里叶变换和反傅里叶变换,G为万有引力常数,f为波数,c为地壳厚度变化的校正因子；H为地壳厚度，κ为波速比。附图4中

(3)利用求得的莫霍面上下密度差Δρ_Moho和地壳密度随波速比的变化率计算出模型理论重力异常；模型理论重力异常的计算具体为：将步骤(2)中计算获得的密度参数反代入步骤(2)的公式(2)、(3)中计算出Δg_Moho和Δg_Crust，计算所述模型理论重力异常Δg_Model公式为：

Δg_Model＝Δg_Moho+Δg_Crust (4)

其中：Δg_Moho是莫霍面重力异常，Δg_Crust为地壳重力异常。

(4)将布格重力异常与所述模型理论重力异常作差，求得异常残差，即异常噪音，采用似然估计法分别计算出异常残差的均值μ与方差σ²，再计算重力似然函数值；

假设异常残差服从高斯分布，因而采用下式计算μ与σ²：

其中：是待反演地震台站所在点上的重力异常残差，n代表计算窗口内的网格点数。

(5)以接收函数H-κ叠加图相同的范围和步长选取H和κ值，重复步骤(3)和(4)，形成重力反演的H-κ似然图(即附图4中的MLE(maximum likelihood estimation))，重复步骤(3)时继续使用步骤(2)中第一次计算的密度参数Δρ_Moho和重复步骤(4)时继续使用步骤(4)中第一次计算的噪声参数μ与σ²；

在与接收函数H-κ叠加图中H和κ相同的参数范围内，变动所述待反演地震台站的H和κ，当似然值最大时，拟合出的模型重力异常与布格重力异常最接近，故此时的H和κ为最佳估计值，重力反演的H-κ似然图采用如下算法：

其中：L(μ,σ²)表示最大似然值；μ与σ²分别表示异常残差的均值和方差，是待反演地震台站所在点上的重力异常残差；n代表计算窗口内的网格点数。

(6)将重力反演的H-κ似然图与接收函数H-k叠加图均归一化，通过两矩阵点乘实现相互加权，得到联合反演的H-κ图，从联合反演的H-κ图中拾取最佳估计的H和κ值。

在本实施例所述方法中，不断更新并迭代研究区域内所有待反演地震台站的H和κ，以避免由于地壳厚度H和波速比κ初值误差而引起的联合反演结果不准确的情况。

对上述方法进行模型测试：

步骤一、建立一个平均深度为35km的地壳厚度模型，一个平均波速比为1.8的波速比模型，网格大小均为21×21，网格间距为15km。设置Δρ_Moho＝0.5g/cm³，利用如下公式正演重力异常：

Δg_B＝Δg_Moho+Δg_Crust

其中,Δg_B、Δg_Moho及Δg_Crust分别代表布格重力异常、莫霍面重力异常及地壳重力异常；Δρ_Moho为莫霍面上下密度差,为地壳密度随波速比的变化率,为研究区地壳厚度的平均值,为研究区波速比的平均值,F{·}和F^-1{·}分别表示傅里叶变换和反傅里叶变换,G为万有引力常数,f为波数,c为地壳厚度变化的校正因子；H为地壳厚度，κ为波速比。

步骤二、以区域中心的地壳厚度和波速比为参考建立两种地壳模型，分别为简单模型和复杂模型，其地壳厚度均为35km，然后正演接收函数；

步骤三、模型区域内随机产生65个点位作为台站位置，其中有一个位于区域中心，假设为需要联合反演的台站，提取这些点上的地壳厚度与波速比，以此作为初始值；

步骤四、将65个台站的地壳厚度与波速比网格化成与重力一样的网格大小，即网格大小为21×21，网格间距为15km；

步骤五、使用上述网格化后布格重力异常及H、κ,用线性迭代算法计算出Δρ_Moho与计算出Δρ_Moho＝0.52g/cm³，

步骤六、利用求得的密度参数Δρ_Moho与计算出模型理论重力异常；

步骤七、将布格重力异常与模型重力异常作差求得异常噪音，采用似然估计法分别计算出μ＝0与σ²＝4.85，将这两个参数带入以下步骤中计算重力似然值；

步骤八、以接收函数H-κ叠加相同的范围和步长选取H和κ值，重复步骤五和六，形成重力反演的H-κ似然图，如图2，重复步骤五时继续使用第一次计算的密度参数，重复步骤六时继续使用第一次计算的噪声参数；

步骤九、分别将两种地壳模型的接收函数H-κ叠加图与重力反演的H-κ似然图归一化并点乘，得到联合反演的H-κ图，如图3，可见无论是对应于接收函数数据数量和质量问题的简单模型，还是对应于存在沉积层和壳内低速层的复杂模型，与重力联合后均得到了一个相对稳定的值，从中拾取最佳估计的H和κ值分别为35km/1.85与36km/1.84，与模型设定值的误差在可接受范围内，因此，模型试验结果证明该方法有效。

Claims (10)

1.一种获取地壳厚度和波速比的方法，所述方法用于降低接收函数H-κ叠加法的不确定性，提高获取地壳厚度H和波速比κ的精度和效率，其特征在于，所述方法基于接收函数和重力进行联合反演，所述方法首先利用布格重力异常反演获取地壳厚度H初值，进而获得波速比κ初值；然后，对选择的任意一个待反演地震台站为中心进行单一台站联合反演，得到单一台站的联合反演的H-κ图，从而获取该地震台站的最佳估计的地壳厚度H和波速比κ值；最后，对研究区域内所有待反演站台进行上述单一台站联合反演过程，以获得研究区域内所有待反演地震台站的地壳厚度H和波速比κ；所述联合反演的H-κ图中地壳厚度H和波速比κ拾取简单；所述方法能够为地壳增厚和地壳变形机制模式的研究提供相对可靠的地壳厚度和波速比数值。

2.根据权利要求1所述方法，其特征在于，所述方法包括如下步骤：

利用布格重力异常反演研究区域内的莫霍面深度并转换成地壳厚度H，提取待反演台站的位置所对应的H作为联合反演的地壳厚度H初值；

以所述地壳厚度H为引导，在接收函数H-κ叠加图中拾取与所有待反演地震台站的地壳厚度相对应的波速比κ，作为联合反演的波速比κ初值；

选择任意一个待反演的地震台站为中心，搜索一定窗口内所有地震台站的地壳厚度H、波速比κ和布格重力异常数据；

对上述选择的地震台站进行单一台站联合反演，使用单一台站联合反演中获得的最佳估计的H和κ值替换中心地震台站的原H和κ值；其中，所述单一台站联合反演基于重力反演的H-κ似然图与接收函数反演的H-κ叠加图；

然后选择下一个待联合反演的台站，重复所述单一台站联合反演的过程，直到研究区域内所有待反演地震台站的H和κ被替换为最佳估计的H和κ值，从而获得研究区域内所有待反演地震台站的地壳厚度H和波速比κ。

3.根据权利要求2所述方法，其特征在于，所述单一台站联合反演的过程具体为：

(1)将以所述待反演地震台站为中心的一定窗口内的所有地震台站的地壳厚度H、波速比κ和布格重力异常数据，分别网格化成相同大小的网格数据；

(2)用步骤(1)获得的网格化后的地壳厚度H、波速比κ和布格重力异常，结合线性迭代算法计算出莫霍面上下密度差Δρ_Moho和地壳密度随波速比的变化率

(3)利用求得的所述莫霍面上下密度差Δρ_Moho和地壳密度随波速比的变化率计算出模型理论重力异常；

(4)将布格重力异常与所述模型理论重力异常作差，求得异常残差，采用似然估计法分别计算出异常残差的均值μ与方差σ²，再计算似然函数值；

(5)以接收函数H-κ叠加图相同的范围和步长选取H和κ值，重复步骤(3)和(4)，形成重力反演的H-κ似然图；

(6)将接收函数反演的H-κ叠加图与重力反演的H-κ似然图相乘，得到联合反演的H-κ图，从联合反演的H-κ图中拾取最佳估计的H和κ值。

4.根据权利要求2所述方法，其特征在于，利用布格重力异常反演研究区域内的莫霍面深度并转换成地壳厚度H，具体为：

用频率域滤波法对研究区域内的布格重力异常作异常分离获得莫霍面重力异常，再应用密度界面反演方法获得研究区域内的莫霍面深度分布，并加上地形换算为地壳厚度分布，提取研究区内所有待反演台站的位置所对应的地壳厚度值作为联合反演的地壳厚度H初值。

5.根据权利要求3所述方法，其特征在于，步骤(2)中，所述布格重力异常由莫霍面重力异常与地壳重力异常组成，所述地壳密度随波速比的变化率的计算使用以下公式：

Δg_B＝Δg_Moho+Δg_Crust (1)

${\mathrm{\Δg}}_{Moho}={\mathrm{\Δ\ρ}}_{Moho}\·F^{-1}\{2\mathrm{\π}G\·F\{H-\stackrel{\‾}{H}\}\·e^{-f\stackrel{\‾}{H}}\}---\left(2\right)$

${\mathrm{\Δg}}_{Crust}=\∂{\mathrm{\ρ}}_{Crust}/\∂k\·F^{-1}\{2\mathrm{\π}G\[\frac{1-e^{-f\stackrel{\‾}{H}}}{f}F\{k-\stackrel{\‾}{k}\}-c\·e^{-f\stackrel{\‾}{H}}\]\}---\left(3\right)$

其中,Δg_B、Δg_Moho及Δg_Crust分别代表布格重力异常、莫霍面重力异常及地壳重力异常；Δρ_Moho为莫霍面上下密度差，为地壳密度随波速比的变化率，为研究区地壳厚度的平均值，为研究区波速比的平均值，F{·}和F^-1{·}分别表示傅里叶变换和反傅里叶变换，G为万有引力常数，f为波数，c为地壳厚度变化的校正因子；H为地壳厚度，κ为波速比。

6.根据权利要求5所述方法，其特征在于，步骤(3)中计算出模型理论重力异常具体为：将步骤(2)中计算获得的密度参数Δρ_Moho与反代入步骤(2)的公式(2)、(3)中计算出Δg_Moho和Δg_Crust，计算所述模型理论重力异常Δg_Model的公式为：

Δg_Model＝Δg_Moho+Δg_Crust (4)

其中：Δg_Moho是莫霍面重力异常，Δg_Crust为地壳重力异常。

7.根据权利要求3所述方法，其特征在于，步骤(4)中，假设异常残差服从高斯分布，因而采用下式计算μ与σ²：

$\mathrm{\μ}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\Δg}}_{B,i}^{noise}---\left(5\right)$

${\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{({\mathrm{\Δg}}_{B,i}^{noise}-\mathrm{\μ})}^2---\left(6\right)$

其中：是待反演地震台站所在点上的重力异常残差，n代表计算窗口内的网格点数。

8.根据权利要求3所述方法，其特征在于，步骤(5)中，在与接收函数H-κ叠加图相同的H和κ参数范围内，变动所述待反演地震台站的H和κ，当其似然值最大时，拟合出的模型重力异常与布格重力异常最接近，故此时的H和κ为最佳估计值，重力反演的H-κ似然图采用如下算法：

$L(\mathrm{\μ},{\mathrm{\σ}}^2)=\underset{i=1}{\overset{n}{\Π}}\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\·e^{\frac{{({\mathrm{\Δg}}_{B,i}^{noise}-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}}={\left(\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}^2}\right)}^{\frac{n}{2}}\·e^{\frac{\underset{i=1}{\overset{n}{\Σ}}{({\mathrm{\Δg}}_{B,i}^{noise}-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}}---\left(7\right)$

其中：L(μ,σ²)表示最大似然值；μ与σ²分别表示异常残差的均值和方差，是待反演地震台站所在点上的重力异常残差；n代表计算窗口内的网格点数。

9.根据权利要求3所述方法，其特征在于，步骤(6)具体为：将重力反演的H-κ似然图与接收函数H-k叠加图均归一化，通过两矩阵点乘实现相互加权，得到联合反演的H-κ图，从联合反演的H-κ图中拾取最佳估计的H和κ值。

10.根据权利要求3所述方法，其特征在于，不断更新并迭代研究区域内所有待反演地震台站的H和κ，以避免由于地壳厚度H和波速比κ初值误差太大而引起的联合反演结果不准确的情况。