基于WPD-ASTFT和SVM重載鐵路鋼軌波磨診斷方法

肖炳環(huán)，劉金朝，牛留斌，邵 奇，羅澤霖，陳仕明

(中國鐵道科學(xué)研究院集團有限公司 基礎(chǔ)設(shè)施檢測研究所，北京 100081)

波磨指鋼軌表面不均勻的波浪形磨耗，是重載鐵路中常見的軌道病害。鋼軌波磨會加劇輪軌間相互作用，加上重載列車本身的軸重大，異常的輪軌力會破壞車輛和軌道部件，造成疲勞傷損，嚴重時會影響行車安全。目前尚無有效消除鋼軌波磨的技術(shù)手段，相比通過鋼軌潤滑調(diào)解輪軌摩擦因數(shù)、采用軌道吸振器等措施，打磨被認為是抑制波磨過快發(fā)展的有效措施[1]。利用波磨小車人工上道查找鋼軌波磨是最直接的辦法，然而，上道作業(yè)受到天窗時間限制，同時受作業(yè)人員操作水平等因素影響，會造成人工查找鋼軌波磨的作業(yè)效率低下。如何及時高效地發(fā)現(xiàn)并評估鋼軌波磨成為國內(nèi)外眾多學(xué)者研究的重要課題。

1999年，文獻[2]采用簧下質(zhì)量加速度信號的方法，通過研究波蘭的鐵路信號總結(jié)并分析不同波磨區(qū)段頻譜；2009年，文獻[3] 提出基于時頻分析技術(shù)對鋼軌波磨進行診斷，利用ARCAP方法估計波磨的波深和頻率；2001年，文獻[4]提出了一種基于鋼軌表面紋理分析的技術(shù)，用于鋼軌波磨的檢測和分類；2015年，文獻[5]采用軸箱加速度測量系統(tǒng)來檢測波磨，并提出了一種基于識別出的特征曲線的自動檢測算法；2017年，文獻[6]提出了一種基于激光攝像技術(shù)的波磨檢測方法，通過多傳感器并行工作對鋼軌輪廓進行高頻率采樣，并將圖像的ROI傳送至上位機提取輪廓數(shù)據(jù)，進而完成鋼軌波磨信息的提?。?020年，文獻[7]提出一種基于參數(shù)優(yōu)化變分模態(tài)分解(Variable mode decomposition, VMD)和平滑偽維格納分布(Smoothed pseudo Wigner-Ville distribution, SPWVD)的軌道波磨辨識方法，信號分解后以包絡(luò)熵為指標(biāo)，采用SPWVD方法對分解后的信號進行時頻分析，確定波磨發(fā)生的位置及波長。

利用圖像檢測方法診斷鋼軌波磨時，受外部環(huán)境因素影響較大，容易造成漏判；利用時頻分析方法判斷數(shù)據(jù)特性雖然可行，但是由于時頻分布的分辨率問題經(jīng)常導(dǎo)致結(jié)果準(zhǔn)確率不高。因此本文提出了自適應(yīng)時頻分析和機器學(xué)習(xí)結(jié)合的方法，通過重載鐵路綜合檢測車采集的軸箱加速度信號診斷鋼軌波磨。先利用小波包分解(Wavelet packet decomposition, WPD)的方法將信號分解為若干子信號，然后對子信號進行自適應(yīng)短時傅里葉變換(Adaptive short-time Fourier transform, ASTFT)獲得信號高分辨率的時頻分布，并提取每個時頻分布的熵值作為分類特征，再對特征數(shù)據(jù)進行降維，最后將所有數(shù)據(jù)投入支持向量機(Support vector machine，SVM)進行訓(xùn)練和分類。上述基于WPD-ASTFT和SVM的診斷重載鋼軌波磨方法通過自適應(yīng)時頻分析可以提高時頻分布分辨率，從而提高診斷結(jié)果準(zhǔn)確率。

1 基礎(chǔ)理論

小波包分解能對原始信號進行全頻段的細致分解，即在低頻段和高頻段繼續(xù)進行分解，可改進小波分解導(dǎo)致的分辨率不能同時兼顧高低頻段的問題，從而可以使原始信號的頻率分辨率更佳[8]。原始信號s(t)經(jīng)過n層小波包分解，得到N=2n個子信號，且：

s(t)=s1(t)+s2(t)+…+sN(t)

(1)

式中：t——時間；

si——分解后的第i個子信號，i=1,2,...,N。

對數(shù)據(jù)采用SVM分類，如圖1所示。

圖1 SVM數(shù)據(jù)分類原理

圖1中，兩類樣本點分布在紅色線的兩側(cè)，該紅色線所在平面稱為最大間隔超平面[9]，可以表示為：

ωTx+b=0

(2)

式中：ω——平面法向量；

b——支持向量到平面的距離；

x——空間上一點，x=(x1,x2,...,xn)。

樣本點到紅色線的距離要最大化，空間上一點x到ωTx+b=0的距離為l，表示為：

(3)

