地球物理学报  2017, Vol. 60 Issue (4): 1571-1583   PDF    
引用本文
高秀鹤, 黄大年. 2017. 基于共轭梯度算法的重力梯度数据三维聚焦反演研究. 地球物理学报, 60(4): 1571-1583, doi: 10.6038/cjg20170429.  
GAO X H, HUANG D N. 2017. Research on 3D focusing inversion of gravity gradient tensor data based on a conjugate gradient algorithm. Chinese J. Geophys. (in Chinese), 60(4): 1571-1583, doi: 10.6038/cjg20170429.  
基于共轭梯度算法的重力梯度数据三维聚焦反演研究
高秀鹤1, 黄大年1,2     
1. 吉林大学地球探测科学与技术学院, 长春 130000;
2. 国土资源部应用地球物理重点实验室, 长春 130000
收稿日期 2016-09-03, 2017-01-09 收修定稿
基金项目: 国家重点研发计划（2016YFC0600505）与吉林大学研究生创新基金资助项目（2016095）联合资助
作者简介: 高秀鹤, 女, 博士研究生在读, 1991年生, 吉林通化人, 主要从事重力与重力梯度数据反演研究.E-mail:842091616@qq.com
通讯作者: , 男, 1958-2017年, 教授, 广西南宁人, 主要从事移动平台探测数据处理与解释及一体化软件平台开发.E-mail:dnhuang@jlu.edu.cn
摘要: 重力梯度数据相对于传统重力数据，能够更细致、准确地描述地球浅部构造和研究矿产资源分布等信息.本文采用共轭梯度算法，在加权密度域求解重力梯度数据三维聚焦反演最优化问题，以恢复地下三维密度分布，目标函数包括数据不拟合函数和最小支撑稳定函数.首先，在推导目标函数对加权密度的一阶导数时，为了得到更合理的计算公式，我们考虑变加权函数中含有密度变量；此外，本文通过密度上下限约束，改善了传统聚焦反演中聚焦因子选取困难的问题.新算法获得的反演结果，对聚焦因子的选择约束较少，相比传统聚焦算法，能够更容易的获得理想结果.将方法应用于理论模型验证其有效性和正确性，并应用本文方法处理文顿盐丘地区的航空全张量重力梯度数据，得到了与已知地质信息匹配的密度分布，表明本文方法具有处理实际数据的能力.
关键词: 重力梯度数据      共轭梯度      聚焦反演      聚焦因子     
Research on 3D focusing inversion of gravity gradient tensor data based on a conjugate gradient algorithm
GAO Xiu-He1, HUANG Da-Nian1,2     
1. College of GeoExploration Science and Technology, Jilin University, Changchun 130000, China;
2. Key Laboratory of Applied Geophysics, The Ministry of Land and Resources, Changchun 130000, China
Abstract: Gravity gradient data can be used to resolve the earth shallow structure in more detail and study the distribution of mineral resources more accurately compared with traditional gravity data. In order to recover the subsurface 3D density distribution we solve the optimization problem about the focusing inversion of gravity gradient tensor data in a space of the weighted density by the conjugate gradient algorithm. The objective function includes a misfit functional between the observed data and the predicted data and a minimum support stabilizing functional. At first, when deriving the gradients of the objective function, we consider the density variable in the variable weighting matrix to obtain a more reasonable formula. At the same time, by performing penalization of the density we successfully overcome the problem that the inversion result is so sensitive to the focusing factor e. Compared with the conventional focusing algorithm we can obtain the focusing results more easily. To investigate the validity and practicability, we apply the method to synthetic gravity gradient anomalies on typical models. Then we apply the inversion method to gravity gradient data from the Vinton salt dome, Louisiana, USA. The results show a good match to the known geologic information.
Key words: Gravity gradient data      Conjugate gradient      Focusing inversion      Focusing factor     
1 引言

