農(nóng)漢彪，曾巧妮

(百色學(xué)院，廣西 百色 533000)

0 引言

輪軌相互作用力監(jiān)測(cè)是軌道交通車輛運(yùn)行安全狀態(tài)中非常重要的方法和途徑。通過輪軌力的監(jiān)測(cè)可以實(shí)時(shí)了解車輛運(yùn)行過程中的減載率，脫軌系數(shù)等。在基于地面的輪軌力檢測(cè)中，通常通過測(cè)量鋼軌的動(dòng)態(tài)響應(yīng)來測(cè)量車輪通過時(shí)的輪軌相互作用力[1]。而在現(xiàn)場(chǎng)實(shí)測(cè)時(shí)，由于鋼軌支撐條件會(huì)存在一定的不確定性，比如懸枕，扣件松弛，以及地基彈性突變等情況，使得輪軌力測(cè)試信號(hào)中會(huì)出現(xiàn)基線漂移的干擾[2]，基線漂移的出現(xiàn)影響了后續(xù)的輪軌力信號(hào)的定量分析以及對(duì)列車運(yùn)行狀態(tài)的安全評(píng)估，所以必須予以去除。

去除基線漂移方法在心電圖信號(hào)處理[3-4]，脈搏信號(hào)處理[5]，眼電圖處理[6]等生物醫(yī)學(xué)領(lǐng)域中廣泛應(yīng)用，在光學(xué)測(cè)試[7]和地震監(jiān)測(cè)[8]中也有了相關(guān)的文獻(xiàn)論述。而針對(duì)輪軌力信號(hào)降噪和基線漂移的消除，目前國(guó)內(nèi)外較少涉及。

文獻(xiàn)[9]中提出了基于分段數(shù)據(jù)中值聚類的算法，該算法將信號(hào)數(shù)據(jù)分成很多分段，取分段數(shù)據(jù)的中值，將中值的時(shí)間間距與幅值間距作為距離測(cè)度進(jìn)行最近距離聚類，得到了能代表基線漂移的基點(diǎn)，擬合基點(diǎn)得到基線漂移干擾。算法效果良好，然而算法中需要設(shè)置不同測(cè)度之間的協(xié)調(diào)因子，且該因子的設(shè)置與數(shù)據(jù)長(zhǎng)度，采樣頻率，數(shù)據(jù)幅值等因數(shù)相關(guān)。文獻(xiàn)[10]利用LabVIEW提供的高級(jí)信號(hào)處理工具包和數(shù)字濾波器設(shè)計(jì)工具包對(duì)實(shí)測(cè)輪軌力應(yīng)變信號(hào)進(jìn)行處理,內(nèi)容包括消除基準(zhǔn)漂移、去除寬帶噪聲和提取信號(hào)特征等。文獻(xiàn)[11]提出一種小波包降噪算法與經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)相關(guān)算法相結(jié)合的數(shù)據(jù)處理方法，有效地消除輪軌力在采集過程中隨機(jī)白噪聲的干擾和抑制模態(tài)的混疊。文獻(xiàn)[12]提出了基于集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)與小波變換相結(jié)合的去噪方法。對(duì)含有基線漂移的分量通過小波變換進(jìn)行分解，將代表基線漂移的趨勢(shì)項(xiàng)置零達(dá)到去除基線漂移的目的。

近年來，基線漂移的去除方法趨向于將原信號(hào)通過不同的算法分解成純凈信號(hào)部分和干擾成分，主要算法有基EMD及其改進(jìn)算法，如EEMD[13], CEEMDAN[14]等，有經(jīng)驗(yàn)小波分解(EWT)[15]，傅里葉分解[16]，變分模態(tài)分解(VMD)[17]等，而信號(hào)的分解總會(huì)存在頻率混疊問題，傳統(tǒng)的方法有FIR或IIR的方法，中值濾波[18]，曲線擬合[19]和數(shù)學(xué)形態(tài)學(xué)濾波[20]等。曲線擬合方法是將信號(hào)中能代表基線漂移走勢(shì)的基點(diǎn)擬合成曲線并從原始信號(hào)中除去。該方法高效、準(zhǔn)確，計(jì)算量小，但是從原始信號(hào)中提取和確定基點(diǎn)比較困難。本文根據(jù)輪軌力信號(hào)短時(shí)分段數(shù)據(jù)的高階統(tǒng)計(jì)量分布和基線漂移的關(guān)系，提出了基于分段數(shù)據(jù)高階統(tǒng)計(jì)量聚類分析的方法獲取基點(diǎn)所在的數(shù)據(jù)分段，以數(shù)據(jù)分段中值作為基點(diǎn)，再擬合即可得到基線漂移。

