基于深度學(xué)習(xí)的2D-3D地下密度異常體反演方法

朱雨菲, 王武星, 曹建玲, 郝梅煊

中國(guó)地震局地震預(yù)測(cè)研究所, 北京 100036

0 引言

重力位場(chǎng)數(shù)據(jù)是地球內(nèi)部密度結(jié)構(gòu)的直接反映,重力勘測(cè)具有快速、經(jīng)濟(jì)、范圍大等優(yōu)勢(shì),是進(jìn)行構(gòu)造劃分和異常圈定快速而有效的手段(馬國(guó)慶, 2013).通過分離出與研究對(duì)象對(duì)應(yīng)的重力異常區(qū)域,根據(jù)地表重力異常范圍、數(shù)值大小等進(jìn)行反演,就能獲得密度異常體的地下分布形態(tài).因此,重力異常反演是研究地下物質(zhì)性質(zhì)結(jié)構(gòu)的重要手段之一.地下密度異常體常指示著地下礦產(chǎn)資源的存在,對(duì)油氣勘探(Constantino et al., 2019; ?zdemir et al., 2022)、成礦預(yù)測(cè)(袁峰等, 2014)等均有重要意義.

當(dāng)?shù)叵麓嬖谂c圍巖密度不一致的異常埋藏體時(shí),會(huì)引起地表局部重力異常(Singh and Biswas, 2016)、重力梯度異常(Hou et al., 2016),可能引起磁異常(王昊, 2020),研究人員通常利用這些地表異常進(jìn)行反演,得到目標(biāo)區(qū)域地下異常體的埋深、位置、范圍等信息.傳統(tǒng)的反演方法大多基于反演理論,引入約束條件,在最小二乘意義上使目標(biāo)函數(shù)達(dá)到最小值而獲得最優(yōu)解,包括希爾伯特變換(Hinojosa and Mickus, 2002)、歐拉反褶積(Thompson, 1982)等.然而,傳統(tǒng)方法對(duì)初始條件敏感,反演重力異常時(shí)參數(shù)眾多,反演矩陣維數(shù)大,計(jì)算過程復(fù)雜(郭良輝等, 2009),難以取得理想的結(jié)果.

早些年,機(jī)器學(xué)習(xí)在地球物理領(lǐng)域的應(yīng)用包括模擬退火算法(張霖斌, 1997)、支持向量機(jī)(張克, 2011)、反向傳播神經(jīng)網(wǎng)絡(luò)(郭文斌等, 2012)、徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)(耿美霞和楊慶節(jié), 2013)等,受限于計(jì)算機(jī)性能、樣本數(shù)量及模型復(fù)雜度等,這些方法的計(jì)算結(jié)果往往與真實(shí)分布存在較明顯差異,并不顯著優(yōu)于傳統(tǒng)方法.隨著計(jì)算機(jī)人工智能(Artificial Intelligence,AI)技術(shù)的發(fā)展,一些學(xué)者開始利用深度學(xué)習(xí)方法解決地球物理問題,當(dāng)前主要聚焦于地震震相拾取(Zhu and Beroza, 2018; Zhou et al. 2021)、地磁反演(王丹鶴, 2021; Xiong et al., 2020)、斷層識(shí)別(常德寬等, 2021)等領(lǐng)域.相對(duì)傳統(tǒng)方法,經(jīng)歷完整的訓(xùn)練過程后,深度學(xué)習(xí)能在極短的時(shí)間內(nèi)獲取更高精度的結(jié)果.

在深度學(xué)習(xí)中,卷積神經(jīng)網(wǎng)絡(luò)常用于提取圖像信息,而重磁及其梯度數(shù)據(jù)與影像數(shù)據(jù)結(jié)構(gòu)相似,可以網(wǎng)格化表示.因此,近些年國(guó)內(nèi)外學(xué)者開始將重磁及其梯度數(shù)據(jù)作為輸入,利用深度卷積神經(jīng)網(wǎng)絡(luò)反演地下結(jié)構(gòu).例如He等(2021)將重力數(shù)據(jù)作為輸入,利用卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks,CNN)反演沉積物基底深度,捕捉復(fù)雜地質(zhì)模型的沉積物-基底界面的詳細(xì)特征;Hu等(2021)利用小批量正演磁異常數(shù)據(jù),結(jié)合深度神經(jīng)網(wǎng)絡(luò),恢復(fù)磁性礦體的物理性質(zhì)及分布,獲得了分布更集中、邊界更清晰的反演結(jié)果.在利用重力異常數(shù)據(jù)反演地下密度異常體方面,王逸宸等(2020)將重力異常作為輸入,利用卷積神經(jīng)網(wǎng)絡(luò)LeNet-5,反演了澳大利亞Kauring地區(qū)異常體密度、線度信息.Yang等(2022)利用理論合成無噪重力異常數(shù)據(jù)作為輸入,對(duì)比了傳統(tǒng)最小二乘方法、全卷積神經(jīng)網(wǎng)絡(luò)(Fully Convolutional Networks,FCN)、卷積神經(jīng)網(wǎng)絡(luò)(CNN)的反演效果,并用實(shí)測(cè)數(shù)據(jù)進(jìn)行測(cè)試,結(jié)果表明了CNN反演地下密度異常體的優(yōu)勢(shì).張志厚等(2021)結(jié)合重力異常與重力梯度異常數(shù)據(jù),利用大量規(guī)則樣本訓(xùn)練深度學(xué)習(xí)模型,得到了邊界十分清晰的地下密度異常體三維反演結(jié)果.Huang等(2021)用大量高復(fù)雜度的樣本訓(xùn)練深度學(xué)習(xí)模型,使其具有反演復(fù)雜樣本清晰邊界的能力.這些研究表明深度學(xué)習(xí)能通過地表異常數(shù)據(jù),有效提取密度異常體的三維信息,反演其地下三維結(jié)構(gòu).

