基于深度學習的二維航空大地電磁數(shù)據(jù)反演

李 賡, 曹飛翔

(河南理工大學物理與電子信息學院, 焦作 454150)

大地電磁測深法是一種天然源頻率域電磁法,其以自然產(chǎn)生的平面電磁波為場源,通過觀測地表正交的電磁場分量獲取地下的電阻率結(jié)構(gòu)[1]。然而,在地形較為復雜的地區(qū),地面大地電磁法施工效率較低,大面積勘探成本急劇上升。目前,航空大地電磁法具有機動靈活、成本低、效率高的優(yōu)點,能夠在復雜環(huán)境中進行大規(guī)模快速勘探,已被應用于油氣勘探、礦產(chǎn)資源勘探等領(lǐng)域[2-4]。

航空大地電磁法利用高空中電離層產(chǎn)生的平面波作為發(fā)射源,通過搭載在直升機上的探測設備采集垂直磁場信號,同時在地面基站采集水平磁場信號,獲得傾子數(shù)據(jù)。經(jīng)反演[5-6]計算,得到地下介質(zhì)的電性結(jié)構(gòu)。許智博等[7]采用非線性共軛梯度算法實現(xiàn)了二維航空大地電磁反演計算,大大加快了目標函數(shù)的收斂速度。蘇揚等[8]針對層狀介質(zhì)模型推導實現(xiàn)了廣義模型約束條件下時間域航空電磁一維反演。殷長春等[9]通過引入雙向約束實現(xiàn)了航空電磁擬空間約束反演。

目前的電磁數(shù)據(jù)反演算法基本都是基于吉洪諾夫正則化的反演算法,這類方法的主要工作是求解目標函數(shù)的偏導數(shù),偏導數(shù)的求解以牛頓法的變種(求解梯度向量、海森矩陣)為主,這會導致反演結(jié)果對反演初始模型依賴高,且易陷入局部極小值。針對傳統(tǒng)電磁反演方法存在的問題,通過引入智能非線性算法來解決電磁反演問題,已被證明有效可行。陳紫靜[10]采用粒子群算法順利完成了大地電磁測深資料反演。王鶴等[11]實現(xiàn)了基于遺傳神經(jīng)網(wǎng)絡的大地電磁數(shù)據(jù)反演。李棟[12]采用人工神經(jīng)網(wǎng)絡開展了航空瞬變電磁擬電阻率反演成像研究工作。邸龍[13]采用人工蜂群算法對不同類型的地質(zhì)開展了大地電磁非線性反演方法的研究。

現(xiàn)基于二維航空大地電磁法正演理論基礎(chǔ),首先,針對地下介質(zhì)中嵌入異常體的地電模型,采用有限元方法計算其傾子數(shù)據(jù),并結(jié)合地下電阻率值標簽,構(gòu)建深度學習樣本數(shù)據(jù)集。接著采用深度學習卷積神經(jīng)網(wǎng)絡進行二維航空大地電磁傾子數(shù)據(jù)反演,針對不同埋深的異常體地電模型,分別采用傳統(tǒng)電磁反演和深度學習反演方法開展二維航空電磁數(shù)據(jù)反演研究。

1 二維航空大地電磁正演理論

1.1 邊值問題

研究的直升機航空大地電磁法是頻率域電磁方法,滿足頻率域麥克斯韋方程組,其表達形式為

(1)

式(1)中：μ為介質(zhì)磁導率；ω為角頻率，rad/s；σ=1/ρ，S/m,表示介質(zhì)電導率；E為電場強度；H為磁場矢量；表示哈密頓運算符。

在二維地電模型中,磁場的垂直分量Hz只存在于TE模式中。因此,二維航空大地電磁法控制方程為

(2)

式(2)中：Ex為電場E在x軸方向上的分量。

式(2)可進一步表示為

(3)

式(3)中:u=Ex；τ表示自由電荷密度，C/m3；τ=1/iωμ；λ表示電導率，S/m，λ=σ。

