黃 亮 吳超華 高小榕

(清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084)

基于自回歸模型和相位斜率指數(shù)的動(dòng)態(tài)腦連接分析

黃 亮 吳超華 高小榕#*

(清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084)

隨著腦網(wǎng)絡(luò)理論的發(fā)展，人們越來越關(guān)注不同腦區(qū)域之間的功能性和效應(yīng)性連接。在常見的腦活動(dòng)無創(chuàng)檢測方法中，腦電圖(EEG)具有較高的時(shí)間分辨率，適于進(jìn)行效應(yīng)性連接分析。提出一種估計(jì)不同通道的EEG信號間的效應(yīng)性連接的方法，即自回歸相位斜率指數(shù)(AR-PSI)。該方法結(jié)合多元自回歸(MVAR)模型對短時(shí)數(shù)據(jù)進(jìn)行譜估計(jì)頻率分辨率高的特點(diǎn)和相位斜率指數(shù)(PSI)對源信號混疊不敏感的特點(diǎn)。與傳統(tǒng)的格蘭杰因果模型相比，它可以有效地摒除由于容積導(dǎo)體效應(yīng)造成的信號混疊所帶來的干擾；與傳統(tǒng)的PSI相比，它在短時(shí)數(shù)據(jù)上能夠更準(zhǔn)確地估計(jì)不同通道間EEG信號的效應(yīng)性連接。首先，分別生成具有強(qiáng)效應(yīng)性連接和相互獨(dú)立的混疊噪聲這兩組信號進(jìn)行模擬實(shí)驗(yàn)，結(jié)果表明AR-PSI方法確實(shí)能夠更有效地檢測出信號中存在的效應(yīng)連接，排除信號線性混疊可能引起的誤檢。然后，應(yīng)用此方法并結(jié)合滑動(dòng)窗技術(shù)，對Stroop實(shí)驗(yàn)記錄到的EEG數(shù)據(jù)進(jìn)行網(wǎng)絡(luò)連接的動(dòng)態(tài)分析，發(fā)現(xiàn)250～500 ms和550～800 ms時(shí)間段內(nèi)腦網(wǎng)絡(luò)連接密度在兩種條件下存在顯著性差異。分析結(jié)果，顯示語義和顏色不一致條件的刺激能夠引起腦網(wǎng)絡(luò)連接密度更迅速地增加，且連接跨度更大。

效應(yīng)性連接；MVAR模型；相位斜率指數(shù)(PSI)；動(dòng)態(tài)分析

引言

近年來，腦網(wǎng)絡(luò)理論[1-8]越來越受到人們的關(guān)注，產(chǎn)生了大量的研究成果。要構(gòu)建一個(gè)腦網(wǎng)絡(luò)，關(guān)鍵步驟之一是如何定義節(jié)點(diǎn)與節(jié)點(diǎn)之間的連接。根據(jù)研究方法的不同，腦連接可以分為結(jié)構(gòu)連接[9]、功能連接[10]和效應(yīng)連接[11]三大類。結(jié)構(gòu)連接描述不同腦區(qū)在物理上的聯(lián)系，目前常用的主要是通過解剖和醫(yī)學(xué)影像技術(shù)觀察和計(jì)算不同大腦區(qū)域的物理連接或者形態(tài)相似性；功能連接表示不同大腦區(qū)域的功能活動(dòng)在時(shí)間上的統(tǒng)計(jì)相關(guān)性，目前主要是根據(jù)腦電圖(EEG)、腦磁圖(MEG)和功能核磁共振(fMRI)測量的信號來進(jìn)行統(tǒng)計(jì)計(jì)算；效應(yīng)連接則反映不同腦區(qū)之間的有向信息傳輸機(jī)制，估計(jì)效應(yīng)連接需要的數(shù)據(jù)與估計(jì)功能連接所用的數(shù)據(jù)基本相同，不同之處主要在于連接的統(tǒng)計(jì)計(jì)算方法上。目前，計(jì)算效應(yīng)連接的方法大多是基于參數(shù)模型的方法，包括結(jié)構(gòu)方程模型、多變量回歸模型、動(dòng)態(tài)因果模型、格蘭杰因果模型等。

格蘭杰因果(Granger-causal, GC)模型[12]是目前比較常用并且比較成熟的估計(jì)變量間因果關(guān)系的方法，在基于腦電信號的效應(yīng)連接分析中也占有一席之地。該方法的基本思想是：如果加入第二個(gè)變量的信息能夠使對第一個(gè)變量的預(yù)測更加準(zhǔn)確，那么就認(rèn)為第二個(gè)變量是第一個(gè)變量的因。但是，在實(shí)際的多通道頭皮腦電信號采集中，由于容積導(dǎo)體效應(yīng)，所觀察到的腦電信號實(shí)際上是多個(gè)活動(dòng)源信號混疊的結(jié)果，GC模型在應(yīng)用于此類信號時(shí)，就會(huì)估計(jì)出虛假的效應(yīng)連接?；诖?，Nolte等提出了相位斜率指數(shù)(phase slope index，PSI)[13]。該方法基于多通道信號的交叉譜相位，利用其相位的斜率，能夠確定不同頻率帶上信號的時(shí)間先后關(guān)系，從而確定信號之間的效應(yīng)連接。后來，Haufe等通過仿真實(shí)驗(yàn)對GC和PSI進(jìn)行了較為嚴(yán)格的評估，證明PSI在避免因信號混疊引起的虛假效應(yīng)連接估計(jì)方面確實(shí)有較強(qiáng)的優(yōu)勢[14]。實(shí)際上，PSI方法計(jì)算的準(zhǔn)確性與多通道信號交叉譜頻率分辨率密切相關(guān)。Nolte等采用一般的非參數(shù)方法估計(jì)多通道信號的交叉譜，當(dāng)數(shù)據(jù)長度較短時(shí)，估計(jì)所得譜的頻率分辨率低，使實(shí)際計(jì)算出的PSI值準(zhǔn)確度降低。為了克服這個(gè)問題，本研究提出了自回歸相位斜率指數(shù)(autoregressive phase slope index，AR-PSI)。該方法通過引入多元自回歸模型(MVAR)[15]，把MVAR和PSI方法的優(yōu)勢結(jié)合起來，能更加準(zhǔn)確有效地估計(jì)EEG信號各個(gè)通道之間的效應(yīng)連接。