地表的重力異常和重力梯度異常均為二維數(shù)據(jù),而地下異常體分布為三維密度結(jié)構(gòu),現(xiàn)有方法多利用二維深度學(xué)習(xí)網(wǎng)絡(luò)來反解這種三維信息(張志厚等,2021; Hu et al., 2021; Huang et al., 2021),且理論反演時(shí)多把樣本的密度異常設(shè)置為固定值(秦朋波和黃大年,2016; 鄭玉君等, 2020).通常,深度學(xué)習(xí)利用二維網(wǎng)絡(luò)提取二維數(shù)據(jù)的抽象信息(Shelhamer et al., 2017),利用三維網(wǎng)絡(luò)處理三維數(shù)據(jù)(Tran et al., 2016).因此,本文在模型結(jié)構(gòu)上做了新的嘗試,提出2D-3D InvNet深度學(xué)習(xí)網(wǎng)絡(luò),采用二維至三維的編-解碼器結(jié)構(gòu),在編碼階段利用二維網(wǎng)絡(luò)結(jié)構(gòu)提取網(wǎng)格數(shù)據(jù)的抽象信息,解碼階段利用三維網(wǎng)絡(luò)結(jié)構(gòu)還原異常體在地下半空間的位置分布及剩余密度.為了使網(wǎng)絡(luò)模型能反演密度不同的地下異常體,本研究不再固定訓(xùn)練數(shù)據(jù)的剩余密度,而是利用隨機(jī)函數(shù),控制剩余密度在一定范圍內(nèi)波動(dòng).實(shí)際實(shí)驗(yàn)發(fā)現(xiàn),若以常用的均方誤差(Mean Square Error,MSE)作為損失函數(shù),無法在數(shù)值上準(zhǔn)確擬合密度多變的地下異常體.本文將地下半空間分為圍巖分布的背景區(qū)與異常體分布的目標(biāo)區(qū),改進(jìn)了深度學(xué)習(xí)模型的損失函數(shù),將目標(biāo)區(qū)與背景區(qū)的誤差值進(jìn)行分離并重新賦予權(quán)重,使最終結(jié)果更關(guān)注于局部區(qū)域密度異常,從而更好的恢復(fù)地下異常體的密度信息.

1 重力異常及重力梯度異常樣本正演

訓(xùn)練深度學(xué)習(xí)網(wǎng)絡(luò)需要大量的樣本,然而目前缺乏大量已知地下三維物質(zhì)分布的樣本,因此,本文采用姚長(zhǎng)利等(2003)提出的幾何格架存儲(chǔ)技術(shù),結(jié)合正演計(jì)算公式模擬地下異常體在地表產(chǎn)生的異常響應(yīng).同時(shí),為每個(gè)正演的樣本添加高斯噪聲,生成大量已知地下密度分布的含噪樣本.最后將樣本按比例隨機(jī)分配為訓(xùn)練樣本和驗(yàn)證樣本,形成訓(xùn)練集和驗(yàn)證集.

1.1 正演計(jì)算方法

姚長(zhǎng)利等(2003)將地下劃分為規(guī)則長(zhǎng)方體網(wǎng)格,利用地面點(diǎn)和網(wǎng)格之間的相對(duì)位置關(guān)系計(jì)算地下密度異常體在地表引起的響應(yīng),地表點(diǎn)(x,y,0)和地下網(wǎng)格(ξi,ηj,ζk)(i,j,k為整數(shù),取1或2)位置關(guān)系如圖1所示.

圖1 地下網(wǎng)格示意圖(x,y,0)為地表觀測(cè)點(diǎn)坐標(biāo),(ξ1,η1,ζ1)和(ξ2,η2,ζ2)分別為長(zhǎng)方體單元網(wǎng)格距原點(diǎn)最近點(diǎn)與最遠(yuǎn)點(diǎn).

