一种未知转速下的RV减速器故障诊断方法

A Fault Diagnosis Method for RV Reducer Under Unknown Rotating Speed

摘要

随着制造业的快速发展， 工业机器人已经成为智能制造的核心执行单元。针对旋转矢量（Rotary vector， RV）减速器在运行过程中转速时变、 瞬时转速不易被读取的问题， 提出一种基于电流信号和振动信号的RV减速器健康状态评估方法。基于电流信号的瞬时频率获取输入转频信息， 计算得出减速器组成部件的故障特征频率； 同步截取平稳振动信号， 再利用粒子寻优的变分模态分解算法对提取出来的振动特征进行分解， 并基于多维评价函数选取故障敏感分量进行包络分析， 以挖掘故障信息进而评估其健康状态。最后， 通过实验数据集验证该方法的有效性。结果表明， 本文所提方法可以获得工业机器人在作业时的关节输入转速， 同时， 多维评价函数可以从RV减速器信号中提取富含冲击信息的故障敏感分量。

Abstract

With the rapid development of the manufacturing industry， industrial robots have become the core execution unit of intelligent manufacturing. Aiming at the problem that the input speed of rotary vector （RV） reducer is time varying and not easy to be read during operation， a health state evaluation method of RV reducer is proposed with the combination of current signal and vibration signal. Based on the instantaneous frequency of the current signal， the input rotation frequency information is obtained. The fault characteristic frequency of the reducer components is calculated， and the stationary vibration signal is intercepted synchronically. The extracted vibration characteristics are decomposed using the variational mode decomposition algorithm of particle optimization， and the fault sensitive component is selected based on the multidimensional evaluation function for envelope analysis. Fault information is mined to evaluate the health status of the fault， and the effectiveness of the method is verified by the experimental data set. The results show that the proposed method can obtain the joint input speed of the industrial robot during operation. The multi-dimensional evaluation function can extract the fault sensitive component rich in shock information from the RV reducer signal.

Graphical abstract

关键词

未知转速 / 工业机器人 / 旋转矢量（Rotary vector， RV）减速器 / 故障诊断 / 多维评价函数

Key words

unknown rotating speed / industrial robots / rotary vector （RV） reducer / fault diagnosis / multidimensional evaluation function

0 引 言

工业机器人是一种复杂的机电系统， 一般由伺服驱动器、 电动机、 减速器、 连杆和控制单元组成。减速器是机器人系统中最昂贵的部件之一， 约占总成本的35%^［1］。作为电机和连杆的中间连接， 减速器降低了电机的输出速度， 并向连杆提供了更大的扭矩。在大多数工业机器人中， 用于反馈控制的传感器数据通常是从执行器侧捕获的， 如电机的旋转角度和输出转矩。然而， 减速器不在反馈控制回路中。 当减速器出现损伤退化时， 控制回路不能有效地对其进行补偿。此外， 减速器还面临着恶劣的工作环境， 如重载和不稳定的转速， 这大大增加了损坏的风险。因此， 减速器的健康预测问题在学术界和工业实践中越来越被重视。

工业机器人常用的减速器有两种， 即谐波减速器和旋转矢量（Rotary vector， RV）减速器。谐波减速器的齿隙为零， 但在重载条件下具有一定的灵活性， 一般在轻载工况下使用。RV减速器不能提供零隙， 但具有更高的刚度， 适合在重载条件下使用。从结构上看， RV减速器是一个复杂的齿轮传动系统。齿轮表面的任何磨损或点蚀都会改变传动系统的动态行为， 从而导致工业机器人的抖动。RV减速器的故障不仅会影响工业机器人的运行精度， 还会在电机反馈控制回路中引入随机干扰。

针对工业机器人的故障诊断方法， 已有文献提出了不同的诊断方法， 如基于物理模型的方法和基于数据驱动的方法。基于物理模型的方法试图构建运动学^［2］、 动力学^［3］或功耗模型^［4］来描述工业机器人在不同健康状态下的物理特性。考虑到工业机器人复杂的结构和不同部件之间的相互作用， 构建一个完善的工业机器人模型是困难的， 甚至是不现实的。

