CN103632050B - 电力***噪声自适应抗差状态估计方法 - Google Patents

Abstract

本发明提出一种电力***噪声自适应抗差状态估计方法，包括以下步骤：S1，获取L个量测断面，其中，L为正整数；S2，对每个量测断面进行量测误差估计以获取误差矢量；以及S3，基于L个量测断面估计得到的误差矢量，运用统计学习方法估计得到广义高斯密度模型GGD的参数，以获取量测噪声的分布类型，并根据噪声分布类型选择对应的最优抗差状态估计模型。本发明的噪声自适应抗差状态估计方法，能够通过统计学习获取噪声的分布规律，并将其与抗差状态估计方法进行在线匹配，从而可实现对各种噪声类型的自适应，即在任何噪声分布类型下，均可以得到更接近于状态变量真值的最优估计结果。

Description

电力***噪声自适应抗差状态估计方法

技术领域

本发明属于电力***调度自动化领域，具体涉及一种电力***噪声自适应抗差状态估计方法。

背景技术

电力***状态估计是能量管理***的基础和核心，是电力***安全、可靠、优质和经济运行必不可少的保障环节。状态估计在国内外的应用已经有四十多年的历史，这期间涌现了各种各样的状态估计方法。

状态估计中的量测量一般包括三个部分，即量测量的真值、相对不良数据及量测噪声。最优状态估计不仅应能可靠辨识不良数据，而且应是相应量测噪声分布类型下的最优估计。传统状态估计方法往往假定噪声服从某种已知的分布，并选择相应的准则函数，选择不同的准则函数可得到不同的状态估计模型，进而得到不同的状态变量估计结果，而特定的准则函数仅在特定的噪声分布下是最优估计。当假定的噪声分布与实际噪声分布不一致时，所得到的状态估计结果往往并不是最优估计。以加权最小二乘法(Weighted LeastSquares，WLS)为例，传统统计学中的大数定律保证了当量测量的数目趋近于无穷大时，任何噪声分布都趋近于高斯分布，高斯分布下的最优估计WLS也在国内外得到了广泛应用。但是，这样做仍然存在两个缺点：（1）WLS不具有抗差性，为此在WLS估计之后往往要用最大正则化残差检验法或估计辨识法等来辨识不良数据，但这些方法在大多数情况下无法有效辨识强相关的不良数据；（2）状态估计可用的量测量的数目毕竟有限，此时大数定律不再成立，噪声的实际分布与高斯分布往往不符合（尤其是电流幅值量测），因而WLS估计结果的性能没有理论上的保证。与WLS一样，其他的传统状态估计方法也存在同样的问题，而由国内学者提出的最小信息损失法（Minimum Information Loss，MIL）对噪声具有较强的适应性，但未给出建模所依赖的噪声统计规律的具体获取方法。

发明内容

本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。

为此，本发明的一个目的在于提出一种电力***噪声自适应抗差状态估计方法（Adaptive Robust State Estimation，ARSE）。该方法能够通过统计学习获取噪声的分布规律，并将其与抗差状态估计方法进行在线匹配，从而可实现对各种噪声类型的自适应。

为了实现上述目的，本发明实施例的电力***噪声自适应抗差状态估计方法，包括以下步骤：S1，获取L个量测断面，其中，L为正整数；S2，对每个所述量测断面进行量测误差估计以获取误差矢量；以及S3，基于所述L个量测断面估计得到的误差矢量，运用统计学习方法估计得到广义高斯密度模型GGD（Generalized Gaussian Density）的参数，以获取量测噪声的分布类型，并根据所述分布类型选择对应的最优抗差状态估计模型。

在本发明的一个实施例中，所述S2包括：S21，令迭代计数器j=1；S22，对于第j个所述量测断面，假设所述第j个量测断面中具有m_j个量测量，基于所述m_j个量测量，运行最大指数绝对值状态估计方法MEAV，辨识出不良数据并剔除，以获取个正常量测量；S23，基于所述个正常量测量运行加权最小二乘法WLS，计算残差矢量S24，通过不良数据的估计辨识法，基于所述WLS得到的残差方程，由残差矢量估计量测方程中的误差矢量，并将所述误差矢量估计值表示为以及S25，若j<L，则令j=j+1，继续执行S22，否则继续执行S3。

