参数估计理论与应用第三章

参数估计理论与应用第三章3、1参数估计得评价准则 参数估计就是通过样本去估计总体得某些数字特征或统计量。任何一个统计量都可作为参数得估计量,但其效果得优劣有所差别。3、1、1无偏性、有效性与相容性

(1)无偏性设样本得总体分布密度函数为p(x;θ),θ就是未知参数。从总体中抽取容量为N得样本x={x1,…,xN

},用样本得估计量来估计θ,如果希望多次估计中,平均得估计值没有偏差,即

则称就是θ得无偏估计量。10/18/2024 例3-1

样本均值就是总体数学期望得无偏估计。 设x1,…,xN就是随机过程{xk}得N个独立观测样本,如果参数θ就是总体得数学期望E[x],即用样本得均值作为θ得估计量,对该估计量取期望值,有 一个无偏估计量在多次估计中将不会产生系统偏差,但并不意味着有偏估计就不好。如果一个有偏估计就是渐进无偏得,即10/18/2024那么她仍然有可能就是一个好得估计。 考虑实随机过程{xk}得相关函数得两种估计量:

假定数据{xk}就是独立观测得,容易验证

式中,Rx(τ)=E[xk+τxk]就是随机数据{xk}得相关函数。 以上二式表明,估计量1(τ)就是无偏得,而2(τ)则就是有偏得。但就是,2(τ)就是渐进无偏得,即10/18/2024渐进无偏估计量2(τ)就是半正定得,而无偏估计量1(τ)却不一定就是半正定得,故2(τ)得使用场合较多。

(2)有效性 如果1

和2

就是两个根据N个独立观测样本得到得无偏估计量,无疑地,对θ得平均偏差较小就是选择得标准之一。例如,如果则

1得值比2得值更密集地聚集在真值θ得附近。通常将方差(或协方差阵)在所有得无偏估计量中达到最小得

称为有效估计量。

例3-2设x1,…,xN就是N个独立观测样本,若被估计参数10/18/2024θ=E[x],则对任何满足都就是θ得无偏估计量。利用不等式可得在估计总体得数学期望时,简单得算术平均比加权平均好。 (3)一致性估计量得精度就是与样本得容量N有关系得。一般说来,总就是认为N越大估计得效果应该越好。如果记依赖样本容量N得估计为N

,当满足10/18/2024则称N就是θ得一致性估计量,或相容估计。

例3-3设总体x具有均匀分布,分布密度为其中,θ1

和θ2

就是未知参数。

总体样本得均值和二阶矩分别为(严格按定义计算)解得10/18/2024 按矩得估计方法,用独立样本得均值和独立样本得二阶矩,分别作为总体均值和二阶矩得估计量,就有 下面说明1

和2

分别就是θ1

和θ2得相容估计。 设y1,…,yN

就是具有同分布得独立观测样本,根据大数定律,有令y=x2,就有10/18/2024于就是3、1、2Fisher信息和Cramer-Rao不等式

通常希望获得有效得参数估计量。但就是,由于不存在导致最小方差无偏估计量得最佳算法,所以通常采用参数无偏估计得Cramer-Rao下限(或CR下界),作为评价参数估计性能得测度。为了简洁叙述这一得评价测度,先定义一个重要得概念。

Fisher信息Fisher信息用J(θ)表示,定义为(3、1、1)10/18/2024大家学习辛苦了，还是要坚持继续保持安静 当考虑N个观测样本X={x1,…,xN},此时,联合条件分布密度函数可表示为 将式(3、1、1)中得p(x|θ)改为p(X|θ)就可给出N个样本变量X得Fisher信息得表达式。

定理(Cramer-Rao不等式)设观测样本X={x1,…,xN

},若参数估计就是真实参数θ得无偏估计,并且条件分布密度函数得p(X|θ)对参数θ得一、二阶偏导数存在,则有(3、1、2) 参数得方差所能达到得下限(称为CR下限),即上式等号成立得充要条件就是10/18/2024其中,函数K(θ)>0,并与样本向量X无关。 当为有偏估计量时,Cramer-Rao不等式为

(3、1、3)

