新型刚柔混联连续体机器人运动建模

Kinematic Modeling of a Novel Rigid⁃Flexible Hybrid Continuum Robots

摘要

针对一种新型软轴拉扭协同驱动的刚柔混联连续体机器人，开展了其运动建模研究。建立了基于分段恒曲率（PCC）并综合考虑多种载荷的运动静力学模型。为求解其高度非线性的逆运动学，构建了一种由牛顿-拉夫逊优化算法（NRBO）优化的反向传播（BP）神经网络模型（NRBO-BP模型）。实验结果表明，单/双段柔性体机器人末端弯曲角度平均误差分别为4.2°和7.1°；基于NRBO-BP模型的轨迹跟踪最大位置误差为2.5 mm，验证了所提方法的准确性与有效性。

Abstract

A study of the kinemetic modeling was carried out for a novel rigid-flexible hybrid continuum robot driven by tension-torsion synergistic actuation. A kinetostatic model was developed based on the piecewise constant curvature（PCC） framework， which comprehensively considered various loads. To solve the highly nonlinear inverse kinematics， a BP neural network model optimized by the Newton-Raphson based optimizer（NRBO）， denoted as NRBO-BP model， was constructed. Experimental results show that the average bending angle errors at the end of the single/dual-segment flexible robots are as 4.2° and 7.1°， respectively. The maximum position error in trajectory tracking based on NRBO-BP model is as 2.5 mm， which verifies the accuracy and effectiveness of the proposed methods.

Graphical abstract

Author summay

董佳祥，男，1997年生，博士研究生。研究方向为微创手术机器人创新设计与精准控制。E-mail： [e-mail]

刘铨权^*（通信作者），男，1984年生，副教授。研究方向为智能机器人与系统集成、医疗机器人机构综合与控制。E-mail： [e-mail]。

0 引言

近年来，受生物结构启发的连续体机器人凭借其高顺应性与灵巧度为微创手术等场景提供了新技术路径。然而，连续体机器人的精确建模面临三大挑战：①无限自由度的超冗余特性使传统D⁃H建模不再适用^［1］；②静力学与运动学的强耦合；③柔性材料的高度非线性、大变形及迟滞效应^［2］。

为应对上述挑战，研究人员发展了多种建模策略，特别是分段恒曲率（PCC）模型^［1］，因其计算高效而在实时控制中得到广泛应用。OUYANG等^［3］基于PCC假设，设计并分析了一款用于微创手术的三段式连续体机器人。MISHRA等^［4］采用PCC假设和D-H参数法对一种模块化连续体机器人进行运动学建模。为提高模型的保真度，YUAN等^［5］提出了一种考虑了摩擦、重力及外部负载的综合静态模型，以提高PCC框架下的模型描述精度。本文所研究的新型刚柔混联机器人具有独特的刚柔耦合结构与软轴驱动方式，因此，笔者在PCC框架下为其建立了一个综合考虑多种载荷的运动静力学模型。

与此同时，对于连续体机器人高度非线性的逆运动学求解，纯解析方法往往难以胜任，为此，基于数据驱动的代理模型为直接学习输入到输出的复杂映射提供了思路。YANG等^［6］提出了一种基于粒子群优化-反向传播（PSO-BP）神经网络的端到端建模方法，直接学习末端坐标到驱动气压的逆向映射。SCHÄFKE等^［7］利用循环神经网络（RNN）预测了一个五自由度关节式软体机器人的动态行为。此外，为规避传统机器学习方法中繁琐的数据采集和离线训练过程，研究人员提出在线学习方法。TAN等^［8］通过两个连续时间归零神经网络分别在线生成驱动信号，实时估计未知机器人雅可比矩阵，仅依赖当前控制输入和实时传感器数据实现对连续体机器人的运动学控制。然而，上述数据驱动方法普遍依赖于大规模、高质量数据集，通过物理实验进行数据采集往往耗时耗力，且难以覆盖整个工作空间。

综上，本文针对软轴拉扭共轴协同驱动的新型刚柔混联连续体机器人，首先采用PCC和D⁃H参数法构建其运动静力学模型；然后，为高效获取训练数据，利用ADAMS多体动力学仿真环境构建数据集；最后，在此数据基础上提出基于牛顿⁃拉夫逊优化算法（NRBO）^［9］优化的反向传播（BP）神经网络模型（简称NRBO⁃BP模型）的逆运动学求解策略，并通过实验对所提建模方法的准确性与有效性进行验证。

1 机器人机构概述与运动学分析

