CN116509400A - 基于希尔伯特黄变换的脑磁图动态功能连接构建方法 - Google Patents

Abstract

本发明涉及基于希尔伯特黄变换法的脑磁图动态脑功能连接构建方法。该方法包括：1)脑磁图数据的预处理和源重构，脑区代表信号的选取；2)脑区信号的希尔伯特黄变换和脑区相干矩阵构建、瞬时脑网络构建；3)网络分析，包括瞬时脑网络分析、平均脑网络分析和聚类脑网络分析。采用希尔伯特黄变换元素共轭相乘相干变换构建的瞬时网络连接，完整保留大脑活动固有的非线性信息，且充分利用脑磁图高时间分辨率的特点，在时间上动态显示大脑活动过程，为了解神经疾病发生产生的影响、阐述认知过程、脑部手术进行术前/术后评估、脑机接口探索等领域提供帮助。

Description

基于希尔伯特黄变换的脑磁图动态功能连接构建方法

技术领域

本发明涉及脑磁图功能网络构建技术领域，尤其涉及一种基于静息态脑磁图数据的动态脑功能连接构建方法。

背景技术

人类大脑由大约10^11个神经元细胞组成，神经元之间由约10^15个神经突触相互连接进而形成一个极为复杂的网络。从传统的细胞生物学角度去了解人脑功能、神经疾病原理、大脑发育与老化机制、认知的产生过程是远远不够的。迅速发展的医学影像分析技术证实了脑网络的存在，越来越多的神经解剖和神经生理研究聚焦于脑网络。然而，很多研究仅专注于网络连接的空间特征，忽略了网络的时间维度信息。

大脑活动是大脑响应外界刺激和适应外界活动的动态变化过程，该变化通常发生在毫秒级尺度。脑磁图是一种非侵入式脑功能信号采集设备,其获取的颅外磁信号可反映神经活动发生的位置和时间过程。脑磁图记录信号时间的精度可达毫秒，与大脑处理信息的时间尺度较为一致，因此，脑磁图为探究动态功能脑网络提供了潜在的有效工具。如能在脑磁图时间分辨率下构建并分析动态功能脑网络，将可能推动神经外科手术评估、神经系统疾病发病机制等研究的发展。

脑磁图动态功能连接方法中,以滑窗法为典型代表，该方法描述连续时间窗内功能连接的波动。时间窗内的功能连接构建，通常以信号间存在线性相关为前提，比如计算信号间幅度的相关性、考虑信息流方向的格兰杰因果关系、信号间的相干性等；或者考虑了大脑活动固有的非线性信息,比如相位耦合的相位锁相值、相位滞后指数等。上述方法中，在限定窗宽的时间尺度下,一定程度上丢失了大脑活动随时间变化的细节信息。同时，大量研究表明，大脑震荡产生的节律可在不同频带上形成脑网络，由于窗宽的限制，动态分析的频率范围将受限于窗周期。

以小波变换为代表的时频分析法很好地解决了基于滑窗法进行动态功能连接构建的不足。针对滑窗法中由于窗宽导致的动态分析频率范围受限的问题，小波变换借助丰富的小波基将信号分解至多频段，在研究过程中，我们不仅可以在单频带中评估功能连接，还可扩展至多频带甚至跨频带。在时间分辨率方面，小波相干支持信号在逐样本级别上的动态功能连接构建。

然而，尽管小波变换优于其他常见的时频变换方法(如短时傅里叶变换等)，但其在信号分解方面依赖于小波基函数的选取。由于基函数长度的限制，采用小波相干法构建功能连接时，无法敏锐地捕获信号的快速波动信息，在一定程度上也将导致信号能量谱的泄露。此外，加之小波变换的非自适应性无法体现信号本身的性质，因此，对于非线性非稳定且时空分辨率高的脑磁图信号来说，小波变换并不是最优的分析方法。综上所述，寻找一种能充分体现脑磁图信号本身特点、且不丢失信号信息和分辨率的方法，将在脑磁图动态功能连接构建上具有重要意义。

发明内容

针对现有方法特别是小波变换在脑磁图动态功能连接构建中存在的非自适应性、信息泄露、动态分辨率不足的问题，本发明提出一种基于希尔伯特黄变换的脑磁图动态脑功能连接构建方法。通过自适应信号分解，以瞬时相干来建立动态功能连接，可消除小波变换中基函数的选择对信号的限制，充分发挥脑磁图高时间分辨率的优势。此方法构建的动态功能连接能在时间上保留大脑活动过程的细节，且原理简单计算方便。为实现上述目的，本发明采用如下技术方案：

