基于ARMAX代理模型的車-橋耦合系統(tǒng)可靠性研究

向活躍，陳緒黎，李永樂

（西南交通大學(xué)土木工程學(xué)院， 四川 成都 610031）

高速鐵路中，橋梁占線路比重通常較大，京滬高鐵橋梁工程占比達(dá)86.5%，橋上行車安全成為最受關(guān)注的問題之一.橋上行車安全常采用耦合振動方法進(jìn)行研究，但受軌道幾何不平順等因素的影響，車-橋耦合系統(tǒng)響應(yīng)具有明顯的隨機(jī)性.已有研究表明，同一列車即使以相同速度在同一橋梁上在行駛，車-橋耦合系統(tǒng)的動力響應(yīng)變異系數(shù)也可多達(dá)30%[1].因此，研究車-橋耦合系統(tǒng)響應(yīng)的可靠性具有重要意義.

車-橋耦合振動的求解方法主要有隨機(jī)振動方法和蒙特卡羅模擬（monte carlo simulation, MCS）法.隨機(jī)振動方法是通過求解系統(tǒng)的統(tǒng)計響應(yīng)以分析系統(tǒng)的可靠性，而MCS法是時域法，通過對大量樣本進(jìn)行統(tǒng)計得到系統(tǒng)響應(yīng)的統(tǒng)計值.

目前，在車-橋耦合系統(tǒng)中應(yīng)用的隨機(jī)振動方法主要有協(xié)方差分析法、虛擬激勵法及概率密度演化法等.晉智斌[2]將軌道不平順激勵轉(zhuǎn)化為白噪聲激勵，提出了車-橋時變系統(tǒng)豎向隨機(jī)振動分析的方差遞推算法.張志超等[3-5]采用虛擬激研究了車-橋耦合系統(tǒng)隨機(jī)振動響應(yīng)，得出了有益的結(jié)論.彭勇波等[6-7]基于概率守恒原理獲取系統(tǒng)響應(yīng)的概率密度信息，從而提出分析車-橋耦合系統(tǒng)隨機(jī)振動的概率密度演化近似方法.隨機(jī)振動方法通常求解效率較高，但通過這種方法得到系統(tǒng)響應(yīng)統(tǒng)計量，并根據(jù)均值和方差進(jìn)行評價，且多數(shù)方法需要對輪軌進(jìn)行線性化處理.

在實(shí)際工程中，常用實(shí)際測量或計算機(jī)模擬得到橋面或軌道不平順，用時域法求解非線性車-橋耦合系統(tǒng)的動力響應(yīng)[8-11].采用直接MCS法進(jìn)行可靠性評價時需要大量樣本，當(dāng)橋梁為大跨度桁梁橋時，計算耗時會進(jìn)一步增大.針對其中關(guān)鍵響應(yīng)指標(biāo)，采用動態(tài)系統(tǒng)代理模型可顯著減少計算時間.陳皓等[12]提出了一種基于貝葉斯算法的NARX （nonlinear autoregressive with exogenous input）模型，能較精確地預(yù)測列車脫軌系數(shù).龐學(xué)苗[13]以一種改進(jìn)的NARX神經(jīng)網(wǎng)絡(luò)，提高輪軌相互作用力的預(yù)測精度.Han等[14]用一種NARX-ANN （帶外生輸入的非線性自回歸人工神經(jīng)網(wǎng)絡(luò)）代理模型分析了軌道高低不平順作用下車-橋耦合系統(tǒng)隨機(jī)振動的豎向響應(yīng).Xiang等[15]用一種子集分裂模擬結(jié)合NARX代理模型的方法分析了小概率條件下車橋耦合豎向系統(tǒng)的隨機(jī)振動.通過代理模型的方法，可顯著減少每個樣本的計算時間，提高計算效率，但采用NARX-ANN訓(xùn)練樣本時所需時間仍然較長，此外所有的時間滯后項(xiàng)數(shù)也較多，單個樣本的計算效率有待進(jìn)一步提高.

