CN103761424B - 基于二代小波和独立分量分析肌电信号降噪与去混迭方法 - Google Patents

Abstract

本发明涉及一种基于二代小波变换和独立分量分析相结合的肌电信号降噪与去混迭方法。目前有的信号去混迭方法运用于肌电信号时显现出一些不足。本发明首先采用二代小波降噪算法先滤除信号中的噪声。二代小波降噪算法运用于肌电信号，使降噪后的信号变得平滑，波形当中不必要的震荡得到抑制，肌电信号特征更加明显。再对降噪后的含混信号进行ICA分离，可快速有效的去除信号中的混迭成分。利用本发明对肌电信号做以上预处理，可大大去除信号中的干扰，为后续的肌电信号的特征提取，动作识别等研究提供了便捷。

Description

基于二代小波和独立分量分析肌电信号降噪与去混迭方法

技术领域

本发明属于多通道信号去混迭领域，涉及一种应用于肌电信号的消噪与去混迭方法。

背景技术

肌电信号包含了丰富的肌肉运动的信息，它能够反映肢体的运动模式。目前，肌电信号已经被广泛应用于临床诊断、康复工程、运动医学等领域。表面肌电信号(sEMG)是一种在皮肤表面募集到的，伴随着肢体运动导致肌肉收缩而产生的生物电信号。基于肌电信号的人机接口通过对表面肌电信号的处理和模式识别可识别出人体不同的动作模式，控制外部环境设备工作，如基于肌电信号的仿生假手的控制，基于肌电信号的跌倒辨识。

基于肌电信号的人机接口研究中，要通过肌电信号对多个动作实现模式识别，往往需要用多个电极同时进行多路肌电信号的采集。由于人体本身存在天然电流，各神经元之间又相互连接，以及采样电极相邻通道之间存在耦合串扰等因素，使得多通道肌电信号采集到的信号存在相邻肌电信号的混迭成分。为了得到纯粹的肌电信号，进行进一步研究分析，我们需要一种有效的方法对肌电信号进行去混迭处理。

罗志增等针对脑电信号中的噪声及夹杂的眼电、心电等伪迹，采用小波分解重构消噪和基于熵估计的RADICAL算法进行消除。对3路脑电观测信号进行小波消噪和白化处理后，通过雅可比旋转矩阵分别对其中的两两组合信号用RADICAL算法进行分离，得出最优分离矩阵完成盲源分离。张洪渊等利用随机变量概率密度函数非参数估计的核函数法,对混合信号的概率密度函数及其导数进行估计,可以直接应用于所有以非线性函数代替评价函数的信号源盲分离算法。以上方法都是不错的信号盲分离算法，但是基于sEMG是一种非线性非平稳信号，以及人体活动时各肌肉的协同工作的特点，使得以上方法运用于sEMG时显现出一些不足。

sEMG的有用能量主要分布在10-500Hz之间，非常微弱(其幅值仅为μV级)，极易受到噪声的干扰，在对sEMG作去混迭处理前必须滤除其中的噪声。以改善肌电信号的盲源分离效果。

传统的小波降噪方法对信号进行多尺度分解，采用一个小波基函数去逼近不同尺度上的信号，在小波基函数与逼近信号的差异处会产生细节信息，这些细节信号在阈值处理时会被当作噪声被滤除，使降噪后的信号丢失部分有用信息。这样，使得常规的小波降噪方法无法准确的刻画非平稳信号。

综上所述，已有的信号消噪与去混迭方法在运用于肌电信号时都存在不足之处。本发明提出了一种基于二代小波变换和独立分量分析(ICA)相结合的肌电信号降噪与去混迭方法。

发明内容

本发明的目的就是针对现有技术的不足，提出一种基于二代小波变换和独立分量分析(ICA)相结合的肌电信号降噪与去混迭方法。第二代小波变换继承了传统小波的多分辨率的特性；不依赖傅立叶变换；小波变换后的系数是整数；基于多项式内插的思想，所有的运算都在空间域进行，从而摆脱了对频域的依赖。利用二代小波进行信号降噪能滤除高频噪声并保留有用信息，还能较好地保留sEMG的边缘特性。

