劉付山，曾志平,2，王衛(wèi)東,2

(1.中南大學 土木工程學院，湖南 長沙 410075；2.中南大學 重載鐵路工程結(jié)構(gòu)教育部重點實驗室，湖南 長沙 410075)

軌道隨機不平順是由不同波長、不同幅值、不同相位的不同不平順波疊加而成的隨線路里程變化的復雜隨機過程[1]。軌道隨機不平順常通過線路實地測量獲得，并采用功率譜密度函數(shù)表示其統(tǒng)計特征，我國及英國、日本、德國等都測定了各自的軌道不平順譜密度和相關(guān)函數(shù)。然而，在非線性車輛-軌道耦合系統(tǒng)中，功率譜密度函數(shù)無法直接輸入動力系統(tǒng)，最有效的辦法是以軌道隨機不平順的時域樣本作為輸入。翟婉明等[2]提出一種根據(jù)軌道隨機不平順功率譜求出頻譜的幅值和隨機相位，再通過傅里葉逆變換得到軌道隨機不平順時域樣本的數(shù)值模擬方法。除此以外，三角級數(shù)法也是采用較多的一種軌道隨機不平順時域樣本數(shù)值模擬方法[3]。眾所周知，隨機分析計算效率隨著所需樣本數(shù)增加而增加，為了更加高效的研究軌道隨機不平順激勵下的車輛-軌道耦合系統(tǒng)隨機振動特性，需要用更少的隨機樣本來反映軌道隨機不平順的統(tǒng)計特性。無論是常用的三角級數(shù)法還是傅里葉逆變換法，在模擬軌道隨機不平順時都需要大量的隨機變量，而分析所需的隨機樣本數(shù)與隨機變量數(shù)正相關(guān)。因而用更少的隨機變量來描述軌道隨機不平順，可極大提高的車輛-軌道耦合系統(tǒng)隨機振動分析的效率及難度。

文獻[4]中，通過隨機諧和函數(shù)法生成軌道隨機不平順時域樣本，并通過頻率與相位分別隨機以減少所需的隨機變量數(shù)，其中隨機頻率數(shù)為50個，考慮軌道高低、方向、水平、軌距不平順情況下，各種不平順相互獨立，隨機頻率和隨機相位共計需400個隨機變量，所需隨機變量依然很多，且單個隨機樣本中頻率分量減少必然會導致模擬的軌道不平順與實際偏差增大，并不能很好地反映不平順的時域特征，該方法可用于線性車輛-軌道耦合系統(tǒng)中軌道隨機不平順的模擬，但不適用于非線性車輛-軌道耦合動力學系統(tǒng)。

本文基于譜表示方法，將譜表示方法中的標準正交隨機變量表達為基本隨機變量的正交函數(shù)形式，從而實現(xiàn)用較少隨機變量來描述軌道隨機不平順的目的。減少隨機變量數(shù)的同時，保留更多的頻率及相位分量，可更準確反映軌道隨機不平順統(tǒng)計特性以及時域特征，提高非線性車輛-軌道耦合系統(tǒng)隨機振動計算效率。

1 隨機過程的譜表示-隨機函數(shù)模擬

一維、單變量、功率譜密度函數(shù)為S(ω)的隨機過程η(t)可以表示為[5]

(1)

式中：Ut(ω)和Vt(ω)為隨機過程η(t)的譜過程，其增量dUt(ω)和dVt(ω)滿足隨機過程譜表示的基本條件，即

其中，E[·]表示數(shù)學期望；ω、ω′為隨機過程中的頻率；S(·)為雙邊功率譜密度函數(shù)。

對式(1)進行離散，并引入一組標準正交隨機變量{Xk,Yk}，當頻率間隔Δω足夠小時，式(1)可表述為

(2)

式中：dUt(ω)≈ΔUt(ω)=σkXk; dVt(ω)≈ΔVt(ω)=σkYk;σk=[2S(ωk)Δω]1/2，S為軌道不平順功率譜密度函數(shù)；Nω為頻率離散點數(shù)；頻率分量均勻離散情況下，ωk=ωl+k·Δω，Δω=(ωu-ωl)/Nω，ωu、ωl為頻率離散點上、下限值。

