采用小波變換和高斯過(guò)程的肌電信號(hào)模型預(yù)測(cè)

邵輝， 蘇芳茵， 程海波

(華僑大學(xué) 信息與科學(xué)工程學(xué)院， 福建 廈門(mén) 361021)

?

采用小波變換和高斯過(guò)程的肌電信號(hào)模型預(yù)測(cè)

邵輝， 蘇芳茵， 程海波

(華僑大學(xué) 信息與科學(xué)工程學(xué)院， 福建 廈門(mén) 361021)

根據(jù)表面肌電信號(hào)的生物電信號(hào)特點(diǎn)，采用小波變換和高斯過(guò)程建模的方法對(duì)表面肌電信號(hào)進(jìn)行建模和預(yù)測(cè).對(duì)非線性的表面肌電信號(hào)利用擬合能力強(qiáng)大的高斯過(guò)程進(jìn)行建模，預(yù)測(cè)效果較好，但所需運(yùn)算時(shí)間長(zhǎng).針對(duì)其運(yùn)算時(shí)間長(zhǎng)的缺點(diǎn)進(jìn)行改進(jìn)，將預(yù)處理后的表面肌電信號(hào)小波分解，對(duì)分解后的系數(shù)高斯建模，然后重構(gòu).實(shí)驗(yàn)結(jié)果表明：該改進(jìn)方法在響應(yīng)時(shí)間和預(yù)測(cè)誤差方面效果明顯.

表面肌電信號(hào)； 高斯過(guò)程； 小波變換； 模型預(yù)測(cè)

表面肌電信號(hào)(surface electromyographic signal,sEMG)是人體在運(yùn)動(dòng)時(shí)骨骼肌產(chǎn)生的電信號(hào)在皮膚表面處的疊加，是一種信號(hào)微弱、噪聲復(fù)雜的非線性生物電信號(hào)[1].表面肌電信號(hào)作為智能手臂和假肢的控制源，其采集和處理要求很高.為了更好地控制機(jī)械臂，需要對(duì)肌電信號(hào)進(jìn)行建模研究.目前,對(duì)其進(jìn)行建模預(yù)測(cè)的研究相對(duì)較少，而常見(jiàn)的非線性建模方法在實(shí)際應(yīng)用中都或多或少存在問(wèn)題.高斯過(guò)程(Gaussian process,GP)于20世紀(jì)90年代進(jìn)入機(jī)器學(xué)習(xí)領(lǐng)域[2-4]，已成為非線性機(jī)器學(xué)習(xí)領(lǐng)域的研究熱點(diǎn)[5-7].它可以很好地避免這些問(wèn)題，并且具有學(xué)習(xí)參數(shù)少，對(duì)處理小樣本、高維數(shù)、非線性等復(fù)雜問(wèn)題有良好的適應(yīng)性，且泛化能力強(qiáng)的優(yōu)點(diǎn).高斯過(guò)程模型應(yīng)用非常廣泛.Williams等[8]率先將高斯過(guò)程應(yīng)用于機(jī)器學(xué)習(xí)領(lǐng)域；Murray-Smith等[9]分析比較了基于蒙特卡羅方法的高斯過(guò)程先驗(yàn)?zāi)Ｐ团c多模型方法；Kocijan[10]率先提出了一種將高斯過(guò)程模型中的方差引入控制信號(hào)的預(yù)測(cè)控制方法.高斯過(guò)程模型善于處理小樣本、高維數(shù)、非線性等復(fù)雜問(wèn)題.因此,本文首次將高斯過(guò)程模型應(yīng)用到肌電信號(hào)的建模,并引入了小波變換(wavelet transform，WT).在高斯過(guò)程回歸前,對(duì)樣本集進(jìn)行小波分解處理，從而減少模型學(xué)習(xí)的時(shí)間以提高預(yù)測(cè)速度.

1 表面肌電信號(hào)高斯過(guò)程模型

1.1 表面肌電信號(hào)描述

