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

Abstract

本发明公开了一种基于混合范数和互相关约束的三维重力场反演方法，提出了球坐标系基于深度加权、光滑约束、聚焦约束和拟合残差互相关系数约束的反演方法。使得传统球坐标系下大型三维重力场反演深度分辨率低和多解性强的问题得到解决，并将反演迭代次数和计算时间相对于传统方法减少一半，避免了传统反演方法中大型雅克比矩阵的存储，也避免了求解大型线性方程组的问题，有利于大尺度精细剖分下重磁位场高分辨率反演计算。本发明所提出的反演方法具有较高的普适性，既可用于油气资源勘查，水文环境监测，又可用于卫星重力场大尺度三维密度成像、研究地球、月球或其它类地行星壳幔结构、构造运动及动力学机制等。

Description

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

技术领域

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

背景技术

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

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

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

发明内容

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

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

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

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

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

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

和M_λ段，共计

个球壳棱柱体单元，三个方向上剖分间隔分别为Δr、

和Δλ；重力场观测面是距离三维场源一定高度的曲面，观测点个数为

且观测点中心位置与场球壳棱柱体单元中心点位置一一对应；

S2、将球坐标系下三维场源球壳棱柱体单元在观测点上所产生的重力异常分量g_z和重力梯度分量g_zz用矩阵与向量乘积的形式来表示，具体如下：

其中，K_z是球坐标系下三维重力异常分量正演核矩阵，K_zz是球坐标系下三维重力梯度分量正演核矩阵，ρ是三维球壳棱柱体单元的密度向量；

S3、计算球坐标系下三维重力异常分量正演核矩阵K_z和球坐标系下三维重力梯度分量正演核矩阵K_zz，计算公式如下：

其中，r、

和λ分别是观测点的高度分量、纬度分量和经度分量；r′、

和λ′分别是三维场源球壳棱柱体单元的深度分量、纬度分量和经度分量；G是重力常数；

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

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

其中，Φ(ρ)为反演目标函数，Φ_d为数据目标函数，Φ_m为混合范数模型目标函数，Φ_ms为光滑模型目标函数，Φ_mf为聚焦模型目标函数，W_d为观测数据误差矩阵，W_s为光滑模型加权矩阵，W_f为聚焦模型加权矩阵，W_f,i为聚焦模型加权矩阵W_f的第i列向量，β为正则化因子，用于平衡数据目标函数和混合范数模型目标函数之间的权重，α为光滑约束与聚焦约束之间权重因子，ρ^ref为参考密度模型向量，ρ_i为第i个球壳棱柱体单元密度值。

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

S5、对反演目标函数最小化求极值，即求ρ_i的偏导数，并令其等于0，具体如下：

其中，上标T表示矩阵的转置，sign是一个符号函数，当W_f,iρ_i>0，sign(W_f,iρ_i)＝1，当W_f,iρ_i<0，sign(W_f,iρ_i)＝–1；

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

ρ＝(ρ_-i+e_iρ_i) (6)

其中，e_i表示第i个基向量，规模为1行1列，值为1，ρ_–i为将第i个密度ρ_i赋值为0后的向量，表示为：

ρ_-i＝(ρ₁,…,ρ_i-1,0,ρ_i+1,…,ρ_M)^T (7)

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

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

S8、根据第i个球壳棱柱体单元的密度ρ_i的迭代形式，遍历所有球壳棱柱体单元，令i＝i+1，进入下一次迭代，直至满足迭代终止条件。

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

其中ε是迭代终止标准，取值10^-5，max为求两个数最大值的函数。

优选的，在步骤S4中，所述观测数据误差矩阵W_d表达如下：

其中，σ_z,i和σ_zz,i分别为第i个观测重力异常g_z,i和观测重力梯度g_zz,i的误差标准差，符号diag表示对角线矩阵。

优选的，在步骤S4中，所述光滑模型加权矩阵W_s由光滑矩阵D_s、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω构成，具体表示如下：

W_s＝D_s(D_p+D_ω) (11)

所述光滑矩阵D_s表示为：

所述深度加权矩阵D_p表示为：

