基于擴展的低階多元廣義線性模型的腦節點識別方法

楊雅倩，唐紹婷

(1.數學、信息與行為教育部重點實驗室(北京航空航天大學)，北京 100191；2.北京航空航天大學 大數據科學與腦機智能高精尖創新中心，北京 100191；3.北京航空航天大學 數學與系統科學學院，北京 100191)(*通信作者電子郵箱tangshaoting@buaa.edu.cn)

0 引言

功能性磁共振成像(functional Magnetic Resonance Imaging, fMRI)技術是一種通過檢測血管中血氧水平變化來測量大腦活動的技術，它因其無創性和高空間分辨率而被廣泛運用于大腦活動的研究當中，現已有多篇相關文章，其中一個重要研究方向是比較不同條件下的大腦反應。具體而言，在許多心理fMRI實驗中，研究人員對兩次實驗中的大腦反應進行比較，旨在找出不同環境下對某一刺激具有不同反應的大腦區域。

一種常見的fMRI實驗數據分析方法是廣義線性模型(General Linear Model, GLM)[1]，它通過血液動力學響應函數(Hemodynamic Response Function, HRF)反映不同大腦節點對不同類型刺激的響應情況。在此基礎上又發展了參數模型、非參模型、半參模型等多種HRF估測方法，其中參數方法主要利用有限模型參數來刻畫不同HRF之間的差異特征，具體包括泊松模型[2]和經典血液動力學響應函數模型(Canonical HRF, Canonical)[3]等；非參數方法通常將HRF表示為一系列基函數線性組合的形式，因此在衡量不同大腦區域HRF函數特征時更加靈活，其主要包括平滑有限脈沖響應(Smooth Finite Impulse Response, SFIR)[4]和正則化和廣義交叉驗證(Tikhonov-regularization and Generalized-Cross-Validation, Tik-GCV)[5]；此外，基于大腦活動的群體一般性和個體特異性，Zhang等[6-7]建立了更為靈活的半參數模型。

上述基于GLM框架的HRF計算模型都是單節點方法，即一次只對一個大腦節點的fMRI時間序列進行分析。由于空間上相鄰的節點往往具有相似的fMRI數據，因此將大腦的空間信息并入HRF計算中將提升模型效果。基于同一區塊大腦節點享有共同的函數形狀而僅僅在振幅上有所不同這一假設，Vincent等[8]在大腦區塊貝葉斯模型[9-10]的基礎上提出HRF振幅空間先驗；Chaari等[11]進一步提出了同一大腦區塊中隨空間變化的HRF振幅先驗，并開發了聯合分割檢測估計程序。同時，fMRI研究向多個體、多維度數據發展，例如Degras等[12]提出了多個體功能磁共振成像HRF估計的貝葉斯模型，Zhang等[13-14]建立了適用于fMRI時間序列中復雜時間空間相關性的貝葉斯方法。

現有的基于GLM框架的模型在進行大腦反應比較實驗時，通常需要提取HRF的某些低維特征(如高度等)，利用假設檢驗對這些低維特征進行比較，進而獲取兩次實驗中反應不同的大腦節點。然而這種方法只能比較HRF某一特征(如高度)的差異而忽略了其他特征(如函數形狀)的差異，導致估測結果具有很大的不確定性。而假設檢驗本質上是對每一個節點進行分析，因此忽略了fMRI數據的空間特性，導致被識別節點的準確度較低。針對這一缺陷，Zhang等[15]提出了綜合所有腦節點空間信息的低階多元廣義線性模型(Low-Rank Multivariate General Linear Model, LRMGLM)，該模型利用空間矩陣靈活描述了HRF的變化，識別節點準確率更高，但fMRI數據的低信噪比和高變異性導致其計算復雜、計算效率低下。此外，由于LRMGLM只能對同一實驗的大腦反應進行比較，因此在比較不同實驗的大腦反應時，需要額外設置一組相同的刺激作為參考項，適用范圍有限。

對此，本文提出一種基于擴展的LRMGLM(Extended LRMGLM, ELRMGLM)的腦節點識別方法。該方法建立了可同時處理兩次實驗所有節點數據的ELRMGLM模型，該模型通過將血液動力學響應函數(HRF)的特征矩陣轉化為兩個低階矩陣相乘的形式，在實現腦節點靈活比較的同時簡化參數運算。模型參數利用基于fMRI數據時空特性的優化函數和迭代算法進行估測，同時為降低模型對fMRI數據高變異性和低信噪比的敏感度，開發了基于K-means的快速選擇策略來實現兩次實驗中反應不同大腦區域的快速選擇。

1 本文模型建立

(1)

其中:

以往基于GLM框架的方法一次只對一個大腦節點或一次實驗數據進行單獨處理，本文提出一個可以同時處理兩次實驗所有節點數據的聯合模型，以更多的時空信息來減少fMRI數據中噪聲的干擾。由于不同個體大腦節點的HRF形狀不同，首先利用B-樣條插值對其進行擬合：

