本發(fā)明屬于圖像處理領(lǐng)域，具體涉及一種基于PCA與K均值聚類的動(dòng)態(tài)腦功能連接模式分解方法，尤其涉及使用主成分分析(Principal,Component Aanlysis,PCA)和K均值聚類算法對(duì)默認(rèn)模式網(wǎng)絡(luò)、執(zhí)行控制網(wǎng)絡(luò)以及感覺網(wǎng)絡(luò)的動(dòng)態(tài)腦功能連接模式分解。

背景技術(shù)：

fMRI是研究腦活動(dòng)、腦功能的主要的無(wú)創(chuàng)方法之一，具有毫米級(jí)的空間分辨率，已成為神經(jīng)科學(xué)探索人類大腦神經(jīng)機(jī)制的重要工具。fMRI基于血氧水平依賴對(duì)比度(BOLD)間接測(cè)量神經(jīng)元活動(dòng)，它通過測(cè)量由神經(jīng)活動(dòng)引起的腦血流和腦血氧等成分變化而造成的磁共振信號(hào)變化來反應(yīng)腦活動(dòng)。BOLD-fMRI方法的發(fā)展對(duì)腦認(rèn)知功能的研究具有突破性的進(jìn)展。

靜息態(tài)fMRI刻畫了大腦在沒有任務(wù)時(shí)的神經(jīng)基準(zhǔn)活動(dòng)，具有非常重要的生理意義。研究表明在靜息態(tài)下，仍存在BOLD信號(hào)的波動(dòng)，在這些自發(fā)振蕩的信號(hào)中，存在某種特定的功能連接，構(gòu)成了靜息態(tài)功能連接網(wǎng)絡(luò)，反映了靜息態(tài)下人腦的活動(dòng)。有研究表明，神經(jīng)精神疾病如阿爾茨海默、輕度認(rèn)知受損(mild cognitive impairment)、抑郁癥(depression)以及精神分裂癥(schizophrenia)等的靜息態(tài)功能 網(wǎng)絡(luò)連接存在異常，因此研究靜息態(tài)功能網(wǎng)絡(luò)連接的異?？梢詾樯窠?jīng)精神疾病的研究提供有效和可靠的生物學(xué)指標(biāo)。靜息態(tài)功能連接網(wǎng)絡(luò)主要包含默認(rèn)模式網(wǎng)絡(luò)、感覺網(wǎng)絡(luò)(聽覺、視覺和運(yùn)動(dòng))和執(zhí)行控制網(wǎng)絡(luò)等，本發(fā)明主要針對(duì)以上三個(gè)網(wǎng)絡(luò)進(jìn)行研究。

目前大多研究只關(guān)注靜態(tài)功能連接，對(duì)于靜息態(tài)腦功能成像靜止性的假設(shè)可能忽略了重要的動(dòng)態(tài)信息。事實(shí)上，大腦本身是一個(gè)復(fù)雜的動(dòng)態(tài)系統(tǒng)，是高度變化的，并且大腦的功能連接具有瞬時(shí)變化特性，其包含更豐富的組織信息。研究大腦功能連接的時(shí)間動(dòng)態(tài)變化信息，有助于對(duì)大腦功能組織結(jié)構(gòu)更綜合和全面的認(rèn)識(shí)。

技術(shù)實(shí)現(xiàn)要素：

為了達(dá)到上述目的，本發(fā)明的目的在于提供一種基于PCA與K均值聚類的動(dòng)態(tài)腦功能連接模式分解方法，一方面能夠捕捉不同基本連接模式之間的轉(zhuǎn)變，其轉(zhuǎn)變可能與神經(jīng)認(rèn)知的變化有關(guān)；另一方面，對(duì)于僅存在于部分連接模式中的，采用靜態(tài)功能連接難以預(yù)測(cè)，因此，動(dòng)態(tài)功能連接模式分解也為生理病理學(xué)機(jī)制研究提供新的思路。

為了達(dá)到上述目的，本發(fā)明的技術(shù)方案為：

基于PCA與K均值聚類的動(dòng)態(tài)功能連接模式分解方法，具體步驟如下：

