基于下肢肌张力感知的痉挛特征提取算法研究

2022-01-11汪步云魏壮壮许德章

兵器装备工程学报 2021年12期

汪步云，吴臣，魏壮壮，张振，许德章

(1.安徽工程大学机械工程学院，安徽芜湖 241000； 2.安徽工程大学人工智能学院，安徽芜湖 241000； 3.芜湖安普机器人产业技术研究院有限公司，安徽芜湖 241007)

1 引言

从被动辅助行走训练到人体下肢运动功能复原这一过程，对脑卒中患者的康复疗效尤为重要；外骨骼机器人辅助下肢运动障碍患者开展步态训练已经成为一种有效的康复临床治疗手段[1]。在此治疗过程中，针对脑卒中患者，较少考虑由于训练强度与康复状况的不良匹配而诱发的痉挛及其带来的二次伤害[2]。当前，在步态康复训练中，融入痉挛在线的客观量化评价是康复医学上的迫切需求。由于痉挛发生存在较大的随机性，这将给患者带来极大的生理和心理压力，不利于患者接受康复治疗。因此，在线量化评定康复过程中伴随出现的病理性痉挛是康复临床的迫切需求，预防痉挛及防范痉挛带来二次损伤亦是康复机器人亟需解决的关键问题。

病理性痉挛是脑卒中患者在步态康复时最为常见的并发症[3]，在临床上往往表现为肌张力异常增高，外在特征主要为肌肉间断性或连续性的抽搐并伴随不同严重程度的疼痛[4]。针对痉挛检测及其特征提取，国内外一些学者和研究机构开展了相关研究。中国科学技术大学的朱晓斐等[5]采用表面肌电与惯性传感器相结合的方法，开展了人体上肢的肌痉挛量化评估，在临床应用中引入了客观量化的评估参数，但在交互力激励下的肌痉挛评定有待进一步开展。张敏等[6]设计了一种线阵CCD(charge coupled device)数字差分式痉挛传感器，仅适用于微张力的测量。Niels Buchhold等[7]开发了一种应变式传感器用于痉挛评测，但临床应用的实验验证未见报道。Hu Baohua等[8]针对上肢痉挛在临床上不易量化评定的问题，在时域内对获取的痉挛信号采用了数学变换的方法，提升了临床客观评定痉挛的应用效果。Heung HL等[9]针对脑卒中患者，设计了柔性康复驱动器，引入了客观量化的痉挛在线评定，有效提升了安全性。上述研究取得了有意义的研究成果，但对于运动功能的恢复过程，预警病理性痉挛、评估痉挛发生的具体类型则较少涉及。

本文依据病理性痉挛在康复临床中的评定机理，设计了痉挛传感器，通过获取与运动关联度较高的下肢肌肉区域的肌张力信息，给出时域与频域的信号特征表达；针对病理性痉挛中阵挛和痛性痉挛信号不易量化与特征阈值随机出现的问题，提出BFT算法，解析特征信息及痉挛映射关系，开展实验研究，验证了算法对痉挛特征评定的可行性。

2 下肢病理性痉挛的感知机理

2.1 下肢运动的病理性痉挛评定

脑卒中患者的下肢相关肌群在运动过程中，由病理性痉挛带来的肌肉收缩并非随意，即痉挛性肌张力信号具有一定规律，肌张力信息与病理性痉挛特征是关联的，这为病理性痉挛的特征信息感知提供了前提条件[10]。因此，检测下肢特定肌肉区域的肌张力信息可以获得病理性痉挛的相关特征信息。根据临床经验总结，得出患者在运动功能恢复的这一过程中，病理性痉挛较大概率发生在下肢小腿肌群中的伸肌与屈肌群，即小腿的胫骨、腓肠、腓骨等区域。因此，应该重点检测下肢小腿后侧比目鱼肌、大腿前侧腓肠等小腿相关肌群。

下肢肌肉的病理性痉挛，主要包括强直、阵挛、痛性痉挛[11]。针对上述痉挛在临床上的表现，痉挛性肌张力幅值突然增加是强直的外在表现；肌张力的频率有较大变动为阵挛表现；而痛挛则是上述两者特征的同时兼备，并且发生的随机性增大。