由于坐標(biāo)系統(tǒng)、參數(shù)定義方式等的不同,網(wǎng)格化的正演數(shù)值模擬有多種形式.綜合眾多正演模擬的文獻(xiàn)(Li and Chouteau, 1998; Pilkington, 2014; 陳召曦等, 2012; 陳閆等, 2014),本文采用圖1所示坐標(biāo)系(原點(diǎn)在地表),及與該坐標(biāo)系定義相同的正演公式(Li and Chouteau, 1998)進(jìn)行正演計(jì)算.

(1a)

(1b)

(2)

為了檢驗(yàn)正演計(jì)算的正確性,我們劃定了250 m×250 m×60 m的地下半空間,設(shè)置50×50的地面測(cè)網(wǎng),按測(cè)網(wǎng)大小將地下半空間劃分為網(wǎng)格,單元網(wǎng)格大小為5 m×5 m×3 m.假設(shè)在網(wǎng)格坐標(biāo)(20,20,13)處存在由10×2×2個(gè)單元網(wǎng)格組成、剩余密度為1.5 g·cm-3的異常體.利用公式(1)和公式(2) 計(jì)算其在地表產(chǎn)生的重力異常和重力梯度異常,獲得的正演結(jié)果如圖2所示,通過對(duì)比分析(舒晴等,2015; Rogers et al., 2009)可見是正確的.

圖2 重力異常及重力梯度異常正演數(shù)值模擬結(jié)果其中,Δg為重力異常,Uxx,Uxy,Uxz,Uyy,Uyz,Uzz為重力梯度異常的六個(gè)分量.

1.2 樣本正演及樣本集設(shè)置

深度學(xué)習(xí)一般將樣本劃分為訓(xùn)練集、驗(yàn)證集和測(cè)試集,訓(xùn)練集是訓(xùn)練樣本的集合,直接參與網(wǎng)絡(luò)參數(shù)的調(diào)整.驗(yàn)證集是驗(yàn)證樣本的集合,用于評(píng)估模型在驗(yàn)證樣本上的效果及參與模型超參數(shù)的選取.測(cè)試集只在模型訓(xùn)練完成后使用,不參與訓(xùn)練過程,僅用于模型訓(xùn)練結(jié)果的評(píng)價(jià).設(shè)地下半空間圍巖區(qū)域的剩余密度為0,每個(gè)地下半空間樣本含1～3個(gè)異常塊體,每個(gè)異常塊體為邊長(zhǎng)占3～12個(gè)網(wǎng)格單元的長(zhǎng)方體或階梯型異常體,剩余密度取1～2.5 g·cm-3之間的隨機(jī)浮點(diǎn)值.本文的訓(xùn)練集含樣本30000條,驗(yàn)證集含樣本6000條.

為評(píng)估網(wǎng)絡(luò)模型在不同復(fù)雜度樣本上的表現(xiàn),參考他人(陳閆等, 2014; 張志厚等, 2021; Yang et al., 2022)的研究,我們定義了十余條形態(tài)各異的理論測(cè)試樣本,包括正方體、長(zhǎng)方體、階梯型樣本以及這些基礎(chǔ)元素的組合.通過繪制剖切面和三維立體圖,可視化展示部分樣本來分析評(píng)估反演效果.表1給出了參與反演結(jié)果可視化展示的典型樣本.每個(gè)樣本含七層信息,為1個(gè)重力異常和6個(gè)重力梯度異常分量,樣本標(biāo)簽為引起當(dāng)前地表響應(yīng)的地下密度異常(網(wǎng)格)分布.為模擬真實(shí)測(cè)量數(shù)據(jù),提高網(wǎng)絡(luò)模型的抗噪能力,每個(gè)樣本都添加了標(biāo)準(zhǔn)差為樣本數(shù)據(jù)最大值10%的高斯噪聲(圖3).

表1 測(cè)試樣本

圖3 樣本加噪處理示意圖(a) 地下半空間異常體位置示意圖; (b)、(c) 分別為正演生成的樣本及加噪后的樣本,其中堆疊的圖片從前到后分別為Uxx,Uxy,Uxz,Uyy,Uyz,Uzz,Δg.

2 2D-3D InvNet模型訓(xùn)練

2.1 深度學(xué)習(xí)反演算法

通常來說,反演問題都是欠定問題,其解具有不唯一性(Tikhonov and Arsenin, 1977).傳統(tǒng)反演手段中,人們選擇添加約束條件(Portniaguine and Zhdanov, 1999; Zhdanov, 2009)來抑制這種不唯一性,從而得到和真實(shí)情況相符合的反演結(jié)果.而深度學(xué)習(xí)利用一種端到端的結(jié)構(gòu),由數(shù)據(jù)驅(qū)動(dòng)模型訓(xùn)練,不需要添加約束條件,就能根據(jù)以往學(xué)習(xí)經(jīng)驗(yàn),自適應(yīng)地得出符合樣本分布的結(jié)果.

