一种基于混合范数和互相关约束的三维重力场反演方法

技术领域

本发明涉及地球物理及勘探技术领域，具体涉及一种基于混合范数和互相关约束的三维重力场反演方法。

背景技术

重力场反演是重力勘探资料处理解释中的重要环节之一，可以有效刻画地下三维物性分布，从而推断异常体赋存状况，被广泛应用于矿产资源勘探和岩石圈密度结构研究。利用重力异常反演地下三维密度结构(亦称密度成像)面临着诸多困难。一方面三维重力场反演问题高度欠定，这将导致反演解的不稳定和非唯一；另一方面由于重力场是体积场这一自然属性，导致重力反演的多解性强以及深度分辨率低。

针对三维密度反演非唯一性严重及深度分辨率低的问题，目前行之有效的改进措施主要有两种。一种是最大限度地挖掘观测数据空间本身的潜力，例如在模型目标函数中加入改善“趋肤效应”的深度加权函数、采用不同测度函数、通过计算构造指数作为自约束、相关成像法。另一种措施是多数据联合反演，例如重力异常和梯度张量的联合、重力和磁测数据的联合、重力和地震数据的联合。不同数据间的相互融合能有效弥补重力数据深度分辨能力的不足，提高反演解释的可靠性。然而上述成熟的重力场约束反演方法都是基于直角坐标系下提出的，其在球坐标系下的适用性还未得到广泛验证。其次，当前可用于月球重力场联合反演的数据较少。因此，有必要利用合理数学手段充分挖掘观测数据本身信息。

另一方面，21世纪以来，卫星重力测量技术的发展不仅革命性地推动了地球重力场的研究，也推动了利用重力场探测月球乃至其它类地行星内部结构和动力学研究。随着高阶次月球卫星重力场模型和地形球谐模型的发布，高分辨率重力场数据对传统球坐标系下三维重力场正反演方法提出了挑战。传统的球坐标系下重力场正反演方法已难以使用海量高分辨率数据的基本要求，且对于大尺度薄层球壳模型，其反演深度分辨率低和多解性强的问题愈发严重，亟需开发一套高效高精度反演方法。

发明内容

为解决现有技术中存在的问题，本发明提供了一种基于混合范数和互相关约束的三维重力场反演方法，是一种球坐标系基于深度加权、光滑约束、聚焦约束和拟合残差互相关系数约束的反演方法，避免了传统反演方法中大型雅克比矩阵的存储，也避免了求解大型线性方程组的问题，并将迭代次数和计算时间减小一半，有利于大尺度精细剖分下重磁位场高分辨率反演计算，解决了上述背景技术中提到的问题。

为实现上述目的，本发明提供如下技术方案：一种基于混合范数和互相关约束的三维重力场反演方法，包括如下步骤：

步骤一、在球坐标系下基于光滑和聚焦的混合范数约束构建反演目标函数；

步骤二、对单个单元体密度求偏导，遍历所有球壳棱柱体单元，直至迭代终止。

优选的，所述步骤一具体包括如下步骤：

S1、将球坐标系下三维场源在深度方向、纬度方向和经度方向上分别剖分为等间隔的M

S2、将球坐标系下三维场源球壳棱柱体单元在观测点上所产生的重力异常分量g

其中，K

S3、计算球坐标系下三维重力异常分量正演核矩阵K

其中，r、

l和cosψ为观测点到球壳棱柱体单元内任一点的距离及其方位角的余弦，表达如下：

S4、基于光滑约束和聚焦约束的混合范数模型目标函数，并结合数据目标函数，构建得到反演目标函数，如下：

其中，Φ(ρ)为反演目标函数，Φ

优选的，所述步骤二具体包括如下步骤：

S5、对反演目标函数最小化求极值，即求ρ

其中，上标T表示矩阵的转置，sign是一个符号函数，当W

S6、根据公式(5)，将ρ写成如下形式：

