螺旋灸法机器人的系统设计及控制方法

System Design and Control Method of Spiral Moxibustion Robot

摘要

雷火灸属于艾灸的一种， 具有药力猛、 渗透力强的特点， 螺旋灸法作为雷火灸的一种手法， 对青少年近视、 干眼症和过敏性鼻炎等疾病有着良好的治疗效果。但是， 利用人工进行螺旋灸法存在精确度不高、 温度控制不准的问题， 且目前可替代人工进行艾灸治疗的设备都不能实现螺旋灸法。针对以上问题， 设计了一种控制方法来模拟中医的螺旋艾灸方法， 并将其应用于艾灸机器人。首先， 使用CoppeliaSim机器人仿真平台对机械臂进行螺旋灸法的轨迹规划仿真； 然后， 采用改进D-H参数法建立机器人的运动学方程， 使用正运动学分析和三次多项式轨迹规划来控制真实艾灸机器人进行螺旋灸法的运动； 最后， 通过真实机械臂测试得到了30次螺旋灸法的实际轨迹误差， 均在误差允许的范围之内。因此， 所提螺旋灸法控制方法的运动精度有效地满足了大多数艾灸应用的要求。

Abstract

Thunder-fire moxibustion belongs to a kind of moxibustion. The efficacy of thunder-fire moxibustion is closely related to the moxibustion methods. Different ailments require different methods. Spiral moxibustion， as a thunder-fire moxibustion method， has a great therapeutic effect on many ailments. In addition， manually spiral moxibustion has the problem of poor accuracy and temperature control， but the current equipment that can replace manual moxibustion treatment can not realize the spiral moxibustion method. In this study， a control method is designed to simulate the spiral moxibustion method of traditional Chinese medicine， which was applied to moxibustion robots. CoppeliaSim is used to simulate the trajectory planning of the robotic arm for spiral moxibustion and the Modified D-H parametric method is used to establish the kinematic equations of the robot. Then the positive kinematics analysis and cubic polynomial trajectory planning are used to control the real moxibustion robot. Finally， through real robotic arm testing， the actual trajectory error of 30 spiral moxibustion methods is obtained， which is within the allowable range of error. Therefore， the motion accuracy of the spiral moxibustion control method meets the requirements of most moxibustion applications effectively.

Graphical abstract

关键词

雷火灸 / 螺旋灸法 / 艾灸机器人 / 轨迹规划 / 运动学分析

Key words

thunder-fire moxibustion / spiral moxibustion / moxibustion robot / trajectory planning / kinematic analysis

0 引　言

艾灸已被广泛应用于疾病的预防和治疗中， 其利用艾草燃烧产生的热力和药效， 通过对人体特定穴位的刺激， 达到预防疾病、 调整身体机能的目的^［1］。其中， 雷火灸作为一种新型的艾灸疗法， 使用独特的艾条^［2］， 在燃烧时产生更强的热辐射力、 远近红外线和药力因子。

雷火灸的疗效与艾灸手法密切相关， 不同的疾病使用不同的治疗手法^［3］。随着中医的成熟和规范化， 国家中医药管理局于2013年发布了《中医医疗技术手册》， 其中说明了雷火灸的基本操作方法不仅包括常见的雀啄灸、 回旋灸和往复灸等基本灸法， 还包括螺旋灸法^［4］。同时， 艾灸的热变化对镇痛有显著效果， 可以在治疗过程中起到关键效果^［5］。因此， 本文模拟了二维螺旋灸法和三维螺旋灸法， 其中， 三维螺旋灸法在治疗过程中的波动温度可以更好地为穴位区域提供动态和反复的热刺激。

近年来， 艾灸技术与机器人技术相结合， 艾灸机器人得到了快速发展。王新庆等^［6］采用工控机+STM32的分层系统架构设计了一款智能艾灸机器人， 采用空间插补的方式实现回旋灸、 循径灸、 雀啄灸的艾灸手法设计； 马蓓蓓等^［7］采用主从双边结构， 设计了一种艾灸机器人系统， 可以实现腧穴的温和灸； 白善明等^［8］使用机器人操作系统下的MoveIt平台对机械臂运动学控制进行快速配置， 并引入笛卡尔轨迹规划和多项式插值等方法， 实现了摆尾灸、 回旋灸、 雀啄灸的艾灸手法设计。这些艾灸机器人都做到了精准化、 标准化、 简便化施灸， 证明了艾灸机器人代替传统人工艾灸的可行性。

然而， 上述自动化艾灸机器人均只能实现简单的轨迹规划灸法， 如回旋灸、 雀啄灸等， 并不能实现更为复杂的螺旋灸法， 这也导致了目前已有艾灸机器人对青少年近视、 干眼症、 过敏性鼻炎等疾病不具备良好的治疗效果， 限制了艾灸机器人在临床中的推广和应用。

