一种基于贝叶斯的儿童脑活动异常病灶区成像系统

1.本发明属于自闭症儿童脑病灶区重构技术领域，更具体地，涉及一种基于贝叶斯的儿童脑活动异常病灶区成像系统。


背景技术：

2.大脑作为中枢神经系统的最重要部分，它控制着人的高级思维，如认知、学习以及社交等，自闭症儿童的教育与脑功能障碍干预治疗首先需要弄清自闭症儿童认知学习过程与脑功能障碍的大脑运行机理。脑电源成像是探究大脑活动机理、理解各种思维和感知脑内活动动态变化过程的重要手段，对揭示自闭症儿童大脑运行机理有重要的意义。
3.脑电源成像逆问题是通过采样脑电磁数据反演脑内神经元活动信息的过程，是脑电源成像的核心。近年来可用于自闭症儿童认知学习以及功能障碍脑电源成像技术国内外研究进展显著，但现存脑源成像技术大多基于最小范数的方法，自由度有限，且人脑源活动的复杂度强、分辨率需求高的特点导致自闭症儿童认知学习与功能障碍脑电源成像精准性低下。如何实现脑电磁异常采样数据的自动识别化，快速和高精度的重构自闭症儿童脑病灶区活动及其位置仍然是一个充满挑战的课题。


技术实现要素：

