基于深度學習的重力異常與重力梯度異常聯合反演

張志厚， 廖曉龍， 曹云勇， 侯振隆， 范祥泰， 徐正宣,，路潤琪， 馮濤， 姚禹， 石澤玉

1 西南交通大學地球科學與環境工程學院， 成都 611756 2 西南交通大學高速鐵路線路工程教育部重點實驗室， 成都 610031 3 中鐵二院成都地勘巖土工程有限責任公司， 成都 610000 4 東北大學深部金屬礦山安全開采教育部重點實驗室， 沈陽 110819 5 東北大學資源與土木工程學院， 沈陽 110819

0 引言

重力物性反演是典型的線性反問題，存在多解性，設計合理的反演算法是關鍵步驟之一.為獲得穩定、準確的反演結果，一些正則化方法被引入到重力數據的線性反演中，并取得了較好的計算結果.如：基于地質異常體最小體積約束方法及其改進方法(Last and Kubik, 1983; Guillen and Menichetti, 1984; Barbosa and Silva, 1994; Li and Oldengburg, 1998; Silva and Barbosa, 2006; Dias et al., 2009, 2011)和基于最小梯度約束的“聚焦反演”(Portniaguine and Zhdanov, 1999; Zhdanov et al., 2004; Zhdanov, 2009)等.此外，為了抑制反演的“趨膚效應”(Fedi et al., 2005)，將深度加權函數引入到目標函數中進行反演(Li and Oldengburg, 1998; Chasseriau and Chouteau, 2003; Commer, 2011)的方法在一定程度上提高了深度分辨率，以及加入其他地質先驗信息或密度模型的重力反演方法(Shamsipour et al., 2010)也取得了較好的效果.隨著FTG(Full Tensor Gravity Gradiometry)系統的問世與成功應用，以上諸多方法也被應用在重力梯度數據的反演中(Martinez et al., 2013; Geng et al., 2014; Lu and Qian, 2015; 高秀鶴等，2017).由于重力梯度數據能夠提供更豐富的地質信息，因此利用重力及重力梯度數據聯合反演不僅能夠更加準確的勾畫出地下異常體的物性大小與空間分布特征，而且聯合反演也是約束多解性的一種重要手段(Gallardo and Meju, 2003; Fregoso and Gallardo, 2009; Zhdanov et al., 2012).于是重力及重力梯度數據的聯合反演方法與技術研究備受諸多學者關注(Wu et al., 2013; Qin et al., 2016; 秦朋波和黃大年，2016).但線性類反演方法通常容易陷入局部極小，并且計算效率低、穩定性差(徐世浙等，2009).在實際應用中，只有選擇良好的初始模型，線性類反演方法才可能獲得較為可靠的結果，然而初始模型應包含有關地質體屬性的先驗信息，對于空間復雜的地質情況，通常很難對其進行可靠的估算.

大數據時代的到來，人工智能的興起，聚類(Sun and Li, 2015；Liu and Jin, 2020)、決策樹(戴詩華等, 2010；Barak et al,. 2020)、支持向量機(Li et al., 2020a)、貝葉斯(Guo et al., 2011; 殷長春等, 2014; Jiang et al., 2016)及人工神經網絡(Artificial Neural Network, ANN)(Al-Garni, 2009, 2013, 2015； Neyamadpour et al., 2010)等機器學習方法被引入到地球物理數據的處理中.利用機器學習的方法進行反演無需依賴初始模型，也不需要對靈敏度矩陣進行計算和存儲，并且有著良好的全局尋優性能，此類方法已逐步成為地球物理數據處理的研究熱點(張旗和周永章, 2017).尤其是機器學習的一個分支——深度學習(Deep Learning，DL)(Schmidhuber, 2015; Lecun et al., 2015)或稱深度機器學習(Deep Machine Learning)，在不適定反問題方面取得了引領性的進展(Girshick et al., 2014; Russakovsky et al., 2015; Krizhevsky et al., 2017).DL是一種深度的神經網絡，能夠用非監督式或半監督式的特征學習和分層特征提取高效算法來代替手工獲取特征(王昊等，2020).DL相比傳統的ANN具有更多的隱含層和更強的復雜函數逼近性能(Hornik et al., 1989; Lecun et al., 2015)，DL與地球物理反演方法具有類似的優化算法(Kim and Nakata, 2018),極有助于不適定反問題的求解(Jin et al., 2017; Liu et al., 2020)，也將在地球物理的海量數據反演中取得巨大的優勢(王昊等, 2020).