1 方法和實(shí)驗(yàn)

1.1 AR-PSI方法

AR-PSI方法實(shí)際上是AR模型和PSI方法的結(jié)合，它充分利用了AR模型譜估計(jì)頻率分辨率高的優(yōu)勢和PSI方法對源信號混疊不敏感的特點(diǎn)，能夠更準(zhǔn)確地估計(jì)多通道信號之間的效應(yīng)連接，下面將詳細(xì)論述AR-PSI方法的原理和步驟。

1.1.1 MVAR模型擬合

MVAR是AR模型在多維變量上的擴(kuò)展。通常，MVAR模型可以用方程表示為

(1)

式中：xn∈m表示m維平穩(wěn)時(shí)間序列；A1,A2,…,Ap∈m×m是MVAR模型的系數(shù)矩陣;un表示均值為零、協(xié)方差矩陣為Σ∈m×m的不相關(guān)噪聲向量；v∈m是非零均值時(shí)間序列的截距項(xiàng)。

為了方便計(jì)算，通常在對時(shí)間序列進(jìn)行去均值處理后，把截距項(xiàng)設(shè)為零。因此，MVAR模型實(shí)際需要估計(jì)的未知參數(shù)就只是Ak和Σ。

有文獻(xiàn)報(bào)道了多種MVAR模型的參數(shù)估計(jì)方法(即模型擬合方法)，包括最小二乘類方法(如MLS、ARFIT[16]及LWR[17])、點(diǎn)陣類方法(如Vieira-Morf[18])、狀態(tài)空間模型類方法(Kalmanfiltering[19])等，這些方法的性能可參考文獻(xiàn)[20]。在這些方法中，MLS方法(多維最小二乘法)最為簡單，而且性能不差，所以本研究采用此方法進(jìn)行MVAR模型的擬合。

在MVAR模型擬合之前，需要指定模型的階數(shù)p。階數(shù)的選擇與模型的擬合效果有密切關(guān)系：階數(shù)選得過小，不足以充分利用觀測數(shù)據(jù)的信息進(jìn)行準(zhǔn)確擬合；階數(shù)選得過大，需要估計(jì)的未知量增多，容易出現(xiàn)過擬合，且估計(jì)所需計(jì)算量也大大增加。一般情況下，可以使用一種或多種準(zhǔn)則來選擇模型階數(shù)，常用的準(zhǔn)則包括SBC(Schwarz-Bayescriterion)準(zhǔn)則、AIC(Akaikeinformationcriterio)準(zhǔn)則、FPE(Akaike′sfinalpredictionerror)準(zhǔn)則和HQ(Hannan-Quinncriterio)準(zhǔn)則等[15]。對于小樣本數(shù)據(jù)來說，AIC準(zhǔn)則和FPE準(zhǔn)則具有更好的效果[15]。本研究的實(shí)驗(yàn)數(shù)據(jù)樣本量有限，所以選用AIC準(zhǔn)則[21]進(jìn)行階數(shù)選擇，該準(zhǔn)則的計(jì)算公式為

(2)

式中，Σ(p)表示p階模型擬合誤差的協(xié)方差矩陣，N表示用于模型擬合的總數(shù)據(jù)點(diǎn)數(shù)。

MVAR模型擬合完成之后，需要對該模型進(jìn)行統(tǒng)計(jì)上的檢驗(yàn)，以確保擬合的結(jié)果是可信的。通常，檢驗(yàn)內(nèi)容包括3個(gè)方面：擬合度檢驗(yàn)、穩(wěn)定性檢驗(yàn)和一致性檢驗(yàn)。

擬合度檢驗(yàn)旨在檢驗(yàn)擬合系數(shù)矩陣是否完全描述了原始數(shù)據(jù)的內(nèi)部關(guān)聯(lián)結(jié)構(gòu)，該檢驗(yàn)等效為檢驗(yàn)擬合系數(shù)矩陣用于原始數(shù)據(jù)上所得殘差是否是白噪聲。一般的做法是：檢驗(yàn)一定遲滯h以內(nèi)的殘差自相關(guān)系數(shù)是否足夠小，以確保在一定置信水平上不能拒絕殘差是白噪聲的假設(shè)。在多變量自回歸模型的擬合驗(yàn)證中，采用混成檢驗(yàn)法(portmanteautests)[22]得到的結(jié)果更加可信?；斐蓹z驗(yàn)法定義了一系列的檢驗(yàn)統(tǒng)計(jì)量，這些統(tǒng)計(jì)量包括Box-Pierce(BPP)、Ljung-Box(LBP)、Li-Mcleod(LMP)等，它們在殘差是白噪聲的零假設(shè)下近似服從自由度為m2(h-p)的χ2分布。選用其中一個(gè)統(tǒng)計(jì)量LMP，表示為

(3)

式中，Cl表示遲滯l時(shí)的殘差自相關(guān)系數(shù)。

穩(wěn)定性檢驗(yàn)即檢驗(yàn)系統(tǒng)是否能夠收斂，可以通過判斷系數(shù)矩陣構(gòu)成的穩(wěn)定性矩陣的特征值是否都在單位圓內(nèi)來實(shí)現(xiàn)。穩(wěn)定性矩陣的結(jié)構(gòu)為

(4)

穩(wěn)定指數(shù)可定義為

(5)

式中，λmax是穩(wěn)定性矩陣A絕對值最大的特征值。當(dāng)SI<0時(shí)，該MVAR模型是穩(wěn)定的。

一致性檢驗(yàn)即判斷真實(shí)數(shù)據(jù)擬合結(jié)果和擬合得到的模型重新生成的仿真數(shù)據(jù)擬合結(jié)果的相似程度。Ding等[23]定義了一致性百分比指數(shù)的計(jì)算式，即

(6)

式中，Rr和Rs分別表示真實(shí)數(shù)據(jù)和由擬合模型產(chǎn)生的模擬數(shù)據(jù)的自相關(guān)矩陣。

PC值越接近于100%，表明擬合模型越能夠產(chǎn)生和原始數(shù)據(jù)具有相同關(guān)聯(lián)結(jié)構(gòu)的數(shù)據(jù)。

1.1.2 相干譜估計(jì)

一旦擬合的模型通過上述檢驗(yàn)，就可以在此基礎(chǔ)上進(jìn)行相干譜估計(jì)。對式(1)的兩邊同時(shí)進(jìn)行傅里葉變換，把MVAR模型轉(zhuǎn)換為頻域形式，可得

