阅读说明：本技术 基于最小地质特征单元约束的重力反演方法 (Gravity inversion method based on minimum geological feature unit constraint ) 是由 宋明水 王金铎 于会臻 王千军 陈学国 毕丽飞 李守济 王有涛 相鹏 王圣柱 于 2019-03-05 设计创作，主要内容包括：本发明提供一种基于最小地质特征单元约束的重力反演方法,包括：步骤1,建立最小地质特征单元模型集合；步骤2,建立基于最小地质特征单元模型约束的重力反演目标函数；步骤3,建立密度初始模型；步骤4,求取重力正演核函数；步骤5,计算最小地质特征单元权重系数；步骤6,基于最小地质特征单元及权重系数恢复模型；步骤7,输出最终反演模型。该基于最小地质特征单元约束的重力反演方法将研究靶区的地质认识通过结构特征提取方式获取最小地质特征单元模型集合,并将其施加至重力反演过程中,以保证反演结果来自于最小地质特征单元模型的组合,从而完善现有模型约束方法,提高解释成果的分辨率及可靠性。(The invention provides a gravity inversion method based on minimum geological feature unit constraint, which comprises the following steps: step 1, establishing a minimum geological feature unit model set; step 2, establishing a gravity inversion target function based on minimum geological feature unit model constraint; step 3, establishing a density initial model; step 4, solving a gravity forward kernel function; step 5, calculating a minimum geological feature unit weight coefficient; step 6, recovering the model based on the minimum geological feature unit and the weight coefficient; and 7, outputting the final inversion model. According to the gravity inversion method based on the minimum geological feature unit constraint, geological knowledge of a research target area is subjected to a structural feature extraction mode to obtain a minimum geological feature unit model set, and the minimum geological feature unit model set is applied to a gravity inversion process to ensure that an inversion result comes from the combination of the minimum geological feature unit models, so that the existing model constraint method is perfected, and the resolution and reliability of an interpretation result are improved.)

基于最小地质特征单元约束的重力反演方法

技术领域

本发明涉及勘探地球物理领域，特别是涉及到一种基于最小地质特征单元约束的重力反演方法。

背景技术

重力勘探作为一种经典地球物理勘探方法，已广泛应用于深部地质探测和矿产、水、石油、地热等资源勘查。近年来，随着重力仪器精度的提高，其应用范围已不再局限于早期的区域勘探，而是渗入到了资源勘查的各个环节。

重力反演是有效恢复地下密度分布、确定勘探目标的重要技术，国内外学者已取得了大量研究成果，应用效果也逐步提高。在众多研究方向中，如何施加更有效的模型约束是重点攻关内容之一。

重力模型约束项可通过限制模型密度分布来缩小求解空间，以增加反演结果可靠性。通常采用以下两类方式来施加模型约束：在已知资料丰富的地区，可采用井约束扩展、限定物性取值范围的方式来对反演结果进行控制；而对于已知资料较少、勘探程度低的地区，采用的方案包括施加模型光滑约束可获取构造连续模型、施加聚焦约束可改善孤立场源体的反演分辨率等技术手段。虽然上述两种手段在一定程度上可降低重力反演的多解性，但是前者太依赖于已知资料的精确度及丰富度，仅能保证测井资料附近的反演结果准确度，而后者又过于宽泛，无法因地制宜地分析研究区内的地质资料、解释成果，丢失了研究区的地质结构特征，反演结果分辨率、可靠性难以得到保证，降低了重力勘探技术的应用水平。因此，对于高精度重力反演来说，具有针对性强、应用灵活的模型约束方法仍然需要开展进一步研究。为此我们发明了一种新的基于最小地质特征单元约束的重力反演方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种完善现有模型约束方法，提高解释成果的分辨率及可靠性的基于最小地质特征单元约束的重力反演方法。

本发明的目的可通过如下技术措施来实现：基于最小地质特征单元约束的重力反演方法，该基于最小地质特征单元约束的重力反演方法包括：步骤1，建立最小地质特征单元模型集合；步骤2，建立基于最小地质特征单元模型约束的重力反演目标函数；步骤3，建立密度初始模型；步骤4，求取重力正演核函数；步骤5，计算最小地质特征单元权重系数；步骤6，基于最小地质特征单元及权重系数恢复模型；步骤7，输出最终反演模型。

本发明的目的还可通过如下技术措施来实现：

在步骤1中，以地质解释模型为基础样本数据，通过人工或自动构建的方式，求解可表征原始样本数据的最小地质特征单元模型集合。