目前，深度學習在地球物理反演中的應用大都是在電磁數據及地震數據處理領域.最新的研究如：王鶴等(2018)利用遺傳神經網絡實現了大地電磁反演，并驗證了該算法的可行性和有效性；Puzyrev(2019)利用全卷積神經網絡(Fully Convolutional Neural Networks，FCN)對井中垂直磁偶極子源發射地面接收的電磁數據進行二維反演；Puzyrev和Swidinsky(2020)將DL應用于海洋頻率域可控電磁數據的一維反演，其結果表明，在淺層存在異常體時，深度學習反演相比傳統反演，其不會影響深層的目標探測，具有較高的分辨率；Noh等(2020)實現了頻率域航空電磁數據卷積神經網絡(Convolutional Neural Networks，CNN)一維反演，結果表明CNN在合成數據量較少時其反演結果仍然具有較高的分辨率；Liu等(2020)提出了一種基于平滑約束和深度加權的直流電阻率FCN反演方法，有效的提高了解的分辨率；Li等(2020b)利用FCN直接從地震記錄數據重建速度模型，其“端到端”的反演思路驗證了DL在地球物理反演方法中的有效性；奚先和黃江清(2018)實現了散射波場的CNN反演成像方法，該方法具有良好、穩定的反演能力和泛化能力.此外，基于DL的地震全波形(Moseley et al., 2019; Zhang and Alkhalifah, 2019; Sun and Demanet, 2020)反演研究也取得了良好的進展.這些研究為重磁反演帶來了更為廣泛的應用前景.

基于此，本文提出一種基于DL的重力及重力梯度數據聯合反演的方法.反演流程包括三個階段：樣本數據集生成，網絡模型訓練和利用網絡模型進行預測.前兩個階段計算耗時較長，因為必須在足夠大的訓練集上訓練網絡模型，才能使網絡模型具有泛化新數據的能力，但前兩個階段僅執行一次.反演預測較為迅速，通常在幾十毫秒內完成.針對反演第一階段，本文首先提出了一種重力及重力梯度空間域快速正演方法，以此來快速構建樣本數據集；對于DL學習的網絡結構，本文基于U-net網絡結構(Russakovsky et al., 2015)建立重力數據(輸入)到剩余密度模型(輸出)的非線性映射網絡結構，即三維重力場反演的“端到端”網絡結構，本文稱為GraInvNet (Gravity Inverse Net)網絡結構；然后利用樣本數據集對網絡進行監督學習(訓練)，最后利用GraInvNet網絡對理論數據和實測重力及重力梯度數據進行反演，以此來檢驗文中方法的可行性和有效性.

1 重力及重力梯度空間域快速正演

利用地下半空間模型的重力及重力梯度正演計算生成樣本數據集，以及擴大樣本數據集是減少DL反演過擬合的有效方法.因此，本文提出一種快速的重力異常及重力梯度異常的正演計算方法，以此來高效構建樣本數據集，從而提高方法的實用性.

1.1 正演計算表達式

通常重力及重力梯度的正演計算是將地下半空間剖分成若干個規則長方體模型單元(如圖1a所示)，然后計算每一個長方體模型單元對觀測點的異常(圖1b)，再將所有長方體模型對觀測點的異常求和，即為地下半空間整個模型體對觀測點的異常.

圖1 地下半空間剖分模型單元示意圖(a) 剖分單元示意圖； (b) 長方體模型單元.Fig.1 Sketch map of partition model of underground half-space subdivision(a) Sketch map of subdivision element; (b) Prism model element.

單個長方體模型單元相對觀測點P(x,y,z)的重力異常及重力梯度異常正演計算理論表達式如下(Li and Chouteau, 1998; 舒晴等, 2015)：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

1.2 等效存儲網格點幾何格架技術

姚長利等(2003)提出了幾何格架技術進行了重力異常的空間域快速正演計算，取得較好的效果.本文在此基礎上提出一種網格點幾何格架技術，此技術能夠對已有的正演算法進行加速.以重力異常的正演計算為例進行說明，公式(1)可以寫成如下表達式：

(8)

其中，

