基于CEEMDAN-DFA與FCM聚類算法的大地電磁強噪聲識別與抑制

史維，嚴良俊，謝興兵，周磊

1.油氣資源與勘探技術(shù)教育部重點實驗室(長江大學)，湖北 武漢 430100 2.非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心(長江大學)，湖北 武漢 430100 3.長江大學工程技術(shù)學院，湖北 荊州 434020

大地電磁測深(magnetotelluric sounding，簡稱MT) 法是一種通過研究大地對天然交變電磁場的頻率響應(yīng)來獲取地下不同深度介質(zhì)電性特征分布的電磁勘探方法。MT法作為一種重要的地球物理勘探方法，現(xiàn)已被廣泛應(yīng)用于油氣田勘探、深部找礦、地熱資源調(diào)查、地震監(jiān)測、深部地質(zhì)構(gòu)造研究等眾多領(lǐng)域[1-3]。MT法以天然場為場源，天然電磁場具有能量弱、幅度變化大、頻帶范圍寬等特點，野外實測信號常常容易受到各類噪聲的干擾[4，5]。隨著人類社會工業(yè)化進程的加快，噪聲也越來越復(fù)雜，特別是強干擾地區(qū)(如礦山、城區(qū))附近，噪聲的幅值往往是天然信號的幾倍甚至是幾個數(shù)量級，導(dǎo)致采集數(shù)據(jù)的質(zhì)量明顯下降，并且很難將噪聲與天然信號分離，致使后續(xù)的處理結(jié)果不能真實地反映地下介質(zhì)的電性分布特征，從而影響整個勘探效果。因此，MT信號中的噪聲抑制問題長期以來都是廣大國內(nèi)外學者長期關(guān)注和研究的熱點之一。

近年來，如小波變換[6，7]、經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)[8，9]、數(shù)學形態(tài)濾波[10]、壓縮感知重構(gòu)算法[11]、局域均值分解[12]等現(xiàn)代信號處理方法均被引入到MT信號噪聲壓制領(lǐng)域，對于MT數(shù)據(jù)質(zhì)量的改善都取得了一定的效果。但是上述方法均是對含噪的MT數(shù)據(jù)時間片段進行整體濾波處理，缺少噪聲識別環(huán)節(jié)，在壓制噪聲的同時會損失一部分緩慢變化的低頻信息，導(dǎo)致過處理現(xiàn)象的產(chǎn)生。鑒于此，筆者嘗試對含強噪聲的MT數(shù)據(jù)運用CEEMDAN(complete ensemble empirical mode decomposition with adaptive noise，自適應(yīng)噪聲的完備經(jīng)驗?zāi)B(tài)分解)-DFA(detrended fluctuation analysis，去趨勢波動分析)與模糊聚類算法相結(jié)合的技術(shù)進行噪聲識別與噪聲抑制來解決該問題，以期更好地改善實測數(shù)據(jù)質(zhì)量。

1 CEEMDAN-DFA濾波方法

1.1 CEEMDAN

CEEMDAN是在EMD和集合經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)基礎(chǔ)上發(fā)展而來的。EMD不需要選擇基函數(shù)，能夠根據(jù)信號的時間特征自適應(yīng)地將復(fù)雜信號分解為有限個頻率由高到低分布的本征模態(tài)函數(shù)(intrinsic mode function，IMF)，是一種分析和處理非線性、非平穩(wěn)信號的有效方法[13]。然而，當信號中存在間歇干擾和噪聲時，EMD容易產(chǎn)生模態(tài)混疊現(xiàn)象，為了解決該問題，2009年WU和HUANG提出了EEMD方法[14]。EEMD雖然在一定程度上解決了模態(tài)混疊問題，但其弊端在于信號在分解過程中不能完全消除殘余噪聲、重構(gòu)信號誤差較大。在EEMD的基礎(chǔ)上，CEEMDAN被TORRES等提出[15]，該方法在EMD分解的每一階段添加自適應(yīng)高斯白噪聲，計算唯一的余量來獲取各個模態(tài)分量，相較于EEMD，它具有更好的分解完備性、重構(gòu)誤差幾乎為零、計算量小等優(yōu)點。信號x(n)經(jīng)CEEMDAN分解后，可以被精確地分解為k個IMF分量和1個殘差r(n)，即：

(1)

式中：xi(n)為IMF分量；i=1,2,…,k。

