一种基于Stacking的帕金森震颤检测方法

摘要

本发明属于信号处理以及模式识别领域，具体涉及一种基于Stacking的帕金森震颤检测方法，该方法包括获取原始信号，对获取的原始信号进行预处理；提取预处理后信号的特征，提取的特征包括信号的时域特征、时频特征以及非线性特征；对提取的特征进行筛选，得到特征子集；将特征子集输入到训练好的Stacking模型中进行帕金森震颤检测；根据检测结果判断待检测人的状态；本发明中通过计算多种特征，并对特征进行筛选，结合提出的新的、泛化能力强的、准确率高的帕金森震颤检测方法，能够准确检测帕金森震颤，为帕金森病诊断提供参考。

说明书

技术领域

本发明属于信号处理以及模式识别领域，具体涉及一种基于Stacking的帕金森震颤检测方法。

背景技术

帕金森病是一种常见的中老年疾病，是一种具有进行性的神经退行性疾病。帕金森病的致死率是比较低的，但是它的致残率很高，对于老年人的身体健康有着严重的威胁。目前，对于帕金森病的致病原因尚在猜测之中，没有定论。有研究表明可能的原因有遗传因素、环境因素、人口老龄化等。帕金森病的临床症状主要包括两个方面：运动相关和非运动相关症状。运动相关的症状主要有静息震颤、肌肉强直、步态异常和运动迟缓等。非运动相关的症状包括抑郁、睡眠障碍和便秘等。据统计有超过70％的PD患者以静息性震颤为首发症状。

目前，采用帕金森病评定量表(Unified Parkinson Disease Rating Scale，UPDRS)对帕金森病进行诊断。UPDRS是一种分级量表，根据严重程度分为六级：无症状、单侧发病、双侧疾病没有平衡障碍、中度双侧疾病、严重残疾、只能坐轮椅或者卧床。在UPDRS中静息震颤也是一个重要的指标。研究一种帕金森震颤检测方法，提供一种客观的震颤检测检测指标。将其作为医生诊断方法的一个重要辅助手段。

利用不同的特征采集震颤信号，提取包括时域、频域以及其他的重要的特征进行分析。主要研究信号源包括局部场电位信号(Local field potential，LFP)、肢体的表面肌电信号(Electromyography，EMG)、脑部电信号(Electrocardiography，EEG)、肢体运动参数的加速度信号(Acceleration，ACC)。在进行震颤信号采集过程中George Rigas通过四个三轴加速度传感器，记录四个不同位置的加速度信号，采用通过离散傅里叶变换，提取与静息震颤相关的频域特征，根据这些特征采用人工神经网络和SVM对静息震颤进行分类，分类准确率为76％。Ada Zhang针对加速度传感器信号，利用多层感知机和随机森林对帕金森震颤进检测最高准确率为88.7％。Pedrosa等人利用利用光电传感器获取的加速度传感器信号对震颤进行评级准确率为92.8％。对于加速度传感器信号学者们利用通过谱估计、奇异值分解等方法对震颤进行分类。也有通过计算信号的时域、频域特征以及统计特征用于帕金森震颤分级评价的。使用人工神经网络、支持向量机、随机森林等对帕金森震颤进行检测。但是上述方法中，存在震颤信号分类的准确率低，检测效果差的缺点。

发明内容

为解决以上现有技术存在的问题，本发明提出了一种基于Stacking的帕金森震颤检测方法，该方法包括：获取原始信号，对获取的原始信号进行预处理；提取预处理后信号的特征，提取的特征包括信号的时域特征、时频特征以及非线性特征；对提取的特征进行筛选，得到特征子集；将特征子集输入到训练好的Stacking模型中进行帕金森震颤检测；根据检测结果判断待检测人的状态。

优选的，对原始信号进行预处理包括：采用小波变换的方法对原始信号进行降噪处理；对原始信号进行小波分解的层数为3层，小波函数为‘db4’，设置的阈值为启发式阈值和软阈值。

优选的，所述时域特征包括：hjorth参数、平均振幅、最大振幅以及均方根值；所述频域特征包括频谱均值、频谱方差、频谱的一阶导数、信号总功率以及频谱峰值功率；所述非线性特征包括样本熵。

进一步的，提取时域特征的具体过程包括：根据预处理后的特征计算信号的hjorth参数，所述hjorth参数包括hjorth的活动性、移动性以及复杂性；hjorth参数的活动性为信号的方差；hjorth参数的移动性为信号一阶导数方差与信号方差的均方根值；hjorth参数的复杂性为信号的一阶导数的移动性与信号的移动性的比值；计算平均振幅和最大振幅的过程包括采用envelope函数计算信号的上包络线env

进一步的，提取频域特征的具体过程包括：采用fft函数计算信号的频谱；通过mean函数对得到的频谱信号进行计算，得到频谱均值；使用var函数对频谱信号进行计算，得到频谱方差；频谱一阶导数为后一点频谱信号减去前一点频谱信号之差的和除以频谱点数后减1；信号总功率为频谱信号中每个点的平方之和除以频谱点数；计算频谱信号的最大值，并获取最大值对应频点位置，选取对应频点位置以及该位置的左边2个点和右边2个点，计算选取的5个点的平均峰值功率，将该峰值功率作为频谱的峰值功率。

进一步的，计算非线性特征的具体过程包括：根据预处理后的信号构建m维空间向量；计算空间向量中两个点的之间的距离；设置偏差阈值r，根据偏差阈值对对计算得到空间向量中的各个点的距离进行筛选，计算筛选出的点与总距离的比值，并根据该比值对所有筛选出的点作平均计算，得到m维空间向量的平均值；对空间向量的维度加1，并重复上述过程，得到m+1维空间向量的平均值和；根据m维空间向量的平均值和m+1维空间向量的平均值使用SampEn函数计算信号的样本熵。

优选的，对提取的特征进行筛选的过程包括：

步骤1：计算每个特征在分类算法中的准确率，将各个特征按准确率由高到低排序；

步骤2：设置特征子集，将设置的特征子集进行初始化，初始化后的特征子集为空集；

步骤3：将排序后的特征按顺序依次加入到特征子集中；

步骤4：根据加入到特征子集中的各个特征的算法准确率判断算法性能是否提升，若算法性能提升，则保留特征子集中的特征，若算法性能不提升，则删除特征子集中的特征；

步骤5：设置阈值，获取特征子集中择特的征数量，判断特征数量与设置的阈值大小，若大于设置的阈值，则将特征子集中的特征作为筛选结果，否则重复步骤3-步骤5，直到所有的特征筛选完为止。

优选的，训练Stacking模型的过程包括：

步骤1：获取原始数据，并提取原始数据的特征，得到m个特征；根据m个特征获取样本总量，样本总量为N*m；其中，N表示信号的长度，m表示特征数量；

步骤2：将样本总量按照7:2:1的比列进行划分，得到对应的N

步骤3：采用N

步骤4：将N

步骤5：将预测结果作为训练集，训练一个LR模型；

步骤6：训练得到的KNN模型、SVM模型、RF模型和Adaboost模型作为第一层，LR模型作为第二层，得到训练好的stacking模型；

步骤7：将N

本发明的优点：

本发明通过对原始信号进行小波变换处理，降低了原始信号的噪声；本发明通过提取信号的多个特征，并对提取的特征进行筛选，将筛选后的特征输入到Stacking模型中进行帕金森震颤检测，提高了检测的效率和准确率。

附图说明

图1为本发明中帕金森震颤检测方法流程图；

图2为本发明中原始加速度信号与小波降噪后信号对比图；

图3为本发明的前向选择流程图；

图4为本发明的基于Stacking的帕金森震颤检测模型图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种基于Stacking的帕金森震颤检测方法，如图1所示，该方法包括：获取原始信号，对获取的原始信号进行预处理；提取预处理后信号的特征，提取的特征包括信号的时域特征、时频特征以及非线性特征；对提取的特征进行筛选，得到特征子集；将特征子集输入到训练好的Stacking模型中进行帕金森震颤检测；根据检测结果判断待检测人的状态。

优选的，获取的原始信号包括采用加速度计采集患者的震颤信号，根据采集的震颤信号得到锁相DBS数据。数据由信号和触发时间序列组成，分别对应于主震颤轴加速度计数据和刺激触发器的数据。触发器是发送到刺激装置的信号。在记录数据过程中，关闭标准高频DBS，将一个三轴加速度计固定在患者最颤抖的手食指掌指关节上。加速度计信号被生物计量学K800放大器放大，来自主要震颤轴的信号被1902隔离前置放大器在线滤波，然后被记录和处理使用1401放大器和Spike2软件进行记录。对于在线滤波器，使用了一个两极数字巴特沃思滤波器，其截止频率为患者震颤频率周围的2Hz。平均震颤频率为4.22±0.25Hz。在无刺激时，三轴加速度计在震颤频率处的频谱峰最高的信号被定义为主震颤轴。后续通过主成分分析表明所选的数据为主要震颤轴的数据。将数据分别标注为震颤与非震颤状态。

