基于小波變換的功能磁共振圖像時間序列分步去噪

李微微　梅雪　周宇

摘要：

功能磁共振圖像（fMRI）數(shù)據(jù)中反映大腦神經(jīng)活動的感興趣信號常受到結(jié)構(gòu)噪聲和隨機噪聲的影響。為消除上述噪聲對分析激活體素的影響，對經(jīng)過SPM標準預(yù)處理的體素時間序列進行Activelets小波變換，并在得到尺度系數(shù)及細節(jié)系數(shù)后，針對兩類噪聲的不同特點進行分步去噪。第一步，在受結(jié)構(gòu)噪聲影響的尺度系數(shù)上，選用獨立成分分（ICA）析去識別并消除結(jié)構(gòu)噪聲源；第二步，提出一種改進的空域相關(guān)去噪算法在細節(jié)系數(shù)上對信號進行處理。值得注意的是，該算法利用鄰域體素之間的相似性，判定所處位置的細節(jié)系數(shù)反映噪聲還是神經(jīng)活動。實驗結(jié)果表明，經(jīng)過這兩步處理的數(shù)據(jù)可有效消除噪聲的影響，其中框架位移減少了1.5mm，尖峰百分比減少了2%，此外由去噪后的信號獲得的腦激活圖中一些明顯的偽激活區(qū)得到抑制。

關(guān)鍵詞：

功能磁共振圖像；去噪；結(jié)構(gòu)噪聲；隨機噪聲；小波變換

中圖分類號：

TP391.4

文獻標志碼：A

Abstract：

The neural activity signal of interest is often influenced by structural noise and random noise in functional Magnetic Resonance Imaging （fMRI） data. In order to eliminate noise effects in the analysis of activate voxels， the time series of voxels preprocessed by Statistical Parametric Mapping （SPM） were transformed by Activelets wavelet. After getting scale coefficient and detail coefficient， the two kinds of noise denoised were eliminated separately according to their corresponding characteristics. Firstly， the Independent Component Analysis （ICA） was used to identify and eliminate the structural noise sources. Secondly， an improved algorithm for spatial correlation was presented on the detail coefficient. In particular， in the improved algorithm， the voxel similarity in the neighborhood was used to determine whether the detail coefficient reflected the noise or the neural activity. Experimental results show that the processing of data effectively eliminate the effect of noise； specifically， the frame displacement decreased by 1.5mm and the percentage of spikes decreased by 2%； in addition， the false activation regions are obviously restrained in the spatial map got by denoised signals.

英文關(guān)鍵詞Key words：

functional Magnetic Resonance Imaging （fMRI）； denoising； structural noise； random noise； wavelet transform

0引言

功能磁共振圖像（functional Magnetic Resonance Imaging， fMRI）兼顧空間及時間分辨率，且具有無入侵性的特點，為研究大腦在一些特定腦區(qū)的反應(yīng)提供了強有力的工具，進而可以幫助研究者更好地理解感覺、認知和運動等信息在人腦中的活動[1]。然而fMRI信號常常受到來源眾多的噪聲影響，導(dǎo)致大量感興趣信號的丟失。例如隨機噪聲不僅會降低信號的信噪比，而且會影響最終血氧依賴水平的敏感性與特異性；而由于機器過熱或者頭動產(chǎn)生的結(jié)構(gòu)噪聲會產(chǎn)生大量的偽激活區(qū)[2]。