1.2 DFA

在EMD和其改進的去噪算法中，如何判斷某個IMF分量是信號主導(dǎo)還是噪聲主導(dǎo)，是影響算法去噪效果的關(guān)鍵因素。通過人工判定的方式選擇，需要一定的先驗知識，且缺乏自適應(yīng)性。DFA利用分形特性對信號的復(fù)雜性進行量化，提供了一個定量參數(shù)(標度指數(shù)α)來表示非平穩(wěn)時間序列的自相關(guān)特性，是一種被成功運用于檢測非平穩(wěn)信號的長程相關(guān)屬性的方法[16]。其計算步驟如下：

2)將y(k)不重疊地分為n個長度相等的窗口,對每一窗口，根據(jù)最小二乘原理，利用多項式擬合局部趨勢項yn(k)來計算波動均方根F(n)：

(2)

標度指數(shù)α可以作為識別噪聲IMF的判別依據(jù)，具有可靠的評判標準[16]。當α=0.5時，說明時間序列不相關(guān)，如白噪聲；當0<α<0.5時，表明時間序列呈反向相關(guān)，即小的波動后緊接著會出現(xiàn)大的波動，反之亦然；當0.5<α<1時，說明信號呈現(xiàn)相關(guān)性；當α≥1時，相關(guān)性不具有冪律性。因此，當α<0.5時，IMF的自相關(guān)性差，可以認為它是噪聲；當α>0.5時，IMF具有長程相關(guān)性，即為有用信號。

1.3 仿真試驗

MERT等[17]利用EMD-DFA對不同信噪比下的合成信號和真實信號進行去噪處理，效果優(yōu)于軟、硬小波閾值法。在此基礎(chǔ)上，筆者采用CEEMDAN-DFA方法對含噪的MT信號進行濾波，選取小波分析中常用的piecewise-regular信號為原始信號，其采樣率為4096Hz，時間長度為1s，添加信噪比(signal-to-noise ratio，SNR)為10dB的高斯白噪聲得到含噪信號。

圖1為含噪信號經(jīng)CEEMDAN分解后所得的IMFs和殘差(residual)，但從圖中無法判斷哪些IMFs是屬于信號分量還是噪聲分量。計算每一階IMF(IMF1～ IMF12)對應(yīng)的標度指數(shù)α，結(jié)果分別為0.386、0.219、0.214、0.308、0.531、0.900、1.450、1.760、1.948、1.986、1.996、2.012?；谥暗姆治?，將計算結(jié)果α>0.5的IMF判定為有用信號，CEEMDAN-DFA濾波結(jié)果應(yīng)為IMF5～IMF12之和，其結(jié)果如圖2(c)所示，此外圖2還給出了小波閾值濾波(db6)和EMD-DFA濾波結(jié)果，分別如圖2(d)、圖2(e)所示，可以看出，相較于小波閾值濾波(db6)和EMD-DFA濾波，經(jīng)CEEMDAN-DFA方法濾波所得的結(jié)果曲線毛刺明顯減少，形態(tài)更加清晰、光滑，更好地還原了原始信號特征。

圖1 含噪信號CEEMDAN分解過程Fig.1 Decomposition of the noise signal by CEEMDAN

圖2 SNR為10dB時3種方法濾波效果Fig.2 The filtering effect by three denosing methods with SNR=10dB

為了更加綜合、客觀地評價上述3種方法所獲得的去噪效果，筆者對原始piecewise-regular信號加入不同強度的噪聲，再分別利用上述3種方法對含噪信號進行去噪處理，并選用SNR、歸一化相關(guān)系數(shù)、均方根誤差3個指標參數(shù)對去噪效果進行評價，對比結(jié)果如表1所示。其中，在小波閾值方法中選擇sym6和db6小波函數(shù)，分解層數(shù)為5層，閾值方式選用sqtwolog，閾值函數(shù)選用軟閾值。分析表1可知，當含噪信號的SNR為20dB時，3種方法的去噪效果相差無幾；然而，隨著SNR的降低，經(jīng)CEEMDAN-DFA方法的歸一化相關(guān)系數(shù)更接近1，均方根誤差最小。由此可見，在不同噪聲強度背景下，上述3種方法中CEEMDAN-DFA方法的去噪效果是最好、最穩(wěn)定的。

表1 不同噪聲強度背景下3種方法去噪效果對比

