李紀永, 李舜酩, 田國成, 陳曉紅
(1.南京航空航天大學能源與動力學院 南京,210016) (2.山東中實易通集團有限公司 濟南,250000)(3.南京航空航天大學理學院 南京,210016)
?
轉子共頻相關故障源源數估計與子帶盲分離*
李紀永1, 李舜酩1, 田國成2, 陳曉紅3
(1.南京航空航天大學能源與動力學院 南京,210016) (2.山東中實易通集團有限公司 濟南,250000)(3.南京航空航天大學理學院 南京,210016)
針對轉子異常振動產生含交叉頻率的響應,其共頻相關故障源不滿足統計獨立要求,提出利用非負矩陣法在頻域中計算故障源個數,不計及源信號和混合系統特性,可以正確估計出故障源數目或源數上限。提出利用小波包分解故障信號,選擇互信息較小的子帶進行重構,剔除共頻信號并進行盲分離,得到獨立非相關的源信號,保留了故障信息。理論及實驗結果證明了所提出方法的有效性。
相關源; 非負矩陣分解; 互信息; 子帶分解; 盲源分離
航空發動機在研制和生產試車中經常出現異常振動現象,非相關源及存在交叉頻率的相關源振動信號混合在一起,不滿足源信號統計獨立性條件[1-3]。傳統的源數估計方法無法正確識別源信號數目,而基于標準獨立量分析方法無法正確分離出此類信號[4-5]。文獻[6]利用功率譜密度方法估計相關源數,在理論上估計出上下界,但不能確定具體數目。文獻[4]指出利用高通濾波器過濾相關信號,然后利用獨立量分析分離。基于高通濾波后的信號往往比原信號更獨立這一特點,在實際應用中并不理想,指出自適應濾波方法需要源信號的先驗概率分布。
筆者利用非負矩陣分解 (non-negative matrix factorization,簡稱NMF)方法[7]分解觀測信號的功率譜密度系數,進而估計相關源數。非負矩陣對源信號無統計性要求,其純加性的和對稀疏性的描述特性使得描述數據更加清晰。利用子帶盲分離算法[1](sub-band decomposition independent component analysis,簡稱SDICA)對觀測信號分離,得到觀測信號子帶,計算子帶互信息[8](mutual information,簡稱MI),剔除互信息最大的子帶進行信號重構。該算法魯棒性強,是分離相關源振動的有效方法。
NMF數學表達形式為
Vm×n≈Wm×rHr×n
(1)
其中:V,W及H均為非負矩陣。
以V及W×H的歐幾里德距離平方(square of eculidian distance,簡稱SED)為目標函數,約束條件為H的每一列歸一化,表達式為
(2)
文獻[7]提出了乘性迭代規則,證明其單調下降收斂特性。迭代規則為
(3)
式(1)形式上與盲信號分離線性瞬時混合模型相似,其表達式為
x(t)m×n=Am×ps(t)p×n
(4)
其中:x(t)為觀測信號;s(t)為源信號;A為混合矩陣。
實際上所觀測振動信號不滿足非負條件,故將其進行傅里葉變換,得到
x(w)m×q=Am×rs(w)r×q=Am×rBr×ps(w)p×q
(5)
令Am×p=Am×rBr×p,Br×p為幅值系數矩陣。由于V,A及B均為正矩陣,滿足非負矩陣分解條件。利用乘性迭代規則單調下降收斂特性,求解源信號與降維信號歐式距離首先降為0時的r即為源數。
含有同頻的各觀測信號是相關的[9],基于獨立分量分析的盲分離技術不能直接將相關信號分離,筆者利用小波包分解,將相關的子頻帶信號即共頻信號剔除,重構后進行基于獨立量分析的盲分離,即可得到分離信號。
源信號x(t)的第i個源信號si(t)可以表示為
si(t)=si,1(t)+si,2(t)+si,3(t)+…+si,L(t)
(6)
其中:si,1(t),si,2(t),si,3(t),…,si,L(t)為源信號si(t)的各個子帶分量。
信號通過多分辨率小波包分解,相當于通過一個迭代的帶通濾波器Ti[10-11]。源信號矩陣信號s經過這個濾波器得到子帶分量為
si(t)=Ti[s(t)]
(7)
子帶盲分離的數學模型為
xi(t)=Ti[As(t)]=ATi[s(t)]=Asi(t)
(8)
由式(8)可知,混合矩陣通過濾波器未發生變化,觀測信號子帶分量是混合矩陣與源信號子帶分量的乘積,可利用觀測信號子帶盲分離進行相關源分析。利用互信息值大小表征子帶信號之間的相關性。當互信息為1時,子帶信號完全相關;互信息為0時,子帶信號不相關。為編程方便,利用高階累積量表達互信息值,表達式為
(9)
選取互信息較小的子帶重構信號,利用快速獨立量分析[12]進行盲分離,即可得到剔除同頻分量的故障源信號。
構造一組相關振動信號
(10)
其中:s0為s1,s2,s3,s4的共頻信號。
圖1為未添加共頻信號時信號波形。首先對觀測信號進行源數識別,當傳感器數目大于源數時(這里設源數為4,傳感器數為5),其混合矩陣A1為對系數矩陣進行非負矩陣分解,得到源數與目標函數關系如圖2,3所示。當傳感器數目大于源數時,可正確估計信號源數;當傳感器數目小于源數時,可估計信號源數的下限。

