CN105549106B - 一种重力多界面反演方法 - Google Patents

Abstract

本发明公开了一种重力多界面反演方法，该方法包括以下步骤：获取观测点的信息、重力观测值、初始模型范围、地下介质信息、待反演层位信息、及模型泛函各部分的正则化权重；对地下介质进行网格剖分；计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值；根据重力观测值与重力正演值，计算失配泛函F₁(x)；计算全变差函数F₂(x)；计算模型界面反演偏差函数F₃(x)；根据计算结果、模型泛函各部分的正则化权重λ₂～λ₃与失配泛函F₁(x)构成目标函数：F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)；采用迭代方法使目标函数最小。通过上述方式，本发明能够具有较好的反演精度和效率。

Description

一种重力多界面反演方法

技术领域

本发明涉及地球物理勘探技术领域，特别是涉及一种重力多界面反演方法，其为重力勘探中的密度界面反演方法，可以得到地下多层界面的起伏分布，且具有较高的运算效率。

背景技术

地质勘探主要应用的地球物理勘探方法包括重力勘探、磁法勘探、电法勘探和地震勘探等。重力勘探通过利用重力仪观测地下物质密度差异引起的重力异常以查明地下的地质构造和岩性异常体。相较于地震等勘探方法而言，重力勘探具有横向分辨率好、成本低廉等优点。重力勘探分为三个大的环节：野外重磁异常采集、重磁异常处理和重磁异常解释。

从地质角度，重力反演问题的目标在于寻找、研究或推断金属或非金属矿体和研究地质构造。从地球物理角度，重力反演问题的目标在于确定地质体的几何和物性参数，确定物性分界面的深度和起伏，以及密度的分布等等。

总体来看，当前重力反演方法研究主要可分为两大类。第一类是密度分布研究，着重于物性(即密度)在地下半空间中的分布情况的研究；第二类是界面反演研究，关注点在于明显的密度分界面，一般用于处理密度已知(或密度变化规律已知)的情形下密度分界面的起伏变化情况。

重力界面反演方法的研究起步较早，自20世纪60年代起，Bott(1960)、Cordell etal.(1968)等即对重力异常在空间域的正反演问题进行了初步研究。这个时期的计算方法都是基于空间域的，通过长方体模型进行拟合，逐次逼近消去剩余值，无法利用先验信息。随着计算机技术的发展，快速傅里叶变换技术(FFT)投入使用，Parker(1973)借助FFT，给出了频率域计算地下单层物质重力异常的正演公式，较以往的离散模型方法速度提高了至少一个数量级(冯锐等,1986)，极大地提高了计算效率。随后由Oldenburg(1974)加以发展，建立了反演问题迭代求解方法，后人称为Parker-Oldenburg反演法，目前已经成为了地球物理领域最常用的重力反演方法之一，具有理论严谨、适应性强和快速有效的特点。除此以外，也有其他科研人员提出了一些反演方法，如Spector et al.(1970)的功率谱分析方法、Chavez et al.(1985)的Parker频谱展开线性反演法和Pederson(1977)的广义矩阵反演法等。这些方法里有一些是对Parker-Oldenburg方法的改进或引申，另外一些方法则常常用于某些特定情况，比如功率谱方法适用于对物质界面的大体预估。不过总体而言，以上方法都存在同样的局限性：假定地下物质为单一界面模型，如此通过重力异常的反演只能得到单一界面的深度变化情况。针对这一问题，王万银等(1993)给出了双界面反演方法，实质采用了单层地壳双界面模型，方法上与Parker-Oldenburg方法非常类似，需要额外考虑的是“调整深度”的合理选取以加快收敛速率。该方法理论较为简易，也可以引入控制点进行约束(黄松,2010)，但模型仅针对单层物质的顶底界面，应用范围有限，对诸如俯冲带地区“双莫霍”的反演问题无能为力，且只能反演两个界面，而且密度、顶界面、底界面三者的反演进行中，对其中一项的反演总是相当于假定其他两个参数为已知确定值，使用上造成了一定的局限性。为避开这种将未知界面“假定为已知确定值”的情况，一些基于双界面模型发展的方法(冯旭亮等，2014)干脆以地形作为双界面模型的顶界面，实质又退化为了单界面反演问题，且因模型本身受限，关注点转变为沉积盆地基底研究，方法无法用于深层界面反演。其他的一些非线性全局最优化方法(朱自强等，1995；Snopek，2005；于鹏等，2007；李倩等，2010；明圆圆等，2012)采用非线性全局最优化方法，包括模拟退火算法、遗传算法和人工鱼群算法等，这类算法为避免线性/非线性方程组大规模矩阵求解问题而采用了零阶算法，其优点是对初始模型的依赖程度低，但参数取值范围和选取参数存在随机性，采用启发式算法导致每次的反演结果变动很大，且求解效率不高。

