CN116661014A - 一种用于任意变密度界面的反演方法 - Google Patents

Abstract

本发明公开了一种用于任意变密度界面的反演方法，具体涉及重力勘探技术领域。本发明基于三维棱柱体剖分的方式，在三维空间网格中构建三维密度界面模型，通过对密度界面正演函数进行泰勒级数展开并保留一次项，近似得到密度界面模型正演的线性函数并带入至正则化目标函数中，通过极小化目标函数得到密度界面反演迭代方程，利用线性共轭梯度算法求解密度界面反演模型，通过多次对密度界面进行迭代修正，反演得到研究区的密度界面和密度模型。本发明实现了在复杂地质构造情况下密度界面分布的准确获取，基于三维重力反演的快速算法，实现了重力密度界面的批量快速反演，为复杂地质构造的勘探开发奠定了基础。

Description

一种用于任意变密度界面的反演方法

技术领域

本发明属于重力勘探技术领域，具体涉及一种用于任意变密度界面的反演方法。

背景技术

重力密度界面与基底起伏、圈层空间展布、构造隆起、莫霍面等地质构造相对应，在研究区域构造、圈定油气远景区以及壳幔结构等方面具有重要的理论和实际意义。

频率域反演方法出现于20世纪70年代，该方法因计算速度快的优势得到了迅速发展与应用，空间域密度界面反演方法出现较早，最早可追溯到20世纪50年代，随后得到了广泛研究和使用，并形成了多种反演方法，典型的方法包括经验公式法、直接迭代法、脊回归法、正则化方法、压缩质面法、级数法和样条函数法等。其中直接迭代法、脊回归法和正则化方法是研究和使用最多的方法，其他方法的相关研究较少。

无论是在频率域正演还是在空间域正演，目前密度界面模型所采用的方法基本上都是利用棱柱体剖分模型构建密度界面模型，其通过将空间域在水平方向进行剖分，然后用棱柱体顶面或底面的埋深在垂直方向上表示表示密度界面的深度，该剖分方式在进行变密度正演时，需要将密度变化归纳为函数，然后再根据函数推导正演公式。这就导致必须将已知的密度变化用有限的函数进行表达，例如，采用只随深度变换的线性密度函数、二次多项式密度函数、三次多项式密度函数、抛物线密度函数、双曲线密度函数、指数密度函数等函数进行表达，或采用基于多项式的二维密度变化和三位密度变化函数等函数进行表达。

上述函数整体光滑，能够较好的拟合简单地质构造的密度分布，但是，受构造作用的影响、地层形态变化较大、地层缺失等密度突变情况的影响，实际地质构造往往很复杂，上述函数难以准确表达地质构造的密度分布。因此，对于大多数地质构造而言，现有的密度变化规律并不能准确表达地质构造实际的密度分布特征，亟需提出一种用于任意变密度界面的反演方法，实现对复杂地质构造密度分布的准确获取。

发明内容

本发明为了解决现有密度变化函数难以准确描述复杂地质构造密度分布特征的问题，提出了一种用于任意变密度界面的反演方法，本发明实现了在复杂地质构造情况下密度界面分布的准确获取，基于三维重力反演进行快速计算，有利于重力密度界面的批量快速反演。

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

一种用于任意变密度界面的反演方法，具体包括如下步骤：

步骤1，获取研究区的密度属性信息，在三维空间网格中构建三维密度界面模型，利用三维棱柱体正演确定研究区三维密度界面模型的重力场；

步骤2，基于研究区三维密度界面模型的正演重力场，确定密度界面正则化反演的目标函数，建立密度界面反演计算模型；

步骤3，基于密度界面反演计算模型进行密度界面反演，得到研究区的密度界面和密度模型。

优选地，所述步骤1中，利用三维棱柱体模型表示研究区的密度界面，根据研究区的密度属性信息确定密度界面上方区域和下方区域的密度属性信息。

优选地，所述步骤1中，基于三维棱柱体模型重力正演确定重力场为：

g＝g(h) (1)

式中，g为地表重力值，h为密度界面深度。

优选地，所述步骤2中，根据研究区的重力场，确定密度界面正则化反演的目标函数，如公式(2)所示：

式中，为密度界面正则化反演的目标函数，W_d为重力数据加权矩阵，g^obs为重力观测值，g(h)为重力场函数，λ为正则化参数，W_h为密度界面加权矩阵，h_pre为密度界面参考数据；