因此， 本文设计了一种能模拟传统螺旋灸法的控制方法， 并开发了一种可替代人工进行螺旋灸法的机器人。艾灸机器人采用工控机＋STC51单片机的分层结构， 可以简单有效地控制机械臂的运动、 除灰滤烟模块、 温度检测与反馈等系统。因此， 本文使用CoppeliaSim机器人仿真软件进行轨迹规划， 并使用了改进D-H参数法、 正运动学分析以及三次多项式轨迹规划等控制方法， 不仅可以让艾灸机器人复刻常用的灸法， 而且可以完成更加复杂的轨迹规划， 如平面螺旋灸和三维螺旋灸。除灰滤烟模块可以减少艾烟中有害气体对医生和患者的影响， 使艾灸理疗更加绿色安全。温度监测与反馈系统可以对艾灸穴位的皮肤表面温度进行实时监测， 由此调节机械臂轨迹运行的高度， 使患者穴位的皮肤表面温度一直保持在最佳范围内， 提高了艾灸机器人操作过程的安全性、 智能性和友好性。

1 艾灸机器人控制系统设计

艾灸机器人采用工控机＋STC51单片机的分层结构， 使用工控机搭载上位机软件系统， 完成机械臂轨迹规划、 定时清灰和温度数据监测等功能， 使用STC51单片机作为下位机系统， 负责采集并发送温度数据等功能。艾灸机器人的整体结构如图 1 所示， 核心是六自由度机械臂， 可以像人类手臂一样做出灵活的动作， 适合于几乎任何轨迹的工作^［9］。

艾灸机器人整体搭配了一体机、 急停按钮、 吹气泵和烟雾净化器等装置。同时， 还为机械臂末端设计了灸疗器具， 用于放置艾柱， 底部有温度传感器和红外线， 用来监测人体皮肤表面的温度和辅助定位穴位。

1.1 灸疗器具设计

为了使灸疗过程更加方便、 高效和安全， 设计了灸疗器具， 安装在机械臂末端， 用于艾条的放置。灸疗器具由内筒和外筒构成， 如图 2（a） 所示， 集排烟、 除灰、 防烫、 温度监测、 红外线辅助穴位定位等多种功能于一体， 各部分均可拆卸， 能够克服灸疗过程中容易烫伤的问题， 减少艾烟气味对医患及诊疗环境的刺激和污染。

内筒为圆柱形筒体， 由筒体、 隔灰网和顶部的通风孔构成， 筒体内部放置艾条。隔灰网在筒体的底部， 由铁丝网构成， 用于拦截掉落的艾灰。同时， 艾灸机器人搭配了吹气泵和网络继电器等装置， 实现了定时清灰功能， 可有效阻止艾灰落到皮肤上造成烫伤。艾烟中有200多种成分， 其中， 单环芳香烃、 甲醛、 多环芳香烃等对人体有害， 艾烟浓度较低时对人体的影响较小， 但长期暴露于高浓度艾烟中会对人体产生一定的损伤^［10］。因此， 筒体顶部设置了通风孔， 通过管道与烟雾净化器连接， 灸疗过程中可以通过烟雾净化器吸收艾烟， 尽可能减少散发至空气中的艾烟， 提高环境空气质量。

外筒整体为圆柱形筒体， 侧面有用于安装电控组件的长方体放置槽， 筒体底部安装红外测温模块和红外线定位装置， 如图 2（b） 所示。其中， 电控组件由STC89C51单片机和无线传输模块构成。

1.2 控制系统设计

艾灸机器人进行灸疗的流程如图 3 所示。首先， 选择艾灸穴位点， 将灸疗器具底部的中心对准需要灸疗的穴位， 并记录该穴位的坐标。

然后， 选择一种或几种艾灸治疗手法， 并在上位机程序中设置对应的参数， 如灸疗的时间、 回旋灸的半径、 雀啄灸的高度、 平面螺旋灸的选择半径范围或三维螺旋灸的旋转高度等， 即可规定机械臂轨迹规划的运动范围。对于回旋灸、 雀啄灸等简单灸法， 只需要重复进行圆形、 线性轨迹运动至施灸时间结束即可， 而对于螺旋灸法等复杂的轨迹运行， 机械臂会在进行一轮完整的轨迹运动之后， 从设置的螺旋灸法最大半径处再次移动至穴位中心点， 开启新一轮的螺旋灸法运动， 如此重复至灸疗时间结束。

与此同时， 为了让机械臂在整个施灸过程中与人体保持在安全的距离范围内， 设计了温度监测与反馈系统。通过灸疗器具底部的红外温度传感器， 对艾灸穴位的皮肤表面温度进行实时监测， 并通过无线传输模块将温度数据上传到控制系统进行运算， 求解出温度变化的范围和速度。若艾灸温度不符合适宜范围， 则改变运行的位置高度， 对艾灸机器人的手法轨迹进行重新规划， 从而确保艾灸过程中， 患者穴位的皮肤表面温度一直保持在最佳范围内， 有效防止了燃烧的艾条烫伤人体。与此同时， 机械臂设置了碰撞保护功能， 只要检测到机械臂末端与物体发生微小碰撞， 机械臂就会紧急停止运动， 并且恢复至初始位置， 大大提高了艾灸机器人的安全性。

2 螺旋灸法的控制方法

