CN112826512B - 一种癫痫棘波自动检测与峰值定位方法 - Google Patents

Abstract

本发明公开了一种癫痫棘波自动检测与峰值定位方法，其步骤包括：1）将每一脑磁图数据拆分成多个大小相同的数据分段；2）构建深度学习网络模型，包括单通道编码模块，深度特征传递模块以及全局解码模块；其中，单通道编码模块用于对数据片段中的单通道数据进行深度特征提取，得到多种不同尺寸的深度特征并发送给深度特征传递模块进行数据降维；全局编码模块用于从尺寸最小的深度特征开始，对深度特征进行卷积与上采样计算，并将得到的深度特征与大小对应一致的深度特征进行拼接，直至处理完尺寸最大的深度特征；3）利用各数据分段训练深度学习网络模型；4）将待处理的脑磁图数据输入训练后的深度学习网络模型进行棘波检测和峰值定位。

Description

一种癫痫棘波自动检测与峰值定位方法

技术领域

本发明属于生物特征识别领域中的脑磁图信号识别领域，具体为一种基于深度学习的癫痫棘波自动检测及峰值定位方法，不仅能够判断脑磁图数据中是否存在棘波，还能确定棘波峰值点的时刻，帮助临床医生减少工作量，更加准确的定位致痫灶。

背景技术

癫痫是一种常见的神经性脑疾病，其特征主要为大脑神经元过度放电异常导致的中枢神经***功能失常。对于癫痫的诊断与治疗，其关键在于确定异常放电的大脑区域，即致痫灶。通过对致痫灶进行针对性的药物抑制或外科手术切除，是治疗癫痫疾病的有效方式。

脑磁图（Magnetoencephalography，MEG）是一种无创性检测大脑神经元放电活动的技术，其通过超过300个传感器同时记录脑组织内的神经活动。相比脑电图，脑磁图具有更好的信噪比，并且能够抵抗复杂脑组织导致的信号失真。癫痫因其异常放电产生的棘波在脑磁图具有十分明显的特异性，棘波作为一种生物标志物被广泛接受。它可以用于判断癫痫的异常放电时间，识别癫痫的发作情况，这使得癫痫患者的脑磁图棘波分析具有十分重要的临床意义。

然而，癫痫异常放电具有发作性，短暂性、反复性的特点，这些特点导致一位癫痫患者的脑磁图检查一般需要持续1.5小时。并且在临床上，棘波检测一般是由临床医生人工进行主观标注，由于脑磁图具有高时间分辨率和高采集通道数的特点，临床医生的主观检查费时费力，并且高度依赖与临床医生的经验。一位经验丰富的临床医生，对一次脑磁图数据进行完全检查，需要观察300多个通道，约6000s，1000Hz采样率的脑磁图信号，耗时在2小时左右，对于资历较浅的医生，耗时可能更长。在脑磁图的长采集时间和高诊断专业要求的共同影响下，催生出了许多自动化的、高准确度的棘波检测方法。

目前，已经有多种自动化的棘波检测算法，通过对脑磁图数据进行时域分析、频域分析、时频分析，并进行相关特征提取，从而实现固定时间片段的有无棘波检测。也有部分研究使用深度学习网络，直接从脑磁图数据入手，进行棘波检测。但是，目前已有的算法还只能满足于检测固定时间片段内是否包含棘波，并不能准确定位棘波能量最高的时刻点，即棘波峰值点。对于致痫灶的定位往往需要使用溯源算法，使用长时间段的棘波数据段并不能得到一个准确的致痫灶的定位，若使用了棘波峰值时刻点，针对峰值周围一小段数据进行溯源定位，能够得到定位更加准确的致痫灶。因此，对于相比于时间段内的棘波检测，棘波的峰值定位具有更大的实现意义。

发明内容

为了改善目前已有的算法在峰值检测上的不足，本发明提供了一种基于深度学习的癫痫棘波检测与峰值定位方法，该方法能有效的从脑磁图的原始数据中识别出癫痫棘波放电事件并且准确定位棘波峰值的时刻点。

本发明提出的深度学习网络通过对脑磁图每个通道的信号分别进行编码，获取信号的深度特征，再对深度特征进行特征融合，共同解码。在解码步骤中，将相同尺寸大小的编码数据与解码数据进行连接，提高网络对特征的复用。最终网络将输出脑磁图信号在各个时间点上的峰值点概率，通过后续的阈值算法对棘波进行判断以及准确定位棘波峰值点的位置。

