低采樣率下經驗模態分解性能提升研究

黎　恒，李　智，莫　瑋

(1.西安電子科技大學 機電工程學院,西安　710071; 2.桂林航天工業學院 自動化系,桂林　541004; 3.桂林電子科技大學 電子工程與自動化學院,桂林　541004)

?

低采樣率下經驗模態分解性能提升研究

黎恒1，李智2,3，莫瑋1

(1.西安電子科技大學 機電工程學院,西安710071; 2.桂林航天工業學院 自動化系,桂林541004; 3.桂林電子科技大學 電子工程與自動化學院,桂林541004)

經驗模態分解(EMD)使用信號極值點的位置和取值信息進行分解，對采樣率有較高的要求。針對EMD在低采樣率下性能降低的現象，提出一種基于B樣條擬合的信號局部均值計算方法。首先提取信號極值點出現的時刻作為尺度，然后通過對極值點時刻進行重新采樣構造B樣條節點，最后應用B樣條最小二乘擬合方法直接計算局部均值。與EMD方法相比，該方法不需要信號極值點的準確位置和取值，因此不容易受到低采樣率的影響。對平穩信號和非平穩信號的仿真結果表明，該方法能在接近奈奎斯特頻率的低采樣率下獲得較高的性能。與基于插值的解決方案相比，該方法的分離性能更好。

經驗模態分解；B樣條擬合；低采樣率；信號分解；時頻分析

希爾伯特-黃(HHT)變換[1]作為一種時間-頻率分析方法，已經成為一種分析非線性和非平穩信號的標準方法。HHT的核心是經驗模態分解(EMD)。通過EMD將輸入信號分解成為一系列單模態分量，這些分量被稱為本征模態函數(IMF)。在HHT中，對IMF進行希爾伯特變換即得到對應的時間-頻率-幅度分布。與小波變換等方法不同，EMD方法不需要定義任何基函數，具有很強的自適應性。因此EMD在振動分析、機械系統、故障檢測等領域得到了廣泛的應用[2-7]。盡管如此，EMD方法也存在若干需要解決的問題，如低采樣率誤差[8]、模態混疊[9]和邊緣效應[10]等。EMD方法的基本原理是通過一種篩分算法，不斷地從輸入信號中去除均值分量，從而得到關于時間軸對稱的IMF。在篩分過程中，算法使用信號的極大值和極小值點產生均值曲線。因此，信號極值點的計算十分重要。在實際應用中，信號采樣率會受到設備和實驗條件的限制，往往不可能遠大于信號的奈奎斯特頻率。當采樣率較低時，極值點的位置和取值將會出現誤差，這些誤差直接影響著均值曲線的精度，進而影響分解結果的精度。

針對這個問題，現有改進方法主要專注于對信號進行插值，然后對極值點進行定位。文獻[11-12]給出基于sinc插值的重構方法，對信號進行插值，再進行EMD分解。文獻[13]對信號包絡使用升余弦插值，提高信號的等效采樣率。這些方法在一定程度上提升了EMD在低采樣率的性能。但是，sinc插值和升余弦插值做插值計算時收斂慢，需要較多的采樣點，并且存在嚴重的邊界失真；另一方面，這些方法對線性信號有很好的效果，而對非平穩信號可能產生額外誤差；再者，對間斷信號插值時通常會產生震蕩，表現為較大的上沖或下沖，為接下來的EMD分解帶來困難；最后，在采樣率接近信號的奈奎斯特頻率時，往往存在較大的誤差。文獻[14]使用二代小波插值對信號極值點進行定位，提高了極值點預測的精度。這些基于插值的方法對極值點位置和取值信息恢復的精度依賴性較強，并且其分解性能往往受限于EMD本身的精度。

