基于运动模式划分的工业机器人健康监测方法

Health Monitoring Method of Industrial Robots Based on Motion Pattern Division

摘要

由于工业机器人存在多运动模式耦合的问题， 传统的健康监测方案需要在每个关节处单独安装传感器， 难以满足实际工业现场需求。本文以6关节工业机器人为研究对象， 基于振动信号研究了多运动模式切换场景下工业机器人的健康监测方法。首先， 通过跳变点算法实现机器人运动模式划分， 获取不同运动模式对应的信号区间。其次， 对不同运动模式的信号分别提取监测指标。最后， 基于控制图法实现工业机器人不同关节的健康监测。在工业机器人退化实验数据中验证了本文所提方法， 表明本文所提方法能够在仅使用两个振动传感器的条件下实现机器人6个关节的健康状态监测。本文所提运动模式划分算法在对大量历史退化数据进行分析时， 所需运行时间更短、 单次精度更高且重复性更好。本文提出的方法能够在使用少量传感器的条件下， 有效避免运动模式耦合和采样信号的差异可能导致的监测结果误判， 使得监测结果更加精准可靠， 适用于实际工业现场的工业机器人健康状态监测。

Abstract

Due to the problem of multi-motion mode coupling of industrial robots， the traditional health monitoring scheme needs to install sensors at each joint separately， which is difficult to meet the needs of actual industrial field. A 6-joints industrial robot is taken as the research object， and the health monitoring method of industrial robot in multi-motion mode switching scenario is studied based on vibration signals. Firstly， robot motion modes are divided by jump point algorithm， and signal intervals corresponding to different motion modes are obtained. Secondly， the monitoring indexes are extracted from the signals of different motion modes. Finally， the health monitoring of different joints of industrial robots is realized based on the control chart method. The proposed method is verified by the experimental data of industrial robot degradation. It shows that the proposed method can realize the health status monitoring of 6 joints of the robot with only two vibration sensors. The motion pattern partitioning algorithm proposed in this paper costs shorter running time， has higher single precision and better repeatability when analyzing a large number of historical degraded data. The proposed method can effectively avoid the misjudgment of monitoring results which may be caused by the difference of motion mode coupling and sampling signal under the condition using only a few sensors， making the monitoring results more accurate and reliable， and more suitable for the health status monitoring of industrial robots in actual industrial sites.

Graphical abstract

关键词

工业机器人 / 健康监测 / 运动模式划分 / 控制图法

Key words

industrial robot / health monitoring / motion mode division / control chart method

0 引 言

随着劳动力的日益短缺及劳动力成本的不断增长， 大量的劳动密集型生产企业需要转型升级。工业机器人由于具有通用性强、 高效可靠、 重复精度高等优点， 在传统制造业等劳动密集型产业的转型升级中发挥着重要的作用， 在企业提高生产效率、 改进产品质量和一致性、 降低企业生产及劳动力成本等方面起到关键作用^［1］。工业机器人作为一种结构精密复杂的机电一体化系统， 集机械工程技术、 电子工程技术、 信息传感器技术以及计算机科学技术于一体， 一旦发生故障， 不仅需要较高的维修保养费用， 甚至会严重威胁生产过程的安全。因此， 亟待开展工业机器人的健康监测研究， 进而全面掌控其健康动态， 以保障生产质量， 提高生产效率^［2］。

针对工业机器人机械传动部件健康监测问题的研究， 主要参考依据包括机器人运动过程中的动力学参数、 电流信号和振动信号等。基于动力学参数的健康监测研究方面， Sneider等^［3］通过建立工业机器人基于观测器的非线性数学模型来获取其运动过程中的信息， 并运用参数估计方法通过与已定参数进行比较对工业机器人进行健康监测。Endo等^［4］基于控制系统的动态模型进行实时健康监测， 通过比较工业机器人控制器的传递特性和各种故障类型的故障模型函数来识别故障。

基于电流信号的健康监测研究方面， Borgi等^［5］基于工业机器人控制器中的三相电流数据， 经过特征提取与特征筛选后， 建立特征与机器人精度误差之间的多元映射关系。Chen等^［6］从6轴工业机器人的控制器中采集反馈电流信号， 并利用连续小波变换将其转换为时频图像， 将去噪的时频图输入紧凑型卷积Transformer中实现工业机器人的复合故障诊断。Jiang等^［7］将最大相关峰度反卷积和完全集成经验模态分解相结合， 对原始电流信号进行预处理， 然后， 从中提取退化指标， 并利用改进的BP神经网络对减速器剩余寿命进行预测。

