曹小玲,陳清禮,蔣 濤

(1.長(zhǎng)江大學(xué)油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,湖北武漢430100;2.長(zhǎng)江大學(xué)信息與數(shù)學(xué)學(xué)院,湖北荊州434023;3.長(zhǎng)江大學(xué)電工電子國(guó)家級(jí)實(shí)驗(yàn)教學(xué)示范中心,湖北荊州434023)

大地電磁測(cè)深(magnetotelluric sounding,MT)是研究地球深部構(gòu)造的重要方法之一,該方法利用天然電磁場(chǎng)作為場(chǎng)源因而易受到各種干擾的影響,其中一種重要的干擾類型為工頻干擾。工頻干擾是大地電磁信號(hào)采集過(guò)程中最為普遍的干擾之一,近年大地電磁信號(hào)中的工頻干擾噪聲愈加嚴(yán)重,極大影響了大地電磁勘探的效果。工頻干擾噪聲的信號(hào)能量強(qiáng),一般的工頻干擾噪聲能量為正常大地電磁信號(hào)的幾倍,有的工頻干擾噪聲能量甚至比正常大地電磁信號(hào)高幾個(gè)數(shù)量級(jí)。在高壓輸電線附近,工頻干擾噪聲信號(hào)能量尤其巨大,常常使得大地電磁有效信號(hào)完全淹沒在工頻干擾噪聲中,嚴(yán)重影響后續(xù)的數(shù)據(jù)處理和反演解釋[1-3]。隨著國(guó)內(nèi)外對(duì)大地電磁信號(hào)噪聲壓制的重視和持續(xù)研究,一些學(xué)者針對(duì)廣泛存在的工頻干擾噪聲問題提出了解決方案。蔡劍華等[4]提出先對(duì)包含工頻干擾噪聲的MT信號(hào)進(jìn)行傅里葉變換,得到其實(shí)部和虛部后,再用小波閾值去噪法分別對(duì)實(shí)部及虛部進(jìn)行去噪處理,最后將去噪后的實(shí)部和虛部聯(lián)合起來(lái)進(jìn)行傅里葉逆變換得到去噪后的信號(hào)。TANG等[5]提出先對(duì)采集的大地電磁信號(hào)進(jìn)行傅里葉變換,然后設(shè)計(jì)與干擾信號(hào)相匹配而對(duì)有用信號(hào)不敏感的冗余字典原子,再結(jié)合改進(jìn)的正交匹配追蹤(improved orthogonal matching pursuit,IOMP)算法分離出頻域信號(hào)中的工頻干擾成分,最后將處理后的頻域信號(hào)進(jìn)行傅里葉逆變換得到去噪后的信號(hào)。上述兩種方法壓制大地電磁信號(hào)工頻干擾噪聲效果較好,但均為基于傅里葉變換及其反(逆)變換的方法,故往往受限于傅里葉變換的優(yōu)勢(shì)和劣勢(shì)。CAO等[6]和曹小玲等[7]提出基于盲源分離的大地電磁工頻干擾噪聲壓制方法,將受到工頻干擾噪聲的大地電磁信號(hào)進(jìn)行離散小波變換(discrete wavelet transform,DWT)和集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)處理,而后再進(jìn)行盲源分離以壓制噪聲,該方法采用了基于DWT的處理模型并引入了自適應(yīng)權(quán)重因子,不僅有效降低了盲源分離算法恢復(fù)信號(hào)幅值的不確定性,而且在處理含有工頻干擾噪聲且信噪比較低的MT數(shù)據(jù)時(shí)仍然效果良好。該方法采用小波變換方法構(gòu)造多通道信號(hào),因小波變換基函數(shù)的選擇和分解層次的選擇過(guò)于依賴于人工經(jīng)驗(yàn)及試驗(yàn)分析,故方法的自適應(yīng)性不強(qiáng)。我們?cè)谘芯恐幸恢痹噲D避免采用小波變換方法,探索不過(guò)分依賴人工經(jīng)驗(yàn)及試驗(yàn)分析的信號(hào)分解方法和工頻干擾噪聲壓制方法。