式中η(θ)为估计偏差,即η(θ)=E[]-θ,并假定b(θ)就是可微分得。 对于多个参数得情况,记θ={θ1,…,θp},则用矩阵J(θ)表示Fisher信息,其元素Jij(θ)定义为(3、1、4)10/18/2024且Cramer-Rao不等式变为矩阵不等式:(3、1、5)

上式表示无偏估计量得协方差矩阵cov()与逆Fisher信息阵之差就是一半正定矩阵。 Fisher信息就是描述从观测数据中得到得θ得“信息”测度,她给出利用观测数据估计参数θ得方差下界。但就是,满足这一下界得估计量有得时候可能不存在。3、2基于统计分布得参数估计方法 参数估计量得优劣取决于所采用得评价准则(或代价函数)和估计算法。现在介绍已知总体统计分布得两种最有效得参数估计方法:Bayes估计和最大似然估计。10/18/20243、2、1Bayes估计

在参数估计中,估计误差θ-通常不为零。因此,除了采用前面介绍得无偏、有效和相容估计作为评价准则外,还可以利用估计误差得变化范围作为参数估计得测度,这种测度叫做代价函数,用符号C(,θ)表示。常用得代价函数有绝对型、二次型和均匀型三种。OOO∆/2∆/2绝对型二次型均匀型10/18/2024 本节仅介绍最常用得二次型代价函数,即 当总体得分布密度函数p(X|θ)已知时,利用X={x1,…,xN

}进行参数估计,通常就是采用代价函数得期望值作为评价参数估计量效果得测度,并称之为风险函数。使风险函数最小得参数估计叫做Bayes估计;基于二次型风险函数最小得估计称为最小均方误差(minimummeansquareerror,MMSE)估计。二次型风险函数定义为(3、2、1) 根据条件概率公式,有10/18/2024其中,p(θ|x1,…,xN

)就是给定N个观测样本X={x1,…,xN}条件下θ得后验分布密度函数。于就是,式(3、2、1)可以写成(3、2、2) 为使风险函数RM

MSE最小,对上式取得偏导,并令其结果为零,便得到由于p(x1,…,xN

)就是非负得,因此,∂RM

MSE

/∂=0,等价于上式中[·]=0。故有10/18/2024(3、2、3)

注意,在式(3、2、3)中,利用了以下事实: 由此可得出重要得结论:未知参数θ得MMSE估计就是给定样本X条件下θ得条件均值。

例3-4某一随机参量x服从高斯N(mx,Cx)分布,用仪器可测量其线性组合y

,即(1)式中,y－N维,k－N×M维,x－M维,e－N维。10/18/2024其中,测量误差e服从高斯N(0,Ce)分布;k为给定得常数阵。假设 (ⅰ)e与x

独立;(ⅱ)e与x

相关,互协方差函数为Cxe

。 试分别求出两种情况下得MMSE估计xˆ(y)和估计误差x

(y)得协方差Rx

(y)。

解如果将x

看作未知参数,那么,根据上面讨论,x得MMSE估计就是给定观测样本{y1,…,yN}时x得条件均值。因此,可利用公式(1、4、16)和(1、4、17)[pp、29]

(2)(3)来求解。10/18/2024 对式(1)两边取均值,得到(4)

将式(1)和(3)代入有关定义式,得(5)(6)(7)10/18/2024

(i)当e与x

互相独立,Cxe=0。将式(4)~(7)代入式(2)和(3),得到xˆ(y)得估计及协方差Rx

(y)

(ii)当e

与x

相关,只需注意Cxe

≠0即可。 这个问题留给读者解决。请构造一组数据,在Matlab平台上仿真这两种得估计结果。3、2、2最大似然估计

最大似然估计(maximumlikelihoodestimate,ML估计)得基本思路就是:在给定参数θ条件下,将观测样本xK10/18/2024联合条件概率密度函数p(x|θ)视为真实参数θ得函数,即似然函数L(x,θ)

(包含未知参数θ得可能性函数),然后利用容量为N得观测样本x={x1,…,xN

},求出使L(x,θ)达到最大化得参数作为θ={θ1,…,θp}得估计值。在数学上,通常将未知参数θ得最大似然估计量记为式中Θ就是参数θ得值域。故ML估计量ML就就是p(x|θ)得全局极大点。 由于对数函数就是严格单调得,故L(x,θ)得极大点与lnL(x,θ)得极大点就是一致得。通常,将lnL(x,θ)称为对数似然函数。于就是,ML估计量ML可由(3、2、4)10/18/2024确定。如果x1,…,xN就是N个独立得观测样本,则对数似然函数可写作(3、2、5)

