龍雨涵 王 彬 薛 潔 杜 芬 劉 輝 熊 新

(1. 昆明理工大學(xué)信息工程與自動化學(xué)院，云南昆明 650500; 2. 云南警官學(xué)院信息網(wǎng)絡(luò)安全學(xué)院，云南昆明 650500)

1 引言

靜息態(tài)功能磁共振成像(Resting State-funcational Magnetic Resonance Imaging, rs-fMRI)是一種有效研究人腦功能和結(jié)構(gòu)的非介入檢測技術(shù)。該技術(shù)通過對血氧平衡依賴(blood oxygen level-dependent, BOLD)信號的分析來研究人腦的功能機制和神經(jīng)活動。研究人員一般可以通過對不同腦區(qū)之間BOLD信號相關(guān)性的分析來度量腦區(qū)之間的功能連接，從而構(gòu)建人腦功能網(wǎng)絡(luò)[1-3]。

目前大多數(shù)研究采用線性相關(guān)分析方法來構(gòu)建人腦功能網(wǎng)絡(luò)。Geerligs等人[4]用線性Pearson相關(guān)分析方法構(gòu)建腦區(qū)之間的功能連接，用于研究功能網(wǎng)絡(luò)特性隨年齡的變化。Ye等人[5]利用Pearson相關(guān)分析方法構(gòu)建人腦功能網(wǎng)絡(luò)，發(fā)現(xiàn)重度抑郁癥患者相對于正常人有更高的局部效率和模塊化結(jié)構(gòu)。Khan等人[6]用偏相關(guān)分析方法構(gòu)建青少年自閉癥患者的腦功能連接關(guān)系，發(fā)現(xiàn)在大腦感覺運動區(qū)域有著更強的功能連接。

但是由于大腦是一個復(fù)雜的系統(tǒng)，僅僅采用線性分析方法難以觀測其內(nèi)在特性，尤其是不同腦區(qū)之間相互關(guān)系隨時間變化的動態(tài)特性。已有研究表明人腦有著非常明顯的非線性特征，而血液動力學(xué)響應(yīng)函數(shù)也具有非線性特點[7-9]，僅通過線性相關(guān)分析方法來構(gòu)建大腦功能網(wǎng)絡(luò)并無法捕獲兩個腦區(qū)之間的非線性關(guān)系。所以在構(gòu)建動態(tài)功能腦網(wǎng)絡(luò)時，應(yīng)該更充分的考慮其非線性屬性。

很多研究人員已經(jīng)開始嘗試利用非線性的分析方法來構(gòu)建人腦功能網(wǎng)絡(luò)并展開相關(guān)研究。Han等人[10]通過非線性斯皮爾曼(Spearman)等級相關(guān)分析方法研究不同年齡與不同性別人群間腦功能網(wǎng)絡(luò)差異。董澤芹等人[11]用肯德爾(Kendall)相關(guān)分析方法來探究癲癇患者與正常人的腦功能網(wǎng)絡(luò)差異。牛寶東等人[12]基于希爾伯特黃變換(Hilbert Huang Transform，HHT)處理非線性腦信號，提高癲癇病的識別率。Salvador等人[13]采用基于相干分析的互信息方法研究發(fā)現(xiàn)功能連接主要存在于低頻上。蔣進航[14]通過實驗發(fā)現(xiàn)最大信息系數(shù)(Maximal Information Coefficient, MIC)適合應(yīng)用于腦網(wǎng)絡(luò)的相關(guān)性度量分析。蘇龍飛[15]利用MIC從rs-fMRI數(shù)據(jù)中提取非線性功能連接特征用于精神分裂癥的診斷。但是上述研究是將在數(shù)據(jù)采集期間內(nèi)的BOLD信號視為靜態(tài)，通過計算其平均值來分析和構(gòu)建一個腦區(qū)網(wǎng)絡(luò)，并沒有考慮BOLD信號隨時間變化的影響，因此不足以全面理解和描述人腦網(wǎng)絡(luò)的動態(tài)特性。也有一些人員開始針對大腦網(wǎng)絡(luò)的動態(tài)特性展開研究，以BOLD時間序列為對象來進行腦區(qū)之間的相關(guān)性分析，從而確定腦區(qū)之間的功能連接網(wǎng)絡(luò)，但是一般都是采用Pearson線性相關(guān)系數(shù)的分析方法構(gòu)建功能連接網(wǎng)絡(luò)。Yu等人[16]利用時變?nèi)四X功能連接的動態(tài)圖論理論去識別精神分裂癥的拓撲異常。Wee等人[17]利用功能動態(tài)網(wǎng)絡(luò)的稀疏特性用于輕度認知障礙的早期診斷。馬灑灑等人[18]將滑動窗口應(yīng)用于功能腦網(wǎng)絡(luò)的BOLD信號時間序列分析，從多維數(shù)據(jù)分析的角度研究人腦網(wǎng)絡(luò)的動態(tài)特征。

