近年来，由于核磁共振成像(Magnetic Resonance Imaging, MRI)独特的软组织分辨能力，MRI介入下经皮穿刺手术变得越来越流行。与计算机化X线体层照相术 (Computerized Tomography, CT) 相比，MRI有更好的软组织对比效果，且无电离辐射^[1]，这使得它非常适合神经外科手术。

核磁共振引导的机器人系统已经得到了广泛的研究，并已经被证实可以用于前列腺活检^[2-3] 、乳腺活检^[4-5]、胸腹部干预^[6-7]和大脑穿刺^[8-9]等。然而，在MRI约束空间中的穿刺行为和强磁场环境仍然对进行穿刺机器人辅助手术提出了准确性和安全性的挑战。

为将医生工作与核磁共振空间分离，Patel 协同其小组成员使用超声波电机、机械手及其电子器械搭建了一个简单的机器人辅助系统，将消融器放进脑中进行MRI引导下的消融实验，该机器人系统的平均平移和旋转误差精度为(1.39 ± 0.64) mm和(1.27 ± 0.56)°^[8]。

同时，为了应对核磁共振室的强磁场问题，如何拟定传动系统也是实现MRI设备兼容的一个重要问题。Johns Hopkins University实验室、Automation and Interventional Medicine实验室都曾经使用超声波电机来设计MRI兼容手术机器人，但稳定性和精度并不太理想^[10-12]，而且基于超声波电机的机器人，其信噪比(Signal Noise Ratio) 相对较高，高于10%。因此，Franco等^[13]开发了一种没有任何金属材料的气动驱动机器人，其信噪比始终低于5%。后续还有很多学者开发了气动驱动器或液压驱动器。Dong等^[14]设计了一款三缸液压马达，并使用该马达设计了液压驱动的手术机器人，同时结合了一种新型的磁共振无线跟踪技术，使得它能够在适当的磁共振跟踪序列下提供实时连续的机器人仪器三维定位^[9]。Farimani等^[15]提出了一种气动步进执行器，该执行器可以代替超声波电机作为磁共振兼容机器人的执行器。但这些驱动器的尺寸太大，无法直接应用于内径核磁共振引导活检。Frishman等^[16]开发了一种带有力反馈的液压驱动MRI内孔引导活检装置，可以通过观察反馈力来区分针是否穿透膜。然而，该设备在手术过程中不能移动，它只有一个自由度(Degree of Freedom, DOF)。

因此，基于流体驱动的核磁共振兼容设备有重要的研究价值，本文重点研究了一种基于连通器原理的液压驱动穿刺手术机器人的动力学模型，以及其 ${H}_{\infty } $ 控制方法，并通过仿真和实物平台验证了其精度。

1 实验平台简介 1.1 总体方案概述

MRI设备兼容穿刺手术机器人整体方案设计如图1所示。其具有的二自由度远端运动中心结构，如图2所示，可以在仅造成微小创口的情况下在人体内进行较大范围的运动。同时，本文采用了非金属材质进行搭建并采用液压的驱动方式，使该机器人可以兼容MRI设备。其主体框架主要由玻璃纤维、塑料轴承以及基于光敏树脂的3D打印件构成，并由陶瓷制成的螺丝、螺母作为紧固。此外，采用了光纤以及亚克力材质的直线光栅尺制成编码器，作为该机器人的位置反馈系统。

此六自由度机器人可以被分为4个独立的功能模块：

(1) 液压驱动模块：由位于MRI设备室外的滚珠丝杆驱动模组作为执行器，利用液压的方式进行机器人本体的控制；

(2) 穿刺点定位模块：用于将穿刺针定位在目标位置，其可以沿着坐标轴的x轴、y轴、z轴运动；

(3) 穿刺针位姿调整模块：用于调整穿刺针的位姿，以通过头颅的小孔（通常由外科医生在将患者置于扫描仪之前钻出）对患者实施穿刺手术；

(4) 穿刺针介入模块：用于调整探针刺入患者脑内的深度。

2 机器人的数学模型 2.1 机器人系统的正逆运动学模型

图3 描述了机器人每个轴的名称、位置、参考坐标系，以及相应的参数。将机器人的位姿矢量定义为：L = $ \left({x, y,{\textit{z}} , {\textit{θ}}, {\textit{φ}} }\right) $ ，其中x为机器人在x轴相对于初始位置的位移；y为机器人在y轴相对于初始位置的位移；z为机器人在z轴相对于初始位置的位移； $ {\textit{θ}} $ 为pitch轴旋转，是穿刺针与图中xz平面的夹角，以穿刺针垂直于xy平面时为零参考点； $ {\textit{φ}} $ 为roll轴旋转，是穿刺针与图中yz平面的夹角，以穿刺针垂直于xy平面为零参考点。定义机器人关节变量广义坐标为V = $ \left(a,b,c,d,e,f\right) $ ，其中， $ a,b,c,d,e $ 依次为驱动端 $ {x, y, {\textit{z}}, {\textit{θ}}, {\textit{φ}} } $ 轴的活塞的移动量，f为控制穿刺针的刺入深度的活塞的移动量。机器人相对于参考坐标系初始坐标为 $ {P}_{0}\text{=}\left({x}_{0},{y}_{0},{\textit{z}}_{0}\right) $ ，其中 $ {x}_{0} $ =123 mm, $ {y}_{0} $ =446.5 mm, $ {\textit{z}}_{0} $ =164 mm。机器人pitch轴和roll轴的力臂分别表示为 $ {l}_{\text{1}} $ 和 $ {l}_{\text{2}} $ ，其中 $ {l}_{\text{1}}$ =40 mm ， $ {l}_{\text{2}}$ =38 mm 。