4.针对现有技术的至少一个缺陷或改进需求，本发明提供了一种基于贝叶斯的儿童脑活动异常病灶区成像系统，可以极大地提高病灶区定位和重构的精准性，特别适用于自闭症儿童的儿童脑电活动异常病灶区成像。
5.为实现上述目的，本发明提供了一种基于贝叶斯的儿童脑活动异常病灶区成像系统，包括：
6.信号处理模块，用于获取脑电磁采样数据，对脑电磁采样数据进行预处理；
7.脑活动异常识别模块，用于从预处理后的脑电磁采样数据中识别出异常波和非异常波；
8.噪声估计模块，用于构建非异常波的第一采样数据生成模型，采用变分贝叶斯方法对第一采样数据生成模型进行求解，获取数据采集的背景干扰噪声的分布信息；
9.脑电源活动重构模块，用于构建异常波的第二采样数据生成模型，背景干扰噪声的分布信息是第二采样数据生成模型中的参数，采用经验贝叶斯方法，对第二采样数据生成模型进行求解，获取脑电源活动时间序列；
10.成像模块，用于根据脑电源活动时间序列进行成像。
11.优选的，所述信号处理模块包括：
12.移除采样偏移模块，用于移除脑电磁采样数据中的采样偏移；
13.去干扰模块，用于采用独立成分分析方法和空间投影方法对移除采样偏移后的信号进行去噪；
14.滤波模块，用于对去噪处理后的信号进行滤波。
15.优选的，所述第一采样数据生成模型为：
16.y
con
＝fu+e，其中，y
con
是非异常波，f是未知的m
×
k的加权矩阵，u是未知的k
×
t
k
的因子参数，t
k
为非异常波中的采样点数目，u＝[u1,u2,...,u
k
]
t
，e为非异常波分解后残余信号，m为采集点的通道数目，设k＜＜m。
[0017]
优选的，所述采用变分贝叶斯方法对第一采样数据生成模型进行求解包括步骤：
[0018]
设因子参数u的先验分布为标准高斯分布：
[0019][0020]
数据采集的背景干扰噪声e的分布为：
[0021][0022]
其中，λ＝diag(λ1,λ2,...,λ
m
)是对角矩阵，λ为e的先验方差，
[0023]
加权矩阵f设为高斯分布：
[0024][0025]
其中，f
m,k
为加权矩阵f中的值，α
k
为加权矩阵f的先验方差，1≤k≤k，基于变分贝叶斯估计对第一采样数据生成模型进行演变，可得加权矩阵f和因子参数u的更新方式为：
[0026][0027][0028]
ψ＝r
uu
+α
[0029][0030][0031]
其中r
yy
为采样数据协方差矩阵，r
uu
为估计因子协方差矩阵，r
uy
为估计因子和采样数据相乘，α＝diag(α1,α2...,α
k
)；
[0032]
按照加权矩阵f和因子参数u的更新方式进行迭代循环直到加权矩阵f和因子参数u的值不再改变，根据此时加权矩阵f和因子参数u的值获取降维的和背景干扰噪声的协方差矩阵σ
e
为：
[0033][0034]
优选的，所述第二采样数据生成模型为：
[0035]
y
act
＝ls+ε，
[0036]
其中，y
act
是异常波，y
act
＝[y(t1),y(t2),...,y(t
k
)]，t
k
为异常波的采样点数目，为获取的已知的单位源活动产生的脑外采样强度,n是待求解源数目，s为脑电源活动时间序列，ε为背景干扰噪声，位置n的单位源活动产生的脑外采样强度为l
n
＝[l1,l2,...,l
m
]，n＜n，m为采集点的通道数目，m＜＜n，s＝[s1,s2,...,s
n
]
t
。
[0037]
优选的，所述采用经验贝叶斯方法对第二采样数据生成模型进行求解包括步骤：
[0038]
设脑电源活动时间序列s在空间和时间上相互独立且服从高斯分布，采用经验贝叶斯方法将第二采样数据生成模型转换为第二采样数据概率生成模型；
[0039]
构建基于凸函数平面边界的代价函数，通过对代价函数进行优化来对第二采样数据概率生成模型进行求解。
[0040]
优选的，所述第二采样数据概率生成模型为：
[0041]
脑电源活动时间序列s在空间和时间上相互独立且服从高斯分布：
[0042][0043]
为位置n的先验高斯分布的方差，则条件概率分布p(y|s)为，
[0044][0045]
σ
e
是噪声估计模块获取的背景干扰噪声的协方差矩阵，基于经验贝叶斯模型可得：
[0046][0047]
后验概率分布p(s|y
act
)为：
[0048][0049]
γ
‑1为超参数，和γ
‑1分别为t时刻采样源活动的后验概率均值和方差。
[0050]
优选的，超参数γ
‑1的求解通过最大化边缘分布p(y
act
|γ)，边缘分布对数模型logp(y
act
|γ)为：
[0051][0052]
为异常波协方差矩阵模型，是第n个源活动点的构建基于γ
‑1的log|σ
y
|边界函数为，
[0053][0054]
tr(*)是计算矩阵的迹；
[0055]
基于凸函数平面边界的代价函数为：
[0056][0057]
λ
n
与λ0为辅助变量；
[0058]
由第二采样数据概率生成模型可得后验分布p(s|y
act
)的均值与方差如下：
[0059][0060][0061]
对代价函数求超参数γ
n
和辅助参数λ
n
倒数，可得：
[0062][0063][0064]
使为零，超参数γ
n
和辅助参数λ
n
的公式如下：
[0065][0066][0067]
通过循环迭代γ
‑1和γ
n
，直至边缘分布对数模型收敛，此时的即为求解脑电源活动时间序列。
[0068]
优选的，所述识别出异常波和非异常波包括步骤：
[0069]
根据预处理后的电磁采样信号的振幅识别第一候选异常波；
[0070]
读取预先存储的标准异常波，将第一候选异常波的局部和标准异常波进行匹配，根据匹配结果来识别出异常波。
[0071]
优选的，所述根据脑电源活动时间序列进行成像包括步骤：
[0072]
获取预先存储的模型参数，根据模型参数和脑电源活动时间序列进行成像。
[0073]
总体而言，本发明与现有技术相比，具有有益效果：通过在儿童脑活动异常病灶区成像系统中采用基于变分贝叶斯的脑活动背景干扰噪声提取方法，以及与采样异常脑电磁数据结合并基于经验贝叶斯模型对儿童异常脑活动病灶区进行重构和定位，由于贝叶斯估计的方法是基于概率模型的匹配，可自适应根据采样数据特征来估计脑活动异常区的脑源活动参数，极大提高了病灶区定位和源活动重构的精准性。特别适用于自闭症儿童脑活动异常病灶区成像。
附图说明
[0074]
图1是本发明实施例的成像系统的组成示意图；
[0075]
图2是本发明实施例的成像系统的工作原理图。
具体实施方式
[0076]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0077]
如图1所示，本发明实施例的基于贝叶斯的儿童脑活动异常病灶区成像系统包括信号采集处理模块、脑活动异常识别模块、噪声估计模块、脑电源活动重构模块和成像模块。本系统可适用于各类通过脑电磁活动异常病灶区成像来识别疾病的场景中，尤其适用于儿童自闭症。
[0078]
信号处理模块，用于获取电磁采样数据，并对电磁采样数据进行预处理。
[0079]
下面说明信号处理模块的优选实现方式。
[0080]
信号采集：
[0081]
输入脑电磁采样数据，y
raw
＝[y
raw
(1),y
raw
(2),...,y
raw
(t)],t是原始信号采样点数目，y
raw
(t)＝[y1(t),y2(t),
…
,y
m
(t)]
t
，这里m为采集点的通道数目，t表示转置矩阵。需要对输入数据做初步筛查，保证采样数据的有效性以及完整性，如数据格式转换以及参数确认，包括采样数据信息等。
[0082]
信号预处理：
[0083]
信号采集处理模块包括移除采样偏移模块、去干扰模块和滤波模块，对输入数据做预处理，包括移除采样偏移，去燥，滤波等预处理。
[0084]
移除采样偏移的公式如下：
[0085][0086]
去干扰需要根据采样数据特征来进行。
[0087]
优选地，基于独立成分分析(ica)或空间投影等方式，如基于ica的去干扰(包括眼电等)，公式如下：
[0088][0089]
其中a＝[a1,a2,...,a
n
]为分解加权矩阵，x＝[x1,x2,...,x
n
]
t
为独立成分，设x
m
包含显著干扰或者噪声的特征，则去燥后采样数据y
new
如下：
[0090][0091]
优选地，利用基于投影的去干扰方法对经过基于ica的去干扰处理后的信号进一
步进行处理。基于投影的去干扰方法则需要找出投影矩阵，其可以基于脑结构空间信息等特征，假设投影矩阵为w，其去干扰公式为：
[0092]
y
clean
＝w
t
y
new
ꢀꢀ
(4)
[0093]
滤波:对采样数据的滤波去除非脑活动频段的采样数据，主要应用基于傅立叶或者小波变换的时频滤波，保留感兴趣频段数据。通过预处理后的数据保证了后续异常识别精准性。
[0094]
脑活动异常识别模块，用于从预处理后的脑电磁采样数据中识别出异常波和非异常波。
[0095]
优选的，脑活动异常识别是基于采样脑电数据的振幅以及局部波形来识别。
[0096]
根据振幅初选：对采样脑电数据的振幅进行统计分析，设置阈值筛选潜在异常采样波形；此处阈值的设置为粗略选取，优选地，满足多个采样数据片段作为第一候选异常波，进入下一个阶段的识别。
[0097]
波形匹配:对初选后的片段进一步精选，基于预先存储在波形库中的异常波形态对初选后对波形进行匹配，通过调整采样频率实现窗口大小的切换，进一步地，选取200ms到100ms之间的异常波识别。在一实施例中优选采用匹配度大于80％的脑活动采样数据作为异常活动的阈值。此步骤同时对采样异常波进行标记，同时对异常波出现前的相同时间尺度的非异常波进行二次标记作为噪声估计模块的输入。
[0098]
噪声估计模块：用于构建非异常波的第一采样数据生成模型，采用变分贝叶斯方法对第一采样数据生成模型进行求解，获取数据采集的背景干扰噪声，即采用基于变分贝叶斯的方法对非异常脑电分析估计背景干扰噪声统计分布信息。
[0099]
非异常波也即平稳波。通过对平稳波进行变分贝叶斯估计求解概率分布信息，作为背景干扰噪声的协方差矩阵。具体过程如下，假设非异常波为y
con
，变分贝叶斯因子分析假设采样数据表示为：
[0100]
y
con
＝fu+e
ꢀꢀ
(5)
[0101]
这里f是未知的m
×
k的加权矩阵，u＝[u1,u2,...,u
k
]
t
是未知的k
×
t
k
的因子参数，t
k
为非异常波中的采样点数目，设k＜＜m，即k为预设的小于m的任意正整数，e为非异常波分解后残余信号。因子参数的先验分布设为标准高斯分布如下：
[0102][0103]
同时传感器噪声的分布如下：
[0104][0105]
λ＝diag(λ1,λ2,...,λ
m
)是对角矩阵，λ
m
为第m个通道e的方差，λ为e的先验方差。加权矩阵同样满足高斯分布如下：
[0106][0107]
α
k
为加权矩阵f的先验方差，1≤k≤k。
[0108]
基于变分贝叶斯估计可得未知参数的更新规则如下：
[0109][0110][0111]
ψ＝r
uu
+α
ꢀꢀ
(11)
[0112][0113][0114]
其中，为m
×
k的加权矩阵的所有元素的均值，为k
×
t
k
的因子参数矩阵中所有元素的后验概率分布均值，r
yy
为采样数据协方差矩阵，r
uu
为估计因子协方差矩阵，r
uy
为估计因子和采样数据相乘，α＝diag(α1,α2...,α
k
)。迭代循环公式9至13直到和的值不在改变，根据此时和的值可得降维的干扰和背景噪声的协方差矩阵信息σ
e
描述为：
[0115][0116]
公式14即为基于非异常脑活动采样数据估计的背景干扰信息，其用于脑源活动重构模块脑异常活动的病灶区输入。
[0117]
脑源活动重构模块用于构建异常波的第二采样数据生成模型，背景干扰噪声是第二采样数据生成模型中的一个参数，采用经验贝叶斯方法，对第二采样数据生成模型进行求解，获取脑源活动时间序列。
[0118]
优选的，采用经验贝叶斯方法对第二采样数据生成模型进行求解包括步骤：设脑电源活动时间序列s在空间和时间上相互独立且服从高斯分布，采用经验贝叶斯方法将第二采样数据生成模型转换为第二采样数据概率生成模型；构建基于凸函数平面边界的代价函数，通过对代价函数进行优化来对第二采样数据概率生成模型进行求解。
[0119]
下面具体说明脑电源活动重构模块的优选实现方式。
[0120]
异常波的第二采样数据生成模型构建：
[0121]
设y
act
＝[y(t1),y(t2),...,y(t
k
)]为识别的异常脑电采样数据，t
k
为异常波的采样点数目，和非异常波的采样点数目相同。y(t
k
)∈r
m*1
为t
k
时刻的脑电采样数据，m为通道数目。脑源活动空间包含n个源活动点，位置i的单位源活动产生脑外采样强度为l
n
＝[l1,l2,...,l
m
]，l
n
是第n个源活动的导向矩阵，为已知参数，可通过源重构正问题求解，也可通过开源软件fieldtrip，nutmeg，spm等生成。异常波的第二采样数据生成模型如下：
[0122]
y
act
＝ls+ε
ꢀꢀ
(15)
[0123]
s＝[s1,s2,...,s
n
]
t
为源活动的时间序列，s
n
为位置n的源活动时间序列，ε为噪声与干扰数据，其统计特征信息由噪声估计模块公式14求解。设脑源活动在空间和时间上相互独立且服从高斯分布如下：
[0124]
[0125]
为位置n的先验高斯分布的方差γ＝diag(γ1,
…
,γ
n
)，则条件概率分布p(y|s)为，
[0126][0127]
σ
e
是s4中估计的噪声与干扰的协方差矩阵。基于经验贝叶斯模型可得：
[0128][0129]
采样数据p(y
act
)的计算困难且在经验贝叶斯估计过程中仅起规约化作用，对估计结果并不造成影响。由于先验分布p(s|γ)与条件分布p(y
act
|s)均为高斯分布，后验概率分布p(s|y
act
)为
[0130][0131]
其中，是源活动的后验概率分布的均值，γ
‑1为后验分布等方差。
[0132]
公式16，17，18和19构成第二采样数据生成模型的第二采样数据概率生成模型，下面求解模型参数。
[0133]
s5.2：基于凸函数平面边界的代价函数优化
[0134]
超参数γ
‑1的求解通过最大化边缘分布p(y
act
|γ)，边缘分布对数模型logp(y
act
|γ)为：
[0135][0136]
为异常波采样数据协方差矩阵模型，公式21可通过期望最大化算法(em算法)优化，但log|∑
y
|的收敛速度极慢，由于log|∑
y
|是γ
‑1的凸函数，构建基于γ
‑1的log|∑
y
|边界函数，
[0137][0138]
优化求解
‑
log|∑
y
|的最大值即边界代价函数的最小值求解，如下所示，
[0139][0140]
λ
n
与λ0为辅助变量，tr(*)是计算矩阵的迹，基于凸函数边界的代价函数始终满足：
[0141][0142]
logp(y
act
|γ)的最大值优化即为代价函数的最小值求解。
[0143]
s5.3:复杂病灶区重构
[0144]
通过以上分析可得，脑电源的求解首先假设源活动服从相互独立的高斯分布，采用经验贝叶斯估计，求解源活动的值将公式16、17和19带入公式18中求解后验分布p(s|y
act
)的均值与方差如下：
[0145][0146][0147]
为对角矩阵。公式20中的均值即待求解源活动时间序列，其求解需超参数γ已知。超参数γ需通过求解边缘分布的最大值来求取，构建基于凸函数边界的代价函数模型，超参数γ求解即转化为代价函数的优化。对函数求超参数γ
n
和辅助参数λ
n
倒数，可得：
[0148][0149]
使公式25为零，超参数γ
n
和辅助参数λ
n
的公式如下：
[0150]
[0151]
通过循环迭代公式24和公式26直至代价函数公式20收敛，即为求解脑电源的活动时间序列。
[0152]
成像模块，用于根据脑电源活动重构模块获取的脑电源活动时间序列进行成像脑电源活动输出与成像。
[0153]
获取脑电源活动时间序列后，需要对结果进行展示，进一步地，展示结果基于现存模型参数，如脑空间分割后源位置以及脑结构信息，本专利采用mri模版，渲染源活动与通用脑模型上，运行结果如图2所示。
[0154]
必须说明的是，上述任一实施例中，方法并不必然按照序号顺序依次执行，只要从执行逻辑中不能推定必然按某一顺序执行，则意味着可以以其他任何可能的顺序执行。
[0155]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。