CN111096745B - 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法 - Google Patents

基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法 Download PDF

Info

Publication number
CN111096745B
CN111096745B CN202010001431.5A CN202010001431A CN111096745B CN 111096745 B CN111096745 B CN 111096745B CN 202010001431 A CN202010001431 A CN 202010001431A CN 111096745 B CN111096745 B CN 111096745B
Authority
CN
China
Prior art keywords
steady
alpha
data
gamma
evoked response
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.)
Active
Application number
CN202010001431.5A
Other languages
English (en)
Other versions
CN111096745A (zh
Inventor
胡南
陈婷婷
曲铭雯
孙兵
王加俊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou University
Original Assignee
Suzhou University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Suzhou University filed Critical Suzhou University
Priority to CN202010001431.5A priority Critical patent/CN111096745B/zh
Publication of CN111096745A publication Critical patent/CN111096745A/zh
Priority to PCT/CN2020/105690 priority patent/WO2021135205A1/zh
Priority to US17/285,087 priority patent/US20220125385A1/en
Application granted granted Critical
Publication of CN111096745B publication Critical patent/CN111096745B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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]
    • A61B5/372—Analysis of electroencephalograms
    • A61B5/374—Detecting the frequency distribution of signals, e.g. detecting delta, theta, alpha, beta or gamma waves
    • 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/25—Bioelectric electrodes therefor
    • A61B5/279—Bioelectric electrodes therefor specially adapted for particular uses
    • A61B5/291—Bioelectric electrodes therefor specially adapted for particular uses for 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/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]
    • A61B5/377—Electroencephalography [EEG] using evoked responses
    • A61B5/38—Acoustic or auditory stimuli
    • 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/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • 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/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
    • 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
    • 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
    • A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier 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
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H30/00—ICT specially adapted for the handling or processing of medical images
    • G16H30/20—ICT specially adapted for the handling or processing of medical images for handling medical images, e.g. DICOM, HL7 or PACS
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G06N20/10—Machine learning using kernel methods, e.g. support vector machines [SVM]
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Artificial Intelligence (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Surgery (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Biophysics (AREA)
  • Veterinary Medicine (AREA)
  • Psychiatry (AREA)
  • Signal Processing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Psychology (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Fuzzy Systems (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Mathematical Optimization (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Epidemiology (AREA)

Abstract

本发明公开了一种基于稀疏贝叶斯学习的稳态诱发响应(SSR)脑源定位方法。本发明首先将稳态诱发响应记录分为多个数据段，通过快速傅里叶变换提取各段稳态诱发响应记录数据段的频域信息并构造数据矩阵。然后设置迭代自动停止条件以及稀疏支撑向量和自发脑电‑电噪声联合功率向量的初始值。而后迭代更新信号的后验均值与后验协方差并由此更新稀疏支撑向量和自发脑电‑电噪声联合功率向量。最后，当迭代结束时利用最新的稀疏支撑向量给出源定位结果。本发明在频域上对稳态诱发响应源定位问题进行建模，结合多段数据中信号的联合稀疏性，在稀疏贝叶斯学习框架下给出了适用于各种稳态诱发响应的脑源定位方法。

Description

基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
技术领域
本发明涉及稳态诱发响应领域，具体涉及一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)脑源定位方法。
背景技术
稳态诱发响应(SSR)是由周期性刺激信号诱发产生的脑电图(EEG)信号成分，具有与周期性刺激信号频率相同的正弦信号形式。与瞬态诱发电位相比，稳态诱发响应具有易于使用傅立叶分析在频域上分辨信号和噪声的优点。目前对稳态诱发响应的研究主要集中在听觉、视觉和体感等领域，在医学领域有重要的潜在应用价值。如听觉稳态响应由周期性的听觉刺激产生，可用于听力测试，麻醉监测和神经学评估，而视觉稳态诱发响应可能与被试的注意力运作方式有关。
从EEG信号起源的角度阐明认知过程的基本机制一直是亟待解决的问题。在头皮上通过电极采集到的脑电记录可以反映大脑内部的神经元活动，对脑电记录和神经元电生理活动关系的合理建模成为解决脑源定位问题的关键。经过几十年的研究，已经有很多关于脑源定位的研究。就源模型的选择而言，通常更倾向于分布式源模型。与等效电流偶极子模型相比，分布式源模型不对源的数量作假设，而且可以得到更好的定位精度。当使用分布式源模型时，利用观测到的EEG数据对激活神经元的空间分布进行估计是一个严重欠定的逆问题。为约束解空间，需要引入合理的先验假设。传统的线性分布式方法通常使用固定并且已知的先验假设，而人为主观设定的先验假设很大程度上影响了源定位结果的准确性。为了得到从数据角度考虑下更为合理的先验条件，越来越多的研究将贝叶斯方法引入脑源定位问题，在贝叶斯框架中以先验分布的形式嵌入先验假设，并通过贝叶斯推理过程确定合适的先验。
虽然稳态诱发响应在频域上的特征是稀疏的，可以通过傅里叶分析在频域上对稳态诱发响应源定位问题进行简化，然而目前还没有研究针对稳态诱发响应特征提出的源定位方法。基于稳态诱发响应在频域上的特征，在贝叶斯框架下建立一个适用于稳态诱发响应脑源定位方法将具有十分重要的意义。
目前尚无已公开专利应用解决基于稀疏贝叶斯学习方法的稳态诱发响应脑源定位问题，但是有一些论文提出了有关贝叶斯方法应用于瞬态诱发电位的脑源定位问题，如：
文献(Wipf and Nagarajan,A unified Bayesian framework for MEG/EEGsource imaging)分析并扩展了适用于源定位问题的几大类贝叶斯方法，包括经验贝叶斯方法，标准MAP估计和多元变分贝叶斯近似，结合现有方法进行改进，提出了可应用于脑电图/脑磁图源成像的统一贝叶斯框架。
文献(Saha,et al.Evaluating the Performance of BSBL Methodology forEEG Source Localization On a Realistic Head Model)和(Costa,et al.Bayesian EEGSource Localization Using a StructuredSparsity Prior)提出为了提高源定位算法的性能，还需要考虑结构稀疏性，并在贝叶斯框架中分别引入了内部块体结构信息和多元Bernoulli-Laplacian结构稀疏先验。
综合当前研究现状，贝叶斯方法应用于脑源定位问题具有很好的优势，但应用于稳态诱发响应脑源定位还需要提出切实可行的方法。
传统技术存在以下技术问题：
1、需要人为设定先验假设或经验参数，很大程度上限制并影响了源定位的结果。
2、未充分利用稳态诱发响应的频域信息进行源定位分析。
3、未充分利用数据段间的联合稀疏性。
4、低信噪比情况下，源定位结果空间分辨率不足。
发明内容
本发明要解决的技术问题是提供一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)脑源定位方法，(1)如何在频域上对稳态诱发响应脑源定位问题进行建模；(2)如何在不需要预设经验参数或变量的前提下，完全基于采集到的脑电记录数据对源的位置进行估计；(3)如何充分利用多段稳态诱发响应脑电数据之间的联合稀疏性，利用多段稳态诱发响应数据进行源定位；(4)如何在低信噪比条件下得到较为精确的稳态诱发响应脑源定位结果。
为了解决上述技术问题，本发明提供了一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，包括：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量
Figure BDA0002353647190000031
并将L个数据段对应的
Figure BDA0002353647190000032
整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中
Figure BDA0002353647190000033
其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：迭代阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵
Figure BDA0002353647190000034
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或
Figure BDA0002353647190000041
若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
在其中一个实施例中，对数据进行预处理以包括：数据的基线校正与叠加平均。
在其中一个实施例中，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为
Figure BDA0002353647190000042
其中
Figure BDA0002353647190000043
(·)Η表示共轭转置，⊙表示Hadamard积，而||·||F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
Figure BDA0002353647190000044
在其中一个实施例中，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure BDA0002353647190000045
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1)
Figure BDA0002353647190000046
其中Λ＝diag(α)；
(2)
Figure BDA0002353647190000047
在其中一个实施例中，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1)
Figure BDA0002353647190000048
其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2)
Figure BDA0002353647190000051
其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
在其中一个实施例中，头模型和电极分布对应的导程场矩阵
Figure BDA0002353647190000056
的计算包括：
(1)使用边界元方法(BEM)或有限元法(FEM)计算由头皮、头骨和皮质组成的3层真实头模型；
(2)在皮质表面上以N个均匀分布的三角网格表示分布式偶极子的空间位置；
(3)结合M导联的脑电采集电极在头皮上的空间信息，通过OpenMEEG软件计算出基于真实头模型的导程场矩阵
Figure BDA0002353647190000052
其规模为M×N，后续可以将该矩阵保存以重复利用。
一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，包括：计算机，所述计算机被编程以便执行如下步骤：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量
Figure BDA0002353647190000053
并将L个数据段对应的
Figure BDA0002353647190000054
整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中
Figure BDA0002353647190000055
其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵
Figure BDA0002353647190000061
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或
Figure BDA0002353647190000062
若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
在其中一个实施例中，对数据进行预处理以包括：数据的基线校正与叠加平均。
在其中一个实施例中，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为
Figure BDA0002353647190000063
其中
Figure BDA0002353647190000064
⊙表示Hadamard积，而||·｜｜F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
Figure BDA0002353647190000065
在其中一个实施例中，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵L，旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1)
Figure BDA0002353647190000066
其中Λ＝diag(α)；
(2)
Figure BDA0002353647190000067
在其中一个实施例中，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1)
Figure BDA0002353647190000071
其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2)
Figure BDA0002353647190000072
其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
在其中一个实施例中，头模型和电极分布对应的导程场矩阵
Figure BDA0002353647190000074
的计算包括：
(1)使用边界元方法(BEM)或有限元法(FEM)计算由头皮、头骨、脑脊液和皮质组成的4层真实头模型；
(2)在皮质表面上以N个均匀分布的三角网格表示分布式偶极子的空间位置；
(3)结合M导联的脑电采集电极在头皮上的空间信息，通过OpenMEEG软件计算出基于真实头模型的导程场矩阵
Figure BDA0002353647190000073
其规模为M×N，后续可以将该矩阵保存以重复利用。
本发明的有益效果：
(1)本发明根据稳态诱发响应的频域特性在频域上对稳态诱发响应源定位问题进行建模；
(2)本发明利用了多数据段间信号的联合稀疏性，得到的稳态诱发响应源定位结果具有更好的一般性；
(3)本发明中的稀疏支撑向量和自发脑电-电噪声功率的估计是在数据驱动的迭代过程中学习得到的，而不需要根据经验设置任何先验条件；
(4)本发明可以在低信噪比的条件下联合多段数据对稳态诱发响应进行源定位。
附图说明
图1是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法的流程图。
图2是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源分布的三种情况示意图。
图3是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为1时的定位效果对比图。
图4是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为2时的定位效果对比图。
图5是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为3时的定位效果对比图。
图6是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的SSR-SBL方法取得的40-Hz听觉稳态诱发响应源定位效果图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
本发明旨在根据稳态诱发响应(SSR)频域特性，建立一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应源定位算法。首先稳态诱发响应记录被分为多个数据段，通过快速傅里叶变换(FFT)，可以提取各段稳态诱发响应记录数据的频域信息。在稀疏贝叶斯学习框架下，在频域上对稳态诱发响应源定位问题进行建模。结合多段数据中稳态诱发响应源信号具有的联合稀疏性，利用稀疏支撑向量来约束解空间，最终通过期望最大化(EM)算法实现了由数据驱动学习稀疏支撑向量的迭代过程，从而得到了适用于稳态诱发响应脑源定位问题的稀疏贝叶斯学习算法。
可以理解，基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法和基于同样的发明构思的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置中所述计算机被编程以便执行的步骤一样，下面只详述基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法。
本发明的技术方案由图1所示的流程图示出：
一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法包括：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段。经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L。
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量
Figure BDA0002353647190000091
并将L个数据段对应的
Figure BDA0002353647190000092
整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中
Figure BDA0002353647190000093
其中(·)Τ表示转置。若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X。
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure BDA0002353647190000094
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或
Figure BDA0002353647190000095
若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
其中：
1、对数据进行预处理以包括：数据的基线校正与叠加平均。
2、对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为
Figure BDA0002353647190000101
其中
Figure BDA0002353647190000102
(·)Η表示共轭转置，⊙表示Hadamard积，而｜｜·｜｜F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
Figure BDA0002353647190000103
3、利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure BDA0002353647190000104
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1)
Figure BDA0002353647190000105
其中Λ＝diag(α)；
(2)
Figure BDA0002353647190000106
4、根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1)
Figure BDA0002353647190000107
其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2)
Figure BDA0002353647190000108
其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
附：本发明涉及稀疏贝叶斯学习算法的理论推导过程
(1)信号模型
稳态诱发响应(SSR)由周期性刺激信号诱发产生，具有正弦信号的形式，其频率与刺激信号的频率相同。假设某种稳态诱发响应由频率为f0的刺激信号诱发产生，在皮层上产生了K个正弦形式的激发电流，即s1(t),s2(t),...,sK(t)，且在一定时间内，这些信号相对稳定，具有如下形式：
Figure BDA0002353647190000111
其中αk和
Figure BDA0002353647190000112
分别表示sk(t)的幅度和相位。根据欧拉公式，sk(t)又可以分解为两个部分：
Figure BDA0002353647190000113
对于更长时间的刺激和脑电图记录，本发明可以假设振幅在短时间内基本保持恒定，但超过该时间片就可能发生变化。在此假设下，如果将整个记录分为多个片段，则在第l个数据片段中的第K个源信号可以表示为：
Figure BDA0002353647190000114
其中
Figure BDA0002353647190000115
和
Figure BDA0002353647190000116
分别表示对应于复分量
Figure BDA0002353647190000117
和
Figure BDA0002353647190000118
的权重。
头皮上采集到的脑电记录是神经元活动通过真实的头模型的传导过程产生的。经典的头模型通常包括具有不同电导率的三层，包括皮质，头骨和头皮。这种信号传输过程可以通过导前场矩阵(LFM)L0对偶极子电流活动的加权进行建模，其中LFM可通过有限元法(FEM)或边界元法(BEM)求出。如果在EEG采集中使用了M个电极，则当某种类型的周期性刺激诱发稳态诱发响应时，第l个数据段中的EEG信号xl(t)可以以向量形式表示为：
Figure BDA0002353647190000119
其中L0(rk)是位于rk位置的第k个源的LFM分量，ψk是包含3个元素的方向向量，用来表示第k个源的电流方向。ξl(t)和nl(t)分别表示第l个数据段中包含的自发脑电和电噪声。实际上自发脑电也可以由神经元活动通过LFM的映射来表示。但在定位稳态诱发响应相关的神经元活动的问题中，自发脑电被视为噪声干扰。
如果仅关注皮层表面的神经元活动，则由于锥体皮质神经元垂直于皮质表面，公式(4)表示的模型可以进一步简化为：
xl(t)＝Lsl(t)+ξl(t)+nl(t)， (5)
其中，源的方向向量ψ1,ψ2,…,ψK由ψ(r1),ψ(r2),…,ψ(rK)表示，LFM的第k列对应于第k个源的分量为
Figure BDA0002353647190000121
表示K个源产生的稳态诱发响应信号波形，符号(·)T表示转置。将(3)式代入(5)式得到：
Figure BDA0002353647190000122
其中
Figure BDA0002353647190000123
对采集到的EEG数据进行时频分析，比如对第l段数据进行快速傅里叶变换(FFT)，可以提取数据在刺激频率f0处的分量。根据(6)式，EEG记录在f0处的分量
Figure BDA0002353647190000124
可以表示为：
Figure BDA0002353647190000125
其中
Figure BDA0002353647190000126
和
Figure BDA0002353647190000127
分别表示ξl(t)和nl(t)在f0频点的分量。假设
Figure BDA0002353647190000128
在数据段之间可以是可变的或不变的，而
Figure BDA0002353647190000129
始终服从具有固定方差
Figure BDA00023536471900001210
的复数高斯分布，即
Figure BDA00023536471900001211
假定自发EEG分量
Figure BDA00023536471900001212
满足复数高斯分布，为了使模型更接近实际，认为方差
Figure BDA00023536471900001213
随数据段变化，即：
Figure BDA00023536471900001214
其中
Figure BDA00023536471900001215
服从伽马分布，即使得
Figure BDA00023536471900001216
服从K分布。
(2)层级贝叶斯模型
为解决公式(7)呈现的稳态诱发响应源定位问题，首先需要建立合适的贝叶斯网络。在(7)中，假设
Figure BDA00023536471900001217
和
Figure BDA00023536471900001218
分别满足方差为变量和方差固定不变的高斯分布。由于
Figure BDA00023536471900001219
和
Figure BDA00023536471900001220
相互独立，则有：
Figure BDA00023536471900001221
假设
Figure BDA00023536471900001222
且
Figure BDA00023536471900001223
给定，则
Figure BDA00023536471900001224
满足均值为
Figure BDA00023536471900001225
方差为γl的高斯分布，即
Figure BDA00023536471900001226
的条件分布可以表示为：
Figure BDA0002353647190000131
其中(·)Η表示共轭转置。联合L段数据对应的
Figure BDA0002353647190000132
可以得到联合多数据段结构信息的矩阵X。又因为当m≠n时，
Figure BDA0002353647190000133
与
Figure BDA0002353647190000134
无关，X的条件分布可以表示为：
Figure BDA0002353647190000135
这里
Figure BDA0002353647190000136
γ＝[γ1,γ2,…,γL]T。
由于稳态诱发响应激活源的位置目前还是未知的，因此应该在整个皮层范围内进行搜索。如果将整个皮层均匀划分为N个网格，每个网格对应的源的位置分别用
Figure BDA0002353647190000137
表示，则可以将扩展源空间后得到的大小为M×N的LFM表示为
Figure BDA0002353647190000138
通常情况下，在周期性刺激信号激发稳态诱发响应的过程中，位于皮质上的神经元只有少部分会被激活，即K＜N。所以本发明定义了N×L的行稀疏矩阵W，以表示分布在整个皮层上的N个源位置对应的活动强度，使得
Figure BDA0002353647190000139
当稳态诱发响应激发产生的某个源恰好位于第n个源位置，即
Figure BDA00023536471900001310
则第n个源位置对应于W的分量W[n,:]与第k个稳态诱发响应激活源对应的强度相对，即W[n,:]＝S[k,:]，否则，W[n,:]的元素都应该为0。只要
Figure BDA00023536471900001311
确定了，
Figure BDA00023536471900001312
就可以预先确定，而(10)中X的条件分布可以转化为：
Figure BDA00023536471900001313
其中wl是W的第l列。
由于对于某种确定类型的稳态诱发响应，激活源的位置在数据片段之间不会有很大变化，也就导致了W的行稀疏性。通过建模可以证明W列之间的联合稀疏性，即W的列服从相同的复数高斯分布，即：
Figure BDA0002353647190000141
其中N×1的向量α是公共方差向量,α中的数值较大的量表示稳态诱发响应激活源更有可能位于相应位置，当l≠l′，wl|α与wl′|α相互独立，故W的条件分布可以表示为：
Figure BDA0002353647190000142
其中wn,l表示W的第(n,l)个元素，可以发现，α控制了W的行稀疏性，因此，如果对α进行正确的估计，就可以给出稳态诱发响应激活源的位置。
在(11)和(13)中，γ和α被视为随机变量向量。在下面的分析中，它们将被替换为参数，即将(11)中的p(X|W,γ)和(13)中的p(W|α)分别替换为p(X|W；γ)和p(W；α)，参数γ和α将是贝叶斯推断的目标。
A.单个被试条件下的贝叶斯网络
对于单个被试，联合概率分布p(X,W；γ,α)可以由条件分布p(X|W；γ)和p(W；α)的乘积表示：
p(X,W；γ,α)＝p(X|W；γ)p(W；α)， (14)
其中W在贝叶斯推断中被视为隐变量矩阵。
B.多被试条件下的贝叶斯网络
为利用某种确定类型的稳态诱发响应源的共同位置，通常对多个被试进行稳态诱发响应引出实验。稀疏贝叶斯学习通过多个被试的数据组合来恢复常见源的分布。若共有P个被试参加了稳态诱发响应引出实验，则在f0频点提取的数据矩阵可以扩展并表示为
Figure BDA0002353647190000151
其中第p个被试的数据矩阵Xp由Lp个列向量组成。同样，本发明也定义了
Figure BDA0002353647190000152
和
Figure BDA0002353647190000153
多被试条件下稳态诱发响应激活源的联合位置信息可以由W的行稀疏性进行建模，即W1,W2,…WP共享的稀疏性。
Figure BDA0002353647190000154
其中α用于控制W1,W2,…WP之间相同的稀疏性。当p1≠p2，
Figure BDA0002353647190000155
与
Figure BDA0002353647190000156
相互独立，故
Figure BDA0002353647190000157
现在可以得到p(X,W；γ,α)＝p(X|W；γ)p(W；α)，与单个被试条件下得到的(14)相同。因此，相同的SBL方法可以应用于单被试场景和多被试场景。
(3)稀疏贝叶斯学习
基于上面建立的贝叶斯模型，可以使用期望最大化(EM)算法来推导参数γ和α。应用EM的前提是推导该模型中隐变量矩阵的后验分布，即p(W|X；γ,α)。
根据贝叶斯规则可以得到：
Figure BDA0002353647190000158
其中
Figure BDA0002353647190000159
Λ＝diag(α)。又因为当l≠l′时，wl|xl；γ,α与wl′|xl′；γ,α相互独立，可以得到：
Figure BDA00023536471900001510
在EM的过程中，通过最大化条件期望Ep(W|X；γ,α){log p(X,W；γ,α)}可以迭代地学习参数向量γ和α。又因为log p(X,W；γ,α)＝log p(X|W；γ)+log p(W；α)，而log p(X|W；γ)与α无关，所以在推导α的过程中，只需要最大化：
Q(α)＝Ep(W|X；γ,α){log p(W；α)}。 (18)
对Q(α)中的αn求偏导，可以得到α的更新值：
Figure BDA0002353647190000161
其中(·)*表示共轭运算，Σl[n,n]表示Σl的第(n,n)个元素，μl[n]为μl的第n个元素。类似地，由于log p(W；α)与γ无关，γ可由最大化Q(γ)的过程推导，其中
Q(γ)＝Ep(W|X；γ,α){log p(X|W；γ)}， (20)
通过Q(γ)中的γl求偏导，可以得到γ的更新值：
Figure BDA0002353647190000162
其中Re(·)和tr(·)分别表示取实部和求迹运算。
在本发明提出的SBL算法中，γ和α会不断地迭代更新。迭代的初始化为条件可以表示为：
Figure BDA0002353647190000163
其中
Figure BDA0002353647190000164
⊙表示Hadamard积，而||·||F表示Frobenius范数。每次迭代过程中计算的复杂度主要由计算Σl时的求逆运算产生。由于N＞M，Σl的更新值可以替代为：
Figure BDA0002353647190000165
基于稀疏贝叶斯学习的稳态诱发响应源定位算法性能验证仿真实验
A.计算头模型和导程场矩阵
在稳态诱发响应源定位实验中，本发明使用了解剖模板MNI/ICBM152以实现对大脑的无偏标准估计。通过分离MRI模板包含的白质和灰质间的界面，获得了皮层表面的空间信息。在Brainstorm软件的帮助下，本发明使用了边界元方法(BEM)计算由头皮、头骨和皮质组成的3层真实头模型。通过下采样，本发明得到了1500个均匀分布在皮质表面的三角网格以表示偶极子的空间位置，假设偶极子的电流方向垂直于皮质表面，结合63导联的脑电采集空间分布信息，通过OpenMEEG软件可以计算出基于真实头模型的导程场矩阵。
B.稳态诱发响应头皮记录仿真设计：
为合成此头模型下的稳态诱发响应头皮记录并测试稳态诱发响应源定位方法的性能，本发明需要根据已经得到的1500个皮层表面网格的空间位置设计稳态诱发响应激活源的分布，并将这些源的源活动设计为稳态诱发响应的信号形式。如图2所示，本发明共设计了3种源分布场景，分对应于1个源，2个源，3个源的情况。编号为1的源位于颞上回附近，编号为2的源位于额中回附近，编号为3的源位于颞叶上方的中央后回，仿真添加的源均位于左半脑。
指定了稳态诱发响应激活源的位置后，本发明设计激活源信号的形式为40Hz的正弦信号，即f0＝40。为模拟自发EEG干扰，本发明还在皮层上的其他区域添加了一定能量的噪声。位于皮层表面的源信号与导程场矩阵相乘后，就可以得到仿真合成相应的稳态诱发响应头皮记录。通过数值仿真，在每种稳态诱发响应激活源分布条件下，本发明都合成了10段(即，L＝10)长度为10秒的稳态诱发响应头皮记录，用于验证稳态诱发响应源定位算法的有效性。假设短时间内稳态诱发响应源信号和自发脑电的功率不变，本发明通过设置信噪比满足一个均值固定的高斯分布，实现同一数据段稳态诱发响应源信号的强度和自发脑电功率为一个确定的值，而不同数据段中，稳态诱发响应源信号的强度和脑电噪声的强度不同。
C.稳态诱发响应源定位结果：
为与其他SBL方法区分，本发明用SSR-SBL表示本发明中提出的基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)源定位算法。为证明SSR-SBL算法的性能，本发明将另外3种线性分布式源定位方法(包括MNE，sLORETA和LCMV波束形成)用于稳态诱发响应源定位以作比较。图3，图4和图5分别给出了当稳态诱发响应仿真源的个数分别为1、2和3时各种方法获得的定位结果。定位结果以灰度图的形式投射在皮层上，深色表示定位得到的稳态诱发响应源强度相对较大，浅灰色表示强度相对较小的稳态诱发响应源。本发明在这些图中仅显示了占定位结果总能量的前80％的结果。
如图3所示，所有种方法中，LCMV方法取得的定位结果最平滑。仅有一个激活源的情况下，相比于LCMV方法，MCE和sLORETA方法能够取得更好的定位精度。就MNE和sLORETA方法而言，当SNR＝10dB时，可以大致确定源的位置在颞叶，但是不能准确给出源的个数。当SNR＝0dB时，定位结果显示的源规模远大于实际尺寸，这时很难确定源的准确位置。当SNR＝-5dB时，这3种对比算法得到的定位结果向后半脑偏移，定位结果很不理想。对比下，本发明提出的SSR-SBL算法在相同条件下，可以实现无定位误差的定位效果。
除SSR-SBL算法外，随着仿真源数量的增加，相同信噪比条件下，其他3种源定位方法取得的定位效果越来越差。比如在稳态诱发响应仿真源的个数为2的情况下，当SNR＝0dB时，无论是MCE还是sLORETA方法取得的定位结果都无法分辨源的个数，更不用说源的尺寸。在稳态诱发响应仿真源的个数为2的情况下，当SNR＝10dB时，相比于MCE，sLORETA方法取得的定位结果更加平滑，很难区分其结果显示的稳态诱发响应源的激活区域为2个还是3个。
综合以上结果及分析，SSR-SBL算法区别于其他线性定位方法的主要特征是，SSR-SBL算法取得的定位结果具有很强的稀疏性，定位结果中每个分离的区域都可以解释为一个可能的源。目前的仿真结果可以证明在信噪比大于-5dB时，SSR-SBL源定位方法可以实现稳态诱发响应激活源位置的无偏估计。
SSR-SBL源定位方法应用于听觉稳态诱发响应源定位
作为一种典型的稳态诱发响应，听觉稳态诱发响应(ASSR)有广泛的临床应用前景，例如听觉稳态诱发响应在听力检测，意识评估和麻醉监测领域都有应用。探究听觉稳态诱发响应起源对解释其临床应用具有重要意义。目前很多研究认为听觉稳态诱发响应主要起源于初级听觉皮层和脑干，但尚无确切的结论。通过本发明中提出的SSR-SBL源定位算法，可以探究听觉稳态诱发响应的皮质起源。这里我们首先将SSR-SBL算法应用于40-Hz听觉稳态诱发响应源定位实验。
本发明收集了9名听力正常的健康志愿者，包括5例男性和4例女性，年龄21±1.89岁，无听力损失和神经系统疾病史。在获得所有受试者书面同意的前提下，对被试实施39.1Hz的正弦振幅调制音听力刺激，用于诱发40-Hz听觉稳态诱发响应，即f0＝39.1。具体来说，诱发40-Hz听觉稳态诱发响应的实验条件为：被试双耳均接受载波频率为500Hz的39.1Hz调幅音作为听力刺激，调制程度为100％，声音刺激的强度为70dB。由于在仿真实验中算法SSR-SBL源定位算法的有效性已经得到验证，因此本发明认为SSR-SBL源定位算法得到的40-Hz听觉稳态诱发响应源定位结果相对可靠。本发明在40-Hz听觉稳态诱发响应引出实验中采用了标准的64通道电极分布(除地电极外，有效电极数为63，即M＝63)和BrainVision Recorder软件记录听觉稳态诱发响应，记录时的采样频率为1000Hz，同时使用了通带频率为1-300Hz(12dB/倍频程)的带通滤波器和50Hz陷波滤波器来抑制噪声干扰。每个被试的脑电采集时长为200s，首先将200s的脑电数据等分为10段长20s的相邻但不重叠的数据段，即L＝10。在40-Hz听觉稳态诱发响应源定位实验之前，可以通过预处理，提高数据的信噪比。
预处理过程包括对每个被试对应的10段EEG数据进行去基线校正、分帧、去除伪迹、滤波和时域平均：第一步，将EEG记录的前5秒数据用于基线校正；第二步，将基线校正后的EEG记录中每1023个数据点分成一个数据帧，每个数据帧的长度恰好是声音刺激信号的40个整周期，且数据帧之间相邻但不重叠。第三步，移除包含数据点幅度超过40uV的数据帧，以去除伪迹。第四步，将数据帧通过中心频率为40Hz，带宽为10Hz的带通滤波器进行滤波。最后，保留10个数据段中前10个数据帧。每个数据段进行5次时域平均得到待处理数据段xl(t),l＝1,2,…L。预处理后，通过对xl(t)进行FFT，得到数据段在f0频点对应的分量
Figure BDA0002353647190000201
其中
Figure BDA0002353647190000202
向量的大小为63×1。在针对单个被试的SSR-SBL模型下，将单个被试得到的10段数据对应的
Figure BDA0002353647190000203
扩展成大小为63×10的输入矩阵X，其中
Figure BDA0002353647190000204
在针对多被试的SSR-SBL模型下，将9个被试对应的63×10的输入矩阵Xp,p＝1,2,…P结果进行再扩展，得到大小为63×90的多被试联合输入矩阵X，其中
Figure BDA0002353647190000205
在已知L＝10,P＝9且导程场矩阵L已确定的情况下，可以通过SSR-SBL算法给出的迭代过程学习得到40-Hz听觉稳态诱发响应的源分布结果，本发明在迭代过程中设置迭代误差阈值ε＝10-3，最大迭代次数Niter＝1000。
SSR-SBL算法取得的定位结果如图6所示。多被试联合定位结果显示，双耳刺激下得到的40-Hz听觉稳态诱发响应的激活神经元主要分布在左右脑的颞叶，在前额也有一定的分布。相比于单个被试的定位结果，多被试联合定位结果更加稳定，可以假设它是单个被试的定位结果的综合，避免了在单个被试的定位结果中出现的噪声过大、激活源能量估计过小的情况。
以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (4)

