最大估计算法

最大估计算法（精选7篇）

最大估计算法 篇1

0 引 言

第三代合作伙伴计划(Third Generation Partnership Program,3GPP)于2004年年底开始了通用移动通信(Universal Mobile Telecommunication System,UMTS)技术的长期演进(Long Term Evolution,LTE)项目,目的是提高UMTS通信标准以跟进未来移动通信技术要求。在UMTS Release 8的空中接口标准中,单载波频分多址接入(SC-FDMA)和正交频分多址接入(OFDMA)分别作为物理层上行和下行的关键技术之一被采用[1,2]。OFDM是一种数字多载波调制技术,其原理是将数据串并变换后调制到一组等小间隔的正交子载波上,这一过程可通过离散傅里叶变换(IFFT)实现。OFDMA尤其适用于频率选择性信道和高数据速率传输。

目前,针对OFDM系统的信道估计方法已得到广泛研究[3]。其中,基于导频的传统LS信道估计方法因其低复杂度且不需要任何信道的先验信息,实现简单而被广泛应用。但LS算法的性能较差,受噪声影响较大[4]。为了提高估计精度,文献[5]提出了一种降低时域信道冲击响应采样率的方法,以减少噪声干扰。本文基于LS算法和降采样思想,为3GPP LTE下行通信提出了一种改进的信道估计方法。首先,对接收到的时域信号做一次LS估计,选择其中能量较集中的几个径抽头,而剩余能量低的径主要是噪声成分,将这些信道抽头置0,就可在一定程度上抑制噪声干扰;然后将处理后得到的时域信道矩阵变换至频域,采用低复杂度的频域均衡来检测信号。

1 系统描述

在LTE物理层规范中,除非特殊说明,各种域的时域大小表示为时间单位Ts=1/(1 500×2 048)的倍数,一个无线子帧的长度可以表示为Tf=307 200×Ts=10 ms。LTE支持两种类型的无线帧结构:类型1,适用于FDD模式;类型2,适用于TDD模式[6]。本文是基于类型1进行的仿真研究。

图1给出了LTE FDD模式下的OFDM调制原理框图。映射到星座图后的串行输入复信号ai(i=0,1,…,N-1) 经过串/并变换后调制到N个正交子载波上,各子载波间间隔为Δf,其后对这个N点的符号序列进行IFFT变换来实现OFDM调制,这样就得到一个时域的OFDM符号。为了抑制OFDM符号间干扰(ISI)及码间干扰(ICI),信号传送前须将每个OFDM符号的后面NCP点的数据复制并添加到该OFDM符号的前面,这一部分被称为循环前缀(CP)。循环前缀的长度选取要大于多径信道的最大时延,这样就可有效抑制因多径时延造成的ISI。添加了循环前缀的传输序列s(k)经过无线信道,到达接收端,经过采样频率为1/Ts的下采样后得到序列r(k),接着依次进行去循环前缀,OFDM解调(FFT变换),信道估计和信号检测。表1给出了系统的调制参数。

图2为LTE时频域帧结构示意图。当采用常规CP时,一个时隙(0.5 ms)符号包含7个OFDM符号。由表1可知,每个OFDM符号中,N个子载波中只有中间Nm个子载波用来调制有效数据,剩余对应于中间直流数据以及两边用做保护带数据的子载波并未得到有效调制。如图2所示,每个LTE子帧中时域上插有两列导频符号:第1参考符号和第2参考符号,它们分别位于每个时隙的第1个OFDM符号和倒数第3个OFDM符号中。在频域,每6个子载波插入一个参考符号,这个数值是在信道估计性能和RS开销之间求取平衡的结果。

在接收端,时域下的接收信号可以表示为:

undefined

其中:h是L×1的时域信道冲击响应向量;FL是N×L的傅里叶矩阵,如式(2)所示;A是N×N的对角矩阵,对角上的元素是待发送的频域信号;FH是N×N反傅里叶矩阵;n则是一个N×1向量,代表加性高斯白噪声。

undefined

其中:

undefined

2 算法描述

2.1 最小二乘法信道估计

式(1)中接收到的时域信号还可以表示为:

undefined

其中:

undefined

式中:Ad和Ap为两个N×N对角矩阵,其对角元素分别是映射到星座图的待传输复信号和导频信号。假设待传输的复信号为统一分布于16QAM星座图上的复信号,那么可以做式(7)所示的假设:

undefined

用LS准则估计出信道时域冲击响应向量h:

undefined

将式(7)代入式(8)可得:

undefined

2.2 最大径算选择算法

将式(9)得到的时域信道响应向量undefined经过FFT变换到频域,就可以得到最终用于频域信号检测的信道响应。显然,这样得到的频域信道响应的性能好坏取决于时域的LS估计。由于大部分时域信道冲击响应径抽头都会受到高斯白噪声的影响,那么考虑将一些功率弱的径置0,而只选取功率最大的J条径,就可以起到抑制高斯白噪声的作用,因为被置0的一些弱径上的功率成分很大程度来自噪声干扰。本文就是基于这一思想提出改进算法。图3给出了运用最大径选择算法的系统框图。对于径数J的确定要视仿真环境而定。这里选取3GPP空时信道多径数的两倍进行仿真[7,8,9,10]。

用n0,n1,…,nJ-1表示J条被选取的功率最大径,则有:

undefined

再进行式(11)的运算就可以得到频域的信道响应向量,并用于最后的频域信号检测。

undefined

3 仿 真

本文利用20 MHz带宽进行仿真,表2和表3列出了仿真参数。

对基于LS的降采样信道估计算法和改进算法同时进行仿真以做比较。图4和图5分别为不同车速下两种算法仿真结果。结果显示针对LTE OFDM系统提出的改进的信道估计算法在一定程度上提高了估计性能,降低了误码率,且在高速环境下性能提高更明显。同时,由于该算法将某些径置0使得后面的傅里叶矩阵相乘运算更简单,从而降低了运算复杂度。

4 结 语

基于对LTE OFDM系统标准的研究,本文在传统LS估计方法的基础上提出了基于最大径选择的改进算法。仿真结果表明基于最大径选择的LS信道估计方法通过进一步消除噪声干扰,在保证低运算复杂度的前提下提搞了估计性能,降低了系统误码率。

参考文献

[1]3GPP,TR 25.814.Physical Layer Aspects for EvolvedUTRA[EB/OL].Available:www.3gpp.org.

[2]沈嘉,索士强,全海洋,等.3GPP长期演进(LTE)技术原理与系统设计[M].北京:人民邮电出版社,2007.

[3]J.J.van de Beek.On Channel Estimation in OFDM Systems[A].Vehicular Technology Conference[C].1995,2:815-819.

[4]Sinem Coleri,Mustafa Ergen,Anuj Puri,et al.Channel Esti-mation Techniques Based on Pilot Arrangement in OFDMSystems[J].IEEE Trans.on Broadcasting,2002,48(3):223-229.

[5]Andrea Ancora,Calogero Bona,Dirk T M.Slock.Down-sam-ple Impulse Response Least-squares Channel Estimation forLTE OFDMA[A].IEEE ICASSP[C].2007.

[6]3GPP TS 36.211 V8.4.0.Physical Channels and Modula-tion.2008.

[7]3GPP TR 25.996(V8.0.0).Spatial Channel Model forMultiple Input Multiple Output(MIMO)Simulation.2003.

[8]I I Gyu Kim,Youngnam Han,Young Hoon Kim,et al.Transmit Diversity and Multiplexing Methods for 3G-LTEDownlink Control Cha-nnels[A].Vehicular TechnologyConference[C].2006:1-4.

[9]Greenstein L,Erceg V,Yeh Y S,et al.A New Path-Gain/Delay-Spread Propagation Model for Digital Cellular Cha-nnels[J].IEEE Trans.on Vehicular Technology,1997,46(2):477-485.

[10]胡宏林,徐景.3GPP LTE无线链路关键技术[M].北京:电子工业出版社,2008.

最大估计算法 篇2

在许多科学和工程领域(如雷达,通讯,声纳等),线性调频信号(简称LFM,又称chirp)得到了广泛的应用,对它的参数检测和估计是信号处理领域中十分关心的问题。对LFM信号的参数检测和估计可以采用解线调、Radon-Wigner变换和分数阶Fourier变换、最大似然估计(MLE)等方法。其中ML估计是一种渐进无偏的估计,当采样点数很大时,对参数的估计方差接近于Cramer-Rao下界,可以说最大似然估计是一种最有效的估计方法,但由于它需要进行二维搜索,存储量和计算量都很大,因此实时应用受到了限制。本文提出了一种基于时域解线调的最大似然估计的改进算法,提高了估计效率和精度,实现了对LFM信号的快速检测和精确估计。

1 最大似然估计法

最大似然估计(MLE)已被证明是一种渐进最优估计,在有限样本的情况下它具有最优的估计性能[1],因此在信号参数估计中,得到了非常广泛的应用,其基本原理如下:

设接收到的信号为:

x(n)=s(n)+w(n), 0≤n≤N-1, (1)

式中,$s(n)=a\exp [\text{j}2\text{π}(f{}_{0}n+\frac{1}{2}kn{}^2)]$为线性调频信号;a,f0,k为线性调频信号的幅度、初始频率和调频斜率;w(n)为零均值的复高斯白噪声,方差为σ2。

式(1)中f0,k的最大似然估计[2]为:

$f{}_0,k=\underset{f{}_0,k}{{\displaystyle \arg \max }}\left|{\displaystyle \sum _{n=0}^{{\rm N}-1}x}{}_n\text{e}{}^{-\text{j}2\text{π}(f{}_{0}n+\frac{1}{2}kn{}^2)}\right|$。 (2)

