CN106155981A - 一种次同步振荡参数检测方法 - Google Patents

Abstract

本发明涉及一种次同步振荡参数检测方法，所述方法包括：对电流信号进行截断，获取长度为N的采样数据；建立长度为N的四项Rife‑Vincent(III)窗函数，并根据所述四项Rife‑Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值；在次同步频段内筛选所述离散频点的幅值的极大值，并利用所述离散频点的幅值的极大值对应的谱线及谱线相邻的两根谱线采用插值修正算法对所述电流信号的幅值、频率和相位进行修正；本发明提供的方法能够利用四项Rife‑Vincent(III)窗三谱线插值FFT的方式有效的抑制频谱泄露。

Description

一种次同步振荡参数检测方法

技术领域

本发明涉及电力系统运行与控制领域，具体涉及一种次同步振荡参数检测方法。

背景技术

为提高长距离输电线路的输送能力，电网结构广泛采用串联补偿和高压直流输电技术，这使得电力系统面临着较大的次同步振荡的风险。随着柔性交流输电系统的逐步应用以及新能源的大规模并网运行，电力系统的次同步振荡风险愈加普遍和严重，严重威胁电力系统的安全和稳定运行。对次同步振荡模态参数的准确检测和分析有利于掌握电力系统的实时运行状况，同时也是将各种抑制和控制次同步振荡的方法付诸实施的基础。

目前，次同步振荡参数识别最常用的方式是快速傅里叶变换(FFT)，然而在非同步采样或对信号进行非整数周期截断的情况下，FFT存在栅栏效应和频谱泄露现象，导致次同步振荡参数(幅值、频率、相位)计算结果不准确。利用性能优良的窗函数可以抑制长范围泄露引起的误差，通过插值算法进行修正可以对短范围泄露进行修正。

常用的窗函数有Hanning窗、Blackman窗、Blackman-Harris窗、Nuttall窗等，不同窗函数的频谱旁瓣特性对于频谱泄露的抑制作用各有特点，插值算法中应用广泛且准确度较高的为双谱线插值修正方法，双谱线插值算法利用峰值频点附近幅值最大和次最大的两条谱线，由于该算法没有充分利用频点附近泄露谱线所包含的信息，因而导致修正误差较大。

发明内容

本发明提供一种次同步振荡参数检测方法，其目的是利用四项Rife-Vincent(III)窗三谱线插值FFT的方式有效的抑制频谱泄露。

本发明的目的是采用下述技术方案实现的：

一种次同步振荡参数检测方法，其改进之处在于，包括：

对电流信号进行截断，获取长度为N的采样数据；

建立长度为N的四项Rife-Vincent(III)窗函数，并根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值；

在次同步频段内筛选所述离散频点的幅值的极大值，并利用所述离散频点的幅值的极大值对应的谱线及谱线相邻的两根谱线采用插值修正算法对所述电流信号的幅值、频率和相位进行修正。

优选的，所述四项Rife-Vincent(III)窗函数w(n)的时域表达式为：

