基于自适应本征维数估计流形学习的相空间重构降噪方法

马婧华，汤宝平，宋　涛

(重庆大学机械传动国家重点实验室，重庆　400044)

第一作者马婧华女，博士生，1989年生

基于自适应本征维数估计流形学习的相空间重构降噪方法

马婧华，汤宝平，宋涛

(重庆大学机械传动国家重点实验室，重庆400044)

摘要：针对实际工程领域振动信号噪声干扰大、具有强烈非线性等问题，提出了基于自适应本征维数估计流形学习的相空间重构降噪方法。利用相空间重构将一维含噪时间序列重构到高维相空间；基于极大似然估计法 (maximum likelihood estimate, MLE) 估计相空间中每个样本点的本征维数并使用自适应加权平均法计算全局本征维数；采用局部切空间排列 (Local tangent space Alignment, LTSA) 流形学习方法将含噪信号从高维相空间投影到有用信号的本征维空间中，剔除分布在高维空间中的噪声后，重构回一维时间序列。通过Lorenz仿真实验和风电机组振动信号降噪实例，证实了该方法具有良好的非线性降噪性能。

关键词：非线性降噪；流形学习；本征维数估计；极大似然估计；自适应加权

基金项目：国家自然科学基金资助项目(51275546，51375514);高等学校博士学科点专项科研基金资助(20130191130001)

收稿日期：2014-03-24修改稿收到日期：2014-05-29

通信作者汤宝平男，博士，教授，博士生导师，1971年生

中图分类号:TN911;TH133

文献标志码：A

DOI:10.13465/j.cnki.jvs.2015.11.006

Abstract:Aiming at the problem that the actual engineering vibration signals are interfered by strong noise with strong nonlinear characteristic, a phase space reconstruction method based on adaptive intrinsic dimension estimation manifold learning was proposed. Firstly, one-dimensional time series containing noise were reconstructed into a high dimensional phase space with the phase space reconstruction method. Secondly, the intrinsic dimension of each sample point in the phase space was estimated based on the maximum likelihood estimate (MLE), the adaptive weighted average method was used to calculate the global intrinsic dimension. At last, the manifold learning algorithm and the local tangent space alignment ( LTSA) were employed to project the signal containing noise from the high-dimensional phase space into the intrinsic dimensional space of useful signals. After eliminating the noise distributing in the high-dimensional space, the signals were reconstructed back into one-dimensional time series. Lorenz simulation and an example of vibration signals’ noise reduction for a wind power generator unit showed that the proposed method has a good performance of nonlinear noise reduction.

Phase space reconstruction method for vibration signal’s noise reduction based on adaptive intrinsic dimension estimation manifold learning

MAJing-hua,TANGBao-ping,SONGTao(The State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing 400044, China)

Key words:nonlinear noise reduction; manifold learning; intrinsic dimension estimation; maximum likelihood estimate; adaptive weighted

随着科技的进步，航空发动机、风电机组等机械设备的故障机理日益复杂，其振动信号具有噪声干扰大、非平稳非线性等特征。良好的降噪效果是保证机械设备故障诊断准确性的重要前提[1]。

传统的线性滤波方法不能有效消除非线性信号中的噪声；而非线性小波降噪方法受到参数取值困难的限制[2]；Takens等[3]提出了相空间重构的延时坐标法，奠定了相空间重构技术的基础，为非线性信号分析提供了新的思路。近年来，结合相空间重构和流形学习进行非线性降噪[4]取得了良好的效果，如张振跃等[5]和Chang等[6]在主成分分析的基础上，引入流形学习的局部平滑思想，提出了加权主成分分析局部平滑算法和鲁棒加权局部线性平滑嵌入算法；杨建宏等[7]研究了基于局部切空间排列(LTSA)的非线性时间序列降噪方法；王广斌等[8]提出分形维与LTSA结合的非线性降噪方法，成功应用于带有断齿的齿轮箱振动分析；余成义[9]提出基于相空间重构与LTSA的降噪方法。

LTSA算法[10]直接给出了从低维流形到高维流形之间的主流形重构的显式表达，且保证了映射变化误差最小，但存在本征维数人为设置的不确定性。因此，本文提出了基于自适应本征维数估计流形学习的相空间重构降噪方法。首先通过基于最大似然估计的自适应加权法得到本征维数，然后将重构得到的高维相空间数据约简为低维局部切空间数据，最后采用均值重构得到一维时间序列。

