基于廣義子空間法的腦電誘發(fā)電位單導(dǎo)少次提取

王永軒 邱天爽 劉 蓉

（大連理工大學(xué)，電子信息與電氣工程學(xué)部，大連 116024）

引言

腦電誘發(fā)電位（evoked potential，EP）是當(dāng)對(duì)人的感官加以特定的刺激時(shí)，中樞神經(jīng)系統(tǒng)所產(chǎn)生的有特定規(guī)律的生物電信號(hào)，它反映了相應(yīng)的感覺通路及皮層區(qū)域的神經(jīng)電活動(dòng)情況，可以客觀地評(píng)價(jià)人的感覺通路是否正常，是臨床醫(yī)學(xué)診斷神經(jīng)系統(tǒng)損傷及病變的重要手段之一［1］。因?yàn)檎T發(fā)電位是和自發(fā)腦電同時(shí)觀測(cè)到的，所以必須先通過(guò)一定的方法提取出來(lái)。目前臨床醫(yī)學(xué)通常使用的方法是疊加平均法［1］，也就是對(duì)感官進(jìn)行多次刺激后，將每次記錄的結(jié)果進(jìn)行疊加平均。當(dāng)刺激次數(shù)很大時(shí)，每次記錄中的EP信號(hào)波形相差會(huì)變大，從而導(dǎo)致結(jié)果不可靠，另外當(dāng)超出受試者的耐受力時(shí)會(huì)產(chǎn)生暴發(fā)偽跡，所以希望盡量減少刺激的次數(shù)，從而實(shí)現(xiàn)少次提取。

二十多年來(lái)，已經(jīng)有很多研究者提出了許多不同的方法來(lái)提取EP信號(hào)，比如各種濾波法［2］、神經(jīng)網(wǎng)絡(luò)法［3］、小波分析法［4-5］、獨(dú)立分量分析法［6］等。這些方法都有各自的特點(diǎn)，它們?cè)诮鉀Q各自的模型問(wèn)題中都得到了很好的結(jié)果，但是這些方法所建立的模型應(yīng)用于誘發(fā)電位單導(dǎo)少次提取問(wèn)題中時(shí)有一定的局限性，所以無(wú)法真正應(yīng)用到臨床。比如獨(dú)立分量分析法只能處理多導(dǎo)腦電信號(hào)，同時(shí)還要做很多假設(shè)，包括信號(hào)源必須是獨(dú)立的，信號(hào)源數(shù)目要不多于觀測(cè)信號(hào)路數(shù)，觀測(cè)信號(hào)是信號(hào)源的線性混合，混和矩陣是列滿秩的等等［7］。這些假設(shè)與實(shí)際并不完全相符，所以首先需要建立一個(gè)與實(shí)際問(wèn)題盡可能相符的模型，然后再解決這個(gè)模型問(wèn)題。

在所有的信號(hào)提取或信號(hào)增強(qiáng)問(wèn)題中，人們都希望盡可能地去除或降低噪聲而同時(shí)不破壞信號(hào)的完整性。只要運(yùn)用得當(dāng)，子空間法通常都可以很好地實(shí)現(xiàn)降低噪聲與保留信號(hào)之間的折衷。因而子空間法不僅在腦電信號(hào)分析中，而且在功率譜估計(jì)［8］、系統(tǒng)辯識(shí)［9］和語(yǔ)音處理［10-11］等信號(hào)分析領(lǐng)域中都有著廣泛的應(yīng)用。例如Ephraim等人使用子空間法實(shí)現(xiàn)在白噪聲背景中提取語(yǔ)音信號(hào)［10］，而Rezayee等人進(jìn)一步應(yīng)用在有色噪聲條件下同時(shí)并不需要對(duì)觀測(cè)信號(hào)進(jìn)行預(yù)白化［11］?；贓P信號(hào)與自發(fā)腦電相互獨(dú)立這一性質(zhì)，本研究提出了一種廣義子空間法實(shí)現(xiàn)EP信號(hào)的單導(dǎo)少次提取。

1 廣義子空間法

在單次刺激的情況下，EP信號(hào)是一段有限長(zhǎng)的隨機(jī)序列，令其為s（n），（n=1，2，…，N），令觀測(cè)信號(hào)為x（n），（n=1，2，…，N），其中的自發(fā)腦電為v（n），（n=1，2，…，N）。將上述序列表示為向量，分別為s，x，v，則有x=s+v。