由于重力场函数g(h)为非线性函数，对重力场函数g(h)进行泰勒级数展开并保留一次项，得到：

其中，采用离散三维网格重力反演时雅可比矩阵A_i中元素A_jk的计算公式为：

式中，i为重力反演的计算次数，h_i为第i次重力反演计算所对应的密度界面深度，A_i为密度界面深度h_i所对应的雅可比矩阵，g_j为第j个观测点，h_k为第k个界面深度点，G_jk为界面深度点h_k所在三维棱柱体模型对观测点g_j的正演系数，m_k为三维棱柱体模型的密度值，Δh_k为三维棱柱体模型的纵向长度；

将泰勒级数展开后的重力场函数g(h)带入到密度界面正则化反演的目标函数中，并对密度界面深度h进行求导，得到：

式中，T为转置矩阵；

将密度界面正则化反演的目标函数的取值设置为零，得到密度界面反演迭代公式为：

将公式(6)两遍同时加上得到：

通过对公式(7)的左侧矩阵取逆，得到密度界面反演计算模型，如公式(8)所示：

优选地，所述步骤3中，利用线性共轭梯度算法求解密度界面反演计算模型，通过对密度界面进行多次迭代修正，反演得到研究区的密度界面和密度模型，具体包括以下步骤：

步骤3.1，反演参数初始化；

输入重力观测数据g^obs，选取密度界面参考数据h_pre，设置密度界面深度的初始值h₀、密度界面加权矩阵W_h、重力数据加权矩阵W_d和最大修正次数Nmax，预设修正精度和密度界面反演计算模型中已知界面之间的密度分布；

步骤3.2，采用离散三维网格进行重力反演，根据密度界面深度h_i以及密度界面深度h_i所对应的密度模型，计算密度界面深度h_i所对应的雅克比矩阵A_i；

步骤3.3，基于线性共轭梯度法求解密度界面反演计算模型，得到修正后密度界面深度h_i+1，确定密度界面深度h_i+1所对应的密度界面；

步骤3.4，计算密度界面深度的修正量Δh，当密度界面深度的修正量Δh小于预设的修正精度或修正次数已达到预设的最大修正次数时，则进入步骤3.5中，否则，则更新修正次数i，继续返回步骤3.2中进行修正；

步骤3.5，输出密度界面和密度模型。

优选地，所述步骤3.3中，基于线性共轭梯度法求解密度界面反演计算模型，获取修正后密度界面深度所对应密度界面，具体包括以下步骤：

步骤3.3.1，迭代计算参数初始化；

预设密度界面反演计算模型的梯度精度值ε，设置线性共轭梯度计算参数Q和b，设置密度界面反演计算模型的初始梯度r₀、迭代计算的初始搜索方向p₀和迭代更新初值m₀，r₀＝b-Qm₀，p₀＝r₀，m₀＝0.001；

步骤3.3.2，根据密度界面反演计算模型的梯度和迭代搜索方向，计算迭代搜索方向的步长，如公式(9)所示：

式中，α_k为第k次迭代计算时迭代搜索方向的步长，r_k为第k次迭代计算时密度界面反演计算模型的梯度值，p_k为第k次迭代计算的搜索方向；

步骤3.3.3，更新迭代计算结果和密度界面反演计算模型的梯度为：

m_k+1＝m_k+α_kp_k (10)

r_k+1＝r_k-α_kQp_k (11)

式中，m_k+1为更新后的迭代计算结果，m_k为第k次迭代计算的结果，r_k+1为更新后密度界面反演计算模型的梯度值；

步骤3.3.3，若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂小于预设梯度精度值ε，则根据第k次迭代计算的结果m_k，确定修正后密度界面深度h_i+1，如公式(12)所示：

h_i+1＝m_k+h_pre (12)

计算得到修正后密度界面深度h_i+1后，进入步骤3.3.5中；

若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂不小于预设梯度精度值ε，则进入步骤3.3.4中；

步骤3.3.4，根据更新后的密度界面反演计算模型梯度值，更新迭代计算的搜索方向和迭代计算次数后，返回步骤3.3.2中继续进行迭代计算；

更新后迭代计算的搜索方向为：

式中，p_k+1为更新后的迭代计算搜索方向；

步骤3.3.5，结束迭代计算，输出修正后的密度界面深度h_i+1以及密度界面深度h_i+1所对应的密度界面。

本发明所带来的有益技术效果：

