楊凱，唐衛(wèi)東，劉誠，賀景龍，姚川

(1.中國地質調查局 西安礦產資源調查中心，陜西 西安 710000； 2.中國地質大學(武漢) 地球物理與空間信息學院，湖北 武漢 430074)

0 引言

大地電磁測深是一種以天然電磁場為場源的地球物理勘探方法，其觀測頻帶較寬，有效信號微弱，隨著人文噪聲的日益增加，大地電磁的“凈土”已經(jīng)基本消失。其中，由測區(qū)附近用電設備的充放電所引起的方波噪聲是一種常見的高強度噪聲，主要集中在電道，可產生高于正常信號十幾倍到幾十倍的干擾[1-2]，這類噪聲影響幅度大，影響頻點多，極大干擾了阻抗估計的結果。為了壓制方波干擾，許多學者做了大量研究:嚴家斌[3]提出了基于小波變換的迭代回歸噪聲改正方法，對含方波干擾的脈沖類噪聲進行了壓制，有效改善了數(shù)據(jù)質量；蔡建華等[4]將經(jīng)驗模態(tài)分解(EMD)應用于方波噪聲處理，突出了有用信號；湯井田、李晉等[5]提出了一種基于數(shù)學形態(tài)濾波的大地電磁去噪方法，對安徽廬樅礦集區(qū)的方波、脈沖、三角波等噪聲進行了處理，有效抑制了大尺度干擾和基線漂移；湯井田、劉祥等[6]討論了不同仿真方波噪聲對測深曲線的影響及遠參考法對其的去噪效果，研究表明遠參考在一定條件下可以消除方波干擾；王輝、魏文博等[7]利用同步大地電磁時間序列信號之間的關系，用合成的無噪數(shù)據(jù)段代替含噪數(shù)據(jù)段，成功去除了大于窗口長度的方波噪聲，精度較高，但該方法需要一段無明顯干擾的遠參考數(shù)據(jù)；湯井田、李廣等[8]通過字典學習提取人文干擾特征，利用構建的冗余字典分離了AMT數(shù)據(jù)中的仿真方波噪聲；李晉、張賢等[9]利用變分模態(tài)分解(VMD)和匹配追蹤(MP)對模擬方波噪聲進行了處理，明顯改善了數(shù)據(jù)質量。

近年來，深度學習已經(jīng)在地球物理部分方法的應用中取得了不錯的效果[10]，不少學者將其成功引入電磁、地震數(shù)據(jù)處理以及重磁反演等地球物理領域[11-15]。深度學習方法適用范圍廣，泛化能力強，在大地電磁方面的應用主要集中在時間序列處理上，Manoj和Nagarajan[16]最早提出了利用人工神經(jīng)網(wǎng)絡自動執(zhí)行大地電磁時間序列編輯的方法，該方法提高了工作效率減少了人工編輯的主觀因素。本文所應用的長短時記憶網(wǎng)絡(long short-term memory, LSTM)是一種特殊的循環(huán)神經(jīng)網(wǎng)絡，該神經(jīng)網(wǎng)絡在語音識別、自然語言處理以及時間序列相關的領域有較為廣泛的應用，近年來也被逐漸引入到地球物理領域，例如:許滔滔等[17]將LSTM網(wǎng)絡應用于工頻干擾壓制，有效去除了工頻干擾;王斯昊等[18]用LSTM網(wǎng)絡去除了時間序列的階躍信號;汪凱翔等[19]利用LSTM網(wǎng)絡對地電場數(shù)據(jù)進行了處理，去除了測試集中不同種類的噪聲。以上方法都是以時間序列的低頻、大尺度噪聲或者信號輪廓為網(wǎng)絡輸出值。與前人不同的是，本文使用大地電磁信號本身作為網(wǎng)絡輸出，通過選取標準大地電磁時間序列隨機添加仿真方波噪聲作為網(wǎng)絡訓練輸入，以無噪大地電磁時間序列作為網(wǎng)絡目標輸出，讓網(wǎng)絡存儲和學習大地電磁信號本身的特征，從含噪時間序列中自動提取符合大地電磁信號特征的序列，從而實現(xiàn)抑制方波噪聲的目的。

