參數(shù)估計理論與應(yīng)用(第三章 )

第三章參數(shù)估計理論與應(yīng)用3.1參數(shù)估計的評價準(zhǔn)則3.2基于統(tǒng)計分布的參數(shù)估計方法3.3基于模型的參數(shù)最小二乘估計

本章小結(jié)2/1/2023第三章參數(shù)估計理論與應(yīng)用

在許多情況下，觀測數(shù)據(jù)所服從的概率模型已知的，而模型的未知部分是以未知參數(shù)形式出現(xiàn)的。 參數(shù)估計的基礎(chǔ)是優(yōu)化理論，即被估計的參數(shù)應(yīng)該在某種準(zhǔn)則下是最優(yōu)的，以及任何獲得最優(yōu)的估計。 非參數(shù)估計方法不假定觀測數(shù)據(jù)服從某種特定的概率模型。例如，頻域上的譜估計與譜線擬合就是典型的非參數(shù)估計方法。觀測到的狀態(tài)狀態(tài)控制x(t)y(t)u(t)v(t)w(t)觀量噪聲設(shè)備噪聲設(shè)備（模型結(jié)構(gòu)已知、參數(shù)未知）測量裝置圖3-1

系統(tǒng)辨識中的參數(shù)估計問題2/1/2023第三章參數(shù)估計理論與應(yīng)用3.1參數(shù)估計的評價準(zhǔn)則

參數(shù)估計是通過樣本去估計總體的某些數(shù)字特征或統(tǒng)計量。任何一個統(tǒng)計量都可作為參數(shù)的估計量，但其效果的優(yōu)劣有所差別。3.1.1無偏性、有效性與相容性

（1）無偏性設(shè)樣本的總體分布密度函數(shù)為p(x;θ),θ是未知參數(shù)。從總體中抽取容量為N的樣本x={x1,…,xN

},用樣本的估計量來估計θ，如果希望多次估計中，平均的估計值沒有偏差，即

則稱是θ的無偏估計量。2/1/2023第三章參數(shù)估計理論與應(yīng)用

例3-1

樣本均值是總體數(shù)學(xué)期望的無偏估計。 設(shè)x1,…,xN

是隨機過程{xk}的N個獨立觀測樣本，如果參數(shù)θ是總體的數(shù)學(xué)期望E[x]，即用樣本的均值作為θ的估計量，對該估計量取期望值，有

一個無偏估計量在多次估計中將不會產(chǎn)生系統(tǒng)偏差，但并不意味著有偏估計就不好。如果一個有偏估計是漸進(jìn)無偏的，即2/1/2023第三章參數(shù)估計理論與應(yīng)用那么它仍然有可能是一個好的估計。 考慮實隨機過程{xk}的相關(guān)函數(shù)的兩種估計量：

假定數(shù)據(jù){xk}是獨立觀測的，容易驗證

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

xk]

是隨機數(shù)據(jù){xk}的相關(guān)函數(shù)。 以上二式表明，估計量1(τ)是無偏的，而2(τ)則是有偏的。但是，2(τ)是漸進(jìn)無偏的，即2/1/2023第三章參數(shù)估計理論與應(yīng)用漸進(jìn)無偏估計量2(τ)是半正定的，而無偏估計量1(τ)卻不一定是半正定的，故2(τ)的使用場合較多。

（2）有效性 如果1

和2

是兩個根據(jù)N個獨立觀測樣本得到的無偏估計量，無疑地，對θ的平均偏差較小是選擇的標(biāo)準(zhǔn)之一。例如，如果則

1的值比2

的值更密集地聚集在真值θ的附近。通常將方差（或協(xié)方差陣）在所有的無偏估計量中達(dá)到最小的

稱為有效估計量。

例3-2

設(shè)x1,…,xN

是N個獨立觀測樣本，若被估計參數(shù)2/1/2023第三章參數(shù)估計理論與應(yīng)用θ=E[x]，則對任何滿足都是θ的無偏估計量。利用不等式可得在估計總體的數(shù)學(xué)期望時，簡單的算術(shù)平均比加權(quán)平均好。

（3）一致性

估計量的精度是與樣本的容量N

有關(guān)系的。一般說來，總是認(rèn)為N

越大估計的效果應(yīng)該越好。如果記依賴樣本容量N的估計為N

，當(dāng)滿足2/1/2023第三章參數(shù)估計理論與應(yīng)用則稱N是θ的一致性估計量，或相容估計。

例3-3

設(shè)總體x具有均勻分布，分布密度為其中，θ1

和θ2

是未知參數(shù)。

總體樣本的均值和二階矩分別為（嚴(yán)格按定義計算）解得2/1/2023第三章參數(shù)估計理論與應(yīng)用 按矩的估計方法，用獨立樣本的均值和獨立樣本的二階矩，分別作為總體均值和二階矩的估計量，就有 下面說明1

和2

分別是θ1

和θ2

的相容估計。 設(shè)y1,…,yN

是具有同分布的獨立觀測樣本，根據(jù)大數(shù)定律，有令y=x2,就有2/1/2023第三章參數(shù)估計理論與應(yīng)用于是3.1.2Fisher信息和Cramer-Rao不等式

通常希望獲得有效的參數(shù)估計量。但是，由于不存在導(dǎo)致最小方差無偏估計量的最佳算法，所以通常采用參數(shù)無偏估計的Cramer-Rao下限(或CR下界),作為評價參數(shù)估計性能的測度。為了簡潔敘述這一的評價測度，先定義一個重要的概念。

