基于非負Tucker 3分解的稀疏分量分析在故障信號提取中的應用

王海軍 許飛云

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

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

1 Tucker 3分解模型

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

圖1 Tucker 3分解模型

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

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

對式(2)中的A(n)分別逐個求偏導，便可得到模矩陣 A(1)，A(2)，A(3)的計算式.對 A(n)進行交替迭代計算，并增加對A(n)非負約束，可得到張量核G.非負約束對計算陷于局部過擬合起到了較大的抑制作用［13］.但是，這種迭代計算方式會產生巨大的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］)，將各堆疊的矩陣重新排列成一行.根據牛頓-高斯梯度下降迭代法，得到矩陣M的更新方程為

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

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

3 SCA的二次特征處理

3.1 信號的稀疏處理

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

式中，σ為 Y與 ^Y之間的協方差.在這里，cexp(-μ(Y-G?{A})2與離散信號的窗函數作用類似.

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

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

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

3.2 混合矩陣估計

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

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

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

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

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

4 齒輪箱信號分解驗證

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

圖3 齒輪箱結構圖

取雙譜2個正頻率軸的頻率點數為64，將分別采集到的信號加噪后組成255組包含3種故障狀態的數據，即構成一個Ω×Ω×S的張量，其中Ω=64，S=255.取NTD后的張量核為32×32×64，設定NTF分解因子維數為323，迭代中的收斂誤差為

將本文的迭代收斂誤差與傳統的NTF進行比較，如圖4所示.

圖4 迭代收斂誤差比較

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

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

圖5 初始3種故障狀態的頻譜圖

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

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

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

表1 齒輪箱故障特征的稀疏值個數

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

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

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

5 結語

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

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］.系統仿真學報，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)