1 LSTM網(wǎng)絡簡介

長短時記憶網(wǎng)絡(LSTM)是一種特殊的循環(huán)神經(jīng)網(wǎng)絡，由Hochreiter和Schmidhuber于1997年提出[20]，主要是為了解決長序列訓練過程中的梯度消失和梯度爆炸問題，被廣泛應用于時間序列相關領域。

一個標準的LSTM網(wǎng)絡神經(jīng)元包含了3個門(輸入門、輸出門、遺忘門)和1個記憶細胞(圖1)，可以歸納為以下幾式：

圖1 LSTM網(wǎng)絡神經(jīng)元基本結構

ft=σ(Wf·[ht-1,xt]+bf)，

(1)

it=σ(Wi·[ht-1,xt]+bi)，

(2)

ot=σ(Wo·[ht-1,xt]+bo)，

(3)

gt=tanh(Wg·[ht-1,xt]+bg)，

(4)

ct=gt⊙it+ct-1⊙ft,

(5)

ht=tanh(ct)⊙ot,

(6)

式中：ht-1為t-1時刻即上一個神經(jīng)元的隱藏層；xt為t時刻的特征向量；σ為sigmoid激活函數(shù)；it、ft、ot分別為t時刻輸入門、遺忘門和輸出門的狀態(tài)；gt為t時刻記憶細胞的候選值；ct為t時刻記憶細胞的狀態(tài)，也作為t+1時刻即下一個神經(jīng)元的初始記憶細胞；ht為t時刻的隱藏層，也為t+1時刻即下一個神經(jīng)元的初始隱藏層；Wi、Wf、Wo、Wg為傳播權重矩陣；bi、bf、bo、bg為偏置向量；⊙表示矩陣元素對應相乘。

總的來說，單向LSTM神經(jīng)元內部有3個處理步驟：

第一步為選擇忘記階段，主要是對上一個節(jié)點傳進來的輸入進行選擇性忘記。根據(jù)上一個神經(jīng)元所傳遞的隱藏層ht-1和本神經(jīng)元的輸入xt生成遺忘門ft，來控制上一個神經(jīng)元所傳遞的記憶細胞ct-1中哪些需要留，哪些需要“忘”(式(1)、式(5))。

第二步為選擇記憶階段，主要是通過上一個神經(jīng)元的隱藏層和本次的輸入xt生成輸入門it，再與t時刻記憶細胞的候選值gt作用來決定需要記住的信息(式(2)、式(4)、式(5))。

第三部為輸出階段。通過前兩個階段的“遺忘”和“記憶”，共同決定了記憶細胞ct的最新狀態(tài)，而后ct通過激活函數(shù)tanh的放縮，由輸出門ot控制生成新的隱藏層狀態(tài)ht(式(4)、式(5)、式(6))。

本文使用雙向LSTM網(wǎng)絡。雙向LSTM網(wǎng)絡就是在序列正向處理的基礎上將序列逆向再處理一次，這樣神經(jīng)元不僅能獲取“過去”時刻的序列信息，也能獲取“未來”時刻的信息，能更好地記錄其上下文的關系，從而取得更佳的學習效果。

2 LSTM網(wǎng)絡搭建

2.1 網(wǎng)絡設計

本網(wǎng)絡主要由輸入層、隱藏層和全連接層構成，通過實際測試和查閱文獻將隱藏層層數(shù)設置為2，每層128個節(jié)點，全連接層輸出維度設為1(如圖2所示)，損失函數(shù)使用均方誤差損失函數(shù)(MSEloss)，優(yōu)化器使用適應性矩估計優(yōu)化器(adaptive moment estimation，Adam)，該優(yōu)化器具有收斂速度快、調參方便等優(yōu)點，適合解決含大規(guī)模數(shù)據(jù)和參數(shù)的優(yōu)化問題[21]。需要注意的是，網(wǎng)絡結構的設計和參數(shù)的選擇帶有一定的經(jīng)驗性和主觀性，需要在參考前人結果的基礎上不斷調整試驗，針對需要解決的問題選擇最合適的組合。

