基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
阅读：28发布：2020-05-08
专利汇可以提供基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法专利检索，专利查询，专利分析的服务。并且本 发明 公开了一种基于稀疏贝叶斯学习的稳态诱发响应(SSR)脑源 定位 方法。本发明首先将稳态诱发响应记录分为多个数据段，通过快速傅里叶变换提取各段稳态诱发响应记录数据段的频域信息并构造数据矩阵。然后设置 迭代 自动停止条件以及稀疏 支撑 向量和自发脑电-电噪声联合功率向量的初始值。而后迭代更新 信号 的后验均值与后验协方差并由此更新稀疏支撑向量和自发脑电-电噪声联合功率向量。最后，当迭代结束时利用最新的稀疏支撑向量给出源定位结果。本发明在频域上对稳态诱发响应源定位问题进行建模，结合多段数据中信号的联合稀疏性，在稀疏贝叶斯学习 框架 下给出了适用于各种稳态诱发响应的脑源定位方法。，下面是基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法专利的具体信息内容。
1.一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，包括：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L。
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量 并将L个数据段对应的 整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中 (·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横
向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或 若不满足则令
α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
2.如权利要求1所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，对数据进行预处理以包括：数据的基线校正与叠加平均。
3.如权利要求1所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为 其中 (·)Η表示
共轭转置，⊙表示Hadamard积，而||·||F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
l＝1,2,…,L。
4.如权利要求1所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1) 其中Λ＝diag(α)；
(2)
5.如权利要求1所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，其特征在于，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1) n＝1,2,…,N，其中αnew[n]为αnew的第n个元
*
素，(·) 表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2) l＝1,2,…,L，其中γnew[l]为
γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
6.一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，包括：计算机，所述计算机被编程以便执行如下步骤：
步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量 并将L个数据段对应的 整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中 其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩
阵横向排列整合成一个矩阵X；
步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或 若不满足则令
α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
7.根据权利要求6所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，对数据进行预处理以包括：数据的基线校正与叠加平均。
8.根据权利要求6所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，对迭代程序中的各变量的初始化包括：
(1)稀疏支撑向量α初始化为 其中 ⊙表示
Hadamard积，而||·||F表示Frobenius范数；
(2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
l＝1,2,…,L。
9.根据权利要求6所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
(1) 其中Λ＝diag(α)；
(2)
10.根据权利要求6所述的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，其特征在于，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
(1) n＝1,2,…,N，其中αnew[n]为αnew的第n个元
素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
(2) l＝1,2,…,L，其中γnew[l]为
γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
技术领域
[0001] 本发明涉及稳态诱发响应领域，具体涉及一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)脑源定位方法。
背景技术
[0002] 稳态诱发响应(SSR)是由周期性刺激信号诱发产生的脑电图(EEG)信号成分，具有与周期性刺激信号频率相同的正弦信号形式。与瞬态诱发电位相比，稳态诱发响应具有易于使用傅立叶分析在频域上分辨信号和噪声的优点。目前对稳态诱发响应的研究主要集中在听觉、视觉和体感等领域，在医学领域有重要的潜在应用价值。如听觉稳态响应由周期性的听觉刺激产生，可用于听力测试，麻醉监测和神经学评估，而视觉稳态诱发响应可能与被试的注意力运作方式有关。
[0003] 从EEG信号起源的角度阐明认知过程的基本机制一直是亟待解决的问题。在头皮上通过电极采集到的脑电记录可以反映大脑内部的神经元活动，对脑电记录和神经元电生理活动关系的合理建模成为解决脑源定位问题的关键。经过几十年的研究，已经有很多关于脑源定位的研究。就源模型的选择而言，通常更倾向于分布式源模型。与等效电流偶极子模型相比，分布式源模型不对源的数量作假设，而且可以得到更好的定位精度。当使用分布式源模型时，利用观测到的EEG数据对激活神经元的空间分布进行估计是一个严重欠定的逆问题。为约束解空间，需要引入合理的先验假设。传统的线性分布式方法通常使用固定并且已知的先验假设，而人为主观设定的先验假设很大程度上影响了源定位结果的准确性。为了得到从数据角度考虑下更为合理的先验条件，越来越多的研究将贝叶斯方法引入脑源定位问题，在贝叶斯框架中以先验分布的形式嵌入先验假设，并通过贝叶斯推理过程确定合适的先验。
[0004] 虽然稳态诱发响应在频域上的特征是稀疏的，可以通过傅里叶分析在频域上对稳态诱发响应源定位问题进行简化，然而目前还没有研究针对稳态诱发响应特征提出的源定位方法。基于稳态诱发响应在频域上的特征，在贝叶斯框架下建立一个适用于稳态诱发响应脑源定位方法将具有十分重要的意义。
[0006] 文献(Wipf and Nagarajan,A unified Bayesian framework for MEG/EEG source imaging)分析并扩展了适用于源定位问题的几大类贝叶斯方法，包括经验贝叶斯方法，标准MAP估计和多元变分贝叶斯近似，结合现有方法进行改进，提出了可应用于脑电图/脑磁图源成像的统一贝叶斯框架。
[0007] 文献(Saha,et al.Evaluating the Performance of BSBL Methodology for EEG Source Localization On a Realistic Head Model)和(Costa,et al.Bayesian EEG Source Localization Using a StructuredSparsity Prior)提出为了提高源定位算法的性能，还需要考虑结构稀疏性，并在贝叶斯框架中分别引入了内部块体结构信息和多元Bernoulli-Laplacian结构稀疏先验。
[0008] 综合当前研究现状，贝叶斯方法应用于脑源定位问题具有很好的优势，但应用于稳态诱发响应脑源定位还需要提出切实可行的方法。
[0009] 传统技术存在以下技术问题：
[0010] 1、需要人为设定先验假设或经验参数，很大程度上限制并影响了源定位的结果。
[0011] 2、未充分利用稳态诱发响应的频域信息进行源定位分析。
[0012] 3、未充分利用数据段间的联合稀疏性。
发明内容
[0014] 本发明要解决的技术问题是提供一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)脑源定位方法，(1)如何在频域上对稳态诱发响应脑源定位问题进行建模；(2)如何在不需要预设经验参数或变量的前提下，完全基于采集到的脑电记录数据对源的位置进行估计；(3)如何充分利用多段稳态诱发响应脑电数据之间的联合稀疏性，利用多段稳态诱发响应数据进行源定位；(4)如何在低信噪比条件下得到较为精确的稳态诱发响应脑源定位结果。
[0015] 为了解决上述技术问题，本发明提供了一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，包括：
[0016] 步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
[0017] 步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量 并将L个数据段对应的 整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中 其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
[0019] 步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
[0020] 步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
[0021] 步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或 若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
[0023] 在其中一个实施例中，对迭代程序中的各变量的初始化包括：
[0024] (1)稀疏支撑向量α初始化为 其中 (·)Η表示共轭转置，⊙表示Hadamard积，而||·||F表示Frobenius范数；
[0025] (2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
[0026] 在其中一个实施例中，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
[0027] (1) 其中Λ＝diag(α)；
[0028] (2)
[0029] 在其中一个实施例中，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
[0030] (1) 其中αnew[n]为αnew的第n个元素，(·)*表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
[0031] (2) 其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
[0032] 在其中一个实施例中，头模型和电极分布对应的导程场矩阵 的计算包括：
[0033] (1)使用边界元方法(BEM)或有限元法(FEM)计算由头皮、头骨和皮质组成的3层真实头模型；
[0034] (2)在皮质表面上以N个均匀分布的三角网格表示分布式偶极子的空间位置；
[0036] 一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置，包括：计算机，所述计算机被编程以便执行如下步骤：
[0037] 步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段；经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L；
[0038] 步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量 并将L个数据段对应的 整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中 其中(·)Τ表示转置；若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X；
[0039] 步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
[0040] 步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
[0041] 步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
[0042] 步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或 若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
[0043] 在其中一个实施例中，对数据进行预处理以包括：数据的基线校正与叠加平均。
[0044] 在其中一个实施例中，对迭代程序中的各变量的初始化包括：
[0045] (1)稀疏支撑向量α初始化为 其中 ⊙表示Hadamard积，而||·｜｜F表示Frobenius范数；
[0046] (2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
[0047] 在其中一个实施例中，利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场导程矩阵L，旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
[0048] (1) 其中Λ＝diag(α)；
[0049] (2)
[0050] 在其中一个实施例中，根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
[0051] (1) 其中αnew[n]为αnew的第*
n个元素，(·) 表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
n个元素，(·) 表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
[0052] (2) 其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
[0053] 在其中一个实施例中，头模型和电极分布对应的导程场矩阵 的计算包括：
[0054] (1)使用边界元方法(BEM)或有限元法(FEM)计算由头皮、头骨、脑脊液和皮质组成的4层真实头模型；
[0055] (2)在皮质表面上以N个均匀分布的三角网格表示分布式偶极子的空间位置；
[0056] (3)结合M导联的脑电采集电极在头皮上的空间信息，通过OpenMEEG软件计算出基于真实头模型的导程场矩阵 其规模为M×N，后续可以将该矩阵保存以重复利用。
[0057] 本发明的有益效果：
[0058] (1)本发明根据稳态诱发响应的频域特性在频域上对稳态诱发响应源定位问题进行建模；
[0059] (2)本发明利用了多数据段间信号的联合稀疏性，得到的稳态诱发响应源定位结果具有更好的一般性；
[0060] (3)本发明中的稀疏支撑向量和自发脑电-电噪声功率的估计是在数据驱动的迭代过程中学习得到的，而不需要根据经验设置任何先验条件；
[0062] 图1是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法的流程图。
[0063] 图2是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源分布的三种情况示意图。
[0064] 图3是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为1时的定位效果对比图。
[0065] 图4是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为2时的定位效果对比图。
[0066] 图5是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的仿真稳态诱发响应源个数为3时的定位效果对比图。
[0067] 图6是本发明基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法中的SSR-SBL方法取得的40-Hz听觉稳态诱发响应源定位效果图。
具体实施方式
[0068] 下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
[0069] 本发明旨在根据稳态诱发响应(SSR)频域特性，建立一种基于稀疏贝叶斯学习(SBL)的稳态诱发响应源定位算法。首先稳态诱发响应记录被分为多个数据段，通过快速傅里叶变换(FFT)，可以提取各段稳态诱发响应记录数据的频域信息。在稀疏贝叶斯学习框架下，在频域上对稳态诱发响应源定位问题进行建模。结合多段数据中稳态诱发响应源信号具有的联合稀疏性，利用稀疏支撑向量来约束解空间，最终通过期望最大化(EM)算法实现了由数据驱动学习稀疏支撑向量的迭代过程，从而得到了适用于稳态诱发响应脑源定位问题的稀疏贝叶斯学习算法。
[0070] 可以理解，基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法和基于同样的发明构思的基于稀疏贝叶斯学习的稳态诱发响应脑源定位装置中所述计算机被编程以便执行的步骤一样，下面只详述基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法。
[0071] 本发明的技术方案由图1所示的流程图示出：
[0072] 一种基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法包括：
[0073] 步骤1.将M个电极采集到稳态诱发响应的头皮记录分为L段。经过对L段数据的预处理以提高信噪比，得到L个数据段xl(t),l＝1,2,…L。
[0074] 步骤2.对各数据段对应的xl(t)进行快速傅里叶变换，即FFT，提取xl(t)在刺激频点f0处的复数分量 并将L个数据段对应的 整合为联合稳态诱发响应多数据段结构信息的矩阵X，其中 其中(·)Τ表示转置。若有多个被试的数据，将他们的数据矩阵横向排列整合成一个矩阵X。
[0075] 步骤3.设置迭代程序的参数：误差阈值ε和最大迭代次数Niter，并对迭代程序中的各变量初始化：稀疏支撑向量α初始化为αinit，自发脑电-电噪声联合功率向量γ初始化为γinit；
[0076] 步骤4.利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算第l段数据对应的源信号的后验均值向量μl和后验协方差矩阵Σl；
[0077] 步骤5.根据μl、Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew；
[0078] 步骤6.判断是否满足迭代停止条件：迭代次数n≥Niter或 若不满足则令α＝αnew，γ＝γnew，回到步骤4继续执行迭代；否则结束迭代，输出稀疏支撑向量α，得到源定位结果。
[0079] 其中：
[0080] 1、对数据进行预处理以包括：数据的基线校正与叠加平均。
[0081] 2、对迭代程序中的各变量的初始化包括：
[0082] (1)稀疏支撑向量α初始化为 其中 (·)Η表示共轭转置，⊙表示Hadamard积，而｜｜·｜｜F表示Frobenius范数；
[0083] (2)第l段数据对应的自发脑电-电噪声联合功率γl初始化为
[0084] 3、利用稳态诱发响应数据信息矩阵X、头模型和电极分布对应的导程场矩阵 旧的稀疏支撑向量α和自发脑电-电噪声联合功率向量γ计算μl和Σl包括：
[0085] (1) 其中Λ＝diag(α)；
[0086] (2)
[0087] 4、根据μl和Σl计算稀疏支撑向量α和自发脑电-电噪声联合功率向量γ的更新值αnew和γnew包括：
[0088] (1) 其中αnew[n]为αnew的第*
n个元素，(·) 表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
n个元素，(·) 表示共轭运算，Σl[n,n]为Σl的第(n,n)个元素，μl[n]为μl的第n个元素；
[0089] (2) 其中γnew[l]为γnew的第l个元素，Re(·)和tr(·)分别为取实部和求迹运算。
[0090] 附：本发明涉及稀疏贝叶斯学习算法的理论推导过程
[0091] (1)信号模型
[0092] 稳态诱发响应(SSR)由周期性刺激信号诱发产生，具有正弦信号的形式，其频率与刺激信号的频率相同。假设某种稳态诱发响应由频率为f0的刺激信号诱发产生，在皮层上产生了K个正弦形式的激发电流，即s1(t),s2(t),...,sK(t)，且在一定时间内，这些信号相对稳定，具有如下形式：
[0093]
[0094] 其中αk和 分别表示sk(t)的幅度和相位。根据欧拉公式，sk(t)又可以分解为两个部分：
[0095]
[0096] 对于更长时间的刺激和脑电图记录，本发明可以假设振幅在短时间内基本保持恒定，但超过该时间片就可能发生变化。在此假设下，如果将整个记录分为多个片段，则在第l个数据片段中的第K个源信号可以表示为：
[0097]
[0098] 其中 和 分别表示对应于复分量 和 的权重。
[0099] 头皮上采集到的脑电记录是神经元活动通过真实的头模型的传导过程产生的。经典的头模型通常包括具有不同电导率的三层，包括皮质，头骨和头皮。这种信号传输过程可以通过导前场矩阵(LFM)L0对偶极子电流活动的加权进行建模，其中LFM可通过有限元法(FEM)或边界元法(BEM)求出。如果在EEG采集中使用了M个电极，则当某种类型的周期性刺激诱发稳态诱发响应时，第l个数据段中的EEG信号xl(t)可以以向量形式表示为：
[0100]
[0101] 其中L0(rk)是位于rk位置的第k个源的LFM分量，ψk是包含3个元素的方向向量，用来表示第k个源的电流方向。ξl(t)和nl(t)分别表示第l个数据段中包含的自发脑电和电噪声。实际上自发脑电也可以由神经元活动通过LFM的映射来表示。但在定位稳态诱发响应相关的神经元活动的问题中，自发脑电被视为噪声干扰。
[0102] 如果仅关注皮层表面的神经元活动，则由于锥体皮质神经元垂直于皮质表面，公式(4)表示的模型可以进一步简化为：
[0103] xl(t)＝Lsl(t)+ξl(t)+nl(t)， (5)
[0104] 其中，源的方向向量ψ1,ψ2,…,ψK由ψ(r1),ψ(r2),…,ψ(rK)表示，LFM的第k列对应于第k个源的分量为 表示KT
个源产生的稳态诱发响应信号波形，符号(·) 表示转置。将(3)式代入(5)式得到：
个源产生的稳态诱发响应信号波形，符号(·) 表示转置。将(3)式代入(5)式得到：
[0105]
[0106] 其中
[0107] 对采集到的EEG数据进行时频分析，比如对第l段数据进行快速傅里叶变换(FFT)，可以提取数据在刺激频率f0处的分量。根据(6)式，EEG记录在f0处的分量 可以表示为：
[0108]
[0109] 其中 和 分别表示ξl(t)和nl(t)在f0频点的分量。假设 在数据段之间可以是可变的或不变的，而 始终服从具有固定方差 的复数高斯分布，即 假定自发EEG分量 满足复数高斯分布，为了使模型更接近实际，认为方差 随数据段变化，即：
[0110]
[0112] (2)层级贝叶斯模型
[0113] 为解决公式(7)呈现的稳态诱发响应源定位问题，首先需要建立合适的贝叶斯网络。在(7)中，假设 和 分别满足方差为变量和方差固定不变的高斯分布。由于 和 相互独立，则有： 假设 且 给定，则 满足均值为方差为γl的高斯分布，即 的条件分布可以表示为：
[0114]
[0115] 其中(·)Η表示共轭转置。联合L段数据对应的 可以得到联合多数据段结构信息的矩阵X。又因为当m≠n时， 与 无关，X的条件分布可以表示
为：
为：
[0116]
[0117] 这里 γ＝[γ1,γ2,…,γL]T。
[0118] 由于稳态诱发响应激活源的位置目前还是未知的，因此应该在整个皮层范围内进行搜索。如果将整个皮层均匀划分为N个网格，每个网格对应的源的位置分别用 表示，则可以将扩展源空间后得到的大小为M×N的LFM表示为通常情况下，在周期性刺激信号激发稳态诱发响应的过程中，位于皮质上的神经元只有少部分会被激活，即K＜N。所以本发明定义了N×L的行稀疏矩阵W，以表示分布在整个皮层上的N个源位置对应的活动强度，使得 当稳态诱发响应激发产生的某个源恰好位
于第n个源位置，即 则第n个源位置对应于W的分量W[n,:]与第k
个稳态诱发响应激活源对应的强度相对，即W[n,:]＝S[k,:]，否则，W[n,:]的元素都应该为
0。只要 确定了， 就可以预先确定，而(10)中X的条件分布可以转化为：
于第n个源位置，即 则第n个源位置对应于W的分量W[n,:]与第k
个稳态诱发响应激活源对应的强度相对，即W[n,:]＝S[k,:]，否则，W[n,:]的元素都应该为
0。只要 确定了， 就可以预先确定，而(10)中X的条件分布可以转化为：
[0119]
[0120] 其中wl是W的第l列。
[0121] 由于对于某种确定类型的稳态诱发响应，激活源的位置在数据片段之间不会有很大变化，也就导致了W的行稀疏性。通过建模可以证明W列之间的联合稀疏性，即W的列服从相同的复数高斯分布，即：
[0122]
[0123] 其中N×1的向量α是公共方差向量,α中的数值较大的量表示稳态诱发响应激活源更有可能位于相应位置，当l≠l′，wl|α与wl′|α相互独立，故W的条件分布可以表示为：
[0124]
[0125] 其中wn,l表示W的第(n,l)个元素，可以发现，α控制了W的行稀疏性，因此，如果对α进行正确的估计，就可以给出稳态诱发响应激活源的位置。
[0126] 在(11)和(13)中，γ和α被视为随机变量向量。在下面的分析中，它们将被替换为参数，即将(11)中的p(X|W,γ)和(13)中的p(W|α)分别替换为p(X|W；γ)和p(W；α)，参数γ和α将是贝叶斯推断的目标。
[0127] A.单个被试条件下的贝叶斯网络
[0128] 对于单个被试，联合概率分布p(X,W；γ,α)可以由条件分布p(X|W；γ)和p(W；α)的乘积表示：
[0129] p(X,W；γ,α)＝p(X|W；γ)p(W；α)， (14)
[0130] 其中W在贝叶斯推断中被视为隐变量矩阵。
[0131] B.多被试条件下的贝叶斯网络
[0132] 为利用某种确定类型的稳态诱发响应源的共同位置，通常对多个被试进行稳态诱发响应引出实验。稀疏贝叶斯学习通过多个被试的数据组合来恢复常见源的分布。若共有P个被试参加了稳态诱发响应引出实验，则在f0频点提取的数据矩阵可以扩展并表示为其中第p个被试的数据矩阵Xp由Lp个列向量组成。同样，本发明也定义了 和 多被试条件下稳态诱发响应激活源的联合
位置信息可以由W的行稀疏性进行建模，即W1,W2,…WP共享的稀疏性。
位置信息可以由W的行稀疏性进行建模，即W1,W2,…WP共享的稀疏性。
[0133]
[0134] 其中α用于控制W1,W2,…WP之间相同的稀疏性。当p1≠p2， 与相互独立，故 现在可以得到p(X,W；γ,α)＝p(X|W；γ)p
(W；α)，与单个被试条件下得到的(14)相同。因此，相同的SBL方法可以应用于单被试场景和多被试场景。
(W；α)，与单个被试条件下得到的(14)相同。因此，相同的SBL方法可以应用于单被试场景和多被试场景。
[0135] (3)稀疏贝叶斯学习
[0136] 基于上面建立的贝叶斯模型，可以使用期望最大化(EM)算法来推导参数γ和α。应用EM的前提是推导该模型中隐变量矩阵的后验分布，即p(W|X；γ,α)。
[0137] 根据贝叶斯规则可以得到：
[0138]
[0139] 其中 Λ＝diag(α)。又因为当l≠l′时，wl|xl；γ,α与wl′|xl′；γ,α相互独立，可以得到：
[0140]
[0141] 在EM的过程中，通过最大化条件期望Ep(W|X；γ,α){log p(X,W；γ,α)}可以迭代地学习参数向量γ和α。又因为log p(X,W；γ,α)＝log p(X|W；γ)+log p(W；α)，而log p(X|W；γ)与α无关，所以在推导α的过程中，只需要最大化：
[0142] Q(α)＝Ep(W|X；γ,α){log p(W；α)}。 (18)
[0143] 对Q(α)中的αn求偏导，可以得到α的更新值：
[0144]
[0145] 其中(·)*表示共轭运算，Σl[n,n]表示Σl的第(n,n)个元素，μl[n]为μl的第n个元素。类似地，由于log p(W；α)与γ无关，γ可由最大化Q(γ)的过程推导，其中[0146] Q(γ)＝Ep(W|X；γ,α){log p(X|W；γ)}， (20)
[0147] 通过Q(γ)中的γl求偏导，可以得到γ的更新值：
[0148]
[0149] 其中Re(·)和tr(·)分别表示取实部和求迹运算。
[0150] 在本发明提出的SBL算法中，γ和α会不断地迭代更新。迭代的初始化为条件可以表示为：
[0151]
[0152] 其中 ⊙表示Hadamard积，而||·||F表示Frobenius范数。每次迭代过程中计算的复杂度主要由计算Σl时的求逆运算产生。由于N＞M，Σl的更新值可以替代为：
[0153] 基于稀疏贝叶斯学习的稳态诱发响应源定位算法性能验证仿真实验
[0154] A.计算头模型和导程场矩阵
[0155] 在稳态诱发响应源定位实验中，本发明使用了解剖模板MNI/ICBM152以实现对大脑的无偏标准估计。通过分离MRI模板包含的白质和灰质间的界面，获得了皮层表面的空间信息。在Brainstorm软件的帮助下，本发明使用了边界元方法(BEM)计算由头皮、头骨和皮质组成的3层真实头模型。通过下采样，本发明得到了1500个均匀分布在皮质表面的三角网格以表示偶极子的空间位置，假设偶极子的电流方向垂直于皮质表面，结合63导联的脑电采集空间分布信息，通过OpenMEEG软件可以计算出基于真实头模型的导程场矩阵。
[0156] B.稳态诱发响应头皮记录仿真设计：
[0157] 为合成此头模型下的稳态诱发响应头皮记录并测试稳态诱发响应源定位方法的性能，本发明需要根据已经得到的1500个皮层表面网格的空间位置设计稳态诱发响应激活源的分布，并将这些源的源活动设计为稳态诱发响应的信号形式。如图2所示，本发明共设计了3种源分布场景，分对应于1个源，2个源，3个源的情况。编号为1的源位于颞上回附近，编号为2的源位于额中回附近，编号为3的源位于颞叶上方的中央后回，仿真添加的源均位于左半脑。
[0158] 指定了稳态诱发响应激活源的位置后，本发明设计激活源信号的形式为40Hz的正弦信号，即f0＝40。为模拟自发EEG干扰，本发明还在皮层上的其他区域添加了一定能量的噪声。位于皮层表面的源信号与导程场矩阵相乘后，就可以得到仿真合成相应的稳态诱发响应头皮记录。通过数值仿真，在每种稳态诱发响应激活源分布条件下，本发明都合成了10段(即，L＝10)长度为10秒的稳态诱发响应头皮记录，用于验证稳态诱发响应源定位算法的有效性。假设短时间内稳态诱发响应源信号和自发脑电的功率不变，本发明通过设置信噪比满足一个均值固定的高斯分布，实现同一数据段稳态诱发响应源信号的强度和自发脑电功率为一个确定的值，而不同数据段中，稳态诱发响应源信号的强度和脑电噪声的强度不同。
[0159] C.稳态诱发响应源定位结果：
[0160] 为与其他SBL方法区分，本发明用SSR-SBL表示本发明中提出的基于稀疏贝叶斯学习(SBL)的稳态诱发响应(SSR)源定位算法。为证明SSR-SBL算法的性能，本发明将另外3种线性分布式源定位方法(包括MNE，sLORETA和LCMV波束形成)用于稳态诱发响应源定位以作比较。图3，图4和图5分别给出了当稳态诱发响应仿真源的个数分别为1、2和3时各种方法获得的定位结果。定位结果以灰度图的形式投射在皮层上，深色表示定位得到的稳态诱发响应源强度相对较大，浅灰色表示强度相对较小的稳态诱发响应源。本发明在这些图中仅显示了占定位结果总能量的前80％的结果。
[0161] 如图3所示，所有种方法中，LCMV方法取得的定位结果最平滑。仅有一个激活源的情况下，相比于LCMV方法，MCE和sLORETA方法能够取得更好的定位精度。就MNE和sLORETA方法而言，当SNR＝10dB时，可以大致确定源的位置在颞叶，但是不能准确给出源的个数。当SNR＝0dB时，定位结果显示的源规模远大于实际尺寸，这时很难确定源的准确位置。当SNR＝-5dB时，这3种对比算法得到的定位结果向后半脑偏移，定位结果很不理想。对比下，本发明提出的SSR-SBL算法在相同条件下，可以实现无定位误差的定位效果。
[0162] 除SSR-SBL算法外，随着仿真源数量的增加，相同信噪比条件下，其他3种源定位方法取得的定位效果越来越差。比如在稳态诱发响应仿真源的个数为2的情况下，当SNR＝0dB时，无论是MCE还是sLORETA方法取得的定位结果都无法分辨源的个数，更不用说源的尺寸。在稳态诱发响应仿真源的个数为2的情况下，当SNR＝10dB时，相比于MCE，sLORETA方法取得的定位结果更加平滑，很难区分其结果显示的稳态诱发响应源的激活区域为2个还是3个。
[0163] 综合以上结果及分析，SSR-SBL算法区别于其他线性定位方法的主要特征是，SSR-SBL算法取得的定位结果具有很强的稀疏性，定位结果中每个分离的区域都可以解释为一个可能的源。目前的仿真结果可以证明在信噪比大于-5dB时，SSR-SBL源定位方法可以实现稳态诱发响应激活源位置的无偏估计。
[0164] SSR-SBL源定位方法应用于听觉稳态诱发响应源定位
[0165] 作为一种典型的稳态诱发响应，听觉稳态诱发响应(ASSR)有广泛的临床应用前景，例如听觉稳态诱发响应在听力检测，意识评估和麻醉监测领域都有应用。探究听觉稳态诱发响应起源对解释其临床应用具有重要意义。目前很多研究认为听觉稳态诱发响应主要起源于初级听觉皮层和脑干，但尚无确切的结论。通过本发明中提出的SSR-SBL源定位算法，可以探究听觉稳态诱发响应的皮质起源。这里我们首先将SSR-SBL算法应用于40-Hz听觉稳态诱发响应源定位实验。
[0166] 本发明收集了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听觉稳态诱发响应源定位实验之前，可以通过预处理，提高数据的信噪比。
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听觉稳态诱发响应源定位实验之前，可以通过预处理，提高数据的信噪比。
[0167] 预处理过程包括对每个被试对应的10段EEG数据进行去基线校正、分帧、去除伪迹、滤波和时域平均：第一步，将EEG记录的前5秒数据用于基线校正；第二步，将基线校正后的EEG记录中每1023个数据点分成一个数据帧，每个数据帧的长度恰好是声音刺激信号的40个整周期，且数据帧之间相邻但不重叠。第三步，移除包含数据点幅度超过40uV的数据帧，以去除伪迹。第四步，将数据帧通过中心频率为40Hz，带宽为10Hz的带通滤波器进行滤波。最后，保留10个数据段中前10个数据帧。每个数据段进行5次时域平均得到待处理数据段xl(t),l＝1,2,…L。预处理后，通过对xl(t)进行FFT，得到数据段在f0频点对应的分量 其中 向量的大小为63×1。在针对单个被试的SSR-SBL模型下，将单个被试得到的10段数据对应的 扩展成大小为63×10的输入矩阵X，其中 在针对多被试的SSR-
SBL模型下，将9个被试对应的63×10的输入矩阵Xp,p＝1,2,…P结果进行再扩展，得到大小为63×90的多被试联合输入矩阵X，其中 在已知L＝10,P＝9且导程场
矩阵L已确定的情况下，可以通过SSR-SBL算法给出的迭代过程学习得到40-Hz听觉稳态诱发响应的源分布结果，本发明在迭代过程中设置迭代误差阈值ε＝10-3，最大迭代次数Niter＝1000。
SBL模型下，将9个被试对应的63×10的输入矩阵Xp,p＝1,2,…P结果进行再扩展，得到大小为63×90的多被试联合输入矩阵X，其中 在已知L＝10,P＝9且导程场
矩阵L已确定的情况下，可以通过SSR-SBL算法给出的迭代过程学习得到40-Hz听觉稳态诱发响应的源分布结果，本发明在迭代过程中设置迭代误差阈值ε＝10-3，最大迭代次数Niter＝1000。
[0168] SSR-SBL算法取得的定位结果如图6所示。多被试联合定位结果显示，双耳刺激下得到的40-Hz听觉稳态诱发响应的激活神经元主要分布在左右脑的颞叶，在前额也有一定的分布。相比于单个被试的定位结果，多被试联合定位结果更加稳定，可以假设它是单个被试的定位结果的综合，避免了在单个被试的定位结果中出现的噪声过大、激活源能量估计过小的情况。
[0169] 以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 基于消息传递算法的SCMA译码方法 | 2020-05-11 | 567 |
| 一种基于麦克风阵列的电子交警系统 | 2020-05-12 | 479 |
| 一种加权多尺度字典学习框架的行星轮轴承故障辨识方法 | 2020-05-08 | 11 |
| 一种基于压缩感知的波前探测系统 | 2020-05-08 | 42 |
| 一种穿钢系统的信道估计方法 | 2020-05-12 | 894 |
| 基于星型网络结构的三维场景数据获取方法及系统 | 2020-05-08 | 632 |
| 在大规模MIMO中基于光滑范数的压缩感知信道估计方法 | 2020-05-08 | 51 |
| 一种室外无人机的自主定位方法 | 2020-05-08 | 120 |
| 电力系统扰动信号重构方法、预处理方法、终端及系统 | 2020-05-12 | 432 |
| 一种目标检测方法、装置、稀疏阵列和成像设备 | 2020-05-08 | 918 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
QQ群二维码
意见反馈
意见反馈