一种经补偿电力系统故障工频分量快速提取算法的制作方法

本发明属于电力系统数字信号处理技术领域，尤其是一种经补偿电力系统故障工频分量快速提取算法。

背景技术：

电力系统信号处理主要涉及对工频分量的提取，例如电力系统中系统状态估计、稳定控制、故障分析及继电保护等领域中都需要对工频分量进行计算及提取。尤其是继电保护领域更是要求快速准确提取工频分量，缩短保护动作时间，从而满足现代电力系统对继电保护快速可靠性的要求。

目前在电力系统继电保护中，广泛采用的方法有全波傅氏算法、半波傅氏算法、最小二乘算法、卡尔曼滤波算法等。全波傅氏算法是基于周期函数模型推导出来的。当电力系统发生故障时，故障信号不是周期函数。此时，由于衰减直流分量和分数次谐波的存在，全波傅氏算法计算的工频量幅值和相位有较大的误差。另外，傅氏算法至少需要一个周波(20ms)的采样数据。虽然半波傅氏算法仅需半个周波的采样数据，但是其受衰减直流分量和分数次谐波的影响更大，并且半波傅氏算法无法滤除偶数次谐波，故它计算得到的工频量幅值和相位的误差更大。最小二乘算法对数据窗长的要求非常灵活，滤波性能良好，精度也很高，但是它很难选择算法的噪声模型。最小二乘算法在噪声模型的选取上不具有普遍性：不同的系统中，噪声模型可能不同，即使同一系统下，发生不同故障时，噪声模型的选取也可能不同。卡尔曼滤波算法受随机干扰和测量噪声的影响较小，但是需要复杂的矩阵运算，并且算法的滤波参数选择基于对含噪声信号的统计特征分析，求取困难。因此有必要研究一种新的能够快速准确提取故障特征量的算法，来提高保护的动作特性，以满足电力系统发展的需要。

技术实现要素：

本发明的目的在于提供一种设计合理且能够快速且准确的得到故障后的稳态分量特征的经补偿电力系统故障工频分量快速提取算法。

本发明解决其现实问题是采取以下技术方案实现的：

一种经补偿电力系统故障工频分量快速提取算法，包括以下步骤：

步骤1、通过信息获取或数据采集得到电力系统工频分量的原始信号；

步骤2、对电力系统工频分量的原始信号进行低通滤波处理；

步骤3、若步骤2中滤波后得到的电力系统工频分量的信号是连续信号，则须对其进行采样，得到采样值x，设采样率为fs，采样率满足采样定理；

步骤4、利用步骤3中得到的采样值x构造样本矩阵r；

步骤5、利用svd方法确定样本矩阵r的有效秩p；

步骤6、利用多元线性回归方法求取prony算法对应ar模型参数

α1,α2,…αp；

步骤7、对求取到的ar模型参数进行补偿得到更加精确的ar模型参数；

步骤8、求解特征值方程1+α1z^-1+α2z^-2+…+αpz^-p＝0的根zi(i＝1,2,…,p)；

步骤9、构造矩阵z及

矩阵按以下形式构成：

其中，

则

取

步骤10、计算

步骤11、计算频率fi＝im(lnzi)fs/2π，1≤i≤p

步骤12、计算衰减因子αi＝re(lnzi)fs，1≤i≤p

步骤13、计算相位θi＝angle(zi)，l≤i≤p

步骤14、计算幅值ai＝|bi|，1≤i≤p

步骤15、在求得的所有频率分量中必然包含有工频分量，从而求得工频分量。

而且，所述步骤4的具体步骤包括：

(1)取n个连续采样点x(i),i＝1,2,...,n且n>2p；其中，n为采样点个数，p为信号中含有的频率分量个数所对应的prony算法模型阶数；

(2)定义样本函数其中pe＞＞p；pe＜＜n；0≤i≤pe；0≤j≤pe进而构建样本矩阵r：

而且，所述步骤6的具体步骤包括：

(1)定义传统prony算法测量数据模型为公式1：

式中：n＝0,1,...,n-1；n为采样点数；为第n个采样点的估计值；bl与zl均是复数；p是模型阶数；

公式1是一个常系数线性差分方程的齐次解，公式1对应的常系数线性差分方程为公式2：

式中：p≤n≤n-1；αm(1≤m≤p)为其对应特征多项式的系数；

则公式2对应特征多项式为公式3：

