CN1775172A - 基于独立成份空间相关性消除功能磁共振数据噪声的方法 - Google Patents

Abstract

本发明公开了一种基于独立成份空间相关性消除功能磁共振数据噪声的方法，先对图像的灰质区域和脑脊液区域进行图像分割，分别对其进行主成份分析、频谱分析，确定并消除主成分分量中的随机噪声成分，重建消除随机噪声后的灰质数据和脑脊液数据，对重建后数据进行独立成份分解，利用独立成份分别构造矩阵，再进行典型相关分后，排序，将这些成分置零后得到一组新的独立成份，重建消除各种噪声成份后的数据，重复以上步骤直至消除各层灰质数据与脑脊液数据中相关最大的各种噪声成分。本发明能够有效地消除多层功能磁共振数据中生理噪声，能够消除与脑脊液区域数据中相关最大的其它低频噪声成份，较好地实现了多层功能磁共振数据中噪声的消除。

Description

基于独立成份空间相关性消除功能磁共振数据噪声的方法

技术领域

本发明涉及功能磁共振成像技术中数据处理领域，尤其涉及一种基于独立成份空间相关性消除功能磁共振数据噪声的方法。

背景技术

基于血氧水平依赖(Blood Oxygenation Level Dependent，BOLD)对比度机制，功能磁共振成像(functional MRI，fMRI)能够实时地检测与大脑皮层神经活动相关联的局部血流动力学改变。但是，基于BOLD效应的磁共振信号改变非常小，在磁场强度为1.5T时只有0.5～2％。此外，在序列图像采集过程中周期呼吸和心脏运动导致的脑组织脉动以及其它一些复杂的生理运动都会导致许多生理相关的MRI信号改变。这些生理相关的MRI信号改变构成了功能磁共振数据的生理噪声。除了与心脏和呼吸相关的生理噪声之外，功能磁共振数据中还包含源于射频线圈和受检者的随机噪声以及其它低频噪声，这些低频噪声主要是由于扫描器的不稳定、运动校正后残留的运动关联效应、长期的生理漂移等因素造成的。功能磁共振数据中这些噪声成份与神经活动无关，不仅降低了对功能磁共振数据统计分析的显著性，而且直接影响功能磁共振的敏感性及空间定位的可靠性。

功能磁共振数据中的噪声主要是由于与心脏和呼吸相关的非神经性的生理效应引起的。除了与心脏和呼吸相关的生理噪声之外，功能磁共振数据中还包含源于射频线圈和受检者的随机噪声以及其它低频噪声，这些低频噪声主要是由于扫描器的不稳定、运动校正后残留的运动关联效应、长期的生理漂移等因素造成的。因此，功能磁共振数据模型可以表示为：

y_i(t)＝s_i(t)+p_i(t)+r_i(t)                                  (1)

其中，y_i(t)表示每个象素的时间过程，s_i(t)表示真正与神经活动相关的功能信号改变，p_i(t)表示非神经性的生理效应引起的生理的噪声及各种低频噪声，r_i(t)表示系统的随机噪声。

为了减少由于生理运动导致的功能磁共振信号改变，已经提出了许多有效的方法，如采用心电门控扫描、导航回波、数字滤波、回溯性校正及基于k-空间的估计和校正技术等。这些方法或者要求实时采集或监测生理运动信息，或者要求特定的扫描序列，都是以增加数据采集负担或数据后处理复杂性为代价的。此外，为了捕捉生理运动周期，这些方法要求功能磁共振数据的采集频率高于生理运动频率。一些学者最近提出采用成分分析方法对功能磁共振数据进行成分分解，然后通过对分离成分进行频谱估计确定并消除相应生理噪声成分。虽然这种方法不需要额外采集生理运动信息，数据后处理较为简单。但是，为了估计并确定生理噪声成分，该方法仍然要求功能数据的采集频率高于生理运动频率。

为了覆盖足够的脑空间，功能磁共振常采用多层扫描方式采集数据。由于时间序列图象的采样率低于生理运动频率，导致高频生理噪声混叠到低频成分，生理噪声不再保持原来的波形，无法基于频谱分析方法确定多层功能磁共振数据的生理噪声成分。

发明内容