本发明的技术方案为。

一种基于深度学习的脑磁图棘波检测与棘波峰值点定位方法，其步骤包括：

步骤一：对脑磁图进行数据预处理：去除脑磁图数据中的大量伪迹信号，并将数据拆分成大小相同的数据分段。首先使用带通滤波器与工频陷波滤波器进行滤波，再使用ICA对眼跳信号与心跳信号进行伪迹去除，随后对数据进行标准化，使其服从标准正态分布。最后对脑磁图数据进行数据拆分与分割，基于大脑的不同脑区，将脑磁图数据拆分成8个脑区，每个脑区上39个脑磁图通道，同时在时间上，对脑磁图数据以296毫秒的时间窗，148毫秒的重叠进行数据分割，得到脑磁图数据片段。在步骤1）中，最终将得到多个大小为[39*296]的脑磁图数据片段。

步骤二：构建深度学***均，输出给全局解码模块。全局编码模块则首先对尺寸大小最小的深度特征进行卷积与上采样计算，并将得到的深度特征与特征传递模块中尺寸大小一致的深度特征进行特征拼接，重复上述计算，最终输出数据大小与网络输入数据在时间维度上大小相同的一维矩阵。因此网络输出矩阵中每一个值都可以与原始信号（即输入的数据分段）上的每一个时刻点进行一一对应，矩阵中的每一个值代表该值所对应时刻点的网络预测为棘波峰值点的概率。

步骤三：训练深度学***滑重构；二是数据增强，对于每一个脑磁图数据片段，对其时间窗进行随机左右平移，平移时间不超过50毫秒（原数据片段为固定的时间点，即每148ms进行一次长度为296ms的数据段截取），平移9次则数据增强倍数为9倍。三是损失函数的选择，本发明方法使用如下损失函数：

。

其中

为网络输出矩阵中索引值为

的矩阵值，即每一个对应的时刻点预测为棘波峰值的概率，

为重构后的标签矩阵中索引值为

的矩阵值，即每一个对应的时刻点重构后的标签值。

步骤四：对一待处理的脑磁图数据进行棘波峰值检测与峰值定位。将预处理后的脑磁图数据片段，输入训练好的深度学习网络中，得到输入的脑磁图数据片段对应的棘波峰值预测矩阵，再将其输入至检测定位模块，计算输出数据的峰度，检测定位模块针对每一个脑磁图片段的峰度使用阈值法，判断该脑磁图片段是否包含棘波。最后，在包含有棘波的片段中选取概率最高的时刻点作为预测棘波峰值的具体时刻点(

。

在步骤四中，本发明提出了针对峰值定位的评价指标，峰值偏差b_p，其计算公式为：

。

即根据样本i的预测的峰值点

与样本i的实际的峰值点

的误差绝对值，计算得到样本i的峰值偏差

，单位为毫秒。

本发明所提出的癫痫棘波检测与峰值定位方法，能够有效的脑磁图的各个数据片段中识别出棘波，并且准确定位棘波峰值点在片段中的时刻点。结果表明，本发明提出的棘波检测方法具有与现有方法相当的性能。同时，首次提出了峰值检测的评估指标，能够准确地预测棘波峰值点的时刻，在测试集中，峰值预测的峰值偏差平均值为14毫秒，标准差为17毫秒。

和现有的技术相比，本发明的有益效果是。

1. 本发明直接从数据本身着手，能够准确的预测棘波以及棘波峰值点的位置，为后续癫痫致痫灶的溯源定位提供了精确的指标。

2. 本发明使用了目标检测的思路，将棘波检测与峰值定位融合至同一个网络任务中，提高了模型的实用性。

3. 本发明使用的网络结合了单通道与全局通道的信息，并且将不同尺寸大小的特征从解码器传递给编码器，提高了特征的复用，使网络能够更加准确地预测棘波峰值时刻点。

附图说明

图1为本发明方法的流程图。

图2为本发明方法中预处理模块。

图3为本发明方法中深度学习网络框架。

图4为本发明方法中深度学习网络中的单通道卷积模块。

图5为本发明方法中深度学习网络中的卷积模块。

图6为本发明方法中深度学习网络中的深度特征传递模块。

图7为本发明方法中深度学习网络中的全局解码模块。

图8为本发明方法中标签重构后的标签值。

具体实施方式

下面将结合附图对本发明方法的实施例做进一步说明。

如图1所示，为本发明实施例的棘波自动检测与峰值定位流程图，主要包含以下步骤。

步骤一：数据预处理。

此步骤对输入的脑磁图原始数据进行预处理，对信号进行滤波，伪迹去除，标准化和分段操作，使数据能够更好的放入深度学习网络中进行训练。

如图2所示，为数据预处理模块，将输入的脑磁图数据进行以下操作的预处理，主要包含以下步骤。

1) 使用1-100Hz带通滤波器与50Hz陷波滤波器，除去各个通道上的低频、高频和工频噪声。

2) 使用独立成分分析（ICA），除去信号中与眼跳信号、心跳信号的伪迹。

3) 对各个通道上的信号进行数据标准化，使数据服从标准正态分布。

4)对脑磁图数据进行拆分与分割，根据脑磁图通道的定位信息，将306个脑磁图通道拆分成8个脑区，每个脑区39个通道。对于不足39通道的脑区，使用全零通道数据进行补全。

5) 对脑磁图数据进行训练集和测试集的拆分，训练集 : 测试集比例为8 : 2。

步骤二：构建深度学习网络。

此步骤为训练的深度学习网络，通过输入训练集脑磁图数据分段进行训练，预测脑磁图数据分段中每个时刻点为棘波峰值的概率。

如图3所示，为本发明的深度学习网络框架，主要包含以下算法模块：单通道编码模块301，深度特征传递模块302和全局解码模块303。以下将分别介绍各个模块的细节。

1. 单通道编码模块301。

如图4所示，为单通道编码模块的基本结构。主要包括脑磁图数据输入层，4个卷积模块，3个最大池化层和4个输出层。

脑磁图数据输入层401：输入脑磁图片段数据中的单通道数据，数据大小为[1*296]。

如图5所示，为卷积模块（包括第一卷积模块402，第二卷积模块405，第三卷积模块408，第四卷积模块411）的基本结构，数据通过卷积模块，进行一维卷积，提取深度特征。卷积模块包括了第一卷积层501，第一激活函数502，第一批标准化层503，第二卷积层504，第二激活函数505，第二批标准化层506。

第一和第二卷积层使用了一维卷积核，对输入数据进行卷积计算，填充方式为相同零填充，即卷积计算输出数据大小与输入数据大小一致，数据两端使用零进行自动填充。每个卷积模块中定义的卷积核大小和卷积核数量根据卷积模块的在网络中的深度而定。第一卷积模块402中两个卷积层的卷积核大小都为[1*20]，卷积核数量都为16。第二卷积模块405中两个卷积层的卷积核大小都为[1*10]，卷积核都数量为32。第三卷积模块408中两个卷积层的卷积核大小都为[1*5]，卷积核数量都为64。第四卷积模块411中两个卷积层的卷积核大小都为[1*3]，卷积核数量都为128。

第一和第二激活函数为线性整流函数（Rectified Linear Unit, ReLU），用于增加网络的非线性拟合能力。第一和第二批标准化层使每一次计算的输入都服从标准正态分布，提高网络的收敛速度，防止产生梯度***的现象。

最大池化层（包括第一最大池化层404，第二最大池化层407，第三最大池化层410）：为了减少网络的参数以及防止网络产生过拟合，使用了池化大小为2，步长为2的最大池化层。使输出的数据大小只有输入的数据大小一半。

输出层（包括第一输出层403，第二输出层406，第三输出层409，第四输出层412）：单通道编码模块一共有四个尺寸大小不同的数据输出，第一输出层403输出矩阵大小由第一卷积模块402决定，大小为[16*296]。第二输出层406输出矩阵大小由第二卷积模块405决定，大小为[32*148]。第三输出层409输出矩阵大小由第三卷积模块408决定，大小为[64*74]。第四输出层412输出矩阵大小由第四卷积模块411决定，大小为[128*37]。

需要注意的是，单通道编码模块会对脑磁图数据分段中的39个通道都做上述的计算，最终将得到39个通道，每个通道有4种不同尺寸大小的深度特征输出。

2. 深度特征传递模块302。

深度特征传递模块将单通道编码模块的输出层在脑磁图通道维度上进行了全局池化平均后，传递给全局解码模块303。

如图6所示，为深度特征传递模块302的基本结构，共包括了4个全局池化层（第一全局池化层601，第二全局池化层603，第三全局池化层605，第四全局池化层607）和4个全局输出层（第一全局输出层602，第二全局输出层604，第三全局输出层606，第四全局输出层608）。

