曹超凡, 李明亮, 蔣雙云, 張廣濤, 李中梁, 盧 娜
(1. 鄭州大學 水利與交通學院,鄭州 450001; 2. 華電電力科學研究院有限公司,杭州 310030;3. 潤電能源科學技術有限公司,鄭州 450052)
水電機組實際運行過程中各種工況轉換頻繁,加之受到水力、電氣、機械等多方面綜合因素影響,可能產生多種故障,因此需開展機組狀態趨勢預測,以掌握機組實時運行狀態,及時發現機組運行異常和早期故障,并據此制定合理的檢修計劃,避免發生重大事故[1-2]。
水電機組狀態趨勢預測是利用機組歷史監測數據進行學習構建趨勢預測模型,并基于當前機組運行數據實現未來機組運行狀態的預測[3-4]。在實際運行中,水電機組故障通常在振動信號上有所體現,提取水電機組振動信號的重要特征進行機組狀態趨勢預測,可以有效反應機組真實運行狀態,并實現故障預警[5-6]。
目前專家學者針對旋轉機械開展故障預警研究已經取得一定進展。鹿衛國等[7]假定水電機組正常狀態下歷史振動數據服從某一概率分布,將偏離概率分布的數據歸類為故障狀態實現故障預警。安學利等[8]利用最小二乘支持向量機構建三維曲面預警模型,通過對機組實時運行振動數據進行狀態評估實現故障預警。劉濤等[9]利用多元狀態評估技術對電廠風機觀測向量進行最優估計,從而得到估計向量,通過計算觀測向量與估計向量之間的偏離度,判斷是否超過預警閾值,從而實現風機的實時故障預警。然而,現有的預警方法多存在預警指標特征信息單一,所構建預警指標難以完全表征機組狀態,早期預警困難的問題。而利用振動信號中的多元特征構建綜合預警指標,可以更加有效地判斷機組健康狀態。為此,劉東等[10]通過計算水電機組振動信號時域特征與健康值之間的相對差值作為時域劣化指標,計算實測信號頻域特征向量與健康聚類中心之間的歐氏距離作為頻域劣化指標,實現了水電機組劣化在線評估。Mao等[11]提出一種具有交替最小化方案的算法進行張量表示和無監督特征提取,并利用提取的特征構建健康指標來進行地鐵車輪劣化狀態評估。Liu等[12]利用深度卷積神經網絡提取軸承退化數據與正常數據之間的特征距離構建演化狀態指標,然后利用無監督聚類方法劃分狀態階段,實現軸承早期故障演化監測。盡管如此,利用綜合指標實現機組早期故障預警的研究仍然較少,且當前研究中大多使用單傳感器振動數據作為模型輸入,其所提取的特征信息遠不如多源數據的全面,并且容易遺漏重要特征。
針對上述問題,并考慮到遺傳規劃(genetic program,GP)可以更有效地構造表征機組狀態的高級特征,進而獲得更為可靠的預警結果,且目前其在故障預警領域應用較少,本文提出了基于集成多傳感器遺傳規劃(integrated multi-sensor genetic programming, IMSGP)與權重歐式距離指標(weighted euclidean distance index, WEDI)的水電機組故障預警方法。本文主要工作為:
(1) 融合多傳感器信息,引入GP算法,構建了IMSGP故障預警模型,可更加全面地反映機組運行狀態。
(2) 綜合考慮時頻特征信息構建多元原始預警特征集,利用復合檢測指數(composite detection index,CDI)進行特征選擇,剔除無關特征向量,提高了特征敏感性,降低了后續算法的復雜性。
(3) 結合主成分分析(principal component analysis,PCA)[13]與歐式距離[14]構建WEDI綜合預警指標,以偏離健康狀態權重歐式距離系數3倍偏差為預警閾值,實現水電機組早期故障預警。
(4) 利用所提出的方法對水電機組實測信號進行處理,結果表明,所構建的狀態趨勢預測模型能夠及時發現水電機組狀態變化,預警時間比明顯出現癥狀的時間提前6天左右。
遺傳規劃是遺傳算法與計算機編程相結合的一種進化方法[15],具有利用候選解表示的動態結構和全局搜索的特點,但與其他進化算法不同的是,GP采用更靈活的層次模型(樹)來表示解空間,其結構和大小可自適應調整,更適合于表示復雜問題[16]。GP中的個體通常由函數集和終端集組成。GP的函數集可分為一般函數(如+,-,×,÷等)、邏輯函數(如if,and,or等)和其他數學函數(如exp,sin,cos等)。終止符集通常為變量和特征等。如圖1所示,該樹模型是GP的一個個體,數學上可以表示為[8×(x+5)]-(3÷y),其中+、-、×和÷包含在函數集中,終端集包含x、y、8、3和5。在水電機組故障特征構造中,該方程可作為所構造的高級特征,其中,x和y為水電機組信號的兩個原始特征。由此可以看出,GP算法具有很好的可視化和模型解釋能力。此外,GP算法樹深度可以根據實際問題進行調整,模型復雜程度與問題復雜程度相關。

