牛留斌，劉金朝，李　谷，祖宏林，楊　飛

(1.中國鐵道科學研究院 基礎設施檢測研究所，北京　100081； 2.中國鐵道科學研究院 機車車輛研究所，北京　100081)

輪軌力是運營車輛對軌道狀態(tài)的動態(tài)響應，它是評價車輛運行安全的重要參數。任何影響車輛運行安全的軌道缺陷都會在輪軌力檢測數據中有所反映，通過分析輪軌力檢測數據、定位查找線路缺陷、消除車輛和軌道部件疲勞損傷的根源成為一種新興的軌道檢測技術[1]，安裝于軌道檢查車上的輪軌力檢測系統為獲取實時輪軌力檢測數據提供了技術手段。

輪軌力受到車輛及其運行速度、軌道空間狀態(tài)、輪軌作用形式等多方面隨機因素的影響，輪軌力的函數形式不能通過單一確定性變量表達，而應從其幅值和頻率的分布特性等方面進行研究。功率譜分析方法是研究輪軌力的有效工具之一，輪軌力功率譜從幅值和頻率2個方面反映了輪軌力在時域和頻域的分布特征和規(guī)律。在既往輪軌力及其功率譜分布規(guī)律的研究中，輪軌力數據主要是通過車輛—軌道耦合模型的數值計算[2]、成熟商業(yè)仿真軟件仿真[3]、檢測系統實測[4]或者直接在特定的軌道斷面上布置測試傳感器[5]得到，但耦合模型或仿真模型中仿真參數的設定和軌道輸入與實際軌道狀態(tài)的關聯程度均會制約著輪軌力計算和仿真的精度，在軌道特定位置安裝傳感器實測不同速度條件下輪軌力的數據量也會影響輪軌力功率譜分布特征的研究。截至目前，國內外還未建立類似于軌道不平順譜的具有普適意義的輪軌力功率譜。

本文從有砟軌道上不同速度條件下的輪軌力檢測數據時域波形重復性入手，在分析輪軌力功率譜特征分布的基礎上，利用非線性最小二乘法(Nonlinear Least Squares)中的Levenberg-Marquardt算法，提出了一種結構簡單、適用性強的輪軌力功率譜有理式函數擬合模型，并給出了有砟軌道輪軌力功率譜擬合結果，以驗證所提擬合模型的實用性。

1　輪軌力檢測數據的時域分析

本文采用的輪軌力檢測數據來源于SY998799型安全綜合檢測車上配備的輪軌力檢測系統，該檢測車最高運行速度為160 km·h-1。輪軌力檢測系統安裝有中國鐵道科學研究院研制的“單周期雙橋路正弦余弦合成法”測力輪對，具有實時連續(xù)檢測功能。輪軌力檢測數據的采樣頻率為500 Hz。

圖1為檢測車3次以140 km·h-1速度通過中國鐵道科學研究院試驗場環(huán)形鐵道某段200 m長有砟軌道區(qū)段時的輪軌力檢測數據波形對比圖。從圖1可以看出：3組輪軌力檢測數據的波形吻合良好，變化趨勢一致，具有相似特性，反映了該區(qū)段相同軌道狀態(tài)引起的輪軌力響應是相對穩(wěn)定，也就是在車輛、運行速度等因素相對固定的條件下，輪軌力主要受到軌道狀態(tài)的影響，相同區(qū)段軌道狀態(tài)引起的輪軌力波形具有較高的相似性。

圖1相同速度條件下相同區(qū)段上3組輪軌力檢測數據的波形對比

圖2為檢測車以不同的速度通過環(huán)形鐵道某段100 m長有砟軌道區(qū)段時的輪軌力檢測數據波形對比圖。從圖2可以看出：檢測車以不同速度經過同一軌道時輪軌力的波形具有很強的相似性，特別是在輪軌力波形峰值處；如圖2(a)中在軌道短波不平順引起的輪軌垂向力雖然峰值大小不相等，但輪軌垂向力峰值波形顯示出同樣的沖擊特征，圖2(b)中相同的橫向軌道不平順引起的輪軌橫向力波動特征很明顯，這說明了相對于固定的軌道狀態(tài)，檢測車對軌道的輪軌力響應特征是相對固定的，檢測車運行速度的變化并不改變輪軌力波形的形狀，僅會引起輪軌力波形峰值的變化。

