陳江城,張小棟,李睿,石強(qiáng)勇,王賀

(西安交通大學(xué)現(xiàn)代設(shè)計(jì)及轉(zhuǎn)子軸承系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室, 710049, 西安)

?

利用表面肌電信號(hào)的下肢動(dòng)態(tài)關(guān)節(jié)力矩預(yù)測(cè)模型

陳江城,張小棟,李睿,石強(qiáng)勇,王賀

(西安交通大學(xué)現(xiàn)代設(shè)計(jì)及轉(zhuǎn)子軸承系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室, 710049, 西安)

為實(shí)現(xiàn)表面肌電信號(hào)的下肢關(guān)節(jié)力矩動(dòng)態(tài)解碼,建立了從表面肌電信號(hào)到關(guān)節(jié)力矩輸出的人體下肢運(yùn)動(dòng)系統(tǒng)正向生物力學(xué)模型。首先,從幅值和頻率兩個(gè)角度建立表面肌電信號(hào)到骨骼肌激活程度模型;其次,根據(jù)肌絲滑移理論,構(gòu)建反映骨骼肌生理結(jié)構(gòu)和微觀力學(xué)特性的肌肉力模型,同時(shí)確定活動(dòng)肌肉拉力線方向及力作用點(diǎn)位移矢量,將骨骼肌力轉(zhuǎn)換到關(guān)節(jié)力矩;最后,以牛頓-歐拉逆動(dòng)力學(xué)方法獲得關(guān)節(jié)力矩作為準(zhǔn)確值,給出正向生物力學(xué)模型參數(shù)動(dòng)態(tài)標(biāo)定方法。在模型基礎(chǔ)上,對(duì)4名對(duì)象進(jìn)行隨意步態(tài)下膝關(guān)節(jié)屈伸動(dòng)態(tài)力矩預(yù)測(cè)試驗(yàn),結(jié)果表明:所建模型對(duì)步態(tài)行走下的膝關(guān)節(jié)動(dòng)態(tài)關(guān)節(jié)力矩具有很好的動(dòng)態(tài)跟蹤性能,最大絕對(duì)誤差為(11.0±1.32) N·m,平均殘差為(4.43±0.698) N·m,預(yù)測(cè)值與準(zhǔn)確值之間的平均線性相關(guān)系數(shù)為0.927±0.042,驗(yàn)證了該方法的正確性和有效性;可為康復(fù)訓(xùn)練機(jī)器人人機(jī)協(xié)同過程中的力學(xué)交互模式研究提供接口。

表面肌電信號(hào);關(guān)節(jié)力矩預(yù)測(cè);肌肉模型;正向生物力學(xué)

運(yùn)動(dòng)訓(xùn)練對(duì)中風(fēng)偏癱患者的運(yùn)動(dòng)功能恢復(fù)十分重要,康復(fù)訓(xùn)練機(jī)器人具有比人工醫(yī)療師輔助訓(xùn)練更多的優(yōu)點(diǎn)而被研究者和康復(fù)工作者所認(rèn)可,目前世界上典型的下肢康復(fù)訓(xùn)練機(jī)器人有LOPES、LOKOMAT、Gait Trainer和Haptic Walker等[1]?？祻?fù)訓(xùn)練過程一般分為被動(dòng)階段和主動(dòng)階段,主動(dòng)階段注重患者的主動(dòng)參與和自身能力,選擇性地給患者提供助力[2-3]。實(shí)際上,康復(fù)機(jī)器人輔助的根本作用是,在患者自身力量無(wú)法提供目標(biāo)運(yùn)動(dòng)所需力的時(shí)候,給患者額外提供部分需要的力,通過訓(xùn)練逐漸擺脫機(jī)器人的幫助,達(dá)到人進(jìn)機(jī)退,同時(shí)應(yīng)保證人、機(jī)力交互的柔順性[4]。因此,康復(fù)機(jī)器人需要對(duì)患者實(shí)際運(yùn)動(dòng)能力進(jìn)行實(shí)時(shí)感知,但是目前由于缺少對(duì)患者運(yùn)動(dòng)能力的直接感知手段,從力交互機(jī)理出發(fā)設(shè)計(jì)的康復(fù)控制策略很少,制約康復(fù)機(jī)器人的進(jìn)一步提高。

骨骼肌是人體運(yùn)動(dòng)的主要?jiǎng)恿υ?在神經(jīng)系統(tǒng)調(diào)控下收縮,牽引骨骼和關(guān)節(jié)運(yùn)動(dòng),因此通過運(yùn)動(dòng)中骨骼肌收縮力檢測(cè)可以感知人體的運(yùn)動(dòng)能力。目前,直接測(cè)量體內(nèi)肌肉力大小是不現(xiàn)實(shí)的,表面肌電(sEMG)蘊(yùn)含肌肉活動(dòng)強(qiáng)度、活動(dòng)時(shí)間等信息,成為體表無(wú)創(chuàng)檢測(cè)肌肉活動(dòng)的重要方法,被廣泛用于康復(fù)領(lǐng)域。文獻(xiàn)[5]將表面肌電信號(hào)的強(qiáng)度與關(guān)節(jié)力矩等效成一階線性關(guān)系,通過表面肌電信號(hào)估計(jì)關(guān)節(jié)力矩,并控制助行外骨骼機(jī)器人,增強(qiáng)了穿戴者的行走能力。文獻(xiàn)[6]通過Hill肌肉模型建立表面肌電信號(hào)與關(guān)節(jié)力矩的關(guān)系,用于控制外骨骼機(jī)器人。文獻(xiàn)[7]通過收集特定肌肉上的肌電信號(hào)并進(jìn)行處理,當(dāng)被處理的信號(hào)超過一定閾值時(shí)觸發(fā)輔助,作為給患者的獎(jiǎng)勵(lì),相同的肌電生物反饋方法也被用于上肢的康復(fù)訓(xùn)練[8-9]。另外,文獻(xiàn)[10]通過神經(jīng)網(wǎng)絡(luò)學(xué)習(xí),建立了表面肌電信號(hào)與關(guān)節(jié)力矩的關(guān)系,用于上肢假肢控制;文獻(xiàn)[11-12]都基于Hill肌肉模型建立了幅值大小到關(guān)節(jié)力矩的映射模型。