圖1 某個體樹形GP結構表達Fig.1 Tree structure representation of some individual in GP
GP算法的流程圖如圖2所示,主要包含:創建初始個體、個體適應度評價、選擇、交叉、變異、判斷是否滿足終止條件,結束。其中選擇指GP算法中個體的選擇與復制操作,其遵循生物界的“優勝劣汰”原則,目的是保留種群中優秀的個體,淘汰劣質個體,根據適應度函數值,選擇復制上一代優秀個體,保證下一代個體更優秀;交叉指選擇父代中的兩個個體,相互替代他們的部分結構,從而產生兩個新的個體,該操作是為了使不同個體的信息相互交換,豐富種群的多樣性,根據適應度值與被選定概率成正比的關系,隨機選擇父代個體,也可以通過隨機選擇父代樹節點作為交叉節點,以交叉點為根的整個子樹為交叉段,兩個子樹相互交換,產生新的樹結構;變異指GP個體的變異,包含對終止符集和函數符集進行變異操作,變異位置是隨機選擇父代樹的任一節點,對于原位置任意替換,當變異點為函數符集,則替換前后運算數目元素相同。

圖2 遺傳規劃算法流程圖Fig.2 Genetic programming algorithm
對于預處理后的多元水電機組信號,需提取其基礎特征信息作為后續IMSGP模型的輸入。本文選取了15個時域特征和13個頻域特征進行基礎特征信息提取。
時域特征分別為算術平均值、平均幅值、方差、標準差、方根幅值、均方根值、最大值、最小值、峰值、峰值因子值、峰度值、波形因子值、脈沖因子值、邊際因子值和偏度值。時域特征計算式如表1所示,其中,v(i)為信號樣本點,N為樣本點數。

表1 時域特征參數Tab.1 Time domain feature parameters
頻域特征考慮了信號的頻率特性,通過對信號進行快速傅里葉變換(fast Fourier transform, FFT)[17]得到的頻譜幅值進行計算,提取出信號的頻率特性,13個頻域特征為:平均頻率F1、標準頻率的偏差F2、頻率中心F3、均方根頻率F4、主波段的位置變化特征F5和F6、頻譜的集中或分散程度特征F7~F13。計算公式如表2所示,其中,s(i)為原始信號的頻譜,M為譜線數,f(j)為第j條譜線的頻率值,Aj為第j條譜線的幅值。

