基于灰箱模型的垂向車-橋系統響應極值分布研究

唐 平，李永樂，向活躍

(西南交通大學 土木工程學院，成都 610031)

隨著我國軌道交通建設的快速發展，高速列車以高速舒適的優勢在全國范圍內得到大規模地建設。在我國高鐵大發展的同時，在營高鐵的可靠性也越來越多地被研究。由于車-橋系統相對復雜的非線性和軌道激勵的隨機性，使得這一研究遇到很多障礙。

隨著隨機振動理論的發展和成熟，學者們提出了很多隨機性方法，并被應用到了車-橋耦合振動中，如協方差分析法[1]、演變隨機響應法[2-3]以及虛擬激勵法[4-7]等，但當對象為非線性系統或需要計算小概率失效事件時，以上方法還有很多缺陷。除了上述隨機振動方法，最直觀通用的還是Monte Carlo模擬(Monte Carlo simulation，MCS)，此模擬方法不受系統線性或非線性影響，但得到的結果變異系數與模擬次數的開方成反比，在失效概率很小時，需要的樣本數目很大。因此，提高每個樣本計算效率或減少樣本數目是非常必要的。在減少非線性結構小失效概率動力可靠性計算所需MCS樣本量方面，比較有效的是子集模擬以及在其框架上衍生的一些方法[8-11]，其本質還是基于條件概率的MCS。

而在提高樣本計算效率方面，有學者提出橋梁代理模型[12]；對于垂向車-橋系統，也有學者采用非線性神經網絡模型[13]；還有學者針對輪軌幾何接觸點的計算[14]、輪軌力的預測[15]提出過一系列人工神經網絡模型。以上方法都將重心放在簡化車-橋模型的某個點上，使用各類代理模型簡化計算。大部分代理模型雖然能在一定程度上提高計算效率，但由于未對計算框架進行簡化，效率提高有限。

本文針對垂向車-橋系統，采用車輛結構的灰箱模型進行參數優化，根據無軌道不平順時近似輪軌力作為車輛灰箱模型輸入，以及考慮輪軌不平順時原始車-橋系統車輛響應作為輸出，對車輛結構灰箱模型參數進行最佳擬合，提出一種高效且高精度的計算框架。由于車輛結構已知，考慮到MATLAB軟件的系統辨識工具箱中灰箱辨識模型能夠方便地建立車輛結構模型，并進行參數優化，為了方便高效，本工作采用了MATLAB軟件中的灰箱模型來進行車輛參數優化。

1 垂向車-橋系統耦合模型

1.1 車輛模型

車輛模型采用CRH型動車，8車編組，每輛車包括一個車體，兩個轉向架，分別考慮其沉浮與點頭自由度，則整個車輛模型為6自由度動力系統，車輛垂向系統模型如圖1所示[16]。

圖1 列車的質量彈簧阻尼模型Fig.1 Mass-spring-damping model of the train

輪軌作用力簡化為直接作用于轉向架的外荷載，此作用力的表達式為：

豎向荷載

(1)

轉動荷載

(2)

式中：j為轉向架號；i為輪對號；Kp，Cp為一系懸掛剛度和阻尼；lw為同一轉向架兩輪對距離之半。

根據達朗貝爾原理，可以列出6個自由度上的運動方程，其具體表達式可以參考毛云霄等的研究，最終組集成矩陣形式表達為

(3)

式中：Mv，Cv，Kv分別為車輛的質量、阻尼、剛度矩陣；Fv為外荷載，其表達式為Fv=[0 0Fbz1Fbφ1Fbz2Fbφ2]T；uv為車輛位移。

1.2 橋梁模型

橋梁采用高鐵中廣泛使用的高架簡支梁橋，橋梁每跨長32 m，共50跨1 600 m。墩高10 m，橋寬12 m，高3.05 m，橋梁截面如圖2所示。

