針對人體站立平衡姿態保持過程中, 下肢主要肌肉的肌力變化分布的最優求解問題。本研究將人體下肢運動肌肉骨骼簡化為具有3關節和9塊肌肉的平面物理模型, 并在此基礎上建立了用于冗余肌力優化求解的數學模型。分別利用粒子群優化(POS)單目標和多目標算法進行最優化求解。數值計算的結果表明多目標優化可以更合理地得到9組肌力的分布及變化規律。最后, 通過對仿真結果的分析, 定性地分析了被動運動下人體恢復站立平衡過程中各肌肉群的運動協調規律。
引用本文： 王洪瑞, 鄭輝, 劉琨. 站立平衡調節的肌力優化求解與分析. 生物醫學工程學雜志, 2015, 32(1): 59-66. doi: 10.7507/1001-5515.20150011 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
人體站立平衡機能一直是生命科學研究的熱點和難點，其平衡過程是在小腦的控制支配下，通過神經系統傳遞生物電信號刺激肌肉收縮，牽拉肌肉所附著的骨骼繞關節旋轉，改變骨骼間的相對位置，從而協調人體恢復站立平衡。生活中由于意外傷害或者中風等疾病導致肌肉組織功能障礙，引起的平衡能力失調，導致了大量的傷害事故。因此，預測恢復站立平衡過程中肌力的分布變化能夠更好地理解其調節機制，為醫療康復起到實際指導作用。采用植入性的儀器測量人體肌肉作用力，不僅技術復雜，成本高，而且會對人體造成創傷，不適合進行正常人體肌力的研究[1]。因此，人體肌肉骨骼模型為間接求解肌力變化分布提供了新思路。由于關節力矩是在多塊肌肉的共同作用下產生的，故一般多采用約束優化方法求解冗余肌力的分布情況。針對不同問題的肌力變化規律，諸多學者采用了不同的優化算法，Neptune[2]、Raikova等[3]分別采用模擬退火(simulated annealing, SA)算法和拉格朗日乘子法對腳踏板曲柄運動過程中下肢肌力分布進行了優化求解，Seth等[4]利用神經骨肌跟蹤(neuromusculoskeletal tracking, NMT)法分析了縱跳過程中下肢肌力的變化規律。然而，上述采用的算法均只是單個目標函數的肌力優化求解。
迄今為止，利用智能優化算法求解多目標問題被廣泛應用于軌跡規劃、系統工程和控制論等領域。Santana-Quintero等[5]針對帶有約束的多目標問題采用粗集理論與差分進化結合的算法予以求解，Bensghaier等[6]采用遺傳算法(genetic algorithm，GA)優化求解了拇指與食指捏起物品過程中肌力的分布與變化。多目標粒子群優化(multi-objective particle swarm optimization, MOPSO)算法同樣備受諸多學者的青睞，如Tan等[7]利用其求解工程約束問題，Carvalho等[8]應用該算法測試了多目標問題中支配解的收斂性和多樣性。Yang[9]利用改進后的粒子群算法解決了多目標離散優化問題。但是，目前該算法還沒有被用于預測人體恢復站立平衡過程中下肢肌力的分布與變化規律。因此，本文在建立的3-DOF平面運動仿真簡化模型的基礎上采用混合罰函數化簡約束，然后利用粒子群優化(particle swarm optimization, PSO)算法分別求解恢復站立平衡過程中單目標函數和多目標函數的下肢9組肌力最優值，進而分析各組肌力在該過程中的變化規律和發揮的作用。
1 模型的分析與建立
1.1 下肢骨骼肌肉簡化模型
本文旨在預測被動運動下人體站立姿態由失衡至恢復平衡過程中肌肉力動態變化的最優值。研究中所施加的外力具體指水平運動平臺對足底提供后向水平加速度[10]。根據此外力的作用方式和人體下肢生理結構，可近似認為腳掌始終保持與平臺平穩接觸且兩腿的運動關于矢狀面左右對稱。建立二維平面下肢骨架3自由度、4剛體(包含軀干、大腿、小腿和腳)和9塊肌肉的運動簡化物理模型如圖 1所示。參考人體解剖學數據和文獻，選取與下肢運動相關的9組主要肌肉，包括臀大肌(gluteus maximus，GM)、髂肌(iliopsoas，IL)、股直肌(rectus femoris，RF)、股側肌群(vastus group，VA)、股二頭肌長頭(biceps femoris long，BFL)、股二頭肌短頭(biceps femoris short，BFS)、腓腸肌(gastrocnemius，GA)、比目魚肌(soleus，SO)、脛骨前肌(tibialis anterior，TA)，如圖 1中肌肉直線1~9所示。
圖1
下肢骨骼肌肉簡化模型
1:臀大肌(GM); 2:髂肌(IL); 3:股直肌(RF); 4:股側肌群(VA); 5:股二頭肌長頭(BFL); 6:股二頭肌短頭(BFS); 7:腓腸肌(GA); 8:比目魚肌(SO); 9:脛骨前肌(TA)
Figure1. Simplified muscular-skeletal model in lower extremity1: gluteus maximus (GM); 2: iliopsoas (IL); 3: rectus femoris (RF); 4: vastus group (VA); 5: biceps femoris long (BFL); 6: biceps femoris short (BFS); 7: gastrocnemius (GA); 8: soleus (SO); 9: tibialis anterior (TA)
1.2 動力學建模
建立拉格朗日-歐拉動力學方程如下：
| $ \boldsymbol{M}\left(\boldsymbol{\theta} \right)\boldsymbol{\ddot{\theta }}=\boldsymbol{B}\left(\boldsymbol{\theta}, \boldsymbol{\dot{\theta }} \right){{{\boldsymbol{\dot{\theta }}}}^{2}}+\boldsymbol{G}\left(\boldsymbol{\theta} \right)+{{\boldsymbol{\tau} }_{d}}\left(\boldsymbol{\theta}, \boldsymbol{\dot{\theta }} \right)+{{\boldsymbol{\tau} }_{m}}, $ |
| $ {{\boldsymbol{\tau }}_{m}}=\boldsymbol{R}\left(\boldsymbol{\theta} \right){{\boldsymbol{F}}_{m}}, $ |
其中θ=[θ1 θ2 θ3]T、= 1, 2, 3]T、=[1, 2, 3]T分別為髖關節、膝關節、踝關節變化的角位移、角速度和角加速度矢量(3×1)；M(θ)為系統慣性質量矩陣(3×3)；B(θ，) 2為Coriolis力和離心力的作用項(3×1)；G(θ)為重力矢量(3×1)；τ d(θ, )為系統所受的外力矩矢量(3×1)；τ m=[τHτKτA]T為肌力矢量F m產生的關節力矩矢量(3×1)；R(θ)為各肌肉的力臂矩陣矢量(3×9)。
2 多目標粒子群優化算法
2.1 多目標優化問題
多目標優化問題的研究對象是多個目標函數在給定區域上的最優化問題。對于帶有約束的多目標優化問題可定義為：
| $ s.t.\text{ }\begin{matrix} {{h}_{1}}\left(x \right)=0, \ \ \ \ \ \ i=1, 2, \cdots, p \\ {{g}_{j}}\left(x \right)\le 0, \ \ \ \ \ \ j=1, 2, \cdots, q \\ \end{matrix}, $ |
其中f(x)為目標函數，hi(x)為等式約束矢量；gj(x)為不等式約束矢量；x=(x1, x2, …, xm)∈Rm稱為m維決策向量，將滿足所有約束條件的解空間S稱為式(3)的可行域。當n=1時，式(3)為單目標優化問題；當n＞1時，式(3)為多目標優化問題[11]。
多目標優化問題中子目標之間通過決策變量相互制約，優化一個子目標可能會導致其他子目標性能降低，因此很難客觀地評價多目標問題解的優劣性。多目標優化問題的求解不同于單目標優化問題。單目標優化問題的最優值只有一個，多目標優化問題的解不是惟一的，而是存在一個最優解集合，集合中的元素稱為Pareto最優或非劣最優[12]。
2.2 目標函數及約束條件
人體恢復站立平衡過程中下肢肌力的優化求解，不僅要求有效的優化算法，還要選取適合人體生理要求的優化目標，并且肌力還需滿足人體相關運動的運動學和動力學約束條件。基于本研究內容，選取的目標函數需要滿足兩個條件：①能夠清晰地反應該過程各階段活躍肌肉的數量，②能夠較好地顯示出該過程中各肌力的分布及變化規律。參考國內外相關文獻，選取常用的三種肌力優化求解目標函數為
| $ {{f}_{1}}={{\sum\limits_{i}{\left({{F}_{i}}/{{S}_{i}} \right)}}^{2}}, $ |
| $ {{f}_{2}}=\sum\limits_{i}{{{F}_{i}}}, $ |
| $ {{f}_{3}}=\max \left({{F}_{i}}/{{S}_{i}} \right), $ |
其中Si為第i塊肌肉的生理截面積，Fi為第i塊肌肉產生的肌力。
通常選取所消耗的能量最小或肌肉疲勞程度最小作為目標函數。Crownins1ield等[13]提出了以肌肉應力平方和(見式(4))最小來度量肌肉疲勞度，這也是最常用的評價標準。Raikova等以式(4)為目標函數優化求解了腳踏板曲柄運動過程中下肢肌力的分布情況。式(5)表示肌力總和，它能夠預測出當活躍的肌肉群數量最小時各肌力的分布變化。Yoshioka等[14]以此為目標函數預測出人體由坐姿到站立過程中下肢肌力的變化。式(6)表示某組肌肉的最大應力，Bensghaier等將其作為目標函數之一優化求解了拇指與食指捏起物品過程中肌力的分布與變化規律。
此外，根據動力學方程(見式(1))和力矩平衡原理(見式(2))，肌力最優解還需滿足以下約束條件，即
| $ {{h}_{1}}={{F}_{1}}{{l}_{1}}-{{F}_{2}}{{l}_{2}}+{{F}_{3}}{{l}_{3H}}-{{F}_{4}}{{l}_{4H}}-{{\tau }_{H}}=0, $ |
| $ {{h}_{2}}={{F}_{3}}{{l}_{3K}}-{{F}_{4}}{{l}_{4K}}-{{F}_{5}}{{l}_{5}}+{{F}_{6}}{{l}_{6}}-{{F}_{7}}{{l}_{7K}}-{{\tau }_{K}}=0, $ |
| $ {{h}_{3}}={{F}_{7}}{{l}_{7A}}+{{F}_{8}}{{l}_{8}}-{{F}_{9}}{{l}_{9}}-{{\tau }_{A}}=0, $ |
| $ {{g}_{i}}:0\le {{F}_{i}}\le {{F}_{i\max }}\left(i=1, 2, \cdots, 9 \right), $ |
其中τH、τK、τA分別為髖關節、膝關節和踝關節對應的關節凈力矩(見圖 1)；li為肌力Fi對應的力臂；Fimax為各肌肉組產生肌力的最大值。
2.3 算法步驟
MOPSO算法步驟如下：
(1)初始化系統參數，包括：慣性權重系數wmin、wmax，學習因子c1、c2，最大迭代次數kmax，罰因子r0，收斂精度ε。
(2)利用混合罰函數化簡約束，構造如式(11)所示的增廣目標函數。
| $ F\left(\boldsymbol{x} \right)=f\left(\boldsymbol{x} \right)+{{r}_{0}}\sum\limits_{j=1}^{m}{\frac{1}{{{g}_{j}}\left(\boldsymbol{x} \right)}+\frac{1}{\sqrt{{{r}_{0}}}}}\sum\limits_{i=1}^{l}{{{\left[{{h}_{i}}\left(\boldsymbol{x} \right)\right]}^{2}}}, $ |
其中x為所求的肌力矢量，x=[F1 F1…F9]；f(x)為目標函數；k為迭代次數。
(3)初始化粒子個數N，隨機賦予每個粒子初始速度vi, j(0)和初始位置xi, j(0)，保證每個粒子均為混合罰函數的內點。
(4)計算初始粒子的增廣目標函數值，并將其值作為粒子的局部最優值pbest。所有pbest中最好的作為本代粒子的全局最優值gbest。
(5)k=k+1，生成下一代粒子。
(6)根據式(12)更新慣性權重系數，由式(13)、(14)計算新一代粒子的速度和位置。
| $w = {w_{\max }} - \frac{{{w_{\max }} - {w_{\min }}}}{{{k_{\max }}}}k, $ |
| $ \begin{align} & {{v}_{i, j}}\left(k+1 \right)=w{{v}_{i, j}}\left(k \right)+{{c}_{1}}{{r}_{1}}\left(k \right)\left[{{p}_{i, j}}-{{x}_{i, j}}\left(k \right)\right] \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{c}_{2}}{{r}_{2}}\left(k \right)\left[{{p}_{g, j}}-{{x}_{i, j}}\left(k \right)\right], \end{align} $ |
| $ {{x}_{i, j}}\left(k+1 \right)-{{x}_{i, j}}\left(k \right)+{{v}_{i, j}}\left(k+1 \right), $ |
其中r1(k)、r2(k)是[0, 1]上的隨機數，在每次迭代時均隨機生成；pi, j、pg, j分別對應上一代粒子的pbest和gbest；k為當前的迭代次數。
(7)計算粒子更新后的增廣目標函數值，并與上一代pbest進行比較，得到本代pbest，將本代最優作為本代gbest。
(8)搜索全部粒子，確定并篩選Pareto最優值。
(9)若滿足停止條件，則停止搜索，輸出全局最優位置，否則返回步驟(5)。
3 仿真結果及分析
本文輸入的肌肉參數包括肌肉的生理截面積S及各肌力對應的力臂l，這些數據取自文獻[3, 15]，具體數值如表 1所示。平臺傳動時人體恢復站立平衡過程的動作姿態如圖 2(a)所示，圖 2(b)給出了該過程中髖關節、膝關節和踝關節扭矩變化數據[16]。經過多次仿真實驗確定，優化算法的各參數分別取wmin=0.4, wmax=0.9，c1=c2=2, kmax=300，r0=100，ε=1×10-4, N=200。
圖2
恢復站立平衡人體姿態及關節扭矩變化曲線
(a)人體姿態變化；(b)踝、膝、髖關節扭矩變化曲線
Figure2. Postures and joint torque curves during standing balance(a) human postures; (b) the torque curves of ankle, knee and hip
表1
下肢各肌肉的生理截面積及力臂
Table1.
Physiological cross-sectional area and moment arm values of muscles in lower extremity
3.1 單目標粒子群優化算法
對于前文給出的三種不同目標函數 f1、f2、f3，取其線性加權和作為評價函數如式(15)所示，即
| $ \Phi \left({{f}_{1}}, {{f}_{2}}, {{f}_{3}} \right)={{\lambda }_{1}}{{f}_{1}}+{{\lambda }_{2}}{{f}_{2}}+{{\lambda }_{3}}{{f}_{3}}, $ |
其中λi≥0(i=1, 2, 3)表示權系數，且。不同的權系數經單目標粒子群優化(single obje-ctive particle swarm optimization, SPSO)算法求解出的肌力變化也不同，本文采用的4組權值如表 2所示。利用PSO算法分別對不同權值組成的目標函數進行優化仿真，得到各目標函數的變化分布(見圖 3)和與之對應的下肢各肌力的F-t變化曲線(見圖 4)。
表2
單目標粒子群優化算法權值
Table2.
Weights of SPSO
圖3
各目標函數值變化分布圖
Figure3.
Distributions of objective functions
圖4
SPSO肌力優化結果
Figure4.
Muscular forces of SPSO
站立平衡恢復過程主要包括緩沖(0~0.4 s)和直立(0.4~1.2 s)兩個階段。緩沖階段下肢各關節均會產生一定程度的彎曲用以消減平臺傳動所帶來的沖擊力。直立階段是當人體適應平臺傳動速度后各關節在相關肌肉的協同作用下恢復正常狀態，即恢復站立平衡。由圖 4可以看出，該過程中活躍的肌肉為GM、RF、VA、GA、SO和TA，但顯然能夠發現采用線性加權和的形式確定目標函數對權值的選取依賴性較大，采用不同的目標函數得到的肌力結果存在較明顯的差異。以f1為目標函數得到的仿真結果(見SPSO1圖)能夠較好地顯示各肌力分布和變化規律，但該結果顯示處于活躍的肌肉主要是GM和VA。故其不能有效地明確各階段肌肉激活數量。如SPSO2圖所示，當以f2為目標函數時，各個肌力的變化范圍較小，由于優化求解肌力和最小時得到的解分別對應各肌力的最小值。這就使得其求解的肌力大小與其他情況出現數量級上的差異，不能有效地預測人體恢復站立平衡過程中下肢肌力的分布與變化規律，文獻[6]也將f2作為目標函數之一用于預測手指捏起物體過程中手部肌力的變化，當其獨自作為目標函數得到的優化結果同樣不符合引證文獻[17]中的數據范圍。以f3為目標函數得到的仿真結果(見SPSO3圖)能夠較好地顯示出各階段活躍肌肉的數量，但由SPSO3圖顯然能夠發現RF、VA、GA和SO的肌力變化幾乎一致，不能有效地預測下肢各肌肉的肌力分布和恢復站立平衡過程中發揮的作用。綜上所述，對于人體恢復站立平衡這一動態過程而言，同時以f1和f3作為預測肌力變化的目標函數，能夠清晰地顯示該過程中激活肌肉的數量，并有效地預測出下肢各肌肉肌力的分布和變化規律，更符合人體運動學和生理學的實際情況。
3.2 多目標粒子群優化算法
人體站立平衡調節是在多關節和多肌肉群協調作用下完成的一種生理機能。由于關節和肌肉的不可見性，導致預測該調節過程中下肢肌力的變化存在諸多不確定性因素。所以，采用單目標作為評價函數就存在一定的局限，基于上節對單目標優化結果的分析，本節同時以f1和f3為目標函數，利用PSO算法求解該多目標問題。
由圖 2可知，在t=0.5 s附近時髖關節、膝關節和踝關節的扭矩變化最大，此時下肢各肌肉也處于較活躍的階段。利用MOPSO算法優化仿真得到Pareto最優曲線如圖 5(a)所示。由圖 5(a)可以看出，采用MOPSO算法得到一個非支配解集，其趨向于由f1和f3決定的Pareto最優邊界。在該邊界選取三組最優值SolutionⅠ、SolutionⅡ和SolutionⅢ，對應坐標分別為(449.9, 381.8)、(472.2, 148.9)和(1858.1, 56.1)。此時下肢9組肌力的分布變化如圖 5(b)所示。
圖5
t=0.5 s時Pareto最優曲線及相應肌力變化分布
a)Pareto最優曲線；(b)肌力變化直方圖
Figure5. Pareto optimal curve and muscular forces of t=0.5 s(a) the Pareto optimal curve; (b) the histogram of muscular forces
當t=0.5 s時，如圖 5(b)所示，主要活動的肌肉是GM、RF、VA和SO，這與利用單目標優化得到的仿真結果基本相同。同時優化多目標(f1，f3)得到的非支配解SolutionⅠ和SolutionⅢ趨向于分別以f1和f3為目標函數得到的肌力變化，而SolutionⅡ能夠合理地取到二者同時最小的Pareto最優值，而非單一的絕對最優值。這也是多目標優化方法突出的優點，有效地改善了單目標優化結果的單一性和絕對性。人體恢復站立平衡過程中經MOPSO算法優化求解所得Pareto最優曲線和肌力分布變化規律如圖 6所示。
圖6
Pareto最優曲線及肌力變化分布
(a)Pareto最優曲線；(b)MOPSO肌力變化分布
Figure6. Pareto optimal curves and muscular forces(a) the Pareto optimal curve; (b) the muscular forces of MOPSO
由圖 6(a)可知，人體恢復站立平衡過程的每個時刻均會得到一組相似的非支配解集，其均趨向于由f1和f3決定的Pareto最優邊界，選取每一時刻得到的最優值如圖 6(a)中“○”標記所示。經MOPSO算法優化求解得到的下肢肌力分布及變化規律如圖 6(b)所示。由圖 6(b)可以看出，以f1和f3同時為目標函數，采用MOPSO優化求解人體恢復站立平衡過程中下肢各肌力變化明顯，該過程主要活動的肌肉為TA、GM、RF、VA、GA和SO，其他肌肉群肌力變化不明顯。本文采用八通道有線表面肌電儀采集平臺低速傳動時試驗者該6塊肌肉的肌電信號，經放大、帶通濾波后得到的表面肌電圖(electro-myography，EMG)如圖 7所示。
由圖 7可以看出，采集的該6塊肌肉的肌電信號與經MOPSO算法優化求解得到的相應肌力變化趨勢基本一致。結合圖 6和圖 7顯然能夠看出，緩沖階段主要活動的肌肉為VA、GA、SO和TA，其中VA、SO和TA在人體受到傳動平臺瞬時作用時肌力迅速增加，GA和TA迅速收縮分別調節膝關節和踝關節發生一定程度的彎曲用于消減所受到的平臺沖擊力。與單目標優化結果相似，伸膝肌群VA和伸踝肌群SO也做退讓性收縮以對抗膝關節和踝關節的彎曲，但比較圖 6(b)和圖 4的肌力優化結果可知，采用MOPSO算法得到的VA和SO肌力變化，在緩沖階段0.3~0.4 s時，VA變化迅速減小、SO變化趨于平緩，二者均為直立階段做好準備，這也更合乎人體生理調節機制。IL增加的幅度比較小且緩慢。因此，緩沖階段主要由膝關節和踝關節完成，而髖關節基本未參與緩沖動作。
圖7
活躍肌肉的肌力及EMG變化
Figure7.
Muscular force and EMG of the active muscles
直立階段主要活動的肌肉為GM、RF、VA、GA和SO，其他肌力變化不明顯。GM和RF的肌力在該階段逐漸增大，分別調節髖關節和膝關節恢復正常狀態。VA同時也迅速增大，在t=0.6 s附近達到峰值，做退讓性收縮以恢復膝關節的彎曲狀態。SO的肌力迅速增大，調節踝關節恢復正常狀態，其在t=0.55 s時達到峰值后迅速減弱，待降至350~400 N時肌力變化趨于平緩，保證其它關節平衡調節時的穩定性。在單目標優化的直立階段，GA活動減弱，肌力變化不明顯，而由圖 6(b)可以看出，在t=0.6 s時，GA的肌力迅速增加，配合SO調節踝關節恢復正常狀態。這幾組肌肉群在人體恢復站立平衡后均逐漸減小直至恢復常態。經分析可知，采用MOPSO算法可以清晰的得到人體恢復站立平衡過程中各階段下肢活動肌肉的數量和相應肌力的分布及變化規律。
目前肌肉力理論計算使用最廣泛的方法主要分為兩類：①基于反向動力學的靜態優化方法；②動態優化方法。本文采用的模型和方法屬于前者統籌的范圍，由于所建模型的純力學特性，也就導致其應用存在一定的局限性。它能夠較好的用于宏觀的肌力求解問題，但不能夠精確得獲得肌肉生理和力學特性模型。而且使用基于反向動力學的靜態優化方法進行肌肉力的預測時，人體運動參數參與到整個計算過程，所以，肌肉力的預測值對運動測量參數比較敏感。Anderson等[18]用優化的方法研究人體完整步態，比較了靜態優化方法和動態優化方法預測得到的肌力，得出結論，兩種方法等價，可以交互使用。并且進一步指出，反向動力學分析方法適用于運動學參數易于測量的常規運動。對于那些可能對人體造成嚴重傷害，無法獲知準確運動學參數的運動，采用前向動力學模擬的方法。所以，本文提出的模型和方法適用于諸如步態、跳躍等常規運動的肌力理論預測。另外，如何將該力學模型與肌肉的生物模型結合一起來預測肌力的變化有待進一步的研究學習。
4 結語
本文主要研究了被動運動下人體恢復站立平衡過程中下肢主要肌肉的肌力分布與變化情況。首先，建立了人體下肢3-DOF平面運動的仿真簡化肌骨模型，以及其數學模型；其次，分析了多目標最優化問題中子目標函數的選取方法，并求解出該過程中肌力的變化規律；最后，分析所得到的優化結果，討論了各肌肉群在恢復站立平衡過程中發揮的作用及各關節與對應肌肉群的協調機制。
該模型和算法能夠應用于其他運動過程的肌力求解，為冗余肌力優化求解問題提供了一種新方法。預測人體恢復站立平衡過程中肌力的變化，不僅能夠對人體姿態平衡醫療康復起到一定的指導作用，而且能為進一步研究人體神經對肌肉的協調控制、建立完整的神經—肌肉模型提供理論依據。
引言
人體站立平衡機能一直是生命科學研究的熱點和難點，其平衡過程是在小腦的控制支配下，通過神經系統傳遞生物電信號刺激肌肉收縮，牽拉肌肉所附著的骨骼繞關節旋轉，改變骨骼間的相對位置，從而協調人體恢復站立平衡。生活中由于意外傷害或者中風等疾病導致肌肉組織功能障礙，引起的平衡能力失調，導致了大量的傷害事故。因此，預測恢復站立平衡過程中肌力的分布變化能夠更好地理解其調節機制，為醫療康復起到實際指導作用。采用植入性的儀器測量人體肌肉作用力，不僅技術復雜，成本高，而且會對人體造成創傷，不適合進行正常人體肌力的研究[1]。因此，人體肌肉骨骼模型為間接求解肌力變化分布提供了新思路。由于關節力矩是在多塊肌肉的共同作用下產生的，故一般多采用約束優化方法求解冗余肌力的分布情況。針對不同問題的肌力變化規律，諸多學者采用了不同的優化算法，Neptune[2]、Raikova等[3]分別采用模擬退火(simulated annealing, SA)算法和拉格朗日乘子法對腳踏板曲柄運動過程中下肢肌力分布進行了優化求解，Seth等[4]利用神經骨肌跟蹤(neuromusculoskeletal tracking, NMT)法分析了縱跳過程中下肢肌力的變化規律。然而，上述采用的算法均只是單個目標函數的肌力優化求解。
迄今為止，利用智能優化算法求解多目標問題被廣泛應用于軌跡規劃、系統工程和控制論等領域。Santana-Quintero等[5]針對帶有約束的多目標問題采用粗集理論與差分進化結合的算法予以求解，Bensghaier等[6]采用遺傳算法(genetic algorithm，GA)優化求解了拇指與食指捏起物品過程中肌力的分布與變化。多目標粒子群優化(multi-objective particle swarm optimization, MOPSO)算法同樣備受諸多學者的青睞，如Tan等[7]利用其求解工程約束問題，Carvalho等[8]應用該算法測試了多目標問題中支配解的收斂性和多樣性。Yang[9]利用改進后的粒子群算法解決了多目標離散優化問題。但是，目前該算法還沒有被用于預測人體恢復站立平衡過程中下肢肌力的分布與變化規律。因此，本文在建立的3-DOF平面運動仿真簡化模型的基礎上采用混合罰函數化簡約束，然后利用粒子群優化(particle swarm optimization, PSO)算法分別求解恢復站立平衡過程中單目標函數和多目標函數的下肢9組肌力最優值，進而分析各組肌力在該過程中的變化規律和發揮的作用。
1 模型的分析與建立
1.1 下肢骨骼肌肉簡化模型
本文旨在預測被動運動下人體站立姿態由失衡至恢復平衡過程中肌肉力動態變化的最優值。研究中所施加的外力具體指水平運動平臺對足底提供后向水平加速度[10]。根據此外力的作用方式和人體下肢生理結構，可近似認為腳掌始終保持與平臺平穩接觸且兩腿的運動關于矢狀面左右對稱。建立二維平面下肢骨架3自由度、4剛體(包含軀干、大腿、小腿和腳)和9塊肌肉的運動簡化物理模型如圖 1所示。參考人體解剖學數據和文獻，選取與下肢運動相關的9組主要肌肉，包括臀大肌(gluteus maximus，GM)、髂肌(iliopsoas，IL)、股直肌(rectus femoris，RF)、股側肌群(vastus group，VA)、股二頭肌長頭(biceps femoris long，BFL)、股二頭肌短頭(biceps femoris short，BFS)、腓腸肌(gastrocnemius，GA)、比目魚肌(soleus，SO)、脛骨前肌(tibialis anterior，TA)，如圖 1中肌肉直線1~9所示。
圖1
下肢骨骼肌肉簡化模型
1:臀大肌(GM); 2:髂肌(IL); 3:股直肌(RF); 4:股側肌群(VA); 5:股二頭肌長頭(BFL); 6:股二頭肌短頭(BFS); 7:腓腸肌(GA); 8:比目魚肌(SO); 9:脛骨前肌(TA)
Figure1. Simplified muscular-skeletal model in lower extremity1: gluteus maximus (GM); 2: iliopsoas (IL); 3: rectus femoris (RF); 4: vastus group (VA); 5: biceps femoris long (BFL); 6: biceps femoris short (BFS); 7: gastrocnemius (GA); 8: soleus (SO); 9: tibialis anterior (TA)
1.2 動力學建模
建立拉格朗日-歐拉動力學方程如下：
| $ \boldsymbol{M}\left(\boldsymbol{\theta} \right)\boldsymbol{\ddot{\theta }}=\boldsymbol{B}\left(\boldsymbol{\theta}, \boldsymbol{\dot{\theta }} \right){{{\boldsymbol{\dot{\theta }}}}^{2}}+\boldsymbol{G}\left(\boldsymbol{\theta} \right)+{{\boldsymbol{\tau} }_{d}}\left(\boldsymbol{\theta}, \boldsymbol{\dot{\theta }} \right)+{{\boldsymbol{\tau} }_{m}}, $ |
| $ {{\boldsymbol{\tau }}_{m}}=\boldsymbol{R}\left(\boldsymbol{\theta} \right){{\boldsymbol{F}}_{m}}, $ |
其中θ=[θ1 θ2 θ3]T、= 1, 2, 3]T、=[1, 2, 3]T分別為髖關節、膝關節、踝關節變化的角位移、角速度和角加速度矢量(3×1)；M(θ)為系統慣性質量矩陣(3×3)；B(θ，) 2為Coriolis力和離心力的作用項(3×1)；G(θ)為重力矢量(3×1)；τ d(θ, )為系統所受的外力矩矢量(3×1)；τ m=[τHτKτA]T為肌力矢量F m產生的關節力矩矢量(3×1)；R(θ)為各肌肉的力臂矩陣矢量(3×9)。
2 多目標粒子群優化算法
2.1 多目標優化問題
多目標優化問題的研究對象是多個目標函數在給定區域上的最優化問題。對于帶有約束的多目標優化問題可定義為：
| $ s.t.\text{ }\begin{matrix} {{h}_{1}}\left(x \right)=0, \ \ \ \ \ \ i=1, 2, \cdots, p \\ {{g}_{j}}\left(x \right)\le 0, \ \ \ \ \ \ j=1, 2, \cdots, q \\ \end{matrix}, $ |
其中f(x)為目標函數，hi(x)為等式約束矢量；gj(x)為不等式約束矢量；x=(x1, x2, …, xm)∈Rm稱為m維決策向量，將滿足所有約束條件的解空間S稱為式(3)的可行域。當n=1時，式(3)為單目標優化問題；當n＞1時，式(3)為多目標優化問題[11]。
多目標優化問題中子目標之間通過決策變量相互制約，優化一個子目標可能會導致其他子目標性能降低，因此很難客觀地評價多目標問題解的優劣性。多目標優化問題的求解不同于單目標優化問題。單目標優化問題的最優值只有一個，多目標優化問題的解不是惟一的，而是存在一個最優解集合，集合中的元素稱為Pareto最優或非劣最優[12]。
2.2 目標函數及約束條件
人體恢復站立平衡過程中下肢肌力的優化求解，不僅要求有效的優化算法，還要選取適合人體生理要求的優化目標，并且肌力還需滿足人體相關運動的運動學和動力學約束條件。基于本研究內容，選取的目標函數需要滿足兩個條件：①能夠清晰地反應該過程各階段活躍肌肉的數量，②能夠較好地顯示出該過程中各肌力的分布及變化規律。參考國內外相關文獻，選取常用的三種肌力優化求解目標函數為
| $ {{f}_{1}}={{\sum\limits_{i}{\left({{F}_{i}}/{{S}_{i}} \right)}}^{2}}, $ |
| $ {{f}_{2}}=\sum\limits_{i}{{{F}_{i}}}, $ |
| $ {{f}_{3}}=\max \left({{F}_{i}}/{{S}_{i}} \right), $ |
其中Si為第i塊肌肉的生理截面積，Fi為第i塊肌肉產生的肌力。
通常選取所消耗的能量最小或肌肉疲勞程度最小作為目標函數。Crownins1ield等[13]提出了以肌肉應力平方和(見式(4))最小來度量肌肉疲勞度，這也是最常用的評價標準。Raikova等以式(4)為目標函數優化求解了腳踏板曲柄運動過程中下肢肌力的分布情況。式(5)表示肌力總和，它能夠預測出當活躍的肌肉群數量最小時各肌力的分布變化。Yoshioka等[14]以此為目標函數預測出人體由坐姿到站立過程中下肢肌力的變化。式(6)表示某組肌肉的最大應力，Bensghaier等將其作為目標函數之一優化求解了拇指與食指捏起物品過程中肌力的分布與變化規律。
此外，根據動力學方程(見式(1))和力矩平衡原理(見式(2))，肌力最優解還需滿足以下約束條件，即
| $ {{h}_{1}}={{F}_{1}}{{l}_{1}}-{{F}_{2}}{{l}_{2}}+{{F}_{3}}{{l}_{3H}}-{{F}_{4}}{{l}_{4H}}-{{\tau }_{H}}=0, $ |
| $ {{h}_{2}}={{F}_{3}}{{l}_{3K}}-{{F}_{4}}{{l}_{4K}}-{{F}_{5}}{{l}_{5}}+{{F}_{6}}{{l}_{6}}-{{F}_{7}}{{l}_{7K}}-{{\tau }_{K}}=0, $ |
| $ {{h}_{3}}={{F}_{7}}{{l}_{7A}}+{{F}_{8}}{{l}_{8}}-{{F}_{9}}{{l}_{9}}-{{\tau }_{A}}=0, $ |
| $ {{g}_{i}}:0\le {{F}_{i}}\le {{F}_{i\max }}\left(i=1, 2, \cdots, 9 \right), $ |
其中τH、τK、τA分別為髖關節、膝關節和踝關節對應的關節凈力矩(見圖 1)；li為肌力Fi對應的力臂；Fimax為各肌肉組產生肌力的最大值。
2.3 算法步驟
MOPSO算法步驟如下：
(1)初始化系統參數，包括：慣性權重系數wmin、wmax，學習因子c1、c2，最大迭代次數kmax，罰因子r0，收斂精度ε。
(2)利用混合罰函數化簡約束，構造如式(11)所示的增廣目標函數。
| $ F\left(\boldsymbol{x} \right)=f\left(\boldsymbol{x} \right)+{{r}_{0}}\sum\limits_{j=1}^{m}{\frac{1}{{{g}_{j}}\left(\boldsymbol{x} \right)}+\frac{1}{\sqrt{{{r}_{0}}}}}\sum\limits_{i=1}^{l}{{{\left[{{h}_{i}}\left(\boldsymbol{x} \right)\right]}^{2}}}, $ |
其中x為所求的肌力矢量，x=[F1 F1…F9]；f(x)為目標函數；k為迭代次數。
(3)初始化粒子個數N，隨機賦予每個粒子初始速度vi, j(0)和初始位置xi, j(0)，保證每個粒子均為混合罰函數的內點。
(4)計算初始粒子的增廣目標函數值，并將其值作為粒子的局部最優值pbest。所有pbest中最好的作為本代粒子的全局最優值gbest。
(5)k=k+1，生成下一代粒子。
(6)根據式(12)更新慣性權重系數，由式(13)、(14)計算新一代粒子的速度和位置。
| $w = {w_{\max }} - \frac{{{w_{\max }} - {w_{\min }}}}{{{k_{\max }}}}k, $ |
| $ \begin{align} & {{v}_{i, j}}\left(k+1 \right)=w{{v}_{i, j}}\left(k \right)+{{c}_{1}}{{r}_{1}}\left(k \right)\left[{{p}_{i, j}}-{{x}_{i, j}}\left(k \right)\right] \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{c}_{2}}{{r}_{2}}\left(k \right)\left[{{p}_{g, j}}-{{x}_{i, j}}\left(k \right)\right], \end{align} $ |
| $ {{x}_{i, j}}\left(k+1 \right)-{{x}_{i, j}}\left(k \right)+{{v}_{i, j}}\left(k+1 \right), $ |
其中r1(k)、r2(k)是[0, 1]上的隨機數，在每次迭代時均隨機生成；pi, j、pg, j分別對應上一代粒子的pbest和gbest；k為當前的迭代次數。
(7)計算粒子更新后的增廣目標函數值，并與上一代pbest進行比較，得到本代pbest，將本代最優作為本代gbest。
(8)搜索全部粒子，確定并篩選Pareto最優值。
(9)若滿足停止條件，則停止搜索，輸出全局最優位置，否則返回步驟(5)。
3 仿真結果及分析
本文輸入的肌肉參數包括肌肉的生理截面積S及各肌力對應的力臂l，這些數據取自文獻[3, 15]，具體數值如表 1所示。平臺傳動時人體恢復站立平衡過程的動作姿態如圖 2(a)所示，圖 2(b)給出了該過程中髖關節、膝關節和踝關節扭矩變化數據[16]。經過多次仿真實驗確定，優化算法的各參數分別取wmin=0.4, wmax=0.9，c1=c2=2, kmax=300，r0=100，ε=1×10-4, N=200。
圖2
恢復站立平衡人體姿態及關節扭矩變化曲線
(a)人體姿態變化；(b)踝、膝、髖關節扭矩變化曲線
Figure2. Postures and joint torque curves during standing balance(a) human postures; (b) the torque curves of ankle, knee and hip
表1
下肢各肌肉的生理截面積及力臂
Table1.
Physiological cross-sectional area and moment arm values of muscles in lower extremity
3.1 單目標粒子群優化算法
對于前文給出的三種不同目標函數 f1、f2、f3，取其線性加權和作為評價函數如式(15)所示，即
| $ \Phi \left({{f}_{1}}, {{f}_{2}}, {{f}_{3}} \right)={{\lambda }_{1}}{{f}_{1}}+{{\lambda }_{2}}{{f}_{2}}+{{\lambda }_{3}}{{f}_{3}}, $ |
其中λi≥0(i=1, 2, 3)表示權系數，且。不同的權系數經單目標粒子群優化(single obje-ctive particle swarm optimization, SPSO)算法求解出的肌力變化也不同，本文采用的4組權值如表 2所示。利用PSO算法分別對不同權值組成的目標函數進行優化仿真，得到各目標函數的變化分布(見圖 3)和與之對應的下肢各肌力的F-t變化曲線(見圖 4)。
表2
單目標粒子群優化算法權值
Table2.
Weights of SPSO
圖3
各目標函數值變化分布圖
Figure3.
Distributions of objective functions
圖4
SPSO肌力優化結果
Figure4.
Muscular forces of SPSO
站立平衡恢復過程主要包括緩沖(0~0.4 s)和直立(0.4~1.2 s)兩個階段。緩沖階段下肢各關節均會產生一定程度的彎曲用以消減平臺傳動所帶來的沖擊力。直立階段是當人體適應平臺傳動速度后各關節在相關肌肉的協同作用下恢復正常狀態，即恢復站立平衡。由圖 4可以看出，該過程中活躍的肌肉為GM、RF、VA、GA、SO和TA，但顯然能夠發現采用線性加權和的形式確定目標函數對權值的選取依賴性較大，采用不同的目標函數得到的肌力結果存在較明顯的差異。以f1為目標函數得到的仿真結果(見SPSO1圖)能夠較好地顯示各肌力分布和變化規律，但該結果顯示處于活躍的肌肉主要是GM和VA。故其不能有效地明確各階段肌肉激活數量。如SPSO2圖所示，當以f2為目標函數時，各個肌力的變化范圍較小，由于優化求解肌力和最小時得到的解分別對應各肌力的最小值。這就使得其求解的肌力大小與其他情況出現數量級上的差異，不能有效地預測人體恢復站立平衡過程中下肢肌力的分布與變化規律，文獻[6]也將f2作為目標函數之一用于預測手指捏起物體過程中手部肌力的變化，當其獨自作為目標函數得到的優化結果同樣不符合引證文獻[17]中的數據范圍。以f3為目標函數得到的仿真結果(見SPSO3圖)能夠較好地顯示出各階段活躍肌肉的數量，但由SPSO3圖顯然能夠發現RF、VA、GA和SO的肌力變化幾乎一致，不能有效地預測下肢各肌肉的肌力分布和恢復站立平衡過程中發揮的作用。綜上所述，對于人體恢復站立平衡這一動態過程而言，同時以f1和f3作為預測肌力變化的目標函數，能夠清晰地顯示該過程中激活肌肉的數量，并有效地預測出下肢各肌肉肌力的分布和變化規律，更符合人體運動學和生理學的實際情況。
3.2 多目標粒子群優化算法
人體站立平衡調節是在多關節和多肌肉群協調作用下完成的一種生理機能。由于關節和肌肉的不可見性，導致預測該調節過程中下肢肌力的變化存在諸多不確定性因素。所以，采用單目標作為評價函數就存在一定的局限，基于上節對單目標優化結果的分析，本節同時以f1和f3為目標函數，利用PSO算法求解該多目標問題。
由圖 2可知，在t=0.5 s附近時髖關節、膝關節和踝關節的扭矩變化最大，此時下肢各肌肉也處于較活躍的階段。利用MOPSO算法優化仿真得到Pareto最優曲線如圖 5(a)所示。由圖 5(a)可以看出，采用MOPSO算法得到一個非支配解集，其趨向于由f1和f3決定的Pareto最優邊界。在該邊界選取三組最優值SolutionⅠ、SolutionⅡ和SolutionⅢ，對應坐標分別為(449.9, 381.8)、(472.2, 148.9)和(1858.1, 56.1)。此時下肢9組肌力的分布變化如圖 5(b)所示。
圖5
t=0.5 s時Pareto最優曲線及相應肌力變化分布
a)Pareto最優曲線；(b)肌力變化直方圖
Figure5. Pareto optimal curve and muscular forces of t=0.5 s(a) the Pareto optimal curve; (b) the histogram of muscular forces
當t=0.5 s時，如圖 5(b)所示，主要活動的肌肉是GM、RF、VA和SO，這與利用單目標優化得到的仿真結果基本相同。同時優化多目標(f1，f3)得到的非支配解SolutionⅠ和SolutionⅢ趨向于分別以f1和f3為目標函數得到的肌力變化，而SolutionⅡ能夠合理地取到二者同時最小的Pareto最優值，而非單一的絕對最優值。這也是多目標優化方法突出的優點，有效地改善了單目標優化結果的單一性和絕對性。人體恢復站立平衡過程中經MOPSO算法優化求解所得Pareto最優曲線和肌力分布變化規律如圖 6所示。
圖6
Pareto最優曲線及肌力變化分布
(a)Pareto最優曲線；(b)MOPSO肌力變化分布
Figure6. Pareto optimal curves and muscular forces(a) the Pareto optimal curve; (b) the muscular forces of MOPSO
由圖 6(a)可知，人體恢復站立平衡過程的每個時刻均會得到一組相似的非支配解集，其均趨向于由f1和f3決定的Pareto最優邊界，選取每一時刻得到的最優值如圖 6(a)中“○”標記所示。經MOPSO算法優化求解得到的下肢肌力分布及變化規律如圖 6(b)所示。由圖 6(b)可以看出，以f1和f3同時為目標函數，采用MOPSO優化求解人體恢復站立平衡過程中下肢各肌力變化明顯，該過程主要活動的肌肉為TA、GM、RF、VA、GA和SO，其他肌肉群肌力變化不明顯。本文采用八通道有線表面肌電儀采集平臺低速傳動時試驗者該6塊肌肉的肌電信號，經放大、帶通濾波后得到的表面肌電圖(electro-myography，EMG)如圖 7所示。
由圖 7可以看出，采集的該6塊肌肉的肌電信號與經MOPSO算法優化求解得到的相應肌力變化趨勢基本一致。結合圖 6和圖 7顯然能夠看出，緩沖階段主要活動的肌肉為VA、GA、SO和TA，其中VA、SO和TA在人體受到傳動平臺瞬時作用時肌力迅速增加，GA和TA迅速收縮分別調節膝關節和踝關節發生一定程度的彎曲用于消減所受到的平臺沖擊力。與單目標優化結果相似，伸膝肌群VA和伸踝肌群SO也做退讓性收縮以對抗膝關節和踝關節的彎曲，但比較圖 6(b)和圖 4的肌力優化結果可知，采用MOPSO算法得到的VA和SO肌力變化，在緩沖階段0.3~0.4 s時，VA變化迅速減小、SO變化趨于平緩，二者均為直立階段做好準備，這也更合乎人體生理調節機制。IL增加的幅度比較小且緩慢。因此，緩沖階段主要由膝關節和踝關節完成，而髖關節基本未參與緩沖動作。
圖7
活躍肌肉的肌力及EMG變化
Figure7.
Muscular force and EMG of the active muscles
直立階段主要活動的肌肉為GM、RF、VA、GA和SO，其他肌力變化不明顯。GM和RF的肌力在該階段逐漸增大，分別調節髖關節和膝關節恢復正常狀態。VA同時也迅速增大，在t=0.6 s附近達到峰值，做退讓性收縮以恢復膝關節的彎曲狀態。SO的肌力迅速增大，調節踝關節恢復正常狀態，其在t=0.55 s時達到峰值后迅速減弱，待降至350~400 N時肌力變化趨于平緩，保證其它關節平衡調節時的穩定性。在單目標優化的直立階段，GA活動減弱，肌力變化不明顯，而由圖 6(b)可以看出，在t=0.6 s時，GA的肌力迅速增加，配合SO調節踝關節恢復正常狀態。這幾組肌肉群在人體恢復站立平衡后均逐漸減小直至恢復常態。經分析可知，采用MOPSO算法可以清晰的得到人體恢復站立平衡過程中各階段下肢活動肌肉的數量和相應肌力的分布及變化規律。
目前肌肉力理論計算使用最廣泛的方法主要分為兩類：①基于反向動力學的靜態優化方法；②動態優化方法。本文采用的模型和方法屬于前者統籌的范圍，由于所建模型的純力學特性，也就導致其應用存在一定的局限性。它能夠較好的用于宏觀的肌力求解問題，但不能夠精確得獲得肌肉生理和力學特性模型。而且使用基于反向動力學的靜態優化方法進行肌肉力的預測時，人體運動參數參與到整個計算過程，所以，肌肉力的預測值對運動測量參數比較敏感。Anderson等[18]用優化的方法研究人體完整步態，比較了靜態優化方法和動態優化方法預測得到的肌力，得出結論，兩種方法等價，可以交互使用。并且進一步指出，反向動力學分析方法適用于運動學參數易于測量的常規運動。對于那些可能對人體造成嚴重傷害，無法獲知準確運動學參數的運動，采用前向動力學模擬的方法。所以，本文提出的模型和方法適用于諸如步態、跳躍等常規運動的肌力理論預測。另外，如何將該力學模型與肌肉的生物模型結合一起來預測肌力的變化有待進一步的研究學習。
4 結語
本文主要研究了被動運動下人體恢復站立平衡過程中下肢主要肌肉的肌力分布與變化情況。首先，建立了人體下肢3-DOF平面運動的仿真簡化肌骨模型，以及其數學模型；其次，分析了多目標最優化問題中子目標函數的選取方法，并求解出該過程中肌力的變化規律；最后，分析所得到的優化結果，討論了各肌肉群在恢復站立平衡過程中發揮的作用及各關節與對應肌肉群的協調機制。
該模型和算法能夠應用于其他運動過程的肌力求解，為冗余肌力優化求解問題提供了一種新方法。預測人體恢復站立平衡過程中肌力的變化，不僅能夠對人體姿態平衡醫療康復起到一定的指導作用，而且能為進一步研究人體神經對肌肉的協調控制、建立完整的神經—肌肉模型提供理論依據。