傳統(tǒng)的數(shù)據(jù)處理流程包含如高斯濾波、空間平滑等去噪的處理步驟，但是這些去噪方法往往達不到理想的效果[3]。SPM（Statistical Parametric Mapping）團隊中Maldjian等[4]也提出了一種簡單的方法——混合效應(yīng)建模。雖然目前該方法得到了廣泛應(yīng)用，但該方法更加側(cè)重于去分析群體激活模式的信息而不是個體體素的響應(yīng)。就目前而言，研究者Feis等[5]最常用的數(shù)據(jù)驅(qū)動去噪算法主要是基于獨立成分分析（Independent Component Analysis， ICA），該算法通過線性分解fMRI數(shù)據(jù)得到空間獨立的各個成分。其在去噪的過程中最重要的是判定各個獨立成分是代表著噪聲還是代表著感興趣信號。此時去噪的效果取決于能否準確判定各個獨立成分，這并不是件容易的事，因為有時候各個成分中不僅包含結(jié)構(gòu)噪聲信號，而且包含神經(jīng)信號[6]。此外，由于獨立成分分析算法數(shù)學(xué)模型的局限性導(dǎo)致對隨機噪聲的估計效果不佳，僅能對結(jié)構(gòu)噪聲作出估計。因此如何針對結(jié)構(gòu)和隨機兩類噪聲設(shè)計不同的去噪策略，成為目前研究的熱點。

小波分解為解決上述問題提供了新的思路，在對數(shù)據(jù)進行小波處理時，小波函數(shù)有類似多尺度導(dǎo)數(shù)算子的特性，因此類似邊界等的奇異點可以用少量的系數(shù)表示為特征，這個特性對于信號去噪非常有用。而最近Fikret Isik Karahanoglu等[7]的研究也表明，fMRI信號經(jīng)過小波變換后結(jié)構(gòu)噪聲的能量主要集中在尺度系數(shù)上，而隨機噪聲主要集中在細節(jié)系數(shù)上。由此，使用小波分解可有效地處理這兩類噪聲。

本文首先使用Activelets小波基分解fMRI信號，并獲得信號的稀疏表示；然后分別對尺度系數(shù)及細節(jié)系數(shù)進行分步處理；最后重構(gòu)信號，通過SPM獲得腦激活圖。結(jié)果表明，本文方法無需先驗知識，針對經(jīng)過標準SPM數(shù)據(jù)預(yù)處理步驟后的數(shù)據(jù)，能夠根據(jù)數(shù)值特征自動消除噪聲影響，提高定位由fMRI數(shù)據(jù)獲得的腦激活區(qū)的準確性。

1本文算法

本文算法基于原始fMRI信號模型，在小波域?qū)υ糵MRI信號進行濾波，具體流程見圖1?？傮w而言，由于考慮到兩類噪聲在小波域的不同特點：結(jié)構(gòu)噪聲的能量主要集中在尺度系數(shù)上，而隨機噪聲主要集中在細節(jié)系數(shù)，本文通過選擇Activelets小波基分解體素時間序列后，分別對這兩類系數(shù)進行處理。

1.1基于Activelets小波基的fMRI信號分解

合適的小波基可以很好地將信號用小波系數(shù)稀疏表示，也決定著去噪的結(jié)果?；谘簞恿W(xué)為一個穩(wěn)定的線性系統(tǒng)的假設(shè)上，Khalidov 等[8]提出了Activelets小波基。本文使用Activelets小波基對fMRI時間序列進行三層分解與重構(gòu)，得到新的尺度系數(shù)和細節(jié)系數(shù)。值得注意的是，Activelets小波基ψ對應(yīng)的細節(jié)系數(shù)部分保證了血氧動力學(xué)中小波系數(shù)快速消退。具體而言，fMRI單體素時間序列可以寫成由尺度系數(shù)及細節(jié)系數(shù)構(gòu)成的如下形式：

x（t）=∑kcj（k）2j/2φ（2jt-k）+dj（k）2j/2ψ（2jt-k）（1）

其中，鑒于尺度空間的嵌套性，Vi+1∈Vi，其尺度函數(shù)φi+1（t）為：

φi+1（t）=∑khi[k]φi（t-2ik）（2）

其中：hi為尺度濾波器。同樣的，因為Vi⊕⊥Wi=Vi-1，Activelets小波函數(shù)ψi+1（t）為：