本发明提供了一种基于独立成份空间相关性消除功能磁共振数据噪声的方法，可以有效地消除多层功能磁共振数据中生理噪声，解决现有的消除相应生理噪声成分必须以增加数据后处理复杂性为代价或者必须要求功能数据的采集频率高于生理运动频率等苛刻条件的问题。

本发明采用的技术方案是：

一种基于独立成份空间相关性消除功能磁共振数据噪声的方法，采用下列步骤：

1)设功能磁共振体数据序列的时间采样点数为N，每个体数据包含S层图像，分别以第i层图像的灰质区域和脑脊液区域为模板，构造该层时间序列图像对应的灰质数据矩阵Y^g和脑脊液数据矩阵Y^c，矩阵Y^g和Y^c的行表达了像素的时间过程；

2)对灰质数据矩阵Y^g进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前p个特征值对应的特征向量的累积方差贡献率大于90％选取灰质数据矩阵Y^g的主成份分量；

3)对脑脊液数据矩阵Y^c进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前q个特征值对应的特征向量的累积方差贡献率大于90％选取脑脊液数据矩阵Y^c的主成份分量；

4)采用傅立叶频谱分析方法分别对灰质数据矩阵Y^g和脑脊液数据矩阵Y^c的主成分分量进行频谱分析，根据随机噪声判定准则确定并消除灰质数据矩阵Y^g和脑脊液数据矩阵Y^c主成分分量中的随机噪声成分；

5)分别利用步骤4)保留的主成分分量重建消除随机噪声后的灰质数据

和脑脊液数据

6)对灰质数据 进行独立成份分解，得到分离矩阵W_g及m个统计上相互独立的成份c_i ^g(t)，i＝1，2，Λm；

7)对脑脊液数据

进行独立成份分解，得到分离矩阵W_c及n个统计上相互独立的成份c_j ^c(t)，j＝1，2，Λn；

8)利用灰质数据的m个独立成份c_i ^g(t)构造矩阵X_m×N，利用脑脊液数据的n个独立成份c_j ^c(t)构造矩阵Y_n×N；

9)对矩阵X_m×N和Y_n×N进行典型相关分析，得到c_i ^g(t)的加权系数矩阵A_m×r、c_j ^c(t)加权系数矩阵B_n×r及包含r个元素的行向量 r为m和n中较小的一个数，

的第k个元素为矩阵X_m×N和Y_n×N的各行分别以矩阵A_m×r和B_n×r的第k列元素为加权系数加权后得到的典型变量的典型相关系数的平方，k＝1，2，Λr；

10)计算行向量 的各元素对矩阵A_m×r各行元素的加权和，得到表达灰质数据中各独立成份c_i ^g(t)与脑脊液数据中各独立成份c_j ^c(t)相关测度的列向量

11)对列向量

的各元素从大到小排序，那么，排序后的列向量 的前r个元素对应的灰质数据的独立成分c_i ^g(t)为与脑脊液数据中各独立成份c_j ^c(t)相关最大的噪声成分，将这些成分置零后得到一组新的独立成份

12)用分离矩阵W_g的逆矩阵乘以由独立成份 构成的成份矩阵，重建消除与脑脊液数据中各独立成份c_j ^c(t)相关最大的各种噪声成份后的灰质数据；

14)重复步骤1)到12)，直至消除功能磁共振时间序列数据各层灰质数据中的与脑脊液数据中相关最大的各种噪声成分。

与现有的技术相比，本发明具有如下优点：

本发明提出基于独立成份空间相关性消除功能磁共振数据噪声的方法基于生理噪声成分时间过程在大脑灰质区域以及脑脊液区域的空间相关性，通过对脑脊液区域和灰质区域数据进行独立成份分解，利用典型相关分析方法分析脑脊液区域和灰质区域独立成份时间过程的空间相关性，能够识别出灰质区域数据中与脑脊液区域数据中相关最大的噪声成份，本发明提出的方法不仅能够有效地消除多层功能磁共振数据中生理噪声，而且能够消除与脑脊液区域数据中相关最大的其它低频噪声成份，较好地实现了多层功能磁共振数据中噪声的消除。