在f0、k均未知的情况下,要得到两参数的最大似然估计,不得不采用全域搜索,其计算量很大,仿真实验证明,对含噪声LFM信号直接采用最大似然估计,计算量很大,而且估计精度也要低于本文提出的最大似然估计的改进算法。

2 延时相关解线调

所谓“解线调”,顾名思义就是解除信号的线性调制,使之变为具有单一频率的信号。更进一步的理解是指:在本地产生一个不含噪声的且初始频率为0的LFM信号,与接收到的待估计的LFM信号相乘,以此来消除接收信号中的线性调制项[3]。

所谓“延时相关”,就是信号延迟某一时间τ后与原始信号的关联程度。设信号s(t)延迟时间τ后的表达式为s(t+τ),则定义二者的瞬时自相关函数为:

R(t,τ)=s(t+τ)s*(t), (3)

而“延时相关解线调”就是利用信号的瞬时自相关函数来解除信号的线性调制规律。下面举例说明这一算法的应用。

设接收的LFM信号为:

$s(t)=\exp (\text{j}2\text{π}(f{}_{0}t+\frac{1}{2}kt{}^2))+n(t)‚0\le t\le {\rm T}$。 (4)

则当信号延迟时间τ后,信号的瞬时自相关函数为:

R(t,τ)=s(t+τ)s*(t)=exp(j2π(f0τ+kτt+

kτ2/2))+n′(t), 0≤t≤T, (5)

式中,n′(t)为新的噪声。若延时τ固定,则R(t,τ)又可以被看作频率为kτ、噪声为n′(t)的正弦信号。

从解线调的角度来看,在τ固定时,LFM信号的瞬时自相关函数R(t,τ)还可以看作是2个含噪声的LFM信号相乘,即R(t,τ)=s1(t)s*(t),其中:

s1(t)=s(t+τ)=

exp(j2π(f0(t+τ)+k(t+τ)2/2))+n1(t)=

exp{j2π[(f0+kτ)t+kt2/2]}exp[j2π(f0τ+kτ2/2)]+n1(t)=

exp{j2π(f1t+kt2/2)}exp[j2π(f0τ+kτ2/2)]+

n1(t),f1=f0+kτ。 (6)

比较式(4)与式(6)可知:LFM信号s(t)和信号s1(t)中包含的那一部分LFM信号,二者具有相同的调频斜率,不同的初始频率,但它们的初始频率又具有f1=f0+kτ的对应关系,故很容易想到用信号s(t)去解调信号s1(t),然后用正弦频率估计算法(本文中采用经典功率谱估计法)从R(t,τ)中估计出f1-f0,从而估计出调频斜率$\widehat{k}$,即

$\widehat{k}=\frac{f{}_1-f{}_0}{\tau }$, (7)

这就是延时相关解线调的应用[4]。

3 最大似然估计改进算法

对于接收到的含噪声LFM信号,首先利用延时相关解线调,将其转化为含噪声的正弦信号,该正弦信号的频率为kτ(k为LFM信号的调频斜率,τ为信号延时时间)。采用经典功率谱估计中的韦尔奇法[5],去估计该含噪声的正弦信号频率。韦尔奇(Welch)法对巴特利法进行了2方面的修正:一是选择适当的窗函数w(n),并在周期图计算前直接加进去,这样得到的每一段的周期图为:

$\stackrel{⌢}{{\rm P}}{}_{{\rm P}ER}^i(k)=\frac{1}{{\rm M}U}\left|{\displaystyle \sum _{n=0}^{{\rm M}-1}x}{}^i(n)w(n)W{}_{{\rm M}}^{-kn}\right|{}^2$, (8)

式中,$U=\frac{1}{{\rm M}}{\displaystyle \sum _{n=0}^{{\rm M}-1}w}{}^2(n)$为归一化因子。加窗的优点:一是使得无论对于什么样的窗函数均可以使谱估计为非负值;二是在分段时,可以使各段之间有重叠,这样会使方差减小。

韦尔奇法估计出正弦信号的频率$\widehat{f}{}_1$,进而得到LFM信号调频斜率的粗略估计值$\widehat{k}{}_1=\widehat{f}{}_1/\tau$,将$\widehat{k}{}_1$作为最大似然估计的初始值,对LFM信号的调频斜率进行最大似然估计,得到LFM信号调频斜率的精确估计值$\widehat{k}$。

构造一个调频斜率为$\widehat{k}$、初始频率为0的LFM信号$\exp (\text{j}2\text{π}\widehat{k}t{}^2)$,用该LFM信号去解调原含噪声LFM信号,得到一个频率为原LFM信号初始频率的含噪声正弦信号,再利用韦尔奇法估计该信号频率$\widehat{f}{}_2$,将$\widehat{f}{}_2$作为最大似然估计的初始值,对LFM信号的初始频率进行最大似然估计,得到LFM信号初始频率的精确估计值$\widehat{f}{}_0$。

基于延时相关解线调的最大似然估计改进算法的参数估计步骤如下:

① 根据式(5)对LFM信号进行延时相关解线调;

② 用韦尔奇法估计解线调正弦信号的频率$\widehat{f}{}_1$,进而得到LFM信号调频斜率的粗略估计值$\widehat{k}{}_1=\widehat{f}{}_1/\tau ；$

③ 将$\widehat{k}{}_1$作为确定初始值,进行最大似然估计,得到LFM信号调频斜率的精确估计值$\widehat{k}$;

④ 构造一LFM信号$\exp (\text{j}2\text{π}\widehat{k}t{}^2)$,对原LFM信号解调频;

⑤ 用韦尔奇法估计解调频信号频率$\widehat{f}{}_2$;

⑥ 将$\widehat{f}{}_2$作为确定初始值,进行最大似然估计,得到LFM信号初始频率的精确估计值$\widehat{f}{}_0$。

4 仿真实验

考虑一个淹没在复高斯白噪声中的单分量LFM信号,初始频率为5 Hz,调频斜率为40 Hz,幅度为1,观测数据长度为1 000,采样频率为1 000 Hz。

仿真硬件平台:处理器为Inter(R) Pentium (R)4,2.66 GHz;内存为512 MB。

仿真软件平台:Windows XP Professional;Matlab 2006(a)。

以下分别采用3种方法对LFM信号的参数进行估计:方法1:采用延时相关解线调结合经典功率谱估计法;方法2:传统最大似然估计法,即对含噪声LFM信号直接采用最大似然估计;方法3:采用本文提出的基于延时相关解线调的最大似然估计改进算法。采用每一种方法对信号做了100次的仿真实验,所求估计值为100次的统计平均值。3种估计方法的仿真实验结果对比如表1所示。

通过表1的实验结果可以看出:① 本文提出的基于延时相关解线调的最大似然估计改进算法的估计精度要明显优于其他2种方法;该算法的均方误差比延时相关解线调结合经典功率谱估计法的均方误差要低几个数量级,由于确定初始值的设定,其均方误差较传统最大似然估计法也减小很多;② 由于该算法是在一维空间内对LFM信号两参数分别进行最大似然估计,其运算速度比传统最大似然估计法的二维搜索快很多,因此该方法具有很好的工程应用前景。

5 结束语

本文在研究传统最大似然估计、延时相关解线调及经典功率谱估计的基础上,提出了一种基于延时相关解线调的最大似然估计的改进算法,该算法结合经典功率谱估计的直观性和实时性特点以及最大似然估计高精度的特点,估计精度得到进一步提高,同时减小了运算量,因此具有很好的工程应用前景。

摘要：为实现含噪声LFM信号参数的快速检测和精确估计,提出了一种基于延时相关解线调的最大似然估计改进算法,即首先在时域内进行延时相关解线调,然后对解线调后含噪声信号进行经典功率谱估计,得到调频斜率的粗略估计,将此估计值作为初始值,再进行最大似然估计,得到调频斜率的精确估计值,用此精确估计值对原LFM信号进行解线调,再以同样的思路可以得到LFM信号初始频率的最大似然精确估计值。仿真实验证明了该算法的有效性。

关键词：线性调频,解线调,最大似然估计,谱估计

参考文献

[1]KAY S M.Modern Spectral Estimation:Theory and Application[M].EnglewCliffs,NJ:Prentice-Hall,1988:243-270.

[2]ABATZOGLOU T J.Fast Maximum Likelihood Joint Estimation of Frequency and Frequency Rate[J].IEEE Trans.on AES,1986,22(6):708-715.

[3]杨辉.线性调频雷达信号检测与参数估计[D].第二炮兵工程学院硕士学位论文,2006:51-53.

[4]李英祥,肖先赐.低信噪比下线性调频信号检测和参数估计[J].系统工程与电子技术,2004,24(8):43-45.

最大估计算法 篇3

在OFDM系统中, 符号定时和载波频偏的估计是决定其性能的关键问题之一。利用基于循环前缀最大似然算法估计, 是实现符号定时和载波频偏的估计的一种基本的方法[1], 之后又有很多文献报道了对这种方法的改进[2,3,4]。上述最大似然算法对于频偏的估计是以符号同步为前提, 所以当符号同步估计出现偏差时, 频偏估计受到较大的影响。

在最大似然算法的基础上, 通过引入观察窗和与其对应的频偏估计标志位, 根据频偏估计标志位对观察窗口进行相应的移动, 使得当符号同步估计出现偏差时, 仍然更加准确地估计频率偏差。

1 频偏估计和定时同步算法[5,6]

最大似然算法利用循环前缀重复有用数据的特点按照最大似然准则估计出定时与载波频率偏移。OFDM符号循环前缀的结构如图1所示。