綜上所述，由于人腦的復(fù)雜性、實時變化特性以及BOLD信號的非線性特點，無論是采用線性分析方法展開的動態(tài)特性研究，還是用非線性分析方法展開的平均時間的分析方法都不足以表現(xiàn)腦區(qū)之間的真實動態(tài)屬性。因此本文以rs-fMRI原始數(shù)據(jù)提取到的BLOD信號時間序列為對象，討論了構(gòu)建全腦腦區(qū)功能網(wǎng)絡(luò)過程中的非線性相關(guān)分析方法及閾值確定等關(guān)鍵技術(shù)，采用基于滑動窗口的動態(tài)數(shù)據(jù)分析技術(shù)，構(gòu)建得到動態(tài)特征矩陣并進行特征降維，將其映射到二維空間內(nèi)可視化，在此基礎(chǔ)上對腦功能網(wǎng)絡(luò)的狀態(tài)及系統(tǒng)的動態(tài)特性展開進一步研究。將基于3種非線性相關(guān)分析方法斯皮爾曼相關(guān)分析(Spearman)、肯德爾相關(guān)分析(Kendall)和最大信息系數(shù)(Maximal Information Coefficient, MIC)的腦網(wǎng)絡(luò)動態(tài)特性的實驗對比結(jié)果顯示：在閾值參數(shù)選擇合理的前提下，無論對BOLD信號展開哪種非線性分析方法，所得到的腦網(wǎng)絡(luò)狀態(tài)及分類結(jié)果均顯現(xiàn)出相似的規(guī)律性，從而為基于非線性相關(guān)分析方法的動態(tài)腦網(wǎng)絡(luò)構(gòu)建方法提供了一種可行的解決方案，為下一步展開腦網(wǎng)絡(luò)的動態(tài)特性分析和演化過程研究奠定了基礎(chǔ)。

2 基于非線性相關(guān)分析方法的人腦動態(tài)特征矩陣構(gòu)建

2.1 數(shù)據(jù)來源及數(shù)據(jù)預(yù)處理方法說明

靜息態(tài)fMRI數(shù)據(jù)的采集使用美國通用電氣(GE)公司制造的3.0Tesla磁共振成像儀，使用頭部正交線圈，數(shù)據(jù)采集過程中要求被試者平躺、閉眼，保持清醒，全身禁止不動，盡量不進行特別的思維活動。采用梯度回波-回波平面成像序列，重復(fù)時間TR為1500 ms，回波時間TE為30 ms，采集矩陣為64×64，視野為FOV=256 mm×265 mm，掃描層后4 mm，層間距0 mm，翻轉(zhuǎn)角Filp angle為60°，共采集200個時間點。

本文使用Python/FSL Resting State Pipeline平臺[19]對rs-fMRI原始數(shù)據(jù)進行預(yù)處理，采用自動解剖圖集(Automated Anatomical Labeling，AAL)[20]將全腦分為90個腦區(qū)。對數(shù)據(jù)進行預(yù)處理的過程主要由以下7個步驟：(1)去除前4個時間點；(2)時間層校正；(3)頭動校正；(4)顱骨去除；(5)空間標準化；(6)帶通濾波；(7)提取腦區(qū)的平均BOLD 信號時間序列[18]。預(yù)處理后得到90個腦區(qū)在196采樣點上的BOLD信號時間序列。

