彭国民， 刘展

1 齐鲁工业大学(山东省科学院)，山东省计算中心(国家超级计算济南中心)，济南 250014 2 中国石油大学(华东)地球科学与技术学院， 青岛 266580

0 引言

重力勘探已广泛应用于矿产、油气等资源勘探中，重力反演作为重力勘探的重要定量解释手段，能够得到地下密度的空间分布(Nabighian，2005；姚长利等，2007；刘展等，2011；孟小红等，2012；王万银等，2014)，并结合已有的其他资料进行后续的地质和地球物理解释.众所周知，重力异常反演具有多解性和不稳定性，目前的三维重力异常反演方法主要针对反演目标函数中的模型项采用有效的正则化技术来获得可靠、稳定且具有一定特点的地质目标体模型(刘银萍等，2013；高秀鹤和黄大年，2017；李泽林等，2019；周少伟等，2021；李芳等，2021).

从观测数据拟合角度，三维重力反演问题通常构建为一最小二乘优化问题，通过极小化观测数据和由当前模型正演得到的模拟数据之间的误差，得到能够拟合观测数据的地下密度模型.传统的三维重力反演方法假设数据残差(观测数据与模拟数据之差)服从高斯分布(Li and Oldenburg，1998)，需要可靠地估计观测误差.此外，野外采集的重力数据的观测误差可能较大，甚至含有极端异常值.因此，观测误差的大小和可靠估计是传统的三维重力反演方法面临的挑战之一.

尽管基于Tikhonov正则化的三维重力异常反演能够获得稳定解，但是当地下密度异常体存在尖锐边界时，基于L2范数正则化的三维重力异常反演得到的密度异常体的边界比较光滑，这给后续的地质和地球物理解释带来一定的困难.针对此问题，众多学者提出了基于非L2范数正则化的非光滑反演方法.大多数非光滑反演方法采用某一函数作用于物性模型，通过极小化目标函数来获得具有尖锐边界的模型解.Last和Kubik(1983)首先基于最小体积约束思想利用最小支撑函数实现了重力异常聚焦反演.Portniaguine和Zhdanov(1999)把最小支撑函数拓展应用到模型梯度中，开发了基于最小支撑梯度函数的三维重力异常聚焦反演.Farquharson和Oldenburg(1998)提出了基于广义的模型度量的非线性反演方法，能够利用一定的表征模型复杂度的函数产生稀疏解.Sun和Li(2014)、Li 等(2018)基于物性模型的Lp范数实现了重磁异常稀疏反演方法.一些学者分别利用模型的总变差函数(Bertete-Aguirre et al., 2002)、柯西范数(Pilkington, 2009；彭国民等，2018b)、模型指数变化性质(Zhao et al., 2016)、L0拟范数(Meng, 2018)、最小熵函数(Rezaie, 2019)实现了非光滑重磁异常反演问题.此外，在重磁异常反演中可以利用已知的物性信息获得具有尖锐边界的模型解.Krahenbuhl和Li(2006)利用已知的物性信息开发了二进制反演方法来改善成像结果.Liu和Jin(2020)实现了基于改进的模糊c均值聚类算法的重力异常反演方法.

相比于高斯分布，基于Tsallis熵最大化原理推导得出的q-高斯分布已应用于复杂系统和交叉学科领域中，且在拟合实测数据方面具有更灵活的优点(Da Silva et al., 2020).基于最小熵正则化的地球物理反演可以获得异常体与围岩之间的尖锐边界，已被应用于大地电磁和重力数据反演中(Ramos et al., 1999；Rezaie, 2019).为此，本文提出了一种基于q-高斯分布和零阶最小熵正则化的三维重力异常聚焦反演方法，该方法不需要估计数据误差，且对较大数据误差的敏感性较弱，同时能够获得聚焦反演结果，且不需要确定合理的聚焦参数，并将本文方法应用于理论模型和实际资料以验证本文方法的有效性.

1 重力正演

重力正演采用点元法，首先将地下剖分成若干个小长方体单元，每个剖分单元赋予一剩余密度值，然后计算每个剖分单元产生的重力异常，根据重力场的叠加原理，将每个剖分单元产生的重力异常进行叠加，得到观测点位置处的重力异常(Blakely，1995)，即：

Gm=d,

(1)

其中G为重力灵敏度矩阵，m为密度模型向量，d为重力异常数据向量.

