0 引言

重力勘探已广泛应用于区域地质调查、矿产勘查、常规油气及页岩气等非常规油气勘探等领域(曾华霖，1999；Nabighian et al., 2005；刘康和郝天珧，2014；柳建新等，2016).三维重力反演作为重力勘探的重要定量解释手段可以得到地下的密度分布，从而为后续的地质解释服务.

众所周知，位场数据反演通常为病态的，解存在非唯一性和不稳定性.为了减小可行解的范围，需要在反演过程中加入不同类型的约束条件及先验信息，包括光滑约束(Li and Oldenburg, 1998)、稀疏约束(Pilkington, 2009; Marchetti et al., 2014)、聚焦反演(Portniaguine and Zhdanov, 1999, 2002；秦朋波和黄大年，2016；高秀鹤和黄大年，2017)、定性约束和定量约束相结合的双重约束机制(刘展等，2011)、相关成像(郭良辉等，2009；孟小红等，2012).尽管约束条件和先验信息的使用对位场数据反演非常重要，但是有些约束条件和先验信息难以与位场数据反演使用的优化方法相融合，从而限制了约束条件和先验信息作用的发挥(耿美霞等，2016).

随着航空重力梯度仪等高精度多分量采集仪器逐渐应用于重力勘探，野外采集的重力数据明显增加，且待反演模型的精细化剖分也使得待反演密度参数增多，这些都给三维重力反演的计算效率带来了挑战.为了提高计算效率，许多学者分别从正演模拟(Caratori Tontini et al., 2009)、灵敏度矩阵存储(Li and Oldenburg, 2003；姚长利等，2003; Martin et al., 2013)、优化算法(Pilkington, 1997; Abedi et al., 2013; Meng et al., 2017; Rezaie et al., 2017)、反问题构建空间(Pilkington, 2009; Marchetti et al., 2014)、并行计算(Moorkamp et al., 2010; Chen et al., 2012; Liu et al., 2012)等方面开展了研究.

近年来，基于近端思想(Rockafellar, 1976)的近端目标函数优化(Proximal Objective Function Optimization, POFO)方法已被应用于地球物理勘探领域，在求解大规模线性地球物理反问题中展现出了良好的应用前景.Dai等(2016)基于POFO方法实现了地震波阻抗反演，与传统的迭代重加权最小二乘算法相比，在保持相同反演精度的同时提高了反演计算效率.Peng等(2017)将POFO方法应用于三维重力数据反演中，取得了较好的反演效果.相比传统的梯度类优化算法，POFO方法有一个明显的优势，即在迭代过程中逐一计算每个未知参数，有较低的计算复杂度和较高的计算效率.

在复杂地质构造背景下，不同岩性单元之间可能会发生物性突变，产生尖锐边界.考虑到POFO方法容易与以求和形式表示的约束相结合，本文将能够以求和形式表示的长尾巴柯西分布加入到反演目标函数中，对待反演模型参数施加稀疏约束，产生块状效果.在Peng等(2017)的工作基础上，本文将POFO方法拓展到快速近端目标函数(Fast Proximal Objective Function, FPOF)优化方法以加速收敛，实现了基于柯西分布约束和FPOF优化的三维重力反演方法.

1 方法原理 1.1 目标函数构建

对于密度物性分布反演，将待反演区域剖分成若干个小长方体单元，每个剖分单元赋予一密度差.根据重力正演公式(Blakely, 1995)，观测点i处的重力异常为

(1)

其中M、N分别为模型剖分单元数和重力数据观测点数，m_j为剖分单元j内的剩余密度，G_ij为核函数，与剖分单元和观测点的位置有关，d_i为观测点i处的重力观测值.将式(1)写成矩阵形式为

(2)

其中G为灵敏度矩阵，m为模型向量，s为观测数据向量.

本文构建的三维重力反演目标函数为

(3)

其中‖·‖₂为L₂范数.假定不同观测数据的噪音是不相关的，W_d为对角数据加权矩阵，其对角线元素为观测数据噪音的标准差的倒数.W为深度加权矩阵，其具体形式为W=diag{1/z₁, 1/z₂, …, 1/z_M}，其中z_j为第j个剖分单元中心的深度.m_ref为参考模型，包含已知的密度测井曲线、岩芯测量等先验信息.σ为柯西分布参数，控制着解的稀疏程度.λ、γ为正则化参数，分别控制着Tikhonov正则化项和柯西分布项对目标函数的贡献.

1.2 FPOF方法

FPOF优化方法的数学理论基础为数学分析里的Lipschitz连续，具有Lipschitz连续的函数的定义及其相关性质可参考Dai等(2016).注意到式(3)由三项构成，即数据拟合项、模型正则化项和柯西分布约束项，首先令

(4)

由具有Lipschitz连续梯度的凸函数的性质可得