圖2 橋梁橫截面圖(mm)Fig.2 Bridge cross-section (mm)

橋梁采用商業有限元軟件ANSYS建模，并導出橋梁有限元模型的質量與剛度矩陣，再結合橋梁的第一階豎彎與第一階扭轉頻率，采用瑞利阻尼計算公式求得橋梁模型阻尼矩陣。外荷載矩陣需要根據車輛輪對不同時刻位于橋梁各單元的位置以及作用力大小，進行橋梁單元內插值計算，轉換為所在橋梁單元的節點荷載。最終，橋梁模型的運動方程可以采用矩陣形式表達為

(4)

式中：Mb，Cb，Kb分別為橋梁的質量、阻尼、剛度矩陣；Fb為外荷載矢量，其合力與作用在車輛轉向架上的力等大反向。

2 軌道不平順的模擬

由于本工作側重于提出近似的計算方法，而軌道不平順的選取對方法的說明并無影響。本工作中軌道不平順采用德國高速垂向低干擾譜，其表達式為

(5)

式中：Ω為空間頻率，rad/m；Av=4.032×10-7rad·m；Ωc=0.824 6 rad/m；Ωr=0.020 6 rad/m。

軌道不平順的模擬采用白噪聲濾波法，其主要做法是根據需要模擬的功率譜的特點，設計相應的濾波器，將白噪聲作為輸入，輸出與所需模擬的功率譜相吻合的隨機樣本。由于德國低干譜給出的是空間域功率譜函數，需要將其轉化為時間域功率密度譜，式(5)轉化為時間域功率密度譜后表達式為

(6)

根據文獻[17]，可以推導出所需白噪聲濾波器的微分方程為

(7)

(8)

對于式(7)和式(8)所示的線性時不變系統，寫成矩陣形式為

(9)

(10)

由式(9)和式(10)，得到此線性時不變系統的A，B，C，D4個矩陣。即

(11)

此時，只要輸入白噪聲序列r1，即可通過式(9)和式(10)得到軌道不平順時程y。通過這種方法得到的軌道不平順包含了所有的頻率成分，需要再次濾波，去除不需要的低頻和高頻部分。文獻[18]表明，在軌道不平順激勵作用下，車體和輪對的垂向、橫向振動加速度主要頻率范圍為1～3 Hz和3～120 Hz。本工作中取0.4～80.0 m的波長段，并將不屬于此波長段對應頻率的頻率成分濾除。通過濾波后得到的軌道垂向不平順樣本以及樣本功率譜與理論譜的對比，如圖3所示。從圖3可知：通過白噪聲濾波法生成并通過濾波后的軌道不平順樣本功率譜與理論譜吻合較好，驗證了生成的軌道不平順樣本的正確性。

圖3 軌道不平順模擬Fig.3 Simulation of track irregularity

3 車輛灰箱代理模型

3.1 灰箱模型理論

灰箱系統辨識作為系統識別的一種方法，與黑箱系統辨識不同，此辨識方法需要事先了解系統的結構，但是對于其結構參數值不清楚，根據系統輸入和輸出，采用最小化誤差的原則，對結構的參數進行最佳估值[19]。

建立灰箱系統辨識模型時，需要將系統的運動微分方程改寫成狀態方程，并定義狀態初值，狀態方程具有以下形式

(12)

式中：A，B，C，D為由各個系統參數參數化表示的矩陣；K為噪聲矩陣，Ke(t)可以為由系統參數參數化表示的標量；x0為狀態初值向量。

在灰箱辨識模型結構建立之后，以定義的系統參數值初值開始，采用子空間Gauss-Newton最小二乘搜索，Levenberg-Marquardt最小二乘搜索，自適應的子空間Gauss-Newton最小二乘搜索，最速下降法等優化算法[20]，迭代估計系統參數值，以將模型估計的輸出值與需要辨識的輸出值的誤差最小化。

3.2 軌道不平順對橋梁和車輛響應影響

由于高鐵橋梁主梁剛度都比較大，因此，橋梁變形主要來自于車輛重力，軌道不平順導致的輪軌力對橋梁變形影響不大。在此，分別計算了8車編組動車過橋時，有軌道不平順作用及無軌道不平順作用下，第一輛車第一個輪對處的橋梁響應，如圖4所示。從圖4可知：有軌道不平順作用時，橋梁的位移和速度與無軌道不平順作用時大體吻合，只有局部地方有很小的波動。而這些橋梁位移和速度小波動將影響車輛加速度。

圖4 第一輛車第一個輪對處的橋梁響應對比Fig.4 Comparison of bridge displacement and velocity at the front wheel of the first vehicle

為了探究橋梁位移與速度差異對于車輛加速度響應的影響。定義兩種計算方式為方式①和方式②：方式①為直接以無軌道不平順時橋梁位移疊加軌道不平順，再計算輸入車輛系統作用力，最后求解車輛方程的車體響應；方式②為直接解車橋運動方程得到車體響應。分別采用隨機生成的400不平順樣本按方式①和方式②計算。400個樣本中第一輛車車體加速度最大值結果與相同軌道不平順作用下直接解車橋運動方程的結果比較，如圖5所示。

圖5 400個樣本車體加速度最大值比較Fig.5 Comparison of vehicle maximum acceleration of 400 sample

由圖5可知，直接以無軌道不平順時橋梁變形作為有軌道不平順作用下橋梁變形值時(方式①)，計算得到的車體加速值總體偏小。因此，如果能將此誤差均值盡量縮小，就能在計算車輛加速度響應時，直接采用無軌道不平順時橋梁變形值作為有軌道不平順作用下橋梁變形值。這樣做的好處在于，可以在無軌道不平順條件下，計算橋梁變形時程，只需計算一次橋梁變形。這樣，就可以確定每個時步各輪對處橋梁位移，這部分也只需要計算一次。在考慮軌道不平順后，直接在此基礎上疊加軌道不平順。而求解原始車橋運動方程時(方式②)，每個樣本每一時步都要計算橋梁位移值，并且都要進行插值計算各輪對處橋梁位移。因此，如果能直接采用無軌道不平順時橋梁變形值作為有軌道不平順作用下橋梁變形值，當求解多個樣本時，整個計算時間將會大大減少。

3.3 車輛灰箱模型建立流程

對于垂向車橋系統，本應是參數確定系統，但由于要將橋梁位移近似，為了使此近似對車輛響應影響最小，可以細微調整車輛參數。而灰箱系統辨識可以很好實現這個目標。在進行車輛系統灰箱系統辨識時，按方式①計算得到作用在轉向架上的力，以此作為車輛灰箱系統的輸入。再將相同軌道不平順作用下，求解原始車橋運動方程(方式②)得到的車輛響應作為系統輸出值，根據此輸入輸出，對車輛參數進行辨識。

為了灰箱模型建立與辨識的高效可靠，本工作采用了MATLAB軟件中系統辨識工具箱里的灰箱系統辨識組件。建立車輛灰箱系統辨識模型時，由于垂向車橋系統的車輛模型考慮6個自由度，因此，狀態量一共選取了12個，即6個自由度方向的位移與速度，并將式(3)改寫成狀態方程。并將車輛狀態方程編寫為C MEX文件，以導入MATLAB軟件生成車輛灰箱模型。在進行灰箱模型參數優化時，本工作采用的是MATLAB軟件內置的Levenberg-Marquardt最小二乘搜索算法。

在進行灰箱系統識別時，代表樣本的選擇值得一提。根據圖5所示的計算結果，方式①計算得到的車體加速度響應最大值普遍偏小。因此，優化車輛參數的目的是將此種方式計算的響應最大值整體增大。要使所有樣本的最大值誤差均值接近0，則應該選擇中間誤差的樣本作為代表樣本。為此，在選擇代表樣本時，可以先采用100個樣本，按方式①和方式②分別計算車體加速度最大值。并將誤差按大小排列，取中間誤差樣本作為代表樣本。