2.2 基于非線性相關(guān)分析的全腦動態(tài)特征矩陣的構(gòu)建

在文獻[18]的研究基礎(chǔ)之上，以90個腦區(qū)BOLD信號的采樣時間軸作為基準，采用滑動窗口技術(shù)依次遍歷整個時間序列，每一個小窗口內(nèi)的BOLD信號時間序列形成當前采樣點的狀態(tài)觀測窗口[18]。本文滑動窗口的大小設(shè)置為n，窗口以步長l從左向右依次移動，即可將整體的BOLD信號時間序列劃分為g(g=196-n+l)個狀態(tài)觀測窗口。

以上述步驟得到的g個狀態(tài)觀測窗口為基礎(chǔ)，針對腦網(wǎng)絡(luò)在時間上的同步性和相關(guān)性展開分析和研究。本文分別選擇Spearman、Kendall和MIC 3種非線性相關(guān)分析方法構(gòu)建每個滑動窗口內(nèi)腦區(qū)之間的動態(tài)功能連接，每個觀測窗口內(nèi)均可以得到一個90×90的狀態(tài)觀測矩陣，該狀態(tài)觀測矩陣是沿對角線對稱的矩陣，在不考慮自相關(guān)的情況下，將腦區(qū)之間相關(guān)系數(shù)作為腦網(wǎng)絡(luò)的實時特征提取出來，一共有4005個特征。將相關(guān)系數(shù)矩陣對角線以上的特征依次取出首尾連接可得到1×4005維的腦網(wǎng)絡(luò)狀態(tài)觀測向量。通過以上方法每一個靜息態(tài)fMRI樣本在采樣區(qū)間內(nèi)將得到g個向量，即得到g×4005矩陣G。如果用矩陣G的行向量代表某一掃描時間區(qū)間內(nèi)腦區(qū)的功能連接特征，那么該矩陣能夠反映全腦網(wǎng)絡(luò)在全部數(shù)據(jù)采集時間上的動態(tài)變化過程，下文將該矩陣稱為動態(tài)特征矩陣。構(gòu)建滑動窗口和人腦動態(tài)特征矩陣的原理如下圖1所示。

圖1 滑動窗口和動態(tài)特征矩陣G的構(gòu)建原理圖Fig.1 The schematic diagram of constructing sliding window and dynamic feature matrix G

2.3 基于非線性相關(guān)分析的狀態(tài)矩陣構(gòu)建方法

在2.2中所構(gòu)建的觀測窗口的基礎(chǔ)上，本文分別采用三種相關(guān)分析方法Kendall、Spearman、MIC確定腦區(qū)i和腦區(qū)j之間的功能連接，從而構(gòu)建狀態(tài)矩陣，這三種非線性相關(guān)分析方法的原理及算法如下。

2.3.1 斯皮爾曼相關(guān)分析(Spearman)

Spearman相關(guān)系數(shù)與變量的分布和樣本容量都無關(guān)，只要兩個變量的觀測都是成對的等級評定資料，就可以用Spearman相關(guān)分析法進行研究[21]。在一個窗口內(nèi)兩個腦區(qū)的BOLD信號是成對的等級評定資料，所以可用Spearman來測定腦區(qū)之間的相關(guān)強度。

(1)

其中SRij∈[-1,1]。當一個腦區(qū)的BOLD信號隨著另一個腦區(qū)單調(diào)遞增的時候SRij=1，反之SRij=-1。

2.3.2 Kendall相關(guān)分析

Kendall相關(guān)分析法(Kendall coefficient concordance)是衡量等級變量相關(guān)程度的一個統(tǒng)計量[21]。在腦網(wǎng)絡(luò)研究中根據(jù)兩個腦區(qū)之間BOLD信號序?qū)Φ囊恢滦詠砼袛嗥湎嚓P(guān)性。

定義2 對于兩個時間序列長度為n的腦區(qū)，它們之間的Kendall系數(shù)KRij定義為：

(2)

其中，P表示一致的序?qū)€數(shù)。而且KRij∈[-1,1]。

