统计参数图和多尺度特征提取用于事件相关fMRI分析的比较

支联合 王华东 张 洁

1(周口师范学院物理与电子工程系，周口 466001)

2(周口师范学院计算机科学系，周口 466001)

3(北京中医药大学东直门医院放射科，北京 100700)

引言

基于血氧水平依赖性(blood oxygenation level dependence，BOLD)对比度机制成像的功能磁共振成像技术(functional magnetic resonance imaging，fMRI)，凭其非侵入性、可重复性和良好的时空分辨率等优点，成为探索大脑功能的理想工具［1］。一般地，BOLD-fMRI实验设计主要有组块型(blockdesigned paradigm)和事件相关型(event-related paradigm)两种［2］。其中，事件相关型设计方法可以有效地避免实验中被试者的期待、疲劳等因素的干扰，因此获得了广泛的应用。但是，由于刺激时间短，大脑对这种设计类型的响应信号(paradigm responsive signal，PRS)的幅度变化较小，而数据里又参杂有随机噪声、基线漂移等成分［3－5］。因此，分析事件相关fMRI数据的方法需要有高的灵敏度。

目前，较为权威的分析方法是 SPM(statistical parametric mapping)［5］。SPM 用实验刺激函数与血流动力学函数卷积的方法模拟PRS，用空间高斯平滑的方法去除每层数据里的高频噪声，用基于离散余弦变换的高通滤波方法，去除每个体素时间序列数据里的基线漂移。因此，SPM能较为全面地去除干扰成分，提高分析的灵敏度。

小波变换是现代信号分析技术之一，由于具有许多优良的特性而被广泛用于包括生物医学信号处理在内的信号分析领域［6－9］。文献［10］根据离散小波变换(discrete wavelet transform，DWT)具有多尺度分析的功能，能把一个复合信号里不同时频特征的成分分离到不同的小波尺度里，而fMRI数据里包含的PRS具有特定的时频属性，因此借助小波变换必能把它分离到特定的小波尺度里的事实，提出了基于DWT的特征提取方法来分析组块型fMRI数据。为了叙述方便，文中称这种方法为多尺度特征提取(multiscale feature extraction，MFE)。通过与SPM2的分析结果比较，发现新的方法显示出了较好的灵敏度［10］。然而，本研究发现，这种新的分析方法虽然是为组块型fMRI数据设计的，却更适合分析事件相关型fMRI数据。在以下的内容中，首先叙述视觉试验方法，然后是 MFE方法以及 MFE和SPM2对实验数据的分析结果，最后是讨论和结论。

1 材料和方法

1.1 视觉fMRI实验

3名健康被试者分别参加了视觉事件相关型和组块型fMRI试验，被试者均知情同意。对于视觉事件相关型设计，试验的任务刺激为棋盘格闪光刺激(持续2 s，频率8 Hz)，静息刺激为空白屏幕。任务刺激分别在第82、122、162和202 s时给予被试者，而其余时间则给予静息刺激。1.5 T Siemens扫描机获取了全脑120帧BOLD/EPI图像，每帧图像的大小为64体素×64体素×20体素，每个体素大小为3.4 mm×3.4 mm×6 mm(层间距1.2 mm)。扫描的时间分辨率 TR=2 s。对于视觉组块型设计，除了任务刺激的持续时间为10 s外，其余条件与事件相关试验相同。

在进一步分析前，这些视觉fMRI数据都经过了空间校准、标准化和空间高斯平滑处理，操作由SPM2软件完成。平滑时涉及到的半高宽为9 mm×9 mm×12 mm。另外，为了有效地去除脑外体素，每一个图像中的体素如果其值小于该图像中最大值的1/5，则该体素的值设置为0。

1.2 方法

1.2.1 离散小波变换

为理解上的方便，在叙述 MFE之前简要叙述DWT的基本原理。DWT分为离散小波分解和重构两个相反的过程，前者把时域信号变换成小波系数，后者再把小波系数逆变换为时域信号。对于能量有限的离散信号序列sj，离散小波分解可表示为