D_p＝diag{D_p,1,…,D_p,i,…,D_p,M}

所述拟合残差互相关系数矩阵D_ω表示为：

其中，Δg_z,j和Δg_zz,j分别是第j个观测上重力异常和重力梯度分量的拟合残差，K_z(j,i)和K_zz(j,i)分别表示第i个球壳棱柱体单元在第j个观测点所产生的重力异常和重力梯度分量核矩阵元素，ω_i表示拟合残差互相关系数矩阵D_ω第i个对角线元素的倒数。

优选的，在步骤S4中，所述聚焦模型加权矩阵W_f由聚焦矩阵D_f、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω构成，具体表示如下：

W_f＝D_f(D_p+D_ω) (15)

所述聚焦矩阵D_f表达式如下：

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

本发明的有益效果是：

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

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

3)反演过程中无需存储大型核矩阵K_z和K_zz，也无需存储光滑和聚焦约束矩阵W_s和W_f，只需将它们分解成M个列向量，其中M为球壳棱柱体剖分个数，有效解决了高分辨率反演的存储难题；

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_r、

和M_λ段，共计

个球壳棱柱体单元，三个方向上剖分间隔分别为Δr、

和Δλ；重力场观测面是距离三维场源一定高度的曲面，观测点个数为

且观测点中心位置与场球壳棱柱体单元中心点位置一一对应；

1.2、将球坐标系下三维场源球壳棱柱体单元在观测点上所产生的重力异常分量g_z和重力梯度分量g_zz用矩阵与向量乘积的形式来表示，具体如下：

其中，K_z是球坐标系下三维重力异常分量正演核矩阵，K_zz是球坐标系下三维重力梯度分量正演核矩阵，矩阵大小为N*M，ρ是三维球壳棱柱体单元的密度向量，向量大小为M；g_z和g_zz分别是观测点上正演计算所得重力异常和重力梯度向量，向量大小为N。

1.3、在整个反演过程中，只需计算和存储核矩阵K_z和K_zz一次，因此提前计算并保存。

计算球坐标系下三维重力异常分量正演核矩阵K_z和球坐标系下三维重力梯度分量正演核矩阵K_zz，计算公式如下：

其中，r、

和λ分别是观测点的高度分量、纬度分量和经度分量；r′、

和λ′分别是三维场源球壳棱柱体单元的深度分量、纬度分量和经度分量；G是重力常数；

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

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

其中，Φ(ρ)为反演目标函数，Φ_d为数据目标函数，Φ_m为混合范数模型目标函数，Φ_ms为光滑模型目标函数，Φ_mf为聚焦模型目标函数，W_d为观测数据误差矩阵，W_s为光滑模型加权矩阵，W_f为聚焦模型加权矩阵，W_f,i为聚焦模型加权矩阵W_f的第i列向量，β为正则化因子，用于平衡数据目标函数和混合范数模型目标函数之间的权重，α为光滑约束与聚焦约束之间权重因子，取值范围[0,1]，ρ^ref为参考密度模型向量，ρ_i为第i个球壳棱柱体单元密度值。

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

进一步的，观测数据误差矩阵W_d表达如下：

其中，σ_z,i和σ_zz,i分别为第i个观测重力异常g_z,i和观测重力梯度g_zz,i的误差标准差，符号diag表示对角线矩阵。

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

进一步的，所述光滑模型加权矩阵W_s由光滑矩阵D_s、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω三部分构成，具体表示如下：

W_s＝D_s(D_p+D_ω) (11)

所述光滑矩阵D_s表示为：

其中，符号

表示偏导数。

所述深度加权矩阵D_p表示为：

其中，D_p,i表示深度加权矩阵D_p的第i个对角线元素，K_z,i表示核矩阵K _z的第i列元素所组成的数组，K_zz,i表示核矩阵K _zz的第i列元素所组成的数组。

所述拟合残差互相关系数矩阵D_ω表示为：

其中，Δg_z,j和Δg_zz,j分别是第j个观测上重力异常和重力梯度分量的拟合残差，K_z(j,i)和K_zz(j,i)分别表示第i个球壳棱柱体单元在第j个观测点所产生的重力异常和重力梯度分量核矩阵元素，ω_i表示拟合残差互相关系数矩阵D_ω第i个对角线元素的倒数。

