CN105486921A - 凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法 - Google Patents

Abstract

本发明公开了一种凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法，包括：步骤1，互感器采样电力系统的信号获得离散化信号，对离散化信号加Kaiser三阶互卷积窗截断，得截断信号，对截断信号进行离散傅里叶变换；步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α；步骤3，计算各位置α对应的γ值；步骤4，对各位置α及其对应的γ值进行多项式拟合得α和γ的函数关系；步骤5，利用三峰谱线修正法计算n次谐波下真实谱线k₀的幅值和相位。本发明可有效抑制频谱泄露现象，从而可极大提高谐波和间谐波的检测精度。

Description

凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法

技术领域

本发明涉及频率、谐波与间谐波参数测量技术领域，具体涉及一种凯撒(Kaiser)三阶互卷积窗三谱线插值的谐波与间谐波检测方法。

背景技术

电力系统中非线性负荷的大量增加，特别是电力电子装置的广泛应用，使电网中谐波和间谐波成分增加，引起配网中电机、变压器等设备损耗增加，噪声增加，同时也会加剧电力电缆的老化。对谐波和间谐波分量的准确分析将有利于电能质量的分析和治理。快速傅里叶变换由于其成熟的理论和易于电网嵌入式实现的原因已成为目前电网谐波分析最主要的方法，但直接采用快速傅里叶变换对电网谐波进行分析时由于非整周期截取和非同步采样会产生频谱泄露和栅栏效应，严重影响信号参数测量结果的准确性，以至于无法达到谐波测量的国标要求。经研究发现加窗后再进行快速傅里叶分析可有效减少由于非整周期截取造成的频谱泄露，而与不同窗函数相适应的插值算法可以有效解决栅栏效应带来的测量误差。不同的窗函数其时域表达式和频谱特性不同，相应的计算得出的插值和参数修正公式也不同，考虑到电网谐波和间歇波参数测量的时效性要求，应在窗函数的精确度和计算量之间进行权衡。

目前现有的一些谐波分析方法正是针对上述问题所提出的。例如，基于插值FFT算法的谐波参数估计、基于FFT算法的分次谐波测量与分析、基于FFT的高精度窗函数双谱线插值FFT的电力谐波分析方法、矩形卷积窗法、三角卷积窗法等。

然而上述窗函数的频谱特性不够理想，在实际应用中对频谱泄露现象的抑制效果有限，因此其在谐波和间谐波检测过程中对信号参数的测量误差较大难以达到高精度的分析测量；在仿真分析中可以看出需要多达二十个工频周波的采样数据才能实现较为准确的测量，在实时性方面不满足要求。

发明内容

针对现有技术中存在的不足，本发明提出了一种凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法。本发明在电力系统谐波和间谐波分析中，采用恒定采样频率，仅通过十个工频周期的采样数据就能实现快速和高精度的谐波与间谐波参数测量。

本发明通过对十个工频周期的采样数据加Kaiser三阶互卷积窗后的FFT数据中真实频率点附近的三条最大谱线进行修正，通过最小均方算法简化修正公式，从而进一步分析谐波和间谐波参数。

为解决上述技术问题，本发明采用如下的技术方案：

步骤1，互感器采样电力系统的信号获得离散化信号，对离散化信号加Kaiser三阶互卷积窗截断，得截断信号，对截断信号进行离散傅里叶变换；

步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α，其中，n为当前谐波次数，f₀为基波频率；

步骤3，计算各位置α对应的γ值，具体为：在离散傅里叶变换后的截断信号中找出距位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，对应的幅值分别为y₁、y₂和y₃，计算

步骤4，对各位置α及其对应的γ值进行多项式拟合得α和γ的函数关系；

步骤5，将(n·f₀-6)～(n·f₀+6)范围内幅值最大的谱线的位置α及其对应的y₁、y₂和y₃代入A＝N^-1(y₁+2y₂+y₃)·(b₀+b₂α²+…+b_2lα^2l)，所得A值即真实谱线的幅值；l为根据实际情况人为设定的正整数，系数b₀、b₂……b_2l采用如下方法获得：

计算各位置α对应的基于各位置α及其对应的A、y₁、y₂和y₃，采用多项式拟合函数拟合A＝N^-1(y₁+2y₂+y₃)·(b₀+b₂α²+…+b_2lα^2l)中的系数b₀、b₂……b_2l；

步骤6，采用离散傅里叶变换后的连续信号和离散谱线计算真实谱线的相位；

步骤7，重复步骤2～6直至所有谐波次数下真实谱线的幅值和相角计算完毕。

上述Kaiser三阶互卷积窗w(n)为：