(7)

式中，H(f)被稱為系統(tǒng)轉(zhuǎn)移方程，其定義如下：

(8)

據(jù)此，很容易得到其譜矩陣，有

(9)

式中，H*(f)是H(f)的共軛轉(zhuǎn)置。

以此方法得到的譜分辨率高，不受計(jì)算樣本數(shù)據(jù)長短的限制。傳統(tǒng)的非參數(shù)方法為保證譜的頻率分辨率需要更長的數(shù)據(jù)，所以當(dāng)數(shù)據(jù)長度有限時(shí)，用MVAR模型的方法比用傳統(tǒng)非參數(shù)方法能夠得到更精細(xì)的譜分析結(jié)果，有利于更精確地計(jì)算PSI。

1.1.3PSI估計(jì)

PSI基于相干譜相位的斜率，定義如下：

(10)

歸一化的相干譜定義如下：

(11)

通過適當(dāng)?shù)淖儞Q，式(10)可以寫成如下便于理解的形式，有

(12)

當(dāng)該相位譜的頻率分辨率足夠高時(shí)，有

(13)

因此，PSI實(shí)際表示的是相干譜相位斜率的加權(quán)平均，權(quán)系數(shù)為各個(gè)頻率上的歸一化相干譜幅度。

從PSI的定義來看，PSI取得的是相干譜的虛部，而相干譜的虛部信息并不會(huì)因?yàn)樵葱盘栔g的混疊而發(fā)生變化[24]。所以，PSI能夠有效地避免源信號混疊所帶來的效應(yīng)連接誤估計(jì)。從式(13)可以看出，δf越小(頻率分辨率越高)，等式左右越趨近于相等，PSI值越準(zhǔn)確。在傳統(tǒng)的PSI估計(jì)方法中，運(yùn)用非參譜估計(jì)方法得到的譜頻率分辨率受信號長度限制，而筆者采用MVAR模型進(jìn)行譜估計(jì)，大大提高了相干譜的頻率分辨率，能夠更準(zhǔn)確地估計(jì)出PSI。

1.2 仿真實(shí)驗(yàn)

為了證實(shí)AR-PSI方法在測量不同信號之間的效應(yīng)連接方面的優(yōu)越性，首先將該方法在一些仿真數(shù)據(jù)上進(jìn)行測試。仿真數(shù)據(jù)包括兩組二維的時(shí)間序列，每組的采樣點(diǎn)數(shù)均為2 000點(diǎn)，采樣率為100 Hz。其中一組數(shù)據(jù)是用如下所示的一階MVAR模型生成的，即

(14)

另一組是相互獨(dú)立的白噪聲和粉紅噪聲的線性混疊，即

(15)

對仿真信號先進(jìn)行分段，每段長度為400點(diǎn)，然后分別用GC方法、PSI方法和AR-PSI方法對信號間的效應(yīng)連接進(jìn)行估計(jì)。按照上述方法生成2000組仿真信號，求取3種指標(biāo)的平均值和方差，如圖1(c)、(d)所示。

從圖1所示的結(jié)果中可以看出：對于有強(qiáng)效應(yīng)連接的兩個(gè)信號(MVAR模型生成的仿真信號)，3種方法都比較好地估計(jì)出了信號間的效應(yīng)連接，但是AR-PSI方法估計(jì)結(jié)果的均值顯著大于其他兩種方法(配對t檢驗(yàn)，P<0.5)，而其估計(jì)結(jié)果的方差也顯著小于其他兩種方法(F檢驗(yàn)，P<0.5)。結(jié)果表明，對于實(shí)際存在效應(yīng)連接的信號，AR-PSI方法能夠估計(jì)出更大且更穩(wěn)定的結(jié)果。而對于線性混疊的獨(dú)立噪聲信號，GC方法誤檢出了效應(yīng)連接，其結(jié)果顯著大于0(t檢驗(yàn)，P<0.5)，傳統(tǒng)PSI方法和AR-PSI方法估計(jì)結(jié)果的均值沒有顯著性差異(配對t檢驗(yàn)，P=0.28)，但是，AR-PSI方法估計(jì)結(jié)果的方差卻顯著小于傳統(tǒng)PSI方法的方差(F檢驗(yàn)，P<0.5)。結(jié)果表明，AR-PSI方法能夠最有效地消除效應(yīng)連接的誤檢，且比傳統(tǒng)PSI方法更穩(wěn)定。綜上所述，在參與比較的3種方法中，AR-PSI方法能夠更有效地檢測出信號中存在的效應(yīng)連接，排除信號線性混疊可能引起的誤檢。

圖1 數(shù)據(jù)仿真實(shí)驗(yàn)。 (a)二維一階AR模型生成的信號；(b) 相互獨(dú)立的白噪聲與粉紅噪聲混疊的信號；(c) GC方法、PSI方法和AR-PSI方法分別對(a)中信號估計(jì)所得效應(yīng)連接的均值和方差；(d) GC方法、PSI方法和AR-PSI方法分別對(b)中信號估計(jì)所得效應(yīng)連接的均值和方差Fig.1 Simulation results. (a) Signal generated by 2 dimensional 1 order AR model; (b)Mixture signal of independent white and pink noise; (c) Effective connectivity estimated by Granger Causality, PSI, AR-PSI using data (a). The height of the bar indicates the mean and the error bar indicates the standard deviation; (d) Effective connectivity estimated by Granger Causality, PSI, AR-PSI using data (b). The height of the bar indicates the mean and the error bar indicates the standard deviation

1.3 真實(shí)EEG信號的動(dòng)態(tài)腦連接分析

仿真數(shù)據(jù)上的實(shí)驗(yàn)，證明了AR-PSI方法相對GC方法及傳統(tǒng)PSI方法在估計(jì)腦效應(yīng)連接上的優(yōu)越性。接下來，結(jié)合滑動(dòng)窗技術(shù)和AR-PSI方法，利用真實(shí)EEG信號進(jìn)行動(dòng)態(tài)腦連接分析。

1.3.1 數(shù)據(jù)采集