独立分量分析(ICA)是根据从传感器监测到的观测信号来估计未知的源信号和混合矩阵的。它从信号高阶统计特性出发，主要目的按是为了提取信号的独立成分，与已有的主分量分析(PCA)方法相比，它不仅可去除分量间的低阶相关性，而且还能够去除分量间的高阶相关性，使之相互独立。因此，ICA是一种有效的盲源分离方法，相对其他方法能更全面地体现数据间的本质结构。

所以，二代小波变换和独立分量分析(ICA)相结合的方法可以有效的降低肌电信号中的噪声并去除相邻通道信号的混迭成分。

为了实现以上目的，本发明方法主要包括以下步骤：

步骤一，对仪器采集到的含噪声和线性混迭的M路观测信号，X＝[x₁,x₂,…,x_M]^T进行二代小波分解。

对每一路信号x[n]进行二代小波分解的具体步骤如下：

①分裂：将信号序列x[n]分成两个互不相交的子集，通常按奇样本x_odd[n]、偶样本x_even[n]分成两个等长的序列。

②预测：采用一个与数据结构无关的预测算子P，保持偶样本不变，由偶样本估计奇样本x_odd[n]＝P(x_even[n])。如果信号具有局部相关结构，偶样本集合和奇样本集合一定是高度相关的，可以以预测算子P用一个子集估计另一个子集。原值与预测值之差作为小波系数d[n]，表示为d[n]＝x_odd[n]-P(x_even[n])。

③修正：为了保持存在于原始信号中的低频尺度特性，引入修正算子U，以得到下一尺度上的尺度系数：c[n]＝x_even[n]-U(d[n])。

上面这三步构成一个提升过程，对输出c[n]重复以上提升过程，可构成一个完整的离散小波变换。得到表示原始数据的低频成分的尺度系数c_j[n]和表示原始数据的高频成分的小波系数d_j[n]。j为二代小波分解的层级，最高层数记为L。

步骤二，对步骤一所得到的小波系数d_j[n]进行阈值处理得到消噪后的高频系数d'_j[n]，阈值函数如下：

$d_j^{\&prime;}\&lsqb;n\&rsqb;=\left\{\begin{array}{c}{d}_j\&lsqb;n\&rsqb;-\mathrm{\&alpha;}\mathrm{\&tau;}\left(\frac{2+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}{1+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}\right),|d_j\&lsqb;n\&rsqb;|\&GreaterEqual;\mathrm{\&tau;}\\ 0,|d_j\&lsqb;n\&rsqb;|<\mathrm{\&tau;}\end{array}\right.$

式中α∈[0,1]为一可变参数。当α＝0时，阈值函数等效于原始信号。阈值σ＝1/0.6745med(|x|)，只要选取合适的参数α就可以得到平滑的降噪信号，抑制波形当中不必要的振荡。

步骤三,将步骤二所得的d'_j[n]与步骤一所得的c_j[n]进行小波重构得到消噪后的信号x'[n]。M路信号都经过降噪处理后得到消噪后的肌电信号X'＝[x₁',x₂',…,x_M']^T。

步骤四,对步骤三得到的X'进行去均值处理：

X”＝X'-E(X')

其中，去均值后的观测数据X”的协方差矩阵可以分解为C_x＝E[X”X”^T]＝UλU^T，其中U为C_x的特征向量矩阵，λ为C_x的特征值矩阵，λ＝Diag(λ₁,λ₂,...,λ_N)。

对X”进行白化处理：

Z＝BX”

其中，白化矩阵B＝λ^-1/2U^T。去均值和白化过程使信号成为零均值和具有单位方差且各分量互不相关的矢量Z＝[z₁,z₂,…,z_M]^T。

步骤五,解混矩阵W的第i列记为解混向量w_i。将步骤四得到的Z带入FastICA算法的迭代式得：

