CN101447676A - 一种电力系统低频振荡分析方法 - Google Patents

Abstract

本发明提出了一种电力系统低频振荡分析方法，其使用高性能的数字滤波器对原始数据滤波；用短时傅里叶算法识别滤波后数据的频谱及主导振荡；再使用Prony算法分析振荡模态、识别主导振荡。滤波操作解决了Prony算法易受高频噪声影响的问题。根据Prony算法识别的4种基本参数：振幅、频率、相位、衰减因子，推算各振荡模态对应的其他参数，包括：阻尼比、累计振荡能量；通过计算拟合百分比和拟合信噪比参数，判定当前Prony分析的运算效果，客观标定当前分析模态是否能够准确反映原数据信息。在Prony分析得到的各振荡模态中，不仅仅依据振幅识别主导振荡，而是综合考虑振荡模态的振荡类型和累积振荡能量。

Description

一种电力系统低频振荡分析方法

技术领域

本发明提出了一种电力系统低频振荡分析方法，属于动态电力系统的理论和分析领域。

背景技术

1、低频振荡简介

系统稳定问题，即电力系统在受到各种扰动后是否保持稳定运行，是动态电力系统需要研究解决的核心问题。系统稳定问题包括暂态稳定、动态稳定和静态稳定三个大类。电力系统受小扰动时发电机转子间由于阻尼不足而引起的持续低频振荡属于系统动态稳定分析的范畴。

这种振荡现象的振荡频率一般为0.2～2.5Hz。由于发生低频振荡时，输电线路上的功率值也会发生相应振荡，所以又将其称为功率振荡或机电振荡。

低频振荡有两种振荡模式：系统域间振荡(interarea mode)和区域内振荡(plant mode)。发生系统域间振荡时，参与振荡的机组间电气距离较大，振荡频率较低(如0.2～0.5Hz)；发生区域内振荡时，参与振荡的机组间电气距离较小，振荡频率较高(如1Hz以上)。系统域间振荡比区域内振荡具有更大的危害。

2、我国电网特征及实施低频振荡监测分析的必要性

我国电网覆盖面积大，结构薄弱，各种一次能源的分布和负荷的密度极不均匀，而电源又往往远离负荷中心，单位装机容量分摊到的标准输电线长度比发达国家的少得多。

由于低频振荡或功率振荡通常出现在长距离、重负荷输电线上，在采用现代快速、高顶值倍数励磁系统的条件下更容易产生，所以我国电网上发生低频振荡的可能性很大，电力系统的安全稳定运行正面临巨大的挑战。

从国内外已经发生的几次低频振荡事故来看，这种事故对电网危害严重，大大制约了电网的输电能力。通过运行经验的长期积累，发现产生域间低频振荡时，系统网络中各节点的相角摆动有一个由弱到强的过程，即为事故发展的过程。通过对低频振荡进行有效监测，可以在事故发生的初期将其捕捉，实时跟踪分析它的发展状况，及时识别振荡中心，切除扰动源，从而将低频振荡对电网的不良影响降到最低。但是对于区域内发生的低频振荡，由于其往往具有起振快，持续短，平息快的特点，所以这种振荡现象仅能被监测到，而来不及采取应对措施。对于这种振荡，应当对实时数据进行记录，然后分析识别数据中的振荡信息，并以提取的信息为基础，改善系统阻尼特性，最终实现对低频振荡现象的抑制。

3、现有的振荡模式分析方法

传统的低频振荡模式分析方法是：1)根据电网结构和电力元件参数建立微分方程组；2)在系统不同稳定点对微分方程进行线性化并形成状态空间表达式；3)计算特征矩阵的特征根，并由特征根计算系统的低频振荡模式。然而随着电力系统规模的日益扩大，该方法由于“维数灾”的缘故，在实际应用中面临越来越多的困难。

目前较为普遍使用的振荡模式辨识方法是改进Prony法。这种方法能够有效识别各振荡模态的相关参数，但是由于其对噪声的抑制能力较差，使得分析结果的精确程度受到了置疑。

发明内容

本发明提出了一种电力系统低频振荡分析方法，其特征在于包括以下步骤：

1)使用数字滤波器对载入的电力系统有功功率数据进行滤波，可以选择“从前往后滤波方式”或“从后往前滤波方式”，以避开原始数据中较大的无意义尖峰值，达到更好的滤波效果，所述带通滤波器提取原始数据中振荡频率位于0.2～2.5Hz之间的振荡分量；

2)手动截取滤波后数据曲线中的一段进行后续分析；

