基于數(shù)據(jù)特性分析的風(fēng)電機(jī)組葉片結(jié)冰辨識(shí)

董興輝，張勁草，李 佳，高 迪，楊志凌

（1.華北電力大學(xué) 能源動(dòng)力與機(jī)械工程學(xué)院，北京 102206；2.北京市城市照明管理中心，北京 100078；3.中國電力科學(xué)研究院有限公司，北京 100192）

0 引言

在一定氣象條件下，風(fēng)電機(jī)組葉片會(huì)出現(xiàn)結(jié)冰現(xiàn)象，既帶來安全生產(chǎn)問題，又影響機(jī)組功率輸出。針對(duì)葉片結(jié)冰條件，以及結(jié)冰對(duì)風(fēng)電機(jī)組運(yùn)行參數(shù)的影響，有關(guān)學(xué)者開展了比較深入的研究，按照研究內(nèi)容大致歸納為數(shù)值模擬與數(shù)值仿真兩種。

針對(duì)葉片霧凇、雨凇結(jié)冰，文獻(xiàn)[1]進(jìn)行了不同程度的模擬，得出雨凇對(duì)葉片表面粗糙度、葉片氣動(dòng)性能、風(fēng)機(jī)轉(zhuǎn)矩以及輸出功率的影響均遠(yuǎn)大于霧凇的結(jié)果。將三維葉片簡化為二維[2]，雖然可以較好地模擬霧凇覆冰，但受雨凇覆冰明顯改變?nèi)~片形狀，激發(fā)振動(dòng)等相應(yīng)復(fù)雜變化，模擬計(jì)算還有待進(jìn)一步優(yōu)化。采用歐拉法模擬旋轉(zhuǎn)葉片覆冰過程[3]，但模型的準(zhǔn)確性缺乏試驗(yàn)驗(yàn)證。也有基于多參考坐標(biāo)系（MRF）法改進(jìn)歐拉模型用于模擬旋轉(zhuǎn)風(fēng)機(jī)結(jié)冰[4]，還有研究覆冰對(duì)葉片不同位置、不同攻角的氣動(dòng)性能影響，但缺乏試驗(yàn)驗(yàn)證。文獻(xiàn)[5]研究了葉片覆冰下槳距角的變化，得出不同功率限定值下的“槳距角-風(fēng)速”曲線。文獻(xiàn)[6]研究了葉片在不同覆冰厚度下的模態(tài)參數(shù)變化。

針對(duì)葉片結(jié)冰檢測(cè)方法的研究也分為兩類：氣象觀測(cè)檢測(cè)與基于數(shù)據(jù)的分析檢測(cè)。

氣象觀測(cè)檢測(cè)是通過分析大氣變化與葉片結(jié)冰的關(guān)系，間接判定葉片是否結(jié)冰；或者在機(jī)組機(jī)艙上安裝專門的傳感器，通過測(cè)量結(jié)冰引起的葉片物理特性變化，判斷葉片是否結(jié)冰。上述方法不僅增加了機(jī)械復(fù)雜度，而且增加維護(hù)成本[7]。

基于數(shù)據(jù)分析的檢測(cè)是根據(jù)專家經(jīng)驗(yàn)知識(shí)提取數(shù)據(jù)特征，利用隨機(jī)森林算法對(duì)葉片結(jié)冰進(jìn)行識(shí)別[8]。但該檢測(cè)方法難以全面客觀地反應(yīng)葉片結(jié)冰特性。文獻(xiàn)[9]基于深度自編碼網(wǎng)絡(luò)（DNN）學(xué)習(xí)SCADA故障特征，訓(xùn)練幾個(gè)基分類器，然后通過集成學(xué)習(xí)組合構(gòu)建覆冰檢測(cè)模型。也有基于葉片振動(dòng)監(jiān)測(cè)葉片覆冰狀態(tài)[10]，根據(jù)覆冰概率與氣象參數(shù)的關(guān)系，通過建立葉片覆冰概率矩陣來評(píng)估機(jī)組出力損失[11]。構(gòu)建葉片結(jié)冰前后風(fēng)功率曲線，取得冰凍功率折減系數(shù)[12]，對(duì)比理想功率曲線和溫度偏離預(yù)測(cè)葉片結(jié)冰，基于環(huán)境溫度、濕度、功率損耗建立結(jié)冰檢測(cè)模型[13]。隨著人工智能、大數(shù)據(jù)分析技術(shù)的發(fā)展，挖掘SCADA數(shù)據(jù)中可以表征風(fēng)機(jī)葉片結(jié)冰的特征量[14]，利用粒子群算法（PSO）優(yōu)化SVR回歸模型，建立結(jié)冰檢測(cè)模型。

