一种源成像的逆问题分析方法、装置及电子设备

技术领域

本发明涉及生理信号处理技术领域，具体涉及一种源成像的逆问题分析方法、装置及电子设备。

背景技术

脑电/脑磁(EEG/MEG)等神经电生理信号具备时间分辨率高、无创、采集方便等优点，使我们可以通过颅外测得的电磁信号在毫秒级的时间尺度上探索颅内大脑神经电活动成为可能，因此，它们在临床和生活中是一种非常重要的脑信号探测/成像手段。但是EEG/MEG在颅外探测神经活动电磁信号，空间分辨率相对较低，为精准获取大脑皮层上的电活动，通常会采用源成像(M/ESI)的方法，包括正问题建模和逆问题求解，即首先将大脑离散化并将皮层空间神经元电活动建模为电流偶极子(简称偶极子)模型(这个模型也叫源模型)，进而构建皮层空间偶极子到头皮空间传感器(即头皮脑电设备的脑电电极/脑磁的传感器)的正向传导模型(以上过程称为正问题建模)，最后利用已知的M/EEG信号及构建的正向传导模型完成皮层空间电活动的估计(即逆问题求解)。源成像通过结合被试的头脑模型，一定程度上提升了M/EEG的空间分辨率，在脑机接口，癫痫灶定位等方面具有重要研究价值。

现有逆问题求解无法同时准确估计偶极子激活范围和激活强度，导致无法估计激活范围内皮层偶极子电活动传播的路径，而这个传播路径可以帮助我们找到更为准确的致痫灶/解码脑活动等。因此，同时估计偶极子的激活范围以及激活强度，精准还原皮层电活动，对进一步提升癫痫灶定位精度/脑活动解码精度具有重要意义。

发明内容

有鉴于此，本发明实施例提供了一种源成像的逆问题分析方法、装置及电子设备，以解决现有的逆问题求解方法无法同时自动准确估计偶极子激活范围及强度的缺陷。

根据第一方面，本发明实施例提供了一种源成像的逆问题分析方法，包括：

构建真实头脑模型和源模型；

基于电极/传感器的结构、所述真实头脑模型、所述源模型、以及所述电极/传感器与所述真实头脑模型的位置关系，求解得到传导矩阵；

按照空间距离对所述源模型中的偶极子进行区域化，并构建组间稀疏约束模型；

构建偶极子强度广义全变分稀疏约束模型；

将所述组间稀疏约束模型与所述偶极子强度广义全变分稀疏约束模型融合，并基于所述传导矩阵，构建逆问题求解模型；

基于所述逆问题求解模型，确定所述源模型中偶极子的激活范围及强度。

一些可选的具体实施方式中，所述组间稀疏约束模型为

一些可选的具体实施方式中，所述按照空间距离对所述源模型中的偶极子进行区域化，包括：

S201：获取当前未被区域化的所述偶极子，并从中选取一个所述偶极子作为种子点；

S202：选取所述种子点的n阶邻域内偶极子与所述种子点一起组成一个组；其中，n＝1,2,…,N，1≤N≤N

重复步骤S201和S202，直至所述偶极子区域化完成。

一些可选的具体实施方式中，所述偶极子强度广义全变分稀疏约束模型为

一些可选的具体实施方式中，所述逆问题求解模型为

其中，

一些可选的具体实施方式中，所述体现源空间偶极子激活强度具有变分特性的线性变换矩阵D的表达式为：

其中，

一些可选的具体实施方式中，所述基于所述逆问题求解模型，确定所述源模型中偶极子的激活范围及强度，包括：

将所述逆问题求解模型变形为：

s.t.z＝Dx-p,g＝Ep,h＝Gx；其中，z、g和h为定义的潜在变量，G为遍历所述源空间所有偶极子组的辅助稀疏矩阵，其行数为

引入缩放后的增广拉格朗日算子，基于交替更新的方式对变形后的所述逆问题求解模型进行求解，确定所述源模型中激活的偶极子及其强度。

一些可选的具体实施方式中，所述交替更新的规则为：

