楊雅倩，唐紹婷

(1.數(shù)學(xué)、信息與行為教育部重點(diǎn)實(shí)驗(yàn)室(北京航空航天大學(xué))，北京 100191；2.北京航空航天大學(xué) 大數(shù)據(jù)科學(xué)與腦機(jī)智能高精尖創(chuàng)新中心，北京 100191；3.北京航空航天大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，北京 100191)(*通信作者電子郵箱tangshaoting@buaa.edu.cn)

0 引言

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

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

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

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

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

1 本文模型建立

(1)

其中:

以往基于GLM框架的方法一次只對一個(gè)大腦節(jié)點(diǎn)或一次實(shí)驗(yàn)數(shù)據(jù)進(jìn)行單獨(dú)處理，本文提出一個(gè)可以同時(shí)處理兩次實(shí)驗(yàn)所有節(jié)點(diǎn)數(shù)據(jù)的聯(lián)合模型，以更多的時(shí)空信息來減少fMRI數(shù)據(jù)中噪聲的干擾。由于不同個(gè)體大腦節(jié)點(diǎn)的HRF形狀不同，首先利用B-樣條插值對其進(jìn)行擬合：

(2)

(3)

(4)

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

2 參數(shù)估計(jì)

2.1 優(yōu)化函數(shù)

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

q=1,2}

則有如下代價(jià)函數(shù)：

SSE(Θ)=

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

(5)

同時(shí)由于空間上相鄰的大腦節(jié)點(diǎn)通常有相似的fMRI時(shí)間序列和HRF函數(shù)，因此得到矩陣Vk(q)上的空間平滑懲罰項(xiàng)：

(6)

τPS(Vk(q))]

(7)

2.2 迭代算法

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

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

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

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

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

由于本文的研究重點(diǎn)是比較兩次實(shí)驗(yàn)中反應(yīng)不同的大腦節(jié)點(diǎn)而非計(jì)算HRF的具體值，即只需要選擇能清楚區(qū)分反應(yīng)相同和反應(yīng)不同腦節(jié)點(diǎn)的懲罰參數(shù)組合，因此可利用聚類方法加快懲罰參數(shù)和大腦節(jié)點(diǎn)的選擇過程，將輪廓系數(shù)作為選擇標(biāo)準(zhǔn)，在提高模型對懲罰參數(shù)容忍度的同時(shí)保證所選節(jié)點(diǎn)的準(zhǔn)確性，具體過程如下：

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

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

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

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

聚類方法提高了算法對不同懲罰參數(shù)組合的容忍度，高輪廓系數(shù)保證了所選節(jié)點(diǎn)的準(zhǔn)確性，通過采用基于K-means的快速選擇策略，利用聚類加快懲罰參數(shù)和大腦節(jié)點(diǎn)的選擇過程，可以在保證準(zhǔn)確性的同時(shí)快速找到兩次實(shí)驗(yàn)中反應(yīng)不同的大腦區(qū)域。

4 實(shí)驗(yàn)分析

4.1 實(shí)驗(yàn)數(shù)據(jù)

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

值得注意的是，本文采用的實(shí)驗(yàn)設(shè)置并不完全遵循提出的ELRMGLM，盡管如此，實(shí)驗(yàn)結(jié)果將表明ELRMGLM方法能夠超越更適合該實(shí)驗(yàn)數(shù)據(jù)的單節(jié)點(diǎn)分析方法和LRMGLM。

4.2 實(shí)驗(yàn)結(jié)果與分析

本文將提出的ELRMGLM運(yùn)用于實(shí)驗(yàn)數(shù)據(jù)，為保持運(yùn)算簡潔，令P=2來刻畫不同個(gè)體大腦節(jié)點(diǎn)HRF在振幅和延遲上的差異。圖1展示了兩次實(shí)驗(yàn)中對第二種刺激反應(yīng)不同的大腦區(qū)域，其中深色區(qū)域表示被選中節(jié)點(diǎn)，白色區(qū)域表示未選節(jié)點(diǎn)，X、Y、Z分別表示大腦節(jié)點(diǎn)的三維坐標(biāo)?？梢钥吹剑珽LRMGLM方法選中了大腦中心偏左的4×9×9網(wǎng)格，準(zhǔn)確度較高。

圖1 被選中節(jié)點(diǎn)圖像Fig. 1 Image of selected voxels

圖2 三組實(shí)驗(yàn)數(shù)據(jù)的ROC曲線Fig. 2 ROC curves of three experimental datasets

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

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

下面以第一種刺激作為參考項(xiàng)，對ELRMGLM和LRMGLM方法進(jìn)行比較，其結(jié)果如表1所示?？梢钥吹皆谌M實(shí)驗(yàn)中，ELRMGLM對兩次實(shí)驗(yàn)中反應(yīng)不同的大腦節(jié)點(diǎn)識(shí)別準(zhǔn)確度略高于LRMGLM。在計(jì)算時(shí)間上，ELRMGLM的平均迭代次數(shù)少于100而LRMGLM需要迭代上萬次；同時(shí)ELRMGLM的單次迭代時(shí)間為8 s，遠(yuǎn)小于LRMGLM的單次迭代時(shí)間60 s。不難算出，ELRMGLM的計(jì)算時(shí)間是LRMGLM的1/750，算法效率大幅提高。

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

5 結(jié)語

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