工作記憶是高級認知功能的重要基礎。本文結合腦電圖（EEG）和功能近紅外成像（fNIRS）的時空優勢研究工作記憶的神經血管耦合機制。在數據分析中，提取EEG數據中不同試次的時間序列與血氧動力學響應函數（HRF）卷積后的矩陣和fNIRS的血氧變化矩陣作為耦合的特征。然后，使用典型相關分析（CCA）計算兩種模態特征間的交叉關聯。結果表明，CCA算法能夠提取出相關成分試次間相近的變化趨勢，并發現額極區和背外側前額葉的fNIRS激活與EEG數據的delta、theta和alpha節律相關。本研究揭示了工作記憶下的神經血管耦合機制，為EEG數據和fNIRS數據融合提供了新方法。
引用本文： 劉文正, 張昊, 楊柳, 古悅. 基于腦電圖與功能近紅外成像的工作記憶神經血管耦合分析. 生物醫學工程學雜志, 2022, 39(2): 228-236, 247. doi: 10.7507/1001-5515.202108048 復制
引言
近年來，認知任務的研究越來越多，高級認知功能已經成為認知神經學、教育學和臨床醫學等諸多學科的重要話題。工作記憶是指人的大腦在執行認知任務過程中，用于大腦的信息加工和短時存儲過程的一個容量有限的系統[1]。工作記憶是高級認知功能的重要基礎，對于認知科學的研究具有重要意義[2]。隨著腦功能成像技術的不斷發展，在工作記憶時人類大腦的變化可以通過腦成像設備檢測，例如腦電圖（electroencephalography，EEG）和功能近紅外成像（functional near-infrared spectroscopy，fNIRS）。EEG數據是通過電極將神經元自發、節律性的運動電位記錄下來而得到的數據，具有很高的時間分辨率[3]。但是，EEG設備無法深入到大腦皮層以下采集信號，空間分辨率較低。有研究使用EEG設備進行溯源分析，發現音樂會降低記憶學習的效率[4]。在工作記憶任務中，學習障礙兒童的EEG數據具有更多的delta和theta節律，以及gamma節律的減少[5]。fNIRS技術利用血液對近紅外光的散射性測量出氧合血紅蛋白（oxygenated hemoglobins，HbO）和脫氧血紅蛋白（deoxygenated hemoglobins，HbR）濃度，具有價格低、無創性和易攜帶等特點，受到腦科學研究的高度重視[6]。基于fNIRS技術，有研究表明適度的運動有利于老年人進行工作記憶任務[7]。在工作記憶任務時，使用fNIRS技術發現左撇子和右撇子有不同的大腦皮層間功能連接[8]。
由于不同的采集設備具有不同的特性，因此單一采集設備的研究會具有一定局限性。基于神經血管耦合機制，融合多種設備來達到一種設備無法比擬的時空分辨率，可以從不同角度研究大腦功能[9]。在國內外神經血管耦合的研究中，EEG數據和功能磁共振成像（functional magnetic resonance imaging，fMRI）數據融合的居多。有研究使用EEG數據與fMRI數據探究睡眠和覺醒狀態下的神經血管耦合穩定性[10]。事件相關電位（event related potentials，ERP）間的振幅調制（amplitude modulations，AMs）與fMRI大腦激活也具有關聯性[11]。雖然，有很多學者研究EEG數據與fMRI數據融合，但是在融合EEG數據與fMRI數據時存在很多問題。例如，fMRI采集過程中要求受試者不能動，否則會對采集的數據造成偽影，影響分析的結果；fMRI釋放的磁場會影響EEG數據與fMRI數據聯合采集的效果。所以，如何解決兩者融合過程中出現的問題也成為了相關研究的熱點。fNIRS是一種新興的大腦信號測量設備，它和fMRI都是基于腦功能血氧動力學的測量方法。相比于fMRI，fNIRS可以在運動的環境下使用，且與EEG設備沒有不良干擾。因此，EEG數據和fNIRS數據的融合具有很大的研究意義和應用前景。在無創腦機接口（brain computer interface，BCI）研究領域，有研究證實EEG數據和fNIRS數據的多模態融合提升了BCI的性能[12]。
在已有的融合方法中，基于數據驅動的方法不需要對復雜的神經元種群和神經血管耦合動力學建模[13]。基于數據驅動的方法有偏最小二乘法（partial least squares，PLS）和典型相關分析（canonical correlation analysis，CCA）等一系列經典算法。因為CCA能夠更好地應用到多模態、多受試者的數據[14]，所以本文使用CCA計算EEG數據試次間的AMs和fNIRS數據的最大相關性，并提取出有關聯的EEG數據時間分量和fNIRS數據空間分量，有效促進了兩者的融合。綜上所述，使用CCA算法融合EEG數據和fNIRS數據，可以從時間和空間兩個方面研究大腦的神經活動。
本文采用工作記憶中經典范式n-back，通過采集設備獲得同步EEG數據和fNIRS數據。首先對采集到的兩種模態的數據進行預處理和特征提取。對特征提取后的矩陣，使用CCA計算模態間相關系數以及對應的EEG數據時間分量和fNIRS數據空間分量。通過對EEG數據時間分量和fNIRS數據空間分量的分析，得到了與工作記憶相關的結果[15]。
1 EEG數據和fNIRS數據融合理論
EEG數據與fNIRS數據融合模型如圖1所示。EEG數據矩陣經過預處理后的行是單個試次的采樣點，列是試次間的變化。然后，將預處理后的矩陣與血氧動力學響應函數（hemodynamic response function，HRF）卷積得到特征矩陣X。HRF是描述fNIRS對神經元瞬間刺激的響應函數。當大腦的某些區域進行神經活動時，氧氣消耗量增加并且刺激血管擴張以增加血流量來提供更多的HbO，達到給神經元活動提供氧氣的目的。在神經元刺激后，HbO濃度會增加，HbR濃度會下降，這種現象形成了血氧動力學響應特征Ri(t)，其計算公式如式（1）所示[16]：
圖1
EEG數據與fNIRS數據融合模型
Figure1.
The fusion model of EEG data and fNIRS data
 |