Martins et al.(2011)和冯旭亮等(2014)针对沉积盆地基底反演问题，采用过全变差泛函与控制点信息相结合的方式进行研究。但需要指出，Martins et al.(2011)和冯旭亮等(2014)的方法针对的仅是单一界面反演问题，而且反演对象为浅层物质，无法对深层物质进行反演。

综上所述，现有的地球物理勘探方法难以处理更为复杂的多界面反演情况，并且运算效率不高。

发明内容

本发明主要解决的技术问题是提供一种重力多界面反演方法，可以得到地下多层界面的起伏分布，且具有较高的运算效率和反演精度。

为解决上述技术问题，本发明采用的一个技术方案是：提供一种重力多界面反演方法，该方法包括以下步骤：

(1)获取观测点的信息、重力观测值、初始模型范围、地下介质信息(地下介质的相关数据)、待反演层位信息、与模型泛函各部分的正则化权重；

(2)对地下介质进行网格剖分；

(3)计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值；

(4)根据重力观测值与重力正演值，计算失配泛函F₁(x)，其中，x为待求的各柱状体底界面深度；

(5)计算全变差函数F₂(x)；

(6)计算模型界面反演偏差函数F₃(x)；

(7)根据步骤(5)～(6)的计算结果构成模型泛函，并根据模型泛函各部分的正则化权重λ₂～λ₃与失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)；

(8)采用迭代方法使目标函数最小，期间根据迭代的结果与给定的不等约束，不断修改模型以使反演结果趋向真实值。

其中，步骤(3)进一步包括利用长方体重力正演公式计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值。

其中，步骤(1)进一步包括：

获取用于约束反演模型的信息，其包括控制点信息、模型各层上下界信息、失配泛函协方差矩阵以及模型界面反演偏差矩阵，其中，控制点信息包括控制点位置信息；

步骤(4)和步骤(5)之间还包括：步骤(45)存在控制点信息时，计算控制点信息函数F₄(x)；

步骤(7)包括：根据步骤(45)、(5)以及(6)计算结果构成模型泛函，并根据模型泛函各部分的正则化权重λ₂～λ₄与失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)+λ₄F₄(x)。

其中，步骤(8)具体包括：不断重复步骤(3)～(8)，直至达到迭代退出条件为止，最终求得各层界面深度反演结果，其中，迭代退出条件包括迭代次数达到预设的次数以及迭代的结果误差小于预设的误差。

其中，观测点的信息包括观测点数目和各观测点位置(x,y,z)，其中，用于表示x值的X轴和用于表示y值的Y轴位于同一平面，用于表示z值的Z轴向下为正，X轴、Y轴以及Z轴三者构成右手直角坐标系；

重力观测值包括重力异常观测值和/或重力梯度观测值，当进行模型实验研究时，由标准模型计算出重力异常正演值和/或重力梯度正演值，若使用重力梯度观测值，则重力梯度观测值至少包括梯度张量的其中一个分量；

地下介质的相关数据包括地下介质层数、各层介质的密度和底界面初始深度；

待反演层位信息包括要反演的层位编号，层位编号为由浅部到深部排列或自定义的编号。

其中，重力梯度张量的形式为：

其中G为万有引力常数，ρ为密度，v为体积微元，和分别为观测点位置和场源位置，算子▽作用于观测点，在X-Y-Z坐标系下的重力梯度为

其中一个分量的形式为：

其中，(x¹,x²,x³)＝(x,y,z)，为偏微分符号，为在坐标x上的偏微分。

其中，步骤(4)中，失配泛函形式为

其中，f₁(x)的取值如下：

当无重力梯度观测值时：f₁(x)＝grav_obs(x)-grav_cal(x)，

当有重力梯度观测值时：

