陳繼朋，陳惠賢，楊小龍，劉曉娟，喬宇，馬利強

1. 蘭州理工大學(xué) 機電工程學(xué)院，甘肅 蘭州 730050；2. 晉中學(xué)院 機械學(xué)院，山西 晉中 030600

引言

醫(yī)用機械臂的任務(wù)是在放療過程中承載腫瘤患者，并以等中心點為目標對腫瘤治療靶區(qū)進行放療前期的擺位和放療過程中的位姿調(diào)整。與傳統(tǒng)治療床相比較，醫(yī)用機械臂具有靈活性高、工作范圍大等優(yōu)點[1]。由于腫瘤區(qū)域形狀多為復(fù)雜曲面，即醫(yī)用機械臂末端軌跡為不連續(xù)的直線或者弧線，這對機械臂的運動軌跡精度提出了更高的要求，要求醫(yī)用機械臂走出沿腫瘤區(qū)域形狀的復(fù)雜軌跡[2-3]。減少機械臂的運動誤差，提高計劃靶區(qū)（Planning Target Volume，PTV）與放療靶區(qū)（Radiotherapy Target Volume，CTV）之間的重合度，進而提高放療精度。已有研究對床板末端旋轉(zhuǎn)關(guān)節(jié)搭載力矩傳感器，根據(jù)牛頓-歐拉方程建立末端動力學(xué)方程，對運動位移誤差進行補償[4-7]。對醫(yī)用機械臂動力學(xué)研究可以解決機械臂系統(tǒng)的控制問題，實現(xiàn)預(yù)期軌跡運動和動態(tài)性能，同時為醫(yī)用機械臂的最優(yōu)化設(shè)計提供理論基礎(chǔ)[8-11]。在提升放療精度方面，腫瘤患者的擺位、圖像引導(dǎo)（Image Guide Radiotherapy，IGRT）、醫(yī)用加速器放療劑量控制等方面對放療精度都有不同程度的影響[12-17]。其中對腫瘤位置精度影響最大、研究成本較低的是醫(yī)用機械臂。

本文針對腫瘤患者載荷大負載特性，基于ADMAS 軟件對虛擬樣機進行了仿真，分析了患者重心位置變化對醫(yī)用機械臂的動力學(xué)影響。研究結(jié)果為后期醫(yī)用機械臂運動控制和位置補償提供依據(jù)，以提高醫(yī)用機械臂的運動及動力學(xué)性能，為腫瘤精準放療提供保證。醫(yī)用機械臂的結(jié)構(gòu)簡圖如圖1 所示。

圖1 醫(yī)用機械臂的結(jié)構(gòu)簡圖

1 醫(yī)用機械臂模型及治療環(huán)境

考慮精準放療所需的設(shè)備相關(guān)性能要求和約束條件，本團隊自主設(shè)計了面向重離子精準放療的醫(yī)用機械臂，同時醫(yī)用機械臂及治療環(huán)境下的模擬示意圖如圖2 所示。

圖2 醫(yī)用機械臂

2 醫(yī)用機械臂的運動學(xué)分析

2.1 基于MATLAB/robotic toolbox 機器人工具箱的建模

創(chuàng)建機器人的兩個最重要的函數(shù)是：Link 和SerialLink。對其中Link 函數(shù)的參數(shù)作簡單介紹：theta 為關(guān)節(jié)角度；d 為連桿偏移量；a 為連桿長度；alpha 為連桿扭角；sigma 表示旋轉(zhuǎn)關(guān)節(jié)為0，移動關(guān)節(jié)為1；mdh 標準的D&H 為0，否則為1；offset 表示關(guān)節(jié)變量偏移量；isprismatic 測試是否為移動關(guān)節(jié)。醫(yī)用機械臂的D-H 參數(shù)，見表1。為了減少運算量，將整體結(jié)構(gòu)縮小十倍，通過相關(guān)函數(shù)可以建立醫(yī)用機械臂的簡化模型。醫(yī)用機械臂的MATLAB 模型及編程部分程序，見圖3。

表1 醫(yī)用機械臂的D-H參數(shù)

圖3 醫(yī)用機械臂的MATLAB模型及編程部分程序

