CN105956565A - 一种考虑量测信号丢失的动态振荡信号参数辨识方法 - Google Patents

Abstract

本发明公开了一种考虑量测信号丢失的动态振荡信号参数辨识方法。该方法首先采用二项分布建立了量测信号丢失的数学模型；然后，在此基础上设计出了适用于量测信号丢失情形下动态振荡信号参数辨识的方法，即改进的扩展卡尔曼滤波算法。该算法在设计时不仅考虑了量测信号丢失的情况，而且引入了系统噪声和量测噪声协方差矩阵的动态计算方法，可以有效的避免因噪声协方差矩阵设置不当而引起的算法失效，提高了算法的效率。该算法因考虑了实际工程背景，且简单方便，具有较高的工程应用价值。

Description

一种考虑量测信号丢失的动态振荡信号参数辨识方法

技术领域

本发明涉及一种考虑量测信号丢失的动态振荡信号参数辨识方法，属于信号分析与参数辨识技术领域。

背景技术

近年来，随着现代电网规模不断扩大，电网互联程度的日益提高，系统受到大、小扰动后产生的动态振荡已经成为制约电网安全稳定运行的最主要的因素之一。由于这些动态振荡信号可以提供关于电力系统运行模式的重要信息，所以发现并准确掌握这些振荡信号特征对于电力系统安全稳定运行具有重要意义。

鉴于动态振荡信号辨识的重要性，研究人员提出许多辨识方法，如最大似然法、普龙尼法，卡尔曼滤波方法等。但是，这些方法在设计时未考虑量测信号发生数据丢包情况。因此，在量测信号发生丢包时，这些方法无法有效的对动态振荡信号的参数进行有效的辨识。所以，研究量测信号丢包的动态振荡信号参数辨识方法具有重要现实意义和工程应用价值。

发明内容

发明目的：针对现有技术中存在的问题，为了有效地解决量测信号丢包时动态振荡信号的参数辨识，克服已有辨识算法的缺点，本发明设计了一种考虑量测信号丢失的动态振荡信号参数辨识方法，有效的实现了量测信号丢包情况下的动态振荡信号参数辨识。

技术方案：一种考虑量测信号丢失的动态振荡信号参数辨识方法，该方法在计算机中是依次按照如下步骤实现的：

(1)、获取量测信号丢失的离散状态空间表达式，模型公式为：

$\left\{\begin{array}{c}x\left(k+1\right)=f\left(x\left(k\right)\right)+w\left(k\right)\\ y\left(k\right)=\mathrm{\Ξ}\left(k\right)h\left(x\left(k\right)\right)+v\left(k\right)\end{array}\right.$

式中，x(k)表示k时刻的状态向量，y(k)表示k时刻的输出量测向量，f(·)和h(·)是对应于具体问题中的非线性函数，w(k)和v(k)分别是系统噪声和量测噪声，Ξ(k)是符合二项分布的白噪声随机序列，即，Ξ(k)＝diag{γ₁(k),γ₂(k),…,γ_m(k)}，γ_i(k)(i＝1...m)是不相关的随机变量，且与w(k)和v(k)，x(0)不相关。

(2)、γ_i(k)的取值为0或者1，且满足如下计算公式：

prob(γ_i(k)＝0)＝1-μ_i(k)

prob(γ_i(k)＝1)＝μ_i(k)

式中，prob(γ_i(k)＝0)＝1-μ_i(k)表示的是第i个量测量的丢失率。

(3)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的初始协方差矩阵Q₀和R₀，整体算法迭代次数最大值S；

(4)、获取量测数据丢包情况下的量测数据y(k)；

(5)、由已知条件，计算k时刻的系统噪声和量测噪声所满足的动态协方差矩阵Q(k)和R(k)，计算公式为：

$s\left(k\right)=y\left(k\right)-\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)h\left(\tilde{x}\right(k\left)\right)$

$C_{vk}=\frac{1}{N}\underset{i=k-N+1}{\overset{k}{\Σ}}s\left(k\right)s{\left(k\right)}^T$

Q(k)＝G(k)C_vkK(k)^T

$R\left(k\right)=C_{vk}+\left(\stackrel{\‾}{\mathrm{\Ξ}}\right(k)\·H(k\left)\right)\·\hat{P}\left(k\right)\·{\left(\stackrel{\‾}{\mathrm{\Ξ}}\right(k)\·H(k\left)\right)}^T$