表2 頻域特征參數Tab.2 Frequency domain feature parameters
從多傳感器信號中能夠獲得大量的特征信息,但其中不乏冗余的無用特征,為了提高計算信號特征的性能,降低后續GP的計算復雜度,本文利用復合檢測指數CDI作為特征選擇的敏感特征評價指標,將進行多元原始信號特征集特征篩選。CDI的計算步驟為:
步驟1計算檢測指數DI值
假設Ni,Nj分別是由旋轉機械狀態i和狀態j的信號計算得到的時頻域特征值,且服從正態分布。則狀態i和狀態j下某一特征值的DI[18]值為
(1)
式中:σi,σj為Ni,Nj的標準差值;μi,μj為Ni,Nj的平均值。
步驟2計算CDI值
由式(1)可知,故障樣本的類間間距越大,類內間距越小時,對于故障分類越有利。因此,DI值越大,越容易識別出對應的兩種不同狀態樣本的差別,據此,本文構造CDI為
(2)
式中:n為故障類型的種類數; DIi,j為狀態i和狀態j對應特征的DI值。
可見,CDI綜合考慮到不同故障類型,同一特征值計算得到的DI平均值與最小值。CDI越大,表明該特征對所有故障類型的敏感性越高。
考慮到將多傳感器數據輸入GP后,在各設定優化參數不變的前提下重復運行GP模型,其輸出的結果會有差別,利用其進行后續的預警指標構建得到的預警結果有優有劣,即多傳感器遺傳規劃(multi-sensor genetic programming, MSGP)運行結果不穩定,魯棒性不強,因此本文引入集成化思想,構建了集成多傳感器遺傳規劃IMSGP模型,這里的集成是指通過重復運行MSGP得到高維高級特征集合,即假設集成k次,就會得到k個互相獨立的一維高級特征,也就是k維高級特征。這樣,所得的k維高級特征就可以輸入到PCA算法進行降維處理,獲得低維敏感特征,提高模型的泛化性能及水電機組故障預警的效率。
GP中目標函數值反映了種群個體的特性,在優化程序運行時,常常通過目標函數評價來尋找最佳個體。而當訓練樣本數量較少時,GP構造的特征在訓練集上容易獲得較好的分類性能,但在測試集上可能會產生過擬合的現象。為了解決這一問題,本文構建了一種融合分類精度和樣本距離的遺傳規劃適應度函數,其計算公式如式(3)所示,利用最小-最大歸一化方法將Fitness變換到 [0,1]范圍內。
(3)
式中:Acc為構造特征的分類準確率,通過k-fold交叉驗證計算得出,在本文計算中取5折[19]; Dist為距離度量,其原理是使樣本的類內距離最小,類間距離最大,以此來降低同類樣本的差異性同時提高構造特征的鑒別能力。Dist計算公式如式(4)~式(5)所示,其計算基礎是兩個樣本特征的JS(jensen-shannon)散度值[20],同類樣本之間的距離計算公式如式(4)所示,類間樣本之間的距離計算公式如式(5)所示,最后受sigmod函數的啟發,利用式(6)對距離進行歸一化處理,保證式(3)處于[0,1]范圍內。
(4)
(5)
(6)
式中:Din為類內距離;Dout為類間距離;Ni,Nj為水電機組運行狀態i,j下的特征個數;X為IMSGP構造的高級特征。
為及時發現水電機組早期故障,本文結合PCA算法與歐式距離構建一種權重歐式距離故障預警指標,其流程如圖3所示,具體步驟為:

圖3 權重歐式距離計算流程Fig.3 Weighted Euclidean distance calculation procedure
步驟1設定PCA算法貢獻率閾值,將貢獻率達到閾值的前z個指標提取出來,從而將高維故障特征集合降到z維度,并記錄降維特征的貢獻率。
步驟2計算健康狀態下歷史信號降維特征的聚類中心,并將其作為健康基準值。
(7)
式中:z為降維故障特征集合的維度;n為健康狀態下歷史信號樣本;X為降維后特征值;H為健康基準值。
步驟3計算水電機組待測信號特征與健康基準值之間的歐式距離。
步驟4將步驟3得到的歐氏距離與步驟1中相應的貢獻率相乘并累加,得到最終的故障預警指標。
(8)
式中:N為待測信號樣本數;w為z維故障特征對應的貢獻率,WED即為所求權重歐式距離。
本文所提出的水電機組故障預警方法流程如圖4所示。具體步驟為:

圖4 基于IMSGP-WEDI的水電機組故障預警方法流程Fig.4 Fault early warning process of hydropower unit based on IMSGP-WEDI
步驟1獲取水電機組歷史正常狀態信號和待測信號多傳感器數據,利用小波閾值降噪[21]方法對振動信號進行預處理。
步驟2提取預處理后信號的時頻域特征,并進行CDI特征選擇。
步驟3已經CDI選擇后的特征參數為輸入變量,以式(3)為目標函數,進行集成多傳感器遺傳規劃,得到由多個高級特征表達的高維特征向量。
步驟4設定PCA降維閾值,利用PCA將高維特征向量映射到低維空間,得到代表水電機組運行狀態的融合特征向量。
步驟5以正常運行狀態下信號數據融合特征向量的均值作為聚類中心,并計算正常狀態數據融合特征與聚類中心的權重歐式距離標準差。
步驟6以三倍標準差為預警閾值,將第一次超過報警閾值所對應的數據定義為預警樣本,并在此時觸發報警。
水電機組的結構如圖5所示,主要包括一個發電機、一個水輪機、上下機架、三個導軸承和一個推力軸承,發電機與水輪機的軸通過法蘭連接。

圖5 水電機組結構圖Fig.5 Hydropower unit structure
在運行過程中,工作人員發現電站3號機組出現嚴重振動,經專業人員排查,得知該故障是由水力不平衡因素引起。數據來源于安裝在上導軸承、推力軸承和水導軸承上各個方向的9個振擺數據和7個工況數據(包含水頭、導葉開度等)共16個傳感器數據。機組停機檢修前,先后經歷了正常、故障預警、故障三種狀態。
本文選取三種不同狀態下水電機組的信號樣本進行訓練,采樣頻率為458 Hz,采樣點數為2 048,每類數據包含20個樣本,其中訓練樣本每類狀態下取5個,共15個,測試樣本每類狀態下取15個,共45個,預處理后的信號樣本經過特征選擇后輸入IMSGP模型中運行20次。振動信號樣本時序波形(以水導軸承x向擺度傳感器、上導軸承x向擺度傳感器和上導軸承y向擺度傳感器為例)如圖6所示。

圖6 水電機組部分傳感器時域波形圖Fig.6 Partial time series waveforms of hydropower generating unit
根據機組狀態將信號數據劃分為歷史數據和待測數據。
(1) 歷史信號數據
采集水電機組健康狀態下多傳感器信號作為歷史數據計算預警閾值,共有253個數據樣本,每個樣本包含4 096個數據點。
(2) 待測信號數據
采集故障發生前后機組多傳感器信號作為待測數據,采集時間為2015-08-15—2015-08-28,共有375個數據樣本,每個樣本包含4 096個數據點。
采用小波閾值降噪方法對振動信號數據進行預處理。小波閾值降噪方法的參數設置為:小波分解層數設定為3層,閾值選擇方法為“sqtwolog”方法,閾值處理方法為軟閾值方法,閾值縮放方法為“mln”方法,小波基函數為“DB4”小波。
IMSGP參數設置如表3所示,計算16個傳感器數據的時頻域特征,得到448維特征集合,各個特征對應的CDI值如表4所示,選擇每個傳感器中的CDI值較大的前6個特征,得到96維特征,其對應的CDI值如表5所示。

表3 IMSGP參數設置Tab.3 IMSGP parameter setting

表4 高維特征集的復合檢測指數值Tab.4 The CDI of high-dimensional feature set

表5 特征選擇后對應特征的CDI值
將所得96維時頻特征作為遺傳規劃的函數集,訓練IMSGP模型,其所得高維特征集合如表6所示。

表6 IMSGP所得歷史數據的高級特征向量集
將待測信號樣本輸入到訓練完成的IMSGP模型中,得到待測數據樣本的高維特征集合,如表7所示。

表7 IMSGP 所得待測數據的高級特征向量集
設定PCA降維特征貢獻率為0.95,利用PCA算法對待測信號樣本的高級特征進行降維處理,降維特征貢獻率如圖7所示。由圖7可知,前兩個特征的貢獻率依次為0.747 0和0.235 8,累計貢獻率為0.982 8,超過所設定閾值0.95,因此,將原來的20維度特征向量降低到2維,其特征值如表8所示,記錄這兩個特征的貢獻率值作為后續權重歐式距離的計算權重。根據權重歐氏距離計算方法,選擇表8中前50組正常狀態下融合特征的平均值作為健康狀態基準值,計算待測樣本與健康基準值的權重歐氏距離并作為故障預警指標值。