(1)、對(duì)被試進(jìn)行靜息態(tài)磁共振數(shù)據(jù)采集，并進(jìn)行預(yù)處理，磁共振數(shù)據(jù)的預(yù)處理目的是提高腦功能圖像的信噪比，并將被試圖像與標(biāo)準(zhǔn)模板進(jìn)行配準(zhǔn)變換；

(2)、預(yù)處理之后，提取DMN、ECN以及SN網(wǎng)絡(luò)的時(shí)間序列， 采用種子點(diǎn)法提取DMN23個(gè)子區(qū)的時(shí)間序列、ECN12個(gè)子區(qū)的時(shí)間序列、SN13個(gè)子區(qū)的時(shí)間序列，種子點(diǎn)區(qū)域選取半徑為6mm的球形區(qū)域，選取種子點(diǎn)之后，并對(duì)其中所有體素的信號(hào)進(jìn)行加權(quán)平均，進(jìn)而得到種子點(diǎn)的時(shí)間序列曲線；

(3)、基于PCA進(jìn)行動(dòng)態(tài)功能連接模式分解，首先通過采用滑動(dòng)窗方法計(jì)算每個(gè)被試所有腦區(qū)兩兩之間的皮爾森相關(guān)系數(shù)，進(jìn)而得到功能連接矩陣，功能連接矩陣的大小為M*N，其中M為功能連接對(duì)的個(gè)數(shù)，N為滑動(dòng)窗的個(gè)數(shù)與被試個(gè)數(shù)的乘積；其次，采用PCA方法對(duì)功能連接矩陣進(jìn)行模式分解，得到矩陣特征向量，按從大到小順序保留最上面的k個(gè)特征向量，k個(gè)特征向量表示動(dòng)態(tài)功能連接的k個(gè)主要連接模式；

(4)、為研究某些連接模式可能是準(zhǔn)穩(wěn)定的，即它們隨時(shí)間重復(fù)出現(xiàn)以及存在于眾多被試中，基于K均值聚類方法對(duì)相應(yīng)的功能連接矩陣進(jìn)行聚類，每個(gè)聚類中心代表存在的一種連接模式，實(shí)現(xiàn)基于K均值聚類的動(dòng)態(tài)功能連接模式分解。

本發(fā)明的創(chuàng)新點(diǎn)在于：提出基于PCA與K均值聚類的動(dòng)態(tài)功能連接模式分解方法，一方面能夠捕捉不同基本連接模式之間的轉(zhuǎn)變；另一方面，基本連接模式能夠預(yù)測(cè)一些神經(jīng)精神疾病(如阿爾茲海默癥、精神分裂癥、抑郁癥等)，為生理病理學(xué)機(jī)制研究提供新的思路。

附圖說明

圖1是靜息態(tài)功能圖像預(yù)處理流程圖。

圖2是頭顱分割前后對(duì)比圖。

圖3是功能像配準(zhǔn)圖，左側(cè)為標(biāo)準(zhǔn)模板，右側(cè)為配準(zhǔn)后的圖像。

圖4是基于PCA的動(dòng)態(tài)功能連接模式分解流程圖。

圖5是對(duì)DMN、ECN以及SN網(wǎng)絡(luò)的PCA模式分解結(jié)果。

圖6是對(duì)三個(gè)網(wǎng)絡(luò)總體進(jìn)行PCA模式分解的結(jié)果圖。

圖7是基于K均值聚類的動(dòng)態(tài)功能連接模式分解流程圖。

圖8是對(duì)三個(gè)網(wǎng)絡(luò)總體進(jìn)行K均值聚類模式分解的結(jié)果圖。

圖9是其中3個(gè)被試的狀態(tài)分配圖。

圖10是各腦區(qū)種子點(diǎn)坐標(biāo)。

圖11是所有被試的平均轉(zhuǎn)移矩陣。

具體實(shí)施方式

下面結(jié)合附圖對(duì)本發(fā)明做詳細(xì)敘述。

本發(fā)明基于PCA與K均值聚類對(duì)靜息態(tài)動(dòng)態(tài)功能連接進(jìn)行功能連接模式分解。