式中k代表第k时刻迭代，N是动态估计窗口值，其是一个常数，G(k)是k时刻的卡尔曼滤波增益，是k时刻的状态估计协方差，表示非线性函数h(·)在处的雅克比矩阵，其中h(·)对应具体问题输出方程中的非线性函数，上标T表示矩阵的装置。系数的计算公式如下：

$\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)=E\left\{\mathrm{\Ξ}\left(k\right)\right\}=diag\left\{\begin{array}{cccc}{\mathrm{\μ}}_1\left(k\right)& {\mathrm{\μ}}_2\left(k\right)& \mathrm{...}& {\mathrm{\μ}}_m\left(k\right)\end{array}\right\}$

(6)、由已经得到的k-1时刻的状态估计值和状态估计误差协方差，利用本发明算法的预测步，得到k时刻的状态预测值和状态预测误差协方差，计算公式为：

$\tilde{x}\left(k\right)=f\left(\hat{x}\right(k-1\left)\right)$

$\tilde{P}\left(k\right)=F_{k-1}\hat{P}(k-1)F_{k-1}^T+Q(k-1)$

式中，表示k时刻的状态预测值，f(·)表示非线性函数，表示k-1时刻的状态估计向量。表示k时刻的状态预测误差协方差，表示非线性函数f(x)在处的雅克比矩阵，表示k-1时刻的状态估计误差协方差，上标T表示转置，Q(k-1)是系统噪声k-1时刻所满足的动态协方差矩阵。

(7)、由已知条件计算k时刻的中间迭代变量M(k)，其计算公式为：

式中ο符号表示Hadamard乘积，其定义为[AοB]_ij＝[A_ij×B_ij]，式中系数的取值计算规则为：

$\widetilde{\mathrm{\Ξ}}\left(k\right)=diag\{{\widetilde{\mathrm{\μ}}}_1\left(k\right),{\widetilde{\mathrm{\μ}}}_2\left(k\right),...,{\widetilde{\mathrm{\μ}}}_m\left(k\right)\}$

${\widetilde{\mathrm{\μ}}}_i\left(k\right)={\mathrm{\μ}}_i\left(k\right)(1-{\mathrm{\μ}}_i(k\left)\right),(i=1,2,...m)$

(8)、利用k时刻的中间迭代值M(k)和k时刻的预测误差协方差计算k时刻的最优滤波增益，计算步骤为：

$G\left(k\right)=\tilde{P}\left(k\right)H^T\left(k\right)\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)M^{-1}\left(K\right)$

式中，G(k)表示k时刻的最优滤波增益，表示k时刻的状态预测误差协方差，上标T表示转置，表示非线性函数h(x)在处的雅克比矩阵。

(9)、利用k时刻的最优滤波增益和k时刻的预测误差协方差计算k时刻的估计误差协方差，计算步骤为：

$\hat{P}\left(k\right)=(I-G(k\left)\stackrel{\‾}{\mathrm{\Ξ}}\right(k\left)H\right(k\left)\right)\tilde{P}\left(k\right)$

(10)、结合已得到的k时刻的最优滤波增益和k时刻的状态预测值计算k时刻的状态估计值，计算步骤为：

$\hat{x}\left(k\right)=\tilde{x}\left(k\right)+G\left(k\right)\[y\left(k\right)-\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)h\left(\tilde{x}\right(k\left)\right)\]$

(11)、按照上述步骤，进行多次迭代辨识，若k≤S，则迭代继续，若k＞S，则迭代结束，输出辨识结果。

附图说明

图1为本发明实施例的方法流程图；

图2实施例的动态振荡信号；

图3为实施例采用本发明方法以及EKF算法的信号频率辨识结果对比；

图4为实施例采用本发明方法以及EKF算法的信号阻尼因子辨识结果对比；

图5为实施例采用本发明方法信号频率和阻尼因子的辨识误差。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，动态振荡信号参数辨识方法，其包含如下步骤：

(1)、获取量测信号丢失的离散状态空间表达式；

(2)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的初始协方差矩阵Q₀和R₀，整体算法迭代次数最大值S；

(3)、获取量测数据丢包情况下的量测数据y(k)；

(4)、由已知条件，计算k时刻的系统噪声和量测噪声所满足的协方差矩阵Q(k)和R(k)；

