趙洪山，郭 偉，邵 玲，張興科

（1.華北電力大學(xué) 電氣與電子工程學(xué)院，河北 保定 071003；2.泰州供電公司，江蘇 泰州 225300）

0 引言

風(fēng)力發(fā)電是當(dāng)前發(fā)展最為迅速的一種綠色能源，據(jù)中國風(fēng)能協(xié)會(huì)的數(shù)據(jù)統(tǒng)計(jì)，2012年底我國風(fēng)電并網(wǎng)總裝機(jī)容量60.83 GW，躍居世界第一；發(fā)電量為100.4 TW·h，已經(jīng)超過核電的98.2 TW·h，成為繼煤電和水電之后的第三大主力電源。風(fēng)電場(chǎng)一般坐落在廣闊的偏遠(yuǎn)地區(qū)，受到惡劣的自然環(huán)境因素影響，再加上風(fēng)機(jī)自身制造工藝和技術(shù)發(fā)展的不完善，這些因素都會(huì)增加風(fēng)電機(jī)組發(fā)生故障的風(fēng)險(xiǎn)，且維修時(shí)往往很困難，造成風(fēng)電場(chǎng)后期的運(yùn)行維修費(fèi)用居高不下。據(jù)統(tǒng)計(jì)，對(duì)于一臺(tái)設(shè)計(jì)壽命為20 a的750 kW風(fēng)機(jī)，它的運(yùn)行和維修費(fèi)用將占整個(gè)發(fā)電成本的 25%～30%，占其投資費(fèi)用的 75%～90%[1]。

作為風(fēng)機(jī)傳動(dòng)系統(tǒng)的關(guān)鍵部件，齒輪箱的任務(wù)是將風(fēng)輪在風(fēng)力作用下所產(chǎn)生的動(dòng)力傳遞給發(fā)電機(jī)。雖然風(fēng)機(jī)齒輪箱的制造工藝已較為成熟，故障率不高，然而一旦故障其修復(fù)過程很復(fù)雜，是造成風(fēng)電機(jī)組停機(jī)時(shí)間最長(zhǎng)的故障之一[2-3]。齒輪箱的故障一般是由輪齒損壞和軸承磨損造成的，在其故障發(fā)生之前，會(huì)有一段漸變的發(fā)展過程，在這個(gè)過程中會(huì)出現(xiàn)一些故障征兆信號(hào)。如果能提前識(shí)別出這些故障征兆，盡早采取措施，就可以避免演變?yōu)閲?yán)重故障。

振動(dòng)分析是一種有效的狀態(tài)檢測(cè)方法，尤其對(duì)于旋轉(zhuǎn)機(jī)械設(shè)備[4]。如廣泛應(yīng)用于齒輪箱、軸承等風(fēng)機(jī)部件振動(dòng)故障檢測(cè)的頻譜分析[5]、小波變換[6]和Hilbert-Huang 變換[7-8]等時(shí)頻分析技術(shù)，其中 Hilbert-Huang變換在處理旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)上要比前2種方法更為有效，但其缺點(diǎn)是耗時(shí)長(zhǎng)[9]。另外，還有一些基于模型的故障診斷算法，如隱Markov模型[10-11]、支持向量機(jī)[12]和神經(jīng)網(wǎng)絡(luò)算法[13]等，它們一般都需要從振動(dòng)信號(hào)中提取特征，并對(duì)數(shù)學(xué)模型進(jìn)行訓(xùn)練，以實(shí)現(xiàn)對(duì)振動(dòng)趨勢(shì)的預(yù)測(cè)和故障識(shí)別。在過去的二十多年里，子空間識(shí)別方法應(yīng)用于振動(dòng)信號(hào)故障分析的發(fā)展也很快[14-15]，特別是對(duì)建筑物[16]和旋轉(zhuǎn)設(shè)備[17]等都有較多的應(yīng)用，但目前用于風(fēng)電機(jī)組振動(dòng)信號(hào)故障預(yù)測(cè)的研究還較少。該方法的特點(diǎn)是，直接在時(shí)域里分析數(shù)據(jù)建立模型，并識(shí)別相應(yīng)的參數(shù)，不僅具有良好的數(shù)值穩(wěn)定性和簡(jiǎn)易性，而且狀態(tài)空間方程的形式非常適合于預(yù)測(cè)、濾波、估計(jì)和控制[18]。

