高速列車乘坐舒適度評價的概率密度演化方法

劉章軍， 張傳勇， 齊東春

(1. 防災減災湖北省重點實驗室(三峽大學)，湖北 宜昌 443002； 2. 三峽大學 土木與建筑學院，湖北 宜昌 443002)

科學技術的不斷進步和社會經濟的快速發展，為我國高速鐵路的發展奠定了堅實基礎。隨著列車運行速度的不斷提升，使得軌道不平順的激勵作用被放大，影響旅客乘坐的平穩性和舒適性。近年來，國內外學者在列車舒適度方面進行了大量研究，取得了一系列成果[1-8]。

軌道不平順作為車輛-軌道耦合系統隨機振動的主要外部激勵源，其獲取方法一般有兩種：線路實測和數值模擬。由于線路實測工作量大、成本高、對設備要求高而受到限制，使得數值模擬被廣泛應用。目前，常用的數值模擬方法有：逆傅里葉變換法[9]、三角級數法[10]、白噪聲濾波法和二次濾波法[10]。為了在概率密度層次上描述軌道不平順隨機過程，本文采用譜表示-隨機函數方法[11-13]，利用隨機函數的約束形式能夠有效地降低隨機過程的隨機度，從而實現用一個基本隨機變量即可表達軌道不平順隨機過程的目的。

從狹義上講，舒適度評價是指列車運行過程中由于車體振動對乘客影響的評價。陳祥等基于因子分析與層次分析方法建立了高速列車乘坐舒適度評價模型，通過確定旅客乘坐需求、生理和心理感受等因素的影響權重進行綜合評價；王海涌等利用粗糙集和模糊分析相結合的方法對城際高速鐵路舒適性進行了分析；國云鵬等基于建筑學振動煩惱率對高速列車的平穩性進行了分析；王海涌等基于ANP和VPRS理論，確定高速列車舒適性綜合評價指標的權重。上述研究為高速列車乘坐舒適度評價提供了不同的方法。在隨機動力系統分析方面，Li等[14-17]提出的概率密度演化理論可以獲取隨機響應的精細概率信息，在線性與非線性多自由度系統的隨機動力反應分析、體系可靠度計算以及振動控制等方面取得了系統性研究進展。在本質上，軌道不平順模擬的譜表示-隨機函數方法與概率密度演化理論具有一致性，兩者的結合可為車輛-軌道耦合系統的隨機動力響應提供高效精確分析方法，實現旅客乘坐舒適度的精細化計算，進而為高速列車乘坐舒適度評價提供新的途徑和方法。

1 譜表示-隨機函數模擬方法

對于零均值的平穩隨機過程X(t)， 假定其單邊的功率譜密度為SX(ω)，ω為圓頻率(rad/s)，則平穩隨機過程模擬的譜表示為

(1)

式中： Δω為頻率離散步長；N為頻率截斷項數；ω1和ωN分別為下截止頻率和上截止頻率，存在如下的關系

(2)

ωk=ω1+(k-1)Δω

(3)

在式(1)中， {Xk,Yk}(k=1,2,…,N)為一組標準的正交隨機變量，滿足如下的基本條件

E[Xk]=E[Yk]=0,E[XkYk]=0

(4a)

E[XkXl]=E[YkYl]=δkl

(4b)

式中：δkl為Kronecker-delta符號，E[·]為數學期望。

在式(1)中，如果令

(5)

式中：φk(k=1,2,…,N)為在(0,2π)上均勻分布且相互獨立的隨機變量。顯然，式(5)滿足式(4)的基本條件。

將式(5)代入式(1)中，即可得到經典的譜表示模擬方法[18]

(6)

一般地，式(1)稱為基于正交隨機變量的譜表示，式(6)稱為基于隨機相位角的譜表示。由于基于隨機相位角的譜表示可以采用FFT算法，因此，在實際模擬中得到廣泛地應用。

從上述推導過程中，可以發現基于隨機相位角的譜表示是基于正交隨機變量譜表示的一種特例，即通過在基于正交隨機變量的譜表示中定義式(5)得到基于隨機相位角的譜表示。顯然，只要能夠滿足基本條件式(4)，還可以定義標準正交隨機變量Xk和Yk的不同隨機函數形式，從而得到不同類型的譜表示模擬公式。

事實上，式(5)可以理解為是對標準正交隨機變量Xk和Yk施加的一種約束形式，從而使得式(1)中隨機變量的數量2N減少為式(6)中的N。可見，通過施加某種約束條件可以有效地降低隨機過程模擬的隨機度(隨機變量的數量)。

