基于狀態(tài)空間法的軌道不平順與車(chē)體橫向加速度關(guān)聯(lián)模型

牛留斌,劉金朝,曲建軍,尹 峰

(中國(guó)鐵道科學(xué)研究院集團(tuán)有限公司 基礎(chǔ)設(shè)施檢測(cè)研究所，北京 100081)

軌道不平順是引起車(chē)輛振動(dòng)的主要來(lái)源，軌道的平順性好壞直接影響車(chē)輛高速運(yùn)行平穩(wěn)性和舒適性。車(chē)體振動(dòng)加速度是衡量車(chē)體振動(dòng)程度的量度，是軌道不平順輸入特征及車(chē)輛動(dòng)力學(xué)響應(yīng)特性的綜合反映。兩者之間關(guān)聯(lián)的研究有助于深入理解車(chē)輛的傳遞特性、找出車(chē)輛敏感波長(zhǎng)、避免車(chē)輛振動(dòng)加劇。同時(shí)，結(jié)合車(chē)輛振動(dòng)和軌道不平順綜合評(píng)價(jià)軌道質(zhì)量狀態(tài)，有利于全面分析軌道病害，科學(xué)指導(dǎo)工務(wù)養(yǎng)護(hù)維修。因此，國(guó)內(nèi)外對(duì)軌道不平順與車(chē)體振動(dòng)之間關(guān)聯(lián)展開(kāi)了大量的研究。文獻(xiàn)[1-2]采用傳遞函數(shù)或短時(shí)傅里葉變換小波分析技術(shù)研究軌道不平順與車(chē)輛振動(dòng)響應(yīng)信號(hào)之間關(guān)聯(lián)；文獻(xiàn)[3-4]通過(guò)車(chē)輛-軌道耦合動(dòng)力學(xué)理論及頻率分析方法研究軌道不平順與車(chē)體振動(dòng)之間的傳遞特性。以上研究多是采用信號(hào)處理方法(如相干分析等)或車(chē)輛-軌道動(dòng)力學(xué)理論借助頻譜分析研究軌道不平順與車(chē)體振動(dòng)加速度之間的關(guān)聯(lián)，側(cè)重于找出引起車(chē)體強(qiáng)烈振動(dòng)的單項(xiàng)軌道不平順敏感波長(zhǎng)定量分布，不能給出軌道不平順與車(chē)體振動(dòng)之間的顯式傳遞特性。車(chē)體振動(dòng)是由軌道不平順引起的綜合響應(yīng)，在近似的情況下，單項(xiàng)軌道不平順引起的車(chē)體振動(dòng)起到主要因素，如文獻(xiàn)[5]利用系統(tǒng)辨識(shí)方法建立ARX模型研究軌道長(zhǎng)波高低不平順與車(chē)體垂向加速度之間傳遞關(guān)系。但相對(duì)于軌道長(zhǎng)波高低不平順和車(chē)體垂向加速度之間近似線性的傳遞關(guān)系，軌道不平順與經(jīng)過(guò)車(chē)輛一、二系懸掛傳遞后的車(chē)體橫向加速度之間存在著較為復(fù)雜的非線性傳遞特性，單項(xiàng)軌道不平順對(duì)車(chē)體橫向加速度的影響隨著頻段而不同，這使得兩者的傳遞關(guān)系更為復(fù)雜，文獻(xiàn)[5]中的ARX模型不能有效地表征兩者之間的傳遞關(guān)系。

本文基于系統(tǒng)辨識(shí)理論，以我國(guó)某列高速綜合檢測(cè)列車(chē)在某高速鐵路上的實(shí)測(cè)軌道不平順及車(chē)體橫向加速度數(shù)據(jù)為基礎(chǔ)，采用多輸入單輸出狀態(tài)空間(State-space)模型建立軌道不平順與車(chē)體橫向加速度之間的關(guān)聯(lián)模型研究?jī)烧咧g的傳遞關(guān)系。利用狀態(tài)空間模型的傳遞函數(shù)及實(shí)測(cè)軌道不平順對(duì)車(chē)體橫向加速度進(jìn)行預(yù)測(cè)，并與文獻(xiàn)[5]中ARX模型預(yù)測(cè)車(chē)體加速度結(jié)果進(jìn)行了對(duì)比。

