基于Copula函數(shù)的高速列車轉(zhuǎn)向架故障特征提取

金煒東， 呂乾勇， 孫永奎

(西南交通大學(xué)電氣工程學(xué)院，四川 成都 610031)

由列車轉(zhuǎn)向架故障引起的振動嚴(yán)重影響列車的舒適性和安全性［1］.傳統(tǒng)列車故障診斷的方法有小波變換、聚合經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)、功率譜、能量譜等，小波變換［2］和 EEMD［3］方法是較為經(jīng)典的方法.小波變換方法是通過對故障信號進(jìn)行分解得到各個頻段的信號，然后提取各個頻段的信號，用熵值、能量等作為特征診斷故障.EEMD方法是通過對原始信號進(jìn)行分解得到各個本征模態(tài)函數(shù)(intrinsic mode function，IMF)，這些本征模態(tài)函數(shù)表征了原始信號在各個頻段下的分量.文獻(xiàn)［4］中對原始信號作EEMD變換得到各個階段的IMF，通過對每個IMF進(jìn)行單獨分析，提取各階段IMF的熵值作為特征.文獻(xiàn)［5］中對原信號進(jìn)行EEMD變換，選取與原始信號相關(guān)性較大的IMF，從選出的IMF中提取Teager能量譜作為特征用于軸承故障診斷.此類方法雖然能對分解得到的各個階段的IMF進(jìn)行分析，但由于只對各階段IMF進(jìn)行了單獨分析，不能捕捉到IMF之間的依賴關(guān)系.

Copula函數(shù)作為研究隨機(jī)變量之間關(guān)聯(lián)性的方法，近年來被廣泛應(yīng)用于金融市場分析［6］、水文分析［7］、紋理圖像識別［8-9］等領(lǐng)域.采用 Copula 函數(shù)的關(guān)聯(lián)性研究列車故障診斷方法能克服傳統(tǒng)方法對各個IMF進(jìn)行單獨研究的不足.

相對其他分析相關(guān)性的方法，Copula函數(shù)具有兩個明顯的優(yōu)點［7］.首先，靈活且求解簡單，構(gòu)造多維聯(lián)合分布函數(shù)時可以分為兩個步驟完成:構(gòu)建邊緣分布以及選擇合適的Copula函數(shù)構(gòu)建聯(lián)合分布;其次，不要求相同的邊緣分布，使得每種邊緣分布函數(shù)保留了自己分布的特點，因此，在轉(zhuǎn)換中不會造成信息失真.

本文提出的Copula函數(shù)和EEMD相結(jié)合的方法提取故障信號的特征，通過Copula函數(shù)建立起EEMD變換，得到IMF之間的關(guān)聯(lián)性，彌補了IMF進(jìn)行單獨分析的不足.

1 Copula函數(shù)

Copula函數(shù)又稱為連接函數(shù)，它可以構(gòu)建多個隨機(jī)變量之間的聯(lián)合分布函數(shù)，反映了各隨機(jī)變量之間的相關(guān)性.

1.1 Sklar定理

令 F(·，…，·)為具有邊緣分布 F1(·)，F(xiàn)2(·)，…，F(xiàn)n(·)的聯(lián)合分布函數(shù)，則存在一個Copula 函數(shù) C(·，…，·)，滿足［10］

若 F1(·)，F(xiàn)2(·)，…，F(xiàn)n(·)連續(xù)，則 C(·，…，·)唯一確定，其中C(·，…，·)為 Copula函數(shù).

1.2 常用Copula函數(shù)介紹

常用的 Copula函數(shù)有3種類型:橢圓形(Gaussian Copula函數(shù)、t-Copula函數(shù))、二次型和Archimedean型.

二維Gaussian Copula函數(shù)的表達(dá)形式為

式中:ρ∈［－1，1］為函數(shù)的相依系數(shù);Φ、Φ－1分別為標(biāo)準(zhǔn)正態(tài)分布及其反函數(shù).

二維Gaussian Copula函數(shù)的分布密度如圖1所示.Gaussian Copula函數(shù)能描述隨機(jī)變量之間正負(fù)兩方面的相關(guān)性，而且具有參數(shù)估計的計算量小等優(yōu)點，所以被用于特征提?。?］.

圖1 二維Gaussian Copula函數(shù)的分布密度Fig.1 The distribution density function diagram of the two-dimensional Gaussian Copula function