ML估计量ML只要能够求出来,总就是比较好得估计,她具有以下性质:最大似然估计就是有效和一致估计;对于大得N,ML估计量ML服从高斯分布,并且就是无偏得,方差可达CR下界。

例3-5设样本x={x1,…,xN

}服从高斯分布N(m,σ),则其对数似然函数为10/18/2024分别求lnL关于m和σ2

得偏导,并令她们等于零,得到解得显然有 可见,均值得ML估计ML

就是无偏得,而方差得ML估计

ML就是有偏得。但若将ML

·N/(N-1)作为新得估计量,则该估计就是无偏得。10/18/2024

计算L(x,θ)得相对于m得二阶偏导数,有由式(3、1、1)得Fisher信息:Cramer-Rao不等式为等号成立得充要条件就是 事实上,我们有10/18/2024因此,只要取K(m)=N/σ2,ML估计ML就可达CR下界σ2/N。这表明ML估计ML就是一有效估计量。

例3-6(二元阵最大似然测向系统)设二元阵布置在x轴上,两个基元坐标分别为x1

和x2,如图3-2所示。如果取x1=0,则x2=d,d为两传感器得位置间隔。假设信号为平面波,入射角为θ,则传感器1相对于传感器2得信号时延τ为(3、2、6)式中,c为声速。我们得问题就是如何利用二元阵中两个输入过程得时差τ来测定目标得方位角θ。θxx2=dx1=0图3-2二元阵测向系统的几何关系10/18/2024 解

设两传感器得零均值接收过程可分别表示为其中,si为单频平面波信号,wi(i=1,2)为零均值高斯噪声,二者互相独立。 如果采用图3-3所示得时延补偿方法,则单频平面波信号得归一化声程补偿(或指向)向量v在所考虑得二元阵中可表示为

下面,我们来推导信号得协方差矩阵和噪声得协方差矩阵,以便于求出观测样本得似v*图3-3声程补偿系统x2x11exp(-jωnτ)∑10/18/2024函数。记输入信号和输入噪声得傅立叶系数为设信号和噪声得功率谱分别为S(ωn)和N(ωn),那么,由公式(1、4、6)[pp、26,ωn=2πn/T)]

信号和噪声得协方差矩阵可分别表示为(3、2、7)10/18/2024于就是,观测样本得似然函数可表示为

(3、2、8)式中,X(1)=[X1

(1),X2(1)]T

,…,X(TW)=[X1(TW),X2(TW)]T就是传感器得接收过程{x=[x1,x2]T}得傅立叶系数阵;T就是过程得持续时间(采样数据得长度),W就是接收过程得带宽。 容易验证,行列式|Cw

+Cs

|与时延τ无关。于就是,ML估计就就是选择τ,使lnp(X|τ)最大,也即使式(3、2、8)得指数函数(3、2、9)10/18/2024最大。下面,我们从式(3、2、9)出发,推导时延参数τ得最大似然估计得等效形式。为此,首先引进下列求逆公式(3、2、10)式中,A为n×n非奇异矩阵;g为n×1列向量。证明留给请读者课外练习【利用恒等式g(1+gHA-1g)=(A+gHg)A-1g)】。 利用求逆公式,可知<g=[1exp(jωnτ)]>10/18/2024将上式代入式(3、2、9),略去与τ无关得量T/N(ωn)。因此,选择τ使式(3、2、9)最大,等价于使下式(3、2、11)

最大。现引入记号在此将X(ωn,τ)视为某时间函数x(t,τ)在时间(t-T,t)内 得傅立叶系数。将上述替换量代入式(3、2、11)后,再应用 周期函数得Parseval公式,就有10/18/2024略去无关紧要得常数项1/2,计算z(x,τ)得结构如图3-4所示。调节时延τ,使输出z(x,τ)达到最大,相应得时延就就是 真实时延得ML估计ML。 根据ML估计得传递性,由式(3、2、6)可得真实方位得ML估计

(3、2、12)

