基于改進的EMD方法提取車輛-軌道垂向耦合系統動態特性

陳雙喜，林建輝，陳建政

(西南交通大學 牽引動力國家重點實驗室，成都 610031)

輪軌系統的激勵可以分為脈沖型、諧波型和動力型，由車輛和軌道兩方面的因素造成。車輛方面較為單一，主要是車輪擦傷、車輪踏面幾何不平順及車輪偏心。而軌道方面的因素比較復雜，大體包括三個方面［1］:①軌道幾何不平順，如高低、水平、方向不平順以及軌面的波浪形磨損;②鋼軌接頭不良;③軌下基礎缺陷。軌面波浪型磨損是鋼軌沿縱向表面出現的周期性波浪性不平順，一般波長在25 mm～1 500 mm，波深0.1 mm～2 mm，其產生機理至今不明［2-4］。軌道的幾何不平順包括靜態不平順和動態不平順，一般用多諧波函數表示［5］。由于軌道不平順、波浪形磨損等因素的影響，高速列車運行時輪軌間將產生接觸振動，甚至導致高頻共振，對車輛和軌道使用壽命影響極大。通過動力學仿真計算出各種不同激勵模型下的動力學響應，國內外學者對此進行了大量的研究，但是對復雜激勵模型下車輛-軌道耦合系統動力學響應的動態特性提取主要以頻譜分析為主，研究并不深入。而此方面的研究有助于對車輛-軌道系統動力學特性及其成因進行深入了解，具有重要的研究意義。

本文提出利用一種改進的經驗模態分解方法對車輛-軌道垂向耦合系統的振動響應進行經驗模態分解，提取耦合系統的各種動態特性，并且對垂向輪軌接觸力、轉向架和車體的加速度本征函數(IMF)進行了分析和比較，取得了較好的效果。

1 改進的EMD算法

經驗模式分解(EMD)是美國華裔科學家Huang提出的新的信號時頻分析方法。它將信號自適應地分解為多個本征函數(IMF)及一個余項，從而反映信號的內部特點。滿足這樣兩個條件的信號稱為本征函數(IMF):① 在整個數據序列中，極值點的數量(包括極大值點和極小值點)與過零點的數量必須相等，或最多相差不多于一個;②在任一時間點上，信號局部極大值確定的上包絡線和局部極小值確定的下包絡線的均值為零。EMD方法對非平穩、非線性信號處理具有較高的效率，但其局部均值的精度較低，計算速度較慢。此外EMD方法在對非平穩信號進行分解時，在數據的兩端還會產生發散現象，并且這種發散的結果會逐漸向內“污染”整個數據序列而使所得分解結果嚴重失真(端點效應)。因此，有兩個問題需要解決:一是如何進一步提高局部均值的求解精度，二是如何有效地消除因邊界不連續而產生的邊界效應。

在EMD方法中以局部極大值與極小值的包絡線的均值代替信號局部均值并不是唯一的方法，其他求解方法有［6，7］:自適應時變濾波法(ATVFD)、極值域均值模式分解法(EMMD)、改進的極值域均值模式分解法(IEMMD)等。本文以改進的極值域均值模式分解方法來提高局部均值的求解精度。其步驟如下:

① 求出原始數據x(t)中所有局部極值點組成極值點序列e(ti)，其中i=1，2，3…k。按照公式(1)計算出兩相鄰極值點間的局部均值序列mi:

② 假設mi在原始數據x(tii)與x(tii+1)(1≤ii≤k-1)之間。按式(2)求得mi對應的時間tγi:

③ 用兩個相鄰的局部均值mi(tγ1)和mi+1(tγ2)加權平均求ti+1處極值點的局部均值m(ti+1)，即:

公式(3)中的k(ti)和k(ti+1)是通過相似梯形得到的加權系數，即:

求得極值點處的局部均值之后，就可以用這些點來擬合數據的局部均值曲線，進而分解出IMF。

目前已經提出了一些抑制端點效應的方法［6-9］:如采用神經網絡對數據延拓法、在端點處按照端點數據變化的“平衡位置”附加兩條平行線段的方法、極值點延拓法、基于AR模型的時間序列線性預測方法等。邊界波形匹配預測法［10］能最大限度地維護原始信號在端部的變化趨勢，相對與其它方法而言，其構造的虛擬極值點更加合理，對端點效應的抑制效果更加理想。考慮到計算的精度與成本和適應性，本文以邊界波形匹配預測法來抑制端點效應。其基本步驟如下:

① 分別確定數據x(t)左右端的兩個相鄰極值點。假設左端兩個極值點為e0和e1，右端極值點為en-1和en，左端點到e0這段波形記作w0，長度為l0，右端點到en-1這段波形記作wn，長度為ln。

② 計算整個數據段每一個極值點附近與w0(wn)長度相同的一段波wi與w0(wn)的匹配度mati。

s(j)表示波w0(wn)所有數據點，u(j)表示波形匹配計算中wi平移后的波形，N為序列長度。

③記mat=min(mati)，如果mat＜αl(其中α是一個常數，l為波的長度)則選取wi作為原始數據的延拓，否則選取原始數據兩端的兩個相鄰的極值點的均值為兩端的極值點。

2 車輛－軌道耦合動力學模型

車輛-軌道垂向耦合動力學模型如圖1所示，包括車輛模型、輪軌耦合模型、軌道模型。其中輪軌耦合模型采用赫茲非線性接觸理論實現車輛系統與軌道系統的耦合，其計算方法見文獻［1］。

圖1 車輛-軌道耦合系統模型Fig.1 Model of vehicle-track coupling system

2.1 車輛模型

本文建立10自由度客車垂向振動模型，包括車體和轉向架的沉浮與點頭、輪對的沉浮。車輛系統的振動響應微分方程為:

2.2 軌道模型

對于軌道系統的模型，本文選取彈性支承塊式無砟軌道作為研究對象。文獻［11］計算得知，鋼軌的剪切變形和轉動慣量對鋼軌的固有頻率有影響，模態越高，影響越明顯。對于高速高頻激勵響應，只有鋼軌采用Timoshenko梁模型，考慮鋼軌斷面的剪切變形，才能真實反映輪軌系統的高頻振動。為了滿足高速高頻振動要求，本文選取鋼軌模態階數為NM=90，鋼軌長度為120 m。于是，鋼軌模型微分方程［11］為:

式中，A為鋼軌斷面面積;mr為單位長度鋼軌質量;E為鋼軌彈性模量;I為鋼軌截面慣量;k為鋼軌截面影響系數;Φr(x，t)為鋼軌軸線的轉角;zr(x，t)為鋼軌的振動位移;G為鋼軌剪切模量;Fj為輪軌作用力;Ri為支承塊反作用力。

彈性支承塊的振動微分方程［1］為:

式中，Msi為彈性支承塊質量;Zsi為支承塊位移;Kbi、Kpi為塊下膠墊和扣件的剛度;Cbi、Cpi為塊下膠墊和扣件的阻尼。

2.3 軌面激勵模型

對于高速列車，整個車輛-軌道耦合系統振動的原因可能是車輪缺陷也可能是軌道不平順。鋼軌波浪形磨損是鋼軌沿縱向表面出現的周期性的類似波浪形狀的不平順現象。Zarembski和 Grassie［2，3］按照波長確定機理和損傷機理對其進行了分類，指出低頻力P2共振導致鋼軌塑性彎曲、滾動接觸疲勞和波谷縱向振動磨損，從而分別形成波長為500 mm～1 500 mm、150 mm～450 mm、150 mm～450 mm的波磨。

本文以諧波波磨激勵模型進行計算，諧波激勵模型如文獻［1］。

單一諧波激勵模型:

組合諧波激勵模型:

公式中a為波的深度;L為波長;n為波數。

根據上述的波磨特點，本文的諧波激勵模型假設為兩種不同波長(波長分別為L1=300 mm、L2=1 500 mm)的鋼軌面波浪形磨損疊加在軌道垂向諧波不平順(波長L3=6 m)上，即激勵為三種不同波長的諧波組合。

2.4 系統動力學特性的分析

為了研究車輛-軌道耦合系統動力學的時域特性，需要在時域上對由方程(5)～(9)組成的高階微分方程組進行動力學計算。對于這種大型復雜非線性動力學高階微分方程組，目前只能采用直接數值積分法，本文采用翟婉明［1］提出的新型快速顯式積分求解此方程組。為了研究不同車速下特別是高速情況下，波浪形磨損和軌道不平順對車輛的振動響應的影響，本文選取的速度分別為 350 km/h、300 km/h、250 km/h、200 km/h。通過計算得到了在激勵模型下，系統的動力學響應特性如輪軌接觸力、輪對加速度、轉向架加速度、車體加速度等，典型的計算結果如輪軌垂向接觸力如圖2所示。