由于目前基于插值的EMD改進方法存在上述弊端，本文提出一種基于B樣條擬合的均值計算方法。該方法根據輸入信號的極值點分布情況進行節點選取并通過擬合得到局部均值。由于不需要對信號插值進行極值點定位，因此該方法不存在極值點定位不準確帶來的誤差。本工作首先研究采樣率對EMD方法的影響，然后針對這一問題提出基于B樣條擬合的改進算法，并給出實現的原理和步驟，最后通過仿真實驗來驗證該方法的有效性和優越性。

1　采樣率對EMD分解結果的影響

EMD方法以信號的極值點為特征，將輸入信號分解成一組包絡關于時間軸對稱的IMF分量[1]。EMD計算步驟中的均值信號通過求其上下包絡曲線的均值得到。由于上下包絡曲線是通過對信號的極值點進行插值得到的，極值點的準確性直接影響著分解性能。隨著采樣率的降低，極值點的位置和取值將會偏離真實情況。文獻[8]深入研究了采樣率對EMD分解性能的影響并給出相關理論推導，這里不再贅述。為直觀理解，給出實驗進行說明。考察由10 Hz和3 Hz頻率組成的雙音信號s(t)=cos(2π×10t)+cos(2π×3t)。分別使用高采樣率f1=1 000 sps和低采樣率f2=25 sps對其進行采樣得到離散時間序列，然后進行EMD分解。分解結果中第一階IMF對應著頻率較高的信號分量，即cos(2π×10t)。該分量的比較結果如圖1(a)和圖1(b)所示，可以明顯看出采樣率對分解結果有明顯的影響。圖1(c)為低采樣率下高頻信號分量，對應著預期準確的分解結果。通過對比可以直觀地看出圖1(b)和圖1(c)有明顯的不同，說明在低采樣率下EMD分解帶來了較大的誤差。造成該現象的原因是在低采樣率下，極值點的出現時刻和取值出現了波動，導致上/下包絡出現了變化。為進一步說明低采樣率帶來的分解誤差，考慮輸入信號x(t)=cos(2π×10t)+cos(2π×10ft)。f為頻率比，當f<1時cos(2πft)對應著低頻信號分量。在不同采樣率下對x(t)進行EMD分解，并計算分解結果第一階IMF與真實結果的均方誤差。實驗中f取0.1～0.4，即對頻率比不同的四組雙音信號進行分解。采樣率范圍設置為25～1 000 sps。圖2為均方誤差對采樣率的變化曲線。從圖中可以得到以下結論。首先在低采樣率下，特別是采樣率小于100 sps時(或采樣率/信號頻率<10)，分解結果有較大誤差。其次，不同頻率比的雙音信號具有相似的誤差曲線。最后，當采樣率/信號頻率>10時，分解結果的誤差趨于平緩。通過該實驗可以得到一個經驗法則，即采樣率應至少取有用信號頻率的10倍才能保證分解精度。

下面從另一個角度說明EMD在低采樣率下導致的問題，即在低采樣率下IMF定義存在的問題。由于希爾伯特變換只能處理窄帶信號，需要先對輸入信號進行分解。EMD分解即是將輸入信號分解成一系列局部窄帶分量的方法。HUANG[1]認為單模態信號應與窄帶信號的特征類似，其包絡是關于時間軸對稱的。然而他忽略了一個重要的事實，在低采樣率下單模態信號的包絡不一定關于時間軸對稱(圖1(c))。EMD方法強制要求IMF包絡的對稱性，在低采樣率下將不可避免地產生誤差。因此在這種情況下，使用信號上下包絡產生均值是不恰當的。

圖1　不同采樣率下EMD分解結果Fig.1　Decomposition under different sampling rates

圖2　不同采樣率下分解第一階結果與真實值的均方誤差Fig.2　Mean squared error of the first IMF under various sampling rates

2　基于B樣條擬合的均值計算方法