對(duì)表面肌電信號(hào)的處理方法是肌電信號(hào)分析和應(yīng)用的重要部分，其處理方法的選擇很大程度上要考慮肌電信號(hào)的產(chǎn)生機(jī)理.表面肌電信號(hào)的數(shù)學(xué)模型主要有4種：線型系統(tǒng)模型、雙極型模型、集中參數(shù)模型和非平穩(wěn)模型.前3種模型均屬于穩(wěn)態(tài)分析，即當(dāng)肌肉用力恒定的平穩(wěn)情況；當(dāng)肌肉用力不均時(shí)，肌電信號(hào)變化不平穩(wěn)，可以用非平穩(wěn)模型表征肌電信號(hào)，即

(1)

圖1 非平穩(wěn)模型Fig.1 Non-stationary model

式(1)中：c(t)為肌肉收縮程度的調(diào)制信號(hào)；m(t)為高斯噪聲(零均值單位方差)載波信號(hào).表面肌電信號(hào)非平穩(wěn)模型,如圖1所示.

1.2 高斯過(guò)程建模

高斯過(guò)程是在貝葉斯學(xué)習(xí)理論基礎(chǔ)上,發(fā)展起來(lái)的機(jī)器學(xué)習(xí)方法.它描述了任意有限變量集合的聯(lián)合密度函數(shù),屬于高斯分布的隨機(jī)過(guò)程.對(duì)于任意有限隨機(jī)變量x1,…,xn，其中,n≥1,且為任意整數(shù)，與其對(duì)應(yīng)的過(guò)程狀態(tài)f(x1),f(x2),…,f(xn)的聯(lián)合概率分布服從n維高斯分布.其具體過(guò)程可由均值m(x)與協(xié)方差函數(shù)k(x,x′)確定,即

(2)

給定訓(xùn)練D={(xi,yi)|i=1,2,…,N}，其中，xi∈Rd為輸入量，yi∈R是輸出向量.則輸出向量y由零均值函數(shù)m(x)和一個(gè)正定的協(xié)方差函數(shù)k(x,x′)的高斯先驗(yàn)分布組成，即y～N(m,p).

假設(shè)輸入新的向量x*∈Rd，KN為訓(xùn)練樣本的協(xié)方差矩陣，則其相應(yīng)的待預(yù)測(cè)數(shù)據(jù)的高斯分布函數(shù)為y～N(m*,p*)，其中

(3)

協(xié)方差函數(shù)的選擇沒(méi)有統(tǒng)一的指導(dǎo)性原則[11]，只需滿(mǎn)足對(duì)稱(chēng)性和半正定性.目前,常用的協(xié)方差函數(shù)為平方指數(shù)協(xié)方差函數(shù)，也稱(chēng)為徑向基函數(shù)(RBF)或高斯函數(shù)，其形式為

(4)

式(4)中：x和x′為訓(xùn)練集或測(cè)試集中的某個(gè)變量組合；σf,l,σn為超參數(shù)；δ為符號(hào)函數(shù).一般常用極大似然法求解最優(yōu)超參數(shù)[12].

2 小波變換結(jié)合高斯模型

高斯過(guò)程模型較適用于處理小樣本、高維數(shù)的問(wèn)題，而表面肌電信號(hào)的樣本數(shù)量多，其運(yùn)算量與訓(xùn)練集樣本數(shù)成正比.因此,在訓(xùn)練集樣本數(shù)較大時(shí)，計(jì)算逆矩陣所需要的時(shí)間和空間將成倍增大.考慮到這種情況，提出了小波變換結(jié)合高斯模型的預(yù)測(cè)方法.

小波變換的實(shí)質(zhì)是將信號(hào)分解為一系列小波函數(shù)的疊加，而這些小波函數(shù)都是由一個(gè)母小波函數(shù)經(jīng)過(guò)平移和尺度伸縮得來(lái)的.對(duì)信號(hào)逐步進(jìn)行時(shí)間頻率多尺度細(xì)化，最終達(dá)到高頻處時(shí)間細(xì)分、低頻處頻率細(xì)分，能自動(dòng)適應(yīng)時(shí)頻信號(hào)分析的要求，從而可聚焦到信號(hào)的任意細(xì)節(jié).

通過(guò)多分辨率離散小波變換將肌電信號(hào)y(t)分解,即