w_i(k+1)＝E{Z(w_i(k)^T)³Z}-3γE{ZZ^T}w_i(k)

k表示迭代次数。随机产生w_i(0)，且||w_i(0)||₂＝1。γ为一大于零的可变参数。

步骤六，为保证每次提取出来的都是之前从未提取过的独立分量，在计算w_i且i≥2时，需要添加正交化步骤，并把步骤五中得到的w_i(k+1)减去已提取的独立分量，即按下式处理：

$w_i(k+1)\&LeftArrow;w_i(k+1)-\underset{j=1}{\overset{i-1}{\&Sigma;}}\left(w_i^T\right(k+1\left)w_j\right)w_j,i\&GreaterEqual;2$

步骤七，为确保分离出的独立分量具有单位能量，需要对步骤六得到的w_i(k+1)进行归一化处理：

$w_i(k+1)\&LeftArrow;\frac{w_i(k+1)}{\left|\right|w_i(k+1)\left|\right|}$

步骤八,如果|w_i(k+1)^Tw_i(k)|收敛于1,则算法收敛,得到一个解混向量w_i，i＝i+1,否则返回步骤五，直至算法收敛。

步骤九,重复步骤五，六，七，八得到解混矩阵W＝[w₁,w₂,...,w_k]。可求出解混信号：

Y＝WZ

即Y＝[y₁,y₂,...,y_k]为最终降噪去混迭的信号。

本发明的有益效果：在基于肌电信号的人体动作模式识别的研究中，往往需要用多个电极同时进行多路信号的采集。由于人体本身存在天然电流，各神经元之间又相互连接，以及采样电极相邻通道之间存在耦合串扰等因素，使得多通道肌电信号采集到的信号存在相邻肌电信号的混迭成分。利用本发明对含混信号进行ICA分离，可快速有效的去除信号中的混迭成分。同时，肌电信号本身具有非线性、非平稳，非常微弱极易受到噪声干扰的特点，噪声的加入会影响模式识别的结果。本发明在对sEMG作去混迭处理前采用二代小波降噪算法先滤除信号中的噪声。二代小波降噪算法运用于肌电信号，使降噪后的信号变得平滑，波形当中不必要的震荡得到抑制，肌电信号特征更加明显。利用本发明对肌电信号做以上预处理，可大大去除信号中的干扰，以提高模式识别的准确度。

附图说明

图1为本发明的实施流程图；

图2为本发明实施例的3路原始含有噪声和混迭的肌电信号图；

图3为本发明实施例经过二代小波降噪的结果图；

图4为本发明实施例的降噪去混迭结果图。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，本实施例包括如下步骤：

步骤一，对仪器采集到的含噪线性混迭的M路观测信号X＝[x₁,x₂,…,x_M]^T进行二代小波分解，在本实施例中M＝3，肌电信号如图2所示。

对每一路信号x[n]进行二代小波分解的具体步骤如下：

①分裂：将信号序列x[n]分成两个互不相交的子集，通常按奇样本x_odd[n]、偶样本x_even[n]分成两个等长的序列。

②预测：采用一个与数据结构无关的预测算子P，保持偶样本不变，由偶样本估计奇样本x_odd[n]＝P(x_even[n])。如果信号具有局部相关结构，偶样本集合和奇样本集合一定是高度相关的，可以以预测算子P用一个子集估计另一个子集。原值与预测值之差作为小波系数d[n]，表示为d[n]＝x_odd[n]-P(x_even[n])。

③修正：为了保持存在于原始信号中的低频尺度特性，引入修正算子U，以得到下一尺度上的尺度系数：c[n]＝x_even[n]-U(d[n])。

上面这三步构成一个提升过程，对输出c[n]重复以上提升过程，可构成一个完整的离散小波变换。得到表示原始数据的低频成分的尺度系数c_j[n]和表示原始数据的高频成分的小波系数d_j[n]。j为二代小波分解的层级，最高层数记为L，在本实施例中L＝4。