1.1 机器人机构概述

本研究建模与控制对象为一种新型刚柔混联连续体机器人，其系统构成如图1所示。该机器人创新性地利用了软轴的拉扭协同驱动方式，其机械手本体由近端柔性段、远端柔性段及末端刚性关节串联构成。该系统共配置6根软轴，通过其牵拉运动控制双柔性段的弯曲形态，利用其扭转运动驱动末端刚性关节，共实现7个操作自由度（DOF 2~8）。此外，整个后置驱动单元被安装于一个独立的平移模块上，为机械手提供第8个自由度（DOF 1），即整体的轴向往返运动。

1.2 机器人运动学分析

本节采用PCC与D⁃H参数法建立机械手的正向运动学模型。此外，考虑到系统高度非线性特性，为求解机械手的逆运动学，构建了NRBO-BP模型，以建立末端坐标与软轴驱动量之间的映射关系。

首先，假设柔性连续体的弹性骨架由一系列具有恒定曲率相互正切的圆弧单元组成，且这一系列圆弧单元均发生平面变形。分隔盘过线孔坐标见图2a，其中，1、2和3代表远端段布线路径，4、5和6代表近端段布线路径，且该序号与电机轴ID一一对应。机械手局部坐标系下过线孔坐标可以表示为

\left.\begin{array}{l}{\sigma }_{\mathrm{1}i}=-\mathrm{\pi }/\mathrm{6}+\left(i-\mathrm{1}\right)\beta \\ {\sigma }_{\mathrm{2}i}=\mathrm{\pi }/\mathrm{6}-\left(i-\mathrm{1}\right)\beta \end{array}\right\}（1）

i=\mathrm{1},\mathrm{2},\mathrm{3};\beta =\mathrm{2}\pi /\mathrm{3}

式中：{\sigma }_{\mathrm{1}i}、{\sigma }_{\mathrm{2}i}分别为近、远端段柔性关节中第i根软轴的角度坐标。

本文采用分段恒曲率方法将连续的柔性段离散为一系列虚拟的节。其中，任意两个相邻的分隔盘及其之间的弹性骨架构成一节柔性关节。柔性关节坐标系定义如图2b所示，局部相切的每节柔性关节亦可采用该种定义方式。通过({\kappa }_i,{\varphi }_i)参数化弹性骨架，其中，{\kappa }_i为第i节弯曲曲率，{\varphi }_i为第i节弯曲平面与参考轴x的夹角（即方位角）。采用弧几何法描述机械手位姿，第i节弧上任意一点在s处关于该节弧局部坐标系的位置可描述为

{\mathit{P}}_s=\left[\begin{array}{c}\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_i}{{\kappa }_i}(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}({\kappa }_{i}s))\\ \frac{\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_i}{{\kappa }_i}(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}({\kappa }_{i}s))\\ \frac{\mathrm{1}}{{\kappa }_i}\mathrm{s}\mathrm{i}\mathrm{n}({\kappa }_{i}s)\end{array}\right]（2）

式中：({\kappa }_{i}s)表示第i节柔性关节绕y轴旋转，简化为{\vartheta }_i。

弯曲曲率{\kappa }_i和方位角{\varphi }_i可描述为

\left.\begin{array}{l}{\kappa }_i=\sqrt[]{{\kappa }_{ix}^{\mathrm{2}}+{\kappa }_{iy}^{\mathrm{2}}}\\ {\varphi }_i=\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{2}({\kappa }_{iy},{\kappa }_{ix})\end{array}\right\}（3）

式中：{\kappa }_{ix}、{\kappa }_{iy}分别为第i节柔性关节局部坐标系下沿x、y轴的弯曲曲率分量。

通过ZYZ欧拉角组方式，考虑弹性骨架发生扭转情形，则第i节坐标系相对于第i-1节坐标系的齐次变换矩阵可以描述为

{}_i{}^{i-\mathrm{1}}\mathit{T}(s)=\left[\begin{array}{cc}{\mathit{R}}_z({\varphi }_i){\mathit{R}}_y({\vartheta }_i){\mathit{R}}_z({\epsilon }_i-{\varphi }_i)& {\mathit{P}}_s\\ \mathrm{0}& \mathrm{1}\end{array}\right]（4）

式中：{\epsilon }_i为第i节柔性关节绕z轴额外旋转角；{\mathit{R}}_z、{\mathit{R}}_y分别为绕z\mathrm{、}y轴的旋转矩阵。

