朱莎莎，楊 光，林雯雯

(1.華北水利水電大學(xué)水利學(xué)院，河南 鄭州 450046；2.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

1 概述

20世紀(jì)末開始，可持續(xù)發(fā)展一直是我們國家堅(jiān)守的戰(zhàn)略。隨著經(jīng)濟(jì)、社會(huì)持續(xù)穩(wěn)定且快速發(fā)展，對(duì)能源的需求量愈來愈大。但我國能源結(jié)構(gòu)復(fù)雜，能源利用占比主要以煤炭為主，水能資源位居第二，清潔能源占比仍較少。由于化石能源對(duì)于環(huán)境的污染巨大，導(dǎo)致溫室效應(yīng)加劇，不利于可持續(xù)發(fā)展，而水能具有清潔高效、在電網(wǎng)中能調(diào)峰、調(diào)頻、填谷等優(yōu)點(diǎn)，國家能源局在“十四五”規(guī)劃提出，要加快抽水蓄能電站建設(shè)和新型儲(chǔ)能技術(shù)規(guī)?；瘧?yīng)用。隨著電子科技和城鎮(zhèn)化的發(fā)展，電力系統(tǒng)的高峰低谷差值愈來愈大。而電力作為一種特殊的、不可儲(chǔ)存的產(chǎn)品，面臨著十分嚴(yán)峻的考驗(yàn)[1]。抽水蓄能電站將不可儲(chǔ)存的電能轉(zhuǎn)化為水流的勢(shì)能儲(chǔ)存起來，是世界上最重要的大規(guī)模儲(chǔ)能手段之一，與傳統(tǒng)電站相比，更能適應(yīng)頻繁改變需求量的電力市場(chǎng)。但水輪機(jī)頻繁改變運(yùn)行工況以適應(yīng)電力市場(chǎng)日益增長的需求，在偏離設(shè)計(jì)工況條件下，會(huì)歷經(jīng)壓力脈動(dòng)幅值較高的區(qū)域運(yùn)行[2]。抽水蓄能電站水泵水輪機(jī)甩負(fù)荷過程中產(chǎn)生的壓力脈動(dòng)信號(hào)組成成分復(fù)雜，而信號(hào)采集過程中又會(huì)再次受到各種干擾信號(hào)的影響，導(dǎo)致采集到的信號(hào)難以直觀且準(zhǔn)確地表明有用的信息。對(duì)水泵水輪機(jī)的壓力脈動(dòng)信號(hào)進(jìn)行數(shù)字信號(hào)處理，剝離出其中干擾的噪音或者不相關(guān)信息，提取表征水輪機(jī)狀態(tài)的特征參量進(jìn)行信號(hào)分析，可實(shí)時(shí)監(jiān)測(cè)壓力脈動(dòng)信號(hào)；且目前設(shè)備維修的方式已經(jīng)逐步轉(zhuǎn)變?yōu)闋顟B(tài)檢修，而不是過去的定期維修[3]，準(zhǔn)確判斷水輪機(jī)運(yùn)行工況，可以減少事故的發(fā)生，還能提高經(jīng)濟(jì)性和檢修效率。

開始國內(nèi)外學(xué)者運(yùn)用傅里葉分析處理信號(hào)，但學(xué)者們發(fā)現(xiàn)傅里葉分析適用范圍有限，后續(xù)又提出了短時(shí)傅里葉變換、小波變換等方法，雖然小波變換能處理非線性信號(hào)，但一直無法脫離小波基的束縛。而水輪機(jī)水道系統(tǒng)中的壓力脈動(dòng)產(chǎn)生原因復(fù)雜，為了準(zhǔn)確地從壓力脈動(dòng)信號(hào)中提取出有用的信號(hào)，并及時(shí)判斷機(jī)組運(yùn)行狀態(tài)，國內(nèi)外展開了許多研究。目前較為常用的非線性非平穩(wěn)信號(hào)處理方法包括：經(jīng)驗(yàn)?zāi)B(tài)分解、變分模態(tài)分解、自適應(yīng)迭代濾波算法等[4]。1998年，NordenE.Huang等人提出了Hilbert-Huang變換(HHT)[5]，這是第一種真正意義上可以處理非線性非平穩(wěn)信號(hào)的信號(hào)處理方法。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)為HHT變換第一步，其時(shí)頻分辨率較高且自適應(yīng)性較好，但可能會(huì)出現(xiàn)模式混疊、端點(diǎn)效應(yīng)和停止條件等問題[6]，適合于分析低頻非平穩(wěn)的水輪機(jī)水力壓力脈動(dòng)信號(hào)。楊世錫[7]等在2004年提出基于高次樣條插值的EMD新算法。2008年NordenE.Huang等人再次提出EEMD方法，在整個(gè)時(shí)頻空間中加入白噪音，再進(jìn)行足夠多的獨(dú)立測(cè)試，促進(jìn)抗混分解，唯一穩(wěn)定的即是所需的信號(hào)[8]。2014年Dragomiretskiy等[6]提出變分模態(tài)分解(VMD)，與EMD相比具有更好的自適應(yīng)能力，能夠更加精確有效地揭示整個(gè)機(jī)組振動(dòng)的物理學(xué)特性。2013年Gilles等提出經(jīng)驗(yàn)小波變換(EWT)代替EMD分解。2017年于曉東等[9]運(yùn)用VMD-PE的方法去除壓力脈動(dòng)信號(hào)中的噪聲，與EEMD相比，去噪效果良好，但自適應(yīng)性較差，且參數(shù)K需提前確定。2018年Xiao-li WANG等[10]基于低噪聲會(huì)對(duì)HHT結(jié)果造成很大干擾的問題，提出運(yùn)用HHT變換和中位值濾波法聯(lián)合對(duì)水輪機(jī)尾水管脈動(dòng)信號(hào)進(jìn)行分析，結(jié)果表明，該方法能有效揭示信號(hào)的時(shí)頻特征。2020年趙志爐等[11]提出了ALIF-PE算法，該方法采用自適應(yīng)迭代濾波算法(ALIF)對(duì)信號(hào)進(jìn)行分解，與VMD-PE相比的抗混疊表現(xiàn)更好。

目前的研究大多數(shù)都是對(duì)抽水蓄能電站壓力脈動(dòng)數(shù)據(jù)運(yùn)用濾波算法進(jìn)行處理，研究內(nèi)容大多為理論研究，有關(guān)工程實(shí)際應(yīng)用方面的內(nèi)容較少。因此，本論文運(yùn)用VB.NET語言編寫處理水道系統(tǒng)中的壓力脈動(dòng)實(shí)測(cè)數(shù)據(jù)的濾波軟件，使抽水蓄能電站甩負(fù)荷數(shù)據(jù)處理更加方便快捷，具有較高的理論價(jià)值和研究意義。

2 濾波軟件程序編程

濾波器可以分為模擬濾波器和數(shù)字濾波器。模擬濾波器可以直接處理時(shí)間連續(xù)、幅值也連續(xù)的模擬信號(hào)。但模擬濾波器受限于元件，更適合處理高頻信號(hào)，對(duì)于低通濾波要求硬件有較大的電容或電感，經(jīng)濟(jì)開銷大。而數(shù)字濾波器通過計(jì)算機(jī)A/D轉(zhuǎn)換器實(shí)現(xiàn)模擬信號(hào)與數(shù)字信號(hào)的轉(zhuǎn)換，運(yùn)用數(shù)字信號(hào)處理器對(duì)轉(zhuǎn)換后的數(shù)字信號(hào)進(jìn)行處理，再輸出相對(duì)應(yīng)的模擬信號(hào)，基本運(yùn)行過程如圖1所示。不同于模擬濾波器，數(shù)字濾波器增加功能就是增加程序，不受元件誤差的影響，經(jīng)濟(jì)性好，繼承了模擬濾波器的優(yōu)點(diǎn)，可以同時(shí)處理高頻和低頻信號(hào)。

圖1 數(shù)字濾波器運(yùn)行流程圖

本文利用Visual Studio 2017編寫濾波軟件，該軟件采用可視化(Visual)程序設(shè)計(jì)，它能讓程序設(shè)計(jì)人員通過拼圖的方式來構(gòu)建各種界面，不再需要為設(shè)計(jì)界面而編寫大量代碼，節(jié)省了程序設(shè)計(jì)人員的時(shí)間與精力。

2.1 程序編制

濾波算法的濾波過程代碼偏長，且需要多次調(diào)用，因此定義函數(shù)過程，即Function過程，就可以在濾波軟件中調(diào)用這個(gè)函數(shù)過程了。調(diào)用時(shí)通常需要將實(shí)際參數(shù)傳遞給形式參數(shù)，函數(shù)過程利用這些參數(shù)進(jìn)行計(jì)算，然后將結(jié)果返回。

2.1.1中位值濾波法

中位值濾波法即中值濾波法，是一種經(jīng)典的非線性信號(hào)處理技術(shù)，能有效濾除偶發(fā)性非周期性波動(dòng)引發(fā)的脈沖信號(hào)，常用于濾除溫度、液位等慢變量信號(hào)，但不適宜處理流量、壓力等快變量[12]。中位值濾波法的Function過程中的形式參數(shù)為字符串File、十進(jìn)制數(shù)val、二維數(shù)組value2(，)，其中File代表導(dǎo)入的數(shù)據(jù)文件；val代表連續(xù)采樣次數(shù)；value2(，)代表存儲(chǔ)濾波后數(shù)據(jù)的二維數(shù)組，以子傳遞。中位值濾波法的Function過程流程圖如圖2所示。當(dāng)連續(xù)采樣次數(shù)val值取n時(shí)，代碼的運(yùn)行原理分為7步：①判斷連續(xù)采樣次數(shù)值是否為奇數(shù)，若非奇數(shù)，則在原來基礎(chǔ)上加1，val為n=n+1，此時(shí)假定最開始的n為奇數(shù)。②i作為計(jì)步器，每讀取一行數(shù)據(jù)時(shí)，i=i+1，當(dāng)i=n時(shí)，表明已經(jīng)讀取了n行數(shù)據(jù)了，可以進(jìn)行濾波。③采用冒泡排序法，將n個(gè)時(shí)間點(diǎn)對(duì)應(yīng)的值從小到大排列。④選取最中間的值，將n個(gè)時(shí)間點(diǎn)的中值，即第(n+1)/2個(gè)時(shí)間點(diǎn)與之配對(duì)，形成一組數(shù)據(jù)，并賦值給value2。⑤再次讀取第i+1行到第i+n行的數(shù)據(jù)賦值給數(shù)組time和value。⑥重新回到步驟3，開始循環(huán)。⑦一直到全部數(shù)據(jù)錄入，并處理完畢。

圖2 中位值濾波法流程圖

2.1.2算數(shù)值平均濾波法

算術(shù)平均值濾波法，也稱算術(shù)平均濾波法，適用于對(duì)存在周期性干擾進(jìn)行濾波。這種信號(hào)的特點(diǎn)是信號(hào)本身在某一數(shù)值范圍附近上下波動(dòng)，如測(cè)量流量、液位時(shí)，但該方法不適用于受到嚴(yán)重脈沖信號(hào)干擾的場(chǎng)合。算數(shù)值平均濾波法的Function過程中的形式參數(shù)其定義與中位值濾波法相同，代碼的運(yùn)行原理不再贅述。其過程流程圖如圖3所示。

圖3 算術(shù)平均值濾波法流程圖

2.1.3遞推平均濾波法

遞推平均濾波法是一種低通濾波器。其原理與算數(shù)平均值濾波法類似，與其不同的是，遞推平均濾波法，每得到一個(gè)運(yùn)算結(jié)果時(shí)間更短，可以應(yīng)用于要求較高的實(shí)時(shí)系統(tǒng)中。遞推平均濾波法的Function過程中的形式參數(shù)其定義與中位值濾波法、算數(shù)平均濾波法相同。代碼的運(yùn)行原理不再贅述。其過程流程圖如圖4所示。

圖4 遞推平均濾波法流程圖

2.1.4經(jīng)驗(yàn)?zāi)B(tài)分解

EMD方法本質(zhì)是通過數(shù)據(jù)特征時(shí)間尺度來獲得固有波動(dòng)模式，然后分解數(shù)據(jù)。其Function過程命名為E_M_D，形式參數(shù)為單精度數(shù)組xy(，)，單精度數(shù)組per_error(，)，整型max_iterations，整型base_extremums，單精度數(shù)組I_M_F。其中xy(，)代表導(dǎo)入的全部濾波數(shù)據(jù)，以值傳遞；per_error(，)代表允許誤差，以值傳遞；max_iterations代表最大迭代次數(shù)，以值傳遞；base_extremums代表基準(zhǔn)極值個(gè)數(shù)，以值傳遞；I_M_F(，)用來存放IMF函數(shù)和最后的殘值，以子傳遞。其過程流程圖如圖5所示。EMD方法Function過程原理如下：①將輸入的全部濾波數(shù)據(jù)賦值給二維數(shù)組h(，)和re_value(，)，以便于接下來IMF和殘值的運(yùn)算；②K=1，求出數(shù)據(jù)的極大值序列和極小值序列，并分別賦值給二維數(shù)組xy1(，)和xy2(，)；③利用xy1(，)和xy2(，)求出上包絡(luò)線序列和下包絡(luò)線序列，并賦值給up_envelopes(，)和down_envelopes(，)；④利用up_envelopes(，)和down_envelopes(，)求出上下包絡(luò)線均值序列并賦值給aver_envelopes(，);使h(，)減去aver_envelopes(，)，并賦值給h(，)，使之臨時(shí)存放原數(shù)據(jù)減去IMF；⑤判斷h(，)是否滿足約束條件，若滿足繼續(xù)下一條，若不滿足跳轉(zhuǎn)到步驟2，K值加1；⑥將h(，)賦值給存放IMF函數(shù)的數(shù)組I_M_F(，);⑦將re_value(，)減去I_M_F(，)，并賦值給re_value(，)，并判斷殘值是否為單調(diào)函數(shù)，若不為單調(diào)函數(shù)則繼續(xù)迭代，直到最大迭代次數(shù)或者殘值稱為單調(diào)函數(shù)；⑧若殘值為單調(diào)函數(shù)，則將此賦值給I_M_F(，)，至此，EMD方法完成。

圖5 EMD方法流程圖

2.2 代碼運(yùn)行過程

(1)導(dǎo)入文件。單擊“瀏覽”文本框，則該文本框的click事件啟動(dòng)，清除濾波前后波形數(shù)據(jù)文件，運(yùn)用StreamReader讀取原數(shù)據(jù)文件。

(2)圖像顯示。該過程引用次數(shù)多，故而定義一個(gè)Function過程，命名為draw()。在導(dǎo)入文件、濾波后都會(huì)將圖像顯示在Chart1中。濾波前數(shù)據(jù)波形為Series(0)，濾波后數(shù)據(jù)波形為Series(1)，數(shù)據(jù)的時(shí)間由x軸表示，壓力值由y軸表示。

(3)數(shù)據(jù)處理。選取“濾波算法”單選框，其對(duì)應(yīng)的控制參數(shù)“允許最大偏差值”的文本框顯現(xiàn)，其余不相干的參數(shù)的Visible變?yōu)镕alse，其界面顯示如圖6所示。

圖6 限幅濾波法界面

(4)單擊“導(dǎo)出”文本框，輸出濾波后數(shù)據(jù)文件，新建一個(gè)csv文件，在數(shù)據(jù)處理的Function過程中，添加一段append代碼，在不清空已錄入數(shù)據(jù)的基礎(chǔ)上，錄入下一行數(shù)據(jù)。

3 數(shù)據(jù)分析與分析

3.1 試驗(yàn)數(shù)據(jù)分析

某抽水蓄能電站采用混流式水泵水輪機(jī)，其壓力脈動(dòng)主要來源見表1[13]。

表1 水泵水輪機(jī)壓力脈動(dòng)主要來源表

(1)運(yùn)用中位值濾波法進(jìn)行分析時(shí)，連續(xù)采樣次數(shù)分別選取5、50、100次，結(jié)果如圖7所示。由濾波后的波形可以看出，在連續(xù)采樣次數(shù)達(dá)到50次后，大約在1405s處出現(xiàn)的尖峰脈沖干擾信號(hào)被過濾，曲線隨著連續(xù)采樣次數(shù)的增多，平滑度越高，但靈敏度越低。例如，連續(xù)采樣次數(shù)為5和50次的圖中在大約580s處出現(xiàn)的局部極小值，在連續(xù)采樣次數(shù)為100次的圖中被過濾了，但也有可能與波形顯示框精度較低有關(guān)。

圖7 中位值濾波法

(2)運(yùn)用算數(shù)平均濾波法進(jìn)行分析時(shí)，連續(xù)采樣次數(shù)分別選取5、50、100次，結(jié)果如圖8所示。由濾波后波形可以看出，曲線隨著連續(xù)采樣次數(shù)增多平滑度越高；大約在1405s處出現(xiàn)的尖峰信號(hào)幅值越來越小，直到連續(xù)采樣次數(shù)達(dá)到100次時(shí)，尖峰信號(hào)完全被過濾；大約在580s處出現(xiàn)的局部極小值并沒有完全消失，說明濾波精度較好。

圖8 算數(shù)平均濾波法

(3)運(yùn)用遞推平均濾波法進(jìn)行分析時(shí)，連續(xù)采樣次數(shù)分別選取5、50、200次，結(jié)果如圖9所示。由濾波后的波形可以看出，曲線平滑度隨著連續(xù)采樣次數(shù)的增多而愈來愈高，連續(xù)采樣次數(shù)達(dá)到200次時(shí)，鋸齒狀波形基本消失，260～470s和1440～1540s鋸齒狀波形較明顯；大約在1405s處出現(xiàn)的尖峰信號(hào)幅值越來越小，直到連續(xù)采樣次數(shù)達(dá)到200次時(shí)，尖峰信號(hào)仍然沒有被完全被過濾，還存在一定幅值；大約第505s處和第605s處出現(xiàn)的局部極小值并沒有隨著采樣次數(shù)的增多而消失。

圖9 遞推平均濾波法

(4)運(yùn)用EMD方法進(jìn)行分析時(shí)，先對(duì)仿真信號(hào)進(jìn)行EMD分解。仿真信號(hào)由一個(gè)幅值為10Pa、頻率為2Hz信號(hào)和一個(gè)幅值為5Pa、頻率為10Hz信號(hào)組成，闕值取0.001，IMF允許最大迭代次數(shù)取5次，殘值允許極值點(diǎn)個(gè)數(shù)取10個(gè)，極值延拓基準(zhǔn)極值個(gè)數(shù)取5個(gè)，分解出的IMF如圖10所示，可以分析出；IMF最大階數(shù)為2，EMD方法分解出IMF1為10Hz的信號(hào)，IMF2為2Hz的信號(hào)，排列順序是由頻率從高到低排列。但是在第4s處信號(hào)分解，存在一定問題，這是程序需要改進(jìn)的地方。

圖10 仿真信號(hào)EMD分解

對(duì)實(shí)測(cè)信號(hào)進(jìn)行分析，闕值取0.001，IMF允許最大迭代次數(shù)取5次，殘值允許極值點(diǎn)個(gè)數(shù)取10個(gè)，極值延拓基準(zhǔn)極值個(gè)數(shù)取5個(gè)。濾波前波形與分解結(jié)果如圖11所示。

圖11 實(shí)測(cè)信號(hào)EMD分解

由分解出來的IMF可知，最大IMF階數(shù)為8，只顯示具有一定物理意義的單分量信號(hào)IMF1-IMF4，剩下的為殘值，不具任何意義。該分解過程只是將復(fù)雜的尾水管進(jìn)口處壓力脈動(dòng)信號(hào)在一定特征尺度下分離，得到的單分量信號(hào)只具有一定的物理意義，但無法得出能夠體現(xiàn)水輪機(jī)壓力脈動(dòng)信號(hào)的特征[14]，需要進(jìn)行傅里葉變換或者希爾伯特變換。該程序目前無法提取出各IMF分量，因此，仍然需要進(jìn)一步改進(jìn)。

3.2 單機(jī)甩負(fù)荷試驗(yàn)濾波結(jié)果分析

工程實(shí)際中為了機(jī)組的安全，施工人員比較關(guān)注尾水管進(jìn)口最大/小壓力、以及最大轉(zhuǎn)速的上升率等。即需要對(duì)抽水蓄能電站的濾波結(jié)果進(jìn)行修正，確保工程的安全性。

