脑电信号的非线性动力学分析①

刘桂珍，于 影，李宪芝，闻邦椿

(1.佳木斯大学机械工程学院，黑龙江 佳木斯154007;2.东北大学机械工程与自动化学院，辽宁 沈阳110004)

0 引 言

脑电图是脑神经细胞电生理活动在大脑皮层的总体反映，它蕴含大量的生理和疾病信息，在临床医学中常作为某些脑疾病的诊断依据和提供辅助治疗的手段.在传统研究方法中，脑电信号被视为一种准平稳信号，通过谱分析、正交变换、逆滤波、神经网络等方法来分析脑电信号［1］，虽取得一定的成效，但仍有许多不足.自1985 年Babloyantz［2］首次运用非线性动力学理论研究EEG(electroencephalogram)信号，得出人类在睡眠时EEG 信号处于混沌状态伊始，脑电信号的研究迅速进入到非线性动力学时代.EEG 信号起源于一个高度的非线性系统，基于人体自身生理信息的复杂程度和生命现象的多重耦合［3］，利用非线性理论定量描述出脑电信号的复杂动力系统的演变信息，对于揭示生命体征更好的造福人类，具有重要意义.

本文基于非线性动力学理论，对LILEY 建立的EEG 的动力学模型进行了相关维数和Lyapunov指数进行了计算，得出脑电信号具有混沌的结论.

1 EEG 动力学方程

LILEY 等人以人脑生理结构和医学解剖为基础，用非线性微分方程定义了脑神经元集群处于兴奋和抑制两种状态时的数学模型，其微分方程如方程组1［4］所示.

式(1)［5～15］中的1 ～2 个方程分别表述的是脑神经元集群处于兴奋型e 和抑制型i 两种状态时的平均膜电位he和hi随时间变化的函数关系;第3 ～6 个方程表述的是突触活动与时间的变化关系，其中Ijk(j，k=e，i)表示前者j 对后者k 的神经元作用;S 描述的是对应e，i 两类神经元细胞膜电位在等价平均点时放电率的函数关系;A 和B，a 和b 分别表示对应e，i 两类神经元集群的突触后膜电位峰值的振幅和突触速率常数;τe，τi是两类神经元集群的细胞膜时间常数;her，hir和heeq，hieq分别表示e，i 两类神经元集群的静止和平衡电位;pee，pie和pie，pii代表e，i 两类神经元输入.Nee，Nie和Nei，Nii分别表示e，i 两类神经元接受其他两类神经元集群传递信号的总数.

式(2)中的q=e，i，是对式(1)中的函数S 量化关系的补充，hq表示两种神经元集群的平均细胞膜电位he和hi;qmax=emax，imax，代表两种神经元集群的最大阈值emax和imax;θq和sq分别表示为对应于e，i 最大阈值时的电位值和标准偏差.

其中，

图1 兴奋型神经元集群平均膜电位he

2 基本原理

2.1 相关维数

相关维数为非整数维数是数值识别混沌运动的一种判据.

本文采用Grassberger 和Procaccia 在1983 年提出的针对时间序列计算吸引子维数的方法［10］.其目的是求出距离接近重构相空间中的点占其重构空间中的比重.其计算方法如下［11］:

(1)对一个离散时间序列{x(t)}，先以嵌入维数m 和时间延迟τ 重构，重构后空间中的一个点可表示为

(2)计算相关积分C(r).设重构空间中有N个点，r 为两两配对点对之间任意两个向量之间的距离，C(r)是满足|Xi－Xj|＜r 的点对数目的概率，其计算公式如下:

式(4)中，M 表示在m 维空间中向量点的个数，H 是阶跃函数，|Xi－Xj|表示相空间中任意两个向量点之间的距离.

(3)做logr－logC(r)二维曲线.

(4)计算相关维数D2.将(3)所做曲中近似线性部分拟合成直线，求该直线的斜率，即为相关维数D2的估计值.用公式表示即为:

图2 兴奋型神经元集群细胞膜电位的相关维数

图3 兴奋型输入pee 作参数的EEG 分岔图

2.2 Lyapunov 指数

混沌运动的另一种判据是Lyapunov 指数为正数.

求解最大的Lyapunov 指数的方法如下:

(1)在m 维相空间重构中找到n 个向量点，得到的向量集记为

(2)计算该向量两点间最近的距离记为d1:

(3)取时间步长Δt，建立一个新向量，找出与新向量最近的向量距离d2，定义该时间段内系统的指数增长率λ1，则

(5)增加嵌入维数m，重复步骤(2)～(4)，当L1随嵌入维数m 的改变趋于稳定，即为所求Lyapunov 指数.

图4 兴奋型输入Pee 作参数时LLE 变化图

3 数值仿真

对模型(1)作降维处理，得到10 个一阶常微分方程.模型(1)中参数分别取a=0.49ms－1，b=0.592ms－1，A=0.81mV，B =4.85mV，emax=imax=0.5ms－1，heeq=45mV，her=hei=－70mV，hieq=－90mV，Nee=Nei=3034，Nie=Nii=536，pie=pii=0ms－1，se=si=5mV，pei=11.9ms－1，θe=θi=－50mV，τe=9ms，τi=39ms［7］.

当pee=12.9ms1时，选取迭代步长0.01，时间长度1000，迭代点数为100000，为减少初值对系统的影响，舍弃前10000 个点.he是大脑皮层某局部区域兴奋性神经元集群的平均膜电位，它与EEG信号成正比［9］.图1 是he的信号波形.

4 结 论

(1)模型(1)的动力学行为主要是由控制参数pee和pei决定的［6，7］.

(2)图2 表明，对he信号做关联积分，当嵌入维数m=8 时，求出关联维数D2=1.7，关联维数的非整数，证明脑电信号呈现混沌状态.

(3)由数值仿真图3 可见，当改变兴奋性神经元集群兴奋性输入参数pee值时，EEG 经历了稳定周期→两倍周期→倍周期分岔→混沌→半周期分岔→稳定的变化历程，结合图3 和图4，pee ∈［13.1ms－1，13.6ms－1］附近，出现周期窗.

(4)图4 表明，最大Lyapunov 指数LLE 值出现在pee∈(12.64，14.36)内，且为正数，证明脑电信号呈现混沌状态.

以上分析说明，EEG 的非线性特征与输入不同的电信号有关.在一定参数区间内，大脑呈现混沌状态.

［1］ 孟欣，欧阳楷.脑电信号的几个非线性动力学分析方法［J］.北京生物医学工程，1997，16(3):135－139.

［2］ Bab loyantz A，Salazar JM，N icolisc Eviclence of Chaotic Dynma Ics of Brain Activity During The Sleep Cycle［J］.Physics letters A，1985，111(3):152－156.

［3］ 缪经良.振动与波利用技术的新进展［M］.沈阳:东北大学出版社，2000.

［4］ LILEY D T J，CADUSCH P J，DAFILIS M P.A Spatially Continuous Field Theory of Electrocortical Activity［J］.Computation in Neural Systems，2002，13(1):67－113.

［5］ 刘桂珍;刘丽群;于影;闻邦椿.非线性振动理论在脑电信号分析中的应用［J］.中国工程机械学报，2012，2.

［6］ 李春胜王宏，赵海滨等.EEG 模型的混沌动力学特性研究［J］.系统仿真学报，2010，22(10):2262－2264.

［7］ 王兴元;谭贵霖.Liley 模型的模拟EEG 信号的非线性预测和分析［J］.计算物理，2007，(09).

［8］ 王兴元，骆超，谭贵霖.EEG 动力学模型中混沌现象的研［J］.生物物理学报，2005，(09).

［9］ 向学勤;庞全;范影乐;薛凌云.EEG 动力学模型中随机共振现象的仿真研究［J］.系统仿真学报，2008，(08).

［10］ KORN H，FAURE P.Is There Chaos in the Brain?II.Experimental Evidence and Related Models［J］.Comptes Rendus Biologies(S1631－0691)，2003，326(9):787－840.

［11］ GRASSBERGERP.Characterization of strange.attrators［J］.Physial Review Letters，1983，50(5):346－349.

［12］ 黄小娜;张军英.;癫痫脑电信号的非线性分析方法［J］.河西学院学报，2010，(04).

［13］ 徐琳，许百华.非线性动力学脑电信号分析方法的研究与应用［J］.心理科学，2005，(05).

［14］ 李春胜，王宏，赵海滨，等.脑电混沌维数复杂度连续检测方法的研究［J］.中国生物医学工程学报，2009，(12).

［15］ 袁野.癫痫脑电的非线性方法分析［D］.长春:吉林大学，2009.