CN102222911A - 基于ar模型和卡尔曼滤波的电力系统间谐波估计方法 - Google Patents

Abstract

基于AR模型和卡尔曼滤波的电力系统间谐波估计方法，它涉及一种间谐波估计方法。本发明为解决现有的快速傅里叶变换方法实时性差以及现有的Burg算法存在谱分析性能低、抗干扰能力低的问题。步骤A、采集电力系统的信号数据；步骤B、由采集电力系统信号的数据建立AR模型，然后利用AR模型对电力系统信号进行分析；步骤C、建立状态方程和观测方程；步骤D、利用最终预测误差准则确定AR模型阶数；步骤E、利用自适应卡尔曼滤波方法对AR模型参数在线估计；步骤F、由AR模型进行功率谱估计：将最佳估计值

代入到AR模型功率谱密度公式中进行功率谱估计。本发明用于对电力系统间谐波的功率谱进行估计。

Description

基于AR模型和卡尔曼滤波的电力系统间谐波估计方法

技术领域

本发明涉及一种电力系统间谐波估计方法，具体涉及一种基于AR模型和卡尔曼滤波的电力系统间谐波估计方法。

背景技术

随着电力系统的迅速发展，大量的电力电子装置、非线性波动负荷以及各种变频调速装置产生的间谐波对电力系统的污染日益严重，并且威胁到电网中各种电气设备的安全经济运行。因此，对电力系统间谐波的分析也成为国内外学者广泛关注的热点，并且具有非常重要的实际意义。

常用的间谐波检测方法有快速傅里叶变换(FFT)、特征分解法、人工神经网络模型、Burg算法等。

快速傅里叶变换(FFT)：其计算量较小，实现较为简单而被广泛使用。但是间谐波分量很难确定其波形周期，且其幅值一般小于基波和谐波，传统的FFT算法就难以实现精确的同步采样，从而引起频谱泄露和栅栏效应等误差。加窗插值FFT虽然对传统FFT算法的效果有一定的改善，但是其分析窗的宽度一般要几十个信号周期，当间谐波与谐波频率相近时，要能准确地检测出间谐波分量，分析窗的宽度还需要进一步增加，则其检测实时性就也会大大降低。

特征分解法：包括Pisarenko(皮萨连科)谐波分析法、MUSIC(多信号分类法)和Prony(波朗尼)算法等，用于电网间谐波分析时，虽然结果较为精确，但需要进行自相关矩阵估计，运算量较大，且对噪声比较敏感，算法的实用性不强。

人工神经网络模型：其神经网络的训练可能要花很长时间，有时根本就不会收敛，且其算法复杂度较高，计算量太大，严重制约间谐波检测的实时性。

Burg算法(伯格算法)：其分辨率高，数值稳定性好，计算效率也较高，因此经常被应用于电力系统间谐波分析中。但由于Burg算法采用的Levinson(莱文森)递推公式这一强约束条件，使其出现谱峰偏移和频谱分裂现象，

综上，现有的快速傅里叶变换方法存在分析窗的宽度大导致实时性差；现有的Burg算法存在谱峰偏移和频谱分裂现象导致谱分析性能低、抗干扰能力低。

发明内容

本发明的目的是为解决现有的快速傅里叶变换方法实时性差以及现有的Burg算法存在谱分析性能低、抗干扰能力低的问题，进而提供了一种基于AR模型和卡尔曼滤波的电力系统间谐波估计方法。

本发明为解决上述技术问题采取的技术方案是：所述电力系统间谐波估计方法由以下步骤实现的：

步骤A、采集电力系统的信号数据，信号数据如下：

Y(n)＝[y(1)，y(2)，…y(N)]

其中N为信号数据的个数；

步骤B、由采集电力系统信号的数据建立AR模型，然后利用AR模型对电力系统信号进行分析；