為了求解式(3),還必須給出邊界條件?，F(xiàn)假設求解區(qū)域如圖1所示,研究對象為天然電磁場,場源在距離地面足夠遠的高空中,電磁場以平面波的形式垂直入射地下介質(zhì)。

圖1 TE模式求解區(qū)域Fig.1 Solution domain for TE polarization

TE模式下的邊界條件如下。

(1)上邊界AB離地面足夠遠,使異常場值在AB上為0,該處u=1。

(2)下邊界CD以下為均勻介質(zhì),其邊界條件為u/n+ku=0。其中：n表示外法線；k為復波數(shù)，為CD以下介質(zhì)的電導率。

(3)左右邊界AD、BC處的邊界條件為u/n=0。

綜上所述,邊值問題可總結(jié)為

(4)

式(4)中:Ω表示求解區(qū)域;AB、CD、AD、BC分別表示求解區(qū)域的上、下、左、右邊界。

1.2 有限元求解方法

與上述邊值問題等價的變分問題為

(5)

式(5)中:Ω為待求解區(qū)域面積；Г為Ω邊界。

有限單元法在求解上述變分問題時,首先將待求解區(qū)域離散為一系列互不重合的四邊形單元；然后在各離散單元內(nèi)用插值函數(shù)與單元中節(jié)點上未知電場值的乘積作為近似解,將近似解代入邊值問題后,在每個單元上可得一個復系數(shù)線性方程組。向該復系數(shù)線性方程組中加入求解邊界條件(定解條件),假設為u=b。最終可得到一個大型稀疏對稱復系數(shù)線性方程組,該大型稀疏對稱復稀疏線性方程組表示為

(6)

式(6)中：Ke和ue分別表示單元系數(shù)矩陣和向量；Ne表示單元總個數(shù)。采用直接求解法求解該方程組即可得到求解區(qū)域內(nèi)各節(jié)點電場值組成的向量u。

1.3 傾子響應計算公式

研究的二維航空大地電磁數(shù)據(jù)為傾子數(shù)據(jù)。當?shù)叵陆橘|(zhì)電阻率在水平方向上發(fā)生變化時,磁場的垂直分量Hz≠0,此時,Hz與Hx和Hy之間的復系數(shù)關(guān)系為

(7)

式(7)中：T=[Tzx,Tzy]稱為傾子。

TE模式下式(7)可進一步表示為

Hz=TzyHy

(8)

最終二維傾子響應計算公式為

(9)

2 傳統(tǒng)電磁反演與深度學習反演理論

2.1 傳統(tǒng)電磁反演

利用傳統(tǒng)電磁反演方法進行反演時,首先構(gòu)建出帶有約束項的觀測數(shù)據(jù)目標函數(shù),接著采用反演算法求解該目標函數(shù)梯度的極小值,最終獲得最優(yōu)的反演地電模型。傳統(tǒng)電磁反演方法將帶有約束項的觀測數(shù)據(jù)目標函數(shù)定義為

Φ(m)=[d-F(m)]TV-1[d-F(m)]+

λmTLTLm

(10)

式(10)中:d=[d1,d2,…,dN]為N維傾子數(shù)據(jù)向量;m=[m1,m2,…,mM]為M維模型參數(shù)向量;F(·)表示從模型向量到傾子數(shù)據(jù)向量映射的正演算子;V為一個用來描述觀測數(shù)據(jù)誤差的正定矩陣；λ為正則化因子；L表示二階差分矩陣。

傳統(tǒng)電磁反演方法的最終目標就是采用合適的反演算法求解上述目標函數(shù)梯度的極小值?，F(xiàn)采用非線性性共軛梯度(nonlinear conjugate gradients, NLCG)反演算法求解目標函數(shù)梯度的極小值。利用目標函數(shù)的梯度信息對模型向量進行修正,從而實現(xiàn)搜索目標函數(shù)極小值的迭代過程。