B樣條是信號分析的重要工具，其計算方法包括B樣條插值、B樣條擬合和平滑B樣條等。一些學者使用B樣條對EMD進行了改進。文獻[15]使用極值點的值做二項平均作為B樣條控制系數進行均值計算，文獻[16]使用平滑B樣條生成上/下包絡。與傳統EMD方法類似，這些方法都使用了極值點的位置和取值，在低采樣率下性能會受到影響。與已有方法不同，本文方法利用B樣條擬合的低通濾波器性質，能很好地抑制低采樣率導致的分解性能降低。

2.1B樣條擬合

定義域[a,b]內的p階B樣條擬合的曲線g(t)表示為：

(1)

式中：cj為待求解的控制系數，Bj,p為p階B樣條函數。對[a,b]內給定的非遞減節點u={u0,u1,u2,…,un+p}，Bj,p使用Cox-de Boor遞歸公式計算[16]：

(2)

(3)

當p>1時，式(3)的分子和分母中都存在節點，因此Bj,p(t)乃至g(t)是關于節點的非線性函數。對輸入信號y(t)進行B樣條擬合可以表示為：

y(t)=g(t)+ε(t)

(4)

式中：ε(t)為擬合誤差。B樣條擬合就是求解關于控制系數cj的最小二乘問題。即求出cj，使Q最小

(5)

B樣條是非線性函數，同時也是擬合非線性和非平穩信號的重要工具[16]。研究表明[17],B樣條擬合具有低通濾波器的性質。濾波器的局部截止頻率約為1/2m，其中fs為信號的采樣頻率，m為節點的局部距離。以等間隔節點為例，當節點間隔為2時，B樣條擬合相當于截至頻率為1/4的濾波器；節點間隔為1.5時，B樣條擬合相當于截至頻率為1/3的濾波器。該研究還指出，B樣條階數p趨向于正無窮時，該濾波器趨近于理想低通濾波器。

2.2算法實現

EMD方法認為，極值點的分布包含信號局部尺度信息。研究表明，當信號中高頻模態幅度較明顯時，極值點的分布往往和該高頻模態分量的尺度一致[17]。通過前面分析可知，B樣條擬合有濾波器的性質，其截至頻率可控。在信號分解中，均值信號通常是頻率較低的分量。因此，將極值點分布特征和B樣條擬合相結合，可以計算出均值信號。本文方法的實質是使用B樣條擬合濾波器實現高頻分量和低頻分量的分離。其基本步驟是：使用信號極值點出現時刻作為初始特征，然后對這些時刻進行重新采樣得到時間分布密度較低的節點，最后使用該節點進行B樣條擬合得到均值。本文關于均值信號的計算方法按照如下步驟進行：

(1)找出y(t)所有極值點出現的時刻，記為vi，i=1,2,3，…；

(2)在每個vi和vi+1之間以(vi+1-vi)/10為步長等間隔插入新的時刻，得到的序列記為ri，i=1,2,3,…；

(3)以k>10(本文取14)作為間隔對ri進行采樣，得到的序列記為si，i=1,2,3，…；

(4)使用si作為節點進行p階(本文取p=14)B樣條擬合，并將擬合曲線作為均值曲線m(t)；

計算出均值曲線后，剩下的處理步驟和EMD方法一致。步驟(3)中，k的取值決定了ri較vi的“稀疏”程度。以極值點等間隔的單音信號為例，設vi的時間間隔為τ，通過計算得到ri的時間間隔為τk/10。在這種情況下，以ri作為節點的B樣條濾波器截至頻率為10/2kτ。注意到該單音信號的頻率為1/2τ>10/2kτ，因此當k>10時，B樣條濾波器可以濾除該單音信號，得到頻率較低的均值分量。在低采樣率下，極值點出現時刻會出現微小波動，但是這些微小波動并不會引起si間隔的顯著變化，因此不會引起B樣條擬合濾波器截至頻率的顯著變化。

3　仿真實驗

由于實際信號往往存在噪音、失真等現象，會對分解結果帶來額外的影響，不利于性能評價。本文使用仿真信號來進行分析算法在低采樣率下的性能，包括平穩信號和非平穩信號。

3.1平穩信號的分離性能