其中，grav_obs(x)为重力观测值，grav_cal(x)为重力正演值，gravgra_obs(x)为重力梯度观测值，gravgra_cal(x)为重力梯度正演值，C_ddg为失配泛函协方差矩阵，M为观测点数目。

其中，步骤(45)中，控制点信息函数为

其中

f₄(x)＝k-Wx

其中，W是一个B×N的稀疏矩阵，B为控制点数目，N为待求柱状体个数，稀疏矩阵里面元素非0即1，保证在相应控制点i上的柱状体底界面深度x_i趋近于k里的对应元素值，而n₄为归一化因子。

其中，步骤(5)中，全变差函数为

其中

其中，为全变差函数核，L表示x里空间相邻的元素对个数，x_i与x_j空间上相邻；n₂为归一化因子；β为全变差函数参数，为正数，其取值范围为(0,0.5)。

其中，步骤(6)中，模型界面反演偏差函数为

其中

f₃(x)＝x-x_pred

其中，N为待求底界面深度的地下网格个数，x_pred是上一步的反演结果；C_RR为模型界面反演偏差矩阵；n₃为归一化因子。

其中，步骤(7)中，分别调整模型泛函各部分的正则化权重λ₂～λ₄，使其范围是[0,1]；

或者对模型泛函的各个部分采用相同的正则化权重，采用Zhdanov(2011)的自适应方法来确定正则化权重。

其中，步骤(8)中的不等约束包括各层界面深度范围约束0≤Δx≤Δx_max及各层界面大小范围约束s_min≤s≤s_max，采用条件判断或对数界限法应用约束。

其中，对数界限法使用参数的范围是(0,1)。

本发明的有益效果是：区别于现有技术的情况，本发明提供一种重力多界面反演方法，该方法具体为：首先(1)获取观测点的信息、重力观测值、初始模型范围、地下介质的相关数据、待反演层位信息与模型泛函各部分的正则化权重；然后(2)对地下介质进行网格剖分；进而(3)利用长方体重力正演公式计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值；再而(4)根据重力观测值与重力正演值，计算失配泛函F₁(x)，其中，x为待求的各柱状体底界面深度；进一步(5)计算全变差函数F₂(x)；再进一步(6)计算模型界面反演偏差函数F₃(x)；又进一步(7)根据步骤(5)～(6)的计算结果构成模型泛函，并根据模型泛函各部分的正则化权重λ₂～λ₃与失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)；最后(8)采用迭代方法使目标函数最小，期间根据迭代的结果与给定的不等约束，不断修改模型以使反演结果趋向真实值。因此，本发明结合全变差与迭代偏差约束的重力多界面反演方法，具有较好的反演精度和效率。

附图说明

图1为本发明提供的一种重力多界面反演方法的流程图；

图2为第一实施例的真实模型、初始模型及控制点信息；

图3为第一实施例的界面反演结果；

图4为第一实施例的历次迭代失配泛函情况图；

图5为第一实施例的历次迭代模型泛函情况图；

图6为第二实施例的真实模型、初始模型及控制点信息；

图7为第二实施例的界面反演结果；

图8为第二实施例的历次迭代失配泛函情况图；

图9为第二实施例的历次迭代模型泛函情况图。

具体实施方式

请参阅图1，本发明的多界面重力反演方法包括以下步骤：

步骤S1：获取观测点的信息、重力观测值、初始模型范围、地下介质的相关数据、待反演层位信息与模型泛函各部分的正则化权重。

其中，步骤S1还进一步获取用于约束反演模型的信息，其包括控制点信息、模型各层上下界信息、失配泛函协方差矩阵以及模型界面反演偏差矩阵，其中，控制点信息包括控制点位置信息。因此，可以减少反演多解性，增强反演的准确性。

进一步地，步骤S1还可以反演控制信息，其为反演过程中的一些参数信息，如全变差函数参数、对数界限法使用参数、迭代次数上限、迭代停止误差等，所需要参数的数目与给定的模型其他信息数目与迭代方法有关。值得注意的是，前文所述的模型泛函各部分的正则化权重也可以理解为属于反演控制信息。

