葛詩琪，孟祥輝，張浩威，李靚珠

（長安大學(xué) 電子與控制工程學(xué)院，陜西 西安 710064）

0 引 言

大腦不同功能的實現(xiàn)需要多個腦區(qū)的協(xié)同作用，這些腦區(qū)在空間上可能相隔較遠(yuǎn)，但在激活模式上卻存在著高度的相關(guān)性。功能磁共振成像（Functional Magnetic Resonance Imagining, fMRI）是一種基于血氧水平依賴（BOLD）的成像技術(shù)，為研究人腦功能提供了一種重要方法[1]。其中靜息狀態(tài)功能磁共振成像（RS-fMRI）由于能夠反映人腦的自發(fā)神經(jīng)活動，已成為研究熱點[2]。此外，多種腦功能障礙表現(xiàn)為腦區(qū)間功能連接之間的顯著變化。故對腦區(qū)之間的功能連接（Functional Connectivity, FC）進(jìn)行研究，有利于促進(jìn)對大腦活動本質(zhì)的探索[3]。

傳統(tǒng)的基于RS-fMRI的功能網(wǎng)絡(luò)研究通常假設(shè)系統(tǒng)在一段時間之內(nèi)處于平穩(wěn)狀態(tài)。雖然RS-fMRI是受試者在無特定任務(wù)的情況下得到的磁共振掃描，然而事實上大腦活動是一個復(fù)雜的時變過程[4]。Hutchison等人的研究表明，靜息態(tài)功能連接網(wǎng)絡(luò)隨時間變化具有波動性[5]。因此研究腦功能網(wǎng)絡(luò)的動態(tài)特性具有重要意義。目前最常用的動態(tài)功能連接分析方法是滑動窗口法（Sliding Window Analysis）[6]。Sako?lu等人利用滑動窗口和獨立成分分析發(fā)現(xiàn)動態(tài)連接與受試者的臨床指標(biāo)相關(guān)[7]?；瑒哟翱诜椒檠芯看竽X內(nèi)部狀態(tài)變化提供了一種有效的方法，但也存在明顯的局限性，窗口選擇過大無法準(zhǔn)確捕捉功能連接的變化，窗口過小會引入不真實的波動。研究表明，30～60 s的窗口能夠提供較為穩(wěn)健的動態(tài)連接網(wǎng)絡(luò)估計[6]。

本文提出一種基于張量分解的動態(tài)功能連接網(wǎng)絡(luò)分析方法。首先根據(jù)靜息態(tài)fMRI不同維度的信息構(gòu)建張量模型，并根據(jù)不同模態(tài)的特點選擇約束項進(jìn)行優(yōu)化求解。然后利用動態(tài)檢測法對時間改變點進(jìn)行檢測，構(gòu)建動態(tài)網(wǎng)絡(luò)。為了驗證模型的有效性，在仿真數(shù)據(jù)和真實數(shù)據(jù)上將本文模型獲得的實驗結(jié)果與滑動窗口方法獲得的結(jié)果進(jìn)行了比較。仿真實驗結(jié)果表明，本文算法能夠克服滑動窗口方法對于窗口選擇的限制，顯著降低了實驗結(jié)果的相對誤差（最多降低38.4%）。真實數(shù)據(jù)上的實驗結(jié)果表明，靜息狀態(tài)下不同腦區(qū)之間存在功能連接，主要集中在SM、AUD、DMN、VIS和SN之間。

1 基于張量分解的動態(tài)功能連接網(wǎng)絡(luò)

1.1 張量分解模型

張量（tensor）是對高維數(shù)組的總稱[8]。一維數(shù)組稱為向量，二維數(shù)組稱為矩陣，三維及高維數(shù)組稱為高階張量。CP分解是一種常用的張量分解方法，可以將高階張量分解為一系列秩一張量的線性組合。三階張量T∈RI×J×K的CP分解可以表示為：

式中：“°”表示張量的外積；R為張量的秩。根據(jù)張量的標(biāo)準(zhǔn)矩陣表示形式，可通過下列目標(biāo)函數(shù)對公式（1）求解：

式中： ⊙表示 Khatri-Rao積；||·||F表示Frobenius范數(shù)；A=[a1,a2, ...,aR]∈RI×R，B=[b1,b2, ...,bR]∈RJ×R，C=[c1,c2, ...,cR]∈RK×R，為因子矩陣；T(k)表示張量以第k模態(tài)展開的矩陣。

1.2 基于靜息態(tài)fMRI的稀疏張量分解

