重力梯度张量的预条件共轭梯度法反演

陈少华，朱自强，鲁光银，曹书锦

(中南大学 地球科学与信息物理学院，湖南 长沙，410083)

摘要：针对重力梯度张量反演中的问题，提出基于预条件共轭梯度法的重力梯度张量反演。通过在目标函数中加入粗糙度对模型进行约束以避免反演参数远多于采集点数的欠定问题不稳定，并在目标函数中添加深度加权矩阵对核函数进行补偿，以避免核函数随着深度的增大而快速衰减的问题。分别反演、比较各重力梯度张量分量和联合5个独立分量，并将重力梯度张量5个独立联合反演应用于Y型岩脉。研究结果表明：联合反演效果明显优于单一分量的反演效果，且能较好地与原始模型相吻合，证明了本文算法的有效性。

关键字：重力梯度张量；深度加权函数；预条件共轭梯度法；反演

中图分类号：P631          文献标志码：A         文章编号：1672-7207(2013)02-0619-07

Inversion of gravity gradient tensor based on preconditioned conjugate gradient

CHEN Shaohua, ZHU Ziqiang, LU Guangyin, CAO Shujin

(School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Abstract: The 3d inversion based on preconditioned conjugate gradient was proposed for Tensor Gradient gravity data. Roughness matrix was added in the objective function to avoid the unstable defect of inversion problem caused by the number of inverse parameter far than observational point’s. Depth weighting function was integrated into the objective function to compensate kernel function avoid diminishing rapidly with depth. The inverse tests with every single component and joint five independent measured components of Tensor Gradient gravity data of the synthesized simple model and Y-type dyke show that the joint inversion is better than the others. The inversion result of the last test with joint five independent measured components of Tensor Gradient gravity data of the synthesized Y-type dyke coincides well with the real model. The results demonstrate the effect of the new algorithm.

Key words: gravity gradient tensor; depth weighting function; preconditioned conjugate gradient; inversion

重力勘探方法在区域地质调查、成矿远景区划分评价、水文及工程地质勘察等方面发挥重要作用。常规重力勘探只观测重力位的铅垂一阶导数，而重力梯度勘探则观测重力位的二阶导数，因此，相比传统的重力测量，重力梯度测量有更高的分辨率^[1]，更能够清晰地反映场源体边界形态，且其张量数据较为丰富，能够提高异常体的解释精度。作为一种前沿性的重力测量方法，国外学者不仅对重力梯度测量仪器进行了大量开发，还对采集的数据进行了大量研究。Zhang等^[2]利用欧拉反褶积法解释了重力全张量数据； Zhdanov等^[3]利用重加权正则化共轭梯度法，采用聚焦正则化对目标函数进行约束，数值模拟表明该算法对尖锐边界有较强的反演效果；Wilson等^[4]针对巨大的航空重力梯度测量数据，采用并行算法，利用共轭梯度法对其进行反演，该算法在Vredefort得到了具体应用。此外，国外对重力张量数据的去噪、谱分析等方面的研究也取得了一定的成果^[5-7]。在国内，姚长利等^[8]针对三维物性反演中存在的困难和问题，提出三维物性反演的随机子域方法技术；郭良辉等^[9]提出了基于重力梯度数据异常分离的三维相关成像方法提高成像分辨率,数值模拟结果表面该方法能够使重力梯度数据具有良好的横向及纵向分辨率。与国外相比，国内对这一方面的研究较少，目前还没有开发出重力梯度张量测量仪器，对重力梯度张量数据的解释也较少。相较于非线性反演方法，线性反演算法出现的时间较长，应用较为广泛。传统线性反演算法中雅可比矩阵的存储需占用大量内存，计算该矩阵与其转置相乘将耗费大量时间且占用大量内存，而共轭梯度法只要计算两次雅可比矩阵与向量的乘积^[10-11]，计算时间将大量减少且所需内存小，有利于大规模计算。针对雅可比矩阵的条件数较大而产生病态解的问题，首先对方程组进行预条件处理^[12-13]，通过减小该矩阵的条件数来增强线性方程组的稳定性。本文作者在目标函数中加入粗糙度函数及深度加权函数，利用预条件共轭梯度法对目标函数进行反演，数值模拟结果证明了反演的有效性；针对重力梯度张量采集的数据有5个分量是独立的，充分利用重力梯度张量数据的丰富性，利用5个独立分量进行反演，取得了较好的效果。

1  正演

首先，在笛卡儿坐标系下(图1)，设地面上某点O为坐标原点。根据万有引力定律，一个剩余密度为ρ、体积为V的地质体，其在外部空间任意一定产生的引力位为：

              (1)

式中：r为观测点P(x,y,z)到异常点Q(ξ,η,ζ)的矢量；G为万有引力常量。

    (2)

式中：g，g_x，g_y和g_z均为笛卡儿坐标系(x,y,z)的函数；g为U(r)的一次导数；U(r)的二阶导数可表示为：

     (3)

由于在地球外部，位场满足拉普拉斯方程G_xx+G_yy+G_zz=0，以及重力场的无旋性：G_xy=G_yx，G_xz=G_zx，G_yz=G_zy，因此，在重力梯度张量中只有5个分量是独立的。

图1  笛卡儿坐标系下的长方体模型及地下空间剖分

Fig.1  Rectangular model and subsurface division in Cartesian coordinates

将探测区域的地下空间离散成一系列大小相等的长方体^[14]，赋予每个小长方体1个固定的剩余密度，在探测区域建立笛卡儿坐标系(图1)，x轴指向正北方向，y轴指向东方向，z轴垂直向下，则地面任意观测点各梯度分量的异常值可表示为：

      (4)

式中：代表某重力梯度分量中第i个长方体的几何架构，称为位置函数；ρ_i为第i个长方体的剩余密度；M为所划分长方体的个数。

则由地下空间离散的所有长方体在地面所有观测点产生的重力梯度异常值满足线性关系：

     (5)

式中：是N×1维观测数据向量，是M×1维长方体单元的剩余密度；G_αβ的元素G_ij表示第j个长方体单元在地面上第i个观测点产生的响应。

2  反演

2.1  粗糙度约束

重力梯度张量的反演便是求解线性方程组(5)。由于采集的数据向量d远少于要反演的未知数m，因此该反演是个欠定问题，方程组的解不唯一且不稳定，必须对目标函数加上适当的约束来缩小解的范围。本文中通过引入粗糙度函数^[15]，来对模型空间进行约束，目标函数可表示为：

          (6)

式中：，为数据拟合泛函；为模型稳定泛函(也称稳定器)；λ为拉格朗日乘子。

定义模型粗糙度R为m(r)在x，y和z方向一阶偏微分的平方和，即

          (7)

2.2  深度加权

与其他位场数据一样,重力梯度数据没有深度分辨率。从式(1)和(3)可知：地面重力梯度值与场源到观测点的立方成反比，即其核函数随着距离的增大快速衰减，而构造极小化模型目标函数时，重构模型的核矩阵是线性的，这样导致模型重构时靠近地表的网格单元比深部的网格单元有更大的权重，反演的构造多集中于地表，不能反映地质体的真实深度情况，因此，必须引进1个合适的深度加权因子来抵消核函数随着深度增大的衰减量，Li等^[16]在3D磁化率成像中引入深度加权：

              (8)

式中：z为块体单元中心点埋深；z₀和β为常数。β一般情况下取2~4；z₀取决于块体单元的尺寸以及观测面的高度。

将式(8)及最小化模型构造代入式(7)：

         (9)

式中：α_s，α_x，α_y和α_z为模型目标函数中各项的加权因子。

对式(9)按照模型区域的网格剖分进行离散化，用有限差分代替微分算子，写成如下矩阵形式：

               (10)

每一个矩阵分量都可以写成2个独立的矩阵与1个系数矩阵的乘积，即

(i= s，x，y，z)         (11)

式中：R_i为有限差分算子；Z为离散化的深度加权对角线矩阵。将式(10)代入式(6)：

   (12)

式(12)可写成方程组：

            (13)

则可用代表雅可比矩阵；，则式(13)可写成：

Am=b                 (14)

当对采集的5个独立分量数据进行联合反演时，可用代表雅可比矩阵，为5个独立分量采集数据的组合，将A及b代入式(14)即可得到反演的方程组。

2.3  预条件共轭梯度法

共轭梯度法是一种通过逐步迭代求解线性化或非线性化方程的算法。该方法求解的目标函数为二次方程，目标函数的解是m=A^-1b。这与式(14)的解是等价的，因此，可用共轭梯度法求解式(14)。与仅考虑梯度方向的最速下降法相比，共轭梯度法结合了梯度方向和共轭方向，利用初始点的梯度构造一组共轭方向，沿着这组方向而不是负梯度方向去搜索目标函数极小点，根据共轭方向的性质，共轭梯度法具有二次终止性^[17]。

由于模型参数的规模较大，在实际计算中，收敛的次数越少越好，而方程的收敛速度由雅可比矩阵的条件数决定的。位场中雅可比矩阵的条件数通常巨大，同时，van Decar等^[18]指出在目标函数中增加粗糙度约束将减缓收敛速度，为了避免过多的收敛次数及增强解的稳定性，将常规的方程A^tAm=A^tb改进成：

SA^tAm=SA^tb                (15)

因此，预条件共轭梯度法的计算步骤：

(1) 给定初始值ε，令m₀=0；k=0；r₀=A^tb；

(2) 构造搜索方向，令p_k=z_k-1+β_kp_k-1；其中，当k=1时，β=0，p_k=z_k-1；当k＞1时，；

(3) m_k=m_k-1+α_kp_k，q_k=Ap_k，求步长 ；

(4) 计算r_k=r_k-1-α_kA^tq_k，若||r_k||₂＜ε，则停止计算，否则返回第(2)步。

其中：S为预条件因子。若S≈(A^tA)^-1，则SA^tA≈I(其中，I为单位矩阵)，那么预条件因子S就可以大大改善式(14)的条件数，使得特征值集中分布在对角线附近，有效地提高迭代速度，这便是预条件共轭梯度法。但是，这样求解S是非常困难的，在实际应用中，用简单的深度加权函数即式(8)作为预条件因子^[12]，能很好地改善式(14)的条件数，实际模型反演结果证明了这一点。

3  数值模拟

3.1  单分量反演

将场源空间划分为20×20×10=4 000个单元格，每个单元格沿x，y和z轴方向的长度均为50 m。如图2所示，地下空间存在2个异常体，2个异常体的长×宽×高均为200 m×200 m×200 m，异常体顶板埋深均为125 m，底板为325 m，异常体的剩余密度为1.0 t/m³。总共有20×20=400个采集数据，数据采集点位于表层，沿x和y轴间距都为50 m。

图2  正演模型和G_xx，G_xy，G_xz，G_yy，G_yz和G_zz的反演结果

Fig.2 Forward model and inversion results of G_xx, G_xy, G_xz, G_yy, G_yz and G_zz

定义绝对拟合差：

反演参数的设置中，取α_x=α_y=α_z=1，α_x=0.00 005，z₀=20，β=4，收敛终止的条件是上一次迭代的绝对拟合差大于后一次的绝对拟合差。从以上各个单分量反演结果可知：每个分量基本能够反演出异常地质体的位置和形态，反演出来的剩余密度也都不同程度地接近真实值，但都出现了大小不等的负值，其中以G_xx和G_xz出现的负值较大，而G_xy，G_yy，G_zz和G_yz反演出来的负值相对较小；同时可知G_xx，G_xz，G_xy和G_yy分量的反演结果较为发散，在深度方向上的延伸较大，不能完全反演出异常点的位置；而G_yz及G_zz分量的反演结果集中于异常体周围，且准确地反映了目标体的深度。

为了各个分量的比较，定义相对拟合差R_ms：

图3所示为各个分量在不同迭代次数下的相对拟合差。从图3可知：G_xz分量的反演非常不稳定，随着迭代的进行R_ms上下波动比较厉害，G_xx分量也出现了小范围的上下波动，这可能是以上2个分量反演结果中出现较大负值的原因；G_xy及G_yy分量基本上随着迭代的进行R_ms逐渐变小，但减小的速度较慢，最终其R_ms仍然较大；G_yz及G_zz分量不仅随着迭代的进行R_ms快速减小，而且最终的R_ms也比较小。

图3  各个分量在不同迭代次数下的相对拟合差

Fig.3  Relative misfit of all single component in different numbers of iteration

3.2  全张量反演

虽然理论上可以证明同一方法不同分量并不能改变由位场性质决定的多解性^[19]，这由于引起位场多解性的一个重要原因是数据的不完整和不精确，而观测得到的数据总是受限于一定网格距和一定的精度。尽管如此，同一高度不同分量的观测数据在完整性上仍有所补充，并减少了误差引起的多解性。重力梯度张量相比于重力勘探，不仅有更高的分辨率，还在于其更为丰富的数据量，充分利用5个独立分量的信息将有利于提高反演的精度。

本文采用5个独立分量G_xx，G_xy，G_xz，G_yy及G_zz进行联合反演，并对反演数据加入10%的噪声干扰，检验反演的抗噪性能。取α_x=α_y=α_z=1，α_s=0.000 05，z₀=10，β=4，虽然随着迭代的进行，数据的绝对拟合差能不断地降低，但可能导致过拟合而引起多余构造，故本次收敛中取最多的迭代次数为30次，观测数据的正演模型见图2(a)。

与没有加入高斯噪声的图2中各个单分量相比较：图4(a)所示的反演的异常体不仅集中，而且和异常体的真实深度吻合，验证了同一种方法的不同分量的联合反演可以减少多解性，提高解的精度。图4(b)中加入10%的高斯噪声后反演的异常形态仍较为明显，只是在异常体的周围出现一些小异常，但可以忽略不计，表面5个分量联合反演有较好的抗噪性。

图4  无噪声干扰和加入10%噪声干扰的5个独立分量的联合反演

Fig.4  Joint-inversion of five independent components and without noise and with 10% Gaussian noise

3.3  合成模型反演

现设计1个Y型岩脉，该模型曾被姚长利等^{[8-9, 20]}应用于位场相关领域的分析。图5(a)和(b)所示分别为深度120 m的深度切片和图5(a)中A-A′的断面切片，图5(c)所示为120 m深度剖面的反演结果，图5(d)所示为图5(a)中A-A′断面切片的反演结果，白色方框为真实异常体的位置。异常区Ⅰ(颜色较浅)的异常体其剩余密度为0.8 t/m³ ，深度方向延伸不大；异常区Ⅱ(颜色较深)的异常体其剩余密度为1.0 t/m³，深度方向延伸大，两岩脉的延伸方向相反。将地下空间剖分为20×20×10=4 000个长方体，每个长方体的大小都为50 m×50 m×50 m。正演计算的观测点位于地表处，共有20×20=400个观测点，在观测数据中加入3%的高斯噪声，采用5个独立分量对地下异常体进行联合反演。

反演参数的设置中，取α_x=α_y=α_z=1，α_x=0.000 05，z₀=20，β=4，最多迭代次数为30次，加入先验参数范围约束，当解小于0时令其为0，当解大于1.1时令其为1.1。整个反演过程的迭代次数为30次，消耗的时间为21.742 587 s，相对拟合差为0.026 371。

图5(c)所示较好地反映了异常体的形态和位置，左边岩脉的形态和位置基本吻合，右边岩脉沿y轴方向完全吻合，在x轴方向稍微扩大了一点，但基本上与异常体的真实形态一致；图5(d)中基本重构出了模型的异常形态，其异常走向与Y型岩脉一致；结合图5(c)和5(d)可发现：反演算法在浅部取得了较清晰的结果，而在深部反演的分辨率则有所下降，表明虽然基于深度加权的重力梯度张量反演能够使物性的参数值大致能够分布在相应的位置，但重构模型的分辨率不可避免地随着深度的增加而降低。以上结果表明：基于粗糙度约束和深度加权的重力梯度张量反演，能够在一定程度上较好地重构模型的异常形态。

图5  合成模型切片图及异常体的反演结果切片图

Fig.5  Slice maps of the synthesized model and inversion results of anomalous mass

4  结论

(1) 提出了基于预条件共轭梯度法的重力梯度张量反演，收敛速度快，耗时较短，可用于数据量较大的反演计算。

(2) 相对于单分量反演，充分利用5个独立分量的信息进行联合反演，能够较准确地反演出异常体的形态及剩余密度特征。

(3) 仅用理想模型数据对本文的算法进行了验证，需要进一步采用现场采集的数据进行验证。

参考文献：

[1] 曾华霖. 重力场与重力勘探[M]. 北京: 地质出版社, 2005: 68-71.

ZENG Hualin. gravitational field and gravitational prospecting[M]. Beijing: Geological Publishing House, 2005: 68-71.

[2] Zhang C, Mushayandebvu M F, Reid A, et al. Euler deconvolution of gravity tensor gradient data[J]. Geophysics, 2000, 65(2): 512-520.

[3] Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925-937.

[4] Wilson G, Cuma M, Zhdanov M S. Massively parallel 3D inversion of gravity and gravity gradiometry data[J]. Preview, 2011, 152(6): 29-34.

[5] Lyrio L S, Tenorio L, Li Y G. Efficient automatic denoising of gravity gradiometry data[J]. Geophysics, 2004, 69(3): 772-782.

[6] While J, Jackson A, Smit D, et al. Spectral analysis of gravity gradiometry profiles[J]. Geophysics, 2006, 71(1): J11-J22.

[7] Afif H. Understanding gravity gradients: Atutorial[J]. The Leading Edge, 2006, 25(8): 942-947.

[8] 姚长利, 郑满元, 张聿文. 重磁异常三维物性反演随机子域法方法技术[J]. 地球物理学报, 2007, 50(5): 1576-1583.

YAO Changli, ZHENG Manyuan, ZHANG Yuwen. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces[J]. Chinese Journal of Geophysics, 2007, 50(5): 1576-1583.

[9] 郭良辉, 孟晓红, 石磊, 等. 重力和重力梯度数据三维相关成像[J]. 地球物理学报, 2009, 52(4): 1098-1106.

GUO Lianghui, MENG Xiaohong, SHI Lei, et al. 3-D correlation imaging for gravity and gravity gradiometry data[J]. Chinese Journal of Geophysics, 2009, 52(4): 1098-1106.

[10] 吴小平, 汪彤彤. 利用共轭梯度方法的电阻率三维反演的若干问题研究[J]. 地震地质, 2001, 23(2): 321-327.

WU Xiaoping, WANG Tongtong. Study on some problems for 3D resistivity inversion using conjugate gradient[J]. Seismology and Geology, 2001, 23(2): 321-327.

[11] 吴小平, 徐国明. 利用共轭梯度法的电阻率三维反演研究[J]. 地球物理学报, 2000, 43(3): 420-427.

WU Xiaoping, XU Guoming. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420-427.

[12] Pilkington M. 3-D magnetic imaging using conjugate gradient[J]. Geophysics, 1997, 62(4): 1132-1142.

[13] Pilkington M. 3D magnetic data-space inversion with sparseness constraints[J]. Geophysics, 2009, 74(1): L7-L15.

[14] 骆遥, 姚长利, 薛典军, 等. 2.5D地质体重磁异常无解析奇点正演计算研究[J]. 石油地球物理勘探, 2009, 44(4): 487-493.

LUO Yao, YAO Changli, XUE Dianjun, et al. Analytic singularity-free forward computation of gravity and magnetic anomaly for 2.5-D geologic body[J]. Oil Geophsical Prospecting, 2009, 44(4): 487-493.

[15] Constable S C, Parker R L, Constable C G. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 1987, 52(3): 289-300.

[16] Li Y G, Oldenburg D W. 3-D inversion of magnetic data[J]. Geophysics, 1996, 61(2): 394-408.

[17] 姚姚. 地球物理反演基本理论与应用方法[M]. 武汉: 中国地质大学出版社, 2002: 57-59.

YAO Yao. Basic theory and application of geophysical inversion methods[M]. Wuhan: China University of Geosciences Press, 2002: 57-59.

[18] van Decar J, Snieder R. Obtaining smooth solution to large, linear, inverse problems[J]. Geophysics, 1994, 59(5): 818-829.

[19] 管志宁. 地磁场与磁力勘探[M]. 北京: 地质出版社, 2005: 199-204.

GUAN Zhining. Magnetic field and magnetic exploration[M]. Beijing: Geological Publishing House, 2005: 199-204.

[20] Li Y G, Oldenburg D W. 3-D inversion of gravity data[J]. Geophysics, 1998, 63: 103-119.

(编辑  杨幼平)

收稿日期：2012-03-12；修回日期：2012-05-01

基金项目：国家自然科学基金资助项目(41174061)；中南大学自由探索计划项目(2011QNZT011)

通信作者：鲁光银(1976-)，男，湖北宜昌人，教授，从事地震及重磁相关领域的研究；电话：13975894898；E-mail：csulgy@163.com