1 問題的描述

子空間方法是一種利用已知數(shù)據(jù)建立多變量的線性狀態(tài)空間模型的黑盒子識(shí)別方法，非常適合于振動(dòng)信號(hào)的建模分析。隨機(jī)子空間的線性狀態(tài)空間模型描述如下：

其中，Xk∈Rn和Yk∈Rl分別為在離散時(shí)刻k的狀態(tài)量和系統(tǒng)的輸出量；wk∈Rn和vk∈Rl分別為系統(tǒng)噪聲和測(cè)量噪聲，一般為不可觀的向量；A∈Rn×n為系統(tǒng)矩陣，描述系統(tǒng)的動(dòng)態(tài)行為；C∈Rl×n為系統(tǒng)輸出矩陣。

本文利用子空間方法進(jìn)行基于振動(dòng)信號(hào)的風(fēng)電機(jī)組齒輪箱故障預(yù)測(cè)。其基本思路為：首先，建立形如式（1）所示的齒輪箱隨機(jī)狀態(tài)空間模型；然后，利用正常運(yùn)行振動(dòng)監(jiān)測(cè)數(shù)據(jù)估計(jì)線性模型的參數(shù)矩陣A和C，并計(jì)算出系統(tǒng)穩(wěn)態(tài)運(yùn)行時(shí)矩陣A的特征值，將其作為齒輪箱線性動(dòng)態(tài)系統(tǒng)的參考特征值。當(dāng)齒輪箱穩(wěn)態(tài)運(yùn)行時(shí)，由實(shí)時(shí)振動(dòng)數(shù)據(jù)計(jì)算出的特征參數(shù)，與系統(tǒng)矩陣A的參考特征值誤差很?。欢?dāng)齒輪箱處于異常狀態(tài)時(shí)，求得的系統(tǒng)矩陣A的特征值會(huì)偏離參考特征值，這樣就可以實(shí)現(xiàn)對(duì)齒輪箱異常狀態(tài)的識(shí)別。當(dāng)特征值較多時(shí)，為避免對(duì)每個(gè)特征值都進(jìn)行比較，定義了均方根誤差這個(gè)總體評(píng)價(jià)指標(biāo)，通過該指標(biāo)可以從數(shù)值上直觀識(shí)別出齒輪箱的故障狀態(tài)。

2 隨機(jī)子空間識(shí)別算法

2.1 正交投影

定義由量測(cè)量組成的分塊Hankel矩陣Y：

其中，｛yk｝k=1，…，i+j+N為系統(tǒng)的 i+j+N 個(gè)測(cè)量輸出；Yp和Yf均為分塊Hankel矩陣，其下標(biāo)分別表示歷史測(cè)量值和未來測(cè)量值；i和j+1分別為歷史測(cè)量和未來測(cè)量形成分塊矩陣的行數(shù)，一般不小于系統(tǒng)模型最大階數(shù) n，min｛i，j+1｝≥n；矩陣 Yp的維數(shù)為 i×N，矩陣Yf的維數(shù)為（j+1）×N，N?max｛i，j+1｝。

對(duì)Y進(jìn)行重新分塊，如式（3）所示。

其中，Y+p為Yf的第1行移到了Yp的末行之后的矩陣；Y-f為Yf去掉了第1行后的矩陣。

將Yf正交投影到Y(jié)p的空間上，則正交投影Pm可以定義為：

其中，表示求Moore-Penrose逆。

同樣對(duì)于式（3），將 Y-f正交投影到Y(jié)+p的空間上，則有正交投影Pm-1為：

2.2 奇異值分解

對(duì)投影矩陣Pm進(jìn)行奇異值分解，則有：