式中:H和G分别为低通和高通滤波器，由小波函数和与之对应的尺度函数得到;D0为二进下采样算子，用于采集滤波器输出的偶数次项;sj+1为sj中的低频成分，dj+1为高频成分;j=0时，sj表示原始信号序列。

显然，如果连续分解 J次，则s0中J+1个不同时频特征的序列成分 d1，d2，…，dJ和 sJ会被一一抽取出来，依次分离到J+1个不同的小波尺度里，这就是小波变换的多尺度分解。

这样，对fMRI一个体素的时间序列数据进行小波分解后，如果知道了 PRS存在的小波尺度，比如假如存在于dj里，则只需保留尺度dj里的小波系数不变，而将其他尺度里的系数设置为0，然后进行小波重构，就能获取隐藏在这一序列数据里的PRS。

1.2.2 多尺度特征提取

文献［10］视一个体素的时间序列数据 X里包含的PRS为X的特征，视小波分解这一序列数据后PRS存在的小波尺度为特征小波尺度(feature wavelet scales，FWS)，然后基于 DWT来提取PRS。

步骤 1:确定 FWS。先由 SPM2软件模拟出PRS，再比较模拟的 PRS的频谱以及离散小波多尺度分解PRS后获得的各个小波尺度重构信号的频谱，由此确定 PRS经过小波分解后存在的小波尺度，即 FWS。

步骤2:特征提取。先对X进行离散小波多尺度分解，再保留FWS尺度里的系数不变，而设置其他尺度里的系数为0，后进行小波重构，得到的信号记为eX。

步骤3:激活检查。为简单计，激活检测通过相关分析eX和参考信号 ePRS实现，ePRS为对模拟的PRS施行同步骤2一样的小波变换操作得到的，相关系数r通过Fisher’s Z变换转化为Z-值，其中ePRS为数据序列的长度。

为对比MFE分析性能的好坏，采用SPM2在相同条件下分析视觉数据。其中，SPM2用到的高通滤波的截止周期为该软件缺省的128 s，分析每个被试最后得到的 t-值，通过 SPM2自带的程序转化为Z-值。两种方法得到的Z-值阈值化，显著性水平选为P=0.001，且经过高斯随机场理论的修正，修正方法通过SPM2完成。

图1显示了基于组块型fMRI实验设计模拟的PRS经过频谱分析后得到的频谱(左图)，图中阴影部分为SPM2用到的截止周期128 s的高通滤波的范围，以及 MFE涉及到的 ePRS的频谱(右图)。MFE用到的小波为 sym3，FWS为 d3和 d4。显然，模拟的 PRS包含0.075、0.05、0.025 Hz和直流共4个频率成分。而SPM2通过高通滤波去掉了直流成分，MFE则通过选取尺度d3和d4为特征尺度的方法主要保留了0.075、0.05、0.025 Hz的成分。

图1 SPM2和MFE分析组块型fMRI数据涉及到的PRS的频谱(纵轴为相对谱密度)Fig.1 The related frequency components of PRS used by SPM2 and MFE for the block designed fMRI data(The vertical axis is the relative spectral density)

图2显示了基于事件相关型fMRI实验设计模拟的PRS经过频谱分析后得到的频谱(左图)，图中阴影部分为SPM2用到的截止周期为128 s的高通滤波的范围，以及 MFE涉及到的 ePRS的频谱(右图)。MFE用到的小波为 sym2，FWS为 d2、d3和d4。显然，模拟的PRS包含7种频率成分，SPM2通过高通滤波去掉了其中的直流成分，而MFE则通过选取d2、d3和d4为特征尺度的方法也去掉了直流成分。

图2 SPM2和MFE分析事件相关型fMRI数据涉及到的PRS的频谱(纵轴为相对谱密度)Fig.2 The related frequency components of PRS used by SPM2 and MFE for the event-related fMRI data(The vertical axis is the relative spectral density)

2 结果

