賈子文, 顧煜炯
(華北電力大學 能源動力與機械工程學院, 北京 102206)
風電機組的安全穩(wěn)定高效運行是風力發(fā)電行業(yè)的基本要求。但是風電場選址大多在偏遠地區(qū),惡劣的環(huán)境和非穩(wěn)定負荷變化等因素對機組的正常運行帶來較大影響[1]。隨著風電機組運行時間的增加,其維護費用呈逐年上升的趨勢。當機組服役時間超過20年時,其維護成本占能源成本的15%左右[2]。風力發(fā)電行業(yè)前期投入大,資本回籠慢,過高的運維開銷是風電行業(yè)發(fā)展的重大隱患。因此,對風電機組設(shè)備運行狀態(tài)進行實時監(jiān)測、及時準確診斷機組故障是避免嚴重故障發(fā)生和降低事后維修費用的有效途徑。
風電機組作為大型戶外發(fā)電設(shè)備,風速、溫度、雨水和沙塵等自然因素均會對機組運行狀態(tài)造成影響。運行工況的時變性增大了診斷風電機組故障的難度。在既定工況下通過單純的數(shù)據(jù)特征分析已不能滿足風電機組故障診斷的要求[3]。盧錦玲等[4]對極限學習機方法進行改進,將風電機組主軸軸承振動信號作為模型輸入,研究了主軸軸承的故障診斷問題。Kusiak等[5]根據(jù)實時風速對機組運行工況進行劃分,確定了監(jiān)測各個工況振動數(shù)據(jù)的閾值。向玲等[6]將小波包與固有時間尺度分解方法相結(jié)合,通過振動信號主要特征分量的關(guān)聯(lián)維數(shù)計算實現(xiàn)了風電機組齒輪箱的故障診斷。但以上研究大多僅對機組某一類型的運行參數(shù)進行分析,單一數(shù)據(jù)特征有時并不能充分反映機組運行工況的時變性與多樣性。李靜立等[7]將模態(tài)參數(shù)識別與階次分析應用到風電機組齒輪箱狀態(tài)監(jiān)測領(lǐng)域,完成了在動態(tài)環(huán)境下齒輪箱運行狀態(tài)的識別工作。Zhang等[8]和王雅琳等[9]采用多塊核主元分析(MBKPCA)方法對運行數(shù)據(jù)進行分塊處理,實現(xiàn)了設(shè)備故障發(fā)生位置的診斷。但是以上研究對于運行數(shù)據(jù)劃分的主觀性過強,缺乏對變量之間關(guān)聯(lián)性的分析。郭鵬等[10]采用非線性狀態(tài)估計方法對風電機組齒輪箱油溫趨勢進行分析,提高了齒輪箱診斷結(jié)果的精度,雖然考慮了機組在確定工況下運行過程的組成,但是在進行樣本數(shù)據(jù)遴選過程中沒有形成樣本數(shù)據(jù)與運行過程的對應關(guān)系,對數(shù)據(jù)特征沒有進行明確說明。
針對以上問題,筆者提出一種基于因子分析改進的MBKPCA方法,可實現(xiàn)風電機組的狀態(tài)監(jiān)測與診斷。采用對應分析方法對風電機組的樣本數(shù)據(jù)與變量進行關(guān)聯(lián)性分析,確定機組在正常狀態(tài)下各運行過程對變量的影響程度,明確數(shù)據(jù)分塊份數(shù)及每個子塊的實際意義,可提高各子塊分析結(jié)果的解釋能力;采用因子分析方法,通過對機組各運行狀態(tài)下數(shù)據(jù)特征貢獻率的計算,確定正常情況下的典型數(shù)據(jù)并作為最終診斷模型的樣本,提高MBKPCA方法的準確性。
核主元分析法(KPCA)[11]是在傳統(tǒng)主元分析法(PCA)基礎(chǔ)上引入核函數(shù)概念的方法,通過非線性函數(shù)將輸入的低維數(shù)據(jù)映射到高維空間中進行處理,用核函數(shù)運算代替非線性變換后進行特征空間內(nèi)積運算,可減少計算量。假設(shè)X={X1,X2,…,Xn}作為KPCA的變量輸入序列,其中每個變量有m個樣本??紤]到機組運行數(shù)據(jù)的非線性特點,采用非線性映射Φ:RM→F,將原始輸入數(shù)據(jù)在高維特征空間中展開。其中,特征空間F的協(xié)方差矩陣為:
(1)
式中:n為參與分析的變量序列長度。
(2)
式中:αi為核矩陣K的特征向量中第i個值。
引入核函數(shù)Kjk=K(Xj,Xk)=Φ(Xj)Φ(Xk),表示核矩陣K中第j行、第k列的值,K的特征值可表示為:
λα=(1/n)Kα
(3)
(4)
式中:E為單位矩陣。
通過主元分析中的累計方差貢獻率 (取80%),確定主元個數(shù)p,即原始樣本的最終特征表示為1個p維向量:
(5)
新觀測樣本XT的得分向量為:
se=[se1,se2,…,sep]
(6)
(7)
對于新的觀測數(shù)據(jù),可通過監(jiān)控統(tǒng)計量T2和SPE進行監(jiān)測,其表達式如下:
(8)
(9)
式中:T2為觀測數(shù)據(jù)特征與樣本數(shù)據(jù)特征的重合程度;SPE為觀測數(shù)據(jù)與樣本數(shù)據(jù)的偏差程度;Λ為每個主元特征值形成的對角矩陣。
(10)
MBKPCA通過對輸入數(shù)據(jù)進行分塊的方式來解釋設(shè)備各組成之間存在的固有聯(lián)系,在一定程度上提高了最終分析結(jié)果的精度,但是在整個分塊過程中存在以下問題:
(1) 通常直接指定或依據(jù)變量采集位置與設(shè)備結(jié)構(gòu)之間的關(guān)系給出數(shù)據(jù)分塊份數(shù)。分塊劃分過程缺乏依據(jù),各子塊表示的實際含義不明確,主觀性較強。
(2) 傳統(tǒng)的MBKPCA方法缺乏對典型樣本選取的過程。因為樣本數(shù)據(jù)的特征過于分散會影響子塊數(shù)據(jù)的特性,弱化各子塊數(shù)據(jù)特征的差異性。分塊策略雖然削減了原始變量核矩陣的尺寸,但對每個子塊核矩陣的分析仍屬于高維空間運算。原始數(shù)據(jù)量較大時,MBKPCA運算效率也會受到影響。
針對以上問題,筆者采用對應分析方法建立各變量與設(shè)備運行工況之間的關(guān)聯(lián)機制,通過樣本數(shù)據(jù)與各個變量的對應分布明確分塊份數(shù),提高各個子塊中數(shù)據(jù)的解釋能力;應用因子分析方法,合理提取各子塊中具有典型特征的數(shù)據(jù),提升子塊間的特征差異,減少樣本總數(shù),提高運算效率。
1.2.1 基于對應分析的分塊原則
對應分析方法[13]主要通過對定性變量列聯(lián)表解析,揭示變量之間、變量和樣品(指工況)之間的關(guān)系。此方法適用于數(shù)量型變量和品質(zhì)型變量,并且可通過圖形直觀反映變量與樣品之間的關(guān)系,便于后續(xù)分析結(jié)果的揭示與推理。結(jié)合MBKPCA方法,可觀察變量與設(shè)備運行過程的分布情況,以確定分塊份數(shù),并根據(jù)子塊中變量的內(nèi)容賦予子塊一定的含義,具體過程如下:
(1) 根據(jù)輸入變量結(jié)構(gòu),構(gòu)建輸入變量矩陣X=(xij)n×m,計算規(guī)格化概率矩陣P=(pij)n×m。
(11)
(2) 計算過渡矩陣Z=(zij)n×m,zij可表示為:
(12)
(13)
(14)
(4) 繪制變量點與樣本點分布圖。載荷矩陣U和V中各元素取值范圍相同,元素含義相近,可將載荷矩陣中的數(shù)據(jù)看成二維點繪制在同一坐標平面,以各點對應的因子載荷值作為坐標值,形成對應分布圖。
(5) 在分布圖中,載荷矩陣U的點表示變量的分布,V中的點表示樣本的分布,即點分布表示變量與樣本之間的對應關(guān)系。通過對橫縱坐標區(qū)間的選取,可提取變量與樣本集中的點集。其中,點集表示被測對象當前狀態(tài)下的某一種子狀態(tài),點集中樣本表示此子狀態(tài)的表征數(shù)據(jù),變量表示影響這一子狀態(tài)樣本數(shù)據(jù)的主要因素。例如,輸入風電機組正常數(shù)據(jù),則點集表示設(shè)備在機組正常情況下的運行過程。其中,點集中樣本數(shù)據(jù)對應的是MBKPCA中子塊的輸入數(shù)據(jù),區(qū)域中變量數(shù)據(jù)表示的是該子塊對應運行過程中的主導變量,點集個數(shù)即為分塊份數(shù)。可進一步綜合各子集中變量的物理意義,明確各子集的含義,即子塊的物理意義。
1.2.2 基于因子分析的典型樣本提取
通過對原始數(shù)據(jù)的對應分析,確定了MBKPCA分塊份數(shù)和各運行過程的主導變量。但只說明輸入樣本數(shù)據(jù)與各運行過程的分布特性,對數(shù)據(jù)特征與運行過程特性之間的關(guān)系沒有做出明確解釋。如上所述,在已確定的子塊中過于分散的數(shù)據(jù)特征分布會削弱子塊對新數(shù)據(jù)的分析能力。在高維空間中對過多的數(shù)據(jù)進行分析,會增加整個算法的運算負擔。因此,筆者運用因子分析[14]方法,提取每個子塊中的典型數(shù)據(jù)作為最終樣本,以提高子塊的解釋能力與運算效率。具體過程如下:
(1) 對子塊中的數(shù)據(jù)變量進行標準化處理。
(2) 建立變量間相關(guān)系數(shù)矩陣,對其進行取樣適切性量數(shù)(KMO)和巴特利特球度檢驗,驗證所選數(shù)據(jù)進行因子分析的可行性。
(3) 采用主元分析法計算相關(guān)矩陣的特征值和特征向量,通過參數(shù)累計方差貢獻率確定相關(guān)矩陣的主要成分,并對主要成分進行因子旋轉(zhuǎn),提高主要成分間的差異性。
(4) 利用回歸方法計算子塊中數(shù)據(jù)在各主要成分下的得分,提取主要成分得分較高的數(shù)據(jù)作為樣本數(shù)據(jù)。
重復以上步驟,完成所有子塊數(shù)據(jù)的選取,形成最終的MBKPCA輸入樣本數(shù)據(jù)。
以1.5 MW雙饋風電機組為研究對象建立故障診斷模型。根據(jù)風電機組結(jié)構(gòu)及運行特點,確定機組正常情況下的模型輸入變量;通過對應分析確定MBKPCA分塊份數(shù);通過因子分析方法確定樣本集;確定機組正常情況下監(jiān)測統(tǒng)計量的控制上限,作為衡量機組運行是否正常的標準。
按照風電機組結(jié)構(gòu)和功能進行系統(tǒng)邊界劃分,可分為變槳系統(tǒng)、偏航系統(tǒng)、傳動系統(tǒng)和發(fā)電機系統(tǒng)。目前絕大部分機組配備數(shù)據(jù)采集與監(jiān)視控制系統(tǒng)(SCADA),對機組運行狀態(tài)進行實時監(jiān)測。篩選機組SCADA主要運行參數(shù),選擇以下參數(shù)作為診斷模型的輸入:變槳系統(tǒng)主要通過調(diào)節(jié)槳葉實現(xiàn)機組最大風能跟蹤和防止機組過載等功能,主要運行參數(shù)包括槳葉角度和風速;偏航系統(tǒng)主要根據(jù)實時風向改變機艙位置,保證機艙對風準確,主要運行參數(shù)為機艙位置和風向角度;傳動系統(tǒng)主要完成轉(zhuǎn)軸升速和傳遞轉(zhuǎn)矩的任務,運行參數(shù)以油溫數(shù)據(jù)為主;發(fā)電機系統(tǒng)將旋轉(zhuǎn)機械能轉(zhuǎn)化為電能,主要參數(shù)為有功功率。
通過現(xiàn)場調(diào)研,發(fā)現(xiàn)SCADA系統(tǒng)中機組傳動鏈測點較少。為滿足風電機組傳動系統(tǒng)狀態(tài)監(jiān)測與診斷的要求,在機組傳動鏈特定位置增加振動測點,具體測點位置及描述如表1所示。
表1 測點位置及描述Tab.1 Arrangement and description of measuring points
對傳動鏈主要組成部件進行結(jié)構(gòu)及運行特點分析,進一步提取振動數(shù)據(jù)中的有效信息并作為模型輸入變量。
齒輪箱結(jié)構(gòu)緊湊,各測點間信號耦合性較強,特征數(shù)據(jù)易受其他信號干擾,選擇無量綱因子[15]范疇中的波形裕度和偏態(tài)因子,其變量數(shù)值只與設(shè)備結(jié)構(gòu)本身有關(guān),與運行工況無關(guān)?;趫D形學的計算方法可有效避免其他信號干擾對結(jié)果帶來的影響。
發(fā)電機前、后軸承發(fā)生故障時,信號特征十分微弱,不易觀察,選擇J散度和KL散度[16]作為軸承部件的輸入變量。
為使模型輸入變量能全面反映機組的運行情況,從機組能效角度將各系統(tǒng)效率指標也作為模型的輸入變量,主要包括風輪損失、齒輪箱損失和發(fā)電機損失[17]。各損失均與主軸轉(zhuǎn)速或主軸轉(zhuǎn)速的倍數(shù)有關(guān),所以在最終的變量信息統(tǒng)計中,未獨立給出主軸轉(zhuǎn)速。模型輸入變量如表2所示。
表2 模型輸入變量Tab.2 Input variables of model
基于數(shù)據(jù)本身特性和數(shù)據(jù)來源的差異,信號采樣頻率以所有變量中采樣頻率最低的數(shù)據(jù)為準。其中,由SCADA系統(tǒng)獲得的機組各類運行數(shù)據(jù)均為該數(shù)據(jù)10 min的平均值,其采樣頻率最低。對于振動參數(shù),為保證數(shù)據(jù)分析的同步性,計算振動每個采樣時刻數(shù)據(jù)的無量綱因子和散度指標,并取其10 min的平均值作為模型輸入數(shù)據(jù)。
采集風速為11~12 m/s時機組的700組正常數(shù)據(jù)作為輸入變量,表3給出了協(xié)方差矩陣R的前4個特征值、方差貢獻率和累計方差貢獻率。
表3 協(xié)方差矩陣數(shù)據(jù)統(tǒng)計Tab.3 Covariance matrix data statistics
由表3可知,輸入變量前2個特征值的累計方差貢獻率已達到81.1%,說明這2個特征值已經(jīng)能夠表述所有變量的絕大部分信息,故提取前2個特征值,即確定公共因子個數(shù)為2。構(gòu)建R型載荷矩陣U和Q型載荷矩陣V,并根據(jù)載荷矩陣中的數(shù)據(jù)繪制輸入樣本與變量的對應分布,如圖1所示。
圖1 輸入樣本與變量的對應分布Fig.1 Corresponding distribution of input samples and variables
由于第1個公共因子的方差貢獻率達62.2%,遠高于第2個公共因子,因此圖1中以橫坐標為分類標準。變量與樣本大致可劃分為3個點集,即將原始數(shù)據(jù)劃分為3個子塊。每個子塊中主導變量分別為:子塊1變量,包括槳葉角、機艙位置、風向角度、風速和風輪損失;子塊2變量,包括波形裕度、油溫、偏態(tài)因子和齒輪箱損失;子塊3變量,包括J散度、KL散度、有功功率和發(fā)電機損失。
子塊1中的變量主要涉及風電機組能量來源和風能捕獲環(huán)節(jié),將此子塊命名為“能量子塊”;子塊2中的變量主要與機械能傳遞環(huán)節(jié)相關(guān),將此子塊命名為“傳動子塊”;子塊3中的變量主要與機組發(fā)電環(huán)節(jié)有關(guān),將此子塊命名為“發(fā)電子塊”。通過機組正常工況下數(shù)據(jù)點集分布,發(fā)現(xiàn)機組整個正常運行狀態(tài)由多個運行過程(子塊)組成,且不同運行過程對變量類型的影響也存在差異。因此,數(shù)據(jù)的對應分析不僅對每個子塊的物理意義作出解釋,同時說明其運行狀態(tài)特性,使后續(xù)分析結(jié)果的意義更為明確。
為提高各子塊對機組運行過程的解釋能力,并改善模型的運算效率,對每個子塊中的數(shù)據(jù)進行因子分析,尋找具有機組運行過程典型特征的數(shù)據(jù)作為最終樣本。建立每個子塊中各變量間的相關(guān)系數(shù)矩陣,并進行相關(guān)性分析,結(jié)果如表4所示。
表4 KMO和巴特利特球度檢驗Tab.4 KMO and bartlett's test
各子塊的KMO值均大于0.7,說明各變量之間關(guān)聯(lián)性較強。各子塊的巴特利特球狀檢驗概率均為0.00,小于顯著性水平0.05,說明樣本充足,可進行因子分析。以傳動子塊為例,計算該子塊相關(guān)系數(shù)矩陣的特征值、方差貢獻率和累計方差貢獻率,結(jié)果如表5所示。
表5 各成分原有指標總方差情況Tab.5 Original indicator and total variance of various components
表5中前2個成分的累計方差貢獻率達89.00%,由累計方差貢獻率要求(大于80%)可知,這2個成分可表示矩陣中數(shù)據(jù)的大部分信息,故確定主要成分個數(shù)為2。計算主要成分載荷矩陣,并對其進行旋轉(zhuǎn)計算。筆者采用Kaiser方差最大旋轉(zhuǎn)方法計算主要成分與變量間的載荷值,結(jié)果如表6所示。
表6 主要成分載荷值Tab.6 Load value of main components
變量間的載荷越高,其在對應成分中包含的信息就越多。成分1中的主要參數(shù)與齒輪箱振動相關(guān),此成分稱為“振動成分”。成分2中的主要參數(shù)與齒輪箱油溫和損失有關(guān),該成分稱為“溫度成分”。采用回歸方法計算各成分在樣本數(shù)據(jù)中的得分,結(jié)果如表7所示。
表7 各成分得分Tab.7 Component score
樣本得分越高,其在對應成分中的解釋能力越強,樣本的數(shù)據(jù)特征越具有代表性。因此,按照一定比例提取所有主要成分中的高分樣本數(shù)據(jù),形成傳動子塊的典型樣本。
同理,對其他2個子塊中的數(shù)據(jù)做相同分析,得出各子塊的典型樣本,形成最終參與MBKPCA分析的輸入樣本。樣本比例可依照子塊中各主要成分方差貢獻率進行選取。經(jīng)因子分析提取出的風電機組樣本與變量的對應分布如圖2所示。
圖2 因子分析樣本與變量的對應分布Fig.2 Distribution of samples and variables by factor analysis
為驗證風電機組故障診斷模型的有效性,以國內(nèi)某1.5 MW雙饋風電機組齒輪箱磨損故障為例進行分析。該機組在2013年2月3日出現(xiàn)解列停機情況,通過查找文件傳輸協(xié)議(FTP)服務器記錄,
圖3 故障診斷流程圖
Fig.3 Flow chart of fault diagnosis
風電機組齒輪箱油溫過高故障,經(jīng)現(xiàn)場勘查,發(fā)現(xiàn)齒輪箱一級太陽輪磨損嚴重。提取停機前7天的歷史數(shù)據(jù),利用所提出的方法對3個子塊單元進行監(jiān)測統(tǒng)計量T2和SPE分析,結(jié)果如圖4所示。
(a) 能量子塊統(tǒng)計量結(jié)果
(b) 傳動子塊統(tǒng)計量結(jié)果
(c) 發(fā)電子塊統(tǒng)計量結(jié)果圖4 改進的MBKPCA方法的故障檢測圖Fig.4 Fault detection plot by improved MBKPCA method
由圖4(b)可知,從第600點開始傳動子塊的2個統(tǒng)計量顯著上升,超過此子塊的控制上限,并且數(shù)據(jù)具有較好的識別持續(xù)性。由圖5可知,傳動子塊的統(tǒng)計值明顯高于其他子塊,且超出閾值,說明機組的傳動部分異常。傳動子塊主要變量為齒輪箱數(shù)據(jù),因此診斷為風電機組齒輪箱出現(xiàn)了故障。由圖4(a)可知,從600點開始能量子塊也出現(xiàn)了統(tǒng)計量數(shù)據(jù)上升的現(xiàn)象,這是因為主軸與風輪依靠法蘭盤連接,屬于剛性連接,為齒輪箱振動能量的傳遞提供了路徑。當齒輪箱振動較大時,振動能量可通過主軸傳遞到風輪,造成風輪數(shù)據(jù)特征偏離正常運行狀態(tài)。但風輪不是故障源頭,所以偏離程度較小,如圖5所示。齒輪箱輸出軸與發(fā)電機軸之間靠柔性聯(lián)軸器連接,柔性聯(lián)軸器不僅補償了兩軸平衡性偏差和角度誤差,較強的阻尼性可有效減弱由齒輪箱傳遞的振動能量,因此齒輪箱故障引起的振動對發(fā)電機運行無太大影響,發(fā)電子塊的統(tǒng)計值最低。綜上說明所提出的方法能夠有效診斷機組故障,并明確故障發(fā)生的具體位置。
(a) 各子塊SPE的統(tǒng)計量
(b) 各子塊T2的統(tǒng)計量圖5 改進的MBKPCA方法的監(jiān)測統(tǒng)計量Fig.5 Monitoring statistics by improved MBKPCA method
應用傳統(tǒng)MBKPCA方法,樣本數(shù)據(jù)分塊按照機組組成結(jié)構(gòu)進行劃分,分別為偏航子塊、變槳子塊、傳動子塊和發(fā)電子塊4部分。對同一組故障數(shù)據(jù)進行分析,結(jié)果如圖6所示。
(a) 偏航子塊統(tǒng)計量結(jié)果
(b) 變槳子塊統(tǒng)計量結(jié)果
(c) 傳動子塊統(tǒng)計量結(jié)果
(d) 發(fā)電子塊統(tǒng)計量結(jié)果圖6 傳統(tǒng)MBKPCA方法的故障監(jiān)測圖Fig.6 Fault detection plot by original MBKPCA method
由圖6可知,在650點處傳動子塊統(tǒng)計值超過控制上限,在時間上明顯滯后于改進的MBKPCA方法,且數(shù)據(jù)波動較大。由于傳統(tǒng)MBKPCA方法缺少對數(shù)據(jù)之間關(guān)聯(lián)性的分析,將能量子塊拆分為變槳子塊和偏航子塊,使統(tǒng)計量的趨勢變化相比改進的MBKPCA方法不明顯,2子塊的統(tǒng)計量變化與傳動子塊未形成對應關(guān)系。由圖7可知,相比改進的MBKPCA方法,傳統(tǒng)MBKPCA方法各子塊的統(tǒng)計量差異性較小,偏航子塊、變槳子塊和傳動子塊的監(jiān)測統(tǒng)計量相近,且均接近控制上限,無法明確機組故障發(fā)生的確切部位。這說明改進的MBKPCA方法能更及時準確地對風電機組故障作出診斷。
(a) 各子塊SPE的統(tǒng)計量
(b) 各子塊T2的統(tǒng)計量圖7 傳統(tǒng)MBKPCA方法的監(jiān)測統(tǒng)計量Fig.7 Monitoring statistics by original MBKPCA method
(1) 深度挖掘SCADA數(shù)據(jù)與振動數(shù)據(jù)信息,引入機組效率、無量綱因子變量和散度指標,提高了數(shù)據(jù)特征對機組故障的識別能力。
(2) 通過對應分析對傳統(tǒng)MBKPCA方法進行改進,明確了機組正常工況下各運行過程的特征,賦予了各子塊實際的物理意義,建立了數(shù)據(jù)與變量之間的關(guān)聯(lián)體系,增強了分塊劃分的客觀性。
(3) 采用因子分析方法分析數(shù)據(jù)與各運行過程之間的影響程度,提取對每個子塊影響程度較大的數(shù)據(jù)作為改進MBKPCA方法的輸入樣本,提高了診斷模型的精度與解釋能力。