本发明提出了一种用于任意变密度界面的反演方法，本发明方法不同于传统密度界面反演方法必须将界面密度变化归纳为函数并推导正演公式后，基于正演推导进行密度界面反演。本发明方法以三维空间网格为基础，将三维密度界面模型在三维空间网格内划分为多个三维棱柱体，通过将各三维棱柱体赋予单独的密度属性，既将三维空间网格内界面的密度变化视为基于函数的空间分布，又能够基于已知的密度属性信息约束插值的空间分布，实现对任意变密度界面的反演，且适用于密度界面出现穿插的情况，有利于准确获取复杂地质构造情况下的密度界面分布，为复杂地质构造的勘探开发提供了依据。

本发明基于三维重力快速反演方法，提高了三维密度界面反演的计算效率，反演速度快且计算成本低，有利于重力密度界面的批量快速反演，适于在复杂地质构造解释过程中推广应用。

附图说明

图1为本发明基于三维棱柱体模型表示密度界面的示意图。

图2为本发明中雅可比矩阵的示意图。

图3为本实施例三维密度界面模型中各层界面的示意图。

图4为本实施例三维密度界面模型中各密度模型的示意图。

图5为本实施例中密度界面处的密度分布。图5中(a)为界面1处的密度分布，图5中(b)为界面2处的密度分布。

图6为本实施例中重力正演所确定重力场的示意图。

图7为本实施例的迭代反演结果。图7中(a)为三维密度界面模型第1次修正后的重力正演场，图7中(b)为三维密度界面模型第1次修正后的重力正演场，图7中(c)为第1次修正后的密度界面，图7中(d)为第1次修正后密度界面与理论界面之间的差值，图7中(e)为三维密度界面模型第3次修正后的重力正演场，图7中(f)为三维密度界面模型第3次修正后的重力正演场，图7中(g)为第3次修正后的密度界面，图7中(h)为第3次修正后密度界面与理论界面之间的差值，图7中(i)为三维密度界面模型第10次修正后的重力正演场，图7中(j)为三维密度界面模型第10次修正后的重力正演场，图7中(k)为第10次修正后的密度界面，图7中(l)为第10次修正后密度界面与理论界面之间的差值，图7中(m)为三维密度界面模型第17次修正后的重力正演场，图7中(n)为三维密度界面模型第17次修正后的重力正演场，图7中(o)为第17次修正后的密度界面，图7中(p)为第17次修正后密度界面与理论界面之间的差值。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

实施例1

本发明提出了一种用于任意变密度界面的反演方法，具体包括如下步骤：

步骤1，基于研究区的地震资料和测井资料获取研究区的密度属性信息，在三维空间网格中构建三维密度界面模型，根据研究区的密度属性信息设置密度界面上方区域和下方区域的密度属性，并利用三维棱柱体正演确定研究区三维密度界面模型的重力场。

由于三维棱柱体模型中密度属性为已知信息，密度界面深度为未知量，因此重力正演是密度界面的函数，利用三维棱柱体正演确定研究区的重力场，如图1所示，确定重力场如公式(1)所示：

g＝g(h) (1)

式中，g为地表重力值，h为密度界面深度。

步骤2，基于研究区三维密度界面模型的重力场，确定密度界面正则化反演的目标函数，建立密度界面反演计算模型。

根据研究区的重力场，确定密度界面正则化反演的目标函数，如公式(2)所示：

式中，为密度界面正则化反演的目标函数，W_d为重力数据加权矩阵，g^obs为重力观测值，g(h)为重力场函数，λ为正则化参数，W_h为密度界面加权矩阵，h_pre为密度界面参考数据。

由于重力场函数g(h)为非线性函数，无法直接展开为线性函数，所以对重力场函数g(h)进行泰勒级数展开并保留一次项，得到近似线性函数：

其中，采用离散三维网格重力反演时雅可比矩阵A_i中元素A_jk的计算公式为：

式中，i为重力反演的计算次数，h_i为第i次重力反演计算所对应的密度界面深度，A_i为密度界面深度h_i所对应的雅可比矩阵，g_j为第j个观测点，h_k为第k个界面深度点，G_jk为界面深度点h_k所在三维棱柱体模型对观测点g_j的正演系数，如图2所示，即正演系数G_jk为单位密度正演值，m_k为三维棱柱体模型的密度值，Δh_k为三维棱柱体模型的纵向长度。

将泰勒级数展开后的重力场函数g(h)带入到密度界面正则化反演的目标函数中，并对密度界面深度h进行求导，得到：

式中，T为转置矩阵。