(5)

其中η为矩阵(W_dG)^T(W_dG)的最大特征值.

传统的梯度类优化算法的迭代公式为

(6)

其中α_k为迭代步长，迭代公式(6)与求解下面的二次问题是等价的，

(7)

由式(5)和式(7)可看出，式(7)中的二次函数是f(m)在点m=m_k-1的近似，且.因此，式(7)变为

(8)

经过一定的数学运算，式(8)等价于

(9)

其中

(10)

同理，由式(3)的第二项可推导出

(11)

其中τ为矩阵W^TW的最大特征值，且

(12)

由式(3)中前两项的推导过程可得，最小化式(3)中的目标函数等价于

(13)

将式(13)中的L₂范数写成加和形式，可得

(14)

注意到式(14)中求和号里的每一项均为非负数，因此，在第k次迭代中，最小化式(3)中的目标函数进一步等价于最小化式(14)中的每一项，即

(15)

令x=m_i，a=a_k(i)，b=b_k(i)，ψ=，ξ= ，经过一定的数学运算，式(15)等价于

(16)

其中

(17)

柯西分布约束项的引入使得式(17)为非线性，本文采用重加权策略求解式(16).令式(17)关于变量x的导数为零，得到

(18)

其中

(19)

因此，迭代求解公式为

(20)

其中q^k-1=q(x^k-1).

由式(10)和式(12)可看出，计算a_k和b_k利用的只是前一次迭代结果m_k-1.为了加速收敛，本文利用前两次迭代结果m_k和m_k-1计算a_k和b_k，得到快速近端目标函数(FPOF)优化方法.a_k和b_k的计算公式为

(21)

(22)

其中h_k=，且θ_k=.

由上述推导过程可看出FPOF优化方法的一个明显特点是在每次迭代过程中逐一计算每个未知参数m_k(i)，有较低的计算复杂度和较高的计算效率.

基于FPOF优化方法的三维重力反演流程如下：

(1) 令k=1、h₀=m₀和θ₀=1，分别计算矩阵(W_dG)^T(W_dG)和W^TW的最大特征值η和τ；

(2) 计算a_k=h_k-1+(W_dG)^TW_d(d_obs-Gh_k-1)和b_k=h_k-1+W^TW(m_ref-h_k-1)；

(3) 计算m_k=arg min{η(m_i-a_k(i))²+λτ(m_i-b_k(i))²+γln(1+)}；

(4) 当迭代次数k达到指定的最大迭代次数k_max或满足下面条件时，终止迭代过程，

其中ε为容忍误差.

(5) 计算和更新h_k=m_k+(m_k-m_k-1)；

(6) 令k=k+1，然后返回第(2)步.

值得注意的是，从目标函数中数据拟合项的推导过程可看出，本文提出的方法在地球物理反演中有很好的应用前景，可应用于正演响应算子是线性且与物性参数无关的地球物理反问题中.

1.3 正则化参数的选取

本文构建的目标函数包含两个正则化参数：Tikhonov正则化参数λ和柯西分布正则化参数γ，正则化参数的选取对重磁反演的收敛性、效率及反演结果的可靠性是至关重要的.在整个反演过程中，本文将柯西分布正则化参数γ设置为一固定值，数值试验表明γ取值在区间[10^-5, 10^-4]能产生较好的反演结果.对于Tikhonov正则化参数λ，当λ较大时，模型约束项主导反演过程，而当λ较小时，反演过程主要拟合观测数据，此时有可能拟合观测数据中的噪音，使得反演结果中含有假象.为此，本文采用呈现递减规律变化的Tikhonov正则化参数λ，即

(23)

其中初始值λ₀赋予为一较大值或利用下式由初始模型m₀给定，

(24)

相比L曲线准则和GCV准则，该方法比较简单，有更高的计算效率.

2 理论模型测试

为了验证本文提出方法的有效性，本文设计了含有两个长方体的模型(图 1a)，模型大小为2000 m×1500 m×1000 m，两个长方体的大小均为300 m×300 m×200 m，其中一个异常体的中心位置为(650，750，400)，剩余密度为1.0 g·cm^-3，另一个异常体的中心位置为(1350，750，400)，剩余密度为0.6 g·cm^-3.重力异常观测点数为39拆29，测点间距为50 m.由于实际观测数据不可避免带有误差，在由理论模型产生的重力数据中加入了3%高斯噪音(图 1b).

将模型剖分成40×30×40个小长方体，每个小长方体的大小为50 m×50 m×25 m.图 2为理论模型正演响应的相对拟合误差(百分比)，从图中可看出数据拟合效果比较好.图 3为由本文方法反演得到的反演结果，图中黑色矩形框代表异常体的真实位置.图 3a为深度为400 m处的水平切片图，而图 3b为沿x方向的反演结果剖面图(y=750 m)，可以看出，反演得到的异常体位置与理论模型比较一致，且反演结果也比较聚焦，水平切片上的聚焦效果较为明显.