机器人运动学问题包括机器人正运动学和机器人逆运动学，正运动学问题为给定机器人各关节变量，求解机器人末端的位置姿态；逆运动学问题为已知机器人末端的位置姿态，求解该位姿对应的关节变量。

根据上文的描述以及机器人结构，容易推得机器人的正运动学表达式为

$ \left\{\begin{split}&\left[ \begin{array}{c}{x}\\ {y}\\ {z}\end{array} \right]=\left[ \begin{array}{cccc}\text{1}& \text{0}& \text{0}& \dfrac{e}{{l}_{2}}\\ \text{0}& \text{1}& \text{0}& \dfrac{d}{{l}_{1}}\\ \text{0}& \text{0}& \text{1}& -\text{cos}\left(\text{arc}\text{sin}\left(\dfrac{{d}}{{{l}}_{{1}}}\right) \right) \text{cos}\left(\text{arc}\text{sin}\left(\dfrac{{e}}{{{l}}_{{2}}}\right) \right) \end{array}\right]\times\\ &\qquad\;\;\; \left[\begin{array}{c}{a}\\{b}\\ {c}\\ {f}\end{array}\right]\text+\left[\begin{array}{c}{x}_{\text{0}}\\ {y}_{\text{0}}\\ {z}_{\text{0}}\end{array}\right]\\ &{\textit{θ}}\text=\text{a}\text{r}\text{c}\text{sin}\left(\dfrac{{d}}{{{l}}_{\text{1}}}\right) \\ &{\textit{φ}}\text=\text{a}\text{r}\text{c}\text{sin}\left(\dfrac{{e}}{{l}_{\text{2}}}\right) \end{split}\right. $

(1)

将机器人的运动分为2个阶段：第一个阶段为机器人定位阶段，该阶段机器人的f 为0，不控制穿刺针的刺入深度，仅将远端控制中心(Remote Center of Motion, RCM)点定位到人体颅骨的某一位置，其逆运动学表达式为

$ \left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right] \left[\begin{array}{c}x\\ y\\ {{z}}\end{array}\right]-\left[\begin{array}{c}{x}_{0}\\ {y}_{0}\\ {{z}}_{0}\end{array}\right] $

(2)

第二个阶段为穿刺针的刺入阶段，此时a, b, c均为0，表示机器人的RCM点不移动，仅通过穿刺针刺入深度和刺入角度的调节来控制针头末端的位置，其逆运动学表达式为

$\left\{\begin{split}d=&\mathrm{sin}\;\theta \times (x-{x}_{0}) +\mathrm{sin}\;\varphi \times (y-{y}_{0}) -\text{c}\text{o}\text{s}\;\theta \text{c}\text{o}\text{s}\;\varphi \times \\& ({\textit{z}}-{{\textit{z}}}_{0}) \\ e=&{l}_{1} \mathrm{sin}\;\theta \\ f=&{l}_{2} \mathrm{sin}\;\varphi \end{split}\right. $

(3)

2.2 机器人系统的动力学模型 2.2.1 液压驱动系统的等效模型

首先规定，在这一由基于活塞−缸体的液压连通器构成的液压驱动系统中，位于核磁共振室外的滚珠丝杆驱动模组的丝杆滑块，称为机器人驱动端；位于机器人上的导向轴承盒，称为机器人的从动端，如图4所示。其次，系统的输入形式为驱动端速度。这是由于采用了滚珠丝杆系统作为驱动，选用了型号为CBX1204-100ZC的滚珠丝杆滑台导轨线性滑台，其行程为10 cm，螺距为4 mm。57HS09型号步进电机，扭矩为0.9 N·m，根据滚珠丝杆扭矩计算理论驱动力，如式(4)所示，其中 $ \eta $ 为正效率， $ \tau $ 为驱动扭矩， $ I $ 为丝杆的螺距。

$ {F=}\frac{\text{2}\text{π}{\textit{×η×τ}}}{I}\text=\frac{\text{2}\text{π}\text{×0.95×900}}{\text{4}}\text{=1\;343} \;{\rm{{N}}} $

(4)

理论驱动力为1343 N，实际机器人系统质量为2.5 kg，则加速度可达537.2 m/s²。在本模型中可以认为若驱动端滑块的速度在1 m/s的范围内，则其被视为一个可以突变的量。因此在本模型中，不以驱动端的力作为系统输入，而是以驱动端速度作为系统输入。

将该液压驱动系统等效为弹簧−阻尼模型，其中 $ {{x}}_{\text{robot}} $ 为机器人从动端移动的距离， ${x}_{\text{push}}$ 和 ${x}_{\text{pull}}$ 分别为驱动端滑块在压缩流体和抽吸流体的过程中移动的距离，理论上有 ${x}_{\text{push}}={x}_{\text{pull}}$ ， $ {f}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{t}} $ 为机器人从动端受到的干摩擦阻尼， $ {f}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{t}} $ 可表示为

