CN113740915A - 一种球坐标系重力和接收函数联合反演地壳结构参数的方法 - Google Patents

Abstract

本发明属于地球物理学反演研究技术领域，公开了一种球坐标系重力和接收函数联合反演地壳结构参数的方法，反演过程中同时拟合重力和接收函数数据，通过联合反演算法实现重力和接收函数的互补作用，减少单一数据体反演的多解性。并且，在联合反演中考虑到了地球曲率的影响，引入了球坐标系下Tesseroid单元体的正演方法；本方法考虑到重力在横向的高分辨率和接收函数在台站附近深度方向的高分辨率，从而获取更加精确的地壳结构参数。

Description

一种球坐标系重力和接收函数联合反演地壳结构参数的方法

技术领域

本发明属于地球物理学反演研究技术领域，具体涉及一种球坐标系重力和接收函数联合反演地壳结构参数的方法。

背景技术

地壳结构参数是研究壳幔结构，壳内组分的重要方法，对于研究地球深部动力学、区域构造和壳内组分迁移等有重要的参考意义，其中地壳厚度(或莫霍面深度)和地壳速度比(或泊松比)是地球物理反演中最常获取的两个参数。重力学和地震学是研究地壳结构参数常用的两种方法。

重力学对地球深度的界面异常有很好的揭示意义，前人多次采用重力学的方法来研究深度的莫霍面起伏，一般而言重力数据在横向上有较好的分辨率，可以获得高精度的横向壳幔结构，同时，在大尺度的地壳结构研究中，考虑到地球的实际曲率，将地球近似为球坐标的相关理论算法日渐成熟。但是，由于重力数据随着深度快速衰减，其深度方向的分辨率较差。在地震学领域，接收函数常被用来研究来获取地壳厚度和地壳平均速度比，接收函数在地震台站的下方的分辨率较高，但是接收函数数据易受到浅层沉积层和复杂介质的干扰，且受限于地震台阵分布的稀疏离散不规则，常规的接收函数结果没有横向的物性约束。

因此，有必要设计一种新的联合反演算法，综合考虑到重力数据在横向上的高分辨率和接收函数在纵向的高分辨率，同时，在联合反演的过程中考虑到地球曲率的影响，将地球近似为球体。通过联合反演算法实现重力和接收函数的互补作用，减少单一数据体反演的多解性，获取更加精确的地壳结构参数。

发明内容

针对现有技术存在的上述不足，本发明的目的在于提供一种球坐标系重力和接收函数联合反演地壳结构参数的方法。

为实现以上目的，本发明采用如下技术方案：

一种球坐标系重力和接收函数联合反演地壳结构参数的方法，将规则网格下的球坐标下重力和接收函数进行联合反演，求取相关的地壳结构参数，反演过程主要依赖于接收函数中的转换波：Ps和其多次反射波：PpPs和PsPs+PpSs相对于直达P波的到时信息和重力异常数据，将非线性反演简化为线性反演问题，其方程简化为：

Gm＝d (1)

其中G为数据相对于模型的偏导数矩阵，m为当前模型的扰动矩阵，d为数据的残差矩阵；

在实际的反演中加入正则化的参数，联合反演的线性系统方程可以表示为：

其中，联合反演的矩阵扰动模型m包含莫霍面的起伏扰动ΔH和地壳中平均波速比的扰动Δκ；G矩阵的第一列包含了重力数据相对于莫霍面模型的偏导数，上角标SB代表球坐标下重力数据，下角标H代表莫霍面模型，μ为重力数据在联合反演中的权重；G矩阵的第二列到第四列分别代表了接收函数对于矩阵的扰动，上角标t_Ps，t_PpPs和t_PsPs+PpSs分别代表了转换波Ps，多次反射波PpPs，PsPs+PpSs相对于直达P波的到时差，下角标H为接收函数到时相对于莫霍面模型的偏导数，下角标κ为接收函数到时相对于速度比模型的偏导数，γ₁,γ₂,γ₃分别代表转换波Ps和多次反射波PpPs，PsPs+PpSs在联合反演中的权重；G矩阵的最后三列代表正则化参数，其中矩阵L表示采用L₁范数的平滑原则，即在矩阵中的每一个点和不同方向的相邻点采用一阶平滑模型；ω_H为求解莫霍面模型的平滑权重，I为单位矩阵，λ_H和λ_κ分别为莫霍面模型和速度比模型的阻尼权重参数；该线性联合反演系统可以通过阻尼最小二乘求解。