在本发明的一个实施例中，所述S3包括：S31，对所述L个量测断面，获取与噪声矢量的所有维都对应的相同量测量，假设符合要求的共有G（G≤L）个相同量测量，将所述相同量测量放入一个集合，且表示为{α₁,α₂,L,α_G}，其中，α_i表示第i个正常量测量，令迭代计数器i=1；S32，对于α_i，在{ε₁,ε₂,…,ε_L}中找出与之对应的所有维，然后通过式 $\left\{\begin{array}{c}{\widehat{\mathrm{\&mu;}}}_i=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\&ap;0\\ {\widehat{\mathrm{\&tau;}}}_i=f_r^{-1}\left(\right[\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{|{\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i|}^r]/\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^{2r}})\\ {\widehat{\mathrm{\&sigma;}}}_i=\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^2}\&times;\sqrt{\frac{\mathrm{\&Gamma;}(1/{\widehat{\mathrm{\&tau;}}}_i)}{\mathrm{\&Gamma;}(3/{\widehat{\mathrm{\&tau;}}}_i)}}\end{array}\right.$ 估计GGD的三个参数，估计结果表示为其中，和分别是由第i个量测噪声对应的GGD中参数μ_i，τ_i和σ_i的估计值，μ_i为位置参数，τ_i为形状参数，σ_i为核宽度参数，j(1≤j≤L)是量测断面编号，L是量测断面的总数目，是的估计值，ε_ij代表量测误差ε_i在第j个量测断面中的数值；S33，若i<G，则令i=i+1，继续执行S32，否则继续执行S34；S34，对估计值取均值，作为GGD中参数的最终估计值，将最终估计值表示为S35，将S34得到的形状参数的平均值进行取整，以得到对应的取整值以及S36，根据所述形状参数的取整值选择对应的最优抗差状态估计模型，并取当前量测断面进行状态估计计算，以得到最优的状态变量估计值。