1.一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，包括：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量
Figure FDA0002941347770000011
并将L个数据段对应的
Figure FDA0002941347770000012
整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中
Figure FDA0002941347770000013
(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure FDA0002941347770000018
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或
Figure FDA0002941347770000014
若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果；
利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure FDA0002941347770000015
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1)
Figure FDA0002941347770000016
其中Λ＝diag(α)；
(2)
Figure FDA0002941347770000017
根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1)
Figure FDA0002941347770000021
其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2)
Figure FDA0002941347770000022
其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算；
其中，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为
Figure FDA0002941347770000023
其中
Figure FDA0002941347770000024
(·)Η表示共轭转置，⊙表示Hadamard积，而||·||F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
Figure FDA0002941347770000025
2.如权利要求1所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，对数据进行预处理包括：数据的基线校正与叠加平均。
3.一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，包括：计算机，所述计算机被编程以便执行如下步骤：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量
Figure FDA0002941347770000026
并将L个数据段对应的
Figure FDA0002941347770000027
整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中
Figure FDA0002941347770000028
其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure FDA0002941347770000031
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或
Figure FDA0002941347770000032
若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果；
利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵
Figure FDA0002941347770000033
旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1)
Figure FDA0002941347770000034
其中Λ＝diag(α)；
(2)
Figure FDA0002941347770000035
根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1)
Figure FDA0002941347770000036
其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2)
Figure FDA0002941347770000037
其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算；
其中，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为
Figure FDA0002941347770000041
其中
Figure FDA0002941347770000042
(·)Η表示共轭转置，⊙表示Hadamard积，而||·||F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
Figure FDA0002941347770000043
4.根据权利要求3所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，对数据进行预处理包括：数据的基线校正与叠加平均。
CN202010001431.5A 2020-01-02 2020-01-02 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法 Active CN111096745B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN202010001431.5A CN111096745B (zh) 2020-01-02 2020-01-02 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
PCT/CN2020/105690 WO2021135205A1 (zh) 2020-01-02 2020-07-30 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
US17/285,087 US20220125385A1 (en) 2020-01-02 2020-07-30 Sbl-based ssr brain source localization method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010001431.5A CN111096745B (zh) 2020-01-02 2020-01-02 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法