2 基于模糊C均值聚類算法的噪聲識別

模糊聚類算法是一種典型的無監(jiān)督學習方法，它以模糊集合論為數(shù)學基礎(chǔ)，其基本思想是將待分類的對象按照數(shù)據(jù)的特征分成若干類，使得劃分為同一類的對象之間的相似度最大，而不同類之間對象的相似度最小[18]。模糊C均值(fuzzy C-means，F(xiàn)CM)聚類是模糊聚類中應(yīng)用最為廣泛的一種算法，該算法理論完備，聚類效果良好，已在圖像處理[19]、機械故障診斷[20，21]等領(lǐng)域得到廣泛應(yīng)用。

特征參數(shù)的選定是數(shù)據(jù)聚類分析的首要問題，選擇哪些特征參數(shù)作為聚類特征來使用，會直接影響到聚類效果。通過分析大量的實測MT數(shù)據(jù)可知，常見的強干擾有方波噪聲、脈沖噪聲、三角波噪聲等，相較于隨機的天然電磁場信號，這些噪聲具有能量強、幅值大、形態(tài)特征明顯等特點。鑒于此，筆者嘗試從時間序列的復(fù)雜度與能量大小2個角度出發(fā)，選取模糊熵與短時能量2個特征參數(shù)作為聚類特征，并利用FCM聚類算法對實測MT數(shù)據(jù)中受強噪聲干擾的部分進行識別。

2.1 模糊熵

模糊熵是陳偉婷等[22]在近似熵和樣本熵概念的基礎(chǔ)上提出的一種新方法，這3者具有類似的物理意義，都能衡量時間序列維數(shù)變化時產(chǎn)生新模式的概率大小。概率越大，則時間序列的復(fù)雜度越大，其熵值也越大；反之亦然。然而，相較于近似熵和樣本熵，模糊熵具有更好的連續(xù)性、一致性和數(shù)據(jù)長度獨立性，且不受基線漂移影響的特點[23]。

假定長度為N的時間序列為{u(i)：1≤i≤N}，模糊熵的計算流程如下：

1)對于給定的參數(shù)m，按順序?qū)r間序列進行相空間重構(gòu)得到m維向量：

(3)

(4)

(5)

4)定義函數(shù)：

(6)

5)令m=m+1，重復(fù)步驟1)～4)可得φm+1(n,r)：

(7)

6)模糊熵FE定義為：

FE(m,n,r,N)=lnφm(n,r)-lnφm+1(n,r)

(8)

2.2 短時能量

能量分析用來反映信號能量隨時間變化的情況，短時能量為一個短段(幀)的能量。目前，短時能量被作為一種識別語音信號與非語音信號的重要參數(shù)在語音信號處理中已得到了廣泛應(yīng)用[24]。鑒于MT信號與語音信號都屬于非線性、非平穩(wěn)信號，且天然MT信號和強噪聲在能量特征上有顯著區(qū)別。因此，選取短時能量作為模糊聚類算法的特征參數(shù)去識別MT有用信號和強噪聲在一定程度上是可行的。

對MT時間序列進行不重疊分幀，第i幀短時能量E(i)公式為：

(9)

式中：ti(n)為第i幀時間序列，1≤i≤M；M為總幀數(shù)；L為每幀的長度。

2.3 FCM聚類算法

FCM聚類是傳統(tǒng)硬聚類算法的改進，其基本思想是計算從所有數(shù)據(jù)到每個聚類中心的歐氏距離及模糊隸屬度的加權(quán)和所確定的目標函數(shù)，通過反復(fù)修改聚類中心矩陣和隸屬度矩陣使目標函數(shù)最小化，以達到將具有相似特征的數(shù)據(jù)對象聚為一類的目的[20]。

FCM的目標函數(shù)J(U,V)為：

(10)

從式(10)可以看出，目標函數(shù)J(U,V)等于加權(quán)類內(nèi)距離平方和，其值越小，表明樣本離某個聚類中心越近。因此，F(xiàn)CM聚類算法的實質(zhì)是使目標函數(shù)J(U,V)最小化來求解隸屬度矩陣U和聚類中心矩陣V，具體步驟如下：

1)給定聚類個數(shù)c和加權(quán)指數(shù)m，初始化隸屬度矩陣U=[uij]c×n，設(shè)定精度ε，令迭代次數(shù)l=0。