Fisher信息Fisher信息用J（θ）表示，定義為（3.1.1）2/1/2023第三章參數(shù)估計理論與應(yīng)用 當(dāng)考慮N個觀測樣本X={x1,…,xN

},此時，聯(lián)合條件分布密度函數(shù)可表示為 將式（3.1.1）中的p(x|θ)改為p(X|θ)就可給出N個樣本變量X的Fisher信息的表達(dá)式。

定理（Cramer-Rao不等式）設(shè)觀測樣本X={x1,…,xN

},若參數(shù)估計是真實參數(shù)θ

的無偏估計，并且條件分布密度函數(shù)的p(X|θ)對參數(shù)θ

的一、二階偏導(dǎo)數(shù)存在，則有（3.1.2） 參數(shù)的方差所能達(dá)到的下限（稱為CR下限），即上式等號成立的充要條件是2/1/2023第三章參數(shù)估計理論與應(yīng)用其中,函數(shù)K(θ)>0，并與樣本向量X無關(guān)。 當(dāng)為有偏估計量時，Cramer-Rao

不等式為

（3.1.3）

式中η(θ)為估計偏差，即η(θ)=E[]-θ，并假定b(θ)是可微分的。 對于多個參數(shù)的情況，記θ={θ1，…，θp}，則用矩陣J(θ)表示Fisher信息，其元素Jij(θ)定義為（3.1.4）2/1/2023第三章參數(shù)估計理論與應(yīng)用且Cramer-Rao不等式變?yōu)榫仃嚥坏仁剑海?.1.5）

上式表示無偏估計量的協(xié)方差矩陣cov()與逆Fisher信息陣之差是一半正定矩陣。

Fisher信息是描述從觀測數(shù)據(jù)中得到的θ的“信息”測度，它給出利用觀測數(shù)據(jù)估計參數(shù)θ的方差下界。但是，滿足這一下界的估計量有的時候可能不存在。3.2基于統(tǒng)計分布的參數(shù)估計方法 參數(shù)估計量的優(yōu)劣取決于所采用的評價準(zhǔn)則（或代價函數(shù)）和估計算法?，F(xiàn)在介紹已知總體統(tǒng)計分布的兩種最有效的參數(shù)估計方法：Bayes

估計和最大似然估計。2/1/2023第三章參數(shù)估計理論與應(yīng)用3.2.1Bayes

估計

在參數(shù)估計中，估計誤差θ-通常不為零。因此，除了采用前面介紹的無偏、有效和相容估計作為評價準(zhǔn)則外，還可以利用估計誤差的變化范圍作為參數(shù)估計的測度，這種測度叫做代價函數(shù)，用符號C(,θ)表示。常用的代價函數(shù)有絕對型、二次型和均勻型三種。OOO?/2?/2絕對型二次型均勻型2/1/2023第三章參數(shù)估計理論與應(yīng)用 本節(jié)僅介紹最常用的二次型代價函數(shù)，即 當(dāng)總體的分布密度函數(shù)p(X|θ)已知時，利用X={x1,…,xN

}進(jìn)行參數(shù)估計，通常是采用代價函數(shù)的期望值作為評價參數(shù)估計量效果的測度，并稱之為風(fēng)險函數(shù)。使風(fēng)險函數(shù)最小的參數(shù)估計叫做Bayes

估計；基于二次型風(fēng)險函數(shù)最小的估計稱為最小均方誤差（minimummeansquareerror,MMSE）估計。二次型風(fēng)險函數(shù)定義為（3.2.1） 根據(jù)條件概率公式，有2/1/2023第三章參數(shù)估計理論與應(yīng)用其中，p(θ

|x1,…,xN

)是給定N個觀測樣本X={x1,…,xN

}條件下θ

的后驗分布密度函數(shù)。于是，式（3.2.1）可以寫成（3.2.2） 為使風(fēng)險函數(shù)RM

MSE

最小，對上式取的偏導(dǎo)，并令其結(jié)果為零，便得到由于p(x1,…,xN

)是非負(fù)的，因此，?RM

MSE

/?=0,等價于上式中[·]=0。故有2/1/2023第三章參數(shù)估計理論與應(yīng)用（3.2.3）

注意，在式（3.2.3）中，利用了以下事實： 由此可得出重要的結(jié)論：未知參數(shù)θ

的MMSE估計是給定樣本X條件下θ的條件均值。

例3-4

某一隨機參量x服從高斯N(mx,Cx)分布，用儀器可測量其線性組合y

，即（1）式中，y－N維，k－N×M維，x－M維，e－N維。2/1/2023第三章參數(shù)估計理論與應(yīng)用其中，測量誤差e服從高斯N(0,Ce)分布；k為給定的常數(shù)陣。假設(shè)

(ⅰ)e與x

獨立；

(ⅱ)e與x

相關(guān)，互協(xié)方差函數(shù)為Cxe

。 試分別求出兩種情況下的MMSE估計x?(y)和估計誤差x

(y)的協(xié)方差Rx(y)。

解如果將x

看作未知參數(shù)，那么，根據(jù)上面討論,x的MMSE估計是給定觀測樣本{y1,…,yN}時x的條件均值。因此，可利用公式（1.4.16）和（1.4.17）[pp.29]

