震顫是肢體一種非自主和重復性的擺動動作，在基于功能性電刺激 (FES) 的震顫抑制系統中可將其視為周期性擾動信號。因此，利用重復控制器調節施加在相應肌肉上的FES強度和時間，從而產生與震顫運動反相的肌肉力矩，是震顫抑制的一種可行手段。目前，大多數基于FES的重復控制系統都假設震顫是固定的單一頻率信號，但實際上震顫可能是多頻率信號，且震顫頻率也會隨時間變化。本文從頻率角度對意向性震顫患者數據進行分析，提出了內模切換的自適應重復控制器來抑制不同頻率的腕關節震顫。仿真和試驗結果表明，所提出的并聯多內模和串聯高階內模切換的自適應重復控制對多頻率和單頻率震顫信號可實現平均84.98%的抑制效果，比傳統單內模重復控制器和基于濾波的反饋控制器在震顫抑制性能上均有明顯提高。因此本文所提基于FES的自適應重復控制方法可以更為有效地解決意向性震顫患者腕部震顫問題，為患者后續運動功能障礙的康復提供有效的技術支撐。
引用本文： 張贊, 劉艷紅, 楚冰, 霍本巖, OWENS David Howard. 基于功能性電刺激的腕關節震顫抑制自適應重復控制研究. 生物醫學工程學雜志, 2023, 40(4): 663-675. doi: 10.7507/1001-5515.202202008 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
0 引言
意向性震顫也稱運動性震顫，是一種非自主、有節律、交替性的肌肉收縮和放松，通常發生在有目標意圖的上肢運動中，且當肢體接近目標物時震顫幅度顯著增大[1]。因此，震顫患者無法完成精準的上肢和手部運動，如飲食、穿衣、書寫等，嚴重影響其生活質量[2]。傳統的震顫抑制手段主要有藥物治療和外科手術兩種。藥物的長期服用會產生藥品依賴性，且會帶來副作用[3]。侵入性手術治療，如立體定向丘腦切開術[4]、立體定向蒼白球切開術[5]和深腦部電刺激[6]等，具有一定的手術風險和并發癥，且費用昂貴。
隨著機電一體化技術的發展，輔助式震顫抑制手段受到越來越多的重視[7-8]。其中，功能性電刺激（functional electrical stimulation，FES）因其具有便攜性、低成本、非創傷性等特點在康復醫療領域應用非常廣泛[9-11]。FES是一種人工電刺激，可以直接激活運動神經元使肌肉收縮，達到運動功能恢復和重建的目的。利用FES抑制震顫是通過設計閉環反饋控制系統產生適當的電脈沖信號刺激相關肌肉，使之產生與震顫運動相反的肌肉力矩，從而實現震顫幅度的減小。Prochazka等[12]最早提出了基于濾波的反饋控制器，并通過試驗對三種類型的單頻震顫分別達到73%、62%和38%的震顫抑制效果。Zhang等[13]使用基于FES的神經振蕩器結合比例積分微分（proportional integral derivative，PID）控制器，仿真實現90%的震顫抑制率。Taheri等[14]利用高通濾波器和反演法估計震顫力矩，在手臂關節模擬器上評估所提出的算法，取得了97.5%～99.2%震顫抑制。張偉等[15]提出了前饋控制與線性二次型控制相結合的復合控制策略，通過仿真實現了多自由度震顫幅值約94%的抑制效果。
由于震顫是肢體一種節律性、交替性擺動動作，是肌肉非自主、重復性的收縮和放松，故在震顫抑制控制系統中可將其視為周期性擾動信號。因此，基于內模原理的重復控制器比上述傳統基于濾波的反饋控制器在震顫抑制效果方面更顯著[16-18]。然而，現有基于重復控制的震顫抑制研究中通常假設震顫為固定的單一頻率信號，但實際患者的震顫信號可能具有多頻率特點[19-20]，且峰值頻率也不是固定的。并聯多內模重復控制對于多頻率周期性干擾信號具有很好的抑制性能[20-22]。Zhang等[20]針對多頻震顫信號設計并聯多內模重復控制方案，試驗結果表明可實現平均81.23%的震顫抑制。Kurniawan等[21]針對伺服電機多頻干擾的有效抑制設計了并聯多內模重復控制器，并通過仿真驗證了控制器的有效性。文獻[22]利用多內模重復控制方法實現了所有特征諧波的提取和跟蹤，提高了有源濾波器的動態響應特性。但以上研究并沒有考慮峰值頻率變化對控制性能的影響，由于患者震顫信號頻率的個體差異性和時變性等特點，本文將并聯多內模結構和串聯高階內模結構結合，提出了針對不同震顫頻率特點的內模切換自適應重復控制（adaptive repetitive control，A-RC）策略。與以往基于重復控制的震顫抑制相關研究相比，本文的主要創新之處如下：
（1）提出了基于哈默斯坦（Hammerstein）模型結構的手腕肌骨模型及其快速參數辨識方法，并期望通過搭建的手腕震顫抑制實驗平臺驗證模型的有效性。
（2）針對意向性患者震顫信號頻率特性的個體差異性和時變性，結合線性化控制設計了并聯多內模和串聯高階內模切換的A-RC控制器，分析了閉環系統穩定性，提出基于頻率修正逆的重復控制補償器構造方法。
（3）利用搭建的手腕震顫抑制實驗平臺，通過健康參與者的誘發性震顫抑制試驗和意向性患者的手腕震顫抑制試驗，以期驗證所提A-RC方法的可行性和有效性。
1 方法
1.1 手腕肌骨模型
為了設計刺激強度和刺激時間精準調節的A-RC控制器，首先需要建立手腕肌骨模型。手腕肌骨模型由肌肉模型和骨骼模型組成。非線性Hammerstein模型由于其結構簡單，便于震顫抑制控制器的設計，被廣泛應用于基于FES的肌肉建模中[23]。Hammerstein模型由靜態非線性招募特性和線性動態模型組成，前者表現肌肉的激活特性，后者表現肌肉對電刺激信號的收縮響應。而骨骼模型的建立通常描述為線性剛體動力學（rigid body dynamics，RBD）模型。
本文主要研究手腕屈伸運動時的震顫抑制。腕部屈伸運動主要由橈側腕伸肌（extensor carpi radialis， ECR）和橈側腕屈肌（flexor carpi radialis，FCR）驅動橈骨實現。因而，根據Hammerstein模型結構，構造了FCR和ECR的非線性肌肉招募特性fFCR(uFCR)、fECR(uECR)和線性肌肉收縮動態特性GFCR(z)、GECR(z) ，以及骨骼線性剛體動力學模型GRBD(z)。電刺激手腕肌骨模型如圖1所示，其中uFCR(k)和uECR(k)是電刺激輸入信號，y(k)為手腕屈伸運動角度輸出信號。wFCR(k)和wECR(k)為招募特性的輸出，
為作用在骨骼系統上的肌肉總力矩。d(k)可以看作是由于震顫引起的干擾信號。
圖1
電刺激手腕肌骨模型
Figure1.
Electrically stimulated wrist musculoskeletal model
為簡化模型結構和參數辨識過程，將uFCR(k)、uECR(k)分別等效為如式(1)、式(2)所示的單一輸入
：
 |
 |
其中，uFCR,ZO和uECR,ZO是FCR和ECR肌肉招募特性的死區閾值，根據試驗設置為50 μs。由于典型的招募特性包括初始死區、單調遞增區和飽和區，進而等效非線性招募特性輸出w(k)與輸入u(k)可以通過如式(3)所示的靜態映射關系表示：
 |