Publications (2)

Publication Number Publication Date
CN111096745A CN111096745A (zh) 2020-05-05
CN111096745B true CN111096745B (zh) 2021-03-30

Family

ID=70426947

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010001431.5A Active CN111096745B (zh) 2020-01-02 2020-01-02 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法

Country Status (3)

Country Link
US (1) US20220125385A1 (zh)
CN (1) CN111096745B (zh)
WO (1) WO2021135205A1 (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111096745B (zh) * 2020-01-02 2021-03-30 苏州大学 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
CN112401905B (zh) * 2020-11-11 2021-07-30 东南大学 一种基于源定位和脑网络的自然动作脑电识别方法
CN112674773B (zh) * 2020-12-22 2021-12-24 北京航空航天大学 基于Tucker分解和ripple时间窗的脑磁图源定位方法和装置
CN114052668B (zh) * 2022-01-17 2022-06-17 北京航空航天大学杭州创新研究院 一种基于脑磁图数据的脑功能分析方法
CN116491960B (zh) * 2023-06-28 2023-09-19 南昌大学第一附属医院 脑瞬态监测设备、电子设备及存储介质

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140323899A1 (en) * 2006-12-22 2014-10-30 Neuro-Insight Pty. Ltd. Psychological Evaluation and Methods of Use
CN104688222B (zh) * 2015-04-03 2017-04-26 清华大学 基于脑电信号的音色合成装置
CN107703477B (zh) * 2017-09-11 2021-02-05 电子科技大学 基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法
WO2019182915A1 (en) * 2018-03-19 2019-09-26 The Board Of Trustees Of The Leland Stanford Junior University Treatment of depression
CN108802683B (zh) * 2018-05-30 2021-04-27 东南大学 一种基于稀疏贝叶斯学习的源定位方法
CN110251083A (zh) * 2019-06-20 2019-09-20 深圳大学 一种癫痫病灶的定位数据的处理方法、系统和存储介质
CN111096745B (zh) * 2020-01-02 2021-03-30 苏州大学 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法

Also Published As

Publication number Publication date
US20220125385A1 (en) 2022-04-28
CN111096745A (zh) 2020-05-05
WO2021135205A1 (zh) 2021-07-08

Similar Documents

Publication Publication Date Title
CN111096745B (zh) 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
Colclough et al. A symmetric multivariate leakage correction for MEG connectomes
Safieddine et al. Removal of muscle artifact from EEG data: comparison between stochastic (ICA and CCA) and deterministic (EMD and wavelet-based) approaches
Castaño-Candamil et al. Solving the EEG inverse problem based on space–time–frequency structured sparsity constraints
Dähne et al. SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters
Anemüller et al. Complex independent component analysis of frequency-domain electroencephalographic data
US12016634B2 (en) Methods and apparatus for electromagnetic source imaging using deep neural networks
CN110090017B (zh) 一种基于lstm的脑电信号源定位方法
Ramírez et al. Spectral signal space projection algorithm for frequency domain MEG and EEG denoising, whitening, and source imaging
CN109199376B (zh) 基于oa-wmne脑源成像的运动想象脑电信号的解码方法
US20160051161A1 (en) Method for locating a brain activity associated with a task
Hansen et al. Unmixing oscillatory brain activity by EEG source localization and empirical mode decomposition
Das et al. Neuro-current response functions: A unified approach to MEG source analysis under the continuous stimuli paradigm
Limpiti et al. A spatiotemporal framework for estimating trial-to-trial amplitude variation in event-related MEG/EEG
CN113995422B (zh) 基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法及系统
Ferdowsi et al. Multi layer spectral decomposition technique for ERD estimation in EEG μ rhythms: An EEG–fMRI study
Chen et al. Solving the inverse problem for auditory steady-state response via sparse Bayesian learning
Pezard et al. Why bother to spatially embed EEG? Comments on Pritchard et al., Psychophysiology, 33, 362–368, 1996
Maki et al. Enhancing event-related potentials based on maximum a posteriori estimation with a spatial correlation prior
Li et al. ICA on sensor or source data: a comparison study in deriving resting state networks from EEG
Woon et al. Can we learn anything from single-channel unaveraged MEG data?
Cabrera Implementation of Machine Learning Algorithm to Exploit Information from Multimodal fMRI/EEG Fused Image Data
Aydın Extraction of auditory evoked potentials from ongoing EEG
Limpiti et al. A spatiotemporal framework for MEG/EEG evoked response amplitude and latency variability estimation
Nuanprasert et al. Fst-Filter: A flexible spatio-temporal filter for biomedical multichannel data denoising

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
GR01 Patent grant
GR01 Patent grant