机器人辅助微创外科手术作为目前研究的热点领域，随着机器人技术的发展也变得更加成熟。机器人辅助微创外科手术相比较于传统的微创外科手术具有独特的优势，例如，从根本上改变了医生的手术方法和手术环境，有效减少了操作医生的疲劳和手术中的误操作，从而大大提高了微创外科手术的效果^[1]。现在微创外科手术机器人的结构已经逐渐成熟，其中较为有名的商业化的手术机器人系统有AESOP系统^[2]、Zeus系统^[3]和Da Vinci手术系统^[4]等。

由于机械臂的尺寸和结构对自身的运动性能影响很大，因此如何通过优化几何参数和拓扑结构得到较优的运动性能在机械臂的设计过程中具有重要的意义。邓宗全等^[5]基于变密度法对摇臂悬架进行了结构拓扑优化设计。贾世元等^[6]提出了姿态可操作度的指标，并利用该指标对机械臂的尺寸进行了优化。丁渊明等^[7]提出一种基于工作空间和能量消耗的综合指标对机械臂的结构进行了优化。岳龙旺等^[8]将条件数的倒数作为灵活度指标实现对双四连杆机构的参数优化。杨世强等^[9]基于对模型的有限元分析，对关键构件进行了减重设计和结构优化。杜志江等^[10]利用条件数和可操作度建立了综合性能指标，并基于该指标实现了对远心点调整机构的杆件参数的优化，王宏民等^[11]基于全域性能指标对主操作手结构进行了优化。赵凯等^[12]以可操作度为性能指标，对机械臂被动关节参数进行了优化，提高了整体的灵活性。Surbhi Gupta等^[13]基于手术中可能的切口位置数量和空间可达性，实现了对三自由度串联手术机械手杆长参数的优化。

结合实际手术机械臂的运动要求，以机械臂运动的灵活性和平稳性为首要要求，建立能够反映该要求的评价指标。通过对手术机械臂远心点调整机构运动学性能的分析，建立了基于全域空间内的条件数均值和条件数波动值的手术机械臂的综合灵巧度评价指标，实现对手术机械臂远心点调整机构的尺寸优化。基于对手术机械臂的有限元分析，对被动关节进行了结构优化。

1 通用手术机械臂结构

优化的微创外科手术通用机械臂如图 1所示，分为主动关节与被动关节，主动关节主要用于手术过程中调整腹腔镜或者微器械末端的位姿，被动关节主要用于术前调整远心点位姿。机械臂末端可以安装腹腔镜或者微器械，对应的成为持镜臂和持械臂。一般的手术系统由一条持镜臂和两条持械臂组成，如图 2所示。在手术过程中，持镜臂实时为医生提供病人腹腔内的手术视野，持械臂通过医生控制的主手实现手术过程中对病灶区域的操作。

双平行四边形机构可以从机械角度保证手术过程中远心点的固定位置(如图 3所示)，防止远心点偏移术前规划位置发生手术事故。

2 远心点调整机构尺寸优化

远心点调整机构对机械臂的灵活摆位有重要的影响，因此针对远心点调整机构的尺寸参数进行优化。

2.1 远心点调整机构运动学分析

根据通用手术机械臂的结构分析，建立被动关节部分的运动学坐标系，如图 4所示。机构初始的连杆参数如表 1所示。

根据图 4中的运动学坐标系及表 1中的D-H参数，可以得到各个关节之间的齐次转换矩阵：

$ {}^0{\mathit{\boldsymbol{T}}_1} = {\left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&{{a_1}}\\ 0&0&0&1 \end{array}} \right]^1} $

