一、EMD降噪1.1 EMD的基本原理经验模态分解Empirical Mode Decomposition, EMD是将复杂信号分解为有限个固有模态函数Intrinsic Mode Function, IMF和残差信号的方法如图1所示。图1 EMD与IMF各IMF分量包含了原始信号的不同时间尺度的局部特征信号IMF信号仅包含一种振荡模式即单一的瞬时频率IMF需要满足以下两个条件1极值点和过零点的数目应相等或最多相差一个2局部最大值和局部最小值的上下包络线均值为零。EMD分解的步骤如下1包络线获取确定原始信号的极大值序列和极小值序列采用三次样条曲线对极值点进行拟合形成上下包络线、此处常用希尔伯特变换Hilbert Transform, HT求取包络线计算上下包络线的平均值2残余信号获取原始信号与平均包络信号做差得到残余信号3IMF条件检验检验是否满足IMF条件当满足时即获得一个固有模态函数若不满足则重复步骤1和24残差信号获取原始信号与平均包络信号做差得到残余信号将残差信号作为原始信号重复以上过程不断得到各固有模态函数5分解终止标志当残差信号为单调函数时终止分解此时原始信号可表示为1.2 EMD降噪的实现过程根据EMD分解过程不难看出在进行信号处理时只需要对信号分解的IMFs和残差信号进行约束然后采用分解得到的IMFs对信号重构即可获得降噪的信号。基于EMD的降噪过程包括以下步骤1EMD分解得到一系列IMFs和残差信号2IMFs筛选筛选包含有用信息的IMFs3信号重构将筛选的IMFs叠加获得降噪后的信号。不难发现基于EMD的降噪过程中最关键的是IMFs的筛选目前常用的方法包括频谱筛选法、峭度/峰值因子筛选法、相关系数筛选法和自适应法该部分暂不进行展开等。频谱筛选法计算各IMF的频谱根据需要选取频段筛选IMFs但这种方法基于研究人员清晰知道自己所需的频段峭度/峰值因子筛选法计算各IMF的峭度值或峰值因子一般认为峭度值越大信号中包含较多的毛刺该类信号通常与噪声相关重构时进行去除相关系数筛选法计算IMFs和原始信号的相关系数设置合理阈值筛选需要的IMF重构信号。一般地相关系数介于-1和1之间相关系数绝对值越大信号相关性越强。1.3 EMD的不足采用EMD分解的方法一般存在末端效应和模态混叠的问题。末端效应发生在信号的头部或者尾部通常由于拟合包络线时末端发散导致进而扭曲IMF末端的波形模态混叠在同一个IMF分量中,存在尺度分布范围很宽却又各不相同的信号在不同的IMF分量中,存在着尺度相近的信号。模态混叠主要是因为在求取包络线的过程中信号中存在异常事件影响极值点的选取使极值点分布不均匀从而导致求取的包络为异常事件的局部包络和真实信号包络的组合。经该包络计算出的均值再筛选出的IMF分量就包含了信号的固有模式和异常事件或者包含了相邻特征时间尺度的固有模式从而产生了模态混叠现象。1.4 EMD降噪的python实现​​​​​​​# todo: 1.定义函数对信号进行EMD去噪。def denoised_with_EMD(data): 基于EMD分解峭度准则筛选IMF分量重构得到去噪信号 筛选规则保留超额峭度小于0的IMF分量进行信号重构 Parameters ---------- data : np.ndarray, 1D 原始一维振动时域信号 Returns ------- imf_select : list 满足峭度条件的IMF分量集合 order_select : list 被选中的IMF分量序号从1开始计数 signal_denoised : np.ndarray IMF分量叠加重构后的去噪信号 # 1.定义变量存储筛选后符合条件的IMF分量 imf_select [] # 2.定义变量存储被选中IMF对应的序号 order_select [] # 3.EMD分解 emd EMD() imfs emd.emd(data) # imfs: (n_imf, n_samples), 最后一行为残差趋势项 # 4.遍历所有分量【不包含最后残差项】(len(imfs)-1) for i in range(len(imfs)-1): # 4.1 计算当前IMF的超额峭度并打印 print(kurtosis(imfs[i])) # 4.2 筛选条件超额峭度 0扁峰类平缓周期分量 if kurtosis(imfs[i]) 0: imf_select.append(imfs[i]) order_select.append(i1) else: # 峭度0尖峰分量冲击/噪声直接舍弃 continue # imf_select[:len(imf_select)]等价于imf_select此处简化为imf_select signal_denoised np.sum(imf_select[:len(imf_select)], axis0) return imf_select, order_select, signal_denoised二、VMD降噪2.1 VMD的基本原理变分模态分解Variational Mode Decomposition, VMD假设信号由一系列具有特定中心频率、有限带宽的子信号IMFs组成如图2所示。因此VMD问题最终即是求解中心频率与带宽限制的问题其在经典wiener滤波的基础上通过变分问题进行求解最终得到模态函数。VMD的数学原理和推导过程较为复杂本文不进行详细展开具体过程可参考文献[1]。图2 VMD与IMF实际上变分问题最终是求解泛函的极值问题最终可以转换成优化求解的问题VMD过程可简单描述如下1变分问题构造/优化模型构造以下从优化模型的角度对变分问题进行简单介绍假设信号可分解为K个有限带宽的模态分量组成每个IMF中心频率为则优化模型可表示如下目标函数通过希尔伯特变换并采用算子调制得到基带计算解调梯度的平方范数并估计每个模分量的带宽最终得到如下目标函数约束条件模态和等于输入信号表示如下设计变量模态分量每个IMF中心频率为。2变分问题求解以上优化模型的求解可采用诸多优化求解器包括二次惩罚项、拉格朗日算子、增广拉格朗日函数等各求解过程较为复杂此处不展开讨论感兴趣可参考文献[2]。在使用VMD分解时需要设定两个重要参数即分解层数K和惩罚因子α优化求解时引入的二阶惩罚系数两个分解参数对结果影响较大。当分解层数K小于待分解信号的有效成分的个数时出现欠分解导致模态混叠当分解层数K大于待分解信号中的有用成分的个数时出现过分解产生虚假分量为了获取最佳的K值可使用峭度最大原理、能量差值原则和优化算法进行计算。惩罚系数 α 决定着IMF分量的带宽惩罚系数越小 各 IMF 分量的带宽越大过大的带宽会使得某些分量包含其他分量信号α值越大各IMF分量的带宽越小过小的带宽是使得被分解的信号中某些信号丢失该系数常见取值范围为1000~3000。2.2 VMD降噪的实现过程与EMD降噪过程类似VMD降噪过程如下1VMD分解得到各模态分量IMF2IMFs筛选筛选包含有用信息的IMFs3信号重构将筛选的IMFs叠加获得降噪后的信号。2.3 VMD的优点1可自定义模态个数2通过VMD方法分解出来的IMF都具有独立的中心频率并且在频域上表现出稀疏性的特征具备稀疏研究的特质3在对IMF求解过程中通过镜像延拓的方式避免了类似EMD分解中出现的端点效应4通过合理的K值选取可避免模态混叠。2.4 VMD降噪的python实现# todo: 2.定义函数对信号进行VMD分解 def denoised_with_VMD(data): 基于VMD变分模态分解峭度准则筛选模态分量重构得到去噪信号 筛选规则保留超额峭度小于0的模态分量用于信号重构 Parameters ---------- data : np.ndarray, 1D 原始一维振动时域信号 Returns ------- imf_select : list 满足峭度条件的VMD模态分量集合 order_select : list 被选中的VMD模态分量序号从1开始计数 signal_denoised : np.ndarray 筛选后模态分量叠加重构的去噪信号 # 1.定义变量存储筛选后符合峭度条件的VMD模态分量 imf_select [] # 2.定义变量存储被选中模态对应的序号 order_select [] # 3.VMD分解,alpha为惩罚因子K为分解层数tol为优化终止条件 u, u_hat, omega VMD(data, alpha2000, tau0.0, K4, DC0, init1, tol1e-7) # 4.遍历所有VMD模态分量u为二维数组每行是一个模态分量 for i in range(len(u)): # 4.1 筛选条件超额峭度 0保留扁峰平缓周期分量 if kurtosis(u[i]) 0: imf_select.append(u[i]) order_select.append(i1) # 4.2 峭度≥0代表尖峰分量冲击/噪声舍弃该模态 else: continue # 5.将选中的模态分量叠加求和重构信号 # 6.imf_select[:len(imf_select)]为冗余切片简化为imf_select signal_denoised np.sum(imf_select[:len(imf_select)], axis0) return imf_select, order_select, signal_denoised# todo: 3.测试代码 if __name__ __main__: # 1.生成数据 n 500 # 生成500个点的信号 t np.linspace(0, 30*np.pi, n, endpointFalse) s np.cos(0.1 * np.pi * t) np.sin(0.3 * np.pi * t) np.cos(0.5 * np.pi * t) np.sin(0.7 * np.pi * t) # 原始信号 r 0.5 * np.random.randn(n) # VMD噪声 # r 0.1 * np.random.randn(n) # EMD噪声 y s r # 加噪声 imf_select, order_select, signal_denoised denoised_with_VMD(y) n len(imf_select) 2 fig plt.figure(figsize(16, 12)) plt.subplots_adjust(hspace0.5) font {family: Times New Roman, size: 16, weight: normal, } matplotlib.rc(font, **font) plt.subplot(n, 1, 1) plt.plot(t, y, linewidth1.5, colorb) plt.title(原始模拟信号, fontnamemicrosoft yahei, fontsize16) for i in range(2, n): plt.subplot(n, 1, i) plt.plot(t, imf_select[i - 2], linewidth1.5, colorr) plt.title(IMF%d % order_select[i - 2], fontnamemicrosoft yahei, fontsize16) plt.subplot(n, 1, n) plt.plot(t, signal_denoised, linewidth1.5, colorb) plt.title(降噪信号, fontnamemicrosoft yahei, fontsize16) plt.show()EMD降噪信号VMD降噪信号