基于Matlab的人體膝關節動力學仿真平臺設計

唐 剛,程榮山,胡 雄,張博一,吳 鋼

(上海海事大學 物流工程學院,上海 201306)

基于Matlab的人體膝關節動力學仿真平臺設計

唐 剛,程榮山,胡 雄,張博一,吳 鋼

(上海海事大學 物流工程學院,上海 201306)

在對小腿運動學和膝關節動力學分析的基礎上,開發了一款基于Matlab的各種步態下人體膝關節動力學仿真平臺. 基于反向動力學原理,編寫了Matlab計算程序,設計了仿真平臺的圖形用戶界面(GUI),通過讀取運動捕捉系統采集的下肢運動信息及三維測力臺的支反力,實現了各種步態下膝關節關節力、關節力矩的快速計算. 以正常步態為例,通過對一名健康人體右側膝關節分析,并與Visual 3D軟件結果進行對比,驗證了該平臺的可靠性,從而為人工膝關節置換和膝關節康復運動評定提供了一種新的分析工具.

膝關節動力學;Matlab;批量讀取;快速計算;步態

隨著計算機技術的迭代更新,應用計算機對人體下肢膝關節進行建模和動態仿真變得越來越重要[1-7],如下肢膝關節置換前的力學分析、患者下肢膝關節的康復運動評定、膝關節風濕病所引起的非正常行走的分析等. 膝關節是下肢運動的最主要關節之一,對各種步態運動的研究起著非常顯著的作用. 同時,膝關節又承載著人體的絕大部分重量,由于它的受力活動頻繁,膝關節更易受到傷害. 因此,研究膝關節的動力學參數對人體生物力學、醫療等領域也具有重要意義[8-10]. 目前,國外成熟的運動分析軟件主要有Visual 3D步態分析軟件等. 雖然這些軟件精度高,但是需要人工手動建模,建模速度慢,且需要專業的建模知識,耗費的時間較長,大大減慢了數據處理的速度,因此,設計一個可以快速自動建模,且精度在許可范圍內相對較高的仿真平臺是本文要解決的問題.

本文基于反向動力學原理,利用Matlab編寫了人體膝關節的動力學程序,設計了仿真平臺的圖形用戶界面(GUI),通過讀取運動捕捉系統采集的下肢運動信息及三維測力臺的支反力數據,實現了膝關節動力學參數的輸出.

1 試驗數據采集

1.1 試驗對象與試驗設備

試驗對象為一名健康男性,身高1.80 m,體重67.0 kg,無關節或肌肉損傷病史.

試驗采用OPTOTRAK CERTUS運動捕捉系統、兩臺AMTI測力平臺同步測量. 運動捕捉系統主要由位置傳感器、控制中心和Marker點組成,采樣頻率為128 Hz. AMTI測力平臺主要由兩塊測力板組成,采樣頻率為512 Hz. 試驗前首先將Marker點粘貼在加熱彎曲的有機玻璃板上,然后將帶有Marker點的有機玻璃板貼附在試驗對象大腿、小腿和足部的側前部,并用彈性繃帶固定. 試驗過程中按要求完成5次符合試驗要求的步態(如圖1所示).

圖1 步態測量Fig.1 Gait measurement

1.2 試驗測量

由于運動捕捉系統捕捉的Marker點要用全局坐標系表示,需要建立統一的全局坐標系. 本文則以測力板的某直角交點為圓心,兩直角邊分別為x軸方向和y軸方向,豎直向上為z軸方向. 全局坐標系確定后,對試驗對象進行正常步態測量. 試驗要求: 試驗前對試驗對象進行適應性訓練直至按節拍行走,試驗中則要求試驗對象的行走節奏與節拍器相同及步態的步長保持穩定. 為保證測力板上的支反力準確,試驗過程中僅單腳接觸測力板. 試驗測量中,運動捕捉系統每次采集640幀數據,測力板每次采集2 560幀數據,數據采集結束后,運動捕捉系統采集的下肢Marker點的空間坐標原始數據和測力板采集的支反力保存為相應名稱的xls格式文件.