k、k+1分别为更新次数的序号，ρ为拉格朗日惩罚算子，u，v与m为所述缩放后的增广拉格朗日算子，

t

t

一些可选的具体实施方式中，所述交替更新的结束条件为：

一些可选的具体实施方式中，所述基于交替更新的方式对变形后的所述逆问题求解模型进行求解，包括：

基于更新后的x、p、z、z、h

若所述惩罚项的权重满足所述预设条件或者当前的求解次数已经达到所述第二次数阈值，则完成对变形后的所述逆问题求解模型的求解；

若所述惩罚项的权重不满足所述预设条件且当前的求解次数尚未达到所述第二次数阈值，则再次更新x、p、z、z、h

一些可选的具体实施方式中，所述惩罚项的权重所需满足的所述预设条件为：

其中，THR1为第二预设阈值，THR2为第三预设阈值，THR3为第四预设阈值，

根据第二方面，本发明实施例提供了一种源成像的逆问题分析装置，包括：

第一模型构建模块，用于构建真实头脑模型和源模型；

计算模块，用于基于电极/传感器的结构、所述真实头脑模型、所述源模型以及所述电极/传感器和所述真实头脑模型的位置关系，求解得到传导矩阵；

第二模型构建模块，用于按照空间距离对所述源模型中的偶极子进行区域化，并构建组间稀疏约束模型；

第三模型构建模块，用于构建偶极子强度广义全变分稀疏约束模型；

第四模型构建模块，用于将所述组间稀疏约束模型与所述偶极子强度广义全变分稀疏约束模型融合，并基于所述传导矩阵，构建逆问题求解模型；

确定模块，用于基于所述逆问题求解模型，确定所述源模型中偶极子激活范围及强度。

根据第三方面，本发明实施例提供了一种电子设备，包括：

存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器用于存储计算机程序，所述计算机程序被所述处理器执行时，实现上述第一方面所述的任一种源成像的逆问题分析方法。

根据第四方面，本发明实施例提供了一种计算机可读存储介质，所述计算机可读存储介质用于存储计算机程序，所述计算机程序被处理器执行时，实现上述第一方面所述的任一种源成像的逆问题分析方法。

本发明实施例提出了基于组变分稀疏约束的源成像的逆问题分析方法、装置及电子设备，本发明实施例首先考虑神经元放电的集群特性，按照偶极子的空间距离对偶极子进行区域化，并考虑组间稀疏性，刻画皮层电活动的结构稀疏先验。另外，通过偶极子强度广义全变分刻画皮层上激活区域边界处的突变性与激活区域内的平滑特性，实现偶极子激活范围及强度的准确估计。最终，以重加权的求解方式实现自动的激活范围及强度的准确估计。

附图说明

通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：

图1为本发明实施例提供的一种源成像的逆问题分析方法的流程示意图；

图2为本发明实施例提供的一种源成像的逆问题分析方法的过程示意图；

图3为本发明实施例提供的源成像的逆问题分析方法与现有两种方法的对比结果示意图；

图4为本发明实施例提供的一种源成像的逆问题分析装置的结构示意图；

图5为本发明实施例提供的一种电子设备的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。在以下各实施例的描述中，“多个”的含义是两个以上，除非另有明确具体的限定。

源成像的逆问题分析方法的目标是准确求解皮层空间(即源空间)偶极子的激活范围(位置)和强度。现有方法可分为基于等效电流偶极子模型的方法与基于分布式电流偶极子模型的方法。其中，基于等效电流偶极子模型的方法是在已知激活偶极子数量的条件下固定其位置、方向以及强度中的某些参数，调整其他未固定的参数直至其产生的M/EEG与测量到的M/EEG吻合。由此可见，等效电流偶极子模型是非线性模型，随着偶极子数量的增多，其求解变得复杂。一般情况下，此种方法只能实现一个或几个偶极子的位置及强度估计。同时，其需要将偶极子的数量作为先验知识，这在实际中是无法准确得知的。相比之下，分布式电流偶极子模型是皮层电活动与头皮脑电的线性模型。无需将偶极子数量作为先验知识，更符合实际情况。