真實(shí)EEG數(shù)據(jù)是在漢字語義顏色Stroop實(shí)驗(yàn)[25]中采集到的數(shù)據(jù)。該實(shí)驗(yàn)包含語義和顏色一致、語義和顏色不一致以及語義和顏色無關(guān)3種刺激條件，每一例受試在3種條件下各采集得到120個(gè)試次的EEG數(shù)據(jù)。實(shí)驗(yàn)階段共采集了15例受試，其中第14例受試信號受眼電污染嚴(yán)重，信號極不平穩(wěn)，所以沒有進(jìn)行分析。在實(shí)驗(yàn)時(shí)，受試位于密閉的屏蔽室，采用Synamps2 (Neuroscan Inc.)腦電系統(tǒng)，記錄60個(gè)通道的腦電數(shù)據(jù)，電極位置如圖2所示。采樣率為1 000 Hz，以頭頂為參考。

圖2 電極位置及分區(qū)Fig.2 Electrodes position and subarea

為了方便后續(xù)的分析，對這60個(gè)電極進(jìn)行分區(qū)，并用不同的色彩標(biāo)識。其中，F(xiàn)P1、FPZ、FP2、AF3、AF4這5個(gè)電極定義為F區(qū)域(前額葉區(qū))，F(xiàn)7、F5、F3、F1、FZ、FC1、FC3、FC5、FT7這9個(gè)電極定義為LF區(qū)域(左側(cè)額葉區(qū))，F(xiàn)2、F4、F6、F8、FT8、FC6、FC4、FC2、FCZ這9個(gè)電極定義為RF區(qū)域(右側(cè)額葉區(qū))，T7、C5、C3、C1、CZ、CP1、CP3、CP5、TP7這9個(gè)電極定義為LC區(qū)域(左側(cè)中央?yún)^(qū))，C2、C4、C6、T8、TP8、CP6、CP4、CP2、CPZ這9個(gè)電極定義為RC區(qū)域(右側(cè)中央?yún)^(qū))，P7、P5、P3、P1、PZ、PO3、PO5、PO7這8個(gè)電極定義為LP區(qū)域(左側(cè)頂葉區(qū))，P2、P4、P6、P8、PO8、PO6、PO4、POZ這8個(gè)電極定義為RP區(qū)域(右側(cè)頂葉區(qū))，O1、OZ、O2這3個(gè)電極定義為O區(qū)域(枕葉區(qū))。*該分區(qū)只是為了方便描述而對電極進(jìn)行的粗略分區(qū)，它們與解剖上的大腦分區(qū)并不一定存在嚴(yán)格的一一對應(yīng)關(guān)系。

1.3.2 數(shù)據(jù)預(yù)處理

在進(jìn)行分析前，實(shí)驗(yàn)采集的原始EEG數(shù)據(jù)都需要經(jīng)過特定的預(yù)處理過程,主要分為以下兩大步驟。

1)第一大步驟的預(yù)處理過程。預(yù)處理的目的在于獲取可用于腦連接分析的數(shù)據(jù)源：首先，進(jìn)行1～40 Hz的帶通濾波，消除漂移噪聲和高頻噪聲；其次，降采樣到200 Hz，其目的在于增強(qiáng)數(shù)據(jù)的平穩(wěn)性以適于進(jìn)行MVAR模型擬合；再次，截取刺激前1 200 ms到刺激后1 000 ms的數(shù)據(jù)用作分析。經(jīng)過這3項(xiàng)處理，每一例受試和每一種條件都分別得到一個(gè)60×440×120(通道數(shù)×數(shù)據(jù)點(diǎn)數(shù)×試次數(shù))的EEG數(shù)據(jù)矩陣，這便是后續(xù)進(jìn)行腦連接分析的數(shù)據(jù)源。

2)第二大步驟的預(yù)處理過程[23]。采用AR-PSI方法進(jìn)行腦連接分析，一個(gè)很重要的步驟就是MVAR模型的擬合，而MVAR模型擬合必須保證數(shù)據(jù)源在時(shí)間上的平穩(wěn)性。實(shí)際上，帶有ERP的EEG信號卻包含非平穩(wěn)成分，因此筆者通過第二大步驟的預(yù)處理過程來提高數(shù)據(jù)的平穩(wěn)性。首先，對上一步驟獲得的數(shù)據(jù)源進(jìn)行歸一化處理，即減去其時(shí)域平均值，并除以其時(shí)域標(biāo)準(zhǔn)差，保證每一試次和每一通道的信號幅度近似一致；其次，去除信號中的ERP成分(非平穩(wěn)成分)，即用上述處理后的信號減去其總體平均。經(jīng)過這兩步處理后，還可以用矩形窗截取小段數(shù)據(jù)，因?yàn)閿?shù)據(jù)越短，其平穩(wěn)性越好。在有多個(gè)試次的情況下，即使是只有小段數(shù)據(jù)，MVAR模型也能夠較好地?cái)M合其模型參數(shù)。經(jīng)過試驗(yàn)，窗長為250 ms時(shí)，數(shù)據(jù)能夠比較好地進(jìn)行MVAR模型的擬合。為了連續(xù)地觀察不同時(shí)間段內(nèi)腦連接的變化情況，還引入了滑動(dòng)窗技術(shù)，窗長為250 ms(即50個(gè)數(shù)據(jù)點(diǎn))，窗滑動(dòng)步長為50 ms(即10個(gè)數(shù)據(jù)點(diǎn))。

1.3.3 腦連接估計(jì)

經(jīng)過以上預(yù)處理之后，每個(gè)短窗中的數(shù)據(jù)都可用于MVAR模型的擬合。首先，根據(jù)AIC準(zhǔn)則選擇一個(gè)合適的模型階數(shù)。通過設(shè)定模型階數(shù)依次為1～20，根據(jù)式(2)計(jì)算出每一階數(shù)對應(yīng)的AIC值，可得到一條單調(diào)下降的曲線。該曲線顯示，隨著模型階數(shù)的增長，AIC值逐漸下降，并且其下降速率也逐漸減小，經(jīng)過某個(gè)比較明顯的拐點(diǎn)之后，AIC值就趨于平穩(wěn)。對于每一段窗數(shù)據(jù)，都能夠得到類似的曲線，其拐點(diǎn)基本上都位于階數(shù)10左右。為了保持對所有數(shù)據(jù)處理方法的一致，方便之后做不同條件和不同時(shí)間下的結(jié)果對比，對所有窗中的數(shù)據(jù)都采用60個(gè)變量的10階MVAR模型進(jìn)行擬合。