由于两段柔性关节共采用10个分隔盘，此时，将远端段柔性关节末端分隔盘的局部坐标系定义为x_{\mathrm{10}}{y}_{\mathrm{10}}{z}_{\mathrm{10}}，刚性关节坐标系定义如图3所示，其中关节特征点定义为M₁、M₂、M₃和M₄，m_{\mathrm{0}}\mathrm{、}{m}_{\mathrm{1}}\mathrm{、}{m}_{\mathrm{2}}为连杆偏距，{\theta }_{\mathrm{R}\mathrm{1}}、{\theta }_{\mathrm{R}\mathrm{2}}为关节角。改进D-H参数如表1所示，其中，{\alpha }_{i-\mathrm{1}}为连杆扭转角（即关节轴z_{i-\mathrm{1}}与z_i之间的夹角），a_{i-\mathrm{1}}为连杆长度（即沿公法线方向，两相邻关节轴z_{i-\mathrm{1}}与z_i之间的距离），d_i为连杆偏距（即沿关节轴z_i方向，两相邻公法线之间的距离），{\theta }_i为关节角（即绕关节轴z_i，两相邻连杆之间的转动角度）。

在刚性、柔性关节连接处通过螺钉固连方式实现了刚柔耦合处位置和姿态的连续性，则刚性关节夹钳处相对于机械手基坐标系的齐次变换矩阵可以描述为

{}_{\mathrm{14}}{}^{\mathrm{0}}\mathit{T}={}_{\mathrm{10}}{}^{\mathrm{0}}\mathit{T}{}_{\mathrm{11}}{}^{\mathrm{10}}\mathit{T}{}_{\mathrm{13}}{}^{\mathrm{11}}\mathit{T}{}_{\mathrm{14}}{}^{\mathrm{13}}\mathit{T}={}_{\mathrm{10}}{}^{\mathrm{0}}\mathit{T}\left[\begin{array}{cccc}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& m_{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\cdot\left[\begin{array}{cccc}{c}_{\mathrm{1}}{c}_{\mathrm{2}}& -c_{\mathrm{1}}{s}_{\mathrm{2}}& s_{\mathrm{1}}& m_{\mathrm{1}}{s}_{\mathrm{1}}\\ s_{\mathrm{2}}& c_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ -c_{\mathrm{2}}{s}_{\mathrm{1}}& s_{\mathrm{1}}{s}_{\mathrm{2}}& c_{\mathrm{1}}& m_{\mathrm{1}}{c}_{\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& m_{\mathrm{0}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（5）

c_{\mathrm{1}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{R}\mathrm{1}}　 s_{\mathrm{1}}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{R}\mathrm{1}}　 c_{\mathrm{2}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{R}\mathrm{2}}　 s_{\mathrm{2}}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{R}\mathrm{2}}

软轴拉扭共轴协同驱动的新型刚柔混联连续体机器人坐标系定义如图4所示，其中，刚性关节处绿、红、黑色线段分别与连杆偏距m₀、m₁、m₂对应。

为求解系统高度非线性的逆运动学，本研究构建了NRBO-BP模型。该模型旨在学习从机器人末端空间坐标到6根软轴牵拉位移的映射关系。为训练该模型，利用ADAMS仿真环境生成了覆盖机器人完整工作空间的3000组数据点，其分布如图5所示。基于BP神经网络的软轴驱动预测末端位置原理如图6所示。

输入(X_{k\mathrm{1}},Y_{k\mathrm{2}},Z_{k\mathrm{3}},k=\mathrm{1},\mathrm{2},\cdots ,\mathrm{3000})为机械手末端位置数据，输出(L_{k\mathrm{1}},L_{k\mathrm{2}},L_{k\mathrm{3}},L_{k\mathrm{4}},L_{k\mathrm{5}},L_{k\mathrm{6}},k=\mathrm{1,2},\dots ,\mathrm{3000})分别为6根软轴的位移驱动量，其中，k为数据点序号。设输入向量为\mathit{X}=(X_{k\mathrm{1}},Y_{k\mathrm{2}},Z_{k\mathrm{3}}{)}^{\mathrm{T}}，输入层到隐藏层的线性组合描述为

{\mathit{Z}}^{(\mathrm{1})}={\mathit{W}}^{(\mathrm{1})}\mathit{X}+{\mathit{b}}^{(\mathrm{1})}（6）

式中：{\mathit{W}}^{(\mathrm{1})}为输入层到隐藏层的权重；{\mathit{b}}^{(\mathrm{1})}为隐藏层的偏置。

采用tanh激活函数，隐藏层输出可以描述为

{\mathit{A}}^{(\mathrm{1})}=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}({\mathit{Z}}^{(\mathrm{1})})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{Z}}^{(\mathrm{1})})-\mathrm{e}\mathrm{x}\mathrm{p}(-{\mathit{Z}}^{(\mathrm{1})})}{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{Z}}^{(\mathrm{1})})+\mathrm{e}\mathrm{x}\mathrm{p}(-{\mathit{Z}}^{(\mathrm{1})})}（7）