圖2 LSTM網(wǎng)絡結構示意

2.2 網(wǎng)絡訓練

網(wǎng)絡訓練主要包括數(shù)據(jù)集構建、數(shù)據(jù)集歸一化(或標準化)、網(wǎng)絡訓練和驗證等幾個步驟。

2.2.1 數(shù)據(jù)集構建

選取2020年10月12日在陜西省寧強縣某地用MTU-5A大地電磁儀采集的、無明顯人文干擾且阻抗估計穩(wěn)定的數(shù)據(jù)段Ex通道序列，將其視為無噪聲的大地電磁時間序列，長度2 400個點，采樣率2 400 Hz(圖3)。首先，隨機生成100組方波信號，因為待處理的實測方波噪聲主頻在24 Hz左右，故模擬方波噪聲頻率隨機分布在23～25 Hz之間，振幅在該段無噪時間序列最大值與最小值之差的0.1～8倍之間隨機取值，相位隨機。而后，將所有模擬方波噪聲各自疊加在無噪大地電磁時間序列上，合成100組仿真含噪信號，并在每組仿真信號中隨機截取64組長度 1 200個點的信號，共產生6 400組信號作為訓練集；選取與訓練集Ex對應的Ey通道數(shù)據(jù)，長度2 400個點，采樣率2 400 Hz(圖4)，將該段數(shù)據(jù)采用與合成訓練集同樣的方式生成20組仿真含方波噪聲時間序列，而后每組信號隨機截取64組長度為1 200個點的信號，共1 280組信號作為驗證集，以驗證網(wǎng)絡的實際處理能力。

圖3 訓練集原始時間序列

圖4 驗證集原始時間序列

2.2.2 數(shù)據(jù)歸一化

數(shù)據(jù)集構建好后還要進行數(shù)據(jù)歸一化(或標準化)，歸一化后數(shù)據(jù)可以提高網(wǎng)絡的收斂速度和網(wǎng)絡精度，根據(jù)大地電磁時間序列數(shù)據(jù)的特征，將其歸一化至1～-1之間。

對于網(wǎng)絡輸入數(shù)據(jù)：

xnorm=(x-xmean)/xmax,

(7)

式中：x為添加模擬噪聲的大地電磁時間序列數(shù)據(jù)，xmean為其平均值，xmax為其最大值，xnorm為歸一化值。

對于網(wǎng)絡的目標輸出數(shù)據(jù)：

ynorm=y/xmax,

(8)

式中：y為不含噪聲的大地電磁時間序列數(shù)據(jù)，xmax為含噪數(shù)據(jù)的最大值，ynorm為歸一化值。

對所有輸入網(wǎng)絡的數(shù)據(jù)都要進行歸一化，包括訓練集、驗證集，由于網(wǎng)絡的理想輸出為歸一化值，輸出后還要進行反歸一化：

ypred=xmaxynorm_pred,

(9)

式中：ypred為反歸一化值，即為實際無噪大地電磁數(shù)據(jù)；ynorm_pred為網(wǎng)絡輸出的歸一化預測值；xmax為網(wǎng)絡輸入的含噪大地電磁信號的最大值。

2.2.3 網(wǎng)絡的訓練和驗證

將bitch_size(一次訓練所選取的樣本數(shù))設為32，epoch(使用訓練集中的全部樣本訓練一次即為1個epoch)數(shù)設置為3 000，以保證網(wǎng)絡的收斂。

評價網(wǎng)絡訓練效果的主要有網(wǎng)絡損失曲線和網(wǎng)絡精度曲線。本文選取網(wǎng)絡理想輸出時間序列和實際輸出時間序列的歸一化互相關系數(shù)(normalized cross correlation,NCC)作為檢驗網(wǎng)絡精度的參數(shù)，具體計算如下：

