基于UNet++卷積神經網絡的重力異常三維密度反演

李柏森,魯寶亮,2,3*,安國強,巨鵬,2,3,朱武,3,4,蘇子旺

1 長安大學地質工程與測繪學院,西安 710054 2 海洋油氣勘探國家工程研究中心,北京 100028 3 長安大學西部礦產資源與地質工程教育部重點實驗室,西安 710054 4 自然資源部生態地質與災害防控重點實驗室,西安 710054

0 引言

三維物性結構是近年來地球物理反演的熱點.而重磁方法作為勘探地球物理研究中的重要方向,是研究地殼結構與密度、磁性特征的主要方法,也是地球深部探測向全三維、大深度發展的重要方向.故通過建立合理的地球內部結構及物性賦存狀態的地質模型,進行重磁三維物性反演問題研究,對解決地球內部物理性質、結構、構造以及能源礦產分布等問題具有重要的科學研究意義.

利用重力異常進行三維物性反演主要采用基于直立長方體單元剖分的正則化反演(Aghaee et al.,2021; Chasseriau and Chouteau,2003; Codd et al.,2021; Commer,2011; Dramsch et al.,2021; Geng et al.,2020; Guo et al.,2021; Hou et al.,2016; Li and Oldenburg,1996,1998,2003; Meng et al.,2016; Nagihara and Hall,2001; Pilkington,2006; Sun and Li,2014; Sun et al.,2018; Wang et al.,2017; Zhdanov et al.,2004; 劉銀萍等,2013; 秦朋波和黃大年,2016; 姚長利等,2007).該方法首先將地下場源區域規則剖分為長方體單元,其次建立目標函數(含引入模型約束項),通過求解目標函數來反演確定剖分單元的密度.其中,模型約束有深度加權約束(Commer,2011; Li and Oldenburg,1996,1998,2003)、模型特征約束,其中分為光滑約束(Li and Oldenburg,1998,2003)、聚焦約束(Commer,2011)、稀疏約束(Pilkington,2006)、物性范圍約束(Commer,2011; Li and Oldenburg,2003)等,這些模型約束條件的加入對正則化問題的解有很大的改善,但是模型約束中正則化參數以及經驗選取的約束因子等對反演結果帶來了較大的影響和不確定性.同時在求解最優解的過程中也經常容易陷入局部極值狀態,無法得到全局最優解.

被證明可以擬合任意函數以及可以進行大數據處理的深度學習逐漸被應用于地球物理領域(胡琪鑫和徐亞,2022;劉俊等,2022;歐炳霖等,2023;王玥天等,2023;朱振宇等,2023).而大數據時代的到來以及近幾年硬件技術的發展,使得神經網絡迎來了新一輪研究和應用熱潮.對地球物理來說,深度學習可以不依賴先驗模型以及約束項、不需要存儲模型矩陣、極強的非線性復雜函數逼近性能正是實際處理當中所需要的(Yu and Ma,2021; 張志厚等,2021a,2022).

目前深度學習神經網絡在地震、測井、電法或者巖性分類方面的應用較為廣泛,而這些成果也為深度學習在重磁領域應用提供了思路.國內外有很多學者使用不同的網絡結構,例如BP神經網絡、卷積神經網絡(Convolutional Neural Networks,CNN)、全卷積神經網絡(Fully Convolutional Neural Networks,FCN)等,對數據網格化(馬國慶等,2021,2022)、二維(耿美霞和楊慶節,2013)或三維物性(Codd et al.,2021; Elawadi et al.,2001; Guo et al.,2021; Huang et al.,2021,2022; Lv et al.,2023; Wang et al.,2021; Yang et al.,2022; Zhang L Z et al.,2022; Zhang S et al.,2022; 管志寧等,1998; 張志厚等,2021a,b)、界面(He et al.,2021; Santos et al.,2015)以及特征位置(Poulton et al.,1992)等進行了反演方法研究,均取得了良好的效果.