圖2 不同速度下的輪軌力Fig.2 Vertical contact force at different speed

圖2中，V-F表示輪軌垂向力，在軌面波浪型磨損與軌道垂向諧波不平順的激勵下，輪軌垂向力V-F呈現準周期波動(由于赫茲非線性接觸力存在，輪軌力并不是絕對的周期變化)，隨著車速的提高，輪軌力不斷增大。同時可見，輪軌力由這樣幾部分組成:①靜態載荷;②波浪性磨損引起的動載荷;③軌道垂向不平順引起的動載荷。計算和識別它們對輪軌力的貢獻能了解不同波長波磨對輪軌接觸力的影響。而前述的改進的EMD方法具有自適應的特性，能對復雜信號(輪軌力、加速度等)進行分解，系統各個部分的動力學指標均能通過這種方法進行分析。

3 車輛－軌道耦合動力學模型

3.1 輪軌接觸力特性的提取與分析

在前述的激勵模型下，輪軌垂向力及其分解如圖3所示，V-F表示輪軌垂向力，V-F c1～c3表示其本征函數，V-F r3表示剩余分量，即靜態輪軌力。車速為350km/h。改進的EMD方法將垂向輪軌接觸力完全分解:V-F c1的幅值為28.25 kN;V-F c2的幅值為12.22 kN;V-F c3的幅值為4.96 kN，即高頻力成分占到動態輪軌力的62%，中頻輪軌力占到動態輪軌力的27%，低頻輪軌力僅僅占動態輪軌力的1%。同理，計算得到速度為 300 km/h、250 km/h、200 km/h情況下動態輪軌力的成分組成:它們的高頻成分分別占到動態輪軌力的50%、43%、46%;它們的中頻成分分別占到動態輪軌力的39%、48%、43%。在四種不同車速下，波長分別為L1=300 mm、L2=1 500 mm的鋼軌面波浪形磨損引起的輪軌力分別占到動態輪軌力的89%、89%、91%、88%，是動態輪軌力的主要成分。

圖3 垂向輪軌力分解(350 km/h)Fig.3 Decomposition of vertical contact force

圖4為不同車速下垂向輪軌力本征函數的比較結果。對于高頻輪軌力成分V-F c1，車速越高，幅值越大，從350 km/h ～200 km/h，其數值分別為 28.25 kN、18.97 kN、13.64 kN、9.02 kN。對于中頻輪軌力成分 V-F c2，其值不隨車速線性變化，在車速350 km/h時，其幅值為12.22 kN;車速為300 km/h時，其幅值為14.93 kN;車速為 250 km/h 時，其幅值為 14.85 kN;車速為200 km/h時，其幅值為8.40 kN。對于低頻力，其幅值隨著速度降低而降低，從350 km/h～200 km/h，其數值分別為 4.96 kN、3.84 kN、2.59 kN、2.2 kN。

圖4 垂向輪軌力本征函數(IMF)Fig.4 Intrinsic mode functions of vertical contact force

由此可見，波磨引起的輪軌力是動態輪軌力的主要成分，速度越高，高頻動態輪軌力所占比例越大，中、低頻輪軌力所占比例越小。

3.2 轉向架動態特性的提取與分析

在前述的激勵模型下，不同車速下轉向架的垂向加速度響應本征函數的比較及其分解結果如圖5所示。

圖5 轉向架加速度響應本征函數(IMF)Fig.5 Intrinsic mode functions of bogie acceleration response

Acc c1～c3為其加速度響應的本征函數，Acc r3為其剩余分量。對于其低頻加速度響應成分Acc c3，車速越高，幅值越大，從200 km/h～350 km/h，其數值分別為0.46 m·s-2、0.67 m·s-2、0.85 m·s-2、1.04 m·s-2。中頻加速度響應成分Acc c2隨速度非線性變化，最大幅值為1.46 m·s-2，出現在250 km/h;其次幅值為1.31 m·s-2，出現在300 km/h;再次幅值為0.99 m·s-2，出現在350 km/h;最小的幅值為0.91 m·s-2，出現在200 km/h。高頻加速度響應成分Acc c1也隨速度非線性變化，最大幅值為0.53 m·s-2，出現在300 km/h;最小的幅值為0.39 m·s-2，出現在200 km/h。