其中，观测点的信息包括观测点数目和各观测点位置(x,y,z)，其中，用于表示x值的X轴和用于表示y值的Y轴位于同一平面，用于表示z值的Z轴向下为正，X轴、Y轴以及Z轴三者构成右手直角坐标系。因此，本实施例结合了观测点深度位置。现有技术的反演方法通常假定观测点位于同一水平面，但实际勘探中存在地形影响，观测点位置大多起伏不均，导致这类界面反演方法要求导入的重力观测值(或重力梯度观测值)要经过各种改正和曲化平等处理。而引入这类处理有导入误差的可能。本实施例的观测点位置结合深度信息，符合实际观测情况，导入的重力异常值不需要进行广义地形改正(即狭义地形改正+中间层改正)和曲化平等操作，提高了反演工作效率。

重力观测值包括重力异常观测值和/或重力梯度观测值，当进行模型实验研究时，由标准模型计算出重力异常正演值和/或重力梯度正演值，若使用重力梯度观测值，则重力梯度观测值至少包括梯度张量的其中一个分量。也就是说，可以使用重力异常观测值，也可以使用重力梯度观测值，或者重力异常观测值和重力梯度观测值两个都用也可以。重力梯度观测值用梯度张量的任一分量均可，且分量数量不限(用几个分量均可)。

其中，梯度张量的形式为：

其中G为万有引力常数，ρ为密度，v为体积微元，比如在直角坐标系xyz下，dv就相当于dxdydz，也就是x、y和z方向各自一个很微小的长度dx、dy、dz构成的小立方体。和分别为观测点位置和场源位置，算子▽作用于观测点，在X-Y-Z坐标系下的重力梯度为

其中一个分量的形式为：

其中，(x¹,x²,x³)＝(x,y,z)，为偏微分符号，为在坐标x上的偏微分。如公式中k＝1，l＝2时，意味着对的部分，分别就x轴和y轴求取偏微分，这样是一个二阶的偏微分。

地下介质的相关数据包括地下介质层数、各层介质的密度和底界面初始深度。实际地形中，同一层位介质的密度可具有横向不均匀性，且不需由特定函数约束。各层介质的密度值根据引入重力异常数据性质的不同而进行给定，如对于海水层，在使用空间异常时密度给定为1.03g/cm³，使用布格异常时应给定为2.67g/cm³。如需考虑沉积层随深度的密度变化，可采用反演时将沉积层沿深度向细分为若干层位的方法来解决，或在反演之前采用其他方法进行层位剥离或重力改正。

待反演层位信息包括要反演的层位编号，层位编号为由浅部到深部排列或自定义的编号。可选择反演所有层位的底界面，或只反演其中某些层位的底界面，后者情形中假定其他不反演的界面层位深度为真实值。

步骤S2：对地下介质进行网格剖分。进行网格建模时，采用层状介质模型，每一层的底界面是其下层的顶界面，待求取数据为感兴趣层位的底界面深度。所有层位均由模型最左端到最右端延伸。模型界面重合时，认为两层界面重合区域的底界面深度相同，也即这两层地层中下层的厚度为0。

步骤S3：利用长方体重力正演公式计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值。

将各柱状体对观测点的重力效应并进行累加，计算得到重力正演值后，针对实际情况可以决定是否除去正常场。如读入数据为生产实测重力异常数据，则通常应减去正常场。如正常场已知，则直接在各观测点的重力正演值减去正常场；反之，可将各观测点的重力正演值取平均作为正常场值，在重力正演值中减去。如已知正常场的变化趋势，则可基于此求取正常场变化函数，将正常场从各观测点的重力正演值中去除。如读入数据为根据标准模型计算出的重力正演值，则可不考虑正常场问题。

当梯度信息可获取到时，可利用长方体重力梯度正演公式对初始模型结果进行梯度正演计算，并类似地去除正常场，用于后续失配泛函计算。

步骤S4：根据重力观测值与重力正演值，计算失配泛函F₁(x)，其中，x为待求的各柱状体底界面深度。

其中，失配泛函形式为：

其中，f₁(x)的取值如下：

当无重力梯度观测值时：f₁(x)＝grav_obs(x)-grav_cal(x)，

其中，grav_obs(x)为重力观测值，grav_cal(x)为重力正演值，gravgra_obs(x)为重力梯度观测值，gravgra_cal(x)为重力梯度正演值，C_ddg为失配泛函协方差矩阵，M为观测点数目。本实施例中的失配泛函进行了归一化处理，避免了观测点数目不同对反演参数选择及反演效果造成的影响。进一步地，本实施例重力梯度数据的引入提高了本实施例的方法的可扩展性。

