陳道云，孫守光，李 強

（北京交通大學　機械與電子控制工程學院，北京100044）

基于經(jīng)驗模態(tài)分解及時域積分的軌道高低不平順獲得方法

陳道云，孫守光，李 強

（北京交通大學機械與電子控制工程學院，北京100044）

分析了信號中的趨勢項產(chǎn)生原因，并推導了時域積分過程中趨勢項的產(chǎn)生過程；利用經(jīng)驗模態(tài)分解的方法提取并去除了實測列車軸箱垂向加速度信號的趨勢項，然后對去趨勢項的加速度信號在時域內(nèi)連續(xù)兩次積分得到軌道高低不平順數(shù)據(jù)；分別用經(jīng)驗模態(tài)分解法和最小二乘法去除了積分過程中產(chǎn)生的趨勢項，并將兩種方法的處理結(jié)果進行了對比分析。分析結(jié)果表明，利用最小二乘法去除信號及積分過程中產(chǎn)生的趨勢項存在較大的計算誤差，而利用經(jīng)驗模態(tài)分解的方法則可以有效去除信號及積分過程產(chǎn)生的趨勢項，從而得到軌道高低不平順數(shù)據(jù)。

加速度；經(jīng)驗模態(tài)分解；最小二乘法原理；趨勢項；時域積分；軌道高低不平順

軌道不平順對機車車輛系統(tǒng)是一種外部激擾，是產(chǎn)生機車車輛系統(tǒng)振動的主要根源。按其對車輛激擾作用的方向可分為垂向軌道不平順和橫向軌道不平順。左右軌垂向、橫向不平順又可組合成常說的高低、水平、軌向及軌距4種不平順。

目前軌道不平順的獲取主要依靠軌檢車，系統(tǒng)很復雜。為了簡單地獲得軌道不平順數(shù)據(jù)，本文采用了一種基于經(jīng)驗模態(tài)分解及時域內(nèi)軸箱加速度積分的方法，該方法易于實現(xiàn)且精度較高，比較適合對軌道進行簡易評估。

對于加速度信號的積分，之前有大量的科研工作者對其進行了不同形式的研究，而對于趨勢項的消除方法，普遍采用的是最小二乘法原理［1-2］，這種方法對處理平穩(wěn)且數(shù)據(jù)量小的信號具有很大的優(yōu)勢。對于本文中研究的非平穩(wěn)隨機的加速度信號，采用這種方法得到的結(jié)果往往是不準確的。為此，本文利用經(jīng)驗模態(tài)分解技術對加速度信號及時域積分產(chǎn)生的趨勢項進行了去除，得到了軌道高低不平順數(shù)據(jù)。

1　趨勢項

1.1信號本身的趨勢項

實際采集的加速度振動信號，由于放大器隨著溫度變化產(chǎn)生的零點漂移、傳感器頻率范圍外低頻性能的不穩(wěn)定以及傳感器周圍的環(huán)境干擾，往往會偏離基線，甚至偏離基線的大小還會隨著時間變化。偏離基線隨時間變化的過程稱為信號的趨勢項，它的存在會直接影響到信號的正確性，因此必須將其去除。

1.2時域積分及其趨勢項

在時域內(nèi)對軸箱加速度信號兩次積分是獲得軌道高低不平順的最直接方法。數(shù)值積分的方法有很多，如梯形公式、Sim pson公式、Newton-Cotes公式、復化梯形公式、復化Sim pson公式等［3］。為了提高積分的精度且方便運算，本文采用復化Sim pson公式對加速度信號進行兩次積分，分別得到速度信號和位移信號：

式中

然而，直接通過時域積分的方法獲得的軌道高低不平順數(shù)據(jù)是不準確的，這種偏差的原因正是由于趨勢項的存在。

假設加速度信號中含有微小直流分量，即

一次積分得：

二次積分得：

式中δ，η分別為ε一次、二次積分后產(chǎn)生的積分常量。

由式（3）～式（5）可以看出，在一次、二次積分后分別產(chǎn)生了一次、二次趨勢項，因此應予以去除。

1.3趨勢項的去除

（1）基于經(jīng)驗模態(tài)分解的趨勢項去除

1998年美國宇航局的N orden E H uang等［4］提出了一種新的信號處理方法，該方法主要由兩部分組成：首先采用經(jīng)驗模態(tài)分解（E m pirical M ode Deco m position，以下簡稱E M D）將原始時間信號分解成一組本征模函數(shù)（Intrinsic M ode Function，以下簡稱I M F），然后再對每一個I M F進行Hilbert變換。