1 算法思想與實(shí)現(xiàn)

鋼軌的模型通常采用連續(xù)彈性基礎(chǔ)梁進(jìn)行表示，鋼軌垂直位移與移動(dòng)載荷的大小，速度，位置之間關(guān)系表示為[21]：

(1)

其中：EI為鋼軌的抗彎剛度，k為地基的彈性系數(shù)，m為單位長(zhǎng)度鋼軌質(zhì)量，P為鋼軌上載荷的大小，x為載荷作用點(diǎn)與觀測(cè)點(diǎn)之間的距離，v為載荷移動(dòng)的速度，δ(·)為狄拉克函數(shù)。

靜態(tài)或準(zhǔn)靜態(tài)時(shí)，在多個(gè)載荷的同時(shí)作用下，鋼軌上某觀測(cè)點(diǎn)的剪力為：

(2)

在實(shí)測(cè)信號(hào)中，測(cè)試系統(tǒng)采集到信號(hào)除了鋼軌載荷引起的響應(yīng)之外，還包含了基線漂移成分和其他的隨機(jī)干擾，即采集到的信號(hào)f(t)可以表示為：

f(t)=Q(t)+BW(t)+n(t)

(3)

其中：Q(t)為車輪載荷所引起的剪力測(cè)量值；BW(t)為基線漂移信號(hào)，n(t)為白噪聲干擾信號(hào)，當(dāng)所有|xi|>3πL/2時(shí)，Q(t)降為峰值的1%以內(nèi)，此時(shí)BW(t)+n(t)Q(t)，稱無載荷區(qū)間即：

f(t)≈BW(t)+n(t)

(4)

根據(jù)隨機(jī)噪聲的對(duì)稱性和基線漂移信號(hào)的緩變性，當(dāng)所有|xi|>3πL/2時(shí)，短時(shí)間段內(nèi)測(cè)試信號(hào)的中值即可表示為基線漂移成分的幅值，稱為基點(diǎn)。而當(dāng)存在|xi|<3πL/2時(shí)，測(cè)量值f(t)包含鋼軌對(duì)載荷的響應(yīng)，其數(shù)據(jù)段的中值不能表示為基點(diǎn)。擬合所有基點(diǎn)即可達(dá)到信號(hào)的基線漂移成分。而前提是正確篩選出能被視為基點(diǎn)的數(shù)據(jù)段中值數(shù)據(jù)。

對(duì)于f(t)，由于Q(t)，BW(t)，n(t)相互獨(dú)立，則有f(t)的方差：

D(f(t))=D(Q(t)+BW(t)+n(t))=

D(Q(t))+D(BW(t))+D(n(t))

(5)

BW(t)為緩變信號(hào)，在短時(shí)間內(nèi)BW(t)為可視為常數(shù)，即D(BW(t))=0，而白噪聲的方差D(n(t))可視為常數(shù)C;從而有當(dāng)所有|xi|>3πL/2時(shí)：

D(f(t))=0+0+C=C

(6)

而在|xi|<3πL/2時(shí)：

D(f(t))=D(Q(t))+0+C=D(Q(t))+C

(7)

另外，同樣可以看到，在短時(shí)間內(nèi)信號(hào)的峭度值：

(8)

其中μ為f(t)的數(shù)學(xué)期望，σ為f(t)的標(biāo)準(zhǔn)差，由于Q(t)，BW(t)，n(t)相互獨(dú)立，有：

K(f(t))=K(Q(t))+K(BW(t))+K(n(t))

(9)

短時(shí)間內(nèi)BW(t)是均勻分布，即K(BW(t))=0，n(t)為正態(tài)分布，K(n(t))=3；從而有當(dāng)所有|xi|>3πL/2時(shí)：

K(f(t))=0+0+3=3

(10)

而在|xi|<3πL/2時(shí)：

K(f(t))=K(Q(t))+0+3=K(Q(t))+3

(11)

聯(lián)合方差與峭度，無載荷期間內(nèi)數(shù)據(jù)段的D(f(t))，和K(f(t))將主要分布在以方差和峭度為坐標(biāo)軸的點(diǎn)(0，3)和點(diǎn)(C,3)附近，而承載期間的坐標(biāo)位置將遠(yuǎn)離這兩點(diǎn)坐標(biāo)。根據(jù)以上分析，可以通過對(duì)分段數(shù)據(jù)的方差和峭度進(jìn)行基于密度的聚類分析來篩選出無載荷期間的數(shù)據(jù)段。

基于密度的聚類分析可以將具有相同或相近特征的對(duì)象樣本之間的分為同一類，即在某類別任意樣本周圍一定空間范圍內(nèi)一定有同類別的樣本存在，能有效區(qū)分有載荷期間數(shù)據(jù)段和無載荷期間數(shù)據(jù)段?；诿芏鹊木垲惙治龅某Ｓ盟惴ㄓ蠨BSCAN，OPTICS和DENCLUE等。DBSCAN需要設(shè)定點(diǎn)數(shù)閾值minPts和鄰域半徑，全局密度，DENCLUE需要設(shè)定密度閥值和參數(shù)，這兩種算法的效果對(duì)參數(shù)非常敏感，需要確定合適的參數(shù)才能得到理想的效果。OPTICS是DBSCAN的一種改進(jìn)方案，也要設(shè)定點(diǎn)數(shù)閾值minPts和鄰域半徑，但聚類結(jié)果對(duì)參數(shù)不敏感。

基于以上分析，應(yīng)用基于密度的聚類分析可以有效找到信號(hào)中基點(diǎn)，從而實(shí)現(xiàn)基線漂移的消除。該方法的具體過程和步驟如下：

1)對(duì)測(cè)量數(shù)據(jù)進(jìn)行分段。前后數(shù)據(jù)段之間可以部分重疊，以加大數(shù)據(jù)段數(shù)量和保證數(shù)據(jù)的連續(xù)性；

2)計(jì)算每數(shù)據(jù)段方差和峭度值，并分別進(jìn)行歸一化處理；

3)針對(duì)歸一化方差和峭度兩個(gè)維度進(jìn)行基于密度的聚類分析；

4)篩選出聚類結(jié)果中處于零點(diǎn)附近的一個(gè)分類或幾個(gè)分類；

5)對(duì)選用分類所對(duì)應(yīng)的數(shù)據(jù)段取中值，即為基點(diǎn)；

6)曲線擬合所有基點(diǎn)得到基線漂移干擾。

分段時(shí)，數(shù)據(jù)段的大小需要根據(jù)信號(hào)自身的成分來確定。涉及到信號(hào)數(shù)據(jù)的采樣頻率，通過列車的轉(zhuǎn)向架間距，軸距和運(yùn)行速度等。分塊大小的最優(yōu)設(shè)置是能將信號(hào)的負(fù)荷狀態(tài)和非負(fù)荷狀態(tài)完全分開，這是理想的情況，而實(shí)際上由于信號(hào)數(shù)據(jù)初始采集相位未知且列車運(yùn)行速度并不恒定原因，無法在未識(shí)別壓頭時(shí)刻的情況下實(shí)現(xiàn)理想分段。而總會(huì)出現(xiàn)部分分段橫跨負(fù)荷狀態(tài)和非負(fù)荷狀態(tài)兩種狀態(tài)的情況。根據(jù)數(shù)據(jù)試驗(yàn)結(jié)果一般采用的分大小可為：

(12)