振动信号作为常规的无损监测技术， 被广泛应用于关节减速器的状态监测与诊断。振动测量和分析被认为是旋转机械最有效的状态诊断方法^［8］。Jaber等^［9］基于振动信号提取方差作为监测指标， 通过对不同故障的工业机器人执行100次拾取动作， 计算不同故障下方差的均值， 再据此使用控制图方法对PUMA 560第 3 关节进行健康监测。Vagaš等^［10］通过加速度传感器采集机器人满载运行时的振动信号， 并计算得到信号的频谱， 通过工业机器人振动的测量和评估振动监测的方法对工业机器人进行健康监测。Bosen等^［11］针对工业机器人关节减速器在不同工况下电流幅值差异大的问题， 提出了基于维纳过程的多任务阶段退化模型， 并根据不同任务阶段设置失效阈值的方法， 较好地描述了谐波减速器在多任务运行条件下的退化过程。Kim等^［12］提出了一种基于相位的时域平均方法， 该方法能够从原始振动信号中提取出含有故障信息的残差信号， 并根据从残差信号中提出的特征来评估机器人减速器的健康状态。Zhi等^［13］利用联合小波区域相关阈值去噪算法对原始振动信号进行去噪， 然后， 利用CNN提取特征， 最后， 通过LSTM网络识别故障类型， 完成了谐波减速器的故障检测。

除了利用振动和电流信号对关节减速器进行状态监测以外， 温度监测^［14］、 红外热成像^［15］和磨粒分析^［16］等都被用来监测关节减速器健康状态。

通过文献调研可以发现， 基于外部振动信号的工业机器人健康监测是通过安装于机器人腕关节处的加速度传感器来采集机器人运行过程中的振动信号数据， 然后， 从信号数据中提取退化指标来反映机器人的健康退化， 研究对象为仅包含单个运动模式的工业机器人。

然而， 对于实际工业场景， 工业机器人在执行特定作业任务过程中往往存在多运动模式的切换。对于传统基于振动信号的监测方法， 一方面， 考虑到对多个腕关节分别布置加速度传感器成本较高； 另一方面， 考虑到多运动模式的耦合以及单次振动信号数据采样样本内所包含运动模式的差异性， 如果直接针对单次振动信号数据采样样本提取退化指标， 将无法实现精确的健康监测功能， 甚至会造成健康状态的误判。

图 1 为工业机器人退化实验过程中两次采集的振动信号样本对比。可以看出， 相同的采样方案设定下（采样时长、 采样间隔、 采样频率等参数均相同）， 由于包含的运动模式不同， 导致不同时刻两次采样所采集的振动信号数据存在明显的差异， 故针对全部样本计算健康指标并进行健康监测不具有实际参考意义。

本文基于振动信号研究了在仅使用少量传感器的条件下， 实现多运动模式切换场景下的工业机器人机械传动部件健康监测方法。首先， 通过联合多物理源信号实现工业机器人运动模式划分， 获取不同运动模式对应的信号区间。其次， 针对运动模式信号分别提取退化指标。最后， 基于控制图法实现了工业机器人不同机械传动部件的健康监测， 以提高工业机器人健康监测的可靠性和精确性。

1 跳变点识别算法

1.1 跳变点识别理论

对于信号x_{\mathrm{1}:n}=(x_{\mathrm{1}},\cdots ,x_n)， 若其中存在单个跳变点x_k， 则从统计学角度满足

（1）

式中： S\left(\bullet \right)为表征信号统计属性的函数， 如二次累加和函数为S(x_{\mathrm{1}:n})={\displaystyle \sum _{i=\mathrm{1}}^n}{‖x_i-{\overline{x}}_{\mathrm{1}:n}‖}_{\mathrm{2}}^{\mathrm{2}}， 其中， {\overline{x}}_{\mathrm{1}:n}表示数据序列x_{\mathrm{1}:n}的均值。

多跳变点问题为单个跳变点问题的延伸。若有m个跳变点以及它们的位置{\tau }_{\mathrm{1}:m}=\left({\tau }_{\mathrm{1}},\dots ,{\tau }_m\right)， 每个跳变点位置都是1到n-\mathrm{1}之间的整数， 定义{\tau }_{\mathrm{0}}=\mathrm{0}和{\tau }_{m+\mathrm{1}}=n， 并且假设跳变点位置按照单调递增排序， 即当时， 。因此， m个跳变点把信号分成m+1段， 第i段信号表示为x_{{\tau }_{i-\mathrm{1}}+\mathrm{1}:{\tau }_i}， 则可获得多跳变点问题的目标损失函数为