式中，自變量t是時間，exp是以e為底的指數函數，i = 1代表HbO模型，i = 2代表HbR模型。對于HbO模型，C = 8log2、τ1,1 = 7.8、w1,1 = 5.95、τ1,2 = 12.16、w1,2 = 16.35和r1=0.37；對于HbR模型，C = 8log2、τ2,1 = 8.5、w2,1 = 7.3、τ2,2 = 21.4、w2,2 = 20.31和r2 = 0.331。
fNIRS數據矩陣經過預處理和特征提取后，得到特征矩陣Y。特征矩陣Y的行是所有通道，列是在時間上的采樣。特征矩陣X的試次間變化和Y的時間上的采樣相互對應，都是隨時間變化的數據采樣。特征矩陣X和Y的融合模型如式（2）所示：
 |
式中，
，
，
，
，
，
，其中R代表實數集。TX、TY是矩陣X和Y的時間維度（數據的采樣點），VX、VY是矩陣X和Y的特征變量，DX、DY是矩陣X和Y的最小秩。相關矩陣AX和AY中第i列向量為AiX和AiY（i=1，2，
，D），且同一相關矩陣AX或AY的列之間不相關；不同相關矩陣AX和AY的列之間相關，且相乘等于相關系數。相關矩陣AX和AY可由CCA計算得到。
CCA算法整體思想：對于X矩陣，存在投影向量aT將X投影到一維；對于Y矩陣，存在投影向量bT將Y投影到一維；投影后的一維向量分別是X′和Y′，那么可以表示成X′ = aTX和Y′ = bTY，并計算矩陣X′和Y′間的最大相關性[17]。
如式（3）所示，CCA的優化目標是最大化相關系數λ，并算出對應的投影向量a和b。
 '/> |
其中，SXY是矩陣X和Y的互協方差矩陣，SXX是矩陣X的協方差矩陣，SYY是矩陣Y的協方差矩陣。公式（3）的限定條件是aTSXXa = 1，bTSYYb = 1。然后，使用拉格朗日特征分解方法優化公式（3）得到的特征方程組如式（4）所示：
 |
解特征方程組（4）得到最大相關系數λ1和第1對投影向量a1和b1。
當得到前k ? 1對投影向量ai和bi（i = 1，
，k ? 1）后，第k對投影向量的限定條件，就如式（5）所示：
 |
然后，再次用拉格朗日方法優化公式（3）得到相關系數λk和第k對投影向量ak和bk。aiX和biY是X和Y的第i對典型變量Ai。如式（6）所示，不同相關矩陣Ai之間相關，且相乘等于相關系數λi。
 |
如式（7）所示，因為模型與最小二乘法原理相同，所以可以用最小二乘法擬合出模態的相關成分Ci。
 |