完成模型擬合之后，采用多變量混成檢驗(yàn)方法，檢驗(yàn)由各個(gè)窗數(shù)據(jù)的擬合模型得到的剩余噪聲是否為白噪聲。在剩余噪聲為白噪聲的原假設(shè)條件下，應(yīng)用替代數(shù)據(jù)法構(gòu)造隨機(jī)的1 000例剩余噪聲樣本，并依據(jù)式(3)分別計(jì)算每一個(gè)樣本的Li-McLeod統(tǒng)計(jì)量，就可以得到該統(tǒng)計(jì)量的一個(gè)分布，在此分布的基礎(chǔ)上計(jì)算P值。經(jīng)過計(jì)算，對于各段窗數(shù)據(jù)MVAR擬合模型的統(tǒng)計(jì)檢驗(yàn)結(jié)果都有P>0.05，因而可判定在統(tǒng)計(jì)意義上這些擬合模型的剩余噪聲均為白噪聲。

同時(shí)，根據(jù)式(5)、(6)驗(yàn)證各擬合模型的穩(wěn)定性和一致性。穩(wěn)定性矩陣A絕對值的最大特征值的對數(shù)值均小于0，說明估計(jì)得到的模型都是穩(wěn)定的;模型一致性的計(jì)算結(jié)果也均在70%左右，可以認(rèn)為具有良好的一致性。綜合其檢驗(yàn)結(jié)果，可以判定各個(gè)窗數(shù)據(jù)擬合得到的MVAR模型均能夠很好地反映原始數(shù)據(jù)的內(nèi)在結(jié)構(gòu)，這些擬合結(jié)果是可信的。

圖3 兩種條件及不同時(shí)間段AR-PSI值統(tǒng)計(jì)直方圖Fig.3 Histogram of AR-PSI value in all time windows and conditions

1.3.4 腦連接分析

AR-PSI描述的是不同信號之間的因果聯(lián)系，所以其值有正有負(fù)，而得到的也是一個(gè)個(gè)有向網(wǎng)絡(luò)?；趫D理論的復(fù)雜網(wǎng)絡(luò)有許多分析指標(biāo)，選取網(wǎng)絡(luò)連接密度，其定義式為

(16)

式中，E表示網(wǎng)絡(luò)的邊數(shù)，V表示網(wǎng)絡(luò)的節(jié)點(diǎn)數(shù)。

網(wǎng)絡(luò)連接密度與網(wǎng)絡(luò)平均度意義等同，可反映當(dāng)前腦網(wǎng)絡(luò)的連接數(shù)，間接反映當(dāng)前腦活動(dòng)的強(qiáng)弱。

統(tǒng)計(jì)各個(gè)時(shí)間段所構(gòu)建的平均腦網(wǎng)絡(luò)的連接密度，并采用t檢驗(yàn)的方法對兩種條件之間進(jìn)行統(tǒng)計(jì)檢驗(yàn)，找出存在顯著性差異的時(shí)間段。然后，對于存在顯著性差異的時(shí)間段的腦網(wǎng)絡(luò)，根據(jù)前面所劃分的8個(gè)區(qū)域，統(tǒng)計(jì)這些區(qū)域之間的連接方向和連接數(shù)，并利用重采樣的方法對兩種條件之間進(jìn)行統(tǒng)計(jì)檢驗(yàn)。最后，計(jì)算腦網(wǎng)絡(luò)各區(qū)域的出度(本區(qū)域節(jié)點(diǎn)是其他區(qū)域節(jié)點(diǎn)因的連接總數(shù))和入度(其他區(qū)域節(jié)點(diǎn)是本區(qū)域節(jié)點(diǎn)因的連接總數(shù))，進(jìn)一步分析腦網(wǎng)絡(luò)中的信息流向。

2 結(jié)果

為了驗(yàn)證Stroop實(shí)驗(yàn)數(shù)據(jù)的可靠性，給出了FCZ通道平均后的ERP波形，如圖4(b)所示。可以發(fā)現(xiàn)，兩種條件下的ERP波形在400～600 ms這個(gè)時(shí)間段內(nèi)存在顯著性差異(逐點(diǎn)配對t檢驗(yàn)，P<0.05)，且不一致條件下的波形比一致條件下的波形小，與已有文獻(xiàn)的研究結(jié)果相符[26，28]，表明所采集到的數(shù)據(jù)是可靠的。

圖4 腦連接動(dòng)態(tài)分析。 (a) 兩種條件下有顯著性差異時(shí)間段的腦網(wǎng)絡(luò)連接密度；(b) 兩種條件下FCZ電極的受試平均ERP波形(綠色覆蓋部分表示兩種條件下ERP波形存在顯著性差異的時(shí)間段)；(c) “語義與顏色”一致條件和不一致條件ERP差值在顯著性差異的時(shí)間段內(nèi)均值的空間分布；(d) “語義與顏色”一致條件和不一致條件ERP配對t檢驗(yàn)P值的空間分布Fig.4 Dynamic brain connectivity analysis. (a) Connectivity density with significant difference under two conditions; (b) Grand averaged ERP at electrode FCZ under two conditions (The region covered by green rectangle is the period with significant difference under two conditions); (c) Topographic map of the mean of ERP difference between two conditions in the periods with significant difference; (d) Topographic map of P value computed by pairwise t test

兩種條件下的腦網(wǎng)絡(luò)連接密度在250～500 ms (配對t檢驗(yàn)，P=0.047)及550～800 ms (配對t檢驗(yàn)，P=0.011)時(shí)間段有顯著性差異，見圖4(a)。刺激發(fā)生之后，兩種條件下的腦網(wǎng)絡(luò)連接密度均呈上升趨勢。從兩個(gè)有顯著性差異的時(shí)間段來看，在250～500 ms時(shí)間段內(nèi)，不一致條件下的腦網(wǎng)絡(luò)連接密度顯著高于一致條件下的腦網(wǎng)絡(luò)連接密度，而在550～800 ms時(shí)間段內(nèi)，一致條件下的腦網(wǎng)絡(luò)連接密度反而顯著高于一致條件的腦網(wǎng)絡(luò)連接密度。