2.2 醫(yī)用機械臂的軌跡規(guī)劃

軌跡規(guī)劃既可以在笛卡爾空間坐標系下進行，也可以在關(guān)節(jié)空間坐標系進行。關(guān)節(jié)空間的軌跡規(guī)劃是以關(guān)節(jié)角度的函數(shù)來描述醫(yī)用機械臂軌跡的，也就是說醫(yī)用機械臂末端執(zhí)行器的運動軌跡是由關(guān)節(jié)變量直接確定的。在笛卡爾空間坐標系中，直接設(shè)定末端執(zhí)行器的始末位置，進而求得各關(guān)節(jié)角度的關(guān)節(jié)角位移。

本文針對腫瘤患者放療結(jié)束后，醫(yī)用機械臂從擺位位置運動到停止工作的位置進行研究。數(shù)值仿真模型仿真起點坐標為（270，0，113.8），仿真終點坐標為（170，-100，73.8）。在滿足醫(yī)用機械臂行程范圍且與其他設(shè)備沒有干涉的要求下選擇一條行程最短、能量消耗最小和機械臂空間變換最少的路徑。根據(jù)設(shè)計參數(shù)，移動關(guān)節(jié)的最大速度為v=75 mm/s?？紤]到已經(jīng)將整體結(jié)構(gòu)在MATLAB 中縮小十倍，故將速度等比例縮小。將v=75 mm/s 帶入MATLAB 仿真程序，假設(shè)整個過程運行位移為S。

取仿真時間為20 s，為了將仿真軌跡盡量平滑，將仿真步數(shù)設(shè)置為step=50，得到醫(yī)用機械臂末端軌跡如圖4 所示。

圖4 醫(yī)用機械臂末端運行軌跡

2.3 MATLAB逆向運動求解

圖4 所示為醫(yī)用機械臂的末端運行軌跡，調(diào)用robotic toolbox 中求逆函數(shù)ikine。

T1=transl(270,0,113.8)；

T2=transl(170,100,73.8)；t=50；

T=ctraj(T1,T2,t)；

Q=robot.ikine(T)。

函數(shù)進行機器人的軌跡規(guī)劃，其中T1 為末端執(zhí)行器的初始值，T2 為末端執(zhí)行器的終點值，t 為采樣步數(shù)，即得到各關(guān)節(jié)角度的變化量。打開MATLAB 左下方的工作欄，找到各個關(guān)節(jié)角度變量保存為txt 文本。截取其中關(guān)節(jié)1、2、3 的位移變化曲線如圖5 所示。

圖5 部分關(guān)位移仿真步數(shù)變化曲線

3 基于ADAMS的動力學(xué)分析

3.1 醫(yī)用機械臂的動力學(xué)模型

目前，機器人動力學(xué)的研究方法很多，如拉格朗日方法、牛頓-歐拉方法、高斯方法、凱恩方法、羅伯遜-魏登堡法、旋量方法等。其中以牛頓-歐拉法和拉格朗日法運用較多。從計算量方面考慮，使用牛頓-歐拉方法，將所有桿件的速度、加速度、慣性矩陣、質(zhì)心位置、力和力矩等都表示在各桿的坐標系中，從而使計算更加簡單，使計算關(guān)節(jié)驅(qū)動力矩的時間不僅與機器人關(guān)節(jié)數(shù)成線性比例，而且與醫(yī)用機械臂的構(gòu)型無關(guān)。

首先從連桿1 到連桿n 向外遞推計算各連桿的速度和加速度，計算出每個連桿的慣性力和力矩。第二步，從連桿n 到連桿1 向內(nèi)遞推計算各連桿內(nèi)部相互作用力和力矩、關(guān)節(jié)驅(qū)動力和力矩，見圖6。其中，xi、yi、zi分別代表第i 個關(guān)節(jié)的三個方向坐標。mig 代表第i 個連桿的重力，ni作用于連桿i 質(zhì)心上的力矩，ipi+1 為i 連桿的質(zhì)心位置，rci為i 連桿質(zhì)心相對關(guān)節(jié)i的旋轉(zhuǎn)半徑，Oi為關(guān)節(jié)i的坐標中心，fi作用于連桿i質(zhì)心上的力。

