李 浩 鄭祥豪 張宇寧 ,2) 李金偉

*(華北電力大學電站能量傳遞轉化與系統(tǒng)教育部重點實驗室，北京 102206)

?(華北電力大學能源動力與機械工程學院，北京 102206)

**(中國水利水電科學研究院，北京 100048)

為了提高可再生能源發(fā)電的并網規(guī)模，更好地利用可再生能源，能夠用于削峰填谷、減輕電網波動的儲能技術成為當下研究的熱點。抽水蓄能是當前負荷容量最大、技術最成熟的大規(guī)模電能存儲方式，除了能夠調峰填谷、調頻調相外，還具備黑啟動、緊急事故備用等多種功能[1]。水泵水輪機是抽水蓄能電站實現(xiàn)能量轉換的核心設備，其特點是可以根據(jù)實際需求發(fā)揮水泵或水輪機兩種相反的功能。水泵水輪機運行水頭和轉速高、工況頻繁轉換的工作特性，使其內部經常出現(xiàn)不穩(wěn)定的流態(tài)。這些不穩(wěn)定流態(tài)會使機組產生異常壓力脈動和振動等現(xiàn)象，嚴重威脅抽水蓄能電站的運行穩(wěn)定性和安全性[2-5]。其中，在無葉區(qū)發(fā)生的由動靜干涉現(xiàn)象所產生的異常壓力脈動被認為是危害最大的水力不穩(wěn)定現(xiàn)象。因此，通過實時監(jiān)測和準確識別水泵水輪機的無葉區(qū)流態(tài)，防止產生異常壓力脈動現(xiàn)象，對于保證機組穩(wěn)定和安全運行是十分必要的。

由于目前技術手段的限制，原型水泵水輪機無葉區(qū)的流態(tài)無法被直接觀測到，需要通過間接方法進行無葉區(qū)流態(tài)監(jiān)測。無葉區(qū)的流態(tài)通?？梢苑从吃诰哂胁煌卣黝l率分量的壓力脈動信號中，因此可以通過分析無葉區(qū)壓力脈動信號的頻譜信息，從中提取流態(tài)特征并進行流態(tài)監(jiān)測與識別。目前，許多學者針對水輪機領域的信號特征提取和識別方法進行了研究。薛延剛等[6]使用經驗模態(tài)分解（empirical mode decomposition,EMD）分析了水輪機空化信號，并利用希爾伯特變換（Hilbert transform,HT）將分解結果轉換為二維時頻圖且提取了空化信號的特征，提升了空化信號中高頻部分的識別準確率。張飛等[7]使用EMD 分析了機組甩負荷過程中蝸殼進口壓力脈動信號，并計算各模態(tài)分量的相關性系數(shù)作為特征，識別出了信號中的趨勢項和脈動項。Lu 等[8]使用EMD 對混流式水輪機不同負荷下的尾水管壓力脈動信號進行了分析，并計算了各模態(tài)分量的能量指標作為特征向量，實現(xiàn)了尾水管中不同強度渦帶的特征提取。Zhang 等[9]使用集合經驗模態(tài)分解（ensemble empirical mode decomposition,EEMD）和希爾伯特–黃變換（Hilbert–Huang transform,HHT）對水泵水輪機飛逸暫態(tài)過程中的無葉區(qū)壓力脈動信號進行分析，發(fā)現(xiàn)信號主要特征頻率來源于動靜干涉現(xiàn)象。Zheng 等[10]基于EEMD 和排列熵對水泵水輪機發(fā)電工況不同負荷下的水導軸承處的軸心軌跡進行了提取，并闡述了其與機組內部流態(tài)的關聯(lián)性。蔣文君等[11]使用EEMD 和近似熵提取了水電機組轉子的不同故障特征，并在概率神經網絡中成功對上述故障進行識別。王衛(wèi)玉等[12]提出基于自適應噪聲完備經驗模態(tài)分解（complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN）和樣本熵的水電機組故障信息特征提取方法，并將提取到的故障特征向量輸入至經過參數(shù)優(yōu)化后的支持向量機中實現(xiàn)了分類識別。