Q(xi,yj,zk)=-G0[xiln(yj+rijk)+yjln(xi+rijk)

(9)

與物性無關，本文稱為網格點幾何格架函數.網格點幾何格架函數與文獻(姚長利等, 2003)類似的幾何格架函數：

S(xi,yj,zk)=

(10)

存在如下關系：

(11)

分析得出：(1)幾何格架S與網格點幾何格架Q的區別在于，前者是觀測點相對于長方體模型單元的概念，后者是觀測點相對于地下半空間的模型剖分網格點的概念；(2)幾何格架S與網格點幾何格架Q都具有平移等效性，重力異常Δg及重力梯度張量Uzz的幾何格架S具有互換對稱性；(3)幾何格架S值是長方體模型單元8個頂點的網格點幾何格架Q的代數和；(4)直接計算不同幾何格架S的值存在8次(圖2)計算同一網格點幾何格架Q值的過程(邊界除外).

圖2 長方體模型單元共用網格點幾何格架示意圖Fig.2 Sketch map of geometric grid of common grid points of prism model elements

基于以上分析，本文對部分觀測點網格點幾何格架Q的值進行計算、存儲，其余觀測點網格點幾何格架的值通過平移等效性進行調用，而不需要計算與存儲，從而通過公式 (8)完成正演計算.

1.3 計算步驟

下面給出基于網格點幾何格架的重力及重力梯度異常的空間域快速正演算法具體步驟.

步驟1：將地下半空間長方體模型單元的剖分與計算點相對應(圖3所示)，即計算點與地下半空間模型單元的網格點一一對應.并將計算點位置換算成實際里程，對長方體模型單元的剩余密度進行賦值；

步驟2：分別計算地下半空間內所有網格點對長方形觀測區4個頂點的重力異常或重力梯度異常的網格點幾何格架值(分別記為Q1，Q2，Q3及Q4.圖3所示)，并存儲待后續計算調用；

步驟3：根據計算點與任意單一長方體模型單元8個頂點之間的相對位置關系，利用平移等效性調取其相對應的網格點幾何格架值，再依據公式 (8)計算該長方體模型單元對計算點的重力異常或重力梯度異常；

步驟4：利用步驟3計算地下半空間每一個長方體模型單元對某一計算點的重力異常或重力梯度異常，最后求和得到該計算點總的重力異常或重力梯度異常；

步驟5：對步驟4遍歷循環觀測區，最終獲得地下半空間剩余密度體在觀測區的重力異常或重力梯度異常.

圖3 計算點與網格點一一對應關系Fig.3 The corresponding relation between calculation points and grid points

2 深度學習反演

2.1 方法

目前DL網格結構主要有CNN、深度置信網絡(Deep Belief Network, DBN)、循環神經網絡(Recurrent Neural Networks, RNN)及生成對抗網絡(Generate Adversarial Network, GAN)等，其余網絡結構都是在此基礎上改進而來.而應用最為廣泛的是CNN，CNN在大型圖像處理方面取得了一系列新的突破(Krizhevsky et al., 2017)，其網絡結構具備較強泛化能力和較好的遷移能力，以及健壯性也會更加優質(LeCun et al., 2015)，尤其是FCN作為特殊的一種CNN其在醫學圖像語義分割方面(Russakovsky et al., 2015)顯現出超越性的優勢.

本文將觀測平面上的多元重力數據視為二維圖像的像素數據，作為FCN網絡模型的輸入端，將已知的地下模型作為訓練標簽，并在搭建的FCN網絡模型上進行監督學習(訓練)，最終通過訓練完成的FCN網絡模型對重力數據進行反演.這一點本文反演方法與圖像語義分割類似，但也存在差別.圖像語義分割(圖4a、b)的輸入與輸出的空間維度相同，本文反演方法的輸入與輸出的空間維度(圖4c、d)不同，即輸出端的剩余密度為三維數據(圖4d)(地下半空間異常體物性值)，輸入端重力場響應為二維數據(圖4c).因此，本文FCN網絡結構的輸出端采用多通道二維數據來連接分類模型，即將地下半空間每一深度層的物性值看作為自然圖像的二維像素數據，多通道的通道數為反演深度的網格單元數.

圖4 FCN端對端示意圖(a)與(b)分別為圖像分割的輸入與輸出示意圖，(c)與(d)分別為本文方法的輸入與輸出示意圖.Fig.4 Sketch map of FCN end-to-end(a) and (b) are the input and output schematic of image segmentation, (c) and (d) are the input and output schematic of the method proposed in this paper respectively.

地下半空間存在的密度異常體會在觀測面上引起重力場響應，不同異常體引起的重力場響應也表現出不同特征(圖4c、d).即具有一定的空間相關性和局部存在性，此特征也可以從極為發展的位場邊緣識別技術(王萬銀, 2010；Wang et al., 2009, 2010)得到印證.總之，位場數據反演的輸入和輸出有兩個特點：(1)空間的對應性，(2)局部的存在性.而FCN方法作為一種特殊的DL方法，其利用卷積算子代替矩陣相乘，重點學習輸入圖像與輸出標簽的局部性和空間性(Shi et al., 2019).因此，利用FCN方法對重力數據進行反演是切實可行的.

相同大小的異常體位于兩個鄰近垂直位置時，其引起的重力異常變化較小，FCN的網絡結構需要具有足夠的深度才有可能將其區分.然而，網格結構加深則會造成訓練難度增加，諸如反向傳播梯度消失、計算量增加等.因此，有必要將以高頻信息為主的重力梯度異常作為以低頻信息為主的重力異常的補充，共同作為FCN的輸入端，可以增加同大小異常體位于兩鄰近垂直位置時的重力場響應差異，從而進一步提高FCN的預測精度.換言之，在不增加感受野的前提下，提高位置相似異常體的重力響應的辨識度來進行FCN網絡訓練，以此提高其預測的精度.本文反演的目標是通過FCN的方法得到重力異常和重力梯度異常到異常體剩余密度的映射函數：

ρ=F(Δg,Uxx,Uxy,Uxz,Uyy,Uyz,Uzz).

(12)

2.2 網絡結構

Long等(2015)提出了FCN方法，該方法是將CNN方法的全連接層修改為卷積層，實現了端到端的學習.FCN方法相比傳統的CNN方法，其優點一是節點超參數大大減少從而提高計算效率，二是避免全連接層數據一維化而丟棄信息，三是“跳躍結構”(Russakovsky et al., 2015)融合了低中高階特征信息，進一步優化了輸出結果，以此可以高效進行密集型像素的預測.本文基于FCN方法的U-Net架構設計重力場數據反演的網絡結構(GraInvNet)，如圖5所示.GraInvNet是U-Net的一種簡化與改進，簡化了網絡結構的深度及卷積通道數，改進了輸出端的標簽分類(圖4d)，以及當輸入與輸出在不同空間維度時，利用第一次卷積進行裁剪或填充，以便跳躍鏈接時編碼層與解碼層的組合.GraInvNet有15層采用卷積運算(寬度為3，跨度為1)，3層采用最大池化運算(寬度為2，跨度為2)，其感受野為80(計算方法參照Luo 等(2017)文獻)，這對于本文重力場輸入數據空間維度的反演任務來說，非線性逼近程度已足夠大.

2.3 正則化

在DL中，泛化能力是指已完成訓練的網絡結構對未泄露數據的預測性能(Zhang et al., 2018; Chan et al., 2017).在實際情況中，我們通常將樣本數據集按照一定比例分成兩部分，大部分用來訓練，小部分未泄露數據進行測試，最后通過測試誤差來評價方法的泛化能力.決定DL反演成敗的兩個關鍵因素是(1)降低訓練誤差和(2)減少訓練誤差與測試誤差之間的差距.訓練誤差較大對應為欠擬合，可以更改FCN網格結構或訓練算法來完善，訓練誤差和測試誤差之間的較大差距對應過擬合.在數據量充足的情況下正則化方案是解決過擬合的一種有效方法.

Dropout是深度學習模型訓練中最簡單、最常用的正則化技術(Hinton et al., 2012b).Dropout是指在每一輪訓練中按一定比例隨機丟棄部分網絡神經單元之間的鏈接(Srivastava et al., 2014; Krizhevsky et al., 2017)，避免訓練數據在所有節點進行訓練，該策略減少了網絡單元之間復雜的共適應關系，同時也對網絡單元起一定的平均作用，故而可減少過擬合.此外，本文對于GraInvNet網絡結構，不僅在卷積層后使用Dropout，還在卷積層間使用批處理歸一化的正則化方案.該方案的優勢是在FCN網絡訓練期間性能更穩定(Ioffe and Szegedy, 2015; He et al., 2016; Puzyrev, 2019).

2.4 損失函數

損失函數是DL網絡訓練過程的指揮棒，其通過正向傳播的輸出結果與樣本標簽值之間的誤差反向傳播來修改網絡單元節點權值及其他參量，本文指通過預測地下半空間長方體模型密度與已知模型密度的誤差進行調參.常用損失誤差為L1、L2范數，其表達式分別為

圖5 GraInvNet網絡結構由3個下采樣層和3個上采樣層組成，下采樣和左側卷積層對應編碼階段，上采樣和右側卷積層對應解碼階段，圖中黑色數字表示特征圖的通道數，Conv為卷積，BN為批歸一化，ReLU為激活函數，Max pooling為最大池化，Upsampling為上采樣， Skip-Connection為跳躍鏈接.Fig.5 GraInvNet architectureIt is composed of three down-sampling layers and three up-sampling layers. The down-sampling layer and the left convolution layer correspond to the encoding stage, the up-sampling layer and the right convolution layer correspond to the decoding stage. In the figure, the black number indicates the number of channels in the feature map, Conv is convolution, BN is batch normalization, ReLU is activation function.

(13)

GraInvNet網絡模型訓練的實質就是損失函數數學期望極小化的過程.傳統地球物理反演的實質是目標函數極小化的過程.其實兩者的二范數表達式類似，區別在于傳統反演目標函數是預測數據的正演響應與實測數據的二范數或其正則化改進，且該目標函數極小化過程的計算成本極為昂貴，而DL損失函數極小化過程效率較高.本文采用歸一化數據的L2范數作為損失函數.

2.5 學習算法的選擇

強大的學習算法是提高深度學習反演性能必不可少的步驟，通過GraInvNet網絡學習算法一步步的迭代使損失函數極小化，達到全局的最優值.

隨機梯度下降法(Stochastic Gradient Descent，SGD)是一種高效的優化方法，是DL中最有效的核心算法之一(Hinton et al., 2012a; Deng et al., 2013; Hinton and Salakhutdinov, 2006; Krizhevsky et al., 2017).然而SGD選取一個樣本數據作為其梯度方向時，其方向具有隨機性，極有可能不能較快收斂到最優解.此外，SGD算法的學習率是全局統一的，而待優變量對于損失函數的依賴卻各不相同，因此其取值不當易造成部分變量收斂速度慢及部分變量的收斂過程產生震蕩.針對SGD固定的學習率問題，AdaGrad (Adaptive Gradient)算法將每一個變量采用不同的學習率，自適應進行學習率的調整.但該算法存在一個平方梯度迅速積累而造成學習速度下降的問題，從而阻止訓練模型達到局部極小值.RMSProp (Kingma and Ba, 2014)是AdaGrad的改進算法，其采用平方梯度累積指數衰減的動態平均算法，使得模型訓練得以繼續.考慮到動量法具有對梯度算法加速的優點，本文采用Kingma和Ba(2014)提出的ADAM算法(Adaptive Moment Estimation)，它是動量法與RMSProp算法的結合，適合大數據量超參數的優化問題，并具備計算效率高收斂速度快的優勢.

2.6 網絡規格及調參

激活函數選用CNN以及其他DL類型最常用的線性整流函數(ReLU，見公式(14))及其改進函數(Xu et al., 2015)，其主要優勢是稀疏性，能夠減少梯度消失的可能性，還具有非常好的收斂性能和較低的計算成本.輸出選擇S型函數(見公式(15)).此外，鑒于在圖形特征提取中最大值池化相比平均值池化具有一定優勢(Boureau et al., 2010)，本文GraInvNet網絡選用最大值池化方案.

fReLU(t)=max(t,0),

(14)

fsigmoid(t)=1/(1+e-t).

(15)

在調參方面需要考慮到訓練輪數、樣本批量數以及Dorpout的比例等因素對模型性能的影響.以上因素需兼顧計算效率和精度，從而優化GraInvNet網絡結構，最終的訓練完成的GraInvNet具備計算的準確性，同時也需對新樣本具備一定的泛化性.

3 理論模型試驗

3.1 樣本數據集的構建及訓練

將地下半空間劃分為40×40×20個長方體模型單元，如圖1a所示，其中長方體模型單元大小為100 m×100 m×50 m，模擬觀測數據區大小為4 km×4 km，計算的單通道(重力異常或某種分量重力梯度異常)正演數據為41×41個.樣本數據集的構建首先是將多個長方體模型單元組合為一定規模的長方體模型塊體或臺階模型，如8個長方體模型單元組合為2×2×2長方體模型塊體，64個長方體模型單元組合為4×4×4長方體模型塊體等；隨后將不同大小的長方體模型塊體分別設置在地下半空間不同位置，再利用本文第1節提供的快速正演方法獲取其重力異常及重力梯度異常數據，其中剩余密度大小為1.0 g·cm3；最后對歸一化后正演數據作為樣本數據集中的輸入數據，不同位置處的密度歸一化后作為標簽.樣本數據集總共包含了14696對不同組合模型塊體正演數據，單個樣本數據輸入端包含7個通道，即Δg，Uxx，Uxy，Uxz，Uyy，Uyz和Uzz異常，數據規格為41×41×7，輸出端為40×40×20，即地下半空間長方體模型單元的個數.

由于重力異常及重力梯度異常的正演為線性關系(式(1)—(7))，故可利用已正演的樣本數據進行不同密度大小的樣本數據擴充，如本文剩余密度大小為-1.0 g·cm3的樣本數據是直接擴充而獲取的，加上正演數據集，共計29392個樣本數據.將樣本數據劃分為訓練數據、測試數據以及驗證數據，其對應比為18∶1∶1，即訓練數據集包含26452個樣本數據，而驗證和測試集各有1470個樣本數據.

訓練輪數為100，每次訓練后進行一次驗證，圖6所示為本文所提GraInvNet網絡的訓練誤差與驗證誤差，可以看出，模型訓練良好，未出現過欠擬合.

圖6 GraInvNet訓練誤差及驗證誤差Fig.6 Training error and validation error of GraInvNet

圖7 組合模型透視圖Fig.7 Perspective view of combination model

3.2 組合模型一

長方體組合模型一由4個不同大小長方體模型組成，如圖7所示，長方體剩余密度都為1.0 g·cm-3.其余參數見表1.利用文中已訓練完成的GraInvNet網絡對組合模型一分別進行重力異常反演、重力異常與部分重力梯度異常聯合反演以及重力異常與全部重力梯度異常聯合反演的三種反演方案進行對比，其結果如圖8a、b及c所示.可以看出在無噪聲情況下，重力異常反演結果(圖8a)能夠“聚焦”于真實模型，表明本文GraInvNet網絡不僅可以有效的對重力異常進行反演，而且還說明了對單一塊體模型所構建的樣本數據集進行訓練，可對類似塊體的組合模型進行預測，驗證了FCN的卷積核具有局部性和權值共享的特點和優勢，從而表明了本文GraInvNet網絡所具備的泛化性能.但同時也可以看出單獨重力異常反演結果的邊界不夠清晰，而重力異常與部分重力梯度異常組合的多維度信息作為GraInvNet網絡的輸入端進行反演，其結果(圖8b)擬補了重力異常反演結果邊界模糊的劣勢，重力異常與全部重力梯度異常的聯合反演結果(圖8c)邊界更為清晰，精度更高，與理論模型輪廓保持高度一致，說明了基于DL的多維度聯合反演精度相比單一的重力異常反演結果精度高.

表1 長方體模型參數Table 1 Cuboid model parameters

為了更加貼合實際資料處理，對重力異常及重力梯度異常分別加入其最大值10%的隨機噪聲，并對其進行反演，結果如圖8d、e及f所示.可以看出，含噪聲重力異常反演結果(圖8d)雖然可以反映出異常體的位置，但相比無噪聲情況下重力異常反演結果(圖8a)，其效果較差；含噪聲的重力異常與部分重力梯度異常的聯合反演結果(圖8e)對比無噪聲情況下的反演結果(圖8b)同樣受到一定的影響，但聯合反演相比單獨重力異常反演結果其在邊界處的影響更小，聯合反演的結果可以準確的反映出異常體的輪廓.含噪聲重力異常與全部重力梯度異常的聯合反演結果(圖8f)基本上與無噪聲情況下的反演結果(圖8c)一致.以上說明了基于DL的多維度聯合反演具有一定抗噪聲性能，計算穩定性較高.

為了進一步對比三種反演方案的優劣，給出了坐標點(0 km，0 km)至坐標點(4 km，4 km)剖面的切塊圖(圖9)，以及z=0.40～0.45 km反演結果水平切塊圖(圖10).可以看出無噪聲情況下重力異常反演結果縱向分辨率較差(圖9a)，橫向分辨率相對較好，但也存在邊界的模糊性(圖10a，編號為4的異常體).含噪聲情況下重力異常反演結果橫縱向分辨率都變差，橫向的邊界不夠清晰(圖10d)，異常體的縱向埋藏深度出現較大的偏差(圖9d，編號為1和2的異常體).這是由于重力異常為低頻信息，同姿態同大小的異常體埋深極為相近時，其在觀測面引起的重力異常變化較小，一方面，隨機噪聲的大小極有可能大于其垂向相鄰異常體在觀測面上重力異常的變化，從而造成GraInvNet網絡預測結果存在深度上的偏差；另一方面如果采用DL去準確區分觀測重力異常差異較小的垂向相鄰異常體時，則需要更深的網絡結構和更高的計算成本；此外，為了防止過擬合則需要更廣泛和更充足的樣本數據集，從而進一步增加訓練難度和成本.因而，為了提高DL的反演的精度，極為有必要增加同姿態同大小異常體位于相近位置時重力場響應的差異.而多維度重力梯度異常和重力異常能夠多角度增加異常體響應的辨識度，其全方位的重力場響應數據作為DL的訓練樣本數據，在一定程度上會減少“盲人摸象”的困境，進一步溯源求真.如圖9b、e與圖10b、e所示的部分重力梯度異常聯合反演結果相比重力異常的反演結果在縱向上的精度大大提高，重力異常與全部重力梯度異常的聯合反演結果(圖9c、f與圖10c、f)更為精確，橫縱分辨率高.

圖8 反演結果三維透視圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果，(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果.Fig.8 3D perspective views of inversion result(a), (b) and (c) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the absence of noise, respectively; (d), (e) and (f) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the presence of noise, respectively.

圖9 反演結果切塊圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果；(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果；切塊線首尾坐標分別為(0 km，0 km)和(4 km，4 km)，白色虛線框為模型體的邊界位置.Fig.9 Block diagram of inversion result(a), (b) and (c) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the absence of noise, respectively; (d), (e) and (f) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the presence of noise, respectively; The first and last coordinates of the cutting line are (0 km，0 km) and (4 km，4 km)，the white dotted box indicates the boundary position of the model.

圖10 z=0.40～0.45 km反演結果切塊圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果；(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結果；白色虛線框為模型體的邊界位置.Fig.10 Block diagram of inversion result at z=0.40～0.45 km(a), (b) and (c) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the absence of noise, respectively; (d), (e) and (f) are the inversion results of Δg, Δg & Uxz & Uyz & Uzz and Δg & all components in the presence of noise, respectively；the white dotted box indicates the boundary position of the model.

為了進一步說明聯合反演結果與理論異常的擬合程度，文中分別將無噪聲與含噪聲數據的反演結果進行正演，其正演重力異常結果與理論異常的殘差值如圖11所示，重力異常與重力梯度異常殘差的標準差與隨機噪聲標準差見表2.無噪聲情況下重力異常差值的最大最小值分別為0.912、0.015 mGal，含噪聲情況下重力異常差值的最大最小值分別為1.613、0.008 mGal，經對比分析各分量的方差及標準差，表明文中所述反演方法的精度高.

3.3 組合模型二

組合模型二為四層臺階與長方體的組合模型，如圖12所示，其中四層臺階模型的每層臺階大小相同都為1.4 km×0.6 km×0.15 km，每層臺階的中心坐標分別為(2.0 km,0.9 km,0.275 km)、(2.0 km,1.2 km,0.425 km)、(2.0 km,1.5 km,0.575 km)及(2.0 km,1.8 km,0.725 km)，長方體大小為0.8 km×0.5 km×0.3 km，長方體頂面埋深為0.2 km，中心點坐標為(2.0 km,3.05 km,0.35 km)，四層臺階與長方體剩余密度都為1.0 g·cm-3.組合模型二相比組合模型一的變化是加入新類型模型——臺階，以及樣本數據集中未包含的長方體模型.利用本文的GraInvNet網絡對組合模型二的重力異常及全部重力梯度異常進行聯合反演，其結果如圖13所示，可以看出臺階模型的反演結果與真實模型高度吻合，其幅值在深部處略有差別(圖13c),長方體模型在深部的幅值也略有差別(圖13e及f)，其余部分也與理論模型高度吻合.表明GraInvNet網絡可以對與樣本模型相似組合模型的異常體進行精確定位，同時也可以準確地預測出模型值的大小，表明其反演能力很好.為了更加直觀地了解反演結果中異常體的位置、形狀和物性值，圖13a、b、c與d 分別給出不同方向的反演切塊圖，可以看出反演結果與真實模型的邊界及形狀吻合，并且反演結果的幅值大小也與已知模型接近.

圖11 重力理論異常與預測結果正演的殘差(a)無噪聲情況; (b)含噪聲情況.Fig.11 The residuals of theoretical gravity anomaly and forward value of prediction result(a) Free noise situation; (b) Noise situation.

表2 理論異常最大最小值、反演結果擬合殘差的標準差及隨機噪聲標準差Table 2 The maximum and minimum value of theoretical anomaly, the standard deviation of the fitting residuals of inversion result and random noise

圖12 組合模型透視圖Fig.12 Perspective view of combination model

同時為了更加貼合實際資料處理，對重力異常及重力梯度異常分別加入其最大值10%的隨機噪聲，再對其進行反演，結果如圖14所示.可以看出，除去異常體部分邊界位置的幅值略有差異(圖14c、e及f)，其余部分也與理論模型高度吻合，表明GraInvNet網絡反演的穩定性較高，具有一定的抗噪聲能力.

為了進一步說明聯合反演結果與理論異常的擬合程度，文中分別將無噪聲與含噪聲數據的反演結果進行正演，其正演重力異常結果與理論異常的殘差值如圖15所示，重力異常與重力梯度異常殘差的標準差與隨機噪聲標準差見表3.無噪聲情況下重力異常差值的最大最小值分別為0.985、0.012 mGal，含噪聲情況下重力異常差值的最大最小值分別為2.975、-0.864 mGal，經對比分析各分量的標準差與理論異常的大小，表明文中所述反演方法的精度高.

表3 理論異常最大最小值、反演結果擬合殘差的標準差及隨機噪聲標準差Table 3 The maximum and minimum value of theoretical anomaly, the standard deviation of the fitting residuals of inversion result and random noise

圖13 無噪聲反演結果(a)、(b)、(c)及(d)分別為y=1.3～1.4 km、x=1.9～2.0 km、x=2.5～2.6 km及z=0.35～0.4 km的反演結果切塊圖，白色虛線框為模型體的邊界位置，(e)、(f)為剩余密度大于0.1 g·cm-3的三維透視圖.Fig.13 Inversion result in the absence of noise(a), (b), (c) and (d) are block diagrams of inversion results at y=1.3～1.4 km, x=1.9～2.0 km, x=2.5～2.6 km and z=0.35～0.4 km respectively, the white dotted box indicates the boundary position of the model, (e) and (f) are the three-dimensional perspectives, where residual density is greater than 0.1 g·cm-3.

圖14 含噪聲反演結果圖(a)、(b)、(c)及(d)分別為y=1.3～1.4 km、x=1.9～2.0 km、x=2.5～2.6 km及z=0.35～0.4 km的反演結果切塊圖，白色虛線框為模型體的邊界位置，(e)、(f)為剩余密度大于0.1 g·cm-3的三維透視圖.Fig.14 Inversion result in the presence of noise(a), (b), (c) and (d) are block diagrams of inversion results at y=1.3～1.4 km、x=1.9～2.0 km、x=2.5～2.6 km and z=0.35～0.4 km respectively, the white dotted box indicates the boundary position of the model, (e) and (f) are the three-dimensional perspectives, where residual density is greater than 0.1 g·cm-3.

圖15 重力理論異常與預測結果正演的殘差(a)無噪聲情況; (b)含噪聲情況.Fig.15 The residuals of theoretical gravity anomaly and forward value of prediction result(a) Free noise situation; (b) Noise situation.

4 實際數據反演

采用本文的方法對美國路易斯安那州西南部Vinton鹽丘的地面重力異常和航空重力梯度異常數據進行反演.鹽丘構造的研究對于尋找油氣區帶及封閉圈具有非常重要的作用(賈承造等, 2003).Vinton鹽丘是一個位于新生代和中生帶晚期頁巖砂巖中的淺層刺穿鹽丘，是由一定規模的鹽核和其頂部的鹽丘帽組成，鹽丘帽是由石灰巖、石膏、硬石膏和黏土組成(Coker et al., 2007).有較多鉆孔穿過該鹽丘，鹽丘的成分、物性大小、空間分布狀態資料較為詳細(Ennen, 2012)，因此，該地區的重力數據也被諸多學者用于檢驗重力數據處理方法的效果(Yuan et al., 2013;秦朋波和黃大年, 2016).

鹽丘帽埋深大約位于地下200～650 m之間，密度為2.75 g·cm-3，周圍圍巖密度約為2.2 g·cm-3.因此，選擇2.2 g·cm-3的密度對85 m高空處的FTG數據進行地形改正，以及對高頻噪聲進行抑制，其結果如圖16a、b、c、d、e、f所示，同時采用插值分割法(徐世浙等, 2009)對地面重力異常數據進行局部重力異常的提取，其結果如圖16g所示，以上預處理結果可以認為是由0.55 g·cm-3的剩余密度所產生.圖中觀測區大小為4.0 km×4.0 km，數據點距為100 m×100 m，共41×41個點，將地下半空間劃分為40×40×20個長方體模型單元，其大小為100 m×100 m×50 m.

根據已知情況，并利用本文所提的GraInvNet網絡對圖16所示的數據進行預測，該過程對反演密度進行了約束，其范圍為：-0.2～0.55 g·cm-3.最終反演結果如圖17所示，從圖17a、b可以看出巖丘帽中心深度大約為500 m，深度范圍在300～700 m之間，東西向長約為1.2 km，南北向長約為1.1 km，這一結果與已知的地質信息較為吻合(Thompson and Eichelberger, 1928).從三維透視圖(圖17c、d)可以看出高密度的物性分布更為“聚焦”，近似為不規則橢球體，低密度至高密度過渡帶范圍較小，約為100 m間隔，表明密度較大值占反演結果的主要分布范圍.

5 結論與討論

本文提出一種基于DL的重力異常及重力梯度異常聯合反演的新方法，該方法對應的DL網絡結構稱為GraInvNet.本文首先提出了一種快速正演計算方法對GraInvNet網絡所需要的樣本數據集進行高效構建，從而提高DL反演的實用性；然后，將樣本數據輸入到GraInvNet網絡，以此來學習從正演數據到模型的復雜映射關系.在訓練階段，根據損失函數值與驗證值的變化來優化網絡結構、正則化方案等；最后，通過已訓練完成的GraInvNet網絡進行反演預測，此階段可以達到高效的計算速度，通常幾十毫秒就可完成.

組合模型一為樣本數據集相同模型的組合，重力異常的反演結果能夠“聚焦”于真實模型，且采用單一塊體模型所構建的樣本數據集對GraInvNet網絡進行訓練，該網絡可對類似塊體的組合模型進行預測，表明基于DL的重力異常反演所具有的可行性，以及本文GraInvNet網絡所具有的較強泛化性.隨著多維度重力梯度異常信息的加入，GraInvNet網絡的反演結果精度顯著提升，含噪數據試驗表明GraInvNet網絡具有較強的抗噪聲性能.組合模型二為樣本數據集的類型模型，GraInvNet網絡的反演結果與理論模型高度一致，表明本文反演方法具有較強的泛化性能和魯棒性.將文中的方法應用于實測的航空重力梯度異常和地面重力異常的聯合反演，其結果與已知的地質情況吻合，表明本文的反演方法可以應用于實際數據.

圖16 Vinton巖丘實測重力數據(a)、(b)、(c)、(d)、(e)、(f)及(g)分別對應Uxx、Uxy、Uxz、ΔUyy、Uyz、Uzz與Δg分量.Fig.16 The measured gravity data of Vinton salt dome(a), (b), (c), (d), (e), (f) and (g) correspond to Uxx, Uxy, Uxz, ΔUyy, Uyz, Uzz and Δg components respectively.

圖17 Vinton巖丘實測數據反演結果(a)與(b)分別為南北向1.6～1.7 km與東西向1.9～2.0 km的反演結果切片圖；(c)為反演密度分別為0.45、0.2 g·cm-3的包絡曲面圖；(d)為反演結果在深度方向上切片圖，切片深度為250 m變化至650 m，間隔100 m.Fig.17 Inversion results of the measured data of Vinton salt dome(a) and (b) are slice maps of inversion results in north-south direction of 1.6～1.7 km and east-west direction of 1.9～2.0 km respectively; (c) is the envelope surface with inversion density of 0.45 and 0.2 g·cm-3 respectively; (d) is slice map of inversion results in depth direction with slice depth from 250 m to 650 m and interval of 100 m.

基于DL的重力反演實質上是利用FCN方法對地下半空間的剩余密度體的位置及物性大小進行分類預測，其“端到端”的網絡結構也可以應用到其他地球物理反演中，如磁及磁梯度反演、電磁反演及面波波速反演等.另一方面，由于地球物理反演的復雜性，DL反演方法需具有的普適性在于如何構建一個足量的、多樣性的、代表性的高質量樣本數據集，然而，這勢必會引起網絡復雜度與計算量增加.因此，在下一步工作中，我們將進一步增加樣本的多樣性來提升復雜地質模型的反演效果，同時對網絡結構與節點超參數進一步優化來提高訓練效率和反演精度，從而使得基于DL學習的重力反演方法更加普適化與智能化.

致謝感謝Bell Geospace公司提供Vinton巖丘實測重力數據.非常感謝匿名審稿專家對論文提出的寶貴修改建議.