CoppeliaSim是一款基于分布式控制架构， 具有集成开发环境的机器人仿真器， 包含快速可自定义的动力学计算、 任何类型机制的运动学计算、 碰撞检测等功能^{［11‐12］}。在对艾灸机器人进行螺旋灸的轨迹规划和控制前， 需要完成机械臂的基本建模和配置： 首先， 创建机械臂的STL模型， 并将其导入到CoppeliaSim进行特征提取及关节添加， 同时对动力学特性及碰撞特性进行配置， 防止机械臂自身连杆之间碰撞； 然后， 为机械臂配置运动学求解器， 运动学算法是机械臂各类算法的核心， 它包含了正运动学求解和逆运动学求解等算法， 本次使用逆运动学求解工具IK来进行机械臂的运动学解算； 之后， 在CoppeliaSim中完成螺旋灸的轨迹规划仿真， 并导出机械臂的实时关节变化数据； 最后， 根据正运动学和三次多项式轨迹规划算法， 将各个关节角度变化转化为真实机械臂的位姿变化， 使真实机械臂完成螺旋灸的运动轨迹。

2.1 基于艾灸手法的螺旋灸轨迹规划仿真

平面螺旋灸和三维螺旋灸的轨迹规划参考阿基米德螺线， 如图 4 所示。阿基米德螺线是一个点匀速离开一个固定点的同时， 又以固定的角速度绕该固定点转动而产生的轨迹， 它的的极坐标方程式为

r=a+b\theta , （1）

式中： r为与极点的距离； θ为与极轴的角度， 极轴就是在平面直角坐标系中的x轴正方向； a和b均为实数。并且半径随着角度变化而线性变化， 连续两条螺旋线之间的间距是均匀的， 为Δr=2πb^［13］。这就保证了艾灸过程中火头运行的轨迹是均匀的， 即温度变化也是均匀的。

基于阿基米德螺线的特点， 可以在CoppeliaSim中规划平面螺旋和三维螺旋的运行轨迹。首先， 在CoppeliaSim中建模的六自由机械臂基础上， 在机械臂末端关节的Z轴方向上添加一个旋转关节和两个平移关节， 分别控制机械臂末端进行固定角度旋转和水平、 竖直方向平移， 如图 5 所示。

使用Reflex运动库同时移动多个关节， 设置参数围绕旋转关节以固定的弧度运动， 同时向平移关节1沿机械臂末端X轴正方向以固定的速度进行移动， 即可实现顺时针平面螺旋运动， 运动仿真轨迹如图 6（a） 所示。此外， 在围绕旋转关节和平移关节1按上述移动时， 再向平移关节2沿机械臂末端Z轴正方向以固定的速度进行移动， 即可实现顺时针三维螺旋运动， 运动仿真轨迹如图 6（b） 所示。

根据中国国家中医药管理局发布的《中医医疗技术手册》， 可以得到螺旋灸法的具体参数。如图 7 所示， 艾灸火头对准应灸部位中心点， 一般火头距离皮肤2~3 cm， 做顺时针方向螺旋式旋转， 旋转直径为1~5 cm^［4］。设置平面螺旋和三维螺旋运动仿真轨迹的参数， 使其旋转直径和旋转高度符合螺旋灸法所规定的参数。

围绕旋转关节运动的固定弧度影响着仿真轨迹的拟合度以及仿真后导出机械臂关节角度的数据量， 为此设置了多种弧度参数进行对比。结果表示， 当旋转弧度设置为1/12π时， 三维螺旋的运行轨迹已经基本光滑， 并且仿真后导出的关节角度数据的数量合适。因此， 旋转弧度设置为1/12π， 此时能确定平面螺旋和三维螺旋运动仿真轨迹的各个参数， 能够将螺旋灸的仿真关节轨迹和实时关节角度变化数据导出。

2.2 正运动学分析

机器人运动学是从几何角度研究机器人的位置、 速度和加速度随时间变化规律的科学， 主要在机器人的工作空间与关节空间中讨论， 即分析机器人的关节变量与末端位置姿态之间的关系^{［14‐15］}。

本文所涉及的逆运动学求解使用CoppeliaSim中求解工具IK， 因此只讨论正运动学分析， 即根据机器人各关节变量的变化来计算机器人末端的位置姿态^［16］。本文采用改进的Denavit-Hartenberg矩阵构建螺旋灸机器人的运动学模型^{［17‐19］}， 机械臂的连杆坐标如图 8 所示。

在图 8 中， （X₀，Y₀，Z₀）为机器人基础坐标系， （X_i，Y_i，Z_i ）（i=1，2，…，6）为机器人每个关节的坐标系。根据机器人运动模型， 可以获得机器人的改进D-H参数， 如表 1 所示。

表 1 中， i为关节的标号， α_i-1为连杆的旋转角度， a_i-1为连杆的长度， d_i 为连杆的偏移距离， θ_i 为关节角度。

通过定义连杆坐标系和相应的连杆参数， 可以直接建立运动学方程^［20］。相邻连杆坐标系之间的变换可以通过坐标系相对位置和姿态的旋转和平移来实现， 用齐次变换矩阵表示为