通過上述研究發(fā)現(xiàn)，EMD[13]及其改進算法（EEMD[14]和CEEMDAN[15]等）是目前水輪機領域主流的信號特征提取方法。它們均能夠自適應地將目標信號分解成包含不同時間尺度的本征模態(tài)函數(shù)（intrinsic mode function,IMF）。然而，這些方法在處理復雜信號時會發(fā)生模態(tài)混疊、產生虛假分量等現(xiàn)象。模態(tài)混疊即相近的特征時間尺度分布在不同的IMF 中，或者一個IMF 中包含了差異極大的特征時間尺度，會嚴重降低分解精度。另外，目前的研究對于特征提取后的模式識別往往采用的是淺層的機器學習算法，遇到復雜信號特征時難以準確識別，迫切需要采用深度學習算法進一步提升識別準確率。針對以上問題，本文提出了一種創(chuàng)新的原型水泵水輪機無葉區(qū)流態(tài)的特征提取和識別方法。該方法能有效避免模態(tài)混疊現(xiàn)象，并實現(xiàn)復雜流態(tài)的智能識別，為水泵水輪機無葉區(qū)流態(tài)的監(jiān)測提供了一種解決方案。

1 研究方法

本節(jié)介紹本文所使用的相關方法的具體原理，以及所提出方法的具體步驟。

1.1 經驗小波變換

經驗小波變換（empirical wavelet transform,EWT）[16]是由Gilles 在2013 年提出的一種信號模態(tài)分解方法，它結合了小波變換的理論完備性和EMD 自適應分解的優(yōu)勢，通過對信號的傅里葉頻譜分割，找出“有意義”的模態(tài)分量，實現(xiàn)模態(tài)分解。EWT 的主要步驟如下[16]。

（1）對信號x(t) 進行傅里葉變換，得到其傅里葉頻譜x?(ω)，并將頻率ω歸一化至 [ 0,π] 范圍內。

（2）利用N+1 個邊界ωi把傅里葉頻譜x?(ω) 劃分為N個連續(xù)區(qū)間Λi，如圖1 所示，其中ω0=0，ωN=π，剩余N -1 個 邊界由 頻譜劃分方法確定[16]。常用的頻譜劃分方法有極大值法[17]、極大極小值法以及尺度空間變換法[18]等。之后，以每個邊界ωi為中心，定義一個寬度為2τi的過渡區(qū)（圖1 中的陰影部分）。

（3）構 造每個連續(xù)區(qū)間Λi上的經驗小波，其定義為每個Λi上的帶通濾波器。根據(jù)Littlewood–Paley 和Meyer 小波結構[19]，定義經驗尺度函數(shù)φn(ω) 和經驗小波ψn(ω) 的表達式。

（4）計算EWT。按照經典的小波變換形式定義經驗小波變換的細節(jié)相關系數(shù)(n,t) 和近似系數(shù)(0,t) 。(n,t) 由經驗小波函數(shù)與信號的內積得到[16]

(0,t)由尺度函數(shù)與信號內積得到[16]

式中，〈,〉 表示內積運算，“ ”表示復共軛，“^”表示傅里葉變換，“∨”表示傅里葉逆變換。

（5）重構信號，并獲得信號各模態(tài)的表達式。由相關系數(shù)(n,t) 和近似系數(shù)(0,t) 組成的重構信號的表達式為[16]

信號x(t) 的經驗模態(tài)可根據(jù)式（3）形式被定義為[16]

根據(jù)上述介紹可知，執(zhí)行EWT 過程中的關鍵步驟是確定頻譜劃分方法。本文選擇頻譜劃分方法為極大極小值法，即首先計算局部極大值，然后將兩個局部極大值之間最小的局部極小值設置為模態(tài)邊界。使用極大極小值法時需要預先確定模態(tài)個數(shù)N。為保證良好的分解效果，N的值需要恰當選取。如果N值過小，則會導致無法把各個模態(tài)完全分開，造成模態(tài)混疊；如果N值過大，則會導致過度分解，產生虛假模態(tài)，降低計算效率。在本文中，首先將N值范圍設置為5～15，經過多次測試后N值最終確定為9。

1.2 散布熵

散布熵（dispersion entropy,DE）[20]是近年來提出的一種較為新穎的用于評估非線性時間序列復雜性和不規(guī)律性的方法，對比經典的樣本熵[21]和排列熵[22]具有計算速度快、受突變信號影響小的優(yōu)點，其具體計算步驟如下[20]。

（1）使用正態(tài)分布函數(shù)將長度為L的非線性時間序列x={x1,x2,x3,...,xL}進行正態(tài)累計分布積分，得到y(tǒng)={y1,y2,y3,...,yL}，其中0＜y ＜1。

（2）選取類別個數(shù)c，按照式（6）對yj進行線性變換，得到[20]

