楊友濤,　劉國祥

(1. 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院, 四川 成都 611756; 2. 西南交通大學(xué)高速鐵路運(yùn)營安全空間信息技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室, 四川 成都 611756; 3. 西南交通大學(xué)土木工程學(xué)院, 四川 成都 610031)

軌道不平順是指軌道的實(shí)際位置、幾何尺寸相對設(shè)計(jì)位置和設(shè)計(jì)尺寸的偏差,是由很多隨機(jī)因素共同作用的結(jié)果.主要包括線路施工過程中的偏差、鋼軌的初始平直性、道床、路基等構(gòu)筑物的不均勻沉降及殘余變形積累,以及機(jī)車車輛等荷載的不規(guī)則變化的動(dòng)力作用[1]等因素.由于這些不確定性因素的共同作用,使得軌道不平順具有隨機(jī)性,主要采用均值、方差等統(tǒng)計(jì)數(shù)據(jù)表示其特性.目前,我國軌道幾何狀態(tài)評價(jià)主要采用局部幅值超限方法和軌道質(zhì)量指數(shù)方法,這兩種方法從空間域(線路里程)幅值方面進(jìn)行評價(jià),具有一定的局限性.

有學(xué)者利用信號的時(shí)頻分析方法對軌道動(dòng)態(tài)不平順數(shù)據(jù)進(jìn)行分析,文獻(xiàn)[2-3]中利用小波分析理論對軌道不平順局部高頻災(zāi)害進(jìn)行檢測，提取不同波長范圍的軌道不平順.文獻(xiàn)[4-5]中利用小波和Wigner-Ville相結(jié)合的方法對軌道不平順特征進(jìn)行提取與分析.文獻(xiàn)[6]中利用短時(shí)傅里葉變換、魏格納變換和Gabor變換等時(shí)頻分析方法分析了軌道實(shí)測振動(dòng)加速度信號,提取了長波不平順信號.短時(shí)傅立葉變換采用固定時(shí)間窗口的函數(shù)進(jìn)行計(jì)算,不能根據(jù)信號的特點(diǎn)來調(diào)整時(shí)間分辨率和頻率分辨率,沒有自適應(yīng)能力.小波變換克服了短時(shí)傅立葉變換時(shí)間窗口大小不隨頻率變化的缺點(diǎn),可根據(jù)信號頻率的高低自適應(yīng)調(diào)整時(shí)頻窗口的大小,但仍受Heisenberg不確定原理限制,不能精確表示頻率隨時(shí)間的變化.另外,各個(gè)小波基函數(shù)的適用范圍不一樣,選擇適合待分析信號的小波基函數(shù)是小波分析的難點(diǎn),使小波分析不具有自適應(yīng)性.

與以上建立在先驗(yàn)基函數(shù)上的時(shí)頻分析方法不同,經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, EMD)是以信號本身的時(shí)間尺度特征為基礎(chǔ)的自適應(yīng)分解方法,更適合非線性、非平穩(wěn)信號分析.經(jīng)驗(yàn)?zāi)B(tài)分解與希爾伯特變換相結(jié)合,更能揭示隱藏于信號中的內(nèi)在物理信息[7].針對經(jīng)驗(yàn)?zāi)B(tài)分解僅能處理一元數(shù)據(jù)的不足, Rehman提出適合于多通道信號分析的多元經(jīng)驗(yàn)?zāi)B(tài)分解(multivariate empirical mode decomposition, MEMD)[8].多元經(jīng)驗(yàn)?zāi)B(tài)分解把多元數(shù)據(jù)分解為具有相同數(shù)目的本征模態(tài)函數(shù),可以減少各層模態(tài)混疊和尺度不對齊的問題,具有自適應(yīng)性[9].本文利用多元經(jīng)驗(yàn)?zāi)B(tài)分解(intrinsic mode functions, IMF)和希爾伯特變換方法,對軌道幾何不平順的時(shí)頻特征分布進(jìn)行分析,為線路養(yǎng)護(hù)維修、列車安全運(yùn)行提供技術(shù)支持.

1　希爾伯特-黃變換方法

1.1　多元經(jīng)驗(yàn)?zāi)B(tài)分解

經(jīng)驗(yàn)?zāi)B(tài)分解是以數(shù)據(jù)驅(qū)動(dòng)方式把非線性、非平穩(wěn)信號分解為多個(gè)幅值-頻率調(diào)制的本征模態(tài)函數(shù)以及殘余分量[7],其表達(dá)式為