{}_i{}^{i-\mathrm{1}}\mathit{T}=R_x({\alpha }_{i-\mathrm{1}})D_x(a_{i-\mathrm{1}})R_z({\theta }_i)D_z(d_i),（2）

\begin{array}{c}{}_{}{}^{i-\mathrm{1}}\mathit{T}=\\ \left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i& -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i& \mathrm{0}& a_{i-\mathrm{1}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& -d_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& d_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\end{array}（3）

式中： {}_i{}^{i-\mathrm{1}}\mathit{T}为连杆i相对于连杆i-1的均匀变换矩阵； R_x （）和R_z （）为坐标系围绕X轴和Z轴的旋转函数； D_x （）和D_z （）为坐标系围绕X轴与Z轴的平移函数。将表 1 中各连杆参数代入式（2）可以得到各连杆的变换矩阵为

{}_{\mathrm{1}}{}^{\mathrm{0}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{1}}& -s_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ s_{\mathrm{1}}& c_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{240.5}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

{}_{\mathrm{2}}{}^{\mathrm{1}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{2}}& -s_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}\\ s_{\mathrm{2}}& c_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

{}_{\mathrm{3}}{}^{\mathrm{2}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{3}}& -s_{\mathrm{3}}& \mathrm{0}& \mathrm{256}\\ s_{\mathrm{3}}& c_{\mathrm{3}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

{}_{\mathrm{4}}{}^{\mathrm{3}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{4}}& -s_{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{210}\\ s_{\mathrm{4}}& c_{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

{}_{\mathrm{5}}{}^{\mathrm{4}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{5}}& -s_{\mathrm{5}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -s_{\mathrm{5}}& -c_{\mathrm{5}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],

{}_{\mathrm{6}}{}^{\mathrm{5}}\mathit{T}=\left[\begin{array}{cccc}{c}_{\mathrm{6}}& -s_{\mathrm{6}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{144}\\ s_{\mathrm{6}}& c_{\mathrm{6}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]，（4）

式中： c_i为\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i的缩写； s_i为\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i的缩写； T 为对应的变化矩阵； {}_{\mathrm{1}}{}^{\mathrm{0}}\mathit{T}， {}_{\mathrm{2}}{}^{\mathrm{1}}\mathit{T}， {}_{\mathrm{3}}{}^{\mathrm{2}}\mathit{T}， {}_{\mathrm{4}}{}^{\mathrm{3}}\mathit{T}， {}_{\mathrm{5}}{}^{\mathrm{4}}\mathit{T}， {}_{\mathrm{6}}{}^{\mathrm{5}}\mathit{T}分别为每个关节的变化矩阵。将这些连杆变换矩阵连乘就能得到机械臂末端相对于基坐标系位置姿态的变换矩阵{}_{\mathrm{6}}{}^{\mathrm{0}}\mathit{T}， 如式（5）所示。

{}_{\mathrm{6}}{}^{\mathrm{0}}\mathit{T}={}_{\mathrm{1}}{}^{\mathrm{0}}\mathit{T}{}_{\mathrm{2}}{}^{\mathrm{1}}\mathit{T}{}_{\mathrm{3}}{}^{\mathrm{2}}\mathit{T}{}_{\mathrm{4}}{}^{\mathrm{3}}\mathit{T}{}_{\mathrm{5}}{}^{\mathrm{4}}\mathit{T}{}_{\mathrm{6}}{}^{\mathrm{5}}\mathit{T}=\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],（5）

式中： {}_{\mathrm{1}}{}^{\mathrm{0}}\mathit{T}， {}_{\mathrm{2}}{}^{\mathrm{1}}\mathit{T}， {}_{\mathrm{3}}{}^{\mathrm{2}}\mathit{T}， {}_{\mathrm{4}}{}^{\mathrm{3}}\mathit{T}， {}_{\mathrm{5}}{}^{\mathrm{4}}\mathit{T}， {}_{\mathrm{6}}{}^{\mathrm{5}}\mathit{T}分别为每个关节的变化矩阵。将{}_{\mathrm{1}}{}^{\mathrm{0}}\mathit{T}， {}_{\mathrm{2}}{}^{\mathrm{1}}\mathit{T}， {}_{\mathrm{3}}{}^{\mathrm{2}}\mathit{T}， {}_{\mathrm{4}}{}^{\mathrm{3}}\mathit{T}， {}_{\mathrm{5}}{}^{\mathrm{4}}\mathit{T}， {}_{\mathrm{6}}{}^{\mathrm{5}}\mathit{T}代入式（5）， 可以得到

\left\{\begin{array}{l}{n}_x=c_{\mathrm{6}}(c_{\mathrm{5}}(c_{\mathrm{1}}{c}_{\mathrm{23}}{c}_{\mathrm{4}}+s_{\mathrm{1}}{s}_{\mathrm{4}})-c_{\mathrm{1}}{s}_{\mathrm{23}}{s}_{\mathrm{5}})+(c_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{c}_{\mathrm{23}}{s}_{\mathrm{4}})s_{\mathrm{6}},\\ n_y=c_{\mathrm{6}}(c_{\mathrm{5}}(c_{\mathrm{23}}{c}_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{s}_{\mathrm{4}})-s_{\mathrm{1}}{s}_{\mathrm{23}}{s}_{\mathrm{5}})+(-c_{\mathrm{1}}{c}_{\mathrm{4}}-c_{\mathrm{23}}{s}_{\mathrm{1}}{s}_{\mathrm{4}})s_{\mathrm{6}},\\ n_z=c_{\mathrm{6}}(c_{\mathrm{4}}{c}_{\mathrm{5}}{s}_{\mathrm{23}}+c_{\mathrm{23}}{s}_{\mathrm{5}})-s_{\mathrm{23}}{s}_{\mathrm{4}}{s}_{\mathrm{6}},\\ o_x=c_{\mathrm{6}}(c_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{c}_{\mathrm{23}}{s}_{\mathrm{4}})-(c_{\mathrm{5}}(c_{\mathrm{1}}{c}_{\mathrm{23}}{c}_{\mathrm{4}}+s_{\mathrm{1}}{s}_{\mathrm{4}})-c_{\mathrm{1}}{s}_{\mathrm{23}}{s}_{\mathrm{5}})s_{\mathrm{6}},\\ o_y=c_{\mathrm{6}}(-c_{\mathrm{1}}{c}_{\mathrm{4}}-c_{\mathrm{23}}{s}_{\mathrm{1}}{s}_{\mathrm{4}})-(c_{\mathrm{5}}(c_{\mathrm{23}}{c}_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{s}_{\mathrm{4}})-s_{\mathrm{1}}{s}_{\mathrm{23}}{s}_{\mathrm{5}})s_{\mathrm{6}},\\ o_z=-c_{\mathrm{6}}{s}_{\mathrm{23}}{s}_{\mathrm{4}}-(c_{\mathrm{4}}{c}_{\mathrm{5}}{s}_{\mathrm{23}}+c_{\mathrm{23}}{s}_{\mathrm{5}})s_{\mathrm{6}},\\ a_x=c_{\mathrm{1}}{c}_{\mathrm{5}}{s}_{\mathrm{23}}+(c_{\mathrm{1}}{c}_{\mathrm{23}}{c}_{\mathrm{4}}+s_{\mathrm{1}}{s}_{\mathrm{4}})s_{\mathrm{5}},\\ a_y=c_{\mathrm{5}}{s}_{\mathrm{1}}{s}_{\mathrm{23}}+(c_{\mathrm{23}}{c}_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{s}_{\mathrm{4}})s_{\mathrm{5}},\\ a_z=-c_{\mathrm{23}}{c}_{\mathrm{5}}+c_{\mathrm{4}}{s}_{\mathrm{23}}{s}_{\mathrm{5}},\\ p_x=\mathrm{256}{c}_{\mathrm{1}}{c}_{\mathrm{2}}+\mathrm{210}{c}_{\mathrm{1}}{s}_{\mathrm{23}}-\mathrm{144}(c_{\mathrm{1}}{c}_{\mathrm{5}}{s}_{\mathrm{23}}+(c_{\mathrm{1}}{c}_{\mathrm{23}}{c}_{\mathrm{4}}+s_{\mathrm{1}}{s}_{\mathrm{4}})s_{\mathrm{5}}),\\ p_y=\mathrm{256}{c}_{\mathrm{2}}{s}_{\mathrm{1}}+\mathrm{210}{s}_{\mathrm{1}}{s}_{\mathrm{23}}+\mathrm{144}(-c_{\mathrm{5}}{s}_{\mathrm{1}}{s}_{\mathrm{23}}-(c_{\mathrm{23}}{c}_{\mathrm{4}}{s}_{\mathrm{1}}-c_{\mathrm{1}}{s}_{\mathrm{4}})s_{\mathrm{5}}),\\ p_z=\mathrm{240.5}+\mathrm{256}{s}_{\mathrm{2}}-\mathrm{210}{c}_{\mathrm{23}}+\mathrm{144}(c_{\mathrm{23}}{c}_{\mathrm{5}}-c_{\mathrm{4}}{s}_{\mathrm{23}}{s}_{\mathrm{5}}),\end{array}\right.（6）

式中： 为了简化表示， 令\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_i+{\theta }_j)=c_{ij}， \mathrm{s}\mathrm{i}\mathrm{n}({\theta }_i+{\theta }_j)=s_{ij}。

2.3 三次多项式轨迹规划

机械臂轨迹规划的目的是确保机械臂安全、 正确、 稳定地运行， 这就需要在运行过程中实施关节空间轨迹规划。轨迹规划可以通过调整每个关节角度， 使每个关节变量的位置和速度函数曲线在运行过程中是连续和平滑的。本文使用三次多项式插补， 它是指关节空间之间多节点满足二阶导数连续地插补， 可以保证机器人的起始、 停止速度为零^［21］。

已知初始和终止时刻的关节位置、 速度， 可以使用三阶多项式来确保关节角度和角速度的平滑连续性， 即

\theta (t)=a_{\mathrm{0}}+a_{\mathrm{1}}t+a_{\mathrm{2}}{t}^{\mathrm{2}}+a_{\mathrm{3}}{t}^{\mathrm{3}}\mathrm{。}（7）

初始时刻和终止时刻的关节角度和角速度的4个约束为

\theta (\mathrm{0})={\theta }_{\mathrm{0}}， \theta (t_f)={\theta }_f， \dot{\theta }(\mathrm{0})=\mathrm{0}， \dot{\theta }(t_f)=\mathrm{0},（8）

式中： θ₀为初始时刻的关节角度； θ_f 为终止时刻的关节角度， 初始时刻和终止时刻的关节角速度均为0。通过求解式（7）的一阶和二阶导数， 可以得到角速度函数为

\dot{\theta }(t)=a_{\mathrm{1}}+\mathrm{2}{a}_{\mathrm{2}}t+\mathrm{3}{a}_{\mathrm{3}}{t}^{\mathrm{2}}\mathrm{。}（9）

结合式（7）~（9）， 三阶多项式系数为

\left\{\begin{array}{l}{a}_{\mathrm{0}}={\theta }_{\mathrm{0}},\\ a_{\mathrm{1}}=\mathrm{0},\\ a_{\mathrm{2}}=\frac{\mathrm{3}}{t_f^{\mathrm{2}}}({\theta }_f-{\theta }_{\mathrm{0}}),\\ a_{\mathrm{3}}=-\frac{\mathrm{2}}{t_f^{\mathrm{3}}}({\theta }_f-{\theta }_{\mathrm{0}})\mathrm{。}\end{array}\right.（10）

式（10）可用于计算从任何起始关节角位置到所需结束位置的三次多项式， 但该解仅适用于起始和结束关节角速度均为零的情况。

在机械臂轨迹规划过程中， 使用三次多项式插值定义起始关节角度、 目标关节角度以及关节移动速度为机械臂生成一系列的关节角度， 可以确保机械臂从起始位置平滑地移动到终止位置， 使机械臂在运动过程中避免出现突兀的动作， 减少震动， 提高运动效率和精度。使用绘图函数生成机器人在一定时间段内的轨迹变化过程中关节角度和角速度随时间变化的曲线， 如图 9 所示。

图 9（a） 反映了每个关节随时间变化的运动过程， 图9（b）反映了每个关节角速度随时间的变化， 它们满足三阶多项式插值的起始速度和结束速度均为零的要求。关节角度和角速度曲线平滑， 没有拐点、 中断、 跳跃等现象， 表明机器人运行平稳。

3 测试与结果分析

为了验证艾灸机器人螺旋灸控制方法的可行性， 在轨迹规划过程中收集了真实六自由度机械臂的关节角度误差， 并将其与模拟关节角度进行了比较。此外， 还设计了一个测试平台， 用于评估螺旋灸法实际操作过程中控制方法的准确性。

3.1 可行性分析

当艾灸机器人运行平面螺旋灸法时， 收集了机械臂关节角度的实时数据， 使用最小二乘法拟合轨迹曲线， 并与模拟曲线进行比较， 如图 10（a） 所示。为了更直观地比较， 从Z轴的正方向观察平面螺旋灸轨迹， 可以将其转换为二维平面曲线， 如图 10（b） 所示， 发现实际和模拟的轨迹非常接近。

同时， 处理了各关节角度的数据， 并计算了角度误差， 结果如表 2 所示。研究发现， 关节4和关节6的弧度总是为零， 因此在计算中可以忽略关节4和关节6。计算得出平面螺旋艾灸方法的角度误差保持在4.8×10^-4 rad以下， 在误差允许范围内。

当艾灸机器人执行三维螺旋灸法时， 将实际轨迹与模拟轨迹进行比较， 如图 11（a） 所示。将其转换为二维平面曲线， 从Z轴的正方向观察三维螺旋灸的轨迹， 如图 11（b） 所示。可以看出机械臂末端的整体运动轨迹与模拟轨迹密切相关。

此外， 计算了各关节的弧度误差， 如表 3 所示。结果表明， 三维螺旋灸法的弧度误差保持在4.2×10^-4 rad以下， 在可接受的范围内。

3.2 精确度分析

在实验室环境中构建的测试平台如图 12 所示， 主要包括1个带控制模块的艾灸机器人、 1台戴尔笔记本电脑和1个红外激光指示器。红外激光指示器取代了机械臂末端的灸疗器具， 为观察螺旋灸法的控制方法准确性提供了更直观的方法。

测试实验在桌面上有限的规划区域内进行， 使艾灸机器人的机械臂末端沿着一条标记的螺旋路径进行运动， 以测试它是否能完成计划的螺旋轨迹。螺旋灸法的精确度测试过程如下： 首先， 在机械臂的工作区域内手动选择穴位中心， 由系统记录坐标； 然后， 根据预设的参数启动螺旋灸法； 接着， 通过观察机械臂末端的红外激光指示器来直观地验证螺旋灸法的准确性， 完成一个螺旋轨迹后， 机械臂返回穴位中心， 重复该测试过程； 最后， 经过多次重复试验， 可以计算出本文设计的螺旋灸法的轨迹误差。

图 13 和图 14 为执行平面螺旋灸法和三维螺旋灸法的机械臂运行时的照片， 其中红外激光指示器始终指向标记的螺旋轨迹， 说明本文提出的控制方法可以完成螺旋灸法的整个过程。

在本次路径规划实验中， 设置合适的参数后， 分别对平面和三维螺旋灸法进行了30次测试， 并计算了最大误差， 结果如表 4 所示。

由表 4 可以看出， 该控制方法有效地满足了大多数螺旋灸法应用的精度要求。

在螺旋灸法的测试过程中观察到误差的原因分析如下：

1） 在仿真软件中构建的机械臂模型存在微小的尺寸差异， 包括各关节之间的距离和各个臂长等；

2） 真实的机械臂在调用仿真软件中各个关节角度进行轨迹规划时， 会出现微小的角度误差。

4 结　论

针对人工艾灸存在的不准确、 不安全、 负担重等缺点， 基于六自由度机械臂设计了艾灸机器人， 包括温度监测与反馈、 定时清灰、 净化艾烟等功能， 能够做到精准化、 安全化、 标准化、 简便化施灸。同时， 为了提高艾灸机器人可操作手法的多样性， 根据《中医医疗技术手册》规定的螺旋灸法基本操作方法， 设计了一种能模拟平面螺旋灸法和三维螺旋灸法的运动规划和机械臂控制方法。利用CoppeliaSim对艾灸机器人进行轨迹规划， 使用正运动学分析和三次多项式轨迹规划算法保证了艾灸机器人进行螺旋灸法运动时的精确度。最后， 设计了一个测试平台， 可以直观地评估螺旋灸法实际操作过程中轨迹的准确性。通过多次测试并计算误差， 验证了艾灸机器人实际运行螺旋灸法时的可行性和精确度。与现有的艾灸机器人相比， 丰富了雷火灸的艾灸手法， 对更多的疾病有着良好的治疗效果， 因此在艾灸理疗店和医院有着广阔的应用前景。

参考文献

[1] YAO Z， GUO X. The Application of moxibustion in the concept of “preventive treatment”［J］. MEDS Chinese Medicine， 2022， 4（7）： 6-12.

[2] 孙凯， 赵希瑞， 伏敏睿， 等. 基于现代文献的赵氏雷火灸病谱研究［J］. 湖南中医杂志， 2024， 40（10）： 89-94.

[3] SUN Kai， ZHAO Xirui， FU Minrui， et al. Study on disease spectrum of Zhao’s Leihuo moxibustion based on modern literature［J］. Hunan Journal of Traditional Chinese Medicine， 2024， 40（10）： 89-94. （in Chinese）

[4] YAO Y， ZHOU L， CHEN F Q， et al. The effect and safety of thunder-fire moxibustion for low back pain： a meta-analysis of randomized controlled trials［J］. Evidence-Based Complementary and Alternative Medicine， 2022（1）： 6114417.

[5] 王国强. 中医医疗技术手册（2013普及版）［EB/OL］. ［2015-12-31］.

[6] 马镓钰， 杨志新， 陈建华， 等. 基于红外热成像技术探索“相对穴” 内关穴透灸外关穴温度特征［J］. 针灸临床杂志， 2024， 40（6）： 37-42.

[7] MA Jiayu， YANG Zhixin， CHEN Jianhua， et al. Exploring temperature characteristics of ‘relative acupoints’ during penetration moxibustion of PC6-TE5 based on infrared thermal imaging technology［J］. Journal of Clinical Acupuncture and Moxibustion， 2024， 40（6）： 37-42. （in Chinese）

[8] 王新庆， 徐振兴， 李承蒿， 等. 艾灸机器人控制系统设计及灸法实验研究［J］. 制造业自动化， 2024， 46（11）： 30-37.

[9] WANG Xinqing， XU Zhenxing， LI Chenghao， et al. Control system design of moxibustion robot and experimental study of moxibustion method［J］. Manufacturing Automation， 2024， 46（11）： 30-37. （in Chinese）

[10] 马蓓蓓， 胡志刚， 时鹏， 等. 艾灸机器人系统设计与实现［J］. 计算机工程， 2024， 50（2）： 214-223.

[11] MA Beibei， HU Zhigang， SHI Peng， et al. Design and implementation of moxibustion robot system［J］. Computer Engineering， 2024， 50（2）： 214-223. （in Chinese）

[12] 白善明， 崔娟， 张志东， 等. 基于六自由度机械臂的艾灸机器人手法控制算法［J］. 科学技术与工程， 2023， 23（6）： 2403-2410.

[13] BAI Shanming， CUI Juan， ZHANG Zhidong， et al. Manipulative control algorithm of moxibustion robot based on 6-DOF manipulator［J］. Science Technology and Engineering， 2023， 23（6）： 2403-2410. （in Chinese）

[14] 侍昊， 金鼎森， 俞言骐， 等. 六自由度机械臂若干基本问题的研究［J］. 机械管理开发， 2024， 39（4）： 19-22.

[15] SHI Hao， JIN Dingsen， YU Yanqi， et al. Study on some basic problems of six-degree-of-freedom robotic arm［J］. Mechanical Management and Development， 2024， 39（4）： 19-22. （in Chinese）

[16] KUANG J J， TIAN Z Y， WANG D， et al. Online detection of the volatile intermediates evolved from moxa floss pyrolysis with photoionization mass spectrometry［J］. Journal of Analytical and Applied Pyrolysis， 2022， 168： 105737.

[17] 朱玉权， 张世超， 蔡武. 基于CoppeliaSim的UR5机器人建模与运动学分析［J］. 舰船电子对抗， 2024， 47（5）： 90-94.

[18] ZHU Yuquan， ZHANG Shichao， CAI Wu. Modeling and kinematic analysis of UR5 robot based on CoppeliaSim［J］. Shipboard Electronic Countermeasure， 2024， 47（5）： 90-94. （in Chinese）

[19] 齐伟， 陈雪艳， 迎春， 等. 基于CoppeliaSim的蛇形机器人建模与运动仿真［J］. 计算机仿真， 2024， 41（10）： 443-447.

[20] QI Wei， CHEN Xueyan， YING Chun， et al. Snake robot modeling and motion simulation based on CoppeliaSim［J］. Computer Simulation， 2024， 41（10）： 443-447. （in Chinese）

[21] DRAELOS M. Time-optimal spiral trajectories with closed-form solutions［J］. IEEE Robotics and Automation Letters， 2023， 8（4）： 2213-2220.

[22] 覃洪汉， 黄立仁， 郭志宏， 等. 基于六自由度机械臂的带电作业机械手模型设计研究［J］. 机械设计与制造工程， 2024， 53（9）： 67-72.

[23] QIN Honghan， HUANG Liren， GUO Zhihong， et al. Research on model design of live working manipulator based on six degrees of freedom manipulator［J］. Machine Design and Manufacturing Engineering， 2024， 53（9）： 67-72. （in Chinese）

[24] 陈大光， 刘昱鑫， 刘兴德， 等. 六自由度机械臂运动分析及轨迹规划［J］. 吉林化工学院学报， 2023， 40（1）： 80-86.

[25] CHEN Daguang， LIU Yuxin， LIU Xingde， et al. Motion analysis and trajectory planning of 6-dof manipulator［J］. Journal of Jilin Institute of Chemical Technology， 2023， 40（1）： 80-86. （in Chinese）

[26] 潘明章， 王硕， 李劲， 等. 一种6自由度机械臂的结构设计和运动学分析［J］. 机械传动， 2024， 48（6）： 50-57.

[27] PAN Mingzhang， WANG Shuo， LI Jin， et al. Mechanism design and kinematics analysis of a 6-DOF manipulator［J］. Journal of Mechanical Transmission， 2024， 48（6）： 50-57. （in Chinese）

[28] 王维镝， 汪建松， 刘志刚， 等. 基于宠物猫智能养育箱平台和六自由度机械臂设计与试验［J］. 农业开发与装备， 2024（10）： 137-140.

[29] WANG Weidi， WANG Jiansong， LIU Zhigang， et al. Design and experiment of intelligent incubator platform and six-degree-of-freedom manipulator based on pet cat［J］. Agricultural Development & Equipments， 2024（10）： 137-140. （in Chinese）

[30] 温树仁， 李艳杰. 六自由度机械臂建模与仿真［J］. 上海电气技术， 2024， 17（1）： 67-72.

[31] WEN Shuren， LI Yanjie. Modeling and simulation of 6-DOF manipulator［J］. Journal of Shanghai Electric Technology， 2024， 17（1）： 67-72. （in Chinese）

[32] 孙泽昆， 陆永华， 朱赟， 等. 基于机器视觉的装车臂自动对接系统研究［J］. 测试技术学报， 2024， 38（4）： 427-434.

[33] SUN Zekun， LU Yonghua， ZHU Yun， et al. Research on automatic docking system of LNG land loading arm based on machine vision［J］. Journal of Test and Measurement Technology， 2024， 38（4）： 427-434. （in Chinese）

[34] 陈新秀. 六自由度工业机械臂运动学分析与避障轨迹规划［J］. 自动化应用， 2023（19）： 18-20.

[35] CHEN Xinxiu. Kinematics analysis and obstacle avoidance trajectory planning of six degrees of freedom industrial manipulator［J］. Automation Application， 2023（19）： 18-20. （in Chinese）

[36] 张益峰. 六自由度机械臂轨迹规划研究［J］. 信息技术与信息化， 2023（6）： 79-82.

[37] ZHANG Yifeng. Research on trajectory planning of six-degree-of-freedom manipulator［J］. Information Technology and Informatization， 2023（6）： 79-82. （in Chinese）

基金资助

山西省中医药强省专项(204-110457)

山西省针灸医院智能艾灸机器人项目(HCZB2021110)