基于CEEMD-排列熵的循環策略信號提取方法*

尹昱東 明 勇 邊 羽

(①西安交通大學機械工程學院，陜西 西安 710049;②成都理工大學地球物理學院，四川 成都 610009;③成都開放大學，四川 成都 610051;④內蒙古大學電子信息工程學院，內蒙古 呼和浩特 010000)

目前，非線性科學已迅速發展涵蓋到了各個學科領域，其中，尤以混沌時間序列的研究最為突出，占據了很大的比重[1]，混沌預測是非線性科學領域的一個熱門課題，已經應用到短時交通流、深部巖體變形、風電場短期風速預測、短期電力負荷預測及海雜波中雷達目標信號提取[2]等多個方面，由于混沌時間序列具有較強的非平穩性、突變性等特點，如何建立高精度的逼近混沌系統的預測模型本身就是一項難題，若是實測信號受到噪聲的干擾，勢必導致混沌預測難上加難，干擾與混沌時間序列本質特征相互耦合和纏繞，掩蓋了系統的內在動態特性，這為系統狀態分析和參數預測帶來了困難，并且混沌信號具極強的似噪聲性，其頻譜是連續的寬譜，導致混沌信號和噪聲會在相同或相近頻率段重合，致使傳統的線性去噪方法和頻譜分析方法難以將其分離[3]。

目前常用的非線性非平穩信號降噪方法是Huang N E等[4]提出的經驗模態分解(empirical mode decomposition，EMD)及其相關延伸改進算法，已被廣泛應用于多個領域，但尚存在以下問題：①EMD由于異常事件的干擾，使分解得到的固有模態分量(intrinsic mode function，IMF)中含有不同尺度的信號，或是相近尺度信號分布在不同固有模態分量中，即模式混疊，增加了信號重構誤差；②總體平均經驗模態分解[5](ensemble empirical mode decomposition，EEMD)通過多次疊加和抵消白噪聲減小模式混疊，但是EEMD預處理過程耗時過長，時效性差；③EEMD固有模態分量中有用信號的特征提取方法在噪聲未知情況下難以取得理想效果[6]。針對以上問題，國內外學者進行了大量研究，文獻[7]對輸入信號進行B樣條最小二乘擬合，改善時間尺度分布，弱化異常事件對信號的影響，達到抑制EMD模態混疊的作用，但是對于信號提取方法的處理方面尚未提及；文獻[8]在EMD閾值處理時，根據IMF間功率譜密度關系確定噪聲強度，然后采用硬閾值處理方法，忽略了硬閾值帶來的不連續性問題；文獻[9]在使用EEMD對慣性導航信號降噪時，提出了改進的區間閾值降噪方法，合理設置調節因子能兼顧軟硬閾值函數的優勢，降噪效果比較依賴于調節因子的選擇，并且沒有改善EEMD耗時的問題；Yeh J R等[10]在EEMD基礎上提出了補充的總體平均經驗模態分解方法(complementary ensemble empirical mode decomposition，CEEMD)，通過添加相反白噪聲抵消信號中的部分噪聲，抑制了噪聲引起的模態混疊的同時減小EEMD疊加白噪聲引起的重構誤差；在閾值處理方面，文獻[11]利用CEEMD-排列熵確定噪聲含量較多的IMF分量，然后采用小波閾值降噪方法對含有較多噪聲的IMF分量進行降噪處理以保留有用信息，小波閾值降噪效果較為依賴閾值選取的準確性。

為解決經驗模態分解及其延伸算法在信號特征提取存在的問題，在抑制噪聲的同時，更好保留有用信號，本文提出了基于CEEMD-排列熵的梯度投影稀疏重構降噪方法，原始信號經CEEMD處理，對信號疊加相反白噪聲抑制白噪聲引起重構誤差的同時簡化了計算方法，有效解決了EMD模式混疊和EEMD耗時過長的問題，對分解得到的本征模態分量通過計算排列熵確定噪聲分量和信號分量，考慮到信號中噪聲先驗知識未知，提出了基于奇異值分解的循環策略信號提取方法，利用噪聲的排列熵較大這一特性，合理設置奇異值有效重構階次的選取，該方法無需信號和噪聲的先驗知識，在抑制噪聲的同時，可以較好保留有用信號。

1 算法基本原理

1.1 CEEMD算法

EMD分解是基于以下假設：信號至少有兩個極值點：一個最大值和一個最小值；數據的局部時域特性是由極值間隔來確定；若數據缺乏極值但是含有拐點，則可以通過對數據微分一次或多次求取極值點，然后再通過積分獲取分解結果。實際采集的振動信號數據較長，在時域范圍內含有多個局部極值點，滿足上述假設條件，所以EMD已被廣泛應用于非線性振動信號處理。實測的振動信號通過EMD篩分過程分解為一系列頻率由高到低排列的固有模態分量和一個余項。