綜上可以看出,通過肌電感知肌肉力和關(guān)節(jié)力矩是一種可行的方法。然而研究表明,產(chǎn)生相同的肌肉力在不同的運(yùn)動(dòng)任務(wù)下肌電活動(dòng)模式是不一樣的,而且不是簡(jiǎn)單的線性關(guān)系,若要保證預(yù)測(cè)模型的精度及在任意運(yùn)動(dòng)下都有效,所建立的模型必須能反映運(yùn)動(dòng)系統(tǒng)的解剖結(jié)構(gòu)和生理特點(diǎn)[11]。Hill模型是在離體肌肉實(shí)驗(yàn)基礎(chǔ)上提出的[13],沒有考慮肌肉收縮的微觀特性,也不能描述活體肌肉完整生物力學(xué)過程。文獻(xiàn)[14]基于肌絲滑移理論,從微觀角度研究肌肉的力學(xué)特性,從頻率角度分析表面肌電信號(hào)與骨骼肌力之間的關(guān)系,重點(diǎn)討論了基于表面肌電的肌力預(yù)測(cè)結(jié)果,沒有對(duì)關(guān)節(jié)力矩的預(yù)測(cè)進(jìn)行研究。

本文以提高下肢康復(fù)機(jī)器人的運(yùn)動(dòng)感知能力為目標(biāo),考慮下肢系統(tǒng)運(yùn)動(dòng)生物力學(xué)模型的完整性和動(dòng)態(tài)力矩預(yù)測(cè)的準(zhǔn)確性,對(duì)骨骼肌生物電調(diào)控(激勵(lì))、收縮力學(xué)模型以及肌肉力到關(guān)節(jié)力矩等過程進(jìn)行深入研究,從運(yùn)動(dòng)機(jī)理出發(fā),構(gòu)建了下肢運(yùn)動(dòng)系統(tǒng)正向生物力學(xué)全過程模型,實(shí)現(xiàn)了下肢運(yùn)動(dòng)力學(xué)信息的表面肌電信號(hào)動(dòng)態(tài)解碼,為下肢康復(fù)機(jī)器人力交互機(jī)理及輔助控制策略研究提供合適的人機(jī)交互接口。

1 正向生物力學(xué)模型

建立正向生物力學(xué)模型的最終目標(biāo)是:通過輸入原始表面肌電信號(hào)和關(guān)節(jié)運(yùn)動(dòng)學(xué)數(shù)據(jù)得到關(guān)節(jié)力矩的大小,且使模型能夠適用于較寬的運(yùn)動(dòng)任務(wù)(除發(fā)生痙攣等意外)和任意下肢關(guān)節(jié)。人體下肢運(yùn)動(dòng)系統(tǒng)主要包括神經(jīng)、肌肉和骨骼關(guān)節(jié)3個(gè)部分,分別發(fā)揮控制、驅(qū)動(dòng)和執(zhí)行的作用。根據(jù)其力學(xué)過程,本文分別建立肌肉生物電激活模型、肌肉收縮力學(xué)模型及肌肉功能模型,3者按照力學(xué)關(guān)系聯(lián)合形成完整的正向生物力學(xué)模型,如圖1所示。其中,肌肉生物電激活模型的作用是將表面肌電信號(hào)sEMG表征為肌肉激活程度a(t),肌肉收縮生物力學(xué)模型描述肌肉在激活程度a(t)下的輸出力Fm(t)的情況,肌肉功能模型作用是將相關(guān)聯(lián)肌肉力合成得到關(guān)節(jié)力矩Tj(t)輸出。

圖1 正向生物力學(xué)模型建立過程

1.1 肌肉生物電激活模型

人體骨骼肌由大量肌纖維組成,肌纖維的活動(dòng)受到神經(jīng)系統(tǒng)控制,一個(gè)神經(jīng)元和受其支配的所有肌纖維構(gòu)成一個(gè)運(yùn)動(dòng)單位。肌纖維的收縮是集力學(xué)、生理學(xué)和生物化學(xué)的一個(gè)復(fù)雜過程。在神經(jīng)沖動(dòng)的激勵(lì)下,肌細(xì)胞膜激發(fā)動(dòng)作電位,一方面疊加形成皮膚表面肌電信號(hào),另一方面興奮整個(gè)細(xì)胞,經(jīng)過興奮收縮耦聯(lián)引起肌絲滑移。因此,表面肌電信號(hào)表征了肌肉的激活強(qiáng)度,其按照以下3個(gè)步驟形成肌肉生物電激活模型。

(1)對(duì)于單根肌纖維,當(dāng)動(dòng)作電位幅值超過門限值,會(huì)引起肌漿網(wǎng)Ca2+通道開啟,使Ca2+進(jìn)入肌漿,當(dāng)動(dòng)作電位幅值低于門限值,Ca2+又迅速回流,肌漿內(nèi)的Ca2+與肌纖維中的肌鈣蛋白作用引起肌絲滑移,Ca2+濃度決定著肌纖維收縮強(qiáng)度。肌纖維激發(fā)動(dòng)作電位存在“全或無(wú)”的特點(diǎn),即只要激發(fā),其波形幅值變化固定,所以肌漿內(nèi)的Ca2+濃度是通過動(dòng)作電位發(fā)放率調(diào)控的。文獻(xiàn)[15]從肌纖維的生理結(jié)構(gòu)和運(yùn)動(dòng)機(jī)理出發(fā),建立了肌漿內(nèi)平均Ca2+濃度μ與動(dòng)作電位發(fā)放率的數(shù)學(xué)模型,并通過Sigmoid函數(shù)擬合