但是之前使用的卷積神經網絡最后通常為全連接層,并且無法提取融合不同層次的數據特征,使得應用場景及準確率受限; 其次,對重力數據來說,通常在正負半區均有分布,且絕對值相對較大,所以之前使用的tanh、sigmoid或者ReLU會存在梯度消失問題; 此外學習率通常需要隨著訓練進行動態調整,而傳統使用的學習率衰減方法會使得陷入局部最優無法跳出,無法進行全局尋優.因此,不同的網絡結構、激活函數以及學習率對反演具有重要影響,在應用過程中如何選取仍需進一步開展研究.

在前人的研究基礎上,本文研究了基于改進的UNet++網絡的重力異常三維密度深度學習反演方法,使用了對輸入數據絕對值較大時梯度更為穩定、輸入值為負時輸出不為0的LeakyReLU激活函數,加入Bacth Normalization層對每層數據進行歸一化操作,使其分布與輸入數據相同,以增加網絡訓練的穩定性以及收斂速度,并采用了基于余弦退火的學習率變更策略以提高網絡的全局最優化能力,避免陷入局部極小值.反演流程主要包括模型數據集及標簽集的生成、深度神經網絡的構建及訓練以及模型預測.通過組合模型的試驗,與UNet、FCN以及正則化方法對比,驗證了本文方法具有較強的抗噪性及泛化能力,并通過對南澳大利亞Olympic Dam多金屬礦區的實測重力數據反演,進一步證明了本文方法的正確性以及有效性.

1 方法原理

1.1 UNet++網絡結構

神經網絡通過大量數據集和對應標簽的學習,來逼近數據集與標簽之間的非線性關系.對數據來說,通常需要提取不同層次的特征,才能更全面地學習逼近.對于傳統卷積神經網絡來說,最后通常為全連接層,就使得應用范圍受限.而隨著全卷積神經網絡的提出,通過卷積層和池化層對輸入數據進行處理和下采樣,學習不同層次的特征,而后通過反卷積層將其恢復至原始數據大小,這一系列也被稱為編碼-解碼架構.

UNet是最常用的編碼器-解碼器網絡,也是一種被使用較多的“端到端”網絡,已廣泛應用于各種場景(Ronneberger et al.,2015).過去大多數研究都將其作為主干,并對不同的任務進行了一些更改.UNet++是最具代表性的基于UNet的分割架構之一(Zhou et al.,2018).為了減少來自編碼器和解碼器子網絡的特征圖之間的語義差距,UNet++使用了一系列嵌套和密集的跳過路徑(圖1),而不僅僅是編碼器和解碼器網絡之間的連接.與UNet降采樣到最底層才返回不同,UNet++每進行一次降采樣,就將其上采樣回去并進行結果預測,最后再將不同尺度特征的預測結果進行融合輸出預測結果,從而使得可以提取不同大小層次的特征.相較于UNet,UNet++填滿了原本長鏈接部分的空白,通過抓取不同層次的特征進行整合疊加,從而實現對不同大小的目標的識別,提高網絡精度.

本文基于UNet++網絡架構搭建了如圖1所示網絡結構,網絡共有37層,其中包含18個卷積層(卷積核尺寸3×3,跨度1)、4個最大池化層(寬度2,跨度2)、9個反卷積層(卷積核尺寸2×2,跨度2)、4個跳躍鏈接層以及輸入輸出層,相對于CNN來說未使用全連接層,既極大減少了節點參數,也減少了二維數據一維化導致的數據特征丟失,提高了計算效率以及反演精度.

為了提高訓練網絡的穩定性、增加神經網絡的泛化能力,通常還會加入Dropout層或Batch Normalization層.Dropout層是指在訓練過程中以某個概率暫時地使某些神經元停止工作,使得每次訓練都類似于訓練不同的網絡結構,進而提高網絡的泛化能力.但全卷積神經網絡在執行過程中深層本身鏈接較少,使用Dropout放棄部分鏈接可能使得有效信息無法提取到.在訓練過程中,每個batch的訓練都會使得網絡的參數進行更新.在神經網絡較深時,后面數據的學習會受到前面數據分布的影響,而隨著深度增加,受影響程度會越來越大,此時就會使得每層學習的數據分布不同,網絡收斂速度變慢,泛化能力變差.本文通過加入Batch Normalization層對每層數據進行歸一化操作,使其分布與輸入數據相同,通過防止梯度爆炸或消失增加網絡訓練的穩定性、提高收斂速度(Ioffe and Szegedy,2015).

