CN108596879A

CN108596879A - 一种基于希尔伯特黄变换的fMRI时频域动态网络构建方法

Info

Publication number: CN108596879A
Application number: CN201810281273.6A
Authority: CN
Inventors: 张馨; 杨春兰; 吴水才
Original assignee: Beijing University of Technology
Current assignee: Beijing University of Technology
Priority date: 2018-04-02
Filing date: 2018-04-02
Publication date: 2018-09-28
Anticipated expiration: 2038-04-02
Also published as: CN108596879B

Abstract

本发明公开一种基于希尔伯特黄变换的fMRI时频域动态网络构建方法，包括：输入原始图像；对原始图像进行预处理；提取基于全脑强功能连接的脑区及其对应的时间序列；时间序列后处理；fMRI的HHT的时频域动态网络构建；在时频域上对网络的变化模型进行分析研究。本发明采用希尔伯特黄变换算法可以根据数据本身自适应生成基函数，避免提前选取不合适的基函数产生的不准确结果。

Description

一种基于希尔伯特黄变换的fMRI时频域动态网络构建方法

技术领域

本发明属于医学图像处理领域，涉及一种静息态下fMRI(functional MagneticResonance Image,功能磁共振成像)的时频域动态网络算法,主要采用希尔伯特黄变换算法进行时频域动态网络构建。

背景技术

目前为止，大部分基于fMRI的功能连接研究假设，远距离大脑区域间信号在统计学上的相互依赖模式是固定不变的，比如相关、协方差和不同区域时间序列的相互作用等，这一状态贯穿记录静息态实验的整个周期。在这一假设下，研究大规模尺度下的大脑功能已经取得了卓越的成果，这些结果描述了复杂时空的平均现象。

然而，人类的大脑连接最有可能具有动态性和时间依赖性，并与正在进行的活动有关。因此，有人提出研究功能连接随时间变化的动态特征可能会更好的发现大脑网络的基本特性。

功能连接中检测动态性最常用的方法是滑窗分析算法。滑窗分析的主要限制是使用了固定窗长。窗的大小影响着时间序列分析的时间尺度，理想情况下，窗长要足够长以适应fMRI测得的BOLD(Blood oxygenation level dependent，血氧水平依赖)信号相对低频的的部分，同时，窗长还要足够短以便捕捉网络连接中转瞬即逝的变化。然而，神经元的相关频率和连接变化的合适时间尺度仍是待解决的问题。

时频分析可用来评估两个信号的相干和相位延迟。基于小波变换的时频域相干分析可以提供多分辨率算法进行时频分析，摆脱了时间滑窗分析中固定窗宽的限制。小波变换中，有效分析窗的大小随信号频率本身的时间尺度而改变：高频用短窗分析，低频采用长窗分析。通过相干提供的多重时间尺度，小波变换用于动态网络构建可提供更加灵活的探索分析。例如，可以呈现出脑网络所分布的主要频域范围，也可以在特定频率段研究网络中节点间连接强度随时间的改变和与相位的关系。在时频域动态网络构建中，目前均采用小波变换算法。基于小波变换的相干方法虽然可以同时捕获高频和低频信息，但需要提前选取基函数，细节信息的捕获不够精确。基函数选取对数据分析的结果影响较大，然而基函数的选取一般要凭经验确定。BOLD信号是时变非平稳信号，小波变换分析时变非平稳信号不是最优选择。因此，构建fMRI时频域动态网络时需要一种自适应选取基函数的算法来代替小波变换。

发明内容

针对小波相干在fMRI时频域动态网络构建中的一些不足，如在时频域分析方法中不是最适合分析时变非平稳信号的方法,不能自适应选取基函数和不能精确捕获动态网络细节信息，本发明提出采用希尔伯特黄变换(Hilbert-Huang transform，HHT)取代小波变换构建时频域动态网络。HHT可以根据数据特征自适应生成相应的基函数，避免主观选取基函数导致结果不准确。HHT相对小波变换更适合分析时变非平稳信号，理论上可以精确提取动态网络的细节信息，捕获更精确的网络特征。HHT在算法参数设置上更简洁。

