CN103245831A - 一种基于广义卡尔曼滤波的谐波辨识方法 - Google Patents

Abstract

本发明涉及的是一种谐波辨识方法领域，具体涉及一种用于电力系统谐波辨识或者控制系统中的谐波辨识方法。本发明包括：建立响应信号非线性状态空间模型；提取k时刻状态向量估计值

提取出k时刻的谐波幅值和相位。本发明所提供的一种基于广义卡尔曼滤波的谐波辨识方法，能够快速、实时地得到简谐波激励响应信号中谐波成分的准确信息，从而直接为谐波成分的消除提供数据依据。本发明所提供的一种基于广义卡尔曼滤波的谐波辨识方法可以方便地应用于电力系统谐波辨识或者控制系统的谐波辨识中。

Description

一种基于广义卡尔曼滤波的谐波辨识方法

技术领域

本发明涉及的是一种谐波辨识方法领域，具体涉及一种用于电力系统谐波辨识或者控制系统中的谐波辨识方法。

背景技术

在电力系统中，谐波是负荷非线性所产生的电压电流的畸变，谐波消耗无功储备、影响电能计量，增加电能传输和利用的损耗，使设备过热，加剧设备的老化并缩短使用其使用寿命。产生的谐波谐振会造成电器元件的故障损毁，严重危害着电气环境。在信号干扰方面，谐波会产生电磁干扰，降低通信系统的通讯质量，造成重要元件和敏感元件的误动作，扰乱自动控制装置正常运行，使可靠性下降。谐波已成为电力公害。为了消除谐波成分，需要首先将谐波成分的信息辨识出来，以作为谐波消除的依据。

从查到的专利和文献的情况来看，已有使用卡尔曼滤波作为谐波辨识方法的例子，但他们将信号状态空间假设为简单的线性情况，使用线性形式的信号模型，以标准卡尔曼滤波作辨识算法，没有考虑实际信号的非线性特性。这样的辨识方法在卡尔曼滤波计算后还需要进行解算才能得到作为辨识结果的各次谐波幅值和相位，实时性不佳。

期刊文章“基于扩展卡尔曼滤波器的异步电机转速辨识”(李剑飞，尹泉，万淑芸.基于扩展卡尔曼滤波器的异步电机转速辨识[J].电工技术学报.2002,10,17(5).40-44.)采用扩展卡尔曼滤波来辨识异步电机的感应电流的方法来间接地辨识电机转速，所解决的实际问题和采用的方法及创新点均与本发明的谐波辨识存在巨大差异，是基于异步电机模型进行的。而硕士学位论文“基于卡尔曼滤波算法的动态谐波状态估计技术研究”(祝石厚.基于卡尔曼滤波算法的动态谐波状态估计技术研究[D].重庆大学2008)是以线性模型和线性卡尔曼滤波作为研究对象的，与本发明差异巨大。在“中国搜索专利数据库”（http://www.soopat.com/）和“中国知识产权数据库”(http://www.cnipr.com/)中检索到几篇相关专利。发明专利“一种基于卡尔曼滤波器的电液伺服系统波形再现控制方法”和发明专利“一种基于卡尔曼滤波器的电液伺服系统随机振动控制方法”都是卡尔曼滤波器在电液伺服系统中的应用，但是并没有将卡尔曼滤波作为辨识算法进行谐波辨识，与本发明有本质的区别。发明专利“系统模型未知的自适应卡尔曼滤波方法”是一种特殊的卡尔曼滤波方法，没有涉及在谐波辨识问题中的具体应用，与本发明的特点没有冲突。

“Power system harmonic analysis using the Kalman filter”(Karen Kennedy,GordonLightbody,Robert Yacamini.Power system harmonic analysis using the Kalman filter.PowerEngineering Society General Meeting,2003,IEEE.13-17July2003.Toronto,Canada.)也使用了扩展卡尔曼滤波来进行辨识，但响应信号模型的状态向量不是各次谐波的幅值和相位，与本专利有较大区别。“Self-Tuning of Kalman Filters for Harmonic Computation”(José A.RosendoMacías,Senior Member,IEEE,and Antonio Gómez Expósito,Fellow,IEEE.Self-Tuning ofKalman Filters for Harmonic Computation.IEEE TRANSACTIONS ON POWER DELIVERY,VOL.21,NO.1,JANUARY2006.501-503.)响应信号的模型是线性状态空间的形式，使用标准卡尔曼滤波作辨识算法，与本专利有较大差异。

发明内容

本发明的目的在于提供一种实现快速、实时地得到响应信号谐波成分信息谐波辨识方法。

本发明的目的是这样实现的：

本发明包括以下步骤：

(1)建立响应信号非线性状态空间模型：

k时刻的响应信号为n个频率为基频整数倍的谐波成分之和为：