图1中, d为符号的定时同步位置, 即OFDM符号的起始位置, 但对于系统来说这个位置d是未知的。观察2N+L个连续样值$r\left(n\right)$, 其中这些样值包含一个完整的N+L个样值的OFDM符号, N为OFDM符号的长度, L为循环前缀的长度。

定义2个集合

${\rm I}=\left\{d,\cdots ,d+L-1\right\}$和${{\rm I}}^{\prime }=\left\{d+{\rm N},\cdots ,d+{\rm N}+L-1\right\}$,

式中, I为第n个符号的循环前缀所对应的位置序列, 集合I′和集合I中元素对应的样值是对应相同的, 因此存在如下的相关特性:

$E\left\{r\left(k\right)r{}^{*}\left(k+m\right)\right\}=\{\begin{array}{cc}\sigma {}_{\text{s}}^2+\sigma {}_n^2‚& \\ \sigma {}_{\text{s}}^2\exp \left(-\text{j}2\text{π}\epsilon \right)‚& \\ 0‚& \end{array}\begin{array}{cc}m=0；& \\ m={\rm N}；& \\ \text{其}\text{他}；& \end{array}$

式中, $\sigma {}_{\text{s}}^2=E\left[\left|s\left(n\right)\right|{}^2\right]‚\sigma {}_n^2=E\left[\left|\eta \left(n\right)\right|{}^2\right]$, 分别表示有用信号和加性高斯白噪声的能量, d和ε是要估计的符号定时同步位置和载波频率偏差。

由文献[5]可以得到, d和ε的联合最大似然估计为:

$\overline{d}{}_{\text{Μ}\text{L}}=\underset{d}{{\displaystyle \arg \max }}\left[\left|\gamma \left(d\right)-\rho \text{ϕ}\left(d\right)\right|\right]$, (1)

$\overline{\epsilon }{}_{\text{Μ}\text{L}}\left(d\right)=-\frac{1}{2\text{π}}\angle \gamma \left(\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}\right)$, (2)

式 (1) 和式 (2) 中, $\gamma \left(d\right)$等于2个长度为Ng, 相隔N个采样点的连续采样序列的互相关值:

$\begin{array}{l}\gamma \left(d\right)={\displaystyle \sum _{n=d{}_{}}^{d+L-1}r}\left(n\right)r{}^{*}\left(n+{\rm N}\right)‚\\ \text{ϕ}\left(d\right)\text{为}\text{接}\text{收}\text{信}\text{号}\text{的}\text{能}\text{量}\text{值}:\end{array}$

$\text{ϕ}\left(d\right)=\frac{1}{2}{\displaystyle \sum _{n=d{}_{}}^{d+L-1}\left(\left|r\left(n\right)\right|{}^2+\left|r\left(n+{\rm N}\right)\right|{}^2\right)}$。

ρ为相关系数的幅度:

$\begin{array}{l}\rho =\left|\frac{Er\left(n\right)r{}^{*}\left(n+{\rm N}\right)}{\sqrt{E\left|r\left(n\right)\right|{}^{2}E\left|r\left(n+{\rm N}\right)\right|{}^2}}\right|=\\ \\ \left|\frac{\sigma {}_{\text{s}}^2\text{e}{}^{-\text{j}2\text{π}\epsilon }}{\sigma {}_{\text{s}}^2+\sigma {}_n^2}\right|=\frac{\sigma {}_{\text{s}}^2}{\sigma {}_{\text{s}}^2+\sigma {}_n^2}=\frac{S{\rm N}R}{S{\rm N}R+1}‚\end{array}$

式中, SNR为信噪比。

2对最大似然算法的改进

考虑3个连续符号的定时, 根据式 (1) 依次可以得到第n、n-1和n-2的3个符号起始位置估计$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}\left(n\right)$、$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$, 该算法流程如图2所示。

定义F为频偏估计标志位, F=-1表示对于当前符号频偏估计观察窗要前向选取, F=1表示频偏估计观察窗要向后选取。

当$|\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-1)|=|\overline{d}{}_{\text{Μ}\text{L}}(n)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-2)|$时,

当$|\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-1)|\ne |\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-2)|$时,

$|\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-t{}_0)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)|=\min \{|\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-1)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)|,|\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n-2)-\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)|\}$, 式中$t{}_0\in \left\{1,2\right\}$。

采用对$\stackrel{⌢}{\epsilon }{}_{{\rm M}L}$在观察窗内取平均的方法, 设观察窗的窗口长度为Nc+1。

$\stackrel{⌢}{\epsilon }{}_{\text{Μ}\text{L}}=\frac{1}{{\rm N}{}_c}{\displaystyle \sum _{i=\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)+F*{\rm N}{}_c}^{\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)}\stackrel{⌢}{\epsilon }}{}_{\text{Μ}\text{L}}\left(i\right)$。 (5)

频偏估计标志位F的符号取决定着观察窗的方向, 当$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$都比$\widehat{d}{}_{\text{Μ}\text{L}}\left(n\right)$小时, 说明当前符号的估计起始位置可能偏后, 即真实的起始位置在当前估计的起始位置之前, 观察窗的方向要向前, 目的是为了将真实的起始位置包含在观察窗之内。以便在观察窗内进行均值运算时, 尽可能多使真实值附近的频偏估计值参与到均值运算中。

当$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$都比$\widehat{d}{}_{\text{Μ}\text{L}}\left(n\right)$大时, 说明当前符号的估计起始位置可能偏前, 真实的起始位置在当前估计的起始位置之后, 采用向后取观察窗的方式。

以上均为$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$对$\widehat{d}{}_{\text{Μ}\text{L}}\left(n\right)$的作用都是同向的, 下面要讨论的是$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$对$\widehat{d}{}_{\text{Μ}\text{L}}\left(n\right)$的作用异向。

当$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-1\right)$和$\widehat{d}{}_{\text{Μ}\text{L}}\left(n-2\right)$对$\widehat{d}{}_{\text{Μ}\text{L}}\left(n\right)$的作用的作用方向不同, 在前2个符号的估计起始位置中, 以接近当前符号估计的起始位置为判定的依据。

通过和当前符号起始位置的比较来确定观察窗口的方向, 当前符号的估计起始位置在为判定依据的符号估计起始位置之前, 则选择向后选取观察窗;否则, 选择向前选取观察窗。

在对符号进行处理时, 该算法选取的观察窗为3个符号长度, 满足观察窗的长度大于N+2L的要求。缓冲单元中存储3个长度的OFDM符号, 假设当前符号为第n个, 将接收到的第n、n-1和n-2个符号同时进行处理。根据式 (1) 做最大似然估计运算, 在最大似然运算之后得到当前符号的定时信息, 并传输到符号定时单元, 实现对当前符号的定时, 确定其起始位置。根据式 (3) 和式 (4) , 可以得到符号的频偏标志位F, 根据式 (2) , 估计当前符号的频偏, 得到的符号频偏信息传输到频偏补偿单元。

最后, 根据F频偏标志位, 确定观察区间移动的方向, 对频率偏差进行取平均运算。符号估计的起始位置和真实起始位置关系如图3所示, 其中e表示预设的频偏值。

假设当前符号的起始位置$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}\left(n\right)$估计的不准确 (图中用d2表示) , 相对真实起始位置$d^{\prime }\left(n\right)$向后 (图中用d1表示) 。

如图3中 (a) 所示, 当频偏估计标志位F=-1时, 因为$\stackrel{⌢}{\epsilon }{}_{\text{Μ}\text{L}}$估计曲线在真实起始位置$d^{\prime }\left(n\right)$处, 最接近频偏估计真实值, 所以$\left|e-\frac{1}{{\rm N}{}_c}{\displaystyle \sum _{i=\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)-{\rm N}{}_c}^{\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)}\epsilon }\left(i\right)\right|<\left|e-\epsilon \left(\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)\right)\right|$, 其中e为频偏真实值。如图3 (b) 所示, 当$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)$相对真实起始位置d′ (n) 向前时, 频偏估计标志位$F=1‚\left|\epsilon -\frac{1}{{\rm N}{}_c}{\displaystyle \sum _{i=\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n){}_{}}^{\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)+{\rm N}{}_c}\epsilon }(i)\right|<\left|\epsilon -\epsilon \left(\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}(n)\right)\right|$成立。

该文提出对于得到的符号起始位置$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}$附近区域的所有$\stackrel{⌢}{\epsilon }{}_{\text{Μ}\text{L}}\left(i\right)$取平均, 这样即使$\stackrel{⌢}{d}{}_{\text{Μ}\text{L}}$存在这符号同步的偏差, 能够使得$\stackrel{⌢}{\epsilon }{}_{\text{Μ}\text{L}}$更接近于真实的频率偏差ε。该算法与原算法仿真比较如图4所示, 图中仿真环境为, 采用3径信道, 信道增益分别是0.8、0.2和0.05, 对应的延时样值为64、128和192。图4中频偏标志位F=1, 向后取观察区间, 图5为频偏标志位F=-1。 由图4和图5可知, 通过仿真得到, 当符号同步出现误差, 无论在F=1还是在F=-1时, 该算法能够减小由于符号起始位置的估计不准确而造成的频偏估计的误差, 从而更好地对频偏进行估计, 实现载波的同步。

3结束语

研究了最大似然算法中频率偏差估计, 通过对连续3个符号进行处理, 得到频偏估计标志位, 确定观察窗的移动方向, 对最大似然估计中得到的频偏估计值在观察窗内进行取平均运算, 从而即使当符号同步估计出现误差, 仍然更加准确地估计频率偏差。