1 關(guān)聯(lián)模型構(gòu)建

本文構(gòu)建的軌道不平順與車(chē)體橫向加速度之間的狀態(tài)空間模型是系統(tǒng)辨識(shí)理論中的方法之一，是在不使用系統(tǒng)結(jié)構(gòu)本身固有參數(shù)和內(nèi)部復(fù)雜傳遞關(guān)系的基礎(chǔ)上，通過(guò)實(shí)測(cè)得到的輸入輸出數(shù)據(jù)來(lái)辨識(shí)系統(tǒng)的傳遞特性。首先利用相干分析法找出不同頻段上軌道不平順與車(chē)體橫向加速度之間的相干特性，為不同頻段的狀態(tài)空間模型選擇合適的輸入輸出波長(zhǎng)；再根據(jù)高速綜合檢測(cè)列車(chē)實(shí)測(cè)數(shù)據(jù)選擇關(guān)聯(lián)模型的合理階數(shù)及相應(yīng)的模型參數(shù)；最后利用殘差及相關(guān)系數(shù)對(duì)關(guān)聯(lián)模型進(jìn)行驗(yàn)證。

1.1 分頻段方法

車(chē)體橫向加速度是由不同類型的軌道不平順引起的，同一類型的軌道不平順的不同波長(zhǎng)成分將影響著車(chē)體橫向加速度的不同頻段[6]。為了簡(jiǎn)化模型輸入變量，在軌道不平順有效測(cè)試波長(zhǎng)范圍內(nèi)，根據(jù)軌道不平順與車(chē)體橫向加速度之間的相干特性找出不同頻段上與車(chē)體橫向加速度相干較強(qiáng)的模型輸入，將狀態(tài)空間傳遞模型分為低頻和中高頻兩個(gè)頻段。相干函數(shù)[7]能夠表征模型輸入與輸出之間的相干特性，可以反映軌道不平順u和車(chē)體橫向加速度y之間的在不同頻段上的依賴程度。相干函數(shù)Cuy為

( 1 )

式中：ω為空間頻率，m-1；Guy(ω)是軌道不平順輸入u和車(chē)體橫向加速度輸出y的互功率譜，mm·m2/s2；Guu(ω)為模型輸入軌道不平順u的自功率譜，mm2·m；Gyy(ω)為模型輸出車(chē)體橫向加速度y的自功率譜，m3/s4。自功率譜及互功率譜的分析方法參考文獻(xiàn)[7]。本文采用平均周期圖譜法(Welch韋爾奇)進(jìn)行功率譜計(jì)算，段內(nèi)傅里葉變換點(diǎn)數(shù)是4 096，加漢明窗(Hamming)，段與段間數(shù)據(jù)樣本重疊50%。

相干函數(shù)的取值范圍在0～1.00之間，相干函數(shù)值越接近1.00，說(shuō)明模型輸入與輸出之間的依賴性越強(qiáng)。通常認(rèn)為相干函數(shù)值大于0.80時(shí)，輸出與輸入之間存在有較強(qiáng)的相干關(guān)系。

1.2 狀態(tài)空間模型

在狀態(tài)空間模型中，軌道不平順輸入u和車(chē)體橫向加速度輸出y之間的關(guān)系表述[8-9]為

( 2 )

式中：u(t)為t時(shí)刻的軌道不平順輸入，mm；y(t)為t時(shí)刻車(chē)體橫向加速度輸出，m/s2；x(t)為狀態(tài)空間模型的狀態(tài)變量；v1(t)和v2(t)為t時(shí)刻模型誤差值；A、B、C、D為狀態(tài)模型待估參數(shù)。在本文中通過(guò)子空間算法得到模型待估參數(shù)，假定空間狀態(tài)模型的階數(shù)n、模型輸入輸出的個(gè)數(shù)分別為m和r，則參數(shù)矩陣A為n×n矩陣、B為n×m矩陣、C為r×n矩陣、D為r×m矩陣。

