NMF的有監督算法在瞬變電磁信號降噪中的應用

熊松齡,曾慶寧,龍 超,王師琦,祁瀟瀟,鄭展恒

(1.桂林電子科技大學信息與通信學院,廣西桂林541004;2.桂林電子科技大學認知無線電與信息處理教育部重點實驗室,廣西桂林541004)

瞬變電磁法是一種重要的地球物理勘探方法。該方法根據電磁感應原理,利用不接地回線或接地線源向地下發送一次場,在其激發下,地下地質體中激勵起二次場,通過分析、處理二次場信號,達到探測地下地質體的目的[1]。瞬變電磁法具有較高的探測和分辨能力,已被廣泛應用于煤礦、金屬礦和油氣田等地下目標的勘探工作[2-6]。

在野外作業時,用瞬變電磁法采集到的信號會受到各種噪聲的干擾,比如:天電噪聲、工頻干擾等,特別是在瞬變電磁信號的晚期,噪聲嚴重時,甚至會淹沒干凈的瞬變電磁信號,這些噪聲嚴重影響了瞬變電磁信號的后續處理工作。對此,在瞬變電磁信號處理前期,人們研究了各種噪聲的特性,提出了多種降噪方法,大致可分為3類:傳統方法、子空間法和基于深度學習的方法。傳統方法主要有小波變換[7-8]、Hilbert-Huang變換[9]和卡爾曼濾波[10]等。利用小波變換來對瞬變電磁信號降噪時,存在小波選擇困難的問題,也就是缺乏自適應性。Hilbert-Huang變換存在端點效應等問題,降噪效果不夠理想。對于卡爾曼濾波來說,這種適用于線性系統的降噪方法對于按指數規律下降的非線性瞬變電磁信號的處理能力有限。子空間法主要有奇異值分解[11]等。相比于傳統方法,雖然奇異值分解法降噪效果有所提高,但一直難以確定有效的奇異值,其處理步驟也比較復雜。近年來隨著深度學習的發展和應用,基于深度學習的方法在瞬變電磁信號處理領域受到較高的關注,主要有堆棧式降噪自編碼器[12]等。堆棧式降噪自編碼器的降噪效果優于傳統方法,但該模型必須建立在特定的地質條件上,適用范圍受到較大限制。一般來說,基于深度學習的方法都需要大量的數據作支撐,模型結構復雜,算法存在難以收斂的風險。

非負矩陣分解(NMF)的有監督算法是一種字典學習方法,在信源分離和語音增強等方面有著廣泛的應用。MOHAMMADIHA等[13]將此法用于估計噪聲的功率譜密度,并與維納濾波相結合進行語音增強。張龍[14]將較長的房間沖激響應分解為兩個較短響應的卷積,提出了一種兩步序貫處理的非負矩陣分解模型的語音增強算法。張星[15]采用對數譜估計方法對時頻譜圖中的語音存在概率進行估計,對非負矩陣分解得到的字典進行自適應補償,提高了語音增強的性能。本文首次將NMF的有監督算法引入瞬變電磁信號的降噪處理中,先研究了瞬變電磁信號模型和NMF的有監督算法,然后介紹了采用該算法對瞬變電磁信號進行降噪處理的技術流程,最后給出了該算法對模擬信號和實際瞬變電磁信號的實際應用。

1 基本理論

1.1 信號與噪聲模型

瞬變電磁法中,采用中心回線法采集時,接收線圈在時刻t測得的感應電動勢表達式為[16]:

(1)

其中核函數為:

(2)

(3)

實際采集的瞬變電磁信號混有噪聲。本文使用加性噪聲模型,即:

y(t)=v(t)+n(t)

(4)

式中:y(t)為含噪的瞬變電磁信號;v(t)為純凈的瞬變電磁信號;n(t)為加性噪聲。

1.2 非負矩陣分解