ARMAX （autoregression and moving average with exogenous）模型作為多項(xiàng)式線性回歸模型[16]，與傳統(tǒng)的ARMA （autoregression and moving average）相比，既融合了AR （autoregression）模型和MA（moving average）模型的特點(diǎn)，線性組合現(xiàn)在的干擾值、過去的觀測值和干擾值，又引入與輸出序列相關(guān)的另一序列（即輸入序列）作為回歸項(xiàng)，提高模型的準(zhǔn)確性.

本文采用ARMAX代理模型研究三維列車過橋的行車安全性.從車輛模型、橋梁模型、輪軌相互作用力及軌道幾何不平順模擬等方面介紹車-橋耦合振動模型.回顧ARMAX模型的原理，給出基于ARMAX代理模型的車-橋耦合系統(tǒng)隨機(jī)振動響應(yīng)分析的算法步驟.討論ARMAX代理模型預(yù)測列車橫向和豎向加速度響應(yīng)的效率和精度，分析車-橋耦合系統(tǒng)的可靠性.

1 車-橋耦合振動模型

1.1 車輛模型

列車可視為四軸二系懸掛質(zhì)量彈簧阻尼車輛[17]，車體各部件由彈簧元件和阻尼元件連接，包含1個車體、2個轉(zhuǎn)向架和4個輪對共7個剛體部件.車體和轉(zhuǎn)向架有橫移、浮沉、側(cè)滾、點(diǎn)頭及搖頭5個自由度，輪對有橫移和搖頭2個自由度，整節(jié)車輛共23個自由度.基于達(dá)朗貝爾原理可建立車輛運(yùn)動方程，見式（1）.

式中：、和yv分別為車輛的加速度、速度和位移向量；Fbv為列車受到來自橋梁輪軌的相互作用力矩陣；Mv、Cv和Kv分別為車輛質(zhì)量、阻尼和剛度矩陣.

分析中采用德國ICE3高速列車參數(shù)，車速取200 km/h，列車為8車編組.為討論方便，8車均采用動車參數(shù)進(jìn)行分析.

1.2 橋梁模型

橋梁采用跨度為32 m的預(yù)應(yīng)力混凝土箱形截面橋[14]，截面面積為9.031 m2，豎向慣性矩為10.939 m4，橫向慣性矩為85.920 m4.

在商業(yè)軟件中，采用空間梁單元建立橋梁單元，并通過模態(tài)分析得到橋梁結(jié)構(gòu)的自振特性，同時導(dǎo)出橋梁模型的質(zhì)量矩陣和剛度矩陣，結(jié)合瑞利阻尼計算方法，得到橋梁模型的阻尼矩陣.將橋梁模型的質(zhì)量、剛度和阻尼矩陣導(dǎo)入MATLAB中可得到橋梁運(yùn)動方程，見式（2）.

本研究中，為減弱非平穩(wěn)特性的影響，橋梁設(shè)置為50跨，均為簡支梁，共長1 600 m，墩高為10 m，墩梁均為混凝土結(jié)構(gòu)，阻尼比為2%.

1.3 輪軌相互作用

輪軌接觸力包括法向力和蠕滑力，法向力公式中的部分系數(shù)是蠕滑力的函數(shù)，蠕滑力公式中的蠕滑系數(shù)又是法向力的函數(shù)，所以法向力與蠕滑力互為隱函數(shù)，分析中輪軌接觸法向力采用輪對的豎向和側(cè)滾振動方程進(jìn)行迭代求解.

輪軌接觸法向力為

式中：NL、NR分別為左、右輪法向力（L表示左，R表示右，余同）；δL、δR分別為輪軌左、右接觸角；θw為側(cè)滾角位移；系數(shù)A、B為蠕滑力的函數(shù)，見式（4）.

式中：mw為輪對質(zhì)量；為輪對橫移加速度；Fypw為一系懸掛系統(tǒng)對輪對的橫向力；Fzpw為一系懸掛系統(tǒng)對輪對的垂向力；yw、ψw分別為橫移位移和搖頭角位移；TyL、TyR（TxL、TxR）分別為左、右輪的橫向（縱向）蠕滑力.

蠕滑力采用Kalker線性理論進(jìn)行計算，不計自旋蠕滑影響，見式（5）.

式中：Tx、Ty分別為縱向、橫向蠕滑力；Mz為旋轉(zhuǎn)蠕滑力矩；f11、f22、f23為蠕滑系數(shù)，由法向力、輪軌接觸幾何參數(shù)和Kalker系數(shù)等進(jìn)行計算，可參考文獻(xiàn)[18]；ξx、ξy分別為縱向、橫向蠕滑率，見式（6）.

式中：λ為踏面斜度；r0為輪對中心滾動圓半徑；ya為軌道中心線橫移位移；θc為軌道水平不平順位移；b為輪對兩滾動圓間距之半；為輪對搖頭角速度；v為輪對中心速度；為輪對橫移速度；為輪對合成轉(zhuǎn)角速度.

輪軌接觸力的計算參數(shù)會隨軌道不平順變化而變化，具有一定非線性輪軌接觸，而代理模型對這類系統(tǒng)參數(shù)變化問題適用性相對較差，為減弱系統(tǒng)參數(shù)的非線性和時變特性的影響，將式（6）的部分參數(shù)進(jìn)行線性化處理，即：滾動圓半徑r0近似取為0.43；踏面斜度 λ取為0.04；b取為0.75.而該簡化對列車響應(yīng)的極值影響相對較小.

由于未考慮輪對的豎向自由度，并假定輪對和軌道為剛體，且與軌道密貼，因此，豎向的軌道不平順作用會直接傳遞到轉(zhuǎn)向架，可寫成式（7）形式[19].

轉(zhuǎn)向架浮沉振動荷載：

轉(zhuǎn)向架側(cè)滾振動荷載：

轉(zhuǎn)向架搖頭振動荷載：

式中：m= (j?1) × 2 + l，控制輪對求和，當(dāng)j= 1時表示前轉(zhuǎn)向架，i= 1，2，當(dāng)j= 2時表示后轉(zhuǎn)向架，i=3，4；zwi和為輪對中心豎向移位和速度；θwi和 θ˙wi分別為輪對中心側(cè)滾移位和速度；K1y、K1z、K1θ分別為輪對橫向、豎向、側(cè)滾剛度；C1y、C1z、C1θ分別為輪對軸箱懸掛橫向、豎向、側(cè)滾阻尼；l1為輪對軸距之半；h4為轉(zhuǎn)向架構(gòu)架中心到車軸中心線的高度.

在得到輪對的法向力和蠕滑力后，根據(jù)輪軌接觸幾何關(guān)系，可得到車輛和橋梁的相互作用力.

1.4 軌道幾何不平順的模擬

采用AR法模擬軌道幾何不平順，該方法是一種線性自回歸模型，也是全極點(diǎn)模型[20].AR模型作為一種線性濾波器，能近似將白噪聲信號轉(zhuǎn)化為具有指定譜特征的離散隨機(jī)過程.AR模型可用差分方程表示，見式（8）.

式中：W(t)為自回歸序列；p為模型階數(shù)，取值為4；A(k)為自回歸系數(shù)；k為1 ～p中的任意一階模型階數(shù)；t為時間；Δt為時間步長；U(t)為均值取0、方差取σ2的白噪聲序列.

基于維納-辛欽定理可得到軌道不平順功率密度譜的自相關(guān)函數(shù)，見式（9）.

式中：Sw(f)為軌道不平順功率密度譜函數(shù)；f為頻率；Rw(dΔt)為Sw(f)的自相關(guān)函數(shù)；d= 0，1，···，p.

再根據(jù)尤拉-沃克方程建立自相關(guān)函數(shù)與AR模型自回歸系數(shù)A(k)的關(guān)系式，見式（10）

聯(lián)立式（9）和式（10a）求得自回歸系數(shù)A(k)，并代入式（10b）中求得R0；根據(jù)式（11）求得白噪聲序列U(t).

式中：n(t)為滿足標(biāo)準(zhǔn)正態(tài)分布的白噪聲.

生成軌道幾何不平順時，取空間步長為0.2 m，Sw(f)采用德國低干擾譜.

針對式（6）中被簡化的參數(shù)，采用與文獻(xiàn)[21-22]相同的車輛參數(shù)，對比在正弦波軌道不平順[21]和采用德國低干擾譜作為隨機(jī)軌道不平順輸入時（按式（8）進(jìn)行模擬），簡化前、后車輛的加速度響應(yīng)如表1所示.由于文獻(xiàn)[21-22]中車輛模型與本文模型有一定差異（文獻(xiàn)中車輛為35個自由度，并考慮軌道的影響），但考慮在激勵相同時，計算結(jié)果基本吻合（見表1）；雖然在正弦激勵下，簡化后的車體橫向加速度明顯變大，但考慮在隨機(jī)激勵作用時，模型參數(shù)簡化對車體加速度響應(yīng)影響有限，后續(xù)研究均采用隨機(jī)激勵.

表1 加速度響應(yīng)結(jié)果對比Tab.1 Comparison of acceleration responsesm/s2

2 動態(tài)系統(tǒng)的ARMAX代理模型

代理模型類似于一個黑箱，在不了解復(fù)雜系統(tǒng)內(nèi)部詳細(xì)結(jié)構(gòu)的情況下，通過輸入和輸出數(shù)據(jù)來訓(xùn)練代理模型，再次輸入數(shù)據(jù)便可預(yù)測動態(tài)響應(yīng).該模型是針對某一個指標(biāo)，通過在時間軸上迭代來計算響應(yīng)，計算效率會明顯提高.

2.1 ARMAX模型

ARMAX模型是帶外源輸入的自回歸滑動平均模型，也是一種線性輸入-輸出多項(xiàng)式模型，可當(dāng)做一個黑箱模型.單輸入和單輸出的ARMAX 模型結(jié)構(gòu)[24]如圖1所示.圖中：p(t?nk)為輸入信號，nk為輸入的延遲階數(shù)；v(t)為輸出信號；e(t)為均值為0，方差不為0的白噪聲信號；Z(q)、H(q)與G(q)分別為模型輸出、輸入和噪聲項(xiàng)，見式（12）.

圖1 ARMAX模型結(jié)構(gòu)Fig.1 ARMAX model structure

式中：nz為自回歸項(xiàng)階數(shù)；nh為回歸項(xiàng)階數(shù)；ng為噪聲項(xiàng)階數(shù)；[z1z2···znz]、[h1h2···hnh]和[g1g2···gng]為模型參數(shù)；q為后向移位算子，可由式（13）求得.

式中：n取值為nz、nh+ 1、ng.

因此，ARMAX模型可表達(dá)為

式中：Z(q)v(t)為自回歸項(xiàng)，即AR部分；G(q)e(t)為滑動平均項(xiàng)，即MA部分；H(q)p(t?nk)為外源輸入（也可稱為回歸項(xiàng)），即X部分.

聯(lián)立式（12）～ （14），即可得到ARMAX模型的差分方程為

分析中，延遲階數(shù)通過試算法得到.在試算過程中，可觀察系數(shù)A、B、C隨延遲系數(shù)的變化情況，若延遲系數(shù)非常小，則可刪除相應(yīng)的項(xiàng).

2.2 分析過程

對車-橋耦合系統(tǒng)模型而言，列車所受荷載除軌道幾何不平順外，還有橋梁位移.但在分析中，橋梁位移與軌道不平順有一定關(guān)系，考慮計算的效率，針對軌道幾何不平順為0時列車過橋時各個輪對處的橋梁響應(yīng)為輸入，與軌道幾何不平順疊加，形成等效的軌道不平順.這個過程相當(dāng)于僅進(jìn)行了一次全過程迭代，但橋梁的剛度通常較大，一次迭代的橋梁響應(yīng)已與實(shí)際曲線較為吻合.具體分析時，采用以下步驟：

步驟1軌道幾何不平順為0時，提取列車過橋時每個輪對對應(yīng)處的振動位移和加速度.

步驟2任意生成一條隨機(jī)的軌道不平順，進(jìn)行車-橋耦合振動計算，得到車輛的響應(yīng)，將軌道幾何不平順與步驟1中的橋梁響應(yīng)疊加，得到輸入和輸出數(shù)據(jù).

