基于多通道图谱结构特征挖掘的癫痫发作预测方法

本发明涉及医疗图像处理，尤其涉及一种基于多通道图谱结构特征挖掘的癫痫发作预测方法。

背景技术：

1、癫痫是一种由大脑的大量神经元异常同步放电导致的中枢神经系统疾病，病因包括但不限于遗传、药物、脑损伤、代谢紊乱以及心理因素，发作时的常见症状为意识障碍、感官异常以及不自主的抽搐痉挛。癫痫的发作不定时、迅速且容易反复，患者可能会因为突然的抽搐和意识障碍导致肢体损伤甚至死亡，此外癫痫的污名也会带来心理负担。大部分患者可以通过药物或者病灶切除手术获得良好的疗效，但是仍然有30%左右患者为难治性癫痫。

2、癫痫发作预测对于尚未获得良好疗效的癫痫患者具有重要意义。通过癫痫发作预测，可以提醒患者提前采取自我保护措施，通知家属或者医护人员提供救助，还可以联动经颅磁/电刺激、脑起搏器等装置抑制甚至避免即将到来的癫痫发作，最大程度降低癫痫导致的伤害。

3、癫痫的发作可以分为4个时期：发作前期、发作期、发作后期和发作间歇期。发作前期也称为癫痫先兆期，一般指癫痫发作前的数十分钟，发作间歇期是指未受癫痫影响的正常时期。癫痫发作预测的本质是区分发作前期和发作间歇期。

4、相比于颅内脑电、核磁共振成像、脑磁图、正电子发射断层扫描等神经成像方式，脑电图（electroencephalogram，eeg）因为其廉价、便携、安全、时间分辨率高等优点，在癫痫发作预测领域使用最为广泛。基于脑电图的癫痫发作预测方法可以分为两大类：传统机器学习方法以及基于神经网络的方法。

5、传统的机器学习可以根据其设计的特征细分为两大类：线性特征方法和非线性特征方法。线性特征方法主要使用脑电的时域、频域、时频域、空间域特征，非线性特征方法主要使用李雅普诺夫指数或者熵度量。值得一提的是，脑电作为一种非线性的多源混合信号，单一的一种特征往往没有包含分类所需的足够的判别性信息，因此这些方法普遍使用了多种特征进行组合，这导致系统的冗余和复杂。

6、神经网络方法使用卷积神经网络、长短期记忆网络、transformer（转换器网络）以及图神经网络等结构，还可以结合自编码器、残差连接等模块，实现端到端的自动特征提取和分类。凭借强大拟合能力，这类方法往往能够在已有数据上获得优秀的性能，但是这类方法缺乏可解释性，并且过拟合的普遍存在以及脑电的个体间差异性与时变性导致其泛化性能不佳；增加数据量是提高泛化性能的有效方式，但是在癫痫发作预测等医疗相关领域，数据往往是稀少而宝贵的。

7、在癫痫预测或者癫痫检测时，对患者进行身份识别这一需求在现有的工作中没有得到关注。使用同一个系统框架实现癫痫预测以及身份识别，可以为患者提供个性化的服务，例如：基于身份识别结果，加载针对该患者优化的系统运行参数；为阻断癫痫发作，实施专为该患者设定的经颅磁/电刺激、脑起搏器治疗参数；通知该患者对应的紧急联系人。

8、此外，如果未能阻止发作，处于癫痫状态的患者很可能无法与他人有效沟通。身份识别可以帮助医护人员迅速获取患者的病史信息（例如过敏史、癫痫细分类型、现行治疗方案等），这对于采取正确的急救措施至关重要。

技术实现思路

1、本发明针对现有技术存在的上述问题，提供了一种基于多通道图谱结构特征的癫痫发作预测法，其以脑电信号的时频图等高线的轮廓特征为分类特征实现自动预测癫痫发作，进而还可以基于该等高线的轮廓特征对预测出癫痫将会发作（癫痫发作前期）的对象进行患者身份的自动识别。

2、本发明提供的技术方案为：

3、基于脑电多通道图谱结构特征挖掘的癫痫自动预测方法，包括以下步骤：

4、步骤1，对脑电采集设备输入的原始多通道脑电信号进行信号预处理；

5、步骤2，对信号预处理后的每个通道的脑电信号进行连续小波变换，得到变换后的每个通道的脑电数据，再使用非负tucker（塔克）分解对所有通道的脑电数据构成的三阶张量进行降维，对降维后的时频图进行归一化处理得到降维的多通道时频域图谱；

