PDE信號修補(bǔ)方法及在EMD端點(diǎn)效應(yīng)處理中的應(yīng)用

尹愛軍，王 璇

（重慶大學(xué) 機(jī)械傳動國家重點(diǎn)實驗室，重慶大學(xué) 機(jī)械學(xué)院測試中心，重慶 400044）

PDE是一個含有多元未知數(shù)函數(shù)及其偏導(dǎo)數(shù)的方程。PDE是在討論自然現(xiàn)象（特別是物理現(xiàn)象）的過程中逐步建立起來的，所以也稱為數(shù)學(xué)物理方程。PDE最新的研究與應(yīng)用已擴(kuò)展到經(jīng)濟(jì)、金融預(yù)測、圖像處理、振動信號處理等領(lǐng)域［1－4］。

二階線性PDE是偏微分方程理論中較為成熟的部分，并成為其他分支借鑒的典范。熱傳導(dǎo)方程見式（1）是最簡單、最基礎(chǔ)的拋物型PDE，可用于描述熱傳導(dǎo)、擴(kuò)散等物理現(xiàn)象，得到了非常廣泛的應(yīng)用。

熱傳導(dǎo)方程初值問題的解實際上就是一個高斯濾波器，即當(dāng)f（x，t）=0時，利用傅立葉積分法可得式（1）的解為：

其中，

因此，熱傳導(dǎo)方程在信號濾波降噪中得到了廣泛應(yīng)用，特別是在二維圖像信號的濾波降噪中。

Hilbert－Huang變換是由 Huang等［5］提出的一種新的非平穩(wěn)信號分析方法。HHT包含兩個過程：首先對信號進(jìn)行非線性的自適應(yīng)分解，稱為經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition，EMD），它把復(fù)雜信號分解成有限個本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF）之和；然后對各個IMF進(jìn)行Hilbert變換，研究信號的時頻能量分布。

如果定義y（t）為原信號序列，imf（t）為經(jīng)EMD得到的本征模態(tài)函數(shù)（IMF），r（t）為余量，則原始信號可以表示為所有的IMF及余量之和：

對于每一個IMF，應(yīng)滿足：① 在整個數(shù)據(jù)段內(nèi)，極值點(diǎn)的個數(shù)和零交叉點(diǎn)的個數(shù)必須相等或相差最多不能超過一個。② 在任何一點(diǎn)，由局部極大值點(diǎn)形成的包絡(luò)線和由局部極小值點(diǎn)形成的包絡(luò)線的平均值為零。

HHT分析的核心是EMD分解。EMD分解過程中，每一個IMF都需要多次篩選過程來實現(xiàn)，而每一次篩選過程必須找到局部極大值構(gòu)成的上包絡(luò)和局部極小值構(gòu)成的下包絡(luò)，并用三次樣條插值來分別計算出上下包絡(luò)的插值及其平均值。由于信號的端點(diǎn)一般不會同時是局部極值點(diǎn)，因此利用EMD進(jìn)行分解時，因包絡(luò)線不準(zhǔn)確會引起誤差，并且這種誤差會隨著IMF分解層數(shù)的增加而向內(nèi)傳播，繼而“污染”整個數(shù)據(jù)序列，使得最后的分解結(jié)果失去意義。為此，眾多研究人員相繼提出了相關(guān)的端點(diǎn)處理方法［6］，如端點(diǎn)鏡像外延［7］、多項式外延［8］、神經(jīng)網(wǎng)絡(luò)處理［9］、支持向量機(jī)預(yù)測［10］、相似極值延拓［11］等方法。另有研究人員則根據(jù)EMD方法的基本思想，提出了局部均值分解改進(jìn)方法［12］。

本文針對振動信號的內(nèi)在性質(zhì)，分析基于插補(bǔ)等方法的信號修補(bǔ)方法，提出了PDE信號修補(bǔ)新方法，并將這一方法成功用于EMD分解過程中的端點(diǎn)效應(yīng)處理。

1 基于PDE的波動信號修補(bǔ)

1.1 信號修補(bǔ)方法

