基于協整分析的風力機多工況監測與故障診斷

汪千程 蘇 春， 文澤軍

1.東南大學機械工程學院，南京，2111892.湖南科技大學機械設備健康維護湖南省重點實驗室，湘潭，411201

0 引言

風能儲量豐富、分布廣泛、污染小、可再生，因此其開發與利用受到廣泛關注。風力機是風電系統的核心裝備，常年工作于酷暑、寒冬、風沙等極端環境，其可靠運行面臨挑戰。

有效的狀態監測與故障診斷對減少風力機的非計劃停機、避免災難性后果至關重要，現有的監測與診斷方法大致可分為數學模型方法、專家知識方法和數據驅動方法[1-2]。日趨復雜的工業裝備導致精確的數學模型和完備的專家知識往往難以獲取，裝備中普遍安裝的監控與數據采集(supervisory control and data acquisition，SCADA)系統使得數據驅動方法逐漸成為狀態監測與故障診斷領域的重要發展方向。

QIN[3]采用主成分分析(principal component analysis，PCA)模型表征裝備的正常與異常狀態。NAMDARI等[4]采用支持向量機(support vector machine，SVM)模型完成裝備早期的故障診斷。ZHAO等[5]提出一種稀疏相異度算法，以實現無先驗故障信息的裝備早期故障診斷。上述方法大多假定變量具有平穩性，難以描述非平穩變量之間的關系。

若一組非平穩的被測信號存在協整關系，利用協整分析(cointegration analysis，CA)可以消除信號的非線性趨勢和環境因素等影響，獲得協整殘差模型，其殘差序列可用于表征變量之間長期均衡的動態關系[6]。當被檢測對象的殘差偏離平穩狀態時，表明被檢測對象的狀態出現異常。

陳前等[7]將協整理論用于非平穩工程系統的狀態監測與故障診斷。ZHAO等[8]利用協整分析設置殘差閾值，完成了風力機狀態監測與齒輪箱故障預警。針對時變轉速下軸承的損傷檢測，TABRIZI等[9]提出一種基于協整的故障檢測方法。為了消除環境因素和運行不確定性的影響，ZHANG等[10]采用協整分析實現風力機的狀態監測。HU等[11]針對非平穩故障問題，提出一種雙重協整分析的故障診斷策略。現有的協整分析研究主要關注裝備正常和故障兩種工況，但在工程實際中，裝備通常存在多種運行工況，若不加以區分可能會造成故障誤報[12]，如受紊流、風向變化和相鄰風力機干擾等因素影響，風力機常處于偏航狀態，協整分析時的偏航狀態會引起殘差幅值超過正常閾值、造成故障誤報。DAO等[13]將偏航狀態定義為“異常運行狀態”，采用協整分析監測風力機狀態和齒輪箱故障，但在確定殘差預警閾值時，該研究未能有效區分風力機偏航和齒輪箱故障兩種狀態，影響故障診斷的準確性。

本文提出一種基于協整分析的裝備多工況監測與故障診斷方法，以SCADA系統的運行數據為基礎，采用隨機森林(random forest，RF)特征選擇算法提取關鍵特征變量；通過對關鍵特征序列的協整分析建立協整殘差模型，獲得最優殘差序列；基于殘差序列，采用概率圖(probability plot，PP)確定每種工況對應的殘差區間和預警閾值，實現狀態監測與故障預警。

1 問題描述與分析步驟

在系統和環境等多重隨機因素的共同作用下，裝備狀態常呈現非平穩、多工況等特征[16]。協整分析適合處理非平穩時間序列長期均衡的動態關系，可用于裝備狀態監測與故障診斷。現有的協整分析模型存在以下不足[8-11，13，15]：①關鍵特征變量的選取多依靠經驗或只針對某一類信號，適用范圍有限，可能會導致關鍵特征缺失，影響模型精度；②殘差閾值的設定主要基于統計過程控制(statistical process control，SPC)方法，難以有效劃分多重工況，易造成故障誤報，降低故障診斷的準確性。

為了解決上述問題，本文提出一種基于協整分析的裝備多工況監測與故障診斷方法，該方法的主要步驟包括[15]：

(1)構建訓練集。以SCADA系統收集的特征變量歷史數據構建訓練集，通過預處理剔除異常狀態數據，得到正常運行狀態下特征變量的時間序列{xj}，j=1，2，…，J。

(2)提取關鍵特征變量。對序列{xj}采用隨機森林特征選擇算法，提取目標變量的關鍵特征序列{xj′}，j′=1，2，…，J′；J′