(10)

式中：N表示數(shù)據(jù)長度；f(n) 和g(n)為2組離散序列；NCC值在-1～1之間，-1代表2組序列相位相反，0代表2組序列正交，1代表2組序列完全相同，NCC越接近于1說明2組序列相似性越高。

由LSTM網(wǎng)絡的學習曲線(圖5)可知，前500次epoch損失曲線急劇下降，精度曲線快速上升；訓練500次以后損失曲線緩慢下降，精度曲線緩慢上升，訓練集損失(train_loss)略小于驗證集(validation_loss)，訓練集精度(train_NCC)略大于驗證集(validation_NCC)。當訓練至1 500次epoch時，訓練集的平均NCC可達0.999 98，驗證集的平均NCC可達0.999 03，說明網(wǎng)絡很好地從含方波噪聲的序列中提取出了有效大地電磁時間序列。訓練1 500次以后網(wǎng)絡趨于穩(wěn)定，驗證集損失幾乎不再減小，精度不再明顯增加，故取epoch為1 500的網(wǎng)絡為最終模型，用來進行下一步去噪測試。

圖5 LSTM網(wǎng)絡學習曲線

3 去噪測試

3.1 仿真含噪信號去噪

選取陜西省寧強縣某地2020年10月16日使用MTU-5A大地電磁儀所采集的無明顯人文干擾且阻抗估計穩(wěn)定的Ex通道時間序列數(shù)據(jù)段作為測試信號，數(shù)據(jù)長24 000個點，采樣率2 400 Hz，由于數(shù)據(jù)過長故截取其中1 s(2 400個點)進行展示(圖6)。給測試信號疊加一主頻24 Hz，振幅為該段測試信號最大值與最小值之差1.5倍的方波噪聲作為仿真含噪信號(圖6c)，用訓練好的LSTM網(wǎng)絡進行去噪試驗。由于天然大地電磁信號的復雜性，在實際使用LSTM網(wǎng)絡提取一次大地電磁信號后，殘余噪聲里還含有部分低頻有效信號，可以使用網(wǎng)絡進行多次提取，一定程度上可以減小信號損失，但同時也會引入更多的噪聲，應該具體情況具體分析。本文所有去噪均為網(wǎng)絡提取一次的結果。

由圖6a、b可以看出原始時間序列片段無明顯人文干擾，頻譜均勻分布，無明顯干擾頻段，疊加仿真噪聲后時間序列基本被方波所淹沒，幾乎無法分辨出大地電磁信號(圖6c)，頻譜呈現(xiàn)出明顯的方波各級諧波頻譜特征，且離主頻越遠幅值越小(圖6d)。經(jīng)LSTM網(wǎng)絡去噪后的時間序列與原始時間序列的NCC高達0.971 8，較好地還原了信號，頻譜也恢復了原始信號的特征，達到了去噪的目的(圖6e、f)。

使用最小二乘法進行阻抗估計。仿真信號未去噪前阻抗結果極不穩(wěn)定，電阻率曲線在265 Hz、115 Hz和22.5 Hz頻點有較大偏離(圖7a)，阻抗相位曲線在229 Hz、159 Hz、115 Hz和22.5 Hz頻點處有較大偏離(圖7b)；這些頻點基本都在方波噪聲的各次諧波頻率附近，且隨著頻率升高干擾逐漸減小。LSTM網(wǎng)絡去噪后，幾處方波噪聲影響的頻點數(shù)據(jù)均得到了較好改善，除了離噪聲主頻較近的22.5 Hz頻點附近有能量損失外，其余頻點幾乎與原始數(shù)據(jù)重合(圖7a、b)，說明了本文方法可以有效抑制方波噪聲并改善阻抗估計質量。

圖7 仿真信號去噪前后電阻率、阻抗結果對比

3.2 實測含噪信號去噪