$s\left(k\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i\left(k\right)\sin ({\mathrm{i\ω}}_{1}k+{\mathrm{\φ}}_i\left(k\right))$

式中k为时刻序列；n为谐波阶次；ω₁为基波圆频率ω₁=2πf₁/f_s，f₁、f_s为基波和系统采样频率；iω₁为第i次谐波的圆频率；φ_i(k)为k时刻第i次谐波的相角；C_i(k)为k时刻第i次谐波的幅值；

状态空间状态向量x(k)=[x₁(k) x₂(k) ... x_2n-1(k) x_2n(k)]^T的2n个元素为：

x₁(k)=C₁(k)，…x_2i-1(k)=C_i(k)，…x_2n-1(k)=C_n(k)，

x₂(k)=φ₁(k)，…x_2i(k)=φ_i(k)，…x_2n(k)=φ_n(k)，

噪声影响下包含n次谐波的响应信号状态方程和测量方程为：

x(k+1)=F(k,x(k))+w(k)=x(k)+w(k)，

z(k)=H(k,x(k))+v(k)

=(x₁(k)sin(ω₁k)cos(x₂(k))+x₁(k)cos(ω₁k)sin(x₂(k)))+…

(x_2i-1(k)sin(ω₁k)cos(x_2i(k))+x_2i-1(k)cos(ω₁k)sin(x_2i(k)))+…，

(x_2n-1(k)sin(ω₁k)cos(x_2n(k))+x_2n-1(k)cos(ω₁k)sin(x_2n(k)))+v(k)

式中z(k)为测量值，z(k)=s(k)，1×1矩阵；F(k,x(k))为k时刻的非线性状态转移函数；H(k,x(k))为k时刻的非线性测量函数；w(k)为状态转移过程噪声；v(k)为测量噪声；

(2)提取k时刻状态向量估计值

提取状态空间方程的雅克比矩阵F(k)和H(k)，状态转移函数F(k,x(k))的雅克比矩阵F(k)为单位矩阵，测量函数H(k,x(k))的雅克比矩阵H(k)为维数为1×2n的矩阵：

$H\left(k\right)={\left[\begin{array}{c}\sin \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)\\ -x_1\left(k\right)\sin \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)+x_1\left(k\right)\cos \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)\\ -x_{2i-1}\left(k\right)\sin \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)+x_{2i-1}\left(k\right)\cos \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)\\ -x_{2n-1}\left(k\right)\sin \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)+x_{2n-1}\left(k\right)\cos \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)\end{array}\right]}^T,$

提取k时刻状态向量估计值

提取k时刻一步预测状态向量

提取k时刻一步预测均方误差P^-(k) P^-(k)=A(k-1)P(k-1)A^T(k-1)+Q(k-1)，

提取k时刻卡尔曼滤波增益K(k) K(k)=P-(k)B^T(k)[B(k)P^-(k)B^T(k)+R(k)]^-1，

提取k时刻状态向量估计值

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)[z\left(k\right)-B\left(k\right){\hat{x}}^{-}\left(k\right)],$

提取k时刻估计均方误差P(k) P(k)=[I-K(k)B(k)]P^-(k)，

式中Q(k)为过程噪声w(k)的相关矩阵，R(k)为测量噪声v(k)的相关矩阵；

(3)提取出k时刻的谐波幅值和相位：

k时刻第i次谐波的幅值为：

C_i(k)=x_2i-1(k)；

相位为：

φ_i(k)=x_2i(k)。

响应信号为电网电流、电压或控制系统的反馈信号。

本发明的有益效果在于：

本发明所提供的一种基于广义卡尔曼滤波的谐波辨识方法，能够更快速、实时地得到简谐波激励响应信号中谐波成分的准确信息，从而直接为谐波成分的消除提供数据依据。本发明所提供的一种基于广义卡尔曼滤波的谐波辨识方法可以方便地应用于电力系统谐波辨识或者控制系统的谐波辨识中。

附图说明

图1为辨识方法的原理示意图；

图2为辨识方法对电液伺服振动台正弦输入信号为4sin(2π×5t)m/s²时的台面加速度响应信号的幅值辨识结果；

图3为辨识方法对电液伺服振动台正弦输入信号为4sin(2π×5t)m/s²时的台面加速度响应信号的相位辨识结果。

具体实施方式

参照附图对本发明作更详细的描述：

基于广义卡尔曼滤波的谐波辨识方法主要包括以下三个关键步骤：

(1)建立响应信号非线性状态空间模型：

k时刻的响应信号表示为n个频率为基频整数倍的谐波成分之和可以表示为：

