基于非負(fù)Tucker 3分解的稀疏分量分析在故障信號(hào)提取中的應(yīng)用

王海軍 許飛云

(東南大學(xué)機(jī)械工程學(xué)院,南京211189)

Tucker 3分解模型是Tucker[1]針對(duì)多維數(shù)據(jù)降解而提出的一種高效數(shù)學(xué)分解模型.隨著計(jì)算機(jī)技術(shù)的發(fā)展,在交替最小二乘算法的基礎(chǔ)上,Paatero等[2]建立了三維CANDECOMP/PARAFAC(CP)并行計(jì)算模型,CP模型對(duì)非負(fù)分解算法的發(fā)展起到了很大的推動(dòng)作用.從實(shí)體數(shù)據(jù)的有效性出發(fā),Lee等[3]提出并證實(shí)了非負(fù)矩陣分解方法對(duì)圖像的局部特征具有良好的解釋性.自此,非負(fù)矩陣分解方法在盲信號(hào)處理、圖像特征提取、神經(jīng)系統(tǒng)學(xué)和化學(xué)計(jì)量學(xué)等領(lǐng)域中得到了廣泛的應(yīng)用.目前,對(duì)三維目標(biāo)特征的局部化處理和分類仍然是研究中的熱點(diǎn)和難點(diǎn).這樣便引出了對(duì)算法的稀疏性控制優(yōu)化等方面的研究.在特征稀疏性處理上比較典型的研究主要包括利用拉格朗日優(yōu)化、增加并行因子補(bǔ)償項(xiàng)和稀疏性控制項(xiàng)等方法來(lái)表達(dá)固有分量的特征[4-9].在較小規(guī)模數(shù)據(jù)的計(jì)算上,這些方法能達(dá)到預(yù)期的目的.但是,對(duì)于高維大規(guī)模數(shù)據(jù)的計(jì)算往往會(huì)帶來(lái)收斂速度慢和局部過(guò)擬合問(wèn)題,而這些問(wèn)題會(huì)直接影響到計(jì)算的精度和二次特征的有效表達(dá).因此,針對(duì)過(guò)擬合和特征稀疏性的問(wèn)題,本文提出了非負(fù)Tucker 3分解(NTD)結(jié)合稀疏分量分析(SCA)的方法來(lái)提高分解的效率以及二次信號(hào)特征的稀疏性;同時(shí),對(duì)Tucker 3的分解因子進(jìn)行非負(fù)約束和一次性更新計(jì)算以提高計(jì)算的精度和效率.這樣提取出的二次特征在稀疏性和精度方面也會(huì)得到有效的控制.因而,該方法的研究從理論上來(lái)說(shuō)顯得很有必要,在實(shí)踐中也具有重要意義.

1 Tucker 3分解模型

作為一種成功的張量分解方法,Tucker提出的分解方法可以簡(jiǎn)單地概括為3種模型,即Tucker 1,Tucker 2和Tucker 3分解模型.其中,Tucker 3分解模型又是前2種方法研究的延續(xù)和發(fā)展,在非負(fù)張量分解(NTF)的特征局部化處理與應(yīng)用中具有十分重要的意義[10-12].根據(jù)NTF的優(yōu)越性,本文主要以Tucker 3分解模型作為研究對(duì)象.

s.t. ‖A(n)‖=1

(1)

圖1 Tucker 3分解模型

對(duì)式(1)的最佳近似進(jìn)行分解,可轉(zhuǎn)化為求解以下最優(yōu)化方程:

(2)

對(duì)式(2)中的A(n)分別逐個(gè)求偏導(dǎo),便可得到模矩陣A(1),A(2),A(3)的計(jì)算式.對(duì)A(n)進(jìn)行交替迭代計(jì)算，并增加對(duì)A(n)非負(fù)約束，可得到張量核G.非負(fù)約束對(duì)計(jì)算陷于局部過(guò)擬合起到了較大的抑制作用[13].但是,這種迭代計(jì)算方式會(huì)產(chǎn)生巨大的Jacobian矩陣，同時(shí)也帶來(lái)了收斂慢和效率低的問(wèn)題[12].因此,研究一種更高效合理的計(jì)算方法是本文的研究目的之一.

2 牛頓-高斯梯度下降的迭代方法

對(duì)分解因子A(1),A(2),A(3)以及G進(jìn)行重新組合,得到一個(gè)新的矩陣M=[A(1)T,A(2)T,…,A(N)T,vec(G)].其中,vec(·)表示展開(kāi)堆疊的張量G(關(guān)于張量的展開(kāi),可參見(jiàn)文獻(xiàn)[10]),將各堆疊的矩陣重新排列成一行.根據(jù)牛頓-高斯梯度下降迭代法,得到矩陣M的更新方程為