步驟3試算代理模型延遲階數(shù)，預(yù)測響應(yīng)，基于均方根值和最大絕對誤差，判斷預(yù)測準(zhǔn)確性，直到確定代理模型參數(shù)為止.

步驟4隨機(jī)生成軌道幾何不平順，與步驟1中的橋梁響應(yīng)疊加，預(yù)測響應(yīng)，直到所需數(shù)量.尋找響應(yīng)絕對最大值所對應(yīng)的樣本，以相應(yīng)的軌道不平順為輸入，進(jìn)行車-橋耦合振動計算，驗(yàn)證代理模型預(yù)測的準(zhǔn)確性.

步驟5若代理模型預(yù)測準(zhǔn)確性差，則以最大值所對應(yīng)的輸入和輸出重新訓(xùn)練代理模型，重新進(jìn)行步驟4的計算，并將兩次預(yù)測結(jié)果進(jìn)行平均值處理.

步驟6統(tǒng)計所有預(yù)測樣本的信息，形成超越概率曲線，評價行車安全.

不同的車輛響應(yīng)輸入有所差異.對于豎向響應(yīng)，主要由豎向的不平順決定，如車體豎向加速度，主要和式（7）有關(guān).因此，采用代理模型計算豎向加速度時，在生成軌道幾何不平順后，按式（7a）計算豎向的荷載，得到相應(yīng)的輸入用于訓(xùn)練代理模型和預(yù)測響應(yīng).對于橫向加速度，主要受橫向和豎向軌道幾何不平順的影響，但輪軌橫向力（見式（3））不能直接計算得到，因此可直接將4個輪對處的橫向軌道不平順和式（7b）作為輸入.

對于步驟5的處理方式，文獻(xiàn)[14]直接采用最大值樣本的處理方式，但本文中的ARMAX代理默認(rèn)的基函數(shù)為線性函數(shù)，采用最大值時仍有一定誤差，通過對比發(fā)現(xiàn)取平均值的精度會更好.當(dāng)然，若任意產(chǎn)生的樣本即為所模擬樣本集中的最大值，則需要再生成一個響應(yīng)最大值較小的組樣本重新訓(xùn)練代理模型.

分析中，取一條樣本中絕對最大值進(jìn)行統(tǒng)計分析，假如有1000個樣本，最大值樣本是絕對最大值所對應(yīng)的時程曲線.

3 結(jié)果分析

為討論ARMAX代理模型的準(zhǔn)確性和計算效率，針對車-橋耦合振動模型，采用直接蒙特卡羅法計算了1 000個樣本，每個樣本的計算時間約5 min.由計算結(jié)果的統(tǒng)計特征觀察到，在車輛參數(shù)相同，8輛車編組時，各車輛的車體豎向加速度av和橫向加速度ah絕對最大值的均值、1%和99%分位處的響應(yīng)差異均較小（見圖2），因此，后續(xù)主要討論第1輛車的響應(yīng).

圖2 8輛車的車體加速度響應(yīng)Fig.2 Acceleration response of 8 car bodies of the train

圖2的加速度響應(yīng)值較小，這在類似分析算例中也有相似的結(jié)果[25-26].這是因?yàn)榉治鲋杏嬎丬囁龠h(yuǎn)小于其最高運(yùn)營速度，加之橋梁剛度較大，因此車輛響應(yīng)較小.由于軌道不平順對橋梁響應(yīng)影響較小，因此后續(xù)不再討論橋梁響應(yīng).

針對第一輛車的車體豎向加速度進(jìn)行代理模型訓(xùn)練，建立代理模型過程中，通過試算得nz= 4，nh=5，ng= 0，nk= 0，此時訓(xùn)練的代理模型準(zhǔn)確性高.相比NARX等其他代理模型，此時ARMAX代理模型的延遲階數(shù)明顯較少，因此，計算效率會明顯提高，結(jié)果顯示每個樣本的計算時間不超過0.2 s.豎向系統(tǒng)不需進(jìn)行2.2節(jié)中的步驟5便可得到預(yù)測結(jié)果.圖3給出了第1輛車的車體豎向加速度響應(yīng)的理論值和代理模型的預(yù)測值對比.圖中：橫軸為目標(biāo)值，由式（1）和式（2）進(jìn)行直接MCS得到；縱軸為預(yù)測值，由代理模型預(yù)測得到.由圖3可見：ARMAX模型對車-橋耦合豎向系統(tǒng)車輛加速響應(yīng)的預(yù)測具有良好精度和效率.

