特征提取是基于P300的腦機接口(BCI)系統中非常關鍵的步驟。獨立分量分析(ICA)算法是效果較好的P300特征提取方法,但目前常用的ICA迭代方法收斂性能均不理想。提出一種基于量子粒子群優化(QPSO)和ICA算法的P300特征提取方法。該方法利用量子計算在計算速度上的優勢,加快了ICA算法的全局收斂,達到了快速有效提取P300的目的。實驗針對BCI Competition Ⅱ dataset Ⅱb和BCI Competition Ⅲ dataset Ⅱ兩組公共數據集進行測試,提取出的P300特征送入線性分類器,系統識別正確率在15次疊加平均情況下達94.4%。實驗結果表明,本文方法用于P300特征提取,在保證提取效果的同時,計算速度更快,為在線BCI系統的進一步研究提供了實驗基礎。
引用本文： 黃璐, 王宏. 基于量子粒子群優化和獨立分量分析的腦電特征提取. 生物醫學工程學雜志, 2014, 31(3): 502-505. doi: 10.7507/1001-5515.20140093 復制
引言
腦機接口(brain-computer interface，BCI)是一種能夠不依賴于外周神經和肌肉組織的交流和控制通路，通過這一通路可以達到利用人腦直接控制外界的目的[1-2]。P300是常用于BCI系統的特征電位，因其無需訓練、無損測量等優點，目前越來越受到學者們的關注[3-5]。當靶刺激出現后300 ms左右，P300出現在被試的頭皮頂區。一般情況下，靶刺激出現的概率越小，誘發的P300電位越明顯[6-7]。
P300是一種誘發電位，可以看成是與自發腦電、眼電偽差等相互獨立的成分，所以，獨立分量分析(independent component analysis，ICA)技術用于P300提取是比較合適的。目前，較常用的ICA迭代方法有梯度法、擬牛頓法，等等[8]。這些方法的求解結果受初始值影響較大，如果初始值選擇不當，很容易收斂于局部最優解，因此導致ICA分離效果不好，分離出的P300電位不明顯。另外，這些算法均是從一點開始搜索，計算時間較長，使用這些算法會導致系統運行緩慢，很難滿足在線BCI系統的實時性要求。
針對上述問題，本文提出一種結合量子粒子群優化(quantum particle swarm optimizer，QPSO)和ICA算法的P300特征提取方法，通過模擬鳥類群體行為的智能活動，多點并行搜索，實現ICA算法的全局收斂，保證了P300特征提取的精確性。并從量子計算角度優化搜索策略，使ICA算法快速收斂，提高了P300特征提取的時間性能。
1 獨立分量分析
ICA技術是一種在源信號和傳輸通道參數均未知的情況下，參照源信號的統計特性，僅根據觀測信號即可求得各獨立源信號的方法[9]。設m個相互獨立的源信號s(t)=(s1(t),s2(t),…,sm(t))T經過線性混合，即
| $x\left( t \right)=As\left( t \right),$ |
得到n維觀測信號x(t)=(x1(t),x2(t),…,xn(t))T，其中A為n×m維線性混合矩陣。則存在一個解混矩陣W，使得
| $y\left( t \right)=x\left( t \right),$ |
ICA需要解決的問題是在混合矩陣A未知的情況下找到W的最佳估計，使得估計信號y(t)是對源信號s(t)的有效逼近，即
| $s\left( t \right)=\bar{W}x\left( t \right),$ |
式中是A的逆矩陣，y(t)=(y1(t),y2(t),…,ym(t))T，且y(t)中各分量是相互獨立的。
目前，已提出很多不同的目標函數可以度量估計信號y(t)中各分量之間的獨立性，本文采用負熵最大法則，建立如下目標函數
| $\begin{align} & Maximize:~C\left( y \right)=\sum\limits_{i=1}^{m}{J\left( {{y}_{i}} \right)}, \\ & Subject\text{ }to:~\left| \left| W \right| \right|=0 \\ \end{align}$ |
其中 yi=wix(t)，J(yi)≈ρ (E {Gi (yi)}-E {Gi (v)}) 2，v是均值為0、方差為1的高斯變量，ρ是正常數，E{·}是均值函數，Gi (·)是非二次函數。
2 量子粒子群優化算法
粒子群優化(particle swarm optimizer，PSO)算法是一種群體智能優化算法，源于對鳥類群體活動的簡單模擬[10]。PSO算法可以有效地優化多種功能。在PSO算法中，優化問題的解決方案可以看作是由幾個粒子構成的群體在一個D維搜索空間中飛行，在每次迭代中，每一個粒子根據個體最優位置和全局最優位置這兩個因素調整其自身運動方向。標準PSO算法不能保證全局收斂性，它經常因為迭代早熟而陷入局部最優解。
2004年，Sun等[11]從量子力學的角度出發，提出了一種新的PSO算法，即QPSO算法。在QPSO算法中，演化方程不再需要速度矢量，形式更簡單，更容易控制，因此QPSO算法的搜索能力比PSO算法更強。QPSO算法具體的迭代公式為
| $p=\frac{({{\alpha }_{1}}\times {{p}_{i}}+{{\alpha }_{2}}\times {{p}_{g}})}{({{\alpha }_{1}}+{{\alpha }_{2}})},$ |
| $mbest=\frac{1}{M}\sum\limits_{i=1}^{M}{{{p}_{i}},}~$ |
| $b\left( k+1 \right)=p\pm \beta \times |mbest-b\left( k \right)|\times ln(\frac{1}{\mu })~,$ |
式中 M為種群大小，pi為第i個粒子的個體歷史最優位置，pg為粒子群的全局最優位置，p為粒子的局部最優位置，mbest為粒子群中所有粒子的局部最優位置的平均值,即所有粒子的重心，α1、α2和μ是在0~1之間均勻分布的隨機數，b(k)是k次迭代后粒子的位置，β為收縮擴展因子，它控制算法的收斂速度。一般情況下，k次迭代后，β可以用如下公式表示，即
| $\beta =0.5+0.5\times \frac{Ite{{r}_{max}}-k}{Ite{{r}_{max}}},$ |
式中 Itermax為最大迭代次數，k為粒子的當前迭代次數。
3 基于QPSO的ICA算法
ICA實際上是一個優化問題，包含兩個因素：目標函數和尋優算法。在ICA的尋優環節中引入QPSO算法，就構成了基于QPSO的ICA算法。目標函數為公式(4)所示，迭代完成后，粒子群的全局最優位置即為解混矩陣W的最優解。具體步驟如下：
(1)預處理觀測信號x(t)，包括中心化和白化，得到z(t)；
(2)初始化粒子群向量空間W(0)和所有相關參數；
(3)根據公式(4)計算各粒子的適應度值，其中，y(t)=Wz(t)；
(4)對于每一個粒子，比較其當前位置與其個體歷史最優位置pi處的適應度值，如果當前位置的適應度比pi處的適應度高，則用當前位置替換pi；
(5)計算粒子群中所有粒子的個體歷史最優位置，將所有粒子群的全局最優位置賦給pg；
(6)根據公式(5)、(6)和(7)計算得到mbest、 p和粒子位置向量b；
(7)判斷粒子適應度是否滿足收斂條件或者迭代次數是否達到最大值。如果不是，回到步驟(3)，否則轉到步驟(8)；
(8)得到混合矩陣和相應的估計信號y(t)。
4 實驗及結果分析
本實驗采用BCI Competition Ⅱ的dataset Ⅱb[12]和BCI Competition Ⅲ的dataset Ⅱ[13]兩組公共數據集對提出的方法進行驗證。這兩個數據集的腦電圖(electroencephyalogram,EEG)采集均基于Farwell等[14]于1988年提出的speller實驗范式，數據采集平臺為BCI 2000實驗平臺。
實驗過程中，被試的任務是盯住如圖 1所示的屏幕上方提供的單詞中某一個字符，即目標靶。實驗開始后，字符行/列隨機閃爍，每行每列均閃爍1次，共計12次，稱為1個trial，其中包含2次目標靶刺激，10次非靶刺激。針對每1個字符，這樣的trial共重復15次。BCI Competition Ⅱ dataset Ⅱb 數據采集自一個被試Subject Ⅱ，BCI Competition Ⅲ dataset Ⅱ數據采集自兩個被試Subject ⅢA和Subject ⅢB，采樣頻率為240 Hz，導聯位置采用國際10-20標準。
本實驗基于MATLAB R2011b運行平臺，對觀測EEG數據進行0.1~12 Hz濾波，濾波器階數設為24階。對訓練集EEG數據，采用基于QPSO的ICA算法，以公式(4)為目標函數做ICA計算。然后，將計算出的解混矩陣直接左乘經過預處理后的測試集EEG數據，得到獨立分量。將相同刺激行/列的獨立分量進行0~650 ms分段并疊加平均，順次連接在一起，選擇刺激后250~400 ms內振幅最大，且在頭皮頂部電極位置Cz、C1、C2處振幅最大或次大者保留，其余分量置零，這樣就保留了大部分含有P300信息的獨立分量。將保留下來的獨立分量反算到頭皮導聯空間。圖 2給出了提取的Subject ⅢA的15次疊加平均P300波形。
圖1
P300拼寫器刺激界面
Figure1.
The user displagy in P300 speller
圖2
Subject ⅢA平均P300波形
Figure2.
Average P300 waveform of Subject ⅢA
采用性能指標PI衡量算法的分離性能[15]，即
| $\begin{align} & PI=\frac{1}{2M}\sum\limits_{i=1}^{M}{\left\{ \left( \sum\limits_{j=1}^{M}{\frac{|{{c}_{ij}}{{|}^{2}}}{\underset{k}{\mathop{max}}\,|{{c}_{ik}}{{|}^{2}}}-1} \right) \right.}+ \\ & \left. \left( \sum\limits_{j=1}^{M}{\frac{|{{c}_{ij}}{{|}^{2}}}{\underset{k}{\mathop{max}}\,|{{c}_{ik}}{{|}^{2}}}-1} \right) \right\} \\ \end{align}$ |
式中 C=WA= (cij),M為變量的個數，求得的PI 指數越小，說明分離效果越好。當PI=0時，表示信號被完全分離。為了驗證算法的有效性，分別與Infomax ICA算法和基于PSO的ICA (PSO_ICA)算法進行了EEG分離結果對比，如表 1所示，表中的迭代次數和迭代時間針對一個字符EEG數據而言。
表1
分離性能分析
Table1.
Analysis of the separation performance
從表 1可以看出，采用PSO_ICA算法進行P300特征提取，由于PSO算法本身的多點并行搜索特性，其迭代次數和迭代時間比Infomax ICA算法明顯降低。但是，由于標準PSO算法不能保證全局收斂，因此，PI指標并不十分理想。而本文采用基于QPSO 的ICA 算法進行P300提取，因為結合了量子計算在計算速度上的優勢，迭代次數和迭代時間較前兩種方法更為減少，達到最優。同時，采用本文算法與采用Infomax ICA算法進行EEG分離所得的PI性能指標接近,說明本文算法較好地分離了實驗所使用EEG數據集中的各獨立分量，適合P300電位的提取。
在此基礎上，提取Cz通道EEG在250~400 ms時間窗內的峰振幅、波形面積和50%面積潛伏期，構成三維特征向量送入線性分類器進行分類。針對兩個數據集中所有被試的平均識別正確率如表 2所示。
表2
識別正確率
Table2.
Recognition accuracy
從表 2可以看出，應用本文算法提取出的P300特征進行BCI分類，系統識別正確率基本上與應用Infomax ICA算法的分類效果相同，進一步驗證了本文算法進行P300特征提取的有效性。
5 結論
在基于P300的BCI系統中，P300的有效提取是順利進行下一步分類的保證。因為用于BCI分類的EEG數據量龐大，且P300振幅遠遠低于自發腦電，目前的P300特征提取算法均耗時較多，很難滿足在線BCI系統的需求。因此，如何快速有效地提取P300成為重要的研究課題。本文采用基于QPSO和ICA的P300特征提取方法，利用QPSO算法的全局收斂性和量子計算的高效性，快速實現了EEG數據的獨立分量分解，為后續基于P300的BCI系統分類奠定了良好的基礎。實驗采用BCI Competition Ⅱ和Ⅲ的兩組公共數據集進行了方法驗證，并從PI、迭代次數和迭代時間等指標考察了算法性能。實驗結果表明，基于QPSO和ICA的特征提取方法收斂速度更快，提取P300效果理想，為在線BCI系統的進一步研究提供了更好的方法。
引言
腦機接口(brain-computer interface，BCI)是一種能夠不依賴于外周神經和肌肉組織的交流和控制通路，通過這一通路可以達到利用人腦直接控制外界的目的[1-2]。P300是常用于BCI系統的特征電位，因其無需訓練、無損測量等優點，目前越來越受到學者們的關注[3-5]。當靶刺激出現后300 ms左右，P300出現在被試的頭皮頂區。一般情況下，靶刺激出現的概率越小，誘發的P300電位越明顯[6-7]。
P300是一種誘發電位，可以看成是與自發腦電、眼電偽差等相互獨立的成分，所以，獨立分量分析(independent component analysis，ICA)技術用于P300提取是比較合適的。目前，較常用的ICA迭代方法有梯度法、擬牛頓法，等等[8]。這些方法的求解結果受初始值影響較大，如果初始值選擇不當，很容易收斂于局部最優解，因此導致ICA分離效果不好，分離出的P300電位不明顯。另外，這些算法均是從一點開始搜索，計算時間較長，使用這些算法會導致系統運行緩慢，很難滿足在線BCI系統的實時性要求。
針對上述問題，本文提出一種結合量子粒子群優化(quantum particle swarm optimizer，QPSO)和ICA算法的P300特征提取方法，通過模擬鳥類群體行為的智能活動，多點并行搜索，實現ICA算法的全局收斂，保證了P300特征提取的精確性。并從量子計算角度優化搜索策略，使ICA算法快速收斂，提高了P300特征提取的時間性能。
1 獨立分量分析
ICA技術是一種在源信號和傳輸通道參數均未知的情況下，參照源信號的統計特性，僅根據觀測信號即可求得各獨立源信號的方法[9]。設m個相互獨立的源信號s(t)=(s1(t),s2(t),…,sm(t))T經過線性混合，即
| $x\left( t \right)=As\left( t \right),$ |
得到n維觀測信號x(t)=(x1(t),x2(t),…,xn(t))T，其中A為n×m維線性混合矩陣。則存在一個解混矩陣W，使得
| $y\left( t \right)=x\left( t \right),$ |
ICA需要解決的問題是在混合矩陣A未知的情況下找到W的最佳估計，使得估計信號y(t)是對源信號s(t)的有效逼近，即
| $s\left( t \right)=\bar{W}x\left( t \right),$ |
式中是A的逆矩陣，y(t)=(y1(t),y2(t),…,ym(t))T，且y(t)中各分量是相互獨立的。
目前，已提出很多不同的目標函數可以度量估計信號y(t)中各分量之間的獨立性，本文采用負熵最大法則，建立如下目標函數
| $\begin{align} & Maximize:~C\left( y \right)=\sum\limits_{i=1}^{m}{J\left( {{y}_{i}} \right)}, \\ & Subject\text{ }to:~\left| \left| W \right| \right|=0 \\ \end{align}$ |
其中 yi=wix(t)，J(yi)≈ρ (E {Gi (yi)}-E {Gi (v)}) 2，v是均值為0、方差為1的高斯變量，ρ是正常數，E{·}是均值函數，Gi (·)是非二次函數。
2 量子粒子群優化算法
粒子群優化(particle swarm optimizer，PSO)算法是一種群體智能優化算法，源于對鳥類群體活動的簡單模擬[10]。PSO算法可以有效地優化多種功能。在PSO算法中，優化問題的解決方案可以看作是由幾個粒子構成的群體在一個D維搜索空間中飛行，在每次迭代中，每一個粒子根據個體最優位置和全局最優位置這兩個因素調整其自身運動方向。標準PSO算法不能保證全局收斂性，它經常因為迭代早熟而陷入局部最優解。
2004年，Sun等[11]從量子力學的角度出發，提出了一種新的PSO算法，即QPSO算法。在QPSO算法中，演化方程不再需要速度矢量，形式更簡單，更容易控制，因此QPSO算法的搜索能力比PSO算法更強。QPSO算法具體的迭代公式為
| $p=\frac{({{\alpha }_{1}}\times {{p}_{i}}+{{\alpha }_{2}}\times {{p}_{g}})}{({{\alpha }_{1}}+{{\alpha }_{2}})},$ |
| $mbest=\frac{1}{M}\sum\limits_{i=1}^{M}{{{p}_{i}},}~$ |
| $b\left( k+1 \right)=p\pm \beta \times |mbest-b\left( k \right)|\times ln(\frac{1}{\mu })~,$ |
式中 M為種群大小，pi為第i個粒子的個體歷史最優位置，pg為粒子群的全局最優位置，p為粒子的局部最優位置，mbest為粒子群中所有粒子的局部最優位置的平均值,即所有粒子的重心，α1、α2和μ是在0~1之間均勻分布的隨機數，b(k)是k次迭代后粒子的位置，β為收縮擴展因子，它控制算法的收斂速度。一般情況下，k次迭代后，β可以用如下公式表示，即
| $\beta =0.5+0.5\times \frac{Ite{{r}_{max}}-k}{Ite{{r}_{max}}},$ |
式中 Itermax為最大迭代次數，k為粒子的當前迭代次數。
3 基于QPSO的ICA算法
ICA實際上是一個優化問題，包含兩個因素：目標函數和尋優算法。在ICA的尋優環節中引入QPSO算法，就構成了基于QPSO的ICA算法。目標函數為公式(4)所示，迭代完成后，粒子群的全局最優位置即為解混矩陣W的最優解。具體步驟如下：
(1)預處理觀測信號x(t)，包括中心化和白化，得到z(t)；
(2)初始化粒子群向量空間W(0)和所有相關參數；
(3)根據公式(4)計算各粒子的適應度值，其中，y(t)=Wz(t)；
(4)對于每一個粒子，比較其當前位置與其個體歷史最優位置pi處的適應度值，如果當前位置的適應度比pi處的適應度高，則用當前位置替換pi；
(5)計算粒子群中所有粒子的個體歷史最優位置，將所有粒子群的全局最優位置賦給pg；
(6)根據公式(5)、(6)和(7)計算得到mbest、 p和粒子位置向量b；
(7)判斷粒子適應度是否滿足收斂條件或者迭代次數是否達到最大值。如果不是，回到步驟(3)，否則轉到步驟(8)；
(8)得到混合矩陣和相應的估計信號y(t)。
4 實驗及結果分析
本實驗采用BCI Competition Ⅱ的dataset Ⅱb[12]和BCI Competition Ⅲ的dataset Ⅱ[13]兩組公共數據集對提出的方法進行驗證。這兩個數據集的腦電圖(electroencephyalogram,EEG)采集均基于Farwell等[14]于1988年提出的speller實驗范式，數據采集平臺為BCI 2000實驗平臺。
實驗過程中，被試的任務是盯住如圖 1所示的屏幕上方提供的單詞中某一個字符，即目標靶。實驗開始后，字符行/列隨機閃爍，每行每列均閃爍1次，共計12次，稱為1個trial，其中包含2次目標靶刺激，10次非靶刺激。針對每1個字符，這樣的trial共重復15次。BCI Competition Ⅱ dataset Ⅱb 數據采集自一個被試Subject Ⅱ，BCI Competition Ⅲ dataset Ⅱ數據采集自兩個被試Subject ⅢA和Subject ⅢB，采樣頻率為240 Hz，導聯位置采用國際10-20標準。
本實驗基于MATLAB R2011b運行平臺，對觀測EEG數據進行0.1~12 Hz濾波，濾波器階數設為24階。對訓練集EEG數據，采用基于QPSO的ICA算法，以公式(4)為目標函數做ICA計算。然后，將計算出的解混矩陣直接左乘經過預處理后的測試集EEG數據，得到獨立分量。將相同刺激行/列的獨立分量進行0~650 ms分段并疊加平均，順次連接在一起，選擇刺激后250~400 ms內振幅最大，且在頭皮頂部電極位置Cz、C1、C2處振幅最大或次大者保留，其余分量置零，這樣就保留了大部分含有P300信息的獨立分量。將保留下來的獨立分量反算到頭皮導聯空間。圖 2給出了提取的Subject ⅢA的15次疊加平均P300波形。
圖1
P300拼寫器刺激界面
Figure1.
The user displagy in P300 speller
圖2
Subject ⅢA平均P300波形
Figure2.
Average P300 waveform of Subject ⅢA
采用性能指標PI衡量算法的分離性能[15]，即
| $\begin{align} & PI=\frac{1}{2M}\sum\limits_{i=1}^{M}{\left\{ \left( \sum\limits_{j=1}^{M}{\frac{|{{c}_{ij}}{{|}^{2}}}{\underset{k}{\mathop{max}}\,|{{c}_{ik}}{{|}^{2}}}-1} \right) \right.}+ \\ & \left. \left( \sum\limits_{j=1}^{M}{\frac{|{{c}_{ij}}{{|}^{2}}}{\underset{k}{\mathop{max}}\,|{{c}_{ik}}{{|}^{2}}}-1} \right) \right\} \\ \end{align}$ |
式中 C=WA= (cij),M為變量的個數，求得的PI 指數越小，說明分離效果越好。當PI=0時，表示信號被完全分離。為了驗證算法的有效性，分別與Infomax ICA算法和基于PSO的ICA (PSO_ICA)算法進行了EEG分離結果對比，如表 1所示，表中的迭代次數和迭代時間針對一個字符EEG數據而言。
表1
分離性能分析
Table1.
Analysis of the separation performance
從表 1可以看出，采用PSO_ICA算法進行P300特征提取，由于PSO算法本身的多點并行搜索特性，其迭代次數和迭代時間比Infomax ICA算法明顯降低。但是，由于標準PSO算法不能保證全局收斂，因此，PI指標并不十分理想。而本文采用基于QPSO 的ICA 算法進行P300提取，因為結合了量子計算在計算速度上的優勢，迭代次數和迭代時間較前兩種方法更為減少，達到最優。同時，采用本文算法與采用Infomax ICA算法進行EEG分離所得的PI性能指標接近,說明本文算法較好地分離了實驗所使用EEG數據集中的各獨立分量，適合P300電位的提取。
在此基礎上，提取Cz通道EEG在250~400 ms時間窗內的峰振幅、波形面積和50%面積潛伏期，構成三維特征向量送入線性分類器進行分類。針對兩個數據集中所有被試的平均識別正確率如表 2所示。
表2
識別正確率
Table2.
Recognition accuracy
從表 2可以看出，應用本文算法提取出的P300特征進行BCI分類，系統識別正確率基本上與應用Infomax ICA算法的分類效果相同，進一步驗證了本文算法進行P300特征提取的有效性。
5 結論
在基于P300的BCI系統中，P300的有效提取是順利進行下一步分類的保證。因為用于BCI分類的EEG數據量龐大，且P300振幅遠遠低于自發腦電，目前的P300特征提取算法均耗時較多，很難滿足在線BCI系統的需求。因此，如何快速有效地提取P300成為重要的研究課題。本文采用基于QPSO和ICA的P300特征提取方法，利用QPSO算法的全局收斂性和量子計算的高效性，快速實現了EEG數據的獨立分量分解，為后續基于P300的BCI系統分類奠定了良好的基礎。實驗采用BCI Competition Ⅱ和Ⅲ的兩組公共數據集進行了方法驗證，并從PI、迭代次數和迭代時間等指標考察了算法性能。實驗結果表明，基于QPSO和ICA的特征提取方法收斂速度更快，提取P300效果理想，為在線BCI系統的進一步研究提供了更好的方法。
