摘要
当齿轮发生故障时,频谱中出现以啮合频率及其高阶谐波频率为中心、以齿轮旋转频率为间隔的多阶调制边频簇现象。为了自动聚焦故障边频成分,提出一种惩罚回归的故障边频簇提取方法,通过自适应稀疏群套索回归自数据驱动策略确定惩罚系数大小,在线更新频谱权重,以此找到故障边频簇。在稀疏群套索回归获得的各边频权重系数基础上,提出一种新稀疏群套索边带指标对齿轮传动系统进行健康监测,实现齿轮传动系统早期故障预警与定位。结果分析表明,所提出的方法可以实现更准确的齿轮早期故障预警与故障定位。
Abstract
Gear faults manifest in the frequency spectrum as multi-order modulation sideband clusters phenomenon center on the meshing frequency and its higher-order harmonics, spaced by the gear rotation frequency. In order to automatically focus the fault side frequency components, a method of fault side frequency cluster extraction with penalty regression was proposed. Adaptive sparse group lasso regression self-data-driven strategy was used to determine the penalty coefficient size and update the spectrum weight online to find the fault sideband clusters. Based on the sideband weight coefficients obtained from sparse group lasso regression, a new index called the sparse group lasso sidebands indicator was proposed for the health monitoring of gear transmission systems, enabling the early fault warning and location of gear transmission systems. Results analysis show that the proposed method can provide more accurate gear early fault detection and fault location.
在现代工业生产领域,齿轮传动系统的高效稳定运行对于确保工业设施的持续畅顺运转以及提升生产效率至关重要。然而,齿轮在恶劣工作环境下不可避免地会遭遇各类故障,例如磨损、裂纹和断齿等,这些故障不仅会导致生产中断和维修成本的增加,还可能引发严重的安全事故。因此,对齿轮传动系统进行健康监测尤为重要[1]。
在齿轮传动系统中,故障会导致振动信号的幅值和频率发生变化。每当齿轮旋转到故障位置时,都会产生冲击,齿轮故障部位只在特定位置参与啮合,因此冲击通常以齿轮旋转频率为周期产生,冲击频率低于齿轮啮合频率,频率较高的啮合频率为载波信号,由频率较低的冲击频率进行调制,形成新的调制信号[2-3]。频谱中出现以啮合频率及其高阶谐波频率为中心、以齿轮旋转频率为间隔的多阶调制边频簇现象。通过分析频谱上的边频成分,可以获得丰富的齿轮状态信息。另外,边频成分的幅值可以反映故障的程度[4]。针对边频特征提取,学者们提出了多种基于故障机理的健康指标:Yan等[5]通过融合频谱幅值构建复合健康指标来检测早期故障;Li等[6]提出自聚焦谱方法提取故障分量,并创建了WHI指数,用于对变速器进行在线监测。另一种普遍使用的是数据驱动的智能诊断方法,如Ha等[7]提出了基于领域知识的样本数据生成方法,通过从健康数据中生成故障数据,增加了故障样本数据。
齿轮传动系统发生不同故障时,频谱中会出现不同的边频成分,因此提取出故障齿轮的边频成分是进行齿轮健康监测的核心。传统方法中,基于经验知识的健康指标构建方法具有较强的物理可解释性,但是对复杂工况下不同单元之间的差异鲁棒性差;而数据驱动方法虽能适应不同工况条件,但需依赖大量带标签数据且可解释性差。为了解决上述两种方法存在的问题,引入稀疏群套索(sparse group lasso,SGL)回归进行故障边频成分选择。通过正则化和分组特征选择,SGL回归能够在保持模型简洁的同时,提高分类的准确性和组结构可解释性[8-9]。通过对边频成分成组应用SGL回归确定惩罚回归参数,在达到早期故障预警点后采用自适应稀疏群套索(adaptive sparse group lasso,ASGL)回归,使用原始特征成分构造一种稀疏群套索边带指标(sparse group lasso sidebands indicator,SSI),进行健康监测。所提出的自数据驱动方法融合经验知识和齿轮箱独有的动态特性,并基于故障机理原理,通过仿真数据和齿轮加速实验数据验证了方法的有效性。
1 SGL回归方法介绍
SGL方法在预先定义的变量组基础上引入l1、l2范数惩罚项选取相关变量[8]。该方法能够同时选择组级变量以及组内变量,为预设组双级选择方法,具体公式为:
(1)
其中:n为特征总个数;X(l)是X的子矩阵,其列对应于l组中的预测变量;m为组数;β(l)是该组的权重向量;pl是β(l)的长度;||·||2为l2范数,|·|1为l1范数;α、λ均是正则化参数。 α用来调整模型组内稀疏性,即控制非零组内非零系数的数量。λ用来调整组级稀疏性,即控制具有至少一个非零系数的组的数量。
2 ASGL回归参数确定及模型拟合算法
在进行SGL回归之前,需要确定α和λ两个参数。通常情况下,期望组内系数稀疏时选取α=0.95,不期望组内系数稀疏时选取α=0.05[9]。下面给出最优λ值的求解过程。
对于第k组,求使得该组全部系数为0的最小λ值:
(2)
其中:nk为组内特征数量;S(·)为坐标软阈值运算符,
(3)
r(-k)为Y与除第k组外所有组的拟合值的残差,。上述方程中恒为正,为了使得上述方程有解,令,需满足,因此计算k=1,2,···,K时的所有值得到一个序列,搜索该序列中全部非重复元素。非重复元素按降序排列,并加入元素0,最终序列表示为[λ1,λ2,···,λM],并形成M-1个区间[λm,λm+1],其中m=1,2,···,M-1且λM=0。
设置预期保留的组数为gexp。考虑到随着惩罚系数的增加,将保留更少的组数。因此按照λ区间的递减顺序,不断搜索最小λ值,使该组系数为0。达到组数gexp后,迭代步骤将停止。第m步的搜索过程按如下方式执行。
对于第k组,找到满足的索引pk,m,并计算满足该条件的元素数量nk,m,上述方程(2)化简为一元二次方程,可以求得对应参数如下:
(4)
(5)
(6)
(7)
如果δ≥0,求解,可以求得:
(8)
迭代所有组并找到满足上述条件的解。每个区间中的解被附加到集合Λ中。当集合Λ中的元素个数超过gexp时,停止迭代步骤。若此时Λ中非零元素个数为gexp,则输出λm+1作为最优λ值;若Λ中非零元素个数大于gexp,则输出此时集合Λ中所有元素的第gexp+1大的值作为最优λ值。
获得最优值后,进行模型拟合过程:
1)使用最小二乘估计初始化系数。
2)迭代所有组,计算组系数是否全为0。k组的计算过程如下:若,则,否则使用式(9)所示梯度下降法进行计算[10]。
(9)
其中,t为迭代步长,迭代终止条件为或迭代次数i≥Imax。
3)为了使组内产生稀疏性,对满足的特征置0,达到组内特征选择的效果。其中,j是组内该特征的序号,是Y的部分残差减去除X(k)j的所有其他协变量拟合。
使用以上方法,最终获得回归系数β。
3 齿轮箱故障诊断流程方法
基于齿轮故障时的边频簇特性,利用边频簇幅值构建指标SSI以反映齿轮系统的健康状况并进行故障定位。具体过程如下:
3.1 构建NSI
获得振动信号的频谱,并计算出各啮合对的啮合频率和各齿轮轴的转频。首先求边频成分的理论位置:
(10)
其中,fi,j,k表示当前时刻tk的频带中心理论位置,c是啮合频率,s是转轴转频,i是啮合频率的谐波阶次,j是齿轮所在轴转频的谐波阶次。
因啮合阶次频率成分一直较大,先将啮合阶次频率成分去除,再将同一啮合阶次、同一转轴转频下的边频簇求幅值,构造反映该齿轮健康状态的单齿特征指标SI(single indicator)。SI指标构建示意图如图1所示。
(11)
(12)
(13)
其中,SIj,k表示当前时刻tk第j个特征的SI指标值,y是频谱的幅值,τ是边频成分的频率波动范围,xτ是频谱的分辨率,Q是该齿轮所在轴转频阶次的最高谐波阶次, 是向下取整函数。