式中，round 為取整函數(shù)，j=1,2,3,...,L。

（3）選擇嵌入維數(shù)m和時間延遲d，按照式（7）構造嵌入向量[20]

式中，i=1,2,3,...,L-(m-1)d。

（5）計算每種散布模式的概率[20]

（6）由香農熵定義，計算散布熵[20]

其歸一化形式為[20]

嵌入維數(shù)m，類別個數(shù)c和時間延遲d是散布熵計算過程中的關鍵參數(shù)。根據(jù)文獻[20]，d一般取1；c和m取值越大，散布熵的結果越可靠，但同時應滿足cm的值小于L。按照以上原則，本文選取嵌入維數(shù)m=4，類別個數(shù)c=6，時間延遲d=1 。

1.3 卷積神經網絡

卷積神經網絡（convolutional neural network,CNN）是深度學習的代表算法之一，其特點是含有卷積計算且具有深層結構。CNN 主要由若干個輸入層、隱含層和輸出層組成，其中隱含層還分為卷積層、池化層和全連接層。本文使用的CNN 包含一個一維輸入層，兩個一維卷積層，一個全局平均池化層，一個全連接層和一個分類輸出層，具體結構如圖2 所示。

圖2 本文所使用CNN 結構Fig.2 CNN structure adopted in the present paper

1.4 EWT–DE–CNN 特征提取和流態(tài)識別方法

基于以上理論，本文提出了一套EWT–DE–CNN 方法，用于原型水泵水輪機無葉區(qū)壓力脈動的特征提取和流態(tài)識別。首先使用經驗小波變換將壓力脈動信號分解成若干分量，然后計算各分量的散布熵并將計算結果與信號所屬流態(tài)組成特征–標簽對，最后將特征–標簽對輸入卷積神經網絡訓練得到了智能流態(tài)識別模型，具體流程如圖3 所示。

圖3 EWT–DE–CNN 方法流程圖Fig.3 Flow diagram of the EWT–DE–CNN method

2 壓力脈動信號采集實驗

本節(jié)簡要介紹了壓力脈動信號采集實驗相關情況，包括機組信息、實驗概況以及各流態(tài)信號特征，同時對涉及到的水泵水輪機相關參數(shù)進行了無量綱化，以方便后續(xù)討論。

2.1 無量綱參數(shù)

各無量綱參數(shù)右上角用上標“′”標識，具體定義如下。

發(fā)電工況下無量綱負荷（P′）定義為

式中，P表示實際負荷，Prated表示額定輸出功率；當水泵水輪機處于空轉工況時，P′=0 。

抽水工況下無量綱輸入功率（）定義為

式中，Ps表示實際輸入功率，Ps-max表示最大輸入功率，負號表示抽水工況。

無量綱頻率（f′）定義為

式中，f和fn分別表示特征頻率和轉輪轉動頻率。

壓力脈動信號時域圖及其快速傅里葉變換后頻域圖的幅值（和Af）分別定義為

式中，At和Af分別表示其時域和頻域的原始幅值，ρ表示液體密度，g表示重力加速度，H表示水頭或者揚程。

2.2 機組參數(shù)及實驗概況

本文選擇了國內某抽水蓄能電站的水泵水輪機組進行壓力脈動信號采集實驗，機組的詳細工作參數(shù)如表1 所示。

表1 機組工作參數(shù)Table 1 Working parameters of the unit

該機組分別在空轉、發(fā)電和抽水工況進行了實驗，在實驗中使用了滿量程精度為±0.5%的MPM480 型壓力變送器測量機組無葉區(qū)中的壓力脈動信號。信號采樣頻率為1000 Hz，信號長度為2000 個數(shù)據(jù)點。

2.3 流態(tài)特征介紹

實驗共采集了4 種水輪機無葉區(qū)壓力脈動信號，分別對應4 種流態(tài)信息。流態(tài)1 通常發(fā)生在發(fā)電工況P′=0 時的空轉啟動階段。流態(tài)2 通常發(fā)生在發(fā)電工況下P′為24.50%～40.83%的低負荷區(qū)。流態(tài)3 發(fā)生在發(fā)電工況下P′為57.17%～98.00%的高負荷區(qū)。流態(tài)4 發(fā)生在=100% 附近抽水工況。圖4～圖7 分別展示了上述4 種流態(tài)的典型壓力脈動信號樣本。

圖4 P ′=0 時的無葉區(qū)壓力脈動信號Fig.4 Pressure fluctuation signal in the vaneless area atP′=0