2)計算樣本的c個聚類中心vi：

(11)

3)更新隸屬度矩陣U：

(12)

4)當滿足‖Ul+1-Ul‖<ε時，迭代停止，否則重新執(zhí)行步驟2)、3)直到滿足精度為止。

2.4 仿真試驗

先選取內(nèi)蒙古錫林浩特某工區(qū)幾乎未受干擾的MT信號中的30個時間片段，然后再從受人為干擾嚴重影響的湖南某工區(qū)選取90個受典型強噪聲干擾的時間片段(含矩形波噪聲、含脈沖噪聲及含三角波噪聲干擾的時間片段各30個)，總共120個時間片段組成測試樣本集。其中每個時間片段均包含150個采樣點，從每種類型的信號中隨機選取一個片段，波形如圖3所示。

圖3 各種實測MT信號時間片段Fig.3 Time domain waveform segments of various measured MT signals

計算120個樣本的模糊熵及歸一化短時能量，參數(shù)數(shù)值分布如圖4所示。由模糊熵的定義可知，模糊熵的大小能反映時間序列的復(fù)雜程度。從圖4可以看出，無干擾MT信號的模糊熵要大于受到強噪聲干擾的3類MT信號的模糊熵，能很好地將無干擾和受強干擾的MT信號區(qū)分開，其原因是：無干擾的MT信號屬于非線性、非平穩(wěn)的隨機信號，波形無規(guī)則程度較高，相對較復(fù)雜，故模糊熵較大，而受嚴重干擾的MT信號由于受到強噪聲的影響，形態(tài)比較規(guī)則，復(fù)雜度較低，故模糊熵較小。由于實測MT數(shù)據(jù)幅值容易受到地域、基線漂移等因素的影響，從而導(dǎo)致數(shù)據(jù)幅值的數(shù)量級差異巨大，使得短時能量無法客觀地反映信號能量隨時間變化的規(guī)律。為避免上述因素所帶來的影響，在計算短時能量之前需要對數(shù)據(jù)進行預(yù)處理，筆者采用去趨勢項、標準差標準化的預(yù)處理手段使數(shù)據(jù)在同一尺度范圍下進行計算分析，經(jīng)預(yù)處理后計算的歸一化短時能量結(jié)果如圖4(b)所示。從圖4(b)中可以看出，雖然受強干擾的3類MT信號曲線跳變、交叉現(xiàn)象比較明顯，但與無干擾的MT信號并沒有混疊，這表明歸一化短時能量對無干擾和受強干擾的MT信號也有較好的區(qū)分度。

圖4 特征參數(shù)數(shù)值分布Fig.4 Value distribution of characteristic parameters

圖5 FCM模糊聚類結(jié)果Fig.5 Fuzzy clustering results by FCM

將模糊熵和歸一化短時能量作為特征參數(shù)，并利用FCM算法進行模糊聚類，結(jié)果如圖5所示，可以看出，30個未受干擾的MT信號樣本分為一類，其他90個受強噪聲干擾的MT信號樣本則分為另一類，顯現(xiàn)了良好的聚類效果。

3 實測數(shù)據(jù)處理

實測數(shù)據(jù)來源于湖南某地區(qū)受干擾嚴重的TYS2-31A測點的采集資料，采集儀器為加拿大Phoenix公司的MTU-5A。數(shù)據(jù)采樣率為高、中、低3種頻率，分別為2400、150、15Hz，采集數(shù)據(jù)分別保存放在TS3、TS4、TS5這3個文件中，每個文件均包含5個分量(Ex、Ey、Hx、Hy、Hz)。考慮到數(shù)據(jù)量巨大，為了便于解釋說明與對比，僅選取該測點TS5文件的電場分量Ex和磁場分量Hy中受到矩形波、三角波干擾的2個時間片段，并利用傳統(tǒng)的整體濾波和筆者所提方法分別進行處理，對其他噪聲的處理方法也與之類似。

