瞬變電磁長短時記憶網絡深度學習實時反演方法

范濤，薛國強，李萍，燕斌，鮑亮，宋金秋，任笑，李澤林

1 中煤科工西安研究院(集團)有限公司，西安 710077 2 中國科學院地質與地球物理研究所 中國科學院礦產資源研究重點實驗室，北京 100029 3 西安電子科技大學計算機科學與技術學院，西安 710071

0 引言

瞬變電磁法作為一種低成本的地球物理勘探方法，多年以來一直被用于資源與環境領域(Guo et al., 2020；Fan et al., 2020；Di et al., 2020a).瞬變電磁實測數據一維反演方法，如共軛梯度法、高斯-牛頓法和最小二乘法等線性化反演(孫懷鳳等，2019；Maurya et al., 2019；陳衛營等,2017；Li，2016)，也取得了一定的應用效果.實際上，反演問題是非線性、多極值、多參數的，線性化近似在一定程度上會丟失細節信息，尤其是進行較為精細的地質結構探查時，如果初始模型選擇不當，極易陷入局部最優解，而非全局最優解，導致線性化反演的結果往往與實際情況相差甚遠.

由于地下介質的復雜性，高維反演問題尚未妥善解決(Xue et al., 2020；Luan et al., 2020；Yin et al., 2015)，三維反演方法受限于復雜的三維正演算法，數據量巨大，占用資源過多，反演運算速度緩慢，很難真正投入實際應用(Cox et al., 2012；Oldenburg et al., 2013).因此，粒子群優化算法、差分蟻群算法、模擬退火算法、人工魚群優化方法、遺傳算法和貝葉斯方法等優化方法得到研究和應用(Di et al., 2020b；Li et al., 2021, 2019a；程久龍等，2014)，這些方法對初始模型依賴較小，且不需要計算雅克比矩陣，但仍需要大量正演計算，收斂速度較慢，算法參數的要求比較嚴格，參數選擇不當甚至會使反演出現早熟或是不收斂的情況，反演結果隨機性也較大.

神經網絡技術也被應用于瞬變電磁反演問題，出現了監督下降法(SDM)、反向傳播(BP)算法、深度置信網絡算法和剪枝貝葉斯神經網絡(PBNN)等(秦善強，2019；Jiang et al., 2016；王秀臣，2006；李創社等，2001)，無需經典反演中的正演計算過程，也無需求解雅克比矩陣，可大大提高反演速度和效率.隨著技術的發展，深度神經網絡被廣泛應用于各個領域，相比于淺層神經網絡，其具有更強的表達能力.目前，在電磁反演領域內，全連接網絡、卷積神經網絡與遺傳神經網絡是常用的方法(Wu et al., 2021；程久龍等，2020；Li et al., 2019b；王鶴等，2018).然而，目前該研究領域內還存在著數據集構造較為簡單，地層情況豐富度較低，模型設計簡單對數據時序性利用不足等問題.同時，深度神經網絡更容易出現過擬合、收斂速度慢、陷入局部最優等情況.針對上述難點，需要對深度神經網絡進行更好的設計，如使用各種基于動量的優化器、加入dropout、使用Layer Normalization機制(Ba et al., 2016)等.

由瞬變電磁數據可以得出其反演問題屬于時間序列問題，可使用深度學習中專門解決時序問題的循環神經網絡進行處理.長短時記憶網絡(Long Short-Term Memory Network，LSTM)是循環神經網絡(Recurrent Neural Network，RNN)中的一種特殊模型，它具備RNN的遞歸屬性(Visin et al., 2015)，同時也是RNN的一種改進模型.它具有獨特的記憶、遺忘模式，可以靈活地適應網絡學習任務的時序特征(Hochreiter and Schmidhuber, 1997).LSTM解決了RNN在長序列訓練過程中的梯度消失和梯度爆炸問題，能夠充分利用歷史信息，建模信號的時間依賴關系(Chung et al., 2014).至今，LSTM已經在各種序列任務取得了巨大成功，例如，機器翻譯、語音識別、文本分類等(Ma et al., 2020；Merity et al., 2017；Zhang et al., 2016；Cho et al., 2014).

