趙蕾,白雪梅,胡超
(長春理工大學 電子信息工程學院,長春 130022)
酒精是親脂性物質,長期飲用可引起多種神經損傷和精神障礙,包括依賴、戒斷綜合征及精神性疾病癥狀。與短時大量飲酒造成的急性酒精中毒癥狀不同,長期酗酒引發的神經損傷造成的酒精使用障礙不易被察覺,極易給從事高風險工作的人群帶來安全隱患。因此如何更準確、更快速的檢測出長期酗酒所造成的神經損傷的研究是很有必要的。
腦電圖(Electroencephalogram,EEG)信號可以作為長期酗酒造成的腦神經損傷的診斷依據。診斷從患者觀察和數據收集開始,去除腦電噪聲,提取有效特征作為輸入信號送到各分類算法模型中,達到智能輔助診斷的目的[1-2]。而腦電圖的非線性、非平穩性、非周期性的特征使得腦電信號的去噪與特征提取面臨很多困難。由于腦電信號是包含大量腦電信息的多電極數據,且極易受到各種噪聲的影響,信號能量分布分散,而酒精腦電相對于正常腦電變化更復雜,如何在預留信號原始信息的基礎上盡量減少數據分析的計算量,在去除龐大的腦電圖數據中的噪聲的同時又能保留有效數據和主要特征是研究的目的。
EEG數據來源于紐約州立大學健康中心Henri Begleiter神經動力學實驗室。數據集包含了放置在受試者頭皮上的64個電極的測量數據,這些電極在頭皮處的采樣頻率為256 Hz,采樣時間為1 s。
受試者被分為兩組,有長期酗酒史的酒精組和無飲酒史的對照組。受試者受到90張來自Snodgrass和Vanderwart圖片集合的物體圖片的刺激,在實驗過程中兩張圖片刺激交替出現,每張圖片刺激持續時間為300 ms,受試者需要判定兩個刺激是否相同,以此獲得被測試者的事件相關電位(ERP)[3]。每個受試者都進行了120次試驗。原始64通道酒精腦電數據示例如圖1所示。

圖1 64通道原始酒精腦電數據
在腦電信號的研究應用中,為了描述大腦不同區域對刺激的反映,需要對大腦不同區域多個通道進行數據采集[4]。多通道的腦電數據樣本給研究和識別應用提供了大量的信息,但是也給后續的數據處理和識別工作增加了工作量。而在酗酒人群和正常對照組的腦電分析研究中,不同通道的腦電信號之間可能存在相關性,如果分別對每個通道單獨分析,得到的結果是孤立的,盲目的減少指標又會丟失有用信息,可能得到錯誤的結論。因此,需要找到一種在減少分析數據的指標同時最多的保留原始數據信息的方法,以達到對數據進行降維[5]的同時也能對數據進行全面分析的目的。
主成分分析法(PCA)就是這樣一種把多個特征映射成少數幾個綜合特征的統計分析降維方法,它是模式識別中數據縮減、圖像壓縮和特征提取的重要技術。
PCA的主要思想是將n維特征映射到k維上,這個k維是全新的正交特征也叫做主成分,就是以原有n維特征為基礎提取出來的k維特征。PCA的原理就是從原始的空間中順序地找一組相互正交的坐標軸,新的坐標軸的選擇與數據本身是密切相關的。第一個新坐標軸選擇是原始數據中方差最大的方向,第二個新坐標軸選取的是與第一個坐標軸正交的平面中方差最大的方向,第三個軸是與第一、二個軸正交的平面中方差最大的方向以此類推,得到了n個這樣的坐標軸。大部分方差都是包含在前面的k個坐標軸中,后面的坐標軸方差幾乎為0,因此,后面的坐標軸可以被忽略只保留前面k個坐標軸。通過計算數據矩陣的協方差矩陣,得到協方差矩陣的特征值與特征向量,選擇特征值最大(即方差最大)的k個特征對應的特征向量組成的矩陣,從而實現數據降維的目的。
PCA算法的主要步驟如下:
(1)樣本集為m維隨機向量組成的樣本集D={x1,x2,...,xm},其中n個樣本構造樣本矩陣,Di={xi1,xi2,...,xim}T,i=1,2,...,n對樣本陣進行標準變化:

得到標準化矩陣A。
(2)計算樣本的協方差矩陣AAT;
(3)對協方差矩陣AAT矩陣做特征值分解,得到協方差矩陣的特征值λ1,λ2,...,λn(i=1,2,...,n)。
(4)特征值所對應的特征向量為b1,b2,...,bn;得到投影矩陣B=(b1,b2,...,bn)T。
(5)令Y=(y1,y2,...,yn)T,Y=BTA,且有BTAATB=Λ,其中,y1,y2,...,yn彼此不相關,稱y1,y2,...,yn分別為第1、第2、...、第i個主分量,的值稱為主分量的貢獻率。前k個主分量的累積貢獻率定義為:

選取前k個主分量代替原始數據做分析,這樣就可以達到降低原始數據維數的目的。
本研究中的分析變量是64個電極通道,將采集到的腦電信號經過簡單濾波后處理成如下矩陣:

首先采用PCA算法對64通道腦電數據進行降維處理。可得到主成分貢獻率如圖2所示。

圖2 各主成分貢獻率
主成分累積貢獻率如表1所示。

表1 前k個主成分累積貢獻率(%)
一般累計貢獻率達到80%~90%之間[6],則說明這些主成分保留了原始數據的主要信息,因此,保留前十五個主成分進行分析,將數據維數從64維降到了15維,實現了數據壓縮降維。降維后前15個主成分系數如圖3所示。

圖3 前十五個主成分系數
對比圖1和圖3,經過降維后的數據與原始數據特征一致,數據主要特征被很好的保留下來,而數據維數減少,后續的數據處理工作量也隨之減小。
經過PCA降維,酒精腦電數據中的部分冗余特征被去除,但是腦電信號是一種非線性、非周期性的非平穩信號,極易受到噪聲干擾,因此經過降維提取的主成分包含噪聲信號。腦電信號中包含的主要噪聲就是腦電信號偽跡,腦電信號偽跡去除的最大困難在于人體自身偽跡信號的去除(如眼電偽跡等)[7]。自身偽跡與環境噪聲不同,它易與腦電信號重疊,無法通過常規的濾波手段去除,現在常用的腦電信號去噪方法普遍需要同時采集偽跡信號數據作為參考進行去噪的工作,而許多腦電實驗無法使用眼電電極等檢測手段進行信號檢測作為參考。如何在沒有參考數據的情況下去除噪聲并盡可能多的保留原始腦電數據的特征是腦電信號研究要解決的一個問題。
目前有許多信號領域內經常應用的去噪方法例如快速獨立成分分析法(FastICA)等方法被應用于腦電去噪的領域,但是由于腦電信號的特殊性,這些去噪方法并不完美。
取降維后得到的第一個主成分如圖4所示。

圖4 第一主成分
因所用腦電數據無先驗參考信號,所以考慮使用FastICA對腦電數據進行去噪處理。以第一主成分為例,應用FastICA對酒精腦電進行去噪,得到去噪前后效果對比如圖5所示。

圖5 FastICA去噪前后效果對比
從圖5所示的結果可以看出,針對本文采用的實驗數據,由于采樣點數較少,FastICA雖然去除了腦電中的偽跡噪聲,但是原有數據特征也損失嚴重,去噪后的腦電信號只保留了較少的數據信息。且這種方法對初始值的選擇比較敏感,計算量大,分選時間長,對信噪比較低的信號估計精度也偏低。
奇異值分解(SVD)[8]是一種用于信號去噪的非線性濾波方法。這種去噪技術是空間算法的一種,經過SVD可以將帶有噪聲信號的向量空間分解純凈信號主導的子空間和噪聲信號主導的子空間,通過去除噪聲空間中的帶噪聲信號分量來估計純凈信號,適用于無先驗參考的腦電信號的去噪。
有一采集到的含噪聲原始信號X={x(i)},i=1,2,3,...,N,其中N是信號的長度,首先構造Hankel矩陣,形式如下:

式中,m=N-n+1,1<n<N,n的取值范圍一般為N/20~N/2。實際應用中一般令m≥n,存在A∈Rm×n。對A做奇異值分解,即一定存在正交矩陣U∈Rm×m和正交矩陣V∈Rn×n,使得

式中,Σ∈Rm×n為對角矩陣,Σ=diag(σ1,σ2,...,σp),p=min(m,n)。對角線元素σ1,σ2,...,σp稱為A的奇異值并按降序排列,其中包含了A的秩的特性信息。分別將U和V的列向量組成正交基空間,那么A可以認為是由多個分量的線性疊加得到的,每一個分量用Di代表,且每一個分量都可以用對應的奇異值表示出來,描述如下:

式中,ui=Rm×1,vi∈Rn×1,i=1,2,...,p。每個奇異值對應著一個不同的信號分量,并且奇異值滿足σ1≥σ2≥...≥σr≥σr+1≥...≥σp≥0,前面r個奇異值代表了有用的信號分量,后面p-r個奇異值則對應著噪聲信號分量,r即為需要保留的矩陣的有效秩的階次,從而得到重構的去噪后的純凈信號。
依然以酒精腦電信號經過PCA降維后得到的第一個主成分為例,構造Hankel矩陣。

對Hankel矩陣Yi進行奇異值分解,得到奇異值分布如圖6所示。

圖6 第一個主成分矩陣的奇異值
想要得到的重構信號既盡可能地保留原始數據的特征,又要去除噪聲干擾,就要合理地保留奇異值。因此,如何確定有效秩的階次是SVD的關鍵。選取不同有效秩階次對信號去噪的效果是存在著明顯的差異的,階次選擇過大會導致信號去噪不徹底,去噪后的信號中仍然存在噪聲;階次選擇過小,會導致信號信息缺失,信號特征提取不全。
對于這個問題,參考文獻[9]提出一種方法,利用SVD本身的奇異性檢測能力達到SVD信號去噪中有效階次確定的目的,然而經過實驗得到分量信號D2如圖7所示。

圖7 奇異值序列分量信號D2
實驗所得腦電信號數據沒有明顯的頻率特征和周期特點,很難準確的判斷信號分量中奇異點的位置,這種方法對于采樣信號頻率小、采樣點數少的腦電信號數據集并不適用。
在參考文獻[10]中研究了信號與噪聲的奇異值特征,文中指出正常信號的奇異值有突變而噪聲信號的奇異值相對平穩。文獻中提出了奇異值差分譜的概念,它由奇異值序列的前向差分組成,可以分辨出復雜信號奇異值的突變狀態,奇異值的有效秩的階次選擇可以通過差分譜的峰值來實現。

圖8 奇異值差分譜
通過圖8可以看出奇異值差分譜的最大峰值位于第二個坐標,所以有效秩的階次選擇為2,按照差分譜選擇奇異值的原則,保留前兩個奇異值重構后的信號如圖9所示。

圖9 奇異值差分譜法重構信號
通過圖像可以看出,奇異值差分譜法克服了有限采樣點數和噪聲帶來的問題,在去除噪聲的同時比FastICA法保留了更多原始數據的特征,但是數據幅值壓縮過大,波形細節丟失的問題仍然很嚴重。
視覺刺激產生事件相關電位(ERP)的能量分布比較集中,過大的奇異值數列將導致信號分量的能量分散,使原始信號數據特征細節丟失。因此,重新構造Hankel矩陣,為了使奇異值分解后得到的奇異值序列盡可能的小,令行數為255,列數為2,得到Hankel矩陣如下:

對Hankel矩陣進行奇異值分解,選擇第一個奇異值為有效奇異值,令其它奇異值為零,得到重構奇異值矩陣,然后重構信號序列獲得重構信號圖10。

圖10 重構第一個主成分序列
圖10與未經處理的主成分圖4相比,部分噪聲被去除。與FastICA和奇異值差譜法去噪相比波形的更多細節得以保留。這樣去噪的弊端是重構后的波形仍然含有噪聲。繼續對重構后的信號序列進行奇異值去噪處理。將第一次重構后的序列作為第二次SVD的輸入,如此反復,進行奇異值分解的迭代過程直到得到純凈信號,具體流程如圖11所示。

圖11 SVD迭代流程示意圖
如圖12所示,經過四次SVD處理去噪后的主成分相比于經過一次SVD處理波形平滑了很多,很明顯更多的噪聲被去除了,而第一主成分的特征信息也被很好的保留下來。數據處理后的圖像表明,對于本實驗的腦電數據,經過四次SVD處理去噪后已經得到比較理想的效果,繼續重復SVD處理重構信號序列,圖像波形基本無變化。

圖12 多次SVD后重構第一個主成分序列
同樣,對其它14個主成分也進行多次重復SVD處理重構主成分信號序列。
傳統的非線性信號處理方法對采樣時長短、采樣點數少的酒精腦電信號數據進行去噪時,不能保留原始數據結構特點,信息損失較多。針對本文采用的腦電信號的數據特點,充分利用了PCA方法可以提取主成分的特點和奇異值良好的數學特性,將PCA與SVD多次迭代方法相結合,對實驗數據進行降維去噪處理,成功將64通道的腦電信號降至15維,經過改進的SVD多次迭代處理,對于非線性、非周期性的酒精腦電信號去噪效果明顯,計算量減少,且原有的信號特征也被更好的保留了下來。