CN110096998A - 一种用于机加工表面形貌特征提取的分解方法 - Google Patents

Abstract

本发明公开一种用于机加工表面形貌特征提取的分解方法，包括以下步骤：(1)利用粗糙度轮廓仪测量出零件难加工位置的表面数据信号；(2)利用小波包分解将上述表面数据信号进行阈值降噪处理并获得重构信号(3)对现存的自适应噪声的完整集合经验模态分解(CEEMDAN)进行预处理改进，通过选取合理的白噪声上限频率以减少白噪声幅值和集合次数，利用改进的CEEMDAN对上述获得的重构信号进行分解，得到一系列的固有模态函数IMF；(4)对上述获得的IMF进行相关性分析获得频率的特征，剔除相关性较差的IMF。

Description

一种用于机加工表面形貌特征提取的分解方法

技术领域

本发明涉及机加工中加工误差的几何测量与评定，特别是涉及一种用于机加工表面形貌特征提取的分解方法。

背景技术

机械加工过程中会产生振动、环境噪声以及人工测量误差等干扰，这些干扰会使得采集到的表面形貌信号发生调频和调幅。因此，能否从机械振动、背景噪声和测量误差中准确的提取出真实信号，并确定真实信号的频率特征，对确定加工表面误差和预测表面性能具有重要的意义[1]。在实际的机械加工过程中，产生的信号大多是非平稳和非线性的，而传统的时域和频域分析方法都是假设信号是平稳的和线性的，故而不再适用。

经验模态分解(EMD)是由纯数据驱动的一种非线性非平稳信号的自适应分解方法，不需要像小波变换那样提前预置基函数，但是会因为噪声和不连续信号的存在产生模态混叠的问题。模态混叠最早由Huang通过对含有间断信号的EMD分解时发现。模态混叠是指在一个固有模态函数(IMF)中包含差异极大的时间尺度；或者相近的特征时间尺度分布在不同的IMF中，表现为相邻两个IMF波形混叠，相互影响，难以分辨[2]。针对EMD中存在的模态混叠问题，各国学者纷纷对其进行研究，Wu和Huang[3]提出了集合经验模态分解(EEMD)，是一种噪声辅助数据分析(NADA)方法。EEMD的原理[4]就是：将白噪声和不同尺度的信号组分均匀的填充在整个时频空间，不同尺度的信号会自动投射在由白噪声建立的参考尺度上，使得EMD变成了自适应的二元滤波器组。EEMD在很大程度上减轻了模态混叠现象，但并没有完全消除，而且又出现了新的问题[5]。重构后的信号包含了残余噪声；不同的噪声与信号集合产生不同数量的模态，使得最后的平均化比较困难；分解过程中需要大量的迭代次数，时间成本比较高。文献[5]中提出了自适应噪声的完整集合经验模态分解(CEEMDAN)是对EEMD的重要改进，实现了可忽略的重建误差并且解决了针对信号加噪声的不同集合产生不同模态的问题。尽管如此，CEEMDAN仍然有一些方面值得改进：模态中包含了一些残余噪声；在分解的早期阶段，信号信息比在EEMD中显示的晚，并伴有一些虚假模态；分解过程中的迭代次数依旧很多，计算成本较高[6]。

针对现存的分解方法很有必要设计一种将小波包[7]和经验模态分解相结合的方法来更加精确的对机加工表面形貌进行特征分解与提取。

参考文献：

[1]FP.Zhang,JB.Yang,TG.Zhang,Y.Yan，Surface topography separationbased on wavelet reconstruction and empirical mode decomposition,10^th ASMEInternational Manufacturing Science and Engineering Conference.2015.WOS:000365146800027.

[2]N.E.Huang,Z.H.Wu,Areview on Hilbert–Huang transform:method and itsapplications to geophysical studies,Rev.Geophys.46(2)(2008).

[3]Z.H.Wu,N.E.Huang,Ensemble empirical mode decomposition:a noise-assisted data analysis method,Adv.Adapt.Data Anal.1(01)(2009)1–41.

[4]Wu Z.Ensemble empirical mode decomposition:a noise assisted dataanalysis method[J].Adv Adapt Data Anal,2009,1.

[5]E.Torres,Marcelo A.Colominas,Gasto′n Schlotthauer,et al.ACOMPLETE ENSEMBLE EMPIRICAL MODE DECOMPOSITION WITH ADAPTIVE NOISE[C].Proceedings of the IEEE International Conference on Acoustics,Speech,andSignal Processing,ICASSP 2011,May 22-27,2011,Prague Congress Center,Prague,Czech Republic.IEEE,2011.

[6]Colominas M A,Schlotthauer,Gastón,Torres,María E.Improved completeensemble EMD:A suitable tool for biomedical signal processing[J].BiomedicalSignal Processing and Control,2014,14:19-29.