将密度界面正则化反演的目标函数的取值设置为零，得到密度界面反演迭代公式为：

将公式(6)两遍同时加上整理得到：

通过对公式(7)的左侧矩阵取逆，整理得到密度界面反演计算模型，如公式(8)所示：

步骤3，基于密度界面反演计算模型进行密度界面反演，利用线性共轭梯度算法求解密度界面反演计算模型，通过对密度界面进行多次迭代修正，反演得到研究区的密度界面和密度模型，具体包括以下步骤：

步骤3.1，反演参数初始化；

输入重力观测数据g^obs，选取密度界面参考数据h_pre，设置密度界面深度的初始值h₀、密度界面加权矩阵W_h、重力数据加权矩阵W_d和最大修正次数Nmax，预设修正精度和密度界面反演计算模型中已知界面之间的密度分布。

步骤3.2，采用离散三维网格进行重力反演，根据密度界面深度h_i以及密度界面深度h_i所对应的密度模型，计算密度界面深度h_i所对应的雅克比矩阵A_i。

步骤3.3，基于线性共轭梯度法求解密度界面反演计算模型，得到修正后密度界面深度h_i+1，确定密度界面深度h_i+1所对应的密度界面。

步骤3.3中具体包括以下步骤：

步骤3.3.1，迭代计算参数初始化；

预设密度界面反演计算模型的梯度精度值ε，设置线性共轭梯度计算参数Q和b，设置密度界面反演计算模型的初始梯度r₀、迭代计算的初始搜索方向p₀和迭代更新初值m₀，r₀＝b-Qm₀，p₀＝r₀，m₀＝0.001。

步骤3.3.2，根据密度界面反演计算模型的梯度和迭代搜索方向，计算迭代搜索方向的步长，如公式(9)所示：

式中，α_k为第k次迭代计算时迭代搜索方向的步长，r_k为第k次迭代计算时密度界面反演计算模型的梯度值，p_k为第k次迭代计算的搜索方向。

步骤3.3.3，更新迭代计算结果和密度界面反演计算模型的梯度为：

m_k+1＝m_k+α_kp_k (10)

r_k+1＝r_k-α_kQp_k (11)

式中，m_k+1为更新后的迭代计算结果，m_k为第k次迭代计算的结果，r_k+1为更新后密度界面反演计算模型的梯度值。

步骤3.3.3，若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂小于预设梯度精度值ε，则根据第k次迭代计算的结果m_k，确定修正后密度界面深度h_i+1，如公式(12)所示：

h_i+1＝m_k+h_pre (12)

计算得到修正后密度界面深度h_i+1后，进入步骤3.3.5中。

若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂不小于预设梯度精度值ε，则进入步骤3.3.4中。

步骤3.3.4，根据更新后的密度界面反演计算模型梯度值，更新迭代计算的搜索方向和迭代计算次数，迭代计算搜索方向更新如公式(13)所示：

式中，p_k+1为更新后的迭代计算搜索方向。

更新迭代计算的搜索方向和迭代计算次数k后，返回步骤3.3.2中继续进行迭代计算。

步骤3.3.5，结束迭代计算，输出修正后的密度界面深度h_i+1以及密度界面深度h_i+1所对应的密度界面。

步骤3.4，计算密度界面深度的修正量Δh，当密度界面深度的修正量Δh小于预设的修正精度或修正次数已达到预设的最大修正次数时，则进入步骤3.5中，否则，则更新修正次数i，继续返回步骤3.2中进行修正。

步骤3.5，输出密度界面和密度模型。

实施例2

为了验证本发明提出的一种用于任意变密度界面的反演方法的可行性以及应对复杂模型的能力，本实施例采用实施例1中所述的用于任意变密度界面的反演方法进行密度界面反演，具体包括以下步骤：

步骤1，获取研究区的密度属性信息，建立由双界面控制的三维密度界面模型，如图3和图4所示，三维密度界面模型中密度界面存在穿插的情况。将三维密度界面模型剖分为91×81×50个尺寸为100m×100m×50m三维棱柱体模型，本实施例中三维空间网格中三维棱柱体模型的原点坐标为(50，50，0)，设置纵坐标向下为负。

本实施例中三维密度界面模型由三层密度模型组成，第一层密度模型的空间范围是由地表界面和界面1控制，地表界面为sur0_h(x,y)＝0的平面，界面1的界面深度sur1_h(x,y)为：

sur1_h(x,y)＝-(18+(sin(x/2000×π)+sin(y/1000×π))×3)×40 (14)