在步骤1中，将最小地质特征单元集合定义为局部空间范围内的地质体特征集合，即工区内的地质模型由该集合的加权组合而成；建立最小地质特征单元模型集合的方式包括以下两种：一是人工构建方式，根据已知认识手动勾绘地质单元，包括不同大小的长方体、球体、弧形、高斯源，地质单元范围不超过待反演区域的范围；二是自动构建方式，在地质模型中采用滑动窗口的方式将其组合为样本数据集合；之后利用基于稀疏约束的求解方法，通过字典学习方法，获得最小地质特征单元模型集合，其优化目标函数如下：

其中，

为地质模型样本集合，由研究区已有地质模型剖面

滑动窗口获取；为待学习的最小地质特征单元；

为对应的稀疏权重；N_x、N_z分别为原始地质模型的宽和高；n_x、n_z分别为所取窗口宽和高，l为地质模型样本个数，k为最小地质特征单元的个数；λ为正则化参数。

采用ADMM(Alternating Direction Method of Multipliers，交替方向乘子法)算法对公式(1)进行求解，获得同时满足稀疏条件的Γ和与之对应的D，之后将D矩阵中不同列数据转换为n_x×n_z大小并赋予至待反演模型网格的中心位置，此时完成了最小地质特征单元集合的构建；为表征不同规模的地质体，采用多尺度窗口搜索地质模型S的方式，以获取更为表征能力更完备的最小地质特征单元集合D。

在步骤2中，反演目标函数构建的原则为不仅保证重力拟合残差最小，而且追求利用尽可能少的最小地质特征单元表征反演模型两个条件，以保证将地质先验约束中的结构特征加入到反演过程，进一步缩小反演解空间。

在步骤2中，重力正演计算公式为：

b＝AM

(2)

其中：为重力观测数据；为重力正演核函数；

为待求解的密度模型向量；其中，q为观测数据个数，p为模型个数；

假设所构建的最小地质特征单元集合可表征待求解密度模型，即最小地质特征单元与对应权重系数向量卷积求和得到地质模型：

其中，*代表二维卷积，

为待求解的最小地质特征单元集合D所对应的稀疏加权系数集合；

公式(3)利用如下矩阵形式表示：

其中，C_D＝[C_D1,C_D2,…,C_Dk]；

为代表第k个最小地质特征单元构建生产的卷积算子矩阵；

代表第k个稀疏加权系数向量；

则重力正演公式(2)重新表示如下：

则重力反演目标函数如下：

其中C_b为观测数据b的协方差矩阵；

设

此时，原重力反演问题变为在最小地质特征单元集合约束下的稀疏权重系数

求解问题。

在步骤5中，求解满足稀疏条件的最小地质特征单元权重系数，以保证用最少的最小地质特征单元来表征反演模型。

在步骤5中，公式(7)中的

项不可微，求解困难；为获得最优稀疏解，采用Majorization-Minimization(优化最小化)框架对公式(7) 进行转换求解，利用一个更容易的优化函数来逐步逼近

推导得出稀疏求解公式如下：

其中，

公式中的

为对

的阈值函数，

为

的一阶导数，即保证

中部分结果小于某一阈值时置0；由于求解稀疏权重系数M的过程中会出现大量0值，但公式(9)中

位于分母，会出现奇异值；为此，将公式中的求逆项

替换为：

此时，结合公式(8)和(10)便可利用迭代循环的方式求得稀疏解

在步骤3中，在笛卡尔坐标系下沿x、z坐标轴分别将模型空间划分成为N_x、N_z个矩形网格单元；间距为△x、△z，初始模型均采用均匀半空间模型。

在步骤4中，重力异常正演公式采用对任意多边形棱柱进行计算。

在步骤6中，对获得的最小地质特征单元权重系数与最小地质特征单元模型进行卷积来获得实际反演的密度模型。

在步骤7中，对恢复的模型进行正演计算，再与实测重力数据进行拟合，根据拟合残差大小、最大迭代次数这些反演终止条件判别是否终止计算，满足条件则输出该密度模型作为最终反演结果，否则返回步骤1。

本发明中的基于最小地质特征单元约束的重力反演方法，将研究靶区的地质认识通过结构特征提取方式获取最小地质特征单元模型集合，并将其施加至重力反演过程中，以保证反演结果来自于最小地质特征单元模型的组合，从而完善现有模型约束方法，提高解释成果的分辨率及可靠性。与现有技术相比，本发明的创新点主要在于：

(1)模型约束来自于局部结构特征，而传统重力反演模型约束要么过于宽泛，对模型整体分布特征进行控制，如采用光滑约束、聚焦约束，难以有效缩小求解空间；要么直接采用硬约束，如密度取值范围约束、测井约束，属于点对点的约束，扩展推广能力低。此为本发明与传统重力反演的重要区别之一。