分段數(shù)據(jù)的方差和峭度都是非負(fù)數(shù)值，但具有不同的數(shù)值范圍，基于密度的聚類需要考慮數(shù)據(jù)點(diǎn)之間的距離測(cè)度或相似度測(cè)度，為了均衡方差和峭度在距離測(cè)度上的貢獻(xiàn)，需要對(duì)方差和峭度進(jìn)行歸一化處理。歸一化時(shí)可以將最大方差和最大峭度置為1，最小值置為零，其他進(jìn)行線性化處理。

在聚類分析時(shí)，選用對(duì)參數(shù)不敏感的OPTICS算法，OPTICS所需要點(diǎn)數(shù)閾值minPts和鄰域半徑ε兩個(gè)參數(shù)可以通過自適應(yīng)的方式來確定。鄰域半徑ε可以取所有數(shù)據(jù)樣點(diǎn)兩兩間距測(cè)度密度分布函數(shù)的最快下降梯度點(diǎn)所對(duì)應(yīng)的距離測(cè)度，minPts可以取使得當(dāng)前ε下所有分類間距總和最小的閾值。

2 仿真數(shù)據(jù)分析

為了定量分析基線消除方法的效果，預(yù)先準(zhǔn)備一組零漂移輪軌垂直力信號(hào)和一組已知基線漂移信號(hào)。假設(shè)列車為兩軸轉(zhuǎn)向架，車輛定距為7.67 m，轉(zhuǎn)向架軸距為1.75 m，車輛長(zhǎng)度為11.986 m，列車低速經(jīng)過觀測(cè)點(diǎn)的車速為4 km/h，信號(hào)采樣頻率為1 000 Hz，數(shù)據(jù)長(zhǎng)度為60 000點(diǎn)，各輪軌力信號(hào)峰值隨機(jī)分布在1.4～1.7之間零漂移輪軌垂直力信號(hào)圖1所示，數(shù)據(jù)本身具有很低的頻率成分，該頻率成分與列車轉(zhuǎn)向架間距，軸距與列車速度都有關(guān)系?；€漂移信號(hào)為：

圖1 零漂移鋼軌垂直力信號(hào)

BW(t)=0.15*(sin(0.15πt)+sin(0.09πt+0.2π))

(14)

采用與力信號(hào)相同的采樣間隔?；€及與力信號(hào)疊加后的結(jié)果如圖2所示。

圖2 基線及與力信號(hào)疊加的結(jié)果

常用于定量分析濾波算法性能評(píng)估參數(shù)主要有均方誤差(MSE)和信噪比(SNR)等[22]。對(duì)于性能評(píng)價(jià)，MSE值越小越好，SNR越大越好。它們的定義分別為：

(15)

(16)

根據(jù)分段大小的分析，仿真數(shù)據(jù)可以采用分段大小為500數(shù)據(jù)點(diǎn)對(duì)含有基線漂移的信號(hào)數(shù)據(jù)進(jìn)行分段，分段窗口步進(jìn)長(zhǎng)度為250。分段后的數(shù)據(jù)分段中值如圖3所示，大部分分段數(shù)據(jù)的中值和已知的基線漂移重合，另一部分偏離較遠(yuǎn)。

圖3 含基線信號(hào)分段取中值的結(jié)果

隨后對(duì)所有數(shù)據(jù)分段統(tǒng)計(jì)數(shù)據(jù)分段的方差和峭度，得到的方差和峭度數(shù)據(jù)如圖4所示。

方差數(shù)值在0～1范圍之內(nèi)，而峭度數(shù)值在0～80范圍之內(nèi)，分別對(duì)方差和峭度進(jìn)行線性歸一化，然后對(duì)歸一化后的數(shù)據(jù)進(jìn)行聚類分析得到如圖5的分類結(jié)果。

圖5 聚類分析結(jié)果

圖5中，數(shù)據(jù)點(diǎn)被分成了兩個(gè)分類和部分未分類樣點(diǎn)。根據(jù)輪軌力數(shù)據(jù)的特點(diǎn)，選用分類1作為無載荷期間的數(shù)據(jù)段。對(duì)分類1所對(duì)應(yīng)的數(shù)據(jù)分段中值即為信號(hào)數(shù)據(jù)的基點(diǎn)?；c(diǎn)數(shù)據(jù)與目標(biāo)基線漂移的位置關(guān)系如圖6所示。