CEEMD是在EMD基礎上的改進，同時可以有效減小EEMD疊加的白噪聲無法完全抵消引起的重構誤差及運算時間過長的問題，具體步驟如下：

步驟1：對含噪原始信號疊加一組白噪聲序列n1(t)，將加噪信號經EMD處理得到一組固有模態分量c1i和余項r1。

步驟3：重復步驟1-2，得到n組cni、rn、c-ni、r-n。

步驟4：對步驟3的多組IMF分量進行平均組合得到分解后最終的IMF。

(1)

(2)

其中:cj(t)表示CEEMD處理最終得到的第j個IMF分量，cni表示信號經EMD處理后得到的n組中的的第j個IMF分量、rn(t)表示CEEMD處理最終得到的余項，ri表示信號經EMD處理后得到的n組中的的第i個余項。

1.2 排列熵

排列熵(permutation entropy, PE)是度量時間序列非邏輯性和復雜度的一種有效方法，PE具有計算簡單、抗干擾能力強等優點，對非線性數據具有較強的魯棒性，其計算方法如下：

對時間序列{x(i),i=1,2,…,N}進行相空間重構

(3)

其中:m為嵌入維數，為時延。

于是，時間序列{x(i),i=1,2,…,N}排列熵以Shannon熵形式表示為：

在記錄詞集中，詞匯可出現的位置有：標題、關鍵詞、成因和答案等，出現在不同的位置對檢索結果的貢獻值不同，所以其重要性也不同。例如，出現在標題中的詞匯的重要性要比出現在答案中的詞匯重要。

(4)

在Pg=1/m時，得到Hp(m)的最大值ln(m!)，因此，可以通過ln(m!)將Hp(m)進行如下歸一化處理：

Hp=Hp(m)/ln(m!)

(5)

Hp在[0，1]之間取值，表示序列的隨機程度，并且其值越大，表示序列隨機性越強，值越小，則時間序列越規則。

文獻[12]列出了噪聲和幾種常見信號的排列熵，如表1所示，可以看出，噪聲的排列熵較大，這與噪聲的隨機不規則特性吻合，正弦信號、調幅信號以及調幅調頻信號等規則性信號的排列熵較小，而間歇性信號相比正弦信號更加不規則，所以排列熵較大，參考文獻[12]的論證結果，本文在對固有模態分量進行檢測時選取排列熵為0.6。

表1 幾種信號的排列熵

1.3 基于奇異值分解的循環策略

在噪聲信號先驗知識未知的情況下，對于固有模態分量中信號的有效提取是個難題。根據論文1.2節排列熵分析，本文提出了基于奇異值分解的循環策略對有用信號進行提取，具體步驟如下：

(1)對含噪信號進行CEEMD處理，得到一系列本征模態分量ci(t)和余項rn(t)。

(2)對ci(t)進行排列熵檢測，若PEi>0.6，則認為第i個IMF為噪聲主導，若PEi≤0.6，則認為第i個IMF為信號主導。

(3)對噪聲主導IMFs進行奇異值分解，由奇異值分解基本原理可知[13]：信號經奇異值分解，得到1-d個奇異值分量，前面較大的k個奇異值由純信號和噪聲共同貢獻，其對應的奇異矢量所張成的空間即為加噪信號子空間，后面相對較小d-k個奇異值完全由噪聲貢獻，將其剔除，可以減少含噪信號中的噪聲，奇異值分解降噪關鍵在于k值得選取，可以做個假設：若重構前k個奇異值分量得到的信號排列熵不大于0.6，則認為奇異值分解有效重構階次為k，由EMD的原理可知，第一IMF包含了絕大部分噪聲，會出現只重構第一個奇異值分量仍不能滿足其排列熵檢測的情況[14]，因此需要設計一個循環策略，提取第一個奇異值分量繼續進行奇異值分解，直至滿足排列熵檢測，循環過程如圖1所示。

圖1中PE[SVDk]表示前k個奇異值分量的排列熵，max(k)表示滿足PE[SVDk]≤0.6時k的最大值，目的是為了盡可能減少降噪過度。

(4)經步驟(3)處理后，可得到一系列信號主導分量，最后經過平滑濾波，去除信號中的毛刺、尖銳部分，得到降噪后的信號，考慮到Savitzky-Golay濾波是一種滑動窗口中心點的有效擬合方法[15]，本文選用SG濾波方法對信號進行平滑處理。

2 仿真和應用分析

2.1 仿真驗證分析

對Lorenz方程

(6)

取參數σ=10,r=28,b=8/3時，是一種混沌信號[16]，積分步長為0.01，采用四階龍格-庫塔算法運算，取1 000個樣本點，干凈信號和含噪信號如圖2所示。

