廖曉龍,張志厚,姚禹,路潤琪,范祥泰,曹云勇2,馮濤2,石澤玉
(1.西南交通大學地球科學與環境工程學院,四川成都,611756;2.中鐵二院地勘巖土工程有限責任公司,四川成都,610031)
大地電磁測深是一種利用天然交變電磁場研究地球電性結構的地球物理勘探方法,具有勘探深度大、不受高阻屏蔽、對低阻層反應靈敏等優點,已經廣泛運用于礦產油氣勘探、地熱資源勘查和環境工程等領域[1]。反演是大地電磁數據處理過程中的一類重要技術,已有的大地電磁反演方法如OCCAM 反演[2]、快速松弛反演(RRI)[3]、非線性共軛梯度反演(NLCG)[4]及其改進的REBOCC 反演[5]、SBI 反演[6]、自適應正則化反演[7]、最小二乘正則化反演[8-9]、三維快速松弛反演[10]等,都在大地電磁資料處理與解釋中發揮著重要作用,但此類方法對初始模型的依賴性較強,當目標為多峰值的復雜函數時,反演容易陷入局部極小,在一定程度上影響反演的精確性和穩定性。近年來,一些基于神經網絡模型的非線性反演方法在地球物理數據處理中受到研究者的極大關注[11],此類方法不依賴初始模型的選取,無需計算靈敏度矩陣,具有良好的全局尋優性能,如:徐海浪等[12]采用神經網絡實現了直流電法二維電阻率反演;李創社等[13]將BP 神經網絡應用于瞬變電磁勘探;耿美霞等[14]采用RBF 神經網絡反演二維重力密度分布;王鶴等[15]搭建了BP 神經網絡模型并用遺傳算法進行優化,實現了遺傳神經網絡的大地電磁二維反演;戴前偉等[16]將BP神經網絡與粒子群算法結合,應用于電阻率成像非線性反演。以上研究表明神經網絡能夠以一定精度逼近地球物理參數與介質參數模型之間的數學關系,但這些研究均是基于淺層神經網絡模型,在一些更加復雜的非線性特征映射問題上逼近能力仍有不足。針對這一問題,HINTON 等[17]提出了深度神經網絡(deep neural networks,DNN)模型,DNN 在特征提取和建模上有明顯的優勢,能夠從原始的輸入數據中挖掘出更深層次的特征,具有超強的復雜函數逼近性能。各種類型的深度網絡在斷層檢測[18-19]、地震相分類[20]、數字巖石物理[21-22]等方面展現出廣闊的應用前景。 卷積神經網絡(convolutional neural networks,CNN)作為深度神經網絡的一種類型,其局部連接、權值共享及池化操作等特性可有效降低網絡復雜度,提高網絡模型的魯棒性和容錯能力[23]。目前,CNN 在地球物理領域的應用已經取得一些研究成果,如:林年添等[24]將CNN 運用于地震油氣儲層的預測;奚先等[25]研究了基于CNN 的地震偏移剖面中散射體的定位和成像;何旭等[26]利用CNN 識別測井相;BI 等[27]通過CNN 對微震波形進行自動識別;PUZYREV[28]利用CNN實現了特定裝置的可控源電磁的二維反演。以上研究表明CNN 在地球物理領域具有較大應用潛力,但這些研究大部分集中在地震方面,在電磁方面運用尚少。為此,本文作者提出一種基于CNN 的大地電磁反演方法,此方法完全非線性,不會產生傳統線性反演方法出現的問題,通過給定觀測值(視電阻率、相位)就能立刻對模型參數進行反演。其步驟為:首先闡述CNN 的基本原理,介紹利用CNN 實現大地電磁反演的基本步驟,然后對理論模型進行反演測試,并討論不同輸入分量以及深度對反演效果的影響,最后通過實測結果驗證此方法的有效性。
20世紀60年代,HUBEL等[29]在研究貓腦皮層中用于方向選擇和局部敏感的神經元時,發現其獨特的網絡結構可以有效地降低反饋神經網絡的復雜性,繼而提出了CNN。近年來,隨著遷移學習理論[30]在CNN 上的成功應用,CNN 的應用領域得到了進一步擴展[31-32]。目前,CNN已成為當前最受關注的研究熱點之一[33]。
CNN 具體結構如圖1所示,由輸入層、卷積層、池化層、全連接層及輸出層組成,卷積層和池化層一般取值若干并且交替設置,即1個卷積層連接1個池化層,池化層后再連接1個卷積層。以下介紹CNN的各組成部分。