(3)關鍵特征的增廣迪基-富勒(augmented Dicky-Fuller，ADF)檢驗[15]。確定關鍵特征序列的單整階數，即判斷序列是否為一階單整序列。

(4)關鍵特征的Johansen協整檢驗。確定關鍵特征的協整關系數目r并計算協整系數矩陣βr。

(5)構建測試集。以SCADA系統收集的關鍵特征運行數據構建測試集，保留異常狀態數據，得到運行狀態下關鍵特征的時間序列{zj′}。

(6)計算協整殘差模型。基于協整系數矩陣βr以及構建的測試序列{zj′}，計算協整殘差模型ξ。

(7)協整殘差模型的ADF檢驗。判斷協整殘差模型中各殘差序列的平穩性，即是否為平穩時間序列。

(8)選取最優殘差序列。以最大特征值對應的殘差模型為最佳監測模型，其殘差序列為最優殘差序列。

(9)確定各工況的殘差區間和預警閾值。利用概率圖分析最優殘差序列，劃分多工況狀態的殘差區間，并確定每種工況對應的預警閾值。

(10)狀態監測與故障診斷。觀測裝備運行中的殘差是否超過設定的閾值，完成裝備狀態監測與故障預警。

2 基于隨機森林的關鍵特征提取

隨機森林是裝袋算法的一個拓展體，采用重抽樣技術從原始數據集中隨機抽取數據構建決策樹，在生成決策樹的過程中引入隨機特征選擇策略，從海量特征中提取關鍵特征，達到降維和提升模型性能的目的[16]。

2.1 隨機森林特征重要性度量方法

隨機森林特征重要性的度量主要依據特征隨機置換前后生成的袋外數據和袋外數據誤差。袋外數據是指利用重抽樣技術建立決策樹時沒有參與決策樹構建的數據，這部分數據可以用于對決策樹的性能進行評估。袋外數據誤差為將袋外數據代入所建立的決策樹而計算得到的模型預測誤差率。特征置換后的袋外數據誤差通常會增大，而特征重要性越高，隨機森林袋外數據誤差的變化值越大[16]。因此，通過分析特征置換前后袋外數據誤差的變化，可評估每個特征對目標變量的重要性。

定義1 第k棵決策樹的袋外數據

對于訓練集D，通過bootstrap方法生成數據集Dk，建立第k棵決策樹，則未參與第k棵決策樹構建的數據D-Dk稱為第k棵決策樹的袋外數據Ok。

定義2 基于袋外數據誤差的特征j重要性度量

將袋外數據分別代入各自的決策樹，計算袋外數據誤差，然后對袋外數據中的特征j進行隨機置換，計算新的袋外數據誤差。決策樹在特征j置換前后的袋外數據誤差的平均變化量即為特征j的重要性。

2.2 隨機森林特征重要性計算步驟

假設隨機森林中的決策樹總數為K，原始數據集有J個特征。對于特征j，基于袋外數據誤差分析的特征重要性計算步驟如下：

(1)計算第k棵決策樹對應的袋外數據Ok和袋外數據誤差Ek。

(2)保持其他特征不變，對Ok中的特征j進行隨機序列置換，計算新的袋外數據Okj和袋外數據誤差Ekj。

(3)重復步驟(1)、步驟(2)，得到{Ek|k=1，2，…，K}和{Ekj|k=1，2，…，K；j=1，2，…，J}。

(4)計算決策樹在特征j置換前后袋外數據誤差的平均變化量[17]：

(1)

式中，V(j)為特征j對于目標變量的重要性。

依次計算各特征對目標變量的重要性，按由大到小排序得到關鍵特征變量。

3 基于協整分析的建模

3.1 單整與協整的定義

若一個時間序列在成為平穩序列之前需經過d(d>0)次差分，則稱該序列為d階單整，記作I(d)。如果序列{xit}(i=1,2,…,n)為d階單整，且存在一個向量η=(η1，η2，…，ηn)，使得η·Xt～I(d-b)，其中，b>0，Xt=({x1t},{x2t},…,{xnt})，則稱Xt為d-b階協整，記作Xt～CI(d-b)，η為協整系數[8]。