本文利用其中的E M D方法提取并去除信號中的趨勢項。算法如下：

（1）找出信號x（t）的所有極大點和極小點，用3次樣條曲線分別擬合為原數(shù)據(jù)序列的上包絡線U（t）和下包絡線L（t），上、下包絡線的均值

將原數(shù)據(jù)序列x（t）減去m1（t）可得到一個去掉低頻的新數(shù)據(jù)序列h1（t），其式為

判斷h1（t）是否為I M F的條件有兩個：

①對于一列數(shù)據(jù)，極值點和過零點數(shù)目必須相等或至多相差一點；

②在任意點，由局部極大點和極小點構(gòu)成的兩條包絡線的平均值為零。

如果h1（t）不滿足上述條件，那么將h1（t）看成，則

重復以上過程m次，直到所得的h1m（t）滿足I M F所必需的條件，此時h1m（t）就是第一個I M F，記為I M F1（t），它表示信號中的最高頻成分。

（2）用x（t）減去I M F1（t）得到一個去掉高頻成分的新數(shù)據(jù)序列r1（t），其式為

將r1（t）看作是原數(shù)據(jù)序列重復步驟（1），可得到一系列I M Fi（t）和最后一個不可分解的序列rn（t）。由此，原數(shù)據(jù)序列x（t）可表示為一組I M F分量和一個殘余分量的和：

E M D分解的收斂特性表明，分解得到的殘余信號分量rn（t）為單調(diào)函數(shù)，其中包含了測試信號中頻率最低的成分，其周期大于采樣信號的長度，所以rn（t）即為測試信號中包含的趨勢項［5］。

（2）基于最小二乘法原理的趨勢項去除

信號測試的采樣頻率記為fs，總采樣時間記為T，則采樣時間間隔Δt=1/fs，實測的加速度信號記為a（ti），其中：

構(gòu)造一個多項式函數(shù)

其中Pk為多項式系數(shù)；φ為所有次數(shù)不超過n（N≤T/ Δt-1）的多項式構(gòu)成的函數(shù)類，使得函數(shù)λn（ti）與離散數(shù)據(jù)a（ti）的誤差平方和最小，即：

若想確定擬合趨勢項，只需確定某一組Pk的值，使得E最小。由多元函數(shù)求極值的必要條件：

其中，j=0，1，2，…，n因而可以得到：

用矩陣表示為

可以證明，方程組的系數(shù)矩陣是一個對稱正定矩陣，故存在惟一解。

2　加速度信號測試及去除趨勢項處理

2.1測試情況

2014年11月12日對尚未開通的滬昆高鐵線路上運行的CRH3 A型動車組的軸箱垂向加速度等參數(shù)進行了線路實測，測試區(qū)間為南昌西至玉山南。

圖1　CRH3 A型動車組

測試選用IC壓電式單向加速度傳感器，其靈敏度為25.1 m V/g，頻率范圍為0.7～11 000 Hz，量程為100g，布置在1-2L號軸箱與一系減振器之間。采用德國IM C數(shù)據(jù)采集系統(tǒng)對加速度信號進行采集，如圖2所示。

2.2加速度信號去除趨勢項

在車輛運行平穩(wěn)的區(qū)段，截取1 s時長的加速度測試數(shù)據(jù)，觀察加速度值隨時間的變化情況，如圖3所示。

圖2　加速度傳感器安裝位置及數(shù)據(jù)采集系統(tǒng)

圖3　含有趨勢項的加速度信號

現(xiàn)有的消除趨勢項的方法大多需要預先假定信號中趨勢項的類型，這就需要對測試信號中包含趨勢項的特征具有一定的預先驗證知識，因而不適用于處理隨機信號。

本文利用E M D算法，通過M atlab編程對上述加速度信號進行分解，所得的I M F及其對應的功率譜（PSD）分析如圖4所示。

對圖4進行分析可以發(fā)現(xiàn)：

（1）軸箱加速度信號被成功分解為I M F1～I M F10及Trend共11個分量。I M F1頻率最高，波長最短，在隨后的分量中，頻率逐漸變低，波長變長，其中Trend的頻率最低。各個I M F分別以不同的分辨率體現(xiàn)了信號中該頻段范圍內(nèi)的振動模態(tài)。