(5)

式(5)中：φ(t)為尺度函數(shù);ψ(t)為小波函數(shù).則尺度系數(shù)和小波系數(shù)的表達(dá)式為

(6)

(7)

將小波分解產(chǎn)生的尺度系數(shù)和小波系數(shù)分別進(jìn)行高斯過(guò)程訓(xùn)練，則其相應(yīng)的待預(yù)測(cè)數(shù)據(jù)的高斯分布函數(shù)為cj～N(mc,pc),dj～N(md,pd).

(a) GP預(yù)測(cè)模型 (b) WT-GP預(yù)測(cè)組合模型 圖2 2種預(yù)測(cè)模型的對(duì)比Fig.2 Comparison of two forecast models

高斯過(guò)程模型(GP)和小波變換結(jié)合高斯過(guò)程(WT-GP)模型在實(shí)現(xiàn)過(guò)程上不同，具體的模型結(jié)構(gòu),如圖2所示.WT-GP預(yù)測(cè)模型選用db5小波，對(duì)表面肌電信號(hào)進(jìn)行三級(jí)分解.db5小波具有較好的正則性，引入的光滑誤差不易被察覺(jué)，使信號(hào)重構(gòu)過(guò)程較光滑.采用小波變換結(jié)合高斯過(guò)程模型對(duì)表面肌電信號(hào)進(jìn)行建模預(yù)測(cè).首先，將樣本數(shù)據(jù)進(jìn)行小波三級(jí)分解，第一級(jí)可以得到一組低頻系數(shù)和一組高頻系數(shù)，繼續(xù)分解低頻部分，以此類(lèi)推，最終可以得到一組低頻系數(shù)和三組高頻系數(shù)，其樣本數(shù)遠(yuǎn)小于原來(lái)的樣本數(shù)據(jù);然后,對(duì)四組系數(shù)分別進(jìn)行高斯過(guò)程回歸建模；最后,將每個(gè)模型預(yù)測(cè)的系數(shù)重組得到最終的預(yù)測(cè)結(jié)果.

小波變換有多尺度的分析能力，將信號(hào)按照不同的尺度分解，降低信號(hào)非線性程度.由于非線性程度降低，預(yù)測(cè)信號(hào)的誤差減小，因此，提高了預(yù)測(cè)的精確度.而在小波分解后，低頻系數(shù)和高頻系數(shù)的維度明顯低于信號(hào)本身，所以對(duì)其分別求逆的運(yùn)算量也低于對(duì)信號(hào)本身求逆.即使低頻和高頻系數(shù)維度相加等于信號(hào)自身，但是求逆過(guò)程的計(jì)算難度確隨著矩陣維數(shù)成倍增長(zhǎng)，故運(yùn)算時(shí)間明顯減少.由此可以得出，引入小波后的高斯過(guò)程模型預(yù)測(cè)回歸在誤差和運(yùn)算時(shí)間上均有優(yōu)化.

3 實(shí)驗(yàn)仿真

3.1 肌電信號(hào)采集方案

實(shí)驗(yàn)采集設(shè)備：表面電極貼、PC機(jī)、PCL818L多功能數(shù)據(jù)采集卡、肌電傳感器.實(shí)驗(yàn)采用的傳感器是日本司機(jī)工工程技術(shù)公司的TYE-1000M雙極型表面電極，其中,包含傳感器和放大器，配合使用一次性銀/氯化銀按鈕式電極貼片.

圖3 現(xiàn)場(chǎng)信號(hào)采集實(shí)驗(yàn)圖Fig.3 Field signal acquisition experimental picture

采集系統(tǒng)的構(gòu)建：實(shí)驗(yàn)平臺(tái)采用MATLAB/Simulink/xPC Target的實(shí)時(shí)環(huán)境建立整個(gè)肌電信號(hào)采集系統(tǒng)，采樣頻率1 000 Hz，信號(hào)主要能量集中在10～500 Hz.現(xiàn)場(chǎng)信號(hào)采集實(shí)驗(yàn),如圖3所示.實(shí)驗(yàn)方案：實(shí)驗(yàn)前,上臂經(jīng)過(guò)清洗處理，減少皮膚阻抗，且實(shí)驗(yàn)之前未進(jìn)行激烈運(yùn)動(dòng);實(shí)驗(yàn)時(shí)，將電極貼片黏在肱二頭肌和肱三頭肌上，受測(cè)者坐在椅子上，手成握拳狀，手臂自然下垂，采集兩種動(dòng)作前臂屈伸和前臂旋轉(zhuǎn)各10次.