图3显示了SPM2和MFE分析视觉组块型fMRI试验3个被试一层数据(Talairach坐标 z=0 mm)的激活Z-图。从图上可知，SPM2检测到的激活脑区都是公认的视觉区，为双侧枕叶，主要是舌回和距状回，但MFE检测到的这些视觉区显然要比SPM2的大。进一步分析检测结果，SPM2检测3个被试视觉区体素个数分别为426、534和238，MFE分别为 871、1 230和 600，因此前者只是后者的48.9%、43.4%和39.7%。就检测到的视觉激活区的平均Z值而言，SPM2分别为7.36、7.19和6.43，而MFE分别为10.84、10.3和8.41。但应该注意的是，除视觉区外，MFE在非视觉区也检测到激活。

图3 SPM2和MFE分析视觉组块型fMRI试验3个被试一层数据(Talairach z=0mm)的激活 Z-图(显著性水平为P=0.001，经过高斯随机场理论修正)Fig.3 The active Z-maps(Talairach coordinate z=0mm)generated by SPM2 and MFE to the visual blocked paradigm fMRI data for three subjects.The significantlevelisP =0.001，corrected based on Gaussian random field.

图4显示了SPM2和MFE分析视觉事件相关型fMRI试验3个被试一层数据(Talairach坐标轴向z=0 mm)的激活 Z-图。显然，与分析组块型 fMRI数据不同，SPM2和MFE两种方法分析事件相关型数据检测到的激活区都在公认的视觉区里。也就是说，如果以SPM2的激活区为标准的话，MFE在处理这种类型的数据时表现出了良好的特异性。但非常明显的是，检测到的激活体素个数SPM2远少于MFE。进一步分析3个被试的检测结果可知，SPM2检测到的激活体素个数分别是94、18和55，MFE分别是 191、292和 375，前者只是后者的49.2%、6.2%和14.7%。而就检测到的激活区的平均 Z值而言，SPM2分别为 6.75、5.69和 5.91，MFE分别为8.08、6.37和6.61。因此，综合两种类型的数据处理结果可知，同SPM2相比，MFE在分析事件相关型fMRI数据时同时显示出了极高的灵敏度和良好的特异性。

图4 SPM2和MFE分析视觉事件相关型fMRI试验3个被试一层数据(Talairach z=0mm)的激活Z-图(显著性水平为P=0.001，经过高斯随机场理论修正)Fig.4 The active Z-maps(Talairach coordinate z=0mm)generated by SPM2 and MFE to the visual eventrelated paradigm fMRI data for three subjects(The significantlevelisP =0.001，corrected based on Gaussian random field)

3 讨论

灵敏性和特异性是考察新方法的重要指标。就本研究而言，灵敏性是指新方法正确地检测激活体素的能力，特异性则是指新方法区分激活体素和非激活体素的能力，也就是正确地检测非激活体素的能力。因为SPM2是目前国际上权威的处理fMRI数据的方法，因此把它的处理结果作为“金标准”是可以接受的。显然，无论是处理组块型还是事件相关型fMRI数据，MFE都比 SPM2灵敏，这归结于两种方法去除干扰成分的能力不同。SPM2虽然去除了每个体素时间序列数据里的低频成分，但仅此而已。与SPM2不同，MFE不但去除了体素时间序列数据里的低频成分，而且还去除了高频噪声成分，因为它在处理时舍弃了高频成分存在的小波尺度。而MFE在处理视觉事件相关型数据时，表现出了比处理组块型数据更高的灵敏性和更好的特异性，这可归结为:一方面，视觉事件相关型数据里 PRS信号的变化幅度弱小，信噪比较低，因此受噪声的影响比较严重，对去噪更灵敏;另一方面，MFE在处理视觉事件相关型数据和选择参考信号时保留了更多的小波尺度，因为它只舍弃了d1尺度里的小波系数，而在处理组块型数据时，除d1尺度外还舍弃了d2尺度里的小波系数。这也是我们强烈推荐MFE处理事件相关型数据的理由和根据。