第一层模型的密度ρ1(x,y,z)的空间分布为：

ρ1(x,y,z)＝-(x+y-8700)/85000+exp(z/200)/1500-0.65 (15)

式中，(x,y,z)为三维密度界面模型中点的空间位置。

第二层密度模型的空间范围由界面1和界面2控制，界面2的界面深度sur2_h(x,y)为：

第二层模型的密度ρ1(x,y,z)的空间分布为：

ρ2(x,y,z)＝(x+y-8700)/85000+exp(z/333)/4000-0.45 (17)

第三层密度模型的空间范围是由界面2和模型空间底界面控制，第三层模型的密度ρ3(x,y,z)＝0。

从模型中提取密度界面上方的密度分布，如图5所示，从图5中可以看出，两个界面上方的密度分布均为起伏变化，整体呈现西南低、东北高的趋势。由于界面2穿透界面1，导致界面1出现部分空白区域。，面2穿透界面1的区域密度出现突变，整体密度变化很复杂。

利用三维棱柱体正演确定研究区的重力场，研究区的正演重力场如图6所示，重力测点位于地表界面上且与三维棱柱体的水平中心点重合。为了更加符合实际重力测量，减弱模型边界效应的影响，本实施例中将重力正演场各边界均缩减10个测点后用于反演计算。

步骤2，基于研究区的正演重力场，确定密度界面正则化反演的目标函数，建立密度界面反演计算模型为：

步骤3，基于密度界面反演计算模型进行密度界面反演，利用线性共轭梯度算法求解密度界面反演计算模型，通过对密度界面进行多次迭代修正，得到研究区的密度界面和密度模型。

本实施例中密度界面反演过程中，将三维密度界面模型中的界面2设置为目标界面。设置界面2的密度界面深度初始值h₀＝-1383，即将界面2的密度界面深度设置为理论界面深度的平均值。对于理论模型实验，由于密度分布、固定界面等要素已知且准确，目标界面的反演结果唯一，因此本实施例中无需设置参考密度界面，密度界面加权矩阵W_h设置为深度加权矩阵，设置重力数据加权矩阵W_d＝1、最大修正次数Nmax＝17。

根据反演流程，最终经过17次界面修正，得到三维密度界面模型中界面2的反演结果，如图7所示。从迭代反演过程可以看出，随着修正次数的增加，界面反演误差和修正模型的拟合重力场误差逐渐减小，反演结果趋于稳定。最终密度反演界面误差整体很小，仅在反演区域的边界处有出现较大的误差，个别区域最大误差为168m，在反演区域的中部反演误差基本小于50m，小于单个三维棱柱体的纵向长度。

通过上述实验验证了本发明所构建的密度界面反演计算模型能够实现对任意变密度界面的反演，同时也进一步证明了在已知信息丰富且准确的情况下，采用本发明方法进行三维密度界面反演能够取得优异的反演效果。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (6)

1.一种用于任意变密度界面的反演方法，其特征在于，具体包括如下步骤：

步骤1，获取研究区的密度属性信息，在三维空间网格中构建三维密度界面模型，利用三维棱柱体正演确定研究区三维密度界面模型的重力场；

步骤2，基于研究区三维密度界面模型的正演重力场，确定密度界面正则化反演的目标函数，建立密度界面反演计算模型；

步骤3，基于密度界面反演计算模型进行密度界面反演，得到研究区的密度界面和密度模型。

2.根据权利要求1所述的用于任意变密度界面的反演方法，其特征在于，所述步骤1中，利用三维棱柱体模型表示研究区的密度界面，根据研究区的密度属性信息确定密度界面上方区域和下方区域的密度属性信息。

3.根据权利要求1所述的用于任意变密度界面的反演方法，其特征在于，所述步骤1中，基于三维棱柱体模型重力正演确定重力场为：

g＝g(h) (1)

式中，g为地表重力值，h为密度界面深度。

4.根据权利要求3所述的用于任意变密度界面的反演方法，其特征在于，所述步骤2中，根据研究区的重力场，确定密度界面正则化反演的目标函数，如公式(2)所示：

式中，为密度界面正则化反演的目标函数，W_d为重力数据加权矩阵，g^obs为重力观测值，g(h)为重力场函数，λ为正则化参数，W_h为密度界面加权矩阵，h_pre为密度界面参考数据；

由于重力场函数g(h)为非线性函数，对重力场函数g(h)进行泰勒级数展开并保留一次项，得到：