采集后的肌電信號(hào)含有噪聲，需要進(jìn)行濾波.文中選擇常用的切比雪夫I型濾波器，其通帶下限頻率為20 Hz，上限頻率為510 Hz；阻帶下限為10 Hz，阻帶上限為520 Hz；通帶衰減小于1 dB，阻帶衰減大于60 dB.由此，可確定切比雪夫?yàn)V波器階數(shù)為4.

3.2 高斯過(guò)程模型預(yù)測(cè)結(jié)果

由理論分析可知：訓(xùn)練集和測(cè)試集是兩個(gè)不同的集合，但這兩個(gè)集合是由同一個(gè)動(dòng)作產(chǎn)生，因此，這兩個(gè)集合在本質(zhì)上是屬于同一個(gè)樣本點(diǎn)集，其預(yù)測(cè)效果會(huì)比較好.如果測(cè)試集由其他動(dòng)作產(chǎn)生，其效果會(huì)有所不同.

采集2種動(dòng)作的肌電信號(hào)，從中分別抽取1 500組和1 000組樣本數(shù)據(jù)，選擇前一種動(dòng)作的500組數(shù)據(jù)作為訓(xùn)練樣本，剩下的1 000組數(shù)據(jù)作為測(cè)試樣本；而后一種動(dòng)作的1 000組僅作為測(cè)試樣本.如此，兩個(gè)測(cè)試集中，一個(gè)和樣本集屬于同一集合，一個(gè)屬于不同集合.整個(gè)學(xué)習(xí)訓(xùn)練的過(guò)程如上敘述，分別對(duì)1 000組測(cè)試樣本進(jìn)行預(yù)測(cè)，結(jié)果如圖4，5所示.

由圖4可知：不管測(cè)試樣本與訓(xùn)練樣本是否來(lái)自同一動(dòng)作數(shù)集，其預(yù)測(cè)值與實(shí)際值的曲線都幾乎相似.圖5中:誤差曲線和方差曲線的精確度數(shù)量級(jí)很小.由此可以證明，高斯過(guò)程對(duì)非線性回歸有很強(qiáng)大的學(xué)習(xí)能力和泛化能力.

由圖5可知：屬于同一集合測(cè)試樣本的預(yù)測(cè)誤差曲線優(yōu)于不同集合測(cè)試樣本的預(yù)測(cè)效果；同一集合的預(yù)測(cè)方差明顯低于不同集合的方差.這說(shuō)明訓(xùn)練后的高斯過(guò)程模型對(duì)同一集合的預(yù)測(cè)具有更大的可信度.在高斯過(guò)程建模時(shí)，根據(jù)訓(xùn)練集500組數(shù)據(jù)求取最優(yōu)化參數(shù).因此，它對(duì)訓(xùn)練集內(nèi)1 000組數(shù)據(jù)預(yù)測(cè)時(shí)，準(zhǔn)確度高于訓(xùn)練集外部的1 000組數(shù)據(jù).由此可知：高斯過(guò)程模型對(duì)訓(xùn)練集外點(diǎn)預(yù)測(cè)的可信度已經(jīng)降低，訓(xùn)練集內(nèi)的預(yù)測(cè)效果高于訓(xùn)練集外的預(yù)測(cè)效果.

(a) 同一數(shù)集預(yù)測(cè)效果 (b) 不同數(shù)集預(yù)測(cè)效果圖4 實(shí)際信號(hào)與預(yù)測(cè)信號(hào)對(duì)比圖Fig.4 Actual signals and forecast comparison chart

(a) 同一集合誤差曲線 (b) 不同集合誤差曲線