全张量重力梯度 (Full Tensor Gravity Gradiometry, FTG) 仪器的问世和成功使用，使利用重力梯度数据进行地质解释成为可能，与其相关的处理方法和解释技术迅速发展 (耿美霞等，2016；秦朋波和黄大年，2016；叶周润等，2016).地球物理反演问题存在多解性，为了得到唯一解，我们基于一些法则引入约束条件，最常用的是模型光滑度最大约束.Last和Kubik (1983)在重力反演中加入异常体体积最小的约束条件，最小二乘解得到了聚焦的密度分布，从此，体积最小的约束条件被广泛应用于地球物理反演 (Guillen and Menichetti, 1984; Barbosa and Silva, 1994; Silva and Barbosa, 2006; Silva Dias et al., 2009, 2011)，并得到较好的反演图像.Portniaguine和Zhdanov (1999)引入最小梯度约束条件，提出聚焦反演方法，并基于共轭梯度算法迭代求解.Portniaguine和Zhdanov (2002)提出在加权密度参数mw空间利用重加权共轭梯度算法求解紧致反演，再转化成密度m，与直接计算密度m相比，收敛速度更快，共轭梯度算法需要对目标函数求导数，求导过程中将变加权函数We视为常数.Zhdanov (2009)直接反演密度参数m，对目标函数求梯度时，考虑了We中含有变量m，视其为变量求导，进一步发展了重加权共轭梯度算法在反演地球物理数据中的应用.虽然，理论上聚焦因子是一个趋近于0的极小值，但是，利用共轭梯度求解目标函数时，聚焦因子的选取都是非常复杂的问题，太大会导致解不聚焦，太小会使解不稳定.

本文在处理重力梯度数据恢复地下密度分布的过程中，在加权密度参数域求解，并在推导目标函数导数过程中，考虑变加权函数We中含有密度变量m.此外，通过设置密度上下限，使聚焦因子可以选择任意非零极小值，避免了聚焦因子的不合理选择对反演问题的影响.

2 方法理论和计算过程 2.1 方法理论

考虑线性地球物理反演问题，公式为

(1)

其中，M维向量m是地质模型参数，本文为密度向量；N维向量d是观测异常数据，本文为重力梯度数据；N×M维矩阵A是正演模型算子.

由于观测数据中不可避免的噪声干扰，反演问题总是病态的，问题的解不稳定，吉洪诺夫正则化方法作为一种常规的处理手段，有效的解决了这个问题.我们可以把问题 (1) 转化为求解目标函数最优解的问题.目标函数为

(2)

其中ϕd是数据不拟合函数，ϕd是模型目标函数，α代表正则化参数.数据不拟合函数定义为

(3)

模型目标函数为最小支撑稳定函数 (Last and Kubic, 1983)，公式为

(4)

以上公式可以表示为伪二次表达式，公式为

(5)

(6)

参数加权矩阵We取决于m，所以必须使用迭代的计算方法.由于不同参数对观测的贡献是不同的，数据对不同参数的积分灵敏度也是变化的.为了消除这种差异，引入模型加权矩阵和数据加权矩阵，公式为

(7)

(8)

引入模型加权矩阵和数据加权矩阵后，公式 (2) 改写为

(9)

将目标函数 (9) 转换到加权密度参数域中求解mw，通过以下变量替换实现，公式为

(10)

(11)

(12)

为了简便，我们取目标函数的一半作为待求解的目标函数.新的目标函数为

(13)

采用共轭梯度算法求最优解mw，使得目标函数 (13) 最小.共轭梯度算法的实现过程，必须计算目标函数ϕ对加权密度mw的一阶导数.目标函数ϕ对密度m的一阶导数以及加权密度mw对密度m的一阶导数容易计算，公式为

(14)

(15)

利用公式 (14) 和公式 (15) 的耦合关系，能够得到目标函数ϕ对加权密度mw的一阶导数为

(16)

在处理位场数据时，传统的求导方法，不考虑We含有密度变量m的影响，视We为常量，因此将公式 (9) 视为伪二次函数，这样能够直接求得目标函数ϕ对mw的一阶导数为 (Portniaguine and Zhdanov, 2002):

(17)

对比新公式 (16) 和传统公式 (17) 的异同：新公式多出e－2We－2，对于某一密度为m的块体，此项为：

(18)

当某一块体密度m趋近于0时，公式 (18) 趋近于1，则对于此块体由公式 (17) 得到的梯度与公式 (16) 相同.当某一块体密度值m不为0时，本文推导出的传统公式 (17) 与公式 (16) 不同，本文考虑了We中含有变量m的影响，在求导过程中视其为变量，因此更合理.

