0 引言

癫痫是一种严重危害人类安全健康的常见慢性脑部疾病，全世界大约0.5%的人患有癫痫，使运动感觉、自主神经、意识和精神状态等反复出现不同程度的障碍。癫痫发作时表现为大脑神经元兴性增高以及过度同步化放电，并导致短暂性中枢神经系统功能失常^[1-2]，患者承受了巨大的痛苦和危险，极大影响了患者及其家庭的生活质量。

脑电图(ElectroEncephaloGram,EEG)是在头皮上通过电极记录的脑细胞群自发性、节律性的电活动轨迹^[3]。当大脑出现病理性或功能性改变时，EEG就会发生相应变化^[4-5]。此外，脑电图具有时间分辨率高、价格低廉、获取容易且可进行实时监测等优点，日益成为癫痫临床诊断的重要辅助工具，在临床医学和脑科学的研究中都起着不可替代的重要作用^[5-7]。癫痫诊断目前大都由医生通过视觉检测患者的脑电图，根据经验诊断完成的。人工视觉检测存在费时、效率低、且缺乏标准的制约等不利因素，因此，癫痫脑电的自动识别与诊断研究具有重大的临床应用价值。

针对这一问题，众多学者致力于癫痫EEG的自动检测技术领域的研究，并提出很多自动检测方法。其中最常见的分析方法主要有时域分析、频域分析、复杂度分析等^[8-16]。这些自动方法诊断主要是对脑电信号进行分析，提取癫痫疾病的特征向量，以实现对健康人群与癫痫患者的区分以及癫痫不同时期的区分。大部分的参考文献基于对不同的EEG通道进行谱分析^[8-10]。Mporas等^[11]通过使用6阶自回归(AutoRegression,AR)模型估计癫痫EEG的功率谱密度并以此为特征进行分类分析，准确率高达近90%。此外，样本熵、多尺度熵、模糊熵、相空间分析等非线性混沌算法也被用来提取脑电特征，用于癫痫的自动诊断研究^[12-16]。

近年来，基于复杂网络的非线性时间序列的动力学分析受到了广泛关注^[17-20]。如Zhang等^[20]在相空间内将准周期时间序列的每一周期作为网络的一个节点，每两个准周期间的距离与所选取的阈值共同决定节点之间是否存在连边，以此将心电信号转化为复杂网络用于心电信号的分析。Gao等^[21]将一维时间序列重构到高维的相空间中，由此构成对应原时间序列的复杂网络，并将此相空间重构方法应用于识别与分析多相流的流形。Wang等^[22]通过相干性分析将多通道Alzheimer症脑电映射为复杂脑网络，通过网络特征探究Alzheimer脑功能的异常。

2008年Lacasa等^[23]提出了一种全新的一维时间序列映射为复杂网络的方法：可视图(Visibility Graph,VG)算法。其思想是将离散时间序列中的数据点作为网络节点，数据点之间满足可视化准则的连线定义为网络连接。可视图算法已成功应用于大脑神经障碍分析、心率分析、飞机湍射流分析、海洋潮汐记录分析、金融经济分析^[24-29]。随后，Luque等^[30]提出了水平可视图(Horizontal Visibility Graph,HVG)算法。可视图与水平可视图建网方法具有算法简单易实现、计算效率高等优点。但是这两种算法对不同类别的时间序列的识别能力尚存在问题，且对噪声的鲁棒性差。基于此，周婷婷等^[31]在此基础上提出了有限穿越可视图(Limited Penetrable Visibility Graph，LPVG)算法，并将其应用于油气水三相流电导波动序列的识别中。Pei等^[32]采用LPVG算法成功地区分开了针刺前、中、后的脑电信号。本文采用功率谱密度分析和有限穿越可视图(LPVG)算法对癫痫发作前、中、后阶段的脑电信号进行分析，旨在获取单通道的癫痫脑电序列的谱特征以及EEG序列映射成复杂网络的网络特征，用于癫痫脑电信号的自动识别与诊断，为癫痫的病理研究和临床诊断提供参考。

