近年来，微创手术机器人由于创伤小、恢复快、并发症少等优点得到广泛应用，在全世界范围内，许多公司及科研机构开展了手术机器人的研制工作[1]．医生操作主手实时控制从机械臂运动的过程中，机器人逆运动学的求解效率和准确性直接影响手术操作的效果．通常采用Denavit-Hartenberg(D-H)参数法对机械臂进行建模[2]，但是该方法几何意义不明显．微创手术机器人的手术器械的关节不满足Pieper准则，通常采用数值迭代和析配消元的方法求解，这影响了主从控制的求解效率和精度．Liu等[3]使基于迭代雅克比矩阵求解磁驱动可操纵血管导管系统的运动学，采用阻尼最小二乘法增强逆解的数值稳定性．夏进军等[4]提出一种改进的加速隐式近似线性规划方法求解逆运动学，该方法避免动力学奇异，并且有效减少迭代次数，但是存在线性化误差．于凌涛等[5]在归一化基系下逆解位姿矩阵，给出了主动关节及手术器械的唯一解，逆解结果复杂，计算量大．唐奥林[6]采用两步近似求解方法和误差补偿求解从机械臂的逆解，通过简单计算即可得到逆运动学的显示解析表达式，但是所得解为逆运动学的近似解．Brockett[7]将旋量指数积法引入到机器人建模中，在机器人运动学、动力学等领域取得了较好的效果[8]．采用旋量理论对机器人进行运动学建模，将机器人的逆运动学求解问题转换为若干Paden-Kahan子问题[9]求解．求逆解具有明确的几何意义，将机器人运动学的分析简化，求解过程简单．所得逆解为解析解，没有累积误差，因而被广泛应用在串联机器人[10]、冗余机器人、并联机器人[11-12]等复杂系统的运动学分析中．赵杰等[13]提出了RTS运动链并联机器人的子问题逆解方法．孔民秀等[14]提出了RRSR运动支链的子问题逆解，并应用到相应的运动支链中．钱东海等[15]基于消元理论和Paden-Kahan子问题结合的方法解决了后三个关节交于一点而前三个关节轴线没有交点的六自由度机器人的逆解问题，具有一定的通用性．陈庆诚等[16]针对“钱江一号”的结构提出了一种新的逆解子问题并给出了逆解方法，完成了六自由度机器人的逆解，并扩展了逆解子问题的适用范围．Leoro等[17]针对相邻两旋转关节平行，并与另一旋转关节垂直的结构提出了针对拟人手臂的新子问题的求解方法，推导了适用于6R机器人的运动学逆解．Xu等[18]针对冗余机械臂提出了三种新型子问题，采用几何约束和代数约束进行逆运动学求解，可应用于多关节冗余机械臂的逆解问题．ZHAO等[19]针对相邻三关节平行的结构提出了新子问题，采用消元理论和几何约束求解出8组逆解，逆解下机械臂末端的位置误差和姿态误差均满足机器人的精度要求．WANG等[20]提出了针对任意三个旋转关节的逆解方法，针对平行与相交等特殊位姿构型均明确说明，可应用于任意的三个回转关节的逆解．本研究提出了三个关节轴线空间垂直的PRR类型的Paden-Kahan子问题，并给出了该子问题逆运动学的求解方法．该子问题描述为一个点绕空间垂直的两个轴线旋转，再沿直线方向移动至指定位置．以主动关节轴线汇交于一点、手术器械关节轴线空间垂直的华鹊III微创手术机器人为实验平台，通过任意给定运动量的实例计算，验证该子问题的运动学逆解方法的有效性．1 机器人运动学模型约定某一符号表示某一个点或者是某一个值，黑斜体的符号表示三维位置矢量、齐次位置矢量或矩阵．对于转动关节，其轴线的方向表示为ω=[ωxωyωz]T，关节位置表示为r=[rxryrz]T，转动关节的单位运动旋量表示为ξ=ωr×ω．(1)对于移动关节，其轴线的方向表示为v=[vxvyvz]T，移动关节的单位运动旋量表示为ξ=03×1v．(2)根据欧拉定理，指数映射eω^θ与旋转矩阵R∈SO(3)是等价的．刚体的螺旋运动采用指数坐标形式表示为eξ^θ=I3×3vθ03×11    (ω=0);eω^θ(I3×3-eω^θ)[ω×(r×ω)]+θωωTv03×11    (ω≠0), (3)式中：θ为广义关节角度；ω^为ω的反对称阵，ω^=0-ωzωyωz0-ωx-ωyωx0；ξ^为旋量坐标ξ的瞬时运动关节旋量．刚体上的工具坐标系相对于惯性坐标系的初始位姿变换为gst(0)，刚体变换运动后的工具坐标系相对惯性坐标系的位姿为gst(θ)=eξ^θgst(0)．(4)2 PRR子问题运动学逆解新子问题的几何关系如图1所示，图中：v1为第一移动关节的方向向量；ω2和ω3为第二和第三回转关节的方向向量；v1''为过点q与向量v1平行的向量．点p2和p3分别位于旋量ξ2和ξ3的轴线上，其位置是任意选取的．新的子问题中三个运动旋量ξ1,ξ2,ξ3空间垂直不交叉，点p绕旋量ξ3的运动轨迹以圆C3表示，旋转角度θ3至点c1；点c1绕旋量ξ2的运动轨迹以圆C2表示，旋转角度θ2至点c2；点c2沿旋量ξ1运动至点q，满足10.13245/j.hust.230209.F001图1新子问题的几何关系示意图q=eξ^1θ1eξ^2θ2eξ^3θ3p．(5)2.1　构造投影点由于第一个关节是移动关节，无法采用消元法消除第一个变量，而且p1位置的选择与计算无关．三个运动旋量在空间中两两垂直，因此无须建立公垂线．令u31=p-r3，易知：r3=p3+ω3ω3T(p-p3);u31=p-r3;u31ω3=0. (6)令u11=q-r2，易知：r2=p2+ω2ω2T(q-p2);u11=q-r2;u11ω2=0. (7)令u12=n-r2=u11-v1v1Tu11．2.2　构造几何关系并求逆解点c1在运动轨迹C3上，令u32=c1-r3，易得u32=u31．(8)由图1(b)可知，fc1=ω2ω2T(r2-r3)，在△fc1r3中，由勾股定理可得fr3=u312-fc12．(9)由此可以计算关节角u32=fc1±fr3v1,θ3=arccosu31u32u31u32. (10)由上式可得θ3∈[0,π]，但不确定正负．因此取u31×u32=λ3ω3，若λ30，则为正值；若λ30，则为负值．求解θ3后，可以得到点c1的位置为c1=eξ^θ3p，点c1和点c2都位于绕旋量ξ2的运动轨迹C2上，令u21=c1-r2，u22=c2-r2，可知u21=u22．(11)由图1(c)可知，在△nc2r2中，满足nc2=u212-u122．(12)由此可以计算关节角u22=u12±nc2v1,θ2=arccosu21u22u21u22. (13)同理，由上式可解得θ2∈[0,π]，但不确定正负．因此，取u21×u22=λ2ω2，若λ20，则为正值；若λ20，则为负值．求解θ2后，可以得到点c2的位置为c2=eξ^θ2c1，可得θ1=(q-c2)v1．(14)2.3　逆解个数的判断与取舍该子问题逆解的个数取决于各运动轨迹交点的个数，可能存在4个、2个、1个解或无解的情况．若部分关节有角度范围限制，则在求解过程中，对当前关节角度求解之后判断角度是否满足约束．若满足，则保留；若不满足，则舍去．最后，应当对所求逆解进行正运动学求解，判断其正解的位姿是否与期望位姿相同．3 逆解算法的实例验证微创手术串联机器人具有10个自由度，结构简图如图2所示，图中：ωi为第i个旋转关节的方向向量；ri为第i个旋转关节的参考位置；vi为第i个移动关节的方向向量；{S}和{T}分别为基座坐标系和工具坐标系；θi为第i个关节的广义运动量；d1，a2，a3，α4，d5和a9均为机械臂结构参数常量．前4个自由度为被动关节，用于术前调整远心点的位置和主动关节介入人体的初始姿态，在术中被动关节保持静止，因此被动关节的位移值和角度值是已知的．主动关节包括3个自由度，分别为绕远心点的俯仰、翻滚和沿直线方向位移．手术器械包括4个自由度，其中小爪的开合由主操作手单独控制，此处不予考虑，因此该手术器械具有3个自由度，分别为绕自身器械杆的旋转、腕部俯仰、小爪翻滚．关节7和关节8同轴，且关节7为移动关节，关节8为回转关节，基于互换关节顺序不影响末端位姿的前提，可以互换两个关节的顺序[21]．由于互换关节7和关节8的顺序后不影响末端位姿，为了便于求解关节的逆解，此处将关节7和关节8顺序进行调换．10.13245/j.hust.230209.F002图2微创手术机器人结构简图3.1　微创手术机器人正运动学建立微创手术机器人各运动关节的旋量坐标：ω1=ω8=000;ω2=ω3=ω4=001;ω5=cos α40-sin α4;ω6=ω9=0-10;ω7=00-1;ω10=100. (15)确定各关节的轴线点：r2=00d1;r3=a20d1;r3=a2+a30d1;r5=r6=r7=r9=a2+a3+d5cos α40d1-d5sin α4;r10=a2+a3+d5cos α40d1-d5sin α4-a9. (16)手术器械末端的运动学正解为gst(θ)=∏i=110eξ^iθigst(0)，(17)式中gst(0)为初始状态下工具坐标系的位姿，是已知矩阵．根据给定的期望位姿求解各关节的角度值，因此gst(θ)是已知的变换矩阵．3.2　微创手术机器人逆解过程微创手术机器人的被动关节在术中保持静止，其关节角度值可由编码器读取，θ1~θ4是已知数值，微创手术串联机器人的逆解实质上是RRR-PRR六自由度机器人逆解．式(17)可变换为∏i=14eξ^iθi-1gst(θ)gst-1(0)=∏i=510eξ^iθi．(18)令g1=∏i=14eξ^iθi-1gst(θ)gst-1(0)，式(18)变换为g1∏i=810eξ^iθi-1=∏i=57eξ^iθi．(19)关节5~7均为回转关节，且轴线汇交于一点，在式(19)两边均右乘汇交点的齐次坐标矢量，并根据位置不变原则，可得∏i=810eξ^iθig1-1p5=p5，(20)式中p5=[r5T1]T．令q1=g1-1p5，若关节8、关节9、关节10的构型配置符合子问题型式，则该式求解与所提子问题求解相同．求解θ8，θ9和θ10后，须求解三轴汇交于一点的θ5，θ6和θ7．由式(19)两边同时右乘关节7轴线上不同于点p7的任意点p7'的齐次坐标p7，可得g1∏i=810eξ^iθi-1p7=∏i=57eξ^iθip7．(21)令q2=g1∏i=810eξ^iθi-1p7，根据位置不变原则，式(21)可变换为q2=eξ^5θ5eξ^6θ6p7．(22)式(22)可根据文献[22]中所提方法求解，此处不再赘述．求解θ5和θ6后，代入式(19)中，可得∏i=56eξ^iθi-1g1∏i=810eξ^iθi-1=eξ^7θ7．(23)将式(23)两边右乘关节6轴线上不同于点p7的任意点p6'的齐次坐标p6，可得∏i=56eξ^iθi-1g1∏i=810eξ^iθi-1p6=eξ^7θ7p6．(24)令q3=∏i=56eξ^iθi-1g1∏i=810eξ^iθi-1p6，式(24)可根据文献[22]中所提方法求解，此处不再赘述．至此，主动关节及手术器械的所有关节逆解均可求得．3.3　微创手术机器人逆解验证微创手术机器人在手术中由RRR-PRR六个自由度运动决定末端位姿，前三个自由度汇交于空间不动点，后三个自由度空间垂直．通过实例计算验证本研究提出的子问题逆解求解方法的有效性和可行性．对微创手术机器人各关节建立坐标系，如图2所示．图2中各参数值为d1=1 000 mm，a2=500 mm，a3=500 mm，α4=π/4，d5=800 mm，a9=10 mm．各关节运动范围为θ5∈[-π/2，π/2]，θ6∈[-π/3，π/3]，θ7∈[-π，π]，θ8∈[100，350] mm，θ9∈(-π/2，π/2)，θ10∈(-π/2，π/2)．由于关节1~4为直接读取编码器数值，在此略去其运动范围．微创手术机器人初始状态下各关节旋量值为ξ1=[0，0，0，0，0，1]T，ξ2=[0，0，1，0，0，0]T，ξ3=[0，0，1，0，-500，0]T，ξ4=[0，0，1，0，-1 000，1]T，ξ5=[0.707 11，0，-0.707 11，0，1 414.213 56，0]T，ξ6=[0，-1，0，434.31，0，-1 565.685 42]T，ξ7=[0，0，-1，0，1 565.685 42，0]T，ξ8=[0，0，0，0，0，-1]T，ξ9=[0，-1，0，434.314 58，0，-1 565.685 42]T，ξ10=[1，0，0，0，424.314 58，0]T． 惯性坐标系设置在基座上，工具坐标系固定在小爪轴线上，姿态与惯性坐标系相同．在机器人初始状态(θ=[0,0,0,0,0,0,0,0,0,0]T)下，工具坐标系的初始位姿为gst(0)=1001 565.685 420100.000 00001424.314 580001.000 00．给定各关节一组运动函数，从中取若干点，进行正运动学计算，以计算得到的位姿为期望位姿，应用提出的机器人逆解算法计算该位姿下机器人关节的角度．首先计算出手术器械PRR关节的逆解，从中筛选出满足约束要求的逆解，进而对主动关节求解逆解．给定术中静止不动关节的运动量为θ1=200，θ2=π9，θ3=-π18，θ4=-π18．给定各主动关节函数为：θ5=π2sinπ25k;θ6=π3sinπ25k;θ7=πsinπ25k;θ8=9(k-1)2+105;θ9=4π9cosπ25k;θ10=-4π9cosπ25k,式中k=1,2,⋯,51．给定运动量与逆解的对比如图3所示，图中：X为采样点；Y1为关节角度值；Y2为关节移动值．各曲线代表各关节的运动函数，θ8曲线对应Y2，其余曲线对应Y1；曲线上标记“о”代表给定的运动量，曲线上标记“+”代表逆解所得运动量，可知给定运动量与逆解运动量的值是一致的．图4为给定运动量与逆解的末端点之间的位置误差(Er)．由图4可知：给定的运动量与逆解的末端位置之间的最大误差为4.875 4×10-14 m，平均误差为4.839×10-15 m，满足微创手术机器人高精度的控制要求．通过实例计算，验证了所提PRR子问题，即一个点绕空间垂直的两个轴线旋转，再沿直线方向移动至指定位置子问题逆解方法的有效性．10.13245/j.hust.230209.F003图3给定运动量与逆解的对比10.13245/j.hust.230209.F004图4给定运动量与逆解的末端点之间的位置误差4 结语针对一种前三个旋转关节相交于一点、后三个关节空间垂直RRR-PRR的六自由度串联机器人的逆解问题，给出了一类PRR三轴空间垂直子问题的求解方法，末端点和参考点向各轴线的垂直平面进行投影，建立几何约束进行求解，表明逆解方法具有明确的几何意义．该方法解决了不满足Pieper准则的六自由度机器人逆解问题．由于机器人特殊的构型，参考点选择为轴线交点，计算过程更加简便．提出的逆解算法在逆解范围约束下，能够迅速筛选出满足要求、且与期望位姿一致的机器人逆解．逆解算法仿真实例中，给定运动量与逆解的末端位置之间的最大位置误差为4.875 4×10-14 m，满足机器人的高精度要求．