基于克里金插值的EEG二维成像最优化分析

蔡 靖， 胡纪锋， 刘锋华， 蒙坚发

（吉林大学仪器科学与电气工程学院，长春 130061）

0 引 言

脑电信号（Electroencephalography EEG）是脑神经细胞群产生的微弱的非平稳伪随机生物电信号，含有丰富的大脑活动状态信息。随着近年信息技术的发展，运用脑电信号进行脑电成像分析并应用于临床研究越来越成为研究的重点，脑电成像可以直观显示大脑活动状态，对于脑血管疾病、脑积水、脑肿瘤等病灶的定位［1-3］以及癫痫、中风等疾病发作时脑电成像特点［4-5］的临床医学研究有很重要的意义。

Lee等［6］开发了一台带有IBM PC AT的电脑脑电图成像系统，采用反距离加权插值的方法，使每个插值点的权重与其最近的4个点的距离成反比，实现了脑电成像。Yang等［7］借助混合效应多项式回归模型的统计分析方法，成功实现了用高密度脑电图来描述次最大肌肉收缩时脑电成像中相应动态源的变化。Paul等［8］对于癫痫病症进行了研究，结果表明，在大多数情况下使用k-最近邻点插值法可以有效提高灵敏度，并且误差也较小。Taran等［9］使用IMF1的径向基核函数对不同运动图像任务的脑电信号进行分类，该方法与现有的方法相比，在分类准确度，灵敏性等方面展现了更好的性能。

在以上研究的基础上，本文对基于高斯模型变异函数普通克里金插值算法进行优化分析，在运用数据拟合理论的基础上建立了变异函数各参量之间的函数关系，并结合模型的均方误差给出了优化的普通克里金变异函数模型取值规范。

1 脑电成像算法分析

1.1 普通克里金插值

克里金插值是Matheron首先提出［10-12］的一种对有限区域内部的区域化变量进行无偏最优估计的空间插值算法。

克里金插值法中的普通克里金法具有下列约束条件：

式中：x为位置坐标；h为距离；Z（x＋h）、Z（x）分别为x＋h与x位置处的实测值；γ（h）为距离h对应的半变异函数。

对于需要插值获得的单个待测样本点x0，由n个样本点的实测值Z（xi）可得到待测点Z*（x0）的估计值：

式中，λ为权重系数。则估计值的方差为：

式中：hi为x0到xi的距离；L为拉格朗日乘数。

对于多个待测点，通常采用估计值方差的平均值来衡量插值效果。估计值方差的平均值为：

式中：N为待测点数；n为已知实际测量点的个数。估计值方差的平均值作为脑电插值成像精度的评判标准，当估计值方差的平均值＝0时即认为是最优的脑电成像。

1.2 半变异函数模型与优化关系确定

式（5）中，半变异函数γ（h）表征了插值区域的空间变异结构，一个插值区域的半变异函数定义为：

实际运用中为了简化计算，通常选用下列半变异函数模型。

高斯模型：

指数模型：

球型模型：

式中：C0为块金值，是由样本误差和短距离的变异性引起的半变异函数的偏差；ɑ为相关尺度，表示当样本之间的距离大于等于此距离时，各样本相互独立。C0＋C1＝var［Z（x）］，称为基台值。

可以得出C0与ɑ的取值决定了方差均值的大小。经过比较，本文选用高斯模型，在高斯模型下使

的C0与ɑ的值即为高斯模型的最优参数。

1.3 相关尺度a与块金值最优值获取

对于n个样本点，以任意两个样本点之间的距离h作为横坐标，由式（6）得出的C（h）值作为纵坐标，进行二次曲线拟合，即可得到取样区域的近似半变异函数，拟合结果如图1所示。

图1 二次拟合曲线

拟合多项式如下：

联立求解式（7）、（11）式可得ɑ与C0的关系：

对于不同的距离h，ɑ和C0有相似的关系，取所有h的平均值h0代入上式，得到最优的ɑ和C0的关系（见图2）。

图2 ɑ-C0关系曲线

联立式（10）、（12）得：

即可求出C0与ɑ的最优值。

由式（13）得出的ɑ与C0参数进行优化普通克里金插值，并用此插值结果进行最优化脑电成像。