由此可見，轉向架加速度響應的低頻與中頻成分幅值大于高頻成分，速度越高，低頻加速度響應幅值越大，中、高頻加速度響應呈非線性變化。

3.3 車體動態特性的提取與分析

在前述的激勵模型下，不同車速下車體垂向加速度響應本征函數的比較結果及其分解如圖6所示，Acc c1～c3為其加速度響應的本征函數，Acc r3為其剩余分量。高頻加速度響應成分Acc c1隨速度非線性變化，最大幅值為0.022 m·s-2，出現在250 km/h;其次幅值為0.019 m·s-2，出現在300 km/h;再次幅值為0.18 m·s-2，出現在 200km/h;最小的幅值為0.016 m·s-2，出現在350 km/h。對于其低頻加速度響應成分也隨速度非線性變化，當車速為250 km/h，其Acc c2與之對應且其幅值最大為0.026 m·s-2。其它車速時Acc c3與之對應:在車速為200 km/h，其幅值為0.015 5 m·s-2;在車速為 300 km/h，其幅值為0.014 6 m·s-2;在車速為 350 km/h，幅值為 0.015 m·s-2。

圖6 車體加速度響應本征函數(IMF)Fig.6 Intrinsic mode functions of car body acceleration response

由此可見，由于一系懸掛和二系懸掛的隔離，車體的加速度響應及其分量隨速度非線性變化，且其幅值變化的范圍很小，小于乘車舒適度標準0.13 g，因此車體具有很好的舒適性。對于不同車體，也可以通過此方法研究其設計參數的合理性、乘坐舒適度。

4 結論

本文對經典的EMD方法進行了適當的改進:①以改進的極值域均值模式分解方法(IEMMD)提高局部均值的求解精度;② 以邊界波形匹配預測法抑制端點效應。同時運用車輛-軌道耦合動力學，計算了系統在波浪形磨損和垂向諧波不平順組合情況下的各種動力學響應(輪軌力、輪對和車體的加速度響應)。然后運用改進的EMD方法對系統的動力學響應進行了經驗模式分解，詳細分析了各個響應分量及其變化規律。計算結果表明這種改進的EMD分解方法自適應地將信號分解成本征模態，能有效地提取車輛-軌道系統的動力學特性。

［1］翟婉明.車輛軌道耦合動力學［M］.北京:中國鐵道出社，1997:20-60.

［2］Zarembski A M.Type of rail corrugation［J］.Rail Gazette International.1989，85(8):12-13.

［3］ Grassie S L ，Kalousek J.Rail corrugation:characteristics，cases and treatments［J］.Proceeding of the Institution of Mechanical Engeersrs，1993，207:57-68.

［4］張繼業，金學松，張衛華.高頻力輪軌相互作用下鋼軌的波磨［J］.摩擦學學報，2003，23(2):128-131.

［5］張 潔，陳春俊，林建輝.基于小波變換的軌道不平順信號分析［J］.西南交通大學學報，2004，39(4):469-471.

［6］余 泊.自適應時頻分析方法及其在故障診斷中的應用［D］.大連:大連理工大學，1998.

［7］蓋 強.局域波時頻分析方法的理論研究與應用［D］.大連:大連理工大學，2001.

［8］ Zhao J P，et al.Mirror extending and circular spline function for empirical mode decomposition method［J］.Journal of Zhejiang University，2001，2(3):247-252.

［9］鄧擁軍.EMD方法及Hilbert變換中便捷問題的處理［J］.科學通報，2001，46(3):257-263.

［10］邵晨曦，王 劍，范金鋒，等.一種自適應的EMD端點延拓方法［J］.電子學報，2007，35(10):1944-1946.

［11］徐志勝，翟婉明，王開云.基于Timoshenko梁模型的車輛軌道耦合振動分析［J］.西南交通大學學報，2003，38(1):22-27.

［12］雷曉燕，張 斌，劉慶杰.列車-軌道系統豎向動力分析的車輛軌道單元模型［J］.振動與沖擊，2010，29(3):168-173.

［13］宋 穎，杜彥良，孫寶臣.壓電傳感技術在輪軌力實時監測中的應用探討［J］.振動與沖擊，2010，29(1):228-232.

［14］ Gullers P，Andersson L，Lundén R.High-frequency wheel-rail contact forces-field measurements and influence of corrugation［C］.Proceedings of the Seventh International Conference on Contact Mechanics and Wear of Rail/Wheel Systems，Brisbane，Australia，September 2-27，2008，265，1472-1478.