一階線性系統隨機共振在轉子軸故障診斷中的應用研究

冷永剛，田祥友

（1．天津大學 機械工程學院，天津 300072；2．天津大學 機構理論與裝備設計教育部重點實驗室，天津 300072）

大型汽輪機、水輪機、燃氣輪機等旋轉動力機械常工作在高溫、高壓和高轉速等較惡劣的工作環境下，其正常運行對于保證工廠的安全生產意義重大。因此需要對旋轉機械進行狀態監測，提取振動信號并進行分析以判別是否存在早期微弱故障。但在工程實際中，信號測點往往離故障源較遠，使測得的信號容易被強噪聲污染，加上早期故障信號本身就比較微弱，這些因素使轉子軸的早期微弱故障診斷較為困難。

目前常用的轉子軸類故障診斷方法有：基于諧波窗、EEMD濾波的轉子軸心軌跡提純方法［1－2］，利用小波或小波包對故障信號進行分解和重構［3］方法，以及基于非線性系統隨機共振檢測微弱信號的振動分析等［4－6］。這些檢測方法各有優點，在一定程度上能反映出轉子軸的故障特征，但同時也存在一些局限性。如轉子軸心軌跡提純法步驟較繁瑣，測試數據時對硬件設備要求較高；小波或小波包降噪時如何選擇閾值和進行閾值優化比較困難；利用非線性系統的隨機共振原理降噪時存在系統參數調節等復雜問題。本文以一階線性系統隨機共振方法為基礎，通過模型簡化，提出一種基于線性系統調參廣義隨機共振的轉子軸類早期微弱故障診斷方法。

1 一階線性系統的調參廣義隨機共振現象

很多學者［7－10］關于線性系統隨機共振現象的研究都是通過復雜的數學模型推導與仿真得出系統輸出存在傳統的或是廣義的隨機共振現象，但這樣的系統中信號與噪聲的作用方式比較復雜，很難應用于實際故障信號的處理分析中。本文僅針對周期信號加高斯白噪聲作用的一階線性系統中出現的廣義隨機共振現象［11］能否用于轉子軸類早期微弱故障信號的檢測進行討論。

一階線性系統動力學模型

式中：a是系統參數，輸入信號 sn（t）＝A sin（2πf0t）＋ξ（t），A sin（2πf0t）是頻率為 f0、幅值為 A、相位為 0的正弦信號，ξ（t）是功率譜密度為2D的噪聲，其中ξ（t）是均值為0，方差為1的高斯白噪聲。

求解系統模型（1），得系統輸出功率譜密度函數為

為比較系統輸出端信號與噪聲的大小關系，引入輸出信噪比。這里信噪比的定義采用信號處理及通信等工程實際中常用的信噪比定義，即系統輸出端信號總功率與噪聲總功率之比［12－13］

式中：S（f）為輸出端信號的功率譜密度函數，SN（f）為輸出端噪聲功率譜密度函數。

由輸出信噪比表達式（4）可以看出，系統參數a不變時，SNRout隨噪聲強度D單調遞減變化，表明周期信號與加性噪聲驅動的一階線性系統中不會出現傳統意義上的隨機共振現象。但仔細觀察式（4）可知，分母中時等號成立），故在噪聲強度D不變的情況下，系統參數a取2πf0時，系統輸出信噪比取得最大值SNRmax＝當輸入信號的頻率與系統參數取得某種協同匹配關系時，系統輸出能夠達到共振狀態。這種輸出信噪比隨系統參數非單調變化的現象稱之為調參廣義隨機共振［11］。

下面通過數值模擬驗證系統的上述特性，設定一組仿真參數：A＝0．1，f0＝0．01 Hz，采樣頻率 fs＝5 Hz，計算數據點數N＝4 000。用四階Runge-Kutta方法（以下同）對方程（1）進行數值計算，觀察輸出信噪比分別隨參數D、a的變化趨勢，計算結果如圖1。算法中的輸出信噪比定義為

式（5）是給定計算數據點數條件下的信噪比，其中PS表示輸出功率譜中信號總功率，PN表示輸出功率譜中噪聲總功率。

圖1 輸出與輸入信噪比分別隨參數D、a的變化趨勢Fig．1 The variation tendency of SNRout and SNRin by the parameter D＆a．Parameter

