一种电网频率估计方法与流程

本发明涉及电力系统运行与控制技术领域，特别是一种电网频率估计方法。

背景技术：

频率是电力系统运行最基本的参数，准确快速的频率估计对电力系统运行、监测和控制具有重要的应用价值。目前电力系统频率估计方法以快速傅里叶变换(fastfouriertransform，fft)为主，并对频谱谱线进行两点或三点插值，从而获得信号的频率估计值。

现有技术的不足有：

(1)电力系统频率总是处于波动之中，无法始终保持采样频率为50hz(基波频率，国外为60hz)的整数倍，对采样数据进行周期延拓势必在数据窗的首端、末端出现突变，产生频率泄漏，使得fft结果不能真实反映信号的频谱；

(2)fft为频域分析方法，待分析信号必须满足一定的长度，fft分析结果为一定长度信号频率的平均值，这也意味着基于fft的频率估计算法存在一定的时延。

技术实现要素：

有鉴于现有技术的上述缺陷，本发明的目的就是提供一种电网频率估计方法，该方法为时域估计法，不借助数学优化算法或迭代过程即可估计电网频率，所需的采样点数目少、运算速度快、响应时间短，适合短时间窗条件下的频率估计，亦可推广用于正弦信号的频率估计。

本发明的目的是通过这样的技术方案实现的，一种电网频率估计方法，它包括有：

s1：任意给定一组正弦信号采样序列x(n)，扫描序列x(n)，将异常值用其左右两个采样点的平均值替代；

s2：将x(n)划分为一段段子序列，子序列的长度为n；

s3：任选一段子序列，采用fft得出待估计信号的基波频率初始值；

s4：根据ω＝f/fs计算圆周频率ω，其中f为步骤s2中获得的基波频率初始值；

s5：根据ωk＝π/3计算采样间隔k，根据采样间隔k，生成a(n)和b(n)，其中，

a(n)＝x(n)+x(n-k)(1)；

b(n)＝[x(n-2k)+x(n-k)+x(n)+x(n+k)]/2(2)；

s6：计算子序列频率估计值

s7：重复步骤s5和s6，直到估计出所有子序列的频率。

进一步，所述步骤s6中计算子序列频率估计值公式如下：

其中，n为子序列的长度，子序列的长度n应大于计算采样间隔k。

进一步，所述步骤s6中计算子序列频率估计值还包括有建立频率与采样点之间的方程组，建立过程如下：

设一理想电压信号采样表达式为：

式中，a、f、fs和θ分别表示电压的幅值、频率、采样率和初相，ω为圆周频率,n＝1,2,…,n；

考察s(n)相邻两个采样点的和有：

s(n+1)+s(n)＝acos[ω(n+1)+θ]+acos[ωn+θ]

＝2acos(ωn+θ+ω/2)cos(ω/2)(5)；

s(n)+s(n-1)＝acos[ωn+θ]+acos[ω(n-1)+θ]

＝2acos(ωn+θ-ω/2)cos(ω/2)(6)；

s(n-1)+s(n-2)＝acos[ω(n-1)+θ]+acos[ω(n-2)+θ]

＝2acos(ωn+θ-3ω/2)cos(ω/2)(7)；

由式(5)至(7)可知，正弦信号相邻采样点之和组成的信号仍为同频率的正弦信号，仅仅幅值和相位与原信号不同；由三角函数和差化积公式可得：

s(n+1)+s(n)+s(n-1)+s(n-2)

＝2×2acos(ωn+θ-ω/2)cos(ω/2)cosω

＝2×[s(n)+s(n-1)]cosω(8)；

即，

式(9)表明理想正弦信号的频率可通过四个相邻的采样点准确给出，若四个采样点不相邻，则该式为

或

s(n+k)+s(n)＝2cosωk×[s(n)+s(n-k)]-[s(n-2k)+s(n-k)](11)；

根据式(10)或(11)即可建立频率与采样点之间的方程组。

进一步，根据建立的频率与采样点之间的方程组计算子序列频率估计值具体过程如下：

由于式(10)或(11)满足ax＝b的形式，故可根据最小二乘法的原理定义误差e，如下：

其中，a(n)＝x(n)+x(n-k)，

和b(n)＝[x(n-2k)+x(n-k)+x(n)+x(n+k)]/2以及n为数据窗的长度；

误差e是圆周频率ω的一元函数，对ω求导并令其为零得：

式中，为频率测量值。