隐藏层到输出层的线性组合为

{\mathit{Z}}^{(\mathrm{2})}={\mathit{W}}^{(\mathrm{2})}{\mathit{A}}^{(\mathrm{1})}+{\mathit{b}}^{(\mathrm{2})}（8）

式中：{\mathit{W}}^{(\mathrm{2})}为隐藏层到输出层的权重；{\mathit{b}}^{(\mathrm{2})}为输出层的偏置。

通过对输出层应用线性激活函数purelin，此时，输入到输出的映射关系表示为

\left.\begin{array}{l}{L}_{k\mathrm{1}}={\displaystyle \sum _{n=\mathrm{1}}^{N_h}}{W}_{\mathrm{1}n}^{\left(\mathrm{2}\right)}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}({\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{3}}}{W}_{nm}^{(\mathrm{1})}{X}_k+b_n^{(\mathrm{1})})+b_{\mathrm{1}}^{(\mathrm{2})}\\ L_{k\mathrm{2}}={\displaystyle \sum _{n=\mathrm{1}}^{N_h}}{W}_{\mathrm{2}n}^{\left(\mathrm{2}\right)}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}({\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{3}}}{W}_{nm}^{(\mathrm{1})}{X}_k+b_n^{(\mathrm{1})})+b_{\mathrm{2}}^{(\mathrm{2})}\\ ⋮\\ L_{k\mathrm{6}}={\displaystyle \sum _{n=\mathrm{1}}^{N_h}}{W}_{\mathrm{6}n}^{\left(\mathrm{2}\right)}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}({\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{3}}}{W}_{nm}^{(\mathrm{1})}{X}_k+b_n^{(\mathrm{1})})+b_{\mathrm{6}}^{(\mathrm{2})}\end{array}\right\}（9）

式中：N_h为隐藏层节点个数；n为隐藏层节点的索引；m、p（p=1，2，…，6）分别为输入层和输出层节点的索引；W_{nm}^{(\mathrm{1})}为连接输入层第m个节点与隐藏层第n个节点的权重；W_{pn}^{(\mathrm{2})}为连接隐藏层第n个节点与输出层第p个节点的权重；b_n^{(\mathrm{1})}为隐藏层第n个节点的偏置值；b_p^{(\mathrm{2})}为输出层第p个节点的偏置值。

损失函数F_loss是模型预测与真值之差的量化指标，本文将其定义为

F_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}N}{\displaystyle \sum _{k=\mathrm{1}}^N}(Y_k-{\widehat{Y}}_k{)}^{\mathrm{2}}（10）

式中：Y_k为第k个样本的真值；{\widehat{Y}}_k为第k个样本的预测值；N为样本总数。

为克服传统梯度下降法在优化BP网络参数\mathit{\Theta }=({\mathit{W}}^{(\mathrm{1})},{\mathit{b}}^{(\mathrm{1})},{\mathit{W}}^{(\mathrm{2})},{\mathit{b}}^{(\mathrm{2})})时易陷入局部最优的缺陷，引入新型元启发算法NRBO。该算法将牛顿法的二阶收敛性与群体智能结合，优化 BP 神经网络的权重和偏置参数\mathit{\Theta }，其中解空间\Theta \in {\mathbf{R}}^D（D为网络参数维度），其差分扰动项可以写为