圖1 CNN示意圖Fig.1 Schematic diagram of CNN
1.1.1 卷積層
卷積層由多個特征面組成,每個特征面由多個神經元組成,它的每一個神經元通過卷積核與上一層特征面的局部區域相連。卷積層通過卷積操作提取輸入的不同特征,第一層卷積層提取低級特征如邊緣、線條、角落,更高層的卷積層提取更高級的特征[23]。卷積過程如圖2所示。前一層的特征圖與一個可學習的卷積核(卷積核實際是一個權值矩陣[34])進行卷積,再加上偏置,然后通過激活函數激活,最后得到輸出特征圖。每一個特征輸出圖可以組合多個特征卷積圖,卷積計算公式如下:

式中:為卷積層l第j個通道的凈激活輸出,它是通過對前一層輸出的特征圖進行卷積求和與偏置操作后得到的;為卷積層l第j個通道的輸出;f(·)為激活函數;Mj為特征圖子集;為卷積核;為卷積層的偏置項;“*”為卷積符號。

圖2 卷積示意圖Fig.2 Schematic diagram of convolution
1.1.2 池化層
池化層在卷積層之后,同樣由多個特征面構成,它的每一個特征面與卷積層的特征面相對應,且池化層的神經元與其輸入層的局部接受域相連,不同神經元的局部接受域不重疊。池化層通過降低特征面的分辨率來獲得空間不變性特征,降低網絡的規模,同時起到提取二次特征的作用,常用的池化方法有平均池化(mean-pooling)、最大池化(max-pooling)等[35],池化的過程可由以下公式描述:

式中:為池化層l第j個通道的凈激活輸出,通過前一層的特征圖池化加權、偏置后得到;β為池化權重系數;為池化層的偏置項;pool(·)為池化函數。
1.1.3 全連接層
在CNN 結構中,全連接層位于卷積層與池化層之后,它的每一個神經元與前一層神經元完全連接,可以整合卷積層或池化層中具有類別區分的局部信息,全連接層的計算過程如下:

式中:ul為全連接層l的凈激活輸出,由前一層的輸出加權求和偏置后得到;ωl為全連接網絡的權值系數;bl為全連接層的偏置項。
1.1.4 激活函數
激活函數在神經網絡信息傳遞過程中具有十分重要的作用,若不采用激活函數,則每一層的輸入將是上一層輸出的線性函數,網絡的逼近能力將非常有限。因此,常采用一些非線性函數作為激活函數,這樣,神經網絡就能擁有更強大的非線性映射能力。常用激活函數有sigmoid 函數、tanh函數和ReLU函數,數學形式如下:

sigmoid 函數能夠將輸入的連續實值變換為0和1之間的輸出,tanh函數能夠將變量映射到-1到1 之間。對于ReLU 函數,當輸入大于0 時,輸出等于輸入,否則輸出為0。與sigmoid 函數和tanh函數相比,ReLU 函數能夠解決梯度消失等問題,還能加快收斂速度,因此,在目前的CNN結構中,常用ReLU函數作為激活函數。
首次將CNN運用于大地電磁反演。目前CNN的主要應用集中在圖像識別、目標檢測等領域[36-37],處理的目標通常是圖像,網絡的輸入為像素點。在地球物理領域,采集的某一個觀測數據理論上可以看作1個像素點,因此,本文將大地電磁視電阻率及相位看作不同顏色通道的像素點,作為卷積神經網絡的輸入,將地電模型參數作為網絡的輸出,將CNN 運用于大地電磁反演。實現流程如圖3所示,具體實施步驟如下。
1)CNN 搭建及初始化。依據待解決問題搭建CNN 基本框架,采用小隨機數(均值為0,方差為0.01 的高斯分布)對網絡的權值與偏置進行初始賦值。
2)樣本數據集構建及輸入。通過大地電磁二維正演獲取樣本數據集,正演方法采用有限單元法[38]。為提高計算效率與精度,對網格進行雙二次插值,并將整個網格劃分為2個區域:目標區域和網格外延區域。目標區域為地質體的賦存區域,同時也是數據采集區域,以均勻網格剖分;網格外延區域的網格步長按二倍遞增。在目標區域內,地質體按一定步長進行移動,以獲取不同位置的視電阻率及相位,并記錄相應地電模型參數。采用零-均值歸一化方法對數據集進行標準化,將標準化后的視電阻率及相位作為CNN 雙通道輸入,相應模型地電參數作為輸出。

