一种基于离散傅里叶变换的频域插值电力谐波分析方法

摘要

本发明公开了一种基于离散傅里叶变换的频域插值电力谐波分析方法，包括以下步骤，(1)？时域连续信号离散化，(2)？根据测量要求，选取合适的组合余弦窗对离散采样序列加权，(3)？对序列进行快速傅里叶变换（4）以谱峰搜索的形式寻找所测谐波附近最大的两条谱线，并计算他们的谱值比，（5）根据牛顿插值方法获得频率偏差，（6）根据求得的频率偏差，进一步求解该次谐波的幅值、频率和初相位。本发明的基于查表的快速插值与窗函数无关，具有普适性，可与任何窗函数结合进行谐波分析，同时，该方法无复杂运算，计算量较小，易于嵌入式系统中实现。

说明书

技术领域

本发明涉及一种基于离散傅里叶变换的频域插值电力谐波分析方法，属于电力系统谐波 测量领域。

背景技术

随着电网智能化的发展，由于越来越多的分布式电源并网运行，使得谐波污染越来越严 重，也使谐波问题更加复杂。准确的谐波测量是实施管理和控制电网谐波的前提，因此，谐 波测量技术受到了广泛关注。

目前，研究者已提出较多谐波参数估计的方法。在这些方法中，基于DFT的频谱分析方 法，是当前应用最广泛的电力系统谐波参数估计方法。DFT方法物理意义明确，并有FFT提 高信号估计的实时性，故被IEC(InternationalElectrotechnicalCommission,IEC)标准所推荐。

DFT方法假定信号是周期性的。但由于实际电网负荷变动导致基波频率波动和间谐波的 影响，使得整周期采样不可实现。即使采用数字锁相环技术等硬件同步，也会有一个或几个 采样周期的延时，不能完全实现同步采样。因此，常采用加窗插值方法对FFT分析结果进行 频谱校正。旁瓣性能较好的窗函数可以抑制信号由于非整周期截断造成的频谱泄漏，插值方 法则能够消除栅栏效应误差。

加窗插值FFT方法在一定程度上提高了谐波分析的精度。目前，已有较多的窗函数和插 值方法被提出。但是，在这些加窗插值方法中，只有基于Rife-Vincent(I)窗的插值方法，在 频率偏差与谱值比之间有解析解，其他形式的加窗插值方法只能通过求解高阶方程、迭代等 方式，求得频率偏差的近似解，求解高阶方程或迭代运算，大大增加了谐波分析过程中的 计算量，影响了相关方法在实际中的应用。

发明内容

本发明提出一种基于离散傅里叶变换的频域插值电力谐波分析方法，该方法独立于窗 函数项数和类型，根据查表和牛顿插值方法来获得非整周期采样条件下的频率偏差，进而求 得谐波幅值、频率和初始相位，该插值方法具有较高的谐波分析准确度，同时降低谐波分析 过程的计算量，易于在嵌入式系统中实现。

为达到上述目的，本发明采用的技术方案如下：

一种基于离散傅里叶变换的频域插值电力谐波分析方法，包括以下步骤：

(1)时域连续信号离散化：将电力系统中含有谐波成分的电压/电流时域连续信号x(t) 通过低通滤波器后，以固定的采样频率f_s进行采样，得到离散采样序列x(n)，其中，n表示 第n个采样点；

(2)根据测量要求，选取合适的组合余弦窗对离散采样序列x(n)加权，得到N点长 序列h_N(n)，其中，N表示采样点数；

(3)对序列h_N(n)进行快速傅里叶变换，然后对其计算结果取模求绝对值，获得包括 基波在内的各次谐波幅值；

(4)以谱峰搜索的形式寻找所测谐波附近最大的两条谱线，并计算谱值比γ：

γ＝次最大谱线值/最大谱线值；

(5)根据牛顿插值方法获得频率偏差；

(6)根据所述步骤(5)求得的频率偏差，进一步求解该次谐波的幅值、频率和初相位。 前述的步骤(3)中，各次谐波幅值的求解方法如下：

3-1)对h_N(n)进行快速傅里叶变换，根据傅里叶变换频域卷积定理，第r次谐波的离散 频谱H(λ)为：

$H\left(\lambda \right)=\frac{A_r}{2}\underset{r=1}{\overset{R}{\Sigma }}[G(\lambda -{\lambda }_r)e^{-j[\pi (\lambda -{\lambda }_r)-{\phi }_r]}+G(\lambda +{\lambda }_r)e^{-j[\pi (\lambda +{\lambda }_r)+{\phi }_r]}]---\left(2\right)$

式中：G(λ)为余弦窗函数g(n)的频域函数；A_r和φ_r分别为第r次谐波的幅值和初始相 位，R为谐波最高次数，λ_r为第r次谐波频率被频率分辨率归一化后的值，λ为被频率分辨 率归一化的频率；

3-2)在不计谐波间相互干扰及负频次谐波影响情况下，式(2)可整理为：

$H\left(\lambda \right)=\frac{A_r}{2}G\left(\sigma \right)e^{-j(\pi \sigma -{\phi }_r)}---\left(3\right)$

式中，σ＝λ-λ_r为频率偏差，根据式(3)，第r次谐波幅值A_r表示为：

$A_r=\frac{2\left|H\left(\lambda \right)\right|}{\left|G\left(\sigma \right)\right|}---\left(4\right).$

前述的步骤(5)根据牛顿插值方法获得频率偏差，包括以下步骤：

5-1)在组合余弦窗从中心位置到半个归一化频率[0,0.5]区间内，将频率均分为P等份， 则均分点处对应的频率值λ_i等于i＝0,1...P，计算均分点处频率λ_i和λ_i-1处的主 瓣幅值，计算方法如下：

组合余弦窗的离散时域表达式g(n)为：

$g\left(n\right)=\underset{m=0}{\overset{M-1}{\Sigma }}{(-1)}^m{a}_{m}cos\left(\frac{2\pi }{N}mn\right),n=0,1,2...N-1---\left(5\right)$

其中，M为窗函数的项数，a_m为窗函数系数，

式(5)经离散傅里叶变换，得到余弦窗频域函数G(λ)：

$G\left(\lambda \right)=\underset{m=0}{\overset{M-1}{\Sigma }}{(-1)}^m\frac{a_m}{2}[D(\lambda -m)+D(\lambda +m)]---\left(6\right)$

式中，λ为被频率分辨率归一化的频率，

把λ_i和λ_i-1的值作为自变量λ代入式(6)求绝对值，即可获得对应的主瓣幅值；

5-2)计算λ_i-1处主瓣幅值与λ_i处主瓣幅值的比值，该比值为一组离散数据，用ζ_i表 示：

${\zeta }_i=\frac{\left|G({\lambda }_i-1)\right|}{\left|G\left({\lambda }_i\right)\right|}=\frac{\left|G(i·\frac{0.5}{P}-1)\right|}{\left|G(i·\frac{0.5}{P})\right|},i=0,1.P.---\left(7\right);$

5-3)以ζ_i为自变量，λ_i为因变量，构建牛顿插值多项式f(x)如下：

f(x)＝f[x₀]+f[x₀,x₁](x-x₀)

+f[x₀,x₁,x₂](x-x₀)(x-x₁)+···(8)

+f[x₀,x₁,···,x_i](x-x₀)···(x-x_i-1)

式中，

$f[x_0,x_1,...,x_i]=\frac{f[x_0,x_1···,x_{i-2},x_i]-f[x_0,x_1...,x_{i-1}]}{x_i-x_{i-1}}$

其中，x₀,x₁,···,x_i分别对应ζ₀,ζ₁,···,ζ_i的计算值，f(x₀),f(x₁),···,f(x_i)分别对应 λ₀,λ₁,···,λ_i的数值；

5-4)利用构造的牛顿插值多项式，计算谱值比为γ时对应的频率偏差σ。

前述的步骤5-1)中，P≤10。

前述的步骤(6)中，谐波的幅值、频率和初相位的求解方法为：

将所述步骤(5)求得的σ作为自变量λ代入式(6)求绝对值，即可获得G(σ)的值，H(λ) 为最接近第r次谐波离散频谱H(λ)的幅值，也为已知，故根据式求得第r次 谐波幅值A_r；

设频域轴上的第p根谱线最接近第r次谐波谱线，则第r次谐波频率f_r可表示为：

$f_r={\mathrm{rf}}_0=\frac{(p-1-\sigma )}{T_{s}N}---\left(9\right)$

其中，f₀为基波频率，T_s为数据采集系统的采样周期；

最后，第r次谐波的初相位φ_r为：

φ_r＝φ_p+πσ(10)

式中，φ_p为频域轴上最接近第r次谐波的第p根谱线的初相角。

本发明的优点是：1、具有普遍适用性，在本发明的公式推导过程中并没有指定谐波次 数，故该方法可用于任意次数谐波的参数求解；2、方法较为简单，无复杂运算，计算量小， 易于在嵌入式系统中实现。

附图说明

图1本发明的实施例中余弦窗主瓣[0,0.5]区间内，i取不同值时的λ_i和λ_i-1位置关 系；

图2本发明的实施例中四项Blackman-Harris窗归一化频率[0,0.5]区间内，均分为10份条 件下，由本发明方法计算得到频率偏差与理论值之间的误差。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的说明。

本发明的基于离散傅里叶变换的频域插值电力谐波分析方法包括以下步骤：

(1)时域连续信号离散化：将电力系统中含有谐波成分的电压/电流时域连续信号x(t) 通过低通滤波器后，以固定的采样频率f_s进行采样，得到离散采样序列x(n)，其中，n表示 第n个采样点。

非整周期采样下，一段电网采样信号，其时域离散采样序列x(n)可表示为：

$x\left(n\right)=\underset{r=1}{\overset{R}{\Sigma }}{A}_r\sin (2{\mathrm{\pi rf}}_0{\mathrm{nT}}_s+{\phi }_r),n=0,1,...,N-1---\left(1\right)$

式中：R为谐波最高次数，f₀为基波频率，A_r和φ_r分别为第r次谐波的幅值和初始相位，T_s为数据采集系统的采样周期，采样频率1/T_s满足奈奎斯特抽样定理。

(2)根据测量要求，选取合适的组合余弦窗对离散采样序列x(n)加权，得到N点长序 列h_N(n)，其中，N表示采样点数。

用长度为N的余弦窗函数g(n)对离散采样序列x(n)加权，得到序列h_N(n)， h_N(n)＝x(n)·g(n)。

(3)对序列h_N(n)进行快速傅里叶变换(fastFouriertransform,FFT)，然后对其计算结果 取模求绝对值，获得包括基波在内的各次谐波幅值。

对h_N(n)进行FFT变换，根据傅里叶变换频域卷积定理，第r次谐波的离散频谱H(λ)为：

$H\left(\lambda \right)=\frac{A_r}{2}\underset{r=1}{\overset{R}{\Sigma }}[G(\lambda -{\lambda }_r)e^{-j[\pi (\lambda -{\lambda }_r)-{\phi }_r]}+G(\lambda +{\lambda }_r)e^{-j[\pi (\lambda +{\lambda }_r)-{\phi }_r]}]---\left(2\right)$

式中：G(λ)为余弦窗函数g(n)的频域函数；λ_r为第r次谐波频率被频率分辨率归一化 后的值，λ为被频率分辨率归一化的频率，

记频率分辨率为F₀，F₀＝1/NT_s，

λ_r＝rf₀/F₀。

在不计谐波间相互干扰及负频次谐波影响情况下，式(2)可整理为：

$H\left(\lambda \right)=\frac{A_r}{2}G\left(\sigma \right)e^{-j(\pi \sigma -{\phi }_r)}---\left(3\right)$

式中，σ＝λ-λ_r为频率偏差。根据式(3)，第r次谐波幅值A_r表示为：

$A_r=\frac{2\left|H\left(\lambda \right)\right|}{\left|G\left(\sigma \right)\right|}---\left(4\right).$

(4)以谱峰搜索的形式寻找所测谐波附近最大的两条谱线，并计算他们的谱值比γ，

γ＝次最大谱线值/最大谱线值。

(5)根据牛顿插值方法获得频率偏差σ，

由于在数据采样时频率分辨率满足各次谐波频率处于不同的频率间隔中，故频率偏差 小于半个频率分辨率，将频率用频率分辨率归一化。在组合余弦窗从中心位置到半个归一化 频率[0,0.5]区间内，将频率均分为P等份(P≤10)，则均分点处对应的频率值λ_i等于i＝0,1...P。

频域采样时，如果一个采样点落在λ_i处，由于相邻采样点之间的差值为1，则在λ_i减 去1处，也会有一个采样点，记为ρ_i，则ρ_i＝λ_i-1。λ_i、ρ_i及余弦窗频域函数G(λ)(λ为 频率分辨率归一化频率)之间的关系如图1所示。

计算ρ_i和λ_i处的主瓣幅值G(ρ_i)、G(λ_i)，求两者的比值并将其存储在内存中。

这里，计算主瓣幅值的方法如下：

组合余弦窗的离散时域表达式g(n)为：

$g\left(n\right)=\underset{m=0}{\overset{M-1}{\Sigma }}{(-1)}^m{a}_{m}cos\left(\frac{2\pi }{N}mn\right),n=0,1,2...N-1---\left(5\right)$

其中，M为窗函数的项数，a_m为窗函数系数。

式(5)经离散傅里叶变换，得到余弦窗频域函数G(λ)：

$G\left(\lambda \right)=\underset{m=0}{\overset{M-1}{\Sigma }}{(-1)}^m\frac{a_m}{2}[D(\lambda -m)+D(\lambda +m)]---\left(6\right)$

式中，λ为被频率分辨率归一化的频率，

把λ_i和ρ_i的值作为自变量λ代入式(6)求绝对值，即可获得对应频率处的主瓣幅值。

计算ρ_i对应的主瓣幅值与λ_i处主瓣幅值的比值，该比值为一组离散数据，用ζ_i表示：

${\zeta }_i=\frac{\left|G({\lambda }_i-1)\right|}{\left|G\left({\lambda }_i\right)\right|}=\frac{\left|G(i·\frac{0.5}{P}-1)\right|}{\left|G(i·\frac{0.5}{P})\right|},i=0,1.P.---\left(7\right)$

以ζ_i为自变量，λ_i为因变量，构建牛顿插值多项式f(x)如下：

f(x)＝f[x₀]+f[x₀,x₁](x-x₀)

+f[x₀,x₁,x₂](x-x₀)(x-x₁)+···(8)

+f[x₀,x₁,···,x_i](x-x₀)···(x-x_i-1)

式中，

$f[x_0,x_1,...,x_i]=\frac{f[x_0,x_1···,x_{i-2},x_i]-f[x_0,x_1...,x_{i-1}]}{x_i-x_{i-1}}$

其中，x₀,x₁,···,x_i分别对应ζ₀,ζ₁,···,ζ_i的计算值，f(x₀),f(x₁),···,f(x_i)分别对应 λ₀,λ₁,···,λ_i的数值。

假设被测信号是第r次谐波，非整周期采样下，频域内接近第r次谐波的次最大谱线与 最大谱线之比为γ。显然，γ的值介于数组ζ_i中最小值和最大值之间，故可以用上述构造的 牛顿插值多项式，计算谱值比为γ时对应的频率偏差σ。

(6)根据求得的频率偏差，进一步求解该次谐波的幅值、频率和初相位。

在步骤(3)，已经得到第r次谐波幅值A_r表示为式(4)：

式中的σ为频率偏差，由步骤(5)求得，则将σ作为自变量λ代入式(6)求绝对值，即可获 得|G(σ)|的值，|H(λ)|为最接近第r次谐波离散频谱H(λ)的幅值，也为已知，故可以根据(4) 求得第r次谐波幅值A_r。

设频域轴上的第p根谱线最接近第r次谐波谱线。则第r次谐波频率f_r可表示为：

$f_r={\mathrm{rf}}_0=\frac{(p-1-\sigma )}{T_{s}N}---\left(9\right)$

最后，由式(3)可知，第r次谐波的初相位φ_r为：

φ_r＝φ_p+πσ(10)

式中，φ_p为频域轴上最接近第r次谐波的第p根谱线的初相角。

以上即为该频域插值方法求解电力谐波参数的过程。

为了进一步说明本发明的具体实施方式，以四项Blackman-Harris窗为例，选用该窗对采 样数据加权，然后用本发明所提频域插值方法进行谐波分析。

将四项Blackman-Harris窗在归一化频率[0,0.5]区间内，均分为10份，则各均分点处频率 为λ_i＝0.05·i(i＝0,1,···,10)，分别计算λ_i和λ_i-1处的主瓣幅值；然后计算i取不同值时， λ_i-1与λ_i对应主瓣幅值的比值ζ_i。其主瓣幅值和比值如表1所示。将λ_i和ζ_i以内存形式存 入内存。以ζ_i为自变量，λ_i为因变量，编写牛顿插值程序。

表1为四项Blackman-Harris窗归一化频率[0,0.5]区间内，均分为10份时，各均分点处频 率为λ_i＝0.05·i(i＝0,1,···,10)与(λ_i-1)处的主瓣幅值，及两者主瓣幅值之比。

表1

从表1可知，λ_i和ζ_i都是单调递增的，故可以用插值方法求解自变量的某个内插值点对 应的函数值。以电力谐波为例，频域内某次谐波附近，次最大谱线与最大谱线的幅值比为 0.8849，把该值代入以ζ_i为自变量λ_i为因变量的牛顿插值多项式程序，求得频率偏差为 0.3400。然后，根据该频率偏差就可以求得该次谐波的幅值、频率和初相位。

实际中，任意次数谐波附近的次最大谱线与最大谱线之比在[0.6805,1]之间任意取值时， 不考虑随机噪声和间谐波等影响，及四项Blackman-Harris窗在归一化频率[0,0.5]区间内均分 为10等份条件下，由本发明方法计算得到频率偏差与理论值之间的误差如图2所示，稳定 在10^-13数量级，测量误差非常小。