钢轨纵向波磨是指钢轨投入使用后，轨顶面出现的沿纵向分布的周期性类似波浪形状的不平顺现象，是轨道损伤的一种主要形式^[1].如图 1所示，当车辆行驶至路基不匀的轨道或上下坡时，置于车体底部的激光位移传感器由于车体刚性不能实时跟踪轨道起伏的变化，导致测量的短周期波磨数据中含有轨道起伏引起的长周期趋势项，从而影响波磨测量的精度.因此，提取和消除信号中存在的趋势项在钢轨波磨检测领域扮演着一个重要的角色^[2].

趋势项通常被认为是包含信号全局改变信息的加性平滑成分^[3]，是信号的“骨架”，和信号中其他成分相比具有典型的低复杂度、高幅特点.传统的趋势项去除方法普遍采用最小二乘原理，只能处理平稳且数据量小的信号，且需预先假定信号中趋势项的类型^[4].而实际采集到的趋势项数据多为非平稳的随机信号，采用这种方法得到的结果是不准确的.

聚合经验模态分解(EEMD, ensemble empirical mode decomposition)以其自适应的多尺度分解特性^[5-6]，特别适合分析非线性、非平稳信号，目前已有的EEMD/EMD去趋算法主要从低频本征模式函数(IMF, Int Rinsicmode Function)分量中选取连续或间断的几个组合重构作为信号的趋势，如通过线性规划来选取相应IMF的EMD-SBP^[7]；通过分析相邻IMF边际谱的相关性找到信号与趋势项的分界点，然后选取后几个连续IMF进行组合的EMD-HHT等^[8].这些方法都为EMD去趋提供了很好的思路，取得了较好的去趋效果.但忽视了一些潜在的问题：低频IMF分量中除了占据较大比例能量的趋势项成分外，还可能含有白噪声与信号的低频成分.这三种成分混合在一起，并没有完全隔离开来，因此直接选取低频IMF分量组合趋势是不精确的.

针对上述问题，本文提出一种EEMD-SVD-PE的IMF分量重构方法，其中，奇异值分解是一种信号的特征提取的技术^[9].将IMF分量组成多维特征矩阵投影到奇异值分解(SVD, Singular Value Decomposition)得到的能量特征向量矩阵空间重构出新的奇异值分量，此时奇异值分量按照能量高低进行分布，同时引入的排列组合熵算法是一种衡量一维序列复杂度的平均熵参数^[10]，能够有效的提取出复杂度低的奇异值分量即钢轨波磨趋势项.它解决了传统EMD直接以低频IMF分量组合趋势时存在的信号混杂问题，有效提升了去趋效果.

1 基于PE选取低复杂度奇异值分量重构趋势的EEMD-SVD信号去趋方法

由于钢轨波磨信息中的趋势项具有典型的低复杂度、高幅的特点，表明在提取钢轨波磨趋势项时应提取能量高且复杂度低的分解信号，本文提出的基于PE选取低复杂度奇异值分量重构趋势的EEMD-SVD信号去趋方法，详细步骤如下：

1.1 高斯平滑滤波与EEMD分解

对长度为M的一维含噪时间序列s(t)，首先对其进行高斯平滑，滤除部分随机脉冲噪声对后续分解的影响，得到平滑后的序列$\tilde s$(t), 然后对其进行噪声辅助分析的EEMD分解^[11-12]，其分解过程如下：

1) 在目标数据中加入白噪声序列；

2) 将加入白噪声的序列分解为IMF；

3) 每次加入不同的白噪声序列，反复重复步骤1)、步骤2) N次；

4) 把分解得到的各个IMF对应求平均作为最终的结果.

得到$\tilde s\left( t \right) = \sum\limits_{i = 1}^N {{c_i}\left( t \right) + r\left( t \right)} $.若将余项r(t)也作为IMF分量，则$\tilde s\left( t \right) = \sum\limits_{i = 1}^{N + 1} {{c_i}\left( t \right)} $.

1.2 SVD构建奇异值分量

