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

第三章参数估计理论与应用3.1参数估计的评价准则3.2基于统计分布的参数估计方法3.3基于模型的参数最小二乘估计

本章小结2/1/2023第三章参数估计理论与应用

在许多情况下，观测数据所服从的概率模型已知的，而模型的未知部分是以未知参数形式出现的。 参数估计的基础是优化理论，即被估计的参数应该在某种准则下是最优的，以及任何获得最优的估计。 非参数估计方法不假定观测数据服从某种特定的概率模型。例如，频域上的谱估计与谱线拟合就是典型的非参数估计方法。观测到的状态状态控制x(t)y(t)u(t)v(t)w(t)观量噪声设备噪声设备（模型结构已知、参数未知）测量装置图3-1

系统辨识中的参数估计问题2/1/2023第三章参数估计理论与应用3.1参数估计的评价准则

参数估计是通过样本去估计总体的某些数字特征或统计量。任何一个统计量都可作为参数的估计量，但其效果的优劣有所差别。3.1.1无偏性、有效性与相容性

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

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

则称是θ的无偏估计量。2/1/2023第三章参数估计理论与应用

例3-1

样本均值是总体数学期望的无偏估计。 设x1,…,xN

是随机过程{xk}的N个独立观测样本，如果参数θ是总体的数学期望E[x]，即用样本的均值作为θ的估计量，对该估计量取期望值，有

一个无偏估计量在多次估计中将不会产生系统偏差，但并不意味着有偏估计就不好。如果一个有偏估计是渐进无偏的，即2/1/2023第三章参数估计理论与应用那么它仍然有可能是一个好的估计。 考虑实随机过程{xk}的相关函数的两种估计量：

假定数据{xk}是独立观测的，容易验证

式中，Rx(τ)=E[xk+τ

xk]

是随机数据{xk}的相关函数。 以上二式表明，估计量1(τ)是无偏的，而2(τ)则是有偏的。但是，2(τ)是渐进无偏的，即2/1/2023第三章参数估计理论与应用渐进无偏估计量2(τ)是半正定的，而无偏估计量1(τ)却不一定是半正定的，故2(τ)的使用场合较多。

（2）有效性 如果1

和2

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

1的值比2

的值更密集地聚集在真值θ的附近。通常将方差（或协方差阵）在所有的无偏估计量中达到最小的

称为有效估计量。

例3-2

设x1,…,xN

是N个独立观测样本，若被估计参数2/1/2023第三章参数估计理论与应用θ=E[x]，则对任何满足都是θ的无偏估计量。利用不等式可得在估计总体的数学期望时，简单的算术平均比加权平均好。

（3）一致性

估计量的精度是与样本的容量N

有关系的。一般说来，总是认为N

越大估计的效果应该越好。如果记依赖样本容量N的估计为N

，当满足2/1/2023第三章参数估计理论与应用则称N是θ的一致性估计量，或相容估计。

例3-3

设总体x具有均匀分布，分布密度为其中，θ1

和θ2

是未知参数。

总体样本的均值和二阶矩分别为（严格按定义计算）解得2/1/2023第三章参数估计理论与应用 按矩的估计方法，用独立样本的均值和独立样本的二阶矩，分别作为总体均值和二阶矩的估计量，就有 下面说明1

和2

分别是θ1

和θ2

的相容估计。 设y1,…,yN

是具有同分布的独立观测样本，根据大数定律，有令y=x2,就有2/1/2023第三章参数估计理论与应用于是3.1.2Fisher信息和Cramer-Rao不等式

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

Fisher信息Fisher信息用J（θ）表示，定义为（3.1.1）2/1/2023第三章参数估计理论与应用 当考虑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下限），即上式等号成立的充要条件是2/1/2023第三章参数估计理论与应用其中,函数K(θ)>0，并与样本向量X无关。 当为有偏估计量时，Cramer-Rao

不等式为

（3.1.3）

式中η(θ)为估计偏差，即η(θ)=E[]-θ，并假定b(θ)是可微分的。 对于多个参数的情况，记θ={θ1，…，θp}，则用矩阵J(θ)表示Fisher信息，其元素Jij(θ)定义为（3.1.4）2/1/2023第三章参数估计理论与应用且Cramer-Rao不等式变为矩阵不等式：（3.1.5）

上式表示无偏估计量的协方差矩阵cov()与逆Fisher信息阵之差是一半正定矩阵。

