結合相角重構法的軌道不平順隨機過程的數值模擬

徐 寧， 任尊松, 李 響

(北京交通大學 機械與電子控制工程學院， 北京 100044)

車輛運用過程中,由于軌道和車輛不斷地相互作用,會使軌面逐漸產生不均勻磨耗, 輪對踏面出現擦傷、不圓等不良損壞，這些破壞又極易形成反作用于車輛-軌道系統的扁疤、表面不圓和波紋形激勵。與此同時，由于鋼軌扣件、軌枕、道床和路基等部件彈性不均,各部件之間存在間隙以及長時間使用的破壞損耗使軌道系統的幾何狀態逐漸發生改變，產生軌道不平順。研究發現軌道不平順是一個隨機過程[1], 是引起車輛-軌道系統隨機振動的主要外界激勵，直接關系到行車安全和車輛運行的舒適性。

為便于在車輛動態模擬與仿真以及進行軌道系統病害診斷等領域應用，德國、美國等通過對軌道不平順進行實測和分析，構造了形式較為統一的譜密度函數[2]。中國鐵道科學研究院以及長沙鐵道學院等研究單位也曾對國內一些個別線路的軌道不平順進行了測量和研究，相應地建立起一定標準的譜密度函數的表達式[3]。這些研究得到的譜密度函數均在頻域范圍內給出，但在車輛運行穩定性、平穩性和安全性分析等方面都直接需要軌道不平順的時間-歷程作為輸入激勵。因此，通過數值模擬獲得軌道譜的時域序列，是車輛和軌道動力學研究中不可回避的問題。很多研究人員也對此開展了研究工作，總結了軌道不平順數值模擬的幾種主要方法，包括白噪聲濾波法[6]、三角級數法[5]、二次濾波法[4]、小波變換法[7]和逆傅里葉變換法[8]。這些方法較為常用，但同時也存在一定的不足，二次濾波法需依據軌道不平順譜密度函數的不同，設計出合理的濾波器；小波變換法獲得隨機時域序列與小波基的選取有關；白噪聲濾波法以及三角級數法需將模擬的不平順隨機序列看作是平穩高斯過程。且通過這幾種方法獲得時域序列經FFT的變換譜與原軌道譜存在一定的差距。逆傅里葉變換法[8]計算速度快，且模擬得到的時域序列經FFT的變換譜與原不平順譜較為一致。但該方法得到的時域隨機序列的正態性與否不確定，樣本數字特征也不確定，這與文獻[9-10]中實測得到的線路不平順為具有一定偏度和峭度范圍的非正態隨機過程并不能完全相符。與此同時，當車輛在處于勻速狀態下，不平順激勵在時間域還是在空間域都是平穩過程。但變速條件下運行時，不平順在時間域上已是非平穩的隨機過程[14]，作為激勵輸入到車輛和軌道系統之中，獲得的各部件響應也較勻速條件下有所差別。為在不同運行工況下，獲得更為準確的車輛系統響應特征，進而分析系統頻域傳遞等特征，獲取變速條件下的不平順隨機序列十分必要。

鑒于此，本文在逆傅里葉變換方法基礎之上，依據實測線路不平順峭度和峰度等數字特征的變化范圍，采用相角重構方法獲取了具有一定峭度值和峰度值的平穩隨機序列。由于變速過程中軌道不平順譜為慢變的時變功率譜，將變速過程看作多個短時段內的勻速過程，使時變功率譜轉化為各時段上的平穩功率譜，并通過逆傅里葉變換法分別得到各時段對應的隨機時域序列，利用相角重構方法實現各序列的連接，得到加速和變加速過程的非平穩軌道不平順時域序列。

1 具有一定偏度值和峭度值的隨機時域序列獲取

1.1 由逆傅里葉變換法生成隨機時域序列

根據逆傅里葉變換方法[8]進行軌道譜的數值模擬，其實現步驟如下:

Step1首先將不平順單邊功率譜密度函數轉化為雙邊譜密度函數。進而將得到的空間雙邊不平順功率譜S(Ω)轉化為時間功率譜S(f)。設車輛運行的速度為v，軌道不平順的最大波長為λmax，最小波長為λmin，則模擬的頻率范圍最大值為fmax=v/λmin，最小值為fmin=v/λmax。根據采樣定理，時域采樣間隔應滿足Δt≤0.5fmax。假定總的模擬時間為T，那么時域采樣點數為T/Δt，一般需保證其個數Nr為2的整數次冪，頻率采樣間隔為Δf=1/(NrΔt)，有效頻段內采樣點數為Nf=(fmax-fmin)/Δf，由周期圖法估計出的功率譜為偶對稱序列，在Nr/2-Nf～Nr/2+Nf頻率區間內的采樣點記為零，得到在1～Nr/2上功率譜的Nf個離散采樣點值S(f=kΔf)，k=0，1，…，Nr/2及在Nr/2+1～Nr-1上離散采樣點值S(f=kΔf)，k=Nr/2+1，…，Nr。

Step2根據周期圖法，可獲得各頻率采樣點的頻譜模值為

X(k)=DFT[x(n)]=

( 1 )

Step3所求時域序列x(n)為一隨機過程，要求其FFT得到的復序列X(k)的相位需具有隨機性。設ψn是均值為零的獨立相位序列，那么ψn應為復數，且滿足|ψn|=1,即

ψn=cosβn+i sinβn=exp(iβn)

( 2 )

式中：βn為在[0,2π]上服從均勻分布的隨機相角。且X(k)的虛部關于Nr/2奇對稱，實部關于Nr/2偶對稱，因此只需求出k在0～Nr/2間各頻率采樣點處的頻譜值X(k)，其表達式為

( 3 )

Step4將式( 3 )中得到的序列X(k)進行傅里葉逆變換，可獲得相應的不平順時域序列。

( 4 )

1.2 通過重構相角生成具有一定統計特性的隨機時域序列

對均值為零的隨機時域序列，根據偏度和峭度的定義，其表達式為

( 5 )

( 6 )

對于隨機位移-時間歷程x(t)，其各階原點矩均可以通過積分逼近，即

( 7 )

將式( 7 )代入到式( 5 )和式( 6 )之中，按照文獻[11]方法，經過化簡，序列的偏度和峭度可表達為關于功率譜幅值和隨機相角的形式

( 8 )

( 9 )

式中：Ak、Aj、Am、An為功率譜中各頻率點的幅值；φk、φj、φm、φn是[0,2π]上均勻分布的隨機相角，經證明[11]，為不影響對功率譜的精確模擬，可使得式( 8 )和式( 9 )中Ak和Aj等各幅值不變，通過調整隨機相位角中一部分進而改變生成隨機時域序列的峭度和偏度值，也使得相位角的整體隨機性得以保證，進而確保重構相角后得到的時域序列仍為隨機序列。

例如對式( 8 )，取下標為j=2k的各組相位角φk和φj(如{φ1，φ2}，{φ2，φ4}，…)符合關系φj=2φk，那么cos(φj-2φk)的值恒為1，偏度值S逐漸增大；或符合關系φj=2φk+π，使cos(φj-2φk)的值恒為-1，偏度值S逐漸減小。

1.3 改進相角重構方式獲取隨機不平順序列

逆傅里葉變換方法應用過程中，生成的隨機相角是偽隨機數，以至于得到的隨機序列雖然均值接近于零，但往往并不是高斯序列，峰度值和峭度值大小也不確定。且軌道譜本身為頻變譜甚至為折線譜，而非文獻[11]中所采用的平直功率譜，所以僅采用文獻中絕對正向和絕對反向這樣兩種相角重構方式重構相角，往往會造成找尋效果比較差；甚至某些情形下，無法獲得目標參數下的隨機時域序列，鑒于此，以文獻[11]中方法為基礎，并改進原有方法，可獲得兼具目標偏度值S*和峭度值K*的隨機序列，具體實現過程步驟:

Step1根據偏度S和峭度K表達式中第一個三角函數項，在已有的Nr/2個[0,2π]上均勻分布的隨機相位角中找出所有下標滿足關系j=2k和m=3n(0≤j,k,m,n≤Nr/2-1)的各組相角組合{φj,φk}和{φm,φn}。數據長度較大時，滿足上述條件的相角很多，其集合可記作{φj1,φk1;φj2,φk2;…}和{φm1,φn1;φm2,φn2;…}。