传统的基于模型的研究大多忽略了减速器的灵活性^［5‐7］。他们假设工业机器人是由刚性部件组成的， 这些刚性部件的连杆位置通过减速器的特定传动比与电机的转动同步连接。在这些研究中， 由于减速器故障引起的转矩波动只能作为外部关节转矩的扰动来处理。McIntyre等^［8］采用了类似的策略， 将不同类型的关节故障简化为动态模型中的随机转矩波动， 通过监测转矩波动的变化来检测故障。Capisani等^［9］提出了用于工业机器人故障诊断的高阶滑模观测器， 包括用于检测执行器故障的未知输入观测器和用于检测传感器故障的广义观测器。

一些研究考虑到减速器的灵活性， 假设电机和连杆通过一个柔性关节连接^{［10‐12］}。在这种情况下， 电机侧输出与连杆侧输入之间的位置、 速度和转矩误差可以在动态模型中表述。Bittencourt等^［13］构建了一个扩展的摩擦模型来描述转速、 载荷、 温度和磨损对RV减速器摩擦的影响， 并在此基础上开发了等速水平下的磨损估计器。

基于物理模型的研究主要关注机器人系统的整体动力学行为， 将复杂的减速器结构粗略地建模为电机和连杆之间的简单关节。Qian等^［14］分别建立了行星轮、 摆线轮和曲轴滚针轴承失效模式的极限状态函数， 并将其耦合为RV减速器的多失效模型进行时变可靠性分析。

由于RV减速器在结构上较为复杂， 具有特殊的结构和动态特性， 这也导致其动态响应特性复杂， 在RV减速器的振动信号分析过程中存在3点问题。1） 机器人在执行特定作业任务过程中， 不易直接读取各个关节的转速信息， 导致无法计算关节组件的特征频率及其故障频率。2） 机器人在执行特定作业任务过程中， 运动模式需要经历加速、 匀速、 减速3个阶段， 即需要考虑转速对振动信号分析的影响。3） 频谱分析结果中频率成分相对复杂， 需要考虑如何从减速器故障信号中提炼出故障敏感分量， 进而挖掘出故障信息。

针对现阶段的技术特点， 本文提出一种未知转速下的减速器故障诊断方法， 以克服RV减速器在运行过程中输入转速不易读取的问题， 从而对其健康状态进行监测， 以提高设备的可靠性和精确性。

1 RV减速器故障诊断方法

1.1 诊断方法介绍

本文所述的故障诊断方法， 需要对数据进行采集和处理， 利用采集到的数据进行健康状态评估， 详细步骤如图 1 所示。

首先， 通过对采集到的RV减速器的振动信号与电流信号进行时频变换， 获得时频矩阵， 并同步采集RV减速器运行过程中的振动信号\mathit{x}和电流信号\mathit{s}。\mathit{x}\in {\mathbf{R}}^{n\times \mathrm{1}},\mathit{s}\in {\mathbf{R}}^{n\times \mathrm{1}}， 每组信号包含n个数据点， 且n=T{f}_{\mathrm{s}}， 其中T表示采样时长， f_{\mathrm{s}}表示采样频率。同时， 计算电流信号的解析信号。 z\left(t\right)为信号s(t)的解析信号， 此解析信号是由希尔伯特变换得到的， 即z\left(t\right)=s\left(t\right)+\mathrm{j}H\left[s\left(t\right)\right]， 其中， H\left[s\left(t\right)\right]表示对信号s\left(t\right)做的希尔伯特变换。

然后， 利用脊线提取算法获取电流信号的瞬时频率， 通过观察时频分布矩阵当中第一个时间点对应的峰值\mathit{T}\mathit{F}\left(t_{\mathrm{1}},f_{\mathrm{1}}\right)， 找出局部峰值对应的频率{\mathit{f}}_{\mathrm{1}}=[f_{\mathrm{11}},f_{\mathrm{12}},f_{\mathrm{13}},\cdots ]， 即为瞬时频率曲线提取的起始位置。同时， 同步截取转速平稳时对应的振动信号。其次， 对平稳振动信号进行自适应变分模态分解。最后， 利用多维评价函数筛选故障敏感分量， 并进行包络解调。

1.2 诊断方法的理论基础

可将工业机器人加速和减速阶段的速度近似为线性变化， 对于处理线性调频信号以及求取此类信号的瞬时频率和时频分布， 线性调频小波变换是一种较有效的时频分析方法^［15］。对采集到的信号进行分析， 根据线性调频小波变换的定义计算窗函数

