CN114692680A - 脑电信号特征处理方法及装置 - Google Patents
脑电信号特征处理方法及装置 Download PDFInfo
- Publication number
- CN114692680A CN114692680A CN202210257465.XA CN202210257465A CN114692680A CN 114692680 A CN114692680 A CN 114692680A CN 202210257465 A CN202210257465 A CN 202210257465A CN 114692680 A CN114692680 A CN 114692680A
- Authority
- CN
- China
- Prior art keywords
- matrix
- modal
- hand
- data
- sliding window
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 18
- 239000011159 matrix material Substances 0.000 claims abstract description 348
- 238000000034 method Methods 0.000 claims abstract description 119
- 230000008859 change Effects 0.000 claims abstract description 91
- 239000013598 vector Substances 0.000 claims abstract description 59
- 230000003190 augmentative effect Effects 0.000 claims abstract description 28
- 210000004556 brain Anatomy 0.000 claims abstract description 25
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 230000009466 transformation Effects 0.000 claims description 46
- 230000009467 reduction Effects 0.000 claims description 40
- 238000000354 decomposition reaction Methods 0.000 claims description 26
- 238000012545 processing Methods 0.000 claims description 17
- 238000012360 testing method Methods 0.000 claims description 14
- 230000003416 augmentation Effects 0.000 claims description 13
- 230000002087 whitening effect Effects 0.000 claims description 12
- 238000004590 computer program Methods 0.000 claims description 11
- 238000013507 mapping Methods 0.000 claims description 4
- 238000002360 preparation method Methods 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 3
- 239000000284 extract Substances 0.000 abstract description 8
- 230000008569 process Effects 0.000 description 18
- 230000007177 brain activity Effects 0.000 description 16
- 238000000537 electroencephalography Methods 0.000 description 16
- 230000006870 function Effects 0.000 description 15
- 238000010586 diagram Methods 0.000 description 12
- 238000000605 extraction Methods 0.000 description 10
- 238000001228 spectrum Methods 0.000 description 10
- 230000000694 effects Effects 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 7
- 230000001965 increasing effect Effects 0.000 description 7
- 230000033001 locomotion Effects 0.000 description 7
- 230000033764 rhythmic process Effects 0.000 description 6
- 238000012549 training Methods 0.000 description 6
- 238000005312 nonlinear dynamic Methods 0.000 description 5
- 230000000638 stimulation Effects 0.000 description 5
- 230000003925 brain function Effects 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 238000013528 artificial neural network Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 230000001427 coherent effect Effects 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 230000006399 behavior Effects 0.000 description 2
- 238000013527 convolutional neural network Methods 0.000 description 2
- 230000001054 cortical effect Effects 0.000 description 2
- 238000013016 damping Methods 0.000 description 2
- 238000013135 deep learning Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000008846 dynamic interplay Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000002599 functional magnetic resonance imaging Methods 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 230000001537 neural effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 210000000225 synapse Anatomy 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 206010008190 Cerebrovascular accident Diseases 0.000 description 1
- 240000002627 Cordeauxia edulis Species 0.000 description 1
- 241000764238 Isis Species 0.000 description 1
- 208000006011 Stroke Diseases 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 235000000332 black box Nutrition 0.000 description 1
- 210000003710 cerebral cortex Anatomy 0.000 description 1
- 230000002490 cerebral effect Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005183 dynamical system Methods 0.000 description 1
- 238000002566 electrocorticography Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 210000003128 head Anatomy 0.000 description 1
- 230000001939 inductive effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000004007 neuromodulation Effects 0.000 description 1
- 210000002569 neuron Anatomy 0.000 description 1
- 230000007996 neuronal plasticity Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 210000004761 scalp Anatomy 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems
- A61B5/7267—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems involving training the classification device
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Theoretical Computer Science (AREA)
- Pathology (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- General Health & Medical Sciences (AREA)
- Mathematical Physics (AREA)
- Signal Processing (AREA)
- Animal Behavior & Ethology (AREA)
- Psychiatry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Biophysics (AREA)
- General Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Physiology (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Fuzzy Systems (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Evolutionary Computation (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Psychology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明涉及一种脑电信号特征处理方法及装置，所述方法包括获取待测患者在预设时间段内的脑电信号，根据脑电信号得到脑电数据矩阵；以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵，并计算模态特征；选取与任务相关的模态向量，将模态向量沿时间进行拼接，得到模态变化信息矩阵并进行空间滤波，得到模态共空间模式；选取区分度最大的模态共空间模式滤波器，分别计算空间滤波后的左、右手运动想象对应模态共空间模式数据的方差信息作为最终特征进行分类。本发明可以基于较少的EEG时间序列数据提取与运动想象相关的全局动态特性，获得大脑模态变化的参数化精确描述，从而提高分类精度，减少脑机接口的控制延时。
Description
技术领域
本发明属于脑电信号检测技术领域，具体涉及一种脑电信号特征处理方法及装置。
背景技术
从脑电信号中提取任务相关的大脑活动特征有助于准确解码患者运动意图，进行康复运动。然而，由于脑电信号包含丰富信息并且是典型的非线性、非平稳信号，因此从脑电信号中提取有效的大脑活动特征是脑机接口中的最具挑战性的任务。现有的基于运动想象任务的脑电信号特征提取方法包括基于时域信息、频域信息、时频域和空间信息的方法。
其中，时域或频域方法，例如自回归(AR)模型方法、快速傅立叶变换(FFT)方法，能够分别从时域或频域表征大脑活动。时频方法则分析频率分量随时间的变化。短时傅立叶变换(STFT)使用固定长度的窗口来分析信号的瞬时频谱，给出了一种时频谱估计方法，而小波变换(WT)则通过部署多分辨率时间窗口提供了一种动态时频方法，被广泛用于脑电信号研究。尽管如此，基于单通道信息的时频谱方法无法反映脑电信号中不同通道之间的相互作用，且往往忽略大脑不同功能区域之间的相对相位关系，使得时频谱估计容易受噪声干扰，从而影响解码准确度和解码效率。
基于空间信息的方法也被广泛用于脑电特征提取。其中，共空间模式(CommonSpatial Pattern，CSP)及其变种算法广泛应用于运动想象任务的特征提取。CSP方法生成的空间滤波器通过同时最小化某一类别脑电信号在投影方向的方差和最大化另一类别的方差实现对不同任务的分类。共空间频谱模式(Common Spatial-Spectral Patterns，CSSP)在CSP空间滤波器的基础上进行了优化，通过对EEG信号插入时间延迟τ后再进行滤波，实现对频域特征的提取。共稀疏光谱空间模式(Common Sparse Spectral SpatialPattern，CSSSP)进一步改进了CSSP方法，不同于CSSP为每个通道计算各自的频谱，CSSSP计算所有通道共用的频谱模式，提高了特征的鲁棒性。子带共空间模式(Sub-band CommonSpatial Pattern,SBCSP)方法则在CSP和CSSP的基础上使用Gabor滤波器库，将脑电信号分解为多个子带。然后，在每个子带上使用CSP空间滤波，并融合各个子带分别提取的特征进行分类，实现了更高的分类精度。然而，SBCSP忽略了从不同子带获取的CSP特征的潜在关联性。滤波器组共空间模式(Filter Bank Common Spatial Pattern，FBCSP)方法通过计算SBCSP中来自多个子带的CSP特征的互信息，以便选择其中最具辨识性的子带特征用于分类。FBCSP相较于前面的几种方法取得了更好的表现，在实际场景中获得了更广泛的应用。
尽管基于CSP的方法在运动想象任务中取得了不错的效果，但其使用的信息仍比较粗糙，无法直接反映运动想象相关的脑功能连接变化，从而在跨个体分类任务中表现不理想。此外，传统的CSP方法基于投影数据的二阶统计量(方差信息)，需要基于较长的脑电数据段来获得准确的统计信息，因此在实际脑机接口应用中具有较长的控制延时。例如，基于3秒时间窗的CSP来说，需要待测患者在运动想象任务进行3秒后才可能给出分类结果，实时性无法满足在神经康复中快速诱导神经可塑性的需求。
随着深度学习的发展，可以将脑电信号直接输入到端到端的神经网络，自动提取特征并对信号进行分类。其中，卷积神经网络(CNN)常被用于提取空间或者时频谱特征，长短期记忆(LSTM)网络常被用于提取时间相关特征。深度学习方法在很多脑电分类任务中已经取得了很好的效果，但训练复杂的网络模型大数据依赖特性使得在有限数据下模型的泛化能力不佳。此外，神经网络与生俱来的黑盒特性，使得所使用的分类特征不透明，可解释性差，无法揭示与任务相关的大脑活动机制。
EEG特征提取的主要困难在于，EEG是一种非线性、非平稳信号，包含丰富的大脑活动信息，全面刻画大脑活动过程中与任务相关的时间、空间、频率、相位变化信息变得极为困难。Koopman算子理论及动态模态分解(DMD)方法，将复杂的非线性过程转变为观测空间内的线性动态过程，为应用成熟的线性系统方法研究复杂非线性动态提供了新的思路。该方法多应用于高维流体动力学研究，相关方法在脑电活动方面的研究仍比较少。
综上所述，现有的脑电信号特征处理方法存在较多缺陷。
发明内容
有鉴于此，本发明的目的在于克服现有技术的不足，提供一种脑电信号特征处理方法及装置，以解决现有技术中脑电信号特征处理方法存在较多缺陷的问题。
为实现以上目的，本发明采用如下技术方案：
一种脑电信号特征处理方法，包括：
获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
进一步的，所述获取待测患者在第一预设时间段内的脑电信号，包括：
获取左右手运动想象发生前第一预设时间点至发生后第二预设时间点的脑电信号。
进一步的，所述以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵，包括：
以预设的滑窗窗长、步长以及时域宽度确定滑窗窗口的个数；
以所述滑窗窗长、步长及滑窗窗口个数对所述脑电数据矩阵进行截取，得到多个沿时间采样的滑窗数据；
对每个滑窗数据进行时延堆叠，得到滑窗增广数据矩阵。
进一步的，所述对多个所述滑窗增广数据矩阵进行动态模态分解，提取模态特征，包括：
根据所述滑窗增广数据矩阵，获取第一观测数据矩阵和第二观测数据矩阵；其中，所述第二观测数据矩阵中的每一列状态值为所述第一观测数据矩阵中相应列的下一时刻的状态值；
对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵；
计算所述降维系统矩阵的特征值和特征向量，将所述降维系统矩阵的特征值确定为线性系统矩阵的特征值，将所述降维矩阵的特征向量映射回高维空间得到对应的线性系统矩阵的特征向量；
根据所述线性系统矩阵的特征值和特征向量，得到模态变化信息矩阵；其中，线性系统矩阵的特征值和特征向量表征所述线性系统矩阵的特征频率和对应不变子空间。
进一步的，所述对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵，包括：
确定低维空间的维数；其中所述维数表示奇异值个数；
选取所述第一观测数据矩阵和第二观测矩阵中截止至奇异值个数的奇异值和对应的奇异向量，构成第一降维变换矩阵和第二降维变换矩阵；
根据所述第一降维变换矩阵和第二降维变换矩阵，计算得到降维系统矩阵。
进一步的，所述选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵，包括：
从所述模态变化信息矩阵中选取关键模态对应的模态特征向量，计算所述模态特征向量中各元素的幅值和相位，根据所述幅值和相位构造该滑动窗口的模态幅相特征向量；其中，所述关键模态包括一个以上；
将一个以上的关键模态对应的模态幅相特征向量沿空间方向进行拼接，得到预设滑动窗口内的多模态特征向量；
依次对各滑动窗口计算对应的多模态特征向量，并沿时间方向进行拼接，得到模态变化信息矩阵。
进一步的，所述模态变化信息矩阵包括左手运动想象下的左手模态变化信息矩阵和右手运动想象下的右手模态变化信息矩阵；所述采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵，包括：
计算每次试验的左手模态变化信息矩阵和右手模态变化信息矩阵的归一化样本协方差矩阵，得到左手模态协方差矩阵和右手模态协方差矩阵；
对多次试验的左手模态协方差矩阵和右手模态协方差矩阵分别进行平均计算，得到左手平均模态协方差矩阵和右手平均模态协方差矩阵；
将所述左手平均模态协方差矩阵和右手平均模态协方差矩阵相加，得到和模态协方差矩阵，并对所述和模态协方差矩阵进行特征分解，以定义白化矩阵；
利用所述白化矩阵分别对左手平均模态协方差矩阵和右手平均模态协方差矩阵进行相似变换，得到左手协方差相似矩阵和右手协方差相似矩阵；
求解所述左手协方差相似矩阵和右手协方差相似矩阵的共同特征向量，以左手协方差相似矩阵对应特征值从大到小顺序，将对应的所述特征向量依次排列组成共同特征向量矩阵；
根据所述和模态协方差矩阵对应的白化矩阵和所述共同特征向量矩阵构造投影变换矩阵，即空间滤波器。
进一步的，所述投影变换矩阵，可同时对角化左手模态协方差矩阵、右手模态协方差矩阵，且其对应特征值之和为1。
本申请实施例提供一种脑电信号特征处理装置，包括：
获取模块，用于获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
数据准备模块，用于以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
模态分解模块，用于对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
空间滤波模块，用于采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
分类模块，用于应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
本申请实施例提供一种计算机设备，包括：存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行上述任一项脑电信号特征处理方法的步骤。
本申请实施例还提供一种计算机存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行上述任一项脑电信号特征处理方法的步骤。
本发明采用以上技术方案，能够达到的有益效果包括：
本发明提供一种脑电信号特征处理方法及装置，本申请将DMD方法应用于脑电模态分析，揭示在运动想象过程中脑功能连接的变化情况，进而用于运动想象分类任务。与现有的脑电特征提取方法相比，基于DMD的方法有诸多优点：首先，DMD方法基于脑电通道之间的动态交互关系来提取脑电模式特征，相较于传统逐通道单独提取特征的方法考虑了通道之间的动态约束关系，使得提取的特征更为稳定；其次，DMD方法基于数据驱动方法提取脑电活动中频率成分，相较于基于固定脑电节律范围的方法，可以更准确反映不同个体间主要频率成分的个体差异性；第三，每一个表征时空关联信息的脑电模态，其对应的时空相干模式反映了大脑不同区域之间在该节律的功能连接关系，即包括脑活动的相对强度(幅值信息)，也包括各区域之间大脑活动的时序差异(相位信息)；另外，DMD方法的本质上是一种多变量自回归模型，相较于傅里叶变换等非参数方法，可以基于较少的数据获得动态过程的参数化精确描述，从而减少脑机接口的控制延时；最后，依据Koopman算子理论，DMD方法将脑电信号作为非线性脑活动的高维观察变量，获得的脑模态从线性观测空间中表征了大脑的非线性动力学行为。
除此之外，本申请所提出的DMD-CSP方法在较短的延迟下可以同时提高个体内和个体间的分类性能。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明脑电信号特征处理方法的步骤示意图；
图2为本发明脑电信号特征处理方法的流程示意图；
图3为本发明脑电信号特征处理装置的结构示意图；
图4为本发明脑电信号特征处理方法的步骤示意图；
图5为本发明提供的模态数量为11时动态模态频谱图；
图6为本发明提供的模态数量为11时快速傅里叶变换频谱图；
图7为本发明提供的脑电信号特征处理装置的结构示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。
根据Hebb学习理论，神经元突触在经历多次重复刺激后其连接性得到加强，且神经可塑性与刺激延时密切相关。因此，研究与运动相关的皮层突触之间的相互作用并降低神经调控的延迟时间对于脑机接口在神经康复中的应用至关重要。
脑机接口(BCI)作为脑卒中康复过程中的重要方法，与传统的康复方法不同，无需让患者进行身体部位执行实际运动，通过解码患者的运动意图生成控制信号，与患者进行基于运动想象的双向交互：一方面通过输出指令控制外设，另一方面通过调节大脑皮层运动区域内的神经活动可塑性，促进患者的康复水平。不同于采集成本较高的皮层电位图(ECoG)和功能性磁共振成像(fMRI)技术，相对经济便携的采集过程使得脑电图(EEG)成为BCI研究中最常用的信息获取方法，这是因为EEG具有很高的时间分辨率，脑电信号尤其适合于大脑动态活动和功能网络的研究。
根据Koopman算子理论，一个有限维的非线性动态系统可以表示为无限维观测空间下的线性过程。动态模态分解(DMD)方法提供了上述无限维观测空间的一种有限维近似。根据线性系统理论，对于离散的线性系统xk+1＝Axk，其动态进化过程可分解为：
xk＝ΨΛkB0
其中，Ψ是矩阵A的特征向量，Λ是矩阵A的特征值，
是参与度系数，即初始状态在特征空间上的投影。如果系统的特征值为实数，那么它所对应的特征向量是指数级增长或衰退的；相反，如果得到的特征值是复数，那么它对应的模态则存在振荡。因此可以发现，矩阵A的特征向量决定了系统的模态，而A的特征值决定了模态的频率以及增长率。
本申请中将N个脑电测量通道视为当前分析时刻大脑活动状态的N个观测函数，DMD的目标则是在观测函数空间找到一个满足于以下近似关系的线性系统：
根据Koopman算子理论，通过选择合适的标量可观测状态函数空间，可以将有限维的非线性动力系统转化为无限维的线性动力系统问题。该线性系统定义为：
其中，
是Koopman算子，g：Rn→C为观测函数。
大脑活动是一个非线性动态系统，脑电信号是通过采集设备从头皮上测量而得。因此脑电信号被视为运动想象过程中大脑状态函数的观测值，即g(x(t))＝x(t)，
需要说明的是，上述的工作就是在确定以原始脑电信号作为观测函数以后，通过DMD算法来近似Koopman算子的特征值和特征函数。
由于EEG信号的强非线性特征，对应的Koopman算子具有丰富的特征谱，基于线性“全状态”观测函数的DMD方法，即仅依赖于当前脑电通道数据构造的观测空间，很难确保是在Koopman算子作用下的不变空间，从而难以准确地逼近Koopman算子的特征值和特征函数。而Takens时延嵌套理论表明，一个动态系统可以用足够多的历史信息精确描述。因此时延嵌套可以有效增广观测函数空间，即利用过去历史的测量值来增广观测函数，产生一个更为丰富的新的观测函数空间
其中△t是延时时间。
下面结合附图介绍本申请实施例中提供的一个具体的脑电信号特征处理方法及装置。
如图1所示，本申请实施例中提供的脑电信号特征处理方法，包括：
S101，获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
S102，以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
S103，对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
S104，采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
S105，应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
脑电信号特征处理方法的工作原理为：如图2所示，患者分别进行左手运动想象和右手运动想象，记录待测患者运动想象时间段T内(想象开始前0.3秒至开始后1.3秒)的N通道脑电信号，得到N×T的脑电数据矩阵。从该N×T的脑电数据矩阵沿时序方向应用滑动窗口，依次截取窗长为L的脑电数据。具体的，设滑窗步长为d，则可得到
个沿时间采样的滑动窗口。对每个滑窗数据进行时延堆叠，得N′×L′的滑窗增广数据矩阵。对每个滑窗增广数据矩阵进行DMD分析，提取出数据模态，并从中选择与运动想象活动相关的q个频率成分所对应模态的幅值和相位作为特征，这些特征反映了瞬时脑电模态及对应不同节律脑活动的功能连接特性。由于不同的滑动窗口表征了各脑功能连接随时间的变化情况，因此将随时间推移的不同滑窗对应的模态特征进行拼接，得到2qN×p的模态变化信息矩阵。之后使用共空间模式对该模态变化信息矩阵进行空间滤波，选取左右手特征值差异最大的投影作为投影变换矩阵即空间滤波器，提高分类准确度。在应用空间滤波器的时候是将待解码模态变化信息矩阵采用空间滤波器进行滤波，得到模态共空间模式数据，根据模态共空间模式数据得到滤波后的模态变化信息沿时间方向的方差，最后应用线性判别分析进行分类。
本申请提供的技术方案提取的特征刻画了核心频率时空关联的全脑功能连接随时间变化关系，同时包含时间、空间、频率、相位信息，提升运动想象任务的分类性能。
一些实施例中，所述获取患者在第一预设时间段内的脑电信号，包括：
获取左右手运动想象发生前第一预设时间点至发生后第二预设时间点的脑电信号。
具体的，如图3所示，选取运动想象发生前0.3秒至发生后1.3秒的预处理后的EEG数据作为样本(左手运动想象与右手运动想象)，得到N通道脑电信号数据矩阵N×T。由于DMD方法需要对一定时长的脑电信号进行分析才能得到模态特征，因此需要对该1.6秒的数据滑窗采样成不同子段，用各个窗内子段的特征来表征该时刻的系统状态。
一些实施例中，所述以滑动窗口对所述脑电数据矩阵进行截取，得到多个滑窗数据矩阵，包括：
以预设的滑窗窗长、步长以及时域宽度确定滑窗窗口的个数；
以所述滑窗窗长、步长及滑窗窗口个数对所述脑电数据矩阵进行截取，得到多个沿时间采样的滑窗数据；
对每个滑窗数据进行时延堆叠，得到多个滑窗增广数据矩阵。
例如：滑窗窗长为L，步长为d，在时域宽度T内提取的滑窗个数为
得到p个N×L的滑窗数据矩阵。滑窗N×L可以表示为：
在实际BCI应用中使用较少的数据可以提升MI运动意图解码的实时性并减少控制延时，因此本申请提供的方法仅使用较少的数据片段来进行MI分类。例如选取窗长为0.6秒，选取步长为0.0625秒，则一共产生17个滑动窗口。
为使观测数据更好的满足Koopman算子作用下的空间不变性，根据Takens延时嵌套理论，当延时嵌套维度s需满足s≥2N+1时，其中N是原状态维度，动态过程可以由历史数据线性表出，即满足如下线性延时嵌套关系：
因此，我们对原来的滑动窗口内的数据进行最大延时为s的堆叠。
具体的，在本申请中对每一个N×L滑窗数据矩阵按嵌入维度s进行时延堆叠，得到对应滑窗增广数据矩阵，所述滑窗增广数据矩阵为如下Hankel矩阵：
一些实施例中，所述对多个所述滑窗增广数据矩阵进行动态模态分解，提取模态特征，包括：
根据所述滑窗增广数据矩阵，获取第一观测数据矩阵和第二观测数据矩阵；其中，所述第二观测数据矩阵中的每一列状态值为所述第一观测数据矩阵中相应列的下一时刻的状态值；
对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵；
计算所述降维系统矩阵的特征值和特征向量，将所述降维系统矩阵的特征值确定为线性系统矩阵的特征值，将所述降维矩阵的特征向量映射回高维空间得到对应的线性系统矩阵的特征向量；
根据所述线性系统矩阵的特征值和特征向量，得到模态变化信息矩阵；其中，线性系统矩阵的特征值和特征向量表征所述线性系统矩阵的特征频率和对应不变子空间。
优选的，所述对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵，包括：
确定低维空间的维数；其中所述维数表示奇异值个数；
选取所述第一观测数据矩阵和第二观测矩阵中截止至奇异值个数的奇异值和对应的奇异向量，构成第一降维变换矩阵和第二降维变换矩阵；
根据所述第一降维变换矩阵和第二降维变换矩阵，计算得到降维系统矩阵。
一些实施例中，所述选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵，包括：
从所述模态变化信息矩阵中选取关键模态对应的模态特征向量，计算所述模态特征向量中各元素的幅值和相位，根据所述幅值和相位构造该滑动窗口的模态幅相特征向量；其中，所述关键模态包括一个以上；
将一个以上的关键模态对应的模态幅相特征向量沿空间方向进行拼接，得到预设滑动窗口内的多模态特征向量；
依次对各滑动窗口计算对应的多模态特征向量，并沿时间方向进行拼接，得到模态变化信息矩阵。
具体的，类似的构造下一时刻的Hankel矩阵得到两个增广观测数据矩阵：
即，第二观测矩阵X′的每一列是第一观测矩阵X的每一列的下一时刻的状态值。
根据Koopman算子理论，满足离散的线性系统X′＝AX，其中线性系统矩阵A描述了在观测空间中系统状态从X到状态X′的线性演化过程，其最直接的近似解为
其中
是矩阵X的广义逆矩阵。
根据线性系统理论，线性系统矩阵A的特征值和特征向量给出了该线性系统的特征频率和对应不变子空间，称为DMD模态。当观测空间维度比较高时，求解A矩阵的问题为不适定(ill-posed)问题。因此DMD方法首先对数据降维，通过降维后空间上的相似系统矩阵
间接计算矩阵A的特征向量和特征值，从而提取系统中支配该线性过程的主要模态。
对第一观测矩阵X进行奇异值个数为r的截断奇异值分解，得到
X≈UΣV*
其中，*表示共轭转置，提取其前r个奇异值和对应奇异向量组成矩阵
将矩阵X的秩由N减少到r。参数r的选择至关重要，r过小或者r过大都会造成动态信息的损失或者引入冗余信息，使得模态信息不准确。
构造低维度空间下系统状态的的线性演化算子(也称降维系统矩阵)：
即
为r维空间下的线性系统矩阵，也就是降维系统矩阵。
计算降维系统矩阵
的特征值和特征向量，有
其中，矩阵Λ为对角阵，对角元素为
的特征值，其对应的特征向量为矩阵M中的列向量。根据相似矩阵的性质，矩阵
的特征值同时也是矩阵A的特征值，而矩阵A对应特征向量矩阵Ψ可以通过将特征向量M映射回高维空间中的对应向量得到：
Ψ[N,r]＝Λ-1YV∑-1M
可以理解的是，通过上述DMD方法可以得到观察空间线性系统矩阵A(近似Koopman算子)的离散谱，其中每个模态包含特征值和特征向量，特征值描述了对应频率成分的频率及稳定性，特征向量描述了空间相干(coherence)关系，以及如何重构原始信号。
通过特征值矩阵Λ的对角元素λ得到对应的模态的阻尼系数和特征频率：
g＝real(ω)
其中，ω为该模态对应的角频率，Δt为样本点间的采样间隔，g为计算得到的阻尼系数，f为模态对应的特征频率。当我们得到线性系统的系统矩阵A的特征值和特征向量后，可以用动态模态分解得到的模态重构原信号：
其中，r为模态数量，Ψ为由ψk构成的模态变化信息矩阵，Ω为由ωk构成的矩阵，b为每个模态的系数。
由此得到系统的DMD模态，如图4所示，对应于不同的特征频率，模态振型可能是实向量，也可能是复向量。当它为实模态时，代表各空间位置点总是以相同或者相反的相位振动，各点振动形成空间上的驻波，各点相位差为0度或180度，模态振型向量不同的值代表各点的相对幅度和方向；而当该模态为复模态时，模态振型向量为复数，系统内各观测点存在相位差，各点振动形成空间上的行波，幅值表示各点的振幅大小，相位反映振动发生的先后顺序。不论实模态或复模态，在同一模态内，各个位置点均按照对应特征频率ωk振动。
因此，本申请中使用时延堆叠的动态模态分解，得到模态变化信息矩阵：
Ψ[N×s,r]＝UM＝XV∑-1M
取上述模态特征向量矩阵Ψ[N×s,r]的后N行作为最终的模态变化信息矩阵Ψ[N,r]。
需要说明的是，本申请中计算的到的DMD模态由：模态的特征值、模态振型向量、参与度因子三元组组成。模态的特征值反映了模态的自然频率及稳定性；模态的振型反映了空间相干性；参与度因子反映了各个模态在指定信号中的占比。
对于每个脑电滑窗中的数据，应用本算法提取新的模态特征。因此对于一个信号样本片段中的一系列滑窗可以生成一组包含多通道之间功能连通性变化的特征。由于这种新的特征包含各个通道在所关心频率能量和相位随时间的变化，因此这种全新的模态特征向量包含脑连接的时空变化特征。最后采用CSP方法作为空间滤波器以提取每个样本中脑功能连接变化的方差。
DMD复模态反映了跨通道的全脑动态特性，因为多通道脑电信号在计算DMD的过程中被视为一个相互耦合的多变量系统。EEG信号可以被分解为具有对应DMD模态形状的离散谱，其中DMD模态是复向量，代表对应节律的相干模式(功能连接)。其幅值反映了特定频率下不同通道的幅度差异，其相位描述了脑电中不同通道的连接性。
需要说明的是，在动态模态分解步骤中，当选择不同数量的奇异值时，我们可以分解出不同的频率成分。当选择数量较少的奇异值时，分解出的频率成分对于数据矩阵更为重要，影响更大，而当选择奇异值数量逐渐增多时，新出现的频率成分对于数据矩阵来说相对不重要。因此，我们可以先不考虑信号的还原程度，从而确定在μ波段内各个频率之间的重要性差异。如图5所示，模态数量为3时，只有直流成分和频率为11Hz的成分被提取出来，因此可以发现11Hz的频率成分是最重要的，因此我们接下来的模态分析主要针对11Hz，该频率也和有关于运动想象的大脑活动研究结论相符。
而当模态数量提高为11时，可以发现频率为8Hz和频率为11Hz的信号都被提取出来，如图6所示。由此可见，8Hz频率重要性低于11Hz，但都较8-12Hz内的其他频率更为重要。因此我们为了保证信号的还原程度，选择模态数量依然为51，但是在提取信息时，只着重关注8Hz与11Hz的信息，避免其他频率信息对于结果产生无关的干扰。
综上，通过逐步在动态模态分解中增加奇异值选取数量的方法，可以预分析得到该脑电任务最为关键的频率及其对应的模态。从r个模态中将这些最为关键的模态(q个)选取出来，用于后续共空间模式空间滤波，可以有效提高特征的代表性和鲁棒性。其中，模态个数与奇异值个数相同。
一些实施例中，所述模态变化信息矩阵包括左手运动想象下的左手模态变化信息矩阵和右手运动想象下的右手模态变化信息矩阵；所述采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；包括：
计算每次试验的左手模态变化信息矩阵和右手模态变化信息矩阵的归一化样本协方差矩阵，得到左手模态协方差矩阵和右手模态协方差矩阵；
对多次试验的左手模态协方差矩阵和右手模态协方差矩阵分别进行平均计算，得到左手平均模态协方差矩阵和右手平均模态协方差矩阵；
将所述左手平均模态协方差矩阵和右手平均模态协方差矩阵相加，得到和模态协方差矩阵，并对所述和模态协方差矩阵进行特征分解，以定义白化矩阵；
利用所述白化矩阵分别对左手平均模态协方差矩阵和右手平均模态协方差矩阵进行相似变换，得到左手协方差相似矩阵和右手协方差相似矩阵；
求解所述左手协方差相似矩阵和右手协方差相似矩阵的共同特征向量，依左手协方差相似矩阵对应特征值从大到小顺序，将对应的所述特征向量依次排列组成共同特征向量矩阵；
根据所述和模态协方差矩阵对应的白化矩阵和所述共同特征向量矩阵构造投影变换矩阵，即空间滤波器。
具体的，经动态模态分解得到模态变化信息矩阵Ψ[N,r]，从r个系统模态中选取q个与运动想象活动相关的关键频率成分对应的模态，分别计算其各元素的幅值和相位(复数)，得到每个滑窗下的关键复模态向量的模值和相位模态变化信息矩阵：
将幅值和相位矩阵按列向量堆叠为该滑动窗口下的系统幅相特征向量：
将p个滑窗下的系统幅相特征向量组合，得到该运动想象时间段内的时频幅相模态变化信息矩阵，即模态变化信息矩阵为
使用共空间模式对其进行空间滤波，其原理本质上是要找到一个投影变换矩阵W，使得
PL、PR分别为左手模态变化信息矩阵和右手模态变化信息矩阵，共空间模式即通过找到一个线性投影变换，在投影空间内同时最大化不同类别之间方差差异，从而有利于提高分类的性能。
计算该投影变换矩阵的方法为：
首先对每一个Trial计算左右手运动想象模态变化信息矩阵的归一化协方差矩阵：
对左右手所有Trial计算出的协方差矩阵分别进行平均，得到左手平均模态协方差矩阵
右手平均模态协方差矩阵
定义和模态协方差矩阵为
对其进行特征分解，得到
R＝UDUT
其中U的列为为矩阵的单位特征向量，D为对角阵，其元素为矩阵对应的特征值。
定义白化变换矩阵：
用白化矩阵对协方差矩阵进行投影，得到
其中，SL为左手协方差相似矩阵，SR为右手协方差相似矩阵。
可以证明，上述矩阵SL与SR有相同的特征向量，也就是共同特征向量Q，并且特征值对应的对角阵之和为单位阵，即对应的特征值之和为1：
SL＝QDLQT
SR＝QDRQT
DL+DR＝I
因此，共同特征向量B中对应于SL的最大特征值是对应于SR的最小特征值，即当左手协方差相似矩阵SL的特征值最大时，右手协方差相似矩阵SR的特征值最小，由此达到对两类任务的最优分离的目的。其中，所述投影变换矩阵，可同时对角化左手模态协方差矩阵、右手模态协方差矩阵，且其对应特征值之和为1。可以理解的是，共同特征向量矩阵可以同时对角化两个白化处理以后的协方差矩阵，且其对应特征值之和为1。
投影变换矩阵由此定义为：
W＝QTF
其中，W即为空间滤波器，将待解码模态变化信息矩阵输入到空间滤波器中进行滤波，得到模态共空间模式数据，具体如下：
将模态变化信息矩阵投影得到最优特征空间：
ZL＝WPL
ZR＝WPR
其中，ZL为左手共空间模式数据，ZR为右手共空间模式数据。
算法生成的CSP模态变化信息矩阵，其信息不是等效的。特征信息主要集中在模态变化信息矩阵的头部和尾部，而中间的特征信息不明显可以忽略，所以选取m行和后m行数据作为CSP特征提取的模态变化信息矩阵。
在该特征空间中，特征的方差区分性较好，取共空间模式数据的前m行和后m行提取方差对数FL和FR作为最终分类特征：
FL＝ln Var(ZL)
FR＝ln Var(ZR)
其中，FL为左手共空间模式数据的方差，FR为右手共空间模式数据的方差。
最后，将上述提取到的最终特征送入线性判别分类器LDA，从而对左手右手运动想象进行分类。LDA的方法是在投影空间中寻找一种最优投影以使得同类样本之间的距离较小、不同类样本之间的距离较大。在确定分类器模型参数后，相较于复杂神经网络，LDA对新的脑电信号具有更快的判别速度，可以提高BCI的实时性能。
本方法融合Koopman算子理论和Takens嵌入理论，将DMD方法应用于脑电模态分析，揭示在运动想象过程中脑功能连接的变化情况，进而用于运动想象分类任务。与现有的脑电特征提取方法相比，基于DMD的方法有诸多优点：首先，DMD方法基于脑电通道之间的动态交互关系来提取脑电模式特征，相较于传统逐通道单独提取特征的方法考虑了通道之间的动态约束关系，使得提取的特征更为稳定；其次，DMD方法基于数据驱动方法提取脑电活动中频率成分，相较于基于固定脑电节律范围的方法，可以更准确反映不同个体间主要频率成分的个体差异性；第三，每一个表征时空关联信息的脑电模态，其对应的时空相干模式反映了大脑不同区域之间在该节律的功能连接关系，即包括脑活动的相对强度(幅值信息)，也包括各区域之间大脑活动的时序差异(相位信息)；另外，DMD方法的本质上是一种多变量自回归模型，相较于傅里叶变换等非参数方法，可以基于较少的数据获得动态过程的参数化精确描述，从而减少脑机接口的控制延时；最后，依据Koopman算子理论，DMD方法将脑电信号作为非线性脑活动的高维观察变量，获得的脑模态从线性观测空间中表征了大脑的非线性动力学行为。本方法在Physionet EEG运动图像和BCI竞赛IV 2a数据集上进行了验证。结果表明，所提出的DMD-CSP方法在较短的延迟下可以同时提高个体内和个体间的分类性能。
本申请使用了Physionet EEG运动想象数据集(eegmmidb)以及BCI CompetitionIV 2a数据集对所提DMD-CSP算法的有效性进行验证。
在eegmmidb数据集中，每个待测患者有45次试验，包括1：1的左手运动想象和右手运动想象。由于样本数量较少，通过采取滑动窗口处理进行数据集增强。对于每次运动想象任务，在运动图像视觉刺激信号发出后，使用偏移量为10个采样点的滑动窗口获得5个长度为1.3秒的数据段，每个待测患者共生成225个样本。将数据集分为训练集和测试集，其中36次试验的数据为训练集和9次试验的数据为测试集；因此一共有180个训练样本和45个测试样本。对每个训练样本，采用DMD-CSP特征提取方法生成时空模态变化信息矩阵，并根据CSP抽取方差差别最大的前后各3个投影方向，共提取6个特征。最后将该特征送入线性判别分类器中进行分类，每个待测患者的最终分类准确率为运行10次实验的平均分类准确率，不同特征提取方法的对比结果见表1。
表1：eegmmidb数据集上待测患者内不同分类方法的对比结果
可以看出，所提算法产生的平均分类准确率最高。DMD-CSP方法对于待测患者2(增加16.6％)、待测患者3(增加21.06％)待测患者7(增加24.67％)和待测患者9(增加16.14％)的分类性能具有显著提升。
对于BCI Competition IV 2a数据集，每个待测患者有包含144次试验的训练集和包含144次试验的测试集。同样适用刺激前0.3秒到刺激后1秒的数据段被用于训练分类器，不同分类方法的分类比较结果见表2。DMD-CSP方法的平均分类精度高于相应的CSP方法。
值得注意的是BCI Competition IV 2a数据集的平均精度低于eegmmidb数据集。其主要原因在于该数据集中的EEG通道数量较少(22通道)，可见DMD-CSP更适合于多通道脑电测量数据。其原因在于DMD-CSP方法在模态空间中提取全脑的时空特征，基于Koopman理论，更多的通道意味着在更高维的观测空间中研究大脑的非线性动态变化。
表2：BCI Competition IV 2a数据集上待测患者内不同分类方法的对比结果
跨待测患者分类性能，具体如下：
为了进一步验证DMD-CSP方法在跨待测患者MI分类任务中的表现性能，本部分实验采用留一交叉验证，从数据集选取9名待测患者数据，采用8个待测患者的数据作为训练集，剩余1名待测患者的EEG信号作为测试集。表3列出了本算法在eegmmidb数据集上于其他方法的对比结果。结果表明，该方法总体上优于现有方法。具体地，相较于CSP平均精度提高了8.79％，相较于FBCSP提高了4.51％。与其他方法相比，该方法在受试者8的性能上有了显著的提高(提高了10.8％)。表4显示了在BCI Competition IV 2a数据集中测试的比较结果。该方法仍然可以提高平均分类精度，特别是受试者9的性能(提高6.46.％)。
表3：eegmmidb数据集上跨待测患者不同分类方法的对比结果
表4：BCI Competition IV 2a数据集上跨待测患者不同分类方法的对比结果
如图7所示，本申请提供一种脑电信号特征处理装置，包括：
获取模块701，用于获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
数据准备模块702，用于以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
模态分解模块703，用于对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
空间滤波模块704，用于采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
分类模块705，用于应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
本申请提供的脑电信号特征处理装置的工作原理为，获取模块701获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；数据准备模块702以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；模态分解模块703对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；空间滤波模块704采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；分类模块705应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
综上所述，本发明提供一种脑电信号特征处理方法及装置，本发明可以基于较少的数据获得动态过程的参数化精确描述，从而减少脑机接口的控制延时，能够从较短的EEG时间序列片段中提取与运动想象相关的全局的动力学特性进行分类。
可以理解的是，上述提供的方法实施例与上述的装置实施例对应，相应的具体内容可以相互参考，在此不再赘述。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令方法的制造品，该指令方法实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
Claims (10)
1.一种脑电信号特征处理方法，其特征在于，包括：
获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
2.根据权利要求1所述的方法，其特征在于，所述获取待测患者在第一预设时间段内的脑电信号，包括：
获取左右手运动想象发生前第一预设时间点至发生后第二预设时间点的脑电信号。
3.根据权利要求1或2所述的方法，其特征在于，所述以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵，包括：
以预设的滑窗窗长、步长以及时域宽度确定滑窗窗口的个数；
以所述滑窗窗长、步长及滑窗窗口个数对所述脑电数据矩阵进行截取，得到多个沿时间采样的滑窗数据；
对每个滑窗数据进行时延堆叠，得到滑窗增广数据矩阵。
4.根据权利要求1所述的方法，其特征在于，所述对多个所述滑窗增广数据矩阵进行动态模态分解，提取模态特征，包括：
根据所述滑窗增广数据矩阵，获取第一观测数据矩阵和第二观测数据矩阵；其中，所述第二观测数据矩阵中的每一列状态值为所述第一观测数据矩阵中相应列的下一时刻的状态值；
对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵；
计算所述降维系统矩阵的特征值和特征向量，将所述降维系统矩阵的特征值确定为线性系统矩阵的特征值，将所述降维矩阵的特征向量映射回高维空间得到对应的线性系统矩阵的特征向量；
根据所述线性系统矩阵的特征值和特征向量，得到模态变化信息矩阵；其中，线性系统矩阵的特征值和特征向量表征所述线性系统矩阵的特征频率和对应不变子空间。
5.根据权利要求4所述的方法，其特征在于，对所述第一观测矩阵和第二观测矩阵进行数据降维处理，得到对应降维系统矩阵，包括：
确定低维空间的维数；其中所述维数表示奇异值个数；
选取所述第一观测数据矩阵和第二观测矩阵中截止至奇异值个数的奇异值和对应的奇异向量，构成第一降维变换矩阵和第二降维变换矩阵；
根据所述第一降维变换矩阵和第二降维变换矩阵，计算得到降维系统矩阵。
6.根据权利要求4所述的方法，其特征在于，所述选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵，包括：
从所述模态变化信息矩阵中选取关键模态对应的模态特征向量，计算所述模态特征向量中各元素的幅值和相位，根据所述幅值和相位构造该滑动窗口的模态幅相特征向量；其中，所述关键模态包括一个以上；
将一个以上的关键模态对应的模态幅相特征向量沿空间方向进行拼接，得到预设滑动窗口内的多模态特征向量；
依次对各滑动窗口计算对应的多模态特征向量，并沿时间方向进行拼接，得到模态变化信息矩阵。
7.根据权利要求6所述的方法，其特征在于，所述模态变化信息矩阵包括左手运动想象下的左手模态变化信息矩阵和右手运动想象下的右手模态变化信息矩阵；所述采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵，包括：
计算每次试验的左手模态变化信息矩阵和右手模态变化信息矩阵的归一化样本协方差矩阵，得到左手模态协方差矩阵和右手模态协方差矩阵；
对多次试验的左手模态协方差矩阵和右手模态协方差矩阵分别进行平均计算，得到左手平均模态协方差矩阵和右手平均模态协方差矩阵；
将所述左手平均模态协方差矩阵和右手平均模态协方差矩阵相加，得到和模态协方差矩阵，并对所述和模态协方差矩阵进行特征分解，以定义白化矩阵；
利用所述白化矩阵分别对左手平均模态协方差矩阵和右手平均模态协方差矩阵进行相似变换，得到左手协方差相似矩阵和右手协方差相似矩阵；
求解所述左手协方差相似矩阵和右手协方差相似矩阵的共同特征向量，以左手协方差相似矩阵对应特征值从大到小顺序，将对应的所述特征向量依次排列组成共同特征向量矩阵；
根据所述和模态协方差矩阵对应的白化矩阵和所述共同特征向量矩阵构造投影变换矩阵，即空间滤波器。
8.根据权利要求7所述的方法，其特征在于，
所述投影变换矩阵，可同时对角化左手模态协方差矩阵、右手模态协方差矩阵，且其对应特征值之和为1。
9.一种脑电信号特征处理装置，其特征在于，包括：
获取模块，用于获取待测患者在预设时间段内的脑电信号，根据所述脑电信号得到脑电数据矩阵；所述脑电信号根据待测患者进行的左手和右手运动想象生成；
数据准备模块，用于以滑动窗口对所述脑电数据矩阵进行截取并进行时延堆叠，得到滑窗增广数据矩阵；
模态分解模块，用于对多个随时间变化的所述滑窗增广数据矩阵分别进行动态模态分解，提取模态特征，选取任务相关的模态向量，并按照时间推移对各滑窗增广数据矩阵对应的模态向量进行拼接，得到模态变化信息矩阵；
空间滤波模块，用于采用共空间模式方法计算所述模态变化信息矩阵的模态共空间模式，并得到相应投影变换矩阵；所述投影变换矩阵即为空间滤波器；
分类模块，用于应用所述空间滤波器依次对待解码模态变化信息矩阵进行空间滤波，得到模态共空间模式数据，并根据所述模态共空间数据计算滤波后的模态变化信息沿时间方向的方差，根据所述方差确定为最终分类特征并对所述最终分类特征进行分类。
10.一种计算机设备，包括：存储器和处理器；
所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述权利要求1至8任一项所述的脑电信号特征处理方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210257465.XA CN114692680A (zh) | 2022-03-16 | 2022-03-16 | 脑电信号特征处理方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210257465.XA CN114692680A (zh) | 2022-03-16 | 2022-03-16 | 脑电信号特征处理方法及装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114692680A true CN114692680A (zh) | 2022-07-01 |
Family
ID=82138638
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210257465.XA Pending CN114692680A (zh) | 2022-03-16 | 2022-03-16 | 脑电信号特征处理方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114692680A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116341616A (zh) * | 2022-11-11 | 2023-06-27 | 南京工程学院 | 基于矩阵重构二维卷积网络的电力负载信息获取方法 |
| CN116491960A (zh) * | 2023-06-28 | 2023-07-28 | 南昌大学第一附属医院 | 脑瞬态监测设备、电子设备及存储介质 |
| CN116942184A (zh) * | 2023-07-24 | 2023-10-27 | 山东睿芯半导体科技有限公司 | 一种脑电图生物特征核验方法、装置、芯片及终端 |
| CN117455013A (zh) * | 2023-11-10 | 2024-01-26 | 无锡鸣石峻致医疗科技有限公司 | 一种训练样本数据生成方法、系统、电子设备及介质 |
-
2022
- 2022-03-16 CN CN202210257465.XA patent/CN114692680A/zh active Pending
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116341616A (zh) * | 2022-11-11 | 2023-06-27 | 南京工程学院 | 基于矩阵重构二维卷积网络的电力负载信息获取方法 |
| CN116341616B (zh) * | 2022-11-11 | 2023-10-17 | 南京工程学院 | 基于矩阵重构二维卷积网络的电力负载信息获取方法 |
| CN116491960A (zh) * | 2023-06-28 | 2023-07-28 | 南昌大学第一附属医院 | 脑瞬态监测设备、电子设备及存储介质 |
| CN116491960B (zh) * | 2023-06-28 | 2023-09-19 | 南昌大学第一附属医院 | 脑瞬态监测设备、电子设备及存储介质 |
| CN116942184A (zh) * | 2023-07-24 | 2023-10-27 | 山东睿芯半导体科技有限公司 | 一种脑电图生物特征核验方法、装置、芯片及终端 |
| CN116942184B (zh) * | 2023-07-24 | 2024-04-26 | 山东睿芯半导体科技有限公司 | 一种脑电图生物特征核验方法、装置、芯片及终端 |
| CN117455013A (zh) * | 2023-11-10 | 2024-01-26 | 无锡鸣石峻致医疗科技有限公司 | 一种训练样本数据生成方法、系统、电子设备及介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114692680A (zh) | 脑电信号特征处理方法及装置 | |
| Lakshmi et al. | Survey on EEG signal processing methods | |
| Ince et al. | Classification of single trial motor imagery EEG recordings with subject adapted non-dyadic arbitrary time–frequency tilings | |
| Ince et al. | Adapting subject specific motor imagery EEG patterns in space–time–frequency for a brain computer interface | |
| Zhang et al. | Bayesian learning for spatial filtering in an EEG-based brain–computer interface | |
| Singh et al. | Small sample motor imagery classification using regularized Riemannian features | |
| CN110781945A (zh) | 一种融合多种特征的脑电信号情感识别方法及系统 | |
| CN114533086A (zh) | 一种基于空域特征时频变换的运动想象脑电解码方法 | |
| Katthi et al. | Deep correlation analysis for audio-EEG decoding | |
| CN113967022B (zh) | 一种基于个体自适应的运动想象脑电特征表征方法 | |
| Ramos-Aguilar et al. | Analysis of EEG signal processing techniques based on spectrograms | |
| CN111310656A (zh) | 基于多线性主成分分析的单次运动想象脑电信号识别方法 | |
| KR102443961B1 (ko) | 뇌파를 활용한 동작 상상 분류 장치 및 방법 | |
| CN115414051A (zh) | 一种脑电信号自适应窗口的情绪分类识别方法 | |
| Khare et al. | Multiclass sleep stage classification using artificial intelligence based time-frequency distribution and CNN | |
| CN112869743B (zh) | 一种考虑认知分心的运动起始意图神经解析方法 | |
| Ziehe | Blind source separation based on joint diagonalization of matrices with applications in biomedical signal processing | |
| Elgharabawy et al. | Decoding of finger movement using kinematic model classification and regression model switching | |
| CN117503157A (zh) | 一种基于sgcrnn模型的脑电信号情绪识别方法 | |
| Seha et al. | A new training approach for deep learning in EEG biometrics using triplet loss and EMG-driven additive data augmentation | |
| Mardiansyah et al. | Multivariate eeg signal using pca and cnn in post-stroke classification | |
| Wang | Temporally local maximum signal fraction analysis for artifact removal from biomedical signals | |
| CN113159205A (zh) | 一种基于最佳通道的稀疏时频块共空间模式特征提取方法 | |
| Ince et al. | ECoG based brain computer interface with subset selection | |
| Zhang et al. | Research on feature extraction algorithm commonly used in brain-computer interface technology |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |