電力系統(tǒng)強噪聲測量數(shù)據(jù)快速魯棒回歸平滑降噪

靳宗帥，張恒旭

（電網(wǎng)智能化調(diào)度與控制教育部重點實驗室（山東大學(xué)）,山東省濟南市 250061）

0 引言

隨著大規(guī)模風(fēng)電/光伏等分布式發(fā)電、柔性直流、電動汽車/儲能等大功率互動性多元電力電子設(shè)備的接入［1-2］,電網(wǎng)呈現(xiàn)有源化、電力電子化,信號形態(tài)、故障形態(tài)發(fā)生顯著變化,動態(tài)行為更加復(fù)雜［3］。以同步相量測量單元（synchrophasor measurement unit,PMU）為終端的廣域測量系統(tǒng)（wide area measurement system,WAMS）［4］目前主要覆蓋高壓輸電網(wǎng)；隨著配電網(wǎng)測量技術(shù)逐漸成熟［5-9］,其提供的中低壓電網(wǎng)頻率、幅值、相位測量數(shù)據(jù)既蘊含高壓主網(wǎng)響應(yīng)信息,又保留了局部響應(yīng)信息,作為WAMS 的補充數(shù)據(jù)源,實現(xiàn)了電力系統(tǒng)多電壓等級全景式同步監(jiān)測。但測量點與負荷的電氣距離較近,測量數(shù)據(jù)受噪聲污染嚴重,如何快速降噪是亟待解決的問題。

受噪聲污染的信號經(jīng)互感器接入測量裝置,然后由模數(shù)轉(zhuǎn)換模塊轉(zhuǎn)換為離散采樣序列,其所包含的噪聲被測量算法的頻域卷積過程或時域最小二乘估計過程引入到頻率、幅值、相位測量結(jié)果中,使測量結(jié)果出現(xiàn)偏差,其中較小的偏差可視為測量數(shù)據(jù)的背景噪聲,較大的偏差或壞數(shù)據(jù)可視為測量數(shù)據(jù)的隨機脈沖噪聲［8-9］。均值濾波、局部回歸平滑等低通線性濾波方法適合過濾背景噪聲,對隨機脈沖噪聲缺乏魯棒性［10］。以中值濾波為代表的時域幾何成形濾波器對離群數(shù)據(jù)具有天然敏感性,對偏離正常數(shù)據(jù)序列的脈沖噪聲有很強的抑制作用［11］；但是當(dāng)測量數(shù)據(jù)快速變化時,容易出現(xiàn)平頂畸變且難以進行頻域補償［12］。加權(quán)中值濾波［13］、自適應(yīng)加權(quán)中值濾波［14］等改進中值濾波方法并沒有突破中值濾波的基本原理,所以濾波畸變問題無法徹底解決。魯棒局部回歸平滑（robust local regression smoothing,RLRS）方法是一種基于離群數(shù)據(jù)檢測與迭代加權(quán)修正的濾波方法,通過迭代調(diào)整權(quán)重矩陣,提高對隨機脈沖噪聲的魯棒性,整體呈非線性濾波特性,是對隨機脈沖噪聲與背景噪聲均具有良好降噪能力的常用降噪算法,應(yīng)用范圍覆蓋電力系統(tǒng)網(wǎng)絡(luò)攻擊診斷［15］、計算機視覺檢測［16］、多角度遙感圖像魯棒辨識［17］、腦磁信號壞數(shù)據(jù)剔除與修復(fù)［18］、頻域信號提?。?9-21］、地形數(shù)據(jù)修復(fù)［22］等多個領(lǐng)域。但由于RLRS 方法的局部降噪過程每次迭代都要執(zhí)行回歸平滑、脈沖噪聲檢測、魯棒平滑修復(fù)3 個步驟,每次更新擬合權(quán)重魯棒修正系數(shù)后都要重新計算擬合參數(shù)的加權(quán)最小二乘解,其矩陣乘法與求逆運算過程消耗大量計算資源,降噪效率低。同時,RLRS降噪數(shù)據(jù)在階躍跳變附近會發(fā)生畸變,一定程度上會失去部分時域原始信息。此外,RLRS 方法的局部加權(quán)修正系數(shù)由局部平滑殘余決定,對降噪數(shù)據(jù)序列造成的頻域畸變具有隨機性,補償困難。

針對上述問題,本文首先簡述RLRS 方法,然后提出快速魯棒回歸平滑（fast robust regression smoothing,FRRS）方法,最后通過仿真與實測數(shù)據(jù)分析驗證所提方法的降噪性能,并通過算法計算復(fù)雜度理論分析與算法耗時測試闡釋所提方法的計算速度優(yōu)勢。所提FRRS 快速降噪方法的主要特點為:1）實現(xiàn)了全局快速回歸平滑、隨機脈沖噪聲自適應(yīng)檢測、局部平滑畸變快速修復(fù),極大降低了計算復(fù)雜性；2）對背景噪聲與隨機脈沖噪聲均具有良好的降噪能力；3）有效辨識擾動階躍突變引起的偽脈沖噪聲,提升了降噪過程的擾動響應(yīng)原始信息保留能力。

1 RLRS 方法

設(shè)局部數(shù)據(jù)序列為D=[(xn,yn)],n=1,2,…,N,其中,xn為第n個數(shù)據(jù)的位置參數(shù),yn為第n個數(shù)據(jù)的數(shù)值,N為序列窗口長度且通常為奇數(shù)。

局部回歸平滑本質(zhì)上是一種多項式擬合平滑過程,可表示為Y=Xα。其中,X∈RN×(M+1)為多項式位置參數(shù)矩陣,其第n行、第m列元素為xn,m-1,M為多項式階數(shù)；Y∈RN×1為局部數(shù)據(jù)矩陣,其第n行元素為yn；α∈R(M+1)×1為多項式擬合參數(shù)矩陣,其第m行為第m-1 階擬合參數(shù)αm-1。利用最小二乘法由式（1）求解α得到αLS,進而由式（2）得到平滑后的數(shù)據(jù)序列YLS。

最小二乘的本質(zhì)是使總體偏差平方和最小,因此局部回歸平滑能夠有效抑制覆蓋整個時域的背景噪聲,但對離群異常數(shù)據(jù)不敏感,導(dǎo)致隨機脈沖噪聲平滑殘差較大,進而引起脈沖噪聲附近降噪畸變。而RLRS 本質(zhì)上則是一種通過自適應(yīng)迭代調(diào)整擬合權(quán)重實現(xiàn)離群異常數(shù)據(jù)魯棒平滑的多項式擬合平滑過程,既保留了最小二乘過程對背景噪聲良好的抑制能力,又通過調(diào)整擬合權(quán)重提高了對隨機脈沖噪聲的魯棒性。RLRS 方法運用加權(quán)最小二乘法求解α得到αWLS,即

式中:Δ∈RN×N為魯棒修正系數(shù)對角陣,其第n行、第n列元素為第n點擬合權(quán)重魯棒修正系數(shù)δn且初始值為1；W∈RN×N為擬合權(quán)重對角陣,其第n行、第n列元素為第n點擬合權(quán)重wn,通常取三次核函數(shù)值,如式（4）所示。

式中:xs為濾波位置；d(xs)為xs到序列窗口內(nèi)其他點的最遠距離,即

利用αWLS由式（6）計算平滑后的數(shù)據(jù)序列YWLS,并由式（7）計算平滑殘余R∈RN×1,其第n行元素為rn。

取R的中值λ,將6λ作為離群值判斷門檻,由式（8）更新Δ,若平滑殘余均小于6λ則停止迭代并輸出YWLS,否則由式（3）重新計算αWLS與YWLS。RLRS 算法迭代執(zhí)行次數(shù)表示為TR。

由于每次迭代更新擬合權(quán)重魯棒修正系數(shù)后都要重新計算αWLS,其矩陣乘法與求逆運算過程消耗大量計算資源,降噪效率低。有限長度窗口的局部平滑殘余中值受局部幾何特征偶然性影響較大,易出現(xiàn)中值偏大情況,導(dǎo)致離群值判斷門檻可靠性較低,難以表征背景噪聲強度的正常水平,引起幅值較小的隨機脈沖噪聲被誤判為背景噪聲。另外,RLRS 方法無法辨識擾動階躍突變引起的偽脈沖噪聲,導(dǎo)致降噪數(shù)據(jù)在階躍跳變附近發(fā)生畸變,一定程度上會失去部分時域原始信息。由于擬合權(quán)重魯棒修正系數(shù)與平滑殘差及其中值有關(guān),頻域畸變呈現(xiàn)隨機性,導(dǎo)致振蕩序列分析結(jié)果補償困難。本文基于上述分析,提出一種FRRS 方法,在保持良好的背景噪聲與隨機脈沖噪聲降噪能力的同時,降低計算復(fù)雜性,并提升擾動階躍突變辨識能力,提高降噪畸變可補償性。

2 FRRS 方法

所提FRRS 方法分為3 個部分:局部回歸平滑驅(qū)動的快速滑動濾波、全局平滑殘差驅(qū)動的離群數(shù)據(jù)自適應(yīng)檢測和離群數(shù)據(jù)局部平滑畸變快速修復(fù)。

2.1 局部回歸平滑驅(qū)動的快速滑動濾波

根據(jù)式（1）與式（2）,局部回歸平 滑形式可寫成:

式中:ΨLS∈RN×N為局部回歸平滑參數(shù)矩陣,計算公式如式（10）所示。

由于窗口中點的低通濾波特性最佳,即能夠達到最小幅頻波動、最大高頻衰減、零相位畸變,故取濾波位置xs為窗口中點,對應(yīng)的局部回歸平滑參數(shù)為φLS∈R1×N,即濾波位置對應(yīng)的ΨLS行參數(shù)。

設(shè)測量數(shù)據(jù)序列Z∈R1×(L+N)為:

式中:S∈RN×1為以l為中點的局部數(shù)據(jù)序列。

2.2 全局平滑殘差驅(qū)動的離群數(shù)據(jù)自適應(yīng)檢測

計算全局平滑殘差RZ∈R1×L,

RZ=[rZ,l]l=0,1,…,L-1 （15）

式中:rZ,l為第l點平滑殘差值,如式（16）所示。

根據(jù)文獻［23］諧波（間諧波）自適應(yīng)檢測門檻值,建立自適應(yīng)離群值判斷門檻τ,即

式中:μ和σ分別為RZ所含背景噪聲的均值與標(biāo)準(zhǔn)差估計值；η為魯棒調(diào)節(jié)系數(shù)。

假設(shè)背景噪聲為高斯白噪聲,η=6 時離群值誤判概率僅為2×10-9,故本文取η=6。τ的迭代估計過程可簡述為:首先,將RZ的所有元素視為背景噪聲,估計μ與σ,進而得到τ；然后,剔除RZ中高于τ的元素并重新估計τ,直到不存在高于τ的元素為止,得到最終的τ。迭代執(zhí)行次數(shù)表示為Tτ。取RZ中高于τ的極大值點為離群數(shù)據(jù)位置點,表示為:

U=[uh]h=1,2,…,H（18）

式中:uh為第h組離群數(shù)據(jù)位置；H為組數(shù)。

2.3 離群數(shù)據(jù)局部平滑畸變快速修復(fù)

根據(jù)式（3）和式（6）,RLRS 形式可寫成式（19）,其中,ΨWLS∈RN×N為RLRS 參數(shù)矩陣且由式（20）計算。

平滑畸變修復(fù)過程,以第h組離群數(shù)據(jù)為例,設(shè)S0∈R1×N0為非階躍跳變離群數(shù)據(jù)魯棒回歸平滑結(jié)果,可由式（24）計算得到,其中,SOutlier∈RN×1為以uh為中點的局部數(shù)據(jù)序列,如式（25）所示。

2.4 FRRS 方法整體流程

所提FRRS 方法的整體流程如圖1 所示。首先,通過局部回歸平滑實現(xiàn)針對測量數(shù)據(jù)序列背景噪聲的降噪處理；然后,根據(jù)全局平滑殘差自適應(yīng)檢測離群數(shù)據(jù)；最后,針對離群數(shù)據(jù)進行魯棒回歸平滑降噪并滑動修復(fù)離群數(shù)據(jù)造成的局部平滑畸變。

圖1 所提FRRS 方法整體流程圖Fig.1 Overall flow chart of the proposed FRRS method

3 算法分析與驗證

3.1 降噪性能分析

電力系統(tǒng)頻率測量曲線特征主要包括:1）負荷與發(fā)電有功功率嚴格平衡情況下的恒定值；2）準(zhǔn)穩(wěn)態(tài)情況下的緩慢斜坡變化；3）大擾動后快速斜坡變化及振蕩過程。相位由頻率偏移量的積分與系統(tǒng)物理參數(shù)的變化共同決定,相位測量曲線特征主要包括:1）額定頻率情況下的恒定值；2）頻率偏移量恒定情況下的斜坡變化；3）頻率斜坡變化情況下的二次函數(shù)變化；4）頻率振蕩情況下的振蕩變化；5）系統(tǒng)物理參數(shù)突變引起的階躍變化。幅值測量曲線特征主要包括:1）系統(tǒng)物理參數(shù)突變引起的階躍變化；2）負荷水平波動引起的幅值緩慢波動；3）大擾動引起的幅值快速斜坡變化與振蕩變化。

綜上,測量數(shù)據(jù)曲線特征應(yīng)包括常數(shù)值、斜坡、二次函數(shù)、振蕩、階躍等。其中,常數(shù)值、斜坡、二次函數(shù)等曲線特征與回歸平滑數(shù)學(xué)模型完全一致,降噪后測量曲線所包含的原始信息理論上能夠完全被保留,回歸平滑階數(shù)應(yīng)不小于2 階。而振蕩過程與回歸平滑數(shù)學(xué)模型并不一致,測量曲線所包含的振蕩信息經(jīng)降噪后將發(fā)生畸變,具體畸變程度由平滑過程的頻率響應(yīng)特性決定。設(shè)Pω為頻率ω處的頻率響應(yīng),計算方法為Pω=φLSΩω,其中,Ωω∈CN×1為頻率響應(yīng)計算因子矩陣,其第n行參數(shù)為e-jωn?？衫肞ω對后續(xù)振蕩分析結(jié)果進行補償。

本文采用噪聲強度抑制百分比β作為降噪能力量化指標(biāo),表示被濾除噪聲強度占原噪聲強度的百分比,β越大則降噪效果越好。背景噪聲降噪能力由平滑過程頻率響應(yīng)決定,當(dāng)回歸平滑階數(shù)M=2時,背景噪聲的β理論值如圖2 所示,β隨著局部窗長N增大而增大,當(dāng)N≥23 時β可達90%以上。而隨機脈沖噪聲抑制能力與魯棒修正系數(shù)直接相關(guān),零權(quán)強制噪聲抑制過程理論上可實現(xiàn)隨機脈沖噪聲100%降噪。

圖2 所提FRRS 方法背景噪聲強度抑制百分比理論值Fig.2 Theoretical value of background noise intensity suppression percentage of the proposed FRRS method

3.2 仿真驗證

3.2.1 仿真模型

根據(jù)前述電力系統(tǒng)測量數(shù)據(jù)曲線特征分析,測量數(shù)據(jù)曲線仿真模型設(shè)置如表1 所示。其中,a0、a1、a2分別為常數(shù)項參數(shù)、斜坡斜率、二次項參數(shù)；Ts為采樣間隔；f為振蕩頻率。

表1 測量數(shù)據(jù)曲線仿真模型及參數(shù)Table 1 Simulation models and parameters of measurement data curves

式中:AIPN為脈沖噪聲幅度；tARR,i和TIPN,i分別為第i個隨機脈沖噪聲的到達時間與持續(xù)時間。本文設(shè)置AIPN為5～10 之間的隨機數(shù),設(shè)置tARR,i隨機分布且覆蓋率為1%,設(shè)置TIPN,i=0.2 s。

擾動階躍跳變模型如式（31）所示:

式中:ASTEP為階躍幅度；tSTEP為階躍發(fā)生時間。本文設(shè)置ASTEP=10,tSTEP=30 s。

3.2.2 仿真結(jié)果

降噪?yún)?shù)設(shè)置為N=23,N0=5,M=2。噪聲強度抑制百分比如表2 所示,降噪結(jié)果如圖3所示。

表2 噪聲強度抑制百分比結(jié)果Table 2 Results of noise intensity suppression percentage

圖3 仿真數(shù)據(jù)降噪結(jié)果Fig.3 Denoising results of simulation data

由表2 和圖3 結(jié)果可以得出:所提FRRS 方法的背景噪聲強度抑制百分比能夠達到89.43%,與理論值90.18%相差不大,比RLRS 方法提高約3%；所提FRRS 方法的隨機脈沖噪聲強度抑制百分比達到98.31%,比RLRS 方法提高約1.2%；2 種方法的脈沖噪聲抑制結(jié)果均無明顯畸變,所提FRRS 方法能夠完整地保留階躍跳變特征并抑制與階躍跳變同時發(fā)生的脈沖噪聲,而RLRS 方法的降噪結(jié)果在階躍跳變附近發(fā)生畸變,一定程度上失去了部分時域原始信息。需要說明的是,在實際應(yīng)用所述FRRS 方法時,零權(quán)擬合長度參數(shù)N0需要預(yù)先設(shè)定為不小于隨機脈沖噪聲持續(xù)時間的數(shù)值,因此,測量數(shù)據(jù)的隨機脈沖噪聲持續(xù)時間統(tǒng)計特征應(yīng)為先驗已知信息。通過設(shè)置較大的N0參數(shù),本文應(yīng)用所述FRRS 方法分析了連續(xù)一個月的測量數(shù)據(jù)。分析結(jié)果表明,頻率、幅值、相位測量數(shù)據(jù)的隨機脈沖噪聲持續(xù)時間分別不大于0.4、0.7、0.6 s。因此,在實際應(yīng)用所提方法時可預(yù)先設(shè)定N0參數(shù),使N0Ts不小于0.7 s。

3.3 實測數(shù)據(jù)驗證

本節(jié)選取一組實測數(shù)據(jù)驗證所提方法的有效性,數(shù)據(jù)源為輕型廣域測量系統(tǒng)（wide area measurement system light,WAMS Light）［5］。圖4 是 華 東 電 網(wǎng) 與“兩華”同步電網(wǎng)某直流輸電閉鎖前后WAMS Light在上海、武漢、濟南獲得的測量數(shù)據(jù)及其降噪結(jié)果,閉鎖擾動發(fā)生后,華東電網(wǎng)有功功率缺失導(dǎo)致頻率、電壓幅值快速下降然后緩慢爬升,擾動后約13 s 上海測點電壓發(fā)生階躍突增；而“兩華”同步電網(wǎng)失去有功負荷后進入衰減振蕩過程。由圖4 可以看出,所提FRRS 方法與傳統(tǒng)RLRS 方法均能有效抑制背景噪聲與高強度隨機脈沖噪聲。尤其從圖4（a）與（b）的局部放大區(qū)域可以看出,快速變化數(shù)據(jù)的FRRS 降噪畸變明顯小于傳統(tǒng)RLRS 降噪畸變,所以FRRS 方法比傳統(tǒng)RLRS 方法能夠更加準(zhǔn)確地保留擾動階躍突變信息與擾動后快速斜坡變化信息。

圖4 實測數(shù)據(jù)降噪結(jié)果Fig.4 Denoising results of field measurement data

3.4 計算復(fù)雜度分析

傳統(tǒng)RLRS 方法與本文所提FRRS 方法的計算復(fù)雜度分析結(jié)果如表3 所示。傳統(tǒng)RLRS 方法每次迭代需要進行5 次矩陣相乘與1 次矩陣求逆運算；而所提FRRS 方法每個滑動窗口僅需要1 次矩陣相乘運算,若存在離群數(shù)據(jù),每組離群數(shù)據(jù)局部平滑畸變快速修復(fù)過程也僅需要1 次用于抑制離群數(shù)據(jù)平滑畸變的矩陣乘法運算并重新平滑離群數(shù)據(jù)局部N點數(shù)據(jù)。因此,所提FRRS 方法計算復(fù)雜度遠小于傳統(tǒng)RLRS 方法。2 種算法耗時測試結(jié)果如圖5 所示,測試環(huán)境為MATLAB 2016b,處理器為Intel Core i7-7700,降噪?yún)?shù)仍設(shè)置為N=23,N0=5,M=2,數(shù)據(jù)時間長度從2 h 到24 h,取100 次重復(fù)測試結(jié)果的均值。所提FRRS 方法的降噪耗時明顯低于傳統(tǒng)RLRS 方法,降噪速度提高了約20 倍,且隨著數(shù)據(jù)長度增大,降噪速度提升倍數(shù)沒有明顯變化,進一步驗證了所提FRRS 方法在計算效率方面的優(yōu)勢。

表3 計算復(fù)雜度分析結(jié)果Table 3 Analysis results of computational complexity

圖5 2 種算法耗時測試結(jié)果Fig.5 Test results of time consuming for two methods

4 結(jié)語

針對電力系統(tǒng)測量數(shù)據(jù)受強噪聲污染問題,本文提出了一種FRRS 降噪方法,相比傳統(tǒng)RLRS 方法,所提方法不僅保持了良好的背景噪聲與高強度隨機脈沖噪聲降噪能力,還能夠有效保留擾動階躍突變原始信息。

由于所述FRRS 方法實現(xiàn)了全局快速回歸平滑、隨機脈沖噪聲自適應(yīng)檢測、局部平滑畸變快速修復(fù),極大降低了計算復(fù)雜性,具有在線降噪應(yīng)用潛力。需要指出,所提FRRS 方法的零權(quán)擬合長度參數(shù)選擇依賴隨機脈沖噪聲持續(xù)時間統(tǒng)計特征先驗信息,隨著對電力系統(tǒng)測量數(shù)據(jù)噪聲特征統(tǒng)計規(guī)律的不斷認知,可逐漸優(yōu)化零權(quán)擬合長度參數(shù)選擇,或者進一步提出能夠自適應(yīng)調(diào)整零權(quán)擬合長度參數(shù)的FRRS 改進方法。