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

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

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

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

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

1 基本理論

1.1 信號(hào)與噪聲模型

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

(1)

其中核函數(shù)為:

(2)

(3)

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

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

(4)

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

1.2 非負(fù)矩陣分解

非負(fù)矩陣分解(NMF)理論實(shí)質(zhì)上是利用非負(fù)約束條件來(lái)獲取數(shù)據(jù)表示的一種方法。NMF的理論問(wèn)題可以描述為:對(duì)于任意給定的一個(gè)非負(fù)矩陣Y,NMF算法能夠找到一個(gè)非負(fù)矩陣W和一個(gè)非負(fù)矩陣H,滿(mǎn)足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)化問(wèn)題:

(6)

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

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

(7)

(8)

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

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

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

(9)

(10)

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

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

(11)

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

降噪模型為:

(12)

(13)

(14)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖6 實(shí)測(cè)信號(hào)的時(shí)域波形

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

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

5 結(jié)論

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

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

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

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