1. 引言

脑性瘫痪(Cerebral Palsy, CP)，简称脑瘫，是自受孕开始至婴儿期由于非进行性脑损伤和发育缺陷所导致的永久性中枢神经障碍综合征[1]，通常伴有神经肌肉系统功能异常[2]，如肌肉活动异常或协同异常引起运动障碍。考虑到，表面肌电(surface electromyography, sEMG)信号，携带有丰富的神经肌肉控制信息和肌肉激活信息，能反映神经肌肉运动状态。因此，sEMG信号已被广泛应用于脑瘫患者的步态分析、日常活动监测等研究中。研究者Bernhard通过下肢胫骨前肌、比目鱼肌和腓肠肌的sEMG信号来评估19名CP患儿步态过程中的踝关节运动[3]。Josephine的研究表明10名CP患儿下肢股四头肌、腓肠肌、胫骨前肌的最大随意收缩在行走、跑步、跳跃等运动任务中，与正常发育的同龄人相比显示出显著较低的EMG水平[4]。可见，sEMG能够提供实时的肌肉活动信息，帮助研究者评估运动能力、监测康复过程。然而，目前的研究CP患儿数目较少，需进一步扩充样本数量。由于CP患儿sEMG信号采集受患者数及伦理审批的限制，特别是面临不同运动障碍分级的患者数的短缺，快速获取大规模样本量用于CP患儿研究仍面临挑战。因此，能准确生成sEMG信号及其特征信号是解决患者样本量不足的核心。

目前关于肌电信号的生成研究主要有以下两种方法。第一种，从生理角度出发。早期的一些研究主要基于sEMG信号产生的原理，通过运动单位数量和运动单位动作电位的特定叠加合成sEMG信号。2001年，Dario等人构建了融入肌肉、脂肪和皮肤层参数的体积导体模型以实现对EMG信号的生成和检测[5]。进一步，D. Farina等人引入多层圆柱体描述体积导体建立sEMG信号生成模型[6]。第二种，基于对抗神经网络生成EMG信号[7]。起初，Fei Zhu等人提出了一种将双向长短期记忆网络(Bi-LSTM)与卷积神经网络(CNN)相结合的生成对抗网络(GAN)模型，用于心电图数据的生成[8]。随后，2020年，Hazra等人提出了一种新颖的SynSigGAN模型能够自动生成与实际EMG信号极为相似的时序数据[9]，钱尧佳提出了一种基于时间加权生成对抗网络(TWGAN-GP)的sEMG特征样本生成模型[10]。上述研究已完成对正常人的生理电信号的生成，然而还缺少针对患者sEMG信号的生成。为深入探索CP患儿肌肉活动异常或协同异常引起的运动障碍，研究者们常常分析sEMG信号的时频域特征[11]。Campanini等人使用sEMG信号的包络信号来评估步态中电极位置不同所带来的影响[12]。为了识别复杂任务中肌肉激活的准确性，Francesco通过使用sEMG信号的线性包络峰值来确定肌肉的激活次数[13]，针对CP患儿，李飞基于sEMG信号的包络信号提取CP患儿的肌肉协同模式，用于评估患者的运动障碍程度并揭示了患者病理机制[14]。可见，肌电包络信号作为sEMG信号的关键特征之一，能够直观地反映肌肉的激活时间及激活强度，因而成为研究肌肉功能及其障碍的重要指标。

基于此，本文以肌肉协同模型为基础，提出了一种新型的方法用于CP患儿肌电包络信号的生成，旨在解决CP患儿样本不足的问题。具体而言，采集健康成人和CP患儿步态过程双下肢各8块肌肉的sEMG信号，通过非负矩阵算法提取受试者的协同结构矩阵和激活系数矩阵，以健康成人的激活系数矩阵为参考，通过引入时间延迟变量获得一系列更新后的激活系数矩阵，再与CP患儿的协同结构矩阵相乘，进而得到全新的肌电包络信号。最终，通过对比生成信号与真实CP患儿步态中各肌肉激活时间，进一步验证本研究提出的包络信号生成模型的有效性。

2. 信号采集

2.1. 实验对象