非負矩陣分解(NMF)理論實質上是利用非負約束條件來獲取數據表示的一種方法。NMF的理論問題可以描述為:對于任意給定的一個非負矩陣Y,NMF算法能夠找到一個非負矩陣W和一個非負矩陣H,滿足Y≈WH,從而將一個非負矩陣分解成兩個非負矩陣的乘積[17]。

給定一個非負矩陣Y,分解的算法有很多種[18-22]。本文使用的是基于歐式距離的NMF算法,其目標函數為:

(5)

式中:‖·‖F為Frobenius范數。

最小化上述目標函數可定義為如下優化問題:

(6)

式中:Wia和Hbj分別為非負矩陣W和H的元素。

解上述優化問題,使用如下乘性迭代規則,可使該算法收斂,W和H達到某個局部最優點。

(7)

(8)

2 NMF的有監督算法

將非負矩陣分解(NMF)的有監督算法應用于瞬變電磁信號的降噪處理,該算法可大致分為訓練、降噪和求平均3步。

在訓練階段,得到表示純凈的瞬變電磁信號矩陣Vt和噪聲信號矩陣Nt后,分別對其進行非負矩陣分解,可通過最小化如下目標函數完成:

(9)

(10)

利用公式(7)和公式(8)求解該優化問題后,分別得到表示純凈的瞬變電磁信號矩陣Vt和噪聲信號矩陣Nt特征的原子字典Wv和Wn。然后,將Wv和Wn拼接成總字典W=[WvWn],矩陣W包含了表示純凈的瞬變電磁信號矩陣Vt和噪聲信號矩陣Nt的特征。

在降噪階段,利用在訓練階段得到的總字典W對電磁信號進行處理,即在固定總字典W的情況下,將表示含噪的瞬變電磁信號的矩陣Yd進行非負矩陣分解,公式如下:

(11)

式中:arg表示取(11)式范數最小值的Hd。

降噪模型為:

(12)

(13)

(14)

實驗表明,由前兩步處理得到的降噪后的瞬變電磁信號不夠理想。本文通過對前兩步進行多次循環,再求平均來提高降噪效果。顯然,該求平均值的平均次數是一個關鍵參數。本文對此進行了實驗,實驗是在10dB的高斯白噪聲環境下進行,結果如圖1所示。從圖1可以看出,在平均次數小于50次時,輸出信噪比隨著平均次數的增加呈現出上升趨勢;在平均次數大于50次時,輸出信噪比不再隨平均次數的增加而增加。因此,在保證較低的時間、空間復雜度和較高質量的降噪效果的前提下,本文選擇的平均次數為50次。

圖1 平均次數與輸出信噪比的關系曲線

NMF的有監督算法的實現步驟(圖2)總結如下:

2)對含噪的瞬變電磁信號y、純凈的瞬變電磁信號v和噪聲信號n分別進行短時傅里葉變換,得到各自的幅度譜和相位譜,取其幅度譜,得到表示相應信號的非負矩陣Y、V和N;

3)將表示純凈的瞬變電磁信號的非負矩陣V和噪聲信號的非負矩陣N進行非負矩陣分解,分別得到表示各自信號特征的原子字典Wv和Wn,然后將其拼接成總字典W=[WvWn];

4)在固定總字典W的前提下,對表示含噪信號的非負矩陣Y進行非負矩陣分解,得到系數矩陣H;

圖2 NMF的有監督算法流程

3 實驗仿真及對比

3.1 實驗說明

通過實驗仿真及對比的方式來驗證非負矩陣分解(NMF)有監督算法的有效性。該方法適用于數據采集后的預處理階段,采用的信號模型為加性噪聲模型。在本實驗中,加入的噪聲為高斯白噪聲和天電噪聲。天電噪聲由Alpha穩定分布[23-24]模擬。純凈的瞬變電磁信號使用公式(1)模擬,發射天線St的有效面積取10m2,發射電流I取5A,接收天線的有效面積Sr取100m2,視電阻率ρ取大地平均值50Ω·m。另外,為了與本文采用的實測數據的個數一致(本文的實測數據由重慶地質儀器廠的ATEM-Ⅱ型瞬變電磁儀采集,密集均勻采樣,樣本數為1000個),模擬信號由采樣頻率為50kHz的均勻采樣獲得,采樣時間為20ms,采樣點數為1000個,取其前40個數據為早中期數據,后960個數據為晚期數據。為了保證實驗的可信度,噪聲訓練數據與加在純凈瞬變電磁信號上的噪聲為同一類型,只是數值不同。輸入信噪比分別為-5,0,5,10,15dB。