进一步地，在联合反演系统公式(4)中，第一列球坐标下重力异常相对于莫霍面的偏导数矩阵，能够通过公式(1)球坐标系下Tesseroid单元体的正演方程的有限差分算法近似求解：

其中

为观测点

的重力异常，

分别为球坐标系下重力观测点的经度，纬度和半径，G为万有引力常量，ρ为Tesseroid单元体的密度，Tesseroid单元体在球坐标下经度，纬度和半径的范围为(λ₁,λ₂),

(r₁,r₂)，假定

为单元体内的任意一个场源点，其中λ′∈(λ₁,λ₂)，

r′∈(r₁,r₂)，l为重力观测点P到场源点Q的欧几里得距离，其中

进一步地，在联合反演系统公式(4)中，对于接收函数转换波和多次波到时相对于莫霍面或者速度比的偏导数(第二列到第四列的偏导数)，能够依据公式(2)进行推导：

其中，t_Ps及

分别为转换波和多次波相对于直达P波的到时差，H为地壳厚度，其中H＝H_moho+topo，H_moho为台站下方的莫霍面起伏，topo为台站相对于海平面的地形起伏)，κ＝V_p/V_s，p是射线参数；由于上地壳对于V_p不敏感，因此，能够假定V_p为定值(Zhuand Kanamori，2000)，κ的大小随着V_s的变化而变化，根据公式(2)至(4)能够求解台站下方的地壳厚度和平均波速比信息。

进一步地，联合反演开始需要给定初始的莫霍面模型和速度比模型，反演的目标是通过多次迭代同时拟合接收函数到时信息和球坐标系下重力异常数据，联合反演总的目标方程为

其中，RMS_RF和RMS_Gra分别为接收函数到时和重力异常的拟合均方根残差，反演共计有M行和N列个观测点，i,j代表当前点所在的行和列(i∈M,j∈N)，T和SB分别代表接收函数到时和球坐标系重力异常，T中的上角标表示转换波或多次反射波(Ps，PpPs，PsPs+PpSs)；T和B的下角标中，cal和real分别代表模型正演计算获取的接收函数到时信息(或球坐标系下重力异常)和真实提取的接收函数到时信息(或球坐标系下重力异常)；φ₁,φ₂,φ₃分别为不同波形的目标权重系数，其中φ₁+φ₂+φ₃＝1。

进一步地，在实际的数据中，转换波Ps的波形质量一般高于多次反射波，因此反演方程中的权重一般选取为γ₁＞γ₂＞γ₃；在本申请的一种优选实施方式中，联合反演中默认为γ₁＝1，γ₂＝0.5，γ₃＝0.25，在实际反演中，可以通过L-curve曲线分析来确定重力和接收函数相关的权重系数。

进一步地，在实际的接收函数中，转换波的振相比后续的多次反射波更加清晰，所以转换波比多次波保留更多的目标权重，一般而言φ₁＞φ₂＞φ₃；在本申请的一种优选实施方式中，测试中选择的目标系数φ₁、φ₂、Φ₃分别为0.6，0.3，0.1。

进一步地，所述球坐标系重力和接收函数联合反演地壳结构参数的方法，具体流程如下：

1)根据重力数据，得到重力预处理提取莫霍面的重力异常；根据地震三分量数据，进行接收函数提取和远震波形重建构建虚拟地震台网；

2)根据莫霍面的重力异常和初始地壳模型参数(H，K)得到重力异常残差:d^SB；根据虚拟地震台网和初始地壳模型参数(H，K)得到到时残差:

3)求解扰动矩阵：Gm＝d；

4)更新模型空间：m^k+1＝m^k+m；

5)判断RMS_RF和RMS_Gra是否小于阈值，“是”即可输出地壳结构参数(H，K)；“否”即按照更新后的模型：m^k+1，正演球坐标系重力异常和正演接收函数转换振相到时，选取迭代次数K＝K+1，更新数据残差矩阵：d^k＝d^k-1-d^real，继续求解扰动矩阵：Gm＝d，并重复后续动作，直到输出地壳结构参数(H，K)。