进一步，还包括有对噪声的抑制，抑制过程如下：

首先考虑加性高斯白噪声对频率估计的影响，假设该噪声序列w(n)的均值为零、方差为σ²，则

x(n)＝s(n)+w(n)＝acos(ωn+θ)+w(n)(14)；

将x(n-2k)、x(n-k)、x(n)和x(n+k)代入a(n)和b(n)，化简式(13)的分子和分母：

由此可见，当数据窗长度n较大时，这些求和项同式(15)和(16)等号右边的前两项系数n-3k相比可以忽略不计，故式(13)可近似为：

式(17)给出圆周频率的估计值与真实值ω之间的近似关系表达式，其中a²/σ²表示测量信号的信噪比snr。

进一步，还包括有对谐波的处理，处理过程如下：

考虑谐波的信号模型为：

式中，ai和θi分别为第i次谐波的幅值和初始相位

将x(n-2k)、x(n-k)、x(n)和x(n+k)代入a(n)和b(n)，化简式(13)的分子和分母得：

考虑谐波影响的频率估计表达式为：

简化式(21)有：

由此可见，谐波与基波幅值之比的平方影响频率估计的精度。

进一步，给定数据窗长度n，满足条件的k值有n/3个，因此总可能存在某个k值，使得：

cos3ωk×[cos3ωk+1]＝0(23)

或，cos5ωk×[cos5ωk+1]＝0(24)。

由于采用了上述技术方案，本发明具有如下的优点：本发明的目的是提出一种快速准确的电力系统频率估计方法，该方法引入采样点间隔参数k，采用四个间距为k的采样点建立电压采样点与频率之间的数学表达式，根据最小二乘原理得出频率估计的表达式，分别考虑白噪声和谐波的影响，推导出频率估计的近似显性表达式。

本发明通过正弦信号采样值之间的三角函数关系推导出频率估计表达式，避免了现有fft方法频率泄漏问题和对数据窗长度的要求，能对短时间范围内依靠四个采样点即可完成频率估计，计算速度快，测量的响应时间短，且具有最佳的白噪声和三次谐波抑制能力。

采样点间隔参数影响频率估计精度，当信号的圆周频率与采样点间隔参数的乘积取π/3时，算法具有最佳的白噪声和三次谐波抑制能力。仿真得出的最佳采样点间隔和理论值一致，该方法为时域估计法，不借助数学优化算法或迭代过程即可估计电网频率，所需的采样点数目少、运算速度快、响应时间短，适合短时间窗条件下的频率估计，亦可推广用于正弦信号的频率估计。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。

附图说明

本发明的附图说明如下：

图1为电网频率估计方法的流程示意图。

图2为当k＝1、snr＝10db时，cos(2ωn-ωk+2θ)的仿真求和结果。

图3为当k＝1、snr＝10db时，cos(ωn-ωk/2+θ)的仿真求和结果。

图4为当k＝1、snr＝10db时，w(n)+w(n-k)的仿真求和结果。

图5为当k＝1、snr＝10db时，w(n+k)+w(n)+w(n-k)+w(n-2k)的仿真求和结果。

图6为采样点间隔k对不同噪声的抑制性能。

图7为采样点间隔k对不同谐波干扰的抑制性能。

图8为当采样间隔点k＝4，snr＝30db时，不同数据窗长度对频率估计结果的影响。

图9为当采样间隔点k＝4，snr＝10db时，不同数据窗长度对频率估计结果的影响。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

实施例，如图1至图9所示；一种电网频率估计方法，它包括有：

s1：任意给定一组正弦信号采样序列x(n)，扫描序列x(n)，将异常值用其左右两个采样点的平均值替代；

s2：将x(n)划分为一段段子序列，子序列的长度为n；

s3：任选一段子序列，采用fft得出待估计信号的基波频率初始值；

s4：根据ω＝f/fs计算圆周频率ω，其中f为步骤s2中获得的基波频率初始值；

s5：根据ωk＝π/3计算采样间隔k，根据采样间隔k，生成a(n)和b(n)，其中，

a(n)＝x(n)+x(n-k)(1)；

b(n)＝[x(n-2k)+x(n-k)+x(n)+x(n+k)]/2(2)；

s6：计算子序列频率估计值

s7：重复步骤s5和s6，直到估计出所有子序列的频率。

所述步骤s6中计算子序列频率估计值公式如下：

其中，n为子序列的长度，子序列的长度n应大于计算采样间隔k。