NLCG算法通過初始化反演模型m=m0,設置最大迭代次數(shù)為I,在每次迭代時,計算目標函數(shù)、梯度,從而計算出搜索方向Pl,采用線搜索的方法確定搜索方向Pl上的最優(yōu)步長α。

梯度計算公式為

gl=-2A(mref)TV-1[d-F(ml)]+2λLTLml

(11)

式(11)中：A為雅可比矩陣。

最優(yōu)步長計算公式為

(12)

模型向量通過如下方式修正:

(13)

式(13)中:Pl為搜索方向；α為搜索步長。

式(13)中搜索方向通過如下方式修正:

(14)

式(14)中：Cl為預調(diào)節(jié)矩陣；βl利用Polak-Ribiere方法計算，即

(15)

由于NLCG反演算法已經(jīng)是很成熟的方法,Rodi等[14]已針對其進行了詳細的介紹,這里不做贅述。

2.2 深度學習反演

2.2.1 卷積神經(jīng)網(wǎng)絡結(jié)構(gòu)

卷積神經(jīng)網(wǎng)絡(convolutional neural networks, CNN)是深度學習的代表算法之一,可實現(xiàn)卷積計算。目前,卷積神經(jīng)網(wǎng)絡已被廣泛應用于人像識別、手勢識別等領(lǐng)域[15]。卷積神經(jīng)網(wǎng)絡結(jié)構(gòu)如圖2所示。

卷積和池化層是卷積神經(jīng)網(wǎng)絡結(jié)構(gòu)的核心組成部分。一般在卷積和池化層之后還需加入一層全連接層,卷積層可實現(xiàn)輸入數(shù)據(jù)和卷積核的卷積操作。池化層的作用是降采樣,可實現(xiàn)特征降維。卷積層和池化層實現(xiàn)了輸入數(shù)據(jù)有效特征的提取,將經(jīng)過特征提取后的有效特征輸入全連接層,最終輸出預測值。

圖2 卷積神經(jīng)網(wǎng)絡示意圖Fig.2 The diagram of convolutional neural network

2.2.2 反演原理

基于卷積神經(jīng)網(wǎng)絡的航空大地電磁數(shù)據(jù)反演采用有監(jiān)督的機器學習技術(shù),即把觀測點電磁響應數(shù)據(jù)作為訓練數(shù)據(jù),同時將相對應的地下電性結(jié)構(gòu)作為標簽。研究目標是把數(shù)據(jù)與標簽之間的復雜關(guān)系利用卷積神經(jīng)網(wǎng)絡進行參數(shù)化,從而利用觀測的電磁響應數(shù)據(jù)獲得地下電性結(jié)構(gòu)電阻率值。

用數(shù)學公式可以這樣描述這個過程:給定一組電磁訓練數(shù)據(jù)XTrain(傾子實部和虛部數(shù)據(jù)),地下電性結(jié)構(gòu)數(shù)據(jù)YTrain。通過損失函數(shù)L去調(diào)整參數(shù)θ,使得輸入數(shù)據(jù)集為XTrain,輸出數(shù)據(jù)集為YTrain,最終獲得參數(shù)化模型f(θ;XTrain)。神經(jīng)網(wǎng)絡參數(shù)更新表達式為

θ*=argminL[YTrain-f(θ;XTrain)]

(6)

假設該參數(shù)化的模型可近似表示電磁響應數(shù)據(jù)與地電結(jié)構(gòu)標簽之間的關(guān)系,只要訓練和預測數(shù)據(jù)具有相同的分布且以相似的方式采樣,其對應的地下電性結(jié)構(gòu)便可由勘探數(shù)據(jù)Xpredict獲得,該預測值的表達式為

Ypredict=f(θ*;Xpredict)

(7)

換句話說,只要訓練數(shù)據(jù)和預測數(shù)據(jù)的勘探目標和環(huán)境相似,且數(shù)據(jù)采集和處理方式相同,就可以用訓練好的神經(jīng)網(wǎng)絡模型來獲得對應的地下電性結(jié)構(gòu)。

3 算例

3.1 卷積神經(jīng)網(wǎng)絡反演算例

