动力学模态解耦的工业机器人加工颤振规避方法

Chatter Avoidance Method of Industrial Robotic Machining Based on Dynamics Mode Decoupling

摘要

针对机器人加工过程的颤振规避问题，提出了一种动力学模态解耦的工业机器人加工颤振规避方法。建立了某一确定位置下不同姿态的机器人动力学微分方程；分别对再生型颤振和振型耦合型颤振进行稳定性分析，将机器人可达姿态判定为稳定姿态或颤振姿态；分别推导机器人再生型颤振和振型耦合型颤振的稳定性判据，得到不发生再生型颤振和振型耦合型颤振的姿态。依据这些稳定性判据，选择稳定姿态和颤振姿态，开展了机器人锤击实验和机器人加工状态下的颤振规避实验。实验结果验证了所提方法的有效性。

Abstract

To address the issue of chatter avoidance in robotic machining processes， a method was proposed to avoid chatters in industrial robot machining based on dynamics modal decoupling. Firstly， dynamics differential equations were established for the robot at various orientations within a specific location. Secondly， stability analyses were conducted separately for regenerative chatter and mode-coupling chatter， allowing the reachable robot orientations to be classified as either stable or chatter-prone. Thirdly， stability criteria for both regenerative chatter and mode-coupling chatter were derived to identify orientations where neither type of chatter occured. Based on these stability criteria， stable and chatter orientations were selected， and hammering experiments， as well as robotic machining chatter avoidance tests， were carried out.Experimental results validate the effectiveness of the proposed method.

Graphical abstract

关键词

机器人加工颤振 / 颤振规避 / 模态解耦 / 再生型颤振 / 振型耦合型颤振

Key words

robotic machining chatter / chatter avoidance / modal decoupling / regenerative chatter / mode coupling chatter

0 引言

近年来，工业机器人智能化加工逐渐成为高质量制造的新趋势，其中接触式作业工业机器人主要应用于打磨、抛光等复杂曲面制造的关键工序中^［1］。与数控机床相比，工业机器人结构刚度较低，加工过程中易发生颤振现象，导致加工质量变差、加工效率降低，甚至机器人或工件损坏^［2］。

针对工业机器人加工过程中的颤振问题，通常可以在机器人加工之前预测颤振是否发生。叶松涛等^［3］提出一种预测任意姿态下机器人铣削加工的刀尖频响函数（FRF）的方法，实现机器人铣削加工的颤振预测。WU等^［4］提出一种基于机器人模态方向性的振动反馈模态耦合机制来预测机器人铣削中的低频颤振。YANG等^［5］提出一种基于五点高斯求积的铣削颤振预测方法，提高颤振预测的收敛速度。目前，工业机器人加工过程中颤振预测研究主要基于建立的工业机器人动力学模型，然而还存在无法准确建立工业机器人动力学模型和较难准确识别工业机器人系统动力学参数的局限。

为了提高工业机器人加工过程的稳定性，WANG等^［6］建立了一种基于信息熵的变分模式分解支持向量机（VMD-SVM）模型来检测机器人铣削中的颤振。YANG等^［7］提出一种在可变机器人姿态和切削参数的情况下检测机器人铣削中颤振的方法。梁志强等^［8］提出一种基于变分模态滤波、包络滤波和注意力机制网络辨识相结合的重载机器人铣削系统颤振辨识方法。ZHAO等^［9］提出了一种基于磁流变阻尼器的变刚度动力吸振器以实现颤振吸收。TAKAHASHI等^［10］提出一种各向异性结构设计方法，以达到无限动态刚度，实现钻孔作业颤振规避。当前针对工业机器人加工过程中颤振检测、监测、规避和抑制的研究取得了一定效果，但现有方法通常没有建立精确的机器人动力学模型和充分考虑机器人系统动力学参数的耦合作用，并且在颤振分析中通常仅针对某一种颤振类型。在实际机器人加工过程中，再生型颤振和振型耦合型颤振可能同时发生。对于颤振规避，当前研究分析过程中因未充分考虑机器人系统动力学参数耦合的影响、研究颤振类型单一，一定程度上影响了机器人颤振稳定性分析的准确性和适用性。为此，本文提出了一种动力学模态解耦的工业机器人加工颤振规避方法。

1 机器人动力学微分方程建立

1.1 机器人运动学建模

研究对象为一种3T2R构型五自由度工业机器人（图1），机器人机构简图见图2，各参数的含义及取值同文献［11］，其中L_i（i=1，2，…，6）为各杆件几何尺寸。移动关节1和2采用伺服电机M₁和M₂并联驱动与差速驱动的复合方式，实现机器人末端沿水平面内两个方向的平动及绕转台的转动；移动关节3由伺服电机M₃驱动，实现机器人沿竖直方向的平动；转动关节4和5分别由呈45°夹角布置的伺服电机M₄和M₅驱动，实现机器人末端加工姿态的调整。机器人加工姿态调整机构如图1b所示，机器人末端第4、5转动关节两轴线呈45°夹角布置，末端执行器作业刀具的作业点与两个转动关节轴线的交点重合，在4、5两个转动关节转动调整加工姿态的过程中，末端执行器作业刀具的作业点的空间位置始终保持不变，该机器人能够实现末端执行器在上半个完整球面内所有加工姿态可达。

机器人的移动关节1~3为低功率交流伺服电机、联轴器、滚珠丝杠副构成的关节传动链结构，具有较强的弱刚性特征，开展接触式加工作业过程中发生颤振风险较大，严重时甚至造成机器人无法完成作业任务。

根据前期机器人理论研究^［11］建立图3所示的D-H参考坐标系，将图2中机器人结构参数e和刀补L_{\mathrm{p}}（图3中O₅点在Z₄方向坐标值）均视为0，机器人连杆参数见表1。其中，{\phi }_{\mathrm{1}}\sim {\phi }_{\mathrm{5}}分别为J_{\mathrm{1}}~J_{\mathrm{5}}轴转角；\alpha为转台转角；{\theta }_i为关节变量；d_i为关节偏距；a_i为连杆长度。

根据建立的D-H参考坐标系和机器人结构参数，可得到如下相邻坐标系间齐次变换矩阵：

{\mathit{T}}_{\mathrm{1}}^{\mathrm{0}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & -\mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{0}& (s_x{\phi }_{\mathrm{1}}+s_x{\phi }_{\mathrm{2}})/\mathrm{4}\mathrm{\pi }\\ \mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{c}\mathrm{o}\mathrm{s}\alpha & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（1）

式中：s_x为滚珠丝杠{\mathrm{S}}_{\mathrm{1}}、{\mathrm{S}}_{\mathrm{2}}的螺距。

通过机器人D-H参考坐标系O_i相对于O_{i-\mathrm{1}}的齐次变换矩阵{\mathit{A}}_i，可得机器人末端执行器作业刀具坐标系O_{\mathrm{5}}相对于基坐标系O_{\mathrm{0}}的位姿矩阵：

{\mathit{T}}_{\mathrm{5}}^{\mathrm{0}}={\mathit{A}}_{\mathrm{1}}{\mathit{A}}_{\mathrm{2}}{\mathit{A}}_{\mathrm{3}}{\mathit{A}}_{\mathrm{4}}{\mathit{A}}_{\mathrm{5}}=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（2）

式中：[n_x{n}_y{n}_z{]}^{\mathrm{T}}、[o_x{o}_y{o}_z{]}^{\mathrm{T}}、[a_x{a}_y{a}_z{]}^{\mathrm{T}}为机器人姿态矩阵^［11］；[p_x{p}_y{p}_z{]}^{\mathrm{T}}为机器人末端作业点的空间位置。

根据相邻坐标系间齐次变换矩阵计算可得

\left.\begin{array}{l}{p}_x=s_x({\phi }_{\mathrm{1}}+{\phi }_{\mathrm{2}})/(\mathrm{4}\mathrm{\pi })+L_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \\ p_y=L_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \\ p_{\mathrm{z}}=L_{\mathrm{01}}+L_{\mathrm{2}}+L_{\mathrm{4}}+\sqrt[]{\mathrm{2}}{L}_{\mathrm{5}}+s_z{\phi }_{\mathrm{3}}/(\mathrm{2}\mathrm{\pi })\end{array}\right\}（3）

\alpha =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}((X_{\mathrm{2}}-X_{\mathrm{1}})/L_{\mathrm{1}})（4）

式中：s_z为滚珠丝杠{\mathrm{S}}_{\mathrm{3}}的螺距；X_{\mathrm{1}}、X_{\mathrm{2}}分别为螺母在滚珠丝杠{\mathrm{S}}_{\mathrm{1}}、{\mathrm{S}}_{\mathrm{2}}上的行程；转台转角\alpha可达范围为\left[-\mathrm{\pi }/\mathrm{4},\mathrm{\pi }/\mathrm{4}\right]；L_{\mathrm{1}}~L_{\mathrm{5}}均为机器人结构尺寸^［11］；L_{\mathrm{01}}为L_{\mathrm{3}}到丝杠 {\mathrm{S}}_{\mathrm{1}}在竖直方向的距离。