Fisher信息是描述从观测数据中得到的θ的“信息”测度，它给出利用观测数据估计参数θ的方差下界。但是，满足这一下界的估计量有的时候可能不存在。3.2基于统计分布的参数估计方法 参数估计量的优劣取决于所采用的评价准则（或代价函数）和估计算法。现在介绍已知总体统计分布的两种最有效的参数估计方法：Bayes

估计和最大似然估计。2/1/2023第三章参数估计理论与应用3.2.1Bayes

估计

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

}进行参数估计，通常是采用代价函数的期望值作为评价参数估计量效果的测度，并称之为风险函数。使风险函数最小的参数估计叫做Bayes

估计；基于二次型风险函数最小的估计称为最小均方误差（minimummeansquareerror,MMSE）估计。二次型风险函数定义为（3.2.1） 根据条件概率公式，有2/1/2023第三章参数估计理论与应用其中，p(θ

|x1,…,xN

)是给定N个观测样本X={x1,…,xN

}条件下θ

的后验分布密度函数。于是，式（3.2.1）可以写成（3.2.2） 为使风险函数RM

MSE

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

)是非负的，因此，∂RM

MSE

/∂=0,等价于上式中[·]=0。故有2/1/2023第三章参数估计理论与应用（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维。2/1/2023第三章参数估计理论与应用其中，测量误差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）来求解。2/1/2023第三章参数估计理论与应用 对式（1）两边取均值，得到（4）

将式（1）和（3）代入有关定义式，得（5）（6）（7）2/1/2023第三章参数估计理论与应用

（i）当e与x

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

（ii）当e

与x

相关，只需注意Cxe

≠0即可。 这个问题留给读者解决。请构造一组数据，在Matlab

平台上仿真这两种的估计结果。3.2.2最大似然估计

最大似然估计（maximumlikelihoodestimate,ML估计）的基本思路是：在给定参数θ条件下，将观测样本xK2/1/2023第三章参数估计理论与应用联合条件概率密度函数p(x|θ)视为真实参数θ

的函数，即似然函数L(x，θ)

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

}，求出使L(x，θ)达到最大化的参数作为θ={θ1，…，θp}的估计值。在数学上，通常将未知参数θ的最大似然估计量记为式中Θ是参数θ的值域。故ML估计量ML就是p(x|θ)的全局极大点。 由于对数函数是严格单调的，故L(x，θ)的极大点与ln

L(x，θ)的极大点是一致的。通常，将ln

L(x，θ)称为对数似然函数。于是，ML估计量ML可由（3.2.4）2/1/2023第三章参数估计理论与应用确定。如果x1,…,xN

是N个独立的观测样本，则对数似然函数可写作（3.2.5）

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

例3-5

设样本x={x1,…,xN

}服从高斯分布N（m，σ），则其对数似然函数为2/1/2023第三章参数估计理论与应用分别求lnL

关于m和σ2

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

是无偏的，而方差的ML估计

ML是有偏的。但若将ML

·N/(N-1)作为新的估计量，则该估计是无偏的。2/1/2023第三章参数估计理论与应用

计算L(x，θ)的相对于m

的二阶偏导数，有由式（3.1.1）得Fisher信息：Cramer-Rao不等式为等号成立的充要条件是 事实上，我们有2/1/2023第三章参数估计理论与应用因此，只要取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

二元阵测向系统的几何关系2/1/2023第三章参数估计理论与应用

解

设两传感器的零均值接收过程可分别表示为其中，si

为单频平面波信号，wi(i=1，2)为零均值高斯噪声，二者互相独立。

如果采用图3-3所示的时延补偿方法，则单频平面波信号的归一化声程补偿（或指向）向量v

在所考虑的二元阵中可表示为

下面，我们来推导信号的协方差矩阵和噪声的协方差矩阵，以便于求出观测样本的似v*图3-3声程补偿系统x2x11exp(-jωnτ)∑2/1/2023第三章参数估计理论与应用函数。记输入信号和输入噪声的傅立叶系数为设信号和噪声的功率谱分别为S(ωn)和N(ωn),那么，由公式（1.4.6）[pp.26,ωn=2πn/T）]

信号和噪声的协方差矩阵可分别表示为（3.2.7）2/1/2023第三章参数估计理论与应用于是，观测样本的似然函数可表示为

（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估计就是选择τ，使ln

p(X|τ)最大,也即使式（3.2.8）的指数函数（3.2.9）2/1/2023第三章参数估计理论与应用最大。下面，我们从式（3.2.9）出发，推导时延参数τ