一种基于希尔伯特黄变换法的脑磁图动态功能连接构建方法，包括以下步骤：

步骤1，获取被试的脑磁图数据及相应被试的MRIT1图像。

步骤2，对原始脑磁图信号进行预处理，采用被试MRIT1图像构建头模型和源模型。

步骤2.1，对原始脑磁图数据进行检查，去除由于传感器故障，导致出现>800fT的高振幅或无振幅的通道，并用邻近通道的平均值/加权平均值来代替；采用0.5HZ高通滤波去除基线漂移，50Hz陷波器去除工频干扰，采用独立成分分析法去除信号中含有的肌电，眼电，心电这些干扰信号；对数据以等时长为间隔进行分割。

步骤2.2，对MRIT1图像进行分割提取脑表面，并建立头的体积传导模型；对MRIT1图像进行切片，构成基于皮质表面信息描述偶极子位置和方向的源模型。

步骤3，基于Beamformer法，对经过预处理的脑磁图信号、头模型、源模型进行源重构，创建通道信号与源的映射关系。

步骤4，根据解剖标记模板将大脑分为116个脑区，每个脑区有数量不等的经过源映射的通道信号。

步骤5，对各个脑区提取功率最大的信号作为脑区代表信号。

步骤6，对脑区代表信号进行希尔伯特黄变换。

步骤6.1，对脑区代表信号进行经验模态分解(EMD)。根据公式X(t)＝MF₁(t)+IMF₂(t)+…+MF_N(t)+(t)，每个代表信号X(t)将自适应分解得到数量不一的本征模态函数IMF和一个残差函数r(t)，本发明中使用EMD需要满足的条件与其他人描述一致，通过计算各个本征模态函数与脑区代表信号的相关性，一致选取相关性较强的前几个本征模态函数进行下一步计算。

步骤6.2，将选取得到的本征模态函数进行希尔伯特变换。以g(t)代表当前进行希尔伯特变换的IMF(t)，希尔伯特变换公式为g(t)的希尔伯特分析变换值z(t)为/>(α(t)代表瞬时幅度，/>代表瞬时相位)。

步骤7，脑磁图动态功能连接构建。

步骤7.1，对经过希尔伯特黄变换的本征模态函数信号矩阵进行元素共轭相乘构建交叉谱/>其中/>和/>表示信号X₁(t)和信号X₂(t)经过希尔伯特黄变换后形成的本征模态函数信号矩阵。

步骤7.2，对交叉谱进行平滑和归一化其中E表示平滑，E′为E的导数，得到相干矩阵。可选地，计算交叉谱各时间点的相位值得到信号间的相位锁相矩阵。

步骤7.3，对任意两两感兴趣区求其相干矩阵/相位锁相矩阵，取出同一时间点的相干值/相位锁相值组成N×N的瞬时相干网络/瞬时相位锁相网络，将各个时间点的瞬时相干网络/瞬时相位锁相网络按时间顺序排列，即得到大脑动态功能连接网络。

与现有技术小波变换相比，本发明的效果是：

我们模拟了一对在时间和频率上具有动态一致性的时间序列，以研究与小波变换相比，我们的方法在脑磁图信号中的有效性。每个时间序列包含两个不同的频率部分，我们将这对时间序列模拟信号命名为1和2。实验中引入了显著性检验，以确保我们的连接性是有效的。在图7(a)中，两个模拟信号具有相似的功率谱分布；在图7(b)中，图谱分为两部分，左列是小波变换的结果，模拟信号1和2的小波谱为前两行所示，第三行表示它们的交叉谱，第四行是95％置信区间的结果。希尔伯特谱在右栏显示，其结果与小波结果相似。比较图7(b)的左右列，其结果显示基于小波变换和希尔伯特黄变换的频谱在时间-频率区域具有相似的能量分布，但希尔伯特黄变换比小波变换提供了更精细的能量细节和更清晰的动态变换。

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

