基于獨(dú)立成分分析的時(shí)間自相關(guān)方法在功能磁共振激活區(qū)自適應(yīng)提取中的應(yīng)用

白保東 劉 健 郭紅宇

(沈陽工業(yè)大學(xué)生物醫(yī)學(xué)與電磁工程研究所， 沈陽 110870)

基于獨(dú)立成分分析的時(shí)間自相關(guān)方法在功能磁共振激活區(qū)自適應(yīng)提取中的應(yīng)用

白保東*劉 健 郭紅宇

(沈陽工業(yè)大學(xué)生物醫(yī)學(xué)與電磁工程研究所， 沈陽 110870)

傳統(tǒng)基于ICA的激活區(qū)檢測手段是將分離后的獨(dú)立成分與參考信號做相關(guān)性分析。實(shí)際問題中，不同區(qū)域的腦血流動力學(xué)響應(yīng)情況不同，因此往往得不到標(biāo)準(zhǔn)的參考信號。針對此類問題，提出時(shí)間自相關(guān)方法(TSC)與ICA方法結(jié)合，在不需要參考信號的情況下，通過檢測體素點(diǎn)各周期的時(shí)間序列相關(guān)性，對fMRI數(shù)據(jù)進(jìn)行激活區(qū)提取。應(yīng)用5-鄰域ICA方法對fMRI數(shù)據(jù)逐點(diǎn)處理，然后應(yīng)用時(shí)間自相關(guān)算法檢測各時(shí)間序列周期間的相關(guān)性，選擇最大的自相關(guān)系數(shù)作為該體素點(diǎn)的信號值。再通過Z變換將相關(guān)系數(shù)分布轉(zhuǎn)換為服從N(0,1)的Z分布，提取出具有顯著性差異(a=0.05)的激活區(qū)。將自相關(guān)算法應(yīng)用于仿真數(shù)據(jù)和12組雙手握拳運(yùn)動的真實(shí)fMRI數(shù)據(jù)的處理，結(jié)果表明該方法能夠準(zhǔn)確提取出仿真數(shù)據(jù)中的激活區(qū)。對真實(shí)數(shù)據(jù)的處理，該方法在空間準(zhǔn)確性上與GLM方法無顯著性差別(0.465 3±0.136 8 vs 0.490 5±0.134 1)，在時(shí)間準(zhǔn)確性上顯著優(yōu)于GLM方法 (0.636 4±0.011 1 vs 0.369 2±0.010 9)，具有良好的腦功能激活區(qū)檢測及空間定位能力。

獨(dú)立分量分析(ICA)；fMRI；時(shí)間自相關(guān)；假設(shè)檢驗(yàn); 自適應(yīng)閾值

引言

獨(dú)立分量分析(independent component analysis，ICA)是信號處理領(lǐng)域在20世紀(jì)90年代后期發(fā)展起來的一項(xiàng)新處理方法[1]，可以在源信號和傳輸通道參數(shù)未知的情況下把觀測信號分解成若干個(gè)互相獨(dú)立的成分，現(xiàn)已成為信號處理、神經(jīng)網(wǎng)絡(luò)(神經(jīng)計(jì)算)、統(tǒng)計(jì)學(xué)的研究熱點(diǎn)[2-3]。國際上應(yīng)用最廣的是FastICA算法和Infomax算法，其中Infomax算法一種基于信息最大化傳輸?shù)膯螌忧胺答伾窠?jīng)網(wǎng)絡(luò)算法，它的基本思想是利用了似然函數(shù)極大法估計(jì)中要估計(jì)x的獨(dú)立性最大等價(jià)于估計(jì)x的某一非線性函數(shù)獨(dú)立性的最大化法的準(zhǔn)則[4]。

以血氧水平依賴測量為原理的功能磁共振成像技術(shù) (blood oxygenation level dependent-functional magnetic resonance imaging，BOLD-fMRI)是20世紀(jì)90年代在MRI基礎(chǔ)上發(fā)展起來的一種新的腦功能成像技術(shù)，既保留普通MRI的解剖學(xué)成像特點(diǎn)，又可同時(shí)獲得生理學(xué)信息[5]。BOLD-fMRI技術(shù)的出現(xiàn)，為研究腦卒中后康復(fù)的機(jī)制、評價(jià)和判斷預(yù)后提供了嶄新的思路。

