孙玉波，李海涛，舒智林，于洋，韩建达，梁思泉，于宁波

（1.南开大学 人工智能学院, 天津 300350; 2.南开大学 天津市智能机器人技术重点实验室, 天津 300350; 3.天津市环湖医院 神经外科, 天津 300350; 4.天津市环湖医院 康复医学科, 天津 300350）

帕金森病(Parkinson’s disease, PD)是一种常见的神经退行性疾病[1]。我国是世界帕金森第一大国，且病人数量随人口老龄化持续上升。由于脑部黑质多巴胺能神经元的进行性缺失，导致患者出现运动症状如步态障碍和震颤，非运动症状如抑郁和焦虑等，严重影响患者的生活和自理能力[2-5]。

脑深部电刺激(deep brain stimulation, DBS)已被证明是一种有效的帕金森病治疗手段[6-8]。术后DBS神经调控参数的调整需要量化、持续、精确的评估，对患者运动症状的动态测量与日常生活中的跟踪测量也有助于评估患者的病情进展情况[9-10]。但是，目前这主要依赖于医生的经验和判断，迫切需要客观、量化的测评手段和分析方法。

在人体步行过程中足跟击地时，腓肠肌的激活是保证脚踝稳定的重要因素，所以人体步行时腓肠肌表面肌电信号(surface electromyogram signal, sEMG)具有周期性特征[11]。根据腓肠激活时刻的检测，可以对足跟击地时刻进行识别，并提取步态周期时间信息。步态障碍是帕金森患者的典型运动症状[12]，步态量化测评对于帕金森病筛查、鉴别诊断和治疗评价都具有非常重要的意义。有研究表明患者下肢腓肠肌肌电信号变化与患病程度密切相关[13-14]，因此帕金森患者的腓肠肌表面肌电信号分析对患者步态量化测评也具有重要意义[15-16]。

DBS神经调控治疗，是对患者大脑直接施加电刺激，这可能对与步态直接相关的腓肠肌的神经指令造成影响。在步态研究的前向动力学方法中，神经指令作为输入，直接确定了肌肉激活的幅度。这可以通过从肌电信号中提取的神经激活度与肌肉激活度进行量化[17]。

帕金森患者脑部黑质多巴胺能神经元的进行性缺失，会导致基底神经节输出过度。基底神经节下行调节脑干活动，脑干是负责调节人体运动节律的重要部分，因此，运动节律失调是帕金森患者步态障碍的一个主要特征[18]。为了调整基底神经节的输出，大部分针对帕金森患者的DBS电极刺激靶点为基底神经节中的丘脑底核[8]。针对帕金森的病理特征和DBS神经调控治疗的特点，本文提出与步态节律密切相关的腓肠肌神经激活度与肌肉激活度变异性指标。

足底压力可以充分反映步态。地面反作用力是人在行走过程中受到的唯一外力。在支撑相，足底与地面相互作用，地面反作用力对于步态有着根本性的作用，该信息可以通过力传感器采集[19-21]。本文通过足底压力进行步态在时域的准确识别，提取在临床广泛应用的步态周期时间变异性指标，对所提基于表面肌电信号的变异性指标进行验证。

本文针对帕金森患者脑深部电刺激神经调控治疗步态测评需求，分析步态与帕金森病理特点，提出基于外侧腓肠肌(lateral gastrocnemius,LG)表面肌电信号的步态分析方法和量化测评指标，建立一个肌电和足底压力信号同步采集与无线传输系统，并开展注册临床研究进行有效性验证，为帕金森患者提供了一种轻量化的、基于外侧腓肠肌sEMG的步态量化测评方案。本文的工作特色和创新具体如下：

1) 针对步态与帕金森病理特点，本文提出了基于外侧腓肠肌表面肌电信号的步相识别方法与步态变异性指标提取方法。采用平滑先验算法去除肌电信号中的低频趋势并保留帕金森临床症状相关信息，通过小波变换与激活程度识别，检测足跟击地时刻、进行步相识别，进而对肌电信号进行划分，计算神经激活度、肌肉激活度、及其变异性指标。