由于低頻輪軌力與測力輪對的靜輪重及軌道的長波不平順有關，且不影響輪軌力分布特征的描述，因此在下文中濾除了輪軌力檢測數據中1 Hz及以下頻段的低頻數據成分。圖3為檢測車運行速度為100 km·h-1和通過環(huán)形鐵道某段100 m長有砟軌道區(qū)段時25萬個輪軌力樣本數據分布的頻數直方圖及其擬合曲線。從圖3可以看出：在該速度條件下，輪軌力檢測數據分布的頻數直方圖呈現正態(tài)分布特征，其中輪軌垂向力的頻數直方圖擬合曲線服從0均值正態(tài)分布，而輪軌橫向力的則服從正偏右拖尾正態(tài)分布。

圖2　不同速度條件下輪軌力檢測數據的波形對比

圖3100 km·h-1速度條件下輪軌力檢測數據分布的直方圖及擬合曲線

采用標準GB/T 4882—2001《數據的統計處理和解釋 正態(tài)性檢驗》[6]中的非參數愛潑斯—普利(Epps-Pully)檢驗方法對輪軌力檢測數據是否符合正態(tài)分布進行假設檢驗。對于n組輪軌力檢測數據x的檢驗統計量TEP為

(1)

其中，

根據文獻[6]查表得到由顯著性水平α(如0.05)和樣本量n確定的檢驗值Tn，如果輪軌力檢測數據的檢驗統計量TEP大于Tn，則表示輪軌力檢測數據不符合正態(tài)分布；否則，表明輪軌力檢測數據符合正態(tài)分布。

本文采用的輪軌力檢測數據經過Epps-Pully檢驗后、在顯著性水平α=0.05的條件下，輪軌力檢測數據的檢驗統計量TEP小于檢驗值Tn，即輪軌力檢測數據符合正態(tài)分布，與圖3中頻數直方圖擬合曲線呈正態(tài)分布的結果一致。

圖4為不同速度條件下輪軌力檢測數據正態(tài)分布時的概率密度曲線。從圖4可以看出：隨著檢測車運行速度的增加，輪軌力檢測數據正態(tài)分布概率密度曲線的頂峰越平緩，即曲線的峰值逐步減??；圖2和圖4都顯示了相對于相同的軌道狀態(tài)，不同速度條件下輪軌力檢測數據的波形在時域上具有相似性。

功率譜分析是建立在數據平穩(wěn)性假設的基礎上，因此進行輪軌力檢測數據平穩(wěn)性檢驗的目的是檢驗該信號是否屬于平穩(wěn)隨機過程。

平穩(wěn)性檢驗的方法[7]有逆序檢驗法、游程(輪次)檢驗法、特征根檢驗法等，本文采用游程檢驗法進行輪軌力檢測數據的平穩(wěn)性檢驗。游程檢驗的步驟為：在保持隨機序列原有順序的情況下，將輪軌力檢測數據分成相互排斥的2組符號序列，比如大于序列中位數的數據標記為“+”，小于中位數的標記為“-”，從而得到1組符號序列，將N1和N2分別記為符號序列內標記為“+”和“-”的個數，N為兩者之和，而游程檢驗法中的游程定義為具有相同元素符號的序列，即把夾在異類元素符號之間的同類元素定義為1個游程。統計上述輪軌力檢測數據符號序列中的游程數r。按照文獻[7],平穩(wěn)隨機過程的游程統計量Z應近似服從標準正態(tài)分布N(0,1)，即Z～N(0,1)，且

圖4不同速度條件下輪軌力檢測數據的概率密度分布曲線

(2)

其中，