2 基于Matlab GUI人體膝關節動力學仿真平臺的實現

2.1 平臺運算的主要參數

平臺運算的主要參數包括小腿運動學參數(小腿轉動角、小腿角速度、小腿角加速度)和膝關節動力學參數(膝關節力、膝關節力矩). 各參數的數學符號和在本試驗中的物理意義如表1所示.

表1 小腿運動學參數和膝關節動力學參數Table 1 Kinematics parameters for shank and dynamics parameters for knee joint

2.2 平臺的算法設計

本平臺的算法設計基于反向動力學原理,引入空間向量坐標編寫Matlab程序求解小腿的運動學參數,利用小腿運動學參數進而求得膝關節動力學參數[11]. 仿真平臺算法設計主要圍繞膝關節動力學算法展開.

為了獲得膝關節動力學參數首先對小腿進行運動學研究,主要包括對小腿角速度和角加速度的分析. 小腿段對應踝關節和膝關節,小腿段轉動角即為局部坐標系下小腿圍繞x軸、y軸、z軸的旋轉角度. 每一幀對應一個小腿轉動角度值,相鄰兩幀的角度變化量dθ與兩幀時間間隔dt的商即為該時間內的小腿平均角速度ω,角加速度αi的計算與角速度基本相同. 小腿角速度ωi、角加速度αi涉及的算法以及對應的Matlab的程序(圖2)如下所述.

圖2 小腿角速度和角加速度的Matlab程序圖Fig.2 Matlab program diagram of shank angular velocity and angular acceleration

踝關節中心的向量表示為

R1=(RLA+RMA)/2

(1)

其中:RLA為踝外側Marker點坐標;RMA為踝內側Marker點坐標.

膝關節中心的向量表示為

R2=(RLK+RMK)/2

(2)

其中:RLK為膝外側Marker點坐標;RMK為膝內側Marker點坐標.

在小腿的局部坐標系中,za軸由踝關節中心指向膝關節中心:

za=R1-R2

(3)

ya軸垂直于由踝關節中心、RLK、RMK三點確定的平面:

ya=za×(RMK-RLK)

(4)

根據右手螺旋定則進一步確定xa軸:

xa=ya×za

(5)

由于每一幀的坐標已知,小腿繞x軸旋轉α角,旋轉后x軸方向的坐標值不變,y軸和z軸方向的坐標值變化相當于在yOz平面內作正α角旋轉:

(6)

小腿繞y軸旋轉β角,旋轉后y軸方向的坐標值不變,z軸和x軸方向的坐標值變化相當于在zOx平面內作正β角旋轉:

(7)

小腿繞z軸旋轉γ角,旋轉后z軸方向的坐標值不變,x軸和y軸方向的坐標值變化相當于在xOy平面內作正γ角旋轉:

(8)

綜上,小腿的旋轉角為

(9)

小腿角速度ωi、角加速度αi的公式為

(10)

利用上文獲得的小腿運動學分析參數對膝關節進行動力學研究,主要包括對膝關節關節力和關節力矩的分析. 通常在人體步態運動研究中會采用多剛體動力學方法將人體簡化為待研究的關節鏈,即將關節鏈的各段近似為剛體,從而將人體簡化成為具有有限個自由度的多剛體模型[12-15],因此,可以通過多剛體模型計算出與人體實際運動等效的膝關節的關節力和關節力矩. 基于牛頓動力學和歐拉動力學方程,在全局坐標系下膝關節的關節力(F)和關節力矩(M)涉及的算法以及對應的Matlab的程序(圖3)如下所述.

(11)

其中:mF和mS分別為足部、小腿的質量;aF和aS分別為足部、小腿的全局坐標系下的質心加速度;FD為測力臺測得的足底反力;MG為足部和小腿的重力產生的力矩;M慣性力為小腿和足部的慣性力產生的慣性力矩;MJ為小腿轉動角加速度和小腿質心相對膝關節中心坐標產生的力矩;MD為足底反力產生的力矩.