由于本发明提出的方法数据处理简单，不需要实时采集或监测生理运动信息，可以对早期采集的功能磁共振数据进行离线处理，因此是一种较为适用的功能磁共振数据处理方法，非常适于临床医院及研究单位利用功能磁共振进行脑科学研究。

附图说明

图1为听觉功能实验刺激时序图

图2-1为相关系数阈值T＝0.60时对实验数据的第30层图像直接进行相关分析检测到的可能的激活像素

图2-2为相关系数阈值T＝0.25时对第30层图像直接进行相关分析检测到的可能的激活像素

图2-3为采用本发明提出的方法消除灰质数据中生理相关的各种噪声后，相关系数阈值T＝0.60时检测到的激活像素

图2-4为将图2-3的检测结果映射到实验数据第30层图像上的结果

图3-1为灰质数据中的10个独立成份

图3-2为脑脊液数据中的5个独立成份

图4-1为某一激活像素消除噪声之前的时间过程

图4-2为某一激活像素消除噪声之后的时间过程

图5-1为消除噪声前后激活像素时间过程与功能刺激参考函数的相关系数比较

图5-2为消除噪声前后激活像素时间过程的方差比较

图6-1到图6-6为采用本发明提出的噪声消除方法对实验数据进行处理之后，采用相关分析法得到的从第28到第33层连续6层图像的听觉功能激活映射图。

具体实施方式

一种基于独立成份空间相关性消除功能磁共振数据噪声的方法，采用下列步骤：

1)设功能磁共振体数据序列的时间采样点数为N，每个体数据包含S层图像，分别以第i层图像的灰质区域和脑脊液区域为模板，构造该层时间序列图像对应的灰质数据矩阵Y^g和脑脊液数据矩阵Y^c，矩阵Y^g和Y^c的行表达了像素的时间过程；

2)对灰质数据矩阵Y^g进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前p个特征值对应的特征向量的累积方差贡献率大于90％选取灰质数据矩阵Y^g的主成份分量；

3)对脑脊液数据矩阵Y^c进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前q个特征值对应的特征向量的累积方差贡献率大于90％选取脑脊液数据矩阵Y^c的主成份分量；

4)采用傅立叶频谱分析方法分别对灰质数据矩阵Y^g和脑脊液数据矩阵Y^c的主成分分量进行频谱分析，根据随机噪声判定准则确定并消除灰质数据矩阵Y^g和脑脊液数据矩阵Y^c主成分分量中的随机噪声成分；

5)分别利用步骤4)保留的主成分分量重建消除随机噪声后的灰质数据

和脑脊液数据

6)对灰质数据 进行独立成份分解，得到分离矩阵W_g及m个统计上相互独立的成份c_i ^g(t)，i＝1，2，Λm；

7)对脑脊液数据

进行独立成份分解，得到分离矩阵W_c及n个统计上相互独立的成份c_j ^c(t)，j＝1，2，Λn；

8)利用灰质数据的m个独立成份c_i ^g(t)构造矩阵X_m×N，利用脑脊液数据的n个独立成份c_j ^c(t)构造矩阵Y_n×N；

9)对矩阵X_m×N和Y_n×N进行典型相关分析，得到c_i ^g(t)的加权系数矩阵A_m×r、c_j ^c(t)加权系数矩阵B_n×r及包含r个元素的行向量

r为m和n中较小的一个数，的第k个元素为矩阵X_m×N和Y_n×N的各行分别以矩阵A_m×r和B_n×r的第k列元素为加权系数加权后得到的典型变量的典型相关系数的平方，k＝1，2，Λr；

10)计算行向量 的各元素对矩阵A_m×r各行元素的加权和，得到表达灰质数据中各独立成份c_i ^g(t)与脑脊液数据中各独立成份c_j ^c(t)相关测度的列向量

11)对列向量

的各元素从大到小排序，那么，排序后的列向量

的前r个元素对应的灰质数据的独立成分c_i ^g(t)为与脑脊液数据中各独立成份c_j ^c(t)相关最大的噪声成分，将这些成分置零后得到一组新的独立成份

12)用分离矩阵W_g的逆矩阵乘以由独立成份

构成的成份矩阵，重建消除与脑脊液数据中各独立成份c^c _j(t)相关最大的各种噪声成份后的灰质数据；