(1)

式中:fa為動(dòng)作電位發(fā)放率;a1、b1和c1為擬合系數(shù)。另外,μ與肌纖維激活程度β的關(guān)系可以表示為

(2)

式中:K1=8.72,表示Ca2+與肌鈣蛋白弱作用時(shí)兩者分離速率與結(jié)合速率之比;K2=0.194,表示Ca2+與肌鈣蛋白強(qiáng)作用時(shí)兩者分離速率與結(jié)合速率之比。

以往研究發(fā)現(xiàn),肌纖維的動(dòng)作電位發(fā)放率在8～50 Hz之間[16],而表面肌電信號(hào)的有效頻率范圍是在20～500 Hz之間。本文通過線性模型建立兩者之間的映射關(guān)系,以表面肌電信號(hào)短時(shí)傅里葉變換后得到中值頻率fm表示表面肌電信號(hào)頻率,線性變換得到肌纖維動(dòng)作電位發(fā)放率,即

(3)

式中:fi(t)為Ai(t)時(shí)刻短時(shí)傅里葉變換的各頻率;Ai(t)為各頻率點(diǎn)對(duì)應(yīng)的幅值;N為頻率點(diǎn)數(shù),由頻率分辨率,即短時(shí)傅里葉變換所用窗函數(shù)寬度決定。至此,式(1)～式(3)描述了表面肌電信號(hào)頻率與單根肌纖維激活程度的關(guān)系。

(2)對(duì)于整個(gè)肌肉,其激活程度還與肌纖維的募集數(shù)量有關(guān),肌纖維募集越多,就有更多的肌膜動(dòng)作電位疊加形成表面肌電信號(hào),因此可用表面肌電信號(hào)的幅值反映肌纖維募集數(shù)量。對(duì)于原始表面肌電信號(hào),通過零滯后的四階巴斯沃特濾波器進(jìn)行濾波,以達(dá)到消除運(yùn)動(dòng)偽跡的目的,截止頻率為20 Hz。在此基礎(chǔ)上,進(jìn)行全波整流和4 Hz低通濾波得到時(shí)域信號(hào)R(t),表示表面肌電信號(hào)活動(dòng)強(qiáng)度。用相同的方法對(duì)每塊肌肉在最大主動(dòng)收縮力下測(cè)得的表面肌電信號(hào)進(jìn)行處理,得到時(shí)域信號(hào)RMVC(t),取其中最大值作為標(biāo)準(zhǔn),對(duì)所有試驗(yàn)下處理后的肌電信號(hào)進(jìn)行歸一化處理,轉(zhuǎn)化為肌纖維募集程度,即

(4)

(3)任意骨骼肌的激活程度可以表示為肌纖維激活程度與肌纖維募集程度的乘積

(5)

利用上述模型,由某對(duì)象股直肌表面肌電信號(hào)計(jì)算肌肉激活程度,如圖2所示。從圖可知,同時(shí)考慮幅值與頻率兩個(gè)因素時(shí)獲得的肌肉激活程度a與單獨(dú)考慮幅值時(shí)的肌肉激活程度u相比,變化趨勢(shì)微小,但是幅值上變化較為明顯。激活程度的改變將影響后續(xù)模型參數(shù)的準(zhǔn)確性。

圖2 肌肉激活程度變化曲線

1.2 肌肉收縮力學(xué)模型

目前,對(duì)肌肉生物力學(xué)研究多數(shù)是在Hill提出的肌肉力學(xué)模型基礎(chǔ)上進(jìn)行的,其肌肉力方程為

(6)

式中:Fmmax為肌肉最大等長(zhǎng)收縮力;l為肌肉纖維長(zhǎng)度;v為肌肉收縮速度;f(l)、f(v)和fp(l)為歸一化的肌纖維力學(xué)特性曲線,是通過離體實(shí)驗(yàn)測(cè)量得到的,沒有考慮微觀特性,因此無(wú)法有效反映活體肌肉生理參數(shù)對(duì)生物力學(xué)的影響。

從微觀上看,骨骼肌收縮是由于大量肌球蛋白分子馬達(dá)對(duì)細(xì)肌絲集體做功引起粗細(xì)肌絲相對(duì)滑移導(dǎo)致的?；诩〗z滑移理論,文獻(xiàn)[17]提出了肌肉橫橋動(dòng)力學(xué)模型,推出了肌絲滑動(dòng)速率方程和肌絲橫截面上的主動(dòng)力力學(xué)方程。文獻(xiàn)[18-19]通過二態(tài)布朗棘輪模型及統(tǒng)計(jì)動(dòng)力學(xué)方法,研究了分子馬達(dá)的集體運(yùn)行特性。文獻(xiàn)[20]針對(duì)肌纖維結(jié)構(gòu)特點(diǎn)和肌球蛋白馬達(dá)運(yùn)行機(jī)制,利用非平衡統(tǒng)計(jì)力學(xué)的方法,構(gòu)建了肌節(jié)主動(dòng)收縮力學(xué)模型,該模型可以用于定量討論負(fù)載力、Ca2+及能量因子三磷酸腺苷濃度等對(duì)收縮力學(xué)特性的影響。

肌節(jié)是肌肉收縮最小功能單元,骨骼肌是由大量肌節(jié)串并聯(lián)而成,而肌肉整體收縮力大小與肌節(jié)并聯(lián)數(shù)目相關(guān),任意骨骼肌的主動(dòng)收縮力學(xué)模型可以表示為