1.2 雅可比矩阵求解

雅可比矩阵是描述关节线速度和角速度与相应末端执行器线速度和角速度之间关系的映射矩阵，它取决于机器人自身的结构与位姿。雅可比矩阵的求解表达式为

\left.\begin{array}{l}{\mathit{J}}_{Pi}=\left[\begin{array}{c}{\mathit{z}}_{i-\mathrm{1}}\\ \mathrm{0}\end{array}\right]\\ {\mathit{J}}_{Oi}=\left[\begin{array}{c}{\mathit{z}}_{i-\mathrm{1}}\times ({\mathit{p}}_{\mathrm{e}}-{\mathit{p}}_{i-\mathrm{1}})\\ {\mathit{z}}_{i-\mathrm{1}}\end{array}\right]\end{array}\right\}（5）

{\mathit{z}}_{i-\mathrm{1}}={\mathit{R}}_{i-\mathrm{1}}^{\mathrm{0}}{\mathit{z}}_{\mathrm{0}}（6）

i=\mathrm{1,2},\cdots ,\mathrm{5}

式中： J_Pi 、 J_Oi 分别为移动关节和转动关节的雅可比矩阵的第i列；{\mathit{p}}_{\mathrm{e}}为机器人末端执行器位置向量，通过变换矩阵{\mathit{T}}_{\mathrm{e}}^{\mathrm{0}}的第4列的前3个元素给出；{\mathit{p}}_{i-\mathrm{1}}由变换矩阵{\mathit{T}}_{i-\mathrm{1}}^{\mathrm{0}}第4列的前3个元素构成；{\mathit{z}}_{i-\mathrm{1}}为关节变量的函数；{\mathit{R}}_{i-\mathrm{1}}^{\mathrm{0}}为坐标系{O_i}_{-\mathrm{1}}相对于O_{\mathrm{0}}的旋转矩阵。

根据图2所示机器人机构简图和图3所示D-H参考坐标系可知，关节1、关节2并联驱动连杆1沿X₀方向平动且绕Z₀方向转动，不属于单纯的平动或转动，存在耦合关系，可依据此运动的物理意义对平动和转动的耦合关系进行等效处理，将其解耦表示，在此将关节1等效处理为移动关节、关节2等效处理为转动关节。根据式（5）求解{\mathit{J}}_{\mathrm{1}}\sim {\mathit{J}}_{\mathrm{5}}列，可得到该机器人雅可比矩阵：

\mathit{J}=\left[\begin{array}{cccc}\begin{array}{cc}{\mathit{J}}_{\mathrm{1}}& {\mathit{J}}_{\mathrm{2}}\end{array}& {\mathit{J}}_{\mathrm{3}}& {\mathit{J}}_{\mathrm{4}}& {\mathit{J}}_{\mathrm{5}}\end{array}\right]（7）

1.3 动力学微分方程建立

在机器人加工作业中，易发生一种比较强烈的振动即颤振，属于自激振动的一种。TLUSTY等^［12］和MERRITT^［13］分别指出再生型颤振和振型耦合型颤振是引起加工颤振的主要因素。建立机器人动力学微分方程：

\mathit{M}\ddot{\mathit{x}}+\mathit{C}\dot{\mathit{x}}+\mathit{K}\mathit{x}=\mathit{F}（8）

式中： M 、 K 、 C 分别为质量矩阵、刚度矩阵和阻尼矩阵； x 为位移向量； F 为外力向量。

1）计算质量矩阵。根据机器人的质量特征和几何特征，建立拉格朗日函数：

T_i={\displaystyle \sum _{j=\mathrm{1}}^n}{D}_{ij}{\ddot{\theta }}_{ij}+{\displaystyle \sum _{j=\mathrm{1}}^n}{\displaystyle \sum _{k=\mathrm{1}}^n}{D}_{ijk}{\dot{\theta }}_j{\dot{\theta }}_k+D_i（9）

i=\mathrm{1,2},\dots ,\mathrm{5}

式中：T_i为机器人加工中关节i所有外力矩之和；{D_i}_j为关节i和j处的耦合惯量，当i=j时表示关节i处的有效惯量；n为机器人关节个数；带有{\dot{\theta }}_j{\dot{\theta }}_k的项表示科里奥利加速度，乘以相应的惯量即为科里奥利力；D_i为关节i处的重力。

由式（9）可进一步得出质量矩阵：

{\mathit{M}}_{\mathit{\theta }}=\left[\begin{array}{cccc}\begin{array}{cc}{D}_{\mathrm{11}}& D_{\mathrm{12}}\end{array}& D_{\mathrm{13}}& D_{\mathrm{14}}& D_{\mathrm{15}}\\ \begin{array}{cc}{D}_{\mathrm{21}}& D_{\mathrm{22}}\end{array}& D_{\mathrm{23}}& D_{\mathrm{24}}& D_{\mathrm{25}}\\ \begin{array}{cc}{D}_{\mathrm{31}}& D_{\mathrm{32}}\end{array}& D_{\mathrm{33}}& D_{\mathrm{34}}& D_{\mathrm{35}}\\ \begin{array}{cc}\begin{array}{c}{D}_{\mathrm{41}}\\ D_{\mathrm{51}}\end{array}& \begin{array}{c}{D}_{\mathrm{42}}\\ D_{\mathrm{52}}\end{array}\end{array}& \begin{array}{c}{D}_{\mathrm{43}}\\ D_{\mathrm{53}}\end{array}& \begin{array}{c}{D}_{\mathrm{44}}\\ D_{\mathrm{54}}\end{array}& \begin{array}{c}{D}_{\mathrm{45}}\\ D_{\mathrm{55}}\end{array}\end{array}\right]（10）

D_{ij}={\displaystyle \sum _{p=max(i,j)}^n}trace({\mathit{U}}_{pj}{\mathit{J}}_p{\mathit{U}}_{pi}^T)（11）

{\mathit{U}}_{ij}=\frac{\partial {\mathit{T}}_i^{\mathrm{0}}}{\partial {\theta }_j}=\frac{\partial ({\mathit{A}}_{\mathrm{1}}{\mathit{A}}_{\mathrm{2}}\cdots {\mathit{A}}_j\cdots {\mathit{A}}_i)}{\partial {\theta }_j}（12）

j\le i

式中：{\mathit{U}}_{pj}为速度传递矩阵；{\mathit{J}}_p为连杆惯性张量；trace（·）用于求矩阵的迹。

由式（11）、式（12）可知，D_{ij}的值与齐次变换矩阵 A_i 和关节变量θ_j 有关。于是，改变机器人关节变量可以得到对应机器人位姿的质量矩阵。

2）计算刚度矩阵。根据雅可比矩阵和胡克定律，建立机器人静刚度模型：

\mathit{F}={\mathit{J}}^{-\mathrm{T}}{\mathit{K}}_{\theta }{\mathit{J}}^{-\mathrm{1}}\mathit{X}（13）

式中： {\mathit{K}}_{\theta }为机器人关节刚度矩阵；\mathit{X}为机器人末端在外力\mathit{F}下的形变。

由式（13）可得

