CN1333003A - 大脑头皮电位信息图形的生成方法 - Google Patents

Abstract

本发明公开了一种利用计算机系统实现大脑头皮电位图形的生成方法。它涉及一种电位信息图形的生成方法。其实现步骤是从头部不同导联电极上采集的脑电原始数据出发,按照非线性动力学理论,运用相空间重建技术构成不同导联之间的头皮电位信息传输矩阵,然后将各导联的信息传输量和时间序列复杂度数据绘制成头皮拓扑分布图－脑电信息图(BIM),以一种全新的角度连续反映脑电活动的时空特征,观察人脑功能状态的变化过程。本发明可运用于诸如癫痫病灶定位和诱发响应等的研究中的头皮脑电诊断,还可用于大脑处于不同生理状态下或老年痴呆病患者的EEG信号分析。

Description

大脑头皮电位信息图形的生成方法

                        技术领域

本发明涉及一种图形的生成方法，尤其涉及一种电位信息图形的生成方法。

                        背景技术

目前临床使用的脑电地形图(Brain Electric Activity Mapping——BEAM)是基于头皮脑电计算机分析的一种诊断技术。其原理是从头部不同导联电极上采集脑电信号，经滤除各种干扰和A/D变换后作快速富氏变换得到不同频段的信号功率谱，再用二维插值方法显示为等值线地形图。BEAM不仅用于临床诊断脑部器质性病变，而且还用于观察脑功能的变化，已成为诊断大脑疾病的常规手段。但由于其信号分析的理论基础是从时域到频域富氏变换的分频段功率谱，所含数据量相对较少，数据选取存在一定主观性，在诸如癫痫病灶定位和诱发响应等研究中难以全面反映脑电动态变化过程。

                          发明内容

本发明为克服现有技术的不足，提供一种大脑头皮电位信息传输的图形显示方法，它从头部不同导联电极上采集的脑电原始数据出发，按照非线性动力学理论，运用相空间重建技术构成不同导联之间的头皮电位信息传输矩阵，然后将各导联的信息传输量和时间序列复杂度数据绘制成头皮拓扑分布图——脑电信息图(Brain Information Mapping——BIM)，以探索一种新的观察大脑功能活动的诊断方法，该信息图以时空演化方式显示头皮上不同导联之间信息传输情况和信息传输时间序列复杂度等非线性动力学参数变化，可以直观地反映脑电信息传输分布模式在不同时相中的变化进程，以一种全新的角度连续反映脑电活动的时空特征，观察人脑功能状态的变化过程。

为达到上述技术效果，本发明采用下列方法：(1)从头部不同导联电极上采集脑电信号波形；(2)将上述信号作A/D变换并滤除各种干扰后将其存储于存储器中；(3)选择脑电中无显著伪差、无脑电图机关机干扰的数据段，存储于存储器中；(4)将编制的计算机程序安装到计算机系统后，启动该计算机程序，利用该计算机程序对存储器中的每个导联的脑电信号重构一个三维嵌入空间的向量集；(5)利用上述计算机程序建立信息传输矩阵，分析上述向量集在相空间的相互信息传输情况；(6)利用上述计算机程序，将时间序列信号粗粒化为(0，1)序列并用复杂度Kc描述事物随机性的复杂度，用复杂度C1和C2反映事物结构性的复杂度；(7)在计算机程序中，采用球面二维插值方法，将各导联之间的信息传输时间序列复杂度数据绘制成头皮拓扑分布图——脑电信息图(BIM)；(8)将得到的脑电信息图(BIM)再输出并存盘。

所述的三维嵌入空间的向量集的构成采用下列方法：