無(wú)論是經(jīng)典濾波器還是現(xiàn)代濾波器，最終需要確定的都是由一組系數(shù)構(gòu)成的序列，輸出序列為輸入序列與該系數(shù)序列的卷積，區(qū)別在確定這個(gè)系數(shù)序列所需的前提條件和準(zhǔn)則不一樣。因?yàn)樾蛄锌梢杂孟蛄縼?lái)表示，所以可以將它們統(tǒng)稱為向量濾波器。而當(dāng)待處理的信號(hào)為已記錄的有限長(zhǎng)的數(shù)據(jù)時(shí)，濾波器可以用一個(gè)矩陣來(lái)實(shí)現(xiàn)，稱之為矩陣濾波器。如果用H表示矩陣濾波器，對(duì)于輸入信號(hào)x，則有輸出信號(hào)為y=Hx。y與期望信號(hào)s的誤差為e=y-s=（H-I）s+Hv，均方誤差為［10］

式（1）中誤差自相關(guān)矩陣Re可表示為

其中Rs和Rv分別為信號(hào)和噪聲的自相關(guān)矩陣，Rsv和Rvs分別為信號(hào)和噪聲的互相關(guān)矩陣。由于EP信號(hào)與自發(fā)腦電是相互獨(dú)立的，并且自發(fā)腦電在短時(shí)間內(nèi)可以近似看作是零均值的平穩(wěn)有色噪聲，所以有Rsv=Rvs=0N×N，因此式（2）可簡(jiǎn)化為

現(xiàn)尋求矩陣Hopt，使E（e2）達(dá)到最小，這樣就可以得到EP信號(hào)在最小均方誤差意義下的最佳估計(jì)=Hoptx。為了得到Hopt，需要先確定濾波矩陣H的基本結(jié)構(gòu)。

若令T=［η1，η2，…，ηN］，由雙正交關(guān)系，可知［ξ1，ξ2，…，ξN］=T-T，則上述投影、加權(quán)和重構(gòu)的過(guò)程可表達(dá)為=T-TDTTx。實(shí)際上信號(hào)和噪聲所在的子空間不可能完全與線性無(wú)關(guān)，總會(huì)有相交的部分，這時(shí)可令D=diag（d），d=［d1，d2，…，dN］T，0≤di≤1。這種方法即為廣義子空間法，由此可得濾波矩陣H的表達(dá)示為

其中T為投影矩陣，且為非奇異陣，其列向量是N維歐氏空間中的一組基；D為系數(shù)加權(quán)矩陣，且為對(duì)角陣，對(duì)投影系數(shù)進(jìn)行加權(quán)；T-T為重構(gòu)矩陣，其列向量是T的列向量的對(duì)偶基，從而組成雙正交基。當(dāng)T是正交陣時(shí)，有T=T-T，則T的列向量就是一組標(biāo)準(zhǔn)正交基。將式（6）帶入式（3）中可得

在式（7）中同時(shí)出現(xiàn)了TTRsT和TTRvT，這相當(dāng)于是用T同時(shí)對(duì)Rs和Rv進(jìn)行合同變換。假設(shè)Rv是正定陣，現(xiàn)尋求投影矩陣Topt，使之可以實(shí)現(xiàn)Rs和Rv的同時(shí)合同對(duì)角化。由于Rs和Rv均為對(duì)稱陣，所以可通過(guò)正交相似變換實(shí)現(xiàn)對(duì)角化，即

令

顯然Rp是對(duì)稱陣，又是非負(fù)定的，所以也可以通過(guò)正交相似變換實(shí)現(xiàn)對(duì)角化，且對(duì)角元素非負(fù)，即

定義投影矩陣Topt為

則Topt可以將Rs和Rv同時(shí)合同對(duì)角化為［12］

將式（7）、（12）代入式（1）中，可得

將E（e2）對(duì)di求偏導(dǎo)并置零，可得

由此可得系數(shù)加權(quán)矩陣Dopt為

在以上算法過(guò)程中，需要知道噪聲自相關(guān)矩陣Rv和信號(hào)自相關(guān)矩陣Rs。Rv可由刺激前記錄的觀測(cè)信號(hào)（不包含誘發(fā)電位成分的純凈的自發(fā)腦電）計(jì)算得到。因?yàn)镽sv=Rvs=0N×N，所以有Rs=Rx-Rv，而觀測(cè)信號(hào)自相關(guān)矩陣Rx可由觀測(cè)信號(hào)計(jì)算得到。

2 仿真實(shí)驗(yàn)與分析