M←M-H-1r

(3)

式中,H為海森矩陣;r為梯度矩陣，計(jì)算公式為

(4)

3 SCA的二次特征處理

3.1 信號(hào)的稀疏處理

高斯-笛卡爾密度核函數(shù)是由Khoromskij等[14]針對(duì)三維PARAFAC分解因子計(jì)算而提出，并應(yīng)用于冗余化學(xué)原子庫(kù)信號(hào)稀疏化的方法.高斯-笛卡爾密度核函數(shù)不僅能處理離散化信號(hào),而且還具有濾波降噪的作用.因此,本文主要根據(jù)NTD各因子間叉積的特點(diǎn),結(jié)合高斯-笛卡爾積,建立聯(lián)合核函數(shù)為

Φ(Y)=c(Y-G×{A})exp(-μ(Y-G?{A})2)

(5)

假如Φ(Y)為冗余完備庫(kù),Y為其觀測(cè)信號(hào)，則令Φ(Y)為一高斯原子,與Y信號(hào)長(zhǎng)度相同,均進(jìn)行歸一化處理.兩者間的最優(yōu)化核函數(shù)求解可轉(zhuǎn)化為解決以下映射內(nèi)積的優(yōu)化問(wèn)題:

(6)

式中,〈Y,Φ(Y0)〉表示Y與Φ(Y0)的內(nèi)積;Φ(Y0)為Y與Φ(Y)間的最佳原子庫(kù);α為充分考慮信號(hào)長(zhǎng)度和Φ(Y0)損失等原因的常數(shù).則信號(hào)最終分解為2部分:一部分為最佳高斯核函數(shù);另一部分為分配后的殘余信號(hào).其數(shù)學(xué)表達(dá)式為

〈Y,Φ(Y)〉=〈Y,Φ(Y0)〉Φ(Y0)+rY

(7)

式中,〈Y,Φ(Y0)〉Φ(Y0)表示觀測(cè)信號(hào)在Φ(Y0)上的最佳映射;rY為映射后的殘余信號(hào).實(shí)際上,高斯核函數(shù)與NTD后各因子參數(shù)直接相關(guān).而根據(jù)交替迭代計(jì)算方式,核張量為G=Y×{A}T.因此,模矩陣的計(jì)算直接影響著高斯核函數(shù)的質(zhì)量.這也充分說(shuō)明了更新算法在迭代計(jì)算中很重要.

3.2 混合矩陣估計(jì)

混合矩陣是稀疏分量處理的重要組成部分,其精確度直接決定著信號(hào)分離的結(jié)果.假設(shè)稀疏化處理后的混合信號(hào)由Y1,Y2,…,YN子信號(hào)組成,將各分量進(jìn)行分層處理,令yn=Yn(:,:,i)∈Yn,1≤n≤N.其中，y表示m×N的觀測(cè)矩陣.SCA類似獨(dú)立分量分析[15],主要解決以下線性信號(hào)分解問(wèn)題：

Y=HS

(8)

式中,S為n×N稀疏源信號(hào)矩陣，N為信號(hào)樣本；H=(hi,j)(i=1,2,…,m;j=1,2,…,n)為未知的混合矩陣.為了保證在信號(hào)特征依然稀疏的情況下能得到最佳混合矩陣,演化成解的最優(yōu)化問(wèn)題，即

min‖Y-HS‖

(9)

當(dāng)信號(hào)損失足夠小時(shí),對(duì)方程(9)求最優(yōu)解,得到

H=YS?

(10)

式中,?表示偽逆.估算出H矩陣后,用最小交替迭代二乘法對(duì)稀疏源信號(hào)矩陣進(jìn)行逼近計(jì)算.整個(gè)信號(hào)特征提取的流程如圖2所示.

圖2 混合信號(hào)的SCA處理流程

4 齒輪箱信號(hào)分解驗(yàn)證

為了驗(yàn)證該算法的稀疏性和可靠性,采用東南大學(xué)故障診斷研究所的3種齒輪箱的故障數(shù)據(jù).實(shí)驗(yàn)設(shè)備包括3套齒輪箱-電機(jī)系統(tǒng),將位于中間的單級(jí)齒輪箱作為測(cè)試對(duì)象.齒輪箱的內(nèi)部結(jié)構(gòu)原理如圖3所示.3個(gè)齒輪箱中的主動(dòng)輪分別設(shè)置為正常、齒面點(diǎn)蝕和均勻磨損3種故障狀態(tài).齒輪箱的輸入軸通過(guò)剛性聯(lián)軸器與電機(jī)相連,轉(zhuǎn)速可由Siemens MicroMaster420控制器進(jìn)行調(diào)節(jié).主動(dòng)輪與從動(dòng)輪間的傳動(dòng)比為31∶46.在電機(jī)轉(zhuǎn)速約為4 000 r/min的情況下,通過(guò)分別安裝在齒輪箱上垂直和水平方向上的壓電傳感器采集振動(dòng)信號(hào).如圖3所示,傳感器靈敏度為100 mV/g,誤差范圍為±3 dB,采樣頻率為3 838 Hz.分別對(duì)3種故障狀態(tài)每種采集20組振動(dòng)信號(hào),每組長(zhǎng)度為4 096點(diǎn).齒輪的嚙合頻率和滾動(dòng)軸承外圈通過(guò)頻率分別為310和99.7 Hz.采樣頻率約為10 kHz.