步骤B1、采集电力系统信号的输入向量由y(n-1)、…、y(n-p)构成，所述y(n-1)、…、y(n-p)为电力系统信号的采样序列，由这些输入向量与p个参数a_k进行加权运算，建立AR模型方程，其方程如下：

$y\left(n\right)=\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{a}_{k}y(n-k)+e\left(n\right)---\left(1\right)$

其中e(n)为白噪声序列，其方差为σ²；p为AR模型的阶次；a_k为AR模型参数，k＝1、2、…、p；y(n)为电力系统信号的输出；

由(1)式可以看出，只需要求得AR模型的系数a_k，便可以建立起电力系统信号的状态空间模型，由维纳滤波器的基本原理可得，AR模型参数a_k的最优估计等价于输入向量为y(n-1)、…、y(n-p)，y(n)的维纳滤波器在时刻n的最优解

步骤B2、求得最优估计误差：

令Y(n)为n时刻滤波器的输入向量，

则Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，令

为AR参数最优解，则

将Y(n)和

代入(1)式可得到方程如下：

$y\left(n\right)=Y\left(n\right){\hat{A}}_0\left(n\right)+{\mathrm{\ϵ}}_0\left(n\right)---\left(2\right)$

其中ε₀(n)为输出y(n)的最优估计误差；

步骤C、建立状态方程和观测方程：

步骤C1、根据AR模型的参数向量建立状态方程，将AR参数估计的过程看作非平稳过程，引入过程噪声v(n)，

状态方程如下：

A(n+1)＝A(n)+v(n)     (3)

其中，A(n)＝[a₁(n)，a₂(n)，…，a_p(n)]^T，状态方程的转移矩阵为单位阵，v(n)是零均值平稳随机过程，v(n)的相关矩阵为Q(n)＝E[v(n)v^T(n)]＝qI，其中q为v(n)的方差，I为单位阵；

步骤C2、由步骤B2中的公式(2)建立观测方程，观测方程如下：

y(n)＝Y(n)A(n)+ε(n)  (4)

其中，Y(n)为观测矩阵，Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，观测噪声ε(n)均值为零，方差为

观测方程中的观测噪声{ε(n)}与状态方程中的过程噪声{v(n)}相互独立；

步骤D：利用最终预测误差准则确定AR模型阶数：

步骤D1：建立最终预测误差准则方程(5)和(6)，方程如下：

$S_p\left(N\right)=\underset{i=p+1}{\overset{N}{\mathrm{\Σ}}}(y\left(i\right)-{\hat{A}}_0\left(i\right)Y^T\left(i\right)){(y\left(i\right)-{\hat{A}}_0\left(i\right)Y^T\left(i\right))}^T---\left(5\right)$

其中S_p(N)为各个估计值的残差平方和，p为AR模型的阶次，y(i)为采样值，Y(i)＝[y(i-1)，y(i-2)，…y(i-p)]，为AR参数最优解；

$\mathrm{PPE}\left(p\right)=\frac{N+1+p}{N-1-p}\left(\frac{S_p\left(N\right)}{N}\right)---\left(6\right)$

其中N为信号数据的个数，p为AR模型的阶次；

将步骤一中采集的数据Y(n)＝[y(1)，y(2)，…y(N)]依次代入公式(5)中，经过运算得到估计值的残差平方和S_p(N)；

步骤D2：将步骤D1得到的估计值的残差平方和S_p(N)代入到公式(6)中，

依次计算PPE(1)、PPE(2)、……PPE(N/3)的值；

选取PPE(1)、PPE(2)、……PPE(N/3)中的最小值，所述最小值min{PPE(p)}为PPE的最优解从而得到AR模型的最佳阶数

步骤E：KF训练AR模型参数：

KF训练是利用自适应卡尔曼滤波方法对AR模型参数在线估计：根据AR模型的参数向量作为卡尔曼滤波的状态向量，将步骤A中的电力系统的信号数据y(1)，y(2)，…y(N)代入自适应卡尔曼滤波算法进行迭代运算，便可以得出AR模型参数a_k的无偏估计值；