痉挛产生的内在机理和病理性特征具有不确定性，表现形式多样，评价方法也不统一。对于患者的病理性痉挛症状，不同的医师会得出不一致的结论。由于病理性痉挛特征阈值随机出现，痉挛性肌张力信号在时域与频域上信号特征存在交替变化的特点，阵挛和痛性痉挛信号不易量化。因而，在痉挛传感器设计基础上，匹配BFT特征提取算法，给出阵挛与痛性痉挛发生时肌张力在频率谱上的特征信号表达，通过实验探求特征阈值表达的有效性，优化评定病理性痉挛方法，将匹配BFT算法的传感器与外骨骼机器人相结合使用，在外骨骼式下肢康复机器人的步态训练中融入客观量化的痉挛在线评价。据此，本文总体思路与评定过程如图1所示。

2.2 痉挛传感器设计

依据上述检测机理，在人体下肢佩戴痉挛传感器时，下肢对应检测区域的肌张力可转换为传感器机构与人体肌肉接触力作用在传感器内部的压敏单元上的压力比值。通过配置特征提取算法，进而转化为传感器获取的肌张力大小及其变化的频率。图2表示痉挛传感器装配图和实物图，展示了机械结构、处理电路及压敏电阻模块，传感器材料组成主要为铝，其总质量较轻，约为248 g。主要技术参数如表1所列。

图1 下肢肌肉痉挛的评定过程框图Fig.1 Assessment process of lower limb muscle spasm

1～15依次为上调节螺钉、上压紧端盖、处理电路、下调节螺钉、固定螺钉、压敏电阻、压敏电阻安装凸台、顶柱板预警弹簧、导向轴、顶柱、滑动导套、传感器外壳、衬套、下压紧端盖

图2 痉挛传感器Fig.2 Spasm sensor

表1 痉挛传感器主要技术参数Table 1 Main technology parameters of spasm sensor

痉挛传感器电路设计实物图和压力传感器如图2所示，痉挛传感器硬件电路原理设计简图如图3所示，STM32具有丰富的外设功能和数据处理能力[12]。以STM32F103RCT6处理器构成的最小单片机系统为核心，扩展了包括以NRF24L01为核心的无线数据通讯模块、CP2102为核心的串口通讯模块，以及压敏电路等，通过采集压敏电阻的模拟量，解析接触力幅值及频率，以转换为痉挛性肌张力信息。痉挛传感器硬件电路的主要功能模块及接口如表2所列，压敏电阻采用FSR402电阻式压力传感器，其实物与性能测试曲线如图4所示。

图3 痉挛传感器硬件电路设计原理图Fig.3 Schematic diagram of hardware circuit design for spasm sensor

表2 痉挛传感器硬件电路功能模块及接口Table 2 Description of function module of hardware circuit of spasm sensor

2.3 痉挛传感器的稳定性判别

假定频率ω为剪切频率ωc(其值大于零)时，相频特性曲线距离-180°线的相位差γ为相位裕度。表达式为

γ=180°+φ(ωc)

(1)

当ω是相位交界频率wg时，开环幅频特性|G(jω)H(jω)|的倒数，叫做幅值裕度Kg。表达式为

(2)

Bode图一般用对数坐标轴表达，幅值裕度为Kg:

Kg=-20lg|G(jω)H(jω)|

(3)

若开环稳定系统，当Kg和γ为正数时，闭环系统亦是稳定的。用Matlab软件绘制出痉挛传感器开环控制系统Bode曲线，如图5所示。

图4 FSR402传感器实物(上)及其性能曲线(下)Fig.4 FSR402 sensor and performance curve

图5 传感器开环控制系统Bode曲线Fig.5 Bode diagram of open-loop control system for the sensor

从图中可知，Kg=53.6 dB，γ=11.8 deg(均大于零),因此传感器闭环控制系统也是稳定的。

肌张力信号的采集，其采集流程如图6所示。传感器采集下肢肌肉痉挛信号传递给上位机，信号经放大后传递给主机模块，经处理后将数据通过无线通讯的方式传递给上位机模块，综合解析后显示出痉挛的等级。