\psi \left(t,t_{\mathrm{0}},\alpha ,\sigma \right)={\omega }_{\left(\sigma \right)}(t-t_{\mathrm{0}})\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{j}\frac{\alpha }{\mathrm{2}}(t-t_{\mathrm{0}}{)}^{\mathrm{2}}\right),（1）

式中： t_{\mathrm{0}},\alpha \in \mathbf{R}分别表示时间和调频率， 其中， t_{\mathrm{0}}表示时间窗在滑动时时间所在的中心位置。当调频率参数\alpha =\mathrm{0}时， 线性调频小波变换就退变为短时傅里叶变换； 当\alpha \ne \mathrm{0}时， 时间宽度不变， 而被截取部分的频带宽度得到降低， 由此可以改善时频分布的聚集性和能量集中程度。w\in L^{\mathrm{2}}\left(R\right)定义了一个非负对称的标准化实窗， 通常是高斯窗， 且高斯窗的表达式为w_{\left(\sigma \right)}\left(t\right)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}\sigma }\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{t}{\sigma }\right)}^{\mathrm{2}}\right)。

从线性调频小波变换的定义式分析可得， 线性调频小波变换的定义可以被解释为由解析信号与窗函数\psi \left(t,t_{\mathrm{0}},\alpha ,\sigma \right)相乘， 再对其乘积求短时傅里叶变换， 即实质为信号加调频窗的傅里叶变换。调频窗是一种窗函数， 其特殊之处在于其频率会随时间而变化， 即呈现调频的特性。在时域信号上使用调频窗后， 窗口的频率会随着时间而变化， 这意味着对于信号的不同时间段， 窗口的频谱特性也会有所不同， 在时频域上可准确地反映信号的频率变化情况。因此， 线性调频小波变换可以计算最终的时频矩阵， 即

\mathit{T}\mathit{F}\left(t_{\mathrm{0}},\omega ,\alpha ,\sigma \right)=A(t_{\mathrm{0}}){\displaystyle \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}}\overline{z}(t){\omega }_{\left(\sigma \right)}(-t_{\mathrm{0}})\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{j}\omega t)\mathrm{d}t,（2）

式中： \overline{z}(t)为解析信号z\left(t\right)在时频平面经过旋转、 平移变化后的信号， 其具体计算公式为

\left\{\begin{array}{l}\overline{z}\left(t\right)=z\left(t\right){\Phi }_{\alpha }^R\left(t\right){\Phi }_{\alpha }^M\left(t,t_{\mathrm{0}}\right),\\ {\Phi }_{\alpha }^R\left(t\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{j}{\displaystyle \sum _{k=\mathrm{2}}^{n+\mathrm{1}}}\frac{\mathrm{1}}{k}{\alpha }_{k-\mathrm{1}}{t}^k\right),\\ {\Phi }_{\alpha }^M\left(t,t_{\mathrm{0}}\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{j}{\displaystyle \sum _{k=\mathrm{2}}^{n+\mathrm{1}}}\frac{\mathrm{1}}{k}{\alpha }_{k-\mathrm{1}}{t}_{\mathrm{0}}^{k-\mathrm{1}}t\right),\\ A\left(t_{\mathrm{0}}\right)=\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{j}\alpha t_{\mathrm{0}}^{\mathrm{2}}/\mathrm{2}),\end{array}\right.（3）

式中： {\Phi }_{\alpha }^R\left(t\right)表示频率旋转算子； {\Phi }_{\alpha }^M\left(t,t_{\mathrm{0}}\right)代表频率平移算子。旋转算子是将解析信号z\left(t\right)在时频平面旋转\theta度， \mathrm{t}\mathrm{a}\mathrm{n}\theta =-\alpha； 平移算子是将t_{\mathrm{0}}处的频率分量向上平移\alpha t_{\mathrm{0}}；A\left(t_{\mathrm{0}}\right)是一个复数， 且模值\left|A\left(t_{\mathrm{0}}\right)\right|=\mathrm{1}。