F(n)=\underset{\left\{\tau :\mathrm{0}={\tau }_{\mathrm{0}},...,{\tau }_{m+\mathrm{1}}=n\right\}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{\displaystyle \sum _{i=\mathrm{1}}^{m+\mathrm{1}}}\left[S\left(x_{{\tau }_{i-\mathrm{1}}+\mathrm{1}:{\tau }_i}\right)\right]+\beta f(m)\right\},（2）

式中： \beta f\left(m\right)为防止过拟合的惩罚项。

对多跳变点目标损失函数做进一步变换， 即

\begin{array}{c}F(n)=\underset{\left\{\tau :\mathrm{0}={\tau }_{\mathrm{0}},\dots ,{\tau }_{m+\mathrm{1}}=n\right\}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{\displaystyle \sum _{i=\mathrm{1}}^{m+\mathrm{1}}}\left[S\left(x_{{\tau }_{i-\mathrm{1}}+\mathrm{1}:{\tau }_i}\right)+\beta \right]\right\}=\\ \underset{t}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\underset{\left\{\tau :\mathrm{0}={\tau }_{\mathrm{0}},\dots ,{\tau }_m=t\right\}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{i=\mathrm{1}}^m}\left[S\left(x_{{\tau }_{i-\mathrm{1}}+\mathrm{1}:{\tau }_i}\right)+\beta \right]+S\left(x_{t+\mathrm{1}:n}\right)+\beta \right\}=\\ \underset{t}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{F(t)+S\left(x_{t+\mathrm{1}:n}\right)+\beta \right\}\mathrm{。}\end{array}（3）

通过对t=\mathrm{1,2},\cdots ,n进行遍历， 即可通过递推关系， 利用动态规划算法求解跳变点位置集合\mathit{\tau }。

1.2 改进的跳变点识别算法

为了提高程序运行效率， 假设信号在分段内均值变化， 方差恒定。由于实际信号中往往存在随机噪声， 即

x(t_n)=s(t_n)+z_n，n=\mathrm{1,2},\cdots ,N，（4）

式中： s(t_n)为纯净信号； z_n为随机噪声。

假设随机噪声z_n服从具有零均值和给定方差的正态概率分布P\left(z_n\right)， 即

P\left(z_n\right)\mathrm{d}{z}_n=\frac{\mathrm{1}}{{\sigma }_n\sqrt[]{\mathrm{2}\mathrm{\pi }}}{\mathrm{e}}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{z_n}{{\sigma }_n}\right)}^{\mathrm{2}}\mathrm{d}{z}_n}\mathrm{。}（5）

如果分段内信号的均值为\lambda， 且保持不变， 则分段内信号的似然函数L_n为

L_n=\frac{\mathrm{1}}{{\sigma }_n\sqrt[]{\mathrm{2}\mathrm{\pi }}}{\mathrm{e}}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{x_n-\lambda }{{\sigma }_n}\right)}^{\mathrm{2}}}\mathrm{。}（6）

由于分段信号之间相互独立， 故信号的全局似然函数为

L^{\left(k\right)}=\prod _n{L}_n=\frac{{\left(\mathrm{2}\mathrm{\pi }\right)}^{-\frac{N_k}{\mathrm{2}}}}{\prod _m{\sigma }_m}{\mathrm{e}}^{-\frac{\mathrm{1}}{\mathrm{2}}\sum _n{\left(\frac{x_n-\lambda }{{\sigma }_n}\right)}^{\mathrm{2}}}\mathrm{。}（7）

为了方便方程化简和表达， 定义

\begin{array}{c}{a}_k=\frac{\mathrm{1}}{\mathrm{2}}\sum _n\frac{\mathrm{1}}{{\sigma }_n^{\mathrm{2}}},\\ b_k=-\sum _n\frac{x_n}{{\sigma }_n^{\mathrm{2}}}\mathrm{。}\end{array}（8）

通过最大化信号的全局似然函数消除\lambda的影响， 即最大化

-\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \sum _n^{}}\left(\frac{x_n-\lambda }{{\sigma }_n}\right)}^{\mathrm{2}},（9）

从而

{\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\sum _n\frac{x_n}{{\sigma }_n^{\mathrm{2}}}/\sum _{n^{\text{'}}}\frac{\mathrm{1}}{{\sigma }_{n^{\text{'}}}^{\mathrm{2}}}=-b_k/\mathrm{2}{a}_k\mathrm{。}（10）

将式（10）代入到方程（7）中， 可以得到信号的全局对数似然函数