的最大似然估计的等效形式。为此，首先引进下列求逆公式（3.2.10）式中，A为n×n非奇异矩阵；g为n×1列向量。证明留给请读者课外练习【利用恒等式g(1+gHA-1g)=(A+gHg)A-1g）】。 利用求逆公式，可知<g=[1exp(jωnτ)]>2/1/2023第三章参数估计理论与应用将上式代入式（3.2.9），略去与τ无关的量T/N(ωn)。因此，选择τ使式（3.2.9）最大，等价于使下式（3.2.11）

最大。现引入记号在此将X(ωn，τ)视为某时间函数x(t,τ)在时间（t-T，t）内 的傅立叶系数。将上述替换量代入式（3.2.11）后，再应用 周期函数的Parseval

公式，就有2/1/2023第三章参数估计理论与应用略去无关紧要的常数项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ωτ)2/1/2023第三章参数估计理论与应用 二元阵最大似然测向系统与二元阵似然比检测系统具有完全相同的结构。这是因为：在H1

情况下，p(X|τ)等价于p(X

|H1)，后者也可看作是时延参数τ的函数；而在H0

情况下，p(X

|H0)与τ无关。因此，选取τ

使似然函数最大，也就是使似然比

p(X|H1)/

p(X|H0)最大。由此可见，检测问题与参数估计问题是密切相关。 顺便指出，可用测向测距近似公式（3.2.13） 构成最大似然联合测向测距系统。其中，di

表示第i

个传感器与“基准”传感器位置的间距；D

表示目标与“基准”传感器位置之间的距离。2/1/2023第三章参数估计理论与应用3.3基于模型的参数最小二乘估计

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

在一般的最小二乘问题中，线性系统的参数化模型可以表示为(3.3.1)

其中，u=[u1,…,up]T

是模型的输入向量，f1,…,fn

是u的已知函数，也可以是未知输入的观测数据；θ1,…,θn

是待估计2/1/2023第三章参数估计理论与应用的参数，又称为回归系数；y是系统的输出。

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

yi]

或[fk

(ui)

yi],i=1,2,…,N,k=1,2,…,n；N≥n}以构成训练数据。将数据对代入方程(3.3.1)，可以获得一组线性方程： 用矩阵表示方法，将上式写成更简洁的形式，即2/1/2023第三章参数估计理论与应用

(3.3.2)其中

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

(3.3.3)2/1/2023第三章参数估计理论与应用

参数θ的最小二乘估计LS

，就是使目标函数

(3.3.4)

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

对θ求导数，有

J

极小化的条件是一般均假设ΦTΦ非奇异，于是，LS有唯一的解：2/1/2023第三章参数估计理论与应用

(3.3.5)

式中Φ+表示Φ的伪逆。 上述表示误差向量对整体平方误差有相同权重。可以进一步扩展，令每个误差项有不同的权重。设W为所需的权值矩阵，它是对称和正定的，则加权的目标函数为(3.3.6)

按上述求极小值的方法，可得加权的最小二乘估计量：(3.3.7)显然，当W选为单位矩阵时，WLS

=LS。

例3-7

考虑最简单的一维线性模型（静态的），即只有一个控制变量u的情形，这时模型的形式是2/1/2023第三章参数估计理论与应用求未知参数θ0

和θ1的LS估计量。

解实际过程输出是模型的输出加上一随机误差项，即观测数据对[ui,yi]的结构应为式中，ei

称为模型的残差或观测噪声，一般认为是零均值、相互独立的随机序列，并具有相同的方差σ2。将上式写成矩阵形式：

根据式(3.3.5)，可得LS估计量：2/1/2023第三章参数估计理论与应用 如果进一步假定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。2/1/2023第三章参数估计理论与应用

例5-8

（可线性化的非线性静态模型——曲线回归）假设有一个非线性模型的输出为其中，x1，x2

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

解

上式两边经简单的代数运算，再同时取自然对数，可转化为一个线性模型：这说明变换后的输出ln(y-1-1)

可以显式地表达为以

lnx1和

x2为输入、以lna，b和

c

为参数的线性模型。因此，就可以按变换后的线性模型用最小二乘法来估计变换后的未知参数，然后，再根据变换后的估计参数计算出原参数。2/1/2023第三章参数估计理论与应用 判定输入x-输出y之间的关系能否用一个线性模型来描述的标准，通常用互相关系数的大小来衡量：(3.3.8)ρxy

的绝对值越大，表示变量之间的线性关系越密切，因而线性回归的效果就越好。

例3-9

设某一结构参数n，m和d已知的离散线性系统，其差分方程的形式为：

