CN109085426A - 一种电力系统工频分量快速提取方法 - Google Patents

Abstract

本发明涉及一种电力系统工频分量快速提取方法，其技术特点在于：包括以下步骤：步骤1、采集电力系统工频分量的原始信号；步骤2、进行低通滤波处理；步骤3、得到采样值x，步骤4、构造样本矩阵R；步骤5、确定样本矩阵R的有效秩p；步骤6、求取Prony算法对应AR模型参数。本发明能够在半个周波(10ms)内(理论上可以时间更短)准确提取电力系统故障工频分量。本发明能够在半个周波(10ms)内(理论上可以时间更短)准确提取电力系统故障工频分量。

Description

一种电力系统工频分量快速提取方法

技术领域

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

背景技术

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

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

发明内容

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

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

一种电力系统工频分量快速提取方法，包括以下步骤：

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

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

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

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

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

步骤6、利用基于拉格朗日乘子法函数求偏导法求取Prony算法对应AR模型参数α₁,α₂,…α_p；

步骤7、求解差分方程1+α₁z^-1+α₂z^-2+…+α_pz^-p＝0的根z_i(i＝1,2,…,p)；

步骤8、构造矩阵Z及

矩阵按以下形式构成：

则

取

步骤9、计算

步骤10、计算频率f_i＝Im(ln z_i)f_s/2π，1≤i≤p

步骤11、计算衰减因子α_i＝Re(ln z_i)f_s，1≤i≤p

步骤12、计算相位θ_i＝angle(z_i)，1≤i≤p

步骤13、计算幅值A_i＝|b_i|，1≤i≤p

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

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

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

(2)定义样本函数其中0≤i≤p_e；0≤j≤p_e，进而构建样本矩阵R：

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

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

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

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

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

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

式中：α₀＝1；

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

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

(2)利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

①首先找到对应于工频分量的特征值方程的根如公式5所示：

由公式3可知，它们还满足如下如公式6所示条件：

②利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

目标函数为公式7：

拉格朗日函数为公式8：

式中：λ₁、λ₂为拉格朗日乘子；

要使目标函数公式7取最小值，必须拉格朗日函数对所有变量及乘子的偏导数为0，即有公式9：

故可得下面一组方程公式10:

定义公式11：其中0≤i≤p；0≤j≤p

则公式10和公式7可以合并写成如公式12所示形式：

③用上述公式12求解得到AR模型参数α_m(1≤m≤p)，从而得到各分量幅值和相位。

而且，所述步骤6第(2)步的第③步的具体方法为：

分析方程可得，λ₁与λ₂互为共轭，又z₁与z₂互为共轭，故有公式13:

又公式12左边系数矩阵最后两行对应方程为公式14：

上述公式14等效于如下公式15所示方程组：

故可用上述公式15替代公式14，改进后的矩阵方程为：

由上式可解得AR模型参数α_m(1≤m≤p)。

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

1、本发明公开了一种新型电力系统工频分量快速提取方法，传统算法受故障信号中噪声及谐波分量的影响，难以在短时间窗下提取出工频分量。本发明充分利用故障信号中包含有稳定工频分量这一特征，将传统AR模型参数求取的最小二乘问题变为有等式约束条件的最小二乘问题，提出基于拉格朗日乘子法的新型Prony算法，该算法有效提高了计算精度和效率，仿真结果表明该算法能够在半个周波(10ms)内准确提取电力系统故障工频分量。

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

附图说明

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

具体实施方式

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

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

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

在电力系统中，特别是故障后的暂态信号，除了含有工频分量外，还有大量整次谐波、非整次谐波和衰减直流分量，故在本实施例中，假设原始信号下式所示：

x(t)＝360sin(2π×34t)e^-3t+500sin(2π×50t+π/3)+150sin(2π×150t)e^-20t

+20sin(2π×1785t+π/6)e^-600t+10sin(2π×2000t+π/7)e^-500t+210e^-25t+δ₁

再加入信噪比为60的高斯白噪声。

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

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

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

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

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

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

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

(2)定义样本函数其中0≤i≤p_e；0≤j≤p_e进而构建样本矩阵R：

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

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

在本实施例中，用SVD法确定R的有效秩为p＝9；

步骤6、利用基于拉格朗日乘子法函数求偏导法求取Prony算法对应AR模型参数α₁,α₂,…α₉；

在本实施例中，利用基于拉格朗日乘子法函数求偏导法求取Prony算法对应AR模型参数α₁,α₂,…α₉；

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

充分利用故障信号中包含有稳定工频分量这一特征，将传统AR模型参数求取的最小二乘问题变为有等式约束条件的最小二乘问题，即拉格朗日乘子法求偏导法，其具体步骤包括：

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

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

上式是一个常系数线性差分方程的齐次解，有：

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

则式(2)对应特征多项式为：

式中：α₀＝1；

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

式中，ε(i)，p≤i≤N-1表示拟合误差；

