CN103353550A - 一种测量电力系统信号频率及谐波参数的方法 - Google Patents

Abstract

本发明涉及一种测量电力系统信号频率及谐波参数的方法，通过对采样数据加7项Rife-Vincent（I）窗的FFT数据中真实频率点附近的2根最大的谱线进行修正，从而得到谐波的频率，幅值和相位，所述的7项Rife-Vincent（I）窗W(n)为：n=0,1,...,N-1其中：b₀=1;b₁=12/7;b₂=15/14;b₃=10/21;b₄=1/7;b₅=2/77;b₆=1/462;本发明实现简单，能够在数字化变电站谐波分析应用中达到很高的精度，计算量比一般的余弦窗小。

Description

一种测量电力系统信号频率及谐波参数的方法

技术领域

本发明涉及测量频率及谐波参数的方法，尤其是涉及一种基于7项Rife-Vincent窗的双谱线插值FFT的测量频率及谐波参数的方法。 

背景技术

在数字化变电站中，互感器采集的电压电流信号，传送到合并单元（进行初步的数据处理和数据打包），之后由合并单元传送到电能质量监测装置（这个装置完成电能质量各项指标的分析）。电力系统中的电压电流频率是不断变化的，而IEC61850-9-2标准规定采样频率为恒定值，一般为10kHz。采用常规的谐波软件算法会发生严重的频谱泄露和栅栏效应。频谱泄漏和栅栏现象使检测出谐波的幅值、相角和频率均存在较大的误差，不能够满足检测精度的要求。利用插值算法可以解决栅栏问题，但是不能消除因频谱泄漏现象而导致的测量误差；利用加窗后再进行DFT分析的算法可以减少频谱泄漏误差，但是需要权衡窗函数复杂度和计算量。 

针对这个问题，国内学者作了相应的研究，并提出了一些改进方法：如基于FFT的高精度谐波检测方法、基于Nuttall窗和5项Rife-Vincent（I）窗双谱线插值FFT的电力谐波分析方法、矩形卷积窗算法等。 

然而由于以上窗函数的频谱旁瓣特性不够理想，对频谱泄露的抑制作用有限，试验中对信号的频率、幅值和相位的分析误差较大，且迭代公式的次数高、项数多，计算时间较长。因此，高精度且快速的频率测量和谐波分析难以实现。 

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种测量电力系统信号频率及谐波参数的方法，是基于7项Rife-Vincent（I）窗的双谱 线插值FFT的测量频率及谐波参数的方法，本方法在数字化变电站谐波分析中固定采样频率的情况下能减少修正公式的计算量，提高计算精度，实现快速、高精度的频率及谐波测量。 

本发明的上述技术问题主要是通过下述技术方案得以解决的： 

一种测量电力系统信号频率及谐波参数的方法，其特征在于，基于以下公式： 

7项Rife-Vincent（I）窗W(n)为： 

$w\left(n\right)=\underset{i=0}{\overset{K}{\mathrm{\Σ}}}{(-1)}^i{b}_i\cos \left(\frac{2\mathrm{\π}}{N}\mathrm{in}\right),n=\mathrm{0,1},...,N-1$ 式一； 

其中：b₀=1;b₁=12/7;b₂=15/14;b₃=10/21;b₄=1/7;b₅=2/77;b₆=1/462; 

包括以下步骤： 

步骤1，互感器采样电力系统信号，传输到合并单元打包处理，通过数字化信息网络将整理后的数据得x(n)传输到电能质量监测仪； 

步骤2，电能质量监测仪采用7项Rife-Vincent（I）窗对复杂电力信号进行处理后采用FFT插值计算，从而得到基波和谐波的频率、幅值和相位； 

步骤3，有数字化变电站将步骤2得到的基波和谐波的频率、幅值和相位显示在控制层。 

在上述的一种测量电力系统信号频率及谐波参数的方法，所述的步骤2中，具体包括以下子步骤： 

步骤2.1：电能质量监测仪对电力系统信号进行加7项Rife-Vincent窗W(n)截断，得到：x_m(n)=x(n)·W(n)n=0...N-1，N为数据窗的长度；对截断的信号进行离散傅里叶变换，如式二所示，得： 

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

