一种检测帕金森病症状严重程度的脑电数据分析方法

技术领域

本发明涉及一种脑电数据分析方法，尤其涉及一种检测帕金森病症状严重程度的脑电数据分析方法。

背景技术

深部脑刺激(Deep Brain Stimulation，DBS)是一种治疗帕金森疾病的有效治疗方式。传统的刺激方式为持续性周期刺激，刺激参数固定，但是帕金森病患者的症状会随时间变化而变化，取决于患者认知和运动状态，服用多巴胺能药物后症状也会有所改善，因此持续性周期刺激可能会产生副作用，如语言障碍。

根据帕金森病患者症状实时调整脑刺激参数的系统称之为闭环神经刺激器(Closed-loop Nerve Stimulator，CNS)或者自适应深部脑刺激器(adaptive DBS，aDBS)，只在需要时进行刺激，既能减少刺激带来的副作用，又可以延长电池使用寿命。

基底神经节局部场电位中beta(下文称之为病理性beta)活性增强是帕金森病的典型特征，研究表明，beta活性与运动迟缓或僵硬症状之间存在相关性：beta频带能量越高，症状越明显；beta频带能量越低，症状改善，呈正相关。所以目前已有闭环神经刺激方法多以此特征为反馈控制信号。

但是以beta活性作为反馈控制信号会存在一些问题。第一，不同帕金森病患者的病理性beta频带会有所不同，不一定能准确找到。第二，某些患者脑电信号中可能会出现多条beta频率子带，难以确定哪一个是病理性频带。第三，大脑beta频带与运动相关，运动或者运动想象会使大脑中beta活性增强，影响反馈结果。第四，帕金森病静息性震颤症状与beta活性无关。因此，以beta活性作为闭环神经刺激的反馈控制信号是不合适的。

发明内容

针对上述问题，本发明的目的是实现一种基于帕金森病患者症状的闭环神经刺激方法，并应用于CNS或者aDBS这种植入式脉冲刺激器(Implantable Pulse Generator，IPG)中，这就要求反馈算法满足如下条件：

第一，算法所用反馈参数可以准确描述出帕金森病患者症状严重程度。

第二，算法每次计算数据量不宜过大，算法本身不能过于复杂。数据量太大会占用过多内存，通常IPG内存都有限。而算法过于复杂会导致计算时间过长影响系统实时性，算法的复杂程度与计算时间成正比。

通过分析帕金森病患者丘脑底核(Subthalamic Nucleus，STN)脑电信号发现，帕金森病患者在服用多巴胺能药物后，药效生效时，患者症状改善，脑电信号中病理性beta能量降低，同时伴随低于10Hz的低频段能量上升。而当药物失去效果后，信号中病理性beta活性重新升高，同时伴随低频段能量下降。因此，除了病理性beta频带能量外，低频段能量变化也与帕金森病症状之间存在相关性：低频段能量越低，症状越明显；低频段能量越高，症状改善，呈负相关。

因此，本发明通过计算帕金森病患者脑电信号中3-6Hz频带能量来描述患者症状严重程度，并以此为依据实时调整刺激参数，实现反馈式闭环刺激。该刺激方法具体实现步骤如下：

S1、针对实时采集的帕金森病患者STN脑电信号进行处理，需要每采集1秒处理一次，防止采样率过高时信号占用内存过大，影响计算效率；

S2、对S1中所采数据进行3-6Hz带通滤波；频带选择3-6Hz是因为：

第一、任何带通滤波器都会存在滚降现象，即滤波器并不能够将所设截止频率外的所有频率完全衰减掉，比如低端截止频率设为2Hz的时候，信号中依然出现1Hz的能量，因此本发明低端截止频率设为3Hz，是为了防止信号中频率为0Hz成分的直流分量衰减不完全，直流分量通常较大，影响计算结果。

第二、某些帕金森病患者STN脑电信号中会出现持续的alpha8-13Hz能量，虽然目前尚不清楚原因，但筛选频带时需避开此频带，同时考虑到滚降现象存在，因此高端截止频率设为6Hz。

S3、对滤波后的信号进行降采样处理，将采样率降为250Hz；降采样的目的是为了减少数据量，提高计算效率，该步骤一定要在滤波后进行，否则会引起高频混叠。

S4、计算降采样后信号的线长ll，表示信号能量，其中x为信号，N为信号中样本点数量，计算公式如下：