xH0(t)z(x,τ)∑x1(t)x2(t)H0(ω)(·)2图3-4二元阵最大似然测向系统exp(-jωτ)10/18/2024 二元阵最大似然测向系统与二元阵似然比检测系统具有完全相同得结构。这就是因为:在H1情况下,p(X|τ)等价于p(X|H1),后者也可看作就是时延参数τ得函数;而在H0

情况下,p(X|H0)与τ无关。因此,选取τ使似然函数最大,也就就是使似然比

p(X|H1)/p(X|H0)最大。由此可见,检测问题与参数估计问题就是密切相关。 顺便指出,可用测向测距近似公式(3、2、13) 构成最大似然联合测向测距系统。其中,di表示第i个传感器与“基准”传感器位置得间距;D表示目标与“基准”传感器位置之间得距离。10/18/20243、3基于模型得参数最小二乘估计

最小二乘法(Leastsquaremethod,LS)就是一种不需要任何先验知识得参数估计方法。在被测系统得静态(稳态)模型和动态模型得参数辨识中,最小二乘法就是最常用得参数估计方法,在测控技术领域获得了广泛得应用。3、3、1最小二乘估计器及其统计特性

在一般得最小二乘问题中,线性系统得参数化模型可以表示为(3、3、1)

其中,u=[u1,…,up]T就是模型得输入向量,f1,…,fn就是u得已知函数,也可以就是未知输入得观测数据;θ1,…,θn就是待估计10/18/2024得参数,又称为回归系数;y就是系统得输出。

当f1,…,fn就是u得稳态响应状态或就是实测得确定性变量,且y就是系统得稳态输出,则称式(3、3、1)就是描述线性系统得静态模型;当y就是u得动态响应或瞬态观测数据,那末式(3、3、1)就就是描述线性系统得动态模型。 为了估计未知参数θi,必须做实验来获得数据对{[uiyi]或[fk

(ui)yi],i=1,2,…,N,k=1,2,…,n;N≥n}以构成训练数据。将数据对代入方程(3、3、1),可以获得一组线性方程: 用矩阵表示方法,将上式写成更简洁得形式,即10/18/2024

(3、3、2)其中

为了唯一地识别出未知参数,通常要求N>n,即数据对得数目多于拟合参数得数目。满足所有N个方程得精确解就是不可能得,因为观测数据难免受到噪声得污染,或者描述系统得参数化数学模型不够精确。故必须考虑随机噪声或建模误差,在方程(3、3、2)中引入随机误差向量e,得到

(3、3、3)10/18/2024 参数θ得最小二乘估计LS

,就就是使目标函数

(3、3、4)

达到最小值得参数估计。为此,通常都采用求极值得方法。 将式(3、3、4)展开后,得到

对θ求导数,有 J

极小化得条件就是一般均假设ΦTΦ非奇异,于就是,LS有唯一得解:10/18/2024

(3、3、5)

式中Φ+表示Φ得伪逆。 上述表示误差向量对整体平方误差有相同权重。可以进一步扩展,令每个误差项有不同得权重。设W为所需得权值矩阵,她就是对称和正定得,则加权得目标函数为(3、3、6) 按上述求极小值得方法,可得加权得最小二乘估计量:(3、3、7)显然,当W选为单位矩阵时,WLS

=LS。

例3-7考虑最简单得一维线性模型(静态得),即只有一个控制变量u得情形,这时模型得形式就是10/18/2024求未知参数θ0

和θ1得LS估计量。

解实际过程输出就是模型得输出加上一随机误差项,即观测数据对[ui,yi]得结构应为式中,ei称为模型得残差或观测噪声,一般认为就是零均值、相互独立得随机序列,并具有相同得方差σ2。将上式写成矩阵形式:

根据式(3、3、5),可得LS估计量:10/18/2024 如果进一步假定ei得分布就是正态得,则容易验证,方差σ2

得ML估计量就是 作为练习,请读者在Matlab平台上输入以下数据和函数:x=[12345];y=[1、31、82、22、93、5];

[p,s]=polyfit(x,y,1) %生成拟合一次多项式运行结果就是:p=[0、550、69],s=0、1643。即y=0、55x+0、69标准差为0、1643。10/18/2024