(2)将原重力反演问题转换为稀疏约束下的最小地质特征单元权重系数求解，再由最小地质特征单元与其卷积求和获得反演结果，追求的是用最少的最小地质特征单元来恢复密度模型，而传统重力反演直接求解密度模型，虽然求解简单，但难以提取能反映密度模型分布的关键特征，反演结果可解释性差。此为本发明与传统重力反演的又一重要的区别。

本发明将研究靶区的地质认识有效地定量化融合至重力反演过程中，通过求解最优的最小地质特征单元模型和与之对应稀疏权重系数的方式来获得最终密度反演模型，与常规重力反演方法相比，本发明着重强调反演结果与先验地质模型结构特征的匹配，反演结果分辨率及可靠性得到了有效提高。

附图说明

图1为本发明的基于最小地质特征单元约束的重力反演方法的一具体实施例的流程图；

图2为本发明的一具体实施例中设计的密度模型及重力异常正演结果图；

图3为本发明的一具体实施例中常规重力反演结果图；

图4为本发明的一具体实施例中已知地质模型样本集图；

图5为本发明的一具体实施例中计算得到的四种最小地质特征单元模型集合图；

图6为本发明的一具体实施例中基于图5(a)最小地质特征单元模型集合所反演得到的稀疏权重系数；

图7为本发明的一具体实施例中基于图5(b)最小地质特征单元模型集合所反演得到的稀疏权重系数；

图8为本发明的一具体实施例中基于图5(c)最小地质特征单元模型集合所反演得到的稀疏权重系数；

图9为本发明的一具体实施例中基于图5(d)最小地质特征单元模型集合所反演得到的稀疏权重系数；

图10为本发明的一具体实施例中为反演得到的密度模型图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

如图1所示，本发明一种基于最小地质特征单元约束的重力反演方法一实施例流程图，该实施例具体步骤如下：

(1)建立最小地质特征单元模型集合。以地质解释模型为基础样本数据，通过人工或自动构建的方式，求解可表征原始样本数据的最小地质特征单元模型集合；将最小地质特征单元集合定义为局部空间范围内的地质体特征集合，即工区内的地质模型可由该集合的加权组合而成。建立最小地质特征单元模型集合的方式包括以下两种：一是人工构建方式。可根据已知认识手动勾绘地质单元，如不同大小的长方体、球体、弧形、高斯源等，一般地质单元范围不超过待反演区域的范围；二是自动构建方式。在地质模型中采用滑动窗口的方式将其组合为样本数据集合，之后利用基于稀疏约束的求解方法，通过字典学习方法，获得最小地质特征单元模型集合，其优化目标函数如下：

其中，

为地质模型样本集合，由研究区已有地质模型剖面

滑动窗口获取；

为待学习的最小地质特征单元；

为对应的稀疏权重。N_x、N_z分别为原始地质模型的宽和高；n_x、n_z分别为所取窗口宽和高(即最小地质特征单元模型大小)，l为地质模型样本个数， k为最小地质特征单元的个数；λ为正则化参数。

采用Stephen Boyd(2011)总结的ADMM(Alternating Direction Method ofMultipliers，交替方向乘子法)算法对公式(1)进行求解，获得同时满足稀疏条件的Γ和与之对应的D，之后将D矩阵中不同列数据转换为n_x×n_z大小并赋予至待反演模型网格(大小为N_x×N_z)的中心位置，此时完成了最小地质特征单元集合的构建。为表征不同规模的地质体，采用多尺度窗口搜索地质模型S的方式，以获取更为表征能力更完备的最小地质特征单元集合D。

(2)建立基于最小地质特征单元模型约束的重力反演目标函数。反演目标函数构建的原则为不仅保证重力拟合残差最小，而且追求利用尽可能少的最小地质特征单元表征反演模型两个条件，以保证将地质先验约束中的结构特征加入到反演过程，进一步缩小反演解空间；重力正演计算公式为：

b＝AM

(2)

其中：

为重力观测数据；为重力正演核函数；

为待求解的密度模型向量。其中，q为观测数据个数，p为模型个数。

假设所构建的最小地质特征单元集合可表征待求解密度模型，即最小地质特征单元与对应权重系数向量卷积求和得到地质模型：

其中，*代表二维卷积，为待求解的最小地质特征单元集合D所对应的稀疏加权系数集合。

公式(3)可利用如下矩阵形式表示：

其中，C_D＝[C_D1,C_D2,…,C_Dk]。