同樣由圖1仿真曲線可以看出，線性系統的輸出信噪比隨輸入噪聲強度D單調遞減變化，隨系統本身參數a非單調變化，存在a值（a＝0．064，約等于2πf0）使輸出信噪比達到極值。圖中虛線表示相應情況下的輸入信噪比，其定義為在給定數值分析長度條件下，信號總功率與噪聲總功率之比。

2 一階線性系統調參廣義隨機共振的微弱信號檢測

由前述分析知，在給定輸入信號與噪聲時，系統取參數a＝2πf0輸出信噪比達到極大值。為了滿足實際工程應用，文獻［11］分析了系統參數a＝2πf0固定不變情況下，采樣頻率與信號頻率比值對輸出頻譜的影響，并給出了檢測特征信號可辨識性好的比值范圍。這表明，實際應用中，使系統輸出信噪比最大的參數a值并不一定就能使輸出特征信號達到理想的可辨識性，因為信噪比大，只說明系統輸出的信號能量與噪聲能量比值大，而特征信號的可辨識性與信號幅值、噪聲大小以及噪聲的分布有關。因此需要進一步分析，實際應用時應如何調節a值才能使特征信號的辨別性最優。下面仿真說明高斯白噪聲條件下，輸出頻譜上特征信號可辨識性與參數a的變化關系。

根據文獻［11］比較輸出頻譜圖上特征信號可辨識程度的方法，即提取輸出譜圖上信號頻率處的峰值h，將其與整個譜圖上除去信號頻率處譜峰后其余譜值中最大的一個值hr作比較，比值越大表明特征信號的可辨識性越好，并將這個比值定義為辨別率r＝h／hr．，一般只有r＞1時特征信號才容易識別出來。在采樣頻率與信號頻率比值不變的情況下，改變參數a值進行仿真，計算不同a值下輸出譜圖上特征信號的辨別率r。仿真參數同上：A＝0．1，f0＝0．01 Hz，D＝0．8，采樣頻率fs＝5 Hz，計算數據點數 N＝4 000，令參數 a與2πf0的比值在0．1～10之間變化，結果如圖2所示。

由圖2可以看出，保持采樣頻率與信號頻率比值fs／f0＝500不變的情況下，輸出譜圖特征信號的辨別率r與參數a的取值有很大關系，隨a／2πf0比值的變化趨勢是先增大至某一極大值，并幾乎保持一段水平值，然后再減小，當a／2πf0＝4～9時存在一個水平極大值。由此可見選取一個合適的參數a值對特征信號頻率的判別有著直接的影響，實際信號檢測時應盡量選取使辨別率達到最佳的a值。圖3比較了參數a取兩個不同值時輸出譜圖中特征信號的可辨識性，其余參數同圖2。

圖2 系統輸出譜圖上特征信號辨別率r隨 a／2πf0變化曲線Fig．2 The characteristic signal recognition r changes with a／2πf0 in system output spectrum

由圖3可見，圖（b）比圖（a）的特征信號辨別性明顯要好。雖然圖（a）中f0＝0．01 Hz處的峰值已經是整個譜圖上的最大值，但與其附近峰值高度相接近，容易造成誤判，而圖（b）在f0＝0．01 Hz處的譜值明顯高于附近其他峰值，很容易確定這就是所尋找的特征信號頻率成分。

圖3 不同參數a值下系統輸出頻譜Fig．3 The output spectrum of different parameter a

這里需要指出的是，圖2得到的辨別率r隨比值a／2πf0變化關系是在采樣頻率fs與信號頻率f0比值保持恒定的情況下得到的，如果改變采樣頻率與信號頻率的比值，那么r隨比值a／2πf0的變化關系曲線也會發生變化。為綜合分析r隨比值a／2πf0和fs／f0的變化關系，圖4給出了三者變化關系的三維圖。從圖4可見，對應于每一組較優的fs／f0值，均可以找到合適的a／2πf0值，使特征信號辨別率r達到最大。另外，應用此方法處理含噪微弱信號時，可不必局限在小頻率范圍，只要采樣頻率合適，可以用于檢測任意頻率大小信號［11］。

圖4 特征信號辨別率r隨比值a／2πf0和 fs／f0變化三維圖Fig．4 The 3D graph of characteristic signal recognition r changes with a／2πf0 and fs／f0

根據以上分析，應用線性系統調參廣義隨機共振檢測實際微弱信號時的步驟可總結為：① 根據已知條件估算故障信號特征頻率f0的可能范圍，參考特征信號辨別率隨比值 a／2πf0和 fs／f0變化的三維圖，確定合適的采樣頻率fs和參數a使信號可識別性最佳。②若輸出譜圖上信號特征頻率不明顯，可通過微調參數a的值進行修正或重新確定一組采樣頻率fs和參數a采樣，再次識別計算。