其中，C_ddg为对角阵形式，可采用Zhdanov(2002)灵敏度加权法计算或等效地采用观测数据偏差倒数的平方来进行赋值，也可采用单位矩阵代替。

定义J为f₁(x)的Jacobian矩阵，其大小为M×N，其中第i行第j列对应着第i个观测点的重力观测值和第j个柱状体的底界面。因此，推导了多界面情形下的重力反演失配泛函Jacobian矩阵。为简便，以下基于无梯度数据情形进行说明，此时J的形式为

其中，

f_1i(x)＝grav_obs,i(x)-grav_cal,i(x)

Δg_ij代表了第j个柱状体对第i个观测点的重力效应，可表示为

其中

x、y、z右下角标记中i表示该值对应观测点坐标位置，j表示对应柱状体坐标位置，a、b和c为1指代柱状体各坐标轴方向的较小值，为2表示较大值。由上可得

其中

且

这里nx指一层上的柱状体个数。

步骤S5：计算全变差函数F₂(x)。

本步骤之前，还包括存在控制点信息时，计算控制点信息函数F₄(x)。

其中，控制点信息函数为：

其中，f₄(x)＝k-Wx

其中，W是一个B×N的稀疏矩阵，B为控制点数目，N为待求柱状体个数，稀疏矩阵里面元素非0即1，保证在相应控制点i上的柱状体底界面深度x_i趋近于k里的对应元素值，而n₄为归一化因子。本实施例的控制点信息函数进行了归一化处理，避免了控制点数目不同对反演参数选择及反演效果造成的影响。

本步骤中，全变差函数为：

其中，定义为1范数，其形式为：

其中，为全变差函数核，是全变差函数核心部分，L表示x里空间相邻的元素对个数，x_i与x_j空间上相邻；n₂为归一化因子；β为全变差函数参数，为正数，其取值范围为(0,0.5)。因此，本实施例的全变差函数进行了归一化处理，避免了地下网格数不同对反演参数选择及反演效果造成的影响。

步骤S6：计算模型界面反演偏差函数F₃(x)。

本步骤中，模型界面反演偏差函数为

其中，f₃(x)＝x-x_pred；

其中，N为待求底界面深度的地下网格个数，换言之即为待求向量x的元素个数。x_pred是上一步的反演结果；C_RR为模型界面反演偏差矩阵；n₃为归一化因子。

其中，C_RR可采用Zhdanov(2002)灵敏度加权，也可等效地采用当前模型与真实模型的界面深度大致偏差平方的倒数来加权。模型协方差矩阵应在一定的迭代次数范围内保持稳定性，如每5次迭代更新一次模型协方差矩阵，以保证较好的迭代收敛性。

本实施例中，模型界面反演偏差函数的引入相较于Gallardo et al.(2005)的反演初始模型偏差度函数，本实施例使用的模型界面反演偏差函数每次针对于上一步的迭代结果而不是初始模型，从而降低了反演对初始模型的依赖，同时达到了深度加权效果，改善了重力趋肤效应，具有更好的反演效果。同时，本实施例中的模型界面反演偏差函数进行了归一化处理，避免了地下网格数不同对反演参数选择及反演效果造成的影响。

步骤S7：根据步骤(5)～(6)的计算结果构成模型泛函，并根据所述模型泛函各部分的正则化权重λ₂～λ₃与所述失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)。

更具体地，存在控制点信息时，将F₂(x)、F₃(x)以及F₄(x)计算结果构成模型泛函，并根据模型泛函各部分的正则化权重λ₂～λ₄与失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)+λ₄F₄(x)。

本步骤中，分别调整模型泛函各部分的正则化权重λ₂～λ₄，使其范围是[0,1]。

或者对模型泛函的各个部分采用相同的正则化权重，采用Zhdanov(2011)的自适应方法来确定正则化权重。此时取λ₂＝λ₃＝λ₄＝α，第一次迭代时α₀＝0，第一次迭代后取

α_k＝α₁q^k-1,k＝1,2,...,itermax；0＜q＜1

其中itermax为最大迭代次数。q控制迭代速率，可通过实验获得，对反演结果影响不大。

步骤S8：采用迭代方法使所述目标函数最小，期间根据迭代的结果与给定的不等约束，不断修改模型以使反演结果趋向真实值。