$s\left(k\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i\left(k\right)\sin ({\mathrm{i\ω}}_{1}k+\mathrm{\φ}\left(k\right))$

式中k——时刻序列；

n——谐波阶次；

ω₁——基波圆频率ω₁=2πf₁/f_s（rad/s），f₁、f_s为基波和系统采样频率（Hz）；

iω₁——第i次谐波的圆频率（rad/s）；

φ_i(k)——k时刻第i次谐波的相角（rad）；

C_i(k)——k时刻第i次谐波的幅值（m/s²）。

状态空间状态向量x(k)=[x₁(k) x₂(k) ... x_2n-1(k) x_2n(k)]^T的2n个元素为

x₁(k)=C₁(k)，…x_2i-1(k)=C_i(k)，…x_2n-1(k)=C_n(k)

x₂(k)=φ₁(k)，…x_2i(k)=φ_i(k)，…x_2n(k)=φ_n(k)

噪声影响下包含n次谐波的响应信号状态方程和测量方程可表示为

x(k+1)=F(k,x(k))+w(k)=x(k)+w(k)

z(k)=H(k,x(k))+v(k)

=(x₁(k)sin(ω₁k)cos(x₂(k))+x₁(k)cos(ω₁k)sin(x₂(k)))+…

(x_2i-1(k)sin(ω₁k)cos(x_2i(k))+x_2i-1(k)cos(ω₁k)sin(x_2i(k)))+…

(x_2n-1(k)sin(ω₁k)cos(x_2n(k))+x_2n-1(k)cos(ω₁k)sin(x_2n(k)))+v(k)

式中z(k)——测量值，z(k)=s(k)，1×1矩阵；

F(k,x(k))——k时刻的非线性状态转移函数；

H(k,x(k))——k时刻的非线性测量函数；

w(k)——为状态转移过程噪声；

v(k)——测量噪声。

由于正弦振动试验频率即基频信号的频率ω₁已知，且已确定了响应信号需要辨识的谐波波次n，所以可以得到响应信号状态方程和测量方程。响应信号状态方程为线性形式，而测量方程为非线性形式。

(2)应用广义卡尔曼滤波来递推计算k时刻状态向量估计值

首先求取状态空间方程的雅克比矩阵F(k)和H(k)，使得卡尔曼滤波在线性的模型下进行。状态转移函数F(k,x(k))的雅克比矩阵F(k)为单位矩阵，测量函数H(k,x(k))的雅克比矩阵H(k)为维数为1×2n的矩阵。

$H\left(k\right)={\left[\begin{array}{c}\sin \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)\\ -x_1\left(k\right)\sin \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)+x_1\left(k\right)\cos \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)\\ -x_{2i-1}\left(k\right)\sin \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)+x_{2i-1}\left(k\right)\cos \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)\\ -x_{2n-1}\left(k\right)\sin \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)+x_{2n-1}\left(k\right)\cos \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)\end{array}\right]}^T$

然后应用如下的广义卡尔曼滤波方程来递推计算k时刻状态向量估计值

a)计算k时刻一步预测状态向量

b)计算k时刻一步预测均方误差P^-(k) P^-(k)=A(k-1)P(k-1)A^T(k-1)+Q(k-1)

c）计算k时刻卡尔曼滤波增益K(k) K(k)=P^-(k)B^T(k)[B(k)P^-(k)B^T(k)+R(k)]^-1

d）计算k时刻状态向量估计值

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)[z\left(k\right)-B\left(k\right){\hat{x}}^{-}\left(k\right)]$

e）计算k时刻估计均方误差P(k) P(k)=[I-K(k)B(k)]P^-(k)

式中Q(k)为过程噪声w(k)的相关矩阵，R(k)为测量噪声v(k)的相关矩阵。

(3)通过分解

提取出k时刻的谐波幅值和相位：

广义卡尔曼滤波计算出的状态向量直接对应着谐波的幅值和相位，因此在经过广义卡尔曼滤波算法估计出状态向量后，可以直接进行分解得到各次谐波的幅值和相位的值。辨识方法的原理示意如图1所示。k时刻第i次谐波的幅值和相位分别为

C_i(k)=x_2i-1(k)

φ_i(k)=x_2i(k)

分解

过程是分解

为奇数下标元素x_2i-1(k)序列和偶数下标元素序列x_2i(k)。将x_2i-1(k)序列提取出作为幅值辨识结果，而将x_2i(k)序列提取出作为相位辨识结果。

实施例：

电液伺服振动台控制系统的正弦输入信号为4sin(2π×5t)m/s²，此时振动台台面加速度响应信号出现了非常严重的失真现象。用本发明所提供的一种基于广义卡尔曼滤波的谐波辨识方法对该加速度响应信号作谐波辨识，能得到如图2所示的各次谐波幅值辨识结果和如图3所示的各次谐波相位辨识结果。

