基于ESN的下肢表面肌電信號(hào)預(yù)測(cè)算法

熊 鳴

(北京信息科技大學(xué) 自動(dòng)化學(xué)院，北京 100192)

0 引言

隨著殘疾人和老年人口不斷增加，依靠科技創(chuàng)新來保障及改善他們的健康已經(jīng)成為當(dāng)前許多國家的戰(zhàn)略需求。研發(fā)先進(jìn)的功能康復(fù)機(jī)器人，對(duì)老年健康服務(wù)與助殘公益事業(yè)具有十分重要的社會(huì)意義[1]。相比于專業(yè)康復(fù)師，康復(fù)機(jī)器人有執(zhí)行效果始終保持一致的優(yōu)點(diǎn)，缺點(diǎn)是無法根據(jù)訓(xùn)練者的實(shí)際情況更改訓(xùn)練計(jì)劃。人的中樞神經(jīng)系統(tǒng)中運(yùn)動(dòng)神經(jīng)發(fā)出的動(dòng)作電位通過神經(jīng)纖維傳導(dǎo)到肌肉，肌肉因此產(chǎn)生伸縮，肌電信號(hào)(electromyography，EMG)就是由許多微弱動(dòng)作電位組合的序列總和，其中蘊(yùn)涵了很多與人體動(dòng)作有關(guān)聯(lián)的信息[2]。如果康復(fù)機(jī)器人能夠通過傳感器直接獲知訓(xùn)練者肌電信號(hào)，就可以適時(shí)地對(duì)訓(xùn)練計(jì)劃做出改變以滿足訓(xùn)練者要求。

Borui Li[3]設(shè)計(jì)了可穿戴下肢和踝關(guān)節(jié)EMG裝置，利用EMG信號(hào)來提示使用者從坐到站的趨勢(shì)，誤判率和漏判率都低于5%。Jimson Ngeo等[4-5]設(shè)計(jì)了輔助運(yùn)動(dòng)用的外骨骼裝置，通過采集健康手指的EMG信號(hào)，預(yù)測(cè)出健康手指的運(yùn)動(dòng)趨勢(shì)以及力量，并控制輔助殘肢的外骨骼，取得了較為理想的效果。Tanvir等[6]采集膝關(guān)節(jié)的EMG信號(hào)，并利用自適應(yīng)神經(jīng)模糊網(wǎng)絡(luò)推理系統(tǒng)來估算膝蓋的彎曲角度。傳統(tǒng)的EMG信號(hào)處理都采用了神經(jīng)網(wǎng)絡(luò)來進(jìn)行，需要確定神經(jīng)網(wǎng)絡(luò)的規(guī)模；而網(wǎng)絡(luò)規(guī)模的選擇只能依靠研究者的經(jīng)驗(yàn)，數(shù)據(jù)預(yù)測(cè)結(jié)果存在著不確定性。Wang等[7]利用二進(jìn)制粒子群優(yōu)化方法 (binary particle swarm optimization,BPSO) 優(yōu)化回聲狀態(tài)網(wǎng)絡(luò)(echo state network，ESN)的輸出層連接,提出了一種利用改進(jìn)PSO優(yōu)化ESN網(wǎng)絡(luò)參數(shù)，并用優(yōu)化后的ESN對(duì)EMG信號(hào)進(jìn)行建模以及預(yù)測(cè)的方法。張杰爍等[8]用基于遞歸最小二乘法的回聲狀態(tài)網(wǎng)絡(luò)心電信號(hào)降噪算法來消除心電信號(hào)中的復(fù)雜噪聲。該算法得到含噪心電數(shù)據(jù)中非線性的且具有區(qū)分度的深層次特征，并利用這些特征自動(dòng)分離心電信號(hào)與噪聲。王磊等[9]提出基于靈敏度分析的模塊化回聲狀態(tài)網(wǎng)絡(luò)修剪算法，利用靈敏度大小判斷子儲(chǔ)備池模塊的貢獻(xiàn)度,刪除靈敏度低的子模塊，修剪之后的網(wǎng)絡(luò)具有較好的泛化能力和魯棒性。