圖6 選用基點(diǎn)與原始信號(hào)

為了對(duì)比分析，利用移動(dòng)中值濾波，形態(tài)學(xué)濾波，(BEADS，baseline estimation and denoising with sparsity)[23]，CEEMDAN，EWT，VMD，零相位濾波器以及本文提出的方法分別對(duì)圖2中含有基線漂移的力信號(hào)進(jìn)行處理。在這些方法可分為兩類，一類是CEEMDAN，EWT，VMD和零相位濾波器，另一類是移動(dòng)中值濾波，形態(tài)學(xué)濾波，BEADS和本文的方法。第一類主要是信號(hào)分解與重構(gòu)的方法和數(shù)字濾波的方法，這些方法會(huì)將所有的測(cè)量數(shù)據(jù)值進(jìn)行計(jì)算，這使得在有載荷期間的基線漂移的提取會(huì)受到載荷大小的影響。各種方法通過優(yōu)化設(shè)置參數(shù)得到最優(yōu)的結(jié)果如圖7所示。所得結(jié)果與目標(biāo)基線存在較大的誤差，不適合后續(xù)的定量分析。

圖7 第一類方法提取基線漂移結(jié)果

第二類方法中，設(shè)置移動(dòng)中值濾波的移動(dòng)窗口大小為700數(shù)據(jù)點(diǎn)，移動(dòng)步進(jìn)為一個(gè)數(shù)據(jù)點(diǎn)，然后再經(jīng)過1 Hz低通濾波器得到最優(yōu)的平滑基線；數(shù)學(xué)形態(tài)學(xué)濾波采用長(zhǎng)度為350數(shù)據(jù)點(diǎn)的橫條型結(jié)構(gòu)分別進(jìn)行開閉運(yùn)算然后取兩部分和的一半，再經(jīng)過1 Hz低通濾波器得到最優(yōu)的平滑基線；依據(jù)文獻(xiàn)[23]BEADS方法采用截止頻率為0.8 Hz，階數(shù)為1，不對(duì)稱參數(shù)為12，正則化參數(shù)為0.01。由于BEADS方法對(duì)任何信號(hào)提取的基線首尾均為零，需要先對(duì)原信號(hào)的首尾置零提取趨勢(shì)線再進(jìn)行提取，才能得到最佳的評(píng)估效果。不同方法提取的基線漂移與目標(biāo)基線的誤差如圖8所示

圖8 提取的基線漂移與目標(biāo)基線的誤差曲線

從基線漂移誤差曲線看到，移動(dòng)中值濾波和形態(tài)學(xué)濾波具有形同的變化規(guī)律，即會(huì)在觀測(cè)點(diǎn)負(fù)載時(shí)會(huì)有較大誤差而其他時(shí)刻的誤差趨于零，BEADS方法在信號(hào)的兩頭出現(xiàn)較大的誤差，而中間部分誤差較小。而本文的方法提取基線的誤差相對(duì)其他3種方法都要小。原始帶基線漂移的信號(hào)在經(jīng)過不同方法去除基線前后的相似測(cè)度參數(shù)如表1所示。

表1 不同方法去除基線漂移的效果

從表1的數(shù)據(jù)結(jié)果表明，本文的方法對(duì)于仿真數(shù)據(jù)的消除基線漂移的兩種評(píng)價(jià)指標(biāo)均優(yōu)于形態(tài)濾波等其他3種方法。本文方法的均方誤差MSE僅約為次好的BEADS算法的0.47%，而信噪比SNR則高出23 dB。

3 實(shí)測(cè)數(shù)據(jù)分析

實(shí)測(cè)數(shù)據(jù)中，除了基線漂移的干擾之外還會(huì)有其他噪聲的存在，如工頻干擾和隨機(jī)干擾等。如圖9所示為一組數(shù)據(jù)現(xiàn)場(chǎng)實(shí)測(cè)輪軌力信號(hào)。

圖9 實(shí)測(cè)輪軌力信號(hào)