c₁位于IMF最高频段，所占据的频带宽度最宽，拥有的噪声能量最多，所以将其舍去，其余N个长度为M的IMF分量组成一个多维特征信号矩阵X_M×N=(c₂^T, c₃^T, …, c_N+1^T).将信号矩阵去均值

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;{{\bar X}_{{\rm{M}} \times {\rm{N}}}} = \left( {{{\mathit{\boldsymbol{\bar x}}}_1},{\mathit{\boldsymbol{x}}_2}, \cdots ,{{\mathit{\boldsymbol{\bar x}}}_{\rm{N}}}} \right) = \\ \left( {{x_1} - E\left( {{\mathit{\boldsymbol{x}}_1}} \right),{\mathit{\boldsymbol{x}}_2} - E\left( {{\mathit{\boldsymbol{x}}_2}} \right), \cdots ,{\mathit{\boldsymbol{x}}_N} - E\left( {{\mathit{\boldsymbol{x}}_{\rm{N}}}} \right)} \right). \end{array} $

(1)

同时，希望通过信号的能量对现有信号矩阵进行重构，信号的协方差矩阵的特征值往往可以用来表示信号强度^[13]，所以将原始信号的协方差矩阵C_N×N=E((X^T)_N×M(X)_M×N)作为转换矩阵，通过奇异值分解，

$ {\mathit{\boldsymbol{C}}_{{\rm{N}} \times {\rm{N}}}} = {\mathit{\boldsymbol{U}}_{{\rm{N}} \times {\rm{N}}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{\rm{N}}}_{ \times {\rm{N}}}\mathit{\boldsymbol{U}}_{{\rm{N}} \times {\rm{N}}}^{\rm{T}}. $

(2)

其中: Λ=diag(λ₁, λ₂, …, λ_N)为对角矩阵，λ₁≥λ₂≥…≥λ_N为C的特征值，U=(u₁, u₂, …, u_N)为特征值对应的特征向量组成的正交矩阵，去除特征值为0的特征向量项及原始数据向量，实现数据降维.降维后的奇异值分量个数K(K≤N).重构后的奇异值分量信号为

$ {\mathit{\boldsymbol{P}}_{{\rm{M}} \times {\rm{K}}}} = {\left( {\mathit{\boldsymbol{\bar X}}} \right)_{{\rm{M}} \times {\rm{K}}}}{\mathit{\boldsymbol{U}}_{{\rm{K}} \times {\rm{K}}}} = \left( {{\mathit{\boldsymbol{p}}_1},{\mathit{\boldsymbol{p}}_2}, \cdots ,{\mathit{\boldsymbol{p}}_{\rm{K}}}} \right). $

(3)

降维重构的原始信号为

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{\tilde X = }}{\mathit{\boldsymbol{P}}_{{\rm{M}} \times {\rm{K}}}}\mathit{\boldsymbol{U}}_{{\rm{K}} \times {\rm{K}}}^{\rm{T}} + {\mathit{\boldsymbol{E}}_{{\rm{M}} \times {\rm{K}}}}.\\ {\rm{式中:}}{\mathit{\boldsymbol{E}}_{{\rm{M}} \times {\rm{K}}}} = {\left( {\begin{array}{*{20}{c}} {E\left( {{x_1}} \right)}&{E\left( {{x_2}} \right)}& \cdots &{E\left( {{x_{\rm{K}}}} \right)}\\ {E\left( {{x_1}} \right)}&{E\left( {{x_2}} \right)}& \cdots &{E\left( {{x_{\rm{K}}}} \right)}\\ \vdots & \vdots & \vdots & \vdots \\ {E\left( {{x_1}} \right)}&{E\left( {{x_2}} \right)}& \cdots &{E\left( {{x_{\rm{K}}}} \right)} \end{array}} \right)_{{\rm{M}} \times {\rm{K}}}} \end{array} $

(4)

为M个样本的均值扩展矩阵.记为P_M×K=(p₁, p₂, … p_K)

1.3 基于排列组合熵选取低复杂度奇异值分量

第i个奇异值分量p_i的排列组合熵为

$ H\left( {{\mathit{\boldsymbol{p}}_i}} \right) = - \sum\limits_{{\pi _j} \in {S_{\rm{D}}}} {{\mathit{\boldsymbol{p}}_i}\left( {{\pi _j}} \right){\rm{ln}}\left( {{\mathit{\boldsymbol{p}}_i}\left( {{\pi _j}} \right)} \right).} $

(5)

式中:π_j是奇异值分量的任意一种排列组合方式，对各个奇异值分量p_i(i=1, 2, …, K)的排列组合熵按从大到小排列得到矢量Η=(h₁, h₂, …, h_K), h₁≥h₂≥…≥h_K，由排列组合熵的定义可以看出：排列组合熵的大小反应序列的不规则程度，排列组合熵值越大，序列越不规则，值越小越规则.将得到的K个奇异值分量的排列组合熵进行k-means聚类，将聚类中的排列组合熵最小的一类作为选定的奇异值分量.由此可得所有奇异值分量的趋势项选取权系数矢量w=(w₁, w₂, …, w_K).

