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

曹小玲,陳清禮,蔣 濤

(1.長江大學油氣資源與勘探技術教育部重點實驗室,湖北武漢430100;2.長江大學信息與數學學院,湖北荊州434023;3.長江大學電工電子國家級實驗教學示范中心,湖北荊州434023)

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

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

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

大地電磁天然場源的特點[11-12]決定了大地電磁信號非常容易受到噪聲的干擾,且這些噪聲形態各異、種類繁多、強弱不一。根據噪聲產生的機理劃分,大地電磁測深信號包含的噪聲包括工頻干擾噪聲、儀器噪聲、人文噪聲、地質噪聲等。工頻干擾噪聲一般指電力傳輸系統、人工電磁場、通信過程及有線廣播等產生的噪聲。這些電磁噪聲呈現非平面波的形態,且與觀測點距離較近,因此不符合大地電磁測深對場源的要求。

因工頻干擾噪聲很大一部分是人為導致的,故也有學者認為應將其歸為人文噪聲。如文獻[1]認為大地電磁測深中的噪聲包括場源噪聲、地質噪聲和人文噪聲3個類別,其中人文噪聲主要是指“人類活動所產生的干擾性電磁波以及人類的活動產生的其它電磁噪聲”[1],這里的“人文噪聲”包括了工頻干擾噪聲,即認為工頻干擾噪聲是人文噪聲中的一部分。

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

圖1 包含工頻干擾噪聲的實測大地電磁信號

2 方法原理

2.1 盲源分離基本原理

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

圖2 盲源分離流程

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

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

(1)

根據王書明等[14-15]對不同地區實測大地電磁信號的性質特征分析可知,大地電磁信號符合盲源分離對源信號的統計性質要求。

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

圖3 ICA方法的算法流程

2.2 EMD方法基本原理

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

1) 假設原始信號為x(t),首先我們找出x(t)中所有的局部極值點(包括極大值點和極小值點),用三次樣條函數連接所有的局部極大值點構成上包絡線,用三次樣條函數連接所有的局部極小值點構成下包絡線。

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

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

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

(2)

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

本文沒有選擇EEMD而選擇EMD處理的原因如下:一是因為EEMD是在信號中添加白噪聲后再進行EMD處理,添加白噪聲可能對原有的工頻干擾噪聲產生影響;二是因為大地電磁信號普遍數據量龐大,采用EEMD處理會極大地增加計算量。

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

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

圖4 基于經驗模態分解和盲源分離的大地電磁工頻干擾噪聲壓制方法流程

該方法主要優勢在于EMD模型、PCA處理及盲源分離方法的應用,受文獻[18]的啟發,我們利用衰減率確定源信號的數目和待處理的主成分的個數,解決了欠定盲源分離中源信號個數無法確定的問題,實現了單通道觀測信號的源信號的分離和提取。

PCA處理在盲源分離處理領域是一種極其有效的處理手段,我們對EMD分解得到的IMF分量進行PCA處理[19],可以有效降低信號數據的維數,通過尋找信號中能量最大的分量得到信號的主要特征量,即起到“化繁為簡”的作用。另外,PCA處理作為一種統計分析方法,還能實現對信號的去相關處理。因為它能夠使得變換后產生的新分量正交或不相關,故變換后的分量能量更集中、性質更穩定[20]。PCA處理的優勢在于利用某種變換將數據原有的大量特征變換為幾個主要特征,而這些特征包含了原始數據的最主要的特征信息。因此,當我們僅僅提取前面的部分特征時,可以既減小特征的個數又保留原有信號最主要的特征信息。同時,由于PCA處理后的信號不相關,所以噪聲子空間和源信號子空間相互分離,進而可以達到初步壓制部分噪聲的目的,即初步壓制部分工頻干擾噪聲。

PCA處理雖然能夠去除信號的相關性,但在高階累積量的定義下信號可能仍然具有相關性。為使所得信號具有的統計獨立性最大,我們又進行了ICA處理,選取有用分離向量與混合矩陣相乘重構得到分離后的獨立分量信號,進而得到消噪信號。

3 仿真模擬實驗

如圖5a所示,利用蒙特卡洛方法對仿真信號中的原始信號隨機構造3組信號(幅值為400)。因為工頻干擾噪聲一般為50Hz的周期信號及其諧波,故構造工頻干擾噪聲信號(噪聲信號)可表示為S=800sin(2π×50t)。現實情況下強烈的工頻干擾噪聲主要來自電力傳輸系統,其頻率一般為50Hz,從方便分析的角度考慮,本文設工頻干擾噪聲為50Hz。因工頻干擾噪聲信號幅值一般是有效電磁信號幅值的數倍,所以這里幅值應設較大值,即為原始信號的數倍(圖5b)。圖5c為工頻干擾噪聲加入原始信號后得到的含噪信號。