(2)利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

①首先找到对应于工频分量的特征值方程的根：

由式(3)可知，它们还满足下面条件：

②利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

目标函数为：

拉格朗日函数为：

式中：λ₁、λ₂为拉格朗日乘子。

要使目标函数式(7)取最小值，必须拉格朗日函数对所有变量及乘子的偏导数为0，即有

故可得下面一组方程:

定义：

则式(10)和式(7)可以合并写成如下形式：

在本实施例中，

③用式(12)求解得到AR模型参数α_m(1≤m≤p)，从而得到各分量幅值和相位。(具体推导略)：

在本实施例中，求解上述方程组得α₁,α₂,…α₉；

在本实施例中，采用新方法求解AR参数存在问题及解决措施为：

用式(12)求解AR模型参数时，存在如下两个问题：

a)上式左边的系数矩阵的元素数量级一般相差悬殊，即和与rr(i,j)(0≤i≤9,0≤j≤9)数量级一般相差悬殊，这样导致方程病态，此时方程的解误差非常大。

b)系数矩阵不是实数矩阵，求解比较困难。

为了克服上面的问题，这里采用了以下两种解决办法：

a)采取平衡措施：先求出的数量级，不妨设为10^k，然后将和变换为和相应的将λ₁和λ₂代换为λ₁/10^k和λ₂/10^k，变换后原矩阵方程仍成立。这样就解决了方程病态的问题。

b)分析方程可得，λ₁与λ₂互为共轭。又z₁与z₂互为共轭，故有:

又式(12)左边系数矩阵最后两行对应方程为

式(14)等效于下面方程组：

故可用式(15)替代式(14)，改进后的矩阵方程为

由上式可解得AR模型参数α₁,α₂,…α₉，这样做保证了分析结果中工频分量的计算精度；

步骤7、求解差分方程1+α₁z^-1+α₂z^-2+…+α_pz^-p＝0的跟z_i(i＝1,2,…,p)

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

步骤8、构造矩阵Z及

矩阵按以下形式构成：

其中

则

在本实施例中，

矩阵按以下形式构成：

其中

取

步骤9、计算

步骤10、计算频率f_i＝Im(ln z_i)f_s/2π，1≤i≤9

步骤11、计算衰减因子α_i＝Re(ln z_i)f_s，1≤i≤9

步骤12、计算相位θ_i＝angle(z_i)，1≤i≤9

步骤13、计算幅值A_i＝|b_i|，1≤i≤9

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

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

表1工频分量分析结果

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

Claims (4)

1.一种电力系统工频分量快速提取方法，其特征在于：包括以下步骤：

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

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

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

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

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

步骤6、利用基于拉格朗日乘子法函数求偏导法求取Prony算法对应AR模型参数α₁,α₂,…α_p；

步骤7、求解特征值方程1+α₁z^-1+α₂z^-2+…+α_pz^-p＝0的根z_i(i＝1,2,…,p)；

步骤8、构造矩阵Z及

矩阵按以下形式构成：

其中

则

取

步骤9、计算

步骤10、计算频率f_i＝Im(lnz_i)fs/2π，1≤i≤p

步骤11、计算衰减因子α_i＝Re(lnz_i)f_s，1≤i≤p

步骤12、计算相位θ_i＝angle(z_i)，1≤i≤p

步骤13、计算幅值A_i＝|b_i|，1≤i≤p

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

2.根据权利要求1所述的一种电力系统工频分量快速提取方法，其特征在于：所述步骤4的具体步骤包括：

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

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

3.根据权利要求1所述的一种电力系统工频分量快速提取方法，其特征在于：所述步骤6的具体步骤包括：

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

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

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

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

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

式中：α₀＝1；

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

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

(2)利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

①首先找到对应于工频分量的特征值方程的根如公式5所示：

由公式3可知，它们还满足如下如公式6所示条件：

②利用拉格朗日乘子法函数求偏导法求取PRONY算法对应AR模型参数：

目标函数为公式7：

拉格朗日函数为公式8：

式中：λ₁、λ₂为拉格朗日乘子；

要使目标函数公式7取最小值，必须拉格朗日函数对所有变量及乘子的偏导数为0，即有公式9：

故可得下面一组方程公式10:

定义公式11：

则公式10和公式7可以合并写成如公式12所示形式：

③用上述公式12求解得到AR模型参数α_m(1≤m≤p)，从而得到各分量幅值和相位。

4.根据权利要求3所述的一种电力系统工频分量快速提取方法，其特征在于：所述步骤6第(2)步的第③步的具体方法为：

分析方程可得，λ₁与λ₂互为共轭，又z₁与z₂互为共轭，故有公式13:

又公式12左边系数矩阵最后两行对应方程为公式14：

上述公式14等效于如下公式15所示方程组：

故可用上述公式15替代公式14，改进后的矩阵方程为：

由上式可解得AR模型参数α₁,α₂,…α_p。