為了對靜息態(tài)fMRI數(shù)據(jù)中不同模態(tài)的信息進(jìn)行融合分析，將所有個體的fMRI數(shù)據(jù)堆疊為一個三階張量X∈RI×R×K，其中I表示個體數(shù)量，J表示ROI節(jié)點個數(shù)，K表示時間序列長度。

張量分解完成了對高維數(shù)據(jù)的降維，分解得到的3個矢量分別表示個體、ROI、時間在有效成分中所占的權(quán)重。為了使分解的結(jié)果包含更多有意義的信息，需要根據(jù)不同維度信息對數(shù)據(jù)進(jìn)行必要的約束。對A矩陣添加L2,1正則約束項以避免不同個體之間的相互影響；對B矩陣進(jìn)行正交約束以保證不同ROI之間的線性無關(guān)性；對矩陣C添加L1正則項實現(xiàn)稀疏性。添加約束之后的模型如下：

式中：A∈RI×R，air表示個體i在成分r中所占的權(quán)重；B∈RJ×R，bjr表示第j個ROI在成分r中所占的權(quán)重；C∈RK×R，ckr表示第k個時間點在成分r中所占的權(quán)重。λ1和λ3不僅可以控制重建功能性連接網(wǎng)絡(luò)的特征提取，而且還確定了有效成分的稀疏性和大小，分解的過程中通過調(diào)節(jié)懲罰參數(shù)λ2實現(xiàn)正交性。

1.3 優(yōu)化求解

使用交替最小二乘（Alternating Least Squares, ALS）算法對式（3）求解，其主要思想是在每次迭代中固定其中2個因子矩陣來更新另一個因子矩陣，重復(fù)該步驟直至算法收斂?？傻镁仃嘇、B和C的更新公式分別為：

式中：sgn(·)為符號函數(shù)；Σ∈RR×R是對角矩陣，該矩陣中第r個對角元素為1/||ar||2。當(dāng)目標(biāo)函數(shù)滿足|et+1-et|≤10-8時終止迭代，其中et=||Tt-T||2/||T||2。

2 動態(tài)功能連接網(wǎng)絡(luò)

2.1 時間改變點檢測

靜息態(tài)fMRI數(shù)據(jù)隨時間而波動，故實驗的一個重點就是對時間改變點進(jìn)行檢測，并根據(jù)時間改變點對狀態(tài)進(jìn)行劃分。假設(shè)時間權(quán)重矩陣C可以被劃分為τ+1個狀態(tài)，時間改變點分別位于{t1,t2, ...,tτ}，使用公式（7）對連續(xù)2個時間點的距離進(jìn)行度量[9]。

式中：C(t,k)表示時間點t在第k個有效成分中所占的權(quán)重；R表示有效成分的總數(shù)。通過式（7）確定d(t)是否為時間改變點。

式中，T為時間點個數(shù)，由于距離d(t)服從高斯分布，故選擇μ+2σ作為閾值。

2.2 功能連接網(wǎng)絡(luò)構(gòu)建

在得到時間改變點之后，提取各狀態(tài)的主導(dǎo)ROI，計算各狀態(tài)下的動態(tài)功能連接網(wǎng)絡(luò)?；谙∈鑿埩緾P分解的動態(tài)功能連接網(wǎng)絡(luò)構(gòu)建主要分為以下步驟。

（1）張量分解：對原始數(shù)據(jù)構(gòu)成的張量模型根據(jù)式（3）進(jìn)行分解，得到關(guān)于個體、ROI及時間的矩陣A、B、C。

（2）時間改變點檢測：根據(jù)式（8）基于時間權(quán)重矩陣C確定時間改變點{t1,t2, ...,tτ}，并進(jìn)行狀態(tài)劃分。

（3）提取各狀態(tài)的主導(dǎo)ROI：2個時間改變點(ti,ti+1)之間為一段穩(wěn)定的狀態(tài)。通過矩陣C求出每個穩(wěn)定狀態(tài)下的主要成分，提取B矩陣中的對應(yīng)列，并保留其中權(quán)重較大的部分，提取出其中主要的ROI。

（4）構(gòu)建dFC網(wǎng)絡(luò)：根據(jù)原始數(shù)據(jù)的連接網(wǎng)絡(luò)及步驟（3）得到的主導(dǎo)ROI的權(quán)重矩陣，建立每個狀態(tài)下的動態(tài)連接網(wǎng)絡(luò)。

3 實驗及分析

3.1 模擬實驗

為了評估算法的性能，首先使用模擬數(shù)據(jù)集進(jìn)行驗證，并與滑動窗口方法進(jìn)行比較。使用SimTB[10]工具箱進(jìn)行BOLD時間序列信號生成。產(chǎn)生的信號包括每個受試者的60 ROI，每個ROI包含120個時間點，重復(fù)時間為2 s。時間序列被分為4個不同的狀態(tài)，每個狀態(tài)包含8～16個不同的ROI，共100個受試者。用二進(jìn)制表示信號的激活狀態(tài)，1表示激活，0表示未激活[11]。功能區(qū)劃分及激活狀態(tài)如圖1所示。