2) 帕金森患者的运动障碍表现极易受到仪器和场景等带来的心理暗示的影响，为客观测评带来挑战。为此，本文采用无线通讯、穿戴式、轻量化的肌电和足底压力采集设备，建立数据同步采集方法，设计基于在临床广泛应用的10 m步行测试的实验方案，构建面向帕金森患者的步态量化测评系统和范式。

3) 取得伦理许可，针对接受DBS治疗的帕金森患者开展注册临床研究，基于时间同步的足底压力信号对所提方法和指标的有效性进行了验证。证明了本文所提步态量化测评方法的有效性和适用于帕金森患者的步态量化测评系统的可行性，为帕金森神经调控治疗提供重要临床手段。

1 系统设计与实验范式

1.1 系统设计

人体在步行过程中受到的唯一外力是地面反作用力，而人的运动由肌肉驱动，因此肌电信号和足底压力中包含着步态的关键信息。基于帕金森的病理特点，病人步态障碍极易受到声音、视觉等提示的影响。为了减小对病人的心理暗示和干扰，实现对患者自然真实、客观量化的步态测评，本文构建了一个采用无线通信、穿戴式、轻量化的肌电和足底压力信号同步采集系统。系统设计示意图如图1所示。

图1 系统示意图Fig.1 System diagram

本文采用Delsys trigno avanti设备(delsys incorporated, USA)记录下肢外侧腓肠肌表面肌电信号，采样频率为2 000 Hz。采用无线压力鞋垫(ZIGUN AVOIN company, Finland)采集患者在行走过程中的足底压力，采样频率为100 Hz，最大测量值为1000 N, 分辨率为0.1 N/cm，压力鞋垫质地柔软且厚度不超过1.2 mm，可以方便、舒适地放到病人的鞋中。

两种信号采样频率不同，并且是由两台独立的计算机分别进行采集。为此，需设计一个能够实现信号同步的打标软件模块，如图1所示。图1中左侧虚线箭头表示的均为无线连接通信，而中间和右侧实线箭头表示的均为有线连接通信。

同步软件单独在一台计算机中运行，程序为Windows Form窗体程序，在Microsoft Visual Studio 2019中基于C#语言开发。通过串口分别连接到表面肌电和足底压力设备对应的计算机。按对应的波特率（肌电：110；压力鞋垫：115 200）和数据位（肌电：5；压力鞋垫：8）向对应设备发送一个字节数据，作为同步信号。

1.2 实验范式

基于在临床广泛应用的10 m步行测试，设计测评实验范式。受试者在起点处静止站立，在听到开始指令后向前自然行走，抵达终点后停止。通过无线足底压力测量单元和表面肌电测量单元对受试者的步态数据进行采集。

步行过程可分为3个部分：启动、步行和停止。对于神经功能障碍患者，其启动和停止部分的步态易受到多种因素影响。因此，本文将每次测试的前两个步态周期（启动）和最后两个步态周期（停止）的数据剔除，仅采用中间步行部分的数据进行分析，如图2所示，从左到右，步行过程划分为3个部分：启动、步行和停止。

图2 基于标准临床10 m步行测试的实验范式Fig.2 Experimental paradigm based on the standard clinical 10-meter walking test

2 信号处理与特征提取

通过建立的肌电和足底压力信号同步采集系统，在临床10 m步行测试实验过程中同步采集数据，通过无线传输至计算机，随后分别对两种信号进行分析与特征提取，流程如图3所示。

图3 步态变异性指标提取流程Fig.3 Diagram of gait variability features extraction

对时间对齐后的表面肌电信号进行滤波等预处理，然后基于小波变换进行足跟击地时刻识别，提取步态周期时间，随后对肌电数据进行划分并提取神经激活度、肌肉激活度及其变异性指标。基于足底压力数据,在时域对步态进行准确识别，进而计算步态周期的变异性系数，验证肌电信号变异性指标的有效性。

2.1 外侧腓肠肌sEMG信号分析与特征提取