$ {f}_{\text{robot}}\text=\left\{\begin{array}{c}{K}\left({x}_{\text{push}}-{x}_{\text{robot}}\right) {,}\; \;{\dot{x}}_{\text{robot}}\leqslant 0\\ {f}_{\text{slide}}\text{s}\text{g}\text{n}\left({\dot{x}}_{\text{robot}}\right) \text{,}\; \;{\dot{x}}_{\text{robot}} > 0\end{array}\right. $

(5)

式中：sgn为符号函数， $ {f}_{\mathrm{s}\mathrm{l}\mathrm{i}\mathrm{d}\mathrm{e}} $ 为机器人的滑动摩擦力， $ {f}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{c}} $ 为机器人的静摩擦力。图4 的液压驱动系统等效模型所示的机器人的动力学模型可以表示为

$ \left({m}_{\text{liquid}}\text+{m}_{\text{robot}}\right) {\ddot{x}}_{\text{robot}} = {N}{K}(x-{x}_{\text{robot}}) -{N}{C}_{\text{liquid}}\left(\dot{x} - {\dot{x}}_{\text{robot}}\right) - \text{}{f}_{\text{robot}} $

(6)

式中： ${m}_{\text{liquid}}$ 为机器人驱动系统中液体的质量， ${m}_{\text{robot}}$ 为机器人自身的质量，K为流体等效弹簧刚度， $ {C}_{\text{liquid}} $ 为液体的阻尼系数， $ {N} $ 为驱动机器人在某个方向上运动所需要的连通器的对数，以两个活塞−缸体系统为一对连通器。依此可以构造驱动端位移 $ {x} $ 与机器人实际位移 ${x}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{t}}$ 的关系。

对于机器人中部分机构，比如z轴在运动的过程中还要受重力的影响，因此，该模型并不是可以完全照搬在机器人的每个轴上，还需要针对具体的运动方向对式(4)进行微调整：

