基于非負(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ì)象.

圖1 Tucker 3分解模型

式中，×1表示矩陣形式的模數(shù)積；{A}表示所有模矩陣的Kronecker積；Y為Y被分解后的近似值.

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

對(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的更新方程為

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

為了避免海森矩陣H出現(xiàn)極值為零而影響計(jì)算的效率，取H的近似矩陣，以保證極值不為零，令=H+uI，其中0＜u?1.

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與 ^Y之間的協(xié)方差.在這里，cexp(-μ(Y-G?{A})2與離散信號(hào)的窗函數(shù)作用類似.

假如Φ(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)題:

式中，〈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，Φ(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)題:

式中，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)題，即

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

式中，?表示偽逆.估算出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)速約為4000 r/min的情況下，通過(guò)分別安裝在齒輪箱上垂直和水平方向上的壓電傳感器采集振動(dòng)信號(hào).如圖3所示，傳感器靈敏度為 100 mV/g，誤差范圍為±3 dB，采樣頻率為3838 Hz.分別對(duì)3種故障狀態(tài)每種采集20組振動(dòng)信號(hào)，每組長(zhǎng)度為4096點(diǎn).齒輪的嚙合頻率和滾動(dòng)軸承外圈通過(guò)頻率分別為310和99.7 Hz.采樣頻率約為10 kHz.

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

取雙譜2個(gè)正頻率軸的頻率點(diǎn)數(shù)為64，將分別采集到的信號(hào)加噪后組成255組包含3種故障狀態(tài)的數(shù)據(jù)，即構(gòu)成一個(gè)Ω×Ω×S的張量，其中Ω=64，S=255.取NTD后的張量核為32×32×64，設(shè)定NTF分解因子維數(shù)為323，迭代中的收斂誤差為

將本文的迭代收斂誤差與傳統(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ì)算精度和效率的目的.

References)

［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］//10th IEEE International Conference on Computer Vision.Beijing，2005:50-57.

［5］Morup M，Hansen L K，Arnfred S M.Algorithms for sparse nonnegative Tucker decompositions［J］.Neural Computation，2008，20(8):2112-2131.

［6］Cichocki A，Zdunek R，Phan A H，et al.Alternating least squares and related algorithms for NMF and SCA problems in nonnegative matrix and tensor factorizations［M］.Chichester，UK:John Wiley ＆ Sons，Ltd，2009:203-266.

［7］Peng S，Xu F，Jia M，et al.Sparseness-controlled nonnegative tensor factorization and its application in machinery fault diagnosis［J］.Journal of Southeast University:English Edition，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］.Pattern Recognition，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］//IEEE International Conference on Acoustics，Speech and Signal Processing.Kyoto，Japan，2012:2425-2428.

［10］Kolda T G.Multilinear operators for higher-order decompositions［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］.Frontiers of Mathematics in China，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］.Multidimensional Systems and Signal Processing，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］.Journal of Computational Physics，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］.Journal of System Simulation，2009(7):1826-1828.(in Chinese)