[7]Hao Q,Zhang X,Wang Y,et al.A novel rail defect detection methodbased on undecimated lifting wavelet packet transform and Shannon entropy-improved adaptive line enhancer[J].Journal of Sound and Vibration,2018,425:208-220.

发明内容

本发明的目的是为了克服现有技术中的不足，提供一种用于机加工表面形貌特征提取的分解方法，该方法将小波包分解和改进的经验模态分解相结合对获得的信号进行降噪重构处理，再进行表面形貌特征分解与提取，为追溯表面性能和评估表面误差做准备。

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

一种用于机加工表面形貌特征提取的分解方法，包括以下步骤：

(1)利用粗糙度轮廓仪测量出零件难加工位置的表面数据信号；

(2)利用小波包分解将上述表面数据信号进行阈值降噪处理并获得重构信号，具体如下：

(201)在小波库中选取不同的正交小波基，得到小波包系数序列{u_k}，同时给定代价函数M，计算不同小波基的代价函数，选取最优小波包树；

(202)通过半硬阈值判断方法对上述获得的最优小波包树分解系数进行阈值判断，剔除噪声信息获得重构信号。

(3)对现存的自适应噪声的完整集合经验模态分解(CEEMDAN)进行预处理改进，利用改进版CEEMDAN对上述获得的重构信号进行分解，得到一系列的固有模态函数IMF；

(4)对上述获得的IMF进行相关性分析获得频率的特征，剔除相关性系数小于0.95的IMF。

进一步的，步骤(202)中半硬阈值的判断函数如下：

式中w_j,k表示原小波系数，表示重构后的估计小波系数，th为阈值，j,k分别为分解尺度和数据长度。

进一步的，步骤(3)中预处理改进步骤如下：

a.给定原始信号x(t)以及采样频率fs，预测上限频率fu_i＝(i+1)·fs,i＝1,2...,n，n为15-20，i+1代表每次频率放大倍数，n代表最大的放大倍数；

b.利用fu_i重新采样并利用三次样条插值法获得插值信号x′_i(t)；

c.对上述获得的插值信号进行EEMD分解，添加白噪声的上限频率为fu_iHz，白噪声幅度为0.2std(x)，std为标准差，集合实验次数为100；

d.计算每个IMF和原始信号的相关系数，计算不同上限频率的相对均方根误差RRMSE值，t为时间，N为最大的时间长度，imf_max为与原始信号相关系数最大的固有模态函数；

e.选取上述最大的RRMSE值所对应的频率为添加白噪声的上限频率fp，得到带限白噪声，将带限白噪声添加到所对应插值信号x′(t)中进行CEEMDAN。

与现有技术相比，本发明的技术方案所带来的有益效果是：

本发明方法前期选取最优小波包树进行半硬阈值判断的降噪处理并重构，可以有效的剔除原始信号中存在的背景噪声，后续的改进版CEEMDAN经过预处理改进，加入最优的带限白噪声，有效的减小所添加白噪声的幅值和集合实验次数，减少运行成本提高计算效率，而且还使得表面形貌中的不同频率特征能够被很好的提取出来，避免了残余噪声和虚假模态的产生，与此同时该方法还具有很好的鲁棒性和通用性。

附图说明

图1是本发明方法的流程示意图。

图2是所获得模拟信号的组成成分。

图3最优小波包树搜索流程图。

图4(a)为硬阈值降噪效果图；图4(b)为软阈值降噪效果图；图4(c)为半硬阈值降噪效果图。

图5白噪声上限频率的选取结果图。

图6改进CEEMDAN算法流程示意图.

图7改进CEEMDAN的分解图。

具体实施方式

以下结合附图和具体实施例对本发明作进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示为本发明的整体流程示意图，具体实施步骤如下：

1、获取表面数据信号

为了验证本发明方法的合理性和正确性，在该步骤选取模拟信号来充当实际的机加工表面信号。所选取的模拟信号如下：

x₂＝sin(2·pi·f₂·t) (2)

x＝x₁+x₂ (3)

上式中x₁为一个间断信号，x₂为一个正弦信号，pi代表圆周率π，在这里信号频率分别设为f₁＝100Hz，f₂＝40Hz，时间t为1秒，并在信号x上加入一些背景随机噪声得到x′，获得的信号如下图2所示。

对于上述获得的模拟数据信号按照图1所示的流程路线进行表面形貌的分解与特征提取。2.小波包分解降噪处理

小波分解是一种多尺度分解的时-频局域变换，可分析包含不同尺度的信号，但在某些情况下高频端频率分辨率低，而小波包分解可以弥补这个缺点，小波包库包含了丰富的小波包基，不同的小波包基具有不同的性质，反映不同的信号特性，因此选取一个合理的最优小波包树是有必要性的。