(1)

式中:x(t)為信號;cm(t)為IMF;m為IMF層數(shù);r(t)為剩余分量;t為時(shí)間.

本征模態(tài)函數(shù)包含不同的頻率成份,反映了信號的不同波動(dòng)情況.每個(gè)本征模態(tài)函數(shù)必須滿足條件[7]:(1) 在整個(gè)數(shù)據(jù)范圍內(nèi),極值點(diǎn)和過0點(diǎn)的個(gè)數(shù)相等或相差1個(gè);(2) 由極大值和極小值計(jì)算的包絡(luò)值之和為0.

經(jīng)驗(yàn)?zāi)B(tài)分解通過循環(huán)篩分的方法把滿足上述兩個(gè)條件的IMF分解出來.在一元經(jīng)驗(yàn)?zāi)B(tài)分解中,通過對信號的極大值和極小值分別內(nèi)插計(jì)算獲取上、下包絡(luò)線,然后對其取平均得到局部均值[7],對于多元數(shù)據(jù),不能直接定義極值.

為解決這個(gè)問題,多元經(jīng)驗(yàn)?zāi)B(tài)分解把多元信號X(t)=(x1(t),x2(t),…,xn(t))沿不同方向投影到n維空間上,計(jì)算這些投影序列的極值并進(jìn)行內(nèi)插,得到投影序列的包絡(luò)線,這些包絡(luò)線的平均值即為多元信號的局部均值[8].多元經(jīng)驗(yàn)?zāi)B(tài)分解的計(jì)算過程如下:

(1) 在n-1維超球面上選擇合適的均勻采樣點(diǎn)集θk,獲得n維空間的方向向量集vθk,k為方向向量個(gè)數(shù),k=1,…K;

(2) 計(jì)算多元信號X(t)在每個(gè)方向向量vθk上的投影向量Pθk(t);

(3) 找出所有投影向量集中,每個(gè)投影向量極大值所對應(yīng)的時(shí)刻ti,θk,i為極大值點(diǎn)的序號;

(4) 利用樣條函數(shù)對極值點(diǎn)(ti,θk,X(ti,θk))進(jìn)行插值,獲得K條包絡(luò)曲線eθk(t);

(5) 計(jì)算K個(gè)方向向量的包絡(luò)線均值

(2)

(6) 提取高頻細(xì)節(jié)信號d(t),d(t)=X(t)-m(t).如果d(t)滿足多元IMF的停止準(zhǔn)則,記C1(t)=d(t)為第 1層多元IMF,將X(t)-d(t)作為(2)中輸入信號,繼續(xù)提取多元IMF;否則,將d(t)作為輸入信號繼續(xù)迭代計(jì)算.

經(jīng)過多元經(jīng)驗(yàn)?zāi)B(tài)分解后,多元信號分解為不同尺度的多元IMF和多元剩余分量之和的形式為

(3)

式中Cm(t)=(c1，m(t),c2，m(t),…,cn，m(t))為n元IMF;R(t)=(r1(t),r2(t),…,rn(t))為n元剩余分量.

1.2　希爾伯特譜分析

通常,瞬時(shí)頻率和瞬時(shí)幅值用于表征非平穩(wěn)、非線性信號的短暫特性.但是,對于多分量組成的非平穩(wěn)信號,瞬時(shí)頻率沒有明顯的物理意義.所以,首先利用經(jīng)驗(yàn)?zāi)B(tài)分解將多分量信號分解為一系列單分量的IMF,然后利用希爾伯特變換計(jì)算其瞬時(shí)頻率和瞬時(shí)幅值.對于單一成分的本征模態(tài)函數(shù)cj,m(t),其希爾伯特變換為

(4)

式中:P為柯西主值；j為多元變量序號.

相應(yīng)的解析信號形式為

Zj，m(t)=cj,m(t)+iH(cj,m(t))=aj,m(t)eiφj,m(t).

(5)

信號瞬時(shí)幅值aj,m(t)和瞬時(shí)相位φj,m(t)的計(jì)算式為

(6)

通過對瞬時(shí)相位求導(dǎo),可得瞬時(shí)頻率為

(7)