仿真實(shí)驗(yàn)中EP信號(hào)的來(lái)源為通過(guò)疊加平均法得到的貓的EP信號(hào)，采樣頻率為1kHz，單次EP信號(hào)存在時(shí)長(zhǎng)為512ms。使用兩個(gè)周期的EP信號(hào)（相當(dāng)于兩次刺激情況下的記錄）令其為ss，即ss=［s（1），s（2），…，s（512），s（1），s（2），…，s（512）］T，將ss疊加上不同強(qiáng)度的有色噪聲（對(duì)應(yīng)不同的初始信噪比條件）做為觀測(cè)信號(hào)xx=［x（1），x（2），…，x（1024）］T。觀測(cè)信號(hào)自相關(guān)矩陣Rx由樣本自相關(guān)矩陣近似代替，具體方法是首先由xx得到xi=［x（i），x（i+1），…，x（i+511）］T，（i=1，2，…，512），然后形成Hankel矩陣X=［x1，x2，…，x512］512×512，最后估計(jì)出觀測(cè)信號(hào)的自相關(guān)矩陣Rx≈XTX/512。噪聲自相關(guān)矩陣Rv可由不包含在xx中的噪聲數(shù)據(jù)段計(jì)算得到。仿真實(shí)驗(yàn)使用的有色噪聲的AR模型為［13］

其中w（n）為白噪聲。式（16）的AR模型產(chǎn)生的有色噪聲的有效頻帶為0～30Hz，峰值頻率為13Hz，符合自發(fā)腦電的典型特征。而EP信號(hào)的有效頻帶為0～20Hz，峰值頻率為10Hz，有色噪聲在頻帶上完全覆蓋了EP信號(hào)。在仿真實(shí)驗(yàn)中，初始信噪比（SNRi）由0dB逐漸降低到-10dB，衡量實(shí)驗(yàn)結(jié)果的指標(biāo)為估計(jì)結(jié)果的歸一化均方差（σ）、估計(jì)結(jié)果的信噪比（SNRo）和估計(jì)結(jié)果與純凈EP信號(hào)的相關(guān)系數(shù)（r）。

圖1和圖2分別為在SNRi=0dB和SNRi=-10dB條件下不同算法的估計(jì)結(jié)果。從實(shí)驗(yàn)結(jié)果可以看出，廣義子空間法較好地抑制了噪聲的干擾，使信噪比得到了很大程度的提高，在SNRi為-10dB時(shí)仍有較好的效果。在本文算法仿真實(shí)驗(yàn)中，使用了兩個(gè)周期的EP信號(hào)（相當(dāng)于兩次刺激的情況），但估計(jì)結(jié)果比兩次疊加平均的結(jié)果要好得多。

圖1 不同算法在SNRi=0dB條件下估計(jì)結(jié)果的比較。（a）有色噪聲；（b）純凈EP信號(hào)；（c）帶噪信號(hào)；（d）子空間法估計(jì)結(jié)果；（e）兩次疊加平均法估計(jì)結(jié)果；（f）子空間疊加平均法估計(jì)結(jié)果Fig.1 The comparison of the results estimated by different methods for EP signal（SNRi=0dB）.（a）Colored noise；（b）Clean EP signal；（c）Noisy signal；（d）EP signal estimated by subspace method；（e）EP signal estimated by double average method；（f）EP signal estimated by subspace average method

圖3、圖4和圖5分別為隨著SNRi的降低，不同算法運(yùn)行20次平均得到的σ、SNRo和r的變化情況。由圖可見EP信號(hào)估計(jì)效果受制于信噪比。另外EP信號(hào)與噪聲的相關(guān)程度也很重要，這是因?yàn)楫?dāng)相關(guān)程度較大時(shí)，就不能忽略信號(hào)和噪聲的互相關(guān)，盡管EP信號(hào)與噪聲的產(chǎn)生機(jī)理是獨(dú)立的，但在時(shí)域表現(xiàn)上還會(huì)有一定的相關(guān)性。還有一個(gè)因素是噪聲自相關(guān)矩陣和觀測(cè)信號(hào)自相關(guān)矩陣的估計(jì)質(zhì)量，更合理的估計(jì)方法會(huì)得到更好的效果。

3 結(jié)論

圖2 不同算法在SNRi=-10dB條件下估計(jì)結(jié)果的比較。（a）有色噪聲；（b）純凈EP信號(hào)；（c）帶噪信號(hào)；（d）子空間法估計(jì)結(jié)果；（e）兩次疊加平均法估計(jì)結(jié)果；（f）子空間疊加平均法估計(jì)結(jié)果Fig.2 The comparison of the results estimated by different methods for EP signal（SNRi=-10dB）.（a）Colored noise；（b）Clean EP signal；（c）Noisy signal；（d）EP signal estimated by subspace method；（e）EP signal estimated by double average method；（f）EP signal estimated by subspace average method