15)重复步骤1)到12)，直至消除功能磁共振时间序列数据各层灰质数据中的与脑脊液数据中相关最大的各种噪声成分。

本实施例还可以通过计算功能实验设计时序与上述步骤1)所述灰质数据矩阵Y^g中每行所表达像素的时间过程的相关系数，用相关系数大于阈值T的像素重新构造灰质数据的矩阵。上述随机噪声判定准则是：如果某一成分傅立叶频谱的平均功率大于或等于频谱功率的标准偏差，则认为该成分为随机噪声成分。

我们通过对实际的听觉功能磁共振数据进行处理与分析来验证提出的方法。实验数据由伦敦大学Wellcom神经成像实验室提供。听觉实验包含静息和激发两种状态，静息期不施加任何刺激，激发期以每分钟60个双音节英文单词的速度施加听觉刺激。这样的实验共重复8次，每次实验中，对静息期和激发期分别进行6次采样，每次采样得到64幅连续的横轴位断层图象，扫描矩阵为64×64×64，单次扫描采样重复时间TR＝7s，整个实验持续6分钟。具体实施方式如下：

1)利用目前较为通用的功能磁共振统计分析软件SPM对实验数据进行运动校正、空间标准化等预处理后，得到包含96个79×95×68体数据的数据集，并对数据集中第一个体数据进行图象分割，得到灰质(GM)、白质(WM)及脑脊液(CSF)，并以分割的灰质和脑脊液图像作为模板，分别构造灰质数据集和脑脊液数据集。

在实验中，我们首先以听觉功能实验的刺激时序作为参考函数ref(t)，如图1所示，利用(7)计算ref(t)与灰质数据集中每一像素时间过程y_i ^g(t)的相关系数。

$d=\frac{\mathrm{\Σ}(y_i^g\left(t\right)-{\stackrel{\‾}{a}}_i)(\mathrm{ref}\left(t\right)-\stackrel{\‾}{b})}{\sqrt{\mathrm{\Σ}{(y_i^g\left(t\right)-{\stackrel{\‾}{a}}_i)}^2\mathrm{\Σ}{(\mathrm{ref}\left(t\right)-\stackrel{\‾}{b})}^2}}---\left(7\right)$

其中 a_i、 b分别表示灰质数据中的第i个像素时间过程y_i ^g(t)的平均值以及参考函数ref(t)的时间过程的平均值。

如果相关系数大于设定的阈值T，则认为该象素为可能的激活象素。然后由这些可能的激活像素重新构造灰质数据集。经过这样的初筛，处理的数据大量减少，一方面数据处理效率相应提高，另一方面不仅保证了灰质数据集的主成分中包含真正与神经活动相关的功能信号成份，而且对灰质数据集独立成分估计的准确度也提高。

相关系数阈值T可以根据统计检验的显著性水平确定。这里，我们选取T＝0.60，相应的Z分位数大于6，统计检验的显著性水平α小于0.01。如图2-1中所示高亮区域显示了阈值T＝0.60时实验数据第30层图像中27个可能的激活象素。由于噪声的影响，直接利用相关分析方法只检测到少数可能的激活象素。为了增加可能的激活象素的样本，我们降低阈值T。图2-2中高亮区域显示了阈值T＝0.25时实验数据第30层图像中262个可能的激活象素。比较图2-1和图2-2，由于相关系数阈值降低，检测到的可能的激活象素明显增多，但是伪激活象素也明显增加。这种结果正是由于功能磁共振数据中包含各种噪声成份导致的。功能磁共振空间定位的可靠性降低。所以，为了消除功能磁共振数据中包含的各种噪声成份，我们用图2-2中262个可能的激活象素重新构造灰质数据集。

2)对重新构造的灰质数据集和脑脊液数据集分别采用奇异分解算法(SVD)进行主成分分析，并根据保留的特征值对应的特征向量的累积方差贡献率大于90％选取主成份分量，对灰质数据和脑脊液数据的主成份利用傅立叶变换进行频谱分析，如果某一成分整个频谱的平均功率大于或等于整个频谱功率的标准偏差，则认为该成份为系统随机噪声成份，并将该成份置零，消除随机噪声成份。经过这样的处理，灰质数据集中保留10个主成分，脑脊液数据集中保留5个主成分。