在給定顯著性水平為α的情況下，游程統計量Z落入標準正態(tài)分布N(0,1)的水平分別為1-α/2和α/2時，其由分位數u1-α/2，uα/2組成的置信區(qū)間為[u1-α/2,uα/2]，則認為所檢驗的輪軌力檢測數據屬于平穩(wěn)隨機過程；如果游程統計量Z落入置信區(qū)間[u1-α/2,uα/2]之外，則認為輪軌力檢測數據不屬于平穩(wěn)隨機過程。

利用游程檢驗法對輪軌力檢測數據進行平穩(wěn)性檢驗，在顯著性水平α=0.05的條件下，對應的置信區(qū)間為[-1.96，1.96]，當由式(2)計算出輪軌力檢測數據游程統計量Z滿足|Z|<1.96時，即認為本文采用的輪軌力檢測數據屬于平穩(wěn)隨機過程。

2　輪軌力檢測數據的頻域分析

本文采用修正周期圖譜法[8](Welch功率譜分析方法)得到輪軌力檢測數據的功率譜密度曲線，Welch方法先對輪軌力檢測數據分段，因相鄰段內數據的重疊率為50%，傅里葉變換長度選為4 096。為了減少功率譜密度能量的泄露，在做功率譜分析時需對輪軌力檢測數據進行加窗處理，本文采用加漢明窗(Hamming)進行處理。

圖5為100 km·h-1速度條件下輪軌力檢測數據的功率譜曲線，該數據在時域上對應的頻數直方圖即為圖3。

從圖5可以看出：輪軌力功率譜曲線包含有周期性諧波成分和隨機成分；周期性諧波成分對應輪軌力功率譜曲線中的高頻尖峰處，如輪軌力功率譜曲線在10，20和30 Hz等處出現規(guī)律性的尖峰值，對應的周期性諧波成分的波長分別為2.78，5.56和8.34 m等，這個周期性諧波的波長是與測力車輪的周長相關，是由車輪的不圓順引起的；隨機成分對應輪軌力功率譜中的低頻成分，其變化平緩、特征明顯，顯示了輪軌力檢測數據在頻域上的變化規(guī)律。

不同速度條件下輪軌力檢測數據的功率譜曲線如圖6所示。從圖6可以看出：在軌道狀態(tài)相對穩(wěn)定的條件下，由于速度的差異造成輪軌力功率譜曲線中周期性諧波成分的峰值出現的頻率不一樣，但是各周期性諧波成分對應的波長是一致的，即與測力車輪的周長相關；輪軌力功率譜曲線中隨機成分的波形具有很強的相似性，可以通過相同的曲線模型對其進行擬合，并且在輪軌力功率譜曲線擬合前應先消除周期性諧波成分的影響。

圖5100 km·h-1速度條件下輪軌力檢測數據的功率譜曲線

圖6不同速度條件下輪軌力檢測數據的功率譜密度曲線

3　輪軌力功率譜曲線擬合模型

相對于輪軌力功率譜分析，輪軌力功率譜模型的研究尚處于起步階段，國內外鮮有實用的輪軌力功率譜模型報道。本文提出采用有理式函數得到輪軌力功率譜曲線的擬合模型(Rational Wheel-rail Force Model，RFM)。有理式函數定義為2個多項式函數相除，因此首先定義m階多項式Fm為

Fm=F(f,p,m)=fm+p(1)fm-1+

p(2)fm-2+…+p(m)

(3)

式中：f為頻率，Hz；p為由多項式函數Fm中待定系數向量,共有m個值，即p(1),p(2),…，p(m)。

經過多次擬合分析，將RFM模型的表達式定義為3階多項式函數F3與4階多項式函數F4的比值，即

(4)

式中：S(f，c)為輪軌力功率譜的曲線擬合模型；p3和p4分別為3階多項式函數F3與4階多項式函數F4中的待定系數向量，它們的長度分別為3×1和4×1；A為RFM模型的放大系數；c為RFM模型待定系數組成的向量，即由p3,p4,A組成的向量，它的長度為8×1。

由式(4)可以看出：RFM模型的待估參數少，結構簡單。

本文采用非線性最小二乘法中的Levenberg-Marquardt算法估計RFM模型的待定參數c。

