唐貴基，孫軍科，，王曉龍，伍小林，周福成，崔彥亭，吳 韜，胥佳瑞

（1. 華北電力大學(xué) 機械工程系 河北省電力機械裝備健康維護與失效預(yù)防重點實驗室，河北 保定 071003；2. 中國大唐集團科學(xué)技術(shù)研究院有限公司 火力發(fā)電技術(shù)研究院，北京 100040）

0 引言

我國電力工業(yè)發(fā)展迅速，根據(jù)中國電力企業(yè)聯(lián)合會發(fā)布的《2019—2020年度全國電力供需形勢分析預(yù)測報告》［1］統(tǒng)計結(jié)果，2019年全年全社會的用電量為7.23×1012kW·h，相比上年增長了4.5%。截至2019年底，全國全口徑發(fā)電裝機容量為2.01×109kW，其中火電裝機容量為1.19×109kW，水電裝機容量為3.6×108kW，核電裝機容量為4.874×107kW，并網(wǎng)風(fēng)電裝機容量為2.1×108kW，并網(wǎng)太陽能發(fā)電裝機容量為2.0×108kW，占比分別為59.2%、17.9%、2.4%、10.4%、9.9%，由此可見在未來一段時間內(nèi)，火力發(fā)電仍將占據(jù)我國電力供應(yīng)的主導(dǎo)地位。汽輪發(fā)電機組、風(fēng)機、給水泵等旋轉(zhuǎn)機械作為火力發(fā)電廠的主要設(shè)備對電力生產(chǎn)的重要性不言而喻，其轉(zhuǎn)子更是重中之重，轉(zhuǎn)子承受蒸汽對工作葉片的回轉(zhuǎn)力并作為重要的負載和傳動結(jié)構(gòu)，它的運行狀態(tài)將直接影響整個發(fā)電系統(tǒng)，因此，對轉(zhuǎn)子進行故障診斷具有重要的現(xiàn)實意義與經(jīng)濟價值［2］。

常見的轉(zhuǎn)子故障類型包括動靜碰摩、油膜渦動、油膜振蕩、裂紋、不平衡、不對中等［3-4］。大量研究表明轉(zhuǎn)子故障振動信號通常為非平穩(wěn)多分量耦合信號［5-7］，由于時頻分析方法能同時提供振動信號的時域和頻域信息，因而在轉(zhuǎn)子故障診斷領(lǐng)域中得到了廣泛的應(yīng)用，探究行之有效的時頻分析方法成為國內(nèi)外學(xué)者的研究熱點。希爾伯特-黃變換HHT（Hilbert-Huang Transform）是近年來信號處理領(lǐng)域常用的時頻分析方法之一，它由經(jīng)驗?zāi)B(tài)分解EMD（Empirical Mode Decomposition）和希爾伯特變換HT（Hilbert Transform）兩部分組成［8-9］。文獻［10］利用HHT 清晰地表征了橫向裂紋轉(zhuǎn)子的扭振所形成的相位調(diào)制現(xiàn)象，為轉(zhuǎn)子裂紋故障的預(yù)測和診斷提供了有效的方法。文獻［11］基于HHT 的思想和Teager 能量算子在一定程度上改善了HT 存在的端點效應(yīng)，提出了Teager-Huang 時頻分析方法用于轉(zhuǎn)子碰摩故障診斷，但所使用的信號分解方法為EMD，其存在模態(tài)混疊、包絡(luò)擬合時易出現(xiàn)端點效應(yīng)等問題，從而影響分析效果［12］。文獻［13］針對局部特征尺度分解LCD（Local Characteristic-scale Decomposition）方法的分解精度不足這一缺陷，提出了基于分段多項式的局部特征尺度分解PPLCD（Piecewise Polynomial based Local Characteristic-scale Decomposition）方法，并成功應(yīng)用于轉(zhuǎn)子碰摩故障診斷。文獻［14］針對轉(zhuǎn)子疲勞裂紋、動靜碰摩和轉(zhuǎn)子-軸承系統(tǒng)不平衡3 類故障振動信號特征提出了基于局部均值分解LMD（Local Mean Decomposition）的時頻分析方法。文獻［15］將互相垂直方向上的振動信號組成復(fù)數(shù)信號，利用復(fù)局部均值分解CLMD（Complex Local Mean Decomposition）對合成復(fù)數(shù)信號進行自適應(yīng)分解，實驗結(jié)果表明該方法可準(zhǔn)確診斷滾動軸承微弱故障和復(fù)合故障。EMD、LCD 和LMD 方法的基本原理均是基于極值點擬合的遞歸分離方法，當(dāng)待分解的信號中存在噪聲、間斷信號等異常事件干擾時，其分解所得的本征模態(tài)函數(shù)IMF（Intrinsic Mode Function）分量波形易發(fā)生畸變，從而出現(xiàn)模態(tài)混疊現(xiàn)象，如何避免或解決模態(tài)混疊現(xiàn)象是目前信號處理領(lǐng)域的一個重要研究方向。