1 葉片結(jié)冰的有關(guān)參量變化特性

1.1 葉片結(jié)冰的必要條件

葉片結(jié)冰與季節(jié)和氣象條件密切相關(guān)。盡管不同地域存在結(jié)冰時(shí)節(jié)、結(jié)冰期長短等差異，但在結(jié)冰條件上，存在著溫度、風(fēng)速、液態(tài)水含量等相近特性。葉片結(jié)冰的主要影響因素如表1所示。

表1 葉片結(jié)冰的主要影響因素Table1 Main influencing factors of blade icing

仿真與實(shí)測(cè)數(shù)據(jù)統(tǒng)計(jì)表明，葉片結(jié)冰溫度發(fā)生 在0～-11℃，風(fēng) 速 在3.8～12m/s。

1.2 葉片結(jié)冰SCADA參數(shù)變化特性

圖1為 風(fēng)電 場(chǎng)一 在20170115T00：00-05：40時(shí)間內(nèi)的SCADA系統(tǒng)數(shù)據(jù)（采樣頻率為10min，共計(jì)35個(gè)數(shù)據(jù)點(diǎn)）?，F(xiàn)場(chǎng)運(yùn)維記錄顯示，05：00出現(xiàn)機(jī)組葉片結(jié)冰。

圖1 葉片結(jié)冰前后部分SCADA參數(shù)變化特性Fig.1 Variation characteristics of SCADA parameters before and after blade icing

機(jī)組在不同工況下的運(yùn)行參數(shù)值如表2所示。其中，風(fēng)速V、葉片轉(zhuǎn)速n、功率P、槳距角 β可直接從SCADA中獲取。

表2 不同工況下機(jī)組參數(shù)特性Table2 Parameter characteristic of the WT under different working conditions

根據(jù)風(fēng)速值，對(duì) 比表2，在 序列點(diǎn)[1～13]時(shí) 段（圖1）機(jī)組處于工況Ⅰ區(qū)，V，β，風(fēng)能利用系數(shù)Cp值都處于正常狀態(tài)；在序列點(diǎn)[13～29]時(shí)段，風(fēng)速降低，機(jī)組處于Cpmax控制狀態(tài)的工況Ⅱ區(qū)，V隨風(fēng)速減小而降低屬于正常，各項(xiàng)參數(shù)值保持在正常范圍內(nèi)；在序列點(diǎn)[29～33]時(shí)段，風(fēng)速由5.7m/s持續(xù)增大到10m/s，接近額定風(fēng)速，此時(shí)機(jī)組運(yùn)行于Ⅱ區(qū)-Ⅲ區(qū)邊界工況區(qū)間，Cp還維持在Cpmax，β也在0°附近，V理應(yīng)逐漸增大到滿發(fā)狀態(tài)，但從29點(diǎn)開始，Cp自0.43快速下降，V雖然在增大，但輸出功率達(dá)不到（滿發(fā)）狀態(tài)Ⅲ，說明機(jī)組在轉(zhuǎn)換風(fēng)能方面出現(xiàn)異常。即從04：40時(shí)刻起葉片結(jié)冰。

