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

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

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

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

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

1 CEEMDAN-DFA濾波方法

1.1 CEEMDAN

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

(1)

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

1.2 DFA

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

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

(2)

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

1.3 仿真試驗

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

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

圖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、歸一化相關系數、均方根誤差3個指標參數對去噪效果進行評價，對比結果如表1所示。其中，在小波閾值方法中選擇sym6和db6小波函數，分解層數為5層，閾值方式選用sqtwolog，閾值函數選用軟閾值。分析表1可知，當含噪信號的SNR為20dB時，3種方法的去噪效果相差無幾；然而，隨著SNR的降低，經CEEMDAN-DFA方法的歸一化相關系數更接近1，均方根誤差最小。由此可見，在不同噪聲強度背景下，上述3種方法中CEEMDAN-DFA方法的去噪效果是最好、最穩定的。

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

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

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

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

2.1 模糊熵

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

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

1)對于給定的參數m，按順序對時間序列進行相空間重構得到m維向量：

(3)

(4)

(5)

4)定義函數：

(6)

5)令m=m+1，重復步驟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 短時能量

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

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

(9)

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

2.3 FCM聚類算法

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

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

(10)

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

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

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

(11)

3)更新隸屬度矩陣U：

(12)

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

2.4 仿真試驗

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

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

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

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

圖5 FCM模糊聚類結果Fig.5 Fuzzy clustering results by FCM

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

3 實測數據處理

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

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

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

圖6 實測數據傳統整體濾波方法效果Fig.6 The effect of traditional integral filtering method for measured data

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

圖7 實測數據傳統整體濾波方法與CEEMDAN-DFA與FCM聚類算法的噪聲識別及濾波結果對比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為經上述2種方法處理后所得結果的時頻分布對比圖，可以看出，經整體濾波處理后的結果在低頻部分幅值出現快速下降的現象，說明該方法損失了有用的低頻信息；而經筆者所提方法處理后的結果則較好地保持原始信號中低頻部分能量分布的特征，表明該方法能更多地保留反映低頻信息的細節成分。

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

4 結論

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

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

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