时频脊线提取技术可以通过在时频图像中的时间维度作“切片”处理， 从而获得瞬时频率估计结果^［16］， 本文采用了基于代价函数自适应搜索脊线算法。首先， 利用代价函数自适应搜寻脊线的瞬时坐标， 目的是让这些坐标使得代价函数值最小。第i个谐波分量的计算公式为

\tilde{f}\left(t_k\right)=\mathrm{m}\mathrm{i}\mathrm{n}\{-{[\mathit{T}\mathit{F}\left(t_k,f\left(t_k\right)\right)]}^{\mathrm{2}}+

\lambda {[f\left(t_k-\tilde{f}\left(t_{k-\mathrm{1}}\right)\right)]}^{\mathrm{2}}\},（4）

\tilde{f}\left(t_{k-\mathrm{1}}\right)-p_{\mathrm{0}}\le f\left(t_k\right)\le \tilde{f}\left(t_{k-\mathrm{1}}\right)+p_{\mathrm{1}},（5）

式中： \tilde{f}\left(t_{\mathrm{1}}\right)取值为f_{\mathrm{1}i}；\tilde{f}\left(t_k\right)为在t_k时刻搜索的频率； \tilde{f}\left(t_{k-\mathrm{1}}\right)-p_{\mathrm{0}}和\tilde{f}\left(t_{k-\mathrm{1}}\right)+p_{\mathrm{1}}分别为脊线频率搜索的下限和上限； p_{\mathrm{0}},p_{\mathrm{1}}值需提前设置； k=\mathrm{1,2},\cdots ,N。 连接各个瞬时脊线坐标点， 最终提取出目标脊线， 以此来估计振动信号的瞬时频率 F。在得到瞬时频率之后， 可以计算电机转频。根据瞬时频率估计结果截取转速平稳段{\mathit{F}}_{\mathrm{1}}， 根据f_{\mathrm{m}}=\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}({\mathit{F}}_{\mathrm{1}})可以计算减速器驱动电机的电流频率， 根据驱动电机的工作原理， 电机转频和电流频率之间存在如式（6）的关系。

f_{\mathrm{s}}=f_{\mathrm{m}}/p,（6）

式中： f_{\mathrm{s}}表示电机转频； f_{\mathrm{m}}表示电流频率； p表示电机极对数。

利用电机转频计算各部件特征频率。为方便各组件的齿数和转频的表示， 首先， 对其进行了标号： 1-太阳轮； 2-行星轮； 3-前端行星架； 4-曲柄轴； 5-摆线轮； 6-针轮； 7-后端行星架。继而计算减速器各组件的特征频率， 分别为

\begin{array}{c}{f}_{\mathrm{2}}=\frac{z_{\mathrm{1}}{z}_{\mathrm{6}}}{\left(z_{\mathrm{5}}-z_{\mathrm{6}}\right)\left(z_{\mathrm{1}}+z_{\mathrm{2}}{z}_{\mathrm{6}}\right)}{f}_{\mathrm{1}},\\ f_{\mathrm{4}}=\frac{z_{\mathrm{1}}{z}_{\mathrm{6}}}{\left(z_{\mathrm{5}}-z_{\mathrm{6}}\right)\left(z_{\mathrm{1}}+z_{\mathrm{2}}{z}_{\mathrm{6}}\right)}{f}_{\mathrm{1}},\\ f_{\mathrm{5}}=\frac{z_{\mathrm{1}}}{z_{\mathrm{1}}+z_{\mathrm{2}}{z}_{\mathrm{6}}}{f}_{\mathrm{1}},f_{\mathrm{6}}=\mathrm{0},\\ f_{\mathrm{7}}=\frac{z_{\mathrm{1}}}{z_{\mathrm{1}}+z_{\mathrm{2}}{z}_{\mathrm{6}}}{f}_{\mathrm{1}},f_{\mathrm{1}\mathrm{c}}=\frac{z_{\mathrm{1}}{z}_{\mathrm{2}}{z}_{\mathrm{6}}}{z_{\mathrm{1}}+z_{\mathrm{2}}{z}_{\mathrm{6}}}{f}_{\mathrm{1}},\\ f_{\mathrm{1}\mathrm{f}}=\frac{f_{\mathrm{1}\mathrm{c}}}{z_{\mathrm{1}}}{N}_{\mathrm{2}},f_{\mathrm{2}\mathrm{f}}=\frac{f_{\mathrm{1}\mathrm{c}}}{z_{\mathrm{2}}}\end{array}，