\left[\begin{array}{c}x\\ y\\ z\\ {\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]=\left[\begin{array}{ccc}{J}_{\mathrm{11}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{i\mathrm{1}}{F}_i& \dots & J_{\mathrm{15}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{i\mathrm{5}}{F}_i\\ ⋮& & ⋮\\ J_{\mathrm{61}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{i\mathrm{1}}{F}_i& \cdots & J_{\mathrm{65}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}}}{J}_{i\mathrm{5}}{F}_i\end{array}\right]\left[\begin{array}{c}{K}_{\theta \mathrm{1}}^{-\mathrm{1}}\\ ⋮\\ K_{\theta \mathrm{5}}^{-\mathrm{1}}\end{array}\right]（14）

式中：x、y、z、{\omega }_x、{\omega }_y和{\omega }_z为 X 的分量。

设计机器人关节刚度辨识实验，如图4所示。测量不同位姿下的机器人末端形变与外力，根据式（14）和最小二乘法可求解出该机器人关节刚度矩阵^［14］：

{\mathit{K}}_{\theta }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(k_{{\theta }_{\mathrm{1}}},k_{{\theta }_{\mathrm{2}}},k_{{\theta }_{\mathrm{3}}},k_{{\theta }_{\mathrm{4}}},k_{{\theta }_{\mathrm{5}}})（15）

3）计算阻尼矩阵。机器人的关节阻尼系数c_{\theta i}可以通过关节的阻尼比{\xi }_i和临界阻尼系数{\overline{c}}_i得到，其中阻尼比可以通过对数衰减法计算：

\left.\begin{array}{l}{c_{\theta }}_i={\xi }_i{\overline{c}}_i\\ {\overline{c}}_i=\mathrm{4}\mathrm{\pi }{m}_i{\omega }_{\mathrm{d}}\end{array}\right\}　　 i=\mathrm{1,2},\cdots ,\mathrm{5}（16）

式中：m_i 为机器人连杆i的质量；ω_d为连杆i的自然频率。

由式（16）计算关节阻尼系数c_{{\theta }_{\mathrm{1}}}\sim c_{{\theta }_{\mathrm{5}}}，进而得到由关节阻尼系数构成的阻尼矩阵：

{\mathit{C}}_{\theta }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(c_{{\theta }_{\mathrm{1}}},c_{{\theta }_{\mathrm{2}}},c_{{\theta }_{\mathrm{3}}},c_{{\theta }_{\mathrm{4}}},{c_{\theta }}_{{}_{\mathrm{5}}})（17）

上述质量矩阵、刚度矩阵和阻尼矩阵都是在机器人关节空间中求得，在关节空间中，每个关节是相互独立的，上述矩阵均为5×5对角矩阵。

由上述可知，通过拉格朗日动力学方程、关节刚度辨识实验和对数衰减模型所开展的工业机器人系统动力学参数质量矩阵、刚度矩阵和阻尼矩阵的推导计算和实验计算，可以获取较为精确的工业机器人系统动力学参数，可提高所建立的动力学微分方程的准确性。

2 机器人动力学模态解耦

2.1 末端笛卡儿空间机器人系统动力学参数矩阵

通过雅可比矩阵 J，关节空间转换至末端笛卡儿空间的机器人动力学微分方程为

{\mathit{M}}_{\mathrm{D}}\ddot{\mathit{x}}+{\mathit{C}}_{\mathrm{D}}\dot{\mathit{x}}+{\mathit{K}}_{\mathrm{D}}\mathit{x}={\mathit{F}}_{\mathrm{D}}（18）

\left.\begin{array}{l}{\mathit{M}}_{\mathrm{D}}={\mathit{J}}^{-\mathrm{T}}{\mathit{M}}_{\theta }{\mathit{J}}^{-\mathrm{1}}\\ {\mathit{K}}_{\mathrm{D}}={\mathit{J}}^{-\mathrm{T}}{\mathit{K}}_{\theta }{\mathit{J}}^{-\mathrm{1}}\\ {\mathit{C}}_{\mathrm{D}}={\mathit{J}}^{-\mathrm{T}}{\mathit{C}}_{\theta }{\mathit{J}}^{-\mathrm{1}}\end{array}\right\}（19）

式中：{\mathit{M}}_{\mathrm{D}}、{\mathit{K}}_{\mathrm{D}}、{\mathit{C}}_{\mathrm{D}}和{\mathit{F}}_{\mathrm{D}}分别为末端笛卡儿空间的质量矩阵、刚度矩阵、阻尼矩阵和外力向量；{\mathit{J}}^{-\mathrm{1}}为雅可比矩阵的逆。

由式（7）可知，雅可比矩阵取决于机器人自身的结构和位姿，末端笛卡儿空间的质量矩阵{\mathit{M}}_{\mathrm{D}}、刚度矩阵{\mathit{K}}_{\mathrm{D}}和阻尼矩阵{\mathit{C}}_{\mathrm{D}}均与机器人本身的结构和位姿有关，并且均为对称矩阵，但并非对角阵。由式（18）可知，由于机器人动力学微分方程中{\mathit{M}}_{\mathrm{D}}、{\mathit{K}}_{\mathrm{D}}、{\mathit{C}}_{\mathrm{D}}均为非对角阵，{\mathit{F}}_{\mathrm{D}}每个分量对应自由度方向上的机器人动力学微分方程会受非对角线元素的耦合作用，无法建立每个自由度方向上独立的机器人动力学微分方程，若忽略非对角线元素则会造成机器人动力学模型信息的缺失，降低机器人动力学微分方程的准确性。为充分考虑机器人系统动力学参数的耦合作用，需要对末端笛卡儿空间的质量矩阵{\mathit{M}}_{\mathrm{D}}、刚度矩阵{\mathit{K}}_{\mathrm{D}}和阻尼矩阵{\mathit{C}}_{\mathrm{D}}进行解耦。

2.2 末端笛卡儿空间机器人系统动力学参数矩阵模态解耦

本文研究中，假设影响加工质量的主要因素是作业刀具的线位移形变^［15］，不考虑加工过程中力矩引起的作业刀具扭转形变，只考虑由外力引起的线位移。于是，对{\mathit{M}}_{\mathrm{D}}、{\mathit{K}}_{\mathrm{D}}、{\mathit{C}}_{\mathrm{D}}矩阵进行划分，只考虑第1~3行和第1~3列构成的维度为\mathrm{3}\times \mathrm{3}的质量-线位移子矩阵{\mathit{M}}_x、刚度-线位移子矩阵{\mathit{K}}_x、阻尼-线位移子矩阵{\mathit{C}}_x。

由于末端笛卡儿空间中的质量-线位移子矩阵、刚度-线位移子矩阵和阻尼-线位移子矩阵均为非对角阵，非对角线元素即为耦合项，为了解耦质量-线位移子矩阵中的惯性耦合项和刚度-线位移子矩阵中出现的弹性耦合项，将质量-线位移子矩阵、刚度-线位移子矩阵和阻尼-线位移子矩阵由末端笛卡儿坐标系O_{\mathrm{D}}{X}_{\mathrm{D}}{Y}_{\mathrm{D}}{Z}_{\mathrm{D}}转换至解耦坐标系O_{\mathrm{P}}{X}_{\mathrm{P}}{Y}_{\mathrm{P}}{Z}_{\mathrm{P}}中，如图5所示。质量-线位移子矩阵、刚度-线位移子矩阵转换至解耦坐标系后均为对角阵，不再存在耦合项，阻尼-线位移子矩阵为非对角阵，但通常采用忽略非对角元素的处理方式^［16］，于是，将原本相互耦合的三自由度问题，解耦转化为3个单自由度问题。

同时，将质量-线位移子矩阵和刚度-线位移子矩阵转换至解耦坐标系（即同时将该两个子矩阵对角化处理），需要计算对应模态矩阵，通过求解具有特征值{\omega }_i^{\mathrm{2}}的特征多项式，并由特征向量{\mathit{\varphi }}_i（也称主振型）建立模态矩阵\mathit{\Phi }：

({\mathit{K}}_x-{{\omega }_i}^{\mathrm{2}}{\mathit{M}}_x)\mathit{\Phi }=\mathrm{0}（20）

\mathit{\Phi }=\left[\begin{array}{ccc}{\mathit{\varphi }}_{\mathrm{1}}& {\mathit{\varphi }}_{\mathrm{2}}& {\mathit{\varphi }}_{\mathrm{3}}\end{array}\right]（21）

式中：第i个特征值{\omega }_i^{\mathrm{2}}的算术平方根{\omega }_i^{}即为第i阶固有频率。

根据所求解的模态矩阵，解耦末端笛卡儿空间的质量-线位移子矩阵和刚度-线位移子矩阵：

\left.\begin{array}{l}{\mathit{M}}_{\mathrm{P}}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{M}}_x\mathit{\Phi }\\ {\mathit{K}}_{\mathrm{P}}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{K}}_x\mathit{\Phi }\\ {\mathit{C}}_{\mathrm{P}}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{C}}_x\mathit{\Phi }\end{array}\right\}（22）

{\mathit{M}}_{\mathrm{P}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{M}}_{\mathrm{P}x},{\mathit{M}}_{\mathrm{P}y},{\mathit{M}}_{\mathrm{P}z})

{\mathit{K}}_{\mathrm{P}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{K}}_{\mathrm{P}x},{\mathit{K}}_{\mathrm{P}y},{\mathit{K}}_{\mathrm{P}z})

式中：{\mathit{M}}_{\mathrm{P}}、{\mathit{K}}_{\mathrm{P}}分别为解耦后在解耦坐标系的质量-线位移子矩阵和刚度-位移子矩阵；{\mathit{C}}_{\mathrm{P}}一般为非对角阵，忽略{\mathit{C}}_{\mathrm{P}}中的非对角元素，即{\mathit{C}}_{\mathrm{P}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{C}}_{\mathrm{P}\mathrm{x}},{\mathit{C}}_{\mathrm{P}\mathrm{y}},{\mathit{C}}_{\mathrm{P}\mathrm{z}})。

解耦后的机器人系统动力学微分方程为

{\mathit{M}}_{\mathrm{P}}\ddot{\mathit{x}}+{\mathit{C}}_{\mathrm{P}}\dot{\mathit{x}}+{\mathit{K}}_{\mathrm{P}}\mathit{x}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{F}}_{\mathrm{D}}\mathit{\Phi }（23）

由{\mathit{M}}_{\mathrm{P}}\mathrm{、}{\mathit{K}}_{\mathrm{P}}\mathrm{、}{\mathit{C}}_{\mathrm{P}}为对角矩阵可知，在质量矩阵中，不存在惯性耦合项，某个自由度所对应运动方向上产生的加速度只在该方向上引起惯性力；在刚度矩阵中，不存在弹性耦合项，某个自由度所对应运动方向上产生的位移只在该方向上引起弹性恢复力。于是，将一个原来相互耦合的三自由度机器人颤振问题，经解耦后转化为3个单自由度系统的机器人颤振问题，使得解耦后3个自由度所对应运动方向上的动力学微分方程稳定性分析互不影响，有效降低了机器人颤振问题分析的复杂程度。