(5)、由已经得到的k-1时刻的状态估计值和状态估计误差协方差，利用本发明算法的预测步，得到k时刻的状态预测值和状态预测误差协方差；

(6)、由已知条件计算k时刻的中间迭代变量M(k)；

(7)、利用k时刻的中间迭代值M(k)和k时刻的预测误差协方差计算k时刻的最优滤波增益；

(8)、利用k时刻的最优滤波增益和k时刻的预测误差协方差计算k时刻的估计误差协方差；

(9)、结合已得到的k时刻的最优滤波增益和k时刻的状态预测值计算k时刻的状态估计值；

(10)、按照上述步骤，进行多次迭代辨识，若k≤S，则迭代继续，若k＞S，则迭代结束，输出辨识结果。

一般情况下动态振荡信号可以表示为多个指数衰减的正弦信号的之和，可以描述为如下形式：

$y\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{A}_i{e}^{-{\mathrm{\δ}}_{i}t}\cos \left(w_{i}t+{\mathrm{\φ}}_i\right)+n\left(t\right)$

式中，A_i,δ_i,w_i,φ_i是实数的未知参数，n(t)是一个零均值白噪声。其中，δ_i称为动态信号的阻尼因子，w_i是动态信号的频率，其中w_i，δ_i为待估参数。经过推理可以得到动态信号的状态变量分量中包含待估参数的离散状态空间模型。考虑由N个指数衰减的正弦信号总和组成的动态信号，其4N个状态变量形式可以表达如下：

$x_{4i-3,k}=A_i{e}^{-{\mathrm{\δ}}_i\frac{k}{f_s}}cos\left(w_i\frac{k}{f_s}\right)$

$x_{4i-2,k}=A_i{e}^{-{\mathrm{\δ}}_i\frac{k}{f_s}}sin\left(w_i\frac{k}{f_s}\right)$

x_4i-1,k＝w_i

x_4i,k＝δ_i

式中i代表这些变量和参数是属于动态信号的第i个衰减正弦信号。k代表时刻，f_s代表采样频率。根据推理可得到k+1时刻的状态分量：

$x_{4i-3,k+1}=e^{-\frac{x_{4i,k}}{f_s}}\[x_{4i-3,k}cos\left(\frac{x_{4i-1,k}}{f_s}\right)-x_{4i-2,k}sin\left(\frac{x_{4i-1,k}}{f_s}\right)\]+w_{4i-3,k}$

$x_{4i-2,k+1}=e^{-\frac{x_{4i,k}}{f_s}}\[x_{4i-3,k}sin\left(\frac{x_{4i-1,k}}{f_s}\right)+x_{4i-2,k}cos\left(\frac{x_{4i-1,k}}{f_s}\right)\]+w_{4i-2,k}$

x_4i-1,k+1＝x_4i-1,k+w_4i-1,k

x_4i,k+1＝x_4i,k+w_4i,k

则考虑量测数据丢包的输出方程为：

$y_k=\mathrm{\Ξ}\left(k\right)\left(\underset{i=1}{\overset{N}{\Σ}}{k}_{2i-1}{x}_{4i-3,k}+k_{2i}{x}_{4i-2,k}\right)+n_k$

式中，k_2i-1＝cos(φ_i)，k_2i＝-sin(φ_i)，n_k为均值为零的白噪声，Ξ(k)为对应维度的二项分布随机序列。

所以，考虑量测数据丢包的动态振荡信号的状态空间模型一般可以表示为：