图6 肌张力信号的采集流程框图Fig.6 Collection flow chart of muscle tension signal

2.4 痉挛性肌张力信号的预处理

通过检测得到痉挛传感器的初始信号含有不同噪声，例如，痉挛传感器的压敏单元转换噪声、生理因素引起的噪声等。需要对信息滤波处理，以便为后续特征提取提供有效测试数据。移动均值滤波法是线性滤波法的一种，对于肌张力这类不平稳且具有随机性生理信号的平滑滤波非常适用。该方法计算过程简便，对原始信号的信息量能够较好保存。在分析肌张力信号特点规律和对比生理信号预处理滤波方法基础上，本文应用移动均值滤波算法对原始肌张力信号展开预处理[13]，其具体表达如下：

(4)

式(4)中，M(s)是第i个信号值的平均值，T是邻域范围。对典型的痉挛肌张力信号，运用移动均值滤波法展开处理，其结果如图7所示。通过对比分析得出，选取窗口值T=15，如图7(c)所示，信号更平滑，滤波效果更好。由于痉挛对时效性有较高要求，所以肌张力信息可以使用移动均值滤波得到整体性的保存，同时也保留了痉挛性肌张力信息中的频率跳变，可以更准确地提取出痉挛发生前后变化及特征阈值。

图7 各窗口下的移动均值滤波曲线Fig.7 Filtering diagram of moving mean under each window

3 下肢痉挛性肌张力的时域特征表达

在康复临床应用领域，肌张力等生理信号特征的时域表达，得到了更加广泛的应用。因此，针对强直、阵挛和痛性痉挛的3种临床分类，给出3种痉挛在时域上的表达与发生时的特征阈值。

通过检测脑卒中患者下肢腓肠内肌等部位的肌张力信号，实时获取患者下肢痉挛性肌张力实验数据，经移动均值滤波后，提取痉挛前后阶段的肌张力信号特征，综合解析后给出痉挛评定，然后求得时域的特征表达。同时，利用统计学信息，如最值、均值、绝对积分平均值(integral absolute mean value，IAV)等，对比分析痉挛前后2个阶段肌张力信号，提取反映痉挛性肌张力信号在时域表达的统计学特征[14]。

3.1 强直下的痉挛时域特征提取

强直作为病理性痉挛的一种，其外在症状为下肢痉挛区域肌肉强烈且连续性性的收缩，并且伴随着痛感，同时肌张力的幅值变化快速升高，并保持较高值。强直的出现具有不确定性，故人为施加交互力激励模拟强直的产生。从传感器信息获取的角度看，稳定的矩形波为信号的具体表现，在跳变点处信号会出现突跃。据此，在强直将要出现时，提取肌张力信号幅值T0(N0个采样点数据的均值，N0为0～500 ms时刻测得的信号数据)；在强直产生过后，信号会发生突然跃升，并且信号在突然跃升后其幅值保持在一定的平稳状态，再次提取这一信号幅值Tt(N1个采样点的均值，N1为500～2 400 ms时刻测得的信号数据)，记二者的差值为Tq，计算式如下：

(5)

对痉挛性肌张力采集到的数据，进行处理后如图8所示。患者下肢肌肉皆为正常的状况下，信号值在小于500 ms范围内基本平稳，当信号值大于500 ms时，肌张力值会发生突然跃升，并且保持在1 600 g左右，利用式(5)，得到Tq的均值为1 360 g，若超过此临界阈值，即判定发生强直性痉挛。进一步的，将Tq的平均值看作是强直的特征阈值，并设置为患者出现强直的警界值。

图8 强直下的肌张力信号幅值曲线Fig.8 Amplitude of muscular tension signal under spasticity

测试者及其下肢各区域肌肉的不相同，会使痉挛性肌张力特征阈值的提取存在偶然性的误差。征集了五名测试者参加了强直的时域特征提取实验，测试者的身高和体重各不相同，在同一条件下开展实验，总结得出下肢肌肉强直信号幅值差Tq的均值作为特征阈值，即强直发生的预警值，如表3所示。

表3 下肢肌肉强直预警值

