电力系统强噪声测量数据快速鲁棒回归平滑降噪

靳宗帅，张恒旭

（电网智能化调度与控制教育部重点实验室（山东大学）,山东省济南市 250061）

0 引言

随着大规模风电/光伏等分布式发电、柔性直流、电动汽车/储能等大功率互动性多元电力电子设备的接入［1-2］,电网呈现有源化、电力电子化,信号形态、故障形态发生显著变化,动态行为更加复杂［3］。以同步相量测量单元（synchrophasor measurement unit,PMU）为终端的广域测量系统（wide area measurement system,WAMS）［4］目前主要覆盖高压输电网；随着配电网测量技术逐渐成熟［5-9］,其提供的中低压电网频率、幅值、相位测量数据既蕴含高压主网响应信息,又保留了局部响应信息,作为WAMS 的补充数据源,实现了电力系统多电压等级全景式同步监测。但测量点与负荷的电气距离较近,测量数据受噪声污染严重,如何快速降噪是亟待解决的问题。

受噪声污染的信号经互感器接入测量装置,然后由模数转换模块转换为离散采样序列,其所包含的噪声被测量算法的频域卷积过程或时域最小二乘估计过程引入到频率、幅值、相位测量结果中,使测量结果出现偏差,其中较小的偏差可视为测量数据的背景噪声,较大的偏差或坏数据可视为测量数据的随机脉冲噪声［8-9］。均值滤波、局部回归平滑等低通线性滤波方法适合过滤背景噪声,对随机脉冲噪声缺乏鲁棒性［10］。以中值滤波为代表的时域几何成形滤波器对离群数据具有天然敏感性,对偏离正常数据序列的脉冲噪声有很强的抑制作用［11］；但是当测量数据快速变化时,容易出现平顶畸变且难以进行频域补偿［12］。加权中值滤波［13］、自适应加权中值滤波［14］等改进中值滤波方法并没有突破中值滤波的基本原理,所以滤波畸变问题无法彻底解决。鲁棒局部回归平滑（robust local regression smoothing,RLRS）方法是一种基于离群数据检测与迭代加权修正的滤波方法,通过迭代调整权重矩阵,提高对随机脉冲噪声的鲁棒性,整体呈非线性滤波特性,是对随机脉冲噪声与背景噪声均具有良好降噪能力的常用降噪算法,应用范围覆盖电力系统网络攻击诊断［15］、计算机视觉检测［16］、多角度遥感图像鲁棒辨识［17］、脑磁信号坏数据剔除与修复［18］、频域信号提取［19-21］、地形数据修复［22］等多个领域。但由于RLRS 方法的局部降噪过程每次迭代都要执行回归平滑、脉冲噪声检测、鲁棒平滑修复3 个步骤,每次更新拟合权重鲁棒修正系数后都要重新计算拟合参数的加权最小二乘解,其矩阵乘法与求逆运算过程消耗大量计算资源,降噪效率低。同时,RLRS降噪数据在阶跃跳变附近会发生畸变,一定程度上会失去部分时域原始信息。此外,RLRS 方法的局部加权修正系数由局部平滑残余决定,对降噪数据序列造成的频域畸变具有随机性,补偿困难。

针对上述问题,本文首先简述RLRS 方法,然后提出快速鲁棒回归平滑（fast robust regression smoothing,FRRS）方法,最后通过仿真与实测数据分析验证所提方法的降噪性能,并通过算法计算复杂度理论分析与算法耗时测试阐释所提方法的计算速度优势。所提FRRS 快速降噪方法的主要特点为:1）实现了全局快速回归平滑、随机脉冲噪声自适应检测、局部平滑畸变快速修复,极大降低了计算复杂性；2）对背景噪声与随机脉冲噪声均具有良好的降噪能力；3）有效辨识扰动阶跃突变引起的伪脉冲噪声,提升了降噪过程的扰动响应原始信息保留能力。

1 RLRS 方法

设局部数据序列为D=[(xn,yn)],n=1,2,…,N,其中,xn为第n个数据的位置参数,yn为第n个数据的数值,N为序列窗口长度且通常为奇数。

局部回归平滑本质上是一种多项式拟合平滑过程,可表示为Y=Xα。其中,X∈RN×(M+1)为多项式位置参数矩阵,其第n行、第m列元素为xn,m-1,M为多项式阶数；Y∈RN×1为局部数据矩阵,其第n行元素为yn；α∈R(M+1)×1为多项式拟合参数矩阵,其第m行为第m-1 阶拟合参数αm-1。利用最小二乘法由式（1）求解α得到αLS,进而由式（2）得到平滑后的数据序列YLS。

最小二乘的本质是使总体偏差平方和最小,因此局部回归平滑能够有效抑制覆盖整个时域的背景噪声,但对离群异常数据不敏感,导致随机脉冲噪声平滑残差较大,进而引起脉冲噪声附近降噪畸变。而RLRS 本质上则是一种通过自适应迭代调整拟合权重实现离群异常数据鲁棒平滑的多项式拟合平滑过程,既保留了最小二乘过程对背景噪声良好的抑制能力,又通过调整拟合权重提高了对随机脉冲噪声的鲁棒性。RLRS 方法运用加权最小二乘法求解α得到αWLS,即

式中:Δ∈RN×N为鲁棒修正系数对角阵,其第n行、第n列元素为第n点拟合权重鲁棒修正系数δn且初始值为1；W∈RN×N为拟合权重对角阵,其第n行、第n列元素为第n点拟合权重wn,通常取三次核函数值,如式（4）所示。

式中:xs为滤波位置；d(xs)为xs到序列窗口内其他点的最远距离,即

利用αWLS由式（6）计算平滑后的数据序列YWLS,并由式（7）计算平滑残余R∈RN×1,其第n行元素为rn。

取R的中值λ,将6λ作为离群值判断门槛,由式（8）更新Δ,若平滑残余均小于6λ则停止迭代并输出YWLS,否则由式（3）重新计算αWLS与YWLS。RLRS 算法迭代执行次数表示为TR。

由于每次迭代更新拟合权重鲁棒修正系数后都要重新计算αWLS,其矩阵乘法与求逆运算过程消耗大量计算资源,降噪效率低。有限长度窗口的局部平滑残余中值受局部几何特征偶然性影响较大,易出现中值偏大情况,导致离群值判断门槛可靠性较低,难以表征背景噪声强度的正常水平,引起幅值较小的随机脉冲噪声被误判为背景噪声。另外,RLRS 方法无法辨识扰动阶跃突变引起的伪脉冲噪声,导致降噪数据在阶跃跳变附近发生畸变,一定程度上会失去部分时域原始信息。由于拟合权重鲁棒修正系数与平滑残差及其中值有关,频域畸变呈现随机性,导致振荡序列分析结果补偿困难。本文基于上述分析,提出一种FRRS 方法,在保持良好的背景噪声与随机脉冲噪声降噪能力的同时,降低计算复杂性,并提升扰动阶跃突变辨识能力,提高降噪畸变可补偿性。

2 FRRS 方法

所提FRRS 方法分为3 个部分:局部回归平滑驱动的快速滑动滤波、全局平滑残差驱动的离群数据自适应检测和离群数据局部平滑畸变快速修复。

2.1 局部回归平滑驱动的快速滑动滤波

根据式（1）与式（2）,局部回归平 滑形式可写成:

式中:ΨLS∈RN×N为局部回归平滑参数矩阵,计算公式如式（10）所示。

由于窗口中点的低通滤波特性最佳,即能够达到最小幅频波动、最大高频衰减、零相位畸变,故取滤波位置xs为窗口中点,对应的局部回归平滑参数为φLS∈R1×N,即滤波位置对应的ΨLS行参数。

设测量数据序列Z∈R1×(L+N)为:

式中:S∈RN×1为以l为中点的局部数据序列。

2.2 全局平滑残差驱动的离群数据自适应检测

计算全局平滑残差RZ∈R1×L,

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

式中:rZ,l为第l点平滑残差值,如式（16）所示。

根据文献［23］谐波（间谐波）自适应检测门槛值,建立自适应离群值判断门槛τ,即

式中:μ和σ分别为RZ所含背景噪声的均值与标准差估计值；η为鲁棒调节系数。

假设背景噪声为高斯白噪声,η=6 时离群值误判概率仅为2×10-9,故本文取η=6。τ的迭代估计过程可简述为:首先,将RZ的所有元素视为背景噪声,估计μ与σ,进而得到τ；然后,剔除RZ中高于τ的元素并重新估计τ,直到不存在高于τ的元素为止,得到最终的τ。迭代执行次数表示为Tτ。取RZ中高于τ的极大值点为离群数据位置点,表示为:

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