为实现上述目的，本发明采用如下技术方案：首先对原始图像进行常规fMRI预处理，主要为时间层矫正、去除头动和伪影、配准到结构像、标准化、平滑、滤波、去除生理噪声等；然后进行群组独立成分分析(Group independent component analysis，GICA)，提取强连接的大脑区域及其对应的时间序列进行后续网络构建；最后对时间序列进行后处理进一步降低噪声对网络动态性的干扰。将经过后处理的时间序列由HHT中的经验模态分解(Empirical mode decomposition,EMD)算法分解，得到本征模态函数(Intrinsic modefunction,IMF)。将得到的IMF采用Hilbert变换映射到时频域上，得到每个时频点的HHT系数。接着，对HHT系数进行相干计算。最后对所有被试的每个时频点对应的相干连接矩阵进行K-means聚类得到5个不同的动态网络连接状态。

本发明采取以下步骤：

步骤1，输入原始图像。

步骤2，对原始图像进行预处理。

对fMRI数据进行预处理，主要为时间层矫正、去除头动和伪影、配准到结构像、标准化、平滑、滤波、去除生理噪声等。

步骤3，在全脑范围提取基于强功能连接的脑区以及与其对应的时间序列。

使用GICA软件筛选出大脑连接较强的区域及对应的时间序列，以便进行后续动态网络连接构建。首先，使用主成分分析(Principal component analysis，PCA)提取150个主成分，接着再次降维，进行第二次PCA分析选出100个主成分。使用Infomax ICA算法提取100个独立成分，GICA进行后重建，具体的参数设置如图1。然后，在100个独立成分中人工选取具有功能活动的大脑区域，排除噪声成分，得到具有强连接的脑区及其时间序列。上述选取依据为激活区域在灰质皮层并且时间序列的频率范围集中在低频区域的部分，如图2。

步骤4，基于HHT构建fMRI时频域动态网络。

本发明采用的方法如下：

步骤4.1计算时频域上的HHT系数

HHT主要包括两部分：1)执行EMD算法；2)Hilbert谱分析(Hilbert spectrumanalysis,HSA)。EMD是为了获得IMF,它具有自适应性、正交性、完备性等特性。使用EMD须满足如下两个条件：1)信号极值点的数量与零点数相等或相差1；2)由极大值定义的上包络和由极小值定义的下包络的信号局部均值为0。EMD过程如下：

1)对输入信号求取所有的极大值点和极小值点。

2)对极大值点和极小值点采用三次样条进行拟合，求上、下包络的曲线，计算均值函数，进而求出信号和其均值的差值h。

3)考察h是否满足IMF条件，如果满足，则把h作为第1个IMF；否则，对其进行前两步操作，直到第k步满足IMF条件，然后求得第1个IMF，求出原信号与IMF的差值r。

4)把差值r作为待分解的信号，直到r为单调信号或者只存在一个极点为止，得到的表达式如下：

其中，S(t)为原信号，C_i(t)表示第i次筛选得到的IMF分量，N为筛选次数，R_n(t)为最终的剩余分量。执行EMD过程后，进行Hilbert谱分析，下式对每个IMF分量进行Hilbert谱变换：

解析信号为：

进而，由式(4)和式(5)分别求得瞬时幅值和瞬时相位：

首先，将后处理得到的时间序列由HHT中的EMD算法进行分解，得到IMFs。去除残差，进行希尔伯特变换，将IMFs变换到时频域上，得到每个时频点的HHT系数。

步骤4.2计算时频域相干矩阵

对HHT系数进行相干计算，首先要构建交叉希尔伯特黄变换，它是希尔伯特黄变换后系数间的元素共轭相乘，公式如下：

W^xy＝W^x*W^y (6)

W^x和W^y是时间序列X和Y经过HHT得到的，点乘代表元素共轭乘法。上面的W^x、W^y和W^xy测量结果还需要由信号频谱标准化，目的是相干性估计不会偏向于有更多能量的信号部分。另外，在标准化度量上引入了平滑函数，以避免团块化。这种标准化和平滑的测量被称为相干变换，其定义如下：

最后，对每个被试的每个时频点都构建一个矩阵，可以得到大小为被试个体数×频率点×时间点个矩阵，这就构建出了时频域动态网络，从而方便进行网络的动态特性研究。