經(jīng)驗(yàn)?zāi)B(tài)分解[8-9](empirical mode decomposition,EMD)顯著特點(diǎn)為:根據(jù)信號(hào)特性通過(guò)迭代的方式自適應(yīng)地獲取基函數(shù)和分解層次,即無(wú)需事先人為給定基函數(shù)和分解層次,無(wú)需依賴人工經(jīng)驗(yàn)及試驗(yàn)分析,因其具有良好的自適應(yīng)性故被廣泛應(yīng)用于噪聲壓制。盲源分離(blind source separation,BSS)是近年來(lái)發(fā)展起來(lái)的一種現(xiàn)代信號(hào)處理技術(shù),已在通信信號(hào)處理、機(jī)械故障識(shí)別、語(yǔ)音信號(hào)去噪、地震信號(hào)去噪等領(lǐng)域獲得廣泛應(yīng)用,前期我們將其應(yīng)用于大地電磁信號(hào)去噪,取得了較好的效果[10]。為了結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解和盲源分離的優(yōu)勢(shì)并彌補(bǔ)二者的不足,本文提出一種基于單通道盲源分離的大地電磁工頻干擾噪聲壓制方法,主要利用工頻干擾噪聲的特點(diǎn)、EMD及盲源分離的優(yōu)良特性,首先將含有工頻干擾噪聲的大地電磁信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和主成分分析(principal component analysis,PCA),然后利用衰減率確定有效的固有模態(tài)函數(shù)(intrinsic mode function,IMF)的個(gè)數(shù)(即源信號(hào)個(gè)數(shù),也即有效主成分的個(gè)數(shù)),最后利用盲源分離中的獨(dú)立分量分析(independent component analysis,ICA)方法對(duì)有效主成分進(jìn)行處理,以壓制工頻干擾噪聲。本文方法根據(jù)觀測(cè)信號(hào)確定源信號(hào)的個(gè)數(shù),解決了在欠定盲源分離情況下源信號(hào)個(gè)數(shù)無(wú)法確定的問題。

1 大地電磁工頻干擾噪聲特征

大地電磁天然場(chǎng)源的特點(diǎn)[11-12]決定了大地電磁信號(hào)非常容易受到噪聲的干擾,且這些噪聲形態(tài)各異、種類繁多、強(qiáng)弱不一。根據(jù)噪聲產(chǎn)生的機(jī)理劃分,大地電磁測(cè)深信號(hào)包含的噪聲包括工頻干擾噪聲、儀器噪聲、人文噪聲、地質(zhì)噪聲等。工頻干擾噪聲一般指電力傳輸系統(tǒng)、人工電磁場(chǎng)、通信過(guò)程及有線廣播等產(chǎn)生的噪聲。這些電磁噪聲呈現(xiàn)非平面波的形態(tài),且與觀測(cè)點(diǎn)距離較近,因此不符合大地電磁測(cè)深對(duì)場(chǎng)源的要求。

因工頻干擾噪聲很大一部分是人為導(dǎo)致的,故也有學(xué)者認(rèn)為應(yīng)將其歸為人文噪聲。如文獻(xiàn)[1]認(rèn)為大地電磁測(cè)深中的噪聲包括場(chǎng)源噪聲、地質(zhì)噪聲和人文噪聲3個(gè)類別,其中人文噪聲主要是指“人類活動(dòng)所產(chǎn)生的干擾性電磁波以及人類的活動(dòng)產(chǎn)生的其它電磁噪聲”[1],這里的“人文噪聲”包括了工頻干擾噪聲,即認(rèn)為工頻干擾噪聲是人文噪聲中的一部分。

工頻干擾噪聲主要特征為:①等頻率(50Hz)且等振幅,通常導(dǎo)致原始的大地電磁有效信號(hào)幾乎完全被淹沒而不可見;②形態(tài)較為規(guī)則,一般呈現(xiàn)出類似于正弦曲線的形狀;③幅值通常很大,因?yàn)榕c原始有效信號(hào)幅值對(duì)比強(qiáng)烈,故常常使得原始有效信號(hào)無(wú)法被識(shí)別,因此對(duì)電磁勘探等實(shí)際生產(chǎn)的破壞性非常大;④諧波干擾(一般為奇次諧波干擾)嚴(yán)重,影響后期的數(shù)據(jù)處理和反演工作。本文從實(shí)際大地電磁信號(hào)中挑選出了一段包含工頻干擾噪聲的實(shí)測(cè)大地電磁信號(hào)(圖1),該測(cè)點(diǎn)地址位于31:08.487(N),114:06.865(E)。一般而言,工頻干擾噪聲可能同時(shí)出現(xiàn)在所有電道和磁道中,如圖1所示Ex、Ey、Hx、Hy、Hz各道信號(hào)中都出現(xiàn)了工頻干擾噪聲,原始有效信號(hào)幾乎被完全覆蓋。