圖5 3組含噪信號構建a 隨機信號; b 50Hz正弦信號的工頻干擾噪聲; c 含噪信號

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

圖6 對圖5c所示的含噪信號進行EMD分解(a)和PCA(b)的結果

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

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

(3)

式中:vmax為衰減率最大值;k為衰減率最大值的序號。

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

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

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

本文進行了多次實驗,如將50,150,250Hz的正弦波分別添加或一起添加到原始信號中,實驗結果均證明利用上述方法可以確定源信號的數目。因為原始信號是隨機構造的且上述方法相關理論知識較為成熟,故采用該方法確定的源信號數目有效。

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

圖7 對含噪信號進行ICA處理的結果a 原始信號; b 工頻干擾信號; c 加噪后信號; d 對加噪后的信號進行ICA處理后的結果

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

圖8 對含噪信號進行去噪后的結果a 原始信號; b 加噪后信號; c 去噪后的信號

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

(4)

表2 不同R值下3種方法得到的相關系數

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

4 實際大地電磁數據處理

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

圖9 測點信號的視電阻率和相位曲線a 高頻結果(MTH); b 低頻結果(MTL)

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

圖10 測點信號的電場和磁場時間序列a 原始信號; b 添加工頻干擾噪聲后的信號

采用本文方法進行工頻干擾噪聲壓制,得到了整個測點的工頻干擾噪聲壓制結果,本文限于篇幅,僅選取其中1000個采樣點展示處理結果(圖11)。圖11的各個子圖中,從上到下依次為原始信號、添加了工頻干擾噪聲的信號、提取的工頻干擾噪聲、壓制工頻干擾噪聲后的信號(即消噪信號)。可以看出,原始信號添加工頻干擾噪聲后被工頻干擾噪聲完全淹沒,幾乎看不出任何原始信號的信息。但經過本文方法處理后,原始信號的概貌基本得到了恢復。

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

表3列出了該測點(整個時間序列)5道信號的提取噪聲與添加噪聲的相關系數、消噪信號與原始信號的相關系數,可以看出,本文方法對大地電磁中工頻干擾噪聲的提取效率非常高,5道相關系數均達到了0.96,接近于1。從表3還可以看出:消噪信號與原始信號的相關系數雖然也較高,但無法達到0.93。這是因為大地電磁信號屬于不穩定的隨機信號,添加大幅值的工頻干擾噪聲后,對原始信號能量及幅值均有較大影響,故目前還難以完全恢復出原始信號。

表3 5道信號的相關系數

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

圖12a為含工頻干擾噪聲大地電磁信號的視電阻率曲線和相位曲線,可以發現,視電阻率曲線在50Hz(橫坐標值約為1.7)附近存在畸變,在5Hz(橫坐標值約為0.7)附近存在較大畸變,阻抗相位曲線在50Hz和5Hz附近均存在較大畸變。圖12b為本文方法處理結果,本文方法使視電阻率曲線和相位曲線在50Hz附近的偏移獲得了較大的改善。因為添加的工頻干擾噪聲幅值較大,采用本文方法能將原始測點信號的視電阻率曲線和相位曲線恢復至現在的程度,因此說明本文方法在壓制工頻干擾噪聲方面具有較好的效果。不足的是,本文方法無法同時使5Hz附近曲線平滑,如何對5Hz附近信號進行有效處理,是今后的研究方向。

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

5 結論

本文提出基于經驗模態分解和盲源分離的大地電磁工頻干擾噪聲壓制方法,擺脫了小波分析等經典方法對小波基函數和分解層次的限定和人工經驗的制約,充分利用了經驗模態分解和盲源分離的優良特性。該方法在利用衰減率確定出源信號個數的同時,利用PCA處理EMD分解后的IMF,而后采用ICA處理PCA提取的若干有效主成分,最終得到了壓制工頻干擾噪聲之后的去噪信號。該方法在工頻干擾噪聲的幅值為原始信號幅值的數倍的情況下能在一定程度上壓制工頻干擾噪聲,魯棒性較強。實驗研究發現,工頻干擾噪聲信號與原始信號的幅值比為4～7,本文算法的去噪性能最佳,未來將著眼于拓展算法的應用范圍,并結合人工智能技術,將此方法進一步推向智能化。

致謝:感謝長江大學物探重點實驗室為本研究提供的野外實際大地電磁觀測數據。