運(yùn) 行 在[Vmin，Vmax]，[nmin，nmax]各 工 況 區(qū) 間 上 的 機(jī)組 氣 動(dòng) 特 性 參 數(shù) λmin與 λmax，Cpmin與Cpmax， 可 由 式（1），（2）計(jì) 算。

式中：λ為葉尖速比；R為葉片半徑。

1.3 葉片結(jié)冰時(shí)的塔架振動(dòng)特性

風(fēng)力機(jī)葉片表面結(jié)冰是隨機(jī)過程，呈非均勻分布，必然導(dǎo)致葉輪質(zhì)量不平衡，機(jī)組塔架產(chǎn)生新的振動(dòng)特性。圖2為風(fēng)電場(chǎng)二某機(jī)組在20170319 T00：00-23：59，葉片結(jié)冰前后塔架振動(dòng)數(shù)據(jù)時(shí)序圖（采樣頻率為1min，共計(jì)1440個(gè)數(shù)據(jù)點(diǎn)）。顯然，葉 片 結(jié) 冰 前（0～450）后（450～1200）直 至 停 機(jī)（1200點(diǎn)之后），塔架振動(dòng)特征明顯不同。

圖2 風(fēng)電機(jī)組塔架時(shí)頻振動(dòng)Fig.2 WT tower time-frequency vibration

1.4 機(jī)械測(cè)風(fēng)儀測(cè)值變化特性

結(jié)冰增大了測(cè)風(fēng)儀轉(zhuǎn)軸的摩擦力，致使其所測(cè)風(fēng)速值偏小，開始出現(xiàn)測(cè)值誤差。相對(duì)于機(jī)械式風(fēng)速儀在低溫環(huán)境中易受雨雪、凍雨等結(jié)冰影響，超聲波風(fēng)速儀則不受低溫環(huán)境影響，有著更高的可靠性和準(zhǔn)確性。因此，同一臺(tái)機(jī)組上的兩臺(tái)測(cè)風(fēng)儀在結(jié)冰前后所測(cè)風(fēng)速值之差就會(huì)發(fā)生變化。設(shè)：

2 葉片結(jié)冰辨識(shí)邏輯

結(jié)冰改變?nèi)~片表面翼型，直接導(dǎo)致其氣動(dòng)性能下降，減小風(fēng)能利用率，影響輸出功率。葉片結(jié)冰呈無規(guī)則形狀，必然引起葉輪不平衡，導(dǎo)致塔架產(chǎn)生新的振動(dòng)特性。結(jié)冰期間的機(jī)械式風(fēng)速儀也因結(jié)冰阻尼增大，測(cè)風(fēng)精度下降。綜合葉片結(jié)冰狀態(tài)特性和結(jié)冰天氣條件，形成葉片結(jié)冰辨識(shí)邏輯如圖3所示，葉片結(jié)冰辨識(shí)流程如圖4所示。

圖3 葉片結(jié)冰與辨識(shí)邏輯Fig.3 Blade icing and identification logic

圖4 葉片結(jié)冰辨識(shí)流程圖Fig.4 Blade icing identification process

3 葉片結(jié)冰數(shù)值特性辨識(shí)

遴選與葉輪轉(zhuǎn)換特性、振動(dòng)特性、輸出功率特性有關(guān)的參數(shù)，將非結(jié)冰機(jī)組歷史數(shù)據(jù)劃分為訓(xùn)練集和測(cè)試集兩部分，運(yùn)用支持向量機(jī)算法建立參數(shù)回歸模型，計(jì)算參數(shù)歷史回歸值。

依據(jù)皮爾遜積矩相關(guān)系數(shù)法的相關(guān)性分析，獲得機(jī)組葉片結(jié)冰狀態(tài)的伴隨參數(shù)。相關(guān)系數(shù)r是描述變量之間線性相關(guān)程度的量。

式中：Xi，Yi為變量樣本；，樣本平均值。

r的 取 值 為[-1，1]。當(dāng){r|-1≤r＜0}時(shí)，負(fù) 相 關(guān)；當(dāng){r|0＜r≤1}時(shí)，正相 關(guān)；當(dāng)r=0時(shí)，不 相 關(guān)。r的數(shù)值越大，表明變量之間相關(guān)性越強(qiáng)。相關(guān)度分布如表3所示。