圖1 包含工頻干擾噪聲的實(shí)測(cè)大地電磁信號(hào)

2 方法原理

2.1 盲源分離基本原理

盲源分離[13]的流程如圖2所示,其中,s(t)=[s1(t),s2(t),…,sn(t)]T是n維未知源信號(hào)向量,A是未知混合系統(tǒng),x(t)=[x1(t),x2(t),…,xm(t)]T是m維觀測(cè)信號(hào)向量,它們是源信號(hào)向量受噪聲向量n(t)=[n1(t),n2(t),…,nm(t)]T干擾而形成的。x(t)被分離系統(tǒng)W處理后得到輸出信號(hào)y(t),它是源信號(hào)向量s(t)的估計(jì)。

圖2 盲源分離流程

盲源分離的目的是在源信號(hào)s(t)和混合系統(tǒng)A及噪聲n(t)均未知的情況下,利用觀測(cè)信號(hào)向量x(t)調(diào)整分離系統(tǒng)W(尋求分離矩陣W),使得輸出y(t)是源信號(hào)s(t)的估計(jì),即:

y(t)=W(x(t))?s(t)

(1)

根據(jù)王書明等[14-15]對(duì)不同地區(qū)實(shí)測(cè)大地電磁信號(hào)的性質(zhì)特征分析可知,大地電磁信號(hào)符合盲源分離對(duì)源信號(hào)的統(tǒng)計(jì)性質(zhì)要求。

獨(dú)立分量分析(ICA)方法是近年發(fā)展的一種有效盲源分離技術(shù),其應(yīng)用范圍極為廣泛[16-17]。從線性變換和線性空間角度來(lái)看,源信號(hào)為相互獨(dú)立的非高斯信號(hào),可以看作線性空間的基信號(hào),而觀測(cè)信號(hào)為源信號(hào)的線性組合,ICA方法是在源信號(hào)和線性變換均不可知的情況下,從觀測(cè)的混合信號(hào)中估計(jì)數(shù)據(jù)空間的基本結(jié)構(gòu)或者信源信號(hào)[7]。ICA方法的算法流程[7]如圖3所示,在信源信號(hào)s(k)中各分量相互獨(dú)立的假設(shè)下,由觀察信號(hào)x(k)通過(guò)解混系統(tǒng)B將信源信號(hào)分離,使輸出y(k)逼近s(k)。ICA方法中應(yīng)用較多的是FastICA算法,本文也采用此算法。

圖3 ICA方法的算法流程

2.2 EMD方法基本原理

EMD方法將相對(duì)復(fù)雜的原始信號(hào)分解成有限個(gè)相對(duì)簡(jiǎn)單且具有不同特征尺度的數(shù)據(jù)序列,即固有模態(tài)函數(shù)(IMF)分量,反映了原始信號(hào)的本質(zhì)和真實(shí)信息,因此我們對(duì)其進(jìn)行操作和處理就等同于對(duì)原始信號(hào)進(jìn)行操作和處理。EMD方法的基本步驟如下。

1) 假設(shè)原始信號(hào)為x(t),首先我們找出x(t)中所有的局部極值點(diǎn)(包括極大值點(diǎn)和極小值點(diǎn)),用三次樣條函數(shù)連接所有的局部極大值點(diǎn)構(gòu)成上包絡(luò)線,用三次樣條函數(shù)連接所有的局部極小值點(diǎn)構(gòu)成下包絡(luò)線。

2) 求出上包絡(luò)線與下包絡(luò)線的均值(設(shè)為m1),我們?cè)O(shè)信號(hào)x(t)與m1的差為h1,h1=x(t)-m1,如果h1滿足IMF的條件,那么h1就是原始信號(hào)的第一個(gè)IMF分量。