兩種條件下腦網(wǎng)絡(luò)連接密度存在顯著性差異，其時(shí)間段的網(wǎng)絡(luò)連接如圖5所示。圖中節(jié)點(diǎn)表示的是EEG信號采集的各個(gè)通道，每個(gè)節(jié)點(diǎn)上都標(biāo)有對應(yīng)通道的標(biāo)號；連接表明兩個(gè)節(jié)點(diǎn)區(qū)域存在效應(yīng)連接，連接的一端加粗指示連接的方向性，即未加粗的一端影響加粗的一端。節(jié)點(diǎn)按照前后左右腦上的電極通道分別排布在圓圈的上下左右，區(qū)域與電極通道的對應(yīng)關(guān)系用不同的色彩標(biāo)出。從圖5可以直觀地看出，在連接方向上，不一致條件下主要是從LF區(qū)域(左側(cè)額葉區(qū))指向RP區(qū)域(右側(cè)頂葉區(qū))，而一致條件下主要是從RF區(qū)域(右側(cè)額葉區(qū))指向LC區(qū)域(左側(cè)中央?yún)^(qū))、RP區(qū)域(右側(cè)頂葉區(qū))指向RF區(qū)域(右側(cè)額葉區(qū))和LP區(qū)域(左側(cè)頂葉區(qū))。

圖5 腦網(wǎng)絡(luò)圖。 (a) 250～500 ms時(shí)間段內(nèi)“語義與顏色”不一致條件下構(gòu)建的腦網(wǎng)絡(luò)圖，網(wǎng)絡(luò)密度顯著大于對應(yīng)時(shí)間段內(nèi)“語義與顏色”一致條件下構(gòu)建的腦網(wǎng)絡(luò)圖；(b) 550～800 ms時(shí)間段內(nèi)“語義與顏色”一致條件下構(gòu)建的腦網(wǎng)絡(luò)圖，網(wǎng)絡(luò)密度顯著大于對應(yīng)時(shí)間段內(nèi)“語義與顏色”不一致條件下構(gòu)建的腦網(wǎng)絡(luò)圖Fig.5 Brain network. (a) Brain network under incongruent condition at 250-500ms, whose density is larger than that under congruent condition at the corresponding period. (b) Brain network constructed under congruent condition at 550-800ms, whose density is larger than that under incongruent condition at the corresponding period

對區(qū)域之間的連接方向和連接數(shù)進(jìn)行統(tǒng)計(jì)，得到兩種條件下存在顯著性差異的腦區(qū)域連接分布，如圖6所示。在250～500 ms時(shí)間段內(nèi)，不一致條件下LF區(qū)域(左側(cè)額葉區(qū))指向RP區(qū)域(右側(cè)頂葉區(qū))的連接數(shù)遠(yuǎn)遠(yuǎn)多于其他區(qū)域之間的連接數(shù)(見圖6(a))；在550～800 ms時(shí)間段內(nèi)，一致條件下RP區(qū)域(右側(cè)頂葉區(qū))指向RF區(qū)域(右側(cè)額葉區(qū))以及RP區(qū)域(右側(cè)頂葉區(qū))指向LP區(qū)域(左側(cè)頂葉區(qū))相對于其他區(qū)域之間存在較多的連接(見圖6(b))。

圖6 腦區(qū)域連接分布。 (a)250～500 ms時(shí)間段兩種條件下腦連接數(shù)具有顯著性差異的腦區(qū)域連接分布；(b)550～800 ms時(shí)間段兩種條件下腦連接數(shù)具有顯著性差異的腦區(qū)域連接分布Fig.6 Brain connectivity between different areas. (a) Brain connectivity between different areas with significant difference under two conditions in 250～500 ms; (b) Brain connectivity between different areas with significant difference under two conditions in 550～800 ms

表1 250～500 ms及550～800 ms時(shí)間段不同條件下各腦區(qū)域的出度和入度統(tǒng)計(jì)

Tab.1 Out-degree and in-degree of different brain areas under two condition in 250～500 ms and 550～800 ms

區(qū)域時(shí)間段250～500ms550～800ms一致Out/In不一致Out/In一致Out/In不一致Out/InF11/56/3524/133/12LF22/1663/942/1332/23RF14/412/2548/5254/0LC53/1222/2517/6814/39RC19/2615/274/6317/58LP20/1839/727/3935/15RP15/2820/4681/814/35O0/457/1030/1713/0

計(jì)算腦網(wǎng)絡(luò)各區(qū)域的出度和入度，得到表1所示結(jié)果。首先重點(diǎn)關(guān)注250～500 ms的不一致條件，結(jié)果顯示從LF區(qū)域(左側(cè)額葉區(qū))信息流出最多，而RP區(qū)域(右側(cè)頂葉區(qū))信息流入最多，相對于該時(shí)間段內(nèi)的一致條件下，LF區(qū)域(左側(cè)額葉區(qū))與RP區(qū)域(右側(cè)頂葉區(qū))有更強(qiáng)的效應(yīng)連接，文獻(xiàn)[29]也得出左額葉腦區(qū)與頂葉腦區(qū)在不一致條件下有更強(qiáng)相干性的結(jié)論。然后，再考察550～800 ms的一致條件，發(fā)現(xiàn)信息流出最多的區(qū)域是RP區(qū)域(右側(cè)頂葉區(qū))，信息流入最多的是LC區(qū)域(左側(cè)中央?yún)^(qū))，RC區(qū)域(右側(cè)中央?yún)^(qū))的信息流入量次之。兩者對比可以發(fā)現(xiàn)，不一致條件下大腦前后方向的大跨度信息流動(dòng)偏多，而一致條件下小跨度信息流動(dòng)偏多，并且傾向于左右方向之間的流動(dòng)。文獻(xiàn)[30-31]報(bào)道稱，Stroop效應(yīng)主要是前扣帶皮層(ACC)活動(dòng)差異引起的，不一致條件下ACC區(qū)域活動(dòng)增強(qiáng)，并同時(shí)影響腦前額葉和其他區(qū)域，因而使得腦前后連接增強(qiáng)，與本研究結(jié)果有一定的相似性。

3 討論