图1SI特征指标构建示意图
Fig.1Schematic diagram of the construction of SI characteristic indicator
由于不同齿轮的传动特性不同,最终构造的不同齿轮的SI指标在幅值数量级上也会存在显著差异。为了方便对不同齿轮健康状态演变过程的比较,将同一时刻构造的所有SI指标进行归一化,得到归一化的单个边频指标(normalized single-sideband indicator,NSI),以将其统一到同一数量级。tk时刻,NSI具体表达式为:
(14)
其中,μk为时刻tk的不同SI特征指标的均值,σk为时刻tk的不同SI特征指标的标准差。
3.2 组惩罚回归生成实时惩罚回归系数βj,k
步骤1:首先生成NSIH。当一个齿轮传动系统在健康状态时(即对应的标签为“0”),采集一段历史振动信号并计算出NSI值,记作NSIH。
(15)
(16)
其中,kH为健康状态的预设时间对应的数据点数,N为齿轮个数,Ni为第i个齿轮啮合频率的最高次谐波阶数。
步骤2:生成初始的NSIF。对一个齿轮传动系统的健康状态进行持续监测,采集另一段实时振动信号生成NSI,其健康状态未知,假设为“异常”状态(即对应的标签为“1”)。
(17)
其中,kW为异常状态所包含时间对应的数据点数。
步骤3:计算初始惩罚回归系数βj,k(表示当前时刻tk第j个特征的实时回归系数)。首先将NSIH和NSIF的全部特征进行分组,同一齿轮对应的特征为一组,同一齿轮啮合频率的不同倍频处的SI对应组内不同特征。使用前文提到的SGL回归模型拟合算法分离NSIH和NSIF获得初始的实时惩罚回归系数β。在tk时刻,将Y=(0,0,···,0,1,1,···,1)T(“0”的个数为kH,“1”的个数为kW)、、、、m=N分别代入式(1),最终可以求得实时回归系数βj,k=βj,1(矩阵β第j行第1列)
步骤4:更新NSIF。在下一监测时刻,采集一段新的实时振动信号生成新的NSI作为NSIF。
步骤5:计算实时惩罚回归系数βj,k。使用步骤3的SGL回归分离NSIH和新产生的NSIF获得实时惩罚回归系数βj,k。
3.3 生成新健康监测指标SSI
早期时刻λ对生成的SSI影响较小,为加快运算速度,将惩罚参数λ设置为定值并获得各个特征实时回归系数βj,k。由于对NSI进行标准化等原因,会出现少量的负系数。系数为正,说明该特征是反映齿轮故障的较重要特征;系数为负,说明此特征不是反映齿轮退化的明显特征,但负系数不代表与故障发生趋势相反,因此将负系数置零。
(18)
将当前时刻tk之前的所有标准化后的回归系数取平均数作为指标融合的权重系数wj,k:
(19)
最后获得当前时刻tk的SSI健康指标:
(20)
3.4 早期故障预警及故障齿轮定位
以前一段时间的SSI值为基准,计算μ+3σ作为早期报警阈值,利用连续触发预警机制[11]进行故障预警。故障发生后使用ASGL回归方法获得最优λ值,并使用前文模型拟合算法获得各特征权重系数,分析权重系数中非零的Q个特征对应的齿轮,计算各个齿轮对应的归一化概率pL:
(21)
其中,KL为权重系数非零的Q个特征中第L个齿轮对应的个数。根据概率大小进行故障定位,其中每个齿轮对应的归一化概率即为该齿轮发生故障的概率。
4 仿真验证
首先使用仿真信号来验证所提出的方法。齿轮箱结构和各齿轮特征频率分别如图2和表1所示。

