一種高速列車蛇行分岔圖和頻率的快速算法

姜成英，羅仁

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

0 引言

蛇行運動的穩(wěn)定性是鐵道車輛系統(tǒng)動力學(xué)中最為重要的研究內(nèi)容之一，圍繞該課題，學(xué)者們多從蛇行運動的發(fā)生原因和抑制等方面展開研究。李曉峰等針對動車組在運營期間存在主頻1～3Hz的晃車問題，從車輪踏面外形狀態(tài)、車輪鏇修質(zhì)量、鋼軌廓形打磨精度等方面進行了調(diào)查分析，給出了經(jīng)濟性鏇修、線路適應(yīng)性的改進方案[1]。李凡松等基于實測車輪外形和鋼軌廓形，建立剛?cè)狁詈蟿恿W(xué)模型，再現(xiàn)了動車組車體異常彈性振動現(xiàn)象。分析表明：異常的輪軌匹配關(guān)系導(dǎo)致轉(zhuǎn)向架的蛇行運動頻率與車體一階菱形模態(tài)頻率接近，進而引發(fā)了車體的異常振動；前后轉(zhuǎn)向架的反相位蛇行運動則不會激發(fā)車體的一階菱形模態(tài)[2]。山長雄亮等通過滾動試驗臺研究了輪對不同初始激擾下的蛇行運動發(fā)生條件。研究表明：當車速介于非線性臨界速度和線性臨界速度之間時，亞臨界Hopf分岔的發(fā)生與否取決于輪對初始激擾幅值，車速越高，發(fā)生蛇行運動所需的激擾幅值越小[3-4]。劉宏友等基于根軌跡法、極限環(huán)法、多種標準等12種方法，計算了某型高速貨車轉(zhuǎn)向架在空載工況下的蛇行失穩(wěn)臨界速度。結(jié)果表明：降速法、UIC 518法、TSI L84法的計算結(jié)果最為保守[5]。鄔平波等以時間為慢變參數(shù)，給車輛系統(tǒng)施加小激擾并緩慢提高車速，可以計算得到車輛系統(tǒng)的線性臨界速度；讓車輛系統(tǒng)在足夠高的速度下發(fā)生周期蛇行運動，然后緩慢降低車速并讓車輛在光滑軌道上運行，可以計算得到非線性臨界速度[6]，后者也被稱為降速法。

綜上所述，蛇行運動臨界速度、分岔圖幅值、蛇行頻率均是高速列車的重要參數(shù)，快速精確地計算這些參數(shù)對高速列車優(yōu)化設(shè)計具有重要意義。本文針對降速法計算得到的輪對橫向位移時間歷程曲線存在毛刺、一個周期多個峰值、整體偏移等結(jié)果，提出了采用滑移窗口取峰谷值點的方法計算位移的包絡(luò)線，從而得到蛇行運動分岔圖；并根據(jù)同相位零點迭代搜索計算蛇行頻率。

1 降速法數(shù)據(jù)處理

在車輛高速運行時，蛇行運動具有復(fù)雜的非線性，在不同速度下的蛇行運動可能是單一周期運動、擬周期運動或者混沌運動。大量仿真分析得到的輪對橫移曲線亦存在多種形態(tài)：部分工況在高速下出現(xiàn)不等幅的擬周期振蕩；部分工況出現(xiàn)幅值調(diào)制現(xiàn)象；部分工況存在毛刺和局部凸起。這些曲線加大了幅值包絡(luò)處理的難度。

常規(guī)的包絡(luò)算法會將所有的極值點篩選出來，時域信號局部微小的突起或凹陷容易導(dǎo)致上下包絡(luò)線的異常抖動，局部毛刺也會被判斷為包絡(luò)線，從而得不到正確的蛇行分岔圖和蛇行頻率。圖1(a)為Simpack后處理中通過93號濾波器的極值算法對某工況輪對橫移曲線的處理結(jié)果。圖中可見，該算法將局部突起處的極值誤判為上下包絡(luò)點，導(dǎo)致蛇行分岔圖和蛇行頻率異常波動。圖1(b)為另一工況經(jīng)包絡(luò)處理得到的蛇行運動分岔圖和頻率，可見，曲線存在嚴重的異常波動，蛇行特征并不明顯。

圖1 常規(guī)包絡(luò)算法異常計算結(jié)果

2 蛇行運動幅值的處理

降速法的初始激擾會導(dǎo)致初始時段的輪對橫移極限環(huán)幅值較大或存在劇烈波動。為了消除初始激擾的影響，將初始車速設(shè)置為高于需要分析的車速范圍，然后以所研究車速范圍內(nèi)的穩(wěn)定蛇行運動量為研究對象。由于采樣和計算的誤差，少數(shù)工況的車速中可能存在重復(fù)值，因而，在加載數(shù)據(jù)之后需檢驗車速的單調(diào)性。部分工況的輪對平衡位置偏離了軌道中心線，需要去掉輪對橫移的趨勢項。