$ {\mathit{\boldsymbol{T}}_2} = \left[ {\begin{array}{*{20}{c}} {\cos {\theta _1}}&{ - \sin {\theta _1}}&0&{{a_2}\cos {\theta _1}}\\ {\sin {\theta _1}}&{\cos {\theta _1}}&0&{{a_2}\sin {\theta _1}}\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right] $

$ {}^2{\mathit{\boldsymbol{T}}_3} = {\left[ {\begin{array}{*{20}{c}} {\cos {\theta _2}}&{ - \sin {\theta _2}}&0&{{a_3}\cos {\theta _2}}\\ {\sin {\theta _2}}&{\cos {\theta _2}}&0&{{a_3}\sin {\theta _2}}\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right]^3} $

$ {\mathit{\boldsymbol{T}}_n} = \left[ {\begin{array}{*{20}{c}} {\sin {p_1}\cos {\theta _3}}&0&{\sin {\theta _3} - \cos {p_1}\cos {\theta _3}}&0\\ {\sin {p_1}\sin {\theta _3}}&0&{ - \cos {\theta _3} - \cos {p_1}\sin {\theta _3}}&0\\ {\cos {p_1}}&0&{\sin {p_1}}&{{a_4}}\\ 0&0&0&1 \end{array}} \right] $

式中p₁=π/4。

联立关节之间的转换矩阵，即可得到远心点坐标系x_ny_nz_n在基坐标系x₀y₀z₀下的位姿矩阵:

$ {}^0{\mathit{\boldsymbol{T}}_n} = {}^0{\mathit{\boldsymbol{T}}_n}{}^1{\mathit{\boldsymbol{T}}_n}{}^2{\mathit{\boldsymbol{T}}_n}{}^3{\mathit{\boldsymbol{T}}_n} $

(1)

各关节的运动范围为

$ \left\{ {\begin{array}{*{20}{c}} {{a_1} \in \left[ {800,1400} \right]{\rm{mm}}}\\ {{\theta _1} = \left[ { - 5{\rm{ \mathsf{ π} }}/6,5{\rm{ \mathsf{ π} }}/6} \right]{\rm{rad}}}\\ {{\theta _2} = \left[ { - 5{\rm{ \mathsf{ π} }}/6,5{\rm{ \mathsf{ π} }}/6} \right]{\rm{rad}}}\\ {{\theta _3} = \left[ { - 2{\rm{ \mathsf{ π} }}/3,2{\rm{ \mathsf{ π} }}/3} \right]{\rm{rad}}} \end{array}} \right. $

由于第一被动关节对机械臂的运动性能不存在影响，因此优化过程不考虑该关节。为验证运动学模型的正确性，在Matlab中搭建Simulink模型，如图 5所示。以被动关节Simmechanics模型的输出为实际值，运动学模型的输出为理论值。设置仿真时间为10 s，运行仿真可以得到关节的输入角度和远心点位置误差，如图 6所示。仿真结果显示建立的被动关节的运动学模型是正确的。

2.2 建立综合灵巧度评价指标

根据远心点调整机构的正运动学模型，可以得到远心点在基座标系下的坐标：

$ \left[ {\begin{array}{*{20}{c}} {{p_x}}\\ {{p_y}}\\ {{p_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{a_2}\cos {\theta _1} + {a_3}\cos \left( {{\theta _1} + {\theta _2}} \right) + {a_4}\sin {p_1}\cos \left( {{\theta _1} + {\theta _2} + {\theta _3}} \right)}\\ {{a_2}\sin {\theta _1} + {a_3}\sin \left( {{\theta _1} + {\theta _2}} \right) + {a_4}\sin {p_1}\sin \left( {{\theta _1} + {\theta _2} + {\theta _3}} \right)}\\ {{a_1} - {a_4}\sin {p_1}} \end{array}} \right] $

式中: a₁、a₂、a₃、a₄为远心点机构几何尺寸，θ₁、θ₂、θ₃为关节转角。

由于第一被动关节对机构的灵巧度没有影响，因此不考虑第一被动关节。雅可比矩阵可以通过对θ₁、θ₂、θ₃求导得到：

$ \mathit{\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {p_x}}}{{\partial {\theta _1}}}}&{\frac{{\partial {p_x}}}{{\partial {\theta _2}}}}&{\frac{{\partial {p_x}}}{{\partial {\theta _3}}}}\\ {\frac{{\partial {p_y}}}{{\partial {\theta _1}}}}&{\frac{{\partial {p_y}}}{{\partial {\theta _2}}}}&{\frac{{\partial {p_y}}}{{\partial {\theta _3}}}}\\ {\frac{{\partial {p_z}}}{{\partial {\theta _1}}}}&{\frac{{\partial {p_z}}}{{\partial {\theta _2}}}}&{\frac{{\partial {p_z}}}{{\partial {\theta _3}}}} \end{array}} \right] $