摘要：深入研究了基于循环前缀 (Cyclic-Prefix, CP) 的最大似然算法 (ML) 算法。针对这个问题, 对基于CP的ML算法做了如下改进, 改进的算法同时对连续3个符号进行ML估计。符号定时同步时, 当前符号同步信息能够通过前面2个符号的同步信息得到校正。在进行频偏估计时, 通过使用观察窗, 在窗口内进行连续频偏值的均值运算, 对频偏估计值误差进行矫正。

关键词：频偏估计,正交频分复用,循环前缀

参考文献

[1]LANDSTROMD, WILSONS K, VANJ J, et al.Symbol Time Offset Estimation in Coherent OFDM Systems[J].IEEE Transaction on Communication, 2002, 50 (4) :545-549.

[2]陆震, 王勇.OFDM系统中的一种符号定时和载波频率同步算法[J].华中科技大学学报 (自然科学版) , 2007, 35 (6) :63-75.

[3]王卫兵, 朱秋萍.频率选择信道下OFDM帧同步及性能分析[J].北京邮电大学学报, 2008, 31 (5) :10:29-38.

[4]RAMASUBRAMANIAN K, BAUM K.An OFDM timing recoveryscheme with inherent delay-spread estimation[C]//Global Telecommunications Conference.IEEE, 5:3111-3115.

[5]汪裕民.OFDM关键技术与应用[M].北京:机械工业出版社, 2006:86-94.

基于最大熵谱估计的水文周期分析 篇4

水文时间序列是水文系统在气象、流域下垫面和人类活动等多种因素共同作用下的输出,其中隐含有随机成分、周期成分和趋势成分。水文学的一个重要研究途径就是利用现有的分析技术对水文时间序列进行描述,以探讨水文系统的演变规律。目前对水文序列主周期的确定,常采用周期图法、方差分析等方法,然而这些方法都存在某些缺陷,如不能处理周期的位相突变、分辨率不高或方差性能不好等[1]。

Burg于1967年提出了最大熵谱估计法,这种方法是以自回归模型为基础的一种参数谱估计。最大熵谱估计克服了经典谱估计方法需要主观假定实测以外数据的缺点,使模型的主观偏见最小。尤其是该方法适用于短序列,但得到的主要周期比较精确,从而广泛应用于通信、化学、经济学等领域。本文将最大熵谱估计法应用于水文周期特性分析,根据岷江紫坪铺水文站1937～2004年68年的年径流实测资料,对紫坪铺站年径流序列的周期变化进行了定量分析和研究,对揭示岷江中上游水文周期的变化规律具有一定的探讨意义。

1 最大熵谱估计的原理及步骤

1.1 最大熵谱估计

熵(Entropy)这个名词来源于统计热力学,用来表示系统无序或混乱程度。Shannon将熵的概念扩展到信息科学领域,定义自信息的数学期望为信息熵[2],即信源{x1,x2,…,xn}的平衡信息量为:

${\rm H}(x)=E[-\lg {\rm P}(x{}_i)]=-{\displaystyle \sum _{i=1}^{{\rm N}}{\rm P}}(x{}_i)\lg {\rm P}(x{}_i)(1)$

式中,P(xi)表示第i个状态发生的概率。

熵的单位是bit,反映了变量的随机性,也是表征随机变量统计特性的一个特征参数,是总体的平均不确定性的度量,随机性越强的序列其熵值也就越大。

熵谱是以熵的概念为基础进行的谱估计,其外推思想是,在观察时间之内的估计值等于观察值;在观察时间之外的取值不做任何假定,即保持最随机、最不确定性,也就是使得熵为最大,从而得到一种新的非线性谱估计法,即最大熵谱法(Maximum Entropy Method)。

设g(λ)是平稳序列的谱密度,$g(\lambda )>0(-\frac{1}{2}\le \lambda \le \frac{1}{2})$,称

${\rm H}={\displaystyle {\int }_{-1/2}^{1/2}\text{l}}\text{g}g(\lambda )\text{d}\lambda (2)$

为此序列的谱熵。而式(2)是定义在$[-\frac{1}{2},\frac{1}{2}]$的连续函数空间上的泛函,以H达到极大为准则来估计其功率谱g(λ),就称为最大熵谱估计。

1.2 Burg算法

工程上常采用Burg提出的递推算法(即Burg算法)来估算序列的参数,其基本思想是按最小二乘法原理估计AR模型中的自回归系数时,使向前预报误差与向后预报误差之和最小[3]。

AR模型可用差分方程表示为:

$x(n)=-{\displaystyle \sum _{r=1}^{p}a}{}_{r}x(n-k)+u(n)(3)$

式中:u(n)为白噪声序列;p为模型的阶数;ar为模型的参数(r=1,2,…,p)。

Burg递推算法具体步骤如下。

(1)计算初始平均功率P0、向前预测误差bp(n)及向后预测误差ep(n),其中

$\begin{array}{l}{\rm P}{}_0=\frac{1}{{\rm N}}{\displaystyle \sum _{n=1}^{{\rm N}}|}x{}_n|{}^2\\ b{}_0(n)=e{}_0(n)=x(n)\\ e{}_r(n)=e{}_{r-1}(n)+k{}_{r}b{}_{r-1}(n-1)\\ b{}_r(n)=b{}_{r-1}(n-1)+k{}_{r}e{}_{r-1}(n)\end{array}$

(2)计算反射系数kr。

$k{}_r=\frac{-2{\displaystyle \sum _{n=r}^{{\rm N}-1}e}{}_{r-1}(n)b{}_{r-1}(n-1)}{{\displaystyle \sum _{n=r}^{{\rm N}-1}|}e{}_{r-1}(n)|{}^2+{\displaystyle \sum _{n=r}^{n-1}|}b{}_{r-1}(n-1)|{}^2}$

(3)根据Levinson-Durbin递推算法求出阶次为r时的AR模型参数。

$\begin{array}{l}a{}_{rk}=a{}_{r-1,k}+k{}_{r}a{}_{r-1,r-k},a{}_{rr}=k{}_r\\ \sigma {}_r^2=(1-|k{}_r|{}^2)\sigma {}_{r-1}^2\\ \sigma {}_0^2=R{}_x(0)=\frac{1}{{\rm N}}{\displaystyle \sum _{n=0}^{{\rm N}-1}|}x(n)|{}^2\end{array}$

(4)将计算出的ark代入下式,直到r等于模型所需阶数,PH(ω)即为所求功率谱密度。

${\rm P}{}_{{\rm H}}(\omega )=\frac{\sigma {}_{r0}^2}{\left|1-{\displaystyle \sum _{r=1}^{r{}_0}a}{}_{r{}_{0}r}\text{e}{}^{-j\omega r}\right|{}^2}(4)$

1.3 模型阶数的确定

AR模型的阶数r一般事先是不知道的,需要在递推的过程中确定。在使用Levison-Durbin递推算法时,可以给出由低阶到高阶的每一组参数,且模型的最小预测误差功率Pmin(相当于白噪声序列的方差σ2)是递减的。从理论上讲,当预测误差功率P达到指定的期望值,或者是不再发生变化时,这时的阶数即是应选的最佳阶数[3]。

因为预测误差功率P是单调下降的,因此,该值降到多少才合适,往往不好选择。为此,有几个不同的准则被提出,其中较为常见的两个如下。

(1)最终预测误差准则(FPE)。

设{x1,x2,…,xn}是均值为零的随机序列,则对于一个阶次为r的预测误差滤波器的最终预测误差FPE(r)定义为:

$F{\rm P}E(r)={\rm P}{}_r\frac{{\rm N}+(r+1)}{{\rm N}-(r+1)}(5)$

(2)信息论准则(AIC)。

信息论准则是建立在作为预测误差滤波器阶次r的函数,预测误差对数似然极小值之上,对于一个阶次为r的预测误差滤波器,定义AIC(r)为:

$A{\rm I}C(r)={\rm N}\ln {\rm P}{}_r+2r(6)$

式中:N为水文时间序列的样本数,当阶数r由1增加时,FPE(r)和AIC(r)都将在某一个r处取得极小值。将此时的r定为最合适的阶数。

在实际运用时发现,当数据较短时,它们给出的阶次偏低,且二者给出的结果基本上是一致的,即:

$\underset{{\rm N}\to \infty }{{\displaystyle \lim }}\{\text{l}\text{g}F{\rm P}E(r)\}=A{\rm I}C(r)$

2 岷江紫坪铺站年径流周期分析

2.1 流域及测站简介

岷江干流全长735 km,流域面积14万km2,属长江上游的一级支流,发源于川、甘两省交界的岷山山脉南麓的贡嘎岭和郎加岭。都江堰以上为上游,以漂木、水力发电为主;都江堰市至乐山段为中游,流经成都平原地区,与沱江水系及众多人工河网一起组成都江堰灌区;乐山以下为下游,以航运为主。岷江有大小支流90余条,上游有黑水河、杂谷脑河;中游有都江堰灌区的黑石河、金马河、江安河、走马河、柏条河、蒲阳河等;下游有青衣江、大渡河、马边河、越溪河等。

紫坪铺水文站位于岷江上游,设立于1936年8月,控制流域面积22 664 km2。根据紫坪铺水文站资料,岷江实测多年平均流量476 m3/s,年径流总量为150.11亿m3,最大月均流量一般出现在7月份,最小月均流量多出现在2月份,最大月平均流量为最小月平均流量的5～7倍,年际变差系数为0.11。由于紫坪铺站期间历经多次停测和迁移,像1953,1988,1989年等年份无实测流量资料。本文由降雨径流关系差补出这些年份的径流资料。

年径流序列的周期分析主要是研究年径流随时间变化而呈现出的周期性变化,本文根据岷江紫坪铺站1937～2004年共68年的年径流资料,采用最大熵谱估计方法进行周期分析。