取雙譜2個(gè)正頻率軸的頻率點(diǎn)數(shù)為64,將分別采集到的信號(hào)加噪后組成255組包含3種故障

圖3 齒輪箱結(jié)構(gòu)圖

狀態(tài)的數(shù)據(jù),即構(gòu)成一個(gè)Ω×Ω×S的張量,其中Ω=64,S=255.取NTD后的張量核為32×32×64,設(shè)定NTF分解因子維數(shù)為323，迭代中的收斂誤差為

(11)

將本文的迭代收斂誤差與傳統(tǒng)的NTF進(jìn)行比較,如圖4所示.

圖4 迭代收斂誤差比較

設(shè)定迭代停止誤差為10-3，在訓(xùn)練張量維數(shù)相同的情況下,NTD達(dá)到目標(biāo)精度所需的迭代步數(shù)約為150,收斂效率明顯高于NTF.另一方面,在迭代誤差計(jì)算過(guò)程中,NTD的收斂誤差較平滑,從而說(shuō)明了該算法具有良好的健壯性.因此,從收斂效率和健壯性兩方面看,本文算法均優(yōu)于NTF.

為了讓特征信息更加容易識(shí)別,本實(shí)驗(yàn)需對(duì)故障信號(hào)進(jìn)行時(shí)頻變換.隨機(jī)選取3種故障數(shù)據(jù),進(jìn)行快速FFT變換后得到的初始狀態(tài)頻譜圖如圖5所示.

圖5 初始3種故障狀態(tài)的頻譜圖

由圖5可見(jiàn),3種狀態(tài)對(duì)應(yīng)的振動(dòng)頻率分布在99.7 Hz倍頻時(shí)的概率較大.對(duì)于正常狀態(tài)和均勻磨損狀態(tài),初始信號(hào)的二次特征并不容易判斷識(shí)別.在電機(jī)高速旋轉(zhuǎn)情況下,點(diǎn)蝕狀態(tài)振動(dòng)明顯,頻譜特征相對(duì)容易判斷.但是,如果在噪聲干擾下,頻譜特征分布不均勻,稀疏性差,信號(hào)的優(yōu)勢(shì)頻率并不突出.類似地,正常狀態(tài)和均勻磨損狀態(tài)的信號(hào)特征稀疏性更需要改進(jìn).對(duì)此,本實(shí)驗(yàn)將采用上面提出的SCA與NTD相結(jié)合(SCA_NTD)的方法提取信號(hào)的二次故障特征，其頻域特征如圖6所示.

圖6 SCA_NTD提取出的齒輪故障頻譜圖

與初始信號(hào)的二次特征相比,SCA_NTD提取出的特征頻率能滿足周期性的特點(diǎn),也與齒輪減速箱嚙合頻率和軸承通過(guò)外圈頻率相符合:點(diǎn)蝕狀態(tài)的振幅值對(duì)應(yīng)于基頻99.7和310 Hz的多倍頻;均勻狀態(tài)對(duì)應(yīng)基頻為99.7 Hz的倍頻，與齒輪嚙合頻率也相符.另外,特征信號(hào)的稀疏性比較好,容易觀測(cè),易于判斷.在此,將特征值小于10-6近似作為特征信號(hào)的稀疏值.初始狀態(tài)與SCA_NTD方法處理后的特征稀疏值個(gè)數(shù)的結(jié)果如表1所示.不難發(fā)現(xiàn),處理后的特征稀疏值個(gè)數(shù)較多,從而說(shuō)明了SCA_NTD處理后的盲信號(hào)具有良好的稀疏性.

表1 齒輪箱故障特征的稀疏值個(gè)數(shù)

為了證明SCA_NTD的可靠性,將其與經(jīng)典的交替最小二乘法的NTD(ALS_NTD)和非負(fù)張量分解NTF算法進(jìn)行比較.計(jì)算精度為Ac=(1-Et)×100%.實(shí)驗(yàn)結(jié)果如表2所示.