在進行曲線擬合前，為了使后續(xù)的曲線擬合結果能更真實地反映輪軌力功率譜的信息狀態(tài)，應先對輪軌力功率譜曲線進行平滑處理，消除初始輪軌力功率譜中的周期性諧波成分和奇異點對擬合結果的影響。

3.1　平滑處理方法

本文采用五點三次平滑算法對輪軌力檢測數據功率譜曲線(初始曲線)y進行平滑處理，得到用于擬合的輪軌力功率譜平滑曲線Y。五點三次平滑是一種局部加權平均的方法，算法如下。

Yi=

(5)

式中：i為輪軌力功率譜的點數，i=1，2，…，q，它的值與Welch譜分析方法中選取的傅里葉變換點數有關；yi為第i個初始輪軌力功率譜數據；Yi為第i個平滑處理后的輪軌力功率譜數據。

為了得到輪軌力功率譜平滑曲線擬合數據，以上五點三次平滑過程需要進行多次迭代。

3.2　擬合算法

通過采用非線性最小二乘擬合方法中的Levenberg-Marquardt算法對輪軌力功率譜平滑曲線Y進行擬合。該算法[9-10]是對高斯—牛頓擬合方法進行了修正的非線性優(yōu)化方法，既具有高斯—牛頓法的局部收斂性，又具有梯度下降法的全局特性，它通過自適應調整阻尼因子λ達到收斂。

首先設置擬合曲線系數向量c的初始值c0(在RFM模型中，可選8個隨機數據作為c0的初始值)，初始權重向量h0，目標誤差ξ，迭代次數M，初始阻尼因子λ0，阻尼調節(jié)因子ν(本文選擇整數10)等控制參數。則當輪軌力功率譜曲線擬合到第k步時，輪軌力功率譜為Yk,擬合系數為ck，權重向量為hk，阻尼因子為λk。當k≥m時，停止迭代，輸出擬合系數ck；當k

步驟1：計算第k步輪軌力功率譜擬合誤差εk，即

εk=‖Yk-S(f,ck)‖

(6)

步驟2：計算系數矩陣Ck的雅克比(Jacobi)矩陣(切映射)Jk，構造增量正規(guī)方程為

(7)

其中，

式中：δ為增量正規(guī)方程的根；I為單位矩陣。

步驟3：求解增量正規(guī)方程的根δk，判斷如下。

(1)若‖yk-S(k,ck)‖≥εk， 則阻尼因子為λk+1=νλk， 重新求解增量正規(guī)方程的根δk+1， 且令ck+1=ck+δk+1后返回步驟1， 并重新調整權重向量hk， 使得靠近擬合曲線的輪軌力功率譜權重大，遠離擬合曲線的輪軌力功率譜數據權重小，即調整權重hk有利于消除離散邊緣數據點對擬合曲線的干擾，提高曲線擬合質量。權重向量的調整方法如下。

定義經第k次迭代后的輪軌力功率譜Yk與擬合的輪軌力功率譜曲線S(f，c)之間的向量差為Dk，其加權平方和為Lk，即有

Dk=Yk-S(f,ck)

(8)

(9)

式中：hk，l為第l個權重向量hk的值。

銷售渠道的轉型升級。采用內部轉化的方法，將傳統的發(fā)行渠道進行揚棄，改造提升傳統出版渠道，逐步提高傳統渠道轉化為數字渠道的比例；獨立自主地建構數字產品銷售渠道，建立健全個人用戶、機構用戶客戶關系管理系統，持續(xù)推動數字出版的市場化運營和產業(yè)化發(fā)展。

(10)

其中，

式中：hk+1，l為第l個調整后權重向量hk+1的值；t為修正差值變量。

權重向量調整為hk+1后，輪軌力功率譜Yk根據各點權重的變化，生成新的輪軌力功率譜Yk+1，進行新一輪的迭代計算。

(2)若‖yk-S(k,ck)‖<εk， 且‖δk‖≥ξ， 則調整阻尼因子， 令λk+1=λk/ν，εk+1=εk， 返回步驟2，重新構造增量正規(guī)方程并求解。