本文開展基于LSTM的瞬變電磁反演算法研究，在提高反演精度的前提下進一步實現實時成像.首先介紹本文采用的LSTM方法理論和在此基礎上實現的瞬變電磁反演算法，然后將本文算法應用于一維和三維數值模擬數據的反演中以檢驗本文算法的有效性，并對比本文方法與傳統反演方法的計算耗時.

1 長短時記憶網絡(LSTM)

LSTM網絡由輸入層、輸出層和隱藏層構成，單元內部結構如圖1所示.

圖1 LSTM單元內部結構Fig.1 Internal structure of the LSTM unit

LSTM單元的工作流程如下：每一個時刻，LSTM單元通過3個門接受當前狀態xt和上一時刻LSTM的隱藏狀態ht-1這兩類外部信息的輸入.此外，每一個門還接收一個內部信息輸入，即記憶單元的狀態Ct-1.接收輸入信息后，每一個門將對不同來源的輸入進行運算，并且由其邏輯單元決定其是否激活.輸入門的輸入經過非線性函數的變換后與遺忘門處理過的記憶單元狀態進行疊加，形成新的記憶單元Ct.最終，記憶單元狀態Ct通過非線性函數的運算和輸出門的工臺控制形成LSTM單元的輸出ht.

各變量之間的計算公式見(1)—(5)式：

it=σ(Wxfxt+Whfht-1+bf),

(1)

ft=σ(Wxfxt+Whfht-1+bf),

(2)

Ct=ftCt-1+ittanh(Wxcxc+Whfht-1+bc),(3)

ot=σ(Wxoxt+Whoht-1+bo),

(4)

ht=ottanh(Ct),

(5)

式中：Wxc、Wxi、Wxf、Wxo為連接輸入信號xt的權重矩陣；Whc、Whi、Whf、Who為連接隱含輸出信號ht的權重矩陣；Wci、Wcf、Wco為連接神經元激活函數輸出矢量Ct和門函數的對角矩陣；bi、bc、bf、bo為偏置向量；σ為激活函數，通常為tanh或sigmoid函數.

算法的基本思路為：先將LSTM網絡按照時間順序展開為一個深層的網絡，然后使用經典的誤差反向傳播(back propagation，BP)算法對展開后的網絡進行訓練，其示意圖如圖2所示.

圖2 LSTM網絡重復模式Fig.2 LSTM network duplication mode

2 序列到序列(Seq2seq)模型

2.1 Seq2seq模型簡述

序列到序列(Seq2seq)模型為自然語言處理(natural language processing，NLP)領域的經典模型，其中最常用的Encoder-Decoder模型如圖3所示.encoder接收每個時間步的輸入信息(X1,X2,X3,X4)，通過encoder將其編碼為一個表示向量，稱為context vector(圖中用C表示)，然后交給decoder來進行解碼.decoder根據context vector與上個時間步的輸出信息產生每個時間步的輸出(Y1,Y2,〈EOS〉)，其中，〈g〉代表輸入序列的結束與輸出序列的開始，〈EOS〉代表輸出序列的結束，是一個事先定義在輸出字典里的標志，循環網絡在產生此標志后，不再進行下一時間步的操作.

圖3 Encoder-Decoder結構(a)Encoder結構(b) Decoder結構Fig.3 Encoder-Decoder structure (a) Encoder structure (b) Decoder structure

2.2 free running與teacher forcing

在NLP領域的機器翻譯問題中，decoder部分每一時間步的輸出代表著詞袋中每一個詞的預測出現概率，通常會將概率最大的詞作為翻譯的結果，并將相應的詞嵌入向量作為下一時間步的輸入之一，這種結構通常被稱為free running，此時Seq2seq內部結構如圖4所示.

圖4 Free running機制圖Fig.4 The free running mechanism diagram

為了提升訓練速度，常常使用teacher forcing(Williams and Zipser, 1989)機制，在decoder端，使用真實的目標文本作為每一時間步的輸入.每一步根據當前正確的輸出詞和上一步的隱狀態，來預測下一步的輸出詞，Seq2seq的內部結構如圖5所示.