3)使用短时傅里叶分析算法分析被截取数据段的幅频特性，识别主导振荡，标志主导振荡的振幅和频率；

4)使用改进的Prony分析算法识别被截取数据段中的各个振荡模态，识别每个振荡模态的参数，包括：振幅、频率、相角、衰减因子、阻尼比和累计振荡能量。根据“振幅”和“累计振荡能量”的数值，按照由大到小的顺序对各振荡模态排序，根据识别的参数信息判定各模态的振荡类型，依据累计振荡能量和振荡类型标志主导振荡；

5)将Prony分析得到的各个振荡模态叠加，拟合滤波后的有功功率数据，根据拟合信噪比和拟合百分比判定拟合效果，判定Prony分析结果是否可信。

本发明的方法其特征还在于所述振荡类型包括域间振荡、区域内振荡和高频分量。

离线分析算法流程如附图1所示。

本发明具有以下优点：

1、低频振荡离线分析算法中使用带通滤波器滤除原始信号中的高频噪声后，再进行Prony分析，解决了Prony算法易受高频噪声影响的问题；

2、短时傅里叶分析算法能够准确识别主导振荡频率，Prony算法能够识别主要的振荡模态及它们的特性参数，用这两种方式同时标定主导振荡更为准确；

3、能够利用振荡参数分析Prony算法识别的振荡模态的类型；

4、按振幅及累积振荡能量递减的顺序，对Prony算法识别的振荡模态排序，能够判定各振荡模态在现阶段所占能量比重；

5、应用拟合百分比和拟合信噪比参数客观判定Prony分析是否有效。

附图说明

下面结合附图对本发明进一步说明。

图1是依据本发明的电力系统低频振荡离线分析算法流程；

图2是本发明的方法中的数字滤波器滤波过程示意图；

图3是短时傅里叶分析算法流程图；

图4是改进Prony分析算法流程图。

具体实施方式

以下更具体地介绍本发明的内容。

振荡模式识别不但能够了解单条被监测线路的固有振荡特性，对分析振荡源位置和振荡原因也具有积极意义。

本发明的方法整合传统的低频振荡模式分析方法，通过各算法相互配合，结果相互映证，有效识别低频振荡模式，准确标定主导振荡及其参数。

离线分析算法流程如附图1所示，该方法包括以下步骤：

步骤1：载入告警录波数据。

将在线监测算法获得的打包数据解包，将数据按采样顺序保存在离线分析算法的“原始数据缓冲区”中。缓冲区的长度取决于原始数据采样率，采样速率越高，“原始数据缓冲区”的长度越长。“原始数据缓冲区”共能容纳1分钟的采样数据。

步骤2：应用数字滤波器提取录波数据中振荡频率位于0.2～2.5Hz频率区间的振荡分量。

如附图2所示，使用数字滤波器对“原始数据缓冲区”中的采样值进行滤波的过程，就是使用滤波器窗口由前至后依次截取“原始数据缓冲区”中的一段连续数据，依次得到连续的滤波输出值。由于滤波器窗口具有一定长度，所以并不能对所有原始数据进行滤波。如果滤波器窗口长度为L，L为奇数，那么“原始数据缓冲区”头部和尾部的各(L-1)/2的采样值将不能被滤波。滤波数据结果将依次存入离线分析算法的“滤波后数据缓冲区”。该缓冲区的长度比“原始数据缓冲区”的长度小(L-1)。

步骤3：从滤波后数据缓冲区中截取需要分析的数据段。

所截取数据段的长度将对后续的“短时傅里叶分析算法”和“改进Prony分析算法”产生影响：数据长度越短，由“短时傅里叶分析算法”得到的幅频谱的频率分辨率越低；数据长度越长，执行“改进Prony分析算法”耗时越长；如果数据长度过长，以至超过“短时傅里叶分析算法”的时窗宽度，则仅能对部分所选数据进行频谱分析，同时还可能导致Prony分析算法拟合失败。所以截取数据长度以包含两个完整周波为宜。

步骤4：应用“短时傅里叶分析算法”分析所截取数据段的幅频特性，识别主导振荡，标定主导振荡振幅和频率。

1)基本原理

短时傅里叶分析算法是对离散傅里叶分析算法的发展和提高。应用它在时域也具有分辨率的特性，可以通过不断平移时间窗，连续观察被测数据在不同时刻具有的频谱特征，更可了解任意频率分量延时间轴正向的能量变化情况。信号x[n]的短时傅里叶变换定义为：

$X[n,\mathrm{\λ})=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}x[n+m]w\left[m\right]e^{-\mathrm{j\λm}}---\left(1\right)$

