CN107688120A - 基于模糊熵的含噪信号处理方法及迭代奇异谱软阈值去噪方法 - Google Patents

Abstract

本发明公开了一种基于模糊熵的含噪信号处理方法及迭代奇异谱(SSA)软阈值去噪方法。该方法适用于含噪信号，假设长度N的含噪信号x_in＝{x₁,x₂,…,x_N}，并假设其中的加性白噪声与信号不相关，利用原始信号x_in构造d维矢量并且定义相似度以及模糊概率同样方法构造(d+1)维矢量及其相应的相似度和模糊概率模糊熵则定义为：对利用已知信号分解方法获得的分量，其所有分量的奇异谱分布定义为模糊熵谱。本发明利用混沌理论中量化***复杂度的模糊熵来表征噪声平面，为含噪信号的处理提供了更加有效的路径；其基于模糊熵谱的迭代奇异谱(SSA‑IST)软阈值去噪方法，去噪性能优于传统的截断奇异谱方法，及小波变换和经验模态分解去噪方法。

Description

基于模糊熵的含噪信号处理方法及迭代奇异谱软阈值去噪 方法

技术领域

本发明涉及信号去噪和滤波，具体涉及一种基于模糊熵的含噪信号处理方法及迭代奇异谱软阈值去噪方法。

背景技术

奇异谱分析(SSA)是一种结合了经典时间序列分析、线性代数、多变量统计学、动力***的高级信号处理方法。SSA的目的在于将待分析信号分解为多个具有物理意义的分量之和，比如趋势、震荡分量和噪声。根据这些分量，研究人员提出了不同的去噪、变点检测、缺失值插补、同步检测、特征提取和预测算法。为了滤除噪声，该方法首先计算延时协方差矩阵的特征值并按降序排列，这些特征值可能会从某个阶数开始形成一个较为平坦的拖尾状平面，即所谓的“噪声平面”。对应于奇异谱上最初陡峭部分较大特征值的分量构成了信号子空间的基底，而对应于噪声平面部分的分量则认为是白噪声而被舍弃。高的特征值对应于信号中基本震荡成分，最大的一个奇异值通常与缓变趋势相关。这种截断奇异谱的方法已被广泛用于语音、超声、多普勒雷达信号，生物医学信号(脑电，心电，肌电)和机械信号以及高光谱图像的去噪。为了获得较为满意的去噪效果，以前的研究和发明侧重于如何确定或查找噪声平面的最佳阶数。

截断SSA方法的最初出发点是任何信号都有一个固有的“噪声平面”，但是许多信噪比(SNR)较低的信号，特别是工程现场记录的信号，其奇异谱往往是平滑的指数函数形状，并不存在一个清晰地噪声平面。另外，奇异谱算法本身仅解决了如果确定一个信号和噪声混合测试矩阵的最佳低秩表达问题，并没有对如何获得低秩信号矩阵的最佳估计给出回答。因此，人们不应期望SSA去噪能获得最佳性能。相反，奇异谱采用保留部分分量，舍弃部分分量的二元法。形同滤除了信号的高频分量，并不能显著提高信噪比。

模糊熵是混沌理论中用来表征***复杂度的一个混沌不变量，在此，我们提出模糊熵谱的概念，并利用该谱来表征信号噪声平面，根据含噪信号模糊熵噪声谱发明了一种迭代SSA阈值去噪方法。

发明内容

本发明的目的在于：针对现有技术中存在的上述不足和欠缺，提供一种基于模糊熵的含噪信号处理方法谱及迭代奇异谱软阈值去噪方法。

本发明是通过以下技术方案实现的：

本发明的基于模糊熵的含噪信号处理方法是：对于含噪信号，假设长度N的含噪信号x_in＝{x₁,x₂,…,x_N}，并假设其中的加性白噪声与信号不相关，即：这里代表无噪信号，n为白噪声，利用原始信号x_in构造一系列d维矢量(1≤e≤N-d+1)，并且定义相似度以及模糊概率同样方法构造(d+1)维矢量及其相应的相似度和模糊概率模糊熵则定义为：这里d和r分别设置为2和0.2，N为待去噪信号长度，对利用任何信号分解方法获得的分量，其所有分量的奇异谱分布定义为模糊熵谱。

本发明的基于上述模糊熵谱的迭代奇异谱软阈值去噪方法，包括如下步骤：