\mathit{\rho }=r_{\mathrm{1}}({\mathit{P}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-{\mathit{X}}_k)+r_{\mathrm{2}}({\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{1}}-{\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{2}})（11）

式中：r₁、r₂为随机数； {\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{1}}、{\mathit{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{2}}为从种群中随机选择的两个不同个体；{\mathit{X}}_k为个体当前位置；{\mathit{P}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}为当前全局最优解。

由此通过多部构造形成基础更新位置{\mathit{X}}_{\mathrm{u}\mathrm{p}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{e}}。为防止算法过早收敛于局部最优，本文设置预设概率P_{\mathrm{T}\mathrm{A}\mathrm{O}}=\mathrm{0.6}以激活陷阱规避算子，该算子的扰动形式可以概括为

{\mathit{X}}_{\mathrm{T}\mathrm{A}\mathrm{O}}={\mathit{X}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}+{\theta }_{\mathrm{1}}\left(u_{\mathrm{1}}{\mathit{P}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-u_{\mathrm{2}}{\mathit{X}}_k\right)+

{\theta }_{\mathrm{2}}\delta (u_{\mathrm{1}}\overline{\mathit{P}}-u_{\mathrm{2}}{\mathit{X}}_k)（12）

式中：{\mathit{X}}_{\mathrm{T}\mathrm{A}\mathrm{O}}为经过陷阱规避算子更新后的新位置向量；{\mathit{X}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}为扰动的基础向量；\delta为动态衰减因子；{\theta }_{\mathrm{1}}、{\theta }_{\mathrm{2}}为随机方向系数；u_{\mathrm{1}}、u_{\mathrm{2}}为随机加权系数；\overline{\mathit{P}}为所有个体的平均位置向量。

整体训练过程如下：

1） 初始化。将BP网络参数\mathit{\Theta }作为优化变量。在预设的参数边界内随机生成初始种群\{{\mathit{X}}_i{\}}_{i=\mathrm{1}}^N，并基于式（10）评估适应度。

2） 迭代优化。达到最大迭代次数前，对种群中每个个体{\mathit{X}}_k：①计算动态因子和差分扰动向量；②应用牛顿-拉夫逊搜索规则生成基础更新位置{\mathit{X}}_{\mathrm{u}\mathrm{p}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{e}}；③根据预设概率P_{\mathrm{T}\mathrm{A}\mathrm{O}}决定是否执行陷阱规避算子，得到最终候选解{\mathit{X}}_{\mathrm{n}\mathrm{e}\mathrm{w}}；④对{\mathit{X}}_{\mathrm{n}\mathrm{e}\mathrm{w}}进行边界处理，并评估其适应度；⑤根据适应度更新个体最优和全局最优解{\mathit{P}}_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}（即最优网络参数{\mathit{\Theta }}^{\mathrm{*}}）。

3） 输出。返回优化后的网络参数{\mathit{\Theta }}^{\mathrm{*}}，构建NRBO-BP模型。

2 刚柔混联连续体机器人驱动力-关节运动模型构建

刚柔混联连续体机器人由一系列相同相邻单元构成，分析某段单元静力学时，通过迭代方式即可建立整只机械手的静力学模型，如图7所示。为与机械手物理空间中坐标系一致，则重力方向沿-y方向。本节将从软轴上的驱动力、与外部环境的交互力、重力以及弹性力这4个方面的综合因素进行分析。

结合式（1）所示的过线孔坐标定义，软轴作用在分隔盘上的驱动力可以表示为

（13）

式中：{}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{f}\mathrm{s},i}为i-1节坐标系下软轴对第i个分隔盘所施加的驱动力；{}^{i-\mathrm{1}}{\mathit{P}}_{\mathrm{f}\mathrm{s},i+\mathrm{1}}为第i+1个分隔盘上的布线坐标在i-1节坐标系下的位置向量；{\mathit{F}}_{\mathrm{f}\mathrm{s}\_i}为沿软轴第i节方向上的驱动力；N_{\mathrm{s}\mathrm{d}}为分隔盘总数。

考虑到分隔盘上布线孔径处粉末冶金含油轴承的良好润滑特性，假设忽略软轴与分隔盘之间的摩擦，则垂直于分隔盘的力为0，当i\ne N_{\mathrm{s}\mathrm{d}}或i\ne N_{\mathrm{s}\mathrm{d}}/\mathrm{2}时，i-1节坐标系下软轴对第i个分隔盘所施加的驱动力为

{}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{f}\mathrm{s},i}{=}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{f}\mathrm{s},i}-({(}^{i-\mathrm{1}}{\mathit{n}}_i{)}^{\mathrm{T}}{·}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{f}\mathrm{s},i}){\cdot }^{i-\mathrm{1}}{\mathit{n}}_i（14）

式中：{}^{i-\mathrm{1}}{\mathit{n}}_i为第i个分隔盘在i-1节坐标系下的单位法向量，可通过齐次变换矩阵{}_i{}^{i-\mathrm{1}}\mathit{T}的第三列获得。

软轴对分隔盘的驱动力所产生的力矩描述为

{}^{i-\mathrm{1}}{\mathit{M}}_{\mathrm{f}\mathrm{s},i}{=}^{i-\mathrm{1}}{\mathit{P}}_{\mathrm{f}\mathrm{s},i}{\times }^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{f}\mathrm{s},i}（15）