与现有技术相比，本发明反演过程中同时拟合重力和接收函数数据，通过联合反演算法实现重力和接收函数的互补作用，减少单一数据体反演的多解性。并且，在联合反演中考虑到了地球曲率的影响，引入了球坐标系下Tesseroid单元体的正演方法；且考虑到重力在横向的高分辨率和接收函数在台站附近深度方向的高分辨率，从而获取更加精确的地壳结构参数。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为重力和接收函数反演地壳结构信息流程图；

图2.(a)球坐标系莫霍面的起伏模型；(b)球坐标系地壳平均的波速比模型；其中三角形为虚拟地震台站的分布位置；

图3.(a)和(b)为虚拟台站(经度103.5°E，纬度32.5°N)壳幔的速度结构和台站下方接收函数的正演模型；图(c)为图2.a中莫霍面起伏所引起的重力异常；

图4.(a)和(b)分别反演中接收函数和重力异常拟合曲线，其中三角形为单独采用接收函数数据结果，圆点为采用重力和接收函数联合反演结果；(c)和(d)为联合反演迭代稳定后第六次的莫霍面和地壳波速比的反演结果；(e)和(f)分别为联合反演结果和图2中莫霍面模型和速度比模型的残差；其中图(c)到(f)中的黑色三角形为虚拟台站。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。

实施例1

一种球坐标系重力和接收函数联合反演地壳结构参数的方法，将规则网格下的球坐标下重力和接收函数进行联合反演，求取相关的地壳结构参数，反演过程主要依赖于接收函数中的转换波：Ps和其多次反射波：PpPs和PsPs+PpSs相对于直达P波的到时信息和重力异常数据，将非线性反演简化为线性反演问题，其方程简化为：

Gm＝d (6)

其中G为数据相对于模型的偏导数矩阵，m为当前模型的扰动矩阵，d为数据的残差矩阵；

在实际的反演中加入正则化的参数，联合反演的线性系统方程可以表示为：

其中，联合反演的矩阵扰动模型m包含莫霍面的起伏扰动ΔH和地壳中平均波速比的扰动Δκ；G矩阵的第一列包含了重力数据相对于莫霍面模型的偏导数，上角标SB代表球坐标下重力数据，下角标H代表莫霍面模型，μ为重力数据在联合反演中的权重；G矩阵的第二列到第四列分别代表了接收函数对于矩阵的扰动，上角标t_Ps，t_PpPs和t_PsPs+PpSs分别代表了转换波Ps，多次反射波PpPs，PsPs+PpSs相对于直达P波的到时差，下角标H为接收函数到时相对于莫霍面模型的偏导数，下角标κ为接收函数到时相对于速度比模型的偏导数，γ₁,γ₂,γ₃分别代表转换波Ps和多次反射波PpPs，PsPs+PpSs在联合反演中的权重，在实际的数据中，转换波Ps的波形质量一般高于多次反射波，因此反演方程中的权重一般选取为γ₁＞γ₂＞γ₃；在本实施例的联合反演中默认为γ₁＝1，γ₂＝0.5，γ₃＝0.25，在实际反演中，可以通过L-curve曲线分析来确定重力和接收函数相关的权重系数；G矩阵的最后三列代表正则化参数，其中矩阵L表示采用L₁范数的平滑原则，即在矩阵中的每一个点和不同方向的相邻点采用一阶平滑模型；ω_H为求解莫霍面模型的平滑权重，I为单位矩阵，λ_H和λ_κ分别为莫霍面模型和速度比模型的阻尼权重参数；该线性联合反演系统可以通过阻尼最小二乘求解。

进一步地，在联合反演系统公式(4)中，第一列球坐标下重力异常相对于莫霍面的偏导数矩阵，能够通过公式(1)球坐标系下Tesseroid单元体的正演方程的有限差分算法近似求解：

其中

为观测点

的重力异常，

分别为球坐标系下重力观测点的经度，纬度和半径，G为万有引力常量，ρ为Tesseroid单元体的密度，Tesseroid单元体在球坐标下经度，纬度和半径的范围为(λ₁,λ₂),