图2齿轮箱结构示意图
Fig.2Schematic diagram of gearbox structure
表1模拟齿轮箱转频及啮合频率
Tab.1 Rotation frequency and meshing frequency of the simulated gearbox

振动仿真信号生成公式如下:
(22)
(23)
(24)
(25)
(26)
(27)
其中,ω(t)为幅值为1的高斯白噪声。仿真信号和加速寿命实验采用相同的参数值,采样频率为12 800 Hz。
仿真信号在0~150 s无边频成分;150~300 s出现以f2的2倍频为中心、fb为间隔的边频成分且该边频成分逐渐增大;300~400 s出现以f2的2、4倍频为中心,fb为间隔的边频成分且该边频成分逐渐增大;400~500 s出现以f2的2、4、5倍频为中心,fb为间隔的边频成分且该边频成分逐渐增大。
初始参数设置如下:gexp=3,α=0.05,λ=0.01,利用20 s的数据计算频谱,每2 s更新一次频谱,前30次频谱认为健康状态,故障窗长设置为10次频谱计算,迭代步长为1次频谱更新。回归前要进行特征分组,本仿真中,每一阶啮合频率周围计算的转轴转频最高阶次设置为3,每一阶啮合频率对应边频个数为2×3=6。由于采样定理的限制,频谱可以测量的极限为5 000 Hz左右(仿真的信号采样频率为12 800 Hz),因此将f1的最高谐波设置为4、f2的最高次谐波设置为10,所以4个齿轮的全部特征个数为2×4+2×10=28。分组方法为同一齿轮对应的特征为一组,共4组。
按照上述分组方式,画出齿轮Ⅲ啮合频率2倍频的边频成分的幅值变化,如图3所示。分析得出,仿真信号可以正确反映故障齿轮的退化信号。

