基于Copula函数的高速列车转向架故障特征提取

金炜东， 吕乾勇， 孙永奎

(西南交通大学电气工程学院，四川 成都 610031)

由列车转向架故障引起的振动严重影响列车的舒适性和安全性［1］.传统列车故障诊断的方法有小波变换、聚合经验模态分解(ensemble empirical mode decomposition，EEMD)、功率谱、能量谱等，小波变换［2］和 EEMD［3］方法是较为经典的方法.小波变换方法是通过对故障信号进行分解得到各个频段的信号，然后提取各个频段的信号，用熵值、能量等作为特征诊断故障.EEMD方法是通过对原始信号进行分解得到各个本征模态函数(intrinsic mode function，IMF)，这些本征模态函数表征了原始信号在各个频段下的分量.文献［4］中对原始信号作EEMD变换得到各个阶段的IMF，通过对每个IMF进行单独分析，提取各阶段IMF的熵值作为特征.文献［5］中对原信号进行EEMD变换，选取与原始信号相关性较大的IMF，从选出的IMF中提取Teager能量谱作为特征用于轴承故障诊断.此类方法虽然能对分解得到的各个阶段的IMF进行分析，但由于只对各阶段IMF进行了单独分析，不能捕捉到IMF之间的依赖关系.

Copula函数作为研究随机变量之间关联性的方法，近年来被广泛应用于金融市场分析［6］、水文分析［7］、纹理图像识别［8-9］等领域.采用 Copula 函数的关联性研究列车故障诊断方法能克服传统方法对各个IMF进行单独研究的不足.

相对其他分析相关性的方法，Copula函数具有两个明显的优点［7］.首先，灵活且求解简单，构造多维联合分布函数时可以分为两个步骤完成:构建边缘分布以及选择合适的Copula函数构建联合分布;其次，不要求相同的边缘分布，使得每种边缘分布函数保留了自己分布的特点，因此，在转换中不会造成信息失真.

本文提出的Copula函数和EEMD相结合的方法提取故障信号的特征，通过Copula函数建立起EEMD变换，得到IMF之间的关联性，弥补了IMF进行单独分析的不足.

1 Copula函数

Copula函数又称为连接函数，它可以构建多个随机变量之间的联合分布函数，反映了各随机变量之间的相关性.

1.1 Sklar定理

令 F(·，…，·)为具有边缘分布 F1(·)，F2(·)，…，Fn(·)的联合分布函数，则存在一个Copula 函数 C(·，…，·)，满足［10］

若 F1(·)，F2(·)，…，Fn(·)连续，则 C(·，…，·)唯一确定，其中C(·，…，·)为 Copula函数.

1.2 常用Copula函数介绍

常用的 Copula函数有3种类型:椭圆形(Gaussian Copula函数、t-Copula函数)、二次型和Archimedean型.

二维Gaussian Copula函数的表达形式为

式中:ρ∈［－1，1］为函数的相依系数;Φ、Φ－1分别为标准正态分布及其反函数.

二维Gaussian Copula函数的分布密度如图1所示.Gaussian Copula函数能描述随机变量之间正负两方面的相关性，而且具有参数估计的计算量小等优点，所以被用于特征提取［9］.

图1 二维Gaussian Copula函数的分布密度Fig.1 The distribution density function diagram of the two-dimensional Gaussian Copula function

2 聚合经验模态分解

EEMD是EMD的改进方法，因为EMD存在严重的模态混叠问题，使得IMF中含有较宽的尺度信号，或者同一尺度的信号包含在不同的IMF中.EEMD在EMD的基础上，通过对原始信号加入高斯白噪声，避免了信号的间断问题，通过对加入高斯白噪声的原始信号分解的结果多次求平均，使白噪声相互抵消，能抑制模态混叠.

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

(1)初始化聚合次数N和高斯白噪声的幅值，并使计算次数m=1.

(2)在计算m次时加入的高斯白噪声，步骤为:

①按给定幅值添加白噪声序列到信号中，则有

其中，k为加入白噪声的复制系数，nm(t)为第m次添加的白噪声序列;

②利用EMD将加入白噪声后的信号xm(t)分解为一组IMFs;

③当m＜N(N为经验模态分解的聚合次数)时，重复①和②，每次加入不同的白噪声信号，并使m=m+1.

(3)计算N次分解出的各个IMFs均值

其中:ci，m为由第m次经验模态分解得到的第i个IMF.

(4)选取每个IMF在N次分解的均值作为最终的本征模态函数，有

(5)各个本征模态函数所包含的频率成分按照从高到低的顺序进行排列.

3 Copula函数与EEMD结合的特征提取方法

3.1 选择 IMF