本步骤中，具体包括：不断重复步骤S3～S8，直至达到迭代退出条件为止，最终求得各层界面深度反演结果，其中，迭代退出条件包括迭代次数达到预设的次数以及迭代的结果误差小于预设的误差。

其中，不等约束包括各层界面深度范围约束0≤Δx≤Δx_max及各层界面大小范围约束s_min≤s≤s_max，采用条件判断或对数界限法应用约束。对数界限法使用参数的范围是(0,1)。

本步骤中，可采用的迭代方法包括Gauss-Newton法、Marquardt法和LSQR法等。一般地，以含控制点信息的情形为例，迭代公式为

(A₁+λ₂A₂+λ₃A₃+λ₄A₄)h＝-(g₁+λ₂g₂+λ₃g₃+λ₄g₄)

其中

这里对于函数省去了自变量标识。

对于每次迭代的结果，应进行判断和修正，避免出现不合理的情况。这种修正包括两类：

a.当不清楚各层界面深度上下界且地下界面有重合可能时，对于每次迭代结果进行判断，如果出现了下层界面高于上层界面的情况，强制将下层界面值赋为上层界面值；如出现上层界面穿透了下层的已知深度(不反演)的界面，则将上层界面深度强制赋为该已知深度界面值。

b.当已知各层界面上下界且地下界面有重合可能时，可同时结合对数界限法进行控制。采用公式为

A＝A+λ(X^TX),g＝g+λX

这里λ为对数象限法使用系数，一般不超过1，可以通过实验给定。A和g分别为采用类Gauss-Newton方法迭代时的近似二阶导数项和近似一阶导数项(即A₁+λ₂A₂+λ₃A₃+λ₄A₄和g₁+λ₂g₂+λ₃g₃+λ₄g₄)，d和u为各层界面深度下限和上限。

迭代停止条件可包括三类：

a.到达先前指定的迭代次数；

b.目标函数迭代步长与目标函数之商小于先前设定的阈值；

c.重力拟合残差或目标函数残差小于先前设定的阈值。

迭代停止条件包括了绝对残差、相对残差，也包括了迭代次数上限，避免为达到过低的阈值限而进行的不必要迭代。

本发明较原有技术有较大改进，使其完全符合地质勘探科研及生产的需要并产生了好用及实用的效果。使用该方法提高了多界面反演的效率和精度，改善了反演效果，节省了勘探成本，降低了勘探风险。

本发明还结合以上的方法提供了以下两个实际应用的实施例。第一个实施例请一并参阅图1-图5。

第一实施例采用的标准地球物理模型含四层界面，如图2所示，是对俯冲地区的简单模拟，右侧洋壳俯冲至左侧过渡壳之下。实线部分为真实模型，虚线部分为初始模型。模型信息见表1。整个模型底界深度约束在20km。因所有界面均从模型最左端延至模型最右端，因此层位1的底界面实际包括了过渡壳莫霍部分和俯冲洋壳顶部，层位2底界面左侧为20km底界，右侧向上直接与俯冲洋壳顶层重叠，层位3底界面左侧为20km底界，右侧为洋壳莫霍，层位4底界面恒定为20km。20km以下不填充物质。为简便，各层y方向范围定义为-10km至10km。基于标准模型计算重力正演值，作为反演的输入重力数据，不引入梯度信息。控制点用星号表示。观测点在x轴上均匀分布于-300km-400km处，间隔为0.5km，同时满足y＝0,z＝0。正则化权重采用Zhdanov(2011)的自适应方式，选取q＝0.9，全变差函数参数β为0.001。

表1

由于洋壳顶层界面不参与反演，这里将层位1底界面从40km到100km的值“锁住”，其中，底界面在x轴方向上，使之保持初始模型原状态，反演中若出现这部分值改动的情形，强制将其设定回初始模型值，以此避免因这部分值的不必要变动导致的反演减缓或结果偏离。初始模型与真实模型待反演界面的偏差，层位1和层位3大致为0.3km和0.7km左右，因此取反演控制信息设定迭代次数上限为100，当目标函数相对变化(即)小于5％时迭代停止。