图3齿轮Ⅲ啮合频率2倍频的幅值变化
Fig.3Amplitude variation of 2 times the engagement frequency of gear Ⅲ
SSI显示的初始退化时刻为162 s,如图4所示,与实际情况基本一致。该指标可以用于早期故障预警。

图4SSI实时监测仿真齿轮箱
Fig.4SSI monitoring on the simulated gearbox
使用SGL回归获得齿轮Ⅲ啮合频率2倍频特征指标回归系数随时间的变化,如图5所示。回归系数变为非零对应时刻为162 s,与故障发生时刻150 s基本一致。对齿轮Ⅲ的啮合频率4倍频和5倍频进行特征分析,回归系数变为零的时间和真实故障时间基本一致。

图5齿轮Ⅲ啮合频率2倍频的实时回归系数
Fig.5Gear real-time regression coefficient of 2 times the engagement frequency of gear Ⅲ
除上述外的其他边频对应回归系数始终为0。因此齿轮Ⅲ的故障概率为1,其余齿轮故障概率为0,可以进行故障定位。
5 齿轮加速寿命实验验证
为了验证提出方法对齿轮健康监测的效果,使用齿轮加速寿命实验台的数据进行验证。采集信号为测试齿轮箱振动加速度信号(共6通道,输入轴X、Y、Z三个方向,输出轴X、Y、Z三个方向);采样频率为12.8 kHz、采样时长2.56 s、采样间隔2 min;输入转频为40 Hz,齿轮箱结构示意图如图2所示。各啮合频率及各转轴转频与仿真案例相同。设备全寿命为110 h,实际的失效形式为齿轮Ⅲ的四齿连续断齿,齿轮Ⅲ与齿轮Ⅳ各齿面均有点蚀。因故障齿轮距离输入轴较远,故障信号微弱,为验证本方法对于微弱故障特征的提取效果,使用输入轴Y方向(此方向为加载方向)的振动加速度信号。
本实验中该齿轮传动系统的全寿命振动时域波形如图6所示。通过时域振动信号无法看出早期故障发生时刻,只能发现在设备寿命末期振动加速度幅值显著增大,难以进行健康状态监测。