\mathrm{l}\mathrm{o}\mathrm{g}{L}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}=-\frac{\mathrm{1}}{\mathrm{2}}\sum _n{\left(\frac{x_n+\frac{b_k}{\mathrm{2}{a}_k}}{{\sigma }_n}\right)}^{\mathrm{2}}=

-\frac{\mathrm{1}}{\mathrm{2}}\left[\sum _n\frac{x_{{}_n}^{\mathrm{2}}}{{\sigma }_{{}_n}^{\mathrm{2}}}+\frac{b_k}{a_k}\sum _n\frac{x_{{}_n}^{}}{{\sigma }_{{}_n}^{\mathrm{2}}}+\frac{b_{{}_k}^{\mathrm{2}}}{\mathrm{4}{a}_{{}_k}^{\mathrm{2}}}\sum _n\frac{\mathrm{1}}{{\sigma }_{{}_n}^{\mathrm{2}}}\right]\mathrm{。}（11）

由于第一项\frac{x_n^{\mathrm{2}}}{{\sigma }_n^{\mathrm{2}}}是与模型无关的常数， 故将其忽略。 然后， 将式（8）代入式（11）， 由于最小化损失函数等于最大化对数似然函数， 因此， 式（2）的目标函数可简化为

F(n)=\mathrm{l}\mathrm{o}\mathrm{g}{L}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}=\frac{b_k^{\mathrm{2}}}{\mathrm{4}{a}_k},（12）

式中： F(n)为简化后的目标函数； \mathrm{l}\mathrm{o}\mathrm{g}{L}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(k\right)}为分段信号的全局对数似然函数。

在假设分段内信号方差与时间无关的情况下， {\sigma }_n可以作为一个常数， 而不必在每个数据分段中指定{\sigma }_n。改进后的跳变点识别算法为

输入： 信号\mathit{s}， 惩罚项系数\beta；

初始化： 信号长度N=\mathrm{l}\mathrm{e}\mathrm{n}(\mathit{s})； 归一化\widehat{\mathit{s}}=\frac{\mathit{s}-{\mu }_s}{{\sigma }_s}；

for R=\mathrm{1},\cdots ,N do

1） 根据式（12）计算目标函数 A(r)=F(r)-\beta +\left\{\begin{array}{ll}\mathrm{0},& r=\mathrm{1},\\ best(r-\mathrm{1}),& r=\mathrm{2},\cdots ,R\mathrm{。}\end{array}\right.

2） 最大化目标函数求解跳变点位置： r^{opt}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}(A(r))；

3） 保存跳变点位置： best(r)=A(r^{opt})， last(r)=r^{opt}；

end for

c{p}_{\mathrm{1}}=last(N)；

while c{p}_i>\mathrm{1} do

4） 获取分段信号最后一个跳变点： c{p}_i=last(c{p}_{i-\mathrm{1}}-\mathrm{1})；

5） i=i+\mathrm{1}；

end while

输出跳变点集合： cp。

1.3 跳变点识别结果评价指标

为了定量评价跳变点识别算法结果的准确性， 引入了3个评价指标： 精确率（Positive predictive value， PPV）、 召回率（True positive rate， TPR）和F1^［17］。假设目标跳变点集合为{\left\{{\tau }_{ti}\right\}}_{i=\mathrm{1}}^n， 跳变点算法获得的跳变点集合为{\left\{{\tau }_{pi}\right\}}_{i=\mathrm{1}}^m， 根据算法获取的跳变点与两个集合间的所属关系计算相应条件下的跳变点数量， 分别用真正例TP、 假反例FN、 假正例FP和真反例TN来表示， 具体定义过程见表 1。

精确率表示跳变点算法获得的跳变点集中真正属于目标跳变点集的跳变点个数， 占跳变点集总个数的比例， 计算表达式为

PPV=\frac{TP}{TP+FP}\mathrm{。}（13）

召回率表示跳变点算法获得的跳变点集中真正属于目标跳变点集的跳变点个数， 占目标跳变点集总个数的比例， 计算表达式为

TPR=\frac{TP}{TP+FN}\mathrm{。}（14）

F1是准确率和召回率的调和平均数， 计算表达式为

F\mathrm{1}=\frac{\mathrm{2}\cdot PPV\cdot TPR}{PPV+TPR}\mathrm{。}（15）

2 基于运动模式划分的工业机器人健康监测方法流程

基于运动模式划分的工业机器人健康监测方法流程如图 2 所示。

下面详细阐述该方法的实施步骤。