根据实验测出的强直特征，取痉挛出现时的阈值Tq为1 200 g。由式(5)求得强直时域的特征表达式为

F=KaTa，Ka>1，Ta=1 200 g

(6)

式(6)中：Ka为比例系数表征强直严重程度的，F表示肌张力的幅值。

3.2 阵挛在时域的特征提取及分析

阵挛如图9所示,阵挛发生时，肌张力的信号值Fs(t)从小于500 g上升到500～1 500 g的变化范围内，采集数据的时间设为20 s，肌张力信号的均值为1 252 g，上限值为 1 738 g，下限值为352 g，具有周期性特征和明显的规律性。

图9 阵挛性肌张力信号曲线Fig.9 Signals of clonic muscle tension

由以上肌张力信号的特点，可以得出阵挛信号在时域的表达为

F=FZ+Asin(ωt+θ),ω≥6,FZ>500 g

(7)

式(7)中：FZ为阵挛发生时初始幅值，A为阵挛的幅值，ω为阵挛信号频率，θ为初相(阵挛发生时)，F为阵挛下的肌张力信号。幅值A为|F|-FZ,ω与时间有关，当ω≥6时，且Fz>500 g时发生阵性痉挛。

3.3 痛性痉挛在时域的特征提取及分析

痛挛数据采集处理后如图10所示。下肢肌肉张力表现为痛性痉挛时，痉挛信号值Fs(t)在500～2 000 g范围内波动，肌张力上限值和下限值分别为2 092 g、110 g，幅值平均值约为1 202 g。在500 g时呈现出跃升状态，大约在7 000 ms开始出现周期性波动。在2 000～3 000 ms、4 000～5 000 ms两个时间段内，未出现异常状况，表现为间歇性的痉挛状况。由此得出，痛性痉挛表现出并不是很明显的规律性，并且其周期性特征存在间断性。

综合来看，导致痛性痉挛的直接原因是强直和阵挛的相结合。由此可得，痛挛信号的特征表达式为

(8)

若要判定痛性痉挛的发生，需达到强直与阵挛信号在时域的临界阈值。

图10 痛性痉挛时肌张力信号曲线Fig.10 Signals of muscle tension during algospasm

3.4 病理性痉挛在时域特征表达

对受测者的下肢腓肠肌施加外部激励，在同一条件下获取病理性痉挛的特征，获取了相关测试数据，完成了数据的相关处理，包括标准差、绝对积分平均值(IAV)、均方根和方差等，实验测得数据如表4所示。

表4 3种不同类型痉挛时域特征的统计数据

由表4可得到：病理性痉挛的3类肌张力的平均值、信号采样值大小接近，不好区分。标准差、方差的值越小，其离散度就越低，实验数据的大小接近，则其稳定性好。另外，与阵挛、痛挛相比较，强直的标准差、方差小得多，相差近10倍，由此看出强直下的肌张力信号数据采集的样点具有较低的离散度。肌张力信号变化的强弱度主要是依据IAV、均方根所表示。对照痛挛、强直、阵性痉挛肌张力信号的3种统计量，均方根和IAV在数值上并不存在很明显的差异，不能完全适用于病理性痉挛的信号特征的提取要求。

在时域角度，病理性痉挛特征应用统计方法提取分析，简便快捷，对肌张力信号的时效性要求比较有利。仅用时域统计法，痉挛性肌张力信号在频率维度及其时间随机性的特征信息则难以提取，且存在着较大约束，不利于准确有效地判定阵挛和痛性痉挛。

4 下肢痉挛性肌张力的频域特征表达

针对病理性痉挛中阵挛和痛性痉挛信号不易量化与特征阈值随机出现的问题，在频域引入痉挛性肌张力信号，观测频谱和分析功率谱，揭示痉挛性肌张力信号的特征分量，进而在频域维度分解为一系列基本信号，剖析痉挛信号频率、幅度和相位的组合，提取反映病理性痉挛的特征阈值。

设获取的病理性肌张力信号为周期信号Fs(t)，将肌张力信号表示为傅里叶级数，基本形式如下：

(9)

式中，Ω1为角频率(Hz)，a0为直流分量，an、bn分别为幅值分量。