表8 經PCA降維后的高級特征值Tab.8 High-level feature values after PCA

圖7 高級特征對PCA算法的累計貢獻率Fig.7 Cumulative contribution rate of high-level features to PCA
經過IMSGP特征提取與PCA降維后,將得到的權重歐氏距離數據減去其均值進行歸一化處理,使預警指標從0開始,計算健康狀態下253組數據的融合特征與健康基準值之間權重歐氏距離的標準差,通過計算,得到正常狀態下253組數據的融合特征與健康基準值之間的標準差為0.315 8,三倍標準差為0.947 3,即為故障預警指標的閾值[22]。
將待測樣本對應的故障預警指標值與預警閾值繪于同一張圖上,如圖8所示,由圖8可知,在08-15—08-20,機組預警指標一直穩定在較低水平,而且波動微乎其微;在08-22左右機組振動預警指標首次突破閾值,可認為在該時段內機組開始產生故障,將發出預警。通過與運行管理人員溝通交流得知,在08-15—08-20,機組振動及各項參數值均處于正常水平,無明顯波動;根據電站的事后分析報告,發現機組在08-28開機過程中,上機架和水輪機蝸殼及尾水管等處異常聲音明顯,通過檢修發現是由水力不平衡故障引起的。可以看出,利用所提出的預警方法得到的結果表明機組早在08-22前后就出現了故障趨勢,因此本文提出的故障預警指標能夠有效反映水電機組的運行狀態并實現機組故障預警,對保證機組和電站安全具有實際應用意義。

圖8 水電機組故障預警指標計算結果Fig.8 Calculation results of fault early warning indicators
為進一步證明本文方法的有效性,本文針對圖4所示方法流程中的多傳感器數據輸入以及IMSGP步驟設置兩組對比試驗,即:
(1) 單傳感器試驗
單傳感器試驗與本文方法不同的是:在2.5節流程的步驟1中利用單傳感器數據輸入代替多傳感器輸入,相應的步驟3中以集成單傳感器GP代替IMSGP,即通過重復運行單傳感器GP模型得到高維高級特征集合。其余步驟與本文方法相同。
(2) 無集成試驗
無集成試驗與本文方法不同的是:在2.5節流程的步驟3中以MSGP代替IMSGP,另外由于只運行MSGP模型一次,得到的一維高級特征可直接進行故障預警指標計算,無需進行PCA降維,所以不進行步驟4的操作。其余步驟與本文方法相同。
如圖9所示,圖9(a)、圖9(b)、圖9(c)分別為三個單傳感器數據所得試驗結果,可以看出圖9(a)、圖9(b)兩圖結果基本無法反映機組狀態從正常到故障的劣化過程,預警閾值無法起到作用,圖9(c)所示結果基本與本文方法的結果接近,在08-23左右首次突破閾值,但其早期故障征兆不如本文故障預警指標明顯,說明多傳感器輸入可以考慮到更多有效特征;圖9(d)、圖9(e)、圖9(f)為無集成情況下試驗結果,因為GP算法本身穩定性不強,可以看出,圖9(d)、圖9(e)所對應的兩組試驗無法得到有效的結果,只有圖9(f)所示結果與本文結果相近,說明不集成時結果優劣難以預測,而通過集成可以得到更加穩定和準確的結果。

圖9 單傳感器和無集成對比試驗Fig.9 Signal sensor and no integrated comparison experiment
本文將集成多傳感器遺傳規劃方法應用到水電機組狀態趨勢預測中,結合PCA算法與歐氏距離構建了水電機組故障預警指標。通過對水電機組多傳感器數據進行試驗,分別利用歷史數據與待測數據進行模型的訓練與測試,試驗結果表明本文提出的基于IMSGP與WEDI的水電機組故障預警方法能夠有效發現早期故障,并進行及時預警。