ψi+1（t）=∑kgi[k]φi（t-2ik）（3）

其中：hi及gi的濾波器組值的設(shè)定是根據(jù)血液動力學(xué)響應(yīng)函數(shù)演變而來，詳見文獻[9]。鑒于它們的濾波器性質(zhì)，其分解與重構(gòu)仍基于經(jīng)典的Mallat快速分解算法[10]。

1.2基于尺度系數(shù)的ICA分解與重構(gòu)系數(shù)

為消除結(jié)構(gòu)噪聲的能量，提高ICA分解信號的準確性，本文對單腦區(qū)分別使用ICA算法。假設(shè)單一腦區(qū)內(nèi)有k個體素，即有同時記錄的k個fMRI信號X（t）={x1（t），x2（t），…，xk（t）}，其是由N個未知的獨立源信號S（t）={s1（t），s2（t），…，sN（t）}線性疊加而成，公式如下：

X（t）=MS（t）（4）

其中：M為未知混合矩陣，ICA算法的目的就是從觀測信號X（t）中估計出混合矩陣M及源信號S（t）。因結(jié)構(gòu)噪聲的能量主要集中在尺度系數(shù)上，其單個源信號可有小波變化后的尺度系數(shù)表出，公式如下：

si=∑Kk=1ci（k）φk（5）

為簡化公式，現(xiàn)表示為如下形式：

si=Csi·Φ（6）

其中：Csi=（ci（k））Kk=1是具有k個小波系數(shù)的行向量；如果Φ選擇的是一組正交基，那么便存在源信號si和其小波系數(shù)Csi=si·Φ-1（Csi=si·ΦH）一一對應(yīng)的關(guān)系，因此其信號模型如式（7）。

S=Cs·Φ（7）

其中Cs為：

Cs=cs1cs2csk=[Cs（1），Cs（2），…，Cs（K）]（8）

因為尺度系數(shù)占據(jù)了信號的大部分能量，可用其近似原始信號，公式如下：

X≈CXΦ（9）

其中CX為尺度系數(shù)的小波系數(shù)矩陣。至此稀疏近似混合觀測信號矩陣X可用盲源分離方法的混合矩陣M及小波分解的尺度系數(shù)表出，公式如下：

CX≈MCs（10）

本文使用FastICA算法[11]由CX估計Cs，通過結(jié)構(gòu)噪聲特性從Cs中識別并消除其對數(shù)據(jù)的影響。結(jié)構(gòu)噪聲常表現(xiàn)出如下三個特性：1）高頻時間序列經(jīng)過傅里葉變換后有50%的能量高于0.1Hz；2）尖峰，有一個或者多個尖峰信號；3）鋸齒波模式，即劇烈且有規(guī)律地上下震蕩變化的時間序列[12]。而其在小波域的特征主要表現(xiàn)為異常大的系數(shù)，因此若某一源信號的32個尺度系數(shù)中，有某一系數(shù)的數(shù)值大于或者小于均值五倍方差，或者小于均值五倍方差，即將其視為噪聲。傳統(tǒng)方法一般將噪聲獨立成分中的全部系數(shù)置0，但是本文通過觀察發(fā)現(xiàn)，通過傳統(tǒng)預(yù)處理的信號，該成分中仍然包含著有用信息，因此本文選擇只將部分超過閾值的系數(shù)置0。其公式如下：

s（i， j）=0，cs（i， j）>閾值或cs（i， j）<閾值

cs（i， j），其他（11）

通過這種方式獲得的新的數(shù)據(jù)sC^s表示著消除結(jié)構(gòu)噪聲影響的尺度系數(shù)。

1.3改進的空域相關(guān)去噪

