基于独立向量分析的高密度表面肌电信号消噪方法
阅读说明：本技术 基于独立向量分析的高密度表面肌电信号消噪方法 (High-density surface electromyographic signal denoising method based on independent vector analysis ) 是由 陈勋 王堃 吴乐 张旭 陈香 于 2019-10-15 设计创作，主要内容包括：本发明公开了一种基于独立向量分析的高密度表面肌电信号消噪方法,其步骤包括：1、首先将高密度表面肌电信号通过延时构造多个数据集；2、利用独立向量分析法进行联合盲源分离,得到每个数据集相对应的源信号矩阵和解混矩阵；3、选择高密度表面肌电信号相对应的源信号矩阵和解混矩阵；4、计算源信号矩阵中各个独立源成分肌肉收缩段和静息段的强度比；5、设定消噪阈值,将强度比低于消噪阈值的独立源成分置零；6、盲源分离逆变换得到消除噪声后的干净高密度表面肌电信号。本发明基于独立向量分析,通过设定消噪阈值,在去除工频干扰和高斯白噪声的同时,保证了肌电信息在处理过程中尽可能不丢失。(The invention discloses a high-density surface electromyographic signal denoising method based on independent vector analysis, which comprises the following steps: 1. firstly, constructing a plurality of data sets by delaying high-density surface electromyographic signals; 2. performing combined blind source separation by using an independent vector analysis method to obtain a source signal matrix and a demixing matrix corresponding to each data set; 3. selecting a source signal matrix and a demixing matrix corresponding to the high-density surface myoelectric signals; 4. calculating the intensity ratio of each independent source component muscle contraction section to each independent source component muscle rest section in the source signal matrix; 5. setting a noise elimination threshold value, and setting the independent source components with the intensity ratio lower than the noise elimination threshold value to be zero; 6. and (4) carrying out blind source separation inverse transformation to obtain a clean high-density surface myoelectric signal after noise elimination. The method is based on independent vector analysis, and guarantees that the electromyographic information is not lost as much as possible in the processing process while removing power frequency interference and Gaussian white noise by setting a noise elimination threshold.)
技术领域
本发明属于肌电信号处理技术领域，具体涉及一种基于独立向量分析，从高密度表面肌电信号中自动识别工频干扰和高斯白噪声的方法，主要应用于肌电信号的预处理。
背景技术
表面肌电信号(Surface Electromyography,sEMG)是通过附着在人体皮肤表面的电极记录下来的、能够反映神经肌肉系统相关活动信息的微弱电生理信号。sEMG信号已经被广泛应用于神经生理学、临床医学、运动医学和康复医学等领域。高密度表面肌电信号(High Density sEMG,HD-sEMG)技术是单通道表面肌电信号技术的一种时空变体，其具有大量紧密相间的电极通道，在较宽的肌肉区域内提供了肌电强度的时空分布的信息。这些信息能够改善肌力估计、肌电分解和肌电假肢控制等。然而，肌电信号作为微弱的电生理信号，经常受到如工频干扰、高斯白噪声和运动伪迹等噪声的干扰，影响后续对肌电分析的准确性。
在过去数十年中，研究人员已经提出多种方法用于消除肌电信号中的噪声。一种常用的方法是利用工频陷波器消除50Hz的工频干扰，但同时也使肌电信号的有效信息受到较大的影响，因为肌电信号的频率范围主要分布在10Hz～500Hz之间，在滤波的同时也对肌电信号造成不可挽回的损失。小波变换具有良好的时频分析能力，具有传统去噪方法不可比拟的优越性，在信号去噪中被广泛应用。然而小波去噪的效果与小波基函数的选择直接相关，同时在信号分析的过程中，一旦小波基函数被选定后将无法改变，因此小波变换的自适应比较差。经验模态分解将信号分解为多个固有模态函数的叠加，每个固有模态函数具有不同尺度的特征，具有很好的自适应性，但存在着模态混合的缺点。集合经验模态分解在信号分解之前将白噪声添加到肌电信号中，虽然克服了模态混合的缺点，但同时信号也会受到残留噪声的污染。
此外，盲源分离方法也常被用于生物医学信号的噪声消除，盲源分离的目的是从观测信号中分离出不同的信号源并得到混合矩阵。其中最常用的盲源分离方法是独立成分分析(Independent Component Analysis,ICA)和典型相关分析(Canonical CorrelationAnalysis,CCA)。盲源分离方法应用于肌电信号中噪声消除的原理是，该方法从多通道的观测信号中提取出一系列独立源成分，然后通过将这些源成分中代表噪声的源成分置零，最后通过逆盲源分离得到消噪后的肌电信号。CCA利用二阶统计量将肌电信号分解成互不相关而自相关性最大的典型变量。由于高斯白噪声的频带宽，与肌电信号相比具有相当低的自相关系数，所以利用二阶统计特性可以很好地分离高斯白噪声成分与肌电成分。但仅仅要求典型变量之间互不相关并不能完全分离出工频干扰分量，因此，CCA方法存在肌电信号分量还混叠有工频干扰的问题。ICA利用高阶统计量将肌电信号分解成相互独立的分量，其可以较好地将工频干扰分量和肌电信号分量分离，但分离出的工频干扰分量仍含有较多的肌电信息。此外，由于ICA无法将多个高斯成分彼此之间分离开来，因此直接将分离出来的高斯成分置零并不能达到消除高斯白噪声的目的。
发明内容
本发明为了克服现有技术的不足之处，提出一种基于独立向量分析的高密度表面肌电信号消噪方法，以期能完全去除工频干扰和高斯白噪声对肌电信号的影响，同时尽可能地保留肌电分量的信息不丢失，从而能提高肌电信号的分析准确性，为肌电去噪提供了一种新思路。
本发明为解决技术问题，采用如下技术方案：
本发明一种基于独立向量分析的高密度表面肌电信号消噪方法的特点是按如下步骤进行：
步骤一：由肌电测量设备采集并记录t时刻N通道的高密度表面肌电信号矩阵，记为：X(t)＝[x1(t),x2(t),...,xn(t),...,xN(t)]T，xn(t)为t时刻第n通道的肌电信号，若t∈[0,t1]，则高密度表面肌电信号矩阵X(t)为不含肌电活动的静息段；若t∈(t1,t2]，则高密度表面肌电信号矩阵X(t)为含有肌电活动的肌肉收缩段，T为矩阵的转置，t1为肌肉开始收缩时刻，t2为采集结束时刻，1≤n≤N,0≤t≤t2；
步骤二：将所述高密度表面肌电信号矩阵X(t)通过延时k-1个点构造K个数据集X＝[X[1],X[2],…,X[k]…,X[K]]，其中，
步骤三：利用独立向量分析法对所述K个数据集X进行联合盲源分离，得到所述高密度表面肌电信号矩阵X(t)的t时刻源信号矩阵S(t)以及逆混合矩阵W*，其中，S(t)＝[s1(t),s2(t),...,sn(t),...,sN(t)]T，sn(t)表示t时刻第n个源信号，并有：S(t)＝W*X(t)；
步骤四：利用式(1)计算所述源信号矩阵S(t)中第n个独立源信号sn(t)的肌肉收缩段和静息段的强度比IRn：
式(1)中，cn(i)表示t时刻第n个源信号sn(t)的肌肉收缩段的第i个采样点；qn(i)表示t时刻第n个源信号sn(t)静息段的第i个采样点，m1＝t1×fs，m2＝(t2-t1)×fs，fs为信号的采样频率；
步骤五：设置消噪阈值，将强度比IRn小于消噪阈值的第n个源信号sn(t)判断为高斯白噪声或工频干扰成分，并将第n个源信号sn(t)置零；将强度比IRn大于消噪阈值的第n个源信号sn(t)判断为肌电信息成分，并保留第n个源信号sn(t)；
步骤六：重复步骤四和步骤五，得到噪声分量置零后的t时刻源信号矩阵记为
式(3)中，(W*)-1为逆混合矩阵W*的逆矩阵。
本发明所述的高密度表面肌电信号消噪方法的特点也在于，所述步骤三是按如下过程进行：
步骤3.1、利用多元高斯分布的独立向量分析法对所述K个数据集X进行联合盲源分离，得到K个初始解混矩阵
步骤3.2、利用初始解混矩阵对多元拉普拉斯分布的独立向量分析法进行初始化，并对所述K个数据集X进行联合盲源分离，得到K个解混矩阵W＝[W[1],W[2],...,W[k]...,W[K]]，W[k]为通过多元拉普拉斯分布的独立向量分析法得到的第k个数据集X[k]的解混矩阵；
步骤3.3、利用式(3)得到第k个数据集X[k]的待估计的源信号矩阵
式(3)中，k＝1,2,...,K；当k＝1时，为所述高密度表面肌电信号矩阵X(t)通过独立向量分析法得到的t时刻源信号矩阵，记为S(t)＝[s1(t),s2(t),...,sn(t),...,sN(t)]T，其中sn(t)为t时刻第n个源信号；W[1]为所述高密度表面肌电信号矩阵X(t)的逆混合矩阵，记为W*；X[1]为所述高密度表面肌电信号矩阵X(t)。
本发明与CCA和ICA方法相比，不仅能去除工频干扰和高斯白噪声对肌电的影响，同时能最大程度地保留肌电信息，具体的有益效果体现在：
1、本发明步骤二和步骤三中，通过延时得到了多个数据集，并利用独立向量分析同时对多个数据集的肌电信号进行联合盲源分离，这一做法与ICA和CCA相比，能够为盲信号分离提供更多的源信号相关信息，促使工频干扰源、高斯白噪声源和肌电源更为精确地分离到不同的独立分量中去，因而在随后的步骤四和步骤五中，通过阈值自动地去除噪声分量时，肌电信号的成分损失更少。
2、本发明步骤三中，独立向量分析对多数据集信号进行联合盲源分离时，分离出的源信号在同一数据集内相互独立，在不同的数据集中，相对应的源信号有着最大的相关性。与ICA方法相比，其考虑了高斯白噪声的二阶统计特性，使得其对高斯白噪声的分离更加有效，同时其分离出的工频干扰分量也含有更少的肌电信息，减少了去噪过程中肌电信号信息的丢失。与CCA方法相比，其分离得到的工频干扰源信号与肌电信号源在统计上相互独立，使得其对于工频干扰源的分离更加彻底，解决了CCA方法中工频干扰去除不彻底的问题。
3、本发明步骤三中，相对于现有的基于ICA和CCA的盲源分离方法，独立向量分析方法利用对源分量向量为多元高斯分布和多元拉普拉斯分布的假设，将二阶统计量和高阶统计量统一在联合盲源分离的数学模型中，充分结合了ICA和CCA方法的优势，同时充分考虑了肌电信号与噪声信号各自的统计特性，能够很好地分离高斯白噪声源成分、工频干扰源成分和肌电信号源成分，更好地解决肌电信号中的噪声消除问题。
附图说明
图1为本发明方法流程图；
图2a为仿真高密度表面肌电信号前20个通道前2秒的示意图；
图2b为被两种噪声污染的高密度表面肌电信号前20个通道前2秒的示意图；
图3为通过本发明方法得到的源信号矩阵前20个通道前2秒的示意图；
图4为通过本发明方法得到的源信号强度比示意图；
图5为通过本发明方法消噪后的源信号前2秒的示意图；
图6为通过本发明方法消噪后的高密度表面肌电信号前20个通道前2秒的示意图
图7为本发明方法和CCA、ICA方法，其去噪性能RRMSE的比较图；
图8a为本发明方法和CCA、ICA方法处理仿真肌电，其去噪性能信噪比的比较图
图8b为本发明方法和CCA、ICA方法处理仿真肌电，其去噪后信号的波形比较图；
图9a为48通道的真实高密度表面肌电信号前20个通道前2秒的示意图；
图9b为通过本发明方法去噪后的48通道表面肌电信号前20个通道前2秒的示意图；
图10a为本发明方法和CCA、ICA方法处理真实肌电，其去噪性能信噪比的比较图；
图10b为本发明方法和CCA、ICA方法处理真实肌电，其去噪后信号的波形比较图；
图11a为去噪前信号的频谱图；
图11b为ICA方法去除噪声前后的频谱局部对比图；
图11c为CCA方法去除噪声前后的频谱局部对比图；
图11d为IVA方法去除噪声前后的频谱局部对比图。
具体实施方式
本实施例中，如图1所示，一种基于独立向量分析的高密度表面肌电信号消噪方法是：首先将高密度表面肌电信号通过延时得到多个数据集，再用独立向量分析的方法对多个数据集进行联合盲源分离；然后得到高密度表面肌电信号的源信号矩阵和逆混合矩阵；再计算源信号矩阵中每个独立源成分肌肉收缩段和静息段的强度比，将低于消噪阈值的源信号置零，高于消噪阈值的源信号保留；最后盲源分离逆变换得到消除噪声后的高密度表面肌电信号。
下面分别通过模拟肌电信号与真实肌电信号为例，结合附图来说明具体的实施方式。
1模拟肌电信号
在这一部分，将说明两个实例，第一个实例是将介绍本发明的具体实施方式，第二个实例是将本发明与ICA(FastICA)和CCA处理方法相比较。
(1)实例一
步骤一：由肌电测量设备采集并记录t时刻N通道的高密度表面肌电信号矩阵，记为：X(t)＝[x1(t),x2(t),...,xn(t),...,xN(t)]T，xn(t)为t时刻第n通道的肌电信号，若t∈[0,t1]，则高密度表面肌电信号矩阵X(t)为不含肌电活动的静息段；若t∈(t1,t2]，则高密度表面肌电信号矩阵X(t)为含有肌电活动的肌肉收缩段，T为矩阵的转置，t1为肌肉开始收缩时刻，t2为采集结束时刻，1≤n≤N,0≤t≤t2；
本实施例中，由肌电测量设备采集并记录t时刻N＝64通道的高密度表面肌电信号矩阵，记为：X(t)＝[x1(t),x2(t),...,xn(t),...,x64(t)]T，xn(t)为t时刻第n通道的肌电信号，若t∈[0,0.5]，则高密度表面肌电信号矩阵X(t)为不含肌电活动的静息段；若t∈(0.5,5.5]，则高密度表面肌电信号矩阵X(t)为含有肌电活动的肌肉收缩段，T为矩阵的转置，t1＝0.5为肌肉开始收缩时刻，t2＝5.5为采集结束时刻，1≤n≤64,0≤t≤5.5；如图2b所示的X(t)＝XEMG(t)+XWGN(t)+XPLI(t)，其中，XEMG(t)＝[xEMG1(t),xEMG2(t),...,xEMG64(t)]T表示64通道的仿真高密度表面肌电信号矩阵，由图2a可以清楚地看出，t∈[0,0.5]信号幅值为零，为不含肌电活动的静息段，0.5秒为肌肉开始收缩时刻，信号幅值逐渐增大，t∈(0.5,5.5]为含有肌电活动的肌肉收缩段，XWGN(t)＝[xWGN1(t),xWGN2(t),...,xWGN64(t)]T表示64通道的仿真高斯白噪声信号，XPLI(t)＝[xPLI1(t),xPLI2(t),...,xPLIch(t)]T,15≤ch≤20，表示仿真的ch通道的工频干扰信号；X(t)的具体仿真过程为：首先添加64通道的高斯白噪声XWGN(t)，使得添加高斯白噪声后，信号X(t)的64通道的信噪比服从N(20,3)的正态分布；然后添加工频干扰XPLI(t)，随机从64通道中选择ch个通道添加工频干扰，使得信号X(t)中添加工频干扰的通道服从信噪比为N(7,3)的正态分布。
步骤二：将高密度表面肌电信号矩阵X(t)通过延时k-1个点构造K个数据集X＝[X[1],X[2],…,X[k]…,X[K]]，其中，
本实施例中，取K＝21，将高密度表面肌电信号矩阵X(t)通过延时k-1个点构造21个数据集X＝[X[1],X[2],…,X[k]…,X[21]]，其中，
步骤三：利用独立向量分析法对K个数据集X进行联合盲源分离，得到高密度表面肌电信号矩阵X(t)的t时刻源信号矩阵S(t)以及逆混合矩阵W*，其中，S(t)＝[s1(t),s2(t),...,sn(t),...,sN(t)]T，sn(t)表示t时刻第n个源信号，并有：S(t)＝W*X(t)；其具体步骤为：
步骤3.1、利用多元高斯分布的独立向量分析法对K个数据集X进行联合盲源分离，得到K个初始解混矩阵
步骤3.2、利用初始解混矩阵
步骤3.3、利用式(3)得到第k个数据集X[k]的待估计的源信号矩阵
式(3)中，k＝1,2,...,K；当k＝1时，为高密度表面肌电信号矩阵X(t)通过独立向量分析法得到的t时刻源信号矩阵，记为S(t)＝[s1(t),s2(t),...,sn(t),...,sN(t)]T，其中sn(t)为t时刻第n个源信号；W[1]为高密度表面肌电信号矩阵X(t)的逆混合矩阵，记为W*；X[1]为高密度表面肌电信号矩阵X(t)。
本实施例中，利用独立向量分析法对21个数据集X进行联合盲源分离，得到高密度表面肌电信号矩阵X(t)的t时刻源信号矩阵S(t)，如图3所示，由于多元高斯分布的独立向量分析法考虑了数据集间相对应的源信号的二阶相关性，所以其分离出的源信号按照自相关系数从大到小自动排列，从图3可以看出自相关系数最大的工频干扰成分为分离出的第一个源信号，分布在后面的独立源信号为自相关系数最低的高斯白噪声成分，介于两者之间的独立源信号为肌电信号成分。又由于多元拉普拉斯分布的独立向量分析法利用了高阶统计量，其分离出来的源信号在同一个数据集内相互独立，所以其可以将工频干扰成分与肌电信号成分完全分离，从图3中可以看出，只有第一个源信号为工频干扰成分，其它肌电信号成分都不含有工频干扰成分。
步骤四：利用式(1)计算源信号矩阵S(t)中第n个独立源信号sn(t)的肌肉收缩段和静息段的强度比IRn：
式(1)中，cn(i)表示t时刻第n个源信号sn(t)的肌肉收缩段的第i个采样点；qn(i)表示t时刻第n个源信号sn(t)静息段的第i个采样点，m1＝t1×fs，m2＝(t2-t1)×fs，fs为信号的采样频率；
本实施例中，利用式(1)计算出源信号矩阵S(t)中64个独立源信号的肌肉收缩段和静息段的强度比如图4所示，图中各个独立源信号的IR值按降序排列。IR值越大说明源信号肌肉收缩段所含的肌电活动越多，如图3所示的第2个源信号至第11个源信号，可以明显看出信号静息段只含少量噪声，肌肉收缩段含有大量肌电活动；IR值越小说明源信号所含噪声成分更多，因为高斯白噪声和工频干扰在肌肉收缩段和静息段的信号幅值较为稳定，没有明显的变化，如图3所示的第一个工频干扰源信号以及后面的高斯白噪声成分，两种噪声在静息段和肌肉收缩段的幅值都很稳定，所以计算出的IR值会在1附近。结合图3和图4所示，可以看出噪声分量的IR值都非常低，所以设定消噪阈值能够在去掉噪声分量的同时，保留大部分的肌电分量，尽可能保证肌电信息不损失。
步骤五：设置消噪阈值，将强度比IRn小于消噪阈值的第n个源信号sn(t)判断为高斯白噪声或工频干扰成分，并将第n个源信号sn(t)置零；将强度比IRn大于消噪阈值的第n个源信号sn(t)判断为肌电信息成分，并保留第n个源信号sn(t)；
本实施例中，设置消噪阈值为1.1，得到消噪后的源信号矩阵如图5所示，可以看出工频干扰和大部分高斯白噪声都被消除，保留下来的第1到第10个源信号明显含有丰富的肌电活动，而第11到第17个源信号虽然也含有白噪声，但同时也含有较多的肌电信息成分，去除这些源信号，肌电信号的信噪比反而会下降，所以保留这些源信号可以尽可能地保证肌电信息不损失；
步骤六：重复步骤四和步骤五，得到噪声分量置零后的t时刻源信号矩阵记为
式(3)中，(W*)-1为逆混合矩阵W*的逆矩阵。
本实施例中，得到消噪后的肌电信号
(2)实例二
为了量化评估本发明的效果，为此将本发明方法(IVA)与ICA(Fast ICA)和CCA两种盲源分离算法相比较。选择相对均方根误差(RRMSE)和信噪比(SNR)两个性能指标作为评价指标。相对均方根误差的定义如下：
信噪比的定义为：
SNRn表示第n个通道的信噪比，xn(i)表示X(t)中第n个通道的第i个样本点，qn(i)表示X(t)中第n个通道的静息段的第i个样本点，n＝1,2,...,N，m1＝t1*fs，m＝t2*fs，fs为信号的采样频率。
RRMSE的值越小越好，越小表示去噪后的信号和原始的干净肌电信号差别越小；SNR的值越大越好，越大表示噪声的功率越小，噪声去除的越干净。
按照实例1中的步骤，重复仿真实验30次，对30次仿真实验的RRMSE取平均，画出了三种方法的RRMSE对比图，如图7所示，从图中可以看出，ICA和CCA消除噪声的能力相当，但都明显弱于IVA方法，且IVA方法的RRMSE值随着数据集数量的增加呈缓慢下降的趋势，也就是说，IVA方法的消噪能力随着数据集数量的增多而缓慢增强。
图8a为一次仿真实验中三种方法消除噪声前后通道的信噪比情况,从图中可以看出IVA方法能够大幅提高所有通道的信噪比，而ICA和CCA方法在提升低信噪比通道的信噪比的同时，也降低了高信噪比通道的信噪比。图8b为第38通道去除噪声前后EMG信号的对比图，可以看出IVA算法能够很好去除高斯白噪声和工频干扰。而CCA不能完全去除工频干扰，ICA对高斯白噪声的分离不彻底，两者的缺点使得盲源分离逆变换后的信号仍然有明显的工频干扰和高斯白噪声干扰，也因此使高信噪比通道的信噪比降低。
2真实肌电信号
使用真实肌电数据作为实验对象，用IVA算法进行处理，并与ICA和CCA算法进行比较，评判发明方法的去噪效果。图9a是一段48通道的真实肌电信号，前0.5秒为静息段，后5秒为肌肉活动段。从图9a中可以看出该段信号被高斯白噪声和工频干扰的严重污染，同时可以观察到每个通道受干扰的程度不相同，所以消除工频干扰，提高通道的信噪比是十分有必要的。
按照本发明方法的步骤，对上述真实肌电信号进行处理，消噪后重建的肌电信号如图9b所示。通过对比图9a和图9b，可以发现本发明能将工频干扰去除得非常干净，同时能完整地保留肌电信号的信息。
此外，使用CCA和ICA方法对上述的真实肌电信号进行处理，得到去除噪声前后通道的信噪比如图10a所示。通过对比，可以发现CCA方法与ICA方法去除噪声的能力相当，ICA方法对于某些通道的信噪比提升并不明显，而本发明能够稳定地提升所有通道的信噪比。图10b为48通道的表面肌电信号第22通道去除噪声前后的信号在时域上的对比图，可以发现三种方法都能有效地去除工频干扰。但是CCA方法不能将工频干扰成分完全去除，ICA方法在去除工频干扰的同时，引入了少量的干扰，使得信噪比下降。
图11a为48通道的表面肌电信号第22通道去除噪声前的频谱图。通过频谱图可以看出，信号在50Hz受到严重的工频干扰；图11b为ICA方法去除噪声前后的频谱局部对比图，图11c为CCA方法去除噪声前后的频谱局部对比图，图11d为IVA方法去除噪声前后的频谱局部对比图；通过对比可知，经过本发明方法消除噪声后的信号在50Hz频率处的能量最低，去除工频干扰的效果最好，而ICA和CCA的方法还含有不同程度的工频干扰。此外，在其它频率附近，经过本发明方法消除噪声后的信号与原始信号的曲线最为接近，说明本发明方法在去除工频干扰的同时，没有损失肌电信号的其它信息，而ICA和CCA方法均在不同程度上损失了肌电信息。
综上所述，本发明能够解决高密度表面肌电信号中工频干扰和高斯白噪声的消除问题，而且不损失任何肌电信息。该方法适用于肌电分解和肌电假肢控制等肌电信号的预处理，与CCA和ICA方法相比，能够取得更好的消噪效果，对进一步研究神经肌肉系统的活动具有重要意义。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种磁场发生器支撑架