共轭梯度算法求解目标函数最优解的过程中，聚焦因子e是常量，为了使求导公式简洁、计算过程简便，我们进一步修改目标函数 (13)，写成如下形式为

(19)

则目标函数对加权密度mw的一阶导数 (公式16) 改写为

(20)

将求导公式 (20) 带入计算过程，寻找使目标函数 (19) 最小的解mw，再利用公式 (12) 转换求出密度分布m.

2.2 密度上下限的选取

理论上，重力位在三个方向上的一阶导数是Vx、Vy、Vz，全张量重力梯度是重力位V的二阶导数，描述的是Vx、Vy、Vz在三个方向上的变化情况，与均匀的背景围岩密度无关.所以，在反演地下密度异常的过程中，我们将观测数据看成是由异常体引起的，与背景围岩无关，一般将背景密度视为0.如果异常体的剩余密度mc大于0，将反演范围设为0 < m < mc；如果异常体的剩余密度mc < 0，将反演密度设为mc < m < 0.而聚焦反演无论是在加权密度域求解还是直接在密度域求解，都涉及到We的计算和使用.其中，聚焦因子e的选取非常麻烦，理论上e为趋近于0的极小值，但是，在实际计算过程中，e太小会在某一块体密度趋于0时，We中出现奇异值，使结果不稳定，因此聚焦因子的合理设置对反演结果影响非常大.本文通过密度上下限的重新设定成功解决了这一问题.例如剩余密度mc>0的正演模型，我们将其设置成背景异常密度为1 g·cm-3 (非0)、异常体密度差为1+mc的模型，此时异常体剩余密度仍然是mc，反演时强制约束密度范围为1 < m < (1+mc).由于密度模型中密度最小值为1，不存在趋近于0的极小值，即使聚焦因子e为趋近于0的极小值，We也不会出现奇异值，能得到稳定的反演结果.这降低了方法对聚焦因子选取的依赖性，模型试验也证明e的选取非常自由，同时反演的稳定性和准确性有明显提高.对于实际测量数据，由于观测异常值是仅由异常体的剩余密度引起的，我们要在同一观测系统下, 假设背景异常密度均匀且为1 g·cm-3，模拟一组全张量重力梯度异常数据与实际测量数据相加，作为待反演数据.当然，最后要在反演结果中减去背景异常1 g·cm-3，才能得到异常体的剩余密度.

2.3 反演计算过程

首先，设k=0，m0=0，mw0=0，e=10-10计算We=diag ((m02+e2)－1/2)，Aw=Wd AWm－1 We－1, Q=We-2 AwT Aw +α0e2 I, q=We-2 AwT dw.共轭梯度算法的初始搜索方向选择目标函数的最速下降方向d0=－ f0=－ Qmw+ q，初始搜索步长为

if mk(i, 1)>up_bound (密度上限)

mk(i, 1)=up_bound

end

if mk(i, 1) < low_bound (密度下限)

mk(i, 1)=low_bound

end

更新We、mwk、Aw：We=diag ((mk2+ε)－1/2)

计算ϕ对mwk的一阶导数: fk=Qmwk－ q