在文獻[13]中，Xu等[13]基于Witkin的研究基礎(chǔ)上提出了空域相關(guān)濾波算法，該算法中最重要的就是信號在經(jīng)過小波變換后，在各個尺度上具有極強的相關(guān)性，特別是處于信號邊緣的位置的相關(guān)性得到加強，而集中在小尺度上的噪聲系數(shù)的相關(guān)性得到抑制。而fMRI具有基本的圖像鄰域像素相關(guān)特性，表現(xiàn)在fMRI數(shù)據(jù)上即為空間相鄰體素間時間序列具有極強的相似性，本文通過比較鄰域信號在同一尺度上的小波系數(shù)，來抑制隨機噪聲。

首先定義空域相關(guān)：

Corrin（j，k）=Wfn（j，k）·Wfn（i，k）（12）

其中：Wfn（j，k）為處于空間j坐標位置k處的時間序列第n層的小波變換系數(shù)。為了使相關(guān)系數(shù)與小波系數(shù)具有可比性，將Corrin（j，k）歸一化到Wfn（j，k）上去，定義歸一化相關(guān)系數(shù)：

NewCorrin（j，k）=Corrin（j，k）pw（j）pCorrin（j）； k=1，2，…，N（13）

其中：

PCorrin（j）=∑Nk=1Corrin（j，k）2（14）

Pw（j）=∑Nk=1Wfn（j，k）2（15）

本文取空間3×3×3的空間立方體為鄰域空間，通過計算目標點與周邊26個鄰域體素各個尺度的空域相關(guān)值來鑒別信號的重要邊緣，即單個體素（除去各個腦區(qū)邊緣）的時間序列每一層小波變換系數(shù)均對應(yīng)有一組1×26的向量Ln，定義如下：

Ln（j，k）={NewCorr1n（j，k），NewCorr2n（j，k），…，NewCorr26n（j，k）}（16）

通過比較Ln均值與Wf（j，k）的大小，可以判斷Wf（j，k）是噪聲還是信號邊緣：若Wf（j，k）代表著信號的邊緣，儲存Wf（j，k）的位置和大小；若代表噪聲，則將該點置零。其公式如下：

NewWf（j，k）=Wf（j，k），mean（Ln）>Wf（j，k）

0，其他（17）

本文進行三層小波分解，因此通過在三個細節(jié)尺度上的計算，將得到的新系數(shù)與上文得到的尺度系數(shù)相結(jié)合得到完整的新小波系數(shù)，顯然該系數(shù)中不僅去除了結(jié)構(gòu)噪聲，而且去除了大部分隨機噪聲。最后重構(gòu)信號，可得到去噪后的數(shù)據(jù)。

2實驗及結(jié)果

2.1實驗數(shù)據(jù)

本實驗中具體掃描參數(shù)設(shè)置如下：重復(fù)時間（Time of Repetition， TR）1.5s，回波時間（Time of Echo， TE）27ms，視野（Field of View， FOV）24cm，獲得矩陣64×64，翻轉(zhuǎn)角70°，體素大小為3.75mm×3.75mm×4mm，層厚4mm，層間距1mm，29層，至下而上獲取圖像。聽力刺激任務(wù)由兩個8min的run組成，聲音是由電腦呈現(xiàn)系統(tǒng)通過內(nèi)嵌在核磁共振相容的衰減30dB的耳機發(fā)出。標準刺激是500Hz的音色，目標刺激是1000Hz的音色，新奇刺激是由非重復(fù)的隨機的電子噪聲構(gòu)成，比如掃地音、汽笛哨子音。目標刺激和新奇刺激出現(xiàn)的概率都為0.1，標準刺激出現(xiàn)的概率為0.8，刺激時長（duration）為200ms，刺激間歇隨機等概率出現(xiàn)1000ms，1500ms或者2000ms的間隔。所有刺激都為80dB，高于聽力的標準閾值。

2.2實驗結(jié)果

為了客觀評價本文去噪算法，本文分別通過比較并分析去噪前后數(shù)據(jù)的體素時間序列、時間簇分析曲線[14]、頭動指標[15]，以及腦區(qū)激活圖作為去噪質(zhì)量的評價標準。