(7)

(8)

另外,根據(jù)文獻(xiàn)[21]對(duì)肌節(jié)主動(dòng)力和肌節(jié)長(zhǎng)度關(guān)系的實(shí)驗(yàn)研究表明,主動(dòng)力隨著肌節(jié)長(zhǎng)度增加先增后減,本文用高斯函數(shù)模擬粗細(xì)肌絲的重疊程度α,即

(9)

式中:lopti為肌肉最大等長(zhǎng)收縮時(shí)肌肉纖維長(zhǎng)度。

除主動(dòng)收縮力外,肌肉組織本身還存在被動(dòng)拉伸的特性,采用剛度阻尼模型表示為

(10)

式中:km表示肌肉剛度;γm為肌肉阻尼;lo為肌纖維初始長(zhǎng)度;Δlm為肌肉伸長(zhǎng)量。

因此,肌肉收縮力學(xué)模型為

Fm=Fa+Fp

(11)

從肌肉收縮力學(xué)模型可知,除激活程度a外,肌肉截面積Sm、最優(yōu)長(zhǎng)度lopti、剛度km及阻尼γm等結(jié)構(gòu)和物理參數(shù)都影響肌肉收縮力的大小。

1.3 肌肉功能模型

根據(jù)肌肉結(jié)構(gòu)及力學(xué)原理,肌肉收縮力在關(guān)節(jié)上產(chǎn)生的轉(zhuǎn)動(dòng)力矩為

T=r×Fcosφ

(12)

式中:r表示關(guān)節(jié)轉(zhuǎn)動(dòng)中心到力作用點(diǎn)的位移矢量;F表示肌肉力矢量;φ表示肌肉羽狀角。因此,要將肌肉力轉(zhuǎn)換到關(guān)節(jié)力矩的大小,需要進(jìn)一步確定肌肉拉力線方向矢量I和位移矢量r。

肌肉附著點(diǎn)位置,包括起點(diǎn)、代起點(diǎn)、止點(diǎn)和代止點(diǎn),是確定肌肉長(zhǎng)度、肌肉伸縮速度、肌肉拉力線和肌肉力作用點(diǎn)的基礎(chǔ)。文獻(xiàn)[22]通過大量尸體標(biāo)本檢測(cè)確定人體下肢肌肉附著點(diǎn)位置和人體下肢骨形參數(shù),并以3個(gè)關(guān)節(jié)中心為坐標(biāo)基準(zhǔn)點(diǎn),建立骨盆、大腿、小腿局部坐標(biāo)系,給出了用于計(jì)算活體肌肉附著點(diǎn)三維坐標(biāo)的線性回歸方程。本文基于文獻(xiàn)[22]研究結(jié)論,給出下肢位姿動(dòng)態(tài)變化過程中的肌肉拉力線和力位移矢量的確定方法,具體步驟如下。

第1步,利用光學(xué)運(yùn)動(dòng)捕捉系統(tǒng)記錄人體下肢關(guān)鍵點(diǎn)的位置,確定髖、膝和踝關(guān)節(jié)中心的位置坐標(biāo)。對(duì)于膝關(guān)節(jié)和踝關(guān)節(jié),以關(guān)節(jié)內(nèi)外側(cè)標(biāo)記點(diǎn)位置的中心作為關(guān)節(jié)中心。采用文獻(xiàn)[23]提出的方法確定髖關(guān)節(jié)中心:首先,以兩髂前上棘中心為基準(zhǔn)點(diǎn),根據(jù)左右髂前上棘和髂后上棘位置坐標(biāo),建立骨盆坐標(biāo)系,如圖3所示,其坐標(biāo)軸方向確定方法如下

(13)

(14)

U=V×W

(15)

其次,計(jì)算髖關(guān)節(jié)中心在骨盆坐標(biāo)系中的位置,并轉(zhuǎn)換至全局坐標(biāo)系下。同時(shí),根據(jù)文獻(xiàn)[24]給出的線性回歸方程確定髖關(guān)節(jié)中心在骨盆坐標(biāo)系中的位置,方程如下

(16)

式中:Pd代表骨盆深度,為兩髂前上棘中點(diǎn)與兩髂后上棘中點(diǎn)之間的距離,mm;Pw代表骨盆寬度,為兩髂前上棘之間的距離,mm。

圖3 骨盆坐標(biāo)系及與髖關(guān)節(jié)中心的關(guān)系

第2步,從骨盆、大腿、小腿分別選擇3個(gè)標(biāo)記點(diǎn),以關(guān)節(jié)中心為基準(zhǔn)點(diǎn)建立各肢體的局部坐標(biāo)系,然后根據(jù)文獻(xiàn)[22]給出的肌肉附著點(diǎn)三維坐標(biāo)計(jì)算線性回歸方程,計(jì)算附著點(diǎn)在局部坐標(biāo)系中的位置,同時(shí)通過坐標(biāo)變換,將附著點(diǎn)位置轉(zhuǎn)換到全局坐標(biāo)系下表示。

第3步,利用肌肉附著點(diǎn)以及關(guān)節(jié)中心坐標(biāo),確定在全局坐標(biāo)下的肌肉拉力線方向矢量和力位移矢量r,同時(shí)計(jì)算肌肉長(zhǎng)度lm、肌肉纖維長(zhǎng)度l以及肌肉伸縮速度v。

在確定所有肌肉力方向矢量和力相對(duì)于關(guān)節(jié)中心的位移矢量后,分別選出對(duì)下肢各關(guān)節(jié)運(yùn)動(dòng)起作用的肌肉,計(jì)算每塊肌肉對(duì)關(guān)節(jié)的力矩貢獻(xiàn)并相加,得到該關(guān)節(jié)力矩,即