重力異常及重力梯度異常數(shù)據(jù)是結(jié)構(gòu)化的網(wǎng)格數(shù)據(jù),相鄰格網(wǎng)之間具有一定的關(guān)聯(lián),單獨(dú)考慮強(qiáng)度信息不足以實(shí)現(xiàn)精確的反演,還需要考慮相鄰網(wǎng)格的輔助信息.為了同時(shí)利用強(qiáng)度與空間信息,深度學(xué)習(xí)采用卷積神經(jīng)網(wǎng)絡(luò),如FCN(Shelhamer et al., 2017)、InceptionNet(Szegedy et al., 2016).這些工作多是在二維、三維數(shù)據(jù)的基礎(chǔ)上,獲得同維度或降維輸出.對(duì)于密度異常體反演問題,需要將二維輸入反解為三維輸出.為了更好地提取地表數(shù)據(jù)信息來恢復(fù)地下三維形態(tài),我們采用二維卷積神經(jīng)網(wǎng)絡(luò)對(duì)地面數(shù)據(jù)進(jìn)行編碼,獲得數(shù)據(jù)的抽象特征;利用三維卷積神經(jīng)網(wǎng)絡(luò)對(duì)提取的特征進(jìn)行解碼,恢復(fù)異常體的幾何位置及密度信息.

2.2 網(wǎng)絡(luò)結(jié)構(gòu)及參數(shù)設(shè)置

為了提取地表信息反演地下三維結(jié)構(gòu),本文構(gòu)建了2D-3D InvNet網(wǎng)絡(luò),將U形網(wǎng)絡(luò)結(jié)構(gòu)UNet(Ronneberger et al., 2015)和3D UNet (?i?ek et al., 2016)進(jìn)行有效結(jié)合.UNet是一種二維語義分割網(wǎng)絡(luò),廣泛應(yīng)用于遙感影像分割(Ye et al., 2019)和醫(yī)學(xué)影像分割(Xu et al., 2021)等領(lǐng)域,特點(diǎn)是能完成逐像素分類.由于三維立體數(shù)據(jù)十分常見,對(duì)數(shù)據(jù)進(jìn)行切片處理的方式冗余且低效,3D UNet將UNet結(jié)構(gòu)擴(kuò)展至三維,直接對(duì)三維數(shù)據(jù)進(jìn)行處理及分割,得到三維結(jié)果(?i?ek et al., 2016).地下密度異常體反演由二維輸入獲得三維結(jié)果,2D-3D InvNet網(wǎng)絡(luò)采用U形框架,左半邊使用與UNet相同的結(jié)構(gòu)及2D卷積,對(duì)二維輸入數(shù)據(jù)進(jìn)行編碼,獲得圖像的高階抽象信息;右半邊采用與3D UNet相同的結(jié)構(gòu)及3D卷積,對(duì)抽象信息進(jìn)行解碼,獲取地下三維反演結(jié)果.編碼過程含3個(gè)卷積塊,為填補(bǔ)底層信息以提高結(jié)果精度,將編碼過程每個(gè)卷積塊中的最后一層運(yùn)算結(jié)果和同級(jí)別的解碼過程進(jìn)行跳躍連結(jié).在跳躍連接前,將同級(jí)別編碼結(jié)果擴(kuò)展至三維;即如圖4所示,同級(jí)別的編碼結(jié)果與解碼塊進(jìn)行跳躍連接前,要先在新的維度重復(fù)堆疊,使深度與解碼塊輸入數(shù)據(jù)一致,其中深度為數(shù)據(jù)的第三維度.

圖4 2D-3D InvNet網(wǎng)絡(luò)

除編碼器和解碼器的連接部分外,每個(gè)卷積塊各含2次卷積,下文除特別說明外,所有操作步長(zhǎng)均為1,卷積過程默認(rèn)含填充(padding).編碼過程中選用3×3的卷積核,解碼過程中選用3×3×3的卷積核.塊與塊之間分別由上、下采樣進(jìn)行連接,下采樣過程為2×2的最大池化,上采樣為卷積核大小2×2×2、步長(zhǎng)2的反卷積.編-解碼器的連接部分將數(shù)據(jù)形式擴(kuò)充至三維,深度為1.此外,網(wǎng)絡(luò)還在每次卷積操作后添加了批歸一化(batch norm,BN)和激活函數(shù)(Relu),具體設(shè)置見圖4.

2.2 數(shù)據(jù)標(biāo)準(zhǔn)化處理