McKeown等把腦神經(jīng)活動信號和非自然信號看成是相互獨(dú)立的信號，最早將ICA應(yīng)用在fMRI領(lǐng)域[6]。近幾年來，學(xué)者們紛紛探索ICA 算法來分離FMRI 信號,目前使用的ICA 算法有Orth2Infomax 算法、Fast ICA算法、空間獨(dú)立分量分析算法、Group-ICA算法、Orth-ExtBS 算法、牛頓型算法等。

基于時(shí)間序列中自相似特性的時(shí)間自相關(guān)(temporal self-correlation，TSC)方法是一種激活區(qū)提取方法。該方法通過fMRI圖像中每個(gè)體素在時(shí)間序列不同周期(相鄰一個(gè)靜息和一個(gè)刺激為一個(gè)周期)之間求相關(guān)系數(shù)，計(jì)算出TSC值，并與特定臨界值比較，判斷該體素點(diǎn)是否激活[7]。在fMRI應(yīng)用中，傳統(tǒng)基于ICA的激活區(qū)檢測方法是將分離后的獨(dú)立成分與參考信號進(jìn)行相關(guān)分析。參考信號是通過刺激周期與血流動力學(xué)相應(yīng)函數(shù)(HRF)做卷積得到的。在實(shí)際問題中，不同受試者、不同區(qū)域、不同狀態(tài)下的腦血流動力學(xué)響應(yīng)情況也會不同[8]，因此往往得不到準(zhǔn)確地參考信號，而全部采用標(biāo)準(zhǔn)參考信號將會造成極大的個(gè)體誤差。本研究提出時(shí)間自相關(guān)方法(TSC)與ICA方法結(jié)合的ICA-TSC方法，在不需要參考信號的情況下，通過檢測體素點(diǎn)各周期的時(shí)間序列相關(guān)性，對fMRI數(shù)據(jù)進(jìn)行激活信息檢測。

采用Z變換方法，將相關(guān)系數(shù)轉(zhuǎn)換成服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)的Z參數(shù)，然后根據(jù)檢驗(yàn)的顯著性水平(α=0.05)閾值化統(tǒng)計(jì)參數(shù)圖，以判別體素的激發(fā)與未激發(fā)。判斷fMRI結(jié)果的空間準(zhǔn)確性沒有一個(gè)完全正確的標(biāo)準(zhǔn)，而SPM是當(dāng)前國際上主流的fMRI數(shù)據(jù)分析軟件之一，它基于像素分析法，以腦功能的像素為基本單位，逐像素地對數(shù)據(jù)做一般線性模型(general linear model，GLM)的統(tǒng)計(jì)分析，得到某個(gè)顯著性水平下的腦激活圖，并能以此為基礎(chǔ)，給出某些(或全部) 激活區(qū)的生理信號時(shí)間變化曲線，即時(shí)間序列[9]。通過SPM軟件得到群體統(tǒng)計(jì)結(jié)果，其正確性被國際學(xué)者所廣泛接受和肯定。筆者將以此作為參照標(biāo)準(zhǔn)來檢驗(yàn)新算法的空間準(zhǔn)確性，同時(shí)以SPM軟件得到的時(shí)域解釋變量為參照標(biāo)準(zhǔn)，檢驗(yàn)新算法的時(shí)間準(zhǔn)確性。

1 方法

1.1ICA算法流程

標(biāo)準(zhǔn)ICA學(xué)模型為

x=As

(1)

混合矩陣A和源信號向量s都是未知的，但只要設(shè)法使分離矩陣B恢復(fù)的各個(gè)輸出成分之間相互統(tǒng)計(jì)獨(dú)立，就可以實(shí)現(xiàn)信號源的盲分離。所以，ICA實(shí)際是在某一獨(dú)立性判據(jù)意義下進(jìn)行的尋優(yōu)計(jì)算，求解ICA可歸結(jié)為兩步：第一步是建立目標(biāo)函數(shù)，第二步是尋求求解該目標(biāo)函數(shù)的算法。

目前，常用的幾類衡量信號之間獨(dú)立性的目標(biāo)函數(shù)有：非高斯性目標(biāo)函數(shù)、互信息目標(biāo)函數(shù)、極大似然估計(jì)目標(biāo)函數(shù)、信息極大化目標(biāo)函數(shù)。