为代表第k个最小地质特征单元构建生产的卷积算子矩阵。

代表第k个稀疏加权系数向量。

则重力正演公式(2)可重新表示如下：

则重力反演目标函数如下：

其中C_b为观测数据b的协方差矩阵。

设

此时，原重力反演问题变为在最小地质特征单元集合约束下的稀疏权重系数

求解问题。

(3)建立密度初始模型；在笛卡尔坐标系下沿x、z坐标轴分别将模型空间划分成为N_x、N_z个矩形网格单元。间距为△x、△z，初始模型均采用均匀半空间模型。

(4)重力正演核函数的求取；重力异常正演公式采用Singh(2002) 提出的对任意多边形棱柱进行计算。

(5)计算最小地质特征单元权重系数。求解满足稀疏条件的最小地质特征单元权重系数，以保证用最少的最小地质特征单元来表征反演模型；公式(7)中的项不可微，求解困难。为获得最优稀疏解，采用 Majorization-Minimization(优化最小化)框架对公式(7)进行转换求解，利用一个更容易的优化函数来逐步逼近推导得出稀疏求解公式如下。

其中，

公式中的

为对

的阈值函数，

为

的一阶导数，即保证

中部分结果小于某一阈值时置0。由于求解稀疏权重系数M的过程中会出现大量0值，但公式(9)中

位于分母，会出现奇异值。为此，将公式中的求逆项替换为：

此时，结合公式(8)和(10)便可利用迭代循环的方式求得稀疏解

(6)基于最小地质特征单元及权重系数恢复模型。对获得的最小地质特征单元权重系数与最小地质特征单元模型进行卷积来获得实际反演的密度模型；在一实施例中，将求解得到的最小地质特征单元权重系数

与最小地质特征单元模型集合D进行卷积并求和，以获得密度模型，具体方式可参照公式(4)。

(7)最终反演模型输出。对恢复的模型进行正演计算，再与实测重力数据进行拟合，根据拟合残差大小、最大迭代次数等其它反演终止条件判别是否终止计算，满足条件则输出该密度模型作为最终反演结果，否则返回步骤(1)。该基于地质特征单元约束的重力反演方法可更好的融入地质先验认识，有效提升重力反演的可靠性。

以下结合具体的实施例对本发明做进一步说明。

实验分析充分考虑到了观测数据噪声、背景密度干扰等影响重力反演计算精度的问题，同时，为验证方法对网格剖分误差的鲁棒性，反演网格相比正演网格在x、y方向各移动了5m、3m。

图2为本发明设计的密度模型及重力异常正演结果图。如图中所示，设计一个二维密度模型组合来验证反演效果，地下半空间剖分x方向范围为50～2950m，间距为100m，网格数为29个；z方向范围为0～290m，间距为 10m，网格数为29个。重力观测点范围为100～2900m，间距为100m，观测点共29个，观测点高点为0m。共有3个长方体密度模型，左上角块体剩余密度0.4g/cm³，右上角块体剩余密度0.4g/cm³，下部剩余密度为-0.2g/cm³。为保证方法的实用性，设置为了均值为0.01高斯密度扰动，并对该密度模型进行低通滤波，将其作为背景密度；同时在重力异常正演结果中加入了为原重力异常最大幅值2％的随机噪声。

图3为常规重力反演结果，可看出反演结果分辨率低，且密度值也与真实模型的密度值有很大差别。

图4为已知地质模型样本集示意图。选择不同块体分布的模型构建地质模型样本集合。图中展示中了其中两个模型示意。

图5为本发明计算得到的最小地质特征单元模型集合。选用滑动窗口 x、z方向网格大小为11×7，间距与模型间距相同。为方便示意，图5(a)～图5(d)分别展示了其中4个最小地质特征单元，可见获得的最小地质特征单元模型集合提取出了已知地质模型中的不同块体结构特征。

图6为本发明重力反演结果图。其中，图6-9分别为反演得到的与图5 (a)～图5(d)种最小地质特征单元模型集合所对应的稀疏权重系数；图10为本发明密度模型反演结果。从图6-9中可看出稀疏权重系数只有在第2、4个块体最小地质模型对应的结果中出现较大幅值，有效的利用了地质结构关键特征，将其与最小地质特征单元集合进行卷积求和有效地恢复出了地下密度模型。相比常规技术，本发明密度反演结果(图10)不仅保持已知地质模型的结构特征、提高了反演分辨率，而且与真实地质模型空间分布与取值范围非常接近，证明了本发明方法的精确性和有效性。

以上所述的具体实施方式，对本发明的目的、技术方案和效果进行了进一步的说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。