本文采用Matlab中的findpeaks函數(shù)計算包絡(luò)線，通過指定合理的參數(shù)閾值，有效地剔除了包絡(luò)線中的異常點。設(shè)置3個參數(shù)：MinPeakDistance、MinPeakProminence和MinPeakHeight，前兩者用于確定包絡(luò)線，后者用于確定最小的包絡(luò)線幅值(低于此幅值認為包絡(luò)線=0)。MinPeakDistance定義了峰值點間的最小間隔n。從最高的峰值點開始計算，并在該峰值點前后各n個點之外尋找次高峰值點，以此類推；通過指定MinPeakDistance可剔除相鄰且峰值較小的點，減小分岔圖中幅值的波動，本文將該參數(shù)設(shè)置為128(與Simpack軟件中的采樣頻率相同)，即大約每隔2s確定1個峰值點。MinPeakProminence和MinPeakHeight則是根據(jù)峰的凸出量、峰點相對零點的高度來篩選峰點。峰的凸出量、高度如圖2中的實線所示?？梢姡ㄟ^指定MinPeakProminence便可以剔除圖1(a)上包絡(luò)線中的小突起。本文將該參數(shù)設(shè)置為0.4mm。通過findpeaks函數(shù)得到包絡(luò)線的峰谷值之后，再選擇峰值>0、谷值<0的點，進行下述去毛刺處理。

圖2 峰值計算方法示意圖

如果輪對橫移曲線中存在毛刺，即使通過設(shè)置最小峰值間隔處理得到的峰谷值中也必會包含這些毛刺。本文通過滑移窗并根據(jù)第一四分位數(shù)Q1、第三四分位數(shù)Q3以及四分位距IQR來進行剔除：窗長設(shè)置為16個點(即在128Hz采樣頻率下的時窗長為0.125s)，滑移的步長與窗長相同；將幅值Q3+IQR的點認定為是離群點，予以剔除。通常最后一幀的長度會小于窗長，本文是將其自后向前補齊至窗長后再行判斷。此外，還應(yīng)去除與上一幀重復(fù)的點。

如果不存在峰谷值以及峰谷值長度不足1幀，該工況未發(fā)生蛇行運動，將最大車速作為臨界速度。否則，根據(jù)蛇行運動幅值=0.2mm來確定臨界速度。

由于上下包絡(luò)線中峰谷值點對應(yīng)的速度不一致，本文將去除毛刺后的峰谷值點均以1m/s為間隔，從6m/s線性插值到166m/s，取上下包絡(luò)的均值得到輪對蛇行運動分岔圖，并將小于臨界速度的極限環(huán)幅值設(shè)為0。蛇行運動幅值處理流程如圖3所示。

圖3 蛇行運動幅值處理的流程圖

3 蛇行運動頻率的處理

在較短的時間間隔下，可以將輪對的蛇行運動看做是周期信號，將相距n個周期同相位點的時間差記為Δt、速度差記為Δv，降速法中的車輛減速度為a，則蛇行運動頻率f為

(1)

在輪對橫移幅值隨時間變化的曲線中，比較易得的同相位點有極值點和零點。由于部分工況極值點附近的幅值波動較大(尤其是發(fā)生擬周期和混沌等復(fù)雜的蛇行運動時)，相鄰的極值點間會存在相位差，從而導(dǎo)致計算的蛇行頻率不準確。因而，本文采用零點來計算蛇行頻率。限于離散采樣，很難在數(shù)據(jù)點中找到幅值為0的點，根據(jù)零點定理，對樣本點進行迭代搜索，找到相鄰異號的兩個點，并將兩點連線與x軸的交點定為零點；再根據(jù)式(1)計算近似蛇行頻率，并將該頻率記做Δv/2處的頻率。為了減弱毛刺的干擾，采用和幅值異常檢測相同的策略來剔除頻率中的毛刺；采用幅值處理中的方法對車速進行插值，進而得到蛇行運動頻率隨車速的變化曲線。也可以采用短時傅里葉變換對這種周期信號進行處理，但要想得到更精確的頻率信息(即足夠的頻域分辨率)，時窗就要足夠?qū)挘瑥亩鴮?dǎo)致時域內(nèi)的分辨率降低。采用極值法、傅里葉變換、零點法計算得到的蛇行運動頻率對比如圖4所示。圖1(a)中的突起會導(dǎo)致極值法計算的頻率曲線劇烈波動；傅里葉變換受限于時窗寬，頻域的分辨率較低，故曲線呈階梯狀變化；而零點法計算的頻率曲線最為光滑，頻率信息也更為完整。

