基于VFFRLS濾波器的腦電時(shí)序隨機(jī)噪聲去除方法

連召洋 段立娟* 陳軍成 喬元華

1(北京工業(yè)大學(xué)信息學(xué)部 北京 100124) 2(可信計(jì)算北京市重點(diǎn)實(shí)驗(yàn)室 北京 100124) 3(信息安全等級(jí)保護(hù)關(guān)鍵技術(shù)國(guó)家工程實(shí)驗(yàn)室 北京 100124) 4(北京工業(yè)大學(xué)應(yīng)用數(shù)理學(xué)院 北京 100124)

0 引 言

顱外的腦電信號(hào)(Electro Encephalography,EEG)中包含了大量的生理和病理信息,但也可能包含有大量的偽跡[1],這些偽跡噪聲嚴(yán)重影響了采集到的腦電信號(hào)[2]。近年來(lái)出現(xiàn)了許多偽跡去除方法,例如:Phadikar等[3]提出通過(guò)啟發(fā)式優(yōu)化閾值的小波變換的方法自動(dòng)去除眼電偽跡;Pushpendra等[4]采用Fourier分解的方法去除心電中的基線漂移和工頻干擾偽跡;Yasoda等[5]采用WICA(Wavelet Independent Component Analysis)算法識(shí)別眼電和電機(jī)運(yùn)動(dòng)偽跡;Ahmed等[6]采用RLS(Recursive Least Square)自適應(yīng)濾波算法去除工頻干擾偽跡;文獻(xiàn)[7]采用ICA和小波算法去除腦點(diǎn)信號(hào)中的肌電偽跡;Malhotra等[8]采用小波濾波算法去除心電中的基線漂移偽跡;文獻(xiàn)[9]分別采用Wiener和Kalman算法去除心電信號(hào)中的噪聲;Altay等[10]采用寬帶陷波器去除心電信號(hào)中的偽跡;文獻(xiàn)[11]采用子代自適應(yīng)濾波器去除心電信號(hào)中的偽跡。但現(xiàn)在大部分腦電偽跡去除方法都是集中在眼電偽跡、心電偽跡、電機(jī)運(yùn)動(dòng)偽跡、基線漂移偽跡和工頻干擾等偽跡的去除,但對(duì)腦電信號(hào)中較強(qiáng)時(shí)序隨機(jī)噪聲的偽跡去除研究較少。

自適應(yīng)濾波器可以根據(jù)輸入信號(hào)自動(dòng)不斷調(diào)整參數(shù)的數(shù)字濾波器,能夠解決非自適應(yīng)濾波器參數(shù)固定的缺陷。腦電是一種隨機(jī)信號(hào),而加入腦電信號(hào)中的時(shí)序隨機(jī)噪聲也是復(fù)雜多變的,事先不知道模型中所需要設(shè)定的參數(shù),所以采用自適應(yīng)濾器來(lái)自適應(yīng)調(diào)節(jié)參數(shù)可以達(dá)到偽跡去除的目的。

常用的自適應(yīng)濾波器分為基礎(chǔ)的LMS(Least Mean Square)和RLS(Recursive Least Squares)及其改進(jìn)算法。在文獻(xiàn)[12]中對(duì)比了基礎(chǔ)的LMS和RLS自適應(yīng)算法在有源濾波器的自適應(yīng)控制中的效果。

LMS的步長(zhǎng)是固定的,濾波器通過(guò)計(jì)算每一次輸入值與目標(biāo)值的梯度來(lái)到得自相關(guān)矩陣和互相關(guān)向量,從而來(lái)估計(jì)噪聲[13]。歸一化的LMS算法NLMS(Normalized LMS)[14]是LMS的改進(jìn),其步長(zhǎng)是一個(gè)隨時(shí)間變化的量,與上個(gè)時(shí)刻輸入信號(hào)的自身能量值有關(guān),輸入能量大的,調(diào)整一個(gè)小的步長(zhǎng),輸入能量小的,調(diào)整一個(gè)大的步長(zhǎng),從而保證輸入信號(hào)的大小對(duì)算法的收斂速度影響不太大。系數(shù)比例自適應(yīng)算法PNLMS(Proportionate NLMS)[15]是一種在處理稀疏信號(hào)時(shí)可以快速收斂的算法,其收斂速度比NLMS要快,但是處理非稀疏的信號(hào)時(shí),其收斂速度比NLMS要慢。IPNLMS(Improved PNLMS)[16]是一種改進(jìn)的PNLMS算法,改善了PNLMS的非稀疏信號(hào)情況下的收斂速度,但它需要預(yù)先確定一些行為參數(shù),給使用帶來(lái)不便。UPNLMS(Under weigh of IPNLMS and NLMS)[17]是一種將IPNLMS與NLMS相融合的算法,采用兩次估計(jì)和兩次更新的方法,第一次估計(jì)后采用IPNLMS算法進(jìn)行更新,第二次估計(jì)后,采用NLMS算法進(jìn)行更新,該算法有較快收斂速度。