1、首先對(duì)采集的原始靜息態(tài)磁共振數(shù)據(jù)進(jìn)行預(yù)處理，由于磁共振掃描過程中各種各樣的噪聲的影響，個(gè)體自身存在尺度和位置上的差異，非常有必要在分析數(shù)據(jù)之前對(duì)數(shù)據(jù)做一定的預(yù)處理。在整個(gè)的實(shí)驗(yàn)的數(shù)據(jù)獲取中，主要的噪聲信息來源包括：(1)物理頭動(dòng)；(2)圖像內(nèi)層間掃描時(shí)間差別；(3)外在磁場(chǎng)的不均勻性。腦功能圖像預(yù)處理是在保留腦功能圖像細(xì)節(jié)的同時(shí)，使用腦功能圖像與標(biāo)準(zhǔn)模板進(jìn)行仿射配準(zhǔn)變換方式的預(yù)處理，并提高腦功能圖像的信噪比。

功能磁共振數(shù)據(jù)的預(yù)處理選取了Linux系統(tǒng)Ubuntu12.04下的AFNI和FSL軟件，并通過編寫B(tài)atch文件實(shí)現(xiàn)數(shù)據(jù)的批處理。預(yù)處理 流程見圖1，主要包含以下幾個(gè)方面：

1)結(jié)構(gòu)像頭顱分割

采集得到的結(jié)構(gòu)像通常包含頭顱信息，需要進(jìn)行頭顱分割，從而消除顱骨部位引入的偽影對(duì)后續(xù)數(shù)據(jù)分析的影響。基于FSL軟件提供的3drefit、3dresample和fast segment實(shí)現(xiàn)頭部顱骨和大腦內(nèi)部組織結(jié)構(gòu)的分割，分割結(jié)果如圖2所示。

2)時(shí)間對(duì)齊

血液動(dòng)力學(xué)函數(shù)表明血液對(duì)刺激的響應(yīng)有一定的時(shí)間延遲，由于一個(gè)TR期間采集到全腦圖像，這導(dǎo)致每一層的圖像并非在同一時(shí)刻采集，而是發(fā)生在整個(gè)掃描時(shí)間段內(nèi)，時(shí)間校正就是通過類似于插值的方法對(duì)每一層圖像進(jìn)行層時(shí)間處理，使得一個(gè)TR周期內(nèi)各層圖像近似在同一時(shí)刻獲取。輸入每幀圖像的各掃描參數(shù)進(jìn)行時(shí)間校正，包括掃描時(shí)間、掃描層數(shù)、第一層到最后一層的時(shí)間間隔以及掃描順序，以中間層為參考進(jìn)行時(shí)間校正，將每幀圖像校正為幾乎是同一時(shí)間點(diǎn)采集的各層圖像。

3)頭動(dòng)校正

為避免掃描過程中發(fā)生微小頭動(dòng)，對(duì)掃描質(zhì)量及結(jié)果造成影響，以每個(gè)被試的第一幀圖像為基準(zhǔn)，將其余的所有圖像和基準(zhǔn)圖像對(duì)齊。進(jìn)行頭動(dòng)校正時(shí)一般將被試的大腦看作一個(gè)剛體，因此在fMRI實(shí)驗(yàn)中被試頭部的運(yùn)動(dòng)可以近似成一種剛體運(yùn)動(dòng)，即只有平移變換與旋轉(zhuǎn)變換的組合。通過AFNI的3dvolreg函數(shù)使得其余的所有圖像與參考圖像配準(zhǔn)，如果頭動(dòng)超過一個(gè)體素則去除該被試。

4)空間平滑

對(duì)圖像進(jìn)行高斯平滑，主要目的是提高圖像信噪比以及確保圖像數(shù)據(jù)具有隨機(jī)高斯場(chǎng)性質(zhì)，這對(duì)于需要做平均以及統(tǒng)計(jì)結(jié)果非常重要?？臻g平滑采用高斯函數(shù)進(jìn)行高斯平滑，能有效地消弱隨機(jī)噪聲對(duì)fMRI信號(hào)的影響，提高數(shù)據(jù)的信噪比。三維高斯函數(shù)是比較常用的空間平滑方法，其半高全寬決定了空間平滑的力度，本文選用半高全寬(Full Width at Half Maximum,FWHM)為6mm的高斯核函數(shù)進(jìn)行數(shù)據(jù)平滑。

5)時(shí)域帶通濾波