3 轉子軸早期微弱故障識別分析

在旋轉機械中軸彎曲是一種常見故障，轉子軸彎曲是指轉軸各截面幾何中心連線與旋轉軸線不重合，從而使轉子產生偏心質量，進而引起不平衡振動。實驗在滑動軸承故障試驗臺上進行，轉子軸中心偏離旋轉軸線0．38 mm，在試驗臺遠離軸承座向外0．5 m處布置一個加速度傳感器，一方面模擬工程實際中由于空間限制無法在軸承座附近布置傳感器采集信號的情況，另一方面使故障信號強度進一步衰減，模擬早期微弱故障。

為獲得這一早期微弱故障特征信號，利用一階線性系統調參共振方法進行處理時，首先根據軸的轉頻和故障類型，判定彎曲故障信號頻率f0值，實驗中轉子軸的轉速約為1 700 r／min，則故障信號頻率 f0約在28 Hz左右。利用NI數據采集儀和Lab VIEW采集軟件進行數據采集，依據特征信號可辨識性最好的采樣頻率取值原則［11］，設置采樣頻率為 fs＝1 kHz，采樣點數N＝5 000，實測信號的時域波形及低頻頻譜如圖5所示。可知，信號譜干擾成分較多，很難通過直接觀察的方法從頻譜圖上準確判斷出故障特征。參考圖4取線性系統參數 a＝1．1×2πf0，將實測信號輸入式（1）系統進行處理，得到圖6時頻域結果。

圖5 實測轉子軸彎曲故障信號的時域圖及頻域圖Fig．5 The time domain and frequency domain of observed rotor shaft bending fault signal

圖6 轉子軸彎曲故障信號經過線性系統處理后的時域及頻域圖Fig．6 The time domain and frequency domain of rotor shaft bending fault signal after being processed in linear system

圖7 轉子軸故障信號經過非線性雙穩系統處理后的時域及頻域圖（參數 a＝1，b＝10）Fig．7 The time domain and frequency domain of rotor shaft fault signal after being processed in non-linear bistable system．（parameter a＝1，b＝10）

由圖6可以看到，經過處理后輸出頻譜圖上27．8 Hz處有一明顯譜峰，且遠高于周圍其他譜峰，由此識別出滑動軸承轉子存在“早期微弱”彎曲不平衡故障。

作為比較，以下采用非線性雙穩隨機共振方法處理該微弱故障信號［5－6］，非線性雙穩隨機共振檢測模型是

式中：a，b為系統參數，sn（t）為待測含噪信號。將上述實驗采集數據代入模型（6）進行處理。由于圖5顯示信號幅值很低，能量不足以引起粒子在兩勢阱間的躍遷，所以首先對采集數據進行幅值放大處理，將原始含噪信號幅值同比增大 200倍，即 sn′（t）＝200×sn（t）。其次，根據產生隨機共振的小參數特性，按照文獻［5－6］變尺度方法對原始數據進行二次采樣，取二次采樣頻率 fsr＝5 Hz。令參數 a＝1，b＝10，將經過上述處理的數據代入式（6），用四階Runge-Kutta法計算，并進行尺度還原得到系統輸出時域圖和低頻頻譜如圖7所示。

由圖7可以看出，在相應參數條件下系統為欠共振狀態，頻譜圖上不易看出明顯的特征信號。嘗試改變系統參數a＝1，b＝50，重新計算得到系統輸出時域圖和低頻頻譜如圖8所示。經過參數的再次調節，圖8的頻譜可以分辨出故障特征信號成分。

圖8 轉子軸故障信號經過非線性雙穩系統處理后的時域及頻域圖（參數a＝1，b＝50）Fig．8 The time domain and frequency domain of rotor shaft fault signal after being processed in non-linear bistable system．（parameter a＝1，b＝50）

從兩種檢測方法的對比中可以看出，利用非線性系統隨機共振法檢測微弱信號時，多參數的協調選擇是一個難點，這些參數調節一般需要根據經驗選取，多次嘗試后才能得到較理想的結果。相比較用一階線性系統的調參隨機共振法檢測，采樣頻率選定后只需要調節系統參數a，就可以得到有效結果，因此該方法在實際應用上更簡便有效。

