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

李微微　梅雪　周宇

摘要：

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

關鍵詞：

功能磁共振圖像；去噪；結構噪聲；隨機噪聲；小波變換

中圖分類號：

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.

英文關鍵詞Key words：

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

0引言

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

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

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

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

1本文算法

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

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

合適的小波基可以很好地將信號用小波系數稀疏表示，也決定著去噪的結果。基于血液動力學為一個穩定的線性系統的假設上，Khalidov 等[8]提出了Activelets小波基。本文使用Activelets小波基對fMRI時間序列進行三層分解與重構，得到新的尺度系數和細節系數。值得注意的是，Activelets小波基ψ對應的細節系數部分保證了血氧動力學中小波系數快速消退。具體而言，fMRI單體素時間序列可以寫成由尺度系數及細節系數構成的如下形式：

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

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

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

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

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

其中：hi及gi的濾波器組值的設定是根據血液動力學響應函數演變而來，詳見文獻[9]。鑒于它們的濾波器性質，其分解與重構仍基于經典的Mallat快速分解算法[10]。

1.2基于尺度系數的ICA分解與重構系數

為消除結構噪聲的能量，提高ICA分解信號的準確性，本文對單腦區分別使用ICA算法。假設單一腦區內有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）。因結構噪聲的能量主要集中在尺度系數上，其單個源信號可有小波變化后的尺度系數表出，公式如下：

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

為簡化公式，現表示為如下形式：

si=Csi·Φ（6）

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

S=Cs·Φ（7）

其中Cs為：

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

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

X≈CXΦ（9）

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

CX≈MCs（10）

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

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

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

通過這種方式獲得的新的數據sC^s表示著消除結構噪聲影響的尺度系數。

1.3改進的空域相關去噪

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

首先定義空域相關：

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

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

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個鄰域體素各個尺度的空域相關值來鑒別信號的重要邊緣，即單個體素（除去各個腦區邊緣）的時間序列每一層小波變換系數均對應有一組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）

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

2實驗及結果

2.1實驗數據

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

2.2實驗結果

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

2.2.1體素時間序列分析

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

2.2.2時間簇分析曲線

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

2.2.3頭動指標分析

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

2.2.4腦區激活點分析

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

3結語

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

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

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

參考文獻：

[1]

薛紹偉，唐一源， 李健，等.一種基于fMRI數據的腦功能網絡構建方法[J].計算機應用研究，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.