靜息態(tài)fMRI信號(hào)的低頻波動(dòng)反映了自發(fā)的神經(jīng)活動(dòng)，因此采用帶通為0.01-0.1Hz的帶通濾波器去除與呼吸、心跳等有關(guān)的生理噪聲。

6)去線性

由于機(jī)器的長(zhǎng)期工作而升溫或者被試的不適應(yīng)，隨著時(shí)間的積累會(huì)存在一個(gè)線性趨勢(shì)。線性趨勢(shì)也屬于噪聲引起的，因此需進(jìn)行去除。

7)圖像分割

為了去除腦脊液(Cerebro-Spinal Fluid,CSF)、白質(zhì)(White Matter,WM)等冗余信號(hào)，需要對(duì)結(jié)構(gòu)像進(jìn)行分割，利用分割得到的信息制作腦脊液與白質(zhì)模板；

8)冗余去除

去除白質(zhì)、腦脊液、全腦信號(hào)以及頭動(dòng)偽跡等冗余信號(hào)。

2、采用種子點(diǎn)方法在標(biāo)準(zhǔn)MNI空間中提取被試腦區(qū)時(shí)間序列，包含DMN(23個(gè)子區(qū))、ECN(12個(gè)子區(qū))、SN(13個(gè)子區(qū))的時(shí) 間序列。為了準(zhǔn)確定位各腦區(qū)，首先將靜息態(tài)磁共振圖像配準(zhǔn)至MNI標(biāo)準(zhǔn)空間；為提高配準(zhǔn)精度，采用兩步配準(zhǔn)：首先將功能像配準(zhǔn)到結(jié)構(gòu)像，其次將結(jié)構(gòu)像配準(zhǔn)到標(biāo)準(zhǔn)空間，從而利用所得到的變換矩陣將功能像配準(zhǔn)到標(biāo)準(zhǔn)空間并將體素重采樣為3mm×3mm×3mm；配準(zhǔn)模板采用MNI模板，該模板是由加拿大蒙特利爾神經(jīng)科學(xué)研究所(Montreal Neurological Institute,MNI)研發(fā)而成；整個(gè)配準(zhǔn)過程使用FMRIB提供的線性配準(zhǔn)工具實(shí)現(xiàn)，配準(zhǔn)后圖像如圖3所示?；趂MRI的研究的種子點(diǎn)區(qū)域選取以各子腦區(qū)的坐標(biāo)為圓心，半徑為6mm的球形區(qū)域，種子點(diǎn)區(qū)域包含33個(gè)體素，各腦區(qū)種子點(diǎn)坐標(biāo)如圖10所示。選取種子點(diǎn)之后，分別對(duì)每一個(gè)被試提取所有種子點(diǎn)區(qū)域內(nèi)全部體素的一維信號(hào)，并對(duì)其求加權(quán)平均作為種子點(diǎn)的時(shí)間序列曲線；整個(gè)基于種子點(diǎn)的時(shí)間序列提取過程基于FSL軟件實(shí)現(xiàn)。

3、接下來基于PCA方法進(jìn)行動(dòng)態(tài)功能連接模式分解。PCA的一般過程是：

1)去除平均值：means；

2)計(jì)算A的協(xié)方差矩陣Σ；

3)計(jì)算Σ的特征向量和特征值；

4)將特征值從大到小排序；

5)保留最上面的k個(gè)特征向量(這k個(gè)特征向量保證了數(shù)據(jù)映射到特征值最大的特征向量的方向時(shí)，數(shù)據(jù)間的累積方差最大，數(shù)據(jù)映射到第二大的特征向量時(shí)，數(shù)據(jù)間的累積方差次之)。

整個(gè)PCA動(dòng)態(tài)功能連接模式分解流程如圖4所示。首先采用滑動(dòng) 窗方法計(jì)算各腦區(qū)時(shí)間序列兩兩之間的相關(guān)系數(shù)，進(jìn)而構(gòu)造基于PCA分析的功能連接矩陣。本文采用的滑動(dòng)窗窗長(zhǎng)為90s，滑動(dòng)步長(zhǎng)為1TR＝645ms，相關(guān)系數(shù)采用皮爾森相關(guān)系數(shù)。相關(guān)系數(shù)ρ考察兩個(gè)變量的相關(guān)程度，取值范圍在-1～1之間，其中，1表示變量完全正相關(guān)，0表示無(wú)關(guān)，-1表示完全負(fù)相關(guān)。