其中，umax是最大電刺激輸入脈寬，為避免過強電刺激輸入造成肌肉收縮時的不適感，本文最大脈寬設置為300 μs [24]。
和 
分別為FCR和ECR肌肉的非線性招募特性，
和 
分別為 
和 
肌肉的招募特性死區閾值。
為等效的橈側肌肉非線性招募特性。化簡后的非線性招募特性f (u(k))在工作區間 
是連續且遞增的，如式(4)所示：
 |
其中，
以及 
分別為FCR和ECR 肌肉需要識別的非線性標量參數；e為指數函數。
由于相似的肌群具有相似的線性收縮特性[25]，因而FCR和ECR動態特性近似為 
，可用GL(z)表示。人體骨骼動力學可以建模為剛體動力學GRBD(z)，其中阻尼和彈性函數都可以看作線性函數[26]。故，等效的線性動態肌骨模型 
如式(5)所示：
 |
其中，
和 
分別為等效線性肌骨模型數學表達式的分子和分母，
和
分別是
的階數和待辨識的參數，
和
分別是
的階數和待辨識的參數。
1.2 控制器設計
由于腕關節肌骨模型具有固有的非線性，首先設計線性化控制器 
，在此基礎上設計A-RC控制器C(z)以抑制不同類型的腕部震顫，提高手腕對參考軌跡的跟蹤性能。基于FES的閉環震顫抑制系統控制框圖如圖2所示，其中r(k)為指定的跟蹤參考位置，e(k)為手腕實際運動位置與參考位置的偏差，
為控制器C(z)的輸出。
圖2
基于FES的閉環震顫抑制控制框圖
Figure2.
FES based closed-loop tremor suppression control diagram
1.2.1 線性化控制器設計
由于電刺激肌骨模型具有復雜的非線性因素，且非線性主要表現在靜態招募特性，因而采用線性化控制器補償肌骨模型中的招募非線性，線性化控制器 
如式(6)所示：
 |
1.2.2 A-RC控制器設計
重復控制是基于內模原理的一種控制方法。內模原理指出，如果在穩定的閉環系統中包含外部周期性信號的數學模型，系統就能夠實現對周期性信號的無靜差跟蹤或抑制。本文根據患者震顫數據的分析，設計了內模切換的A-RC控制器。經過線性化之后，等效的閉環控制系統框圖如圖3所示，H(z)是提高閉環系統穩定性和動態性能的重復控制補償器。自適應調節器根據傅里葉變換得到的震顫數據頻域特性調節內模結構，實現并聯多內模或串聯高階內模的自適應選擇，即根據震顫信號頻域分析結果，調節內模結構的數學模型，以實現對不同震顫患者最佳的震顫抑制效果。
圖3
等效的閉環震顫抑制系統框圖
Figure3.
Equivalent closed-loop tremor suppression system diagram
其中，A-RC控制器C(z)如式(7)所示：
 |
其中，補償器H(z)如式(8)所示：
 |
H(z)包含n個零點，其位置由系數c1, c2, ···, cn確定，并在原點引入n ? m個極點。自適應內模結構 
由m條并聯支路構成，Wq為第q條內模支路增益，q = 1, 2, ···, m。其中，第q條高階內模支路 
如式(9)所示：
 |
其中，
為一個周期的采樣次數，
為震顫信號的周期，
為采樣周期，
（p = 1, 2, ···, n）是第p階支路增益。
1.2.3 穩定性分析
定理1：在A-RC控制器的作用下，選擇控制增益 
、
和補償器 
使得閉環控制系統特征多項式的根均在單位圓內，如式(10)所示：
 |
那么，如圖3所示的閉環控制系統是漸近穩定的，且能夠完全抑制 
周期的擾動。尤其是，選擇 
，則當 
時，
需要滿足如式（11）所示：
 |
證明：根據圖(3)可得誤差 
、參考輸入 
和擾動 
之間的關系如式(12)所示：
 |
將式(7)帶入式(12)，則閉環系統的特征多項式如式(13)所示：
 |
當特征多項式的根都在單位圓內時，則閉環反饋控制系統是漸近穩定的。特別地，當 
時，式(13)可寫為如式(14)、式(15)所示：
 |
 |
為了有效抑制震顫，需要選擇重復控制器的內模結構 
在震顫信號基波和諧波處增益無窮大，即 
，令 
從而得到離散系統的頻率響應，其中，
表示虛數單位，
為角頻率，那么如式(16)所示：
 |
故當選擇 
，
時，有 
，這時第 
條內模支路增益 
選擇滿足 
。
證明完畢。
1.2.4 補償器H(z)設計
在實際中，由于手腕肌骨模型參數的不確定性和模型單位圓外可能存在零點等因素，補償器H(z)無法通過 
得到。因此，本文采用頻率修正逆重復控制方法來獲得補償器H(z)的階數和參數[20, 27]。即適當選擇式（8）中的參數c=[c1, c2,…, cn]T來最小化以下成本函數J，從而實現 
的頻率響應。成本函數J如式(17)所示：
 |
其中，
為一個周期的采樣次數，i表示虛數單位，
是從零到奈奎斯特的一組離散頻率，
為采樣周期，*表示復共軛。
2 仿真與試驗驗證
在試驗驗證環節，招募4名健康參與者和4名意向性震顫患者。健康參與者編號為No.1～No.4，3男1女，平均年齡為（28.25 ± 8.54）歲。震顫患者來自鄭州大學第五附屬醫院康復科，編號為No.5～No.8，3男1女，平均年齡為（43.25 ± 10.90）歲。試驗中8位受試者均簽署書面知情同意書，且已通過鄭州大學倫理審查，倫理審批編號為No.ZZURIB2019-004。基于搭建的手腕震顫抑制實驗平臺，每位受試者首先進行手腕肌骨模型參數辨識試驗，然后根據辨識的模型參數設計相應控制器，進而開展手腕震顫數據采集和震顫抑制試驗。
2.1 實驗平臺
基于FES的腕部震顫抑制實驗平臺如圖4所示。震顫抑制實驗平臺由輪椅改裝而成，其中手臂支架是固定的，前端的U型手掌支撐夾板運動范圍被限制在± 80°之內，這符合手腕屈伸的運動范圍。試驗時受試者舒適、放松地坐在實驗平臺的輪椅上，右前臂放于手臂支架，右手放在U型手掌支撐夾板中間。受試者的前臂用綁帶固定，以防止產生除手腕屈伸平面外其他任何自由度的運動。支架外側設計扇形角度標尺，方便受試者跟蹤指定的參考位置。實時的手腕屈伸運動角度數據由編碼器（E6B2-CWZ6C，深圳布瑞特科技有限公司，中國）測量，其中編碼器的軸與U型手掌支撐夾板的轉軸通過直徑比為1∶1的齒輪連接。仿真系統平臺dSPACE（Microlabbox 1202, dSPACE Inc.，德國）實現手腕角度數據和控制信號的實時傳輸，為控制算法的開發和驗證提供支持。最終的控制信號由電刺激儀（4 CHANNEL STIMULATOR, Odstock Medical Ltd.，英國）輸出。本文電刺激信號的頻率、電流幅值和最大脈寬設置分別為40 Hz、20 mA和300 μs，并且可以根據受試者對電刺激信號強弱的不同反應，利用電刺激儀上每個通道的電位器調節信號的大小。
圖4
腕部震顫抑制實驗平臺
Figure4.
Experimental platform for wrist tremor suppression
在臨床試驗之前，對健康參與者進行控制方案的可行性測試是必要的步驟。因此，在實驗平臺上嵌入直流電機（VDM08SGN24-60-1800/JB5G12T，北京微特微電機制造有限公司，中國），利用電機來模擬震顫信號。其原理是U型手掌支撐夾板與電機同軸，當電機產生輸出轉矩時，可直接帶動U型手掌支撐夾板動作。因此，通過軟件設置震顫頻率和幅度，電機驅動器（ESCON50/5，maxon motor Inc.，匈牙利）驅動電機運轉即可輸出相應的轉速和角度。這樣，健康參與者的手腕會產生不自主的屈伸運動，即得到誘發性震顫信號。
2.2 手腕肌骨模型參數辨識
手腕肌骨模型的參數辨識基于輸入/輸出數據集［u, y］，利用迭代最小二乘法獲得，整個辨識過程大約2 min左右。輸入信號u(k)為0.2～4.0 Hz、脈沖寬度為300 μs的不同頻率正弦信號。y(k)為編碼器測量的手腕屈伸運動的角度數據。在參數辨識過程中，受試者手腕運動是由輸入的電刺激信號驅動，而不是自主運動產生。為了計算模型的擬合度，定義最佳擬合比 
如式(18)所示：
 |