圖3 隨著SNRi的降低，不同算法（子空間疊加平均法、子空間法、兩次疊加平均法）得到的歸一化均方差σ的變化情況。Fig.3 The relationship between normalized mean squared error σwith SNRifor different algorithms（subspace average method，subspace method and double average method）

本研究提出了一種廣義子空間法用于腦電誘發(fā)電位單導(dǎo)少次提取，其核心是首先利用投影矩陣將信號(hào)和噪聲同時(shí)投影到系數(shù)空間，再根據(jù)觀測(cè)信號(hào)和噪聲的自相關(guān)矩陣得到系數(shù)加權(quán)矩陣，對(duì)系數(shù)進(jìn)行分離估計(jì)出信號(hào)的投影系數(shù)，最后利用重構(gòu)矩陣進(jìn)行重構(gòu)得到期望的EP信號(hào)。該算法較好地抑制了自發(fā)腦電的干擾，使EP信號(hào)的信噪比獲得了較大程度的提高。因?yàn)樾盘?hào)和噪聲的頻譜重疊，所以頻譜分離法難以適用，而廣義子空間法是根據(jù)觀測(cè)信號(hào)構(gòu)造投影矩陣，信號(hào)和噪聲的投影系數(shù)可以同時(shí)非零，這樣使算法的應(yīng)用條件更加廣泛。此外，與疊加平均法相比，該算法最少只需兩次刺激得到的觀測(cè)信號(hào)即可實(shí)現(xiàn)，大大減少了信號(hào)提取所需的周期。理論分析和試驗(yàn)結(jié)果均表明了該算法的有效性。我們將在下一步的研究工作中繼續(xù)完善該算法，并將其應(yīng)用于實(shí)際測(cè)試EP信號(hào)提取問(wèn)題中。

圖4 隨著SNRi的降低，不同算法（子空間疊加平均法、子空間法、兩次疊加平均法）得到的SNRo的變化情況Fig.4 The relationship between SNRowith SNRifor different algorithms（subspace average method，subspace method and double average method）

圖5 隨著SNRi的降低，不同算法（子空間疊加平均法、子空間法、兩次疊加平均法）得到的相關(guān)系數(shù)r的變化情況Fig.5 The relationship between correlation coefficient r with SNRifor different algorithms（subspace average method，subspace method and double average method）

［1］潘映輔.臨床誘發(fā)電位學(xué)（第2版）［M］.北京：人民衛(wèi)生出版社，2000.

［2］Stefanos DG，Perttu OR，Mika PT，et al.Single-trial dynamical estimation of event-related potentials：A kalman filter-based approach［J］.IEEE Transactions on Biomedical Engineering，2005，52（8）：1397-1406.

［3］Qiu Wei，Kenneth SM，F(xiàn)rancis HY，et al.Adaptive filtering of evoked potentials with radial basis function neutral network prefilter［J］.IEEE Transactions on Biomedical Engineering，2002，49（3）：225-231.

［4］Wang Zhisong，Alexander M，David A，et al.Single-trial evoked potential estimation using wavelets［J］.Computers in Biology and Medicine，2007，37（4）：463-473.

［5］Quiroga R，Garcia H.Single-trial event-related potentials with wavelet denoising［J］.Clin Neurophysiol，2003，114：376-390.

［6］畢曉輝，邱天爽，朱勇，等.基于經(jīng)驗(yàn)?zāi)Ｊ椒纸夂酮?dú)立分量分析的單道少次EP信號(hào)提?。跩］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2008，27（6）：817-821.

［7］Aapo H，Erkki O.Independent component analysis：algorithms and applications［J］.Neural Networks，2000，13（4-5）：411-430.

［8］Kumaresan R，Tufts DW.Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise［J］.IEEE Transaction on Acoustics，Speech and Signal Processing，1982，30（6）：833-840.

［9］Overschee PV，Moor BD.Subspace Identification for Linear Systems［M］.Boston：Kluwer Academic Publishers，1996.

［10］Ephraim Y，Trees HLV.A signal subspace approach for speech enhancement［J］.IEEE Transactions on Speech and Audio Processing，1995，3（4）：251-266.

［11］Rezayee A，Gazor S.An adaptive KLT approach for speech enhancement［J］.IEEE Transactions on Speech and Audio Processing，2001，9（2）：87-95.

［12］夏璇.二個(gè)矩陣同時(shí)對(duì)角化［J］.南昌航空工業(yè)學(xué)院學(xué)報(bào)（自然科學(xué)版），2003，17（3）：26-32.

［13］Yu Xiaohu，He Zhenya，Zhang Yisheng.Time-varying adaptive filters for evoked potential estimation［J］.IEEE Transactions on Biomedical Engineering，1994，41（11）：1062-1071.