2.2.1體素時間序列分析

圖2顯示了同一被試者的左側(cè)杏仁核腦區(qū)和右側(cè)海馬腦區(qū)去噪前后的時間序列曲線。可以明顯看出，經(jīng)過本文算法去噪后的時間序列相較原始時間序列，其峰值變化平緩，大部分的峰值已被置零。比較圖2（a）和圖2（b），可以看出本文算法針對不同腦區(qū)信號，有效抑制了偽激活，都具有良好的去噪效果。

2.2.2時間簇分析曲線

時間簇分析是一種對噪聲異常敏感的分析方法，圖3（a）為未經(jīng)去噪處理數(shù)據(jù)的時間簇分析曲線，圖3（b）為使用本文算法去噪后的時間簇曲線。從圖3（a）中可以看出，在掃描的前幾個時間序列中，由于被試并未適應(yīng)實驗環(huán)境，頭動較為嚴重。而由于頭動結(jié)構(gòu)噪聲的影響，原始信號的時間簇分析曲線能量主要集中在前幾個時間點，大腦響應(yīng)的刺激不明顯。而圖3（b）所示的腦區(qū)刺激，雖然在幾個時間點仍然受到結(jié)構(gòu)噪聲影響，但尖峰信號明顯減弱，其對大腦響應(yīng)刺激的觀測要明顯優(yōu)于處理前所示，這也說明了使用基于小波變換的分步去噪算法的優(yōu)越性。

2.2.3頭動指標分析

頭動是引起fMRI噪聲的一個重要因素，本文分別選用三個頭動指標去衡量去噪前后被試者的頭動情況，三個頭動指標分別為框架位移、DVARS（D為時間序列的導(dǎo)數(shù)， VARS為體素的均方根偏差）、尖峰百分比。其中一般要求被試的框架位移要小于2mm，框架位移越大，其信號的可靠性越低；而尖峰百分比衡量的是信號中尖峰信號出現(xiàn)的個數(shù)，其百分比越低，信號的質(zhì)量越好。表1為原始信號、文獻[16]方法、文獻[17]方法及基于本文算法的頭動指標。從表1中可以看出，與原始信號相比，三種去噪方法均能有效去除頭動影響，而綜合比較三個指標，本文算法具有較好的性能，特別是框架位移相較原始信號得到了大大降低。

2.2.4腦區(qū)激活點分析

圖4為腦區(qū)激活的平面及三維示意圖，其中圖4（a）為由未經(jīng)去噪的數(shù)據(jù)得到的激活腦區(qū)圖，從圖中可以明顯看出，激活的體素大部分集中在大腦的額葉部分，但是由于灰質(zhì)邊界的影響，這樣大范圍占據(jù)整個腦區(qū)的激活現(xiàn)象是不合理的，并且額葉主要負責思維，演算及于個體的需求和感情相關(guān)，與實驗設(shè)計的聽覺刺激并不吻合，因此該部分激活很有可能是由噪聲引起，因此該激活圖并不能很好地反映大腦神經(jīng)活動。與圖4（a）不同的是，經(jīng)過本文分步去噪算法后的圖4（b），腦區(qū)中包括顳上回、中回、下回前部區(qū)域的聽覺皮層均反映出一定程度的激活，這與實驗設(shè)計高度吻合，激活體素也正常分布在不同的腦區(qū)。

3結(jié)語

fMRI數(shù)據(jù)處理的重點就是捕獲“感興趣”信號，然而這類信號不僅能量只占全腦信號的一小部分，還常受到各類噪聲的影響；甚至在經(jīng)過傳統(tǒng)的一系類標準處理后，這些噪聲仍會影響數(shù)據(jù)分析的準確性。本文提出的基于小波變換的分步處理結(jié)構(gòu)噪聲及隨機噪聲的算法，在一定程度上減少了噪聲的影響，并且根據(jù)處理的數(shù)據(jù)分析得到的腦區(qū)激活圖，準確地反映了實驗設(shè)計對聽覺腦區(qū)的刺激。