其中，y為實測輸出，
為模型模擬輸出。y的均值用 
表示。一般來說 
值越高，參數辨識的精度越高，但模型階數也越高。在本文中，
和 
的值［即 
中 
和 
的階數］從4～10進行選取。以受試者No.1為例，其最佳擬合率如圖5所示。相對低階模型的最佳擬合率為62.06%（
和 
），相對高階模型的最佳擬合率最高為67.94%（
和 
）。當 
和 
在4～10之間取值時，最佳擬合率變化并不明顯。因此，為了減少辨識計算量、簡化控制器設計，本文選取了 
。
圖5
受試者No.1手腕模型最佳擬合率
Figure5.
The best fit rate of participant No.1’s wrist model
2.3 震顫信號分析
利用搭建的實驗平臺采集意向性震顫患者No.5～No.8手腕屈伸運動時伴隨的震顫數據，如表1所示。
表1
患者No.5～No.8數據
Table1.
Data of Patients No.5～No.8
具體的采集方法為：讓意向性震顫患者右手通過屈伸手腕自主地跟蹤一個指定的參考軌跡，指定的參考軌跡類似于一些簡單的日常任務，例如，屈伸手腕以達到一個特定的角度。患者需要通過自主屈伸手腕跟蹤參考軌跡，產生的運動跟蹤軌跡為黑色虛線。由于震顫的存在，實際手腕運動的跟蹤軌跡會偏離指定的參考軌跡。如圖6所示，通過隨機選擇的患者No.5和No.7的跟蹤軌跡數據時頻圖可以看到，自主運動的頻率范圍低于1 Hz，且在20 s和40 s手腕有自主屈伸動作時，功率譜密度增強。意向性震顫頻率范圍主要在2～5 Hz之間[1, 28]，本文采集的意向性震顫頻率在2～3 Hz之間，且震顫頻率隨時間會有稍許的變化，如患者No.5震顫信號時頻圖40 s之后，震顫頻率明顯降低。
圖6
震顫患者No.5和No.7震顫信號頻率分析
Figure6.
Frequency analysis of tremor signals in patient No.5 and No.7
通過對患者No.5和No.7跟蹤軌跡信號分析可得，患者No.5為單頻率震顫信號，而患者No.7為多頻率震顫信號。對于單頻率或多頻率的震顫信號，通過自適應調節器選擇切換不同的重復控制器內模結構，即通過傅里葉變換得到最大振幅對應的震顫頻率（frequency corresponding to peak amplitude，FPA），如患者No.5的FPA為2.2 Hz，患者No.7的FPA為2.3 Hz。若其余頻率振幅均小于FPA的80% （患者No.5沒有超過FPV的80%的震顫頻率，故歸類為單頻率震顫信號），則選擇單一支路的高階內模結構（圖3中q = 1）；若震顫頻率幅值中有高于80%的FPA頻率（患者No.7在2.5 Hz處的震顫幅值超過FPV的80%，故歸類為多頻率震顫信號），則選擇并聯內模結構（圖3中q > 1），因為震顫頻率的振幅越大，對患者的影響越大。
2.4 仿真測試
在本小節中，使用受試者No.1的手腕模型驗證提出的A-RC方法對周期性震顫抑制的有效性，并將并聯多內模結構和串聯高階內模結構切換的A-RC方法與單周期重復控制方法（single periodic repetitive control，SP-RC）[16]和傳統的基于高通濾波器的比例積分控制方法（proportional integral control algorithm based on high pass filter，PI-HPF）[12]進行比較。根據震顫信號的單頻率或多頻率特點，將仿真測試分為以下兩個環節。
2.4.1 仿真環節1結果與分析
仿真測試環節1中，震顫信號d(k)采用頻率為2 Hz、幅值為1的正弦信號。采樣周期TS為0.005 s，因此A-RC控制器q = 1，選擇單一支路的串聯高階內模結構。W1 = 1，內模延時周期N1 = 100，p = 2，K1 = 2，
。對于單周期重復控制器，圖3中q = 1，W1 = 1，N1 = 100，p = 1且K1 = 1。重復控制器補償器H(z)由頻率修正逆算法得到。采用臨界比例度法確定了PI-HPF調節器的參數（P = 36.82，I = 1.4），高通巴特沃斯濾波器階數為6，截止頻率為1.2 Hz。設定參考輸入r(k)為零，輸出y(k)在干擾d(k)的作用下跟蹤參考信號。
根據采集的患者震顫數據頻譜分析可得，震顫信號的峰值頻率隨時間會產生稍許變化，故在仿真環節進行蒙特卡洛實驗，以驗證所提控制器對于頻率變化的周期性震顫信號的抑制效果。令震顫信號d(k)在峰值頻率±5%的范圍內隨機波動（即1.9～2.1 Hz隨機變化），控制器參數不變，進行100次仿真實驗，如圖7所示。圖7左側圖為仿真環節1中不同控制器下一次仿真測試的跟蹤誤差結果，PI-HPF方法的控制輸入相對于兩個重復控制器（A-RC和SP-RC）小一些，但是震顫抑制效果最差，因而驗證了對于周期性擾動信號重復控制器的優越性。對比A-RC和SP-RC兩種重復控制方法，在震顫峰值頻率有所變化的情況下（設定2 Hz，實際2.05 Hz），本文所提的A-RC控制器跟蹤誤差信號最小，震顫抑制效果最佳。
圖7
不同控制器下的跟蹤誤差和控制量
Figure7.
Tracking error and control input under different controllers
通過以下兩個指標來量化各控制方法在震顫抑制方面的性能。
（1）均方根誤差（root mean square error，RMSE）的計算公式如式(19)所示：
 |
其中，
為輸出信號，
為參考信號，n為測試的總樣本數。
（2）震顫抑制率（tremor suppression rate，TSR）如式(20)所示：
 |