第一全局池化层601：输入39个通道在单通道编码模块中第一输出层403的深度特征，输入数据大小为[39*16*296]，并在脑磁图通道维度上进行全局池化计算，第一全局输出层602输出数据大小为[16*296]。

第二全局池化层603：输入39个通道在单通道编码模块中第二输出层406的深度特征，输入数据大小为[39*32*148]，并在脑磁图通道维度上进行全局池化计算，第二全局输出层604输出数据大小为[32*148]。

第三全局池化层605：输入39个通道在单通道编码模块中第三输出层409的深度特征，输入数据大小为[39*64*74]，并在脑磁图通道维度上进行全局池化计算，第三全局输出层606输出数据大小为[64*74]。

全局池化层607：输入39个通道在单通道编码模块中第四输出层412的深度特征，输入数据大小为[39*128*37]，并在脑磁图通道维度上进行全局池化计算，第四全局输出层608输出数据大小为[128*37]。

3. 全局解码模块303。

全局解码模块303将深度特征传递模块302的4个输出作为输入，通过卷积计算与上采样，最终得到大小为[1*296]的数据，每一个数据点代表原始脑磁图数据段上对应的时刻点为棘波峰值的概率。

如图7所示，为全局解码模块303的基本结构，包括了4个输入层，3个卷积模块，3个上采样层和1个输出层。

全局输入层701：输入深度特征传递模块302中的第四全局输出层608输出的深度特征。

上采样层（第一上采样层702，第二上采样层705，第三上采样层708）：为了使深度特征传递模块中各个尺寸大小的数据能够与全局解码模块的数据进行拼接，对全局输入层701和卷积模块704、707计算后的数据进行了2倍上采样。第一上采样层702对第四全局输出层608输出的深度特征数据进行上采样后，得到的数据大小为[128*74]，第二上采样层705输出的数据大小为[64*148]，第三上采样层708输出的数据大小为[32*296]。

输入层（第一输入层703，第二输入层706，第三输入层709）：将上采样层和深度特征传递模块的输出的数据进行拼接，实现特征的复用。第一输入层703拼接了第三全局输出层606与第一上采样层702输出的数据，最终数据大小为[192*74]。第二输入层706拼接了第二全局输出层604与第二上采样层705的数据，最终数据大小为[96*148]。第三输入层709拼接了第一全局输出层602与第三上采样层708的数据，最终数据大小为[48*296]。

卷积模块（即第五卷积模块704，第六卷积模块707，第七卷积模块710）：此处卷积模块的结构与单通道编码模块301中的卷积模块结构相近。第五卷积模块704中两个卷积层的卷积核大小都为[1*5]，卷积核数量都为64。第六卷积模块707中两个卷积层的卷积核大小都为[1*10]，卷积核数量都为32。需要注意的是，第七卷积模块704中，两个卷积层的卷积核大小都为[1*20]，但其中第一卷积层的卷积核数量为16，第二卷积层的卷积核数量为1。

输出层711：将第七卷积模块710的输出结果经过Sigmoid激活函数进行计算后，作为整个网络的最后输出结果，大小为[1*296]。

步骤三：训练深度学习网络。

在训练深度学***滑重构，主要步骤为。

1. 用ipeak代表峰值所在时刻点在标签矩阵中的索引，

为标签矩阵中的索引值。

2. 计算标准差为0.01，均值为

的正态分布函数，将函数最大值取为1，其函数表达式为：

。

3. 计算标准差为0.10，均值为

的正态分布函数，将函数最大值取为0.5，其函数表达式为：

。

4. 取

为：

。

作为新的标签矩阵上索引值为

的标签值。

图8展示了原始标签与重构后的标签。

对于没有棘波的脑磁图数据片段，标签值为全0，即大小为[1*296]的零矩阵。

对于每个脑磁图数据片段，通过时间窗随机平移实现了9倍的数据增强，其中时间窗最大平移时间为50ms。

通过深度学***时，停止训练，防止模型过拟合。

在损失函数的选取上，本发明使用了如下损失函数：

。

目的是让模型输出与重构后标签更加接近的预测值。对于存在棘波的脑磁图数据片段，预测结果将呈现出单峰，而对于不存在棘波的脑磁图数据片段，预测结果全都趋近于0。这样能够为后续的棘波预测与峰值定位带来更好的评判依据。

步骤四：棘波预测与峰值定位

此步骤对深度学习网络输出结果进行进一步处理，输出该脑磁图数据片段中是否包含棘波以及棘波峰值点的时刻点。

