基于无迹卡尔曼滤波的频率与谐波检测方法与流程

本发明涉及电力系统电能质量检测领域，特别是一种基于无迹卡尔曼滤波的频率与谐波检测方法。

背景技术：

电力系统中频率的测量是其他各项电能质量参数测量的基础。以谐波测量为例，各次谐波的频率与基波的频率都是倍数关系，如果要准确测量各次谐波的幅值和相位，基波频率必须为已知，即需要提前测量。

测频方法有过零点法，检测波形两个过零点之间的时间间隔来计算频率，但谐波严重时过零点检测不准确导致频率测量不准确，同时器件零点漂移的问题也会给测量带来误差。插值周期法对过零点附近进行插值处理，降低了过零点检测难度，但采样不同步会导致线性化误差。傅里叶变换法既可以用来测频率也可以用来测谐波，但是这个方法需要获得一个周期的数据之后才能计算，实时性不好，且存在因采样周期与波形周期不严格同步导致的频谱泄露和栅栏效应的问题。传统的卡尔曼滤波测谐波的算法，状态变量大都选取为不能从状态变量中直接得到幅值和相位，需要进一步计算。

技术实现要素：

本发明的目的在于提供一种基于无迹卡尔曼滤波的频率与谐波检测方法。

实现本发明目的的技术方案为：一种基于无迹卡尔曼滤波的频率与谐波检测方法，包括以下步骤：

步骤1，建立含谐波的信号模型：

其中，y是观测信号，r是谐波次数，m是总谐波次数，f是基波频率，ar是r次谐波幅值，是r次谐波相位，e是协方差为r的零均值高斯白噪声；

步骤2，选取观测状态变量：

其中，ar,k是k时刻r次谐波的幅值，是k时刻r次谐波的相位，fk是k时刻基波频率；

步骤3，建立系统动态方程和量测方程：

其中，

ηk和ek是过程噪声和观测噪声，方差分别为qk和rk，ts为采样间隔；

步骤4，进行无迹卡尔曼滤波，得到状态变量的估计值：

其中，是算法得到的k时刻r次谐波的幅值，是算法得到的k时刻r次谐波的相位，r＝1,2,…,m，是算法得到的k时刻基波频率。

与现有技术相比，本发明的显著优点为：本发明的检测方法在建立谐波检测模型时兼顾了频率检测问题，不再需要单独的测频算法，避免了遇到现有测频方法存在的一些问题；本发明的检测方法可以直接测得各次谐波幅值和相位信息，无需再进一步计算；本发明的检测方法避免了傅里叶变换法采样不同步导致的频谱泄露和栅栏效应的问题。

附图说明

图1是本发明使用的算法流程图。

图2是本发明实施例中基波频率检测结果图。

图3是本发明实施例中各次谐波幅值检测结果图。

图4是本发明实施例中各次谐波相位检测结果图。

具体实施方式

如图1所示，一种基于无迹卡尔曼滤波的频率与谐波检测方法，包括以下步骤：

步骤1，建立含谐波的信号模型：

其中，y是观测信号，r是谐波次数，m是总谐波次数，f是基波频率，ar是r次谐波幅值，是r次谐波相位，e是协方差为r的零均值高斯白噪声；

步骤2，选取观测状态变量：

每一特定次谐波的幅值和相位均对应两个分量，基波频率对应一个分量，其中，ar,k是k时刻r次谐波的幅值，是k时刻r次谐波的相位，fk是k时刻基波频率；

步骤3，建立系统动态方程和量测方程：

其中，

ηk和ek是过程噪声和观测噪声，方差分别为qk和rk，ts为采样间隔；

步骤4，进行无迹卡尔曼滤波，得到状态变量的估计值：

其中，是算法得到的k时刻r次谐波的幅值，是算法得到的k时刻r次谐波的相位，r＝1,2,…,m，是算法得到的k时刻基波频率。这些幅值、相位和频率信息作为检测方法的最终结果。

本发明通过选取状态变量可一次性直接测得频率及各次谐波的幅值和相位，不再需要单独的测频算法，也无需进一步计算幅值和相位，避免了遇到上述频率和谐波检测方法存在的问题。使用无迹卡尔曼滤波不需要严格控制每个波形周期内采样点的个数，即没有上述傅里叶变换法采样不同步导致的频谱泄露和栅栏效应的问题。

下面结合实施例对本发明进行详细说明。

实施例

本部分以检测系统输出波形中3、5次谐波分量为例，详细说明本发明的实施方式。

步骤1，根据算法设定的模型，在仿真中预先生成一个含3、5次谐波分量的表达式：

其中，a1＝1，a3＝0.23，a5＝0.13，ω＝2πf＝2π×50.1；

然后在这些幅值和相位及频率中加入一定的高斯白噪声信号，并离散化，用来模拟系统状态：

其中，ηk～n(0,qk)，qk＝0.001*[0.04,0.001,0.01,0.0004,0.002,0.0001,0.01]^t；

最后在上述i信号中加入信噪比为35db的高斯白噪声，并离散化，用来模拟互感器采集到的信号：

其中，ts＝0.0002，ek用来对应加入的高斯白噪声信号；

步骤2，选取观测状态变量：

上述状态变量包含基波和3、5次谐波的幅值相位信息，也包含基波频率信息；

步骤3，建立系统动态方程和量测方程：

其中，

步骤4，进行无迹卡尔曼滤波，主要分为两步，预测、更新，具体如下：

下面计算中的rk＝0.02；

状态方程预测

sigma点的计算：

表示(n+λ)pk-1平方根的第i列，n＝7为维数，λ＝α²(n+k)-n，α＝1为比例修正因子，κ＝2为高阶状态分布的先验知识；

预测：

观测方程预测

sigma点的计算：

表示(n+λ)pk|k-1平方根的第i列；

预测：

状态更新

算法计算结果如图2、图3、图4所示，分别为基波频率检测结果、各次谐波幅值检测结果和各次谐波相位检测结果。

本发明的检测方法在建立谐波检测模型时兼顾了频率检测问题，不再需要单独的测频算法，避免了遇到现有测频方法存在的一些问题；本发明的检测方法可以直接测得各次谐波幅值和相位信息，无需再进一步计算；本发明的检测方法避免了傅里叶变换法采样不同步导致的频谱泄露和栅栏效应的问题。

实施例及附图是对本发明的功能进行描述的，并不是用来限定本发明的，任何在本发明的基础上所作的等效变换，均包括在本发明的专利保护范围内。