k次沿mwk搜索方向: dk=－ fk+βdk－1(β=

k次沿mwk搜索步长:

end.

3 理论模型试验

将上述方法运用于合成的理论数据进行试验.地下三维空间剖分为23×16×10=3680个大小为100 m×100 m×100 m的规则立方体单元，背景剩余密度为0，两个密度差为1 g·cm-3的规则立方体大小均400 m×400 m×300 m，埋深均为200 m，相距500 m，如图 1所示.23×16=368个地面观测点位于每个密度单元的中心，地面观测网格为100 m×100 m，在理论观测数据中加入平均值为0，标准差为理论观测数据幅值5%的高斯噪声，如图 2所示.全张量重力梯度数据中每个单分量的幅值不同，所以噪声标准差也不同，如表 1所示.假设在同一地下3维空间的3680个立方体单元密度均为1 g·cm-3，地面观测无噪声数据如图 3，将两组数据对应相加组成的数据如图 4所示.图 4数据作为我们在反演过程中处理的数据.

图 1 密度差为1 g·cm-3的组合立方体理论模型图 Fig. 1 Synthetic model comprising two prisms with density contrast 1 g·cm-3
表 1 全张量重力梯度理论观测数据各分量或分量组合中高斯噪声标准差 Table 1 Standard deviations of Gaussian noise in synthetic gravity gradient data
图 2 理论模型正演全张量重力梯度地面观测数据 (含5%高斯噪声) (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 2 Surface gravity gradients calculated from synthetic model with 5% Gaussian noise
图 3 密度差为1 g·cm-3的背景密度引起的地面观测异常 (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 3 Surface anomalies from the background homogeneous host rocks with density of 1 g·cm-3
图 4 图 2和图 3中两组数据对应加和所得数据 (即反演过程中待处理数据) (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 4 Sum of two groups of data (Fig. 1 and 3) for inversion

将图 4中的数据带入目标函数 (19).利用共轭梯度算法求解重加权聚焦反演，利用新推导出的公式 (20) 求目标函数的梯度，正则化因子α的选择方法是利用Portniaguine和Zhdanov在2002年采用的方法，(k>1)，本文q=0.9.由正演模型可知，密度参数上限为2 g·am-3，下限为1 g·cm-3.聚焦因子e首先选择10-10，分别处理Txx、Tyy、Txy、Txz、Tyz、Txx| Txy| Tyy、Txz| Tyz| Tzz、5C(Txx| Txy| Tyy| Txz| Tyz)、6C(5C+ Tzz)，对于某组数据，当预测数据和观测数据达到拟合要求时迭代停止.给出3个水平分量Txx| Txy| Tyy联合反演结果如图 5c所示，反演结果收敛，密度分布聚焦，边界清晰，能够较好的恢复理论模型的位置.其他分量或分量组合的反演结果类似，不一一展示.反演模型与理论模型的拟合差如表 2中第3行所示.由6个分量单独反演得到的地下密度分布中，反演拟合差由大到小依次为Txy> Txz> Txx> Tyy> Tyz> Tzz，由单分量Tzz恢复的地下异常体模型与理论模型的拟合差最小，说明对于本实验数据，由Tzz分量恢复的地下密度分布最好.由4种组合分量反演得到的地下密度分布中，拟合差大小依次为6 C >3 C(Txz| Tyz| Tzz)>5 C>3 C(Txx| Txy| Tyy)，由3个水平分量的组合恢复的地下异常体与理论模型的拟合差最小，说明观测数据量的增加不一定会提高反演的准确性.另外，我们注意到，在单分量反演中，单独垂直分量反演效果大体上优于水平分量，但是3个水平分量的组合优于3个垂直分量的组合.

图 5 处理3 C(Txx| Txy| Tyy) 恢复的模型 (参数e取不同值) (a) e=1；(b) e=10-3；(c) e=10-10；(d) e=10-15. Fig. 5 Recovered models by 3 C (Txx| Tyy| Txy) for different e values
表 2 参数e取不同值处理全张量重力梯度的不同分量和分量组合得到的反演模型与理论模型的均方根误差 (单位：g·cm-3) Table 2 Root-mean-square errors (RMSE) between inversion models in Fig. 5 and theoretical model in Fig. 1 for four different e values

我们以e=10-10的计算为标准，对比不同聚焦因子对反演结果的影响.当e=1、10-3和10-15时，处理3个水平分量Txx| Txy|Tyy，反演过程中的反演次数与e=10-10保持相同，反演结果如图 5a、5b、5d所示.图 5b和图 5d中的密度分布与图 5c极为相似，异常体收敛聚焦；图 5a中的反演结果还未达到完全收敛.其他分量或分量组合的反演结果与其类似，不一一展示.反演模型与理论模型的拟合差如表 2第1、2和4行所示.对任意一组数据，拟合差由大到小依次为1>10-3=10-10=10-15，说明当聚焦因子e=1时，收敛速度慢，密度分布还未达到完全收敛.当聚焦因子取10-3、10-10和10-15时，相同的反演次数得到的结果相同，说明收敛速度相等.这也说明了，聚焦因子选择对反演结果影响小，任意大于0的极小值，都能够得到稳定收敛的结果.克服了以往聚焦反演，聚焦因子的选择对反演结果的影响.

4 讨论

为了充分验证方法的实用性和可靠性，我们在相同的观测系统下设置模型试验，分析讨论背景异常密度不同对反演结果的影响，以及异常体密度差不同对反演结果的影响.首先，选择异常体密度差为1 g·cm-3，背景异常密度为1 g·cm-3作为对照组研究背景异常密度不同对反演结果的影响，观测数据选择处理3个水平分量的组合 (Txx| Txy| Tyy).基于实际地质情况分别设置背景异常密度为0、0.1 g·cm-3、0.5 g·cm-3、5 g·cm-3、10 g·cm-3(实际地质情况下，背景异常密度一般不会大于10 g·cm-3).正则化参数α的选择同上，聚焦因子e=10－10，在相同的计算条件下迭代次数相同，对比反演结果 (图 6).反演模型和理论模型的均方根误差如表 3.当背景异常密度为0时，由于e太小反演结果不收敛且无收敛趋势，即使继续增加迭代次数反演结果也无法收敛；当背景异常密度为0.1 g·cm-3、5 g·cm-3、10 g·cm-3时，反演结果未达到完全收敛，单具备收敛趋势，这说明收敛速度较慢，继续增加迭代次数反演结果能够完全收敛；当背景异常密度为0.5 g·cm-3、1 g·cm-3时，反演结果基本达到收敛；背景异常密度为0.5 g·cm-3时，收敛速度最快.我们得到两点结论：聚焦因子e取很小值时，常规的背景异常密度为0，反演结果不收敛，加适当的背景异常可使结果收敛；为了节约计算时间，提高效率，建议选择背景密度异常为0.5 g·cm-3.

图 6 相同计算条件下，设置不同背景异常密度反演Txx| Txy| Tyy数据得到密度分布模型，背景异常密度分别为 (a) 0；(b) 0.1 g·cm-3；(c) 0.5 g·cm-3；(d) 1 g·cm-3；(e) 5 g·cm-3；(f) 10 g·cm-3(异常体密度差为1 g·cm-3，e=10-10) Fig. 6 The recovered models by inverting Txx| Txy| Tyy for the different density of background host rocks
表 3 图 6中反演模型与真实模型的均方根误差 Table 3 Root-mean-square errors (RMSE) between inversion models in Fig. 6 and actual model for different densities of background host rocks (density contrast is 1 g·cm-3, e=10-10)

最后，选择异常体密度差为1 g·cm-3，背景异常密度为1 g·cm-3作为对照组研究异常体密度差不同对反演结果的影响，观测数据选择处理3个水平分量的组合 (Txx| Txy| Tyy).基于实际地质情况分别选择异常体密度差为0.2 g·cm-3、0.5 g·cm-3、0.8 g·cm-3、1.2 g·cm-3、1.5 g·cm-3，在相同计算条件下，对比反演结果 (图 7).反演模型和理论模型的均方根误差如表 4.结果表明，异常体密度差在合理范围内，反演结果均能达到收敛；异常体密度差越大，达到收敛所需的迭代次数越少，证明收敛速度越快.

图 7 设置不同异常体密度差反演Txx| Txy| Tyy数据得到的密度分布模型，异常体密度差分别为 (a) 0.2 g·cm-3；(b) 0.5 g·cm-3；(c) 0.8 g·cm-3；(d) 1 g·cm-3；(e) 1.2 g·cm-3；(f) 1.5 g·cm-3(背景异常密度差为1 g·cm-3，e=10-10) Fig. 7 Recovered models by inverting Txx| Txy| Tyy for different density contrasts. (a) 0.2 g·cm-3; (b) 0.5 g·cm-3; (c) 0.8 g·cm-3; (d) 1 g·cm-3; (e) 1.2 g·cm-3; (f) 1.5 g·cm-3 (background density is 1 g·cm-3, e=10-10)
表 4 图 7中反演模型与真实模型的均方根误差 Table 4 Root-mean-square errors (RMSE) (in g·cm-3) between inversion models in Fig. 6 and actual model in Fig. 1 for four different density contrasts (background density is 1 g·cm-3, e=10-10)
5 实际数据应用

一个多世纪以来，美国墨西哥湾沿岸地区一直是石油和天然气的大量产出地，沿岸大多数的储层都是在盐丘侧翼、背斜、地层尖灭等处被找到的 (Branson, 1991).文顿盐丘位于美国路易斯安纳州东南部墨西哥湾沿岸，作为一个典型的墨西哥湾海岸盐丘，这是第一个在侧翼找到油气的盐丘油田 (Ennen, 2012).盐丘坐落在几乎完整的新生代和中生代沉积层顶部，侵入第三纪中新世弗莱明地层、渐新世的维克斯堡地层及始新世的杰克逊地层，大部分油气来源于第三纪中新世弗莱明地层和渐新世的维克斯堡地层 (Branson, 1991)，且多出现在盐丘东部、东北部、北部，有学者认为源岩是始新世的杰克逊地层 (Swanson and Karlsen, 2009).

文顿盐丘包含一个大规模的盐核和定义非常好的岩盖延伸覆盖在顶部 (Thompson et al., 1928; Coker et al., 2007).岩盖在盐核顶部，埋深相对较浅，故本地区有很多钻井穿过岩盖，岩盖成分、展布范围、埋深、厚度、密度等信息都相对完善.岩盖主要由石灰岩、石膏、硬石膏组成；岩盖展布由东向西为1500 m，由南到北为1100 m；岩盖顶部埋深约为160 m，底部埋深约为360 m，密度约为2.75 g·cm-3，周围沉积物以及盐核密度都约为2.2 g·cm-3(Ennen, 2012)，所以，岩盖的密度差为0.55 g·cm-3，能够引起重力及重力梯度数据观测异常.2008年，Bell Geospace公司在这个地区搜集了航空全张量重力梯度数据，很多学者对这组数据进行了分析和处理 (Oliveira and Barbosa, 2013; Geng et al., 2014; Hou et al., 2016; Qin et al., 2016).为了消除地形影响，得到由岩盖引起的观测数据异常，我们用2.2 g·cm-3的密度对原始数据做了地形改正.设x、y、z的正方向分别代表指东方向、北方向、垂直向下的方向，截取4 km×4 km的子区域作为我们的研究的观测区域 (图 8)，观测网格为100 m×100 m，共有41×41个观测点.地下区域划分成nx×ny×nz=41×41×20个规则立方体，大小为100 m×100 m×50 m，所以，反演深度为1000 m.

图 8 美国文顿盐丘地区实际观测全张量重力梯度数据 (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 8 Measured gravity gradients over the Vinton dome, USA

同样，为了消除聚焦因子选择对反演结果的影响，我们假设背景异常密度为1 g·cm-3，正演计算得到一个无噪声的理论观测异常数据 (图 9)，将其与实测观测异常 (图 8) 相加，这组理论观测异常与实测观测异常的组合数据 (图 10) 作为反演过程中的观测数据.反演得到的地下三维密度体，背景异常应为1 g·cm-3，异常体的剩余密度应为1.55 g·cm-3.所以密度约束上下限为1 g·cm-3-1.55 g·cm-3.最后，在反演结果中再减去背景密度1 g·cm-3，得到异常体的剩余密度.

图 9 背景异常密度差为1 g·cm-3得到的理论数据 (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 9 Synthetic data in the same observation system with Fig. 8
图 10 文顿盐丘地区实际数据与图 9中理论数据对应加和，用做反演中的观测数据 (a) Txx；(b) Txy；(c) Txz；(d) Tyy；(e) Tyz；(f) Tzz. Fig. 10 Measured data over the Vinton dome added with synthetic data in Fig. 8 used in inversion

对比图 10与图 8，实际数据中加入合成数据后，噪声水平降低，将带有合成数据的实际数据带入目标函数对反演可能是有益的.由理论模型部分，我们已经知道单一的Tzz分量以及3个水平分量的组合 (Txx| Txy| Tyy) 能够得到令人满意的效果，所以，在实际数据部分，我们分别利用Tzz数据以及3个水平分量组合反演地下密度分布.反演过程中正则化因子α的选择与理论模型部分相同.聚焦因子e取10-10.最后在密度分布中减去人为加入的背景异常得到异常体剩余密度.处理Tzz数据得到的结果如图 11所示，处理3个水平分量组合得到的结果如图 12所示，两个图很相似，初步说明结果可信.

图 11 处理文顿盐丘地区实际Tzz数据恢复的密度分布 (a) 水平切片z=0.35 km；(b) 南北方向切片y=2.1 km；(c) 东西方向切片x=2.1 km；(d) 去除密度小于0.2 g·cm-3的块体的三维立体图俯视图；(e) 去除密度小于0.2 g·cm-3的块体的三维立体图仰视图. Fig. 11 Recovered density distribution by inverting Tzz data over Vinton salt dome (a) Horizontal slice of z=0.35 km; (b) NS-trending slice of y=2.1 km; (c) EW-trending slice of x=2.1 km; (d) Top view of 3D stereogram with density contrast less than 0.2 g·cm-3 removed; (e) Bottom view of 3D stereogram with density contrast less than 0.2 g·cm-3 removed.
图 12 处理文顿盐丘地区实际Txx| Tyy| Txy数据恢复的密度分布 (a) 水平切片z=0.35 km；(b) 南北方向切片y=2.1 km；(c) 东西方向切片x=2.1 km；(d) 去除密度小于0.2 g·cm-3的块体的三维立体图俯视图；(e) 去除密度小于0.2 g·cm-3的块体的三维立体图仰视图. Fig. 12 Recovered density distribution by inverting Txx| Tyy| Txy data over Vinton salt dome (a) Horizontal slice of z=0.35 km; (b) NS slice of y=2.1 km; (c) EW slice of x=2.1 km; (d) Top view of 3D stereogram with density contrast less than 0.2 g·cm-3 removed; (e) Bottom view of 3D stereogram with density contrast less than 0.2 g.cm-3 removed.

我们的反演结果 (图 11、图 12) 显示，异常密度体集中分布，边界清晰，大体图像为一个西北部有一缺口的楔状体.南部及东南部埋深最浅，埋深200 m，中间埋深250 m，西北部埋深较深，为300 m；东南部异常体较厚，平均厚度为200 m，西北部异常体比较薄，平均厚度为100 m.南部边界变化陡，北部变化平缓.东西方向延伸1300 m, 南北方向延伸1100 m.将本文结果与前人研究成果对比，如表 5.基于地质资料建模得到的岩盖埋深较浅，厚度较小；而基于密度反演理论得到的岩盖模型埋深较深，且厚度较大.本文基于密度反演理论的方法能够得到更接近地质资料的结果，体现了方法处理实际数据的优越性.

表 5 前人方法与本文方法估计或计算的文顿盐丘埋深及厚度结果比较 Table 5 Comparison of results by this paper and previous work
6 结论

本文基于正则化反演理论，利用聚焦反演方法处理全张量重力梯度数据以恢复地下三维密度分布.应用共轭梯度算法在加权密度域求解目标函数最优化问题，在考虑变加权函数W e中含有密度变量的情况下，重新推导了目标函数的梯度，得到了与传统梯度不同的新的公式，新公式更合理.另外，在正演过程中，通过在研究的反演空间强制加背景异常，并将其产生的观测异常与异常体观测异常融合，共同进入反演程序，达到改变密度参数上下限的目的，从而使聚焦因子的选择更加灵活.最后，我们讨论了所加背景异常密度的大小对反演结果的影响，以及异常体密度差的大小对反演结果的影响，结果证明，背景异常为0.5 g·cm-3，收敛速度最快；异常体密度差在合理范围内，模型都能够被合理恢复，说明加背景异常是可行的.本文方法处理文顿盐丘地区实际重力梯度数据，反演计算出岩盖大体位于地下200 m到450 m，与前人研究结果相符，证明本文方法具有处理实际数据的能力.

致谢

感谢Bell Geospace允许使用其在文顿盐丘地区获得的航空全张量重力梯度数据.

参考文献
Barbosa V C F, Silva J B C. 1994. Generalized compact gravity inversion. Geophysics, 59(1): 57-68. DOI:10.1190/1.1443534
Branson R B. 1991. Productive trends and production history, south Louisiana and adjacent offshore. New Orleans Geological Society, 61-70.
Coker M O, Bhattacharya J P, Marfurt K J. 2007. Fracture patterns within mudstones on the flanks of a salt dome:syneresis or slumping?. Gulf Coast Association of Geological Societies Transactions, 57: 125-137.
Ennen C. 2012. Mapping gas-charged fault blocks around the Vinton Salt Dome, Louisiana using gravity gradiometry data[Master's thesis]. Houston:University of Houston.
Geng M X, Huang D N, Yang Q J, et al. 2014. 3D inversion of airborne gravity-gradiometry data using cokriging. Geophysics, 79(4): G37-G47. DOI:10.1190/geo2013-0393.1
Geng M X, Huang D N, Yu P, et al. 2016. Three-dimensional constrained inversion of full tensor gradiometer data based on cokriging method. Chinese J. Geophys., 59(5): 1849-1860. DOI:10.6038/cjg20160528
Guillen A, Menichetti V. 1984. Gravity and magnetic inversion with minimization of a specific functional. Geophysics, 49(8): 1354-1360. DOI:10.1190/1.1441761
Hou Z L, Huang D N, Wei X H. 2016. Fast inversion of probability tomography with gravity gradiometry data based on hybrid parallel programming. Journal of Applied Geophysics, 124: 27-38. DOI:10.1016/j.jappgeo.2015.11.009
Last B J, Kubik K. 1983. Compact gravity inversion. Geophysics, 48(6): 713-721. DOI:10.1190/1.1441501
Oliveira V C Jr, Barbosa V C F. 2013. 3-D radial gravity gradient inversion. Geophysical Journal International, 195(2): 883-902. DOI:10.1093/gji/ggt307
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596
Portniaguine O, Zhdanov M S. 2002. 3-D magnetic inversion with data compression and image focusing. Geophysics, 67(5): 1532-1541. DOI:10.1190/1.1512749
Qin P B, Huang D N. 2016. Integrated gravity and gravity gradient data focusing inversion. Chinese J. Geophys., 59(6): 2203-2224. DOI:10.6038/cjg20160624
Qin P B, Huang D N, Yuan Y, et al. 2016. Integrated gravity and gravity gradient 3D inversion using the non-linear conjugate gradient. Journal of Applied Geophysics, 126: 52-73. DOI:10.1016/j.jappgeo.2016.01.013
Silva Dias F J S, Barbosa V C F, Silva J B C. 2009. 3D gravity inversion through an adaptive-learning procedure. Geophysics, 74(3): I9-I21. DOI:10.1190/1.3092775
Silva Dias F J S, Barbosa V C F, Silva J B C. 2011. Adaptive learning 3D gravity inversion for salt-body imaging. Geophysics, 76(3): I49-I57. DOI:10.1190/1.3555078
Silva J B C, Barbosa V C F. 2006. Interactive gravity inversion. Geophysics, 71(1): J1-J9. DOI:10.1190/1.2168010
Swanson S M, Karlsen A W. 2009. USGS assessment of undiscovered oil and gas resources for the Oligocene Frio and Anahuac formations, onshore Gulf of Mexico basin, USA. Search and Discovery.
Ye Z R, Liu L T, Liang X H, et al. 2016. Formula of Moho inversion in the spectral domain using vertical gravity gradient and its application. Chinese J. Geophys., 59(2): 476-483. DOI:10.6038/cjg20160207
Zhdanov M S. 2009. New advances in regularized inversion of gravity and electromagnetic data. Geophysical Prospecting, 57(4): 463-478. DOI:10.1111/gpr.2009.57.issue-4
Zhdanov M S, Ellis R, Mukherjee S. 2004. Three-dimensional regularized focusing inversion of gravity gradient tensor component data. Geophysics, 69(4): 925-937. DOI:10.1190/1.1778236
耿美霞, 黄大年, 于平, 等. 2016. 基于协同克里金方法的重力梯度全张量三维约束反演. 地球物理学报, 59(5): 1849–1860. DOI:10.6038/cjg20160528
秦朋波, 黄大年. 2016. 重力和重力梯度数据联合聚焦反演方法. 地球物理学报, 59(6): 2203–2224. DOI:10.6038/cjg20160624
叶周润, 柳林涛, 梁星辉, 等. 2016. 垂直重力梯度反演Moho面的频谱域公式及其应用. 地球物理学报, 59(2): 476–483. DOI:10.6038/cjg20160207