对获取的数据进行预处理包括：使用小波变换方法对原始信号进行降噪处理，可以有效地提高信号的信噪比，小波降噪的主要参数：设置小波分解层数为3层，选择小波函数为‘db4’，阈值选用启发式阈值和软阈值。原始加速度传感器信号及其去噪后的信号如图2所示。

降噪后的加速度传感器信号与原始信号相比，信噪比为8.91dB，均方根误差为0.0016。从图中可以看出，原始信号上局部混杂较为严重，而通过降噪处理之后，混杂部分得到有效消除，并且降噪后的信号能够保留原始信号的特点，减少失真，对后续的帕金森震颤检测有一定帮助。

提取预处理后信号的特征，提取的特征包括信号的时域特征、时频特征以及非线性特征。时域特征主要包括hjorth参数、平均振幅、最大振幅、以及均方根值。频域特征包括频谱的均值、频谱方差、频谱的一阶导数、以及信号的总能量、频谱的峰值功率。非线性特征计算信号的样本熵。

时域特征

计算时域特征的过程包括：信号的Hjorth参数描述了信号在时域中的统计特性。它的参数包括活动性(Activity)，移动性(Mobility)和复杂性(Complexity)。其中活动性能表明信号的方差，移动性能测量平均频率。另外在给定周期时间内，复杂性能表明频率的变化。y(t)为原始信号,var为计算数据方差，计算公式为：

Activity＝var(y(t))

其中，Activity表示活动性参数，var表示数据方差，y(t)表示原始信号，Complexity表示负杂性参数，

在计算加速度信号的振奋过程中，通过计算信号的平均幅度和最大振幅得到信号的振幅。为了计算加速度信号的平均振幅，提取了总加速度的包络线。信号在时间窗内的平均幅值定义为上包络线与下包络线之差的平均值定义如下：

其中，Mean_A表示平均振幅，N表示总采样点数，env

最大振幅的计算公式如下：