本研究信号采集的方案包括4名健康成人(AD组，年龄23.0 ± 2岁)，6名脑瘫患儿(CP组，年龄2.32 ± 1岁)共10名受试者。AD组实验对象来自实验室同学，CP组实验对象招募于上海养志康复医院。表1和表2分别展示了健康成人和脑瘫患儿的相关信息。根据研究的需要，对脑瘫患儿的选取标准为：1) 临床医生已确诊其为脑瘫患儿；2) 具有明显的步态异常症状，但能够独立或在他人辅助下完成迈步行走动作；3) 除先天性脑瘫外没有患有其他对步态运动有障碍的疾病。GMFCS级别由三位临床医生进行评估，并对三位医生所评定的级别结果取平均。本研究采用的实验方案经评估后获得康复大学的伦理审查委员会的批准(KFDX: NO. 2021-035)。

Table 1. AD group basic information

表1. AD组基本信息

Table 2. CP group basic information

表2. CP组基本信息

2.2. 实验方案

步态运动是指由中枢神经系统控制下肢肌肉驱动髋、膝、踝、足趾的一系列连续活动，激活的肌肉多且复杂。考虑到脑瘫患儿的肌肉相较于成人较小，因此选取主要参与步态运动的8块尺寸偏大的肌肉作为目标肌肉。分别为：胫骨前肌(TA)、比目鱼肌(SOL)、腓肠肌(LG)、股直肌(RF)、股外侧肌(Vlat)、半腱肌(ST)、股二头肌(BF)和阔筋膜张肌(TFL)。传感器放置位置如图1所示，其中，蓝色标记为sEMG信号传感器，红色标记为加速度信号(acceleration signal, acc)传感器。实验中受试者被要求以合适的速度步行，左右腿各完成约30个迈步动作。信号采集设备为实验室自行研发的多通道表面肌电采集系统，最多可同步采集16个通道的sEMG信号和2个三轴acc信号，表面肌电电极的采样频率为1000 Hz，acc信号传感器采样频率为200 Hz。

Figure 1. Schematic diagram of the placement of surface electromyography electrodes and acceleration sensors

图1. 表面肌电电极与加速度传感器放置位置示意图

3. 数据处理与分析

Figure 2. Signal processing flow chart

图2. 信号处理流程图

本文对实验采集10名受试者双下肢各8块肌肉的sEMG信号和acc信号进行了信号处理，具体整体流程如图2。

3.1. 信号预处理

3.1.1. sEMG信号预处理

首先需要对信号进行去噪处理。由于人体sEMG信号的能量主要集中在20~500 Hz，而运动伪迹的噪声频率通常集中在0~20 Hz范围内，故可通过低截止频率为20 Hz、高截止频率为500 Hz的巴特沃斯6阶带通滤波器滤除噪声。同时，采用50 Hz的陷波滤波器来消除工频干扰。最后，每通道的sEMG信号通过截止频率为5 Hz，阶数为1的低通滤波器获得sEMG包络信号。以健康成人受试者AD1的比目鱼肌为例，预处理后的sEMG如图3(c)所示：

Figure 3. Preprocessing of Soleus muscle sEMG signals and acc signals in AD1

图3. AD1比目鱼肌sEMG信号和acc信号的预处理过程

3.1.2. acc信号预处理

一个完整的步态周期可以定义为第一次脚跟着地到第二次脚跟着地。脚跟着地时小腿动作的迅速变化引起加速度变化，加速度传感器竖直方向加速度(y轴acc)达到峰值，因此可通过检测加速度峰值点来分割步态周期[13]。将检测到的加速度信号峰值点映射回sEMG来实现步态周期的分割，加速度处理结果如上图3(d)所示。

3.2. 肌肉协同提取

肌肉协同的概念最早由Bernstein提出并由D’avella等人完善，指的是人体任何运动都是由多块肌肉共同协同配合才能完成，并不是任何一单块肌肉作用[15]。通过非负矩阵分解算法(Non-Negative Matrix Factorization, NMF)可将sEMG包络信号( $V_{m\times t}$ )分解为协同结构矩阵( $W_{m\times n}$ )与激活系数矩阵( $H_{n\times t}$ )，如公式(1)。协同结构矩阵 $W_{m\times n}$ 反映不同肌肉之间的协同关系，描述每个协同模式中各肌肉的相对贡献；而激活系数矩阵 $H_{n\times t}$ 则表示协同结构模式在时间上的激活强度，展现肌肉在整个时间序列内的动态激活情况。该分解过程可以有效地描述运动过程中肌肉群的收缩行为，便于对运动过程中产生的肌电信号进行分析。

$V_{m\times t}=W_{m\times n}\times H_{n\times t}\times Error$ (1)