（2）軸箱加速度信號的頻譜比較豐富，大部分信號分量集中在1 000 H z以下。I M F1～I M F6屬于測試信號的優(yōu)勢頻段，體現(xiàn)了軸箱垂向振動的主要時頻特征。I M F7～I M F10屬于測試信號中包含的低頻成分。余量Trend頻率很低且所占能量很大，其為單調(diào)函數(shù)，周期大于采樣信號的長度，因此屬于低頻趨勢項，須將其去除。

3　軌道高低不平順數(shù)據(jù)的獲取

在時域內(nèi)，對去除趨勢項的軸箱加速度信號連續(xù)兩次積分便得到軌道的高低不平順，如圖5所示。

可見時域內(nèi)積分獲得的位移激擾含有明顯的趨勢項，這是因為在去除了加速度信號的趨勢項之后，新的加速度信號中不可避免地殘留有微小的直流分量，在積分過程中，這種微小的直流分量被逐漸放大，以至于偏離正確值，因而需對其進行去除趨勢項的處理?，F(xiàn)分別利用E M D方法和最小二乘法，編寫M atlab程序?qū)Ψe分得到的位移信號去除趨勢項。

3.1E M D方法去除積分趨勢項

利用E M D原理編寫M atlab程序?qū)ι鲜鑫灰菩盘栠M行分解，所得的I M F 分量及其對應的功率譜（PSD）分析如圖6所示。

對圖6進行分析，發(fā)現(xiàn)經(jīng)E M D分解得到的4個分量中，Trend的頻率最低且所占能量較大，其為單調(diào)函數(shù)，周期大于采樣信號的長度，因此為趨勢項，須將其去除。去除趨勢項后的位移信號如圖7所示。

3.2最小二乘法去除積分趨勢項

利用前面提到的最小二乘法也可以去除由時域積分產(chǎn)生的趨勢項，去除趨勢項后的位移數(shù)據(jù)如圖8所示。

3.3兩種去除趨勢項方法的對比

首先將E M D得到的趨勢項和最小二乘法擬合的二次趨勢項進行對比，如圖9所示。

可以看出，兩條擬合趨勢曲線的大致走向是一致的，相比最小二乘法擬合趨勢項，E M D分解得到的趨勢項波動更明顯一些，這與E M D分解的本質(zhì)特征有關，因為這種分解完全依據(jù)數(shù)據(jù)自身的時間尺度特征來進行信號分解，具有自適應特征，因而提取趨勢項的精度非常高，而利用最小二乘法擬合的趨勢項則是一條標準的二次函數(shù)曲線，與E M D相比，精度偏低。

然后將兩種方法去除趨勢項后得到的位移數(shù)據(jù)進行對比，如圖10所示。

觀察圖10可以發(fā)現(xiàn)，通過最小二乘法去除趨勢項得到的位移信號波動幅度要明顯大于通過E M D方法去除趨勢項得到的位移信號，但是二者的波形走勢卻基本相同。導致前者波動幅度大的原因是利用最小二乘法原理去除趨勢項本身存在一定缺陷。

圖4　I M F及相應功率譜

圖5　含趨勢項的位移信號

3.4最小二乘法去除積分趨勢項的缺陷

在時域內(nèi)對積分后的數(shù)據(jù)去除趨勢項時，當截取起點相同而長度不同的幾組數(shù)據(jù)時，會導致在數(shù)據(jù)相同部分得到的位移修正曲線隨著截取數(shù)據(jù)長度的增加而逐漸發(fā)生偏移。這是因為當截取的數(shù)據(jù)長度不同時，利用最小二乘法得到的擬合多項式系數(shù)便會不同，從而導致擬合多項式的值也發(fā)生變化。

圖6　I M F分量及相應功率譜

圖7　EMD去除趨勢項得到的位移

圖8　最小二乘法原理去除趨勢項得到的位移

圖9　兩種方法得到的趨勢項對比

圖10　兩種方法去除趨勢項后得到的位移對比

例如若要去除一階趨勢項，則根據(jù)式（15）得擬合多項式系數(shù)

可見，當截取的數(shù)據(jù)點數(shù)不同時，T/Δt的值便會發(fā)生變化，因而擬合的一次式系數(shù)P0和P1也不同，在起點相同終點不同的各組積分區(qū)間內(nèi)，數(shù)據(jù)相同部分得到的位移修正曲線便會不同。這也是利用最小二乘法去除趨勢項的一個很大的缺陷，這種缺陷正是算法本身的缺陷造成的。

為了用試驗數(shù)據(jù)驗證這種積分誤差，分別取兩個時間區(qū)間［0，0.2］，［0，1］對加速度數(shù)據(jù)兩次積分并去除二階趨勢項，得到兩條位移曲線，為了便于觀察，取相同時間區(qū)間［0，0.2］作圖如圖11。