遞歸最小二乘算法RLS[18-19]充分利用過(guò)去的信息,使每個(gè)時(shí)刻的期望信號(hào)與濾波器輸出的估計(jì)信號(hào)的平方和最小,提高了收斂速度,克服了LMS收斂速度過(guò)慢的問(wèn)題?？焖贆M向遞歸最小二乘算法FTRLS(Fast Transversal RLS)[20]是一種RLS的快速實(shí)現(xiàn)方法,在時(shí)序上同時(shí)求解前向和后向更新,外加橫向的聯(lián)合估計(jì)過(guò)程和輔助濾波過(guò)程,這四個(gè)過(guò)程并行并且交換變量,從而實(shí)現(xiàn)快速RLS的目的。

對(duì)于自適應(yīng)濾器,遺忘因子的大小與算法的收斂速度、穩(wěn)定性有密切的關(guān)系,為此本文提出采用可變遺忘因子的自適應(yīng)濾波器VFFRLS(Variable Forgetting Factor RLS)來(lái)去除腦電信號(hào)中的時(shí)序隨機(jī)噪聲,并對(duì)比了不同算法在運(yùn)動(dòng)想象腦電信號(hào)中混入較強(qiáng)的時(shí)序隨機(jī)噪聲后的去除效果。與其他7種自適應(yīng)濾波方法相比,VFFRLS可變遺忘因子的自適應(yīng)濾波算法在運(yùn)動(dòng)想象的腦電信號(hào)中的時(shí)序隨機(jī)噪聲的去除效果最好,在6個(gè)電極中每個(gè)樣本重復(fù)試驗(yàn)20次,VFFRLS算法的平均根均方差RMSE都最小,SNR都最大。

1 自適應(yīng)濾波器的腦電時(shí)序隨機(jī)噪聲去除框架

圖1中的SEEG,noise表示含有時(shí)序隨機(jī)強(qiáng)噪聲的混合腦電信號(hào),Noise表示參考的隨機(jī)噪聲信號(hào),圖1中的自適應(yīng)濾波器為各種經(jīng)典的自適應(yīng)濾波器及其改進(jìn)算法,包含LMS、NLMS、PNLMS、IPNLMS、UPNLMS、RLS、FTRLS、VFFRLS自適應(yīng)濾波器,Noise,est是經(jīng)過(guò)自適應(yīng)濾波器計(jì)算出來(lái)的噪聲信號(hào),SEEG,est是計(jì)算出來(lái)不含隨機(jī)噪聲的腦電信號(hào)。

圖1 自適應(yīng)濾波算法的腦電時(shí)序隨機(jī)噪聲去除框架圖

2 基于VFFRLS的腦電時(shí)序隨機(jī)噪聲去除方法

VFFRLS[21]是在權(quán)重更新過(guò)程中引入可變的遺忘因子。遺忘因子的大小和算法的收斂速度、穩(wěn)定性有密切的關(guān)系。當(dāng)遺忘因子較小時(shí),距離當(dāng)前時(shí)間點(diǎn)近的數(shù)據(jù)起作用比較大,距離當(dāng)前時(shí)間點(diǎn)遠(yuǎn)的數(shù)據(jù)作用會(huì)被減弱,從而提高算法的收斂速度,但穩(wěn)定性相對(duì)較差。當(dāng)遺忘因子為1時(shí),就是普通的最小二乘法。因此為兼顧穩(wěn)定性和跟蹤速度,找到合適的遺忘因子非常重要。在算法開(kāi)始的時(shí)候,先驗(yàn)知識(shí)少,誤差大,可以設(shè)定一個(gè)較小的遺忘因子,提高收斂速度,隨著時(shí)間的推移,逐漸增大遺忘因子,提高算法的穩(wěn)定性。

基于該思想的腦電的時(shí)序隨機(jī)強(qiáng)噪聲去除的具體步驟如下:

(1) 將混合有時(shí)序隨機(jī)強(qiáng)噪聲的腦電信號(hào)和參考噪聲信號(hào)放進(jìn)VFFRLS自適應(yīng)濾波器。

(2) 設(shè)置VFFRLS的初始參數(shù)。初始化遺忘因子λ0、δ、w0、y0、M、P0等參數(shù),其中λ0=0.99,δ=0.001,M=3,l=M-1=2,P0=δ·M,y0、w0是1×3維的0向量。

(3) 進(jìn)行VFFRLS自適應(yīng)濾波。

首先,生成輸入的混合信號(hào)序列片段SEEG,noiset和參考噪聲片段xt。

SEEG,noiset=[SEEG,noise(t),…,SEEG,noise(t+l)]

