基于RBF和TFT的自适应宽频振荡监测方法

项 颂, 苏 鹏, 吴晓丹, 王 杨, 蒋小龙

(1. 国网内蒙古东部电力有限公司, 内蒙古 呼和浩特 010010; 2. 四川大学电气工程学院, 四川 成都 610065)

1 引言

近年来,以风电和光伏为代表的大量可再生能源通过电力电子设备接入电网,逐步形成了以高比例可再生能源和高比例电力电子设备为特征的“双高”电力系统[1,2]。而电力电子设备带来的弱抗扰、低惯性以及强耦合等特点已严重影响了电力系统的安全稳定运行,导致宽频振荡事故频发,如美国德州电网的双馈风电场发生了多起20 Hz次同步振荡(SubSynchronous Oscillation, SSO)事故[3];我国河北沽源地区出现了多次3～10 Hz次同步振荡现象[4];法国-西班牙联网工程发生了1.4 kHz的振荡事故;我国新疆哈密地区直驱风电机组发生了20～80 Hz的次/超同步振荡事故[5];云南鲁西柔直工程发生了1.27 kHz的高频振荡事故。若振荡事故没有得到有效抑制,可能造成设备损毁、可再生能源机组脱网,甚至传递至用户侧造成大范围的停电事故。随着中国“双碳”目标的提出,可再生能源渗透率将进一步提高,电网抵御宽频振荡的能力将进一步受到挑战。

目前,国内外针对电力系统宽频振荡的分析方法以建模分析为主,如阻抗建模法和状态空间法。尽管该类方法具有坚实的理论基础,但其分析精度在实际应用中极大地受制于系统模型的准确性。基于监测数据,分析人员可以从全局的角度进行振荡风险评估、振荡事故反演、振荡源定位等[6-8]。文献[9]提出,在系统模型不明确的情况下,构建宽频振荡的广域监测系统是分析宽频振荡产生机理与传播特性的有效手段。同时,构建宽频振荡监测系统也是实现振荡实时预警、控制、保护的硬件基础。

当前的广域测量系统以相量测量单元(Phasor Measurement Unit, PMU)为基础,只关注工频分量[9],难以满足宽频振荡监测需求,为此,可以利用子站PMU的波形记录功能或改进广域监测系统(Wide-Area Measurement System, WAMS)动态监测算法将同步测量能力扩大到宽频范围。目前宽频振荡监测的常用方法包括傅里叶变换算法、小波变换算法、希尔伯特黄变换算法(Hilbert-Huang Transform, HHT)、卡尔曼滤波算法、Prony算法等。文献[10-12]均采用加窗的插值傅里叶变换(Interpolated Discrete Fourier Transform, IpDFT)结合滤波器进行宽频振荡参数估计,对稳态信号的估计比较准确,但是难以处理非平稳信号。文献[13,14]采用了小波算法处理时变信号,但是在遇到多频率分量时存在分辨不高的问题。文献[15,16]使用了HHT算法,但是并没有完全解决端点效应及模态混叠的问题。文献[17-19]利用卡尔曼滤波算法对谐波进行跟踪,但是所需参数较多,计算量较大。文献[20,21]将形态滤波和小波算法及Prony算法结合,提高了一定的抗噪性能,但没有解决模型阶数难以确定及计算量较大的问题。除此之外,上述算法在处理强时变特性的信号时均会出现较大误差,甚至会完全失效,因而无法直接应用于复杂多变的振荡信号中。

泰勒傅里叶变换(Taylor-Fourier Transform, TFT)算法因采用动态相量模型而具有更为优越的动态性能[22],因此本文将TFT算法用于宽频振荡监测中,但TFT算法易受到噪声干扰,常需要增加数据窗长[23]或者进行滤波预处理[24]。经作者大量测试得知,常规滤波算法难以解决TFT算法的抗噪问题或以大幅提升计算时间为代价,部分滤波算法甚至会降低TFT算法的精度,因此增加窗长是更为有效的抗噪方案。而TFT算法在噪声较小时并不需要较长的窗长,故本文提出根据噪声强度自适应调整窗长的方案。径向基(Radial Basis Function, RBF)神经网络具有较强输入和输出的映射能力,可以高精度逼近任意非线性函数,且有拓扑结构紧凑、收敛速度快等优点,因此本文将其应用于噪声强度估计。最后,本文通过大量仿真测试及河北沽源和新疆哈密的实测数据进行了验证,结果表明文中方法具有较好的可行性和准确性。