${\mathrm{\ρ}}_{X,Y}=\frac{\mathrm{cov}(X,Y)}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}=\frac{E\[(X-{\mathrm{\μ}}_X)(Y-{\mathrm{\μ}}_Y)\]}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y}$

式中：ρ(X,Y)——變量X和Y的相關(guān)系數(shù)；cov(X,Y)——變量X和Y的協(xié)方差；σ_X——變量X的標(biāo)準(zhǔn)方差；σ_Y——變量Y的標(biāo)準(zhǔn)方差；μ_X——變量X的均值；μ_Y——變量Y的均值。

功能連接矩陣的具體構(gòu)造過程如下，首先計(jì)算每個(gè)被試DMN、ECN和SN網(wǎng)絡(luò)全部48個(gè)子區(qū)兩兩之間的動(dòng)態(tài)功能連接，得到W個(gè)48×48的功能連接矩陣，其中W為滑動(dòng)窗的個(gè)數(shù)。其次，取每個(gè)功能連接矩陣的上三角進(jìn)行展開，并將展開后的矩陣以時(shí)間順序進(jìn)行連接，得到每個(gè)被試的功能連接矩陣C_S(S＝1,2,...,s)，為減少被試個(gè)體之間的差異，對(duì)每個(gè)被試的C_S減去C_S每行的平均值，即得到最后，將各個(gè)被試的去均值功能連接矩陣進(jìn)行連接，得到用于PCA分析的功能連接矩陣

得到功能連接矩陣后，對(duì)功能連接矩陣進(jìn)行PCA分析，得到分解后的連接模式，這些連接模式盡可能多地功能連接矩陣的信息，為基本連接模式。首先計(jì)算矩陣A的協(xié)方差矩陣AA^T的特征值和特征向量：AA^T＝U∧U^T，其中U的各列向量為正交特征向量，∧的對(duì)角元素為對(duì) 應(yīng)于特征向量的特征值。這些特征向量即為矩陣的主成分，表示功能連接矩陣的主要特征，為基本連接模式，因此稱特征向量為“特征連接”。

正交矩陣U的大小為(N²-N)/2×(N²-N)/2，其中N表示腦區(qū)的數(shù)量，即得到(N²-N)/2個(gè)特征連接。然而，對(duì)應(yīng)的特征值大的特征向量表示的動(dòng)態(tài)功能連接矩陣的累積方差最大，因此，少數(shù)的特征連接即可有效地代表矩陣的特征。將特征值按從大到小排列，取前10個(gè)特征值對(duì)應(yīng)的特征向量來表示功能連接模式分解結(jié)果。

整個(gè)基于PCA的動(dòng)態(tài)功能連接模式分解流程使用MATLAB編寫程序?qū)崿F(xiàn)。對(duì)DMN、ECN以及SN網(wǎng)絡(luò)的PCA模式分解結(jié)果如圖5，對(duì)三個(gè)網(wǎng)絡(luò)總體進(jìn)行模式分解的結(jié)果如圖6。

4、為研究某些連接模式可能是準(zhǔn)穩(wěn)定的，即它們隨時(shí)間重復(fù)出現(xiàn)以及存在于眾多被試中，基于K均值聚類方法對(duì)相應(yīng)的功能連接矩陣進(jìn)行聚類分析。

整個(gè)K均值聚類的動(dòng)態(tài)功能連接模式分解流程如圖7所示。首先采用滑動(dòng)窗方法計(jì)算各腦區(qū)時(shí)間序列兩兩之間的相關(guān)系數(shù)，進(jìn)而構(gòu)造基于K均值聚類分析的功能連接矩陣。功能連接矩陣的具體構(gòu)造過程如下，首先計(jì)算每個(gè)被試DMN、ECN和SN網(wǎng)絡(luò)全部48個(gè)子區(qū)兩兩之間的動(dòng)態(tài)功能連接，得到W個(gè)48×48的功能連接矩陣，其中W為滑動(dòng)窗的個(gè)數(shù)。其次，取每個(gè)功能連接矩陣的上三角進(jìn)行展開，并將展開后的矩陣以窗的時(shí)間順序進(jìn)行連接。最后，將各個(gè)被試的功能連接矩陣進(jìn)行連接，得到用于K均值聚類分析的功能連接矩陣B，B的 大小為P×Q，其中P為窗的個(gè)數(shù)與被試個(gè)數(shù)的乘積，Q為功能連接對(duì)的個(gè)數(shù)。