ρ＝(ρ

其中，e

ρ

S7、根据公式(6)所述的密度表示方式，将公式(5)中第i个球壳棱柱体单元的密度ρ

其中，上标(k)表示第k次迭代的结果，上标(k+1)表示第k+1次的迭代结果；

S8、根据第i个球壳棱柱体单元的密度ρ

优选的，所述步骤S8中的迭代终止条件为：

其中ε是迭代终止标准，取值10

优选的，在步骤S4中，所述观测数据误差矩阵W

其中，σ

优选的，在步骤S4中，所述光滑模型加权矩阵W

W

所述光滑矩阵D

所述深度加权矩阵D

D

所述拟合残差互相关系数矩阵D

其中，Δg

优选的，在步骤S4中，所述聚焦模型加权矩阵W

W

所述聚焦矩阵D

其中，符号exp表示自然底数e的指数函数，δ为一个接近0的小数，取值为0.01。

本发明的有益效果是：

1)本发明提出了一种新的球坐标系下三维重力场高效反演方法，通过构建基于光滑、聚焦约束的混合范数正则化反演目标函数，并将深度加权矩阵引入到模型目标函数，有效解决了传统球坐标系下三维重力场反演深度分辨率低的问题；

2)通过重力异常分量与重力梯度分量的联合反演，以及拟合残差互相关系数矩阵的构建，将反演迭代次数减小一半，计算时间减少一半，有效解决了传统球坐标系三维重力场反演多解性强的问题；

3)反演过程中无需存储大型核矩阵K

4)本发明在求反演目标函数最小化的过程中，采用对单个球壳棱柱体单元密度求偏导数，使得避开了传统求解大型线性方程组的问题，仅需对单一单元体密度进行迭代，然后遍历所有单元体，同时也避免了存储大型雅克比矩阵，减少了大量矩阵乘积运算及存储。

附图说明

图1为本发明反演方法步骤流程图；

图2为实施例2中简单台阶模型的深度断面图，展示的是纬度为0.125°所在的深度断面真实密度分布，该台阶由两个不同倾角的台阶组成；

图3为实施例2中简单台阶模型在10km高度观测面所产生的重力异常和重力梯度分量图，(a)为重力异常，(b)为重力梯度；

图4为实施例2中简单台阶模型的反演对比结果图，(a)为本发明方法的反演结果，(b)为传统方法反演结果；

图5为实施例3中大尺度复杂异常体模型的深度断面图，展示的是纬度为0.5°所在深度断面的真实密度分布，该复杂模型由倾斜台阶、正负异常体层叠、正负异常体相间等不同系列组成；

图6为实施例3中大尺度复杂模型在10km高度观测面上所产生的的重力异常和重力梯度分量图，(a)为重力异常，(b)为重力梯度；

图7为实施例3中大尺度复杂模型的反演对比结果图，(a)为真实模型，(b)为本发明方法的反演结果，(c)为传统方法反演结果；

图8为本发明方法技术路线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

一种基于混合范数和互相关约束的三维重力场反演方法，如图1所示，包括如下步骤：

步骤1、在球坐标系下基于光滑和聚焦的混合范数约束构建反演目标函数

1.1、对勘探尺度而言，重力场反演计算常在直角坐标系下进行，然而当研究区域较大或者研究整个行星岩石圈时，行星曲率的影响不得不考虑，因此需要在球坐标系下进行。在球坐标系下，地下三维场源通常被剖分球壳棱柱单元体。

将球坐标系下三维场源在深度方向、纬度方向和经度方向上分别剖分为等间隔的M

1.2、将球坐标系下三维场源球壳棱柱体单元在观测点上所产生的重力异常分量g

其中，K

1.3、在整个反演过程中，只需计算和存储核矩阵K

计算球坐标系下三维重力异常分量正演核矩阵K

其中，r、

l和cosψ为观测点到球壳棱柱体单元内任一点的距离及其方位角的余弦，表达如下：

1.4、正演是反演的核心，根据吉洪诺夫正则化理论，基于光滑约束和聚焦约束的混合范数模型目标函数，并结合数据目标函数，构建得到反演目标函数，如下：