2 泰勒傅里叶变换

电力系统宽频振荡信号的时域模型为:

(1)

式中,Am(t)、fm、φm(t)分别为第m个信号分量的峰值、频率和初相位;N(t)为系统的噪声;M为信号分量总数;t为时间。

以频率fs对x(t)进行采样,经采样后的离散相量模型可表示为:

(2)

pm(n)=Am(nTs)ejφm(nTs)

(3)

式中,ωm为信号的角频率;pm(n)为第m个信号分量的相量,n=-Nw,…,0,…,Nw;Ts为采样周期。

将相量pm(n)在零点进行k阶泰勒级数展开:

(4)

式中,k为泰勒展开的阶数。

根据欧拉公式可以将式(2)写成矩阵形式:

(5)

式中,X为宽频振荡信号的采样序列;B为泰勒傅里叶的基向量矩阵;P为相量的各阶导数及其共轭项所组成的列向量,表达如下:

X=[x(-Nw) …x(0) …x(Nw)]T

(6)

P=[ph1…phm…phM]T

(7)

B=[B1…Bm…BM]

(8)

phm=(pm(k)(0),…,pm(0),pm(0)*,…,pm(k)(0)*)

(9)

(10)

式中,*为共轭运算符。

根据矩阵理论,利用最小二乘法便可求解P。

(11)

式中,H为Hermitian运算符;W为权重对角矩阵,用于降低高频干扰的影响。宽频振荡信号的峰值、相位和频率计算为:

(12)

式中,fm0为第m个分量的频率初始值;abs(·)为取模运算符;angle(·)为求角运算符。

3 TFT算法的影响因素

3.1 初始频率

由基向量矩阵B的构成可知,TFT算法极其依赖初始代入频率,错判频率个数甚至会导致TFT算法直接失效。相反,如果代入频率接近实际频率,TFT算法的辨识精度就会得到保证。因此,本文采用了离散傅里叶变换(Discrete Fourier Transform, DFT)算法来确定初始频率,为了降低频谱泄露以及栅栏效应的影响,使用了窗函数及插值算法,窗函数选择了通带平稳、阻带衰减快的Hanning窗,插值算法选择了三谱线插值算法[25]。

考虑到加窗插值算法无法完全消除栅栏效应和频谱泄露的影响,初始频率的准确性仍可能较低,因此将由式(9)得到的频率回代到式(10)中迭代更新B矩阵并重新用式(9)计算相应参数。理论上来说,迭代的次数越多,TFT算法越精确,但是计算效率也会相应降低。

3.2 窗口长度

TFT算法在计算向量P时所采用的加权最小二乘法在窗口较长时计算量大、运算时间长,在实际应用时需根据具体情况尽可能地缩短所需数据窗长。然而向量P中包含了较多的高阶导数项,高阶导数项对噪声比较敏感,在噪声较大时需要足够多的数据才能保证其准确性,因此需要同时考虑噪声强度及计算量来确定TFT算法的数据窗长,故提出了基于噪声强度的窗长自适应方法,步骤如下:

(1)利用RBF神经网络算法拟合信号。

(2)将采样波形与拟合波形做序列减法,得到噪声序列,然后估计信噪比。

(3)根据信噪比计算数据窗长。

3.2.1 拟合去噪

人工神经网络相比传统拟合算法具有对系统模型的低依赖性、强大自适应能力以及容错性高等优点,且非常适合处理复杂的非线性问题,可以逼近任意函数,因而本文采用神经网络算法对信号进行拟合。在各类神经网络算法中,RBF神经网络具有结构简单、速度快且函数拟合精度高的优点,非常适用于拟合去噪[26]。

