基于深度學(xué)習(xí)的重力異常與重力梯度異常聯(lián)合反演

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

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

0 引言

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

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

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

基于此，本文提出一種基于DL的重力及重力梯度數(shù)據(jù)聯(lián)合反演的方法.反演流程包括三個(gè)階段：樣本數(shù)據(jù)集生成，網(wǎng)絡(luò)模型訓(xùn)練和利用網(wǎng)絡(luò)模型進(jìn)行預(yù)測.前兩個(gè)階段計(jì)算耗時(shí)較長，因?yàn)楸仨氃谧銐虼蟮挠?xùn)練集上訓(xùn)練網(wǎng)絡(luò)模型，才能使網(wǎng)絡(luò)模型具有泛化新數(shù)據(jù)的能力，但前兩個(gè)階段僅執(zhí)行一次.反演預(yù)測較為迅速，通常在幾十毫秒內(nèi)完成.針對(duì)反演第一階段，本文首先提出了一種重力及重力梯度空間域快速正演方法，以此來快速構(gòu)建樣本數(shù)據(jù)集；對(duì)于DL學(xué)習(xí)的網(wǎng)絡(luò)結(jié)構(gòu)，本文基于U-net網(wǎng)絡(luò)結(jié)構(gòu)(Russakovsky et al., 2015)建立重力數(shù)據(jù)(輸入)到剩余密度模型(輸出)的非線性映射網(wǎng)絡(luò)結(jié)構(gòu)，即三維重力場反演的“端到端”網(wǎng)絡(luò)結(jié)構(gòu)，本文稱為GraInvNet (Gravity Inverse Net)網(wǎng)絡(luò)結(jié)構(gòu)；然后利用樣本數(shù)據(jù)集對(duì)網(wǎng)絡(luò)進(jìn)行監(jiān)督學(xué)習(xí)(訓(xùn)練)，最后利用GraInvNet網(wǎng)絡(luò)對(duì)理論數(shù)據(jù)和實(shí)測重力及重力梯度數(shù)據(jù)進(jìn)行反演，以此來檢驗(yàn)文中方法的可行性和有效性.

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

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

1.1 正演計(jì)算表達(dá)式

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

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

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

(1)

(2)

(3)

(4)

(5)

(6)

(7)

1.2 等效存儲(chǔ)網(wǎng)格點(diǎn)幾何格架技術(shù)

姚長利等(2003)提出了幾何格架技術(shù)進(jìn)行了重力異常的空間域快速正演計(jì)算，取得較好的效果.本文在此基礎(chǔ)上提出一種網(wǎng)格點(diǎn)幾何格架技術(shù)，此技術(shù)能夠?qū)σ延械恼菟惴ㄟM(jìn)行加速.以重力異常的正演計(jì)算為例進(jìn)行說明，公式(1)可以寫成如下表達(dá)式：

(8)

其中，

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

(9)

與物性無關(guān)，本文稱為網(wǎng)格點(diǎn)幾何格架函數(shù).網(wǎng)格點(diǎn)幾何格架函數(shù)與文獻(xiàn)(姚長利等, 2003)類似的幾何格架函數(shù)：

S(xi,yj,zk)=

(10)

存在如下關(guān)系：

(11)

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

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

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

1.3 計(jì)算步驟

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

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

步驟2：分別計(jì)算地下半空間內(nèi)所有網(wǎng)格點(diǎn)對(duì)長方形觀測區(qū)4個(gè)頂點(diǎn)的重力異?；蛑亓μ荻犬惓５木W(wǎng)格點(diǎn)幾何格架值(分別記為Q1，Q2，Q3及Q4.圖3所示)，并存儲(chǔ)待后續(xù)計(jì)算調(diào)用；

步驟3：根據(jù)計(jì)算點(diǎn)與任意單一長方體模型單元8個(gè)頂點(diǎn)之間的相對(duì)位置關(guān)系，利用平移等效性調(diào)取其相對(duì)應(yīng)的網(wǎng)格點(diǎn)幾何格架值，再依據(jù)公式 (8)計(jì)算該長方體模型單元對(duì)計(jì)算點(diǎn)的重力異常或重力梯度異常；

步驟4：利用步驟3計(jì)算地下半空間每一個(gè)長方體模型單元對(duì)某一計(jì)算點(diǎn)的重力異?；蛑亓μ荻犬惓?，最后求和得到該計(jì)算點(diǎn)總的重力異?；蛑亓μ荻犬惓?；

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

圖3 計(jì)算點(diǎn)與網(wǎng)格點(diǎn)一一對(duì)應(yīng)關(guān)系Fig.3 The corresponding relation between calculation points and grid points

2 深度學(xué)習(xí)反演

2.1 方法

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

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