基于分布式电流偶极子模型的方法假设大量的偶极子均匀分布在皮层空间中，根据正向传导模型构建头皮电极/传感器与偶极子的线性方程并求解出皮层空间的电活动。由于头皮电极/传感器的数量(≈10

请参阅图1和图2，本发明实施例提供了一种源成像的逆问题分析方法，包括：

S101：构建真实头脑模型和源模型；

具体的，请参阅图2，可以根据头皮脑电设备的电极/传感器结构及磁共振成像(Magnetic Resonance Imaging，MRI)得到的图像，基于fieldtrip构建个性化的真实头脑模型；使用freesurfer等软件，根据MRI图像，构建源模型；该源模型可以是分布式源模型，分布式源模型指皮层放电源(电流偶极子，简称“偶极子”)均匀分布在三角网格化的皮层上。同时，我们假设一个三角网格内有一个偶极子且偶极子位于三角网格的重心，方向为三角网格的法向；

S102：基于电极/传感器的结构、所述真实头脑模型、所述源模型、以及所述电极/传感器与所述真实头脑模型的位置关系，求解得到传导矩阵；

S103：按照空间距离对所述源模型中的偶极子进行区域化，并构建组间稀疏约束模型；这里，基于神经元的集群放电特性，认为相邻偶极子更有可能同时激活，因此按照空间距离对偶极子进行区域化。

S104：构建偶极子强度广义全变分稀疏约束模型；

S105：将所述组间稀疏约束模型与所述偶极子强度广义全变分稀疏约束模型融合，并基于所述传导矩阵，构建逆问题求解模型；

S106：基于所述逆问题求解模型，确定所述源模型中偶极子的激活范围及强度。具体的，可以以重加权的方式求解所述逆问题求解模型，自动确定所述源模型中偶极子的激活范围及强度。

本发明实施例提出了基于组变分稀疏约束的源成像的逆问题分析方法，该方法首先考虑神经元放电的集群特性，按照偶极子的空间距离对偶极子进行区域化，并考虑组间稀疏性，刻画皮层电活动的结构稀疏先验。另外，通过偶极子强度广义全变分刻画皮层上激活区域边界处的突变性与激活区域内的平滑特性，实现偶极子激活范围及强度的准确估计。最终，以重加权的求解方式实现自动的激活范围及强度估计。

一些具体的实施方式中，所述组间稀疏约束模型为

本发明实施例中，在按照空间距离对所述源模型中的偶极子进行区域化时，会将源模型中的偶极子划分为I个组。

一些具体的实施方式中，所述按照空间距离对所述源模型中的偶极子进行区域化包括：

S201：获取当前未被区域化的所述偶极子，并从中选取一个所述偶极子作为种子点；

S202：选取所述种子点的n阶邻域内偶极子与所述种子点一起组成一个组；其中，n＝1,2,…,N，1≤N≤N

重复步骤S201和S202，直至所述偶极子区域化完成。

本发明实施例中，以偶极子的邻域来表征偶极子的空间距离。当N为10时，即种子点以及种子点的一阶邻域偶极子、二阶邻域偶极子、三阶邻域偶极子、四阶邻域偶极子、五阶邻域偶极子、六阶邻域偶极子、七阶邻域偶极子、八阶邻域偶极子、九阶邻域偶极子和十阶邻域偶极子为一个偶极子分组。

一些具体的实施方式中，所述偶极子强度广义全变分稀疏约束模型为

可选的，为同时保证偶极子强度一阶变分与二阶变分特性，α可以设定为0.5。

一些具体的实施方式中，所述逆问题求解模型为

其中，/>

可选的，δ＝10

一些具体的实施方式中，所述体现源空间偶极子激活强度具有变分特性的线性变换矩阵D的表达式为：

/>

其中，

一些具体的实施方式中，所述基于所述逆问题求解模型，确定所述源模型中偶极子的激活范围及强度，包括：

将所述逆问题求解模型变形为：

s.t.z＝Dx-p,g＝Ep,h＝Gx；其中，z、g和h为定义的潜在变量，G为遍历所述源空间所有偶极子组的辅助稀疏矩阵，其行数为

引入缩放后的增广拉格朗日算子，基于交替更新的方式对变形后的所述逆问题求解模型进行求解，确定所述源模型中激活的偶极子及其强度。

本发明实施例中，利用交替方向乘子法(Alternating Direction Method ofMultipliers，ADMM)方法求解逆问题求解模型，以得到偶极子的激活情况。

其他的可替换实施方式中，也可以采用深度学习网络求解上述的逆问题求解模型。

其中一些具体的实施方式中，所述交替更新的规则为：

k、k+1分别为更新次数的序号，ρ为拉格朗日惩罚算子，u，v与m为所述缩放后的增广拉格朗日算子；

t

t

其他可选的具体实施方式中，上述交替更新的规则也可以用深度神经网络替代，例如FISTA(一种快速的迭代阈值收缩算法)网络等。

一些具体的实施方式中，所述交替更新的结束条件为：

例如，所述交替更新的结束条件为：

一些具体的实施方式中，所述基于交替更新的方式对变形后的所述逆问题求解模型进行求解，包括：

基于更新后的x、p、z、z、h

若所述惩罚项的权重满足所述预设条件或者当前的求解次数已经达到所述第二次数阈值，则完成对变形后的所述逆问题求解模型的求解；

若所述惩罚项的权重不满足所述预设条件且当前的求解次数尚未达到所述第二次数阈值，则再次更新x、p、z、z、h

其中一些具体的实施方式中，所述惩罚项的权重所需满足的条件为：

其中，THR1为第二预设阈值，THR2为第三预设阈值，THR3为第四预设阈值，

本发明实施例中，基于求解出的

例如，THR1、THR2和THR3均可以为10

本文选择一类经典且与本发明实施例提供的逆问题求解方法强相关的逆问题求解方法，包括：最小一范数融合变分方法(TV-l

为验证所提方法的准确性，我们基于获取的头脑模型源模型，分别在无噪声、信噪比为30dB、20dB、10dB以及信噪比为5dB条件下分别构建激活范围相同的50个仿真头皮脑电。其中，仿真的激活范围为种子点偶极子的5阶邻域(如权利要求3所述)，激活范围内偶极子激活强度呈以种子点偶极子为中心，均值为0，标准差为

偶极子范围定位误差(Dipole Localization Error,DLE)及相对均方根误差(Normalized Root Mean Square Error,NMSE)作为方法性能的评价指标分别评估所述方法估计偶极子激活范围及强度的准确性。其中，DLE的计算方式为：

其中，E代表DLE，N

NMSE的计算方式为：

其中，η代表NMSE，

相应地，请参考图4，本发明实施例提供一种源成像的逆问题分析装置，包括：

第一模型构建模块301，用于构建真实头脑模型和源模型；

计算模块302，用于基于电极/传感器的结构、所述真实头脑模型、所述源模型、以及所述电极/传感器与所述真实头脑模型的位置关系，求解得到传导矩阵

第二模型构建模块303，用于按照空间距离对所述源模型中的偶极子进行区域化，并构建组间稀疏约束模型；

第三模型构建模块304，用于构建偶极子强度广义全变分稀疏约束模型；

第四模型构建模块305，用于将所述组间稀疏约束模型与所述偶极子强度广义全变分稀疏约束模型融合，并基于所述传导矩阵，构建逆问题求解模型；

确定模块306，用于基于所述逆问题求解模型，确定所述源模型中偶极子的激活范围及强度。

本发明实施例提出了基于组变分稀疏约束的源成像的逆问题分析装置，该装置首先考虑神经元放电的集群特性，按照偶极子的空间距离对偶极子进行区域化，并考虑组间稀疏性，刻画皮层电活动的结构稀疏先验。另外，通过偶极子强度广义全变分刻画皮层上激活区域边界处的突变性与激活区域内的平滑特性，实现偶极子激活范围及强度的准确估计。最终，以重加权的求解方式实现自动的激活范围及强度估计。

一些具体的实施方式中，所述组间稀疏约束模型为

一些具体的实施方式中，所述第二模型构建模块303包括：

种子点选取单元，用于获取当前未被区域化的所述偶极子，并从中选取一个所述偶极子作为种子点；

分组单元，用于选取所述种子点的n阶邻域内偶极子与所述种子点一起组成一个组；其中，n＝1,2,…,N，1≤N≤N

控制单元，用于控制种子点选取单元和分组单元重复执行，直至所述偶极子区域化完成。

一些具体的实施方式中，所述偶极子强度广义全变分稀疏约束模型为

一些具体的实施方式中，所述逆问题求解模型为

一些具体的实施方式中，所述体现源空间偶极子激活强度具有变分特性的线性变换矩阵D的表达式为：

其中，

一些具体的实施方式中，所述确定模块306包括：

变形单元，用于将所述逆问题求解模型变形为：

，

s.t.z＝Dx-p,g＝Ep,h＝Gx；其中，z、g和h为定义的潜在变量，G为遍历所述源空间所有偶极子组的辅助稀疏矩阵，其行数为

求解单元，用于引入缩放后的增广拉格朗日算子，基于交替更新的方式对变形后的所述逆问题求解模型进行求解，确定所述源模型中激活的偶极子及其强度。

一些具体的实施方式中，所述交替更新的规则为：

k、k+1分别为更新次数的序号，ρ为拉格朗日惩罚算子，u，v与m为所述缩放后的增广拉格朗日算子，

t

一些可选的具体实施方式中，所述交替更新的结束条件为：

一些具体的实施方式中，所述求解单元，具体用于：

基于更新后的x、p、z、z、h

若所述惩罚项的权重满足所述预设条件或者当前的求解次数已经达到所述第二次数阈值，则完成对变形后的所述逆问题求解模型的求解；

若所述惩罚项的权重不满足所述预设条件且当前的求解次数尚未达到所述第二次数阈值，则再次更新x、p、z、z、h

一些具体的实施方式中，所述惩罚项的权重所需满足的条件为：

其中，THR1为第二预设阈值，THR2为第三预设阈值，THR3为第四预设阈值，

本发明实施例是与上述方法实施例基于相同的发明构思的装置实施例，因此具体的技术细节和对应的技术效果请参阅上述方法实施例，此处不再赘述。

本发明实施例还提供了一种电子设备，如图5所示，可以包括处理器41和存储器42，其中处理器41和存储器42可以通过总线或者其他方式互相通信连接，图5中以通过总线连接为例。

处理器41可以为中央处理器(Central Processing Unit，CPU)。处理器41还可以为其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等芯片，或者上述各类芯片的组合。

存储器42作为一种非暂态计算机可读存储介质，可用于存储非暂态软件程序、非暂态计算机可执行程序以及模块，如本发明实施例中的源成像的逆问题分析方法对应的程序指令/模块(例如，图4所示的第一模型构建模块301、计算模块302、第二模型构建模块303、第三模型构建模块304、第四模型构建模块305和确定模块306)。处理器41通过运行存储在存储器42中的非暂态软件程序、指令以及模块，从而执行处理器的各种功能应用以及数据处理，即实现上述方法实施例中的源成像的逆问题分析方法。

存储器42可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储处理器41所创建的数据等。此外，存储器42可以包括高速随机存取存储器，还可以包括非暂态存储器，例如至少一个磁盘存储器件、闪存器件、或其他非暂态固态存储器件。在一些实施例中，存储器42可选包括相对于处理器41远程设置的存储器，这些远程存储器可以通过网络连接至处理器41。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。

所述一个或者多个模块存储在所述存储器42中，当被所述处理器41执行时，执行如图1-2所示实施例中的源成像的逆问题分析方法。

上述电子设备的具体细节可以对应参阅图1至图2所示的实施例中对应的相关描述和效果进行理解，此处不再赘述。

相应地，本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质用于存储计算机程序，所述计算机程序被处理器执行时，实现上述源成像的逆问题分析方法实施例的各个过程，且能达到相同的技术效果，为避免重复，这里不再赘述。

计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。