RBF神经网络是一种三层前向网络,由输入层、隐含层及输出层组成,结构如图1所示。

图1 RBF神经网络结构Fig.1 Structure of RBF neural network

图1中,X为输入层;xi为输入向量;R为隐含层;ri为神经元;wj为权重向量;F为输出层空间。RBF神经网络采用径向基函数作为激活函数,从而直接将输入向量空间映射到隐含层,因此隐含层是非线性的,输出层则为线性的。

(13)

式中,f(x)为输出层的输出,由隐含层神经元输出的加权求和可得;φ(‖x-xj‖)为隐含层的激活函数,多为高斯函数,如下所示:

(14)

式中,σ为方差;x0为径向基函数的中心点。

令输出层的输出为输入向量的期望响应,记为:

FT=[f1f2…fu]

(15)

同时,令:

φij=φ(‖xi-xj‖)

WT=[w1w2…wv]

(16)

则式(13)可以改写为:

(17)

即:

ΦW=F

(18)

如果Φ可逆,则可得:

W=Φ-1F

(19)

本文采用正则化RBF网络,其隐含层节点数与输入样本数相同,且各输入样本均为径向基函数的中心点,因此无需再另外设计算法。方差σ计算为:

(20)

式中,dmax为中心点之间的最大距离;b为隐含层节点数。

由于RBF神经网络会带来运算量的增加,降低实时性,因此用于拟合的窗不宜过长,文中采用1个基波周期的数据用于噪声估计。

3.2.2 噪声强度估计

将原始波形与经RBF神经网络拟合的波形做差即可得到噪声数据,信噪比便可根据式(21)得到,单位为分贝(dB)。

(21)

式中,PS和PN分别为信号和噪声的功率。

3.2.3 确定窗长

信噪比越高,所需数据越少。对于采样率为12 800 Hz的系统,通过测试发现,当信噪比为30 dB及以上时,TFT算法只需要5个基波周期即可保证高精度的参数估计。当信噪比较低时,窗长应根据实际情况适当增加,但仍不建议超过10个基波周期以保证计算效率和动态性能。因此,窗长Win采用分段处理(遇小数则四舍五入):

(22)

文中所需的窗长一直保持在5～10个基波周期,在信噪比较低时,该方法会在一定程度上降低算法的动态实时性,但是削弱了噪声对参数辨识结果的影响。因此,本文方法兼顾了动态实时性和准确性,在实际应用中,窗长函数可根据应用场景和需求进行适当调整。

3.3 算法实现流程

本文算法的整体流程如图2所示,首先确定泰勒展开的阶数以及TFT的迭代次数,然后对采样序列进行快速估计和噪声强度估计分别得到初始频率和窗口长度,最后基于TFT算法通过多次迭代计算出频率、幅值和相位等参数。

图2 算法流程图Fig.2 Flowchart of algorithm

4 仿真结果及分析

4.1 性能评价标准

本文采用IEEE Std. C37.118.1-2011[27]标准中以评估相量测量算法性能的总体相量误差(Total Vector Error, TVE)和频率误差(Frequency Error, FE)作为评价标准,TVE和FE的定义如下:

(23)

(24)

4.2 仿真信号测试

以次同步振荡为例(基波存在频率偏移)验证本文算法的效果,仿真信号如式(25)所示,并加入了30 dB的高斯白噪声,采样率为12 800 Hz,泰勒傅里叶阶数为2阶,共迭代3次,波形如图3所示。

图3 含次同步振荡波形图Fig.3 Waveform plot with subsynchronous oscillation

(25)

30 dB的信噪比对应的窗长为5个基波周期,拟合结果如图4所示,工频分量的总体相量误差百分比和频率误差分别为0.16%和0.003 Hz,SSO的则为1.77%和0.33 Hz。

图4 拟合对比图Fig.4 Comparison between waveforms