步骤E1：设定{Y(i)，i＝p+1，p+2，…，N}为电力系统信号的输入向量，对应的电力系统信号的输入向量Y(p+1)，Y(p+2)，……，Y(N)的输出为y(p+1)，y(p+2)，……，y(N)，公式如下：

$\left\{\begin{array}{c}Y(p+1)\\ Y(p+1)\\ .\\ .\\ .\\ Y\left(N\right)\end{array}\right\}=\left\{\begin{array}{cccc}y\left(p\right)& y(p-1)& ...& y\left(1\right)\\ y(p+1)& y\left(p\right)& ...& y\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ y(N-1)& y(N-2)& ...& y(N-p)\end{array}\right\}---\left(7\right)$

步骤E2：令AR参数最优解的初始值协方差矩阵的初始值P(1，0)＝P(0)＝cI，其中C为很小的正数，取c＝2，I为单位阵，当n＝1，2，…N时，按照以下步骤进行迭代运算：

假设n＝1时，将P(0)、Y(N)和步骤C2中观测噪声ε(n)的方差

代入卡尔曼增益矩阵K_n公式进行计算：

$K_n=P(n-1)Y^T\left(n\right){(Y\left(n\right)P(n-1)Y^T\left(n\right)+{\mathrm{\σ}}_{\mathrm{\ϵ}}^2)}^{-1}---\left(8\right)$

从而计算出卡尔曼增益矩阵K_n；

步骤E3：将第n次的电力系统信号的输出y(n)与第n次电力系统信号的估计值

做差，得出卡尔曼新息β(n)：

$\mathrm{\β}\left(n\right)=y\left(n\right)-Y^T\left(n\right)\hat{A}(n-1)---\left(9\right)$

步骤E4：将新息和卡尔曼增益K_n的加权结果与模型参数相叠加，得到经过修正以后的模型参数

$\hat{A}\left(n\right)=\hat{A}(n-1)+K_n\mathrm{\β}\left(n\right)---\left(10\right)$

步骤E5：对预测协方差矩阵估计P(n，n-1)进行修正，得到协方差矩阵估计P(n)：

P(n)＝P(n，n-1)-K_nY^T(n)P(n，n-1)    (11)

步骤E6：将协方差矩阵估计P(n)与步骤C2中状态方程的相关矩阵Q(n)相叠加，得到预测协方差矩阵估计P(n，n+1)：

P(n，n+1)＝P(n)+Q(n)＝P(n)+qI       (12)

步骤E7：由AR模型进行功率谱估计：

将n＝1、n＝2、n＝3……n＝N代入步骤E2～步骤E6的自适应卡尔曼滤波算法进行多次迭代运算，当迭代次数趋于无限大时，由卡尔曼滤波算法得到的AR模型参数将趋于最优解，所求的AR模型参数

将会收敛于其最佳估计值

步骤F：将步骤E7求得的最佳估计值

代入到AR模型功率谱密度公式(13)中，公式(13)如下：

$p_y\left(f\right)=\frac{{{\mathrm{\σ}}_{\mathrm{\ϵ}}}^2}{{|1+\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{\hat{a}}_k{e}^{-2\mathrm{\πjk}\left(\frac{f_s}{f}\right)}|}^2}---\left(13\right)$

其中

为白噪声序列e(n)的方差；ω为圆频率；通过公式(13)计算出所需要的信号的功率谱p_y(ω)，得到各次谐波和间谐波的频率信息。

本发明具有以下有益效果：

本发明的电力系统间谐波估计方法与现有快速傅里叶变换方法相比，避免了FFT算法的栅栏效应，并且利用较短的数据即可得到较好的谱估计，大大减少了谱估计运算的数据量，可广泛应用于电力系统间谐波分析中；

