CN112183477A - 一种基于持续同调的脑电信号持续特征提取方法 - Google Patents

Abstract

本发明公开一种基于持续同调的脑电信号持续特征提取方法，包括以下步骤：S1、获取大脑成像数据并进行预处理，得到EEG时间序列信号，确定大脑的网络节点；S2、将网络节点映射到有限度量空间构造邻接矩阵，度量空间采用皮尔逊相关性度量，将各网络节点间皮尔逊相关系数的倒数作为网络节点间连接的权重；本发明基于持续同调性的全尺度脑网络分析模型，并对模型中的每一数据处理步骤涉及到的算法和参数做了分析，将模型应用于精神分裂症WM持续拓扑特征分析中，在获取持续拓扑特征的基础上，本发明提出的模型具有稳定性和对抗噪性，可以为精神类疾病医学影像分析提供稳定生物参考指标。

Description

一种基于持续同调的脑电信号持续特征提取方法

技术领域

本发明属于脑电技术领域，尤其涉及一种基于持续同调的脑电信号持续特征提取方法。

背景技术

宇航员面临着我们大多数人永远都不会经历的精神和心理挑战，自从美国女宇航员莉萨·诺瓦克事件后，NASA一直致力于减少太空飞行期间任何精神健康问题的风险问题的研究。因为在载人航天面临的重大挑战中，航天员需面临的心理学和精神病学的挑战是最大的制约因素。为此，本发明将从宇航员的视角，运用拓扑数据分析中的持续同调理论，以精神分裂症数据为例，分析其复杂脑网络的持续性特征，以便能够为提前预防和检测宇航员精神疾病提供有效的脑网络指标。

发明内容

本发明提供一种基于持续同调的脑电信号持续特征提取方法，旨在解决上述存在的问题。

本发明是这样实现的，一种基于持续同调的脑电信号持续特征提取方法，包括以下步骤：

S1、获取大脑成像数据并进行预处理，得到EEG时间序列信号，确定大脑的网络节点；

S2、将网络节点映射到有限度量空间构造邻接矩阵，度量空间采用皮尔逊相关性度量，将各网络节点间皮尔逊相关系数的倒数作为网络节点间连接的权重；

S3、通过邻接矩阵构造边权矩阵，选择过滤阈值和过滤算法对边权矩阵进行过滤，构造脑网络复形；

S4、可视化表示过滤过程中脑网络的拓扑特征变化，并获取持续拓扑特征；

S5、分析持续拓扑特征的稳定性，并输出脑网络的持续不变特征。

进一步的，在步骤S4中，通过条形码或者持续性图可视化表示过滤过程中脑网络的拓扑特征变化。

进一步的，在过滤过程中，计算p维贝蒂数间隔[ε_birth，ε_death]；其中，ε_birth为单纯复形中p维孔开始的时间，ε_death为单纯复形中p维孔消失的时间。

进一步的，当通过持续性图可视化表示过滤过程中脑网络的拓扑特征变化时，持续性图中横坐标表示持续特征出现的时间ε_birth，纵坐标表示持续特征消失的的时间ε_death，将过滤过程中求得的间隔集合[ε_birth，ε_death]作为持续性图中点的坐标，将集合中所有的间隔对表示在坐标中。

进一步的，当通过条形码可视化表示过滤过程中脑网络的拓扑特征变化时，条形码中横坐标表示过滤阈值ε，条形码的长度表示[ε_birth，ε_death]的长度。

进一步的，在步骤S5中，采用Bottleneck距离和Wasserstein距离来分析持续拓扑特征的稳定性。

进一步的，两个图X和Y之间的p阶Wassertein距离定义为：

其中，p∈|1，∞)，

是从X到Y的映射；

当p＝∞时，距离d是二维空间的度量，以上公式表示为：

其中，Bottleneck距离为W_∞[d_∞]。