在得到代表樣本后，按方式①計算作用在車輛轉向架上的力，作為灰箱辨識系統的輸入，按方式②計算得到的車體相應作為輸出，進行系統參數識別。部分參數識別前后變化，如表1所示。由1表可知：車輛參數在辨識前后有了細微調整，從而使得優化后的車輛模型能最小化近似計算的誤差。系統辨識結果如圖6所示。由6圖可知：車體加速度相應的模擬值與原始動力方程求解值吻合很好，這也驗證了建立的灰箱系統辨識模型正確性。

表1 車輛部分參數識別前后比較Tab.1 Comparison of partial vehicle parameters before and after identification

圖6 灰箱辨識模型模擬值與直接求解運動方程結果比較Fig.6 Comparison between simulated values of grey-box identification model and solution of motion equation

在得到車輛灰箱系統后，便可使用其進行垂向車橋系統車輛動力響應計算。

4 計算流程

在車輛的灰箱辨識模型參數優化后，便可利用其計算車輛的加速度響應。計算步驟總結如下：

步驟1在無軌道不平順作用下，計算橋梁位移與速度時程；

步驟2根據步驟1中的橋梁位移與速度時程，插值計算每個時步車輛各輪對處橋梁位移與速度值；

步驟3白噪聲濾波法生成軌道不平順時程，并疊加步驟2中的橋梁位移與速度，得到車輛各輪對處的等效不平順時程，并計算作用在車輛轉向架上的力；

步驟4將步驟3中得到的力時程輸入車輛系統的灰箱模型，得到車輛的響應值。

5 結果分析

按照第4章的計算流程，計算了1 200個樣本(包括圖5中400個樣本的軌道不平順)的第一輛車車體加速度響應最大值。此最大值與根據原始的車橋運動方程計算得到的最大值比較，如圖7所示。由圖7可知：圖5中的400條不平順樣本結果總體偏小的趨勢已被明顯改善，額外的800條隨機樣本下的結果也與原始系統計算值吻合很好，最大誤差不超過5%，并且，誤差均值接近0。在相同計算設備條件(計算機處理器：Intel(R)Core(TM)i5-3210M CPU @2.50 GHz 2.5 GHz，內存：8 GB)下，求解原始車橋運動方程共用時106 670 s，而基于灰箱辨識模型計算方法共用時15 707 s，用時僅為前者的1/7左右，大大提高了求解車輛響應的效率。

圖7 1 200個樣本車體加速度最大值比較Fig.7 Comparison of vehicle maximum acceleration of 1 200 sample

根據此1 200個樣本計算結果，對車體加速度響應極值分布進行了研究，使用本文框架計算的極值分布圖與原始車-橋方程計算的極值分布，如圖8所示。由圖8可知：本文所提方法得到的響應極值分布與求解原始車-橋運動方程得到的極值分布吻合很好，驗證了使用本文所提方法研究垂向車-橋系統響應極值分布的可行性。

圖8 1 200個樣本車體加速度最大值分布比較Fig.8 Comparison of vehicle maximum acceleration distribution of 1 200 sample

6 結 論

本文所提的基于灰箱系統計算車體動力響應的方法，利用橋梁的近似響應時程，結合車輛灰箱辨識模型，簡化了垂向車-橋系統中計算車輛動力響應的過程，并使用1 200個樣本驗證了其高效性與準確性。結果表明：車輛灰箱模型計算誤差較小，相同樣本量的計算用時較少，為求解原始車-橋運動方程的1/7左右，且車體加速度響應極值分布曲線與求解系統原始運動方程結果吻合很好。這對于小概率事件進行模擬時能夠凸顯出其高效的優勢。