一种基于功率谱划分的复杂度谱脑电预测和诊断方法

摘要

本发明提出了一种基于功率谱划分的复杂度谱脑电预测和诊断方法，属于脑电信号（EEG）分析和脑精神疾病预测和诊断领域，本发明主要包括：一种基于功率谱划分的脑电信号复杂度谱定义、分析和提取方法和一种非线性罗辑斯蒂（Logistic）复杂度谱参考模型构建方法。首先对脑电信号定义了一个基于功率谱划分的复杂度谱，并给出其计算方法，然后利用本发明对此映射生成的数据序列进行计算，据此建立一个脑电信号的分析复杂度谱参考模型，并分析各结构谱线序列的大小、数量和分布所反映的物理生物意义，画出该映射基于功率谱划分的复杂度谱参考空间分布模型。本发明可以对脑精神疾病作出预测和诊断。

说明书

技术领域

本发明属于脑电信号（EEG）分析和脑精神疾病预测和诊断领域，涉及一种脑电信号处 理中的分析、建模和计算方案，尤其是一种能够对脑电信号进行参考模型化非线性处理以及 对脑电信号的复杂度细微结构进行分析的方法。

背景技术

目前，脑科学的研究是科学研究领域的一个热点，脑电信号处理是研究大脑的主要手段 之一。脑电信号的处理现在共有两类方法，一类是线性分析方法，另一类是非线性方法。线 性方法属于传统信息分析方法；非线性方法属于现代信息处理方法。线性方法主要有时域分 析、频域分析、时频分析等，非线性方法主要包括非线性动力学方法、信息熵、复杂度等。 越来越多的研究证据表明：大脑是一个非常复杂的非线性系统，脑电信号是这个系统的一个 输出。虽然线性方法在脑电的研究领域有着不可替代的作用，也没有证据表明它有被取代的 可能，但是在研究中确实发现存在许多不可避免地存在缺陷。例如，大脑中的非线性关联信 息量在分析中丢失，大脑复杂性得不到充分的研究等。在此种情况下，出现了脑电的非线性 分析方法，并且快速发展起来，已经取得许多重要成果。其中，复杂度的方法就是发展起来 比较快的一种方法。

20世纪60年代，柯尔莫哥洛夫（Kolmogorov）等人把复杂性定义为：一个系统的复杂 程度与该系统的空间结构或时间序列所表示的变化行为的最小描述有关，一般称作柯尔莫哥 洛夫复杂性。由于柯尔莫哥洛夫复杂性不可计算，1976年A.Lempel和J.Ziv提出一种度量 符号序列复杂性的简单算法，称为Lempel-Ziv复杂度。直到1987年，才由卡斯帕（Kaspar） 和舒斯特（Schuster）提出了该算法的计算机实现方法。解幸幸、李舒等也在《Lempel-Ziv 复杂度在非线性检测中的应用研究》一文中认为：某事物的算法复杂度等于产生该事物的图 形结构或符号序列的最短程序长度与该图形结构或符号序列本身大小之比的极限。

脑电信号是一种复杂的非线性信号，复杂度能够客观有效地测度和评价了信号的复杂性 和混乱程度。1991年，吴祥宝和徐京华首次将脑电复杂度分析引入脑电信号的研究中，并且 通过计算已知的混沌系统（主要是Lorenz和Rossler）产生的序列的LZ复杂度，与脑电信号 的LZ复杂度对比发现脑电信号计算的LZ复杂度远远大于混沌系统序列的复杂度。顾凡及等 用LZ复杂度研究了8中精神状态下13名受试者的脑电信号，通过研究他们认为LZ复杂度 的数值与脑电信号的不同部位、和不同的精神状态有着显著的关系。

1998年陈芳、顾凡及、徐京华等提出C₀复杂度的概念和算法，其主要思想就是：把信号 分解成规则成分和不规则成分两部分，C₀复杂度就定义为不规则成分在原信号里所占的比例。 随后在2004年，蔡志杰、顾凡及、沈恩华对该算法作了进一步的改进，主要是改进了划分规 则部分与不规则部分的方法,采用幅值平方的功率谱代替了幅度谱,然后再计算功率谱的平均 值,以此作为分界线来区分规则部分与不规则部分，进而计算出复杂度值。本文即是在此工作 的基础上的进一步改进和拓展。

