本文提出一种不需人工干预的EMG信号自动分解算法。算法的学习阶段使用初始的一小段EMG记录以估计不同运动单位的动作电位模板。分解阶段根据运动单位发放统计特性,采有基于最小错误率的贝叶斯准则,将EMG信号分解为构成它的各运动单位动作电位序列,这里重点讨论了具有迭加动作电位波形的EMG信号分解问题。最后利用合成的模拟EMG信号及真实EMG信号对算法进行了检验,表明其具有很高的分解正确率.
分类号: R318.04; R540.41
STUDY ON AUTOMATIC DECOMPOSITION OF EMG SIGNALS
WITH OVERLAPPING MUAP WAVEFORM
Yang Jihai, Yang Hongning, Zhou Ping, Zhang jinsong, Sun Luyao
(University of Science and Technology of China, Hefei, Anhui, 230026)
ABSTRACTA scheme of EMG signal decomposition was reported no need of human supervision. A learning stage was first applied at the beginning of recording to estimate typical spike shapes of different motor units (MUs). As for the decomposition stage, a method was devoloped, with special consideration of temporal spike overlap. It minimized the probability of error, taking into accout the statistical properties of discharges of MUs. The method was tested on a real recording as well as synthetic data.
Key words: Automatic decomposition; MUAP waveform overlapping
0前言
肌肉随意收缩过程中,运动单位(MU)在中枢神经系统的控制下产生兴奋脉冲,即运动单位动作电位(MUAP),运动单位动作电位序列(MUAPT)在空间和时间的叠加结果就形成了肌电信号(EMG)。而EMG信号的分解则是其形成的求逆过程,即根据EMG信号形成的机理,将其还原为构成它的各运动单位动作电位序列MUAPT的过程[1]。EMG信号分解的研究无论是对神经肌肉系统控制特性的基础研究,还是对临床诊断的应用研究都具有十分重要的意义。
1分解方法
1.1MUAP的自学习
EMG信号首先被分成许多活动段(ACS)和非活动段[2],活动段包含有激活的MUAP信息,而非活动段仅由噪声和背景肌电活动组成。区分的方法是采用移动窗峰-峰值检测:取1ms长度的移动窗,若窗内EMG信号的峰-峰值超过某一阈值Vt1,判断为一个含有MUAP的ACS的开始;而当峰-峰值下降到另一阈值Vt2时,则判断该ACS结束。
MUAP的自学习过程包括对MU个数和模板的确定、各MU平均发放率(FR)的估计。采用相关法[3]对MUAP进行分类,以建立各MUAP模板。设当前检测到的含有MUAP的ACS波形序列为xi(i=1,2,…,N),已有一组MUAP模板集合为{S1,…,Sj,…SM},模板Sj中的MUAP序列为Sji(i=1,2,…,N),令它们的均值分别为:
定义信号xi与Sji的相关系数为:
Cj越大表明波形xi与模板Sji相关性越好,不过由于相关系数对波值不敏感,实际判决时还增加一个条件,即只有同时满足Cj≥0.95及波峰值差≤36.6μV时,才把xi归入Sji中。如果xi与所有的模板均不能满足上述两个条件,则把xi作为一个新的模板存入模板集合。
处理中力求不遗漏MU并准确地建立MUAP模板,采取的方法是适当提高相关系数Cj的判决门限。不过这种办法可能导致属于同一MU的MUAP被划归为几个类,主要是由于受到较多干扰的那些MUAP引起的,由于干扰的随机性,这类MU所包含的MUAP数一般很少。为此,我们在处理中把所含MUAP个数小于3的那些模板去除,这样做的结果既可以把各个MU中有所变异的MUAP组成的模板去掉,也可把多个MUAP叠加波形所形成的模板删除。后者所以能删除是因为MU发放的随机性,在较短的学习EMG信号中,多个MUAP重叠波形出现概率较低的原因。
通过上述处理得到的各模板内的MUAP具有高度相似性,从而由它们经迭加平均产生的模板Sj也比较准确。
建立模板采用遗忘因子算法[4]:
Sj(new)=[(M-1)*Sj(old)+X](2)
其中,M为已检测出的属于该MU的MUAP数,Sj模板序列,X为信号序列。
计算结果模板集合中包含的模板个数即是募集到的MU个数。
为了估计平均FR,仍采用相关法求出ACS信号序列X与属于MUi的模板Si的相关系数Ci,i=1,2,…,N。若存在Ci,它满足所有Ci≥0.85(1<i<N),波峰值差≤36.6μV,并且是所有Ci中的最大者,则判断X属于MUi;若不存在满足上述两个条件的Ci,判断X是噪声,拒绝分类。考虑到尽可能精确地得到各MU的平均FR,希望分类中尽量不遗漏任何一个MUAP,阈值就取得略小一些(0.85)。学习EMG信号采用一秒时间长初始段的EMG信号,当它们处理完后,各MU所包含的MUAP个数就作为各自的平均FR估计。
1.2EMG信号的分解
对EMG信号的分解实际就是将上述检测到的ACS进行判断并分类,确定它是否含有MUAP,含有几个MUAP,各MUAP分别属于哪个MU。分解步骤是,对于各ACS,首先采用相关法,将它们与各模板Si进行相关比较,如果相关系数C≥0.97且波峰值差≤36.6μV,则判断ACS只含有属于该MUi的一个MUAP。若不满足,则继续下一步,用更准确的基于Bayes决策理论的最大后验概率判别准则进行分解,这里包括对具有叠加MUAP波形的分解问题,其分解算法如下:
设X为检测到的ACS波形的特征矢量(本文采用波形的采样点),ACS的起始时间为T1,结束时间为T2,间隔I=(T1,T2)。为简化起见,这里仅考虑有两个MUAP叠加的波形。定义:
Pk=P(I中有MUk的发放|X,A)k=1,2,……,N
Pij=P(I中有MUi与MUj的发放|X,A)i,j=1,2,……,N, j<i
其中A是间隔I中有一个或两个MUAP出现的概率事件,Pk是I中属于MUk的MUAP出现的概率,Pij是属于MUi和MUj两类MUAP迭加形成的波形出现的概率。
令Hi(t)表示MUi在时刻t产生了一个MUAP,f(Hi(t))表示其概率密度函数。这样可将分类问题转化为对一些似然函数的比较。
对于Pij,可以表示为:
利用Bayes准则,得:
设各运动单位的发放是相互独立的,可得:
f(Hi(t1),Hj(t2))=f(Hi(t1))*f(Hj(t2))
由于时间长度I很短,f(Hi(t1)),f(Hj(t2))在I中变化很小,可作如下近似:
f(Hi(t1),=f(Hi(t), f(Hj(t2)=f(Hj(t))(5)
式中t为I中的任一时刻,于是:
这里,积分项中把A从概率密度函数内删除,因为两类MUAP出现在时刻t1和t2,t1.t2∈I即意味着事件A的发生。
同样可以得到:
略去公共项f(X,A),实际只需比较如下似然函数:
