基于張力格林樣條的EMD均值曲線插值方法

戴豪民，許愛強

（海軍航空工程學院 飛行器檢測與應用研究所，山東 煙臺264001）

0 引 言

經驗模式分解 （empirical mode decomposition，EMD）是一種自適應信號分解方法，適用于分析非線性、非平穩信號［1］。經 典EMD 方 法 采 用 三 次 樣 條 插 值 （cubic spline interpolation，CSI）擬合信號的上、下包絡線，從而獲得信號的均值曲線。實驗結果表明，CSI方法擬合的包絡線經常會出現過沖現象，造成均值曲線的誤差較大。目前，針對CSI方法中出現的曲線過沖現象，已經提出了多種改進方法。文獻 ［2］指出三次樣條曲線是二階光滑的 （二階連續可導），過沖的主要原因是曲線的 “柔性”不夠，會在臨近間隔大且缺少約束的插值點之間出現劇烈振蕩 （如圖1所示），遂提出分段冪函數插值法，通過犧牲曲線的光滑性使包絡線更具 “柔性”。其冪指數可以取任意大于1的整數或小數，但隨著冪指數的增大，計算時間也顯著增加。通過實驗驗證，當冪指數β＝2.5時，包絡線同時具備較好的光滑性和柔性，但此時包絡線仍然只有一階光滑性 （一階連續可導）。文獻 ［3］指出過沖的原因是由于CSI方法在待擬合的兩個鄰近插值點之間不具備單調性，提出分段三次 厄 米 特 插 值 法 （piecewise cubic hermite interpolation，PCHI），它能夠使待擬合的兩個鄰近插值點之間具備單調性，較好克服包絡插值過程中出現的過沖現象，但包絡線只有一階光滑性，只能保證各段曲線在連接點處的連續性，不能保證插值曲線在這些點上的光滑性。另外，文獻 ［4］提出B樣條插值法，該方法是為了便于對EMD 分解進行理論分析而引入的一種替代方法，插值特性并沒有顯著提高［5］。

圖1 三次樣條、無張力格林樣條和張力格林樣條對離散數據點插值結果

本文針對EMD 的包絡插值問題，提出采用張力格林樣條 插 值 方 法 （Green’s spline interpolation in tension，GSIT），在不犧牲包絡線光滑性的前提下，減輕包絡線的過沖問題，另外，采用極值中點的一次插值也能進一步降低信號分解的誤差。通過對仿真和實際信號的分析驗證了該方法的有效性，其能夠使EMD 分解更加精確。

1 格林樣條插值原理

1.1 無張力格林樣條插值

格林樣條插值又叫雙調和樣條插值，它是一種基于雙調和算子的全局插值法［6］。在一組不規則的數據點之間，通常可以利用三次樣條插值找出通過這些數據點的光滑曲線。從彈性力學的角度上看，就是求解彈性細桿在不同位置上受到點壓力作用下的彎曲變形問題［7］，這樣受位于xj點的垂直集中力δ（x－xj）作用的彈性細桿，其位移函數g（x）滿足以下雙調和方程

式中：4——雙調和算子，xj——彈性細桿的受力位置，δ（x）——Delta函數。式 （1）的特解是

當對位于xj（j＝1，2，…，N）的N 個數據點wj進行插值時，需要求解的雙調和方程就成為

式中：αj是加權系數，w（xj）是內插曲線在位置xj處的位移。式 （3）和式 （4）的特解w（x）可以表示成各個點壓力的格林函數的線性組合，即

加權系數αj可以通過解下面線性方程組獲得

確定αj后，就可以利用式 （5）計算任一點x處的插值函數值w（x）。

1.2 張力格林樣條插值

無張力格林樣條只受到垂直方向力的作用，水平方向并未受到力的作用。由圖1可見，無張力格林樣條插值會在一些數據點間出現大幅度的振蕩，尤其是數據點間隔較大時，由于水平方向未受力的作用，這種變化更為明顯。通常情況下，這種大幅振蕩現象在三次樣條插值過程中也會出現。為避免出現這種欠合理的現象，可以通過在格林樣條中引入張力參數，從而獲得更為真實的曲線形狀［8］。

當彈性細桿同時受到水平張力p ＞0作用時，其位移函數gp（x）滿足以下雙調和方程

式中：2——調和算子，該方程的一個特解是

通過N 個數據點的受張力p 的位移函數wp（x）為

與式 （6）類似，式 （9）中的加權系數αj可以通過以下線性方程組得到

1.3 計算實例

2 實驗分析

為了驗證格林樣條插值在經驗模式分解中的效果，本文選取EMD 中最為常用的兩種插值方法——三次樣條插值和分段三次厄米特插值與其進行對比，采用正交指數 （index of orthogonality，IO）和能量相對誤差 （energy relative error，ERE）評價各種插值方法的性能

其中，s（t）表示信號，ci（t）表示EMD 分解得到的第i個imf分量。IO 指標越接近于0，說明EMD 分解得到的各個imf分量之間越接近于正交。ERE 指標越接近于0，說明EMD分解能量損失越低。

2.1 仿真信號分析

考察以下多分量信號：s（t）＝sin（20πt）＋sin（50πt）＋sin（140πt）t∈［0，1）求取其極大值后，采用3種方法擬合的該信號上包絡線，如圖2所示 （截取信號中0.16～0.23段表示，張力參數取0.5），從圖2中A 點可以看出，CSI方法和GSIT 方法擬合得到的包絡線在此處存在過沖現象，過沖程度GSIT 方法要輕于CSI方法，而PCHI方法則不存在過沖問題；在B 點，由于PCHI方法要保證插值曲線的單調性，3種方法均出現過沖現象，其過沖程度PCHI方法最輕，GSIT 方法次之，CSI方法最差。從信號擬合的包絡線來看，PCHI方法只存在一階導數，其光滑性最差，后續的分解情況會發現，雖然PCHI方法過沖現象最不明顯，但其分解誤差卻是最大的，這主要是由于包絡線的光滑性差造成的；GSIT 方法存在四階導數，在保證包絡線光滑性的情況下，一定程度上克服了三次樣條插值過程中出現的大幅度的振蕩現象，從而降低插值誤差。

圖2 三次樣條插值、分段三次厄米特插值以及張力格林樣條插值擬合仿真信號上包絡結果

從表1可以看出，GSIT 方法無論在IO 指標還是ERE指標都要優于CSI和PCHI方法；PCHI方法雖然能在插值過程中的在一定程度上克服過沖問題，但其插值曲線光滑性差的特點使得分解誤差較大。

表1 3種插值方法對仿真信號進行EMD分解的IO 和ERE指標比較

2.2 直接擬合均值曲線的格林樣條插值方法

EMD 算法的關鍵是求取均值曲線，均值曲線生成的不準確將導致信號的EMD 分解結果不準確。經典經驗模態分解采用三次樣條插值方法進行上、下包絡擬和，來求取信號的均值曲線。但是這種方法求取均值曲線需要上、下包絡兩次擬合，加之擬合曲線和實際包絡會存在偏差，兩次擬合更會加大這種偏差［10］，隨著迭代的進行，這種偏差會不斷積累，最終會使信號分解不能得到正確的結果。

為了克服上述這些缺點，本文提出一種基于極值中點的格林樣條插值方法來直接擬合均值曲線。該方法通過以下三步來實現：

（1）采用平行延拓法［11］分別得到信號x（t）左右兩端的中點。

（2）找出信號x（t）上所有相鄰極值點連線的中點，與信號兩端的中點組成信號的中點序列。

（3）用張力格林樣條函數擬合中點序列作為均值曲線。

這種擬合均值曲線的方法，由于利用了信號的所有極值點，其擬合精度也會更高。為了驗證該方法的有效性，以上述仿真信號為例，比較一下兩次包絡插值擬合的均值曲線和本文提出的方法 （張力參數取0.5）的差異，其差異性仍采用IO 和ERE指標。

從圖3可以看出，通過上下包絡兩次插值求出的均值曲線與直接張力格林樣條插值形成的均值曲線，差異性還是很明顯的。通過均值曲線穿過中點的情況來看，后者的誤 差 還 是 很 小 的，其IO 和 ERE 指 標 僅 為0.0181和0.0223。

2.3 實際信號分析

圖4是電機軸承加速度數據，軸承型號是6326－C3，采樣頻率fs＝25600 Hz，采樣時間1S。

圖5分別是CSI、PCHI方法和本文方法擬合的均值曲線 （截取信號中0.022s～0.225s 段表示，張力參數取0.6）。CSI方法在部分相鄰極大值點間出現劇烈振蕩現象，其均值曲線與極值中點偏差較大。PCHI方法曲線光滑性較差，擬合的效果一般。GSIT 方法直接對極值中點進行插值，擬合的均值曲線與實際信號均值曲線最為接近。從表2可以看出，PCHI方法的IO 指標高達0.2561，分解出的固有模態函數正交性很差，這說明了PCHI方法為了保持相鄰極大值點 （或極小值點）間插值曲線的單調性，失去了信號波動的本質，從而導致分解誤差較大。

圖3 三次樣條插值、分段三次厄米特插值和本文方法擬合仿真信號均值均線的結果

圖4 電機軸承加速度時域波形

圖5 三次樣條插值、分段三次厄米特插值和本文方法擬合實際信號均值均線的結果

表2 3種插值方法對實際信號進行EMD分解的IO 和ERE指標比較

3 關于張力參數的選取問題

格林函數的加權系數αj是通過求解數據點的線性方程組獲得，線性方程組的穩定性很大程度上會受到數據分布的影響。當數據的最大與最小距離值之比很大時，線性系統往往會不穩定［12］。如果數據點分布比較均勻時，線性系統會趨于穩定。通過實驗發現，張力參數的取值對線性系統的穩定性有很大影響。當張力參數取值增大時，線性系統特征值的衰減速率明顯變慢，系統會更加穩定，但t的取值不能無限增大，當t＝1時，格林函數插值就蛻變為線性插值。所以張力參數的取值一般設定為0.5≤t≤0.7，這樣既可以滿足分解誤差要求，又可以兼顧系統穩定性要求。圖6是張力參數t＝0.1和t＝0.7時特征值的衰減曲線，其特征值的最小幅度相差兩個數量級。

圖6 特征值衰減速度在張力參數為t＝0.1和t＝0.7的比較

4 結束語

針對現有CSI方法存在的問題，本文提出了一種基于張力格林樣條的均值曲線插值方法，它在保證插值曲線光滑性的基礎上，可以在一定程度上克服CSI方法出現的過沖現象，同時基于極值中點的一次插值也能進一步降低信號分解的誤差。仿真信號和實際信號的分析結果也說明了該方法的有效性。但是，對于較長的復雜數據，會碰到大量非規則分布的極值點，由于極值點數增加，對應加權參數的線性方程組求解需要較長時間。如果引入分段格林樣條插值［13］，不僅可以減少計算時間，同時也會改變極值點分布的均衡性，有利于提高插值結果的穩定性，這也是后續工作的一個重點。

［1］Huang N E，Shen Z，Long S R.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non－stationary time series analysis ［J］.Proceeding of Royal Society，1998，454 （1971）：903－995.

［2］ZHONG Youming，JIN Tao，QIN Shuren.New envelope algorithm for Hilbert－Huang transform ［J］.Journal of Data Acquisition ＆Processing，2005，20 （1）：13－17 （in Chinese）.［鐘佑明，金濤，秦樹人.希爾伯特黃變換中的一種新包絡線算法 ［J］.數據采集與處理，2005，20 （1）：13－17.］

［3］LONG Sisheng，ZHANG Tiebao，LONG Feng.Causes and solutions of overshoot and end swing in Hilbert－Huang transform ［J］.Acta Seismologica Sinica，2005，27 （5）：561－568 （in Chinese）.［龍思勝，張鐵寶，龍峰.希爾伯特—黃變換中擬合過沖和端點飛翼 的 原 因 及 解 決 辦 法 ［J］.地 震 學 報，2005，27 （5）：561－568.］

［4］Chen Qiuhui.A B－spline approach for empirical mode decompositions［J］.Advances in Computational Mathematics，2006，24 （1－4）：171－195.

［5］ZHU Weifang，ZHAO Heming，CHEN Xiaoping.A leastlength constrained envelope approach for EMD ［J］.Acta Electronica Sinica，2009，40 （9）：1909－1912 （in Chinese）. ［朱偉芳，趙鶴鳴，陳小平.一種最小長度約束的EMD 包絡擬合方法 ［J］.電子學報，2009，40 （9）：1909－1912.］

［6］Wessel P.Interpolation using ageneralized Green’s function for a spherical surface spline in tension ［J］.Computer ＆Geosciences，2009，36 （6）：1247－1254.

［7］XIA Jizhuang.An application of biharmonic spline interpolation to integration of logging and seismic data［J］.Oil ＆Gas Geology，2010，31 （2）：244－248 （in Chinese）.［夏吉莊.雙調和樣條內插方法在測井和地震資料整合中的應用 ［J］.石油與天然氣地質，2010，31 （2）：244－248.］

［8］XU Feng，HE Rizheng，ZHANG Guibin.Green’s functionbased interpolator and its application ［J］.Progress in Geophysics，2013，28 （4）：1721－1728 （in Chinese）. ［許豐，賀日政，張貴賓.格林樣條插值算法及其應用 ［J］.地球物理學進展，2013，28 （4）：1721－1728.］

［9］Wessel P，Becker JM.Interpolation using ageneralized Green’s function for a spherical surface spline in tension ［J］.Geophysical Journal International，2008，174 （1）：21－28.

［10］ZHU Zheng，WANG Yilin，CAI Ping.A method of obtaining the mean curve for empirical mode decomposition ［J］.Journal of Harbin Engineering University，2011，32 （7）：872－876 （in Chinese）.［朱正，王逸林，蔡平.EMD 中的求取均值曲線方法 ［J］.哈爾濱工程大學學報，2011，32 （7）：872－876.］

［11］Wu Z，Huang NE.Ensemble empirical mode decomposition：A noise assisted data analysis method ［J］.Advances in Adaptive Data Analysis，2009，1 （1）：1－41.

［12］Wessel P，Bercovici D.Interpolation with splines in tension：A Green’s function approach ［J］.Mathematical Geology，1998，30 （1）：77－94.

［13］DENG Xingsheng，TANG Zhong’an.Study on two dimensional spline interpolation based on moving Green’s function［J］.Journal of Geodesy and Geodynamic，2011，31 （6）：69－72 （in Chinese）. ［鄧興升，湯仲安.移動格林基函數樣條二維插值算法研究 ［J］.大地測量與地球動力學，2011，31 （6）：69－72.］