圖4 不同計算方法下的蛇行運動頻率對比

4 典型工況的計算結(jié)果

建立某時速350 km高速列車拖車非線性動力學(xué)模型，考慮一二系懸掛中的彈簧、阻尼和止擋非線性，輪軌接觸采用Hertz彈簧，輪軌蠕滑力采用Fastsim計算。根據(jù)實際測量的鋼軌型面和車輪型面，組合得到比較惡劣的輪軌匹配狀態(tài)，采用降速法計算車輛系統(tǒng)從170m/s降速到0m/s過程中的蛇行運動響應(yīng)，對一位輪對橫向位移數(shù)據(jù)進行處理來計算蛇行運動分岔圖和蛇行頻率。由于采用了較惡劣的輪軌匹配關(guān)系，車輛的8個車輪采用相同的廓形，左右鋼軌廓形也恒定，某些工況計算的臨界速度很低，這與實際運營狀態(tài)不符，這樣設(shè)置主要是為了研究高速列車的蛇行運動規(guī)律，而不是模擬線路運營狀態(tài)。

經(jīng)過大量動力學(xué)仿真發(fā)現(xiàn)，該車固有的蛇行運動主要包括兩類：一類是車輛在較低的速度下就開始發(fā)生蛇行運動，輪對橫移幅值相對較大，蛇行頻率較低；另一類是車輛在較高的速度下才開始蛇行運動，輪對橫移幅值較小，蛇行頻率較高。前者主要是由于輪軌匹配等效錐度低于正常運用范圍所導(dǎo)致，后者主要是由于車輪磨耗導(dǎo)致輪軌等效錐度過大引起的。前者當?shù)刃уF度接近0時容易出現(xiàn)輪對橫移具有趨勢項的現(xiàn)象；后者的輪對橫移曲線容易在毛刺、峰值點附近出現(xiàn)局部凹凸現(xiàn)象以及在高速下發(fā)生擬周期和混沌等復(fù)雜蛇行運動。

采用本文方法對大量動力學(xué)仿真結(jié)果進行處理，均能得到比較光滑和準確的蛇行運動近似分岔圖和蛇行運動頻率曲線，下面給出3個典型結(jié)果。

對于低頻大幅晃車的蛇行運動處理結(jié)果如圖5所示。由圖可見，本文算法能夠較好地提取不同車速下的極限環(huán)幅值，得到比較光滑的蛇行運動分岔圖(即圖中的極限環(huán)幅值曲線)和蛇行頻率圖。某些工況會發(fā)生輪對橫移幅值在不同車速下的不等幅振蕩，縮小滑移窗口也可以得到更加細微的蛇行分岔圖。

圖5 低頻大幅晃車的蛇行運動計算結(jié)果

對高頻小幅抖車工況的蛇行運動處理結(jié)果如圖6所示。由圖可見，極限環(huán)幅值曲線很好地貼合了輪對橫移包絡(luò)線，臨界速度的取值也能滿足工程精度要求，蛇行運動頻率曲線比較光滑。

圖6 高頻小幅抖車工況的蛇行運動計算結(jié)果

輪對橫移幅值曲線中存在毛刺的工況處理結(jié)果如圖7所示。由圖可見，極限環(huán)幅值曲線和蛇行頻率曲線比較光滑，基本不受局部毛刺的影響，說明本文的算法可以有效地消除毛刺的影響。

由于在車輛優(yōu)化設(shè)計和性能預(yù)測時，主要關(guān)注蛇行運動臨界速度、蛇行幅值范圍和頻率范圍，故圖5-圖7的計算結(jié)果已經(jīng)足夠滿足工程應(yīng)用需要。

圖7 對橫移曲線具有毛刺工況的處理結(jié)果

5 結(jié)語

1)在計算峰谷值的滑移窗函數(shù)中設(shè)置合適的極值最小橫向間隔和極值最小高度間隔，并采用離群點剔除方法，可以有效消除時間歷程曲線中的局部凹凸和毛刺的影響，得到比較光滑的蛇行運動分岔圖。

2)極值法、短時傅里葉變換計算的頻率曲線容易出現(xiàn)劇烈波動、分辨率較低等問題，而零點法計算的頻率曲線較為光滑且更加準確。

3)在車輛優(yōu)化設(shè)計和性能預(yù)測的仿真分析中，主要關(guān)注蛇行運動臨界速度、幅值和頻率范圍，本文方法的計算結(jié)果能夠滿足工程應(yīng)用需要。