式中:uh为第h组离群数据位置；H为组数。

2.3 离群数据局部平滑畸变快速修复

根据式（3）和式（6）,RLRS 形式可写成式（19）,其中,ΨWLS∈RN×N为RLRS 参数矩阵且由式（20）计算。

平滑畸变修复过程,以第h组离群数据为例,设S0∈R1×N0为非阶跃跳变离群数据鲁棒回归平滑结果,可由式（24）计算得到,其中,SOutlier∈RN×1为以uh为中点的局部数据序列,如式（25）所示。

2.4 FRRS 方法整体流程

所提FRRS 方法的整体流程如图1 所示。首先,通过局部回归平滑实现针对测量数据序列背景噪声的降噪处理；然后,根据全局平滑残差自适应检测离群数据；最后,针对离群数据进行鲁棒回归平滑降噪并滑动修复离群数据造成的局部平滑畸变。

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

3 算法分析与验证

3.1 降噪性能分析

电力系统频率测量曲线特征主要包括:1）负荷与发电有功功率严格平衡情况下的恒定值；2）准稳态情况下的缓慢斜坡变化；3）大扰动后快速斜坡变化及振荡过程。相位由频率偏移量的积分与系统物理参数的变化共同决定,相位测量曲线特征主要包括:1）额定频率情况下的恒定值；2）频率偏移量恒定情况下的斜坡变化；3）频率斜坡变化情况下的二次函数变化；4）频率振荡情况下的振荡变化；5）系统物理参数突变引起的阶跃变化。幅值测量曲线特征主要包括:1）系统物理参数突变引起的阶跃变化；2）负荷水平波动引起的幅值缓慢波动；3）大扰动引起的幅值快速斜坡变化与振荡变化。

综上,测量数据曲线特征应包括常数值、斜坡、二次函数、振荡、阶跃等。其中,常数值、斜坡、二次函数等曲线特征与回归平滑数学模型完全一致,降噪后测量曲线所包含的原始信息理论上能够完全被保留,回归平滑阶数应不小于2 阶。而振荡过程与回归平滑数学模型并不一致,测量曲线所包含的振荡信息经降噪后将发生畸变,具体畸变程度由平滑过程的频率响应特性决定。设Pω为频率ω处的频率响应,计算方法为Pω=φLSΩω,其中,Ωω∈CN×1为频率响应计算因子矩阵,其第n行参数为e-jωn。可利用Pω对后续振荡分析结果进行补偿。

本文采用噪声强度抑制百分比β作为降噪能力量化指标,表示被滤除噪声强度占原噪声强度的百分比,β越大则降噪效果越好。背景噪声降噪能力由平滑过程频率响应决定,当回归平滑阶数M=2时,背景噪声的β理论值如图2 所示,β随着局部窗长N增大而增大,当N≥23 时β可达90%以上。而随机脉冲噪声抑制能力与鲁棒修正系数直接相关,零权强制噪声抑制过程理论上可实现随机脉冲噪声100%降噪。

图2 所提FRRS 方法背景噪声强度抑制百分比理论值Fig.2 Theoretical value of background noise intensity suppression percentage of the proposed FRRS method

3.2 仿真验证

3.2.1 仿真模型

根据前述电力系统测量数据曲线特征分析,测量数据曲线仿真模型设置如表1 所示。其中,a0、a1、a2分别为常数项参数、斜坡斜率、二次项参数；Ts为采样间隔；f为振荡频率。

表1 测量数据曲线仿真模型及参数Table 1 Simulation models and parameters of measurement data curves

式中:AIPN为脉冲噪声幅度；tARR,i和TIPN,i分别为第i个随机脉冲噪声的到达时间与持续时间。本文设置AIPN为5～10 之间的随机数,设置tARR,i随机分布且覆盖率为1%,设置TIPN,i=0.2 s。

扰动阶跃跳变模型如式（31）所示:

式中:ASTEP为阶跃幅度；tSTEP为阶跃发生时间。本文设置ASTEP=10,tSTEP=30 s。

3.2.2 仿真结果

降噪参数设置为N=23,N0=5,M=2。噪声强度抑制百分比如表2 所示,降噪结果如图3所示。

表2 噪声强度抑制百分比结果Table 2 Results of noise intensity suppression percentage

图3 仿真数据降噪结果Fig.3 Denoising results of simulation data