Infomax算法采用信息極大化的目標(biāo)函數(shù)。它的特點(diǎn)是在輸出y之后逐分量地引入一個(gè)非線性函數(shù)ri=gi(yi)來代替對高階統(tǒng)計(jì)量的估計(jì)，在給定合適的gi(yi)后，通過調(diào)節(jié)解混矩陣B,使輸出r=[r1,r2,…,rM]T的總熵量H(r，B)極大,則其目標(biāo)函數(shù)為ε=H(r，B)。以下是對H(r，B)的推導(dǎo)：

由熵(信號中所含有的平均信息量)的定義有

(2)

熵的單位依對數(shù)函數(shù)log的底數(shù)而定：當(dāng)a=2時(shí)其單位為比特(bit)，當(dāng)a=e時(shí)單位為奈特(nat)，而當(dāng)a=10時(shí)單位則為哈特(hat)，其中ri=gi(yi)

(3)

將式(3)代入式(2)，有

(4)

式(4)就是所求的目標(biāo)函數(shù)ε。

Infomax算法對于非線性函數(shù)g(·)的選擇要求不是很苛刻，某些單調(diào)增長函數(shù)都可以用，如sigmoid函數(shù)、tanh函數(shù)等。

選擇sigmoid函數(shù)，則有

(5)

y=Bx+B0

(6)

將式(5)中的各式代入式(4)，并對兩邊求導(dǎo)，有

(7)

為避免矩陣求逆問題，用自然梯度來代替常規(guī)梯度，則式(7)變?yōu)?/p>

(8)

將式(8)代入式(6)，有

(9)

式中，μ為步長，I為元素均為1的向量。

判斷Infomax算法是否收斂的終止條件為：迭代前的B和迭代后的B同指一個(gè)方向，即它們的內(nèi)積接近或等于1[3]。

算法的主要步驟如下：

步驟1：初始化解混矩陣B(隨機(jī)數(shù))。

步驟2：取第i個(gè)樣本矢量xi(i=1,2,…,m)。

1)計(jì)算解向量yi=Bxi和網(wǎng)絡(luò)輸出，有

(10)

2)計(jì)算權(quán)值增量，有

3)更新權(quán)值B=B+ΔB

步驟3：判斷是否達(dá)到收斂條件，“是”表示結(jié)束；否則回到步驟2。

1.2激活區(qū)檢測方法

應(yīng)用ICA判斷fMRI圖像上的激活區(qū)，可以采用時(shí)域ICA(temporal ICA，tICA)或空域ICA(spatial ICA，sICA)。無論是tICA還是sICA，在成分分離算法上沒有明顯區(qū)別，主要區(qū)別在于排列為二維數(shù)組的方式[10]。

fMRI數(shù)據(jù)反映的是三維空間結(jié)構(gòu)和時(shí)間序列組成的四維數(shù)據(jù)，tICA是將同一時(shí)間點(diǎn)上的各空間數(shù)據(jù)排成一列，同一體素的不同時(shí)間數(shù)據(jù)排成一行，sICA則相反。sICA是在一定程度上犧牲了時(shí)間的獨(dú)立性來換取空間獨(dú)立性，而tICA 則是犧牲空間獨(dú)立性來換取時(shí)間的獨(dú)立性。

在fMRI數(shù)據(jù)分析中，感興趣的效應(yīng)主要是與運(yùn)動任務(wù)有關(guān)的信號。與感興趣效應(yīng)無關(guān)的其他因素主要有兩個(gè)來源：一個(gè)是生理噪聲，主要由心跳、呼吸和血流等引起；一個(gè)是熱噪聲和成像噪聲，由電磁干擾、受試者的頭部運(yùn)動等因素引起。通過合理設(shè)計(jì)實(shí)驗(yàn)，可以保證任務(wù)頻率與生理噪聲的頻率不同，則生理噪聲所產(chǎn)生的信號與任務(wù)所產(chǎn)生的信號在時(shí)域內(nèi)可以看作是獨(dú)立的[11]。由于生理噪聲對腦的影響是全局性的，所以它們所產(chǎn)生的干擾效應(yīng)與任務(wù)效應(yīng)在空間上的獨(dú)立性廣受質(zhì)疑。同樣，熱噪聲和成像噪聲也不可避免地會影響到全腦的fMRI信號?？偠灾?，干擾效應(yīng)與任務(wù)效應(yīng)在空間上可能是相關(guān)的，但它們的時(shí)間特性卻仍然可能是相互獨(dú)立的。因此，假設(shè)fMRI數(shù)據(jù)中可能的信號源是時(shí)間獨(dú)立而不是空間獨(dú)立更為合理，并且相比較sICA而言，tICA的物理意義更為明確[12]。

