本發(fā)明涉及大地電磁信號(hào)去噪技術(shù)領(lǐng)域，具體涉及一種基于組合濾波的大地電磁信號(hào)去躁方法。

背景技術(shù)：

由于大地電磁測(cè)深（Magnetotelluric，MT）使用天然場(chǎng)源，所以大地電磁測(cè)深信號(hào)極易受到噪聲的污染，如場(chǎng)源噪聲、地質(zhì)噪聲、人文噪聲。此外作業(yè)現(xiàn)場(chǎng)，風(fēng)會(huì)導(dǎo)致電磁探頭和信號(hào)傳輸線的擺動(dòng)，機(jī)器的振動(dòng)以及樹(shù)木的晃動(dòng)等導(dǎo)致地表微震，這些都會(huì)引起干擾噪聲，使得所測(cè)得的信號(hào)基線發(fā)生飄移。噪聲嚴(yán)重地影響了阻抗響應(yīng)參數(shù)的穩(wěn)定估計(jì), 尤其是近年來(lái)，隨著國(guó)民經(jīng)濟(jì)和現(xiàn)代化建設(shè)的發(fā)展，電磁背景噪聲(如無(wú)線電、工廠、鐵路、電網(wǎng)、礦場(chǎng)等)日益嚴(yán)重，所以如何有效地消除噪聲的干擾，提高采集數(shù)據(jù)質(zhì)量，是MT數(shù)據(jù)處理的首要任務(wù)?，F(xiàn)有技術(shù)中，對(duì)于影響大地電磁信號(hào)的噪聲主要是以壓制的方式進(jìn)行處理，噪聲去除不徹底，仍然影響大地電磁測(cè)探的精準(zhǔn)度。

經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）是一種完全自適應(yīng)的信號(hào)分解方法，它將信號(hào)分解為含有不同時(shí)間尺度且滿足以下兩個(gè)定義條件的一組（intrinsic mode function，IMF）： (1) 對(duì)于一列數(shù)據(jù)極值點(diǎn)和過(guò)零點(diǎn)數(shù)目必須相等或至多相差一點(diǎn)；(2) 在任意點(diǎn)，由局部極大點(diǎn)和極小點(diǎn)構(gòu)成的兩條包絡(luò)線的平均值為零。每個(gè)IMF可以被認(rèn)為是信號(hào)中固有的一個(gè)模態(tài)函數(shù)。信號(hào)經(jīng)過(guò)EMD分解后的各階模態(tài)函數(shù)能夠完全重構(gòu)，幾乎沒(méi)有能量損失。

而數(shù)學(xué)形態(tài)學(xué)基本思想是用具有一定形態(tài)的結(jié)構(gòu)元素去度量和提取信號(hào)中的對(duì)應(yīng)形狀，結(jié)構(gòu)元素的形態(tài)決定了所提出的形態(tài)信息。其基本運(yùn)算有膨脹、腐蝕、開(kāi)啟和閉合4種，其中腐蝕運(yùn)算的定義 ；膨脹運(yùn)算定義為；關(guān)閉和開(kāi)啟操作可分別列式如：，；廣義的數(shù)學(xué)形態(tài)閉-開(kāi)和開(kāi)-閉操作可表示為 ，； 最后廣義的數(shù)學(xué)形態(tài)濾波可寫成： 。

將經(jīng)驗(yàn)?zāi)B(tài)分解與數(shù)學(xué)形態(tài)學(xué)的組合，充分利用EMD 的多尺度分解及其可重構(gòu)特性和數(shù)學(xué)形態(tài)學(xué)濾波方法的優(yōu)點(diǎn)，為大地電磁信號(hào)的噪聲識(shí)別、去噪，改善數(shù)據(jù)質(zhì)量，提供了新途徑。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)上述現(xiàn)有技術(shù)存在的不足，本發(fā)明的目的是提供一種能有效去噪、改善數(shù)據(jù)質(zhì)量的基于組合濾波的大地電磁信號(hào)去躁方法。