由協整的定義可知，如果一組變量具有保持線性關系的共同趨勢，那么協整分析可以建立這種線性關系，即協整分析可以對兩個或多個單整階數均為一階的非平穩序列進行線性組合，使其變為一個平穩序列。例如，{xit}為一階單整的非平穩序列，如果它們滿足協整關系，則存在一組系數a1，a2，…，al，使得φ=a1x1t+a2x2t+…+alxlt成立，其中，φ為協整殘差。

3.2 協整分析建模步驟

對于非平穩時間序列，通常將其轉變為平穩序列再開展相關研究。時間序列單位根檢驗也稱為時間序列平穩性檢驗，常用方法是增廣迪基-富勒檢驗。若非平穩時間序列存在單位根，則采用差分方法得到平穩序列。

3.2.1關鍵特征的ADF單位根檢驗

考慮下列回歸模型：

xt=ρxt-1+et

(2)

式中，x0=0；t為觀測時刻；ρ為系數；et～(0，σ2)為均值是0、標準差是σ的白噪聲序列隨機誤差。

對于時間序列{xt}，若|ρ|<1，則{xt}是平穩的；若|ρ|=1，則{xt}是非平穩的，{xt}的單整階數為1并含有一個單位根。

變量yt的ADF檢驗方程為

Δyt=γyt-1+θ1Δyt-1+θ2Δyt-2+…+

θp-1Δyt-p+1+et

(3)

γ=ρ-1

其中，Δyt為yt的差分；{yt-1}為{yt}的一階滯后序列；θi為系數，i=1，2，…，p-1；p為滯后階數，通常采用赤池信息準則(Akaike information criterion，AIC)確定；Δyt-i為yt的i階差分；et為隨機誤差[18]。

使用最小二乘法對式(3)中的系數γ進行估計，通過檢驗系數γ的t統計量，從而判斷序列yt是否為一階單整序列。檢驗假設H0：γ=0(yt非平穩)；H1：γ<0(yt平穩)。計算系數γ的t檢驗統計量：

(4)

式中，σγ為γ的標準差。

由ADF分布臨界值表查得在給定顯著性水平下對應γ的臨界值。若tγ小于臨界值，則拒絕H0假設，即原序列不存在單位根，是平穩時間序列；否則原序列不平穩，則必須對其差分進一步檢驗單位根，建立相應的檢驗模型再次開展假設檢驗，若差分序列的t檢驗統計量小于臨界值，則拒絕H0假設，說明差分序列是平穩時間序列，原序列為一階單整序列。

3.2.2關鍵特征的Johansen協整檢驗

假設有一階單整分量的向量時間序列{yt}，它的p階向量自回歸(vector auto regression，VAR)可以表示為

yt=A1yt-1+A2yt-2+…+Apyt-p+ξt

(5)

式中，yt=(y1t，y2t，…，ymt)T是由m個一階單整非平穩變量構成的m×1維向量，t=1，2，…，T；yt-1，yt-2，…，yt-p是與yt對應的p個m×1維滯后項；A1，A2，…，Ap為m×m維參數矩陣；ξt是m×1維殘差向量。

式(5)兩邊同時減去yt-1完成差分變換，可得到對應于p階向量自回歸的誤差修正模型[11]：

(6)

(7)

式中，Em為單位矩陣；α為調整系數矩陣；β為協整系數矩陣；ξt～N(0，Ω)為平穩的殘差向量；Ω為ξt的方差。

α和β均為m×r維矩陣，其秩均為r。

整合式(6)、式(7)，可得

(8)

為確定協整關系數目r和協整系數矩陣β，建立兩個輔助回歸方程，并完成系數的最小二乘估計[19]：

(9)

(10)

利用殘差矩陣R0t、R1t構建積矩陣Swz(w，z=0，1)：

(11)

建立并求解特征值方程，可得協整系數矩陣β：

(12)

(13)

式中，λ為特征方程的解；βi(i=1，2，…，m)為λi的特征向量且λ1≥λ2≥…≥λm。

如果由r個最大的特征值給出了協整向量，則剩余m-r個非協整組合的特征值λr+1，λr+2，…，λm為0。以特征值的跡統計量Qr來檢驗變量的協整關系，給出Johansen檢驗假設H0：存在r(r=0，1，…，m-1)個協整關系。跡統計量Qr為[9]

(14)

按順序依次檢驗統計量Qr的顯著性，Qr大于檢驗臨界值時，拒絕H0假設，繼續按順序開展統計量顯著性檢驗，直到檢驗假設H0不被拒絕為止。確定協整關系數目r后，即可得到協整系數矩陣βr，由此建立協整殘差模型[20]：

(15)