1自适应加权本征维数估计

在无信息丢失的前提下，能够表示数据自由变量的最少变量个数，称为数据集合的本征维数，它是数据集的固有属性。信号的本征维数决定着相空间内相点邻域内的分布性质，因此，可以通过分析相点邻域所张成的子空间来确定。在重构相空间中，有用信号以吸引子形式存在，实质上是分布在以本征维数为大小的特征空间里，而噪声则分布在整个高维相空间中。因此，在使用流形学习进行维数约简时，如果约简目标维数选择过大，则容易将噪声分布方向的数据包含进去，难以很好地消除噪声；而约简目标维数选择过小，则可能将有用信号当作噪声消除掉。因此，本征维数的估计至关重要。

1.1极大似然法本征维数估计

一般来说，数据间的关系可以充分的反映数据的局部几何特征。极大似然估计法[11]就是通过建立近邻间距离的似然函数，来得到本征维数的极大似然估计的。最大似然法是先假定一个小球面Sx(r)，r为半径，当r足够小时，认为密度采样函数f(x)近似为常量，假定观察样本服从齐次泊松分布，对非齐次过程N(t,x)，0≤t≤r进行考察，它表示的是以数据点x为中心，t距离内覆盖的样本的数目，则：

(1)

式中:I{xi∈Sx(r)}为示性函数，用泊松公式近似该二项式得到：

λ(t)=f(x)V(d)dtd-1

(2)

式中:V(d)为d维空间单位球的体积。令θ=lg(f(x))后，N(t)的似然函数形式如下：

(3)

(4)

(5)

式中:Tk(xi)为样本点xi与其第k个近邻点的欧氏距离。

1.2本征维数的自适应加权计算

对每个样本点，分别计算基于极大似然估计的本征维数，取其平均值即可得到全局本征维数：

(6)

取所有样本点的本征维数的平均值作为全局本征维数有其缺点，当数据集样本足够大且k选择合适的情况下，可以得到比较好的估计维数。然而，k值过小，估计值明显高于固有维数，当k值过大时，又有明显的负偏差。当样本数比较小时很难选择合适的k值。主要由于上述算法中直接取各个样本点估计维数的平均，而没有区分每个数据点的贡献，容易导致冗余信息放大、重要信息湮没，非常不利于后续降维操作。故本文对其算法进行改进，采用自适应加权平均算法：考虑每个数据点与其近邻点之间的关系，通过加权的方式修正数据点的贡献。

自适应权重计算：令Dk为近邻距离矩阵，Ak(x)为点x距最近邻的距离，则点x的权重可表示为：

(7)

修正后的本征维数为：

(8)

直接取平均值得方法问题在于位于稠密区的点给出过高的估计维数，相反稀疏区的点则给出过低的估计。而采用自适应权重后在一定程度上弥补了这一缺陷，原因是：

k较小时，大多数点位于稠密区，wk(x)→exp(1)=e，减少了正偏差；

k较大时，使得部分点变成稀疏点，wk(x)→exp(0)=1，对数据集维数贡献没有影响，可以削弱了负偏差的激增。

自适应加权估计根据数据点的局部分布，通过加权的方式修正了每个点对固有维数的贡献，在一定程度上提高了固有维数的估计准确率。

2基于LTSA的主流形识别

局部切空间排列算法的原理是在每个样本点处由其k个最近邻构造切空间，计算k个最近邻在切空间中一个正交坐标系下的坐标，然后对所有数据点处的局部坐标系进行平移缩放及旋转并排列，得到样本点的全局坐标。排列时最大可能地保留每个样本点处局部坐标系下的几何结构信息。基于LTSA的主流形识别基本过程如下：

(1)样本邻域选择：

对数据样本集X=[x1,x2,….xn]中的每个样本点，用欧式距离或最大范数确定k个距离最小的近邻点。设xi的k近邻矩阵为Xi=[xi1,xi2,….xik]。

(2)局部线性拟合

(3)全局坐标构建