2 反演方法

2.1 q-高斯分布和零阶最小熵正则化

Tsallis熵是统计物理学中的一种信息度量方式，能够对具有不规则形状的非可加性系统给出较好的物理解释，其数学定义式为(Da Silva et al., 2020)：

(2)

其中k>0，q为一实参数，p(x)为概率密度函数.在一定的约束条件下通过最大化Tsallis熵能够获得q-高斯分布，该分布的一个重要统计性质是更加接近于数据残差满足的分布特点.q-高斯分布函数可表示为(Da Silva et al., 2020)：

(3)

其中1

基于零阶最小熵正则化思想，本文采用如下函数(Ramos et al., 1999)作为反演目标函数的稳定函数，即：

(4)

(5)

其中：

为一对角最小熵矩阵，起聚焦反演作用.由于参数Q是由模型m和参数δ得到，而参数δ取值为10-15能够得到较好的反演结果，并且具有较好的稳定性，因此，该对角矩阵不需要事先确定合理的聚焦参数.

2.2 目标函数构建

传统的三维重力反演方法假设数据残差服从高斯分布，基于最小模型范数正则化的三维重力光滑反演目标函数为：

(6)

其中‖·‖2为L2范数，dobs为重力观测数据，Wd为数据加权矩阵，当不同观测数据的误差互不相关时，该矩阵为一对角矩阵，其对角线元素为观测噪声的标准差的倒数，mref为参考密度模型，Wm为一对角深度加权矩阵，其对角线元素为1/zj，zj为第j个小长方体单元中心的深度(Peng et al., 2018a)，λ为正则化因子，控制着数据拟合项和模型拟合项对目标函数的相对贡献.由目标函数(6)可看出，传统的三维重力光滑反演需要可靠地估计观测误差以确定数据加权矩阵Wd，此外，基于最小模型范数正则化的光滑反演得到的异常体边界比较平滑.

为了避免确定数据加权矩阵以及获得具有尖锐边界的反演模型，本文构建基于q-高斯分布和零阶最小熵正则化的三维重力聚焦反演目标函数，即：

(7)

2.3 优化求解

(8)

首先从变换后目标函数(8)的梯度角度分析该方法对数据误差的敏感性，目标函数(8)的梯度为：

(9)

下面利用高斯牛顿(GN)算法对变换后的目标函数(8)进行最优化求解，其法方程为：

(10)

(11)

(12)

其中I为单位矩阵.

(13)

由于构建的反演目标函数(7)中涉及到参数q，本文将在后面的理论模型测试和实际资料应用中探讨q值的选择对反演结果的影响，即反演模型的不确定性分析.

2.4 正则化因子选取

正则化因子的选择对于三维重力反演的收敛性及反演模型的可靠性是非常重要的.选择的正则化因子应该能使得反演模型在一定程度上拟合观测数据，且同时有最小的结构复杂度(Farquharson and Oldenburg，2004).本文采用一种呈现递减规律变化的正则化因子选取方法，即：

(14)

3 理论模型测试

3.1 模型1

为了测试本文反演方法的有效性，设计含有两个大小相同的规则异常体的模型(图1)，模型大小为1000 m×600 m×500 m，两个异常体的大小均为200 m×200 m×100 m，其中心位置分别为(300 m,300 m,150 m)、(700 m,300 m,210 m).该模型的背景密度值为0.0 g·cm-3，两个异常体的剩余密度均为1.0 g·cm-3.在地表观测重力数据，其观测点数为39×29=1131，观测点间距在x、y方向上分别为25 m、20 m.由于野外采集的重力数据不可避免地带有观测误差，在生成的正演模拟数据中加入均值为0、标准差分别为模拟数据大小的3%、5%、10%的高斯干扰，得到含有不同强度干扰的重力数据(图2).

图1 理论模型示意图Fig.1 Sketch map of synthetic model

将该模型剖分成40×30×25=30000个剖分单元，每个剖分单元的大小为25 m×20 m×20 m.本文将密度模型的最小值、最大值分别设置为0.0 g·cm-3、1.0 g·cm-3.将本文反演方法应用于这三种带有干扰的重力数据，以探讨本文反演方法对含有不同强度干扰的重力数据的适用性，并与传统的三维重力光滑反演方法进行对比.为了定量评价反演结果，本文采用如下两种统计度量，即模型拟合度mRMS和Pearson相关系数R，其定义式分别为：