圖3 豎向加速度的目標(biāo)值和預(yù)測值Fig.3 Target value vs predicted value of vertical acceleration

列車橫向系統(tǒng)為非線性系統(tǒng)，而ARMAX代理模型是線性模型，在模擬橫向系統(tǒng)時會有一定誤差.圖4給出了橫向加速度響應(yīng)的最大絕對值.建立代理模型過程中，通過試算得到na= 50，nb= 50，nc=10，nk= 0.由圖4可知：初始樣本訓(xùn)練得到的代理模型預(yù)測值總體較目標(biāo)值偏大，而最大值樣本所得到的代理模型預(yù)測值較總體目標(biāo)值偏小，取兩者平均值后準(zhǔn)確性明顯提高.

圖4 橫向加速度的目標(biāo)值和預(yù)測值Fig.4 Target value vs predicted value of lateral acceleration

為進(jìn)一步討論代理模型預(yù)測車體橫向加速度響應(yīng)的準(zhǔn)確性，圖5（a）給出了車體橫向加速度目標(biāo)時程和預(yù)測時程的對比曲線，圖5（b）給出了目標(biāo)值和預(yù)測值的最大誤差時程曲線.由圖5（b）可見：最大誤差的幅值小于0.03 m/s2，所以利用ARMAX代理模型預(yù)測非線性的橫向系統(tǒng)響應(yīng)時仍有一定準(zhǔn)確性.同樣，車體橫向加速度每條樣本的計算時間小于0.2 s，代理模型的精度和效率是可接受的.

圖5 車體橫向加速度時程曲線Fig.5 Lateral acceleration time-history curves of car body

鐵路線路有固定的養(yǎng)護(hù)周期，每天運(yùn)行列車的對數(shù)是確定值，因此，評價行車安全的概率水平也應(yīng)確定，這表明從概率角度考慮行車安全可行.為進(jìn)一步研究行車安全性，圖6（a）、（b）給出車體豎向和橫向加速度的超越概率曲線，當(dāng)超越概率為10?3時，列車加速度響應(yīng)仍較小，這在文獻(xiàn)[18]中也有類似的結(jié)果，同時加速度響應(yīng)較小也與計算速度較低有一定關(guān)系，這表明行車安全具有較高的安全儲備.此外，在代理模型基礎(chǔ)上，按MCS統(tǒng)計超越概率時，在保證變差系數(shù)為0.1的條件下，模擬的次數(shù)為超越概率倒數(shù)100倍，由圖6可見：在超越概率為10?2時，預(yù)測值與目標(biāo)值已較為吻合.這表明基于ARMAX代理可近似車橋系統(tǒng)響應(yīng).

圖6 車體加速度超越概率Fig.6 Exceedance probability of car body acceleration

4 結(jié) 論

本文采用MATLAB編寫車-橋耦合系統(tǒng)的分析程序，提出一種ARMAX代理模型近似得到車-橋耦合系統(tǒng)中的車輛豎向和橫向響應(yīng)，并與MCS計算結(jié)果進(jìn)行對比，探究車-橋耦合系統(tǒng)可靠性分析的效率和精度，得出如下結(jié)論：

1）ARMAX代理模型預(yù)測車輛豎向和橫向加速度響應(yīng)的時間短，對車-橋耦合系統(tǒng)進(jìn)行可靠性分析的效率較高；

2）就車-橋耦合系統(tǒng)車輛加速度響應(yīng)而言，豎向響應(yīng)的預(yù)測結(jié)果與MCS結(jié)果吻合較好，橫向響應(yīng)的預(yù)測精度雖然略低于豎向，但平均處理后的超越概率曲線與直接MCS得到的結(jié)果吻合較好.