其中，Φ(ρ)为反演目标函数，Φ

在上述球坐标系下基于混合范数的反演目标函数中，观测数据误差将会对反演结果造成影响，因此正确评估每个观测数据的误差水平，对反演结果的好坏至关重要。

进一步的，观测数据误差矩阵W

其中，σ

模型目标函数的作用同样十分重要，其决定这反演结果的总体趋势。由于球坐标系下重力场反演多解性强和深度分辨率低的问题严重，因此本申请提出光滑模型加权矩阵和聚焦模型加权矩阵，使两者的优势互补，以提高反演深度分辨率。

进一步的，所述光滑模型加权矩阵W

W

所述光滑矩阵D

其中，符号

所述深度加权矩阵D

其中，D

所述拟合残差互相关系数矩阵D

其中，Δg

进一步，所述聚焦模型加权矩阵W

W

所述聚焦矩阵D

其中，符号exp表示自然底数e的指数函数，δ为一个接近0的小数，取值为0.01。

步骤2、对单个单元体密度求偏导，遍历所有球壳棱柱体单元，直至迭代终止。

2.1、传统的最优化方法是对整个密度向量ρ求导数并令其等于0，这样最终就得到了一个线性方程组，而求解该大型线性方程组将是一个巨大的问题。

而本发明采用新的最优化方法：对反演目标函数最小化求极值，即求ρ

其中，上标T表示矩阵的转置，sign是一个符号函数，当W

2.2、根据公式(5)，将ρ写成如下形式：