對數(shù)據(jù)進行分類時的分類決策函數(shù)h(x)為：

h(x)=sign(ωTx+b)

(4)

sign(x)為sigmoid函數(shù)。將預(yù)測樣本帶入分類決策函數(shù)中即可得到分類結(jié)果[10]。

在已有的基礎(chǔ)理論上提出WPD-ASTFT波磨診斷方法，含有多頻率成分的信號在進行短時傅里葉變換時，在一個時刻只能選擇某一固定窗長，這就使得信號的時頻分布分辨率不足，因此，先利用WPD把信號中主頻相差較大的頻率分解到不同的子信號，這樣不同的子信號在同一時刻可以選擇不同的窗長進行短時傅里葉變換，然后對不同的子信號進行自適應(yīng)短時傅里葉變換，通過計算獲得使信號時頻分布聚集性最好的窗長，稱為最優(yōu)窗長。利用Renyi熵刻畫時頻分布的聚集性，熵值越小時頻聚集性越好。同時窗長與Renyi熵值對應(yīng)關(guān)系作為粒子群優(yōu)化算法(Particle swarm optimization，PSO)求最優(yōu)窗長時的適應(yīng)度函數(shù)。

2 自適應(yīng)短時傅里葉變換和特征提取

2.1 自適應(yīng)短時傅里葉變換

短時傅里葉變換(Short time Fourier transform, STFT)作為傳統(tǒng)時頻分析工具，其窗長大小對時頻分布的分辨率影響較大，為了避免人工調(diào)試窗長并且保證信號做短時傅里葉變換后具有較好分辨率的時頻分布，結(jié)合Renyi熵和粒子群優(yōu)化算法提出了自適應(yīng)短時傅里葉變換。對每個子信號進行自適應(yīng)短時傅里葉變換，在對應(yīng)的最優(yōu)窗長下計算時頻分布的Renyi熵[11]。

設(shè)某信號包含n個信息y1,y2,y3,…,yn，每一個信息出現(xiàn)的概率是p(y1)，p(y2)，p(y3),…,p(yn)，由它們組成的一個系統(tǒng)S為：

(5)

則這個系統(tǒng)的Renyi熵H(y)為：

(6)

其中，q>0且q≠1。

本文用Renyi熵值作為評判標(biāo)準(zhǔn)，評價時頻表示集中度。為使信號做短時傅里葉變換后時頻分辨率最佳，需要獲取該子信號做短時傅里葉變換時的最佳窗長。因此計算信號在不同窗長下做短時傅里葉變換得到的時頻分布STFT(t,f)的Renyi熵Eζ，計算公式如下：

(7)

式中：l——計算系數(shù)，為大于0的常數(shù)；

α——時間；

β——頻率；

T——信號持續(xù)時間。

粒子群優(yōu)化算法相比于遺傳算法簡單易行，收斂速度快，參數(shù)調(diào)整方便[12]。若D維空間中有P個粒子，粒子i的位置xi=(xi1,xi2,...,xiD),將xi代入適應(yīng)函數(shù)f(x)求適應(yīng)值。粒子i的速度vi=(vi1,vi2,...,viD)。粒子i的第d維速度公式更新為：

(8)

式中：pbestid——粒子i個體經(jīng)歷過的最好位置；

gbestd——種群所經(jīng)歷過的最好的位置；

w——慣性權(quán)重，非負數(shù)，調(diào)節(jié)對解空間的搜索范圍；

c1、c2——均為加速常數(shù)，調(diào)節(jié)學(xué)習(xí)最大步長，通常c1=c2=2；

r1、r2——均為隨機函數(shù)，取值范圍為[0,1]，以增加搜索隨機性；

k——迭代次數(shù)。

粒子i的第d維位置公式更新為：

(9)

子信號的不同窗長參數(shù)設(shè)為粒子，做短時傅里葉變換后得到的時頻分布的Renyi熵值為適應(yīng)度函數(shù)，計算適應(yīng)度函數(shù)最小時粒子的位置和Renyi熵值。信號經(jīng)過小波包分解后得到N個子信號。每個子信號主頻有所差異，經(jīng)過粒子群優(yōu)化過程確定做短時傅里葉變換的最優(yōu)窗長，對每個子信號si(t)分別和高斯窗函數(shù)g做短時傅里葉變換[13]：

(10)

式中：g(ξ-t)——最優(yōu)窗函數(shù)；

f——頻率。

經(jīng)過自適應(yīng)短時傅里葉變換后，N個子信號在對應(yīng)的最優(yōu)窗長下獲得Renyi熵值最小的時頻分布，同時，在變換過程中Renyi熵值既可以作為分辨率評判標(biāo)準(zhǔn)，又可以在分類過程中作為數(shù)據(jù)的特征。

2.2 特征指標(biāo)提取及其降維方法

為了提取兩類數(shù)據(jù)Z列熵值差異較大的值作為數(shù)據(jù)特征，并實現(xiàn)降維，提出均值特征降維法，根據(jù)兩類數(shù)據(jù)Z列熵值的平均值差異對數(shù)據(jù)進行降維處理，具體計算步驟如下：