ξt即為所求的協整殘差模型。

若殘差模型對應的r個殘差序列均為平穩時間序列，則都可用于描述系統的穩定性。這r個殘差序列中，最大特征值對應的模型包含系統信息最多，為最佳監測模型，稱為最優殘差序列。

3.3 協整殘差區間和預警閾值的確定

3.3.1基于SPC的殘差閾值確定方法

殘差處于控制限之內表明裝備運行正常；殘差超出控制限，則認為裝備運行異常。目前，利用最優殘差序列確定裝備狀態殘差區間和預警閾值的方法主要為：①以v±3σ為殘差閾值的上下限即方法1，其中，v、σ分別為殘差序列的均值和標準差；②以裝備正常運行狀態下殘差的最大幅值、最小幅值確定殘差閾值的上下限即方法2[8-11，13，15]。

上述兩種方法均可用于裝備狀態監測和提前預警，但也存在不足：方法1殘差閾值區間過大，提前預警效果較差，甚至會造成報警延誤；方法2沒有考慮裝備運行中可能存在的多種工況，具有局限性。為此，本文基于概率圖，提出考慮裝備多重工況以確定殘差閾值的第3種方法。

3.3.2基于概率圖的殘差閾值確定方法

概率圖是根據樣本的真實數據以及指定理論分布的累積概率所繪制的散點圖，能直觀檢測樣本數據是否符合某一概率分布。由統計學理論可知：若樣本總體服從正態分布，則它在散點圖上近似呈一條直線[21]。

裝備正常運行時，殘差序列多服從正態分布，具體表現為：均值附近的概率較大；離均值越遠，點越少、出現的概率越小。裝備出現異常時，系統會偏離正常狀態，輸出的殘差幅值也會偏離正常值，概率密度曲線將出現偏移，正常狀態下發生概率較小的數據點會大量出現[12]。因此，通過概率圖擬合得到的最優殘差序列，繪制多重工況概率圖，根據不同工況下殘差序列小概率點對應的殘差最小幅值、最大幅值，可以確定各工況對應的殘差分布區間及預警閾值，實現裝備多工況狀態監測與故障診斷。

4 案例分析

本節以某風電場3 MW直驅式風力機為研究對象，利用其2014年6月至12月的SCADA數據和故障維修數據進行分析，已知該風力機每隔10 min獲取一條狀態數據。根據零部件類型，故障可分為電機故障、控制系統故障、空冷系統故障、勵磁故障和饋電故障，故障統計如表1所示。

表1 故障次數及維修時間

由表1可知，電機故障導致的維修時間最長，因此，本文主要考慮風力機正常、偏航和電機故障3種狀態，其中，正常和偏航狀態下的電機狀況良好，但利用協整分析進行狀態監測與故障診斷時，偏航狀態會引起殘差幅值超過正常閾值，造成電機故障誤報。為排除這一干擾，本文提出一種裝備多工況監測與故障診斷方法，具體分析過程如下。

已知該風力機SCADA系統采集風速、轉速、功率、環境溫度等42個特征變量，利用隨機森林特征選擇算法提取反映電機故障的關鍵特征變量。由文獻[22]可知，轉速是能反映電機故障與否的有效特征變量。對SCADA數據進行預處理，篩除異常狀態數據，得到正常運行狀態下各特征變量的時間序列，以轉速作為隨機森林的輸出變量，將其余41個特征作為隨機森林的輸入變量，計算得到最重要的前4個特征變量：風速、功率、無功功率和定子溫度2。繪制上述關鍵特征的時間序列，如圖1所示。

(a)風力機風速時間序列 (b)風力機轉速時間序列 (c)風力機功率時間序列

(d)風力機無功功率時間序列 (e)風力機定子溫度2時間序列圖1 風力機正常運行狀態的SCADA數據Fig.1 SCADA data under normal operation of wind turbine

由圖1可知，風力機運行環境的變化使其特征變量的變化為典型的非平穩隨機變化。ADF檢驗表明上述5個特征變量均為非平穩時間序列，一階差分后可得到平穩時間序列，即原序列均為一階單整序列。通過AIC準則確定模型的滯后階數為5。利用Johansen協整檢驗確定協整關系的數目，如表2所示。

風力機的風速、轉速、功率、無功功率和定子溫度2之間存在4個協整關系，計算可以得到4個協整向量，由4列的協整向量組成的協整系數矩陣為

(16)

將協整系數矩陣β代入式(15)，并結合SCADA系統收集的關鍵特征運行數據，計算得到協整殘差模型：