回聲狀態(tài)網(wǎng)絡(luò)[10]是一類新型回歸神經(jīng)網(wǎng)絡(luò)，在混沌和非線性動(dòng)態(tài)系統(tǒng)建模、辨識(shí)和控制等方面取得了成功應(yīng)用。如在經(jīng)典的Mackey-Glass(MG)混沌時(shí)間序列的前向第84步預(yù)測(cè)，其預(yù)測(cè)精度相比常見的方法提高了約2個(gè)數(shù)量級(jí)。本文采用回聲狀態(tài)網(wǎng)絡(luò)來對(duì)EMG信號(hào)進(jìn)行建模以及預(yù)測(cè)。

1 EMG檢測(cè)裝置

EMG檢測(cè)裝置主要由電極、檢測(cè)電路和信號(hào)調(diào)理電路組成。檢測(cè)下肢股二頭肌EMG信號(hào)的裝置如圖1所示。

其中電極1與電極2分別貼在大腿肌后群股二頭肌上，保持一定間距，保證兩個(gè)電極在大腿做屈伸運(yùn)動(dòng)時(shí)，能夠產(chǎn)生電壓差。由于檢測(cè)裝置后端采用精密整流電路，故兩個(gè)電極相對(duì)位置不影響信號(hào)輸出特性。檢測(cè)電路采用差分放大電路，差分放大器采用AD8221，該芯片具有高增益和高共模抑制比特性，并且增益可調(diào)。

當(dāng)大腿進(jìn)行屈伸運(yùn)動(dòng)時(shí)，所產(chǎn)生的差分信號(hào)經(jīng)過放大器之后成為交流信號(hào)，經(jīng)整流后進(jìn)入濾波電路。濾波電路由工頻陷波電路組成以濾除工頻干擾。為了降低EMG檢測(cè)裝置輸出阻抗，在濾波電路之后采用同向比例放大電路。

本文檢測(cè)股二頭肌EMG信號(hào)，故此EMG參考電位為遠(yuǎn)離大腿肌后群的股四頭肌，相當(dāng)于電極1和2在大腿的后部，參考電極在大腿的前部。參考電極與差分放大電路參考地線相連。電極與檢測(cè)裝置之間采用導(dǎo)聯(lián)線進(jìn)行連接。

2 回聲狀態(tài)網(wǎng)絡(luò)

ESN由被稱為儲(chǔ)備池的內(nèi)部網(wǎng)絡(luò)與簡單的輸入輸出網(wǎng)絡(luò)組成，典型ESN拓?fù)浣Y(jié)構(gòu)如圖2所示。圖中左側(cè)為M個(gè)輸入節(jié)點(diǎn)；中間為漏積分神經(jīng)元組成的儲(chǔ)備池，由稀疏矩陣以及稀疏矩陣所對(duì)應(yīng)的連接權(quán)值組成，其中稀疏矩陣為一L×L的矩陣；圖中右側(cè)為N個(gè)輸出節(jié)點(diǎn)。其中連接輸入權(quán)重矩陣記為Win，儲(chǔ)備池稀疏矩陣記為W，輸出連接權(quán)重矩陣記為Wout，輸出反饋連接權(quán)重矩陣記為Wback。當(dāng)ESN作為波形發(fā)生器和模式發(fā)生器時(shí)，Wback是必需的矩陣。對(duì)于本文來說，ESN作為系統(tǒng)建模使用，故此Wback在本文中默認(rèn)為零矩陣。對(duì)ESN來說，Wout是通過訓(xùn)練得到的矩陣。

回聲狀態(tài)網(wǎng)絡(luò)可表示為

x(t+1)=(1-α)×x(t)+α×f(Win×

u(t)+W×x(t))

(1)

(2)

通過式(2)可以得到Wout的最優(yōu)估計(jì)。實(shí)際計(jì)算中，由于xTx比較容易接近奇異，將會(huì)導(dǎo)致Wout病態(tài)解，因此常利用嶺回歸的方法來替換式(2)的偽逆方法：

(3)

式中λ為正則項(xiàng)系數(shù)，用其作為噪聲改善xTx可能的奇異性，以保證Wout有合適解。由于x初始狀態(tài)為零矩陣，為了消除零矩陣初始狀態(tài)的影響，參與式(3)計(jì)算中的x從n1+1時(shí)刻后開始參與計(jì)算。即[0，n1]為序列初始化，這時(shí)間段內(nèi)x不參與Wout計(jì)算。ESN訓(xùn)練序列長度為n2，其中[0，n1]?[0，n2]。

3 ESN參數(shù)優(yōu)化方法

3.1 參數(shù)優(yōu)化目標(biāo)