其中，U1=[u1，…，un]和 V1= [v1，…，vn]為酋矩陣，分別包含了 n 個(gè)左右奇異值向量；S1=diag｛σ1，…，σn｝為對(duì)角矩陣，σi為 S1的第 i個(gè)奇異值；U0、V0和 S0為零矩陣。

由文獻(xiàn)[19]可知，投影矩陣Pm也可以寫成可觀矩陣Γm和卡爾曼濾波狀態(tài)序列的乘積，即：

同樣，投影矩陣Pm-1也可以寫成如下形式：

其中，Γm-1為Γm去掉最后一行CAm-1后的矩陣。

由式（6）、（7），則可得到可觀矩陣 Γm和狀態(tài)量X^m的表達(dá)式：

2.3 估計(jì)系統(tǒng)矩陣A和輸出矩陣C

其中，ρw和 ρv為殘差向量，與不相關(guān)。

利用最小二乘法，定義目標(biāo)函數(shù)為殘差，當(dāng)目標(biāo)函數(shù)最小時(shí)，則可以估計(jì)出系統(tǒng)矩陣A和輸出矩陣C為：

至此，就由已知的觀測(cè)量Y，得到了隨機(jī)線性狀態(tài)空間方程（1）。

3 基于振動(dòng)數(shù)據(jù)的風(fēng)機(jī)齒輪箱故障預(yù)測(cè)

3.1 齒輪箱線性狀態(tài)空間模型的確立

對(duì)于風(fēng)電機(jī)組齒輪箱，由式（4）—（13）可以獲得其隨機(jī)子空間模型。在該過程中，由于利用投影矩陣的定義式（4）和（5）求其值，涉及大維數(shù)的矩陣運(yùn)算，計(jì)算量很大。為了避免由定義直接計(jì)算投影矩陣Pm，下面給出了簡(jiǎn)化計(jì)算方法。

首先，對(duì)量測(cè)量 Hankel矩陣式（2）和（3）分別進(jìn)行LQ分解，得：

然后，計(jì)算出 Yp、Yf、Y+p和 Y-f，分別代入投影矩陣的定義式（4）和（5），進(jìn)行推導(dǎo)得到：

利用式（16）、（17）取代式（4）、（5）計(jì)算投影矩陣，也可以避免矩陣的求逆運(yùn)算。

由于Q1是行正交矩陣，所以矩陣L21的奇異值分解與Pm的奇異值分解具有相同的可觀矩陣Γm，證明如下。

投影矩陣為：

其中，Q1為正交陣。

投影矩陣又等于可觀矩陣和狀態(tài)序列的乘積：

若直接對(duì)投影矩陣進(jìn)行奇異值分解，則：

其中，U、V為正交陣；S為對(duì)角陣。

聯(lián)立式（18）—（20）得式（21）：

令可觀矩陣為：

對(duì)式（21）等號(hào)兩邊右乘 QT1，則式（21）轉(zhuǎn)化為：

令 U1=U、VT1=VTQT1，則式（23）變?yōu)椋?/p>

顯然U1是正交陣，由于：

則V1是正交陣；式（24）是關(guān)于L21的奇異值分解。令可觀矩陣為：則式（22）與式（27）相等，即可觀矩陣是相等的。

對(duì)L21進(jìn)行奇異值分解，則有：

因此，利用式（28）可以得到與式（9）相同形式的可觀矩陣Γm，而對(duì)于狀態(tài)量，則有：

將可觀矩陣Γm上移1行即可獲得可觀矩陣Γm-1，再利用投影矩陣式（17），可以計(jì)算出式（11）的狀態(tài)量m+1。因此，根據(jù)2.3節(jié)可以得到齒輪箱隨機(jī)狀態(tài)空間模型為：