式中： f=\frac{n}{\mathrm{60}}=\frac{\mathrm{2}\mathrm{\pi }\omega }{\mathrm{60}}， 表示回转频率， Hz； n表示转频， r/min； \omega表示转动角速度， rad/s； N_{\mathrm{2}}表示行星轮个数； 太阳轮回转频率等于输入转频， f_{\mathrm{1}}=f_{\mathrm{s}}；f_{\mathrm{1}\mathrm{c}}表示一级传动啮合频率； f_{\mathrm{1}\mathrm{f}}表示太阳轮故障频率； f_{\mathrm{2}\mathrm{f}}表示行星轮故障频率。

为了从减速器故障信号中提炼出高质量的故障敏感分量， 利用变分模态分解^{［17‐19］}算法提取故障信息。在利用变分模态分解算法时， 需要考虑分解参数的设置^［20］。首先， 同步截取转速平稳时对应的振动信号X_{\mathrm{1}}， 初始化粒子群算法的参数， 包括寻优最大迭代次数K_{\mathrm{m}\mathrm{a}\mathrm{x}}、 种群规模\left({\alpha }_{\mathrm{1}},K_{\mathrm{1}}\right)和初始粒子的迁移速度\left(v_{\alpha },v_K\right)。这里设置K_{\mathrm{m}\mathrm{a}\mathrm{x}},{\alpha }_{\mathrm{1}}=\mathrm{100},K_{\mathrm{1}}=\mathrm{1},v_{\alpha }=v_K=\mathrm{1}。然后， 使用变分模态分解方法和粒子群优化算法中可用的初始化参数来计算分解分量的适应度函数f_K。

根据式（7）计算种群中所有粒子的当前位置和速度， 如果当前适应度函数取值小于历史最小适应度函数取值， 则种群中粒子的当前位置和速度根据式（7）和式（8）进行更新。当K=K_{\mathrm{m}\mathrm{a}\mathrm{x}}时， 寻优过程终止， 输出最优适应度取值f_{\mathrm{1}:K}^{\mathrm{m}\mathrm{i}\mathrm{n}}以及变分模态分解的最优参数组合\left(\widehat{\alpha },\widehat{K}\right)。

{\widehat{\alpha }}_{\mathrm{i}\mathrm{d}}=\widehat{\alpha },{\widehat{K}}_{\mathrm{i}\mathrm{d}}=\widehat{K},（7）

\left[\begin{array}{c}\widehat{\alpha }\\ \widehat{K}\end{array}\right]=\delta \bullet \left[\begin{array}{c}{\alpha }_{\mathrm{i}\mathrm{d}}\\ K_{\mathrm{i}\mathrm{d}}\end{array}\right]+C_{\mathrm{1}}\bullet Rand\left(\mathrm{0,1}\right)\bullet \left({\widehat{\alpha }}_{\mathrm{i}\mathrm{d}}-{\alpha }_{\mathrm{i}\mathrm{d}}\right)+C_{\mathrm{2}}\bullet Rand(\mathrm{0,1})\bullet ({\widehat{K}}_{\mathrm{i}\mathrm{d}}-K_{\mathrm{i}\mathrm{d}}),（8）

式中： \delta表示惯性因数； C_{\mathrm{1}}和C_{\mathrm{2}}控制粒子的移动速度； Rand（0，1）表示生成0~1之间的随机数； {\alpha }_{\mathrm{i}\mathrm{d}}和K_{\mathrm{i}\mathrm{d}}表示种群的当前位置； {\widehat{\alpha }}_{\mathrm{i}\mathrm{d}}和{\widehat{K}}_{\mathrm{i}\mathrm{d}}记录了历史最优位置； 设置\delta =\mathrm{1},C_{\mathrm{1}}=C_{\mathrm{2}}=\mathrm{0.5}。

使用最优参数组合对平稳振动信号进行变分模态分解， 获得多个模态分量{\left\{u_k\right\}}_{k=\mathrm{1}}^{\widehat{K}}。

利用多维评价函数筛选故障敏感分量， 并进行包络解调。基于模态分量集合{\left\{u_k\right\}}_{k=\mathrm{1}}^{\widehat{K}}计算峭度指标（S）、 包络谱峭度指标（V）和包络熵指标（E）， 公式为