3) 如果h1不是IMF分量,我們就將h1看成原始信號(hào)進(jìn)行處理,再根據(jù)上述步驟,得到h11=h1-m11,其中,m11是h1信號(hào)的上、下包絡(luò)線均值。按照上面的操作反復(fù)處理k次得到h1k,如果h1k滿足IMF的條件,h1k被稱為IMF,且有h1k=h1(k-1)-m1k,c1=h1k為原始信號(hào)的第一個(gè)IMF分量,代表x(t)最高頻率的分量。

4) 從x(t)中分離c1,得到r1=x(t)-c1,將r1視為原始信號(hào)重復(fù)上述過(guò)程,可以得到信號(hào)x(t)的第二個(gè)IMF分量c2。重復(fù)上述過(guò)程n次,得到n個(gè)IMF分量,依次為c1,c2,…,cn,且r2=r1-c2,…,rn=rn-1-cn。當(dāng)rn為一個(gè)極小的常量或單調(diào)函數(shù)時(shí),就無(wú)法提取IMF,此時(shí)我們終止分解過(guò)程,可以得到:

(2)

其中,rn稱為殘差,它反映了信號(hào)x(t)的集中趨勢(shì)。IMF分量c1,c2,…,cn包含了信號(hào)的特征尺度信息,也包含了信號(hào)的特征頻率信息,因此利用c1,c2,…,cn就可以了解原始信號(hào)完整的特征尺度信息和特征頻率信息。

本文沒有選擇EEMD而選擇EMD處理的原因如下:一是因?yàn)镋EMD是在信號(hào)中添加白噪聲后再進(jìn)行EMD處理,添加白噪聲可能對(duì)原有的工頻干擾噪聲產(chǎn)生影響;二是因?yàn)榇蟮仉姶判盘?hào)普遍數(shù)據(jù)量龐大,采用EEMD處理會(huì)極大地增加計(jì)算量。

2.3 基于單通道盲源分離的大地電磁工頻干擾噪聲壓制方法

大地電磁信號(hào)通常是由多道信號(hào)數(shù)據(jù)構(gòu)成,由于每道信號(hào)數(shù)據(jù)包含的噪聲形式和噪聲數(shù)量不盡相同,因此我們需要對(duì)每道信號(hào)數(shù)據(jù)分別進(jìn)行去噪處理,即需要進(jìn)行單通道大地電磁工頻干擾噪聲壓制。盲源分離一般要求觀測(cè)信號(hào)的數(shù)目不少于源信號(hào)的數(shù)目,而單通道大地電磁工頻干擾噪聲壓制問題屬于欠定盲源分離問題,即觀測(cè)信號(hào)的數(shù)目少于源信號(hào)的數(shù)目,所以必須構(gòu)造其它觀測(cè)信號(hào),增加觀測(cè)信號(hào)的個(gè)數(shù),使得觀測(cè)信號(hào)的數(shù)目不少于源信號(hào)的數(shù)目,才能應(yīng)用盲源分離方法。本文提出的基于經(jīng)驗(yàn)?zāi)B(tài)分解和盲源分離的大地電磁工頻干擾噪聲壓制方法,可以實(shí)現(xiàn)單通道大地電磁工頻干擾噪聲壓制。該方法的流程如圖4所示,首先采用EMD將觀測(cè)信號(hào)分解為若干個(gè)IMF分量(這里假設(shè)為M個(gè)),即實(shí)現(xiàn)觀測(cè)信號(hào)的升維,以滿足盲源分離對(duì)源信號(hào)數(shù)目的要求;然后將IMF分量組成矩陣,應(yīng)用奇異值分解方法求取矩陣特征值,根據(jù)特征值比求衰減率,根據(jù)衰減率確定源信號(hào)的數(shù)目和有效主成分的個(gè)數(shù);再對(duì)獲得的IMF分量進(jìn)行PCA分析,確定并選擇起主要作用的主成分分量;為進(jìn)一步壓制噪聲和減少主成分信號(hào)之間的相關(guān)噪聲,最后對(duì)選取的若干個(gè)有效主成分分量進(jìn)行ICA處理,獲得壓制工頻干擾噪聲之后的信號(hào)。

圖4 基于經(jīng)驗(yàn)?zāi)B(tài)分解和盲源分離的大地電磁工頻干擾噪聲壓制方法流程

