李思平,劉彩云,熊杰,田慧瀟,王方
(1.長江大學(xué) 電子信息學(xué)院,湖北 荊州 434023;2.長江大學(xué) 信息與數(shù)學(xué)學(xué)院,湖北 荊州 434023)
大地電磁測深是以天然交變電磁場為場源來研究地球內(nèi)部電性結(jié)構(gòu)的一種重要的地球物理勘探方法,其反演問題一直是地球物理工作者研究的重點(diǎn)[1]。常用的反演方法,如OCCAM反演[2]、快速松弛RRI反演[3]、自適應(yīng)正則化反演[4]等在大地電磁數(shù)據(jù)處理中占有重要席位,然而這些方法都受非線性和非唯一性的限制,在計(jì)算過程中要依賴先驗(yàn)信息和初始模型的選擇,且有陷入局部極小的風(fēng)險(xiǎn)。隨著地球物理反演方法的發(fā)展,差分進(jìn)化法(DE)[5]、蟻群算法(ACO)[6]和人工神經(jīng)網(wǎng)絡(luò)[7]等全局尋優(yōu)的反演方法進(jìn)入人們的視野,這些方法不依賴于初始模型,但存在反演時(shí)間長、計(jì)算量大的問題,且針對二維反演或三維反演,效果還有待改進(jìn)。
與以上反演方法相比,基于深度學(xué)習(xí)(DL)的完全非線性反演方法無需計(jì)算偏導(dǎo)數(shù)和靈敏度矩陣,對初始模型依賴性小且具有良好的全局尋優(yōu)性能,因此在地球物理數(shù)據(jù)處理中受到極大關(guān)注。例如:Noh等[8]采用深度神經(jīng)網(wǎng)絡(luò)(DNN)對一維的頻率域航空電磁數(shù)據(jù)進(jìn)行了反演,相比于高斯-牛頓反演方法,盡管合成的訓(xùn)練集數(shù)據(jù)量較少,其結(jié)果仍然具有較高的分辨率;Liu等[9]提出了一種基于深度殘差網(wǎng)絡(luò)(ResNet18)的一維音頻大地電磁反演方法;Wang等[10]實(shí)現(xiàn)了基于深度置信網(wǎng)絡(luò)(DBN)的大地電磁測數(shù)據(jù)的二維非線性反演;廖曉龍等[11]提出一種基于卷積神經(jīng)網(wǎng)絡(luò)(CNN)的大地電磁反演方法;范振宇[12]以深度殘差網(wǎng)絡(luò)(ResNet50)為基礎(chǔ),實(shí)現(xiàn)了大地電磁反演并解決了實(shí)測問題。上述反演方法均克服了傳統(tǒng)方法的不足且在一定程度上提高了大地電磁反演的計(jì)算效率,但也存在一定的問題,例如在二維反演方面,有的網(wǎng)絡(luò)結(jié)構(gòu)簡單[11],受到網(wǎng)絡(luò)層數(shù)的限制導(dǎo)致表達(dá)能力有限,反演準(zhǔn)確性、探測深度和縱向分辨率有待進(jìn)一步提高,并且理論模型反演的視電阻率數(shù)據(jù)都是“學(xué)習(xí)”過的,對于未“學(xué)習(xí)”過的數(shù)據(jù),其反演效果(即泛化能力)如何值得進(jìn)一步研究;有的可以解決大尺度實(shí)測數(shù)據(jù)反演問題,但“以分類替代回歸解決實(shí)際問題,無法完全提取數(shù)據(jù)的特征與標(biāo)定電阻率值精細(xì)擬合”[12]。
因此,本文提出一種改進(jìn)的殘差網(wǎng)絡(luò)(iResNet),通過對多種二維地電模型進(jìn)行網(wǎng)絡(luò)學(xué)習(xí)訓(xùn)練和反演測試,并將反演結(jié)果與CNN網(wǎng)絡(luò)進(jìn)行對比,驗(yàn)證該網(wǎng)絡(luò)對大地電磁非線性反演的可行性和準(zhǔn)確性,評估網(wǎng)絡(luò)的泛化能力和抗噪聲能力,最終應(yīng)用于實(shí)測數(shù)據(jù)解決實(shí)際問題。
卷積神經(jīng)網(wǎng)絡(luò)(CNN)是深度學(xué)習(xí)的代表算法之一, 可實(shí)現(xiàn)卷積計(jì)算?;A(chǔ)的卷積神經(jīng)網(wǎng)絡(luò)由卷積層、池化層、全連接層三種結(jié)構(gòu)組成,網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示。
圖1 卷積神經(jīng)網(wǎng)絡(luò)[13]
對于卷積神經(jīng)網(wǎng)絡(luò)而言,網(wǎng)絡(luò)越深準(zhǔn)確率越高,然而過深的網(wǎng)絡(luò)不僅容易出現(xiàn)梯度消失和梯度爆炸的問題,還會(huì)出現(xiàn)“退化現(xiàn)象”,即當(dāng)網(wǎng)絡(luò)達(dá)到一定的深度后繼續(xù)增加層數(shù),模型的準(zhǔn)確率卻會(huì)下降。于是He等學(xué)者提出了殘差網(wǎng)絡(luò)[14]。殘差網(wǎng)絡(luò)不僅可以有效地解決神經(jīng)網(wǎng)絡(luò)中梯度爆炸和梯度消失的問題,還通過引入“快捷連接”來解決 “退化現(xiàn)象”?!翱旖葸B接”能將輸入x添加到經(jīng)過兩個(gè)權(quán)重層之后的輸出中,如圖2所示,因此輸出H(x)等于輸入的映射F(x)加上輸入x,該模型即為殘差模塊。
圖2 殘差網(wǎng)絡(luò)的一般結(jié)構(gòu)Fig.2 General structure of ResNet
殘差網(wǎng)絡(luò)采用的不是直接學(xué)習(xí)堆疊網(wǎng)絡(luò)的潛在映射H(x),而是增加一個(gè)恒等映射來擬合殘差映射F(x)=H(x)-x,F(x)包含2個(gè)權(quán)值層和1個(gè)ReLU(linear rectification function),易于優(yōu)化。具體在構(gòu)造ResNet時(shí),F(x)可以根據(jù)實(shí)際要求采用不同的殘差塊。如圖3展示的兩種典型的殘差結(jié)構(gòu),基于基本塊(residual block)和基于瓶頸塊(bottleneck block),用于構(gòu)建基于深度卷積網(wǎng)絡(luò)的ResNet。其中,n、i分別表示經(jīng)過1×1或2×2卷積后的特征圖數(shù)量,也就是通道數(shù)量。圖3a與圖3b的不同在于后者具有3個(gè)卷積層,當(dāng)輸入與輸出的維度不相同時(shí)需要用圖3b來調(diào)整維度。
圖3 兩種典型的殘差結(jié)構(gòu)
目前,CNN 的主要應(yīng)用集中在圖像識別、目標(biāo)檢測等領(lǐng)域[15],處理的目標(biāo)通常是圖像。在地球物理領(lǐng)域,利用CNN實(shí)現(xiàn)大地電磁二維反演時(shí),以二維視電阻率作為輸入數(shù)據(jù),二維地電模型參數(shù)作為網(wǎng)絡(luò)的輸出,核心思想是用卷積神經(jīng)網(wǎng)絡(luò)建立從輸入(視電阻率)到輸出(地電模型)之間的映射關(guān)系,即通過求解以下優(yōu)化問題來實(shí)現(xiàn)。
用數(shù)學(xué)公式可以這樣描述這個(gè)過程:給定一組視電阻率數(shù)據(jù)d和地電模型m,通過損失函數(shù)L去調(diào)整參數(shù)θ,最終獲得參數(shù)化模型Net(d,θ)。神經(jīng)網(wǎng)絡(luò)參數(shù)更新表達(dá)式為[16]
(1)
式中:θ表示網(wǎng)絡(luò)中需要更新的權(quán)重和偏差;Net(·)卷積神經(jīng)網(wǎng)絡(luò)表示從視電阻率數(shù)據(jù)d到預(yù)測地電模型mpred之間的映射;L是MSE損失函數(shù)。上式的意義在于優(yōu)化θ,使得預(yù)測地電模型mpred和真實(shí)地電模型m之間的誤差最小。
本文的目標(biāo)函數(shù)設(shè)置為
(2)
如式(2)所示,本文目標(biāo)函數(shù)由兩部分組成:損失函數(shù)和正則化項(xiàng),其中,損失函數(shù)用來均衡模型的預(yù)測結(jié)果與真實(shí)值之間的誤差;正則化項(xiàng)用來避免模型過擬合。具體來說,損失函數(shù)為理論模型與反演結(jié)果的均方誤差,對應(yīng)目標(biāo)函數(shù)第一項(xiàng);正則化項(xiàng)為網(wǎng)絡(luò)權(quán)重約束項(xiàng),對應(yīng)目標(biāo)函數(shù)第二項(xiàng)。如此設(shè)置,能使神經(jīng)網(wǎng)絡(luò)不僅擁有一定的學(xué)習(xí)能力,還能增強(qiáng)其泛化能力。
反演是由觀測數(shù)據(jù)反推模型參數(shù)的映射,具體實(shí)現(xiàn)過程如圖4所示,主要包括4個(gè)步驟:構(gòu)建數(shù)據(jù)集、搭建網(wǎng)絡(luò)、訓(xùn)練網(wǎng)絡(luò)與數(shù)據(jù)反演。
圖4 反演過程
1)利用有限單元法計(jì)算不同地電模型的視電阻率,將視電阻率和模型對組成數(shù)據(jù)集,并按一定的比例分為訓(xùn)練集和驗(yàn)證集。
2)設(shè)計(jì)iResNet超參數(shù)建立網(wǎng)絡(luò)模型,超參數(shù)主要包括:網(wǎng)絡(luò)模型結(jié)構(gòu)、卷積層、殘差模塊的窗口和步長等。
3)分別利用訓(xùn)練集和驗(yàn)證集對網(wǎng)絡(luò)進(jìn)行訓(xùn)練和測試,不斷優(yōu)化網(wǎng)絡(luò)參數(shù)。
4)將未“學(xué)習(xí)”過的視電阻率數(shù)據(jù)輸入到已訓(xùn)練好的iResNet網(wǎng)絡(luò),輸出對應(yīng)的地電模型,實(shí)現(xiàn)反演。
通過圖3所示的殘差結(jié)構(gòu),本文借鑒經(jīng)典的殘差網(wǎng)絡(luò)ResNet18,改進(jìn)得到一種新的反演網(wǎng)絡(luò)iResNet。該網(wǎng)絡(luò)主要由4個(gè)殘差塊和2個(gè)全連接層組成,其中每個(gè)殘差塊由4個(gè)卷積核和2個(gè)快捷連接組成,在Block2、3、4中由于輸入和輸出的數(shù)據(jù)維度不同,因而第一個(gè)快捷連接的方式存在差異。網(wǎng)絡(luò)具體參數(shù)如圖5所示[9]。
圖5 iResNet網(wǎng)絡(luò)結(jié)構(gòu)[9]
在二維大地電磁TM極化模式下,設(shè)計(jì)了5×5(500 m×500 m)、6×6(600 m×600 m)、8×4(800 m×400 m)、4×8(400 m×800 m)、5×10(500 m×1 000 m)、10×5(1 000 m×500 m)共7個(gè)形狀規(guī)則且單一的模型,以及由3個(gè)4×8堆疊的不規(guī)則梯形,并針對每種模型均設(shè)置4個(gè)阻值:高阻3 000、2 500 Ω·m,低阻500、200 Ω·m,部分模型如圖6所示。
圖6 模型示意
在正演過程中,有限單元法目標(biāo)區(qū)域網(wǎng)格規(guī)模設(shè)置為32×57,網(wǎng)格間距為100 m。考慮到計(jì)算量的問題,異常體和圍巖的電阻率值不變,只通過改變異常體位置獲得用于網(wǎng)絡(luò)訓(xùn)練和反演測試的二維數(shù)據(jù),其中異常體移動(dòng)位置遍布整個(gè)網(wǎng)格,每次移動(dòng)間距為200 m,共生成8 987組數(shù)據(jù)構(gòu)成數(shù)據(jù)集。計(jì)算時(shí),選擇53個(gè)測點(diǎn),并記錄27個(gè)頻點(diǎn)的視電阻率。將采集范圍內(nèi)所有測點(diǎn)和頻點(diǎn)電阻率作為網(wǎng)絡(luò)輸入,則輸入為27×53的二維數(shù)據(jù)矩陣;將地電模型參數(shù)作為網(wǎng)絡(luò)輸出,則輸出層神經(jīng)元個(gè)數(shù)為32×57。
對上文中獲得的數(shù)據(jù)集進(jìn)行歸一化處理。數(shù)據(jù)的歸一化操作可有效提升網(wǎng)絡(luò)訓(xùn)練效率,避免梯度爆炸的產(chǎn)生。采用線性函數(shù)歸一化方法實(shí)現(xiàn)樣本數(shù)據(jù)的歸一化,將數(shù)據(jù)歸一化到[0,1],其表達(dá)公式為
(3)
在TM模式下,對上文3.1節(jié)中的模型進(jìn)行反演測試。通常,訓(xùn)練集用于訓(xùn)練網(wǎng)絡(luò)的權(quán)重與偏差,而驗(yàn)證集用于驗(yàn)證訓(xùn)練后網(wǎng)絡(luò)的準(zhǔn)確性。在本文中,對于iResNet網(wǎng)絡(luò),一共有8 987組數(shù)據(jù)樣本,隨機(jī)抽取90%(8 088組)作為訓(xùn)練集,剩下10%(899組)作為驗(yàn)證集。表1為訓(xùn)練網(wǎng)絡(luò)時(shí)的參數(shù)設(shè)置。
表1 iResNet網(wǎng)絡(luò)參數(shù)設(shè)置
通過調(diào)整參數(shù)進(jìn)行迭代訓(xùn)練。圖7給出了隨機(jī)抽取的幾組反演結(jié)果,可以看出:iResNet網(wǎng)絡(luò)反演精度較高,能準(zhǔn)確反演出模型的具體位置、形狀以及電阻率值,且反演精度受探測深度影響較小,具有較好的縱向分辨率。
圖7 部分驗(yàn)證集樣本的反演結(jié)果
圖8為廖曉龍等[11]用CNN網(wǎng)絡(luò)得到的反演結(jié)果,其中用于反演的低阻異常體電阻率為100 Ω·m,形狀為5×5(1 km×1 km),背景電阻率為1 000 Ω·m;地塹模型凹陷部分電阻率為100 Ω·m,長度為6格(1.2 km),寬度為11格(2.1 km);基底電阻率為1 000 Ω·m。在保持背景電阻率、異常體形狀及電阻率均與上述參數(shù)相同的情況下,本文使用iResNet網(wǎng)絡(luò)的反演結(jié)果如圖9所示。對比來看,CNN網(wǎng)絡(luò)得到的反演模型隨著所處位置的加深,其形狀逐漸模糊且邊緣電阻率值與真實(shí)值相差較大,而本文使用的iResNet網(wǎng)絡(luò)得到的反演模型其形狀以及電阻率值不會(huì)受到深度的影響,在反演準(zhǔn)確性、探測深度和縱向分辨率方面均優(yōu)于CNN網(wǎng)絡(luò)。
圖8 CNN的反演結(jié)果[11]
圖9 iResNet的反演結(jié)果
由于訓(xùn)練集和驗(yàn)證集中的模型形狀和電阻率值都是已“學(xué)習(xí)”過的,所以驗(yàn)證集的反演結(jié)果只能用來評估網(wǎng)絡(luò)的質(zhì)量,而不能評估網(wǎng)絡(luò)的泛化能力。因此本文設(shè)計(jì)了多個(gè)未“學(xué)習(xí)”過的模型來評估iResNet網(wǎng)絡(luò)的泛化能力。
3.5.1 實(shí)驗(yàn)一
采用“學(xué)習(xí)”過的電阻率參數(shù)(3 000、2 500、500、200 Ω·m)搭配未“學(xué)習(xí)”過的形狀(6×12、7×7、4個(gè)4×8堆疊的不規(guī)則梯形)構(gòu)成測試集來進(jìn)行評估。隨機(jī)抽取幾組反演結(jié)果,根據(jù)圖10a~d中單一規(guī)則模型的反演結(jié)果來看,反演得到的模型其位置與真實(shí)模型基本一致,模型四周少量單元格的電阻率值與實(shí)際值存在差異,但均可較好地確定模型的整體形狀;圖10e~f中,4個(gè)4×8堆疊的不規(guī)則梯形反演結(jié)果的位置及形狀幾乎與真實(shí)模型一致,少許單元格的電阻率值略低于實(shí)際值。
圖10 實(shí)驗(yàn)一的反演結(jié)果
3.5.2 實(shí)驗(yàn)二
采用“學(xué)習(xí)”過的形狀(6×6、4×8、3個(gè)4×8堆疊的不規(guī)則梯形)搭配未“學(xué)習(xí)”過的電阻率參數(shù)(2 800、2 300、700 Ω·m)構(gòu)成測試集來進(jìn)行評估。
圖11為隨機(jī)抽取的幾組反演結(jié)果,可見反演得到的模型位置、形狀以及電阻率值均與真實(shí)模型基本一致。
圖11 實(shí)驗(yàn)二的反演結(jié)果
3.5.3 實(shí)驗(yàn)三
為了進(jìn)一步評估本方法的泛化能力,實(shí)驗(yàn)三采用未“學(xué)習(xí)”過的形狀(7×7、4個(gè)4×8堆疊的不規(guī)則梯形)搭配未“學(xué)習(xí)”過的電阻率參數(shù)(2 300、700 Ω·m)構(gòu)成測試集來進(jìn)行評估。反演結(jié)果(圖12)表明, 反演得到的模型位置以及形狀與真實(shí)模型一致,電阻率數(shù)值與真實(shí)值接近,僅少許單元格的電阻率值低于實(shí)際值。由此說明本文iResNet網(wǎng)絡(luò)可以用來解決大地電磁反演問題,并且具有較好的泛化能力。
圖12 實(shí)驗(yàn)三的反演結(jié)果
在真實(shí)環(huán)境中,大地電磁通常受噪聲等因素的影響,因此將高斯白噪聲添加到視電阻率數(shù)據(jù)中,以分析網(wǎng)絡(luò)的抗噪能力。抗噪能力實(shí)驗(yàn)采用前文3.1節(jié)中的7種模型,正演時(shí)在視電阻率數(shù)據(jù)中分別加入3%和5%的高斯白噪聲(圖13),構(gòu)成數(shù)據(jù)集用以訓(xùn)練網(wǎng)絡(luò),并用實(shí)驗(yàn)一中設(shè)計(jì)的未“學(xué)習(xí)”過的模型來評估網(wǎng)絡(luò)的抗噪能力。
由圖14、圖15可以看出,無論是加入3%還是5%的高斯白噪聲,反演結(jié)果受噪聲的影響均較小,反演得到的模型形狀、位置都接近于真實(shí)模型,僅模型邊緣較少單元格的電阻率值略低于實(shí)際值,說明本文所提網(wǎng)絡(luò)具有較好的抗噪聲能力。
圖14 加入3%高斯白噪聲后的反演結(jié)果
以非洲南部實(shí)測大地電磁數(shù)據(jù)[12]為例,探討驗(yàn)證iResNet網(wǎng)絡(luò)的可行性與反演結(jié)果的可靠性。圖16所示為非洲南部的剛果克拉通(Congo craton)和達(dá)馬拉—杭濟(jì)—喬貝帶(Damara—Ghanzi—Chobe Belt)的詳細(xì)地質(zhì)圖。前人認(rèn)為在ETO007~ETO009間可能存在達(dá)馬拉造山帶與剛果克拉通南部間的分界線。本文重點(diǎn)研究ETO009及之后的6個(gè)測點(diǎn),即ETO009~ETO014[17]。這6個(gè)測點(diǎn)分布于Damara Northern Zone、Damara Northern Margin Zone、Northern Damara basement_Kamanjab Inlier3塊區(qū)域,前人將這3塊區(qū)域統(tǒng)稱為南部剛果克拉通(Southern Congo craton)的Northern zone[18]。
本文選擇以下27個(gè)頻點(diǎn)進(jìn)行研究:0.0092、0.0134、0.0183、0.0269、0.037、0.054、0.073、0.107、0.146、0.215、0.293、0.43、0.59、0.86、1.17、1.72、2.34、3.4、4.7、6.9、9.4、13.7、18.8、27.5、40、66、97 Hz。從圖17中的視電阻率數(shù)值可以看出,研究區(qū)域多為高阻。
圖17 Northern zone的視電阻率
圖18為iResNet網(wǎng)絡(luò)對Northern zone的反演結(jié)果,可以看出:該區(qū)域整體電阻率值較高,大約為104Ω·m,紅色弧線以上的區(qū)域電阻率值較低,約為103Ω·m,紅色弧線以下區(qū)域電阻率值高達(dá)104.5Ω·m?,F(xiàn)有研究結(jié)果表明,Kamanjab Inlier上層為元古代沉積物覆蓋層,內(nèi)層主要為花崗巖和片麻巖,因此呈現(xiàn)出上層電阻率值低,內(nèi)層電阻率值高(超過104Ω·m)的結(jié)果,且在ETO011和ETO012站點(diǎn)之間可以追蹤到內(nèi)層的表面邊界[20]。整體而言,本文反演結(jié)果與現(xiàn)有研究結(jié)果基本一致:iResNet網(wǎng)絡(luò)可以大體反映出地下物性的分布范圍與趨勢,并劃分電阻率值的變化邊界。
圖18 iResNet對Northern zone反演結(jié)果
為驗(yàn)證反演結(jié)果的準(zhǔn)確性,引用Khoza等[18]采用非線性共軛梯度法對非洲南部的KIM+ETO測線進(jìn)行反演的結(jié)果,圖19中紅色方框即為本文重點(diǎn)研究的區(qū)域。通過對比圖18與圖19,不難看出iResNet反演結(jié)果與正則化反演結(jié)果相似,但對于電阻率值極低的地質(zhì)構(gòu)造,iResNet則無法反演出非常準(zhǔn)確的結(jié)果。相較于前人使用殘差網(wǎng)絡(luò)對非洲南部大尺度實(shí)測數(shù)據(jù)進(jìn)行反演的結(jié)果[12]而言,本文所使用的iResNet在小尺度上擁有更加精細(xì)的反演結(jié)果,該反演結(jié)果也更加接近正則化反演結(jié)果。總體來說,本文方法可以有效用于大地電磁重建地下電阻率結(jié)構(gòu)。
圖19 KIM+ETO測線正則化反演結(jié)果[18]
本文提出一種基于改進(jìn)殘差網(wǎng)絡(luò)iResNet的大地電磁反演方法。實(shí)驗(yàn)結(jié)果表明,該方法可以快速準(zhǔn)確反演出地電模型的位置、形狀和電阻率值,具有較好的泛化能力和抗噪聲能力,在實(shí)際資料處理中具有一定的實(shí)用性,能用于大地電磁二維反演問題。但該方法也存在一定的局限性,例如本文的目標(biāo)函數(shù)中重點(diǎn)約束的是模型形狀,由此導(dǎo)致網(wǎng)絡(luò)會(huì)犧牲模型電阻率值的恢復(fù)以保證模型形狀的恢復(fù),因此在今后的工作中,可以進(jìn)一步修改目標(biāo)函數(shù)使其能對模型的電阻率值進(jìn)行約束,從而提高反演的精度。如今,殘差網(wǎng)絡(luò)應(yīng)用較為廣泛,但面對更加復(fù)雜的三維反演工作,如何進(jìn)行網(wǎng)絡(luò)結(jié)構(gòu)的優(yōu)化和數(shù)據(jù)的選取也值得進(jìn)一步深入研究。