$w\left(n\right)=\frac{I_0\[{\mathrm{\β}}_1{(1-{\[\left(n_1-{\mathrm{\α}}_1\right)/{\mathrm{\α}}_1\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_1\right)}*\frac{I_0\[{\mathrm{\β}}_2{(1-{\[\left(n_2-{\mathrm{\α}}_2\right)/{\mathrm{\α}}_2\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_2\right)}*\frac{I_0\[{\mathrm{\β}}_3{(1-{\[\left(n_3-{\mathrm{\α}}_3\right)/{\mathrm{\α}}_3\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_3\right)}$

其中，I₀(·)表示第一类零阶修正Bessel函数；β_i表示Kaiser窗的形状参数，为经验值，可在5～14范围内取值；n_i为经验值；α_i＝n_i/2，i＝1,2,3。

与现有技术相比，本发明具有以下优点和有益效果：

1、可有效抑制频谱泄露现象，从而可极大提高谐波和间谐波的检测精度。

2、易于工程实现，可在实际应用中获得高精度的谐波和简谐波检测结果。

3、便于嵌入式系统实现，可连续、长期对被测信号进行检测。

附图说明

图1是实施例1中谐波信号示意图；

图2是实施例1中窗函数的时域和频谱特性图，其中，图(a)为窗函数时域特性图，图(b)为窗函数的频谱特性图；

图3是实施例1中截断信号示意图；

图4是实施例1中间谐波信号图；

图5是实施例2中截断信号示意图。

具体实施方式

本发明的具体步骤如下：

步骤1，采样电力系统的信号x(t)，并获得信号x(t)的离散化的信号x(n)，对信号x(n)加Kaiser三阶互卷积窗，得到截断信号x_m(n)，对截断信号x_m(n)进行离散傅里叶变换。

互感器以固定频率f_s对电力系统的信号进行采样，信号可以为电压信号或电流信号；将采样信号传输到合并单元按时段进行打包，通过数字化信息网络将信号x(n)传输到电能质量检测仪。

对信号x(n)加Kaiser三阶互卷积窗进行截断，得到截断信号x_m(n)：

x_m(n)＝x(n)·w(n)(1)

Kaiser三阶互卷积窗函数w(n)如下：

$w\left(n\right)=\frac{I_0\[{\mathrm{\β}}_1{(1-{\[\left(n_1-{\mathrm{\α}}_1\right)/{\mathrm{\α}}_1\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_1\right)}*\frac{I_0\[{\mathrm{\β}}_2{(1-{\[\left(n_2-{\mathrm{\α}}_2\right)/{\mathrm{\α}}_2\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_2\right)}*\frac{I_0\[{\mathrm{\β}}_3{(1-{\[\left(n_3-{\mathrm{\α}}_3\right)/{\mathrm{\α}}_3\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_3\right)}$

(2)

式(2)中：

n表示采样点编号，n＝0,1,…,N-1，N表示窗长度；

I₀(·)表示第一类零阶修正Bessel函数，

β_i表示Kaiser窗的形状参数，为经验值，可在5～14范围内取值；

n_i为经验值，一般可根据IEC61044标准取n₁、n₂、n₃分别为512、512、1024；当n₁＝n₂＝n₃、β₁＝β₂＝β₃时，w(n)为三阶自卷积窗函数，否则，w(n)为三阶互卷积窗函数；

α_i＝n_i/2，i＝1,2,3。

对截断信号x_m(n)进行离散傅里叶变换，经离散傅里叶变换即可得截断信号x_m(n)中所有谱线函数。离散傅里叶变换如下：

$X_m\left(k\right)=\frac{N}{2}\underset{n=0}{\overset{N-1}{\Σ}}(A_n/2)\×\[e^{{\mathrm{j\θ}}_n}W\left(\mathrm{\π}\right(k\·\mathrm{\Δ}f-n\·f_0)/f_s)\]---\left(3\right)$

式(3)中：

X_m(k)表示截断信号x_m(n)的离散傅里叶变换值；

A_n和θ_n分别表示第n次谐波的幅值和初始相位；

W(π(k·Δf-n·f₀)/f_s)表示Kaiser三阶互卷积窗函数w(π(k·Δf-n·f₀)/f_s)的离散傅里叶变换值，k为自变量，表示谱线序号；f₀为基波频率，f_s为采样频率，Δf为离散频率间隔，且

步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α，在离散傅里叶变换后的截断信号中找出距离位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，并假设k₁和k₂分布于位置α两侧，峰值谱线k₁、k₂、k₃的幅值分别记为y₁、y₂和y₃。

k₁、k₂、k₃同时表示峰值谱线的位置。上述n表示谐波次数，基波频率f₀＝(α+k_i+0.5)Δf，k_i可以为k₁、k₂或k₃。

步骤3，计算各位置α对应的γ值，经多项式拟合得α和γ值的函数关系。

计算各位置α对应的γ值具体为：

在离散傅里叶变换后的截断信号中找出距位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，对应的幅值分别为y₁、y₂和y₃，计算

$\mathrm{\γ}=\frac{y_3-y_1}{y_2}=\frac{\left|W\left(\frac{2\mathrm{\π}(1-\mathrm{\α})}{N}\right)\right|-\left|W\left(\frac{2\mathrm{\π}(-1-\mathrm{\α})}{N}\right)\right|}{\left|W\left(\frac{2\mathrm{\π}(-\mathrm{\α})}{N}\right)\right|}---\left(4\right)$

α＝g^-1(γ)(5)

式(4)～(5)中：

g^-1(γ)表示公式(4)所示函数的反函数；

α表示真实谱线k₀的位置；

$W\left(\frac{2\mathrm{\π}(1-\mathrm{\α})}{N}\right),W\left(\frac{2\mathrm{\π}(-1-\mathrm{\α})}{N}\right),W\left(\frac{2\mathrm{\π}(-\mathrm{\α})}{N}\right)$ 分别表示峰值谱线k₃、k₁、k₂的Kaiser三阶互卷积窗函数 $w\left(\frac{2\mathrm{\π}(1-\mathrm{\α})}{N}\right),w\left(\frac{2\mathrm{\π}(-1-\mathrm{\α})}{N}\right),w\left(\frac{2\mathrm{\π}(-\mathrm{\α})}{N}\right)$ 的离散傅里叶变换值；

N表示窗长度。

步骤4，利用三峰谱线修正法计算n次谐波下真实谱线k₀的幅值和相位。

幅值A计算公式如下：

A＝N^-1(y₁+2y₂+y₃)·(b₀+b₂α²+…+b_2lα^2l)(6)

式(6)中：

l为根据实际情况人为设定的正整数，一般l取值越大，计算精度越高，同时计算效率降低，具体实施时，需同时考虑计算精度和计算效率两个因素人为设定l值。

b₀、b₂……b_2l为α的系数，通过多项式拟合得到，具体为：

计算各位置α对应的基于各位置α及其对应的A、y₁、y₂和y₃，采用多项式拟合函数拟合公式(6)中的b₀、b₂……b_2l。

将(n·f₀-6)～(n·f₀+6)范围内幅值最大的谱线的位置α及其对应的y₁、y₂和y₃代入公式(6)，所得A值即真实谱线的幅值。(n·f₀-6)～(n·f₀+6)范围内幅值最大的谱线的位置α可作为真实谱线位置。

相位θ计算公式如下：

θ＝arg[X_m(k₁·Δf)]+π/2-arg[W(2π·(k₁-k₀)/N)](7)

式(7)中：

k₀表示截断信号经傅里叶变换后得到的连续谱线的包络最大值位置，即真实谱线位置；

k₁是截断信号经离散傅里叶变换后的离散谱线中最大值谱线位置，即峰值谱线k₁的位置。

重复步骤2～4直至所有谐波和间歇波的幅值和相角计算完毕。

实施例1

谐波检测

选取含高次谐波的电压信号x(t)，见图1，下面将以图1所示电压信号为例说明本发明的谐波检测方法。

步骤1，固定采样频率f_s＝10000Hz采样谐波信号x(t)，经离散化得信号x(n)，对信号x(n)加Kaiser三阶互卷积窗w(n)，得截断信号x_m(n)，并对截断信号x_m(n)进行离散傅里叶变换。

本实施例中，Kaiser三阶互卷积窗w(n)的形状参数β₁＝β₂＝β₃＝8，窗长度参数n₁＝n₂＝512、n₃＝1024，进行三阶互卷积后补零得到长度N＝2048的窗函数。窗函数的时频特性见图2，截断信号x_m(n)见图3。

步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α，在离散傅里叶变换后的截断信号中找出距离位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，并假设k₁和k₂分布于位置α两侧，峰值谱线k₁、k₂、k₃的幅值分别记为y₁、y₂和y₃。

步骤3，计算各位置α对应的γ值，经多项式拟合得α和γ值的函数关系。

本实施例所得α和γ值的函数关系如下：

α＝-0.0031·γ⁷+0.0170·γ⁵-0.1120·γ³+1.4983·γ(8)

步骤4，利用三峰谱线修正法计算真实谱线的幅值和相位。

采用公式(9)计算真实谱线的幅值：

A＝N^-1·(y₁+2·y₂+y₃)·(0.0077·α⁴+0.0856·α²+0.5102)·10^-4(9)

采用公式(10)相位计算公式为：

θ＝arg[X_m(k₁·Δf)]-π·(α-(-1)ⁿ·0.5)(10)

其中，n表示谐波次数。

步骤5，根据前面的步骤可以得到基波频率f₁＝f₀，重复步骤2～4直至所有谐波和间歇波的幅值和相角计算完毕。

本实施例中信号频率为50.5赫兹，基波频率计算相对误差为4.48×10^-10，幅值和相位及误差分析见表1。

表1谐波信号的误差分析表

实施例2：间谐波检测

选取含间谐波的电压信号x(t)，见图4，下面将以图4所示电压信号为例说明本发明的间谐波检测方法。

步骤1，以固定采样频率f_s＝10000Hz采样谐波信号x(t)得x(n)，对信号x(n)加Kaiser三阶互卷积窗，得到截断后号x_m(n)，β₁＝β₂＝β₃＝8。

本实施例中，窗函数长度n₁＝n₂＝512、n₃＝1024，进行三阶卷积后补零得到长度N＝2048的窗函数，加窗后信号，即截断信号x_m(n)见图5。

步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α，在离散傅里叶变换后的截断信号中找出距离位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，并假设k₁和k₂分布于位置α两侧，峰值谱线k₁、k₂、k₃的幅值分别记为y₁、y₂和y₃。

步骤3，计算各位置α对应的γ值，经多项式拟合得α和γ值的函数关系。

本实施例所得α和γ值的函数关系如下：

α＝-0.0031·γ⁷+0.0170·γ⁵-0.1120·γ³+1.4983·γ(11)

步骤4，利用三峰谱线修正法计算真实谱线的幅值和相位。

采用公式(12)计算真实谱线的幅值：

A＝N^-1·(y₁+2·y₂+y₃)·(0.0077·α⁴+0.0856·α²+0.5102)·10^-4(12)

采用公式(13)相位计算公式为：

θ＝arg[X_m(k₁·Δf)]-π·(α-(-1)ⁿ·0.5)(13)

其中，n表示谐波次数。

步骤5，根据前面的步骤可以得到基波频率f₁＝f₀，重复步骤2～4直至所有谐波和间歇波的幅值和相角计算完毕。

本实施的信号频率以及幅值和相位误差分析如表2所示。

表2间谐波信号的误差分析表

结果表明采用Kaiser三阶互卷积窗三谱线插值FFT的谐波与间谐波检测方法有着较高的计算精度，满足国家标准。

Claims (2)

1.凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法，其特征是，包括：

步骤1，互感器采样电力系统的信号获得离散化信号，对离散化信号加Kaiser三阶互卷积窗截断，得截断信号，对截断信号进行离散傅里叶变换；

步骤2，在(n·f₀-6)～(n·f₀+6)范围按预设步长取一系列值，记为位置α，其中，n为当前谐波次数，f₀为基波频率；

步骤3，计算各位置α对应的γ值，具体为：在离散傅里叶变换后的截断信号中找出距位置α最近的三条峰值谱线，按位置顺序顺次记为峰值谱线k₁、k₂、k₃，对应的幅值分别为y₁、y₂和y₃，计算 $\mathrm{\γ}=\frac{y_3-y_1}{y_2};$

步骤4，对各位置α及其对应的γ值进行多项式拟合得α和γ的函数关系；

步骤5，将(n·f₀-6)～(n·f₀+6)范围内幅值最大的谱线的位置α及其对应的y₁、y₂和y₃代入A＝N^-1(y₁+2y₂+y₃)·(b₀+b₂α²+...+b_2lα^2l)，所得A值即真实谱线的幅值；l为根据实际情况人为设定的正整数，系数b₀、b₂……b_2l采用如下方法获得：

计算各位置α对应的基于各位置α及其对应的A、y₁、y₂和y₃，采用多项式拟合函数拟合A＝N^-1(y₁+2y₂+y₃)·(b₀+b₂α²+...+b_2lα^2l)中的系数b₀、b₂……b_2l；

步骤6，采用离散傅里叶变换后的连续信号和离散谱线计算真实谱线的相位；

步骤7，重复步骤2～6直至所有谐波次数下真实谱线的幅值和相角计算完毕。

2.如权利要求1所述的凯撒三阶互卷积窗三谱线插值的谐波与间谐波检测方法，其特征是：

所述的Kaiser三阶互卷积窗w(n)为：

$w\left(n\right)=\frac{I_0\[{\mathrm{\β}}_1{(1-{\[(n_1-{\mathrm{\α}}_1)/{\mathrm{\α}}_1\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_1\right)}*\frac{I_0\[{\mathrm{\β}}_2{(1-{\[(n_2-{\mathrm{\α}}_2)/{\mathrm{\α}}_2\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_2\right)}*\frac{I_0\[{\mathrm{\β}}_3{(1-{\[(n_3-{\mathrm{\α}}_3)/{\mathrm{\α}}_3\]}^2)}^{\frac{1}{2}}\]}{I_0\left({\mathrm{\β}}_3\right)}$

其中，I₀(·)表示第一类零阶修正Bessel函数；β_i表示Kaiser窗的形状参数，为经验值，可在5～14范围内取值；n_i为经验值；α_i＝n_i/2，i＝1,2,3。