在對比實驗中,本文選取小波變換、Hilbert-Huang變換(HHT)和奇異值分解(SVD)3種方法與本文方法進行對比。在小波變換中,本文使用的是sym8小波,硬閾值函數,啟發式閾值,5層分解。在HHT中,本文通過固有模態函數(IMF)分量的能量密度和平均周期的乘積確定有用的IMF分量。在SVD中,本文采用趙學智等[25]提出的奇異值曲率譜方法確定有效的奇異值。

3.2 對比實驗

首先,對加入天電噪聲的瞬變電磁信號進行降噪處理。天電噪聲是瞬變電磁信號噪聲的主要來源之一,其在波形上表現為突出的毛刺。圖3a是純凈的瞬變電磁時域信號,圖3b是加入了15dB天電噪聲的瞬變電磁信號。圖4是分別使用HHT、SVD、小波變換和本文方法降噪后的信號。

圖3 純凈信號與含噪信號a 純凈的瞬變電磁時域信號; b 加入了15dB高斯白噪聲的瞬變電磁信號

由圖4可見,小波變換對天電噪聲的消除基本沒有作用,經HHT處理后仍有較多毛刺,但經SVD和本文方法降噪處理后,天電噪聲基本被消除,看不出有明顯的毛刺,降噪效果較好。

圖4 在15dB天電噪聲環境下采用不同方法處理后的瞬變電磁信號a HHT; b SVD; c小波變換; d 本文方法

對比在不同噪聲環境下分別采用HHT、SVD、小波變換和本文方法的輸出信噪比,結果如表1和表2所示。一般來說,輸出信噪比越大,混在信號中的噪聲就越小。由表2可見,在天電噪聲的環境下,經小波變換處理后,輸出信噪比與輸入信噪比基本相同,這表明了小波變換對天電噪聲基本沒有抑制能力。HHT對高斯白噪聲和天電噪聲消除效果大致相同,但輸出信噪比不高。在輸入信噪比為10dB和15dB的情況下,相比于其它噪聲環境,采用SVD的輸出信噪比提升幅度不大,此方法對兩種噪聲抑制效果相當,但其輸出信噪比要比小波變換和HHT的輸出信噪比高,與本文方法相比,其輸出信噪比相對較低。

表1 高斯白噪聲采用不同方法處理后的輸出信噪比

表2 天電噪聲下采用不同方法處理后的輸出信噪比

經本文方法降噪處理后,無論在何種噪聲環境下,其輸出信噪比都有大幅度的提升,而且要比采用SVD處理的高十幾個分貝。

表3和表4分別為在高斯白噪聲環境下和天電噪聲情況下采用4種方法處理后的均方根誤差。均方根誤差體現了原始信號和降噪之后的估計信號間的差異程度。降噪后的均方根誤差越小,表示原始有效信號和降噪之后的估計信號之間的差異越小,兩者越接近,表明降噪效果越好。由表3可以看出,在高斯白噪聲環境下,采用HHT處理后的均方根誤差最大,采用小波變換、SVD和本文方法基本上依次減小,說明HHT對高斯白噪聲的消除能力最差,本文方法的降噪能力最高。由表4可以看出,在天電噪聲的環境下,小波變換對天電噪聲的消除能力最差,經本文方法降噪處理后的均方根誤差依然最小,降噪效果最好。

表3 高斯白噪聲下采用不同方法處理后的均方根誤差

表4 天電噪聲下采用不同方法處理后的均方根誤差

