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

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

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

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

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

1 基礎(chǔ)理論

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

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

(1)

式中：t——時(shí)間；

si——分解后的第i個(gè)子信號(hào)，i=1,2,...,N。

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

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

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

ωTx+b=0

(2)

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

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

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

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

(3)

對(duì)數(shù)據(jù)進(jìn)行分類(lèi)時(shí)的分類(lèi)決策函數(shù)h(x)為：

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

(4)

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

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

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

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

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

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

(5)

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

(6)

其中，q>0且q≠1。

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

(7)

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

α——時(shí)間；

β——頻率；

T——信號(hào)持續(xù)時(shí)間。

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

(8)

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

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

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

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

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

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

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

(9)

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

(10)

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

f——頻率。

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

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

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

(1) 計(jì)算平均值。

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

(11)

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

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

(12)

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

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

(2) 計(jì)算熵值差值。

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

(13)

(14)

3 WPD-ASTFT和SVM診斷步驟

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

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

(2) 各單元信號(hào)進(jìn)行小波包分解，得到若干個(gè)子信號(hào)；

(3) 每個(gè)單元各個(gè)子信號(hào)進(jìn)行自適應(yīng)短時(shí)傅里葉變換，并計(jì)算時(shí)頻分布對(duì)應(yīng)的Renyi熵值；

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

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

4 試驗(yàn)驗(yàn)證

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

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

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

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

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

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

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

表2 計(jì)算所得Renyi熵值

圖3 正常和波磨區(qū)段子信號(hào)對(duì)應(yīng)時(shí)頻分布

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

4.2 分類(lèi)結(jié)果

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

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

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

5 結(jié)束語(yǔ)

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