式中：α0＝1；

将问题转化为求公式2的的差分方程的系数，而经推导可得求取该系数实为求ar模型参数，即求解下列矩阵方程公式4：

式中，ε(i)表示拟合误差；

(2)利用多元线性回归方法求取prony算法对应ar模型参数α1,α2,…：αp；

记

x＝(x(p),x(p+1),…,x(n-1))^t；f＝(ε(p),…,ε(n-1))^t；a＝(α1,α2,…,αp)^t。

则公式4可写为：ya＝x+f

由多元线性回归知识可得，a的估计值为

a＝(y^ty)^-1y^tx。

由此计算得到ar模型参数α1,α2,…αp。

而且，所述步骤7的具体步骤方法为：在求出ar模型参数后，在ar参数α1,α2,…αp中找出任意两个参数作为未知数，其余参数作为已知数，重新求取选择的两个未知参数，进而对ar参数α1,α2,…αp进行了补偿，最终确保分析结果中出现工频分量；

而且，所述步骤7的求取未知参数的方法为：

首先找到特征值方程1+α1z^-1+α2z^-2+…+αpz^-p＝0对应于工频分量的的根z1、z2，如下公式5所示：

由公式3可得如下公式6：

上式两边同时除以(z-z1)(z-z2)，得到如下公式7：

对上面左式利用长除法，可建立关于两个未知参数的方程组，进而可以求得两个未知参数：

对(z-z1)(z-z2)展开，可得到如下所示公式8：

(z-z1)(z-z2)＝z²-(z1+z2)z+z1z2

其中z1+z2＝e^j100πt+e^-j100πt＝2cos(100πt)，z1z2＝e^j100πte^-j100πt＝1；

设k＝2cos(100πt)，

则有公式9：(z-z1)(z-z2)＝z²-kz+1

取p＝3为例进行详细说明如下(p取其他值时，方法同下)，

由公式6可得，(z-z1)(z-z2)(z-z3)＝z³+a1z²+a2z+a3，用z²-kz+1去除z³+a1z²+a2z+a3，进行长除法运算，如下所示：

由于z²-kz+1能整除z³+a1z²+a2z+a3，固可得下面方程组公式15：

在ar参数α1,α2,…αp中找出任意两个参数作为未知数，其余参数作为已知数，进而由上面方程组可解得未知数。

本发明的优点和有益效果：

1、本发明公开了一种新型经补偿电力系统故障工频分量快速提取算法，传统算法受故障信号中噪声及谐波分量的影响，难以在短时间窗下提取出工频分量。本发明提出经补偿的prony算法用于电力系统工频分量快速提取，本发明充分利用故障信号中包含有稳定工频分量这一特征，对利用多元线性回归方法求取得到的prony算法对应的ar模型参数进行了有效补偿，该算法有效提高了计算精度和效率，仿真结果表明该算法能够在四分之一个周波(5ms)甚至更短时间内准确提取电力系统故障工频分量。

2、本发明能够提供一种不受整次谐波、非整次谐波和直流分量影响的工频分量提取方法；容许使用任意长度的数据窗(采样率确定时，至少应满足采样定理)，并能解决频率泄露问题；能够在四分之一个周波(5ms)甚至更短时间内准确提取电力系统故障工频分量。

附图说明

图1是本发明的处理流程图。

具体实施方式

以下结合附图对本发明实施例作进一步详述：

一种经补偿电力系统故障工频分量快速提取算法，如图1所示，包括以下步骤：

步骤1、通过信息获取或数据采集得到电力系统工频分量的原始信号；

在电力系统工频分量常用提取算法傅氏算法受衰减直流分量的影响很大，故假设原始信号如式(1)所示，再加入信噪比为60的高斯白噪声；

步骤2、对电力系统工频分量的原始信号进行低通滤波处理；

在本实施例中，低通滤波器的截止频率设为1000hz。

步骤3、若步骤2中滤波后得到的电力系统工频分量的信号是连续信号，则须对其进行采样，得到采样值x，设采样率为fs，采样率满足采样定理；

在本实施例中，采样时选6khz作为采样频率，取5ms数据窗，则取到30个连续点x(i),i＝1,2,...,30。

步骤4、利用步骤3中得到的采样值x构造样本矩阵r；

所述步骤4的具体步骤包括：

