專利名稱：：基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法
技術(shù)領(lǐng)域：
：本發(fā)明涉及功能磁共振成像技術(shù)中數(shù)據(jù)處理領(lǐng)域，尤其涉及一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法。
背景技術(shù)：
：基于血氧水平依賴(BloodOxygenationLevelDependent，BOLD)對比度機(jī)制，功能磁共振成像(functionalMRI，fMRI)能夠?qū)崟r地檢測與大腦皮層神經(jīng)活動相關(guān)聯(lián)的局部血流動力學(xué)改變。但是，基于BOLD效應(yīng)的磁共振信號改變非常小，在磁場強(qiáng)度為1.5T時只有0.5～2％。此外，在序列圖像采集過程中周期呼吸和心臟運(yùn)動導(dǎo)致的腦組織脈動以及其它一些復(fù)雜的生理運(yùn)動都會導(dǎo)致許多生理相關(guān)的MRI信號改變。這些生理相關(guān)的MRI信號改變構(gòu)成了功能磁共振數(shù)據(jù)的生理噪聲。除了與心臟和呼吸相關(guān)的生理噪聲之外，功能磁共振數(shù)據(jù)中還包含源于射頻線圈和受檢者的隨機(jī)噪聲以及其它低頻噪聲，這些低頻噪聲主要是由于掃描器的不穩(wěn)定、運(yùn)動校正后殘留的運(yùn)動關(guān)聯(lián)效應(yīng)、長期的生理漂移等因素造成的。功能磁共振數(shù)據(jù)中這些噪聲成份與神經(jīng)活動無關(guān)，不僅降低了對功能磁共振數(shù)據(jù)統(tǒng)計分析的顯著性，而且直接影響功能磁共振的敏感性及空間定位的可靠性。功能磁共振數(shù)據(jù)中的噪聲主要是由于與心臟和呼吸相關(guān)的非神經(jīng)性的生理效應(yīng)引起的。除了與心臟和呼吸相關(guān)的生理噪聲之外，功能磁共振數(shù)據(jù)中還包含源于射頻線圈和受檢者的隨機(jī)噪聲以及其它低頻噪聲，這些低頻噪聲主要是由于掃描器的不穩(wěn)定、運(yùn)動校正后殘留的運(yùn)動關(guān)聯(lián)效應(yīng)、長期的生理漂移等因素造成的。因此，功能磁共振數(shù)據(jù)模型可以表示為yi(t)＝si(t)+pi(t)+ri(t)(1)其中，yi(t)表示每個象素的時間過程，si(t)表示真正與神經(jīng)活動相關(guān)的功能信號改變，pi(t)表示非神經(jīng)性的生理效應(yīng)引起的生理的噪聲及各種低頻噪聲，ri(t)表示系統(tǒng)的隨機(jī)噪聲。為了減少由于生理運(yùn)動導(dǎo)致的功能磁共振信號改變，已經(jīng)提出了許多有效的方法，如采用心電門控掃描、導(dǎo)航回波、數(shù)字濾波、回溯性校正及基于k-空間的估計和校正技術(shù)等。這些方法或者要求實(shí)時采集或監(jiān)測生理運(yùn)動信息，或者要求特定的掃描序列，都是以增加數(shù)據(jù)采集負(fù)擔(dān)或數(shù)據(jù)后處理復(fù)雜性為代價的。此外，為了捕捉生理運(yùn)動周期，這些方法要求功能磁共振數(shù)據(jù)的采集頻率高于生理運(yùn)動頻率。一些學(xué)者最近提出采用成分分析方法對功能磁共振數(shù)據(jù)進(jìn)行成分分解，然后通過對分離成分進(jìn)行頻譜估計確定并消除相應(yīng)生理噪聲成分。雖然這種方法不需要額外采集生理運(yùn)動信息，數(shù)據(jù)后處理較為簡單。但是，為了估計并確定生理噪聲成分，該方法仍然要求功能數(shù)據(jù)的采集頻率高于生理運(yùn)動頻率。為了覆蓋足夠的腦空間，功能磁共振常采用多層掃描方式采集數(shù)據(jù)。由于時間序列圖象的采樣率低于生理運(yùn)動頻率，導(dǎo)致高頻生理噪聲混疊到低頻成分，生理噪聲不再保持原來的波形，無法基于頻譜分析方法確定多層功能磁共振數(shù)據(jù)的生理噪聲成分。
發(fā)明內(nèi)容本發(fā)明提供了一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，可以有效地消除多層功能磁共振數(shù)據(jù)中生理噪聲，解決現(xiàn)有的消除相應(yīng)生理噪聲成分必須以增加數(shù)據(jù)后處理復(fù)雜性為代價或者必須要求功能數(shù)據(jù)的采集頻率高于生理運(yùn)動頻率等苛刻條件的問題。本發(fā)明采用的技術(shù)方案是一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，采用下列步驟1)設(shè)功能磁共振體數(shù)據(jù)序列的時間采樣點(diǎn)數(shù)為N，每個體數(shù)據(jù)包含S層圖像，分別以第i層圖像的灰質(zhì)區(qū)域和腦脊液區(qū)域?yàn)槟０?，?gòu)造該層時間序列圖像對應(yīng)的灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc，矩陣Yg和Yc的行表達(dá)了像素的時間過程；2)對灰質(zhì)數(shù)據(jù)矩陣Yg進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前p個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取灰質(zhì)數(shù)據(jù)矩陣Yg的主成份分量；3)對腦脊液數(shù)據(jù)矩陣Yc進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前q個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取腦脊液數(shù)據(jù)矩陣Yc的主成份分量；4)采用傅立葉頻譜分析方法分別對灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc的主成分分量進(jìn)行頻譜分析，根據(jù)隨機(jī)噪聲判定準(zhǔn)則確定并消除灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc主成分分量中的隨機(jī)噪聲成分；5)分別利用步驟4)保留的主成分分量重建消除隨機(jī)噪聲后的灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)6)對灰質(zhì)數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wg及m個統(tǒng)計上相互獨(dú)立的成份cig(t)，i＝1，2，Λm；7)對腦脊液數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wc及n個統(tǒng)計上相互獨(dú)立的成份cjc(t)，j＝1，2，Λn；8)利用灰質(zhì)數(shù)據(jù)的m個獨(dú)立成份cig(t)構(gòu)造矩陣Xm×N，利用腦脊液數(shù)據(jù)的n個獨(dú)立成份cjc(t)構(gòu)造矩陣Yn×N；9)對矩陣Xm×N和Yn×N進(jìn)行典型相關(guān)分析，得到cig(t)的加權(quán)系數(shù)矩陣Am×r、cjc(t)加權(quán)系數(shù)矩陣Bn×r及包含r個元素的行向量r為m和n中較小的一個數(shù)，的第k個元素為矩陣Xm×N和Yn×N的各行分別以矩陣Am×r和Bn×r的第k列元素為加權(quán)系數(shù)加權(quán)后得到的典型變量的典型相關(guān)系數(shù)的平方，k＝1，2，Λr；10)計算行向量的各元素對矩陣Am×r各行元素的加權(quán)和，得到表達(dá)灰質(zhì)數(shù)據(jù)中各獨(dú)立成份cig(t)與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)測度的列向量11)對列向量的各元素從大到小排序，那么，排序后的列向量的前r個元素對應(yīng)的灰質(zhì)數(shù)據(jù)的獨(dú)立成分cig(t)為與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的噪聲成分，將這些成分置零后得到一組新的獨(dú)立成份12)用分離矩陣Wg的逆矩陣乘以由獨(dú)立成份構(gòu)成的成份矩陣，重建消除與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的各種噪聲成份后的灰質(zhì)數(shù)據(jù)；14)重復(fù)步驟1)到12)，直至消除功能磁共振時間序列數(shù)據(jù)各層灰質(zhì)數(shù)據(jù)中的與腦脊液數(shù)據(jù)中相關(guān)最大的各種噪聲成分。與現(xiàn)有的技術(shù)相比，本發(fā)明具有如下優(yōu)點(diǎn)本發(fā)明提出基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法基于生理噪聲成分時間過程在大腦灰質(zhì)區(qū)域以及腦脊液區(qū)域的空間相關(guān)性，通過對腦脊液區(qū)域和灰質(zhì)區(qū)域數(shù)據(jù)進(jìn)行獨(dú)立成份分解，利用典型相關(guān)分析方法分析腦脊液區(qū)域和灰質(zhì)區(qū)域獨(dú)立成份時間過程的空間相關(guān)性，能夠識別出灰質(zhì)區(qū)域數(shù)據(jù)中與腦脊液區(qū)域數(shù)據(jù)中相關(guān)最大的噪聲成份，本發(fā)明提出的方法不僅能夠有效地消除多層功能磁共振數(shù)據(jù)中生理噪聲，而且能夠消除與腦脊液區(qū)域數(shù)據(jù)中相關(guān)最大的其它低頻噪聲成份，較好地實(shí)現(xiàn)了多層功能磁共振數(shù)據(jù)中噪聲的消除。由于本發(fā)明提出的方法數(shù)據(jù)處理簡單，不需要實(shí)時采集或監(jiān)測生理運(yùn)動信息，可以對早期采集的功能磁共振數(shù)據(jù)進(jìn)行離線處理，因此是一種較為適用的功能磁共振數(shù)據(jù)處理方法，非常適于臨床醫(yī)院及研究單位利用功能磁共振進(jìn)行腦科學(xué)研究。圖1為聽覺功能實(shí)驗(yàn)刺激時序2-1為相關(guān)系數(shù)閾值T＝0.60時對實(shí)驗(yàn)數(shù)據(jù)的第30層圖像直接進(jìn)行相關(guān)分析檢測到的可能的激活像素圖2-2為相關(guān)系數(shù)閾值T＝0.25時對第30層圖像直接進(jìn)行相關(guān)分析檢測到的可能的激活像素圖2-3為采用本發(fā)明提出的方法消除灰質(zhì)數(shù)據(jù)中生理相關(guān)的各種噪聲后，相關(guān)系數(shù)閾值T＝0.60時檢測到的激活像素圖2-4為將圖2-3的檢測結(jié)果映射到實(shí)驗(yàn)數(shù)據(jù)第30層圖像上的結(jié)果圖3-1為灰質(zhì)數(shù)據(jù)中的10個獨(dú)立成份圖3-2為腦脊液數(shù)據(jù)中的5個獨(dú)立成份圖4-1為某一激活像素消除噪聲之前的時間過程圖4-2為某一激活像素消除噪聲之后的時間過程圖5-1為消除噪聲前后激活像素時間過程與功能刺激參考函數(shù)的相關(guān)系數(shù)比較圖5-2為消除噪聲前后激活像素時間過程的方差比較圖6-1到圖6-6為采用本發(fā)明提出的噪聲消除方法對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理之后，采用相關(guān)分析法得到的從第28到第33層連續(xù)6層圖像的聽覺功能激活映射圖。具體實(shí)施例方式一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，采用下列步驟1)設(shè)功能磁共振體數(shù)據(jù)序列的時間采樣點(diǎn)數(shù)為N，每個體數(shù)據(jù)包含S層圖像，分別以第i層圖像的灰質(zhì)區(qū)域和腦脊液區(qū)域?yàn)槟０?，?gòu)造該層時間序列圖像對應(yīng)的灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc，矩陣Yg和Yc的行表達(dá)了像素的時間過程；2)對灰質(zhì)數(shù)據(jù)矩陣Yg進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前p個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取灰質(zhì)數(shù)據(jù)矩陣Yg的主成份分量；3)對腦脊液數(shù)據(jù)矩陣Yc進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前q個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取腦脊液數(shù)據(jù)矩陣Yc的主成份分量；4)采用傅立葉頻譜分析方法分別對灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc的主成分分量進(jìn)行頻譜分析，根據(jù)隨機(jī)噪聲判定準(zhǔn)則確定并消除灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc主成分分量中的隨機(jī)噪聲成分；5)分別利用步驟4)保留的主成分分量重建消除隨機(jī)噪聲后的灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)6)對灰質(zhì)數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wg及m個統(tǒng)計上相互獨(dú)立的成份cig(t)，i＝1，2，Λm；7)對腦脊液數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wc及n個統(tǒng)計上相互獨(dú)立的成份cjc(t)，j＝1，2，Λn；8)利用灰質(zhì)數(shù)據(jù)的m個獨(dú)立成份cig(t)構(gòu)造矩陣Xm×N，利用腦脊液數(shù)據(jù)的n個獨(dú)立成份cjc(t)構(gòu)造矩陣Yn×N；9)對矩陣Xm×N和Yn×N進(jìn)行典型相關(guān)分析，得到cig(t)的加權(quán)系數(shù)矩陣Am×r、cjc(t)加權(quán)系數(shù)矩陣Bn×r及包含r個元素的行向量r為m和n中較小的一個數(shù)，的第k個元素為矩陣Xm×N和Yn×N的各行分別以矩陣Am×r和Bn×r的第k列元素為加權(quán)系數(shù)加權(quán)后得到的典型變量的典型相關(guān)系數(shù)的平方，k＝1，2，Λr；10)計算行向量的各元素對矩陣Am×r各行元素的加權(quán)和，得到表達(dá)灰質(zhì)數(shù)據(jù)中各獨(dú)立成份cig(t)與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)測度的列向量11)對列向量的各元素從大到小排序，那么，排序后的列向量的前r個元素對應(yīng)的灰質(zhì)數(shù)據(jù)的獨(dú)立成分cig(t)為與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的噪聲成分，將這些成分置零后得到一組新的獨(dú)立成份12)用分離矩陣Wg的逆矩陣乘以由獨(dú)立成份構(gòu)成的成份矩陣，重建消除與腦脊液數(shù)據(jù)中各獨(dú)立成份ccj(t)相關(guān)最大的各種噪聲成份后的灰質(zhì)數(shù)據(jù)；15)重復(fù)步驟1)到12)，直至消除功能磁共振時間序列數(shù)據(jù)各層灰質(zhì)數(shù)據(jù)中的與腦脊液數(shù)據(jù)中相關(guān)最大的各種噪聲成分。本實(shí)施例還可以通過計算功能實(shí)驗(yàn)設(shè)計時序與上述步驟1)所述灰質(zhì)數(shù)據(jù)矩陣Yg中每行所表達(dá)像素的時間過程的相關(guān)系數(shù)，用相關(guān)系數(shù)大于閾值T的像素重新構(gòu)造灰質(zhì)數(shù)據(jù)的矩陣。上述隨機(jī)噪聲判定準(zhǔn)則是如果某一成分傅立葉頻譜的平均功率大于或等于頻譜功率的標(biāo)準(zhǔn)偏差，則認(rèn)為該成分為隨機(jī)噪聲成分。我們通過對實(shí)際的聽覺功能磁共振數(shù)據(jù)進(jìn)行處理與分析來驗(yàn)證提出的方法。實(shí)驗(yàn)數(shù)據(jù)由倫敦大學(xué)Wellcom神經(jīng)成像實(shí)驗(yàn)室提供。聽覺實(shí)驗(yàn)包含靜息和激發(fā)兩種狀態(tài)，靜息期不施加任何刺激，激發(fā)期以每分鐘60個雙音節(jié)英文單詞的速度施加聽覺刺激。這樣的實(shí)驗(yàn)共重復(fù)8次，每次實(shí)驗(yàn)中，對靜息期和激發(fā)期分別進(jìn)行6次采樣，每次采樣得到64幅連續(xù)的橫軸位斷層圖象，掃描矩陣為64×64×64，單次掃描采樣重復(fù)時間TR＝7s，整個實(shí)驗(yàn)持續(xù)6分鐘。具體實(shí)施方式如下1)利用目前較為通用的功能磁共振統(tǒng)計分析軟件SPM對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行運(yùn)動校正、空間標(biāo)準(zhǔn)化等預(yù)處理后，得到包含96個79×95×68體數(shù)據(jù)的數(shù)據(jù)集，并對數(shù)據(jù)集中第一個體數(shù)據(jù)進(jìn)行圖象分割，得到灰質(zhì)(GM)、白質(zhì)(WM)及腦脊液(CSF)，并以分割的灰質(zhì)和腦脊液圖像作為模板，分別構(gòu)造灰質(zhì)數(shù)據(jù)集和腦脊液數(shù)據(jù)集。在實(shí)驗(yàn)中，我們首先以聽覺功能實(shí)驗(yàn)的刺激時序作為參考函數(shù)ref(t)，如圖1所示，利用(7)計算ref(t)與灰質(zhì)數(shù)據(jù)集中每一像素時間過程yig(t)的相關(guān)系數(shù)。d=Σ(yig(t)-a‾i)(ref(t)-b‾)Σ(yig(t)-a‾i)2Σ(ref(t)-b‾)2---(7)]]>其中ai、b分別表示灰質(zhì)數(shù)據(jù)中的第i個像素時間過程yig(t)的平均值以及參考函數(shù)ref(t)的時間過程的平均值。如果相關(guān)系數(shù)大于設(shè)定的閾值T，則認(rèn)為該象素為可能的激活象素。然后由這些可能的激活像素重新構(gòu)造灰質(zhì)數(shù)據(jù)集。經(jīng)過這樣的初篩，處理的數(shù)據(jù)大量減少，一方面數(shù)據(jù)處理效率相應(yīng)提高，另一方面不僅保證了灰質(zhì)數(shù)據(jù)集的主成分中包含真正與神經(jīng)活動相關(guān)的功能信號成份，而且對灰質(zhì)數(shù)據(jù)集獨(dú)立成分估計的準(zhǔn)確度也提高。相關(guān)系數(shù)閾值T可以根據(jù)統(tǒng)計檢驗(yàn)的顯著性水平確定。這里，我們選取T＝0.60，相應(yīng)的Z分位數(shù)大于6，統(tǒng)計檢驗(yàn)的顯著性水平α小于0.01。如圖2-1中所示高亮區(qū)域顯示了閾值T＝0.60時實(shí)驗(yàn)數(shù)據(jù)第30層圖像中27個可能的激活象素。由于噪聲的影響，直接利用相關(guān)分析方法只檢測到少數(shù)可能的激活象素。為了增加可能的激活象素的樣本，我們降低閾值T。圖2-2中高亮區(qū)域顯示了閾值T＝0.25時實(shí)驗(yàn)數(shù)據(jù)第30層圖像中262個可能的激活象素。比較圖2-1和圖2-2，由于相關(guān)系數(shù)閾值降低，檢測到的可能的激活象素明顯增多，但是偽激活象素也明顯增加。這種結(jié)果正是由于功能磁共振數(shù)據(jù)中包含各種噪聲成份導(dǎo)致的。功能磁共振空間定位的可靠性降低。所以，為了消除功能磁共振數(shù)據(jù)中包含的各種噪聲成份，我們用圖2-2中262個可能的激活象素重新構(gòu)造灰質(zhì)數(shù)據(jù)集。2)對重新構(gòu)造的灰質(zhì)數(shù)據(jù)集和腦脊液數(shù)據(jù)集分別采用奇異分解算法(SVD)進(jìn)行主成分分析，并根據(jù)保留的特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取主成份分量，對灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)的主成份利用傅立葉變換進(jìn)行頻譜分析，如果某一成分整個頻譜的平均功率大于或等于整個頻譜功率的標(biāo)準(zhǔn)偏差，則認(rèn)為該成份為系統(tǒng)隨機(jī)噪聲成份，并將該成份置零，消除隨機(jī)噪聲成份。經(jīng)過這樣的處理，灰質(zhì)數(shù)據(jù)集中保留10個主成分，腦脊液數(shù)據(jù)集中保留5個主成分。3)分別利用步驟2)得到的主成分分量重建灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)通過計算功能實(shí)驗(yàn)設(shè)計時序與每一圖象灰質(zhì)區(qū)域象素的時間過程的相關(guān)系數(shù)，用相關(guān)系數(shù)大于閾值T的像素重新構(gòu)造灰質(zhì)數(shù)據(jù)的矩陣。4)對消除隨機(jī)噪聲后的灰質(zhì)數(shù)據(jù)采用Hyvarinen的快速定點(diǎn)獨(dú)立成分分析算法(fastICA)進(jìn)行獨(dú)立成分分解，結(jié)果如圖3-1所示。5)對消除隨機(jī)噪聲后的腦脊液數(shù)據(jù)分別采用與4)同樣的算法進(jìn)行獨(dú)立成分分解，結(jié)果如圖3-2所示。6)為了確定灰質(zhì)數(shù)據(jù)中的生理相關(guān)的噪聲成分，我們分別利用圖3-1、圖3-2所示的獨(dú)立成份構(gòu)造一個10×96的矩陣X和一個5×96矩陣Y。7)對矩陣X和Y進(jìn)行典型相關(guān)分析，可得到灰質(zhì)數(shù)據(jù)中獨(dú)立成份cig(t)的加權(quán)系數(shù)矩陣A5×5、腦脊液數(shù)據(jù)中獨(dú)立成份cjc(t)加權(quán)系數(shù)矩陣B10×5及包含5個元素的行向量行向量的第k個元素為矩陣X和Y的各行分別以矩陣Am×r和Bn×r的第k列元素加權(quán)后的典型相關(guān)系數(shù)的平方。8)計算行向量的各元素對矩陣Am×r各行元素的加權(quán)和，得到表達(dá)灰質(zhì)數(shù)據(jù)中各獨(dú)立成份cig(t)與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)的相關(guān)測度的列向量對列向量的各元素從大到小排序，結(jié)果如表1-1所示表1灰質(zhì)和腦脊液中獨(dú)立成份時間過程典型相關(guān)分析結(jié)果(第30層)<tablesid="table1"num="001"><tablewidth="745">M1.18330.96250.84890.82550.81670.73130.60810.53890.52380.2226IC17105648932</table></tables>從表1的排序結(jié)果可見，圖3-1所示的灰質(zhì)數(shù)據(jù)集中10個獨(dú)立成份與圖3-2所示的腦脊液數(shù)據(jù)集中的5個獨(dú)立成分相關(guān)的相關(guān)測度M從大到小一依次為第1、7、10、5、6、4、8、9、3、2。圖3-1所示的第2個獨(dú)立成分與圖3-2所示的腦脊液數(shù)據(jù)集中的5個獨(dú)立成分相關(guān)的相關(guān)測度最小，該獨(dú)立成分與圖2所示的聽覺功能實(shí)驗(yàn)刺激時序基本一致，反映了大腦灰質(zhì)區(qū)域?qū)β犛X功能刺激響應(yīng)的時間過程。9)根據(jù)前面介紹的典型相關(guān)分析方法，我們選取表1前5個獨(dú)立成分作為與腦脊液數(shù)據(jù)中5個獨(dú)立成份相關(guān)最大的生理相關(guān)的噪聲成分。將這些成分置零后，10)根據(jù)(8)重建灰質(zhì)數(shù)據(jù)集Yg~=W-1C0g,]]>其中，W-1表示分離矩陣W的逆矩陣，C0g表示將噪聲成份置零后的灰質(zhì)成份矩陣。圖4-1、圖4-2中實(shí)線分別表示同一激活像素在消除噪聲之前和消除噪聲之后的時間過程，虛線表示反映聽覺功能實(shí)驗(yàn)刺激時序的參考函數(shù)的時間過程?？梢姡肼曋?，激活像素時間過程的波形明顯改善，聽覺功能實(shí)驗(yàn)的靜息與激發(fā)兩種狀態(tài)對比也更加明顯。為了進(jìn)一步定量說明噪聲的消除效果，我們分別計算像素的時間過程在噪聲消除之前和之后與參考函數(shù)的相關(guān)系數(shù)以及像素的時間過程在噪聲消除之前和之后的方差。圖6-1中虛線顯示了重構(gòu)的灰質(zhì)數(shù)據(jù)集中262個可能的激活像素的時間過程在消除噪聲之前與參考函數(shù)的相關(guān)系數(shù)，實(shí)線顯示了消除噪聲之后與參考波形的相關(guān)系數(shù)，可見，消除噪聲之后相關(guān)系數(shù)明顯變大，之前相關(guān)系數(shù)最大值為0.6715，平均值為0.4261；之后相關(guān)系數(shù)最大值為0.7547，平均值為0.6097，相關(guān)系數(shù)平均值增加43.09％。6-2中虛線和實(shí)線分別顯示了重構(gòu)的灰質(zhì)數(shù)據(jù)集中262個像素的時間過程在消除噪聲之前和之后的均方差。可見，消除噪聲之后均方差明顯變小，之前均方差為5541.2，之后均方差為3110.6，均方差減小43.86％。圖2-3的高亮區(qū)域顯示了消除噪聲之后與參考函數(shù)的相關(guān)系數(shù)大于0.60的激活像素。與圖2-1相比，雖然選取同樣的相關(guān)系數(shù)閾值，但是在對功能磁共振數(shù)據(jù)統(tǒng)計分析的顯著性保持不變的情況下，檢測到的激活像素個數(shù)增加到160個，功能磁共振的敏感性提高。與圖2-2相比，偽激活像素個數(shù)也明顯減少，功能磁共振空間定位的可靠性提高。圖2-4為將圖2-3映射到實(shí)驗(yàn)數(shù)據(jù)集第30層圖像上的聽覺功能激活映射圖。表2為功能磁共振數(shù)據(jù)消除噪聲前后與參考函數(shù)相關(guān)系數(shù)及方差改變比較。Num_B、Num_A分別表示消除噪聲之前和之后各數(shù)據(jù)層可能的激活像素。消除噪聲之前的相關(guān)系數(shù)閾值為T＝0.25，消除噪聲之后的相關(guān)系數(shù)閾值為T＝0.60。dmax_B、dmax_A分別表示消除噪聲前、后激活像素的最大相關(guān)系數(shù)，dmean_B、dmean_A分別表示消除噪聲前、后激活像素的相關(guān)系數(shù)平均值。Δd％表示消除噪聲之后相關(guān)系數(shù)平均值改變百分比。Cmean_B、Cmean_A分別表示消除噪聲前、后激活像素的均方差，Δv％表示消除噪聲前后激活像素的均方差百分比改變。從表2可見，采用本文提出的方法消除噪聲功能磁共振數(shù)據(jù)的噪聲后，激活像素的時間過程與參考函數(shù)的最大相關(guān)系數(shù)、平均相關(guān)系數(shù)都一致增大，激活像素時間過程的方差都一致減小。檢測到的激活像素數(shù)量明顯增多，偽激活像素數(shù)量也明顯減少。表2功能磁共振數(shù)據(jù)消除噪聲前后與參考函數(shù)相關(guān)系數(shù)及方差改變比較圖6-1到圖6-6分別顯示了先采用提出的噪聲消除方法對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理之后，然后采用相關(guān)分析法得到的從第28到第33層連續(xù)6層圖像的聽覺功能激活映射圖。從檢測結(jié)果可見大腦左右兩側(cè)的上顳葉及大腦邊緣葉顯示較顯著的激活，激活區(qū)域檢測結(jié)果與經(jīng)驗(yàn)的聽覺中樞功能解剖區(qū)域基本相符；大腦左右兩側(cè)的語言溝回區(qū)也顯示少量的激活信號，這是由于該功能實(shí)驗(yàn)采用雙音節(jié)英文單詞作為聽覺刺激時，受檢者的語言理解區(qū)功能也有響應(yīng)的緣故，該結(jié)果進(jìn)一步說明了本文提出消除功能磁共振數(shù)據(jù)的噪聲方法的有效性及可靠性。本發(fā)明的原理隨機(jī)噪聲判定準(zhǔn)則是如果某一成分傅立葉頻譜的平均功率大于或等于頻譜功率的標(biāo)準(zhǔn)偏差，則認(rèn)為該成分為隨機(jī)噪聲成分。主成分分析(PrincipalComponentAnalysis，PCA)是一種基于二階統(tǒng)計量的信號分解技術(shù)，其目的是將信號分解為若干統(tǒng)計意義上相互正交(或不相關(guān))的成分。獨(dú)立成分分析(IndependentComponentAnalysis，ICA)是一種基于高階統(tǒng)計量一種信號分解技術(shù)，其目的是將信號分解為若干統(tǒng)計意義上相互獨(dú)立的成分。由于主成份分析能夠較好地分離并消除隨機(jī)噪聲，獨(dú)立成份分析能夠較好地分離并消除生理噪聲及各種低頻噪聲，因此可以先對功能磁共振數(shù)據(jù)進(jìn)行主成份分析，利用傅立葉變換對每個成份進(jìn)行頻譜分析，如果某一成分傅立葉頻譜的平均功率大于或等于頻譜功率的標(biāo)準(zhǔn)偏差，則認(rèn)為該成份為隨機(jī)噪聲成份，并將該成份置零，然后重建消除隨機(jī)噪聲成份功能磁共振數(shù)據(jù)，再對消除隨機(jī)噪聲的功能磁共振數(shù)據(jù)進(jìn)行獨(dú)立成分分析，以消除生理噪聲及各種低頻噪聲。假設(shè)采用主成份分析消除系統(tǒng)隨機(jī)噪聲以后，功能磁共振數(shù)據(jù)中的功能信號、生理噪聲以及各種低頻噪聲是線性混合的，功能磁共振數(shù)據(jù)模型可進(jìn)一步用矩陣表示為Y＝AC(2)其中，A表示未知的線性混合矩陣，C表示包含功能信號和各種噪聲成分的矩陣。為了估計功能磁共振數(shù)據(jù)中的功能信號和噪聲成分，可以尋求一個滿足(3)的分離矩陣WC＝WY(3)當(dāng)分離矩陣W為混合矩陣A的逆矩陣時，功能磁共振數(shù)據(jù)中的功能信號和噪聲成分將能夠被分離出來。分離矩陣W可以采用Hyvrinen提出的基于信息論的快速定點(diǎn)算法迭代估計。利用獨(dú)立成份分析消除功能磁共振數(shù)據(jù)中各種生理相關(guān)噪聲的關(guān)鍵是如何確定生理相關(guān)的各種噪聲成份。在通常情況下，人體的心臟脈動頻率為0.75～1.5Hz，呼吸頻率為0.1～0.5Hz。根據(jù)Nayquist采樣定理，當(dāng)功能磁共振時間序列圖象的采樣率高于2×1.5Hz時，通過對分離出的成分進(jìn)行頻譜分析，能夠確定心臟和呼吸相關(guān)的生理噪聲成分。但是，為了覆蓋足夠的腦空間，功能磁共振常采用多層掃描方式采集數(shù)據(jù)，掃描重復(fù)時間(TR)通常大于或等于2000ms。這樣，每一層時間序列圖象的采樣率小于或等于0.5Hz，低于人體生理運(yùn)動頻率，導(dǎo)致高頻生理噪聲混疊到低頻成分，生理噪聲不再保持原來的波形，無法基于頻譜分析方法確定多層功能磁共振數(shù)據(jù)中生理噪聲成分。功能磁共振數(shù)據(jù)中的噪聲主要是由心臟和呼吸相關(guān)的非神經(jīng)性的生理效應(yīng)引起的。雖然周期性的心臟脈動和呼吸引起的非神經(jīng)性的生理效應(yīng)對不同的腦區(qū)域影響不同，但是大腦灰質(zhì)區(qū)域以及腦脊液區(qū)域同時受周期的心臟和呼吸噪聲影響。由于腦脊液數(shù)據(jù)中不可能有與大腦皮層神經(jīng)活動相關(guān)聯(lián)的功能信號改變，可以認(rèn)為腦脊液中的信號只包含生理相關(guān)的噪聲、各種低頻噪聲及系統(tǒng)隨機(jī)噪聲。因此，我們可以基于生理噪聲成分時間過程在大腦灰質(zhì)區(qū)域以及腦脊液區(qū)域的空間相關(guān)性，通過腦脊液數(shù)據(jù)中生理噪聲成分估計并確定灰質(zhì)數(shù)據(jù)中生理噪聲成分，最終消除多層功能磁共振數(shù)據(jù)的生理噪聲。因此，采用主成份分析方法分別消除灰質(zhì)數(shù)據(jù)、腦脊液數(shù)據(jù)中的隨機(jī)噪聲后，腦脊液區(qū)域每個象素的時間過程yc(t)可以用(4)表示yc(t)＝p(t)(t＝1，2，Λ，N)(4)灰質(zhì)區(qū)域每個象素的時間過程yg(t)可以用(5)表示yg(t)＝s(t)+p(t)(t＝1，2，Λ，N)(5)采用獨(dú)立成份分析方法分別對(4)、(5)表示的灰質(zhì)數(shù)據(jù)Yg、腦脊液數(shù)據(jù)Yc進(jìn)行獨(dú)立成份分解。假設(shè)灰質(zhì)數(shù)據(jù)中包含m個統(tǒng)計上相互獨(dú)立的成分cig(t)，腦脊液數(shù)據(jù)中包含n個統(tǒng)計上相互獨(dú)立的成分cjc(t)。為了確定灰質(zhì)數(shù)據(jù)中的生理相關(guān)的噪聲成分，我們利用灰質(zhì)數(shù)據(jù)的m個獨(dú)立成份cig(t)構(gòu)造矩陣Xm×N，利用腦脊液數(shù)據(jù)的n個獨(dú)立成份cjc(t)構(gòu)造矩陣Ym×N。通過對矩陣Xm×N和Ym×N進(jìn)行典型相關(guān)分析，可得到cig(t)的加權(quán)系數(shù)矩陣Am×r和cjc(t)的加權(quán)系數(shù)矩陣Bn×r及包含r個元素的行向量r＝min(m，n)。的第k個元素為矩陣Xm×N和Ym×N的各行分別以矩陣Am×r和Bn×r的第k列元素為加權(quán)系數(shù)進(jìn)行加權(quán)后所得典型變量的典型相關(guān)系數(shù)的平方。顯然，典型相關(guān)系數(shù)的平方越大，腦脊液中各成份cig(t)與灰質(zhì)中各成分cjc(t)的相關(guān)性越大；灰質(zhì)中成分cig(t)的加權(quán)系數(shù)越大，該成分與腦脊液中成份cjc(t)的相關(guān)性也越大。因此，可以通過計算行向量的各元素對矩陣Am×r各行元素的加權(quán)和，得到表達(dá)灰質(zhì)數(shù)據(jù)中各獨(dú)立成份cig(t)與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)測度的列向量Mρ=A×λTρ---(6)]]>對列向量的各元素從大到小排序，那么，排序后的列向量的前r個元素對應(yīng)的獨(dú)立成分cig(t)為與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的噪聲成分，將這些成分置零，就可以消除灰質(zhì)數(shù)據(jù)中與腦脊液數(shù)據(jù)中各種噪聲成份相關(guān)最大的各種噪聲成份。對真實(shí)的功能磁共振數(shù)據(jù)試驗(yàn)結(jié)果證明了這種方法是很有效的。權(quán)利要求1.一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，其特征在于采用下列步驟1)設(shè)功能磁共振體數(shù)據(jù)序列的時間采樣點(diǎn)數(shù)為N，每個體數(shù)據(jù)包含S層圖像，分別以第i層圖像的灰質(zhì)區(qū)域和腦脊液區(qū)域?yàn)槟０澹瑯?gòu)造該層時間序列圖像對應(yīng)的灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc，矩陣Yg和Yc的行表達(dá)了像素的時間過程；2)對灰質(zhì)數(shù)據(jù)矩陣Yg進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前p個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取灰質(zhì)數(shù)據(jù)矩陣Yg的主成份分量；3)對腦脊液數(shù)據(jù)矩陣Yc進(jìn)行主成份分析，得到N個按降序排列的特征根及其相應(yīng)的特征向量，根據(jù)前q個特征值對應(yīng)的特征向量的累積方差貢獻(xiàn)率大于90％選取腦脊液數(shù)據(jù)矩陣Yc的主成份分量；4)采用傅立葉頻譜分析方法分別對灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc的主成分分量進(jìn)行頻譜分析，根據(jù)隨機(jī)噪聲判定準(zhǔn)則確定并消除灰質(zhì)數(shù)據(jù)矩陣Yg和腦脊液數(shù)據(jù)矩陣Yc主成分分量中的隨機(jī)噪聲成分；5)分別利用步驟4)保留的主成分分量重建消除隨機(jī)噪聲后的灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)6)對灰質(zhì)數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wg及m個統(tǒng)計上相互獨(dú)立的成份cig(t)，i＝1，2，Λm；7)對腦脊液數(shù)據(jù)進(jìn)行獨(dú)立成份分解，得到分離矩陣Wc及n個統(tǒng)計上相互獨(dú)立的成份cjc(t)，j＝1，2，Λn；8)利用灰質(zhì)數(shù)據(jù)的m個獨(dú)立成份cig(t)構(gòu)造矩陣Xm×N，利用腦脊液數(shù)據(jù)的n個獨(dú)立成份cjc(t)構(gòu)造矩陣Yn×N；9)對矩陣Xm×N和Yn×N進(jìn)行典型相關(guān)分析，得到cig(t)的加權(quán)系數(shù)矩陣Am×r、cjc(t)加權(quán)系數(shù)矩陣Bn×r及包含r個元素的行向量r為m和n中較小的一個數(shù)，的第k個元素為矩陣Xm×N和Yn×N的各行分別以矩陣Am×r和Bn×r的第k列元素為加權(quán)系數(shù)加權(quán)后得到的典型變量的典型相關(guān)系數(shù)的平方，k＝1，2，Λr；10)計算行向量的各元素對矩陣Am×r各行元素的加權(quán)和，得到表達(dá)灰質(zhì)數(shù)據(jù)中各獨(dú)立成份cig(t)與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)測度的列向量11)對列向量的各元素從大到小排序，那么，排序后的列向量的前r個元素對應(yīng)的灰質(zhì)數(shù)據(jù)的獨(dú)立成分cig(t)為與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的噪聲成分，將這些成分置零后得到一組新的獨(dú)立成份12)用分離矩陣Wg的逆矩陣乘以由獨(dú)立成份構(gòu)成的成份矩陣，重建消除與腦脊液數(shù)據(jù)中各獨(dú)立成份cjc(t)相關(guān)最大的各種噪聲成份后的灰質(zhì)數(shù)據(jù)；13)重復(fù)步驟1)到12)，直至消除功能磁共振時間序列數(shù)據(jù)各層灰質(zhì)數(shù)據(jù)中的與腦脊液數(shù)據(jù)中相關(guān)最大的各種噪聲成分。2.根據(jù)權(quán)利要求1所述的基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，其特征在于通過計算功能實(shí)驗(yàn)設(shè)計時序與上述步驟1)所述灰質(zhì)數(shù)據(jù)矩陣Yg中每行所表達(dá)像素的時間過程的相關(guān)系數(shù)，用相關(guān)系數(shù)大于閾值T的像素重新構(gòu)造灰質(zhì)數(shù)據(jù)的矩陣。3.根據(jù)權(quán)利要求1所述的基于獨(dú)立成份空間相關(guān)性消除多層功能磁共振數(shù)據(jù)數(shù)據(jù)噪聲的方法，其特征在于隨機(jī)噪聲判定準(zhǔn)則是如果某一成分傅立葉頻譜的平均功率大于或等于頻譜功率的標(biāo)準(zhǔn)偏差，則認(rèn)為該成分為隨機(jī)噪聲成分。全文摘要本發(fā)明公開了一種基于獨(dú)立成份空間相關(guān)性消除功能磁共振數(shù)據(jù)噪聲的方法，先對圖像的灰質(zhì)區(qū)域和腦脊液區(qū)域進(jìn)行圖像分割，分別對其進(jìn)行主成份分析、頻譜分析，確定并消除主成分分量中的隨機(jī)噪聲成分，重建消除隨機(jī)噪聲后的灰質(zhì)數(shù)據(jù)和腦脊液數(shù)據(jù)，對重建后數(shù)據(jù)進(jìn)行獨(dú)立成份分解，利用獨(dú)立成份分別構(gòu)造矩陣，再進(jìn)行典型相關(guān)分后，排序，將這些成分置零后得到一組新的獨(dú)立成份，重建消除各種噪聲成份后的數(shù)據(jù)，重復(fù)以上步驟直至消除各層灰質(zhì)數(shù)據(jù)與腦脊液數(shù)據(jù)中相關(guān)最大的各種噪聲成分。本發(fā)明能夠有效地消除多層功能磁共振數(shù)據(jù)中生理噪聲，能夠消除與腦脊液區(qū)域數(shù)據(jù)中相關(guān)最大的其它低頻噪聲成份，較好地實(shí)現(xiàn)了多層功能磁共振數(shù)據(jù)中噪聲的消除。文檔編號G01R33/48GK1775172SQ20051012261公開日2006年5月24日申請日期2005年11月29日優(yōu)先權(quán)日2005年11月29日發(fā)明者王世杰,羅立民,李松毅申請人:東南大學(xué)