空間狀態(tài)模型輸入u和輸出y之間的傳遞函數(shù)為

G(ω)=C(ωIn-A)-1B+D

( 3 )

式中：In為n階單位矩陣。

由實(shí)測(cè)數(shù)據(jù)訓(xùn)練模型參數(shù)A、B、C、D，根據(jù)式(3)得到模型傳遞函數(shù)，再由模型傳遞函數(shù)及實(shí)測(cè)軌道不平順數(shù)據(jù)預(yù)測(cè)相應(yīng)的輸出，公式為

( 4 )

1.3 模型參數(shù)優(yōu)化方法

狀態(tài)空間模型參數(shù)A、B、C、D是通過(guò)實(shí)測(cè)軌道不平順數(shù)據(jù)與車(chē)體橫向加速度數(shù)據(jù)訓(xùn)練估計(jì)得到的，選擇合理的模型階數(shù)有利于優(yōu)化模型參數(shù)的數(shù)量，提高系統(tǒng)辨識(shí)精度，減少模型誤差。

定義模型的損失函數(shù)V為

( 5 )

由式(5)可知，損失函數(shù)是待估參數(shù)A、B、C、D的函數(shù)，最優(yōu)化預(yù)測(cè)模型的合理階數(shù)n應(yīng)使模型損失函數(shù)V最小。獲得模型最小損失函數(shù)的判別方法有FPE(Final Prediction Error)準(zhǔn)則或AIC(Akaike Information Criterion)準(zhǔn)則[8]，其表達(dá)式分別為

( 6 )

( 7 )

式中：ΔFPE和ΔAIC分別為FPE準(zhǔn)則和AIC準(zhǔn)則計(jì)算出的最小誤差值；d為模型中待估參數(shù)個(gè)數(shù)；N為估計(jì)模型參數(shù)所使用的實(shí)測(cè)數(shù)據(jù)樣本個(gè)數(shù)；V為由式( 5 )定義的模型損失函數(shù)。

通過(guò)FPE準(zhǔn)則或AIC準(zhǔn)則，可以得到合理階數(shù)的狀態(tài)空間模型及對(duì)應(yīng)的模型參數(shù)。

1.4 模型精度驗(yàn)證

為了確保所構(gòu)建的狀態(tài)空間模型能夠應(yīng)用于實(shí)際預(yù)測(cè)并保持較高的預(yù)測(cè)精度，一般采用模型預(yù)測(cè)與實(shí)測(cè)數(shù)據(jù)之間的殘差及相關(guān)系數(shù)來(lái)檢驗(yàn)狀態(tài)空間模型的精度。

( 8 )

( 9 )

相關(guān)系數(shù)ρ的絕對(duì)值為0～1.00；當(dāng)相關(guān)系數(shù)ρ的絕對(duì)值大于0.80時(shí)，認(rèn)為兩組數(shù)據(jù)之間具有較強(qiáng)的線性相關(guān)性。

2 數(shù)值試驗(yàn)

車(chē)體橫向振動(dòng)不但與軌道不平順有關(guān)，還與車(chē)輛系統(tǒng)的動(dòng)力性能、車(chē)輪踏面圓順度等因素有關(guān)，所建立模型反應(yīng)了軌道不平順與車(chē)體橫向加速度之間的傳遞關(guān)系。由于不同的檢測(cè)車(chē)輛在動(dòng)力性能上的差異，利用不同檢測(cè)車(chē)輛上的實(shí)測(cè)數(shù)據(jù)建立模型時(shí)得到的傳遞函數(shù)會(huì)有差異。

車(chē)輪踏面圓順度會(huì)引起車(chē)輛周期性振動(dòng)，基波的波長(zhǎng)與車(chē)輪輪徑有關(guān)，不大于3 m，在建構(gòu)模型時(shí)會(huì)將該波長(zhǎng)成分濾除。車(chē)輪圓順度引起的二階及以上諧波成分在總成分中的比例很少，對(duì)模型的建構(gòu)沒(méi)有影響。