例5-8(可线性化得非线性静态模型——曲线回归)假设有一个非线性模型得输出为其中,x1,x2

为确定性输入变量,a,b和c为待估计参数。

解上式两边经简单得代数运算,再同时取自然对数,可转化为一个线性模型:这说明变换后得输出ln(y-1-1)可以显式地表达为以lnx1和x2为输入、以lna,b和c为参数得线性模型。因此,就可以按变换后得线性模型用最小二乘法来估计变换后得未知参数,然后,再根据变换后得估计参数计算出原参数。10/18/2024 判定输入x-输出y之间得关系能否用一个线性模型来描述得标准,通常用互相关系数得大小来衡量:(3、3、8)ρxy得绝对值越大,表示变量之间得线性关系越密切,因而线性回归得效果就越好。 例3-9设某一结构参数n,m和d已知得离散线性系统,其差分方程得形式为:

(3、3、9)10/18/2024其中,e(k)为噪声,φ(k)为输入-输出观测向量,θ为未知参数向量,且要求根据N次数据对{[y(i),u(i)],i=1,2,…,N;N≥n+m+1}来估计对未知参数θ。

解将式(3、3、9)改写成矩阵形式,得到将数据写成下标形式,就有这样,未知参数向量θ可按式(3、3、5)进行估计。 10/18/2024

考虑如下单输入-单输出系统:用Matlab中rarx函数进行系统辨识,程序如下:

A=[1-1、50、7]; %a0=1,a1=-1、5,a2=0、7 B=[00、30、20、5]; %b0=0,b1=0、3,b2=0、2,b3=0、5 th0=arx2th(A,B,1,1); %实际系统得ARX模型 e=randn(200,1);u=idinput(200,‘prbs’);%高斯噪声和伪随机信号 y=idsim([ue],th0);z=[yu]; %模型仿真;输入-输出信号[z] na=2;nb=3;nk=1 %ARX模型得阶次 [thm,yhat]=rarx(z,[nanbnk],'ng',0、1);%根据[z]进行ARX模型参数辨识 plot(y,'-');grid %作图,实际系统得输出曲线 holdon plot(yhat,':') %作图,辨识系统得输出曲线

参数辨识结果thm:aˆ1=-1、3798,aˆ2=0、7039,bˆ1=0、3007, bˆ2=0、1170,bˆ3=0、4243。10/18/2024 应当指出,要求观测数据容量N≥n+m+1就是为了保证ΦTΦ非奇异,降低过程噪声序列{e(k)}得影响,从而提高参数估计得精度。不论{e(k)}就是何种形式得噪声序列,式(3、3、5)总就是成立得。换言之,噪声性质仅影响LS估计得统计特性。 下面介绍LS估计得统计特性。如果观测噪声或建模误差序列{e(k)}具有零均值和相同得方差,即则LS估计量LS就是无偏、有效和相容得,并具估计误差得协方差为σ

2(ΦTΦ)-1。 对于动态控制系统得辨识,输入信号u(t)必须满足持续激励条件,也即输入信号u(t)得频谱必须包含足够丰富(Sufficientrich)得频率成分,以保证充分激励受控对象得10/18/2024所有振型,从而使观测数据载有动态系统得主要信息。LS估计在满足持续激励条件时,就是渐进无偏得,也称为估计得一致性。 式(3、3、9),也称为CAR模型(即受控得AR模型),可以写成更简洁得形式(3、3、10)式中q表示时间算子,d为整数,表示系统得滞后量;A(·),B(·)分别为q－1

得降次幂多项式。 CAR模型满足一致估计(或相容估计)得条件为:

(1){e(k)}就是白噪声序列; (2){u(k)}得均值和协方差有界;且满足(m+1)阶持续激励条件(或正定条件):10/18/2024

(3)u(k)和e(k)相互独立。通常,u(k)都采用伪随机二元序列。 只要选择恰当得模型阶次或最小二乘多项式阶次(参见taylor、m,Matlab),最小二乘法总就是可以很好地拟合数据,但就是,如果观测数据波动较大,将严重影响参数估计得准确性。对此,可采用数据预处理和数字滤波得方法加以解决。 检验模型准确性得最简单方法就是准备另外一组输入-输出数据对,称为检验数据集,在参数估计时不用,待模型建立后,用这组数据对来验证所得模型得普适性或泛化能力。10/18/2024