基于二維航空大地電磁法正演理論基礎(chǔ),采用有限單元法計算不同埋深的異常體地電模型傾子響應值,神經(jīng)網(wǎng)絡樣本數(shù)據(jù)為傾子響應實虛部,標簽為各異常體模型地下電性結(jié)構(gòu)電阻率值。

3.1.1 樣本數(shù)據(jù)集生成

首先創(chuàng)建一個初始異常體地電模型,如圖3所示。該模型中異常體位于400 m埋深的地下介質(zhì)中,大小為400 m×400 m,電阻率為10 Ω·m,地下介質(zhì)電阻率為100 Ω·m。

然后將初始地電模型中異常體按照50 m間隔依次向右和向下移動,其構(gòu)建過程如圖4所示,共生成1 728組異常體地電模型。具體參數(shù)如下:所有地電模型的航空測線高度均位于地面以上100 m處,在y=-2 000～2 000 m的觀測范圍內(nèi),按照100 m間隔選擇41個觀測點采集磁場垂直分量。在y=-3 000 m處采集磁場水平分量。在1～1 000 Hz范圍內(nèi),對數(shù)等間隔選擇30個計算頻率。針對這些地電模型采用有限單元法計算其傾子響應數(shù)據(jù)。

圖3 異常體模型示意圖Fig.3 The schematic diagram of anomalous body model

圖4 不同位置異常體構(gòu)建示意圖Fig.4 The construction schematic diagram of anomalous bodies at different locations

最后針對各地電模型設置不同的地下電性結(jié)構(gòu)電阻率值標簽,在設置各地電模型標簽時,在y=-2 000～2 000 m,z=0～2 000 m的觀測范圍,采用四邊形網(wǎng)格按照50 m間隔剖分反演網(wǎng)格,最終各地電模型標簽由3 200個反演網(wǎng)格構(gòu)成,共生成1 728個地電模型標簽。

3.1.2 數(shù)據(jù)預處理

針對3.1.1節(jié)中獲取的1 728組樣本數(shù)據(jù),采用數(shù)據(jù)歸一化方法進行歸一化處理,數(shù)據(jù)的歸一化操作可有效提升神經(jīng)網(wǎng)絡訓練效率,避免梯度爆炸的產(chǎn)生。采用線性函數(shù)歸一化方法實現(xiàn)樣本數(shù)據(jù)的歸一化,將樣本集數(shù)據(jù)歸一化到0.1～0.9,其表達公式為

(8)

式(8)中：X為原始數(shù)據(jù)；Xmin表示樣本數(shù)據(jù)最小值；Xmax表示樣本數(shù)據(jù)最大值。按照0.8∶0.1∶0.1的比例劃分樣本集為訓練集、驗證集和測試集。

3.1.3 卷積神經(jīng)網(wǎng)絡模型搭建

通過在卷積和池化層之后拼接一層全連接層,完成卷積神經(jīng)網(wǎng)絡模型的搭建。其中卷積層激活函數(shù)為relu,卷積核大小為32×3×3,池化層采用最大池化方法,在池化層后接一層flatten層實現(xiàn)一維數(shù)據(jù)的轉(zhuǎn)化。采用adam優(yōu)化算法和Dropout正則化技術(shù)對模型進行訓練。最終建立的卷積神經(jīng)網(wǎng)絡模型參數(shù)如圖5所示。

3.1.4 神經(jīng)網(wǎng)絡模型訓練與測試

在進行卷積神經(jīng)網(wǎng)絡模型訓練時,將神經(jīng)網(wǎng)絡的最大訓練輪迭代次數(shù)設置為100。卷積神經(jīng)網(wǎng)絡訓練和驗證誤差隨迭代輪數(shù)變化趨勢如圖6所示,最終模型訓練誤差和驗證誤差分別為2.11×10-2、2.48×10-2。從測試集隨機選擇7組數(shù)據(jù)用于卷積神經(jīng)網(wǎng)絡模型測試,7組測試數(shù)據(jù)測試誤差如圖7所示,其測試誤差平均值為2.44×10-2。