基于這一認識，如果對標準正交隨機變量施加更強的約束條件，則可能極大地降低譜表示模擬公式中的隨機變量數量。為此，將Xk和Yk(k=1,2,…,N)表達為概率分布已知的基本隨機變量Θ的正交函數形式，即定義Xk=gk(Θ)，Yk=hk(Θ)， 其中gk(·)和hk(·)為確定性的正交函數式。在本文中，標準正交隨機變量的隨機函數定義為

(7a)

(7b)

顯然，式(7)定義的隨機函數形式使得基于正交隨機變量的譜表示模擬公式(1)中的隨機變量數量2N減少為1。從而克服Monte Carlo模擬方法需要成千上萬個隨機變量的局限性，為復雜工程結構的隨機動力反應分析提供基礎。

2 軌道不平順模擬

功率譜密度是表征平穩隨機過程的軌道不平順最常用統計函數。工程中常用空間頻率fm(cycle/m)或ωm(rad/m)來表達軌道不平順的功率譜密度。在軌道不平順模擬時，則通常采用時間頻率ft(cycle/s)或ωt(rad/s)來表達功率譜密度。它們之間的關系如下：

ωm=2πfm,ωt=2πft

(8a)

ft=vfm,ωt=vωm

(8b)

式中：v為列車運行速度(m/s)。順便指出，時間頻率的單位也可以表示為

cycle/s=Hz, rad/s=1/s=s-1

(9)

本文選用我國高速試驗線譜作為軌道不平順功率譜密度(七參數譜)，其表達式為[19]

(10)

式中：Sm(fm)為軌道不平順功率譜密度， 其單位為mm2/(cycle/m)；fm為空間頻率(cycle/m)；A～G為軌道不平順譜的特征參數，可根據軌道不平順類型選取不同的參數值。

在利用上節的譜表示-隨機函數方法模擬軌道不平順時，需要將式(10)的空間域功率譜密度轉換為時域功率譜密度。根據能量守恒原理，存在如下的關系

Sm(fm)dfm=St(ft)dft

(11a)

St(ωt)dωt=St(ft)dft

(11b)

于是，根據式(8)、式(10)及式(11)，得到軌道不平順的時域功率譜密度

(12)

或

(13)

式中：St(ft)和St(ωt)都是時域功率譜密度，單位分別為mm2/Hz和mm2/(rad/s)；ft和ωt都是時間頻率，單位分別為Hz和rad/s。

本文僅以軌道高低不平順模擬為例，對于式(10)中的軌道不平順功率譜密度，其左高低不平順特征參數A～G的取值見文獻[19]。考慮列車運行速度v=200 km/h、250 km/h、300 km/h三種工況，即v=55.56 m/s、69.44 m/s、83.33 m/s。由于在計算列車車體的有效加速度時，加權函數的頻率范圍為0.5～80 Hz，因此，空間波長取0.5～200 m。 于是，可計算下截止頻率ω1=2πv/200， 上截止頻率ωN=2πv/0.5，均能夠包含加權函數的頻率范圍。

在模擬軌道不平順隨機過程時，首先需將基本隨機變量Θ在區間(-π,π)上均勻離散化，本文建議的離散代表點計算公式為：

θi=-π+(i-0.5)π/72,i=1,2,3,…,144

(14)

式中：θi為離散代表點值， 代表點的數量nsel=144。

圖1為譜表示-隨機函數方法模擬得到的時域代表性樣本曲線，從圖中可知，生成的代表性樣本具有軌道不平順的平穩性特征。圖2和圖3分別為144條代表性樣本集合的平均功率譜與目標功率譜在時域和空間域上的比較。可見，代表性樣本集合的平均功率譜與目標功率譜具有很高的一致性，表明了本文方法模擬軌道不平順的有效性。

圖1 軌道不平順代表性樣本Fig.1 Representative sample of track irregularity

圖2 樣本功率譜與目標譜在時域的對比Fig.2 Comparison of power spectrum and target spectrum in time domain

圖3 樣本功率譜與目標譜在空間域的對比Fig.3 Comparison of power spectrum and target spectrum in space domain

3 車輛-軌道耦合系統的概率密度演化方法