(1)、SSA分解：噪声信号x_in＝{x₁,x₂,…,x_N}嵌入为m×d的Hankel矩阵，SVD分解Hankel矩阵为d个秩1矩阵之和，并重构为d个信号分量(c＝1，…d)；

(2)、计算模糊熵谱：按上述模糊熵谱定义计算SSA分量的模糊熵谱和原始噪声信号的模糊熵值；

(3)、设置软阈值：第一分量x₁不做任何阈值去噪，对模糊熵值小于噪声信号模糊熵值的分量x₂到x_k，设置较小阈值对剩余的x_k+1到x_d分量，设置较大阈值这里σ_c为分量x_c方差；

(4)、软阈值去噪：除第一分量的所有分量x_c(c＝2,…,d)，利用相应的阈值进行去噪，即每一分量的某噪声信号数值绝对值如小于该分量阈值，输出为零，如信号数值大于阈值，则输出为信号数值减去阈值，剩余情况输出为信号数值加上去阈值，软阈值去噪后所有分量(c＝2,…,d)与x₁之和为第一次估计信号估计噪声为

(5)、估计的方差，用估计信号作为输入信号，重复步骤(1)到(4)；

(6)、比较连续迭代获得的噪声方差，如噪声方差不再明显减小或达到指定迭代次数，迭代停止，否则重复(1)到(5)；

(7)、去噪信号为具有最小噪声方差或指定迭代次数情况下的趋势x₁和所有去噪后分量(c＝2,…,d)之和。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明利用混沌理论中量化***复杂度的模糊熵来表征噪声平面，模糊熵能够定量地表示每个奇异谱分量相对于白噪声和原始噪声信号的噪声水平，根据模糊熵谱和SSA的滤波器特性，通过四个合成信号和两个实验信号的仿真表明，迭代奇异谱软阈值去噪方法的去噪性能明显优于传统的截断奇异谱方法，稍优于小波变换和经验模态分解去噪方法。

2、无论信号性质和噪声水平，无论传统的“噪声平面”是否存在，本发明提出的模糊熵谱均能准确提供含噪信号各分量的相对噪声水平，这为不同信号去噪或改进其他信号去噪方法提供了重要基础。本发明可以广泛应用于移动装备、助听器、可穿戴装备、医疗器械或生物医学信号、机械信号、雷达信号的去噪。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是SNR分别为0dB(左)和15dB(右)的piecewise-regular信号的模糊熵谱(实线)及奇异谱(点线)，虚线和点划线分别为无噪及含噪piecewise-regular信号的模糊熵值，矩形为高斯白噪声模糊熵的95％置信带；

图2是SNR分别为0dB(左)和15dB(右)的Riemann信号的模糊熵谱(实线)及奇异谱(点线)，虚线和点划线分别为无噪及含噪Riemann信号的模糊熵值，矩形为高斯白噪声模糊熵的95％置信带；

图3是本发明的原理图；

图4(a)一含噪piecewise-regular信号样本；

图4(b)至图4(d)分别是SSA-IST算法迭代一次(b)，四次(c)，和14次(d)的去噪信号；

图4(e)是SNR改进曲线；

图5是SSA-IST(实线)，截断SSA(虚线)，WT(点划线)和EMD(点线)四种算法滤波piecewise-regular(a)，Riemann(b)，blocks(c)，和sineoneoverx(d)信号后的信噪比。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书(包括任何附加权利要求、摘要和附图)中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

模糊熵是混沌理论中用来表征***复杂度的一个混沌不变量，在此，我们首先提出模糊熵谱的概念并利用该谱来获得真实意义上的噪声平面。无论平坦的奇异谱是否存在，模糊熵谱都能够揭示每个信号分量的相对噪声水平并鉴别该分量是信号还是噪声占优。

当SSA分解一信号为其多组成分量时，基于SSA分量的模糊熵谱特性，我们发明了一种奇异谱迭代软阈值(SSA-IST)去噪方法：为了能够滤除每个分量中的噪声，根据这些SSA分量的模糊熵谱，我们将其划分为信号占优或噪声占优的两组，我们提出两个不同的阈值公式来分别过滤信号或噪声占优的分量。我们通过四个合成信号和两个实验信号(语音和肌电信号)的去噪实验来表明该方法的有效性及相对截断奇异谱方法的显著性能改善。