拓扑数据分析(TDA)是一个数据分析、代数拓扑、计算几何、计算机科学、统计等多领域相关的一个领域。TDA的主要目标是利用几何学和拓扑学的思想和成果来开发研究数据的定性特征一个工具。为了达到这个目标，TDA需要精确的定义定性特性，还有在具体实践应用中的计算工具，以及保证这些特性稳定、健壮性的理论。解决这三个问题的一种方法就是TDA中的持久同源性(PH)。

对大脑成像数据进行处理和分析时，通过生成表示节点之间连接强度的矩阵来构建脑网络。对于解剖学连通性数据，连通性矩阵可以表示大脑区域之间纤维束的概率，而在功能连通性数据中，连通性矩阵可以只是来自不同体素/传感器的时程之间的相关矩阵。然后通常将阈值应用于此矩阵以生成邻接矩阵，阈值的选择在网络构建中起着重要作用，因为它影响连接的密度并影响网络拓扑。同样，在比较网络时，阈值设置也很重要。

与现有技术相比，本发明的有益效果是：本发明基于持续同调性的全尺度脑网络分析模型，并对模型中的每一数据处理步骤涉及到的算法和参数做了分析，研究了全尺度脑网络构造中的节点、边权矩阵构造、过滤阈值的选择等关键性问题。将模型应用于精神分裂症WM持续拓扑特征分析中，在获取持续拓扑特征的基础上，实验结果还表明本文提出的模型具有稳定性和对抗噪性，可以为精神类疾病医学影像分析提供稳定生物参考指标。本发明不但可以解决构建复杂脑网络时黄金阈值选择难的问题，而且根据尺度的变化可动态观测每一个有效尺度下网络拓扑结构的变化。

附图说明

图1为本发明中基于PH的全尺度脑网络分析模型框架；

图2为本发明编码阶段不同连接密度下邻接矩阵，(a)为健康被试，(b)为精分病人；

图3为本发明编码阶段不同频段及全频段邻接矩阵，(a)为健康被试，(b)为精分病人；

图4为健康被试编码阶段持续图(Fs＝20)；

图5为健康被试编码阶段持续图(Fs＝100)；

图6为健康被试编码阶段持续图(Fs＝1000)；

图7为健康被试编码阶段条形码；

图8为精分病人编码阶段条形码；

图9为精分病人编码阶段持续图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明的描述中，需要理解的是，术语“长度”、“宽度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。

实施例

根据持续同源数据分析流程，结合脑电信号处理设计的全尺度脑网络分析模型如图1所示。

模型的输入是脑电时间序列信号，选取合适的度量空间后，脑电时序信号中的点就称为点云，在该空间构建点云的邻接矩阵，然后通过计算持续同调性获取持续拓扑特征，最后通过持续特征的稳定性分析确定网络的持续不变特征。

模型的输入是经过预处理的EEG时间序列信号。

在一定的有限度量空间下，为EEG信号各通道数据(即点云)构造连接矩阵，度量空间可以选择距离度量，相关性度量或者同步性度量。连接矩阵可以是无向无权，也可以是无向加权网。由于依赖的技术背景和系统稳定性问题，基于持续同调的复杂脑网络分析不讨论有向加权网。

如果输入的数据是数值型的，一般采用欧拉距离，当然也有其他距离如闵科夫斯基距离作为度量空间。对于本模型数据模型，因为输入主要是脑电信号，所以根据以往的研究可以采用相关性或者同步性度量节点间的距离。

选择确定的度量空间后，就可以构造网络节点间的邻接矩阵。如果构造的是无向无权网，网络节点可以直接作为一维的单纯形，而如果是无向加权网，则需要通过增加或者减少权重进行过滤构造一维单纯复形。

构造大脑网络嵌套复形的过程也就是计算持续同调的过程。持续同调分两部分，即同源性和持续性。同源在群论称作同调，它是拓扑集合分类的工具，可以度量一个单纯复形的特定结构；持续性是指给定一个，在所有可能的值下计算哪些结构是持续存在的，即获得持续拓扑特征。特征中能够保持时间长的就是有用的特征，而寿命较短的可能是噪声，这个过程就称为持续同调。构造复形的关键步骤就是选择合适的过滤算法和过滤阈值。