(1)取n个连续采样点x(i),i＝1,2,...,n且n>2p；其中，n为采样点个数，p为信号中含有的频率分量个数所对应的prony算法模型阶数；

(2)定义样本函数其中pe＞＞p；pe＜＜n；0≤i≤pe；0≤j≤pe进而构建样本矩阵r：

在本实施例中，取30个连续点(对应5ms的数据窗)x(i),i＝0,1,...,29,构造样本矩阵r，形式如下：

步骤5、利用svd方法确定样本矩阵r的有效秩p(所求有效秩p极为prony算法模型阶数)；

在本实施例中，用svd法确定r的有效秩为p＝3；

步骤6、利用多元线性回归方法求取prony算法对应ar模型参数α1,α2,…αp；

在本实施例中，利用多元线性回归方法求取prony算法对应ar模型参数α1,α2,α3；

所述步骤6的具体方法为：

记

x＝(x(3),x(4),…,x(29))^t；a＝(α1,α2,α3)^t。

由多元线性回归知识可得，a的估计值为

a＝(y^ty)^-1y^tx。

由此计算得到ar模型参数α1,α2,α3。

步骤7、对求取到的ar模型参数进行补偿得到更加精确的ar模型参数；

在本实施例中，由于电力系统发生故障后，故障信号中包含有稳定的工频分量，故我们用prony算法分析故障信号时，分析结果中应出现50hz频率成分。为此，充分利用故障信号中包含有稳定的工频分量的特点，在用前面方法求取ar模型参数后对其进行补偿，使分析结果中始终有工频分量。

所述步骤7的具体方法为：在求出ar模型参数后，在ar参数α1,α2,…α3中找出任意两个参数作为未知数，其余参数作为已知数，重新求取选择的两个未知参数，进而对ar参数α1,α2,…α3进行了补偿，最终确保分析结果中出现工频分量；

所述步骤7的求取未知参数的方法为：

首先找到特征值方程1+α1z^-1+α2z^-2+α3z^-3＝0对应于工频分量的的根z1、z2，如下公式5所示：

其中

由特征多项式可得如下方程

上式两边同时除以(z-z1)(z-z2)，得到

对上面左式利用长除法，可建立关于两个未知参数的方程组，进而可以求得两个未知参数。具体如下：

对(z-z1)(z-z2)展开，

(z-z1)(z-z2)＝z²-(z1+z2)z+z1z2(13)

其中z1+z2＝e^j100πt+e^-j100πt＝2cos(100πt)，z1z2＝e^j100πte^-j100πt＝1。

设k＝2cos(100πt)，

则有(z-z1)(z-z2)＝z²-kz+1(14)

由式(11)可得，(z-z1)(z-z2)(z-z3)＝z³+α1z²+α2z+α3，用z²-kz+1去除z³+α1z²+α2z+α3，进行长除法运算，如下所示：

由于z²-kz+1能整除z³+α1z²+α2z+α3，固可得下面方程组

在ar参数α1,α2,α3中找出任意两个参数作为未知数，其余参数作为已知数，进而由上面方程组可解得未知数。这样我们就对原ar模型参数进行了补偿。用补偿后的ar模型参数进行计算，能够保证分析结果中肯定包含频率为50hz、衰减因子为0的分量，即工频分量。

步骤8、求解差分方程1+α1z^-1+α2z^-2+…αpz^-p＝0的根zi(i＝1,2,…,p)；

在本实施例中，求解差分方程1+α1z^-1+α2z^-2+α3z^-3＝0的根zi(i＝1,2,…,3)

步骤9、构造矩阵z及

在本实施例中，

矩阵按以下形式构成：

其中，在本实施例中，

其中，

则

在本实施例中，

取

步骤10、计算

步骤11、计算频率fi＝im(lnzi)fs/2π，1≤i≤p

步骤12、计算衰减因子αi＝re(lnzi)fs，1≤i≤p

步骤13、计算相位θi＝angle(zi)，1≤i≤p

步骤14、计算幅值ai＝|bi|，1≤i≤p

步骤15、在求得的所有频率分量中必然包含有工频分量，从而求得工频分量。

本实施例中，计算所得工频分量具体如下：

表1是信噪比为60的高斯白噪声添加，数据窗为10ms时，新型prony算法分析结果表。

需要强调的是，本发明所述实施例是说明性的，而不是限定性的，因此本发明包括并不限于具体实施方式中所述实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。