基于稀疏張量分解的動態功能網絡分析

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

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

0 引 言

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

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

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

1 基于張量分解的動態功能連接網絡

1.1 張量分解模型

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

式中：“°”表示張量的外積；R為張量的秩。根據張量的標準矩陣表示形式，可通過下列目標函數對公式（1）求解：

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

1.2 基于靜息態fMRI的稀疏張量分解

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

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

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

1.3 優化求解

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

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

2 動態功能連接網絡

2.1 時間改變點檢測

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

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

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

2.2 功能連接網絡構建

在得到時間改變點之后，提取各狀態的主導ROI，計算各狀態下的動態功能連接網絡?；谙∈鑿埩緾P分解的動態功能連接網絡構建主要分為以下步驟。

（1）張量分解：對原始數據構成的張量模型根據式（3）進行分解，得到關于個體、ROI及時間的矩陣A、B、C。

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

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

（4）構建dFC網絡：根據原始數據的連接網絡及步驟（3）得到的主導ROI的權重矩陣，建立每個狀態下的動態連接網絡。

3 實驗及分析

3.1 模擬實驗

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

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

圖1 四個功能區和四個功能區的激活狀態

表1顯示了滑動窗口法（SW）和本文方法（STensor）在模擬數據集上的相對誤差。可以看出，四種狀態下本文算法對連接矩陣的擬合均優于滑動窗口方法，相對誤差均明顯降低。在狀態1中，當滑動窗口的窗口長度取40時，相對誤差最多降低了38.4%。

表1 相對誤差

3.2 動態功能連接網絡分析

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

圖2分別為四種狀態下每種功能區域之間的連接數，連接數越大表示連接強度越高。從圖2中可以看出靜息狀態下大腦并不是靜息的，各個腦區之間存在著許多自發性的神經活動，尤其是視覺網絡（VIS）、聽覺網絡（AUD）、默認網絡（DMN）等之間存在著高強度的功能連接，其中DMN的連通性強于其他腦區，這與文獻[13]提出的DMN自身具有相對較高的連通性的結論一致。文獻[14]指出DMN和VIS之間有著顯著的腦行為關聯。

圖2 dFC網絡

圖2中4個狀態的dFC網絡中，較高強度的功能連接均集中在SM、AUD、DMN、VIS和SN幾個腦區中，說明大部分功能連接所對應的腦區具有相似性，且在神經解剖學上有重疊部分[2]。

4 結 語

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