设全局坐标为T={t1,t2,…,tn}∈Rd，Ti={ti1,ti2,…,tik}∈Rd是ti的不包括自身在内的k个近邻，每个样本点局部坐标矩阵排列起来得到全局坐标系，并且全局坐标反映局部坐标矩阵几何结构，其中：

(9)

(10)

式中:Ei=(εi1,εi2,…,εik),e为全1向量，长度为k，则

Ei=Ti×(I-e×eT/k)-LiΘi∈Rd

(11)

(4)主流形重构

得到每个数据点xi对应的低维全局坐标ti后，需要由低维坐标重构出高维空间中数据的主流形，计算如式(12)所示。

(12)

3基于自适应本征维数估计流形学习的相空间重构降噪流程

基于自适应本征维数估计流形学习的相空间重构降噪方法步骤如下：

(1)重构相空间

对一维时间序列为S=(s1,s2,…,sN)，运用互信息法和cao法[12]确定相空间重构最佳维数m和时间延迟τ，得到重构后的相空间矩阵为X=[x1,x2,…,xN-(m-1)τ]∈Rm×[N-(m-1)τ]。

(2)主流形识别

使用本文的自适应加权极大似然法估计出信号本征维数d，并使用LTSA将高维相空间降维到本征维数空间，设低维主流形为T=[t1,…,tN-(m-1)τ]∈Rd。

(3)主流形重构

(4)一维信号反求

根据相空间重构的方法反求一维时间序列。由于采用时间延迟的方式来重构相空间，所以一维时间序列中的某个元素可能会在相空间数据矩阵中的多个位置出现。从相空间的数据矩阵恢复到一维时间序列如式(13)所示。

(13)

式中:{Ii(j,k)}为时间序列中第i个元素在相空间数据矩阵中满足条件k+(j-1)τ=i的所有元素的下标集合，k∈[1,N-(m-1)τ]，j∈[1,m]，Ci为{Ii(j,k)}中元素的个数。

基于自适应本征维数估计流形学习和相空间重构的信号降噪流程见图1。

4仿真信号降噪实验

洛伦兹系统是典型的非线性动力系统：

(14)

式中：取d=10，b=8/3，r=28构造洛伦兹信号，采样频率为100 Hz，用四阶龙格库塔法计算，取x前20 s信号作为试验信号。洛伦兹信号的时域波形和二维相图见图2。

对洛伦兹信号添加高斯白噪声，信噪比为10 dB，带噪洛伦兹信号的的时域波形和二维相图见图3。

图2　洛伦兹信号时域波形和二维相图Fig.2 The time-domain waveform and two-dimensional phase diagram of Lorentz signal

图3　带噪洛伦兹信号时域波形和二维相图Fig.3 The time-domain waveform and two-dimensional phase diagram of Lorentz signal with noise

对添加白噪声后的洛伦兹信号进行相空间重构，使用互信息法计算相空间重构的时间延迟τ=1。利用已知的时延τ使用cao方法求嵌入维数，取嵌入位数m=13。自适应极大似然法估计高维相空间本征维数，其结果为d=3.257，取整为3。进一步进行降维，再重构主流形最后反求一维时间信号，其结果见图4。

图4　流形学习降噪后的洛伦兹信号时域波形和二维相图Fig.4 The time-domain waveform and two-dimensional phase diagram of Lorentz signal after de-noising with manifold learning

从洛伦兹信号时域波形来看，经过本文方法降噪后，基本能够还原原始信号，噪声强度已非常低，二维相图略有变形，说明还是存在少量噪声成分。为验证本文方法的有效性，这里同时选取常用的小波降噪方法与本文方法进行对比，选用“db8”小波，作5层分解，降噪效果见图5。

图5　小波降噪后的洛伦兹信号时域波形和二维相图Fig.5 The time-domain waveform and two-dimensional phase diagram of Lorentz signal after de-noising with wavelet

从时域波形来看，小波降噪与本文方法效果接近，都能很好的去处噪声，但从二维相位图来看，小波降噪后的数据有较大失真。由于小波降噪参数不好确定，难以根据信号非线性特性自适应降噪，效果不及本文方法。

为分析本征维数对降噪效果的影响，同样以上述添加白噪声洛伦兹信号为例，如果LTSA算法的目标维数不取3，而分别取d=4和d=5，其它参数不变，降噪后的洛伦兹信号的二维相位图(见图6)。可以看出，含噪洛伦兹信号噪声没能很好去除，如果目标维数不是本征维数，基于流形学习的非线性降噪效果将大打折扣，效果很不理想。