依據(jù)算法，首先獲取數(shù)據(jù)峰值寬度以確定數(shù)據(jù)分段的大小，接著對(duì)實(shí)測(cè)信號(hào)進(jìn)行數(shù)據(jù)分段，統(tǒng)計(jì)分段數(shù)據(jù)的方差和峭度，基于方差和峭度統(tǒng)計(jì)數(shù)據(jù)段之間距離，按距離分布函數(shù)獲取最優(yōu)聚類領(lǐng)域半徑，然后進(jìn)行基于OPTICS算法進(jìn)行聚類分析，得到聚類結(jié)果如圖10所示，選用分類1所對(duì)應(yīng)的數(shù)據(jù)段，以該分類數(shù)據(jù)段的中值作為基點(diǎn)，擬合所有基點(diǎn)得到基線漂移成分，如圖11所示。

圖10 實(shí)測(cè)信號(hào)分段數(shù)據(jù)高階統(tǒng)計(jì)量聚類結(jié)果

圖11 實(shí)測(cè)輪軌信號(hào)基線漂移的基點(diǎn)的選取

實(shí)測(cè)信號(hào)沒有預(yù)知的基線漂移或純凈的信號(hào)數(shù)據(jù)，無法進(jìn)行量化的評(píng)估。為了驗(yàn)證本文方法的效果，采用移動(dòng)中值濾波，形態(tài)學(xué)濾波，BEADS方法分別對(duì)實(shí)測(cè)信號(hào)進(jìn)行處理，對(duì)比的方法通過調(diào)整相應(yīng)的參數(shù)使得處理結(jié)果達(dá)到最佳效果。圖12為不同方法消除實(shí)測(cè)輪軌垂直力信號(hào)基線漂移的效果。其中最上方的曲線為實(shí)測(cè)的原始信號(hào)，往下分別為移動(dòng)中值濾波加低通濾波的方法，形態(tài)學(xué)濾波加低通濾波的方法，BEADS方法以及本文方法去除基線漂移后的信號(hào)波形圖。

圖12 實(shí)測(cè)輪軌力信號(hào)基線漂移的消除效果比較

從圖上看，本文的方法和BEADS方法于實(shí)測(cè)數(shù)據(jù)中的基線漂移的消除效果明顯優(yōu)于前兩種方法。本文的方法和BEADS方法都能有效地去除基線漂移的干擾，僅從視覺上看，兩者無法區(qū)分兩者優(yōu)劣情況。從計(jì)算量上，由于本文的方法采用分段的方法，使得整個(gè)計(jì)算過程的計(jì)算量相比與移動(dòng)中值濾波、形態(tài)學(xué)濾波和BEADS方法大大的減小，提高了計(jì)算效率。

4 結(jié)束語

本文針對(duì)輪軌垂向力信號(hào)的基線漂移消除問題，提出了基于分段數(shù)據(jù)高階統(tǒng)計(jì)量聚類分析的方法來確定信號(hào)基點(diǎn)所在的數(shù)據(jù)段，統(tǒng)計(jì)選用數(shù)據(jù)段中值作為基點(diǎn)，曲線擬合基點(diǎn)即可得到信號(hào)的基線漂移干擾。文中通過仿真數(shù)據(jù)和現(xiàn)場(chǎng)實(shí)測(cè)信號(hào)從定量和定性兩種方式對(duì)所提出的方法進(jìn)行了驗(yàn)證，并與現(xiàn)有常用的基線漂移去除方法進(jìn)行了橫向比較。驗(yàn)證結(jié)果都表明相比于其他方法，本文的方法在定量上具有更高的信噪比和更小的均方誤差，在對(duì)實(shí)測(cè)信號(hào)的定性分析上同樣具有很好的效果，而且計(jì)算量更小，速度更快。本文方法的參數(shù)設(shè)置可以根據(jù)信號(hào)的特性自適應(yīng)設(shè)置，不需要過多的人為參與。該方法同樣適用于其他包含有基線偏移的測(cè)試信號(hào)的處理，比如心電信號(hào)，光譜信號(hào)，脈搏信號(hào)等。