本发明的电力系统间谐波估计方法采用迭代运算，适合在微处理器上应用，本发明的电力系统间谐波估计方法具有分辨率高、数值稳定性好和计算效率高优点；

本发明的电力系统间谐波估计方法采用AR模型作为描述电网信号的工具，与现有的Burg算法相比，不会出现谱峰偏移和频谱分裂现象，谱分析性能强、抗干扰能力强。

附图说明

图1是本发明电力系统间谐波估计方法的流程图，图2是采用快速傅里叶变换方法短采样时间情况下谱估计的检测结果走势图，图3是采用Burg算法短采样时间情况下谱估计的检测结果走势图，图4是采用本发明电力系统间谐波估计方法的短采样时间情况下谱估计的检测结果走势图。

具体实施方式

具体实施方式一：本实施方式的基于AR模型和卡尔曼滤波的电力系统间谐波估计方法是由以下步骤实现的：

步骤A、采集电力系统的信号数据，信号数据如下：

Y(n)＝[y(1)，y(2)，…y(N)]

其中N为信号数据的个数；

步骤B、由采集电力系统信号的数据建立AR模型，然后利用AR模型对电力系统信号进行分析；

步骤B1、采集电力系统信号的输入向量由y(n-1)、…、y(n-p)构成，所述y(n-1)、…、y(n-p)为电力系统信号的采样序列，由这些输入向量与p个参数a_k进行加权运算，建立AR模型方程，其方程如下：

$y\left(n\right)=\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{a}_{k}y(n-k)+e\left(n\right)---\left(1\right)$

其中e(n)为白噪声序列，其方差为σ²；p为AR模型的阶次；a_k为AR模型参数，k＝1、2、…、p；y(n)为电力系统信号的输出；

由(1)式可以看出，只需要求得AR模型的系数a_k，便可以建立起电力系统信号的状态空间模型，由维纳滤波器的基本原理可得，AR模型参数a_k的最优估计等价于输入向量为y(n-1)、…、y(n-p)，y(n)的维纳滤波器在时刻n的最优解

步骤B2、求得最优估计误差：

令Y(n)为n时刻滤波器的输入向量，

则Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，令

为AR参数最优解，则

将Y(n)和

代入(1)式可得到方程如下：

$y\left(n\right)=Y\left(n\right){\hat{A}}_0\left(n\right)+{\mathrm{\ϵ}}_0\left(n\right)---\left(2\right)$

其中ε₀(n)为输出y(n)的最优估计误差；

步骤C、建立状态方程和观测方程：

步骤C1、根据AR模型的参数向量建立状态方程，将AR参数估计的过程看作非平稳过程，引入过程噪声v(n)，

状态方程如下：

A(n+1)＝A(n)+v(n)        (3)

其中，A(n)＝[a₁(n)，a₂(n)，…，a_p(n)]^T，状态方程的转移矩阵为单位阵，v(n)是零均值平稳随机过程，v(n)的相关矩阵为Q(n)＝E[v(n)v^T(n)]＝qI，其中q为v(n)的方差，I为单位阵；

步骤C2、由步骤B2中的公式(2)建立观测方程，观测方程如下：

y(n)＝Y(n)A(n)+ε(n)     (4)

其中，Y(n)为观测矩阵，Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，观测噪声ε(n)均值为零，方差为

观测方程中的观测噪声{ε(n)}与状态方程中的过程噪声{v(n)}相互独立；

步骤D：利用最终预测误差准则确定AR模型阶数：

步骤D1：建立最终预测误差准则方程(5)和(6)，方程如下：

$S_p\left(N\right)=\underset{i=p+1}{\overset{N}{\mathrm{\Σ}}}(y\left(i\right)-{\hat{A}}_0\left(i\right)-Y^T\left(i\right)){(y\left(i\right)-{\hat{A}}_0\left(i\right)Y^T\left(i\right))}^T---\left(5\right)$