1 癫痫脑电的采集与预处理

研究所需癫痫病例由河北省唐山市工人医院神经内科提供。目前，临床实验获取癫痫病例已达到70余人。本研究在已获取的病例群中筛选特性相同或相近的5例中度异常病例进行算法研究。5例病例均为EP大发作，具体症状表现为右侧顶中后颞持续大量中高幅尖波尖慢。所有研究病例均为自愿参加，且已签署知情协议。

临床采集脑电实验中，电极帽电极位置按照国际导联10-20系统安放，其中，耳垂作为电极的参考地。实验设备的采样频率为256 Hz，硬件滤波器的为0.5~100 Hz。

实验过程中，对实验病例连续采集记录脑电信号10 min。从每位实验病例的每通道EEG信号中截取发作前(R1) 、中(R2) 、后(R3) 这三段长度为8 s的数据(数据长度为256×8=2048个采样点)进行分析。采用0~30 Hz的带通滤波器分别对每段数据进行预处理，并采用小波变换提取EEG信号4个子频带：delta(0~4 Hz)、theta(4~8 Hz)、alpha(8~15 Hz)以及beta(15~30 Hz)。图 1为theta频带下发作前、中、后三种状态下的脑电信号示意图，可以发现相比发作前、后，癫痫发作中的脑电信号出现异常。

2 分析方法

本文将采用功率谱密度分析和有限穿越可视图算法分别从能量和网络特征角度对癫痫发作前、中、后三个状态下的单通道脑电信号进行分析，具体过程如图 2所示。

2.1 功率谱密度算法

本文采用线性自回归(AR)模型估计法来估计EEG信号的功率谱密度。

AR模型如式(1) 所示：

$x\left( n \right)=-\sum\limits_{k=1}^{p}{{{a}_{k}}x\left( n-k \right)}+u\left( n \right)$

(1)

其中：x(n)为待分析时间序列，u(n)为外加白噪声序列，p是AR模型的阶数，本文中取p=10，a_k为模型参数，由Burg算法估计得到，其基本原理是使前向预测与后向预测的均方误差之和最小。

因此，估计出AR模型参数后，可以得到序列x(n)的功率谱密度，如式(2) 所示：

${{{\hat{P}}}_{\text{BURG}}}\left( f \right)\hat{e}_{P}^{2}/1+\sum\limits_{k=1}^{p}{{{a}_{k}}{{\text{e}}^{-\text{j}2\pi fk}}}$

(2)

其中$\hat{e}_{P}^{2}$是最小的预测均方误差之和^[33-34]。

2.2 有限穿越可视图算法

可视图(VG)算法就是通过可视图准则将一维时间序列映射为复杂网络。其思想是将离散时间序列中的数据点作为网络节点，数据点之间满足可视化准则的连线定义为网络连接^[23]。

可视化准则如下：长度为N的时间序列{X_i|1,2,…,N}，其中任意两点X_a与X_c之间存在连接等价于这两点间任意一点X_b，满足式(3) ：

${{X}_{b}}\le {{X}_{e}}+\left( {{X}_{a}}-{{X}_{e}} \right)\cdot \left( c-b \right)/\left( c-a \right)$

(3)

图 3(a)中用黑色垂直的直方条表示时间序列的数据点，同时作为网络的节点，直方条的高度表示对应的幅值，根据上述可视化准则将满足条件的点连接起来，用实线表示，构成可视图网络。VG网络具有如下几点性质：每一个点至少和其前后两个邻点相连接；网络连接具有无向性；横轴和纵轴坐标尺度变化或者经过仿射变换后，可视性保持不变。此外，VG算法能够继承原始时间序列的非线性动力学特征，即应用VG算法可将周期时间序列转化成规则网络，随机时间序列转化成随机网络，分形时间序列转化成标度网络。

