具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明实施例的基于贝叶斯的儿童脑活动异常病灶区成像系统包括信号采集处理模块、脑活动异常识别模块、噪声估计模块、脑电源活动重构模块和成像模块。本系统可适用于各类通过脑电磁活动异常病灶区成像来识别疾病的场景中，尤其适用于儿童自闭症。

信号处理模块，用于获取电磁采样数据，并对电磁采样数据进行预处理。

下面说明信号处理模块的优选实现方式。

信号采集：

输入脑电磁采样数据，Y_raw＝[y_raw(1),y_raw(2),...,y_raw(T)],T是原始信号采样点数目，y_raw(t)＝[y₁(t),y₂(t),…,y_M(t)]^T，这里M为采集点的通道数目，T表示转置矩阵。需要对输入数据做初步筛查，保证采样数据的有效性以及完整性，如数据格式转换以及参数确认，包括采样数据信息等。

信号预处理：

信号采集处理模块包括移除采样偏移模块、去干扰模块和滤波模块，对输入数据做预处理，包括移除采样偏移，去燥，滤波等预处理。

移除采样偏移的公式如下：

去干扰需要根据采样数据特征来进行。

优选地，基于独立成分分析(ICA)或空间投影等方式，如基于ICA的去干扰(包括眼电等)，公式如下：

其中A＝[a₁,a₂,...,a_N]为分解加权矩阵，X＝[x₁,x₂,...,x_N]^T为独立成分，设x_m包含显著干扰或者噪声的特征，则去燥后采样数据Y_new如下：

优选地，利用基于投影的去干扰方法对经过基于ICA的去干扰处理后的信号进一步进行处理。基于投影的去干扰方法则需要找出投影矩阵，其可以基于脑结构空间信息等特征，假设投影矩阵为W，其去干扰公式为：

Y_clean＝W^TY_new (4)

滤波:对采样数据的滤波去除非脑活动频段的采样数据，主要应用基于傅立叶或者小波变换的时频滤波，保留感兴趣频段数据。通过预处理后的数据保证了后续异常识别精准性。

脑活动异常识别模块，用于从预处理后的脑电磁采样数据中识别出异常波和非异常波。

优选的，脑活动异常识别是基于采样脑电数据的振幅以及局部波形来识别。

根据振幅初选：对采样脑电数据的振幅进行统计分析，设置阈值筛选潜在异常采样波形；此处阈值的设置为粗略选取，优选地，满足多个采样数据片段作为第一候选异常波，进入下一个阶段的识别。

波形匹配:对初选后的片段进一步精选，基于预先存储在波形库中的异常波形态对初选后对波形进行匹配，通过调整采样频率实现窗口大小的切换，进一步地，选取200ms到100ms之间的异常波识别。在一实施例中优选采用匹配度大于80％的脑活动采样数据作为异常活动的阈值。此步骤同时对采样异常波进行标记，同时对异常波出现前的相同时间尺度的非异常波进行二次标记作为噪声估计模块的输入。

噪声估计模块：用于构建非异常波的第一采样数据生成模型，采用变分贝叶斯方法对第一采样数据生成模型进行求解，获取数据采集的背景干扰噪声，即采用基于变分贝叶斯的方法对非异常脑电分析估计背景干扰噪声统计分布信息。

非异常波也即平稳波。通过对平稳波进行变分贝叶斯估计求解概率分布信息，作为背景干扰噪声的协方差矩阵。具体过程如下，假设非异常波为Y_con，变分贝叶斯因子分析假设采样数据表示为：

Y_con＝FU+E (5)

这里F是未知的M×K的加权矩阵，U＝[u₁,u₂,...,u_K]^T是未知的K×t_K的因子参数，t_K为非异常波中的采样点数目，设K＜＜M，即K为预设的小于M的任意正整数，E为非异常波分解后残余信号。因子参数的先验分布设为标准高斯分布如下：

同时传感器噪声的分布如下：

Λ＝diag(λ₁,λ₂,...,λ_M)是对角矩阵，λ_m为第m个通道E的方差，Λ为E的先验方差。加权矩阵同样满足高斯分布如下：

α_k为加权矩阵F的先验方差，1≤k≤K。

基于变分贝叶斯估计可得未知参数的更新规则如下：

Ψ＝R_uu+α (11)

其中，为M×K的加权矩阵的所有元素的均值，为K×t_K的因子参数矩阵中所有元素的后验概率分布均值，R_yy为采样数据协方差矩阵，R_uu为估计因子协方差矩阵，R_uy为估计因子和采样数据相乘，α＝diag(α₁,α₂...,α_k)。迭代循环公式9至13直到和的值不在改变，根据此时和的值可得降维的干扰和背景噪声的协方差矩阵信息Σ_E描述为：

