基于小波变换的功能磁共振图像时间序列分步去噪

李微微　梅雪　周宇

摘要：

功能磁共振图像（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.