圖1 仿真信號時間歷程Fig.1 Simulation signal time sierious

圖2 超定情況下最小目標函數值Fig.2 Min objective function value in over-determination case

圖3 欠定情況下源數與最小目標函數值Fig.3 Min objective function value in under-determination case

隨機生成一個混合矩陣,對源信號進行混合,利用FastICA方法分離,恢復出的源信號如圖4所示,其Amari差值計算為0.218 2[8]。由y1的頻譜可以看出共頻信號依然存在,并未分出,分離失敗。

圖4 FastICA恢復信號時間歷程與頻譜圖Fig.4 FastICA recovery signal time series and frequency spectrum
利用小波包子帶分解,分解層數為1時相當于采用高通濾波器提取其高頻信號進行盲分離,其Amari差值計算為0.11,亦未分離出共頻信號,采用5層小波包子帶分解,計算各層子帶間互信息值,提出互信息最大的子帶對信號進行重構,然后進行盲分離,恢復信號及其頻譜圖如圖5所示。可以看出,共頻信號被分離出來,計算出的Amari差值為0.001 8,較好地分離出了共頻信號,并保留了原始信號特征信息。
在單轉子實驗臺測取4 kr/min時基礎松動-碰摩-偏心耦合故障,實驗裝置如圖6所示。對測得的信號進行去均值并進行相關降噪處理,利用非負矩陣法計算源數,計算得源數為3,如圖7所示。

圖5 子帶分解恢復信號時間歷程與頻譜圖Fig.5 Sub-band decomposition blind source separation time series and frequency spectrum

圖6 轉子實驗臺Fig.6 Rotor test rig

圖7 實驗信號源數與最小目標函數Fig.7 Experimental signal source number and min objective function value
利用自相關方法對觀測的故障信號進行降噪,然后直接進行盲源分離,所得頻譜如圖8所示。各路信號同頻成分多,主要包括50,250,500 Hz等,說明各傳感器測得的源信號相關性強,利用相關方法信號降噪并對降噪后的信號進行小波包子帶分解,分解層數為6,計算其互信息值,選取適當閾值去除同頻分子帶分量,然后對不相關的子帶進行重構,經過盲源分離所得的頻譜結果如圖9所示。可明顯看出,信號噪聲降低,頻率成分簡潔,s1信號工頻顯著,為不平衡故障,s2信號倍頻成分豐富,噪聲成分小,可判斷為碰摩故障,s3信號倍頻成分豐富, 但噪聲較多,可認為由基礎松動振動不斷沖擊底座造成。此方向信號較強,四倍頻被剔除掉,說明各源信號存在4倍頻,以上推測符合碰摩故障特征。

