適用復(fù)雜劣化趨勢的軌道不平順魯棒建模方法

楊雅琴，徐 鵬*，李 曄，孫全欣

(1.北京交通大學(xué)交通運(yùn)輸學(xué)院，北京100044；2.中國鐵路南昌局集團(tuán)有限公司，南昌330002)

0 引 言

各國鐵路均利用運(yùn)營線路上鋼軌在三維空間中的投影與其理想位置的偏差量化軌道不平順狀態(tài)，包括軌距、水平、多種弦長的高低和軌向等參數(shù)，通過峰值和軌道質(zhì)量指數(shù)(Track Quality Index，TQI)評(píng)價(jià)軌道局部和整體不平順狀態(tài)[1].為保證行車安全、平穩(wěn)和舒適，鐵路工務(wù)部門開展軌道維修作業(yè)，將軌道狀態(tài)控制在允許范圍內(nèi).掌握軌道不平順劣化規(guī)律，預(yù)測其未來發(fā)展?fàn)顟B(tài)，是合理安排維修作業(yè)計(jì)劃，實(shí)現(xiàn)鐵路線路預(yù)測性維修的基礎(chǔ)，對(duì)鐵路智能化運(yùn)營管理具有重要意義.國內(nèi)外學(xué)者圍繞預(yù)測軌道不平順劣化開展了大量研究.Andrade 等[2]通過條件自回歸(CAR)建立了維修作業(yè)后高低值和劣化速率函數(shù)，據(jù)此建立了軌道高低不平順分層貝葉斯劣化模型；曲建軍等[3]將灰色GM(1,1)與基于殘差修正的改進(jìn)馬爾可夫鏈相結(jié)合，建立改進(jìn)灰色—馬爾可夫鏈組合預(yù)測模型描述兩次維修作業(yè)之間的軌道不平順劣化過程.

因維修作業(yè)明顯改變軌道不平順狀態(tài)，學(xué)者以單一函數(shù)或者組合函數(shù)(線性、指數(shù)等)建立兩次維修作業(yè)間的軌道不平順劣化模型.通過假設(shè)維修作業(yè)前軌道不平順的初始狀態(tài)和劣化率服從特定函數(shù)，將軌道不平順劣化模型的時(shí)間域拓展到多次維修之間.上述假設(shè)在路基支撐剛度一致和自然氣候條件變化不大的區(qū)段是成立的.但我國地域廣闊，南北溫差和降雨差異較大，南方和西南地區(qū)土壤中含有膨脹土，軌道劣化復(fù)雜多變；另外，維修作業(yè)(如大機(jī)搗固、換軌等)會(huì)改變軌道剛度.所以，上述假設(shè)在實(shí)際運(yùn)營線路的大部分區(qū)段是很難滿足的.因此，根據(jù)維修作業(yè)在時(shí)間維度上對(duì)軌道不平順劣化過程進(jìn)行分段是深入研究軌道劣化規(guī)律的前提.中國國家鐵路集團(tuán)公司(簡稱國鐵集團(tuán))從2016年在全國各鐵路局集團(tuán)公司推廣應(yīng)用《鐵路工務(wù)安全生產(chǎn)管理信息系統(tǒng)》，積累了近3～4年的維修作業(yè)地點(diǎn)和日期信息.我國鐵路上，檢測列車運(yùn)行條件下軌道不平順設(shè)備包括高速鐵路上的綜合檢測列車和綜合巡檢車，以及普速鐵路上的軌道檢測車(以下統(tǒng)稱為動(dòng)軌檢車)，將其檢測數(shù)據(jù)統(tǒng)稱為動(dòng)態(tài)檢測數(shù)據(jù).動(dòng)態(tài)檢測數(shù)據(jù)是計(jì)算TQI 的數(shù)據(jù)基礎(chǔ)，也是鐵路部門感知軌道狀態(tài)的重要依據(jù).動(dòng)軌檢車在檢測過程中，不可避免地受到環(huán)境或檢測系統(tǒng)故障影響，產(chǎn)生異常檢測數(shù)據(jù)，使計(jì)算得出的TQI異常，這將對(duì)軌道不平順劣化過程中維修作業(yè)的判別造成干擾.