對于數(shù)字信號分析系統(tǒng)來講，在將模擬信號轉(zhuǎn)換為數(shù)字信號的過程中，總會存在各種各樣的干擾，或者硬件設(shè)備本身的性能指標(biāo)等因素，使得數(shù)字信號可能會存在失真。如采樣頻率太高，而分析系統(tǒng)處理速度過慢，從而引起的分析信號不連續(xù)、數(shù)據(jù)掉點(diǎn)等等。如圖1所示，圖1（a）是理想信號，圖1（b）表示了信號不連續(xù)情況。因此，需要對這一類信號進(jìn)行修補(bǔ)，復(fù)原信號。

對于一個一般實際系統(tǒng)，總可以將其表示為如式（5）所示的系統(tǒng)形式。實際系統(tǒng)內(nèi)部一般存在儲能元件和耗能元件，因此系統(tǒng)的輸出具有連續(xù)性，或者局部光滑性。利用這一特點(diǎn)，我們可以設(shè)計多種修補(bǔ)方式。

圖1 原始信號與缺陷信號Fig.1 Original signal and distortion signal

拉格朗日插值、牛頓插值、埃爾米特插值、樣條插值等方法即是利用實際信號的連續(xù)性和光滑性所實現(xiàn)的典型信號修補(bǔ)方法。式（6）為拉格朗日多項式插值表達(dá)式。顯然，若數(shù)據(jù)量小，缺陷信號的時間很短，則可以直接利用中值插值、線性插值等方法。

1.2 PDE信號修補(bǔ)原理

熱傳導(dǎo)方程表示了一種傳熱過程（或者擴(kuò)散過程），且這種傳導(dǎo)過程是連續(xù)光滑的。并且我們注意到，熱傳導(dǎo)方程初值問題的解就是一個高斯濾波器。因此，對于上述畸變信號同樣可以用PDE方法進(jìn)行修補(bǔ)。

為簡化計算過程，令式（1）中f（x，t）=0，a=1，即有：

設(shè)原始振動信號為y0（x），畸變信號位置范圍為［x1，x2］，于是可以得到如下修補(bǔ)公式：

即對于畸變?nèi)毕菪盘?，根?jù)熱傳導(dǎo)過程中的連續(xù)及光滑性，逐漸重構(gòu)出原始信號；而對于正常信號，則不作處理。

式（8）是直接根據(jù)高斯濾波實現(xiàn)的線性PDE修補(bǔ)方法，這一方法沒有考慮信號之間的變化，從而可能會影響修補(bǔ)信號的瞬態(tài)特性。為此，可以采用式（8）的非線性 PDE 修補(bǔ)方式［1，13－14］。

其中k（·）表示一個特定函數(shù)，如：

其中p、K為調(diào)制系數(shù)。grad（·）為梯度函數(shù)，典型的有前向差分［式（11a）］、后向差分［式（11b）］、中心差分［式（11c）］等。

在這一修補(bǔ)方式下，可以針對信號梯度設(shè)置相應(yīng)的調(diào)制參數(shù)，使修補(bǔ)后的信號即滿足整體連續(xù)光滑、又滿足瞬態(tài)沖擊特性。

若定義如下模板函數(shù)：

x1，x2表示修補(bǔ)范圍。則式（9）可以表示為：

1.3 仿真試驗

設(shè)由兩個正弦信號組成的理想信號，如式（14）。設(shè)置采樣頻率為1 000 Hz，取連續(xù)1 024個數(shù)據(jù)點(diǎn)（y1，y2，…，y1024），如圖 1（a）所示?；冃盘枮閥50，y51，…，y100，如圖 1（b）所示。

圖2（a）為拉格朗日多項式插值修補(bǔ)效果；圖 2（b）為三次樣條插值修補(bǔ)效果；圖2（c）為采用式（8）所示的線性PDE修補(bǔ)效果；圖 2（d）為采用式（9）所示的非線性PDE修補(bǔ)效果。其中式（8）、式（9）的迭代次數(shù)均為1 000次。對于非線性 PDE修補(bǔ)，選擇式［10（a）］作為調(diào)制函數(shù)，并令p=－0.01，梯度函數(shù)為式（11a）。從圖可以看出，上述方法均可實現(xiàn)信號的修補(bǔ)，非線性PDE因考慮了信號內(nèi)部梯度，修補(bǔ)效果比線性PDE好。

圖2 信號修補(bǔ)Fig.2 Signal restored

2 EMD端點(diǎn)效應(yīng)的PDE處理

2.1 PDE端點(diǎn)效應(yīng)處理