过滤阈值的选择非常重要，一般的方法是通过选择不同的构造复形，找到有效的结果所对应的。如果过小，那么构造的复形就可能是原始的点云，或者点云的几条边；如果过大，可能的结果就是原始点云构成一个巨大的超维复形。

针对不同实际的应用，需要构造不同类型的单纯复形，这些复形结构具有不同的属性，有的容易用数学描述，有的容易计算，有的简单但计算效率低。最常见构造单纯复形的算法有：

ech complexu算法，Vietoris–Rips complex算法，alpha complex算法，witness complex算法。其中Vietoris–Rips complex算法和witness complex算法是基于图的过滤而不是基于覆盖的过滤，所以非常适合应用于基于图论的复杂脑网络复形构造，并且对处理高维数据有很好的性能。尤其是Vietoris–Rips complex算法从计算的角度来看，描述和实际操作相对容易，所以应用较多，本实验也主要是采用Vietoris–Ripscomplex算法过滤。

随着过滤阈值的变化，Rips复形的拓扑特征会发生变化。过滤过程中网络的拓扑变化可以使用条形码或者持续性图可视化表示。

在过滤过程中主要是计算p维贝蒂数间隔[ε_birth，ε_death]，其中，ε_birth就是单纯复形中p维孔开始的时间，而ε_death是它消失的时间，同时它们也是条形码中一个条形码的开始和结束点，这些间隔用图表示出来就是持续条形码。与条形码等效的就是持续图。

持续图中，横坐标表示持续特征出现的时间，即ε_birth，纵坐标是持续特征消失的的时间ε_death，将过滤过程中求得的间隔集合[ε_birth，ε_death]作为持续图中点的坐标，将集合中所有的间隔对表示在坐标中就绘制了持续图。

在条形码中，横坐标表示过滤阈值ε，用[ε_birth，ε_death]的长度表示条形码的长度，条形码长度大的表示持续拓扑特征，条形码长度很短或者只是一个点的表示噪声，相对应的在持续图中，离对角线远的点表示持续特征，离对角线很近的点表示噪声。

拓扑特征的稳定性分析也就是条形码的统计分析问题，是目前快速发展的一个研究方向，也就是开发相应的统计方法使用持续图比较对结果做出相应的解释。在网络匹配问题中，为每个网络创建一个持久性图，并对持久性图进行比较，来获得网络的相似性。

解决条形码的统计分析问题需要使用稳定的度量空间来度量持久图之间的距离。如果对数据集的一个小扰动只会在该指标之前的持久性关系图中造成一个小的变动，就可以说明该指标是稳定的。通常，有两个稳定性度量标准(Bottleneck距离和Wasserstein距离)用于度量图之间的距离Bottleneck距离和Wasserstein距离。

定义：两个图X和Y之间的p阶Wassertein距离定义为：

其中，p∈[1，∞)，：X→Y是从X到Y的映射。当p＝∞时，距离d是二维空间的度量,以上公式表示为:

W_∞[d_∞]就是Bottleneck距离。

Bottleneck距离度量的是两个图对应匹配点之间的最大距离，也就是说只输出最大离群值之间的距离，因此它可以捕获持续图大的变化。而Wasserstein距离度量的是两个图对应匹配点之间的总距离，因此可以提供持久性图之间相似性的总体量化。它对持续图小的变化比较敏感。通俗的说Bottleneck距离捕捉匹配点极端情况的距离最大值，Wasserstein距离侧重于全部对应匹配点的平均表现。本发明选取这两个距离来分析所求的持续拓扑特征的稳定性。

试验例

实验数据及预处理