通過希爾伯特變換計(jì)算瞬時(shí)頻率,有時(shí)會(huì)出現(xiàn)較大的波動(dòng),甚至出現(xiàn)負(fù)頻率的情況.同時(shí),希爾伯特變換具有端點(diǎn)效應(yīng),即在信號兩端的瞬時(shí)頻率具有能量泄漏的問題[10].文獻(xiàn)[10]中在希爾伯特變換求取瞬時(shí)頻率的基礎(chǔ)上,提出標(biāo)準(zhǔn)化希爾伯特變換,采用經(jīng)驗(yàn)調(diào)頻-調(diào)幅分解方法,減弱了IMF中的調(diào)幅分量,僅保留了調(diào)頻分量,從而使瞬時(shí)頻率的計(jì)算不再受Bedrosian等式的限制.文中軌道不平順的瞬時(shí)頻率采用標(biāo)準(zhǔn)化希爾伯特變換計(jì)算.

通過對每一個(gè)IMF進(jìn)行希爾伯特變換,信號xj(t)表示為

(8).

用式(8)可在三維坐標(biāo)圖中,將瞬時(shí)幅值和瞬時(shí)頻率表示為時(shí)間的函數(shù),即在時(shí)間和頻率聯(lián)合平面內(nèi)表示幅值的變化.這種幅值隨時(shí)間和頻率(f)變化的分布函數(shù)為信號的希爾伯特譜,其表達(dá)式為

(9)

2　軌道幾何不平順數(shù)據(jù)分析

軌道幾何不平順用軌距、水平、左高低、右高低、左軌向、右軌向和扭曲7個(gè)參數(shù)來描述,其表示為隨里程變化的隨機(jī)函數(shù).

本文選取軌道檢查車在滬寧城際上的軌道動(dòng)態(tài)檢查數(shù)據(jù)進(jìn)行分析,檢測速度為252 km/h,樣本數(shù)據(jù)里程為上行K104+000.00～K104+999.75,軌檢車數(shù)據(jù)采樣間隔為0.25 m.

2.1　軌道不平順多尺度分解

對預(yù)處理后的軌道不平順信號進(jìn)行多元經(jīng)驗(yàn)?zāi)B(tài)分解,得到13層IMF(C1～C13)和1個(gè)趨勢項(xiàng),分解結(jié)果見圖1、2.

由圖1、2可以看出,隨著IMF序號的增加,其空間變化尺度逐漸增大;7個(gè)軌道不平順參數(shù)中,相同序號的IMF具有相近的空間變化尺度.根據(jù)多元經(jīng)驗(yàn)?zāi)B(tài)分解這一特點(diǎn),按信號內(nèi)在規(guī)律分解出多元序列相近尺度的分量.軌距和右軌向不平順的第2層和第3層IMF在里程K104+485.00前后有大的幅值異常值,見圖示.在數(shù)據(jù)預(yù)處理階段采用3倍中誤差準(zhǔn)則并沒有探測出該里程位置的異常值,但是經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解后有異常值出現(xiàn).

2.2　軌道不平順頻率分析

對軌道不平順分解的各層IMF進(jìn)行標(biāo)準(zhǔn)化希爾伯特變換,計(jì)算各個(gè)IMF的瞬時(shí)頻率和瞬時(shí)幅值.

為直觀表示軌道不平順各IMF瞬時(shí)頻率的變化,計(jì)算其平均頻率.計(jì)算平均頻率有兩種方法:

(1) 根據(jù)極值點(diǎn)或過0點(diǎn)個(gè)數(shù)計(jì)算平均頻率[11-12];

(2) 利用信號的瞬時(shí)頻率計(jì)算平均頻率[13-14].

按極值點(diǎn)和過0點(diǎn)個(gè)數(shù)確定平均頻率的方法,在高階IMF中,由于極值點(diǎn)和過0點(diǎn)位置的不均勻性,平均頻率的計(jì)算結(jié)果誤差較大.

文獻(xiàn)[14]中利用瞬時(shí)頻率的算術(shù)平均值和能量加權(quán)均值表示平均頻率,能量權(quán)重平均頻率可寫為

(10)

式中:L為時(shí)間長度.

根據(jù)瞬時(shí)頻率所占的權(quán)比來計(jì)算平均頻率的幾種方法,更能直觀的反映每層IMF的頻率變化情況[13-14].本文采用瞬時(shí)頻率能量權(quán)重的方法,即按式(10)計(jì)算軌道不平順各IMF的平均頻率及平均頻率標(biāo)準(zhǔn)差.與式(10)相對應(yīng),平均頻率標(biāo)準(zhǔn)差也考慮能量權(quán)重,則加權(quán)標(biāo)準(zhǔn)差計(jì)算式為