步骤二，对步骤一所得到的小波系数d_j[n]进行阈值处理得到消噪后的高频系数d'_j[n]，阈值函数如下：

$d_j^{\&prime;}\&lsqb;n\&rsqb;=\left\{\begin{array}{c}{d}_j\&lsqb;n\&rsqb;-\mathrm{\&alpha;}\mathrm{\&tau;}\left(\frac{2+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}{1+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}\right),|d_j\&lsqb;n\&rsqb;|\&GreaterEqual;\mathrm{\&tau;}\\ 0,|d_j\&lsqb;n\&rsqb;|<\mathrm{\&tau;}\end{array}\right.$

式中α∈[0,1]为一可变参数。当α＝0时，阈值函数等效于原始信号。阈值σ＝1/0.6745med(|x|)，只要选取合适的参数α就可以得到平滑的降噪信号，抑制波形当中不必要的振荡。在本实施例中α＝0.5，N＝2020。

步骤三,将步骤二所得的d'_j[n]与步骤一所得的c_j[n]进行小波重构得到消噪后的信号x'[n]。M路信号都经过降噪处理后得到X'＝[x₁',x₂',…,x_M']^T，如图3所示。

步骤四,对步骤三得到的X'进行去均值处理:

X”＝X'-E(X')

其中，去均值后的观测数据X”的协方差矩阵可以分解为C_x＝E[X”X”^T]＝UλU^T，其中U为C_x的特征向量矩阵，λ为C_x的特征值矩阵，λ＝Diag(λ₁,λ₂,...,λ_N)。

对X”进行白化处理：

z＝BX”

其中，白化矩阵B＝λ^-1/2U^T。去均值和白化过程使信号成为零均值和具有单位方差且各分量互不相关的矢量Z＝[z₁,z₂,…,z_M]^T。

步骤五,解混矩阵W的第i列记为解混向量w_i。将步骤四得到的Z带入FastICA算法的迭代式得：

w_i(k+1)＝E{Z(w_i(k)^T)³Z}-3γE{ZZ^T}w_i(k)

k表示迭代次数。随机产生w_i(0)，且||w_i(0)||₂＝1。γ为一大于零的可变参数。在本实施例中γ＝0.8。

步骤六，为保证每次提取出来的都是之前从未提取过的独立分量，在计算w_i且i≥2时，需要添加正交化步骤，并把步骤五中得到的w_i(k+1)减去已提取的独立分量，即按下式处理：

$w_i(k+1)\&LeftArrow;w_i(k+1)-\underset{j=1}{\overset{i-1}{\&Sigma;}}\left(w_i^T\right(k+1\left)w_j\right)w_j,i\&GreaterEqual;2$

步骤七,为确保分离出的独立分量具有单位能量，需要对步骤六得到的w_i(k+1)进行归一化处理：

$w_i(k+1)\&LeftArrow;\frac{w_i(k+1)}{\left|\right|w_i(k+1)\left|\right|}$

步骤八,如果|w_i(k+1)^Tw_i(k)|收敛于1，则算法收敛,得到一个解混向量w_i，i＝i+1,否则返回步骤五，直至算法收敛。

步骤九,重复步骤五，六，七，八得到解混矩阵W＝[w₁,w₂,...,w_k]。可求出解混信号：

Y＝WZ

即Y＝[y₁,y₂,...,y_k]为最终降噪去混迭的信号，如图4所示。

Claims (1)

1.基于二代小波和独立分量分析肌电信号降噪与去混迭方法，其特征在于该方法包括如下步骤：

步骤一，对仪器采集到的含噪声和线性混迭的M路观测信号，X＝[x₁,x₂,…,x_M]^T进行二代小波分解；

对每一路信号x[n]进行二代小波分解的具体步骤如下：

①分裂：将信号序列x[n]分成两个互不相交的子集，按奇样本x_odd[n]、偶样本x_even[n]分成两个等长的序列；

