CN111337993A - 一种基于变密度变深度约束的重力密度界面反演方法 - Google Patents

Abstract

本发明公开了一种基于变密度变深度约束的重力密度界面反演方法，包括：对获取的研究区内重力异常数据通过二维傅里叶变换转换，在频率域中进行正反演计算，得出修正后的重力异常场反演计算递推公式；基于全球地壳模型和实测地震资料信息获取研究区内密度界面的密度、深度信息情况；以得到的已知密度、深度信息数据作为约束条件，利用得出的修正后重力异常场反演计算递推公式，在频率域中进行密度界面反演计算。本发明通过同时考虑密度和深度变化的频率域重力密度界面反演方法，计算速度快，还能方便快捷的对大量数据进行处理。且本发明以已知的先验信息为约束条件，进行密度界面反演，能极大的提高重力方法反演的准确性，与实际情况相贴近。

Description

一种基于变密度变深度约束的重力密度界面反演方法

技术领域

本发明涉及地球物理重力勘探技术领域，特别是涉及一种基于变密度变深度约束的重力密度界面反演方法。

背景技术

随着地球物理数据处理及解释水平的不断提高，重力勘探方法已经成为地球物理技术方法中不可缺少的一环。重力勘探方法是根据地下物质之间的密度差异来研究地球重力场变化特征及其与矿产、构造等地质问题的一种方法。造成地球重力场变化(重力异常)的因素包括从地表附近至地球深处的物质密度的不均匀分布，重力异常是地下各物性界面起伏和岩性不均匀变化的综合体现，通过重力异常的反演，可以获取地下密度界面的起伏。

当地质界面的上、下层的物质的密度存在差异时，就构成了密度分界。密度分界根据上、下层密度差异的地球物理含义的不同，可分别与地层界面、地壳界面(Moho面)、岩石圈界面等地质地球物理界面相对应。因此研究密度界面的分布形态对于Moho面起伏、岩石圈厚度、区域构造等研究具有十分重要的意义。

在地球物理方法进行密度界面的方法研究过程中，地震和重力方法是目前最为常用的方法和手段。地震方法对于密度界面的精细研究可取得较好的应用效果，但地震方法成本高，且多为剖面形式布设，无法快速获得大面积范围内的密度界面起伏特征。另外，在存在屏蔽层的区域，地震方法无法获得部分层位的信号。较地震方法而言，重力勘探具有经济、勘探深度大以及可快速获得面积性信息等优点，特别是随着卫星重力精度的提高，在研究地质体的横向不均匀性等方面，具有独特的优势。因此，对重力异常数据进行密度界面反演是了解地球内部结构的一种重要方法。

目前，主流的重力密度界面反演方法可分为空间域和频率域两类。在空间域方法中，典型的方法有经验公式法、直接迭代法、压缩质面法和非线性反演方法。频率域方法则以Parker、Oldenburg所提出Parker-Olden burg频率域密度界面迭代反演方法为代表。密度界面反演中，密度和深度是非常重要的参数，其取值将直接影响反演界面结果的准确性和可靠性。

而现有的重力密度界面反演方法存在如下缺陷及不足之处：

现有频率域中的重力密度界面反演方法仅采用密度和深度分别为一个不变的常数进行反演，由于地质界面的深度和密度在不同区域均不一样，前人针对于地质界面的密度横向或纵向深度变化情况仅采用一常数进行密度界面反演，不足以满足对实际情况的计算和拟合，特别是计算范围面积超过4°×4°时，这种常数计算所得界面与真实情况相差更远，如果用由此得出的界面开展地质解释将会导致错误的认识。

因此，本申请就是在上述现有技术的基础上，创设一种基于变密度变深度约束的重力密度界面反演方法，使其能同时考虑密度和深度变化的重力密度界面反演，为重力勘探方法更为可靠的获取地球内部圈层结构奠定坚实基础。

发明内容

本发明要解决的技术问题是提供一种基于变密度变深度约束的重力密度界面反演方法，使其能同时考虑密度和深度变化的重力密度界面反演，为重力勘探方法更为可靠的获取地球内部圈层结构奠定坚实基础，从而克服现有的重力密度界面反演方法的不足。

为解决上述技术问题，本发明提供一种基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述方法包括如下步骤：