鄰域tICA采用逐點(diǎn)處理法，選擇大腦切片中某一體素點(diǎn)與N-1個(gè)相鄰體素點(diǎn)的時(shí)間序列共同構(gòu)成N維觀測向量。ICA分解后，得到N個(gè)獨(dú)立時(shí)間成分。引入TSC方法，計(jì)算該體素點(diǎn)的所有獨(dú)立時(shí)間成分的自相關(guān)性。選擇自相關(guān)性最大的成分進(jìn)行假設(shè)檢驗(yàn)，判斷該點(diǎn)是否被激活。

1.3TSC方法

設(shè)每個(gè)時(shí)間序列有m個(gè)周期，每個(gè)周期有n個(gè)時(shí)間點(diǎn)，則數(shù)據(jù)可用矩陣來表示，有：

(11)

式中，矩陣元xij表示對于某一確定位置的體素點(diǎn)在第i個(gè)周期中第j幅圖像內(nèi)該體素點(diǎn)的強(qiáng)度，則有

(12)

用來描述時(shí)間序列的第i個(gè)周期。

對于一個(gè)時(shí)間序列t=[t1,t2,…,tm]，對其中的所有周期兩兩求相關(guān)系數(shù)，有

(13)

式中，corr(ti,tj)是兩個(gè)不同周期數(shù)據(jù)的協(xié)方差，D(ti)和D(tj)分別為ti和tj的方差。

TSC值被定義為cij的數(shù)學(xué)期望，有

(14)

這種通過計(jì)算體素點(diǎn)在時(shí)間序列中各周期間的相關(guān)系數(shù)來判定該體素點(diǎn)是否激活的方法，稱為時(shí)間自相關(guān)方法。對于激活的體素點(diǎn)，其各周期間具有較強(qiáng)的相關(guān)性。將相關(guān)系數(shù)最大的成分組合在一個(gè)二維數(shù)組里，構(gòu)成由最大相關(guān)系數(shù)組成的統(tǒng)計(jì)參數(shù)圖，這種成分叫做一致任務(wù)相關(guān)(consistently task-related，CTR)成分[6]。

對Tij做Z變換

(15)

式中，μ為所有體素點(diǎn)強(qiáng)度的均值，σ為標(biāo)準(zhǔn)差，zij服從標(biāo)準(zhǔn)正態(tài)分布N(0,1)[13]。

接著通過單尾z檢驗(yàn)，計(jì)算zij相應(yīng)的p-值，以檢測Tij的顯著性。若Tij不顯著，說明該體素各周期時(shí)間序列之間不相關(guān)，接受原假設(shè)H0；若Tij顯著，說明該體素各周期時(shí)間序列之間相關(guān)，拒絕原假設(shè)H0。

1.4ICA-TSC方法準(zhǔn)確性分析

群體統(tǒng)計(jì)結(jié)果是對所有GLM方法所得個(gè)體結(jié)果的二次統(tǒng)計(jì)，黑色激活表示在所有個(gè)體結(jié)果中該體素點(diǎn)被激活的次數(shù)高于某一閾值。以SPM群體統(tǒng)計(jì)結(jié)果作為參照標(biāo)準(zhǔn)，分別與ICA-TSC方法和GLM方法的個(gè)體結(jié)果做比較，將各自得到的相關(guān)系數(shù)(R,R0)作為準(zhǔn)確性指標(biāo)。應(yīng)用z檢驗(yàn)方法，檢驗(yàn)ICA-TSC方法與GLM方法的可信度是否具有顯著性差異，若沒有，是否極其相似。通常的顯著性差異檢驗(yàn)，注重“棄真”錯(cuò)誤概率。不同于顯著性差異檢驗(yàn)，衡量兩個(gè)數(shù)據(jù)集相似程度的假設(shè)檢驗(yàn)更注重“取偽”錯(cuò)誤的概率。例如，設(shè)“取偽”錯(cuò)誤概率

(16)

若z>k，則ICA-TSC方法具有極高的可信度。

