金煒東， 呂乾勇， 孫永奎

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

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

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

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

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

1 Copula函數(shù)

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

1.1 Sklar定理

令 F(·，…，·)為具有邊緣分布 F1(·)，F(xiàn)2(·)，…，F(xiàn)n(·)的聯(lián)合分布函數(shù)，則存在一個(gè)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ù)估計(jì)的計(jì)算量小等優(yōu)點(diǎn)，所以被用于特征提?。?］.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 選擇 IMF

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

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

在Chi-plot中，點(diǎn)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)的計(jì)算公式為

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

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

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

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

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

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

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

式中:

α為尺度參數(shù);

β為形狀參數(shù)，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實(shí)驗(yàn)結(jié)果及分析

4.1 數(shù)據(jù)來(lái)源

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

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

4.2 實(shí)驗(yàn)結(jié)果及分析

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

對(duì)車體前部橫向加速度通道，對(duì)4種工況下的信號(hào)使用Chi-plot對(duì)IMF2和IMF3進(jìn)行相關(guān)性檢驗(yà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ù)的點(diǎn)都落在水平控制線以外，IMF2與IMF3之間存在較強(qiáng)相關(guān)性.根據(jù)(λni，χni)偏離水平線 χ=0的程度可以判斷，正常、抗蛇形減振器失效和空氣彈簧失效故障信號(hào)經(jīng)EEMD所得的IMF2與IMF3之間的相關(guān)性較大，橫向減振器失效故障信號(hào)經(jīng)EEMD所得的IMF2與IMF3之間的相關(guān)性相對(duì)其他3種工況較小，但仍存在一定的相關(guān)性.

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

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

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

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

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

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

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

實(shí)驗(yàn)數(shù)據(jù)的特征隨機(jī)二等分為訓(xùn)練樣本和測(cè)試樣本，選擇高斯核函數(shù)訓(xùn)練支持向量機(jī)(support vector machine，SVM)，用交叉驗(yàn)證搜索法［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

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

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

圖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個(gè)部位橫向加速度通道識(shí)別結(jié)果Tab.1 The recognition results of train's lateral acceleration channel signals of three positions %

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

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

［1］王新銳，丁勇，李國(guó)順.鐵路火車可靠性實(shí)驗(yàn)方法及評(píng)價(jià)標(biāo)準(zhǔn)的研究［J］.中國(guó)鐵道科學(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］.振動(dòng)與沖擊，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é)報(bào)，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ù)在多變量水文分析計(jì)算中的應(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òng)的小波域紋理特征提取研究［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.