得到功能連接矩陣后，采用K均值聚類方法對(duì)功能連接矩陣進(jìn)行行聚類分析，具體為：

1)從n個(gè)數(shù)據(jù)對(duì)象中任意選擇k個(gè)對(duì)象作為初始聚類中心M_i(i＝1,2,...,k)。

2)根據(jù)每個(gè)聚類對(duì)象的初值(中心對(duì)象)，計(jì)算每個(gè)對(duì)象X_k(k＝1,2,...,n)與這些中心對(duì)象的距離D_k＝|M_i-X_k|，并根據(jù)最小距離重新對(duì)相應(yīng)的對(duì)象進(jìn)行劃分，形成類簇C_i(i＝1,2,...N_i)，N_i表示類簇C_i中包含的對(duì)象個(gè)數(shù)。

3)計(jì)算每個(gè)聚類的均值作為更新后的聚類中心：

$M_i=\underset{X_k\∈C_i}{\Σ}{X}_k/N_i$

4)重復(fù)步驟(2)和步驟(3)，直到每個(gè)聚類不再發(fā)生變化為止。

整個(gè)基于K均值聚類法的動(dòng)態(tài)功能連接模式分解流程使用MATLAB編寫程序?qū)崿F(xiàn)。每個(gè)聚類中心代表存在的一種連接模式，即一種“狀態(tài)”，本文選取的聚類個(gè)數(shù)為7，對(duì)DMN、ECN以及SN三個(gè)網(wǎng)絡(luò)總體進(jìn)行K均值聚類模式分解的結(jié)果如圖8。

接下來根據(jù)聚類結(jié)果可以得到每個(gè)被試以時(shí)間為函數(shù)的狀態(tài)分配圖。圖9所示為其中3個(gè)被試的狀態(tài)分配圖，可以看出，功能連接在長(zhǎng)時(shí)間內(nèi)趨于單一狀態(tài)，盡管在狀態(tài)轉(zhuǎn)換過渡時(shí)的短時(shí)間內(nèi)也通常會(huì)出現(xiàn)其它狀態(tài)。

除了狀態(tài)轉(zhuǎn)移圖，也可采用馬爾可夫鏈描述其轉(zhuǎn)換行為，將狀態(tài)轉(zhuǎn)換過程看作馬爾科夫鏈。馬爾可夫鏈具體為：設(shè){X(t),t∈T}是隨機(jī)過程，其參數(shù)集T＝{0,1,2,...}，狀態(tài)空間S為可數(shù)集，設(shè)S＝{1,2,...}，若對(duì)任意的m≥1及任意的i₀,i₁,...,i_m,j∈S，當(dāng) 時(shí)，有

P{X(m+1)＝j(luò)|X(0)＝i₀,X(1)＝i₁,...,X(m)＝i_m}

＝P{X(m+1)＝j(luò)|X(m)＝i_m}

則稱{X(t),t∈T}為離散時(shí)間馬爾可夫鏈。通過馬爾可夫鏈得到每個(gè)被試的轉(zhuǎn)移矩陣，轉(zhuǎn)移矩陣中的元素表示由一種狀態(tài)轉(zhuǎn)移到另一種狀態(tài)的概率。然后將所有被試的轉(zhuǎn)移矩陣進(jìn)行平均，得到平均轉(zhuǎn)移矩陣，其矩陣的每一行元素乘以被試個(gè)數(shù)然后相加，結(jié)果為整數(shù)。

綜上所述，本發(fā)明目的在于提出基于PCA與K均值聚類的動(dòng)態(tài)功能連接模式方法，不僅能夠捕捉不同基本連接模式之間的轉(zhuǎn)變，并且為臨床上精神分裂癥、阿爾茲海默癥、抑郁癥等神經(jīng)精神疾病的進(jìn)一步研究和防治提供策略。