2 脑电信号采集系统

如图3所示为脑电信号采集系统结构图，主要由采集电极、前置滤波放大、模数转换、微处理器、蓝牙通信和上位机信号处理等单元组成。

图3 系统结构图

采集电极安放位置参考国际10-20系统电极放置法，8个采集电极Fp1、Fp2、Fz、Cz、P3、P4、O1、O2同时采集脑电信号，信号经由前置滤波放大单元滤去高频段的噪声，并将信号放大。模数转换单元将8通道的模拟脑电信号转换为数字信号，微处理器单元负责接收此数字信号并将其通过蓝牙通信传输到上位机中进行信号的处理。

3 实验及结果分析

3.1 原始脑电信号

由采集系统采集Fp1、Fp2、Fz、Cz、P3、P4、O1、O2点的相对于耳垂的原始电位值脑电信号如图4所示。

3.2 测量电极坐标

为进行脑电2D成像，以国际10-20电极安放标准中的CZ为极点，以CZ点到眉间的弧线为极轴建立极坐标系，极径为电极安放点到原点的大脑弧面距离，极角为电极安放点与极点相连线段在极点所在平面的投影与极轴的夹角。

将人的大脑近似为一个球面，假设某两个电极Fp1与Fz的球面坐标分别为r1、θ1、ϕ1与r1、θ2、ϕ2，那么两者之间的球面距离，即弧长l可以用如下方程组求解：

建立平面直角坐标系，将极坐标点转换为直角坐标点并平移到第1象限内。即可得到电极安放位置坐标如图5所示。

图5 测量电极坐标图

3.3 半变异函数模型选择

普通克里金插值中变异函数模型的选取对插值的无偏性具有重要意义。图6～8分别为变异函数模型为高斯模型、指数模型、球形模型时脑电成像结果及对应的误差分布图。模型待测点估计值方差分布情况如图9所示，不同半变异函数模型插值点误差曲线如图10所示。

图6 高斯模型脑电成像及误差图

图7 指数模型脑电成像及误差图

图8 球形模型脑电成像及误差图

由图可见，高斯模型每个待测点的估计值误差均在零附近，且平均误差明显低于其他2种模型，因此本文进行普通克里金插值优化分析时选用的半变异函数模型为高斯模型。

3.4 最优值确定

当相关尺度ɑ与块金值C0取值不同时，由式（4）可得平均方差变化曲线如图11所示。

图9 方差散点图

图10 不同半变异函数模型插值点误差曲线

图11 平均方差

由式（13）可得参数最优解。图12中曲面与曲线交点即反映了当平均误差趋近于0时的相关尺度ɑ与块金值C0的最优值。

图12 交点图

3.5 结果分析

图13为在正常静默状态下，相关尺度ɑ与块金值C0取位于最优解附近的值时测试者的2D脑电成像情况及对应的误差分布图。

图13 脑电成像及误差图

均方差值如表1所示。

以上实验的基础上采用逐次逼近的方法，可得均方差趋近于0时的最优的ɑ和C0值。

图14为身体健康的成年男性测试者在ɑ和C0取最优值即均方差趋近于0时所绘制的最优2D脑电成像与误差分布图：

在均方差趋近于零的约束条件下，对高斯模型半变异函数的参数ɑ和C0的函数关系进行最优解求取，实现了基于高斯模型半变异函数的普通克里金插值的优化分析。实验结果表明，采用逐步逼近的方法平均方差趋近于零。

表1 均方差值表

图14 最优脑电成像及误差图

4 结 语

本文提出了一种基于克里金插值的EEG二维成像最优化分析方法。通过仿真比较，使用高斯半变异函数模型的EEG二维成像的方差均值最小。通过拟合高斯半变异函数中的协方差函数，得到相关尺度ɑ与块金值C0的函数关系。方差均值是误差评估标准，当方差均值接近零时，ɑ与C0将取得最优值。通过将ɑ与C0的最优值应用于EEG二维成像，实现了最优EEG二维成像。该方法提高了EEG二维成像的准确性，使得分析EEG信号并从EEG中获取生理和疾病信息变得更加容易。