${V^{\prime }}_{m\times t}=W_{m\times n}\times H_{n\times t}$ (2)

式中， $m$ 为所测肌肉数目； $t$ 为样本的时间序列长度，即采样数； $n$ 为肌肉协同数目；Error为分解误差。 ${V^{\prime }}_{m\times t}$ 表示 $W_{m\times n}$ 和 $H_{n\times t}$ 相乘得到的重构肌电包络矩阵。

提取肌肉协同的过程中，确定肌肉协同数目的方法最常用的就是通过变异比例系数(Variability Accounted For, VAF)来衡量。为了确定最小的肌肉协同数目，通常计算每一个 $n$ 值下由分解得到的 $W_{}$ 和 $H$ 矩阵重构的肌电包络信号 ${V^{\prime }}_{}$ ，然后根据 $V^{}{}_{}$ 与 ${V^{\prime }}_{}$ 的数据变异度(Variability accounted for, VAF)来确定肌肉协同数目[16]，其公式如下：

$VAF=1-\frac{V-V^{\prime }}{V^2}$ (3)

式中，VAF的取值范围在[0, 1]，VAF的值越大则重构矩阵与原始矩阵之间的误差越小(即 $V$ 和 $V^{\prime }$ 越接近)。从肌肉协同数目 $n=1$ 开始，依次增加进行矩阵分解，直到当前选定的肌肉协同数目使得VAF ≥ 95%则停止计算[17]，此时的 $n$ 即为选定的肌肉协同数目。

1) CP患儿肌肉协同提取

对6名CP患儿步态分割后的包络信号通过公式(1)可分别提取出2、3和4个肌肉协同模式。如图4所示，不同患儿的协同结构、激活系数以及协同数目存在差异。GMFCS I级患者为CP5，其协同数目为4个；II级患者为CP2，CP6其协同数目分别有2个和3个；III级患者为CP1，CP3，CP4其协同数目有3个。可见较少的协同数目倾向于出现在运动功能评级较低的患者中，即在运动功能严重受限的患儿中更常见。本文与Tang的研究[14]取得了一致的结果。

Figure 4. Muscle synergy extracted from 6 children with cerebral palsy

图4. 对6名脑瘫患儿提取的肌肉协同

研究发现，协同数目较少的患儿通常表现出更显著的肌肉共激活模式，这表明在运动功能受损的情况下，患儿的肌肉可能更难以独立控制，从而倾向于同时激活多块肌肉来完成运动任务。这一现象揭示了协同数目与运动功能评级之间的相关性，表明在部分脑瘫患儿中，运动功能越严重，肌肉间的协同关系越趋向于简单化，并体现出肌肉协同数量的减少和共激活现象的增加[18]。

2) 健康成人肌肉协同提取

将健康成人步态分割后的包络信号通过公式(1)可提取肌肉协同。值得注意是健康成人协同数目的确定。考虑到，肌肉协同募集方式的变化会导致肌肉激活与组合模式的异常，进而影响正常运动子任务发生的时刻与强度。因此，本文将脑瘫患儿肌肉协同数目映射到健康成人，以获得正常下行神经元募集步态子任务发生的时刻与强度。因此，健康成人协同数目 $n_{\begin{array}{c}TD\\ \end{array}}$ 是由脑瘫患儿的协同数目确定，即 $n_{\begin{array}{c}TD\\ \end{array}}$ 为2~4。以脑瘫患儿协同数目 $n_{\begin{array}{c}CP\\ \end{array}}$ = 2为例，为获得相应协同数目下的激活系数，公式(1)中协同数目设定为2，即 $n_{\begin{array}{c}TD\\ \end{array}}$ = 2，将健康成人步态分割后的包络信号通过NMF可提取数目为2的 $W_{m\times n}$ 及相应的 $H_{n\times t}$ 。为获得协同数目为2~4的健康成人协同模式模板，将4名健康成人的 $W_{m\times n}$ 及 $H_{n\times t}$ 取平均为 $\overline{W_{m\times n}}$ 及 $\overline{H_{n\times t}}$ (字母上面加横线，表示平均)。图5展示了 $n_{\begin{array}{c}TD\\ \end{array}}$ 。分别设置为2，3，4时，得到的协同结构即募集系数的平均模板，分别定义为 $\overline{H_{n=2}}，\overline{H_{n=3}}，\overline{H_{n=4}}$ 。

Figure 5. Average muscle synergy template for healthy adults