(3.3.9)2/1/2023第三章参数估计理论与应用其中，e(k)为噪声，φ(k)为输入-输出观测向量，θ为未知参数向量，且要求根据N次数据对{[y(i),u(i)],i=1,2,…,N；N≥n+m+1}来估计对未知参数θ。

解将式(3.3.9)改写成矩阵形式，得到将数据写成下标形式，就有这样，未知参数向量θ可按式(3.3.5)进行估计。 2/1/2023第三章参数估计理论与应用

考虑如下单输入-单输出系统：用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,[na

nbnk],'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。2/1/2023第三章参数估计理论与应用 应当指出，要求观测数据容量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）的频率成分，以保证充分激励受控对象的2/1/2023第三章参数估计理论与应用所有振型，从而使观测数据载有动态系统的主要信息。LS估计在满足持续激励条件时，是渐进无偏的，也称为估计的一致性。 式(3.3.9)，也称为CAR模型（即受控的AR模型），可以写成更简洁的形式

(3.3.10)式中q

表示时间算子，d为整数，表示系统的滞后量；A(·),B(·)分别为q－1

的降次幂多项式。

CAR模型满足一致估计（或相容估计）的条件为：

（1）{e(k)}是白噪声序列； （2）{u(k)}的均值和协方差有界；且满足（m+1）阶持续激励条件（或正定条件）：2/1/2023第三章参数估计理论与应用

（3）u(k)

和e(k)

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

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

来描述。通常，当n

或（n,m）增大时J(n)

或J(n,m)就会减小；而当n(或n,m)大于模型的真实阶次

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

均是非奇异方阵，则

(3.3.11)

下面介绍最小二乘递推算法。为了简化符号，以下推导均用代替LS。 设ΦN是时刻N为止的观测数据，N+1时刻θ的LS估计为2/1/2023第三章参数估计理论与应用式中[参见式(3.3.2)和(3.3.9)]于是有(3.3.12)

令PN=[ΦN

TΦN]-1

，由求逆公式(3.3.11)知

(3.3.13)

定义增益向量为2/1/2023第三章参数估计理论与应用将式(3.3.13)和(3.3.14)代入(3.3.12)，得到(3.3.15)

上式表明，新的估计量N+1

等于前一时刻的估计量N

与修正项KN+1(yN+1-φTN+1N)

之和，这是一切递推公式的共同特征。如果令代表基于前一时刻的估计量N

对N+1时刻的预测。那么，递推估计提供的新息2/1/2023第三章参数估计理论与应用就是预测误差或拟合误差。因此，修正量的大小与新息成正比，而各校正分量的权，由增益向量决定。 在启动上述递推算法时，必须知道初值0

和P0

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

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

ΦNTΦN

通常是正定的。显然，当N

趋于无穷大，ΦNTΦN

/N

接近于非奇异的常数阵。于是，有 可见，式(3.3.16)中自适应增益向量KN+1

随着每次迭代而递减，这意味着递推运算过程将逐渐收敛于参数空间的最优点。事实上，在白噪声或低噪声条件下，递推最小二乘2/1/2023第三章参数估计理论与应用估计是一种简便而又有效的算法。这种递推算法在递推过程中虽然没有保存全部先前的数据，但所有先前数据的影响却一直在起作用，故称为无限增长记忆的递推最小二乘估计。3.3.3卡尔曼滤波器的递推算法（状态估计） 与参数估计器不同，卡尔曼滤波器主要是解决如何从被噪声污染的观测数据中估计出已知动态系统模型的状态，而不是动态系统模型的未知参数。然而，仅从算法上看，这二者是非常相似的。为了便于比较二者的异同，我们在此不加证明地列出卡尔曼滤波器算法。 设离散定常系统的状态方程和观测方程分别为2/1/2023第三章参数估计理论与应用式中，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，称为平滑问题。 下面主要介绍滤波问题。2/1/2023第三章参数估计理论与应用滤波算法：预测算法：滤波增益：滤波误差的协方差：预测误差的协方差：初始条件：2/1/2023第三章参数估计理论与应用

例3-10

考虑某一动态系统表示的空间导航问题，其中加速度为白噪声。被白噪声污染之后观测其位置。因此，过程的状态方程为（1）式中，X(t)=[x1(t),x2(t)]T,表示过程的位置和速度；而w(t)是一均值为0、方差为1的高斯白噪声。观测方程为（2）其中v(t)是一均值为0、方差为10的高斯白噪声。

解

