楊艷麗，李 瑤，谷金曄，李欣蕓，陳俊杰

(太原理工大學 a.信息與計算機學院，b.軟件學院，太原 030024)

利用功能性磁共振成像技術，從復雜網絡分析角度研究腦疾病的發(fā)病機理，已成為腦科學領域的研究熱點之一[1]。該方法主要包括：節(jié)點的定義，連接的定義以及網絡的分析。節(jié)點定義最常用的方法是利用腦圖譜劃分節(jié)點；然而這種節(jié)點的定義方法過于依賴先驗的腦模板，不同模板建立起來的網絡可能會存在顯著差異，這將嚴重影響功能連接網絡的分析結果[2]。ALLEN在研究中使用組獨立成分分析(group independent component analysis，group ICA)提取到的獨立成分作為功能的網絡節(jié)點[3]。

傳統(tǒng)的靜息態(tài)功能連接假定功能連接在時間上是靜止的，功能連接網絡在掃描時間內穩(wěn)定不變。這種方法的主要問題是忽略了掃描時間內可能發(fā)生的神經活動或相互作用[4]。然而，最近的研究已經表明大腦的功能連接承載著豐富的時間信息，無論是任務態(tài)還是靜息態(tài)，功能連接都會隨著神經活動的變化而變化[5]。CHEN使用滑動窗口的方法來劃分整個RS-fMRI時間序列，從而構建出了高階功能連接網絡[6]。這種動態(tài)功能連接方法，在傳統(tǒng)的功能連接網絡上添加了時變因素，充分考慮了時間特性。

目前腦網絡的分析方法有很多，其中最常用的是使用網絡拓撲指標來描述腦網絡特征[7]。然而這種方法受網絡規(guī)模和稀疏度的影響較大，且容易忽視網絡中一些有用的拓撲結構。DU先對腦網絡進行閾值化處理，然后在無權圖上使用頻繁子圖挖掘和判別性特征選擇方法[8]。但閾值化處理時，無論是單閾值還是多閾值，都會造成網絡信息的丟失。KONG et al通過構建腦網絡的加權圖(邊以一定概率存在)，進行加權圖的頻繁子圖挖掘，并將無權圖上的判別性特征選擇方法推廣到了加權圖上[9]。

1 實驗材料及方法

1.1 數據采集及預處理

本文的所有數據都來自ADNI(Alzheimer's disease neuroimaging initiative)數據庫。ADNI數據庫自建立以來，為阿爾茲海默癥(Alzheimer's disease，AD)相關疾病及早期診斷方法的研究提供了豐富的數據支持[11]，得到了廣泛的使用。本實驗數據包括30名正常被試和29名阿爾茲海默癥患者，被試基本信息如表1所示。表中數值是平均值±標準差；NC表示正常對照組；NMSE表示簡易精神狀況檢查量表；a為雙樣本T檢驗；b為皮爾遜卡方檢驗。

表1 被試基本信息Table 1 Basic information for participants

原始數據預處理使用SPM8(http://www.fil.ion.ucl.ac.uk/spm)進行。主要步驟包括1) 時間片校正：校正圖像在獲取時間上的差異；2) 頭動校正：刪除頭動大于3 mm或轉動大于3度的被試數據，丟棄的被試數據不包含在最終的被試中；3) 空間標準化：對圖像進行12維度的優(yōu)化仿射變換，將其標準化到3 mm體素的MNI(montreal neurological institute)標準空間中；4) 平滑：消除不同被試腦結構之間的細微差別、提高信噪比，在10 mm半高全寬高斯空間上進行平滑。

1.2 方法框架

基于組獨立成分分析的加權高階功能連接網絡的靜息態(tài)功能性磁共振成像數據分類方法包括以下4個步驟。

1) 數據采集及預處理。

2) 組獨立成分分析。

組獨立成分分析的主要步驟包括：數據降維、空間獨立成分估計、數據反重構、篩選腦網絡成分[12]。

3) 高階功能連接網絡的構建,主要包括兩個過程:

a.選擇長度固定的滑動窗口，對獨立成分的時間序列進行時間窗劃分，并計算每個時間窗下的皮爾遜相關系數矩陣，得到低階功能連接網絡；

b.堆棧所有的低階功能連接網絡，再在其上計算皮爾遜相關系數矩陣，得到高階功能連接網絡；

4) 加權圖的頻繁子圖挖掘及特征選擇,主要包括兩個過程：

a.對構建的全連接加權的腦網絡進行頻繁子圖挖掘，得到加權的頻繁子圖模式；

b.計算基于動態(tài)規(guī)劃的判別分數函數的相關統(tǒng)計指標。

1.3 組獨立成分分析

組獨立組分分析使用GIFT(http://mialab.mrn.org/software/gift)工具箱進行。本文采用信息極大化(Infomax)算法進行獨立成分的估計和分離，該算法的核心思想是通過最大化輸入和輸出的互信息等價地達到輸出各分量間互信息的最小化[13]。首先使用GIFT軟件估計獨立成分的個數，并結合先驗的研究信息[14]，將獨立成分的個數設為54；接著標準化獨立成分以消除被試個體差異的影響；隨后為了確保獨立成分的穩(wěn)定性和可靠性，在ICASSO軟件上重復了20次Infomax算法；最后使用數據反重構方法獲取被試獨立成分的空間分布和時間序列。

使用組獨立成分分析提取的獨立成分既包括本文感興趣的腦網絡成分，也包括其他無關成分或含噪聲較多的成分。因此需要使用先驗模板匹配方法對這些獨立成分進行篩選，并結合人工肉眼檢查方法對這些成分進行進一步的確認。腦網絡成分的篩選標準包括：與先驗模板匹配的多重回歸系數較大；主要激活區(qū)域分布在灰質區(qū)域。最終去除了32個無關成分或含噪聲較多的成分后保留了22個腦網絡成分，并確認了這些獨立成分分別屬于聽覺網絡、感覺運動網絡、視覺網絡、默認網絡、注意網絡和額葉網絡。

1.4 高階功能連接網絡

在構建高階功能連接網絡之前，還需要對被試獨立成分的時間序列進行后處理，包括去線性漂移和低通濾波，以降低低頻漂移和高頻生理噪聲的影響。接著在每個獨立成分的時間序列上，選擇長度固定的滑動窗口，按照一定步長進行時間窗滑動。假設時間序列的長度為K，滑動窗口的長度為w，滑動窗口每次移動的步長為s，則可以將整個時間序列(維數為K×N，N為獨立成分的個數)劃分為l個相互重疊的時間窗口(維數為w×H).其中l(wèi)的計算公式如下：

l=[(K-w)/s+1] .

(1)

接著在每個時間窗下，計算獨立成分兩兩之間的皮爾遜(Person)相關系數，就可以得到l個低階功能連接網絡(維數為N×N)，它反映了所有獨立成分之間連接強度隨時間的變化情況。其中皮爾遜相關系數的計算公式如下：

(2)

式中：ri,j表示獨立成分i與獨立成分j時間序列的皮爾遜相關系數，cov(i,j)表示兩個獨立成分時間序列的協(xié)方差，σi和σj分別表示兩個獨立成分i與j時間序列的標準差。

最后在新生成的相關時間序列上，再次計算獨立成分兩兩之間的皮爾遜相關系數，就得到了高階功能連接網絡，維度為

(3)

所謂高階體現在矩陣中的任一元素都能反映了一條功能連接對另一條功能連接的影響，體現了多達4個獨立成分之間的交互作用。

1.5 加權圖的頻繁子圖挖掘

在高階功能連接網絡中，功能連接之間的皮爾遜相關系數可以視為功能連接之間交互的概率。因此依據概率的定義，需要將負值的皮爾遜相關系數按零值處理。本文在加權的高階功能連接網絡上進行頻繁子圖挖掘，頻繁子圖挖掘的目的就是在整個網絡中找出出現概率較高的連接模式(即子圖)[15]。

加權圖的頻繁子圖挖掘算法如下：

1) 遍歷加權的高階功能腦網絡，得到所有的1-子圖模式，按照公式(3)計算期望支持度。