表3 相關(guān)度劃分Table3 Divide the correlation

對(duì)于給定的輸入-輸出樣本數(shù)據(jù)集 （xi，yi|i=1，2，…，N），其 中，xi為 第i個(gè)d維 輸 入 樣 本 變 量，yi為第i個(gè)輸出標(biāo)量，N為訓(xùn)練樣本個(gè)數(shù)。依照優(yōu)化目標(biāo)函數(shù)：

式 中：w，b為 超 平 面 參 數(shù)；C為 懲 罰 系 數(shù)；ξi，為松弛變量。

對(duì)于約束條件，SVR模型的回歸函數(shù)為

式 中：f（x）為 輸 入 向 量x對(duì) 應(yīng) 的 回 歸 輸 出；αi，為

拉格朗日乘子。

式中：σ為核寬參數(shù)。

考慮到風(fēng)電機(jī)組SCADA數(shù)據(jù)具有多樣本、多特征的特點(diǎn)，選擇徑向基核函數(shù)作為SVR核函數(shù)，構(gòu)造支持向量機(jī)回歸函數(shù)。

式（5）中的C和式（7）中的 σ取值，決定著SVR模型的學(xué)習(xí)能力和預(yù)測(cè)能力。為了得到具有最佳泛化性能的SVR模型，本文采用PSO算法對(duì)模型中的參數(shù)（C，σ）進(jìn)行優(yōu)化，從中選出最佳的C與σ。基于SVR進(jìn)行參數(shù)偏離回歸預(yù)測(cè)如圖5所示。

圖5 基于SVR算法的參數(shù)偏離回歸辨識(shí)Fig.5 Regression Identification of Parameter Deviation Based on SVR Algorithm

利用回歸函數(shù)可計(jì)算出正常狀態(tài)下的參數(shù)預(yù)測(cè)值。參數(shù)值偏離計(jì)算如下：

當(dāng) δ*偏 離 滿 足 條 件：δP＞ηP＆δA＞ηA＆δCp＞ηCp（η*為偏離閾值），則可判定葉片出現(xiàn)結(jié)冰現(xiàn)象。

4 實(shí)測(cè)數(shù)據(jù)驗(yàn)證

4.1 SCADA相關(guān)參數(shù)辨識(shí)

選取寧夏某風(fēng)電場(chǎng)結(jié)冰機(jī)組數(shù)據(jù)驗(yàn)證上述模型的正確性。機(jī)組主要參數(shù)：葉輪直徑為110m，額定功率為2100kW，切入風(fēng)速為3m/s，額定風(fēng)速為10.5m/s，SCADA數(shù)據(jù)采集頻率為1min。

4.1.1確定預(yù)測(cè)輸入?yún)⒘?/p>

隨機(jī)取#A50機(jī)組連續(xù)5d的歷史數(shù)據(jù)作為樣本，依據(jù)皮爾遜積矩相關(guān)系數(shù)法，取｜r｜≥0.6（中、高相關(guān)程度）作為標(biāo)準(zhǔn)選取相關(guān)參數(shù)，計(jì)算各監(jiān)測(cè)參數(shù)與有功功率的相關(guān)系數(shù)，得到變流器無功功率均值、發(fā)電量均值、風(fēng)速1（機(jī)械式測(cè)風(fēng)儀）均值、風(fēng)速2（超聲波式測(cè)風(fēng)儀）均值、30s風(fēng)速均值、5min風(fēng)速均值、輪轂轉(zhuǎn)速均值、發(fā)電機(jī)轉(zhuǎn)速均值、轉(zhuǎn) 矩 反 饋 平 均 值、電 網(wǎng) 電 壓AB（BC，CA）均值、電網(wǎng)A（B，C）相電流均值、電網(wǎng)線電壓L1L2（L2L3，L3L1）均 值、電 網(wǎng) 出 口 線 電 流L1（L2，L3）均 值 和 塔 架 振 動(dòng) 加 速 度 均 值（波 段5，6，7，8，9，11，12）與有功功率輸出量相關(guān)的輸入?yún)⒘俊?/p>