式二中，X_m(k)：x_m(n)的离散傅里叶变换值；A_n,θ_n：第n次谐波的幅 值和初始相位；W(π(k·Δf-n·f₀)/f_s)：窗函数W(n)的离散傅里叶变换值；Δf：离散频率间隔，且

7项Rife-Vincent（I）窗具有优良的旁瓣特性，旁瓣电平为-101.1分贝，旁瓣衰减速率为96分贝/倍频程（优于目前广泛使用的窗函数），旁瓣特性比较如表1。 

7项Rife-Vincent（I）窗函数及其幅频特性如图1a和图1b所示。 

步骤2.2：在步骤2.1中，基波峰值频率f₀=k₀·Δf很难正好位于离散谱线频点上，即k₀不是整数；设峰值点左右两侧的谱线分别为第k₁和k₂条，这两条谱线应是峰值点附近幅值最大和次最大的谱线；显然地，k₁<k₀<k₂=k₁+1；在离散频谱中找到这两条谱线，确定k₁和k₂，这两条谱线幅值分别是y₁和y₂；寻找到相应的k₁和k₂，y₁和y₂的值，即能够在非同步的情况下求得真实的频谱；得到式三： 

$\frac{y_1}{y_2}=\left|\frac{W(2\mathrm{\&pi;}(k_1-k_0)/N)}{W(2\mathrm{\&pi;}(k_2-k_0)/N)}\right|$ 式三； 

步骤2.3：对于给定的窗函数，将目前的双谱线插值法进行优化，计算曲线拟合的系数；由式三能够计算出唯一未知量k₀，由于0≤k₀-k₁≤1，定义参数

α=k₀-k₁-0.5，得到式四： 

$\mathrm{\&beta;}=\left|\frac{W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}+0.5)/N-W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5)/N}{W\left(2\mathrm{\&pi;}(-\mathrm{\&alpha;}+0.5\right)/N+W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5)/N}\right|$ 式四； 

且当N较大时，上式一般可以简化为式五： 

α=g^-1(β)                 式五； 

多项式逼近是一种近似计算复杂连续函数值的数值方法，在此可采用多项式逼近方法计算该式，通过控制多项式逼近的次数，可以有效地控制逼近的精度，而且，随着硬件乘法器在微处理器中的广泛应用，多项式逼近的计算公式易于采用程序代码实现； 

步骤2.4：为了克服单峰谱线修正算法易受到频谱泄漏和噪声干扰影响的缺点，次强谱线的信息也能够用于幅值修正，即：直接对k₁和k₂两根谱线幅值进行加权平均，从而计算出实际的峰值点的幅值，其相应幅值计算公式如式六所示： 

$A=\frac{2(y_2+y_1)}{\left|W\left(\frac{2\mathrm{\&pi;}(0.5-\mathrm{\&alpha;})}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5}{N}\right)\right|}$ 式六； 

定义两根谱线采用的权重与其各自的幅值成正比，对于实系数窗函数采用多项式逼近求出函数的近似计算公式，结果中将不含奇次项，这样，双峰谱线修正算法的计算公式就可改写为式七所示： 

A=N^-1(y₁+y₂)(b₀+b₂β²+...+b_2lβ^2l)            式七； 

相位的计算公式如式八所示 