對(duì)于ESN來說如何獲得一個(gè)好的儲(chǔ)備池是一個(gè)至關(guān)重要的問題[10]。雖然相關(guān)研究很多，但是沒有一個(gè)較為系統(tǒng)全面的方法。常用方法是利用訓(xùn)練殘差來調(diào)整儲(chǔ)備池權(quán)值，并進(jìn)行相應(yīng)的迭代運(yùn)算以得到最優(yōu)參數(shù)。源于鳥類捕食行為模擬的PSO沒有變異以及交叉環(huán)節(jié)，容易陷入?yún)^(qū)域最優(yōu)解。ESN參數(shù)有4個(gè)：1)譜半徑ρ；2)儲(chǔ)備池規(guī)模L；3)儲(chǔ)備池尺度因子α；4)儲(chǔ)備池稀疏度。儲(chǔ)備池規(guī)模與擬合數(shù)據(jù)類型有關(guān)系，在本文中作為固定參數(shù)。儲(chǔ)備池稀疏度調(diào)整范圍較小，也作為固定參數(shù)使用。將譜半徑以及尺度因子作為待優(yōu)化參數(shù)，利用PSO對(duì)這兩個(gè)參數(shù)所組成的二維空間進(jìn)行搜索。PSO算法實(shí)質(zhì)是在空間搜尋全局最優(yōu)解[12]，對(duì)搜索的過程和方向進(jìn)行了一定約束。如果能夠在搜索方向上有較為明確規(guī)定，則PSO計(jì)算量會(huì)大為減小，可以縮短優(yōu)化時(shí)間。

某個(gè)信號(hào)序列ESN最優(yōu)輸出連接權(quán)重問題可以轉(zhuǎn)化為一個(gè)約束優(yōu)化問題。相對(duì)誤差函數(shù)(RES)定義為

(4)

利用PSO來求解式(4)最優(yōu)化問題。式(4)中所使用序列為[n3+1，n2]，該序列稱為測(cè)試序列，測(cè)試序列長度小于訓(xùn)練序列長度。

式(4)是對(duì)傳統(tǒng)均方根誤差(RMSE)的改進(jìn)。RMSE定義如下：

R=

(5)

RMSE是時(shí)間序列預(yù)測(cè)等回歸分析中常用的評(píng)價(jià)指標(biāo)，著重點(diǎn)在于絕對(duì)誤差，式(4)對(duì)RMSE進(jìn)行了修改。相對(duì)誤差越小則認(rèn)為預(yù)測(cè)輸出向量與期望輸出向量之間歐式距離越小，兩條曲線擬合度越高。

3.2 PSO粒子位置和速度的更新機(jī)制

PSO算法采用速度-位置模型[11]。PSO算法關(guān)鍵在于粒子位置更新，位置更新依靠更新速度。其數(shù)學(xué)描述如下。

設(shè)在一個(gè)N維解空間內(nèi)，由M個(gè)粒子組成一個(gè)群組，第i個(gè)粒子位置和速度可表示為：pid=(pi1,pi2,…,pin,…,piN)和vid=(vi1,vi2,…,vin,…,viN)，其中i=1,2,…，M；第i個(gè)粒子區(qū)域最優(yōu)值為Pbestid=(xi1best,xi2best,…,xinbest,…,xiNbest)，全局最優(yōu)值為Gbestd=(x1best,x2best,…,xnbest,…,xNbest)。

粒子位置更新公式為

(6)

(7)

式中：I為總迭代次數(shù)；i為當(dāng)前迭代次數(shù)；ωmax=1和ωmin=0.1分別為慣性權(quán)值的初始值和終值。

為了加快PSO優(yōu)化速度，對(duì)PSO更新速度方法進(jìn)行了改進(jìn)，在基本PSO更新速度基礎(chǔ)上增加一個(gè)約束條件。處理方式如式(8)所示：

pgbest2=(M+1-j)/M×pgbest+

(j-1)/M×psfbest

(8)

式中：pgbest2為當(dāng)前迭代最優(yōu)解位置；pgbest為RES當(dāng)前迭代最優(yōu)解位置；psfbest為絕對(duì)誤差函數(shù)(absolute error，AE)當(dāng)前迭代最優(yōu)解位置；M為粒子群組的個(gè)數(shù)；j為AE最優(yōu)解所對(duì)應(yīng)位置在本次RES序列排序中的序號(hào)。AE定義如下：

(9)