綜上，本文基于最小描述長度準(zhǔn)則探索建立一套動(dòng)態(tài)檢測數(shù)據(jù)驅(qū)動(dòng)的軌道不平順劣化自適應(yīng)分段建模方法(Minimum-Description-Length-Based Rail Track Deterioration Adaptive Segmentation Framework，MDL-RTDAS). MDLRTDAS 能夠克服檢測數(shù)據(jù)異常影響，靈敏感知劣化趨勢的變化，自動(dòng)識(shí)別軌道不平順劣化過程中維修作業(yè)發(fā)生的日期.本文將為正在開展的“通過挖掘高鐵歷史檢測數(shù)據(jù)，掌握高鐵開通十年來劣化規(guī)律，指導(dǎo)高鐵預(yù)測性維修”奠定可靠的數(shù)據(jù)挖掘基礎(chǔ).

1 軌道不平順劣化自適應(yīng)分段建模方法

1.1 最小描述長度準(zhǔn)則

將軌道劣化過程中維修作業(yè)發(fā)生的日期稱為維修作業(yè)點(diǎn).軌道不平順劣化趨勢會(huì)在維修作業(yè)點(diǎn)處發(fā)生突變.因此，維修作業(yè)點(diǎn)的識(shí)別問題類似于應(yīng)用統(tǒng)計(jì)學(xué)中時(shí)間序列突變點(diǎn)的識(shí)別.基于最小描述長度(Minimum Description Length，MDL)準(zhǔn)則的突變點(diǎn)識(shí)別方法是一種不需要突變點(diǎn)個(gè)數(shù)和位置等先驗(yàn)知識(shí)，依據(jù)時(shí)間序列特征進(jìn)行識(shí)別的方法.Davis[4]基于MDL準(zhǔn)則和AR模型，建立了一套自動(dòng)識(shí)別非平穩(wěn)時(shí)間序列突變點(diǎn)的模型.LU Q.Q.等[5]利用MDL準(zhǔn)則估計(jì)具有自相關(guān)性和周期性的天氣時(shí)間序列中的突變點(diǎn).MDL 準(zhǔn)則是Jorma Rissanen 提出的最優(yōu)編碼理論[6]，基本原理是采用某種模型對(duì)數(shù)據(jù)進(jìn)行編碼，用編碼后的數(shù)據(jù)和模型表示原始數(shù)據(jù)，編碼后占用存儲(chǔ)空間(稱為總描述長度)最小的編碼模型，即為最優(yōu)編碼模型.基于MDL準(zhǔn)則識(shí)別時(shí)間序列突變點(diǎn)的核心思想是將識(shí)別問題轉(zhuǎn)化為模型選擇問題，即通過搜索能夠最完整保留原始數(shù)據(jù)包含信息同時(shí)最大程度將其壓縮的最優(yōu)擬合模型，估計(jì)突變點(diǎn)的個(gè)數(shù)和位置.

依據(jù)MDL 準(zhǔn)則，總描述長度分為兩部分：模型的編碼長度和對(duì)應(yīng)殘差的編碼長度.對(duì)于長度為l的時(shí)間序列y={y1,y2,…,yl} ，其中，y1,…,yl為時(shí)間序列y中的元素；選擇模型擬合y，令表示y的擬合值，則對(duì)應(yīng)殘差為以Ly表示y的編碼長度，則Ly可分解為

1.2 MDL-RTDAS

MDL-RTDAS 通過識(shí)別軌道劣化過程中維修作業(yè)導(dǎo)致的劣化狀態(tài)突變，自動(dòng)將劣化過程分段，使用線性回歸模型擬合每一段劣化過程，據(jù)此，建立軌道不平順劣化自適應(yīng)分段模型.這種方法不依賴維修作業(yè)記錄，不假設(shè)檢測數(shù)據(jù)有效性，完全由動(dòng)態(tài)檢測數(shù)據(jù)驅(qū)動(dòng)，能夠適應(yīng)各種自然條件下的軌道劣化建模.