实验数据采用精神分裂症工作记忆EEG信号，实验数据来自北京市回龙观医院，经过重参考、分段、去除眼电、肌电等伪迹等预处理，网络节点规模为60个，分五个波段θ(4～7Hz)、α(7～14Hz)、β1(14～20Hz)、β2(20～30Hz)、γ(30～40Hz)。实验分三个阶段，每个被试每个频段每个阶段下，选取20个平稳脑电信号拼接形成最终的EEG时间序列。

基于PH的邻接矩阵构建

本实验将网络节点映射到有限度量空间构造邻接矩阵，度量空间采用的是皮尔逊相关性度量。这样将各节点间皮尔逊相关系数的倒数作为节点间连的权重，得到一个60*60无向加权网。作为模型的输入。具体为：如果节点i和j连接边的权重为ω，则i和j之间的距离为1/ω,否则的话，我们定义节点i和j的距离为两个节点间所有距离的最小值。精神分裂症工作记忆编码阶段不同稀疏度下构造的动态邻接矩阵如图2所示。

图3为无阈值条件下构造的边权矩阵，从图3不难发现：

(1)不同的连接密度(阈值)，脑网络构造结果不同；

(2)网络连接密度较小时(20％左右)，健康被试和精分病人脑网络差异很大；

(3)从网络连接密度50％左右开始，健康被试和精分病人的连接矩阵变化逐渐变小。

也就是说健康被试和精分病人在工作记忆编码阶段的连接矩阵是有“形状”的，同时从图3无阈值邻接矩阵也可以直观的看出脑网络是有形状的，通过计算持续同调性来找到其内在的形状(持续拓扑特征)。

构造精神分裂症WM数据的单纯复形

实验采用的是由斯坦福大学拓扑计算小组基于PLEX库开发的软件包JavaPlex。

构造复形需要确定四个参数

(1)由边权矩阵构造的点云坐标文件(.txt)；

(2)最大过滤值ε；

(3)最大维数ε_max；

(4)过滤步骤数(Filatration steps,简称Fs)。

下面根据实验情况分别讨论如何确定这几个参数可以达到最佳实验效果。

由边权矩阵构造点云坐标文件，首先将邻接矩阵转换为边权矩阵，矩阵中每一行为“ij ω_ij”，然后用网络嵌入方法(Network Embedding)中的经典算法多维缩放(MultiDimensional Scalling简称MDS)的一个变种算法等度量映射(Isometric FeatureMapping简称ISOMAP)算法，将高维邻接矩阵中两两节点之间的距离在低维上(通常是2维，但是可以是任意维)找到一组新的样本点(即点云)，使降维后两点间的距离与它们在高维上的距离相等。相比MDS算法，ISOMAP算法既能够保留非线性数据的本质几何结构，又能够保持全局的结构信息。

最大过滤值ε_max，按照以上方法构造边权矩阵后，分别取每个阶段各个节点之间距离的最大值作为过滤最大值。如果ε过大会构造出超维复形，现有技术无法可视化和判断其稳定性。

表1编码阶段最大过滤值

Table1 Max filtration value

最大维数和过滤步骤数，最大维数初值设定为3，也就是计算维数0(dim0)，维数1(dim1),维数2(dim2)和维数3(dim3)四个维度下的持续拓扑特征。

过滤步骤数也就是在确定最大过滤值ε后确定过滤步长，根据文献一般Fs＝20，本实验取Fs＝20，100，1000三种情况分别求持续拓扑特征，来确定模型中最优Fs。

首先根据以上选择的过滤算法和最大过滤值ε_max，分别选取过滤步骤数Fs＝20，100，1000三个数，程序在计算机配置为CPU：Interl(R)core(TM)i7-6700；内存32G；操作系统WindowsX64位机上运行。比较三种情况下构造的复形数、运行时间、dim0持续特征数(节点)、dim1持续特征数(线段)、dim2持续特征数(三角形)、dim3持续特征数(四面体)，实验结果如表2所示。

表2三种过滤步骤数下实验结果比较