本发明提供了基于脑磁图的动态功能连接构建方法。在基于经过源重构过程的信号基础上，采用希尔伯特黄变换的方法对信号进行时频域转换，在信号上可对脑区进行时域分析、频域分析、时频域分析、相位分析。利用希尔伯特黄变换进行元素共轭相乘相干计算得到瞬时网络，在网络上可进行瞬时网络分析、借助滑窗法进行平均功能网络分析、通过聚类法进行聚类网络分析。根据聚类结果还可进行脑区状态分析等。本发明在网络构建过程中，充分利用了脑磁图高时间分辨率的特点，采用瞬时相干的理念，保留了脑磁图信号间的非线性信息，构建的瞬时功能网络，为揭示脑区间时间进程上的变化提供窗口。

附图说明

图1为本发明数据预处理流程示意图；

图2为希尔伯特黄变换及瞬时功能连接构建示意图；

图3为瞬时功能连接按时间顺序排列构成动态网络示意图；

图4(a)表示健康被试和癫痫患者各IMF分量的希尔伯特加权频率分布情况，图4(b)表示健康被试和癫痫患者各脑区的频谱图；

图5(a)为本发明在真实癫痫患者脑磁图数据上构建的瞬时网络，图5(b)为平均网络，图5(c)聚类网络对比示意图；

图6为本发明在聚类结果上进行的状态转换分析示意图；

图7(a)为本发明模拟信号频谱示意图；

图7(b)为信号经小波变换和希尔伯特黄变换后绘制的信号谱图和信号交叉谱图。

具体实施方式

下面结合附图和具体实施案例对本发明进一步说明，但并不以此作为本发明申请保护范围的限定。

一种基于希尔伯特黄变换法的脑磁图动态脑功能连接构建方法，包括以下步骤：

如图1所示，为本发明的预处理过程，包括以下步骤：

步骤1：数据准备阶段，需要获得干净的脑磁图通道数据，头模型，源模型。在本实施例中，需要对采集到的脑磁图通道数据和MRIT1图像进行处理。

脑磁图通道数据的预处理过程包括：

1)去除损坏通道：将没有采集到信号的通道、信号幅度过大的异常通道用临近通道信号幅度的平均值或加权平均值进行代替。

2)滤波：将脑磁图数据通过截止频率为0.5Hz的高通滤波器去除基线漂移，50Hz陷波器去除工频干扰。可选地，根据信号分析频率范围选用相应的低通滤波器，根据数据分辨率要求合理进行下采样操作。

3)去噪：采用独立成分分析法对脑磁图数据中含有的心电、眼电、肌电干扰信号进行去除。

4)信号分割：根据所需，将信号分割成2s一段

MRIT1图像预处理过程包括：

1)头模型创建：将MRI分割为灰质、白质、脑脊液等不同的组织类型，基于不同组织，创建头骨及其内容物的二进制掩码，即创建头模型。

2)源模型创建：将MRI进行切片，使用的厚度分辨率为1mm，维度为256×256×256，对MRI图形进体积处理，创建皮质网格，即源模型。

步骤2：源重构阶段，创建通道信号与源的映射关系。在本实施例中，采用Beamformer法，创建脑磁图通道信号与源位置的映射关系。

步骤3：感兴趣脑区信号提取阶段，提取相应脑区代表信号。在本实施例中，采用116脑区的解剖标记模板对大脑进行分区，每个区域内有数量不等的通道信号，采用最大功率法将脑区内功率最大的通道信号作为脑区的代表信号。

如图2所示，为本发明瞬时功能网络构建过程，包括以下步骤：

步骤4：对提取的感兴趣脑区信号进行分割，分割的信号片段进行希尔伯特变换的第一步，经验模态分解(EMD)。信号经经验模态分解后得到数量不一的本征模态函数(IMF)，经计算本征模态函数与原信号的相关关系，在本实施例中，选取前5个本征模态函数进行希尔伯特黄变换的第二步，希尔伯特变换。

步骤5：感兴趣脑区信号X₁(t)和信号X₂(t)经希尔伯特黄变换得到信号矩阵和/>进行元素共轭相乘构建交叉谱/>对交叉谱进行平滑和归一化得到相干矩阵，相干矩阵的横坐标表示时间，纵坐标表示本征模态函数的编号。本实施例中采用默认脑网络的22个脑区作为感兴趣脑区，故可得到484个相干矩阵。

步骤6：选取484个相干矩阵中同一时间点的相干值，构成22×22的脑区瞬时功能连接网络。

可选地，计算交叉谱各时间点的相位值得到信号间的相位锁相矩阵，提取同一时间点的相位锁相值得到22×22的脑区瞬时相位锁相连接矩阵。

如图3所示，为本发明动态网络的直观体现。将各个时间点的瞬时功能连接网络按时间顺序排列，即可动态显示网络变化。