2.2 资料的标准化处理

在建立自回归模型时先需对原始数据进行标准化处理,即

$Y{}_t=\frac{X{}_t-\overline{X}}{\sigma }(7)$

式中:Xr,Yt分别为t时刻的原始和标准化后的年径流量,$\overline{X}$、σ分别为年径流量的多年平均值和均方差。

2.3 模型阶数的确定及计算结果

本文采用FPE准则确定自回归最佳阶数r0,按式(5)计算得到FPE(r)的关系图,如图1所示,最大试验阶数为30,选取FPE(r)最小值对应的阶数即为最佳阶数。

可见,当r=2时,FPE(r)最小,从而确定最佳阶数r0=2。将r0代入Levinson-Durbin递推算法求出AR模型的参数,用式(4)计算最大熵谱。

在实际计算中,通常采用式(4)的离散形式:

${\rm P}{}_{{\rm H}}(\omega )=\frac{\sigma {}_{r{}_0}^2}{\left[1-{\displaystyle \sum _{r=1}^{r{}_0}a}{}_{r{}_{0}r}\cos (\frac{\text{π}lr}{m})\right]{}^2+\left[{\displaystyle \sum _{r=1}^{r{}_0}a}{}_{r{}_{0}r}\sin (\frac{\text{π}lr}{m})\right]{}^2}(8)$

式中:频率取$\omega {}_l=\frac{2\text{π}l}{2m}(l=0,1,\cdots ,m)‚m$为选取的最大阶数。在序列样本量不大时,m通常取$\frac{n}{2}(n$为序列样本总数)。此时m对应的周期为${\rm T}{}_l=\frac{2m}{l}$。

式(8)的具体推导过程可见参考文献[4]。

计算出最佳阶数r0=2时的熵谱,得到的计算结果见表1,并绘制出最大熵谱谱密度图,如图2所示,该谱图呈单峰型。如果谱密度曲线有尖锐的峰点,其对应的周期就是序列存在的显著周期。由图2可知,在周期为4年处存在最大值,说明在岷江紫坪铺站1937～2004年的年径流时间序列中,存在4年左右的变化准周期。

3 周期的检验

3.1 方法概述

在根据最大熵谱估计提取出水文时间序列的周期后,还需对周期进行检验,因为序列的白噪声可能导致出现虚假的峰值。周期的检验可以采用周期图方法中常用的Fisher检验来判定出现的峰值是否比白噪声的情形要显著得多[7]。

对水文时间序列{xt}进行Fourier级数展开:

$x{}_t={\displaystyle \sum _{i=1}^{s}c}{}_i\cos (2\text{π}f{}_i+\phi {}_i)+\epsilon {}_t(9)$

式中:s=[N/2],N为水文时间序列的样本数,{φi}是在[-π,π]内均匀分布的独立随机变量,{εt}为白噪声序列。

对于水文时间序列{xt},定义IN(f)为:

${\rm I}{}_{{\rm N}}(f)=\frac{1}{{\rm N}}\left|{\displaystyle \sum _{t=1}^{{\rm N}}x}(t)\text{e}{}^{-j2\text{π}ft}\right|{}^2(-\frac{1}{2}\le f\le \frac{1}{2})(10)$

为方便计算,对于实测的水文时间序列,IN(f)的频率一般采用离散形式,记${\rm I}{}_j={\rm I}{}_{{\rm N}}(f{}_j),f{}_j=\frac{j}{{\rm N}},j=0,1,\cdots ,[\frac{{\rm N}}{2}]$。对此,Fisher提出一种max(Ij)检验方法,基于统计量:

$g=\frac{\max {}_{1\le j\le s}({\rm I}{}_j)}{{\displaystyle {\sum }_{j=1}^s{\rm I}}{}_j}(11)$

式(11)被称为Fisher统计量[5]。

对式(9)建立如下零假设:

$\text{Η}{}_0:c{}_i=0,i=1,2,\cdots ,k$

Fisher证明了在H0假设下统计量g的分布为:

$\begin{array}{l}{\rm P}[g>{\rm Z}]=s(1-{\rm Z}){}^{s-1}-\frac{s(s-1)}{2}(1-2{\rm Z}){}^{s-1}\\ +\cdots +(-1){}^a\frac{s!}{a!(s-a)!}(1-a{\rm Z}){}^{s-1}(12)\end{array}$

式中,a为小于$\frac{1}{{\rm Z}}$的最大整数。

根据选择的显著性水平α和Fisher检验表[5]选定Zα,使P(g>Zα)=α。然后根据式(11)计算出统计量g。若计算出的统计量g大于Zα,则表明max(Ij)在显著性水平α时是显著的,即序列的白噪声对周期的影响不显著,水文时间序列{xt}存在隐含周期。

3.2 检验结果分析

根据最大熵谱估计法,计算得出岷江上游紫坪铺站年径流序列存在4年左右的变化准周期,由式(11)Fisher检验计算出统计量$\stackrel{⌢}{g}{}_1=0.237$,取显著性水平 α=0.01,s=[N/2]=34,由Fisher检验表差值得出Z0.01=Z(34,1)=0.208。由于$\stackrel{⌢}{g}{}_1>{\rm Z}{}_{0.01}$,因此接受该序列有一个隐含周期的假设,表明在显著性水平为0.01时,max(Ij)是统计显著的,即该水文时间序列存在4年的显著周期。

4 结 语

(1)本文应用最大熵谱估计法,根据岷江上游紫坪铺水文站的年径流实测资料,进行周期变化规律分析。通过分析和检验,表明岷江上游紫坪铺站年径流序列具有4年左右的变化准周期。

(2)在研究中同时也发现,最大熵谱估计虽然对于提取短序列的周期有一定的实用性,但是也存在峰值漂移的现象,不能准确、客观的反映出水文时间序列在时域和频域方面的变化特征[6],因此其应用范围受到一定的限制。

(3)实际中的水文时间序列,由于受到天文、气象、气候、下垫面和人类活动等众多因素的影响,往往呈现出高度复杂动态的非线性空间过程,因此在今后的研究中还需进一步建立熵与水文过程之间的联系,从而更好的揭示出水文水资源的变化过程和规律。

摘要：最大熵谱估计法是以AR模型为基础的一种参数谱估计方法。该方法通过频谱分析,可以计算出全局性的主要周期,从而反映出水文时间序列的整体特性。介绍了最大熵谱估计法的基本概念和理论,并论述了用Burg算法作为求解AR模型参数的方法。针对岷江紫坪铺水文站19372004年的年径流实测资料,将最大熵谱估计理论应用于水文时间序列周期的提取上,并采用Fisher检验方法对提取出的周期进行了检验。通过分析发现,岷江紫坪铺站年径流序列存在4年左右的变化准周期,对揭示岷江中上游水文周期的变化规律具有一定的探讨意义。在研究中同时也发现,最大熵谱估计也存在一定的缺陷,在今后的研究中还需进一步建立熵与水文过程之间的联系,从而更好地揭示出水文时间序列的变化规律。

关键词：年径流量,最大熵谱法,岷江,水文周期

参考文献

[1]丁晶,邓育仁.随机水文学[M].成都:成都科技大学出版社,1988.

[2]周萌清.信息理论基础[M].北京:北京航空航天大学出版社,1993.

[3]邹鲲,袁俊泉,龚享铱.MATLAB 6.x信号处理[M].北京:清华大学出版社,2002.

[4]魏凤英.现代气候统计诊断与预测技术[M].北京:气象出版社,1999.

[5]顾岚.时间序列分析在经济中的应用[M].北京:中国统计出版社,1994.

[6]易淑珍,王钊.水文时间系列周期分析方法探讨[J].水文,2005,25(4):26-29.

最大估计算法 篇5

准确的线路参数是电力系统中能量管理系统(EMS),如潮流计算、状态估计等各种高级应用的基础。目前EMS中线路参数往往直接采用设计参数,很少使用实测参数。实际上,线路改建、运行环境变化、对补偿电容器的组数掌握不准确[1]等原因,会导致实际参数偏离设计参数,进而影响EMS中的状态估计、不良数据辨识及其他安全分析模块。鉴于此,估计线路实际参数成为电力系统参数辨识的研究方向之一[2]。

传统线路参数辨识方法主要分为两大类:残差灵敏度分析法[3,4]和增广状态估计法[5,6,7]。此外,随着电网中相量测量单元(PMU)布点的增加,将单条线路与全网解耦,直接利用PMU量测获取电网参数的方法也引起了广泛关注[8,9,10,11,12]。其中,残差灵敏度法是在常规状态估计结束后,利用残差与错误参数之间的联系通过多次迭代对错误参数进行修正进而获得参数;该方法不影响已有的状态估计程序,但依赖于现有状态估计方法的准确度。而基于PMU量测的直接计算法利用待辨识线路两端的量测,根据线路模型,直接辨识得到参数,简单易行;然而,分析表明,该类方法中电阻参数的估计结果对电压幅值量测的相对灵敏度极高,若电压幅值有微小误差都会导致电阻参数估计结果的不准确[13];进而意味着该类方法对量测误差要求很高,因此,该类方法目前具有一定的局限性。

增广状态估计法将待辨识参数作为状态量,与原有节点状态量一起进行状态估计。传统的增广状态估计法(如基于加权最小二乘(WLS)的增广状态估计)存在数值稳定性显著变差的问题[5],直接限制了该方法的实用性。针对这一问题,文献[5]提出基于Tabu搜索的参数估计,该方法具有全局寻优能力,并改进了增广矩阵法的数值稳定性,但是Tabu法的搜索方向如何选取不得而知。文献[6,7]用可疑支路潮流补偿量代替支路参数进行参数估计,一定程度上改善了参数估计的数值稳定性。

