（计算机软件与理论专业论文）脑电时间序列非线性高性能计算与仿真.pdf

摘要 摘要 人脑是复杂的非线性动力学系统，脑科学研究已成为2 1 世纪最重要的研究 热点之一。自上世纪2 0 年代脑电( e e g ) 被发现以来，人类便开始利用脑电对 大脑进行无创伤性研究，从而脑电在许多领域都起到了重要作用。传统的方法 从频谱和统计学的角度研究脑电时间序列，能够对脑电的部分特征给出详细的 解释，为脑电在这些领域中的有效应用做出了一定贡献。然而，这些方法无法 对大脑的非线性动力学性质进行深入有效的分析。本文讨论脑电时间序列的非 线性动力学高性能计算与分析方法。 非线性时间序列动力学分析分为相图重构和几何特征值刻画两个重要步 骤。两个步骤现存较多分析方法，本文通过分析比较这些方法后选择了较佳的 且适合脑电分析的方法。然而这些方法的均存在时空复杂度较高的缺点，给对 实际系统分析带来困难，本文对其进行了优化、改进和并行化，较大程度提升 了其效率，并提高了计算结果正确率。本文还设计了e e g 非线性动力学高性能 计算与分析系统，可以高速分析实际脑电数据。本文最后对睡眠脑电数据和癫 痫脑电数据进行了计算，并对睡眠阶段进行了仿真预测，通过评估说明了结果 的正确性。 关键词：脑电非线性分析优化并行化仿真预测 i i a b s t r a c t _ _ - 一一一 ab s t r a c t t 1 1 er e s e a r c ho nt h eh 啪a n b r a i n ，w h i c hi so n eo ft h em o s tc o m p l e xn o n i i n e a r d y l l a m l c ss y s t e m s ，h a sb e c o m eo n eo ft h eh o t t e s ta 诧弱a n dw i l lb e 向l ld e v e l o p e di n t 1 1 e21s t c e n t u r y - a r e rt h ed i s c o v e 巧o fe l e c t r o e n c e p h a l o g m m ( e e g ) i n19 2 9b y h a i l sb e 唱e r ，e e gh a sb e e nu t i l i z e da san o n i n v a s i v em e m o d t os t u d yt h eb r a i na n d p l a y s 觚 i m p o n a n tr o l ei n m a n y a r e a s t r a d i t i o n a l m e t h o d s ， s u c h嬲 舭q u e n c y - s p e c a 1 1 ds t a t i s t i c s ，c a n tc o r r e c t l ye x p l a i nt h en o n l i n e a rd y n a n l i c a l c h a r a c t e r i s t i c so fh 啪a 1 1b r a i n t h i st h e s i sd i s c u s s e sh i 曲巾e r f o 衄a n c en o n l i n e a r d y n a m i cc o m p u t a t i o na n da i l a l y s i sm e t h o do fe e g r e c o n s t i t m i o no fp h a s eg r 印h i c sa 1 1 dc o m p u t a t i o no fg e o m e t r i c a le i g e n v a l u e s a r el m p o n a n tp r o c e s s e si nn o n l i n e a rd y n 锄i c a la n a l y s i s t h e r ea r em a n yd i f f e r e n t m e t h o d smt h et v v r op r o c e s s e s ，a n dw ec h o s et h eb e s tm e t h o d sw h i c h a r ea l s of i tf o r e e g sa n a i y s i sb yc o m p a r i n gt h em e t h o d s i t sd i m c u l tt oa i l a l y z en o i l l i n e a rs y s t e m s u c ha se e gr e a l t i m eb e c a u s et h ea l g o r i t h m so fp h a s e 蓼a p h i c sr e c o n s t i t u t i o na i l d g e o m e t i i i c a le i g e n v a i u ec o m p u t a t i o na r ec o m p l e xo nb o t ht i m ea 1 1 d s p a c e w e o p t i m i z e da n dm o d i f j e dt h ea 】g o r j t h m st or e d u c et h e i rc o m p l e x i t ya n dt 0i m p r o v e t h e i ra c c u r a c y ，a n dm e nw ep a r a l l e l i z e dt h ea l g o r i t h m s w ed e s i g n e d 肌e e g ，s h i 曲p e 怕m l a n c en o n l i n e a rd y l l 锄i ca n a l y s i ss y s t e mw h i c hc a na n a l v z ee e g r e a l - t i m e w e c o m p u t e ds i e e p i n ge e ga n d e p i l e p t i c e e gu s i n ge e g ，s h i g h - p e 墒咖a n c en o l l l i n e a rd y n 锄j ca n a j y s i ss y s t e m ，a n dt h e nf o r e c a s t e dn i 江m s l e 印s t a g e ，w eg o tg o o dr e s u l ta tl a s t k e y w o r d s ：e i e c t r o e n c e p h a l o g r 锄 ( e e g ) ， n o n l i n e a r a 1 1 a l y s i s ，o p t i m i z a t i o 凡 p a r a l l e l i z a t i o n ，s i m u l a t i o na n df o r e c a s t 中国科学技术大学学位论文原创性和授权使用声明 中国科学技术大学学位论文原创性和授权使用声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工 作所取得的成果。除已特别加以标注和致谢的地方外，论文中不包 含任何他人已经发表或撰写过的研究成果。与我一同工作的同志对 本研究所做的贡献均已在论文中作了明确的说明。 本人授权中国科学技术大学拥有学位论文的部分使用权，即： 学校有权按有关规定向国家有关部门或机构送交论文的复印件和电 子版，允许论文被查阅和借阅，可以将学位论文编入有关数据库进 行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论 文。 保密的学位论文在解密后也遵守此规定。 作者签名：遣盐扭一 2 口够年多月石日 蟪矗6 o 刀巧y ，移仞 第一章绪论 1 1 引言 第一章绪论 大脑皮层作为一个整体，其神经元活动所产生的电位变化，可以通过大脑 这个容积导体，反映到大脑表面。在头皮上安置引导电极，通过脑电图仪所记 录到的皮层脑电活动的图形，称为脑电图( e l e c t r o e n c e p h a l o 掣锄，e e g ) 。 脑电活动研究的发展始于1 9 世纪早期。1 8 7 5 年，英国生理学家硒c h 砌 c a t o n ( 1 8 7 5 ) 首先使用验流计从动物的脑子中记录到了微弱的电流。德国精神病 学家h 觚sb e r g e r 被称为“脑波之父”，于1 9 2 9 年首次成功从人的大脑皮质处检 测到活动电流，并首创脑电图这一术语( h b e 玛e r ，1 9 2 9 ) 。此后，b e 唱e r 先后发 表了重要的科技论文达2 0 篇，总的名称都叫“0 ne 1 e c t r o e n c e p h a i o g r 锄o f m a n 口g l o o r ，1 9 6 9 ) ，从生理、心理及i 临床各个方面对脑电活动进行了初步的探讨， 奠定了e e g 研究的基础。1 9 3 4 年，英国的a d r i a i l 和m a t t h e w s ( 1 9 3 4 ) 发明了增 幅器及其同步的记录装置，首次清晰地在纸上描绘出了a 波的形状。1 9 3 5 年， 美国科学家g i b b s ，d a v i s 和l e n n o x 记录下了3 h z 的棘波，g i b b s 还发现了在昏 睡状态下脑电波呈现出显著的徐波化现象。1 9 3 9 年周期分析装置开发成功。在 对脑电活动不断探索的过程中，各种检测和分析技术也随之诞生了，其中诱发 电位检测分析技术的发展将脑电的研究引入了新的领域，诸如闪光刺激诱发、 睡眠诱发等试验方法，进一步提高了人们对大脑的认识和了解。 如今，e e g 已成为诊断大脑疾病，检测大脑功能和分析大脑活动的一个重 要工具，但是由于e e g 信号强烈的非平稳、非线性特性，分析起来十分困难， 实现e e g 研究的突破必须寻找新的方法，因而有关e e g 分析方法的研究正在 不断发展和日益深化。 第一章绪论 1 2e e g 的特点及分析方法 1 2 1e e g 的特点 脑电波是两个电极问的细胞群的自发性、节律性电位差的记录。但是在头 皮上安放电极所采集到的电位变化并不是单一神经细胞的电位变化，而是由许 多神经细胞综合起来的电位活动。记录到的脑电图可以看作是由各种类型的波 叠加而形成的图形或曲线，在时域平面上表示的话，其纵轴为电位活动的振幅 高度( 电压) ，横轴为时间( 电位活动的时间长短) 。通常，在描述波形时一般 使用周期( 频率) 、波幅和相位( 极性) 等基本要素。 1 ) 周期与频率：通常，正弦波的周期是指一个波与下一个波之间的距离， 用时间表示。脑电波的周期稍有不同，一般指一个波底到下一个波底 的距离或者一个波顶到下一个波顶的距离对基线的投影，用时间( m s ) 表示。在临床上一般使用频率来代替周期表示脑电波的特征。频率是 指同一周期的脑波在秒钟内以基线为准出现的次数。 2 ) 波幅：又称振幅或电压，反映了电位差的大小，也就是电压的高低。 正弦波的振幅是指由基线到波顶或者波底的距离。但是在脑电图中通 常用平均振幅来表示，从波顶画一条垂直于基线的直线，与前后两个 波底之间的连接线交叉后，交叉点至波项的距离为波幅。 位相：又称极性、相位，指同一部位在同一导联中所导出的脑波，在 前后不同时间旱波的位置；或者两个不同部位在同一时间( 某一瞬间) 里所导出的脑波的位置关系，即时问关系。脑电波的位相有正相和负 相之分，一般以基线为标准，朝上的波为负相波，朝下的波为j 下相波。 4 ) 脑电波的分类：根据脑电波的频率，通常可以将脑电划分成四个波段： 6 ( 1 3 h z ) 、0 ( 4 7 h z ) 、q ( 8 1 3 h z ) 、b ( 1 4 3 0 h z ) ( 谭郁玲，1 9 9 9 ) 。 通过分析脑电来研究大脑和脑疾病诊断具有如下几点优势： 1 ) 安全：对于表皮脑电波的提取对大脑无任何损伤，只是根据大脑皮层 产生的电活动在表皮出现的电效应进行记录； 2 ) 经济：脑电波的采集无需耗费大量的人力财力，只需具备安静封闭的 环境和并不昂贵的脑电仪就可以采集到非常精准的脑电波； 3 ) 方便：脑电波随时都可以采集，现有的许多便捷式脑电仪可以将电极 2 第一章绪论 和采集设备采用数字无线设备，非常小巧的脑电采集器可以让病人随 身携带，用于病人的长期监护。 4 ) 可靠：根据不同的需求可以选用不同精度的脑电仪，对于常规的临床 监护，可以选择精度相对较低的脑电仪；对于要求较高的科学研究， 如最近新兴的脑一计算机接口( b c i ) ，需要非常精确的脑电，甚至有 的采用植入式电极记录大脑皮层的脑电信号，称之为皮层脑电( e c o g ) ( t n l a l ，t h i n t e r b e r g e r ，gw r i d m a n ，e ta l ，2 0 0 4 ) 。 1 2 2e e g 的常规分析方法简介 对采集到的脑电进行分析和处理，最基本的方法就是由有经验的脑电专家 通过直接观察，识别出脑电图中的不同波形。这种方法需要有丰富的领域知识 和临床经验，而且需要耗费大量的时间和精力。这种方法是脑电发现之初最常 用的方法。随着脑电的应用不断发展，脑电图数量的不断增加，这种人工识别 的方法使得脑电专家不堪重负，因而，这种方法很难对病人进行长时间的全程 监护。 在计算机技术和医学影像技术的不断发展下，许多自动分析和可视化方法 为分析识别脑电信号提供了良好的辅助工具。根据目前的研究成果，可以将脑 电的研究方法分为三大类：时域频域分析、概率统计分析、非线性分析。这些 方法从不同角度研究脑电信号所反映的大脑状态特征。 1 ) 时频分析：这类方法从时域和频域的角度研究脑电信号中的变化，它们 是最初的直接观察法的发展和补充。时域分析方法使用统计学知识直接提取波 形特征，为研究者或医师提供进一步的分析或诊断信息，提取的特征有过零截 点分析、直方图分析、方差分析、相关度分析、峰值检测及波形参数分析、相 干平均等。， 脑电的频域分析方法主要采用功率谱估计和相干分析。功率谱估计方法将 按时间变化的电压信号变换为脑电功率随频率变化的谱图，以便直观地观察到 脑电节律的分布和变化情况。相干分析方法运用相干函数表示两个脑电波的振 幅、频率和相角的相对程度。 近年来，使用时频分析研究e e g 得到了较大的发展。尤其是以c o h e n 为代 表的二次时频分布和各种小波变换方法在e e g 分析中得到了广泛的应用，这些 第一章绪论 应用包括癫痫波的特征分析与监护，脑电的节律分析及诱发电位的检测与估计 盘蟹 1 了= o 2 ) 统计分析：在脑电分析的统计方法中，常用有独立成分分析( i c a ) 、主 成份分析( p c a ) 和概率随机分析。在多通道脑电分析中，i c a 方法假设脑电 的信号源时序上独立或几乎独立于大脑的整个过程，只是不同区域的通道反映 了不同脑功能的信号。 当没有大量数据进行充分分解时，可以用p c a 方法对信号进行线性分解。 p c a 的目标是在相互关联的通道数据空间寻找时序上相互正交的方向，当与主 要活动相关的头皮结构图不正交时，p c a 可将大脑活动状态的一部分组合到一 个或多个主成分中，而不是同i c a 方法那样将他们分离成不同的源。 脑电信号的随机分析假设脑电信号的复杂无序地不可解释性是因为脑电信 号的随机性，是大脑内部活动的无规则性引起的。随机方法通常用来分析系统 的平稳性，通过计算一些时间上的统计均值来讨论不同状态经历的随机过程， 随机过程中的下一个状态都是依赖于前一个状态，并且以一定的概率进行转移。 3 ) 非线性分析：非线性分析方法与前两种方法的出发点完全不同，前两种 方法均假设脑电信号是随机的，脑电信号的难以解释性是由于脑电是由未知系 统无规律运动所产生的结果，即未知系统是不确定的。而非线性分析方法假设 脑电是某个未知的非线性系统的一个变量或多个变量组合的输出信号，脑电信 号所表现出来的随机性是由于附加在上面的其他信号所引起的，脑电信号背后 的系统本身是确定的( h k a n t z ，t s c l u i b e r 1 9 9 7 ) 。对于非线性确定的系统所产 生的信号，用概率统计的方法进行分析显然是不能达到最终目的，只会产生一 些错误的解释。然而，对于一个未知系统的时间序列，在用非线性分析方法分 析之前，必须知道该时间序列是否是确定的，所以，非线性分析方法一般分成 两个阶段：首先检验系统的非线性确定性，然后进行非线性分析。 非线性检验的常用方法是使用替代数据法进行假设检验，替代数据一般是 对分析的时间序列数据通过线性随机变换而获得，替代数据要保证与原信号具 有统计一致性。然后提出一个线性随机过程的无效假设，通过对原时间序列和 所有替代数据集进行数理统计分析，如果在足够高的置信度下拒绝无效假设， 则认为目标系统是非线性确定的，否则就没有足够的证据说明系统的非线性确 定性。运用这种方法，许多文献已经判断出大脑在许多状态下的脑电时间序列 都具有明显的非线性特征( c b i 姗i e ，2 0 0 1 ) ( n k a n n a t h a l ，u r a i e n d r aa c h a r y 钆 4 第一章绪论 c m l i m ，e tm ，2 0 0 5 ) 。 非线性时问序列方法基于动力学系统理论，也就是说系统时间演化定义在 某个相空间中，这些系统会具有确定性混沌特征。因此要研究非线性时间序列 就需要具备刻画系统状态演化的相空间，w k 缸l e y ( 1 9 3 6 ) 的重构理论和融e n s ( 1 9 81 ) 的延迟重构定理提供了从实验时间序列重构系统相空间的理论基础。基 于重构的系统相空间，根据非线性系统的确定性，可以对时间序列进行噪音消 减、非线性预测、探测稳定性、计算复杂度等操作。对于具有非线性特征的系 统，非线性方法能非常有效地解释目标系统。本文主要讨论脑电的非线性高性 能计算与分析方法。 1 2 3 非线性时间序列分析方法简介 首先介绍时间序列的概念。从纵向看，时间序列是指存在于自然科学或社 会科学中的某一变量或指标的数值或观测值，按照其出现时间序列的先后次序， 以相同或者不同的时间间隔排列的一组数值。它是某一现象或若干现象在不同 时刻上的状态所形成的数据，反映的是现象以及现象之间关系的发展变化规律。 从横向看，时间序列也可以是若干相关现象在某一时间点上所处的状态按一定 顺序排列的一组数据，反映的是一定时间、地点条件下各相关现象之间存在的 内在数值联系。因此，从系统的意义上看，时间序列就是某一系统在不同时间 ( 地点、条件等) 的响应。时间序列按所研究系统复杂程度不同，可以分为线 性时间序列和非线性时间序列。显然脑电是时间序列，且人脑是一个非常复杂 的系统，所以脑电是非线性时间序列。 近年来，非线性科学得到了快速发展。实际问题中，面临大量复杂而无法 直接建立解析数学模型的非线性系统，对这些非线性系统可以通过实验或观测 的手段获得非线性时间序列，这些非线性时问序列中蕴含着丰富的系统动力学 信息，如何提取这些信息并应用到实际问题中去解释、分析，甚至控制原来的 复杂系统是非线性时间序列应用的一个重要方面。从2 0 世纪8 0 年代以来，非 线性时间序列分析的相关理论已取得了许多重大成果，在实际复杂系统中也得 到了有效的应用，为解决实际问题中像经济系统、生理系统、电力系统等这些 用传统的确定性和随机性方法都无法很好刻画的复杂系统提供了一种全新的研 究工具。 第一章绪论 1 3 研究内容和本文组织结构 1 3 1 研究内容 1 ) 选择适合分析脑电时间序列非线性时间序列动力学分析方法。 非线性时间序列的动力学分析方法由相空间重构和几何特征值刻画两个步 骤组成，目前两个步骤都存在多种分析方法，本文需要对这些方法进行系统的 分析比较，选择较优的且适合脑电时间序列分析的方法。 2 ) 对非线性时间序列分析方法进行优化改进及并行化。 非线性时间序列动力学分析方法时空复杂度均较高，计算耗时，难以达到 对一些实际非线性系统的快速甚至实时的分析。本文需要对它们进行优化改进 及并行化，以提高其计算效率，同时通过对部分算法进行改进以提高其准确性。 3 ) 设计一套e e g 非线性动力学高性能分析系统。 基于优化改进及并行化后的非线性时间序列分析算法，设计了一套e e g 非 线性动力学高性能分析系统。 4 ) 讨论时间序列线性变换对非线性分析结果的影响。 对这个问题的讨论源于数据的获取。 5 ) 对实际的e e g 数据进行计算和分析。 用非线性动力学高性能分析系统对实际e e g 数据进行了计算，并对结果进 行分析，然后对睡眠阶段进行仿真预测。 1 3 2 本文组织结构 第一章主要介绍脑电分析在脑科学研究中的意义以及脑电波的特点，回顾 现有的脑电分析方法。 第二章分析比较了非线性时间序列分析的各种方法，选择了适合脑电时间 序列分析的较佳方法。 第三章对相空间重构中确定嵌入延迟的c c 算法和确定嵌入维数的 l i a i l g y u ec a o 算法进行了优化改进和并行化，大大提升了它们的效率，并提高 了c c 算法的j 下确性。 第四章对关联维计算算法g p 算法进行了优化改进和并行，很大程度上提 6 第一章绪论 升了算法效率；并设计了e e g 非线性动力学高性能分析系统。 第五章用e e g 非线性动力学高性能分析系统对睡眠脑电、睁闭眼脑电和癫 痫脑电进行了分析，并进行了睡眠阶段预测。 第六章对硕士研究生期间的工作进行总结，列出了本文在脑电非线性时间 序列高性能计算和分析领域的工作，同时对脑电非线性时间序列高性能计算和 分析的研究在未来可开展的工作进行展望。 7 第二章非线性时间序列分析方法的比较 2 1 引言 第二章非线性时间序列分析方法的比较 2 0 世纪下半叶，非线性科学获得了前所未有的蓬勃发展。非线性科学是一 门研究非线性现象共性的基础科学，被誉为2 0 世纪自然科学中的“三大革命之 一”，与量子力学和相对论具有相同的地位( h a 0bl ，1 9 8 4 ) ( h a obl ，1 9 9 0 ) 。这 门学科几乎涉及自然科学和社会科学的各个领域，并不断地改变着人们对现实 世界的许多传统看法。 非线性动力系统理论是研究复杂系统的一种较好的方法，对于那些其演化 过程不可预测但其不可预测性可以预测的复杂系统，即确定性系统内具有内在 随机性的复杂系统，非线性理论是一种有效的研究工具。如果通过观测或实验 手段获得的时间序列是来自于线性系统，则时间序列可以表示成不同幅度的正 弦波的叠加，因此可以用传统的线性方法和f o 嘶e r 变换方法进行处理。但当获 得的时间序列来自于非线性系统时，无论多么高阶次的线性模型都不可能很好 地刻画该动力系统，对这样的时间序列，非线性时间序列分析方法称为重要的 研究工具。对于e e g 的分析，非线性方法与传统的分析方法相比有其自身的优 越性，它往往能够观测到传统方法所不能观测到的结构( c j s t 锄，t c a m v a n w 6 e r k o m ，w s p r i t c h a r d ，19 9 6 ) ，因此非线性分析方法在分析e e g 等非线性时间 序列方面显示出较好的应用前景。 当今非线性时间序列分析方法较多，本章工作首先需要分析比较现有的各 种方法，然后选择合适的方法进行分析及优化。 2 2 非线性时间序列分析方法 非线性时间序列分析可以看作是如图2 1 所示的输入输出过程，输入的是 没有任何背景知识的通过观测或实验手段获得的实际复杂系统的时间序列，输 出的是所研究系统的某些特征。 第二章非线性时间序列分析方法的比较 器裂眷嵩嚣蓊邕得的时间序列 图2 1 非线性时间序列分析的框架模型 下面简要介绍非线性时间序列分析方法的主要步骤，它包括两个主要步骤。 第一个步骤是重构相空间，它基于毗e n s 的嵌入延迟理论，即从观察到的时间 序列重构系统的状态空间轨迹；第二个步骤是刻画系统的几何不变量。下面首 先介绍相空间重构。 2 2 1 相空间重构 相空间重构是分析非线性动力学系统的第一步，它是以瞅e n s ( 1 9 8 1 ) 理论 为基础的。呲e n s 指出，由于非线性系统产生的轨迹经过一定时期的变化后， 会最终做一种有规律的运动，产生一种规则的、有形的轨迹( 混沌吸引子) 。系 统的任一分量的演化都是由与之相互作用着的其他分量所决定的，因此，这些 相关的分量的信息就隐含在任一分量发展的过程之中，这样就可以从某一分量 的时间序列中提取和恢复出动力系统原来的规律。为了重构一个等价的状态空、 间，只需考察一个分量，并将它在某些固定的时间延迟点上的观测量作为新维 处理，即延迟值被看作是新的坐标，它们确定了某个多维状态空间的一点。重 复这个过程并测量相对于不同时间的各延迟量，就可以产生出许多这样的点， 已经证明它可以将吸引子的许多性质保存下来，即用系统的一个观测量可以重 构出原动力系统的模型，并可以初步确定该系统的真实相空间的维数。由此可 见，相空间重构就是从时间序列出发创建一个多维状态的空间，它保持了原系 统的许多几何不变量不变。 设动力系统是由非线性差分方程 z 川= ，1 ( 乙)( 2 1 ) 表示的离散系统，或是由微分方程： 掣：f ( z ( r ) ) ( 2 2 ) = ， i7 j - il ，l 口z 表示的连续系统。其中z t i 或z ( t ) 是系统在时刻n 或t 的状态向量，f ( ) 是向 量值函数。时间序列 x n ) 是观测到的系统某一维输出，即 9 第二章非线性时间序列分析方法的比较 矗= 力( 乙) + ( 2 3 ) “ 戥 x 。= x ( 气+ 玎f ) = j f l 【z ( + 玎f ) 】+ 缈。 ( 2 4 ) 式中，h ( ) 为多元数量值函数；n 为在观测或者观测过程中由于技术手段不 完善或者精度不够引起的测量噪声。 根据f 融e n s 定理，当。= o 时，观察到的时间序列 x n ) 以向量 贾。= ( x 。，x 。一，x 。一( 肿1 ) ，) ( 2 。5 ) 形成m 维空间，只要l n 2 2 d + 1 ，动力系统的几何结构可以完全打开，其中 d 是系统中吸引子的维数，是t 正整数，称为延迟时间间隔。条件m 2 d + 1 是动 力系重构的充分但不必要条件，获得动力系统重构的整数m 叫做嵌入维数。状 态空间中的矗一毛+ 。的演化反映了未知动力系统毛一乏+ 。或z ( f ) 专z o + 1 ) 的演 化，并且状态空间中吸引子的几何特征与原动力系统几何特征等价，这意味着 原动力系统中任何微分或拓扑不变量可以在重构的状态空间中计算。 为了能在重构的相空间中刻画原动力系统的性质，需要j 下确地确定延迟时 间问隔t 和嵌入维数m 。因此如何选择恰当的嵌入时间间隔和嵌入维数就是相 空问重构的主要研究内容。 2 2 2 刻画系统的几何不变量 混沌现象最主要的特征是相空间中存在着奇异吸引子，而描述吸引子的基 本数学量是它的维数。另外在相空间重构过程中，重构参数嵌入时间t 和 嵌入维数m 的选择对重构质量以及后续的分析都有直接影响。因此，吸引子维 数计算在混沌分析中起着重要的作用。 奇异吸引子的维数一方面反映了吸引子结构的复杂性，另一方面也反映了 吸引子的信息量，因此维数的分析和计算不仅对混沌信号的相空间重构有重要 影响，而且也可以用它来刻画吸引子几何结构的特征。对于奇异吸引子，常见 的维数主要有：h a u s d o r f f 维数、信息维数、l y a p u n o v 维数、关联维数( c o 玎e l a t i o n d i m e n s i o n ) ( p g r a s s b e 唱e r ，i p r o c a c c i a ，19 8 3 ，a ) ( p - g r a s s b e 唱e r ，i p r o c a c c i a ，19 8 3 ， b 1 等，其中关联维数是最常用的一种维数度量，尤其对通过观测或实验获得的 时间序列。 在确定性系统中，关联维数就是生成相应复杂系统所必需的独立变量的个 1 0 第二章非线性时间序列分析方法的比较 数，规则的确定性系统有整数关联维数，而混沌系统有非整数关联维数，大于 此关联维数的下一个整数就是系统的独立变量的个数。数据量n 、嵌入延迟t 和嵌入维数m 等参量影响着关联维数的精确估计，因此这些参量的选取显得特 别重要。 2 3 嵌入延迟计算方法选择 t a k e l l s 理论认为，重构的结果不依赖于嵌入延迟，即可以任意选取嵌入延 迟。但是瞅e n s 理论的条件在实际应用中往往得不到满足。同时实践证明，嵌 入延迟1 的精确选取对于重构结果甚至吸引子的刻画起着重要的影响，因此精 确选择t 是必要的。对于嵌入延迟的计算，经常凭借的是经验标准。对于一个 标量的时间序列x n ，基本的想法是选择延迟t ，使得x n 和x n 竹在一定程度上不 相关，但是又不完全不相关，由此可以把它们视为重构的相空间中独立的坐标。 若t 过小，x n 和x n 竹相关性过大以至于不能作为重构的相空间中独立的变量， 在相空间中表现为重构出来的相图集中在主对角线周围，系统的很多状态区分 不出来，这种情况称为冗余( r e d u l l d a n c e ) 。如果t 过大的话，x n 和x m 就会变 成随机的变量，这在相空间中表现为相图过分的分散，没有表现出该有的固定 结构，因而原来系统潜在的确定的动力学性质将会丢失，这种情况称为不相关性 ( i 仃e l e v a i l c e ) 。因此要选择一个比较合适的延迟时间t 。目前，可以使用的方法 非常多，下面对它们进行简单介绍和比较。 2 3 1 嵌入延迟计算方法比较 嵌入延迟t 的选取存在多种方法，比如自相关函数( a u t o c o r r e l a t i o n 矗m c t i o n ) ( l u i sa a ，19 9 5 ) 、互信息( m u t u a li n f o n t l a t i o n ) ( a m f r a s e r ，h l s 埘肌e y ， 1 9 8 6 ) 、重构膨胀( r e e o n s t m “o ne x p a n s i o n ) ( m i c h a e lt r o s e n s t e i n ，j 锄e sj ，c o l l i t l s ， c a r l oj d el u c 如19 9 4 ) 和c c 方法( h s k i m ，r e y k h o l ta n dj d s a l a s ，19 9 9 ) 等。 通过调研对几种方法的特点总结如表2 1 以便于比较。 第二章1 f 线性时间序列分析方法的比较 表2 1 四种嵌入延迟计算方法特点比较 性质白相关函数法互信息法重构膨胀法c c 方法 是否包含时间序列否足是是 非线性特征 足否适用于各种非 否 是否是 线性系统 所需时间序列数据 小大小小 量人小 计算复杂程度简单很复杂( 并有复杂简单简单 窄间划分要求) 足否易于扩展到高 否 是是足 维 抗噪声能力弱弱弱强 足甭能计算时问窗含否否能 口 准确求取时间延迟的方法应该具有下列优点： a ) 概念清晰，在相空间几何学中有明确的意义。 b ) 数学表达式简单易算。 c ) 对数据流要求不大，较小数据量求解结果的可靠性仍较高。 d ) 对含噪数据集合计算结果的鲁棒性较高。 e ) 易于扩展到高维。 通过比较发现c c 方法具有较多优点，因此选择使用该方法计算时间延迟。 2 3 2c c 方法计算时间延迟 1 9 9 9 年，k i m ，e y k h o l t 和s a l a s 提出了c c 方法，该方法应用关联积分能 够同时估计出延迟时间龟和嵌入窗宽k 。设x ( n ) ，n _ 1 ，2 ，n 为时间序列， x ( n ) = x ( n ) ，x ( n 州，x 【n + ( m 一1 ) 司) ( i = l ，2 ，m ) 为相空间中的点。c c 方法的 具体描述如下： 嵌入时间序列的关联积分定义为下式的函数，其中r 0 c ( 删 加壶。s 蝥( ，一i iz 一舢 ( 2 6 ) 式中：m 为嵌入维数；n 为时间序列长度；r 为邻域半径的大小；f 为时间 延迟；m = n - ( m - 1 ) t ；0 ( ) 为h e a v i s i d e 单位函数。 ，fo ，x _ 0 秒( x ) = 1 一一 ( 2 7 ) l l ，z 0 第二章非线性时间序列分析方法的比较 | i i l 是欧几里德距离的某种量度，定义为： 慨( ，z ) y ，( 朋) l j = o 骤。k 机一札r ( 2 8 ) 将时间序列x ( n ) ，n = 1 ，2 ，n ，分成t 个不相交的时间序列，长度为 i n t m ) ，i n t 为取整。对于一般的自然数t ，有 x ( 1 ) ，x ( t + 1 ) ，x ( 2 t + 1 ) ，) x ( 2 ) ，x ( t + 2 ) ，x ( 2 t + 2 ) ， x ( t ) ，x ( t 斗i t ) ，x ( 2 t + t ) ， 然后计算每个子序列的统计量s ( m ，n ，r ，百) 为 s ( 肌，r ，f ) = c f ( 肌，厂，) _ 【c f ( 1 ，f ，r ) 】” ( 2 9 ) ，= l 式中：c i 是第1 个子序列的相关积分。 局部最大间隔可以取s ( ) 的零点或对所有的半径r 相互差别最小的时间点。 选择对应值最大和最小的两个半径r ，定义差量为 丛( 朋，f ) = m a x 降( 删，f ) 】一m i n 陋( m ，f ) 】，f ， ( 2 1 0 ) 根据统计学原理，m 取值在2 到5 之间，r 的取值在以和2 0 之间，n 大 于等于5 0 0 ( n = 3 0 0 0 时比较好) 。o 是时间序列的均方差。得到方程如下： 豆( f ) = 去耄褰s ( 肌，。，力 ( 2 1 1 ) 其中，蜃( f ) 为所有子序列的统计量s ( m ，n ，i ：j ，t ) 的均值，歹( f ) 的第一个极小 值对应第一个局部最大时间f ，s 。o r ( t ) 的最小值对应时间序列独立的第一个整体 最大值时间窗口，即延迟时问窗口。 图2 2 是用c c 算法计算l o r e n z 系统( e n l o r e n z ，1 9 6 3 ) x 变量时间序列得 到的心( f ) ，曲线图。 第二章非线性时间序列分析方法的比较 。口，1 口 蛤 口口8 5 口 图2 2c ，c 算法计算l o r e i l z 系统x 变量时间序列得到心( f ) f 曲线 从结果中可以发现，曩( ，) 的第一个极小值对应第一个局部最大时间下为 1 0 ，即该时问序列嵌入延迟为1 0 ，这符合理论上的结果。 2 4 嵌入维数计算方法选择 关于嵌入维数，f 毗e n s ( 1 9 8 1 ) 、r m a n e ( 1 9 8 1 ) 、t s a u e “1 9 9 1 ) 等先后从理 论上证明了当m 2 2 d + 1 时可获得一个吸引子的嵌入，其中d 是吸引子的分形维 数，但这只是一个充分条件，对观测时间序列选择m 没有帮助。如果仅仅是计 算关联维数，m z d i n g ( 1 9 9 3 ) 证明了对无噪声、无限长的时间序列，只要m 为 大于关联维数d 的最小整数即可，但对长度有限且具噪声的时间序列，m 要比 d 大得多。如果m 选得太小，则吸引子可能折叠以致在某些地方相交，这样在 相交区域的一个小邻域内可能会包含来自吸引予不同部分的点；如果m 选得太 大，理论上是可以的，但在实际应用中，随着m 的增大会大大增加吸引子的几 何不变量( 如关联维数、l y a p u n o v 指数等) 的计算工作量，且噪声和舍入误差 的影响会大大增加。目前，可以使用的方法很多，下面对它们进行简单介绍和 比较。 2 4 1 嵌入维数计算方法比较 嵌入维数m 的选取有许多方法，如根据吸引子计算一些不变量的试算法、 第二章非线性时间序列分析方法的比较 真实矢量场法( d a n i e lt k a p l a n ，l e o ng l a s s ，1 9 9 2 ) 、伪最近邻居法( m b k e 加e l ， r b r o w n ，h d i a b a r b a n e l ，1 9 9 2 ) 和“a n g y u ec a 0 ( 1 9 9 7 ) 提出的改进的伪最近邻 居法。通过调研把几种方法的特点列为表4 2 ，以便于大家分析比较。 表4 2 四种嵌入维数计算方法特点比较 性质试算法 真实矢量场法伪最近邻居法改进的伪最近邻居 法 理论基础 不断测试系统几何特性系统几何特性系统几何特性 计算复杂程度复杂复杂 相对简单相对简单 所需时间序列数据根据具体方法而定大 小小 量人小 足否具有主观性是( 主观选择合适是( 主观判断甲均是( 主观确定闽值， 否 ( 需要靠主观判断 嵌入维数)方向矢指标随嵌入而选择合适闽值非 才能得到结果) 维变化足否收敛)常困难) 是否需要预先计算根据具体方法而定 否是是 嵌入延迟 抗噪声能力 根据具体方法而定较强较强较强 是否易于扩展到高根据具体方法而定 否是是 维 能否区别1 f 线性信 否否否能够 号与完伞随机信号 分析比较表4 2 中四种嵌入维数计算方法特点之后，伪最近邻居法和 l i a i l g ) ，u ec a o 提出的改进的伪最近邻居法具有较多优点，而改进的伪最近邻居 法克服了伪最近邻居法中需要主观选择阈值的缺点，且能够区分非线性时间序 列和完全随机时间序列。通过综合分析比较，最终决定选择l i a n g ”ec a o 提出 的改进的伪最近邻居法来计算嵌入维数。 2 4 2l i a n g y u ec a o 算法计算嵌入维数 l i a i l g y u ec a o 提出的改进的位邻居方法主要步骤介绍如下： 对于输入的时间序列x ( n ) ，n = 1 ，2 ，n ，选定嵌入延迟下( 通过前面的c c 算法计算得到) ，对于某个维数d ，可以重构出向量序列 y i ( d ) ) ，y i ( d ) 表示d 维 重构所得到的第i 个向量，即y i ( d ) = x ( i ) ，x ( i 均，x 【i + ( d 一1 ) 司) 。这里定义： 州= 钤筹 亿 第二章= f 线性时间序列分析方法的比较 这里忆0 是欧几里德距离的某种量度，定义为： i 协( 历) 一y ，) l i = 。燃。k 小一x ， i ( 2 1 3 ) n ( i ，d ) ( 1 n ( i ，d ) 剑一d t ) 是一个整数，在定义的距离1 1 i | 的意义下，它满足 y n ( i ，d ) ( d ) 是y i ( d ) 最近的邻居。这罩需要注意的有： 1 ) 式子( 2 1 0 ) 的分子和分母的n ( i ，d ) 是一样的； 2 ) 如果y n ( i ，d ) ( d ) 等于y i ( d ) ，那么本文取距y i ( d ) 第二近的邻居作为n ( i ，d ) ，也 就是要保证方程( 2 1 0 ) 的分母不为零。 接下来定义a ( i ，d ) 的平均值： 以力2 志荟删) 由上式可见e ( d ) 只依赖于维数d 和延迟t 。为了研究维数从d 到d + 1 出现的 变化，定义： 剧( 炉等 ( 2 1 5 ) 如果时间序列来自某个吸引子，当d 大于某个值d 0 时，e l ( d ) 就会停止变化。 那么d o + 1 便是最小嵌入维数。 图2 3 是用l i a n g y u ec a 0 算法计算l o r e n z 系统( e n l o r e n z ，1 9 6 3 ) x 变量时 间序歹0 得到的e 1 ( d ) d 曲线图。 d 图2 - 3l i a n g y u ec a o 算法计算l o r e n z 系统x 变量时间序列得到e 1 ( d ) d 曲线 从图中可见，当d 大于2 时( 即d 大于等于3 时) ，e 1 ( d ) 值变化很小，可见 该时间序列的嵌入维数为3 。 第二章非线性时间序列分析方法的比较 此外，有必要定义另外一个量，它对于区分确定性的信号和随机信号是有 用的。令 e + ( d ) = 寺k 矿m 出i ( 2 1 6 )、7 一d f 鲁一”。 ”v 一扩“ 、 这里n ( i ，d ) 的意义与上面讨论的一致。同时还定义 鳓= 等 利用e 2 可以把有确定结构的信号与随机信号清楚的区分开来。对于随机的 时间序列，因为将来的值和过去的值互相独立，因此不管d 取任何值，e 2 ( d ) 总 是等于1 ( 当然，实际情况中不免有少许偏差) 。但是，对于确定性的时间序列， e 2 ( d ) 与d 有关，结果，对所有的d 而言，e 2 ( d ) 不是一个恒定值。换句话说，必 定存在某个d ，满足e 2 ( d ) 不等于1 。 2 5 关联维数计算方法选择 2 5 1 关联维数计算算法的选择 g r a s s b e 玛e r 和p r o c a c c i a 提出的关联维数算法( 简称g p 算法) 由于特别 适用于实验观测数据，且算法简洁易于实现，因此获得了广泛的应用。但是g p 算法也是存在一些缺点，如要求时间序列数据量很大且不能含有噪声。而实验 数据或多或少都含有噪声，因此利用g p 算法得到的关联维数