尹昱東 明 勇 邊 羽

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

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

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

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

1 算法基本原理

1.1 CEEMD算法

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

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

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

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

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

(1)

(2)

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

1.2 排列熵

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

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

(3)

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

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

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

(4)

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

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

(5)

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

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

表1 幾種信號(hào)的排列熵

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

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

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

(2)對(duì)ci(t)進(jìn)行排列熵檢測(cè)，若PEi>0.6，則認(rèn)為第i個(gè)IMF為噪聲主導(dǎo)，若PEi≤0.6，則認(rèn)為第i個(gè)IMF為信號(hào)主導(dǎo)。

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

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

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

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

2.1 仿真驗(yàn)證分析

對(duì)Lorenz方程

(6)

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

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

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

(7)

(8)

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

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

表2 4種方法各項(xiàng)指標(biāo)比較

2.2 應(yīng)用分析

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

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

3 結(jié)語(yǔ)

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