本文中軌向不平順是指實(shí)測(cè)左右兩側(cè)軌向不平順的平均值；由于軌道短波不平順引起激振頻率高，對(duì)車(chē)體振動(dòng)加速度影響不大，因此忽略3 m以下軌道不平順對(duì)車(chē)輛橫向加速度的影響，主要考慮中長(zhǎng)波范圍軌道不平順與車(chē)體振動(dòng)之間的關(guān)系，所以在數(shù)據(jù)使用前首先對(duì)實(shí)測(cè)車(chē)體橫向加速度數(shù)據(jù)進(jìn)行了空間頻率為0.35 m-1的低通濾波，保留了大于3 m的波長(zhǎng)成分，軌道檢測(cè)數(shù)據(jù)的采樣間隔為0.25 m，車(chē)輛檢測(cè)速度為200 km/h。

利用構(gòu)建模型的傳遞函數(shù)及實(shí)測(cè)軌道不平順數(shù)據(jù)對(duì)車(chē)體橫向加速度進(jìn)行預(yù)測(cè)，結(jié)果與實(shí)測(cè)數(shù)據(jù)的相關(guān)性對(duì)比是驗(yàn)證模型辨識(shí)精度的方法。所以將實(shí)測(cè)數(shù)據(jù)分為兩個(gè)相互獨(dú)立的部分：一部分用于辨識(shí)訓(xùn)練關(guān)聯(lián)模型的參數(shù)，樣本數(shù)N為10 000組；另一部分軌道不平順數(shù)據(jù)及關(guān)聯(lián)模型用于預(yù)測(cè)車(chē)體橫向加速度，并與實(shí)測(cè)車(chē)體加速度數(shù)據(jù)進(jìn)行對(duì)比。

2.1 狀態(tài)空間模型分頻段相干函數(shù)

高速綜合檢測(cè)列車(chē)實(shí)測(cè)車(chē)體橫向加速度功率譜曲線和車(chē)體橫向加速度與長(zhǎng)波軌向不平順、曲率、超高、水平等的相干函數(shù)曲線，見(jiàn)圖1。由圖1可知，車(chē)體橫向加速度在空間頻域上的能量分布比較集中，主要分布在空間頻域0.010～0.050 m-1范圍內(nèi)。在低于0.005 m-1的頻段時(shí)，車(chē)體橫向加速度與超高及曲率相干函數(shù)值大于0.80；在0.005～0.008 m-1頻段時(shí)，車(chē)體橫向加速度與水平相干函數(shù)值較大；在0.008～0.35 m-1頻段，車(chē)體橫向加速度與長(zhǎng)波軌向不平順相干函數(shù)值較大。

圖1 高速綜合檢測(cè)列車(chē)車(chē)體橫向加速度功率譜及相干函數(shù)

在不同的頻段，不同類型的軌道不平順與車(chē)體橫向加速度的相干特性是不一樣的。分為不同的頻段有助于建立不同類型軌道不平順與車(chē)體橫向加速度之間的狀態(tài)空間模型。由圖1可知，本文分別在小于0.005 m-1頻段和0.005～0.350 m-1頻段，在低頻區(qū)段，狀態(tài)空間模型的輸入為曲率和超高；中高頻區(qū)段，狀態(tài)空間模型的輸入為長(zhǎng)波軌向和水平。低頻模型是對(duì)實(shí)測(cè)數(shù)據(jù)0.005 m-1低通濾波得到的；中高頻數(shù)據(jù)是對(duì)實(shí)測(cè)數(shù)據(jù)頻率為0.005～0.350 m-1帶通濾波得到的。

2.2 低頻狀態(tài)空間模型

由式(7)AIC準(zhǔn)則公式計(jì)算值與狀態(tài)空間低頻車(chē)體橫向加速度關(guān)聯(lián)模型階數(shù)n之間的曲線見(jiàn)圖2，圖2中當(dāng)模型的階數(shù)n選為6時(shí)，對(duì)應(yīng)的模型AIC準(zhǔn)則值最小。所以，在建立中高頻橫向加速度與軌道不平順之間的狀態(tài)空間關(guān)聯(lián)模型時(shí)，選用的階數(shù)n為6，所建模型的傳遞函數(shù)特性曲線見(jiàn)圖3，幅頻特性曲線變化平緩，相頻特性曲線顯示兩者沒(méi)有相位差。