(17)

式中:j代表關(guān)節(jié);N代表作用于該關(guān)節(jié)上的肌肉數(shù)量;φi表示第i塊肌肉的羽狀角。

2 模型的標(biāo)定

從建立的下肢運(yùn)動(dòng)系統(tǒng)正向生物力學(xué)模型看出,模型中包括了反映肌肉生理特性及結(jié)構(gòu)的參數(shù),不同肌肉具有不同的參數(shù)大小,如橫截面積較大的肌肉并聯(lián)更多的肌節(jié),產(chǎn)生的最大肌肉力也更大。因此,在模型應(yīng)用于具體對(duì)象之前,需要對(duì)模型參數(shù)進(jìn)行調(diào)整,即模型標(biāo)定。

2.1 模型標(biāo)定原理

模型標(biāo)定需要一個(gè)標(biāo)準(zhǔn)的參考值?？紤]測(cè)量人體活體肌肉在關(guān)節(jié)上的作用力不現(xiàn)實(shí),本文利用足底力(GRF)和運(yùn)動(dòng)學(xué)數(shù)據(jù),通過牛頓-歐拉逆動(dòng)力學(xué)方程[25]計(jì)算各關(guān)節(jié)的輸出力矩Tjiv,并將其作為模型標(biāo)定的參考值。

模型標(biāo)定的目標(biāo)是通過調(diào)整正向生物力學(xué)模型參數(shù),最終尋找一組合適的參數(shù)值,使模型輸出的關(guān)節(jié)力矩Tj與標(biāo)準(zhǔn)參考值無(wú)限接近。本文采用模型參數(shù)優(yōu)化算法實(shí)現(xiàn)這一目標(biāo),并定義優(yōu)化目標(biāo)函數(shù)如下

(18)

式中:n為采樣點(diǎn)數(shù)。整個(gè)模型的標(biāo)定過程如圖4所示。

圖4 模型標(biāo)定過程

另外,為了使標(biāo)定后模型能夠適用于較寬的運(yùn)動(dòng)任務(wù)條件,優(yōu)化函數(shù)中所選取的樣本應(yīng)盡可能涵蓋復(fù)雜的運(yùn)動(dòng)情況,當(dāng)然這也會(huì)增加模型參數(shù)尋找的時(shí)間。鑒于本文的目的是為康復(fù)機(jī)器人提供患者運(yùn)動(dòng)能力感知的功能,而康復(fù)的主要運(yùn)動(dòng)訓(xùn)練是步態(tài)行走,文中利用不同步行速度下的若干步態(tài)周期數(shù)據(jù)樣本對(duì)模型參數(shù)進(jìn)行標(biāo)定,使標(biāo)定后模型能夠完成步態(tài)行走下的關(guān)節(jié)力矩準(zhǔn)確預(yù)測(cè)。

2.2 模型參數(shù)的優(yōu)化算法

為了實(shí)現(xiàn)目標(biāo)函數(shù)的最小化,采用遺傳優(yōu)化算法尋找模型的最優(yōu)參數(shù)值。模型中優(yōu)化的參數(shù)包括肌肉收縮力學(xué)模型中的高斯函數(shù)標(biāo)準(zhǔn)差σ、肌肉最大力Famax、最優(yōu)長(zhǎng)度lopti以及剛度km和阻尼γm。對(duì)于肌肉的羽狀角φ,文獻(xiàn)[11]研究表明其值變化范圍很小,因此本文采用該文獻(xiàn)提供的羽狀角數(shù)據(jù),不對(duì)其進(jìn)行優(yōu)化。

在對(duì)模型參數(shù)進(jìn)行優(yōu)化時(shí),初始種群的每個(gè)個(gè)體代表模型參數(shù)的一組初始值,進(jìn)化方向?yàn)槭鼓繕?biāo)函數(shù)值減小。另外,本文采用文獻(xiàn)[26]提供的最大肌肉力和肌肉最優(yōu)長(zhǎng)度作為參考初始值,優(yōu)化范圍分別被約束在初始值上下±50%和±20%之間。這樣一方面減小了優(yōu)化搜索空間,加快優(yōu)化速度,另一方面也能防止大幅度偏離實(shí)際情況而與解剖學(xué)知識(shí)不符的參數(shù)值出現(xiàn)。

3 膝關(guān)節(jié)動(dòng)態(tài)力矩預(yù)測(cè)試驗(yàn)

為了驗(yàn)證模型的可靠性和有效性,選取4名男性對(duì)象A、B、C和D(平均年齡(29±5.5)歲,平均身高(172±4.3) cm,平均體質(zhì)量(63.9±5.0) kg)進(jìn)行生物力學(xué)模型標(biāo)定和隨意步態(tài)下的膝關(guān)節(jié)屈伸動(dòng)態(tài)力矩預(yù)測(cè)試驗(yàn)。根據(jù)生理學(xué)先驗(yàn)知識(shí),確定10塊對(duì)膝關(guān)節(jié)運(yùn)動(dòng)有貢獻(xiàn)的主要肌肉,分別是半膜肌(SM)、股二頭肌長(zhǎng)頭(BFL)、股二頭肌短頭(BFS)、股內(nèi)側(cè)肌(VM)、股外側(cè)肌(VL)、股直肌(RF)、腓腸肌內(nèi)側(cè)(MG)、腓腸肌外側(cè)(LG)、闊筋膜張肌(TFL)和縫匠肌(SR)。實(shí)驗(yàn)中,所有受試者完成多種步速(約0.8、1.0和1.2 m/s)行走,記錄相應(yīng)的步態(tài)數(shù)據(jù),試驗(yàn)設(shè)備場(chǎng)景如圖5所示。下肢運(yùn)動(dòng)學(xué)數(shù)據(jù)通過英國(guó)Oxford Metrics Limited公司生產(chǎn)的VICON運(yùn)動(dòng)光學(xué)捕捉系統(tǒng)記錄,包括12臺(tái)高精度攝像機(jī),采樣頻率為50 Hz;足底力信號(hào)利用AMTI測(cè)力板獲取,采樣率為1 000 Hz;表面肌電信號(hào)由16通道多導(dǎo)生理記錄儀記錄,采樣頻率為1 024 Hz。