P. Bonizzi 等人［16］于2014 年提出了一種以奇異譜分析為基礎(chǔ)的信號分析處理方法——奇異譜分解（SSD）方法。SSD 方法不同于EMD 方法通過極值點擬合實現(xiàn)信號分量的提取，其核心思想是通過奇異值分解從頻帶的劃分和篩選角度去實現(xiàn)特定信號分量的重構(gòu)，此重構(gòu)方式使得SSD 能較好地抑制偽分量的產(chǎn)生和模態(tài)混疊現(xiàn)象的出現(xiàn)，表現(xiàn)出了較好的魯棒性。但該方法以迭代殘余分量與原始信號的能量比作為迭代停止條件，不具有穩(wěn)定性，且針對轉(zhuǎn)子類振動信號，該方法在分解過程中的殘余分量與原始信號的能量比不易控制。若能量比閾值設(shè)置過大，則易導(dǎo)致信號分解不充分，出現(xiàn)欠分解的現(xiàn)象，掩蓋包含故障特征信息的分量成分；若能量比閾值設(shè)置過小，則可能會出現(xiàn)過分解的現(xiàn)象，產(chǎn)生過多的虛假分量，導(dǎo)致計算效率低，且會給后續(xù)分析帶來極大的不便。針對該缺陷，本文提出了融合互信息判據(jù)和能量比的循環(huán)迭代協(xié)同控制策略，對SSD 方法進行改進，并與Teager能量算子解調(diào)算法相結(jié)合，提出了基于改進SSD-Teager時頻分析的轉(zhuǎn)子故障診斷方法，并應(yīng)用仿真信號和現(xiàn)場實測信號對該方法的有效性及優(yōu)越性進行驗證。

1 基于改進SSD-Teager時頻分析的轉(zhuǎn)子故障診斷方法

1.1 SSD方法基本原理

SSD 方法通過迭代循環(huán)可將非線性非平穩(wěn)信號分解為按由高頻到低頻順序排列的奇異譜分量SSC（Singular Spectrum Component），具體分解步驟如下［16］。

矩陣X左側(cè)區(qū)域的3行3列部分為軌跡矩陣，將X的右下角3個元素移至矩陣左上角位置，從而得到式（2）所示的新矩陣。

新矩陣中各對角線元素相同且數(shù)量相等，其中矩陣左側(cè)區(qū)域為所構(gòu)建的新軌跡矩陣Y，該構(gòu)造方式增強了原始信號中的振蕩成分，并使得迭代后殘余分量的能量呈現(xiàn)遞減規(guī)律。

2）自適應(yīng)選擇嵌入維數(shù)M。

首先由式（3）計算第j次迭代運算時的殘余分量υj(n)：

然后再估計功率譜密度最大峰值處對應(yīng)的頻率fmax。在第1 次迭代過程中，若fmax/fs（fs為采樣頻率）小于給定閾值，則可將殘余分量視為一個大趨勢項，此時設(shè)定嵌入維數(shù)M=N/3；否則，在j ＞1 的情況下，嵌入維數(shù)M=1.2fs/fmax。