对原始信号进行EEMD变换，设定EEMD的聚合次数以及白噪声幅值，系统能根据信号自适应地确定IMF个数，依次按照频率成分从高到低的顺序排列.

观察变换所得的各IMF的频谱图，根据列车转向架故障时振动信号的频率范围，选择出其中在故障频率范围内的IMF，选出其中关联性最大的两个IMF用作进一步研究.本文采用Chi-plot［12］理论研究其中任意两个的IMF的相关性.

在Chi-plot中，点n的坐标(Xi，Yi)被转换成(λni，χni)，i=1，2，…，n.Chi-plot图在二维坐标中的范围是［－1，1］×［－1，1］.(λni，χni)偏离水平线χ=0的程度反映了X和Y变量的相关程度.(λni，χni)的计算公式为

式中:各参数的具体含义见文献［12］.

在Chi-plot图中设定有两条水平控制线，作为测度变量X和Y相关性的置信区间.文献［12］给出了当置信区间p为不同值时，其他参数的对应取值.当绝大多数点落在两条水平控制线范围内时，说明X和Y相互独立;反之X和Y具有较强的相关性.

作为研究相关性的工具，Chi-plot图直观、定性地刻画两个变量之间的相关性.

3.2 构建Copula函数边缘分布函数

构建边缘分布函数是使用Copula函数构建联合分布函数的重要步骤.高速列车信号经EEMD变换后得到的IMF分量呈现长尾状的非高斯分布，因此高斯分布不能很好地拟合高速列车信号分布的形状.

有两种常用的模型可以拟合高速列车信号:泛化高斯模型(generalized Gaussian distribution，GGD)、核 密 度 估 计 (kernel density estimation，KDE).KDE是非参数估计模型，计算量很大.经实验验证，采用 GGD能对 IMF的边缘分布进行拟合.

GGD密度函数的形式为

式中:

α为尺度参数;

β为形状参数，

当β=2时，GGD为Gaussian分布，

当β=1时，GGD是Laplace分布.

GGD参数计算方法有最大似然估计和Newton-Raphson算法.本文参数获取方式为最大似然估计法.

本文中使用了边缘分布函数的kull-back leibler distance(KLD)作为其中的一个特征，GGD具有解析形式的KLD，可用于模式识别中.

假定GGD概率密度函数分别为fi=p(x;αi，βi)和 fj=p(x;αj，βj)，其 KLD的表达式为

KLD作为一个重要特征，表征了两个边缘分布函数之间的距离，被广泛使用在模式识别问题中［8］.

3.3 Copula函数构建联合分布及特征提取

使用GGD拟合两个IMF的边缘分布函数后，使用Gaussian Copula函数构建其联合分布函数.这是由于Gaussian Copula和t-Copula均能用于描述随机变量之间的正负两方面的相关性，而且用最大似然估计法求解t-Copula参数的计算量远大于Gaussian Copula的计算量.

对Copula函数参数进行估计常用的方法有两种.一种是完全最大似然估计法，即一次性估计出边缘分布函数和Copula函数的参数;另一种是两阶段最大似然估计法，即先估计出边缘分布函数的参数，再使用最大似然估计法求出Copula函数的参数值.本文所用方法为两阶段最大似然估计法，先通过最大似然估计法求得式(5)中的参数，即边缘分布函数的参数，再通过最大似然估计法求出Gaussian Copula函数的参数.

使用Gaussian Copula函数构建得到两个IMF的联合分布函数后，提取联合概率密度函数的均值和方差作为另外的两个特征.

本文特征提取部分计算量主要在估计Copula多维模型参数，该复杂度主要由3部分组成:估计GGD参数的复杂度、估计Copula函数参数的复杂度以及估计相关矩阵R的复杂度.

设L为数据长度，d为特征向量维数，根据文献［13］，整个特征提取过程的复杂度为

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

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

4 实验结果及分析

4.1 数据来源

采用动力学仿真分析的多体动力学分析软件包，针对某型号动车组动车转向架的机械故障进行了仿真实验，得到车体前部横向加速度、车体中部横向加速度和车体后部横向加速度信号.工况分为正常、抗蛇行减振器失效、空气弹簧失效、横向减振器失效，速度设定为 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 实验结果及分析

对仿真得到的高速列车转向架振动信号，使用小波包滤波对原始信号进行预处理.对预处理后的信号进行EEMD变换，设置聚合次数为100，白噪声幅值为0.2，系统自适应选取得到16个IMF.其中前6个IMF的频谱范围在列车转向架故障信号的频率范围内.再使用Chi-plot图分析6个IMF中任意2个IMF的相关性，设定置信区间p的值为0.95，用式(3)、(4)计算并绘制任意两个IMF之间的Chi-plot图，得到IMF2与 IMF3的 Chi-plot图中(λni，χni)偏离 χ=0 程度最远，说明 IMF2 与 IMF3之间具有最大相关性，故使用IMF2以及IMF3作下一步分析.

