基于小波包变换的生物电脑部信号分析方法

摘要

本案涉及一种基于小波包变换的生物电脑部信号分析方法，包括：对采集到的脑部信号进行预处理，得到预处理信号；针对该预处理信号，选择匹配度最高的小波基函数进行小波变换分解，获得感兴趣频带或频率子带的特征信号；通过建立阈值T来对所述特征信号进行状态划分，最终得到感兴趣频带或频率子带的状态编码。本案利用小波包变换来更加精确的对感兴趣的频带或频率子带进行提取，并进行状态判别，联合不同的频带状态对患者状态进行刻画，从而指导对患者状态的分析；本方法具有状态特异性和自适应性。

说明书

技术领域

本发明涉及一种生物电脑部信号的数据处理方法，特别涉及一种基于小波包变换的脑部场电位信号分析方法。

背景技术

患者状态的精确确定对于临床分析十分重要，但是目前很多生理状态都没有客观评价的标准如运动障碍的程度，慢性疼痛的等级，通常对这些状态的确定主要依靠医生经验以及患者自身的体验，具有很大的主观性。对于人体生物电脑信号的研究发现大脑活动的水平跟很多生理状态相关，这些信号可以应用于患者状态的定量评估，并指导患者的治疗。

生物电脑部信号中蕴含着大量的信息，但只有部分频率带的信息与特定状态相关，因此需要对频率带信息进行提取，并得到这些信息与患者状态之间的关系。但是在现有技术中，对这些信息的提取方法比较单一，基本类似“一刀切”的办法，即设定一个功率阈值，超过该阈值是一种状态，未超过的是另一种状态，这种粗略的提取方法显然为后续的状态分析提供了较高的误判信息。

发明内容

针对现有技术存在的不足，本案提供了一种基于小波包变换的生物电脑部信号分析方法。小波包变换提供了一种更加精细的数据分类方法，它能够动态有效的提取窄带信号，并且可以根据需求自适应地选择频带然后根据频率特征选择合适的小波函数得到准确的频带信息。当多个频带与患者状态相关时进行联合表达能够更加精细的对患者状态进行刻画。

为实现上述目的，本案通过以下技术方案实现：

一种基于小波包变换的生物电脑部信号分析方法，其包括：

步骤1)预处理：

对采集到的脑部信号的原始数据进行筛选，挑选出有效信号；

对所述有效信号进行计算得到差分信号；

对所述差分信号进行去工频干扰和基线漂移，最终得到预处理信号；

步骤2)特征提取：

针对所述预处理信号，选择匹配度最高的小波基函数进行小波包变换分解，获得感兴趣频带或频率子带的特征信号；

步骤3)状态判别：

通过建立阈值T来对所述特征信号进行状态划分，最终得到感兴趣频带或频率子带的状态编码。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，在步骤1)中还包括对所述差分信号进行低通滤波。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，在步骤1)中还包括对所述差分信号进行降采样处理。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，步骤2)中，所述小波基函数的选择方法为：

根据感兴趣频带或频率子带的特征信号在对称性、正交性、紧支性方面的要求，粗选出若干个具有相似程度的小波基函数，随后分别代入小波包变换函数中分别计算得到若干个小波包系数的熵，将最小的熵所对应的小波包基函数确立为匹配度最高的小波基函数。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，感兴趣频带或频率子带的特征信号是利用小波包变换后得到的对应的感兴趣频带或频率子带的小波包系数C来进行表征。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，步骤3)中，所述阈值T通过以下方法获得：

对早于当前数据的一段时间内的小波包系数C取中值，得median(|C|)；

设σ＝median(|C|)/0.6745；

将σ代入minimaxi函数得到计算结果t；

将t乘以阈值权重，得到当前数据所对应的阈值T。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，在步骤3)中，状态划分的方法为：

以阈值T为临界，将每个小波包系数C分别与阈值T进行比较，从而将所有C≥T的情形和所有C＜T的情形区分开。

优选的是，所述的基于小波包变换的生物电脑部信号分析方法，其中，所述状态编码采用二进制编码。

本发明的有益效果是：本案利用小波包变换来更加精确的对感兴趣的频带或频率子带进行提取，并进行状态判别，联合不同的频带状态对患者状态进行刻画，从而指导对患者状态的分析；本方法具有状态特异性和自适应性，选用小波包变换更为精细地提取与患者状态相关的生物电脑信号特征，并且提供了解决小波包变换中参数的选择方案，有效的提高的小波包变换的有效性和准确性；并且本案选用阈值模型对患者状态进行判别，对状态进行编码，能够更加精细的对疾病状态进行划分。

附图说明

图1为基于生物电脑信号实现患者状态判别的模型框架图。

图2为患者状态判别算法模型中的预处理示意图。

图3为患者状态判别算法模型中的特征提取示意图。

图4为对提取的特征进行状态判别的示意图。

图5为利用小波包变换提取出的theta和alpha频带活动示意图。

图6为不同小波基函数在提取theta和alpha频带活动的小波包系数熵。