(15)

图3为利用光滑反演方法反演含有不同强度干扰重力数据的反演结果，图4—图6分别为使用不同q值反演含有3%、5%、10%干扰重力数据的反演结果，图中水平切片图的深度均为150 m，纵向剖面图均沿着x方向且y=300 m，黑色实线为异常体的真实边界.表1为使用不同方法反演含有不同强度干扰重力数据的反演结果的mRMS和R值.对于传统的重力光滑反演方法，从图3和表1可看出：①尽管异常体的位置能够反演出来，但是异常体的边界不是很清晰，且反演的密度值也不能较好地反映异常体的真实值；②随着干扰强度变大，反演模型的mRMS值变大，R值变小，反演结果变差.对于本文聚焦反演方法，从图4—图6和表1可看出：①当q取值为1.1、1.5、2.0时，随着干扰强度变大，反演模型仍具有比较清晰的异常体边界，且反演密度值也比较接近于真实值；②当q取值为2.5时，尽管本文方法得到的反演结果相对变差，但是也优于传统的重力光滑反演结果；③在该理论模型测试中，当干扰强度为3%时，q取值为1.5反演结果具有最小的mRMS值和最大的R值，q取值为2.0，次之；当干扰强度为5%时，q取值为1.5反演结果具有最小的mRMS值和最大的R值，q取值为1.1，次之；当干扰强度为10%时，q取值为1.1反演结果具有最小的mRMS值和最大的R值，q取值为1.5，次之.

图2 含有不同强度干扰的重力数据 (a) 3%； (b) 5%； (c) 10%.Fig.2 Gravity data with noise of different level

图3 利用光滑反演方法反演含有不同强度干扰重力数据的反演结果 (a) 3%，平面图； (b) 3%，剖面图； (c) 5%，平面图； (d) 5%，剖面图； (e) 10%，平面图； (f) 10%，剖面图.Fig.3 Inversion results from inversion of gravity data with noise of different level using smooth inversion method (a) 3%, horizontal slice; (b) 3%, vertical section; (c) 5%, horizontal slice; (d) 5%, vertical section; (e) 10%, horizontal slice; (f) 10%, vertical section.

图4 使用不同q值反演含有3%干扰重力数据的反演结果 (a) q=1.1，平面图； (b) q=1.1，剖面图； (c) q=1.5，平面图； (d) q=1.5，剖面图； (e) q=2.0，平面图； (f) q=2.0，剖面图； (g) q=2.5，平面图； (h) q=2.5，剖面图.Fig.4 Inversion results from inversion of gravity data with 3% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

图5 使用不同q值反演含有5%干扰重力数据的反演结果 (a) q=1.1，平面图； (b) q=1.1，剖面图； (c) q=1.5，平面图； (d) q=1.5，剖面图； (e) q=2.0，平面图； (f) q=2.0，剖面图； (g) q=2.5，平面图； (h) q=2.5，剖面图.Fig.5 Inversion results from inversion of gravity data with 5% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

图6 使用不同q值反演含有10%干扰重力数据的反演结果 (a) q=1.1，平面图； (b) q=1.1，剖面图； (c) q=1.5，平面图； (d) q=1.5，剖面图； (e) q=2.0，平面图； (f) q=2.0，剖面图； (g) q=2.5，平面图； (h) q=2.5，剖面图.Fig.6 Inversion results from inversion of gravity data with 10% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

3.2 模型2

设计含有两个方向相反、大小不同的倾斜脉状体的模型(图7)，模型大小为600 m×600 m×300 m，背景密度值为0.0 g·cm-3.两个倾斜脉状体的顶部埋深均为地表以下40 m位置处，左边倾斜脉状体的宽度为80 m，且向下延伸100 m，而右边倾斜脉状体的宽度为100 m，且向下延伸120 m，每个倾斜脉状体的剩余密度均为1.0 g·cm-3.在地表观测重力数据，其观测点数为29×29=841，观测点间距在x、y方向上均为20 m.由于野外采集的重力数据不可避免地带有观测误差，在生成的正演模拟数据中加入均值为0、标准差分别为模拟数据大小的3%、5%、10%的高斯干扰，得到含有不同强度干扰的重力数据(图8).

图7 理论模型示意图Fig.7 Sketch map of synthetic model