(2)

(3)

(4)

相較于LRMGLM，ELRMGLM可以直接處理兩次實驗所有節點的數據信息，不僅比原模型適用范圍更廣，而且可以利用更多的時間空間信息來提高識別的準確度。相較于傳統的單節點模型，ELRMGLM可以通過保持相應Uk(q)不變，對Vk(q)進行比較的方法來實現對HRF靈活而全面的比較。例如，在計算過程中可以令U2(1)=U2(2)，比較V2(1)和V2(2)的估計值來識別兩次實驗中對第二種刺激有不同反應的大腦區域。

2 參數估計

2.1 優化函數

在利用上述模型對數據進行分析時，通常希望估計得到的Y值盡可能地接近真實值，令

q=1,2}

則有如下代價函數：

SSE(Θ)=

由于HRF在時間上連續，為避免過度擬合，得到矩陣Uk(q)上的時間平滑懲罰項：

(5)

同時由于空間上相鄰的大腦節點通常有相似的fMRI時間序列和HRF函數，因此得到矩陣Vk(q)上的空間平滑懲罰項：

(6)

τPS(Vk(q))]

(7)

2.2 迭代算法

1)給定V，找到U、β，最小化

2)給定U、β，找到V，最小化

上述迭代算法的計算效率主要取決于兩個最優子問題的計算效率，由于步驟1)、2)中的目標函數均為二次函數且存在最優解析表達式，因此可對其進行直接求解，具體推導過程在此不多作贅述。當數據量較大、數據維度過高時，直接求解可能耗時較長，此時可利用最速梯度下降法加速運算。

3 基于K-means的快速選擇策略

懲罰參數λ和τ分別控制著HRF的時間平滑性和空間平滑性。在懲罰優化問題中，普通交叉驗證(Ordinary Cross Validation, OCV)和廣義交叉驗證(Generalized Cross Validation, GCV)[18]是選擇懲罰參數的經典方法；在成像數據分析中，文獻[19-20]提出了基于GCV的選擇過程，同時文獻[21]對約束最大似然法(Restricted Maximum Likelihood, REML)進行了研究。由于本文節點數量較多且有兩個懲罰參數，OCV耗時過長而REML不能直接適用，因此提出了基于K-means的快速選擇策略。

由于本文的研究重點是比較兩次實驗中反應不同的大腦節點而非計算HRF的具體值，即只需要選擇能清楚區分反應相同和反應不同腦節點的懲罰參數組合，因此可利用聚類方法加快懲罰參數和大腦節點的選擇過程，將輪廓系數作為選擇標準，在提高模型對懲罰參數容忍度的同時保證所選節點的準確性，具體過程如下：

1)對每個懲罰參數，在e-1～e5范圍內選取大量候選參數值，利用2.2節中提到的迭代算法對帶有不同懲罰參數組合的ELRMGLM進行模型參數估計，得到相應的時間矩陣Uk(q)和空間矩陣Vk(q)。

2)計算每組懲罰參數對應矩陣Vk(q)在兩次實驗中的差值S=|Vk(1)-Vk(2)|，其中k為所要比較的刺激類型。

3)對每組懲罰參數，利用K-means聚類將所有節點的S值分為兩類，選擇輪廓系數最大的懲罰參數組合和聚類結果，其中S均值較大的群組為兩次實驗中對刺激k反應不同的大腦節點集合。

上述快速選擇策略的原理在于兩方面：首先，大部分懲罰參數組合對應的大腦節點S值可以被自然地分為兩部分，一部分在一個較小值附近上下浮動，另一部分的S值則顯著較大，它們分別代表了由生理噪聲引起的低頻漂移和兩次實驗中不同的大腦反應，因此可通過聚類方法加以區分；其次，不同的懲罰參數組合影響了兩組節點S值的差異性以及各組節點S值的穩定性，從而影響了聚類結果的準確性，而輪廓系數是一種評價聚類效果好壞的方式，它衡量了個體相較其他群集與其所屬群集的相似程度，其范圍從-1到1，輪廓系數值越大表明聚類效果越好，因此可將輪廓系數作為懲罰參數的選擇標準。一組輪廓系數較高的懲罰參數組合會使大腦節點的S值具有高差異性和高穩定性，聚類效果較好，從而保證了被選節點的準確性和可靠性。

聚類方法提高了算法對不同懲罰參數組合的容忍度，高輪廓系數保證了所選節點的準確性，通過采用基于K-means的快速選擇策略，利用聚類加快懲罰參數和大腦節點的選擇過程，可以在保證準確性的同時快速找到兩次實驗中反應不同的大腦區域。

4 實驗分析

4.1 實驗數據