選取陜西省旬陽市某地2021年4月27日使用MTU-5A大地電磁儀采集的實測含方波噪聲Ex通道數(shù)據(jù)段進行實際去噪試驗。該數(shù)據(jù)段長2 400個點，采樣率2 400 Hz，由于受附近用電設備的干擾，該段數(shù)據(jù)表現(xiàn)出典型的含方波噪聲大地電磁時間序列特征，噪聲已基本將原始信號淹沒，頻譜被方波各次諧波特征所占據(jù)(圖8a、b)。經(jīng)LSTM網(wǎng)絡去噪后，時間序列更加接近大地電磁原始信號的時頻特征，除高頻和低頻有殘存的方波干擾外，大部分方波干擾被有效抑制(圖8c、d)。

圖8 實測信號去噪前后對比

使用最小二乘法進行阻抗估計。未處理的原始實測信號電阻率及阻抗相位多處頻點偏移較大，特別是中高頻段數(shù)據(jù)嚴重畸變，不連續(xù)，無法進行下一步處理。用本文方法去噪后多處干擾頻點恢復正常，特別是中高頻數(shù)據(jù)變得穩(wěn)定連續(xù)，由于去噪作用的影響，低頻段數(shù)據(jù)有部分能量損失，電阻率值有小幅度波動(圖9a)，但是總體來看電阻率曲線和阻抗相位曲線變得更加穩(wěn)定連續(xù)，阻抗質量得到了較好改善(圖9a、b)。

圖9 實測信號去噪前后電阻率、阻抗對比

為了進一步評價實測信號去噪效果在此引入奈奎斯特圖[22-23]，在無噪聲干擾條件下，大地電磁阻抗的奈奎斯特圖是從低頻到高頻呈順時針展布的連續(xù)光滑曲線，一旦某個頻點受到干擾該頻點將會脫離這種趨勢，所以可以根據(jù)奈奎斯特圖是否更具有這種趨勢來評判去噪效果。去噪前，實測信號阻抗的奈奎斯特圖特征比較混亂，特別是在22.5 Hz、79 Hz、132 Hz、229 Hz等頻點處嚴重脫節(jié)，幾乎無法識別出其順時針旋轉特征(圖10a)。經(jīng)本文方法去噪后，除22.5 Hz較其他頻點偏離較大之外，其余頻點雖然也不是光滑連續(xù)分布，但是基本濾除了大尺度方波對阻抗的干擾，奈奎斯特圖總體趨勢基本符合從低頻到高頻順時針展布的特征(圖10b)，與去噪前相比阻抗結果更加穩(wěn)定和合理。

圖10 實測信號去噪前后的奈奎斯特圖對比

4 結論與建議

與前人將噪聲作為LSTM網(wǎng)絡的期望輸出不同的是，本文采用大地電磁信號本身作為網(wǎng)絡期望輸出，用訓練好的網(wǎng)絡進行了仿真含方波噪聲和實測含方波噪聲的大地電磁時間序列去噪測試，結果表明本文所提方法能有效消除方波噪聲干擾，改善阻抗估計結果，為深度學習在大地電磁時間序列去噪領域的應用提出了新思路。

本文所提方法仍具有局限性：一是本文只討論了方波噪聲的處理，對于其他典型的強干擾噪聲也可以借鑒本文的思路進行處理，但是怎樣讓網(wǎng)絡適應不同類型的噪聲需要研究；二是由于天然大地電磁信號的非平穩(wěn)性和隨機性[24-25]，在選取訓練數(shù)據(jù)時不能很好地囊括所有信號特征，只有當待處理數(shù)據(jù)的信號特征與訓練網(wǎng)絡所用大地電磁信號特征相似時才可以有效去除噪聲，當遇到網(wǎng)絡未學習過的信號時去噪效果就會大打折扣，所以怎樣更合理地構建訓練神經(jīng)網(wǎng)絡所需的數(shù)據(jù)集將是下一步需要解決的問題。

致謝：感謝中國地質調查局西安礦產資源調查中心郝子瓊對本文網(wǎng)絡訓練的支持！