3）重構(gòu)第j個SSC。若在第1 次迭代過程中就檢測到大趨勢項，則計算矩陣X1=σ1u1vT1（其中，σ1為所構(gòu)造的第1 個高斯函數(shù)的帶寬，u1為高斯函數(shù)的譜模型擬合的功率譜密度輪廓，v1為第1 次迭代運算時的殘余分量），對X1進行對角平均得到分量序列時間尺度g()1(n)；否則執(zhí)行下一次j ＞1 時的循環(huán)，使用左特征向量在頻譜[fmax-Δf，fmax+Δf]（Δf為殘余項功率譜密度中主峰的半帶寬）范圍內(nèi)具有突出主頻率的所有特征組和對選取模態(tài)分量的主峰能量貢獻最大的一個特征組，創(chuàng)建子集Ij={i1，i2，…，ip}（is（s=1，2，…，p）為所提取的特征組，p為所有特征組的數(shù)量），然后通過矩陣XIj=Xi1+Xi2+…+Xip的對角平均來獲得相應(yīng)的SSC。

4）迭代停止條件。在分解過程中得到分量序列時間尺度g(j)(n)的同時會產(chǎn)生一個殘余信號分量υ(j+1)(n)=υ(j)(n)-g(j)(n)，由式（4）計算殘余分量與原始信號的歸一化均方值誤差λ(j)NMSE。

式中：m為SSD得到的SSC數(shù)量。

1.2 改進SSD方法

由SSD 方法原理可知，處理信號時其迭代停止條件能量比閾值λth需預(yù)先設(shè)定，該值設(shè)定合理與否將直接影響分解效果以及SSC 數(shù)量，然而受不同設(shè)備、不同工況的影響，λth值通常難以準(zhǔn)確設(shè)定。針對這一問題，本文提出一種融合互信息MI（Mutual Information）判據(jù)和能量比的循環(huán)迭代協(xié)同控制策略，對SSD 的信號處理過程進行改進。該策略的核心思想是通過互信息計算結(jié)果甄別SSD 過程是否產(chǎn)生虛假分量，從而實現(xiàn)循環(huán)迭代過程的有效控制，同時為避免由于無法滿足互信息判定閾值下限而陷入無限死循環(huán)，利用能量比指標(biāo)作為輔助控制條件來強制結(jié)束迭代過程。

互信息引申于信息論中熵的概念，用于判斷2個隨機變量A與B之間的相關(guān)性，能夠定量地表示二者的相互依賴程度，比相關(guān)系數(shù)和互相關(guān)函數(shù)法更能表征不同變量間的相關(guān)程度［17］，互信息表達式為：

對互信息進行歸一化處理后，將其用于判斷分解殘余分量與原始信號的相關(guān)性，當(dāng)最小歸一化互信息值低于設(shè)定閾值時，認(rèn)為殘余分量不再蘊含故障特征信息，原始信號已被充分分解，從而結(jié)束整個運算過程，即可根據(jù)實際分析信號自適應(yīng)地控制SSD 方法的迭代過程。參考相關(guān)文獻［17-19］中的閾值設(shè)定范圍，并通過多組不同類型仿真信號進行分析驗證，綜合對比發(fā)現(xiàn)，當(dāng)互信息判定閾值為0.05時，改進SSD 方法普適性最佳，因此本文中設(shè)置互信息判定閾值為0.05?；诨バ畔⑴袚?jù)和能量比協(xié)同控制策略的改進SSD 方法流程圖見附錄A 圖A1，具體實現(xiàn)步驟如下：

1）初始化i=1，設(shè)置SSD 方法的迭代停止條件能量比閾值=0.1，搜索范圍從0.1 開始至0.001 結(jié)束，搜索步長為0.001；

2）設(shè)定第i次迭代分解的輸入信號xi(t)，加載原始振動信號x(t)，x1(t)=x(t)，通過自適應(yīng)選擇嵌入維數(shù)M構(gòu)造新矩陣X，進而獲取新軌跡矩陣Y，對Y進行SSD 重組，將信號xi(t)分解為mi個SSC（j=1，2，…，mi）和1個殘余分量ri(t)；

6）結(jié)束整個循環(huán)輸出結(jié)果，得到n個SSC。