圖5 Teacher forcing機制圖Fig.5 The teacher forcing mechanism diagram

2.3 Attention機制

隨著輸入序列的不斷增長，使用上述方法的序列到序列模型的表現會越來越差，于是人們提出了注意力(attention)機制，此方法使decoder的每個時間步的輸出可以根據不同的context vector生成，在decoder的不同時間步，產生對輸入序列的不同關注重點.

假設共有n個時間步的輸入，encoder第i個時間步的輸出為ei，decoder第j-1個時間步的輸出為qj-1，則decoder第j個時間步的context vectorCj的計算公式如式(6)到式(8)所示.

(6)

(7)

(8)

Attention機制模型如圖6所示.decoder每個時間步使用不同的Cj.

圖6 Attention機制圖Fig.6 The attention mechanism diagram

3 基于LSTM的瞬變電磁反演

在模型訓練環節，將正演得到視電阻率數據輸入模型進行訓練，并通過超參數優化提升模型性能，之后運用評估指標進行模型評估.最終，得到訓練、優化完成的模型用于預測未知的地層結構，即使用模型完成瞬變電磁實時反演成像.

3.1 模擬數據準備與評價指標

模擬數據主要用來在理論層面上展示算法的可用性.而由于實際地層構成情況十分復雜，則需要根據實際面臨的反演問題對模擬生成的訓練數據進行調整，并隨后在實測數據上進行反演測試.

本文使用基于漢克爾變換和余弦變換的瞬變電磁一維正演程序生成樣本集.在模擬數據實驗中，將正演模擬數據分為train(訓練集)、eval(驗證集)、test(測試集)三部分.在檢測數據實驗中，使用人工指定參數、通過三維數值模擬生成的模型數據作為額外的測試數據.檢測數據主要用來驗證反演算法對復雜模型的適用性.

通過三個指標評價算法的效果，分別是正演曲線點相對誤差(FPRE)、反演折線面積相對誤差(BARE)和反演折線點相對誤差(BPRE).

FPRE是通過在目標地層結構與預測地層結構分別正演出的曲線上相應取出N個點，計算它們的平均相對誤差得出的，它衡量了兩個正演曲線的擬合程度.具體計算公式如式(9)，其中pi為預測地層結構的正演曲線上的第i個點，ti為目標地層結構的正演曲線上的第i個點：

(9)

BARE是通過計算代表目標地層結構與預測地層結構的折線間的相對面積誤差得出的，它通過面積的差值比衡量了兩個地層結構折線的相似程度.具體計算公式如式(10)，其中Stp指的是兩個折線之間圍城的面積，St指的是目標地層結構折線與x軸(x軸為深度，y軸為電阻率)圍成的面積：

(10)

BPRE是通過在代表目標地層結構與預測地層結構的折線上相應取出N個點，計算它們的平均相對誤差得出的，它通過點之間的差值比衡量了兩個地層結構折線的相似程度.具體計算公式如式(11)，其中pi為預測地層結構折線上的第i個點，ti為目標地層結構折線上的第i個點：

(11)

3.2 Seq2seq網絡設計

反演算法整體框架使用Seq2seq模型，包含encoder，decoder，attention三部分，整體結構如圖7所示.

Encoder部分的主要目的是將input的信息用一個向量表示出來，本文所使用的encoder為雙向LSTM網絡.在進行地層預測時，需要綜合各個地層的信息，最終的地層預測結果由前面若干輸入和后面若干輸入共同決定，因此使用雙向LSTM網絡以使得預測結果更加準確.

Encoder的詳細結構如圖8所示.

圖7 反演算法整體框架Fig.7 Overall framework of the inversion algorithm

圖8 Encoder架構圖Fig.8 Encoder Architecture diagram

早期的Seq2seq模型沒有使用attention機制，在輸入序列input較長或者包含信息較多時，會丟失很多信息，attention機制可使decoder在翻譯時能夠依靠包含不同input信息側重點的context vector.本文使用Bahdanau Attention(Bahdanau et al.，2016)，主要計算過程如下：

