向活躍，陳緒黎，李永樂

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

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

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

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

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

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

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

1 車-橋耦合振動模型

1.1 車輛模型

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

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

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

1.2 橋梁模型

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

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

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

1.3 輪軌相互作用

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

輪軌接觸法向力為

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

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

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

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

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

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

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

轉向架浮沉振動荷載：

轉向架側滾振動荷載：

轉向架搖頭振動荷載：

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

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

1.4 軌道幾何不平順的模擬

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

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

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

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

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

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

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

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

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

表1 加速度響應結果對比Tab.1 Comparison of acceleration responsesm/s2

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

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

2.1 ARMAX模型

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

圖1 ARMAX模型結構Fig.1 ARMAX model structure

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

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

因此，ARMAX模型可表達為

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

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

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

2.2 分析過程

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

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

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

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

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

步驟5若代理模型預測準確性差，則以最大值所對應的輸入和輸出重新訓練代理模型，重新進行步驟4的計算，并將兩次預測結果進行平均值處理.

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

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

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

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

3 結果分析

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

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

圖2的加速度響應值較小，這在類似分析算例中也有相似的結果[25-26].這是因為分析中計算車速遠小于其最高運營速度，加之橋梁剛度較大，因此車輛響應較小.由于軌道不平順對橋梁響應影響較小，因此后續(xù)不再討論橋梁響應.

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

圖3 豎向加速度的目標值和預測值Fig.3 Target value vs predicted value of vertical acceleration

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

圖4 橫向加速度的目標值和預測值Fig.4 Target value vs predicted value of lateral acceleration

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

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

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

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

4 結 論

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

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

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