(11)

圖1　軌道幾何不平順多元經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果 (A)Fig.1　MEMD results of irregularities in track geometry (A)

圖2　軌道幾何不平順多元經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果(B)Fig.2　MEMD results of irregularities in track geometry (B)

為減少希爾伯特變換邊界效應(yīng)對平均頻率的影響,在軌道不平順的每個(gè)IMF平均頻率及標(biāo)準(zhǔn)差計(jì)算中剔除了樣本數(shù)據(jù)兩端的部分?jǐn)?shù)據(jù),其計(jì)算結(jié)果如表1所示.表中,相對標(biāo)準(zhǔn)差指每個(gè)IMF的平均頻率的標(biāo)準(zhǔn)差與其平均頻率的比值.

從表1可以看出:(1) 軌道不平順7個(gè)參數(shù)中,經(jīng)多元經(jīng)驗(yàn)?zāi)B(tài)分解得出的同層IMF的平均頻率基本一致,每層IMF中平均頻率相對較差最大的為第13層的IMF,其原因是高階IMF的平均頻率計(jì)算誤差較大造成的;其他各層的平均頻率相對較差很小,說明經(jīng)多元經(jīng)驗(yàn)?zāi)B(tài)分解的軌道不平順各層的IMF頻率具有一致性[9].(2) 每個(gè)軌道不平順參數(shù)的各層IMF的平均頻率逐漸減少,并且其值約為上一層IMF平均頻率的1/2.文獻(xiàn)[9]中通過對白色噪聲進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解,得出該方法具有二進(jìn)濾波特性,軌道不平順數(shù)據(jù)的頻率分布與白色噪聲具有相似性,說明軌道不平順數(shù)據(jù)具有隨機(jī)特性.(3) 不同序列的同層IMF的平均頻率相對標(biāo)準(zhǔn)差差別較小,其中：第1層IMF的相對標(biāo)準(zhǔn)差較大,其值為0.4左右,說明第1層瞬時(shí)頻率波動(dòng)較大,可能是信號中含有高頻噪聲所致;第2～13層IMF的相對標(biāo)準(zhǔn)差均較小,說明多元經(jīng)驗(yàn)?zāi)B(tài)分解方法獲得的軌道不平順的頻率帶寬較窄,瞬時(shí)頻率波動(dòng)較小.

表1　軌道不平順各層IMF的平均頻率和相對標(biāo)準(zhǔn)差Tab.1　Average frequency and relative standard deviation of IMFs of track irregularity

多元經(jīng)驗(yàn)?zāi)B(tài)分解方法是根據(jù)數(shù)據(jù)本身的波動(dòng)情況自適應(yīng)的分解為多層IMF.本文選取的軌道不平順樣本數(shù)據(jù)被分解為13層,其平均波長(空間頻率的倒數(shù))范圍為0.7～761.4 m,說明軌道不平順數(shù)據(jù)是由不同波長組成的多分量隨機(jī)函數(shù).

目前軌道檢查車的最大檢測波長為200 m,由多元經(jīng)驗(yàn)?zāi)B(tài)分解方法得出的最大波長遠(yuǎn)超此范圍,說明利用多元經(jīng)驗(yàn)?zāi)B(tài)分解方法獲得的IMF有虛假波長成份,所以如何根據(jù)軌道不平順數(shù)據(jù)的實(shí)際波動(dòng)情況選擇有效的IMF也是需要解決的問題.

2.3　軌道不平順能量分布

文獻(xiàn)[14]中在分析小波尺度圖的基礎(chǔ)上,提出多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖,即用各層IMF的瞬時(shí)幅值的平方來表示信號不同尺度的能量變化.根據(jù)式(8),信號的多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖可表示為

(12)

多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖在時(shí)間和尺度的聯(lián)合平面內(nèi)表示能量的變化情況,軌道不平順7個(gè)參數(shù)的多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖如圖3所示.

由圖3可以看出,軌距和水平不平順的能量分布較為廣泛,涵蓋在第3～11層IMF范圍內(nèi),軌距的第11層IMF能量所占比重較大,其能量主要分布在低頻部分;軌向和高低4個(gè)不平順參數(shù)的能量主要集中在中間頻段的IMF中,低頻和高頻段范圍的能量幾乎可以忽略不計(jì);扭曲不平順的能量主要分布在前6個(gè)IMF中,即高頻部分.

