一种基于最近邻多元线性回归模型的脑磁源定位方法

1.本发明设计生物医学工程领域，具体涉及一种基于最近邻多元线性回归模型的脑磁源定位方法。


背景技术：

2.脑磁图(meg:magnetoencephalogram)是一种用于探究人脑活动的非侵入式功能成像技术，在近几十年里是神经科学和临床研究中备受关注。相较于一些其他功能成像技术(如脑电图和近红外成像等)，脑磁图以其高时间和空间分辨率的优势，能更精准地实现对大脑活动的定位。脑磁源定位是指从传感器处获得的脑磁信号中重建出潜在的大脑活动。目前常见meg源定位有两大类方法：偶极子拟合方法和分布式源成像方法。偶极子拟合方法假设大脑中的偶极子数量有限，以偶极子代表源活动，然后用最小二乘法拟合估计这些偶极子参数。在分布式源成像方法中，皮质表面被划分为一定数量的三角形网格，每个网格代表一个偶极子。与偶极子拟合方法相比，成像方法事先确定了源的数量，只需要解决一个线性问题。由于未知参数的数量大大超过方程，这个线性问题仍有许多解。应加入适当的先验约束以获得唯一的解。
3.传统的分布式源成像方法通常采用l2范数约束，如最小规范估计(mne)。mne对电流源的l2范数进行约束以获得最小能量的解决方案。在传感器信号层面上，浅层源产生的磁场强度比深层源大，所以mne偏向于浅层源。mne的变种已经出现在现有的一些种类的文献中，包括加权mne(wmne)，以及低分辨率脑电磁断层扫描(loreta)等。虽然基于l2范数的方法计算效率高，但它们的估计源是扩散的，覆盖了实际皮质活动区域以外的大部分大脑皮层，其结果是许多来源不能在空间上被分离。
4.准确估计源范围对研究大脑活动是必要的。估算癫痫致病脑组织和绘制脑功能区图是其重要的应用。为了实现准确的源范围估计，本发明提出了一种基于最近邻的多元线性回归算法来实现源范围估计。图拉普拉斯空间平滑器被用来对源协方差进行建模，源协方差利用凸优化理论不断地被更新。


技术实现要素：

