穩態視覺誘發電位(SSVEP)是由持續的視覺刺激而誘發的節律性腦電信號。SSVEP頻率由固定的視覺刺激頻率及其諧波頻率組成。二維集合經驗模式分解(2D-EEMD)是經典的經驗模式分解算法的改進算法, 將分解拓展到二維方向上。本文首創性地將2D-EEMD應用于SSVEP。分解得到的本征模式函數(IMF)的二維圖像可清晰地觀測到SSVEP頻率。經過噪聲和偽跡濾除的SSVEP主要有效IMF成分投影到頭圖上, 可以反映大腦對視覺刺激的時變趨勢, 以及大腦不同區域的反應程度, 結果顯示枕葉區對于視覺刺激的反應最為強烈。最后本文用短時傅里葉變換(STFT)對2D-EEMD的重構信號進行SSVEP頻率提取, 其識別準確率提高了16%。
引用本文： 楊晨, 黃麗亞, 文念, 楊俊宇. 基于二維集合經驗模式分解的穩態視覺誘發電位目標檢測研究. 生物醫學工程學雜志, 2015, 32(3): 508-513. doi: 10.7507/1001-5515.20150093 復制
引言
腦機接口(brain computer interface, BCI)是通過腦電圖(electroencephalogram，EEG)信號的轉化而建立起來的大腦與外設之間的控制系統，它不依賴于常規的腦信號輸出通路——外周神經及肌肉組織，可以為癱瘓或者神經肌肉受損的殘疾人提供全新的通信控制方式[1-2]。基于EEG的BCI系統憑借其非侵入式檢測所帶來的低傷害性和便捷性，在當今的BCI系統中得到了廣泛的應用。其中，穩態視覺誘發電位(steady state visual evoked potential，SSVEP)在近幾年的研究中更是引起了廣泛的關注。SSVEP是由持續視覺刺激而誘發的節律性EEG信號[3]，SSVEP頻率由固定的視覺刺激頻率及其諧波頻率組成[4]，因此通過研究SSVEP可以準確提取出視覺刺激頻率。將不同的任務和不同的刺激頻率相組合，可以真正實現人機交互。
經驗模式分解(empirical mode decomposition，EMD)是希爾伯特-黃變換(Hilbert-Huang transform，HHT)的預處理過程[5]，EMD分解以信號自身為基礎，通過反復篩選得到的本征模態函數(intrinsic mode function，IMF)是自適應的，適合處理像腦電一類的非平穩非線性信號[6]。集合經驗模式分解(ensemble EMD, EEMD)作為一種新的噪聲輔助分析法，在一定程度上緩解了EMD所帶來的模態混疊問題[7-10]。而二維集合經驗模式分解(two-dimensional EEMD, 2D-EEMD)則在兩個正交方向上分別獨立地對信號進行EEMD分解，然后根據“可比較的小尺度組合準則”對分解結果進行重新組合，得到最終的二維分解結果[11]。2D-EEMD已廣泛應用在海洋颶風、核磁共振成像(nuclear magnetic resonance imaging，MRI)、Lena圖像等圖像處理相關領域[11]。
本文首創性地將2D-EEMD算法應用到腦電處理領域，對原始EEG信號進行2D-EEMD分解，觀測SSVEP視覺刺激頻率及其諧波頻率，以發現有效的主頻率成分。再將結果投射到頭圖上，觀察規律性明暗變化的視覺刺激頻率，以及對視覺刺激反應最強烈的區域。最后將有效成分進行重構信號，并用短時傅立葉變換(short time Fourier transform，STFT)實現SSVEP頻率提取。
1 研究方法
1.1 數據采集
數據由美國加州大學圣迭戈分校的Swartz計算神經研究中心(Swartz Center for Computational Neuroscience，SCCN)的實驗室提供，實驗設備采用Bio-semi系統。5個視力正常的被試者(subjects)舒服地坐在電腦屏幕前方的椅子上，佩戴264個電極(channel)的電極帽進行EEG信號的采集，采樣頻率為2 048 Hz。被試者視線與電腦屏幕保持水平。電腦屏幕中央有一個3 cm×3 cm大小、以固定頻率進行黑白閃爍的方塊，即SSVEP視覺刺激頻率[12-13]。每人每次實驗包括4段(session)，每段包括10個30 s的試驗(trial)，SSVEP刺激頻率從9~13 Hz中隨機選取。本文選取了9 Hz和10 Hz兩個頻率。每30 s的刺激之后休息15 s以消除被試者的視覺疲勞，每段試驗間也有幾分鐘的休息間隔[12-13]。
采集到的SSVEP數據中去除了13個電極粘連不牢固而采集到的失效信號，然后進行了256 Hz的重新采樣。圖 1顯示了9 Hz刺激下，從1號被試者(Subject 01)采集的數據前1 s的二維圖像，橫坐標為時間軸，縱坐標為有效的251個電極。全文統一采用圖 1所示標尺，灰度值在[-10, 10]之間，灰度值代表信號電位的大小，顏色越深電位越小，顏色越淺電位越大。由于各電極之間電位漂移值差距很大，因此在原始的二維圖像上很難觀察到信號隨時間的變化情況。
圖1
1號被試者的1 s采樣數據的二維圖像
Figure1.
2D image of 1 s sample data (Subject 01)
1.2 2D-EEMD
EEMD是為了解決EMD算法存在的模態混疊問題而提出來的[14-15]，每次分解都給原始信號附加不同振幅的高斯白噪聲，只要附加的高斯白噪聲的振幅是有限的，那么在保證總實驗次數足夠大的情況下，利用白噪聲均值為零的特性，對每一個IMF求取集合平均，不僅可消除附加白噪聲所帶來的影響，也解決了模態混疊問題。以EEMD為基礎的2D-EEMD，則是在二維數據或圖像的任意兩個正交方向上分別進行EEMD分解[11]。
對于任意一個二維數據f(m, n)，包含M個空間點，N個時間點，
| $f\left( m,n \right)\left( \begin{matrix} {{f}_{1,1}} & {{f}_{2,1}} & \cdots & {{f}_{M,1}} \\ {{f}_{1,2}} & {{f}_{2,2}} & \cdots & {{f}_{M,2}} \\ \vdots & \vdots & \cdots & \vdots \\ {{f}_{1,N}} & {{f}_{2,N}} & \cdots & {{f}_{M,N}} \\ \end{matrix} \right) $ |
每列代表一個時間序列，對第m列
| $f\left( m,\tilde{\ } \right)=\left( \begin{align} & {{f}_{m,1}} \\ & {{f}_{m,2}} \\ & \vdots \\ & {{f}_{m,N}} \\ \end{align} \right) $ |
進行一維EEMD分解，得到
| $f\left( m,\tilde{\ } \right)=\sum\limits_{j=i}^{\text{J}}{{{c}_{j}}\left( m,\tilde{\ } \right)=}\sum\limits_{j=i}^{\text{J}}{\left( \begin{align} & {{c}_{m,1,j}} \\ & {{c}_{m,2,j}} \\ & \vdots \\ & {{c}_{m,N,j}} \\ \end{align} \right)} $ |
其中j表示第j個IMF分量。對f(m, n)的所有m列進行一維EEMD分解，得到J個矩陣，其中第j個矩陣為：
| $ {{g}_{j}}\left( m,n \right)=\left( \begin{matrix} {{c}_{1,1,j}} & {{c}_{2,1,j}} & \cdots & {{c}_{M,1,j}} \\ {{c}_{1,2,j}} & {{c}_{2,2,j}} & \cdots & {{c}_{M,1,j}} \\ \vdots & \vdots & \cdots & \vdots \\ {{c}_{1,N,j}} & {{c}_{2,N,j}} & \cdots & {{c}_{M,N,j}} \\ \end{matrix} \right) $ |
上述J個矩陣即f(m, n)的一維分解結果，接下來進行二維分解：對于gj(m, n)的每一行數據進行EEMD分解。gj(m, n)第n行數據及其EEMD分解結果分別為
| $ {{g}_{j}}\left(\tilde{\ }, n \right)=\left({{c}_{1, n, j\ \ \ }}{{c}_{2, n, j}}\ \cdots \ \ {{c}_{M, n, j}} \right) $ |
和
| $ \begin{align} & {{g}_{j}}\left( \tilde{\ },n \right)=\sum\limits_{k=1}^{\text{K}}{{{D}_{j,k}}}\left( \tilde{\ },n \right) \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum\limits_{k=1}^{\text{K}}{\left( {{d}_{1,n,j,k}}\ \ {{d}_{2,n,j,k}}\ \ \ \cdots \ \ {{d}_{M,n,j,k}} \right)} \\ \end{align} $ |
對gj(m, n)所有N行數據分別進行EEMD分解，得到K個矩陣，其中第k個矩陣為:
| $ {{h}_{j,k}}\left( m,n \right)=\left( \begin{matrix} {{d}_{1,1,j,k}} & {{d}_{2,1,j,k}} & \cdots & {{d}_{M,1,j,k}} \\ {{d}_{1,2,j,k}} & {{d}_{2,1,j,k}} & \cdots & {{d}_{M,2,j,k}} \\ \vdots & \vdots & \cdots & \vdots \\ {{d}_{1,N,j,k}} & {{d}_{2,N,j,k}} & \cdots & {{d}_{M,N,j,k}} \\ \end{matrix} \right) $ |
對所有J個一維分解結果進行二維分解則得到J×K個矩陣，因此信號最終可以表示為
| $ f\left(m, n \right)=\sum\limits_{i=1}^{J}{\sum\limits_{k=1}^{K}{{{h}_{j, k}}}}\left(m, n \right) $ |
基于上述原理，我們對251(電極)·4(秒)數據引入了2D-EEMD處理。對9 Hz視覺刺激下的SSVEP數據進行2D-EEMD的結果如圖 2所示(J=K=5)，從中尚難以確定SSVEP的有效成分。為了解決該問題，引入文獻[11]提出的“可比較的小尺度組合準則”的重組方法：當在兩個正交方向上應用EEMD時，如果兩個IMF分量間的尺度是可以比較的，那么它們的重組成分會有最佳的物理含義。根據研究，圖 2中的二維分解圖像中，每一列圖像之間都有近似相同的垂直尺度，每一行圖像之間都有近似相同的水平尺度，所以最終重組了K個IMF成分，其中第i個IMF成分Ci可以表示為：
| $ {{c}_{i}}=\sum\limits_{k=1}^{K}{{{h}_{i, k}}+}\sum\limits_{j=i+1}^{K}{{{h}_{j, i}}} $ |
圖2
9 Hz視覺刺激下SSVEP的2D-EEMD分解結果(J=K=5)
對于每一個子圖，橫坐標表示1 s的數據采集時間，縱坐標表示251個電極
Figure2. 2D-EEMD results of 9 Hz-stimulated SSVEP (J=K=5)for each graph,
該方法應用于EEG信號時，存在一個問題：圖 2的最后一列分解成分屬于低頻漂移干擾成分，其幅值變化巨大。如果將其加入重組，將會覆蓋掉SSVEP的特征。因此，本文改進原有重組規則，去除最后一列成分(由于右下角的成分為重組的IMF5唯一組成，所以仍舊保留)，最終重組得到的5個IMF結果如圖 3所示。
圖3
9 Hz視覺刺激下2D-EEMD重組后的IMF1~IMF5和對應的頻譜圖
Figure3.
2D-EEMD results and their power spectra of recombined IMFs (IMF1 to IMF5) at 9 Hz stimulated SSVEP
1.3 STFT
STFT作為傅里葉變換的一種改進算法，通過在原始信號上附加一個隨時間不斷滑動的分析窗，對信號進行加窗截斷，將非平穩信號分解成一系列近似平穩的短時信號，對這些短時信號再應用傅里葉變換可以得到其近似的頻譜分析[16-17]。所以STFT非常適合處理像EEG信號的非平穩信號。
2 結果
2.1 IMFs
圖 3為2D-EEMD分解得到的IMFs二維圖像。從IMF對應的頻譜圖，可以知道IMF1高頻成分較多，而EEG信號主要集中在1~50 Hz范圍，因此IMF1為各種噪聲成分，包括工頻噪聲。IMF2的頻率峰值為18、27、36、45 Hz, 均為SSVEP頻率9 Hz的諧波頻率。IMF3和IMF4中則以9 Hz為主，因此為SSVEP主要成分；且從IMF3和IMF4的圖像中，我們可以看到在1 s時間內，峰值(較淺色)或峰谷(較深色)的個數正好為9個，這也對應了9 Hz的響應頻率，這種空間上的一致性更加突顯了2D-EEMD分解的直觀性。IMF5為剩余信號，頻率集中在1 Hz以內，屬于低頻成分，基于目前的研究基礎認為該分量亦為低頻噪聲。
對10 Hz視覺刺激下的SSVEP二維信號，我們也進行了分解。圖 4顯示了IMF3和IMF4的二維圖像和頻譜圖，同樣證明了這兩個成分為SSVEP的主要響應成分。
圖4
10 Hz視覺刺激下2D-EEMD重組后的IMF3、IMF4和對應的頻譜圖
Figure4.
2D-EEMD results and their power spectra of recombined IMFs (IMF3, IMF4) at 10 Hz stimulated SSVEP
2.2 頭圖
將分解的IMF主要成分在不同時間點上的數據投影到頭圖上，可以觀察該成分隨時間變化的趨勢。圖 5顯示了9 Hz刺激下IMF3投影到頭圖的時變情況。整個時長為1 s，從0 s開始，以1/18 s為步長。我們可以看到9 Hz響應頻率的規律性明暗變化，說明了大腦各區域之間具有高度的響應同步性。同時，對照黑白強度可以發現后腦枕葉區的振幅最強，說明該區域對于視覺刺激的反應最強烈。
圖5
9 Hz視覺刺激下，1 s內的頭圖變化趨勢
Figure5.
Topoplot graphs' changing trend at 9 Hz stimulation SSVEP in 1 second
這種二維分解結果投影到頭圖的意義在于EEG的電極之間位置接近，同一時刻的相鄰電極數據具有相關性，因此這樣的二維分解是有物理意義的，這種分析方法為分析EEG信號提供了一個有效的手段。
2.3 枕葉區與非枕葉區對視覺刺激的反應強度
圖 6比較了各IMF成分在枕葉區和非枕葉區提取的SSVEP頻率及其諧波頻率的幅值。可以發現在枕葉區提取到的SSVEP頻率及其諧波頻率在IMF1到IMF3中都遠大于其他區的頻率，在IMF4中兩個區域的值基本相等，而且可以非常清晰地看到IMF3是最有效的IMF。圖 6的結果也說明了枕葉區在人腦的視覺刺激反應中占據主導地位。
圖6
枕葉區和非枕葉區的反應強度比較
Figure6.
Comparison of reacting intensity between occipital region and other regions
2.4 SSVEP頻率提取準確率
實驗中使用STFT來處理EEG信號。圖 7顯示了5個被試者平均的檢測準確度隨STFT的滑動窗長度變化的曲線。通過比較我們可以發現，隨著窗口長度的增加，準確度增加，而且在同一個窗口長度下，2D-EEMD分解后的IMF的識別準確度也有明顯的增加。例如在窗口長度為0.2 s時，識別準確度從55%(原始數據)增加到71%(IMF1, 2, 3, 4)。
圖7
STFT頻率提取準確度比較
Figure7.
Comparison of frequency accuracy rate by STFT
3 結論
本文主要研究了2D-EEMD對于SSVEP頻率的檢測。2D-EEMD在數據的兩個正交方向上獨立地進行EEMD分解，在“可比較的小尺度組合準則”這一重組準則下，對于二維分解結果進行了有效重組，得到了最終的IMFs成分。從最后的有效分解成分IMF3和IMF4，以及頭圖的規律明暗變化中，我們可以清晰地觀測到SSVEP的響應頻率。將IMF主要成分投影到頭圖上，可以觀察到各個區域對視覺刺激的響應情況和時變情況。最后使用STFT提取反應頻率，可以發現2D-EEMD處理使頻率的識別準確度得到了提升。
引言
腦機接口(brain computer interface, BCI)是通過腦電圖(electroencephalogram，EEG)信號的轉化而建立起來的大腦與外設之間的控制系統，它不依賴于常規的腦信號輸出通路——外周神經及肌肉組織，可以為癱瘓或者神經肌肉受損的殘疾人提供全新的通信控制方式[1-2]。基于EEG的BCI系統憑借其非侵入式檢測所帶來的低傷害性和便捷性，在當今的BCI系統中得到了廣泛的應用。其中，穩態視覺誘發電位(steady state visual evoked potential，SSVEP)在近幾年的研究中更是引起了廣泛的關注。SSVEP是由持續視覺刺激而誘發的節律性EEG信號[3]，SSVEP頻率由固定的視覺刺激頻率及其諧波頻率組成[4]，因此通過研究SSVEP可以準確提取出視覺刺激頻率。將不同的任務和不同的刺激頻率相組合，可以真正實現人機交互。
經驗模式分解(empirical mode decomposition，EMD)是希爾伯特-黃變換(Hilbert-Huang transform，HHT)的預處理過程[5]，EMD分解以信號自身為基礎，通過反復篩選得到的本征模態函數(intrinsic mode function，IMF)是自適應的，適合處理像腦電一類的非平穩非線性信號[6]。集合經驗模式分解(ensemble EMD, EEMD)作為一種新的噪聲輔助分析法，在一定程度上緩解了EMD所帶來的模態混疊問題[7-10]。而二維集合經驗模式分解(two-dimensional EEMD, 2D-EEMD)則在兩個正交方向上分別獨立地對信號進行EEMD分解，然后根據“可比較的小尺度組合準則”對分解結果進行重新組合，得到最終的二維分解結果[11]。2D-EEMD已廣泛應用在海洋颶風、核磁共振成像(nuclear magnetic resonance imaging，MRI)、Lena圖像等圖像處理相關領域[11]。
本文首創性地將2D-EEMD算法應用到腦電處理領域，對原始EEG信號進行2D-EEMD分解，觀測SSVEP視覺刺激頻率及其諧波頻率，以發現有效的主頻率成分。再將結果投射到頭圖上，觀察規律性明暗變化的視覺刺激頻率，以及對視覺刺激反應最強烈的區域。最后將有效成分進行重構信號，并用短時傅立葉變換(short time Fourier transform，STFT)實現SSVEP頻率提取。
1 研究方法
1.1 數據采集
數據由美國加州大學圣迭戈分校的Swartz計算神經研究中心(Swartz Center for Computational Neuroscience，SCCN)的實驗室提供，實驗設備采用Bio-semi系統。5個視力正常的被試者(subjects)舒服地坐在電腦屏幕前方的椅子上，佩戴264個電極(channel)的電極帽進行EEG信號的采集，采樣頻率為2 048 Hz。被試者視線與電腦屏幕保持水平。電腦屏幕中央有一個3 cm×3 cm大小、以固定頻率進行黑白閃爍的方塊，即SSVEP視覺刺激頻率[12-13]。每人每次實驗包括4段(session)，每段包括10個30 s的試驗(trial)，SSVEP刺激頻率從9~13 Hz中隨機選取。本文選取了9 Hz和10 Hz兩個頻率。每30 s的刺激之后休息15 s以消除被試者的視覺疲勞，每段試驗間也有幾分鐘的休息間隔[12-13]。
采集到的SSVEP數據中去除了13個電極粘連不牢固而采集到的失效信號，然后進行了256 Hz的重新采樣。圖 1顯示了9 Hz刺激下，從1號被試者(Subject 01)采集的數據前1 s的二維圖像，橫坐標為時間軸，縱坐標為有效的251個電極。全文統一采用圖 1所示標尺，灰度值在[-10, 10]之間，灰度值代表信號電位的大小，顏色越深電位越小，顏色越淺電位越大。由于各電極之間電位漂移值差距很大，因此在原始的二維圖像上很難觀察到信號隨時間的變化情況。
圖1
1號被試者的1 s采樣數據的二維圖像
Figure1.
2D image of 1 s sample data (Subject 01)
1.2 2D-EEMD
EEMD是為了解決EMD算法存在的模態混疊問題而提出來的[14-15]，每次分解都給原始信號附加不同振幅的高斯白噪聲，只要附加的高斯白噪聲的振幅是有限的，那么在保證總實驗次數足夠大的情況下，利用白噪聲均值為零的特性，對每一個IMF求取集合平均，不僅可消除附加白噪聲所帶來的影響，也解決了模態混疊問題。以EEMD為基礎的2D-EEMD，則是在二維數據或圖像的任意兩個正交方向上分別進行EEMD分解[11]。
對于任意一個二維數據f(m, n)，包含M個空間點，N個時間點，
| $f\left( m,n \right)\left( \begin{matrix} {{f}_{1,1}} & {{f}_{2,1}} & \cdots & {{f}_{M,1}} \\ {{f}_{1,2}} & {{f}_{2,2}} & \cdots & {{f}_{M,2}} \\ \vdots & \vdots & \cdots & \vdots \\ {{f}_{1,N}} & {{f}_{2,N}} & \cdots & {{f}_{M,N}} \\ \end{matrix} \right) $ |
每列代表一個時間序列，對第m列
| $f\left( m,\tilde{\ } \right)=\left( \begin{align} & {{f}_{m,1}} \\ & {{f}_{m,2}} \\ & \vdots \\ & {{f}_{m,N}} \\ \end{align} \right) $ |
進行一維EEMD分解，得到
| $f\left( m,\tilde{\ } \right)=\sum\limits_{j=i}^{\text{J}}{{{c}_{j}}\left( m,\tilde{\ } \right)=}\sum\limits_{j=i}^{\text{J}}{\left( \begin{align} & {{c}_{m,1,j}} \\ & {{c}_{m,2,j}} \\ & \vdots \\ & {{c}_{m,N,j}} \\ \end{align} \right)} $ |
其中j表示第j個IMF分量。對f(m, n)的所有m列進行一維EEMD分解，得到J個矩陣，其中第j個矩陣為：
| $ {{g}_{j}}\left( m,n \right)=\left( \begin{matrix} {{c}_{1,1,j}} & {{c}_{2,1,j}} & \cdots & {{c}_{M,1,j}} \\ {{c}_{1,2,j}} & {{c}_{2,2,j}} & \cdots & {{c}_{M,1,j}} \\ \vdots & \vdots & \cdots & \vdots \\ {{c}_{1,N,j}} & {{c}_{2,N,j}} & \cdots & {{c}_{M,N,j}} \\ \end{matrix} \right) $ |
上述J個矩陣即f(m, n)的一維分解結果，接下來進行二維分解：對于gj(m, n)的每一行數據進行EEMD分解。gj(m, n)第n行數據及其EEMD分解結果分別為
| $ {{g}_{j}}\left(\tilde{\ }, n \right)=\left({{c}_{1, n, j\ \ \ }}{{c}_{2, n, j}}\ \cdots \ \ {{c}_{M, n, j}} \right) $ |
和
| $ \begin{align} & {{g}_{j}}\left( \tilde{\ },n \right)=\sum\limits_{k=1}^{\text{K}}{{{D}_{j,k}}}\left( \tilde{\ },n \right) \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum\limits_{k=1}^{\text{K}}{\left( {{d}_{1,n,j,k}}\ \ {{d}_{2,n,j,k}}\ \ \ \cdots \ \ {{d}_{M,n,j,k}} \right)} \\ \end{align} $ |
對gj(m, n)所有N行數據分別進行EEMD分解，得到K個矩陣，其中第k個矩陣為:
| $ {{h}_{j,k}}\left( m,n \right)=\left( \begin{matrix} {{d}_{1,1,j,k}} & {{d}_{2,1,j,k}} & \cdots & {{d}_{M,1,j,k}} \\ {{d}_{1,2,j,k}} & {{d}_{2,1,j,k}} & \cdots & {{d}_{M,2,j,k}} \\ \vdots & \vdots & \cdots & \vdots \\ {{d}_{1,N,j,k}} & {{d}_{2,N,j,k}} & \cdots & {{d}_{M,N,j,k}} \\ \end{matrix} \right) $ |
對所有J個一維分解結果進行二維分解則得到J×K個矩陣，因此信號最終可以表示為
| $ f\left(m, n \right)=\sum\limits_{i=1}^{J}{\sum\limits_{k=1}^{K}{{{h}_{j, k}}}}\left(m, n \right) $ |
基于上述原理，我們對251(電極)·4(秒)數據引入了2D-EEMD處理。對9 Hz視覺刺激下的SSVEP數據進行2D-EEMD的結果如圖 2所示(J=K=5)，從中尚難以確定SSVEP的有效成分。為了解決該問題，引入文獻[11]提出的“可比較的小尺度組合準則”的重組方法：當在兩個正交方向上應用EEMD時，如果兩個IMF分量間的尺度是可以比較的，那么它們的重組成分會有最佳的物理含義。根據研究，圖 2中的二維分解圖像中，每一列圖像之間都有近似相同的垂直尺度，每一行圖像之間都有近似相同的水平尺度，所以最終重組了K個IMF成分，其中第i個IMF成分Ci可以表示為：
| $ {{c}_{i}}=\sum\limits_{k=1}^{K}{{{h}_{i, k}}+}\sum\limits_{j=i+1}^{K}{{{h}_{j, i}}} $ |
圖2
9 Hz視覺刺激下SSVEP的2D-EEMD分解結果(J=K=5)
對于每一個子圖，橫坐標表示1 s的數據采集時間，縱坐標表示251個電極
Figure2. 2D-EEMD results of 9 Hz-stimulated SSVEP (J=K=5)for each graph,
該方法應用于EEG信號時，存在一個問題：圖 2的最后一列分解成分屬于低頻漂移干擾成分，其幅值變化巨大。如果將其加入重組，將會覆蓋掉SSVEP的特征。因此，本文改進原有重組規則，去除最后一列成分(由于右下角的成分為重組的IMF5唯一組成，所以仍舊保留)，最終重組得到的5個IMF結果如圖 3所示。
圖3
9 Hz視覺刺激下2D-EEMD重組后的IMF1~IMF5和對應的頻譜圖
Figure3.
2D-EEMD results and their power spectra of recombined IMFs (IMF1 to IMF5) at 9 Hz stimulated SSVEP
1.3 STFT
STFT作為傅里葉變換的一種改進算法，通過在原始信號上附加一個隨時間不斷滑動的分析窗，對信號進行加窗截斷，將非平穩信號分解成一系列近似平穩的短時信號，對這些短時信號再應用傅里葉變換可以得到其近似的頻譜分析[16-17]。所以STFT非常適合處理像EEG信號的非平穩信號。
2 結果
2.1 IMFs
圖 3為2D-EEMD分解得到的IMFs二維圖像。從IMF對應的頻譜圖，可以知道IMF1高頻成分較多，而EEG信號主要集中在1~50 Hz范圍，因此IMF1為各種噪聲成分，包括工頻噪聲。IMF2的頻率峰值為18、27、36、45 Hz, 均為SSVEP頻率9 Hz的諧波頻率。IMF3和IMF4中則以9 Hz為主，因此為SSVEP主要成分；且從IMF3和IMF4的圖像中，我們可以看到在1 s時間內，峰值(較淺色)或峰谷(較深色)的個數正好為9個，這也對應了9 Hz的響應頻率，這種空間上的一致性更加突顯了2D-EEMD分解的直觀性。IMF5為剩余信號，頻率集中在1 Hz以內，屬于低頻成分，基于目前的研究基礎認為該分量亦為低頻噪聲。
對10 Hz視覺刺激下的SSVEP二維信號，我們也進行了分解。圖 4顯示了IMF3和IMF4的二維圖像和頻譜圖，同樣證明了這兩個成分為SSVEP的主要響應成分。
圖4
10 Hz視覺刺激下2D-EEMD重組后的IMF3、IMF4和對應的頻譜圖
Figure4.
2D-EEMD results and their power spectra of recombined IMFs (IMF3, IMF4) at 10 Hz stimulated SSVEP
2.2 頭圖
將分解的IMF主要成分在不同時間點上的數據投影到頭圖上，可以觀察該成分隨時間變化的趨勢。圖 5顯示了9 Hz刺激下IMF3投影到頭圖的時變情況。整個時長為1 s，從0 s開始，以1/18 s為步長。我們可以看到9 Hz響應頻率的規律性明暗變化，說明了大腦各區域之間具有高度的響應同步性。同時，對照黑白強度可以發現后腦枕葉區的振幅最強，說明該區域對于視覺刺激的反應最強烈。
圖5
9 Hz視覺刺激下，1 s內的頭圖變化趨勢
Figure5.
Topoplot graphs' changing trend at 9 Hz stimulation SSVEP in 1 second
這種二維分解結果投影到頭圖的意義在于EEG的電極之間位置接近，同一時刻的相鄰電極數據具有相關性，因此這樣的二維分解是有物理意義的，這種分析方法為分析EEG信號提供了一個有效的手段。
2.3 枕葉區與非枕葉區對視覺刺激的反應強度
圖 6比較了各IMF成分在枕葉區和非枕葉區提取的SSVEP頻率及其諧波頻率的幅值。可以發現在枕葉區提取到的SSVEP頻率及其諧波頻率在IMF1到IMF3中都遠大于其他區的頻率，在IMF4中兩個區域的值基本相等，而且可以非常清晰地看到IMF3是最有效的IMF。圖 6的結果也說明了枕葉區在人腦的視覺刺激反應中占據主導地位。
圖6
枕葉區和非枕葉區的反應強度比較
Figure6.
Comparison of reacting intensity between occipital region and other regions
2.4 SSVEP頻率提取準確率
實驗中使用STFT來處理EEG信號。圖 7顯示了5個被試者平均的檢測準確度隨STFT的滑動窗長度變化的曲線。通過比較我們可以發現，隨著窗口長度的增加，準確度增加，而且在同一個窗口長度下，2D-EEMD分解后的IMF的識別準確度也有明顯的增加。例如在窗口長度為0.2 s時，識別準確度從55%(原始數據)增加到71%(IMF1, 2, 3, 4)。
圖7
STFT頻率提取準確度比較
Figure7.
Comparison of frequency accuracy rate by STFT
3 結論
本文主要研究了2D-EEMD對于SSVEP頻率的檢測。2D-EEMD在數據的兩個正交方向上獨立地進行EEMD分解，在“可比較的小尺度組合準則”這一重組準則下，對于二維分解結果進行了有效重組，得到了最終的IMFs成分。從最后的有效分解成分IMF3和IMF4，以及頭圖的規律明暗變化中，我們可以清晰地觀測到SSVEP的響應頻率。將IMF主要成分投影到頭圖上，可以觀察到各個區域對視覺刺激的響應情況和時變情況。最后使用STFT提取反應頻率，可以發現2D-EEMD處理使頻率的識別準確度得到了提升。