本文采用與文獻[22-24]中真實fMRI實驗相同的實驗設計進行分析，該實驗使用了以頻率27.9%、21.1%、50.8%和0.2%隨機出現的4種不同刺激，共有106個受試者，每個受試者的fMRI時間序列包含205次掃描，每次掃描時間為2 s。考慮到LRMGLM只適用于存在參考項的實驗數據，本文采用文獻[15]中使用的實驗數據集。該數據集共含有三組fMRI數據，分布在15×15×15的大腦網格上。在組1中，HRF函數服從經典形式且前兩種刺激的HRF在振幅和延遲上均不同；在組2中，HRF函數形狀改變，其余與組1保持相同；在組3中，改變第二種刺激的振幅使得前兩種刺激的HRF僅僅在延遲時間上有所不同，其余與組2保持相同。此外，該數據集中與第二種刺激有關的HRF參數只針對中心9×9×9的節點網格，其余大腦節點的HRF參數與第一種刺激相同。為比較實際問題中大腦節點在兩次實驗的中不同反應，本文增設含三組相同fMRI數據的實驗二，并令其中與第二種刺激有關的HRF參數只針對中心偏右的5×9×9大腦網格，其余與實驗一保持相同。此時約9.6%的大腦節點在兩次實驗中對第二個刺激反應不同。

值得注意的是，本文采用的實驗設置并不完全遵循提出的ELRMGLM，盡管如此，實驗結果將表明ELRMGLM方法能夠超越更適合該實驗數據的單節點分析方法和LRMGLM。

4.2 實驗結果與分析

本文將提出的ELRMGLM運用于實驗數據，為保持運算簡潔，令P=2來刻畫不同個體大腦節點HRF在振幅和延遲上的差異。圖1展示了兩次實驗中對第二種刺激反應不同的大腦區域，其中深色區域表示被選中節點，白色區域表示未選節點，X、Y、Z分別表示大腦節點的三維坐標。可以看到，ELRMGLM方法選中了大腦中心偏左的4×9×9網格，準確度較高。

圖1 被選中節點圖像Fig. 1 Image of selected voxels

圖2 三組實驗數據的ROC曲線Fig. 2 ROC curves of three experimental datasets

下面將ELRMGLM分別與canonical方法[3]、SFIR方法[4]和Tik-GCV方法[5]進行比較。在進行計算時，先用上述單節點方法對HRF進行估計，再提取得到的HRF低維特征(如高度)，通過t假設檢驗對提取的低維特征進行比較，進而獲得兩次實驗中反應不同的大腦區域。通過改變t檢驗中的P值，得到了顯示不同臨界值對應的真正類率(True Positive Rate, TPR)和負正類率(False Positive Rate, FPR)組合的受試者工作特征曲線(Receiver Operator characteristic Curve, ROC)。為進行比較，本文還通過使用不同的懲罰參數組合畫出了ELRMGLM的ROC曲線，具體如圖2所示。

從圖2可知：Canonical方法在第一組和第三組實驗數據中表現較差，但在第二組實驗數據中表現較好，其TPR和FPR分別達到了約80%和4%。與之相反，SFIR在第一組和第三組數據中表現較好，其TPR和FPR分別達到了約95%、3%和99%、20%，但在第二組數據中表現一般。Tik-GCV在所有單節點分析方法中表現最為穩定，其TPR和FPR在三組數據中均達到99%和20%左右。通過使用第3章基于K-means的快速選擇策略，本文提出的ELRMGLM在三組實驗數據中均表現優異，其TPR和FPR均實現了99%以上和1%以下 (其TPR和FPR分別達到99.73%、0.44%，99.99%、0.09%和99.99%、0.01%)，比以上三種方法的最優結果分別提升了約20%、8%、20%，不僅實現了高敏感度(sensitivity)和高特異度(specificity)，同時在不同的數據集上表現穩定。

下面以第一種刺激作為參考項，對ELRMGLM和LRMGLM方法進行比較，其結果如表1所示。可以看到在三組實驗中，ELRMGLM對兩次實驗中反應不同的大腦節點識別準確度略高于LRMGLM。在計算時間上，ELRMGLM的平均迭代次數少于100而LRMGLM需要迭代上萬次；同時ELRMGLM的單次迭代時間為8 s，遠小于LRMGLM的單次迭代時間60 s。不難算出，ELRMGLM的計算時間是LRMGLM的1/750，算法效率大幅提高。

表1 ELRMGLM與LRGMLM方法比較Tab. 1 Comparison of ELRMGLM and LRGMLM

5 結語

本文提出了一種用于識別兩次實驗中反應不同大腦區域的擴展的低階多元廣義線性模型(ELRMGLM)。該模型同時綜合了兩次實驗的數據信息，通過帶懲罰項的優化函數考慮了fMRI數據的時空特性，并利用K-means聚類提高了模型對參數的容忍度，進而實現了對大腦節點的快速準確識別。通過在三組實驗數據集上的分析，該模型在準確度、計算效率和穩定性方面均高于現有模型。ELRMGLM主要用于群體大腦活動的評估，對個體大腦反應的比較還稍有欠缺，因此如何擴展ELRMGLM使其適用于大腦活動的個體差異將成為下一步的研究重點。