进一步，运用欧拉公式推出指数形式傅里叶级数为：

(10)

式(10)中，F(nΩ1)是傅里叶级数中指数形式的系数，导出可得：

(11)

结合病理性痉挛在临床上的表现，痉挛信号在展开处理后，根据其特点，确定周期信号采用阵挛信号，非周期信号采用强直和痛性痉挛，运用傅里叶数学方法分类处理上述信号，可以使信号在任意有限区间得到拓展或者是周期延拓，以获取各基波、谐波分量的频谱值，从而进一步获得痉挛性肌张力的特征信息。

4.1 运用B4-FFT算法处理痉挛性肌张力信息

运用康复机器人对脑卒中患者实施下肢步态训练是当前康复治疗的重要手段。为防范由于患者下肢在训练过程中随机出现的病理性痉挛及其带来的二次损伤，康复临床对病理性痉挛预警的实时性提出了很高的要求。一般地，降低痉挛传感器信号处理的运算量可有效提高传感器的响应速度，提升痉挛感知的实时性，进一步提升痉挛预警的安全效果。

针对痛挛，运用BFT算法，提取阵挛与痛挛的特征阈值，归纳病理性痉挛的信号规律及特点，算法表述如下：

(12)

已知信号序列x(n)，x(n)表示为实数或者复数，它的DFT变换为X(k)，若X(k)为复数，令DFT的变换核(加权因子)为

WN=e-j2π/N

(13)

相对一个实际痉挛信号，采集样本点数较多时，需要多次完成乘法运算，并且消耗大量的计算资源，严重影响了痉挛检测及预警的时效性。为提高运算速度，实时快速处理痉挛性肌张力数据，且DFT运算具有周期性特点，由式(13)分析可得出：

(14)

以N为周期，可以运用蝶算方式，具体如下：

(15)

对于N=Mr的DFT运算分r级，每级有N/M个蝶算，蝶算的复乘次数为：

(16)

与直接DFT计算相比，α为蝶算的复乘效率，若α越小则表明实时性越好；若N越大，则运算效率越高。

(17)

表5 3种算法的运算结果

通过表4比较得出，在上述3种算法中，DFT直接算法运算最快，但响应速度最慢，即痉挛信号表达的实效性较差，相比B2-FFT与B4-FFT二者的运算量，后者的运算效率比前者节省了50%，由于痉挛性肌张力信号的采样点数较多，随机性大、时效性要求高，故应用B4-FFT算法能分析其所具有的痉挛特征。

4.2 B4-FFT算法提取病理性痉挛的频域特征

结合上述实验所采集到的肌张力数据，根据病理性痉挛的阵挛和痛挛特征，运用B4-FFT算法提取在频域上的特征。

阵挛信号的频域处理结果如图11所示，阵性痉挛发生时，信号频率小于20 Hz，幅值超过40 000，频谱幅值对应的频率主要分布在6～10 Hz间，在Hz时对应的频谱幅值最大，相比于为16～18 Hz所对应的幅值变化更为显著，此时可判断阵挛发生。

图11 B4-FFT算法处理阵挛频谱曲线Fig.11 B4-FFT algorithm processing clonic spectrum

痛性痉挛信号的频域处理结果如图12所示，痛性痉挛发生时，信号频率小于20 Hz，幅值超过20 000；频谱幅值F(ω)对应的频率ω主要分布在2～8 Hz、11 Hz、13 Hz、19 Hz等范围内，在ω为2～8 Hz之间所对应的频谱幅值变化最大，频谱幅值F(ω)达到了50 000上下，相比于ω为11 Hz、13 Hz、19 Hz所对应的幅值则较为明显，频谱幅值F(ω)均在25 000左右，痛性痉挛呈现强直与阵挛信号特征的交替出现，并且在多个频率点上的幅值呈现小范围增高的现象，此时判断痛性痉挛发生。

图12 B4-FFT算法处理后的痛挛频谱曲线Fig.12 Spectrum of pain after B4-FFT algorithm processing

为了能够直观且量化地判定病理性痉挛的3类特征信息，归一化处理所获取的痉挛信号，归一化式子如下：

(18)