其中，本实施例网格剖分时每个网格x向长度为0.5km。待求值为待求层位上各网格底界面深度。本实施例中待求的是层1和层3的底界面，每层有201个网格，共计402个网格，网格编号按从上到下、每层从左到右递增的顺序编号。由此可构成待求向量x＝(x₁,x₂,...,x_201,x₂₀₂,...,x₄₀₂)^T，其中前201个待求值为层位1底界面，后201个待求值为层位3底界面。

其中，本实施例LSQR法迭代，未引入对数界限法限制。经过前文所述的方法和本实施例的实际应用之后请参阅图3-图5，图3为反演结果与真实模型的对比。经过100次迭代得到最终结果。图4为历次迭代重力拟合残差情况，图5为历次迭代控制点残差、全变差函数及界面反演偏差函数情况。可见反演迭代可较快收敛，反演最终重力拟合残差为1.16×10-3mGal，反演精度较好。

第二个实施例请一并参阅图1、图6-图9。

第二实施例采用的标准模型为俯冲带模型，是基于高德章等(2004)对东海地区俯冲带的研究，对俯冲部分所做的简化模型。模型长度350km。模型信息见表2，上地幔部分没有模拟，同时由于该地区俯冲洋壳层1、层2相对于整个模型而言很浅，洋壳俯冲部分建模只涉及了洋壳层3。本实施例只反演层位5和层位6的底界面。为简便，各层y方向范围定义为-10km至10km。基于标准模型计算重力正演值，作为反演的输入重力数据，不引入梯度信息。控制点用星号表示。观测点在x轴上均匀分布于-200km-550km处，间隔为0.5km，同时满足y＝0,z＝0。

模型其他信息采用了控制点信息、待反演层位界面上下界信息和模型界面反演偏差矩阵C_RR。其中待反演层位界面上下界相对于初始模型深度的范围是[-1,1.2]，保证真实模型在此深度范围内。C_RR＝diag(1,1)。反演控制信息设定迭代次数上限为100，当目标函数相对变化(即)小于1％时迭代停止。正则化权重λ₂～λ₄均为0.5。全变差函数参数β为0.001。

表2

图6展示了标准模型、初始模型和采用的控制点，其中，实线为标准模型，虚线为初始模型，星号为控制点。可见初始模型与真实模型的界面起伏偏差大约在±1km左右。为简便，各层y方向范围定义为-10km至10km。

本实施例中，网格间距0.5km，观测点在x轴上均匀分布于-200km..550km处，间隔为0.5km。待求值为待求层位上各网格底界面深度。本例中待求的是层5和层6的底界面，每层有701个网格，共计1402个网格，网格编号按从上到下、每层从左到右递增的顺序编号。由此可构成待求向量x＝(x₁,x₂,...,x_701,x₇₀₂,...,x₁₄₀₂)^T，其中前701个待求值为层位5，后701个待求值为层位6。

其中，本实施例采用LSQR法迭代。由于界面存在相交情形，需要采用对数界限法限制。经测试，λ＝0.2时效果较好。

经过前文所述的方法和本实施例的实际应用之后请参阅图7-图9。图7为反演结果与真实模型的对比，图8为历次迭代重力拟合残差情况，图9为历次迭代控制点残差、全变差函数及界面反演偏差函数情况，从中可见反演收敛速度较快。经过9次迭代收敛，重力拟合残差1.23mGal。

上述数据处理的过程，本专业技术领域分析人员能熟练完成。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (12)

1.一种重力多界面反演方法，其特征在于，该方法包括以下步骤：

(1)获取观测点的信息、重力观测值、初始模型范围、地下介质信息、待反演层位信息、及模型泛函各部分的正则化权重；

(2)对地下介质进行网格剖分；

(3)计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值；

(4)根据所述重力观测值与重力正演值，计算失配泛函F₁(x)，其中，x为待求的各柱状体底界面深度；

(5)计算全变差函数F₂(x)；

(6)计算模型界面反演偏差函数F₃(x)；

(7)根据步骤(5)～(6)的计算结果构成模型泛函，并根据所述模型泛函各部分的正则化权重λ₂～λ₃与所述失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)；

(8)采用迭代方法使所述目标函数最小，期间根据迭代的结果与给定的不等约束，不断修改模型以使反演结果趋向真实值。

2.根据权利要求1所述的方法，所述步骤(3)为利用长方体重力正演公式计算初始模型网格剖分后的各柱状体对观测点的重力效应并进行累加，得到重力正演值。