2 聚合經(jīng)驗?zāi)B(tài)分解

EEMD是EMD的改進(jìn)方法，因為EMD存在嚴(yán)重的模態(tài)混疊問題，使得IMF中含有較寬的尺度信號，或者同一尺度的信號包含在不同的IMF中.EEMD在EMD的基礎(chǔ)上，通過對原始信號加入高斯白噪聲，避免了信號的間斷問題，通過對加入高斯白噪聲的原始信號分解的結(jié)果多次求平均，使白噪聲相互抵消，能抑制模態(tài)混疊.

EEMD 算法的步驟如下［11］:

(1)初始化聚合次數(shù)N和高斯白噪聲的幅值，并使計算次數(shù)m=1.

(2)在計算m次時加入的高斯白噪聲，步驟為:

①按給定幅值添加白噪聲序列到信號中，則有

其中，k為加入白噪聲的復(fù)制系數(shù)，nm(t)為第m次添加的白噪聲序列;

②利用EMD將加入白噪聲后的信號xm(t)分解為一組IMFs;

③當(dāng)m＜N(N為經(jīng)驗?zāi)B(tài)分解的聚合次數(shù))時，重復(fù)①和②，每次加入不同的白噪聲信號，并使m=m+1.

(3)計算N次分解出的各個IMFs均值

其中:ci，m為由第m次經(jīng)驗?zāi)B(tài)分解得到的第i個IMF.

(4)選取每個IMF在N次分解的均值作為最終的本征模態(tài)函數(shù)，有

(5)各個本征模態(tài)函數(shù)所包含的頻率成分按照從高到低的順序進(jìn)行排列.

3 Copula函數(shù)與EEMD結(jié)合的特征提取方法

3.1 選擇 IMF

對原始信號進(jìn)行EEMD變換，設(shè)定EEMD的聚合次數(shù)以及白噪聲幅值，系統(tǒng)能根據(jù)信號自適應(yīng)地確定IMF個數(shù)，依次按照頻率成分從高到低的順序排列.

觀察變換所得的各IMF的頻譜圖，根據(jù)列車轉(zhuǎn)向架故障時振動信號的頻率范圍，選擇出其中在故障頻率范圍內(nèi)的IMF，選出其中關(guān)聯(lián)性最大的兩個IMF用作進(jìn)一步研究.本文采用Chi-plot［12］理論研究其中任意兩個的IMF的相關(guān)性.

在Chi-plot中，點n的坐標(biāo)(Xi，Yi)被轉(zhuǎn)換成(λni，χni)，i=1，2，…，n.Chi-plot圖在二維坐標(biāo)中的范圍是［－1，1］×［－1，1］.(λni，χni)偏離水平線χ=0的程度反映了X和Y變量的相關(guān)程度.(λni，χni)的計算公式為

式中:各參數(shù)的具體含義見文獻(xiàn)［12］.

在Chi-plot圖中設(shè)定有兩條水平控制線，作為測度變量X和Y相關(guān)性的置信區(qū)間.文獻(xiàn)［12］給出了當(dāng)置信區(qū)間p為不同值時，其他參數(shù)的對應(yīng)取值.當(dāng)絕大多數(shù)點落在兩條水平控制線范圍內(nèi)時，說明X和Y相互獨立;反之X和Y具有較強(qiáng)的相關(guān)性.

作為研究相關(guān)性的工具，Chi-plot圖直觀、定性地刻畫兩個變量之間的相關(guān)性.

3.2 構(gòu)建Copula函數(shù)邊緣分布函數(shù)

構(gòu)建邊緣分布函數(shù)是使用Copula函數(shù)構(gòu)建聯(lián)合分布函數(shù)的重要步驟.高速列車信號經(jīng)EEMD變換后得到的IMF分量呈現(xiàn)長尾狀的非高斯分布，因此高斯分布不能很好地擬合高速列車信號分布的形狀.

有兩種常用的模型可以擬合高速列車信號:泛化高斯模型(generalized Gaussian distribution，GGD)、核 密 度 估 計 (kernel density estimation，KDE).KDE是非參數(shù)估計模型，計算量很大.經(jīng)實驗驗證，采用 GGD能對 IMF的邊緣分布進(jìn)行擬合.

GGD密度函數(shù)的形式為

式中:

α為尺度參數(shù);

β為形狀參數(shù)，

當(dāng)β=2時，GGD為Gaussian分布，