圖8 降噪故障信號盲分離源信號頻譜Fig.8 Denoising fault signal blind source separation spectrum

圖9 實驗信號子帶盲分離頻譜Fig.9 Experimental signal sub-band decomposition blind source separation spectrum
1) 利用非負矩陣乘性迭代規則收斂特性進行故障信號降噪并在頻域進行分解,無需考慮源信號及混合系統特性,尋找代價函數最小值即可求得源的個數。當傳感器數大于源數時,可正確估計出源數,當傳感器數目小于源數時,可估計其下限。
2) 利用相關小波包子帶分解故障信號,計算其互信息,選取閾值剔除同頻相關信號,進行信號重構后,保留了故障信號的信息。
3) 對重構后的信號進行盲分離,提取出不相關的信號,與直接進行盲分離的結果進行對比,更易于診斷故障及分析因果關系,從而得出系統振動源故障特性。
[1] 周曉峰,楊世錫,甘春標.相關機械振源的盲源分離方法[J].振動與沖擊,2012,31(14):60-64.
Zhou Xiaofeng, Yang Shixi, Gan Chunbiao. Blind source separation of statistically correlated sources[J]. Journal of Vibration and Shock, 2012,31(14):60-64. (in Chinese)
[2] 杜建建,李舜酩,張袁元,等.平板葉片的相關激勵響應實驗[J].實驗室研究與探索,2011,30(1):14-17.
Du Jianjian, Li Shunming, Zhang Yuanyuan, et al. Research on stimulus-response test of flat blade[J]. Research and Exploration in Laboratory, 2011,30(1):14-17. (in Chinese)
[3] 陳茉莉,李舜酩,溫衛東,等.多源拍振分析方法與實驗[J].振動、測試與診斷,2011,31(2):202-206.
Chen Moli, Li Shunming, Wen Weidong, et al. Analysis and experiment on multi-source beat vibration[J]. Journal of Vibration, Measurement & Diagnosis, 2011,31(2):202-206. (in Chinese)
[4] Ivica K, Damir S. Wavelet packets approach to blind separation of statistically dependent sources[J]. Neuro Computing, 2008,7(1):1642-1655.
[5] Araki S, Makino S, Aichner R, et al. Subband-based blind separation for convolutive mixtures of speech[J]. Trans Funament, 2005,2(88):3593-3603.
[6] 李寧.頻率域盲信號分離理論研究[D].武漢:華中科技大學,2007.
[7] Lee D D, Seung H S. Alg or ithms for non-negative matrix factorization advances in neural information processing[J]. Advances in Neural Information Processing Systems, 2001,13:556-562.
[8] Li Y, Cichocki A, Amari S. Analysis of sparse representation and blind source separation[J]. Neural Computation, 2004,16(6):1193-1234.
[9] 楊世錫,焦衛東,吳昭同.應用JADE盲分離算法分離統計相關源[J].振動工程學報,2003,16(4):498-501.
Yang Shixi, Jiao Weidong, Wu Zhaotong. Application of JADE to separation of statistically correlated sources[J]. Journal of Vibration Engineering, 2003,16(4):498-501. (in Chinese)
[10]Mallat S G. A theory for multiresolution signal decomposition: the wavelet representation[J]. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 1989,11(7):674-693.
[11]Graps A. An introduction to wavelets[J]. Computational Science & Engineering, 1995,2(2):50-61.
[12]Erkki O, Yuan Zhijian. The fast ICA algorithm revisited: convergence analyis[J]. Neural Networks, 2006,17(6):1370-1381.

*航空自然科學基金資助項目(2012ZD52054);國家自然科學基金資助項目(61403193);南航基本科研業務費科研資助項目(NS2014081)
2013-01-28;
2013-03-18
10.16450/j.cnki.issn.1004-6801.2015.01.025
TH132.4
李紀永,男,1985年1月生,博士研究生。主要研究方向為旋轉機械故障信號處理。 E-mail:ljynav@163.com