實際上，上述所提到的這些方法，都是利用信號的連續(xù)性（樣條修補(bǔ)、多項式等）、相似性（鏡像、周期延拓等）等特征，根據(jù)已有信號兩端的性質(zhì)，實施的信號修補(bǔ)方法。也就是說，我們需要利用信號的一端特性，恢復(fù)出原始信號。相對于前面討論的信號中間修補(bǔ)方法，顯然這里的修補(bǔ)條件要“弱”。

同樣，我們可以利用PDE的信號修補(bǔ)性質(zhì)，實現(xiàn)信號的端點(diǎn)擴(kuò)展處理。此時，只需要定義模板函數(shù)式（12）的延托范圍即可。假定原來信號長度為L0，延托之后信號長度為L1，定義：

其中，L0=x2－x1，即實現(xiàn)在原始信號左右兩端同時擴(kuò)展，左端擴(kuò)展長度為Ll=x1，右端擴(kuò)展長度為Lr=L1－x2。

2.2 試驗分析

對式（14）所示的模擬信號，如圖1（a），令其左右兩端各50個數(shù)據(jù)為0，如圖3所示。圖4－圖8分別表示理想模擬信號［圖1（a）］的 EMD分解結(jié)果，及利用不同擴(kuò)展方法恢復(fù)出左右各50個數(shù)據(jù)后的EMD分析結(jié)果。EMD算法為Rilling等［7］提出的改進(jìn)算法。為便于觀察和比較，圖中同時給出了信號兩端的擴(kuò)展數(shù)據(jù)波形。圖9為式（13）實現(xiàn)的非線性PDE擴(kuò)展及其分解，其中調(diào)制函數(shù)為式（10a），p=－0.01，梯度函數(shù)為式（11a），迭代次數(shù)為1 000。

由圖可以看出，端點(diǎn)擴(kuò)展是信號修補(bǔ)的一種形式，因此關(guān)于信號修補(bǔ)的方法同樣適用于端點(diǎn)擴(kuò)展。但因在擴(kuò)展時，只有一端的特征信息可以利用，因此，常規(guī)的信號修補(bǔ)（擴(kuò)展）方法完成擴(kuò)展后，兩端畸變大，如圖5～圖8。圖8所示SVM擴(kuò)展方法，采用多項式核函數(shù)，多項式次數(shù)為5。實際上，SVM擴(kuò)展方法很不理想，本文試驗了多種核函數(shù)，均無法獲得較好的擴(kuò)展效果。對于傳統(tǒng)方法，當(dāng)擴(kuò)展信號較長時，畸變更為明顯。非線性PDE擴(kuò)展端點(diǎn)畸變少，且適用于較長數(shù)據(jù)的擴(kuò)展。PDE擴(kuò)展之后的分解結(jié)果與理想邊界信號分解結(jié)果一致。

圖10為一數(shù)控車床電主軸低速狀態(tài)下的測試振動信號利用非線性PDE端點(diǎn)處理后的EMD分解結(jié)果。由圖可以看出，端點(diǎn)處理較好的保證了EMD分解的有效性，沒有存在發(fā)散等現(xiàn)象。限于篇幅，本文沒有給出其他端點(diǎn)擴(kuò)展方式下的EMD分解結(jié)果。

圖3 端點(diǎn)擴(kuò)展模擬信號Fig.3 Analog signal for extending end

圖4 理想邊界EMD分解Fig.4 EMD decomposition based on ideal end

圖5 拋物線擴(kuò)展EMD分解Fig.5 EMD decomposition based on extending end using parabola

圖6 拉格朗日多項式擴(kuò)展EMD分解Fig.6 EMD decomposition based on extending end using lagrange polynomial

圖7 三次樣條擴(kuò)展EMD分解Fig.7 EMD decomposition based on extending end using cubic spline

圖8 SVM擴(kuò)展EMD分解Fig.8 EMD decomposition based on extending end using SVM

圖9 非線性PDE擴(kuò)展EMD分解Fig.9 EMD decomposition based on extending end using nonlinear PDE

3 結(jié)論