在本发明的实施例中，所述S35具体包括：将S34得到的形状参数的平均值代入式 $\left[\widehat{\mathrm{\&tau;}}\right]=\left\{\begin{array}{cc}1,& \mathrm{if}\widehat{\mathrm{\&tau;}}<1.5\\ 2,& \mathrm{if}1.5\&le;\widehat{\mathrm{\&tau;}}<3\\ \&infin;& \mathrm{if}3<\widehat{\mathrm{\&tau;}}\end{array}\right.$ 求得其取整值。

在本发明的一个实施例中，所述噪声分布类型包括拉普拉斯分布、高斯分布和均匀分布中的一种或多种。

在本发明的一个实施例中，所述抗差状态估计模型包括最大指数绝对值状态估计方法MEAV、指数型目标函数状态估计方法MES（Maximum Exponential Square）和以合格率最大为目标的状态估计方法MNMR（Maximum Normal Measurement Rate）中的一种或多种。

在本发明的一个实施例中，当所述噪声分布类型为拉普拉斯分布时，对应的最优抗差状态估计模型为MEAV；当所述噪声分布类型为高斯分布时，对应的最优抗差状态估计模型为MES；当所述噪声分布类型为均匀分布时，对应的最优抗差状态估计模型为MNMR。

根据本发明实施例的电力***噪声自适应抗差状态估计方法，可获取L个量测断面，并对每个量测断面进行量测误差估计以获取误差矢量，以及基于L个量测断面，可运用统计学习方法估计GGD模型的参数，以获取量测噪声的分布规律，并将其与最优抗差状态估计方法进行在线匹配，从而可实现对各种噪声类型的自适应，即在任何噪声分布类型下，均可以得到更接近于状态变量真值的最优估计结果。

本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中，

图1是根据本发明一个实施例的电力***噪声自适应抗差状态估计方法的流程图；

图2是GGD变化曲线（μ=0,ε=1.0）；

图3是反函数的变化曲线；

图4是一个2节点***算例；

图5是四种状态估计得到的2节点***中节点2电压幅值的估计值变化曲线；

图6是四种状态估计得到的2节点***中节点2电压相角的估计值变化曲线；

图7是四种状态估计得到的IEEE-118节点***中节点3电压幅值的估计值变化曲线；

图8是四种状态估计得到的IEEE-118节点***中节点3电压相角的估计值变化曲线；

图9是四种状态估计得到的IEEE-118节点***中节点92电压幅值的估计值变化曲线；以及

图10是四种状态估计得到的IEEE-118节点***中节点92电压相角的估计值变化曲线。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。相反，本发明的实施例包括落入所附加权利要求书的精神和内涵范围内的所有变化、修改和等同物。

在本发明的描述中，需要理解的是，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性。在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。此外，在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。

下面结合说明书附图详细说明本发明实施例的电力***噪声自适应抗差状态估计方法。

如图1所示，电力***噪声自适应抗差状态估计方法包括以下步骤：

步骤S1，获取L个量测断面，其中，L为正整数。

步骤S2，对每个量测断面进行量测误差估计以获取误差矢量。

在本发明的一个实施例中，发明人对每个量测断面进行量测误差估计以获取误差矢量的具体步骤归纳如下：

步骤S21，令迭代计数器j=1；

步骤S22，对于第j个量测断面，假设第j个量测断面中具有m_j个量测量，基于m_j个量测量，运行MEAV，辨识出不良数据并剔除，以获取个正常量测量。

步骤S23，基于个正常量测量运行WLS，然后计算残差矢量

步骤S24，通过不良数据的估计辨识法，基于WLS得到的残差方程，由残差矢量估计量测方程中的误差矢量，并将误差矢量估计值表示为

步骤S25，若j<L，则令j=j+1，继续执行步骤S22，否则继续执行步骤S3。

步骤S3，基于L个量测断面估计得到的误差矢量，运用统计学习方法估计得到广义高斯密度模型GGD的参数，以获取量测噪声的分布类型，并根据估计得到的噪声分布类型选择对应的最优抗差状态估计模型。

其中，对应于不同量测噪声分布类型下的抗差状态估计模型一般包括最大指数绝对值状态估计方法MEAV、指数型目标函数状态估计方法MES和以合格率最大为目标的状态估计方法MNMR等中的一种或多种。

在本发明的一个实施例中，运用统计学习方法估计得到广义高斯密度模型GGD的参数，以获取量测噪声的分布类型，并根据分布类型选择对应的最优抗差状态估计模型，即步骤S3具体包括：

步骤S31，对L个量测断面，获取与噪声矢量的所有维都对应的相同量测量，假设符合要求的共有G（G≤L）个相同量测量，将其放入一个集合，表示为{α₁,α₂,…,α_G}，其中，α_i表示第i个正常量测量，令迭代计数器i=1。

步骤S32，对于α_i，在{ε₁,ε₂,L,ε_L}中找出与之对应的所有维，然后通过

$\left\{\begin{array}{c}{\widehat{\mathrm{\&mu;}}}_i=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\&ap;0\\ {\widehat{\mathrm{\&tau;}}}_i=f_r^{-1}\left(\right[\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{|{\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i|}^r]/\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^{2r}})\\ {\widehat{\mathrm{\&sigma;}}}_i=\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^2}\&times;\sqrt{\frac{\mathrm{\&Gamma;}(1/{\widehat{\mathrm{\&tau;}}}_i)}{\mathrm{\&Gamma;}(3/{\widehat{\mathrm{\&tau;}}}_i)}}\end{array}\right.---\left(1\right)$

式估计GGD的三个参数，估计结果表示为其中，和分别是由第i个量测噪声对应的GGD中参数μ_i，τ_i和σ_i的估计值，μ_i为位置参数，τ_i为形状参数，σ_i为核宽度参数，j(1≤j≤L)是量测断面编号，L是量测断面的总数目，是的估计值，ε_ij代表量测误差ε_i在第j个量测断面中的数值。

步骤S33，若i<G，则令i=i+1，继续执行步骤S32，否则继续执行步骤S34。

步骤S34，对估计值取均值，作为GGD中参数的最终估计值，将最终估计值表示为

步骤S35，将步骤S34得到的形状参数的平均值进行取整，以得到对应的取整值

在本发明的实施例中，步骤S35具体包括：将步骤S34得到的形状参数的平均值代入

$\left[\widehat{\mathrm{\&tau;}}\right]=\left\{\begin{array}{cc}1,& \mathrm{if}\widehat{\mathrm{\&tau;}}<1.5\\ 2,& \mathrm{if}1.5\&le;\widehat{\mathrm{\&tau;}}<3\\ \&infin;& \mathrm{if}3<\widehat{\mathrm{\&tau;}}\end{array}\right.---\left(2\right)$

式求得其取整值

步骤S36，根据形状参数的取整值选择对应的最优抗差状态估计模型，并对当前量测断面进行状态估计计算，以得到最优的状态变量估计值。

例如，当即噪声分布类型为拉普拉斯分布时，选择对应的状态估计方法MEAV，当即噪声分布类型为高斯分布时，选择对应的状态估计方法MES，当即噪声分布类型为均匀分布时，选择对应的状态估计方法MNMR，并对当前量测断面进行状态估计计算，以得到最优的状态变量估计值。

为了使得本领域的技术人员更加地理解本发明，下面详细介绍ARSE具体实现方式。

以下首先给出一种电力***抗差状态估计一般模型，其既可以从理论上统一多种已有的抗差状态估计，又可以导出新的状态估计方法，之后给出基于统计学习的状态估计量测噪声分布类型的估计方法，最后给出ARSE的实施方式和实施步骤。

（1）电力***抗差状态估计一般模型

一般地，状态估计的量测方程可以表示为

z=h(x)+ε （3）

其中，z∈R^m为量测矢量，可包括支路功率量测、节点注入功率量测和节点电压幅值量测等，PMU量测中还包括相角量测；x∈Rⁿ为包括所有节点电压幅值和相角（参考节点相角除外）的状态变量，n=2N-1，N为网络中节点的总数目；h:Rⁿ→R^m即量测表达式，为状态变量到量测矢量的非线性映射；ε∈R^m为量测误差矢量。

状态估计量测噪声的实际分布类型是未知的。一般地，量测噪声类型可用如下形式的广义高斯密度模型GGD来表示

$G_{\mathrm{\&sigma;}}\left(\mathrm{\&epsiv;}\right)=\frac{\mathrm{\&tau;}}{2\mathrm{\&sigma;\&Gamma;}(1/\mathrm{\&tau;})}\exp (-{\left(\frac{|\mathrm{\&epsiv;}-\mathrm{\&mu;}|}{\mathrm{\&sigma;}}\right)}^{\mathrm{\&tau;}})---\left(4\right)$

其中， 反映了噪声水平，为对应量测量的标准差，Γ(·)为由下式定义的Gamma函数

$\mathrm{\&Gamma;}\left(s\right)={\&Integral;}_0^{\&infin;}{x}^{s-1}{e}^{-x}\mathrm{dx},s>0---\left(5\right)$

其中，τ为形状参数，σ为核宽度参数，μ为位置参数，这三个参数应满足以下条件：

$\left\{\begin{array}{c}{G}_{\mathrm{\&sigma;}}\left(\mathrm{\&epsiv;}\right)\&GreaterEqual;0\\ {\&Integral;}_{-\&infin;}^{\&infin;}{G}_{\mathrm{\&sigma;}}\left(\mathrm{\&epsiv;}\right)\mathrm{d\&epsiv;}=1\\ G_{\mathrm{\&sigma;}}\left(\mathrm{\&epsiv;}\right)=G_1(\mathrm{\&epsiv;}/\mathrm{\&sigma;})/\mathrm{\&sigma;}\end{array}\right.---\left(6\right)$

对以上三个参数，形状参数最为重要，因为该形状参数决定了量测噪声的分布类型；对于状态估计来说，量测误差的均值为0，因而位置参数可视为0，此时GGD的变化曲线如图1所示。

将噪声的一般模型GGD应用于电力***状态估计问题中，并对量测噪声分布类型进行非参数估计，可得如下形式的电力***状态估计模型

$\mathrm{MaxJ}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&tau;}}{2\mathrm{\&sigma;\&Gamma;}(1/\mathrm{\&tau;})}\exp (-{\left(\frac{|z_i-h_i\left(x\right)|}{\mathrm{\&sigma;}}\right)}^{\mathrm{\&tau;}})---\left(7\right)$

模型（7）适应于任何噪声分布类型，它有多种具体的形式。严格的理论证明和实际试验都表明，模型（7）具有很好的抗差性。因此，模型（7）可作为电力***抗差状态估计一般模型。下面来说明模型（7）可从理论上统一已有的几种抗差状态估计方法，并可以导出新的抗差状态估计模型。

当GGD中参数取μ=0,τ=2，即假设量测噪声符合高斯分布时，模型（7）变为如下形式

$\mathrm{MaxJ}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{\sqrt{\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp (-\frac{{(z_i-h_i\left(x\right))}^2}{{\mathrm{\&sigma;}}^2})---\left(8\right)$

模型（8）等价于国内学者提出的指数型目标函数状态估计方法MES。

当GGD中参数取μ=0,τ=1，即假设量测噪声符合拉普拉斯分布时，模型（7）变为如下形式

$\mathrm{MaxJ}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{2\mathrm{\&sigma;}}\exp (-\frac{|z_i-h_i\left(x\right)|}{\mathrm{\&sigma;}})---\left(9\right)$

模型（9）等价于本发明人提出的最大指数绝对值状态估计方法MEAV。

满足式（6）的均匀噪声（μ=0,τ→+∞）的密度函数可表示如下

$G_{\mathrm{\&sigma;}}\left(\mathrm{\&xi;}\right)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\&sigma;}},& \mathrm{if}\left|\mathrm{\&xi;}\right|\&le;\mathrm{\&sigma;}/2\\ 0,& \mathrm{if}\left|\mathrm{\&xi;}\right|>\mathrm{\&sigma;}/2\end{array}\right.\left(10\right)$

当GGD采用式（10）描述的均匀噪声分布类型时，模型（10）变为如下形式

$\mathrm{MaxJ}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}f(z_i-h_i\left(x\right))$