2.3.3 最大信息系數(shù)(MIC)

Reshef等人[22]發(fā)表在《Science》上提出的最大信息系數(shù)MIC屬于最大化的基于信息的非參數(shù)性探索，可用于衡量兩個腦區(qū)之間非線性的強度。

定義3 令腦區(qū)i和腦區(qū)j的n個BOLD信號時間序列為數(shù)據(jù)集D，將數(shù)據(jù)集劃分為x行y列，則兩個腦區(qū)之間的最大信息系數(shù)MRij定義為：

MRij=max{M(D)x,y}

(3)

其中M(D)為D的特征矩陣，MRij∈[0,1]。具體公式定義可參考文獻[22]。

相對于前兩種非線性相關(guān)系數(shù)，最大信息系數(shù)能檢測出兩個腦區(qū)之間更復(fù)雜的非線性關(guān)系[22]。

對于基于滑動窗口技術(shù)而展開的腦區(qū)非線性相關(guān)性分析原理如圖 2 所示，在每一個觀測窗口內(nèi)構(gòu)建狀態(tài)矩陣的具體步驟如下：

(2)分別求取腦區(qū)i和腦區(qū)j的Kendall相關(guān)系數(shù)(KRij)、Spearman相關(guān)系數(shù)(SRij)、MIC相關(guān)系數(shù)(MRij)得到該時刻90×90的腦網(wǎng)絡(luò)狀態(tài)觀測矩陣；

(3)滑動窗口沿采樣時間方向變化，逐一按(1)、(2)中的步驟求取不同時刻的狀態(tài)觀測矩陣；

(4)得到全部采樣時刻的狀態(tài)觀測矩陣。

圖2 腦區(qū)相關(guān)性分析原理圖Fig.2 Schematic diagram of correlation analysis of brain regions

3 腦功能網(wǎng)絡(luò)的構(gòu)建及閾值確定方法

3.1 腦功能網(wǎng)絡(luò)的構(gòu)建準則和閾值范圍確定

利用相關(guān)分析方法構(gòu)建腦功能網(wǎng)絡(luò)時主要有兩方面因素需要考慮：確定腦網(wǎng)絡(luò)的節(jié)點及確定節(jié)點之間連接的邊[23]。本文通過ALL先驗?zāi)０逡呀?jīng)將整個大腦分為90個腦區(qū)，即腦網(wǎng)絡(luò)有90個節(jié)點。而節(jié)點之間的邊即腦區(qū)之間的功能連接則可以通過腦區(qū)之間196個BOLD信號時間序列的相關(guān)系數(shù)來確定。但無論采用上述哪種非線性相關(guān)分析算法構(gòu)建得到的動態(tài)特征矩陣G，在構(gòu)建網(wǎng)絡(luò)的過程中都需要進行閾值處理，其目的是使其符合網(wǎng)絡(luò)特性，只有當相關(guān)強度大小達到某一閾值時，才認為兩個腦區(qū)之間存在功能連接。但是目前對于確定腦功能網(wǎng)絡(luò)的閾值并沒有統(tǒng)一的標準，大多數(shù)研究者都是根據(jù)經(jīng)驗或參考其他文獻進行閾值選取。同時已有研究發(fā)現(xiàn)，在使用相關(guān)分析構(gòu)建人腦功能網(wǎng)絡(luò)時，不同的閾值會對所構(gòu)建的腦網(wǎng)絡(luò)拓撲性質(zhì)產(chǎn)生較大的影響[23]。所以合理選取閾值是根據(jù)動態(tài)特征矩陣G進行腦網(wǎng)絡(luò)構(gòu)建的關(guān)鍵。

研究發(fā)現(xiàn)腦功能網(wǎng)絡(luò)具有小世界性[24-25]和高效性[26]的特點，所以本節(jié)將根據(jù)大腦功能網(wǎng)絡(luò)完整性、小世界性和高效性的要求來討論不同相關(guān)分析方法所對應(yīng)的合理閾值范圍確定方法，上述要求分別可以做如下表述：