圖3 膝關節力和關節力矩的Matlab程序圖Fig.3 Matlab program diagram of knee joint force and knee joint torque

圖4 膝關節動力學仿真平臺的GUIFig.4 GUI of knee dynamics simulation platform

2.3 平臺的GUI設計

人體膝關節動力學仿真平臺界面設計以Matlab為編程語言,平臺界面主要由數據輸入(Import Data)、仿真運行(Simulation)、小腿段運動學參數輸出(Shank)、膝關節動力學參數輸出(Knee joint)、工具箱(Pipeline)、下肢棍棒行走演示(Begin、Pause、Continue)共6個功能模塊構成,如圖4所示.

數據輸入模塊包括手動輸入試驗對象的身高和體重,試驗采樣頻率以及讀取的xls數據單元格位置等. 仿真運行模塊主要用于對平臺算法的Matlab程序的運行. 小腿運動學參數輸出模塊包括本平臺算法輸出的小腿角速度、角加速度曲線圖，以及兩者與Visual 3D輸出的小腿角速度、角加速度對比曲線圖. 膝關節動力學參數輸出模塊包括本平臺算法輸出的膝關節關節力、關節力矩曲線圖，以及兩者與Visual 3D輸出的膝關節關節力、關節力矩對比曲線圖. 工具箱模塊包括優化工具箱(Optimtool)和曲線擬合工具箱(Sftool),分別用于對平臺算法的優化以及平臺輸出曲線的擬合. 下肢棍棒行走演示模塊主要通過Matlab對下肢在一個完整步態周期內的行走過程的再現.

3 試驗結果分析

3.1 平臺結果輸出及分析

將試驗對象的身高與體重、試驗采樣頻率、保存的試驗采集的xls數據單元格位置等參數，通過數據輸入模塊輸入到平臺內(隨機選擇多次采集中的5組比較穩定的、完整的右側下肢步態數據作為本次試驗數據),之后運行仿真模塊,結束后即可獲得試驗對象右小腿的角速度、角加速度和右膝關節的關節力、關節力矩的曲線圖和數值矩陣. 曲線圖則可以直觀地反映膝關節正常步態下的關節力和關節力矩的變化,而數值矩陣則可以直接保存使用或用作下一步計算的原始數據.

圖5和6分別為本平臺和Visual 3D軟件在一個運動周期內沿x軸、y軸、z軸3個方向輸出的膝關節的關節力、關節力矩的對比曲線圖(實線表示本平臺的輸出曲線,點劃線表示Visual 3D軟件的輸出曲線， 5次試驗輸出曲線按照線寬由細到粗加以區分).

(a) x軸方向

(b) y軸方向

(c) z軸方向圖5 膝關節關節力對比曲線圖Fig.5 Contrast curves of knee joint force

(a) x軸方向

(b) y軸方向

(c) z軸方向圖6 膝關節力矩對比曲線圖Fig.6 Contrast curves of knee joint torque

由圖5可以看出,本平臺和Visual 3D 5次試驗的關節力曲線基本吻合.x軸方向上,運動周期在10%左右時,出現了明顯的波谷(圖5(a));y軸方向上,一個運動周期內出現了明顯的波峰和波谷(圖5(b));z軸方向上,一個運動周期內出現了兩個波谷(圖5(c)). 由圖6可以看出,本平臺和Visual 3D 5次試驗的關節力矩曲線基本吻合.x軸方向上,運動周期在5%左右時,出現了小的波谷,運動周期在20%左右時,出現了明顯的波峰(圖6(a));y軸方向上,一個運動周期內出現了一次波峰和兩次波谷(圖6(b));z軸方向上,一個運動周期內出現了兩次波峰和一次波谷(圖6(c)).