完成數(shù)據(jù)獲取后,對(duì)不同信號(hào)進(jìn)行濾波和采樣率同步化處理。另外,文獻(xiàn)[27]進(jìn)行了肌電爆發(fā)時(shí)刻與關(guān)節(jié)力矩產(chǎn)生時(shí)刻之間的時(shí)間差檢測(cè),結(jié)果表明,關(guān)節(jié)力矩產(chǎn)生滯后表面肌電信號(hào)50～80 ms;文獻(xiàn)[11]在利用表面肌電信號(hào)進(jìn)行關(guān)節(jié)力矩預(yù)測(cè)時(shí)將滯后時(shí)間定為40 ms。為此,本文綜合各方面因素,將運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)數(shù)據(jù)都提前50 ms。利用同步化數(shù)據(jù),按照第2、3節(jié)的方法分別對(duì)4名對(duì)象進(jìn)行正向生物力學(xué)模型參數(shù)的標(biāo)定和預(yù)測(cè)。

圖5 試驗(yàn)設(shè)備與場(chǎng)景圖

圖6 對(duì)象A在3個(gè)步態(tài)周期中肌肉激活程度變化曲線

圖7 對(duì)象A在3個(gè)步態(tài)周期中肌肉纖維長(zhǎng)度變化曲線

以對(duì)象A為例,從3種步速下各選取一個(gè)步態(tài)周期數(shù)據(jù)作為優(yōu)化樣本,對(duì)模型進(jìn)行標(biāo)定。首先計(jì)算3個(gè)步態(tài)周期中的肌肉激活程度、肌肉運(yùn)動(dòng)參數(shù)以及關(guān)節(jié)力矩。圖6是對(duì)象A在3個(gè)步態(tài)周期下10塊肌肉的激活程度變化情況。從圖中可以看出,在步態(tài)過程中,每塊肌肉的激活程度變化非常復(fù)雜,但是在行走這一特定任務(wù)中,肌肉激活程度的變化趨勢(shì)基本固定,主要在大小上存在差異。圖7是對(duì)象A在3個(gè)步態(tài)周期下的肌肉纖維長(zhǎng)度變化情況,可以看出,步態(tài)周期中肌纖維長(zhǎng)度變化緩和,周期性很好,其中縫匠肌(SR)起始于骨盆髂前上棘,止于脛骨上端前緣,是人體最長(zhǎng)的肌肉。其次,根據(jù)肌肉激活程度、肌肉運(yùn)動(dòng)學(xué)參數(shù)以及膝關(guān)節(jié)屈伸力矩值,對(duì)模型參數(shù)進(jìn)行優(yōu)化,初始種群規(guī)模為60,最大代數(shù)為400,優(yōu)化結(jié)果誤差如圖8所示,優(yōu)化最大誤差達(dá)到6.25 N·m。此外,表1給出了對(duì)象A經(jīng)過模型標(biāo)定后的肌肉部分參數(shù),包括肌肉最大力Famax、最優(yōu)長(zhǎng)度lopti和高斯函數(shù)中的標(biāo)準(zhǔn)差σ,同時(shí)還給出了文獻(xiàn)[26]基于Hill肌肉力學(xué)模型給出的肌肉參數(shù)。在Delp模型中,通過測(cè)量21名尸體相關(guān)肌肉的橫截面積取均值并乘以單位面積最大張力(61 N/cm2),得到肌肉最大力,因此肌肉之間最大力之比等于肌肉橫截面積之比,而肌肉最優(yōu)長(zhǎng)度通過尸體肌肉直接測(cè)量得到。通過比較,本文優(yōu)化得到的不同肌肉之間最大力之比與Delp模型中不同肌肉之間最大力之比關(guān)系一致,即遵循實(shí)際人體不同肌肉之間的發(fā)達(dá)程度關(guān)系,而肌肉的最優(yōu)長(zhǎng)度相差在1～2 cm之間,最長(zhǎng)的肌肉縫匠肌(SR)相差為8 cm,這可以解釋為是由于個(gè)體之間肌肉發(fā)達(dá)程度、身高差異引起的。

圖8 模型參數(shù)優(yōu)化結(jié)果誤差

肌肉本文模型Famax/Nlopti/mDelp模型Famax/Nlopti/mσSM12390.09610300.0800.126BFL8370.1037200.0890.112BFS3870.0934000.1100.657VM11660.08912950.0970.323VL14360.07518710.0840.263RF7700.0887800.0760.125MG12500.04211150.0510.155LG5610.0504900.0590.245TFL1630.1301550.1600.656SR1060.3981050.4790.570

模型標(biāo)定后,對(duì)所有對(duì)象步行過程中的膝關(guān)節(jié)屈伸方向動(dòng)態(tài)力矩進(jìn)行預(yù)測(cè),圖9給出了對(duì)象A某次步態(tài)下的預(yù)測(cè)結(jié)果,圖中還給出了逆動(dòng)力學(xué)計(jì)算的動(dòng)態(tài)力矩曲線。從圖中可以看出,建立的生物力學(xué)模型對(duì)步態(tài)行走過程中的膝關(guān)節(jié)動(dòng)態(tài)力矩有較好的預(yù)測(cè)能力,能夠快速地捕捉實(shí)際關(guān)節(jié)力矩的變化趨勢(shì),而且預(yù)測(cè)曲線波動(dòng)較小。