根據(jù)式(1)推導(dǎo)線性回歸模型的編碼長度和殘差編碼長度的具體過程如下.

(1)模型編碼長度.

軌道高低不平順VTQI-L是安排維修作業(yè)的主要狀態(tài)指標(biāo)，故本文主要圍繞軌道高低劣化進(jìn)行研究.采用一元線性回歸模型擬合每一段子變化過程，設(shè)某一時(shí)間段內(nèi)進(jìn)行了n次動(dòng)態(tài)檢測，每次檢測日期為t1,t2,…,tn，記其集合為t={t1,t2,…,tn} ，某200 m 線路單元在各次檢測的VTQI-L值為q1,q2,…,qn，記其集合為q={q1,q2,…,qn} ，期間進(jìn)行過m次維修作業(yè)(0＜m ＜n)，第j個(gè)和第j+1個(gè)子過程之間的維修作業(yè)日期為τj，令τ0=t1，τm+1=tn，則第j個(gè)子過程的擬合模型表示為

式中：βj={βj0,βj1} 為第j個(gè)子過程擬合模型的參數(shù)集合，其中，βj0為截距，βj1為斜率；ε為隨機(jī)項(xiàng)，令μe為ε的均值，為ε的方差，則可見，的待估計(jì)參數(shù)包括τj和βj.令表示所有子過程擬合模型的集合，則表示為

令nj=τj-τj-1代表第j個(gè)子過程包含的數(shù)據(jù)點(diǎn)數(shù)，則(τ1,…,τm)與(n1,…,nm+1)蘊(yùn)含的信息是等價(jià)的，式(3)可轉(zhuǎn)換為

依據(jù)MDL準(zhǔn)則：對(duì)于有限整數(shù)I，若I的上限未知，則LI=lbI；若I的上限已知，設(shè)為IU，則LI≈lbIU.通過極大似然估計(jì)擬合N個(gè)觀測值，得到的一元線性回歸模型參數(shù)的編碼長度為由此，可得

(2)模型殘差編碼長度.

對(duì)于Le，MDL 準(zhǔn)則中給出定義：擬合模型的殘差e的編碼長度Le為的負(fù)對(duì)數(shù)似然估計(jì)(以2為底)[6].令表示q的估計(jì)值，則中各子過程擬合模型相互獨(dú)立，可得

將式(8)和式(9)相加，得到軌道不平順劣化分段擬合模型的編碼長度，記為Lq.為簡化運(yùn)算，將對(duì)數(shù)底由2換為e，得到

式(10)為本文解決軌道不平順劣化過程中維修作業(yè)點(diǎn)識(shí)別問題的目標(biāo)方程式，使式(10)取得最小值的一組模型，即為軌道不平順劣化的最優(yōu)分段擬合模型，其對(duì)應(yīng)的分段點(diǎn)則為要識(shí)別的維修作業(yè)點(diǎn).

2 模型求解

在突變點(diǎn)個(gè)數(shù)未知的情況下，求解式(10)最小值的計(jì)算量很大.本文數(shù)據(jù)VTQI-L具有以下兩個(gè)特征：①樣本量小，根據(jù)《普速鐵路線路修理規(guī)則》對(duì)動(dòng)態(tài)檢測周期的規(guī)定及現(xiàn)場調(diào)研，主要繁忙干線1個(gè)月有2次檢測[1]，故研究對(duì)象的樣本量較??；②檢測異常值干擾較大，異常檢測數(shù)據(jù)的存在對(duì)識(shí)別維修作業(yè)點(diǎn)造成干擾，如何克服檢測數(shù)據(jù)異常的干擾是本文設(shè)計(jì)模型求解算法要考慮的關(guān)鍵問題之一.結(jié)合研究數(shù)據(jù)特征，克服異常檢測值的干擾，同時(shí)減小計(jì)算復(fù)雜度，模型求解算法具體步驟如下：

Step 1確定疑似維修作業(yè)點(diǎn).

