靳宗帥,張恒旭
(電網(wǎng)智能化調(diào)度與控制教育部重點(diǎn)實(shí)驗(yàn)室(山東大學(xué)),山東省濟(jì)南市 250061)
隨著大規(guī)模風(fēng)電/光伏等分布式發(fā)電、柔性直流、電動(dòng)汽車/儲(chǔ)能等大功率互動(dòng)性多元電力電子設(shè)備的接入[1-2],電網(wǎng)呈現(xiàn)有源化、電力電子化,信號(hào)形態(tài)、故障形態(tài)發(fā)生顯著變化,動(dòng)態(tài)行為更加復(fù)雜[3]。以同步相量測(cè)量單元(synchrophasor measurement unit,PMU)為終端的廣域測(cè)量系統(tǒng)(wide area measurement system,WAMS)[4]目前主要覆蓋高壓輸電網(wǎng);隨著配電網(wǎng)測(cè)量技術(shù)逐漸成熟[5-9],其提供的中低壓電網(wǎng)頻率、幅值、相位測(cè)量數(shù)據(jù)既蘊(yùn)含高壓主網(wǎng)響應(yīng)信息,又保留了局部響應(yīng)信息,作為WAMS 的補(bǔ)充數(shù)據(jù)源,實(shí)現(xiàn)了電力系統(tǒng)多電壓等級(jí)全景式同步監(jiān)測(cè)。但測(cè)量點(diǎn)與負(fù)荷的電氣距離較近,測(cè)量數(shù)據(jù)受噪聲污染嚴(yán)重,如何快速降噪是亟待解決的問題。
受噪聲污染的信號(hào)經(jīng)互感器接入測(cè)量裝置,然后由模數(shù)轉(zhuǎn)換模塊轉(zhuǎn)換為離散采樣序列,其所包含的噪聲被測(cè)量算法的頻域卷積過程或時(shí)域最小二乘估計(jì)過程引入到頻率、幅值、相位測(cè)量結(jié)果中,使測(cè)量結(jié)果出現(xiàn)偏差,其中較小的偏差可視為測(cè)量數(shù)據(jù)的背景噪聲,較大的偏差或壞數(shù)據(jù)可視為測(cè)量數(shù)據(jù)的隨機(jī)脈沖噪聲[8-9]。均值濾波、局部回歸平滑等低通線性濾波方法適合過濾背景噪聲,對(duì)隨機(jī)脈沖噪聲缺乏魯棒性[10]。以中值濾波為代表的時(shí)域幾何成形濾波器對(duì)離群數(shù)據(jù)具有天然敏感性,對(duì)偏離正常數(shù)據(jù)序列的脈沖噪聲有很強(qiáng)的抑制作用[11];但是當(dāng)測(cè)量數(shù)據(jù)快速變化時(shí),容易出現(xiàn)平頂畸變且難以進(jìn)行頻域補(bǔ)償[12]。加權(quán)中值濾波[13]、自適應(yīng)加權(quán)中值濾波[14]等改進(jìn)中值濾波方法并沒有突破中值濾波的基本原理,所以濾波畸變問題無法徹底解決。魯棒局部回歸平滑(robust local regression smoothing,RLRS)方法是一種基于離群數(shù)據(jù)檢測(cè)與迭代加權(quán)修正的濾波方法,通過迭代調(diào)整權(quán)重矩陣,提高對(duì)隨機(jī)脈沖噪聲的魯棒性,整體呈非線性濾波特性,是對(duì)隨機(jī)脈沖噪聲與背景噪聲均具有良好降噪能力的常用降噪算法,應(yīng)用范圍覆蓋電力系統(tǒng)網(wǎng)絡(luò)攻擊診斷[15]、計(jì)算機(jī)視覺檢測(cè)[16]、多角度遙感圖像魯棒辨識(shí)[17]、腦磁信號(hào)壞數(shù)據(jù)剔除與修復(fù)[18]、頻域信號(hào)提取[19-21]、地形數(shù)據(jù)修復(fù)[22]等多個(gè)領(lǐng)域。但由于RLRS 方法的局部降噪過程每次迭代都要執(zhí)行回歸平滑、脈沖噪聲檢測(cè)、魯棒平滑修復(fù)3 個(gè)步驟,每次更新擬合權(quán)重魯棒修正系數(shù)后都要重新計(jì)算擬合參數(shù)的加權(quán)最小二乘解,其矩陣乘法與求逆運(yùn)算過程消耗大量計(jì)算資源,降噪效率低。同時(shí),RLRS降噪數(shù)據(jù)在階躍跳變附近會(huì)發(fā)生畸變,一定程度上會(huì)失去部分時(shí)域原始信息。此外,RLRS 方法的局部加權(quán)修正系數(shù)由局部平滑殘余決定,對(duì)降噪數(shù)據(jù)序列造成的頻域畸變具有隨機(jī)性,補(bǔ)償困難。
針對(duì)上述問題,本文首先簡(jiǎn)述RLRS 方法,然后提出快速魯棒回歸平滑(fast robust regression smoothing,FRRS)方法,最后通過仿真與實(shí)測(cè)數(shù)據(jù)分析驗(yàn)證所提方法的降噪性能,并通過算法計(jì)算復(fù)雜度理論分析與算法耗時(shí)測(cè)試闡釋所提方法的計(jì)算速度優(yōu)勢(shì)。所提FRRS 快速降噪方法的主要特點(diǎn)為:1)實(shí)現(xiàn)了全局快速回歸平滑、隨機(jī)脈沖噪聲自適應(yīng)檢測(cè)、局部平滑畸變快速修復(fù),極大降低了計(jì)算復(fù)雜性;2)對(duì)背景噪聲與隨機(jī)脈沖噪聲均具有良好的降噪能力;3)有效辨識(shí)擾動(dòng)階躍突變引起的偽脈沖噪聲,提升了降噪過程的擾動(dòng)響應(yīng)原始信息保留能力。
設(shè)局部數(shù)據(jù)序列為D=[(xn,yn)],n=1,2,…,N,其中,xn為第n個(gè)數(shù)據(jù)的位置參數(shù),yn為第n個(gè)數(shù)據(jù)的數(shù)值,N為序列窗口長度且通常為奇數(shù)。
局部回歸平滑本質(zhì)上是一種多項(xiàng)式擬合平滑過程,可表示為Y=Xα。其中,X∈RN×(M+1)為多項(xiàng)式位置參數(shù)矩陣,其第n行、第m列元素為xn,m-1,M為多項(xiàng)式階數(shù);Y∈RN×1為局部數(shù)據(jù)矩陣,其第n行元素為yn;α∈R(M+1)×1為多項(xiàng)式擬合參數(shù)矩陣,其第m行為第m-1 階擬合參數(shù)αm-1。利用最小二乘法由式(1)求解α得到αLS,進(jìn)而由式(2)得到平滑后的數(shù)據(jù)序列YLS。
最小二乘的本質(zhì)是使總體偏差平方和最小,因此局部回歸平滑能夠有效抑制覆蓋整個(gè)時(shí)域的背景噪聲,但對(duì)離群異常數(shù)據(jù)不敏感,導(dǎo)致隨機(jī)脈沖噪聲平滑殘差較大,進(jìn)而引起脈沖噪聲附近降噪畸變。而RLRS 本質(zhì)上則是一種通過自適應(yīng)迭代調(diào)整擬合權(quán)重實(shí)現(xiàn)離群異常數(shù)據(jù)魯棒平滑的多項(xiàng)式擬合平滑過程,既保留了最小二乘過程對(duì)背景噪聲良好的抑制能力,又通過調(diào)整擬合權(quán)重提高了對(duì)隨機(jī)脈沖噪聲的魯棒性。RLRS 方法運(yùn)用加權(quán)最小二乘法求解α得到αWLS,即
式中:Δ∈RN×N為魯棒修正系數(shù)對(duì)角陣,其第n行、第n列元素為第n點(diǎn)擬合權(quán)重魯棒修正系數(shù)δn且初始值為1;W∈RN×N為擬合權(quán)重對(duì)角陣,其第n行、第n列元素為第n點(diǎn)擬合權(quán)重wn,通常取三次核函數(shù)值,如式(4)所示。
式中:xs為濾波位置;d(xs)為xs到序列窗口內(nèi)其他點(diǎn)的最遠(yuǎn)距離,即
利用αWLS由式(6)計(jì)算平滑后的數(shù)據(jù)序列YWLS,并由式(7)計(jì)算平滑殘余R∈RN×1,其第n行元素為rn。
取R的中值λ,將6λ作為離群值判斷門檻,由式(8)更新Δ,若平滑殘余均小于6λ則停止迭代并輸出YWLS,否則由式(3)重新計(jì)算αWLS與YWLS。RLRS 算法迭代執(zhí)行次數(shù)表示為TR。
由于每次迭代更新擬合權(quán)重魯棒修正系數(shù)后都要重新計(jì)算αWLS,其矩陣乘法與求逆運(yùn)算過程消耗大量計(jì)算資源,降噪效率低。有限長度窗口的局部平滑殘余中值受局部幾何特征偶然性影響較大,易出現(xiàn)中值偏大情況,導(dǎo)致離群值判斷門檻可靠性較低,難以表征背景噪聲強(qiáng)度的正常水平,引起幅值較小的隨機(jī)脈沖噪聲被誤判為背景噪聲。另外,RLRS 方法無法辨識(shí)擾動(dòng)階躍突變引起的偽脈沖噪聲,導(dǎo)致降噪數(shù)據(jù)在階躍跳變附近發(fā)生畸變,一定程度上會(huì)失去部分時(shí)域原始信息。由于擬合權(quán)重魯棒修正系數(shù)與平滑殘差及其中值有關(guān),頻域畸變呈現(xiàn)隨機(jī)性,導(dǎo)致振蕩序列分析結(jié)果補(bǔ)償困難。本文基于上述分析,提出一種FRRS 方法,在保持良好的背景噪聲與隨機(jī)脈沖噪聲降噪能力的同時(shí),降低計(jì)算復(fù)雜性,并提升擾動(dòng)階躍突變辨識(shí)能力,提高降噪畸變可補(bǔ)償性。
所提FRRS 方法分為3 個(gè)部分:局部回歸平滑驅(qū)動(dòng)的快速滑動(dòng)濾波、全局平滑殘差驅(qū)動(dòng)的離群數(shù)據(jù)自適應(yīng)檢測(cè)和離群數(shù)據(jù)局部平滑畸變快速修復(fù)。
根據(jù)式(1)與式(2),局部回歸平 滑形式可寫成:
式中:ΨLS∈RN×N為局部回歸平滑參數(shù)矩陣,計(jì)算公式如式(10)所示。
由于窗口中點(diǎn)的低通濾波特性最佳,即能夠達(dá)到最小幅頻波動(dòng)、最大高頻衰減、零相位畸變,故取濾波位置xs為窗口中點(diǎn),對(duì)應(yīng)的局部回歸平滑參數(shù)為φLS∈R1×N,即濾波位置對(duì)應(yīng)的ΨLS行參數(shù)。
設(shè)測(cè)量數(shù)據(jù)序列Z∈R1×(L+N)為:
式中:S∈RN×1為以l為中點(diǎn)的局部數(shù)據(jù)序列。
計(jì)算全局平滑殘差RZ∈R1×L,
RZ=[rZ,l]l=0,1,…,L-1 (15)
式中:rZ,l為第l點(diǎn)平滑殘差值,如式(16)所示。
根據(jù)文獻(xiàn)[23]諧波(間諧波)自適應(yīng)檢測(cè)門檻值,建立自適應(yīng)離群值判斷門檻τ,即
式中:μ和σ分別為RZ所含背景噪聲的均值與標(biāo)準(zhǔn)差估計(jì)值;η為魯棒調(diào)節(jié)系數(shù)。
假設(shè)背景噪聲為高斯白噪聲,η=6 時(shí)離群值誤判概率僅為2×10-9,故本文取η=6。τ的迭代估計(jì)過程可簡(jiǎn)述為:首先,將RZ的所有元素視為背景噪聲,估計(jì)μ與σ,進(jìn)而得到τ;然后,剔除RZ中高于τ的元素并重新估計(jì)τ,直到不存在高于τ的元素為止,得到最終的τ。迭代執(zhí)行次數(shù)表示為Tτ。取RZ中高于τ的極大值點(diǎn)為離群數(shù)據(jù)位置點(diǎn),表示為:
U=[uh]h=1,2,…,H(18)
式中:uh為第h組離群數(shù)據(jù)位置;H為組數(shù)。
根據(jù)式(3)和式(6),RLRS 形式可寫成式(19),其中,ΨWLS∈RN×N為RLRS 參數(shù)矩陣且由式(20)計(jì)算。
平滑畸變修復(fù)過程,以第h組離群數(shù)據(jù)為例,設(shè)S0∈R1×N0為非階躍跳變離群數(shù)據(jù)魯棒回歸平滑結(jié)果,可由式(24)計(jì)算得到,其中,SOutlier∈RN×1為以u(píng)h為中點(diǎn)的局部數(shù)據(jù)序列,如式(25)所示。
所提FRRS 方法的整體流程如圖1 所示。首先,通過局部回歸平滑實(shí)現(xiàn)針對(duì)測(cè)量數(shù)據(jù)序列背景噪聲的降噪處理;然后,根據(jù)全局平滑殘差自適應(yīng)檢測(cè)離群數(shù)據(jù);最后,針對(duì)離群數(shù)據(jù)進(jìn)行魯棒回歸平滑降噪并滑動(dòng)修復(fù)離群數(shù)據(jù)造成的局部平滑畸變。
圖1 所提FRRS 方法整體流程圖Fig.1 Overall flow chart of the proposed FRRS method
電力系統(tǒng)頻率測(cè)量曲線特征主要包括:1)負(fù)荷與發(fā)電有功功率嚴(yán)格平衡情況下的恒定值;2)準(zhǔn)穩(wěn)態(tài)情況下的緩慢斜坡變化;3)大擾動(dòng)后快速斜坡變化及振蕩過程。相位由頻率偏移量的積分與系統(tǒng)物理參數(shù)的變化共同決定,相位測(cè)量曲線特征主要包括:1)額定頻率情況下的恒定值;2)頻率偏移量恒定情況下的斜坡變化;3)頻率斜坡變化情況下的二次函數(shù)變化;4)頻率振蕩情況下的振蕩變化;5)系統(tǒng)物理參數(shù)突變引起的階躍變化。幅值測(cè)量曲線特征主要包括:1)系統(tǒng)物理參數(shù)突變引起的階躍變化;2)負(fù)荷水平波動(dòng)引起的幅值緩慢波動(dòng);3)大擾動(dòng)引起的幅值快速斜坡變化與振蕩變化。
綜上,測(cè)量數(shù)據(jù)曲線特征應(yīng)包括常數(shù)值、斜坡、二次函數(shù)、振蕩、階躍等。其中,常數(shù)值、斜坡、二次函數(shù)等曲線特征與回歸平滑數(shù)學(xué)模型完全一致,降噪后測(cè)量曲線所包含的原始信息理論上能夠完全被保留,回歸平滑階數(shù)應(yīng)不小于2 階。而振蕩過程與回歸平滑數(shù)學(xué)模型并不一致,測(cè)量曲線所包含的振蕩信息經(jīng)降噪后將發(fā)生畸變,具體畸變程度由平滑過程的頻率響應(yīng)特性決定。設(shè)Pω為頻率ω處的頻率響應(yīng),計(jì)算方法為Pω=φLSΩω,其中,Ωω∈CN×1為頻率響應(yīng)計(jì)算因子矩陣,其第n行參數(shù)為e-jωn??衫肞ω對(duì)后續(xù)振蕩分析結(jié)果進(jìn)行補(bǔ)償。
本文采用噪聲強(qiáng)度抑制百分比β作為降噪能力量化指標(biāo),表示被濾除噪聲強(qiáng)度占原噪聲強(qiáng)度的百分比,β越大則降噪效果越好。背景噪聲降噪能力由平滑過程頻率響應(yīng)決定,當(dāng)回歸平滑階數(shù)M=2時(shí),背景噪聲的β理論值如圖2 所示,β隨著局部窗長N增大而增大,當(dāng)N≥23 時(shí)β可達(dá)90%以上。而隨機(jī)脈沖噪聲抑制能力與魯棒修正系數(shù)直接相關(guān),零權(quán)強(qiáng)制噪聲抑制過程理論上可實(shí)現(xiàn)隨機(jī)脈沖噪聲100%降噪。
圖2 所提FRRS 方法背景噪聲強(qiáng)度抑制百分比理論值Fig.2 Theoretical value of background noise intensity suppression percentage of the proposed FRRS method
3.2.1 仿真模型
根據(jù)前述電力系統(tǒng)測(cè)量數(shù)據(jù)曲線特征分析,測(cè)量數(shù)據(jù)曲線仿真模型設(shè)置如表1 所示。其中,a0、a1、a2分別為常數(shù)項(xiàng)參數(shù)、斜坡斜率、二次項(xiàng)參數(shù);Ts為采樣間隔;f為振蕩頻率。
表1 測(cè)量數(shù)據(jù)曲線仿真模型及參數(shù)Table 1 Simulation models and parameters of measurement data curves
式中:AIPN為脈沖噪聲幅度;tARR,i和TIPN,i分別為第i個(gè)隨機(jī)脈沖噪聲的到達(dá)時(shí)間與持續(xù)時(shí)間。本文設(shè)置AIPN為5~10 之間的隨機(jī)數(shù),設(shè)置tARR,i隨機(jī)分布且覆蓋率為1%,設(shè)置TIPN,i=0.2 s。
擾動(dòng)階躍跳變模型如式(31)所示:
式中:ASTEP為階躍幅度;tSTEP為階躍發(fā)生時(shí)間。本文設(shè)置ASTEP=10,tSTEP=30 s。
3.2.2 仿真結(jié)果
降噪?yún)?shù)設(shè)置為N=23,N0=5,M=2。噪聲強(qiáng)度抑制百分比如表2 所示,降噪結(jié)果如圖3所示。
表2 噪聲強(qiáng)度抑制百分比結(jié)果Table 2 Results of noise intensity suppression percentage
圖3 仿真數(shù)據(jù)降噪結(jié)果Fig.3 Denoising results of simulation data
由表2 和圖3 結(jié)果可以得出:所提FRRS 方法的背景噪聲強(qiáng)度抑制百分比能夠達(dá)到89.43%,與理論值90.18%相差不大,比RLRS 方法提高約3%;所提FRRS 方法的隨機(jī)脈沖噪聲強(qiáng)度抑制百分比達(dá)到98.31%,比RLRS 方法提高約1.2%;2 種方法的脈沖噪聲抑制結(jié)果均無明顯畸變,所提FRRS 方法能夠完整地保留階躍跳變特征并抑制與階躍跳變同時(shí)發(fā)生的脈沖噪聲,而RLRS 方法的降噪結(jié)果在階躍跳變附近發(fā)生畸變,一定程度上失去了部分時(shí)域原始信息。需要說明的是,在實(shí)際應(yīng)用所述FRRS 方法時(shí),零權(quán)擬合長度參數(shù)N0需要預(yù)先設(shè)定為不小于隨機(jī)脈沖噪聲持續(xù)時(shí)間的數(shù)值,因此,測(cè)量數(shù)據(jù)的隨機(jī)脈沖噪聲持續(xù)時(shí)間統(tǒng)計(jì)特征應(yīng)為先驗(yàn)已知信息。通過設(shè)置較大的N0參數(shù),本文應(yīng)用所述FRRS 方法分析了連續(xù)一個(gè)月的測(cè)量數(shù)據(jù)。分析結(jié)果表明,頻率、幅值、相位測(cè)量數(shù)據(jù)的隨機(jī)脈沖噪聲持續(xù)時(shí)間分別不大于0.4、0.7、0.6 s。因此,在實(shí)際應(yīng)用所提方法時(shí)可預(yù)先設(shè)定N0參數(shù),使N0Ts不小于0.7 s。
本節(jié)選取一組實(shí)測(cè)數(shù)據(jù)驗(yàn)證所提方法的有效性,數(shù)據(jù)源為輕型廣域測(cè)量系統(tǒng)(wide area measurement system light,WAMS Light)[5]。圖4 是 華 東 電 網(wǎng) 與“兩華”同步電網(wǎng)某直流輸電閉鎖前后WAMS Light在上海、武漢、濟(jì)南獲得的測(cè)量數(shù)據(jù)及其降噪結(jié)果,閉鎖擾動(dòng)發(fā)生后,華東電網(wǎng)有功功率缺失導(dǎo)致頻率、電壓幅值快速下降然后緩慢爬升,擾動(dòng)后約13 s 上海測(cè)點(diǎn)電壓發(fā)生階躍突增;而“兩華”同步電網(wǎng)失去有功負(fù)荷后進(jìn)入衰減振蕩過程。由圖4 可以看出,所提FRRS 方法與傳統(tǒng)RLRS 方法均能有效抑制背景噪聲與高強(qiáng)度隨機(jī)脈沖噪聲。尤其從圖4(a)與(b)的局部放大區(qū)域可以看出,快速變化數(shù)據(jù)的FRRS 降噪畸變明顯小于傳統(tǒng)RLRS 降噪畸變,所以FRRS 方法比傳統(tǒng)RLRS 方法能夠更加準(zhǔn)確地保留擾動(dòng)階躍突變信息與擾動(dòng)后快速斜坡變化信息。
圖4 實(shí)測(cè)數(shù)據(jù)降噪結(jié)果Fig.4 Denoising results of field measurement data
傳統(tǒng)RLRS 方法與本文所提FRRS 方法的計(jì)算復(fù)雜度分析結(jié)果如表3 所示。傳統(tǒng)RLRS 方法每次迭代需要進(jìn)行5 次矩陣相乘與1 次矩陣求逆運(yùn)算;而所提FRRS 方法每個(gè)滑動(dòng)窗口僅需要1 次矩陣相乘運(yùn)算,若存在離群數(shù)據(jù),每組離群數(shù)據(jù)局部平滑畸變快速修復(fù)過程也僅需要1 次用于抑制離群數(shù)據(jù)平滑畸變的矩陣乘法運(yùn)算并重新平滑離群數(shù)據(jù)局部N點(diǎn)數(shù)據(jù)。因此,所提FRRS 方法計(jì)算復(fù)雜度遠(yuǎn)小于傳統(tǒng)RLRS 方法。2 種算法耗時(shí)測(cè)試結(jié)果如圖5 所示,測(cè)試環(huán)境為MATLAB 2016b,處理器為Intel Core i7-7700,降噪?yún)?shù)仍設(shè)置為N=23,N0=5,M=2,數(shù)據(jù)時(shí)間長度從2 h 到24 h,取100 次重復(fù)測(cè)試結(jié)果的均值。所提FRRS 方法的降噪耗時(shí)明顯低于傳統(tǒng)RLRS 方法,降噪速度提高了約20 倍,且隨著數(shù)據(jù)長度增大,降噪速度提升倍數(shù)沒有明顯變化,進(jìn)一步驗(yàn)證了所提FRRS 方法在計(jì)算效率方面的優(yōu)勢(shì)。
表3 計(jì)算復(fù)雜度分析結(jié)果Table 3 Analysis results of computational complexity
圖5 2 種算法耗時(shí)測(cè)試結(jié)果Fig.5 Test results of time consuming for two methods
針對(duì)電力系統(tǒng)測(cè)量數(shù)據(jù)受強(qiáng)噪聲污染問題,本文提出了一種FRRS 降噪方法,相比傳統(tǒng)RLRS 方法,所提方法不僅保持了良好的背景噪聲與高強(qiáng)度隨機(jī)脈沖噪聲降噪能力,還能夠有效保留擾動(dòng)階躍突變?cè)夹畔ⅰ?/p>
由于所述FRRS 方法實(shí)現(xiàn)了全局快速回歸平滑、隨機(jī)脈沖噪聲自適應(yīng)檢測(cè)、局部平滑畸變快速修復(fù),極大降低了計(jì)算復(fù)雜性,具有在線降噪應(yīng)用潛力。需要指出,所提FRRS 方法的零權(quán)擬合長度參數(shù)選擇依賴隨機(jī)脈沖噪聲持續(xù)時(shí)間統(tǒng)計(jì)特征先驗(yàn)信息,隨著對(duì)電力系統(tǒng)測(cè)量數(shù)據(jù)噪聲特征統(tǒng)計(jì)規(guī)律的不斷認(rèn)知,可逐漸優(yōu)化零權(quán)擬合長度參數(shù)選擇,或者進(jìn)一步提出能夠自適應(yīng)調(diào)整零權(quán)擬合長度參數(shù)的FRRS 改進(jìn)方法。