傳統(tǒng)的整體濾波方法首先會對含噪的MT信號進行CEEMDAN-DFA濾波，提取強噪聲整體輪廓，然后從含噪的MT信號中將強噪聲剔除，得到重構(gòu)的有用MT信號，結(jié)果如圖6所示。分析圖6可知，利用CEEMDAN-DFA方法進行整體濾波雖然可以較好地壓制矩形波和三角波干擾所產(chǎn)生的影響，但重構(gòu)后的MT信號僅在基線附近振動，原始信號表現(xiàn)出的整體波動趨勢已不存在，表明原始MT信號中緩慢變化的低頻信息在整體濾波過程中已被剔除，即出現(xiàn)了過處理現(xiàn)象，使得后續(xù)卡尼亞視電阻率計算結(jié)果中會出現(xiàn)低頻信息的缺失，從而導(dǎo)致后續(xù)的解釋結(jié)果無法客觀地反映地下介質(zhì)深部構(gòu)造的電性特征。

筆者在傳統(tǒng)整體濾波方法的基礎(chǔ)上增加了強噪聲干擾識別環(huán)節(jié)，使得濾波更具針對性，處理步驟如下：①對受干擾信號先進行整體濾波，保留濾波結(jié)果；②對原始信號進行不重疊分幀，每一幀數(shù)據(jù)為一個樣本，樣本長度為10s，即150個樣點，并對數(shù)據(jù)進行預(yù)處理；③計算每個樣本的模糊熵和歸一化短時能量；④利用FCM算法對樣本進行分類識別；⑤將識別為受到強干擾的數(shù)據(jù)段所對應(yīng)的整體濾波結(jié)果與識別為未受干擾的MT信號組合作為最終的處理結(jié)果。

圖6 實測數(shù)據(jù)傳統(tǒng)整體濾波方法效果Fig.6 The effect of traditional integral filtering method for measured data

圖7為實測數(shù)據(jù)傳統(tǒng)整體濾波方法與CEEMDAN-DFA與FCM聚類算法的噪聲識別及濾波結(jié)果對比。分析圖7可知，CEEMDAN-DFA與FCM聚類算法可自動將受嚴重干擾和未受干擾的MT信號區(qū)分開，所得結(jié)果不僅能壓制矩形波和三角波干擾的影響，而且能很好地保留原始信號中緩慢變化的成分。

圖7 實測數(shù)據(jù)傳統(tǒng)整體濾波方法與CEEMDAN-DFA與FCM聚類算法的噪聲識別及濾波結(jié)果對比Fig.7 Comparison of the noise recognition and filtering results between the traditional integral filtering method and CEEMDAN-DFA combined with FCM clustering algorithm for measured data

圖8為經(jīng)上述2種方法處理后所得結(jié)果的時頻分布對比圖，可以看出，經(jīng)整體濾波處理后的結(jié)果在低頻部分幅值出現(xiàn)快速下降的現(xiàn)象，說明該方法損失了有用的低頻信息；而經(jīng)筆者所提方法處理后的結(jié)果則較好地保持原始信號中低頻部分能量分布的特征，表明該方法能更多地保留反映低頻信息的細節(jié)成分。

注：lg(Amp)、lg(Fre)表示對幅值和頻率取以10為底對數(shù)。圖8 傳統(tǒng)整體濾波方法與CEEMDAN-DFA與FCM聚類算法的時頻分析結(jié)果對比Fig.8 Comparison of time-frequency analysis results between the traditional integral filtering method and CEEMDAN-DFA combined with FCM clustering algorithm

4 結(jié)論

1)考慮到不同噪聲背景下的去噪效果，選用CEEMDAN-DFA方法對含噪的仿真信號進行濾波處理，并與小波閾值和EMD-DFA的去噪結(jié)果進行了對比分析，結(jié)果表明，CEEMDAN-DFA的濾波效果要整體好于其他2種方法。

2)從時間序列復(fù)雜度和能量特征2個角度出發(fā)，選取模糊熵與短時能量參數(shù)對天然MT信號和受強干擾影響MT信號之間的差異進行定量分析，并在FCM聚類算法中運用上述2個參數(shù)對無干擾和受強干擾影響的MT信號進行聚類識別，獲得了較好的識別效果。

3)實測數(shù)據(jù)處理結(jié)果表明，相較于傳統(tǒng)整體濾波方法，將CEEMDAN-DFA與FCM聚類算法相結(jié)合的方法增加了噪聲識別環(huán)節(jié)，僅對識別為受強噪聲干擾的數(shù)據(jù)進行定向濾波，既達到了去噪的目的又保留了整體濾波處理中損失的低頻成分。該方法對改善強干擾區(qū)的MT數(shù)據(jù)質(zhì)量具有一定的參考價值。