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

(1.桂林電子科技大學(xué)信息與通信學(xué)院,廣西桂林541004;2.桂林電子科技大學(xué)認(rèn)知無線電與信息處理教育部重點(diǎn)實(shí)驗(yàn)室,廣西桂林541004)

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

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

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

1 基本理論

1.1 信號與噪聲模型

瞬變電磁法中,采用中心回線法采集時(shí),接收線圈在時(shí)刻t測得的感應(yīng)電動勢表達(dá)式為[16]:

(1)

其中核函數(shù)為:

(2)

(3)

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

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

(4)

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

1.2 非負(fù)矩陣分解

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

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

(5)

式中:‖·‖F(xiàn)為Frobenius范數(shù)。

最小化上述目標(biāo)函數(shù)可定義為如下優(yōu)化問題:

(6)

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

解上述優(yōu)化問題,使用如下乘性迭代規(guī)則,可使該算法收斂,W和H達(dá)到某個(gè)局部最優(yōu)點(diǎn)。

(7)

(8)

2 NMF的有監(jiān)督算法

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

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

(9)

(10)

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

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

(11)

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

降噪模型為:

(12)

(13)

(14)

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

圖1 平均次數(shù)與輸出信噪比的關(guān)系曲線

NMF的有監(jiān)督算法的實(shí)現(xiàn)步驟(圖2)總結(jié)如下:

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

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

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

圖2 NMF的有監(jiān)督算法流程

3 實(shí)驗(yàn)仿真及對比

3.1 實(shí)驗(yàn)說明

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

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

3.2 對比實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

總之,從輸出信噪比和均方根誤差這兩個(gè)降噪性能評價(jià)指標(biāo)來看,本文方法都要優(yōu)于其它3種方法。

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

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

4 實(shí)測數(shù)據(jù)應(yīng)用

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

圖6 實(shí)測信號的時(shí)域波形

圖7 實(shí)測數(shù)據(jù)降噪前、后的對比結(jié)果a 噪聲分析段; b 降噪后的瞬變電磁信號

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

5 結(jié)論

本文對NMF的有監(jiān)督算法如何應(yīng)用于瞬變電磁信號降噪進(jìn)行了研究,得出如下結(jié)論。

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

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

需要指出的是,本文方法的收斂速度較慢,主要是因?yàn)榉秦?fù)矩陣分解的速度較慢,需要幾千甚至上萬次的迭代計(jì)算,這增加了算法的時(shí)間。另外,本文方法對瞬變電磁信號的早中期數(shù)據(jù)處理效果不佳,這一點(diǎn)有待進(jìn)一步的研究。