2) 如果1-子圖模式的期望支持度大于給定支持度，則把它加入到1-子圖模式集合中。

3) 對每個1-子圖模式，調用模式增長算法，得到新的子圖模式集合，繼續(xù)使用公式(3)計算期望支持度。同樣地，如果新子圖模式的期望支持度大于給定支持度，就把它加入到k-子圖模式集合中(k為新子圖模式的邊數)。

4) 重復上步，直到沒有滿足條件的子圖生成為止。

模式增長算法如下：

1) 將(k-1)-子圖模式所在子搜索空間標記為i.

2) 在1-子圖模式集合中找出標號>i的子圖，將其添加到(k-1)-子圖上，得到k-子圖。如果k-子圖模式連通，則將其加入到k-子圖模式集合中。

給定一個圖集G，一個子圖g的期望支持度定義為：

(3)

式中：Exp(g)表示子圖模式g的期望支持度，n表示加權圖的個數；Pr(g?Gj)表示子圖模式g在加權圖Gj上的概率；k表示子圖模式g的邊數；pj(ek)表示子圖模式g的邊ek在加權圖Gj上的概率。

上述加權圖的頻繁子圖挖掘中挖掘出的子圖有很多，需要從中選出具有判別性的子圖作為特征，本文采用的方法是KONG et al[9]提出的基于統(tǒng)計指標的方法。但這種方法受統(tǒng)計指標的影響較大且最終的分類正確率較低，因此我們提出了一些新的統(tǒng)計指標，并結合無權圖上的判別分數函數(見表2)，通過分類正確率來衡量子圖模式的判別性。

表2 判別分數函數Table 2 Discriminant fraction function

表3 統(tǒng)計指標公式Table 3 Statistical indicator formula

表3中前3種：期望(mean)、方差(variance)和φ-probability，是KONG et al[9]在實驗中使用的方法。受該方法的啟發(fā)，本文提出了兩種新的統(tǒng)計指標——修正的樣本均值(trimmean)和任意階中心距(moment)，來驗證這種方法的合理性。另外，表3中只給出了正樣本中的計算方法，對負樣本類似。

2 實驗結果與分析

2.1 獨立成分分析結果

從圖1中可以看出，注意網絡中包含的成分最多，其次是感覺運動網絡和視覺網絡，然后是默認網絡，額葉網絡和聽覺網絡中包含的成分最少。這說明注意網絡的分布較廣，而聽覺網絡的分布較窄，這與以往基于腦圖譜的全腦功能連接的研究結果基本一致。

圖1 靜息態(tài)腦網絡成分的峰值區(qū)域Fig.1 Peak area of static brain network

2.2 頻繁子圖挖掘結果

圖2中顯示的是作為特征進行分類的20個頻繁子圖模式，不同顏色的節(jié)點分別代表不同的靜息態(tài)腦網路。從中可以發(fā)現，感覺運動網絡、默認網絡和視覺網絡中連接的邊比較密集，這說明這些網絡在患者和正常被試的分類過程中起到的作用比較大。

圖2 頻繁子圖模式Fig.2 Frequent subgraph mode

2.3 分類結果與分析

由于被試人數有限，本文使用基于RBF核函數的支持向量機(support vector machine,SVM)進行分類，并使用留一交叉驗證(leave-one-out cross validation,LOOCV)方法來評估分類器的的泛化性能(具體實驗使用基于MATLAB的LIBSVM工具包進行)。分類器的性能可以量化為正確率(accuracy)，靈敏度(sensitivity)和特異度(specificity).其中，靈敏度可以衡量該方法正確識別患者的能力，特異度衡量的是正確識別正常被試的能力。進一步使用ROC曲線下的面積(the area under receiver operating characteristic curve，AUC)來評價不同方法的性能表現，具體結果如表4所示。

