重力梯度张量的拟BP神经网络反演
郭文斌，朱自强，鲁光银
(中南大学 地球科学与信息物理学院，有色金属成矿预测教育部重点实验室，湖南 长沙，410083)
摘 要：
量是反映重力场空间变化率的参数，比传统的重力异常具有更高的分辨率和更丰富的信息，将改进的BP神经网络算法应用于重力梯度张量的反演中并分析其反演效果。该算法是一种基于RPROP算法的拟BP神经网络反演算法，采用三层神经网络结构，用隐层神经元表示物性单元的密度值，根据RPROP算法自动修改各单元密度值，从而得出场源空间的密度分布。研究结果表明：采用这种算法对重力梯度张量进行反演计算，收敛速度快，对初始模型依赖性小，可准确反映出异常体形态特征和密度特征。
关键词：
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2011)12-3797-07
Quasi-BP neural network inversion of gravity gradient tensor
GUO Wen-bin, ZHU Zi-qiang, LU Guang-yin
(Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,
School of Info-Physics and Geomatics Engineering, Central South University, Changsha 410083, China)
Abstract: Based on the fact that gravity gradient tensor is a parameter which can reflect the spatial variation of gravity field, and that it has a higher resolution compared to the traditional gravity anomaly, a method for interpretation of gravity gradient tensor was proposed. The method is a kind of quasi-BP neural network algorithm which is based on RPROP algorithm. A three-layer network and the hidden layer neurons denote physics value were used. The physics value was automatically modified according to RPROP algorithm, and the physical distribution of field source was gotten. The results show that the method has a fast convergence speed and little dependence on initial model used in the inversion of gravity gradient tensor date, and can reflect the shape and density characters of anomalous body.
Key words: gravity gradient tensor; quasi-BP neural network; RPROP algorithm; inversion
与传统的重力异常相比，重力梯度张量测量的是重力矢量的梯度，能更直接地反映密度异常体的边界，对圈定矿藏范围或地质构造形态方面具有重要意义。此外，梯度测量能够克服运动载体自身加速度的影响，重力梯度张量测量被广泛应用于海洋、航空及卫星重力测量。目前，国外的重力梯度测量技术相对成熟。Bell Geospace公司、ARKeX公司和EDEX公司等都开发出全张量重力梯度测量系统，可以测量重力梯度张量的全部独立分量。超导重力梯度仪的应用使运动载体上的重力测量结果(如航空重力)的分辨率达到地面重力测量结果的分辨率，具有矿产资源勘探的能力。国内外重力梯度张量解释技术也取得了一定成就[1-8]，如：Vasco等[3]利用张量重力的对角线元素对Oklahoma西南地区的航空重力张量进行了反演；Zhang等[4]利用欧拉反褶积法反演了重力全张量；Zhdanov等[5]采用三维正则化收敛方法对重力梯度张量单分量及全张量进行了反演，并用该方法对实测的航空重力梯度张量及海洋重力梯度张量进行反演。另外，国外对重力张量的去噪、谱分析等方面的研究也取得了一定的成果[7-8]。总之，国外重力梯度张量的研究相对成熟，进入了商业应用阶段。而国内在这方面的研究只处于起步阶段，还没有重力梯度张量测量技术及相应的数据解释技术。与线性反演算法相比，非线性反演算法具有更强的全局寻优能力和反演效率，被越来越多地用于地球物理反演。在各种非线性算法中，神经网络方法由于具有独立的学习能力，被成功地用于各种地球物理方法中[7-17]。在重磁异常的数据解释方面，众多研究者[15-17]将神经网络算法和重磁异常反演理论相结合，提出了一些反演效果更好的拟神经网络算法，并在实际的数据解释中取得了很好的效果。在此，本文作者提出一种基于RPROP算法(弹性BP算法)的拟BP神经网络算法，并将其用于重力梯度张量的解释。
1 重力梯度张量
重力梯度张量[18-20]能够反映重力场的空间变化速率，与传统的重力异常相比，具有更高的分辨率，在远景区，其水平分辨率与地震反射波法的水平分辨率相当。从数学角度来讲，重力梯度张量是重力场在各个方向的导数，可用重力势的二次导数表示。在笛卡尔坐标系中，若重力势为V，则重力梯度张量可表示为一个3×3的矩阵：
由于在地球外部，位场满足拉普拉斯方程即
，并满足重力场的无旋性即 
，
，
，因此，重力梯度张量中只有5个分量是独立的。
如图1所示，在右手笛卡尔坐标系中，x轴方向指向正北方向，y轴指向东方向，z轴垂直向下，场源空间存在一剩余密度为
的长方体，其2个角点坐标分别为
和
，则在地面任意观测点(x，y，z)的重力异常及梯度张量的各分量为[19]：
(1)
(2)
(3)
(4)
(5)
(6)
(7)
其中：
G=6.67×10-11 m3/(kg3·s2)，为万有引力常量。式(2)~(7)可统一写为如下形式：
(8)
其中：α和β可为x，y和z；Sαβ称为位置函数。
图1 右手笛卡尔坐标系中的长方体模型
Fig.1 Cuboid model in right-hand Cartesian coordinate
2 拟BP神经网络原理
人工神经网络(Artificial neural network，简称ANN)[21]是通过模拟生物神经元的非线性映射功能对实际问题进行处理的一种方法，具有超强的适应能力和学习能力。在各种神经网络中，误差逆传播网络(Back propagation network，简称BP神经网络)是一种应用比较广泛和成熟的神经网络。典型的BP神经网络是3层前馈阶层网络，由输入层、隐含层和输出层构成，各层之间实行全连接。从理论上讲，神经网络方法可以很好地实现地球物理的反演，但是，在实际应用中，该反演方法过分依赖于用于网络训练的样本集，若样本集种类、数量不足，则难以取得较好的反演效果。为了改变这种情况，管志宁等[17]借鉴神经网络的权值调整方法(梯度下降法)，提出了不需要样本集的拟BP神经网络方法，并在理论模型和实测数据的三维物性反演中取得很好的效果。但是，当数据量大时，梯度下降法收敛速度慢甚至无法收敛，且易陷入局部最小，这使得在实际应用中无法对场源空间进行较精细划分。为克服这个缺点，在BP神经网络中常采用附加动量法、拟牛顿法、弹性BP法(RPROP)和Levenberg-Marquard 法[22]等。在实际应用中，拟牛顿法、Levenberg-Marquard 法收敛速度最快，但需要较大存储量，若数据量过多，则需要较高的硬件设施；附加动量法虽然收敛速度有所提高，但在数据量大时其速度仍不够快；弹性BP算法[23]收敛速度较快，所需存储量小，更适用于大数据量的计算，在此，本文采用弹性BP算法改进拟BP神经网络法，以达到更好的反演效果，并用其对重力梯度张量进行反演。
根据三维物性反演理论，用垂直于x，y和z轴的3组平面将观测面以下的场源空间划分为长方体的堆积，各长方体的体积既可以相等，也可以不等。为了便于计算，将场源空间划分为若干相等的长方体单元。设第j个物性单元的密度为ρi，据式(8)，则观测平面内第i个测点重力梯度张量各分量Gαβ,i的观测值可以表示为：
(9)
设各测点重力张量分界理论值为Ti，则式(9)可表示为：
(10)
其中：i=1，2，…，m；j=1，2，…，n；m为正演计算的总点数；n为划分的物性单元的总数。
本文使用3层神经网络，输入层为各测点实测重力张量的各分量值，隐含层表示各物性单元的密度，输出层为各测点理论值，通过输入与输出的不断比较，来调整隐层各节点的数值，达到反演的目的。假设张量的各分量实测值为fi，理论值为Ti，对实测值和理论值进行S型变换：
(11)
(12)
式中：θi为神经元的阈值；θ0为控制系数。由于S型函数具有非线性特性和可微型，这是实现BP算法的重要一步[17]。定义网络误差为：
(13)
对E求各单元密度值ρj的偏导数：
(14)
由式(13)可得：
(15)
将其代入式(14)，得：
(16)
若依据传统的拟BP神经网络方法，使用梯度下降法对各单元密度值进行调整，则当式(16)计算所得偏导数很小时，则收敛速度很慢，甚至无法收敛。在对重力梯度张量的计算中，偏导数可达到10-7，甚至更小，若依据梯度下降法对个单元进行调整，则调整量很小，这使得计算难以收敛。而RPROP 算法(弹性BP算法)权值调整幅度是独立的，只考虑导数的符号，偏导数的符号决定物性参数的更新方向，若连续2次迭代的偏导数符号不同，则说明上次调整过度权值调整幅度减少，反之，则增加。由于物性参数调整的幅值与偏导数的值无关，这种算法更易收敛。
弹性BP算法对各单元的密度值调整规则如下：
设第j个单元的第k次迭代的密度值为
，其调整值为
，则
(17)
(18)
式中：
是1个独立量，第k次迭代第j个物性单元物密度值的调整幅度仅由决定，调整的方向则由偏导数的符号决定。若2次迭代的偏导数符号相同，则表明调整的幅度不够大，需要增大调整量继续调整，以加快收敛速度；若2个迭代的偏导数符号不同，则表明调整幅度过大，密度值超过最佳值，需要减小调整量反方向调整。这种方法有效地避免了因偏导数过小导致的难以收敛问题，而且当调整量不足时，其调整幅度会持续增大，可以跳过一些局部最小点，可见：与梯度下降法相比，弹性BP法不易陷入局部最小。在数据量较大时，其优势更加明显。
每次迭代的由下式决定：
(19)
式中：0＜η-＜1＜η+；η+和η-分别为用于增加和减少调整量的系数，这两者应相互协调，不宜差别过大，否则会影响网络的收敛效果。调整量的初始值则根据期望的反演精度来确定，初始调整量应小于反演精度，但不应过小。本文使用下列目标函数作为迭代停止的条件：
(20)
式中：ε为任意正数。当满足式(20)或迭代达到一定次数时，则停止迭代。在实际计算中，往往难以选择合适的ε使反演结果达到最优，若ε选取过大，则反演效果较差；若ε选取过小，则式(20)很难被满足，网络难以收敛。本文在反演过程中选取较小的ε，以期获得更好的反演效果。为防止ε过小出现难以收敛的问题，在迭代达到一定次数且目标函数值变化平缓时也认为网络收敛，停止迭代。
为了减少多解性的影响，加入下列约束条件：当迭代一定次数后，令
(21)
在全部5个分量独立分量的联合反演(全张量反演)中，选取5个独立的分量(
)，令：
则有
(22)
式中：G和S分别为5个独立分量的观测值和相应的位置函数组成的矩阵。然后，用式(22)代替式(9)，依据相同的物质参数调整规则即式(16)，(17)及(18)，计算各个物性单元的密度ρ。
3 重力梯度张量反演
将场源空间划分为15×15×5个物性单元，每个物性单元在x和y方向上的长度为40 m，z方向上的长度为60 m，目标体模型如图2所示。该模型为2个长×宽×高都为120 m×160 m×120 m的长方体，两者之间的距离为120 m，其顶面埋深都为120 m，剩余密度都为0.31×103 kg/m3。地面的矩形测网共有15×15即225个观测点，沿x方向和y方向点距都为40 m，将通过正演公式计算其在观测面上的各重力梯度张量作为观测值。
图2 异常体模型的切片图
Fig.2 Slice of anomalous body model
3.1 单分量反演
采用前面提出的反演方法对重力梯度张量的5个独立分量分别进行反演，令ε=0.01，若式(20)成立，则认为算法收敛，迭代停止；若不满足式(20)，但迭代次数达到1 000次且目标函数值变化平缓时，也停止迭代。在各个分量的反演中，目标函数均未达到要求，迭代1 000次后，目标函数值几乎无变化，停止迭代。反演结果如图3所示。在各个分量的反演计算中，所有物性单元的密度初值都为0 t/m3，各个计算耗时都在23 s左右。
从图3可以看出：Vxx和Vyy分量较准确地反映了目标体的密度，但看不出其形态特征。在Vzz分量中，即使在增加约束条件的情况下，仍在一定区域内出现较大的负值，但负值区域与目标体位置比较接近。Vxz分量反映了异常体的大体形态特征；Vyz 分量在横向上较准确地反映了异常体的形态特征，但对其埋深反映得并不准确；Vxy分量则较为准确地反映了目标异常的形态特征，但反演得出的密度远远小于真实密度。各个分量的反演效果都不太理想。
3.2 全张量反演
3.2.1 不含噪声数据反演
在全张量反演中，目标函数仍没达到预期要求，迭代1 000次后，目标函数值几乎无变化，停止迭代。反演结果如图4所示。所有物性单元密度初值都为0 t/m3，整个计算过程耗时92.77 s。
从图4可以看出：全张量反演结果明显优于单分量的反演结果，无论目标体的形态特征还是密度都分别与真实模型的相应特征和密度十分接近。
图3 单分量反演结果
Fig.3 Inversion results of single component
3.2.2 含噪声数据反演
由于观测误差等影响，实测数据往往含有一定的噪声成分，因此，本文分别对含有10%和20%随机噪声的数据进行反演，并分析噪声对该算法的影响。其反演结果如图5所示。从图5可见：所有物性单元仍旧取密度初值0 t/m3，2次反演计算耗时都在92 s左右，与不含噪声的数据反演相比，并未增加计算时间。
图4 使用5个独立分量对不含噪声数据联合反演的结果
Fig.4 Joint inversion results of five independent components without noise
图5 对含噪声数据联合反演的反演结果
Fig.5 Joint inversion result of five independent components with noise
在含10%随机噪声的情况下，在计算过程中，使用约束条件进行约束，仍在背景区域出现了一些较大的负值，但仍准确地反映了目标体的形态及密度特征。噪声主要影响背景区域的密度参数，若对反演结果按约束条件式(21)再次进行修正，则可达到不含噪声数据的反演效果。
在含20%随机噪声的情况下，反演效果比前者的差。除了背景区出现较多负值外，目标体的形态特征(大小及位置)与模型的差距变大。
4 结论
(1) 将RPROP算法与拟BP神经网络算法相结合，提出了一种适于重力梯度张量反演计算的算法。该算法具有收敛速度快、不易陷入局部最小等特点，可用于数据量较大的反演计算。
(2) 本文使用该方法对理论模型进行了重力梯度张量的单分量反演和全张量反演。该方法耗时较少，反演过程快速。在单分量反演过程中，各个分量反演结果都较差，其中，Vxx和Vyy在一定程度上反映了目标体的密度特征，但并未反映出形态特征；而Vxy和Vyz分量较准确地反映了目标异常的形态特征，但并未反映出其密度特征。由于包含了更多的信息，全张量反演(即5个独立分量的联合反演)结果远远优于单分量反演结果，准确地反映出异常体的形态及密度特征。在反演过程中，初始模型的各单元密度初值都为0 t/m3，这表明该方法的全张量反演对初值的依赖性 较小。
参考文献：
[1] Murphy C A. Interpreting FTG gravity data using horizontal tensor components[C]//EGM 2007 International Workshop- Innovation in EM, Gravity and Magnetic methods: New Perspective for Exploration. Capri, Italy, 2007.
[2] Miller H G, Singh V. Potential field tilt: A new concept for location of potential field sources[J]. Journal of Applied Geophysics, 1994, 38(2): 213-217.
[3] Vasco D W, Taylor C. Inversion of airborne gravity gradient data, southwestern Oklahoma[J]. Geophysics, 1991, 56(1): 90-101.
[4] ZHANG Chang-you, Mushayandebvu M F, Reid A B, et al. Euler deconvolution of gravity tensor gradient data[J]. Geophysics, 2000, 65(2): 512-520.
[5] Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2006, 71(4): 925-937.
[6] de Oliveira L J C S, Tenorio L, LI Yao-guo. Efficient automatic denoising of gravity gradiometry data[J]. Geophysics, 2004, 69(3): 772-782.
[7] While J, Jackson A, Smit D, et al. Spectral analysis of gravity gradiometry profiles[J]. Geophysics, 2006, 71(1): J11-J22.
[8] Mickus K L, Hinojosa J H. The complete gravity gradient tensor derived from the vertical component of gravity: A Fourier transform technique[J]. Journal of Applied Geophysics, 2001, 46(2): 159-174.
[9] Baan M, Jutten C. Nueral networks in geophysical applications[J]. Geophysics, 2000, 65(4): 1032-1047.
[10] ZHANG Lin, Poulton M M, Wang T. Borehole electrical resistivity modeling using neural networks[J]. Geophysics, 2002, 67(6): 1790-1791.
[11] Sarzeaud O, LeQuentrec-Lalancette M F, Rouxel D. Optimal interpolation of gravity maps using a modi?ed neural network[J]. Math Geoscience, 2009(41): 379-395.
[12] Spichak V, Popova I. Artificial neural network inversion of magnetotelluric data in termsof three-dimensional earth macroparameters[J]. Geophys J Int, 2000, 42(1): 15-26.
[13] Calderon-Macias C, Sen M K, Stoffa P. Artificial neural networks for parameter estimation in geophysics[J]. Geophysical Prospecting, 2000, 48(1): 21-47.
[14] 朱自强, 程方道, 黄国祥. 同时反演2个三维密度界面的拟神经网BP算法[J]. 石油物探, 1995, 34(1): 76-85.
ZHU Zi-qiang, CHENG Fang-dao, HUANG Guo-xiang. Quasi neural network BP algorithm for simultaneous inversion of two 3-D density interfaces[J]. Geophysical Prospecting for Petroleum, 1995, 34(1): 76-85.
[15] 朱晓军, 申宁华. BP网络在位场解释中的应用[J]. 长春地质学报, 1995, 25(1): 81-86.
ZHU Xiao-jun, SHEN Ning-hua. The application of BP networks to the interpretation of potential field data[J]. Journal of Changchun University of Science and Technology, 1995, 25(1): 81-86.
[16] 张新兵, 王家林. BP、Hopfield神经网络在位场反演中的应用比较[J]. 物探与化探计算技术, 2007, 29(2): 161-166.
ZHANG Xin-bing, WANG Jia-lin. The application comparisons of BPNN and HNN in potential field inversion[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2007, 29(2): 161-166.
[17] 管志宁, 侯俊生, 黄临平, 等. 重磁异常反演的拟BP神经网络方法及其应用[J]. 地球物理学报, 1998, 41(2): 242-251.
GUAN Zhi-ning, HOU Jun-sheng, HUANG Lin-ping, et al．Inversion of gravity and magnetic anomalies using pseduo-BP neural network method and its application[J]. Chinese Journal of Geophysics, 1998, 41(2): 242-251.
[18] Droujinine A, Vasilevsky A, Evans R. Feasibility of using full tensor gradient (FTG) data for detection of local lateral density contrasts during reservoir monitoring[J]. Geophysics, 2007, 72(3): 795-820.
[19] LI Xiong, Chouteau M. Three-dimensional gravity modeling in all space[J]. Surveys in Geophysics, 1998, 19(4): 339-368.
[20] Saad A H. Understanding gravity gradients—A tutorial[J]. The Leading Edge, 2006, 25(8): 942-947.
[21] 王伟. 人工神经网络原理: 入门与应用[M]. 北京: 北京航空航天大学出版社, 1995: 52-62.
WANG Wei. The principle of the neural network: Introduction and application[M]. Beijing: Beijing Aeronautics and Astronautics Press, 1995: 52-62.
[22] 徐海浪, 吴小平. 电阻率二维神经网络反演[J]. 地球物理学报, 2006, 49(2): 584-589.
XU Hai-lang, WU Xiao-ping. 2-D resistivity inversion using the neural network method[J]. Chinese Journal of Geophysics, 2006, 49(2): 584-589.
[23] Marlin R, Heinrich B A. Direct adaptive method for faster back-propagation learning: The RPROP algorithm[C]// Proceedings of the IEEE International Conference on Neural Net-works. Beijing: IEEE Press, 1993: 586-591.
(编辑 陈灿华)
收稿日期：2010-12-06；修回日期：2011-03-10
基金项目：国家自然科学基金资助项目(41174061)；国家高技术研究发展计划(“863计划”)项目(2007AA06Z102)；中南大学自由探索计划项目(2011QNZT011)
通信作者：鲁光银(1976-)，男，湖北宜昌人，博士，副教授，从事重磁正、反演及数据处理研究；电话：13975894898；E-mail：luguangyin2006@126.com
摘要：基于重力梯度张量是反映重力场空间变化率的参数，比传统的重力异常具有更高的分辨率和更丰富的信息，将改进的BP神经网络算法应用于重力梯度张量的反演中并分析其反演效果。该算法是一种基于RPROP算法的拟BP神经网络反演算法，采用三层神经网络结构，用隐层神经元表示物性单元的密度值，根据RPROP算法自动修改各单元密度值，从而得出场源空间的密度分布。研究结果表明：采用这种算法对重力梯度张量进行反演计算，收敛速度快，对初始模型依赖性小，可准确反映出异常体形态特征和密度特征。