為實(shí)現(xiàn)上述目的，本發(fā)明采用的技術(shù)方案如下：

一種基于組合濾波的大地電磁信號(hào)去躁方法，包括如下步驟：

（1）用EMD方法對(duì)原始MT信號(hào)x(t)進(jìn)行分解，得到高階模態(tài)函數(shù)h(t)和低階模態(tài)函數(shù)l(t)；

（2）對(duì)于低階模態(tài)函數(shù)l(t)，首先采用數(shù)學(xué)形態(tài)濾波方法對(duì)其進(jìn)行處理得到消噪后的部分g(t)，然后用l(t)減去g(t)，得到峰值信號(hào)f(t)；再用自適應(yīng)閾值去噪方法對(duì)峰值信號(hào)f(t)進(jìn)行二次分離，得到結(jié)果為f’(t)；最后，g(t) 和 f’(t)相加，其和就為低階模態(tài)函數(shù)l’(t)的去噪結(jié)果；

（3）對(duì)于高階模態(tài)函數(shù)h(t)，采用平滑濾波法來(lái)消除基線漂移，得到濾波結(jié)果h’(t)；

（4）將步驟（2）得到的l’(t)與步驟（3）得到的h’(t)相加，得到的和即為原始MT信號(hào)的去噪結(jié)果。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果：1、EMD分解是自適應(yīng)信號(hào)的，分解出來(lái)的IMF序列是多通帶濾波的結(jié)果，具有多分辨分析的能力，且可以完全重構(gòu)；2、充分結(jié)合經(jīng)驗(yàn)?zāi)B(tài)的多尺度分解及其可重構(gòu)特性和數(shù)學(xué)形態(tài)學(xué)濾波的有效提取有用信號(hào)能力等優(yōu)點(diǎn)，為大地電磁信號(hào)盡可能多的保留了有用信息，改善受了干擾數(shù)據(jù)的質(zhì)量，為進(jìn)一步估算功率譜估計(jì)和阻抗奠定了良好的基礎(chǔ)。

附圖說(shuō)明

圖1為本發(fā)明的工藝流程圖，其中Ⅰ Hilbert-Huang 變換；Ⅱ 數(shù)學(xué)形態(tài)濾波； Ⅲ 基于形態(tài)學(xué)的信號(hào)提出；Ⅳ 自適應(yīng)閾值濾波； Ⅴ 低階IMF 的消噪。

圖2為原始仿真信號(hào)和加噪后的信號(hào)，其中a為原始仿真信號(hào)，b為加噪后的信號(hào)。

圖3加噪仿真MT信號(hào)的EMD分解和對(duì)應(yīng)各階IMF分量的功率譜，其中a為EMD分解圖，b為各階IMF分量的功率譜。

圖4為去噪前后仿真信號(hào)的HHT時(shí)頻能量譜，其中a為去噪前的，b為去噪后的。

圖5為受三角波干擾的MT信號(hào)。

圖6為受三角波干擾Hy分量的EMD分解。

圖7為受三角波干擾Hx（b）的EMD分解。

圖8為三角波干擾抑制后的MT信號(hào)。

圖9為去噪前后大地電磁視電阻率曲線、相位曲線對(duì)比圖； a為 r_xy 曲線; b為r_yx曲線; c為 Φ_xy曲線; d為Φ_yx曲線。

具體實(shí)施方式

現(xiàn)結(jié)合具體實(shí)施例，來(lái)對(duì)本發(fā)明作進(jìn)一步的闡述。以下實(shí)施例僅為本發(fā)明的優(yōu)選方式，并非只用于限制本發(fā)明的保護(hù)范圍，任意與本發(fā)明相同或相似的技術(shù)方案均視為本發(fā)明的保護(hù)范圍。

實(shí)施例一

基于組合濾波的大地電磁信號(hào)去躁方法，包括如下步驟：