圖5 P ′=32.67% 時的無葉區(qū)壓力脈動信號Fig.5 Pressure fluctuation signal in the vaneless area at P′=32.67%

圖6 P ′=89.83% 時的無葉區(qū)壓力脈動信號Fig.6 Pressure fluctuation signal in the vaneless area at P′=89.83%

圖7 =–97.91%時的無葉區(qū)壓力脈動信號Fig.7 Pressure fluctuation signal in the vaneless area at=–97.91%

3 結果與討論

本節(jié)將展示使用本文提出的EWT–DE–CNN 方法分析實測壓力脈動信號結果，并與傳統(tǒng)方法進行對比，討論該方法的有效性和優(yōu)越性。首先，以具有典型壓力脈動信號特征的圖5 中P′=32.67%的壓力脈動信號為例，通過圖片展示了EWT 分解結果。然后，計算了由EWT 分解得到的分量信號的散布熵結果，構建流態(tài)特征向量。最后，利用特征–標簽對訓練了智能流態(tài)識別模型，并進行了流態(tài)識別測試，驗證了其識別的準確性。

3.1 EWT 分析結果

圖8 和 圖9 展示了P′=32.67% 時壓力脈動信號經EWT 分解得到的模態(tài)分量的時域圖和瞬時頻率分布圖。從圖8 中可以看出，壓力脈動信號被分解成9 個分量，且隨著模態(tài)階數(shù)的增加，分量信號的頻率逐漸增高。從圖9 可以看出，9個模態(tài)分量的頻率均勻分布在整個頻率區(qū)間。單個分量的頻率隨時間的變化幅度較小，總體成一條直線，且未出現(xiàn)相交或重疊。圖8 和圖9 的結果說明，EWT 方法較好地按照不同模態(tài)頻率分解了原始信號，沒有出現(xiàn)分解不充分或過度分解的情況。

圖8 P ′=32.67% 時基于EWT 方法的壓力脈動信號模態(tài)分量時域圖Fig.8 Time–domain diagrams of the mode components of pressure fluctuation signal based on EWT method atP′=32.67%

圖9 P ′=32.67% 時基于EWT 方法的壓力脈動信號模態(tài)分量瞬時頻率Fig.9 Instantaneous frequencies of the mode components of pressure fluctuation signal based on EWT method at P′=32.67%

圖10 和圖11 展示了P′=32.67% 時壓力脈動信號經EMD 分解得到的IMF的時域圖和瞬時頻率。從圖10 和圖11 中可以看出，壓力脈動信號被分解成了7 個IMF和一個余項，且頻率由高到低變化。具體觀察各分量可以發(fā)現(xiàn)，高頻分量IMF1在圖10 中出現(xiàn)了兩次明顯的低頻脈動，而在對應的圖11 中IMF1在1 s 和1.8 s 附近出現(xiàn)了兩次較大幅度的頻率驟減。低頻分量IMF4，IMF5和IMF6在圖10 中的波形變化具有一定的相似性，而在對應的圖11 中，IMF4，IMF5和IMF6在低頻區(qū)重疊在一起。圖10 和圖11 的結果表明，EMD方法在分解壓力脈動信號時在高頻分量IMF1中存在部分低頻未被完全分離出去，而低頻分量則被過度分解成了多個分量，說明發(fā)生了模態(tài)混疊現(xiàn)象。

圖10 P ′=32.67% 時基于EMD 方法的壓力脈動信號模態(tài)分量時域圖Fig.10 Time–domain diagrams of the mode components of pressure fluctuation signal based on EMD method atP′=32.67%

圖11 P ′=32.67% 時基于EMD 方法的壓力脈動信號模態(tài)分量瞬時頻率Fig.11 Instantaneous frequencies of the mode components of pressure fluctuation signal based on EMD method at P′=32.67%

3.2 壓力脈動信號特征提取結果

通過式（6）～式（11）步驟計算上述4 種流態(tài)的典型壓力脈動信號，經EWT 或EMD 分解后得到的若干模態(tài)分量的散布熵合并為特征向量，結果如表2 和表3 所示。從表2 中可以看出，4種流態(tài)的特征向量的維度相等，等于EWT 方法的最大模態(tài)個數(shù)N。而在表3 中，特征向量的維數(shù)并不一致，這是由于EMD 方法的特性導致的。EMD 方法自適應地分解信號，產生的模態(tài)數(shù)量不確定，從而導致特征向量維度不統(tǒng)一。特征向量維度不統(tǒng)一對后續(xù)根據(jù)特征向量訓練智能流態(tài)識別模型造成了極大的困難，因此需要對特征向量的維度統(tǒng)一化。經統(tǒng)計，基于EMD 結果的特征向量的維度范圍為7～10，第7 維之后的數(shù)據(jù)將被舍棄。