对车体前部横向加速度通道，对4种工况下的信号使用Chi-plot对IMF2和IMF3进行相关性检验，结果如图4所示.

图4 IMF2与IMF3依赖关系的Chi-plot图Fig.4 The Chi-plot diagram of the dependence between IMF2 and IMF3

由图4可知，4种工况下绝大多数的点都落在水平控制线以外，IMF2与IMF3之间存在较强相关性.根据(λni，χni)偏离水平线 χ=0的程度可以判断，正常、抗蛇形减振器失效和空气弹簧失效故障信号经EEMD所得的IMF2与IMF3之间的相关性较大，横向减振器失效故障信号经EEMD所得的IMF2与IMF3之间的相关性相对其他3种工况较小，但仍存在一定的相关性.

利用GGD对IMF2和IMF3的边缘分布进行拟合，车辆正常信号经过EEMD变换得到IMF2的分布以及使用GGD拟合的结果如图5所示.由图5可知拟合结果较好.

计算两个边缘分布的KLD值作为特征，车体前部横向加速度通道下4种工况的KLD分布如图6所示.

由图6可知，本次实验提取的KLD值作为特征可将横向减振器失效及抗蛇行减振器失效与剩下两种工况加以区分，说明了提取KLD值作为特征值的合理性.

使用Gaussian Copula函数构建IMF2与IMF3的联合分布，并提取联合概率密度函数求均值和方差作为另外两个特征.车体前部横向加速度通道下所得4种工况的IMF2与IMF3联合概率密度函数均值如图7所示.

构建得到的联合概率密度函数的方差的分布如图8所示.

车体前部横向加速度通道提取得到的KLD、联合概率密度函数的均值、方差的三维特征分布如图9所示.

由图9可知，本次实验所提取的3个特征在三维空间分布上直观地将4种工况信号加以区分，说明了本次实验所提特征的有效性.

实验数据的特征随机二等分为训练样本和测试样本，选择高斯核函数训练支持向量机(support vector machine，SVM)，用交叉验证搜索法［14］寻优惩罚参数C和高斯核函数参数γ.

图5 IMF2分布图以及GGD拟合结果Fig.5 The distribution of IMF2 and the result fitted by GGD

图6 4种工况的KLD分布Fig.6 The KLDdistribution of four working conditions

图7 联合概率密度函数均值的分布Fig.7 The distribution of the mean of the joint probability density function

试验共提取得到280个样本，使用50%做训练样本，50%作测试样本.将提取的3个特征使用SVM加以分类，随机选取训练样本和测试样本，进行20次实验.得到车体前部、中部、后部横向加速度通道在4种工况的分类识别结果如表1所示.

从表1可知，车体横向振动工况在车体前部、中部、后部加速度通道上的平均识别率分别为95.97%、97.21%、96.54%，表明 Copula 函数与EEMD结合的方法提取得到的特征能很好地对高速列车转向架故障信号进行表征.文献［4］通过对信号进行EEMD变换后提取IMF熵值作为特征，所得平均识别率为88%.本文利用IMF分量之间关联性进行特征提取的方法，所得平均识别率超过95%，得到了更好的识别效果.

图8 联合概率密度函数方差的分布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个部位横向加速度通道识别结果Tab.1 The recognition results of train's lateral acceleration channel signals of three positions %

5 结束语

针对高速列车转向架振动信号，提出了一种Copula函数与EEMD相结合的特征提取方法，弥补了传统方法对EEMD变换得到的IMF进行单独分析的不足.通过选取信号经过EEMD变换所得的具有最大互相关性的两个 IMF分量，并使用GGD对其分布进行拟合，使用Gaussian Copula函数构建联合分布并提取特征.通过对某型列车的前部、中部以及后部的横向加速度通道实验结果表明，该特征提取方法的识别率均超过95%，能得到较好的识别效果.

［1］王新锐，丁勇，李国顺.铁路火车可靠性实验方法及评价标准的研究［J］.中国铁道科学，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］唐曦凌，梁霖，高慧中，等.结合小波变换和多约束非负矩阵分解的故障特征提取方法［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］秦娜，金炜东，黄进，等.基于EEMD样本熵的高速列车转向架故障特征提取［J］.西南交通大学学报，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］李悦，程希骏.上证指数和恒生指数的Copula尾部相关性分析［J］.系统工程，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函数在多变量水文分析计算中的应用及研究进展［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，FU 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驱动的小波域纹理特征提取研究［D］.成都:电子科技大学，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.