圖2 狀態(tài)空間模型AIC與階數(shù)關(guān)系圖(低頻)

圖3 狀態(tài)空間法傳遞函數(shù)特性(低頻模型)

利用狀態(tài)空間模型傳遞函數(shù)及實(shí)測(cè)曲率、超高數(shù)據(jù)預(yù)測(cè)的車(chē)體橫向加速度與實(shí)測(cè)加速度的對(duì)比圖見(jiàn)圖4。圖4中顯示了直線及曲線線況條件兩者波形吻合良好，其相關(guān)系數(shù)為0.98，直線部分的對(duì)比見(jiàn)圖5。

圖4 狀態(tài)空間法預(yù)測(cè)低頻車(chē)體橫向加速度及對(duì)比

圖5 狀態(tài)空間法預(yù)測(cè)低頻車(chē)體橫向加速度及對(duì)比(直線區(qū)段)

兩者的殘差re服從正態(tài)分布N(-0.001 3，0.004 6)，近似滿足0均值假設(shè)，其正態(tài)概率分布曲線見(jiàn)圖6。

圖6 狀態(tài)空間模型殘差分布圖

2.3 中高頻狀態(tài)空間模型

中高頻狀態(tài)空間模型的階數(shù)是通過(guò)AIC準(zhǔn)則進(jìn)行優(yōu)化定階的，按式(7)AIC準(zhǔn)則公式計(jì)算值與狀態(tài)空間模型階數(shù)n之間的關(guān)系曲線見(jiàn)圖7，圖7中曲線表明，當(dāng)狀態(tài)空間模型的階數(shù)n選為6時(shí)，對(duì)應(yīng)的模型AIC準(zhǔn)則值最小，模型最合理。所以，在建立中高頻橫向加速度與軌道不平順之間的狀態(tài)空間模型時(shí)，選用的階數(shù)n為6。

在0.005～0.350 m-1頻段利用狀態(tài)空間法建立長(zhǎng)波軌向不平順、水平與車(chē)體橫向加速度之間的多輸入單輸出(MISO)關(guān)聯(lián)模型，由式(2)得到模型輸入與輸出之間的傳遞函數(shù)幅頻特性與相頻曲線見(jiàn)圖8，兩者的相頻特性曲線上在0.047、0.062、0.21 m-1處出現(xiàn)了相位由-π到π的突變，在利用模型傳遞函數(shù)預(yù)測(cè)輸出時(shí)將造成該頻段預(yù)測(cè)數(shù)據(jù)的不準(zhǔn)確。

圖8 狀態(tài)空間法傳遞函數(shù)特性(中高頻)

利用6階狀態(tài)空間模型傳遞函數(shù)及實(shí)測(cè)軌道長(zhǎng)波軌向、水平預(yù)測(cè)到的車(chē)體橫向加速度，并與實(shí)測(cè)車(chē)體橫向加速度進(jìn)行了對(duì)比，見(jiàn)圖9，兩者的相關(guān)系數(shù)為0.85，兩者的殘差服從正態(tài)分布N(0.000 3，0.006 3)，滿足殘差近似0均值假設(shè)，方差0.006 3反映了模型殘差偏離零均值的程度較小，模型輸出值與實(shí)測(cè)值很接近。殘差正態(tài)概率分布曲線見(jiàn)圖10。

圖9 中高頻橫向加速度預(yù)測(cè)結(jié)果圖(圓曲線部分)

圖10 中高頻殘差分布曲線

利用n階狀態(tài)空間模型預(yù)測(cè)車(chē)體橫向加速度時(shí)，前n個(gè)初始值設(shè)置為0，在開(kāi)始預(yù)測(cè)第n+1個(gè)值時(shí)，存在較大的誤差，經(jīng)過(guò)20至40步的迭代計(jì)算，預(yù)測(cè)車(chē)體橫向加速度數(shù)據(jù)將逐漸收斂于實(shí)測(cè)波形。