4 結 論

本文研究了一階線性系統調參廣義隨機共振機理，及其用于檢測微弱故障信號的處理方法。周期信號與高斯白噪聲作用下的一階線性系統，通過調節系統參數，在輸出端信噪比能取到極大值，但這并不能保證輸出譜圖上特征信號能達到理想的識別性。為在輸出譜圖上得到更清晰的特征信號，文中討論了信號特征頻率的可辨識性與系統參數、信號頻率、采樣頻率等參數之間的變化關系，給出實際應用時各參數取值的參考圖。隨后滑動軸承試驗臺上轉子軸彎曲故障的實例分析，驗證了一階線性系統調參共振檢測早期微弱故障信號的有效性與簡便性。

［1］張文斌，周曉軍，楊先勇，等．基于諧波窗方法的轉子軸心軌跡提純［J］．振動與沖擊，2009，28（8）：74－77．ZHANG Wen-bin，ZHOU Xiao-jun，YANG Xian-yong，et al．Harmonic window method for purification of axis trace［J］．Journal of Vibration and Shock，2009，28（8）：74－77．

［2］陳仁祥，湯寶平，呂中亮．EEMD濾波的轉子軸心軌跡提純方法［J］．重慶大學學報，2012，35（11）：15－20．CHEN Ren-xiang， TANG Bao-ping， L Zhong-liang． A mthod of rotor orbit purification based on ensemble empirical mode decomposition filter［J］． Journal of Chongqing University，2012，35（11）：15－20．

［3］楊文志，馬文生，任學平．小波包降噪方法在滑動軸承故障診斷中的應用研究［J］．噪聲與振動控制，2009，29（4）：50－53．YANG Wen-zhi， MA Wen-sheng， REN Xue-ping．Application of wavelet packet denoising method in fault diagnosis of bearings［J］．Noise and Vibration Control，2009，29（4）：50－53．

［4］石鵬，冷永剛，范勝波，等．雙穩系統處理微弱沖擊信號的研究［J］．振動與沖擊，2012，31（6）：150－154．SHI Peng，LENG Yong-gang，FAN Sheng-bo，et al．A bistable system for detecting a weak pulse signal［J］．Journal of Vibration and Shock，2012，31（6）：150－154．

［5］Leng Y G，Wang T Y，Guo Y，et al．Engineering signal processing based on bistable stochastic resonance［J］．Mechanical Systems and Signal Processing，2007，21：138－150．

［6］冷永剛．雙穩調參高頻共振機理［J］．物理學報，2011，60（2）：020503＿1－7．LENG Yong-gang．Mechanism of high frequency resonance of parameter-adjusted bistable system ［J］． Acta Physica Sinica，2011，60（2）：020503＿1－7．

［7］Berdichevsky V，Gitterman M．Stochastic resonance in linear systems subject to multiplicative and additive noise［J］．Phys．Rev．E，1999，60（2）：1494－1499．

［8］Gitterman M．Harmonic oscillator with multiplicative noise：Nonmonotonic dependence on the strengthand the rate of dichotomous noise［J］．Phys．Rev．E，2003，67：057103．

［9］張莉，劉立，曹力．過阻尼諧振子的隨機共振［J］．物理學報，2010，59（3）：1494－1498．ZHANG Li，LIU Li，CAO Li．Stochastic resonance in an overdamped harmonic oscillator［J］．Acta Physica Sinica，2010，59（3）：1494－1498．

［10］陸志新，曹力．輸入方波信號的過阻尼諧振子的隨機共振［J］．物理學報，2011，60（11）：110501．LU Zhi-xin，CAO Li．Stochastic resonance of square wave signal in an overdamped harmonic oscillator［J］．Acta Physica Sinica，2011，60（11）：110501．

［11］田祥友，冷永剛，范勝波．一階線性系統的調參隨機共振研究［J］．物理學報，2013，62（2）：020505．TIAN Xiang-you， LENG Yong-gang， FAN Sheng-bo．Parameter-adjusted stochastic resonance of first-order linear system［J］．Acta Physica Sinica，2013，62（2）：020505．

［12］Mitaim S，Kosko B．Adaptive stochastic resonance［J］．Proceedings of the IEEE，1998，86：2152－2183．

［13］ Gingl Z，Vajtai R，Kiss L B．Signal-to-noise ratio gain by stochastic resonance in a bistable system［J］．Chaos，Solitons＆Fractals，2000，11：1929－1932．