$\left\{\begin{array}{c}{x}_{k+1}=f\left(x_k\right)+w\left(k\right)\\ y\left(k\right)=\mathrm{\Ξ}\left(k\right)h\left(x_k\right)+v\left(k\right)\end{array}\right.$

式中，f(·)和h(·)代表可以根据泰勒级数展开进行线性化的非线性函数，w(k)和v(k)是均值为零的高斯白噪声序列，分别满足动态协方差矩阵Q(k)和R(k)。具体而言：

$f\left(x_k\right)=\left[\begin{array}{c}{M}_1\\ M_2\\ .\\ .\\ .\\ M_i\\ .\\ .\\ .\\ M_N\end{array}\right],M_i=\left[\begin{array}{c}{x}_{4i-3,k}\\ x_{4i-2,k}\\ x_{4i-1,k}\\ x_{4i,k}\end{array}\right]$

而函数h(x_k)可以表示为如下形式：

H＝(k₁k₂00…,k_2i-1k_2i00…,k_2N-1k_2N00)

h(x_k)＝Hx_k

至此，考虑量测数据丢包的动态振荡信号状态空间表达式已建立。在此基础上，则可以运用本发明所介绍的方法，进行量测信号丢包情况下的动态振荡信号参数辨识，得到辨识结果。

下面介绍本发明的一个实施例：

考虑动态振荡信号为：

y(k)＝Ξ(k)e^-0.01k sin(0.4k)+n_k,0≤k≤300

式中k是信号采样时刻，n_k是高斯白噪声。如图2所示，该动态信号是由一个指数衰减的正弦信号组成，该动态振荡信号的频率为w＝0.4，阻尼因子为δ＝0.01。在运用本发明所提出的方法进行动态信号参数辨识时，所采用的相关初始参数值为：

R₀＝10^-4

${\hat{x}}_0={\left[\begin{array}{cccc}0.1& 0& 0.3& 0.004\end{array}\right]}^T$

在对本算例进行辨识时，假设量测信号的丢包率u＝0.1，在计算噪声所满足的动态误差协方差矩阵时，动态估计窗口值N＝2。

图1为实施例所用的算法流程图，图2为实施例的动态振荡信号，图3为实施例采用本发明方法以及EKF算法的信号频率w辨识结果对比，图4为实施例采用本发明方法以及EKF算法的信号阻尼因子δ辨识结果对比，图5为实施例采用本发明方法信号频率和阻尼因子的辨识误差。对比结果表明，本发明所提的方法有效的实现了量测信号丢失情形下的动态振荡信号参数辨识，且较EKF算法具有更好的收敛性和精度。

Claims (5)

1.一种考虑量测信号丢失的动态振荡信号参数辨识方法，其特征在于，包含如下步骤：

(1)、获取量测信号丢失的离散状态空间表达式；

(2)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的初始协方差矩阵Q₀和R₀，整体算法迭代次数最大值S；

(3)、获取量测数据丢包情况下的量测数据y(k)；

(4)、由已知条件，计算k时刻的系统噪声和量测噪声所满足的协方差矩阵Q(k)和R(k)；

(5)、由已经得到的k-1时刻的状态估计值和状态估计误差协方差，利用本发明算法的预测步，得到k时刻的状态预测值和状态预测误差协方差；

(6)、由已知条件计算k时刻的中间迭代变量M(k)；

(7)、利用k时刻的中间迭代值M(k)和k时刻的预测误差协方差计算k时刻的最优滤波增益；

(8)、利用k时刻的最优滤波增益和k时刻的预测误差协方差计算k时刻的估计误差协方差；

(9)、结合已得到的k时刻的最优滤波增益和k时刻的状态预测值计算k时刻的状态估计值；

(10)、按照上述步骤，进行多次迭代辨识，若k≤S，则迭代继续，若k＞S，则迭代结束，输出辨识结果。

2.如权利要求1所述的考虑量测信号丢失的动态振荡信号参数辨识方法，其特征在于，获取量测信号丢失的离散状态空间表达式，模型公式为：

$\left\{\begin{array}{c}x(k+1)=f\left(x\right(k\left)\right)+w\left(k\right)\\ y\left(k\right)=\mathrm{\Ξ}\left(k\right)h\left(x\right(k\left)\right)+v\left(k\right)\end{array}\right.$

式中，x(k)表示k时刻的状态向量，y(k)表示k时刻的输出量测向量，f(·)和h(·)是对应于具体问题中的非线性函数，w(k)和v(k)分别是系统噪声和量测噪声，Ξ(k)是符合二项分布的白噪声随机序列，即，Ξ(k)＝diag{γ₁(k),γ₂(k),…,γ_m(k)}，γ_i(k)(i＝1...m)是不相关的随机变量，且与w(k)和v(k)，x(0)不相关。

(2)、γ_i(k)的取值为0或者1，且满足如下计算公式：

prob(γ_i(k)＝0)＝1-μ_i(k)

prob(γ_i(k)＝1)＝μ_i(k)

式中，prob(γ_i(k)＝0)＝1-μ_i(k)表示的是第i个量测量的丢失率。

3.如权利要求2所述的考虑量测信号丢失的动态振荡信号参数辨识方法，其特征在于，由已知条件，计算k时刻的系统噪声和量测噪声所满足的动态协方差矩阵Q(k)和R(k)，计算公式为：

$s\left(k\right)=y\left(k\right)-\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)h\left(\tilde{x}\right(k\left)\right)$

$C_{vk}=\frac{1}{N}\underset{i=k-N+1}{\overset{k}{\Σ}}s\left(k\right)s{\left(k\right)}^T$

Q(k)＝G(k)C_vkG(k)^T

$R\left(k\right)=C_{vk}+\left(\stackrel{\‾}{\mathrm{\Ξ}}\right(k)\·H(k\left)\right)\·\hat{P}\left(k\right)\·{\left(\stackrel{\‾}{\mathrm{\Ξ}}\right(k)\·H(k\left)\right)}^T$

式中k代表第k时刻迭代，N是动态估计窗口值，其是一个常数，G(k)是k时刻的卡尔曼滤波增益，是k时刻的状态估计协方差，表示非线性函数h(·)在处的雅克比矩阵，其中h(·)对应具体问题输出方程中的非线性函数，上标T表示矩阵的装置。系数的计算公式如下：

$\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)=E\left\{\mathrm{\Ξ}\left(k\right)\right\}=diag\left\{\begin{array}{cccc}{\mathrm{\μ}}_1\left(k\right)& {\mathrm{\μ}}_2\left(k\right)& \mathrm{...}& {\mathrm{\μ}}_m\left(k\right)\end{array}\right\}.$