該方法主要優(yōu)勢(shì)在于EMD模型、PCA處理及盲源分離方法的應(yīng)用,受文獻(xiàn)[18]的啟發(fā),我們利用衰減率確定源信號(hào)的數(shù)目和待處理的主成分的個(gè)數(shù),解決了欠定盲源分離中源信號(hào)個(gè)數(shù)無(wú)法確定的問題,實(shí)現(xiàn)了單通道觀測(cè)信號(hào)的源信號(hào)的分離和提取。

PCA處理在盲源分離處理領(lǐng)域是一種極其有效的處理手段,我們對(duì)EMD分解得到的IMF分量進(jìn)行PCA處理[19],可以有效降低信號(hào)數(shù)據(jù)的維數(shù),通過(guò)尋找信號(hào)中能量最大的分量得到信號(hào)的主要特征量,即起到“化繁為簡(jiǎn)”的作用。另外,PCA處理作為一種統(tǒng)計(jì)分析方法,還能實(shí)現(xiàn)對(duì)信號(hào)的去相關(guān)處理。因?yàn)樗軌蚴沟米儞Q后產(chǎn)生的新分量正交或不相關(guān),故變換后的分量能量更集中、性質(zhì)更穩(wěn)定[20]。PCA處理的優(yōu)勢(shì)在于利用某種變換將數(shù)據(jù)原有的大量特征變換為幾個(gè)主要特征,而這些特征包含了原始數(shù)據(jù)的最主要的特征信息。因此,當(dāng)我們僅僅提取前面的部分特征時(shí),可以既減小特征的個(gè)數(shù)又保留原有信號(hào)最主要的特征信息。同時(shí),由于PCA處理后的信號(hào)不相關(guān),所以噪聲子空間和源信號(hào)子空間相互分離,進(jìn)而可以達(dá)到初步壓制部分噪聲的目的,即初步壓制部分工頻干擾噪聲。

PCA處理雖然能夠去除信號(hào)的相關(guān)性,但在高階累積量的定義下信號(hào)可能仍然具有相關(guān)性。為使所得信號(hào)具有的統(tǒng)計(jì)獨(dú)立性最大,我們又進(jìn)行了ICA處理,選取有用分離向量與混合矩陣相乘重構(gòu)得到分離后的獨(dú)立分量信號(hào),進(jìn)而得到消噪信號(hào)。

3 仿真模擬實(shí)驗(yàn)

如圖5a所示,利用蒙特卡洛方法對(duì)仿真信號(hào)中的原始信號(hào)隨機(jī)構(gòu)造3組信號(hào)(幅值為400)。因?yàn)楣ゎl干擾噪聲一般為50Hz的周期信號(hào)及其諧波,故構(gòu)造工頻干擾噪聲信號(hào)(噪聲信號(hào))可表示為S=800sin(2π×50t)?，F(xiàn)實(shí)情況下強(qiáng)烈的工頻干擾噪聲主要來(lái)自電力傳輸系統(tǒng),其頻率一般為50Hz,從方便分析的角度考慮,本文設(shè)工頻干擾噪聲為50Hz。因工頻干擾噪聲信號(hào)幅值一般是有效電磁信號(hào)幅值的數(shù)倍,所以這里幅值應(yīng)設(shè)較大值,即為原始信號(hào)的數(shù)倍(圖5b)。圖5c為工頻干擾噪聲加入原始信號(hào)后得到的含噪信號(hào)。

圖5 3組含噪信號(hào)構(gòu)建a 隨機(jī)信號(hào); b 50Hz正弦信號(hào)的工頻干擾噪聲; c 含噪信號(hào)

對(duì)圖5c所示的含噪信號(hào)進(jìn)行EMD分解,結(jié)果如圖6a所示,對(duì)該分解結(jié)果進(jìn)行主成分分析(PCA),結(jié)果如圖6b所示,受篇幅限制,只展示圖5中左側(cè)信號(hào)的處理結(jié)果。因?yàn)楣ゎl干擾噪聲的幅值和能量均較大,因此它可能被分解到不同階的IMF分量中(出現(xiàn)EMD分解的混頻現(xiàn)象),對(duì)這些IMF分量進(jìn)行主成分分析后發(fā)現(xiàn),工頻干擾噪聲主要包含在前幾個(gè)主要的主成分之中。因此如果需要提取噪聲,我們只需要對(duì)前幾個(gè)主成分進(jìn)行處理。正如我們?cè)讷@得觀測(cè)信號(hào)后無(wú)法確定其源信號(hào)數(shù)目一樣,確定有效主成分的個(gè)數(shù)非常困難。根據(jù)本文方法的思路和獨(dú)立分量分析方法的前提條件,如果我們確定了待處理的有效主成分個(gè)數(shù),也就確定了盲源分離問題中的源信號(hào)數(shù)目。