正演生成的單個(gè)樣本,重力異常與重力梯度異常信息層的量綱不同,數(shù)量級(jí)上的差異可能導(dǎo)致數(shù)值更大的數(shù)據(jù)層在參數(shù)更新過程中占更高的比例.因此,在進(jìn)行運(yùn)算前,我們對(duì)數(shù)據(jù)做了Z-Score標(biāo)準(zhǔn)化處理.除了將所有數(shù)據(jù)各維度限制在同一數(shù)量級(jí)范圍內(nèi),Z-Score能處理數(shù)據(jù)最大和最小值未知或樣本存在離群值的情況,無需提前定義樣本最大和最小值.Z-Score標(biāo)準(zhǔn)化處理方程如下:

(3)

(4)

2.3 核密度估計(jì)

概率密度是一個(gè)統(tǒng)計(jì)學(xué)上的概念,其值等于一段區(qū)間內(nèi)事件發(fā)生的概率除以該區(qū)間的長(zhǎng)度.核密度估計(jì)是一種估計(jì)連續(xù)數(shù)據(jù)概率密度函數(shù)的非參數(shù)方法,能通過有限的樣本分布直方圖,推斷總體的數(shù)據(jù)分布(圖5a).核密度估計(jì)的表達(dá)公式如下(Botev et al., 2010):

通常,假設(shè)存在剩余密度為ρ的異常體,地下半空間中密度的分布應(yīng)該集中在0(背景值)和ρ(目標(biāo)值)處.如圖5b藍(lán)色曲線,剩余密度集中分布在0和1 g·cm-3處.由于核密度分布曲線的縱坐標(biāo)為概率密度,曲線與橫坐標(biāo)之間所包含的面積為此密度在整體密度分布中所占的比例,因此我們可以根據(jù)結(jié)果曲線與理論曲線之間ρ處的峰值差異、曲線包圍面積、聚焦程度等來判斷密度反演結(jié)果的優(yōu)劣.

圖 5 核密度估計(jì)示意圖(a) 密度分布直方圖與核密度曲線的關(guān)系,藍(lán)色條為密度柱狀圖,線為核密度曲線; (b) 反演結(jié)果與真實(shí)分布的關(guān)系,藍(lán)色為標(biāo)簽真值核密度分布,橙色為反演結(jié)果核密度分布. 圖中橫坐標(biāo)為樣本密度,縱坐標(biāo)為概率密度.

2.4 損失函數(shù)

在深度學(xué)習(xí)中,損失函數(shù)用來衡量真實(shí)值與反演結(jié)果之間的差異,訓(xùn)練過程中通過反向傳播計(jì)算梯度并完成參數(shù)更新.本文所涉及的地下密度異常體反演屬于回歸問題,一般用均方誤差作為損失函數(shù).均方誤差表達(dá)式如下:

(6)

WtdMSE=α·MSEtarget+(1-α)·MSEbackground,(7)

其中α為比例系數(shù),取值在0～1之間,MSEtarget和MSEbackground分別為目標(biāo)區(qū)和背景區(qū)均方誤差.取0.1～0.9之間的小數(shù),以不同比例系數(shù)的WtdMSE作為損失函數(shù)訓(xùn)練2D-3D InvNet,同時(shí)用MSE評(píng)估模型在驗(yàn)證樣本上的誤差.結(jié)果顯示隨著α的增大,整體測(cè)試誤差逐漸增大,目標(biāo)區(qū)誤差逐漸減小后在固定值周圍震蕩,拐點(diǎn)約在α=0.3處.為了使最終總體誤差與目標(biāo)區(qū)誤差均處于較低水平,綜合不同比例系數(shù)時(shí)模型的整體表現(xiàn),后續(xù)實(shí)驗(yàn)中設(shè)置α為0.3.

分別利用公式(6)和(7)作為損失函數(shù)訓(xùn)練2D-3D InvNet網(wǎng)絡(luò),對(duì)比MSE與WtdMSE的損失函數(shù)曲線 (圖6) ,發(fā)現(xiàn)利用MSE優(yōu)化模型參數(shù),盡管最終曲線收斂,但是目標(biāo)區(qū)誤差(圖6a)始終處于一個(gè)很高的數(shù)值,此時(shí)背景誤差主導(dǎo)了參數(shù)的更新過程,導(dǎo)致模型對(duì)目標(biāo)區(qū)異常體密度的預(yù)測(cè)并不準(zhǔn)確.反觀WtdMSE,能在損失函數(shù)值穩(wěn)定下降且最終收斂的同時(shí),顯著減少目標(biāo)區(qū)誤差(圖6b).