1.2 激活函數的選擇

激活函數是通過作用在神經元上,將輸入數據非線性地映射到輸出端,進而提高網絡的非線性擬合能力,使得網絡能力更強,可以擬合更復雜的數據,從而表征輸入輸出之間更復雜的函數關系.目前常用的激活函數有sigmoid(乙狀函數)、tanh(雙曲正切函數)、ReLU(Rectified Linear Unit,整流線性單元)以及LeakyReLU(Leaky Rectified Linear Unit,滲漏整流線性單元),其具體函數圖像見圖2.

圖2 激活函數曲線示意圖

神經網絡中參數的訓練更新主要依賴反向傳播算法,其核心內容是結合損失函數及最優化方法,而梯度類算法是目前應用最為廣泛的一種,在訓練過程中需要計算每一次殘差的梯度值進行更新.對重力數據來說,數值在正負半區均有分布,且變化較大,絕對值一般分布在幾十到上百.從函數圖像可以看出,在輸入絕對值大于2時,sigmoid、tanh變化平緩,梯度趨近于0,ReLU在輸入值為負時輸出恒為0,梯度為0,從而使得反向傳播時無法更新或更新較小,增加了收斂難度; 由于LeakyReLU通過解決輸入值為負時輸出恒為0,減少了梯度消失的可能性,所以其具有更好的收斂能力以及更小的計算成本(Xu et al,2015).本文使用表現更好的LeakyReLU作為網絡中的激活函數,并將在下文第2.2節中開展tanh、sigmoid、LeakyReLU對訓練結果影響的實驗研究.

1.3 基于余弦退火的學習率衰減

對神經網絡訓練來說,每輪次的更新步長會對最終訓練結果以及收斂速度產生較大影響.當使用梯度下降算法來優化目標函數的時候,在訓練初期學習率大一些,使得網絡收斂迅速;在訓練后期學習率小一些,使得網絡更好的收斂到最優解;當越來越接近Loss值的全局最小值時,學習率應該變得更小來使得模型盡可能接近這一點.

對于傳統衰減策略來說,通常是指定衰減方法,如線性衰減(Exponential Decay)、指數衰減(Linear Decay)等,使得學習率隨訓練過程逐漸減小從而使得結果收斂.因為目標優化函數可能不是嚴格的凸函數,所以可能存在多個局部最優解,而由于初始權重的隨機性,在訓練時梯度下降算法可能陷入局部最小值.此時可以通過突然提高學習率,來“跳出”局部最小值并找到通向全局最小值的路徑,所以本文還進行了余弦退火(Cosine Annealing)學習率的應用.通過模擬余弦函數,先快速下降,再線性上升,通過熱重啟的辦法使得尋優有可能跳出局部極小值(圖3),進而提高網絡訓練收斂的全局搜索能力(Loshchilov and Hutter,2017),計算公式為

(1)

圖3 不同學習率衰減策略及其應用示意圖

其中,i為退火循環次數,α為退火率,ηt為當前批次學習率,ηmin,ηmax為學習率的最大值與最小值,Tcur,Ti分別為上一次熱重啟后訓練輪次與給定熱重啟間隔輪次.

從圖3a可以看出,線性衰減和指數衰減均為單調減小至某個最小值后保持恒定,使得當陷入局部最小值時無法尋找到全局最優(圖3b),而余弦退火學習率衰減,可以在某一步時增大學習率,進而通過增大步長跳出局部最優,增強了全局尋優能力,訓練結束后保存所有輪次中的最優解.為進一步證明余弦退火的全局尋優能力,將在下文第2.2節中開展線性衰減、指數衰減、余弦退火對訓練結果影響的實驗研究.

1.4 損失函數及學習算法的選擇

損失函數與傳統地球物理反演方法最小化的目標函數相似,通過計算輸出結果與標簽集的誤差,利用反向傳播算法來修改網絡權值等參數,進而達到擬合效果.本文使用平均絕對誤差(MAE)與均方差(MSE)進行衡量:

(2)