(1)大腦功能網(wǎng)絡(luò)應(yīng)該是一個完整的網(wǎng)絡(luò)，90個腦區(qū)中的每一個腦區(qū)應(yīng)該與其他腦區(qū)存在功能連接，所有網(wǎng)絡(luò)中節(jié)點的度應(yīng)該大于0。即ki>0,其中i=1,2,…,90。

為了觀測不同樣本之間的差別，本文采用6個健康人的rs-fMRI數(shù)據(jù)作為實驗樣本(S1,S2,S3,S4,S5,S6)，擬分別使用以上3種不同的非線性相關(guān)分析法計算出長度為196時間序列的相關(guān)系數(shù)矩陣，并用該閾值將相關(guān)系數(shù)矩陣二值化從而構(gòu)建人腦功能網(wǎng)絡(luò)。為了進行對比分析，除了使用以上3種非線性相關(guān)分析之外，還同時采用Pearson進行線性相關(guān)性分析。由于MIC相關(guān)分析方法對較長的數(shù)據(jù)對敏感，而且實驗結(jié)果也顯示，在對本文實驗數(shù)據(jù)中的196個時間序列BOLD信號進行MIC相關(guān)分析時，所得到的相關(guān)系數(shù)矩陣中的所有元素值普遍較小，因此不采用與其他三種相關(guān)分析相同的方法確定閾值范圍，具體方法見下文。如果構(gòu)建的人腦網(wǎng)絡(luò)全部滿足以上3個條件，則表明構(gòu)建該功能網(wǎng)絡(luò)時所用的閾值合理；反之，構(gòu)建腦功能網(wǎng)絡(luò)不符合其中的任意一個條件，則該閾值不合理。

對Pearson、Spearman和Kendall算法來說，相關(guān)系數(shù)值都在-1到1的范圍內(nèi)，在確定閾值范圍時具體步驟如下：

(1)以0為起始點，0.01為步長，依次提高閾值大小。并將閾值用于全時間序列相關(guān)系數(shù)矩陣進行二值化，即矩陣中元素絕對值大于閾值的元素置為1，小于閾值的元素置為0。

(2)將二值化后的矩陣構(gòu)建為人腦功能網(wǎng)絡(luò)，1代表腦區(qū)之間存在功能連接，0代表腦區(qū)之間沒有功能連接。若構(gòu)建的網(wǎng)絡(luò)滿足3.1節(jié)中的三個條件，則步驟(1)中的閾值為合理閾值。

(3)將閾值從0遍歷至1，最終確定閾值范圍。

通過以上3個步驟能夠得到6個樣本的合理閾值范圍如表1所示。

表1 三種相關(guān)分析方法的閾值范圍

由于人腦樣本的個體差異性，每一個樣本得到閾值范圍有所不同，但是通過表1中的實驗結(jié)果可以看到，在使用同種相關(guān)分析方法時，所有樣本的閾值范圍均有一定程度的重合。

3.2 小世界性驗證實驗及最優(yōu)閾值確定方法

大腦功能網(wǎng)絡(luò)是典型的小世界網(wǎng)絡(luò)，為了驗證腦功能網(wǎng)絡(luò)經(jīng)過閾值處理二值化后，腦功能網(wǎng)絡(luò)是否具有小世界性，本文選用σ指標進行判定。

σ是由Humphries等人[28]為衡量小世界性而提出的指標。σ指標的計算公式(4)如下：

σ=λ/γ

(4)

其中，λ是小世界網(wǎng)絡(luò)與隨機網(wǎng)絡(luò)的聚類系數(shù)比值，γ是小世界網(wǎng)絡(luò)與隨機網(wǎng)絡(luò)的平均最短路徑系長度的比值。相對于隨機網(wǎng)絡(luò)，小世界網(wǎng)絡(luò)有著較短的平均最短路徑長度和較高的聚類系數(shù)，所以同等規(guī)模(相同的節(jié)點數(shù)，邊數(shù)和度分布)的小世界網(wǎng)絡(luò)和隨機網(wǎng)絡(luò)，λ大于1，γ約等于1，即當σ大于1時腦功能網(wǎng)絡(luò)具有小世界性，而且σ越大，則腦網(wǎng)絡(luò)的小世界性越強。