對于上述數(shù)據(jù)現(xiàn)象，結(jié)合人腦對信息的加工機(jī)理以及Stroop效應(yīng)的自動(dòng)化模型[32]加以討論。Stroop效應(yīng)的自動(dòng)化模型認(rèn)為：人腦對語義和顏色的加工采用不同的方式，對語義是自動(dòng)加工，對顏色是控制加工，當(dāng)語義和顏色不一致時(shí)，語義信息的自動(dòng)加工會(huì)對顏色信息的控制加工產(chǎn)生干擾，產(chǎn)生Stroop效應(yīng)。于是，可以做這樣的推測：當(dāng)干擾發(fā)生時(shí)，大腦可能需要立刻調(diào)用更多、更廣范圍內(nèi)的資源協(xié)同來對這種干擾進(jìn)行排除，因而在該條件下腦網(wǎng)絡(luò)連接數(shù)迅速增加，并且主要表現(xiàn)在前后腦大跨度的連接增多；而當(dāng)語義和顏色一致時(shí)，對于兩者的處理不會(huì)產(chǎn)生相互干擾，大腦遵循正常的信息處理模式，調(diào)用的資源可能是逐步增加，因而腦網(wǎng)絡(luò)連接數(shù)是持續(xù)緩慢地增加，并且主要表現(xiàn)為小跨度的連接模式。同時(shí)，筆者也注意到，連接性有顯著性差異的時(shí)間段和ERP信號波形有顯著性差異的時(shí)間段并不完全一致。腦電連接性和ERP信號所反映的認(rèn)知神經(jīng)過程并不完全一致，它們可能存在時(shí)間上的先后關(guān)系。但是，根據(jù)本研究的實(shí)驗(yàn)和數(shù)據(jù)并不能做出準(zhǔn)確的推測，因此將在今后對此開展更為細(xì)致的研究。

本研究所提出的AR-PSI方法在估計(jì)效應(yīng)連接的準(zhǔn)確性方面相比GC方法和傳統(tǒng)PSI方法確實(shí)有一定的優(yōu)勢。進(jìn)一步地，如果以該方法為依托能夠針對短時(shí)數(shù)據(jù)有效地構(gòu)建起腦網(wǎng)絡(luò)，將會(huì)有利于進(jìn)行動(dòng)態(tài)分析。但是，該方法在數(shù)據(jù)量很大的情況下，計(jì)算起來比較耗時(shí)。一方面，MVAR模型擬合花費(fèi)大量時(shí)間；另一方面，利用MVAR模型系數(shù)在全頻率段估計(jì)PSI，由于伴隨著大量的矩陣求逆運(yùn)算也相當(dāng)耗時(shí)。所以，該方法在計(jì)算速度上比不上可以直接在時(shí)域進(jìn)行運(yùn)算的GC方法。因此，該方法更適合用于對實(shí)時(shí)性要求不高的離線分析場合。

另外，腦電信號應(yīng)用腦網(wǎng)絡(luò)分析時(shí)，不同的參考選擇對其分析結(jié)果也有著不同的影響[33-34]。雖然AR-PSI方法能夠很好地消除容積導(dǎo)體效應(yīng)的影響，但是卻不能消除參考電極的影響。由于篇幅所限，本研究在應(yīng)用AR-PSI方法進(jìn)行腦網(wǎng)絡(luò)構(gòu)建時(shí)，直接選擇頭頂作為參考，沒有進(jìn)一步考慮其他的參考可能會(huì)給腦網(wǎng)絡(luò)分析的結(jié)果帶來何種影響，這也將會(huì)是后續(xù)研究工作中的一個(gè)重要方向。

4 結(jié)論

本研究利用MVAR模型在多通道信號譜估計(jì)方面頻率分辨率高的優(yōu)勢，并結(jié)合PSI方法，對獨(dú)立源信號混疊的不敏感性提出了AR-PSI來估計(jì)不同通道EEG信號之間的效應(yīng)連接，從理論上分析了它的可行性以及存在的優(yōu)勢。仿真結(jié)果也顯示，相比GC方法，它可以有效地防止由于信號混疊造成的虛警；相比傳統(tǒng)PSI方法，它能更準(zhǔn)確地估計(jì)效應(yīng)連接。在研究中，將AR-PSI應(yīng)用于Stroop實(shí)驗(yàn)采集的EEG數(shù)據(jù)，并基于AR-PSI進(jìn)行了動(dòng)態(tài)腦連接分析。分析結(jié)果顯示，語義和顏色不一致條件的刺激能夠引起腦網(wǎng)絡(luò)連接密度更迅速地增加，并且連接跨度更大，表明大腦在處理具有沖突的信息時(shí)可能會(huì)更快地調(diào)動(dòng)更多、更廣泛的資源。

(致謝 本課題得到清華大學(xué)醫(yī)學(xué)院神經(jīng)工程實(shí)驗(yàn)室老師和同學(xué)的大力支持，分析所用Stroop實(shí)驗(yàn)數(shù)據(jù)來源于清華大學(xué)神經(jīng)工程實(shí)驗(yàn)室高文靜同學(xué)，在此一并表示感謝。)

[1] Sporns O, Tononi G, K?tter R. The human connectome: a structural description of the human brain[J]. PLoS Computational Biology, 2005, 1(4): e42.

[2] 蔣田仔, 劉勇, 李永輝. 腦網(wǎng)絡(luò)：從腦結(jié)構(gòu)到腦功能[J]. 生命科學(xué), 2009(2):181-188.

[3] 梁夏, 王金輝, 賀永. 人腦連接組研究：腦結(jié)構(gòu)網(wǎng)絡(luò)和腦功能網(wǎng)絡(luò)[J]. 科學(xué)通報(bào), 2010(16): 1565-1583.

[4] 孫俊峰, 洪祥飛, 童善保. 復(fù)雜腦網(wǎng)絡(luò)研究進(jìn)展——結(jié)構(gòu)、功能、計(jì)算與應(yīng)用[J]. 復(fù)雜系統(tǒng)與復(fù)雜性科學(xué), 2011, 7(4): 74-90.

[5] Bullmore E, Sporns O. The economy of brain network organization[J]. Nature Reviews Neuroscience, 2012, 13(5): 336-349.

[6] Van Essen DC, Smith SM, Barch DM, et al. The WU-Minn human connectome project: an overview[J]. Neuroimage, 2013, 80: 62-79.

[7] Sporns O. Making sense of brain network data[J]. Nature Methods, 2013, 10(6): 491-493.

[8] Uehara T, Yamasaki T, Okamoto T, et al. Efficiency of a “small-world” brain network depends on consciousness level: a resting-state FMRI study[J]. Cerebral Cortex, 2014, 24(6): 1529-1539.