從表4中可以看出，本文提出的兩種統(tǒng)計指標(修正的樣本均值和任意階中心距)，都可以很好地找到與分類相關的特征，獲得較高的分類正確率；但與φ-probability方法相比，在分類正確率、靈敏度、特異度以及AUC值上均有些差距。雖然在Confidence判別函數上這兩種方法的分類效果比φ-prob方法好，但分類效果最好的還是φ-prob-HSIC方法，其分類正確率達到95.24%。從判別分數函數角度來看，雖然最好的分類效果出現在HSIC判別函數上，但整體來看，還是Gtest判別函數的分類效果較為穩(wěn)定，分類正確率都在90%以上。因此，使用φ-prob統(tǒng)計指標和Gtest判別函數的分類效果較好。

表4 不同特征選擇的分類結果比較Table 4 Comparison of classification results of different feature selection

此外，為了驗證基于獨立成分分析的高階功能連接網絡的優(yōu)越性，本實驗使用相同的被試數據分別構建了基于AAL的低階功能連接網絡、基于AAL的高階功能連接網絡、基于group ICA的低階功能連接網絡和基于group ICA的高階功能連接網絡，并使用加權圖的頻繁子圖挖掘和φ-prob-Gtest特征選擇方法，進行SVM分類，結果如表5所示。

表5 不同方法的分類結果比較Table 5 Comparison of classification results of different methods

從表5中可以看出，使用本文方法時，在group ICA提取的獨立成分的時間序列上構建高階功能連接網絡，并以加權圖的頻繁子圖挖掘出的子圖作為特征進行分類識別，分類正確率達到95.24%，遠高于其他方法。而且高階功能連接腦網絡上的分類效果要遠好于傳統(tǒng)的低階功能連接網絡，這表明高階功能連接腦網絡上包含更多的對分類識別有用的信息。另外本文方法在靈敏度和特異度等指標上也較傳統(tǒng)方法有一定的提升，在一定程度上能夠為臨床診斷提供客觀參考。

2.4 不同參數的影響

本文所使用的高階功能連接網絡會受一些參數的影響，如滑動窗口的長度w和每次移動的步長s等。由公式(1)可知，選擇不同參數時，得到的時間窗數l不同，進而構建出來的高階功能連接網絡也不同。因此，本實驗先在滑動窗口每次移動步長都為1 s的情況下，分別選擇滑動窗口的長度為23，43，63，83，103 s，探究滑動窗口的長度對分類效果的影響，結果如圖3(a)所示。當滑動窗口的長度為63 s時，分類效果最好，且分類曲線呈開口向下的拋物線形狀。這可能是因為：當選擇的滑動窗口的長度過小時，劃分出來的時間窗數很多，其中包含的特征數更多，從中選擇判別性特征的難度將增大；而選擇的滑動窗口的長度過大時，劃分出來的時間窗數很少，其中所包含的時變特征更少，這會大大影響了分類效果。此外，本實驗還在滑動窗口的長度都為63 s時，分別選擇滑動窗口每次移動的步長為1 s，

圖3 不同參數的影響Fig.3 Impact of different parameters

2 s，3 s，4 s，5 s，探究滑動窗口每次移動的步長對分類結果的影響，結果如圖3(b)所示。當滑動窗口每

次移動的步長為1 s時，分類效果最好，且每次移動的步長越長，分類的效果越差。這可能是因為：滑動窗口每次移動的步長越長，得到的時間窗數越少，其中所包含的時變特性就更少，分類正確率也會隨之降低。因此本實驗中選擇的滑動窗口長度為63 s和每次移動步長為1 s.

3 總結

本文在對靜息態(tài)功能性磁共振成像數據分類時，先使用組獨立成分分析方法提取靜息態(tài)腦網絡成分，接著構建高階功能連接網絡，然后使用加權圖的頻繁子圖挖掘和判別性特征選擇方法，最后進行SVM分類。結果表明，該方法具有較好的分類效果，可以為臨床醫(yī)學診斷提供參考依據。然而，該方法構建高階功能連接網絡的過程較為復雜，大大增加了運算量；而且該方法在加權圖上做頻繁子圖挖掘時存在大量的冗余操作，效率較低，只適用于小樣本數據集。因此，如何降低高階功能連接網絡的復雜性和加權圖的頻繁子圖挖掘的冗余度，是今后研究的重點。