所述步骤s6中计算子序列频率估计值还包括有建立频率与采样点之间的方程组，建立过程如下：

设一理想电压信号采样表达式为：

式中，a、f、fs和θ分别表示电压的幅值、频率、采样率和初相，ω为圆周频率,n＝1,2,…,n；

考察s(n)相邻两个采样点的和有：

s(n+1)+s(n)＝acos[ω(n+1)+θ]+acos[ωn+θ]

＝2acos(ωn+θ+ω/2)cos(ω/2)(5)；

s(n)+s(n-1)＝acos[ωn+θ]+acos[ω(n-1)+θ]

＝2acos(ωn+θ-ω/2)cos(ω/2)(6)；

s(n-1)+s(n-2)＝acos[ω(n-1)+θ]+acos[ω(n-2)+θ]

＝2acos(ωn+θ-3ω/2)cos(ω/2)(7)；

由式(5)至(7)可知，正弦信号相邻采样点之和组成的信号仍为同频率的正弦信号，仅仅幅值和相位与原信号不同；由三角函数和差化积公式可得：

s(n+1)+s(n)+s(n-1)+s(n-2)

＝2×2acos(ωn+θ-ω/2)cos(ω/2)cosω

＝2×[s(n)+s(n-1)]cosω(8)；

即，

式(9)表明理想正弦信号的频率可通过四个相邻的采样点准确给出，若四个采样点不相邻，则该式为：

或

s(n+k)+s(n)＝2cosωk×[s(n)+s(n-k)]-[s(n-2k)+s(n-k)](11)；

对数据窗内的一组采样点序列，根据式(10)或(11)即可建立频率与采样点之间的方程组。

然而在实际测量中往往无法得到理想的电压采样序列s(n)，实际采集的电压信号x(n)总存在一定的噪声或谐波干扰，x(n)与s(n)之间存在一定的偏差，用x(n)代替s(n)建立上述方程或方程组，不可避免地引入一定的估计误差。

由于式(10)或(11)满足ax＝b的形式，故可根据最小二乘法的原理定义误差e，如下：

其中，a(n)＝x(n)+x(n-k)，

和b(n)＝[x(n-2k)+x(n-k)+x(n)+x(n+k)]/2以及n为数据窗的长度；

误差e是圆周频率ω的一元函数，对ω求导并令其为零得：

式中，为频率测量值。

误差e包含n＝2k+1～n-k共n-3k项，这是为了兼顾采样点x(n-2k)和x(n+k)在n点数据窗之内始终有值与之对应，理论上式(13)的测量时延为(3k+1)/fs。

a(n)和b(n)根据采样点序列即可得到，但求和的范围由k决定，因此可以推测该值影响频率的估计值当分子和分母的n-3k项求和结束后，由反余弦函数即可得到

还包括有对噪声的抑制，抑制过程如下：

首先考虑加性高斯白噪声对频率估计的影响，假设该噪声序列w(n)的均值为零、方差为σ²，则

x(n)＝s(n)+w(n)＝acos(ωn+θ)+w(n)(14)；

将x(n-2k)、x(n-k)、x(n)和x(n+k)代入a(n)和b(n)，化简式(13)的分子和分母：

从上式可见，当x(n)为理想信号时，式(13)退化为式(10)；考虑白噪声之后，式(13)的分子和分母有一定数目的误差项。

图2至图5分别给出一个周期后当k＝1、snr＝10db时，cos2(ωn-ωk+2θ)、cos(ωn-ωk/2+θ)、w(n)+w(n-k)和w(n+k)+w(n)+w(n-k)+w(n-2k)的仿真求和结果。

由此可见，当数据窗长度n较大时，这些求和项同式(15)和(16)等号右边的前两项系数n-3k相比可以忽略不计，故式(13)可近似为：

式(17)给出圆周频率的估计值与真实值ω之间的近似关系表达式，其中a2/σ²表示测量信号的信噪比snr。

还包括有对谐波的处理，处理过程如下：

考虑谐波的信号模型为：

式中，ai和θi分别为第i次谐波的幅值和初始相位。

通常谐波的幅值比基波幅值小且电力系统的谐波以奇次谐波为主。

将x(n-2k)、x(n-k)、x(n)和x(n+k)代入a(n)和b(n)，化简式(13)的分子和分母得：

采用类似于式(15)和(16)的近似处理，式(19)和(20)第一项之后的误差项均可忽略，则考虑谐波影响的频率估计表达式为：