I.SSA和信号噪声模型

考虑长度N的含噪信号x_in＝{x₁,x₂,…,x_N}，并假设其中的加性白噪声与信号不相关，即：

这里代表无噪信号，n为白噪声。

选择合适窗长度d，x_in可通过SSA的第一步嵌入变换为一轨迹矩阵：

这里m＝N-d+1是多维延迟矢量的数目，H为一Hankel矩阵，即其沿主对角i+j＝const上的元素相等。H可表示为：

这里和N分别代表方程(1)中和n构造的Hankel矩阵。

假设m≥d，轨迹矩阵可通过奇异值分解(SVD)为：

H＝UΣV^T,U≡(u₁,u₂,…,u_m),V≡(v₁,v₂,…,v_d) (4)

这里和分别由均正交的左右奇异矢量u和v构成，Σ＝diag(λ₁≥λ₂≥…≥λ_d≥0).对角元素Σ称为H的奇异值，其集合即奇异谱。根据方程(3)，H的SVD也可表达为：

这里k为满足λ_k＞＞λ_k+1的截止参数。奇异值(λ_k+1,λ_k+2…λ_d)构成了所谓的“噪声平面”。真实信号的Hankel矩阵也可通过SVD表达为：

为了能从x_in中恢复截断SSA算法将噪声平面上的奇异值置零，即(λ_k+1＝λ_k+2＝…＝λ_d＝0)，这等价于获得从下式获得信号估计

最后，通过对约减秩的矩阵上反对角上元素的平均运算，我们可以获得去噪信号

从方程(6)到(8)我们可以发现截断SSA算法移除了噪声空间，但是保留了噪声在信号子空间的投影，这样的算法实际上包含了可能性最高的残存噪声。

II.模糊熵辅助的SSA阈值去噪方法

A.模糊熵

模糊熵是一个鲁棒的量化信号复杂度的测度，它适用于任何非线性非平稳信号。我们重新表达方程(2)中H的第e个矢量序列：

矢量对和之间的距离定义为：

矢量对和之间的相似度通过高斯型的模糊隶属度函数来确定：

这里r为模糊函数边界宽度，给定r，矢量与其所有邻点之间的平均相似度为：

利用模糊数学中模糊概率的思想，矩阵H中所有矢量对相似的概率定义为：

类似地，我们可以利用原始信号x_in构造一系列(d+1)维矢量(1≤e≤N-d)，并且定义相似度以及概率模糊熵则定义为：

这里d和r分别设置为2和0.2。

任何无噪或去噪后的带限信号具有较小的熵值，而含噪信号具有较大熵值，且随噪声水平升高，熵值变大，我们用5000个样本的蒙特卡洛模拟表明，高斯分布和均匀分布白噪声的模糊熵95％置信带的上下限分别为1.66和1.74。与其他多个量化信号复杂度的线性、非线性统计量相比，模糊熵具有更好的单调性、相对一致性，对任何一组通过SSA、WT、EMD或任何其他信号分解方法获得的分量，我们定义这些分量模糊熵的分布即为模糊熵谱模。

B.SSA模糊熵谱

利用SSA分解噪声信号，第一个分量为低频分量，代表了原始信号的主要趋势，而剩余分量则具有明显的震荡和高频特性。

图1为piecewise-regular信号在信噪比0dB(左)和15dB(右)下的模糊熵谱(实线)和奇异谱(点线)，无噪(虚线)及含噪(点划线)信号的模糊熵以及高斯白噪声模糊熵的95％置信带(矩形)也置于图形中供比较。图2提供了Riemann信号的同样信息，两图中的模糊熵排列对应于奇异值从大到小的顺序。

分析图1和图2，我们可以得出如下结论：第一，很多情况下，并无所谓的奇异谱噪声平面存在，比如我们考察的信号中，piecewise-regular信号在两个信噪比下，Riemann信号在15dB下均无平坦的噪声平面存在，但我们提出的模糊熵谱则清晰地勾画出了每个分量的相对噪声水平，提供了另一个“噪声平面”来有效地表征每个分量相对于白噪声、干净信号和含噪信号的噪声水平；第二，在每种情况下，第一个SSA分量的模糊熵均与干净信号的模糊熵非常接近；最后，也是最重要的发现是奇异谱对应的高频分量并不都在白噪声95％置信带内，说明这里高频分量中包含信号成分，也表明截断SSA剔除所有高频分量会导致信息损失，影响去噪性能。