圖3 基于CNN的大地電磁反演流程圖Fig.3 Flow chart of MT inversion based on CNN
3)卷積。將步驟2)中的輸入數據進行卷積處理,即將每一個位置的輸入進行卷積變換映射成新值,將卷積核看成權重記為向量w,輸入數據記為向量x,則該位置卷積結果為y=w'x+b,即向量內積加上偏置,將x變換為y,然后通過激活函數激活即可得到卷積結果,在此過程中,不同的卷積核將實現數據不同特征的提取。
4)池化。采用最大池化法對步驟3)中卷積后的數據進行池化處理,即將池化窗口區域內數據的最大值保留,再通過激活函數激活就可得到結果,通過此步驟可對數據體進行降維及二次特征提取。根據本文待解決問題的復雜程度,對數據再次進行卷積操作,即回到步驟3)。
5)信息正向傳播和誤差反向傳遞。CNN 在學習訓練過程中,通過BP算法進行監督訓練,需要通過信息正向傳播和誤差反向傳遞2個階段,此操作在全連接層中進行。在誤差正向傳播階段,池化后的數據進入全連接層實現信息的非線性變換,到達輸出層;在誤差反向傳播階段,網絡計算出正向輸出與理論值的誤差,再計算所有網絡層的靈敏度及總誤差對所有網絡參數的偏導數,從而逐層計算出誤差,以修正網絡權值與閾值。通過不斷的信息正向傳播和誤差反向傳遞實現網絡訓練,當誤差達到期望值或迭代次數達到預定值時,網絡停止訓練,獲得最佳反演網絡排列及超參數。
6)模型反演。將學習訓練獲得的最佳反演網絡排列及超參數作為網絡的初始賦值,將未知模型響應數據標準化后直接輸入網絡,通過一次正向傳播即可得到輸出。
在二維大地電磁TM極化模式下,建立地下高阻異常體、低阻異常體、高低阻組合模型、地壘以及地塹5種基礎模型,如圖4所示。在正演過程中,有限單元法目標區域網格規模設置為35×32,網格間距為200 m。為減少計算量,控制異常體和圍巖電阻率不變,通過移動異常體的位置來獲取樣本數據集,異常體移動位置遍布整個網格,每次移動間距為400 m。計算時,選擇36個測點,測點間距為200 m,并記錄16個頻點的視電阻率及相位,頻點范圍為10-3~103Hz,按對數等間隔分布。將采集范圍內所有測點和頻點電阻率及相位作為網絡輸入,則輸入為2×36×16的三維數據矩陣;將地電模型參數作為網絡輸出,則輸出層神經元個數為35×32,即1 120個。
網絡參數的選擇影響網絡的穩定性與反演效果,通過對比試驗確定了適合本研究的參數取值,具體參數設置如下:卷積層采用雙卷積結構,第一卷積層含卷積核32 個,第二卷積層含卷積核64個,池化層均采用最大池化法,全連接層設計為2層,神經元節點數均為1 000個;設置網絡訓練次數為200輪,損失函數采用均方誤差函數,訓練過程中使用Adam 優化算法[39]進行優化,并采用RELU激活函數與節點丟棄技術防止梯度彌散和過擬合現象。
本文就大地電磁TM模式下高阻異常體、低阻異常體、高低阻組合模型、地壘和地塹進行反演測試,設定理論模型與反演結果的均方誤差為反演模型評價指標,均方誤差eMSE表達式如下:

式中:n為反演參數個數;mi為歸一化后測試模型電阻率參數理論值;為歸一化反演電阻率。

