最大概率（精选4篇）
最大概率 篇1
1 算法概述
首先建立图像退化模型, 低分辨率图像和原图像之间的关系为YK=HKX+NK, 其中HK为退化矩阵, NK为加性零均值高斯噪声。
根据贝叶斯公式, 在已知序列低分辨率图像Y的情况下, 原图像X的概率表达为:
其中P (X|Y) 为后验概率, P (Y|X) 为高分辨率图像退化为低分辨率序列图像Y时的条件概率, P (X) 和P (y) 分别表示高分辨率图像和退化图像的先验概率。式中的最优估计X代表在已知低分辨率序列图像时最可能的高分辨率图像X。
由于项P (Y) 是对结果X无影响的量, 因此可将上式简化为:
由此得到目标方程:
然后利用有约束的最小二乘法来最小化Lagrange方程, 对方程进行求解。
2 基于特征匹配的运动参数估计
为了获取退化矩阵, 需要对图像序列进行配准, 使所有图像的灰度位于相同坐标系。这里我们引入图像的仿射变换公式, 仿射变换在笛卡尔坐标系上定义了原图像的坐标点 (x, y) 在图像旋转、平移和缩放后对应的坐标点 (x’, y’) 与其自身的关系。若设坐标原点在图像的左上角, 图像围绕坐标远点旋转后再作平移, 旋转的角度为, 水平和垂直平移量分别为h和v, 仿射变换公式即 (x, y) 与 (x’, y’) 的关系如下图所示:
该表达式显然就是我们所期待的求G和F运动关系的式子。但用该式直接求出的值并不是正确的结果, 因为在之前的表达式中忽略了高次残差项。为了得到正确的值, 则应当用结果 (h, v, ) 进行迭代。具体的做法是:求得初始值 (h0, v0, θ0) , 让原始图像F按 (h, v, θ) 进行旋转平移得到F1, 再把它带入式中计算得到 (h1, v1, θ1) , 将 (h0, v0, θ0) 累加到 (h1, v1, θ1) , 再让F按 (h1, v1, θ1) 进行旋转平移得到F2, 如此迭代循环, 就可得到期望精度估计值。由于这个过程实际上是让F不断向G配准的过程, 因此, 这种参数估计方法是基于图像配准的运动参数估计方法。
3 矩阵运算的简化
阅读文献的过程中我们发现, 通常在求解退化矩阵时, 为了使得矩阵的大小在运算过程中能够匹配, 常将高分辨率图像和低分辨率图像的矩阵垒叠成MN×1 (M、N分别为图像的行、列数) 的长矩阵来进行运算。这种思路虽然简单, 但是在进行矩阵的乘除法运算时会消耗大量存储空间, 处理较大的图片时, 会出现内存溢出的情况。为此, 我们在这里提出一种新的方法, 用以解决上述问题。
假设M、N、q分别为图像矩阵的行数、列数以及图像的放大倍数, 那么低分辨率图像Y、高分辨率图像X分别为M×N、qM×qN大小的矩阵。而在式||Y-HKXC||因为矩阵大小不匹配, 无法用M×N、qM×qN大小的矩阵进行运算。只有把原图Y与高分辨率图像转化为列数为1的向量才能进行运算, 而此时求得的退化矩阵是MN×q2MN的超大矩阵, 图片较大时, 会出现内存不够的情况。
在上述问题的描述中, YK、XC矩阵的大小不等是产生问题的根本原因。为此我们可以将Y采用最临近差值放大q倍成高分辨率图像Yl, 以Yl、XC来求解新的运动参数矩阵Hkl, 此时矩阵大小为qM×qN的矩阵, Hkl不再是采样、运动估计矩阵的乘积, 而仅代表运动估计矩阵。
4 实验结果
从双线性差值法与本文算法之间的比较可以看出, 本文算法获得的高分辨率图像具有较好的边缘特性。
参考文献
[1]刘晓天.基于MAP技术的图像超分辨率复原研究与实现[D].国防科学技术大学:软件工程, 2004.
[2]许静.基于MAP技术图像超分辨率重建的研究[D].中国海洋大学:通信与信息系统, 2007.
[3]刘雪婷.基于MAP算法的图像超分辨率重构技术研究[D].山东大学:通信与信息系统, 2008.
最大概率 篇2
关键词：可靠度计算,最大熵拟合,概率密度函数
结构可靠度计算常用的方法有一次二阶矩法[1]，有条件则采用高次高阶矩法以提高计算精度，如二次二阶矩法与二次四阶矩法，这类方法需要计算功能函数的偏导数，计算采用迭代法，计算结果为可靠指标，不是直接的失效概率，而且随着功能函数的非线性程度的增大，计算精度有所减小。结构可靠度的计算也采用蒙特卡罗法，计算量大，一般用于各种可靠度近似方法的校核。当极限状态函数有显式表达式而且随机变量不多的简单情况，亦可采用数值积分法，但当失效区域是一般不规则或是多连通或不连通的独立区域时，数值积分法亦有一定困难。如果可以求得极限状态函数的前几阶矩(一般为四阶矩或五阶矩)，则可采用拟合概率密度函数方法近似计算失效概率。拟合的函数有最大熵密度函数[2]和多项式密度函数[3]，前者物理意义明显，计算较复杂，后者没有直接的物理背景，计算较简单。
本文对极限状态函数求得前几阶矩，再用最大熵法拟合其概率密度函数，近似计算失效概率。
1 近似计算失效概率的步骤
1.1 极限状态函数的前几阶矩
已知极限状态函数为Z=g(X)=g(x1,x2，…，xn)，随机向量X的概率密度函数为f(X)，则Z的各阶原点矩为:
式(1)可通过数值积分求得。f(X)可视为权函数，采用高斯积分法，计算效率和精度较高。如正态分布、指数分布、均匀分布，相应的权函数为e-x2,e-x,1，可分别选高斯·埃尔米特(Gauss-Hermite)、高斯·拉盖尔(Gauss-Laguerre)、高斯·勒让德(Gauss-Legendre)积分点。
显然，这里规定m0=1是合理的。
1.2 用概率分布函数的多项式逼近极限状态函数
根据最大熵原理，极限状态函数的概率密度函数为:
式(2)中a0,a1，…，aN为待定系数;N为矩的最大阶数，一般取4～5阶矩即能获得良好的效果[2]。
这里，z的各阶原点矩也可采用积分形式表示为:
式(3)中待定系数有N+1个，方程有N+1个。
式(3)为非线性方程组，可转化为无约束优化问题，采用优化算法求解:
式(4)中ωk为加权系数，为使式(4)各项的量纲一致，可取ωk=m-k2。
优化算法的初始点可按正态分布给出，即:
式(5)中μz=m1，σz2=m2-m12。于是初始点为:
解出待定系数，就可以得到概率密度函数。
1.3 计算失效概率
对式(2)积分得到极限状态函数的失效概率:
注意，实际问题的解应与用中心点法计算的可靠指标的近似值β≈μz/σz，按正态分布换算的失效概率Φ-1(-μz/σz)相差不会很大。
2 简单算例
例1:已知非线性极限状态方程g=567fr+0.5H2=0。f服从正态分布，μf=0.6，δf=0.131;r服从正态分布，μr=2.18，δr=0.03;H服从对数正态分布，μH=32.8，δH=0.03。试求失效概率[1]。
1)用数值积分求极限状态函数的前4阶原点矩:
2)经计算，待定系数初始值为:
用优化算法求出待定系数:
3)求失效概率。这里的积分下界限取为μz-5σz=-320.870 9。失效概率Pf=0.025 951(换算为可靠指标β=1.944 0)。
作为对比，采用4阶原点矩计算的结果见表1。
例2:已知极限状态方程g=1.8-x1-x2=0，其中，x1,x2∈[0,1]，其联合概率密度函数为f(x1,x2)=2-x1-x2。试求失效概率[1]。
计算方法同例1，结果列于表1中。
注:1)例1的精确解由蒙特卡罗法求得,例2的精确解由直接积分求得;2)带*数为按正态分布逆变换β=Φ-1(1-Pf)换算得到的
3 结语
本文采用数值积分求得极限状态函数的前几阶矩，拟合极限状态函数的概率密度，求解失效概率。计算表明:1)采用四阶矩与五阶矩的计算结果相差不大，所以实用中可以采用四阶矩。2)数值积分的界限一般取为μ±5σ，即可满足工程精度要求。3)与多项式拟合法相比，计算精度变化较小，数值较稳定，所以适用性较好。
参考文献
[1]赵国藩,金伟良,贡金鑫.结构可靠度理论[M].北京:中国建筑工业出版社,2000.
[2]章光,朱维申,白世伟.计算近似失效概率的最大熵密度函数法[J].岩石力学与工程学报,1995(8):18-19.
[3]邓建,李夕兵,古德生.结构可靠性分析的多项式数值逼近法[J].计算力学学报,2002(11):26-30.
[4]李庆扬,王能超,李大义.数值分析[M].武汉:华中工学院出版社,1982.
最大概率 篇3
1 系统模型与定位算法
1.1 定位系统模型
定位系统模型如图1所示, 在LLA (经纬高) 坐标系中, T为地面固定辐射源, 位置矢量为Xt, 辐射源信号波长为λ;卫星位置矢量为Xr;基线长度为d (dλ) , 基线矢量方向由天线2指向天线1, 基线矢量为Xb, ωr为基线转角速度, 那么辐射源来波方向与基线的夹角θ有[7]
理论上讲, 天线1与天线2接收信号的无模糊相位差为
由于干涉仪鉴相器输出范围只能在 (-π, π) 内, 致使干涉仪输出相位差存在模糊, 干涉仪输出的模糊相位差φ与无模糊相位差φr满足[8]:
式 (3) 中n为模糊数, 是由基线长度、辐射源信号波长、基线矢量、定位平台和辐射源位置等共同决定。n的存在导致了测量θ模糊, 进而导致定位模糊, 从而给精确快速定位带来更大困难。
1.2 基于最大后验概率的定位算法
根据定位站与地面的相对位置和干涉仪性能参数, 可以确定目标搜索区域S的大致范围。假设定位站是卫星, 并且在ECF坐标系下卫星位置矢量:E= (x, y, z) , 基线方向矢量:M= (mx, my, mz) 。目前, 利用GPS测量卫星位置矢量和基线方向矢量已经达到了较高的精度, 满足实际应用要求[9]。取目标估计区域内的任意一点T= (Tx, Ty, Tz) ∈S, 干涉仪输出相位差为
式 (4) 中d为干涉仪基线的长度, λ为目标辐射源的波长, n为模糊数, 并且
(Tx, Ty, Tz) 是S内点的空间直角坐标, 大地坐标 (L, B, H) 与空间直角坐标 (Tx, Ty, Tz) 的转换关系为
式 (6) 中, N是当地卯酉圈曲率半径, a为椭球模型的半长轴, e2为第一偏心率平方。对于地面和低空目标来讲, H→0, 式 (6) 简化为
由式 (5) 可知, 在 (Tx, Ty, Tz) 处干涉仪输出相位差的值受x、y、z、mx、my、mz、λ的影响是独立的, 那么
目标在目标估计区域S内概率相等, 即P ( (Tx, Ty, Tz) ) =SΔLΔB ( (Tx, Ty, Tz) ∈S) , 点 (Tx, Ty, Tz) 代表子区域的中心, ΔLΔB代表子区域的面积。若干涉仪输出相位差为φm, 那么后验概率[4]
设干涉仪测量相位差的误差是均值为零的高斯白噪声, 那么
式 (10) 中为真值, σ为测量误差的标准差。因此
将式 (11) 带入式 (9) , 得
式 (12) 的分母积分后是一个常数值, 不失一般性, 设
在忽略噪声的情况下, 设 (T'x, T'y, T'z, ) 为目标位置, 经过N次观测后有
在实际应用中, 测量误差较大时, 式
不便于数据处理。用概率累加的方式代替概率累积, 即
称P' ( (Tx, Ty, Tz) /φN) 为伪后验概率。不难得出, 在与式 (15) 条件相同的情况下
在考虑测量误差的情况下, 由P ( (Tx, Ty, Tz) /φN) 最大值确定的定位点在目标位置附近。多次测量得到一组定位点, 首先确定定位点集中的区域, 然后对该区域内的定位点做加权处理, 便可得到目标的精确位置估计。
1.3 算法定位性能分析
研究无源测向定位时, 还应描述定位误差的分布情况。任一种定位系统对空间不同位置的辐射源的定位精度往往是不同的, 这是由于辐射源处在不同位置时, 即使相同的测量误差在不同位置所交汇的区域是不同的, 对于固定的定位系统, 定位误差是目标位置的函数;而对于移动的定位系统, 定位误差是目标辐射源相对于定位系统的几何关系的函数。通常, 工程上用GDOP (geometrical dilution of precision) 描述目标辐射源位置与定位系统位置 (或运动轨迹) 之间几何关系对定位误差的影响[9]。定义GDOP为
式 (19) 中, Δx、Δy、Δz是定位系统给出的辐射源位置坐标与辐射源真实位置坐标的偏差。从直观上讲, 距离星下点越远, 卫星与目标的连线与地面的夹角越小, 卫星与目标的几何关系越差, 定位误差越大, 可以推想定位误差的分布应为星下点的同心圆, 即定位精度是和目标相对于卫星的位置是有关的。为了简化计算过程, 减少变量的个数, 利用大地坐标代替空间直角坐标。由式 (5) 和式 (7) 可知, 对于位于 (L, B) 的地面上任一辐射源信号到卫星两天线的相位差为
由式 (12) 可知, 以干涉仪测量相位差为定位的基本信息时, 定位性能由相位差的测量准确性确定, 并使得P[ (Tx, Ty, Tz) /φm]取最大值的点为目标位置, 此时
由式 (4) 、式 (5) 和式 (22) 可得
式 (23) 中
分别对L、B求全微分, 分析精度和纬度对定位[10]精度的影响。对L、B分别求全微分, 可得
写成矩阵的形式:
从式 (27) 可以直接看出, 相位差测量值φm越精确、基线波长比 (d/λ) 越大, 算法的定位经度越高, 并且由于E[dφm]=0, 从而E[du]=0, 因此基于最大后验概率的定位算法是无偏的。在大地坐标系中, 定位的协方差矩阵为[11]
算法的GDOP为
3 仿真分析
3.1 仿真场景描述
辐射源参数:频率3 GHz, 位于 (-16.82N, 73.45E, 0 km) 的地面上, 卫星对该辐射源的有效侦测距离为1 300 km。
卫星参数:卫星沿圆轨道运行, 轨道高度为500km, 倾角为45°, 星载旋转干涉仪基线长度0.5 m、转速1 r/s, 数据更新间隔为0.05 s, 相位差测量误差为均方根10°的零均值高斯白噪声。
图2中蓝色粗线表示卫星能够侦测到辐射源信号的轨道。
3.2 算法定位仿真
取经度、纬度间隔为ΔL=ΔB=0.05°, 依据干涉仪输出的相位差, 给出了卫星侦察区域每一点的后验概率, 如图3所示。依据后验概率最大化原理, 确定目标位置为: (-16.8N, 73.45E) , 相比目标真实位置偏北0.02°, 定位仿真结果表明了算法的有效性。利用目标周围伪后验概率值大的特点, 可缩小搜索范围, 细化搜索间隔以提高定位精度。还可通过将多次测量的定位点做加权处理 (比如求其均值) , 来得到更精确的定位精度。
针对侦察区域南纬-17.2°~-16.5°、东经73.3°~73.6°取ΔL=0.01°, 分别采用100、200、500、700、1 000、1 200、1 500、1 800个数据并将其分成20组仿真定位效果, 并将每次仿真得到的20个定位点的均值作为定位最终结果, 计算定位误差和圆概率偏差 (CEP) , 结果见图4。
从图4可以看出, 由于测量误差的不确定性, 导致每次的定位结果波动, 但是从走势仍然可以得出利用数据量越多、定位效果越好的结论, 即采用数据量越大, 定位的误差、CEP越小, 说明定位点更加集中在目标点附近, 定位精度越精确、稳定度越高。但要采集更多的数据, 会不可避免的增加定位时间。从图中可以看出, 采用的数据量超过某一范围后, 定位效果会趋向极限, 在采用1 000组数据后, 已能对目标位置做出精确估计, 耗时50 s。
3.3 算法的定位精度分布仿真
当基线方向与卫星速度方向夹角为100°时, 计算GDOP为图5。
当基线方向与卫星速度方向夹角为340°时, 计算GDOP为图6。
分析与结论:从图中可以看出, 定位的GDOP与卫星位置 (星下点) 和基线方向有重要联系。星下点附近定位误差最小, 基线方向定位误差增加最快, 垂直基线方向定位误差增加最慢。这是由于在基线方向上的点与卫星的连线与地面的夹角最小, 使卫星与目标的几何关系最差, 导致定位误差较大。在每次定位中, 由于基线是连续旋转的, 并且基线旋转一周采集数据次数为20, 时间为1 s, 因此须将20次的定位GDOP加权平均, 才能表现旋转干涉仪定位的GDOP, 如图7所示。
从式 (27) 可以得出, GDOP与相位差测量值φm的误差成正比, 没有过多分析的价值。
当基线波长比为20时, 将20次的定位GDOP加权平均, 得到定位的GDOP为图8。
当提高基线波长比后, 定位精度得到了显著地提高。基线比增大后, 能增加波程差, 弥补卫星与目标连线夹角过小的问题, 改善卫星与目标的几何关系。
保持3.1的仿真场景不变, 当基线长度为2 m时, 探讨此时旋转干涉仪的定位性能, 如图9所示。
从仿真中可以看出, 增大基线波长比能成倍改善定位精度, 当采用500组相位差时定位精度已达到百米级, 用时仅为25 s。在实际应用中, 增大基线波长的代价不大, 因此该算法具有很强的操作性。
4 结论
本文简要介绍了基于最大后验概率的单星无源定位算法, 着重分析了算法定位性能。理论分析和仿真计算表明了该算法具有定位精度高、速度快的优点。依据式 (12) 可以看出, 该算法具有开放的架构, 还可以研究波长等定位所需的测量量的误差带来的定位误差。理论分析和仿真计算表明, 若能增大基线长度, 提高基线波长比, 能成倍改善定位精度;若提高数据更新速度, 还可进一步减少定位时间。
参考文献
[1] 龙宁.单星无源定位原理及精度分析.电讯技术, 2012;52 (6) :17—20
[2] 曲志昱, 司锡才.基于虚拟基线的宽带被动导引头测向方法.弹箭与制导学报, 2007;27 (4) :92—95
[3] 司伟建, 程伟.旋转干涉仪解模糊方法研究及实现.弹箭与制导学报, 2010;30 (3) :199—202
[4] 姜勤波, 杨利锋, 马红光.机载单站多目标无源定位算法.系统工程与电子技术, 2006;28 (7) :946—948
[5] 刘学, 焦淑红.拟蒙特卡罗自适应粒子滤波的机载无源定位算法.西安交通大学学报, 2011;45 (9) :34—39
[6] 胡来招.无源定位.北京:国防工业出版社, 2004:31—79
[7] 刘鲁涛, 司锡才.开环旋转相位干涉仪DOA算法分析.解放军理工大学学报 (自然科学版) , 2011;12 (5) :419—424
[8] 龚文斌, 谢恺, 冯道旺, 等.星载无源定位系统测向定位方法及精度分析.长沙电力学院学报 (自然科学版) , 2004;19 (2) :64—67
[9] Kolanek J, Carlsen E.Precision geolocation system and method using a long baseline interferometer antenna system.USA, US 7286085B2, 2007
[10] 许耀伟一种快速高精度无源定位方法的研究.长沙:国防科技大学, 1998:4—18
最大概率 篇4
随着风力发电逐步向大规模、高集中开发的方向发展，电力系统必然会面临大量可再生能源的接入问题。风力发电的规模性接入，对电力系统运行的各方面都会带来明显影响，电压稳定研究作为系统安全稳定运行的重要部分，将面临更多的不确定因素。传统电压稳定分析都是基于确定性模型，忽略了系统中的随机注入功率的不确定性，以系统最危险工作模式为研究对象进行电压稳定性评估，评估结果过于保守，灵活性不强。为了弥补确定性方法的不足，在电力系统不确定问题研究中引入概率分析方法，这些方法可分为以下3类：蒙特卡洛MC(Monte Carlo）法、解析法以及点估计法。其中，蒙特卡洛法可以得到高精度的随机变量统计信息，但由于计算效率低下，一般仅用作其他概率分析方法的验证标准；解析法通过卷积技术或输出变量的累积量求解概率密度函数，这种方法假定输入和输出随机变量成线性关系，不能反映系统运行的实际情况；点估计法根据已知变量的概率分布求解未知变量的各阶矩信息，计算量小、精度高，是一种理想的概率分析方法。上述3类方法中，只有蒙特卡洛法可直接得出随机变量概率分布，后面2类方法还需要结合其他方法得出随机变量的概率分布。
现有文献主要集中于概率潮流计算问题的研究[1,2,3,4,5]，关于电压稳定问题概率分析的文献较少[6,7,8,9]：文献[6]将点估计法应用于电压稳定概率分析的计算，但并未计及风电功率的影响；文献[7]研究了考虑分布式电源的静态电压稳定概率问题，研究过程忽略了输入变量的相关性，用级数展开求取概率分布时，存在一定误差；文献[8]采用随机响应面法对电压稳定概率问题进行求解，取得了较为理想的分析结果；文献[9]对含风电系统进行了电压稳定概率分析，文中采用Nataf变换处理相关输入变量，蒙特卡洛方法求解概率问题，求解过程较为复杂。
本文提出了一种基于多项式正态变换PNT(Polynomial Normal Transformation）和最大熵估计的含风电系统电压稳定概率分析方法。利用多项式正态变换方法处理输入变量之间的相关性，结合点估计法将电压稳定概率分析问题转化为确定性负荷裕度非线性优化模型的求解，并计算负荷裕度的统计特征，采用最大熵估计方法估计输出变量概率分布函数。文中还分析了风速参数变化、相关系数矩阵变化、风电接入容量变化以及风电场功率因数变化对负荷裕度统计特征的影响。本文方法解决了点估计法无法处理输入变量相关性和得出输出变量概率分布的不足，保持了点估计法计算简便、精度高的优点。算例分析表明，本文所提方法具有令人满意的计算精度，且计算效率高，变换方法简洁灵活，能够得到全面反映系统电压稳定性的概率信息。
1 多项式正态变换
概率问题中包含不同概率分布的随机变量，有些变量之间还具有相关性，分布形式的不同和变量的相关性都会对分析结果产生影响。因此有必要对输入变量进行预处理，以减小这些因素带来的影响。多项式正态变换可以同时考虑多个非正态相关的随机变量，利用统一的变换步骤将其变换到独立正态空间，变换方法简洁有效，具有良好的通用性。
1.1 单变量正态变换原理
多项式正态变换是利用标准正态分布随机变量z的多项式对非正态分布的随机变量x进行求解。随机变量x的r阶多项式正态变换表达式为：
其中，a0、a1、a2、…、ar为多项式变换系数。
多项式阶数r可以取不同值，得到不同的多项式变换公式，其中最常用的有二阶、三阶和五阶多项式。尽管高阶多项式可利用高阶矩信息提高多项式正态变换的近似度，但高阶多项式对近似度的改善程度还不明确，而计算过程却更加复杂[10]，因此阶数的确定需要结合具体计算要求进行选择。除了阶数的确定，正态变换的另一关键在于变换系数ai的求取。ai的求解算法主要包括积矩法、L-矩法、最小二乘法以及Fisher-Cornish级数展开法，文献[11]对这些算法做了详细的数值分析和对比，可作为求解算法的选择依据。结合已有文献结论[10,11]，本文选择三阶多项式作为正态变换多项式，并采用L-矩法对变换系数进行求解。
随机变量x的q阶L-矩λq(q=1,2,3,4）与多项式变换系数ai(i=0,1,2,3）满足如式（2）所示函数关系[11]。
利用前4阶L-矩得出变换系数，即可将随机变量x转换为标准正态变量z进行处理。
本文采用Weibull分布对风速进行刻画，可直接利用Weibull分布参数对前4阶L-矩进行计算[12]，如式（3）所示。
其中，τ3=λ3/λ2和τ4=λ4/λ2分别为L-偏度和L-峰度；kw和c分别为Weibull分布的形状参数和尺度参数。通过式（3）即可完成服从Weibull分布的随机变量多项式正态变换。
1.2 多变量多项式正态变换
多随机变量概率分布模型的建立需同时考虑各变量的边际分布和变量间的相关性因素。由1.1节可知，各边际分布随机变量均可采用单变量正态变换方法进行变换处理，本节解决的是把原始随机变量的相关系数矩阵变换为对应的标准正态变量的相关系数矩阵，实现多变量的多项式正态变换。
假设X=[x1,x2，…，xn]T表示由n个随机变量xi组成的随机向量，对应的相关系数矩阵为,Z=[z1,z2，…，zn]T表示由n个标准正态变量zi组成的随机向量，对应的相关系数矩阵为，则两相关系数矩阵对应元素满足如下关系[11]:
其中，ami和amj(m=0,1,2,3）分别为xi和xj的多项式变换系数；分别为xi和xj的标准差和期望。依次求解RX中各元素对应的等式方程，并选择满足不等式条件的解作为的值，即可得到标准正态变量相关系数矩阵RZ。
结合1.1节和1.2节，可以得到由独立标准正态随机变量S=[s1,s2，…，sn]T求取原始随机变量的完整公式：
其中，Am=[am1,am2，…，am n]T(m=0,1,2,3）为由1.1节得出的n个随机变量xi对应的多项式变换系数；L为标准正态变量相关系数矩阵RZ经过Cholesky分解后所得的下三角矩阵。
通过多项式正态变换，为点估计法处理含相关性的随机变量提供数据基础。
2 点估计法
点估计法是求解非线性函数Y=h(X）在随机变量X=[x1,x2，…，xn]T作用下的统计分布的一种有效方法，具有良好的计算精度和较小的计算复杂度。其原理是首先求解一个r×n的采样点集F1(X）和每个采样点对应的权重系数，然后利用式（6）计算函数变量Y的前2r-1矩信息：
其中，（μ1，μ2，…，xi,k，…，μn）和ωi,k分别为随机变量采样点和对应的权重系数，xi,k为X的某一分量估计值，μl(l≠i）为其他分量的自身期望；r和n分别为随机变量xi的估计点数和变量个数。研究表明当r>3时，采样点集F1(X）和权重系数集ω可能包含复数解，不符合实际情况，因此r的数值通常取为3。
将随机变量xi的一个估计点选为自身期望μi的三点估计法称为2n+1估计方案。该方案的采样点集F1(X）中有n个相同采样点（μ1，μ2，…，μi，…，μn），只需进行一次计算即可完成对这n个采样点的处理，具有很高的计算效率。2n+1估计方案下的位置系数ξi,k和对应权重系数ωi,k计算公式如下[13]:
其中，i=1,2，…，n；λi,3和λi,4分别为随机变量xi的偏度系数和峰度系数；σi和μi分别为xi的标准差和期望。
点估计法适用于输入随机变量相互独立的情况，若输入变量之间存在相关性，可利用第1节的方法进行处理。
3 基于最大熵的概率分布估计
任意随机变量x的熵定义为：
其中，H(x）和f(x）分别为随机变量x的熵和概率密度函数。
利用最大熵原理POME(Principle Of Maximum Entropy）估计随机变量的概率密度函数的基本思路为：给定随机变量x的相关统计信息，建立约束条件下的最大熵优化模型，在候选概率分布解集Φ[f(x）]中，选择使H(x）最大的f(x）作为随机变量x的概率密度函数最优估计。POME数学模型如下[14,15]:
其中，为N+1个已知函数，称为基函数集；μ0M=1、μnM(n=1,2，…，N）为给定的N+1个随机变量x的相关统计信息。POME模型的解可用的函数进行表达：
其中，λn(n=0,1，…，N）为拉格朗日参数。求解模型式（8)—(10）得出拉格朗日参数，即可获得概率密度函数f(x）的最优估计。
本文选择幂函数形式(n=0,1，…，N）作为POME模型的基函数集，则POME模型变为：
其中，μnM(n=0,1，…，N）为随机变量x的各阶几何矩。基于随机变量的各阶矩信息，即可利用该模型对概率密度函数f(x）进行最优估计。
4 含风电系统的电压稳定概率分析方法
4.1 风电功率模型
风速的概率分布可用双参数Weibull分布进行描述。风电场风速vi的概率密度函数如式（13）所示：
其中，kw和c分别为Weibull分布的形状参数和尺度参数。
风电机组输出功率与风速的函数关系为：
其中，vci、vco和vR分别为切入风速、切出风速和额定风速；PR为风电机组额定功率。
4.2 负荷方向随机模型
电压稳定分析中，通常选择各节点负荷作为负荷增长方向，以保证各节点负荷功率因数不变且按同比例增加。由于系统测量、估计等方面的误差，实际的节点负荷并不是常数，存在一定的随机性，可选择正态分布反映其随机特点。
假定节点i的有功负荷增长方向IP L i满足以基态有功PL i为均值、以σLi为标准差的正态分布，则其概率密度函数如下[8]:
对应的无功负荷增长方向IQLi由式（16）表示为：
其中，QL i为节点i的基态无功功率。
4.3 负荷裕度非线性优化模型
负荷裕度为运行人员提供了系统从当前运行点到电压崩溃点距离的直观量度，是评估系统电压稳定性的最有效指标。选择负荷裕度为目标函数，同时计及系统等式和不等式约束，构建计算最大负荷裕度的非线性优化模型如下[16,17,18,19]:
其中，ji表示节点j与节点i直接相连，包括j=i的情况；SB为系统节点集合；SG为发电机节点集合；SR为无功源集合；SL为支路集合；PGi和QRi分别为节点i的电源发出的有功和无功功率；PWi和QWi分别为节点i的风电输出有功、无功功率；PL i和QL i分别为节点i的有功负荷和无功负荷；ρ为负荷增长系数，即负荷裕度；αi j=δi-δj-φi j;Ui和δi分别为节点i的电压幅值和相角；Yij和φij分别为节点导纳矩阵对应位置的导纳元素幅值和相角；Sij为节点i与节点j之间的支路潮流大小；下划线变量和上划线变量分别对应各自变量的下限值和上限值。式中等式约束由含负荷裕度参数的扩展潮流方程组成，不等式约束由系统静态安全约束组成。
4.4 含风电系统电压稳定概率分析步骤
在负荷裕度优化模型中加入风电输出功率和节点负荷的随机特性，并采用概率方法进行求解，即可实现含风电系统电压稳定的概率分析评估。
本文假定风电场风速和节点负荷的概率分布参数已知，并给出相应的风速相关系数矩阵，则含风电系统的电压稳定概率分析具体步骤如下。
a.确定输入随机变量个数n，在独立标准正态空间S中对n个标准正态变量进行采样，并求解方程（7），形成基本采样点集F1(S）和对应权重系数集ω。
b.将Weibull分布参数和风速相关系数矩阵代入式（2)—(4），求得相应的变换矩阵L和多项式变换系数Am=[am1,am2，…，amn]T(m=0,1,2,3）。
c.将基本采样点集F1(S）代入式（5），得到输入随机变量的采样点集F1(X）。
d.将F1(X）代入负荷裕度非线性优化模型进行求解，得出对应的负荷裕度离散点集F2（ρ）。非线性优化模型采用预测-校正原对偶内点法进行求解。
e.将ω和F2（ρ）代入式（6），计算负荷裕度的各阶矩信息。
f.将负荷裕度各阶矩信息代入模型式（12）进行求解，即可得到负荷裕度概率分布的最优估计。
5 算例分析
5.1 系统概况
本文以IEEE 30和IEEE 118节点系统为基础，对含风电场的电力系统电压稳定进行概率分析，程序实现平台为MATLAB R2009a，所用计算机的CPU主频为2.1 GHz，内存为2 GB。风电场按恒功率因数方式运行，各风电场原始相关参数如表1所示，表中接入台数列“/”前后数字分别表示118节点和30节点的标准接入台数，风电场原始相关系数矩阵RW为：
IEEE测试系统数据参见文献[20],IEEE 30节点系统中，风电场分别接入节点9、25、28，系统总负荷为189.2 MW;IEEE 118节点系统中，风电场分别接入节点23、39、114、117，系统总负荷为3668 MW。
负荷增长方式为全网负荷同时增加，各节点负荷按基态功率因数等比例增长。各负荷分量服从以基态负荷为均值、标准差为5%的正态分布。
5.2 基于最大熵的概率分布估计结果
应用基于最大熵的概率分布估计方法对算例系统进行含风电系统的电压稳定裕度概率计算，得出相应的分布曲线，并同时与40 000次蒙特卡洛模拟结果进行对比分析。
图1、2和图3、4分别给出了IEEE 30和118节点系统在最大熵估计和蒙特卡洛模拟2种方法下对应的负荷裕度概率密度分布曲线和累计概率分布曲线，图中负荷裕度为标幺值，后同。从图中可以看出，由最大熵估计方法得出的分布曲线较为合理地拟合了蒙特卡洛仿真结果。对比图1和图3中由最大熵估计方法得出的概率密度曲线可知，图1中的密度曲线发生了畸变，存在“翘尾现象”，而图3中的密度曲线变形较小，仅在其尾部有一定程度的收缩。结合5.1节的算例数据可知，IEEE 30节点系统的风电接入比例大于IEEE 118节点系统，这表明非正态变量比例的增加会对最大熵概率分布估计结果产生影响。观察图2和图4可知，2种方法所得累计概率分布曲线较为接近，较好地反映了负荷裕度的累积分布特性。
表2给出了服从Weibull分布的风电场风速随机变量多项式正态变换系数，相应的变换矩阵L为：
表3给出了基于最大熵的概率分布估计参数，可用对应的估计函数描述负荷裕度的概率密度分布。表4给出了2种方法下不同算例系统的负荷裕度均值（标幺值，后同）、标准差（标幺值，后同）以及2种方法计算结果的相对误差。由表4数据可得，均值的相对误差都小于1%，标准差的相对误差都小于4%，计算精度令人满意，在工程应用的误差要求范围之内。表5给出了2种方法下的计算仿真规模和计算时间。由表5可得，最大熵估计方法的计算时间仅为蒙特卡洛方法的0.14%和0.6%，大幅节省了计算时间且能获得较为满意的计算结果，为电压稳定概率分析方法的工程实际应用实现提供了可能性。
5.3 风电参数影响分析
为了研究风速概率密度函数变化和风速相关系数矩阵变化对电压稳定裕度的影响，本节假设在不同风速分布参数和相关系数矩阵组合下对电压稳定问题进行概率求解和结果分析，得出风电参数变化对电压稳定裕度概率分布所带来的影响。
本节以IEEE 118节点系统为研究对象，分别以表1的风电场原始风速参数和原始相关系数矩阵RW为标准，同方向增加或减小风速分布参数数值和相关系数矩阵中的元素数值，分别形成以下3个新风速参数场景和3个新相关系数矩阵场景：风速参数场景取名为“较小风速场景”、“较大风速场景1”和“较大风速场景2”，各场景的具体参数如表6所示；相关系数矩阵场景取名为“零相关系数矩阵”、“较小相关系数矩阵”和“较大相关系数矩阵”，其中“零相关系数矩阵”用4阶单位矩阵I4×4描述，“较小相关系数矩阵”和“较大相关系数矩阵”分别用RW S和RW L描述，各系数矩阵参数如下：
图5为不同风速参数场景下计算得出的负荷裕度概率密度曲线。观察图中曲线，可以看出风速参数对概率密度分布的影响：随着风速参数的逐渐增大，概率密度曲线朝着扁平化的趋势发展，曲线顶点对应的概率密度数值逐渐减小，负荷裕度数值逐渐增大。这说明风速参数的增加有助于系统负荷裕度的提高，但随着风速参数的增加，负荷裕度随机变量的变化也逐渐变大，风速变化所带来的波动性更加明显。
表7所示为不同风速参数和风速相关系数矩阵组合下的负荷裕度计算结果，通过对比表中数据，验证了图5所得相关结论的正确性。由上述分析可知，在选择风能资源时，并不能只考虑风资源本身的变化因素，还要考虑所接入系统的电压稳定承受能力等影响因素，这样才能将风能资源和所接系统进行合理的匹配。
表7中还给出风速相关系数矩阵变化时的负荷裕度模型计算结果。观察表中数据可知，相关系数变化会对负荷裕度计算结果产生影响：随着风速相关性增加，有些负荷裕度均值逐渐减小，有些负荷裕度均值先减后增，相关系数的变化引起负荷裕度均值的变化；随着相关系数增大，负荷裕度标准差逐渐增大，表明其波动更加明显。因此在负荷裕度概率分析时，不能忽略风速相关性，必须加以考虑，否则分析结果与实际情况偏差较大，不能正确反映负荷裕度变化的概率特性。本文方法所得结果和蒙特卡洛模拟所得结果的相对误差如表7所示。从表中可以看出，采用本文方法所得的负荷裕度结果与蒙特卡洛模拟所得结果较为接近，结果均值误差较小，标准差误差偏大，但均在可接受范围之内，可认为本文方法所得结果正确有效。综上可知，采用基于多项式正态变换的点估计法能较好地处理风速相关性问题，计算结果反映了负荷裕度变化趋势，相对蒙特卡洛模拟方法，计算时间缩短为原来的1%以下，适合于对电力系统概率问题的分析求解。
5.4 风电场接入容量影响分析
随着风电场的发展，风电容量的接入比例会越来越大，本节就不同的风电接入比例对负荷裕度概率模型进行求解，观察容量变化带来的影响。标准风电接入容量为230 MW，考虑不同的标准接入容量倍数，得到不同风电接入比例下的负荷裕度均值和标准差，具体数据如表8所示。可以看出，风电比例越高，负荷裕度的均值和标准差均有所增加，这是因为一方面风电的接入提供了电源支持，使得负荷裕度水平有所提高，但是由于风电固有的随机特性，容量的增加也让波动更加明显，标准差随之增大。因此，如何有效地减小风电波动性带来的影响是大规模利用风电必须考虑的问题，否则大量的风电接入会给系统运行的稳定性造成不利的影响。
5.5 风电场功率因数影响分析
风电场发电控制通常采用恒功率策略，不同的功率因数设定值对负荷裕度有所影响。本节采用最大熵估计法快速求解不同功率因数下的负荷裕度模型，计算结果如表9所示。表中负功率因数表示风电场需要消耗无功。
由表9可以看出，随着功率因数的增加，负荷裕度均值逐渐增加，标准差逐渐减小。这是因为随着功率因数的增加，整个风电场的无功消耗逐渐降低，无功输出逐渐增加。由于风电场能够向系统提供无功，因此系统的电压稳定情况得以改善，电压稳定裕度有所增加，能够承受更大的系统变化，在同样的有功波动情况下，能使电压稳定裕度的均值更大，而标准差更小，电压稳定特性更加优良。表中还给出了不同相关系数矩阵条件下，功率因数变化时的负荷裕度计算结果，不同条件下的计算结果变化趋势一致，验证了前述的原因分析。
6 结论