$w\left(n\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{a}_{m}cos(2\mathrm{\π}mn/N)---\left(1\right)$

式(1)中，n∈[0,N-1]，N为所述采样数据的长度，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158。

优选的，所述根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值，包括：

对所述采样数据进行加窗处理，获取加窗序列x_w(n)，公式为：

x_w(n)＝w(n)·x(n) (2)

式(2)中，w(n)为四项Rife-Vincent(III)窗函数，x(n)为所述采样数据；

对所述加窗序列进行离散傅里叶变换，公式为：

式(3)中，X(k)为所述加窗序列的离散傅里叶变换函数，k∈[0,N-1]，N为所述采样数据的长度，A为所述电流信号的幅值，f为所述电流信号的频率，为所述电流信号的相位，W(·)为窗函数的离散傅立叶变换函数，Δf为频率分辨率，其中，所述频率分辨率的公式为：

Δf＝f_s/N (4)

式(4)中，f_s为采样频率；

令则所述窗函数的离散傅立叶变换函数W(r)的公式为：

$W\left(r\right)=\frac{Nrsin\left(\mathrm{\π}r\right)}{\mathrm{\π}}{e}^{-j\mathrm{\π}r}{e}^{j\frac{\mathrm{\π}}{N}r}\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{a_m}{r^2-m^2}---\left(5\right)$

式(5)中，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158；

忽略负频点处谱峰的旁瓣影响，则将式(3)中所述加窗序列的离散傅里叶变换函数X(k)转换为：

优选的，所述次同步频段的频率范围为10～40Hz。

进一步的，所述离散频点的幅值的极大值对应的谱线为在所述次同步频段内谱线幅值的极大值对应的谱线；

设所述离散频点的幅值的极大值对应的谱线为k_a，其相邻的两根谱线分别为k_a-1和k_a+1，则所述谱线k_a-₁、k_a和k_a+1的幅值分别为y₁＝|X(k_a-₁)|、y₂＝|X(k_a)|和y₃＝|X(k_a+1)|，设所述离散频点的幅值的极大值对应的次同步分量的实际频率为kΔf，第一引入参数δ＝k-k_a，第二引入参数α＝(y₃-y₁)/y₂，其中，所述第一引入参数δ∈(-0.5,0.5)，则根据所述加窗序列的离散傅里叶变换函数将第二引入参数α转换为：

$\mathrm{\α}=\frac{\left|W(1-\mathrm{\δ})\right|-\left|W(-1-\mathrm{\δ})\right|}{W(-\mathrm{\δ})}---\left(7\right)$

式(7)中，W(·)为所述窗函数的离散傅立叶变换函数；

记式(7)中α＝g(δ)，采样多项式拟合的方式确定其反函数δ＝g^-1(α)，获取第一引入参数δ，公式为：

δ＝a₁·α+a₃·α³+…+a_2l+1·α^2l+1 (8)

式(8)中，a₁,a₃,…,a_2l+1为2l+1次逼近多项式的系数，l为多项式拟合系数；

对所述电流信号的频率f进行修正，公式为：

f＝kΔf＝(k_a+δ)Δf (9)

对所述电流信号的幅值A进行修正，公式为：

$A=\frac{2(y_1+2y_2+y_3)}{\left|W(-1-\mathrm{\δ})\right|+2\left|W(-\mathrm{\δ})\right|+\left|W(1-\mathrm{\δ})\right|}---\left(10\right)$

由于所述采样数据的长度N＞＞1，则将式(10)转换为：

A＝N^-1(y₁+2y₂+y₃)v(δ) (11)

式(11)中，采用多项式逼近的方式确定v(δ)的近似计算公式为：

v(δ)＝b₀·δ+b₂·δ²+…+b_2l·δ^2l (12)

式(12)中，b₀,b₂,…,b_2l为2l次逼近多项式的系数；

对所述电流信号的相位进行修正，公式为：

本发明的有益效果：

本发明提供的一种次同步振荡参数检测方法，采用四项Rife-Vincent(III)窗对采样数据进行加窗处理，对频谱泄露引起的误差起到了较好的抑制作用，使用峰值频点附近的三根谱线进行插值修正并推导出插值修正公式，更充分的利用了离散频谱包含的信号信息，从而具有更高的计算精度；综合考虑各类窗函数的主瓣宽度和旁瓣幅值衰减速率，选取四项Rife-Vincent(III)窗函数能有效的减小频谱泄露对谐波测量的影响，对抑制频谱泄露的效果更好。

附图说明

图1是本发明一种次同步振荡参数检测方法的流程图；

图2是本发明实施例中采样数据的信号图；

图3是本发明实施例中对采样数据进行加窗处理的信号图。

具体实施方式

下面结合附图对本发明的具体实施方式作详细说明。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本发明提供的一种次同步振荡参数检测方法，如图1所示，包括：

101.对电流信号进行截断，获取长度为N的采样数据；

其中，所述采样数据的信号图如图2所示；

102.建立长度为N的四项Rife-Vincent(III)窗函数，并根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值；

103.在次同步频段内筛选所述离散频点的幅值的极大值，并利用所述离散频点的幅值的极大值对应的谱线及谱线相邻的两根谱线采用插值修正算法对所述电流信号的幅值、频率和相位进行修正。

其中，所述次同步频段的频率范围为10～40Hz。

具体的，所述四项Rife-Vincent(III)窗函数w(n)的时域表达式为：

$w\left(n\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{a}_{m}cos(2\mathrm{\π}mn/N)---\left(1\right)$

式(1)中，n∈[0,N-1]，N为所述采样数据的长度，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158。

所述步骤102中，根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值，包括：

对所述采样数据进行加窗处理的信号图如图3所示，对所述采样数据进行加窗处理，获取加窗序列x_w(n)，公式为：

x_w(n)＝w(n)·x(n) (2)

式(2)中，w(n)为四项Rife-Vincent(III)窗函数，x(n)为所述采样数据；

对所述加窗序列进行离散傅里叶变换，公式为：

式(3)中，X(k)为所述加窗序列的离散傅里叶变换函数，k∈[0,N-1]，N为所述采样数据的长度，A为所述电流信号的幅值，f为所述电流信号的频率，为所述电流信号的相位，W(·)为窗函数的离散傅立叶变换函数，Δf为频率分辨率，其中，所述频率分辨率的公式为：

Δf＝f_s/N (4)

式(4)中，f_s为采样频率；

令则所述窗函数的离散傅立叶变换函数W(r)的公式为：

$W\left(r\right)=\frac{Nrsin\left(\mathrm{\π}r\right)}{\mathrm{\π}}{e}^{-j\mathrm{\π}r}{e}^{j\frac{\mathrm{\π}}{N}r}\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{a_m}{r^2-m^2}---\left(5\right)$

式(5)中，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158；

忽略负频点处谱峰的旁瓣影响，则将式(3)中所述加窗序列的离散傅里叶变换函数X(k)转换为：

进一步的，所述步骤103中，离散频点的幅值的极大值对应的谱线为在所述次同步频段内谱线幅值的极大值对应的谱线；

设所述离散频点的幅值的极大值对应的谱线为k_a，其相邻的两根谱线分别为k_a-1和k_a+1，则所述谱线k_a-₁、k_a和k_a+1的幅值分别为y₁＝|X(k_a-₁)|、y₂＝|X(k_a)|和y₃＝|X(k_a+1)|，设所述离散频点的幅值的极大值对应的次同步分量的实际频率为kΔf，第一引入参数δ＝k-k_a，第二引入参数α＝(y₃-y₁)/y₂，其中，所述第一引入参数δ∈(-0.5,0.5)，则根据所述加窗序列的离散傅里叶变换函数将第二引入参数α转换为：

$\mathrm{\α}=\frac{\left|W(1-\mathrm{\δ})\right|-\left|W(-1-\mathrm{\δ})\right|}{W(-\mathrm{\δ})}---\left(7\right)$

式(7)中，W(·)为所述窗函数的离散傅立叶变换函数；

记式(7)中α＝g(δ)，调用MATLAB的polyfit(α,δ,m)函数求出拟合多项式g^-1(α)的系数，m为拟合多项式的阶数，采样多项式拟合的方式确定其反函数δ＝g^-1(α)，获取第一引入参数δ，公式为：

δ＝a₁·α+a₃·α³+…+a_2l+1·α^2l+1 (8)

式(8)中，a₁,a₃,…,a_2l+1为2l+1次逼近多项式的系数，l为多项式拟合系数，取值应综合考虑计算复杂程度和结果精确度，此处取3，则δ＝g^-1(α)＝1.035494988·α-0.086114685·α³+0.014890668·α⁵-0.002902830·α⁷；

对所述电流信号的频率f进行修正，公式为：

f＝kΔf＝(k_a+δ)Δf (9)

幅值最大谱线及其左右谱线的幅值均较大，其中含有丰富的频谱泄露信息，通过对y₁、y₂和y₃进行加权平均来对幅值进行修正，对所述电流信号的幅值A进行修正，公式为：

$A=\frac{2(y_1+2y_2+y_3)}{\left|W(-1-\mathrm{\δ})\right|+2\left|W(-\mathrm{\δ})\right|+\left|W(1-\mathrm{\δ})\right|}---\left(10\right)$

由于所述采样数据的长度N＞＞1，则将式(10)转换为：

A＝N^-1(y₁+2y₂+y₃)v(δ) (11)

式(11)中，调用MATLAB的polyfit(δ,v(δ),m)函数求出拟合多项式v(δ)的系数，m为拟合多项式的阶数采用多项式逼近的方式确定v(δ)的近似计算公式为：

v(δ)＝b₀·δ+b₂·δ²+…+b_2l·δ^2l (12)

式(12)中，b₀,b₂,…,b_2l为2l次逼近多项式的系数，l取3，则v(δ)＝0.582078950+0.142388577·δ²+0.018604143·δ⁴+0.001795873·δ⁶；

对所述电流信号的相位进行修正，公式为：

本实施例的数据为包含多次谐波和次同步分量的电流数据，基波频率为50hz，采样频率fs＝4000hz，采样点数N＝4096，电流信号中谐波和次同步分量的参数分别如下表所示：

表1 电流信号中谐波的参数

表2 电流信号中次同步分量的参数

应用本发明所述方法分析得到的次同步分量的误差如下表所示：

表3 次同步分量的误差

结果表明本发明提供的一种次同步振荡参数检测方法，采用四项Rife-Vincent(III)窗对采样数据进行加窗处理，对频谱泄露引起的误差起到了较好的抑制作用。使用峰值频点附近的三根谱线进行插值修正并推导出插值修正公式，更充分的利用了离散频谱包含的信号信息，具有更高的计算精度。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。

Claims (5)

1.一种次同步振荡参数检测方法，其特征在于，所述方法包括：

对电流信号进行截断，获取长度为N的采样数据；

建立长度为N的四项Rife-Vincent(III)窗函数，并根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值；

在次同步频段内筛选所述离散频点的幅值的极大值，并利用所述离散频点的幅值的极大值对应的谱线及谱线相邻的两根谱线采用插值修正算法对所述电流信号的幅值、频率和相位进行修正。

2.如权利要求1所述的方法，其特征在于，所述四项Rife-Vincent(III)窗函数w(n)的时域表达式为：

$w\left(n\right)=\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m{a}_{m}cos(2\mathrm{\π}mn/N)---\left(1\right)$

式(1)中，n∈[0,N-1]，N为所述采样数据的长度，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158。

3.如权利要求1所述的方法，其特征在于，所述根据所述四项Rife-Vincent(III)窗函数确定所述采样数据对应的离散频点的幅值，包括：

对所述采样数据进行加窗处理，获取加窗序列x_w(n)，公式为：

x_w(n)＝w(n)·x(n) (2)

式(2)中，w(n)为四项Rife-Vincent(III)窗函数，x(n)为所述采样数据；

对所述加窗序列进行离散傅里叶变换，公式为：

式(3)中，X(k)为所述加窗序列的离散傅里叶变换函数，k∈[0,N-1]，N为所述采样数据的长度，A为所述电流信号的幅值，f为所述电流信号的频率，为所述电流信号的相位，W(·)为窗函数的离散傅立叶变换函数，Δf为频率分辨率，其中，所述频率分辨率的公式为：

Δf＝f_s/N (4)

式(4)中，f_s为采样频率；

令则所述窗函数的离散傅立叶变换函数W(r)的公式为：

$W\left(r\right)=\frac{Nrsin\left(\mathrm{\π}r\right)}{\mathrm{\π}}{e}^{-j\mathrm{\π}r}{e}^{j\frac{\mathrm{\π}}{N}r}\underset{m=0}{\overset{M-1}{\Σ}}{(-1)}^m\frac{a_m}{r^2-m^2}---\left(5\right)$

式(5)中，M＝4为窗函数的项数，a₀＝1，a₁＝1.43596，a₂＝0.49754，a₃＝0.06158；

忽略负频点处谱峰的旁瓣影响，则将式(3)中所述加窗序列的离散傅里叶变换函数X(k)转换为：

4.如权利要求1所述的方法，其特征在于，所述次同步频段的频率范围为10～40Hz。

5.如权利要求3所述的方法，其特征在于，所述离散频点的幅值的极大值对应的谱线为在所述次同步频段内谱线幅值的极大值对应的谱线；

设所述离散频点的幅值的极大值对应的谱线为k_a，其相邻的两根谱线分别为k_a-1和k_a+1，则所述谱线k_a-1、k_a和k_a+1的幅值分别为y₁＝|X(k_a-1)|、y₂＝|X(k_a)|和y₃＝|X(k_a+1)|，设所述离散频点的幅值的极大值对应的次同步分量的实际频率为kΔf，第一引入参数δ＝k-k_a，第二引入参数α＝(y₃-y₁)/y₂，其中，所述第一引入参数δ∈(-0.5,0.5)，则根据所述加窗序列的离散傅里叶变换函数将第二引入参数α转换为：

$\mathrm{\α}=\frac{\left|W(1-\mathrm{\δ})\right|-\left|W(-1-\mathrm{\δ})\right|}{W(-\mathrm{\δ})}---\left(7\right)$

式(7)中，W(·)为所述窗函数的离散傅立叶变换函数；

记式(7)中α＝g(δ)，采样多项式拟合的方式确定其反函数δ＝g^-1(α)，获取第一引入参数δ，公式为：

δ＝a₁·α+a₃·α³+…+a_2l+1·α^2l+1 (8)

式(8)中，a₁,a₃,…,a_2l+1为2l+1次逼近多项式的系数，l为多项式拟合系数；

对所述电流信号的频率f进行修正，公式为：

f＝kΔf＝(k_a+δ)Δf (9)

对所述电流信号的幅值A进行修正，公式为：

$A=\frac{2\left(y_1+2y_2+y_3\right)}{\left|W\left(-1-\mathrm{\δ}\right)\right|+2\left|W\left(-\mathrm{\δ}\right)\right|+\left|W\left(1-\mathrm{\δ}\right)\right|}---\left(10\right)$

由于所述采样数据的长度N＞＞1，则将式(10)转换为：

A＝N^-1(y₁+2y₂+y₃)v(δ) (11)

式(11)中，采用多项式逼近的方式确定v(δ)的近似计算公式为：

v(δ)＝b₀·δ+b₂·δ²+…+b_2l·δ^2l (12)

式(12)中，b₀,b₂,…,b_2l为2l次逼近多项式的系数；

对所述电流信号的相位进行修正，公式为：