董興輝,張勁草,李 佳,高 迪,楊志凌
(1.華北電力大學 能源動力與機械工程學院,北京 102206;2.北京市城市照明管理中心,北京 100078;3.中國電力科學研究院有限公司,北京 100192)
在一定氣象條件下,風電機組葉片會出現結冰現象,既帶來安全生產問題,又影響機組功率輸出。針對葉片結冰條件,以及結冰對風電機組運行參數的影響,有關學者開展了比較深入的研究,按照研究內容大致歸納為數值模擬與數值仿真兩種。
針對葉片霧凇、雨凇結冰,文獻[1]進行了不同程度的模擬,得出雨凇對葉片表面粗糙度、葉片氣動性能、風機轉矩以及輸出功率的影響均遠大于霧凇的結果。將三維葉片簡化為二維[2],雖然可以較好地模擬霧凇覆冰,但受雨凇覆冰明顯改變葉片形狀,激發振動等相應復雜變化,模擬計算還有待進一步優化。采用歐拉法模擬旋轉葉片覆冰過程[3],但模型的準確性缺乏試驗驗證。也有基于多參考坐標系(MRF)法改進歐拉模型用于模擬旋轉風機結冰[4],還有研究覆冰對葉片不同位置、不同攻角的氣動性能影響,但缺乏試驗驗證。文獻[5]研究了葉片覆冰下槳距角的變化,得出不同功率限定值下的“槳距角-風速”曲線。文獻[6]研究了葉片在不同覆冰厚度下的模態參數變化。
針對葉片結冰檢測方法的研究也分為兩類:氣象觀測檢測與基于數據的分析檢測。
氣象觀測檢測是通過分析大氣變化與葉片結冰的關系,間接判定葉片是否結冰;或者在機組機艙上安裝專門的傳感器,通過測量結冰引起的葉片物理特性變化,判斷葉片是否結冰。上述方法不僅增加了機械復雜度,而且增加維護成本[7]。
基于數據分析的檢測是根據專家經驗知識提取數據特征,利用隨機森林算法對葉片結冰進行識別[8]。但該檢測方法難以全面客觀地反應葉片結冰特性。文獻[9]基于深度自編碼網絡(DNN)學習SCADA故障特征,訓練幾個基分類器,然后通過集成學習組合構建覆冰檢測模型。也有基于葉片振動監測葉片覆冰狀態[10],根據覆冰概率與氣象參數的關系,通過建立葉片覆冰概率矩陣來評估機組出力損失[11]。構建葉片結冰前后風功率曲線,取得冰凍功率折減系數[12],對比理想功率曲線和溫度偏離預測葉片結冰,基于環境溫度、濕度、功率損耗建立結冰檢測模型[13]。隨著人工智能、大數據分析技術的發展,挖掘SCADA數據中可以表征風機葉片結冰的特征量[14],利用粒子群算法(PSO)優化SVR回歸模型,建立結冰檢測模型。
葉片結冰與季節和氣象條件密切相關。盡管不同地域存在結冰時節、結冰期長短等差異,但在結冰條件上,存在著溫度、風速、液態水含量等相近特性。葉片結冰的主要影響因素如表1所示。

表1 葉片結冰的主要影響因素Table1 Main influencing factors of blade icing
仿真與實測數據統計表明,葉片結冰溫度發生 在0~-11℃,風 速 在3.8~12m/s。
圖1為 風電 場一 在20170115T00:00-05:40時間內的SCADA系統數據(采樣頻率為10min,共計35個數據點)。現場運維記錄顯示,05:00出現機組葉片結冰。

圖1 葉片結冰前后部分SCADA參數變化特性Fig.1 Variation characteristics of SCADA parameters before and after blade icing
機組在不同工況下的運行參數值如表2所示。其中,風速V、葉片轉速n、功率P、槳距角 β可直接從SCADA中獲取。

表2 不同工況下機組參數特性Table2 Parameter characteristic of the WT under different working conditions
根據風速值,對 比表2,在 序列點[1~13]時 段(圖1)機組處于工況Ⅰ區,V,β,風能利用系數Cp值都處于正常狀態;在序列點[13~29]時段,風速降低,機組處于Cpmax控制狀態的工況Ⅱ區,V隨風速減小而降低屬于正常,各項參數值保持在正常范圍內;在序列點[29~33]時段,風速由5.7m/s持續增大到10m/s,接近額定風速,此時機組運行于Ⅱ區-Ⅲ區邊界工況區間,Cp還維持在Cpmax,β也在0°附近,V理應逐漸增大到滿發狀態,但從29點開始,Cp自0.43快速下降,V雖然在增大,但輸出功率達不到(滿發)狀態Ⅲ,說明機組在轉換風能方面出現異常。即從04:40時刻起葉片結冰。
運 行 在[Vmin,Vmax],[nmin,nmax]各 工 況 區 間 上 的 機組 氣 動 特 性 參 數 λmin與 λmax,Cpmin與Cpmax, 可 由 式(1),(2)計 算。