首先把连续的状态方程离散化（参见Matlab中c2d函数）2/1/2023第三章参数估计理论与应用不妨设采样间隔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的仿真计算。2/1/2023第三章参数估计理论与应用3.3.4限定记忆的递推最小二乘估计

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

设置为P0，重新进行参数估计，这是因为LS估计能快速地收敛到当前的最优参数。处理这一问题的另一种方法是对过去的数据引入带遗忘因子λ，逐渐削弱它们在参数估计中的作用。为此，采用加权的目标函数：2/1/2023第三章参数估计理论与应用其中

(3.3.17)

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

(3.3.18)

每当取得一个新的数据后，就对沟渠以前的加权矩阵乘以λ，于是N+1时刻的估计量为2/1/2023第三章参数估计理论与应用与前面略有不同，在此令PN=[ΦN

TWNΦN]-1

，则由求逆公式(3.3.11)，可知于是，按前面推导出式(3.3.15)思路，可得式中

遗忘因子λ将老的数据逐渐从“记忆”中去掉，因此这种使用数据信息的方式也叫做“渐消记忆”法，相应的算法称为带遗忘因子的递推最小二乘法，即2/1/2023第三章参数估计理论与应用关于λ的选取通常由经验或实验确定，一般范围为0.95≤λ≤0.99 λ

取得愈小，最新数据的权重就愈大，也就更适合于跟踪大的时变参数，但与此同时，估计器也可能会发生较大的波动 从而加大估计的误差。 例5-10

考虑模型采样300次后变为2/1/2023第三章参数估计理论与应用

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

解用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);%时变模型

end2/1/2023第三章参数估计理论与应用

z=[yu]；na=1;nb=1;nk=1;

%产生输入-输出信号[z]和ARX模型的阶次

[thm,yhat]=rarx(z,[na

nbnk],'ff',0.97);

%带遗忘因子0.97的ARX参数辨识

plot(thm(:,1),‘-’);grid

%作图，时变系统的a(k)输出曲线

holdon

plot(thm(:,2),':')

%作图，时变系统的b(k)输出曲线图3-5带遗忘因子λ的时变系统参数辨识结果2/1/2023第三章参数估计理论与应用 从仿真结果（图3-5）可以看出，采用带遗忘因子λ的系统辨识算法能较好地跟踪系统时变参数的变化。3.3.5广义最小二乘估计 上述推导中，假定模型(3.3.10)中的噪声项e(t)为白噪声。如果e(t)为有色噪声，则被测对象就应当的用CARMA模型（受控的ARMA）来表示，即与CAR模型（基本最小二乘估计的数学模型）比较，此处C(q-1)≠1，而有色噪声e(t)

可以认为是由白噪声ξ(t)通过成型滤波器而产生的输出（以下写成离散数据的形式）：2/1/2023第三章参数估计理论与应用 下面介绍广义最小二乘法，用以改善CARMA模型参数估计的统计特性。 将CARMA模型写成常规的最小二乘结构，即(3.3.20)

因为ξ(k)是不可测量的，所以在初始估计时只能用ξ(k)估计值代替φ(k)中有关ξ(k)的分量。按基本最小二乘法估计出ŷ，然后，计算ξ(k)的估计量(3.3.21)再构造2/1/2023第三章参数估计理论与应用这样，就可以按前面介绍的递推公式(3.3.22)

对CARMA模型进行参数估计。估计的初始值可按下式选取其中，σ2»1。在实际应用中，ξ(k)的估计有两种形式，一种是预测的形式：2/1/2023第三章参数估计理论与应用为了避免这种问题的产生，最常采用的方法是对PN

采 用U-D分解法,下面就介绍这种方法。3.3.6改进数值稳定性的U-D分解法 为了保证在递推过程中始终保持PN

的非负定性，可把矩阵PN

进行U-D分解，即【在此，将PN

写成P(N)】其中，U是对角元全为1的上三角矩阵，D为对角矩阵。这种分解只需要实时修正U，而无需修正P，因而能保证P的非负定性。下面，我们不加证明地将U-D分解基本算法的计算步骤总结如下：

2/1/2023第三章参数估计理论与应用

(1)

设置初值θˆ(0)，U(0)和D(0)，输入初始数据；

(2)

采样当前的输入和输出（假设数据的行数为p）；

(3)

按下式计算f(N)和g(N)

(4)

定义由此可得2/1/2023第三章参数估计理论与应用按上式计算βi(N)； （5）计算

（6）计算参数估计 （7）按下式计算di(N)式中，di是D的对角元素； （8）引入上三角阵