(1) 計算平均值。

計算正常軌道區(qū)段加速度Z列熵值每一列的平均值mi：

(11)

式中：J——正常軌道區(qū)段個數(shù)；

Oj×i——正常軌道熵值矩陣。

(12)

式中：K——波磨軌道區(qū)段個數(shù)；

Ck×i——波磨軌道熵值矩陣。

(2) 計算熵值差值。

正常和含有波磨區(qū)段Z列平均值對應(yīng)相減并取其絕對值：

(13)

(14)

3 WPD-ASTFT和SVM診斷步驟

基于WPD-ASTFT和SVM的重載鐵路鋼軌波磨診斷步驟如下：

(1) 將重載鐵路軸箱垂向加速度信號按50 m劃分為單元；

(2) 各單元信號進行小波包分解，得到若干個子信號；

(3) 每個單元各個子信號進行自適應(yīng)短時傅里葉變換，并計算時頻分布對應(yīng)的Renyi熵值；

(4) 利用均值特征降維法對熵值數(shù)據(jù)降維；

(5) 投入SVM分類器訓(xùn)練和預(yù)測。

4 試驗驗證

試驗采用國內(nèi)某重載鐵路綜合檢測列車采集的軸箱垂向加速度信號。將軸箱加速度信號劃分為50 m一個單元，一個單元軸箱加速度信號經(jīng)過小波包分解被分解成若干個頻帶范圍不同的子信號，并且原始信號可由子信號重構(gòu)：s(t)=s1(t)+…+sN(t)，其中N為子信號個數(shù)。本文通過三層小波包分解得到8個子信號，并將第6、7、8個子信號相加作為最終的第6子信號，這樣既可以節(jié)約計算時間，又可以充分分解軸箱加速度信號。

4.1 信號自適應(yīng)時頻特性分析和特征提取

分別選取一段重載鐵路正常軌道區(qū)段和波磨軌道區(qū)段的軸箱加速度信號，如圖2所示，正常軌道區(qū)段里程為K401+050～K401+100，波磨區(qū)段里程為K387+150～K387+200。對2段信號進行三層小波包分解后得到6個子信號，按照WPD-ASTFT步驟計算子信號最優(yōu)窗長，最后求得時頻分布。

圖2 正常和波磨軌道區(qū)段軸箱加速度

對正常軌道區(qū)段和波磨區(qū)段加速度子信號進行自適應(yīng)短時傅里葉變換，結(jié)果如圖3所示。由圖3可以看出，波磨區(qū)段子信號1對應(yīng)的時頻分布中由波磨造成的軸箱加速度的頻率集中性明顯高于正常軌道區(qū)段。計算正常和波磨區(qū)段WPD-ASTFT后時頻分布對應(yīng)的Renyi熵值，計算結(jié)果如表1所示。

表1 WPD-ASTFT后時頻分布對應(yīng)的Renyi熵值

計算所有區(qū)段經(jīng)過WPD-ASTFT后時頻分布對應(yīng)的Renyi熵值，表2為計算結(jié)果，其中1～75為正常區(qū)段軸箱加速度信號的Renyi熵值，76～100為波磨區(qū)段對應(yīng)的Renyi熵值。

表2 計算所得Renyi熵值

圖3 正常和波磨區(qū)段子信號對應(yīng)時頻分布

表3 正常區(qū)段和波磨區(qū)段Renyi熵值的均值和差值

4.2 分類結(jié)果

將均值特征降維得到的第一特征熵值和第二特征熵值放入SVM分類器中進行訓(xùn)練。表2中序號1～75為正常區(qū)段，標(biāo)簽為‘1’，序號76～100為波磨區(qū)段，標(biāo)簽為‘-1’。序號1～100為訓(xùn)練樣本，將其投入SVM分類器中進行訓(xùn)練，結(jié)果如圖4所示，可以看出兩類數(shù)據(jù)被超平面完美分隔。再把測試數(shù)據(jù)放入分類器進行預(yù)測，測試數(shù)據(jù)同樣保留第一特征熵值和第二特征熵值，所有正常區(qū)段和波磨區(qū)段的訓(xùn)練測試數(shù)據(jù)及分類結(jié)果如圖5所示，可以看出15個區(qū)段準(zhǔn)確預(yù)測14個，準(zhǔn)確率達93.33%，準(zhǔn)確率較高。

圖4 SVM訓(xùn)練結(jié)果

圖5 SVM預(yù)測結(jié)果

5 結(jié)束語

本文結(jié)合時頻分析方法和機器學(xué)習(xí)，提出了基于WPD-ASTFT和SVM的重載鐵路鋼軌波磨診斷方法。將該方法應(yīng)用于診斷重載鐵路鋼軌波磨，通過預(yù)測結(jié)果和實際情況對比，驗證了方法的可靠性，預(yù)測準(zhǔn)確率高達93.33%，可以有效指導(dǎo)養(yǎng)護維修工作，抑制波磨的發(fā)展。