(1)

xt=[Noise(t),…,Noise(t+l)]

(2)

式中:SEEG,noiset和x都是一個(gè)1×(l+1)維的向量,SEEG,noise(t)表示向量SEEG,noise在索引t處的數(shù)據(jù),SEEG,noiset為在第t次迭代時(shí)混合信號(hào)序列片段的向量。

其次,更新遺忘因子λt。

如果迭代次數(shù)1≤t≤λnum-1則更新式(3),如果λnum-1

式中:λnum取90,N為896。遺忘因子對(duì)算法有重要影響,太大穩(wěn)定性好,但跟蹤能力差,太小跟蹤能力強(qiáng),但是穩(wěn)定性差,一般來(lái)說(shuō)遺忘因子的大小在0.9到1.0之間,算法效果會(huì)比較好,在本文中初始值為0.99,開(kāi)始取相對(duì)小的值,提高其跟蹤速度,隨著時(shí)間的增加,逐漸提高遺忘因子,提高其穩(wěn)定性。

再次,估算不含噪聲的腦電信號(hào)片段SEEG,est。

SEEG,est(t+l)=SEEG,noise(t+l)-Noise,est(t+l)

(5)

估算噪聲片段Noise,est。

最后,更新wt、Kt、Pt等其他中間參數(shù)。

式中:SEEG,est是估計(jì)出來(lái)的不含隨機(jī)噪聲的腦電信號(hào),SEEG,est(t+l)中的小括號(hào)表示SEEG,est向量中第t+l索引處的數(shù)據(jù),Pt的下標(biāo)為第t次迭代時(shí)更新的P值,上述公式中其他變量含t的小括號(hào)和下標(biāo)含義也一樣。

如果1

算法1VFFRLS算法更新流程

輸入 混合腦電SEEG,noise和參考噪聲SEEG,org

初始化λ0、δ、w0、y0、M、P0等參數(shù)

for循環(huán)t=1,2,…,N-M+1

生成輸入的混合信號(hào)片段SEEG,noiset和參考噪聲片段xt

ift≤λnum-1,更新λtend

估算不含噪聲的腦電信號(hào)片段SEEG,noise

估算噪聲片段Noise,est

更新wt、Kt、Pt等其他中間參數(shù)

end

輸出 不含噪聲的腦電信號(hào)SEEG,noise

3 基于VFFRLS的腦電時(shí)序隨機(jī)噪聲去除方法

采用BCI(Brain Computer Interface)腦機(jī)接口競(jìng)賽第二屆Ia數(shù)據(jù)集[22],數(shù)據(jù)集的受試者都是健康的。該數(shù)據(jù)集是由圖賓根大學(xué)的Niels Birbaumer提供的。數(shù)據(jù)采集于被試者的大腦皮層,實(shí)驗(yàn)任務(wù)是被試者根據(jù)電腦屏幕光標(biāo)上下移動(dòng)想象運(yùn)動(dòng)的光標(biāo),儀器記錄大腦皮層對(duì)應(yīng)電極的電位值。每個(gè)樣本持續(xù)時(shí)間為6 s,這6 s中,前1 s時(shí)間是休息,中間1.5 s時(shí)間是運(yùn)動(dòng)想象的提示,后3.5 s時(shí)間是信息反饋。其中后3.5 s被6個(gè)電極以256 Hz記錄下來(lái)作為樣本,一共有561個(gè)樣本,每個(gè)樣本中有6個(gè)點(diǎn)電極樣本段,每段有896個(gè)樣本點(diǎn)。

隨機(jī)噪聲的生成:

Noise=a0·frnd(m,n)

(10)

式中:a0為噪聲的幅值系數(shù),m為1,n為896。frnd(m,n)是生成m×n維的0到1之間隨機(jī)數(shù)的函數(shù),Noise是一個(gè)隨機(jī)生成的1×896維的行向量。含有時(shí)序隨機(jī)信號(hào)的混合腦電信號(hào)生成方式如下:

SEEG,noise(i)=SEEG,org(i)+1·a1·Noise(i)+

1i>1·a2·Noise(i-1)+

1i>2·Noise(i-2)

(11)

式中:SEEG,org是每個(gè)樣本中每個(gè)電極的樣本段,SEEG,noise是加了時(shí)序隨機(jī)強(qiáng)噪聲的混合信號(hào),延時(shí)添加不同權(quán)重噪聲體現(xiàn)時(shí)序性,a1和a2為權(quán)重系數(shù),其中01表示當(dāng)i>1時(shí)取1,反之取0。