其中，Yi為由振動(dòng)傳感器采集到的振動(dòng)數(shù)據(jù)，傳感器配置個(gè)數(shù)決定了Yi的維數(shù)；為由2.2節(jié)計(jì)算出的狀態(tài)量；與分別為由式（13）計(jì)算出的齒輪箱振動(dòng)狀態(tài)空間方程的系統(tǒng)矩陣和輸出矩陣，是一個(gè)方陣，其階數(shù)由奇異值矩陣S1確定，的列數(shù)與相同，行數(shù)等于振動(dòng)傳感器的個(gè)數(shù)；wi和vi為噪聲。

3.2 齒輪箱故障預(yù)測(cè)

為了有效利用隨機(jī)子空間模型，需要確定其計(jì)算方法的穩(wěn)定性。判斷原則為，即所求系統(tǒng)矩陣的所有特征值均在單位圓以內(nèi)，則說明該隨機(jī)模型是穩(wěn)定的，反之則不穩(wěn)定。

為了便于比較，引入均方根誤差（RMSE）作為評(píng)價(jià)齒輪箱狀態(tài)的指標(biāo)。由于特征值可能為復(fù)數(shù)，因此定義該指標(biāo)為：

其中，n為系統(tǒng)階數(shù)；λi為由實(shí)時(shí)振動(dòng)數(shù)據(jù)求得的特征值。

3.3 預(yù)警值與門檻值的確定

通過定義預(yù)警值與門檻值這2個(gè)閾值，利用RMSE進(jìn)行判斷，來實(shí)現(xiàn)齒輪箱的故障告警和識(shí)別。

通過對(duì)風(fēng)電機(jī)組齒輪箱振動(dòng)監(jiān)測(cè)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)特性分析可知，實(shí)時(shí)求得的特征值λi散落在參考特征值的附近，誤差較??；并對(duì)每個(gè)實(shí)特征值、每個(gè)復(fù)特征值的實(shí)部和虛部與相應(yīng)的參考特征值的殘差進(jìn)行分析，由分析可知，每個(gè)特征值的殘差皆為正態(tài)分布：

其中，μi為齒輪箱振動(dòng)信號(hào)子空間模型特征值殘差的均值；σi為其標(biāo)準(zhǔn)差。因此，可利用統(tǒng)計(jì)過程控制（SPC）原理[20]定義齒輪箱故障判別的 2個(gè)閾值。

根據(jù)SPC原理，服從正態(tài)分布的特征值的殘差上下限為 si=μi±m(xù)σi（m=2，3）。 利用該限值定義齒輪箱故障判別指標(biāo)RMSE的預(yù)警值和門檻值。m=2時(shí)，R*定義為預(yù)警值，當(dāng)RMSE超過該值，意味著將要出現(xiàn)故障；而m=3時(shí)，R*為故障門檻值，當(dāng)RMSE超過該值時(shí)，則認(rèn)為已經(jīng)故障。

因此，根據(jù)實(shí)時(shí)數(shù)據(jù)不斷計(jì)算系統(tǒng)矩陣的特征值和相應(yīng)的RMSE值，并與所設(shè)閾值進(jìn)行比較，進(jìn)而可以評(píng)估出齒輪箱的運(yùn)行狀態(tài)。

4 算例分析

為了充分驗(yàn)證該方法的有效性，對(duì)一個(gè)實(shí)際發(fā)生齒輪箱斷齒故障的風(fēng)電機(jī)組的振動(dòng)監(jiān)測(cè)數(shù)據(jù)進(jìn)行研究分析。首先，利用該風(fēng)機(jī)齒輪箱故障之前近半年的振動(dòng)監(jiān)測(cè)數(shù)據(jù)建立其隨機(jī)子空間模型，分段計(jì)算系統(tǒng)矩陣 A^k；然后，利用所提方法識(shí)別系統(tǒng)參考特征值λi；定義齒輪箱故障的判別函數(shù)，并確定出相應(yīng)的閾值，進(jìn)而就可以實(shí)現(xiàn)對(duì)齒輪箱的故障預(yù)測(cè)，具體計(jì)算過程如圖1所示。

圖1 基于子空間方法的風(fēng)機(jī)齒輪箱故障預(yù)測(cè)算法Fig.1 Gearbox fault prediction algorithm based on subspace method for wind turbine