表2 Johansen協整檢驗結果

(17)

Johanson協整檢驗是基于最大特征值的統計方法，特征值越大，對應模型包含的信息越多，監測系統運行狀態的能力越好。殘差模型的特征值是從大到小排列的，故選取特征值最大的殘差模型ξ1t進行歸一化處理，得到風力機運行狀態的最佳監測模型：

ξ1t=y1t-0.7794y2t-0.0022y3t+

0.0319y4t+0.0076y5t

(18)

將風力機的上述5個特征變量代入式(18)，構造最優殘差序列并開展ADF檢驗，得到檢驗值為-13.29，小于5%臨界值-1.9416，表明最優殘差序列為平穩時間序列，可用于評估風力機的運行狀態。

為排除因風速未達到切入風速導致的偏航狀態對風力機電機故障診斷的干擾，本文將風力機分為正常、偏航、電機故障3種狀態。通過數據預處理得到3種狀態對應的殘差序列，其中，電機故障狀態殘差序列包含故障數據、偏航數據和正常數據，偏航狀態殘差序列包含偏航數據和正常數據，正常狀態殘差序列只包含正常數據。繪制3種狀態的殘差序列概率圖(圖2)。

圖2 風力機3種狀態的殘差分布概率圖Fig.2 Probability plot ofresidual distribution of 3 states of wind turbine

風力機正常狀態的殘差區間為[-3.8957，3.9020]，偏航狀態的殘差區間為[-4.8234，7.7697]，故障狀態的殘差區間為[-4.8234，43.3767]。當風力機出現偏航和故障時，殘差下限變化不大，而殘差上限會出現明顯偏離，因此，本文通過分析3種狀態的殘差上限完成工況劃分。殘差幅值3.9020是風力機正常狀態與偏航狀態的分界線，殘差幅值7.7697是風力機偏航狀態與電機故障狀態的分界線，從而得到如下的劃分結果：風力機正常狀態的殘差區間為[-4.8234，3.9020]，偏航狀態的殘差區間為(3.9020，7.7697]，電機故障狀態的殘差區間為(7.7697，43.3767]。

取同一臺風力機運行狀態下的SCADA數據來驗證模型、殘差預警閾值的有效性。根據風力機SCADA系統采集的數據可知，在樣本點207附近，風速未達到切入風速，風力機轉速和功率很小，處于偏航狀態；在樣本點362附近，風速達到切入風速，但風力機轉速和功率很小，處于電機故障狀態。根據得到的最優殘差序列，利用基于概率圖的殘差閾值確定方法(方法3)進行分析，并與基于統計過程控制的殘差閾值確定方法(方法1和方法2)進行比較，結果如圖3所示。

圖3 不同殘差閾值確定方法的工況劃分Fig.3 Condition division of different residual threshold determination methods

由圖3可知，根據概率圖確定的電機故障預警閾值為7.7697，風力機在樣本點360.6時，殘差幅值超出方法3確定的閾值上限，此時風力機電機發生故障。而利用SCADA系統直接進行監測，直到樣本點362時系統才會發出電機故障警報。風力機每隔10 min記錄一個樣本點，獲取一條狀態數據，故與傳統方法相比，方法3能提前14 min監測到電機故障。根據殘差序列均值和標準差確定的電機故障預警閾值為13.2728，即方法1在樣本點361.2時發出警報，并能較準確地識別出電機故障，但報警提前期僅為8 min，提前預警能力不如方法3。根據正常狀態下殘差序列最大幅值確定的電機故障預警閾值為3.902，即方法2在樣本點360.3時發出警報，能提前17 min識別出電機故障，但不能有效排除樣本點207時偏航狀態造成的干擾，預警準確度不如方法3。

5 結論

本文基于風力機的SCADA數據，在協整分析的基礎上集成隨機森林和概率圖等理論，針對非平穩多工況的裝備提出一種狀態監測與故障診斷方法。本文在協整分析之前利用隨機森林特征選擇算法、提取關鍵特征，提高了協整殘差的精度。在協整分析之后提出基于概率圖的殘差閾值確定方法，有效實現了多種工況的狀態監測與提前預警，提高了故障診斷的準確性。

后續可結合不同類型工業裝備以及故障類型展開故障診斷方法研究。此外，可根據裝備不同月份、季度故障次數的統計，尋找其季節性特征，并通過調整裝備故障診斷的殘差閾值，開展維修策略優化研究。