图6振动信号时域波形图
Fig.6Vibration signal time domain waveform
使用本方法对齿轮传动系统进行健康监测。将前30个振动信号文件视为正常状态下的信号。使用窗长为30且增量为1的移动窗口更新收集到的信号作为“异常信号”。
本实验中,初始参数设置如下:gexp=2,α=0.05,λ=0.01。将所有特征先进行标准化,对正常状态(取“0”)的NSIH和异常状态(取“1”)的NSIF进行SGL回归获得回归系数;对每一个特征,取当前时刻之前所有回归系数的均值作为该特征的权重系数,得到最终的SSI,并实时更新。前期健康监测阶段,因为自适应更新的λ值的变化对特征系数和指标值影响不明显,不会对早期故障预警准确度产生显著影响。当达到早期故障预警点后,这时λ的变化会对特征系数的值产生较大影响。为了准确进行故障定位,这个时期对λ的值采用ASGL算法进行实时更新,使用实时更新的λ进行特征系数的计算,根据每个齿轮对应特征的非零系数的个数进行故障概率求解,定位故障齿轮。整个阶段求出的新指标SSI特征值如图7所示。

图7SSI实时监测实验齿轮箱
Fig.7SSI real-time monitoring on experimental gearbox
为了对设备的早期故障时间进行判断,下一步进行报警阈值的设定。本实验中,选取SSI的前500个数据点进行计算(如图7所示),得到μ=0.402 3,σ=0.089 6,早期故障报警阈值μ+3σ=0.671 1,报警时刻为5 542 min,说明5 542 min时刻开始发生了早期故障。
为验证所提方法在早期故障预警和定位方面具有更优的性能,与几种方法指标进行了对比。图8~11分别给出了OWIE [12](optimal weight impulse extraction)、GITS [13](Gini index of trend spectrum)、峭度、有效值与负熵指标的时域图。有效值、峭度和负熵几种常见指标只有在临近失效时才出现了明显的上升趋势,难以进行设备的健康监测。在故障早期求出各指标指示的早期故障预警时刻,在早期故障发生后对各指标幅值的趋势性和鲁棒性进行计算,结果如表2所示。通过指标对比,所提出的指标SSI具有更准确地早期故障预警时间,优于其他的指标,明显优于有效值、峭度和负熵三种常用指标。SSI与对比指标鲁棒性相差不多,但在趋势性上有明显优势,更能反映故障的发展情况。

图8OWIE指标监测实验齿轮箱
Fig.8OWIE monitoring on experimental gearbox

图9GITS指标监测实验齿轮箱
Fig.9GITS monitoring on experimental gearbox

图10峭度监测实验齿轮箱
Fig.10Kurtosis monitoring on experimental gearbox

图11有效值与负熵监测实验齿轮箱
Fig.11Root mean square value and negative entropy monitoring on experimental gearbox
表2SSI与其他指标对比
Tab.2 SSI compared with other indicators

如表3所示,SGL回归方法计算得出齿轮Ⅲ发生故障概率为66.67%,齿轮Ⅱ发生故障概率为16.67%,齿轮Ⅰ和齿轮Ⅳ发生故障概率为8.33%。ASGL回归方法计算得出齿轮Ⅲ发生故障概率为80%,齿轮Ⅳ发生故障概率为20%,实际拆机检查结果为齿轮Ⅲ的四齿连续断齿,齿轮Ⅲ与齿轮Ⅳ各齿面均有点蚀,与实际情况一致。相比SGL回归方法,该方法可以正确选择出特征组,更为准确地进行齿轮故障定位,验证了该方法的有效性。
表3故障齿轮定位
Tab.3 Faulty gear location

6 结论
针对齿轮传动系统的健康监测,本文提出了使用ASGL回归方法确定惩罚参数值,并提取出故障边频,提出了SSI健康指标,在故障初期阶段进行故障预警和故障定位。通过仿真数据及加速寿命实验数据对本方法进行了验证,SSI在齿轮健康监测方面表现出色,与其他常用的监测指标相比,趋势性有明显优势,更能反映故障的发展情况,具有更准确的早期故障预警性能。此外,该方法也可以提取出故障边频簇并计算出各齿轮的故障概率从而进行故障定位。