圖9 膝關(guān)節(jié)動(dòng)態(tài)關(guān)節(jié)力矩預(yù)測(cè)結(jié)果

另外,為定量描述模型的預(yù)測(cè)效果,采用最大絕對(duì)誤差、平均殘差以及預(yù)測(cè)值和實(shí)際值的線性相關(guān)程度R這3個(gè)參數(shù)來描述預(yù)測(cè)精度。對(duì)4名對(duì)象、3種不同步速下膝關(guān)節(jié)力矩預(yù)測(cè)結(jié)果進(jìn)行定量統(tǒng)計(jì),每名對(duì)象的描述參數(shù)都是對(duì)自身預(yù)測(cè)結(jié)果統(tǒng)計(jì)的平均值,如表2所示。從表中可以看出,正向生物力學(xué)模型預(yù)測(cè)最大絕對(duì)誤差平均值為(11.0±1.32) N·m,平均殘差為(4.43±0.698) N·m,平均線性相關(guān)系數(shù)為0.927±0.042。

表2 所有對(duì)象預(yù)測(cè)結(jié)果統(tǒng)計(jì)

4 結(jié) 論

康復(fù)訓(xùn)練機(jī)器人能否有效地感知患者自身運(yùn)動(dòng)能力是其提供合理輔助的重要前提?；谏镫姷娜藱C(jī)交互接口技術(shù)是康復(fù)機(jī)器人領(lǐng)域的重要發(fā)展方向。本文在對(duì)人體下肢運(yùn)動(dòng)系統(tǒng),包括神經(jīng)、肌肉、骨骼、關(guān)節(jié),各個(gè)組織的解剖結(jié)構(gòu)和運(yùn)動(dòng)機(jī)理進(jìn)行深入研究的基礎(chǔ)上,建立了完整過程的正向生物力學(xué)模型,實(shí)現(xiàn)了基于表面肌電信號(hào)的下肢關(guān)節(jié)力矩動(dòng)態(tài)預(yù)測(cè)。本文模型根據(jù)表面肌電信號(hào)的產(chǎn)生機(jī)理,從頻率和幅值角度描述了肌肉激活程度與肌纖維激活及募集程度的關(guān)系,不同于過去只通過幅值信息描述肌肉激活程度,本文建立的肌肉生物電激活模型更全面。綜合幅值和頻率與單一的幅值相比,可以更加直接影響肌肉激活程度,進(jìn)一步影響肌肉力大小及優(yōu)化后參數(shù)的準(zhǔn)確性,然而由于人體內(nèi)肌力大小和肌肉參數(shù)目前無(wú)法直接測(cè)量,因此本文對(duì)于參數(shù)的影響尚不能給出具體評(píng)價(jià),有待進(jìn)一步研究。與Hill肌肉模型不同,本文根據(jù)肌絲滑移理論,建立能夠描述生理結(jié)構(gòu)和微觀收縮機(jī)理的肌肉力學(xué)模型,因此更加完整和切合實(shí)際收縮過程。膝關(guān)節(jié)屈伸動(dòng)態(tài)力矩的試驗(yàn)預(yù)測(cè)結(jié)果證明了本文模型的準(zhǔn)確性和可靠性,為康復(fù)機(jī)器人運(yùn)動(dòng)感知能力的提高和人、機(jī)力柔順交互研究提供了接口。同時(shí),該生物力學(xué)模型的研究對(duì)不同運(yùn)動(dòng)的生理機(jī)理、探索運(yùn)動(dòng)康復(fù)規(guī)律、提高運(yùn)動(dòng)員成績(jī)和防止運(yùn)動(dòng)損傷等領(lǐng)域存在應(yīng)用價(jià)值。

[2] LUM P, REINKENSMEYER D, MAHONEY R, et al. Robotic devices for movement therapy after stroke: current status and challenges to clinical acceptance [J]. Topics in Stroke Rehabilitation, 2002, 8(4): 40-53.

[3] EKKELENKAMP R, VENEMAN J, VANDER KOOIJ H. LOPES: selective control of gait functions during the gait rehabilitation of CVA patients [C]∥9th International Conference on Rehabilitation Robotics. Piscataway, NJ, USA: IEEE, 2005: 361-364.

[4] 文忠, 錢晉武, 沈林勇, 等. 基于阻抗控制的步行康復(fù)訓(xùn)練機(jī)器人的軌跡自適應(yīng) [J]. 機(jī)器人, 2011, 33(1): 142-149. WEN Zhong, QIAN Jinwu, SHEN Linyong, et al. Trajectory adaptation for impedance control based walking rehabilitation training robot [J]. Robot, 2011, 33(1): 142-149.

[5] HAYASHI T, KAWAMOTO H, SANKAI Y. Control method of robot suit HAL working as operator’s muscle using biological and dynamical information [C]∥IEEE/RSJ International Conference on Intelligent Robots and Systems. Piscataway, NJ, USA: IEEE, 2005: 3063-3068.

[6] FLEISCHER C, HOMMEL G. A human-exoskeleton interface utilizing electromyography [J]. IEEE Transactions on Robotics, 2008, 24(4): 872-882.

[7] KREBS H I, PALAZZOLO J J, DIPIETRO L, et al. Rehabilitation robotics: performance based progressive robot-assisted therapy [J]. Autonomous Robots, 2003, 15(1): 7-20.

[8] PERRY J C, ROSEN J, BURNS S. Upper-limb powered exoskeleton design [J]. IEEE/ASME Transactions on Mechatronics, 2007, 12(4): 408-417.

[9] KREBS H, VOLPE B, WILLIAMS D, et al. Robot-aided neurorehabilitation: a robot for wrist rehabilitation [J]. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2007, 15(3): 327-335.