3)分别利用步骤2)得到的主成分分量重建灰质数据

和脑脊液数据

通过计算功能实验设计时序与每一图象灰质区域象素的时间过程的相关系数，用相关系数大于阈值T的像素重新构造灰质数据的矩阵。

4)对消除随机噪声后的灰质数据采用Hyvarinen的快速定点独立成分分析算法(fastICA)进行独立成分分解，结果如图3-1所示。

5)对消除随机噪声后的脑脊液数据分别采用与4)同样的算法进行独立成分分解，结果如图3-2所示。

6)为了确定灰质数据中的生理相关的噪声成分，我们分别利用图3-1、图3-2所示的独立成份构造一个10×96的矩阵X和一个5×96矩阵Y。

7)对矩阵X和Y进行典型相关分析，可得到灰质数据中独立成份c_i ^g(t)的加权系数矩阵A_5×5、脑脊液数据中独立成份c_j ^c(t)加权系数矩阵B_10×5及包含5个元素的行向量

行向量 的第k个元素为矩阵X和Y的各行分别以矩阵A_m×r和B_n×r的第k列元素加权后的典型相关系数的平方。

8)计算行向量 的各元素对矩阵A_m×r各行元素的加权和，得到表达灰质数据中各独立成份c_i ^g(t)与脑脊液数据中各独立成份c_j ^c(t)的相关测度的列向量

对列向量 的各元素从大到小排序，结果如表1-1所示：

                       表1灰质和脑脊液中独立成份时间过程典型相关分析结果(第30层)

M	1.1833	0.9625	0.8489	0.8255	0.8167	0.7313	0.6081	0.5389	0.5238	0.2226
											IC	1	7	10	5	6	4	8	9	3	2

从表1的排序结果可见，图3-1所示的灰质数据集中10个独立成份与图3-2所示的脑脊液数据集中的5个独立成分相关的相关测度M从大到小一依次为第1、7、10、5、6、4、8、9、3、2。图3-1所示的第2个独立成分与图3-2所示的脑脊液数据集中的5个独立成分相关的相关测度最小，该独立成分与图2所示的听觉功能实验刺激时序基本一致，反映了大脑灰质区域对听觉功能刺激响应的时间过程。

9)根据前面介绍的典型相关分析方法，我们选取表1前5个独立成分作为与脑脊液数据中5个独立成份相关最大的生理相关的噪声成分。将这些成分置零后，

10)根据(8)重建灰质数据集 $\widetilde{Y^g}=W^{-1}{C}_0^g,$ 其中，W^-1表示分离矩阵W的逆矩阵，C₀ ^g表示将噪声成份置零后的灰质成份矩阵。图4-1、图4-2中实线分别表示同一激活像素在消除噪声之前和消除噪声之后的时间过程，虚线表示反映听觉功能实验刺激时序的参考函数的时间过程。可见，消除噪声之后，激活像素时间过程的波形明显改善，听觉功能实验的静息与激发两种状态对比也更加明显。

为了进一步定量说明噪声的消除效果，我们分别计算像素的时间过程在噪声消除之前和之后与参考函数的相关系数以及像素的时间过程在噪声消除之前和之后的方差。图6-1中虚线显示了重构的灰质数据集中262个可能的激活像素的时间过程在消除噪声之前与参考函数的相关系数，实线显示了消除噪声之后与参考波形的相关系数，可见，消除噪声之后相关系数明显变大，之前相关系数最大值为0.6715，平均值为0.4261；之后相关系数最大值为0.7547，平均值为0.6097，相关系数平均值增加43.09％。6-2中虚线和实线分别显示了重构的灰质数据集中262个像素的时间过程在消除噪声之前和之后的均方差。可见，消除噪声之后均方差明显变小，之前均方差为5541.2，之后均方差为3110.6，均方差减小43.86％。

图2-3的高亮区域显示了消除噪声之后与参考函数的相关系数大于0.60的激活像素。与图2-1相比，虽然选取同样的相关系数阈值，但是在对功能磁共振数据统计分析的显著性保持不变的情况下，检测到的激活像素个数增加到160个，功能磁共振的敏感性提高。与图2-2相比，伪激活像素个数也明显减少，功能磁共振空间定位的可靠性提高。图2-4为将图2-3映射到实验数据集第30层图像上的听觉功能激活映射图。