若仅考虑三次、五次谐波，当数据窗长度n较大时，式(19)等号右边第一项之后所有的正余弦乘积等求和项可忽略不计，则简化式(21)有：

由此可见，谐波与基波幅值之比的平方影响频率估计的精度。

若不考虑谐波分量，则式(22)退化为式(10)。给定数据窗长度n，满足条件的k值有n/3个，因此总可能存在某个k值，使得

cos3ωk×[cos3ωk+1]＝0(23)

或，cos5ωk×[cos5ωk+1]＝0(24)。

若以消除三次谐波为主，则ωk的取值可能为π/6、π/2或π/3。特别地，当ωk取π/3时，式(17)分子和分母中的噪声分量可同时消除；此时式(13)不仅能消除三次谐波，而且能消除白噪声干扰。以60.1hz的2秒电压仿真信号为例，若采样率固定为1440hz，则当k＝4，即相邻采样点的间隔为4，估计算法具有最佳的噪声和三次谐波抑制能力。

本发明的具体频率估计效果比对，给定一组电压采样点序列x(n)，则a(n)和b(n)已知，频率估计由n-3k项求和给出；当x(n)的长度n较大，满足k＜n/3条件的起始点较多，故选择不同的k值将影响算法对噪声和谐波的抑制效果。

①采样点间隔系数

参考fdr算法的计算流程，数据窗长度取8个周期，即192个采样点，信号长度为2秒；数据窗逐点移动，算法依次给出该数据窗对应的频率估计值。

图6和图7为考虑白噪声和谐波的影响，给出当取不同值时，频率估计的相对误差曲线。

图6和图7为频率估计相对误差的绝对值的半对数图，在图6中，k值取4时频率估计误差最小；图7采用现有技术(c.i.chen,g.w.chang,r.c.hongetal.extendedrealmodelofkalmanfilterfortime-varyingharmonicsestimation[j].ieeetrans.onpowerdelivery,2010,25(1):17-26.)的谐波信号，其中五、七次谐波的幅值固定为0.2和0.1，而三次谐波的幅值则线性增加，用它的幅值与基波幅值之比表示不同的谐波干扰水平。图中可见，在不同谐波干扰下，k＝4对谐波抑制的整体效果最好，该仿真结果和理论分析的结果相符合。

电网电压信号的频率随时间变化，因此ω＝2πf/fs在π/12左右变化，最佳的采样点间隔k有可能取3或5。

②数据窗长度n

频率估计公式在推导过程中假设数据窗长度n足够大，以确保n-3k远远大于其他误差项的系数。在实际测量中，n决定算法的响应时间，一般要求算法在确保精度的同时，具有最小的时延。

图8和图9给出当采样间隔点k＝4时，不同数据窗长度对频率估计结果的影响。仿真结果表明，数据窗长度对谐波抑制影响不大，其估计的相对误差稳定在10^-4-10^-3之间；当噪声干扰较强时，若选择稍长的数据窗，则噪声环境下的相对误差可控制在10-⁴以下。当数据窗长度取2～3个周期时，算法具有较好的噪声及谐波抑制效果，这意味着算法在该条件下的响应时间为0.033～0.05秒，小于fdr算法0.1～0.133秒的时延，该算法能在局部范围内估计频率，所需采样点数目小、运算速度快。

本发明通过正弦信号采样值之间的三角函数关系推导出频率估计表达式，避免了现有fft方法频率泄漏问题和对数据窗长度的要求，能对短时间范围内依靠四个采样点即可完成频率估计，计算速度快，测量的响应时间短，且具有最佳的白噪声和三次谐波抑制能力。

本发明的目的是提出一种快速准确的电力系统频率估计方法，该方法引入采样点间隔参数k，采用四个间距为k的采样点建立电压采样点与频率之间的数学表达式，根据最小二乘原理得出频率估计的表达式，分别考虑白噪声和谐波的影响，推导出频率估计的近似显性表达式。采样点间隔参数影响频率估计精度，当信号的圆周频率与采样点间隔参数的乘积取π/3时，算法具有最佳的白噪声和三次谐波抑制能力。仿真得出的最佳采样点间隔和理论值一致，该方法为时域估计法，不借助数学优化算法或迭代过程即可估计电网频率，所需的采样点数目少、运算速度快、响应时间短，适合短时间窗条件下的频率估计，亦可推广用于正弦信号的频率估计。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

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