上述讨论,均假设数学模型得阶次就是已知得。实际上,对于动态系统,模型得阶次很少就是预知得。检验模型阶次就是否合适得一种简单而有效得方法就是:评估不同阶次得理论模型对观测数据得拟合度,用拟合误差函数

来描述。通常,当n或(n,m)增大时J(n)或J(n,m)就会减小;而当n(或n,m)大于模型得真实阶次n0(或n0,m0)时,J得减小就不显著了。由此,可以很方便地用多次实验得方法来确定模型得恰当阶次。 注意,对于J(n,m)形式得拟合误差函数,一般应按正交实验法来确定模型得恰当阶次,以减少实验得次数。10/18/20243、3、2递推最小二乘估计 在测控系统中,被测对象通常可以不断提供新得输入-输出数据。如果希望利用新得信息来改善估计精度,那么就应当采用递推估计算法,这不仅可避免观测数据矩阵Φ得行数得不断“膨胀”,而且可减少参数估计得计算量。 在推导最小二乘递推算法前,先引入一个与式(3、2、10)类似得矩阵求逆定理。设A和I+CA-1B均就是非奇异方阵,则

(3、3、11)

下面介绍最小二乘递推算法。为了简化符号,以下推导均用代替LS。 设ΦN就是时刻N为止得观测数据,N+1时刻θ得LS估计为10/18/2024式中[参见式(3、3、2)和(3、3、9)]于就是有(3、3、12) 令PN=[ΦN

TΦN]-1

,由求逆公式(3、3、11)知

(3、3、13)

定义增益向量为10/18/2024将式(3、3、13)和(3、3、14)代入(3、3、12),得到(3、3、15)

上式表明,新得估计量N+1等于前一时刻得估计量N

与修正项KN+1(yN+1-φTN+1N)之和,这就是一切递推公式得共同特征。如果令代表基于前一时刻得估计量N

对N+1时刻得预测。那么,递推估计提供得新息10/18/2024就就是预测误差或拟合误差。因此,修正量得大小与新息成正比,而各校正分量得权,由增益向量决定。 在启动上述递推算法时,必须知道初值0和P0

,通常令其中,σ2»1。然后从得到得一组数据,按式(3、3、16)开始递推运算。 从物理上看,这种初值选取虽然初始误差较大,但校正得作用也大,因此这种递推算法就是有效得。此外,还可以先取得N>m+n+1组数据,算出10/18/2024作为初值,然后,按式(3、3、16)进行递推运算。 增益向量KN+1在递推运算过程就是怎样变化得？先考察PN=[ΦN

TΦN]-1如果N大于估计参数得数目,而且,输入－输出数据对含有足够得“信息”(满足充分激励条件),则 ΦNTΦN

通常就是正定得。显然,当N趋于无穷大,ΦNTΦN

/N接近于非奇异得常数阵。于就是,有 可见,式(3、3、16)中自适应增益向量KN+1

随着每次迭代而递减,这意味着递推运算过程将逐渐收敛于参数空间得最优点。事实上,在白噪声或低噪声条件下,递推最小二乘10/18/2024估计就是一种简便而又有效得算法。这种递推算法在递推过程中虽然没有保存全部先前得数据,但所有先前数据得影响却一直在起作用,故称为无限增长记忆得递推最小二乘估计。3、3、3卡尔曼滤波器得递推算法(状态估计) 与参数估计器不同,卡尔曼滤波器主要就是解决如何从被噪声污染得观测数据中估计出已知动态系统模型得状态,而不就是动态系统模型得未知参数。然而,仅从算法上看,这二者就是非常相似得。为了便于比较二者得异同,我们在此不加证明地列出卡尔曼滤波器算法。 设离散定常系统得状态方程和观测方程分别为10/18/2024式中,X(k)就是n维状态向量,u(k)就是m维控制向量,w(k)就是

p维过程噪声,y(k)就是r维观测向量,v(k)就是r维观测噪声;

A,B,C,Γ分别就是相应维数得系统矩阵、控制矩阵、观测矩阵和过程噪声权矩阵。噪声得统计特性满足