图7为生成的alpha频带模拟信号，其中心频率为10.5Hz，信噪比为-5dB，包络是周期为4s的梯形波图。

图8为利用算法模型对模拟信号进行实验得到的结果与理论结果的对比图。

图9为利用本案模型对真实慢性疼痛信号进行信号分析提取得到theta和alpha频带的小波包系数并对各节律活动状态进行划分。

图10为利用theta和alpha节律活动的状态对患者状态进行编码。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本案列出一实施例的基于小波包变换的生物电脑部信号分析方法，其包括：

步骤1)预处理：

对采集到的脑部信号的原始数据进行筛选，挑选出有效信号；

对有效信号进行计算得到差分信号；

对差分信号进行去工频干扰和基线漂移，最终得到预处理信号；

步骤2)特征提取：

针对预处理信号，选择匹配度最高的小波基函数进行小波变换分解，获得感兴趣频带或频率子带的特征信号；

步骤3)状态判别：

通过建立阈值T来对特征信号进行状态划分，最终得到感兴趣频带或频率子带的状态编码。

其中，在步骤1)中还包括对差分信号进行低通滤波。

其中，在步骤1)中还包括对差分信号进行降采样处理。

其中，步骤2)中，小波基函数的选择方法为：

根据感兴趣频带或频率子带的特征信号在对称性、正交性、紧支性方面的要求，粗选出若干个具有相似程度的小波基函数，随后分别代入小波包变换函数中分别计算得到若干个小波包系数的熵，将最小的熵所对应的小波基函数确立为匹配度最高的小波包基函数。

其中，感兴趣频带或频率子带的特征信号是利用小波包变换后得到的对应的感兴趣频带或频率子带的小波包系数C来进行表征。

其中，步骤3)中，阈值T通过以下方法获得：

对早于当前数据的一段时间内的小波包系数C取中值，得median(|C|)；

设σ＝median(|C|)/0.6745；

将σ代入minimaxi函数得到计算结果t；

将t乘以阈值权重，得到当前数据所对应的阈值T。

其中，在步骤3)中，状态划分的方法为：

以阈值T为临界，将每个小波包系数C分别与阈值T进行比较，从而将所有C≥T的情形和所有C＜T的情形区分开。

其中，状态编码优选采用二进制编码。选择的感兴趣频带或频率子带的个数将决定编码的位数，例如只有1个感兴趣频带或频率子带时，分为两种状态：大于等于阈值的状态为“1”，而小于阈值的状态为“0”；若有两个核团分别有1个和2个感兴趣频率子带，那么患者状态将有8种组合模式：000、001、010、011、100、101、110和111。本案并不限于二进制编码，其他的编码形式亦可以。

根据实际情况确定感兴趣频带和频率子带，其中实际情况指：例如：帕金森状态与STN核团的beta频段(15-30Hz)密切相关，再如研究发现疼痛状态与PAG和VPL核团神经活动的theta和alpha频带或频率子带相关等。本案可以用于任何频带，示例频带包括δ带、α带、β带和γ带以及这些频带的频率子带。频率子带具有比频带窄的频率宽度，因此频带可由多个频率子带进行定义。

生物电脑部信号的不同频带与不同状态相关联。在以下表1中示出频带划分的一个示例：

表1、频带划分

频带范围(Hz)δ频带1.5-4Hzθ频带4-8Hzα频带8-13Hzβ频带13-30Hzγ频带30-90Hz

然而，频带范围不仅限于上述表格中的范围，可以包括多个频率子带(例如，可由多个频率子带组成)，每个频率子带具有比频带更窄的宽度，频带的频率子带可具有相同或不同的宽度。

所检测的生物电脑部信号的特定频带内的活动可能由于患者的状态发生改变而改变，也可能是患者状态发生改变的原因，因此特定频带内的生物电活动与患者状态密切相关，可用来指示患者状态。

相比一个单独的频带信号的状态，多个核团的多个频率子带的活动状态能够更好的对患者状态进行划分，例如多个频率子带可以划分疼痛的不同等级，从一级到十级。基于一个或多个核团的一个或多个感兴趣频带或频率子带的小波包系数计算阈值划分状态可以对患者状态进行分析。

本实施例以慢性疼痛状态为例，但在其他示例中，除了慢性疼痛之外，该算法模型可应用于其他患者状态的分析，诸如，但不限于，运动障碍等状态判断。

虽然贯穿本案的生物电脑部信号主要以人脑深部场电位信号(LFP)为例，但本案所提供的算法还可用于诸如，但不限于，头皮脑电(EEG)，脑磁图，功能磁共振数据等与患者状态分析相关的生物电脑部信号。

本算法模型主要分为3个模块，其结构示意图如图1所示：预处理模块、特征提取模块和状态判别模块。

算法模型的具体实现方式如下：

为了使结果更加精准，对采集到的原始数据进行预处理分析，如图2所示，主要包括以下步骤：

1)对信号质量不佳的信号进行剔除，包括利用信号判断接收到的信号是否来自于预先选择的大脑核团。