綜上所述,膝關節在一個運動周期中,運動周期從 0到 40%為試驗對象右腿接觸到測力板的過程,因右腿屈伸過程中急劇減速使得膝關節沿3個方向所受關節力和關節力矩出現了明顯的波峰和波谷;運動周期從40%到60%為試驗對象右腿在測力板上保持平衡的過程,膝關節沿3個方向所受關節力、關節力矩無明顯波動;運動周期從60%到100%為試驗對象逐漸離開測力板過程,膝關節沿3個方向所受關節力和關節力矩再一次出現了明顯的波峰和波谷并直至趨于0.

3.2 平臺可靠性驗證

目前處理人體運動測量數據較為成熟的C-Motion公司Visual 3D步態/體態分析軟件,可以利用運動捕捉系統采集人體運動信息作為原始數據用于人體運動學、動力學參數的計算. 利用SPSS 20.0軟件對平臺輸出的關節力、關節力矩與Visual 3D軟件的處理結果進行相關性分析，結果如表2和3所示. 二者的相關程度越高,說明計算結果越趨于一致. 二者的相關程度由相關系數決定,當相關系數取值范圍為0.7～1.0,則二者呈高度相關;當相關系數取值范圍為0.40～0.69,則二者呈中度相關;當相關系數取值為0～0.39,則二者呈低度相關. 由表2可知,二者輸出的關節力在x軸、y軸、z軸3個方向均呈現出高度相關. 由表3可知,二者輸出的關節力矩在x軸和y軸方向呈高度相關,而在z軸方向呈現中度相關,這主要因為二者模型具有一定的差異性. 綜上所述,就5次試驗的相關性和總體趨勢而言,二者的結果基本吻合,從而驗證了本仿真平臺的可靠性.

表2 平臺與Visual 3D的關節力在3個坐標軸方向的相關性分析
Table 2 Correlation analysis of joint force on the platform and Visual 3D in three axes direction

變量第1次第2次第3次第4次第5次xyzxyzxyzxyzxyzPearson相關性0.720**0.975**0.999**0.946**0.997**1.00**0.914**0.994**1.00**0.973**0.974**1.00**0.898**0.998**1.00**顯著性(雙側)000000000000000N100100100100100100100100100100100100100100100

注: **. 在 0.01 水平(雙側)上顯著相關.

表3 平臺與Visual 3D的關節力矩在3個坐標軸方向的相關性分析
Table 3 Correlation analysis of joint torque on the platform and Visual 3D in three axes direction

變量第1次第2次第3次第4次第5次xyzxyzxyzxyzxyzPearson相關性0.979**0.689**0.030**0.988**0.982**0.496**0.992**0.964**0.437**0.993**0.984**0.647**0.992**0.971**0.502**顯著性(雙側)000.767000000000000N100100100100100100100100100100100100100100100

注: **. 在 0.01 水平(雙側)上顯著相關.

3.3 平臺結果的實際應用

本平臺設計主要應用于醫療領域,如人工膝關節的假體設計、膝關節的康復訓練等,其中，針對人工膝關節的假體設計,如通過本平臺輸出膝關節假體的關節力和關節力矩,則可以為膝關節假體的設計提供依據. 而對于那些已經安裝膝關節假體的患者,可以通過采集患者的運動參數輸入本平臺對關節力和關節力矩的計算,判斷其膝關節是否過載. 如果過載較輕,可以通過康復訓練對其調節;如果過載較為嚴重,可以考慮為其更換假體.

4 結 語

本文闡述了基于Matlab所設計的人體膝關節動力學仿真平臺的開發過程,并通過對行走步態下膝關節動力學參數即關節力和關節力矩的分析,驗證了本平臺的可靠性,為人工膝關節置換和膝關節康復運動評定提供了一種新的分析工具. 與目前常用且以C++或Java開發的運動分析軟件相比,本平臺利用Matlab編程,程序代碼比較簡潔,沒有復雜的建模環節,通過快速自動建模有效地解決了目前常用運動分析軟件的人工手動建模、建模速度慢、較高的專業知識要求等帶來的不便,同時對不同人體運動數據批量處理和分析的效率提高起著顯而易見的作用. 但本平臺與目前常用運動分析軟件相比,還存在許多不足,如精度相對較低、缺少對人體其他體段運動學和關節動力學參數的分析等,在以后的工作中還需對本平臺的功能不斷完善以及計算結果做更深一步的分析.

