一种自适应更新的在线运动单位识别方法技术

本发明专利技术公开了一种自适应更新的在线运动单位识别方法，包括：1、采集待测肌肉处的表面肌电信号，保存为离线数据集；2、采用自动化逐步变量剥离算法分解离线数据集的肌电信号，解算每个运动单位对应的解混向量并保存；3、在线处理阶段，实时获取表面肌电信号并对其分窗；4、求解每个分析窗内源信号的源信号；5、采用双线程并行操作的策略提取运动单位的发放序列并保存；6、随着实时数据的输入，执行步骤3到步骤5，待没有新的数据输入时，依次连接不同分析窗的分解结果并输出。本发明专利技术能够在线识别表面肌电信号中的运动单位活动信息，并随时间自适应更新以达到较高的精度和识别的运动单位数目。目。目。

【技术实现步骤摘要】
一种自适应更新的在线运动单位识别方法


[0001]本专利技术属于生物信号处理
，具体涉及一种自适应更新的在线运动单位识别方法，主要应用于智能鲁棒的肌电控制技术和人机接口。

技术介绍

[0002]表面肌电信号(Surface Electromyography,SEMG)是一种能够反映人体在进行肢体运动的过程中下发的神经驱动信息的电生理信号。实验所采集的肌电信号是由多个运动单位(Motor Unit,MU)受中枢调控激活发放产生的运动单位动作电位(Motor Unit Action Potential,MUAP)沿肌纤维传播在检测电极处的时间、空间叠加而成的。通过分析肌电信号的各项特征，能够直接解码用户的运动意图以及运动健康信息。此外，越来越多的学者发现，仅仅使用肌电信号的宏观时频域特征，对肌电控制技术性能的改善十分有限。因此，现有的肌电控制技术不再满足于各项宏观时频域特征的简单应用，转而聚焦于蕴含在神经驱动下由大量的信源产生的丰富活动信息，寻找更符合运动底层生理规律的控制策略。作为能直接获取到神经驱动信息的技术，肌电分解在近十年来被广泛应用于各项研究中。肌电分解可以看作肌电形成的逆过程，将混叠的表面肌电还原为其最基本的组成部分——MUAP序列的形式。MUAP序列中包含了运动单位的发放和波形信息，其中发放信息是最为重要的，直接展示了神经指令下发和传达的时刻，表征为运动单位的0
‑
1发放序列。肌电分解算法驱动着肌电控制技术朝着解码微观神经指令的方向逐步发展。
[0003]在过去的十年内，最具代表性的肌电分解算法包括卷积核抵消算法(Convolution Kernel Compensation,CKC)和逐步变量剥离算法(Progressive FastICA Peel
‑
off,PFP)。这两类算法都是在通识的盲源分离技术的基础上，结合肌电信号所特有的规律上所设计而成的。CKC算法粗略地估计运动单位发放序列和表面肌电信号之间的互相关矢量，然后通过估计的互相关矢量和EMG信号的相关矩阵来计算该运动单位的发放模式。PFP算法通过盲源分离经典算法快速独立成分分析(Fast Independent Component Analysis,FastICA)算法求解出一系列运动单位源信号，提取出其中包括的发放序列，通过带时域约束的FastICA对结果进行检测和修正，并通过剥离策略将已得到的运动单位从原始信号中剥离出来，以寻找到更多的运动单位。两种肌电分解策略在离线的过程中取得了较好的成果，且都有其对应的在线分解版本，为实际应用提供了众多有效实用的可参照方案。
[0004]然而，目前的以肌电分解技术为基础的在线运动单位识别方法存在诸多问题。通常来说在线的运动单位识别方法都是需要一段离线数据集获取分解所需要的解混向量，然后将其直接应用于在线分解的流程中，从而避免了庞杂的计算过程。但是，现有的在线运动单位识别方法在离线初始化解混向量后，在长时的在线数据处理过程中不再对其进行更新修正。这就导致随着采集时间的延长，用户肌肉状态和采集环境的逐渐变化导致解混向量的可用性逐渐降低，所获取的运动单位活动信息也就不再可靠。此外，在用户长时运动过程中，存在新的运动单位募集现象，而这些新募集的运动单位是无法被离线初始化过程所捕获，从而消失于在线分解的结果中。如何在实时识别运动单位的同时更新解混向量集以确
保算法框架的高性能，是一个巨大的困难和挑战。

技术实现思路