圖3為Lorenz含噪信號的CEEMD結果，分解為8個本征模態分量(imf1～imf8)和一個余項ri(t)，高頻部分沒有明顯的信號混疊特征，說明了經過白噪聲多次疊加和平均，克制了信號中的偽極值點，減少了模態混疊。圖4分別為傳統CEEMD方法、CEEMD-小波閾值方法和本文方法的降噪結果，傳統CEEMD方法是直接將噪聲主導的本征模態分量剔除，重構信號主導的本征模態分量，CEEMD-小波閾值方法是利用小波閾值處理提取噪聲主導的本征模態分量中的有用細節，然后和號主導的本征模態分量疊加達到降噪目的，圖5為3種方法降噪后的重構誤差，即是降噪后信號和干凈Lorenz信號的差值，三種方法的重構誤差分別為0.147、0.051 6、0.040 2，通過比較可以看出，傳統CEEMD方法的重構誤差明顯高于其他2種方法，其在去除高頻噪聲的同時，也會使其中的信號丟失，增加了信號失真，所以其重構誤差中明顯包含了有很多有用信號特征；CEEMD-小波閾值方法將閾值以上的信號保留，會使部分噪聲遺留，尤其是對于噪聲主導的信號，效果很不理想，圖4b中信號明顯留有粗糙細節，也印證了這一論述；本文方法是一種多次篩選提取信號的方法，即使是低信噪比的信號，經過多次循環篩選，也可以有效將隱藏其中的信號提取出來，并且相比上述兩種方法，本文方法降噪后信號更加光滑，重構誤差最小，證實了所提方法的優越性。

為了進一步驗證所提方法的優越性，本文通過信噪比(SNR)、均方誤差(MSE)分析上述兩個參數的選取，2個參數的計算公式如下：

(7)

(8)

其中：x′(n)表示去噪后的序列，x(n)為原始時間序列，var(·)表示方差，x′(n)-x(n)表示信號中的剩余噪聲。信噪比反應去噪能力的大小，均方誤差的物理意義是表示去噪后信號和原始信號的平均偏離程度，通常認為，SNR越大，MSE越小，說明去噪效果越好。

表2給出了EEMD、傳統CEEMD方法、CEEMD-小波閾值方法和上述3種方法降噪的各項指標。可以看出，本文方法降噪后信號信噪比最高，均方根誤差最小，體現了所提方法對于低信噪比信號提取的優越性；同時，表2是同步比較3種方法的計算時間，CEEMD是同時疊加相反的白噪聲，所以疊加噪聲總次數和EEMD一致，可以看出CEEMD時效性較EEMD有很大的提升，本文方法由于增加了基于奇異值分解的多次循環提取信號，所以計算耗時有所提升，但是所提循環策略是一種自適應的手段，且奇異值分解運算簡單，因而總耗時增加不多，但降噪后信號信噪比卻比CEEMD-小波閾值降噪方法提升了近1dB。

表2 4種方法各項指標比較

2.2 應用分析

對CEEMD-小波閾值和本文方法降噪后的含噪Lorenz信號利用Volterra預測模型進行預測，預測絕對誤差分別為Perr1=0.002 4、Perr2=2.874 1×10-4，可以看出同一預測模型下本文方法將預測精度提升了一個數量級，且圖6a中的誤差曲線有幾個峰值誤差接近0.1，這種大誤差點很可能導致混沌預測在應用中出現錯誤結論，而6b中預測誤差曲線沒有特別明顯的誤差峰值。

混沌振動特征分析是混沌工程實際應用的關鍵環節，基于實驗室條件下的混沌運動特性分析也是將混沌應用于工程實際的重要前提，為此，以文獻[17]產生的混沌信號為對象，驗證所提方法對實測混沌信號的適用性，試驗實物圖如圖7所示，當激勵增益為1V，激勵頻率為9.8 Hz時，實測時間序列如圖8a所示，嵌入維數為4，延遲時間為16，得到混沌吸引子相圖如圖8c所示，利用所提方法降噪前后實測信號的功率譜和預測結果8e、f和圖9所示，可以看出，降噪后預測絕對誤差大大減小，證實了本文方法對無論仿真信號還是實測混沌信號都是有效的，同時證實了本文方法可為混沌運動分析提供了有效的預處理手段。

3 結語

提出了基于CEEMD-排列熵的循環策略信號提取方法。CEEMD較于EEMD，在抑制白噪聲引起重構誤差的同時簡化了計算方法，對分解得到的本征模態分量通過計算排列熵確定噪聲分量和信號分量，在奇異值分解的基礎上，建立循環策略的信號提取方法，該方法無需信號和噪聲的先驗知識，在抑制噪聲的同時，可以較好地保留有用信號。通過仿真信號和實測混沌信號降噪處理，結果表明所提方法較EEMD時效性有很大的提升，降噪效果較CEEMD-小波閾值提升了近1 dB，且本文方法對較低信噪比信號的有更好的降噪效果，并可為混沌實測信號特征分析提供有效的預處理手段。