由式（23）可知，机器人动力学微分方程会随着{\mathit{M}}_{\mathrm{P}}\mathrm{、}{\mathit{K}}_{\mathrm{P}}\mathrm{、}{\mathit{C}}_{\mathrm{P}}的改变而改变。通过改变机器人位姿进而改变机器人动力学微分方程，便可建立机器人所有可达位姿与机器人动力学微分方程之间的映射关系。

2.3 解耦坐标系姿态确定

由于机器人系统的3个主振型{\mathit{\varphi }}_i(i=\mathrm{1,2},\mathrm{3})关于质量矩阵及刚度矩阵相互正交，则各个主振型必然是线性无关的，即任意一个3维位移向量\mathit{x}都能被唯一表示成主振型的线性组合^［16］，将其写成矩阵形式：

\mathit{x}=\mathit{\Phi }\mathit{\eta }（24）

式中：\mathit{\eta }=\left[\begin{array}{ccc}{\eta }_{\mathrm{1}}& {\eta }_{\mathrm{2}}& {\eta }_{\mathrm{3}}\end{array}\right]为末端笛卡儿坐标系中向量\mathit{x}经过模态矩阵\mathit{\Phi }转换后，在解耦坐标系的新坐标向量。

由上述可知，模态矩阵\mathit{\Phi }可以实现末端笛卡儿坐标系中的位移向量与解耦坐标系的新坐标向量的转换，模态矩阵即为解耦坐标系O_{\mathrm{P}}{X}_{\mathrm{P}}{Y}_{\mathrm{P}}{Z}_{\mathrm{P}}相对于机器人加工系统末端笛卡儿坐标系O_{\mathrm{D}}{X}_{\mathrm{D}}{Y}_{\mathrm{D}}{Z}_{\mathrm{D}}的旋转矩阵{\mathit{R}}_{\mathrm{P}}^{\mathrm{D}}。由图3和图5可知，机器人刀具坐标系O₅经旋转矩阵{\mathit{R}}_{\mathrm{D}}^{\mathrm{5}}可转换至加工系统末端笛卡儿坐标系，并根据作业刀具位姿矩阵{\mathit{T}}_{\mathrm{5}}^{\mathrm{0}}计算出加工系统末端笛卡儿坐标系相对于基坐标系O₀的旋转矩阵{\mathit{R}}_{\mathrm{D}}^{\mathrm{0}}，进而计算出解耦坐标系相对于基坐标系的旋转矩阵{\mathit{R}}_{\mathrm{P}}^{\mathrm{0}}：

\left.\begin{array}{l}{\mathit{R}}_{\mathrm{D}}^{\mathrm{0}}={\mathit{R}}_{\mathrm{5}}^{\mathrm{0}}{\mathit{R}}_{\mathrm{D}}^{\mathrm{5}}\\ {\mathit{R}}_{\mathrm{P}}^{\mathrm{0}}={\mathit{R}}_{\mathrm{5}}^{\mathrm{0}}{\mathit{R}}_{\mathrm{D}}^{\mathrm{5}}{\mathit{R}}_{\mathrm{P}}^{\mathrm{D}}\end{array}\right\}（25）

式中：{\mathit{R}}_{\mathrm{5}}^{\mathrm{0}}为机器人刀具坐标系相对于基坐标系的旋转矩阵。

根据式（25）所示的{\mathit{R}}_{\mathrm{P}}^{\mathrm{0}}可计算出解耦坐标系在基坐标系中的具体姿态，得到确定解耦后3个单自由度系统机器人颤振问题在机器人基坐标系的对应方向，从而确定机器人颤振问题的主方向，可以针对性地对主方向的颤振进行规避。

3 机器人颤振稳定性分析

机器人加工系统的自激特性导致的颤振现象是目前机器人在接触式加工应用中的主要挑战。机器人加工过程中，因前后两次相邻加工轨迹的加工厚度和接触力变化，易产生再生型颤振^［17］，将影响机器人作业精度和质量。多自由度下发生振动且存在彼此耦合时，机器人易产生振型耦合型颤振^［18］，当其发生时，整个机器人结构将剧烈振动。这两类机器人颤振均会影响机器人加工作业稳定性以及加工精度和质量，甚至造成加工作业失效且损坏机器人系统和工件^［19］。

3.1 再生型颤振稳定性分析

在机器人加工作业过程中，加工厚度和切削力因前一加工轨迹留下的波纹和当前加工轨迹产生的波纹之间的相位差而变化。如果前一次加工留下的波纹和当前加工的波纹是同相位的，则动态加工厚度和动态接触力将是恒定的；如果相位不同，则会导致动态加工厚度和动态接触力发生变化，可能导致再生型颤振的发生^［20］。根据加工作业时的工况和受力情况，在解耦后的坐标系中建立机器人加工再生型颤振系统的动力学模型，如图6所示。其中，为简洁表示，使用部分作业刀具代表整个机器人加工系统。

由图6所示机器人加工再生型颤振系统动力学模型和式（23）可得，所研究机器人再生型颤振的动力学微分方程为

{\mathit{M}}_{\mathrm{P}}\ddot{\mathit{x}}(t)+{\mathit{C}}_{\mathrm{P}}\dot{\mathit{x}}(t)+{\mathit{K}}_{\mathrm{P}}\mathit{x}(t)=\mathrm{\Delta }\mathit{F}(t)\mathit{V}（26）

\left.\begin{array}{l}\mathrm{\Delta }\mathit{F}\left(t\right)={\mathit{\Phi }}^{\mathrm{T}}{\mathit{K}}_{\mathrm{c}}\mathit{\Phi }\mathit{b}(y(t-T)-y(t))\\ \mathit{V}={\mathit{R}}_{\mathrm{D}}^{\mathrm{5}}{\mathit{R}}_{\mathrm{P}}^{\mathrm{D}}{\mathit{V}}_{\mathrm{a}}\end{array}\right\}（27）

式中：{\mathit{V}}_{\mathrm{a}}为加速度传感器所获得加速度向量与坐标系O_{\mathrm{5}}{X}_{\mathrm{5}}{Y}_{\mathrm{5}}{Z}_{\mathrm{5}}之间的方向余弦\mathrm{c}\mathrm{o}\mathrm{s}{{\gamma }_{\mathrm{a}}}_x\mathrm{、}\mathrm{c}\mathrm{o}\mathrm{s}{{\gamma }_{\mathrm{a}}}_y\mathrm{、}\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{a}z}构成的对角矩阵，{\mathit{V}}_{\mathrm{a}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{c}\mathrm{o}\mathrm{s}{{\gamma }_{\mathrm{a}}}_x,\mathrm{c}\mathrm{o}\mathrm{s}{{\gamma }_{\mathrm{a}}}_y,\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{\mathrm{a}z})； V 为解耦坐标系中刀具接触动态反作用力\mathrm{\Delta }\mathit{F}(t)与3个自由度振动方向的方向余弦构成的对角矩阵， V =diag（\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_x，\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_y，\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_z）；\mathit{x}(t)为作业刀具振动位移向量；\mathit{b}为加工宽度向量；y(t)、y(t-T)为加工产生的波纹面，其中T为加工周期；{\mathit{K}}_{\mathrm{c}}为单位加工宽度的刚度系数矩阵，解耦坐标系中的单位加工宽度的刚度系数矩阵{\mathit{K}}_{\mathrm{c}\mathrm{P}}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{K}}_{\mathrm{c}}\mathit{\Phi }，为了将解耦坐标系中3个自由度所对应运动方向的微分方程独立求解，忽略{\mathit{K}}_{\mathrm{c}\mathrm{P}}中的非对角线元素。

每个加工姿态都有相对应的质量矩阵、刚度矩阵和阻尼矩阵，于是每个加工姿态都有与其相对应的动力学微分方程，在变换加工姿态时，动力学微分方程也会随之改变。在分析每一个加工姿态下的再生型颤振时，将该姿态相对应的动力学微分方程变换至解耦坐标系中，便可将3个自由度所对应运动方向的微分方程独立求解。

将式（27）代入式（26），根据微分方程的稳定性可知，设系统的特征根为S=\sigma +\mathrm{i}w，若系统特征根的实部均为负值，即当时，系统是稳定的，则可以判定此时机器人系统不发生再生型颤振。反之，若特征根中有一个或多个根具有正实部，即存在\sigma >\mathrm{0}，此时系统不稳定，则可以判定此时机器人系统发生再生型颤振。因此，当\sigma =\mathrm{0}时，可认为系统处于临界稳定状态。此时，临界稳定状态下对应的极限加工宽度为

({\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}{)}_{\mathrm{m}\mathrm{i}\mathrm{n}}=\frac{{\mathit{C}}_{\mathrm{P}}(\mathrm{2}\sqrt[]{{\mathit{K}}_{\mathrm{P}}{\mathit{M}}_{\mathrm{P}}}+{\mathit{C}}_{\mathrm{P}})}{\mathrm{2}{\mathit{K}}_{\mathrm{c}\mathrm{P}}{\mathit{M}}_{\mathrm{P}}\mathit{V}}（28）