将该模型剖分成30×30×15=13500个剖分单元，每个剖分单元的大小为20 m×20 m×20 m.将本文反演方法应用于这三种带有干扰的重力数据，以探讨本文反演方法对含有不同强度干扰的重力数据的适用性，并与传统的三维重力光滑反演方法进行对比.

图9为利用光滑反演方法反演含有不同强度干扰重力数据的反演结果，图10—图12分别为使用不同q值反演含有3%、5%、10%干扰重力数据的反演结果，图中水平切片图的深度均为60 m，纵向剖面图均沿着x方向且y=300 m，黑色实线为倾斜脉状体的真实边界.表2为使用不同方法反演含有不同强度干扰重力数据的反演结果的mRMS和R值.对于传统的重力光滑反演方法，从图9和表2可看出：①尽管两个倾斜脉状体的位置和倾斜方向均能够反演出来，但是反演出的倾斜脉状体的边界不是很清晰，且反演的密度值不能较好地反映倾斜脉状体的真实值；②随着干扰强度变大，反演模型的mRMS值变大，R值变小，反演结果变差.对于本文聚焦反演方法，从图10—图12和表2可看出：①随着干扰强度变大，反演模型也具有比较清晰的异常体边界，且反演密度值也比较接近于真实值；②对于干扰强度为10%的重力数据，本文方法仍能较好地得到倾斜脉状体的形态和密度值，反演结果优于传统的重力光滑反演结果；③在该理论模型测试中，q取值为2.0反演结果具有最小的mRMS值和最大的R值，q取值为1.5，次之.

表1 使用不同方法反演含有不同强度干扰 重力数据的反演结果的mRMS和R值Table 1 The values of mRMS and R of inversion results from inversion of gravity data with noise of different level using different inversion methods

表2 使用不同方法反演含有不同强度干扰 重力数据的反演结果的mRMS和R值Table 2 The values of mRMS and R of inversion results from inversion of gravity data with noise of different level using different inversion methods

图8 不同强度噪声的重力数据 (a) 3%； (b) 5%； (c) 10%.Fig.8 Gravity data with different noise level

图9 利用光滑反演方法反演含有不同强度干扰重力数据的反演结果 (a) 3%，平面图； (b) 3%，剖面图； (c) 5%，平面图； (d) 5%，剖面图； (e) 10%，平面图； (f) 10%，剖面图.Fig.9 Inversion results from inversion of gravity data with noise of different level using smooth inversion method (a) 3%, horizontal slice; (b) 3%, vertical section; (c) 5%, horizontal slice; (d) 5%, vertical section; (e) 10%, horizontal slice; (f) 10%, vertical section.

图10 使用不同q值反演含有3%干扰重力数据的反演结果 (a) q=1.1，平面图； (b) q=1.1，剖面图； (c) q=1.5，平面图； (d) q=1.5，剖面图； (e) q=2.0，平面图； (f) q=2.0，剖面图； (g) q=2.5，平面图； (h) q=2.5，剖面图.Fig.10 Inversion results from inversion of gravity data with 3% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

图11 使用本文方法基于不同q值反演含有5%干扰的重力数据的反演结果 (a) q=1.1，平面图； (b) q=1.1，剖面图； (c) q=1.5，平面图； (d) q=1.5，剖面图； (e) q=2.0，平面图； (f) q=2.0，剖面图； (g) q=2.5，平面图； (h) q=2.5，剖面图.Fig.11 Inversion results from inversion of gravity data with 5% noise using this proposed inversion method with different value of q (a) q=1.1, horizontal slice; (b) q=1.1, vertical section; (c) q=1.5, horizontal slice; (d) q=1.5, vertical section; (e) q=2.0, horizontal slice; (f) q=2.0, vertical section; (g) q=2.5, horizontal slice; (h) q=2.5, vertical section.

图13 San Nicolas矿区的剩余重力异常数据 黑色三角为重力观测点.Fig.13 The residual gravity anomaly of San Nicolas deposit Black triangles indicate gravity observation points.

4 实际资料应用