(c) 同一集合方差曲線 (d) 不同集合方差曲線圖5 預(yù)測(cè)的誤差和方差曲線圖Fig.5 Prediction error and variance of graph

3.3 基于小波的高斯過(guò)程模型預(yù)測(cè)結(jié)果

采用上面的500組作為訓(xùn)練樣本，不同集合的1 000組作為測(cè)試樣本.經(jīng)過(guò)仿真結(jié)果比較可知：在誤差分析和響應(yīng)速度上，WT-GP模型預(yù)測(cè)的效果都有良好地改善.引入小波預(yù)測(cè)后，在誤差分析上有變化.在訓(xùn)練和預(yù)測(cè)前,將信號(hào)進(jìn)行小波三層分解，將小波分解得到的低頻系數(shù)ca1和高頻系數(shù)cd1，cd2，cd3分別訓(xùn)練和預(yù)測(cè)，分別重構(gòu)得到低頻分量和高頻分量.各分量的實(shí)際值和預(yù)測(cè)值,如圖6所示.

將得到的低頻系數(shù)和高頻系數(shù)組合重構(gòu)，得到預(yù)測(cè)后的信號(hào)，如圖7所示.由圖7可知：預(yù)測(cè)信號(hào)和實(shí)際信號(hào)沒(méi)有太大差別，而誤差和方差的數(shù)量級(jí)依然很小.實(shí)驗(yàn)仿真用的數(shù)據(jù)不是同一集合，由此可見(jiàn)，其擬合能力和泛化能力依舊沒(méi)變.比較圖5與圖7中非同一集合預(yù)測(cè)的誤差，圖7的誤差更小一些，和圖5中同一集合預(yù)測(cè)誤差范圍比較接近.

(a) 小波分解低頻分量a3 (b) 預(yù)測(cè)低頻分量a3

(c) 小波分解高頻分量d1 (d) 預(yù)測(cè)高頻分量d1

(e) 小波分解高頻分量d2 (f) 預(yù)測(cè)高頻分量d2

(g) 小波分解高頻分量d3 (h) 預(yù)測(cè)高頻分量d3圖6 各分量的實(shí)際值和預(yù)測(cè)值Fig.6 Each component of actual value and predictive value

(a) 實(shí)際信號(hào) (b) 誤差曲線

(c) 預(yù)測(cè)信號(hào) (d) 方差曲線圖7 引入小波后的預(yù)測(cè)效果Fig.7 Prediction effect after introducing wavelet

兩種預(yù)測(cè)模型運(yùn)算時(shí)間的比較，如表1所示.表1中，運(yùn)行時(shí)間指訓(xùn)練和測(cè)試所共用的時(shí)間.由表1可知：引入小波的高斯模型訓(xùn)練時(shí)間明顯低于沒(méi)有使用小波的情況.不論測(cè)試集合和訓(xùn)練集合是否來(lái)自同一集合，運(yùn)行時(shí)間都會(huì)隨著訓(xùn)練樣本數(shù)的下降而明顯變??；而當(dāng)訓(xùn)練組數(shù)確定時(shí)，同一集合的時(shí)間略小于不同集合的運(yùn)行時(shí)間.經(jīng)過(guò)小波分解后的訓(xùn)練樣本數(shù)減少，運(yùn)行的時(shí)間有很大的變化.

表1 兩種預(yù)測(cè)模型運(yùn)算時(shí)間的比較

4 結(jié)束語(yǔ)

表面肌電信號(hào)是非線性的生物電信號(hào)，對(duì)其進(jìn)行高斯過(guò)程回歸建模預(yù)測(cè)，能很好地?cái)M合其非線性.高斯過(guò)程建模時(shí)，首先,確定核函數(shù)；然后,訓(xùn)練高斯過(guò)程模型，從而得到最優(yōu)超參數(shù)；最后,用訓(xùn)練好的模型進(jìn)行預(yù)測(cè).仿真結(jié)果證明，高斯過(guò)程對(duì)表面肌電信號(hào)有良好的擬合能力和強(qiáng)大的泛化性能.在此基礎(chǔ)上，提出了基于小波的高斯建模，解決了單純使用高斯過(guò)程建模帶來(lái)的運(yùn)行時(shí)間長(zhǎng)和非同一集合預(yù)測(cè)誤差大的缺點(diǎn).從仿真結(jié)果來(lái)看，引入小波后的預(yù)測(cè)效果在減小誤差和運(yùn)算時(shí)間上都有明顯的優(yōu)勢(shì).

