基于自適應本征維數估計流形學習的相空間重構降噪方法

馬婧華，湯寶平，宋　濤

(重慶大學機械傳動國家重點實驗室，重慶　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個最近鄰在切空間中一個正交坐標系下的坐標，然后對所有數據點處的局部坐標系進行平移縮放及旋轉并排列，得到樣本點的全局坐標。排列時最大可能地保留每個樣本點處局部坐標系下的幾何結構信息?；贚TSA的主流形識別基本過程如下：

(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))，噪聲成分得到較好的抑制，并保留了有效的高頻成分?？梢娀诹餍螌W習的信號消噪方法在低頻和高頻振動信號降噪中均能有效的消除了非線性噪聲，有異于簡單的通過濾波實現消噪，具有更好的自適應性。

圖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.