(1)第i個target上下文向量ci會根據input的隱藏狀態加權求和得到：

(12)

(2)其中，對于每個hj的系數αij根據下式得到：

(13)

其中eij=a(si-1,hj)，si-1為輸出序列的上一個隱狀態，a將decoder的上一個隱狀態si-1和encoder的隱狀態hj作為輸入，計算出一個xi,yj對齊的值eij(此值可以衡量輸出序列和輸入序列的對應情況，越對應，權重越高，即注意力高)，再歸一化得到權重.

Attention機制的架構如圖9所示.

圖9 Attention架構圖Fig.9 The attention architecture diagram

Decoder部分每個時間步的輸入中不包括上一時間步的輸出，即本文的算法在free running結構的基礎上進行了適應性修改.此做法是為了緩解過擬合，增強模型的泛化能力.在機器翻譯問題中，decoder每一時間步的輸出是詞的選中概率，算法需要根據此概率選出每一時間步的真正輸出詞匯，然后將其相應的詞嵌入向量作為下一時間步的輸入，從而為模型補充被選詞的信息.而針對瞬變電磁反演問題，decoder的每一個時間步的輸出是實數，它們是LSTM的輸出htot經過后續計算得到的，同時ht也作為隱藏狀態傳入到了下一時間步中，所以decoder的每一個時間步中不需要重復地輸入信息.這樣做一方面減少了網絡的參數量，另一方面避免了從輸入項中再次引入誤差.

Decoder的架構如圖10所示.

圖10 Decoder架構圖Fig.10 The decoder architecture diagram

3.3 損失函數設計

損失函數包含三部分，分別為：電阻率損失、厚度損失和深度損失，公式如式(14)所示.電阻率損失、厚度損失分別針對電阻率、地層厚度這兩個預測目標進行擬合.同時，為了避免每個地層的厚度預測誤差累積導致較深層地層的位置發生較大偏移，加入深度損失，對近地層的厚度預測進行校正：

loss=lossρ+losst+lossd.

(14)

(15)

(16)

(17)

d1=t1,

(18)

di=ti+di-1,i=2,…,L,

(19)

(20)

(21)

3.4 超參數設置

本文所使用的Seq2seq模型的超參數以及訓練過程的超參數如表1所示.

由于針對較復雜的地層，模型訓練需要較多的數據，為了平衡正演樣本生成時間成本和模型預測的準確性，設定訓練、驗證和測試樣本的比例為 0.9∶0.05∶0.05.

表1 超參數設置表Table 1 Over the parameter setting table

4 隨機模擬數據訓練和測試

4.1 三層模型

本節通過使用表2所示的參數生成3個地層的正演模擬數據.

表2 三層模型參數設置表Table 2 The parameter setting table of the three-layer model

本節實驗最終的loss函數下降趨勢如圖11所示.從圖中可看出，loss不斷下降至收斂，且無過擬合現象.這說明該模型對三層地層的數據具備較強的擬合能力.

圖11 三層模型loss下降趨勢曲線圖Fig.11 The loss decay curve of the three-layer model

隨機選取了9個測試數據，本模型對地層分布的反演結果及目標結果圖12所示.由于學習率衰減方式為階梯式衰減，如表1所示，所以loss下降趨勢呈階梯式下降，符合預期.從圖中可看出，反演結果中絕大部分的曲線都完全吻合，少部分地層的厚度和電阻率大小有一些誤差，但均處于可接受范圍內.因此，基于時間序列的反演模型具有學習出三層地層形態的能力.

模型的三個衡量指標結果如表3所示.

表3 三層模型衡量指標值表Table 3 Table of three-layer model measures

4.2 五層模型

本節通過使用表4所示的參數生成5個地層的正演模擬數據.

表4 五層模型參數設置Table 4 Five-layer model parameter settings

