融合整體法和Simulink的人體肌肉力矩動態(tài)仿真計算

王彥鳳，王瑞，蔡玉強(qiáng)，趙欣

(1. 華北理工大學(xué) 機(jī)械工程學(xué)院，河北 唐山 063210；2.唐山市拓又達(dá)科技有限公司，河北 唐山 063000)

0 引言

穿戴式下肢外骨骼康復(fù)機(jī)器人是一種典型的外骨骼助力裝置[1]，融合了傳感、控制、信息獲取、移動計算等機(jī)器人技術(shù), 幫助操作者在無意識控制下完成站立、行走等功能[2]。人體肌肉力量作為康復(fù)機(jī)器人共享控制和自適應(yīng)控制中常用的一種參量，很多學(xué)者對肌肉力量的計算進(jìn)行了研究。婁佐政[3]根據(jù)人體解剖學(xué)建立了人體體段肌肉骨骼模型，然后據(jù)此模型在 MATLAB 開發(fā)平臺上設(shè)計了一套可視化的人體關(guān)節(jié)靜態(tài)力學(xué)評估軟件 Bio SAFE估計人體的關(guān)節(jié)力矩和肌肉力。李海龍等[4-6]通過提取患者下肢肌電信號在時域內(nèi)的特征量來測量下肢肌肉的肌肉力量。張希安等[7-9]使用骨肌仿真技術(shù)，選取 Hill 三元素肌肉模型與Lagrange 動力學(xué)模型相結(jié)合，建立骨肌系統(tǒng)模型計算下肢肌肉力矩。賈惠珍等[9]利用計算機(jī)和高速數(shù)字?jǐn)z像技術(shù)，通過量化計算產(chǎn)生力的速率(RFD)，測量等長收縮所需的時間來測定人體肌肉力量。

目前，針對肌肉力量測量研究中主要存在三方面的問題:(1)只做靜態(tài)動力學(xué)分析。由于人體步態(tài)規(guī)律變化復(fù)雜，單純的靜力學(xué)分析沒有考慮肢體慣性力、離心力等動態(tài)因素；(2)建立數(shù)學(xué)模型時將下肢簡化為多剛體的系統(tǒng)組合，沒有反應(yīng)出肌肉等柔性組織對下肢關(guān)節(jié)力矩計算的影響；(3)測量存在的輔助設(shè)備多、信號分析復(fù)雜、易受環(huán)境變化影響和反饋滯后。為了快速、準(zhǔn)確地求出下肢肌肉力量的變化規(guī)律，該項研究提出了一種融合整體法和 MATLAB /Simulink 的動態(tài)仿真計算方法。整體法主要解決2個以上物體組成的系統(tǒng)力學(xué)問題[10]，使用整體法對人體-外骨骼系統(tǒng)應(yīng)用牛頓歐拉法建立動力學(xué)方程，并加入了肢體的被動粘彈性扭矩，對模型進(jìn)行了完整的主動步態(tài)訓(xùn)練仿真計算，求解出了髖關(guān)節(jié)與膝關(guān)節(jié)的肌肉力矩，驗證了所提方法的正確性。

1 建立人體-外骨骼系統(tǒng)動力學(xué)方程

1.1 人體肌肉骨骼動態(tài)力學(xué)分析

人體下肢可看作剛性體和柔性體的組合。其中骨骼為剛性體，肌肉、關(guān)節(jié)軟骨等為柔性體。肌肉骨骼的動作過程等效如圖1所示。在使用穿戴式外骨骼機(jī)器人進(jìn)行步態(tài)訓(xùn)練時，下肢的角度變化大多集中在矢狀面內(nèi)，并且踝關(guān)節(jié)多數(shù)為被動自由度，這為問題分析提供了便利。在研究中作以下假設(shè)：(1)只考慮人體矢狀面的下肢情況；(2)股骨和脛骨為剛體，體段間由鉸鏈連接，桿件間無間隙無摩擦；(3)將腳與小腿作為一個整體。人體下肢模型的運(yùn)動方程表示為：

(1)

圖1 肌肉骨骼的動作過程

1.2 外骨骼動態(tài)力學(xué)分析

外骨骼依據(jù)人體仿生學(xué)設(shè)計，與人體生理結(jié)構(gòu)高度相似度，所以與人體有相同的運(yùn)動學(xué)特性。外骨骼機(jī)器人產(chǎn)生的驅(qū)動力矩轉(zhuǎn)化為以下兩部分：(1)克服外骨骼自身運(yùn)動中的慣性力、重力、離心力和科里奧利力；(2)補(bǔ)充給人體運(yùn)動，施加給下肢的τext。由此得到外骨骼機(jī)器人的運(yùn)動方程為

(2)

1.3 人體-外骨骼系統(tǒng)動態(tài)力學(xué)分析