總之,從輸出信噪比和均方根誤差這兩個降噪性能評價指標來看,本文方法都要優于其它3種方法。

需要說明的是,本文方法不僅對視電阻率為50Ω·m的含噪的瞬變電磁信號有效,對其它視電阻率的含噪信號同樣有降噪效果。另外,雖然本實驗是采用視電阻率為50Ω·m的純凈的瞬變電磁信號進行訓練,得到表征其特征的原子字典,但利用此原子字典對其它視電阻率的含噪的瞬變電磁信號也起到降噪作用,換句話說,本文方法對視電阻率的大小并不敏感。如圖5所示,對視電阻率為150Ω·m的含噪信號進行降噪,用視電阻率為50Ω·m的純凈的瞬變電磁信號進行訓練,經本文方法處理后,降噪后的信號依然與視電阻率為150Ω·m的純凈信號接近?？梢灶A見,這一重要特性將有利于本文方法的實際應用。

圖5 不同視電阻率與降噪信號的對比結果

4 實測數據應用

采用本文方法對廣東某地區的實測數據進行降噪處理。實驗設備為重慶地質儀器廠的ATEM-Ⅱ型瞬變電磁儀,采用中心回線裝置工作,密集均勻采樣,采樣時間為20ms,測得的數據未經過疊加處理。取10號線第38測點1000道數據中第381～900道作為噪聲分析段,其中第381～400道數據作為早中期數據,第401～900道數據作為晚期數據,并使用NMF的有監督算法進行降噪處理。純凈的瞬變電磁信號使用公式(1)仿真得到,相關參數與仿真實驗一致。由于采集數據的最后一部分瞬變電磁信號比較微弱,噪聲占主導地位,可將這部分數據看作噪聲,則噪聲訓練數據將由實測數據的最后100道數據組成。另外,本實驗對實測數據進行分段處理,其中晚期數據為5段,每100道數據為一段,早中期數據為單獨一小段(早中期數據不需要進行非負矩陣分解處理),以此來解決訓練數據與含噪數據維度不匹配的問題。圖6是在某地測得的瞬變電磁信號,其噪聲分析段信號與降噪后的瞬變電磁信號分別如圖7a和圖7b 所示。

圖6 實測信號的時域波形

圖7 實測數據降噪前、后的對比結果a 噪聲分析段; b 降噪后的瞬變電磁信號

由圖6和圖7a可見,實測數據存在明顯的表現為毛刺的天電噪聲和其它噪聲。由圖7b可以明顯看出,降噪后的波形無明顯的毛刺,比較平滑,噪聲壓制得比較干凈,降噪效果較好。這說明本文方法對實測數據是可行有效的。

5 結論

本文對NMF的有監督算法如何應用于瞬變電磁信號降噪進行了研究,得出如下結論。

1) 將含噪的瞬變電磁信號分為早中期和晚期兩部分,對含噪信號采用NMF的有監督算法處理,步驟重復多次,并將處理結果的晚期數據與含噪信號的原始早中期數據分別累加以求算術平均值,最后拼接出完整的降噪信號。結果表明,本文方法能有效地抑制高斯白噪聲和天電噪聲,處理后的波形比較平滑,噪聲壓制得比較干凈。另外,從輸出信噪比和均方根誤差這兩個降噪性能評價指標看,即使在低信噪比的情況下,采用本文方法也能較大幅度地提高輸出信噪比和降低均方根誤差,能取得較好的降噪效果,提高了數據的質量,有利于瞬變電磁信號的后續處理。

2) 相比于傳統的降噪方法,本文方法是一種基于字典學習的降噪算法,能學習到信號的特征。理論上本文方法能對混有多種噪聲的信號進行降噪處理(如非平穩噪聲),局限性較小,降噪效果也較好。

需要指出的是,本文方法的收斂速度較慢,主要是因為非負矩陣分解的速度較慢,需要幾千甚至上萬次的迭代計算,這增加了算法的時間。另外,本文方法對瞬變電磁信號的早中期數據處理效果不佳,這一點有待進一步的研究。