(1)获取研究区范围内的重力数据，并对重力数据信息进行有效提取，得出由密度界面所反映出的重力异常信息，再将所述重力异常信息的数据分布通过二维傅里叶变换进行转换，在频率域中进行正反演计算，得出修正后的重力异常场的反演计算递推公式：

式中，F[·]表示对括号内对应量的傅里叶变换，Δh为界面的修正量，r为扰动源坐标r在xy平面上的投影，上角标b、b-1、b+1为迭代次数序号，G为万有引力常数，ρ为密度差，ζ为界面S下某点的横坐标值；Δg_inv为实测的布格重力异常与第b次模型

正演获得的重力异常之差，即为模型的改正量的重力异常，a为波数场的模值，n为导数阶数；F_s ^(b)为第b次迭代计算出深度正演出的重力异常，

为第b次的迭代向下延拓算子，Ne为迭代次数；

(2)基于全球地壳模型Crust1.0和实测地震资料信息获取研究区范围内密度界面的密度、深度信息分布情况；

(3)以步骤(2)得到的研究区密度界面已知的密度、深度信息分布数据作为约束条件，利用步骤(1)得出的修正后重力异常场的反演计算递推公式，在频率域中进行密度界面反演计算。

进一步改进，所述步骤(1)中在频率域的正反演计算过程中，还包括：

(a)对重力数据信息在进行频率域离散之前，采用插值法与边界对称延伸的扩边方法对重力数据进行合理的扩充及赋值；

其中，采用插值法与边界对称延伸的扩边方法进行重力数据的扩充和赋值的公式如下：

式中，A为扩边后的重力数据，i为第i列，j为第j行，M为x方向的网格点数，l为待保留频谱成分的带宽，对外围点做线性插值，保持四周边界为该物理场的平均值；

(b)对扩充及赋值后的重力数据采用二维傅里叶变换方法，从时间域转换成频率域中进行计算；

(c)对转换成的频率域的重力数据采用余弦低通滤波方法进行滤波；

其中，余弦低通滤波方法采用如下公式进行：

式中，a为波数场的模值，k_s为高截止波数，k_w为低通波数，一般k_w＝0.75k_s。

进一步改进，所述步骤(1)中在频率域的正反演计算过程中，采用Parker-Oldenburg正反演计算公式，得出单点密度界面下的反演计算公式为：

式中，

为波数场的模值，

重力异常的波谱，

起伏深度的波谱，n为导数阶数，ρ为密度差，z₀为平均深度，

为向下延拓算子，G为万有引力常数；

和，二维连续密度界面h(x,y)在地表z＝z₀平面上，产生的重力异常场数据的正演计算公式为：

式中，z₀为地表观测平面高度；r₀为观测点的空间坐标矢量；

为扰动源坐标r在xy平面上的投影；G为万有引力常数；

为物质层的二维密度分布；F[·]表示对括号内对应量的傅里叶变换；a为波数场的模值；Δg为重力异常，

为向下延拓算子、

为密度界面的二维深度分布、n为导数阶数。

进一步改进，对二维连续密度界面的重力异常场数据的正演计算公式进行级数展开，提取第一项，得到简化的迭代反演表达式如下：

第一次迭代令h⁽⁰⁾＝0，通过上式中的第一项获得

然后再代入式子右端，计算获得

进一步改进，在迭代反演计算中，将反演界面的初值修改为：

式中，F[·]表示对括号内对应量的傅里叶变换；

为密度界面深度的初始模型，ξ为小于1的分配系数，一般为0.8ˉ0.95，是对级数项的补充，z₀为反演界面的平均深度，Δg为重力异常；

反演计算通过双重迭代：一是逐次平移反演层的下界面，通过正演公式计算新模型的重力贡献值；二是计算上界面的修正量，即得递推公式：

式中，上角标b，b-1、b+1为迭代次数序号，Δg_inv为实测的布格重力异常与第b次模型

正演获得的重力异常之差，即为模型的改正量的重力异常，n为导数阶数。

通过上面所述的重力异常数据迭代反演结果，与上一次迭代反演结果的对比误差，判断误差是否满足预定误差ε：

式中，M，N为x,y方向的网格点数，i,j分别为x,y方向的坐标，ε为前后两次迭代的误差；

为第b+1次的密度界面模型深度改正量，

为第b次的密度界面模型深度改正量；

若满足，则反演结束，输出反演结果数据；若不满足，则重复计算，直至满足预定误差为止。