其中，
是輸入信號 
和輸出信號 
的偏差。
為震顫信號對跟蹤角度的影響。
是輸入信號 
和 
之間的偏差。
與 
的比值為震顫抑制性能。
仿真環節1蒙特卡洛實驗的統計結果如圖8左側圖所示。PI-HPF方法受震顫信號頻率變化影響最小，其次是本文所提的A-RC方法，而SP-RC方法最差。但PI-HPF方法的RMSE均值（0.357）卻遠遠大于A-RC和SP-RC方法。SP-RC方法的RMSE均值為0.095，本文所提A-RC方法RMSE均值僅為0.043。RMSE 均值越小，意味著跟蹤誤差越小。在TSR的統計指標中，A-RC方法、SP-RC 方法和PI-HPF方法TSR均值分別為93.90%，85.26%和49.58%，TSR指標越大，震顫抑制效果越好。因此，本文所提的A-RC方法對于單頻率變化的周期性震顫信號具有最佳的抑制效果。
圖8
蒙特卡洛實驗統計結果
Figure8.
Statistical results of Monte Carlo trials
2.4.2 仿真環節2結果與分析
在仿真測試環節2中，震顫信號d(k)采用頻率為2 Hz和2.5 Hz，幅值為1和0.8的疊加正弦信號。A-RC控制器切換到并聯多內模結構，q = 2，W1 = W2 = 0.5，內模延時周期分別為N1 = 100和N2 = 80，K1 = K2 = 0.5，其余參數與仿真環節1一致。對于SP-RC控制器和PI-HPF控制器參數均不變。
在不同控制器下系統跟蹤誤差如圖7仿真環節2所示，控制量和仿真環節1的結果類似，PI-HPF方法的控制量最小，A-RC和SP-RC方法的控制量相差不大。但是SP-RC方法的跟蹤誤差和仿真環節1相比有明顯增加，因為單內模結構只能抑制內模中已有的周期性信號，而對于2.5 Hz的信號無法抑制。相比較而言，本文所提的A-RC方法可以實現多頻率擾動信號的完全抑制。利用蒙特卡洛實驗討論在震顫頻率變化情況下，控制器對于多頻率震顫信號的抑制性能。震顫信號
在峰值頻率±5%的范圍內波動，100次仿真測試的統計結果如圖8仿真環節2所示。從統計結果分析來看，PI-HPF和SP-RC方法的RMSE均值接近，分別為0.408和0.394，這說明在震顫頻率變化的情況下，這兩個方法的穩態跟蹤誤差大，TSR指標僅為48.36%和50.21%，無法實現多頻率震顫信號的抑制。而本文所提的A-RC控制器，利用其并聯多內模結構，在頻率變化的情況下，穩態誤差小，TSR均值可達90.27%，即實現對多頻率震顫信號的良好抑制。
2.5 試驗驗證
試驗驗證也分為兩個階段，第一階段在健康參與者身上開展誘發性震顫抑制試驗，用以評價所提控制方法的可行性。第二階段招募意向性震顫患者進行震顫抑制驗證試驗，進一步驗證所提方法的有效性。
2.5.1 試驗階段1結果與分析
在試驗階段1中，通過震顫抑制實驗平臺產生誘發性震顫。為了防止肌肉疲勞，每個受試者的電刺激時長不超過20 min。4名健康參與者需要完成以下三個測試：
測試1（test 1，T1）：自主地屈伸手腕跟蹤參考軌跡，此測試無誘發性震顫和電刺激輸入。
測試2（test 2，T2）：自主地屈伸手腕跟蹤參考軌跡，但此測試加入誘發性震顫信號。
測試3（test 3，T3）：在具有誘發性震顫和功能性電刺激輸入的情況下，自主屈伸手腕跟蹤指定參考軌跡（每位受試者在同一控制器下進行3次T3試驗）。
為了更好地對比不同控制器下T3的試驗結果，將誘發震顫設置為多頻震顫，頻率為2 Hz和2.5 Hz。受試者No.1在不同控制器作用下其中一次試驗的手腕運動軌跡和電刺激控制輸入信號如圖9所示。
圖9
受試者No.1在不同控制器下的手腕角度位置和電刺激信號
Figure9.
Wrist angular position and electrical stimulation level of participant No.1 under different controllers
如圖9中的T3手腕跟蹤軌跡所示，3種不同的控制器均可以不同程度地實現誘發性震顫的抑制，但本文所提的A-RC方法的T3運動軌跡相對參考軌跡的偏離程度最小，具有最佳的震顫抑制效果。其次是SP-RC控制器，PI-HPF控制器的震顫抑制性能最差。對于電刺激控制輸入，PI-HPF控制器產生的電刺激控制輸入信號相對重復控制器要小一些，但震顫抑制性能差。A-RC控制器和SP-RC控制器的電刺激輸入信號相差不大，結合T3的手腕運動軌跡可得本文所提A-RC方法具有震顫抑制的性能優勢。如圖10所示，受試者No.1在不同控制器下T3的跟蹤軌跡時頻圖也可以得到相同的結論。在不同控制器作用下，多頻率震顫信號的功率譜密度在震顫頻率處均有不同程度的下降。其中，本文所提A-RC方法的功率譜密度下降最為明顯，其次是SP-RC方法，而PI-HPF方法對于周期性震顫信號的抑制作用最差。
圖10
參與者No.1在不同控制器下跟蹤軌跡的時頻圖
Figure10.
The time-frequency diagram of tracking trajectory of participant No.1 under different controllers
試驗環節采用與仿真環節相同的性能指標。如式（19）所示，RMSE值用于量化T3（
）與T1（
）的之間的偏離程度。如式（20）所示，
是T3（
）與T1（
）之間的偏差。
是T2（
）與T1（
）之間的偏差。其余三位健康參與者也得到了相似的試驗結果，統計結果如圖11所示。
圖11
四位健康參與者試驗統計結果
Figure11.
Statistical results of experiment with four participants
四位健康參與者RMSE均值A-RC方法最小為1.34，而SP-RC方法為3.12，PI-HPF方法為3.58，這一試驗結果與仿真結果類似。通過TSR值可以更為直觀地判斷每種控制方法對周期性震顫信號的抑制情況，如圖11所示四位健康參與者在本文所提的A-RC方法下震顫抑制性能顯著提升，四位健康參與者TSR均值為87.59%，比SP-RC方法和PI-HPF方法均有明顯提升。
2.5.2 試驗階段2結果與分析
本階段招募意向性震顫患者開展震顫抑制試驗，定義以下三個測試試驗：
測試4（test 4，T4）：將指定參考軌跡提供給患者，患者不需任何動作。
測試5（test 5，T5）：患者自主地屈伸手腕跟蹤參考軌跡，此測試無FES信號。
測試6（test 6，T6）：在FES作用下患者自主屈伸手腕跟蹤參考軌跡（同一控制器下患者需進行3次T6試驗）。
隨機選擇患者No.5和No.7的T6跟蹤軌跡時頻圖如圖12所示。從時頻圖可以看到，本文所提的A-RC方法，對于單頻率震顫和多頻率震顫均有很好的抑制效果，在震顫頻率處的功率譜密度均有明顯下降，比PI-HPF方法在震顫抑制性能方面有明顯提升。而SP-RC方法對于周期性震顫信號也有一定的抑制作用，但效果相對不佳。患者試驗環節也采用相同的性能指標。如式（19）所示，RMSE值用于量化T6（
）與T4（
）的之間的偏離程度。如式(20)所示，
是T6（
）與T4（
）之間的偏差。
是T5（
）與T4（
）之間的偏差。患者No.5和No.7的量化試驗統計結果如圖13所示。
圖12
患者No.5和No.7在不同控制器下T6跟蹤軌跡的時頻圖
Figure12.
The time-frequency diagram of T6 tracking trajectory of patients No.5 and No.7 under different controllers
圖13
意向性震顫患者試驗統計數據
Figure13.
Statistical experimental results of intention tremor patients
如圖13所示，患者震顫抑制試驗結果與第1階段健康參與者的結果類似。A-RC方法對于單頻率和多頻率，且頻率存在微變的震顫均有很好的抑制作用。A-RC方法實現平均82.53%的單頻震顫和82.21%多頻震顫抑制，而SP-RC方法和PI-HPF方法對于這兩種類型震顫的抑制率僅分別為71.22%、58.31%和67.88%、52.33%。
3 結論
本文提出一種基于FES的A-RC方法來抑制腕關節意向性震顫。首先建立了非線性Hammerstein結構的手腕肌骨模型，提出了自動、快速的參數辨識算法，在搭建的手腕震顫抑制實驗平臺完成了肌骨模型的參數辨識。其次，根據辨識的每位受試者肌骨模型設計相應的線性化控制與反饋A-RC相結合的控制器，并根據患者震顫頻率特點自動選擇重復控制內模結構，分析了閉環系統穩定性，并利用頻率修正逆算法實現了重復控制補償器的構造。最后，通過健康參與者誘發性震顫抑制試驗和意向性患者震顫抑制試驗驗證了本文提出的A-RC策略，其平均震顫抑制率可達84.98%，比傳統單內模重復控制和基于濾波的反饋控制均有明顯提高。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：張贊主要負責數據收集與分析、算法程序設計以及論文編寫；霍本巖主要負責試驗素材收集；劉艷紅、楚冰和Owens David Howard主要負責試驗指導、數據分析指導以及論文審閱修訂。
倫理聲明：本研究通過了鄭州大學倫理委員會的審批（批文編號：No.ZZURIB2019-004）。
0 引言
意向性震顫也稱運動性震顫，是一種非自主、有節律、交替性的肌肉收縮和放松，通常發生在有目標意圖的上肢運動中，且當肢體接近目標物時震顫幅度顯著增大[1]。因此，震顫患者無法完成精準的上肢和手部運動，如飲食、穿衣、書寫等，嚴重影響其生活質量[2]。傳統的震顫抑制手段主要有藥物治療和外科手術兩種。藥物的長期服用會產生藥品依賴性，且會帶來副作用[3]。侵入性手術治療，如立體定向丘腦切開術[4]、立體定向蒼白球切開術[5]和深腦部電刺激[6]等，具有一定的手術風險和并發癥，且費用昂貴。
隨著機電一體化技術的發展，輔助式震顫抑制手段受到越來越多的重視[7-8]。其中，功能性電刺激（functional electrical stimulation，FES）因其具有便攜性、低成本、非創傷性等特點在康復醫療領域應用非常廣泛[9-11]。FES是一種人工電刺激，可以直接激活運動神經元使肌肉收縮，達到運動功能恢復和重建的目的。利用FES抑制震顫是通過設計閉環反饋控制系統產生適當的電脈沖信號刺激相關肌肉，使之產生與震顫運動相反的肌肉力矩，從而實現震顫幅度的減小。Prochazka等[12]最早提出了基于濾波的反饋控制器，并通過試驗對三種類型的單頻震顫分別達到73%、62%和38%的震顫抑制效果。Zhang等[13]使用基于FES的神經振蕩器結合比例積分微分（proportional integral derivative，PID）控制器，仿真實現90%的震顫抑制率。Taheri等[14]利用高通濾波器和反演法估計震顫力矩，在手臂關節模擬器上評估所提出的算法，取得了97.5%～99.2%震顫抑制。張偉等[15]提出了前饋控制與線性二次型控制相結合的復合控制策略，通過仿真實現了多自由度震顫幅值約94%的抑制效果。
由于震顫是肢體一種節律性、交替性擺動動作，是肌肉非自主、重復性的收縮和放松，故在震顫抑制控制系統中可將其視為周期性擾動信號。因此，基于內模原理的重復控制器比上述傳統基于濾波的反饋控制器在震顫抑制效果方面更顯著[16-18]。然而，現有基于重復控制的震顫抑制研究中通常假設震顫為固定的單一頻率信號，但實際患者的震顫信號可能具有多頻率特點[19-20]，且峰值頻率也不是固定的。并聯多內模重復控制對于多頻率周期性干擾信號具有很好的抑制性能[20-22]。Zhang等[20]針對多頻震顫信號設計并聯多內模重復控制方案，試驗結果表明可實現平均81.23%的震顫抑制。Kurniawan等[21]針對伺服電機多頻干擾的有效抑制設計了并聯多內模重復控制器，并通過仿真驗證了控制器的有效性。文獻[22]利用多內模重復控制方法實現了所有特征諧波的提取和跟蹤，提高了有源濾波器的動態響應特性。但以上研究并沒有考慮峰值頻率變化對控制性能的影響，由于患者震顫信號頻率的個體差異性和時變性等特點，本文將并聯多內模結構和串聯高階內模結構結合，提出了針對不同震顫頻率特點的內模切換自適應重復控制（adaptive repetitive control，A-RC）策略。與以往基于重復控制的震顫抑制相關研究相比，本文的主要創新之處如下：
（1）提出了基于哈默斯坦（Hammerstein）模型結構的手腕肌骨模型及其快速參數辨識方法，并期望通過搭建的手腕震顫抑制實驗平臺驗證模型的有效性。
（2）針對意向性患者震顫信號頻率特性的個體差異性和時變性，結合線性化控制設計了并聯多內模和串聯高階內模切換的A-RC控制器，分析了閉環系統穩定性，提出基于頻率修正逆的重復控制補償器構造方法。
（3）利用搭建的手腕震顫抑制實驗平臺，通過健康參與者的誘發性震顫抑制試驗和意向性患者的手腕震顫抑制試驗，以期驗證所提A-RC方法的可行性和有效性。
1 方法
1.1 手腕肌骨模型
為了設計刺激強度和刺激時間精準調節的A-RC控制器，首先需要建立手腕肌骨模型。手腕肌骨模型由肌肉模型和骨骼模型組成。非線性Hammerstein模型由于其結構簡單，便于震顫抑制控制器的設計，被廣泛應用于基于FES的肌肉建模中[23]。Hammerstein模型由靜態非線性招募特性和線性動態模型組成，前者表現肌肉的激活特性，后者表現肌肉對電刺激信號的收縮響應。而骨骼模型的建立通常描述為線性剛體動力學（rigid body dynamics，RBD）模型。
本文主要研究手腕屈伸運動時的震顫抑制。腕部屈伸運動主要由橈側腕伸肌（extensor carpi radialis， ECR）和橈側腕屈肌（flexor carpi radialis，FCR）驅動橈骨實現。因而，根據Hammerstein模型結構，構造了FCR和ECR的非線性肌肉招募特性fFCR(uFCR)、fECR(uECR)和線性肌肉收縮動態特性GFCR(z)、GECR(z) ，以及骨骼線性剛體動力學模型GRBD(z)。電刺激手腕肌骨模型如圖1所示，其中uFCR(k)和uECR(k)是電刺激輸入信號，y(k)為手腕屈伸運動角度輸出信號。wFCR(k)和wECR(k)為招募特性的輸出， 為作用在骨骼系統上的肌肉總力矩。d(k)可以看作是由于震顫引起的干擾信號。
圖1
電刺激手腕肌骨模型
Figure1.
Electrically stimulated wrist musculoskeletal model
為簡化模型結構和參數辨識過程，將uFCR(k)、uECR(k)分別等效為如式(1)、式(2)所示的單一輸入：
|
|
其中，uFCR,ZO和uECR,ZO是FCR和ECR肌肉招募特性的死區閾值，根據試驗設置為50 μs。由于典型的招募特性包括初始死區、單調遞增區和飽和區，進而等效非線性招募特性輸出w(k)與輸入u(k)可以通過如式(3)所示的靜態映射關系表示：
|
其中，umax是最大電刺激輸入脈寬，為避免過強電刺激輸入造成肌肉收縮時的不適感，本文最大脈寬設置為300 μs [24]。 和 分別為FCR和ECR肌肉的非線性招募特性， 和 分別為 和 肌肉的招募特性死區閾值。 為等效的橈側肌肉非線性招募特性。化簡后的非線性招募特性f (u(k))在工作區間 是連續且遞增的，如式(4)所示：
|
其中， 以及 分別為FCR和ECR 肌肉需要識別的非線性標量參數；e為指數函數。
由于相似的肌群具有相似的線性收縮特性[25]，因而FCR和ECR動態特性近似為 ，可用GL(z)表示。人體骨骼動力學可以建模為剛體動力學GRBD(z)，其中阻尼和彈性函數都可以看作線性函數[26]。故，等效的線性動態肌骨模型 如式(5)所示：
|
其中， 和 分別為等效線性肌骨模型數學表達式的分子和分母，和分別是的階數和待辨識的參數，和分別是的階數和待辨識的參數。
1.2 控制器設計
由于腕關節肌骨模型具有固有的非線性，首先設計線性化控制器 ，在此基礎上設計A-RC控制器C(z)以抑制不同類型的腕部震顫，提高手腕對參考軌跡的跟蹤性能。基于FES的閉環震顫抑制系統控制框圖如圖2所示，其中r(k)為指定的跟蹤參考位置，e(k)為手腕實際運動位置與參考位置的偏差， 為控制器C(z)的輸出。
圖2
基于FES的閉環震顫抑制控制框圖
Figure2.
FES based closed-loop tremor suppression control diagram
1.2.1 線性化控制器設計
由于電刺激肌骨模型具有復雜的非線性因素，且非線性主要表現在靜態招募特性，因而采用線性化控制器補償肌骨模型中的招募非線性，線性化控制器 如式(6)所示：
|
1.2.2 A-RC控制器設計
重復控制是基于內模原理的一種控制方法。內模原理指出，如果在穩定的閉環系統中包含外部周期性信號的數學模型，系統就能夠實現對周期性信號的無靜差跟蹤或抑制。本文根據患者震顫數據的分析，設計了內模切換的A-RC控制器。經過線性化之后，等效的閉環控制系統框圖如圖3所示，H(z)是提高閉環系統穩定性和動態性能的重復控制補償器。自適應調節器根據傅里葉變換得到的震顫數據頻域特性調節內模結構，實現并聯多內模或串聯高階內模的自適應選擇，即根據震顫信號頻域分析結果，調節內模結構的數學模型，以實現對不同震顫患者最佳的震顫抑制效果。
圖3
等效的閉環震顫抑制系統框圖
Figure3.
Equivalent closed-loop tremor suppression system diagram
其中，A-RC控制器C(z)如式(7)所示：
|
其中，補償器H(z)如式(8)所示：
|
H(z)包含n個零點，其位置由系數c1, c2, ···, cn確定，并在原點引入n ? m個極點。自適應內模結構 由m條并聯支路構成，Wq為第q條內模支路增益，q = 1, 2, ···, m。其中，第q條高階內模支路 如式(9)所示：
|
其中， 為一個周期的采樣次數， 為震顫信號的周期， 為采樣周期，（p = 1, 2, ···, n）是第p階支路增益。
1.2.3 穩定性分析
定理1：在A-RC控制器的作用下，選擇控制增益 、 和補償器 使得閉環控制系統特征多項式的根均在單位圓內，如式(10)所示：
|
那么，如圖3所示的閉環控制系統是漸近穩定的，且能夠完全抑制 周期的擾動。尤其是，選擇 ，則當 時， 需要滿足如式（11）所示：
|
證明：根據圖(3)可得誤差 、參考輸入 和擾動 之間的關系如式(12)所示：
|
將式(7)帶入式(12)，則閉環系統的特征多項式如式(13)所示：
|
當特征多項式的根都在單位圓內時，則閉環反饋控制系統是漸近穩定的。特別地，當 時，式(13)可寫為如式(14)、式(15)所示：
|
|
為了有效抑制震顫，需要選擇重復控制器的內模結構 在震顫信號基波和諧波處增益無窮大，即 ，令 從而得到離散系統的頻率響應，其中， 表示虛數單位， 為角頻率，那么如式(16)所示：
|
故當選擇 ， 時，有 ，這時第 條內模支路增益 選擇滿足 。
證明完畢。
1.2.4 補償器H(z)設計
在實際中，由于手腕肌骨模型參數的不確定性和模型單位圓外可能存在零點等因素，補償器H(z)無法通過 得到。因此，本文采用頻率修正逆重復控制方法來獲得補償器H(z)的階數和參數[20, 27]。即適當選擇式（8）中的參數c=[c1, c2,…, cn]T來最小化以下成本函數J，從而實現 的頻率響應。成本函數J如式(17)所示：
|
其中， 為一個周期的采樣次數，i表示虛數單位， 是從零到奈奎斯特的一組離散頻率， 為采樣周期，*表示復共軛。
2 仿真與試驗驗證
在試驗驗證環節，招募4名健康參與者和4名意向性震顫患者。健康參與者編號為No.1～No.4，3男1女，平均年齡為（28.25 ± 8.54）歲。震顫患者來自鄭州大學第五附屬醫院康復科，編號為No.5～No.8，3男1女，平均年齡為（43.25 ± 10.90）歲。試驗中8位受試者均簽署書面知情同意書，且已通過鄭州大學倫理審查，倫理審批編號為No.ZZURIB2019-004。基于搭建的手腕震顫抑制實驗平臺，每位受試者首先進行手腕肌骨模型參數辨識試驗，然后根據辨識的模型參數設計相應控制器，進而開展手腕震顫數據采集和震顫抑制試驗。
2.1 實驗平臺
基于FES的腕部震顫抑制實驗平臺如圖4所示。震顫抑制實驗平臺由輪椅改裝而成，其中手臂支架是固定的，前端的U型手掌支撐夾板運動范圍被限制在± 80°之內，這符合手腕屈伸的運動范圍。試驗時受試者舒適、放松地坐在實驗平臺的輪椅上，右前臂放于手臂支架，右手放在U型手掌支撐夾板中間。受試者的前臂用綁帶固定，以防止產生除手腕屈伸平面外其他任何自由度的運動。支架外側設計扇形角度標尺，方便受試者跟蹤指定的參考位置。實時的手腕屈伸運動角度數據由編碼器（E6B2-CWZ6C，深圳布瑞特科技有限公司，中國）測量，其中編碼器的軸與U型手掌支撐夾板的轉軸通過直徑比為1∶1的齒輪連接。仿真系統平臺dSPACE（Microlabbox 1202, dSPACE Inc.，德國）實現手腕角度數據和控制信號的實時傳輸，為控制算法的開發和驗證提供支持。最終的控制信號由電刺激儀（4 CHANNEL STIMULATOR, Odstock Medical Ltd.，英國）輸出。本文電刺激信號的頻率、電流幅值和最大脈寬設置分別為40 Hz、20 mA和300 μs，并且可以根據受試者對電刺激信號強弱的不同反應，利用電刺激儀上每個通道的電位器調節信號的大小。
圖4
腕部震顫抑制實驗平臺
Figure4.
Experimental platform for wrist tremor suppression
在臨床試驗之前，對健康參與者進行控制方案的可行性測試是必要的步驟。因此，在實驗平臺上嵌入直流電機（VDM08SGN24-60-1800/JB5G12T，北京微特微電機制造有限公司，中國），利用電機來模擬震顫信號。其原理是U型手掌支撐夾板與電機同軸，當電機產生輸出轉矩時，可直接帶動U型手掌支撐夾板動作。因此，通過軟件設置震顫頻率和幅度，電機驅動器（ESCON50/5，maxon motor Inc.，匈牙利）驅動電機運轉即可輸出相應的轉速和角度。這樣，健康參與者的手腕會產生不自主的屈伸運動，即得到誘發性震顫信號。
2.2 手腕肌骨模型參數辨識
手腕肌骨模型的參數辨識基于輸入/輸出數據集［u, y］，利用迭代最小二乘法獲得，整個辨識過程大約2 min左右。輸入信號u(k)為0.2～4.0 Hz、脈沖寬度為300 μs的不同頻率正弦信號。y(k)為編碼器測量的手腕屈伸運動的角度數據。在參數辨識過程中，受試者手腕運動是由輸入的電刺激信號驅動，而不是自主運動產生。為了計算模型的擬合度，定義最佳擬合比 如式(18)所示：
|
其中，y為實測輸出， 為模型模擬輸出。y的均值用 表示。一般來說 值越高，參數辨識的精度越高，但模型階數也越高。在本文中， 和 的值［即 中 和 的階數］從4～10進行選取。以受試者No.1為例，其最佳擬合率如圖5所示。相對低階模型的最佳擬合率為62.06%（ 和 ），相對高階模型的最佳擬合率最高為67.94%（ 和 ）。當 和 在4～10之間取值時，最佳擬合率變化并不明顯。因此，為了減少辨識計算量、簡化控制器設計，本文選取了 。
圖5
受試者No.1手腕模型最佳擬合率
Figure5.
The best fit rate of participant No.1’s wrist model
2.3 震顫信號分析
利用搭建的實驗平臺采集意向性震顫患者No.5～No.8手腕屈伸運動時伴隨的震顫數據，如表1所示。
表1
患者No.5～No.8數據
Table1.
Data of Patients No.5～No.8
具體的采集方法為：讓意向性震顫患者右手通過屈伸手腕自主地跟蹤一個指定的參考軌跡，指定的參考軌跡類似于一些簡單的日常任務，例如，屈伸手腕以達到一個特定的角度。患者需要通過自主屈伸手腕跟蹤參考軌跡，產生的運動跟蹤軌跡為黑色虛線。由于震顫的存在，實際手腕運動的跟蹤軌跡會偏離指定的參考軌跡。如圖6所示，通過隨機選擇的患者No.5和No.7的跟蹤軌跡數據時頻圖可以看到，自主運動的頻率范圍低于1 Hz，且在20 s和40 s手腕有自主屈伸動作時，功率譜密度增強。意向性震顫頻率范圍主要在2～5 Hz之間[1, 28]，本文采集的意向性震顫頻率在2～3 Hz之間，且震顫頻率隨時間會有稍許的變化，如患者No.5震顫信號時頻圖40 s之后，震顫頻率明顯降低。
圖6
震顫患者No.5和No.7震顫信號頻率分析
Figure6.
Frequency analysis of tremor signals in patient No.5 and No.7
通過對患者No.5和No.7跟蹤軌跡信號分析可得，患者No.5為單頻率震顫信號，而患者No.7為多頻率震顫信號。對于單頻率或多頻率的震顫信號，通過自適應調節器選擇切換不同的重復控制器內模結構，即通過傅里葉變換得到最大振幅對應的震顫頻率（frequency corresponding to peak amplitude，FPA），如患者No.5的FPA為2.2 Hz，患者No.7的FPA為2.3 Hz。若其余頻率振幅均小于FPA的80% （患者No.5沒有超過FPV的80%的震顫頻率，故歸類為單頻率震顫信號），則選擇單一支路的高階內模結構（圖3中q = 1）；若震顫頻率幅值中有高于80%的FPA頻率（患者No.7在2.5 Hz處的震顫幅值超過FPV的80%，故歸類為多頻率震顫信號），則選擇并聯內模結構（圖3中q > 1），因為震顫頻率的振幅越大，對患者的影響越大。
2.4 仿真測試
在本小節中，使用受試者No.1的手腕模型驗證提出的A-RC方法對周期性震顫抑制的有效性，并將并聯多內模結構和串聯高階內模結構切換的A-RC方法與單周期重復控制方法（single periodic repetitive control，SP-RC）[16]和傳統的基于高通濾波器的比例積分控制方法（proportional integral control algorithm based on high pass filter，PI-HPF）[12]進行比較。根據震顫信號的單頻率或多頻率特點，將仿真測試分為以下兩個環節。
2.4.1 仿真環節1結果與分析
仿真測試環節1中，震顫信號d(k)采用頻率為2 Hz、幅值為1的正弦信號。采樣周期TS為0.005 s，因此A-RC控制器q = 1，選擇單一支路的串聯高階內模結構。W1 = 1，內模延時周期N1 = 100，p = 2，K1 = 2，。對于單周期重復控制器，圖3中q = 1，W1 = 1，N1 = 100，p = 1且K1 = 1。重復控制器補償器H(z)由頻率修正逆算法得到。采用臨界比例度法確定了PI-HPF調節器的參數（P = 36.82，I = 1.4），高通巴特沃斯濾波器階數為6，截止頻率為1.2 Hz。設定參考輸入r(k)為零，輸出y(k)在干擾d(k)的作用下跟蹤參考信號。
根據采集的患者震顫數據頻譜分析可得，震顫信號的峰值頻率隨時間會產生稍許變化，故在仿真環節進行蒙特卡洛實驗，以驗證所提控制器對于頻率變化的周期性震顫信號的抑制效果。令震顫信號d(k)在峰值頻率±5%的范圍內隨機波動（即1.9～2.1 Hz隨機變化），控制器參數不變，進行100次仿真實驗，如圖7所示。圖7左側圖為仿真環節1中不同控制器下一次仿真測試的跟蹤誤差結果，PI-HPF方法的控制輸入相對于兩個重復控制器（A-RC和SP-RC）小一些，但是震顫抑制效果最差，因而驗證了對于周期性擾動信號重復控制器的優越性。對比A-RC和SP-RC兩種重復控制方法，在震顫峰值頻率有所變化的情況下（設定2 Hz，實際2.05 Hz），本文所提的A-RC控制器跟蹤誤差信號最小，震顫抑制效果最佳。
圖7
不同控制器下的跟蹤誤差和控制量
Figure7.
Tracking error and control input under different controllers
通過以下兩個指標來量化各控制方法在震顫抑制方面的性能。
（1）均方根誤差（root mean square error，RMSE）的計算公式如式(19)所示：
|
其中， 為輸出信號， 為參考信號，n為測試的總樣本數。
（2）震顫抑制率（tremor suppression rate，TSR）如式(20)所示：
|
其中， 是輸入信號 和輸出信號 的偏差。 為震顫信號對跟蹤角度的影響。 是輸入信號 和 之間的偏差。 與 的比值為震顫抑制性能。
仿真環節1蒙特卡洛實驗的統計結果如圖8左側圖所示。PI-HPF方法受震顫信號頻率變化影響最小，其次是本文所提的A-RC方法，而SP-RC方法最差。但PI-HPF方法的RMSE均值（0.357）卻遠遠大于A-RC和SP-RC方法。SP-RC方法的RMSE均值為0.095，本文所提A-RC方法RMSE均值僅為0.043。RMSE 均值越小，意味著跟蹤誤差越小。在TSR的統計指標中，A-RC方法、SP-RC 方法和PI-HPF方法TSR均值分別為93.90%，85.26%和49.58%，TSR指標越大，震顫抑制效果越好。因此，本文所提的A-RC方法對于單頻率變化的周期性震顫信號具有最佳的抑制效果。
圖8
蒙特卡洛實驗統計結果
Figure8.
Statistical results of Monte Carlo trials
2.4.2 仿真環節2結果與分析
在仿真測試環節2中，震顫信號d(k)采用頻率為2 Hz和2.5 Hz，幅值為1和0.8的疊加正弦信號。A-RC控制器切換到并聯多內模結構，q = 2，W1 = W2 = 0.5，內模延時周期分別為N1 = 100和N2 = 80，K1 = K2 = 0.5，其余參數與仿真環節1一致。對于SP-RC控制器和PI-HPF控制器參數均不變。
在不同控制器下系統跟蹤誤差如圖7仿真環節2所示，控制量和仿真環節1的結果類似，PI-HPF方法的控制量最小，A-RC和SP-RC方法的控制量相差不大。但是SP-RC方法的跟蹤誤差和仿真環節1相比有明顯增加，因為單內模結構只能抑制內模中已有的周期性信號，而對于2.5 Hz的信號無法抑制。相比較而言，本文所提的A-RC方法可以實現多頻率擾動信號的完全抑制。利用蒙特卡洛實驗討論在震顫頻率變化情況下，控制器對于多頻率震顫信號的抑制性能。震顫信號在峰值頻率±5%的范圍內波動，100次仿真測試的統計結果如圖8仿真環節2所示。從統計結果分析來看，PI-HPF和SP-RC方法的RMSE均值接近，分別為0.408和0.394，這說明在震顫頻率變化的情況下，這兩個方法的穩態跟蹤誤差大，TSR指標僅為48.36%和50.21%，無法實現多頻率震顫信號的抑制。而本文所提的A-RC控制器，利用其并聯多內模結構，在頻率變化的情況下，穩態誤差小，TSR均值可達90.27%，即實現對多頻率震顫信號的良好抑制。
2.5 試驗驗證
試驗驗證也分為兩個階段，第一階段在健康參與者身上開展誘發性震顫抑制試驗，用以評價所提控制方法的可行性。第二階段招募意向性震顫患者進行震顫抑制驗證試驗，進一步驗證所提方法的有效性。
2.5.1 試驗階段1結果與分析
在試驗階段1中，通過震顫抑制實驗平臺產生誘發性震顫。為了防止肌肉疲勞，每個受試者的電刺激時長不超過20 min。4名健康參與者需要完成以下三個測試：
測試1（test 1，T1）：自主地屈伸手腕跟蹤參考軌跡，此測試無誘發性震顫和電刺激輸入。
測試2（test 2，T2）：自主地屈伸手腕跟蹤參考軌跡，但此測試加入誘發性震顫信號。
測試3（test 3，T3）：在具有誘發性震顫和功能性電刺激輸入的情況下，自主屈伸手腕跟蹤指定參考軌跡（每位受試者在同一控制器下進行3次T3試驗）。
為了更好地對比不同控制器下T3的試驗結果，將誘發震顫設置為多頻震顫，頻率為2 Hz和2.5 Hz。受試者No.1在不同控制器作用下其中一次試驗的手腕運動軌跡和電刺激控制輸入信號如圖9所示。
圖9
受試者No.1在不同控制器下的手腕角度位置和電刺激信號
Figure9.
Wrist angular position and electrical stimulation level of participant No.1 under different controllers
如圖9中的T3手腕跟蹤軌跡所示，3種不同的控制器均可以不同程度地實現誘發性震顫的抑制，但本文所提的A-RC方法的T3運動軌跡相對參考軌跡的偏離程度最小，具有最佳的震顫抑制效果。其次是SP-RC控制器，PI-HPF控制器的震顫抑制性能最差。對于電刺激控制輸入，PI-HPF控制器產生的電刺激控制輸入信號相對重復控制器要小一些，但震顫抑制性能差。A-RC控制器和SP-RC控制器的電刺激輸入信號相差不大，結合T3的手腕運動軌跡可得本文所提A-RC方法具有震顫抑制的性能優勢。如圖10所示，受試者No.1在不同控制器下T3的跟蹤軌跡時頻圖也可以得到相同的結論。在不同控制器作用下，多頻率震顫信號的功率譜密度在震顫頻率處均有不同程度的下降。其中，本文所提A-RC方法的功率譜密度下降最為明顯，其次是SP-RC方法，而PI-HPF方法對于周期性震顫信號的抑制作用最差。
圖10
參與者No.1在不同控制器下跟蹤軌跡的時頻圖
Figure10.
The time-frequency diagram of tracking trajectory of participant No.1 under different controllers
試驗環節采用與仿真環節相同的性能指標。如式（19）所示，RMSE值用于量化T3（）與T1（）的之間的偏離程度。如式（20）所示，是T3（）與T1（）之間的偏差。是T2（）與T1（）之間的偏差。其余三位健康參與者也得到了相似的試驗結果，統計結果如圖11所示。
圖11
四位健康參與者試驗統計結果
Figure11.
Statistical results of experiment with four participants
四位健康參與者RMSE均值A-RC方法最小為1.34，而SP-RC方法為3.12，PI-HPF方法為3.58，這一試驗結果與仿真結果類似。通過TSR值可以更為直觀地判斷每種控制方法對周期性震顫信號的抑制情況，如圖11所示四位健康參與者在本文所提的A-RC方法下震顫抑制性能顯著提升，四位健康參與者TSR均值為87.59%，比SP-RC方法和PI-HPF方法均有明顯提升。
2.5.2 試驗階段2結果與分析
本階段招募意向性震顫患者開展震顫抑制試驗，定義以下三個測試試驗：
測試4（test 4，T4）：將指定參考軌跡提供給患者，患者不需任何動作。
測試5（test 5，T5）：患者自主地屈伸手腕跟蹤參考軌跡，此測試無FES信號。
測試6（test 6，T6）：在FES作用下患者自主屈伸手腕跟蹤參考軌跡（同一控制器下患者需進行3次T6試驗）。
隨機選擇患者No.5和No.7的T6跟蹤軌跡時頻圖如圖12所示。從時頻圖可以看到，本文所提的A-RC方法，對于單頻率震顫和多頻率震顫均有很好的抑制效果，在震顫頻率處的功率譜密度均有明顯下降，比PI-HPF方法在震顫抑制性能方面有明顯提升。而SP-RC方法對于周期性震顫信號也有一定的抑制作用，但效果相對不佳。患者試驗環節也采用相同的性能指標。如式（19）所示，RMSE值用于量化T6（）與T4（）的之間的偏離程度。如式(20)所示， 是T6（）與T4（）之間的偏差。 是T5（）與T4（）之間的偏差。患者No.5和No.7的量化試驗統計結果如圖13所示。
圖12
患者No.5和No.7在不同控制器下T6跟蹤軌跡的時頻圖
Figure12.
The time-frequency diagram of T6 tracking trajectory of patients No.5 and No.7 under different controllers
圖13
意向性震顫患者試驗統計數據
Figure13.
Statistical experimental results of intention tremor patients
如圖13所示，患者震顫抑制試驗結果與第1階段健康參與者的結果類似。A-RC方法對于單頻率和多頻率，且頻率存在微變的震顫均有很好的抑制作用。A-RC方法實現平均82.53%的單頻震顫和82.21%多頻震顫抑制，而SP-RC方法和PI-HPF方法對于這兩種類型震顫的抑制率僅分別為71.22%、58.31%和67.88%、52.33%。
3 結論
本文提出一種基于FES的A-RC方法來抑制腕關節意向性震顫。首先建立了非線性Hammerstein結構的手腕肌骨模型，提出了自動、快速的參數辨識算法，在搭建的手腕震顫抑制實驗平臺完成了肌骨模型的參數辨識。其次，根據辨識的每位受試者肌骨模型設計相應的線性化控制與反饋A-RC相結合的控制器，并根據患者震顫頻率特點自動選擇重復控制內模結構，分析了閉環系統穩定性，并利用頻率修正逆算法實現了重復控制補償器的構造。最后，通過健康參與者誘發性震顫抑制試驗和意向性患者震顫抑制試驗驗證了本文提出的A-RC策略，其平均震顫抑制率可達84.98%，比傳統單內模重復控制和基于濾波的反饋控制均有明顯提高。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：張贊主要負責數據收集與分析、算法程序設計以及論文編寫；霍本巖主要負責試驗素材收集；劉艷紅、楚冰和Owens David Howard主要負責試驗指導、數據分析指導以及論文審閱修訂。
倫理聲明：本研究通過了鄭州大學倫理委員會的審批（批文編號：No.ZZURIB2019-004）。