其中S_p(N)为各个估计值的残差平方和，p为AR模型的阶次，y(i)为采样值，Y(i)＝[y(i-1)，y(i-2)，…y(i-p)]，

为AR参数最优解；

$\mathrm{PPE}\left(p\right)=\frac{N+1+p}{N-1-p}\left(\frac{S_p\left(N\right)}{N}\right)---\left(6\right)$

其中N为信号数据的个数，p为AR模型的阶次；

将步骤一中采集的数据Y(n)＝[y(1)，y(2)，…y(N)]依次代入公式(5)中，经过运算得到估计值的残差平方和S_p(N)；

步骤D2：将步骤D1得到的估计值的残差平方和S_p(N)代入到公式(6)中，

依次计算PPE(1)、PPE(2)、……PPE(N/3)的值；

选取PPE(1)、PPE(2)、……PPE(N/3)中的最小值，所述最小值min{PPE(p)}为PPE的最优解

从而得到AR模型的最佳阶数

步骤E：KF训练AR模型参数：

KF训练是利用自适应卡尔曼滤波方法对AR模型参数在线估计：根据AR模型的参数向量作为卡尔曼滤波的状态向量，将步骤A中的电力系统的信号数据y(1)，y(2)，…y(N)代入自适应卡尔曼滤波算法进行迭代运算，便可以得出AR模型参数a_k的无偏估计值；

步骤E1：设定{Y(i)，i＝p+1，p+2，…，N}为电力系统信号的输入向量，对应的电力系统信号的输入向量Y(p+1)，Y(p+2)，……，Y(N)的输出为y(p+1)，y(p+2)，……，y(N)，公式如下：

$\left\{\begin{array}{c}Y(p+1)\\ Y(p+1)\\ .\\ .\\ .\\ Y\left(N\right)\end{array}\right\}=\left\{\begin{array}{cccc}y\left(p\right)& y(p-1)& ...& y\left(1\right)\\ y(p+1)& y\left(p\right)& ...& y\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ y(N-1)& y(N-2)& ...& y(N-p)\end{array}\right\}---\left(7\right)$

步骤E2：令AR参数最优解的初始值

协方差矩阵的初始值P(1，0)＝P(0)＝cI，其中C为很小的正数，取c＝2，I为单位阵，当n＝1，2，…N时，按照以下步骤进行迭代运算：

假设n＝1时，将P(0)、Y(N)和步骤C2中观测噪声ε(n)的方差

代入卡尔曼增益矩阵K_n公式进行计算：

$K_n=P(n-1)Y^T\left(n\right){(Y\left(n\right)P(n-1)Y^T\left(n\right)+{\mathrm{\σ}}_{\mathrm{\ϵ}}^2)}^{-1}---\left(8\right)$

从而计算出卡尔曼增益矩阵K_n；

步骤E3：将第n次的电力系统信号的输出y(n)与第n次电力系统信号的估计值做差，得出卡尔曼新息β(n)：

$\mathrm{\β}\left(n\right)=y\left(n\right)-Y^T\left(n\right)\hat{A}(n-1)---\left(9\right)$

步骤E4：将新息和卡尔曼增益K_n的加权结果与模型参数

相叠加，得到经过修正以后的模型参数

$\hat{A}\left(n\right)=\hat{A}(n-1)+K_n\mathrm{\β}\left(n\right)---\left(10\right)$

步骤E5：对预测协方差矩阵估计P(n，n-1)进行修正，得到协方差矩阵估计P(n)：

P(n)＝P(n，n-1)-K_nY^T(n)P(n，n-1)    (11)

步骤E6：将协方差矩阵估计P(n)与步骤C2中状态方程的相关矩阵Q(n)相叠加，得到预测协方差矩阵估计P(n，n+1)：

P(n，n+1)＝P(n)+Q(n)＝P(n)+qI    (12)