5.本发明要解决的技术问题是：传统分布式源成像方法无法准确定位出源的范围大小。为克服现有技术的不足，本发明提供了一种基于最近邻多元线性回归模型的脑磁源定位方法，快速准确地实现了脑磁源定位，并能较为准确地估计出源的范围大小。
6.本发明所述的基于最近邻多元线性回归模型的脑磁源定位方法包括以下步骤：
7.步骤1：在良好的磁屏蔽环境下，先获得实验开始之前的空房间信号数据计算出噪声协方差矩阵σ
ε
，然后开展实验以获取质量较好的脑磁信号数据b，并对获得的脑磁信号数据b进行预处理。
8.步骤2：首先使用核磁共振仪采集受试者的mri数据。对mri图像数据进行分割，把被试的头部分割为：头皮、颅骨、脑脊液、灰质和白质等。之后创建单壳模型头模型，再将头
模型坐标系和脑磁图传感器阵列坐标系进行配准并进一步计算得出引导场矩阵l。
9.步骤3：将经过预处理之后的脑磁信号数据b、噪声协方差矩阵σ
ε
以及引导场矩阵l输入最近邻多元线性回归模型中进行源定位。
10.步骤3.1：先初始化模型参数，利用被试mri数据根据公式4得到图拉普拉斯空间平滑器g，再将分别从步骤1和步骤2中得到的脑磁信号b、噪声协方差矩阵σ
ε
和引导场矩阵l代入最近邻多元线性回归模型之中。
11.步骤3.2：先使用公式5更新源信号估计值s，再使用公式7和8分别更新参数a和γ，然后根据公式6计算损失函数的大小。
12.步骤3.3：判断步骤3.2计算得到的损失函数的大小是否收敛。若不收敛，重复步骤3.2，若收敛，此时源信号的估计值s即为最终的源估计值。
13.本发明所实现的源定位算法能获得更为准确的源估计，包括准确的定位源位置、范围大小以及源时间波形估计，为脑磁的医学应用提供了坚实的技术基础。
附图说明
14.图1脑磁源定位流程图。
15.图2单壳模型头模型。
16.图3基于最近邻多元线性回归模型脑磁源定位原理流程图。
具体实施方式
17.本发明提出一种基于最近邻多元线性回归模型的脑磁源定位方法，下面结合图1、图2和图3具体阐释算法原理和实现步骤：
18.在介绍具体操作步骤之前简单介绍下脑磁信号数据b、头模型导出的引导场l、源定位结果s以及噪声ε之间的关系：
19.b＝ls+ε (1)
20.其中，采集的脑磁信号表示为nc代表通道数，n
t
表时间采用点的个数，源信号表示为ns表示分布式源模型中所假设的源个数，引导场是神经元活动产生的磁源活动在真实头模型下传导规律的数学表示，噪声表示为
21.如下图1所示的是本发明的具体步骤，概括如下：
22.步骤1：在磁屏蔽房或磁屏蔽桶等磁屏蔽环境下，使用传感器(如超导量子干涉仪(squid)和光泵原子磁力计(quspin))测量空房间信号和脑磁信号数据b。在实验开始前需要测量空房间信号以便计算出噪声协方差矩阵σ
ε
，再开展脑磁实验，获取质量较好的脑磁信号数据b，对获得的脑磁信号数据b进行一些必要的预处理工作，包括识别并剔除坏道以及使用带通滤波器和叠加平均等手段去除噪声。经过上述预处理后，脑磁信号数据才可以用于后续的源定位。
23.步骤2：首先使用核磁共振仪采集受试者的mri数据。然后将受试者的mri数据导入freesurfer软件之中，并使用该软件对mri图像数据进行分割，把被试的头部分割为：头皮、颅骨、脑脊液、灰质和白质等。之后使用matlab软件创建单壳模型头模型,如图2所示，只需
要设定好皮质的电导率即可。构建好合适的头模型之后，将头模型坐标系和传感器阵列坐标系进行配准并进一步计算得出引导场矩阵l。
24.步骤3：将经过预处理之后的脑磁信号数据b、噪声协方差矩阵σ
ε
以及引导场矩阵l输入最近邻多元线性回归模型中进行源定位，算法具体的流程如图3所示。
25.在介绍具体操作步骤3之前，先介绍下最近邻多元线性回归模型的框架及原理：
26.首先介绍最近邻的概念，如图2所示。被试mri数据中的皮质部分可以用三角网格分割，得到一系列点云数据，方便重建。mri分割时一般将距离较近且处于同一组织的源构建一个连接关系，而最近邻的概念就从中产生。在图2中，左边是分割的三角网格，右边是局部放大图，显示了中心蓝色点源和其周围的7个有着连接关系且距离较近的绿色点源，我们将绿色点源称之为蓝色点源的最近邻。
27.大脑各神经元活动之间是存在时间和空间上的联系的，那么最直观的便是源和自己的最近邻之间肯定存在时空联系，本发明拟用最近邻多元线性回归模型来描述这种时空联系，如下：
[0028][0029]
其中，λ是一个介于0到1之间的常数，本发明取0.95，an表征上一个时刻源自身和其最近邻的不同贡献，常见的分布式源模型中含有近万个源，所以an也有近万个，这会导致模型过拟合，所以简化为a，即不同的源对应的an相同，n(n)表示第n个源的最近邻集合，d
n,i
是一个比例系数，表征不同最近邻对源活动的影响程度，一般取d
n,i
＝1/nn，nn指最近邻的个数，w
n,t
是未知扰动，用来修正模型偏差，一般认为在所有源中都一样的量，简化为w
t
，上式可以简化为：
[0030][0031]
其中，f矩阵表征源之间的时空联系，i是单位矩阵,q矩阵是加权邻接矩阵，能表征各源与其最近邻源之间的相关程度，m是邻接矩阵，表征各源之间是否处于相邻状态,diag(h)表示将向量h变成一个对角矩阵的操作，对角元素就是向量h的值.为了更好的求解源定位结果，我们假设a服从一个均值为μ，方差为β的高斯分布，未知扰动w
t
也服从一个高斯分布，均值为0，方差为σw。最常见的是假设方差σw是一个对角矩阵，而为了实现更好的源范围估计，本发明采用图拉普拉斯来构建方差σw，具体如下：
[0032][0033]
其中，g是图拉普拉斯空间平滑器，ns表示源偶极子的个数。我们设置方差矩阵γ＝diag(γi),其中我们假设γi相互独立，gi是g的第i列。
[0034]
有了未知扰动之后，我们可以按照下式迭代求解源估计：
[0035][0036]
其中，表示已知的上个时刻的源估计，因此为了迭代能进行下去需要提前已知0时刻源的值，n
t
表示采样点的个数。然而，上式中的f以及σw并不是已知的，而f和σw分别由参数a和γ完全决定，因此需要将参数a和γ的值求解出来。为此，我们由贝叶斯推理可得损失函数如下：
[0037][0038]
其中，σb＝lσwl
t
+σ
ε
，r＝b-lfs-1
，|
·
|是行列式符号，log表示取对数，trace是求矩阵迹的操作，是一个和参数a和γ均无关的常数。我们将损失函数对参数a求导置0可得a的迭代表达式如下：
[0039][0040]
矩阵和是计算a时的过渡矩阵，无实际意义，和分别是矩阵和的第t列。由于对数和行列式的影响，我们将损失函数对参数γ求导置0得不到参数γ的显示迭代表达式。由于损失函数是参数γ的凸函数，所以我们采用凸分析来解决这个问题，得到γ的迭代表达式如下：
[0041][0042]
上式中是第i个源及其最近邻的集合，(
·
)
i,j
表示矩阵的第i行、第j列元素,(s-fs-1
)i是矩阵s-fs-1
的第i列，||
·
||2表示对矩阵取二范数。
[0043]
接下来结合图3阐述本发明中步骤3具体实现：
[0044]
步骤3.1：先初始化参数a＝0.5和γi＝1，再设置0时刻的源s0为一个0向量，利用被试mri数据根据公式4得到图拉普拉斯空间平滑器g，再将分别从步骤1和步骤2中得到的脑磁信号b、噪声协方差矩阵σ
ε
和引导场矩阵l代入最近邻多元线性回归模型之中。
[0045]
步骤3.2：先使用公式5更新源信号估计值s，再使用公式7和8分别更新参数a和γ，然后根据公式6计算损失函数的大小。
[0046]
步骤3.3：判断步骤3.2计算得到的损失函数的大小是否收敛。若不收敛，重复步骤3.2，若收敛，此时源信号的估计值s即为最终的源估计值。
[0047]
本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。
[0048]
提供以上实施例仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