旅客乘車一般可分為站姿、坐姿和臥姿三種方式，不同乘車方式對于列車各個方向的振動敏感性不同，因而在研究乘坐舒適度時所考慮的側重點也不同。由于旅客在乘坐高速列車時，主要以坐姿為主，坐姿的旅客對列車豎向振動加速度敏感性更為明顯，軌道高低不平順所引起的豎向加速度過大時，會引起更大的噪聲危害，對旅客造成不舒適感。

為此，本文僅以軌道高低不平順作為隨機激勵源，在利用譜表示-隨機函數方法表達隨機激勵源時，其隨機性完全來源于基本隨機變量Θ， 且Θ為在區間(-π,π)上均勻分布。車輛-軌道耦合系統采用豎向振動模型[20]，耦合系統響應的求解，可利用有限元方法來實現。通過將耦合系統分解為上部車輛子系統和下部軌道子系統，采用Newmark-β數值積分方法對兩子系統分別交叉迭代獨立求解[21]。

對于車輛子系統，采用附有二系彈簧阻尼的整車模型，其動力方程為

Muau+Cuvu+Kuxu=Qug+Ful(Θ)

(15)

式中：Mu、Cu和Ku分別為車輛子系統的質量矩陣、阻尼矩陣和剛度矩陣；xu、vu和au分別為車輛子系統的位移、速度和加速度向量， 其中加速度向量au={ac,αc,at1,at2,α1,α2,aw1,aw2,aw3,aw4}T，ac和αc分別為車體的豎向加速度和角加速度，ati和αi(i=1,2)分別為轉向架的豎向加速度和角加速度，awi(i=1,2,3,4)為車輪豎向加速度；Qug為車輛重力向量；Ful(Θ)為輪軌接觸力向量。

對于軌道子系統，其動力方程為

Mlal+Clvl+Klxl=Qlg-Ful(Θ)

(16)

式中：Ml、Cl和Kl分別為軌道子系統的質量矩陣、阻尼矩陣和剛度矩陣；xl、vl和al分別為軌道子系統的位移、速度和加速度向量；Qlg為軌道結構重力向量。

在式(15)和式(16)中，車輛單元采用CRH3型高速列車，軌道單元分別采用我國干線60 kg/m有砟無縫線路和CRTS-Ⅱ型無砟軌道線路，其參數見文獻[21]。

對于上述保守的車輛-軌道耦合隨機動力系統，其解答是存在且唯一的，而且是基本隨機變量Θ和時間t的函數。對于感興趣的具體物理量，如車體豎向加速度ac(隨機過程)，其解答可寫為

ac=Z(Θ,t)

(17)

則車體豎向加速度ac的時間變化率為

(18)

由于耦合動力系統的隨機性完全來源于基本隨機變量Θ，因此，擴展的隨機動力系統(Z,Θ)是一個概率保守系統。根據概率守恒原理的隨機事件描述，可得廣義概率密度演化方程為

(19)

于是，在給定的初始條件和邊界條件下求解式(19)獲得聯合概率密度函數pZΘ(z,θ,t)之后， 即可得到車體豎向加速度Z的概率密度函數

(20)

式中：ΩΘ為Θ的分布空間。

在車輛-軌道耦合系統的概率密度演化數值分析中，基本步驟如下：

(1) 生成軌道不平順隨機過程的代表性樣本集合。基本隨機變量Θ在區間(-π,π)上進行均勻離散化，選取代表性點集{θi,i=1,2,…,144}， 并計算代表性點的賦得概率Pi。 利用譜表示-隨機函數方法模擬生成軌道高低不平順的代表性樣本集合{η(θi,t),i=1,2,…,144}， 同時得到代表性樣本的賦得概率Pi。

(2) 求解確定性的車輛-軌道耦合動力系統。對于給定的軌道高低不平順代表性樣本η(θi,t)，進行確定性的車輛子系統式(15)和軌道子系統式(16)的動力響應計算，其中輪軌接觸力向量Ful(θi)是通過Hertz非線性公式來計算[21]

(21)

式中：G為接觸撓度系數；ywk為第k個車輪位移；ylck和ηk分別為第k個輪軌接觸處的鋼軌位移和軌道不平順值。

在Newmark-β數值積分中， 時間步長取為0.000 2 s，動力響應計算時長取為10 s，收斂精度取為1.0×10-7。

(3) 求解廣義概率密度演化方程。采用有限差分方法求解廣義概率密度演化方程式(19)，得到聯合概率密度函數的數值解答。

(4) 累計求和。對式(20)累計求和得到車體豎向加速度的概率密度函數數值解答。