根據(jù)振動信號的局部連續(xù)性、光滑性特性，可以根據(jù)已有信號實現(xiàn)對缺陷信號的修補(bǔ)。三次樣條修補(bǔ)、多項式修補(bǔ)等方法實際上就是依據(jù)這一原理。熱傳導(dǎo)方程反映了熱傳導(dǎo)這一過程中的連續(xù)性、光滑性等特點(diǎn)，其初值問題的解即為一個高斯濾波過程，因此熱傳導(dǎo)方程可以用于信號的修補(bǔ)。非線性PDE信號修補(bǔ)時，同時還兼顧了信號之間的梯度特征，因此修補(bǔ)后既能保證信號的連續(xù)、光滑性，還能保持其瞬態(tài)性。

EMD方法是一種新的非平穩(wěn)信號分析方法。EMD方法把復(fù)雜信號分解成有限個IMF之和，而這一分解過程中，信號端點(diǎn)的處理是其關(guān)鍵問題之一，直接影響到分解的有效性。常用的端點(diǎn)處理方法實際上就是一種信號的修補(bǔ)方法。但是在端點(diǎn)擴(kuò)展處理時，只能考慮一端對“缺陷”信號的影響，因此利用常用方法處理之后，端點(diǎn)會存在較大的畸變，并且隨著擴(kuò)展信號長度的增加，畸變越來越大。PDE信號端點(diǎn)擴(kuò)展方法是一種新的端點(diǎn)處理方法，可擴(kuò)展復(fù)原出較長時間長度的信號，且處理之后的信號畸變少。

圖10 實際信號PDE擴(kuò)展的EMD分解Fig.10 EMD decomposition of real signal based on extending end using PDE

PDE方法可完成信號的修補(bǔ)，但在這一過程中，仍有一些問題需要進(jìn)一步研究。如怎樣確定PDE修補(bǔ)過程的迭代次數(shù)、怎樣選擇或者設(shè)計調(diào)制函數(shù)、怎樣確定調(diào)制函數(shù)的參數(shù)等問題。

［1］ Tschumperle D，Deriche R.PDE’s based regularization of multi-valued images and applications［J］.IEEE Trans on Pattern Analysis Machine Intelligence，2005 ，27（4）：506－517.

［2］Rudin L，Osher S，F(xiàn)atemi E.Nonlinear total variation based noise removal algorithms［J］.Physica D，1992，60（1 －4）：259－268.

［3］ Weickert J.A review of nonlinear diffusion filtering［J］.Scale-Space Theory in Computer Vision，Lecture Notes in Computer Science，Berlin，1997，1252：3 －28.

［4］吳宏鋼，尹愛軍，秦樹人.基于PDE的振動信號去噪［J］.機(jī)械工程學(xué)報，2009，45（5）：91 －94.

［5］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proc.R.Soc.Lond.A，1998，454：903－995.

［6］胡維平，莫家玲，龔英姬，等.經(jīng)驗?zāi)B(tài)分解中多種邊界處理方法的比較研究［J］.電子與信息學(xué)報，2007，29（6）：1394－1398.

［7］Rilling G，F(xiàn)landrin P，Goncalyes P.On empirical mode decomposition and its algorithms［C］.IEEE-EURASIP Workshop on Nonlinear Signaland Image Processing（NSIP2003），Grado（I），June，2003：8 －11.

［8］劉慧婷，張 旻，程家興.基于多項式擬合算法的EMD端點(diǎn)問題的處理［J］.計算機(jī)工程與應(yīng)用，2004，40（16）：84－86.

［9］鄧擁軍，王 偉，錢成春，等.EMD方法及Hilbert變換中邊界問題的處理［J］.科學(xué)通報，2001，46（3）：257－263.

［10］程軍圣，于德介，楊 宇.基于支持矢量回歸機(jī)的Hilbert_Huang變換端點(diǎn)效應(yīng)問題的處理方法［J］.機(jī)械工程學(xué)報，2006，42（4）：23 －31.

［11］沈 路，周曉軍，張志剛，等.Hilbert-Huang變換中的一種端點(diǎn)延拓方法［J］.振動與沖擊，2009，28（8）：168－171，174.

［12］程軍圣，張 亢，楊 宇，等.局部均值分解與經(jīng)驗?zāi)Ｊ椒纸獾膶Ρ妊芯浚跩］.振動與沖擊，2009，28（5）：13 －16.

［13］ Tschumperle D.PDE's based regularization of multivalued images and applications［D］.Antipolis：University of Nice-Sophia，2002.

［14］ Tschumperl D.Fast Anisotropic smoothing of multi-valued images using curvature-preserving PDE's.International Journal of Computer Vision，2006，68（1）：65－82.