圖6 分別利用MSE和WtdMSE作為損失函數(shù)訓(xùn)練時(shí)的誤差曲線對(duì)比(a) 利用MSE作為損失函數(shù)的誤差曲線; (b) 利用WtdMSE作為損失函數(shù)的誤差曲線;藍(lán)色、橙色和綠色曲線分別為訓(xùn)練集、驗(yàn)證集和目標(biāo)區(qū)的誤差.

圖7和圖8展示了兩種損失函數(shù)訓(xùn)練的模型在具體樣本上的表現(xiàn).利用MSE優(yōu)化的網(wǎng)絡(luò)模型,對(duì)異常體密度的估計(jì)不如WtdMSE準(zhǔn)確.由圖7d和圖8d可見,損失曲線收斂后,采用MSE時(shí)對(duì)目標(biāo)區(qū)密度的估計(jì)不僅峰值上和真實(shí)值有一定的差異,范圍上(曲線包圍的面積)也有較大的差異,而采用WtdMSE時(shí)得到的結(jié)果更優(yōu).

圖7 利用MSE作為損失函數(shù)訓(xùn)練的結(jié)果(a) 地下半空間密度異常的理論分布; (b) 反演結(jié)果; (c) 訓(xùn)練過程的誤差曲線; (d) (a、b)的核密度分布曲線.

圖8 利用WtdMSE優(yōu)化模型訓(xùn)練的結(jié)果(a) 地下半空間密度異常的理論分布; (b) 反演結(jié)果; (c) 訓(xùn)練過程的誤差曲線; (d) (a、b)的核密度分布曲線.

3 訓(xùn)練結(jié)果及理論樣本反演結(jié)果

為展示2D-3D InvNet模型在不同形態(tài)樣本上的表現(xiàn),對(duì)表1設(shè)置的不同復(fù)雜度的測(cè)試樣本進(jìn)行反演,采用WtdMSE作為損失函數(shù).利用異常體中心塊體所在位置對(duì)地下半空間切片,可視化2D-3D InvNet的反演效果.

圖 10 樣本1反演結(jié)果在x、y、z方向的切片切片位置分別為沿x、y、z方向第10,第26和第7個(gè)網(wǎng)格處.

3.1 樣本1——單直立長(zhǎng)方體

樣本1地下半空間內(nèi)存在1個(gè)密度固定的直立長(zhǎng)方體異常體,異常體左上角網(wǎng)格坐標(biāo)為(6,9,2),長(zhǎng)、寬、高分別為8、11、9個(gè)單元網(wǎng)格,密度異常為1.5 g·cm-3(樣本參數(shù)設(shè)置見表1).圖9給出了訓(xùn)練好的2D-3D InvNet網(wǎng)絡(luò)對(duì)該樣本的反演結(jié)果,從圖9a可見反演結(jié)果異常體邊緣清晰,整體分布與真實(shí)地下分布一致,反演的密度值與真實(shí)值相當(dāng).圖9b顯示反演結(jié)果密度異常的概率密度分布與真實(shí)分布接近,二者核密度曲線包圍的面積基本相同.圖10的切片顯示網(wǎng)絡(luò)模型對(duì)異常體在x、y方向的邊界均有十分精準(zhǔn)的識(shí)別,z方向上也精準(zhǔn)地識(shí)別出了異常目標(biāo)上邊緣位置與主要塊體深度,下邊緣位置略有偏差.

圖9 樣本1密度反演結(jié)果(a) 左側(cè)為地下異常體理論分布、右側(cè)為其反演結(jié)果; (b) 理論分布與反演結(jié)果的核密度估計(jì)曲線.

總體上,反演結(jié)果在樣本形態(tài)、邊緣清晰度上均獲得了較好的結(jié)果,且核密度分布聚焦,曲線峰值及包圍面積與真實(shí)情況基本一致,這表明2D-3D InvNet能通過單個(gè)異常體造成的地表重力異常及重力梯度異常,有效識(shí)別地下簡(jiǎn)單的異常埋藏.

3.2 樣本2——雙直立長(zhǎng)方體

自然環(huán)境中的密度異常埋藏分布復(fù)雜,地表異常往往由多個(gè)異常體同時(shí)造成、互相影響.因此,需驗(yàn)證網(wǎng)絡(luò)模型利用較復(fù)雜地表異常識(shí)別并分離地下異常體的能力.樣本2在地下半空間存在2個(gè)直立長(zhǎng)方體密度異常,二者相互分離且密度異常均為1.3 g·cm-3(樣本參數(shù)設(shè)置見表1).圖11給出了2D-3D InvNet網(wǎng)絡(luò)對(duì)樣本2的反演結(jié)果,從圖11a可以看出反演結(jié)果很好地恢復(fù)了地下異常體的位置及形態(tài),且能識(shí)別出兩個(gè)相互分離的獨(dú)立異常體塊.圖11b顯示反演結(jié)果的密度異常分布與真實(shí)情況接近,二者曲線包圍的面積基本相同.從圖12的切片可以看出反演結(jié)果沿三個(gè)方向的密度異常分布基本均勻且邊緣較為清晰.