C.迭代SSA软阈值去噪方法

图3为迭代SSA软阈值去噪算法原理图，根据前述对噪声信号模糊熵谱的考察，模糊熵值大于原始噪声信号特别是那些逼近或在白噪声95％置信带内的SSA分量，需要较大的阈值来滤波，另一方面，模糊熵值小于噪声信号的分量具有更多的信号成分，需要较小的阈值，特别是第一个分量，包含了信号的主要趋势，并且其模糊熵值与干净信号非常接近，因此第一分量应完全保留在信号恢复中，据此，迭代SSA软阈值去噪算法归纳如下：

1.SSA分解：噪声信号x_in＝{x₁,x₂,…,x_N}嵌入为m×d的Hankel矩阵，SVD分解Hankel矩阵为d个秩1矩阵之和，并重构为d个分量(c＝1,…d)；

2.计算模糊熵谱：按模糊熵谱公式(14)计算SSA分量的模糊熵谱和原始噪声信号的模糊熵值；

3.设置阈值：第一分量x₁不做任何阈值去噪，

对模糊熵值小于噪声信号模糊熵值的分量x₂到x_k，设置较小阈值：

这里σ_c为分量x_c方差，N为信号长度，

剩余分量的模糊熵值均不小于噪声信号，设置较大阈值：

4.去噪：除第一分量的所有分量x_c(c＝2,…,d)，利用相应的阈值公式(15)或(16)进行去噪：

这里j＝1,2,…,N，去噪后分量为(c＝2,…,d)；

5.信号恢复：第一次去噪信号为第一分量和其余软阈值去噪分量之和：

6.迭代：通过模糊熵谱一次去噪的信号仍可能有残留噪声，我们通过迭代机制进一步提高信号信噪比。估计噪声的方差，并用估计信号重构嵌入矩阵，重复步骤(1)到(5)；

7.迭代终止：比较连续迭代获得的噪声方差，如噪声方差不再明显减小或达到指定迭代次数，迭代停止，否则重复(1)到(6)，去噪信号为具有最小噪声方差或指定迭代次数情况下的趋势和所有去噪后分量之和。

III.性能评估

除了piecewise-regular和Riemann信号，我们还利用另外两个代表性的合成信号，即blocks和sineoneoverx信号来表明我们去噪算法的性能，两个真实的实验信号，即语言和肌电信号也被用来评估该算法，并和截断SSA，WT，及EMD去噪算法进行比较，四个合成信号样本均采用了四个不同采样频率，生成长度分别为1024，2048，4096，和8192的样本，性能参数为50个噪声样本去噪后信噪比的均值。

A.迭代次数效果

图4以piecewise-regular信号为例，说明了SSA-IST算法中迭代机制的效果，图4(a)为一长度8192信噪比15的piecewise-regular信号样本，图4(b)为SSA-IST去噪一次后信号波形，图4(c)和(d)分别为SSA-IST算法迭代4和14次的输出，从这些时域波形上我们可以看出随迭代次数的增加，噪声被显著压制，图4(e)定量显示了SNR随迭代次数提高的效果，通过14次的迭代，含噪信号SNR从15dB提高到约25dB，这充分表明了SSA-IST算法去噪的有效性。

B.滤波合成信号的性能

图5显示了SSA-IST方法，以及截断SSA，WT和EMD四种去噪方法对四个合成信号的去噪效果，每个子图上的曲线对应于不同样本长度信号去噪后的SNR，按去噪前SNR它们被分左侧的0dB和右侧的15dB，可以看出在所有情况下SSA-IST(实线)去噪后信噪比都明显高于截断SSA(虚线)2到5个dB，另外，在大多数情况下，SSA-IST信噪比改善稍好于或非常接近于EMD(点线)，唯一例外是在sineoneoverx信号长度8192情况下，EMD优于SSA-IST，最后，SSA-IST与WT(点划线)相比，在低信噪比的情况下，SSA-IST去噪的优势更为明显，当SNR为15dB时，WT趋向于和SSA-IST性能相似或稍优，这种趋势在sineoneoverx信号中表现尤为明显。

C.滤波实验信号的性能