當(dāng)β=1時，GGD是Laplace分布.

GGD參數(shù)計算方法有最大似然估計和Newton-Raphson算法.本文參數(shù)獲取方式為最大似然估計法.

本文中使用了邊緣分布函數(shù)的kull-back leibler distance(KLD)作為其中的一個特征，GGD具有解析形式的KLD，可用于模式識別中.

假定GGD概率密度函數(shù)分別為fi=p(x;αi，βi)和 fj=p(x;αj，βj)，其 KLD的表達(dá)式為

KLD作為一個重要特征，表征了兩個邊緣分布函數(shù)之間的距離，被廣泛使用在模式識別問題中［8］.

3.3 Copula函數(shù)構(gòu)建聯(lián)合分布及特征提取

使用GGD擬合兩個IMF的邊緣分布函數(shù)后，使用Gaussian Copula函數(shù)構(gòu)建其聯(lián)合分布函數(shù).這是由于Gaussian Copula和t-Copula均能用于描述隨機(jī)變量之間的正負(fù)兩方面的相關(guān)性，而且用最大似然估計法求解t-Copula參數(shù)的計算量遠(yuǎn)大于Gaussian Copula的計算量.

對Copula函數(shù)參數(shù)進(jìn)行估計常用的方法有兩種.一種是完全最大似然估計法，即一次性估計出邊緣分布函數(shù)和Copula函數(shù)的參數(shù);另一種是兩階段最大似然估計法，即先估計出邊緣分布函數(shù)的參數(shù)，再使用最大似然估計法求出Copula函數(shù)的參數(shù)值.本文所用方法為兩階段最大似然估計法，先通過最大似然估計法求得式(5)中的參數(shù)，即邊緣分布函數(shù)的參數(shù)，再通過最大似然估計法求出Gaussian Copula函數(shù)的參數(shù).

使用Gaussian Copula函數(shù)構(gòu)建得到兩個IMF的聯(lián)合分布函數(shù)后，提取聯(lián)合概率密度函數(shù)的均值和方差作為另外的兩個特征.

本文特征提取部分計算量主要在估計Copula多維模型參數(shù)，該復(fù)雜度主要由3部分組成:估計GGD參數(shù)的復(fù)雜度、估計Copula函數(shù)參數(shù)的復(fù)雜度以及估計相關(guān)矩陣R的復(fù)雜度.

設(shè)L為數(shù)據(jù)長度，d為特征向量維數(shù)，根據(jù)文獻(xiàn)［13］，整個特征提取過程的復(fù)雜度為

本文特征提取方法流程如圖2所示.

圖2 特征提取方法流程Fig.2 Flowchart of feature extraction

4 實驗結(jié)果及分析

4.1 數(shù)據(jù)來源

采用動力學(xué)仿真分析的多體動力學(xué)分析軟件包，針對某型號動車組動車轉(zhuǎn)向架的機(jī)械故障進(jìn)行了仿真實驗，得到車體前部橫向加速度、車體中部橫向加速度和車體后部橫向加速度信號.工況分為正常、抗蛇行減振器失效、空氣彈簧失效、橫向減振器失效，速度設(shè)定為 200 km/h，仿真時間為3.6 min，采樣頻率為243 Hz.車體前部橫向加速度在4種工況下的原始信號如圖3所示.

圖3 車體前部橫向加速度在4種工況下的時域信號Fig.3 The train’s front lateral acceleration channel signals under four working conditions

4.2 實驗結(jié)果及分析

對仿真得到的高速列車轉(zhuǎn)向架振動信號，使用小波包濾波對原始信號進(jìn)行預(yù)處理.對預(yù)處理后的信號進(jìn)行EEMD變換，設(shè)置聚合次數(shù)為100，白噪聲幅值為0.2，系統(tǒng)自適應(yīng)選取得到16個IMF.其中前6個IMF的頻譜范圍在列車轉(zhuǎn)向架故障信號的頻率范圍內(nèi).再使用Chi-plot圖分析6個IMF中任意2個IMF的相關(guān)性，設(shè)定置信區(qū)間p的值為0.95，用式(3)、(4)計算并繪制任意兩個IMF之間的Chi-plot圖，得到IMF2與 IMF3的 Chi-plot圖中(λni，χni)偏離 χ=0 程度最遠(yuǎn)，說明 IMF2 與 IMF3之間具有最大相關(guān)性，故使用IMF2以及IMF3作下一步分析.

