基于混合插值与多项式修正的局部特征尺度分解端点效应抑制方法

张 征， 崔伟成

(1. 鲁东大学 山东 烟台 264025； 2. 海军航空大学 岸防兵学院 山东 烟台 264001)

在基于振动信号分析的旋转机械故障诊断领域，时频分析技术的研究是一个热点[1]。其中，经验模态分解(Empirical Mode Decomposition，EMD)是一种有效的自适应方法，但使用过程中存在过/欠包络、频率混淆、端点效应等问题[2-4]。程军圣、杨宇等在本征时间尺度分解 (Intrinsic Time-scale Decomposition, ITD)[5]的基础上创造性地提出了局部特征尺度分解(Local Characteristic-scale Decomposition，LCD)[6-7]。该方法开拓了自适应时频分析方法的新思路，在旋转机械故障诊断领域得到了广泛的应用[8-10]。

在应用过程中，LCD仍存在端点效应：两侧端点附近的均值曲线出现幅值发散。虽然LCD的端点效应不如EMD严重，但对于低频、短数据序列来说，LCD的端点效应往往不能忍受。

为了完善LCD方法，提高分解质量，进而准确地进行故障诊断，给出了混合插值与多项式修正相结合的LCD端点效应抑制方法。以仿真信号和轴承内圈故障振动信号的分析，验证了方法的有效性。

1 局部特征尺度分解

LCD方法将一个信号分解成若干个内禀尺度分量(Intrinsic Mode Components，ISC)之和。每个ISC分量必须满足2个条件：

(1)任意2个相邻的极值点符号互异。

(2)考虑所有的极值点Xk及对应的时刻τk(k= 1，2，…，M，其中M为极值点的个数)。取两个相距最近且符号相同的极值点(τk,Xk)和(τk+2,Xk+2)，按照式(1)定义一个τk+1时刻的函数值

(1)

则Ak+1与极值点Xk+1的比值关系不变，即

aAk+1+(1-a)Xk+1=0

(2)

式中：a∈(0,1)为一常量，典型地a=0.5。

为了筛选出ISC，LCD方法构造了信号的均值曲线(基线)，将均值曲线不断从原信号中分离，其迭代过程与EMD方法类似，不同之处在于均值曲线的构造方法。

LCD方法用均值点构造均值曲线，均值点的定义为

Lk+1=aAk+1+(1-a)Xk+1

(3)

式中：Lk+1为均值点，在迭代结束时数值为0。

LCD方法采用分段线性方法或三次样条方法由均值点插值均值曲线。经仿真发现，三次样条插值可得到更好的分解结果，本文采用三次样条LCD作为标准LCD。

2 LCD的端点效应抑制方法

端点效应(End-effect)是Huang等在研究EMD的筛分过程时发现的，指的是：在EMD的运算过程中，反复运用三次样条方法根据极值点对上、下包络进行插值，得到的包络曲线在两侧的端点附近出现幅值发散的现象。LCD方法由极值点直接估计均值点，然后插值得到均值曲线，在插值的过程中同样也存在端点效应。端点效应产生的根本原因是端点处的取值不是均值点，使插值得到的均值曲线在端点附近出现误差。为了抑制LCD的端点效应，可预测(或延拓)端点处的均值点取值，也可以将端点处的插值误差控制在一定范围内。本文给出的端点效应抑制方法采用简单的极值点延拓方法，大致估计出端点处的均值点，采用混合插值与多项式修正相结合的插值方法，通过将插值误差控制在端点附近的方法抑制端点效应。

本文方法中，对端点处均值点的估计精度要求不高，而极值点延拓法简单、实用，更为贴近工程，因此采用极值点延拓方法。

LCD方法在筛分的过程中要反复调用插值函数求取均值曲线，要求插值曲线满足：①足够光滑；②在均值点处的值尽量逼近理论值。从实际的效果来看，线性插值的端点效应最小，均值曲线在均值点的值等于局部均值，但线性插值不能逼近真实均值曲线，导致分量失真较大；三次样条插值的端点效应较为明显，但中间段分量失真小；格林样条插值[11]、B样条插值[12]则是在光滑性与均值点误差之间做了取舍，效果介于线性插值与三次样条插值之间，从分解的最终结果来看，较三次样条插值并无明显优势。

基于上述考虑，给出了混合插值的方法：在插值过程中，对端点效应明显的左区间(第一个均值点与第二个均值点之间)、右区间(倒数第二个均值点与最后一个均值点之间)的采用线性方法插值均值曲线；将插值结果与剩余的均值点序列合并，形成新的均值点序列；对新的均值点序列按照三次样条方法插值形成均值曲线。