表2为功能磁共振数据消除噪声前后与参考函数相关系数及方差改变比较。Num_B、Num_A分别表示消除噪声之前和之后各数据层可能的激活像素。消除噪声之前的相关系数阈值为T＝0.25，消除噪声之后的相关系数阈值为T＝0.60。dmax_B、dmax_A分别表示消除噪声前、后激活像素的最大相关系数，dmean_B、dmean_A分别表示消除噪声前、后激活像素的相关系数平均值。Δd％表示消除噪声之后相关系数平均值改变百分比。Cmean_B、Cmean_A分别表示消除噪声前、后激活像素的均方差，Δv％表示消除噪声前后激活像素的均方差百分比改变。

从表2可见，采用本文提出的方法消除噪声功能磁共振数据的噪声后，激活像素的时间过程与参考函数的最大相关系数、平均相关系数都一致增大，激活像素时间过程的方差都一致减小。检测到的激活像素数量明显增多，伪激活像素数量也明显减少。

                表2功能磁共振数据消除噪声前后与参考函数相关系数及方差改变比较

	Num_B	dmax_B	dmax_A	dmean_B	dmean_A	Δd％	Cmax_B	Cmax_A	Δv％	Num_A
											Slice28	169	0.7046	0.7309	0.4090	0.5228	27.82％	4854.8	3179.5	52.69％	55
Slice29	211	0.6907	0.7380	0.4342	0.5416	24.74％	6248.6	4206.3	48.55％	88
											Slice30	262	0.6715	0.7547	0.4261	0.6097	43.09％	5541.2	3110.6	43.8％	160
Slice31	295	0.6867	0.7208	0.4168	0.5347	28.27％	4615.3	2184.6	111.3％	129
											Slice32	274	0.6475	0.7002	0.4007	0.5582	39.31％	3903.9	2420.7	73.40％	137
Slice33	183	0.5705	0.6348	0.3628	0.4557	25.61％	2846.1	1583.1	79.78％	9

图6-1到图6-6分别显示了先采用提出的噪声消除方法对实验数据进行处理之后，然后采用相关分析法得到的从第28到第33层连续6层图像的听觉功能激活映射图。从检测结果可见：大脑左右两侧的上颞叶及大脑边缘叶显示较显著的激活，激活区域检测结果与经验的听觉中枢功能解剖区域基本相符；大脑左右两侧的语言沟回区也显示少量的激活信号，这是由于该功能实验采用双音节英文单词作为听觉刺激时，受检者的语言理解区功能也有响应的缘故，该结果进一步说明了本文提出消除功能磁共振数据的噪声方法的有效性及可靠性。

本发明的原理：

随机噪声判定准则是：如果某一成分傅立叶频谱的平均功率大于或等于频谱功率的标准偏差，则认为该成分为随机噪声成分。

主成分分析(Principal Component Analysis，PCA)是一种基于二阶统计量的信号分解技术，其目的是将信号分解为若干统计意义上相互正交(或不相关)的成分。独立成分分析(Independent Component Analysis，ICA)是一种基于高阶统计量一种信号分解技术，其目的是将信号分解为若干统计意义上相互独立的成分。由于主成份分析能够较好地分离并消除随机噪声，独立成份分析能够较好地分离并消除生理噪声及各种低频噪声，因此可以先对功能磁共振数据进行主成份分析，利用傅立叶变换对每个成份进行频谱分析，如果某一成分傅立叶频谱的平均功率大于或等于频谱功率的标准偏差，则认为该成份为随机噪声成份，并将该成份置零，然后重建消除随机噪声成份功能磁共振数据，再对消除随机噪声的功能磁共振数据进行独立成分分析，以消除生理噪声及各种低频噪声。

假设采用主成份分析消除系统随机噪声以后，功能磁共振数据中的功能信号、生理噪声以及各种低频噪声是线性混合的，功能磁共振数据模型可进一步用矩阵表示为：

Y＝AC                                                   (2)