對車體前部橫向加速度通道，對4種工況下的信號使用Chi-plot對IMF2和IMF3進(jìn)行相關(guān)性檢驗，結(jié)果如圖4所示.

圖4 IMF2與IMF3依賴關(guān)系的Chi-plot圖Fig.4 The Chi-plot diagram of the dependence between IMF2 and IMF3

由圖4可知，4種工況下絕大多數(shù)的點都落在水平控制線以外，IMF2與IMF3之間存在較強(qiáng)相關(guān)性.根據(jù)(λni，χni)偏離水平線 χ=0的程度可以判斷，正常、抗蛇形減振器失效和空氣彈簧失效故障信號經(jīng)EEMD所得的IMF2與IMF3之間的相關(guān)性較大，橫向減振器失效故障信號經(jīng)EEMD所得的IMF2與IMF3之間的相關(guān)性相對其他3種工況較小，但仍存在一定的相關(guān)性.

利用GGD對IMF2和IMF3的邊緣分布進(jìn)行擬合，車輛正常信號經(jīng)過EEMD變換得到IMF2的分布以及使用GGD擬合的結(jié)果如圖5所示.由圖5可知擬合結(jié)果較好.

計算兩個邊緣分布的KLD值作為特征，車體前部橫向加速度通道下4種工況的KLD分布如圖6所示.

由圖6可知，本次實驗提取的KLD值作為特征可將橫向減振器失效及抗蛇行減振器失效與剩下兩種工況加以區(qū)分，說明了提取KLD值作為特征值的合理性.

使用Gaussian Copula函數(shù)構(gòu)建IMF2與IMF3的聯(lián)合分布，并提取聯(lián)合概率密度函數(shù)求均值和方差作為另外兩個特征.車體前部橫向加速度通道下所得4種工況的IMF2與IMF3聯(lián)合概率密度函數(shù)均值如圖7所示.

構(gòu)建得到的聯(lián)合概率密度函數(shù)的方差的分布如圖8所示.

車體前部橫向加速度通道提取得到的KLD、聯(lián)合概率密度函數(shù)的均值、方差的三維特征分布如圖9所示.

由圖9可知，本次實驗所提取的3個特征在三維空間分布上直觀地將4種工況信號加以區(qū)分，說明了本次實驗所提特征的有效性.

實驗數(shù)據(jù)的特征隨機(jī)二等分為訓(xùn)練樣本和測試樣本，選擇高斯核函數(shù)訓(xùn)練支持向量機(jī)(support vector machine，SVM)，用交叉驗證搜索法［14］尋優(yōu)懲罰參數(shù)C和高斯核函數(shù)參數(shù)γ.

圖5 IMF2分布圖以及GGD擬合結(jié)果Fig.5 The distribution of IMF2 and the result fitted by GGD

圖6 4種工況的KLD分布Fig.6 The KLDdistribution of four working conditions

圖7 聯(lián)合概率密度函數(shù)均值的分布Fig.7 The distribution of the mean of the joint probability density function

試驗共提取得到280個樣本，使用50%做訓(xùn)練樣本，50%作測試樣本.將提取的3個特征使用SVM加以分類，隨機(jī)選取訓(xùn)練樣本和測試樣本，進(jìn)行20次實驗.得到車體前部、中部、后部橫向加速度通道在4種工況的分類識別結(jié)果如表1所示.

從表1可知，車體橫向振動工況在車體前部、中部、后部加速度通道上的平均識別率分別為95.97%、97.21%、96.54%，表明 Copula 函數(shù)與EEMD結(jié)合的方法提取得到的特征能很好地對高速列車轉(zhuǎn)向架故障信號進(jìn)行表征.文獻(xiàn)［4］通過對信號進(jìn)行EEMD變換后提取IMF熵值作為特征，所得平均識別率為88%.本文利用IMF分量之間關(guān)聯(lián)性進(jìn)行特征提取的方法，所得平均識別率超過95%，得到了更好的識別效果.

圖8 聯(lián)合概率密度函數(shù)方差的分布Fig.8 The distribution of the variation of the joint probability density function

圖9 4種工況的特征分布Fig.9 The feature distribution of four working conditions

表1 車體3個部位橫向加速度通道識別結(jié)果Tab.1 The recognition results of train's lateral acceleration channel signals of three positions %

5 結(jié)束語