1.4 依据选取的低复杂度奇异值分量重建趋势项

依据趋势项选取权系数矢量来选取相应奇异值分量进行重构，再与扩展均值矩阵结合，得到原始信号的趋势项矩阵，即

$ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{{\rm{M}} \times {\rm{K}}}} = \sum\limits_{i\left| {{w_i} = 1} \right.} {{\mathit{\boldsymbol{p}}_i}{{\left( {{u_i}} \right)}^{\rm{T}}} + {\mathit{\boldsymbol{E}}_{{\rm{M}} \times {\rm{K}}}}.} $

(6)

对趋势项矩阵按列求和，即得原始信号的一维趋势项.流程见图 2所示.

2 仿真信号去趋效果对比分析

为验证本文所提方法的有效性，参照文献[14]引入一种仿真信号进行测试，并分别采用无相位延迟的低通滤波法(LF)、与EMD结合的线性规划法(EMD-SBP)^[7]、小波分解的方法(WD)^[15]和本文提出的基于排列组合熵选取低复杂度奇异值分量重构趋势的EEMD-SVD法(EEMD-SVD-PE)进行对比.同时，使用均方根误差(RMSE)，相对L₂范数来定量分析各种方法的去趋效果^[16]，它的表达式为

$ {\rm{RMSE = }}\sqrt {\sum\limits_{t = 1}^M {\frac{1}{M}{{\left( {T\left( t \right) - \Gamma \left( t \right)} \right)}^2}} } , $

(7)

$ {L_2} = \frac{{\sqrt {\sum\limits_{t = 1}^M {{{\left( {\Gamma \left( t \right) - T\left( t \right)} \right)}^2}} } }}{{\sqrt {\sum\limits_{t = 1}^M {{{\left( {T\left( t \right)} \right)}^2}} } }}. $

(8)

式中:Γ(t)为算法提取出的趋势项，T(t)为原始信号的趋势项，M为测试采样长度.RMSE和相对L₂范数越小，表示提取出的趋势越接近原始信号趋势，去趋效果越好.

仿真信号为s(t)=y(t)+T(t)+n(t)，式中y(t)=A₁sin(2πf₁t)，A₁=0.6，f₁=100 Hz，T(t)=A₂e^－bt+A₃sin(2πf₂t)+A₄cos(2πf₃t)，A₂=4，b=0.02，A₃=2，f₂=1 Hz，A₄=2.5，f₃=2 Hz，n(t)为加性高斯白噪声成分，一段长为1 000的信号s(t)见图 3所示.为验证该方法在不同噪声环境下去趋效果的鲁棒性，分别进行单信噪比和多信噪比测试.

2.1 单信噪比测试

令信噪比(SNR, Signal-Noise Ratio)为5 dB，EEMD分解得到的各层IMF分量与SVD分解得到的奇异值分量分别见图 4(a)、(b)所示，各层奇异值分量的PE权系数见表 1所示.可看出，经过k-means分类后取最小排列组合熵的一类，作为趋势项的奇异值分量为p1、p2和p3.这4种方法提取出的趋势项与原趋势项对比效果见图 5所示，低通滤波法(通带截止频率取2 Hz，阻带截止频率取5 Hz)提取的趋势项信号幅度不足；与EMD结合的线性规划法(实线，选取第7层IMF作为趋势项) IMF成分选取不足，提取出的趋势项与原趋势项(粗点线)在峰值处相差较大；WD(细点线，小波分解后的第一项)在t=0.05与t=0.95信号边界处与原趋势相差较大；本文方法(虚线，选取第1、2、3项奇异值分量重构组合趋势)趋势提取的更为充分.表 2为4种方法提取出的趋势与原趋势的RMSE和相对范数比较，本文方法相对其他三种，RMSE分别约降低了0.183 4，0.254 9和0.055 1，相对范数分别约降低了0.112 6，0.127和0.028 3.

2.2 多信噪比测试

令R_SN=-5~20 dB，步长1 dB，每个R_SN等级执行1 000次试验，每次加入不同的白噪声序列，统计每个等级处的提趋准确率.4种方法的提趋准确率-噪声等级曲线见图 6所示.可以看出，随着信噪比增加，与EMD结合的线性规划法、直接选取小波分解项的WD法和本文方法提趋准确率均逐步增高，最高达到98%；低通滤波法在R_SN=10 dB后准确率逐步下降，原因是低通滤波在复杂信号中效果较差.相比而言，本文所提方法提趋准确率保持稳步上升，总体性能优于其他三种，尤其在低信噪比情况下，准确率最大提高约30%.