对于一个以等时间隔T_s采样的脑电信号V(k)(k＝0，1，2，…，N-1)，根据时间弛豫法，取延迟时间T＝1/512秒，通过延时构成M维向量X(k)＝{V(k)，V(k＋T)，…，V(k＋(M－1)T)}，取M＝3，每个导联的脑电信号V(k)构成一个三维嵌入空间的向量集，并采用下列公式表达： $\left[\begin{array}{ccc}X\left(1\right),& X\left(2\right),& X\left(3\right)\\ X\left(2\right),& X\left(3\right),& X\left(4\right)\\ \…\…& \…\…& \…\…\\ X\left(n\right),& X(n+1),& X(n+2)\end{array}\right]$ 将此三维相空间的每一坐标轴作六等分，得到信息熵H(S)： $H\left(S\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_s\left(S_i\right)\mathrm{lo}{g}_2\[P_s\left(S_i\right)\](N=6^3)$

以S和Q分别代表两个导联的脑电信号所构成三维相空间的相互信息传输格点序列；若把两个三维空间联合起来构成一个六维相空间，同样将每一坐标轴作六等分，再求信息熵H(S，Q)得： $H(S,Q)=-\underset{i=1}{\mathrm{\Σ}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\stackrel{N}{P_{\mathrm{SQ}}\left(S_i{Q}_j\right)}\mathrm{lo}{g}_2\[P_{\mathrm{SQ}}\left(S_i{Q}_j\right)\](N=6^6)$

据此计算两个导联的相互信息传输量MI(S，Q)：

    MI(S，Q)＝H(S)＋H(Q)－H(S，Q)；

所述的信息传输矩阵采用下列方法：

用时间延迟方法构成两个导联的相互信息传输时间序列SQ(τ)：

    SQ(τ)＝MI[S(t)Q(t＋τ)]    τ＝T，2T，…，KT

SQ(τ)代表从S导联到Q导联的信息传输；

而相互信息传输时间序列QS(τ)为：

    QS(τ)＝MI[Q(t)S(t＋τ)]    τ＝T，2T，…，KT

QS(τ)则代表从Q导联到S导联的信息传输，序列长为512点；矩阵中行代表输出信息，列代表接收信息；

所述的复杂度计算，采用下列方法：

对所有(0，1)序列，定义b(n)

       b(n)＝(log₂n)/n

n为数据序列号。再定义C(n)

    limC(n)＝lim(log₂n)/n

    n→∝    n→∝

定义复杂度Kc为：

    Kc＝C(n)/b(n)

对于符号(0，1)序列{S_i}^N _i＝1，定义复杂度C1(n)为：

    C1＝lim[log₂N_s(n)]/n

    n→∝

对于时间符号(0，1)序列{S_i}^N _i＝1，定义复杂度C2(n)为：

    C2＝lim[log₂N_f(n)]/n

    n→∝

C1反映了时间序列在相空间运动轨道的随机性结构复杂度；而C2则反映了时间序列征相空间运动轨道的限制性结构复杂度。

本发明与现有技术相比的有益效果是：由本发明获得的脑电信息图(BIM)，能直观地反映脑电信息传输分布模式在不同时相中的变化进程，能表达长程脑电模式变化，全面反映脑电动态变化过程，便于比较大脑左右两侧图象变化特点以对比了解两侧功能的变化；此外，也能将BIM用于大脑处于不同生理状态下或老年痴呆病患者的EEG信号分析。

                            附图说明

图1是本发明工作流程示意图；

图2(a)是本发明实施方式中，当癫痫发作前的脑电信号波形；

图2(b)是本发明实施方式中，当癫痫小发作的脑电信号波形；

图3(a)是图2(a)的8×8信息传输矩阵；

图3(b)是图2(b)的8×8信息传输矩阵；

图4(a)是图3(a)的信息输出量的脑电信息图；

图4(b)是图3(b)的信息输出量的脑电信息图；

图5(a)是图2(a)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度Kc的脑电信息图；

图5(b)是图2(b)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度Kc的脑电信息图；

图5(c)是图2(a)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度C1的脑电信息图；

图5(d)是图2(b)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度C1的脑电信息图；

图5(e)是图2(b)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度C2的脑电信息图；

图5(f)是图2(b)当痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度C2的脑电信息图。

                        具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述：

图1示出本发明的工作流程，使用脑电采集分析系统，按照国际标准导联10～20系统在受试者头皮上安放F_p1(左额)、F_p2(右额)、C₃(左顶)、C₄(右顶)、T₃(左颞)、T₄(右颞)、O₁(左枕)、O₂(右枕)8导电极，参考电极置于左右耳垂，同时采集各导联EEG信号波形；要求受试者实验前三天内不服药，实验时未发病，清醒闭目，静坐在安静的实验室内。原始脑电信号采集时间为5分钟，A/D采样频率为512Hz，记录患者作深呼吸诱导时的脑电数据。选择脑电中无显著伪差、无脑电图机关机干扰的数据段，长度为2562点，时间为5秒，进行相空间重建和构成信息传输矩阵及计算信息传输量和时间序列复杂度。

对于一个以等时间隔T_s采样的脑电信号V(k)(k＝0，1，2，…，N-1)，可根据时间弛豫法，选取延迟时间T＝1/512秒，通过延时构成M维向量X(k)＝{V(k)，V(k＋T)，…，V(k＋(M－1)T)}，取M＝3；这样对于每个导联的脑电信号V(k)可构成一个三维嵌入空间的向量集： $\left[\begin{array}{ccc}X\left(1\right),& X\left(2\right),& X\left(3\right)\\ X\left(2\right),& X\left(3\right),& X\left(4\right),\\ \…\…& \…\…& \…\…\\ X\left(n\right),& X(n+1),& X(n+2)\end{array}\right]$

将此三维相空间的每一坐标轴作六等分，则其相空间的格点数为6³。以S代表相空间的格点序列S＝S₁S₂……S_N，N＝6³。设向量落在第i个格点的概率为P_S(S_i)，可得到信息熵H(S)： $H\underset{i=1}{\overset{N}{\left(S\right)=\mathrm{\Σ}{P}_s}}\left(S_1\right)\mathrm{lo}{g}_2\[P_s\left(S_i\right)\](N=6^3)----\left(2\right)$

以S和Q分别代表两个导联的脑电信号所构成三维相空间的格点序列。若把两个三维空间联合起来构成一个六维相空间，同样将每一坐标轴作六等分，则其相空间的格点6⁶，其对应的概率分布为P_SQ(S_iQ_j)，再求信息熵H(S，Q)得： $H(S,Q)=-\mathrm{\Σ}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{P_{\mathrm{SQ}}\left(S_i{Q}_j\right)}}\mathrm{lo}{g}_2\[P_{\mathrm{SQ}}\left(S_i{Q}_j\right)\](N=6^6)----\left(3\right)$

据此可以计算两个导联的相互信息传输量MI(S，Q)：

MI(S，Q)＝H(S)＋H(Q)－H(S，Q)                 (4)

同样可以用时间延迟方法构成两个导联的相互信息传输时间序列SQ(τ)：

SQ(τ)＝MI[S(t)Q(t＋τ)] τ＝T，2T，…，KT    (5)

SQ(τ)代表从S导联到Q导联的信息传输。而相互信息传输时间序列QS(τ)为：

QS(τ)＝MI[Q(t)S(t＋τ)] τ＝T，2T，…，KT    (6)

QS(τ)则代表从Q导联到S导联的信息传输。序列长为512点。对于8个导联的脑电数据来说，可以得到8×8列的相互信息传输时间序列，也即组成8×8的信息传输矩阵。矩阵中行(横向)代表输出信息，列(纵向)代表接收信息。利用信息传输矩阵可对各导联脑电信号所构成三维向量集在相空间的相互信息传输情况进行分析。

本发明使用复杂度Kc描述事物随机性的复杂度；使用复杂度C1和C2，以反映事物结构性的复杂度。根据事物的复杂性可用描写该事物所用最短计算机语言的长度来表达的复杂性测度概念，可先将时间序列信号粗粒化为(0，1)序列。粗粒化前需计算序列信号平均值，粗粒化时，当信号值大于平均值时取为1，信号值小于平均值时取为0，构成(0，1)序列。

对所有(0，1)序列，定义b(n)

    b(n)＝(log₂n)/n                               (7)

n为数据序列号。再定义C(n)

    limC(n)＝lim(log₂n)/n    (8)

    n→∝    n→∝

复杂度Kc定义为：

    Kc＝C(n)/b(n)             (9)

对于符号(0，1)序列{S_i}^N _i＝1，其N足够大，则n个字节(S₁，S₂，…，S_n)的可能排列有2ⁿ个；令N_a(n)是长度为n个字节(S₁，S₂，…，S_n)在{S_i}^N _i＝1中允许出现的序列数，定义复杂度C1(n)为：

    C1＝Iim[log₂N_a(n)]/n    (10)

    n→∝

对于时间符号(0，1)序列{S_i}^N _i＝1中，令N_f(n)是长度为n-1个字节(S₁，S₂，…，S_n-1)在{S₁}^N _i=1中出现，但不出现n个字节(S₁，S₂，…，S_n)的禁止序列数，定义复杂度C2(n)为：

    C2＝lim[log₂N_f(n)]/n    (11)

    n→∝

C1反映了时间序列在相空间运动轨道的随机性结构复杂度；而C2则反映了时间序列在相空间运动轨道的限制性结构复杂度。

虽然可以直接利用脑电信息传输矩阵中各阵元时间序列所反映的信息量或复杂度来分析各导联脑电的相互信息传输情况，但借助传统脑电地形图(BEAM)的二维等值线地形图示方法更能形象直观地反映大脑皮层信息传输动态变化，便于比较大脑左右两侧图象变化特点以对比了解两侧功能的变化。本发明采用球面二维插值方法，将各导联之间的信息传输量和时间序列复杂度数据绘制成头皮拓扑分布图——脑电信息图(BIM)。

图2(a)和图2(b)给出在局限性癫痫病患者8个导联处测得的典型脑电信号波形。其中图2(a)为癫痫发作前的EEG信号，图2(b)为癫痫小发作时的EE6信号。可以看出分别对应于C₃(左顶)和T₃(左颞)的第3、5导联是痫性病灶区。

图3(a)和图3(b)是相应的8×8信息传输矩阵。其中图3(a)为癫痫发作前的信息传输矩阵，图3(b)为癫痫小发作时的信息传输矩阵。从图3(a)可以看出癫痫发作前各导联之间的信息传输都很活跃，每个导联既输出又接收大量信息。而在癫痫发作时如图3(b)所示，唯痫性的第3、5导联信息传输较活跃，其它导联信息传输则比癫痫发作前减少很多。特别是痫性第3、5导联之间的信息传输量很大。该病例从脑电图观察诊断为“左侧灶性放电”，信息传输矩阵的变化也表明，癫痫发作时病灶区神经细胞活动增加，一方面极大地激励了痫性导联处信息传输，另一方面其类周期性放电又强烈地抑制了其它部位之间的信息传输。可以看出，从信息传输时间序列(8×8矩阵)方式来观察局限性癫痫病发作前与发作时各导联之间信息传输量的变化，很难直观把握其变化特征。若以脑电信息图(BIM)的方式来显示，则不仅具有很好的直观性，而且有可能发现隐含在信息传输时间序列内部的变化规律。

图4(a)和图4(b)是相应于图3(a)和图3(b)信息传输矩阵中信息输出量的复杂度Kc脑电信息图(BIM)。其中图4(a)为癫痫发作前的BIM图，图4(b)为癫痫小发作时的BIM图。从图4(a)可以看出癫痫发作前大脑皮层各部分之间的信息传输有较好的对称性。而在癫痫发作时如图4(b)所示，这种对称性遭到较大的破坏，尤其在痫性导联的对侧，对称性破坏更为严重。这与上述信息传输矩阵各阵元的变化表象是相吻的，但这种对称性破坏的变化特征却难以在信息传输矩阵中反映出来。

根据本发明原理，还可以作出以某导联为中心与其余各导联信息传输时间序列复杂度为表征的脑电信息图(BIM)来比较癫痫病发作前与发作时大脑皮层信息传输的变化。其中较为有意义的是该导联向其它各导联信息输出时间序列复杂度的脑电信息图(BIM)。图5(a)～图5(f)分别为痫性导联C₃(左顶)向其它各导联信息输出时间序列复杂度Kc、C1、C2的脑电信息图(BIM)。图5(a)、图5(c)和图5(e)分别为癫痫发作前的BIM图，图5(b)、图5(d)和图5(f)分别为癫痫小发作时的BIM图。从各图可以看到在癫痫发作前，三种复杂度分布皆有某种程度的对称性，而癫痫发作时三种复杂度的对称性分布皆受到破坏。其中如图5(b)和图5(d)示出了复杂度Kc和C1在癫痫发作时痫性病灶区，即痫性导联C₃和T₃之间区域的复杂度值明显低于其它区域，而图5(e)示出此时复杂度C2在痫性病灶区，痫性导联的复杂度值却明显高于其它区域。即复杂度Kc和C1的变化趋势相似，而复杂度C2的变化趋势与前二者不尽相同。这一结果可理解为痫性导联处神经细胞因过度重复放电使其脑电的类周期活动增强，导致脑电信息传输时间序列的有序性增加，随机性降低，相应的复杂度Kc和C1减小；同时其相空间运动轨迹结构受到更多的限制，使其结构复杂度C2高于其它区域。癫痫发作前后脑电信息图(BIM)的变化再次说明：三种复杂度从不同方面刻化了时间序列的变化特点；复杂度C1、C2的引入，可比单独使用复杂度Kc更为全面地反映EEG变化的随机性和结构性，尤其复杂度C2的定义具有新内含。

Claims (4)

1.一种利用计算机系统实现大脑头皮电位信息图形的生成方法，包括如下步骤：

(1)从头部不同导联电极上采集脑电信号波形；

(2)将上述信号进行A/D变换并滤除各种干扰后，将其存储于存储器中；

(3)选择脑电中无显著伪差、无脑电图机关机干扰的数据段，存储于存储器中；

其特征在于它还采用以下步骤：

(4)将编制的计算机程序安装到计算机系统后，启动该计算机程序，利用该计算机程序对存储器中的每个导联的脑电信号重构一个三维嵌入空间的向量集；

(5)利用上述计算机程序建立信息传输矩阵，分析上述向量集在相空间的相互信息传输情况；

(6)利用上述计算机程序，将时间序列信号粗粒化为(0，1)序列，并用复杂度Kc描述事物随机性的复杂度，用复杂度C1和C2反映事物结构性的复杂度；

(7)在计算机程序中，采用球面二维插值方法，将各导联之间的信息传输量和时间序列复杂度数据绘制成头皮拓扑分布图——脑电信息图(BIM)；

(8)将得到的脑电信息图(BIM)再输出并存盘。

2.根据权利要求1所述的大脑头皮电位信息图形的生成方法，其特征在于所述的三维嵌入空间的向量集的构成采用下列方法：

对于一个以等时间隔T_s采样的脑电信号V(k)(k＝0，1，2，…，N-1)，根据时间弛豫法，取延迟时间T＝1/512秒，通过延时构成M维向量X(k)＝{V(k)，V(k＋T)，…，V(k＋(M－1)T)}，取M＝3，每个导联的脑电信号V(k)构成一个三维嵌入空间的向量集，并采用下列公式表达： $\left[\begin{array}{ccc}X\left(1\right),& X\left(2\right),& X\left(3\right)\\ X\left(2\right),& X\left(3\right),& X\left(4\right)\\ \…\…& \…\…& \…\…\\ X\left(n\right),& X(n+1),& X(n+2)\end{array}\right]$

将此三维相空间的每一坐标轴作六等分，得到信息熵H(S)： $H\left(S\right)=\underset{i=l}{\overset{N}{\mathrm{\Σ}}}{P}_s\left(S_i\right)\mathrm{lo}{g}_2\[P_s\left(S_i\right)\](N=6^3)$

以S和Q分别代表两个导联的脑电信号所构成三维相空间的格点序列；若把两个三维空间联合起来构成一个六维相空间，同样将每一坐标轴作六等分，再求信息熵H(S，Q)得： $H(S,Q)=-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{SQ}}\left(S_i{Q}_j\right)\mathrm{lo}{g}_2\[P_{\mathrm{SQ}}\left(S_i{Q}_j\right)\](N=6^6)$ 据此计算两个导联的相互信息传输量MI(S，Q)：

        MI(S，Q)＝H(S)＋H(Q)－H(S，Q)。