（2）（3）來求解。2/1/2023第三章參數(shù)估計理論與應(yīng)用 對式（1）兩邊取均值，得到（4）

將式（1）和（3）代入有關(guān)定義式，得（5）（6）（7）2/1/2023第三章參數(shù)估計理論與應(yīng)用

（i）當(dāng)e與x

互相獨立，Cxe=0。將式（4）~（7）代入式（2）和（3），得到x?(y)的估計及協(xié)方差Rx(y)

（ii）當(dāng)e

與x

相關(guān)，只需注意Cxe

≠0即可。 這個問題留給讀者解決。請構(gòu)造一組數(shù)據(jù)，在Matlab

平臺上仿真這兩種的估計結(jié)果。3.2.2最大似然估計

最大似然估計（maximumlikelihoodestimate,ML估計）的基本思路是：在給定參數(shù)θ條件下，將觀測樣本xK2/1/2023第三章參數(shù)估計理論與應(yīng)用聯(lián)合條件概率密度函數(shù)p(x|θ)視為真實參數(shù)θ

的函數(shù)，即似然函數(shù)L(x，θ)

（包含未知參數(shù)θ的可能性函數(shù)），然后利用容量為N的觀測樣本x={x1,…,xN

}，求出使L(x，θ)達(dá)到最大化的參數(shù)作為θ={θ1，…，θp}的估計值。在數(shù)學(xué)上，通常將未知參數(shù)θ的最大似然估計量記為式中Θ是參數(shù)θ的值域。故ML估計量ML就是p(x|θ)的全局極大點。 由于對數(shù)函數(shù)是嚴(yán)格單調(diào)的，故L(x，θ)的極大點與ln

L(x，θ)的極大點是一致的。通常，將ln

L(x，θ)稱為對數(shù)似然函數(shù)。于是，ML估計量ML可由（3.2.4）2/1/2023第三章參數(shù)估計理論與應(yīng)用確定。如果x1,…,xN

是N個獨立的觀測樣本，則對數(shù)似然函數(shù)可寫作（3.2.5）

ML估計量ML只要能夠求出來，總是比較好的估計，它具有以下性質(zhì)：最大似然估計是有效和一致估計；對于大的N，ML估計量ML服從高斯分布，并且是無偏的，方差可達(dá)CR下界。

例3-5

設(shè)樣本x={x1,…,xN

}服從高斯分布N（m，σ），則其對數(shù)似然函數(shù)為2/1/2023第三章參數(shù)估計理論與應(yīng)用分別求lnL

關(guān)于m和σ2

的偏導(dǎo)，并令它們等于零，得到解得顯然有 可見，均值的ML估計ML

是無偏的，而方差的ML估計

ML是有偏的。但若將ML

·N/(N-1)作為新的估計量，則該估計是無偏的。2/1/2023第三章參數(shù)估計理論與應(yīng)用

計算L(x，θ)的相對于m

的二階偏導(dǎo)數(shù)，有由式（3.1.1）得Fisher信息：Cramer-Rao不等式為等號成立的充要條件是 事實上，我們有2/1/2023第三章參數(shù)估計理論與應(yīng)用因此，只要取K(m)=N/σ2，ML估計ML就可達(dá)CR下界σ2/N。這表明ML估計

ML是一有效估計量。

例3-6

（二元陣最大似然測向系統(tǒng)）設(shè)二元陣布置在x軸上，兩個基元坐標(biāo)分別為x1

和x2，如圖3-2所示。如果取x1=0，則x2=d，d為兩傳感器的位置間隔。假設(shè)信號為平面波，入射角為θ，則傳感器1相對于傳感器2的信號時延τ為（3.2.6）式中，c

為聲速。我們的問題是如何利用二元陣中兩個輸入過程的時差τ來測定目標(biāo)的方位角θ。θxx2=dx1=0圖3-2

二元陣測向系統(tǒng)的幾何關(guān)系2/1/2023第三章參數(shù)估計理論與應(yīng)用

解

設(shè)兩傳感器的零均值接收過程可分別表示為其中，si

為單頻平面波信號，wi(i=1，2)為零均值高斯噪聲，二者互相獨立。

如果采用圖3-3所示的時延補償方法，則單頻平面波信號的歸一化聲程補償（或指向）向量v

在所考慮的二元陣中可表示為

下面，我們來推導(dǎo)信號的協(xié)方差矩陣和噪聲的協(xié)方差矩陣，以便于求出觀測樣本的似v*圖3-3聲程補償系統(tǒng)x2x11exp(-jωnτ)∑2/1/2023第三章參數(shù)估計理論與應(yīng)用函數(shù)。記輸入信號和輸入噪聲的傅立葉系數(shù)為設(shè)信號和噪聲的功率譜分別為S(ωn)和N(ωn),那么，由公式（1.4.6）[pp.26,ωn=2πn/T）]

信號和噪聲的協(xié)方差矩陣可分別表示為（3.2.7）2/1/2023第三章參數(shù)估計理論與應(yīng)用于是，觀測樣本的似然函數(shù)可表示為

（3.2.8）式中，X(1)=[X1

(1)，X2(1)]T

，…，X(TW)=[X1(TW)，X2(TW)]T是傳感器的接收過程{x=[x1，x2]T}的傅立葉系數(shù)陣；

T

是過程的持續(xù)時間（采樣數(shù)據(jù)的長度），W

是接收過程的帶寬。 容易驗證，行列式|