式中w[n]是一个窗序列。在短时傅里叶表示中，一维序列x[n]是单个离散变量的函数，它将转换为一个离散的时间变量n和连续的频率变量λ的二维函数。

2)参数选择

利用“短时傅里叶分析算法”计算连续信号的频谱时，应根据实际问题中对混叠、泄漏、频域分辨率和时域分辨率等要求确定：

(1)信号采样间隔；

信号采样间隔就是信号采样率，该参数不能调节，由录波数据的采样率决定，采样率需满足采样定理。

(2)信号截取窗函数及其参数；

短时傅里叶分析算法又称加窗傅里叶分析算法，是通过时窗截取连续无限长采样数据，再对时窗内数据进行频谱分析的方法。

对原信号进行加窗，进行周期卷积运算后，将对原信号中的尖峰和不连续点起平滑作用。这种谱平滑作用，将使一个频率处的分量泄漏到相邻的另一个分量中去，如果两个相邻频率的相位相异，将会导致频谱曲线谱峰的降低，严重时还可导致相邻谱峰的重叠，即不可分辨。

频谱发生泄漏和分辨率的降低，是对信号加窗后对原信号频谱产生的两种主要影响。分辨率主要受窗函数频谱函数W(e^jω)主瓣宽度的影响，而泄漏的程度则取决于W(e^jω)的主瓣和旁瓣的相对幅度。就窗函数自身而言，在频谱曲线中，主瓣的宽度和相对旁瓣的幅度主要取决于窗的长度和窗的形状(拖尾的大小)。

加窗操作就是用有限长窗w[n]乘以时窗内截取的有限长序列x[n]，从而产生有限长序列v[n]＝w[n]x[n]。它在频域上表现为周期卷积，即：

$V\left(e^{\mathrm{j\ω}}\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}X\left(e^{\mathrm{j\θ}}\right)W\left(e^{j(\mathrm{\ω}-\mathrm{\θ})}\right)\mathrm{d\θ}---\left(2\right)$

对加窗序列v[n]＝w[n]x[n]进行DFT变换为：

$\begin{array}{cc}V\left[k\right]=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}v\left[n\right]e^{-j(2\mathrm{\π}/N)\mathrm{kn}}& k=\mathrm{0,1},...,N-1\end{array}---\left(3\right)$

有限长序列v[n]的DFT对应于v[n]的傅里叶变换的等间隔采样，采样间隔为2π/N。

此处的窗函数为凯赛窗。使用凯塞窗的好处在于：这种窗有两个参数β和L，它们可以用于调节主瓣宽度和相对旁瓣幅度。

主瓣宽度Δ_ml定义为窗函数频谱曲线上，在中央的两个过零焦点之间的对称距离。相对旁瓣高度A_sl定义为以dB计的主瓣幅度与最大旁瓣幅度之比。

在用凯塞窗进行谱分析时，首先根据能够接受的频谱泄漏程度确定A_sl，然后根据下式确定β。相对旁瓣幅度基本上与窗的长度无关，只取决于β。

$\mathrm{\&beta;}=\left\{\begin{array}{cc}0,& A_{\mathrm{sl}}<13.26\\ 0.76609{(A_{\mathrm{sl}}-13.26)}^{0.4}+0.09834(A_{\mathrm{sl}}-13.26),& 13.26<A_{\mathrm{sl}}<60\\ 0.12438(A_{\mathrm{sl}}+6.3),& 60<A_{\mathrm{sl}}<120\end{array}\right.---\left(4\right)$

试验表明，主瓣的宽度与窗的长度成反比，主瓣宽度、相对旁瓣幅度和窗长度之间的折中关系可以用下式表示，在给定主瓣宽度Δ_ml后，使用下式可以估算出窗宽L。

$L\&ap;\frac{24\mathrm{\&pi;}(A_{\mathrm{sl}}+12)}{155{\mathrm{\&Delta;}}_{\mathrm{ml}}}+1---\left(5\right)$

(3)对窗宽补零后的实际窗宽，即进行DFT运算的实际点数；

在不对窗口进行补零时，频率分辨率是频域中谱线间的最小间隔，它等于信号的基波频率f₁，f₁越小则分辨率越高。

若对窗口进行补零，以扩大窗宽，那么频率分辨率就与实际窗口宽度T₁(即实际DFT运算点数)成反比。若截取数据长度不变，则补零宽度越大，分辨率越高。可以考虑极限情况：当T₁为无限长时，则相当于对连续信号进行频谱分析，此时频率分辨率极高，谱线分布极密，所得的频域频谱曲线几乎成为连续曲线。

补零宽度取决于信号分析所需要的频率分辨率。若信号的基波频率为f₁，补零宽度等于截取的实际数据宽度的N倍，那么补零后的实际频率分辨率为f₁/N。

需要注意，对截断序列补零，应当在被截取序列两端对称补零，这样才不会对窗口的时窗中心产生影响。

3)算法流程