图6　目标维数分别为d=4(左)和d=5(右)时降噪后的二维相位图Fig.6 The two-dimensional phase diagram after de-noising with target dimension d=4(left) and d=5(right)

对不同信噪比的含噪洛伦兹信号分别用不同的目标维数(d=2、3、4、5、6、7)使用流形学习方法进行降噪处理，降噪后的信噪比见表1。由表1可知，在不同的信噪比下，都是d=3时降噪后的信噪比最高，即目标维数小于本征维数或者大于本征维数，会造成降噪过度或者降噪不够，只有目标维数设置为本征维数时，降噪效果最好。可见，目标维数是基于流形学习的非线性降噪方法的关键，必须正确估计出信号的本征维数，以本征维数作为目标维数，才能获得最佳降噪效果。

表1　嵌入维数和本征维数与信噪比的关系

5应用

对海装风电设备有限公司生产的5MW的风电机组传动系统进行振动测试，测试部位包括主轴承、齿轮箱转矩臂、齿轮箱支撑台架、齿轮箱传动轴和发电机。对于主轴部分，由于转速较低，振动多体现为低频，需使用低频加速度传感器。对于齿轮箱部分及发电机部分，振动多体现为高频，可使用标准加速度传感器进行采集。测试中所用采集设备为NI-9188采集机箱和NI9234采集模块。此次测试中，高速轴转速范围为800～1 200 r/min，齿轮箱传动比约为116，低速轴转速范围为6.9～10.3 r/min(0.12～0.18 Hz)，低频部位转频最低接近0.1 Hz，因此对于主轴承选用频率响应为0.1～2 000 Hz的朗斯振动加速度传感器LC0166T，对齿轮箱传动轴选用频率响应为1～10 000 Hz的SKF振动加速度传感器CMSS-WIND-100。低速部位传感器采用磁力座固定安装，高速部位传感器采用胶水粘结固定安装。风电机组传动系统振动信号测试现场见图7。

图7　风电机组传动系统振动信号测试现场Fig.7 The test bed of wind turbine vibration test

在高速轴转速为1 200 r/min时，测得各测点的振动信号，对主轴承轴向测点的振动信号进行降噪处理，其原始振动信号时域波形及其频谱见图8(a)。

主轴承上的振动频率基本在400 Hz以内，主要频率成分约为20 Hz及其倍频成分，20 Hz来源于齿轮箱输出端的转频，转频的倍频成分在轴向表现比较明显，而在径向表现相对较弱。除了转频及其倍频外，图示还有304 Hz的频率成分，经分析，该频率主要来自于行星轮的齿轮啮合频率，此外，由于测试机组附近的电磁干扰，从频谱图可以看出信号中还有一些噪声成分。

采用传统的小波消噪方法对该信号进行降噪处理。选择小波基为Daubechies小波“db8”，小波分解层数为5层，采用stein无偏风险估计原则估计阈值，软阈值去噪，降噪后的频谱结果见图8(b)。由图8(b)可知，小波消噪主要是带通滤波，通常是滤掉高频成分，但由于该测点中行星轮齿轮啮合频率为304 Hz，处于高频部分被滤掉，低频部分的噪声依然存在，显然没有起到有效的消噪效果。采用本文流形学习方法进行非线性降噪，通过互信息法计算得到的时间延迟τ=1，通过cao方法计算所得嵌入维数为86，基于自适应极大似然估计的得到本征维数为5，降噪后的时域信号和频谱见图8(c)。显然，经本文方法降噪后的频率成分比较清晰，20 Hz及其倍频成分都很明显，304Hz行星轮齿轮啮合频率成分有效的保留下来，噪声水平很低。

图8　主轴承轴向测点振动信号时域波形及相应频谱Fig.8 The time-domain waveform and frequency spectrum of vibration signal from axial direction of main bearing

为了进一步验证本文降噪方法的有效性，这里同样分别采用传统的小波消噪方法和本文流形学习方法对齿轮箱传动轴径向测点振动信号进行降噪处理，其原始振动信号时域波形及频谱见图9(a)。