(3)若‖yk-S(k,ck)‖<εk， 且‖δk‖<ξ， 令ck+1=ck+δk，進行步驟4。

步驟4：停止迭代，并將最終權重向量記為H、最終擬合系數記為C。

迭代步驟結束時，得出平滑輪軌力功率譜Y的擬合曲線S(f,C)，兩者的相似程度用擬合優(yōu)度進行檢驗。

3.3　擬合優(yōu)度

(11)

擬合優(yōu)度R2的取值范圍為0～1.0，R2的值越接近1.0，說明擬合效果越好。

4　實測數據驗證

對圖5所示100 km·h-1速度條件下的輪軌力功率譜初始曲線進行五點三次平滑處理，得到其平滑曲線如圖7所示。從圖7可以看出：平滑處理消除了初始輪軌力功率譜數據中的奇異點及周期性諧波成分，形成規(guī)律性明顯的輪軌力功率譜隨機成分，它反映了軌道狀態(tài)引起的輪軌力特征。

圖7　輪軌力功率譜曲線平滑結果

利用RFM模型對圖7中的平滑輪軌力功率譜曲線進行擬合，得到的輪軌垂向力和橫向力功率譜擬合曲線的擬合優(yōu)度R2分別為0.97，0.95，結果如圖8所示。由圖8可以看出：輪軌力功率譜平滑后，高頻成分波動性較大，如200 Hz以上部分，在數據擬合過程中會影響到整體的擬合優(yōu)度，在這種情況下，可以對輪軌力功率譜進行分段擬合。

圖8　輪軌力功率譜擬合結果

圖8中輪軌力功率譜擬合曲線與其98%置信度上、下限曲線的對比結果如圖9所示。從圖9可以看出：輪軌力功率譜中的隨機成分大部分位于擬合曲線置信度上、下限曲線之間。

圖9輪軌力功率譜擬合曲線與其98%置信度上、下限曲線的對比

同理得到圖6所示不同速度條件下輪軌力功率譜的擬合曲線如圖10所示，且擬合曲線的擬合優(yōu)度R2均大于0.94。從圖10可以看出：不同速度級下的輪軌力功率譜變化趨勢相似，反映了圖6中輪軌力功率譜曲線族中隨機成分的變化規(guī)律。

為進一步對RFM擬合模型進行驗證，利用以上擬合步驟對某條有砟軌道線路上不同速度條件下的輪軌垂向力功率譜曲線進行擬合，得到的擬合曲線(見圖11)擬合優(yōu)度R2均不小于0.93。

對比圖10與圖11可以看出：不同的有砟軌道線路狀態(tài)對應的輪軌力響應狀態(tài)是不同的，如圖11中，該有砟軌道上輪軌垂向力功率譜曲線變化相對平緩，但同一有砟軌道線路的輪軌力功率譜擬合模型在形式上相似，均可以用RFM模型進行擬合。

圖10　不同速度級輪軌力功率譜的擬合曲線

5　結　論

(1)有砟軌道輪軌力檢測數據的時域分析結果表明：時域波形具有相似特性，車速的變化并不改變其波形的形狀，僅會引起其峰值的變化；其在分布上符合正態(tài)分布，且滿足平穩(wěn)隨機過程的檢驗條件；頻域分析結果表明：有砟軌道線路上輪軌力的功率譜曲線包含周期性諧波成分和隨機成分，周期性諧波出現的頻段對應特定的空間波長且與車速有關，而輪軌力功率譜中的隨機成分呈現出相似的變化規(guī)律，可用擬合模型對其進行描述。

(2)根據基于平均周期圖譜分析Welch方法得到不同速度級下輪軌力功率譜的特征，提出了輪軌力功率譜曲線擬合模型(RFM)，主要用于描述輪軌力功率譜中隨機成分的分布規(guī)律；輪軌力功率譜中周期性諧波成分的波長主要對應車輪周長、橋梁單倍跨距等。在擬合模型中，首先采用五點三次平滑算法對輪軌力功率譜初始曲線進行平滑處理，然后采用非線性最小二乘擬合方法中的Levenberg-Marquardt算法進行擬合，最后采用擬合優(yōu)度進行檢驗。