(3)

式中：cas(x)=cos(x)+sin(x)為Hartley正交基函數(shù)；隨機變量Θ1、Θ2相互獨立且服從區(qū)間[0，2π)上均勻分布。直接生成的標準正交變量并不能直接代入式(2)，還需要構(gòu)造隨機映射將標準正交變量映射到標準正交隨機變量{Xk,Yk}，由此將描述隨機過程的隨機變量數(shù)降低為2。

2 軌道隨機不平順的譜表示-隨機函數(shù)模擬方法

以TB/T 3352—2014《高速鐵路無砟軌道不平順譜》[6]為例，軌道不平順譜包括高低、方向、軌距和水平四種軌道不平順譜，計算不平順譜時，將實測得到的軌道不平順按1 024 m劃分單元，剔除異常值和趨勢項，加窗處理后進行傅里葉變換(FFT)計算，對所有單元的軌道不平順譜進行擬合，得到軌道不平順功率譜密度函數(shù)。我國的軌道不平順譜擬合公式為

(4)

式中：f為空間頻率；A、k為擬合系數(shù)。

統(tǒng)計誤差可表示為

(5)

需要指出的是，軌道不平順功率譜密度函數(shù)是在大量實測數(shù)據(jù)基礎(chǔ)上，按一定單元長度組成軌道不平順樣本集，由樣本集計算得到的均值譜。單個軌道不平順分段的隨機誤差達到了100%，隨著分段次數(shù)的增加隨機誤差逐漸減小，如果分段100次則隨機誤差可以降到10%[6]。因而采用單個軌道隨機不平順樣本或單元計算得到的不平順功率譜與現(xiàn)有軌道不平順功率譜對比，并不能說明數(shù)值模擬方法的準確性，需采用多個不平順樣本或分段所組成的樣本集不平順功率譜與現(xiàn)有軌道不平順功率譜對比驗證。

對于多樣本模擬，常采用蒙特卡羅法選擇隨機變量，但其計算效率偏低。擬蒙特卡羅法是通過數(shù)論方法產(chǎn)生高度均勻化的擬隨機數(shù)序列代替蒙特卡羅隨機數(shù)模擬，以減少隨機樣本數(shù)從而提高計算效率，但當隨機變量較多，即高維度情況下，數(shù)論選點法選取的隨機數(shù)列不可避免的會產(chǎn)生叢聚現(xiàn)象[7]。采用譜表示-隨機函數(shù)法模擬生成軌道不平順時域樣本，可減少所需的隨機變量數(shù)，避免高維數(shù)論選點中叢聚現(xiàn)象對計算結(jié)果的影響。

綜上所述，采用譜表示-隨機函數(shù)法生成軌道隨機不平順樣本集過程可表示為

Step1根據(jù)生成的軌道隨機不平順類型確定隨機變量數(shù)N，當需同時生成軌道高低、方向、軌距、水平不平順時，由于各類型不平順間相互獨立，因而需要8個隨機變量(各類型不平順均為2個)。

Step2根據(jù)所需隨機變量數(shù)，采用擬蒙特卡羅法生成隨機變量的低偏差序列Θ，常用的低偏差序列有Corput序列、Halton序列、Sobol序列等[7]。本文采用Sobol序列生成隨機變量[5],采用Matlab直接計算Sobol的命令。

Θ={θ1j,θ2j,…,θNj}j=1,2,…,Npt

式中：N為隨機變量數(shù)；Npt為樣本集合的樣本數(shù)。此時生成的隨機變量序列服從[0,1)上均勻分布, 需將其變換為[0，2π)上均勻分布隨機變量Θ′

Θ′=2πΘ

(6)

高低不平順ηV(t)為

方向不平順ηA(t)為

水平不平順ηC(t)為

軌距不平順ηG(t)為

式中：SV、SA、SC、SG分別為高低、方向、水平、軌距不平順功率譜密度：ωk為軌道隨機不平順中不同的頻率分量，其表達式為