进一步改进，所述步骤(3)中对搜集的研究区域密度界面的密度差约束，是通过将重力数据划分成与密度差信息分布矩阵相同大小的网格区域，对不同区域的重力数据采用与之对应区域的密度差值来进行。

进一步改进，所述步骤(3)中研究区域密度界面的深度约束，是通过二维约束矩阵的方式进行，若研究区域内的深度控制点有d个，则构建含有d个点的约束域U，其二维约束矩阵为：

式中，i和j分别为深度控制点矩阵的x,y方向坐标；

若反演的重力数据点与深度控制点重合，则直接将深度控制点的深度作为该点的重力反演解；反之，若不重合，则将重力反演解的深度作为反演的深度值。

进一步改进，所述步骤(3)中，还包括对深度控制点进行模型深度修改的步骤，其模型深度修改的公式如下：

式中，h为已知约束点深度，Δh为每次迭代反演之后，转换到空间域的深度值，Δh_new为约束之后的模型深度值，i和j分别为已知约束点矩阵的坐标，d为已知约束点的个数，Δh(x_i,y_i)为通过迭代反演获得界面深度值，h(x_i,y_i)为约束深度点的深度值。

采用这样的设计后，本发明至少具有以下优点：

本发明通过同时考虑密度和深度变化的重力密度界面反演方法，在频率域中进行计算反演，具有计算速度快的优势，可方便快捷的对大量数据进行处理。且本发明以已知的先验信息为约束条件，进行密度界面反演，能极大的提高重力方法反演的准确性，与实际情况相贴近。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是本发明基于变密度变深度约束的重力密度界面反演方法中约束条件下的界面反演方法的流程示意图。

图2是本发明基于变密度变深度约束的重力密度界面反演方法中单点p下密度界面S计算重力异常的示意图。

图3是本发明基于变密度变深度约束的重力密度界面反演方法中二维连续界面下三维理论体模型密度界面计算重力异常的示意图。

图4是本发明基于变密度变深度约束的重力密度界面反演方法中下界面、上界面修正移动的示意图。

图5是本发明基于变密度变深度约束的重力密度界面反演方法中界面反演区的密度差分块计算示意图，其中，大框内为研究区域，小网格中数值为该分块中的密度差。

具体实施方式

下面对本发明的具体实施方式做进一步的说明：

参照附图1所示，首先确定研究区大小及分辨率，进行数据准备工作：

利用基于全球地壳模型Crust1.0和实测地震资料信息获取研究区范围内密度界面的密度、深度信息分布情况。

如由Indensity_input.dat:搜集研究区密度界面的密度差二维分布，包含密度差和对应的位置坐标文件。

由Depth_input.dat:搜集研究区密度界面的深度二维分布，包含深度和对应的位置坐标文件。

由Bouguer_input.dat:输入研究区密度界面的重力数据，包含研究区的重力异常数据和对应的位置坐标文件。

其次，对得到的重力数据进行反演计算。计算步骤如下：

(1)对重力数据进行预处理，提取出由密度界面所反映出的重力异常信息。预处理步骤包括：

(a)对输入的重力数据在进行频率域离散之前，采用插值法与边界对称延伸的扩边方法进行合理的扩充及赋值。

其中，采用插值法与边界对称延伸的扩边方法进行重力数据的扩充和赋值的公式如下：

式(1)中，A为扩边后的重力数据，i为第i列，j为第j行，M为x方向的网格点数，l值为待保留频谱成分的带宽，对外围点做线性插值，保持四周边界为该物理场的平均值。对数据进行合理的扩充赋值，能确保正演公式两个先决条件的成立。

(b)对扩充及赋值后的重力数据采用二维傅里叶变换方法，从时间域转换成频率域中进行计算；

(c)对转换成的频率域的重力数据采用余弦低通滤波方法进行滤波；

其中，余弦低通滤波方法采用如下公式进行：

式(2)中，a为波数场的模值，k_s为高截止波数，k_w为低通波数，一般设置k_w＝0.75k_s。

(2)根据重力异常信息数据反演计算初始界面的起伏情况。

通过重力异常的波谱来反演计算界面S的起伏Δh的频率域的界面迭代反演公式如下：

上式(3)中，

为波数场的模值，

重力异常的波谱，

起伏深度的波谱，n为导数阶数，ρ为密度差，z₀为平均深度，