以Spearman算法為例，本實驗中的六個樣本采用Spearman相關(guān)分析后，按照表1中的閾值范圍內(nèi)不同閾值處理所構(gòu)建得到的腦功能網(wǎng)絡(luò)的σ指標值如下圖3所示。

從圖3可以看出，合理閾值范圍內(nèi)不同閾值處理所構(gòu)建得到的腦功能網(wǎng)絡(luò)的σ值基本都大于1，并且隨著閾值的增大，腦網(wǎng)絡(luò)的小世界性越強。為了在這個范圍內(nèi)確定最優(yōu)的閾值，我們給出了一種將腦功能網(wǎng)絡(luò)的拓撲特性與小世界性相結(jié)合的綜合閾值確定方法。

圖3 合理閾值范圍內(nèi)的σ指標變化曲線Fig.3 The variation curve of σ value corresponding to different threshold value

已有研究表明，腦功能網(wǎng)絡(luò)信息傳遞效率與網(wǎng)絡(luò)聚類系數(shù)C和網(wǎng)絡(luò)平均最短路徑長度L密切相關(guān)，并且與C成正比，與L成反比[24]。本文中腦功能網(wǎng)絡(luò)聚類系數(shù)C和平均最短路徑長度L定義如下：

(5)

(6)

(7)

圖值變化曲線Fig.4 The variation curve of value

按照本文方法，最后得到的Spearman相關(guān)分析方法的最優(yōu)閾值為0.46，Kendall相關(guān)分析方法的最優(yōu)閾值為0.32，Pearson相關(guān)分析方法的最優(yōu)閾值為0.48。下文將使用這些最優(yōu)閾值展開動態(tài)特征矩陣的稀疏度實驗。

3.3 稀疏度實驗結(jié)果及對比分析

以上述6個樣本為對象，將滑動窗口大小設(shè)置為n=20，步長l=1，使用四種不同的相關(guān)分析方法分別構(gòu)建人腦動態(tài)特征矩陣G。在使用滑動窗口技術(shù)來研究腦網(wǎng)絡(luò)動態(tài)特性時，最重要的兩個參數(shù)就是窗口的大小和步長。但是目前研究者們對于這兩個參數(shù)的設(shè)定并沒有達成統(tǒng)一的共識。Shirer等人[29]研究發(fā)現(xiàn)人腦認知狀態(tài)的改變至少需要30～60 s；Li等人[30]研究發(fā)現(xiàn)窗口時間長度設(shè)定為20～40 s時，能夠較好的捕捉腦連接的動態(tài)變化。為了減小窗口大小對實驗的影響，本實驗選擇將窗口的時間長度設(shè)定為40 s,即窗口大小n為20。而且大部分研究者根據(jù)經(jīng)驗將窗口步長設(shè)定為1[16,18,31]，所以本實驗窗口的移動步長設(shè)定為1。

并針對動態(tài)特征矩陣G設(shè)定九種閾值，從0.1到0.9，步長為0.1。然后根據(jù)閾值將2.2節(jié)構(gòu)建的動態(tài)特征矩陣G二值化，計算每種閾值處理后人腦動態(tài)特征矩陣G的稀疏度，即計算矩陣G中元素1與元素0的比值。隨著閾值增大，人腦動態(tài)特征矩陣G的稀疏度的變化趨勢如下圖5所示：

表2 稀疏度為時的閾值

由圖5以及表2可以看出，閾值設(shè)定為0.32時，Kendall相關(guān)分析構(gòu)建的人腦動態(tài)特征矩陣的稀疏度在50%左右，閾值分別設(shè)定為0.46和0.48時，對應(yīng)的Spearman相關(guān)分析和Pearson相關(guān)分析構(gòu)建的人腦動態(tài)特征矩陣的稀疏度也在50%左右。所以將MIC的閾值確定為0.55。此時腦網(wǎng)絡(luò)的稀疏度要求可以得到滿足。