如附图3所示。

步骤5：应用“改进Prony算法”分析所截取数据段包含的主要振荡模态，标定各振荡模态的特性参数，包括：振幅、频率、相角、衰减因子、阻尼比、累积振荡能量。

1)基本原理

Prony方法采用的数学模型为一组P个具有任意幅值、相位、频率与衰减因子的指数函数，其离散时间的函数形式为：

$\begin{array}{cc}\hat{x}\left(n\right)=\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{b}_i{z_i}^n,& n=\mathrm{0,1},...,N-1\end{array}---\left(6\right)$

拟合输出

是x(n)的近似值。式(6)中，b和z假定为复数，即：

b_i＝A_iexp(jθ_i)                   (7)

z_i＝exp[(α_i+j2πf_i)Δt]            (8)

式中，A为幅值；θ为相位(单位：弧度)；α为衰减因子；f表示振荡频率(单位：赫兹)；Δt代表采样间隔(单位：秒)。将Prony模型与算法应用于低频振荡监测尤为适宜。

2)应用过程

改进Prony方法可以叙述如下。

(1)利用式(9)计算样本函数r(i，j)，

$\begin{array}{cc}r(i,j)=\underset{n=p}{\overset{N-1}{\mathrm{\&Sigma;}}}x(n-j)x^{*}(n-i),& i,j=\mathrm{0,1},...,p\end{array}---\left(9\right)$

并构造扩展阶的矩阵：

$\mathrm{Re}=\left[\begin{array}{cccc}r\left(\mathrm{1,0}\right)& r\left(\mathrm{1,1}\right)& L& r(1,p_e)\\ r\left(\mathrm{2,0}\right)& r\left(\mathrm{2,1}\right)& L& r(2,p_e)\\ M& M& M& M\\ r(p_{e}0,)& r(p_e,1)& L& r(p_e,p_e)\end{array}\right],p_e?p---\left(10\right)$

(2)用SVD—TLS算法确定矩阵Re的有效秩p以及参数a₁，...，a_p的总体最小二乘估计。

(3)求特征多项式(11)的根z₁，...，z_p，

1+a₁z^-1+...+a_pz^-p＝0                (11)

并利用公式(12)计算出n＝p，...，N-1，其中对于n＝0，1，...，p-1，有 $\hat{x}\left(n\right)=x\left(n\right).$

$\begin{array}{cc}\hat{x}\left(n\right)=-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_i\hat{x}(n-i),& n=p,...,N-1\end{array}---\left(12\right)$

(4)利用式(13)～式(15)计算参数b₁，...，b_p。

$b={\left(Z^{H}Z\right)}^{-1}{Z}^H\hat{x}---\left(13\right)$

$Z^{H}Z=\left[\begin{array}{cccc}{\mathrm{\&gamma;}}_{11}& {\mathrm{\&gamma;}}_{12}& L& {\mathrm{\&gamma;}}_{1p}\\ {\mathrm{\&gamma;}}_{21}& {\mathrm{\&gamma;}}_{22}& L& {\mathrm{\&gamma;}}_{2p}\\ M& M& M& M\\ {\mathrm{\&gamma;}}_{p1}& {\mathrm{\&gamma;}}_{p2}& L& {\mathrm{\&gamma;}}_{\mathrm{pp}}\end{array}\right]---\left(14\right)$

$\mathrm{\&gamma;ij}=\frac{{\left({z_i}^{*}{z}_j\right)}^N-1}{\left({z_i}^{*}{z}_j\right)-1}---\left(15\right)$

(5)用下式计算幅值A_i、相位θ_i、频率f_i和衰减因子α_i：

$\left.\begin{array}{c}{A}_i=\left|b_i\right|\\ {\mathrm{\&theta;}}_i=\arctan [\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right)]/\left(2\mathrm{\&pi;\&Delta;t}\right)\\ {\mathrm{\&alpha;}}_i\text{=ln}\left|z_i\right|/\mathrm{\&Delta;t}\\ f_i=\arctan [\mathrm{Im}\left(z_i\right)/\mathrm{Re}\left(z_i\right)]/\left(2\mathrm{\&pi;\&Delta;t}\right)\end{array}\right\},i=1,...,p---\left(16\right)$