（1）根據(jù)EMD分解法，設(shè)原始的大地電磁信號(hào)為x(t)，首先確定x(t)上的所有極值點(diǎn)，然后將所有極值點(diǎn)用曲線擬合，得到信號(hào)的上、下包絡(luò)曲線， 分別設(shè)為和，則上、下包絡(luò)的平均曲線為：； 用減去后，得到剩余部分設(shè)為，即：。不是一個(gè)平穩(wěn)數(shù)據(jù)序列，為此重復(fù)以上過(guò)程k次，使得滿足IMF兩個(gè)定義條件：當(dāng)hk(t)滿足以下兩個(gè)條件：(1)整個(gè)信號(hào)中，零點(diǎn)數(shù)與極點(diǎn)數(shù)相等或至多相差1；(2)信號(hào)上任意一點(diǎn)，由局部極大值點(diǎn)確定的包絡(luò)線和由局部極小值點(diǎn)確定的包絡(luò)線的均值均為零，即信號(hào)關(guān)于時(shí)間軸局部對(duì)稱，則hk(t)就可定義為原始信號(hào)的第一個(gè)IMF。

即

。

此時(shí)，就是第一個(gè)IMF，C1（t）,[微軟用戶1] 是信號(hào)數(shù)據(jù)中最高頻率成分，用x(t)減去C1（t），得到一個(gè)去掉最高頻成分的新數(shù)據(jù)序列，重復(fù)步驟(1)次[微軟用戶2] ，

。

[微軟用戶3] 得到一系列c_n和最后一個(gè)序列r_n，r_n代表x(t)的均值或趨勢(shì)項(xiàng)。那么原序列x(t)可表示為所有IMF和一個(gè)殘余項(xiàng)的和：

；

將分解得到的模態(tài)函數(shù)中前半低序列號(hào)的IMF定義為低級(jí)模態(tài)函數(shù)l(t)，后半高序列號(hào)的IMF定義為高級(jí)模態(tài)函數(shù)h(t)。那么重寫信號(hào)x(t)的分解形式為：

，

從而將大地電磁信號(hào)的高頻部分和低頻分開(kāi)，然后分別進(jìn)行去噪處理。

（2）對(duì)于低階模態(tài)函數(shù)l(t)，先用數(shù)學(xué)形態(tài)學(xué)濾波方法對(duì)其進(jìn)行去噪處理，得到濾波結(jié)果g(t)，這個(gè)結(jié)果保留了MT 信號(hào)低階模態(tài)函數(shù)的分類信息：

，

上式中，k(t) (t=1,2, ... , M,M<<T )是結(jié)構(gòu)元素， 和 分別代表數(shù)學(xué)形態(tài)學(xué)里面的開(kāi)和閉操作；然后用l(t) 減去 g(t)得到峰值信號(hào)f(t)，

，

然而，通常情況下被分離出來(lái)的f(t)信號(hào)依然保留了一些有用的信息，特別是一些低頻信號(hào)保留了一些深層的地質(zhì)結(jié)構(gòu)信息，丟失這部分信息會(huì)大大減小MT的深層勘測(cè)能力。所以對(duì)f(t)進(jìn)行二次信噪分離是非常必要的，以至于能盡可能多的保留真實(shí)MT信號(hào)的信息。對(duì)f(t)用基于3σ準(zhǔn)則的自適應(yīng)閾值濾波方法進(jìn)行二次分離，獲得去噪后的信號(hào)f’(t)，這就是濾波后的峰值信號(hào)：

，

其中，f’(t)和來(lái)自數(shù)學(xué)形態(tài)濾波后g(t)相加得到去噪后的低階模態(tài)函數(shù)l’(t)。

（3）對(duì)于高階模態(tài)函數(shù)h(t)，采用平滑濾波法來(lái)消除基線漂移，對(duì)于MT信號(hào)來(lái)說(shuō)，基線漂移干擾主要存在于高階模態(tài)函數(shù)。對(duì)高階模態(tài)函數(shù)執(zhí)行多點(diǎn)平滑濾波，基線漂移有望得到抑制，計(jì)算格式如下：