网络输出为[1*296]的矩阵，其中每一个点代表对应时刻点为棘波的概率，由于在训练阶段对数据标签进行了重构，因此理想情况下，存在棘波的脑磁图数据片段经过网络计算后体现出单峰的趋势，而不存在棘波的脑磁图数据片段经过计算后则体现出均匀分布的趋势。因此本发明使用峰度作为评判指标，判断当前脑磁图数据片段是否存在棘波。对于均值为

，方差为

的任意样本X，峰度的定义为样本的四阶中心矩，其计算公式为：

。

定义样本T，令网络输出矩阵为样本T的概率密度函数，T的取值

为网络输出矩阵的索引值，样本T取

的概率值

为网络输出矩阵中索引值为

的矩阵值，因此可以计算样本T的均值：

。

因此样本T的峰度计算公式为：

。

在本发明中，使用2作为阈值，当峰度小于2时，判断为没有棘波，峰度大于2时，判断为存在棘波。对于存在棘波的脑磁图数据片段，选择概率最大的预测点作为棘波峰值点。通过计算预测的峰值点与实际峰值点的偏差值，峰值点的偏差值b_p，评判模型的预测性能。

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

Claims (7)

1.一种癫痫棘波自动检测与峰值定位方法，其步骤包括：

1-1.对脑磁图进行数据预处理：去除脑磁图数据中的伪迹信号，然后将每一脑磁图数据拆分成多个大小相同的数据分段；

1-2.构建深度学***均后，输出给所述全局解码模块；所述全局解码模块用于从尺寸大小最小的深度特征开始，对深度特征进行卷积与上采样计算，并将得到的深度特征与所述深度特征传递模块中尺寸大小对应一致的深度特征进行特征拼接，直至处理完尺寸大小最大的深度特征，最终输出数据大小与网络输入数据在时间维度上大小相同的一维矩阵；其中所述一维矩阵中每一个值与输入数据分段上每一个时刻点一一对应，每一个值的大小代表该值所对应时刻点的网络预测为棘波峰值点的概率；

1-3.利用各所述数据分段训练所述深度学习网络模型；

1-4.对于一待处理的脑磁图数据，将其拆分为多个数据片段后输入训练后的所述深度学习网络模型，得到每一数据片段的预测矩阵；检测定位模块根据该预测矩阵计算各数据片段对应的峰度，将峰度大于设定阈值的数据片段判定为包含棘波，然后在包含有棘波的数据片段中选取概率最高的时刻点作为预测棘波峰值的具体时刻点。

2.如权利要求1所述的方法，其特征在于，对训练所述深度学习网络模型的数据分段的标签进行重构，其方法为：

2-1.根据数据分段中标注的峰值点位置，生成二值化的标签矩阵；用ipeak代表峰值所在时刻点在标签矩阵中的索引，x为标签矩阵中的索引值；

2-2.计算标准差为0.01，均值为ipeak/296的正态分布函数f₁(x)，f₁(x)的最大值取为1；

2-3.计算标准差为0.10，均值为peak/296的正态分布函数f₂(x)，f₂(x)的最大值取为0.5；

2-4.取f(x)＝max(f₁(x),f₂(x))作为标签矩阵上的每一个点的标签值。

3.如权利要求1或2所述的方法，其特征在于，训练所述深度学习网络模型所使用的损失函数为

其中p_t为网络输出矩阵中索引值为t的矩阵值，即时刻点t预测为棘波峰值的概率，l_t为重构后的标签矩阵中索引值为t的矩阵值，即时刻点t重构后的标签值。

4.如权利要求1或2所述的方法，其特征在于，对训练所述深度学***移，平移时间不超过若干毫秒，获得多个数据分段。

5.如权利要求4所述的方法，其特征在于，所述时间窗为296毫秒，平移时间不超过50毫秒。

6.如权利要求1所述的方法，其特征在于，计算数据片段网络输出对应的峰度的方法为：定义样本T，令网络输出矩阵为样本T的概率密度函数，其中网络输出矩阵的索引值t为样本T的取值，网络输出矩阵中索引值为t的矩阵值p(t)即为样本T取t的概率值，其峰度为

7.如权利要求1所述的方法，其特征在于，通过计算预测峰值点与实际峰值点的偏差值b_p，评判所述深度学习网络模型的预测性能；其中，峰值偏差b_p＝abs(peak_groundtruth-peak_predict)；peak_predict为数据片段的预测峰值点，peak_groundtruth为数据片段的实际峰值点。

Images