Donoho提出的小波变换去噪技术使小波变换得到了广泛的应用［11］。利用这种技术去噪，遵循下面这些步骤:对已知信号进行小波分解;依据一定的规则选择适当的阈值，并依据一定的规则对小波系数进行修正;对修正后的小波系数进行重构，得到去噪信号。显然，MFE方法遵循上述步骤，但采取了一个相对简单的方法来对原始的小波系数进行修正。它是根据模拟出的PRS在小波尺度里的能量分布来确定所谓的特征小波尺度，进而确定哪些尺度里的小波系数原封不动地保留，哪些尺度里的系数设置为0。本研究采用的这种方法虽然简单，但却是非常合理的。因为，特征尺度是 PRS存在的尺度，因此里面的系数不能受到损害;而在非特征尺度里，或者完全没有 PRS存在，因而完全是干扰成分，或者PRS的能量存在，但干扰成分很多，因此小波系数必须舍弃。MFE的合理性基于一个非常重要的前提，即 PRS必须已知，进而借助频谱分析也能够确定特征小波尺度。Donoho小波去噪技术的假设是:绝对值小的小波系数与噪声有关，绝对值大的与有用信号有关，至于有用信号确切地存在于哪些小波尺度则不知道，因此相应的阈值选取规则和小波系数修正规则就复杂得多。

任何一种数据分析方法都有自己的使用条件，MFE也不例外。如上所述，MFE需要事先知道大脑响应试验刺激的函数来确定特征小波尺度，所以MFE是一种模型驱动的方法，而非Donoho小波去噪技术那样是一种数据驱动的方法。因此，如果fMRI试验无法事先确定大脑响应试验刺激的函数，则这样的数据无法用MFE分析。另外，因为MFE借助小波变换，提取大脑体素fMRI时间序列数据里的PRS，因此与小波变换有关的限制条件也同样适用于MFE。比如，小波变换需要事先选择最合适的小波，以得到最优的分析结果，MFE也同样如此。经过广泛的尝试才发现，sym3是视觉组块型 fMRI数据的最优小波，而sym2是视觉事件相关型fMRI数据的最优小波，因此寻找普适的小波是本课题组下一步的任务。

4 结论

总之，MFE是一个非常值得推荐的分析事件相关型fMRI数据的方法，这源于它简单的设计、优异的灵敏性和良好的特异性。

［1］Ogawa S，Lee TM，Kay AR，et al.Brain magnetic resonance imaging with contrast dependent on blood oxygenation［J］.Proc Natl Acad Sci，USA，1990，87(24):9868－9872.

［2］包尚联.脑功能成像物理学［M］.郑州:郑州大学出版社，2006:200－207.

［3］Hossien-Zadeh GA， Soltanian-Zadeh H， Ardekani A.Multiresolution fMRI activation detection using translation invariant wavelet transform and statistical analysis based on resampling［J］.IEEE Trans Med Imaging，2003，22(3):302－314.

［4］Tanabe J， MillerD， TregellasJ， etal. Comparison of detrendingmethod foroptimalfMRIpreprocessing ［J］.NeuroImage，2002，15(4):902－907.

［5］Friston KJ， HolmesAP， Worsley KJ， etal. Statistical parametric maps in functional imaging:a general linear approach［J］.Hum Brain Mapping，1995，2(4):189 －210.

［6］Bullmore E，Fadili J，Maxim V，et al.Wavelets and functional magnetic resonance imaging ofthe human brain ［J］.NeuroImage，2004，23(1):S234－S249.

［7］侯建华，熊承义，何翔，等.基于小波统计模型的医学超声图像去噪方法研究［J］.中国生物医学工程学报，2009，28(1):31－36.

［8］孙炜，王耀南.基于模糊小波神经网络的磁共振图像分割方法［J］.中国生物医学工程学报，2006，25(3):267－270.

［9］杨建国.小波分析及其工程应用［M］.北京:机械工业出版社，2005:107－128.

［10］支联合，李玉晓，赵书俊，等.基于离散小波变换的 fMRI数据特征提取［J］.中国医学影像技术，2010，26(6):84－87.

［11］Donoho DL，Johnstone IM.Adapting to unknown smoothness via wavelet shrinkage ［J］.J Am Stat Assoc，1995，90(432):1200－1224.