θ=arg[X(k₁·Δf)]+π/2-arg[W(2π·(k₁-k₀)/N]          式八； 

步骤2.5：根据前面的步骤得到基波频率f₁=f₀，接下来计算第i次谐波的幅值和相角：在(i·f₁-5)～(i·f₁+5)范围中再寻求相应的y₁和y₂，回到步骤2.2，一直到所有谐波计算完毕。 

因此，本发明具有如下优点：通过对采样数据加7项Rife-Vincent（I）窗后的FFT数据中真实频率点附近的2根最大的谱线进行修正（修正公式简单），从而得到谐波的频率，幅值和相位；本发明实现简单，能够在实际 应用中达到很高的精度，计算量比一般的余弦窗小。 

附图说明

图1a是7项Rife-Vincent（I）窗时域图。 

图1b是7项Rife-Vincent（I）窗频域图。 

图2是数字化变电站谐波检测方法示意图。 

图3是谐波信号图。 

图4加窗后的信号图。 

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。 

实施例： 

本发明示意图见图2，其中实际信号的波形示于图3，加窗后的波形图示于图4，谐波幅值和相位误差分析形示于表1. 

以下结合附图对本发明的具体实施方式作详细的说明，但不作为本发明的限定。 

本发明可作为数字化变电站按照IEC61850-9-2标准帧数字采样率对电力信号进行采样的系统中进行谐波测量，目前大部分数字化变电站的采样率为10kHz的采样率，选取一个电压谐波信号x(t)，如图3所示。 

本发明所述的非同步采样的谐波检测方法包含六个步骤： 

步骤a：以固定采样频率f_s=10000Hz采样谐波信号x(t)得：x(n)，对其加7项Rife-Vincent（I）窗得到截断后信号x_m(n)，x_m(n)=x(n)·W(n)n=0...N-1，窗函数的长度N=2046，如图4所示； 

7项Rife-Vincent（I）窗W(n)为： 

$w\left(n\right)=\underset{i=0}{\overset{K}{\mathrm{\&Sigma;}}}{(-1)}^i{b}_i\cos \left(\frac{2\mathrm{\&pi;}}{N}\mathrm{in}\right),n=\mathrm{0,1},...,N-1---\left(1\right)$

其中：b₀=1;b₁=12/7;b₂=15/14;b₃=10/21;b₄=1/7;b₅=2/77;b6₌1/462; 

对截断的信号进行傅里叶变换计算， 

如式（2）所示，得： 

$X_m\left(k\right)=\frac{N}{2}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}(A_n/2)\&times;\left[e^{j{\mathrm{\&theta;}}_n}W(\mathrm{\&pi;}(k\&CenterDot;\mathrm{\&Delta;f}-n\&CenterDot;f_0)/f_s)\right]---\left(2\right);$

得到一个待修正的幅频特性和相频特性； 

步骤b：在45～55Hz内寻找最大的相邻的两根谱线，分别为第k₁和k₂根谱线，相应幅值分别是y₁和y₂。并且得到式（9）。 

$\frac{y_1}{y_2}=\left|\frac{W(2\mathrm{\&pi;}(k_1-k_0)/N)}{W(2\mathrm{\&pi;}(k_2-k_0)/N)}\right|---\left(9\right)$

步骤c：求系数α和β，如式（10）、（11）所示。 

β=(y2-y1)/(y2+y1)        （10） 

α=6.5*β       （11） 

步骤d：各次谐波的幅值和相角计算公式，如（12）、（13）所示。 

A=N^-1·(y₁+y₂)·(0.10392+0.62613*β²+1.99889*β⁴)       （12） 

θ=arg[X(k_i·Δf)]-π·(α-(-1)ⁱ·0.5)...(i=1,2)       （13） 

步骤e：根据前面的步骤可以得到基波频率f₁=f₀，接下来计算第i次谐波的幅值和相角：在(i·f₁-5)～(i·f₁+5)范围中再寻求相应的y₁和y₂，回到步骤b，一直到所有谐波计算完毕。 

步骤f：进行相关的误差分析。 

本实施算例的信号频率为50.5Hz，基波频率计算相对误差为1.46×10^-10。基波、谐波参数和误差分析见表1。 

表1误差分析表 

本算例实施还采用了Nuttall窗、Blackman-harris窗、5项Rife-Vincent（I）窗进行了实验。结果表明采用基于7项Rife-Vincent（I）窗的双谱线插值FFT的频率及谐波测量方法的修正公式项数少，计算量少，计算精度更高。 

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。 

Claims (2)

1.一种测量电力系统信号频率及谐波参数的方法，其特征在于，基于以下公式：

7项Rife-Vincent I窗W(n)为：

$w\left(n\right)=\underset{i=0}{\overset{K}{\mathrm{\&Sigma;}}}{(-1)}^i{b}_i\cos \left(\frac{2\mathrm{\&pi;}}{N}\mathrm{in}\right),n=\mathrm{0,1},...,N-1$ 式一；

其中：b₀=1;b₁=12/7;b₂=15/14;b₃=10/21;b₄=1/7;b₅=2/77;b₆=1/462;

包括以下步骤：

步骤1，互感器采样电力系统信号，传输到合并单元打包处理，通过数字化信息网络将整理后的数据得x(n)传输到电能质量监测仪；

步骤2，电能质量监测仪采用7项Rife-Vincent I窗对复杂电力信号进行处理后采用FFT插值计算，从而得到基波和谐波的频率、幅值和相位；

步骤3，有数字化变电站将步骤2得到的基波和谐波的频率、幅值和相位显示在控制层。

2.根据权利要求1所述的一种测量电力系统信号频率及谐波参数的方法，其特征在于，所述的步骤2中，具体包括以下子步骤：

步骤2.1：电能质量监测仪对电力系统信号进行加7项Rife-Vincent窗W(n)截断，得到：x_m(n)=x(n)·W(n)n=0...N-1，N为数据窗的长度；对截断的信号进行离散傅里叶变换，如式二所示，得：

$X_m\left(k\right)=\frac{N}{2}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}(A_n/2)\&times;\left[e^{j{\mathrm{\&theta;}}_n}W(\mathrm{\&pi;}(k\&CenterDot;\mathrm{\&Delta;f}-n\&CenterDot;f_0)/f_s)\right]$ 式二；

式二中，X_m(k)：x_m(n)的离散傅里叶变换值；A_n,θ_n：第n次谐波的幅值和初始相位；W(π(k·Δf-n·f₀)/f_s)：窗函数W(n)的离散傅里叶变换值；Δf：离散频率间隔，且

步骤2.2：在步骤2.1中，基波峰值频率f₀=k₀·Δf很难正好位于离散谱线频点上，即k₀不是整数；设峰值点左右两侧的谱线分别为第k₁和k₂条，这两条谱线应是峰值点附近幅值最大和次最大的谱线；显然地，k₁<k₀<k₂=k₁+1；在离散频谱中找到这两条谱线，确定k₁和k₂，这两条谱线幅值分别是y₁和y₂；寻找到相应的k₁和k₂，y₁和y₂的值，即能够在非同步的情况下求得真实的频谱；得到式三：

$\frac{y_1}{y_2}=\left|\frac{W(2\mathrm{\&pi;}(k_1-k_0)/N)}{W(2\mathrm{\&pi;}(k_2-k_0)/N)}\right|$ 式三；

步骤2.3：对于给定的窗函数，将目前的双谱线插值法进行优化，计算曲线拟合的系数；由式三能够计算出唯一未知量k₀，由于0≤k₀-k₁≤1，定义参数α=k₀-k₁-0.5，得到式四：

$\mathrm{\&beta;}=\left|\frac{W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}+0.5)/N-W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5)/N}{W\left(2\mathrm{\&pi;}(-\mathrm{\&alpha;}+0.5\right)/N+W(2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5)/N}\right|$ 式四；

且当N较大时，上式一般可以简化为式五：

α=g^-1(β)             式五；

步骤2.4：为了克服单峰谱线修正算法易受到频谱泄漏和噪声干扰影响的缺点，次强谱线的信息也能够用于幅值修正，即：直接对k₁和k₂两根谱线幅值进行加权平均，从而计算出实际的峰值点的幅值，其相应幅值计算公式如式六所示：

$A=\frac{2(y_2+y_1)}{\left|W\left(\frac{2\mathrm{\&pi;}(0.5-\mathrm{\&alpha;})}{N}\right)\right|+\left|W\left(\frac{2\mathrm{\&pi;}(-\mathrm{\&alpha;}-0.5}{N}\right)\right|}$ 式六；

定义两根谱线采用的权重与其各自的幅值成正比，对于实系数窗函数采用多项式逼近求出函数的近似计算公式，结果中将不含奇次项，这样，双峰谱线修正算法的计算公式就可改写为式七所示：

A=N^-1(y₁+y₂)(b₀+b₂β²+...+b_2lβ^2l)               式七；

相位的计算公式如式八所示

θ=arg[X(k₁·Δf)]+π/2-arg[W(2π·(k₁-k₀)/N]         式八；

步骤2.5：根据前面的步骤得到基波频率f₁=f₀，接下来计算第i次谐波的幅值和相角：在(i·f₁-5)～(i·f₁+5)范围中再寻求相应的y₁和y₂，回到步骤2.2，一直到所有谐波计算完毕。

Images