1） 获取不同运动模式的运行时长、 振动信号和电流信号。在布置规划轨迹算法时， 需确定每种运动模式的运行时长。 假设有p种运动模式， 获得的运行时长集合为\mathit{T}=(t_{\mathrm{1}},\dots ,t_i,\dots ,t_p)， 获取的振动信号为{\mathit{y}}_{\mathrm{1}:n}=(y_{\mathrm{1}},\cdots ,y_n)， 电流信号为{\mathit{z}}_{\mathrm{1}:n}=(z_{\mathrm{1}},\cdots ,z_n)。

2） 对采集的振动信号进行降采样。一方面， 是为了减少算法计算量， 另一方面， 是为了增强跳变效应。首先， 将采集的振动信号{\mathit{y}}_{\mathrm{1}:n}=(y_{\mathrm{1}},\cdots ,y_n)用l_{\mathrm{w}\mathrm{i}\mathrm{n}}长度的窗口去切分， 获得m=n/l_{\mathrm{w}\mathrm{i}\mathrm{n}}+\mathrm{1}段信号， 且每段信号表示为{\mathit{y}}_{(k-\mathrm{1})l_{\mathrm{w}\mathrm{i}\mathrm{n}}:k{l}_{\mathrm{w}\mathrm{i}\mathrm{n}}},k=\mathrm{1,2},\cdots ,m， 其中， n为采集的振动信号总长度。然后， 计算每个信号段的均方根值， 作为降采样后的信号表示， 即{\widehat{y}}_k=rms\left({\mathit{y}}_{(k-\mathrm{1})l_{\mathrm{w}\mathrm{i}\mathrm{n}}:k{l}_{\mathrm{w}\mathrm{i}\mathrm{n}}}\right),k=\mathrm{1,2},\dots ,m。