如图4所示，从被试样本中随机选取3名被试，健康被试、左侧颞叶癫痫被试与右侧颞叶癫痫被试各一名。采用希尔伯特黄变换对脑区信号进行处理，4(a)体现出各被试IMF分量的希尔伯特加权频率(HWF)差异，4(b)体现出各被试在不同脑区中的能量差异。

如图5所示，为三种不同时间尺度下的脑功能网络的可视化比较。5(a)表示瞬时脑功能网络，5(b)表示窗宽大小为50个时间点下的平均脑功能网络，5(c)表示瞬时脑功能网络经过聚类后得到的聚类网络之一。

如图6所示，为聚类系数k＝5得到的5个状态在时间上的变化。

综上，本发明提供的基于希尔伯特黄变换法的脑磁图动态脑功能连接构建方法，进而进行的网络构建和分析。包括信号处理阶段的时域分析，频域分析，时频域分析，相位分析，能量分析；在网络分析中的瞬时网络分析，基于滑窗法的平均网络分析，基于聚类法的聚类网络分析，基于聚类的网络的状态分析等。

需要理解的是，以上对本发明的具体实施例进行的描述只是为了说明本发明的技术路线和特点，其目的在于让本领域的技术人员能够理解本发明并实施该技术，但本发明并不限于上述的特定实施方式。凡是利用本发明说明书及附图内容所做的等效流程变换，直接或间接运用在其他相关技术领域，均在本发明的专利保护范围内。

Claims (1)

1.一种基于希尔伯特黄变换法的脑磁图动态功能连接构建方法，其特征在于，包括以下步骤：

步骤1，获取被试的脑磁图数据及相应被试的MRIT1图像；

步骤2，对原始脑磁图信号进行预处理，采用被试MRIT1图像构建头模型和源模型；

步骤2.1，对原始脑磁图数据进行检查，去除由于传感器故障，导致出现>800fT的高振幅或无振幅的通道，并用邻近通道的平均值/加权平均值来代替；采用0.5HZ高通滤波去除基线漂移，50Hz陷波器去除工频干扰，采用独立成分分析法去除信号中含有的肌电，眼电，心电这些干扰信号；对数据以等时长为间隔进行分割；

步骤2.2，对MRIT1图像进行分割提取脑表面，并建立头的体积传导模型；对MRIT1图像进行切片，构成基于皮质表面信息描述偶极子位置和方向的源模型；

步骤3，基于Beamformer法，对经过预处理的脑磁图信号、头模型、源模型进行源重构，创建通道信号与源的映射关系；

步骤4，根据解剖标记模板将大脑分为116个脑区，每个脑区有数量不等的经过源映射的通道信号；

步骤5，对各个脑区提取功率最大的信号作为脑区代表信号；

步骤6，对脑区代表信号进行希尔伯特黄变换；

步骤6.1，对脑区代表信号进行经验模态分解(EMD)；根据公式X(t)＝IMF₁(t)+IMF₂(t)+…+IMF_N(t)+r(t)，每个代表信号X(t)将自适应分解得到数量不一的本征模态函数IMF和一个残差函数r(t)，通过计算各个本征模态函数与脑区代表信号的相关性，选取相关性较强的前几个本征模态函数进行下一步计算；

步骤6.2，将选取得到的本征模态函数进行希尔伯特变换；以g(t)代表当前进行希尔伯特变换的IMF(t)，希尔伯特变换公式为g(t)的希尔伯特分析变换值z(t)为/>其中α(t)代表瞬时幅度，/>代表瞬时相位；

步骤7，脑磁图动态功能连接构建；

步骤7.1，对经过希尔伯特黄变换的本征模态函数信号矩阵进行元素共轭相乘构建交叉谱/>其中/>和/>表示信号X₁(t)和信号X₂(t)经过希尔伯特黄变换后形成的本征模态函数信号矩阵；

步骤7.2，对交叉谱进行平滑和归一化其中E表示平滑，E^′为E的导数，得到相干矩阵；计算交叉谱各时间点的相位值得到信号间的相位锁相矩阵；

步骤7.3，对任意两两感兴趣区求其相干矩阵/相位锁相矩阵，取出同一时间点的相干值/相位锁相值组成N×N的瞬时相干网络/瞬时相位锁相网络，将各个时间点的瞬时相干网络/瞬时相位锁相网络按时间顺序排列，即得到大脑动态功能连接网络。