軌道高低不平順劣化過程中異常點(diǎn)和維修作業(yè)點(diǎn)的特征明顯：異常點(diǎn)明顯偏離當(dāng)前的劣化趨勢，進(jìn)行一次維修作業(yè)后，高低不平順值會(huì)降低到規(guī)定范圍內(nèi)，進(jìn)入下一個(gè)劣化過程.如圖1所示.

圖1 軌道高低不平順劣化過程中的異常檢測數(shù)據(jù)及維修作業(yè)Fig.1 Maintenance operations and anomalous measurements in track longitudinal level irregularity deterioration

將異常點(diǎn)和維修作業(yè)點(diǎn)統(tǒng)稱為疑似維修作業(yè)點(diǎn)，令di=qi-qi-1表示VTQI-L的一階差分，記其集合為d，則疑似維修作業(yè)點(diǎn)的一階差分會(huì)明顯偏高/偏低.設(shè)ωi為衡量某點(diǎn)一階差分值與相鄰點(diǎn)的差異程度，則有

設(shè)一階差分d的上四分位數(shù)為W0，當(dāng)ωi ＞W(wǎng)0時(shí)，認(rèn)為該點(diǎn)對(duì)應(yīng)的檢測日期ti為疑似維修作業(yè)點(diǎn).

Step 2計(jì)算不同分段方式的MDL值.

確定疑似維修作業(yè)點(diǎn)后，設(shè)其個(gè)數(shù)為.假設(shè)從個(gè)疑似維修作業(yè)點(diǎn)中選取r個(gè)，利用動(dòng)態(tài)規(guī)劃方法求得令MDL值最小的r個(gè)疑似維修作業(yè)點(diǎn)的組合方式，記其最小值為Mr|min.考慮到1 個(gè)月內(nèi)同一個(gè)線路單元不會(huì)安排兩次維修作業(yè)[1]，故在組合疑似維修作業(yè)點(diǎn)時(shí)，若存在相鄰兩點(diǎn)的時(shí)間差小于1個(gè)月，則認(rèn)為該種組合方式無效.

Step 3選擇最優(yōu)分段方式.

根據(jù)MDL 準(zhǔn)則，令MDL 最小的為最優(yōu)擬合模型.但實(shí)際分析時(shí)發(fā)現(xiàn)，對(duì)于本文研究對(duì)象，當(dāng)MDL 值最小時(shí)，往往將軌道不平順劣化過程過度分段，不利于劣化過程分析.分析Mr|min隨被選疑似維修作業(yè)點(diǎn)的個(gè)數(shù)r增加的變化過程，認(rèn)為當(dāng)|Mr+1|min-Mr|min|小于閾值δ時(shí)，Mr|min對(duì)應(yīng)的疑似維修作業(yè)點(diǎn)組合方式為最優(yōu).最優(yōu)分段擬合模型的分段點(diǎn)，即為本文要識(shí)別的維修作業(yè)點(diǎn)，未被選中的疑似維修作業(yè)點(diǎn)認(rèn)為是異常點(diǎn).在本文中，δ=10.

3 計(jì)算驗(yàn)證

數(shù)據(jù)來源于2014—2019年昌福高速鐵路下行K21+184～K220+308的原始動(dòng)態(tài)檢測數(shù)據(jù).原始動(dòng)態(tài)檢測數(shù)據(jù)的里程偏差修正和波形對(duì)齊，各項(xiàng)TQI計(jì)算等一系列數(shù)據(jù)預(yù)處理過程依托團(tuán)隊(duì)研發(fā)的《基于動(dòng)態(tài)檢測數(shù)據(jù)的軌道變形分析系統(tǒng)》完成[7].MDL-RTDAS 的關(guān)鍵是劣化過程中維修作業(yè)點(diǎn)的識(shí)別，選取多個(gè)200 m 線路單元區(qū)段，針對(duì)這些區(qū)段劣化過程中的維修作業(yè)點(diǎn)，對(duì)比分析MDLRTDAS 的識(shí)別結(jié)果和人工分析確認(rèn)的結(jié)果，同時(shí)展示Davis等的Auto-PARM[4]在這些區(qū)段上的識(shí)別結(jié)果，據(jù)此來驗(yàn)證MDL-RTDAS的有效性.

