基于放療的六自由度醫用機械臂動力學仿真分析

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

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

引言

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

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

圖1 醫用機械臂的結構簡圖

1 醫用機械臂模型及治療環境

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

圖2 醫用機械臂

2 醫用機械臂的運動學分析

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

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

表1 醫用機械臂的D-H參數

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

2.2 醫用機械臂的軌跡規劃

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

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

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

圖4 醫用機械臂末端運行軌跡

2.3 MATLAB逆向運動求解

圖4 所示為醫用機械臂的末端運行軌跡，調用robotic toolbox 中求逆函數ikine。

T1=transl(270,0,113.8)；

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

T=ctraj(T1,T2,t)；

Q=robot.ikine(T)。

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

圖5 部分關位移仿真步數變化曲線

3 基于ADAMS的動力學分析

3.1 醫用機械臂的動力學模型

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

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

圖6 連桿受力分析

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

連桿i+1 的角速度：

連桿i+1 的角加速度：

連桿i+1 的線速度：

連桿i+1 的線加速度：

作用于連桿i+1 質心上的力：

作用于連桿i+1 質心上的力矩：

（2）向內遞推（i：0 →1）。

作用于連桿i 質心上的力：

作用于連桿i 質心上的力矩：

關節驅動力矩：

建立好機器人動力學模型，可以在已知各關節位移、速度和加速度的情況下，求得各關節所需要的驅動力和力矩，為機器人的動力學仿真提供依據。

3.2 關節剛度模型

線性彈性系統關節的力變形關系：

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

在笛卡爾空間坐標系下變形產生的能量跟關節坐標系下產生的能量相等，根據虛功原理：

根據雅可比矩陣的定義，

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

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

Kc為機器人結構變形剛度，綜上所得

同理可得：

3.3 醫用機械臂的動力學仿真分析

3.3.1 驅動函數的設置和末端載荷的施加

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

圖7 關節1的spline函數

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

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

圖9 醫用機械臂各個關節的位移仿真曲線

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

3.3.2 末端質心位移偏差

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

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

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

從表3～5 可知，在仿真剛開始階段位移偏差幾乎為0，在仿真10～20 s 內，位移偏差逐漸增大。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 加載下產生的最大位移誤差更大。位移誤差會累積傳送到腫瘤位置區域，為醫用機械臂的運動補償提供定量分析。

3.3.3 各關節扭矩曲線仿真

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

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

4 結論

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