圖12 三層模型部分反演結果與目標對比(a) 模型17; (b) 模型24; (c) 模型75; (d) 模型94; (e) 模型131; (f) 模型183; (g) 模型206; (h) 模型277; (i) 模型349.Fig.12 The partial inversion results of the three-layer model are compared with the target(a) Model 17; (b) Model 24; (c) Model 75; (d) Model 94; (e) Model 131; (f) Model 183; (g) Model 206; (h) Model 277; (i) Model 349.

本實驗最終的loss函數下降趨勢如圖13所示.從圖中可看出，loss不斷下降至收斂，僅出現了輕微的過擬合現象.這說明該模型對五層地層的數據仍具備較強的擬合能力.

圖13 五層模型loss下降趨勢曲線圖Fig.13 The loss downward trend curve of the five-layer model

同樣隨機選取了9個測試數據，本模型對地層分布的反演結果及目標結果圖14所示.從圖中可看出，隨著地層數的增加，雖然五層的反演結果相較于三層的結果準確度有所下降，但是Seq2seq模型的反演結果依然能夠較好地預測出地層厚度和電阻率的整體變化趨勢，基本能夠和目標地層情況吻合，說明該模型在五層地層數據上依然比較穩定.且從圖14中可以看出，深度的反演誤差相對于厚度和電阻率的反演誤差較小，從而印證了loss函數中深度loss的設計使得模型在淺層地層厚度有預測偏差時，對后續深層地層的預測不會產生較大的影響.

模型的三個衡量指標結果如表5所示.從三個評價指標中可看出，與三層數據的訓練結果相比較，五層數據的測試效果出現下降，這說明面對復雜的地層分布情況，模型需要擴充容量或進行針對性設計.

圖14 五層模型部分反演結果與目標對比(a) 模型59; (b) 模型87; (c) 模型165; (d) 模型216; (e) 模型269; (f) 模型323; (g) 模型485; (h) 模型540; (i) 模型634.Fig.14 The partial inversion results of the five-layer model are compared with the target(a) Model 59; (b) Model 87; (c) Model 165; (d) Model 216; (e) Model 269; (f) Model 323; (g) Model 485; (h) Model 540; (i) Model 634.

表5 五層模型衡量指標值Table 5 Measurement values of the five-layer model

5 檢測數據驗證

5.1 均勻空間單一低阻異常模型

首先設置一個簡單地電模型檢驗反演模型的準確性，模型設計如圖15a所示，背景為均勻半空間，電阻率為500 Ωm，異常體在模型中心，電阻率為10 Ωm，異常體埋深為50 m，厚度為50 m，發射線框邊長100 m，測點共布置3個，異常體上方1個，兩側各1個.

對圖15a中的測點數據進行反演，可得如圖15(b、c)所示的反演結果.對比(b)(c)兩圖可以看出，OCAMM反演結果盡管能夠反演出測區中部有一低阻異常區，橫向邊界也很清晰，但深度方向的邊界分辨率較低，上邊界大約在30 m，下邊界無法分辨，異常中心的電阻率值也要明顯高于設置值，且異常體兩側的正常測點也表現出了輕微的高阻-低阻-高阻分層特征，與設置的均勻背景不符；而Seq2seq模型反演結果中，低阻異常體的邊界非常清晰，與設置位置完全吻合，異常體的電阻率值與設置值相同，異常體兩側的正常測點的電阻率也與設置的均勻背景值一致，僅異常體下方的電阻率值比設置的背景值略低.

圖15 均勻空間模型示意圖和反演結果(a) 模型示意圖; (b) OCAMM反演結果; (c) Seq2seq模型反演結果.Fig.15 Model schematic diagram and inversion results(a) Model schematic; (b) OCAMM inversion results; (c) Seq2seq model inversion results.

本次測試，OCAMM反演結果用時約16.34 min，Seq2seq模型反演用時小于1 s.