为了做对比分析，本发明还进行了复Morlet小波变换的时频域动态网络构建。首先，采用复Morlet小波对之前提取的脑区的时间序列进行分解，得到时频域上的复Morlet小波系数。接着，对小波系数进行相干计算，在计算相干中还涉及到标准化和平滑，采用信号频谱标准化和移动加权平均法进行平滑。对每个被试的每个时频点都构建出一个矩阵，所以可以得到被试个体数×频率点×时间点个矩阵，这就构建出了时频域动态网络，可以研究网络的动态特性。

步骤5，在时频域上对网络的变化状态进行分析研究。

因为得到了被试个体数×频率点×时间点个矩阵，数据量很大，研究起来不方便，而且有研究发现在动态网络中存在着重复的连接状态。因此，建立区域间可重复性、瞬态的相关连接状态以解决这一问题。可以采用聚类实现不同连接状态研究。聚类算法提供了一个潜在有力的方法来检测被试内部状态的自发改变。对所有被试的每个时频点所对应的相干连接矩阵进行聚类，得到几个不同的动态功能连接状态；在时频域上研究这些状态下的网络连接与各个状态在矩阵上的占有率和相位分布的差异等，这有助于更好的理解大脑网络连接。此外，也可以在特定的频率上精确研究时间序列的动态网络连接状态。

与现有技术相比，本发明的有益效果是：

1.采用自适应生成基函数的时频域分析，HHT算法避免了小波变换需主观选取基函数的缺点，理论上可以提取更加精准的时频域动态网络细节信息，更准确的传达出大脑网络的连接状态。

2.基于HHT的fMRI时频域动态网络构建过程中，参数设定相比小波变换更加简单，可以减少计算量。

附图说明

图1GICA参数设置表示图；

图2ICA成分选取图；

图3复Morlet小波变换聚类后具有重复性的5个时频域动态网络连接状态，其中，状态1为40.76％，状态2为24.69％，状态3为22.19％，状态4为8.37％，状态5为：4.02％；

图4HHT聚类后具有重复性的5个时频域动态网络连接状态，其中，状态1为43.35％，状态2为23.97％，状态3为20.82％，状态4为7.36％，状态5为：4.5％；

图5为本发明流程图。

具体实施方式

下面结合附图对本发明作进一步的详细说明。

如图5所示，本发明包括以下步骤：

步骤1，输入原始图像。

步骤2，对原始图像进行预处理。

采用DPARSF软件对fMRI数据进行预处理，主要为时间层矫正，去除头动和伪影、配准到结构像、标准化、平滑、滤波、去除生理噪声等。

步骤3，使用GICA提取基于全脑强功能连接脑区及时间序列。

图1GICA参数设置的表示图，提供GICA分析的详细参数设定；

步骤4，基于HHT的fMRI时频域动态网络构建。

步骤5，在时频域上对网络的变化状态进行分析研究。

下面通过实验数据比较本发明所述方法和复Morlet小波变换构建时频域动态网络方法效果的优劣。

实验数据来自Human Connectome Project(HCP)数据库，选取40例年龄在26～35岁之间的健康被试(男19例，女21例)。被试在静息状态下使用3T西门子核磁共振扫描仪(3TSiemens Connectome Skyra MRI)梯度回波平面回波成像序列进行扫描，时间为15分钟。

图3为复Morlet小波相干矩阵进行聚类得到的5类时频域动态网络连接状态，饱和度表示连接幅值，状态1-5分别表示聚类后得到的5个重复的时频域动态网络连接状态，每个状态图中还包含所有连接矩阵在此状态上的占有率、频域分布直方图和相位信息。时频域动态网络连接状态总结了大脑功能连接的重复状态。并且在每个状态上，可以观察到脑区间的连接情况、连接强度、连接状态主要集中在哪些频率范围内，同时还可观察相应信息。

图4为HHT相干矩阵进行聚类得到5类时频域动态网络连接状态，饱和度显示连接幅值，状态1-5分别表示聚类后得到的5个重复的时频域动态网络连接状态，每个状态图中还包含所有连接矩阵在此状态上的占有率、频域分布直方图和相位信息。时频域动态网络连接状态总结了大脑功能连接的重复状态。并且在每个状态上可以观察到脑区间的连接情况和连接强度和此连接状态主要集中在哪些频率范围内，同时还呈现出相位的分布情况。与复Morlet小波构建时频域动态网络连接状态相比较，HHT构建的状态中脑区连接状态与频率相位分布与图3复Morlet小波构建的较为相似，两种方法构建的5个状态对应的矩阵占有率也几乎相同。图3、4中矩阵占有率最高的状态的相位分布主要都在0附近，随着占有率的降低，状态中的相位分布逐渐分散，并且状态2和3对应的频率分布直方图也很一致。因此，证明HHT确实适合构建时频域动态网络。