步骤E7：将n＝1、n＝2、n＝3……n＝N代入步骤E2～步骤E6的自适应卡尔曼滤波算法进行多次迭代运算，当迭代次数趋于无限大时，由卡尔曼滤波算法得到的AR模型参数将趋于最优解，所求的AR模型参数

将会收敛于其最佳估计值

步骤F：由AR模型进行功率谱估计：

将步骤E7求得的最佳估计值代入到AR模型功率谱密度公式(13)中，公式(13)如下：

$p_y\left(f\right)=\frac{{{\mathrm{\σ}}_{\mathrm{\ϵ}}}^2}{{|1+\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{\hat{a}}_k{e}^{-2\mathrm{\πjk}\left(\frac{f_s}{f}\right)}|}^2}---\left(13\right)$

其中为白噪声序列e(n)的方差；ω为圆频率；通过公式(13)计算出所需要的信号的功率谱p_y(ω)，得到各次谐波和间谐波的频率信息。

所述AR模型为自回归模型。

本发明方法的实验验证：

应用本发明对电力系统间谐波的功率谱进行估计，结果证明本发明方法的确可行有效，具有分辨率高、数值稳定性好和抗干扰能力低的优点；这里，选取快速傅里叶变换方法和Burg算法与本发明的电力系统间谐波的功率谱方法进行比较，说明本发明方法完全可达到发明目的，如下：

IEC标准(国际电工委员会标准)推荐的采样窗长度为200ms，但是对于实际电网系统的间谐波检测采样时间就过于长，其检测过程的实时性要求较高，间谐波频率出现是随机的，其变化也较快且持续时间短暂，长采样窗会给间谐波的频率估计带来较大的误差。在实际应用中，一般只采样2～3个基频周期，即采样窗长度为40～60ms(50Hz电力系统)。

设信号除50Hz的基波外还加入了五个谐波和间谐波分量，各个分量的频率、幅度和初始相位参数见表1。

表1仿真实验加入的谐波与间谐波信息

为了模拟实际电网环境，加入σ＝0.1的零均值白噪声信号，采样频率为1600Hz，采样时间减少到40ms。则快速傅里叶变换方法、Burg算法和本发明的方法(KF-AR)短采样时间情况下的检测结果如图2～4所示(其中Burg算法和本发明的方法的AR模型阶数均为25)：

采样时间窗长度减小，直接影响快速傅里叶变换方法的频率分辨率，严重降低了其谱估计性能，除了能勉强检测基频和150Hz、250Hz谐波的谱峰，快速傅里叶变换方法几乎完全漏掉了其他三个间谐波频率分量；Burg算法虽然还能较为清晰地显示出基频、150Hz、250Hz和325Hz的谱峰分量，但是分别与基频和325Hz相邻较近的85Hz和350Hz分量却被完全漏掉；而本发明的方法(KF-AR)依然能够提取出全部的6个频率分量，且谱峰较为尖锐，易于检测；具体结果如表2所示：

表2是快速傅里叶变换方法、Burg算法和本发明的方法短采样时间下三种谱估计结果：

注：表中的X表示此频率分量峰值无法检测。

从表中可以看出快速傅里叶变换方法和Burg算法仅能提取部分频谱分量信息，且其检测精度较低，误差分别高达1.67Hz和1.26Hz；而本发明的电力系统间谐波估计方法可以提取出全部的频谱分量，而且精度较高，平均误差仍然保持在0.34Hz的较低水平；本发明的电力系统间谐波估计方法的谱分辨率所受影响远小于其对FFT算法和Burg算法，更适合实际间谐波检测，在满足实时性的同时，还保持着较高的检测精度。

Claims (1)

1.一种基于AR模型和卡尔曼滤波的电力系统间谐波估计方法，其特征在于所述间谐波估计方法由以下步骤实现的：

步骤A、采集电力系统的信号数据，信号数据如下：

Y(n)＝[y(1)，y(2)，…y(N)]

其中N为信号数据的个数；