圖6 對(duì)圖5c所示的含噪信號(hào)進(jìn)行EMD分解(a)和PCA(b)的結(jié)果

我們根據(jù)參考文獻(xiàn)[18]利用衰減率確定源信號(hào)的數(shù)目,進(jìn)而確定待處理的有效主成分個(gè)數(shù)。計(jì)算由IMF分量所構(gòu)成的矩陣的特征值及衰減率,步驟如下:假設(shè)IMF分量所構(gòu)成的矩陣為X=[IIMF1,IIMF2,…,IIMFN],對(duì)其進(jìn)行奇異值分解,得到的特征值為λ1,λ2,…,λM,衰減率為vi=λi/λi+1(i=1,2,…,M-1),待處理的有效主成分的個(gè)數(shù)N計(jì)算公式如下:

令(vmax,k)=max{v1,v2,…,vM-1},則N=k+1

(3)

式中:vmax為衰減率最大值;k為衰減率最大值的序號(hào)。

為了驗(yàn)證此方法有效性,將圖5中3組含噪信號(hào)分別按上述方法進(jìn)行處理,得到的衰減率如表1所示。

表1 圖5中3組含噪信號(hào)處理后得到的衰減率

此處k=1,即v1值最大,故得到的有效主成分個(gè)數(shù)N為2,與理論情況一致。如果觀測(cè)信號(hào)主要為噪聲,則有用信號(hào)能量很弱,通常這種情況下壓制噪聲幾乎不可能,因此要求工頻干擾噪聲幅值不得超過(guò)原始有用信號(hào)的10倍。

本文進(jìn)行了多次實(shí)驗(yàn),如將50,150,250Hz的正弦波分別添加或一起添加到原始信號(hào)中,實(shí)驗(yàn)結(jié)果均證明利用上述方法可以確定源信號(hào)的數(shù)目。因?yàn)樵夹盘?hào)是隨機(jī)構(gòu)造的且上述方法相關(guān)理論知識(shí)較為成熟,故采用該方法確定的源信號(hào)數(shù)目有效。

得到待處理的有效主成分個(gè)數(shù)N后,我們可以只選取PCA處理后的前N個(gè)主成分進(jìn)行ICA處理,前述處理得到N=2,故本文選取2個(gè)主成分(記為p1,p2)進(jìn)行ICA處理。圖7a、圖7b和圖7c分別為原始信號(hào)、工頻干擾信號(hào)和加噪后信號(hào);圖7d為對(duì)加噪后的信號(hào)進(jìn)行ICA處理后的結(jié)果?？梢园l(fā)現(xiàn),在原始信號(hào)被干擾噪聲完全淹沒的情況下,原始信號(hào)波形已被較好地提取出來(lái)(如圖7d中C1所示),但C1的幅值與原始信號(hào)存在差距。根據(jù)觀測(cè)信號(hào)的幅值范圍確定權(quán)重因子[10],并將權(quán)重因子作用于C1,最終得到去噪后的信號(hào)。

圖7 對(duì)含噪信號(hào)進(jìn)行ICA處理的結(jié)果a 原始信號(hào); b 工頻干擾信號(hào); c 加噪后信號(hào); d 對(duì)加噪后的信號(hào)進(jìn)行ICA處理后的結(jié)果

圖8a、圖8b和圖8c分別為原始信號(hào)、加噪后信號(hào)和去噪后的信號(hào)。比較去噪后的信號(hào)與原始信號(hào)可以看出,提取出幅值較大的工頻干擾噪聲后,信號(hào)波形和幅值均得到了較好的恢復(fù)。

圖8 對(duì)含噪信號(hào)進(jìn)行去噪后的結(jié)果a 原始信號(hào); b 加噪后信號(hào); c 去噪后的信號(hào)