(7)

采用譜表示-隨機函數(shù)模擬方法生成樣本數(shù)為300的軌道隨機不平順樣本集，其中軌道不平順功率譜計算方法中單元長1 024 m，采樣間隔0.25 m，單個單元共計4 096個點，因而模擬軌道隨機不平順時頻率離散點數(shù)Nω取4 096。

對生成的軌道不平順樣本集中各樣本，參考文獻[8]中高速鐵路無砟軌道不平順譜的生成方法及軌道不平順譜頻率分布規(guī)律檢驗方法，采用周期圖法計算得到各樣本的不平順譜值，并求得樣本集不平順譜均值，對各頻點不平順譜值Sk,j按下式進行變換

(8)

采用Kolmogorov-Smirnov檢驗方法，對λk,j進行顯著性水平α=0.05、自由度為2的χ2分布檢驗。軌道不平順單樣本情況、單樣本功率譜、樣本均值、樣本標準差、部分頻點功率譜概率分布CDF等評價指標以及樣本集總體功率譜PSD與目標譜對比情況如圖1～圖4所示。

圖1 譜表示-隨機函數(shù)法高低不平順樣本評價指標對比

圖2 譜表示-隨機函數(shù)法方向不平順樣本評價指標對比

圖3 譜表示-隨機函數(shù)法水平不平順樣本評價指標對比

圖4 譜表示-隨機函數(shù)法軌距不平順樣本評價指標對比

從圖1～圖4中可以看出，計算得到的軌道高低、方向、水平、軌距不平順單個樣本功率譜密度值與軌道不平順功率譜存在較大的誤差，而軌道不平順樣本集計算得到的功率譜均值與目標譜吻合較好，其中各種不平順樣本總體均值與目標譜最大相差分別為0.18、0.10、0.07、0.08 mm，標準差最大偏差分別為4.68%、5.77%、6.97%、6.99%，樣本集功率譜最大偏差分別為14.54%、10.23%、10.07%、9.87%，各不平順98.00%、96.33%、98.00%、98.33%的頻點功率譜概率分布滿足自由度為2的χ2分布?？梢钥闯鰳颖炯髦笜伺c目標譜吻合良好，且僅需8個隨機變量便可同時模擬生成軌道高低、方向、水平、軌道不平順時域樣本，同時還具有大量頻率分量。

3 軌道隨機不平順的數(shù)值模擬方法對比分析

目前在工程中應用較多的軌道隨機不平順數(shù)值模擬方法有三角級數(shù)和傅里葉逆變換法。三角級數(shù)法模擬軌道隨機不平順可表示為[2]

(9)

式中：ak為相互獨立的均值為0、標準差為σk的高斯隨機變數(shù)；φk為相互獨立的隨機相位。

傅里葉逆變換法模擬軌道隨機不平順可表示為[1]

η=ifft(X(k))

(10)

式中：X(k)為時間序列，其為復序列，實部關(guān)于Nr/2偶對稱，虛部關(guān)于Nr/2奇對稱。設(shè)模擬總時間為Ts，則時域采樣點數(shù)為Ts/ΔT，一般需要在末尾添0以保證采樣點數(shù)為2的整數(shù)次冪，即為Nr。

X(k)表達式為

k=0，1，…，Nr-1

(11)

式中：ξn為獨立隨機相位序列，ξn=cosφn+isinφn=exp(iφn)，φn服從[0,2π]上的均勻分布；Δf=1/(NrΔT)。

三角級數(shù)法中，ak為服從高斯分布的隨機變量，φk為服從均勻分布的隨機變量，因而三角級數(shù)模擬軌道隨機不平順時所需隨機變量數(shù)為2Nω，且存在高斯分布與均勻分布的兩種類型隨機變量。傅里葉逆變化法中，僅有φn服從正態(tài)分布且相互獨立，而所需隨機變量數(shù)與時域采樣點數(shù)相關(guān)，一般為2的整數(shù)次冪，當采樣間隔0.25 m，模擬1 km軌道隨機不平順時所需隨機變量數(shù)為4 096個。三種方法模擬軌道隨機不平順時域樣本時，不平順頻率點計算方法基本一致，采用相同頻率數(shù)時，具有相同的不平順波長成分，因而文中僅對不平順時域樣本集的幅值進行對比分析。