考慮到機(jī)組相關(guān)變量數(shù)據(jù)具有不同單位和量級(jí)，且指標(biāo)之間的量值差異非常大，為保證預(yù)測(cè)結(jié)果的準(zhǔn)確性，需要對(duì)原始數(shù)據(jù)進(jìn)行歸一化[0，1]處理。

4.1.2選取驗(yàn)證數(shù)據(jù)

運(yùn)行記錄顯示，#A50機(jī)組在20170319出現(xiàn)葉片結(jié)冰現(xiàn)象，并在19：56停機(jī)。由于沒有安裝任何監(jiān)測(cè)葉片結(jié)冰的軟、硬件設(shè)備，無法準(zhǔn)確確定風(fēng)力機(jī)葉片開始結(jié)冰的時(shí)間。鑒于此，本文選取機(jī)組SCADA系 統(tǒng) 在20170319T00：00-19：56的 監(jiān) 測(cè) 數(shù)據(jù)作為分析、驗(yàn)證數(shù)據(jù)。

4.1.3輸出功率數(shù)據(jù)辨識(shí)結(jié)冰

圖6（a）為SCADA記錄的有功功率1440個(gè)數(shù)據(jù)點(diǎn)以及應(yīng)用回歸預(yù)測(cè)計(jì)算出的相應(yīng)功率值，圖6（b）為 其 部 分 時(shí) 段 數(shù) 據(jù)。從 圖6（b）可 知，大 約在460點(diǎn)之后，預(yù)測(cè)值開始偏離SCADA記錄值。

圖6 有功功率預(yù)測(cè)值與原始值對(duì)比Fig.6 Compare the predicted value of active power with the original value

圖7為機(jī)組有功功率的相對(duì)誤差曲線。由圖7可知，自458點(diǎn)之后，功率損失持續(xù)上升，即20170319T07：37前后發(fā)生“可能因葉片結(jié)冰”導(dǎo)致的機(jī)組整體性能“異?！?。

圖7 有功功率400～600序列點(diǎn)相對(duì)誤差Fig.7 Active power400～600sequence point relative error

4.1.4Cp數(shù)據(jù)辨識(shí)結(jié)冰

根據(jù)所選SCADA驗(yàn)證數(shù)據(jù)，計(jì)算機(jī)組Cp時(shí)序 值（圖8）。

對(duì)比表3工況劃分與圖8，#A50機(jī)組處于兩種 狀 態(tài)，[1～342]時(shí) 序 點(diǎn) 間 處 于 第 一 狀 態(tài)，[343～1440]時(shí)序點(diǎn)間處于第二狀態(tài)。

[1～342]時(shí) 序 點(diǎn) 間，風(fēng) 速 為10.5m/s，機(jī) 組 處 于Ⅳ工況區(qū)，葉輪轉(zhuǎn)速、Cp，β正常波動(dòng)變化，機(jī)組正常運(yùn)行。

[343～1440]時(shí) 序 點(diǎn) 間，風(fēng) 速 由10.5m/s逐 步降 低 到4.198m/s，之 后[443～1440]時(shí) 序 點(diǎn)，風(fēng) 速 持續(xù)波動(dòng)在4.198～10m/s，正常狀態(tài)下機(jī)組應(yīng)處于Ⅱ工況區(qū)，即Cpmax階段，β不變，發(fā)電機(jī)轉(zhuǎn)速跟隨風(fēng)速而變化。

圖9為[400～600]數(shù)據(jù)時(shí)序點(diǎn)對(duì)應(yīng)的多個(gè)參數(shù)值。

圖9 時(shí)序400～600點(diǎn)間參數(shù)值Fig.9 Parameter values of sequence points in 400～600interval