其中，采用离散三维网格重力反演时雅可比矩阵A_i中元素A_jk的计算公式为：

式中，i为重力反演的计算次数，h_i为第i次重力反演计算所对应的密度界面深度，A_i为密度界面深度h_i所对应的雅可比矩阵，g_j为第j个观测点，h_k为第k个界面深度点，G_jk为界面深度点h_k所在三维棱柱体模型对观测点g_j的正演系数，m_k为三维棱柱体模型的密度值，Δh_k为三维棱柱体模型的纵向长度；

将泰勒级数展开后的重力场函数g(h)带入到密度界面正则化反演的目标函数中，并对密度界面深度h进行求导，得到：

式中，T为转置矩阵；

将密度界面正则化反演的目标函数的取值设置为零，得到密度界面反演迭代公式为：

将公式(6)两遍同时加上得到：

通过对公式(7)的左侧矩阵取逆，得到密度界面反演计算模型，如公式(8)所示：

5.根据权利要求4所述的用于任意变密度界面的反演方法，其特征在于，所述步骤3中，利用线性共轭梯度算法求解密度界面反演计算模型，通过对密度界面进行多次迭代修正，反演得到研究区的密度界面和密度模型，具体包括以下步骤：

步骤3.1，反演参数初始化；

输入重力观测数据g^obs，选取密度界面参考数据h_pre，设置密度界面深度的初始值h₀、密度界面加权矩阵W_h、重力数据加权矩阵W_d和最大修正次数Nmax，预设修正精度和密度界面反演计算模型中已知界面之间的密度分布；

步骤3.2，采用离散三维网格进行重力反演，根据密度界面深度h_i以及密度界面深度h_i所对应的密度模型，计算密度界面深度h_i所对应的雅克比矩阵A_i；

步骤3.3，基于线性共轭梯度法求解密度界面反演计算模型，得到修正后密度界面深度h_i+1，确定密度界面深度h_i+1所对应的密度界面；

步骤3.4，计算密度界面深度的修正量Δh，当密度界面深度的修正量Δh小于预设的修正精度或修正次数已达到预设的最大修正次数时，则进入步骤3.5中，否则，则更新修正次数i，继续返回步骤3.2中进行修正；

步骤3.5，输出密度界面和密度模型。

6.根据权利要求5所述的用于任意变密度界面的反演方法，其特征在于，所述步骤3.3中，基于线性共轭梯度法求解密度界面反演计算模型，获取修正后密度界面深度所对应密度界面，具体包括以下步骤：

步骤3.3.1，迭代计算参数初始化；

预设密度界面反演计算模型的梯度精度值ε，设置线性共轭梯度计算参数Q和b，设置密度界面反演计算模型的初始梯度r₀、迭代计算的初始搜索方向p₀和迭代更新初值m₀，r₀＝b-Qm₀，p₀＝r₀，m₀＝0.001；

步骤3.3.2，根据密度界面反演计算模型的梯度和迭代搜索方向，计算迭代搜索方向的步长，如公式(9)所示：

式中，α_k为第k次迭代计算时迭代搜索方向的步长，r_k为第k次迭代计算时密度界面反演计算模型的梯度值，p_k为第k次迭代计算的搜索方向；

步骤3.3.3，更新迭代计算结果和密度界面反演计算模型的梯度为：

m_k+1＝m_k+α_kp_k (10)

r_k+1＝r_k-α_kQp_k (11)

式中，m_k+1为更新后的迭代计算结果，m_k为第k次迭代计算的结果，r_k+1为更新后密度界面反演计算模型的梯度值；

步骤3.3.3，若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂小于预设梯度精度值ε，则根据第k次迭代计算的结果m_k，确定修正后密度界面深度h_i+1，如公式(12)所示：

h_i+1＝m_k+h_pre (12)

计算得到修正后密度界面深度h_i+1后，进入步骤3.3.5中；

若更新后密度界面反演计算模型梯度值的二范数||r_k+1||₂不小于预设梯度精度值ε，则进入步骤3.3.4中；

步骤3.3.4，根据更新后的密度界面反演计算模型梯度值，更新迭代计算的搜索方向和迭代计算次数后，返回步骤3.3.2中继续进行迭代计算；

更新后迭代计算的搜索方向为：

式中，p_k+1为更新后的迭代计算搜索方向；

步骤3.3.5，结束迭代计算，输出修正后的密度界面深度h_i+1以及密度界面深度h_i+1所对应的密度界面。