整體法適用于求解系統(tǒng)內(nèi)各物體保持相對靜止或勻速直線運(yùn)動的情況，也適用于相對作勻變速直線運(yùn)動的情況[11]。外骨骼帶動人體進(jìn)行康復(fù)訓(xùn)練時二者相對靜止，而且為了使患者更舒適、安全，外骨骼作勻變速直線運(yùn)動，使用整體法分析人體-外骨骼系統(tǒng)更加便捷合理。人體-外骨骼系統(tǒng)模型如圖2所示。

圖2 人體-外骨骼系統(tǒng)模型

該模型由下肢外骨骼與人體下肢組成，各有5個參數(shù)定義:關(guān)節(jié)長度(L)、體段質(zhì)量(m)、質(zhì)心在連桿平行方向上的位置(a)、慣性矩(Iz)、關(guān)節(jié)角(θ)。系統(tǒng)的牛頓歐拉方程為：

(3)

HR指標(biāo)代表了人類-外骨骼系統(tǒng)的組合，即MHR=MH+MR,VHR=VH+VR,GHR=GH+GR具體如式(4)-式(21)所示,式中角標(biāo)1代表大腿，角標(biāo)2代表小腿。

高技術(shù)產(chǎn)業(yè)是我國經(jīng)濟(jì)的戰(zhàn)略性先導(dǎo)產(chǎn)業(yè)，對推動產(chǎn)業(yè)結(jié)構(gòu)升級和經(jīng)濟(jì)增長方式轉(zhuǎn)變具有重要的推動作用。高技術(shù)產(chǎn)業(yè)的發(fā)展水平關(guān)系到一個國家或地區(qū)的綜合競爭力，近年來我國高技術(shù)產(chǎn)業(yè)得到了較快發(fā)展，為我國的經(jīng)濟(jì)發(fā)展做出了巨大貢獻(xiàn)。但高技術(shù)產(chǎn)業(yè)快速發(fā)展的同時也存在許多問題，如高技術(shù)產(chǎn)業(yè)的區(qū)域發(fā)展不均衡，在技術(shù)水平、技術(shù)創(chuàng)新效率等方面與發(fā)達(dá)國家還存在較大差距。由于高技術(shù)產(chǎn)業(yè)具備技術(shù)和知識集中、資源消耗量少等特點，使得技術(shù)創(chuàng)新效率對產(chǎn)業(yè)發(fā)展起到更為重要的引導(dǎo)作用，也決定了各地區(qū)高技術(shù)產(chǎn)業(yè)的總體發(fā)展水平。因此，如何科學(xué)有效地對區(qū)域高技術(shù)產(chǎn)業(yè)創(chuàng)新效率進(jìn)行評價，對于區(qū)域高技術(shù)產(chǎn)業(yè)的健康發(fā)展具有重要意義。

(4)

MHR,11=JHR1+2L1XHR2cosθ2

(5)

MHR,11=JHR1+2L1XHR2cosθ2

(6)

MHR,12=MHR,21=JHR2+L1XHR2cosθ2

(7)

MHR,22=JHR2

(8)

(9)

(10)

(11)

(12)

GHR,1=g(XHR1sinθ1+XHR2sinθ12)

(13)

GHR,2=gXHR2V(14)

(15)

(16)

(17)

(18)

(19)

(20)

X2=m2a2(21)

式(16)和式(17)使用Riener的雙指數(shù)方程[12]對P(θ)進(jìn)行計算。為了簡潔起見，在式(18)～式(21)中省略了角標(biāo)H、R和HR。人體體段慣性參數(shù)為了符合中國人的特點，參照我國的成年人體生物力學(xué)參數(shù)(GB/T17245-2004)標(biāo)準(zhǔn)[13]，應(yīng)用多元線性回歸方程進(jìn)行計算。

2 建立Simulink仿真模型

Simulink提供了一個動態(tài)系統(tǒng)建模、仿真和綜合分析的集成環(huán)境，是MATLAB最重要的組件之一[14]。人體-外骨骼的動力學(xué)方程結(jié)合Simulink，動態(tài)仿真計算流程如圖3所示。未知量外骨骼機(jī)器人關(guān)節(jié)角度和關(guān)節(jié)驅(qū)動力矩由康復(fù)機(jī)器人自帶的角度編碼器和扭矩傳感器進(jìn)行計算，只需將動力學(xué)方程模塊化嵌入到仿真模型中便可實時動態(tài)計算出患者的關(guān)節(jié)肌肉力矩。建立仿真計算平臺如圖4所示。

圖3 動態(tài)仿真計算流程

圖4 仿真計算平臺

3 下肢主動康復(fù)訓(xùn)練實驗及結(jié)果分析

3.1 肌肉扭矩測量實驗過程

以某公司自主研發(fā)的自平衡下肢外骨骼康復(fù)機(jī)器人為實驗平臺，如圖5所示，該款外骨骼康復(fù)機(jī)器人的髖關(guān)節(jié)和膝關(guān)節(jié)內(nèi)置角度編碼器和扭矩傳感器。