除了在個(gè)別高頻數(shù)據(jù)外，利用空間狀態(tài)模型預(yù)測(cè)的車(chē)體橫向加速度與實(shí)測(cè)車(chē)體橫向加速度吻合良好，兩者之間的相關(guān)系數(shù)較高，存在著較強(qiáng)的線性相關(guān)性。

傳遞函數(shù)能夠反映狀態(tài)空間模型能夠辨識(shí)軌道不平順與車(chē)體橫向加速度之間關(guān)聯(lián)關(guān)系的準(zhǔn)確性，而傳遞函數(shù)的正確性可以通過(guò)實(shí)測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證。

3 模型驗(yàn)證及對(duì)比

利用實(shí)測(cè)軌道檢測(cè)數(shù)據(jù)及狀態(tài)空間關(guān)聯(lián)模型傳遞函數(shù)預(yù)測(cè)車(chē)體橫向加速度，并與相應(yīng)的實(shí)測(cè)車(chē)體橫向加速度結(jié)果進(jìn)行對(duì)比。利用ARX模型建立軌道不平順與車(chē)體橫向加速度關(guān)聯(lián)模型的方法可參考文獻(xiàn)[5]。圖11、圖12分別利用狀態(tài)空間模型傳遞函數(shù)、ARX模型傳遞函數(shù)并與相應(yīng)的實(shí)測(cè)車(chē)體橫向加速度進(jìn)行對(duì)比。圖11、圖12中所選用的樣本數(shù)據(jù)區(qū)段長(zhǎng)度為1 500 m，曲線半徑為12 000 m，共有6 000個(gè)樣本點(diǎn)。

圖11 狀態(tài)空間模型預(yù)測(cè)車(chē)體橫向加速度數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比圖

圖12 ARX模型預(yù)測(cè)車(chē)體橫向加速度數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比圖

由圖12可知，無(wú)論線型是直線或者曲線，狀態(tài)空間模型及ARX模型預(yù)測(cè)的車(chē)體橫向加速度波形與實(shí)測(cè)波形吻合一致，但狀態(tài)空間模型預(yù)測(cè)的低頻、中高頻及整體數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的相關(guān)系數(shù)分別為0.97、0.86、0.94，而ARX模型的相關(guān)系數(shù)分別為0.86、0.80、0.82。

由上述利用關(guān)聯(lián)模型傳遞函數(shù)預(yù)測(cè)車(chē)體橫向加速度結(jié)果與實(shí)測(cè)車(chē)體橫向加速度數(shù)據(jù)對(duì)比可知狀態(tài)空間模型能夠較準(zhǔn)確的辨識(shí)軌道不平順與車(chē)體橫向加速度之間的傳遞關(guān)系，而ARX模型辨識(shí)結(jié)果的精度有待提高。

4 結(jié)論

在軌道不平順與車(chē)體橫向加速度相干分析基礎(chǔ)上確定狀態(tài)空間模型的輸入，分別在小于0.005 m-1和0.005～0.350 m-1兩個(gè)空間頻段上建立軌道不平順與車(chē)體橫向加速度之間的關(guān)聯(lián)模型，利用高速綜合檢測(cè)列車(chē)實(shí)測(cè)數(shù)據(jù)訓(xùn)練優(yōu)化關(guān)聯(lián)模型得到傳遞模型的合理階數(shù)及結(jié)構(gòu)參數(shù)，模型的殘差近似符合零正態(tài)分布的假設(shè)。

對(duì)比實(shí)測(cè)軌道不平順及模型傳遞函數(shù)預(yù)測(cè)車(chē)體橫向加速度數(shù)據(jù)和實(shí)測(cè)車(chē)體橫向加速度結(jié)果驗(yàn)證所構(gòu)建的狀態(tài)空間模型傳遞函數(shù)的準(zhǔn)確性，分析結(jié)果表明合理階數(shù)的狀態(tài)空間模型能夠較好地辨識(shí)軌道不平順與車(chē)體橫向加速度之間的傳遞關(guān)系。