式中:P(i)为幅值阈值λ的所占百分比，幅值Fi(ω)为整个频率段的幅值，λ为该频率段的幅度阈值，λ设定区间为[20 000,100 000]。经过处理后，病理性痉挛的频域特征阈值指标如表6所示。

表6 阵挛与痛性痉挛频域特征阈值指标

综上所述，病理性痉挛信号经过B4-FFT算法处理后，特征阈值设置可总结为：当发生阵挛时，峰值频率分布较为单一，此时可设定多个特征阈值以判定阵挛是否发生；当发生痛性痉挛时，可设定幅值阈值λ来判定。

5 实验研究

实时获取痉挛信息是较为困难的，结合康复临床对病理性痉挛的特征评定，拟采用交互力激励人体下肢特定区域的肌张力，使该区域肌张力信号的特征表现，用于验证痉挛提取算法的有效性。如图13所示：选取健康标准的成年男性(64 kg/ 170 cm)为实验数据采集对象，再次实验采集外骨骼穿戴者的下肢腓肠肌的肌张力信号。在此基础上，获取肌张力信息及其特征阈值，验证痉挛信息提取算法的有效性。

实验设计为2种工况下的获取方式：静立时，手动激励下肢特定肌区以获取肌张力信号；外骨骼机器人步态激励下的肌张力信号。

图13 交互力激励下腓肠肌处肌张力信号提取图Fig.13 Signal extraction of muscle tension at gastrocnemius under interaction force excitation

5.1 交互力激励下的肌张力获取及特征评定

如图13(a)所示，受测者静止站立时，对下肢腓肠肌内侧施加手动激励，此时肌张力信号如图14所示。如图13(b)图13(c)所示，受测者穿戴外骨骼，外骨骼对人体下肢施加交互力激励，所获取的肌张力信息如图14所示。

如图14所示，静立姿态下，信号初始值表现为不平稳状态且肌张力幅值小于1 000 g，此刻放松下肢部分使其处于放松状态，经过手动激励后上升到大于2 500 g，达到了强直发生的阈值，伴随着非常明显的均匀波动，其值最大达到2 600 g以上，依据信号表现的特点和本文提出算法的评定，由此可以判定下肢肌肉出现强直。

图14 静止站立下的手动激励的肌张力信号曲线Fig.14 Signals of muscular tension induced by manual excitation while standing still

如图15所示，受测者穿戴下肢外骨骼康复训练机器人完成助力行走实验，腓肠内肌处于原始状态，此时肌张力信号值维持在1 000 g以内；随着外骨骼机器人助力行走，人体下肢在完成行走步态的同时受到外骨骼交互力激励，肌张力值在500～1 500 g之间，其波动状态处于稳定范围，最大值在3 000 g上下。可以看出：外骨骼施加激励下的肌张力信息已被有效提取。至于肌张力信息是否超过痉挛特征阈值以及评定痉挛类型，则需进一步通过B4-FFT算法处理，提取交互力激励下的肌张力信息中的痉挛信号段及特征信息，同时将病理性痉挛评定与特征阈值相结合。

图15 施加交互力激励下的肌张力信号变化曲线Fig.15 Variation of muscle tension signal under interaction force excitation

5.2 B4-FFT算法提取痉挛特征及实验分析

痉挛信号在临床上的获取具有随机性，很难得到真正的痉挛信号。根据临床痉挛信号的频率，我们通过人为手动的设置来获取相同频率信号。在交互力激励下的肌张力获取基础上，肌张力信息用BFT算法处理，验证其对痉挛特征提取的实际效果。用B4-FFT算法对以上2种实验过程中出现的肌张力信息展开处理，分别得到手动激励与人机交互力激励下的肌张力信息频谱图，如图16、图17所示。如图16所示：在2～20 Hz内有多个超过幅值阈值λ的峰值幅度；在2～10 Hz变化时，肌张力幅值相对比较集中，可据此判定是否发生痛挛；信号值处在16 Hz上下浮动时，其信号幅度存在单一峰值，可判定发生阵挛。综上所示，肌张力信息经B4-FFT算法处理后，可以有效提取病理性痉挛发生的特征阈值并给出评定。