本小節研究本文方法對平穩信號分解性能。考慮下面的雙音信號模型：

x(t;a,f)=cos2πt+acos(2πft+φ)

(6)

式中：a為幅度比，f為頻率比，φ為相位差。當f∈(0,1)時，式(6)中第一項為高頻分量，第二項為低頻分量。由EMD方法的分解原理可知，先分解出來的IMF為頻率較高的分量。為評價分解性能，考察第一階IMF與高頻分量的相似度，使用文獻中常用的公式作為評判標準[18]：

(7)

圖3　不同采樣率下分離性能比較Fig.3　Comparison of separation quantity between EMD and the proposed method under different sampling rates

當分解出的IMF與高頻分量完全相同時，式(7)為零。定義一個閾值，當c(a,f,φ)<0.5時認為分解是成功的。相位差φ通常被認為對分解結果影響不大，本文隨機取5個相位分解結果的平均作為輸出。圖3給出了EMD和本文方法在不同采樣率fs下的分解結果。圖中黑色部分為分解成功的區域，白色部分為分解失敗的區域。圖3(a)和圖3(b)給出EMD方法在采樣率fs=5fN和fs=1.5fN下的分解結果，其中fN為高頻分量的奈奎斯特頻率。從分解結果看，采樣頻率對EMD分離性能有較大影響。當采樣率降低時，代表分解成功的黑色區域面積明顯減小。圖3(c)和圖3(d)給出了本文方法fs=1.5fN和fs=1.1fN下的分解結果。與EMD方法相比，本文方法成功分解的區域面積更大。隨著采樣率的進一步降低，當fs=1.1fN接近奈奎斯特頻率時，本文方法的分解性能并未出現明顯下降。這些實驗結果充分驗證了本文方法性能在低采樣率下的穩定性和優越性。

值得注意的是，本文方法中，k的取值對頻率分離性能有直接影響。回顧本文算法，k=14的情況下，B樣條擬合濾波器的截至頻率約為0.35τ=0.7×0.5τ。因此該算法能分離頻率比f=0.7的信號分量。同理可以得到，當k=12時，f=0.8，當k=11時，f=0.9。圖3(e)給出了分解結果中分解成功/失敗過渡區域的輪廓曲線，即對應C(a,f)=0.5的曲線。使用三組不同的配置以考察頻率分辨率并在結果中使用不同的線條表示：k=11用虛線，k=12用實線，k=14用短劃線。圖中還給出EMD分解結果，用點虛線表示。可以看出，參數k直接決定著SAMD的分解的頻率分辨率，其性能與理論預期一致。

3.2非平穩信號的分離性能

使用FM信號作為非平穩信號的例子來評價本文方法的分解性能。

y(t)=cos(2πft+0.12πft×cos2πt)+

acos(πft+0.06πft×sin2πt)

(8)

式中：t∈[0,1]，FM信號的中心頻率為f=1.5fN。對該信號分別使用EMD、基于sinc插值EMD和本文方法進行分解。為直觀地觀察分解性能，對分解結果求希爾伯特譜。EMD分解結果如圖4(a)所示，可以很明顯地看到低采樣率對EMD分解的影響。希爾伯特譜出現了嚴重的失真，甚至在低頻段出現了多余的雜散分量。圖4(b)給出了基于sinc插值重構的分解結果。從分解結果看，插值確實帶來了性能的提升，但仍存在明顯的失真和邊界問題。本文方法分解結果如圖4(c)所示。可以看出，高頻和低頻分量的時頻分布清晰可辨，失真低，說明本文方法更好地實現了信號的分離。

圖4　分解結果的希爾伯特譜Fig.4　Comparison of Hilbert spectra

3.3精度

