基于LSTM循環(huán)神經(jīng)網(wǎng)絡(luò)的大地電磁方波噪聲抑制

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

(1.中國地質(zhì)調(diào)查局 西安礦產(chǎn)資源調(diào)查中心，陜西 西安 710000； 2.中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

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

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

1 LSTM網(wǎng)絡(luò)簡介

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

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

圖1 LSTM網(wǎng)絡(luò)神經(jīng)元基本結(jié)構(gòu)

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時刻記憶細(xì)胞的候選值；ct為t時刻記憶細(xì)胞的狀態(tài)，也作為t+1時刻即下一個神經(jīng)元的初始記憶細(xì)胞；ht為t時刻的隱藏層，也為t+1時刻即下一個神經(jīng)元的初始隱藏層；Wi、Wf、Wo、Wg為傳播權(quán)重矩陣；bi、bf、bo、bg為偏置向量；⊙表示矩陣元素對應(yīng)相乘。

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

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

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

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

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

2 LSTM網(wǎng)絡(luò)搭建

2.1 網(wǎng)絡(luò)設(shè)計(jì)

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

圖2 LSTM網(wǎng)絡(luò)結(jié)構(gòu)示意

2.2 網(wǎng)絡(luò)訓(xùn)練

網(wǎng)絡(luò)訓(xùn)練主要包括數(shù)據(jù)集構(gòu)建、數(shù)據(jù)集歸一化(或標(biāo)準(zhǔn)化)、網(wǎng)絡(luò)訓(xùn)練和驗(yàn)證等幾個步驟。

2.2.1 數(shù)據(jù)集構(gòu)建

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

圖3 訓(xùn)練集原始時間序列

圖4 驗(yàn)證集原始時間序列

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

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

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

xnorm=(x-xmean)/xmax,

(7)

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

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

ynorm=y/xmax,

(8)

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

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

ypred=xmaxynorm_pred,

(9)

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

2.2.3 網(wǎng)絡(luò)的訓(xùn)練和驗(yàn)證

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

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

(10)

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

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

圖5 LSTM網(wǎng)絡(luò)學(xué)習(xí)曲線

3 去噪測試

3.1 仿真含噪信號去噪

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

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

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

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

3.2 實(shí)測含噪信號去噪

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

圖8 實(shí)測信號去噪前后對比

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

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

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

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

4 結(jié)論與建議

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

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

致謝：感謝中國地質(zhì)調(diào)查局西安礦產(chǎn)資源調(diào)查中心郝子瓊對本文網(wǎng)絡(luò)訓(xùn)練的支持！