经混合插值后，均值曲线在中间区间(第二个均值点与倒数第二个均值点之间)的插值曲线保持了足够的光滑性；左区间、右区间的插值曲线无端点效应，但与实际的均值曲线相比，波形失真较大。采用多项式修正的方法解决线性插值带来的失真问题：①对左区间，在第二个均值点向右沿均值曲线取若干数据，以第一个极值点为约束，采用三次多项式修正均值曲线；②对右区间，在倒数第二个均值点向左沿均值曲线取若干数据，以倒数第一个均值点为约束，采用三次多项式获得修正曲线。

基于混合插值及多项式修正的LCD方法步骤为：

(1)将极值点Xk及对应时刻τk序列进行延拓：

(4)

(5)

得到：Xk、τk(k=0，1，…，M+1)。

(2)计算局部均值点序列L(τk)，(k=1，2，…，M)。

(3)用线性方法插值均值点序列L(τk)，(k=1，2，M-1，M)得到左、右区间的均值曲线l1(t),t∈[τ1,τ2]∪t∈[τM-1,τM]。

(4)将l1(t)作为序列，与中间区间的均值点序列L(τk)，(k=2，,3，…，M-2，M-1)合并，用三次样条方法计算均值曲线，截取中间段的均值曲线l2(t),t∈[τ3,τM-2]。

(5)采用三次多项式插值的方法获得左、右区间的均值曲线l3(t),t∈[τ1,τ2]∪t∈[τM-1,τm]。

(6)将l2(t)、l3(t)合并得到的均值曲线，按照标准方法的迭代过程进行LCD分解。

3 改进算法仿真分析

为了考察混合插值及多项式修正方法对LCD端点效应的抑制能力，考察如下信号：

(6)

式中：x(t)为仿真信号，由调幅调频信号x1(t)、正弦信号x2(t)合成；采样频率为fs=1 024 Hz，仿真时间t∈[0,1]。仿真信号的时域波形见图1。

图1 仿真信号时域波形Fig. 1 Simulated signal in time domain waveform

对仿真信号按照标准LCD方法进行分解，得到两个ISC分量和一个余量，将其时域波形列于图2。

图2 标准LCD分解结果时域波形Fig.2 Results of standard LCD in time domain waveform

从分解结果来看：分量1对应调幅调频信号x1(t)，分量2对应正弦信号x2(t)，分解结果与实际曲线总体趋势一致，但是两个分量的端点处均存在失真。其中，分量2的失真最为严重，第一个极点的幅值只有理论值的一半。残余信号的左、右端点处数值接近1.5，大于信号的幅值。

对仿真信号应用改进LCD进行分解，将结果的时域波形列于图3。

图3 改进LCD分解结果时域波形Fig.3 Results of improved LCD in time domain waveform

从改进LCD的分解结果可以看出：

(1)中间区间的分量保持了原始信号的光滑性，与原始信号非常接近。

(2)左区间、右区间的分量与原始信号在波形上的差异已非常小；

(3)改进LCD方法产生的剩余信号的幅值远小于标准LCD。

为了进一步分析效果，借鉴镜像延拓法[14]、支持向量机(Support Vector Machine，SVM)数据预测法[13]处理LCD的端点效应，并与B样条LCD[12]及本文方法其进行对比。采用正交指数(Index of Orthogonality，IO)、能量相对误差(Energy Relative Error，ERE)和计算时间(Calculate time，CT)评价各种方法的性能。

IO定义[15]为：

(7)

IO作为各ISC分量正交性的判定标准，其值越小，说明分量间的正交性越好，反之亦然。

ERE的定义为：

(8)

ERE指标越接近于0，说明分解能量损失越低。

运用四种方法对仿真信号进行分解，将三种指标列于表1。可以看出：BLCD及镜像延拓LCD的IO及ERE指标均较本文方法差，计算时间基本相当；SVM-LCD方法可得到很好的分解效果，本文方法的效果稍差。但SVM结构复杂、参数不易调整，且运行时间过长，不便于工程应用。本文方法在分解质量及运行时间上的综合优势比较明显。

表1 不同方法对仿真信号进行LCD分解的指标比较

混合插值方法能避免两个端点邻域的插值“过冲”现象，经过多项式方法对线性插值进行修正，整体上保证了插值曲线与原始信号的相似程度。因此，混合插值及多项式修正方法能有效抑制插值带来的端点效应，是一种有效的方法。

4 试验数据分析

采用美国凯斯西储大学电气工程实验室的轴承振动数据进行算法验证[16]。其试验装置构造见图4。

图4 试验装置构造图Fig.4 The basic layout of the test rig