滑動窗口方法中我們測試了兩種不同的窗口長度（l=30，40）。本文算法的稀疏性參數(shù)λ1、λ2和λ3使用網(wǎng)格搜索法確定，分別為0.9、2和0.1×10-1，R設(shè)置為60，此時擬合率相對較高，可達(dá)到85.2%。使用真實連通矩陣Cg與算法估計所得連通矩陣CE的相對誤差作為評價指標(biāo)：

圖1 四個功能區(qū)和四個功能區(qū)的激活狀態(tài)

表1顯示了滑動窗口法（SW）和本文方法（STensor）在模擬數(shù)據(jù)集上的相對誤差?？梢钥闯觯姆N狀態(tài)下本文算法對連接矩陣的擬合均優(yōu)于滑動窗口方法，相對誤差均明顯降低。在狀態(tài)1中，當(dāng)滑動窗口的窗口長度取40時，相對誤差最多降低了38.4%。

表1 相對誤差

3.2 動態(tài)功能連接網(wǎng)絡(luò)分析

本文使用的數(shù)據(jù)來源于費城神經(jīng)發(fā)育隊列（Philadelphia Neurodevelopmental Cohort, PNC）[12]。實驗中選取其中878名受試者的RS-fMRI數(shù)據(jù)，時間序列為124，ROI數(shù)量為264。實驗中三種稀疏性參數(shù)使用網(wǎng)格搜索法進(jìn)行選擇，λ1，λ2和λ3分別取0.1，6×103和0.1×10-3。有效成分個數(shù)R取100。在完成稀疏張量分解之后，首先使用動態(tài)方法檢測得到時間改變點為{3，40，75，90}。由于初始時機器及人體心跳呼吸等噪音造成的波動較大，故舍棄時間改變點3，于是時間序列被分為四段相對穩(wěn)定的狀態(tài)[0，40]，[41，75]，[76，90]和[91，120]。找出每個狀態(tài)的主要成分及主導(dǎo)ROI，計算動態(tài)功能連接矩陣，這些ROI可以被分為不同的功能區(qū)域。

圖2分別為四種狀態(tài)下每種功能區(qū)域之間的連接數(shù)，連接數(shù)越大表示連接強度越高。從圖2中可以看出靜息狀態(tài)下大腦并不是靜息的，各個腦區(qū)之間存在著許多自發(fā)性的神經(jīng)活動，尤其是視覺網(wǎng)絡(luò)（VIS）、聽覺網(wǎng)絡(luò)（AUD）、默認(rèn)網(wǎng)絡(luò)（DMN）等之間存在著高強度的功能連接，其中DMN的連通性強于其他腦區(qū)，這與文獻(xiàn)[13]提出的DMN自身具有相對較高的連通性的結(jié)論一致。文獻(xiàn)[14]指出DMN和VIS之間有著顯著的腦行為關(guān)聯(lián)。

圖2 dFC網(wǎng)絡(luò)

圖2中4個狀態(tài)的dFC網(wǎng)絡(luò)中，較高強度的功能連接均集中在SM、AUD、DMN、VIS和SN幾個腦區(qū)中，說明大部分功能連接所對應(yīng)的腦區(qū)具有相似性，且在神經(jīng)解剖學(xué)上有重疊部分[2]。

4 結(jié) 語

針對靜息態(tài)fMRI數(shù)據(jù)的特性，本文提出了一種基于稀疏張量CP分解的動態(tài)腦功能網(wǎng)絡(luò)分析方法。首先利用動態(tài)連接檢測法對張量分解的結(jié)果進(jìn)行狀態(tài)劃分，該方法打破了滑動窗口方法對于窗口長度選擇的限制，顯著降低了實驗的相對誤差。然后對不同狀態(tài)下的功能連接進(jìn)行建網(wǎng)分析。實驗結(jié)果表明，靜息狀態(tài)下各個腦區(qū)之間依然存在聯(lián)系，這些腦區(qū)可能在解剖學(xué)上相互分離，但由于自發(fā)性的神經(jīng)活動形成了功能連接網(wǎng)絡(luò)。在靜息狀態(tài)下，功能連接所對應(yīng)的腦區(qū)具有相似性，且大部分集中在SM、AUD、DMN、VIS和SN幾個腦區(qū)當(dāng)中，這一結(jié)論也在其他多個文獻(xiàn)中得到了驗證。