Step3繼續按照Step2中的方式，選取后續的多對相角。例如第x次重構的兩組相角分別為{φjx,φkx}和{φmx,φnx}，為減小重構過程中偏度和峭度間的相互影響，當φjx為6的公倍數的情形，跳過該組相角，選取下一組對應相角{φj(x+1),φk(x+1)}。再次計算新生成序列的偏度和峭度值Sxil和Kxil，并按照步驟Step2中的方式進行檢驗，如e(Sxil)和e(Kxil)處于規定的范圍之內，停止重構，獲得滿足要求的隨機序列。如果仍超出規定的范圍，可按照Step2中的方式，確定該次重構的最優方式。并對后一組相角進行類似的重構，直到e(Sxil)和e(Kxil)小于規定限值ξS和ξK為止。

2 變速過程的非平穩隨機時域序列的獲取

勻速運行時，軌道不平順激勵可視為平穩的隨機過程。當車輛處于加速啟動和減速制動等運行工況下，軌道不平順在空間域上是平穩的，但在時域上是非平穩的隨機過程。此時，不平順譜轉變為慢變的時變功率譜。首先，將模擬序列總的時長拆分為多個較短的時間段，在這些較短的時間段上，時變功率譜可近似看做非時變，可求得各時間段上的功率譜密度函數表達式，通過逆傅里葉變換方法可獲得各時段的隨機序列?；诖?，利用相角重構的方法將各段隨機序列首尾相連接，即可獲得變速過程的非平穩隨機時域序列，其具體實現過程如下。

2.1 軌道譜轉化為各短時間段上平穩功率譜

根據時域和空間域的轉化關系式

ft=Ωs

(10)

式中：Ω為空間頻率;s為走行距離。瞬時頻率表達式為

fs=Ωds/dt

(11)

車輛在勻加速情況下

(12)

車輛在變加速度情況下

(13)

式中：v0、a0分別為速度和加速度的初始值；a、u分別為加速度和加速度的變化率。

可得到軌道不平順在勻加速和變加速工況下的時變功率譜，即

(14)

(15)

子段劃分數M的選取，應使得各子段任意頻率點對應的時變功率譜S(f,t)滿足式(16)，且M的取值一般為2的整數次冪。

(16)

式中：fj為功率譜頻域范圍內任一頻率點;timin和timax為第i個子時間段的上、下時刻點;δ為一很小的正數。

設模擬的總時間為T，已劃分時間段數為M，那么各子段的時長Ts0可暫定為T/M，根據時域采樣間隔為Δt，子段時域采樣點數為Ts0/Δt，若滿足2j-1

根據得到各子段時長Ts，計算出各子段中心時刻Ti mid(i=1～M)，依據式(14)和式(15)，將時變功率譜S(f,t)轉化為各子段的中值時刻功率譜S(f,timid)?；诘玫降母髦兄禃r刻功率譜S(f,timid)，并利用1.1中逆傅里葉變換法，即可獲得各子段的隨機時域序列。

2.2 基于相角重構實現各段隨機序列連接

通過2.1節的過程得到了各子段的隨機時域序列，但這些序列在每一段的末尾時刻以及下一段的初始時刻的函數值并不連續，顯然與真實情況不相符。參照相角重構的方法，實現各子段隨機序列的連接按照以下具體步驟:

Step2利用逆傅里葉變換法重新生成隨機序列，經過第j個相角重構后，獲得的該子段序列的初始時刻位移值x(i_j)start_k，按式(17)計算出各x(i_j)start_k與x(i-1)last的絕對誤差。

δi_j_k=|xi_jstart_k-x(i-1)last|

(17)

式中：i為當前子段數;j為該段第j個重構相角;k為該相角第k個重構方向。

Step3將已得到的各δi_j_k與規定限值ξ進行比較，若某一δi_j_k小于規定限值，那么對應的時域序列即為所求，如果k=0～2r的任一重構方式均不能小于規定限制ξ，那么比較各δi_j_k(k=0～2r)，求出其中最小值min(δi_j_k)，并將其對應重構方式作為本次重構的最優方式。選取下一重構相角φj+1，按照Step1的方法進行重構，直到得到δi_j_k小于規定限值ξ為止。

Step4以步驟Step3得到序列的末尾時刻位移值xilast為基準，多次重復Step1～ Step3，實現后續多段隨機序列的連接，直到最后一段為止。

3 數值算例及結果分析

為驗證本文方法的有效性和實用性，結合兩個不同算例，利用重構相角方法對高速鐵路無砟軌道譜進行時域轉化。

3.1 算例1