一台功率为1.47 kW(2马力)的感应电动机驱动转轴，驱动轴上安装一个扭矩传感器和编码器，扭矩通过测功机和电子控制系统作用于驱动轴，轴的驱动端和风扇端均安装了轴承。论文选用驱动端故障轴承数据集IR007_0.mat中的X105_FE_time作为试验数据，该数据集的故障轴承为SKF公司的6205-2RS JEM深沟球轴承，其尺寸参数如表2如示。

表2 滚动轴承6205-2RS的尺寸参数

使用电火花加工技术在该轴承内圈上布置了单点故障，故障直径为0.177 8 mm、深度为0.279 4 mm，电机转速为1 797 r/min。则转频fr为29.95 Hz，滚动轴承内圈点蚀故障特征频率为：

(9)

式中：α为接触角。采用加速度传感器采集振动信号，传感器固定在风扇端垂直方向，采用数据采集卡拾取振动信号，采样频率为12 kHz，本文计算时，采样点数取1 024，时间起点为t=0。

原始信号的时域波形见图5。

图5 原始信号的时域波形Fig.5 The original signal in time domain waveform

对原始信号应用标准LCD进行分解，得到5个ISC分量和1个剩余信号。将峭度值较高的前3个分量的时域波形列于图6。可以看出：①3个分量的时间尺度依次增大，体现了LCD的自适应滤波特性；②3个分量的幅值依次减小；③3个分量均能清晰体现调制现象；④3个分量的右端点处均存在较大的幅值发散现象，端点效应明显。前三点说明LCD的方法是有效的，第四点说明，LCD在短数据集的应用中仍存在端点效应，分解质量有待提高。

图6 标准LCD的ISC分量时域波形Fig.6 The ISCs of standard LCD in time domain waveform

对原始信号应用改进LCD进行分解，得到4个ISC分量和1个剩余信号。将峭度值较高的前3个分量的时域波形列于图7。可以看出：①3个分量的频段、幅值依次减小，且均能清晰体现调制现象；②改进LCD的分解结果中已没有明显的端点效应，而中间段的信号趋势与标准LCD结果一致。说明改进方法避免了端点邻域的插值“过冲”，且能保证了插值曲线与原始信号的相似程度。

图7 改进LCD的ISC分量时域波形Fig.7 The ISCs of improved LCD in time domain waveform

为了进一步对比分解质量，对标准LCD和改进LCD的分量分别求峭度，选取峭度最大的分量作为敏感分量，并分别求取Hilbert包络谱，将低频段(0～1 000 Hz)列于图8。

图8 敏感分量的包络谱Fig.8 Envelope spectrum of sensitive component

可以看出：①在图8(a)所示的标准LCD分解结果中，内圈故障特征频率fi处存在全谱幅值最大的谱线，内圈故障特征频率2倍频与转频的调制边频带(2fi-fr)处、内圈故障特征频率3倍频与转频的调制边频带(3fi-fr)处均存在较为清晰的谱线，说明LCD是有效的；②在图8(b)所示的改进LCD分解结果中，fi处的谱线仍为全谱幅值最大，但相对幅值较标准LCD分解结果突出，2fi-fr处、3fi-fr处的谱线更为清晰，说明了改进LCD的分解质量更高。

为了进一步分析效果，采用BLCD、镜像延拓LCD及SVM-LCD分别对原始信号进行分解。同样求取敏感分量的包络谱，并将结果列于图9。

可以看出：①BLCD与标准LCD的结果差别不大；②镜像延拓LCD分解结果中，fi、3fi-fr处的谱线比较清晰，但2fi-fr处谱线不如图8(b)清晰；③SVM-LCD的结果中，fi、2fi-fr处谱线与图8(b)的清晰程度相当，但3fi-fr处谱线更明显。因此，在试验数据的分析中，本文方法较BLCD及镜像延拓LCD的效果好，较SVM-LCD的效果稍差。考虑到算法复杂性及运行时间，本文方法具有一定的优势。

图9 各改进方法敏感分量的包络谱Fig.8 Envelope spectrum of sensitive component of each improved approach

混合插值方法能避免两个端点邻域的插值“过冲”现象，经过多项式方法对线性插值进行修正，整体上保证了插值曲线与原始信号的相似程度。因此，混合插值及多项式修正方法能有效抑制插值带来的端点效应，是一种有效的方法。

5 结 论

本文为了提高LCD的分解质量，给出了一种端点效应抑制方法：对端点邻域的数据段，采用线性方法插值均值曲线；将端点邻域的均值曲线与中间数据段合并，按照三次样条方法插值形成均值曲线；采用三次多项式修正端点邻域的均值曲线。混合插值方法避免两个端点邻域的插值“过冲”现象，多项式修正方法保证插值曲线的光滑程度。因此，混合插值及多项式修正方法能有效抑制插值带来的端点效应，是一种有效的方法。