j=N_{\mathrm{s}\mathrm{d}}条件下，机械手末端盘处与外部环境的交互力可描述为

{}^{\mathrm{0}}{\mathit{F}}_{\mathrm{e}\mathrm{x}}={}_j^{\mathrm{0}}\mathit{T}\cdot {{}_{}{}^j\mathit{F}}_{\mathrm{e}\mathrm{x}}（16）

式中：{}^{\mathrm{0}}{\mathit{F}}_{\mathrm{e}\mathrm{x}}、{}^j{\mathit{F}}_{\mathrm{e}\mathrm{x}}分别为基坐标系、j坐标系下的外部负载力；{}_j{}^{\mathrm{0}}\mathit{T}为j坐标系相对于基坐标系的齐次变换矩阵。

j=N_{\mathrm{s}\mathrm{d}}条件下考虑到外部力矩，j坐标系下的力矩可描述为

{}^j{\mathit{M}}_{\mathrm{e}\mathrm{x}}={(}_j^{\mathrm{0}}{\mathit{T})}^{-\mathrm{1}}{\cdot }^{\mathrm{0}}{\mathit{M}}_{\mathrm{e}\mathrm{x}}+({}_j{}^{j-\mathrm{1}}\mathit{T}{[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}]}^{\mathrm{T}}){\times }^j{\mathit{F}}_{\mathrm{e}\mathrm{x}}（17）

式中：{}^{\mathrm{0}}{\mathit{M}}_{\mathrm{e}\mathrm{x}}为基坐标系下施加于尖端的力矩向量。

本文讨论的重力包含分隔盘和弹性骨架两部分重力，其中刚性关节重力集中在机械手的末端分隔盘上。在i-1节坐标系下，分隔盘和弹性骨架重力可以描述为

\left.\begin{array}{l}{}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{s}\mathrm{d}\_\mathrm{g}}=m_{\mathrm{s}\mathrm{d}}g\mathit{A}\\ {}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{b}\mathrm{b}\_\mathrm{g}}=m_{\mathrm{b}\mathrm{b}}g\mathit{A}\end{array}\right\}（18）

\mathit{A}=({}_{i-\mathrm{1}}{}^{\mathrm{0}}\mathit{T}{)}^{-\mathrm{1}}{\left[\begin{array}{cccc}\mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}

式中：{}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{s}\mathrm{d}\_\mathrm{g}}为分隔盘重力；m_sd为分隔盘质量；{}^{i-\mathrm{1}}{\mathit{F}}_{\mathrm{b}\mathrm{b}\_\mathrm{g}}为弹性骨架重力；m_bb为弹性骨架质量；g为重力加速度；\mathit{A}为i-1节坐标系下的单位重力列向量。

针对弹性骨架部分，第i节弹性骨架在i-1节坐标系下的质心位置为

{}^{i-\mathrm{1}}{\mathit{P}}_{\mathrm{c}\mathrm{f}\mathrm{s},i}=\left[\begin{array}{c}(\mathrm{1}-\frac{\mathrm{s}\mathrm{i}\mathrm{n}({\vartheta }_i/\mathrm{2})}{{\vartheta }_i/\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}({\vartheta }_i/\mathrm{2}))\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_i\\ (\mathrm{1}-\frac{\mathrm{s}\mathrm{i}\mathrm{n}({\vartheta }_i/\mathrm{2})}{{\vartheta }_i/\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}({\vartheta }_i/\mathrm{2}))\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_i\\ \frac{\mathrm{s}\mathrm{i}\mathrm{n}({\vartheta }_i/\mathrm{2})}{{\vartheta }_i/\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}({\vartheta }_i/\mathrm{2})\end{array}\right]（19）

式中：\frac{\mathrm{s}\mathrm{i}\mathrm{n}({\vartheta }_i/\mathrm{2})}{{\vartheta }_i/\mathrm{2}}为第i节弹性骨架质心位置修正项。

根据分段恒曲率假设，沿弹性骨架方向的弯矩可以简化为其曲率的一维方程^［10］，则第i节弹性骨架（含软轴）的弯矩为

{\mathit{M}}_{\mathrm{b},i}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_i& -\mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_i& \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\varphi }_i& \mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_i& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}\mathrm{0}\\ {\kappa }_{i}EI\\ \mathrm{0}\end{array}\right]（20）

式中：E为弹性骨架（含软轴）的弹性模量；I为弹性骨架（含软轴）的惯性矩。