由于风电机组齿轮箱结构复杂，主轴经过传动系统增速后，其频谱成分较低速轴更为丰富，在10 000 Hz以内均有分布，且噪声成分更为明显。通过小波消噪，处于图9(b)所示的部分有效高频频率成分被滤除，且低频部分的噪声依然存在。通过本文方法降噪后(见图9(c))，噪声成分得到较好的抑制，并保留了有效的高频成分。可见基于流形学习的信号消噪方法在低频和高频振动信号降噪中均能有效的消除了非线性噪声，有异于简单的通过滤波实现消噪，具有更好的自适应性。

图9　齿轮箱传动轴测点径向振动信号时域波形及相应频谱Fig.9 The time-domain waveform and frequency spectrum of vibration signal from radial direction of gearbox shaft

6结论

本文针对风电机组振动信号降噪提出了基于本征维数估计流形学习的非线性降噪法。通过风电机组振动信号降噪实验表明：

(1)结合流形学习和相空间重构进行非线性降噪，更能体现动力学系统主流形的全域信息和整体结构，具有更好的自适应性；

(2)采用极大似然法的本征维数估计法，可以避免目标维数设置不当带来的不完全降噪或过度降噪，并在此基础上进行自适应加权处理，修正数据点的贡献，从而进一步提高风电机组传动系统诊断信号的非线性降噪效果；

(3)与常用的小波降噪法进行对比，结果表明基于本征维数估计流形学习的非线性降噪法具有更好的降噪性能，进一步验证了本方法的有效性。

参考文献

[1]Jiang Y H, Tang B P, Qin Y, et al. Feature extraction method of wind turbine based on adaptive Morlet wavelet and SVD[J]. Renewable Energy, 2011, 36(8): 2146-2153.

[2]Zhong M C, Zhang C L, Cheng L, et al. Vibrationsignal of recoil mechanism de-noising processing based on wavelet analysis[J]. Applied Mechanics and Materials, 2013(401): 1424-1427.

[3]Takens F, Detecting sange attractors in turbulence [R]． Lecture Notes in Math．New York： SpringeL 1981．

[4]Yin J, Hu D, Zhou Z. Noisy manifold learning using neighborhood smoothing embedding[J]. Pattern Recognition Letters, 2008, 29(11): 1613-1620.

[5]Zhang Z, Zha H. Local linear smoothing for nonlinear manifold learning[R]. Technical Report, CSE, Penn State Univ, 2003.

[6]Chang H, Yeung D Y. Robust locally linear embedding[J]. Pattern Eecognition, 2006, 39(6): 1053-1065.

[7]阳建宏,徐金梧,杨德斌,等. 基于主流形识别的非线性时间序列降噪方法及其在故障诊断中的应用[J].机械工程学报,2006,42(8):154-158.

YANG Jian-hong, XU Jin-wu, YANG De-bin, et al. Noise reduction method for nonlinear time series based on principal manifold learning and its application to fault diagnosis[J].Chinese Journal of Mechanical Engineering, 2006, 42(8): 154-158.

[8]王广斌, 黄良沛, 赵先琼.基于分形维和局部切空间均值重构的非线性降噪方法[J]. 电子测量与仪器学报, 2010(8): 699-704.

WANG Guang-bin, HUANG Liang-pei, ZHAO Xian-qiong.Nonlinear noise reduction method based on fractal dimension and the local tangent space mean reconstruction[J]. Journal of Electronic Mesaurement and Instrument, 2010(8): 699-704.

[9]余成义.流形学习在混沌时间序列降噪中的应用[D]. 武汉： 武汉科技大学, 2012.

[10]Xiao H B, Xie X P, Zhou S Q, et al. Featureextraction of diesel engine wear fault based on local tangent space alignment[J]. Applied Mechanics and Materials, 2014(482): 179-182.

[11]Fu B, Wang X, Yang B, et al. Anovel approach for intrinsic dimension estimation based on maximum likelihood estimation[C]//Proceedings of the 2012 International Conference on Cybernetics and Informatics. Springer New York,2013:91-97.

[12]Cao Lianyue. Practical method for determining the minimum embedding dimension of a scalar time series[J]. Physica D : Nonlinear Phenomena, 1997(110)：43-50.