各次谐波的幅值辨识结果和相位辨识结果在2秒后基本都稳定下来，稳定后的幅值和相位波动幅度也很小。

Claims (2)

1.一种基于广义卡尔曼滤波的谐波辨识方法，其特征是，包括以下步骤：

(1)建立响应信号非线性状态空间模型：

k时刻的响应信号为n个频率为基频整数倍的谐波成分之和为：

$s\left(k\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{C}_i\left(k\right)\sin ({\mathrm{i\ω}}_{1}k+{\mathrm{\φ}}_i\left(k\right))$

式中k为时刻序列；n为谐波阶次；ω₁为基波圆频率ω₁=2πf₁/f_s，f₁、f_s为基波和系统采样频率；iω₁为第i次谐波的圆频率；φ_i(k)为k时刻第i次谐波的相角；C_i(k)为k时刻第i次谐波的幅值；

状态空间状态向量x(k)=[x₁(k) x₂(k) ... x_2n-1(k) x_2n(k)]^T的2n个元素为：

x₁(k)=C₁(k)，…x_2i-1(k)=C_i(k)，…x_2n-1(k)=C_n(k)，

x₂(k)=φ₁(k)，…x_2i(k)=φ_i(k)，…x_2n(k)=φ_n(k)，

噪声影响下包含n次谐波的响应信号状态方程和测量方程为：

x(k+1)=F(k,x(k))+w(k)=x(k)+w(k)，

z(k)=H(k,x(k))+v(k)

=(x₁(k)sin(ω₁k)cos(x₂(k))+x₁(k)cos(ω₁k)sin(x₂(k)))+…

(x_2i-1(k)sin(ω₁k)cos(x_2i(k))+x_2i-1(k)cos(ω₁k)sin(x_2i(k)))+…，

(x_2n-1(k)sin(ω₁k)cos(x_2n(k))+x_2n-1(k)cos(ω₁k)sin(x_2n(k)))+v(k)

式中z(k)为测量值，z(k)=s(k)，1×1矩阵；F(k,x(k))为k时刻的非线性状态转移函数；H(k,x(k))为k时刻的非线性测量函数；w(k)为状态转移过程噪声；v(k)为测量噪声；

(2)提取k时刻状态向量估计值

提取状态空间方程的雅克比矩阵F(k)和H(k)，状态转移函数F(k,x(k))的雅克比矩阵F(k)为单位矩阵，测量函数H(k,x(k))的雅克比矩阵H(k)为维数为1×2n的矩阵：

$H\left(k\right)={\left[\begin{array}{c}\sin \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)\\ -x_1\left(k\right)\sin \left({\mathrm{\ω}}_{1}k\right)\sin \left(x_2\left(k\right)\right)+x_1\left(k\right)\cos \left({\mathrm{\ω}}_{1}k\right)\cos \left(x_2\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)\\ -x_{2i-1}\left(k\right)\sin \left({\mathrm{\ω}}_{i}k\right)\sin \left(x_{2i}\left(k\right)\right)+x_{2i-1}\left(k\right)\cos \left({\mathrm{\ω}}_{i}k\right)\cos \left(x_{2i}\left(k\right)\right)\\ \·\·\·\\ \sin \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)+\cos \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)\\ -x_{2n-1}\left(k\right)\sin \left({\mathrm{\ω}}_{n}k\right)\sin \left(x_{2n}\left(k\right)\right)+x_{2n-1}\left(k\right)\cos \left({\mathrm{\ω}}_{n}k\right)\cos \left(x_{2n}\left(k\right)\right)\end{array}\right]}^T,$

提取k时刻状态向量估计值

提取k时刻一步预测状态向量

提取k时刻一步预测均方误差P^-(k) P^-(k)=A(k-1)P(k-1)A^T(k-1)+Q(k-1)，

提取k时刻卡尔曼滤波增益K(k) K(k)=P^-(k)B^T(k)[B(k)P^-(k)B^T(k)+R(k)]^-1，

提取k时刻状态向量估计值

$\hat{x}\left(k\right)={\hat{x}}^{-}\left(k\right)+K\left(k\right)[z\left(k\right)-B\left(k\right){\hat{x}}^{-}\left(k\right)],$

提取k时刻估计均方误差P(k) P(k)=[I-K(k)B(k)]P^-(k)，

式中Q(k)为过程噪声w(k)的相关矩阵，R(k)为测量噪声v(k)的相关矩阵；

(3)提取出k时刻的谐波幅值和相位：

k时刻第i次谐波的幅值为：

C_i(k)=x_2i-1(k)；

相位为：

φ_i(k)=x_2i(k)。

2.根据权利要求1所述的一种基于广义卡尔曼滤波的谐波辨识方法，其特征在于：所述的响应信号为电网电流、电压或控制系统的反馈信号。

Images