1.3 改進SSD-Teager時頻分析方法

Teager 能量算子是一種非線性算子，其通過對振動信號的瞬時值及其微分進行非線性組合放大信號的瞬時能量成分，從而增強信號的沖擊特征，適用于檢測信號中的沖擊成分，且Teager 能量算子的包絡(luò)解調(diào)效果相比HT 具有更高的時間分辨率［20-21］?；诖?，本文將Teager能量算子與改進SSD方法相結(jié)合，提出了一種基于改進SSD-Teager時頻分析方法。

對于連續(xù)時間信號x(t)，Teager 能量算子ψ的定義為［22］：

離散信號x(n)的Teager 能量算子ψ[x(n)]定義為：

由式（9）可知，任一時刻Teager能量算子的計算僅需連續(xù)3 個樣本數(shù)據(jù)即可，計算量較小且計算效率高。

利用Teager 能量算子進行解調(diào)計算，通過式（10）、（11）可以分別得到信號x(t)的瞬時頻率、瞬時幅值。

由式（10）和式（11）可了解信號的時頻分布信息，有利于信號特征的準(zhǔn)確辨識，因此本文將改進SSD 方法與Teager 能量算子解調(diào)算法相結(jié)合，利用改進SSD 方法對多分量信號進行處理，自適應(yīng)得到若干個SSC 后，分別計算各SSC 的瞬時幅值和瞬時頻率，并將所有SSC 的時頻分布特征進行集合匯總，繪制得到多分量信號的Teager 時頻譜，從而獲得其完整的時頻分布特征信息。

1.4 故障診斷流程

改進SSD-Teager 時頻分析方法用于轉(zhuǎn)子故障診斷時的具體流程為：

1）利用改進SSD 方法對傳感器采集的引風(fēng)機轉(zhuǎn)子故障振動信號進行分解處理，得到一系列按由高頻到低頻順序排列的SSC；

2）利用Teager 能量算子解調(diào)求解每個SSC 的瞬時幅值和瞬時頻率，得到各個分量的時頻表示；

3）將所有SSC 的時頻分布進行匯總，得到對應(yīng)的時頻圖，并根據(jù)所得信號的整體時頻特征判斷轉(zhuǎn)子故障狀態(tài)。

2 仿真信號分析

為驗證本文所提方法的分解效果和時頻分析能力，構(gòu)造一含有隨機噪聲的多分量信號x(t)，該信號由1個低頻諧波分量x1(t)、2個頻率成分不同的調(diào)頻信號x2(t)和x3(t)以及隨機噪聲x4(t)疊加而成，具體表達式為：

式中：t∈[0，0.25]s；分析點數(shù)N=1024；randn(1，N)表示生成1 行N列服從均值為0、方差為1 的正態(tài)分布的隨機數(shù)；length(t)為時間t內(nèi)的數(shù)據(jù)長度。合成信號x(t)及各分量信號的時域波形如圖1 所示，信號采樣頻率fs=4096 Hz。

圖2 是利用本文所提改進SSD 方法對仿真信號x(t)進行分析后所得分解結(jié)果，共得到3 個SSC。由圖可見，SSC1、SSC2、SSC3分別對應(yīng)x3(t)、x2(t)、x1(t)這3 個信號分量，可見本文所提改進SSD 方法能準(zhǔn)確地將信號分量按照由高頻到低頻的順序排列，同時能自適應(yīng)地將隨機噪聲信號分量x4(t)視為殘余分量不納入分解結(jié)果中。對改進SSD 方法所得分量進行Teager能量算子解調(diào)得到的時頻分布圖見附錄A 圖A2，由圖可清晰地看出各分量頻率成分以及頻率調(diào)制現(xiàn)象。

同時為了證明本文所提改進SSD方法的優(yōu)越性，將SSD方法和HHT方法與其進行對比。