进一步，所述聚焦模型加权矩阵W_f由聚焦矩阵D_f、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω三部分构成，具体表示如下：

W_f＝D_f(D_p+D_ω) (15)

所述聚焦矩阵D_f表达式如下：

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

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

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

而本发明采用新的最优化方法：对反演目标函数最小化求极值，即求ρ_i的偏导数，并令其等于0，具体如下：

其中，上标T表示矩阵的转置，sign是一个符号函数，当W_f,iρ_i>0，sign(W_f,iρ_i)＝1，当W_f,iρ_i<0，sign(W_f,iρ_i)＝–1；K^T _z,i是核矩阵K_z是第i列向量的转置，K^T _zz,i是核矩阵K_zz是第i列向量的转置，W^T _s,i为光滑模型加权矩阵W_s的第i列向量的转置。

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

ρ＝(ρ_-i+e_iρ_i) (6)

其中，e_i表示第i个基向量，规模为1行1列，值为1，ρ_–i为将第i个密度ρi赋值为0后的向量，表示为：

ρ_-i＝(ρ₁,…,ρ_i-1,0,ρ_i+1,…,ρ_M)^T (7)

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

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

2.4、根据第i个球壳棱柱体单元的密度ρ_i的迭代形式，遍历所有球壳棱柱体单元，令i＝i+1，进入下一次迭代，直至满足迭代终止条件。

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

其中ε是迭代终止标准，取值10^-5，max为求两个数最大值的函数。

在上述球坐标系下基于混合范数和互相关约束的三维重力场高精度反演方法中，不需要存储大型矩阵W_s和W_f，只需存储其对角线元素，不需要求解大型线性方程组，只需对单个球壳棱柱体单元密度值进行更新迭代即可，让于大尺度高分辨率反演成为可能。

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