Table2 Experiment results in three filtration steps

从表中可以得知，三种情况下构造的复形总数不变，第二种情况运行时间比第一种情况多消耗21.14％，但是特征数变化不大，第三种情况运行时间比第二种情况多消耗3.49％，特征数变化较大，因此数据量较大时，权衡时间效率和特征数，建议选择20，也可以选择1000。以下是Fs＝20，100，1000三种情况下的用持续图作为可视化方法展示健康被试持续拓扑特征(图4-图6)。

从三个图中可以清楚的看到dim0区别不大，但是对于dim1和dim2两种情况，虽然Fs＝100和Fs＝1000时特征数变多，但这些特征多分布在对角线附近，也就是持续时间很短的噪声，而间隔为[3.474274,3.709818]始终存在。通过以上分析最终确定过滤步骤选Fs＝20。表中dim3特征数始终是0，所以过滤的最大维数选择2。

精神分裂症WM持续拓扑特征

经以上实验确定过滤参数为(1)最大维数2；(2)各频段健康被试和精分病人最大过滤值ε_max对应取表1中值；(3)过滤步骤Fs＝20时，分别用条形码和持续图可视化求得的持续拓扑特征，以编码阶段为例，展示全频段边权矩阵，α、β1、β2、γ、θ边权矩阵的持续特征。

健康被试和精分病人全频段时持续特征的条形码如图7、8所示。

健康被试全频段边权矩阵持续特征的持续图如图4所示，精分病人的全频段边权矩阵持续特征的持续图如图9所示。

精神分裂症WM持续拓扑特征稳定性分析

实验已经分别提取到健康被试和精分病人在全频段和α、β1、β2、θ、γ五个频段的持续特征，并且用条形码和持续图可视化持续特征，以下采用对持续性图进行比较的方法测度网络相似性。本实验采用两个度量标准：Bottleneck距离和Wasserstein距离来计算。实验数据包采用Python环境下Ripser调用GUDHI包,GUDHI库是一个通用的开源c++库，带有Python接口，用于拓扑数据分析(TDA)。该库提供了最先进的数据结构和算法来构造单纯形和计算持续同源性。Ripser是目前使用Vips过滤计算PH的性能最佳库，调用GUDHI包可以计算Bottleneck距离，GUDHI包3.1.1版本中增加了计算Wasserstein距离的功能。

计算Bottleneck距离的基本思想是基于“瓶颈匹配和相关问题的几何帮助”文献中的思想，需要输入三个参数，两个持续图的间隔数组，一个是计算的精确度参数e，默认情况下e是最小的正双精度数，如果e＝0，结果就是用高昂的算法代价来计算精确距离。本实验e取了两个值，一个是e＝0.01，计算出来是近似值，另一个是e取默认值，也就是双精度数最接近真实值的近似值，我们称之为真实值。计算结果如表3所示。

表3持续图的Bottleneck距离

Table4.3 Bottleneck distance between diagrams

从表中可以得知，近似值和真实值之间的误差都不是很大，除了在dim0的α频段和dim1的γ频段，这两个频段的持续拓扑特征可能存在奇异值。所以最后采用默认值的方式计算持续图间的Bottleneck距离。

计算Wasserstein距离的基本思想是基于“通过最优传输的持久性图均值和聚类的大规模计算”文献中的思想。计算得到的健康被试和精分病人全频段和α、β1、β2、θ、γ各个频段下两个维度的Wasserstein距离如表4所示。

表4持续图的Wasserstein距离

Table4 Wasserstein distance between diagrams

本发明提出全尺度复杂脑网络模型并应用于精神分裂症病人WM数据分析中，经实验分析确定模型中的重要的参数和算法，在边权矩阵和点云文件互换时，运用ISOMAP算法降维后，矩阵由60维基本上降到了34～38维之间，为后续高效率处理数据提供了良好的基础。