式中：{\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}为极限加工宽度矩阵，{\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(b_x,b_y,b_z)，b_x、b_y、b_z为3个自由度所对应运动方向的极限加工宽度。

临界稳定状态下对应的极限加工宽度与加工周期T无关，而主轴转速n=\mathrm{60}/T，因此取任一加工周期和对应主轴转速，所求({\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}{)}_{\mathrm{m}\mathrm{i}\mathrm{n}}一致。

由式（28）可知，{\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}的值与解耦坐标系中的质量矩阵、阻尼矩阵、刚度矩阵、刚度系数向量、解耦坐标系中刀具接触动态反作用力\mathrm{\Delta }\mathit{F}(t)与3个自由度振动方向的方向余弦构成的对角矩阵\mathit{V}有关。因此，可以求出每个加工姿态对应的极限加工宽度，当该姿态下实际作业加工宽度小于极限加工宽度时，可以认为此时机器人系统不发生再生型颤振。在进行实际加工作业之前，可根据所求极限加工宽度，在小于所求极限加工宽度的基础上，综合考虑被加工工件材料特性、加工刀具特性和加工要求（如精度、表面质量和效率）等因素，选择适宜的加工宽度，从而规避机器人加工作业过程中的再生型颤振。

3.2 振型耦合型颤振稳定性分析

在机器人加工作业过程中，当同时考虑两个甚至多个方向的颤振时，根据加工作业时的工况和受力情况，在解耦后的坐标系中建立机器人加工振型耦合型颤振系统的动力学模型，如图7所示。其中，为简洁表示，使用部分作业刀具代表整个机器人加工系统。通常直接计算机器人系统动力学参数含有耦合项的耦合型颤振的稳定性判据较为复杂且难以得到解析解，通过前述研究可将机器人加工系统振型耦合型颤振解耦为不同自由度所对应运动方向的振型耦合型颤振；将非线性耦合的三自由度问题解耦，转化为3个单自由度的线性问题，能够单独规避每个自由度方向的振型耦合型颤振，降低计算复杂度，并可将机器人系统动力学微分方程分解成3个自由度方向上的机器人系统动力学微分方程，来分析振型耦合型颤振，表述形式更为简洁。

根据机器人加工振型耦合型颤振系统的动力学模型（图7），同时考虑解耦坐标系中3个自由度所对应运动方向的振动，建立机器人系统动力学微分方程：

{\mathit{M}}_{\mathrm{P}}\ddot{\mathit{x}}+{\mathit{C}}_{\mathrm{P}}\dot{\mathit{x}}+{\mathit{K}}_{\mathrm{P}}\mathit{x}={\mathit{\Phi }}^{\mathrm{T}}{\mathit{F}}_{\mathrm{D}}\mathit{\Phi }（29）

因阻尼的存在会增加机器人加工系统的稳定性，于是，在振型耦合型颤振分析中，取阻尼系数为0，忽略阻尼的影响，采用该种稳定性较差的情况进行分析。同时，根据“power force”理论^［21］，为假设作用在作业刀具的动态反作用力向量{\mathit{F}}_{\mathrm{D}}正比于位移向量\mathit{x}，设正比系数矩阵为\mathit{Q}，可进一步得到机器人加工系统动力学微分方程：

{\mathit{M}}_{\mathrm{P}}\ddot{\mathit{x}}+{\mathit{K}}_{\mathrm{P}}\mathit{x}={\mathit{\Phi }}^{\mathrm{T}}\mathit{Q}\mathit{\Phi }\mathit{x}（30）

令{\mathit{Q}}_{\mathrm{P}}={\mathit{\Phi }}^{\mathrm{T}}\mathit{Q}\mathit{\Phi }，且{\mathit{M}}_{\mathrm{P}}为非零对角阵，可得

\ddot{\mathit{x}}={\mathit{M}}_{\mathrm{P}}^{-\mathrm{1}}({\mathit{Q}}_{\mathrm{P}}-{\mathit{K}}_{\mathrm{P}})\mathit{x}（31）

由式（31）可知，因需对{\mathit{M}}_{\mathrm{P}}矩阵求逆，对应求逆的时间复杂度为O(n)。若式（31）采用末端笛卡儿空间的动力学参数表达，即\ddot{\mathit{x}}={\mathit{M}}_{\mathrm{D}}^{-\mathrm{1}}({\mathit{Q}}_{\mathrm{D}}-{\mathit{K}}_{\mathrm{D}})\mathit{x}，此时{\mathit{M}}_{\mathrm{D}}为\mathrm{3}\times \mathrm{3}非对角矩阵，对应求逆的时间复杂度为O(n^{\mathrm{3}})。因此，动力学模态解耦后可有效降低计算复杂度。

由式（29）可知，每个加工姿态都有对应的动力学微分方程，在变换加工姿态时，动力学微分方程也会随之改变。由Routh-Hurwitz稳定判据可知，若线性系统特征方程的根均有负的实部，则系统是稳定的，故而系统稳定时可以认为系统没有发生振型耦合型颤振。反之，若特征方程的根中有一个或多个根具有正实部，则系统是不稳定的，此时可认为系统发生了振型耦合型颤振。进一步，依据式（31）可求出每个加工姿态下对应的线性系统特征方程的根，进行实际加工作业之前，可根据所求的线性系统特征方程的根，选择不发生振型耦合型颤振的加工姿态（即系统特征方程的根均有负的实部，所对应的加工姿态），规避加工过程中的振型耦合型颤振。

3.3 加工姿态颤振规避

由式（28）可知，再生型颤振的临界稳定状态下，极限加工宽度与{\mathit{M}}_{\mathrm{P}}\mathrm{、}{\mathit{K}}_{\mathrm{P}}\mathrm{、}{\mathit{C}}_{\mathrm{P}}\mathrm{、}{\mathit{K}}_{\mathrm{c}\mathrm{P}}\mathrm{、}\mathit{V}有关。于是，在某一确定位置下，改变机器人所有可达加工姿态，即可调整机器人3个自由度所对应运动方向上的极限加工宽度，进而可从中选择机器人系统不会发生再生型颤振对应的加工姿态，从而实现机器人加工作业再生型颤振规避。

与此同时，由式（31）可知，振型耦合型颤振的动力学方程是否稳定与{\mathit{M}}_{\mathrm{P}}\mathrm{、}{\mathit{Q}}_{\mathrm{P}}\mathrm{、}{\mathit{K}}_{\mathrm{P}}有关。同样，在某一确定位置下，改变机器人加工姿态，便可得到所有可达加工姿态对应机器人动力学特征方程的根，进而根据稳定性判据选择机器人系统不会发生振型耦合型颤振对应加工姿态，从而实现机器人加工作业振型耦合型颤振规避。

由上所述，提出的动力学模态解耦的工业机器人加工颤振规避方法，既能够分别推导出再生型颤振和振型耦合型颤振的稳定性判据，又能够推导出同时分析两类颤振的稳定性判据。当只考虑再生型颤振或振型耦合型颤振时，可以通过选择满足其对应稳定性判据的加工姿态实现再生型颤振或振型耦合型颤振规避。

所提动力学模态解耦的工业机器人加工颤振规避方法可以分别或同时规避再生型颤振和振型耦合型颤振，将两类颤振规避问题转化为机器人加工姿态选择的求解问题，降低了颤振规避问题分析的复杂度。

4 机器人颤振算例稳定性分析及机器人颤振规避实验验证

4.1 机器人系统动力学参数计算

为了验证所提动力学模态解耦的工业机器人加工颤振规避方法的有效性，需要在不改变机器人末端执行器空间位置的前提下，变换机器人加工姿态开展颤振规避研究。由机器人结构（图2）可知，末端第4、5转动关节两轴线呈45°夹角布置，当作业刀具的作业点与第4、5转动关节轴线交点重合时，仅第4、5转动关节调整加工姿态，即可保持末端执行器作业刀具的作业点空间位置始终不变。于是，保持关节1、2、3的转角{\phi }_{\mathrm{1}}\mathrm{、}{\phi }_{\mathrm{2}}\mathrm{、}{\phi }_{\mathrm{3}}不变，只改变关节4、5的转角{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}，可以保持机器人末端作业点的空间位置不变，对应机器人加工姿态可实现上半个球域范围内灵巧可达。因此，机器人上半球域加工姿态灵巧可达的特性为颤振规避提供了良好的结构基础和条件。在机器人加工作业过程中，可根据实际加工工况，在上半球域内可达姿态中，筛选不发生颤振的加工姿态。

由上文可知，在解耦坐标系O_{\mathrm{P}}{X}_{\mathrm{P}}{Y}_{\mathrm{P}}{Z}_{\mathrm{P}}中，机器人系统动力学微分方程中的参数矩阵会随着姿态的改变而改变。因此，计算机器人所有可达姿态相对应机器人系统动力学微分方程的稳定性判据表达式，便可根据具体的加工工况及参数筛选出机器人加工系统稳定状态下的姿态。仿真分析中的姿态序列及其对应的{\phi }_{\mathrm{4}}、{\phi }_{\mathrm{5}}取值见表2。

根据表2中机器人姿态序列和3T2R构型机器人结构参数，计算该机器人的雅可比矩阵，再根据上文计算出关节空间的质量矩阵和阻尼矩阵，并通过关节刚度辨识实验得到关节空间刚度矩阵的值（表3）。

利用该雅可比矩阵将机器人系统动力学参数矩阵从机器人关节空间转换至机器人末端笛卡儿空间；计算模态矩阵，对末端笛卡儿空间的机器人系统动力学参数矩阵进行解耦使其为对角阵。当{\phi }_{\mathrm{1}}\mathrm{、}{\phi }_{\mathrm{2}}\mathrm{、}{\phi }_{\mathrm{3}}选取两组定值时，解耦坐标系中x轴方向的{\mathit{M}}_{\mathrm{P}x}和{\mathit{K}}_{\mathrm{P}x}的值随着{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}取值改变的变化情况见图8。由图8可知，机器人处于某一确定空间位置时，{\mathit{M}}_{\mathrm{P}x}值和{\mathit{K}}_{\mathrm{P}x}值会随着机器人加工姿态的改变而变化。对比可知，当机器人处于同一姿态而空间位置不同时，两种空间位置所对应的{\mathit{M}}_{\mathrm{P}x}和{\mathit{K}}_{\mathrm{P}x}不同，且{\mathit{M}}_{\mathrm{P}x}值和{\mathit{K}}_{\mathrm{P}x}值关于{\phi }_{\mathrm{4}}=\mathrm{0}和{\phi }_{\mathrm{5}}=\mathrm{0}呈一定规律的对称分布。在机器人所受外力一定的情况下，机器人系统所受惯性力较大时，机器人系统因杆件形变所受力会较小，即当{\mathit{M}}_{\mathrm{P}x}值较大时，{\mathit{K}}_{\mathrm{P}x}值较小。

4.2 颤振算例稳定性分析

4.2.1 再生型颤振算例分析

令{\phi }_{\mathrm{1}}=\mathrm{\pi }/\mathrm{2}rad、{\phi }_{\mathrm{2}}=\mathrm{\pi }/\mathrm{3}rad、{\phi }_{\mathrm{3}}=\mathrm{\pi }/\mathrm{4}rad，以（\mathrm{\pi }/\mathrm{18}）（10°）增量改变{\phi }_{\mathrm{4}}、{\phi }_{\mathrm{5}}的取值，当加工宽度系数{\mathit{K}}_{\mathrm{c}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{2000,2000,2000})、 V =diag（2/3，2/3，1/3）时，便可由式（26）计算出每个姿态下3个自由度所对应运动方向的极限加工宽度b_x、b_y、b_z。假设3个自由度所对应运动方向的进给深度为h_x、h_y、h_z，并将进给深度H用矩阵形式表达，即\mathit{H}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(h_x,h_y,h_z)。为了方便判断机器人系统是否发生再生型颤振，定义综合极限加工宽度b为极限加工宽度{\mathit{B}}_{\mathrm{l}\mathrm{i}\mathrm{m}}的二范数，即b=\sqrt[]{b_x^{\mathrm{2}}+b_y^{\mathrm{2}}+b_z^{\mathrm{2}}}；综合进给深度h为3个自由度所对应运动方向的进给深度 H 的二范数，即h=\sqrt[]{h_x^{\mathrm{2}}+h_y^{\mathrm{2}}+h_z^{\mathrm{2}}}，解耦坐标系中综合极限加工宽度b的变化如图9所示。设定该机器人在实际加工过程中的综合进给深度为0.3 mm，根据上述分析可知，实际加工宽度小于综合极限加工宽度的加工姿态，视为不发生再生型颤振的姿态。对应不发生再生型颤振时机器人作业姿态所对应的{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值如图10所示，图10中红色矩形框内的空洞表示发生再生型颤振的机器人加工姿态所对应的{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值。

4.2.2 振型耦合型颤振算例分析

本文所加工工件为汽车轮毂，其刚度通常在1000~1200 N/mm，考虑机器人加工系统时变特性及作业刀具形变和工件加工形变等综合影响，设置正比系数矩阵\mathit{Q}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{2000},\mathrm{2000},\mathrm{2000})，当{\phi }_{\mathrm{1}}=\mathrm{\pi }/\mathrm{2} rad、{\phi }_{\mathrm{2}}=\mathrm{\pi }/\mathrm{3} rad、{\phi }_{\mathrm{3}}=\mathrm{\pi }/\mathrm{4} rad时，以（\mathrm{\pi }/\mathrm{18}）（10°）增量改变{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值。将计算得到的机器人系统动力学参数代入式（31）中并求解其特征方程的根{\lambda }_i，特征根{\lambda }_i随{\phi }_{\mathrm{4}}、{\phi }_{\mathrm{5}}取值变化的分布情况见图11。由图11可知，特征根{\lambda }_i随{\phi }_{\mathrm{4}}、{\phi }_{\mathrm{5}}取值存在周期性的变化，这是由于关节4和5绕各自轴线旋转，当{\phi }_{\mathrm{5}}为某个定值，{\phi }_{\mathrm{4}}取旋转方向相反且两个旋转角度之和为\mathrm{2}\mathrm{\pi }时，机器人在实际运动中处于相同位姿；{\phi }_{\mathrm{4}}为某个定值时，{\phi }_{\mathrm{5}}的取值存在同样情况。

根据Routh-Hurwitz稳定判据，将特征方程的根均有负的实部时对应的机器人加工姿态视为不发生振型耦合型颤振的姿态。对应不发生振型耦合型颤振的机器人加工姿态所对应的{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值如图12所示，图12中红色矩形框内的空洞表示发生振型耦合型颤振的机器人加工姿态所对应的{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值。

4.2.3 颤振规避与加工姿态选取

从图10和图12中可以看出，当确定机器人末端空间位置（即确定关节1、2、3所对应转角{\phi }_{\mathrm{1}}、{\phi }_{\mathrm{2}}、{\phi }_{\mathrm{3}}的取值）时，将图10和图12上的稳定工作点对应的加工姿态作为加工作业的候选加工姿态，在实际加工中便可从这些候选加工姿态中进行选取。选择机器人加工姿态时，应尽量选择稳定工作点较为密集的区域，这样使得改变加工姿态调整更加平顺。与此同时，选择图10和图12重合区域对应加工姿态，可以同时规避再生型颤振和振型耦合型颤振。两类颤振都不发生的姿态所对应{\phi }_{\mathrm{4}}、{\phi }_{\mathrm{5}}的取值情况如图13所示，可以作为当前位置下机器人系统不发生颤振对应可选加工姿态的数据集，机器人每个可达空间位置都可以计算出相应的可选加工姿态数据集，从而对整个加工过程中的加工姿态进行规划并规避颤振。图13中红色矩形框里面的空洞表示两类颤振都发生的姿态所对应的{\phi }_{\mathrm{4}}\mathrm{、}{\phi }_{\mathrm{5}}的取值。

4.3 机器人模态解耦颤振规避对比

机器人末端执行器前3个关节处于{\phi }_{\mathrm{1}}=\mathrm{90}°、{\phi }_{\mathrm{2}}=\mathrm{60}°、{\phi }_{\mathrm{3}}=\mathrm{45}°位置，在解耦与未解耦情况下，对应的再生型颤振综合极限加工宽度b见图14，并统计极限加工宽度b处于0~1 mm对应的加工姿态数量，如图15所示；振型耦合型颤振所对应的可选加工姿态分布如图16所示。由图14可知，在该空间位置下，同一姿态序列所对应的解耦与未解耦的综合极限加工宽度重合率较低。由图15可知，在综合极限加工宽度为0~1 mm范围内，解耦形式和未解耦形式的综合极限加工宽度大部分都分布在0.2~0.6 mm处，其中解耦形式的综合极限加工宽度分布在0.2~0.4 mm处的加工姿态数量为1886个，占比36.3%，分布在0.4~0.6mm处的加工姿态数量为3067个，占比59.1%；未解耦形式的综合极限加工宽度分布在0.2~0.4 mm处的加工姿态数量为3949个，占比74.1%，分布在0.4~0.6mm处的加工姿态数量为1040个，占比19.5%。机器人末端执行器在前3个关节位于{\phi }_{\mathrm{1}}=\mathrm{90}°、{\phi }_{\mathrm{2}}=\mathrm{60}°、{\phi }_{\mathrm{3}}=\mathrm{45}°位置下不同综合进给深度时，不发生再生型颤振的可选择加工姿态数量如图16所示。

由图15和图16可知，采用未解耦形式进行再生型颤振规避，会导致某些稳定加工姿态被误认为是不稳定加工姿态，或者将某些不稳定姿态误认为是稳定姿态。例如，当综合进给深度为0.5 mm，此时需要保证所选机器人加工姿态对应的综合进给深度小于极限加工宽度，当机器人加工姿态对应的极限加工宽度小于综合进给深度0.5 mm时，该加工姿态为颤振姿态。由图16可知，采用解耦形式分析且综合进给深度0.5 mm小于极限加工宽度的可选加工姿态为1271个；而采用未解耦形式进行再生型颤振分析且综合进给深度0.5 mm小于极限加工宽度的可选加工姿态为661个。该情况下，采用未解耦形式的再生型颤振分析会导致较多的非颤振姿态被认为是颤振姿态而无法选取，降低了加工姿态规划的可选范围。当综合进给深度为0.6 mm时，采用解耦形式所得可选加工姿态为221个，而采用未解耦形式进行再生型颤振分析所得可选加工姿态为341个，前者比后者减少的加工姿态相对应的极限加工宽度均小于0.6 mm，因此，这些加工姿态实际为颤振姿态。如果实际的颤振姿态被误当成稳定姿态、且被错误选择为加工姿态，将引起再生型颤振风险，影响机器人系统的稳定性。

当机器人末端执行器在前3个关节位于相同位置下，针对振型耦合型颤振规避的可选择加工姿态数量如图17所示，与之对应，解耦和未解耦形式颤振姿态数量都为154个。对应机器人系统的特征方程的根中有一个或多个根具有正实部所对应的加工姿态，即认为系统发生振型耦合型颤振所对应的被判定为颤振姿态的加工姿态，占全部可达加工姿态的2.9%。由图17可知，该确定空间位置下，机器人模态解耦与未解耦情况下，不发生振型耦合型颤振的可选加工姿态完全重合，表明在进行振型耦合型颤振稳定性分析时，解耦与未解耦情况的求解结果一致。同时发现，在其他确定空间位置下的振型耦合型颤振稳定性，解耦与未解耦情况的求解结果仍一致。未解耦情况振型耦合型颤振稳定性判据的计算过程考虑了机器人动力学参数的耦合作用，对应解耦与未解耦情况的求解结果一致，表明所采用动力学模态解耦方法不会造成机器人动力学信息的缺失，验证了机器人动力学模态解耦的正确性。

5 实验验证

5.1 机器人空载状态下锤击实验

在图18所示的3T2R构型高灵巧性机器人实验平台^［11］开展颤振姿态与稳定姿态的机器人空载状态下锤击实验，验证所提动力学模态解耦的工业机器人加工颤振规避方法的有效性。实验中机器人振动信号测量系统由振动加速度传感器YD621-2和振动采集系统VK701H+构成，采样频率f为1000 Hz；激振采用力锤74002。如图18所示，将加速度传感器安装在靠近末端执行器的连杆上，以一定的时间间隔使用力锤锤击末端执行器以产生周期性的激振力，采集力信号和加速度振动信号，并对所采集信号进行后处理，得到频率响应函数曲线。

依据前述机器人颤振算例稳定性分析，将实验的机器人可达加工姿态分为稳定姿态和颤振姿态，在保证机器人末端空间位置不变的情况下，选择颤振算例分析中计算得到的3个稳定姿态与3个颤振姿态，如图19所示，以稳定姿态1和颤振姿态1为例进行锤击实验。采集到的振动加速度信号与激振力信号如图20所示，对振动加速度信号与激振力信号进行频域分析，得到稳定姿态1与颤振姿态1对应的频率响应函数曲线，如图21所示。3组稳定姿态和3组颤振姿态锤击实验对应频率响应函数的幅值和均值的实验结果见表4。

由图21可知，对比稳定姿态1和颤振姿态1中频域的振动响应值，稳定姿态1对应频率响应函数的幅值为4.839g，均值为0.064g，而颤振姿态1对应频率响应函数的幅值为77.063g，均值为1.176g；与颤振姿态1相比，稳定姿态1的幅值和均值分别减小了93.7%和94.5%。

由表4实验结果可知，稳定姿态的频率响应函数的幅值和均值均小于颤振姿态，3个稳定姿态的频率响应幅值和频率响应均值的平均值分别为4.758 g和0.378 g；相应的3个颤振姿态对应值分别为 53.585 g和1.070 g；前者与后者相比，分别减小了91.1%和64.7%。实验结果表明，颤振姿态比稳定姿态在特定频率下更易受到外界激励的影响，从而导致机器人颤振的发生，与颤振算例稳定性分析中将该加工姿态判定为颤振姿态的结果相符，验证了所提动力学模态解耦的工业机器人加工颤振规避方法的有效性。

5.2 机器人实验平台加工实验

为进一步验证颤振规避方法在机器人在加工状态下的有效性，在相同机器人实验平台开展颤振姿态与稳定姿态的机器人接触被加工工件下颤振规避实验，如图22所示。实验中机器人振动信号测量系统由振动加速度传感器YD621-2和振动采集系统VK701H+构成，采样频率为1000 Hz，其中加速度传感器安装在靠近末端执行器的连杆上。选择图19所示的稳定姿态和颤振姿态进行打磨加工作业，以稳定姿态1和颤振姿态1为例，实际机器人加工状态下颤振规避实验对应姿态如图23所示。

控制机器人在稳定姿态1和颤振姿态1进行打磨加工作业，利用机器人TwinCAT控制平台和振动采集系统采集机器人关节力矩信号和加速度振动信号，通过机器人雅可比矩阵将机器人关节力矩映射至末端笛卡儿空间六维力/力矩形式，机器人加工实验的稳定姿态1与颤振姿态1所采集的关节力矩信号如图24所示。对所采集信号进行后处理，得到机器人稳定姿态1与颤振姿态1实际加工作业时所对应的机器人单位力作用下加速度信号a_{\mathrm{0}}（称为单位加速度），如图25所示。

由图25可知，对比稳定姿态1和颤振姿态1对应的单位加速度，前2s机器人未接触被加工工件，稳定姿态1和颤振姿态1中单位加速度大小都在0.020g/N附近波动；2~3s为机器人与被加工工件接触阶段，该阶段稳定姿态1和颤振姿态1对应单位加速度的幅值分别为0.032g/N和0.043g/N，均值分别为0.023g/N和0.025g/N；3~10s为机器人稳定打磨加工工件阶段，该阶段稳定姿态1和颤振姿态1对应单位加速度的幅值分别为0.028g/N和0.034g/N，均值分别为0.018g/N和0.022g/N。3组稳定姿态和3组颤振姿态机器人与被加工工件接触阶段对应单位加速度的幅值和均值的实验结果见表5。

由表5可知，稳定姿态的机器人与被加工工件接触阶段对应单位加速度的幅值和均值均小于颤振姿态，3个稳定姿态的机器人与被加工工件接触阶段对应单位加速度的幅值和均值的平均值分别为0.026g/N和0.019g/N；相应的3个颤振姿态对应值分别为0.061g/N和0.046g/N。实验结果表明，在单位力作用情况下，颤振姿态比稳定姿态会产生更大的加速度，从而导致机器人颤振的发生，与颤振算例稳定性分析中将该加工姿态判定为颤振姿态的结果相符，进一步验证了所提方法的有效性。

5.3 颤振稳定性判据正确性验证

因机器人加工过程中无法始终保持设置的加工工况严格一致，且在实际加工过程中无法实时计算两类颤振的稳定性判据，在颤振稳定性判据计算中，忽略机器人加工过程中加工工况的实时变动。将表5所示稳定姿态1和颤振姿态1对应的机器人系统动力学参数矩阵的值和末端笛卡儿空间六维力/力矩代入式（28）和式（31），机器人加工实验稳定姿态和颤振姿态所对应综合极限打磨宽度分别为4.6 mm和3.1 mm，机器人加工实验所设置的打磨宽度为1 mm，所设置打磨宽度小于综合极限打磨宽度，两姿态下机器人均不发生再生型颤振；机器人加工实验稳定姿态对应动力学微分方程系统特征方程的根的实部分别为{\lambda }_{\mathrm{1}}=-\mathrm{33.9}、{\lambda }_{\mathrm{2}}=-\mathrm{0.002}、{\lambda }_{\mathrm{3}}=-\mathrm{0.007}，均为负的实部，该姿态下机器人不发生振型耦合型颤振；机器人加工实验颤振姿态对应动力学微分方程系统特征方程的根的实部分别为{\lambda }_{\mathrm{1}}=-\mathrm{33.9}、{\lambda }_{\mathrm{2}}=-\mathrm{0.015}、{\lambda }_{\mathrm{3}}=\mathrm{0.023}，{\lambda }_{\mathrm{3}}为正的实部，该姿态下机器人发生振型耦合型颤振。

上述稳定性判据计算结果表明，颤振姿态下机器人加工过程中发生了颤振，与机器人加工实验中在单位力作用情况下，颤振姿态比稳定姿态会产生更大的加速度，从而导致机器人颤振发生的实验结果相符，验证了所提方法的再生型颤振和振型耦合型颤振的稳定性判据的正确性。

6 结论

1）针对工业机器人加工作业中的颤振规避问题，提出一种动力学模态解耦的工业机器人加工颤振规避方法。建立工业机器人加工姿态与机器人系统动力学模型参数的映射关系，对机器人动力学微分方程进行模态解耦，降低机器人颤振稳定性分析的复杂程度；确定了机器人颤振问题的主方向，从而可以针对性地对主方向的颤振进行规避；开展再生型颤振和振型耦合型颤振稳定性分析，获得两类颤振的稳定性判据。所提方法可以单一规避再生型颤振或振型耦合型颤振，也可以同时规避这两类颤振，为实现工业机器人加工颤振规避及机器人加工姿态选取提供了一种实现途径。

2）通过机器人颤振算例稳定性分析，验证了机器人动力学微分方程模态解耦方法的有效性。与未解耦形式相比，采用解耦形式的再生型颤振稳定性分析，可更准确地判定稳定姿态或颤振姿态，减小误判风险；对于振型耦合型颤振稳定性分析，解耦和未解耦形式的颤振姿态数量相同，且解耦和未解耦形式对应的不发生振型耦合型颤振的稳定姿态完全重合，验证了解耦过程的正确性。机器人动力学微分方程模态解耦方法可以有效判定加工作业姿态中的稳定姿态与颤振姿态，为机器人加工颤振规避提供了一种理论分析方式。

3）在机器人实验平台开展颤振姿态与稳定姿态的机器人锤击实验和机器人加工状态下的颤振规避实验，验证了所提动力学模态解耦的工业机器人加工颤振规避方法的有效性。

4）在机器人加工作业颤振规避的基础上，后续将进一步研究机器人加工颤振识别，以及机器人作业刚度优化和形变误差补偿，综合减小机器人加工作业过程中的末端形变，提高机器人的加工精度。

参考文献

[1] 李振， 赵欢， 王辉， 等. 机器人磨抛加工接触稳态自适应力跟踪研究［J］. 机械工程学报， 2022， 58（9）：200-209.

[2] LI Zhen， ZHAO Huan， WANG Hui， et al. Research on Contact Steady-state Adaptive Force Tracking of Robot Grinding and Polishing［J］. Journal of Mechanical Engineering， 2022， 58（9）：200-209.

[3] 孙朝阳， 彭芳瑜， 唐小卫， 等. 基于自适应变分模态分解与功率谱熵差的机器人铣削加工颤振类型辨识［J］. 机械工程学报， 2023， 59（9）：90-100.

[4] SUN Zhaoyang， PENG Fangyu， TANG Xiaowei， et al. Robotic Milling Chatter Types Detection Based on Adaptive Variational Mode Decomposition and Difference of Power Spectral Entropy［J］. Journal of Mechanical Engineering， 2023， 59（9）：90-100.

[5] 叶松涛， 严思杰， 李文韬， 等. 面向机器人铣削加工的刀尖动态特性分析与稳定性预测［J］. 机械工程学报， 2022， 58（17）：261-275.

[6] YE Songtao， YAN Sijie， LI Wentao， et al. Analysis of Tool Tip Dynamic Characteristics and Stability Prediction for Robotic Milling Tasks［J］. Journal of Mechanical Engineering， 2022， 58（17）：261-275.

[7] WU Jiawei， TANG Xiaowei， PENG Fangyu， et al. A Novel Mode Coupling Mechanism for Predicting Low-frequency Chatter in Robotic Milling by Providing a Vibration Feedback Perspective［J］. Mechanical Systems and Signal Processing， 2024， 216：111424.

[8] YANG Yun， YUAN Jiawei， Duo TIE， et al. An Efficient and Accurate Chatter Prediction Method of Milling Processes with a Transition Matrix Reduction Scheme［J］. Mechanical Systems and Signal Processing， 2023， 182：109535.

[9] WANG Yu， ZHANG Mingkai， TANG Xiaowei， et al. A kMap Optimized VMD-SVM Model for Milling Chatter Detection with an Industrial Robot［J］. Journal of Intelligent Manufacturing， 2022， 33（5）：1483-1502.

[10] YANG Bin， GUO Kai， ZHOU Qian， et al. Early Chatter Detection in Robotic Milling under Variable Robot Postures and Cutting Parameters［J］. Mechanical Systems and Signal Processing， 2023， 186：109860.

[11] 梁志强， 陈司晨， 杜宇超， 等. 基于变分模态滤波和注意力机制的重载机器人铣削系统颤振辨识方法［J］. 中国机械工程， 2025， 36（5）：1018-1027.

[12] LIANG Zhiqiang， CHEN Sichen， DU Yuchao， et al. Chatter Identification Method for Heavy-duty Robotic Milling Systems Based on Variational Mode Filtering and Attention Mechanism［J］. China Mechanical Engineering， 2025， 36（5）：1018-1027.

[13] ZHAO Wei， LI Bo， TIAN Wei， et al. Magnetorheological Elastomer Absorber-based Chatter Suppression in Robotic Milling［J］. Robotics and Computer-Integrated Manufacturing， 2024， 88：102740.

[14] TAKAHASHI W， SUZUKI N， SHAMOTO E. Development of a Novel Boring Tool with Anisotropic Dynamic Stiffness to Avoid Chatter Vibration in Cutting Part 1： Design of Anisotropic Structure to Attain Infinite Dynamic Stiffness［J］. Precision Engineering， 2021， 68：57-71.

[15] 郭万金. 复杂形状零部件打磨作业机器人研究［D］. 哈尔滨：哈尔滨工业大学， 2017.

[16] GUO Wanjin. Research on Robot Manipulator for Complex Shaped Parts Deburring［D］. Harbin：Harbin Institute of Technology， 2017.

[17] TLUSTY J， POLACEK M. The Stability of the Machine Tool against Self Excited Vibration in Machining［C］∥ Proceedings of International Research in Production Engineering. Pittsburgh， 1963：465-474.

[18] MERRITT H E. Theory of Self-excited Machine-tool Chatter：Contribution to Machine-tool Chatter Research—1［J］. Journal of Engineering for Industry， 1965， 87（4）：447-454.

[19] 郭万金， 李儒， 王利斌， 等. 一种考虑弱刚性关节传动链结构和柔性杆结构的机器人刚度辨识方法［J］. 中国机械工程， 2025， 36（5）：974-985.

[20] GUO Wanjin， LI Ru， WANG Libin， et al. A Robot Stiffness Identification Method Considering Weakly Rigid Joint Transmission Chain Structures and Flexible Rod Structures［J］. China Mechanical Engineering， 2025， 36（5）：974-985.

[21] 陈钦韬， 殷参， 张加波， 等. 面向铣削任务的工业机器人刚度位姿优化［J］. 机器人， 2021， 43（1）：90-100.

[22] CHEN Qintao， YIN Shen， ZHANG Jiabo， et al. Pose Optimization of Industrial Robots Based on Stiffness for Milling Tasks［J］. Robot， 2021， 43（1）：90-100.

[23] 刘延柱， 陈立群， 陈文良. 振动力学［M］. 2版. 北京：高等教育出版社， 2011.

[24] LIU Yanzhu， CHEN Liqun， CHEN Wenliang. Mechanics of Vibration［M］. 2nd ed. Beijing：Higher Education Press， 2011.

[25] ZHANG Yun， ZHOU Tingting， HU Tianliang， et al. Stability Analysis of Milling Chatter in Six-degree-of-freedom Industrial Robots［J］. The International Journal of Advanced Manufacturing Technology， 2023， 127（5）：2861-2880.

[26] JI Yongjian， WANG Liyong， SONG Yue， et al. Investigation of Robotic Milling Chatter Stability Prediction under Different Cutter Orientations by an Updated Full-discretization Method［J］. Journal of Sound and Vibration， 2022， 536：117150.

[27] 廖文和， 郑侃， 孙连军， 等. 大型复杂构件机器人加工稳定性研究进展［J］. 航空学报， 2022， 43（1）：026061.

[28] LIAO Wenhe， ZHENG Kan， SUN Lianjun， et al. Review on Chatter Stability in Robotic Machining for Large Complex Components［J］. Acta Aeronautica et Astronautica Sinica， 2022， 43（1）：026061.

[29] OZAKI N， TANIGUCHI T， HIROGAKI T， et al. Elucidation of the Mechanism of Chatter Mark Formation during End-milling and Inverse Analysis of Chatter Vibration by Two-dimensional Discrete Fourier Transform of Chatter Mark［J］. Precision Engineering， 2023， 82：116-128.

[30] MERCHANT M E. Basic Mechanics of the Metal-cutting Process［J］. Journal of Applied Mechanics， 1944， 11（3）：A168-A175.

基金资助

国家自然科学基金(52275005)

陕西省自然科学基础研究计划(2025JC-QYXQ-027)

陕西省自然科学基础研究计划(2025JC-YBMS-619)

中央高校基本科研业务费专项资金(300102253201)

中国博士后科学基金(2024M760002)

安徽省博士后科研项目(2023B675)

安徽省机器视觉检测与感知重点实验室开放基金(KLMVI-2025-HIT-06)