圖4 基礎地電模型Fig.4 Basic geoelectricity models
將高阻異常體和高低阻組合模型的CNN 反演結果與最小二乘反演結果進行對比,并采用高阻異常體和地壘模型討論輸入分量對反演效果影響,采用低阻異常體和地塹模型探索深度對反演效果的影響。
2.3.1 CNN與基于光滑約束的最小二乘反演對比
選擇高阻異常體和高低阻組合模型的CNN 反演結果與基于光滑約束的最小二乘反演結果進行對比,基準模型示意圖如圖4(a)和(c)所示。用作反演的高阻異常體電阻率為1 000 Ω·m,長度和寬度均為1.0 km,背景電阻率為100 Ω·m,異常體中心點在水平方向位置0 km 處,埋深3.2 km;高低阻組合模型高阻異常電阻率為1 000 Ω·m,低阻異常電阻率為100 Ω·m,長度和寬度均為1.0 km,背景電阻率為500 Ω·m,高阻異常體中心點在水平方向0 km 處,埋深為2.0 km,低阻異常體中心點在水平方向0 km處,埋深為3.8 km。
圖5(a)所示為高阻異常體CNN 反演結果,圖5(b)所示為高阻異常體最小二乘法反演結果,圖6(a)所示為高低阻組合模型CNN反演結果,圖6(b)所示為高低阻組合模型最小二乘法反演結果,其中,紅線代表真實高阻模型位置,藍線代表真實低阻體位置。從圖5和圖6可見:最小二乘法反演結果收斂程度較差,對模型體邊緣反映不夠清晰,而CNN 反演方法能較準確地反映異常體的位置、大小和電阻率,能較清晰地刻畫模型體的邊緣。
2.3.2 反演結果的正演擬斷面圖對比分析
將2.3.1 節中的單一高阻異常體和高低阻組合模型CNN 反演結果進行正演,圖7所示為高阻異常體理論正演擬斷面圖和CNN 反演結果的正演擬斷面圖,圖8所示為高低阻組合模型理論正演擬斷面圖和CNN 反演結果的正演擬斷面圖。從圖7和圖8可見:無論是反演結果的擬合視電阻率還是阻抗相位,都與其理論的正演結果基本吻合,說明了此反演方法的準確性;值得一提的是,在視電阻率對比結果中,這2種模型的反演擬合視電阻率與理論視電阻率在各個頻段均有一定差異;而在相位對比結果中,低頻部分的反演擬合相位與理論相位差異并不大,說明相位在低頻段的反應并不靈敏。

圖5 高阻異常體不同方法反演結果Fig.5 Inversion results of high-resistance anomaly using different methods

圖6 高低阻組合模型不同方法反演結果Fig.6 Inversion results of high-low resistance model using different methods

圖7 高阻異常體理論正演結果與CNN反演結果的正演結果對比Fig.7 Comparison between forward data of theoretical high-resistance anomaly and forward data of CNN inversion results
2.3.3 不同輸入分量反演結果對比
選擇高阻異常體和地壘模型討論不同輸入分量對反演效果的影響。輸入分量包括3種模式:視電阻輸入、相位輸入、視電阻率和相位聯合輸入。基準模型示意圖如圖4(a)和(d)所示。用作反演的高阻異常體電阻率為1 000 Ω·m,長度和寬度均為1.0 km,背景電阻率為100 Ω·m,異常體中心在水平方向位置0 km 處,埋深為2.4 km;地壘模型凸起部分電阻率為1 000 Ω·m,長度為1.2 km,寬度為2.0 km,蓋層電阻率為100 Ω·m,凸起部分中心在水平方向位置0 km處,埋深為2.4 km。
圖9(a)和(d)所示為高阻異常體和地壘模型視電阻率和相位聯合輸入的反演結果,圖9(b)和(e)所示為視電阻率輸入的反演結果,圖9(c)和(f)所示為相位輸入的反演結果。從圖9可以看出:3 種輸入模式均取得良好的反演成像效果,相比于單一輸入反演,電阻率和相位聯合輸入的反演結果對邊界的刻畫更加清晰,更加接近真實模型。表1所示為測試集中10 組數據不同輸入模式下的均方誤差對比,從表1可以看出:高阻異常體模型聯合反演的平均均方誤差為0.024時,視電阻率和相位反演的平均均方誤差分別為0.089和0.090;地壘模型聯合反演的平均均方誤差為0.035時,視電阻率和相位反演的平均均方誤差分別為0.067和0.076。聯合反演的均方誤差是單一反演的30%~50%,表明視電阻率和相位聯合反演具有更好的反演效果。
2.3.4 不同深度反演結果對比
選擇低阻異常體和地塹模型探索CNN 在不同深度的反演效果,基準模型示意圖如圖4(b)和(e)所示。用作反演的低阻異常體電阻率為100 Ω·m,長度和寬度均為1.0 km,背景電阻率為1 000 Ω·m;地塹模型凹陷部分電阻率為100 Ω·m,長度為1.2 km,寬度為2.0 km;基底電阻率為1 000 Ω·m。定義深度[-2,0)km 為淺部區域,[-4,-2)km 為中部區域,[-6,-4)km 為深部區域,若反演模型體中心位于淺部,則將該模型體歸為淺部模型體,以此類推。