由表2 和图3 结果可以得出:所提FRRS 方法的背景噪声强度抑制百分比能够达到89.43%,与理论值90.18%相差不大,比RLRS 方法提高约3%；所提FRRS 方法的随机脉冲噪声强度抑制百分比达到98.31%,比RLRS 方法提高约1.2%；2 种方法的脉冲噪声抑制结果均无明显畸变,所提FRRS 方法能够完整地保留阶跃跳变特征并抑制与阶跃跳变同时发生的脉冲噪声,而RLRS 方法的降噪结果在阶跃跳变附近发生畸变,一定程度上失去了部分时域原始信息。需要说明的是,在实际应用所述FRRS 方法时,零权拟合长度参数N0需要预先设定为不小于随机脉冲噪声持续时间的数值,因此,测量数据的随机脉冲噪声持续时间统计特征应为先验已知信息。通过设置较大的N0参数,本文应用所述FRRS 方法分析了连续一个月的测量数据。分析结果表明,频率、幅值、相位测量数据的随机脉冲噪声持续时间分别不大于0.4、0.7、0.6 s。因此,在实际应用所提方法时可预先设定N0参数,使N0Ts不小于0.7 s。

3.3 实测数据验证

本节选取一组实测数据验证所提方法的有效性,数据源为轻型广域测量系统（wide area measurement system light,WAMS Light）［5］。图4 是 华 东 电 网 与“两华”同步电网某直流输电闭锁前后WAMS Light在上海、武汉、济南获得的测量数据及其降噪结果,闭锁扰动发生后,华东电网有功功率缺失导致频率、电压幅值快速下降然后缓慢爬升,扰动后约13 s 上海测点电压发生阶跃突增；而“两华”同步电网失去有功负荷后进入衰减振荡过程。由图4 可以看出,所提FRRS 方法与传统RLRS 方法均能有效抑制背景噪声与高强度随机脉冲噪声。尤其从图4（a）与（b）的局部放大区域可以看出,快速变化数据的FRRS 降噪畸变明显小于传统RLRS 降噪畸变,所以FRRS 方法比传统RLRS 方法能够更加准确地保留扰动阶跃突变信息与扰动后快速斜坡变化信息。

图4 实测数据降噪结果Fig.4 Denoising results of field measurement data

3.4 计算复杂度分析

传统RLRS 方法与本文所提FRRS 方法的计算复杂度分析结果如表3 所示。传统RLRS 方法每次迭代需要进行5 次矩阵相乘与1 次矩阵求逆运算；而所提FRRS 方法每个滑动窗口仅需要1 次矩阵相乘运算,若存在离群数据,每组离群数据局部平滑畸变快速修复过程也仅需要1 次用于抑制离群数据平滑畸变的矩阵乘法运算并重新平滑离群数据局部N点数据。因此,所提FRRS 方法计算复杂度远小于传统RLRS 方法。2 种算法耗时测试结果如图5 所示,测试环境为MATLAB 2016b,处理器为Intel Core i7-7700,降噪参数仍设置为N=23,N0=5,M=2,数据时间长度从2 h 到24 h,取100 次重复测试结果的均值。所提FRRS 方法的降噪耗时明显低于传统RLRS 方法,降噪速度提高了约20 倍,且随着数据长度增大,降噪速度提升倍数没有明显变化,进一步验证了所提FRRS 方法在计算效率方面的优势。

表3 计算复杂度分析结果Table 3 Analysis results of computational complexity

图5 2 种算法耗时测试结果Fig.5 Test results of time consuming for two methods

4 结语

针对电力系统测量数据受强噪声污染问题,本文提出了一种FRRS 降噪方法,相比传统RLRS 方法,所提方法不仅保持了良好的背景噪声与高强度随机脉冲噪声降噪能力,还能够有效保留扰动阶跃突变原始信息。

由于所述FRRS 方法实现了全局快速回归平滑、随机脉冲噪声自适应检测、局部平滑畸变快速修复,极大降低了计算复杂性,具有在线降噪应用潜力。需要指出,所提FRRS 方法的零权拟合长度参数选择依赖随机脉冲噪声持续时间统计特征先验信息,随着对电力系统测量数据噪声特征统计规律的不断认知,可逐渐优化零权拟合长度参数选择,或者进一步提出能够自适应调整零权拟合长度参数的FRRS 改进方法。