假定到k时刻为止,观测数据为{y(1),…,y(k)},要估计l时刻得状态,就有三种情况:l>k,称为预测问题;l=k,称为滤波问题;l<k,称为平滑问题。 下面主要介绍滤波问题。10/18/2024滤波算法:预测算法:滤波增益:滤波误差得协方差:预测误差得协方差:初始条件:10/18/2024 例3-10考虑某一动态系统表示得空间导航问题,其中加速度为白噪声。被白噪声污染之后观测其位置。因此,过程得状态方程为(1)式中,X(t)=[x1(t),x2(t)]T,表示过程得位置和速度;而w(t)就是一均值为0、方差为1得高斯白噪声。观测方程为(2)其中v(t)就是一均值为0、方差为10得高斯白噪声。

解首先把连续得状态方程离散化(参见Matlab中c2d函数)10/18/2024不妨设采样间隔T=1,则可写出状态方程(1)和观测方程(2)得离散化表达式进一步假设则可按前面介绍得滤波算法进行递推计算:Xˆ0→Xˆ(1|0),P(1|0)→K(1)→Xˆ(1|1);Xˆ(2|1),P(1|1)→P(2|1)→K(2)→Xˆ(2|2)…。 在Matlab中用Kalman函数仿真卡尔曼滤波器得设计。请读者在Matlab平台上完成例3-10得仿真计算。10/18/20243、3、4限定记忆得递推最小二乘估计

以上介绍得递推最小二乘法适用于定常系统得参数估计。对于时变系统,由于参数时变得信息显然更多地蕴藏在当前得观测数据中,而与先前观测数据得关系将逐渐减弱。因此,利用一切观测数据对来决定得自适应增益向量KN+1,显然会削弱递推过程跟踪变化参数得能力。解决这一问题得一种简单方法就是:当我们怀疑观测数据发生显著得变化时,就将当前得PN

设置为P0,重新进行参数估计,这就是因为LS估计能快速地收敛到当前得最优参数。处理这一问题得另一种方法就是对过去得数据引入带遗忘因子λ,逐渐削弱她们在参数估计中得作用。为此,采用加权得目标函数:10/18/2024其中(3、3、17)

而且0<λ≤1,当λ=1,就退化为基本得递推最小二乘法。由式(3、3、7)知,在时刻N,使上述加权得目标函数最小化得估计为(3、3、18)

每当取得一个新得数据后,就对沟渠以前得加权矩阵乘以λ,于就是N+1时刻得估计量为10/18/2024与前面略有不同,在此令PN=[ΦN

TWNΦN]-1

,则由求逆公式(3、3、11),可知于就是,按前面推导出式(3、3、15)思路,可得式中

遗忘因子λ将老得数据逐渐从“记忆”中去掉,因此这种使用数据信息得方式也叫做“渐消记忆”法,相应得算法称为带遗忘因子得递推最小二乘法,即10/18/2024关于λ得选取通常由经验或实验确定,一般范围为0、95≤λ≤0、99 λ取得愈小,最新数据得权重就愈大,也就更适合于跟踪大得时变参数,但与此同时,估计器也可能会发生较大得波动 从而加大估计得误差。 例5-10考虑模型采样300次后变为10/18/2024

试用两组模拟数据,一组不考虑噪声,一组就是带观测噪声得数据,分别用不同得遗忘因子,对时变模型进行参数估计,并讨论估计结果。

解用MATLABrarx函数进行带遗忘因子λ得系统辨识算法。程序如下: e=randn(300,1);u=idinput(300,‘prbs’); %产生高斯噪声和伪随机信号 A1=[10、8];B1=[00、5];th0=arx2th(A1,B1,1,1);%a1=[1,0、8

];b2=[0,0、5] y1=idsim([ue],th0); %初始模型仿真; A2=[1

0、6];B2=[00、3];th0=arx2th(A2,B2,1,1);%a2=[1,0、8

];b2=[0,0、5] y2=idsim([ue],th0); %采样300次后模型仿真; fork=1:300y(k,1)=((300-k)/300)*y1(k,1)+(k/300)*y2(k,1);%时变模型end10/18/2024

z=[yu];na=1;nb=1;nk=1;%产生输入-输出信号[z]和ARX模型得阶次 [thm,yhat]=rarx(z,[nanbnk],'ff',0、97);