②预测：采用一个与数据结构无关的预测算子P，保持偶样本不变，由偶样本估计奇样本x_odd[n]＝P(x_even[n])；如果信号具有局部相关结构，偶样本集合和奇样本集合一定是高度相关的，以预测算子P用一个子集估计另一个子集；原值与预测值之差作为小波系数d[n]，表示为d[n]＝x_odd[n]-P(x_even[n])；

③修正：为了保持存在于原始信号中的低频尺度特性，引入修正算子U，以得到下一尺度上的尺度系数：c[n]＝x_even[n]-U(d[n])；

上面这三步构成一个提升过程，对输出c[n]重复以上提升过程，构成一个完整的离散小波变换；得到表示原始数据的低频成分的尺度系数c_j[n]和表示原始数据的高频成分的小波系数d_j[n]；j为二代小波分解的层级，最高层数记为L；

步骤二，对步骤一所得到的小波系数d_j[n]进行阈值处理得到消噪后的高频系数d'_j[n]，阈值函数如下：

$d_j^{\&prime;}\&lsqb;n\&rsqb;=\left\{\begin{array}{c}{d}_j\&lsqb;n\&rsqb;-\mathrm{\&alpha;}\mathrm{\&tau;}\left(\frac{2+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}{1+{\mathrm{\&tau;}}^{\mathrm{\&alpha;}}}\right)\\ 0,|d_j\&lsqb;n\&rsqb;|<\mathrm{\&tau;}\end{array}\right.,|d_j\&lsqb;n\&rsqb;|\&GreaterEqual;\mathrm{\&tau;}$

式中α∈[0,1]为一可变参数；当α＝0时，阈值函数等效于原始信号；阈值

步骤三,将步骤二所得的d'_j[n]与步骤一所得的c_j[n]进行小波重构得到消噪后的信号x'[n]；M路信号都经过降噪处理后得到消噪后的肌电信号X'＝[x₁',x₂',…,x_M']^T；

步骤四,对步骤三得到的X'进行去均值处理:

X”＝X'-E(X')

其中，去均值后的观测数据X”的协方差矩阵分解为C_x＝E[X”X”^T]＝UλU^T，其中U为C_x的特征向量矩阵，λ为C_x的特征值矩阵，λ＝Diag(λ₁,λ₂,...,λ_N)；

对X”进行白化处理：

Z＝BX”

其中，白化矩阵B＝λ^-1/2U^T；去均值和白化过程使信号成为零均值和具有单位方差且各分量互不相关的矢量Z＝[z₁,z₂,…,z_M]^T；

步骤五,解混矩阵W的第i列记为解混向量w_i；将步骤四得到的Z带入FastICA算法的迭代式得：

w_i(k+1)＝E{Z(w_i(k)^T)³Z}-3γE{ZZ^T}w_i(k)

k表示迭代次数；随机产生w_i(0)，且||w_i(0)||₂＝1；γ为一大于零的可变参数；

步骤六，为保证每次提取出来的都是之前从未提取过的独立分量，在计算w_i且i≥2时，添加正交化步骤，并把步骤五中得到的w_i(k+1)减去已提取的独立分量，即按下式处理：

$w_i(k+1)\&LeftArrow;w_i(k+1)-\underset{j=1}{\overset{i-1}{\&Sigma;}}\left(w_i^T\right(k+1\left)w_j\right)w_j,i\&GreaterEqual;2;$

步骤七，对步骤六得到的w_i(k+1)进行归一化处理：

$w_i(k+1)\&LeftArrow;\frac{w_i(k+1)}{\left|\right|w_i(k+1)\left|\right|}$

步骤八,如果|w_i(k+1)^Tw_i(k)|收敛于1,则算法收敛,得到一个解混向量w_i，i＝i+1,否则返回步骤五，直至算法收敛；

步骤九,重复步骤五，六，七，八得到解混矩阵W＝[w₁,w₂,...,w_k]；可求出解混信号：

Y＝WZ

即Y＝[y₁,y₂,...,y_k]为最终降噪去混迭的信号。