根據文獻[9]的建議，模擬的軌道長度為1 024 m，而文獻[10]中通過正線實測高低不平順數據并進行統計分析，得到不平順序列的偏度值范圍為-1～1之間，峰度值位于零值附近。為不失一般性，本文選取生成隨機序列的目標峭度為0，目標偏度取值分別為0、-0.5和0.5，模擬的運行速度為v為80 m/s(288 km/h)，時域采樣間隔Δt為10-4，軌道不平順波長范圍為3～100 m，規定限值ξS和ξK分別為0.01和0.02。

取高速無砟軌道高低不平順譜為轉化對象，其表達式為分段冪函數，具體可見式(18)，各分段對應的相關參數見表1。

SΩ=A/Ωk

(18)

式中：Ω為空間頻率;A、k為擬合公式系數，1/m。

表1 無砟軌道不平順中位譜擬合公式系數 1/m

高速鐵路無砟軌道譜為分段線性譜,各段間分段點對應空間波長和頻率可見表2。

表2 高速鐵路無砟軌道不平順譜分段點空間頻率及對應波長

圖1、圖2分別為逆傅里葉變換法得到的隨機序列以及不同目標峭度和目標偏度下的隨機時域序列。

從圖1、圖2可以看出，經過本文方法調整得到的具有一定偏度和峭度值的隨機序列在數值范圍上與逆傅里葉變換法直接獲得的隨機序列相近，說明了該方法的可用性。不同偏度和峭度值對應的時域序列的整體走勢有很大的不同，說明了偏度和峭度對于隨機序列確實有一定的影響。

表3為各隨機序列的均值、方差的統計特性。表中顯示重構相角后得到的新隨機序列的均值和方差與原序列幾乎沒有差別，并能保證序列的偏度和峭度較高的精度地接近目標值。隨著偏度值從-0.5～0，再由0逐漸增大到0.5，各隨機序列的中位數也對應地出現了逐漸增大的過程，偏度為0.5序列與偏度為-0.5序列的中位數相差達到1.25，峭度值由基本序列的-0.103左右到達0附近，數值處于均值±2周圍的點的個數多了近8 000。統計特征的對比結果說明偏度和峭度對隨機序列整體特點確有一定的影響。

表3 勻速條件下各序列統計特征及其對比

由圖3可知，重構相角后得到的新隨機序列及逆傅里葉變換法得到初始序列，經FFT的幅頻譜與轉化譜都非常接近，說明本文方法的模擬結果精確性。

3.2 算例2

仍以算例1中的高速鐵路高低不平順為轉化對象，模擬的軌道長度約為2 048 m，以v0=77 m/s(277 km/h)作為車輛在加速過程中的初始速度，將之后的加速過程分別看做是加速度為a=0.2 m/s2的勻加速運動與加速度為a=0.2 m/s2、加速度變化率為u=0.01 m/s3變加速運動兩種情形，并假設二者各子段的初始隨機相位角相同，獲得勻加速過程和變加速過程對應的不平順位移-時間序列見圖4。

由圖4可以看出，勻加速過程和變加速過程隨機不平順序列整體趨勢和數值范圍上比較一致，但各局部極值的大小及最大值和最小值出現的位置均存在差異，且隨著時間的增長，二者在局部上差異更為明顯，說明加速度變化率的對隨機序列的走勢有一定影響。

表4給出勻加速和變加速過程下兩隨機序列的統計特征及其平穩性檢驗結果。

表4 隨機序列的統計特征及平穩性檢驗結果

由表4見，重構相角方式獲得的勻加速和變加速過程的軌道隨機序列是非平穩的，而系統的均值仍然在零值附近，方差與勻速過程軌道隨機序列相近。

由圖5可以看出，由重構相角方法獲得的勻加速過程和變加速過程的隨機序列，經FFT變換得到的功率譜，二者數值上略有差別，且都與原不平順功率譜較為接近，說明本文方法獲得的模擬序列的可用性和準確性。

4 結論

(1) 本文通過二次重構隨機相角的方法可以比較準確的將由逆傅里葉變換法得到的不平順隨機序列轉化具有一定偏度和峭度值的隨機序列，能夠保證序列的均值和方差等統計特征不發生較大改變。偏度和峭度對時域序列的整體走勢，中位值以及均值附近±2周圍的樣點的個數均有不小的影響。

(2) 將加速過程和變加速過程下時變功率譜轉化為一些較短時間段上平穩功率譜，利用逆傅里葉變換法獲得各時間段的隨機時域序列，通過重構相角方法實現各段序列的連接，最終獲得了變速過程下的軌道隨機不平順序列。發現變速過程下不平順序列具有非平穩的特點，其均值和方差等統計特性與勻速過程下較為接近。且變速過程的隨機不平順序列經FFT變換得到的功率譜與原不平順功率譜數值上較為一致。勻加速過程和變加速過程下隨機不平順序列的整體趨勢和數值范圍上比較一致，但各局部極值的大小及最大值和最小值出現的位置均存在差異，且隨著時間的增長，局部上差異越明顯。

參考文獻：

[1] 羅林.軌道隨機干擾函數[J]. 中國鐵道科學,1982,3(1):74-112.

LUO Lin. Track Random Irregularities Function[J]. China Railway Science,1982,3(1):74-112.

[2] 翟婉明.車輛-軌道耦合動力學[M].2版.北京：中國鐵道出版社,2002.

[3] 長沙鐵道學院隨機振動研究室.關于機車車輛-軌道系統隨機激勵函數的研究[J].長沙鐵道學院學報,1985,1(2):1-36.

Research Group on Random Vibration of Changsha Railway Institute. Research on Random Excitation Function for Rolling Stock/Track System[J]. Journal of Changsha Railway University,1985,1(2):1-36.

[4] 王元豐,王穎,王東軍.鐵路軌道不平順模擬的一種新方法[J].鐵道學報,1997,19(6):110-115.

WANG Yuanfen, WANG Ying, WANG Dongjun. A New Simulating Method for the Irregularities of Railway Tracks[J]. Journal of the China Railway Society,1997,19(6):110-115.

[5] 星谷勝,常保琦.隨機振動分析[M].北京:地質出版社,1977:42-62.

[6] 錢雪軍.軌道不平順的時域模擬法[J].鐵道學報,2000,22(4):94-98.

QIAN Xuejun. Track Irregularity Simulation in Time Domain[J]. Journal of the China Railway Society,2000,22(4):94-98.

[7] 何越磊,陳施宇,李再幃.基于小波變換的軌道不平順數值模擬方法[J].城市軌道交通研究, 2014,17(10):59-62.

HE Yuelei, CHEN Shiyu, LI Zaiwei. Track Irregularity Simulation Based on Wavelet Transform[J]. Urban Mass Transit,2014,17(10):59-62.

[8] 陳果,翟婉明. 鐵路軌道不平順隨機過程的數值模擬[J].西南交通大學學報,1999,34(2): 138-142.

CHEN Guo, ZHAI Wanming. Numerical Simulation of the Stochastic Process of Railway Track Irregularities[J]. Journal of Southwest Jiaotong University,1999,34(2):138-142.

[9] 康熊,劉秀波,李紅艷,等.高速鐵路無砟軌道不平順譜[J].中國科學:技術科學,2014,44(7):687 -696.

KANG Xiong, LIU Xiubo, LI Hongyan, et al. PSD of Ballastless Track Irregularity of High Speed Railway[J]. Scientia Sinica:Technologica,2014,44(7):687-696.

[10] 中國鐵道科學研究院基礎設施檢測研究所.2013年JC字第157號 高速線路不平順譜的研究[R]. 北京：中國鐵道科學研究院基礎設施檢測研究所，2013.

[11] 蔣瑜.頻譜可控的超高斯隨機振動環境模擬技術及其應用研究[D].長沙：國防科學技術大學,2005.

[12] 蔣瑜，陶俊勇，王得志，等.一種新的非高斯隨機振動數值模擬方法[J].振動與沖擊, 2012, 31(19):169-173.

JIANG Yu, TAO Junyong, WANG Dezhi,et al. A Noval Approach for Numerial Simulation Method of a Non-Gauss Random Vibration[J]. Journal of Vibration and Shock,2012,31 (19):169-173.

[13] 李茂清，王潔，陳強，等. 基于MATLAB程序的FIR濾波器設計實現[J].電力學報，2008，23(2): 87-90.

LI Maoqing, WANG Jie, CHEN Qiang, et al. The Accomplishment of the Design of FIR Filter Design Based on MATLAB Program[J]. Journal of Electric Power,2008,23(2):87-90.

[14] 張立軍.車輛非平穩行駛動力學及控制研究[D].沈陽：東北大學，2006.