圖8 高低阻組合模型理論正演結果與CNN反演結果的正演結果對比Fig.8 Comparison between forward data of theoretical high-low resistance model and forward data of CNN inversion results

圖9 高阻異常和地壘模型不同輸入分量反演結果Fig.9 Inversion results of different input components of high-resistance anomaly and ground barrier model
從測試集中隨機抽取2種模型在不同深度上的反演結果(由于是隨機抽取,橫向位置與深度均有不同),如圖10所示。圖10(a),(b)和(c)所示分別為低阻異常體淺部、中部和深部的反演結果,圖10(d),(e)和(f)所示分別為地塹模型淺部、中部和深部反演結果,紅線代表真實模型界面。從圖10可以看出:CNN 方法在各個深度處均取得了較好的反演效果,相比于深部,淺部反演模型的效果更好,對異常區域邊界刻畫更為清晰。表2所示為測試集中10 組不同深度反演模型體的均方誤差對比。從表2可知:淺、中、深部反演模型的平均均方誤差依次增加,說明CNN 反演方法對淺部的地質體具有更好的反演效果,隨著深度增加,反演效果變差。
利用CNN 方法對云南省云縣至鳳慶高速公路某隧道的實測電磁參數資料進行反演,采集系統為EM3D,含6個采集通道。經調查,隧道區域附近巖性為前奧陶系石英片巖和絹云母石英片巖。所選測線位于隧道正上方,共布設26 個測點,測點間距10 m,記錄41 個頻點,頻點范圍為2.5~19 983.3 Hz。由于實測數據的頻點、點距、電阻率等與理論模型試驗不同,需要重新構建樣本數據集。首先,根據實際探測區范圍構建網格,點距為10 m,所選頻點與實測數據頻點一致;再根據該地區已知的地質資料和巖石礦物測定結果,確定模型背景電阻率范圍為1 000~1 500 Ω·m,異常體電阻率范圍100~500 Ω·m,通過改變電阻率(按100 Ω·m 等量改變),同時移動異常體的位置(按20 m等間距移動)來獲取樣本數據集;然后,將樣本數據集放入已搭建的網絡(所選參數與2.2節中一致)中進行學習訓練;最后,將該地區實測的電阻率與相位作為CNN 的雙通道輸入,通過一次網絡正向傳播得到輸出結果。

表1 不同輸入分量反演參數的均方誤差對比Table1 Comparison of mean square error of inversion results of different input components

圖10 低阻異常和地塹模型在不同深度上的反演結果Fig.10 Inversion results of low-resistance anomaly and ground ridge model at different depths
圖11(a)所示為實際資料的CNN反演結果,圖11(b)所示為最小二乘反演結果,實線代表隧道位置,虛線代表鉆孔位置。2種反演結果均顯示隧道50 m 處正上方存在低阻異常區域,推測此區域可能存在含水率較高的隱伏破碎帶或溶洞,2種反演方法雖然存在一定差異,但反演出的低阻位置大致相符。鉆孔資料(ZK1)顯示高程2 017~2 050 m為第四系含礫粉質黏土;高程1 867~2 017 m為中等風化云母石英片巖,含水率較高;高程1 823~1 867 m 為強風化石英云母片巖,呈松散狀,含水率高。采用四極法測定鉆孔巖芯電阻率,高程范圍1 823~1 867 m內的電阻率明顯偏低,與反演結果顯示的低阻異常區域吻合,驗證了此方法在實際資料反演中的有效性。

表2 不同深度反演結果均方誤差對比Table2 Comparison of mean square error of inversion results of different depths

圖11 實測資料不同方法反演結果Fig.11 Inversion results of measured data by different methods
1)將卷積神經網絡引入大地電磁反演,通過數據集建立、模型構建、網絡訓練和反演測試等步驟實現了地質模型體的精準定位與成像,其“聚焦”反演效果優于最小二乘法。
2)卷積神經網絡的反演效果與輸入分量和地電模型深度有關,視電阻率和相位聯合反演結果優于單一參量反演結果,淺部模型體的反演結果比深部模型體的優,并且聯合反演的均方誤差是單一反演的30%~50%。
3)CNN 方法在實際資料處理中具有一定實用性,但是目前該方法僅能對簡單的地質體進行反演,其主要原因在于樣本數據集的泛化性不夠,泛化性增加會導致計算量呈幾何數量級增加。下一步將對提高模型的泛化性并同時減少計算量進行研究。