圖6 連桿受力分析

（1）向外遞推（i：0 →n-1）。

連桿i+1 的角速度：

連桿i+1 的角加速度：

連桿i+1 的線速度：

連桿i+1 的線加速度：

作用于連桿i+1 質(zhì)心上的力：

作用于連桿i+1 質(zhì)心上的力矩：

（2）向內(nèi)遞推（i：0 →1）。

作用于連桿i 質(zhì)心上的力：

作用于連桿i 質(zhì)心上的力矩：

關(guān)節(jié)驅(qū)動力矩：

建立好機器人動力學(xué)模型，可以在已知各關(guān)節(jié)位移、速度和加速度的情況下，求得各關(guān)節(jié)所需要的驅(qū)動力和力矩，為機器人的動力學(xué)仿真提供依據(jù)。

3.2 關(guān)節(jié)剛度模型

線性彈性系統(tǒng)關(guān)節(jié)的力變形關(guān)系：

其中，f 參數(shù)為6×1 的關(guān)節(jié)所受的外力或者力矩，Kx參數(shù)為6×6 的剛度矩陣，參數(shù)為6×1 位移偏差或者轉(zhuǎn)角偏差。假定機械臂連桿為剛體，關(guān)節(jié)在旋轉(zhuǎn)方向上為彈性體，關(guān)節(jié)徑向方向剛度很大。

在笛卡爾空間坐標系下變形產(chǎn)生的能量跟關(guān)節(jié)坐標系下產(chǎn)生的能量相等，根據(jù)虛功原理：

根據(jù)雅可比矩陣的定義，

把公式(13)帶入公式(12)可得：

把式(10)和式(11)帶入式(16)可得：

Kc為機器人結(jié)構(gòu)變形剛度，綜上所得

同理可得：

3.3 醫(yī)用機械臂的動力學(xué)仿真分析

3.3.1 驅(qū)動函數(shù)的設(shè)置和末端載荷的施加

對SolidWorks 中的醫(yī)用機械臂三維模型中關(guān)鍵零部件進行結(jié)構(gòu)簡化，把沒有相對運動的零件組合成組件，另存為Parasolid (x_t)格式導(dǎo)入ADAMS 中。在ADAMS 中重新定義各零部件的材料屬性，采用輕質(zhì)鋁合金，確定密度、彈性模量、剛度等參數(shù)。在關(guān)節(jié)模型處添加運動副、固定副、maker 點，得到醫(yī)用機械臂的動力學(xué)模型。利用前一章中Robotics Toolbox 軌跡規(guī)劃得到的各關(guān)節(jié)角位移，在MATLAB 中已經(jīng)將各桿件長度縮小十倍，故應(yīng)將移動關(guān)節(jié)變化數(shù)值放大十倍，轉(zhuǎn)動關(guān)節(jié)角度數(shù)值不變。在ADAMS中通過Data Elements 功能建立各關(guān)節(jié)樣條曲線SPLINE，對應(yīng)的驅(qū)動函數(shù)調(diào)用格式為CUBSPL ( time, 0, SPLINE_number, 0)，最后將六個關(guān)節(jié)樣條驅(qū)動函數(shù)依次施加到相應(yīng)的運動副上。其中關(guān)節(jié)1 的樣條曲線如圖7 所示。

圖7 關(guān)節(jié)1的spline函數(shù)

參考中國正常人的標準體重表，假定腫瘤患者體重為1000 N?？紤]到腫瘤患者的重心位置的不同，在醫(yī)用機械臂床板中心建立一個薄圓槽，沿著床板x 方向距離圓槽兩側(cè)距離60 cm 和40 cm 處分別建立兩個圓槽，方便建立marker 點，床板上從左到右的依次為marker9、marker10和marker11。模擬腫瘤患者重心位置。針對其中任意一個marker 點施加恒定載荷1000 N，仿真時間設(shè)為20 s，采樣步長設(shè)為0.04 s。選取末端床板質(zhì)量中心點作為標記點，記錄各關(guān)節(jié)的位移、速度、加速度變化曲線，各關(guān)節(jié)中心的位移變化曲線如圖8 和圖9 所示。分析對比不同marker 點施加1000 N 時末端床板的位移變化曲線。