圖4 FCN端對(duì)端示意圖(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.

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

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

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

(12)

2.2 網(wǎng)絡(luò)結(jié)構(gòu)

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

2.3 正則化

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

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

2.4 損失函數(shù)

損失函數(shù)是DL網(wǎng)絡(luò)訓(xùn)練過程的指揮棒，其通過正向傳播的輸出結(jié)果與樣本標(biāo)簽值之間的誤差反向傳播來修改網(wǎng)絡(luò)單元節(jié)點(diǎn)權(quán)值及其他參量，本文指通過預(yù)測地下半空間長方體模型密度與已知模型密度的誤差進(jìn)行調(diào)參.常用損失誤差為L1、L2范數(shù)，其表達(dá)式分別為

圖5 GraInvNet網(wǎng)絡(luò)結(jié)構(gòu)由3個(gè)下采樣層和3個(gè)上采樣層組成，下采樣和左側(cè)卷積層對(duì)應(yīng)編碼階段，上采樣和右側(cè)卷積層對(duì)應(yīng)解碼階段，圖中黑色數(shù)字表示特征圖的通道數(shù)，Conv為卷積，BN為批歸一化，ReLU為激活函數(shù)，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網(wǎng)絡(luò)模型訓(xùn)練的實(shí)質(zhì)就是損失函數(shù)數(shù)學(xué)期望極小化的過程.傳統(tǒng)地球物理反演的實(shí)質(zhì)是目標(biāo)函數(shù)極小化的過程.其實(shí)兩者的二范數(shù)表達(dá)式類似，區(qū)別在于傳統(tǒng)反演目標(biāo)函數(shù)是預(yù)測數(shù)據(jù)的正演響應(yīng)與實(shí)測數(shù)據(jù)的二范數(shù)或其正則化改進(jìn)，且該目標(biāo)函數(shù)極小化過程的計(jì)算成本極為昂貴，而DL損失函數(shù)極小化過程效率較高.本文采用歸一化數(shù)據(jù)的L2范數(shù)作為損失函數(shù).

2.5 學(xué)習(xí)算法的選擇

強(qiáng)大的學(xué)習(xí)算法是提高深度學(xué)習(xí)反演性能必不可少的步驟，通過GraInvNet網(wǎng)絡(luò)學(xué)習(xí)算法一步步的迭代使損失函數(shù)極小化，達(dá)到全局的最優(yōu)值.

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

2.6 網(wǎng)絡(luò)規(guī)格及調(diào)參

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

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

(14)

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

(15)

在調(diào)參方面需要考慮到訓(xùn)練輪數(shù)、樣本批量數(shù)以及Dorpout的比例等因素對(duì)模型性能的影響.以上因素需兼顧計(jì)算效率和精度，從而優(yōu)化GraInvNet網(wǎng)絡(luò)結(jié)構(gòu)，最終的訓(xùn)練完成的GraInvNet具備計(jì)算的準(zhǔn)確性，同時(shí)也需對(duì)新樣本具備一定的泛化性.

3 理論模型試驗(yàn)

3.1 樣本數(shù)據(jù)集的構(gòu)建及訓(xùn)練

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

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

訓(xùn)練輪數(shù)為100，每次訓(xùn)練后進(jìn)行一次驗(yàn)證，圖6所示為本文所提GraInvNet網(wǎng)絡(luò)的訓(xùn)練誤差與驗(yàn)證誤差，可以看出，模型訓(xùn)練良好，未出現(xiàn)過欠擬合.

圖6 GraInvNet訓(xùn)練誤差及驗(yàn)證誤差Fig.6 Training error and validation error of GraInvNet

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

3.2 組合模型一

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

表1 長方體模型參數(shù)Table 1 Cuboid model parameters

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

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

圖8 反演結(jié)果三維透視圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果，(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果.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 反演結(jié)果切塊圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果；(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果；切塊線首尾坐標(biāo)分別為(0 km，0 km)和(4 km，4 km)，白色虛線框?yàn)槟Ｐ腕w的邊界位置.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反演結(jié)果切塊圖(a)、(b)及(c)分別為無噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果；(d)、(e)及(f)分別為含噪聲情況下Δg、Δg & Uxz & Uyz & Uzz與Δg &全部分量反演結(jié)果；白色虛線框?yàn)槟Ｐ腕w的邊界位置.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.

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

3.3 組合模型二

組合模型二為四層臺(tái)階與長方體的組合模型，如圖12所示，其中四層臺(tái)階模型的每層臺(tái)階大小相同都為1.4 km×0.6 km×0.15 km，每層臺(tái)階的中心坐標(biāo)分別為(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，中心點(diǎn)坐標(biāo)為(2.0 km,3.05 km,0.35 km)，四層臺(tái)階與長方體剩余密度都為1.0 g·cm-3.組合模型二相比組合模型一的變化是加入新類型模型——臺(tái)階，以及樣本數(shù)據(jù)集中未包含的長方體模型.利用本文的GraInvNet網(wǎng)絡(luò)對(duì)組合模型二的重力異常及全部重力梯度異常進(jìn)行聯(lián)合反演，其結(jié)果如圖13所示，可以看出臺(tái)階模型的反演結(jié)果與真實(shí)模型高度吻合，其幅值在深部處略有差別(圖13c),長方體模型在深部的幅值也略有差別(圖13e及f)，其余部分也與理論模型高度吻合.表明GraInvNet網(wǎng)絡(luò)可以對(duì)與樣本模型相似組合模型的異常體進(jìn)行精確定位，同時(shí)也可以準(zhǔn)確地預(yù)測出模型值的大小，表明其反演能力很好.為了更加直觀地了解反演結(jié)果中異常體的位置、形狀和物性值，圖13a、b、c與d 分別給出不同方向的反演切塊圖，可以看出反演結(jié)果與真實(shí)模型的邊界及形狀吻合，并且反演結(jié)果的幅值大小也與已知模型接近.

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

表2 理論異常最大最小值、反演結(jié)果擬合殘差的標(biāo)準(zhǔn)差及隨機(jī)噪聲標(biāo)準(zhǔn)差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

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

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

表3 理論異常最大最小值、反演結(jié)果擬合殘差的標(biāo)準(zhǔn)差及隨機(jī)噪聲標(biāo)準(zhǔn)差Table 3 The maximum and minimum value of theoretical anomaly, the standard deviation of the fitting residuals of inversion result and random noise

圖13 無噪聲反演結(jié)果(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的反演結(jié)果切塊圖，白色虛線框?yàn)槟Ｐ腕w的邊界位置，(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 含噪聲反演結(jié)果圖(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的反演結(jié)果切塊圖，白色虛線框?yàn)槟Ｐ腕w的邊界位置，(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 重力理論異常與預(yù)測結(jié)果正演的殘差(a)無噪聲情況; (b)含噪聲情況.Fig.15 The residuals of theoretical gravity anomaly and forward value of prediction result(a) Free noise situation; (b) Noise situation.

4 實(shí)際數(shù)據(jù)反演

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

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

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

5 結(jié)論與討論

本文提出一種基于DL的重力異常及重力梯度異常聯(lián)合反演的新方法，該方法對(duì)應(yīng)的DL網(wǎng)絡(luò)結(jié)構(gòu)稱為GraInvNet.本文首先提出了一種快速正演計(jì)算方法對(duì)GraInvNet網(wǎng)絡(luò)所需要的樣本數(shù)據(jù)集進(jìn)行高效構(gòu)建，從而提高DL反演的實(shí)用性；然后，將樣本數(shù)據(jù)輸入到GraInvNet網(wǎng)絡(luò)，以此來學(xué)習(xí)從正演數(shù)據(jù)到模型的復(fù)雜映射關(guān)系.在訓(xùn)練階段，根據(jù)損失函數(shù)值與驗(yàn)證值的變化來優(yōu)化網(wǎng)絡(luò)結(jié)構(gòu)、正則化方案等；最后，通過已訓(xùn)練完成的GraInvNet網(wǎng)絡(luò)進(jìn)行反演預(yù)測，此階段可以達(dá)到高效的計(jì)算速度，通常幾十毫秒就可完成.

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

圖16 Vinton巖丘實(shí)測重力數(shù)據(jù)(a)、(b)、(c)、(d)、(e)、(f)及(g)分別對(duì)應(yīng)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巖丘實(shí)測數(shù)據(jù)反演結(jié)果(a)與(b)分別為南北向1.6～1.7 km與東西向1.9～2.0 km的反演結(jié)果切片圖；(c)為反演密度分別為0.45、0.2 g·cm-3的包絡(luò)曲面圖；(d)為反演結(jié)果在深度方向上切片圖，切片深度為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的重力反演實(shí)質(zhì)上是利用FCN方法對(duì)地下半空間的剩余密度體的位置及物性大小進(jìn)行分類預(yù)測，其“端到端”的網(wǎng)絡(luò)結(jié)構(gòu)也可以應(yīng)用到其他地球物理反演中，如磁及磁梯度反演、電磁反演及面波波速反演等.另一方面，由于地球物理反演的復(fù)雜性，DL反演方法需具有的普適性在于如何構(gòu)建一個(gè)足量的、多樣性的、代表性的高質(zhì)量樣本數(shù)據(jù)集，然而，這勢必會(huì)引起網(wǎng)絡(luò)復(fù)雜度與計(jì)算量增加.因此，在下一步工作中，我們將進(jìn)一步增加樣本的多樣性來提升復(fù)雜地質(zhì)模型的反演效果，同時(shí)對(duì)網(wǎng)絡(luò)結(jié)構(gòu)與節(jié)點(diǎn)超參數(shù)進(jìn)一步優(yōu)化來提高訓(xùn)練效率和反演精度，從而使得基于DL學(xué)習(xí)的重力反演方法更加普適化與智能化.

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