(r₁,r₂)，假定

为单元体内的任意一个场源点，其中λ′∈(λ₁,λ₂)，

r′∈(r₁,r₂)，l为重力观测点P到场源点Q的欧几里得距离，其中

公式(1)为球坐标下的Tesseroid单元体正演公式，该公式不存在数值解，前人通过相关的研究(Lietal.2011；Uiedaetal.，2016；AoweiHaoetal，2019)提出了高斯勒让德的方法，已经可以获取了高精度的解析解，可以为本发明的重力和接收函数联合反演提供支撑。

进一步地，在联合反演系统公式(4)中，对于接收函数转换波和多次波到时相对于莫霍面或者速度比的偏导数(第二列到第四列的偏导数)，能够依据公式(2)进行推导(远震P波在经过台站下方的莫霍面时，由于波阻抗差的存在P波会发生波形转换，接收函数记录到的分别为直达P波，莫霍面的转换波Ps，以及多次的反射波(包含PpPs和PpSs+PSPs)，其中，后续震相与初至P波的到时差可以表示为公式(2)。)：

其中，t_Ps及

分别为转换波和多次波相对于直达P波的到时差，H为地壳厚度，其中H＝H_moho+topo，H_moho为台站下方的莫霍面起伏，topo为台站相对于海平面的地形起伏)，κ＝V_p/V_s，p是射线参数；由于上地壳对于V_p不敏感，因此，能够假定V_p为定值(ZhuandKanamori，2000)，κ的大小随着V_s的变化而变化，根据公式(2)至(4)能够求解台站下方的地壳厚度和平均波速比信息。

进一步地，联合反演开始需要给定初始的莫霍面模型和速度比模型，反演的目标是通过多次迭代同时拟合接收函数到时信息和球坐标系下重力异常数据，联合反演总的目标方程为

其中，RMS_RF和RMS_Gra分别为接收函数到时和重力异常的拟合均方根残差，反演共计有M行和N列个观测点，i,j代表当前点所在的行和列(i∈M,j∈N)，T和SB分别代表接收函数到时和球坐标系重力异常，T中的上角标表示转换波或多次反射波(Ps，PpPs，PsPs+PpSs)；T和B的下角标中，cal和real分别代表模型正演计算获取的接收函数到时信息(或球坐标系下重力异常)和真实提取的接收函数到时信息(或球坐标系下重力异常)；φ₁,φ₂,φ₃分别为不同波形的目标权重系数，其中φ₁+φ₂+φ₃＝1，在实际的接收函数中，转换波的振相比后续的多次反射波更加清晰，所以转换波比多次波保留更多的目标权重，一般而言φ₁＞φ₂＞φ₃；优选的，测试中选择的目标系数Φ₁、Φ₂、Φ₃分别为0.6，0.3，0.1。

如图1所示，所述球坐标系重力和接收函数联合反演地壳结构参数的方法，具体流程如下：

1)根据重力数据，得到重力预处理提取莫霍面的重力异常；根据地震三分量数据，进行接收函数提取和远震波形重建构建虚拟地震台网；

2)根据莫霍面的重力异常和初始地壳模型参数(H，K)得到重力异常残差:d^SB；根据虚拟地震台网和初始地壳模型参数(H，K)得到到时残差:

3)求解扰动矩阵：Gm＝d；

4)更新模型空间：m^k+1＝m^k+m；

5)判断RMS_RF和RMS_Gra是否小于阈值，“是”即可输出地壳结构参数(H，K)；“否”即按照更新后的模型：m^k+1，正演球坐标系重力异常和正演接收函数转换振相到时，选取迭代次数K＝K+1，更新数据残差矩阵：d^k＝d^k-1-d^real，继续求解扰动矩阵：Gm＝d，并重复后续动作，直到输出地壳结构参数(H，K)。

为了验证本发明联合反演算法的有效性，通过建立理论模型来进行测试。

(1)重力和接收函数正演