6、步骤3，以非均匀采样间隔抽取降维的多通道时频域图谱的每个通道的一组等高线，提取每个通道的所有等高线的轮廓特征得到当前通道的时频图特征向量；基于所有通道的时频图特征向量得到最终的脑电信号特征；

7、步骤4，将最终的脑电信号特征输入训练好的袋装树分类器，进行是否处于癫痫发作前期的自动识别。

8、进一步的，步骤4还包括：若袋装树分类器的识别结果为是（即预测出癫痫将会发作），将对应的脑电信号特征输入线性判别分析（linear discriminant analysis，简称为lda）分类器，进行患者身份的自动识别。

9、本发明中，对癫痫预测，将最终的脑电信号特征（记为特征向量）输入已经训练好的预置袋装树分类器，实现自动预测癫痫发作。袋装树是一种集成学习方法，通过构建多个决策树并对它们的预测结果进行组合来提高模型的准确性和稳定性，不易过拟合，具有一定的鲁棒性和抗噪声能力，适合解决癫痫预测面临的低信噪比、低数据量、高时变性问题；而对于患者身份识别，袋装树在预测出癫痫将会发作时，会将同一个特征向量  输入预置的线性判别分析分类器，借助特征共用实现患者身份识别，避免对身份识另行计算特征带来的处理成本，有利于依赖身份识别的下游任务快速开展。由于个体间脑电具有较大的差异性，身份识别不需要袋装树这类繁重的集成学习机制，因此本发明采用在适合多分类任务的简洁高效的线性判别分析。

10、进一步的，步骤1中的信号预处理包括：降采样和滤波，以有利于降低数据大小，减少信号噪声；

11、进一步的，降采样时，将原始多通道脑电信号降采样到指定的采样频率，优选200hz，允许一定范围的偏差。由于是在头皮采集的脑电信号，大脑内部原有的高频成分在头骨和皮肤组织中会遭遇较高的阻抗，采集到的原始多通道脑电信号的有效频率范围大致在0.5hz到60hz区间。采样频率只需要保证有效频率范围的信号不会混叠即可，过高的采样频率反而会增加存储和信号处理的代价。因此本发明选择将原始脑电信号降采样到200hz。这样做还有一个好处：在将本发明部署在不同的脑电采集硬件时，信号处理模块的入口信号获取恒定采样频率的数据，有利于系统的统一设计和稳定运行。

12、进一步的，滤波时采用带通滤波，优选0.5hz~40hz的带通滤波，滤波器优选4阶的巴特沃斯滤波器，波可以抑制有效频率范围以外的噪声分量，因而采用0.5hz~40hz的带通滤波，巴特沃斯滤波器是线性相位滤波器，并且它在通频带内具有平滑的响应，可以较好地避免脑电信号失真。

13、进一步的，连续小波变换时，采用bump小波（邦普小波）作为小波基，以获得良好的震荡尖峰捕获能力；此外，对于小波变换获得的时频图（所有通道的脑电数据构成），截取频率位于0.5 hz ~ 40hz的主要部分进行非负tucker分解，以进一步降低计算量。癫痫前期脑电信号在时域可能不会有显著的变化，在频域进行分析是一个更好的视角。但由于脑电信号的时变特性，其频率成分会随着时间动态变化，本发明采用连续小波变换对每一个通道的脑电信号进行处理，并选择广泛使用的bump小波作为小波基，来获得良好的震荡尖峰捕获能力。

14、目前的脑电采集设备除非是特殊场景使用，都具有多个通道来记录不同脑区的电信号。这些通道记录到的脑电信号不是独立的：由于容积传导效应，当前通道的电极采集到的是颅内不同脑区多个独立源产生信号的叠加，这也就意味着采集到的多通道脑电信号存在很多的相关性和冗余，因而本发明采用非负tucker分解对小波变换获得的时频图的截取部分进行降维处理。在本发明中，将每个通道的二维时频矩阵进行堆叠，很自然地就构成了一个三阶张量。张量是矩阵的推广，是多维数据的天然表示。tucker分解是一种常用的张量分解方法，可以视为奇异值分解的高维推广，具有很好的数据压缩和降维能力，可以提取高维数据中的关键特征。原始的tucker是正交的，非负约束的添加使得tucker分解得到的核张量和模式矩阵都是非负的，更符合物理实际，具有良好的可解释性。因此本发明采用非负tucker分解，使用分解结果中的核张量作为原始数据的降维表示，在保留关键信息的情况下，同时减少数据的时频点数和通道数。非负tucker分解的优化目标为：

15、，

16、，

