應用功能性電刺激(FES)技術, 設計了兩自由度的病理性震顫抑制系統, 能夠實現對人體腕關節內外翻和上下切自主運動過程中伴隨的震顫進行抑制。首先建立了腕關節的肌骨模型, 針對參與腕關節內外翻和上下切運動的四塊主要肌肉(橈側腕長伸肌、橈側腕屈肌、尺側腕伸肌、尺側腕屈肌)進行建模, 同時建立了腕關節的二階機械阻抗動力學模型。然后以腕關節肌骨模型為控制對象設計了FES抑震控制器, 提出了基于逆模型的前饋控制與線性二次型調節器(LQR)最優控制相結合的復合控制策略。在此基礎上仿真了在階躍、正弦以及真實患者數據3種輸入作為自主運動的情況下, 系統對于震顫的抑制, 結果表明該FES抑震控制器能夠對多自由度震顫幅度實現94%以上的抑制。
引用本文： 張偉, 張定國, 劉建榮. 基于功能性電刺激的腕關節多自由度震顫抑制仿真研究. 生物醫學工程學雜志, 2015, 32(2): 423-429. doi: 10.7507/1001-5515.20150076 復制
引言
震顫是指身體某些運動功能區有節律、非自主的抖動，主要是由于神經系統對于肌肉的往復支配引起的肌肉周期性的收縮[1]。根據臨床表現，病理性震顫可以分為靜止性震顫、姿勢性震顫以及運動性震顫[2]。震顫減弱了患者處理日常生活的能力，大大降低了患者的生活質量[3]。目前藥物治療和手術治療是治療震顫的兩種主要方式，藥物治療雖暫時有所好轉但并不能阻止病情發展；手術治療也有局限，比如典型的深部腦刺激技術[4]，只是對某些病癥引起的病理性震顫有效，而且費用高、風險大、副作用明顯。據相關研究表明至少有25%的震顫患者無法通過藥物治療或者手術治療得到治愈[5]。正是由于傳統治療遇到的這些問題和瓶頸，越來越多的工程科研人員加入到病態震顫抑制這個領域，從而使輔助式抑震技術為病理性震顫的治療提供了新的選擇。
功能性電刺激(functional electrical stimulation, FES)是指用弱電流脈沖刺激人體的骨骼肌，肌肉因受到刺激收縮產生力以驅動關節運動，從而實現運動功能恢復和重建的目的，目前被大量應用于康復工程領域。加拿大的Prochazka團隊首次提出將FES作為一種重要的抑制病理性震顫的技術手段，他們利用FES異相刺激目標肌肉以期抑制肘關節、腕關節的病態顫抖[5]。法國的Poignet團隊提出了基于FES技術控制關節剛度來抑制顫抖的方法，設計了一種使用FES控制關節拮抗肌從而增強關節剛度的顫抖抑制系統[6]。筆者前期也開展了FES抑震研究，通過基于神經振蕩器的控制方式分別實現了對腕部和肘部的震顫抑制，取得了良好的效果[7]。經調研發現，目前的國內外相關工作都是針對腕關節或者肘關節單自由度的情況進行的，而且日常生活中的動作都是基于兩到三個甚至更多自由度完成的。本文的主要工作就是在前人的基礎上首先完成腕關節內外翻和上下切兩個自由度的肌骨模型建模，并完成基于FES的腕關節兩自由度病態震顫抑制的控制器設計與抑震系統的仿真。
本文提出的基于FES的抑震控制系統如圖 1所示。系統的期望輸入為自主運動，與反饋信號做差后經過FES控制器產生異相電刺激脈沖作用在兩自由度肌骨模型上從而達到震顫抑制的目的。
圖1
基于FES的腕關節抑震控制系統框圖
Figure1.
System diagram of the wrist joint tremor suppression based on FES
我們建立的肌骨模型包含兩對拮抗肌模型和二自由度的腕關節骨骼動力學，以肌骨模型為控制對象設計了FES控制器。我們提出了反饋加前饋的控制策略，其中反饋控制采用基于狀態空間的線性二次型調節器(linear quadratic regulator, LQR)，前饋控制器采用基于肌肉逆模型[8]，仿真結果表明該復合控制策略能夠有效地對病態震顫進行抑制，簡化了控制系統的復雜程度，降低了向實際抑震系統轉化的難度。
1 肌骨模型的建立
本文建立的腕部肌骨模型主要包括四塊肌肉，它們分別是橈側腕長伸肌(extensor carpi radialis longus, ECRL)、尺側腕伸肌(extensor carpi ulnaris, ECU)、橈側腕屈肌(flexor carpi radialis, FCR)、尺側腕屈肌(flexor carpi ulnaris, FCU)。FCU、FCR主要參與腕關節的內翻，ECRL、ECU主要參與腕關節的外翻，FCR、ECRL主要參與腕關節的上切，FCU、ECU主要參與腕關節的下切[9]。
本文建立的肌骨模型主要由三個部分組成：肌肉的激活，肌肉收縮，骨骼動力特性。第一部分肌肉實際產生的激活程度Act(Nm)由電刺激脈沖的幅值強度來計算；第二部分肌肉收縮特性實際反映關節角度和關節速度與實際關節力矩之間的關系；第三部分主要是關節的機械特性，當有力矩作用在骨骼上時表現為關節的運動。
1.1 FES作用下的肌肉模型
| ${A_{{\rm{ct}}}}=\left\{ \begin{array}{l} 0\;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; {P_{\rm{A}}} \le {T_{{\rm{hres}}}}\\ \frac{{{S_{\rm{f}}} \times \left({{P_{\rm{A}}}-{T_{{\rm{hres}}}}} \right)}}{{{S_{{\rm{at}}}}-{T_{{\rm{hres}}}}}}\; \;\; \;\; {T_{{\rm{hres}}}} < {P_{\rm{A}}} < {S_{{\rm{at}}}}\\ {S_{\rm{f}}}\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{P_{\rm{A}}} \ge {T_{{\rm{hres}}}} \end{array} \right.$ |
肌肉的激活程度Act可以由一個線性分段函數來表述的，參見公式(1)。主要包含三個參數：電脈沖幅值的閾值Thres(mA)，電脈沖的飽和幅值Sat(mA)，比例因子Sf(Nm)。Sf是對肌肉強度以及舒適度的測量，它由肌肉在最佳肌肉長度最大刺激強度下肌肉產生的激活力矩來決定。
肌肉收縮部分主要由關節角度和關節速度兩個輸入，它們反映力矩-角度、力矩-速度的關系。事實上，肌肉力是由肌肉激活程度、長度、速度、最大收縮力共同計算得來的[10]，為了盡可能簡化模型，這一關系同樣可以由關節力矩與角度之間的關系來表示[11]。這里我們采用Hatze′s文章中的肌肉力與纖維長度之間的關系[12]，力矩-角度關系由下面的高斯公式給出：
| ${M_{\rm{a}}}=\exp \left\{ {-{{\left({\frac{{\varphi-{K_1}}}{{{K_2}}}} \right)}^2}} \right\}$ |
公式(2)中，φ是腕關節的角度，K1是Ma最大時所對應的關節角度，K2決定了關節的活動范圍。
力矩-角速度的關系是由Zajac文章[10]中的近似線性力-速度關系fv推導得到公式(3)：
| $\begin{array}{c} {f_{\rm{v}}}=1 + \frac{1}{{{v_{\rm{m}}}}}\frac{{{\rm{d}}l}}{{{\rm{d}}t}} \Rightarrow \\ {M_{\rm{v}}}=1 + \frac{1}{{{v_{\rm{m}}}}}\frac{{{\rm{d}}l}}{{{\rm{d}}\varphi }}\frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}}=1-\frac{1}{{{v_{\rm{m}}}}}r\varphi \left(t \right)=1-{K_3}\varphi \left(t \right) \end{array}$ |
其中，l是肌肉的長度，vm是肌肉的最大收縮速度，r是肌肉力矩，K3=r/vm。
肌肉收縮產生的力矩可以由公式(4)計算：
| ${M_{{\rm{act}}}}={A_{{\rm{ct}}}}{M_{\rm{v}}}{M_{\rm{a}}}$ |
Mv、Ma均為無量綱變量，Mact、Mct單位均為Nm。
模型中相關參數的值如表 1所示，具體可查閱相關文獻[8, 13]以及實際實驗測量。
表1
肌肉模型參數
Table1.
Parameters of the muscle model
1.2 具有耦合關系的兩自由度骨骼模型
為了使模型更貼近真實的腕關節骨骼模型，我們考慮了腕關節內外翻和上下切兩個方向上的耦合關系。Steven等[14]基于大量的實驗數據、數學計算以及相關文獻提出了一個在內外翻和上下切兩自由度上有耦合關系的腕關節動力學模型。他們把腕關節看成一個兩軸不相交的萬向節(圖 2)，腕關節上下切方向以x軸為旋轉軸，旋轉角度為β，腕關節內外翻方向以z軸為旋轉軸，旋轉角度為γ。
圖2
腕關節兩軸不相交的萬向節模型
Figure2.
Diagram of the actual wrist and universal joint model
腕關節動力學方程如下：
| $\left\{ \begin{array}{l} {M_\beta }={M_{inert, \beta }} + {M_{damp, \beta }} + {M_{stiff, \beta }}\\ {M_\gamma }={M_{inert, \gamma }} + {M_{damp, \gamma }} + {M_{stiff, \gamma }} + {M_{grav, \gamma }} \end{array} \right.$ |
公式(5)是一個非常復雜的模型，包含有很多的非線性成分，而且其中的很多成分對總的力矩貢獻是非常小的，對于仿真計算也是一個很大的負擔，因此在盡量不丟失真實特性的基礎上對系統做以下兩個近似處理：
近似一：忽略軸間距離和慣量交叉力矩。此時，實際力矩與近似力矩之間的誤差平均值僅有最大力矩的(0.52±0.58)%，在所有運動中的最大誤差不超過10%。
近似二：腕關節的運動范圍比較小。震顫是一種小幅度的有節律運動，一般在±15°之內，在這種情形下我們可以做近似：sin2γ≈0，cosγ≈1，cos2γ≈1，經過以上兩個近似之后總的累積誤差為(0.80±0.60)%，最大的誤差仍然不超過10%。
經過以上處理，基于震顫的小范圍運動的腕關節二階機械阻抗模型見公式(6)：
| $\begin{array}{c} \left[\begin{array}{l} {M_\beta }\\ {M_\gamma } \end{array} \right]=\left[{\begin{array}{*{20}{c}} {{I_\beta }} & 0\\ 0 & {{I_\gamma }} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{B_{\beta \beta }}} & {{B_{\beta \gamma }}}\\ {{B_{\beta \gamma }}} & {{B_{\gamma \gamma }}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right] + \\ \left[{\begin{array}{*{20}{c}} {{K_{\beta \beta }}} & {{K_{\beta \gamma }}}\\ {{K_{\beta \gamma }}} & {{K_{\gamma \gamma }}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right]-\left[\begin{array}{l} 0\\ mgr \end{array} \right] \end{array}$ |
選擇狀態向量
| $X={\left[{{x_1}, {x_2}, {x_3}, {x_4}} \right]^T}={\left[{\beta, \beta, \gamma, \gamma } \right]^T}$ |
結合公式(6)和公式(7)，得狀態空間方程如下：
| $\left\{ \begin{array}{l} X=AX + BU\\ Y=CX + DU \end{array} \right.$ |
其中：
| $\begin{array}{l} A=\left[{\begin{array}{*{20}{c}} 0 & 1 & 0 & 0\\ {\frac{{{K_{\beta \beta }}}}{{{I_\beta }}}} & {-\frac{{{B_{\beta \beta }}}}{{{I_\beta }}}} & {\frac{{{K_{\beta \gamma }}}}{{{I_\beta }}}} & {\frac{{{B_{\beta \gamma }}}}{{{I_\beta }}}}\\ 0 & 0 & 0 & 1\\ {-\frac{{{K_{\beta \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{B_{\beta \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{K_{\gamma \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{B_{\gamma \gamma }}}}{{{I_\gamma }}}} \end{array}} \right]\\ B=\left[{\begin{array}{*{20}{c}} 0 & 0\\ {\frac{1}{{{I_\beta }}}} & 0\\ 0 & 0\\ 0 & {\frac{1}{{{I_\gamma }}}} \end{array}} \right], \\ C=\left[{\begin{array}{*{20}{c}} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 \end{array}} \right], D=0. \end{array}$ |
表2
骨骼模型參數
Table2.
Parameters of the skeletal model
2 控制器設計
針對前面建立的肌骨模型本文提出了前饋加反饋的控制策略。其中LQR作為反饋控制器用來消除系統穩態誤差，肌肉模型的逆用作前饋控制器對反饋控制器的輸出進行補償，同時提高系統的實時響應性。
2.1 LQR反饋控制器
最優控制是指以最小的損耗實現系統的穩定運行，并且獲得一個最優的控制性能指標。本文中我們希望系統的期望輸入與系統實際輸出的誤差能夠到達一個極小值，同時又不希望對系統的穩定性和動態性能產生比較大的影響。在前面我們已經求出了含有兩個輸入兩個輸出的腕關節狀態方程，傳統的經典控制理論不能夠得到滿意的結果，那么在已知控制對象狀態空間方程的情況下采用LQR實現系統的最優控制將是理想選擇。LQR本質就是求解線性二次最優控制問題，其目標函數是系統狀態和控制輸入的二次型函數，目標就是選擇加權矩陣使得目標函數能夠取極小值。腕關節的骨骼模型是用狀態空間表示的線性定常系統，根據LQR原理選擇目標函數如下：
| $J=\int\limits_0^\infty {\left({{\boldsymbol{X}^T}\boldsymbol{QX} + {\boldsymbol{u}^T}\boldsymbol{Ru}} \right){\rm{d}}t} $ |
公式(9)中，Q為正定實定常對稱陣，狀態向量X的加權矩陣；R為正定實對稱常數矩陣，控制向量U的加權矩陣。根據極值原理，可以求出最優控制律對應的反饋增益矩陣K：
| $u\left(t \right)=-\boldsymbol{K}*x\left(t \right)$ |
此時目標函數J可以取到最小值。公式(10)中的矩陣K可以通過Matlab中的lqr函數進行求解，如下：
| $\boldsymbol{K}=lqr\left({A, B, \boldsymbol{Q}, \boldsymbol{R}} \right)$ |
不同的加權矩陣系統的響應也不相同，Q的值越大，系統到達穩態的時間越短，對于干擾的抵抗能力越強, R的值則與之相反[15]。仿真中我們會根據系統的實際響應來調整Q、R的值來獲得一個最優控制律下對應的反饋增益矩陣K。
2.2 基于模型逆的前饋控制器
我們已經建立了電刺激作用下的肌肉響應模型，由于其復雜的非線性，當我們使用LQR時非線性因素使得系統控制效果很差甚至是不穩定的，因此我們對電刺激作用下的肌肉模型進行求逆，并將其作為前饋控制器放在肌肉模型的前端。前饋控制器與肌肉模型串聯起來就近似于將肌肉模型線性化為單位系統，原來的非線性因素就被前饋控制器給補償掉了，從而可以使LQR控制器取得一個非常好的控制效果，較高的穩態控制精度。逆控制的基本思想就是要用一個來自控制器的信號去驅動對象，而該控制器的設計就是對控制對象本身進行求逆，逆控制是一種簡單的、魯棒的、且精確的控制方法[16]。
對電刺激作用下的肌肉響應模型進行求逆，根據公式(12)、(13)在Simulink中建模作為前饋控制器進行仿真控制，參數取值參見肌肉模型部分表 1。
| $\begin{array}{c} {A_{{\rm{ct}}}}={M_{{\rm{act}}}}*M_{\rm{v}}^{-1} \times M_{\rm{a}}^{-1}\\={M_{{\rm{act}}}} \times \exp \left\{ {{{\left({\frac{{\Phi-{K_1}}}{{{K_2}}}} \right)}^2}} \right\} \times \left\{ {1-{K_3} \times \Phi \left(t \right)} \right\} \end{array}$ |
| ${P_{\rm{A}}}=\left\{ \begin{array}{l} {S_{{\rm{at}}}}\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{A_{{\rm{ct}}}} \ge {S_{\rm{f}}}\\ \frac{{{A_{{\rm{ct}}}} \times \left({{S_{{\rm{at}}}}-{T_{{\rm{hres}}}}} \right)}}{{{S_{\rm{f}}}}} + {T_{{\rm{hres}}}}\; \;\; \;0 < {A_{{\rm{ct}}}}, {S_{\rm{f}}}\\ 0\;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{A_{{\rm{ct}}}} \le 0 \end{array} \right.$ |
3 仿真結果與分析
本文使用Matlab的Simulink進行仿真實驗，根據圖 1在Simulink中建立仿真系統框圖。參考輸入分別代表腕關節內外翻和上下切方向的運動，一般1.5 Hz以下；震顫的影響是通過分別施加在肌骨系統的腕關節內外翻和上下切方向的干擾力矩實現的，根據病理性震顫的特點，通常為3~8 Hz的高頻信號。系統的參考輸入相當于正常人腕關節的自主運動，系統的干擾力矩相當于患者肌骨系統中的病變信號，患者之所以震顫就是由于這些病變信號引起的。仿真目的就是要驗證所設計的系統可以抑制震顫，同時保證不影響自主運動。
依據前面已得到的比較精確的骨骼系統二階阻抗動力學模型，以LQR算法不斷選擇Q、R矩陣計算反饋增益矩陣，跟蹤系統的期望輸入，考察系統對于不同震顫干擾的抑制能力。仿真環境中，首先給定期望輸入為兩個方向的階躍信號：屈伸方向為20°，上下切方向為18°。干擾力矩在屈伸方向為0.2×sin(2×π×5×t), 在上下切方向為0.3×sin (2×π×5×t)。
通過試錯法(trial-and-error)，不斷的調整Q和R矩陣，我們發現在一定范圍內Q選的的越大，系統到達穩態的時間越短，穩態誤差越小，對于震顫的抑震效果越好，同樣減小R也可以得到類似的效果。經過不斷的仿真，我們選擇Q=[300, 0, 0, 0; 0, 0, 0, 0; 0, 0, 300, 0; 0, 0, 0, 0], R=[0.001, 0; 0, 0.001], 計算得到反饋增益矩陣K=[546.3044, 1.1624, 0.1703, 0.0002; 0.2081, 0.0002, 545.8754, 0.9456], 此時系統有比較理想的響應以及震顫抑制效果。圖 3所示是此時系統的輸出，經計算系統對于震顫的抑制效果達到94%，穩態誤差在0.1°以內。
圖3
在干擾力矩作用下系統對階躍信號的響應
Figure3.
System response to step signal with disturbance torques
我們又測試了系統在正弦曲線的參考輸入情況下的性能，也就是腕部兩個方向的自主運動為正弦曲線。我們設定在屈伸方向的期望輸入運動軌跡是5+10×sin(2×π×0.7×t)，上下切方向的期望輸入運動軌跡是5+3×sin(2×π×0.7×t+0.5)，干擾力矩不變。圖 4是仿真結果，5 s以前沒有加入FES控制器，在5 s時加入FES控制器，觀察系統輸出。由圖 4可以觀察到系統對于正弦信號的跟隨接近理想狀況，平均抑制效果達到98%以上。
圖4
參考輸入為正弦信號時的震顫抑制效果
Figure4.
Tremor suppression effect of the reference sinusoidal signal at input
為了使仿真接近真實情況，我們在試驗中采集了真實震顫患者的試驗數據(圖 5)作為系統的期望輸入以及震顫干擾力矩。患者性別男，74歲，患病時間20年，震顫頻率約為4.75 Hz，醫院診斷病情為帕金森。試驗要求患者將一個膠帶從自己左手邊平移到右手邊，整個過程中患者處于放松狀態，采集患者在這個任務中的腕關節角度信號用以仿真。
圖5
帕金森震顫患者實驗數據采集
Figure5.
Acquisition experiment data of a Parkinson patient with tremor
仿真結果如圖 6所示。經計算系統對于震顫的抑制程度為：內外翻方向95%，上下切方向92%。震顫抑制后對于原先的自主運動有微小的偏移，考慮電刺激抑制震顫的本質就是電刺激肌肉產生抵消震顫的力，相當于從外界對腕關節本身施加符合我們要求的力，因而這一微小的偏移完全可以接受。
圖6
震顫患者運動信息作為輸入時的震顫抑制效果
Figure6.
Tremor suppression effect with the movements of tremor patient as input
4 結論
本文針對腕關節首先建立了含有四塊肌肉(ECRL、FCU、FCR、ECU)、兩個自由度(內外翻和上下切)的肌骨模型，對肌肉激活、肌肉收縮、骨骼動力學分別進行了建模，考慮了FES作用以及兩自由度之間的耦合。以建立的肌骨模型為控制對象，針對肌骨模型的特點提出了基于模型逆的前饋控制與LQR最優控制作為反饋控制器的控制策略，并對控制器做了具體的設計與分析。最后在Matlab中建立了完整的基于FES的多自由度腕關節抑制系統，模擬了在震顫干擾力矩作用下系統對于階躍、正弦以及真實患者運動信息的響應。結果表明本文提出的控制策略對不同震顫類型有良好的抑制效果，驗證了本文控制策略的有效性以及可行性。在未來的工作中，將考慮自主運動未知的情況，即自主運動是由人體神經系統給出，仿真中需要完成自主運動與震顫信號的分離，實現對震顫信號的預測，為抑震系統走向實際應用提供支持。
引言
震顫是指身體某些運動功能區有節律、非自主的抖動，主要是由于神經系統對于肌肉的往復支配引起的肌肉周期性的收縮[1]。根據臨床表現，病理性震顫可以分為靜止性震顫、姿勢性震顫以及運動性震顫[2]。震顫減弱了患者處理日常生活的能力，大大降低了患者的生活質量[3]。目前藥物治療和手術治療是治療震顫的兩種主要方式，藥物治療雖暫時有所好轉但并不能阻止病情發展；手術治療也有局限，比如典型的深部腦刺激技術[4]，只是對某些病癥引起的病理性震顫有效，而且費用高、風險大、副作用明顯。據相關研究表明至少有25%的震顫患者無法通過藥物治療或者手術治療得到治愈[5]。正是由于傳統治療遇到的這些問題和瓶頸，越來越多的工程科研人員加入到病態震顫抑制這個領域，從而使輔助式抑震技術為病理性震顫的治療提供了新的選擇。
功能性電刺激(functional electrical stimulation, FES)是指用弱電流脈沖刺激人體的骨骼肌，肌肉因受到刺激收縮產生力以驅動關節運動，從而實現運動功能恢復和重建的目的，目前被大量應用于康復工程領域。加拿大的Prochazka團隊首次提出將FES作為一種重要的抑制病理性震顫的技術手段，他們利用FES異相刺激目標肌肉以期抑制肘關節、腕關節的病態顫抖[5]。法國的Poignet團隊提出了基于FES技術控制關節剛度來抑制顫抖的方法，設計了一種使用FES控制關節拮抗肌從而增強關節剛度的顫抖抑制系統[6]。筆者前期也開展了FES抑震研究，通過基于神經振蕩器的控制方式分別實現了對腕部和肘部的震顫抑制，取得了良好的效果[7]。經調研發現，目前的國內外相關工作都是針對腕關節或者肘關節單自由度的情況進行的，而且日常生活中的動作都是基于兩到三個甚至更多自由度完成的。本文的主要工作就是在前人的基礎上首先完成腕關節內外翻和上下切兩個自由度的肌骨模型建模，并完成基于FES的腕關節兩自由度病態震顫抑制的控制器設計與抑震系統的仿真。
本文提出的基于FES的抑震控制系統如圖 1所示。系統的期望輸入為自主運動，與反饋信號做差后經過FES控制器產生異相電刺激脈沖作用在兩自由度肌骨模型上從而達到震顫抑制的目的。
圖1
基于FES的腕關節抑震控制系統框圖
Figure1.
System diagram of the wrist joint tremor suppression based on FES
我們建立的肌骨模型包含兩對拮抗肌模型和二自由度的腕關節骨骼動力學，以肌骨模型為控制對象設計了FES控制器。我們提出了反饋加前饋的控制策略，其中反饋控制采用基于狀態空間的線性二次型調節器(linear quadratic regulator, LQR)，前饋控制器采用基于肌肉逆模型[8]，仿真結果表明該復合控制策略能夠有效地對病態震顫進行抑制，簡化了控制系統的復雜程度，降低了向實際抑震系統轉化的難度。
1 肌骨模型的建立
本文建立的腕部肌骨模型主要包括四塊肌肉，它們分別是橈側腕長伸肌(extensor carpi radialis longus, ECRL)、尺側腕伸肌(extensor carpi ulnaris, ECU)、橈側腕屈肌(flexor carpi radialis, FCR)、尺側腕屈肌(flexor carpi ulnaris, FCU)。FCU、FCR主要參與腕關節的內翻，ECRL、ECU主要參與腕關節的外翻，FCR、ECRL主要參與腕關節的上切，FCU、ECU主要參與腕關節的下切[9]。
本文建立的肌骨模型主要由三個部分組成：肌肉的激活，肌肉收縮，骨骼動力特性。第一部分肌肉實際產生的激活程度Act(Nm)由電刺激脈沖的幅值強度來計算；第二部分肌肉收縮特性實際反映關節角度和關節速度與實際關節力矩之間的關系；第三部分主要是關節的機械特性，當有力矩作用在骨骼上時表現為關節的運動。
1.1 FES作用下的肌肉模型
| ${A_{{\rm{ct}}}}=\left\{ \begin{array}{l} 0\;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; {P_{\rm{A}}} \le {T_{{\rm{hres}}}}\\ \frac{{{S_{\rm{f}}} \times \left({{P_{\rm{A}}}-{T_{{\rm{hres}}}}} \right)}}{{{S_{{\rm{at}}}}-{T_{{\rm{hres}}}}}}\; \;\; \;\; {T_{{\rm{hres}}}} < {P_{\rm{A}}} < {S_{{\rm{at}}}}\\ {S_{\rm{f}}}\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{P_{\rm{A}}} \ge {T_{{\rm{hres}}}} \end{array} \right.$ |
肌肉的激活程度Act可以由一個線性分段函數來表述的，參見公式(1)。主要包含三個參數：電脈沖幅值的閾值Thres(mA)，電脈沖的飽和幅值Sat(mA)，比例因子Sf(Nm)。Sf是對肌肉強度以及舒適度的測量，它由肌肉在最佳肌肉長度最大刺激強度下肌肉產生的激活力矩來決定。
肌肉收縮部分主要由關節角度和關節速度兩個輸入，它們反映力矩-角度、力矩-速度的關系。事實上，肌肉力是由肌肉激活程度、長度、速度、最大收縮力共同計算得來的[10]，為了盡可能簡化模型，這一關系同樣可以由關節力矩與角度之間的關系來表示[11]。這里我們采用Hatze′s文章中的肌肉力與纖維長度之間的關系[12]，力矩-角度關系由下面的高斯公式給出：
| ${M_{\rm{a}}}=\exp \left\{ {-{{\left({\frac{{\varphi-{K_1}}}{{{K_2}}}} \right)}^2}} \right\}$ |
公式(2)中，φ是腕關節的角度，K1是Ma最大時所對應的關節角度，K2決定了關節的活動范圍。
力矩-角速度的關系是由Zajac文章[10]中的近似線性力-速度關系fv推導得到公式(3)：
| $\begin{array}{c} {f_{\rm{v}}}=1 + \frac{1}{{{v_{\rm{m}}}}}\frac{{{\rm{d}}l}}{{{\rm{d}}t}} \Rightarrow \\ {M_{\rm{v}}}=1 + \frac{1}{{{v_{\rm{m}}}}}\frac{{{\rm{d}}l}}{{{\rm{d}}\varphi }}\frac{{{\rm{d}}\varphi }}{{{\rm{d}}t}}=1-\frac{1}{{{v_{\rm{m}}}}}r\varphi \left(t \right)=1-{K_3}\varphi \left(t \right) \end{array}$ |
其中，l是肌肉的長度，vm是肌肉的最大收縮速度，r是肌肉力矩，K3=r/vm。
肌肉收縮產生的力矩可以由公式(4)計算：
| ${M_{{\rm{act}}}}={A_{{\rm{ct}}}}{M_{\rm{v}}}{M_{\rm{a}}}$ |
Mv、Ma均為無量綱變量，Mact、Mct單位均為Nm。
模型中相關參數的值如表 1所示，具體可查閱相關文獻[8, 13]以及實際實驗測量。
表1
肌肉模型參數
Table1.
Parameters of the muscle model
1.2 具有耦合關系的兩自由度骨骼模型
為了使模型更貼近真實的腕關節骨骼模型，我們考慮了腕關節內外翻和上下切兩個方向上的耦合關系。Steven等[14]基于大量的實驗數據、數學計算以及相關文獻提出了一個在內外翻和上下切兩自由度上有耦合關系的腕關節動力學模型。他們把腕關節看成一個兩軸不相交的萬向節(圖 2)，腕關節上下切方向以x軸為旋轉軸，旋轉角度為β，腕關節內外翻方向以z軸為旋轉軸，旋轉角度為γ。
圖2
腕關節兩軸不相交的萬向節模型
Figure2.
Diagram of the actual wrist and universal joint model
腕關節動力學方程如下：
| $\left\{ \begin{array}{l} {M_\beta }={M_{inert, \beta }} + {M_{damp, \beta }} + {M_{stiff, \beta }}\\ {M_\gamma }={M_{inert, \gamma }} + {M_{damp, \gamma }} + {M_{stiff, \gamma }} + {M_{grav, \gamma }} \end{array} \right.$ |
公式(5)是一個非常復雜的模型，包含有很多的非線性成分，而且其中的很多成分對總的力矩貢獻是非常小的，對于仿真計算也是一個很大的負擔，因此在盡量不丟失真實特性的基礎上對系統做以下兩個近似處理：
近似一：忽略軸間距離和慣量交叉力矩。此時，實際力矩與近似力矩之間的誤差平均值僅有最大力矩的(0.52±0.58)%，在所有運動中的最大誤差不超過10%。
近似二：腕關節的運動范圍比較小。震顫是一種小幅度的有節律運動，一般在±15°之內，在這種情形下我們可以做近似：sin2γ≈0，cosγ≈1，cos2γ≈1，經過以上兩個近似之后總的累積誤差為(0.80±0.60)%，最大的誤差仍然不超過10%。
經過以上處理，基于震顫的小范圍運動的腕關節二階機械阻抗模型見公式(6)：
| $\begin{array}{c} \left[\begin{array}{l} {M_\beta }\\ {M_\gamma } \end{array} \right]=\left[{\begin{array}{*{20}{c}} {{I_\beta }} & 0\\ 0 & {{I_\gamma }} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{B_{\beta \beta }}} & {{B_{\beta \gamma }}}\\ {{B_{\beta \gamma }}} & {{B_{\gamma \gamma }}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right] + \\ \left[{\begin{array}{*{20}{c}} {{K_{\beta \beta }}} & {{K_{\beta \gamma }}}\\ {{K_{\beta \gamma }}} & {{K_{\gamma \gamma }}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} \beta \\ \gamma \end{array}} \right]-\left[\begin{array}{l} 0\\ mgr \end{array} \right] \end{array}$ |
選擇狀態向量
| $X={\left[{{x_1}, {x_2}, {x_3}, {x_4}} \right]^T}={\left[{\beta, \beta, \gamma, \gamma } \right]^T}$ |
結合公式(6)和公式(7)，得狀態空間方程如下：
| $\left\{ \begin{array}{l} X=AX + BU\\ Y=CX + DU \end{array} \right.$ |
其中：
| $\begin{array}{l} A=\left[{\begin{array}{*{20}{c}} 0 & 1 & 0 & 0\\ {\frac{{{K_{\beta \beta }}}}{{{I_\beta }}}} & {-\frac{{{B_{\beta \beta }}}}{{{I_\beta }}}} & {\frac{{{K_{\beta \gamma }}}}{{{I_\beta }}}} & {\frac{{{B_{\beta \gamma }}}}{{{I_\beta }}}}\\ 0 & 0 & 0 & 1\\ {-\frac{{{K_{\beta \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{B_{\beta \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{K_{\gamma \gamma }}}}{{{I_\gamma }}}} & {-\frac{{{B_{\gamma \gamma }}}}{{{I_\gamma }}}} \end{array}} \right]\\ B=\left[{\begin{array}{*{20}{c}} 0 & 0\\ {\frac{1}{{{I_\beta }}}} & 0\\ 0 & 0\\ 0 & {\frac{1}{{{I_\gamma }}}} \end{array}} \right], \\ C=\left[{\begin{array}{*{20}{c}} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 \end{array}} \right], D=0. \end{array}$ |
表2
骨骼模型參數
Table2.
Parameters of the skeletal model
2 控制器設計
針對前面建立的肌骨模型本文提出了前饋加反饋的控制策略。其中LQR作為反饋控制器用來消除系統穩態誤差，肌肉模型的逆用作前饋控制器對反饋控制器的輸出進行補償，同時提高系統的實時響應性。
2.1 LQR反饋控制器
最優控制是指以最小的損耗實現系統的穩定運行，并且獲得一個最優的控制性能指標。本文中我們希望系統的期望輸入與系統實際輸出的誤差能夠到達一個極小值，同時又不希望對系統的穩定性和動態性能產生比較大的影響。在前面我們已經求出了含有兩個輸入兩個輸出的腕關節狀態方程，傳統的經典控制理論不能夠得到滿意的結果，那么在已知控制對象狀態空間方程的情況下采用LQR實現系統的最優控制將是理想選擇。LQR本質就是求解線性二次最優控制問題，其目標函數是系統狀態和控制輸入的二次型函數，目標就是選擇加權矩陣使得目標函數能夠取極小值。腕關節的骨骼模型是用狀態空間表示的線性定常系統，根據LQR原理選擇目標函數如下：
| $J=\int\limits_0^\infty {\left({{\boldsymbol{X}^T}\boldsymbol{QX} + {\boldsymbol{u}^T}\boldsymbol{Ru}} \right){\rm{d}}t} $ |
公式(9)中，Q為正定實定常對稱陣，狀態向量X的加權矩陣；R為正定實對稱常數矩陣，控制向量U的加權矩陣。根據極值原理，可以求出最優控制律對應的反饋增益矩陣K：
| $u\left(t \right)=-\boldsymbol{K}*x\left(t \right)$ |
此時目標函數J可以取到最小值。公式(10)中的矩陣K可以通過Matlab中的lqr函數進行求解，如下：
| $\boldsymbol{K}=lqr\left({A, B, \boldsymbol{Q}, \boldsymbol{R}} \right)$ |
不同的加權矩陣系統的響應也不相同，Q的值越大，系統到達穩態的時間越短，對于干擾的抵抗能力越強, R的值則與之相反[15]。仿真中我們會根據系統的實際響應來調整Q、R的值來獲得一個最優控制律下對應的反饋增益矩陣K。
2.2 基于模型逆的前饋控制器
我們已經建立了電刺激作用下的肌肉響應模型，由于其復雜的非線性，當我們使用LQR時非線性因素使得系統控制效果很差甚至是不穩定的，因此我們對電刺激作用下的肌肉模型進行求逆，并將其作為前饋控制器放在肌肉模型的前端。前饋控制器與肌肉模型串聯起來就近似于將肌肉模型線性化為單位系統，原來的非線性因素就被前饋控制器給補償掉了，從而可以使LQR控制器取得一個非常好的控制效果，較高的穩態控制精度。逆控制的基本思想就是要用一個來自控制器的信號去驅動對象，而該控制器的設計就是對控制對象本身進行求逆，逆控制是一種簡單的、魯棒的、且精確的控制方法[16]。
對電刺激作用下的肌肉響應模型進行求逆，根據公式(12)、(13)在Simulink中建模作為前饋控制器進行仿真控制，參數取值參見肌肉模型部分表 1。
| $\begin{array}{c} {A_{{\rm{ct}}}}={M_{{\rm{act}}}}*M_{\rm{v}}^{-1} \times M_{\rm{a}}^{-1}\\={M_{{\rm{act}}}} \times \exp \left\{ {{{\left({\frac{{\Phi-{K_1}}}{{{K_2}}}} \right)}^2}} \right\} \times \left\{ {1-{K_3} \times \Phi \left(t \right)} \right\} \end{array}$ |
| ${P_{\rm{A}}}=\left\{ \begin{array}{l} {S_{{\rm{at}}}}\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{A_{{\rm{ct}}}} \ge {S_{\rm{f}}}\\ \frac{{{A_{{\rm{ct}}}} \times \left({{S_{{\rm{at}}}}-{T_{{\rm{hres}}}}} \right)}}{{{S_{\rm{f}}}}} + {T_{{\rm{hres}}}}\; \;\; \;0 < {A_{{\rm{ct}}}}, {S_{\rm{f}}}\\ 0\;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;\; \;{A_{{\rm{ct}}}} \le 0 \end{array} \right.$ |
3 仿真結果與分析
本文使用Matlab的Simulink進行仿真實驗，根據圖 1在Simulink中建立仿真系統框圖。參考輸入分別代表腕關節內外翻和上下切方向的運動，一般1.5 Hz以下；震顫的影響是通過分別施加在肌骨系統的腕關節內外翻和上下切方向的干擾力矩實現的，根據病理性震顫的特點，通常為3~8 Hz的高頻信號。系統的參考輸入相當于正常人腕關節的自主運動，系統的干擾力矩相當于患者肌骨系統中的病變信號，患者之所以震顫就是由于這些病變信號引起的。仿真目的就是要驗證所設計的系統可以抑制震顫，同時保證不影響自主運動。
依據前面已得到的比較精確的骨骼系統二階阻抗動力學模型，以LQR算法不斷選擇Q、R矩陣計算反饋增益矩陣，跟蹤系統的期望輸入，考察系統對于不同震顫干擾的抑制能力。仿真環境中，首先給定期望輸入為兩個方向的階躍信號：屈伸方向為20°，上下切方向為18°。干擾力矩在屈伸方向為0.2×sin(2×π×5×t), 在上下切方向為0.3×sin (2×π×5×t)。
通過試錯法(trial-and-error)，不斷的調整Q和R矩陣，我們發現在一定范圍內Q選的的越大，系統到達穩態的時間越短，穩態誤差越小，對于震顫的抑震效果越好，同樣減小R也可以得到類似的效果。經過不斷的仿真，我們選擇Q=[300, 0, 0, 0; 0, 0, 0, 0; 0, 0, 300, 0; 0, 0, 0, 0], R=[0.001, 0; 0, 0.001], 計算得到反饋增益矩陣K=[546.3044, 1.1624, 0.1703, 0.0002; 0.2081, 0.0002, 545.8754, 0.9456], 此時系統有比較理想的響應以及震顫抑制效果。圖 3所示是此時系統的輸出，經計算系統對于震顫的抑制效果達到94%，穩態誤差在0.1°以內。
圖3
在干擾力矩作用下系統對階躍信號的響應
Figure3.
System response to step signal with disturbance torques
我們又測試了系統在正弦曲線的參考輸入情況下的性能，也就是腕部兩個方向的自主運動為正弦曲線。我們設定在屈伸方向的期望輸入運動軌跡是5+10×sin(2×π×0.7×t)，上下切方向的期望輸入運動軌跡是5+3×sin(2×π×0.7×t+0.5)，干擾力矩不變。圖 4是仿真結果，5 s以前沒有加入FES控制器，在5 s時加入FES控制器，觀察系統輸出。由圖 4可以觀察到系統對于正弦信號的跟隨接近理想狀況，平均抑制效果達到98%以上。
圖4
參考輸入為正弦信號時的震顫抑制效果
Figure4.
Tremor suppression effect of the reference sinusoidal signal at input
為了使仿真接近真實情況，我們在試驗中采集了真實震顫患者的試驗數據(圖 5)作為系統的期望輸入以及震顫干擾力矩。患者性別男，74歲，患病時間20年，震顫頻率約為4.75 Hz，醫院診斷病情為帕金森。試驗要求患者將一個膠帶從自己左手邊平移到右手邊，整個過程中患者處于放松狀態，采集患者在這個任務中的腕關節角度信號用以仿真。
圖5
帕金森震顫患者實驗數據采集
Figure5.
Acquisition experiment data of a Parkinson patient with tremor
仿真結果如圖 6所示。經計算系統對于震顫的抑制程度為：內外翻方向95%，上下切方向92%。震顫抑制后對于原先的自主運動有微小的偏移，考慮電刺激抑制震顫的本質就是電刺激肌肉產生抵消震顫的力，相當于從外界對腕關節本身施加符合我們要求的力，因而這一微小的偏移完全可以接受。
圖6
震顫患者運動信息作為輸入時的震顫抑制效果
Figure6.
Tremor suppression effect with the movements of tremor patient as input
4 結論
本文針對腕關節首先建立了含有四塊肌肉(ECRL、FCU、FCR、ECU)、兩個自由度(內外翻和上下切)的肌骨模型，對肌肉激活、肌肉收縮、骨骼動力學分別進行了建模，考慮了FES作用以及兩自由度之間的耦合。以建立的肌骨模型為控制對象，針對肌骨模型的特點提出了基于模型逆的前饋控制與LQR最優控制作為反饋控制器的控制策略，并對控制器做了具體的設計與分析。最后在Matlab中建立了完整的基于FES的多自由度腕關節抑制系統，模擬了在震顫干擾力矩作用下系統對于階躍、正弦以及真實患者運動信息的響應。結果表明本文提出的控制策略對不同震顫類型有良好的抑制效果，驗證了本文控制策略的有效性以及可行性。在未來的工作中，將考慮自主運動未知的情況，即自主運動是由人體神經系統給出，仿真中需要完成自主運動與震顫信號的分離，實現對震顫信號的預測，為抑震系統走向實際應用提供支持。