Cw

+Cs

|與時延τ無關(guān)。于是，ML估計就是選擇τ，使ln

p(X|τ)最大,也即使式（3.2.8）的指數(shù)函數(shù)（3.2.9）2/1/2023第三章參數(shù)估計理論與應(yīng)用最大。下面，我們從式（3.2.9）出發(fā)，推導(dǎo)時延參數(shù)τ

的最大似然估計的等效形式。為此，首先引進(jìn)下列求逆公式（3.2.10）式中，A為n×n非奇異矩陣；g為n×1列向量。證明留給請讀者課外練習(xí)【利用恒等式g(1+gHA-1g)=(A+gHg)A-1g）】。 利用求逆公式，可知<g=[1exp(jωnτ)]>2/1/2023第三章參數(shù)估計理論與應(yīng)用將上式代入式（3.2.9），略去與τ無關(guān)的量T/N(ωn)。因此，選擇τ使式（3.2.9）最大，等價于使下式（3.2.11）

最大?，F(xiàn)引入記號在此將X(ωn，τ)視為某時間函數(shù)x(t,τ)在時間（t-T，t）內(nèi) 的傅立葉系數(shù)。將上述替換量代入式（3.2.11）后，再應(yīng)用 周期函數(shù)的Parseval

公式，就有2/1/2023第三章參數(shù)估計理論與應(yīng)用略去無關(guān)緊要的常數(shù)項1/2，計算z(x,τ)的結(jié)構(gòu)如圖3-4所示。調(diào)節(jié)時延τ，使輸出z(x,τ)達(dá)到最大，相應(yīng)的時延就是 真實時延的ML估計ML。 根據(jù)ML估計的傳遞性，由式（3.2.6）可得真實方位的ML估計

（3.2.12）

xH0(t)z(x,τ)∑x1(t)x2(t)H0(ω)(·)2圖3-4

二元陣最大似然測向系統(tǒng)exp(-jωτ)2/1/2023第三章參數(shù)估計理論與應(yīng)用 二元陣最大似然測向系統(tǒng)與二元陣似然比檢測系統(tǒng)具有完全相同的結(jié)構(gòu)。這是因為：在H1

情況下，p(X|τ)等價于p(X

|H1)，后者也可看作是時延參數(shù)τ的函數(shù)；而在H0

情況下，p(X

|H0)與τ無關(guān)。因此，選取τ

使似然函數(shù)最大，也就是使似然比

p(X|H1)/

p(X|H0)最大。由此可見，檢測問題與參數(shù)估計問題是密切相關(guān)。 順便指出，可用測向測距近似公式（3.2.13） 構(gòu)成最大似然聯(lián)合測向測距系統(tǒng)。其中，di

表示第i

個傳感器與“基準(zhǔn)”傳感器位置的間距；D

表示目標(biāo)與“基準(zhǔn)”傳感器位置之間的距離。2/1/2023第三章參數(shù)估計理論與應(yīng)用3.3基于模型的參數(shù)最小二乘估計

最小二乘法（Leastsquaremethod，LS）是一種不需要任何先驗知識的參數(shù)估計方法。在被測系統(tǒng)的靜態(tài)（穩(wěn)態(tài)）模型和動態(tài)模型的參數(shù)辨識中，最小二乘法是最常用的參數(shù)估計方法，在測控技術(shù)領(lǐng)域獲得了廣泛的應(yīng)用。3.3.1最小二乘估計器及其統(tǒng)計特性

在一般的最小二乘問題中，線性系統(tǒng)的參數(shù)化模型可以表示為(3.3.1)

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

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

是u的已知函數(shù)，也可以是未知輸入的觀測數(shù)據(jù)；θ1,…,θn

是待估計2/1/2023第三章參數(shù)估計理論與應(yīng)用的參數(shù)，又稱為回歸系數(shù)；y是系統(tǒng)的輸出。

當(dāng)f1,…,fn是u的穩(wěn)態(tài)響應(yīng)狀態(tài)或是實測的確定性變量，且y是系統(tǒng)的穩(wěn)態(tài)輸出，則稱式(3.3.1)是描述線性系統(tǒng)的靜態(tài)模型；當(dāng)y是u的動態(tài)響應(yīng)或瞬態(tài)觀測數(shù)據(jù)，那末式(3.3.1)就是描述線性系統(tǒng)的動態(tài)模型。 為了估計未知參數(shù)θi,必須做實驗來獲得數(shù)據(jù)對{[ui

yi]

或[fk

(ui)

yi],i=1,2,…,N,k=1,2,…,n；N≥n}以構(gòu)成訓(xùn)練數(shù)據(jù)。將數(shù)據(jù)對代入方程(3.3.1)，可以獲得一組線性方程： 用矩陣表示方法，將上式寫成更簡潔的形式，即2/1/2023第三章參數(shù)估計理論與應(yīng)用

(3.3.2)其中

為了唯一地識別出未知參數(shù)，通常要求N>n，即數(shù)據(jù)對的數(shù)目多于擬合參數(shù)的數(shù)目。滿足所有N個方程的精確解是不可能的，因為觀測數(shù)據(jù)難免受到噪聲的污染，或者描述系統(tǒng)的參數(shù)化數(shù)學(xué)模型不夠精確。故必須考慮隨機噪聲或建模誤差，在方程(3.3.2)中引入隨機誤差向量e，得到