几个重要性参数通过全方位的实验效果比较，最后确定为最大维数维1，因为健康被试β2、θ和γ三个频段在dim2，精分病人β2、和γ两个频段在dim2都没有找到持续性拓扑特征，所以在稳定性分析时，实验分析了dim0和dim1两个维度。最大过滤步骤数确定为Fs＝20，实验分析在这个过滤步长下，即丢失不了重要的特征，又可以剔除噪声，同时还可以提高时间效率。从确定最优参数后构造的复形数比较可以在全频段和γ频段差异性很大，其它频段差异较小。

在稳定性分析中，持续图间的Bottleneck距离结果(表3)显示，精确度参数e采用默认值实验效果更容易接近真实值，但是当持续图没有奇异值是真实值和差异值差别不大。实验结果(表3)还进一步表明，在α和θ频段Bottleneck距离较小，也就是健康被试和精分病人的持续图大的变化较小，可以选取这两个阶段的持续拓扑特征作为模型的输出。

本发明提出基于持续同调性的全尺度脑网络分析模型，并对模型中的每一数据处理步骤涉及到的算法和参数做了分析，研究了全尺度脑网络构造中的节点、边权矩阵构造、过滤阈值的选择等关键性问题。将模型应用于精神分裂症WM持续拓扑特征分析中，在获取持续拓扑特征的基础上，实验结果还表明本文提出的模型具有稳定性和对抗噪性，可以为精神类疾病医学影像分析提供稳定生物参考指标。本发明不但可以解决构建复杂脑网络时黄金阈值选择难的问题，而且根据尺度的变化可动态观测每一个有效尺度下网络拓扑结构的变化。

以上仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于持续同调的脑电信号持续特征提取方法，其特征在于，包括以下步骤：

S1、获取大脑成像数据并进行预处理，得到EEG时间序列信号，确定大脑的网络节点；

S2、将网络节点映射到有限度量空间构造邻接矩阵，度量空间采用皮尔逊相关性度量，将各网络节点间皮尔逊相关系数的倒数作为网络节点间连接的权重；

S3、通过邻接矩阵构造边权矩阵，选择过滤阈值和过滤算法对边权矩阵进行过滤，构造脑网络复形；

S4、可视化表示过滤过程中脑网络的拓扑特征变化，并获取持续拓扑特征；

S5、分析持续拓扑特征的稳定性，并输出脑网络的持续拓扑特征。

2.根据权利要求1所述的提取方法，其特征在于：在步骤S4中，通过条形码或者持续性图可视化表示过滤过程中脑网络的拓扑特征变化。

3.根据权利要求2所述的提取方法，其特征在于：在过滤过程中，计算p维贝蒂数间隔[ε_birth，ε_death]；其中，ε_birth为单纯复形中p维孔开始的时间，ε_death为单纯复形中p维孔消失的时间。

4.根据权利要求3所述的提取方法，其特征在于：当通过持续性图可视化表示过滤过程中脑网络的拓扑特征变化时，持续性图中横坐标表示持续特征出现的时间ε_birth，纵坐标表示持续特征消失的的时间ε_death，将过滤过程中求得的间隔集合[ε_birth，ε_death]作为持续性图中点的坐标，将集合中所有的间隔对表示在坐标中。

5.根据权利要求3所述的提取方法，其特征在于：当通过条形码可视化表示过滤过程中脑网络的拓扑特征变化时，条形码中横坐标表示过滤阈值ε，条形码的长度表示[ε_birth，ε_death]的长度。

6.根据权利要求1所述的提取方法，其特征在于：在步骤S5中，采用Bottleneck距离和Wasserstein距离来分析持续拓扑特征的稳定性。

7.根据权利要求6所述的提取方法，其特征在于：

两个图X和Y之间的p阶Wassertein距离定义为：

其中，p∈「1，∞)，

X→Y是从X到Y的映射；

当p＝∞时，距离d是二维空间的度量，以上公式表示为：

其中，Bottleneck距离为W_∞[d_∞]。

Images