3） 基于降采样后的信号， 使用跳变点识别算法获取跳变点集合。当工业机器人进行多个运动模式切换时， 即降采样信号中有w个跳变点， 并且位置分别是{\mathit{\tau }}_{\mathrm{1}:w}=({\tau }_{\mathrm{1}},\cdots ,{\tau }_i,\cdots ,{\tau }_w)， 其中， {\tau }_{\mathrm{0}}=\mathrm{0},{\tau }_{w+\mathrm{1}}=m， 将降采样信号分成w+1段， 第i段表示为{\widehat{y}}_{{\tau }_{i-\mathrm{1}}+\mathrm{1}:{\tau }_i}， 利用改进后的跳变点识别算法获得跳变点位置集合{\mathit{\tau }}^{\text{'}}=\left({\tau }_{\mathrm{1}}^{\text{'}},\cdots ,{\tau }_c^{\text{'}}\right)， 下标c表示最终获取的跳变点集合中元素总个数， 。

4） 结合关节电流信号和运动模式切换时振动信号的趋势， 确定最终跳变点作为运动模式的起始点。首先， 将获得的跳变点位置转换到原始振动信号中， {\mathit{\tau }}^{″}=\left({\tau }_{\mathrm{1}}^{″},\cdots ,{\tau }_c^{″}\right),{\tau }_i^{″}={\tau }_i^{\text{'}}\cdot l_{\mathrm{w}\mathrm{i}\mathrm{n}}。 然后， 判断电流信号{\mathit{z}}_{\mathrm{1}:n}=(z_{\mathrm{1}},\cdots ,z_n)在跳变点\pm \delta范围内均值是否存在， 并且跳变点处振动信号趋势（上升沿）是否符合

（16）

5） 根据不同运动模式的运行时长和起始跳变点位置， 确定不同运动模式相应的信号区间。获取到筛选后的跳变点集合， 即获得了每次采样下运动模式切换的起始时刻信息， 结合前期获取的不同运动模式的运行时长， 即可进一步划分每一种运动模式对应的振动信号区间， 为

{\mathit{y}}_i^i={\mathit{y}}_{\mathrm{\Omega }}，

\mathrm{\Omega }={\tau }_{\mathrm{1}}^{\mathrm{*}}+{\displaystyle \sum _{j=\mathrm{1}}^i}{t}_{i-\mathrm{1}}\cdot f_{\mathrm{s}}:{\tau }_{\mathrm{1}}^{\mathrm{*}}+{\displaystyle \sum _{j=\mathrm{1}}^i}{t}_i\cdot f_{\mathrm{s}},（17）

式中： {\mathit{y}}_i^i为第i次采样第i个运动模式对应的振动信号区间； \mathrm{\Omega }为第i个运动模式对应的振动信号区间索引； f_{\mathrm{s}}为采样频率。

6） 对上述步骤获取的每个运动模式对应的振动信号区间{\mathit{y}}_i^i， 分别提取监测指标。为获得适用于健康监测的监测指标， 前期进行了工业机器人关节故障模拟预实验， 通过分析不同指标对通道、 负载、 转速以及运动路径的鲁棒性， 最终选取均方根值指标为默认优选指标。

7） 基于控制图法对监测指标进行监测。经过N次采样后， 获得每个运动模式监测指标序列(I\mathrm{d}{x}_{\mathrm{1}}^i,\cdots ,I\mathrm{d}{x}_i^i,\cdots ,I\mathrm{d}{x}_N^i)， 其中， 控制图法的上下监测边界阈值计算公式为

UC{L}^i={\mu }_{I\mathrm{d}{x}_{\mathrm{1}:\alpha }^i}+\mathrm{3}{\sigma }_{I\mathrm{d}{x}_{\mathrm{1}:\alpha }^i},

LC{L}^i={\mu }_{I\mathrm{d}{x}_{\mathrm{1}:\alpha }^i}-\mathrm{3}{\sigma }_{I\mathrm{d}{x}_{\mathrm{1}:\alpha }^i},（18）

式中： I\mathrm{d}{\mathit{x}}_{\mathrm{1}:\alpha }^i为前次采样第i个运动模式的监测指标序列(I\mathrm{d}{x}_{\mathrm{1}}^i,\cdots ,I\mathrm{d}{x}_i^i,\cdots ,I\mathrm{d}{x}_{\alpha }^i)； UC{L}^i为第i个运动模式的监测上边界； LC{L}^i为第i个运动模式的监测下边界。

3 工业机器人退化实验数据集验证

3.1 实验介绍

为了验证所提方法的有效性， 利用某型号机器人退化实验台上采集的数据进行验证， 此型号工业机器人为典型的6关节工业机器人。为模拟实际工业机器人执行特定作业任务时的场景， 设定工业机器人轨迹规划方案为： 各个关节依次往复摆动进而所有关节联动， 1个周期共有7种运动模式， 其中， 运动模式1~6对应关节1~6往复摆动， 运动模式7对应所有关节联动。

此外， 对工业机器人运行过程中不同关节位置处的振动信号和电流信号进行了采集， 累计共18通道信号， 具体各通道详细信息如表 2 所示。

如图 3 所示， 4个三向振动传感器分别布置在底座、 关节 2 和关节 3 之间的臂上、 关节 4 和关节 5 之间的臂上以及负载重物上， 振动传感器型号为亿恒EV4540， 灵敏度为100 mV/g； 电流传感器型号为维博WBI411N95， 灵敏度为1/6 mV/mA； 数采型号为NET0824。为完整地获取多个工作周期的信号， 采样方案设定为： 单次采样时长为3 min（大于2个工业机器人的工作周期）， 采样间隔为30 min， 采样频率为6 250 Hz。

3.2 实验结果分析

选择某次采集的位置3 y方向振动加速度信号（选择依据： 对所有通道进行跳变点算法预实验， 选取识别结果较为鲁棒的通道）和关节 4 电机U相电流信号（选择依据： 靠近位置 3 的关节电机）， 详细展示算法的实施过程。首先， 根据前期轨迹规划算法获得7种运动模式的运行时长， 分别为9.16 s， 5.85 s， 6.27 s， 10.81 s， 9.63 s， 8.28 s， 7.50 s。其次， 对振动信号进行降采样处理。设置长度为100的滑窗在原始振动信号上滑动， 对滑窗内的信号区间求均方根值获得降采样信号， 原始振动信号和降采样信号对比如图 4 所示， 可以看出降采样信号的跳变点效应得到了改善。

在使用跳变点识别算法之前， 需要确定算法的惩罚项系数\beta， 因此， 本文将探究惩罚系数取值与跳变点识别结果评价指标F1得分之间的关系。设置惩罚系数取值区间为［1， 191］， 取值步长为10（即\beta分别取值为1， 11， 21， …， 191）。此外， 在确定运动模式切换时刻时， 需要从跳变点识别结果中进行筛选。为了精确和鲁棒地获得运动模式间的切换时刻， 需要简化筛选条件的设置， 因此， 人为设置目标跳变点集为（1 614， 2 544， 5 220， 6 149， 8 824， 9 754）。在不同惩罚系数取值下， 计算跳变点算法获得的跳变点集与目标跳变点集的F1得分， 结果见图 5（a）。随着惩罚系数取值增大， F1结果取值不断增大， 自惩罚系数取值为101后， F1结果开始收敛， 因此， 将惩罚项系数\beta设置为101。进一步设置跳变点算法的惩罚系数\beta， 分别为21， 51， 101， 将获得的跳变点位置进行可视化， 如图5（b）~5（d）所示。

可以看出随着惩罚系数的增大， 识别的跳变点数量不断减小， 说明惩罚系数具有防止跳变点识别结果过拟合的作用， 可以通过改变其取值大小来关注不同程度的跳变效应。

设置惩罚项系数\beta为101后， 可获得跳变点位置集合， 如图 5（d） 所示。结合靠近位置3的关节4的电流信号， 通过设置筛选条件获得筛选后的跳变点位置， 如图 6 所示， 筛选后的跳变点集即为当前采样样本内运动模式4的起始位置。由于运动模式运行时长已知， 因此， 可以确定当前采样样本内不同运动模式的切换时刻， 如图 7 所示， 进而获取不同运动模式对应的振动信号区间。针对每一种运动模式对应的信号区间提取监测指标， 然后， 计算监测指标的上下边界阈值（取序列总长度的前20%为默认健康状态）， 如图 8 所示。为了减少运算量和便于可视化， 单次采样样本的间隔是20 d。

从图 8 可以发现， 在位置2测点处运动模式2的监测指标超过监测上边界阈值。为进一步排除传感器漂移、 信号幅值正常波动等因素， 将该测点处获取的运动模式1和运动模式2的监测指标和不同时刻采集的振动信号（选取依据： 根据均方根值指标的幅值选择退化的3个阶段， 对应图中的健康状态1， 2， 3）对比情况展示在图 9 中， 而且不同时刻采集的运动模式信号分别用蓝、 绿、 红表示， 并对应标注在监测指标中。通过对比可以发现， 该测点处运动模式1并没有明显的幅值变化， 因此， 可以排除传感器漂移、 信号幅值正常波动等因素， 确定运动模式2监测指标幅值的变化是由健康状态变化引起的。

为了对大量退化数据进行分析， 需要运动模式划分算法具有 3 个特点： 1） 运行时间短； 2） 单次精度高； 3） 重复性好。因此， 本文针对随机选取的10组采样数据， 选用跳变点识别结果评价指标， 对5种对比方法（包括本文方法）进行了评价。

5种方法分别为： 1） 本文方法（M1）， 设置算法惩罚系数为8； 2） 二元分隔（Binary segmentation algorithm， BSA）方法（M2）^［18］， 设定序列变化模式是均值变化， 惩罚项为贝叶斯信息准则（Bayesian information criterion， BIC）； 3） 贝叶斯在线变化点检测方法（M3）^［19］， 窗长和时间尺度\lambda均为1 000； 4） 递归平均差最大化（Recursive mean difference maximization， RMDM）方法（M4）^［20］， 显著性水平设置为0.95； 5）隐逻辑回归方法（M5）， 跳变点个数设置为6， 多项式次数设置为0。为了降低运算耗时， 在进行本次对比实验之前， 设置降采样窗长为1 000， 此外， 人为设定目标跳变点集为关节4与关节5运行区间起始位置处， 如图 5（d） 所示。

最终计算10次评价结果的均值和方差， 同时， 给出了识别算法的平均运行时间， 详见表 3。

结果显示， 本文方法的F1得分最高， 为0.800， 虽然对比方法M5也取得了较为接近的得分， 但是计算时间是本方法的270倍。进一步可视化了不同跳变点识别方法以及提出方法在不同惩罚系数设置、 不同组别样本上的跳变点识别结果， 见图 10。通过对比图 10（a）~10（c）， 可以发现\beta可以有效地限制提出算法获取的跳变点集的个数。通过对比图 10（d）~10（f）， 可以发现对比方法虽然可以较好地识别到目标跳变点， 但是由于识别结果过多， 不易设置筛选条件。通过对比图 10（g）~10（i）， 可以发现本文方法在不同样本上的识别结果基本一致， 也进一步证明了本文方法在10组随机采集的信号跳变点分析中取得了较低的评价结果， 方差为0.004， 如表 4 所示。

4 结 论

本文研究了基于运动模式划分的工业机器人健康监测方法， 针对多运动模式切换带来的问题， 提出了联合关节振动和关节电机电流信号实现工业机器人运动模式划分的方法， 并针对不同运动模式对应的信号区间提取监测指标， 使用控制图法实现了运动模式的健康监测。实验结果表明， 该方法能够有效地避免运动模式耦合和采样信号的差异可能导致的监测结果误判， 使得监测结果更加精确和可靠。本文针对工业机器人的健康监测， 只需要在位置2和位置3布置1个y方向的单向振动传感器， 便可对工业机器人各关节健康状态进行分别监测， 为后续针对核心组件如减速器、 电机等的研究奠定了数据基础。

参考文献

[1] CUI G， LI B， TIAN W， et al. Dynamic modeling and vibration prediction of an industrial robot in manufacturing［J］. Applied Mathematical Modelling， 2022， 105： 114-136.

[2] 马宝瑞. 基于数字孪生的工业机器人关节轴承健康状况评估技术研究［D］. 武汉： 华中科技大学， 2022.

[3] SNEIDER H, FRANK P M. Observer-based supervision and fault detection in robots using nonlinear and fuzzy logic residual evaluation[J]. IEEE Transactions on Control Systems Technology, 1996, 4(3): 274-282.

[4] ENDO Y， OSHIMA M， MIYAZAKI T， et al. A prototype of fault diagnostic system for robots［C］//Proceedings of the 2nd International Conference on Industrial and Engineering Applications of Artificial Intelligence and Expert Systems. 1989： 26-35.

[5] BORGI T， HIDRI A， NEEF B， et al. Data analytics for predictive maintenance of industrial robots［C］//2017 International Conference on Advanced Systems and Electric Technologies （IC_ASET）. IEEE， 2017： 412-417.

[6] CHEN C， LIU C， WANG T， et al. Compound fault diagnosis for industrial robots based on dual-transformer networks［J］. Journal of Manufacturing Systems， 2023， 6（3）： 163-178.

[7] JIANG L L， YOURUI T， JIA W. Remaining useful life prediction for reducer of industrial robots based on MCSA［C］//2021 Global Reliability and Prognostics and Health Management （PHM-Nanjing）. 2021： 1-7.

[8] SAFIZADEH M S， LATIFI S K. Using multi-sensor data fusion for vibration fault diagnosis of rolling element bearings by accelerometer and load cell［J］. Information Fusion， 2014， 18： 1-8.

[9] JABER A A， BICKER R. Wireless fault detection system for an industrial robot based on statistical control chart［J］. International Journal of Electrical and Computer Engineering， 2017， 7（6）： 3421.

[10] VAGAŠ M， SEMJON J， BALAŽ V， et al. Methodology for the vibration measurement and evaluation on the industrial robot KUKA［C］//2014 23rd International Conference on Robotics in Alpe-Adria-Danube Region （RAAD）. IEEE， 2014： 1-6.

[11] BOSEN M， LULU J， JIA W. Degradation modeling for the reducer of industrial robot under periodic multi-stage task［C］//2022 Global Reliability and Prognostics and Health Management （PHM-Yantai）. 2022： 1-6.

[12] KIM Y， PARK J， NA K， et al. Phase-based time domain averaging （PTDA） for fault detection of a gearbox in an industrial robot using vibration signals［J］. Mechanical Systems and Signal Processing， 2020， 138： 106544.

[13] ZHI Z， LIU L， LIU D， et al. Fault detection of the harmonic reducer based on CNN-LSTM with a novel denoising algorithm［J］. IEEE Sensors Journal， 2021， 22（3）： 2572-2581.

[14] SEO J J， YOON H， HA H， et al. Infrared thermographic diagnosis mechnism for fault detection of ball bearing under dynamic loading conditions［J］. Advanced Materials Research， 2011， 295： 1544-1547.

[15] BAGAVATHIAPPAN S， LAHIRI B B， SARAVANAN T， et al. Infrared thermography for condition monitoring—A review［J］. Infrared Physics & Technology， 2013， 60： 35-55.

[16] PENG P， WANG J. Wear particle classification considering particle overlapping［J］. Wear， 2019， 422： 119-127.

[17] CHICCO D， JURMAN G. The advantages of the Matthews correlation coefficient （MCC） over F1 score and accuracy in binary classification evaluation［J］. Bmc Genomics， 2020， 21（6）： 1-13.

[18] CHEN J， GUPTA A K. Parametric statistical change point analysis： with applications to genetics， medicine， and finance［J］. 2012， 2（1）： 212-215.

[19] ADAMS R P， MACKAY D J. Bayesian online changepoint detection［J］. 2007， 12（2）： 34-41.

[20] BERNAOLA-GALVÁN P， IVANOV P C， AMARAL L A N， et al. Scale invariance in the nonstationarity of human heart rate［J］. Physical Review Letters， 2001， 87（16）： 168-172.

基金资助

国家自然科学基金资助项目(52375121)