由表2可以看出人腦動態(tài)特征矩陣G的稀疏度在50%時，Kendall相關(guān)分析方法對應(yīng)的閾值較為接近0.32，Spearman和Pearson相關(guān)分析方法對應(yīng)的閾值分別較為接近0.46和0.48，而當MIC構(gòu)建的人腦動態(tài)特征矩陣二值化后稀疏度為50%時，閾值在0.55左右，為了得到與其他3種相關(guān)分析方法大致相同的腦功能網(wǎng)絡(luò)規(guī)模，本文將MIC的閾值確定為0.55。

4 腦網(wǎng)絡(luò)狀態(tài)聚類方法及結(jié)果分析

根據(jù)上文方法所得到的腦網(wǎng)絡(luò)在每個時刻的狀態(tài)是一個4005維的矩陣，由于其維數(shù)過高，無法直接觀測其特征，這使得深入展開人腦網(wǎng)絡(luò)動態(tài)特征的識別、分析和研究變得極其困難。針對這一現(xiàn)狀，需要對動態(tài)特征矩陣進行降維處理。文獻[32]的研究表明，通過利用多種降維算法對腦網(wǎng)絡(luò)人腦動態(tài)特征矩陣進行降維，發(fā)現(xiàn)相對于其他降維算法，t分布隨機領(lǐng)域嵌入算法(t-Distributed Stochastic Neighbor Embedding, t-SNE)降維效果較為合理，在2維空間上能夠觀測到狀態(tài)隨時間的變化特征。所以本論文使用t-SNE降維算法對人腦動態(tài)特征矩陣從4005維降至2維，并進一步研究和分析不同的相關(guān)分析方法對狀態(tài)識別的影響。

由上文研究可知，在對人腦動態(tài)特征矩陣G進行降維之前，還需要對該矩陣進行閾值處理。為了研究靜息態(tài)腦網(wǎng)絡(luò)狀態(tài)的動態(tài)變化，在對含有動態(tài)信息的人腦動態(tài)特征矩陣G進行閾值處理時，并不是簡單地對人腦動態(tài)特征矩陣進行二值化。上節(jié)主要從靜態(tài)的角度構(gòu)建人腦功能網(wǎng)絡(luò)，不需要考慮rs-fMRI數(shù)據(jù)的動態(tài)變化，所以將相關(guān)系數(shù)矩陣進行閾值二值化時，只需考慮腦區(qū)之間是否存在功能連接，即相關(guān)系數(shù)矩陣二值化后矩陣中元素1表示在某兩個腦區(qū)之間有功能連接，元素0表示兩個腦區(qū)之間不存在功能連接。為了觀測靜息態(tài)腦功能網(wǎng)絡(luò)的動態(tài)變化，本文將人腦動態(tài)特征矩陣G中正值元素大于閾值T的置為1，負值元素絕對值大于閾值T的置為-1，將元素絕對值小于閾值T的元素置為0。因為元素-1可以表示在一個窗口時間內(nèi)，兩個腦區(qū)之間存在功能連接并有著反向響應(yīng)趨勢，所以應(yīng)該將動態(tài)特征矩陣G中相關(guān)系數(shù)值的正負考慮進去。

動態(tài)特征矩陣G閾值處理的方法如式(8)：

(8)

其中i=1,2,…,177，j=1,2,…,4005。由上文分析可知，Spearman相關(guān)分析的閾值為0.46，Kendall相關(guān)分析的閾值為0.32，Pearson相關(guān)分析的閾值為0.48，MIC相關(guān)分析的閾值為0.55。6個樣本用2.2節(jié)的方法以窗口大小n=20，步長l=1構(gòu)建人腦動態(tài)特征矩陣并進行以上閾值分析后，使用t-SNE降維后的可視化結(jié)果圖如下圖6所示。

