文章编号：1004-2539（2021）04-0029-08DOI：10.16578/j.issn.1004.2539.2021.04.005基于指数窗截取递推最小二乘法的齿轮啮合刚度辨识算法研究
王茂辉李海翔杨平陈娇夏伟
（重庆工商职业学院智能制造与汽车学院，重庆401520）
摘要齿轮在机械传动系统中有着广泛应用，由于齿轮啮合过程中参与啮合的轮齿对数周期变化，因此，齿轮啮合刚度为时变参数，在啮合时会产生啮合振动。当齿轮副出现齿根裂纹时，啮合刚度会减小，齿轮啮合产生的系统振动响应也发生改变，通过振动响应辨识齿轮啮合刚度能够监测齿轮副的健康状态。针对齿轮啮合刚度的时变特征，提出了基于指数窗截取递推最小二乘（Exponential window recursive least square，EWRLS）算法和振动信号瞬时频率的齿轮啮合刚度辨识方法。进行啮合刚度辨识时，EWRLS算法将输入、输出齿轮的转速曲线分别作为辨识输入信号和观测信号，使用指数窗函数进行数据截断，使用递推最小二乘算法估计系统参数。为了计算输入、输出齿轮的转速曲线，使用经验模态分解（Empirical mode decomposition，EMD）方法将振动信号分解为具有不同变化频率的本征模态函数（Intrinsic mode function，IMF），并根据IMF的平均频率重构输入、输出齿轮的特征信号。通过Hilbert变换计算特征信号的瞬时频率曲线，从而获得各齿轮的转速曲线。使用仿真和实测信号对算法进行验证，结果表明，EWRLS算法能够辨识齿轮副的时变啮合刚度。
关键词刚度识别齿轮啮合刚度递推最小二乘法指数窗瞬时频率
Rearch of Gear Meshing Stiffness Identification Algorithm bad on Exponential
Window Interception Recursive Least Square Method
Wang Maohui Li Haixiang Yang Ping Chen Jiao Xia Wei
（School of Intelligent Manufacturing and Automotive，Chongqing Technology and Business Institute，Chongqing401520，China）
Abstract Gear is widely applied in mechanical transmission systems.Gear meshing stiffness is a time varying parameter becau of the periodical changing of the numbers of teeth involved in meshing process，the meshing vibration is generated during gear meshing.When the cracks are existed in roots of gear teeth，the mesh⁃ing stiffness decread，vibration respon of the system due to gear meshing changed as well.Thus，identifying gear meshing stiffness through vibration respon is a method monitoring healthy state of gear pairs.For the time-varying character of meshing stiffness，a gear meshing stiffness identification algorithm is propod bad on exponential window recursive least square method（EWRLS）interception and instantaneous frequency of the vibration signal.Dur
ing the identification of meshing stiffness，EWRLS algorithm takes the speed curves of the input and output gears as the identification input signal and obrvation signal，respectively.In the algorithm，the exponential window is applied to intercept data，and the recursive least square algorithm is applied to esti⁃mate parameters of the system.To calculate speed curve of the input and output gears，the empirical mode de⁃composition（EMD）method is ud to decompo vibration signal into intrinsic mode function（IMF）of differ⁃ent frequency，the IMFs are ud to reconstruct character signals of input and output gears basing on the mean frequency of IMFs.The Hilbert transform is applied to calculate the instantaneous frequency curve of the charac⁃ter signals to obtain the speed curves of gears.The simulated signal and measured signal are ud to validate the algorithm，results show that the EWRLS algorithm can identify time-varying meshing stiffness of gear pair.
Key words Stiffness identification Gear Meshing stiffness Recursive least square method Expo⁃nential window Instantaneous frequency
0引言
齿轮副由于具有能够传递大转矩、结构紧凑、传动效率高等优点，而被广泛应用于机械传动系统之中。
齿轮传动过程中会由于单、双齿交替啮合导致啮合刚度随时间周期变化，从而产生刚度激励。齿轮刚度往往与齿轮副的健康状态相关，轮齿裂纹、齿面磨损、弯曲疲劳、断齿等故障发生时，啮合刚度也会相应减小，减小的程度与故障程度有关，因此，啮合刚度也反映了齿轮故障的程度[1]。齿轮啮合刚度的变化会导致齿轮系统振动响应的变化，通过系统的振动响应可以监测系统状态，从而保证系统安全工作。
一般基于振动信号对系统进行状态监测的方法包括时域和频率两种方法。时域方法的特征指标包括峰-峰值、有效值、平均功率以及峭度等，频域分析方法则将振动信号傅里叶频谱在特定频率处的幅值作为特征指标[2]。时域、频域指标能够很好地反映被测齿轮传动系统的振动状态，但是其幅值往往与齿轮转速、负载以及支承结构参数有关；系统参数改变后，指标数值也随之改变[3]。而刚度参数作为齿轮的固有特征，其指标数值大小仅与齿轮的啮合状态有关，与系统的其他参数无关[4]，因此，齿轮传动系统的刚度指标具有随外界参数变化波动较小、稳定性好等优点，通过振动信号对齿轮啮合刚度进行直接辨识，可以得到齿轮状态的稳定表示。
目前，不同方法所识别的刚度以定刚度为主，识别方法以系统参数识别方法和模态实验为主。系统参数识别方法包括最小二乘算法[5]、扩展卡尔曼滤波算法[6]、遗传算法[7]、自适应辨识算法等[8]。模态实验方法一般先利用实验获得系统参数，再通过实验结果与被测系统特征参数的差异来确定啮合刚度的数值。董冠华等[9]基于模态分析理论对结合部动刚度辨识方法进行研究，建立了动力学特性和系统
固有频率的关系，并通过模态分析得到的固有频率对结合部刚度进行了辨识。王二化等[10]41-44使用粒子群优化算法，应用锤击实验测试结果辨识了轴承、螺母和导轨的径向刚度。Wang等[11]1-24提出了一种基于刚度和阻尼边际曲线的结构动力学非线性模型及参数辨识方法，该方法定义了多个非线性指数用以描述非线性模型的特征，结合支持向量机方法和非线性最小二乘算法对非线性模型进行辨识。
由于转轴的瞬时旋转速度在很大程度上反映了传动系统的工作状态，因此，一些研究人员开展了使用振动信号估计转轴瞬时旋转频率的研究。康德等[12]利用变分模态分解（Variational mode decomposi⁃tion，VMD）变换对振动信号进行降噪重构，然后对重构信号频谱进行Viterbi瞬时频率估计，该方法提高了瞬频的估计精度。陈建新等[13]结合Hilbert变换和傅里叶变换微分性质提出了一种针对声信号的瞬时频率特征提取方法。程为东等[14]采用经验模态分解（Empirical mode decomposition，EMD）降噪方法优化了基于瞬时故障特征频率（Instantaneous fault char⁃acteristic frequency，IFCF）[15]的滚动轴承瞬时转频估计方法。
模态识别方法或其他辨识方法能够准确地计算系统刚度，但所辨识的刚度多为定刚度，且一般都需要先进行实验标定，方法的参数跟踪能力和经济性较差。考虑到指数窗函数能够对数据进行指数加权，减小截断误差，提高识别方法的跟踪能力[16]，因此，本研究对既有递推最小二乘法进行了改进，提出了指数窗递推最小二乘（Exponential window re⁃cursive least square，EWRLS）算法。所提出的
算法使用EMD方法进行信号去噪，利用本征模态函数（Intrinsic mode function，IMF）平均频率提取特征信号，使用Hilbert变换计算一级传动齿轮系统输入、输出齿轮的旋转频率。算法将输入、输出齿轮的旋转角速度作为输入参数，在系统工作状态下完成齿轮的啮合刚度辨识，解决了啮合刚度的跟踪辨识问题，且测试系统中仅使用了单个振动传感器，解决了以往系统参数测试过程中使用设备较多，费用昂贵的问题。
1齿轮传动动力学离散模型
假设主动齿轮和从动齿轮仅存在旋转自由度，典型的齿轮传动动力学模型如图1
所示。
图1齿轮啮合动力学模型
Fig.1Dynamics model of gear meshing
如图1所示，主动齿轮的旋转角度为φ，从动齿轮的旋转角度为θ，主动齿轮节圆半径为r1，从动齿轮节圆半径为r2，齿轮副之间的啮合刚度为K（t）。由于齿轮副在不同啮合角度下参与啮合的轮齿对数不同，因此，啮合刚度为一时变参数。由于齿轮副的啮合阻尼比较小，阻尼比的取值范围一般在0.03~ 0.17范围之内[17]，因此，模型中忽略了啮合阻尼比。图1所示模型的动力学系统微分方程为
éëêùûúJ 100
J 2éëêùûúφθ+éëêùûúK (t )r 21
-K (t )r 1r 2-K (t )r 1r 2K (t )r 22éëêùûúφθ=éëêùû
ú00（1）式中，J 1为主动齿轮的转动惯量；J 2为从动齿轮的转
动惯量。假设初始状态下齿轮静止，令θ和φ
分别为ωΘ和ωΨ，则转角之间的函数关系可写为
J 2ωθ+K (t )tr 22∫ωθd t -K (t )tr 1r 2∫
ωφd t =0
（2）
通过拉氏变换，可以推导系统的传递函数为
ωθ(s )ωφ(s )=
1J 2K (t )r 1r 2
s 2+K (t )J 2r 22
（3）
应用冲击响应不变法，可以将系统传递函数由拉式变换域转换到z 变换域，系统传递函数的两个极点分别为
s 1，2=±
2
（4）
式中，j 表示虚数。
应用部分分式展开法可以将传递函数转化为
ωθ(s )ωφ(s )=12j J 2K (t )r 1s +
r 2+-12j J 2K (t )r 1s -
2
（5）转换后，系统的z 域传递函数为ωθ(z )
ωφ(z )
=
1-2cos(
2T )z -1+z -2
（6）
式中，T 为采样周期。利用z 变换的性质，可以将连续域的传递函数转换为离散差分方程形式，以便于在计算机中通过离散点的数值进行参数辨识。系统的离散模型为
ωθ(n )=p 1ωθ(n -1)-ωθ(n -2)+p 2ωφ(
n -1)（7）其中，
p 1=2cos(2T )（8）p 2
=
K (t )r 1sin (
2T )（9）
通过辨识参数p 1或p 2，就能够获取齿轮副的啮合
动刚度。由于通过参数p 2计算刚度K （t ）需要将正弦函数作级数展开并略去其中的高阶项，计算过程较
复杂且存在一定误差，因此，通过参数p 1的辨识结果p 1计算啮合动刚度K （t ）是较适宜的。啮合动刚度的计算式为
K (t )=J 2(arccos (p 1/2)r 2T
)
2
（10）
不同类型故障对啮合刚度的影响不同，对啮合刚度影响最大的故障类型为齿根裂纹，齿面点蚀、齿形误差等其他故障类型对刚度的影响在2%以内[18]。根据齿轮时变啮合刚度的计算理论[19]，对于一般的斜齿轮，当齿根裂纹长度分别为10mm 、20mm 和30mm ，裂纹深度为3mm 时，含裂纹的轮齿通过啮合区间时齿轮副啮合刚度如图2所示。
当齿根裂纹长度为30mm ，裂纹深度分
别为
1mm 、2mm 和3mm 时，含裂纹的轮齿通过啮合区间时齿轮副啮合刚度如图3所示。
图2
不同裂纹长度下啮合刚度
Fig.2
Meshing stiffness under crack with different length
图3
不同裂纹深度下啮合刚度
Fig.3
Meshing stiffness under
crack with different depth
由图2和图3可以看出，齿根裂纹对啮合刚度的影响较大，不同裂纹长度和深度造成刚度减小的程度不同。不同裂纹长度和深度下齿轮啮合刚度最大减小量如图4所示。
图4
不同裂纹长度和深度下刚度最大减小百分比
Fig.4Maximum percentage reduction of stiffness under different
crack length and depth
图4所示的结果表明，裂纹深度的增加更能够导
致刚度下降，裂纹深度在1mm 以内且裂纹长度在10mm 以下时，最大啮合刚度的下降率在1%以下，当裂纹深度达到3mm 且裂纹长度达到30mm 时，最
大啮合刚度的下降率约为10%。轮齿含有不同深度或不同长度裂纹时，最大刚度的下降量存在差异。因此，通过比较齿轮系统正常和当前状态下的最大啮合刚度差异，就能够估计轮齿是否存在裂纹，同时可以估计齿根裂纹的深度和长度。
2啮合刚度识别算法
2.1
RLS 算法
一般的因果系统可离散为
z (k )=-∑i =1
n
a i y (k -i )+∑j =0
m
b j u (k -j )+v (k )（11）
式中，z （k ）为系统输出量的第k 次观测结果；y （k ）为系统输出的第k 次真值；u （k ）为系统的第k 次输入
值；ν（k ）为均值为0的系统输入随机噪声；m 和n 为系统阶数。定义辨识过程中的第k 次辨识向量为
h (k )=[-y ()k -i u (k -j )]
（12）
式中，i =1，2，…，m ；j=1，2，…，n 。定义待辨识的参数向量为
ξ=[a i b j ]
（13）
当使用N 次观测结果进行最小二乘法参数辨识
时，辨识结果为
ξ=(H T N H N )-1H T
N Z N
桑叶祛斑（14）
式中，辨识矩阵H 和观测向量Z 的定义分别为
H N =éëêêêêêùûúúúúúh (1)h (2)⋮h (N )，Z N =éëêêêêêù
û
úúúú
úz (1)z (2)⋮z (N )（15）
最小二乘法（Least square ，LS ）通过最小化误差的平方和寻找数据的最佳拟合参数，对于多输入
单输出系统，LS 法能计算基于数据的全局最优无偏估计。但是，由于LS 法的计算结果无法跟踪系统参数的变化，因此，无法对具有时变参数的系统进行参数估计，即使使用加权最小二乘（Weighted least square ，WLS ）法，也无法解决辨识过程中因数据饱和、新观测值对参数估计起不到修正作用，导致估
计误差增大的现象。对具有时变参数的系统进行参数估计时，需要使用递推最小二乘（Recursive least square ，RLS ）算法，RLS 算法利用新的观测结果不断修正原有的辨识结果，达到识别系统动态参数的目的。当新的观测结果进入到辨识过程中后，辨识矩阵和观测向量分别为
H N +1=éëêù
ûúH N h (N +1)，Z N +1=éëêùû
úZ N z (N +1)（16）
假设辨识向量中各数值对辨识结果的影响是平均的，即加权矩阵为单位阵，则原观测值的递推辨识矩阵和新观测值加入后的递推辨识矩阵分别为
P N =[H T N IH N ]-1，P N +1=[H T N +1IH N +1]
-1
（17）RLS 算法的辨识计算过程为P N +1=P N -A N +1h (N +1)P N托班
（18）
A N +1=P N h T (N +1)⋅
[w -1(N +1)+h (N +1)P N h T (N +1)]-1
（19）
ξN +1=ξN +A N +1[z (N +1)-h (N +1)ξN ]（20）
式中，A 为辨识增益矩阵，表示对原参数的修正程度，A 越大则辨识结果越能够跟随参数的变化；w 为遗忘因子，表示旧数据对新辨识结果的作用，当
w （N ）为1时，递推最小二乘算法变为一般最小二乘算法，w （N ）小于1时，则表示削弱旧数据的作用。2.2
EWRLS 算法
由于RLS 算法的递推辨识矩阵同样会由于新观测值的进入而不断扩大阶数，最终导致新数据被旧数据
淹没，因此，需要修改新、旧数据的权重，同时对辨识向量进行截断，从而达到辨识时变参数的目的。既有方法一般通过矩形窗对数据进行截断或使用指数窗对旧数据进行截尾。两种方法各有优势，矩形窗可以减小辨识矩阵的阶数，特别适合系统在线辨识；指数窗截取能够获得更好的参数跟踪结果。结合两种方法，本研究提出指数窗截取算法，并应用在递推最小二乘算法之中，称为EWRLS 算法。指数窗截取函数的示意图如图5
所示。
图5
指数窗截取函数
Fig.5
Exponential window interception function
如图5所示，指数窗函数截取数据中N 个数据点，当新观测值进入辨识过程后，指数窗向前移动1位，将最后1位数据舍弃，同时，对所有观测值和辨识向量中的数值进行指数加权。加权后，新观测数据、新输入真值和新输入数据的比重为1，应用指数窗截断算法后，EWRLS 辨识算法的计算式为
A i +N ，i +1=
P i +N ，i h T (N +1)[λ+h (N +1)P i +N ，i h T (N +1)]
-1（21）P i +N ，i +1=1
λ
(P i +N ，i -A i +N ，i +1h (N +1)P i +N ，i )（
22）
ξi+N，i+1=ξi+N，i+A i+N，i+1[z(N+1)-h(N+1)ξi+N，i]
（23）式中，λ为指数函数的指数，一般λ在0~1之间。λ越小意味着旧数据对估计结果的影响降低，新数据的影响加大；λ过小可能会导致噪声干扰影响过大，影响估计精度。因此，λ的取值范围一般在0.9~0.99之间。
佛像摆放2.3瞬时频率计算
对连续信号s（t）进行Hilbert变换，可以求取信号的瞬时相位角，Hilbert变换的频域滤波器为
H(f)={-j f>0
j f<0（24）信号的Hilbert变换通过Hilbert滤波器在频域与信号的频域表达相乘实现，即
s
H(t)=F-1(H(f)S(f))（25）信号s（t）和s H（t）共同构成信号s（t）的解析信号z（t）。将解析信号z（
t）用欧拉方程写作指数形式为z(t)=a(t)e jϕ()t（26）式中，a（t）为解析信号的幅度；ϕ（t）为解析信号的瞬时相位角。信号s（t）的瞬时频率定义为
f(t)=12πdϕ(t)d t（27）3信号特征提取算法
由于实测齿轮传动系统振动信号中必定含有除输入、输出齿轮转动引起的振动以及啮合振动以外的其他振动成分，因此，需要利用信号分解重构算法提取测试信号中与所需物理量有关的信号成分。由于齿轮的理论转动频率以及齿轮的啮合频率都可以通过输入轴转速进行计算，因此，可以使用分解信号所在的频带提取所需的信号成分。结合文献[10]41-44和文献[11]1-24的有关方法，本文中使用EMD算法对信号进行分解和重构。原始信号s（t）经EMD分解后可以得到多个本征模态函数（IMF）为
s(t)=∑i=1n imf i(t)+r(t)（28）式中，r（t）为分解余项。EMD分解过程一般是从高频到低频逐步提取信号的IMF分量，各IMF具有不同频率，频率最低的分量（余项）一般认为是信号的趋势项[20]。因此，可以根据各IMF所在的频带，选择相应的IMF函数。利用式（24）~式（27），可以求取各个IMF函数的瞬时频率为
f
i(t)=1
2π
d{}
arctané
ë
ê
ù
û
ú
imf H(t)
imf(t)
d t（29）
在时间-T内IMF函数的平均频率为
fˉi=1-T∫0-T f i(t)d t（30）
假设在正常工作状态下齿轮传动系统输入齿轮
和输出齿轮的理论旋转角速度分别为ωˉi和ωˉo，齿轮
服务器教程
的啮合频率为f m，则选择含对应频率成分的IMF的方
法为
n
i=in
|
|2πfˉi-ωˉi（31）
n
j=in
|
|2πfˉj-ωˉo（32）
n
k=in
||fˉk-f m（33）
啮合刚度主要影响输出齿轮的角速度，造成输
出齿轮旋转角速度在一定范围内波动，波动频率与啮
合频率一致[21]。利用这一特征，可以构造输入齿轮和
输出齿轮的振动特征信号，即将平均频率接近输入
齿轮转频的IMF作为输入齿轮振动特征信号，平均
后赤壁赋翻译
电脑主板频率接近输出齿轮转频的IMF与平均频率接近齿轮
啮合频率的IMF之和作为输出齿轮振动特征信号。
4齿轮啮合刚度识别
4.1仿真分析
假设主动齿轮转速为1800r/min，齿轮副的物理有头有脑
参数如表1所示，啮合刚度K（t）的表达式为
K(t)=106+2×105sin(2π×300t)（34）
表1齿轮副参数
Tab.1Parameters of gear pair
物理量
从动齿轮转动惯量/（kg∙m2）
主动齿轮节圆半径/mm
从动齿轮节圆半径/mm
符号
J2
r1
r2
参数值
100
160
天上一个月亮
320
使用4阶龙格-库塔积分算法，结合式（1），可
以通过数值积分[22]得到主动齿轮旋转后从动齿轮的
旋转角速度，各齿轮的瞬时角速度曲线如图6
所示。
图6主、从动齿轮旋转角速度曲线
Fig.6Rotating angle speed curve of driving and driven gear
如图6所示，由于主、从动齿轮的节圆半径比为
1∶2，因此，从动齿轮的转速为主动齿轮的一半。对从
动齿轮角速度曲线进行频谱分析，结果如图7所示。