風(fēng)電機(jī)組容量為1.5 MW，齒輪箱傳動(dòng)結(jié)構(gòu)為2級(jí)行星齒輪+1級(jí)平行軸，齒輪箱傳動(dòng)比率為100.48；所安裝的傳感器為加速度傳感器，采樣率為12 kHz，所采集到的振動(dòng)信號(hào)為相應(yīng)的加速度信號(hào)，單位為m／s2，所分析的正常和異常振動(dòng)數(shù)據(jù)如圖2所示。

算法中每一段數(shù)據(jù)窗口的選擇應(yīng)盡可能足夠長(zhǎng)，以滿足子空間方法的參數(shù)識(shí)別有效性。對(duì)該風(fēng)機(jī)齒輪箱的振動(dòng)監(jiān)測(cè)數(shù)據(jù)，窗口長(zhǎng)度選為6000，其中Hankel矩陣 Yp為 7×5994維矩陣，Yf為 7×5994維矩陣。對(duì)形成的測(cè)量矩陣進(jìn)行LQ分解，得到子塊矩陣L21，然后利用式（28），對(duì)矩陣 L21進(jìn)行奇異值分解，從而可以確定齒輪箱的隨機(jī)子空間模型的階數(shù)n為 3，其對(duì)應(yīng)的參考特征值i（i=1，2，3）為：

圖2 正常與異常狀態(tài)下的振動(dòng)信號(hào)Fig.2 Normal and abnormal vibration signals

圖3給出了參考特征值的分布。從圖中可以看出3個(gè)參考特征值都分布在單位圓內(nèi)，說明識(shí)別出的齒輪箱隨機(jī)子空間模型是穩(wěn)定的。利用穩(wěn)態(tài)運(yùn)行時(shí)的振動(dòng)采樣數(shù)據(jù)，代入齒輪箱故障判別模型，實(shí)時(shí)計(jì)算出的特征值與參考特征值之間的誤差非常小，說明該風(fēng)機(jī)齒輪箱的運(yùn)行狀況很好，是健康的。

圖3 齒輪箱正常狀態(tài)下的參考特征值分布Fig.3 Distribution of reference eigenvalues of gearbox operating in normal state

分析圖2中的異常振動(dòng)信號(hào)，其由正常數(shù)據(jù)和故障數(shù)據(jù)構(gòu)成，并從故障前一段時(shí)間的振動(dòng)監(jiān)測(cè)數(shù)據(jù)開始分析，進(jìn)行齒輪箱子空間模型的參數(shù)計(jì)算。通過故障判別指標(biāo)的計(jì)算結(jié)果發(fā)現(xiàn)，其中1個(gè)實(shí)數(shù)特征值變化不大，但另外2個(gè)復(fù)特征值逐漸偏離其各自的參考特征值，而且變化趨勢(shì)逐漸增大，但仍然位于單位圓內(nèi)。圖4給出了各個(gè)特征值的演變過程。

為了更加清楚地說明齒輪箱狀態(tài)的變化過程，利用定義的綜合指標(biāo)分析齒輪箱故障發(fā)生的過程。通過對(duì)長(zhǎng)時(shí)間穩(wěn)態(tài)運(yùn)行的齒輪箱的振動(dòng)數(shù)據(jù)的分析，可以計(jì)算出各個(gè)穩(wěn)態(tài)特征值與相應(yīng)參考特征值殘差的均值μi和標(biāo)準(zhǔn)差σi，并利用SPC原理確定其殘差的上下限為 si=μi±m(xù)σi，具體計(jì)算結(jié)果見表1。 再將si代入式（33）中，則可以得到RMSE指標(biāo)的預(yù)警值為0.0084、門檻值為0.0126。

圖4 齒輪箱異常狀態(tài)的特征值分布圖Fig.4 Distribution of eigenvalues of gearbox operating in abnormal state

表1 特征值殘差上下限Table 1 Upper and lower limits of residual errors of eigenvalues