時(shí)域解釋變量反映實(shí)驗(yàn)設(shè)計(jì)的血流動力學(xué)響應(yīng)曲線。以時(shí)域解釋變量為參照標(biāo)準(zhǔn)，分別與兩種方法所得的峰值點(diǎn)時(shí)間序列做比較，同樣將各自得到的相關(guān)系數(shù)作為準(zhǔn)確性指標(biāo)。應(yīng)用t檢驗(yàn)方法檢驗(yàn)兩種方法誰的時(shí)間準(zhǔn)確性更高，是否有顯著性差異，即

(17)

2 實(shí)驗(yàn)

2.1仿真數(shù)據(jù)分析

為了驗(yàn)證算法的可靠性，構(gòu)造一個(gè)分辨率為64像素×64像素的橫斷面仿真fMRI圖像,共300張。每3層組合，模擬成一個(gè)三維數(shù)據(jù)，則共100個(gè)時(shí)間點(diǎn)。靜息-任務(wù)-靜息相互交替，每個(gè)靜息或任務(wù)模塊為10個(gè)時(shí)間點(diǎn)，即20個(gè)時(shí)間點(diǎn)為一周期，共5個(gè)周期，如圖1所示。靜息態(tài)為像素為100的均勻圖像，4個(gè)白色方塊為預(yù)先設(shè)定的模擬激活區(qū)域或噪聲區(qū)域，如圖2所示。在仿真數(shù)據(jù)中，各區(qū)域信號仿真如圖3所示。

圖1 時(shí)間序列設(shè)計(jì)Fig.1 Time series design

圖2 仿真數(shù)據(jù)位置分布Fig.2 simulation data location distribution

圖3 仿真數(shù)據(jù)中各部分信號的時(shí)間序列波形。(a)激活信號,2%；(b)系統(tǒng)噪聲(高斯白噪聲),5%；(c)生理噪聲(包括心跳、呼吸等)；(d)以上3種信號的混合信號Fig.3 The time series signal waveform of each component of the simulation data. (a): Active signal,2%; (b): System noise(Gaussian white noise),5%; (c): Physiological noise (including heartbeat, breathing, etc.); (d): Mixed signal of the above three signals

區(qū)域A與區(qū)域D中包含激活信號，B與C中全部是噪聲信號。處理層的選擇是任意的，在本實(shí)驗(yàn)中則選擇了第2層作為處理層，處理后得到圖4所示的統(tǒng)計(jì)參數(shù)。像素亮度代表相關(guān)系數(shù)的數(shù)值大小，越亮則表示值越大，反之亦然。

圖4 由相關(guān)系數(shù)組成的統(tǒng)計(jì)參數(shù)Fig.4 The correlation coefficient of statistical parameter diagram

圖5表示相關(guān)系數(shù)的分布情況，從0～1分別是B區(qū)、C區(qū)、D區(qū)、A區(qū)，A區(qū)最大、B區(qū)最小。因此，只要選擇合適的閾值，就可以實(shí)現(xiàn)含有激活成分的區(qū)域的提取，手動選擇T=0.2。

圖5 相關(guān)系數(shù)分布Fig.5 Correlation coefficient distribution

相關(guān)系數(shù)大于T則保留原值，小于T則置零。

圖6 仿真實(shí)驗(yàn)結(jié)果Fig.6 The simulation experiment results

圖7 分離后D區(qū)的激活信號Fig.7 Active signal of the D area after separation

2.2雙手握拳運(yùn)動激活區(qū)提取

2.2.1數(shù)據(jù)來源及類型

材料為健康右利手志愿者12例(8名男性，4名女性)，年齡20～28歲。所有受試者均無神經(jīng)、精神損傷史，受試前無任何不適感，并且未經(jīng)歷過樂器演奏等指法訓(xùn)練。

在功能圖像的采集中，飛利浦1.5 T掃描儀，EPI序列，TR為250 0 ms，TE為50 ms，矩陣為64×64×47，層厚為4 mm，層距為4 mm。FOV為200 mm，Voxel size為4×4×3。實(shí)驗(yàn)設(shè)計(jì)采用block設(shè)計(jì)，任務(wù)為主動雙手握拳運(yùn)動，頻率約1 Hz。先休息，然后運(yùn)動，再休息，反復(fù)5個(gè)周期，block長度為12個(gè)TR時(shí)間(30 s)，共10個(gè)block(120個(gè)采樣時(shí)間點(diǎn))。