步骤(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³，右边台阶密度异常为-1000kg/m³，具体位置及形态如图2所示。将该模型在深度方向、纬度方向和经度方向上分别均匀剖分为10、40和80段，总的模型个数为32000，观测数据个数3200。正演计算该模型在观测面上所产生的重力异常，并加入均值为0标准差为最大值2％的误差，将其作为反演的观测数据，如图3所示。

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

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

表1

反演方法	内存占用(Mb)	迭代次数	计算时间(s)
				本发明所提方法	10.24	619	1613.55
传统反演方法	819.20	1366	3890.79

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

实施例3

上述实施例2通过一个简单倾斜台阶模型证明了本发明所提方法的正确性，为了进一步证明本发明所提出的球坐标系下基于混合范数和互相关约束的三维重力场高精度反演方法的有效性和高效性，接下来本发明给出一个大尺度复杂模型的反演案例。该模型经度范围为-57°到57°，纬度范围为-20°到20°，深度方向半径范围为1638km到1738km，共100km厚，观测面高度为1748km，也即距离该模型外球壳10km高度。异常体可分为6组，涵盖了倾斜台阶、上下正负异常体、上下同性异常体、水平方向正负交替异常体和上下同性水平板，密度异常从-500kg/m³到500kg/m³，具体位置及形态如图5所示。将该模型在深度方向、纬度方向和经度方向上分别均匀剖分为10、40和114段，总的模型个数为45600，观测数据个数4560。正演计算该模型在观测面上所产生的重力异常，并加入均值为0标准差为最大值2％的误差，将其作为反演的观测数据，如图6所示。

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

表2

反演方法	内存占用(Mb)	迭代次数	计算时间(s)
				本发明所提方法	14.59	2613	6671.96
传统反演方法	1663.49	6082	15962.88

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

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

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

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

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

Claims (7)

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

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

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

2.根据权利要求1所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：所述步骤一具体包括如下步骤：

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

和M_λ段，共计

个球壳棱柱体单元，三个方向上剖分间隔分别为Δr、

和Δλ；重力场观测面是距离三维场源一定高度的曲面，观测点个数为

且观测点中心位置与场球壳棱柱体单元中心点位置一一对应；

S2、将球坐标系下三维场源球壳棱柱体单元在观测点上所产生的重力异常分量g_z和重力梯度分量g_zz用矩阵与向量乘积的形式来表示，具体如下：

其中，K_z是球坐标系下三维重力异常分量正演核矩阵，K_zz是球坐标系下三维重力梯度分量正演核矩阵，ρ是三维球壳棱柱体单元的密度向量；

S3、计算球坐标系下三维重力异常分量正演核矩阵K_z和球坐标系下三维重力梯度分量正演核矩阵K_zz，计算公式如下：

其中，r、

和λ分别是观测点的高度分量、纬度分量和经度分量；r′、

和λ′分别是三维场源球壳棱柱体单元的深度分量、纬度分量和经度分量；G是重力常数；

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

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

其中，Φ(ρ)为反演目标函数，Φ_d为数据目标函数，Φ_m为混合范数模型目标函数，Φ_ms为光滑模型目标函数，Φ_mf为聚焦模型目标函数，W_d为观测数据误差矩阵，W_s为光滑模型加权矩阵，W_f为聚焦模型加权矩阵，W_f,i为聚焦模型加权矩阵W_f的第i列向量，β为正则化因子，用于平衡数据目标函数和混合范数模型目标函数之间的权重，α为光滑约束与聚焦约束之间权重因子，ρ^ref为参考密度模型向量，ρ_i为第i个球壳棱柱体单元密度值。

3.根据权利要求1所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：所述步骤二具体包括如下步骤：

S5、对反演目标函数最小化求极值，即求ρ_i的偏导数，并令其等于0，具体如下：

其中，上标T表示矩阵的转置，sign是一个符号函数，当W_f,iρ_i>0，sign(W_f,iρ_i)＝1，当W_f,iρ_i<0，sign(W_f,iρ_i)＝–1；

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

ρ＝(ρ_-i+e_iρ_i) (6)

其中，e_i表示第i个基向量，ρ_–i为将第i个密度ρ_i赋值为0后的向量，表示为：

ρ_-i＝(ρ₁,…,ρ_i-1,0,ρ_i+1,…,ρ_M)^T (7)

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

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

S8、根据第i个球壳棱柱体单元的密度ρ_i的迭代形式，遍历所有球壳棱柱体单元，令i＝i+1，进入下一次迭代，直至满足迭代终止条件。

4.根据权利要求3所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：所述步骤S8中的迭代终止条件为：

其中ε是迭代终止标准，取值10^-5，max为求两个数最大值的函数。

5.根据权利要求2所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：在步骤S4中，所述观测数据误差矩阵W_d表达如下：

其中，σ_z,i和σ_zz,i分别为第i个观测重力异常g_z,i和观测重力梯度g_zz,i的误差标准差，符号diag表示对角线矩阵。

6.根据权利要求2所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：在步骤S4中，所述光滑模型加权矩阵W_s由光滑矩阵D_s、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω构成，具体表示如下：

W_s＝D_s(D_p+D_ω) (11)

所述光滑矩阵D_s表示为：

所述深度加权矩阵D_p表示为：

所述拟合残差互相关系数矩阵D_ω表示为：

D_ω＝diag{1/|ω₁|,1/|ω₂|,…,1/|ω_i|,…,1/|ω_M|}

其中，Δg_z,j和Δg_zz,j分别是第j个观测点上重力异常和重力梯度分量的拟合残差，K_z(j,i)和K_zz(j,i)分别表示第i个球壳棱柱体单元在第j个观测点所产生的重力异常和重力梯度分量核矩阵元素，ω_i表示拟合残差互相关系数矩阵D_ω第i个对角线元素的倒数。

7.根据权利要求1所述的基于混合范数和互相关约束的三维重力场反演方法，其特征在于：在步骤S4中，所述聚焦模型加权矩阵W_f由聚焦矩阵D_f、深度加权矩阵D_p和拟合残差互相关系数矩阵D_ω构成，具体表示如下：

W_f＝D_f(D_p+D_ω) (15)

所述聚焦矩阵D_f表达式如下：

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

Images