小波包定义可以表示为如下的二尺度方程：

其中g_k＝(-1)^kc_1-k，即尺度方程的两个系数具有正交性。尤其当n＝0时，式(4)变成：

式(5)中分别为尺度函数u₀(t)与小波函数u₁(t)的二尺度方程。式(4)(5)中c_k为尺度系数，g_k为小波系数，k为分解尺度，l²为空间，t为时间，u_2n和u_2n+1分别为尺度函数和小波函数。

由式(4)和(5)构造的序列{u_n(t)},(n∈Z₊)，称为由基函数确定的正交小波包。由于尺度函数是由c_k唯一确定的，那么就可以将序列{c_k}看成{u_n(t)}_n∈Z的正交小波包。

2.1最优小波包树的选取

给定信号x(t)，使其在一组给定的小波基下展开，得到一系列的小波包系数，在该序列上定义代价函数M，在小波库的所有小波基中寻找使代价函数最小的基，则这个代价最小的小波包基所构成的分解树就是最优小波包树。此处选取的代价函数M为信息熵。

最优小波包树的选取步骤如下，见图3：

a.将分解树的各节点上写上代表信息熵值的数字，并在树的最低层各节点上标上*号；

b.将最低层的信息熵值作为一个初始值，将N层节点熵值和N-1层熵值进行比较，用*号标定低熵节点，如此一直进行到顶层。

c.*号的节点表示该树杈在该节点的代价函数最小，搜索到此为即可停止。搜索从顶层开始，分别经过每一个节点一直到最低层，最后将带*号的节点选出，则经过修剪后的分解树就是最优的小波包树。

2.2阈值降噪

传统的阈值降噪分为硬阈值降噪和软阈值降噪，通过硬阈值得到的估计小波系数的连续在较差，会引起重构信号的振荡，产生较大的均方差；通过软阈值得到的估计小波系数和原小波系数存在着恒定偏差，使重构信号产生较大的均方差。针对传统的硬阈值和软阈值降噪的缺点与不同，提出一种半硬阈值降噪，阈值函数如下：

在上式中w_j,k表示原小波系数，表示重构后的估计小波系数，th为阈值，j,k分别为分解尺度和数据长度，|w_j,k|≥th时，有随着|w_j,k|的增大，和w_j,k之间的偏差会逐渐减小，使得到的估计小波系数连续性变好，产生较小的误差。上述三种降噪效果图见图4(a)、图4(b)和图4(c)。

3.改进的CEEMDAN

对于噪声信号，极值密度(每单位长度的极值平均数)表示信号的频率，高频噪声信号相比低频噪声信号可以更好的改变原始信号的极值分布。在经验模态分解中，原始信号与相应IMF之间的标准差由下式表示：

在式(7)中，A是添加白噪声的振幅，N是集合的总数，ε是原始信号和相应IMF的标准差。白噪声振幅的减小可以减少噪声残留引起的误差，与此同时，白噪声的频率并非越高越好，过高的频率会引起计算成本的提高，因此要选择合适的上限频率噪声，既可以减少所添加白噪声的振幅，也可以产生带限白噪声提高效率。

3.1白噪声上限频率的确定过程如下：

i.给定原始信号x(t)以及采样频率fs，预测上限频率fu_i＝(i+1)·fs,i＝1,2...,n，n为15-20，i+1代表每次频率放大倍数，n代表最大的放大倍数；

ii.利用fu_i重新采样并利用三次样条插值法获得插值信号x′_i(t)；

iii.对上述获得的插值信号进行EEMD分解，添加白噪声的频率为fu_iHz，噪声幅度为0.2std(x)，std为标准差，集合实验次数为100；

iv.计算每个IMF和原始信号的相关系数，计算不同上限频率的相对均方根误差RRMSE值，imf_max为与原始信号相关系数最大的固有模态函数，t为时间，N为最大的时间长度；

v.选取上述最大的RRMSE值所对应的频率为添加白噪声的上限频率fp。白噪声上限频率的选取结果见图5。

3.2改进CEEMDAN的算法，流程图见图6

i.给定原始信号x[n]，对白噪声信号wⁱ[n]进行预处理，确定其上限频率fp，得到带限白噪声bⁱ[n]；

ii.对信号x[n]按照fp重新采样，并进行三次样条插值得信号x′[n]；

iii.利用EMD对x′[n]+β₀E₁(bⁱ[n]),i＝1,...,I进行分解获得第一个固有模态函数IMF，符号E代表经验模态分解，i代表集合的个数，I代表集合的总个数，代表分解得到的固有模态函数，β为白噪声的幅度；