步骤B、由采集电力系统信号的数据建立AR模型，然后利用AR模型对电力系统信号进行分析；

步骤B1、采集电力系统信号的输入向量由y(n-1)、…、y(n-p)构成，所述y(n-1)、…、y(n-p)为电力系统信号的采样序列，由这些输入向量与p个参数a_k进行加权运算，建立AR模型方程，其方程如下：

$y\left(n\right)=\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{a}_{k}y(n-k)+e\left(n\right)---\left(1\right)$

其中e(n)为白噪声序列，其方差为σ²；p为AR模型的阶次；a_k为AR模型参数，k＝1、2、…、p；y(n)为电力系统信号的输出；

由(1)式可以看出，只需要求得AR模型的系数a_k，便可以建立起电力系统信号的状态空间模型，由维纳滤波器的基本原理可得，AR模型参数a_k的最优估计等价于输入向量为y(n-1)、…、y(n-p)，y(n)的维纳滤波器在时刻n的最优解

步骤B2、求得最优估计误差：

令Y(n)为n时刻滤波器的输入向量，

则Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，令

为AR参数最优解，则

将Y(n)和

代入(1)式可得到方程如下：

$y\left(n\right)=Y\left(n\right){\hat{A}}_0\left(n\right)+{\mathrm{\ϵ}}_0\left(n\right)---\left(2\right)$

其中ε₀(n)为输出y(n)的最优估计误差；

步骤C、建立状态方程和观测方程：

步骤C1、根据AR模型的参数向量建立状态方程，将AR参数估计的过程看作非平稳过程，引入过程噪声v(n)，

状态方程如下：

A(n+1)＝A(n)+v(n)        (3)

其中，A(n)＝[a₁(n)，a₂(n)，…，a_p(n)]^T，状态方程的转移矩阵为单位阵，v(n)是零均值平稳随机过程，v(n)的相关矩阵为Q(n)＝E[v(n)v^T(n)]＝qI，其中q为v(n)的方差，I为单位阵；

步骤C2、由步骤B2中的公式(2)建立观测方程，观测方程如下：

y(n)＝Y(n)A(n)+ε(n)      (4)

其中，Y(n)为观测矩阵，Y(n)＝[y(n-1)，y(n-2)，…y(n-p)]，观测噪声ε(n)均值为零，方差为

观测方程中的观测噪声{ε(n)}与状态方程中的过程噪声{v(n)}相互独立；

步骤D：利用最终预测误差准则确定AR模型阶数：

步骤D1：建立最终预测误差准则方程(5)和(6)，方程如下：

$S_p\left(N\right)=\underset{i=p+1}{\overset{N}{\mathrm{\Σ}}}(y\left(i\right)-{\hat{A}}_0\left(i\right)Y^T\left(i\right)){(y\left(i\right)-{\hat{A}}_0\left(i\right)Y^T\left(i\right))}^T---\left(5\right)$

其中S_p(N)为各个估计值的残差平方和，p为AR模型的阶次，y(i)为采样值，Y(i)＝[y(i-1)，y(i-2)，…y(i-p)]，

为AR参数最优解；

$\mathrm{PPE}\left(p\right)=\frac{N+1+p}{N-1-p}\left(\frac{S_p\left(N\right)}{N}\right)---\left(6\right)$

其中N为信号数据的个数，p为AR模型的阶次；

将步骤一中采集的数据Y(n)＝[y(1)，y(2)，…y(N)]依次代入公式(5)中，经过运算得到估计值的残差平方和S_p(N)；

步骤D2：将步骤D1得到的估计值的残差平方和S_p(N)代入到公式(6)中，

依次计算PPE(1)、PPE(2)、……PPE(N/3)的值；

选取PPE(1)、PPE(2)、……PPE(N/3)中的最小值，所述最小值min{PPE(p)}为PPE的最优解从而得到AR模型的最佳阶数

步骤E：KF训练AR模型参数：