17、其中， 为原始脑电张量（由小波变换获得的时频图的截取部分所构成，即分解前的脑电张量），表示实数域， 表示分解前的时间点数，表示分解前的频点数，表示分解前的通道数量，为分解的核张量， 为不同维度的模式矩阵，n为维度标识，其取值为1,2,3，为单位矩阵，表示f范数，和为模-n乘，用于张量和矩阵之间的乘法，表示核张量的第一个模态和模式矩阵相乘，表示张量“”的第二个模态和模式矩阵相乘。该优化问题可以通过非负交替最小二乘、高阶乘法及其他众多算法求解。本发明中采用非负交替最小二乘法，因为它的内存消耗和迭代次数比较平衡，适合用在类似本发明的数据规模不大的场景。

18、非负tucker分解的秩是影响性能的关键超参数，本发明中，秩取原始张量对应方向上维数的一半。具体地，降维以后的压缩数据为，数据大小将会降为原来的八分之一。其中的即为tucker秩。

19、进一步的，步骤3中，等高线的轮廓特征包括等高线比重特征、时域分布特征和频域分布特征。

20、进一步的，步骤3具体包括：

21、步骤301，采用采样间隔依次线性递减的方式抽取每个通道的时频域图谱的一组等高线，分别收集每条等高线上的所有坐标点构成等值的二维坐标点集，其中，下标为每个通道的等高线标识符；

22、步骤302，将二维坐标点集分别向时间轴和频率轴投影，获得时间轴的一维点集和频率轴的一维点集；

23、由二维坐标点集的元素个数计算等高线比重特征，分别计算一维点集和的均值和方差作为时域分布特征和频域分布特征；

24、拼接每个通道的所有等高线的等高线比重特征、时域分布特征和频域分布特征得到当前通道的时频图特征向量；基于所有通道的时频图特征向量得到最终的脑电信号特征；

25、此外，步骤2中，归一化处理具体为：对执行非负tucker分解后得到的每个通道的时频图（即降维后每个通道的时频图），将功率归一化到 (0, 1] 区间，得到归一化时频图g。这样做可以保持对不同通道的时频图进行等高线采样时，可以使用相同的采样点。

26、进一步的，步骤 301具体为：

27、在 (0, 1] 区间，从0开始，以预置的初始间隔（优选0.05）开始采样，每采完一个点，采样间隔按照指定值（优选值为0.001）进行衰减，对采集到的若干个点构成采样点集合s；

28、对于采样点集合s中的任意一点，定义其对应的高度为x，点表征当前通道的归一化时频图g的值为x的一条等高线；将归一化时频图g中所有值在  范围的点都归属于点所对应的等高线，得到当前等高线的二维坐标点集 ，每一个二维坐标点集 均对应着一个等高线轮廓，其中，表示预置的偏差，优选取值为0.02；这样做的原因是：小波变化计算得到的功率值和当前时刻这个频点的真实功率值是存在偏差的，对等高线的取值进行微小的扩大可以避免遗漏。

29、进一步的，步骤 302中，等高线比重特征、时域分布特征和频域分布特征分别为：

30、等高线 的等高线比重特征为：，代表了当前等高线对应的功率值在所有时频点中的占比；其中，所有等高线的总点数，表示每个通道的时频域图谱的等高线条数， 表示等高线 的二维坐标点集的元素个数；

31、时域分布特征为一维点集的均值和方差；

32、频域分布特征为一维点集的均值和方差。

33、本发明提供的技术方案至少带来如下有益效果：

34、（1）本发明可以通过从脑电采集硬件获取的脑电信号，使用非负tucker分解降维，并采用时频图等高线的轮廓特征作为共享特征，进而实现癫痫的自动预测，以及患者身份的自动识别。相比于现有的技术，本发明存储开销小，可以协助癫痫患者有准备地应对癫痫发作；

35、（2）本发明使用非负tucker对数据降维，降低了后续处理流程的计算量。相比于传统的主成分分析等降维方法，可以更好地探索和保留原始脑电数据在时域、频域、空间域的复杂关联；

36、（3）本发明使用功率等高线的轮廓作为时频域的多通道图谱结构特征。脑电信号在时频域具有丰富的判别性信息，其中最主要的是在一些特定时间点上出现高功率的频率成分，时频图的等高线轮廓特征可以很好地刻画这种时间、频率和能量的耦合关系。由稀疏变稠密的非均匀等高线采样，在减少等高线数量的同时，确保了对高功率区间的细致刻画；

37、（4）本发明直接使用癫痫预测时已获取的特征来进行身份识别，可以简化系统模块，加快身份识别响应速度。