[10]MORITA S, KONDO T, ITO K. Estimation of forearm movement from EMG signal and application to prosthetic hand control [C]∥International Conference on Robotics and Automation. Piscataway, NJ, USA: IEEE, 2001: 3692-3697.

[11]LLOYD D G, BESIER T F. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo [J]. Journal of Biomechanics, 2003, 36(6): 765-776.

[12]BUCHANAN T S, LLOYD D G, MANAL K, et al. Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command [J]. Journal of Applied Biomechanics, 2004, 20(4): 367-395.

[13]HILL A V. The heat of shortening and the dynamic constants of muscle [J]. Proceedings of the Royal Society of London: Series B Biological Sciences, 1938, 126(843): 136-195.

[14]GUO Zhao, FAN Yuanjie, ZHANG Jianjun, et al. A new 4M model-based human-machine interface for lower extremity exoskeleton robot [C]∥The 5th International Conference on Intelligent Robotics and Applications. Montreal, Canada: ICIRA, 2012: 545-554.

[15]YIN Y H, CHEN X. Bioelectrochemical control mechanism with variable-frequency regulation for skeletal muscle contraction-biomechanics of skeletal muscle based on the working mechanism of myosin motors (II) [J]. Science China: Technological Sciences, 2012, 55: 2115-2125.

[16]WANG W, DE STEFANO A, ALLEN R. A simulation model of the surface EMG signal for analysis of muscle activity during the gait cycle [J]. Computers in Biology and Medicine, 2006, 36(6): 601-618.

[17]HUXLEY A F, SIMMONS R M. Proposed mechanism of force generation in striated muscle [J]. Nature, 1971, 233(5321): 533-538.

[18]JULICHER F, AJDARI A, PROST J. Modeling molecular motors [J]. Reviews of Modern Physics, 1997, 69(4): 1269-1281.

[19]LECARPENTIER Y, BLANC F X, QUILLARD J, et al. Statistical mechanics of myosin molecular motors in skeletal muscles [J]. Journal of Theoretical Biology, 2005, 235(3): 381-392.

[20]YIN Yuehong, GUO Zhao. Collective mechanism of molecular motors and a dynamic mechanical model for sarcomere [J]. Science China: Technological Sciences, 2011, 54(8): 2130-2137.

[21]GOLLAPUDI S K, LIN D C. Experimental determination of sarcomere force-length relationship in type-I human skeletal muscle fibers [J]. Journal of Biomechanics, 2009, 42(13): 2011-2016.

[22]單大卯. 人體下肢肌肉功能模型及其應(yīng)用的研究 [D]. 上海: 上海體育學(xué)院, 2003.

[23]CAMOMILLA V, CEREATTI A, VANNOZZI G, et al. An optimized protocol for hip joint center determination using the functional method [J]. Journal of Biomechanics, 2006, 39(6): 1096-1106.

[24]HARRINGTON M E, ZAVATSKY A B, LAWSON S E, et al. Prediction of the hip joint centre in adults, children, and patients with cerebral palsy based on magnetic resonance imaging [J]. Journal of Biomechanics, 2007, 40(3): 595-602.

[25]VAUGHAN C L, DAVIS B L, O’CONNOR J C. Dynamics of human gait [M]. Champaign, Illinois, USA: Human Kinetics Publishers, 1992: 83-106.

[26]ARNOLD E M, WARD S R, LIEBER R L, et al. A model of the lower limb for analysis of human movement [J]. Annals of Biomedical Engineering, 2010, 38(2): 269-279.

[27]HOWATSON G, GLAISTER M, BROUNER J, et al. The reliability of electromechanical delay and torque during isometric and concentric isokinetic contractions [J]. Journal of Electromyography and Kinesiology, 2009, 19(5): 975-979.

(編輯 杜秀杰)

Prediction Model for Dynamic Joint Torque of Lower Limb with Surface EMG

CHEN Jiangcheng,ZHANG Xiaodong,LI Rui,SHI Qiangyong,WANG He

(Key Laboratory of Education Ministry for Modern Design & Rotor-Bearing System, Xi’an Jiaotong University, Xi’an 710049, China)

To achieve the dynamic joint torque decoding from surface electromyography (EMG), the forward biomechanical model of lower limb motor system, which relates the surface EMG and joint torque, is established. The dynamic surface EMG to skeletal muscle activation model is constructed from the perspective of amplitude and frequency. Then the muscle contraction model reflecting physiological structure and micromechanical properties is constructed according to the sliding-filament theory. The force direction and displacement vector of active muscle are determined and the transformation from muscle force to joint moment is realized. The dynamic calibration for the forward biomechanical model using the exact joint torque value obtained with Newton-Eular method is finally put forward. Following the calibration, the flexion/extension (FE) knee joint torque of four objects under different speed walking is predicted. The results show that the forward biomechanical model can capture the general shape and timing of the joint torque, the maximum absolute error is (11.0±1.32) N·m, the mean residual error is (4.43±0.698) N·m, and the linear relationship between predicted and exact knee FE torque reaches 0.927±0.042. This prediction model provides an interface for the study of force interaction pattern in the process of human-machine cooperation in training.

surface EMG; joint torque prediction; muscle model; forward biomechanics

2015-03-08。

陳江城(1987—),男,博士生;張小棟(通信作者),男,教授,博士生導(dǎo)師。

國(guó)家高技術(shù)研究發(fā)展計(jì)劃資助項(xiàng)目(2015AA042301);國(guó)家自然科學(xué)基金重大研究計(jì)劃資助項(xiàng)目(91420301)。

時(shí)間:2015-09-13

10.7652/xjtuxb201512005

TP242

A

0253-987X(2015)12-0026-08

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20150913.1822.002.html