S5、重复5次S1-S4，得到5秒钟的5次每一秒的线长ll，进行中值滤波代表此5秒内3-6Hz的能量p，计算公式如下：p＝median(ll

S6、将能量p与设置阈值thr对比，若p

本发明的有益效果是：

1、本发明基于帕金森病患者服用药物以后低频段能量会上升这一现象，选择3-6Hz频段能量变化来描述患者症状变化情况，既能排除信号中直流分量干扰，又可以排除某些帕金森病患者STN脑电信号中出现的alpha振荡，同时考虑到了滤波器的滚降现象，解决了传统帕金森闭环反馈方法中病理性beta频带难以定位问题。

2、本发明在处理信号过程中，每次只处理1秒钟数据，并进行降采样处理，保证了运算速率，提高了实时性。每隔5秒进行一次中值滤波，滤除信号中波动较大的干扰，提高准确性。此外，整个过程占用内存少、计算简单，适合IPG这种算力有限的系统。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是A帕金森病患者STN脑电时频图和PSD；

图2是A帕金森病患者服用多巴胺能药物后STN脑电信号的时频图和PSD；

图3是B帕金森病患者STN脑电信号时频图；

图4是存在多条β频带的STN脑电时频图；

图5是C帕金森病患者STN脑电信号中的alpha振荡；

图6是实施例一中帕金森病患者1秒钟STN脑电信号带通滤波后波形；

图7是实施例一中降采样后脑电信号波形；

图8是实施例二中帕金森病患者药物代谢周期脑电信号的时频图和反馈算法值变化曲线图；

图9-17是实施例三中9名帕金森病患者药物代谢周期STN脑电信号分析结果。

具体实施方式

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

基底神经节局部场电位中beta活性增强是帕金森病的典型特征，图1是A帕金森病患者一段STN脑电信号的时频分析图和功率谱密度(Power Spectral Density，PSD)，图中，上图为时频图，横坐标表示时间，纵坐标表示频率，颜色表示某一时刻某个频率的能量大小，颜色越深能量越大，可见此段信号中有非常明显的病理性beta能量，频带范围为13-18Hz。截取其中一段数据进行功率谱密度分析，时长一分钟，如图1下图所示，图中横坐标表示频率，纵坐标表示频率所对应的能量，依然可见很明显的beta能量，其峰值频率为15Hz。

图2为A帕金森病患者服用多巴胺能药物后STN脑电信号的时频图和PSD，患者约于18:00服用药物，约18:39药效开始生效，由时频图可知，药效生效后发生如下现象：①13-18Hz频带颜色变浅；②低频段颜色变深。表示信号中病理性beta(13-18Hz)能量下降以及低频段能量上升。截取药物生效后一段信号进行功率谱密度分析，时长一分钟，如图2下图所示，对比图1中PSD可以发现，信号中beta能量下降明显，峰值频率由15Hz变为2.5Hz，可见低频段能量明显上升。

目前已有文献中多以帕金森病患者脑电信号中病理性beta能量变化作为反馈控制信号，但是研究表明会存在一些问题：第一，不同帕金森病患者的病理性beta频段不同，可能与患者症状不同或症状严重程度不同有关，图3是B帕金森病患者STN脑电信号时频图，由图可知，此患者病理性beta范围在18-22Hz之间，与图1中13-18H明显不同。第二，某些患者脑电信号中可能会出现多条beta频带，难以确定哪一条是病理性频带，如图4所示，图中清晰可见两条beta频带，一条低beta、一条高beta。第三，大脑beta频带与运动相关，运动或者运动想象会导致大脑中beta活性增强，影响反馈结果。第四，静息性震颤症状与beta活性无关。

上述几点问题会导致定位病理性beta频带困难，影响闭环刺激效果。通过上文分析可知，服用多巴胺能药物以后，除了病理性beta能量下降，同时也会伴随低频段能量上升，因此本发明以3-6Hz频带能量变化作为反馈控制信号。

闭环刺激方法流程如下：

S1、实时采集帕金森病患者STN脑电信号，每采集1秒处理一次。每1秒处理一次是为了防止采样率过高时信号占用内存过大，影响计算效率。

S2、对S1中所采数据进行3-6Hz带通滤波。频带选择3-6Hz是因为：

第一、任何带通滤波器都会存在滚降现象，即滤波器并不能够将所设截止频率外的所有频率完全衰减掉，比如低端截止频率设为2Hz的时候，信号中依然出现1Hz的能量，因此本发明低端截止频率设为3Hz，是为了防止信号中直流分量(频率为0Hz的成分)衰减不完全，直流分量通常较大，影响计算结果。

第二、某些帕金森病患者STN脑电信号中会出现持续的alpha8-13Hz能量，图5为C帕金森病患者STN脑电时频分析图，由图可知信号中出现了频段为8-12Hz的alpha振荡，即图中箭头所指的频带，虽然目前尚不清楚原因，但筛选频带时需避开此频带，同时考虑到滚降现象存在，因此高端截止频率设为6Hz。

S3、对滤波后的信号进行降采样处理，将采样率降为250Hz。降采样的目的是为了减少数据量，提高计算效率。此步骤一定要在滤波后进行，否则会引起高频混叠。

S4、计算降采样后信号的线长ll，表示信号能量，其中x为信号，N为信号中样本点数量。计算公式如下：

其它能用于表示信号能量的方法亦可，如分形维数(Fractal Dimension，FD)或面积(Area)，但不能过于复杂，否则会影响计算效率。

S5、重复5次S1-S4，得到5秒中每一秒的线长ll，进行中值滤波代表此5秒内3-6Hz的能量p，计算公式如下：p＝median(ll

S6、设置阈值thr，若p

实施例，如图6-9所示：

S1，实时采集帕金森病患者STN脑电信号，采样率为500Hz，每采集1秒处理一次。每1秒处理一次是为了防止采样率过高时信号占用内存过大，影响计算效率。

S2，对S1中所采数据进行3-6Hz带通滤波，滤波后波形如图6所示。

S3，对滤波后的信号进行降采样处理，将采样率降为250Hz，降采样后波形如图7所示，对比图7和图6可知，降采样后信号波形形态并未发生变化，但横坐标采样点数却变为原来的一半，便于后续计算。降采样的目的是为了减少数据量，提高计算效率。此步骤一定要在滤波后进行，否则会引起高频混叠。

S4，计算降采样后信号的线长ll，表示信号能量，其中x为信号，N为信号中样本点数量。计算公式如下：

此实施例中，N＝250，x波形见图7，因此能量ll即为：ll＝(|(-0.458)-(-0.472)|+…+|0.123-0.144|)＝12.85dB，表示此1秒钟信号内3-6Hz能量为12.85dB。

S5，重复5次S1至S4，得到5秒中每一秒的线长ll，进行中值滤波代表此5秒内3-6Hz的能量p，计算公式如下：p＝median(ll

S6，设置阈值thr，若p

实施例二：

持续采集并分析此患者一个药物代谢周期脑电，包括服用药物前、药效生效时以及药效失效后，以验证算法有效性，结果如图8所示。图中，上图为信号时频图，下图为反馈算法值随时间变化曲线，计算方法重复步骤一到步骤五。患者约于10:30睡觉，约于11:15药物生效，约于14:00药效失效。由时频图可知，药效生效后，信号中病理性beta能量下降明显，同时3-6Hz能量上升；而当药物失效后，信号中病理性beta能量又重新上升，同时伴随3-6Hz能量下降。从图8下图可知，当药物生效时，算法值上升并维持在较高水平，而当药物失效后，算法值又重新回到基线，基线是指帕金森病患者症状较差时低频带能量水平。因此，上述反馈算法可以很好地描述出帕金森病患者在药物代谢周期中病理性beta活性的变化趋势，设置阈值thr为20dB。当反馈算法值低于阈值时，信号中病理性beta活性较强，患者症状较差，应开启刺激；当反馈算法值高于阈值时，信号中病理性beta活性下降，患者症状改善，此时应减小或关闭刺激。

实施例三：

为进一步验证上述方法有效性，采集16名帕金森病患者药物代谢周期STN脑电信号并进行时频分析和反馈算法计算，图9-17展示了其中9名患者分析结果。结果显示，药物生效时，患者症状改善，3-6Hz能量上升，反馈算法值上升；药物失效时，患者症状变差，3-6Hz能量下降，反馈算法值下降，因此本发明实现的闭环刺激方法可以针对帕金森病患者症状变化时脑电信号中频带能量变化的典型特征进行描述，根据患者症状严重程度实时调整刺激参数，实现反馈式闭环刺激。