为了探讨本文提出的FPOF优化方法在计算效率上的优势，本文还对该理论模型进行了基于柯西约束和共轭梯度(CG)算法的三维重力反演，具体的基于柯西约束和CG算法的三维重力反演算法可见附录A.基于CG优化方法的反演结果如图 4所示，可以看出，由CG方法得到的异常体位置与理论模型比较一致，且反演结果也比较聚焦.FPOF方法和CG方法均是在4 GB系统内存、3.2 GHz处理器和32位操作系统上运行，表 1给出了两种反演方法的耗时.从图 3、图 4和表 1可看出，尽管FPOF方法和CG方法能够得到相同精度的反演结果，但是相比CG方法，FPOF方法有较少的运行时，这主要是由于FPOF方法是逐一计算每个模型参数，在每次迭代过程涉及的矩阵-向量乘积运算较少.

3 实际资料应用

现将本文提出的方法应用于我国西部某地区的重力数据.研究区目标层是一套第四系湖相沉积地层，以泥质成分为主，其特有的地质特征使得该目标层富集大量生物气，从而气田分布规律等问题一直是该研究区油气地球物理勘探关注的热点问题.研究区气田埋藏相对较浅(埋深＜2000 m)、储量较大，且目标地层具有构造比较平缓、地层压实作用小及孔隙度大等地质地球物理特点.含气层密度与围岩密度的最大差值为0.37 g·cm^-3，使得气藏密度明显低于围岩密度，在野外采集的重力异常数据中气田位置有明显的重力低异常，这为利用重力勘探直接找油气提供了较好的地球物理基础.图 5为研究区的剩余布格重力异常，由已知的钻井资料知，图中的重力低异常是由地下的生物气藏引起的.

将地下剖分成50×30×50个网格，每个网格大小为400 m×500 m×100 m.由于实测数据的噪音的标准差是未知的，本文设置数据加权矩阵为W_d=I，I为单位矩阵.图 6为反演结果三维立体图，从图中可看出气藏能够较好地反演出来，气藏位置比较清晰，且具有一定的聚焦效果.为了更好地分析反演结果，从反演结果中分别提取深度为1500 m处的水平切片(图 7a)、过图 5中重力低异常中心的沿南北方向(x=11 km)(图 7b)和东西方向(y=6.5 km)(图 7c)的剖面，图中的白色虚线框为气藏在水平方向和垂向上的位置.从图 7a可看出，气藏位置在水平方向上比较聚焦，而从图 7b和7c可看出气藏分布的深度范围为1000~2000 m，中心位置在1500 m左右，这与已知的认识(埋深＜2000 m)比较一致.图 8为观测数据相对拟合误差，可看出观测数据拟合比较好.

为了验证气藏位置外反演结果的正确性，提取图 5中井A位置处反演结果，并与测井曲线对比(图 9)，可看出反演结果和测井数据有相同的变化趋势，即地层密度从浅到深逐渐增加.此外，从图 6和图 7可看出气藏位置外的地层密度也是从浅到深逐渐增加，这也与研究区前期研究结果是一致的.因此，由本文方法得到的反演结果与已知的资料有较好的一致性，证明了该方法的实用性.实际数据反演同样是在4 GB系统内存、3.2 GHz处理器和32位操作系统上运行，耗时为103263 ms.

4 结论

本文提出了一种基于柯西分布约束和FPOF优化方法的三维重力反演方法，该方法逐一求解模型参数，有较低的计算复杂度和较高的计算效率，且在目标函数中加入的长尾巴柯西分布约束可得到稀疏解，有助于产生块状效果.理论模型测试表明，相比基于CG优化算法的三维重力反演方法，本文提出的方法在保持相同反演精度的同时具有更高的计算效率.实际资料应用获得了与已知信息一致的反演结果，证明了该方法的实用性.更重要的是，本文提出的方法可扩展到正演算子是线性且与物性参数无关的地球物理反问题中.

附录A 基于柯西约束和CG算法的三维重力反演算法

将本文的反演目标函数即式(3)写为

(A1)

其中为柯西约束项.

柯西约束项C(m)导数为

(A2)

其中为一对角矩阵.

因此，式(A1)的解为

(A3)

由于矩阵P依赖于模型参数m，因此，采用迭代方式求解式(A3)，广泛应用的求解式(A3)的迭代优化算法是基于梯度下降的CG算法.

致谢  本文在研究过程中与中国石油大学(华东)地球科学与技术学院张繁昌教授关于POFO算法进行了有帮助的讨论，两位匿名评审专家为本文提出了中肯的修改意见，在此一并表示感谢.