构建的莫霍面起伏模型和速度比模型分别见图2.a和图2.b，模型空间的大小为6°×4°，虚拟台站在东西向和南北向的网格间距都为0.5°，其中三角形为虚拟台站的分布，在实际的联合反演中，虚拟台站可以通过远震P波波形重建来完成。此外，在合成的数据测试中，假定观测面为海平面，忽略了地形的起伏，即莫霍面深度和地壳厚度相同。图2.a中莫霍面变化的范围为45km至35km，平均值为40km，其整体趋势为自西向东逐渐的依次递减，图2.b的平均速度比模型自西向东依次递增，其变化范围为1.55至1.83(平均值为1.68)。

为了进一步说明如何构建联合反演并且如何确定相关反演的参数，构建了存在两个低速层的地壳结构。以图2中的一个台站(Moho＝39.6Km，K＝1.65)为例，其速度分布和合成的接收函数数据，分别见图3.a，图3.b。从图3.b中横向为不同的射线参数，范围为0.044到0.08(间距为0.003)，纵向的时间轴上可以清晰地看到受到干扰的转换波Ps及其后续的多次波，每一个台站的转换波和多次波相对于直达P波的相对到时可以通过波形的叠加获取；对于莫霍面起伏引起的重力异常，可以通过球坐标系正演公式(1)来计算，模型测试中假定参考界面为40Km，上下界面的剩余密度为0.5g/cm³，图2中莫霍面起伏所引起的重力异常见图3.c，其重力异常的范围为-71mGal至69mGal，整体的变化趋势为自西向高递增。

(2)反演过程

在实际的台站中，由于数据干扰的存在，分别将提取到的理论转换波、多次波的到时信息和重力异常加入一定幅度的随机噪声干扰，对于接收函数三个不同振相加入5％的随机干扰，对于重力数据分别加入5％的高斯噪音干扰，将加入噪音干扰之后的接收函数到时信息和重力异常数据作为后续的反演理论值，此外在反演过程中选取的速度Vp为6.3Km/s，剩余密度为0.5g/cm³，重力的界面参考深度为40Km，L-curve曲线分析得到平滑参数和阻尼参数的最佳点分别为

λ_H＝300，λ_κ＝8000；重力的最优权重μ＝20。图4为联合反演的结果，图4.a和4.b中表明联合反演结果相对于单独的接收函数反演而言，不仅同时拟合接收函数和重力数据，而且可以加快迭代的效率，使反演结果在接收函数和重力数据中得到最佳的折中点；图4.c和图4.d为本发明联合反演结果，其形态和理论模型存在较大的一致性；图4.e和图4.f为联合反演结果和理论模型的残差，对于莫霍面而言，联合反演结果与理论模型的最大差异为1.5km，平均差异0.3km可以很好的恢复模型，对于地壳平均波速比而言，联合反演结果相对于理论结果的最大差异为0.05，平均为0.023。模型测试表明，本发明提出的联合反演算法可以很好的恢复理论模型，通过接收函数和球坐标系重力数据的结合，减少了反演结果的不确定性和多解性。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (9)

1.一种球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，将规则网格下的球坐标下重力和接收函数进行联合反演，求取相关的地壳结构参数，反演过程主要依赖于接收函数中的转换波：Ps和其多次反射波：PpPs和PsPs+PpSs相对于直达P波的到时信息和重力异常数据，将非线性反演简化为线性反演问题，其方程简化为：

Gm＝d (1)

其中G为数据相对于模型的偏导数矩阵，m为当前模型的扰动矩阵，d为数据的残差矩阵；

在实际的反演中加入正则化的参数，联合反演的线性系统方程可以表示为：

其中，联合反演的矩阵扰动模型m包含莫霍面的起伏扰动ΔH和地壳中平均波速比的扰动Δκ；G矩阵的第一列包含了重力数据相对于莫霍面模型的偏导数，上角标SB代表球坐标下重力数据，下角标H代表莫霍面模型，μ为重力数据在联合反演中的权重；G矩阵的第二列到第四列分别代表了接收函数对于矩阵的扰动，上角标t_Ps，t_PpPs和t_PsPs+PpSs分别代表了转换波Ps，多次反射波PpPs，PsPs+PpSs相对于直达P波的到时差，下角标H为接收函数到时相对于莫霍面模型的偏导数，下角标κ为接收函数到时相对于速度比模型的偏导数，γ₁,γ₂,γ₃分别代表转换波Ps和多次反射波PpPs，PsPs+PpSs在联合反演中的权重；G矩阵的最后三列代表正则化参数，其中矩阵L表示采用L₁范数的平滑原则，即在矩阵中的每一个点和不同方向的相邻点采用一阶平滑模型；ω_H为求解莫霍面模型的平滑权重，I为单位矩阵，λ_H和λ_κ分别为莫霍面模型和速度比模型的阻尼权重参数；该线性联合反演系统能够通过阻尼最小二乘求解。