在表面肌电信号预处理过程中，需要去除肌电信号的运动伪迹以及由电极与皮肤的不稳定接触造成的低频趋势。通常做法是对表面肌电信号进行截止频率不超过20 Hz的高通滤波。然而，对帕金森患者来说，肌电信号低频成分中存在着与其临床特征相关的重要信息。因此，在去除低频趋势的同时，要避免相关信息的丢失。为此，本文采用平滑先验算法。

对时间对齐后的表面肌电信号的预处理流程如图4所示，分为3个步骤：通过平滑先验算法去趋势、对肌电信号进行幅值归一化处理、全波整流。

图4 表面肌电信号预处理流程Fig.4 Preprocessing workflow for sEMG signals

假设表面肌电信号的时间序列由有效信号和低频趋势成分两部分组成：

式中：zs是去掉趋势项之后的表面肌电信号序列；zt是表面肌电信号中的低频趋势成分。

对有效信号的估计可以表示为

式中 α为正则化参数。D2的取值为

通过平滑先验算法去趋势后，进行幅值归一化处理。归一化尺度一种可行的选择是表面肌(surface electromyogram, sEMG)信号绝对值中的最大值[22]。但为了排除异常值和采集过程中由于运动而产生的高幅值异常值对激活时刻识别的影响，归一化尺度参数选择为肌电信号50个最大绝对值中的最小值[23]。经过全波整流后，得到去除了基线漂移、低频趋势和高频噪声的肌电信号；再通过低通滤波器得到线性包络。

在人体步行过程中足跟击地时，腓肠肌的激活是保证脚踝稳定的重要因素。外侧腓肠肌表面肌电信号具有与足底压力信号相同的周期性特征，由传感器采集到的原始足底压力信号与同侧下肢的外侧腓肠肌表面肌电信号如图5所示。本文提出了一种基于小波变换和阈值设定的步相识别方法。通过对外侧腓肠肌的激活状态识别，确认足跟击地时刻，进而识别步相。

图5 一名帕金森患者的时间对齐的足底压力原始信号与外侧腓肠肌表面肌电原始信号Fig.5 Synchronized raw data of plantar pressure signal and LG sEMG signal in a PD patient

小波变换可以同时保留时域和频域的信息，有研究表明小波变换在鉴别肌电信号模式方面优于其他时频方法，并且相比传统方法提供了更好的时频分辨率，更加适用于不稳定的肌电信号分析[24-25]。

使用离散小波变换代替均方根方法，对表面肌电信号进行处理，可以更加准确地识别肌肉激活时刻，去除噪声影响。经过离散小波变换，高频噪声分量主要集中在细节系数中，而近似系数主要表示信号的低频信息。在识别激活时刻过程中，主要关注的是肌电信号幅值的低频变化，关注低频信息。在肌肉激活程度识别过程中，离散小波变换方法具有更高的准确度[26]。

基于小波的方法是分析非平稳信号如肌电信号的实用技术。离散小波变换(discrete wavelet transform, DWT)对母小波进行缩放和移位，并将离散时间信号分解成一组信号。DWT使用高通和低通滤波器与降采样方法，将信号分解成多分辨率系数。原始肌电信号S经低通滤波器和高通滤波器，获得第一级的近似系数(cA1)和细节系数(cD1)。为了得到多分辨率子集，对上一级的近似系数进行降采样和滤波，并重复此操作，获得近似系数 cAn，细节系数 cD1，cD2，…，cDn。离散小波变换分解流程如图6所示。

图6 表面肌电信号分解树Fig.6 Tree decomposition of the sEMG signal

相关研究表明基于Haar小波的DWT可以更加有效地对不同动作进行识别，对肌肉活动检测准确度较高[26]。并且基于Haar小波的小波变换方法可以更有效地去除肌电信号中的噪声[27]。因此本文中，对预处理后的外侧腓肠肌表面肌电信号进行离散小波变换，采用Haar小波作为母小波，进行5级分解，获得近似系数cA5。随后通过3 Hz的巴特沃斯低通滤波器，获得近似系数的线性包络，再通过阈值方法，检测激活状态，从而实现对足跟击地时刻的识别。通过外侧腓肠肌表面肌电信号进行足跟击地时刻识别的具体流程与对应信号示例如图7所示。