KF训练是利用自适应卡尔曼滤波方法对AR模型参数在线估计：根据AR模型的参数向量作为卡尔曼滤波的状态向量，将步骤A中的电力系统的信号数据y(1)，y(2)，…y(N)代入自适应卡尔曼滤波算法进行迭代运算，便可以得出AR模型参数a_k的无偏估计值；

步骤E1：设定{Y(i)，i＝p+1，p+2，…，N}为电力系统信号的输入向量，对应的电力系统信号的输入向量Y(p+1)，Y(p+2)，……，Y(N)的输出为y(p+1)，y(p+2)，……，y(N)，公式如下：

$\left\{\begin{array}{c}Y(p+1)\\ Y(p+1)\\ .\\ .\\ .\\ Y\left(N\right)\end{array}\right\}=\left\{\begin{array}{cccc}y\left(p\right)& y(p-1)& ...& y\left(1\right)\\ y(p+1)& y\left(p\right)& ...& y\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ y(N-1)& y(N-2)& ...& y(N-p)\end{array}\right\}---\left(7\right)$

步骤E2：令AR参数最优解的初始值

协方差矩阵的初始值P(1，0)＝P(0)＝cI，其中C为很小的正数，取c＝2，I为单位阵，当n＝1，2，…N时，按照以下步骤进行迭代运算：

假设n＝1时，将P(0)、Y(N)和步骤C2中观测噪声ε(n)的方差

代入卡尔曼增益矩阵K_n公式进行计算：

$K_n=P(n-1)Y^T\left(n\right){(Y\left(n\right)P(n-1)Y^T\left(n\right)+{\mathrm{\σ}}_{\mathrm{\ϵ}}^2)}^{-1}---\left(8\right)$

从而计算出卡尔曼增益矩阵K_n；

步骤E3：将第n次的电力系统信号的输出y(n)与第n次电力系统信号的估计值

做差，得出卡尔曼新息β(n)：

$\mathrm{\β}\left(n\right)=y\left(n\right)-Y^T\left(n\right)\hat{A}(n-1)---\left(9\right)$

步骤E4：将新息和卡尔曼增益K_n的加权结果与模型参数

相叠加，得到经过修正以后的模型参数

$\hat{A}\left(n\right)=\hat{A}(n-1)+K_n\mathrm{\β}\left(n\right)---\left(10\right)$

步骤E5：对预测协方差矩阵估计P(n，n-1)进行修正，得到协方差矩阵估计P(n)：

P(n)＝P(n，n-1)-K_nY^T(n)P(n，n-1)    (11)

步骤E6：将协方差矩阵估计P(n)与步骤C2中状态方程的相关矩阵Q(n)相叠加，得到预测协方差矩阵估计P(n，n+1)：

P(n，n+1)＝P(n)+Q(n)＝P(n)+qI       (12)

步骤E7：将n＝1、n＝2、n＝3……n＝N代入步骤E2～步骤E6的自适应卡尔曼滤波算法进行多次迭代运算，当迭代次数趋于无限大时，由卡尔曼滤波算法得到的AR模型参数将趋于最优解，所求的AR模型参数

将会收敛于其最佳估计值

步骤F：由AR模型进行功率谱估计：

将步骤E7求得的最佳估计值

代入到AR模型功率谱密度公式(13)中，公式(13)如下：

$p_y\left(f\right)=\frac{{{\mathrm{\σ}}_{\mathrm{\ϵ}}}^2}{{|1+\underset{k=1}{\overset{p}{\mathrm{\Σ}}}{\hat{a}}_k{e}^{-2\mathrm{\πjk}\left(\frac{f_s}{f}\right)}|}^2}---\left(13\right)$

其中

为白噪声序列e(n)的方差；ω为圆频率；通过公式(13)计算出所需要的信号的功率谱p_y(ω)，得到各次谐波和间谐波的频率信息。

Images