[1] 孫瑜. 人下肢運動鏈的運動建模仿真分析[J]. 計算機仿真,2015,32(10): 365-368,388.

[2] 于佳彬,郝衛亞,周興龍. 縱跳落地動作地面反作用力計算機仿真方法的研究[J]. 天津體育學院學報,2013,28(6): 497-501.

[3] 劉書朋,司文,嚴壯志,等. 基于AnyBodyTM 技術的人體運動建模方法[J]. 生物醫學工程學進展,2010,31(3): 131-134.

[4] 郝衛亞. 人體運動的生物力學建模與計算機仿真進展[J]. 醫用生物力學,2011,26(2): 97-104.

[5] 唐剛,魏高峰,聶文忠,等. 人體下肢關節坐標系的一種簡單定義方法[J]. 北京生物醫學工程,2009,28(6): 606-609.

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

[7] HOLZBAUR K R,MURRAY W M,DELP S L. A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control [J]. Ann Biomed Eng,2005,33 (6): 829-840.

[8] ALKNER B A ,TESCH P A. Knee extensor and plantar flexor muscle size and function following 90 days of bed rest with or without resistance exercise [J]. Eur J Appl Physiol,2004,93(3): 294-305.

[9] 唐剛,白雪嶺,王洪生,等. 基于關節坐標系的人體運動學仿真[J]. 計算機仿真,2011,28(8): 94-97,241.

[10] FUKASHIRO S,HAY D C,NAGANO A. Biomechanical behavior of muscle-tendon complex during dynamic human movements[J]. Journal of Applied Biomechanics,2006,22 (2): 131-147.

[11] TANG G,WANG C T. A muscle-path-plane method for representing muscle contraction during joint movement [J]. Computer Methods in Biomechanics and Biomedical Engineering,2010,14(1): 59-69.

[12] HATZE H. A comprehensive model for human monition simulation and its application to the take-off phase of the long jump [J]. J Biomechanics,1981,14(3): 135-142.

[13] PANDY M G,ZAJAC F E. Optimal muscular coordination strategies for jumping [J]. J Biomechanics,1991,24(1): 1-10.

[14] MILLER D I. A computer simulation model of the airborne phase of diving [D]. Pennsylvania: Pennsylvania State University,1970.

[15] PASSERELLO C E,HUSTON R L. Human attitude control [J]. J Biomechanics,1971,4(2): 95-102.

Simulation Platform Design for Human Knee Joint Dynamics Based on Matlab

TANGGang,CHENGRong-shan,HUXiong,ZHANGBo-yi,WUGang

(School of Logistics Engineering,Shanghai Maritime University,Shanghai 201306,China)

On the basis of the analysis of shank kinematics and knee joint dynamics,a Matlab-based under various kinds of gait of simulation platform for human knee joint dynamics was designed. Basing on inverse dynamics principle to write the Matlab program,graph user interface (GUI) of simulation platform was designed,through reading lower limb movement information from motion capture system and reaction force from 3D load stage,knee joint force and joint torque for various gaits were calculated rapidly. With a normal gait as an example,the right knee joint of a health subject was analyzed,the reliability of the platform was verified by comparing the results with those generated by Visual 3D software,and the platform also provided new analysis tools for the artificial knee joint replacement and knee rehabilitation exercise evaluation.

knee joint dynamics;Matlab;batch read;fast calculation;gait

1671-0444 (2016)05-0725-07

2016-01-19

國家高技術研究發展計劃(863)資助項目(2013A2041106);國家自然科學基金資助項目(31300783);中國博士后科學基金資助項目(2014M561458);教育部博士點基金聯合資助項目(20123121120004);上海海事大學科研基金資助項目(20130474);上海高校一流學科-管理科學與工程資助項目

唐 剛 (1982—),男,重慶人,副教授,博士,研究方向為人機工程. E-mail: gangtang@shmtu.edu.cn

TP 391

A