仙居抽水蓄能電站位于浙東南中心地帶仙居縣境內(nèi)，電站設(shè)計(jì)安裝4臺(tái)單機(jī)容量375MW的可逆式抽水蓄能機(jī)組，總裝機(jī)容量1500MW，年發(fā)電量25.125億度，年抽水耗電量34.5億度。在連續(xù)采樣次數(shù)為100次的情況下，運(yùn)用中位值濾波法、算數(shù)平均濾波法以及遞推平均濾波法對(duì)機(jī)組實(shí)測(cè)數(shù)據(jù)進(jìn)行濾波，根據(jù)機(jī)組單機(jī)甩負(fù)荷試驗(yàn)情況，進(jìn)行實(shí)測(cè)與濾波后結(jié)果分析，根據(jù)實(shí)測(cè)數(shù)據(jù)與濾波后數(shù)據(jù)進(jìn)行對(duì)比得到最值差值和壓力脈動(dòng)修正值，并比較3種濾波法的結(jié)果。

3.2.1尾水管壓力對(duì)比

不同工況時(shí)尾水管進(jìn)口壓力濾波后波形對(duì)比、甩負(fù)荷試驗(yàn)水位組合表及尾水管進(jìn)口最大/小壓力數(shù)據(jù)如圖12—13及表2—4所示，其中，計(jì)算誤差和壓力脈動(dòng)修正值采用“(濾波后數(shù)據(jù)-試驗(yàn)實(shí)測(cè)值)/甩前凈水頭”確定。

圖12 甩100%負(fù)荷工況下尾水管進(jìn)口壓力的濾波前后圖形

圖13 甩50%負(fù)荷工況下尾水管進(jìn)口壓力的濾波前后圖形

表2 單機(jī)甩負(fù)荷試驗(yàn)工況水位組合表

表3 尾水管進(jìn)口最大/小壓力實(shí)測(cè)值與濾波后極值對(duì)比 單位：m

表4 尾水管進(jìn)口最大/小壓力相對(duì)差值表 單位：%

3.2.2數(shù)據(jù)分析

為確保工程的安全性，需要對(duì)仙居電站可能發(fā)生的極端控制性工況進(jìn)行復(fù)核計(jì)算，并對(duì)計(jì)算結(jié)果進(jìn)行修正，修正按照已有的試驗(yàn)工況按照最不利情況進(jìn)行修正。修正后特征參量取值見表5—6。

表5 各個(gè)特征參數(shù)修正原則 單位：%

表6 尾水管進(jìn)口最大/小壓力修正后特征參數(shù)表 單位：m

通過對(duì)單機(jī)甩負(fù)荷試驗(yàn)與濾波后數(shù)據(jù)對(duì)比分析可知：單機(jī)100%、50%甩負(fù)荷工況尾水管進(jìn)口最大壓力的試驗(yàn)和不同濾波法濾波后數(shù)據(jù)結(jié)果相對(duì)差值絕對(duì)值對(duì)比可知中位值濾波法>遞推平均濾波法>算術(shù)平均值濾波法，最小壓力相對(duì)差值對(duì)比結(jié)果可知算術(shù)平均值濾波法>中位值濾波法>遞推平均濾波法。

4 結(jié)語

本文運(yùn)用VB.NET語言編寫濾波軟件，該軟件將傳統(tǒng)濾波和EMD方法結(jié)合在一起，可以實(shí)現(xiàn)濾波前后波形對(duì)比以及提取濾波后信號(hào)等功能。并且不需要用戶具備深厚的數(shù)學(xué)知識(shí)，只需要簡單輸入數(shù)據(jù)和參數(shù)，即可得到想要的濾波結(jié)果。針對(duì)過濾影響工程安全的最值的問題，發(fā)現(xiàn)中位值濾波法對(duì)最大壓力的過濾程度最好，算數(shù)平均值對(duì)最小壓力的過濾程度最好。進(jìn)行濾波處理時(shí)，傳統(tǒng)濾波法能較好地過濾尖峰脈沖信號(hào)和噪聲，并對(duì)極值點(diǎn)有較好地保留；但隨著控制參量愈來愈大，濾波后波形愈來愈平滑，但靈敏度越來越低。經(jīng)驗(yàn)?zāi)B(tài)分解可以提取出具有物理意義的IMF分量，但需要進(jìn)一步轉(zhuǎn)換才有意義，且有運(yùn)行速度慢、數(shù)據(jù)精度低等缺陷。后續(xù)可以將EMD方法的算法進(jìn)一步優(yōu)化改進(jìn)。