2.2.2fMRI數(shù)據(jù)預(yù)處理

首先對磁共振系統(tǒng)掃描得到的DICOM格式數(shù)據(jù)轉(zhuǎn)換成.img格式數(shù)據(jù)，利用SPM軟件進(jìn)行預(yù)處理工作，包括頭動校正、圖像配準(zhǔn)、空間平滑等。利用SPM軟件包中提供的V=spm_vol([pathname filename])函數(shù)，將預(yù)處理后的fMRI數(shù)據(jù)讀入Matlab。

應(yīng)用ICA-TSC方法，對fMRI數(shù)據(jù)逐點(diǎn)進(jìn)行處理， 將每一層fMRI數(shù)據(jù)轉(zhuǎn)換成由相關(guān)系數(shù)構(gòu)成的統(tǒng)計(jì)參數(shù)，并進(jìn)行Z變換再將其轉(zhuǎn)換為Z分布，圖8所示為受試者第1～40層Z變換前后的處理結(jié)果。

圖8 相關(guān)系數(shù)分布轉(zhuǎn)換為Z分布。(a)(b)表示相關(guān)系數(shù)的分布情況及其灰度直方圖。(c)(d)表示Z分布及其灰度直方圖。Fig. 8 The distribution of the correlation coefficient is converted to Z distribution. (a)(b) Shows the distribution of the correlation coefficient and the gray histogram. (c)(d) Shows the distribution of Z and the gray histogram.

Z變換后，拉大了激活區(qū)域與非激活區(qū)域之間的差距，提高了激活區(qū)提取的靈敏度。根據(jù)檢驗(yàn)的顯著性水平(α=0.05)閾值化Z參數(shù)圖，以判別體素是否激發(fā)，并將激發(fā)體素疊加在腦結(jié)構(gòu)圖像上。

3 結(jié)果

在仿真數(shù)據(jù)實(shí)驗(yàn)中，圖4為利用ICA-TSC方法得到的統(tǒng)計(jì)參數(shù)，各點(diǎn)的灰度值與自相關(guān)系數(shù)成正比。可以看出，對應(yīng)的代表激活信號的A區(qū)域由于沒有噪聲影響，相關(guān)性很高；包含3種混合信號的D區(qū)域，相關(guān)性較弱；生理噪聲C區(qū)域也能表現(xiàn)出微弱的相關(guān)性；而系統(tǒng)噪聲B區(qū)域的相關(guān)性，則幾乎為零。實(shí)驗(yàn)結(jié)果證實(shí)，ICA-TSC方法可以成功提取出A和D區(qū)域(包含激活信息區(qū)域)，并剔除B和C區(qū)域(不包含激活信息區(qū)域)，結(jié)果如圖6所示。其中D區(qū)域分離后的激活信號如圖7所示。

在真實(shí)fMRI數(shù)據(jù)實(shí)驗(yàn)中，圖9為前4名受試者第40層的處理結(jié)果，左側(cè)所示的是所提出的ICA-TSC方法，右側(cè)的是作為對照的GLM方法。接著給出群體統(tǒng)計(jì)結(jié)果(見圖10)和時(shí)域解釋變量(見圖11(a))。應(yīng)用本文2.4節(jié)中描述的方法，得到兩種方法的準(zhǔn)確性指標(biāo)，如表1所示。

表1 兩種方法的準(zhǔn)確性指標(biāo)(均值±標(biāo)準(zhǔn)差)

空間準(zhǔn)確性檢驗(yàn)結(jié)果z=-0.650 3>k，與GLM方法沒有差別，所得結(jié)果在“取偽”概率P<0.01的約束條件下具有統(tǒng)計(jì)意義，證明本ICA-TSC方法的空間提取能力可信度極高。從圖9中可以看出，ICA-TSC方法所得的結(jié)果與GLM方法非常相似。

時(shí)間準(zhǔn)確性檢驗(yàn)結(jié)果t=-6.240 4，遠(yuǎn)小于k=-t0.01=-2.508 3，證明ICA-TSC方法的時(shí)間準(zhǔn)確性遠(yuǎn)高于GLM的時(shí)間準(zhǔn)確性，結(jié)果具有極顯著性差異(“棄真”概率P<0.01)。從圖11中可以看出，與 (b)相比，(c)與(a)的相似度更高。