圖11　不同長度積分區(qū)間對相同時間段積分結(jié)果影響

從圖中可以看到，隨著積分區(qū)間長度的增大，在相同區(qū)間段［0，0.2］的位移值發(fā)生偏移，這正是因為增加了數(shù)據(jù)點數(shù)而導致對前面的擬合帶來了影響。為了消除這種影響，可以采用動態(tài)遞推法［6］，步驟如下：

采用某段長度的加速度信號進行積分和擬合修正，將修正位移結(jié)果第一個點的值作為最終位移的第一個值，將該段長度的數(shù)據(jù)點向后遞推一個采樣點進入下一次循環(huán)處理，各次循環(huán)處理得到的修正位移的第一個值組成最終位移值。

雖然動態(tài)遞推法可以有效消除后面數(shù)據(jù)對前面擬合影響帶來的誤差，但是該方法的關鍵在于找到一段積分結(jié)果準確的數(shù)據(jù)區(qū)間，這樣才能利用動態(tài)遞推的原理消除因不同數(shù)據(jù)區(qū)間長度而引起的積分結(jié)果偏移。對于軸箱加速度信號，如何找到這一個積分結(jié)果準確的基準積分區(qū)間還有待進一步研究。

4　結(jié) 論

（1）趨勢項普遍存在于采集得到的信號以及時域積分過程中，運用E M D法可有效地消除實測信號以及時域積分產(chǎn)生的趨勢項，使通過軸箱加速度信號積分得到的軌道高低不平順數(shù)據(jù)更加準確。

（2）利用最小二乘法在時域內(nèi)對積分后的數(shù)據(jù)去除趨勢項時，當截取起點相同而長度不同的幾組數(shù)據(jù)時，會導致在相同時間區(qū)間部分得到的位移修正曲線隨著截取數(shù)據(jù)長度的增加而逐漸發(fā)生偏移。雖然采用動態(tài)遞推法可以消除這種偏移，但是如何選擇合理的基準積分區(qū)間有待研究。

［1］ 陳為真，汪秉文，胡曉婭.基于時域積分的加速度信號處理［J］.華中科技大學學報，2010，38（1）：2-4.

［2］ 李東文，熊曉燕，李 博.振動加速度信號處理探討［J］.機電工程技術，2008，37（9）：50-52.

［3］ 王兵團，張作泉，趙平福.數(shù)值分析簡明教程［M］.北京：清華大學出版社，北京交通大學出版社，2012.

［4］ H uang N E，Shen Z，Long S R，et al.The e-mpirical mode decomposition and the Hi lbe-rt spectrum for nonl inear and non-stationa-ry time series analysis［J］.Proceedings of the Royal Society of London，Series A，1998，454：903-995.

［5］ Rilling G，F(xiàn)landrin P.O n em pirical m ode deco m position and its algorith ms［C］∥Grado.IE E E-E U R A SIPW orkshop on Nonlinear Signal and Image Processing.N SIP-03，June 2003.

［6］ 馬建倉，劉 琦，程存虎，等.依據(jù)振動加速度信號繪制航空發(fā)動機轉(zhuǎn)子軸心軌跡的方法［J］.航空計算技術，2009，39（4）：10-14.

An Acquired M ethod of Track Vertical Profile Irregularity Based on E M D and Time Domain Integral

C H E N D aoyun，S U N Shouguang，LI Qiang
（School of M echanical，Electronic and Control Engineering，Beijing Jiaotong U niversity，Beijing 100044，China）

This paper analyzed the cause of trend item in signals and derived the generation process of trend item in time do main integral. Then it extracted the trend item of actual measured vertical acceleration signals of the axle box and rem oved it.After that，the new acceleration signals were treated by twice-integration in time do main，w hich finally turned into track vertical profile irregularity data.This paper rem oved the trend item produced in the process of time do main integral by using E m pirical M ode Deco m position（E M D）and least square principle and contrasted the treatment results of two methods.The results showed that big calculation errors will happen by using least square principle to rem ove trend item in signals and process of time do main integral.H owever，E M D can rem ove trend item effectively and guarantee the accuracy of track vertical profile irregularity data.

acceleration；E M D；least square principle；trend term；time do main integral；track vertical profile irregularity

U26011+1

A

10.3969/j.issn.1008-7842.2015.05.03

1008-7842（2015）05-0009-06

陳道云（1988—）男，博士研究生（2015-04-27）