ρ＝(ρ

其中，e

ρ

2.3、根据公式(6)所述的密度表示方式，将公式(5)中第i个球壳棱柱体单元的密度ρ

其中，上标(k)表示第k次迭代的结果，上标(k+1)表示第k+1次的迭代结果。

2.4、根据第i个球壳棱柱体单元的密度ρ

进一步的，所述迭代终止条件为：

其中ε是迭代终止标准，取值10

在上述球坐标系下基于混合范数和互相关约束的三维重力场高精度反演方法中，不需要存储大型矩阵W

利用本发明在上述球坐标系下基于混合范数和互相关约束的三维重力场高精度反演方法时，具体操作步骤为：

步骤(1)：数据读取

首先程序从网格化的grid文件中读取观测重力异常及观测重力梯度数据，并自动获取观测数据个数、起始范围等参数，然后根据三个方向上数据范围和个数，计算出剖分间隔，以及每一个测点的坐标信息等；

步骤(2)：参数设置

然后人为设置观测高度和反演深度的起始范围，观测高度一般是已知的，而反演深度范围一般是未知的，需要根据地质资料以及其它地球物理资料去确定，并给出观测数据的误差均方差；接着设置反演参数，包括混合范数相对权重因子和正则化因子等；

步骤(3)：运行程序

程序按照前述公式进行运算，根据提前设定好的一组正则化因子，画出L曲线，找到曲线上拐点及其所对应的正则化因子大小，并再次反演，得到最终最优的反演结果。

反演验证

为了证明本发明所提出反演方法的正确性，本发明给出一个简单双异常体反演模型算例，该模型包含40*80*10个单元体，观测点数为40*80，异常体为两个倾斜台阶，具体见实施例2。本发明还设计一个大尺度复杂模型，验证本申请所提方法的正确性，具体见实施例3。

实施例2

本实施例使用一个月球尺度的小型台阶模型来验证本发明所提反演方法的正确性和高深度分辨率，具体为两个异常体反演模型算例，该模型包含40*80*10个单元体，观测点数为40*80，异常体为两个倾斜台阶。该模型纬度范围为-5°到5°，经度范围为0°到20°，深度方向半径范围为1638km到1738km，共100km厚，观测面高度为1748km，也即距离该模型外球壳10km高度。异常体为两个倾斜台阶，左边台阶密度异常为1000kg/m

设置混合范数权重因子为0.9，采用L曲线法确定最佳正则化因子。分别对比传统球坐标系下基于光滑约束的重力场三维密度反演与本发明所提出的方法的反演结果，如图4所示。可以看出：传统反演方法根本无法正确看出两个正负密度台阶模型的轮廓，也无法反演出倾斜台阶的倾角和真实位置，整体上呈现一团，且在深度方向上毫无分辨率；而本发明所提出的反演方法所得结果，不仅能正确反演出两倾斜台阶模型的产状，也能真实反映出台阶的倾角和具体位置，且其位置与密度大小与真实模型非常吻合，证明了本申请所提方法的有效性。

此外，本发明所提出的基于混合范数和拟合残差互相关系数约束的球坐标系三维重力场反演方法，相比于传统基于光滑约束的反演方法来说，在内存占用、迭代次数和计算时间方面具有显著优势，简单台阶模型算例使用本发明所提方法与传统方法在内存占用、迭代次数和计算时间上的具体统计如表1所示。

表1

通过表1可以看出，本发明所提方法内存占用小，迭代次数和计算时间为传统方法的一半。

实施例3

上述实施例2通过一个简单倾斜台阶模型证明了本发明所提方法的正确性，为了进一步证明本发明所提出的球坐标系下基于混合范数和互相关约束的三维重力场高精度反演方法的有效性和高效性，接下来本发明给出一个大尺度复杂模型的反演案例。该模型经度范围为-57°到57°，纬度范围为-20°到20°，深度方向半径范围为1638km到1738km，共100km厚，观测面高度为1748km，也即距离该模型外球壳10km高度。异常体可分为6组，涵盖了倾斜台阶、上下正负异常体、上下同性异常体、水平方向正负交替异常体和上下同性水平板，密度异常从-500kg/m

设置混合范数权重因子为0.9，采用L曲线法确定最佳正则化因子。分别对比传统球坐标系下重力场三维密度反演与本发明所提出的方法的反演结果，如图7所示。可以看出：对于该复杂模型，传统反演方法几乎毫无反演能力，无法刻画台阶异常体形态，无法区分正负异常体，在深度方向上毫无分辨能力；而本发明所提出的反演方法，能较好恢复出6组不同异常体的基本形态，台阶的位置和产状与真实模型对应较好，能很好恢复真实异常体所在位置，且台阶异常体浅部的负异常体也能被分辨出来，正负相间的3异常体同样被明显区分开来，展现了本发明所提方法的强大之处，即较高的深度分辨能力。进一步证明了本申请所提方法的正确性，如表2所示。

表2

表2统计了本发明所提方法与传统方法在内存占用、迭代次数和计算时间上的对比，可以看出本发明所提方法具有内存占用小、迭代次数少和计算时间少的优点。

因此，本发明所提方法是一种快速高精度的三维重力场反演方法，提出了球坐标系基于深度加权、光滑约束、聚焦约束和拟合残差互相关系数约束的反演方法。使得传统球坐标系下大型三维重力场反演深度分辨率低和多解性强的问题得到解决，并避免了传统反演方法中大型雅克比矩阵的存储，也避免了求解大型线性方程组的问题，有利于大尺度精细剖分下重磁位场高分辨率反演计算。

本发明技术方案的核心具体有以下两点，本发明方法的技术路线如图8所示，第一、是针对传统球坐标系下三维重力场反演深度分辨率低的问题，提出球坐标系下基于光滑、聚焦混合范数正则化约束和深度加权约束的反演目标函数，通过改变两种不同范数之间的权重因子，以达到光滑和聚焦效果的折中，大幅提高球坐标系重力场反演深度分辨率和准确度。第二、是针对球坐标系下三维重力场反演多解性强的问题，提出重力异常分量和梯度分量联合反演方法以及拟合残差互相关约束方法，通过两种不同数据的相互补充以及互相关系数作为先验信息的加入，以改善球坐标系下三维重力场反演多解性强的问题。

本发明所提出的方法具有较高的普适性，既可用于油气资源勘查，水文环境监测，又可用于卫星重力场大尺度三维密度成像、研究地球、月球或其它类地行星壳幔结构、构造运动及动力学机制等。

尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。