圖9 前4名受試者的第40層處理結(jié)果。(a)ICA-TSC方法。(b)GLM方法(p<0.05)。Fig. 9 The 40th floor processing results of the first four subjects. (a) ICA-TSC method. (b) GLM方法 method (p<0.05).

圖10 SPM群體統(tǒng)計(jì)結(jié)果的玻璃(透明)腦顯示(p<0.005), (a)(b)(c)分別表示玻璃腦的矢狀面、冠狀面和橫截面。將橫截面的第40層作為檢驗(yàn)的參照標(biāo)準(zhǔn)。Fig. 10 A glass (transparent) brain to display the groups statistics results of the SPM, (a)(b)(c) respectively represent sagittal, coronal and cross-section of the glass brain. The 40th floor of the cross-section was the test reference standard.

圖11 時(shí)域解釋變量與兩種方法所得最大峰值點(diǎn)時(shí)間序列的對比。(a)時(shí)域解釋變量；(b)GLM方法所得最相關(guān)時(shí)間序列rmax1=0.434 0； (c) ICA-TSC方法所得最相關(guān)時(shí)間序列rmax2=0.793 5。Fig. 11 Comparison of time domain regressors and the most relevant time series obtained by the two methods. (a)Time domain regressors; (b)The most relevant time series for GLM method rmax1=0.434 0; (c)The most relevant time series for ICA-TSC method rmax2=0.793 5.

4 討論和結(jié)論

從仿真的實(shí)驗(yàn)結(jié)果可以清楚地看出， ICA-TSC算法能夠準(zhǔn)確地提取含有激活成分的區(qū)域。TSC方法由于是通過檢測體素點(diǎn)各周期的時(shí)間序列相關(guān)性，對fMRI數(shù)據(jù)進(jìn)行激活區(qū)提取，因此僅適用于任務(wù)實(shí)驗(yàn)的組塊設(shè)計(jì)，并不適用于事件相關(guān)設(shè)計(jì)和靜息態(tài)。時(shí)間自相關(guān)方法相對參考函數(shù)法更為敏感，體素時(shí)間序列只要滿足特定長度的周期性，就可以被檢測出來。因此，由于血流動力學(xué)響應(yīng)延遲造成的影響就可以被避免，提高了激活區(qū)檢測的靈敏度。Infomax算法中初始值選取的不同，會對結(jié)果造成輕微影響，是算法不穩(wěn)定的因素之一。另外，步長μ也是影響結(jié)果好壞的重要因素。TSC方法中，臨界值的選擇非常關(guān)鍵，應(yīng)用Z變換，將0～1分布的相關(guān)系數(shù)T轉(zhuǎn)換成標(biāo)準(zhǔn)正態(tài)分布N(0,1)的Z系數(shù)，便于應(yīng)用假設(shè)檢驗(yàn)的方法自適應(yīng)地確定閾值，準(zhǔn)確地檢測及定位腦功能激活區(qū)。雙手握拳運(yùn)動的激活區(qū)主要分布于雙側(cè)SI(中央后回, BA 1 2 3), MI(中央前回, BA 4), SMA(額內(nèi)側(cè)回, BA 6)。由于GLM方法是經(jīng)典的模型驅(qū)動方法，具有較高的可信度，基于數(shù)據(jù)驅(qū)動的方法很難超越，因此一般情況下，要檢驗(yàn)新方法與GLM方法的可信度是否相似。實(shí)驗(yàn)結(jié)果表明，ICA-TSC方法與GLM方法的空間結(jié)果相似度極高。在時(shí)間序列上，進(jìn)行了ICA提取后，ICA-TSC方法的周期性要遠(yuǎn)高于GLM方法的周期性。

綜上所述，時(shí)間自相關(guān)方法(TSC)與ICA方法結(jié)合，在不需要參考信號的情況下可以靈活并準(zhǔn)確地對fMRI數(shù)據(jù)進(jìn)行激活區(qū)提取，具有不輸于GLM方法的空間準(zhǔn)確性和優(yōu)于GLM方法的時(shí)間準(zhǔn)確性。

[1] 楊福生, 洪波. 獨(dú)立分量分析的原理與應(yīng)用[M]. 北京：清華大學(xué)出版社, 2006: 184-195.

[2] Jutten C, Harault J. Blind separation of sources, part I:an adaptive algorithm based on neuromimetic architecture[J]. Signal Processing, 1991,24(1):1-10.