由于存在绕z轴的{\epsilon }_i额外旋转角{\epsilon }_i，则第i节弹性骨架（含软轴）的扭矩为

{\mathit{M}}_{\mathrm{t},i}={}_i{}^{i-\mathrm{1}}\mathit{T}{\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \frac{\mathrm{2}IG{\epsilon }_i}{l_i}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}（21）

式中：G为弹性骨架（含软轴）的剪切模量；l_i为第i个分隔盘与第i-1个分隔盘的垂直距离。

综上，静力学平衡方程为

\mathit{E}=\left[\begin{array}{c}{\mathit{M}}_{\mathrm{b},\mathrm{1}}+{\mathit{M}}_{\mathrm{t},\mathrm{1}}-{\displaystyle \sum _{i=\mathrm{1}}^{N_{sd}}}{}_{}{}^{i-\mathrm{1}}{\mathit{M}}_{net,i}\\ ⋮\\ \begin{array}{l}{\mathit{M}}_{\mathrm{b},N_{\mathrm{s}\mathrm{d}}-\mathrm{1}}+{\mathit{M}}_{\mathrm{t},N_{\mathrm{s}\mathrm{d}}-\mathrm{1}}-({}_{}{}^{N_{\mathrm{s}\mathrm{d}}-\mathrm{1}}\mathit{M}{\mathit{M}}_{\mathrm{n}\mathrm{e}\mathrm{t},N_{\mathrm{s}\mathrm{d}}}+{}_{}{}^{N_{\mathrm{s}\mathrm{d}}-\mathrm{2}}\mathit{M}{\mathit{M}}_{\mathrm{n}\mathrm{e}\mathrm{t},N_{\mathrm{s}\mathrm{d}}-\mathrm{1}})\\ {\mathit{M}}_{\mathrm{b},N_{\mathrm{s}\mathrm{d}}}+{\mathit{M}}_{\mathrm{t},N_{\mathrm{s}\mathrm{d}}}-{}_{}{}^{N_{\mathrm{s}\mathrm{d}}-\mathrm{1}}\mathit{M}{\mathit{M}}_{\mathrm{n}\mathrm{e}\mathrm{t},N_{\mathrm{s}\mathrm{d}}}\end{array}\end{array}\right]（22）

式中：{}^{i-\mathrm{1}}{\mathit{M}}_{\mathrm{n}\mathrm{e}\mathrm{t},i}为第i节柔性关节在i-\mathrm{1}节坐标系下的净外部力矩。

通过fsolve函数求解该非线性方程，找到满足该方程组的解，即获得分段恒曲率假设下每一节对应的{\kappa }_{ix}、{\kappa }_{iy}以及弹性骨架的额外旋转角{\epsilon }_i。

3 实验验证

3.1 实验平台概述

为验证所建模型与算法，搭建了图8所示的实验平台。该平台核心是感知系统：1套空间定位精度为±0.2 mm的光学动作捕捉系统（北京度量， NOKOV MARS1.3H）用于获取末端位姿真值；6个拉压力传感器（深圳锐力德， RDF-TM10）串联于驱动链中，用于测量软轴输入牵拉力。所有实验均以机器人因重力产生的自然下垂状态为初始构型。

3.2 基于静力学模型的运动精度实验

本实验基于静力学分析模型，通过对机器人施加不同软轴牵拉力组合以测试其弯曲精度。实验过程中，将主驱动软轴以恒定速率牵拉至目标力值，而非主驱动软轴则通过松弛方式跟随运动，同时保持恒定的2 N预紧力。本文选取了12种典型的单段及双段软轴驱动配置进行测试。具体而言，对于单段机器人，配置涵盖了自然下垂状态、单软轴驱动（仅驱动F₄或F₆）在10、20、30 N三个不同力级下的弯曲，以及双软轴协同驱动（同时驱动F₄和F₅）在不同力值组合下的弯曲。对于双段机器人，配置涵盖了驱动近端、远端以及两段同时运动的组合工况，涉及两根、三根软轴的协同驱动，且最大驱动力达到了50 N。图9给出一种单段驱动配置（F₄=30 N， F₅=F₆=2 N）下的实时力数据。可以看出，当主驱动软轴4的拉力从0逐渐增大至30 N时，软轴5和软轴6上的拉力逐渐到达预设预紧力2 N。