(3.3.3)2/1/2023第三章參數(shù)估計理論與應(yīng)用

參數(shù)θ的最小二乘估計LS

，就是使目標(biāo)函數(shù)

(3.3.4)

達(dá)到最小值的參數(shù)估計。為此，通常都采用求極值的方法。 將式(3.3.4)展開后，得到

對θ求導(dǎo)數(shù)，有

J

極小化的條件是一般均假設(shè)ΦTΦ非奇異，于是，LS有唯一的解：2/1/2023第三章參數(shù)估計理論與應(yīng)用

(3.3.5)

式中Φ+表示Φ的偽逆。 上述表示誤差向量對整體平方誤差有相同權(quán)重?？梢赃M(jìn)一步擴展，令每個誤差項有不同的權(quán)重。設(shè)W為所需的權(quán)值矩陣，它是對稱和正定的，則加權(quán)的目標(biāo)函數(shù)為(3.3.6)

按上述求極小值的方法，可得加權(quán)的最小二乘估計量：(3.3.7)顯然，當(dāng)W選為單位矩陣時，WLS

=LS。

例3-7

考慮最簡單的一維線性模型（靜態(tài)的），即只有一個控制變量u的情形，這時模型的形式是2/1/2023第三章參數(shù)估計理論與應(yīng)用求未知參數(shù)θ0

和θ1的LS估計量。

解實際過程輸出是模型的輸出加上一隨機誤差項，即觀測數(shù)據(jù)對[ui,yi]的結(jié)構(gòu)應(yīng)為式中，ei

稱為模型的殘差或觀測噪聲，一般認(rèn)為是零均值、相互獨立的隨機序列，并具有相同的方差σ2。將上式寫成矩陣形式：

根據(jù)式(3.3.5)，可得LS估計量：2/1/2023第三章參數(shù)估計理論與應(yīng)用 如果進(jìn)一步假定ei

的分布是正態(tài)的，則容易驗證，方差σ2

的ML估計量是 作為練習(xí)，請讀者在Matlab平臺上輸入以下數(shù)據(jù)和函數(shù)：x=[12345];y=[1.31.82.22.93.5];

[p,s]=polyfit(x,y,1) %生成擬合一次多項式運行結(jié)果是：p=[0.550.69]，s=0.1643。即y=0.55x+0.69標(biāo)準(zhǔn)差為0.1643。2/1/2023第三章參數(shù)估計理論與應(yīng)用

例5-8

（可線性化的非線性靜態(tài)模型——曲線回歸）假設(shè)有一個非線性模型的輸出為其中，x1，x2

為確定性輸入變量，a，b和c為待估計參數(shù)。

解

上式兩邊經(jīng)簡單的代數(shù)運算，再同時取自然對數(shù)，可轉(zhuǎn)化為一個線性模型：這說明變換后的輸出ln(y-1-1)

可以顯式地表達(dá)為以

lnx1和

x2為輸入、以lna，b和

c

為參數(shù)的線性模型。因此，就可以按變換后的線性模型用最小二乘法來估計變換后的未知參數(shù)，然后，再根據(jù)變換后的估計參數(shù)計算出原參數(shù)。2/1/2023第三章參數(shù)估計理論與應(yīng)用 判定輸入x-輸出y之間的關(guān)系能否用一個線性模型來描述的標(biāo)準(zhǔn)，通常用互相關(guān)系數(shù)的大小來衡量：(3.3.8)ρxy

的絕對值越大，表示變量之間的線性關(guān)系越密切，因而線性回歸的效果就越好。

例3-9

設(shè)某一結(jié)構(gòu)參數(shù)n，m和d已知的離散線性系統(tǒng)，其差分方程的形式為：

(3.3.9)2/1/2023第三章參數(shù)估計理論與應(yīng)用其中，e(k)為噪聲，φ(k)為輸入-輸出觀測向量，θ為未知參數(shù)向量，且要求根據(jù)N次數(shù)據(jù)對{[y(i),u(i)],i=1,2,…,N；N≥n+m+1}來估計對未知參數(shù)θ。

解將式(3.3.9)改寫成矩陣形式，得到將數(shù)據(jù)寫成下標(biāo)形式，就有這樣，未知參數(shù)向量θ可按式(3.3.5)進(jìn)行估計。 2/1/2023第三章參數(shù)估計理論與應(yīng)用

考慮如下單輸入-單輸出系統(tǒng)：用Matlab中rarx函數(shù)進(jìn)行系統(tǒng)辨識，程序如下：

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); %實際系統(tǒng)的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);%

根據(jù)[z]進(jìn)行ARX模型參數(shù)辨識

plot(y,'-');grid %

作圖，實際系統(tǒng)的輸出曲線

holdon

plot(yhat,':') %

作圖，辨識系統(tǒng)的輸出曲線

參數(shù)辨識結(jié)果thm：

a?1=-1.3798，a?2=0.7039，b?1=0.3007，

b?2=0.1170，b?3=0.4243。2/1/2023第三章參數(shù)估計理論與應(yīng)用 應(yīng)當(dāng)指出，要求觀測數(shù)據(jù)容量N≥n+m+1是為了保證ΦTΦ非奇異，降低過程噪聲序列{e(k)}的影響，從而提高參數(shù)估計的精度。不論{e(k)}是何種形式的噪聲序列，式(3.3.5)總是成立的。換言之，噪聲性質(zhì)僅影響LS估計的統(tǒng)計特性。 下面介紹LS估計的統(tǒng)計特性。如果觀測噪聲或建模誤差序列{e(k)}具有零均值和相同的方差，即則LS估計量LS是無偏、有效和相容的，并具估計誤差的協(xié)方差為σ