2)从记录电极的四个触点得到的四通道信号计算出差分信号；

3)对信号进行50Hz陷波滤波去除工频干扰，由于与患者状态相关的局部场电位主要集中在90Hz以下，因此可以做90Hz低通滤波去除高频信号，然后做3Hz高通滤波去除掉基线漂移，最后信号降采样到采样率为500Hz以减少计算量，提高计算速度。

对预处理后的数据进行特征提取，具体步骤如图3所示的特征提取模块。利用小波包变换分解数据，提取感兴趣频带或频率子带的信号，利用小波包系数表征这些频带信号活动的情况。

特征提取是实时进行的，特征提取模块对收到的信号进行很短的时间窗处理，例如每接收到0.2s信号处理一次，也可以选择其他较短的时间窗。本算法模型具有实时性，其实现方式是每个模块每次只处理很短时间的最新数据，且每个模块的处理窗长均相同，在本例中处理窗长经选择后确认为针对theta(6-9Hz)频带为0.33s，而alpha(9-12Hz)频带为0.29s。同时状态判别模块的先验数据也具有实时性，其方式是当判别模型完成一个时间窗内数据的判别后就将窗内数据纳入到先验数据中，并丢弃原始先验数据中的尾部等长的数据，因此既保持先验信息的长度，又同时进行更新。

在计算窗内利用小波包变换对局部场电位信号进行分解，并提取特征频段，得到对应的小波包系数。在做小波包变换前首先要进行参数选择，包括：确定合适的小波基函数来提取特定的信号节律，以及分解时分解层数，熵，和提取节点等参数的设定。但由于分解层数，熵，和提取节点这些参数的选择在当前状态下是本领域所公知通用的，因此，这些参数相当于就是已知的，唯一需要选择的就是小波基函数。

上述对于小波包变换中小波基函数的选择方法是：结合需要提取的与患者状态相关的节律特征的特性和小波基函数本身的特性，如对称性、正交性和紧支性对小波基函数进行粗筛选得到一组比较可靠的小波基函数组。比较利用这些小波基函数提取特征后的小波包系数的熵，熵越小对应的小波基函数越适合当前的特征提取，从而得到针对该节律特征的最佳小波基函数。图6为利用多组真实信号比较不同小波基函数的结果，统计发现bior3.7小波基在提取alpha频带熵相对于其他具有良好的表现，而对于theta频带rbio3.7较好。

上述做小波包分解提取特征时分解层数和提取节点的确定方法：需要根据特征节律的频段改变原始信号的采样率，使所需节律能够从小波包分解树上完整的提取出来。例如特征节律的频段为9-12Hz，则需要把原始信号重采样使采样率变为384Hz，这样根据小波包分解定律，9-12Hz存在于节点[6 2],因此分解层数为6层，提取节点为[6 2]。

上述小波包变换的参数中熵使用Shannon熵。

对提取的特征进行状态判别，方法如图4所示，状态判别模块基于提取特征的小波包系数利用阈值模型计算阈值，大于等于阈值对应状态“1”，小于阈值对应状态“0”。这种对状态的编码是采用二进制编码，本案也可以选择其他编码形式。

上述阈值模型是基于阈值去噪构建的，其理论为具有模式特征的信号节律的小波包系数要大于等于随机状态下的小波包系数。其具体步骤包括：

1)得到对噪声水平进行估计σ，为median(|C_j,k|)/0.6745,其中C_j,k为小波包系数。

2)选用”minimaxi”阈值估计方法估计阈值t；

3)阈值t在乘以阈值权重得到最终的阈值T，本例中阈值权重设为1。

进一步地，特征判别模块具有预测能力，使用固定长度的先验数据计算阈值作为对当前数据阈值的预测。

为了验证本案方法的有效性实验生成了9-12Hz节律的模拟信号(图8)，其中心频率为10.5Hz的信号并加入-5dB噪声信号，其包络是周期为4s的梯形波，利用算法模型计算得出小波包变换的小波基函数选择bior3.7，分解层数为6层，提取节点为[6 2],并理论上计算开启和关闭刺激的时间。图8中还包括有利用本案的算法模型计算得到的动态阈值和状态划分，并与理论状态进行比较，结果显示算法有着较高的精度和准确度，算法的正确开启率为95.14％，正确关闭率为98.28％。需要说明的是模拟信号的信噪比对于准确率有较大的影响。

在本例中，利用该算法对慢性疼痛患者的场电位信号进行状态划分，感兴趣频带分别为theta(6-9Hz)和alpha频带(9-12Hz)(图9)，其中信号主要来自PVAG和感觉丘脑核团。实验发现数据中同时存在由两种频率段组合的四种模式(图10)：“00”，“01”，“10”，“11”，并且这四种模式的百分比与患者的状态有关。

模拟数据与真实疼痛场电位数据的实验证明了算法模型能够有效实现患者状态的判别，克服了以往单纯的使用功率谱分析作为与患者状态相关的度量，同时也解决了准确提取节律的问题。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。