圖5 外骨骼康復(fù)機(jī)器人

為了驗證所提方法的準(zhǔn)確性，在每位受試者的右腿臀大肌、股內(nèi)側(cè)肌、股直肌、半腱肌、腓腸肌上貼有測量電極片，用以獲得表面肌電信號，通過等距收縮實驗建立表面肌電信號(EMG)數(shù)據(jù)和關(guān)節(jié)扭矩之間的關(guān)系，以此作為對比實驗。選取3名成年男性和2名成年女性作為實驗?zāi)繕?biāo)，其中1名受試者和與他相符合的外骨骼康復(fù)機(jī)器人基本參數(shù)如表1所示。

表1 人體與外骨骼基本參數(shù)

5名受試者因患病程度不同，因此肌肉有力量也不相同。當(dāng)人產(chǎn)生輕微的邁步動作時就會被外骨骼的傳感系統(tǒng)捕捉到，控制系統(tǒng)開始逐步的增加外骨骼的驅(qū)動扭矩，直至人體-外骨骼系統(tǒng)整體關(guān)節(jié)速度達(dá)到預(yù)先的設(shè)定的正常人步態(tài)速度；對于下肢完全沒有力量的患者由醫(yī)護(hù)人員輔助其邁步。實驗過程如圖6所示。受試者每次進(jìn)行10個步幅循環(huán)，測量3次，每次間隔30 min,每次使用之前需要重新標(biāo)定測量位置。角度編碼器與外骨骼扭矩傳感器的測量數(shù)據(jù)實時傳輸至上位機(jī)，由MATLAB計算肌肉扭矩。表面肌電信號測量設(shè)備將電極采集的原始肌電圖信號經(jīng)Bagnoli 8通道系統(tǒng)放大至增益1 000，然后進(jìn)行帶通濾波，所有信號的采樣速率均為1 000 Hz。

圖6 實驗過程

3.2 實驗結(jié)果分析

將5名受試者的外骨骼角度編碼器和扭矩傳感器獲得的數(shù)據(jù)以數(shù)值矩陣格式導(dǎo)入MATLAB工作空間，取仿真時間為3 s、固定步長0.000 1。其中1名受試者的simulink仿真結(jié)果如圖7所示；表面肌電信號測得的模擬量，經(jīng)過等距肌肉收縮實驗換算轉(zhuǎn)換后的肌肉扭矩如圖8所示。

圖7 MATLAB/simulink計算肌肉扭矩

圖8 表面肌電信號測量肌肉扭矩

從圖7與圖8曲線的變化規(guī)律來看，所提方法獲得的下肢肌肉扭矩與表面肌電信號測量出的下肢肌肉扭矩十分相似。為了更加直觀的比較兩種方法,將五位受試者由simulink計算出的肌肉扭矩作為觀測值，表面肌電信號測量的肌肉扭矩作為真值，求出兩種方法獲取的肌肉扭矩的均方根誤差，其計算公式如下：

(22)

式中n為受試者人數(shù)5；Xs,i為simulink計算的第i位受試者在j時刻的肌肉扭矩；XsEMG,i

使用MATLAB對其進(jìn)行計算，髖關(guān)節(jié)肌肉扭矩的RMSE最大值為2.36，膝關(guān)節(jié)肌肉扭矩的RMSE最大值為2.89。兩種方法獲取的髖關(guān)節(jié)肌肉扭矩和膝關(guān)節(jié)肌肉扭矩偏差很小，證明所提方法能夠準(zhǔn)確測量穿戴外骨骼康復(fù)機(jī)器人進(jìn)行主動康復(fù)訓(xùn)練時的人體下肢肌肉扭矩。

4 結(jié)論

(1)提出了一種融合整體法和simulink動態(tài)仿真計算人體肌肉扭矩的方法。將進(jìn)行主動康復(fù)訓(xùn)練的人體下肢與穿戴式外骨骼康復(fù)機(jī)器人看做一個整體，并將該整體簡化為多剛體系統(tǒng)，用牛頓-歐拉法求解動力學(xué)問題，建立動力學(xué)方程。

(2)使用simulink搭建仿真計算平臺，將角度編碼器與扭矩傳感器獲得的數(shù)字量導(dǎo)入MATLAB，進(jìn)行人體肌肉扭矩計算。通過與表面肌電信號實驗結(jié)果進(jìn)行對比，驗證了所提方法的可行性。

(3)該方法不需要在受試者身上安裝傳感器，提高了測量的舒適性；信號的采集、分析和計算也比傳統(tǒng)的方法更簡便、省時。可作為穿戴式外骨骼康復(fù)機(jī)器人康復(fù)效果測評和主動訓(xùn)練時共享控制及自適應(yīng)控制的一種方法。