從測試集中隨機選擇兩組樣本進行卷積神經(jīng)網(wǎng)絡反演結(jié)果分析。表1展示了兩組樣本卷積神經(jīng)網(wǎng)絡反演模型測試結(jié)果。表1中不同樣本對比數(shù)據(jù)來源于異常體的相同位置。從表1中可以看出,卷積神經(jīng)網(wǎng)絡反演電阻率與真實值更為接近,可有效表示真實電阻率。圖8所示為不同位置和埋深的異常體卷積神經(jīng)網(wǎng)絡反演結(jié)果,其對應的模型真實背景電阻率均為100 Ω·m,異常體電阻率為10 Ω·m。由圖8中可以看出，卷積神經(jīng)網(wǎng)絡輸出的地電結(jié)構(gòu)電阻率分布較為均勻,反演結(jié)果中異常體邊界較清晰,異常體的位置和大小與真實異常體較接近。

圖5 卷積神經(jīng)網(wǎng)絡模型參數(shù)Fig.5 Parameters of convolution neural network model

圖6 卷積神經(jīng)網(wǎng)絡反演平均絕對誤差Fig.6 Mean absolute error of inversion of convolution neural network

圖7 7組測試數(shù)據(jù)測試誤差Fig.7 Test error of seven test data

表1 卷積神經(jīng)網(wǎng)絡反演模型測試結(jié)果Table 1 Test result of convolution neural network inversion model

圖8 卷積神經(jīng)網(wǎng)絡反演結(jié)果Fig.8 Inversion result of convolution neural network

3.2 傳統(tǒng)電磁反演與深度學習反演

在反演階段,首先針對一個550 m埋深的異常體模型,分別采用傳統(tǒng)電磁反演方法和卷積神經(jīng)網(wǎng)絡反演方法進行反演測試。

采用傳統(tǒng)電磁反演方法進行傾子數(shù)據(jù)反演時,首先,將反演初始模型背景電阻率設置為與真實的背景電阻率相同的100 Ω·m,其反演結(jié)果如圖9所示,圖中黑色虛線方框表示異常體的真實位置。從圖9可以看出,反演結(jié)果對異常體的橫向邊界反演效果較好,但在垂向上存在明顯偏差。接著,將反演初始模型背景電阻率設置為與真實背景電阻率不同的500 Ω·m時,其反演結(jié)果如圖10所示,從圖10中可以看出,當反演初始模型與真實模型電阻率相差較大時,單獨的傾子數(shù)據(jù)反演結(jié)果不可靠。

圖9 異常體埋深550 m-傳統(tǒng)電磁反演方法反演結(jié)果Fig.9 Inversion result of 550 m buried depth of anomalous body by traditional electromagnetic method

圖10 初始背景電阻率500 Ω·m-傳統(tǒng)電磁反演方法反演結(jié)果Fig.10 Inversion result by traditional electromagnetic method with initial background resistivity-500 Ω·m

圖11 550 m異常體埋深-卷積神經(jīng)網(wǎng)絡反演結(jié)果Fig.11 Inversion result of 550 m buried depth of anomalous body by convolution neural network

采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,其反演結(jié)果如圖11所示,圖中黑色虛線方框表示異常體的真實位置。從圖11中可以看出,反演結(jié)果中異常體的大小和位置與真實情況符合程度較高。對比圖9和圖11,可以得出結(jié)論:相比于傳統(tǒng)電磁反演方法,采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,反演結(jié)果更加準確、可靠。

接著針對一個750 m埋深的異常體模型,分別采用傳統(tǒng)電磁反演和卷積神經(jīng)網(wǎng)絡反演方法進行傾子數(shù)據(jù)反演測試。采用傳統(tǒng)電磁反演方法進行傾子數(shù)據(jù)反演時,將反演初始模型背景電阻率設置為與真實的背景電阻率相同的100 Ω·m,其反演結(jié)果如圖12所示,圖中黑色虛線方框表示異常體的真實位置。從圖12中可以看出:當異常體埋深增加時,傳統(tǒng)電磁反演方法反演結(jié)果中異常體的大小和位置已明顯偏離真實位置,呈現(xiàn)出整體向下偏移趨勢,反演結(jié)果已無法有效反映真實異常體的真實情況。