因?yàn)镋MD算法的卓越性能和PCA分解能量的優(yōu)化操作,所以對(duì)EMD之后的IMF直接進(jìn)行PCA處理,并提取第一個(gè)主成分作為消噪后的信號(hào),也能達(dá)到一定的工頻干擾壓制效果。為了方便對(duì)比分析,本文將這種方法稱為PCA方法。表2為噪聲信號(hào)與原始信號(hào)的幅值之比R分別取5種不同值的情況下,3種方法(小波方法、PCA方法、本文方法)的性能。為了比較提取的恢復(fù)信號(hào)與原始信號(hào)的相似程度,比較去噪效果,我們采用相關(guān)系數(shù)來(lái)度量去噪后信號(hào)與原始信號(hào)(無(wú)噪信號(hào))的相似程度,定義如下:

(4)

表2 不同R值下3種方法得到的相關(guān)系數(shù)

從表2中可以發(fā)現(xiàn),隨著工頻干擾噪聲幅值的不斷增加,采用本文方法提取的去噪信號(hào)與原始信號(hào)相關(guān)系數(shù)始終維持在一個(gè)較高的水平(大于0.84),而采用小波方法得到的結(jié)果始終維持在一個(gè)較低的水平(小于0.18),這說(shuō)明本文方法壓制信號(hào)中的工頻干擾噪聲明顯優(yōu)于小波分析方法。PCA方法在R值較小時(shí),性能較優(yōu)越,但當(dāng)噪聲信號(hào)與原始信號(hào)的幅值之比大于1后,其性能迅速下降,即PCA方法穩(wěn)定性欠佳,而大地電磁中工頻干擾噪聲幅值一般都較大,甚至是原始有效信號(hào)的數(shù)倍,所以PCA方法的實(shí)用性不高。

4 實(shí)際大地電磁數(shù)據(jù)處理

選取西部某地采集的大地電磁測(cè)點(diǎn)信號(hào)作為研究對(duì)象,因該地區(qū)電力線信號(hào)、通訊信號(hào)等均非常微弱,故可以認(rèn)為采集到的大地電磁測(cè)點(diǎn)信號(hào)是不包含工頻干擾噪聲的原始信號(hào)。此測(cè)點(diǎn)信號(hào)的視電阻率曲線和相位曲線如圖9所示,曲線形態(tài)平滑、流暢,也從側(cè)面印證此測(cè)點(diǎn)為未被噪聲污染的理想測(cè)點(diǎn)。

圖9 測(cè)點(diǎn)信號(hào)的視電阻率和相位曲線a 高頻結(jié)果(MTH); b 低頻結(jié)果(MTL)

此測(cè)點(diǎn)原始信號(hào)的電場(chǎng)和磁場(chǎng)時(shí)間序列如圖10a 所示(考慮到篇幅限制這里僅展示1000個(gè)采樣點(diǎn)),對(duì)該測(cè)點(diǎn)的5道信號(hào)分別添加50Hz工頻干擾噪聲(幅值為10000),結(jié)果如圖10b所示,可以看出,添加工頻干擾噪聲后,原始信號(hào)被完全淹沒,幾乎完全不可見。此時(shí)若直接進(jìn)行數(shù)據(jù)處理,視電阻率曲線和阻抗相位曲線均會(huì)在50Hz附近出現(xiàn)畸變現(xiàn)象,進(jìn)而影響后續(xù)反演解釋工作,因此必須壓制工頻干擾噪聲。

圖10 測(cè)點(diǎn)信號(hào)的電場(chǎng)和磁場(chǎng)時(shí)間序列a 原始信號(hào); b 添加工頻干擾噪聲后的信號(hào)

采用本文方法進(jìn)行工頻干擾噪聲壓制,得到了整個(gè)測(cè)點(diǎn)的工頻干擾噪聲壓制結(jié)果,本文限于篇幅,僅選取其中1000個(gè)采樣點(diǎn)展示處理結(jié)果(圖11)。圖11的各個(gè)子圖中,從上到下依次為原始信號(hào)、添加了工頻干擾噪聲的信號(hào)、提取的工頻干擾噪聲、壓制工頻干擾噪聲后的信號(hào)(即消噪信號(hào))?？梢钥闯?原始信號(hào)添加工頻干擾噪聲后被工頻干擾噪聲完全淹沒,幾乎看不出任何原始信號(hào)的信息。但經(jīng)過(guò)本文方法處理后,原始信號(hào)的概貌基本得到了恢復(fù)。

