PDE信號修補方法及在EMD端點效應處理中的應用

尹愛軍，王 璇

（重慶大學 機械傳動國家重點實驗室，重慶大學 機械學院測試中心，重慶 400044）

PDE是一個含有多元未知數函數及其偏導數的方程。PDE是在討論自然現象（特別是物理現象）的過程中逐步建立起來的，所以也稱為數學物理方程。PDE最新的研究與應用已擴展到經濟、金融預測、圖像處理、振動信號處理等領域［1－4］。

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

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

其中，

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

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

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

對于每一個IMF，應滿足：① 在整個數據段內，極值點的個數和零交叉點的個數必須相等或相差最多不能超過一個。② 在任何一點，由局部極大值點形成的包絡線和由局部極小值點形成的包絡線的平均值為零。

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

本文針對振動信號的內在性質，分析基于插補等方法的信號修補方法，提出了PDE信號修補新方法，并將這一方法成功用于EMD分解過程中的端點效應處理。

1 基于PDE的波動信號修補

1.1 信號修補方法

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

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

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

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

1.2 PDE信號修補原理

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

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

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

即對于畸變缺陷信號，根據熱傳導過程中的連續及光滑性，逐漸重構出原始信號；而對于正常信號，則不作處理。

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

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

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

在這一修補方式下，可以針對信號梯度設置相應的調制參數，使修補后的信號即滿足整體連續光滑、又滿足瞬態沖擊特性。

若定義如下模板函數：

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

1.3 仿真試驗

設由兩個正弦信號組成的理想信號，如式（14）。設置采樣頻率為1 000 Hz，取連續1 024個數據點（y1，y2，…，y1024），如圖 1（a）所示。畸變信號為y50，y51，…，y100，如圖 1（b）所示。

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

圖2 信號修補Fig.2 Signal restored

2 EMD端點效應的PDE處理

2.1 PDE端點效應處理

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

同樣，我們可以利用PDE的信號修補性質，實現信號的端點擴展處理。此時，只需要定義模板函數式（12）的延托范圍即可。假定原來信號長度為L0，延托之后信號長度為L1，定義：

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

2.2 試驗分析

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

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

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

圖3 端點擴展模擬信號Fig.3 Analog signal for extending end

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

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

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

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

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

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

3 結論

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

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

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

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

［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，Fatemi 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］.機械工程學報，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］.電子與信息學報，2007，29（6）：1394－1398.

［7］Rilling G，Flandrin 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端點問題的處理［J］.計算機工程與應用，2004，40（16）：84－86.

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

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

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

［12］程軍圣，張 亢，楊 宇，等.局部均值分解與經驗模式分解的對比研究［J］.振動與沖擊，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.