對仿真信號直接進行SSD，所得結(jié)果如圖3 所示。由圖可見，對仿真信號直接進行SSD 共得到5個SSC，這說明SSD 方法出現(xiàn)了過分解的現(xiàn)象，產(chǎn)生了冗余分量，而本文所提改進SSD方法準(zhǔn)確分解出3個有效SSC，自適應(yīng)剔除添加的隨機噪聲成分，且未分解出多余分量。對SSD 方法分解所得分量進行Teager 能量算子解調(diào)，得到的二維時頻分布圖如附錄A 圖A3所示。從圖中的箭頭指示可以看出，過分解現(xiàn)象的存在致使圖中出現(xiàn)了干擾頻率成分，對比圖A2 可以發(fā)現(xiàn)原始SSD-Teager 時頻分析方法的分析效果較本文所提改進SSD-Teager 時頻分析方法差。

圖3 仿真信號的SSD結(jié)果Fig.3 SSD result of simulated signal

利用HHT 方法對仿真信號進行處理，仿真信號經(jīng)EMD 處理后共得到7 個IMF，如圖4 所示。通過分 析 發(fā) 現(xiàn)IMF2、IMF3和IMF4這3 個 分 量 分 別 對 應(yīng)x3(t)、x2(t)和x1(t)，但如圖中虛線標(biāo)記所示，IMF2—IMF4均出現(xiàn)了較大程度的模態(tài)混疊和幅值失真現(xiàn)象，IMF1、IMF5、IMF6和IMF7的存在表明出現(xiàn)了過分解現(xiàn)象。附錄A 圖A4 為對應(yīng)的2 維時頻分布圖，由圖可見由于EMD 方法的不足致使最終的時頻分析結(jié)果較差。

圖4 仿真信號的EMD結(jié)果Fig.4 EMD result of simulated signal

3 工程案例分析

某電廠4 號機組配有2 臺由成都風(fēng)機廠生產(chǎn)制造的雙級動葉可調(diào)軸流式引風(fēng)機，該引風(fēng)機型號為HU27050-22，主要由進氣室、集流室、導(dǎo)葉、擴散器、動葉調(diào)節(jié)驅(qū)動裝置、伺服驅(qū)動裝置等部件組成，引風(fēng)機水平布置，垂直進氣、水平出氣。引風(fēng)機驅(qū)動裝置為湘潭電機廠制造的YXKS1000-8 型電機，其額定功率為6 900 kW，額定轉(zhuǎn)速為745 r／min。

電機-引風(fēng)機軸系結(jié)構(gòu)示意簡圖如圖5 所示，其軸系共有1、2 號2 個單軸承用于支撐電機轉(zhuǎn)子和1個3號組合軸承用于支撐引風(fēng)機轉(zhuǎn)子，2段軸系通過靠背輪由螺栓緊固。鍵相探頭布置在電機側(cè)靠背輪處，電機側(cè)傳感器布置在電機驅(qū)動端2 號軸承蓋的水平和垂直方向上用于測量相應(yīng)方向上的電機振動；引風(fēng)機側(cè)傳感器布置在3 號軸承組前軸承位置，同樣用于測量水平、垂直方向上的電機振動。本次振動測量所使用的傳感器均在有效期內(nèi)，并在使用前均校驗合格，振動數(shù)據(jù)由合肥偉博測控科技有限公司生產(chǎn)的Vib900A 多通道現(xiàn)場振動信號采集系統(tǒng)采集。根據(jù)機組的歷次啟停機振動數(shù)據(jù)分析，在機組升降速過程中存在振動異常情況，采集時電機轉(zhuǎn)速為480 r／min，采樣頻率設(shè)置為1 024 Hz，分析點數(shù)為2048。

圖5 軸系結(jié)構(gòu)示意簡圖Fig.5 Simplified diagram of shafting structure

圖6 為引風(fēng)機側(cè)水平、垂直方向上振動信號的時域波形圖和頻譜圖，圖7 為兩方向振動信號合成的軸心軌跡圖。由圖可見，振動信號時域波形規(guī)律性明顯，呈現(xiàn)出基頻和二倍頻成分疊加的“M”形波形，頻譜圖中存在基頻和二倍頻成分，軸心軌跡大致呈雙環(huán)橢圓狀，初步懷疑軸系存在不對中故障。

圖6 原始振動信號的時域波形和頻譜圖Fig.6 Time-domain waveforms and spectra of original vibration signals