[9] Koch MA, Norris DG, Hund-Georgiadis M. An investigation of functional and anatomical connectivity using magnetic resonance imaging[J]. Neuroimage, 2002, 16(1): 241-250.

[10] Fingelkurts AA, Fingelkurts AA, K?hk?nen S. Functional connectivity in the brain—is it an elusive concept?[J]. Neuroscience & Biobehavioral Reviews, 2005, 28(8): 827-836.

[11] Horwitz B. The elusive concept of brain connectivity[J]. Neuroimage, 2003, 19(2): 466-470.

[12] Handbook of Time Series Analysis: Recent Theoretical Developments and Applications[M]. New York: John Wiley & Sons, 2006.

[13] Nolte G, Ziehe A, Nikulin VV, et al. Robustly estimating the flow direction of information in complex physical systems[J]. Physical Review Letters, 2008, 100(23): 234101.

[14] Haufe S, Nikulin VV, Müller KR, et al. A critical assessment of connectivity measures for EEG data: a simulation study[J]. Neuroimage, 2013, 64(1):120-133.

15] Lütkepohl H. New Introduction to Multiple Time Series Analysis[M]. Berlin: Springer Science & Business Media, 2007.

[16] Schneider T, Neumaier A. Algorithm 808: ARfit—A Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models[J]. ACM Transactions on Mathematical Software (TOMS), 2001, 27(1): 58-65.

[17] Morf M, Vieira A, Lee DTL, et al. Recursive multichannel maximum entropy spectral estimation[J]. IEEE Trans Geosci Electron, 1978, 16(2): 85-94.

[18] Marple Jr SL. Digital Spectral Analysis with Applications[M]. Englewood Cliffs: Prentice-Hall, 1987.

[19] Schl?gl A. The Electroencephalogram and the Adaptive Autoregressive Model: Theory and Applications [M]. Aachen: Shaker, 2000.

[20] Schl?gl A. A comparison of multivariate autoregressive estimators[J]. Signal Processing, 2006, 86(9): 2426-2429.

[21] Akaike H. A new look at the statistical model identification[J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

[22] Hosking JRM. The multivariate portmanteau statistic[J]. Journal of the American Statistical Association, 1980, 75(371): 602-608.

[23] Ding M, Bressler SL, Yang W, et al. Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment[J]. Biological Cybernetics, 2000, 83(1): 35-45.

[24] Nolte G, Bai O, Wheaton L, et al. Identifying true brain interaction from EEG data using the imaginary part of coherency[J]. Clinical Neurophysiology, 2004, 115(10): 2292-2307.

[25] 高文靜, 陳小剛, 高小榕. 基于 Hermite 濾波器的穩(wěn)態(tài)視覺誘發(fā)電位背景標(biāo)記 Stroop 效應(yīng)包絡(luò)分析[J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào), 2014, 31(1): 1-7.

[26] Qiu J, Luo Y, Wang Q, et al. Brain mechanism of Stroop interference effect in Chinese characters[J]. Brain Research, 2006, 1072(1): 186-193.

[27] Liotti M, Woldorff MG, Perez R, et al. An ERP study of the temporal course of the Stroop color-word interference effect[J]. Neuropsychologia, 2000, 38(5): 701-711.

[28] Zysset S, Müller K, Lohmann G, et al. Color-word matching Stroop task: separating interference and response conflict[J]. Neuroimage, 2001, 13(1): 29-36.

[29] Schack B, Chen ACN, Mescha S, et al. Instantaneous EEG coherence analysis during the Stroop task[J]. Clinical Neurophysiology, 1999, 110(8): 1410-1426.

[30] Harrison BJ, Shaw M, Yücel M, et al. Functional connectivity during Stroop task performance.[J]. NeuroImage, 2005, 24(1):181-191.

[31] Hanslmayr S, Past?tter B, B?uml KH, et al. The electrophysiological dynamics of interference during the stroop task[J]. Journal of Cognitive Neuroscience, 2008, 20(2):215-225.

[32] MacLeod CM. Half a century of research on the Stroop effect: an integrative review[J]. Psychological Bulletin, 1991, 109(2): 163-203.

[33] Marzetti L, Nolte G, Perrucci MG, et al. The use of standardized infinity reference in EEG coherency studies[J]. Neuroimage, 2007, 36(1):48-63.

[34] Yun Q, Peng X, Yao D. A comparative study of different references for EEG default mode network: the use of the infinity reference.[J]. Clinical Neurophysiology, 2010, 121(12):1981-1991.

Dynamic Brain Connectivity Analysis Based on Autoregressive Model and Phase Slope Index

Huang Liang Wu Chaohua Gao Xiaorong#*

(DepartmentofBiomedicalEngineering,MedicalSchool,TsinghuaUniversity,Beijing100084,China)

Functional and effective connectivity are important branches in brain network research. Electroencephalogram (EEG) has sufficient temporal resolution to catch the fast brain dynamic changes, so it is suitable for effective connectivity analysis. We proposed a new method to estimate the effective connectivity base on EEG, namely autoregressive phase slope index (AR-PSI). Spectral estimation based on MVAR has high frequency resolution even on short time data. PSI is insensitive to linear mixture of non-interacting sources. AR-PSI combined advantages of the two methods. Compared with conventional Granger causality model, AR-PSI could eliminate the interference caused by volume conduction. Compared with conventional PSI, AR-PSI could get more accurate estimation of effective connectivity with short time data. Experimental data indicated that AR-PSI could exactly detect the effective connectivity between two signals and exclude the false detection. AR-PSI was also applied to dynamic brain connectivity analysis based on EEG recorded in Stroop paradigm. We found that brain connectivity density had significant difference under two conditions in 250～500 ms and 550～800 ms. The results indicated that the incongruent stimulus could make the density of brain network increase more quickly and expand the span of connectivity.

effective connectivity; MVAR model; phase slope index( PSI); dynamic analysis

10.3969/j.issn.0258-8021. 2016. 01.001

2015-06-03， 錄用日期:2015-09-16

國家自然科學(xué)基金重點(diǎn)項(xiàng)目(61431007, 91320202, 91220301)

R318

A

0258-8021(2016) 01-0001-09

# 中國生物醫(yī)學(xué)工程學(xué)會(huì)會(huì)員(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: gxr-dea@tsinghua.edu.cn