图7 帕金森患者外侧腓肠肌sEMG信号激活识别处理过程Fig.7 Schematic diagram of LG sEMG activation recognition process in PD patients

本文中，使用了已经得到广泛认可的肌肉活动检测阈值设定方法[28]，阈值T由式(4)所得。其中，本文取患者在静止站立时长度为2 s的外侧腓肠肌表面肌电信号作为基线，预处理后提取其平均值M和标准差STD。将近似参数包络线中超过阈值T的时刻定义为激活时刻。从非激活向激活过渡的时刻定义为足跟击地时刻。

随后，将经过预处理与低通滤波求出的线性包络e(t)，按照足跟击地时刻进行划分。再通过神经激活模型计算每段信号的肌肉神经激活度u(t)。

2.2 神经激活度与肌肉激活度

肌电信号的线性包络e(t)转化为神经激活度u(t)的过程称为激活动力学。当肌肉纤维被单个动作电位激活时，肌肉就会产生抽动反应。这个响应可以用临界阻尼线性二阶微分系统很好地表示，如式(5)所示。这种类型的响应是根据肌电信号的输入e(t)来确定神经激活度u(t)的基础。

式中M、B和K是定义二阶动力学系统的常数。在实验室中采集到的数据是离散的时间序列信号，离散方程表示为

式中：d是采集设备的延迟，由于所采用设备延迟为48 ms，采样周期0.5 ms，所以d取96； β1和 β2是定义二阶动力学系统的参数。这些参数将e(t)映射到神经激活度u(t)。参数 α、 β1和 β2满足以下条件：

为了表征神经激活度与肌肉激活度之间的非线性关系，引入一个函数：

式中a(t)是计算所得的肌肉激活度。

2.3 足底压力数据处理与步态识别

将时间对齐后的足底压力信号通过一个截止频率为20 Hz的低通滤波器，消除高频噪声。基于在足底压力信号中识别出的周期性特征进行步态周期识别。

在进行步态识别和划分时，考虑到患者体重的特异性以及足底压力信号在摆动相过程中具有低频噪声，设定区分摆动相和支撑相的阈值为体重的5%。记录每个步态周期摆动相的结束时刻，并作为下一个步态周期支撑相的开始时刻，如图8所示，为一名PD患者的足底压力曲线与步态识别示意图。随后可提取出每个步态周期的时间长度。

图8 基于足底压力信号的步态识别Fig.8 Gait recognition based on plantar pressure

在本文中，根据足底压力信号划分步态周期，计算得到步态周期时间变异性指标，与表面肌电信号变异性指标进行交叉验证。

2.4 步态变异性指标

运动节律失调是帕金森患者步态障碍的一个主要特征，并且DBS神经调控可能会对与步态相关的外侧腓肠肌的神经指令造成影响。因此，确定步态变异性作为主要测评指标。

步态周期时间的变异性是公认可以反映帕金森患者运动障碍的评估指标。具体计算方法如式(9)所示，由步态周期时间的平均值和标准差计算得到：

计算肌电信号激活度变异性时，令步态周期个数为n。求得每个步态周期相应的平均肌肉神经激活度u(k)，k∈[1,n]。

式中：t∈[1,mk]；mk为第k个步态周期的样本长度。计算神经激活度变异指标：

对肌肉激活度a(t)计算变异指标的过程与u(t)相同：

3 临床实验与结果分析

3.1 伦理与受试者招募

本研究获得天津市环湖医院伦理委员会批准(2019—35)，并在中国临床试验注册中心注册(ChiCTR1900022715)。入组标准为：

1） 需要重新调整DBS治疗参数；

2） 年龄40岁及以上；

3） 学历小学及以上；

4） 能够自主站立和行走；

5） 无其他会对步态造成影响的肢体或肌肉相关疾病。

本研究中，5名接受了脑深部电刺激手术、并符合入组标准的帕金森患者被招募，其中包括2名男性和3名女性，基本信息如表1所示。所招募的受试者平均年龄61.8岁，平均身高165.0 cm，平均体重67.6 kg。