(3)

机械臂的运动性能评价指标有很多，这些指标也被广泛应用于机构尺寸的优化中。其中应用最广泛的是机构的灵巧度，机构灵巧度是指末端沿任意方向到达指定位姿的运动能力。雅可比矩阵的条件数是一个重要的灵巧度指标，条件数越接近于1，说明在该位姿下的运动性能越优^[14]，定义如下

$ {K_J} = \left\| \mathit{\boldsymbol{J}} \right\| \cdot \left\| {{\mathit{\boldsymbol{J}}^{ - 1}}} \right\| $

(4)

式中J为雅可比矩阵：

$ \left\| \mathit{\boldsymbol{J}} \right\| = \sqrt {{\lambda _{\max }}\left( {{\mathit{\boldsymbol{J}}^{\rm{T}}}\mathit{\boldsymbol{J}}} \right)} = {\sigma _{\max }}\left( \mathit{\boldsymbol{J}} \right) $

(5)

式中: λ_max为雅可比矩阵最大特征值, σ_max为雅可比矩阵最大奇异值。

由于条件数只与机械臂关节的转角和结构尺寸有关，反映了机械臂在指定位姿下的运动能力，但不能反映在全域空间内的运动能力，对此提出了全域空间条件数均值(GCI)这一指标，即对条件数在整个运动空间内取平均值，这一指标可以反映在整个运动空间内运动学灵巧度和位置误差的总体性能:

$ {\rm{GCI}} = \frac{{\int\limits_W {{K_J}{\rm{d}}w} }}{{\int\limits_W {{\rm{d}}w} }} $

(6)

式中: K_J为指定位姿下的条件数，W为全域工作空间。

考虑到GCI指标并没有反映出随着各个位置变化时的空间灵巧度波动情况，对此提出了运动学雅可比矩阵条件数波动性能指标(GCIF)，即在全域工作空间内对雅可比矩阵条件数取方差：

$ {\rm{GCIF}} = \sqrt {\frac{{\int\limits_W {\left( {{K_J} - {\rm{GCI}}} \right){\rm{d}}w} }}{{\int\limits_W {{\rm{d}}w} }}} $

(7)

考虑到以上两个指标对机械臂的运动灵巧度同样重要，因此对两个指标进行加权处理，作为综合灵巧度评价指标CDI：

$ {\rm{CDI}} = {k_1}{\rm{GCI}} + {k_2}{\rm{GCIF}} $

(8)

式中k₁、k₂为对应项加权系数。

2.3 设计优化变量及分析

在优化前利用初始参数对机械臂进行分析，根据实际手术中的摆位角度研究，关节的转角范围为35° ≤θ₁≤85°, -90°≤θ₂≤-55°，-160°≤θ₃≤-135°远心点的工作空间如图 7所示，条件数的分布如图 8所示。

由于运动学性能指标主要由关节角度变化范围和关节杆长所决定的，而常用角度取值范围主要依据手术过程中医生的观察角度及方便术前摆位要求等其他条件所确定，变化范围不大，故此项约束不做研究。三个关节杆长与角度变化范围决定了运动学性能、工作空间及关节内部零件的安装尺寸等，故以工作空间为约束，以三个杆长(a₂，a₃，a₄)为设计变量。考虑到关节2和关节3的结构相似，因此设置l₀=a₂=a₃，l₁=a₄sin p₁，限定杆长范围为140≤l₀≤260 mm，310≤l₁≤400 mm。