iv.当k＝1时计算第一个残余项

v.对r₁[n]+β₁E₂(bⁱ[n]),i＝1,...,I进行分解直到产生第一个模态函数，并定义其为原始信号的第二个IMF，

vi.计算第k个残余项，K为IMF的总个数；

vii.对r_k[n]+β_kE_k+1(bⁱ[n]),i＝1,...,I进行分解直到产生第一个模态函数，并定义为第k+1个IMF，

viii.返回步骤iv进行下一个k，最后的残余项为

改进的CEEMDAN通过对添加白噪声的预处理，减少了个别IMF中因添加白噪声而残留的噪声信号，同时也提高了计算效率；在CEEMDAN中，第一个IMF是直接对x[n]+β₀wⁱ[n]提取的，而第二个IMF是对r₁[n]+β₁E₁(wⁱ[n])提取的，这两个模态之间会产生混叠现象和虚假模态，按照改进CEEMDAN的方法便可消除这种情况。图7便是改进CEEMDAN的分解图。

将式(3)的模拟信号分别用EEMD、CEEMDAN、改进的CEEMDAN进行模态分解结果如下表所示。

表1三种经验模态分解的效率比较

	幅度	集合数量	时间(s)	效率提升
					EEMD	0.2std	200	34.63
CEEMDAN	0.2std	200	32.93	4.9％
					改进CEEMDAN	0.01std	100	25.43	22.8％

4.相关性分析

机加工表面形貌包含着大量的与加工过程相关的有用信息，如切削参数、机床刀具的振动等，同时也包含着噪声等无用信息。有用信息与表面形貌的相关性高，无用信息的相关性就低，因此对于表面形貌分解产生的IMF进行相关性分析，剔除相关性低的无用信息。相关系数的计算过程如下式：

在上式中，ρ_i为每个i的相关系数，cov为协方差,Var为方差，x为信号，imf为固有模态函数，K为IMF的总个数。当相关系数大于0.95表明相关性很高，予以保留；小于0.95表明相关性较低，便直接剔除。

本发明并不限于上文描述的实施方式。以上对具体实施方式的描述旨在描述和说明本发明的技术方案，上述的具体实施方式仅仅是示意性的，并不是限制性的。在不脱离本发明宗旨和权利要求所保护的范围情况下，本领域的普通技术人员在本发明的启示下还可做出很多形式的具体变换，这些均属于本发明的保护范围之内。

Claims (3)

1.一种用于机加工表面形貌特征提取的分解方法，其特征在于，包括以下步骤：

(1)利用粗糙度轮廓仪测量出零件难加工位置的表面数据信号；

(2)利用小波包分解将上述表面数据信号进行阈值降噪处理并获得重构信号，具体如下：

(201)在小波库中选取不同的正交小波基，得到小波包系数序列{u_k}，同时给定代价函数M，计算不同小波基的代价函数，选取最优小波包树；

(202)通过半硬阈值判断方法对上述获得的最优小波包树分解系数进行阈值判断，剔除噪声信息获得重构信号。

(3)对现存的自适应噪声的完整集合经验模态分解(CEEMDAN)进行预处理改进，利用改进版CEEMDAN对上述获得的重构信号进行分解，得到一系列的固有模态函数IMF；

(4)对上述获得的IMF进行相关性分析获得频率的特征，剔除相关性系数小于0.95的IMF。

2.根据权利要求1所述一种用于机加工表面形貌特征提取的分解方法，其特征在于，步骤(202)中半硬阈值的判断函数如下：

式中w_j,k表示原小波系数，表示重构后的估计小波系数，th为阈值，j,k分别为分解尺度和数据长度。

3.根据权利要求1所述一种用于机加工表面形貌特征提取的分解方法，其特征在于，步骤(3)中预处理改进步骤如下：

a.给定原始信号x(t)以及采样频率fs，预测上限频率fu_i＝(i+1)·fs,i＝1,2...,n，n为15-20，i+1代表每次频率放大倍数，n代表最大的放大倍数；

b.利用fu_i重新采样并利用三次样条插值法获得插值信号x′_i(t)；

c.对上述获得的插值信号进行EEMD分解，添加白噪声的上限频率为fu_iHz，白噪声幅度为0.2std(x)，std为标准差，集合实验次数为100；

d.计算每个IMF和原始信号的相关系数，计算不同上限频率的相对均方根误差RRMSE值，t为时间，N为最大的时间长度，imf_max为与原始信号相关系数最大的固有模态函数；

e.选取上述最大的RRMSE值所对应的频率为添加白噪声的上限频率fp，得到带限白噪声，将带限白噪声添加到所对应插值信号x′(t)中进行CEEMDAN。