本文算法針對傳統(tǒng)算法的改進在于，不僅在時間域上考慮了fMRI數(shù)據(jù)的信號特點，而且從fMRI數(shù)據(jù)具有圖像數(shù)據(jù)的特點出發(fā)，利用空間鄰域相關(guān)特性消除了部分隨機噪聲的影響。從實驗結(jié)果來看，同時考慮時間及空間特征可以很好地達到對信號去噪的效果，這也給今后處理fMRI數(shù)據(jù)提供了新的思路；

總言之，基于小波變換的分步去噪方法，可以減少噪聲對fMRI數(shù)據(jù)分析的影響。針對噪聲的特點，對經(jīng)過小波變換得到的尺度系數(shù)及細節(jié)系數(shù)進行分步處理。通過ICA算法分解尺度系數(shù)去除結(jié)構(gòu)噪聲源對數(shù)據(jù)的影響，而細節(jié)系數(shù)經(jīng)過改進的空域相關(guān)算法，其隨機噪聲的影響也得到了削弱。實驗結(jié)果表明，本文算法能提高數(shù)據(jù)反映神經(jīng)影響的能力，從而得到正確的腦區(qū)激活圖。但是本文算法也存在一些不足，譬如基于鄰域運算的空域相關(guān)算法，無法有效處理腦區(qū)邊界區(qū)域，今后將從全腦范圍內(nèi)著手解決該問題。

參考文獻：

[1]

薛紹偉，唐一源， 李健，等.一種基于fMRI數(shù)據(jù)的腦功能網(wǎng)絡(luò)構(gòu)建方法[J].計算機應(yīng)用研究，2010，27（11）：4055-4057.（XUE S W， TANG Y Y， LI J， et al. Method for constructing brain functional networks based on fMRI data [J]. Application Research of Computers， 2010， 27（11）： 4055-4057.）

[2]

SHIRER W R， JIANG H， PRICE C M， et al. Optimization of rsfMRI preprocessing for enhanced signalnoise separation， testretest reliability， and group discrimination [J]. Neuroimage， 2015， 117：67-79.

[3]

BULLMORE E， LONG C， SUCKLING J， et al. Colored noise and computational inference in fMRI time series analysis： resampling methods in time and wavelet domains [J]. Neuroimage， 2001， 13（6）： 86.

[4]

MALDJIAN J A， LAURIENTI P J， KRAFT R A， et al. An automated method for neuroanatomic and cytoarchitectonic atlasbased interrogation of fMRI data sets [J]. Neuroimage， 2003， 19（3）： 1233-1239.

[5]

FEIS R A， SMITH S M， FILIPPINI N， et al. ICAbased artifact removal diminishes scan site differences in multicenter restingstate fMRI [J]. Frontiers in Neuroscience， 2015， 9.

FEIS R A， SMITH S M， FILIPPINI N， et al. ICAbased artifact removal diminishes scan site differences in multicenter restingstate fMRI [EB/OL]. [20151115]. http：//europepmc.org/backend/ptpmcrender.fcgi？accid=PMC4621866&blobtype=pdf.

[6]

BRIGHT M G， MURPHY K. Is fMRI "noise" really noise？ Resting state nuisance regressors remove variance with network structure [J]. Neuroimage， 2015， 114： 158-169.

[7]

KARAHANOGLU F I， CABALLEROGAUDES C， LAZEYRAS F， et al. Total activation： fMRI deconvolution through spatiotemporal regularization [J]. Neuroimage， 2013， 73： 121-134.

[8]

KHALIDOV I， VAN DE VILLE D， FADILI J， et al. Activelets and sparsity： a new way to detect brain activation from fMRI data [J]. Proceedings of SPIE， 2007， 23（3）： 380-388.