強大的學習算法也是收斂至全局最優解必不可少的一環.核心算法之一是隨機梯度下降法(Stochastic Gradient Descent,SGD)(Hinton and Salakhutdinov,2006),但是其收斂速度慢,全局最優化能力較差,所以后續通過加入一階動量或二階動量的方法改進出了SGDM(SGD with Momentum)(Qian,1999)、AdaGrad(Adaptive Gradient)(Duchi et al.,2011).而為了解決AdaGrad二階動量持續累積導致的訓練提前結束問題,不計算全部歷史梯度而只關注一段時間內的梯度平均值,改進得到了RMSprop.而本文采用Kingma和Ba(2017)提出的Adam算法,通過同時使用一階動量以及二階動量,使得收斂速度變快,計算效率提升.

2 數據集生成與網絡訓練

2.1 數據集生成

對深度學習來說,通常需要大量的數據集.本文將地下半空間模型均勻剖分為直立六面體,在其中隨機生成長寬高不等的直立六面體組合及臺階狀、餅狀模型(圖4).通過給定剩余密度值以及六面體頂點坐標,即可利用正演公式(3)、(4)計算其引起的重力異常及其垂向一階導數(Nagy et al.,2000).由于只計算非零密度值的異常體,所以計算速度較快.再將總數據集分為三部分:訓練集、驗證集、測試集,訓練集是指神經網絡用于學習訓練網絡的數據集;驗證集用于確定網絡結構以及調整網絡中的超參數;測試集是用于檢驗網絡最終性能與泛化能力.公式(3)、(4)為

圖4 數據集模型示意圖

