滚动轴承作为旋转机械中不可缺少的零部件，同时也是最易发生故障的零部件之一，其工作状态将直接影响机械设备的运行．研究滚动轴承的故障机理，对其进行故障诊断具有十分重要的意义[1]．目前，基于振动信号分析轴承故障是一种比较完善的诊断方式．然而由于噪声的干扰，直接对原始信号进行分析往往无法获得准确的诊断结果．针对这一问题，对信号实施降噪、特征提取等预处理显得尤为重要[2-5]．由于信号表达的灵活性，稀疏分解在特征重构中独具优势，在图像融合、人脸识别等领域有着广泛的应用．近年来，稀疏分解在故障诊断领域也得到快速发展．文献[6]提出一种基于循环结构学习字典的稀疏分解算法，实现了周期性故障特征的提取；文献[7]结合时频流形分析与稀疏分解实现了滚动轴承故障特征的增强，提升了诊断效果；文献[8]结合非局部平均与自相关脉冲谐波噪声比(AIHN)对K-SVD算法学习到的原子进行增强与选取，然后根据增强后的原子构造字典，基于该字典实现滚动轴承故障特征的提取．基于稀疏分解的局部性故障诊断是以重构信号中的冲击特征实现的[9]．重构精度与字典中的原子有直接关系，同时字典的大小也影响计算效率，这使得稀疏重构的效果与效率存在矛盾．K-SVD算法在一定程度上能够自适应学习故障特征以生成字典，但其准确性基于较为干净的故障信号，而在工程实际中获得这样的信号较为困难．此外，传统的基于残差信号能量阈值的迭代终止准则在实际应用中存在明显的不足．由于信号先验信息的不足，阈值选取须多次尝试或根据普适的经验设定，这导致处理不同信号时缺乏适应性，重构结果也具有随机性[10]．针对上述问题，本研究提出一种双树复小波分解与遗传算法相结合的改进稀疏分解方法，平衡了效率与精度需求，实现了强噪声信号的深度降噪；同时，提出以残差信号包络熵来自适应地选取匹配追踪的迭代次数，克服了稀疏分解应用于故障诊断领域时迭代终止条件设置困难的问题．1 改进的稀疏分解1.1　双树复小波传统的离散小波变换存在小波系数振荡明显、平移不变性差及易发生模态混叠等不足，为解决这些问题，提出了双树复小波变换[11]．定义复小波为ψc(t)=ψh(t)+iψg(t), (1)式中：ψh(t)为实部树小波，ψg(t)为虚部树小波，两者均为实小波；i为虚数单位．双树复小波实部的小波系数dIjRe及尺度系数cIjRe的计算如下dIjRe(k)=2j/2∫-∞∞x(t)ψh(2jt-k)dt; (2)cIjRe(k)=2J/2∫-∞∞x(t)φh(2Jt-k)dt, (3)式中：j为尺度因子，j=1，2，…，J，J为分解层数；k为小波滤波器长度；ψh(•)为实部树小波函数；φh(•)为实部树尺度函数．同理，虚部的小波系数dIjIm和尺度系数cIjIm表示为dIjIm(k)=2j/2∫-∞∞x(t)ψg(2jt-k)dt; (4)cIjIm(k)=2J/2∫-∞∞x(t)φg(2Jt-k)dt, (5)式中：ψg(∙)为虚部树小波函数；φg(∙)为虚部树尺度函数．双树复小波变换输出的小波系数与尺度系数表示为：djc(k)=dIjRe(k)+idIjIm(k); (6)cJc(k)=cIjRe(k)+icIjIm(k). (7)对某一尺度的小波系数保留，其他尺度系数置零以实现特定尺度下小波系数的单支重构，表达式为：dj(t)=2(j-1)/2∑kdjRe(k)ψh(2jt-k)+i∑kdjIm(k)ψg(2jt-k); (8)cJ(t)=2(J-1)/2∑kcJRe(k)φh(2Jt-k)+i∑kcJIm(k)φg(2Jt-k). (9)1.2　基于遗传算法的稀疏重构信号的稀疏分解是从过完备字典中依次匹配与信号最为相似的原子，以实现信号最稀疏表示的过程．Gabor字典是最常用的时频字典，由Gabor函数参数的离散化生成．Gabor函数表示如下gγ(t)=1sgt-usexp[i(ξt+ϕ)], (10)式中：g(t)=exp(-πt2)为高斯窗函数；γ=(s,u，ξ,ϕ)为字典的参数空间，s,u,ξ,ϕ分别为尺度因子、位移因子、频率因子和相位．参数空间的离散化形式为γ=(aj，pajΔu，wa-jΔξ，iΔϕ)，其中a=2，Δu=1/2，Δξ=π，Δϕ=π/6，0j≤log2 L，0≤p≤2-j+1L，0≤w≤2j+1，0≤i≤2，L为信号长度，字典中原子个数为52(Llog2 L+L-1)．基于上述Gabor字典，本研究采用结合遗传算法优化的匹配追踪对故障特征信号进行重构．遗传算法通过选择、交叉等运算模拟生物进化中基因的遗传行为，以适应度函数评价种群中个体的优劣程度，根据概率计算选择部分优秀个体作为新的种群，以此不断进化，最终使种群中包含最优个体[12]．本研究遗传算法中相关设置如下．个体：由生成Gabor原子的参数所组成的集合．种群：在各因子的设置范围内，随机生成多个参数集合，构成种群，其中各参数范围与前述相同．种群大小通常设置为20～100，本研究种群个数设置为50．目标函数：即自适应度函数，为残差信号与个体参数所构造的Gabor原子的内积．选择：目标函数值高的个体遗传至下一代，以保证全局收敛性，其中选择概率设为0.5．交叉：对选择后的个体进行交叉操作，生成同等数量的新个体．终止准则：采用进化代数控制程序的终止，其范围一般为100～500，本研究设置迭代次数为500．1.3　基于残差信号包络熵的终止准则包络谱能有效反映信号的冲击特性．熵是信号不确定性与复杂程度的度量，系统越复杂、越不确定其熵值越大[13]．基于上述理论，提出基于残差信号的包络熵的迭代终止准则．对原始信号进行双树复小波分解后，结合峭度最大准则，选取最优频带，并对该频带计算包络谱．设x(t)为一般信号的时域表示，通过柯西主值对信号进行希尔伯特变换，表示为x̂(t)=1π∫-∞∞x(t)t-τdτ. (11)信号的包络表示为e(t)=x2(t)+x̂2(t). (12)然后，计算信号包络的功率谱即为信号x(t)的包络谱P(t)．将包络谱划分为M块，以每块包络谱占总包络谱能量的比例dm作为熵值计算中事件发生的概率，划分块数越多，计算越准确，但同时考虑到计算效率问题，本研究设定划分块数M=128．因此，信号x(t)的包络熵H可表示为：H=-∑m=1Mdmlg dm;dm=Pm/∑m=1MPm. (13)2 仿真验证为验证DTCWT在信号处理中的有效性，仿真构建一组由冲击信号h1(t)、正弦信号h2(t)及噪声h3(t)构成的混合信号h(t)，混合信号的函数表达式如下：      h1(t)=∑m{exp[(-ζ1/1-ζ12)2πf1(t-τ1-mT0)]sin[2πf1(t-τ1-mT0)]};                          h2(t)=sin(2πf2t);                          h3(t)=rand(0,1);                   h(t)=h1(t)+h2(t)+h3(t). (14)信号h1(t)的仿真参数为：衰减阻尼ζ1=0.06，振荡频率f1=600 Hz，初始时间τ1=0.01 s，冲击周期T0=0.1 s，由冲击周期可计算信号的特征频率f0=10 Hz．信号的h2(t)的仿真参数为正弦频率f2=20 Hz．信号h3(t)为满足方差为1的正态函数．采样频率为10 kHz，仿真时长t=1 s，图1(a)为h(t)的时域波形，其中A为h(t)的时域振幅．10.13245/j.hust.210611.F001图1仿真验证结果对h(t)进行4层双树复小波分解，由分解获得各分量的时域波形如图1(b)所示，其中：D1，D2，D3和D4为高频分量；A4为低频分量．根据各分量的时域波形可知，双树复小波成功地将混合信号h(t)分离，其中：D1，D2和D3主要为噪声信号，D4主要为冲击信号；A4主要为正弦信号．为进一步验证本文提出方法在处理强噪声信号时的优越性，对混合信号h(t)添加-10 dB的高斯白噪声构造高噪声信号，其时域波形如图2所示．由混合信号中冲击的个数设置迭代次数n=10，图3为未经双树复小波处理的信号的稀疏重构结果，图4为处理后信号的稀疏重构结果，其中：A为信号的时域振幅；P为信号频域包络的能量幅值；f为频率．经对比，预处理后的重构信号其时域冲击特征被准确定位，频域倍频现象更加明显．10.13245/j.hust.210611.F002图2高噪声信号10.13245/j.hust.210611.F003图3未经处理信号的稀疏重构结果10.13245/j.hust.210611.F004图4处理后信号的稀疏重构结果3 实验验证3.1　实验设置为了验证本文方法在处理实际信号时的有效性，分别采用不同故障的滚动轴承实验数据进行验证．实验中采用N205EM型圆柱滚子轴承，其主要参数为：节径Dm=38.5 mm，滚子直径Db=7.0 mm，接触角α=0°，滚子数Z=13．其中故障尺寸为0.25 mm×0.30 mm(宽×深)．所用加速度传感器型号为PCBMA352A60，其输出带宽为5 Hz～60 kHz，灵敏度为10 mV/g．实验转速N=1 000 r/min，采样频率fs=50 kHz．根据轴承的故障机理[14]可知，轴承内圈、外圈和滚动体故障特征频率的计算如下：fi=12Z1+DbDmcos αfr; (15)fo=12Z1-DbDmcos αfr; (16)fb=12DmDb1-Db2Dm2cos2 αfr, (17)式中fr为轴承所在轴的旋转频率．经计算，内圈故障特征频率fi=128.08 Hz，外圈故障特征频率fo=88.63 Hz，滚动体故障特征频率fb=44.33 Hz．3.2　外圈故障图5(a)是正常轴承的振动信号，图5(b)是外圈故障轴承振动信号，其中采集信号时间长度t=0.6 s．诊断流程主要分为预处理、稀疏重构以及频域分析．首先对故障振动信号进行双树复小波分解，根据峭度最大准则选取含故障冲击成分的最优分量作为进一步处理的对象，图5(c)为最优分量的时域波形．10.13245/j.hust.210611.F005图5轴承振动信号对选取的最优分量进行稀疏重构，图6为迭代过程中各特征量的变化曲线，图中：Ec为残差信号的包络熵；Ey为残差信号能量．残差信号的包络熵值随着迭代次数的增加，呈现出先增加后保持平稳的趋势，这是由于信号中能量集中的周期性冲击被重构后，残差信号逐步变为无序的噪声信号，因而其包络熵值呈现先增大后平稳的规律，选取转折处的迭代次数(n=287)为迭代的终止条件，可获得最优的重构效果．10.13245/j.hust.210611.F006图6迭代过程中各特征量的变化曲线1—残差信号包络熵；2—残差信号与重构信号包络熵差；3—残差信号能量；4—重构信号包络熵．而残差信号能量则随着迭代次数的增加一直呈现减小的趋势，传统的以残差信号能量作为终止条件设定的准则，由于未知信号的先验知识，使得重构结果存在很大的随机性．同时，遗传算法使得寻优过程中内积计算的效率获得提升．表1对比了传统的匹配追踪(MP)与基于遗传算法改进的匹配追踪(GAMP)字典大小及迭代时间．为方便计算假设信号长度L=1 000．10.13245/j.hust.210611.T001表1计算效率对比参数GAMPMP字典类型实时字典固定字典字典大小50L570 388L寻优一次用时/s0.47171.12迭代的总用时/s2.5n171.12n图7反映了外圈故障的重构信号的时域频域特征，重构信号的时域波形中存在明显的冲击，包络谱中出现外圈故障特征频率及其倍频．为说明本文方法的优越性，与直接对信号进行稀疏重构的结果相对比，结果如图8所示．10.13245/j.hust.210611.F007图7外圈故障诊断结果1—重构信号时域波形；2—原信号时域波形．10.13245/j.hust.210611.F008图8对比实验结果1—重构信号时域波形；2—原信号时域波形．与图7结果相比，图8中重构信号的包络谱图中除了故障频率及其倍频外还存在能量较大的干扰频率，这将妨碍故障信息的正确识别．3.3　内圈及滚动体故障对内圈及滚动体故障信号进行与上节中相同的处理，以验证本文方法在识别内圈和滚动体故障过程中的有效性．图9为内圈故障与滚动体故障轴承的重构信号包络图，包络图中出现明显的倍频现象，其基频分别与轴承内圈故障特征频率fi和轴承滚动体故障特征频率fb相等．10.13245/j.hust.210611.F009图9内圈和滚动体故障诊断结果4 结语稀疏分解在信号的处理上具有灵活、高效的优势，但在实际应用中存在重构不准确、终止条件难以选择等问题．本研究针对传统的稀疏分解故障诊断方法存在的不足，提出了结合双树复小波与遗传算法改进的稀疏分解算法，实现了强噪声信号的深度降噪，提升了计算效率，在处理高噪声、高维度的实际工程振动信号的过程中，更具自适应性与可行性．提出了基于残差信号包络熵的终止准则，解决了传统迭代终止条件误差大、随机性强的问题．这一准则可以自适应选取迭代次数并且无需已知故障信号的先验信息，避免了以经验和试验确定终止条件的不足．最后，通过仿真及实验验证，表明本文提出的方法与传统稀疏分解相比，对冲击特征的检测更为准确，重构信号中包含噪声更少，包络谱中的倍频现象更为明显，能够准确地实现轴承内圈、外圈及滚动体的故障诊断．