$ \left\{\begin{array}{l}\left({m}_{x}\right) {\ddot{x}}_{{{\rm{r}}x}}={2}{K}\left({x}_{x}-{x}_{{{\rm{r}}x}}\right) -{{2}{C}}_{{l}}\left({\dot{x}}_{x}-{\dot{x}}_{{{\rm{r}}x}}\right) -{{f}}_{{{\rm{r}}x}}\\ \left({m}_{{y}}\right) {\ddot{x}}_{{{\rm{r}}y}}={4}{K}\left({x}_{{y}}-{x}_{{{\rm{r}}y}}\right) -{{4}{C}}_{{l}}\left({\dot{x}}_{{y}}-{\dot{x}}_{{{\rm{r}}y}}\right) -{{f}}_{{{\rm{r}}y}}\\ \left({m}_{{{\textit{z}}}}\right) {\ddot{x}}_{{{\rm{r}}{\textit{z}}}}={K}\left({x}_{{{\textit{z}}}}-{x}_{{{\rm{r}}{\textit{z}}}}\right) -{{C}}_{{l}}\left({\dot{x}}_{{{\textit{z}}}}-{\dot{x}}_{{{\rm{r}}{\textit{z}}}}\right) -{{f}}_{{{\rm{r}}{\textit{z}}}}-{m}_{{{\textit{z}}}}g\\ \left({m}_{{p}}\right) {\ddot{x}}_{{{\rm{r}}p}}={K}\left({x}_{{p}}-{x}_{{{\rm{r}}p}}\right) -{{C}}_{{l}}\left({\dot{x}}_{{p}}-{\dot{x}}_{{{\rm{r}}p}}\right) -{{f}}_{{{\rm{r}}p}}\\ \left({m}_{{r}}\right) {\ddot{x}}_{{{\rm{r}}r}}={2}{K}\left({x}_{{{{r}}}}-{x}_{{{\rm{r}}r}}\right) -{{2}{C}}_{{l}}\left({\dot{x}}_{{{{r}}}}-{\dot{x}}_{{{\rm{r}}r}}\right) -{{f}}_{{{\rm{r}}r}}\\ \left({m}_{{d}}\right) {\ddot{x}}_{{{\rm{r}}d}}={K}\left({x}_{{d}}-{x}_{{{\rm{r}}d}}\right) -{{C}}_{{l}}\left({\dot{x}}_{{d}}-{\dot{x}}_{{{\rm{r}}d}}\right) -{{f}}_{{{\rm{r}}d}}\end{array}\right. $

(7)

式中： $ {x}_{{\rm{r}}x} $ ， $ {x}_{{\rm{r}}y} $ ， $ {x}_{{\rm{r}}z} $ ， $ {x}_{{\rm{r}}p} $ ， $ {x}_{{\rm{r}}r} $ ， $ {x}_{{\rm{r}}d} $ 分别为机器人在x，y，z，pitch，roll轴移动距离和刺入深度。而 $ {x}_{x} $ ， $ {x}_{y} $ ， $ {x}_{z} $ ， $ {x}_{p} $ ， $ {x}_{r} $ ， $ {x}_{d} $ 分别为各个轴对应的驱动活塞的移动距离。 ${m}_{x}$ ， ${m}_{y}$ ， ${m}_{z}$ ， ${m}_{p}$ ， ${m}_{r}$ ， ${m}_{d}$ 分别为机器人在沿x，y，z，pitch，roll，deep轴运动时的负载质量，亦即机器人各部分的质量，g为重力加速度。 $ {f}_{{\rm{r}}x} $ ， $ {f}_{{\rm{r}}y} $ ， $ {f}_{{\rm{r}}{\textit{z}}} $ ， $ {f}_{{\rm{r}}p} $ ， $ {f}_{{\rm{r}}r} $ ， $ {f}_{{\rm{r}}d} $ 分别表示机器人在沿x，y，z，pitch，roll，deep轴运动时，所受到的干摩擦阻尼。若K和C前面系数为2，则表示是对称作用器驱动，为1是非对称作用器驱动，为4则表示是两组对称作用器驱动。对称作用器如图4所示，而非对称作用器则为在图4的基础上，除去任意一边的液压连通器。

2.2.2 基于流体力学的运动学模型参数整定

本文机器人系统采用的管路为聚氨酯管，其形变可以忽略不计，因此在后续的建模中不考虑管路因液压而膨胀的情况。

首先进行流体的弹性模量，即流体刚度K的推导，流体的弹性是由流体的可压缩性决定的，流体的压缩公式为

$ {\textit{β}}=V\frac{{\Delta }{P}}{{\Delta }{V}}{=V}\frac{{\Delta }{F}}{{\Delta }{VS}} $

(8)

式中： ${\Delta }{P}$ 为液体压力变化量，V为液体原体积， ${\Delta }{V}$ 为液体体积变化量，S为提供压力的活塞的横截面积。为了找出流体等效弹簧刚度的表达式，式(8)可变化为

$ {\Delta }{F}={(} {\textit{β}}\frac{{{S}}^{{2}}}{{V}}{) }{\Delta }{h}={K} \times {\Delta }{x} $

(9)

式中： ${\Delta }{h} $ 为液面高度变化量， $ {\Delta }{x}$ 为等效弹簧压缩量。

一种描述气、液混合流体的等效流体刚度K的公式为

$ {K}{=  {\textit{β}}}\frac{{{S}}^{{2}}}{{V}}= {\textit{β}}_{{0}}\left[\frac{{1}+\dfrac{{{V}}_{\rm{{air}}}}{{{V}}_{\rm{{liquid}}}}}{{1}+{\left(\dfrac{{p}}{{{p}}_{{0}}}\right) }^{\dfrac{{1}}{\textit{γ}}} \times \dfrac{{{V}}_{\rm{{air}}}}{{{V}}_{\rm{{liquid}}}} \times \dfrac{ {\textit{β}}_{{0}}}{{{\textit{γ}}\text{∙}p}}}\right]\dfrac{{{S}}^{{2}}}{{V}} $

(10)

式中： ${\textit{β}}$ 为混入了空气后，液压驱动系统内所有流体的等效体积模量， ${\;\textit{β}}_{\text{0}}$ 为液体自身的体积模量，p为液压， $ {p}_{0} $ 为标准大气压，S为管道横截面积，V为液压驱动系统容器体积， ${V}_{\text{air}}$ 为容器内气体体积， ${V}_{\text{liquid}}$ 为容器内液体体积， ${\textit{γ}}$ 为比定压热容与空气比定容热容之比。

其次进行流体的阻力系数C的推导。流体在管道中流动时会受到沿程阻力的影响，某段管道上流体产生的总的能量损失 ${h}_{f}$ 用Darcy-Weishach公式表示为

$ {h}_{f}\text={\sum } {h}_{{ \lambda }}+\sum {h}_{\textit{ξ} }$

(11)

$ \displaystyle\sum {h}_{\lambda } $ 是由流体的粘滞力造成的损失，由式(12)描述。

$ {h}_{{ \lambda }}{= \lambda }\frac{L}{{d}}\frac{{v}^{\text{2}}}{\text{2}{g}}\text{} $

(12)

式中： $ \lambda $ 为沿程阻力系数(无量纲) ，L为管道的长度，d为管道的直径，v为管道有效截面上的平均流速。

$\displaystyle\sum {h}_{\xi }$ 是流体的局部碰撞导致的能量损失，由式(13)描述：

$ {h}_{\textit{ξ} }={\textit{ξ} }\frac{{{v}}^{{2}}}{{2}{g}} $

(13)

式中： $ \xi $ 为局部阻力系数(无量纲) ，多由实验测得，对于本文液压驱动系统所使用的管路如图5所示。根据尼古拉兹实验的结果，该种管路对应的局部阻力系数 $ \xi =0.5 $ 。

流体产生的总能量损失与对应摩擦力的关系可由黏性流体流动的伯努利方程得出。

$ {\textit{z}}_{{1}}+\frac{{{p}}_{{1}}}{{\textit{ρ}}g}+{\textit{α}}_{{1}}\frac{{{V}_{{1}}}^{{2}}}{{2}{g}}={\textit{z}}_{{1}}+\frac{{{p}}_{{2}}}{{\textit{ρ}}g}+{\textit{α}}_{{2}}\frac{{{{V}}_{{2}}}^{{2}}}{{2}{g}}+{{h}}_{{f}} $

(14)

式中：由于该液压驱动系统两个截面相同，即 $ {\textit{z}}_{1}={\textit{z}}_{2} $ ， $ {V}_{1}{=V}_{2} $ ， $ {\alpha }_{1}={\alpha }_{2} $ ，得

$ {{h}}_{{f}}=\frac{{{p}}_{{1}}-{{p}}_{{2}}}{{\textit{ρ}}g} $

(15)

$ {\Delta }{p}_{f}={p}_{1}-{p}_{2}={{\textit{ρ}}g}{h}_{f}=\frac{f}{S}{} $

(16)

因此液体流动沿程摩擦力f 可以表示为

$ {f=S\text{∙} \Delta }{{p}}_{{f}}=\mathrm{S}\mathrm{\rho }\mathrm{g}{{h}}_{{f}} $

(17)

由于步进电机最大驱动速度受限在0.2 m/s内，其雷诺数Re < 2000，处于层流区。由尼古拉兹图，沿程阻力系数 $ \lambda $ 在本文液压驱动系统中可以近似表示为

$ { \lambda =}\frac{{64}}{{Re}} $

(18)

由泊肃叶公式得

$ {{h}}_{{f}}=\frac{{32}{{\textit{μ}}l}}{{{\textit{ρ}}g}{{d}}^{{2}}}{v} $

(19)

式中： $ \mu $ 为流体黏度，则f 改写为

$ {f=S{{\textit{ρ}}g}}{{h}}_{{f}}{=S{{\textit{ρ}}g}}\frac{{32}{{\textit{μ}}l}}{{{\textit{ρ}}g}{{d}}^{{2}}}{v=}{8}{\text{π}}{{\textit{μ}}lv} $

(20)

取水的黏度系数为 $0.9\times {10}^{-3} \;\mathrm{P}\mathrm{a} \cdot \mathrm{s}$ ，管道长度为5 m，则等效阻尼系数

$ {C}_{\text{liquid}}\text=\text{8}{\text{π}}{{\textit{μ}}l}\text{=8}\text{×}\text{3.1416}\text{×}\text{0.9} \cdot {\text{10}}^{-3}\text{×}\text{5=0.113\;1} $

(21)

根据式(4)描述的机器人动力学方程式，机器人的状态空间方程可以写为

$ \left\{ \begin{array}{l}\left[\begin{array}{c}{\dot{{x}}}_{{1}}\\ {\dot{{x}}}_{{2}}\\ {\dot{{x}}}_{{3}}\end{array}\right]=\left[\begin{array}{ccc}{0}& {1}& {0}\\ -\dfrac{{2}{K}}{m}& \dfrac{{2}{C}}{m}& \dfrac{{2}{K}}{m}\\ {0}& {0}& {0}\end{array}\right]\left[\begin{array}{c}{{x}}_{{1}}\\ {{x}}_{{2}}\\ {{x}}_{{3}}\end{array}\right]+\\ \qquad\quad\;\;\; \left[\begin{array}{c}{0}\\ -\dfrac{{2}{C}}{m}\\ {1}\end{array}\right]{u}-\left[\begin{array}{c}{0}\\ \dfrac{{1}}{m}\\ {0}\end{array}\right]{f}\\ {y}=\left[\begin{array}{ccc}{1}& {0}& {0}\end{array}\right]\left[\begin{array}{c}{{x}}_{{1}}\\ {{x}}_{{2}}\\ {{x}}_{{3}}\end{array}\right]\end{array} \right. $

(22)

式中：K为液压驱动系统的等效弹簧模型的刚度系数，f 为机器人所受干摩擦阻尼 $ {f}_{\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{t}} $ ，C为液压驱动系统的等效阻尼，m为机器人与流体媒介的质量，即 ${m      = m}_{\text{liquid}}\text{+}{m}_{\text{robot}}$ ， $ {x}_{1} $ 为从动端移动距离， $ {x}_{2} $ 为从动端移动速度， $ {x}_{3} $ 为驱动端移动距离。根据上文分析，取K = 960，C = 0.11， ${m}_{x}$ = 2.0 kg， ${m}_{y}$ = 2.3 kg， ${m}_{z}$ = 2.6 kg， ${m}_{p}$ = 1.0 kg， ${m}_{r}$ = 1.5 kg， ${m}_{d}$ = 0.2 kg； $ {f}_{{\rm{r}}x} $ = 40 N， $ {f}_{{\rm{r}}y} $ = 60 N， $ {f}_{{\rm{r}}z} $ = 40 N， $ {f}_{{\rm{r}}p} $ = 20 N， $ {f}_{{\rm{r}}r} $ = 25 N， $ {f}_{{\rm{r}}d} $ = 15 N。

3 基于 $ {\mathit{H}}_{\mathbf{\infty }} $ 控制理论的机器人控制器设计

$ {H}_{\infty } $ 控制理论广泛应用在解决模型不确定性以及未知外部扰动对系统鲁棒性的影响中，从能量的角度解决了系统的镇定问题。针对忽略管路中的压力损耗导致的系统的模型不确定性问题，以及未建模的干摩擦阻尼带来的不确定死区非线性特征，本文拟采用 $ {H}_{\infty } $ 鲁棒控制器来解决。

由式(22)得出的液压驱动系统状态方程，可转化为

$ \left\{\begin{split}\dot{x}{(}{t}{) }=&{\boldsymbol{A}}{x}{(}{t}{) }+{\boldsymbol{{B}}}_{{1}}{u}{(}{t}{) }+{\boldsymbol{B}}_{{2}}{f}{(}{t}{) }\\ {y}{(}{t}{) }=&{\boldsymbol{C}}{x}{(}{t}{) }\end{split}\right. $

(23)

式中： $ \boldsymbol{A}=\left[\begin{array}{ccc}\text{0}& \text{1}& \text{0}\\ \text{−}\dfrac{\text{2}{K}}{m}& \dfrac{\text{2}{C}}{m}& \dfrac{{2K}}{m}\\ \text{0}& \text{0}& \text{0}\end{array}\right] $ ， $ {\boldsymbol{B}}_{1}=\left[\begin{array}{c}\text{0}\\ \text{−}\dfrac{\text{2}{C}}{m}\\ \text{1}\end{array}\right] $ ， ${\boldsymbol{B}}_{2}= \left[\begin{array}{c}\text{0}\\ \dfrac{\text{1}}{m}\\ \text{0}\end{array}\right]$ ， $ \boldsymbol{C}=\left[\begin{array}{ccc}\text{1}& \text{0}& \text{0}\end{array}\right] $ 。

为了让机器人能够快速、稳定地跟踪某一输入信号，引入新的状态变量 ${\boldsymbol{ \dot{x}}}_{{\rm{e}}}\left(t\right)$ 、 ${\boldsymbol{x}}_{{\rm{e}}}\left(t\right)$ ，用于描述系统跟踪误差。

$ \left\{\begin{split}{\dot{\boldsymbol{x}}}_{\text{e}}\text{(}{t}\text{) }\text=&{y}_{\text{r}}\text{(}{t}\text{)}-{y}\text{(}{t}\text{) }\\ {\boldsymbol{x}}_{\text{e}}\text{(}{t}\text{) }\text=&{\int }_{\text{0}}^{t}{(}{y}_{\text{r}}\text{(}{\tau }{) -}{y}\text{(}{\tau }\text{) ) d}{\tau }\end{split}\right. $

(24)

式中： $ {y}_{\mathrm{r}}\left(t\right) $ 是系统的期望输出，也称参考输出。因此标准线性状态空间表达式可以重写为

$ \left\{\begin{split}{\boldsymbol{\dot{ \varepsilon }}}\text{(}{t}\text{) }\text=&\bar{{\boldsymbol{A}}}{ \varepsilon }\text{(}{t}\text{) }\text+\bar{{{\boldsymbol{B}}}}_{\text{1}}{u}\text{(}{t}\text{) }\text+\bar{{\boldsymbol{B}}}_{\text{2}}{D}\text{(}{t}\text{) }\\ \boldsymbol{y}\text{(}{t}\text{) }\text=&\bar{{\boldsymbol{C}}}{ {\boldsymbol{\varepsilon}} }\text{(}{t}\text{) }\end{split}\right. $

(25)

式中：新的状态变量 ${\boldsymbol{\varepsilon}} \left(t\right) ={[\begin{array}{cc}\begin{array}{cc}{x}_{1}& {x}_{2}\end{array}& \begin{array}{cc}{x}_{3}& {x}_{{\rm{e}}}\end{array}\end{array}]}^{\mathrm{T}} $ ；D(t) 为新的输入向量，即 ${H}_{\infty } $ 控制问题中的扰动项， $ \boldsymbol{D}\left(t\right) ={[\begin{array}{cc}f& {y}_{\mathrm{r}}\left(t\right) \end{array}]}^{\mathrm{T}} $ ； $ \bar{\boldsymbol{A}} $ 、 $ {\bar{\boldsymbol{B}}}_{1} $ 、 $ {\bar{\boldsymbol{B}}}_{2} $ 、 $\bar{\boldsymbol{C}}$ 表达式为

$ \begin{split} &\bar{\boldsymbol{A}}=\left[\begin{array}{cc}\boldsymbol{A}& 0\\ -\boldsymbol{C}& 0\end{array}\right]，{\bar{\boldsymbol{B}}}_{1}=\left[\begin{array}{c}{\boldsymbol{B}}_{1}\\ 0\end{array}\right]，{\bar{\boldsymbol{B}}}_{2}=\left[\begin{array}{cc}{\boldsymbol{B}}_{2}& 0\\ 0& 1\end{array}\right] \\ &\qquad\qquad \bar{\boldsymbol{C}}=\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right] \end{split} $

取控制器的性能指标项为

$ {J=}{\int }_{{0}}^{\text{∞}}\left({{{{\boldsymbol{x}}}}_{{{\rm{e}}}}}^{\rm{{T}}}{\boldsymbol{Q}{{\boldsymbol{x}}}}_{{{\rm{e}}}}+{\boldsymbol{u}}^{\rm{{T}}}{\boldsymbol{R}}{\boldsymbol{u}}-{\textit{γ}}^{{2}}{\boldsymbol{D}}^{{{\rm{T}}}}{\boldsymbol{D}}\right) {\rm{{d}}}{t} $

(26)

式中： ${{\boldsymbol{x}}_{{\rm{e}}}}^{\mathrm{T}}{\boldsymbol{Q}}{\boldsymbol{x}}_{{\rm{e}}}$ 为系统误差的积分的平方，该项表示了外部干扰对系统跟踪性能的影响； ${\boldsymbol{u}}^{\mathrm{T}}{\boldsymbol{R}}\boldsymbol{u}$ 为系统输入的积分的平方，描述了系统输入能量的大小，是为了防止系统解算出一组过大的输入量，导致执行器负载过大。线性矩阵不等式描述为

$ \left[\begin{array}{cccc}\bar {\boldsymbol{A}} \bar{\boldsymbol{X}}\text+{\bar{\boldsymbol{B}}}_{\text{1}}\boldsymbol{V}\text+{\text{(}\bar{\boldsymbol{A}}\bar{\boldsymbol{X}}\text+{\bar{\boldsymbol{B}}}_{\text{1}}\boldsymbol{V}\text{) }}^{\text{T}}& {\bar{\boldsymbol{B}}}_{\text{2}}& \bar{\boldsymbol{X}}{\boldsymbol{M}}^{\text{T}}& {\text{(}\boldsymbol{RV}\text{) }}^{\text{T}}\\ {{\bar{\boldsymbol{B}}}_{\text{2}}}^{\text{T}}& -{\bf{1}}\rho & \text{0}& \text{0}\\ \boldsymbol{M}\bar{\boldsymbol{X}}& \text{0}& \text−{\boldsymbol{Q}}& \text{0}\\ {\boldsymbol{R}}\boldsymbol{V}& \text{0}& \text{0}& \text−{\boldsymbol{R}}\end{array}\right] < 0 $

(27)

式中： ${\bar{\boldsymbol{X}}}={\boldsymbol{X}}^{-1}$ ， ${\boldsymbol{V}={K}\boldsymbol{X}}^{-1}$ ， $ \rho ={\gamma }^{2} $ 。控制器表示为 ${\boldsymbol{u}}= \boldsymbol{V}{\bar{\boldsymbol{X}}}^{-1}x$ ，利用Matlab中LMI工具箱的mincx函数求解，解得控制器输入：

$ \boldsymbol{u}={K{\boldsymbol{x}}}{= [}\begin{array}{cccc} {-17}& {-13.1}& {-1\;535.4}& {8\;672.1} \end{array}{]}{{[}\begin{array}{cc}{{x}}_{{1}}\; {{x}}_{{2}}\;{{x}}_{{3}}\; {{x}}_{{{\rm{e}}}}\end{array}{]}}^{\rm{{T}}} $

4 仿真与实验验证 4.1 模型与控制算法验证的仿真 4.1.1 液压驱动系统数学模型的Simulink仿真验证

为了验证本文第三章提出的液压驱动系统的数学模型的可靠程度，需要对其进行基于Simulink的仿真验证。

根据液压驱动系统的动力学模型(式(7))，以及状态空间描述(式(22))，在Simulink中搭建出系统的框图，将不确定的摩擦阻尼，表示为一个在特定范围内取值的有界函数，作为系统的外部干扰项。通过力传感器测得该阻尼约为20 N，于是本文将边界设为30 N，作为未知干摩擦阻尼的一个上界。

首先验证模型的准确性，本文通过STM32单片机，读取光敏二极管模组的反馈信号，得到机器人实际的运动轨迹，并且与Simulink得到的仿真轨迹相比较。仿真结果如图6所示，其中黄线表示机器人第一次实验测得的运动轨迹；蓝线表示机器人第二次实验测得的运动轨迹；红线表示仿真得出的机器人移动轨迹。图7是在单片机的每个离散的采样周期内得到的与数学模型之间的误差，模型的最大误差不超过0.5 mm，平均误差为0.38 mm。因此可以认为，机器人动力学的仿真结果与机器人实际的运动轨迹基本吻合，该数学模型基本可靠。

4.1.2 $ {\mathit{H}}_{\mathit{\infty }} $ 轨迹跟踪控制Simulink仿真验证

在验证了动力学模型合理的情况下，本文将继续使用Simulink以验证第三章提出的 $ {H}_{\infty } $ 轨迹跟踪控制器的效果。

图8展示了 $ {H}_{\infty } $ 轨迹跟踪控制器的效果，以单位阶跃信号为例，给定机器人目标位置在x轴上距离初始坐标0.1 m。图中reference表示目标位置或者参考输入，H-infinite表示 $ {H}_{\infty } $ 轨迹跟踪控制器作用下的系统响应曲线，PID表示PID控制器作用下的系统响应曲线，其中P = 5，I = 0.5，D = 1。图9为 $ {H}_{\infty } $ 轨迹跟踪控制器的输出，即系统的输入量：滚珠丝杆滑块的速度。可以看到在 $ {H}_{\infty } $ 轨迹跟踪控制器的作用下，系统输入无高频抖动，符合执行器的实际性能。

4.2 机器人平台的运动精度实验 4.2.1 实验准备

本文通过分别测试机器人在x轴、y轴和z轴上的定位精度，来评估机器人的直线定位精度，结果如图10所示。通过编程使得机器人在纸片上对每条直线依次进行穿刺，并留下一系列针孔，再用游标卡尺测量针孔与目标直线之间的间距，以此得到机器人移动距离与期望值之间的误差。本文在机器人上固定一个陀螺仪模块，并利用它来测量俯仰轴和滚转轴的角度精度。规定姿态矢量L = $\left({x},\;{y},\;{\textit{z}},\;{\textit{θ}},\;{\textit{φ}},\;{d}\right) = \left({0}{,}\;{0}{,}\;{0}{,}\;{0}{,}\;{0}{,}\;{0}\right)$ 时，为机器人的起始位姿。

首先，制作了一个放置纸片的定位框架，如图10所示，使纸片边缘到机器人RCM点的距离是完全确定的。然后，根据机器人的Solidworks设计图，在纸片上打印若干条间隔10 mm的等距平行线，并确保第一条线处于RCM点的初始位置。最后，通过一个3D打印件，将陀螺仪模块固定在机器人的末端执行器上，用于传递末端执行器姿态。

4.2.2 实验过程

首先执行物理复位程序，STM32单片机发送复位信号，使得滚珠丝杆滑块带动连通器活塞运动，直至滚珠丝杆滑块碰到限位开关。限位开关的摆放会确保机器人复位完成后，驱动端活塞移动到行程尽头，即内部无液体的初始状态。再由操作人员在电脑上通过上位机以及串口，给单片机发送数据，控制机器人以10 mm的步距移动，直至其工作空间边界。当机器人停止运动后，操控穿刺针进行对纸片的穿刺动作，再抬起穿刺针，发送新的目标位置点，如此往复3次。通过游标卡尺测量穿刺针留下的针孔与纸片上参考直线之间的误差，并记录下实验数据。

4.2.3 实验结果与分析

图11~图14展示了机器人应用了 ${H}_{\infty } $ 轨迹跟踪算法后的实验结果，其中横坐标轴表示机器人移动的距离。

最终得到机器人系统在PID控制器下的平均线性定位误差为0.88 mm，平均角度误差为1.48°，在H_∞控制器下的平均线性定位误差为0.56 mm，平均角度误差为1.03°。

通过观察，发现当目标线定位在10~60 mm时，机器人造成的相邻穿刺点之间的距离增量实际上不是太大，基本保持在0.7 mm以内，小于编码器的分辨率。然而，机器人在10和70 mm两条直线附近的定位精度不是很理想，两个穿刺点之间的距离增量和定位精度都与其他直线有很大的偏差。经过分析，其原因是当机器人运行到工作空间边缘时，气缸内的水已基本排空，剩余的空气为可压缩流体，导致液压驱动系统运行效果不理想。

在角度精度实验研究中，随着目标角度的增大，俯仰轴和滚转轴的误差越来越大。这主要是由于在RCM机构中采用了止转棒轭来驱动。当夹角θ和φ增大时，它们对应的直线长度减小，即相同的直线误差将导致夹角误差增大。

从理论上讲，机器人俯仰轴和滚转轴的角度精度与编码器的线性精度直接相关。由图11至图12可知，机器人的平均线性定位误差分别为0.88 mm和0.56 mm，根据式(1)，理论上对应的角度误差应该分别为1.3°，0.82°，即最终测得的角度定位误差应在这两个值附近。最终在图13与图14中得到验证，实验测得的角度误差分别为1.48°和1.03°，基本合理。

5 结论

本文首先建立了一种核磁共振兼容液压驱动穿刺手术机器人的运动学模型以及动力学模型，并针对系统具有的非线性死区特性，采用 $ {H}_{\infty } $ 控制理论设计系统的状态反馈控制器，将非线性部分视为扰动输入ω。通过仿真测试，即使机器人的传感器分辨率只有0.7 mm，机器人的信号跟踪误差也可以缩小至0.02 mm。

最后，本文展示了搭建出的机器人实物，设计了基于STM32单片机的控制代码，搭建了机器人的电气系统。并且通过实验研究，验证了机器人辅助定位穿刺针的性能。机器人系统在PID控制器下的平均线性定位误差为0.88 mm，平均角度误差为1.48°，在 ${H}_{\infty } $ 控制器下的平均线性定位误差为0.56 mm，平均角度误差为1.03°

本文提出的模型是将液压驱动系统视为线性的弹簧−阻尼模型，叠加上一个未知的、未建模的非线性的干摩擦阻尼的结合。将非线性的干摩擦阻尼视为外部扰动，采用 ${H}_{\infty } $ 控制的方式控制系统完成轨迹跟踪任务。但实际测试中，机器人的摩擦阻尼较大，并不同于普通的小扰动，运用 ${H}_{\infty } $ 控制没有取得特别好的效果。因此未来还将探究关于如何处理不确定非线性的干摩擦阻尼的相关方法。