技术特征：
1.一种基于最近邻多元线性回归模型的脑磁源定位方法，其特征在于，包括以下步骤：步骤1：利用脑磁实验开始之前的空房间信号数据计算出噪声协方差矩阵σ
ε
，获得脑磁实验开始之后的脑磁信号数据b，对数据进行预处理；步骤2：采集受试者的mri数据，利用mri数据构建单壳头模型，将各坐标系进行配准并进一步计算得出引导场矩阵l；步骤3：将经过预处理之后的脑磁信号数据、噪声协方差矩阵以及引导场矩阵输入最近邻多元线性回归模型中进行源定位。2.根据权利要求1所述的脑磁源定位方法，其特征在于，所述步骤1包括：在高度磁屏蔽的环境下，使用传感器测量得到受试者的脑磁信号数据，对获得的脑磁信号数据进行预处理工作，所述预处理工作包括识别并剔除坏道以及去除噪声。3.根据权利要求2所述的脑磁源定位方法，其特征在于，所述步骤2包括：使用核磁共振仪采集受试者的mri数据，将受试者的mri数据导入分割软件中进行分割，把头部皮质分割为：头皮、颅骨、脑脊液、灰质和白质，然后设定上述各皮质的电导率，构建单壳头模型，将头模型坐标系和传感器阵列坐标系下的数据配准到同一坐标系下，并进一步计算得出引导场矩阵。4.根据权利要求1所述的脑磁源定位方法，其特征在于，所述步骤3包括：步骤3.1：先初始化参数a和γ，再设置0时刻的源为一个零向量，利用mri数据根据公式(1)得到图拉普拉斯空间平滑器，其中，g是图拉普拉斯空间平滑器，n
s
表示源的个数，i是单位矩阵,m是邻接矩阵，再将分别从步骤1和步骤2中得到的脑磁信号数据b、噪声协方差矩阵σ
ε
和引导场矩阵l代入最近邻多元线性回归模型之中，步骤3.2：先使用公式(2)更新源信号估计值s，其中，s
t
是t时刻源s的值，表示已知的上个时刻的源估计，σ
w
是算法模型假设的源协方差矩阵，其中假设γ
i
相互独立，g
i
是g的第i列，b
t
是t时刻的脑磁传感器数据，f矩阵表征源之间的时空联系，n
t
表示采样点的个数；再使用公式(3)和(4)分别更新参数a和γ，μ和β是本算法模型假设的参数a的均值和方差，λ是一个介于0到1之间的常数，和
分别是矩阵和的第t列，σ
b
＝lσ
w
l
t
+σ
ε
，trace是求矩阵迹的操作，矩阵和是计算参数a时的过渡矩阵，无实际意义，q矩阵是一个加权邻接矩阵，能表征各源与其最近邻源之间的相关程度；上式中是第i个源及其最近邻的集合，(
·
)
i,j
表示矩阵的第i行、第j列元素,(s-fs-1
)
i
是矩阵s-fs-1
的第i列，||
·
||2表示对矩阵取二范数，n
t
表示采样点的个数；然后根据公式(5)计算损失函数的大小；其中，n
t
表示采样点的个数，σ
b
＝lσ
w
l
t
+σ
ε
，r＝b-lfs-1
，|
·
|是行列式符号，log表示取对数，trace是求矩阵迹的操作，是一个和参数a和γ均无关的常数；步骤3.3：判断步骤3.2计算得到的损失函数的大小是否收敛,若不收敛，重复步骤3.2，若收敛，此时源信号的估计值s即为最终的源估计值。5.根据权利要求4所述的脑磁源定位方法，其特征在于，所述步骤3.1包括：在最近邻多元线性回归模型之中，使用公式(6)中的基于最近邻多元线性回归模型来构建大脑模型并表征出大脑源活动之间的时空关联；其中，λ是一个介于0到1之间的常数，一般取0.95，a
n
表征上一个时刻源自身和其最近邻的不同贡献，n(n)表示第n个源的最近邻集合，d
n,i
是一个比例系数，表征不同最近邻对源活动的影响程度，一般取d
n,i
＝1/n
n
，n
n
指最近邻的个数，w
n,t
是未知扰动，用来修正模型偏差；使用上述公式(1)中的基于图拉普拉斯平滑器构建模型中的源协方差矩阵，表征源活动之间协方差层级的空间相关性；使用基于贝叶斯推理和凸优化理论对模型进行合理的假设并推导出用于模型参数迭代求解的上述公式(3)和(4)。

技术总结
本发明提供一种基于最近邻多元线性回归模型的脑磁源定位方法。方法包括：基于最近邻多元线性回归模型来构建大脑模型并表征出大脑源活动之间的时空关联；基于图拉普拉斯平滑器构建模型中的源协方差矩阵，表征源活动之间协方差层级的空间相关性；基于贝叶斯推理和凸优化理论对模型进行合理的假设并推导出模型参数迭代求解的公式。本发明所实现的源定位算法能获得更为准确的源估计，包括准确的定位源位置、范围大小以及源时间波形估计，为脑磁的医学应用提供了坚实的技术基础。医学应用提供了坚实的技术基础。医学应用提供了坚实的技术基础。


技术研发人员：宁晓琳 李文 高阳
受保护的技术使用者：北京航空航天大学
技术研发日：2023.05.31
技术公布日：2023/8/21