为进一步验证本文方法在复杂工况下的效果,本文参考了IEEE Std. C37.118.1—2011标准中的PMU测试标准,设计了如下仿真算例,且各仿真的噪声水平均从15～60 dB变化。

(1) 幅值动态

在真实情况下,SSO发生的过程往往比较迅猛,本例以衰减指数形式模拟振荡时的幅值动态变化,衰减阻尼为1,信号模型为式(26),仿真结果如图5所示。

(26)

(2) 频率斜坡动态

在实际的振荡事故中,SSO频率可能具有时变性,故以斜坡函数为例进行仿真,信号模型为式(27),SSO分量的频率分量叠加了线性项10t,仿真结果如图6所示。

图6 频率斜坡动态下的误差Fig.6 Errors under frequency ramp test

(27)

(3) 频率调制动态

在实际电力系统中,SSO频率的变化并非一定为线性变化,此处对频率进行正弦调制,调制信号见式(28),频率叠加了正弦项-10sin(2π×0.1t),仿真结果如图7所示。

(28)

(4) 多模式宽频振荡

宽频振荡中可能存在多个信号分量,比如会有谐波和间谐波同时存在的情况,而且各信号分量的动态过程都各有不同,此时的工况愈发复杂,对监测算法的要求也更高。本例的仿真信号为式(29),具体参数见表1,仿真结果如图8所示。

表1 信号参数Tab.1 Parameters of signals

图8 多模式宽频振荡下的误差Fig.8 Error under multi-mode broadband oscillation

(29)

4.3 算法性能讨论

(1)从仿真结果中可以看出,本文算法的动态性能优越,几乎可以应对任何动态变化过程,且在噪声较强时仍具有较高准确性。因TFT算法对泰勒展开项进行了二阶截断,所以当振荡动态变化较快时误差也会有所增大。在实际应用时,如果算力允许,可增加泰勒展开项阶数以进一步增强TFT算法的动态性能。

(2)工频分量的相量误差和频率误差远小于其他信号分量,这是因为工频分量的幅值较大,受噪声影响较小。根据文献[27]中对工频相量的测量标准(最大总体相量误差为1%,最大的频率误差为0.01 Hz),由图5～图8可知,本文算法满足要求。

(3)TFT算法极其依赖初始频率,而振荡信号在振荡初期可能因噪声较大而无法从插值傅里叶算法中识别出来,这有可能导致TFT算法崩溃,因此在使用本算法时可对采样信号进行预处理。

(4)数据窗长对TFT算法也具有较大的影响,在信噪比较低时,若采样率较低,可不限于5个基波周期和10个基波周期的最低限制和最高限制。若需提升动态实时性,也可牺牲准确性降低窗长自适应范围,具体可根据需要做取舍。

(5)需要指出的是,当测量低频振荡时,本文方法与其他基于傅里叶变换的方法相似[10],窗长须有所增加才能保证其准确性。本文当前采用的窗长方案仅适合分析5 Hz及以上的信号。经试验,若需分析更为低频的信号,窗长需要适当增加,如15个基波周期才能达到较高的精确性。而高频振荡并不会影响本文算法的准确度。由于篇幅限制,分析结果不再展示。

(6)本文的测试环境为Matlab R2020b,硬件配置AMD Ryzen 5-4500U CPU@3.00 GHz,运行内存为8 GB DDR4。当信号包含少量谐波或间谐波信号时,本文算法的计算效率较高,比如信号中仅包含次同步振荡分量时,算法运行时间约为0.02 s。但是当信号存在大量的谐波或间谐波时,如本文多模式宽频振荡,此时运行时间约为0.8 s。其主要原因在于信号较多时,TFT算法中的矩阵B的维度会大幅增加,所以导致求解向量P会十分耗时,且RBF算法的拟合时间约为0.01 s。考虑到本文所用的测试环境计算能力较弱,在实际应用时,可选择高性能硬件处理器进一步减少计算时间[24]。

4.4 算法对比