2 材料和方法
2.1 受試者
本研究共有23名成年人參與數據采集，全部為男性（平均年齡22歲，年齡范圍18～26歲）。所有受試者都是右利手，且都身體健康、沒有心理或精神類疾病。每個受試者都在任務前簽署了受試者知情同意書。本研究經北京師范大學研究倫理審查委員會根據赫爾辛基宣言批準。
2.2 n-back范式
本研究采用數字設計的n-back范式，包括3種負荷任務：0-back、1-back和2-back。不同負荷任務的目標刺激不同。0-back負荷時，目標數字刺激是“7”，其他數字刺激是非目標刺激。1-back負荷時，如果當前數字和前1個數字相同則為目標刺激，不同則為非目標刺激。2-back負荷時，如果當前數字和前2個數字相同則為目標刺激，不同則為非目標刺激。以上所述“目標刺激”需要受試者按下鍵盤空格鍵。本范式中每種負荷任務包含10個組塊（block），3種不同的負荷任務按隨機順序顯示。n-back范式任務流程如圖2所示，以1-back為例：首先，屏幕上顯示5 s的任務類型“1-back”、1 s的注視窗口“+”；然后，每2.2 s隨機出現1個數字，每個數字出現0.5 s，接著是1.7 s的注視窗口“+”，每個block共有15個數字；最后，每個block的任務結束后有20 s的休息時間。在每個block中，“目標刺激”出現的概率為三分之一，“非目標刺激”概率為三分之二。因此，本范式共有3×10×15=450個試次。
圖2
n-back范式
Figure2.
n-back paradigm
2.3 EEG數據獲取
使用多通道EEG設備（ActiCHamp，Brain Products.，德國）以1 000 Hz的采樣率記錄連續的任務數據。EEG數據采集電極共32個（AFp1、AFp2、F7、AFF1h、AFF2h、F8、FC5、FC1、FC2、FC6、T7、C3、C4、T8、TP9、CP5、CP1、CP2、CP6、TP10、P7、P3、Pz、P4、P8、AFF5h、O1、Oz、O2、AFF6h、Cz和FCz），其中，Cz是參考電極，FCz是接地電極，所有電極的阻抗保持在5 kΩ以下。
2.4 fNIRS數據獲取
使用fNIRS設備（NIRSport，NIRx Medical Technologies.，美國）以7.812 5 Hz的采樣率記錄連續的任務數據。相鄰的光源和探測器配置為fNIRS通道，使用8個光源和7個探測器在EEG電極（AF7、Fp1、Fpz、Fp2、AF8、AF3、AFz、AF4、F5、F3、F1、Fz、F2、F4、F6）周圍形成了20個通道。20個通道的光源和探測器距離均設置為30 mm[18]。fNIRS設備的光源、探測器和EEG設備的電極被固定在同一個電極帽上。
2.5 EEG數據預處理和特征提取
使用EEG數據處理軟件Analyzer 2.2（Brain Products.，德國）對采集到的EEG數據進行預處理。首先，對EEG數據進行200 Hz降采樣，并對所有EEG電極平均重參考。然后，對EEG數據進行1～40 Hz濾波、50 Hz陷波、偽跡檢測、壞通道替換。本文使用獨立分量分析（independent component analysis，ICA）識別并去除眼電噪聲。最后，將數據按照標記基線校正，取0～800 ms的數據做分析，并對每個試次進行小波去噪[19]。
預處理后的數據標準化為零平均值和單位方差。選取AMs現象明顯的額中央電極（Cz、FC1和FC2）平均值作為EEG特征。接下來分別對3種負荷任務下的EEG數據提取特征。首先，從一般反應中提取出與目標相關的反應，試次選取“目標刺激”，“非目標刺激”置為零。然后，為了與fNIRS的時間進程一致，將EEG數據的每個block后補充5行零。最后，將試次間的變化與HRF卷積，得到每種負荷任務的數據格式為200行 × 160列。
2.6 ERP的AMs
為了證明n-back范式中存在ERP的AMs現象，本文使用預處理后的EEG數據，對所有通道選取ERP關鍵的N200（225 ms）和P300（350 ms）位置進行計算。在3種負荷任務下，將單個block的5個目標提取出來，再對10個block進行疊加平均，得到單個受試者的平均振幅變化。然后，在3種負荷任務下分別疊加平均所有受試者的AMs。
2.7 fNIRS數據預處理和特征提取
根據原始近紅外光強文件，使用修正的比爾 ?朗伯特定律計算出HbO濃度[20]。fNIRS數據的預處理使用MATLAB R2019a（MathWorks., 美國）和柏林BCI工具箱[21]。首先，為了與試次間的間隔保持一致，對fNIRS數據降采樣至2.2 Hz。然后，對HbO濃度數據進行0.01～0.2 Hz帶通濾波和0.1 Hz陷波，去除了心跳、呼吸噪聲和邁耶波[22]。最后，將數據按照標記基線校正，取0～44 s的數據做分析。
預處理后的數據標準化為零平均值和單位方差。在每個試次后，血氧動力學響應與神經元活動相比有延遲。因此，fNIRS特征以2.2 Hz采樣率取任務期15個采樣點和任務期后的5個采樣點，得到每種負荷任務的數據格式為200行 × 20列。
2.8 特征融合
將特征提取后的兩個特征矩陣輸入到CCA中，輸入數據都是200行，EEG數據每行有160個采樣點，fNIRS數據每行20個通道。經過CCA計算后，得到了兩種模態的相關系數和相關矩陣。用最小二乘法擬合相關矩陣得到EEG數據時間分量和fNIRS數據空間分量。為了提高該方法的魯棒性，EEG數據和fNIRS數據輸入CCA使用了留一交叉驗證。因此，CCA執行了23次，并將23個結果疊加平均。EEG數據與fNIRS數據融合分析的步驟如圖3所示。
圖3
EEG數據與fNIRS數據融合分析的步驟
Figure3.
The steps of fusion analysis of EEG data and fNIRS data
3 結果
3.1 ERP的AMs
經過通道間對比后，發現Cz、FC1和FC2存在試次間的AMs，而其它位置AMs不明顯。因此，本文只展示Cz、FC1和FC2的平均AMs結果。如圖4所示，0-back負荷下N200的AMs擬合曲線的斜率為0.06、截距為 ? 0.40；P300的AMs擬合曲線的斜率為 ? 0.01、截距為1.26。1-back負荷下N200的AMs擬合曲線的斜率為0.12、截距為 ? 1.00；P300的AMs擬合曲線的斜率為 ? 0.40、截距為3.08。2-back負荷下N200的AMs擬合曲線的斜率為0.06、截距為 ? 0.12；P300的AMs擬合曲線的斜率為 ? 0.18、截距為2.35。總之，在0-back、1-back和2-back中均觀察到N200的上升趨勢和P300的下降趨勢。
圖4
3種負荷下N200和P300的AMs結果
Figure4.
AMs results of N200 and P300 under three types of loads
3.2 相關系數
在CCA計算后，本文得到20對相關成分。為了進一步提高數據的可靠性，使用統計學分析去篩選相關成分。選取統計學分析的P值小于0.01，得到了3對0-back相關成分、4對1-back相關成分和2對2-back相關成分。統計學分析篩選后的相關系數如表1所示。對于0-back負荷，3對相關成分的相關系數分別是0.85、0.80和0.77。對于1-back負荷，4對相關成分的相關系數分別是0.82、0.80、0.77和0.74。對于2-back負荷，2對相關成分的相關系數分別是0.83、0.80。
表1
EEG數據與fNIRS數據間的相關系數
Table1.
Correlation coefficient between EEG data and fNIRS data
3.3 0-back負荷下EEG數據與fNIRS數據的相關成分
為了觀察EEG數據時間分量和fNIRS數據空間分量的相關變化，將相關矩陣A進行以下處理：首先，對10個block疊加平均；然后，為了獲得所有受試者的平均幅度，把不同受試者的數據進行疊加平均，得到試次間變化如圖5、6和7第一列所示。fNIRS數據空間分量如圖5、6和7第二列所示，右側的標尺代表激活程度。如圖5、6和7第三列所示，由于標準化和卷積的處理，EEG數據的行已經沒有規律；因為EEG數據的平均響應被移除，所以EEG數據時間分量已經沒有ERP峰值特征[23]。為了進一步探究EEG數據時間分量的信息，EEG數據時間分量的功率譜密度（power spectral density，PSD）如圖5、6和7第四列所示。
圖5
0-back負荷下EEG數據與fNIRS數據的相關成分
Figure5.
Correlation components of EEG data and fNIRS data under 0-back loads
圖6
1-back負荷下EEG數據與fNIRS數據的相關成分
Figure6.
Correlation components of EEG data and fNIRS data under 1-back loads
0-back負荷下EEG數據與fNIRS數據的相關成分如圖5所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有顯著的fNIRS激活；EEG數據時間分量的PSD主要包含delta頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和theta頻帶。對于成分3，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta、theta和alpha頻帶。
3.4 1-back負荷下EEG數據與fNIRS數據的相關成分
1-back負荷下EEG數據與fNIRS數據的相關成分如圖6所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和alpha頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含theta頻帶。對于成分3，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和theta頻帶。對于成分4，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活，且右腦激活比左腦激活高；EEG數據時間分量的PSD主要包含delta頻帶。
3.5 2-back負荷下EEG數據與fNIRS數據的相關成分
2-back負荷下EEG數據與fNIRS數據的相關成分如圖7所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有顯著的fNIRS激活；EEG數據時間分量的PSD主要包含delta、theta和alpha頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活，且右腦激活比左腦激活高；EEG數據時間分量的PSD主要包含delta和theta頻帶。
圖7
2-back負荷下EEG數據與fNIRS數據的相關成分
Figure7.
Correlation components of EEG data and fNIRS data under 2-back loads
4 討論
本文使用CCA計算工作記憶期間EEG數據和fNIRS數據的相關性，進而結合EEG數據和fNIRS數據的優勢探究工作記憶的神經血管耦合機制。對于工作記憶，CCA算法可以提取出EEG數據時間分量和fNIRS數據空間分量，從時間和空間兩個方面分析工作記憶期間大腦的變化。
本文發現3種負荷任務在P300處AMs的下降，N200處AMs的上升，這種隨ERP時間序列變化的AMs與fNIRS激活有密切關系。本文使用CCA算法融合EEG數據和fNIRS數據，并使用最小二乘法擬合出EEG數據時間分量和fNIRS數據空間分量。由于大腦在疲勞狀態下會觀察到delta節律[24]。所以，EEG數據時間分量的PSD中delta節律的出現，推斷是受試者任務期間勞累的表現。通常theta節律與長期記憶有關[25]。在EEG數據時間分量的PSD中，theta節律同樣與長期記憶有密切關系。已有研究發現，alpha節律在調節人體注意力方面具有重要作用[26]。在EEG數據時間分量的PSD中，alpha節律的出現可能與任務期間受試者的注意力有關系。在fNIRS數據空間分量中，額極區在3種負荷任務都有fNIRS激活。結果表明，額極區在工作記憶中起關鍵作用，與以往的研究一致[27]。在1-back的成分4和2-back的成分2中，觀察到額極區的fNIRS激活偏側性現象，且2-back的偏側性更為顯著。因此，推測在工作記憶中右腦額極區抑制了左腦額極區的類似功能。
綜合以上EEG數據時間分量和fNIRS數據空間分量的討論，得出額極區和背外側前額葉的fNIRS激活與EEG數據的delta、theta和alpha節律相關。本研究的CCA算法成功地將EEG數據和fNIRS數據的特性融合到一起，在融合算法上具有一定的創新性，并揭示了工作記憶下的神經血管耦合機制。
雖然本文在工作記憶下多模態的耦合取得了一定的成果，但是依然存在一些不足。第一，未充分利用EEG信號的空間信息。受限于當前的CCA方法，本文只采用了EEG信號三個通道的均值。未來將進一步改進擴展當前的CCA方法，實現大空間尺度上的EEG數據和fNIRS數據融合。第二，未考慮不同性別在工作記憶任務下的差異[28]。未來將招募女性受試者，更全面地研究工作記憶任務下的神經血管耦合機制。
5 結論
本文利用n-back范式對大腦神經活動和血氧活動的影響，研究基于EEG數據和fNIRS數據聯合采集下的神經血管耦合機制。本文通過CCA算法發現在n-back任務中，EEG信號delta、theta和alpha節律的出現，相應額極區和背外側前額葉的fNIRS激活，為工作記憶的雙模態時空分析提供了依據，特別在EEG數據和fNIRS數據融合方面提供了新穎的融合方法。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：劉文正是本研究的實驗設計者和實驗研究的執行人，完成數據采集、數據分析、論文寫作；張昊參與本研究的數據分析；楊柳參與本研究的實驗設計；古悅是本研究的構思者及負責人，指導實驗設計、數據分析、論文寫作與修改。
倫理聲明：本研究通過了北京師范大學倫理委員會的審批（批文編號：IRB_A_0009_2021003）。
引言
近年來，認知任務的研究越來越多，高級認知功能已經成為認知神經學、教育學和臨床醫學等諸多學科的重要話題。工作記憶是指人的大腦在執行認知任務過程中，用于大腦的信息加工和短時存儲過程的一個容量有限的系統[1]。工作記憶是高級認知功能的重要基礎，對于認知科學的研究具有重要意義[2]。隨著腦功能成像技術的不斷發展，在工作記憶時人類大腦的變化可以通過腦成像設備檢測，例如腦電圖（electroencephalography，EEG）和功能近紅外成像（functional near-infrared spectroscopy，fNIRS）。EEG數據是通過電極將神經元自發、節律性的運動電位記錄下來而得到的數據，具有很高的時間分辨率[3]。但是，EEG設備無法深入到大腦皮層以下采集信號，空間分辨率較低。有研究使用EEG設備進行溯源分析，發現音樂會降低記憶學習的效率[4]。在工作記憶任務中，學習障礙兒童的EEG數據具有更多的delta和theta節律，以及gamma節律的減少[5]。fNIRS技術利用血液對近紅外光的散射性測量出氧合血紅蛋白（oxygenated hemoglobins，HbO）和脫氧血紅蛋白（deoxygenated hemoglobins，HbR）濃度，具有價格低、無創性和易攜帶等特點，受到腦科學研究的高度重視[6]。基于fNIRS技術，有研究表明適度的運動有利于老年人進行工作記憶任務[7]。在工作記憶任務時，使用fNIRS技術發現左撇子和右撇子有不同的大腦皮層間功能連接[8]。
由于不同的采集設備具有不同的特性，因此單一采集設備的研究會具有一定局限性。基于神經血管耦合機制，融合多種設備來達到一種設備無法比擬的時空分辨率，可以從不同角度研究大腦功能[9]。在國內外神經血管耦合的研究中，EEG數據和功能磁共振成像（functional magnetic resonance imaging，fMRI）數據融合的居多。有研究使用EEG數據與fMRI數據探究睡眠和覺醒狀態下的神經血管耦合穩定性[10]。事件相關電位（event related potentials，ERP）間的振幅調制（amplitude modulations，AMs）與fMRI大腦激活也具有關聯性[11]。雖然，有很多學者研究EEG數據與fMRI數據融合，但是在融合EEG數據與fMRI數據時存在很多問題。例如，fMRI采集過程中要求受試者不能動，否則會對采集的數據造成偽影，影響分析的結果；fMRI釋放的磁場會影響EEG數據與fMRI數據聯合采集的效果。所以，如何解決兩者融合過程中出現的問題也成為了相關研究的熱點。fNIRS是一種新興的大腦信號測量設備，它和fMRI都是基于腦功能血氧動力學的測量方法。相比于fMRI，fNIRS可以在運動的環境下使用，且與EEG設備沒有不良干擾。因此，EEG數據和fNIRS數據的融合具有很大的研究意義和應用前景。在無創腦機接口（brain computer interface，BCI）研究領域，有研究證實EEG數據和fNIRS數據的多模態融合提升了BCI的性能[12]。
在已有的融合方法中，基于數據驅動的方法不需要對復雜的神經元種群和神經血管耦合動力學建模[13]。基于數據驅動的方法有偏最小二乘法（partial least squares，PLS）和典型相關分析（canonical correlation analysis，CCA）等一系列經典算法。因為CCA能夠更好地應用到多模態、多受試者的數據[14]，所以本文使用CCA計算EEG數據試次間的AMs和fNIRS數據的最大相關性，并提取出有關聯的EEG數據時間分量和fNIRS數據空間分量，有效促進了兩者的融合。綜上所述，使用CCA算法融合EEG數據和fNIRS數據，可以從時間和空間兩個方面研究大腦的神經活動。
本文采用工作記憶中經典范式n-back，通過采集設備獲得同步EEG數據和fNIRS數據。首先對采集到的兩種模態的數據進行預處理和特征提取。對特征提取后的矩陣，使用CCA計算模態間相關系數以及對應的EEG數據時間分量和fNIRS數據空間分量。通過對EEG數據時間分量和fNIRS數據空間分量的分析，得到了與工作記憶相關的結果[15]。
1 EEG數據和fNIRS數據融合理論
EEG數據與fNIRS數據融合模型如圖1所示。EEG數據矩陣經過預處理后的行是單個試次的采樣點，列是試次間的變化。然后，將預處理后的矩陣與血氧動力學響應函數（hemodynamic response function，HRF）卷積得到特征矩陣X。HRF是描述fNIRS對神經元瞬間刺激的響應函數。當大腦的某些區域進行神經活動時，氧氣消耗量增加并且刺激血管擴張以增加血流量來提供更多的HbO，達到給神經元活動提供氧氣的目的。在神經元刺激后，HbO濃度會增加，HbR濃度會下降，這種現象形成了血氧動力學響應特征Ri(t)，其計算公式如式（1）所示[16]：
圖1
EEG數據與fNIRS數據融合模型
Figure1.
The fusion model of EEG data and fNIRS data
|
式中，自變量t是時間，exp是以e為底的指數函數，i = 1代表HbO模型，i = 2代表HbR模型。對于HbO模型，C = 8log2、τ1,1 = 7.8、w1,1 = 5.95、τ1,2 = 12.16、w1,2 = 16.35和r1=0.37；對于HbR模型，C = 8log2、τ2,1 = 8.5、w2,1 = 7.3、τ2,2 = 21.4、w2,2 = 20.31和r2 = 0.331。
fNIRS數據矩陣經過預處理和特征提取后，得到特征矩陣Y。特征矩陣Y的行是所有通道，列是在時間上的采樣。特征矩陣X的試次間變化和Y的時間上的采樣相互對應，都是隨時間變化的數據采樣。特征矩陣X和Y的融合模型如式（2）所示：
|
式中，，，，，，，其中R代表實數集。TX、TY是矩陣X和Y的時間維度（數據的采樣點），VX、VY是矩陣X和Y的特征變量，DX、DY是矩陣X和Y的最小秩。相關矩陣AX和AY中第i列向量為AiX和AiY（i=1，2，，D），且同一相關矩陣AX或AY的列之間不相關；不同相關矩陣AX和AY的列之間相關，且相乘等于相關系數。相關矩陣AX和AY可由CCA計算得到。
CCA算法整體思想：對于X矩陣，存在投影向量aT將X投影到一維；對于Y矩陣，存在投影向量bT將Y投影到一維；投影后的一維向量分別是X′和Y′，那么可以表示成X′ = aTX和Y′ = bTY，并計算矩陣X′和Y′間的最大相關性[17]。
如式（3）所示，CCA的優化目標是最大化相關系數λ，并算出對應的投影向量a和b。
| '/> |
其中，SXY是矩陣X和Y的互協方差矩陣，SXX是矩陣X的協方差矩陣，SYY是矩陣Y的協方差矩陣。公式（3）的限定條件是aTSXXa = 1，bTSYYb = 1。然后，使用拉格朗日特征分解方法優化公式（3）得到的特征方程組如式（4）所示：
|
解特征方程組（4）得到最大相關系數λ1和第1對投影向量a1和b1。
當得到前k ? 1對投影向量ai和bi（i = 1，，k ? 1）后，第k對投影向量的限定條件，就如式（5）所示：
|
然后，再次用拉格朗日方法優化公式（3）得到相關系數λk和第k對投影向量ak和bk。aiX和biY是X和Y的第i對典型變量Ai。如式（6）所示，不同相關矩陣Ai之間相關，且相乘等于相關系數λi。
|
如式（7）所示，因為模型與最小二乘法原理相同，所以可以用最小二乘法擬合出模態的相關成分Ci。
|
2 材料和方法
2.1 受試者
本研究共有23名成年人參與數據采集，全部為男性（平均年齡22歲，年齡范圍18～26歲）。所有受試者都是右利手，且都身體健康、沒有心理或精神類疾病。每個受試者都在任務前簽署了受試者知情同意書。本研究經北京師范大學研究倫理審查委員會根據赫爾辛基宣言批準。
2.2 n-back范式
本研究采用數字設計的n-back范式，包括3種負荷任務：0-back、1-back和2-back。不同負荷任務的目標刺激不同。0-back負荷時，目標數字刺激是“7”，其他數字刺激是非目標刺激。1-back負荷時，如果當前數字和前1個數字相同則為目標刺激，不同則為非目標刺激。2-back負荷時，如果當前數字和前2個數字相同則為目標刺激，不同則為非目標刺激。以上所述“目標刺激”需要受試者按下鍵盤空格鍵。本范式中每種負荷任務包含10個組塊（block），3種不同的負荷任務按隨機順序顯示。n-back范式任務流程如圖2所示，以1-back為例：首先，屏幕上顯示5 s的任務類型“1-back”、1 s的注視窗口“+”；然后，每2.2 s隨機出現1個數字，每個數字出現0.5 s，接著是1.7 s的注視窗口“+”，每個block共有15個數字；最后，每個block的任務結束后有20 s的休息時間。在每個block中，“目標刺激”出現的概率為三分之一，“非目標刺激”概率為三分之二。因此，本范式共有3×10×15=450個試次。
圖2
n-back范式
Figure2.
n-back paradigm
2.3 EEG數據獲取
使用多通道EEG設備（ActiCHamp，Brain Products.，德國）以1 000 Hz的采樣率記錄連續的任務數據。EEG數據采集電極共32個（AFp1、AFp2、F7、AFF1h、AFF2h、F8、FC5、FC1、FC2、FC6、T7、C3、C4、T8、TP9、CP5、CP1、CP2、CP6、TP10、P7、P3、Pz、P4、P8、AFF5h、O1、Oz、O2、AFF6h、Cz和FCz），其中，Cz是參考電極，FCz是接地電極，所有電極的阻抗保持在5 kΩ以下。
2.4 fNIRS數據獲取
使用fNIRS設備（NIRSport，NIRx Medical Technologies.，美國）以7.812 5 Hz的采樣率記錄連續的任務數據。相鄰的光源和探測器配置為fNIRS通道，使用8個光源和7個探測器在EEG電極（AF7、Fp1、Fpz、Fp2、AF8、AF3、AFz、AF4、F5、F3、F1、Fz、F2、F4、F6）周圍形成了20個通道。20個通道的光源和探測器距離均設置為30 mm[18]。fNIRS設備的光源、探測器和EEG設備的電極被固定在同一個電極帽上。
2.5 EEG數據預處理和特征提取
使用EEG數據處理軟件Analyzer 2.2（Brain Products.，德國）對采集到的EEG數據進行預處理。首先，對EEG數據進行200 Hz降采樣，并對所有EEG電極平均重參考。然后，對EEG數據進行1～40 Hz濾波、50 Hz陷波、偽跡檢測、壞通道替換。本文使用獨立分量分析（independent component analysis，ICA）識別并去除眼電噪聲。最后，將數據按照標記基線校正，取0～800 ms的數據做分析，并對每個試次進行小波去噪[19]。
預處理后的數據標準化為零平均值和單位方差。選取AMs現象明顯的額中央電極（Cz、FC1和FC2）平均值作為EEG特征。接下來分別對3種負荷任務下的EEG數據提取特征。首先，從一般反應中提取出與目標相關的反應，試次選取“目標刺激”，“非目標刺激”置為零。然后，為了與fNIRS的時間進程一致，將EEG數據的每個block后補充5行零。最后，將試次間的變化與HRF卷積，得到每種負荷任務的數據格式為200行 × 160列。
2.6 ERP的AMs
為了證明n-back范式中存在ERP的AMs現象，本文使用預處理后的EEG數據，對所有通道選取ERP關鍵的N200（225 ms）和P300（350 ms）位置進行計算。在3種負荷任務下，將單個block的5個目標提取出來，再對10個block進行疊加平均，得到單個受試者的平均振幅變化。然后，在3種負荷任務下分別疊加平均所有受試者的AMs。
2.7 fNIRS數據預處理和特征提取
根據原始近紅外光強文件，使用修正的比爾 ?朗伯特定律計算出HbO濃度[20]。fNIRS數據的預處理使用MATLAB R2019a（MathWorks., 美國）和柏林BCI工具箱[21]。首先，為了與試次間的間隔保持一致，對fNIRS數據降采樣至2.2 Hz。然后，對HbO濃度數據進行0.01～0.2 Hz帶通濾波和0.1 Hz陷波，去除了心跳、呼吸噪聲和邁耶波[22]。最后，將數據按照標記基線校正，取0～44 s的數據做分析。
預處理后的數據標準化為零平均值和單位方差。在每個試次后，血氧動力學響應與神經元活動相比有延遲。因此，fNIRS特征以2.2 Hz采樣率取任務期15個采樣點和任務期后的5個采樣點，得到每種負荷任務的數據格式為200行 × 20列。
2.8 特征融合
將特征提取后的兩個特征矩陣輸入到CCA中，輸入數據都是200行，EEG數據每行有160個采樣點，fNIRS數據每行20個通道。經過CCA計算后，得到了兩種模態的相關系數和相關矩陣。用最小二乘法擬合相關矩陣得到EEG數據時間分量和fNIRS數據空間分量。為了提高該方法的魯棒性，EEG數據和fNIRS數據輸入CCA使用了留一交叉驗證。因此，CCA執行了23次，并將23個結果疊加平均。EEG數據與fNIRS數據融合分析的步驟如圖3所示。
圖3
EEG數據與fNIRS數據融合分析的步驟
Figure3.
The steps of fusion analysis of EEG data and fNIRS data
3 結果
3.1 ERP的AMs
經過通道間對比后，發現Cz、FC1和FC2存在試次間的AMs，而其它位置AMs不明顯。因此，本文只展示Cz、FC1和FC2的平均AMs結果。如圖4所示，0-back負荷下N200的AMs擬合曲線的斜率為0.06、截距為 ? 0.40；P300的AMs擬合曲線的斜率為 ? 0.01、截距為1.26。1-back負荷下N200的AMs擬合曲線的斜率為0.12、截距為 ? 1.00；P300的AMs擬合曲線的斜率為 ? 0.40、截距為3.08。2-back負荷下N200的AMs擬合曲線的斜率為0.06、截距為 ? 0.12；P300的AMs擬合曲線的斜率為 ? 0.18、截距為2.35。總之，在0-back、1-back和2-back中均觀察到N200的上升趨勢和P300的下降趨勢。
圖4
3種負荷下N200和P300的AMs結果
Figure4.
AMs results of N200 and P300 under three types of loads
3.2 相關系數
在CCA計算后，本文得到20對相關成分。為了進一步提高數據的可靠性，使用統計學分析去篩選相關成分。選取統計學分析的P值小于0.01，得到了3對0-back相關成分、4對1-back相關成分和2對2-back相關成分。統計學分析篩選后的相關系數如表1所示。對于0-back負荷，3對相關成分的相關系數分別是0.85、0.80和0.77。對于1-back負荷，4對相關成分的相關系數分別是0.82、0.80、0.77和0.74。對于2-back負荷，2對相關成分的相關系數分別是0.83、0.80。
表1
EEG數據與fNIRS數據間的相關系數
Table1.
Correlation coefficient between EEG data and fNIRS data
3.3 0-back負荷下EEG數據與fNIRS數據的相關成分
為了觀察EEG數據時間分量和fNIRS數據空間分量的相關變化，將相關矩陣A進行以下處理：首先，對10個block疊加平均；然后，為了獲得所有受試者的平均幅度，把不同受試者的數據進行疊加平均，得到試次間變化如圖5、6和7第一列所示。fNIRS數據空間分量如圖5、6和7第二列所示，右側的標尺代表激活程度。如圖5、6和7第三列所示，由于標準化和卷積的處理，EEG數據的行已經沒有規律；因為EEG數據的平均響應被移除，所以EEG數據時間分量已經沒有ERP峰值特征[23]。為了進一步探究EEG數據時間分量的信息，EEG數據時間分量的功率譜密度（power spectral density，PSD）如圖5、6和7第四列所示。
圖5
0-back負荷下EEG數據與fNIRS數據的相關成分
Figure5.
Correlation components of EEG data and fNIRS data under 0-back loads
圖6
1-back負荷下EEG數據與fNIRS數據的相關成分
Figure6.
Correlation components of EEG data and fNIRS data under 1-back loads
0-back負荷下EEG數據與fNIRS數據的相關成分如圖5所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有顯著的fNIRS激活；EEG數據時間分量的PSD主要包含delta頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和theta頻帶。對于成分3，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta、theta和alpha頻帶。
3.4 1-back負荷下EEG數據與fNIRS數據的相關成分
1-back負荷下EEG數據與fNIRS數據的相關成分如圖6所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和alpha頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含theta頻帶。對于成分3，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區和背外側前額葉具有微弱的fNIRS激活；EEG數據時間分量的PSD主要包含delta和theta頻帶。對于成分4，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活，且右腦激活比左腦激活高；EEG數據時間分量的PSD主要包含delta頻帶。
3.5 2-back負荷下EEG數據與fNIRS數據的相關成分
2-back負荷下EEG數據與fNIRS數據的相關成分如圖7所示。對于成分1，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有顯著的fNIRS激活；EEG數據時間分量的PSD主要包含delta、theta和alpha頻帶。對于成分2，EEG數據與fNIRS數據在試次間有相近的變化趨勢；在大腦皮層額極區具有微弱的fNIRS激活，且右腦激活比左腦激活高；EEG數據時間分量的PSD主要包含delta和theta頻帶。
圖7
2-back負荷下EEG數據與fNIRS數據的相關成分
Figure7.
Correlation components of EEG data and fNIRS data under 2-back loads
4 討論
本文使用CCA計算工作記憶期間EEG數據和fNIRS數據的相關性，進而結合EEG數據和fNIRS數據的優勢探究工作記憶的神經血管耦合機制。對于工作記憶，CCA算法可以提取出EEG數據時間分量和fNIRS數據空間分量，從時間和空間兩個方面分析工作記憶期間大腦的變化。
本文發現3種負荷任務在P300處AMs的下降，N200處AMs的上升，這種隨ERP時間序列變化的AMs與fNIRS激活有密切關系。本文使用CCA算法融合EEG數據和fNIRS數據，并使用最小二乘法擬合出EEG數據時間分量和fNIRS數據空間分量。由于大腦在疲勞狀態下會觀察到delta節律[24]。所以，EEG數據時間分量的PSD中delta節律的出現，推斷是受試者任務期間勞累的表現。通常theta節律與長期記憶有關[25]。在EEG數據時間分量的PSD中，theta節律同樣與長期記憶有密切關系。已有研究發現，alpha節律在調節人體注意力方面具有重要作用[26]。在EEG數據時間分量的PSD中，alpha節律的出現可能與任務期間受試者的注意力有關系。在fNIRS數據空間分量中，額極區在3種負荷任務都有fNIRS激活。結果表明，額極區在工作記憶中起關鍵作用，與以往的研究一致[27]。在1-back的成分4和2-back的成分2中，觀察到額極區的fNIRS激活偏側性現象，且2-back的偏側性更為顯著。因此，推測在工作記憶中右腦額極區抑制了左腦額極區的類似功能。
綜合以上EEG數據時間分量和fNIRS數據空間分量的討論，得出額極區和背外側前額葉的fNIRS激活與EEG數據的delta、theta和alpha節律相關。本研究的CCA算法成功地將EEG數據和fNIRS數據的特性融合到一起，在融合算法上具有一定的創新性，并揭示了工作記憶下的神經血管耦合機制。
雖然本文在工作記憶下多模態的耦合取得了一定的成果，但是依然存在一些不足。第一，未充分利用EEG信號的空間信息。受限于當前的CCA方法，本文只采用了EEG信號三個通道的均值。未來將進一步改進擴展當前的CCA方法，實現大空間尺度上的EEG數據和fNIRS數據融合。第二，未考慮不同性別在工作記憶任務下的差異[28]。未來將招募女性受試者，更全面地研究工作記憶任務下的神經血管耦合機制。
5 結論
本文利用n-back范式對大腦神經活動和血氧活動的影響，研究基于EEG數據和fNIRS數據聯合采集下的神經血管耦合機制。本文通過CCA算法發現在n-back任務中，EEG信號delta、theta和alpha節律的出現，相應額極區和背外側前額葉的fNIRS激活，為工作記憶的雙模態時空分析提供了依據，特別在EEG數據和fNIRS數據融合方面提供了新穎的融合方法。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：劉文正是本研究的實驗設計者和實驗研究的執行人，完成數據采集、數據分析、論文寫作；張昊參與本研究的數據分析；楊柳參與本研究的實驗設計；古悅是本研究的構思者及負責人，指導實驗設計、數據分析、論文寫作與修改。
倫理聲明：本研究通過了北京師范大學倫理委員會的審批（批文編號：IRB_A_0009_2021003）。
