王海軍 許飛云

(東南大學機械工程學院,南京211189)

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

1 Tucker 3分解模型

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

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

(1)

圖1 Tucker 3分解模型

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

(2)

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

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

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

M←M-H-1r

(3)

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

(4)

3 SCA的二次特征處理

3.1 信號的稀疏處理

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

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

(5)

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

(6)

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

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

(7)

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

3.2 混合矩陣估計

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

Y=HS

(8)

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

min‖Y-HS‖

(9)

當信號損失足夠小時,對方程(9)求最優(yōu)解,得到

H=YS?

(10)

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

圖2 混合信號的SCA處理流程

4 齒輪箱信號分解驗證

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

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

圖3 齒輪箱結構圖

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

(11)

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

圖4 迭代收斂誤差比較

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

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

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

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

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

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

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

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

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

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

5 結語

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

)

[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)仿真學報,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)