2.根据权利要求1所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，在联合反演系统公式(4)中，第一列球坐标下重力异常相对于莫霍面的偏导数矩阵，能够通过公式(1)球坐标系下Tesseroid单元体的正演方程的有限差分算法近似求解：

其中

为观测点

的重力异常，λ

r分别为球坐标系下重力观测点的经度，纬度和半径，G为万有引力常量，ρ为Tesseroid单元体的密度，Tesseroid单元体在球坐标下经度，纬度和半径的范围为(λ₁,λ₂),

(r₁,r₂)，假定

为单元体内的任意一个场源点，其中λ′∈(λ₁,λ₂)，

r′∈(r₁,r₂)，l为重力观测点P到场源点Q的欧几里得距离，其中

3.根据权利要求1所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，在联合反演系统公式(4)中，对于接收函数转换波和多次波到时相对于莫霍面或者速度比的偏导数，能够依据公式(2)进行推导：

其中，t_Ps及

分别为转换波和多次波相对于直达P波的到时差，H为地壳厚度，其中H＝H_moho+topo，H_moho为台站下方的莫霍面起伏，topo为台站相对于海平面的地形起伏)，κ＝V_p/V_s，p是射线参数；由于上地壳对于V_p不敏感，因此，能够假定V_p为定值，κ的大小随着V_s的变化而变化，根据公式(2)至(4)能够求解台站下方的地壳厚度和平均波速比信息。

4.根据权利要求1所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，联合反演开始需要给定初始的莫霍面模型和速度比模型，反演的目标是通过多次迭代同时拟合接收函数到时信息和球坐标系下重力异常数据，联合反演总的目标方程为：

其中，RMS_RF和RMS_Gra分别为接收函数到时和重力异常的拟合均方根残差，反演共计有M行和N列个观测点，i,j代表当前点所在的行和列(i∈M,j∈N)，T和SB分别代表接收函数到时和球坐标系重力异常，T中的上角标表示转换波或多次反射波(Ps，PpPs，PsPs+PpSs)；T和B的下角标中，cal和real分别代表模型正演计算获取的接收函数到时信息(或球坐标系下重力异常)和真实提取的接收函数到时信息(或球坐标系下重力异常)；φ₁,φ₂,φ₃分别为不同波形的目标权重系数，其中φ₁+φ₂+φ₃＝1。

5.根据权利要求1所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，反演方程中的权重一般选取为γ₁＞γ₂＞γ₃。

6.根据权利要求5所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，选取的参数为γ₁＝1，γ₂＝0.5，γ₃＝0.25。

7.根据权利要求4所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，φ₁＞φ₂＞φ₃。

8.根据权利要求7所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，Φ₁、φ₂、φ₃选取的参数分别为0.6、0.3、0.1。

9.根据权利要求4所述球坐标系重力和接收函数联合反演地壳结构参数的方法，其特征在于，具体流程如下：

1)根据重力数据，得到重力预处理提取莫霍面的重力异常；根据地震三分量数据，进行接收函数提取和远震波形重建构建虚拟地震台网；

2)根据莫霍面的重力异常和初始地壳模型参数(H，K)得到重力异常残差:d^SB；根据虚拟地震台网和初始地壳模型参数(H，K)得到到时残差:

3)求解扰动矩阵：Gm＝d；

4)更新模型空间：m^k+1＝m^k+m；

5)判断RMS_RF和RMS_Gra是否小于阈值，“是”即可输出地壳结构参数(H，K)；“否”即按照更新后的模型：m^k+1，正演球坐标系重力异常和正演接收函数转换振相到时，选取迭代次数K＝K+1，更新数据残差矩阵：d^k＝d^k-1-d^real，继续求解扰动矩阵：Gm＝d，并重复后续动作，直到输出地壳结构参数(H，K)。

Images