由圖9可知，[442～600]區(qū)段風(fēng)速從4.198m/s開始持續(xù)增大，機(jī)組理應(yīng)處于Cpmax，但與[400～452]區(qū)段不同的是，圖中自453點(diǎn)開始，Cp出現(xiàn)快速下降，表明機(jī)組于時(shí)序453點(diǎn)即20170319T07：32時(shí)刻起，開始出現(xiàn)氣動(dòng)方面“異?！薄?/p>

4.1.5振動(dòng)特征數(shù)據(jù)辨識(shí)結(jié)冰

將塔架振動(dòng)（振動(dòng)加速度平均值）信號(hào)數(shù)據(jù)應(yīng)用EEMD分解，選取前3個(gè)IMF特征信號(hào)進(jìn)行分析（圖10）。

圖10 EEMD信號(hào)特征分解圖Fig.10 EEMD signal characteristics decomposition diagram

由圖10可知，3種不同頻率信號(hào)在441點(diǎn)附近的幅值都各自發(fā)生突然增大現(xiàn)象，即機(jī)組塔架于20170319T7：20左右，出現(xiàn)塔架振動(dòng)＂異常＂現(xiàn)象。

上述3種模型計(jì)算結(jié)果對(duì)比如表5所示。

表5 3種SCADA參量對(duì)比判定機(jī)組葉片結(jié)冰Table5 Comparison of three SCADA parameters to determine the blade icing of WT

表5中，3種參量計(jì)算結(jié)果表明，#A50機(jī)組大約從450點(diǎn)開始出現(xiàn)運(yùn)行異常，即葉片結(jié)冰。相對(duì)于機(jī)組輸出功率來說，振動(dòng)信號(hào)更加敏感，感知葉片結(jié)冰時(shí)刻更早一些。

4.2 不同類型風(fēng)速測(cè)量儀驗(yàn)證

機(jī)組#A50共安裝了機(jī)械式和超聲波兩種類型風(fēng)速儀，兩種測(cè)速儀實(shí)測(cè)風(fēng)速數(shù)據(jù)如圖11所示。

圖11 兩種風(fēng)速測(cè)量儀的測(cè)量值對(duì)比Fig.11 Comparison of measured values between the two kinds of anemometer

由圖11（a）可知：正常情況下，機(jī)械式測(cè)速儀和超聲波測(cè)速儀所測(cè)得的風(fēng)速值差別不大、趨勢(shì)一致；在序列點(diǎn)[442～527]區(qū)域附近，葉片結(jié)冰之前，機(jī)械式比超聲波式所測(cè)風(fēng)速值偏大，但結(jié)冰之后則相反[圖11（b）]。不考慮測(cè)速儀發(fā)生機(jī)械類故障，則可以推定是受結(jié)冰摩擦，影響到機(jī)械測(cè)量值。這也間接印證了#A50機(jī)組于時(shí)間序列第450點(diǎn)附近出現(xiàn)了結(jié)冰現(xiàn)象。

5 結(jié)論

①與通過觀察機(jī)組功率曲線這種單一檢測(cè)手段相比，本文采用多參數(shù)組合驗(yàn)證方法，可在葉片出現(xiàn)結(jié)冰的30min內(nèi)辨識(shí)出葉片結(jié)冰，便于及時(shí)采取相應(yīng)的運(yùn)維決策。

②對(duì)比分析機(jī)組輸出功率、風(fēng)能利用系數(shù)和塔架振動(dòng)信號(hào)3種參量的異常時(shí)序，機(jī)組塔架振動(dòng)信號(hào)更加靈敏，可以更早反應(yīng)葉片表面出現(xiàn)覆冰現(xiàn)象。

③通過分析葉片結(jié)冰前后SCADA有關(guān)參數(shù)的變化特性而不是通過增加硬件的檢測(cè)方法，避免了因增添設(shè)備而改變?nèi)~片的氣動(dòng)特性、增加維護(hù)維修復(fù)雜度和運(yùn)行成本。