$s.t.\left\{\begin{array}{cc}f(z_i-h_i\left(x\right))=1/{\mathrm{\&sigma;}}_i,& \mathrm{if}|z_i-h_i\left(x\right)|\&le;{\mathrm{\&sigma;}}_i/2\\ f(z_i-h_i\left(x\right))=0,& \mathrm{if}|z_i-h_i\left(x\right)|>{\mathrm{\&sigma;}}_i/2\end{array}\right.---\left(11\right)$

其中，窗宽参数σ_i代表第i个量测量的精度，精度越高，其值越小。

如果认为σ_i/2是第i个量测量的扩展不确定度，满足|z_i-h_i(x)|≤σ_i/2的量测被认为是合格的量测，则模型（11）的目标函数转化为求解最大合格率，此时，模型（11）等价于国内学者提出的以合格率最大为目标的状态估计方法MNMR。

以上分析表明，模型（7）可从理论上统一多种已有的抗差状态估计方法；如果模型（7）的参数取其他的数值（区别于以上三种状态估计模型对应的参数），则由模型（7）可推导得到其他状态估计模型。

（2）噪声分布类型的估计方法

虽然模型（7）具有很好的抗差性，但是其中含有两个未知的参数。在电力***中，如果量测噪声的分布类型已知，即模型（7）中的两个参数已知，则可通过求解模型（7）得到状态变量的最优估计值。但是，实际中状态估计中的量测噪声分布类型是未知的，此时应用模型（7）进行状态估计计算就需要利用大量的量测断面以对其中的两个未知参数进行估计。一旦获得了这两个参数的估计值，模型（7）即可被确定，并对新的量测断面具有适应性，即在新的量测断面下可获得最优的状态变量估计值，从而实现了对量测噪声分布的自适应。

以下给出通过状态估计量测断面估计模型（7）中参数μ、τ及σ的方法。

对于模型（7），噪声的r阶绝对中心矩由下式给出

$\begin{array}{c}{d}_r=E\left[{|\mathrm{\&xi;}-\mathrm{\&mu;}|}^r\right]={\&Integral;}_{-\&infin;}^{+\&infin;}{|\mathrm{\&xi;}-\mathrm{\&mu;}|}^{r}p\left(\mathrm{\&xi;}\right)\mathrm{d\&xi;}\\ ={\mathrm{\&sigma;}}^r\mathrm{\&Gamma;}((r+1)/\mathrm{\&tau;})/\mathrm{\&Gamma;}(1/\mathrm{\&tau;})\end{array}---\left(12\right)$

则GGD模型的2r阶绝对中心矩为

d_2r=σ^2rΓ((2r+1)/τ)/Γ(1/τ) （13）

由式（12）和式（13）可得

$\frac{d_r}{\sqrt{d_{2r}}}=\frac{\mathrm{\&Gamma;}((r+1)/\mathrm{\&tau;})}{\sqrt{\mathrm{\&Gamma;}((2r+1)/\mathrm{\&tau;})\mathrm{\&Gamma;}(1/\mathrm{\&tau;})}}---\left(14\right)$

式（14）的右边是形状参数τ的函数，不妨记为f_r(τ)，则有

$\mathrm{\&tau;}=f_r^{-1}(d_r/\sqrt{d_{2r}})---\left(15\right)$

其中，是f_r(τ)的反函数，可通过数值计算方法获得的数值，如图3所示，图3给出了当r=1和2时，的变化曲线。

对于量测方程（3）中的量测误差矢量ε的第i维ε_i(1≤i≤m)，基于矩匹配法的GGD参数估计公式为

$\left\{\begin{array}{c}{\widehat{\mathrm{\&mu;}}}_i=\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\&ap;0\\ {\widehat{\mathrm{\&tau;}}}_i=f_r^{-1}\left(\right[\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{|{\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i|}^r]/\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^{2r}})\\ {\widehat{\mathrm{\&sigma;}}}_i=\sqrt{\frac{1}{L}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{({\mathrm{\&epsiv;}}_{\mathrm{ij}}-{\widehat{\mathrm{\&mu;}}}_i)}^2}\&times;\sqrt{\frac{\mathrm{\&Gamma;}(1/{\widehat{\mathrm{\&tau;}}}_i)}{\mathrm{\&Gamma;}(3/{\widehat{\mathrm{\&tau;}}}_i)}}\end{array}\right.---\left(1\right)$

其中，和分别是由第i个量测噪声对应的GGD中参数μ_i，τ_i和σ_i的估计值；j(1≤j≤L)是量测断面编号，L是量测断面的总数目，是的估计值，ε_ij代表量测误差ε_i在第j个量测断面中的数值。

通过式（1）可得到量测误差矢量ε的每一维对应的GGD中三个参数的估计值，GGD中三个参数的最终估计值可取ε的所有维得到的这三个参数估计值的平均值。需要说明的是，在本发明的实施例中，由于在一个较长的时期内，量测误差分布类型并不发生改变，因此量测误差分布类型并不需要在ARSE实施的每一次都予以估计。

值得指出的是，式（1）通过大量量测断面中的误差矢量以估计量测噪声的实际分布类型，然而实际量测误差矢量是未知的。一种可行的解决方案是通过状态估计的残差来估计误差矢量。其具体算法为：对每个量测断面，运行MEAV辨识不良数据，并删除不良数据；基于正常量测运行WLS，并计算残差矢量；基于WLS得到的残差方程，利用残差矢量来估计误差矢量，并采用同样的方法估计出所有L个量测断面的误差矢量，然后即可通过式（1）估计量测噪声的分布类型。

（3）ARSE的实现方法

需要说明的是，在理论上说，一旦由上面给出的方法获得了GGD中三个参数的估计值，即可通过将这三个参数带入式（7）得到具体的状态估计模型，从而可实现量测噪声分布类型的自适应，即建立了精确的ARSE模型和方法。但是由式（1）得到的形状参数的估计值经常不是整数，其对应的ARSE模型难以求解。为了解决这个问题，下面提出一种可行的ARSE实现方法。

在本发明的一个实施例中，噪声分布类型一般认为包括拉普拉斯分布、高斯分布和均匀分布等中的一种或多种。具体地，由图2可见，当形状参数τ大于3时，对应的GGD曲线趋近于均匀分布曲线（τ→∞）；当形状参数τ小于1时，对应的GGD曲线趋近于拉普拉斯分布曲线。考虑到高斯分布、拉普拉斯分布和均匀分布是最为常见的三种噪声分布类型，并且它们对应的抗差状态估计方法（分别是MES、MEAV和MNMR）容易求解，因此，可用高斯分布、拉普拉斯分布和均匀分布三种分布类型的某一种来近似表示任何一种具体的噪声分布类型，并从MES、MEAV和MNMR中选择与噪声分布类型相适应的最优状态估计方法，即可近似实现ARSE。

在本发明的实施例中，上面给出的ARSE方法的本质是对形状参数的估计值进行取整，其取整方法为

$\left[\widehat{\mathrm{\&tau;}}\right]=\left\{\begin{array}{cc}1,& \mathrm{if}\widehat{\mathrm{\&tau;}}<{\mathrm{\&beta;}}_1\\ 2,& \mathrm{if}{\mathrm{\&beta;}}_1\&le;\widehat{\mathrm{\&tau;}}<{\mathrm{\&beta;}}_2\\ \&infin;& \mathrm{if}{\mathrm{\&beta;}}_2<\widehat{\mathrm{\&tau;}}\end{array}\right.---\left(16\right)$

其中，代表形状参数估计值的取整值；是式（16）得到的的均值；β₁和β₂是两个门槛参数，推荐值为β₁=1.5,β₂=3。

由以上分析可知，在本发明的实施例中，当GGD分布为拉普拉斯分布时，对应的最优抗差状态估计模型是MEAV；当GGD分布为高斯分布时，对应的最优抗差状态估计模型是MES；当GGD分布为均匀分布时，对应的最优抗差状态估计模型是MNMR；如果将状态估计中的量测噪声分布类型（由GGD中三个参数决定）与抗差状态估计一般模型（7）进行在线匹配，即可实现ARSE。

（4）ARSE的具体实现步骤

以上给出ARSE的具体实现步骤归纳如下：

Step1：（初始化）任取L个连续的断面，用S1,S2,…S_L表示，设j=1；

Step2:（估计量测误差）对于第j个量测断面，假设断面中有m_j个量测量；

Step2.1：（不良数据辨识）基于m_j个量测量，运行MEAV，辨识出不良数据并剔除，假设正常量测的数量为

Step2.2:(残差计算)基于个正常量测运行WLS，然后计算残差矢量 $r_j={[r_{j\_1},r_{j\_2},...,r_{j\_m_j^g}]}^T;$

Step2.3:（噪声估计）通过不良数据的估计法，基于WLS得到的残差方程，由残差矢量估计量测方程中的误差矢量，将误差矢量估计值表示为

Step3:若j<L，则令j=j+1，转Step2；否则转Step4；

Step4:（找出相同的量测量）对所有的L个量测断面，找出与噪声矢量的所有维都对应的相同量测量，假设符合要求的共有G（G≤L）个相同量测量，将其放入一个集合，表示为{α₁,α₂,…,α_G}，其中，α_i表示第i个正常量测量，令迭代计数器i=1；

Step5:（GGD参数估计）对于αi，在{ε₁,ε₂,…,ε_L}中找出与之对应的所有维，然后利用式（1）估计GGD的三个参数，估计结果表示为

Step6:若i<G，则令i=i+1，转Step5，否则转Step7；

Step7:（参数平均）对估计值取均值，作为GGD中参数的最终估计值，将最终估计值表示为

Step8:（形状参数取整）将Step7得到的形状参数的平均值代入式（16）求得其取整值，表示为

Step9:（选择状态估计模型）根据参数在MNMR、MES及MEAV中选择相应的状态估计模型，然后取当前断面进行状态估计计算，得到最优的状态变量估计值；

Step10:END。

为使本领域技术人员更好地理解本发明以及了解本发明相对现有技术的优点，发明人结合具体实施例进行进一步的阐释。

设定利用一个2节点***和IEEE-118节点***来测试所提出的噪声自适应状态估计方法（ARSE）的性能。测试环境为PC机，CPU为Intel(R)Core(TM)i3M370、主频为2.40GHz、内存2.00GB。

1）初步测试

（1）抗差性测试

试验在图4所示的两节点***进行，其中，节点1为参考节点，节点2的电压幅值和相角是待求的状态变量。网络参数和真实潮流在图4中也给出，在所有的测试中，节点2的电压幅值保持为1，相角保持为-π/6。

在试验中，在真实潮流的基础上分别独立叠加三种噪声分布类型（如高斯分布噪声、拉普拉斯分布噪声以及均匀分布噪声）中的一种，噪声的标准差是0.005p.u.，同时在节点1的注入有功量测p₁和支路有功量测p₁₂的真值上迭代50%的误差以仿真不良数据，试验采用全量测。然后分别用状态估计方法MNMR、MES和MEAV进行估计。为了比较在不良的量测噪声下，以上三种状态估计方法的估计性能，定义下面两个均方误差（Mean Square Error，MSE）指标

$\left\{\begin{array}{c}{\mathrm{MSE}.v}_2=\frac{1}{T}{\mathrm{\&Sigma;}}_{i=1}^T{({\hat{v}}_{2.i}-v_{2.\mathrm{true}})}^2\\ {\mathrm{MSE}.\mathrm{\&theta;}}_2=\frac{1}{T}{\mathrm{\&Sigma;}}_{i=1}^T{({\widehat{\mathrm{\&theta;}}}_{2.i}-{\mathrm{\&theta;}}_{2.\mathrm{true}})}^2\end{array}\right.---\left(17\right)$

其中，v_2.true=1和θ_2.true=-π/6分别是节点2电压幅值和相角的真值，和分别为其估计值；T是总试验次数；MSE.v₂和MSE.θ₂分别为节点2电压幅值和相角估计值的MSE指标。

在1000次试验中，在不同的噪声分布类型下（如高斯分布、拉普拉斯分布以及均匀分布），由状态估计方法MNMR、MES和MEAV分别得到的节点2电压幅值和相角估计值的MSE指标如下面表1所示（为表示方便，表1中省略了单位1e-6）。

表1三种状态估计方法在不同噪声下得到的MSE指标

Table1MSE Indicators of MNMR,MES and MEAV on Uniform noise,Laplacenoise and Gaussian noise

由表1可知，由状态估计方法MNMR、MES和MEAV得到的估计结果不受不良数据的影响，由此，证明了它们具有良好的抗差性；由表1同样可知，当噪声分布类型分别是均匀分布、高斯分布和拉普拉斯分布时，由状态估计方法MNMR、MES和MEAV得到的估计结果的MSE指标最小，因而，MNMR、MES和MEAV分别是当噪声分布类型为均匀分布、高斯分布和拉普拉斯分布时的最优抗差状态估计方法，从而证明了进行噪声自适应状态估计的必要性。

（2）GGD参数估计方法性能测试

在本发明的实施例中，GGD中三个参数的估计精度对于ARSE的正确实施至关重要，因此需测试所提出的GGD参数估计方法的性能。

同样考虑图4所示的2节点***，在潮流解的基础上分别独立叠加高斯分布、拉普拉斯分布和均匀分布的噪声（其中噪声的标准差为0.005p.u.）。分别取100、1000和10000个量测断面，在试验中，首先通过状态估计方法MEAV进行估计，以辨识不良数据；之后通过WLS进行估计；当WLS收敛之后，基于估计辨识法通过残差来估计量测误差（即量测噪声）；然后通过公式（1）估计GGD参数，最终通过公式（16）对形状参数取整。三种噪声分布下，GGD参数的估计结果分别如表2、3和4所示。

表2高斯噪声时GGD参数估计结果

Table2Estimation Results Of The Parameters Of GGD When The Noise IsGaussian

表3拉普拉斯噪声时GGD参数估计结果

Table3Estimation Results Of The Parameters Of GGD When The Noise IsLaplace

表4均匀噪声时GGD参数估计结果

Table4Estimation Results of The Parameters Of GGD When The Noise IsUniform

由表2、3和4可知，随着量测断面数目的增加，位置参数估计值的绝对值逐渐减小，并趋近于0，这与理论结果相符；核宽度估计值随着量测断面数目的增加而增大；在试验中，一旦得到的估计值，即可由图3可得形状参数的估计值（试验中取r=1），并得到然后由式（16）求得形状参数的取整估计值根据取整估计值以得到相应的噪声分布类型。由表2、表3和表4可知，随着量测断面数目的增加，形状参数的估计值趋近于其真值，而形状参数的估计值则总是等于其真值，这就证明了所提出的噪声类型估计方法可准确地估计出噪声的正确分布类型，从而为实现噪声自适应状态估计奠定了基础。

2）自适应性测试

（1）2节点***测试结果

下面将测试本发明所提出的ARSE的估计性能。同样考虑图4所示的2节点***。在50次试验中的每一次试验中，分别在潮流的基础上随机叠加三种噪声分布类型（如均匀分布、高斯分布和拉普拉斯分布）中的其中一种，同时在节点1的注入有功量测p₁和支路有功量测p₁₂的真值上迭代50%的误差以仿真不良数据，试验采用全量测。每次试验均包含1000个量测断面，首先由所提出的量测噪声类型分布方法估计噪声的分布类型；然后分别通过状态估计方法MES、MNMR、MEAV以及本发明提出的ARSE进行估计。在50次独立试验中，节点2电压幅值和相角的估计结果分布如图5、图6所示。

由图5和图6可知，在50次试验中，MES、MNMR和MEAV的总体估计性能相当，而在单次试验中，其估计结果不同，具体地说，MES、MNMR和MEAV分别在高斯噪声、均匀噪声与拉普拉斯噪声下得到的估计结果更接近状态变量的真值。而ARSE通过对噪声类型的估计实现了在MES、MNMR和MEAV中的“自动切换”。从总体估计性能上看，由ARSE得到的估计结果更接近于状态变量的真值，从而证明了本发明所提出的噪声自适应状态估计方法ARSE的有效性。

（2）IEEE-118节点***测试结果

进一步在IEEE-118节点***测试ARSE的估计性能。量测数据的产生方法与2节点的测试方法相同，同时随机地叠加2%的一致性不良数据。然后分别通过MES、MNMR、MEAV以及本发明提出的ARSE进行估计。在50次独立试验中，节点3和节点92电压幅值和相角的估计结果分布如图7、图8、图9、图10所示。其中，节点3和节点92电压幅值和相角的真值分别是（0.968，-0.3218）和（0.993，0.0663）。

由图7、图8、图9、图10可知，在50次试验中，MES、MNMR和MEAV的总体估计性能相当，在单次试验中，其估计结果不同；MES、MNMR和MEAV分别在高斯噪声、均匀噪声与拉普拉斯噪声下得到的估计结果更接近状态变量的真值；而ARSE通过对噪声类型的估计实现了在MES、MNMR和MEAV中的“自动切换”；从总体估计性能上看，由ARSE得到的估计结果更接近于状态变量的真值。事实上，其他节点也有类似的估计结果，从而进一步证明了本发明提出的噪声自适应状态估计ARSE的有效性。

综上所述，本发明提出的ARSE能够通过统计学习获取噪声的分布规律，并将其与抗差状态估计方法进行在线匹配，从而可实现对各种噪声类型的自适应，即在任何噪声分布类型下，均可以得到更接近于状态变量真值的最优估计结果。

应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行***执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列（PGA），现场可编程门阵列（FPGA）等。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (5)

1.一种电力***噪声自适应抗差状态估计方法，其特征在于，包括以下步骤：

S1，获取L个量测断面，其中，L为正整数；

S2，对每个所述量测断面进行量测误差估计以获取误差矢量，包括：

S21，令迭代计数器j＝1；

S22，对于第j个所述量测断面，假设所述第j个量测断面中具有m_j个量测量，基于所述m_j个量测量，运行最大指数绝对值状态估计方法MEAV，辨识出不良数据并剔除，以获取个正常量测量；

S23，基于所述个正常量测量运行加权最小二乘法WLS，计算残差矢量

S24，通过不良数据的估计辨识法，基于所述WLS得到的残差方程，由残差矢量估计量测方程中的误差矢量，并将所述误差矢量估计值表示为以及

S25，若j<L，则令j＝j+1，继续执行S22，否则继续执行S3；以及

S3，基于所述L个量测断面估计得到的误差矢量，运用统计学习方法估计得到广义高斯密度模型GGD的参数，以获取量测噪声的分布类型，并根据所述噪声分布类型选择对应的最优抗差状态估计模型，包括：

S31，对所述L个量测断面，获取与噪声矢量的所有维都对应的相同量测量，符合要求的共有G个相同量测量，其中，G≤L，将所述相同量测量放入一个集合，且表示为{α₁,α₂,…,α_G}，其中，α_i表示第i个正常量测量，令迭代计数器i＝1；

S32，对于α_i，在{ε₁,ε₂,…,ε_L}中找出与之对应的所有维，然后通过式估计GGD的三个参数，估计结果表示为其中， 和分别是由第i个量测噪声对应的GGD中参数μ_i，τ_i和σ_i的估计值，μ_i为位置参数，τ_i为形状参数，σ_i为核宽度参数，j是量测断面编号，L是量测断面的总数目，其中，1≤j≤L，是的估计值，其中，d_r为r阶绝对中心矩，d_2r为2r阶绝对中心矩，r＝1，为d_r的估计值，为d_2r的估计值，ε_ij代表量测误差ε_i在第j个量测断面中的数值；

S33，若i<G，则令i＝i+1，继续执行S32，否则继续执行S34；

S34，对估计值取均值，作为GGD中参数的最终估计值，将最终估计值表示为

S35，将S34得到的形状参数的平均值进行取整，以得到对应的取整值以及

S36，根据所述形状参数的取整值选择对应的最优抗差状态估计模型，并取当前量测断面进行状态估计计算，以得到最优的状态变量估计值。

2.如权利要求1所述的电力***噪声自适应抗差状态估计方法，其特征在于，所述S35具体包括：

将S34得到的形状参数的平均值代入式求得其取整值。

3.如权利要求1所述的电力***噪声自适应抗差状态估计方法，其特征在于，所述噪声分布类型包括拉普拉斯分布、高斯分布和均匀分布中的一种或多种。

4.如权利要求1所述的电力***噪声自适应抗差状态估计方法，其特征在于，所述对应的最优抗差状态估计模型包括最大指数绝对值状态估计方法MEAV、指数型目标函数状态估计方法MES和以合格率最大为目标的状态估计方法MNMR中的一种或多种。

5.如权利要求3或4所述的电力***噪声自适应抗差状态估计方法，其特征在于，所述根据所述噪声分布类型选择对应的最优抗差状态估计模型，包括：

当所述噪声分布类型为拉普拉斯分布时，选择所述对应的最优抗差状态估计模型为MEAV；

当所述噪声分布类型为高斯分布时，选择所述对应的最优抗差状态估计模型为MES；

当所述噪声分布类型为均匀分布时，选择所述对应的最优抗差状态估计模型为MNMR。