2(ΦTΦ)-1。 對于動態(tài)控制系統(tǒng)的辨識，輸入信號u(t)必須滿足持續(xù)激勵條件，也即輸入信號u(t)的頻譜必須包含足夠豐富（Sufficientrich）的頻率成分，以保證充分激勵受控對象的2/1/2023第三章參數(shù)估計理論與應(yīng)用所有振型，從而使觀測數(shù)據(jù)載有動態(tài)系統(tǒng)的主要信息。LS估計在滿足持續(xù)激勵條件時，是漸進(jìn)無偏的，也稱為估計的一致性。 式(3.3.9)，也稱為CAR模型（即受控的AR模型），可以寫成更簡潔的形式

(3.3.10)式中q

表示時間算子，d為整數(shù)，表示系統(tǒng)的滯后量；A(·),B(·)分別為q－1

的降次冪多項式。

CAR模型滿足一致估計（或相容估計）的條件為：

（1）{e(k)}是白噪聲序列； （2）{u(k)}的均值和協(xié)方差有界；且滿足（m+1）階持續(xù)激勵條件（或正定條件）：2/1/2023第三章參數(shù)估計理論與應(yīng)用

（3）u(k)

和e(k)

相互獨立。通常，u(k)都采用偽隨機二元序列。 只要選擇恰當(dāng)?shù)哪Ｐ碗A次或最小二乘多項式階次(參見taylor.m,Matlab)，最小二乘法總是可以很好地擬合數(shù)據(jù)，但是，如果觀測數(shù)據(jù)波動較大，將嚴(yán)重影響參數(shù)估計的準(zhǔn)確性。對此，可采用數(shù)據(jù)預(yù)處理和數(shù)字濾波的方法加以解決。 檢驗?zāi)Ｐ蜏?zhǔn)確性的最簡單方法是準(zhǔn)備另外一組輸入-輸出數(shù)據(jù)對，稱為檢驗數(shù)據(jù)集，在參數(shù)估計時不用，待模型建立后，用這組數(shù)據(jù)對來驗證所得模型的普適性或泛化能力。2/1/2023第三章參數(shù)估計理論與應(yīng)用

上述討論，均假設(shè)數(shù)學(xué)模型的階次是已知的。實際上，對于動態(tài)系統(tǒng)，模型的階次很少是預(yù)知的。檢驗?zāi)Ｐ碗A次是否合適的一種簡單而有效的方法是：評估不同階次的理論模型對觀測數(shù)據(jù)的擬合度，用擬合誤差函數(shù)

來描述。通常，當(dāng)n

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

或J(n,m)就會減小；而當(dāng)n(或n,m)大于模型的真實階次

n0(或n0,m0)時，J的減小就不顯著了。由此，可以很方便地用多次實驗的方法來確定模型的恰當(dāng)階次。 注意，對于J(n,m)形式的擬合誤差函數(shù)，一般應(yīng)按正交實驗法來確定模型的恰當(dāng)階次，以減少實驗的次數(shù)。2/1/2023第三章參數(shù)估計理論與應(yīng)用3.3.2遞推最小二乘估計 在測控系統(tǒng)中，被測對象通常可以不斷提供新的輸入-輸出數(shù)據(jù)。如果希望利用新的信息來改善估計精度，那么就應(yīng)當(dāng)采用遞推估計算法，這不僅可避免觀測數(shù)據(jù)矩陣Φ的行數(shù)的不斷“膨脹”，而且可減少參數(shù)估計的計算量。 在推導(dǎo)最小二乘遞推算法前，先引入一個與式(3.2.10)類似的矩陣求逆定理。設(shè)A和I+CA-1B

均是非奇異方陣，則

(3.3.11)

下面介紹最小二乘遞推算法。為了簡化符號，以下推導(dǎo)均用代替LS。 設(shè)ΦN是時刻N為止的觀測數(shù)據(jù)，N+1時刻θ的LS估計為2/1/2023第三章參數(shù)估計理論與應(yīng)用式中[參見式(3.3.2)和(3.3.9)]于是有(3.3.12)

令PN=[ΦN

TΦN]-1

，由求逆公式(3.3.11)知

(3.3.13)

定義增益向量為2/1/2023第三章參數(shù)估計理論與應(yīng)用將式(3.3.13)和(3.3.14)代入(3.3.12)，得到(3.3.15)

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

等于前一時刻的估計量N

與修正項KN+1(yN+1-φTN+1N)

之和，這是一切遞推公式的共同特征。如果令代表基于前一時刻的估計量N

對N+1時刻的預(yù)測。那么，遞推估計提供的新息2/1/2023第三章參數(shù)估計理論與應(yīng)用就是預(yù)測誤差或擬合誤差。因此，修正量的大小與新息成正比，而各校正分量的權(quán)，由增益向量決定。 在啟動上述遞推算法時，必須知道初值0

和P0