为了验证本文提出方法的实用性，现将本文方法应用于墨西哥萨卡特卡斯州圣尼古拉斯(San Nicolas)矿区的重力数据.San Nicolas矿区上侏罗系到下白垩系的长英质和铁镁质火山岩中发育有火山成因块状硫化物，火山基岩被第三系的火山碎屑角砾岩覆盖，其厚度为50～150 m.San Nicolas矿区块状硫化物矿体的最高密度为3.5 g·cm-3，围岩的密度为2.3～2.7 g·cm-3，因此该矿区的密度差最大可达1.2 g·cm-3(Phillips et al.，2001)，这为利用重力勘探方法直接找块状硫化物矿体提供了物性基础.图13为San Nicolas矿区的剩余重力异常图，由钻井信息知大约中间位置处的高重力异常是由地下的块状硫化物矿体引起的.

将该矿区地下空间离散成43×36×20=30960个剖分单元，每个剖分单元的大小为50 m×50 m×50 m，在反演过程中设置密度差的上、下边界约束分别为1.2 g·cm-3、-0.1 g·cm-3(Phillips et al.，2001).图14—图17分别为q取值为1.1、1.5、2.0、2.5时的反演结果，每个图中的子图(a)为反演结果在深度为300 m位置处的水平切片，子图(b)为反演结果沿着图13中的E1-E2线的纵向切片，子图(c)为反演结果沿着图13中的N1-N2线的纵向切片，子图(d)为由反演结果进行正演得到的预测重力异常数据，图中的黑色实线框代表矿体的真实边界.从这些图中可看出，q取值为1.1(对应反演结果图14)或1.5(对应反演结果图15)时，本文方法能够得到与已知信息比较一致的反演模型.在图14a和15a中，中间的密度异常体为San Nicolas矿体，位于东北方向的异常体是由铁镁质火山岩引起的(Phillips et al.，2001)；在图14b、c和图15b、c中，已知的钻井信息表明，块状硫化物矿体的顶部埋深在地表以下150～220 m，矿体向下延伸到地表以下400～450 m，反演的块状硫化物矿体的位置与已知信息比较一致，且与矿体的真实边界也比较吻合；从图14d和图15d可看出，反演结果也能够较好地拟合观测数据.

5 结论

本文提出了一种基于q-高斯分布和零阶最小熵正则化的三维重力聚焦反演方法，该方法不需要估计数据误差，且对较大数据误差的敏感性较弱，同时能够获得聚焦反演结果，也不需要事先确定合理的聚焦参数.理论模型测试表明，相比于传统的光滑反演方法，本文方法能够较好地刻画地质体边界和密度值，且对含有较强干扰的重力数据具有一定的适用性.墨西哥圣尼古拉斯矿区的实测重力数据应用表明，本文方法反演的硫化物矿体与已知资料具有较好的一致性.理论模型和实际资料测试中的反演模型不确定性分析表明，q取值为1.5能够得到令人满意的反演结果.

图14 q=1.1时的反演结果 (a) 深度为300 m处的水平切片； (b) 沿着图13中E1-E2线的纵向切片； (c) 沿着图13中N1-N2线的纵向切片； (d) 预测的重力异常数据.(b)和(c)中的黑色实线框代表矿体真实边界.Fig.14 Inversion results for q=1.1 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

图15 q=1.5时的反演结果 (a) 深度为300 m处的水平切片； (b) 沿着图13中E1-E2线的纵向切片； (c) 沿着图13中N1-N2线的纵向切片； (d) 预测的重力异常数据.(b)和(c)中的黑色实线框代表矿体真实边界.Fig.15 Inversion results for q=1.5 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

图16 q=2.0时的反演结果 (a) 深度为300 m处的水平切片； (b) 沿着图13中E1-E2线的纵向切片； (c) 沿着图13中N1-N2线的纵向切片； (d) 预测的重力异常数据.(b)和(c)中的黑色实线框代表矿体真实边界.Fig.16 Inversion results for q=2.0 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

图17 q=2.5时的反演结果 (a) 深度为300 m处的水平切片； (b) 沿着图13中E1-E2线的纵向切片； (c) 沿着图13中N1-N2线的纵向切片； (d) 预测的重力异常数据.(b)和(c)中的黑色实线框代表矿体真实边界.Fig.17 Inversion results for q=2.5 (a) The horizontal slice at depth=300 m; (b) The vertical slice along line E1-E2 in theFig.13; (c) The vertical slice along line N1-N2 in theFig.13; (d) The predicted gravity data. The black solid lines in the (b) and (c) denote the true boundaries of the ore body.

致谢感谢加拿大英属哥伦比亚大学GIF组的公开实测重力数据，感谢三位评审专家为本文提出了中肯的修改意见.