图5. 健康成人的平均协同模板

3.3. CP患儿包络信号生成

为生成符合CP患儿步态特征的肌电包络信号，本研究通过公式(2)将 $W_{}$ 矩阵与 $H$ 相乘。其中， $W_{}$ 矩阵使用脑瘫患儿的协同结构矩阵以保留其特有的肌肉协同关系，而 $\overline{H_{}}$矩阵则为健康成人的激活系数矩阵均值作为模板，以保证下行运动神经信号对肌肉在整个步态周期内的正常动态激活。同时，通过对健康成人的 $H$ 矩阵施以不同程度的时间延迟来实现动态调整，旨在模拟CP患儿步态中的协同结构被激活的时间序列。

为了更为细化激活时序的变化，分别以整个步态周期的10%、5%、1%为步进参数进行动态调整。健康成人步态中的肌肉激活具有较强的时序规律性，而脑瘫患儿的神经肌肉控制受损，导致其步态中激活时序存在差异。通过对健康成人 $H$ 矩阵施以不同程度的时间延迟，从而更贴近脑瘫患儿步态中的实际肌肉激活特点。

4. CP患儿包络信号生成结果与分析

4.1. CP患儿包络信号生成结果

在完成对6名TD成人以及4名CP患儿的肌肉协同提取之后，分别通过将健康成人动态时间调整后的 $\overline{H_{}}$矩阵与CP患儿的 $W_{}$ 矩阵进行矩阵相乘，以生成CP患儿sEMG包络信号。以CP1患儿为例，采用1%、5%、10%时间延迟其生成结果如图6所示，其中图(a)为相乘后的8块肌肉sEMG包络信号，图(b)为各肌肉sEMG包络信号的激活部分。

Figure 6. The envelope curve and activation time diagram of each muscle after CP1 time shift

图6. CP1时移后的各肌肉包络曲线及激活时间图

4.2. 验证实验

4.2.1. CP患儿的肌肉激活时间

Figure 7. Lower limb muscle activation time during gait in children with CP

图7. CP患儿步态中下肢肌肉激活时间

为进一步验证生成信号的准确性，本文通过对比真实CP患儿步态周期中8块肌肉激活时间的占比进行确认。通过对痉挛性CP患儿步态下肢肌肉激活程度的相关多篇文献进行归纳总结[3] [4] [15] [19]-[22]，获取了实际8块肌肉在整个步态周期中的激活时序，如图7所示。对比健康成人8块肌肉激活时序图[19]，痉挛型CP患儿的步态过程肌肉激活时间相对正常人出现提前或延迟，且激活持续时间更长。研究者表示上述表征可能与脑瘫患儿运动控制困难和能量消耗增加有关[21]。根据步态周期的划分标准，本文同样将一个完整的步态周期划分为7个时相[19]，分别为承重期(10%)，支撑中期(30%)，支撑末期(50%)，预摆动期(60%)，摆动前期(73%)，摆动中期(87%)，摆动末期(100%)。

4.2.2. 验证参数

为评估生成信号的准确性，确保生成的sEMG包络信号能够精确反映CP患儿肌肉激活模式，本文通过生成信号的肌肉激活与真实CP患儿步态中下肢肌肉激活的重合度进行衡量。交并比(Intersection over Union, IoU) [23]为重合度的度量指标，如公式(5)。IoU的取值范围为0至100%，数值越大表示集合的重叠程度越高，通常IoU值大于50%即被认为具有较高的相似性。

4.2.3. 验证结果

采用1%、5%、10%时间延迟的 $\overline{H_{}}$矩阵与CP患儿的 $W$ 矩阵进行矩阵相乘，生成sEMG包络信号的肌肉激活与真实CP患儿步态中下肢肌肉激活的重合度进行计算，以获得生成sEMG包络信号准确性。以患者CP1为例，将健康成人 $\overline{H_{n=3}}$ 矩阵乘以 $W_{CP1}$ 获得的生成信号。应用到所有受试者，计算每块肌肉的IoU均值，以全面评估生成信号的准确性，结果如表3所示。

Table 3. The average IoU of the signals generated by all children

表3. 所有患儿生成信号的IoU均值