在不同噪声等级下4种方法的相对L₂范数统计均值见表 3.由表 3可知，本文方法在不同噪声等级下的结果均优于其他方法，比低通滤波法(LF)平均降低约0.101，比与EMD结合的线性规划法(EMD-SBP)降低约0.115，比WD降低约0.015.

3 EEMD-SVD-PE在钢轨波磨数据去趋中的应用 3.1 数据采集

实验中采用三点弦测法对钢轨波磨进行检测^[17]，激光位移传感器的采样间隔设为2 mm.检测原理见图 7(a).

检测流程如下：

1) 车轴转动带动光电编码器旋转输出触发信号，触发信号发给检测同一侧钢轨的3个1D激光位移传感器，完成同一时刻的钢轨波磨数据采集.

2) 数据经数据采集卡进行数据集中，然后传输给车载综合处理计算机.综合处理计算机利用三点弦测法将含有趋势的波磨信息提取、存储.

本文利用实验搭建的钢轨动态磨耗检测平台在实际50钢轨上进行现场实验见图 7(b)所示.实验中的钢轨由于长时间的使用，路基呈现出较大的起伏，十分适合本文的情况.通过三点弦测法最终得到17 m含有趋势的钢轨波磨数据.同时，为了得到标准的波磨信息，用波磨尺对同一路段的钢轨波磨进行测量，并以测得的数据作为去趋后的标准数据，见图 7(c).采集到的含有趋势项的钢轨波磨数据见图 8.

3.2 钢轨波磨去趋

现场实验采集得到的波磨数据由于趋势项存在，轨道表面高低起伏很大，起伏范围达到±0.6 mm.EEMD分解得到的各层IMF分量与SVD分解得到的奇异值分量分别见图 9(a)、(b)所示，各层奇异值分量的PE权系数见表 4所示.此时选取p1、p3、和p6项奇异值分量来重构趋势，同时也可以看出组成趋势项的奇异值分量并未连续排列.

4种方法提取出的趋势对比见图 10所示，与EMD结合的线性规划法(EMD-SBP，实线)与WD(点线)提趋出的趋势项存在明显的邻近和边界信息丢失与峰值幅度衰减问题，低通滤波法(点划线)在复杂趋势项信号提取中效果较差.相比之下，EEMD-SVD-PE(虚线)提取出的复杂趋势成分更好地反映了信号走势，边界信息、临近信息和峰值信息得到了较好的保留.与其他3种方法相比的结果见表 5，本文方法的RMSE分别降低了0.086 0，0.069 1，0.042 8.应用本文所提方法去趋后的钢轨波磨效果见图 11所示，其中红色为由波磨尺测得的真实波磨值，测得的波磨与真实波磨的趋势相近，且最大误差在±0.2 mm的范围之内，符合实际需要.

实验在matlabR2014a环境下进行，计算机内存为3.39 G，CPU主频为3.2 GHz.统计仿真信号的平均去趋时间开销：低通滤波法(LF)约为0.002 908 s，与EMD结合的线性规划法(EMD-SBP)为0.012 04 s，WD为0.010 528 s，本文提出的EEMD-SVD-PE为0.010 914 s，时间开销没有明显增加.

4 结论

针对现有的EMD去趋方法直接选取低频IMF组合趋势项时存在的IMF中信号成分混杂问题，本文提出了一种基于排列组合熵选取低复杂度奇异值分量重构趋势的EEMD-SVD信号去趋新方法.首先对信号进行高斯平滑去噪和EEMD分解，然后选取第二项至余项IMF分量组成多维特征矩阵投影到SVD得到的能量特征向量矩阵空间重构出新的奇异值分量，最后基于排列组合熵自适应选取相应低复杂度奇异值分量组合成信号趋势项.为了验证所提算法的性能，引入了仿真趋势项信号，并选择了低通滤波法、与EMD结合的线性规划法、小波分解WD进行去趋效果比较.单信噪比和多信噪比的仿真实验，均表明所提算法的去趋效果总体优于其他三种，尤其在低信噪比情况下，提趋准确率最大提高约30%.然后对含趋轨道波磨数据进行的去趋实例应用，也进一步验证了本文方法的有效性.