通過圖6可以看出，4種相關(guān)系數(shù)構(gòu)建的人腦動態(tài)特征矩陣經(jīng)過降維后都能表達靜息態(tài)人腦的狀態(tài)變化過程。上圖中每一個點表示的是人腦動態(tài)特征矩陣G通過降維投影到二維空間內(nèi)后形成的數(shù)據(jù)點，所以每個圖中共有177個數(shù)據(jù)點，數(shù)據(jù)點的顏色代表著窗口次序，如顏色條所示，從藍色到紅色依次顯示的是窗口狀態(tài)從1到177隨時間的狀態(tài)表達。由于t-SNE降維方法是高維空間內(nèi)樣本對和低維空間內(nèi)樣本對相似概率的度量和映射，所以映射到2維空間中表現(xiàn)的并不是腦網(wǎng)絡(luò)狀態(tài)的實際分布，而是樣本間的距離關(guān)系，從圖中可以看到，盡管每次降維實驗高維狀態(tài)投映到二維坐標上的絕對位置可能不一樣，但點之間的相對位置反映的才是狀態(tài)之間的實際關(guān)系。

本文研究主要觀測不同狀態(tài)隨時間的變化趨勢，以及狀態(tài)在哪些點發(fā)生了跳變，降維后，無論哪種分析方法均能將不同時間點的狀態(tài)分為4到6類，這表明了人腦在靜息態(tài)時腦區(qū)之間的交互隨時間是有狀態(tài)變化的，并且從斷點分析可知，多數(shù)圖像中斷點的位置是在相近范圍內(nèi)的，如果考慮人腦樣本在動態(tài)腦網(wǎng)絡(luò)構(gòu)建中的個體差異性，那么雖然同一種樣本使用不同相關(guān)分析方法構(gòu)建的人腦動態(tài)特征矩陣經(jīng)過降維后結(jié)果有差異，但基本的聚簇和斷點范圍是相似的。這可以說明在fMRI數(shù)據(jù)采集過程中，人腦狀態(tài)的變化是具有一定的規(guī)律性的。

圖6 t-SNE降維可視化結(jié)果Fig.6 t-SNE reduced dimensional visualization results

5 結(jié)論

本文以辨別和描述人腦網(wǎng)絡(luò)的狀態(tài)特征為目標，針對靜息態(tài)腦功能磁共振的血氧依賴水平信號中的非線性特征，提出使用非線性相關(guān)分析的方法來構(gòu)建動態(tài)特征功能連接矩陣，重點對構(gòu)建腦區(qū)構(gòu)建全腦狀態(tài)矩陣過程中的不同非線性相關(guān)分析方法及關(guān)鍵技術(shù)展開了研究，并且討論了如何通過閾值控制來排除動態(tài)特征矩陣中相關(guān)程度不夠強的連接，從而確保所構(gòu)建腦網(wǎng)絡(luò)的合理性。

研究發(fā)現(xiàn)，無論采用Spearman、Kendall和MIC三種非線性分析方法中的哪一種，在確保閾值選擇合理的條件下，采用基于滑動窗口的動態(tài)數(shù)據(jù)分析技術(shù)均可完成對網(wǎng)絡(luò)屬性相似的全腦動態(tài)功能網(wǎng)絡(luò)構(gòu)建。本文還實現(xiàn)了對所構(gòu)建的全腦功能網(wǎng)絡(luò)的降維和可視化聚類分析，實驗結(jié)果發(fā)現(xiàn)人腦在靜息態(tài)的情況下，腦網(wǎng)絡(luò)存在自發(fā)的狀態(tài)變化，對于不同的相關(guān)分析方法，動態(tài)特征矩陣經(jīng)過降維后，可視化結(jié)果顯示出的人腦在靜息態(tài)時人腦功能網(wǎng)絡(luò)狀態(tài)的變化具有一定的規(guī)律性，從而驗證了該方法的有效性。

盡管由于在腦網(wǎng)絡(luò)構(gòu)建過程中存在個體差異性，但在不同樣本下相同分析方法在聚類和類簇邊界上相似，而對于相同樣本，不同的相關(guān)分析方法構(gòu)建的動態(tài)特征矩陣可視化降維聚類結(jié)果也只存在很小的差異，上述實驗結(jié)果證明了本文方法的普適性和適用性。

下一步將使用本文方法對腦疾病樣本展開研究，以進一步驗證其有效性，同時對窗口參數(shù)對實驗結(jié)果的影響和狀態(tài)演變規(guī)律展開深入研究。