4.如权利要求3所述的考虑量测信号丢失的动态振荡信号参数辨识方法，其特征在于，由已经得到的k-1时刻的状态估计值和状态估计误差协方差，利用本发明算法的预测步，得到k时刻的状态预测值和状态预测误差协方差，计算公式为：

$\tilde{x}\left(k\right)=f\left(\hat{x}\right(k-1\left)\right)$

$\tilde{P}\left(k\right)=F_{k-1}\hat{P}(k-1)F_{k-1}^T+Q(k-1)$

式中，表示k时刻的状态预测值，f(·)表示非线性函数，表示k-1时刻的状态估计向量；表示k时刻的状态预测误差协方差，表示非线性函数f(x)在处的雅克比矩阵，表示k-1时刻的状态估计误差协方差，上标T表示转置，Q(k-1)是系统噪声k-1时刻所满足的动态协方差矩阵；

由已知条件计算k时刻的中间迭代变量M(k)，其计算公式为：

式中о符号表示Hadamard乘积，其定义为[AоB]_ij＝[A_ij×B_ij]，式中系数的取值计算规则为：

$\widetilde{\mathrm{\Ξ}}\left(k\right)=diag\{{\widetilde{\mathrm{\μ}}}_1\left(k\right),{\widetilde{\mathrm{\μ}}}_2\left(k\right),...,{\widetilde{\mathrm{\μ}}}_m\left(k\right)\}$

${\widetilde{\mathrm{\μ}}}_i\left(k\right)={\mathrm{\μ}}_i\left(k\right)(1-{\mathrm{\μ}}_i(k\left)\right),(i=1,2,...m).$

5.如权利要求4所述的考虑量测信号丢失的动态振荡信号参数辨识方法，其特征在于，利用k时刻的中间迭代值M(k)和k时刻的预测误差协方差计算k时刻的最优滤波增益，计算步骤为：

$G\left(k\right)=\tilde{P}\left(k\right)H^T\left(k\right)\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)M^{-1}\left(K\right)$

式中，G(k)表示k时刻的最优滤波增益，表示k时刻的状态预测误差协方差，上标T表示转置，表示非线性函数h(x)在处的雅克比矩阵；

利用k时刻的最优滤波增益和k时刻的预测误差协方差计算k时刻的估计误差协方差，计算步骤为：

$\hat{P}\left(k\right)=(I-G(k\left)\stackrel{\‾}{\mathrm{\Ξ}}\right(k\left)H\right(k\left)\right)\tilde{P}\left(k\right)$

结合已得到的k时刻的最优滤波增益和k时刻的状态预测值计算k时刻的状态估计值，计算步骤为：

$\hat{x}\left(k\right)=\tilde{x}\left(k\right)+G\left(k\right)\[y\left(k\right)-\stackrel{\‾}{\mathrm{\Ξ}}\left(k\right)h\left(\tilde{x}\right(k\left)\right)\].$