其中，A表示未知的线性混合矩阵，C表示包含功能信号和各种噪声成分的矩阵。为了估计功能磁共振数据中的功能信号和噪声成分，可以寻求一个满足(3)的分离矩阵W：

C＝WY                                                   (3)

当分离矩阵W为混合矩阵A的逆矩阵时，功能磁共振数据中的功能信号和噪声成分将能够被分离出来。分离矩阵W可以采用Hyvrinen提出的基于信息论的快速定点算法迭代估计。

利用独立成份分析消除功能磁共振数据中各种生理相关噪声的关键是如何确定生理相关的各种噪声成份。在通常情况下，人体的心脏脉动频率为0.75～1.5Hz，呼吸频率为0.1～0.5Hz。根据Nayquist采样定理，当功能磁共振时间序列图象的采样率高于2×1.5Hz时，通过对分离出的成分进行频谱分析，能够确定心脏和呼吸相关的生理噪声成分。但是，为了覆盖足够的脑空间，功能磁共振常采用多层扫描方式采集数据，扫描重复时间(TR)通常大于或等于2000ms。这样，每一层时间序列图象的采样率小于或等于0.5Hz，低于人体生理运动频率，导致高频生理噪声混叠到低频成分，生理噪声不再保持原来的波形，无法基于频谱分析方法确定多层功能磁共振数据中生理噪声成分。

功能磁共振数据中的噪声主要是由心脏和呼吸相关的非神经性的生理效应引起的。虽然周期性的心脏脉动和呼吸引起的非神经性的生理效应对不同的脑区域影响不同，但是大脑灰质区域以及脑脊液区域同时受周期的心脏和呼吸噪声影响。由于脑脊液数据中不可能有与大脑皮层神经活动相关联的功能信号改变，可以认为脑脊液中的信号只包含生理相关的噪声、各种低频噪声及系统随机噪声。

因此，我们可以基于生理噪声成分时间过程在大脑灰质区域以及脑脊液区域的空间相关性，通过脑脊液数据中生理噪声成分估计并确定灰质数据中生理噪声成分，最终消除多层功能磁共振数据的生理噪声。

因此，采用主成份分析方法分别消除灰质数据、脑脊液数据中的随机噪声后，脑脊液区域每个象素的时间过程y^c(t)可以用(4)表示：

y^c(t)＝p(t)          (t＝1，2，Λ，N)                        (4)

灰质区域每个象素的时间过程y^g(t)可以用(5)表示：

y^g(t)＝s(t)+p(t)     (t＝1，2，Λ，N)                        (5)

采用独立成份分析方法分别对(4)、(5)表示的灰质数据Y^g、脑脊液数据Y^c进行独立成份分解。假设灰质数据中包含m个统计上相互独立的成分c_i ^g(t)，脑脊液数据中包含n个统计上相互独立的成分c_j ^c(t)。为了确定灰质数据中的生理相关的噪声成分，我们利用灰质数据的m个独立成份c_i ^g(t)构造矩阵X_m×N，利用脑脊液数据的n个独立成份c_j ^c(t)构造矩阵Y_m×N。通过对矩阵X_m×N和Y_m×N进行典型相关分析，可得到c_i ^g(t)的加权系数矩阵A_m×r和c_j ^c(t)的加权系数矩阵B_n×r及包含r个元素的行向量 r＝min(m，n)。 的第k个元素为矩阵X_m×N和Y_m×N的各行分别以矩阵A_m×r和B_n×r的第k列元素为加权系数进行加权后所得典型变量的典型相关系数的平方。显然，典型相关系数的平方越大，脑脊液中各成份c_i ^g(t)与灰质中各成分c_j ^c(t)的相关性越大；灰质中成分c_i ^g(t)的加权系数越大，该成分与脑脊液中成份c_j ^c(t)的相关性也越大。因此，可以通过计算行向量 的各元素对矩阵A_m×r各行元素的加权和，得到表达灰质数据中各独立成份c_i ^g(t)与脑脊液数据中各独立成份c_j ^c(t)相关测度的列向量

$\stackrel{\mathrm{\ρ}}{M}=A\×\stackrel{\mathrm{\ρ}}{{\mathrm{\λ}}^T}---\left(6\right)$