近年来，有关复杂度的研究许多都和混沌时间序列结合在一起，出现许多新的研究成果。 主要分为几方面：(1)侧重关注研究混沌时间序列复杂度的计算方法；(2)侧重关注研究混沌时 间序列复杂度的强度计算；(3)关注研究混沌时间序列复杂度的稳定性分析；(4)研究关注混沌 时间序列复杂度的相关性分析。所有这些研究有一个共同点就是都只关注一个单一指标的复 杂度的计算和分析，没有提出对其细微结构进行进一步的分析研究，没有对多复杂度概念指标 进行研究，没有提出复杂度谱的定义，分析和计算方法。虽然有人提出结构复杂度的概念， 但还是没有对时间序列复杂度的组成、结构成分以及数量分布的进行研究。现有的脑电时间 信号都是具有强的非线性性的信号，仅仅分析计算单一复杂度及其性质不能满足实际存在的 多种脑精神疾病并存、诊断容易失误和误诊的实际情况。因此，提出一种在整个非线性时间 在统一的背景下的从规则周期，经分叉混沌周期，再到随机序列的演化过程的复杂度或复杂 度谱的整个演变过程，并能够计算时间序列复杂度，计算各成分的数量、大小和分布的方法 和分析其性质及物理或者生物意义，对于预测和诊断脑电信号有着很重要的参考意义。

综上所述，现有方法和研究成果虽然比较多，但对于脑精神疾病的诊断和预测研究远远 不够，复杂度的研究虽然给这个研究带来较大的发展，但其单一复杂度指标现状难以满足脑 电信号强非线性的性质的实际需求，尤其对多种脑精神疾病并存的状况，给脑电信号的分析、 诊断和预测带来更大的难度。

发明内容

为解决实际技术中单一复杂度指标能够满足对与多种脑精神疾病的脑电信号的分析、诊断 和预测带来的大概率的失误和误诊的问题，本发明提供一种基于功率谱划分的复杂度谱脑电 预测和诊断方法，这种方法不仅能够对脑电信号进行分析、建模和计算，还能对脑电信号的 进行参考模型化非线性处理以及对脑电信号的复杂度细微结构进行分析计算。

本发明采用如下的技术方案。

一种基于功率谱划分的复杂度谱脑电预测和诊断方法，其特征在于包括以下步骤：

步骤一：基于功率谱划分的复杂度谱模型的定义及计算：本发明在现有技术和研究的基 础上扩展和拓宽C₀复杂度的概念，提出复杂度谱的定义和计算方法，并给出了通过功率谱划 分确定复杂度谱长度的方法。具体方法如下：

(1)基于功率谱划分复杂度谱模型的定义及计算，方法如下：

1)对于给定的脑电信号序列X＝{x(n),n＝0,1,2,...,N-1}，计算其傅立叶（Flourier）变换， 得到功率谱序列，记为F。

2)利用公式计算序列F在f(k)点处的功率谱值。得功率谱值序列: P＝{p(k),k＝0,1,2,...,N-1}。

3)寻找序列P的最大值和最小：P_max=max(P)，P_min=min(P)，计算其差值ΔP=P_max-P_min， 将ΔP平均分成m段(2≤m≤N的正整数,N为序列元素个数)，从下到上依次记为ΔP₁,ΔP₂,…, ΔP_m，记新序列为则按如下规则生成新序列：

4)对于序列{k＝0,1,2,...,N-1，j=1,2,…,m}进行傅立叶（Flourier）逆变换， 得新的时间序列为：

$\widetilde{x^j}\left(n\right)=\underset{k=0}{\overset{N-1}{\Sigma }}\widetilde{f_N^j}\left(k\right)W_N^{\mathrm{nk}},k=\mathrm{0,1,2},...,N-1,j=\mathrm{1,2},...,m---\left(2\right)$

5)定义如下求和比

$C_X\left(j\right)=\frac{\underset{n=0}{\overset{N-1}{\Sigma }}{\left|\widetilde{x^j}\left(n\right)\right|}^2}{\underset{n=0}{\overset{N-1}{\Sigma }}{\left|x\left(n\right)\right|}^2},j=\mathrm{1,2},...,m---\left(3\right)$

为C_X(j)为序列X的第j复杂度；称序列s＝{C_X(j),j＝1,2,...,m}，为脑电时间序列X的基 于功率谱划分的复杂度结构谱。

(2)确定基于功率谱划分的复杂度谱长度m，方法如下：

1)构造关于m的函数

$E_{\mathrm{CX}}\left(m\right)={\left(\underset{i=1}{\overset{m}{\Sigma }}{\left(C_X\left(j\right)\right)}^2\right)}^{\frac{1}{2}}---\left(4\right)$

作为目标函数，m是整数，2≤m≤N，N为序列元素个数。

2)作出目标函数曲线，对每一个m求其线性拟合曲线：

y＝a_mm+c       (5)

其中：a_m是取m时的斜率，c截距。得到一个斜率序列A={a₂,a₃,…,a_m}；

3)计算A中每两个相邻项之间的差值的绝对值序列：

Y={Δa_i＝|a_i+1-a_i|，i=2,3,4,…,m-1}     (6)