表2 基于EWT 方法的不同流態(tài)壓力脈動信號特征向量Table 2 Feature vectors of pressure fluctuation signals of different flow statuses based on EWT method

表3 基于EMD 方法的不同流態(tài)壓力脈動信號特征向量Table 3 Feature vectors of pressure fluctuation signals of different flow statuses based on EMD method

3.3 流態(tài)識別結果

將提取到的特征向量與其所屬的類型標簽組成特征–標簽對，輸入神經網絡。神經網絡根據(jù)輸入的特征預測其類型，并與特征對應的類型標簽對比，計算預測準確率。然后根據(jù)梯度下降算法，神經網絡不斷修改內部的權重系數(shù)，提高預測的準確率。此過程稱為神經網絡的訓練過程。本文使用的CNN 訓練參數(shù)如表4 所示。

表4 CNN 訓練參數(shù)設置Table 4 Setting of training parameters of CNN

將特征向量的流態(tài)類型作為標簽，通過特征提取共得到了900 個特征–標簽對，其中流態(tài)1～4 的特征向量–標簽對數(shù)量分別為30，270，540和60。表5 展示了圖4～圖7 中的壓力脈動信號經EWT–DE 提取到的特征向量與其流態(tài)組成的特征–標簽對。將所有特征–標簽對按照8∶2 的比例隨機化分為訓練集與測試集，其中訓練集用于訓練CNN，測試集用于測試智能流態(tài)識別模型的效果。在訓練集和測試集中，每種流態(tài)的樣本都占一定比例，避免只覆蓋一種或兩種流態(tài)的情況。

表5 基于EWT–DE 方法的不同流態(tài)壓力脈動信號的特征–標簽對Table 5 Feature-label pairs of pressure fluctuation signals of different flow statuses based on EWT–DE method

表6 展示基于EWT–DE–CNN 和EMD–DE–CNN 兩種方法訓練出的流態(tài)識別模型的5 次測試結果，其中基于EWT–DE–CNN 方法的模型5次測試的平均準確率為94.89%，達到了較高水平，并顯著高于EMD–DE–CNN 方法的86.45%，表明該模型對無葉區(qū)流態(tài)的識別效果較為準確。圖12 展示了EWT–DE–CNN 方法流態(tài)識別模型準確率最高的第2 次測試（96.11%）的具體結果，模型對4 種流態(tài)識別準確率分別為80%，96.43%，97.27%和88.89%。

表6 基于EWT–DE–CNN 和EMD–DE–CNN 方法得到的智能流態(tài)識別模型測試的準確率(%)Table 6 Accuracy rate of intelligent identification model of flow status based on EWT–DE–CNN and EMD–DE–CNN methods(%)

圖12 基于EWT–DE–CNN 方法的智能流態(tài)識別模型第2 次測試結果Fig.12 Result of the second test of intelligent identification model of flow status based on EWT–DE–CNN method

4 結論

水泵水輪機由于無葉區(qū)中的異常流動引起的壓力脈動嚴重威脅了機組的安全平穩(wěn)運行。監(jiān)測和識別不同運行工況下無葉空間內的流動狀態(tài)，是保證整個抽水蓄能電站安全穩(wěn)定運行的必要條件。本文基于EWT，DE 和CNN，提出了一種原型水泵水輪機無葉區(qū)流態(tài)的特征提取和識別方法，并使用實測壓力脈動信號進行了驗證，結論如下。

（1）使用EWT 方法對水泵水輪機無葉區(qū)壓力脈動信號進行模態(tài)分解，有效避免了EMD 方法中出現(xiàn)的模態(tài)混疊現(xiàn)象，改善了分解效果。

（2）使用DE 進行特征提取，特征向量的維度由模態(tài)分解結果數(shù)量決定。EWT 方法產生固定個數(shù)的信號分量，解決了EMD 方法中特征向量維度不一致的問題。

（3）本文提出EWT–DE–CNN 方法得到的智能流態(tài)識別模型，有效實現(xiàn)了發(fā)電、抽水、空轉工況下水泵水輪機無葉區(qū)壓力脈動信號的特征提取和識別，識別的平均準確率為94.89%。