尹昱東 明 勇 邊 羽

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

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

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

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

1 算法基本原理

1.1 CEEMD算法

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

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

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

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

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

(1)

(2)

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

1.2 排列熵

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

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

(3)

其中:m為嵌入維數(shù)，為時延。

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

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

(4)

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

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

(5)

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

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

表1 幾種信號的排列熵

1.3 基于奇異值分解的循環(huán)策略

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

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

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

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

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

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

2 仿真和應(yīng)用分析

2.1 仿真驗證分析

對Lorenz方程

(6)

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

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

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

(7)

(8)

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

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

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

2.2 應(yīng)用分析

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

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

3 結(jié)語

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