基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频分析方法
基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频分析方法
【技术领域】
[0001] 本发明涉及一种脑电信号时频分析方法，尤其涉及一种基于多尺度径向基函数和 改进粒子群优化算法的时频分析方法，属于信号分析与处理技术领域。
【背景技术】
[0002] 癫痫是由大脑神经细胞反复超同步放电，引起的自发性、突发性的脑功能紊乱的 神经系统疾病，目前癫痫治疗的主要手段一手术切除病灶和药物治疗存在一些问题，对有 些患者可能带来并发症和一些不良反应，如果能在癫痫发作前进行早期确诊，提前采取保 护措施，就可极大的减少患者受伤的风险，并可对认识癫痫的发病机制和研宄新的治疗方 法起到促进作用。脑电图（Electroencephalogram，EEG)是一种评估脑的状态和疾病的重 要形式。通过脑电图检查并及时分析脑电信号，从癫痫患者的脑电信息中提取出能够反映 大脑功能状态的特征参数，提前对这些人进行临床干预。因此，脑电图检查并提取脑电信 号特征对于癫痫的诊断和研宄有着重要的意义。
[0003] 脑电信号是一种非平稳信号。由于非平稳信号的频率随时间变化，需对其进行时 频分析，即描述信号频谱含量在时间上分布的分析方法，提取其时频特征。脑电信号按频率 可划分为 4 个波段：S波（〇? 5-4HZ)、0 波（4-8HZ)、a波（8-13HZ)、0 波（13-30HZ)。
[0004]目前时频分析方法主要有两类：非参数方法和参数方法。非参数方法如短时傅里 叶变换是基于对信号时频联合分布的非参数表示，根据海森堡测不准原理，该方法的主要 缺点是时间分辨率和频率分辨率不能同时达到最优。参数方法对信号进行建模及参数进行 估计，能够同时得到较高的时频分辨率。
[0005] 时变参数模型方法的主要任务是对其时变参数进行辨识。目前主要采用的基函 数扩展法。基函数扩展法的主要思想将时变参数表示为一组已知基函数的线性加权组合 (陈宇，陈怀海，李赞澄，等.基于时变AR模型和小波变换的时变参数识别[J].国外 电子测量技术，2011，30(7) : 20-23.)，将时变问题转化为关于基函数的时不变参数辨识问 题，通过对时不变参数的辨识进而得到时变参数。目前可供选择的基函数有傅里叶基、径 向基函数（LiY,WeiHL,BillingsSA,etal.Time-varyingmodelidentificationfor time-frequencyfeatureextractionfromEEGdata[J] ?JournalofNeuroscience Methods, 2011，196:151 - 158.)等。每种基函数都有各自逼近特性，如傅里叶基函数和勒让 德多项式可以有效辨识变化缓慢且平滑的时变参数，而径向基函数可以同时辨识平滑及变 化剧烈的时变参数。因而本发明采用径向基函数进行时变参数辨识。
[0006] 脑电信号是一种微弱的生物信号，因其噪声大，能量低及非平稳强等特点，因而很 难提取潜在的具有生物特性的瞬时特征。在此背景下，研宄一种基于时变参数建模的时频 分析方法，并采用径向基函数辨识时变参数能同时获得较高的时频分辨率，能够有效提取 脑电时频特征，对于癫痫脑电信号的准确时频特征提取及辅助癫痫疾病的诊断具有重要意 义。
【发明内容】
[0007] 本发明提供了一种基于时变参数模型的脑电时频分析方法，采用多尺度径向基函 数（Multi-scaleradialbasisfunction，MRBF)对时变参数进行展开，将时变参数估计 问题转化为时不变模型，然后对时不变参数进行辨识。其中，径向基函数具有多尺度和多分 辨率特性，对变化快速和缓慢的时变参数都能有效识别与跟踪，已广泛应用于具有多种动 态特征的时变参数辨识。径向基函数的最优尺度由改进粒子群优化算法（Particleswarm optimization，PSO)决定。PS0算法是一种群体智能算法，能够根据启发性信息一适应度值 在众多粒子中寻到最优解（李秀英，韩志刚.一种基于粒子群优化的非线性系统辨识方法 [J].控制与决策，2011，26(11) =1627-1631.)。本发明通过对癫痫脑电信号和正常脑电信 号进行时频特征提取分析，准确提取了癫痫脑电信号和正常脑电信号的时频特征，能够对 癫痫脑疾病的辅助诊断提供定量的技术支持。
[0008] 为实现上述目的，本发明提供了：
[0009] -种基于多尺度径向基函数和改进粒子群优化算法的非平稳信号时频分析方法， 包括如下步骤：
[0010] 1 ?时变参数建模；
[0011] 2.时变参数展开，时变模型转化为时不变参数模型；
[0012] 3.选择径向基函数的中心及尺度；
[0013] 4?时不变参数估计；
[0014] 5?时变参数估计；
[0015] 6.时频分析，由时变参数估计值求信号时频分布特征。
[0016] 其中，在所述步骤3中，由改进粒子群优化算法选择径向基函数的最优尺度。
[0017] 所述步骤3中，改进粒子群优化算法中速度更新公式的惯性权重为与适应度值有 关的变化值。
[0018] 所述步骤3中，衡量模型拟合效果的相关指数作为粒子的适应度值。
[0019] 本发明所提供的基于多尺度径向基函数和粒子群优化算法的非平稳信号时频分 析方法的优点包括：
[0020] 1.时频分辨率高，可准确辨识强噪声、非平稳、非线性的信号，准确提取其时频特 征；
[0021] 2.能根据信号灵活选择径向基函数尺度，适应性好；
[0022] 3.准确提取癫痫脑电信号与正常脑电信号时频特征，为癫痫脑疾病的辅助诊断提 供定量的技术指导。
【附图说明】
[0023] 图1为根据本发明实施例的时频分析方法的流程示意图。
[0024] 图2为粒子群优化方法选择多尺度径向基函数最优尺度的流程示意图。
[0025] 图3(a)显示了真实癫痫脑电信号。图3(b)_3(c)显示了根据本发明实施例的时 频分析方法与经典时频分析方法对真实癫痫脑电信号时频分析结果对比；其中图3(b)为 短时傅里叶方法的时频分析结果，图3(c)为本发明方法的时频分析结果。
[0026] 图4(a)显示了真实正常人脑电信号。图4(b)_4(c)显示了根据本发明实施例的 时频分析方法与经典时频分析方法对真实正常人脑电信号的时频分析结果对比；其中图 4(b)为短时傅里叶方法的时频分析结果，图4(c)为本发明方法的时频分析结果。
【具体实施方式】
[0027] 下面结合附图和【具体实施方式】对本发明作进一步的详细说明。
[0028] 真实脑电信号的时频分析作为本发明的实施例验证本发明的时频分析效果，并将 本发明与经典时频分析方法进行对比。真实脑电信号的时频分析包括对癫痫脑电信号和正 常人脑电信号的时频分析比较。
[0029] 实验数据来自德国波恩癫痫研宄室临床采集的脑电数据库，采样频率为173. 6Hz， 持续时间为23. 6秒。如图3(a)和图4(a)所示，分别为癫痫脑电信号和正常人脑电信号。
[0030] 图1展示了根据本发明实施例方法的流程图，包括：
[0031] 首先对信号进行时变参数建模，由最终预测误差（Finalpredictionerror,FPE) 模型定阶准则对模型进行定阶（步骤1);接着将时变参数用MRBF展开，表示为时不变参数 形式（步骤2);然后选择MRBF的中心并由PS0算法选择MRBF的最优尺度（步骤3);再对 时不变参数估计，将模型表示为矩阵形式，对时不变参数进行辨识（步骤4)，得到时变参数 估计值；由时不变参数估计值和MRBF估计时变参数（步骤5);最后进行时频分析，由时变 参数计算时频谱值，得到时频分布图（步骤6)。
[0032] 下面具体介绍根据本发明所提供的基于多尺度径向基函数和改进粒子群优化算 法的时频分析方法的具体方法：
[0033] 1.时变参数建模
[0034]时变参数模型参数是随时间变化的，一个p阶时变参数模型如下所示：
[0036] 式中：p为模型阶次，ai (n)(i= 1，2,…，p)为时变参数，e(n)是均值为0,方差为 0 2的尚斯白噪声。
[0037] 模型定阶：本发明采用最终预测误差（FPE)定阶准则选择模型阶数，其表达式如 下：
[0039] 式中，N为数据长度，p为模型阶数，&为预测误差。
[0040] 2.时变参数展开，将时变参数表示为一组多尺度径向基函数的加权线性组合，如 下式所示：
[0042]式中，ci;ni为基函数的时不变权系数，fjt)为高斯径向基函数。M为基函数的维 数。将上式代入式（1)得：
[0044] 式（4)即为时不变参数模型。
[0045] 3.径向基函数中心和尺度的确定
[0046] 本发明采用高斯径向基函数作为基函数扩展时变参数。由于高斯径向基函数具有 中心径向对称、逼近性能强、效果好，学习速度快等优点，且受信号特征的影响较小，所以本 发明选择高斯径向基函数作为基函数对时变参数进行展开，其表达式如下：
[0048]式中c为径向基函数的中心，〇 2为基函数的尺度，它决定了基函数围绕中心点的 距离。X为输入样本。在径向基函数作为基函数时，x(x= 1,2,…，乂?^为采样长度）为离 散采样时间序列。
[0049] 选择合适的径向基函数中心和尺度是准确辨识时变参数的关键。为了使径向基函 数分布到整个时变参数中，以确保径向基函数对时变参数的所有局部进行准确估计，本发 明将径向基函数的中心均匀分布到时变参数中，设第k个径向基函数及径向基函数的中心 分别表不为：
[00 52] 式中，M为径向基函数的维数，x为样本输入，ck为第k个径向基函数的中心，(J2k 为第k个径向基函数的尺度，它决定了该基函数围绕中心点的距离，N为采样数据长度。
[0053] 在本发明中，由改进粒子群优化算法计算径向基函数的最优尺度。
[0054] 图2为改进粒子群优化算法计算径向基函数的最优尺度流程图，具体步骤如下：
[0055] (a)初始化：初始化粒子、粒子个数、适应度值、局部最优粒子、全局最优粒子、粒 子最大值和最小值、速度最大值、迭代次数等。
[0056] 径向基函数的候选尺度为：
[0058] 式中，N为采样数据长度，M为径向基函数的维数，sk为需要进行调节的任意整数， 且最大值在10以内，粒子为Sk构成的向量为ui=[ss2，…，sM]。
[0059] 取随机整数为粒子赋初值，初始适应度值为0,粒子数和迭代次数可根据参数变化 特征来调整。
[0060] (b)计算每个粒子的适应度值
[0061] 粒子的适应度值是帮助找到最优粒子的启发性信息，在本发明中，采用模型拟合 效果的相关指数作为适应度函数，相关指数公式如下：
[0063] 式中，R2为相关指数，N为采样数据长度，；f为数据Y的模型预测输出。相关指数 越大说明模型拟合效果越好。
[0064] 将粒子代入尺度公式（8)，得到径向基函数的尺度，再根据径向基函数中心，得到 多尺度径向基函数。再由多尺度径向基函数展开方法对时变系统的模型参数辨识，得到测 量数据Y的模型预测值最后根据式（9)计算相关指数，得到粒子的适应度值。
[0065] (c)更新局部最优粒子。局部最优粒子为当前循环中适应值最大的粒子，找到适应 度值最大的粒子，该粒子作为局部最优粒子。
[0066] (d)更新全局最优粒子。全局最优粒子是所有循环中适应度值最大的粒子。若局 部最优粒子适应度值比全局最优粒子大，则赋值给全局最优粒子。否则全局最优粒子不变。
[0067] (e)更新每个粒子的速度，按粒子速度更新公式更新每个粒子的速度。
[0068] 粒子速度更新公式如下：
[0070]其中，Cbest为局部最优粒子，gbest为全局最优粒子，#>为第h次循环的第i个粒 子，V广为第h次循环的第i个粒子的速度，Ct为第h次循环的局部最优粒子，gt为第 h次循环后当前全局最优粒子。0为0~1之间的随机数，(^和(32的作用是防止粒子陷入 局部最优：
[0073] 式中，H为最大循环次数，h= 1，2, ???，!!?
[0074] 粒子速度更新公式中，《为惯性权重。在传统粒子群优化算法中，《 -般为常数， 不利于灵活寻找最优粒子。本发明中，采用一种改进粒子群优化算法，令惯性权重值与适应 度值相联系，以便更灵活寻找最优粒子，其表达式如下：
[0076] 式中，0是0~1之间的随机数，Fi是第i个粒子的适应度值，Fbest是全局最优粒 子的适应度值。
[0077] 当速度接近0时，按下式重新为粒子速度赋值：
[0079] 式中，
中的第j个元素，y是常值〇. 1，9是〇~1之间的随机数。
[0080] (f)更新每个粒子的位置。粒子位置更新公式如下：
[0082] (g)返回步骤（b)，重复执行步骤（b)、（c)、（d)、（e)和（f)，直到得到最优粒子。 最后得到全局最优粒子即为最优粒子，进一步根据尺度公式得到径向基函数的最优尺度。
[0083] 4.时不变参数估计。用最小二乘法对时不变参数模型中的时不变参数进行估 计，首先将时不变模型写为矩阵形式。
[0084]将式⑷中参数做如下定义：
[0085]f(t)= f2(t),???, fM(t)],
[0086] Hi(t) = y(t-i)f(t),
[0087]H(t)= H2(t),???, Hp(t)],
[0088] (；=[ci;1,ci;2，…，ci;M]，
[0089] C= [Q,^,-,Cp],
[0090] 则将式（4)写成矩阵形式：
[0091] y(t)=H(t)CT+e(t) (16)
[0092] 式中，T为矩阵的转置。C由最小二乘算法求得，即
[0094]式中，H= [H(1)，H(2)，…，H(N)，]Y= [y(l)，y(2)，...，y(N)]。
[0095] 5.时变参数估计。得到时不变参数估计值后，将时不变参数和多尺度径向基函数 代入式（3)即可得到时变参数的估计值。
[0096]6.时频分析。根据估计的时变参数，应用功率谱公式对信号进行时频分析，其计算 公式如下：
[0098] 式中，fs为采样频率，4⑴为时变参数的估计值，< 为观测误差的方差。
[0099] 求出信号的时频分布图，即得到信号的时频分布特征。
[0100] 本发明方法即基函数扩展法与经典时频分析方法即短时傅里叶变换方法进行比 较，对真实癫痫脑电信号及正常脑电信号时频结果进行分析，揭示脑电信号潜在的瞬时时 变特征。时频分析结果分别如图3(b)-3(c)和图4(b)-4(c)所示。
[0101] 由癫痫脑电信号和正常脑电信号时频分析结果可知，两种方法得到了较为一致 的频率特征分布，但显然短时傅里叶变换方法的时频分辨率低，结果如图3(b)和图4(b)所 示，时频图中频率分布不清晰。而本发明方法的时频分辨率高，结果如图3(c)和图4(c)所 不〇
[0102] 由本发明方法的时频图可清晰得到癫痫脑电信号和正常脑电信号的时频分布特 征。特别地，在癫痫脑电信号的时频分布图中，能够清晰地观察两种频段波：0波和0波。 如图3 (c)所示，0波（4-8Hz)分别在0~6秒之间、9~11秒及14~22秒之间出现；0波 (13-30HZ)出现在1~4秒、及12~22秒之间的时间段。类似地，在正常脑电信号的时频分 布图中，可以观察到两种不同频段的脑电波：S波和a波。如图4(c)所示，S波（〇-4Hz) 在4~5秒、10~12秒、14~16秒及18~22秒之间的时间段出现，a波（8-12Hz)主要 出现在17~20秒及14~15秒之间的时间段。
[0103] 由真实脑电信号的时频分析结果可知，与经典短时傅里叶变换时频方法相比，本 发明方法能够准确得到真实脑电信号的时频分布特征。由癫痫脑电信号和正常脑电信号时 频分析结果对比可知，正常脑电信号在整个时间段上，发生频率较低的脑电活动，而癫痫脑 电信号的脑电频率活动频率明显增大。而且，癫痫患者发作事件均被检测到，被检测到的发 作事件，尽管发作持续时间很短（2s)，本发明方法仍能检测到，并且所提取的脑电特征在 癫痫发作期具有明显的变化。本发明方法准确检测到癫痫发作的脑电信号中存在与正常脑 电信号不同频率的脑电波，而且提供了不同频率的瞬时频率分布信息，这些信息可以为医 生对癫痫脑疾病发作的辅助诊断提供定量的分析工具。
[0104] 本发明所提供的基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频 分析方法，主要是为准确提取脑电信号瞬时时频分布特征，辅助癫痫脑疾病的诊断提出的。 但显然，本说明书中所描述的时频分析方法也适用于其他非平稳信号的时频分析。
[0105] 以上对本发明所提供的基于多尺度径向基函数和改进粒子群优化算法的脑电信 号时频分析方法进行了详细的说明，但显然本发明的范围并不局限于此。在不脱离所附权 利要求书所限定的保护范围的情况下，对上述实施例的各种改变都在本发明的范围之内。
【主权项】
1. 一种基于多尺度径向基函数与改进粒子群优化算法的脑电信号时频分析方法，其特 征在于包括： 步骤1.时变参数建模； 步骤2.时变参数展开，时变模型转化为时不变参数模型； 步骤3.选择径向基函数的中心及尺度； 步骤4.时不变参数估计； 步骤5.时变参数估计； 步骤6.时频分析，由时变参数估计值求信号时频分布特征。2. 如权利要求1所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，由改进粒子群优化算法选择径向基函数的最优尺度。3. 如权利要求2所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，改进粒子群优化算法中速度更新公式的惯性权重为与适应度值有 关的变化值。4. 如权利要求2所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，衡量模型拟合效果的相关指数作为粒子的适应度值。
【专利摘要】本发明提出了一种基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频分析方法。该方法应用时变参数建模法对脑电信号进行时频特征提取分析，引入多尺度径向基函数展开式辨识时变参数。首先建立时变参数模型，然后用基函数扩展法辨识模型时变参数，即将时变参数表示为一组多尺度径向基函数的线性加权组合，将时变参数的辨识问题转化为时不变参数的辨识。径向基函数的最优尺度由粒子群优化算法决定。最后根据时变参数估计值和功率谱密度公式计算脑电信号的时频分布特征。与现有时频分析方法相比，本发明方法能同时获得较高的时间和频率分辨率，能够准确提取脑电信号的时频分布特征，对癫痫脑电信号的应用及辅助癫痫疾病的诊断具有重要意义。
【IPC分类】G06F19/00, G06K9/00
【公开号】CN104899436
【申请号】CN201510284547
【发明人】李阳, 刘青, 王旭东
【申请人】北京航空航天大学
【公开日】2015年9月9日
【申请日】2015年5月29日
【技术领域】
[0001] 本发明涉及一种脑电信号时频分析方法，尤其涉及一种基于多尺度径向基函数和 改进粒子群优化算法的时频分析方法，属于信号分析与处理技术领域。
【背景技术】
[0002] 癫痫是由大脑神经细胞反复超同步放电，引起的自发性、突发性的脑功能紊乱的 神经系统疾病，目前癫痫治疗的主要手段一手术切除病灶和药物治疗存在一些问题，对有 些患者可能带来并发症和一些不良反应，如果能在癫痫发作前进行早期确诊，提前采取保 护措施，就可极大的减少患者受伤的风险，并可对认识癫痫的发病机制和研宄新的治疗方 法起到促进作用。脑电图（Electroencephalogram，EEG)是一种评估脑的状态和疾病的重 要形式。通过脑电图检查并及时分析脑电信号，从癫痫患者的脑电信息中提取出能够反映 大脑功能状态的特征参数，提前对这些人进行临床干预。因此，脑电图检查并提取脑电信 号特征对于癫痫的诊断和研宄有着重要的意义。
[0003] 脑电信号是一种非平稳信号。由于非平稳信号的频率随时间变化，需对其进行时 频分析，即描述信号频谱含量在时间上分布的分析方法，提取其时频特征。脑电信号按频率 可划分为 4 个波段：S波（〇? 5-4HZ)、0 波（4-8HZ)、a波（8-13HZ)、0 波（13-30HZ)。
[0004]目前时频分析方法主要有两类：非参数方法和参数方法。非参数方法如短时傅里 叶变换是基于对信号时频联合分布的非参数表示，根据海森堡测不准原理，该方法的主要 缺点是时间分辨率和频率分辨率不能同时达到最优。参数方法对信号进行建模及参数进行 估计，能够同时得到较高的时频分辨率。
[0005] 时变参数模型方法的主要任务是对其时变参数进行辨识。目前主要采用的基函 数扩展法。基函数扩展法的主要思想将时变参数表示为一组已知基函数的线性加权组合 (陈宇，陈怀海，李赞澄，等.基于时变AR模型和小波变换的时变参数识别[J].国外 电子测量技术，2011，30(7) : 20-23.)，将时变问题转化为关于基函数的时不变参数辨识问 题，通过对时不变参数的辨识进而得到时变参数。目前可供选择的基函数有傅里叶基、径 向基函数（LiY,WeiHL,BillingsSA,etal.Time-varyingmodelidentificationfor time-frequencyfeatureextractionfromEEGdata[J] ?JournalofNeuroscience Methods, 2011，196:151 - 158.)等。每种基函数都有各自逼近特性，如傅里叶基函数和勒让 德多项式可以有效辨识变化缓慢且平滑的时变参数，而径向基函数可以同时辨识平滑及变 化剧烈的时变参数。因而本发明采用径向基函数进行时变参数辨识。
[0006] 脑电信号是一种微弱的生物信号，因其噪声大，能量低及非平稳强等特点，因而很 难提取潜在的具有生物特性的瞬时特征。在此背景下，研宄一种基于时变参数建模的时频 分析方法，并采用径向基函数辨识时变参数能同时获得较高的时频分辨率，能够有效提取 脑电时频特征，对于癫痫脑电信号的准确时频特征提取及辅助癫痫疾病的诊断具有重要意 义。
【发明内容】
[0007] 本发明提供了一种基于时变参数模型的脑电时频分析方法，采用多尺度径向基函 数（Multi-scaleradialbasisfunction，MRBF)对时变参数进行展开，将时变参数估计 问题转化为时不变模型，然后对时不变参数进行辨识。其中，径向基函数具有多尺度和多分 辨率特性，对变化快速和缓慢的时变参数都能有效识别与跟踪，已广泛应用于具有多种动 态特征的时变参数辨识。径向基函数的最优尺度由改进粒子群优化算法（Particleswarm optimization，PSO)决定。PS0算法是一种群体智能算法，能够根据启发性信息一适应度值 在众多粒子中寻到最优解（李秀英，韩志刚.一种基于粒子群优化的非线性系统辨识方法 [J].控制与决策，2011，26(11) =1627-1631.)。本发明通过对癫痫脑电信号和正常脑电信 号进行时频特征提取分析，准确提取了癫痫脑电信号和正常脑电信号的时频特征，能够对 癫痫脑疾病的辅助诊断提供定量的技术支持。
[0008] 为实现上述目的，本发明提供了：
[0009] -种基于多尺度径向基函数和改进粒子群优化算法的非平稳信号时频分析方法， 包括如下步骤：
[0010] 1 ?时变参数建模；
[0011] 2.时变参数展开，时变模型转化为时不变参数模型；
[0012] 3.选择径向基函数的中心及尺度；
[0013] 4?时不变参数估计；
[0014] 5?时变参数估计；
[0015] 6.时频分析，由时变参数估计值求信号时频分布特征。
[0016] 其中，在所述步骤3中，由改进粒子群优化算法选择径向基函数的最优尺度。
[0017] 所述步骤3中，改进粒子群优化算法中速度更新公式的惯性权重为与适应度值有 关的变化值。
[0018] 所述步骤3中，衡量模型拟合效果的相关指数作为粒子的适应度值。
[0019] 本发明所提供的基于多尺度径向基函数和粒子群优化算法的非平稳信号时频分 析方法的优点包括：
[0020] 1.时频分辨率高，可准确辨识强噪声、非平稳、非线性的信号，准确提取其时频特 征；
[0021] 2.能根据信号灵活选择径向基函数尺度，适应性好；
[0022] 3.准确提取癫痫脑电信号与正常脑电信号时频特征，为癫痫脑疾病的辅助诊断提 供定量的技术指导。
【附图说明】
[0023] 图1为根据本发明实施例的时频分析方法的流程示意图。
[0024] 图2为粒子群优化方法选择多尺度径向基函数最优尺度的流程示意图。
[0025] 图3(a)显示了真实癫痫脑电信号。图3(b)_3(c)显示了根据本发明实施例的时 频分析方法与经典时频分析方法对真实癫痫脑电信号时频分析结果对比；其中图3(b)为 短时傅里叶方法的时频分析结果，图3(c)为本发明方法的时频分析结果。
[0026] 图4(a)显示了真实正常人脑电信号。图4(b)_4(c)显示了根据本发明实施例的 时频分析方法与经典时频分析方法对真实正常人脑电信号的时频分析结果对比；其中图 4(b)为短时傅里叶方法的时频分析结果，图4(c)为本发明方法的时频分析结果。
【具体实施方式】
[0027] 下面结合附图和【具体实施方式】对本发明作进一步的详细说明。
[0028] 真实脑电信号的时频分析作为本发明的实施例验证本发明的时频分析效果，并将 本发明与经典时频分析方法进行对比。真实脑电信号的时频分析包括对癫痫脑电信号和正 常人脑电信号的时频分析比较。
[0029] 实验数据来自德国波恩癫痫研宄室临床采集的脑电数据库，采样频率为173. 6Hz， 持续时间为23. 6秒。如图3(a)和图4(a)所示，分别为癫痫脑电信号和正常人脑电信号。
[0030] 图1展示了根据本发明实施例方法的流程图，包括：
[0031] 首先对信号进行时变参数建模，由最终预测误差（Finalpredictionerror,FPE) 模型定阶准则对模型进行定阶（步骤1);接着将时变参数用MRBF展开，表示为时不变参数 形式（步骤2);然后选择MRBF的中心并由PS0算法选择MRBF的最优尺度（步骤3);再对 时不变参数估计，将模型表示为矩阵形式，对时不变参数进行辨识（步骤4)，得到时变参数 估计值；由时不变参数估计值和MRBF估计时变参数（步骤5);最后进行时频分析，由时变 参数计算时频谱值，得到时频分布图（步骤6)。
[0032] 下面具体介绍根据本发明所提供的基于多尺度径向基函数和改进粒子群优化算 法的时频分析方法的具体方法：
[0033] 1.时变参数建模
[0034]时变参数模型参数是随时间变化的，一个p阶时变参数模型如下所示：
[0036] 式中：p为模型阶次，ai (n)(i= 1，2,…，p)为时变参数，e(n)是均值为0,方差为 0 2的尚斯白噪声。
[0037] 模型定阶：本发明采用最终预测误差（FPE)定阶准则选择模型阶数，其表达式如 下：
[0039] 式中，N为数据长度，p为模型阶数，&为预测误差。
[0040] 2.时变参数展开，将时变参数表示为一组多尺度径向基函数的加权线性组合，如 下式所示：
[0042]式中，ci;ni为基函数的时不变权系数，fjt)为高斯径向基函数。M为基函数的维 数。将上式代入式（1)得：
[0044] 式（4)即为时不变参数模型。
[0045] 3.径向基函数中心和尺度的确定
[0046] 本发明采用高斯径向基函数作为基函数扩展时变参数。由于高斯径向基函数具有 中心径向对称、逼近性能强、效果好，学习速度快等优点，且受信号特征的影响较小，所以本 发明选择高斯径向基函数作为基函数对时变参数进行展开，其表达式如下：
[0048]式中c为径向基函数的中心，〇 2为基函数的尺度，它决定了基函数围绕中心点的 距离。X为输入样本。在径向基函数作为基函数时，x(x= 1,2,…，乂?^为采样长度）为离 散采样时间序列。
[0049] 选择合适的径向基函数中心和尺度是准确辨识时变参数的关键。为了使径向基函 数分布到整个时变参数中，以确保径向基函数对时变参数的所有局部进行准确估计，本发 明将径向基函数的中心均匀分布到时变参数中，设第k个径向基函数及径向基函数的中心 分别表不为：
[00 52] 式中，M为径向基函数的维数，x为样本输入，ck为第k个径向基函数的中心，(J2k 为第k个径向基函数的尺度，它决定了该基函数围绕中心点的距离，N为采样数据长度。
[0053] 在本发明中，由改进粒子群优化算法计算径向基函数的最优尺度。
[0054] 图2为改进粒子群优化算法计算径向基函数的最优尺度流程图，具体步骤如下：
[0055] (a)初始化：初始化粒子、粒子个数、适应度值、局部最优粒子、全局最优粒子、粒 子最大值和最小值、速度最大值、迭代次数等。
[0056] 径向基函数的候选尺度为：
[0058] 式中，N为采样数据长度，M为径向基函数的维数，sk为需要进行调节的任意整数， 且最大值在10以内，粒子为Sk构成的向量为ui=[ss2，…，sM]。
[0059] 取随机整数为粒子赋初值，初始适应度值为0,粒子数和迭代次数可根据参数变化 特征来调整。
[0060] (b)计算每个粒子的适应度值
[0061] 粒子的适应度值是帮助找到最优粒子的启发性信息，在本发明中，采用模型拟合 效果的相关指数作为适应度函数，相关指数公式如下：
[0063] 式中，R2为相关指数，N为采样数据长度，；f为数据Y的模型预测输出。相关指数 越大说明模型拟合效果越好。
[0064] 将粒子代入尺度公式（8)，得到径向基函数的尺度，再根据径向基函数中心，得到 多尺度径向基函数。再由多尺度径向基函数展开方法对时变系统的模型参数辨识，得到测 量数据Y的模型预测值最后根据式（9)计算相关指数，得到粒子的适应度值。
[0065] (c)更新局部最优粒子。局部最优粒子为当前循环中适应值最大的粒子，找到适应 度值最大的粒子，该粒子作为局部最优粒子。
[0066] (d)更新全局最优粒子。全局最优粒子是所有循环中适应度值最大的粒子。若局 部最优粒子适应度值比全局最优粒子大，则赋值给全局最优粒子。否则全局最优粒子不变。
[0067] (e)更新每个粒子的速度，按粒子速度更新公式更新每个粒子的速度。
[0068] 粒子速度更新公式如下：
[0070]其中，Cbest为局部最优粒子，gbest为全局最优粒子，#>为第h次循环的第i个粒 子，V广为第h次循环的第i个粒子的速度，Ct为第h次循环的局部最优粒子，gt为第 h次循环后当前全局最优粒子。0为0~1之间的随机数，(^和(32的作用是防止粒子陷入 局部最优：
[0073] 式中，H为最大循环次数，h= 1，2, ???，!!?
[0074] 粒子速度更新公式中，《为惯性权重。在传统粒子群优化算法中，《 -般为常数， 不利于灵活寻找最优粒子。本发明中，采用一种改进粒子群优化算法，令惯性权重值与适应 度值相联系，以便更灵活寻找最优粒子，其表达式如下：
[0076] 式中，0是0~1之间的随机数，Fi是第i个粒子的适应度值，Fbest是全局最优粒 子的适应度值。
[0077] 当速度接近0时，按下式重新为粒子速度赋值：
[0079] 式中，
中的第j个元素，y是常值〇. 1，9是〇~1之间的随机数。
[0080] (f)更新每个粒子的位置。粒子位置更新公式如下：
[0082] (g)返回步骤（b)，重复执行步骤（b)、（c)、（d)、（e)和（f)，直到得到最优粒子。 最后得到全局最优粒子即为最优粒子，进一步根据尺度公式得到径向基函数的最优尺度。
[0083] 4.时不变参数估计。用最小二乘法对时不变参数模型中的时不变参数进行估 计，首先将时不变模型写为矩阵形式。
[0084]将式⑷中参数做如下定义：
[0085]f(t)= f2(t),???, fM(t)],
[0086] Hi(t) = y(t-i)f(t),
[0087]H(t)= H2(t),???, Hp(t)],
[0088] (；=[ci;1,ci;2，…，ci;M]，
[0089] C= [Q,^,-,Cp],
[0090] 则将式（4)写成矩阵形式：
[0091] y(t)=H(t)CT+e(t) (16)
[0092] 式中，T为矩阵的转置。C由最小二乘算法求得，即
[0094]式中，H= [H(1)，H(2)，…，H(N)，]Y= [y(l)，y(2)，...，y(N)]。
[0095] 5.时变参数估计。得到时不变参数估计值后，将时不变参数和多尺度径向基函数 代入式（3)即可得到时变参数的估计值。
[0096]6.时频分析。根据估计的时变参数，应用功率谱公式对信号进行时频分析，其计算 公式如下：
[0098] 式中，fs为采样频率，4⑴为时变参数的估计值，< 为观测误差的方差。
[0099] 求出信号的时频分布图，即得到信号的时频分布特征。
[0100] 本发明方法即基函数扩展法与经典时频分析方法即短时傅里叶变换方法进行比 较，对真实癫痫脑电信号及正常脑电信号时频结果进行分析，揭示脑电信号潜在的瞬时时 变特征。时频分析结果分别如图3(b)-3(c)和图4(b)-4(c)所示。
[0101] 由癫痫脑电信号和正常脑电信号时频分析结果可知，两种方法得到了较为一致 的频率特征分布，但显然短时傅里叶变换方法的时频分辨率低，结果如图3(b)和图4(b)所 示，时频图中频率分布不清晰。而本发明方法的时频分辨率高，结果如图3(c)和图4(c)所 不〇
[0102] 由本发明方法的时频图可清晰得到癫痫脑电信号和正常脑电信号的时频分布特 征。特别地，在癫痫脑电信号的时频分布图中，能够清晰地观察两种频段波：0波和0波。 如图3 (c)所示，0波（4-8Hz)分别在0~6秒之间、9~11秒及14~22秒之间出现；0波 (13-30HZ)出现在1~4秒、及12~22秒之间的时间段。类似地，在正常脑电信号的时频分 布图中，可以观察到两种不同频段的脑电波：S波和a波。如图4(c)所示，S波（〇-4Hz) 在4~5秒、10~12秒、14~16秒及18~22秒之间的时间段出现，a波（8-12Hz)主要 出现在17~20秒及14~15秒之间的时间段。
[0103] 由真实脑电信号的时频分析结果可知，与经典短时傅里叶变换时频方法相比，本 发明方法能够准确得到真实脑电信号的时频分布特征。由癫痫脑电信号和正常脑电信号时 频分析结果对比可知，正常脑电信号在整个时间段上，发生频率较低的脑电活动，而癫痫脑 电信号的脑电频率活动频率明显增大。而且，癫痫患者发作事件均被检测到，被检测到的发 作事件，尽管发作持续时间很短（2s)，本发明方法仍能检测到，并且所提取的脑电特征在 癫痫发作期具有明显的变化。本发明方法准确检测到癫痫发作的脑电信号中存在与正常脑 电信号不同频率的脑电波，而且提供了不同频率的瞬时频率分布信息，这些信息可以为医 生对癫痫脑疾病发作的辅助诊断提供定量的分析工具。
[0104] 本发明所提供的基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频 分析方法，主要是为准确提取脑电信号瞬时时频分布特征，辅助癫痫脑疾病的诊断提出的。 但显然，本说明书中所描述的时频分析方法也适用于其他非平稳信号的时频分析。
[0105] 以上对本发明所提供的基于多尺度径向基函数和改进粒子群优化算法的脑电信 号时频分析方法进行了详细的说明，但显然本发明的范围并不局限于此。在不脱离所附权 利要求书所限定的保护范围的情况下，对上述实施例的各种改变都在本发明的范围之内。
【主权项】
1. 一种基于多尺度径向基函数与改进粒子群优化算法的脑电信号时频分析方法，其特 征在于包括： 步骤1.时变参数建模； 步骤2.时变参数展开，时变模型转化为时不变参数模型； 步骤3.选择径向基函数的中心及尺度； 步骤4.时不变参数估计； 步骤5.时变参数估计； 步骤6.时频分析，由时变参数估计值求信号时频分布特征。2. 如权利要求1所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，由改进粒子群优化算法选择径向基函数的最优尺度。3. 如权利要求2所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，改进粒子群优化算法中速度更新公式的惯性权重为与适应度值有 关的变化值。4. 如权利要求2所述的基于多尺度径向基函数与改进粒子群优化算法的脑电信号时 频分析方法，其特征在于： 所述步骤3中包括，衡量模型拟合效果的相关指数作为粒子的适应度值。
【专利摘要】本发明提出了一种基于多尺度径向基函数和改进粒子群优化算法的脑电信号时频分析方法。该方法应用时变参数建模法对脑电信号进行时频特征提取分析，引入多尺度径向基函数展开式辨识时变参数。首先建立时变参数模型，然后用基函数扩展法辨识模型时变参数，即将时变参数表示为一组多尺度径向基函数的线性加权组合，将时变参数的辨识问题转化为时不变参数的辨识。径向基函数的最优尺度由粒子群优化算法决定。最后根据时变参数估计值和功率谱密度公式计算脑电信号的时频分布特征。与现有时频分析方法相比，本发明方法能同时获得较高的时间和频率分辨率，能够准确提取脑电信号的时频分布特征，对癫痫脑电信号的应用及辅助癫痫疾病的诊断具有重要意义。
【IPC分类】G06F19/00, G06K9/00
【公开号】CN104899436
【申请号】CN201510284547
【发明人】李阳, 刘青, 王旭东
【申请人】北京航空航天大学
【公开日】2015年9月9日
【申请日】2015年5月29日
最新回复(0)