圖5給出了齒輪箱穩(wěn)態(tài)運(yùn)行時(shí)的RMSE指標(biāo)變化過程（圖中每次仿真對(duì)應(yīng)波形上的1個(gè)點(diǎn)，圖6中同），RMSE值均在預(yù)警值以下范圍內(nèi)變化，說明了齒輪箱處于正常運(yùn)行狀態(tài)。

圖5 齒輪箱正常狀態(tài)下的均方根誤差變化曲線Fig.5 Curve of RMSE of gearbox operating in normal state

圖6 齒輪箱故障過程中的均方根誤差變化曲線Fig.6 Curve of RMSE of gearbox operating in abnormal state

圖6為齒輪箱故障過程中的RMSE。從圖6可以看出，對(duì)應(yīng)圖4中的特征值與參考特征值，RMSE出現(xiàn)了越過預(yù)警值的情況，甚至最后越過了門檻值。在剛開始的一段RMSE值變化也不大，均在預(yù)警值以下，但隨著特征值點(diǎn)的發(fā)散，RMSE值越來越大；大約在第609個(gè)點(diǎn)，曲線越過了預(yù)警值0.0084，說明齒輪箱存在異常，有故障的趨勢(shì)；在第669個(gè)點(diǎn)以后，曲線增長(zhǎng)到0.0126附近，并且在門檻值附近波動(dòng)；從第788個(gè)點(diǎn)以后迅速上升，且完全越過了門檻值，并且在越限以后曲線繼續(xù)上升，此時(shí)可以認(rèn)為齒輪箱已經(jīng)故障。故應(yīng)該在第609個(gè)點(diǎn)到第669個(gè)點(diǎn)之間發(fā)出預(yù)警信號(hào)，提醒運(yùn)行人員采取應(yīng)對(duì)措施。

5 結(jié)論與展望

本文提出了一種子空間識(shí)別算法，用于僅有觀測(cè)量已知的風(fēng)機(jī)齒輪故障預(yù)測(cè)。該算法利用了幾何投影的概念，通過投影的行空間可以確定系統(tǒng)的狀態(tài)量i；通過投影的列空間可以確定系統(tǒng)的可觀矩陣Γm；通過對(duì)投影矩陣的奇異值分解可以確定系統(tǒng)的階數(shù)n，進(jìn)而可以確定系統(tǒng)矩陣的維數(shù)，并計(jì)算相應(yīng)的特征值和RMSE指標(biāo)，以實(shí)現(xiàn)齒輪箱運(yùn)行狀態(tài)的評(píng)估。

通過對(duì)振動(dòng)數(shù)據(jù)的仿真分析，可以得出以下結(jié)論：首先，由該隨機(jī)子空間算法識(shí)別出的隨機(jī)狀態(tài)空間模型具有良好的穩(wěn)定性；其次，當(dāng)已知參考特征值后，可以用其作為評(píng)判齒輪箱狀態(tài)的一個(gè)直觀參考；最后，引入了RMSE指標(biāo)對(duì)各個(gè)特征值的分布情況進(jìn)行定量評(píng)價(jià)，當(dāng)齒輪箱處于不同運(yùn)行狀態(tài)時(shí)，通過該故障預(yù)測(cè)算法給出的相應(yīng)RMSE變化趨勢(shì)，可以對(duì)齒輪箱進(jìn)行故障預(yù)測(cè)，給現(xiàn)場(chǎng)監(jiān)測(cè)和制定檢修計(jì)劃提供一定的指導(dǎo)，在提高風(fēng)電機(jī)組的經(jīng)濟(jì)性和可靠性方面具有重要意義。

本文所提方法雖然能夠預(yù)測(cè)出故障的發(fā)生，但在識(shí)別出具體故障方面還需要完善。將繼續(xù)從兩方面對(duì)該算法進(jìn)行深化研究：一是研究子空間算法中特征值與故障類型之間的關(guān)系；二是研究子空間與頻譜分析相結(jié)合的故障預(yù)測(cè)方法，以實(shí)現(xiàn)風(fēng)電機(jī)組齒輪箱的故障識(shí)別。