3.根据权利要求1所述的大脑头皮电位信息图形的生成方法，其特征在于所述的信息传输矩阵采用下列方法：

用时间延迟方法构成两个导联的相互信息传输时间序列SQ(τ)：

       SQ(τ)＝MI[S(t)Q(t＋τ)] τ＝T，2T，…，KT

SQ(τ)代表从S导联到Q导联的信息传输；

而相互信息传输时间序列QS(τ)为：

       QS(τ)＝MI[Q(t)S(t＋τ)] τ＝T，2T，…，KT

QS(τ)则代表从Q导联到S导联的信息传输，序列长为512点；矩阵中行代表输出信息，列代表接收信息。

4.根据权利要求1所述的大脑头皮电位信息图形的生成方法，其特征在于所述的复杂度和时间序列信号粗粒化，采用下列方法：

对所有(0，1)序列，定义b(n)

       b(n)＝(log₂n)/n

n为数据序列号。再定义C(n)

       lim C(n)＝lim(log₂n)/n

        n→∝    n→∝

定义复杂度Kc为：

       Kc＝C(n)/b(n)

对于符号(0，1)序列{S_i}^N _i＝1，定义复杂度C1(n)为：

       C1＝lim[log₂N_s(n)]/n

        n→∝

对于时间符号(0，1)序列{S_i}^N _i＝1，定义复杂度C2(n)为：

       C2＝lim[log₂N_f(n)]/n

        n→∝

C1反映了时间序列在相空间运动轨道的随机性结构复杂度；而C2则反映了时间序列在相空间运动轨道的限制性结构复杂度。

Images