Max_A＝max(env

时频特征

计算时频特征的过程包括：计算信号在频谱上的特征，对于原始信号x(n)，通过快速傅里叶变换得到信号对应的频谱X(k)。并计算信号的总能量，频谱均值频谱方差、信号的峰值功率。

信号的总能量是指频谱系数平方之和。对于长度为n的X(k)，公式如下：

其中，e表示信号的总能量，k表示频谱对，n表示频谱的长度，X(k)表示通过快速傅里叶变换得到信号对应的频谱。

均值是指频谱系数的均值。计算对于长度为n的X(k)，公式如下：

其中，P_mean表示频谱系数的均值。

方差能够表示能量对于平均能量的离散水平。对于长度为n的X(k)，公式如下：

其中，σ表示方差，n表示频谱长度，μ表示频谱均值。

频谱一阶导数是指每一能级能量的第一阶导数的平均值。对于长度为n的X(k)，公式如下：

其中，δ表示频谱一阶导数。

频谱的峰值功率一般是指在一段时间内传感器信号单边功率谱在主频的±0.5Hz附近的功率谱估计，其计算公式为：

其中，P

计算非线性特征

样本熵(SampEn)是通过计算信号中产生新模式的大小来衡量时间序列的复杂性。样本熵的值越小，则说明信号的自似性越大。样本熵越大，则说明信号的越没有规律性。

计算非线性特征的具体过程包括：根据预处理后的信号构建m维空间向量；计算空间向量中两个点的之间的距离；设置偏差阈值r，根据偏差阈值对对计算得到空间向量中的各个点的距离进行筛选，计算筛选出的点与总距离的比值，并根据该比值对所有筛选出的点作平均计算，得到m维空间向量的平均值；对空间向量的维度加1，并重复上述过程，得到m+1维空间向量的平均值和；根据m维空间向量的平均值和m+1维空间向量的平均值使用SampEn函数计算信号的样本熵。对应的参数m取值为2，r取值为信号标准差的0.15倍。

对于一个原始信号长度为N的时间序列表示为{u(i):1

步骤1：构造一个m维空间向量X(1)，X(2)，...X(N-m+1)；即

X(i)＝{u(i),u(i+1),...,u(i+m)}

其中，u(i)表示时间序列第i个点的值，m表示向量空间的维数。

步骤2：定义向量X(i)和X(j)之间的距离，d[X(i),X(j)]为两个向量对应元素差值最大值。

步骤3：对于每一个{i:1

其中，

对所有的1

其中，

步骤4：将m的维度增加1，重复步骤3的过程得到

步骤5：样本熵的计算公式为：

SampEn(N,m,r)＝-ln[φ

其中，SampEn(N,m,r)表示样本熵，N表示信号的长度，m表示空间向量的维数，r表示偏差，

根据上面的计算过程可知样本熵的取值与N，m，r息息相关。其中m取1，r取标准差的0.15倍。

对提取的特征进行筛选的过程包括：特征选择模式识别和机器学习的一个重要过程，通过基于过滤式、包裹式或者嵌入式的特征筛选方法，从计算得到的特征中，选择出适合用于分类的特征，降低数据集的维度的同时，提高算法的性能。筛选特征子集时，首先以单个特征的准确率作为特征重要性评价指标将特征根据重要性排序。

如图3所示，首先设置特征子集为空集，然后根据特征的重要性，从未选择的特征集中选择一个特征加入该特征子集。然后，比较加入该特征之后的特征子集在后续分类算法中的表现，如果相较于该特征加入之前效果有所提升，则保留该特征。反之，则舍弃该特征，直到选择到合适的子集使得分类算法达到最优。

构建一个两层的Stacking模型，用于帕金森震颤检测。根据Stacking特的点，选择不同种类的分类器或者使用同一学习器，而采用不同参数构成，另外在为了防止过拟合，第二层需要采用较为简单的模型。因此在构造该Stacking模型的时候，本文采用四个不同的分类器作为第一层，分别为KNN、SVM、RF以及Adaboost。而采用相对较为经典且容易解释的LR作为第二层。基于Stacking的帕金森震颤检测的模型结构如下图4所示。

训练Stacking模型的过程包括：

步骤1：获取原始数据，并提取原始数据的特征，得到m个特征；根据m个特征获取样本总量，样本总量为N*m；其中，N表示信号的长度，m表示特征数量；

步骤2：将样本总量按照7:2:1的比列进行划分，得到对应的N

步骤3：采用N

步骤4：将N

步骤5：将预测结果作为训练集，训练一个LR模型；

步骤6：训练得到的KNN模型、SVM模型、RF模型和Adaboost模型作为第一层，LR模型作为第二层，得到训练好的stacking模型；

步骤7：将N

采用N

步骤1：采用fitcknn函数对N

步骤2：采用fitcsvm函数对N

步骤3：采用随机森林TreeBagger函数对N

步骤4：采用fitensemble函数对N

训练LR模型模型的过程包括：以N

采用帕金森震颤检测的模型对信号进行处理，得到12种特征，即Hjorth活动性参数(1)，hjorth复杂性(2)，hjorth移动性(3)，信号的峰值(4)，平均振幅(5)，最大振幅(6)，均方根值(6)，频谱的均值(7)，频谱方差(8)，频谱的一阶导数(9)，信号的总能量(10)，频谱的峰值功率(11)，样本熵(12)。

对于LR、SVM、KNN、RF、Adaboost以及提出的基于Stacking的帕金森震颤检测六种分类算法，利用特征筛选方法，采用默认参数设置，得到对应分类方法的分类的最高准确率以及其对应的特征组合如下表1.1所示。

表1.1各分类算法的最高准确率及对应特征组合

在震颤状态与非震颤状态分类中，共使用4386个样本进行五倍交叉验证。分类准确率体现了分类算法的识别率高低，是一个评估分类算法识别能力大小的重要指标，但是为了计算分类算法的泛化能力、验证分类准确率的真实性计算分类算法的精准性，精准度、召回率、以及F1score进行计算如表1.2所示。

表1.2各分类算法的评价指标

震颤检测的目的是区分患者是处于静息震颤状态还是非震颤状态。震颤的准确率表示分类后的真实震颤样本数量和标注的震颤样本数量相比较。当精准度的时候，模型较为容易识别震颤状态，这也有利于实际应用中的震颤检测。结果表明在准确率，召回率以及F1 score三个参数上，基于Stacking的震颤检测算法都优于其他的震颤检测分类算法。说明该算法得到的准确率可靠，具有很高的一致性以及泛化能力。

以上所举实施例，对本发明的目的、技术方案和优点进行了进一步的详细说明，所应理解的是，以上所举实施例仅为本发明的优选实施方式而已，并不用以限制本发明，凡在本发明的精神和原则之内对本发明所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。