，

此處，；

。

（4）現(xiàn)在，就可以構(gòu)造整個(gè)去噪后的大地電磁信號(hào)。來(lái)自低階模態(tài)函數(shù)部分的去噪結(jié)果l’(t) 加上來(lái)自高階模態(tài)函數(shù)部分的去噪結(jié)果h’(t)，其和就被視為受噪MT信號(hào)x(t)的去噪結(jié)果x’(t)，表達(dá)式如下：

。

EMD分解是自適應(yīng)信號(hào)的，分解出來(lái)的IMF序列是多通帶濾波的結(jié)果，具有多分辨分析的能力，且可以完全重構(gòu)，EMD分解為大地電磁信號(hào)的去噪提供了一條良好的途徑。而數(shù)學(xué)形態(tài)學(xué)是由一組形態(tài)學(xué)算子組成，形態(tài)學(xué)的基本運(yùn)算包括腐蝕、膨脹、開(kāi)閉。噪聲通常是在一定范圍內(nèi)作為一個(gè)峰值“峰頂”或“谷低”疊加在信號(hào)里面，形態(tài)學(xué)的開(kāi)操作可以用來(lái)剝離 “峰頂”的噪聲，形態(tài)學(xué)的閉操作可以用來(lái)填充的“谷低”的噪聲，二者都有濾波功能，在實(shí)際應(yīng)用中，將開(kāi)啟和關(guān)閉操作經(jīng)結(jié)合形成形態(tài)學(xué)濾波算法，提取信號(hào)中的有效信息。

為了驗(yàn)證上述組合濾波方法的正確性及其在處理MT信號(hào)方面的優(yōu)勢(shì)，發(fā)明人對(duì)一組仿真數(shù)據(jù)進(jìn)行了處理，仿真信號(hào)被加入了幾種不同的噪聲。在此，仿真信號(hào)被表示為，噪聲是，則這受噪的仿真MT信號(hào)可表示為：

；

圖2a所示為無(wú)噪的原始仿真MT信號(hào)，圖2b顯示了加噪后的信號(hào)(信噪比為：-10dB)，噪聲包含三角波干擾、工頻干擾、高頻噪聲和基線漂移。圖3a是含噪信號(hào)的EMD分解圖，圖3b是對(duì)應(yīng)的被分解后各階IMF分量的頻譜圖。從這些圖中可以看出，信號(hào)被自適應(yīng)地從高頻到低頻一層層的被分解，每一層都有不同的幅度和頻率值。從圖3還可以看出，高頻噪聲和工頻干擾主要存在于低階模態(tài)函數(shù)（主要在高階模態(tài)“imf1”-“imf4”）；基線漂移和工頻干擾的高次諧波存在于高階模態(tài)函數(shù)；三角波由于其豐富的頻率成分，分布在多階模態(tài)函數(shù)上。這些就為消噪處理帶來(lái)了方便。

按照本發(fā)明的方法，對(duì)于低階模態(tài)函數(shù)，先采用數(shù)學(xué)形態(tài)學(xué)濾波方法提取有用信號(hào)，然后再用自適應(yīng)濾波方法進(jìn)行二次分離；對(duì)于高階模態(tài)部分，采用多點(diǎn)平滑進(jìn)行去噪處理，最后兩部分的處理結(jié)果相加得到去噪后的信號(hào)。圖4展示了去噪前后兩種信號(hào)的時(shí)頻分布。從圖4a中也可看出，加噪后，信號(hào)的能量幅值增加了好幾個(gè)數(shù)量級(jí)，且在50Hz和低頻帶的頻率上出現(xiàn)了突出的奇異值。這些都是來(lái)自仿真加噪的三角波干擾、工頻干擾和基線漂移的影響，而在去噪后的的時(shí)頻分布圖4 b中，這些噪聲都得到了有效抑制。