選取雙音信號的分解精度進行考察，比較EMD、二代小波插值[14]和本文方法的性能。本小節實驗的配置與第一章中演示實驗類似，考慮輸入信號x(t)=cos(2π×10t)+cos(2π×10ft)。在不同采樣率下對x(t)進行EMD分解，并計算分解結果第一階IMF與真實結果的均方誤差。實驗中f取0.1到0.4，采樣率范圍設置為25～1 000 sps。對不同雙音信號分解得到的結果進行算術平均并取對數，實驗結果如圖5所示。首先比較采樣率為25 sps到100 sps的分解性能。與EMD相比，基于二代小波插值的方法在該區域的分解精度得到明顯提升，其中當采樣率為25 sps時精度提升了1個數量級。本文方法獲得的精度提升更為明顯，當采樣率為25 sps時精度提升了約1.5個數量級。下面比較100 sps到1 000 sps的分解精度。在該區域插值已無法帶來精度的提升，這是因為插值對極值點位置預測的改善已經很微小，另一方面插值方法帶來的精度提升受限于EMD本身的分解精度。為比較B樣條階數p對分解精度的影響，本文方法分別使用p=14和p=10兩種階數進行分解對比。可以看到本文方法在整個采樣率范圍都可以帶來分解性能的提升。另一方面，提高B樣條階數p能提升分解精度。這是因為p取值越大B樣條濾波器的滾降越大，越接近理想低通濾波器。然而當p的取值過大時，會使式(5)的求解變得更復雜。因此增大p帶來的性能提高是以更多計算時間為代價的。如本文取p=14時，能帶來較明顯的精度改善，計算時間約為EMD方法的4倍。

圖5　不同采樣率下分解精度比較Fig.5　Comparison of accuracy under various sampling rates

4　結　論

本文針對低采樣率下EMD存在的分解不準確的問題，提出一種基于B樣條擬合的均值計算方法。仿真實驗表明通過本文方法能有效地消除低采樣率帶來的影響。本文方法與現有方法的不同之處在于:

(1)該方法的控制系數通過最小二乘方式確定，直接計算均值信號，不需要計算上/包絡。

(2)不需要知道極值點的準確位置或準確取值。該方法只需要知道信號極值點大致出現的時刻。因此，該方法對低采樣率有很好的適應性。但是該方法目前尚不能解決模態混疊問題，在今后的工作中，如何抑制模態混疊是重點研究內容。

[1]HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society of London.Series A,1998,454(1971):903-995.

[2]吳響,錢建生,王海燕.微震信號多尺度非線性特征提取與辨識研究[J].儀器儀表學報,2014,35(5):969-975.

WU Xiang,QIAN Jiansheng,WANG Haiyan.Study on multi-scale nonlinear feature extraction and signal identification for microseismic signal[J].Chinese Journal of Scientific Instrument,2014,35(5):969-975.

[3]王錄雁,王強,張梅軍,等.基于 EMD 的滾動軸承故障灰色診斷方法[J].振動與沖擊,2014,33(3):197-202.

WANG Luyan,WANG Qiang,ZHANG Meijun,et al.A grey fault diagnosis method for rolling bearing based on EMD [J].Journal of Vibration and Shock,2014,33(3):197-202.

[4]YANG Z,LING B W K,BINGHAM C.Trend extraction based on separations of consecutive empirical mode decomposition components in Hilbert marginal spectrum[J].Measurement,2013,46(8):2481-2491.

[5]KOMATY A,BOUDRAA A O,AUGIER B,et al.EMD-based filtering using similarity measure between probability density functions of IMFs[J].IEEE Transactions on Instrumentation and Measurement,2014,63(1):27-34.

[6]GOHARRIZI A Y,SEPEHRI N.Internal leakage detection in hydraulic actuators using empirical mode decomposition and hilbert spectrum[J].IEEE Transactions on Instrumentation and Measurement,2012,61(2):368-378.

[7]LEI Y,LIN J,HE Z,et al.A review on empirical mode decomposition in fault diagnosis of rotating machinery[J].Mechanical Systems and Signal Processing,2013,35(1):108-126.