对列向量

的各元素从大到小排序，那么，排序后的列向量

的前r个元素对应的独立成分c_i ^g(t)为与脑脊液数据中各独立成份c_j ^c(t)相关最大的噪声成分，将这些成分置零，就可以消除灰质数据中与脑脊液数据中各种噪声成份相关最大的各种噪声成份。对真实的功能磁共振数据试验结果证明了这种方法是很有效的。

Claims (3)

1.一种基于独立成份空间相关性消除功能磁共振数据噪声的方法，其特征在于采用下列步骤：

1)设功能磁共振体数据序列的时间采样点数为N，每个体数据包含S层图像，分别以第i层图像的灰质区域和脑脊液区域为模板，构造该层时间序列图像对应的灰质数据矩阵Y^g和脑脊液数据矩阵Y^c，矩阵Y^g和Y^c的行表达了像素的时间过程；

2)对灰质数据矩阵Y^g进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前p个特征值对应的特征向量的累积方差贡献率大于90％选取灰质数据矩阵Y^g的主成份分量；

3)对脑脊液数据矩阵Y^c进行主成份分析，得到N个按降序排列的特征根及其相应的特征向量，根据前q个特征值对应的特征向量的累积方差贡献率大于90％选取脑脊液数据矩阵Y^c的主成份分量；

4)采用傅立叶频谱分析方法分别对灰质数据矩阵Y^g和脑脊液数据矩阵Y^c的主成分分量进行频谱分析，根据随机噪声判定准则确定并消除灰质数据矩阵Y^g和脑脊液数据矩阵Y^c主成分分量中的随机噪声成分；

5)分别利用步骤4)保留的主成分分量重建消除随机噪声后的灰质数据 和脑脊液数据

6)对灰质数据

进行独立成份分解，得到分离矩阵W_g及m个统计上相互独立的成份c_i ^g(t)，i＝1，2，Λm；

7)对脑脊液数据

进行独立成份分解，得到分离矩阵W_c及n个统计上相互独立的成份c_j ^c(t)，j＝1，2，Λn；

8)利用灰质数据的m个独立成份c_i ^g(t)构造矩阵X_m×N，利用脑脊液数据的n个独立成份c_j ^c(t)构造矩阵Y_n×N；

9)对矩阵X_m×N和Y_n×N进行典型相关分析，得到c_i ^g(t)的加权系数矩阵A_m×r、c_j ^c(t)加权系数矩阵B_n×r及包含r个元素的行向量 r为m和n中较小的一个数，

的第k个元素为矩阵X_m×N和Y_n×N的各行分别以矩阵A_m×r和B_n×r的第k列元素为加权系数加权后得到的典型变量的典型相关系数的平方，k＝1，2，Λr；

10)计算行向量

的各元素对矩阵A_m×r各行元素的加权和，得到表达灰质数据中各独立成份c_i ^g(t)与脑脊液数据中各独立成份c_j ^c(t)相关测度的列向量

11)对列向量

的各元素从大到小排序，那么，排序后的列向量

的前r个元素对应的灰质数据的独立成分c_i ^g(t)为与脑脊液数据中各独立成份c_j ^c(t)相关最大的噪声成分，将这些成分置零后得到一组新的独立成份

12)用分离矩阵W_g的逆矩阵乘以由独立成份

构成的成份矩阵，重建消除与脑脊液数据中各独立成份c_j ^c(t)相关最大的各种噪声成份后的灰质数据；

13)重复步骤1)到12)，直至消除功能磁共振时间序列数据各层灰质数据中的与脑脊液数据中相关最大的各种噪声成分。

2.根据权利要求1所述的基于独立成份空间相关性消除功能磁共振数据噪声的方法，其特征在于通过计算功能实验设计时序与上述步骤1)所述灰质数据矩阵Y^g中每行所表达像素的时间过程的相关系数，用相关系数大于阈值T的像素重新构造灰质数据的矩阵。

3.根据权利要求1所述的基于独立成份空间相关性消除多层功能磁共振数据数据噪声的方法，其特征在于随机噪声判定准则是：如果某一成分傅立叶频谱的平均功率大于或等于频谱功率的标准偏差，则认为该成分为随机噪声成分。

Images