由表2中的精確度分布可見(jiàn),相同計(jì)算方法的精度隨著張量維數(shù)的增大而增高.總體上看,SCA_NTD計(jì)算得到的精度要比ALS_NTD和NTF高.隨著張量核維數(shù)的調(diào)整,SCA_NTD的最高精度達(dá)到了97.16%,相比ALS_NTD與NTF的最高精度93.93%和88.81%,優(yōu)勢(shì)明顯.從張量核維數(shù)組合情況可看出,當(dāng)核張量維數(shù)約為張量維數(shù)的一半時(shí),精度最高.因此,根據(jù)這一規(guī)律合理選擇張量核維數(shù),SCA_NTD的可靠性將會(huì)進(jìn)一步得到保證.

表2 SCA_NTD與ALS_NTD在不同張量維數(shù)下的精確度 %

5 結(jié)語(yǔ)

針對(duì)NTD算法提取的二次特征信號(hào)不稀疏問(wèn)題,結(jié)合SCA二次分離的方法得到了更加稀疏的特征信號(hào).在處理SCA的混合矩陣問(wèn)題時(shí),采用了交替迭代計(jì)算稀疏源偽逆矩陣的方法.同時(shí),為了避免NTD在迭代過(guò)程中陷于局部過(guò)擬合而導(dǎo)致誤差增大和效率降低的問(wèn)題,提出了一次更新所有分解因子的方法.實(shí)驗(yàn)結(jié)果表明，SCA_NTD達(dá)到了改善二次特征信號(hào)的稀疏性以及提高了計(jì)算精度和效率的目的.

)

[1] Tucker L R. Some mathematical notes on three-mode factor analysis[J].Psychometrika,1966,31(3): 279-311.

[2] Paatero P,Tapper U. Positive matrix factorization: a non-negative factor model with optimal utilization of error estimates of data values [J].Environmetrics,1994,5(2): 111-126.

[3] Lee D D,Seung H S. Learning the parts of objects by non-negative matrix factorization[J].Nature,1999,401(6755): 788-791.

[4] Hazan T,Polak S,Shashua A. Sparse image coding using a 3D non-negative tensor factorization[C]//10thIEEEInternationalConferenceonComputerVision. Beijing,2005: 50-57.

[5] Morup M,Hansen L K,Arnfred S M. Algorithms for sparse nonnegative Tucker decompositions [J].NeuralComputation,2008,20(8): 2112-2131.

[6] Cichocki A,Zdunek R,Phan A H,et al.AlternatingleastsquaresandrelatedalgorithmsforNMFandSCAproblemsinnonnegativematrixandtensorfactorizations[M]. Chichester,UK: John Wiley & Sons,Ltd,2009: 203-266.

[7] Peng S,Xu F,Jia M,et al. Sparseness-controlled non-negative tensor factorization and its application in machinery fault diagnosis[J].JournalofSoutheastUniversity:EnglishEdition,2009,25(3): 346-350.

[8] Karoui M S,Deville Y,Hosseini S,et al. Blind spatial unmixing of multispectral images: new methods combining sparse component analysis,clustering and non-negativity constraints[J].PatternRecognition,2012,45(12): 4263-4278.

[9] Asaei A,Davies M E,Bourlard H,et al. Computational methods for structured sparse component analysis of convolutive speech mixtures[C]//IEEEInternationalConferenceonAcoustics,SpeechandSignalProcessing. Kyoto,Japan,2012: 2425-2428.

[10] Kolda T G.Multilinearoperatorsforhigher-orderdecompositions[M]. California,USA: Sandia National Laboratories,2006.

[11] Cai X J,Chen Y N,Han D R. Nonnegative tensor factorizations using an alternating direction method[J].FrontiersofMathematicsinChina,2013,8(1): 3-18.

[12] Jiang L L,Yin H Q. Bregman iteration algorithm for sparse nonnegative matrix factorizations via alternating l (1)-norm minimization [J].MultidimensionalSystemsandSignalProcessing,2012,23(3): 315-328.

[13] Albright R,Cox J,Duling D,et al. Algorithms,initializations,and convergence for the nonnegative matrix factorization[R]. Raleigh,USA: Carolina State University,2006.

[14] Khoromskij B,Khoromskaia V,Chinnamsetty S,et al. Tensor decomposition in electronic structure calculations on 3D Cartesian grids[J].JournalofComputationalPhysics,2009,228(16): 5749-5762.

[15] 劉海林,姚楚君.欠定混疊稀疏分量分析的超平面聚類算法 [J].系統(tǒng)仿真學(xué)報(bào),2009(7):1826-1828.

Liu Hailin,Yao Chujun. Hyperplane clustering algorithm of underdetermined mixing sparse component analysis[J].JournalofSystemSimulation,2009(7):1826-1828. (in Chinese)