gz=Gρ{-(x-ξ)ln[r+(y-η)]-(y-η)ln[r+(x-ξ)]

(3)

(4)

數據集生成是將地下半空間模型均勻剖分為64×64×20個立方體,立方體邊長為100 m,在其中隨機生成長寬高不等的直立六面體組合及餅狀、不同厚度的傾斜板狀體模型,訓練集剩余密度在-1.0、-0.8、-0.6、0.6、0.8、1.0 g·cm-3中隨機取值,共生成模型16000組,利用正演公式計算其引起的重力異常生成數據集,其中訓練集、驗證集和測試集比例為80∶15∶5,即分別包含12800、2400、800組數據,每輪次訓練后進行一次驗證.

2.2 激活函數、學習率衰減策略及不同網絡結構的確定

由于不同學習參數會對訓練結果產生一定影響,本文選取了對訓練結果影響較大的三個參數進行實驗,分別是激活函數、學習率更新策略以及網絡結構.同樣使用前文生成的數據集,batchsize固定為64,給定最大訓練輪次200,同時若驗證集精度連續50輪次提升小于1×10-5則終止,訓練結束后保存所有輪次中最優的網絡參數,即以驗證集精度最高的網絡參數作為訓練結果并用于模型實驗.本文實驗硬件平臺為CPU:i9-10900X,GPU:RTX2080Ti×2,RAM:128GB,軟件平臺為Python3.9.12、TensorFlow2.9.0_GPU版.

首先是選取了三種激活函數,分別為LeakyReLU、sigmoid、tanh,使用UNet++網絡以及余弦退火學習率更新策略,不同激活函數的收斂曲線見圖5.從訓練集MSE收斂情況可以看出LeakyReLU同樣訓練輪次下收斂更快,收斂結果也優于sigmoid和tanh;從驗證集曲線來看使用LeakyReLU波動更小,穩定性更高,過擬合可能性要小于其他兩種.

圖5 不同激活函數均方差收斂曲線

其次選取了線性衰減、指數衰減以及余弦退火學習率更新進行實驗,使用UNet++網絡結構,LeakyReLU作為激活函數,最小學習率為1×10-7.MSE收斂曲線見圖6.在訓練剛開始階段三種方法收斂速度接近,但隨著訓練進行指數衰減和線性衰減下降趨于平緩,而由于余弦退火具有熱重啟策略,收斂曲線具有起伏,但收斂情況優于線性衰減和指數衰減.

對于網絡結構實驗,選取了常用的三種網絡模型,分別是UNet++、UNet以及FCN,使用LeakyReLU作為激活函數,以及余弦退火作為學習率更新方法,收斂曲線見圖7.從收斂曲線可以看出UNet++在迭代相同輪次的情況下訓練集和驗證集MSE收斂情況都要優于UNet以及FCN.

圖7 不同網絡結構均方差收斂曲線

綜上可以得出,UNet++網絡模型結合LeakyReLU作為激活函數,在使用余弦退火學習率更新策略的情況下具有較強的全局尋優及防止過擬合能力.

3 模型實驗

3.1 直立六面體組合

模型實驗一為兩個不同大小、密度、位置的直立六面體模型,具體參數及三維示意圖如表1所示,其重力異常如圖8a所示.為了測試噪聲對反演結果的影響,對重力異常添加了其原始數據絕對值最大值的3%、5%、10%的正態分布隨機噪聲,噪聲均值和方差見表2.

表1 直立六面體模型參數

表2 不同噪聲水平及反演方法結果對比

圖8 不同方法在不同噪聲水平數據的反演結果

利用UNet++、FCN、UNet網絡結構及傳統方法進行反演對比,不同的網絡結構均選取LeakyReLU激活函數、余弦退火學習率,即圖7訓練的網絡,傳統方法反演使用正則化方法,加入深度加權、物性約束以及最小長度約束,使用L2范數作為目標函數.反演結果的均方根誤差見表2,三維結果見圖8e—t.結果表明,在不加噪的情況下UNet++與FCN反演結果接近,優于UNet,能夠非常準確地刻畫真實模型的地下空間位置以及相對剩余密度大小,邊界及異常體形態清晰,且明顯優于傳統正則化方法.不同噪聲水平下的反演結果的均方根誤差(表2)及三維幾何形態(圖8),可以看出UNet++均優于其他兩種網絡,也仍然優于傳統方法,說明UNet++受噪聲影響較小,更為穩定,抗噪性較好.

3.2 復雜組合模型

該實驗為三種模型的組合,分別為兩個直立六面體、兩個餅狀模型、一個臺階狀模型,其具體參數及三維示意圖如表3所示,不同噪聲水平的重力異常見圖9a—d,噪聲的均值及方差見表4.

表3 復雜模型中各模型參數

表4 不同噪聲水平及反演方法結果對比

圖9 不同方法在不同噪聲水平數據的反演結果

不同網絡結構及傳統正則化方法反演結果見圖9e—t,均方根誤差見表4.由于該種模型組合未經訓練,所以反演結果均方根誤差均大于3.1節中直立六面體模型.在無噪聲情況下,UNet++網絡反演結果仍然可以較為準確的刻畫異常體形態特征及剩余密度情況,與FCN和UNet結果接近,均優于傳統正則化方法.隨著噪聲水平增加,三種網絡結構對餅狀模型的反演逐漸變差,但UNet++網絡對邊界刻畫明顯優于其他兩種網絡,FCN和傳統正則化方法對餅狀模型較難識別,UNet受噪聲干擾較大,虛假異常較多.

綜合模型實驗結果,從無噪聲的反演結果來看,本文提出的方法模型可以應用于重力異常的三維密度反演,且即使未經過訓練的組合模型也可以具有一定的反演效果,說明了該網絡具有較強的泛化能力,但反演精度會低于經過訓練的模型組合;從添加噪聲后的反演結果來看,即使在添加10%的噪聲后,也可以反演得到異常體模型的空間位置及相對剩余密度情況,且反演結果均優于FCN、UNet以及傳統正則化方法,說明本文所搭建的網絡具有一定的有效性及穩定性.

4 實測數據測試

為進一步檢驗本文所搭建網絡的可行性及實用性,利用Olympic Dam多金屬礦區布格重力數據進行了反演實驗.Olympic Dam多金屬礦區位于澳大利亞南部高勒克拉通省(圖10),是20世紀70年代發現的一個超大型Cu-U-Au-Ag礦床,同時也是典型的IOCG礦床(Iron Oxide-Copper-Gold礦床).礦床產于一個復雜的角礫巖群中,主要礦石為富含赤鐵礦和磁鐵礦的角礫巖(圖11).Olympic Dam多金屬礦床具有明顯的重力異常特征.為了消除背景場,獲得由礦體引起的剩余重力異常,使用二維小波多尺度分解方法(楊文采等,2001)對該區域的布格重力異常數據進行異常的分離,如圖12a所示為重力異常,圖12b為區域重力異常,圖12c為剩余重力異常.

圖10 高勒克拉通省的基底地質圖(Direen and Lyons,2007)

圖11 超巨型Olympic Dam角礫巖群的簡化地質圖及剖面圖(Kirchenbaur et al.,2016)

圖12 重力異常及其分離結果

利用本文所提方法及傳統正則化方法對該區剩余重力異常進行了三維密度反演,兩種方法反演得到的高密度巖體平面位置分布重合度較高,本文方法反演的巖體更為聚焦,邊緣更清晰,而正則化反演更為光滑.從本文方法反演結果(圖13a)來看,高密度巖體平面位置分布與剩余重力異常分布形態相似,縱向上主要分布在300～1400 m的深度范圍,并與Fe元素含量5%等值線(圖11a)吻合度較高.4口鉆井巖心實測密度曲線(圖14)顯示高密度巖體(實測密度大于2.9 g·cm-3)主要分布集中在400～1400 m范圍內;淺部上覆蓋層為低密度,小于2.6 g·cm-3,深部為一套變質的火山巖、沉積巖及侵入巖,密度相對礦體較小,在2.6～2.9 g·cm-3之間(Direen 和 Lyons,2007),反演結果與鉆井RD16W1井密度變化較為一致(圖15),其中該井在455～1316 m深度部分的密度為實測值,其余深度部分密度根據Direen和Lyons(2007)所給數據推測繪制.另外,本文方法反演結果與前人反演結果(Yang et al.,2019; Zhang S et al.,2022)及正則化反演結果(圖13c)的平面形態相似,但最大分布深度略有差別.以上分析表明本文方法反演結果水平位置及垂向深度與目標巖體位置重合度較高.另外,分別對深度學習及正則化方法反演結果進行了重力數據重建(圖16a、d),計算了重建數據與剩余重力異常的均方根誤差(圖16b、e)及殘差(圖16c、f),結果顯示兩種方法重建結果與剩余重力異常形態相似,重建效果均較好,UNet++重建結果均方根誤差為0.6433 mGal,略小于正則化方法重建結果.綜上,說明本文所提方法在實際數據處理中也表現出良好的可行性和實用性.

圖13 反演結果三維透視圖及切片圖

圖14 鉆井巖心實測密度曲線

圖15 過RD16W1井密度剖面(y=6063.68 km)

圖16 UNet++與傳統方法反演結果重建數據及殘差

5 結論

本文提出了基于UNet++網絡的重力異常三維物性反演,采用LeakyReLU激活函數以及加入Batch Normalization層,提高了網絡更新收斂速度和穩定性;引入基于余弦退火的學習率更新策略,明顯提升了網絡的全局最優化能力,進而實現了重力異常數據地下異常體模型的三維物性反演.將本文方法與UNet、FCN以及傳統正則化方法應用于不同噪聲水平模型實驗,結果表明本文搭建的網絡可以準確地反映地下異常體模型的三維物性,且未經訓練的模型也可以較為準確的反映異常體特征,說明網絡具有良好的泛化能力及較強的抗噪性.并將該方法應用于澳大利亞Olympic Dam多金屬礦區的實際重力數據,反演結果與地質情況、鉆井資料以及正則化反演結果吻合度較高,進一步表明本文所提出的方法具有良好的可行性及實用性.

但是深度學習方法通常需要足量、多樣且具有代表性的樣本,而增加樣本量會使得計算數據量及復雜度大為上升,如何通過小批量樣本訓練網絡達到相近結果是本文方法的一個局限性.本文反演未使用先驗模型信息及約束項,僅僅使用了一種網絡架構,而不同的網絡具有不同的優勢,也就會導致反演結果有所不同.在下一步工作中將優化網絡結構及樣本特征性,嘗試不同方法,加入約束項及先驗信息,進而增加網絡收斂能力及應用范圍.另外,在本文研究中發現采用深度學習方法的反演結果出現了垂向分層現象,在其他深度神經網絡重力反演結果中也出現了該現象(Yang et al.,2022; Zhou et al.,2023),具體原因仍需深入研究.

致謝感謝兩位匿名評審專家及編輯部提出的寶貴修改意見和建議,感謝南澳大利亞資源信息網站提供的Olympic Dam多金屬礦區布格重力數據及鉆井資料.