針對高速列車轉(zhuǎn)向架振動信號，提出了一種Copula函數(shù)與EEMD相結(jié)合的特征提取方法，彌補了傳統(tǒng)方法對EEMD變換得到的IMF進(jìn)行單獨分析的不足.通過選取信號經(jīng)過EEMD變換所得的具有最大互相關(guān)性的兩個 IMF分量，并使用GGD對其分布進(jìn)行擬合，使用Gaussian Copula函數(shù)構(gòu)建聯(lián)合分布并提取特征.通過對某型列車的前部、中部以及后部的橫向加速度通道實驗結(jié)果表明，該特征提取方法的識別率均超過95%，能得到較好的識別效果.

［1］王新銳，丁勇，李國順.鐵路火車可靠性實驗方法及評價標(biāo)準(zhǔn)的研究［J］.中國鐵道科學(xué)，2010，31(1):116-122.WANG Xinrui，DING Yong，LI Guoshun.Study on the reliability test method and evaluation specifications of railway freight car［J］.China Railway Science，2010，31(1):116-122.

［2］唐曦凌，梁霖，高慧中，等.結(jié)合小波變換和多約束非負(fù)矩陣分解的故障特征提取方法［J］.振動與沖擊，2013，24(19):7-11.TANG Xiling，LIANG Lin，GAO Huizhong，et al.Fault feature extraction method combining continuous wavelet transformation with multi-consitraint nonnegative matrix factorization［J］.Journal of Vibration and Shock，2013，24(19):7-11.

［3］WU Z H，HUANG N E.Ensemble empirical mode decomposition: a noise-assisted data analysis method［R］. Calverton:Center for Ocean-Land-Atmosphere Studies，2009.

［4］秦娜，金煒東，黃進(jìn)，等.基于EEMD樣本熵的高速列車轉(zhuǎn)向架故障特征提?。跩］.西南交通大學(xué)學(xué)報，2014，49(1):27-32.QIN Na，JIN Weidong，HUANG Jin，et al.Feature extraction of high speed train bogie based on ensemble empirical mode decomposition and sample entropy［J］.Journal of Southwest Jiaotong University，2014，49(1):27-32.

［5］FENG Z P，ZUO M J，HAO R J.Ensemble empirical mode decomposition-based teager energy spectrum for bearing fault diagnosis［J］.J.Vib.Acoust ，2013，135(3):031013-031033.

［6］李悅，程希駿.上證指數(shù)和恒生指數(shù)的Copula尾部相關(guān)性分析［J］.系統(tǒng)工程，2006，24(5):88-92.LI Yue，CHENG Xijun.Tail dependence analysis of SZI＆ HSI based on Copula method［J］. Systems Engineering，2006，24(5):88-92.

［7］郭生練，閆寶偉，肖義，等.Copula函數(shù)在多變量水文分析計算中的應(yīng)用及研究進(jìn)展［J］.水文，2008，28(3):1-7.GUO Shenglian， YAN Baowei， XIAO Yi， et al.Multivariate hydrological analysis and estimation［J］.Journal of China Hydrology，2008，28(3):1-7.

［8］LI Chaorong，LI Jianping，F(xiàn)U Bo.Magnitude-phase of quaternion wavelet transform for texture representation using multilevel Copula［J］.IEEE Signal Processing Letters，2013，20(8):799-802.

［9］LASMAR N E，BERTHOUMIEU Y.Gaussian Copula multivariate modeling for texture image retrieval using wavelet transforms［J］.IEEE Transactions on Image Processing，2012，23(5):2246-2261.

［10］NELSON R B.An introduction to Copulas［M］.New York:Springer，2006:32-105.

［11］ZVOKELJM，ZUPAN S，PREBEL I.Non-liner multivariate and multiscale monitoring and signal denoising strategy using kernel principal component analysis combined with ensemble empiricalmode decomposition method［J］. Mechanical System and Signal Processing，2011，25(7):2631-2653.

［12］FISHER N I，SWITZER P.Chi-plots for assessing dependence［J］.Biometrica，1985，72:253-265.

［13］李朝榮.Copula驅(qū)動的小波域紋理特征提取研究［D］.成都:電子科技大學(xué)，2013.

［14］CHANG Chihchung，LIN Chihjen.LIBSVM:a library for support vector machines［J］.ACM Transactions on Intelligent，Systems and Technology，2011，2(3):1-27.