從式(4)和式(9)可以看出，RES和AE最優(yōu)解位置不可能完全相同，存在最優(yōu)解不一致的可能性。假設(shè)AE和RES最優(yōu)解位置相同，此時(shí)j=1，則pgbest2與pgbest相等，即當(dāng)前迭代最優(yōu)解位置與RES最優(yōu)解位置相同。假設(shè)AE最優(yōu)解位置與RES次最優(yōu)解位置相同，則j=2，以此類推。式(8)對(duì)兩種目標(biāo)函數(shù)所得到的位置進(jìn)行線性組合，以得到當(dāng)前迭代最優(yōu)解位置。從式(8)可以看出最優(yōu)解位置還是以RES最優(yōu)位置為主體，以保證ESN輸出向量與期望輸出向量有更好的擬合度。

由于基本PSO易陷入?yún)^(qū)域最優(yōu)解，故借鑒遺傳算法變異和交叉操作，對(duì)每次迭代中部分粒子位置進(jìn)行變異和交叉工作。每次迭代中，對(duì)所得到的RES數(shù)值從小到大排序，對(duì)于排序中最后30%的粒子進(jìn)行該操作。將最后10%的粒子群組稱為A組，除去這10%之后剩下的20%的粒子群稱為B組。即A與B的并集組成了排序中后30%的粒子群。A組進(jìn)行變異的規(guī)則為：

1)當(dāng)A組RES數(shù)值遠(yuǎn)大于B組數(shù)值；

2)當(dāng)A組RES數(shù)值與B組后半部分RES數(shù)值相同時(shí)。

兩個(gè)條件滿足任一則進(jìn)行變異。當(dāng)不滿足上述條件時(shí)，進(jìn)行交叉操作。交叉操作對(duì)象為一個(gè)隨機(jī)生成位置。本文中算法采用十進(jìn)制編碼方式，交叉操作為將A組粒子位置與隨機(jī)生成位置進(jìn)行加權(quán)平均。

B組進(jìn)行交叉操作，隨機(jī)生成數(shù)如果大于交叉概率0.2時(shí)，進(jìn)行交叉操作。交叉操作為將B組粒子位置與當(dāng)前迭代最優(yōu)解位置進(jìn)行加權(quán)平均。

通過對(duì)A組和B組粒子位置變異與交叉，使得在每次迭代中都有新粒子產(chǎn)生，極大地降低了迭代中陷入?yún)^(qū)域最優(yōu)解的概率。A組變異行為類似于鯰魚效應(yīng)，給下次迭代粒子帶來新位置，使得下一次迭代過程有新的激勵(lì)元素。

3.3 優(yōu)化步驟

本文所使用的優(yōu)化目標(biāo)函數(shù)是RES，故滿足下面兩條中任何一條即可退出優(yōu)化過程：1)當(dāng)A<(n2-n3)×0.052=2.5×10-3×(n2-n3)并且測(cè)試序列中有90%以上的個(gè)體相對(duì)誤差小于5%；2)迭代次數(shù)到達(dá)迭代上限。

優(yōu)化步驟為：1)對(duì)變量進(jìn)行初始化；2)利用式(1)和式(3)得到Wout；3)得到每個(gè)粒子的優(yōu)化目標(biāo)函數(shù)值；4)對(duì)每個(gè)粒子RES和AE進(jìn)行排序，根據(jù)式(8)更新當(dāng)前迭代最優(yōu)解位置；5)對(duì)RES排序中后30%的粒子位置進(jìn)行變異和交叉操作；6)判斷是否滿足退出優(yōu)化條件，不滿足則從第2步開始繼續(xù)優(yōu)化，滿足則退出優(yōu)化過程。

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)結(jié)果

圖3～6為經(jīng)過歸一化處理后的EMG信號(hào)，EMG信號(hào)采樣頻率為280 Hz，S為該采樣數(shù)據(jù)序列的均方差。圖3和圖4為大腿做被動(dòng)屈伸運(yùn)動(dòng)時(shí)EMG信號(hào)。大腿被動(dòng)屈伸運(yùn)動(dòng)主要由外力引起，圖3為左腿EMG信號(hào)，圖4為右腿EMG信號(hào)。圖5和圖6為大腿進(jìn)行自主屈伸時(shí)EMG信號(hào)。大腿屈伸運(yùn)動(dòng)是自主行為，并且該自主行為以右腿蹬踏為主，模擬肢體殘障人士雙腿肌力分布不均的情況。圖5為左腿EMG信號(hào)，圖6為右腿EMG信號(hào)。從圖中可以看出左右腿EMG信號(hào)并不完全相同，主要是由于檢測(cè)裝置以及左右腿肌肉不同所造成。