3.根据权利要求1或2所述的方法，其特征在于，所述步骤(1)进一步包括：

获取用于约束反演模型的信息，其包括控制点信息、模型各层上下界信息、失配泛函协方差矩阵以及模型界面反演偏差矩阵，其中，所述控制点信息包括控制点位置信息；

所述步骤(4)和步骤(5)之间还包括：步骤(45)存在控制点信息时，计算控制点信息函数F₄(x)；

所述步骤(7)包括：根据步骤(45)、(5)以及(6)计算结果构成模型泛函，并根据所述模型泛函各部分的正则化权重λ₂～λ₄与所述失配泛函F₁(x)一起构成目标函数：

F(x)＝F₁(x)+λ₂F₂(x)+λ₃F₃(x)+λ₄F₄(x)。

4.根据权利要求2所述的方法，其特征在于，所述步骤(8)具体包括：不断重复步骤(3)～(8)，直至达到迭代退出条件为止，最终求得各层界面深度反演结果，其中，迭代退出条件包括迭代次数达到预设的次数以及迭代的结果误差小于预设的误差。

5.根据权利要求2所述的方法，其特征在于，所述观测点的信息包括观测点数目和各观测点位置(x,y,z)，其中，用于表示x值的X轴和用于表示y值的Y轴位于同一平面，用于表示z值的Z轴向下为正，所述X轴、Y轴以及Z轴三者构成右手直角坐标系；

所述重力观测值包括重力异常观测值和/或重力梯度观测值，当进行模型实验研究时，由标准模型计算出重力异常正演值和/或重力梯度正演值，若使用重力梯度观测值，则重力梯度观测值至少包括梯度张量的其中一个分量；

所述地下介质的相关数据包括地下介质层数、各层介质的密度和底界面初始深度；

所述待反演层位信息包括要反演的层位编号，所述层位编号为由浅部到深部排列或自定义的编号。

6.根据权利要求5所述的方法，其特征在于：所述梯度张量的形式为：

其中，G为万有引力常数，ρ为密度，v为体积微元，和分别为观测点位置和场源位置，算子作用于观测点，在X-Y-Z坐标系下的重力梯度为

其中一个分量的形式为：

其中，(x¹,x²,x³)＝(x,y,z)，为偏微分符号，为在坐标x上的偏微分。

7.根据权利要求6所述的方法，其特征在于，所述步骤(4)中，失配泛函形式为：

其中，f₁(x)的取值如下：

当无重力梯度观测值时：f₁(x)＝grav_obs(x)-grav_cal(x)，

当有重力梯度观测值时：

其中，grav_obs(x)为重力观测值，grav_cal(x)为重力正演值，gravgra_obs(x)为重力梯度观测值，gravgra_cal(x)为重力梯度正演值，C_ddg为失配泛函协方差矩阵，M为观测点数目。

8.根据权利要求3所述的方法，其特征在于，所述步骤(45)中，控制点信息函数为：

其中，f₄(x)＝k-Wx

其中，所述W是一个B×N的稀疏矩阵，B为控制点数目，N为待求柱状体个数，稀疏矩阵里面元素非0即1，保证在相应控制点i上的柱状体底界面深度x_i趋近于k里的对应元素值，而n₄为归一化因子。

9.根据权利要求2所述的方法，其特征在于：所述步骤(5)中，所述全变差函数为：

其中，定义为1范数，其形式为：

其中，为全变差函数核，L表示x里空间相邻的元素对个数，x_i与x_j空间上相邻；n₂为归一化因子；β为全变差函数参数，为正数，其取值范围为(0,0.5)。

10.根据权利要求2所述的方法，其特征在于：所述步骤(6)中，模型界面反演偏差函数为

其中，f₃(x)＝x-x_pred

其中，N为待求底界面深度的地下网格个数，x_pred是上一步的反演结果；C_RR为模型界面反演偏差矩阵；n₃为归一化因子。

11.根据权利要求1所述的方法，其特征在于：所述步骤(8)中的不等约束包括各层界面深度范围约束0≤Δx≤Δx_max及各层界面大小范围约束s_min≤s≤s_max，采用条件判断或对数界限法应用约束。

12.根据权利要求11所述的方法，其特征在于：所述对数界限法使用参数的范围是(0,1)。