[3] 潘麗麗, 史振威, 唐煥文,等. fMRI 信號盲分離的一種獨(dú)立成分分析算法[J]. 大連理工大學(xué)學(xué)報(bào), 2005,45(4): 407-611.

[4] Esposito F, Formisano E, Seifritz E,etal. Spatial independent component analysis of functional MRI time-series: to what extent do results depend on the algorithm used[J]. Human Brain Mapping, 2002, 16: 146-157.

[5] 韓彤, 崔世民, 劉梅麗, 等. 利用fMRI和雙手交替運(yùn)動模式研究腦腫瘤所致的運(yùn)動功能重組[J]. 中國影像技術(shù), 2005, 9: 1324-1328.

[6] McKeown MJ, Makeig S, Brown GG,etal. Analysis of fMRI data by blind separation into independent spatial components [J]. Human Brain Mapping, 1998, 6: 160-188.

[7] Li Yu, Liu Jingsen. Mechanism and Improvement of Direct Anonymous Attestation Scheme [J]. Journal of Henan University, 2007,37(2): 195-197.

[8] Gossl C, Fahrmeir L. Bayesian Modeling of the Hemodynamic Response Function in BOLD fMRI [J]. Neuro Image, 2001, 14: 140-148.

[9] 吳義根, 李可. SPM軟件包數(shù)據(jù)處理原理簡介——第一部分:基本數(shù)學(xué)原理[J]. 中國醫(yī)學(xué)影像技術(shù), 2004,11:1768-1772.

[10] McKeown MJ. Whole brain function MR imaging activation from a finger-tapping task examined with independent component analysis [J]. American J Neuroradiol, 2000, 21:1629-1635.

[11] Calhoun VD, Adali T, Pearlson GD,etal. A method for making group inference from functional MRI data using independent component analysis [J]. Hum Brain Map, 2001, 14:140-151.

[12] Calhoun VD, Adali T, Pearlson GD,etal. Spatial and temporal independent component analysis of function MRI data containing a pair of task-relaed waveforms[J]. Hum Brain Map, 2001, 13:43-53.

[13] 王世杰，王立功，羅立民. 自適應(yīng)閾值化功能磁共振統(tǒng)計(jì)參數(shù)映射圖[J]. 中國圖像學(xué)報(bào), 2003, 10: 1196-2000.

ApplicationoffMRIActivationAreaAdaptiveExtractionBasedonTemporalSelf-CorrelationMethodofICA

BAI Bao-Dong*LIU Jian GUO Hong-Yu

(InstituteofBiomedicalandElectromagneticEngineering,ShenyangUniversityofTechnology,Shenyang110870,China)

The traditional ICA based activation zone detection is to analyze the correlation between the separated independent component and the reference signal. However, in practical problems, since the differences among the hemodynamic responses of the cerebral regions, the standard reference signal is often not available. Aiming at such problems, in this paper, the method of temporal self-correlation (TSC) combined with infomax-ICA was proposed. This method processed fMRI data point by point with 5-adjacent voxels based ICA, then detected the correlation between each time series period with temporal self-correlation algorithm and selected the maximum autocorrelation coefficient as the signal value of the voxel. After that we conversed correlation coefficient distribution to Z distribution which obeyN(0, 1) by Z-transform, extracted the active regions with significant difference (a=0.05) in the statistical parametric mapping. The algorithm was applied to deal with simulation data and 12 set of real fMRI data of fist movement with both hands. Results: The method can accurately extract the active region of the simulation data. For real data processing, results of this method have a high similarity with GLM method in the spatial domain (0.465 3±0.136 8 vs 0.490 5±0.134 1) and better than GLM method in the temporal domain (0.636 4±0.011 1 vs 0.369 2±0.010 9). These results have statistical significant. Experimental results showed that, this method has good capacities of detection of functional brain activation areas and spatial orientation.

independent component analysis (ICA); fMRI; temporal self-correlation; hypothesis test; adaptive threshold value

10.3969/j.issn.0258-8021. 2014. 02.009

2012-11-10， 錄用日期:2013-12-01

遼寧省自然科學(xué)基金(201102163)；國家自然科學(xué)基金(50907041)；遼寧省教育廳科研項(xiàng)目(201134120)

R318

A

0258-8021(2014) 02-0194-08

*通信作者。E-mail: baibd@sut.edu.cn