圖7 軸心軌跡圖Fig.7 Diagram of axis orbit

為進一步確定該引風(fēng)機轉(zhuǎn)子故障類型，利用本文所提方法進行進一步分析。垂直方向振動信號的改進SSD 結(jié)果如圖8 所示，由圖可見改進SSD 方法共得到4 個SSC。通過Teager 能量算子解調(diào)方法計算各個SSC的瞬時幅值及瞬時頻率，匯總后得到的3維時頻分布圖如附錄A 圖A5 所示。4 個SSC 的瞬時頻率分別為8、16、24、40 Hz，即分別對應(yīng)基頻、二倍頻、三倍頻以及五倍頻分量，并且二倍頻分量為主要分量，結(jié)合軸心軌跡呈雙環(huán)橢圓以及二倍頻能量明顯超過基頻能量等特征［23］可充分說明該引風(fēng)機轉(zhuǎn)子出現(xiàn)了不對中故障。

圖8 實測信號的改進SSD結(jié)果Fig.8 Improved SSD result of measured signal

圖9 為直接利用SSD 處理原始振動信號得到的分解結(jié)果，雖然也準(zhǔn)確提取出了基頻和二倍分量，但高頻諧波分量（三倍頻、五倍頻分量）卻沒有被分解提取出來，出現(xiàn)了欠分解情況。通過Teager 能量算子解調(diào)計算2 個SSC 的瞬時幅值及瞬時頻率，并匯總后得到3 維時頻分布圖如附錄A 圖A6 所示。2 個SSC 的瞬時頻率分別為8 Hz 和16 Hz，即分別對應(yīng)基頻和二倍頻成分，可見由于SSD 方法出現(xiàn)欠分解現(xiàn)象導(dǎo)致最終的時頻分析效果較本文所提方法差。為進一步證明本文方法的優(yōu)越性，將HHT 方法與本文方法進行對比分析。原始振動信號經(jīng)EMD 處理所得結(jié)果如圖10 所示。由圖可見，EMD 共得到6 個IMF，其中IMF2對應(yīng)工頻成分，IMF1分量對應(yīng)二倍頻成分；由虛線標(biāo)記處可以看出包含主要故障特征信號的2 個IMF 分量均出現(xiàn)了模態(tài)混疊和幅值失真現(xiàn)象。附錄A 圖A7 為對應(yīng)的3 維時頻圖，由圖可見，各IMF 分量頻率調(diào)制現(xiàn)象明顯，頻率成分雜亂無章無法準(zhǔn)確識別具體頻率信息，HHT 方法的時頻分析效果明顯較本文方法差。

圖9 實測信號的SSD結(jié)果Fig.9 SSD result of measured signal

圖10 實測信號的EMD結(jié)果Fig.10 EMD result of measured signal

4 結(jié)論

本文結(jié)合SSD 方法和Teager 能量算子解調(diào)各自的優(yōu)點，提出了一種改進SSD-Teager 時頻分析的引風(fēng)機轉(zhuǎn)子故障診斷方法。針對SSD 方法以能量比作為迭代停止條件導(dǎo)致魯棒性欠佳的缺點，本文融合互信息判據(jù)對循環(huán)迭代過程進行協(xié)同控制，從而自動確定SSC 數(shù)量，實現(xiàn)信號自適應(yīng)處理過程。利用改進SSD 方法對原始振動信號進行自適應(yīng)分解得到一系列SSC，再融合Teager 能量算子解調(diào)求取各分量的瞬時振幅和瞬時頻率信息，最后根據(jù)時頻特征判別故障類型。將改進SSD-Teager時頻分析方法應(yīng)用于現(xiàn)場實測引風(fēng)機轉(zhuǎn)子不對中故障振動信號分析，最終準(zhǔn)確地識別出轉(zhuǎn)子故障狀態(tài)，并且與傳統(tǒng)HHT 方法相比優(yōu)勢明顯。綜上所述本文方法具有一定工程實際應(yīng)用價值，可為類似旋轉(zhuǎn)機械故障損傷判定提供參考借鑒。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。