[8]RILLING G,FLANDRIN P.On the influence of sampling on the empirical mode decomposition[C]// IEEE International Conference on Acoustics,Speech and Signal Processing.Toulouse:IEEE,2006:444-447.

[9]肖瑛,殷福亮.解相關 EMD:消除模態混疊的新方法[J].振動與沖擊,2015,34(4):25-29.

XIAO Ying,YIN Fuliang.Decorrelation emd:a new method of eliminating mode mixing[J].Journal of Vibration and Shock,2015,34(4):25-29.

[10]王錄雁,王強,魯冬林,等.EMD自適應三角波匹配延拓算法[J].振動與沖擊,2014,33(4):94-99.

WANG Luyan,WANG Qiang,LU Donglin,et al.A self-adaptive triangular waveform matching extension algorithm for EMD[J].Journal of Vibration and Shock,2014,33(4):94-99.

[11]XU Z,HUANG B,ZHANG F.Improvement of empirical mode decomposition under low sampling rate[J].Signal Processing,2009,89(11):2296-2303.

[12]ROY A,DOHERTY J F.Improved signal analysis performance at low sampling rates using raised cosine empirical mode decomposition[J].Electronics Letters,2010,46(2):176-177.

[13]CHEN Q,HUANG N,RIEMENSCHNEIDER S,et al.A B-spline approach for empirical mode decompositions[J].Advances in Computational Mathematics,2006,24(1/2/3/4):171-195.

[14]趙利強 王建林 于濤.基于改進EMD的輸油管道泄漏信號 特征提取方法研究[J].儀器儀表學報,2013,34(12):2696-2702.

ZHAO Liqiang,WANG Jianlin,YU Tao.Study on the extraction method for oil pipeline leakage signal feature based on improved empirical mode decompositio [J].Chinese Journal of Scientific Instrument,2013,34(12):2696-2702.

[15]KIM D,PARK M,OH H S.Bidimensional statistical empirical mode decomposition[J].Signal Processing Letters,IEEE,2012,19(4):191-194.

[16]UNSER M,ALDROUBI A,EDEN M.B-spline signal processing.I.Theory[J].IEEE Transactions on Signal Processing,1993,41(2):821-833.

[17]UNSER M.Splines:A perfect fit for signal and image processing[J].IEEE Signal Processing Magazine,1999,16(6):22-38.

[18]RILLING G,FLANDRIN P.One or two frequencies? The empirical mode decomposition answers[J].IEEE Transactions on Signal Processing,2008,56(1):85-95.

Performance improvement of EMD under low sampling rates

LI Heng1,LI Zhi2,3,MO Wei1

(1.School of Mechano-Electronic Engineering,Xidian University,Xi’an 710071,China;2.Guilin University of Aerospace Technology,Guilin 541004,China;3.School of Electronic Engineering and Automation,Guilin University of Electronic Technology,Guilin 541004,China)

Empirical mode decomposition (EMD)depends highly on exact locations and values of signal extrema,they need a higher sampling rate.Aiming at improving the performance of EMD under low sampling rates,a local mean estimation method based on B-spline fitting was proposed.Firstly,the occuring instant of extrema was extracted as a time scale.Then,the occuring instant of extrema was re-sampled to construct node points of B-spline.Finally,the local mean was computed directly with B-spline least squares method.Compared with EMD method,the proposed method did not need the exact locations and values of extrema,so the low sampling rates were not easy to affect this technique.The efficiency of this technique was demonstrated using simulated signals.The simulation results showed that the performance of the proposed method is not reduced even the sampling rate is close to Nyquist rate; the proposed method is superior to existing interpolation methods in separation performance.

empirical mode decomposition (EMD); B-spline fitting; low sampling rates; signal decomposition; time-frequency analysis

國家自然科學基金(61361006)

2015-04-07修改稿收到日期：2015-09-16

黎恒 男，博士生，1982年4月生

李智 男，博士，教授，1965年5月生

TH137

A DOI:10.13465/j.cnki.jvs.2016.17.031