采用擬蒙特卡羅法生成低偏差隨機向量序列，通過三角級數(shù)法、傅里葉逆變換法模擬得到300個軌道隨機不平順時域樣本。計算樣本集均值、標準差、功率譜，并檢驗樣本集各頻點軌道譜分布，各項對比情況如圖5和圖6所示。各種方法生成的軌道不平順樣本集不平順幅值范圍情況見表1。

圖5 三角級數(shù)法高低不平順樣本評價指標對比

圖6 傅里葉逆變換法高低不平順樣本評價指標對比

表1 樣本集不平順幅值范圍 mm

可以看出，除了傅里葉逆變換法計算得到高低不平順幅值范圍偏小外，三種方法計算得到的軌道隨機不平順幅值范圍大體相當，最大相差5.51%。三角級數(shù)法計算得到的軌道不平順單個樣本功率譜與目標功率譜存在較大誤差，各方向軌道不平順樣本集與目標譜均值最大相差分別為0.22、0.09、0.11、0.08 mm，標準差最大偏差分別為5.14%、6.78%、4.95%、7.25%，功率譜密度最大偏差分別為15.82%、13.93%、10.78%、9.65%，各不平順97.33%、96.67%、97.33%、98.67%的頻點功率譜概率分布均滿足自由度為2的χ2分布。傅里葉逆變換法計算得到的軌道隨機不平順，無論是單樣本還是樣本集，其功率譜值均能很好滿足目標功率譜，最大偏差僅分別為0.30%、0.24%、0.35%、0.13%，但各頻點功率譜概率分布均不滿足文獻[8]中關(guān)于軌道不平順譜分布的論述。除此以外，傅里葉逆變換法各方向軌道不平順樣本集與目標譜均值最大相差分別為0.28、0.11、0.13、0.07 mm，標準差最大偏差分別為11.38%、8.53%、8.73%、6.57%，與其他兩種方法模擬結(jié)果基本一致。綜上可以看出，譜表示-隨機函數(shù)法在模擬軌道不平順時域樣本時其幅值范圍基本一致。

4 非線性車輛-軌道-橋梁耦合系統(tǒng)隨機振動分析

采用譜表示-隨機函數(shù)方法僅需8個隨機變量就可以描述四種軌道不平順，這為非線性車輛-軌道-橋梁耦合系統(tǒng)隨機振動提供了極大的方便。文獻[9]中指出，由于涉及輪軌空間動態(tài)接觸幾何關(guān)系、輪軌切向蠕滑力和輪軌法向力等多種復雜非線性求解過程，車輛-軌道耦合系統(tǒng)最有效的隨機振動研究方法是將軌道不平順譜轉(zhuǎn)換為時域樣本，采用一段足夠長的時域樣本來代替整個隨機過程，即采用蒙特卡洛法選取足夠長的時域樣本。

分別采用三角級數(shù)法、傅里葉逆變換法、譜表示-隨機函數(shù)法模擬得到的軌道隨機不平順樣本，利用三維非線性車輛-軌道-橋梁耦合動力學模型，計算得到系統(tǒng)動力響應，對不同距離上的動力響應結(jié)果進行統(tǒng)計分析，以輪重減載率、脫軌系數(shù)為例，樣本計算結(jié)果截口概率分布PDF見圖7。

圖7 不同計算法的響應截口概率分布

從圖7中可以看出，在300個樣本集的情況下，三種方法計算得到的動力響應其截口概率分布基本相同，都具有較大的離散性，速度300 km/h、不平順波長2～120 m時，輪重減載率變化范圍為±0.25，脫軌系數(shù)變化范圍為±0.12。由于文中采用中國高速鐵路無砟軌道譜，平順性較其他軌道譜好，因而輪重減載率與脫軌系數(shù)的計算結(jié)果偏小。