另一方面,由于量测误差和电网参数误差同时影响状态估计的精度,因此量测误差通常会对参数估计结果产生一定影响,并且量测误差越大,参数估计结果越差。以往的研究大都采用将多个时间断面所采集的数据求取平均值的方式降低随机量测误差的影响[14],或是人为假设量测噪声标准差小于一定范围[15,16](如小于1%),实际上,以上方式并不能从根本上解决量测误差对参数估计结果影响的问题。

文献[17,18,19,20]提出了一种基于测量不确定度理论,以最大测点正常率(MNMR)为目标函数的抗差状态估计新方法,用理论和大量仿真表明了方法的有效性及优异的“抗差”特性。为此,本文拓展MNMR抗差状态估计,提出了基于MNMR的线路参数增广状态估计新方法,以降低量测误差的影响,使参数估计结果更合理稳定。

1 基于加权最小二乘的增广状态估计方法

增广后电力系统状态估计的量测方程可写成:

式中:z为m维量测向量;x为n维状态向量,一般包括节点电压的幅值和相角;p为待估计参数向量;xa为增广后的状态向量,包括x和p;v为量测误差向量,通常假设量测误差服从均值为0、标准差为σ的正态分布;h为用x和p表示的量测函数向量。

相应的基于加权最小二乘的增广状态估计方法

式中:W为权重系数矩阵。

与求解基于加权最小二乘的状态估计类似,求解式(2)的迭代公式为:

式中:Δxa为xa的修正向量;Δz为z的偏差向量;Hq为增广雅可比矩阵。

其中Ha的第i行第j列的元素可表示为:

式中:hi(x,p)为测点i的量测函数;(x,p)j为第j个状态量。

2 MNMR抗差估计和线路参数估计

2.1 参数估计模型

本文方法是在MNMR抗差状态估计方法[17,18,19,20]的基础上提出的,将待辨识参数向量p增广为待估计量,得到相应的参数估计的数学模型为:

式中:zi为测点i的量测值;Ui为与一定置信概率对应的测点i的扩展不确定度;g(x,p)和l(x,p)分别为等式约束和不等式约束;di (x,p)为测点i在状态估计结果x下的相对偏离。

进而,得到测点i的测点评价函数f(di(x,p))为:

式中:K和λ为2个变量,本文K取2,λ取5。

实际上,f(di (x,p))是理想测点评价函数f'(di(x,p)的近似函数:

当|di(x,p)|≤1时,认为测点i为正常测点,反之为异常测点。因此,目标函数f'(di(x,p))最小即意味着测点的正常率最高。

2.2 约束方程

本文选取的状态量x为“全状态”,包括所有节点的电压幅值和相角,以及注入有功功率和无功功率。

对于实际运行电力系统,真实的系统状态必然满足潮流方程约束及发电机出力上限等必须满足的不等式约束。其中等式约束g(x,p)的具体形式如下:

式中:n1为独立节点数;P(i),Q(i),vi分别为节点i的注入有功功率、无功功率和电压幅值;θij=θi-θj,其中θi和θj分别为节点i和j的相角;Gij和B ij分别为节点i-j间的自导纳(i=j)或互导纳(i≠j)的实部和虚部。

不等式约束l(x,p)的具体形式为:

式中:Pkmax和pk,Qkmax和qk分别为节点k的发电机节点的最大、实际有功出力和无功出力。

系统中存在的联络节点为A,零注入功率测点的量测方程也可以写成:

式中

对于输电线路,待估计参数p是电阻R、电抗X和电纳B;根据线路的实际情况和参数估计值的物理意义引入参数变量的有效性上下限约束,如式(13)所示,并且也可以归纳到不等式约束l(x,p)中。

式(13)中,本文将Rmin,Xmin和Bmin设置为0;Rmax,Xmax和Bmax设置为一个较大的数。

3 参数结果的统计特征提取

一般情况下,随机采样过程中量测误差是不同的,而参数辨识结果又容易受量测误差的影响,从而导致参数辨识结果具有一定的波动性。根据多组量测值的随机采样,可得到大量的参数辨识结果,然后进行统计分析,对比其统计特征更具说服性。

本文分别采用核密度估计方法(kernel density estimation,KDE)[21]和点估计方法,获取参数辨识结果的概率密度分布、置信区间、期望及方差等信息,以更准确地反映参数的特征。

核密度估计不需对数据分布附加任何假定,而是直接从数据样本本身出发研究数据的分布特征,本文采用最常用的高斯核函数。通过核密度估计获得辨识结果的概率密度函数后,可通过式(14)获得包含参数真值p的置信度为1-α的置信区间[,]:

点估计包括矩估计和最大似然估计等,由于矩估计不需要对数据分布进行假设,具有无偏性、有效性与一致性,因此本文采用矩估计获得大量辨识结果样本的期望及方差,分别如式(15)和式(16)所示:

式中:n'为样本总数;Xi为第i个辨识结果样本。

4 仿真算例及分析

本节利用基于加权最小二乘的增广状态估计方法(简称为ASE-WLS方法)和本文提出的方法(简称为ASE-MNMR方法),采用4节点系统和IEEE14节点电力系统进行仿真验证。利用ASE-WLS方法计算时,测点i的权重取值为为该测点量测误差标准差。根据文献[20],当置信概率为99.7%时,与之对应的扩展不确定度取。

4节点系统如附录A图A1所示,将ASE-WLS方法和A SE-MNMR方法进行对比分析,以初步表明方法的有效性。量测数据为所有节点的电压幅值和注入功率,以及所有支路的支路功率,4节点系统的量测和参数信息取自文献[1]。在量测真值上叠加2%的高斯白噪声得到试验用量测生数据[22]。

随机采样1组量测数据,如附录A表A1所示,分别用ASE-WLS方法和ASE-MNMR方法估计线路1-2的参数,估计结果如表1所示。

进一步,取100次随机采样,可以得到线路1-2的参数估计结果如图1所示。相应地,基于核密度估计方法得到的参数估计结果的概率密度函数曲线如图2所示,其参数估计结果的统计特征如表2所示。

图1、图2和表2表明,受量测误差的影响,ASE-MNMR方法所得参数估计的置信区间和方差明显小于ASE-WLS方法,即ASE-MNMR方法所得参数估计结果受量测噪声的影响小。

进一步,用IEEE 14节点系统分析不同水平的量测误差及不同位置的量测误差对本文方法的影响,具体结果见附录B。算例表明,随着量测误差加大,相比ASE-WLS方法,ASE-MNMR方法估计结果的波动性变化不大,受量测误差的影响更小,即本文方法具有更好的抵抗量测误差的能力;此外,相较于其他位置的量测量存在的量测误差,待估计线路本身量测的量测误差对本文方法的影响更大些,这与文献[5]的结论基本一致。

5 结语

本文提出了基于MNMR抗差状态估计的线路参数增广状态估计方法,该方法以正常测点率最大为目标函数,有效继承了MNMR状态估计的“抗差”特性,并利用高斯核密度估计和点估计方法提取线路参数估计结果的统计特征。仿真结果表明,与基于加权最小二乘的增广状态估计方法相比,该方法不易受量测误差影响,参数估计结果合理稳定。但是,该方法对一些数值较小的参数(如线路电纳)的估计效果相对较差,如何提高其估计精度尚需研究。

最大估计算法 篇6

传统的基于CP的最大似然估计同步技术是利用CP来确定OFDM符号的到达时刻，并计算接收端与发射端本地振荡器的频率偏移。这种算法的优点是能够有效提高信道的利用率，因为在传输的过程中没有加入训练序列，能够在非分散信道和高斯信道中发挥良好的作用。但是该算法不适用于时变多径信道，如果FFT窗的起始时刻落在CP内，只会引起相位旋转，而这种相位旋转是可以在FFT之后通过均衡器调整过来的;如果FFT窗的起始时刻落在有效数据内就会引起ISI(码间串扰)和ICI(载波串扰)，影响载波之间的正交性。这些缺点在时变多径信道中能够带来致命的伤害，特别是在延时大于CP的时候。本文提出的新型CP结构允许扩大符号定时的误差，起始时刻可以落在CP内，也可以落在有效数据的头部，这种误差也只会带来相位旋转，不会影响系统性能。因此新CP结构能大大提高OFDM在时变多径信道的性能。

1 传统CP最大似然估计的分析

无线信道中的多径时延效应会给OFDM系统带来码间干扰，增大接受端的误码率，为了减少码间干扰需要加入循环前缀，即假设OFDM符号为位，将其后位复制到符号前面，加入循环前缀后的符号长度为M=N+L。

假设传输信道是慢时变信道，同步过程如图1所示，在接收端观测2N+L个连续的r(n)数据样点，则这些样点中含有一个完整的OFDM符号，但这个OFDM符号的位置是未知的，设其起始点位置为d，即为待估计的符号定时同步位置。

定义两个集合:

其中集合I是第i个符号的循环前缀，它与集合I'中的元素相同，将2N+L个观察点作为一个向量。

由循环前缀的特性知道，集合I和集合I'的元素(即r(k),k∈I∪I')是对应相同的，因此运算特性如下:

式中，σs2=E[|s(n)|2]表示传输过程中有用信息符号的能量，σn2=E[|η(n)|2]表示假设是在高斯信道中传输的加性高斯白噪声的能量。

定义关于d和ε对数似然函数如下(d和ε是要估计的符号定时同步位置和载波频率偏差):

式中对同步参数的估计没有影响，可以忽略。对数似然函数由此可以简化为:

式中的分子项f(r(k),r(k+N))为二维复高斯分布的概率密度函数。

由于r(k)为复高斯随机变量，因此上式中的分母项由两个一维复高斯随机概率密度函数组成，其概率密度的表达式可以写为:

利用式(1)中的相关特性，可以得到概率密度表达式为:

式中，ρ定义为r(k)与r(k+N)之间相关系数的绝对值，运算结果如下:

将式(5)和式(6)的概率密度表达式代入式(4)中，经过一些代数运算处理，对数似然函数可以写为:

其中，为常数，Arg(·)表示求复数的辐角，γ(d)和Φ(d)由式(9)给出:

式(9)中的γ(d)表示OFDM符号中连续L个相距为N的数据对的相关系数之和。式(10)中的Φ(d)表示相关部分的能量。

由于c1,c2为常数，且c2>0，故它们对似然函数的最大化没有影响，因此式(8)可以进一步化简为:

其中等号右边的第一项为γ(d)的加权模值，权值由载波频偏决定;第二项是与频偏无关的能量项，由相关系数ρ加权。

由式(11)可以看出需要对符号定时同步和载波频率偏移进行联合估计，因此要找出似然函数的最大值需要分两步进行:

由式(11)可知，为了使似然函数取得最大值，则载波频率偏移ε应满足使式中的余弦项为1，由此可以得到对ε的最大似然估计:

其中n为整数。由于余弦的周期性，使得有多个频偏值满足要求，这会造成频偏估计的模糊。考虑到一般情况下，载波频偏应该在很小的范围内，因此可以取n=0。此时，从而式(11)中的似然函数可以写为:

此时的似然函数只与符号定时同步参数d有关，因此使式(14)最大化就可以获得对d的估计。之后再将估计得到的d值代入式(13)中，即可得到对频偏ε的估计。最后给出对d和ε的联合最大似然估计公式为:

不过该算法不完全适用于时变多径信道，这是因为在时变多径信道环境下，接收端的OFDM符号延时严重，OFDM符号的相关特性遭到破坏，这样会极大的影响系统性能。

2 改进的CP结构

为了能够克服传统CP在时变多径信道的缺点，对CP结构进行如图2所示的改进。

图2中CP个数与传统CP个数一样，只是把它分为个数相等的两部分。观察symboli,CP由两部分组成，CP前半部分是symnoli-1有效数据的头部的复制，CP后半部分是symboli有效数据的尾部的复制。这种结构的CP能够确保每一个OFDM符号的准确性，因为一个OFDM符号有两个扩展部分，前一个符号头部和该信号尾部的复制。

基于这种新型结构的CP估计符号定时同步和载波频率偏移的计算方法与传统CP结构的计算方法一致。但是，现在观察图3集合I由两个小集合组成，I1和I2，集合I1包含的是CP的第二部分，集合I2包含的是该符号有效数据的头部。集合I'也包含两部分，I'1和I'2，集合I'1包含的是该符号的尾部，集合I'2包含的是下一个符号CP的第一部分。结合图2可以看出集合I和集合I'是数据相同的两部分。这种改进的CP结构，不再是计算从CP第一个样点开始的相关函数，而是从CP的中间点开始计算相关性的大小。这样的话，被检测到的OFDM符号序列样本的开始时间就不是CP的第一个样点，而是CP/2,CP/2的位置对应于相关函数的峰值。

使用这种改进的CP结构，被估计出的起始时间允许落入有效数据区间中，即集合I2中。这样的定时误差在FFT之后只是会引起OFDM的相位旋转，而这种相位旋转是可以通过均衡器来纠正的。

3 仿真

3.1 仿真参数

现在使用MATLAB对该算法进行仿真，仿真的条件如下:采用16QAM对信号进行调制，子载波数N=2048,IFFT和FFT的长度NS=4096，循环前缀个数G=(1/32)NS。

3.2 仿真结果

现在OFDM符号序列的长度是NS+G=4224，利用改进的CP接收到的OFDM数据的起始时刻允许落在第65到192样点之间区间的任意位置，而没有码间串扰。

由图4可以看出当传输时延小于80时，改进的CP结构与传统的CP结构计算的误码率几乎一致，但当传输时延大于80，传统CP的误码率就会明显升高，系统性能下降。改进CP的系统性能提高也很明显。

4 结束语

为了克服传统CP在时变多径信道下性能急剧下降的缺点，本文提出对传统的CP结构进行改进。改进的CP，把CP分为个数相同的两部分，第一部分是复制前一OFDM符号的头部，第二部分是复制该OFDM符号的尾部。这样的改进算法能够使估计接收端信号的起始时刻无码间串扰的落在有效数据区间。这些优点能够提高OFDM的系统性能，特别是延迟大于CP/2的情况下。

参考文献

[1]Chen B,Wang H.Maximum Likelihood Estimation of OFDM Carrier frequency O set[C].IEEE International Conference on Communications,Session A03,2002.

[2]Ninh L V,Vu T A.Depression of Multipath ISI in Symbol Timing Synchronization for OFDM[C]//Proceeding of the 9th Biennial Vietnam Conference on Radio&Electronics(REV04),Nov.27-28,Hanoi,Vietnam,2004:49-53.

[3]Visible Light Communications:Tutorials[S].IEEE 802.15 VLC SG,2008.

[4]赵林靖,李健东,陈晨.基于ML估计的高效OFDM整数倍频偏估计算法[J].西安电子科技大学学报:自然科学版,2006,32(4):559-561,565.

[5]王胜飞.可见光通信OFDM系统同步技术研究[D].上海:华东理工大学,2012.

[6]李伟华,蓓蕾,吴伟陵.OFDM系统定时和频率偏移估计[J].北京邮电大学学报,2006,27(1):36-39.

最大估计算法 篇7

传统的状态估计方法由Schweppe等人于1970年创立[1]。根据求解状态变量所用的目标函数的不同,状态估计器有不同的估计准则,常用的估计准则包括:加权最小二乘(WLS)准则、非二次准则[2]、加权最小绝对值(WLAV)准则[3,4]、最小均方(LMS)准则和最小截取二乘(LTS)准则[5]。当前应用最广泛的是WLS估计准则,其优点是模型简单、计算量小,对理想正态分布的量测量,估计结果具有最优一致且无偏等优良统计特性,缺点是抗差能力差,即估计结果易受不良数据的影响而偏离真值较远。非二次准则、WLAV、LMS和LTS估计器都属于抗差估计器,其显著优点是具有抗差能力,但有计算量大、需要主观确定加权因子等缺点,这也是妨碍其发展的重要原因[6]。

纵观以往方法的基本思想,WLS准则是使残差加权平方和最小,WLAV准则是使残差加权绝对值之和最小,以及其他估计方法,都不是使估计结果的合格率最大。然而,实际应用中评价状态估计效果的标准是合格率。传统状态估计方法的理念与状态估计性能评价标准之间存在不可逾越的鸿沟。

本文提出一种新型电力系统状态估计方法,该方法以合格率最大为目标,可用现代内点法求解[7],计算速度快、收敛性好。

1 方法介绍

1.1 主要思想

电力系统的遥测量方程可表示为:

Zi=hi(x)+eii=1,2,…,m (1)

式中:Zi为系统中测点i的量测值;x为系统状态;hi(x)为测点i的量测函数;ei为测点i的量测误差;m为系统中测点数目。

在实际应用中,对测点i,若

$|e{}_i|=|{\rm Z}{}_i-h{}_i(x)|\le \alpha {}_i(2)$

则称测点i合格,反之则不合格。αi为给定常量,可按电力公司已有标准选取,本文采用值见附录A。

经典最小二乘状态估计方法的基本思想是在已知量测情况下,希望找到某个状态x,使得残差平方和最小。这样,当某个测点量测数据为坏数据时,其残差较大,对最小二乘结果影响也较大,造成所谓“残差污染”。实际应用中,为克服这一问题,往往要通过某些方法[8,9],判断出该量测数据为坏数据,计算时给予较小的权重,来减少坏数据对估计结果的影响(即提高“抗差性”)。这种方法,需要主观确定加权因子,调试和维护较为困难。

对这一问题,本文采用一种全新思路:①任给系统某个状态x。②在这一状态下,对任一测点i,若其是合格测点(满足式(2)),则认为该测点“投票同意”这一状态;反之,认为该测点“投票反对”这一状态。③认为最多测点赞同的状态即为待求的系统状态。显然,在该状态下,系统合格点数最多(即合格率最大)。基于以上思路,本文提出一种新的状态估计方法。

1.2 测点评价函数

根据上文思想,可建立如下测点评价函数:

$g{}_i(e{}_i)=\{\begin{array}{l}0|e{}_i|\le \alpha {}_i\\ 1|e{}_i|>\alpha {}_i\end{array}(3)$

式(3)意义为,当测点i合格(|ei|≤αi)时,gi(ei)为0;反之,gi(ei)为1。gi(ei)形状见图1。

函数gi(ei)的缺点是:在±αi处有间断点,并且不是处处连续可导,不便于实际应用。

为此,建立gi(ei)的近似函数,令

$f{}_i(e{}_i)=\frac{1}{1+\text{e}{}^{-\frac{(e{}_i-\alpha {}_i)c}{\alpha {}_i}}}+\frac{1}{1+\text{e}{}^{\frac{(e{}_i+\alpha {}_i)c}{\alpha {}_i}}}(4)$

式中:c一般为大于3的常数。

fi(ei)形状如图2所示(图中αi=10,c=10)。关于参数取值与算法收敛性关系讨论参见算例。

由图2可见,fi(ei)没有间断点,处处连续可导,且fi(ei)与gi(ei)的特性接近,即当测点i合格时,fi(ei)接近于0;反之,fi(ei)接近于1。因此,本文采用式(4)定义的fi(ei)作为测点评价函数。需要说明的是,测点评价函数也可采用其他形式,只要保证处处连续可导,且当测点i合格时,测点评价函数值接近于0,反之接近于1即可。

1.3 数学模型

设系统中不合格测点数目为k,则${\displaystyle \sum _{i=1}^{m}f}{}_i(e{}_i)\approx k$。而要寻求某一系统状态,使得合格测点数最多,即相当于寻求某一系统状态,使得在该状态下${\displaystyle \sum _{i=1}^{m}f}{}_i(e{}_i)$最小。得到新的状态估计模型如下:

$\{\begin{array}{l}\underset{x}{{\displaystyle \min }}{\displaystyle \sum _{i=1}^{m}f}{}_i(e{}_i)\\ \text{s}.\text{t}.{\rm Z}{}_i=h{}_i(x)+e{}_i\forall i=1,2,\cdots ,m\end{array}(5)$

由于实际运行中,真实系统状态必然满足潮流约束和运行中上下界约束,考虑这些约束后,得到改进的状态估计模型:

$\{\begin{array}{l}\underset{x}{{\displaystyle \min }}{\displaystyle \sum _{i=1}^{m}f}{}_i(e{}_i)\\ \text{s}.\text{t}.{\rm Z}{}_i=h{}_i(x)+e{}_{i}i=1,2,\cdots ,m\\ g(x)=0\\ \underset{¯}{h}\le h{}_i(x)\le \overline{h}i=1,2,\cdots ,m\end{array}(6)$

式中:g(x)=0,代表潮流约束;$\underset{¯}{h}\le h{}_i(x)\le \overline{h}$,代表各测点运行中上下界约束。

1.4 求解算法

式(6)中,将ei代入目标函数中,即得:

$\{\begin{array}{l}\underset{x}{{\displaystyle \min }}{\displaystyle \sum _{i=1}^{m}f}{}_i({\rm Z}{}_i-h{}_i(x))\\ \text{s}.\text{t}.g(x)=0\\ \underset{¯}{h}\le h{}_i(x)\le \overline{h}i=1,2,\cdots ,m\end{array}(7)$

式(7)是一最优潮流问题,可用求解最优潮流的算法求解。由于现代内点法在最优潮流求解中的成功应用[10],最优潮流的收敛性和计算速度均得到保证,因而模型(式(7))也可用现代内点法顺利求解。

2 方法特点

与以往状态估计方法相比,本文方法具有以下特点:

1)估计准确性不易受不良数据影响,具有很强的抗差性。以往状态估计方法易受不良数据影响,且估计值偏离量测值越远,受影响就越大。本文方法中,当测点不合格(不满足式(2))时,无论估计值偏离量测值多远,其反映在目标函数中大小都是1,因而估计准确性不易受不良数据影响,抗差性强。

2)可自动对坏数据点进行识别。本文方法所估计状态为绝大多数测点所赞同的状态。由此推论:若绝大多数点赞同某一状态,则投票反对该状态的测点必为坏量测点。因而,本文方法在进行状态估计时,自动对坏数据点进行了识别。

3)其他特点。①所求得的状态估计结果为潮流解,且满足各种物理约束,更加接近实际;②计算中无需做坏数据校验、可观性校验、权重因子设置,调试和维护极为简单;③求解算法为现代内点法,收敛性好、计算速度快。

3 算例

分别用文献[11]中4节点系统、IEEE 39节点系统、IEEE 118节点系统和某实际输电网的数据进行试验。对IEEE标准系统试验中,通过将潮流计算结果加上2%的高斯噪声,得到不含不良数据的生数据,通过将生数据改变符号、置零或加减量测值20%以上等得到试验用的不良数据。试验中,如无特别说明,WLS方法中测点i的权重取1/σ${}_i^2$,σi为正态分布标准差,本文方法中取αi=3σi。

3.14节点试验系统

4节点系统量测配置见附录B图B1。在正常量测基础上将P12符号取反作为不良数据,表1为试验结果。可以看出:WLS方法受不良数据的影响导致量测P1,P12等估计结果偏离真值较远,本文方法在案例中表现出抗差性,且估计值都逼近真值。

3.2 IEEE39节点系统

采用以下的指标来衡量估计的状态变量(取平衡节点外其余节点的电压和相角)的准确性[12]:

$\begin{array}{l}S{}_1={\displaystyle \sum _i|}x{}_i^{*}-x{}_i|\\ S{}_2=\max |x{}_i^{*}-x{}_i|\end{array}$

式中:x*i为第i个状态变量的估计值;xi为第i个状态变量的真值。

表2给出了试验结果,IEEE 118节点的试验结果见附录C表C1。

可以看出:本文方法随不良数据所占比例增加,S1 和S2指标变化较小,不易受不良数据的影响。

3.3 某实际输电网实际数据

某实际输电网连续5个断面的试验结果见表3。计算时各测点权重取目前该系统实际采用的权重值:500 kV,220 kV和其他电压等级的电压测点权重分别取10,2.236,1.0;有功和无功测点权重分别取1.0和0.2。αi的取值见附录A。从表3中可以看出,本文提出的状态估计方法在合格率方面较传统状态估计方法有了明显提高。

去除不合格测点后再进行状态估计,合格率均为100%,说明本文方法可自动识别不合格测点。

3.4 收敛性试验

试验中,使用IEEE 118节点系统,各次试验使用相同初值,计算模型均采用模型(式(7)),其中c取不同值时,迭代次数和状态估计结果评价指标如表4所示,其他测试结果见附录C表C2和表C3。

从测试结果中可以看出:无论模型中是否含有不良数据,且无论含有多少不良数据,均有迭代次数少、解的合理性好,且迭代次数和合理性指标变化不大,这也说明了本文方法的稳定性和合理性。

3.5 测点合格率标准对估计结果的影响

试验中,使用IEEE 39节点系统,测点中含有2.38%的不良数据,模型(式(7))中参数c取3。表5给出了αi取不同值时的试验结果。

可以看出,式(2)中的αi取0.5倍～6.0倍标准差时,估计结果均比较准确,本文方法对合格率标准不敏感。

4 结语

基于全新理念,建立了以合格率最大为目标的状态估计模型,该方法具有以下特点:①所求得的状态估计结果合格率更高;②估计准确性不易受不良数据影响,具有很强的抗差性;③可自动对坏数据点进行识别;④所求得的状态估计结果为潮流解,且满足各种物理约束,更加接近实际;⑤计算中无需做坏数据检验、可观性校验、权重因子设置,调试和维护极为简单;⑥求解算法为现代内点法,收敛性好,计算速度快。

用4节点系统、IEEE 39节点系统、IEEE 118节点系统和某实际输电网实际数据的算例证明了本文方法的有效性。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

摘要：以合格率最大为目标,考虑潮流约束和运行中上下界约束,提出一种电力系统状态估计方法。该方法可用现代内点法求解,收敛性好。与以往方法相比,文中方法得到的状态估计结果合格率高;估计结果有较强的抗差性;所求的解为潮流解,满足各种物理约束,更加接近实际;计算中无需做坏数据检验、可观性校验,无需主观设定测点权重,大幅减少了调试和维护工作量。IEEE39节点系统、IEEE118节点系统和某实际输电网等试验系统的测试结果验证了该方法的优越性。

关键词：状态估计,合格率,加权最小二乘,抗差估计器,电力系统

参考文献

[1]SCHWEPPE F C,WILDES J.Power system static-state estimation:PartⅠ-Ⅲ.IEEE Trans on Power Apparatus and Systems,1970,89(1):120-135.

[2]BALDICK R,CLEMENTS K A,PI NKO-I AZIGAL Z,et al.Implementing non-quadratic objective functions for state estimation and bad data rejection.IEEE Trans on Power Apparatus and Systems,1997,12(1):376-382.

[3]CELIK M K,ABUR A.Arobust WLAVstate estimator using transformations.IEEE Trans on Power Systems,1992,7(1):106-113.

[4]SI NGH H,ALVARADO F L,LI U W H E.Constrained LAV state estimation using penalty functions.IEEE Trans on Power Systems,1997,12(1):383-388.

[5]MILI L,CHENI AL M G,ROUSSEEUW P J.Robust state estimation of electric power systems.IEEE Trans on Circuits and SystemsⅠ:Fundamental Theory and Applications,1994,41(5):349-358.

[6]李碧君,薛禹胜,顾锦汶,等.电力系统状态估计问题的研究现状和展望.电力系统自动化,1998,22(11):53-60.LI Bijun,XUE Yusheng,GUO Jinwen,et al.Status quo and prospect of power system state estimation.Automation of Electric Power Systems,1998,22(11):53-60.

[7]W CHTER A.An interior point algorithm for large-scale nonlinear optimization with applications in process engineering[D].Pittsburgh,PA,USA:Carnegie Mellon University,2002.

[8]ZHANG B M,LO K L.A recursive measurement error estimationidentification method for bad data analysis in power systemstate estimation.IEEE Trans on Power Systems,1991,6(1):191-197.

[9]ZHANG B M,WANG S Y,XI ANG N D.Alinear recursive bad data identification method with real-time application to power systemstate estimation.IEEE Trans on Power Systems,1992,7(3):1378-1385.

[10]WEI H,SASAKI H,KUBOKAWAJ,et al.Aninterior point nonlinear programmingfor optimal power flowproblems with a novel data structure.IEEE Trans on Power Systems,1998,13(3):870-877.

[11]于尔铿.电力系统状态估计.北京:水利水电出版社,1985.