[1] 張啟忠,席旭剛,羅志增.多重分形分析在肌電信號(hào)模式識(shí)別中的應(yīng)用[J].傳感技術(shù)學(xué)報(bào),2013,26(2):282-288.

[2] RASMUSSEN C E,WILLIAMS C K I.Gaussian processes for machine learning (adaptive computation and machine learning)[M].American:MIT Press,2005:49-51.

[3] KOCIJAN J,GRANCHAROVA A.Application of gaussian processes to the modelling and control in process engineering[M]∥Studies in Computational Intelligence.Germany:Springer Berlin Heidelberg，2014:155-190.

[4] PARK C,HUANG Jianhua,DING Yu.Domain decomposition approach for fast gaussian process regression of large spatial data sets[J].Journal of Machine Learning Research,2011,12(4):1697-1728.

[5] HE Zhikun,LIU Guangbin,ZHAO Xijing,et al.Temperature model for FOG zero-bias using Gaussian process regression[J].Advances in Intelligent Systems and Computing,2012,180(1):37-45.

[6] 李鵬,宋申民,段廣仁.改進(jìn)的平方根UKF及其在交會(huì)對(duì)接中的應(yīng)用[J].電機(jī)與控制學(xué)報(bào),2010,14(11):100-104.

[7] 孫斌,姚海濤,劉婷.基于高斯過(guò)程回歸的短期風(fēng)速預(yù)測(cè)[J].中國(guó)電機(jī)工程學(xué)報(bào),2012,32(29):104-109.

[8] WILLIAMS C K I,RASMUSSEN C E.Gaussian processes for regression[J].Advances in Neural Information Processing Systems Pages,1996,27(6):514-520.

[9] MURRAY-SMITH R,JOHANSEN T A,SHORTEN R.On transient dynamics, off-equilibrium behaviour and identification in blended multiple model structures[C]∥Proceedings of the European Control Conference.Karslruhe：IEEE Press，1999:BA-14.

[11] 張惠澤.基于高斯過(guò)程的pH中和過(guò)程控制研究[D].哈爾濱：哈爾濱工業(yè)大學(xué),2010:26-29.

[12] 陳寶林.最優(yōu)化理論與算法[M].北京：清華大學(xué)出版社,1989:280-288.

(責(zé)任編輯： 黃曉楠 英文審校： 吳逢鐵)

Model Forecasting of EMG Using Wavelet Transformation and Gaussian Process

SHAO Hui, SU Fangyin, CHENG Haibo

(College of Information Science and Engineering, Huaqiao University, Xiamen 361021, China)

According to the characteristics of the surface EMG signal, this paper uses wavelet transform and Gauss process modeling method to model and predict the surface EMG signal. The nonlinear surface EMG signal is used to model the fitting ability of the Gauss process, and the prediction effect is better, but the operation time is longer. To overcome the shortcomings of the long computation time, the wavelet decomposition of the surface EMG signal is processed, and the coefficients of the decomposition are modeled in Guassian. Experimental results show that the improved method has obvious effect on response time and prediction error.

surface electromyogram; Gaussian process; wavelet transform; model prediction

10.11830/ISSN.1000-5013.201606016

2016-01-14

邵輝(1973-)，女，副教授，博士，主要從事機(jī)器人控制、運(yùn)動(dòng)規(guī)劃、智能控制、非線性系統(tǒng)LPV建模的研究.E-mail:shaohuihull@163.com.

福建省科技計(jì)劃項(xiàng)目(2015H0026)； 教育部留學(xué)回國(guó)人員科研啟動(dòng)基金資助項(xiàng)目(Z1534004)； 福建省泉州市科技計(jì)劃項(xiàng)目(2013Z34)

TP 391

A

1000-5013(2016)06-0743-06