为向下延拓算子，G为万有引力常数。

该公式(3)的推导过程为：

参照附图2所示，单点密度界面下某点Q(ξ,η,ζ)处体积元对地表观测点p(x,y,0)所产生的引力场在直角坐标系(Z轴垂直向下)的基本计算公式为：

dΔg＝GKρdv (4)

其中，

G为万有引力常数，ρ为密度差，dv为体积分，Δg为引力异常。

设界面S的平均深度为z₀，将函数K在ζ＝z₀处采用泰勒级数展开，得到

上式(5)中：ξ,η,ζ分别为界面S下某点的横、纵、垂直坐标值；x,y分别为地面观测点横、纵坐标；n为导数阶数，K⁽ⁿ⁾(x-ξ,y-η,z₀)是在ζ＝z₀处K对ζ的n阶导数，z₀为平均深度。

利用上式将K对ζ从(h+Δh)到h积分，其中Δh是S面相对于z₀的距离，在ζ＝z₀以上的Δh为负，其下者为正，且Δh是(ξ,η)的函数。

上式(6)中，Δh为修正的深度值。

假设ρ在S界面以下仅为(ξ,η)的函数，则S界面以下全部物质对于P点产生的重力异常Δg为：

上式(7)中，ρ(ξ,η)为密度差值，Δh(ξ,η)为某点的修正深度值。由上式可知，重力异常Δg为函数K及其ζ＝z₀处的各阶垂向导数或一次垂向积分等函数与Δhⁿ(ξ,η)ρ(ξ,η)在ζ-η平面内相应的褶积式之和。根据褶积定理，令ρ(ξ,η)为常数获得重力异常Δg的波谱

式(8)中，

为波数场的模值，u,v分别为某点坐标ξ,η在x,y方向上的波数，

起伏深度的波谱，ρ为密度差，z₀为平均深度，

为向下延拓算子。

再通过重力异常的波谱来反演计算界面S的起伏Δh的频率域的界面迭代反演公式，则得出：

(3)输入上述搜集到的研究区密度界面的密度差和深度的二维分布数据，根据已知的密度界面的密度差二维分布，将研究区域划分为与密度差分布相适应的区域，在对不同区域的重力数据计算时，采用与之对应区域的密度差值进行计算，如附图5上述。

对于深度约束，可以通过如下方式引入到反演过程，假设研究区域内的深度控制点有d个，则构建含有d个点的约束域U，那么其二维约束矩阵为

式(9)中，M，N为深度控制点x,y方向的网格点数，i,j分别为深度控制点矩阵的x,y方向坐标。

若反演的重力数据点与深度控制点重合，则直接将深度控制点的深度作为该点的重力反演解；反之，若不重合，则将重力反演解的深度作为反演的深度值。

(4)将模型界面深度与已知的深度控制点进行对比，根据深度控制点来修改模型的深度。

在含有深度控制点值区域，重力反演解与深度控制点结果完全重合；在不含深度控制点的区域，所获得的重力反演解的值作为该点的深度值。通过机械方式改变深度控制点区域的深度，仅是对研究区域的密度界面进行点对点的约束，会造成反演结果界面数据的局部突变，且对于恢复重力异常具有较大误差。

为了避免对界面深度进行简单的替换，通过计算已知点位深度与每次迭代反演界面深度的误差算术平均值，然后对界面进行整体移动，从而既能对反演结果进行约束，又能防止界面的局部突变。如利用下式即可得出约束后的模型深度值Δh_new：

式(10)中h为已知约束点深度，Δh为每次迭代反演之后，转换到空间域的深度值，Δh_new为约束之后的模型深度值，i和j分别为已知约束点矩阵的坐标，d为已知约束点的个数，Δh(x_i,y_i)为通过迭代反演获得界面深度值，h(x_i,y_i)为约束深度点的深度值。

(5)对步骤(4)修改后的模型深度通过下式进行正演计算，计算出修改后模型的重力异常值。

二维连续密度界面h(x,y)在地表z＝z₀平面上，如附图3所示，产生的重力异常场的频率域表达式如下：

上式(11)中，z₀为密度界面的平面高度；r₀为观测点的空间坐标矢量；

为扰动源坐标r在xy平面上的投影；G为万有引力常数；

为物质层的二维密度分布；F[·]表示对括号内对应量的傅里叶变换；a为波数场的模值；Δg为重力异常。