4)在Y中寻找项Δa_i，对于序列中l≥i的所有项，都有Δa_l≤Δa_i，i项Δa_i就是序列Y开 始收敛或稳定变化的起始点。从i项后所有序列项就是收敛和稳定变化的区域。此时，序列相 对Δa_i的稳定或收敛的变化精度为：

$\delta =\frac{\Delta a_i}{\max \left\{Y\right\}}---\left(7\right)$

其中：0<δ≤0.05，i即是相对于δ开始收敛点，也是相对于精度δ的最小划分数目。

步骤二：建立非线性罗辑斯蒂（Logistic）映射的复杂度谱参考演化模型：本发明利用步 骤一定义的方法，对经典的非线性罗辑斯蒂映射序列的数据进行计算建立复杂度参考演化模 型，主要包括计算复杂度谱和标定模型，以及复杂度计算窗口的最小长度和稳定复杂度值的 确定方法，具体方法如下：

(1)建立非线性罗辑斯蒂映射的复杂度谱参考演化模型,方法如下：

1)根据罗辑斯蒂（Logistic）映射公式，生成控制参数u∈[0.75,2]的演化区间作为建立 脑电复杂度谱从周期经混沌向随机演化的标准数据，按照步骤一的方法计算其每一时间序列 的复杂度谱：

$s^i=\{C_X^i\left(j\right),j=\mathrm{1,2},...,m\},i\in \left[\mathrm{0.75,2}\right]---\left(8\right)$

sⁱ表示i对应复杂度谱。

2)按照给定的数据步长，逐步加大计算窗口的长度并重复计算S,直到各复杂度谱中的 每个复杂度值收敛或变化稳定为止，此时计算和统计稳定区间的如下过程数据：

i)计算最小数据窗口的长度：

N_min=min{N₁,N₂,N₃,…}        (9)

N_i表示在收敛或变化稳定是各计算窗口长度；

ii)计算S在稳定区域的各个谱中每个的近似平均值：

$\bar{S}=\{\overline{C_X^i\left(j\right)},j=\mathrm{1,2},...,m,i\in [\mathrm{0.75,2}\left]\right\}---\left(10\right)$

iii)计算S中各复杂度值的变化范围：

$\Delta C_X^i\left(j\right)=\max \left\{\widetilde{C_X^i\left(j\right)}\right\}-\min \left\{\widetilde{C_X^i\left(j\right)}\right\}---\left(11\right)$

表示在脑电序列X在控制参数u＝i时第j复杂度在不同长度的窗口计算得到的复杂 度值序列。

3)根据复杂度谱序列最后计算值，并对其校验后在三维坐标图中作出三维空间模型 图，并标出变化范围的边界。作出各序列每个复杂度谱的复杂度值的窗口长度三维模型。

(2)确定最小计算窗口长度和稳定区平均复杂度值，算法如下：

1)脑电序列X在u＝i时第j复杂度在以L/30的步长(L为数据长度)来增加计算窗口长度, 并计算得到的复杂度值序列为：

$\widetilde{C_X^i\left(j\right)}=\{Z_v,Z_{v+d},Z_{v+2d},...\}---\left(12\right)$

v表示起始计算时的数据窗口长度，为1≤v≤L的正整数；d表示每次增加的计算步长， 为1≤d≤L的正整数。

2)计算每相邻两项的中位值：

$t_{v+\mathrm{id}}=\frac{1}{2}(Z_{v+\mathrm{id}}+Z_{v+(i+1)d})---\left(13\right)$

并构成新的序列：

T={t_v+id，i=0,1,2,…}         (14)

3)以新序列T代替重复迭代第2)步10次，取其最后一项作为近似平均复杂度值 t_end；

4)将序列的每一项Z_v+id与t_end做差值：

ΔZ_v+id＝Z_v+id-t_end       (15)

取序列中ΔZ_v+id≥0项的元素下标组成新的序列H_big，取序列中ΔZ_v+id≤0项 的元素下标组成新的序列H_small；

5)计算序列H_big中的最小值：

b_big=min{H_big}        (16)

计算H_small中的最小值：

b_small=min{H_small}          (17)

然后求两者中的最大值：

N₁=max{b_big,b_small}        (18)

6)在中，取下标H≥N₁的所有项组成一个新序列：

G={Z_N,Z_N+d,Z_N+2d,...}            (19)

设G有l项，计算序列G的平均值：

${\bar{Z}}_1=\frac{1}{l}\underset{l=0}{\overset{l-1}{\Sigma }}{Z}_{N+\mathrm{id}}---\left(20\right)$

7)用替代t_end，重复(3)～(5)，直到计算得到的N_i值与上一次计算的N_i-1值相同为止。 记最后一次计算的N_i为N_min，为

这样，最后计算得到的N_min和即为所求的最小计算窗口长度和稳定区平均复杂度值。

步骤三：计算各类脑电信号的复杂度谱序列并提出特征在参考模型中标定，方法如下：