公式14即为基于非异常脑活动采样数据估计的背景干扰信息，其用于脑源活动重构模块脑异常活动的病灶区输入。

脑源活动重构模块用于构建异常波的第二采样数据生成模型，背景干扰噪声是第二采样数据生成模型中的一个参数，采用经验贝叶斯方法，对第二采样数据生成模型进行求解，获取脑源活动时间序列。

优选的，采用经验贝叶斯方法对第二采样数据生成模型进行求解包括步骤：设脑电源活动时间序列S在空间和时间上相互独立且服从高斯分布，采用经验贝叶斯方法将第二采样数据生成模型转换为第二采样数据概率生成模型；构建基于凸函数平面边界的代价函数，通过对代价函数进行优化来对第二采样数据概率生成模型进行求解。

下面具体说明脑电源活动重构模块的优选实现方式。

异常波的第二采样数据生成模型构建：

设Y_act＝[y(t₁),y(t₂),...,y(t_K)]为识别的异常脑电采样数据，t_K为异常波的采样点数目，和非异常波的采样点数目相同。y(t_k)∈R^M*1为t_k时刻的脑电采样数据，M为通道数目。脑源活动空间包含N个源活动点，位置i的单位源活动产生脑外采样强度为L_n＝[l₁,l₂,...,l_M]，L_n是第n个源活动的导向矩阵，为已知参数，可通过源重构正问题求解，也可通过开源软件fieldtrip，nutmeg，spm等生成。异常波的第二采样数据生成模型如下：

Y_act＝LS+ε (15)

S＝[S₁,S₂,...,S_N]^T为源活动的时间序列，S_n为位置n的源活动时间序列，ε为噪声与干扰数据，其统计特征信息由噪声估计模块公式14求解。设脑源活动在空间和时间上相互独立且服从高斯分布如下：

为位置n的先验高斯分布的方差γ＝diag(γ₁,…,γ_n)，则条件概率分布p(Y|S)为，

Σ_E是S4中估计的噪声与干扰的协方差矩阵。基于经验贝叶斯模型可得：

采样数据p(Y_act)的计算困难且在经验贝叶斯估计过程中仅起规约化作用，对估计结果并不造成影响。由于先验分布p(S|γ)与条件分布p(Y_act|S)均为高斯分布，后验概率分布p(S|Y_act)为

其中，是源活动的后验概率分布的均值，Γ^-1为后验分布等方差。

公式16，17，18和19构成第二采样数据生成模型的第二采样数据概率生成模型，下面求解模型参数。

S5.2：基于凸函数平面边界的代价函数优化

超参数Γ^-1的求解通过最大化边缘分布p(Y_act|γ)，边缘分布对数模型logp(Y_act|γ)为：

为异常波采样数据协方差矩阵模型，公式21可通过期望最大化算法(EM算法)优化，但log|∑_y|的收敛速度极慢，由于log|∑_y|是γ^-1的凸函数，构建基于γ^-1的log|∑_y|边界函数，

优化求解-log|∑_y|的最大值即边界代价函数的最小值求解，如下所示，

Λ_n与Λ₀为辅助变量，tr(*)是计算矩阵的迹，基于凸函数边界的代价函数始终满足：

logp(Y_act|γ)的最大值优化即为代价函数的最小值求解。

S5.3:复杂病灶区重构

通过以上分析可得，脑电源的求解首先假设源活动服从相互独立的高斯分布，采用经验贝叶斯估计，求解源活动的值将公式16、17和19带入公式18中求解后验分布p(S|Y_act)的均值与方差如下：

为对角矩阵。公式20中的均值即待求解源活动时间序列，其求解需超参数γ已知。超参数γ需通过求解边缘分布的最大值来求取，构建基于凸函数边界的代价函数模型，超参数γ求解即转化为代价函数的优化。对函数求超参数γ_n和辅助参数Λ_n倒数，可得：

使公式25为零，超参数γ_n和辅助参数Λ_n的公式如下：

通过循环迭代公式24和公式26直至代价函数公式20收敛，即为求解脑电源的活动时间序列。

成像模块，用于根据脑电源活动重构模块获取的脑电源活动时间序列进行成像脑电源活动输出与成像。

获取脑电源活动时间序列后，需要对结果进行展示，进一步地，展示结果基于现存模型参数，如脑空间分割后源位置以及脑结构信息，本专利采用MRI模版，渲染源活动与通用脑模型上，运行结果如图2所示。

必须说明的是，上述任一实施例中，方法并不必然按照序号顺序依次执行，只要从执行逻辑中不能推定必然按某一顺序执行，则意味着可以以其他任何可能的顺序执行。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。