對于車輛-軌道耦合系統(tǒng)動力響應，采用三種軌道不平順譜數(shù)值模擬方法計算得到結(jié)果基本一致，但所需的隨機變量數(shù)卻存在巨大差異。隨機變量數(shù)的減少方便了樣本空間概率賦予，同時也可在一定程度上減少所需的樣本數(shù)，為采用廣義概率密度演化理論對車輛運行安全性進行評價提供了方便。

對于非線性車輛-軌道-橋梁耦合系統(tǒng)的任一動力響應ξ(t)，車輛運行過程中不考慮其他系統(tǒng)隨機的情況下，其隨機性來源于軌道隨機不平順中隨機參數(shù)Θ的隨機性，該系統(tǒng)為一個概率保守隨機系統(tǒng)。根據(jù)隨機系統(tǒng)的概率守恒原理，其應滿足廣義概率密度演化方程[10]

(12)

研究中常通過與大量蒙特卡洛法(MCM)計算結(jié)果對比，對廣義概率密度演化方法(PDEM)計算結(jié)果進行驗證[4,12]。采用三角級數(shù)法計算模擬軌道隨機不平順樣本，結(jié)合MCM法計算得到3 000個樣本動力響應統(tǒng)計情況，并采用譜表述隨機函數(shù)法，結(jié)合廣義概率密度演化理論，計算得到300個樣本下的動力響應概率密度演化結(jié)果，兩種方法車輛運行距離中間點的截口概率密度對比情況見圖8，由于三角級數(shù)法所需的隨機變量數(shù)較多，因而在3 000個樣本情況下，動力響應結(jié)果依然存在一定離散性，但兩種方法計算得到的截口概率密度基本一致，因而采用譜表示-隨機函數(shù)法模擬軌道隨機不平順，結(jié)合廣義概率密度演化理論可極大提高非線性車輛-軌道-橋梁耦合系統(tǒng)隨機振動分析效率，相較于采用三角級數(shù)等方法模擬軌道隨機不平順并結(jié)合MCM法進行車輛-軌道-橋梁耦合系統(tǒng)隨機振動分析，其計算效率至少可提高一個量級。參考“3σ”準則，取E±3σ對應值為動力響應指標最值。車輛運行過程中輪重減載率及脫軌系數(shù)最大值、等概率密度演化曲線見圖9。由此，譜表示-隨機函數(shù)法軌道不平順數(shù)值模擬方法結(jié)合廣義概率密度演化理論可對非線性車輛-軌道-橋梁耦合系統(tǒng)車輛運行安全性進行評價。

圖8 MCM與PDEM截口概率密度對比

圖9 等概率密度演化曲線

5 結(jié)論

軌道不平順譜是描述軌道隨機不平順最常用的統(tǒng)計函數(shù)，本文以軌道隨機不平順時域樣本數(shù)值模擬方法為研究對象，采用譜表示-隨機函數(shù)法分別由2個隨機變量模擬生成軌道高低、水平、方向、軌距不平順集。

對比樣本集與目標譜均值、標準差、軌道譜以及各頻點軌道譜概率分布，驗證了譜表示-隨機函數(shù)法模擬軌道隨機不平順時域樣本準確性；與三角級數(shù)法、傅里葉逆變化法相比，模擬生成軌道不平順幅值范圍基本一致，誤差在5.51%以內(nèi)，其中傅里葉逆變化法生成的不平順樣本功率譜能很好地吻合目標譜，但各頻點功率譜分布并不滿足與文獻[8]所述頻點功率譜分布；采用譜表示-隨機函數(shù)法模擬生成軌道隨機不平順時域樣本集，結(jié)合廣義概率密度演化理論，可用相對較少樣本計算得到車輛運行過程中概率演化情況，相較于采用三角級數(shù)法模擬軌道隨機不平順，并結(jié)合MCM法進行隨機振動分析，計算效率至少提高了1個量級，為非線性車輛-軌道-基礎(chǔ)結(jié)構(gòu)耦合系統(tǒng)隨機振動分析中軌道隨機不平順的模擬提供了新思路。