我们还评估了SSA-IST对两种真实实验信号，即语音和肌电信号的去噪性能，语音信号的采样频率为16kHz，肌电信号的采样频率为1kHz。表1为SSA-IST和其他三种方法对两种实验信号去噪后的信噪比总结，多个样本去噪后信噪比的方差是另一个，特别是针对真实实验信号有效的去噪性能评估方法，四种方法去噪后SNRs的方差也同时列于表1。与合成信号类似，在所有情况下SSA-IST去噪后的信噪比均高于截断SSA法，此外，除了语音信号15dB，SSA-IST去噪后信噪比在其他情况下比WT均高出1–2dB，与SSA-IST相比，EMD是较具竞争力的方法，在不同输入信号信噪比下，EMD对语音信号去噪的信噪比改善都优于SSA-IST，但是对肌电信号，SSA-IST在不同信噪比下都优于EMD，并且在两种信号情况下，SSA-IST去噪SNRs的方差均小于EMD法。

在所测试的六个信号中，几种方法对肌电信号去噪效果最差，其次是Riemann信号，这可能与这两种信号的复杂成分相关，EMG信号中包含了很多棘波，而Riemann信号具有明显的高频成分，波形与1/f噪声极为相似。

表1.四种方法对语音和肌电信号去噪后的信噪比(SNRs)以及信噪比的方差

IV.总结

传统的截断SSA滤波方法采用了保留部分分量，舍弃其他分量的二元方法，在频域等价于低通滤波，损失了高频带的信息，另外，该方法依赖于主观查找很多情况下并不存在的噪声平面，本发明首先提出了用模糊熵谱来替代奇异谱，无论信号性质和噪声水平，模糊熵谱均能准确提供含噪信号各分量的相对噪声水平，这为发明不同信号去噪或其他处理方法提供了重要基础。在模糊熵谱的基础上，我们发明了迭代奇异谱软阈值去噪算法，通过对四种仿真信号和两个实验信号的在不同信噪比下去噪的验证，我们发明的SSA-IST均比截断SSA更能有效地提高含噪信号信噪比，此发明可以广泛应用于移动装备，助听器，可穿戴装备，医疗器械或生物医学信号，机械信号，雷达信号的去噪。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (2)

1.一种基于模糊熵的含噪信号处理方法，其特征在于：对于含噪信号，假设长度N的含噪信号x_in＝{x₁,x₂,…,x_N}，并假设其中的加性白噪声与信号不相关，即：这里代表无噪信号，n为白噪声，利用原始信号x_in构造一系列d维矢量(1≤e≤N-d+1)，并且定义相似度以及模糊概率同样方法构造(d+1)维矢量及其相应的相似度和模糊概率模糊熵则定义为：这里d和r分别设置为2和0.2，N为待去噪信号长度，对利用任何信号分解方法获得的分量，其所有分量的奇异谱分布定义为模糊熵谱。

2.一种基于模糊熵的迭代奇异谱软阈值去噪方法，其特征在于，包括如下步骤：

(1)、SSA分解：噪声信号x_in＝{x₁,x₂,…,x_N}嵌入为m×d的Hankel矩阵，SVD分解Hankel矩阵为d个秩1矩阵之和，并重构为d个信号分量

(2)、计算模糊熵谱：按权利要求1所述模糊熵谱定义计算SSA分量的模糊熵谱和原始噪声信号的模糊熵值；

(3)、设置软阈值：第一分量x₁不做任何阈值去噪，对模糊熵值小于噪声信号模糊熵值的分量x₂到x_k，设置较小阈值对剩余的x_k+1到x_d分量，设置较大阈值这里σ_c为分量x_c的方差；

(4)、软阈值去噪：除第一分量的所有分量x_c(c＝2,…,d)，利用相应的阈值进行去噪，即每一分量的某噪声信号数值绝对值如小于该分量阈值，输出为零，如信号数值大于阈值，则输出为信号数值减去阈值，剩余情况输出为信号数值加上阈值，软阈值去噪后所有分量与x₁之和为第一次估计信号估计噪声为

(5)、估计的方差，用估计信号作为输入信号，重复步骤(1)到(4)；

(6)、比较连续迭代获得的噪声方差，如噪声方差不再明显减小或达到指定迭代次数，迭代停止，否则重复(1)到(5)；

(7)、去噪信号为具有最小噪声方差或指定迭代次数情况下的趋势x₁和所有去噪后分量之和。