以該模型為基礎，通過對正演數據加入不同比例的噪聲，檢測Seq2seq模型反演算法的適用性，結果如圖16所示.圖中可以看出，即使噪聲增加到20%，對與背景值差異較大的異常體邊界的反演結果依然較為準確，但噪聲僅增加2%，均勻背景中就開始出現不均勻的現象，噪聲增加至10%，均勻背景就表現出明顯Q型地層特征，且異常體上部地層電阻率也與兩側背景值出現了肉眼可見的差異，在噪聲進一步增加后，這種表現較為穩定.因此綜合來看，當正演數據中噪聲占比達10%以上時，反演成果中與背景值差異較大的顯著異常仍有較大的可信度，但與背景值差異較小的異常的置信度需要調低.

圖16 加噪聲正演數據的反演結果對比(a) 無噪聲; (b) 噪聲占比2%; (c) 噪聲占比5%; (d) 噪聲占比10%; (e) 噪聲占比15%; (f) 噪聲占比20%.Fig.16 Comparison of inversion results for noisy forward modeling data(a) No noise; (b) Noise accounts for 2%; (c) Noise accounts for 5%; (d) Noise accounts for 10%; (e) Noise accounts for 15%; (f) Noise accounts for 20%.

5.2 復雜地層高低阻異常交錯模型

為進一步驗證Seq2seq模型反演算法對復雜地電模型的反演效果，設計了如圖17a的模型.圖中背景電阻率為300 Ωm，中間含斷層的煤層層位電阻率為700 Ωm，厚度為200 m，左側寬度為11個測點，右側寬度為17個測點，左上角的低阻異常電阻率為10 Ωm，埋深50 m，厚度100 m，寬度為4個測點，左下角的高阻異常電阻率為1000 Ωm，埋深500 m，厚度200 m，寬度為2個測點，右上角的高阻異常電阻率為1000 Ωm，埋深50 m，厚度100 m，寬度為2個測點，右下角的低阻異常電阻率為10 Ωm，埋深400 m，厚度300 m，寬度為4個測點.發射線框邊長400 m.

圖17 復雜地層模型示意圖和反演結果(a) 模型示意圖; (b) OCAMM反演結果; (c) Seq2seq模型反演結果.Fig.17 Model schematic diagram and inversion results(a) Model schematic; (b) OCAMM inversion results; (c) Seq2seq model inversion results.

對圖17a中的測點數據進行反演，可得如圖17b、c所示的反演結果.對比(b)(c)兩圖可以看出，OCAMM反演結果對煤層、斷層和4個異常體均有一定反映，但埋深、厚度和邊界均不準確、不清晰，尤其是異常埋深和下邊界，均與設置值不符，下部兩個異常體的下邊界甚至無法反映，各個部分的電阻率也與設置值差異較大；而Seq2seq模型反演結果中，煤層、斷層和4個異常體的邊界反映均較為清晰，尤其是埋深、厚度、上下界面都與設置值基本吻合，僅煤層厚度存在一定誤差，而較寬的兩個異常體，中間和兩側的厚度及上下邊界也有一定誤差，此外，除低阻異常體外的其他區域，電阻率值較設置值均偏低，但總體來看，反演精度顯著高于OCAMM反演結果，與預設模型吻合度較高.

本次測試，OCAMM反演結果用時約30.58 min，Seq2seq模型反演用時小于1 s.

6 結論

本文基于長短時記憶網絡(LSTM)結合瞬變電磁正演算法，成功實現了可實時成像的高精度瞬變電磁反演.數萬組以上隨機正演模型的反演結果表明，Seq2seq模型的反演算法計算精度較高，但隨著地層層數的增加，會出現一定的過擬合現象，需要擴充樣本庫容量或進行針對性正演模型設計.三維數值模擬模型的反演結果表明，Seq2seq模型的反演算法精度和速度都顯著優于傳統的OCCAM反演方法，尤其是對于復雜模型，Seq2seq模型的反演結果能準確反映多個異常體和地層的邊界位置、電阻率差異，反演速度更是達到實時化水平.本文的反演結果已經證實了基于LSTM網絡的Seq2seq模型相較于傳統反演方法的優越性，但是瞬變電磁三維反演仍是業界一個非常具有挑戰性的問題，基于Seq2seq模型的三維三分量反演算法將是我們未來的重點研究方向.