对全部12种配置完成加载并记录稳定状态后，将实验测得的弯曲角度与静力学模型的仿真预测值进行对比，如图10所示，其中，横坐标代表12种不同的驱动配置。可以看出，单、双段样机仿真数据与实测数据具有良好的一致性。经计算，单段样机在不同配置下弯曲角度的平均误差为4.2°，双段样机在不同配置下弯曲角度的平均误差为7.1°。

3.3 基于NRBO-BP神经网络的位置控制运动精度实验

为评估NRBO-BP神经网络模型的位置控制精度，在机械手任务空间内设计了一个正方形轨迹，其边长为20 mm。首先，在正方形轨迹上均匀地选取12个规划点，将其空间位置数据作为输入传递到NRBO-BP神经网络模型中，从而得到6根软轴的牵拉位移。然后，将获得的拉伸量作为ADAMS平台输入数据，得到机器人运动的仿真轨迹，同时将6根软轴的牵拉拉伸量输入到连续体机器人样机控制系统中得到实测数据。实测轨迹、NRBO-BP神经网络预测轨迹与规划轨迹的对比如图11所示，可以看出，三者具有良好的吻合度。为进一步量化规划轨迹点误差，绘制x、y、z方向上的轨迹点误差，如图12所示。在x、y、z方向上，实测轨迹的最大误差分别为1.8、2.5、2.2 mm，而经过NRBO-BP神经网络预测的轨迹，其最大误差分别降至0.7、2.0、1.9 mm。

4 结论

1） 考虑机器人刚柔耦合特性，综合软轴上的驱动力、与外部环境的交互力、重力及弹性力四个因素，建立了刚柔混联连续体机器人静力学模型，由此构建了驱动力-关节运动模型，单、双段柔性体的弯曲角度平均误差分别为4.2°和7.1°。

2） 针对柔性关节非线性及分布式驱动特点，基于牛顿⁃拉夫逊优化算法（NRBO）优化的反向传播（BP）神经网络构建了从末端空间坐标到软轴驱动量的数据驱动映射模型。正方形轨迹跟踪实验验证了其运动学精度，最大位置误差为2.5 mm。

参考文献

[1] WEBSTER R J， JONES B A. Design and Kinematic Modeling of Constant Curvature Continuum Robots： a Review［J］. International Journal of Robotics Research， 2010， 29（13）： 1661-1683.

[2] 李法民， 郑天江， 沈雯钧， 等. 绳驱动连续体机器人标定方法［J］. 中国机械工程， 2022， 33（2）： 202-208.

[3] LI Famin， ZHENG Tianjiang， SHEN Wenjun， et al. Calibration Method for Cable-driven Conti⁃nuum Robot［J］. China Mechanical Engineering， 2022， 33（2）： 202-208.

[4] OUYANG B， LIU Y， SUN D. Design of a Three-Segment Continuum Robot for Minimally Invasive Surgery［J］. Robotics and Biomimetics， 2016， 3（1）： 2.

[5] MISHRA A， MONDINI A， DEL DOTTORE E， et al. Modular Continuum Manipulator： Analysis and Characterization of Its Basic Module［J］. Biomimetics， 2018， 3（1）： 3.

[6] YUAN H， ZHOU L， XU W. A Comprehensive Static Model of Cable-driven Multi-section Continuum Robots Considering Friction Effect［J］. Mechanism and Machine Theory， 2019， 135： 130-149.

[7] YANG Y， GENG P， HE Y， et al. Inverse Kinematics Modeling of Soft Manipulator Based on Optimized Particle Swarm Optimization-Back Propagation Neural Network［J］. IEEE/ASME Transactions on Mechatronics， 2025，30：1-10.

[8] SCHÄFKE H， LHABICH T， MUHMANN C， et al. Learning-based Nonlinear Model Predictive Control of Articulated Soft Robots Using Recurrent Neural Networks［J］. IEEE Robotics and Automation Letters， 2024， 9（12）： 11609-11616.

[9] TAN N， YU P， ZHONG Z， et al. Data-driven Control for Continuum Robots Based on Discrete Zeroing Neural Networks［J］. IEEE Transactions on Industrial Informatics， 2023， 19（5）： 7088-7098.

[10] SOWMYA R， PREMKUMAR M， JANGIR P. Newton-Raphson-based Optimizer：a New Population-based Metaheuristic Algorithm for Continuous Optimization Problems［J］. Engineering Applications of Artificial Intelligence， 2024， 128： 107532.

[11] RAO P， PEYRON Q， LILGE S， et al. How to Model Tendon-driven Continuum Robots and Benchmark Modelling Performance［J］. Frontiers in Robotics and AI， 2021， 7： 630245.