所選用ESN稀疏矩陣為(500×500)，稀疏度為4%，正則項(xiàng)系數(shù) 為1×10-5，待優(yōu)化參數(shù)為尺度因子和稀疏矩陣譜半徑，訓(xùn)練序列為[0,600]，測(cè)試序列為[500,600]，初始序列為[0,100]。PSO粒子群數(shù)量為10，粒子群速度更新上限為0.1,下限為-0.1。區(qū)域最優(yōu)解學(xué)習(xí)因子以及本次迭代最優(yōu)解學(xué)習(xí)因子均為2，PSO迭代次數(shù)為100。為了驗(yàn)證本文改進(jìn)算法的預(yù)測(cè)效果，對(duì)3種算法進(jìn)行比較:1)利用ESN常用的RMSE作為優(yōu)化目標(biāo)函數(shù)，優(yōu)化采用基本PSO，PSO退出機(jī)制采用本文提出的方法(算法1)。2)利用本文提出的RES作為優(yōu)化目標(biāo)函數(shù)，優(yōu)化采用PSO算法，其中速度-更新算法中慣性權(quán)重為經(jīng)典線性遞減方式如式(7)所示，PSO退出機(jī)制采用本文提出的方法(算法2)。3)利用本文提出的RES作為優(yōu)化目標(biāo)函數(shù)，PSO速度-位置更新辦法采用本文的改進(jìn)算法，退出機(jī)制采用本文提出的方法(算法3)。利用這3種算法對(duì)圖3～6的數(shù)據(jù)進(jìn)行預(yù)測(cè)，結(jié)果如表1～4所示，表中數(shù)據(jù)均為10次預(yù)測(cè)的平均值。

表1 圖3數(shù)據(jù)的10次預(yù)測(cè)平均數(shù)值

表2 圖4數(shù)據(jù)的10次預(yù)測(cè)平均數(shù)值

表3 圖5數(shù)據(jù)的10次預(yù)測(cè)平均數(shù)值

4.2 實(shí)驗(yàn)分析

從表1和表3數(shù)據(jù)對(duì)比來看，3種算法沒有明顯區(qū)別，在不同指標(biāo)中各有優(yōu)勢(shì)。表2和表4數(shù)據(jù)對(duì)比說明，算法3在平均尋優(yōu)速度上具有一定優(yōu)勢(shì)。尤其表4中算法3平均尋優(yōu)時(shí)間是算法2的1/3左右。表2和表4迭代次數(shù)均方差數(shù)據(jù)說明了算法2在10次預(yù)測(cè)中優(yōu)化時(shí)間分布不均衡，因此優(yōu)化時(shí)間最長，這反映了基本PSO在優(yōu)化時(shí)比較容易陷入局部最優(yōu)的特點(diǎn)，同時(shí)也說明本文的改進(jìn)算法能夠在一定程度上避免這種缺陷。

通過數(shù)據(jù)對(duì)比可以看出，當(dāng)數(shù)據(jù)序列均方差較小時(shí)，3種算法效果基本一致，當(dāng)均方差較大時(shí)利用本文算法則可以較快得到符合要求的解。改進(jìn)算法通過對(duì)當(dāng)前全局最優(yōu)解的二次處理以及對(duì)部分粒子進(jìn)行變異和交叉操作之后，提高了解空間的搜索效率。本文算法在各項(xiàng)指標(biāo)中不一定最優(yōu)，但卻是較快找到符合優(yōu)化條件的算法，為預(yù)測(cè)的實(shí)時(shí)性奠定了基礎(chǔ)。

5 結(jié)束語

本文在回聲狀態(tài)網(wǎng)絡(luò)基礎(chǔ)上，使用改進(jìn)PSO算法對(duì)ESN參數(shù)進(jìn)行優(yōu)化，利用PSO找到適合于不同EMG序列的ESN參數(shù)。通過數(shù)據(jù)對(duì)比發(fā)現(xiàn)，ESN在時(shí)間序列預(yù)測(cè)方面具有較好的效果。本文所提出的改進(jìn)PSO算法提高了ESN參數(shù)解空間的搜索效率。