(3)盡管不同的有砟軌道線路對應的輪軌力功率譜分布特征不一樣，但實測數據驗證了本文提出的RFM擬合模型結構簡單，擬合優(yōu)度高，普適性強，該模型能夠較好地描述有砟軌道線路輪軌力功率譜的分布狀態(tài)。

(4)本文采用的輪軌力檢測數據采樣頻率是500 Hz，輪軌力功率譜RFM擬合方法適用頻段范圍為1～250 Hz。由于高頻輪軌力功率譜的分布離散性大，為了提高擬合曲線在全頻段范圍內的擬合優(yōu)度，可以根據需要分頻段利用RFM模型進行曲線擬合。

[1]祖宏林，張志超，汪偉. 輪軌力測量在高速鐵路軌道檢測中的應用研究[J].鐵道機車車輛，2012,32(4):19-24.

(ZU Honglin,ZHANG Zhichao,WANG Wei. Application Research of the Wheel-Rail Interaction Force Measurement on the High-Speed Railway Track Inspection [J].Railway Locomotive & Car, 2012,32(4):19-24. in Chinese)

[2]馮青松，雷曉燕，練松良. 軌道隨機不平順影響下高速鐵路地基動力分析模型[J].振動工程學報，2013,26(6):927-934.

(FENG Qingsong,LEI Xiaoyan,LIAN Songliang. A Model of Dynamic Analysis of Ground for High-Speed Railway with Track Random Irregularities [J].Journal of Vibration Engineering,2013,26(6):927-934. in Chinese)

[3]SUN Yanquan, COLIN Cole, MAKSYM Spiryagin. Study on Track Dynamic Forces Due to Rail Short-Wavelength Dip Defects Using Rail Vehicle-Track Dynamics Simulations [J]. Journal of Mechanical Science & Technology, 2013, 27(3):629-640.

[4]BLANCO-Lorenzo J, SANTAMARIA J, VADILLO E G, et al. Dynamic Comparison of Different Types of Slab Track and Ballasted Track Using a Flexible Track Model [J]. Proceedings of the Institution of Mechanical Engineers Part F：Journal of Rail & Rapid Transit, 2011, 225(6):574-592.

[5]王建西，許玉德，練松良，等.隨機輪軌力作用下鋼軌滾動接觸疲勞裂紋萌生壽命預測仿真[J].鐵道學報,2010,32(3):66-70.

(WANG Jianxi,XU Yude, LIAN Songliang,et al. Simulation of Predicting RCF Crack Initiation Life Rails under Random Wheel-Rail Forces [J].Journal of the China Railway Society, 2010, 32(3):66-70. in Chinese)

[6]全國統計方法應用標準化技術委員會. GB/T 4882—2001數據的統計處理和解釋正態(tài)性檢驗[S].北京：中國標準出版社，2002.

(National Statistical Method Application Standardization Technical Committee .GB/T 4882—2001 Statistical Interpretation of Data-Normality Tests[S].Beijing: Standards Press of China, 2002. in Chinese)

[7]張樹京，齊立新.時間序列分析簡明教材[M].北京：清華大學出版社，2003.

[8]PETRE Stoica,RANDOLPJ Moses. Spectral Analysis of Signals[M].New Jersey：Pearson Prentice Hall,2005.

[9]張鴻燕，耿征.Levenberg-Marquardt算法的一種新解釋[J] .計算機工程與應用，2009，45(19):5-8.

(ZHANG Hongyan, GENG Zheng. Novel Interpretation for Levenberg-Marquardt Algorithm[J]. Computer Engineering and Application,2009,45(19):5-8. in Chinese)

[10]MADSEN K, NIELSEN H B, TINGLEFF O. Methods for Non-Linear Least Square Problems [M]. Lyngby:Informatics and Mathematical Modelling Technical University of Denmark, DTU, 2004.

[11]The Mathworks,Inc. Curve Fitting Toolbox User’s Guide[M].Natick MA：The MathWorks Inc.2010.