水平可视图(HVG)算法如图 3(b)所示。网络的节点定义与网络的连边方式与VG相同，不同的是，网络节点间的连边必须是水平直线，且不能穿越其他的直方条。因此，对于同一列时间序列，水平可视图是可视图的一个子图，且节点的平均度值较小^[30]。

有限穿越可视图(LPVG)算法如图 3(c)所示。定义有限穿越视距为N，当且仅当两个节点所对应的直方条之间的连线穿越其他直方条次数满足n≤N时，两节点之间存在连边^[31]。图 3(c)中实线表示基于VG算法的连边，灰色虚线表示新增有限穿越连边。本文仿真中，有限穿越视距设为N=5。LPVG算法同样具有连接无向性与仿射变换或者坐标轴尺度变换后可视性不变的特点，且LPVG网络具有更强的连接性，长程连接能够反映序列波动趋势。因此，LPVG算法对噪声和干扰具有较强鲁棒性，在含噪情况下仍能较好地继承信号的动力学特性。

2.3 网络参数

复杂网络的拓扑结构可由以下网络参数^[35-37]来表征。

1) 节点度和度分布。

节点度衡量单个节点作为功能网络中多个连接交叉点的重要程度。网络中单个节点的度值等于该节点的连接总数。网络的平均度定义为所有节点度的平均值：

$K=\left\langle {{K}_{i}} \right\rangle =\frac{1}{N}\sum\limits_{i=1}^{N}{{{K}_{i}}}$

(4)

其中：N是网络的节点个数，K_i是节点i的度值，节点度K的分布情况可以用分布函数P(K)表示。

2) 平均路径长度。

平均路径长度是网络中所有点对的距离的平均值，度量了网络的整体连通性与不同节点间的集合程度。平均路径长度定义为某一节点到任一其他节点的平均距离：

$L=\frac{1}{N\left( N-1 \right)}\sum\limits_{i,j,i\ne j}{{{d}_{ij}}}$

(5)

其中：d_ij为节点i和j间的最短距离，N为网络中节点总个数。另一个重要的网络特征参数：聚类系数衡量一个节点最近邻节点之间的连接情况，与平均路径长度呈现相反趋势。

3) 图指数复杂度。

图指数复杂度(Graph Index Complexity,GIC)用于衡量功能网络的复杂度。定义式如下：

${{C}_{{{\lambda }_{\max }}}}\text{=}4m\left( 1-m \right)$

(6)

$\begin{align} & m=\left[ {{\lambda }_{\max }}-2\cos \left( \pi /\left( N+1 \right) \right) \right]/\left[ n-n- \right. \\ & \left. 2\cos \left( \pi /\left( N+1 \right) \right) \right] \\ \end{align}$

(7)

其中λ_max为功能网络邻接矩阵的最大特征值。在无向二值网络中，总有2 cos(π/(n+1) )≤λ_max≤n-1成立，因此，图指数复杂度在区间变化。C_λmax越大表明网络结构越复杂。

3 模型验证

在本章中，用正弦周期序列、随机序列、H=0.5(H为Hurst指数)的布朗分形序列以及Lorenz混沌序列来研究有限穿越可视图(LPVG)、可视图(VG)以及水平可视图(HVG)这三种时间序列建网的可视化算法的性能。此外，对四种序列分别添加10 dB、20 dB、30 dB的高斯白噪声信号。简明起见，四种典型时间序列以及添加10 dB高斯白噪声的时间序列如图 4所示，而其他噪声强度的序列在图 4中不予表示。

分别采用VG、HVG以及LPVG 3种可视图算法对四种典型序列进行建网，可以发现3种可视图算法都具有较好的还原序列的动力学特性。提取3种可视图网络的聚类系数以及平均路径长度，如表 1~2所示。可以发现，随着噪声强度的增大，与VG和HVG相比，LPVG算法提取的聚类系数和平均路径长度，其波动率均为最小，可见LPVG对含噪信号具有较好的适应性，具备更佳的抗干扰能力。