圖8 ADAMS中醫(yī)用機械臂上的三個圓孔槽

圖9 醫(yī)用機械臂各個關(guān)節(jié)的位移仿真曲線

從仿真結(jié)果可知，各關(guān)節(jié)位移變化曲線都比較平緩，沒有出現(xiàn)陡峰和拐點。軌跡規(guī)劃中，醫(yī)用機械臂首先進行水平面的轉(zhuǎn)動，然后才有豎直方向的移動。在5～10 s 運動期間，大臂、小臂及末端床板的位移變化明顯；運動時間在15～20 s 之間，位移變化較為明顯。

3.3.2 末端質(zhì)心位移偏差

僅列出對末端marker11 點施加1000 N 時醫(yī)用機械臂仿真軌跡中位移、速度、加速度的變化曲線，其曲線分別對應(yīng)圖10 a、b 和c。

圖10 Marker11點加載的位移、速度、加速度曲線

由圖10 b 和c 可知，醫(yī)用機械臂在整個運動過程中末端速度、加速度變化平緩，沒有突變，且都符合技術(shù)規(guī)范中對速度、加速度的要求。在ADAMS 中空載運行前一章得到MATLAB 仿真軌跡，分別記錄末端床板在0、2.5、5、10、15、20 s 時的x、y、z 方向的位移，與之前marker9、marker10、marker11 三點施加1000 N 時末端床板的位移曲線進行對比。

從表3～5 可知，在仿真剛開始階段位移偏差幾乎為0，在仿真10～20 s 內(nèi)，位移偏差逐漸增大。Marker9 點y 方向的最大位移誤差為2.078 mm，在x 方向和z 方向的最大位移誤差0.82 mm 和1.07 mm。marker11 點在y 方向的最大位移誤差1.889 mm，在x 方向和z 方向最大位移誤差0.81 mm和0.97 mm。在y 方向的位移偏差比x、z 方向的位移偏差都大，這與患者的重心方向相重合的原因相吻合，對x 方向的位移偏差影響最小。隨著運行時間的增加，運動誤差有增大的趨勢。

表3 Marker9點的位移偏差表

表4 Marker10點的位移偏差表

表5 Marker11點的位移偏差表

Marker9 比marker11 離床板末端位置中心更遠，marker9 加載下產(chǎn)生的最大位移誤差更大。位移誤差會累積傳送到腫瘤位置區(qū)域，為醫(yī)用機械臂的運動補償提供定量分析。

3.3.3 各關(guān)節(jié)扭矩曲線仿真

這里僅列出marker10 和marker11 點加載下的力矩曲線圖，見圖11。從圖11 可知在y 方向的最大力矩大于在x 和z 方向的最大力矩。Marker11 點加載下的力矩在大于marker10 點加載下的力矩。在擺位過程中，應(yīng)該盡量將腫瘤患者擺放在床板等中心點，提高放療精度。

圖11 不同加載點的力矩曲線

4 結(jié)論

本文基于醫(yī)用機械臂特定軌跡，結(jié)合機器人運動學(xué)、動力學(xué)仿真常用的三個軟件SolidWorks、MATLAB 和ADAMS，運用聯(lián)合仿真的方法，對醫(yī)用機械臂剛體多體醫(yī)用機械臂動力學(xué)模型，并進行運動學(xué)、動力學(xué)仿真分析，驗證了腫瘤患者在床板位置中心不同存在位移偏差，得到了醫(yī)用機械臂末端位移偏差曲線和力矩曲線，為RV 減速器的選型和醫(yī)用機械臂的運動補償提供依據(jù)。這種聯(lián)合分析的方法，結(jié)合了SolidWorks 強大的繪圖、MATLAB 的運動學(xué)分析以及ADAMS 的動力學(xué)分析的優(yōu)勢，克服了單獨使用的局限性，提高了在加載狀態(tài)下醫(yī)用機械臂動態(tài)響應(yīng)的準確性。精準的模擬了醫(yī)用機械臂的運動狀態(tài)，為醫(yī)用機械臂的運動控制和位移補償提供了理論依據(jù)。