為了評(píng)價(jià)這種組合濾波方法的良好效果，將處理結(jié)果與小波去噪方法的結(jié)果進(jìn)行了對(duì)比。并采用了一些評(píng)價(jià)參數(shù)如：SNR (信噪比), NSR (噪聲抑制率)和SDR (信號(hào)失真比)，來(lái)對(duì)去噪效果進(jìn)行評(píng)價(jià)，他們的定義如下：

信噪比: ；

噪聲抑制率:；

信號(hào)失真比:；

在此， 是信號(hào)的功率譜， 是噪聲的功率譜，是估算的去噪后的MT信號(hào)，SNR反映了有用信號(hào)被噪聲的污染程度，其值越小，說(shuō)明噪聲干擾越大。NSR 越大，SDR越小，則說(shuō)明去噪效果越好。

表1給出了幾個(gè)參數(shù)的對(duì)比結(jié)果：

從表中可以看出，在不同的SNR條件下，不管是本研究提出的組合方法還是小波軟閾值去噪方法，信號(hào)失真比SDR都不大，且前者有更低的SDR。就噪聲抑制率來(lái)說(shuō)，文中提出的組合方法也優(yōu)于小波軟閾值去噪方法。在組合方法的處理結(jié)果中，在較高的SNR環(huán)境下（10dB)，噪聲很弱，所以NSR也相對(duì)較小其值為0.1703，對(duì)應(yīng)的SDR也很小為0.0201。隨著信噪比的減小，當(dāng)SNR為10dB時(shí)，NSR也相對(duì)增大(從 0.1703 增加到0.7927)，盡管SDR也增大了（從 0.0201增大到0.0409），但是仍舊保持很小的值。這的就說(shuō)明文中提出的組合濾波方法有很好的去噪和信號(hào)細(xì)節(jié)保持能力。

通過(guò)以下實(shí)例來(lái)說(shuō)明本發(fā)明的處理效果。

實(shí)例一：

三角波干擾或由多個(gè)三角波連續(xù)出現(xiàn)構(gòu)成的鋸齒波狀干擾也是電磁測(cè)深作業(yè)中經(jīng)常遇到的干擾。特別是在作業(yè)區(qū)接近礦區(qū)或就在礦區(qū)內(nèi)，礦區(qū)內(nèi)井下存在大功率牽引機(jī)車等設(shè)備，造成大范圍的電磁干擾，成為研究區(qū)重要的干擾源，開(kāi)采時(shí)造成的震動(dòng)干擾和低頻干擾也時(shí)常發(fā)生，三角波干擾尤為常見(jiàn)。一般只影響磁場(chǎng)信號(hào)，在時(shí)間序列中，噪聲曲線突跳明顯，其形態(tài)表現(xiàn)為不規(guī)則的三角波形。通常具有如下表達(dá)式：

；

實(shí)測(cè)點(diǎn)西北方存在多處正在開(kāi)采的大型井下礦山，井下礦石運(yùn)輸采用的是大功率直流電力牽引機(jī)車。據(jù)實(shí)地調(diào)查，單臺(tái)機(jī)車工作電流為84 A，且回路為直接嵌入基巖的鐵軌，可形成大規(guī)模地下游散電流; 此外，大型機(jī)械( 如碎石機(jī)、風(fēng)鉆) 突然啟動(dòng)、關(guān)閉或負(fù)荷改變，均產(chǎn)生地電干擾。研究區(qū)內(nèi)主要地層為高阻的花崗侵入巖層。