[0005]本专利技术为了克服现有在线运动单位识别方法的不足之处，提出一种自适应更新的在线运动单位在线识别方法，以期能够在线识别表面肌电信号中的运动单位活动的同时自适应更新解混向量集，保持较高的精度用以持续追踪相同的运动单位且识别出长时肌肉活动过程中新募集的运动单位，从而达到与离线分解相当的效果。
[0006]本专利技术为解决技术问题，采用如下技术方案：
[0007]本专利技术一种自适应更新的在线运动单位识别方法的特点是按如下步骤进行：
[0008]步骤1、通过肌电采集设备和N通道高密度电极在待测肌肉处采集t时刻肌肉收缩任务的表面肌电信号数据，记为x(t)＝[x1(t),x2(t),
…
,x
I
(t),
…
,x
N
(t)]T
，其中，x
I
(t)表示t时刻第I个通道的肌电数据；从而得到时长为T1的表面肌电信号数据并作为离线数据集；
[0009]步骤2、采用自动化逐步变量剥离算法离线分解所述离线数据集中的表面肌电信号：
[0010]步骤2.1、利用式(1)对所述离线数据集进行扩展，得到扩展信号
[0011][0012]式(1)中，K是信号扩展的延迟因子；
[0013]定义并初始化运动单位发放集合γ及其对应的解混向量集合ε均为空集；
[0014]定义并初始化残差信号
[0015]步骤2.2、定义并初始化剥离次数peel_num＝1；
[0016]步骤2.3、执行盲源分离的FastICA算法，并设置算法执行的总次数为n：
[0017]步骤2.3.1、定义并初始化变量i＝1；
[0018]步骤2.3.2、定义并初始化变量j＝0；
[0019]步骤2.3.3、定义并初始化残差信号中第i个解混向量在第j次的计算结果w
i,j
为一个均值为0，方差为1的高斯信号；
[0020]步骤2.3.4、利用式(2)计算第i个解混向量在第j+1次的计算结果w
i,j+1
：
[0021][0022]式(2)中，G'表示非多项式函数G的一阶导数；G”表示G的二阶导数；E是数学期望；代表第i个解混向量在第j次的计算结果w
i,j
的转置；
[0023]步骤2.3.5、如果则令j+1赋值给j，返回至步骤2.3.3；其中，表示预设的迭代计算次数上限；
[0024]若i<n且则表示得到第i个解混向量w
i
＝w
i,j
，并令i+1赋值给i，返回至步骤2.3.2；
[0025]否则，表示得到n个解混向量所组成的解混矩阵W＝[w1,w2,
…
,w
i
,
…
,w
n
]T
，并利用式(3)得到正交化后的解混矩阵W'：
[0026]W'＝(WW
T
)
‑
1/2
W
ꢀꢀ
(3)
[0027]步骤2.3.6、计算运动单位源信号且Y＝[y1,y2,
…
,y
i
,
…
,y
n
]T
，其中，y
i
代表第i个解混向量w
i
所对应的运动单位源信号；
[0028]步骤2.4：采用图像处理的大津法获取第i个运动单位源信号y本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种自适应更新的在线运动单位识别方法，其特征是按如下步骤进行：步骤1、通过肌电采集设备和N通道高密度电极在待测肌肉处采集t时刻肌肉收缩任务的表面肌电信号数据，记为x(t)＝[x1(t)，x2(t)，...，x
I
(t)，...，x
N
(t)]
T
，其中，x
I
(t)表示t时刻第I个通道的肌电数据；从而得到时长为T1的表面肌电信号数据并作为离线数据集；步骤2、采用自动化逐步变量剥离算法离线分解所述离线数据集中的表面肌电信号：步骤2.1、利用式(1)对所述离线数据集进行扩展，得到扩展信号步骤2.1、利用式(1)对所述离线数据集进行扩展，得到扩展信号式(1)中，K是信号扩展的延迟因子；定义并初始化运动单位发放集合γ及其对应的解混向量集合ε均为空集；定义并初始化残差信号步骤2.2、定义并初始化剥离次数peel_num＝1；步骤2.3、执行盲源分离的FastICA算法，并设置算法执行的总次数为n：步骤2.3.1、定义并初始化变量i＝1；步骤2.3.2、定义并初始化变量j＝0；步骤2.3.3、定义并初始化残差信号中第i个解混向量在第j次的计算结果w
i，j
为一个均值为0，方差为1的高斯信号；步骤2.3.4、利用式(2)计算第i个解混向量在第j+1次的计算结果w
i，j+1
：式(2)中，G
′
表示非多项式函数G的一阶导数；G
″
表示G的二阶导数；E是数学期望；代表第i个解混向量在第j次的计算结果w
i，j
的转置；步骤2.3.5、如果则令j+1赋值给j，返回至步骤2.3.3；其中，表示预设的迭代计算次数上限；若i＜n且则表示得到第i个解混向量w
i
＝w
i，j
，并令i+1赋值给i，返回至步骤2.3.2；否则，表示得到n个解混向量所组成的解混矩阵W＝[w1，w2，...，w
i
，...，w
n
]
T
，并利用式(3)得到正交化后的解混矩阵W
′
：W
′
＝(WW
T
)
‑
1/2
W
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)步骤2.3.6、计算运动单位源信号且Y＝[y1，y2，...，y
i
，...，y
n
]
T
，其中，y
i
代表第i个解混向量w
i
所对应的运动单位源信号；步骤2.4：采用图像处理的大津法获取第i个运动单位源信号y
i
的阈值th
i
：步骤2.4.1、将区间(0，max(y
i
))划分为M个等分，得到M个等分值[th
i，1
，th
i，2
，th
i，3
，...，th
i，m
，...，th
i，M
]，其中，max(y
i
)表示运动单位源信号y
i
的最大值，th
i，m
表示第m个小区间的最大值；步骤2.4.2、定义并初始化变量h＝1；步骤2.4.3、令th
i
＝th
i，h
；求出y
i
中大于th
i
的采样点个数N
0，h
和小于th
i
的采样点个数N
1，h
；并求解个数N
0，h
对应的灰度μ
0，h
作为y
i
中所有大于th
i
的采样点的平均值；求解个数N
1，h
对应的灰度μ
1，h
作为y
i
中所有小于th
i
的采样点的平均值；步骤2.4.4、根据式(4)得到第h次计算的类间方差F
i，h
：
F
i，h
＝N
0，h
×
N
1，h
×
(μ
0，h
‑
μ
1，h
)2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)步骤2.4.5、若h＜M，则令h+1赋值给h，返回至步骤2.4.3；否则，进入步骤2.4.6；步骤2.4.6、找到M个类间方差中的最大值所对应的下标h
max
，并将下标h
max
所对应的作为第i个运动单位源信号y
i
的阈值th
i
；步骤2.5、将y
i
中大于阈值th
i
的波峰峰值位置置为“1”，其位置余置为“0”，从而得到第i个解混向量对应的运动单位发放序列r
i
；步骤2.6、通过盲源分离算法中带时域约束的FastICA算法检测分解结果的可靠性：步骤2.6.1、定义迭代计算的次数为d，并初始化d＝1；定义并初始化第d次迭代的拉格朗日乘子μ
d
＝0、第d次迭代的惩罚因子δ
d
为随机小数；初始化第d次迭代的第i个运动单位源信号y
i，d
＝y
i
、第d次迭代的第i个解混向量w
i，d
＝w
i
；定义第i个运动单位源信号y
i
和待测的运动单位发放序列r
i
之间的第d次迭代的相关性函数为其中，ξ
d
是第d次迭代的相关性下界，并初始化ξ
d
＝1；步骤2.6.2、通过式(5)得到第d+1次迭代的第i个解混向量w
i，d+1
：式(5)中，g
′
d
为第d次迭代的相关性函数g
d
的一阶导数；步骤2.6.3、对第d+1次迭代的第i个解混向量w
i，d+1
进行标准化：将w
i，d+1
/||w
i，d+1
||2赋予w
i，d+1
，其中||w
i，d+1
||2是w
i，d+1
的2
‑
范数；步骤2.6.4、根据式(6)、式(7)和式(8)分别更新其他参数：μ
d+1
＝max{0，μ
d
+δ
d
g
d
}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)δ
d+1
＝α
×
δ
d
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)ξ
d+1
＝β
×
ξ
d
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)式(6)、式(7)和式(8)中，μ
d+1
表示第d+1次迭代的拉格朗日乘子，δ
d+1
表示第d+1次迭代的惩罚因子，ξ
d...

【专利技术属性】
技术研发人员：张旭，赵昊文，宋芮，李心慧，陈香，
申请(专利权)人：中国科学技术大学，
类型：发明
国别省市：