本研究通过IoU评估生成信号和真实CP患儿肌肉激活时间之间的相似性，其平均相似性为64.94%，表明生成信号在总体上能够较好地反映CP患儿的肌肉激活特征，同时在某些肌肉的激活时间匹配上存在一定差异。其中ST的相似性为71.83%，是本次实验中表现最佳的肌肉。ST在步态支撑期中负责膝关节的弯曲，其较高的相似性表明该生成信号能够准确捕捉该肌肉的激活模式。同样地，TFL在步态中起稳定髋关节的作用，而RF在摆动期负责膝关节的伸展[19]，二者的相似性分别为69.39%和69.08%，其较高的相似性可能与这些肌肉在步态中的激活模式相对稳定有关。此外，LG、Vlat和SOL的相似性分别为64.35%、62.06%和60.69%略低于前述肌肉。LG在步态推进期中提供主要动力，Vlat在支撑期中负责膝关节的稳定，而SOL则在步态支撑期中起稳定踝关节的作用。然而TA的相似性为52.84%，由于TA在步态摆动期中主要控制足部下落，负责足背屈，其较低的相似性可能与CP患儿在该肌肉激活模式上存在较大的病理变异有关，这种现象可能源于以下两个方面：一是上行运动神经元损伤导致的选择性运动控制能力下降；二是拮抗肌(跖屈肌群)的痉挛抑制了TA的正常功能。

以上结果证实了即使在病理状态下，运动控制依旧遵循一定程度的模块化组织原则，未来的研究可以针对相似性较低的肌肉进行模型优化，以提高生成信号在所有肌肉上的匹配度。

5. 讨论与分析

目前，针对CP患儿sEMG包络信号的生成研究尚少，现有方法多集中于健康人或通用病理模型的信号生成，难以直接应用于CP患儿的sEMG包络信号生成，本研究首次提出了一种基于时变肌肉协同模型的CP患儿sEMG包络信号生成方法，在解决CP患儿样本不足问题上展现了独特的优势。与现有sEMG信号生成方法相比，基于深度学习的信号处理方法虽能实现信号重构，但难以捕捉肌肉协同的动态特性，且由于其大规模数据训练的需求限制了其在小样本场景中的应用。本文方法通过结合时变肌肉协同模型和非负矩阵分解，生成的信号具有明确的生理意义，通过协同结构矩阵和激活系数矩阵，揭示了CP患儿肌肉激活的肌肉空间协同异常和时间动态病理特征。且在小样本条件下实现了64.94%的平均IoU相似性，显著优于传统方法和部分深度学习方法，为临床数据稀缺场景提供了可行解决方案。

然而，本研究仍存在一定的局限性，生成信号的准确性在复杂病理特征肌肉(如胫骨前肌)上表现较差(IoU = 52.84%)，可能原因在于模型对CP患儿某些肌肉高度变异激活模式的适应性不足，如上行运动神经元损伤导致的选择性控制能力下降，未来可结合深度学习的时序建模能力(如LSTM或Transformer)和多模态数据(如关节运动学参数)进行改进，进一步提升模型对非线性肌肉协同关系的捕捉能力，可以融合多模态数据来同时增强生成信号的时空一致性。本研究首次验证了基于肌肉协同理论的CP患儿sEMG生成框架的可行性，揭示了肌肉协同模型在病理信号生成中的潜力，为病理运动分析提供了新的思路。

6. 结论

本文为解决CP患儿样本不足的问题，提出了一种以时变肌肉协同模型为基础的CP患儿肌电包络信号生成的方法。通过采集健康成人和CP患儿步态过程双下肢各8块肌肉的sEMG信号和加速度信号，利用非负矩阵分解算法提取受试者的协同结构矩阵和激活系数矩阵，并以健康成人的激活系数矩阵为参考，通过引入时间延迟变量获得一系列更新后的激活系数矩阵，随后与CP患儿的协同结构矩阵相乘，进而生成符合CP患儿病理特征的sEMG包络信号。为验证生成信号的准确性，通过交并比IoU评估生成信号和真实CP患儿肌肉激活时间之间的相似性，实验结果表明，其平均相似性为64.94%，其中某些关键肌肉(如RF、ST)的IoU值超过69%，验证了模型的有效性，表明生成的包络信号符合CP患儿的肌肉激活特征。本研究的主要贡献在于提供了一种可解释的小样本CP患儿sEMG包络信号生成框架，突破了传统方法对大批量数据或通用病理模型的依赖，为CP患儿运动康复研究提供了可解释性的数据支持，同时揭示了肌肉协同模型在病理信号生成中的潜力。

NOTES

^*通讯作者。

参考文献