表1 招募病人的基本信息Table 1 Information of PD patients

实验开始前，每位帕金森病患者都被充分告知实验目的和步骤，并签署了同意书。患者的临床特征、所使用的DBS仪器种类、初始DBS刺激参数以及神经调控过程中所使用的刺激频率等详细信息见表2。

表2 病人DBS神经调控初始参数Table 2 Patients’ initial DBS programming parameters

3.2 步态量化测评结果与分析

本文对5名帕金森患者在接受DBS神经调控时的步态变异性进行了量化测评。根据表面肌电信号，计算出外侧腓肠肌肌电信号变异性指标。根据足底压力信号，识别步态周期，计算步态周期时间的变异性指标。在不同DBS神经调控频率下，患者的步态变异性指标有显著不同，具体结果如表3所示。

表3 患者神经调控测试中提取的步态周期时间变异性指标与外侧腓肠肌肌电信号变异性指标Table 3 CV of step time and GL sEMG under different DBS parameters

在由步态周期时间变异性指标与外侧腓肠肌变异性指标构成的二维空间中，特征点变化呈现出一定的趋势，如图9、10所示。图9、10中蓝色直线为最小二乘法线性拟合结果，以表示二者变化趋势，阴影部分为95%置信区间。由此可见，患者的特征点离原点越近，即步态周期时间变异系数与所提肌电变异系数越小，患者步态变异性越小，步态越稳定。所提指标可以为步态量化测评提供参考。

采用皮尔逊相关性系数分析所提出的神经激活度和肌肉激活度变异性指标与步态周期时间变异性之间的关系，相关性系数分别为0.602(p=0.000)和0.591(p=0.000)，具有显著相关性。由此，说明了本文所提根据外侧腓肠肌肌电信号提取神经激活度与肌肉激活度变异性指标的方法的可行性，确切证实了所提变异性指标的有效性。且本文所提的2种指标相较于传统的步态周期时间变异性指标，对患者的状态具有更好的分辨率。如图9、10所示，步态周期时间变异性的变化范围较小，而所提表面肌电信号变异性指标则具有更大的变化范围，可以更好地对步态变异性进行分辨。

图9 外侧腓肠肌神经激活度变异性指标与步态周期时间变异性指标呈线性趋势Fig.9 Linear regression between the CV of step time and LG neural activation

图10 外侧腓肠肌肌肉激活度变异性指标与步态周期时间变异性指标呈线性趋势Fig.10 Linear regression between the CV of step time and LG muscle activation

在患者的连续稳定步态周期中，基于肌电信号提取的神经激活度、肌肉激活度在根据步相划分后，对患者每个步态周期的量化，如图11所示。图11左侧为患者在连续稳定步态过程中的足底压力信号与对应的神经激活度、肌肉激活度曲线，右侧为对应的计算变异性指标的步态周期时间、平均神经激活度和平均肌肉激活度柱状图，与对应提取的变异性指标。相较于步态周期时间长度，所提指标具有更好的区分度。

图11 一名帕金森患者的足底压力信号、神经激活度曲线、肌肉激活度曲线与提取出的变异性指标Fig.11 Plantar pressure signal, neural activation curve and muscle activation curve and the extracted variation index of a PD patient

针对帕金森患者神经调控步态量化测评中步态变异性的评估，所提的基于外侧腓肠肌sEMG信号采集、无线传输与分析的步态量化测评方法，能为医生提供更加有区分度的指标，为临床诊断提供辅助。

4 结束语

本文针对帕金森患者脑深部电刺激神经调控治疗步态量化测评需求，分析帕金森病理与步态特点，提出了一种基于外侧腓肠肌sEMG数据的步态量化测评方法。建立一个肌电和足底压力信号同步采集与无线传输系统，基于在临床广泛应用的步行实验范式，开展临床研究并进行有效性验证。通过时间同步的足底压力信号，验证了所提测评方法和指标在帕金森病患者步态量化测评中的有效性。且所提神经激活度与肌肉激活度变异性指标，相较于传统指标具有更好的区分度。由此，为帕金森神经调控治疗提供了一种重要临床手段。