(1)采集包括抑郁症在内的几种典型的脑精神疾病的脑电信号，按照步骤一和步骤二的 计算方法，计算N_min、和指标；

(2)并将这些指标标定在复杂度谱三维空间模型图和最小窗口三维模型图中，以此确定这 几类脑精神疾病脑电信号指标在这两个模型中的空间位置和分布范围的大小。

步骤四：对于临床采集的脑电信号计算其特征指标，并将其与参考模型的特征空间做相 似度对比，计算其和各类空间之间的距离dis（dis表示在特征空间的欧氏距离，为dis>0的 正实数）,最后依据dis的大小进行排序，取最小值的dis的类作为诊断输出。

本发明的有益效果是，本方法不仅能够对脑电信号进行分析、建模和计算，还能对脑电 信号的进行参考模型化非线性处理以及对脑电信号的复杂度细微结构进行分析计算。同时可 以在降低实际临床诊断中存在的多种脑精神疾病并存而导致的诊断容易失误和误诊的概率， 为医生或研究者为提供一种新的脑电信号分析和计算方法。

附图说明

图1为基于功率谱划分的复杂度谱脑电预测和诊断方法流程图；

图2为计算基于功率谱划分的脑电信号复杂度谱的流程图；

图3为确定功率谱划分m的流程图；

图4为某位抑郁症患者的基于功率谱划分的复杂度谱；

图5为随控制参数u的演化功率谱划分数m变化曲线；

图6为精度变化和划分数m关系曲线；

图7为计算演化参考模型流程图；

图8为计算最小稳定窗口长度和谱线稳定值流程图；

图9为罗辑斯蒂（Logistic）映射复杂度谱线稳定值演化参考模型图(u∈[0.75,2])；

图10为罗辑斯蒂（Logistic）映射复杂度谱线变化幅度值分布图(u∈[0.75,2])；

图11为罗辑斯蒂（Logistic）映射复杂度谱线最小稳定计算窗口长度分布图(u∈[0.75,2])。

具体实施方式

下面结合附图和实施示例对本发明作进一步说明。本实施示例针对在某医院采集的抑郁症 患者脑电信号进行分类和诊断示范过程。

图1是本发明所述的方法流程图，具体包括以下步骤：

步骤一：计算基于功率谱划分的脑电信号复杂度谱，方法如图2所示，包括以下内容：

(1)根据公式(1)～(3),对脑电信号计算基于功率谱划分的复杂度谱序列。图4是抑郁症患者 的功率谱划分为36时的脑电复杂度谱，其中m值是根据(2)步按照精度0.02确定的。

(2)根据图3所示，按照公式(4)～(7)，计算不同精度下的功率谱划分需要的最小复杂度的 长度。本发明已经计算得出了罗辑斯蒂（Logistic）映射在0.75～2之间序列演化时，精度为 0.5～0.01之间需要的最小划分数m，如图5所示。图6为随计算精度变化所需的最小划分数 m变化曲线。

步骤二：建立非线性罗辑斯蒂映射的复杂度谱演化参考模型，方法如下：

(1)根据公式(8)，按照图7的流程计算不同窗口的复杂度值，知道稳定为止。构成复杂 度谱序列。

(2)根据公式(9)～(20)，按照图7的流程计算最小稳定窗口长度和稳定复杂度谱，以及谱线 的变化幅度。

(3)根据计算的稳定复杂度谱值绘制复杂度参考模型图，如图9所示；根据谱线值变化幅 度绘制参考模型变化分布图，如图10所示；根据最小稳定窗口长度绘制参考模型最 小窗口变化分布图，如图11所示。

步骤三：计算各类脑电信号的复杂度谱序列并提出特征在参考模型中标定，方法如下：

(1)利用步骤一计算多种脑疾病脑电信号的复杂度谱，根据复杂图谱距离参考模型的谱线 的欧氏距离排序，将其标定在参考模型里距其最小谱线所在的区域。

(2)逐步积累大量的样本点，最后在参考模型理形成不同脑疾病模型区域。通过这些区域 分析空间分布关系，进而分析其内在联系、演化关系和物理生物意义。

步骤四：对临床脑电信号，计算其基于功率谱划分的复杂度谱，并将其与参考模型空间 的各类疾病区域做相似度对比，计算其和各类空间之间的距离dis,最后依据dis的大小进行 排序，取最小值的dis的类作为诊断输出。例如：图4的复杂度谱序列s₁在演化模型空间中， 经计算得到与癫痫症特征空间、与抑郁症特征空间以及与正常健康人特征空间的综合欧氏距 离分别为0.15、0.05和0.3。对这三个数进行排序，最小者为0.05，它对应着抑郁症特征空 间。据此，输出的诊断结果为：s₁为抑郁症患者的脑电复杂度谱。