，通常令其中，σ2?1。然后從得到的一組數(shù)據(jù)，按式(3.3.16)開始遞推運算。 從物理上看，這種初值選取雖然初始誤差較大，但校正的作用也大，因此這種遞推算法是有效的。此外，還可以先取得N>m+n+1組數(shù)據(jù)，算出2/1/2023第三章參數(shù)估計理論與應(yīng)用作為初值，然后，按式(3.3.16)進(jìn)行遞推運算。 增益向量KN+1在遞推運算過程是怎樣變化的？先考察PN=[ΦN

TΦN]-1如果N大于估計參數(shù)的數(shù)目，而且，輸入－輸出數(shù)據(jù)對含有足夠的“信息”（滿足充分激勵條件），則

ΦNTΦN

通常是正定的。顯然，當(dāng)N

趨于無窮大，ΦNTΦN

/N

接近于非奇異的常數(shù)陣。于是，有 可見，式(3.3.16)中自適應(yīng)增益向量KN+1

隨著每次迭代而遞減，這意味著遞推運算過程將逐漸收斂于參數(shù)空間的最優(yōu)點。事實上，在白噪聲或低噪聲條件下，遞推最小二乘2/1/2023第三章參數(shù)估計理論與應(yīng)用估計是一種簡便而又有效的算法。這種遞推算法在遞推過程中雖然沒有保存全部先前的數(shù)據(jù)，但所有先前數(shù)據(jù)的影響卻一直在起作用，故稱為無限增長記憶的遞推最小二乘估計。3.3.3卡爾曼濾波器的遞推算法（狀態(tài)估計） 與參數(shù)估計器不同，卡爾曼濾波器主要是解決如何從被噪聲污染的觀測數(shù)據(jù)中估計出已知動態(tài)系統(tǒng)模型的狀態(tài)，而不是動態(tài)系統(tǒng)模型的未知參數(shù)。然而，僅從算法上看，這二者是非常相似的。為了便于比較二者的異同，我們在此不加證明地列出卡爾曼濾波器算法。 設(shè)離散定常系統(tǒng)的狀態(tài)方程和觀測方程分別為2/1/2023第三章參數(shù)估計理論與應(yīng)用式中，X(k)是n維狀態(tài)向量，

u(k)是m維控制向量，w(k)是

p維過程噪聲，y(k)是r維觀測向量，v(k)是r維觀測噪聲；

A，B，C，Γ分別是相應(yīng)維數(shù)的系統(tǒng)矩陣、控制矩陣、觀測矩陣和過程噪聲權(quán)矩陣。噪聲的統(tǒng)計特性滿足

假定到k時刻為止，觀測數(shù)據(jù)為{y(1),…,y(k)}，要估計l時刻的狀態(tài)，就有三種情況：l>k，稱為預(yù)測問題；l=k，稱為濾波問題；l<k，稱為平滑問題。 下面主要介紹濾波問題。2/1/2023第三章參數(shù)估計理論與應(yīng)用濾波算法：預(yù)測算法：濾波增益：濾波誤差的協(xié)方差：預(yù)測誤差的協(xié)方差：初始條件：2/1/2023第三章參數(shù)估計理論與應(yīng)用

例3-10

考慮某一動態(tài)系統(tǒng)表示的空間導(dǎo)航問題，其中加速度為白噪聲。被白噪聲污染之后觀測其位置。因此，過程的狀態(tài)方程為（1）式中，X(t)=[x1(t),x2(t)]T,表示過程的位置和速度；而w(t)是一均值為0、方差為1的高斯白噪聲。觀測方程為（2）其中v(t)是一均值為0、方差為10的高斯白噪聲。

解

首先把連續(xù)的狀態(tài)方程離散化（參見Matlab中c2d函數(shù)）2/1/2023第三章參數(shù)估計理論與應(yīng)用不妨設(shè)采樣間隔T=1，則可寫出狀態(tài)方程（1）和觀測方程（2）的離散化表達(dá)式進(jìn)一步假設(shè)則可按前面介紹的濾波算法進(jìn)行遞推計算：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函數(shù)仿真卡爾曼濾波器的設(shè)計。請讀者在Matlab平臺上完成例3-10的仿真計算。2/1/2023第三章參數(shù)估計理論與應(yīng)用3.3.4限定記憶的遞推最小二乘估計

以上介紹的遞推最小二乘法適用于定常系統(tǒng)的參數(shù)估計。對于時變系統(tǒng)，由于參數(shù)時變的信息顯然更多地蘊藏在當(dāng)前的觀測數(shù)據(jù)中，而與先前觀測數(shù)據(jù)的關(guān)系將逐漸減弱。因此，利用一切觀測數(shù)據(jù)對來決定的自適應(yīng)增益向量KN+1，顯然會削弱遞推過程跟蹤變化參數(shù)的能力。解決這一問題的一種簡單方法是：當(dāng)我們懷疑觀測數(shù)據(jù)發(fā)生顯著的變化時，就將當(dāng)前的PN

設(shè)置為P0，重新進(jìn)行參數(shù)估計，這是因為LS估計能快速地收斂到當(dāng)前的最優(yōu)參數(shù)。處理這一問題的另一種方法是對過去的數(shù)據(jù)引入帶遺忘因子λ，逐漸削弱它們在參數(shù)估計中的作用。為此，采用加權(quán)的目標(biāo)函數(shù)：2/1/2023第三章參數(shù)估計理論與應(yīng)用其中

(3.3.17)

而且0<λ≤1，當(dāng)λ=1，就退化為基本的遞推最小二乘法。由式(3.3.7)知，在時刻N，使上述加權(quán)的目標(biāo)函數(shù)最小化的估計為