圖2中每個(gè)子圖的最上面的信號(hào)是原始腦電信號(hào)SEEG,org,是從561個(gè)樣本的對(duì)應(yīng)電極中隨機(jī)抽取的;每個(gè)子圖的中間是生成的混合信號(hào)SEEG,noise,可以看出混合信號(hào)已經(jīng)嚴(yán)重受時(shí)序隨機(jī)強(qiáng)噪聲影響;每個(gè)子圖的下面是經(jīng)過(guò)VFFRLS自適應(yīng)濾波算法估計(jì)出來(lái)的不含時(shí)序隨機(jī)噪聲的腦電信號(hào)SEEG,est,可以看到效果很好,和原始的腦電信號(hào)很接近。

(a) A1

(b) F3圖2 不同電極中原始、混合后、估計(jì)出來(lái)的腦電信號(hào)可視化

4 實(shí)例驗(yàn)證

4.1 評(píng)價(jià)指標(biāo)

(1) 均方根誤差。RMSE(Root Mean Square Error),可以衡量運(yùn)動(dòng)想象腦電中時(shí)序隨機(jī)強(qiáng)噪聲的在時(shí)序上的去除效果,評(píng)價(jià)計(jì)算的純凈信號(hào)和原始的純凈信號(hào)的差別程度,此值越小,說(shuō)明估計(jì)出來(lái)的純凈信號(hào)和原始的純凈信號(hào)越相近,其公式如下:

式中:N為樣本點(diǎn)數(shù),SEEG,org某樣本中某電極導(dǎo)聯(lián)純凈的腦電信號(hào),SEEG,est通過(guò)各種自適應(yīng)濾波器估計(jì)出來(lái)的該電極導(dǎo)聯(lián)的純凈EEG信號(hào)。

(2) 信噪比。SNR是信噪比(Signal Noise Ratio),本來(lái)是指通信系統(tǒng)中一臺(tái)電子設(shè)備中信號(hào)和噪聲的比例。本文中用信噪比表示原始純凈信號(hào)與估計(jì)出來(lái)的純凈EEG信號(hào)的誤差之間的關(guān)系,此值越大越好。

式中:N、SEEG,org、SEEG,est參數(shù)含義和RMSE中的參數(shù)含義一致。

4.2 對(duì)比分析

本文對(duì)比了8種自適應(yīng)濾波器對(duì)混有時(shí)序隨機(jī)噪聲的腦電信號(hào)的噪聲去除效果。本數(shù)據(jù)中有561個(gè)樣本,每個(gè)樣本有6個(gè)電極,每個(gè)方法對(duì)每個(gè)樣本中的每個(gè)電極導(dǎo)聯(lián)信號(hào)做20次實(shí)驗(yàn),所以每個(gè)方法在每個(gè)電極上的RMSE和SNR值都是分別做了561×20=11 220次實(shí)驗(yàn),然后分別對(duì)這些實(shí)驗(yàn)求平均值得到平均的RMSE和SNR值。相當(dāng)于本文做了8×6×(561×20)=538 560次實(shí)驗(yàn)。

在表1和表2中,水平方向代表不同的自適應(yīng)濾波算法,垂直方向代表不同的電極。從表1中可知,在8種自適應(yīng)濾波器中,VFFRLS的在6個(gè)電極中估計(jì)出來(lái)的純凈腦電信號(hào)與原始無(wú)隨機(jī)噪聲的腦電信號(hào)的平均RMSE最小。從表2可知,VFFRLS自適應(yīng)濾波器在6個(gè)電極上估計(jì)出來(lái)的純凈的腦電信號(hào)與原始無(wú)時(shí)序隨機(jī)噪聲的腦電信號(hào)之間的平均SNR都最大。在評(píng)價(jià)指標(biāo)中RMSE越小越好,SNR越大越好。所以本文提出的基于可變遺忘因子VFFRLS的腦電時(shí)序隨機(jī)噪聲去除算法在本對(duì)比實(shí)驗(yàn)中效果最好。

表1 不同自適應(yīng)濾波算法的平均RMSE值比較 單位:μV

表2 不同自適應(yīng)濾波算法的平均SNR值比較 單位:dB

續(xù)表2

5 結(jié) 語(yǔ)

現(xiàn)在大部分腦電偽跡去除方法都是集中在眼電偽跡、心電偽跡、電機(jī)運(yùn)動(dòng)偽跡、基線漂移偽跡和工頻干擾等偽跡的去除,但對(duì)腦電信號(hào)中較強(qiáng)時(shí)序隨機(jī)噪聲的偽跡去除研究較少。本文提出了基于可變遺忘因子VFFRLS的腦電時(shí)序隨機(jī)噪聲的偽跡去除算法,并且對(duì)比了8種主流的自適應(yīng)濾波算法在運(yùn)動(dòng)想象數(shù)據(jù)集上時(shí)序隨機(jī)強(qiáng)噪聲去除的效果,實(shí)驗(yàn)結(jié)果表明,VFFRLS在對(duì)比實(shí)驗(yàn)中效果最好,其RMSE最小,SNR最大。