S_i=\frac{E{\left(u_i-{\mu }_{x_i}\right)}^{\mathrm{4}}}{{{\sigma }_{x_i}}^{\mathrm{4}}},（9）

V_i=\frac{E{\left(a_i-{\mu }_{a_i}\right)}^{\mathrm{4}}}{{{\sigma }_{a_i}}^{\mathrm{4}}},（10）

\left\{\begin{array}{c}{E}_i=-{\displaystyle \sum _{j=\mathrm{1}}^N}{p}_{i,j}\mathrm{l}\mathrm{o}\mathrm{g}{p}_{i,j},\\ p_{i,j}=a_i(j)/{\displaystyle \sum _{j=\mathrm{1}}^N}{a}_i(j),\end{array}\right.（11）

式中： u_i表示分解后的第i个信号分量； a_i表示信号分量u_i经希尔伯特变换后得到的包络信号； p_{i,j}为a_i\left(j\right)的归一化形式； E(\bullet )表示期望函数； \mu和\sigma分别为数据序列的均值和方差。

为了进一步衡量各个指标的重要性， 运用熵权法对上述评价指标进行加权， 最终的多维评价函数构造流程如下。

1） 标准化。根据指标与衡量性能之间的正、 负相关性， 将指标标准化过程分为

{\widehat{S}}_i=\frac{S_i-\mathrm{m}\mathrm{i}\mathrm{n}\left(S\right)}{\mathrm{m}\mathrm{a}\mathrm{x}\left(S\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(S\right)},（12）

{\widehat{V}}_i=\frac{V_i-\mathrm{m}\mathrm{i}\mathrm{n}\left(V\right)}{\mathrm{m}\mathrm{a}\mathrm{x}\left(V\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(V\right)},（13）

{\widehat{E}}_i=\frac{\mathrm{m}\mathrm{a}\mathrm{x}\left(E\right)-E_i}{\mathrm{m}\mathrm{a}\mathrm{x}\left(E\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(E\right)}\mathrm{。}（14）

2） 计算权重。以峭度指标为例， 其余指标计算过程相同， 权重计算公式为

s_i=\frac{{\widehat{S}}_i}{{\displaystyle \sum _{i=\mathrm{1}}^K}{\widehat{S}}_i},（15）

\widehat{s}=-\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^K}{s}_i\mathrm{l}\mathrm{n}{s}_i,（16）

w_s=\frac{\mathrm{1}-\widehat{s}}{\mathrm{3}-\left(\widehat{s}+\widehat{v}+\widehat{e}\right)}\mathrm{。}（17）

3） 多维评价函数。获取权重后利用线性加权求和计算多维评价函数， 即

L_i={\omega }_{\mathrm{s}}{\widehat{S}}_i+{\omega }_{\mathrm{v}}{\widehat{V}}_i+w_{\mathrm{e}}{\widehat{E}}_i,（18）

式中： w_{\mathrm{s}},w_{\mathrm{v}},w_{\mathrm{e}}表示从第i个信号分量获取的峭度指标、 包络谱峭度指标和包络熵指标的权重。计算多维评价函数得分最高的分解分量u_l， 并对其进行包络解调分析获得包络谱e_l， 通过包络谱成分分析， 进而对减速器进行故障诊断。

2 RV减速器摆臂实验台数据集验证

2.1 实验介绍

为了验证诊断方法的可行性， 利用RV减速器实验台数据对本文方法进行验证。实验台主要由伺服电机、 RV减速器和摆臂组成， 伺服电机负责提供动力输入， 通过驱动减速器从而带动摆臂进行摆动， 摆臂末端的质量块用于模拟机器人实际执行任务时关节所受的载荷。实验台使用Net0824数据采集卡进行数据采集， 通过在减速器顶部布置三向振动加速度传感器、 在电机控制柜布置电流传感器、 在电机和减速器侧面布置温度传感器， 来采集减速器的三向振动信号、 电机的三相电流信号以及电机和RV减速器的温度信号。实验设置减速器输出轴运动方式为往复摆动（以竖直方向为0°， 摆动位置区间为［－90°， 90°］）， 通过拆装质量块， 将负载设置为24 kg， 摆动速度设置为45 （°）/s和90 （°）/s。实验模拟了行星轮和太阳轮的不同健康状态， 包括正常和点蚀。每种工况下进行若干组实验， 每组连续采集2 min， 采样频率为6 250 Hz。

2.2 实验数据处理及分析

首先， 利用线性小波变换得到电流信号的时频图。然后， 利用基于代价函数自适应搜寻脊线算法提取出电流信号的瞬时频率。基于电流信号估计瞬时频率结果， 见图 2。该算法可以完整地提取出瞬时频率的变化趋势， 通过分析提取的瞬时频率的变化趋势， 可以推断出工业机器人关节在往复摆动的过程中存在加速、 匀速和减速3个阶段。进一步分析图 2， 可以发现匀速运动阶段关节电机的电流频率为151.4 Hz。

图 2 表明， 通过估计电流信号的瞬时频率， 可以获得工业机器人实际工作时的关节转速信息， 解决了实际工业场景中工业机器人关节转速难以获取的问题， 为机器人关节信号的进一步分析奠定了基础。根据关节电机电流信号的瞬时频率估计结果， 可以计算 RV减速器各组成部分的特征频率， 具体计算结果见表 1。

在得到减速器各部分特征频率后， 同时截取转速稳定段的振动信号， 如图 3 所示。

为提取更加明显的故障特征信息， 需要对转速平稳段对应的振动信号进行变分模态分解， 并利用多维评价函数筛选出故障敏感的分量。

首先， 对变分模态分解算法中的两个参数， 利用粒子群优化算法进行了寻优， 算法参数的初始化设置信息需要对平稳振动信号进行自适应变分模态分解。粒子群优化算法的适应度函数取值随种群进化迭代次数变化的曲线如图 4 所示。可以看出随着算法的不断进化迭代， 适应度函数取值不断减小， 随后自第 15 次迭代后开始收敛， 此时， 最小包络熵取值为8.526 3， 寻优到的最佳参数组合为（9， 4 937）。 然后， 根据寻优结果设定变分模态分解算法中的惩罚因子为4 937和分量个数为9， 对平稳振动信号进行变分模态分解后得到9个信号分量及其频谱结果， 分别如图 5 和图 6 所示。

为了筛选出富含故障信息的IMF， 采用包络熵指标、 峭度指标与包络谱峭度指标对每个IMF中包含故障信息的量进行评价， 运用熵权法求得3个指标的权重分别为 0.49， 0.27， 0.23， 对3个指标加权求和后最终可获得各个IMF的多维评价指标结果， 见图 7。

从图 7 可以看出IMF5得分最高， 说明IMF5中包含故障信息最多。通过对比IMF5时域信号和IMF8时域信号（见图 8）， 可以发现IMF5包含更丰富的冲击信息。说明根据故障信号的脉冲性和循环平稳性而构造的多维评价函数， 在故障敏感分量提取上具有较好的适用性。

对选取的故障敏感分量IMF5进行包络分析， 结果如图 9 所示， 可以发现， 行星轮故障特征频率 （10 Hz） 及其倍频成分较为明显。因此， 基于本文所提方法可以对减速器的故障类别进行辨识。

3 结 论

针对RV减速器在运行过程中输入转速不易被读取的问题， 提出了一种未知转速下减速器故障的诊断方法。一方面， 通过对电流信号的瞬时频率进行估计， 可以获得驱动电机的转速， 根据瞬时频率估计结果获取RV减速器各组件的特征频率。另一方面， 提出了一种RV减速器的故障特征自适应提取方法， 该方法对转速平稳段对应的振动信号进行变分模态分解， 其中分解参数由粒子群优化算法自适应获取。随后， 使用多维评价函数完成对故障敏感分量的筛选， 并进行包络分析以提取故障信息， 有效解决了机器人关节在运行过程中转速信息不易读取和关节振动信号分析中存在变转速等问题。为了精确地提取故障信号成分， 未来研究工作考虑对时频图精细化估计以及时频脊线提取算法进行优化， 并改进多维评价函数指标以筛选富含故障成分的信号成分。

参考文献

[1] Global and china industrial robot speed reducer industry report［R/OL］.（2020-09-01）［2023-01-10］.

[2] ZHANG H W， AHMAD S， LIU G J. Torque estimation for robotic joint with harmonic drive transmission based on position measurements［J］.IEEE Transactions on Robotics， 2015， 31（2）： 322-330.

[3] DE L A， MATTONE R. Actuator failure detection and isolation using generalized momenta ［C］//IEEE International Conference on Robotics and Automation. 2003， 1： 634-639.

[4] SABRY A H， NORDIN F H， SABRY A H， et al. Fault detection and diagnosis of industrial robot based on power consumption modeling［J］. IEEE Transactions on Industrial Electronics， 2020， 67（9）： 7929-7940.

[5] HADDADIN S， DE L A， ALBU-SCHAFFER A. Robot collisions： A survey on detection， isolation， and identification ［J］. IEEE Transactions on Robotics， 2017， 33（6）： 1292-1312.

[6] CACCAVALE F， VILLANI L. Fault diagnosis for industrial robots ［C］//Fault Diagnosis and Fault Tolerance for Mechatronic Systems： Recent Advances. 2003： 85-108.

[7] MA H J， YANG G H. Simultaneous fault diagnosis for robot manipulators with actuator and sensor faults［J］. Information Sciences， 2016， 366： 12-30.

[8] MCINTYRE M L， DIXON W E， DAWSON D M， et al. Fault identification for robot manipulators［J］. IEEE Transactions on Robotics， 2005， 21（5）： 1028-1034.

[9] CAPISANI L M， FERRARA A， FERREIRA D L A， et al. Manipulator fault diagnosis via higher order sliding-mode observers［J］. IEEE Transactions on Industrial Electronics， 2012， 59（10）： 3979-3986.

[10] DE L A， BOOK W J. Springer handbook of robotics［M］. Berlin： Springer， 2016.

[11] QIAO Y T， CAO J Y， HUANG G H， et al. Improved subspace modal identification of industrial robots［J］. Journal of Field Robotics， 2023， 40（6）： 1327-1338.

[12] LIU H， LEI Y G， YANG X， et al. Deflection estimation of industrial robots with flexible joints［J］. Fundamental Research， 2022， 2（3）： 447-455.

[13] BITTENCOURT A C， AXELSSON P. Modeling and experiment design for identification of wear in a robot joint under load and temperature uncertainties based on friction data［J］. IEEE/ASME Transactions on Mechatronics， 2013， 19（5）： 1694-1706.

[14] QIAN H M， LI Y F， HUANG H Z. Time-variant reliability analysis for industrial robot RV reducer under multiple failure modes using Kriging model［J］. Reliability Engineering & System Safety， 2020， 199： 106936.

[15] YU G， ZHOU Y Q. General linear chirplet transform［J］. Mechanical Systems and Signal Processing， 2016， 70： 958-973.

[16] SHI H Y， QIN Y， WANG Y， et al. An improved viterbi algorithm for adaptive instantaneous angular speed estimation and its application into the machine fault diagnosis［J］. IEEE Transactions on Instrumentation and Measurement， 2021， 70（2）： 1-11.

[17] 王宏民， 刘洪超， 薛萍. 变分模态分解在齿轮故障特征提取中的应用［J］. 东北林业大学学报， 2019， 47（8）： 97-102.

[18] WANG Hongmin， LIU Hongchao， XUE Ping. Application of variational mode decomposition in fault feature extraction of gears［J］. Journal of Northeast Forestry University， 2019， 47（8）： 97-102.

[19] KUMAR A， GANDHI C P， ZHOU Y Q， et al. Latest developments in gear defect diagnosis and prognosis： A review ［J］. Measurement， 2020， 158（1）： 107-126.

[20] 杨宗林， 熊继军. 参数优化VMD在爆破振动信号分析中的应用 ［J］. 中北大学学报（自然科学版）， 2020， 41（5）： 467-473.

[21] YANG Zonglin， XIONG Jijun. Application of parameter optimization VMD in blasting vibration signal analysis［J］. Journal of North University of China（Natural Science Edition）， 2020， 41（5）： 467-473.

[22] LI Z X， JIANG Y， HU C， et al. Recent progress on decoupling diagnosis of hybrid failures in gear transmission systems using vibration sensor signal： A review ［J］. Measurement， 2016， 90（2）： 4-19.