(3.3.18)

每當(dāng)取得一個新的數(shù)據(jù)后，就對溝渠以前的加權(quán)矩陣乘以λ，于是N+1時刻的估計量為2/1/2023第三章參數(shù)估計理論與應(yīng)用與前面略有不同，在此令PN=[ΦN

TWNΦN]-1

，則由求逆公式(3.3.11)，可知于是，按前面推導(dǎo)出式(3.3.15)思路，可得式中

遺忘因子λ將老的數(shù)據(jù)逐漸從“記憶”中去掉，因此這種使用數(shù)據(jù)信息的方式也叫做“漸消記憶”法，相應(yīng)的算法稱為帶遺忘因子的遞推最小二乘法，即2/1/2023第三章參數(shù)估計理論與應(yīng)用關(guān)于λ的選取通常由經(jīng)驗或?qū)嶒灤_定，一般范圍為0.95≤λ≤0.99 λ

取得愈小，最新數(shù)據(jù)的權(quán)重就愈大，也就更適合于跟蹤大的時變參數(shù)，但與此同時，估計器也可能會發(fā)生較大的波動 從而加大估計的誤差。 例5-10

考慮模型采樣300次后變?yōu)?/1/2023第三章參數(shù)估計理論與應(yīng)用

試用兩組模擬數(shù)據(jù)，一組不考慮噪聲，一組是帶觀測噪聲的數(shù)據(jù)，分別用不同的遺忘因子，對時變模型進(jìn)行參數(shù)估計，并討論估計結(jié)果。

解用MATLABrarx函數(shù)進(jìn)行帶遺忘因子λ的系統(tǒng)辨識算法。程序如下：

e=randn(300,1);u=idinput(300,‘prbs’); %產(chǎn)生高斯噪聲和偽隨機信號

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第三章參數(shù)估計理論與應(yīng)用

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

%產(chǎn)生輸入-輸出信號[z]和ARX模型的階次

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

nbnk],'ff',0.97);

%帶遺忘因子0.97的ARX參數(shù)辨識

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

%作圖，時變系統(tǒng)的a(k)輸出曲線

holdon

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

%作圖，時變系統(tǒng)的b(k)輸出曲線圖3-5帶遺忘因子λ的時變系統(tǒng)參數(shù)辨識結(jié)果2/1/2023第三章參數(shù)估計理論與應(yīng)用 從仿真結(jié)果（圖3-5）可以看出，采用帶遺忘因子λ的系統(tǒng)辨識算法能較好地跟蹤系統(tǒng)時變參數(shù)的變化。3.3.5廣義最小二乘估計 上述推導(dǎo)中，假定模型(3.3.10)中的噪聲項e(t)為白噪聲。如果e(t)為有色噪聲，則被測對象就應(yīng)當(dāng)?shù)挠肅ARMA模型（受控的ARMA）來表示，即與CAR模型（基本最小二乘估計的數(shù)學(xué)模型）比較，此處C(q-1)≠1，而有色噪聲e(t)

可以認(rèn)為是由白噪聲ξ(t)通過成型濾波器而產(chǎn)生的輸出（以下寫成離散數(shù)據(jù)的形式）：2/1/2023第三章參數(shù)估計理論與應(yīng)用 下面介紹廣義最小二乘法，用以改善CARMA模型參數(shù)估計的統(tǒng)計特性。 將CARMA模型寫成常規(guī)的最小二乘結(jié)構(gòu)，即(3.3.20)

因為ξ(k)是不可測量的，所以在初始估計時只能用ξ(k)估計值代替φ(k)中有關(guān)ξ(k)的分量。按基本最小二乘法估計出?，然后，計算ξ(k)的估計量(3.3.21)再構(gòu)造2/1/2023第三章參數(shù)估計理論與應(yīng)用這樣，就可以按前面介紹的遞推公式(3.3.22)

對CARMA模型進(jìn)行參數(shù)估計。估計的初始值可按下式選取其中，σ2?1。在實際應(yīng)用中，ξ(k)的估計有兩種形式，一種是預(yù)測的形式：2/1/2023第三章參數(shù)估計理論與應(yīng)用為了避免這種問題的產(chǎn)生，最常采用的方法是對PN

采 用U-D分解法,下面就介紹這種方法。3.3.6改進(jìn)數(shù)值穩(wěn)定性的U-D分解法 為了保證在遞推過程中始終保持PN

的非負(fù)定性，可把矩陣PN

進(jìn)行U-D分解，即【在此，將PN

寫成P(N)】其中，U是對角元全為1的上三角矩陣，D為對角矩陣。這種分解只需要實時修正U，而無需修正P，因而能保證P的非負(fù)定性。下面，我們不加證明地將U-D分解基本算法的計算步驟總結(jié)如下：

2/1/2023第三章參數(shù)估計理論與應(yīng)用

(1)

設(shè)置初值θ?(0)，U(0)和D(0)，輸入初始數(shù)據(jù)；

(2)

采樣當(dāng)前的輸入和輸出（假設(shè)數(shù)據(jù)的行數(shù)為p）；

(3)

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

(4)

定義由此可得2/1/2023第三章參數(shù)估計理論與應(yīng)用按上式計算βi(N)； （5）計算

（6）計算參數(shù)估計 （7）按下式計算di(N)式中，di是D的對角元素； （8）引入上三角陣