綜上所述，通過多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖也能確定軌道不平順在不同尺度內(nèi)沿里程方向上能量的突變,可為精確查找軌道病害,進(jìn)行養(yǎng)護(hù)維修提供了基礎(chǔ).

表2為根據(jù)多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖統(tǒng)計(jì)出的各層IMF能量所占的比例.由表2可以看出：

(1) 軌距不平順中，第11層IMF所占能量達(dá)到40.9%,主要原因是在里程K104+500.00～K104+980.00范圍內(nèi)軌距不平順幅值有大的波動(dòng).

(2) 軌向和高低能量分布比較集中,第4～8層 的IMF能量之和比例達(dá)到90%以上,從表1可知，第4～8層的IMF對應(yīng)的空間平均波長為4.37、7.58、14.11、23.96、35.56 m,即能量主要分布在波長為4～36 m的范圍;同時(shí)在每層IMF中左右軌向、左右高低的能量百分比比較一致,說明左右軌道的軌向、高低的能量分布具有相關(guān)性.

(3) 扭曲不平順能量主要集中在高頻部分,其中第4、5層IMF能量所占比例均大于30%；第8～13層IMF所占能量比例之和僅為1.5%,低頻段能量的影響很小.

多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖能清晰的標(biāo)識軌道不平順在不同尺度上的能量分布.軌道不平順樣本數(shù)據(jù)中,除軌距不平順外,其他不平順的能量主要分布在中波區(qū)段,且能量分布相對集中.通過以上分析,能更好的掌握軌道不平順的能量和波長分布規(guī)律,可以有針對性的對軌道進(jìn)行養(yǎng)護(hù)維修,提高鐵路運(yùn)行的安全性.

圖3　軌道不平順多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖Fig.3　MEMD-based scalograms of track irregularity

IMF層號軌距水平左軌向右軌向左高低右高低扭曲C10.1 3.60.80.70.40.24.8C20.35.21.81.61.30.75.5C31.16.74.13.72.72.37.2C45.519.015.215.212.815.030.1C55.417.418.719.017.417.934.0C66.513.520.122.127.823.512.7C78.110.724.427.925.427.44.2C85.26.014.39.210.911.70.9C96.18.60.50.50.80.90.5C1019.15.10.20.10.30.30.1C1140.92.90.00.00.10.10.0C120.70.80.00.00.00.00.0C130.90.50.00.00.00.00.0

3　結(jié)　論

基于多元經(jīng)驗(yàn)?zāi)B(tài)分解方法和希爾伯特變換,對高速鐵路軌道動(dòng)態(tài)不平順進(jìn)行時(shí)頻特征分析和應(yīng)用研究,得出以下主要結(jié)論:

(1) 多元經(jīng)驗(yàn)?zāi)B(tài)分解方法用于分析多元數(shù)據(jù)序列,減少了標(biāo)準(zhǔn)經(jīng)驗(yàn)?zāi)B(tài)分解方法中的模態(tài)混疊和多元序列尺度不對齊的問題,具有很好的自適應(yīng)性和時(shí)頻局部化能力.希爾伯特變換計(jì)算多元序列的瞬時(shí)頻率和瞬時(shí)幅值,能反映多元數(shù)據(jù)序列的局部時(shí)頻分布特征.

(2) 用多元經(jīng)驗(yàn)?zāi)B(tài)分解方法,分解出的軌道幾何不平順各層IMF的頻率帶寬較窄,軌道不平順在不同尺度下的頻率波動(dòng)較小;軌道不平順不同參數(shù)的同層IMF的平均頻率相差較少,具有一致性.

(3) 多元經(jīng)驗(yàn)?zāi)B(tài)分解尺度圖能直觀的反映軌道不平順的能量分布,樣本數(shù)據(jù)的軌距不平順能量主要分布于中長波區(qū)段,波長頻段范圍較廣;軌向和高低能量分布相對集中,主要集中于空間波長4～36 m范圍;扭曲的能量主要集中分布在波長為4.9、7.6 m的兩個(gè)尺度內(nèi).

為更好的提取軌道不平順的時(shí)頻特征,應(yīng)對整條線路軌道不平順數(shù)據(jù)按線路類型、軌道結(jié)構(gòu)以及線下結(jié)構(gòu)等分門別類的進(jìn)行特征分析,同時(shí)比較軌道不平順數(shù)據(jù)在不同檢測時(shí)間的時(shí)頻變化特征.