圖11 采用本文方法得到的5道信號(hào)處理結(jié)果a Ex; b Ey; c Hx; d Hy; e Hz

表3列出了該測(cè)點(diǎn)(整個(gè)時(shí)間序列)5道信號(hào)的提取噪聲與添加噪聲的相關(guān)系數(shù)、消噪信號(hào)與原始信號(hào)的相關(guān)系數(shù),可以看出,本文方法對(duì)大地電磁中工頻干擾噪聲的提取效率非常高,5道相關(guān)系數(shù)均達(dá)到了0.96,接近于1。從表3還可以看出:消噪信號(hào)與原始信號(hào)的相關(guān)系數(shù)雖然也較高,但無(wú)法達(dá)到0.93。這是因?yàn)榇蟮仉姶判盘?hào)屬于不穩(wěn)定的隨機(jī)信號(hào),添加大幅值的工頻干擾噪聲后,對(duì)原始信號(hào)能量及幅值均有較大影響,故目前還難以完全恢復(fù)出原始信號(hào)。

表3 5道信號(hào)的相關(guān)系數(shù)

添加工頻干擾噪聲(幅值以萬(wàn)為計(jì)數(shù)單位)后以及壓制工頻干擾噪聲后得到的視電阻率曲線和相位曲線分別如圖12所示(因MTL在添加噪聲前后幾乎沒有變化,故這里只展示MTH在去噪后的結(jié)果)。

圖12a為含工頻干擾噪聲大地電磁信號(hào)的視電阻率曲線和相位曲線,可以發(fā)現(xiàn),視電阻率曲線在50Hz(橫坐標(biāo)值約為1.7)附近存在畸變,在5Hz(橫坐標(biāo)值約為0.7)附近存在較大畸變,阻抗相位曲線在50Hz和5Hz附近均存在較大畸變。圖12b為本文方法處理結(jié)果,本文方法使視電阻率曲線和相位曲線在50Hz附近的偏移獲得了較大的改善。因?yàn)樘砑拥墓ゎl干擾噪聲幅值較大,采用本文方法能將原始測(cè)點(diǎn)信號(hào)的視電阻率曲線和相位曲線恢復(fù)至現(xiàn)在的程度,因此說(shuō)明本文方法在壓制工頻干擾噪聲方面具有較好的效果。不足的是,本文方法無(wú)法同時(shí)使5Hz附近曲線平滑,如何對(duì)5Hz附近信號(hào)進(jìn)行有效處理,是今后的研究方向。

圖12 采用本文方法去噪前(a)、后(b)的視電阻率曲線和相位曲線

5 結(jié)論

本文提出基于經(jīng)驗(yàn)?zāi)B(tài)分解和盲源分離的大地電磁工頻干擾噪聲壓制方法,擺脫了小波分析等經(jīng)典方法對(duì)小波基函數(shù)和分解層次的限定和人工經(jīng)驗(yàn)的制約,充分利用了經(jīng)驗(yàn)?zāi)B(tài)分解和盲源分離的優(yōu)良特性。該方法在利用衰減率確定出源信號(hào)個(gè)數(shù)的同時(shí),利用PCA處理EMD分解后的IMF,而后采用ICA處理PCA提取的若干有效主成分,最終得到了壓制工頻干擾噪聲之后的去噪信號(hào)。該方法在工頻干擾噪聲的幅值為原始信號(hào)幅值的數(shù)倍的情況下能在一定程度上壓制工頻干擾噪聲,魯棒性較強(qiáng)。實(shí)驗(yàn)研究發(fā)現(xiàn),工頻干擾噪聲信號(hào)與原始信號(hào)的幅值比為4～7,本文算法的去噪性能最佳,未來(lái)將著眼于拓展算法的應(yīng)用范圍,并結(jié)合人工智能技術(shù),將此方法進(jìn)一步推向智能化。

致謝:感謝長(zhǎng)江大學(xué)物探重點(diǎn)實(shí)驗(yàn)室為本研究提供的野外實(shí)際大地電磁觀測(cè)數(shù)據(jù)。