3.1 案例分析

選取3 個(gè)區(qū)段，直觀對(duì)比MDL-RTDAS 與Auto-PARM 的識(shí)別結(jié)果及分段擬合模型，分別計(jì)算兩個(gè)模型擬合殘差的均值(EMR)、標(biāo)準(zhǔn)差(ER-SD)和平均絕對(duì)值(EMAR)，對(duì)比分析擬合效果.圖2(a)、圖3(a)和圖4(a)為本文MDL-RTDAS 的識(shí)別結(jié)果及分段擬合模型，圖2(b)、圖3(b)和圖4(b)為Auto-PARM 的結(jié)果.圖中，實(shí)線代表該區(qū)段的實(shí)際高低不平順值，虛線代表模型計(jì)算出的擬合值，垂直線代表模型識(shí)別出的維修作業(yè)點(diǎn).

(1)區(qū)段1(K84+601～K84+801).

該區(qū)段位于曲線上，根據(jù)其軌道高低不平順劣化過程分析出該區(qū)段于2015-07-24、2017-12-26進(jìn)行了維修作業(yè).MDL-RTDAS識(shí)別出的維修作業(yè)點(diǎn)為2015-07-24、2017-12-26，如圖2(a)所示；Auto-PARM 識(shí)別出的維修作業(yè)點(diǎn)為2015-07-24、2016-01-12、2016-08-12、2017-01-11、2017-12-26、2019-03-21，如圖2(b)所示.通過兩個(gè)模型擬合的殘差計(jì)算得到的EMR、ER-SD、EMAR如表1 所示，可得MDL-RTDAS擬合殘差的均值更接近于0，離散程度更小，平均絕對(duì)值更低.

圖2 K84+601～K84+801 (位于曲線上)Fig.2 K84+601～K84+801(on a curved track)

表1 兩個(gè)模型的殘差對(duì)比分析(區(qū)段1)Table 1 Residual comparison analysis of section one

(2)區(qū)段2(K103+334～K103+534).

該區(qū)段位于直線上，根據(jù)其軌道高低不平順劣化過程分析出該區(qū)段于2016-06-26、2016-11-26、2018-09-11 進(jìn)行了維修作業(yè).MDL-RTDAS 和Auto-PARM 識(shí)別出的維修作業(yè)點(diǎn)均為2016-06-26、2016-11-26、2018-09-11，如圖3(a)、(b)所示.通過殘差計(jì)算得到的EMR、ER-SD、EMAR如表2 所示，可得MDL-RTDAS擬合殘差的均值更接近于0，離散程度更小，平均絕對(duì)值更低.

圖3 K103+334～K103+534 (位于直線上)Fig.3 K103+334～K103+534(on a tangent track)

表2 兩個(gè)模型的殘差對(duì)比分析(區(qū)段2)Table 2 Residual comparison analysis of section two

(3)區(qū)段3(K106+551～K106+751).

該區(qū)段位于直線上，根據(jù)其軌道高低不平順劣化過程分析出該區(qū)段于2016-07-14、2017-09-12、2018-01-26 進(jìn)行了維修作業(yè).MDL-RTDAS 識(shí)別出的維修作業(yè)點(diǎn)為2016-07-14、2017-09-12、2018-01-26，如圖4(a)所示；Auto-PARM 識(shí)別出的維修作業(yè)點(diǎn)為2015-04-13、2016-07-14、2018-01-26，如圖4(b)所示.通過殘差計(jì)算得到的EMR、ER-SD、EMAR如表3所示，可得MDL-RTDAS擬合殘差的均值更接近于0，離散程度更小，平均絕對(duì)值更低.

圖4 K106+551～K106+751 (位于直線上)Fig.4 K106+551～K106+751(on a tangent track)

表3 兩個(gè)模型的殘差對(duì)比分析(區(qū)段3)Table 3 Residual comparison analysis of section three

3.2 綜合對(duì)比