如圖5所示是實(shí)測(cè)的受三角波干擾的MT數(shù)據(jù)。通常在磁道觀測(cè)信號(hào)中出現(xiàn)，且為磁場(chǎng)相關(guān)噪聲，電道極少受其影響。在圖5中，磁道Hx，Hy在記錄時(shí)間段內(nèi)出現(xiàn)明顯的三角狀噪聲，且相關(guān)性較好，噪聲局部能量較強(qiáng)，將正常磁場(chǎng)信號(hào)完全湮滅，而電道信號(hào)（Ex，Ey）卻無(wú)此類波形出現(xiàn)。單個(gè)噪聲幅值為正常信號(hào)幅值的數(shù)倍到數(shù)十倍。Hy信號(hào)的統(tǒng)計(jì)參數(shù)：最大值為6912 nT，最小值為-8588 nT，均值為39.65 nT，方差為5.31×10⁶，能量為4.35×10¹⁰nT²； Hx信號(hào)的統(tǒng)計(jì)參數(shù)：最大值為7951 nT，最小值為-7808 nT，均值為44.50 nT，方差為3.36×10⁶，能量為2.75×10¹⁰ nT²。

將受干擾的信號(hào)進(jìn)行EMD分解如圖6、圖7所示。用本文提出的組合方法對(duì)各階IMF進(jìn)行消噪處理，最后重構(gòu)得到消噪后的信號(hào)。圖8所示為受三角波干擾的大地電磁信號(hào)經(jīng)本文方法消噪后的效果圖。消噪后Hy信號(hào)的統(tǒng)計(jì)參數(shù)：最大值為6205 nT，最小值為-6980 nT，均值為-9.06 nT，方差為4.54×10⁵，能量為3.72×10⁹ nT²， Hy信號(hào)的統(tǒng)計(jì)參數(shù)：最大值為6182 nT，最小值為-7080 nT，均值為0.34 nT，方差為4.01×10⁵，能量為6.03×10⁹ nT²。消噪后鋸齒形狀沒(méi)有了，消噪后信號(hào)的方差減小到近原來(lái)的10％，信號(hào)變得平穩(wěn)。

用去噪前后的數(shù)據(jù)分別計(jì)算電阻率曲線和相位曲線。圖9 顯示了去噪前后視電阻率曲線和相位曲線的對(duì)比。我們用方差，均值來(lái)評(píng)估去噪前后兩參數(shù)曲線的平穩(wěn)性，用曲線相似性參數(shù)(NCC)從整體上評(píng)價(jià)去噪前后兩參數(shù)曲線的相似程度即整體趨勢(shì)。曲線相似性參數(shù)(NCC)定義如：

；

其中f(n)，g(n)分別為兩離散序列，NCC的值在-1~1之間，-1代表變換前后兩數(shù)據(jù)波形反向；0代表兩波形正交；1代表完全相同。

去噪前計(jì)算的ρ_xy曲線的方差是2.21×10⁴, 均值是313.64，去噪后計(jì)算的曲線方差是1.96×10⁴, 均值是312.78，兩曲線的相似性參數(shù)為0.905；去噪前計(jì)算的ρ_yx曲線的方差是1.13×10⁴, 均值是201.58，去噪后計(jì)算的曲線方差是1.11×10⁴, 均值是202.23，兩曲線的相似性參數(shù)為0.913。去噪前計(jì)算的TE 模式下相位曲線的方差是41.3, 均值是42.2，去噪后計(jì)算的曲線方差是32.3, 均值是42.2，兩曲線的相似性參數(shù)為0.917；去噪前計(jì)算的的TM 模式下相位曲線的方差是147.83, 均值是40.68，去噪后計(jì)算的曲線方差是123.34, 均值是40.70，兩曲線的相似性參數(shù)為0.925。

由圖5-9和數(shù)據(jù)的對(duì)比可知，由消噪前后計(jì)算的曲線相似性參數(shù)都較大，即整體趨勢(shì)是一直的，但細(xì)節(jié)發(fā)生了明顯的改善。所有去噪后的曲線方差都減小了很多，約為去噪前的一半，曲線的突變點(diǎn)得到了有效抑制（如ρ_yx曲線的3278 Hz頻點(diǎn)電阻率由原來(lái)的983歐姆米下降到723歐姆米）。所有響應(yīng)曲線的誤差棒平均值都減小了，曲線也更加圓滑、連續(xù)，參數(shù)的估算變得穩(wěn)定。