式中:λ為葉尖速比;R為葉片半徑。
風力機葉片表面結冰是隨機過程,呈非均勻分布,必然導致葉輪質量不平衡,機組塔架產生新的振動特性。圖2為風電場二某機組在20170319 T00:00-23:59,葉片結冰前后塔架振動數據時序圖(采樣頻率為1min,共計1440個數據點)。顯然,葉 片 結 冰 前(0~450)后(450~1200)直 至 停 機(1200點之后),塔架振動特征明顯不同。

圖2 風電機組塔架時頻振動Fig.2 WT tower time-frequency vibration
結冰增大了測風儀轉軸的摩擦力,致使其所測風速值偏小,開始出現測值誤差。相對于機械式風速儀在低溫環境中易受雨雪、凍雨等結冰影響,超聲波風速儀則不受低溫環境影響,有著更高的可靠性和準確性。因此,同一臺機組上的兩臺測風儀在結冰前后所測風速值之差就會發生變化。設:

結冰改變葉片表面翼型,直接導致其氣動性能下降,減小風能利用率,影響輸出功率。葉片結冰呈無規則形狀,必然引起葉輪不平衡,導致塔架產生新的振動特性。結冰期間的機械式風速儀也因結冰阻尼增大,測風精度下降。綜合葉片結冰狀態特性和結冰天氣條件,形成葉片結冰辨識邏輯如圖3所示,葉片結冰辨識流程如圖4所示。

圖3 葉片結冰與辨識邏輯Fig.3 Blade icing and identification logic

圖4 葉片結冰辨識流程圖Fig.4 Blade icing identification process
遴選與葉輪轉換特性、振動特性、輸出功率特性有關的參數,將非結冰機組歷史數據劃分為訓練集和測試集兩部分,運用支持向量機算法建立參數回歸模型,計算參數歷史回歸值。
依據皮爾遜積矩相關系數法的相關性分析,獲得機組葉片結冰狀態的伴隨參數。相關系數r是描述變量之間線性相關程度的量。

式中:Xi,Yi為變量樣本;,樣本平均值。
r的 取 值 為[-1,1]。當{r|-1≤r<0}時,負 相 關;當{r|0<r≤1}時,正相 關;當r=0時,不 相 關。r的數值越大,表明變量之間相關性越強。相關度分布如表3所示。

表3 相關度劃分Table3 Divide the correlation
對于給定的輸入-輸出樣本數據集 (xi,yi|i=1,2,…,N),其 中,xi為 第i個d維 輸 入 樣 本 變 量,yi為第i個輸出標量,N為訓練樣本個數。依照優化目標函數:

式 中:w,b為 超 平 面 參 數;C為 懲 罰 系 數;ξi,為松弛變量。
對于約束條件,SVR模型的回歸函數為

式 中:f(x)為 輸 入 向 量x對 應 的 回 歸 輸 出;αi,為
拉格朗日乘子。

式中:σ為核寬參數。
考慮到風電機組SCADA數據具有多樣本、多特征的特點,選擇徑向基核函數作為SVR核函數,構造支持向量機回歸函數。

式(5)中的C和式(7)中的 σ取值,決定著SVR模型的學習能力和預測能力。為了得到具有最佳泛化性能的SVR模型,本文采用PSO算法對模型中的參數(C,σ)進行優化,從中選出最佳的C與σ。基于SVR進行參數偏離回歸預測如圖5所示。

圖5 基于SVR算法的參數偏離回歸辨識Fig.5 Regression Identification of Parameter Deviation Based on SVR Algorithm
利用回歸函數可計算出正常狀態下的參數預測值。參數值偏離計算如下:

當 δ*偏 離 滿 足 條 件:δP>ηP&δA>ηA&δCp>ηCp(η*為偏離閾值),則可判定葉片出現結冰現象。
選取寧夏某風電場結冰機組數據驗證上述模型的正確性。機組主要參數:葉輪直徑為110m,額定功率為2100kW,切入風速為3m/s,額定風速為10.5m/s,SCADA數據采集頻率為1min。
4.1.1確定預測輸入參量
隨機取#A50機組連續5d的歷史數據作為樣本,依據皮爾遜積矩相關系數法,取|r|≥0.6(中、高相關程度)作為標準選取相關參數,計算各監測參數與有功功率的相關系數,得到變流器無功功率均值、發電量均值、風速1(機械式測風儀)均值、風速2(超聲波式測風儀)均值、30s風速均值、5min風速均值、輪轂轉速均值、發電機轉速均值、轉 矩 反 饋 平 均 值、電 網 電 壓AB(BC,CA)均值、電網A(B,C)相電流均值、電網線電壓L1L2(L2L3,L3L1)均 值、電 網 出 口 線 電 流L1(L2,L3)均 值 和 塔 架 振 動 加 速 度 均 值(波 段5,6,7,8,9,11,12)與有功功率輸出量相關的輸入參量。
考慮到機組相關變量數據具有不同單位和量級,且指標之間的量值差異非常大,為保證預測結果的準確性,需要對原始數據進行歸一化[0,1]處理。
4.1.2選取驗證數據
運行記錄顯示,#A50機組在20170319出現葉片結冰現象,并在19:56停機。由于沒有安裝任何監測葉片結冰的軟、硬件設備,無法準確確定風力機葉片開始結冰的時間。鑒于此,本文選取機組SCADA系 統 在20170319T00:00-19:56的 監 測 數據作為分析、驗證數據。
4.1.3輸出功率數據辨識結冰
圖6(a)為SCADA記錄的有功功率1440個數據點以及應用回歸預測計算出的相應功率值,圖6(b)為 其 部 分 時 段 數 據。從 圖6(b)可 知,大 約在460點之后,預測值開始偏離SCADA記錄值。

圖6 有功功率預測值與原始值對比Fig.6 Compare the predicted value of active power with the original value
圖7為機組有功功率的相對誤差曲線。由圖7可知,自458點之后,功率損失持續上升,即20170319T07:37前后發生“可能因葉片結冰”導致的機組整體性能“異常”。

圖7 有功功率400~600序列點相對誤差Fig.7 Active power400~600sequence point relative error
4.1.4Cp數據辨識結冰
根據所選SCADA驗證數據,計算機組Cp時序 值(圖8)。
對比表3工況劃分與圖8,#A50機組處于兩種 狀 態,[1~342]時 序 點 間 處 于 第 一 狀 態,[343~1440]時序點間處于第二狀態。
[1~342]時 序 點 間,風 速 為10.5m/s,機 組 處 于Ⅳ工況區,葉輪轉速、Cp,β正常波動變化,機組正常運行。
[343~1440]時 序 點 間,風 速 由10.5m/s逐 步降 低 到4.198m/s,之 后[443~1440]時 序 點,風 速 持續波動在4.198~10m/s,正常狀態下機組應處于Ⅱ工況區,即Cpmax階段,β不變,發電機轉速跟隨風速而變化。
圖9為[400~600]數據時序點對應的多個參數值。

圖9 時序400~600點間參數值Fig.9 Parameter values of sequence points in 400~600interval
由圖9可知,[442~600]區段風速從4.198m/s開始持續增大,機組理應處于Cpmax,但與[400~452]區段不同的是,圖中自453點開始,Cp出現快速下降,表明機組于時序453點即20170319T07:32時刻起,開始出現氣動方面“異常”。
4.1.5振動特征數據辨識結冰
將塔架振動(振動加速度平均值)信號數據應用EEMD分解,選取前3個IMF特征信號進行分析(圖10)。

圖10 EEMD信號特征分解圖Fig.10 EEMD signal characteristics decomposition diagram
由圖10可知,3種不同頻率信號在441點附近的幅值都各自發生突然增大現象,即機組塔架于20170319T7:20左右,出現塔架振動"異常"現象。
上述3種模型計算結果對比如表5所示。

表5 3種SCADA參量對比判定機組葉片結冰Table5 Comparison of three SCADA parameters to determine the blade icing of WT
表5中,3種參量計算結果表明,#A50機組大約從450點開始出現運行異常,即葉片結冰。相對于機組輸出功率來說,振動信號更加敏感,感知葉片結冰時刻更早一些。
機組#A50共安裝了機械式和超聲波兩種類型風速儀,兩種測速儀實測風速數據如圖11所示。

圖11 兩種風速測量儀的測量值對比Fig.11 Comparison of measured values between the two kinds of anemometer
由圖11(a)可知:正常情況下,機械式測速儀和超聲波測速儀所測得的風速值差別不大、趨勢一致;在序列點[442~527]區域附近,葉片結冰之前,機械式比超聲波式所測風速值偏大,但結冰之后則相反[圖11(b)]。不考慮測速儀發生機械類故障,則可以推定是受結冰摩擦,影響到機械測量值。這也間接印證了#A50機組于時間序列第450點附近出現了結冰現象。
①與通過觀察機組功率曲線這種單一檢測手段相比,本文采用多參數組合驗證方法,可在葉片出現結冰的30min內辨識出葉片結冰,便于及時采取相應的運維決策。
②對比分析機組輸出功率、風能利用系數和塔架振動信號3種參量的異常時序,機組塔架振動信號更加靈敏,可以更早反應葉片表面出現覆冰現象。
③通過分析葉片結冰前后SCADA有關參數的變化特性而不是通過增加硬件的檢測方法,避免了因增添設備而改變葉片的氣動特性、增加維護維修復雜度和運行成本。