為進(jìn)一步對(duì)比分析，另選若干區(qū)段統(tǒng)計(jì)分析每個(gè)區(qū)段劣化過程中異常檢測的個(gè)數(shù)MF,兩個(gè)模型擬合的殘差，以及識(shí)別結(jié)果的準(zhǔn)確率(PARC)、精確率(PPRC)，結(jié)果如表4所示.

式中：NT為模型正確識(shí)別出的維修作業(yè)點(diǎn)個(gè)數(shù)；NF為模型中錯(cuò)誤識(shí)別的維修作業(yè)點(diǎn)個(gè)數(shù)；Mmanual為人工分析出的維修作業(yè)點(diǎn)個(gè)數(shù).

表4 兩個(gè)模型的識(shí)別準(zhǔn)確度及擬合效果Table 4 Accuracy,precision and fitness of two algorithms

基于兩個(gè)模型識(shí)別的準(zhǔn)確率和精確率，令Fscore表示衡量模型的識(shí)別準(zhǔn)確度，則Fscore=0.5×PARC+0.5×PPRC.根據(jù)表4計(jì)算得到：兩個(gè)模型在各區(qū)段的Fscore，對(duì)比結(jié)果如圖5(a)所示；兩個(gè)模型擬合殘差計(jì)算的EMAR，對(duì)比結(jié)果如圖5(b)所示；MF對(duì)PARC的影響，如圖6所示.

從表4 中EMR和ER-SD可得，MDL-RTDAS 建立的分段擬合模型殘差的均值更接近于0 且離散程度更小；從圖5(a)可得，MDL-RTDAS 的識(shí)別準(zhǔn)確度高于Auto-PARM；從圖5(b)可得，MDLRTDAS 擬合殘差的平均絕對(duì)值低于Auto-PARM；從圖6 可得，隨著區(qū)段中異常檢測數(shù)據(jù)的增加，MDL-RTDAS 的識(shí)別準(zhǔn)確率沒有隨之變大/變小，說明MDL-RTDAS 可克服檢測數(shù)據(jù)異常的影響.綜上，MDL-RTDAS更適用于軌道高低不平順劣化自適應(yīng)分段建模，且具有一定的魯棒性.

圖5 兩個(gè)模型識(shí)別準(zhǔn)確度及平均絕對(duì)擬合誤差對(duì)比情況Fig.5 Comparison between two algorithms'performances

圖6 異常檢測數(shù)據(jù)對(duì)MDL-RTDAS 識(shí)別準(zhǔn)確率的影響Fig.6 Effect of contaminated measurements on PARC

4 結(jié) 論

本文提出將維修作業(yè)導(dǎo)致軌道不平順劣化過程突變的識(shí)別問題轉(zhuǎn)化為模型選擇問題，并基于最小描述長度提出動(dòng)態(tài)檢測數(shù)據(jù)驅(qū)動(dòng)的軌道不平順劣化自適應(yīng)分段建模方法MDL-RTDAS.在難以確保維修作業(yè)記錄完整，檢測數(shù)據(jù)全都有效的情況下，MDL-RTDAS 能夠克服檢測數(shù)據(jù)異常的影響，靈敏感知到軌道劣化過程中劣化趨勢的變化和自動(dòng)識(shí)別軌道不平順復(fù)雜劣化過程中的維修作業(yè)，實(shí)現(xiàn)軌道不平順劣化自適應(yīng)分段建模.將MDL-RTDAS應(yīng)用到2014—2019年昌福高速鐵路下行K21+184～K220+308的軌道劣化分析中，對(duì)其有效性進(jìn)行驗(yàn)證.與同類算法對(duì)比，MDL-RTDAS更加有效.為進(jìn)一步提高M(jìn)DL-RTDAS 的有效性和精度，需要優(yōu)化求解算法.考慮應(yīng)用MDLRTDAS建立的軌道不平順劣化分段擬合模型來探索軌道劣化規(guī)律，輔助工務(wù)部門制定養(yǎng)護(hù)維修計(jì)劃是下一步的研究重點(diǎn).