(6)对修改后模型的重力异常值通过下式进行频率域的反演计算，

上式(12)中，F[·]表示对括号内对应量的傅里叶变换，Δh为界面的修正量，

为扰动源坐标r在xy平面上的投影，上角标b，b-1、b+1为迭代次数序号，G为万有引力常数，ρ为密度差，ζ为界面S下某点的横坐标值；Δg_inv为实测的布格重力异常与第b次模型

正演获得的重力异常之差，即为模型的改正量的重力异常，a为波数场的模值，n为导数阶数；F_s ^(b)为第b次迭代计算出深度正演出的重力异常，

为第b次的迭代向下延拓算子，Ne为迭代次数；

该式(12)的推导过程如下：

根据重力异常场的频率域表达式：

式(11)中，z₀为地表观测平面高度；r₀为观测点的空间坐标矢量；

为扰动源坐标r在xy平面上的投影；G为万有引力常数；

为物质层的二维密度分布；F[·]表示对括号内对应量的傅里叶变换；a为波数场的模值；Δg为重力异常，

为向下延拓算子、

为密度界面的二维深度分布、n为导数阶数。

对于(M,N)的网格点阵

上式(13)的s与t分别为x,y方向上的网格点距，M为行数，N为列数，m₁为x方向的坐标，n₁为y方向的坐标。

对二维连续密度界面的重力异常场进行级数展开后，提取出第一项，可获得迭代反演的表达式：

上式(14)中，F[·]表示对括号内对应量的傅里叶变换；h为起伏界面深度。

第一次迭代令h⁽⁰⁾＝0，通过上述式子中的第一项获得

然后再代入式子右端，计算获得

上述该二维连续密度界面的正反演公式还需要存在两个先决条件：(a)有限域外的界面趋于零；(b)界面起伏是有界可积的。

为改善反演计算的稳定性并加速运算，将反演界面的初值修改如下：

式(15)中，F[·]表示对括号内对应量的傅里叶变换；

为密度界面深度的初始模型，ξ为小于1的分配系数(一般为0.8～0.95)，是对级数项的补充。z₀为反演界面的平均深度。

反演计算过程中包含双重迭代：第一是逐次平移反演层的下界面，通过正演公式计算新模型的重力贡献值；二是计算上界面的修正量，如附图4所示，其递推公式得出下式：

(7)反演结束条件：对比前后两次迭代反演的结果，并计算其误差，通过误差是否满足预定误差来决定反演是否结束，误差计算公式：

式(16)中，M，N为x,y方向的网格点数，i,j分别为x,y方向的坐标，ε为前后两次迭代的误差；

为第b+1次的密度界面模型深度改正量，

为第b次的密度界面模型深度改正量。

若满足，则反演结束，输出反演结果数据；若不满足，则重复计算(4)～(7)步骤，直至满足预定误差为止。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims (9)

1.一种基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述方法包括如下步骤：

(1)获取研究区范围内的重力数据，并对重力数据信息进行有效提取，得出由密度界面所反映出的重力异常信息，再将所述重力异常信息的数据分布通过二维傅里叶变换进行转换，在频率域中进行正反演计算，得出修正后的重力异常场的反演计算递推公式：

式中，F[·]表示对括号内对应量的傅里叶变换，Δh为界面的修正量，

为扰动源坐标r在xy平面上的投影，上角标b、b-1、b+1为迭代次数序号，G为万有引力常数，ρ为密度差，ζ为界面S下某点的横坐标值；Δg_inv为实测的布格重力异常与第b次模型

正演获得的重力异常之差，即为模型的改正量的重力异常，a为波数场的模值，n为导数阶数；

为第b次迭代计算出深度正演出的重力异常，

为第b次的迭代向下延拓算子，Ne为迭代次数；

(2)基于全球地壳模型Crust1.0和实测地震资料信息获取研究区范围内密度界面的密度、深度信息分布情况；

(3)以步骤(2)得到的研究区密度界面已知的密度、深度信息分布数据作为约束条件，利用步骤(1)得出的修正后重力异常场的反演计算递推公式，在频率域中进行密度界面反演计算。

2.根据权利要求1所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述步骤(1)中在频率域的正反演计算过程中，还包括：

(a)对重力数据信息在进行频率域离散之前，采用插值法与边界对称延伸的扩边方法对重力数据进行合理的扩充及赋值；