本文对3种可视图算法进行了时间复杂度和空间复杂度分析，计算结果如表 3所示。时间复杂度是度量算法执行的时间长短；而空间复杂度是度量算法所需存储空间的大小。

由表 3可知，HVG、VG以及LPVG这3种可视图算法的时间复杂度依次增大，表明算法的复杂程度逐渐增大，但计算所需的存储空间均为O(1) 量级。此外，3种可视图算法的计算时间如图 5(a)所示。仿真时间序列长度设置为2500个点。图 5(b)为采用VG与HVG算法对仿真序列的第2000~2500点建网时所耗时间的放大。

由图 5可知，与VG和HVG相比，LPVG算法时间复杂度远大于其他两种可视图算法，且仿真消耗时间也远远大于两种可视图算法。

虽然有限穿越可视图(LPVG)算法的时间复杂度较高，计算时间较长，且在信号识别能力上存在一定的局限性，但是抗噪能力强。在一定程度上，有限穿越可视图以牺牲相邻节点的相关性，来抵抗噪声点的干扰，达到防止重要节点被噪声隔断的目的，能够更好地反映时间序列的动力学特性和波动趋势。如果有限穿越视距N取值小，局部相关性并未全部抹杀。有限穿越可视图在牺牲微观相关性的情况下更能保留信号的宏观信息。

4 实验结果分析

本章将第2章提出的算法应用到癫痫脑电信号发作前、中、后状态的识别上。由于癫痫脑电信号在临床采集过程中含有大量噪声，因此，本章采用功率谱密度(Power Spectrum Density,PSD)以及抗噪能力强的有限穿越可视图(LPVG)算法对脑电信号(EEG)进行分析。

首先将每通道EEG提取为发作前(R1) 、发作中(R2) 、发作后(R3) 三段数据。随后，采用小波分解将三段EEG分解为delta、alpha、theta、beta四个频带，并在每个频带下采用AR Burg参数模型法估计分别19通道脑电信号的功率谱密度。最后，对5名患者的PSD计算结果求平均值，计算结果如图 6所示。

从图 6中可知，在delta频带下(图 6(a))，大脑普遍EEG能量较高，并且在癫痫发作前与发作后，脑区的能量分布有一定的层次性，而在发作期间，整个脑区能量异常增强，各脑区的能量值显著升高；而在theta频带下(图 6(b))，发作期间整个脑区的能量异常减少，脑区之间能量分布的区域性差异显著降低，反映出大脑脑区活动趋于一致性；在alpha频带下(图 6(c))，发作前、中、后三个阶段，脑区能量的变化不明显；在beta频带下(图 6(d))，脑区能量的分布值很低，且不存在明显变化。

为了进一步分析各个脑区能量变化的规律，本文将大脑分为前额区(F)、左颞区(LT)、中心区(CP)、枕区(O)、右颞区(RT)五个脑区。接下来在各个脑区上，对癫痫发作前、中、后(R1，R2，R3) 三个阶段的EEG的功率谱密度(PSD)值进行单因素方差分析，结果如图 7所示。当p＜0.01(图中用“*”注明)，表明状态之间存在显著性差异。图 7中，在delta频带下(图 7(a))，脑区能量的变化与地形图分析的结果一致，发作中各个脑区的能量值均高于发作前。与发作后相比，只有前额区、右颞区的能量值相近，并且发作期间各个脑区之间能量的变化没有发作前、中显著；发作中与发作前、后的能量只在右颞区没有显著性区别。delta频带下，在发作中各脑区能量的异常增强，而theta频带下(图 7(b))，除中心区以外，各脑区能量均低于发作前和发作后。从能量分布的区域性角度看，发作前、后不同脑区的能量存在较明显的区别，而发作中各脑区能量值趋于一致；除中心区以外，其余各脑区发作中与前、后能量存在显著差异。而alpha和beta频带下(图 7(c)、(d))的能量值较低，各脑区能量的变化无明显规律。