本发明的采用希尔伯特黄变换取代小波变换构建时频域动态网络方法，包括：输入原始图像；对原始图像进行预处理；提取基于全脑强功能连接的脑区及其对应的时间序列；时间序列后处理；fMRI的HHT的时频域动态网络构建；在时频域上对网络的变化模型进行分析研究。本发明采用希尔伯特黄变换算法可以根据数据本身自适应生成基函数，避免提前选取不合适的基函数产生的不准确结果。HHT适合于分析时变非平稳信号，理论上可以精确提取动态网络的细节信息，捕获更精确的网络特征。同时在算法参数设置上更简洁。实验结果表明，与复Morlet小波构建时频域动态网络相比较，HHT构建的状态中脑区连接状态与频率相位分布与复Morlet小波构建的较为相似，两种方法构建的5个状态对应的矩阵占有率也几乎相同，证明了HHT确实适合构建时频域动态网络。

Claims

1.一种基于希尔伯特黄变换的fMRI时频域动态网络构建方法，其特征在于包括以下步骤：

步骤1，输入原始图像；

步骤2，对原始图像进行预处理；

使用GICA软件筛选出大脑连接较强的区域及对应的时间序列，以便进行后续动态网络连接构建；首先，使用主成分分析(Principal component analysis，PCA)提取150个主成分，接着再次降维，进行第二次PCA分析选出100个主成分；使用Infomax ICA算法提取100个独立成分，GICA进行后重建，具体的参数设置如图1；然后，在100个独立成分中选取具有功能活动的大脑区域，排除噪声成分，得到具有强连接的脑区及其时间序列；上述选取依据为激活区域在灰质皮层并且时间序列的频率范围集中在低频区域的部分；

步骤4，基于HHT构建fMRI时频域动态网络。

2.如权利要求1所述的基于希尔伯特黄变换的fMRI时频域动态网络构建方法，其特征在于，步骤4包括以下步骤：

步骤4.1计算时频域上的HHT系数

HHT包括两部分：1)执行EMD算法；2)Hilbert谱分析(Hilbert spectrum analysis,HSA)，使用EMD须满足如下两个条件：1)信号极值点的数量与零点数相等或相差1；2)由极大值定义的上包络和由极小值定义的下包络的信号局部均值为0，EMD过程如下：

1)对输入信号求取所有的极大值点和极小值点；

2)对极大值点和极小值点采用三次样条进行拟合，求上、下包络的曲线，计算均值函数，进而求出信号和其均值的差值h；

3)考察h是否满足IMF条件，如果满足，则把h作为第1个IMF；否则，对其进行前两步操作，直到第k步满足IMF条件，然后求得第1个IMF，求出原信号与IMF的差值r；

其中，S(t)为原信号，C_i(t)表示第i次筛选得到的IMF分量，N为筛选次数，R_n(t)为最终的剩余分量；

执行EMD过程后，进行Hilbert谱分析，下式对每个IMF分量进行Hilbert谱变换：

解析信号为：

进而，由式(4)和式(5)分别求得瞬时幅值和瞬时相位：

首先，将后处理得到的时间序列由HHT中的EMD算法进行分解，得到IMFs，去除残差，进行希尔伯特变换，将IMFs变换到时频域上，得到每个时频点的HHT系数；

步骤4.2计算时频域相干矩阵

W^xy＝W^x*W^y (6)

其中，W^x和W^y是时间序列X和Y经过HHT得到的，点乘*代表元素共轭乘法，所述W^x、W^y和W^xy测量结果还需要由信号频谱标准化，另外，在标准化度量上引入了平滑函数，以避免团块化，这种标准化和平滑的测量被称为相干变换，其定义如下：

最后，对每个被试的每个时频点都构建一个矩阵，可以得到大小为被试个体数×频率点×时间点的矩阵，构建出了时频域动态网络。