为了更全面的研究手术机械臂的运动目标性能，采用逐步求解的方法分析每个被动关节杆长对目标性能的灵敏度。通过编写M文件得到全域空间内条件数均值与条件数波动值在变量约束范围内的分布，如图 9、10所示。

从以上仿真结果可以得到：在设计变量约束范围内，随着杆长的变化，全域空间内的条件数均值变化范围为(4, 9)，而条件数波动值变化范围为(6, 10)，两者变化范围在同一数量级。考虑到这两项指标在优化指标中地位同等重要，将其进行线性加权处理，权重因子分别取0.5，即k₁=k₂=0.5。

2.4 优化算法和优化结果

通过编写M文件，分析设计变量对综合目标性能的灵敏度，如图 11所示。从仿真结果中可以看出综合灵巧度指标随着杆3的长度增加呈现先递减后递增的趋势；而随着杆1长度的逐渐增大，目标性能开始阶段一直递减，后半段变成一直递增的趋势。综合灵巧度指标随杆长的变化趋向如图 12所示。

减小设计变量的变化梯度，得到设计变量对综合灵巧度指标的影响关系如图 13所示，从中选择对应最优综合灵巧度指标的参数为最终的优化杆长，得到l₀=185 mm，l₁=370 mm。利用优化后的杆长参数分析远心点的工作空间和条件数的分布，相比于优化前的结果，远心点的工作空间基本不变，全域空间内的条件数分布更加平缓，综合灵巧度指标CDI提高了14.8%，全域条件数均值GCI降低了15.5%，全域条件数波动值GCIF减少了14.2%，如图 14、15所示。

3 结构分析和优化

首先对机械臂的结构进行分析，对整个手术机械臂进行静力学分析、振动模态分析和关键零部件的疲劳分析，从而得到主要零件的强度分布及可优化空间，机械臂静力学应力及形变云图如图 16所示。

从仿真结果中可以发现，整个机构应力值较低且分布很不均匀，存在较大的材料冗余，需要对其优化。

对整个机械臂来说，过渡关节不单要承受整个主动关节的重力以及力矩，还要承受手术过程中关节变化所带来的交变载荷，所以将过渡关节作为关键零部件进行分析。

关键零件的分析结果如图 17所示，由分析结果可知，该过渡关节即使在加大载荷的情况下，疲劳寿命数量级也很大而且安全因子主要集中在转轴的过渡变形部位，基本满足手术要求。

根据最优受力结构对小臂进行参数化建模，以最大应力、变形和振动频率为约束对其进行轻量化设计，分析及优化结果如图 18、19所示。

比较优化前后的指标可以得到：优化后各项指标有明显的提升，优化后被动关节质量减少了21.9%，一阶固有频率减少26.9%，最大应力减小了22%，如表 2所示。

4 结论

1) 借助于运动学雅可比矩阵的奇异值，构造了基于全域空间的条件数均值和条件数波动值加权的综合灵巧度评价指标，并利用该指标对机械臂杆长参数进行优化，使优化后的机械臂灵活性指标达到最优。

2) 基于灵巧度指标得到的最优参数，以质量最小为目标，以最大变形、最大应力、一阶固有频率为约束，对被动关节进行结构优化设计，得到最优的受力结构及尺寸。

3) 优化结果表明，优化后的手术机械臂具有良好的运动性能，且实现手术机械臂轻量化的设计要求，从而提高了机械臂的整体的运动性能。优化后的机械臂很好地满足了微创手术对机械臂机构灵活性的要求。

基于灵巧度指标对主动关节的参数进行优化和加权系数的细分将是下一步的工作。