圖11 樣本2密度反演結(jié)果(a) 左側(cè)為地下異常體理論分布、右側(cè)為其反演結(jié)果; (b) 理論分布與反演結(jié)果的核密度估計(jì)曲線.

圖12 樣本2反演結(jié)果在x、y、z方向的切片切片位置分別為沿x、y、z方向第29,第26和第10個(gè)網(wǎng)格處.

綜上,盡管雙直立長(zhǎng)方體反演結(jié)果的邊緣清晰度略低于單直立長(zhǎng)方體,但能有效識(shí)別相互分離的異常體,且核密度估計(jì)曲線相對(duì)聚焦,曲線峰值及曲線包圍面積與真實(shí)情況基本一致,表明2D-3D InvNet能通過地表復(fù)較雜的重力異常及重力梯度異常,有效反演得到相互分離的異常塊體.

3.3 樣本3——階梯形樣本

以上展示了2D-3D InvNet對(duì)簡(jiǎn)單理論樣本的反演效果.真實(shí)世界地下埋藏的形態(tài)往往非常復(fù)雜,劉銀萍等(2013)和Hou等(2015)曾采用階梯型異常體樣本,驗(yàn)證反演方法的泛化性能.為測(cè)試2D-3D InvNet在復(fù)雜、不規(guī)則異常體樣本上的表現(xiàn),本節(jié)對(duì)階梯型樣本造成的地表異常進(jìn)行反演.階梯型樣本由四個(gè)相互連接的長(zhǎng)方體塊組成,其密度異常為2.0 g·cm-3(樣本參數(shù)設(shè)置見表1).圖13給出了2D-3D InvNet 網(wǎng)絡(luò)對(duì)階梯形樣本的反演結(jié)果,從圖13a可以看出反演結(jié)果較好地恢復(fù)出了異常體的位置和形態(tài).圖13b顯示反演的異常體密度聚焦于1.88 g·cm-3,與樣本的密度真實(shí)值匹配較好,反演結(jié)果的核密度分布曲線相對(duì)聚焦,曲線分布與包圍面積同真實(shí)情況基本一致.圖14的切片顯示反演結(jié)果能較好地?cái)M合目標(biāo)樣本邊界.

圖13 樣本3密度反演結(jié)果(a) 左側(cè)為地下異常體理論分布、右側(cè)為其反演結(jié)果; (b) 理論分布與反演結(jié)果的核密度估計(jì)曲線.

圖14 樣本3反演結(jié)果在x、y、z方向的切片切片位置分別為沿x、y、z方向第29,第27和第10個(gè)網(wǎng)格處.

因此,2D-3D InvNet不僅能通過地表異常反演簡(jiǎn)單地下異常體分布,還能準(zhǔn)確識(shí)別深度變化的復(fù)雜異常體,這進(jìn)一步證明了2D-3D InvNet反演方法的可靠性.

4 2D-3D InvNet反演方法在澳大利亞Kauring重力試驗(yàn)場(chǎng)的應(yīng)用

Kauring航空重力試驗(yàn)場(chǎng)位于西澳大利亞珀斯以東約100 km處(圖15),因豐富的礦產(chǎn)資源而聞名.Kauring試驗(yàn)場(chǎng)內(nèi)已確定具有較大的重力異常和磁異常,澳大利亞政府允許個(gè)人或組織對(duì)此區(qū)域進(jìn)行研究,其重力數(shù)據(jù)可直接從官網(wǎng)獲取.本文所使用的重力梯度分量數(shù)據(jù)來自Fugro公司在2011年7月至2012年2月間,利用FALCON航空重力梯度儀對(duì)Kauring試驗(yàn)場(chǎng)開展的測(cè)量.所有測(cè)量數(shù)據(jù)均基于澳大利亞地心基準(zhǔn)(GDA94)坐標(biāo)系,并做了地形校正.

圖15 西澳大利亞地區(qū)kauring重力試驗(yàn)場(chǎng)位置示意圖

在已經(jīng)獲取的數(shù)據(jù)精度基礎(chǔ)上,為避免數(shù)據(jù)空洞對(duì)反演結(jié)果產(chǎn)生影響,以60 m為單位間隔對(duì)數(shù)據(jù)進(jìn)行了重采樣.部分學(xué)者曾對(duì)此區(qū)域做出研究,坐標(biāo)系正東方向505～505.5 km,正北方向6470.5～6471 km的區(qū)域存在典型的重力異常和重力梯度異常(田宇等,2019).根據(jù)數(shù)據(jù)精度、2D-3D InvNet網(wǎng)絡(luò)輸入,參考其他學(xué)者的研究,取數(shù)據(jù)范圍為東504.1～506.5 km,北6469.3～6471.7 km.以東向?yàn)閤軸正方向,北向?yàn)閥軸正方向,Kauring地區(qū)重力異常及重力梯度異常數(shù)據(jù)如圖16所示.