其中，采用插值法与边界对称延伸的扩边方法进行重力数据的扩充和赋值的公式如下：

式中，A为扩边后的重力数据，i为第i列，j为第j行，M为x方向的网格点数，l为待保留频谱成分的带宽，对外围点做线性插值，保持四周边界为该物理场的平均值；

(b)对扩充及赋值后的重力数据采用二维傅里叶变换方法，从时间域转换成频率域中进行计算；

(c)对转换成的频率域的重力数据采用余弦低通滤波方法进行滤波；

其中，余弦低通滤波方法采用如下公式进行：

式中，a为波数场的模值，k_s为高截止波数，k_w为低通波数，一般k_w＝0.75k_s。

3.根据权利要求2所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述步骤(1)中在频率域的正反演计算过程中，采用Parker-Oldenburg正反演计算公式，得出单点密度界面下的反演计算公式为：

式中，

为波数场的模值，

重力异常的波谱，

起伏深度的波谱，n为导数阶数，ρ为密度差，z₀为平均深度，

为向下延拓算子，G为万有引力常数。

和，二维连续密度界面h(x,y)在地表z＝z₀平面上，产生的重力异常场数据的正演计算公式为：

式中，z₀为地表观测平面高度；r₀为观测点的空间坐标矢量；

为扰动源坐标r在xy平面上的投影；G为万有引力常数；

为物质层的二维密度分布；F[·]表示对括号内对应量的傅里叶变换；a为波数场的模值；Δg为重力异常，

为向下延拓算子、

为密度界面的二维深度分布、n为导数阶数。

4.根据权利要求3所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，对二维连续密度界面的重力异常场数据的正演计算公式进行级数展开，提取第一项，得到简化的迭代反演表达式如下：

第一次迭代令h⁽⁰⁾＝0，通过上式中的第一项获得

然后再代入式子右端，计算获得

5.根据权利要求4所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，在迭代反演计算中，将反演界面的初值修改为：

式中，F[·]表示对括号内对应量的傅里叶变换；

为密度界面深度的初始模型，ξ为小于1的分配系数，一般为0.8～0.95，是对级数项的补充，z₀为反演界面的平均深度，Δg为重力异常；

反演计算通过双重迭代：一是逐次平移反演层的下界面，通过正演公式计算新模型的重力贡献值；二是计算上界面的修正量，即得递推公式：

式中，上角标b，b-1、b+1为迭代次数序号，Δg_inv为实测的布格重力异常与第b次模型

正演获得的重力异常之差，即为模型的改正量的重力异常，n为导数阶数。

6.根据权利要求5所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，重力异常数据迭代反演的误差条件为计算前后两次迭代反演的误差值，通过判断误差是否达到预定误差ε要求，来决定反演是否结束：

式中，M，N为x,y方向的网格点数，i,j分别为x,y方向的坐标，ε为前后两次迭代的误差；

为第b+1次的密度界面模型深度改正量，

为第b次的密度界面模型深度改正量；

若满足，则反演结束，输出反演结果数据；若不满足，则重复计算，直至满足预定误差为止。

7.根据权利要求6所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述步骤(3)中对搜集的研究区域密度界面的密度差约束，是通过将重力数据划分成与密度差信息分布矩阵相同大小的网格区域，对不同区域的重力数据采用与之对应区域的密度差值来进行。

8.根据权利要求书7所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述步骤(3)中研究区域密度界面的深度约束，是通过二维约束矩阵的方式进行，若研究区域内的深度控制点有d个，则构建含有d个点的约束域U，其二维约束矩阵为：

式中，i和j分别为深度控制点矩阵的x,y方向坐标；

若反演的重力数据点与深度控制点重合，则直接将深度控制点的深度作为该点的重力反演解；反之，若不重合，则将重力反演解的深度作为反演的深度值。

9.根据权利要求书8所述的基于变密度变深度约束的重力密度界面反演方法，其特征在于，所述步骤(3)中，还包括对深度控制点进行模型深度修改的步骤，其模型深度修改的公式如下：

式中，h为已知约束点深度，Δh为每次迭代反演之后，转换到空间域的深度值，Δh_new为约束之后的模型深度值，i和j分别为已知约束点矩阵的坐标，d为已知约束点的个数，Δh(x_i,y_i)为通过迭代反演获得界面深度值，h(x_i,y_i)为约束深度点的深度值。

Images