為了使被動平衡康復訓練系統可以制定安全的訓練強度和訓練方式,本論文采用T-S模糊辨識方法建立了一種新的人體站立平衡系統數學模型。該模型以多維運動平臺的運動加速度為模型輸入,人體關節角度為模型輸出。采用人工蜂群尋優算法改進模糊C-均值聚類算法,提高了辨識前件參數的效率。通過試驗,采集了9位健康成年人的被動站立平衡調節數據,用于模型參數訓練和模型結果驗證。采用仿真結果和測量數據的均方差及互相關度,證明了所建模型是準確和合理的。
引用本文： 王洪瑞, 劉琨, 肖金壯, 熊鵬. 一種基于T-S模糊辨識的人體站立平衡系統建模方法. 生物醫學工程學雜志, 2014, 31(6): 1243-1249. doi: 10.7507/1001-5515.20140236 復制
引言
站立平衡功能關系到人們的生存質量，近年來關于人體站立平衡功能的康復設備和訓練方法逐漸成為研究的熱點[1-2]。本課題試圖開發一種被動平衡康復訓練與評測系統，該系統由多維運動平臺提供外力，干擾平臺上訓練者的站立平衡狀態，使訓練者被動地鍛煉自身平衡調節機能[3]。這種康復訓練設備存在的問題是外力的刺激強度難以選擇，如果外力強度過低則難以達到康復治療訓練的效果，而外力強度過大又會使訓練者摔倒或者造成關節組織扭傷，此外患者的身體條件不盡相同也導致系統安全強度需因人而異進行調整。在康復治療必須保證患者安全這一前提下，為了進一步完善該被動訓練系統的訓練方法和訓練策略，需要建立一種可體現測試者個體差異的人體站立平衡調節數學模型，仿真人體主要關節角度在平臺運動下的動態響應過程，最終預測平臺運動對人體刺激的安全強度范圍。
一些研究者提出了基于工程模型的建模方法，將人體的肌肉骨骼簡化為倒立擺，并為其建立控制器和控制算法，模擬人類中樞神經系統(center neural system,CNS)的功能[4]。還有研究者利用一些實驗或醫學理論證明了CNS在功能上可近似等效為比例-積分-微分控制器 (proportion-integration-differentiation controller，PID)、比例-微分控制器(proportion-differentiation controller，PD)、最優控制器 (Optimal controller)或間歇控制器(Intermittent controller)等[5-7]。這類工程模型雖然可以近似地為研究者提供模擬人類運動控制的數學依據，但仿真結果不能夠反映不同人的運動特點，難以應用于工業產品開發或醫療產品設計。
本研究采用T-S模糊辨識方法，基于9位健康成年人的真實測試數據建立數學模型，以平臺的運動加速度數據為輸入信號，人體的主要關節角度動態數據為系統的輸出信號。該方法不需要考慮人體結構的耦合關系，僅將人體的中樞神經系統和人體骨骼肌肉系統視為內部結構未知的黑箱。對T-S模糊模型的前件參數的辨識，采用的是模糊C-均值聚類(fuzzy C-means，FCM)算法，但由于FCM算法對聚類中心初值要求苛刻，聚類結果易陷入局部最小值，導致模型的辨識精度和效率降低[8]，所以本研究采用人工蜂群尋優(artificial bee colony，ABC)算法對聚類原型的初值進行優化。由于ABC具有較好的全局優化能力[9]，克服了聚類算法的缺陷，因此保證了T-S模型前件參數的正確性。
最后本文利用遞歸最小二乘法，求解模型的后件參數，并推導出模型的解析函數。本模型可以有效表征不同受試者個體間的運動特點和差異，并通過計算仿真數據和測量數據的均方差(mean square error，MSE)以及互相關度，證明所建模型的準確性和客觀性。
1 模型辨識
1.1 T-S模糊模型
T-S模糊模型可以表示復雜的系統，以很高的精度逼近定義在緊集上的非線性函數，并且模型的結果用線性方程表示[8]，適合于人體站立平衡系統建模。辨識T-S模糊模型的前件參數采用FCM聚類算法，將數據按照一定的相似性準則劃分成多個子集，以表征系統的不同性征，并對集群的中心進行編碼。因此，基于FCM的T-S模糊模型，可以表示為
| $\begin{align} & {{R}_{i}}:\text{ }IF\text{ }x\text{ }is\text{ }in\text{ }cluster~{{C}_{i}}~with\text{ }({{v}_{i}},{{u}_{i}}), \\ & THEN~{{y}^{i}}={{a}_{0}}^{i}+{{a}_{1}}^{i}{{x}_{1}}+{{a}_{2}}^{i}{{x}_{2}}+\ldots +{{a}_{d}}^{i}{{x}_{d}}, \\ \end{align}$ |
其中 i =1,2,…,R代表聚類中心的序號；向量x=(x1,x2,…,xd)為給定的輸入變量；Ci為模糊變量產生的模糊聚類空間；d為輸入變量維數；θi=(a0i,a1i,a2i,…,adi)代表結論部分的輸出參數；yi為第i個模糊規則的輸出結果。模糊模型的輸出Y是一個加權平均數去模糊化計算方法，可表示為
| $Y=\frac{\sum\limits_{i=1}^{R}{{{w}_{i}}{{y}^{i}}}}{\sum\limits_{i=1}^{R}{{{w}_{i}}~}},$ |
其中wi為第i個模糊規則輸出結果的權值，即
| ${{w}_{i}}=exp-{{\left( \frac{\|x-{{v}_{i}}\|}{{{\sigma }_{i}}} \right)}^{2}}~,$ |
其中σi為常數，取σi=2；vi是輸入變量x的聚類中心坐標位置，即
| ${{v}_{i}}=\frac{\sum\limits_{i=1}^{R}{{{u}_{ij}}{{x}_{i}}}}{\sum\limits_{i=1}^{R}{{{u}_{ij}}~}}~~,$ |
| ${{u}_{ij}}=\frac{1}{\sum\limits_{k=1}^{R}{\left( \frac{\|{{x}_{i}}-{{v}_{j}}{{\|}^{2}}}{\|{{x}_{i}}-{{v}_{k}}{{\|}^{2}}} \right)}}$ |
其中uij為隸屬度值，下標j=1,2,…,N表示數據樣本序號，k=1,2,…,R，且k≠i。
1.2 人工蜂群改進的FCM聚類算法
FCM通過對目標函數的迭代優化來取得對數據集的模糊聚類，目標函數在迭代過程中是遞減的。這種方法在很大程度上依賴于初始聚類中心的選擇，如果初始聚類中心矩陣選擇不合理，會造成錯誤聚類，還會增加系統復雜度，降低算法的效率。
ABC是模擬蜜蜂采蜜方式的仿生式人工智能優化算法，具有引領蜂、跟隨蜂和偵查蜂為三種不同功能的尋優因子。每個食物源的坐標對應樣本的一個可行解(聚類中心)，食物源的收益度對應問題的適應度，最大收益度對應問題的最優效果。這里的適應度用于評價各個聚類中心在優化過程中的價值，評價這個可行解的質量，作為是否需要保留的條件，即
| $fi{{t}_{n}}=\frac{1}{1+objec{{t}_{n}}},$ |
其中objectn是要解決優化問題的目標函數，由具體問題而定。ABC計算步驟如下：
(1)初始化參數：設定尋優問題的食物源個數N，最大循環次數100，最大迭代次數200，循環次數m和迭代次數n的初值均為0。
(2)產生初始解集xij，并且通過式(6)計算初始解集的適應度。
(3)引領蜂通過式(7)做鄰域搜索，產生新解vij，即
| ${{v}_{ij}}={{x}_{ij}}+{{\Phi }_{ij}}({{x}_{ij}}-{{x}_{kj}})\text{ },$ |
并且計算其適應度，如果比上一次xij的適應度大，則用vij取代xij，否則保留原xij；
(4)利用新解xij的適應度，根據式(8)計算概率值Pij，用以選擇食物源，即
| ${{P}_{ij}}=\frac{fi{{t}_{n}}}{\sum\limits_{n=1}^{200}{fi{{t}_{n}}}}$ |
(5)鄰域搜索，產生新解vij，計算其適應度，若大于xij的適應度，則用vij取代xij，否則xij保留。
(6)鄰域搜索計數到達最大循環次數后，判斷是否有要丟掉的解，如果存在，偵查蜂根據式(9)產生一個新解代替原xij。
| ${{x}_{ij}}=xmin+rand\left( 0,1 \right)\times (xmax-{{x}_{min}})$ |
(7)判斷當前的迭代次數是否已經到達最大迭代次數，如果到達則停止迭代，輸出最后一代的食物源作為優化問題的最優解集合，否則轉到步驟(3)，迭代計數加一。
ABC改進FCM算法(ABC-FCM)算法的思想可總結為：利用ABC算法良好的全局搜索能力求得最優解，作為FCM算法的初始聚類中心，再利用FCM算法求得最終聚類中心解。ABC優化目標函數和FCM聚類算法的聚類目標函數均選為式(10)，尋優與聚類使objectn最小。
| $objec{{t}_{n}}=\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{C}{{{u}_{ij}}\|{{x}_{i}}-{{x}_{j}}{{\|}^{2}}}}$ |
ABC-FCM算法即在原ABC優化算法的基礎上加入以下步驟：
(8)根據式(5)計算新隸屬度函數U={uij|j=1,2,…,N; i=1,2,…,R}。
(9)根據式(4)計算新聚類中心V={vi|i=1,2,…,R}。計算相鄰兩代隸屬度之差E，若E＜ε，則尋優停止，否則轉到步驟(8)。
1.3 T-S模型辨識算法
所研究的人體站立平衡系統選取了平臺在水平面的加速度數據為輸入，人體踝、膝、髖關節的角度變換為輸出，組成了二維輸入三維輸出問題。為了便于模型的辨識，本文采取的方法是將三個關節獨立建模，即將原問題視為分別針對踝、膝、髖關節的三個二維輸入一維輸出問題，同時為了進一步提高辨識準確性，將另外兩關節過去時刻的角度數據作為當前辨識計算的輸入，這種方法在實際應用中也易于實現。以辨識踝關節為例介紹具體算法：
(1)T-S模型參數訓練
輸入： 令Z=(A前后(k)，A左右(k)，φH(k-1),φK(k-1),令y=φA(k)，A前后(k)，A左右(k)是運動平臺所在水平面的前后和左右方向的加速度數據；φH(k-1),φK(k-1)人體踝、膝關節前一時刻的角度數據。
輸出： T-S模型前后件參數：{U,V,θ}。
第1步： 設定聚類數R，按照ABC-FCM聚類方法，利用式(4)、(5)，分別求解聚類中心V和隸屬度函數矩陣U。
第2步： 用最小二乘算法求解后件參數。通過目標函數最小化：
| $J=min\frac{1}{N}\left\{ {{(y-Xe\theta e)}^{T}}\Phi i(y-Xe{{\theta }_{e}}) \right\}$ |
其中Xe=[1,X]，Φi=diag(ui1,ui2…um)。因此，模型在本條規則的輸出的參數可以得到：
| $\theta i={{({{x}_{e}}^{T}\Phi i{{x}_{e}})}^{-1}}{{x}_{e}}^{T}{{\Phi }_{iy}}$ |
(2) T-S模型輸出方程
輸入： ① T-S模型結構參數{U,V,θ}；② 受試者在平臺同方向以3 cm幅度運動所獲得的數據樣本Z=(A前后(k)，A左右(k),φH(k-1),φk(k-1))；③ 實際輸出y。
輸出： 預測輸出函數Y=φA(k)和MSE。
第1步： 計算樣本Z和聚類中心向量V的距離。
第2步： 根據式(3)，計算樣本Z對每個本地輸出函數yi的隸屬度，這個隸屬度值wi可以看作這個規則下輸出yi的權值。
第3步： 計算每個規則的輸出相應值yi，利用當前模型的參數θi，根據式(2)可以計算得到輸入為Z的估計輸出Y。
第4步： 利用均方差式(13)，計算該仿真結果與真實數據之間的MSE。
| $MSE=\frac{1}{N}\sum\limits_{j=1}^{N}{{{({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{y}}}_{j}}-{{y}_{j}})}^{2}}}$ |
2 實驗設計
T-S模糊辨識方法需要真實人體數據對模型結構和參數進行訓練，另一方面需要數據對辨識結果進行驗證[8]，所以本節詳細介紹試驗設計與數據采集方法。
2.1 受試者和測試設備
受試者包括9名健康成年人(4女，5男)，年齡(25.2±2.1)歲；身高(169.1±4.9 )cm；質量(62.6±9.7)kg。受試者沒有肌肉或神經運動障礙病史，并且進行了相關臨床檢測[10]。
所用到的測試設備包括：多維運動平臺。四臺高分辨率攝像機。計算機，用于視頻信號記錄和圖像處理。基于ADIS16355陀螺儀傳感器的平臺運動測量裝置，該裝置固定于平臺上，用于采集平臺的運動加速度數值。人體重心壓力坐標采集板[11]。
2.2 實驗方法
測量受試者的身體信息：將人體以踝、膝、髖關節的中心點分為上中下三個體段，A表示踝關節；K表示膝關節；H表示髖關節；L表示下段(小腿)；M表示中段(大腿)；U表示上段(上段包括軀干、手臂和頭部)[11-12]。mL、mM、mU分別指三體段的質量；cL、cM、cU分別指三體段質心到轉動關節的距離；lL、lM、lU分別是三體段的長度；參數Mwhole人體的總質量[13]。
受試者赤腳站立在運動平臺上，雙手自然下垂，放于身體兩側，兩腳并攏，如圖 1所示。測試時要求受試者閉上雙眼，不可抬腳、擺動手臂，僅依靠身體的主要關節來調整平衡姿態。平臺的水平運動數據通過陀螺儀加速度傳感器模塊采集，記錄為A前后(k)，A左右(k)。人體重心壓力坐標采集板將每次試驗過程中人體壓力中心(center of pressure,COP)的坐標記錄下來，作為模型仿真結果的驗證依據。
圖1
受試者的站立位置和姿態
Figure1.
Position and posture of the testee
通過文獻[3]采用的線性分析試驗可知，運動平臺的擾動形式為“水平面內的脈沖式平動”，幅度為2~3 cm，時間帶寬約為200 ms。本文所說外力強度可對應理解為平臺單次脈沖運動的幅度，該擾動幅度不但使人體晃動數據具有較高的信噪比，也保證了受試者的安全。本文將2 cm幅度擾動條件下的人體關節角度動態數據用于T-S模型系統辨識參數訓練；將3 cm幅度擾動條件下的人體關節角度動態數據作為辨識獲得模型預測值的驗證數據。由于T-S模型不需要考慮人體的結構耦合關系，所以平臺在任意方向的運動均可以利用該方法建模。用于參數訓練的人體運動數據越全面，T-S模型的預測精度就越高，但是在臨床使用中，不宜對患者進行過久的測試，可能會造成模型參數訓練不足，降低了模型的準確性。為了盡可能使樣本具有遍歷性，平臺以水平面內的脈沖式多方向運動，各方向間隔45°，如圖 2圓盤內箭頭所示，分別為：前(anterior,A)，后(posterior,P)，左(left,L)，右(right,R)，左前(anterior-left,AL)，左后(posterior-left,PL)，右前(anterior-right,AR)，右后(posterior-right,PR)。
圖2
平臺運動方向和攝像機位置
Figure2.
Platform moving direction and camera position
攝像機與受試者髖關節中心點的高度一致。4臺攝像機的間隔角度為45°，朝向平臺的中心，如圖 2所示。均在人體朝向攝像機方向的肩部、髖關節、膝關節和踝關節貼有4組標識點，用于身體運動位置的視頻提取。平臺每次運動方向所在矢量平面與人體的晃動平面相同，并且總有一臺攝像機的鏡頭方向與人體的晃動平面保持垂直，使對圖像的觀測誤差盡量小。測量開始后，平臺的運動方向按照之前規定的激勵方式進行，給定的平臺運動激勵間隔為10 s。測試過程中，受試者可能會習慣平臺激勵，進而做出提前應激動作，這將破壞數據的客觀性，為此平臺的運動方向順序是隨機的。待8個方向進行完畢后，按同樣方法連續測試三次。
通過圖像處理技術，提取所拍攝畫面中標記點的坐標。將人體簡化為矢量平面內多連桿體段倒立擺，將每幀畫面關節標識點的位置與初始點位置比較，可以得到人體運動狀態。以求解踝關節角度變換數據為例：圖 3為2號攝像機得到的人體運動下肢部分簡圖。2號攝像機正對受試者右側，圖 3中關節處白色圓形為提取得到的標識點，定義體段與鉛垂線的夾角為其相鄰轉動關節的關節角，將相鄰關節標識點的水平運動距離相減可得運動距離stA-st0A-sP=ΔsA，stK-st0K-sP=ΔsK，ΔsAK=ΔsA-ΔsK，其中st0A、stA表示踝關節在單位時間開始和結束的橫坐標，sP表示單位時間內平臺運動的距離，st0K、stK表示膝關節在單位時間開始和結束的橫坐標。踝關節運動角度數據用反三角正弦函數可得，即
| ${{\theta }_{t}}^{A}=\frac{arcsin\left( \frac{\Delta {{s}^{AK}}}{{{l}_{L}}} \right)}{t-{{t}_{0}}}$ |
圖3
踝關節和膝關節標致點運動簡圖
Figure3.
Motion diagram of the mark points on ankle and knee
膝關節和髖關節可用同樣方法獲得動態角度數據集，關節數據分別記做踝：(φA)、膝(φK)、髖(φH)。該計算方法已在Matlab7.0軟件環境編程，可以批量處理試驗獲得的人體數據。
2.3 實驗步驟
(1)試驗前向受試者說明試驗任務，確保受試者正確理解試驗操作并進行熟悉性練習后，簽訂試驗協議。
(2)受試者填寫個人基本資料，測量身高、體段及體重等基本信息后，在身體規定位置，安裝標識點。
(3)受試者按要求于運動平臺上站好。確認無誤并提示受試者后，啟動平臺，開始施加運動刺激。
(4)平臺運動20 s，確認受試者無不良反映，開啟采集設備電源，進行數據采集；若受試者出現不適，應立即停止平臺運動刺激。
(5)按要求方法采集完畢后結束試驗，關閉數據采集系統，輸入平臺停止運動指令，待平臺停止運動過后，令受試者離開平臺。
(6)完成基本數據的計算，編號后存入硬盤。
3 結果與分析
3.1 聚類算法結果與分析
通過對試驗獲得的實際人體角度數據矩陣X=[φA,φK,φH]進行聚類計算，結果如圖 4所示。
圖4
聚類計算結果
Figure4.
Clustering results of calculation
如圖 4(上)所示，*表示聚類的輸入數據集合；△表示ABC優化的聚類初始點；○表示FCM聚類的最終聚類中心點。通過ABC算法的改進，聚類中心分布比較均勻，FCM算法在優化的聚類初始點基礎上進行微小調整即達到最終聚類結果。聚類目標函數的迭代結果比較如圖 4所示，ABC-FCM的初始迭代效果明顯優于FCM，聚類目標在迭代前100次便達到較低的數值，之后趨于平穩。程序設定ABC最大迭代次數為200次，圖 4(下)中實線為ABC-FCM聚類，200次ABC優化迭代后轉為FCM繼續聚類目標函數，結果穩定在2 500左右。實曲線200~250次FCM迭代對應于圖 4(上)中的變化，即FCM算法在ABC聚類初始點基礎上進行的微小調整。虛線為單純利用FCM聚類算法得到的目標函數值，聚類效果較差，250次迭代運算后結果處于3 000左右，從該結果可以看到ABC-FCM的聚類效果更加合理。由于T-S模型辨識的前件參數辨識對準確性有較大的要求，這種改進對模型建立的準確性有較大改善。聚類程序運算過程中，ABC-FCM聚類相對于FCM聚類占用了更多的系統時間，在同樣的計算機和Matlab7.0仿真環境中，ABC-FCM聚類計算占用2.321 s，而FCM聚類計算僅占用0.283 s。可見ABC優化算法對聚類性能的改善是以消耗系統時間為代價的，但是考慮到本課題選用的T-S模糊模型辨識屬于離線辨識方式，并不要求較高的系統辨識速度，因此該問題對于本課題的影響并不明顯。
3.2 辨識結果與分析
隨機選取一位受試者在平臺后向運動過程中記錄的試驗數據建立模型，并且對比仿真運動結果，踝、髖、膝三個關節的實際運動過程與對應關節仿真的結果如圖 5所示。
圖5
關節數據與仿真結果對比
Figure5.
Comparison of joint data with simulation results
如圖 5所示，橫坐標為測試數據樣本序號，該樣本的持續時間從平臺運動開始到1500 ms為止。仿真結果與真實數據之間的MSE分別為：踝關節的MSE=0.017；髖關節的MSE=0.018；膝關節的MSE=0.220。仿真結果的誤差處于較小的范圍，三個關節的運動幅度可以被如實預測，同時各個關節之間的運動協調關系基本與真實人體數據一致。基于同樣的方法，可以求得該受試者在平臺刺激為其它方向時的人體數學模型，并且仿真出平衡調節運動數據；通過計算這些仿真數據和真實測試數據的MSE，可知全部仿真結果的MSE＜0.271。通過式(15)可以計算相應的壓力中心點COP坐標變化曲線[11-12]，
| $\begin{align} & CO{{P}_{前后}}=-[({{m}_{L}}{{c}_{L}}+{{m}_{M}}{{l}_{M}}+{{m}_{U}}lU)\cdot sin({{\varphi }_{A}})+ \\ & ({{m}_{M}}{{c}_{M}}+{{m}_{U}}{{l}_{U}})\cdot sin({{\varphi }_{A}}+{{\varphi }_{K}})+ \\ & {{m}_{U}}{{c}_{U}}\cdot sin({{\varphi }_{A}}+{{\varphi }_{K}}+{{\varphi }_{H}})]/{{M}_{whole}} \\ \end{align}$ |
計算仿真數據和采集實際重心坐標數據的互相關程度。通過對9位受試者的數據仿真，可以得到所有受試者的仿真結果與實際采集得到COP數據的相關度，全部結果在0.641~0.785之間。一般認為，當COP坐標移動到足底與平臺接觸面之外，人體就會失去平衡，現以不同受試者的足底長度為其前向判斷摔倒閾值，用所建立模型預測數據進行判斷，可獲得不同受試者的安全刺激強度范圍。9位受試者安全刺激強度的預測范圍處于5~7.5 cm之間。
需要指出，求解COP的過程存在三種誤差來源：① 本研究所用T-S模糊模型辨識方法受到迭代次數的限制，預測結果不可能零誤差逼近實際測量數據；② 式(15)涉及諸多肢體參數，測量時采用了文獻[13]介紹的方法，由于人為操作的問題同樣會引入測量誤差；③ 人體密度分布不均勻，形狀不規則，各體段質心的位置以現有的測量手段難以獲得，而本研究采用了基于統計數據的比例法，并且將各個體段等效為密度均勻的圓柱體。上述三方面誤差累積，會在一定程度上降低COP結果的準確性，但是較以往的人體運動模型，該方法所建立的人體站立平衡模型能夠較準確地預測人體運動狀態，并且仿真結果足以體現個體間的差異，可針對不同患者制定安全刺激強度。
4 結論
本研究針對目前人體平衡系統建模方法存在的不足，通過T-S模糊模型辨識方法，建立了平臺運動加速度和人體踝、髖、膝三個關節運動的數學模型，并且用人工蜂群算法改進了模型前件參數辨識的運算效率和聚類精度，對9位年齡在22~29歲之間的受試者進行試驗，利用所得的數據訓練模型參數。通過仿真關節角度結果和實際測量數據對比以及求解預測COP與真實記錄COP的相關系數，證實了模型的合理性與實用性。
相對與之前研究者通過簡化身體結構和假設人腦控制器的人體模型建立方法，本文所提出的方法具有明顯優點：首先，方法簡單易行，同時具有更高的準確性，能充分反映人體結構的協調關系；其次，可以用同樣方法辨識多方向的運動模型，人體運動無需局限在一個二維平面內，避免了大量人體運動信息的丟失；最后，由于辨識的輸入是實際人體運動數據，所建立的模型可以準確地反映人類個體之間的差異性以及患者自身在不同訓練階段的康復效果。
引言
站立平衡功能關系到人們的生存質量，近年來關于人體站立平衡功能的康復設備和訓練方法逐漸成為研究的熱點[1-2]。本課題試圖開發一種被動平衡康復訓練與評測系統，該系統由多維運動平臺提供外力，干擾平臺上訓練者的站立平衡狀態，使訓練者被動地鍛煉自身平衡調節機能[3]。這種康復訓練設備存在的問題是外力的刺激強度難以選擇，如果外力強度過低則難以達到康復治療訓練的效果，而外力強度過大又會使訓練者摔倒或者造成關節組織扭傷，此外患者的身體條件不盡相同也導致系統安全強度需因人而異進行調整。在康復治療必須保證患者安全這一前提下，為了進一步完善該被動訓練系統的訓練方法和訓練策略，需要建立一種可體現測試者個體差異的人體站立平衡調節數學模型，仿真人體主要關節角度在平臺運動下的動態響應過程，最終預測平臺運動對人體刺激的安全強度范圍。
一些研究者提出了基于工程模型的建模方法，將人體的肌肉骨骼簡化為倒立擺，并為其建立控制器和控制算法，模擬人類中樞神經系統(center neural system,CNS)的功能[4]。還有研究者利用一些實驗或醫學理論證明了CNS在功能上可近似等效為比例-積分-微分控制器 (proportion-integration-differentiation controller，PID)、比例-微分控制器(proportion-differentiation controller，PD)、最優控制器 (Optimal controller)或間歇控制器(Intermittent controller)等[5-7]。這類工程模型雖然可以近似地為研究者提供模擬人類運動控制的數學依據，但仿真結果不能夠反映不同人的運動特點，難以應用于工業產品開發或醫療產品設計。
本研究采用T-S模糊辨識方法，基于9位健康成年人的真實測試數據建立數學模型，以平臺的運動加速度數據為輸入信號，人體的主要關節角度動態數據為系統的輸出信號。該方法不需要考慮人體結構的耦合關系，僅將人體的中樞神經系統和人體骨骼肌肉系統視為內部結構未知的黑箱。對T-S模糊模型的前件參數的辨識，采用的是模糊C-均值聚類(fuzzy C-means，FCM)算法，但由于FCM算法對聚類中心初值要求苛刻，聚類結果易陷入局部最小值，導致模型的辨識精度和效率降低[8]，所以本研究采用人工蜂群尋優(artificial bee colony，ABC)算法對聚類原型的初值進行優化。由于ABC具有較好的全局優化能力[9]，克服了聚類算法的缺陷，因此保證了T-S模型前件參數的正確性。
最后本文利用遞歸最小二乘法，求解模型的后件參數，并推導出模型的解析函數。本模型可以有效表征不同受試者個體間的運動特點和差異，并通過計算仿真數據和測量數據的均方差(mean square error，MSE)以及互相關度，證明所建模型的準確性和客觀性。
1 模型辨識
1.1 T-S模糊模型
T-S模糊模型可以表示復雜的系統，以很高的精度逼近定義在緊集上的非線性函數，并且模型的結果用線性方程表示[8]，適合于人體站立平衡系統建模。辨識T-S模糊模型的前件參數采用FCM聚類算法，將數據按照一定的相似性準則劃分成多個子集，以表征系統的不同性征，并對集群的中心進行編碼。因此，基于FCM的T-S模糊模型，可以表示為
| $\begin{align} & {{R}_{i}}:\text{ }IF\text{ }x\text{ }is\text{ }in\text{ }cluster~{{C}_{i}}~with\text{ }({{v}_{i}},{{u}_{i}}), \\ & THEN~{{y}^{i}}={{a}_{0}}^{i}+{{a}_{1}}^{i}{{x}_{1}}+{{a}_{2}}^{i}{{x}_{2}}+\ldots +{{a}_{d}}^{i}{{x}_{d}}, \\ \end{align}$ |
其中 i =1,2,…,R代表聚類中心的序號；向量x=(x1,x2,…,xd)為給定的輸入變量；Ci為模糊變量產生的模糊聚類空間；d為輸入變量維數；θi=(a0i,a1i,a2i,…,adi)代表結論部分的輸出參數；yi為第i個模糊規則的輸出結果。模糊模型的輸出Y是一個加權平均數去模糊化計算方法，可表示為
| $Y=\frac{\sum\limits_{i=1}^{R}{{{w}_{i}}{{y}^{i}}}}{\sum\limits_{i=1}^{R}{{{w}_{i}}~}},$ |
其中wi為第i個模糊規則輸出結果的權值，即
| ${{w}_{i}}=exp-{{\left( \frac{\|x-{{v}_{i}}\|}{{{\sigma }_{i}}} \right)}^{2}}~,$ |
其中σi為常數，取σi=2；vi是輸入變量x的聚類中心坐標位置，即
| ${{v}_{i}}=\frac{\sum\limits_{i=1}^{R}{{{u}_{ij}}{{x}_{i}}}}{\sum\limits_{i=1}^{R}{{{u}_{ij}}~}}~~,$ |
| ${{u}_{ij}}=\frac{1}{\sum\limits_{k=1}^{R}{\left( \frac{\|{{x}_{i}}-{{v}_{j}}{{\|}^{2}}}{\|{{x}_{i}}-{{v}_{k}}{{\|}^{2}}} \right)}}$ |
其中uij為隸屬度值，下標j=1,2,…,N表示數據樣本序號，k=1,2,…,R，且k≠i。
1.2 人工蜂群改進的FCM聚類算法
FCM通過對目標函數的迭代優化來取得對數據集的模糊聚類，目標函數在迭代過程中是遞減的。這種方法在很大程度上依賴于初始聚類中心的選擇，如果初始聚類中心矩陣選擇不合理，會造成錯誤聚類，還會增加系統復雜度，降低算法的效率。
ABC是模擬蜜蜂采蜜方式的仿生式人工智能優化算法，具有引領蜂、跟隨蜂和偵查蜂為三種不同功能的尋優因子。每個食物源的坐標對應樣本的一個可行解(聚類中心)，食物源的收益度對應問題的適應度，最大收益度對應問題的最優效果。這里的適應度用于評價各個聚類中心在優化過程中的價值，評價這個可行解的質量，作為是否需要保留的條件，即
| $fi{{t}_{n}}=\frac{1}{1+objec{{t}_{n}}},$ |
其中objectn是要解決優化問題的目標函數，由具體問題而定。ABC計算步驟如下：
(1)初始化參數：設定尋優問題的食物源個數N，最大循環次數100，最大迭代次數200，循環次數m和迭代次數n的初值均為0。
(2)產生初始解集xij，并且通過式(6)計算初始解集的適應度。
(3)引領蜂通過式(7)做鄰域搜索，產生新解vij，即
| ${{v}_{ij}}={{x}_{ij}}+{{\Phi }_{ij}}({{x}_{ij}}-{{x}_{kj}})\text{ },$ |
并且計算其適應度，如果比上一次xij的適應度大，則用vij取代xij，否則保留原xij；
(4)利用新解xij的適應度，根據式(8)計算概率值Pij，用以選擇食物源，即
| ${{P}_{ij}}=\frac{fi{{t}_{n}}}{\sum\limits_{n=1}^{200}{fi{{t}_{n}}}}$ |
(5)鄰域搜索，產生新解vij，計算其適應度，若大于xij的適應度，則用vij取代xij，否則xij保留。
(6)鄰域搜索計數到達最大循環次數后，判斷是否有要丟掉的解，如果存在，偵查蜂根據式(9)產生一個新解代替原xij。
| ${{x}_{ij}}=xmin+rand\left( 0,1 \right)\times (xmax-{{x}_{min}})$ |
(7)判斷當前的迭代次數是否已經到達最大迭代次數，如果到達則停止迭代，輸出最后一代的食物源作為優化問題的最優解集合，否則轉到步驟(3)，迭代計數加一。
ABC改進FCM算法(ABC-FCM)算法的思想可總結為：利用ABC算法良好的全局搜索能力求得最優解，作為FCM算法的初始聚類中心，再利用FCM算法求得最終聚類中心解。ABC優化目標函數和FCM聚類算法的聚類目標函數均選為式(10)，尋優與聚類使objectn最小。
| $objec{{t}_{n}}=\sum\limits_{i=1}^{N}{\sum\limits_{j=1}^{C}{{{u}_{ij}}\|{{x}_{i}}-{{x}_{j}}{{\|}^{2}}}}$ |
ABC-FCM算法即在原ABC優化算法的基礎上加入以下步驟：
(8)根據式(5)計算新隸屬度函數U={uij|j=1,2,…,N; i=1,2,…,R}。
(9)根據式(4)計算新聚類中心V={vi|i=1,2,…,R}。計算相鄰兩代隸屬度之差E，若E＜ε，則尋優停止，否則轉到步驟(8)。
1.3 T-S模型辨識算法
所研究的人體站立平衡系統選取了平臺在水平面的加速度數據為輸入，人體踝、膝、髖關節的角度變換為輸出，組成了二維輸入三維輸出問題。為了便于模型的辨識，本文采取的方法是將三個關節獨立建模，即將原問題視為分別針對踝、膝、髖關節的三個二維輸入一維輸出問題，同時為了進一步提高辨識準確性，將另外兩關節過去時刻的角度數據作為當前辨識計算的輸入，這種方法在實際應用中也易于實現。以辨識踝關節為例介紹具體算法：
(1)T-S模型參數訓練
輸入： 令Z=(A前后(k)，A左右(k)，φH(k-1),φK(k-1),令y=φA(k)，A前后(k)，A左右(k)是運動平臺所在水平面的前后和左右方向的加速度數據；φH(k-1),φK(k-1)人體踝、膝關節前一時刻的角度數據。
輸出： T-S模型前后件參數：{U,V,θ}。
第1步： 設定聚類數R，按照ABC-FCM聚類方法，利用式(4)、(5)，分別求解聚類中心V和隸屬度函數矩陣U。
第2步： 用最小二乘算法求解后件參數。通過目標函數最小化：
| $J=min\frac{1}{N}\left\{ {{(y-Xe\theta e)}^{T}}\Phi i(y-Xe{{\theta }_{e}}) \right\}$ |
其中Xe=[1,X]，Φi=diag(ui1,ui2…um)。因此，模型在本條規則的輸出的參數可以得到：
| $\theta i={{({{x}_{e}}^{T}\Phi i{{x}_{e}})}^{-1}}{{x}_{e}}^{T}{{\Phi }_{iy}}$ |
(2) T-S模型輸出方程
輸入： ① T-S模型結構參數{U,V,θ}；② 受試者在平臺同方向以3 cm幅度運動所獲得的數據樣本Z=(A前后(k)，A左右(k),φH(k-1),φk(k-1))；③ 實際輸出y。
輸出： 預測輸出函數Y=φA(k)和MSE。
第1步： 計算樣本Z和聚類中心向量V的距離。
第2步： 根據式(3)，計算樣本Z對每個本地輸出函數yi的隸屬度，這個隸屬度值wi可以看作這個規則下輸出yi的權值。
第3步： 計算每個規則的輸出相應值yi，利用當前模型的參數θi，根據式(2)可以計算得到輸入為Z的估計輸出Y。
第4步： 利用均方差式(13)，計算該仿真結果與真實數據之間的MSE。
| $MSE=\frac{1}{N}\sum\limits_{j=1}^{N}{{{({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{y}}}_{j}}-{{y}_{j}})}^{2}}}$ |
2 實驗設計
T-S模糊辨識方法需要真實人體數據對模型結構和參數進行訓練，另一方面需要數據對辨識結果進行驗證[8]，所以本節詳細介紹試驗設計與數據采集方法。
2.1 受試者和測試設備
受試者包括9名健康成年人(4女，5男)，年齡(25.2±2.1)歲；身高(169.1±4.9 )cm；質量(62.6±9.7)kg。受試者沒有肌肉或神經運動障礙病史，并且進行了相關臨床檢測[10]。
所用到的測試設備包括：多維運動平臺。四臺高分辨率攝像機。計算機，用于視頻信號記錄和圖像處理。基于ADIS16355陀螺儀傳感器的平臺運動測量裝置，該裝置固定于平臺上，用于采集平臺的運動加速度數值。人體重心壓力坐標采集板[11]。
2.2 實驗方法
測量受試者的身體信息：將人體以踝、膝、髖關節的中心點分為上中下三個體段，A表示踝關節；K表示膝關節；H表示髖關節；L表示下段(小腿)；M表示中段(大腿)；U表示上段(上段包括軀干、手臂和頭部)[11-12]。mL、mM、mU分別指三體段的質量；cL、cM、cU分別指三體段質心到轉動關節的距離；lL、lM、lU分別是三體段的長度；參數Mwhole人體的總質量[13]。
受試者赤腳站立在運動平臺上，雙手自然下垂，放于身體兩側，兩腳并攏，如圖 1所示。測試時要求受試者閉上雙眼，不可抬腳、擺動手臂，僅依靠身體的主要關節來調整平衡姿態。平臺的水平運動數據通過陀螺儀加速度傳感器模塊采集，記錄為A前后(k)，A左右(k)。人體重心壓力坐標采集板將每次試驗過程中人體壓力中心(center of pressure,COP)的坐標記錄下來，作為模型仿真結果的驗證依據。
圖1
受試者的站立位置和姿態
Figure1.
Position and posture of the testee
通過文獻[3]采用的線性分析試驗可知，運動平臺的擾動形式為“水平面內的脈沖式平動”，幅度為2~3 cm，時間帶寬約為200 ms。本文所說外力強度可對應理解為平臺單次脈沖運動的幅度，該擾動幅度不但使人體晃動數據具有較高的信噪比，也保證了受試者的安全。本文將2 cm幅度擾動條件下的人體關節角度動態數據用于T-S模型系統辨識參數訓練；將3 cm幅度擾動條件下的人體關節角度動態數據作為辨識獲得模型預測值的驗證數據。由于T-S模型不需要考慮人體的結構耦合關系，所以平臺在任意方向的運動均可以利用該方法建模。用于參數訓練的人體運動數據越全面，T-S模型的預測精度就越高，但是在臨床使用中，不宜對患者進行過久的測試，可能會造成模型參數訓練不足，降低了模型的準確性。為了盡可能使樣本具有遍歷性，平臺以水平面內的脈沖式多方向運動，各方向間隔45°，如圖 2圓盤內箭頭所示，分別為：前(anterior,A)，后(posterior,P)，左(left,L)，右(right,R)，左前(anterior-left,AL)，左后(posterior-left,PL)，右前(anterior-right,AR)，右后(posterior-right,PR)。
圖2
平臺運動方向和攝像機位置
Figure2.
Platform moving direction and camera position
攝像機與受試者髖關節中心點的高度一致。4臺攝像機的間隔角度為45°，朝向平臺的中心，如圖 2所示。均在人體朝向攝像機方向的肩部、髖關節、膝關節和踝關節貼有4組標識點，用于身體運動位置的視頻提取。平臺每次運動方向所在矢量平面與人體的晃動平面相同，并且總有一臺攝像機的鏡頭方向與人體的晃動平面保持垂直，使對圖像的觀測誤差盡量小。測量開始后，平臺的運動方向按照之前規定的激勵方式進行，給定的平臺運動激勵間隔為10 s。測試過程中，受試者可能會習慣平臺激勵，進而做出提前應激動作，這將破壞數據的客觀性，為此平臺的運動方向順序是隨機的。待8個方向進行完畢后，按同樣方法連續測試三次。
通過圖像處理技術，提取所拍攝畫面中標記點的坐標。將人體簡化為矢量平面內多連桿體段倒立擺，將每幀畫面關節標識點的位置與初始點位置比較，可以得到人體運動狀態。以求解踝關節角度變換數據為例：圖 3為2號攝像機得到的人體運動下肢部分簡圖。2號攝像機正對受試者右側，圖 3中關節處白色圓形為提取得到的標識點，定義體段與鉛垂線的夾角為其相鄰轉動關節的關節角，將相鄰關節標識點的水平運動距離相減可得運動距離stA-st0A-sP=ΔsA，stK-st0K-sP=ΔsK，ΔsAK=ΔsA-ΔsK，其中st0A、stA表示踝關節在單位時間開始和結束的橫坐標，sP表示單位時間內平臺運動的距離，st0K、stK表示膝關節在單位時間開始和結束的橫坐標。踝關節運動角度數據用反三角正弦函數可得，即
| ${{\theta }_{t}}^{A}=\frac{arcsin\left( \frac{\Delta {{s}^{AK}}}{{{l}_{L}}} \right)}{t-{{t}_{0}}}$ |
圖3
踝關節和膝關節標致點運動簡圖
Figure3.
Motion diagram of the mark points on ankle and knee
膝關節和髖關節可用同樣方法獲得動態角度數據集，關節數據分別記做踝：(φA)、膝(φK)、髖(φH)。該計算方法已在Matlab7.0軟件環境編程，可以批量處理試驗獲得的人體數據。
2.3 實驗步驟
(1)試驗前向受試者說明試驗任務，確保受試者正確理解試驗操作并進行熟悉性練習后，簽訂試驗協議。
(2)受試者填寫個人基本資料，測量身高、體段及體重等基本信息后，在身體規定位置，安裝標識點。
(3)受試者按要求于運動平臺上站好。確認無誤并提示受試者后，啟動平臺，開始施加運動刺激。
(4)平臺運動20 s，確認受試者無不良反映，開啟采集設備電源，進行數據采集；若受試者出現不適，應立即停止平臺運動刺激。
(5)按要求方法采集完畢后結束試驗，關閉數據采集系統，輸入平臺停止運動指令，待平臺停止運動過后，令受試者離開平臺。
(6)完成基本數據的計算，編號后存入硬盤。
3 結果與分析
3.1 聚類算法結果與分析
通過對試驗獲得的實際人體角度數據矩陣X=[φA,φK,φH]進行聚類計算，結果如圖 4所示。
圖4
聚類計算結果
Figure4.
Clustering results of calculation
如圖 4(上)所示，*表示聚類的輸入數據集合；△表示ABC優化的聚類初始點；○表示FCM聚類的最終聚類中心點。通過ABC算法的改進，聚類中心分布比較均勻，FCM算法在優化的聚類初始點基礎上進行微小調整即達到最終聚類結果。聚類目標函數的迭代結果比較如圖 4所示，ABC-FCM的初始迭代效果明顯優于FCM，聚類目標在迭代前100次便達到較低的數值，之后趨于平穩。程序設定ABC最大迭代次數為200次，圖 4(下)中實線為ABC-FCM聚類，200次ABC優化迭代后轉為FCM繼續聚類目標函數，結果穩定在2 500左右。實曲線200~250次FCM迭代對應于圖 4(上)中的變化，即FCM算法在ABC聚類初始點基礎上進行的微小調整。虛線為單純利用FCM聚類算法得到的目標函數值，聚類效果較差，250次迭代運算后結果處于3 000左右，從該結果可以看到ABC-FCM的聚類效果更加合理。由于T-S模型辨識的前件參數辨識對準確性有較大的要求，這種改進對模型建立的準確性有較大改善。聚類程序運算過程中，ABC-FCM聚類相對于FCM聚類占用了更多的系統時間，在同樣的計算機和Matlab7.0仿真環境中，ABC-FCM聚類計算占用2.321 s，而FCM聚類計算僅占用0.283 s。可見ABC優化算法對聚類性能的改善是以消耗系統時間為代價的，但是考慮到本課題選用的T-S模糊模型辨識屬于離線辨識方式，并不要求較高的系統辨識速度，因此該問題對于本課題的影響并不明顯。
3.2 辨識結果與分析
隨機選取一位受試者在平臺后向運動過程中記錄的試驗數據建立模型，并且對比仿真運動結果，踝、髖、膝三個關節的實際運動過程與對應關節仿真的結果如圖 5所示。
圖5
關節數據與仿真結果對比
Figure5.
Comparison of joint data with simulation results
如圖 5所示，橫坐標為測試數據樣本序號，該樣本的持續時間從平臺運動開始到1500 ms為止。仿真結果與真實數據之間的MSE分別為：踝關節的MSE=0.017；髖關節的MSE=0.018；膝關節的MSE=0.220。仿真結果的誤差處于較小的范圍，三個關節的運動幅度可以被如實預測，同時各個關節之間的運動協調關系基本與真實人體數據一致。基于同樣的方法，可以求得該受試者在平臺刺激為其它方向時的人體數學模型，并且仿真出平衡調節運動數據；通過計算這些仿真數據和真實測試數據的MSE，可知全部仿真結果的MSE＜0.271。通過式(15)可以計算相應的壓力中心點COP坐標變化曲線[11-12]，
| $\begin{align} & CO{{P}_{前后}}=-[({{m}_{L}}{{c}_{L}}+{{m}_{M}}{{l}_{M}}+{{m}_{U}}lU)\cdot sin({{\varphi }_{A}})+ \\ & ({{m}_{M}}{{c}_{M}}+{{m}_{U}}{{l}_{U}})\cdot sin({{\varphi }_{A}}+{{\varphi }_{K}})+ \\ & {{m}_{U}}{{c}_{U}}\cdot sin({{\varphi }_{A}}+{{\varphi }_{K}}+{{\varphi }_{H}})]/{{M}_{whole}} \\ \end{align}$ |
計算仿真數據和采集實際重心坐標數據的互相關程度。通過對9位受試者的數據仿真，可以得到所有受試者的仿真結果與實際采集得到COP數據的相關度，全部結果在0.641~0.785之間。一般認為，當COP坐標移動到足底與平臺接觸面之外，人體就會失去平衡，現以不同受試者的足底長度為其前向判斷摔倒閾值，用所建立模型預測數據進行判斷，可獲得不同受試者的安全刺激強度范圍。9位受試者安全刺激強度的預測范圍處于5~7.5 cm之間。
需要指出，求解COP的過程存在三種誤差來源：① 本研究所用T-S模糊模型辨識方法受到迭代次數的限制，預測結果不可能零誤差逼近實際測量數據；② 式(15)涉及諸多肢體參數，測量時采用了文獻[13]介紹的方法，由于人為操作的問題同樣會引入測量誤差；③ 人體密度分布不均勻，形狀不規則，各體段質心的位置以現有的測量手段難以獲得，而本研究采用了基于統計數據的比例法，并且將各個體段等效為密度均勻的圓柱體。上述三方面誤差累積，會在一定程度上降低COP結果的準確性，但是較以往的人體運動模型，該方法所建立的人體站立平衡模型能夠較準確地預測人體運動狀態，并且仿真結果足以體現個體間的差異，可針對不同患者制定安全刺激強度。
4 結論
本研究針對目前人體平衡系統建模方法存在的不足，通過T-S模糊模型辨識方法，建立了平臺運動加速度和人體踝、髖、膝三個關節運動的數學模型，并且用人工蜂群算法改進了模型前件參數辨識的運算效率和聚類精度，對9位年齡在22~29歲之間的受試者進行試驗，利用所得的數據訓練模型參數。通過仿真關節角度結果和實際測量數據對比以及求解預測COP與真實記錄COP的相關系數，證實了模型的合理性與實用性。
相對與之前研究者通過簡化身體結構和假設人腦控制器的人體模型建立方法，本文所提出的方法具有明顯優點：首先，方法簡單易行，同時具有更高的準確性，能充分反映人體結構的協調關系；其次，可以用同樣方法辨識多方向的運動模型，人體運動無需局限在一個二維平面內，避免了大量人體運動信息的丟失；最后，由于辨識的輸入是實際人體運動數據，所建立的模型可以準確地反映人類個體之間的差異性以及患者自身在不同訓練階段的康復效果。