圖4～圖7為有砟軌道情況下列車運行速度為200 km/h時車體豎向加速度的概率信息。其中，圖4為車體豎向加速度的均值與標準差；圖5為典型時刻的概率密度函數；圖6為概率密度演化曲面；圖7為等概率密度線。由此可以看出，概率密度演化方法可對系統響應進行精細化分析得到豐富的概率信息，這為基于概率密度函數的、精細化的車輛-軌道隨機動力響應分析奠定了基礎。

圖4 均值與標準差Fig.4 The mean and the standard deviation

圖5 典型時刻的概率密度函數Fig.5 Typical probability density function at certain time of instants

圖6 概率密度演化曲面Fig.6 The probability density function surface

圖7 等概率密度線Fig.7 The contour of the probability density function surface

4 舒適度評價的概率密度演化方法

4.1 評價標準

旅客乘坐舒適性是評價車輛運行品質的重要標準，一般分為舒適度評價和平穩性指標評價。其中，舒適度評價參考標準有ISO2631國際標準和UIC513標準；平穩性指標評價參考標準有Sperling標準和GB 5599—85標準。

本文采用UIC513標準的簡化方法進行舒適度評價，其計算公式為[22]：

(22)

豎向有效加速度值可采用下式進行計算

(23)

(24)

表1 UIC513舒適度等級劃分

4.2 舒適度評價

應用概率密度演化方法，可以方便地計算有效加速度(隨機變量)的概率密度函數和分布函數。當列車運行速度分別為200 km/h、250 km/h和300 km/h時，采用概率密度演化方法得到有效加速度的概率密度函數和分布函數，如圖8和圖9所示。對于有砟軌道線路，列車運行速度分別為200 km/h和250 km/h時，車體有效加速度的95%置信點值分別為0.226 m/s2和0.278 m/s2；對于無砟軌道線路，列車運行速度分別為200 km/h、250 km/h和300 km/h時，車體有效加速度的95%置信點值分別為0.223 m/s2、0.267 m/s2和0.307 m/s2。利用式(22)即可計算相應的舒適度值，對于有砟軌道線路，舒適度值分別為1.356和1.668，均屬于2級舒適度；對于無砟軌道線路，舒適度值分別為1.338、1.602和1.842，也均屬2級舒適度。這與文獻[22]采用實測軌檢數據得到的高速列車舒適度評價結果相一致。此外，CRTS-Ⅱ型無砟軌道線路的乘坐舒適度略優于我國干線60 kg/m有砟軌道線路，且隨列車運行速度的提高，無砟軌道線路的優勢逐漸增大。

圖8 有效加速度的概率密度函數Fig.8 Probability density function of effective acceleration

圖9 有效加速度的分布函數Fig.9 Distribution function of effective acceleration

為進一步說明概率密度演化方法與直方圖統計方法的差異，圖10給出了有砟軌道線路下列車運行速度為200 km/h時，采用概率密度演化方法和直方圖統計方法得到有效加速度的分布函數對比圖。

從圖10可知，直方圖統計方法的分級間隔大小對有效加速度的分布函數影響較大，進而影響舒適度評價結果。采用概率密度演化方法則可得到精確的有效加速度分布函數，表明了本文方法的優越性。

圖10 概率密度演化方法與統計方法對比Fig.10 Comparison of probability density evolution method and statistical method

5 結 論

本文采用平穩隨機過程的譜表示-隨機函數方法模擬了我國干線鐵路軌道不平順，并結合概率密度演化方法和UIC513舒適度評價標準，研究了高速列車在軌道不平順隨機激勵作用下的車體豎向加速度和旅客乘坐舒適度。根據分析結果得出如下結論：

(1) 通過本文建議的譜表示-隨機函數方法，可構造強約束的標準正交隨機變量形式，只需一個基本隨機變量即可精確模擬軌道不平順隨機過程，且模擬生成的代表性樣本集合具有完備的概率集，這一特點使得其與概率密度演化方法具有天然一致性，可以方便地對隨機振動響應進行分析。

(2) 概率密度演化方法可以獲得任意時刻所需系統響應的瞬時概率密度函數及其演化過程，這為基于概率密度函數的、精細化的列車隨機動力響應分析奠定了基礎。

(3) 通過概率密度演化方法和直方圖統計方法的比較，可見，概率密度演化方法對有效加速度進行處理可以獲得精確的分布函數，提高了舒適度評價的精確性。