采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,其反演結(jié)果如圖13所示,從圖13中可以看出:當異常體埋深增加時,卷積神經(jīng)網(wǎng)絡反演結(jié)果依然較好,反演結(jié)果中異常體大小和位置與真實情況非常接近。

最后,針對不同埋深異常體模型的傳統(tǒng)電磁反演結(jié)果與卷積神經(jīng)網(wǎng)絡反演結(jié)果進行定量分析。表2展示了不同埋深異常體模型傳統(tǒng)電磁反演方法反演結(jié)果與卷積神經(jīng)網(wǎng)絡反演結(jié)果。表2中不同樣本對比數(shù)據(jù)來源于異常體的相同位置。從表2中可以看出,采用傳統(tǒng)電磁反演方法進行傾子數(shù)據(jù)反演時,異常體反演電阻率較真實值偏大。相較于傳統(tǒng)電磁反演方法,采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,異常體反演電阻率與真實值偏差更小。

圖12 750 m異常體埋深-傳統(tǒng)電磁反演方法反演結(jié)果Fig.12 Inversion result of 750 m buried depth of anomalous body by traditional electromagnetic method

圖13 750 m異常體埋深-卷積神經(jīng)網(wǎng)絡反演結(jié)果Fig.13 Inversion result of 750 m buried depth of anomalous body by convolution neural network

表2 傳統(tǒng)電磁反演與卷積神經(jīng)網(wǎng)絡反演結(jié)果對比

綜上所述,針對不同埋深的異常體模型,采用傳統(tǒng)電磁反演方法進行傾子數(shù)據(jù)反演時,隨著異常體埋深的增加,傳統(tǒng)電磁方法反演結(jié)果已無法有效反映異常體的真實情況；采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,反演速度快、準確度高,反演結(jié)果與真實情況更加接近。最終結(jié)果表明,采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演,反演效果更好。

4 結(jié)論

首先從麥克斯韋方程組出發(fā),推導了傾子響應計算公式,基于二維航空大地電磁法正演理論基礎(chǔ),針對地下介質(zhì)中嵌入異常體的地電模型,采用有限元方法計算其傾子數(shù)據(jù),并結(jié)合地下電阻率標簽,構(gòu)建了深度學習樣本數(shù)據(jù)集。接著采用深度學習卷積神經(jīng)網(wǎng)絡實現(xiàn)了二維航空大地電磁傾子數(shù)據(jù)反演,針對不同埋深的異常體模型,分別分析了傳統(tǒng)電磁反演和深度學習反演結(jié)果。最終結(jié)論總結(jié)如下。

(1)針對深度學習樣本數(shù)據(jù)集,采用卷積神經(jīng)網(wǎng)絡進行二維航空電磁反演時,反演結(jié)果對異常體邊界刻畫較為清晰,異常體的位置和大小與真實異常體較符合,且反演電阻率值與真實值更為接近,可有效表示真實電阻率值。

(2)針對不同埋深的異常體模型,采用傳統(tǒng)電磁反演方法進行反演時,反演結(jié)果對初始模型依賴度高,當反演初始模型與真實模型電阻率相差較大時,單獨的傾子數(shù)據(jù)反演結(jié)果不可靠。隨著異常體埋深的增加,傳統(tǒng)電磁反演結(jié)果已無法有效反映異常體的真實情況；采用卷積神經(jīng)網(wǎng)絡進行傾子數(shù)據(jù)反演時,反演速度快、準確度高,反演結(jié)果與真實情況更加接近。相較于傳統(tǒng)電磁反演方法,采用卷積神經(jīng)網(wǎng)絡方法進行傾子數(shù)據(jù)反演,反演效果更好。