由于輸入數(shù)據(jù)精度與前文理論模擬時(shí)不同,此處根據(jù)原始網(wǎng)絡(luò)設(shè)計(jì)方法重新劃定單個(gè)網(wǎng)格尺寸為60 m×60 m×40 m,樣本集規(guī)模和劃分,以及樣本訓(xùn)練方式與前文一致.完成模型訓(xùn)練后,將Kauring重力異常和重力梯度異常數(shù)據(jù)投入計(jì)算,反演結(jié)果如圖17所示.最終得到的Kauring地區(qū)地下埋藏異常體的密度異常最大值約為2.2～2.3 g·cm-3,埋藏深度約在320～400m之間.根據(jù)他人(Martinez and Li, 2012; Liu et al., 2015; 田宇等, 2019)的研究結(jié)果,該區(qū)域存在密度異常為2.3～2.4 g·cm-3、深度位于320～400 m之間的地下異常體.本文結(jié)果與其他學(xué)者的反演結(jié)果一致.此外,我們還發(fā)現(xiàn)研究區(qū)東南方向存在密度異常高達(dá)2.7 g·cm-3的異常體塊.

5 結(jié)論與討論

針對(duì)利用地表二維重力異常和重力梯度異常反演地下三維密度異常體的問題,本文提出了2D-3D InvNet深度學(xué)習(xí)反演算法,編碼階段用二維操作提取信息,解碼階段創(chuàng)新性地利用三維操作還原異常體在地下半空間的分布及密度信息.經(jīng)歷了完整的訓(xùn)練過程后,模型能在不到1 s的時(shí)間內(nèi),得到數(shù)據(jù)反演結(jié)果.

利用MSE作為損失函數(shù)反演的密度結(jié)果與真實(shí)值存在較大差異,本文通過改進(jìn)誤差函數(shù),同時(shí)引入核密度估計(jì),評(píng)價(jià)密度反演結(jié)果,提高了地下異常體密度反演的準(zhǔn)確性.將常用的損失函數(shù)MSE改進(jìn)為加權(quán)后的損失函數(shù)WtdMSE,通過比對(duì)訓(xùn)練過程中異常體所在區(qū)域誤差值的變化,發(fā)現(xiàn)在位置和分布形態(tài)準(zhǔn)度相差不大的情況下,WtdMSE對(duì)于目標(biāo)所在區(qū)域最終的計(jì)算誤差值與MSE的相比低至少50%以上,這表明網(wǎng)絡(luò)訓(xùn)練采用WtdMSE作為損失函數(shù)具有優(yōu)越性.

對(duì)測(cè)試樣本進(jìn)行精度評(píng)估,發(fā)現(xiàn)2D-3D InvNet深度學(xué)習(xí)模型的反演結(jié)果,無論是在樣本邊界的清晰度還是密度分布上,均能得到較好的結(jié)果.將此方法應(yīng)用于西澳大利亞Kauring地區(qū),我們成功獲得了地下密度異常體的分布,中心塊體剩余密度反演結(jié)果與他人的研究結(jié)果相符.由此可見,2D-3D InvNet深度學(xué)習(xí)方法穩(wěn)定且具有較強(qiáng)的泛化性能,能在無需規(guī)定限制條件的情況下進(jìn)行快速反演,不僅能很好地反演模擬樣本的密度和形態(tài),且能在實(shí)際礦區(qū)地下密度異常體識(shí)別中應(yīng)用.

“人工智能”已經(jīng)成為當(dāng)前科技以及各領(lǐng)域跨學(xué)科研究的主要方向,地球科學(xué)也是此賽道上的成員之一.重力異常和重力梯度異常數(shù)據(jù)是平面的,可以劃歸為類似圖像的矩陣結(jié)構(gòu),因此跟圖像相關(guān)的深度學(xué)習(xí)算法也可以在此領(lǐng)域應(yīng)用.盡管目前已知地下埋藏的異常數(shù)據(jù)十分稀少,現(xiàn)有方法僅能通過理論數(shù)據(jù)訓(xùn)練模型,但仍能得到較好的反演結(jié)果.深度學(xué)習(xí)或可成為地球物理反演極具前景的發(fā)展方向之一.

致謝西南交通大學(xué)張志厚老師對(duì)相關(guān)問題給出了解答和幫助,中國(guó)地震局第一監(jiān)測(cè)中心劉金釗老師提供了數(shù)據(jù)支持,兩位審稿專家提出了寶貴的修改意見,作者在此表示衷心的感謝.