为兼顾辨识算法的稳定性和准确性,目前工程实际应用中多采用基于傅里叶变换的算法实现宽频量测[10]。文献[10]采用了汉宁窗和双谱线插值算法解决传统傅里叶变换中易出现的频谱泄露和栅栏效应问题。以频率调制动态为例,傅里叶算法窗长固定为0.2 s,仿真设置与式(28)一致,对比结果如图9所示。

图9 算法对比Fig.9 Comparison with other algorithms

从图9对比结果中可以看出,在辨识频率非线性变化的次同步分量时,傅里叶算法辨识结果的误差较大,这是因为傅里叶算法会对窗口内的特征参数进行“平均”,所以该算法更适用于线性变化的信号。而对非线性信号的适应性极差,在辨识稳态工频分量时,其辨识结果也不如本文算法,主要是因为傅里叶变换受噪声和数据窗长的影响更大,而且还存在频谱泄露及栅栏效应等问题。

为展示RBF神经网络的拟合效果,本文将局部加权回归平滑算法(Locally Weighted Regression Smoothing, LWRS)、切比雪夫Ⅱ型(ChebyshevII)数字滤波器以及小波算法与RBF神经网络进行比较,信号模型为:

(30)

其中,仿真信号添加了15 dB的高斯白噪声。另外为方便观察,采样率降低为1 600 Hz。LWRS算法的拟合阶数为2阶、窗口长度为15;切比雪夫Ⅱ型滤波器的通带频率为80 Hz、通带波纹为0.1 dB、阻带频率为81 Hz、阻带衰减为80 dB;小波算法采用了dmey小波,分解层数设为3。拟合误差对比如图10所示。

从拟合结果和误差对比图可以看出,RBF神经网络的拟合误差更小,信噪比估算更准确,而且没有复杂的参数设置过程。综上,在实际应用中,本文算法具有一定的优越性和实用性。

5 实测数据验证

本节分别以河北沽源和新疆哈密地区真实发生的次同步及次/超同步振荡案例对本文方法进行验证。

沽源和哈密系统实测的电流分别如图11和图12所示,辨识结果分别如图13和图14所示。根据RBF拟合结果可知,沽源及哈密的噪声强度分别约为30 dB和55 dB,但由于沽源与哈密的故障录波仪的采样率较低,分别只有3 200 Hz和1 200 Hz,故相应增大了数据窗长(10周期)。

图12 哈密某场站电流波形Fig.12 Current waveform of a wind farm in Hami

从辨识结果中能够得知,沽源只发生了次同步振荡事故,而哈密则同时存在次/超同步分量,而且从频率和幅值变化的过程可以判断出这是一对同源的次超同步分量。根据事故后的分析可知,沽源和哈密振荡事故的发生原理不同,前者为双馈风机与固定串补之间的交互作用,而后者为永磁直驱风机与弱电网之间的交互作用,但无论是哪种振荡事故,短时间内振荡的幅值都会迅速增大,进而可能危害到系统的安全稳定运行。

6 结论

本文提出了一种基于波形数据的宽频振荡监测方法。首先,利用DFT算法得到初步的频率参数,并利用RBF神经网络进行拟合去噪,从而得到噪声强度并自适应确定数据窗长,根据所得窗长和初始频率利用TFT算法迭代提取出相关振荡参数。文中通过动态仿真测试和实际的振荡事故案例进行了验证,从测试结果可以看出本文采用的RBF神经网络算法具有较好的拟合效果,且窗长的自适应处理可以保证TFT算法的计算效率和实时性。从参数辨识结果看本文方法的动态性能优越,在各种工况下均有较强的适应性。

本文方法不局限于振荡监测,还可根据实际需求拓展至其他领域,如谐波检测、电能质量等,并在未来实现基于广域监测系统的宽频监测、稳定分析、风险预警等。RBF神经网络已广泛应用于各工程领域中[28,29],本文所采用的基于RBF神经网络的拟合方法所带来的计算负担完全在可接受范围内,但在振荡分量较多且噪声较大时,TFT算法的计算效率会受到影响,因此在实时监测方面的应用仍需进一步研究。