(6)根据衰减因子和频率用下式求取各振荡模态的阻尼比：

$\mathrm{\&xi;}=\frac{-\frac{\mathrm{\&alpha;}}{f}}{2\mathrm{\&pi;}\sqrt{1+\frac{{\mathrm{\&alpha;}}^2}{4{\mathrm{\&pi;}}^2{f}^2}}}---\left(7\right)$

该式中，ξ是阻尼比，α是衰减因子，f是振荡频率。

(7)求取累积能量

根据振荡模态参数和被分析数据长度得到以各个振荡模态拟合原数据的各条拟合曲线，按与录波数据相同的采样率对各条拟合曲线采样，求取每条曲线上各采样值的平方值并累加，将其作为对应振荡模态的累积能量。

3)算法流程

如附图4所示。

步骤6：根据“频率”和“阻尼比”判断各振荡模态对应的“振荡类型”。

依据表1判断振荡类型：

表1  振荡类型分类列表

步骤7：根据“振幅”及“累积振荡能量”递减的顺序将各振荡模态排序。在这些振荡模态中，若两个模态的“累积振荡能量”相同，则“振幅”较大的模态排位靠前。

步骤8：综合评估“振幅”、“累积振荡能量”和“振荡类型”，得到基于Prony法识别的主导振荡。

识别方法为：在排序后的振荡模态中由前至后获取若干振荡模态，使它们的“累积振荡能量”占总振荡能量的80％以上。在这些振荡模态中剔除“振荡类型”为“直流”或“高频分量”的模态，剩余的即为主导振荡。

步骤9：计算“拟合百分比”和“拟合信噪比”，判断Prony分析所得的振荡模态是否可以客观反映原始数据所包含的振荡信息。

衡量Prony分析结果有两个指标：

(1)拟合信噪比

假设实际数据为y(k)，Prony算法模型输出为

$\begin{array}{cc}\hat{y}\left(k\right),& k=\mathrm{0,1},L,N-1\end{array}---\left(18\right)$

拟合信噪比(Signal/Noise Ratio，SNR)是最常用的指标。其定义为

$\mathrm{SNR}=20\lg \frac{\mathrm{rms}\left[y\left(k\right)\right]}{\mathrm{rms}[y\left(k\right)-\hat{y}\left(k\right)]}---\left(19\right)$

其中rms表示均方根(root mean square)，单位为dB。在物理上，均方根值也称为有效值，它对于每个采样点的计算方法是：先平方、再平均、最后开方。

(2)拟合百分比

其定义为

$\frac{\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}|y\left(k\right)-\hat{y}\left(k\right)|}{\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}|y\left(k\right)-y\left(0\right)|}*100\%---\left(20\right)$

一般认为拟合百分比小于10％以及SNR达到20dB以上时，Prony法分析结果可以被接受，若SNR接近40dB则表明分析结果更加理想。另外，通过比较真实数据曲线和拟合数据曲线，也可以作为一种辅助手段直观判断拟合效果。

此处已经根据特定的示例性实施例对本发明进行了描述。对本领域的技术人员来说在不脱离本发明的范围下进行适当的替换或修改将是显而易见的。示例性的实施例仅仅是例证性的，而不是对本发明的范围的限制，本发明的范围由所附的权利要求定义。

Claims (2)

1、一种电力系统低频振荡分析方法，其特征在于包括以下步骤：

1)使用数字滤波器对载入的电力系统有功功率数据进行滤波，可以选择“从前往后滤波方式”或“从后往前滤波方式”，以避开原始数据中较大的无意义尖峰值，达到更好的滤波效果，所述带通滤波器提取原始数据中振荡频率位于0.2～2.5Hz之间的振荡分量；

2)手动截取滤波后数据曲线中的一段进行后续分析；

3)使用短时傅里叶分析算法分析被截取数据段的幅频特性，识别主导振荡，标志主导振荡的振幅和频率；

4)使用改进的Prony分析算法识别被截取数据段中的各个振荡模态，识别每个振荡模态的参数，包括：振幅、频率、相角、衰减因子、阻尼比和累计振荡能量。根据“振幅”和“累计振荡能量”的数值，按照由大到小的顺序对各振荡模态排序，根据识别的参数信息判定各模态的振荡类型，依据累计振荡能量和振荡类型标志主导振荡；

5)将Prony分析得到的各个振荡模态叠加，拟合滤波后的有功功率数据，根据拟合信噪比和拟合百分比判定拟合效果，判定Prony分析结果是否可信。

2、如权利要求1所述的方法，所述振荡类型包括域间振荡、区域内振荡和高频分量。

Images