接下来，对发作前、中、后能量存在显著变化的theta频带下的第六通道EEG信号C4采用有限穿越可视图(LPVG)进行分析。序列长度选取为2000个数据点。图 8为构造的复杂网络的邻接矩阵，构造的复杂网络能较好地还原脑电序列的时域特性，并且对噪声具有极强的适应能力。由图 8所示可知，采用LPVG算法提取的发作前、中、后三种状态下的EEG序列的邻接矩阵中看出：theta频带下，矩阵中社团连接主要集中在对角线上，其他区域的社团连接面积都较小。发作前、后(图 8(a)、(c))，对角线上社团连接大小随机分布，而发作中(图 8(b))，社团连接面积有所增大，但分布呈现一定的规律性，表明发作中网络的复杂性有所降低。综上，发作发作前、后，网络模块分布无规则，各模块大小有明显区分；而发作中，网络中各模块大小趋于一致，并且小规模的网络各模块明显减少，各模块的独立性有所提高，网络的复杂度降低，相对趋向于规则网络。

图 9为theta频带下网络的节点度和度分布。可以发现，发作中(图 9(b))，网络的度值明显高于发作前、后(图 9(a)、(c))，且度分布的峰值增大且峰值右移，概率最大的度值接近40，而发作前、后概率最大的度值在20左右。而发作后与发作前相比，节点度分布要分散一些，但总体变化不明显。

图 10为采用LPVG对癫痫发作前、中、后的EEG序列在delta、theta、alpha和beta四个频带下分别构建，并提取网络特征参数：平均路径长度和图指数复杂度(GIC)。由图 10(a)可知，在theta和alpha频带下，前、中、后三种状态下区分不明显；在delta频带下，发作中网络的平均最短路径高于发作前、后，表明网络信息传递减慢；在beta频带下，发作中网络的平均路径长度变小，考虑到高频带包含的信息量较少，因此，对分析结果没有明显影响。图 10(b)可知，网络的复杂度随着频带增高而降低，在delta和theta频带下，发作中网络的复杂度要低于前、后，其中在delta频带的衰减最为显著，表明癫痫发作期间，虽然delta频带脑电能量急剧增强，但脑网络的复杂度降低，与EEG信号呈周期性变化相对应；而在alpha和beta频带下，复杂度改变不明显，对应高频带下网络信息量较低。

5 结语

为探究癫痫发作对脑电信号的影响，本文采用功率谱密度和有限穿越可视图算法算法分析癫痫发作前、中、后三种状态下的单通道脑电信号能量分布以及脑电信号的网络特征。研究发现，癫痫发作时，脑电信号能量分布发生明显改变，能量向低频带转移，表现为delta频带能量显著增加，而theta频带能量明显降低。对脑区能量分布的进一步研究发现，在delta频带下，发作前与发作后各脑区能量均有显著性差异；在theta频带下，除中心区外，各脑区均有显著性差异，表明癫痫发作对大部分脑区在低频带的脑电信号影响较大。对第六通道C4在theta频带下的发作前、中、后三种状态下对应的EEG信号进行有限穿越可视图(LPVG)分析，发现三种状态下，LPVG构建的复杂网络具有不同形态的拓扑结构，发作中时更接近于规则网络。与发作前、后相比，在高频带(alpha和beta)下，发作中网络的平均最短路径减小，而在低频带(delta和theta)显著增大，且网络的复杂度明显降低。本文研究结果表明，功率谱密度和有限穿越可视图算法能够有效表征癫痫前、中、后三种状态下脑电信号能量分布和单通道信号可视化后的网络拓扑结构的异常，为癫痫的病理研究和临床诊断提供帮助。