重力梯度張量的擬BP神經網絡反演

郭文斌，朱自強，魯光銀

(中南大學 地球科學與信息物理學院，有色金屬成礦預測教育部重點實驗室，湖南 長沙，410083)

與傳統的重力異常相比，重力梯度張量測量的是重力矢量的梯度，能更直接地反映密度異常體的邊界，對圈定礦藏范圍或地質構造形態方面具有重要意義。此外，梯度測量能夠克服運動載體自身加速度的影響，重力梯度張量測量被廣泛應用于海洋、航空及衛星重力測量。目前，國外的重力梯度測量技術相對成熟。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的矩陣：

由于在地球外部，位場滿足拉普拉斯方程即Vxx+Vyy+Vzz=0，并滿足重力場的無旋性即Vxy=Vyx，Vxz=Vzx，Vzy=Vyz，因此，重力梯度張量中只有5個分量是獨立的。

如圖1所示，在右手笛卡爾坐標系中，x軸方向指向正北方向，y軸指向東方向，z軸垂直向下，場源空間存在一剩余密度為的長方體，其 2個角點坐標分別為 (ξ1,η1,ζ1)和 (ξ2,η2,ζ2)，則在地面任意觀測點(x，y，z)的重力異常及梯度張量的各分量為[19]：

圖1 右手笛卡爾坐標系中的長方體模型Fig.1 Cuboid model in right-hand Cartesian coordinate

G=6.67×10-11m3/(kg3s2)，為萬有引力常量。式(2)～(7)可統一寫為如下形式：

其中：α和β可為x，y和z；Sαβ稱為位置函數。

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的觀測值可以表示為：

設各測點重力張量分界理論值為Ti，則式(9)可表示為：

其中：i=1，2，…，m；j=1，2，…，n；m為正演計算的總點數；n為劃分的物性單元的總數。

本文使用3層神經網絡，輸入層為各測點實測重力張量的各分量值，隱含層表示各物性單元的密度，輸出層為各測點理論值，通過輸入與輸出的不斷比較，來調整隱層各節點的數值，達到反演的目的。假設張量的各分量實測值為fi，理論值為Ti，對實測值和理論值進行S型變換：

式中：θi為神經元的閾值；θ0為控制系數。由于S型函數具有非線性特性和可微型，這是實現 BP算法的重要一步[17]。定義網絡誤差為：

對E求各單元密度值ρj的偏導數：

若依據傳統的擬 BP神經網絡方法，使用梯度下降法對各單元密度值進行調整，則當式(16)計算所得偏導數很小時，則收斂速度很慢，甚至無法收斂。在對重力梯度張量的計算中，偏導數可達到10-7，甚至更小，若依據梯度下降法對個單元進行調整，則調整量很小，這使得計算難以收斂。而RPROP 算法(彈性BP算法)權值調整幅度是獨立的，只考慮導數的符號，偏導數的符號決定物性參數的更新方向，若連續2次迭代的偏導數符號不同，則說明上次調整過度權值調整幅度減少，反之，則增加。由于物性參數調整的幅值與偏導數的值無關，這種算法更易收斂。

彈性BP算法對各單元的密度值調整規則如下：

設第j個單元的第k次迭代的密度值為，其調整值為，則

式中：Δkj是1個獨立量，第k次迭代第j個物性單元物密度值的調整幅度僅由Δkj決定，調整的方向則由偏導數的符號決定。若2次迭代的偏導數符號相同，則表明調整的幅度不夠大，需要增大調整量繼續調整，以加快收斂速度；若2個迭代的偏導數符號不同，則表明調整幅度過大，密度值超過最佳值，需要減小調整量反方向調整。這種方法有效地避免了因偏導數過小導致的難以收斂問題，而且當調整量不足時，其調整幅度會持續增大，可以跳過一些局部最小點，可見：與梯度下降法相比，彈性 BP法不易陷入局部最小。在數據量較大時，其優勢更加明顯。

每次迭代的Δkj由下式決定：

式中：0＜η-＜1＜η+；η+和η-分別為用于增加和減少調整量的系數，這兩者應相互協調，不宜差別過大，否則會影響網絡的收斂效果。調整量的初始值則根據期望的反演精度來確定，初始調整量應小于反演精度，但不應過小。本文使用下列目標函數作為迭代停止的條件：

式中：ε為任意正數。當滿足式(20)或迭代達到一定次數時，則停止迭代。在實際計算中，往往難以選擇合適的ε使反演結果達到最優，若ε選取過大，則反演效果較差；若ε選取過小，則式(20)很難被滿足，網絡難以收斂。本文在反演過程中選取較小的ε，以期獲得更好的反演效果。為防止ε過小出現難以收斂的問題，在迭代達到一定次數且目標函數值變化平緩時也認為網絡收斂，停止迭代。

為了減少多解性的影響，加入下列約束條件：當迭代一定次數后，令

在全部 5個分量獨立分量的聯合反演(全張量反演)中，選取5個獨立的分量(Vxx,Vzz,Vxy,Vxz,Vyz)，令：

式中：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×103kg/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。

圖3 單分量反演結果Fig.3 Inversion results of single component

從圖4可以看出：全張量反演結果明顯優于單分量的反演結果，無論目標體的形態特征還是密度都分別與真實模型的相應特征和密度十分接近。

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 fi 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.