1、第三篇 最優估計理論概 述在科學和技術領域中，經常遇到“估計”問題。所謂“估計”，就是對受到隨機干擾和隨機測量誤差作用的物理系統，按照某種性能指標為最優的原則，從具有隨機誤差的測量數據中提取，信息估計出系統的某些參數狀態變量。這就提出了參數和狀態估計問題。這些被估參數或被估狀態可統稱為被估量。一般，估計問題分兩大類，即參數估計和狀態估計。一、參數估計參數估計屬于曲線擬合問題。例如做完某項試驗之后，得到若干個觀測值 與相應時間 的關系 。我們希望以一條曲線來表示 和 的關系，設izit,1,2,iiz timzt 1 122nnz tx h tx htx ht式中 為已知的時間函數，一般是 的冪

2、函數、指數函數或正余弦函數等等。 為 個未知參數，它們不隨時間而變。 12nh ththt、 、t12nxxx、 、 、n根據 對觀測值 來估計未知參數。按照什么準則來估計這些參數呢？這將是第十章討論的主要問題。m,1,2,iiz timmn；12nxxx、 、 、二、狀態估計設系統的狀態方程和觀測方程分別為 ttttttttH tttxAxBuFwzxv 式中， 為狀態變量，它是隨時間而變的隨機過程， 為控制變量， 為系統噪聲， 為測量噪聲， 為觀測值。現要根據觀測值來估計狀態變量 ，這就是狀態估計問題。卡爾曼濾波是一種最有效的狀態估計方法，將在第十一章討論這個問題。 tx tu tw tv

3、 tz tx 人們希望估計出來的參數或狀態愈接近真值愈好，因此提出了最優估計問題。所謂最優估計，是指在某一確定的準則條件下，從某種統計意義上來說，估計達到最優，顯然，最優估計不是唯一的，它隨著準則不同而不同，因此在估計時，要恰當選擇估計準則。 在自動控制中，為了實現最優控制和自適應控制，遇到許多參數估計或狀態估計問題，促進了估計理論和估計方法的發展。另外，由于電子計算機的迅猛發展和廣泛使用，使得許多復雜的估計問題的解決成為可能，這也促進了估計理論的發展。所以近二十多年來最優估計理論及其應用得到迅速的發展。第十一章 參數估計方法 本章討論參數估計準則和估計方法，根據對被估值統計特性的掌握程度不同

4、，可提出不同的估計準則。依據不同的準則，就有相應的估計方法，即最小方差估計、線性最小方差估計、極大似然估計、極大驗后估計、最小二乘估計等，本章將對這些估計方法進步不同程度的討論。第一節 最小方差估計與線性最小方差估計 一、最小方差估計一、最小方差估計 最小方差準則，要求誤差的方差為最小，它是一種最古典的估計方法，這呼估計方法需要知道被估隨機變量 的概率分布密度 和數學期望 。這種苛刻的先驗條件，使此方法在工程上的應用受到很大限制。這里只以一維隨機變量的估計為例，介紹最小方差估計方法。x p x E x將上式展開，得 2222JExxE xxE xx 設有一維隨機變量 ，它的概率密度 和常數期望

5、 ，都是已知的，求 的估值 。評價估計優劣的準則是 與 的誤差的方差為最小，即x p x xE xmx x xx 22minJExxxxp x dx(11-1)求上式對 的偏導數，令偏導數等于零，得 x 22JxE xx因此 的最小方差估值為 ，估計誤差為xxm 0 xxxxxxxxmE xE xR xE xmmm即 E xE x則 的最優估值為x xxE xxp x dxm(11-2) 如果估值 的數學期望等于 的數學期望，或者估計誤差 的數學期望為零，則最小方差估計是無偏的。因此 的估計是無偏估計。 xx xx所以數學期望 是 的最小方差估計。xmx這種方法可以推廣到多維隨機變量的估值，這

6、里不再敘述。估計誤差 的方差為 x 222xxxExmxmp x dx(11-3)二、線性最小方差估計二、線性最小方差估計 線性最小方差估計就是估計值為觀測值的線性函數，估計誤差的方差為最小。在使用這種方法時，需要知道觀測值和被估值的一、二階矩，即數學期望 和 、方差Varz和Varx及協方差 和 。 E z E xCov xz，Cov zx，的條件來確定系數 和 。ab式中 和 為兩個待定常數。根據估計誤差的方差ab22JExxExazbmin(11-5) 先討論被估值 和觀測值 都是一維隨機變量的情況。線性最小方差估計是把 的估值 表示成 的線性函數，即xzxz x xazb(11-4)

7、求式(11-5)對 和 的偏導數，令偏導數等于零，可求得 和 兩個系數。abba202JExazbzaJExazba (11-7)(11-6)從式(11-7)可得0 xzmamb式中 和 為 和 的數學期望，從此式可得xmzmxzxzbmam(11-8)將式(11-8)代入式(11-6)得0 xzExazmamz把上式改寫成 0 xzzzExma zmzmm展開上式得02xzzxzzzExmzmm E xmaEzmam E zm 式中， 分別為隨機變量 和 的均方根差， 為 與 的相關系數 。于是的估值為xz、xzxzxzCov,/xzxzax z 2Cov,xzxx zxazbmzm(11-

8、10)求上式的數學期望值，可得222220Cov,xzxzxzzxxxCov xzax za ，(11-9)估計誤差為xxx 0 xzzzzxxzzxzzzxE xE xE mE zmmmmm 因此 。所以估計是無偏的。 E xE x 下面討論 和 都是多維隨機變量的估計問題。設 為 維， 為維，已知 和 的一、二階矩，即 xzxxzznq VarVarzCov,Cov,E xE zx zz x、和假定 的估值 是 的線性函數xz x z xbAz式中， 是 維非隨機常數向量， 是 非隨機常數矩陣。nbAnq(11-11)估計誤差方差陣為( )( ) TTJTraceExx zxx zExbA

9、zxbAz(11-12) 估計準則是方差陣J為最小，也可等價為方差陣J的跡 為最小，即 的各分量的方差之和為最小tJx minTtTJTraceEExbAzxbAzxbAzxbAz(11-13) 用函數對矩陣的微分法則（附錄一），求 和 ，令 和 ，聯立求解可得 。tJbtJA0tJb0tJAbA和220tzzJE xbAzbAmmb zzEEbmAmxAx 22220TtTTTTTTJEEEEEEE xbAzxbAzAAxbAz zxbxAzzAzzbzxx(11-15)(11-14)將式(11-14)代入式(11-15)得 000( , )0TTTTTTTTTTEEEEEEzEE z EE

10、EEEEEEEE zAVarzCov x zAzzxzAzzxAzzzxzxzAzzzzxxz因此 1CovVarAxzz、(11-16)將式(11-16)代入式(11-14)，可得 1CovVar EEbxxzzz、根據式(11-16)和式(11-17)求得 代入式(11-11)，得Ab和 11Cov,VarCov,VarzEExxxx zzzzmx zzzm(11-18) 1Cov,VarzEEEExxxmx zzzmmx由式(11-18)可得所以估計是無偏的。估計誤差的方差陣為1Var -CovVarCovJxxzzzx、(11-19)第二節 極大似然法估計與極大驗后法估計 一 、極大似

11、然法估計 極大似然法估計是以觀測值出現的概率為最大作為估計準則的，它是一種常用的參數估計方法。 設 是連續隨機變量，其分布密度為 ，含有 個未知參數 。把 個獨立觀測值 分別代入 中的 ，則得z12,np z n12,n k12,kz zz12,np z z12,1,2,inp zik ，稱函數 為似然函數。當 固定時， 是 的函數。極大似然法的實質就是求出使 達到極大時的 的估值。從式(11-20)可看到 是觀測值 的函數。L12,kz zzL12,n L12,n 12,n 12,n 12,kz zz將所得的 個函數相乘，得k1212121,kkniniL z zzp z ， ；(11-20

12、) 為了便于求出使 達到極大的 ，對式(11-20)取對數，則L12,n 1121lnln,kniLp z (11-21) 由于對數函數是單調增加函數，因此當 取極大值時， 也同時取極大值，將上式分別對 求偏導數，令偏導數等于零，可得下列方程組：LlnL12,n 1nln0ln0LL(11-22) 解上述方程組，可得使 達到極大值的 。按極大似然法確定的 ，使 最有可能出現，并不需要的驗前知識，即不需要知道的概率分布密度和一、二階矩。L12,n 12,n 12,kz zz12,n 例11-1 設有正態分布隨機變量 ，給出 個觀測值 。觀測值相互獨立，試根據這 個觀測值，確定分布密度中的各參數。

13、zk12,kz zzk解: 的分布密度可用下式表示：z221,exp22zmp z m式中的 和 為未知參數。 m現有極大似然法來確定參數 和 。 m作似然函數： 212211,exp22kikizmL z zzm；對上式取對數，可得212212211ln,lnexp221lnln22kikikiizmL z zzmzmkk ；將上式分別對 和 求偏導數，令偏導數等于零，可得m21231lnL0lnlnL10kiikiizmkzm聯立求解可得2211kkiiiizzmmkk，上面介紹了極大似然法的基本概念。現在來討論極大似然法估計參數的問題。 設 為 維隨機變量， 為 維未知參數，假定已知 的

14、條件概率密度 。現在得到 組 的觀測值。觀測值相互獨立。當參數 是何值時， 出現的可能性最大？為此，確定似然函數：/p z xzxmnzzkx12,kz zz 12,/kL z xp zx p zxp zxp z x或 ln,ln/L z xp z x(11-23)(11-24)取極大值的充分條件是L2222ln00LLxx或 因此，用極大似然法時，應先求似然函數 ，然后用微分法求出使似然函數 為極大的的 估值 。LLx x解之，可得 的估值 。x x求出使 為極大的 值，令 Lxln00LLxx或(11-25)下面求似然函數 1,/p xvL z xp z xpx，設有一線性觀測系統zb z

15、,vHzV 式中， 是 維觀測值， 是 維未知參數， 是 維測量誤差。設與 獨立。給出 的統計特性，求 的極大似然估計。zmxnvmvxvx(11-26) 根據不同隨機變量的概率密度變換公式，并考慮到 與 獨立，可得xv 1212221,p x zp x Hzvp x vpx pvpx pvL z xpvpzHxpx令 2,0LpHz xzxxx得上式，可得的 估值 。x x假定噪聲 是正態分布，其均值為零，方差陣為 ，則vTE vvR 121/211,exp -22TmL z xpvvvRR把 代入上式，得vzHx121,exp -2TL z xpczHxzHxRzHx式中 1/212mcR

16、求出 ，使 為最大，也就是使x2,L z xpzHx11min2TJzHxRzHx(11-27)求 對 的偏導數，令偏導數等于零，可得 的估值Jxx x110TT JH R zH R Hxx11TTH R HxH R z111TTxH R HH R z(11-28)二、極大驗后估計二、極大驗后估計 如果給出 維隨機變量 的條件概率分布密度 也稱驗后概率密度，怎樣求 的最優估值 呢？極大驗后估計準則：使 的驗后概率密度 達到最大那個 值為極大驗后估值 。可見，極大驗后估計是已知 求 的最優估值 的一種有效方法。nx/p x zx xx/p x zx x/p x z x 式中 是 的驗前概率密度，

17、 是觀測值 的概率密度，可用計算方法或實驗方法求得。為了計算 需要知道 。在 沒有驗前知識可供利用時，可假定 在很大范圍內變化。 p xx p zz/p z x/p x z p xxx 極大驗后估計是以已知 為前提的。如果只知道 ，可按下式計算 。/p x z/p z x/p x z /p z x p xp x zp z(11-29) 在這種情況下，可把 的驗前概率密度 近似地看作方差陣趨于無限大的正態分布密度x p x 11/211exp -22TxxnpxxmPxmP式中 為 的方差陣， 為 單位陣， ，于是xP PII，nn1P0 /21/2111lnln 2-2lnnTxxxpp xP

18、xmPxmxPxmx(11-30)當 1P0 ln0pxx(11-31) 當缺乏 的驗前概率分布密度時，極大驗后估計與極大似然估計是等同的，現證明如下：x對于極大似然估計，為了求得 的最優估值 ，應令x xln/0pz xx(11-32)對于極大驗后估計，為了求得 的最優估值 ，應令x xln/0px zx(11-33)根據式(11-29)得 ln/ln/ln-lnln/ln/lnln0ppppppppx zz xxzx zz xxzxxxx考慮到 不是 的函數，同時考慮到式(11-31)，可得 p zxln/ln/ppx zz xxx 一般說來，極大似然估計比極大驗后估計應用普遍，這是由于計

19、算似然函數比計算驗后概率密度較為簡單。(11-34) 第三節 最小二乘法估計與加權最小二乘法估計 上面討論的幾種估計方法，分別對被估隨機變量 的概率分布密度 、條件概率密度 、 以及一、二階矩等條件有著不同的要求。假定我們并不掌握上述任何條件，仍要估計隨機變量 的最優估值 ，只有用高斯提出的最小二乘法。x p x/p z x/p x zx x一、最小二乘法估計一、最小二乘法估計 設 次獨立試驗，得到 對觀測值： 。這里 表示時間或其他物理量。現在的任務是：根據這些觀測值，用最優的形式來表示 與 之間的函數關系。m 1122,mmz tz tzt， ，itztm 通常， 的未知函數可用 表示，

20、的類型應根據這 對數據（ 個點）的分布情況或所研究問題的物理性質來確定。 f tz f tmm為了便計算，可采用多項式 21123nnf txx tx tx t(11-35)來表示，也可以用更一般的形式表示： 1 12212nnnf tx h tx htx htf txxx， ， ， ，在式(11-36)中， 為已知的確定性函數，如 的冪函數、正余弦函數和指數函數等。 為 個待定的未知數。 12nh ththt， ，t12nxxx， ， ，n(11-36) 把 對觀測值代入式(11-35)或(11-36)，可得 個方程式。如果 ，即方程數 少于未知參數的數目。則方程的解不確定，不是唯一地確定解

21、出 。當 時，方程數正好與未知參數的數目相等，能唯一地解出 。在這種情況下， 曲線一定通過每一個觀測點 。因為在觀測結果中，不可避免地含有隨機測量誤差，如果曲線通過每一個觀測點，則曲線將包含這些測量誤差，反而不能真實地表達出 與 之間的正確函數關系。所以不應要求 曲線一定通過每一個觀測點。mmmnm12nxxx， ， ，mn12nxxx， ， ， f t( , )(1,2,)iiz timzt f t 一般，試驗次數 ，而且希望 比 大得多，即方程數大于未知參數數目這種情況只能采用數理統計求估值的方法。下面討論這一問題。mnmn 確定了函數 的類型之后，問題就歸結為如何合理地選擇 中的參數 ，

22、使得這一函數在一定意義下比較準確地反映出 與 的函數關系。通常用最小二乘法來選擇這些參數。所謂電波二乘法，就是要求所選擇的 的參數，使得觀測 與對應的函數 的偏差的平方和為最小。設 為觀測值 與對應函數 的偏差的平方和，即 if t f t f t12nxxx， ， ，zt f tizJiz if t 221211,mmiiiniiiJzf tzf x xx t(11-37) 按照 為最小的條件來確定 中的參數 ，將上式分別對 求偏導數，并令它們等于零，可得下列方程組：J f t12nxxx， ， ，12nxxx， ， ，12121211212112121,0,0,0nminixniimini

23、xniiminixniizf x xx tfx xx tzf x xx tfx xx tzf x xx tfx xx t 上述方程組有 個方程， 個未知數，解之可得 的最優估值 。nn12nxxx， ， ，12nxxx， ， ，12nxxx， ， ，(11-38)例11-2 觀測值 和觀測時刻 如表11-1所示。izit245892.012.983.505.025.47itiz表11-1 設 。用最小二乘法確定 和 。 12f txx t1x2x解: 1251215121100 xxiiiiiiifftzxx tzxx t t，把 值分別代入上面兩個方程，經過整理后可得iizt，12125.6

24、3.7966.78574.3846xxxx解上述方程組，可得 的估值：12xx和121.06160.496xx，所以 1.0160.496z tf tt實際上 與 不可能完全一致，這是由于以下幾個原因引起：iz if t 選得不夠準確；ix 存在觀測誤差； 的模型方程 選得不夠確切。z f t當選定 之后，可得觀測值 與 之差 xiz12 ,nif xxx t， ， ，12 ,1,2,iiniezf xxx tim， ， ，(11-39)式(11-39)可寫成12 ,inizf xxx t， ， ，考慮到式(11-36)，上式可寫成 1 1221,2,iiinniizx h tx htx ht

25、eim(11-40)或寫成 11 112211121 12222221 122nnnnmmmnnmmzx h tx htx htezx h tx htx htezx h tx htx hte(11-41)如果設111222mnmzxezxezxezxe 112111222212nnmmnmh ththth ththth ththt12mhhHh(11-43)(11-42) 式中， 。則式(11-41)可寫成下列矩陣形式 121,2,iiinih ththtimh， ，zHxe(11-44)求 對 的偏導數，令偏導數等于零，可得J x20TTTTJH zH HxxH HxH z 下面用矩陣形式不

26、表示最小二乘的公式。殘差的平方可用下式表示：21mTiiJee e(11-45)從式(11-44)得 TJezHxzHxzHx(11-46)(11-47)因而 的估值為x1TTxH HH z為使 能求得估值 ，逆陣 必須存在。J x1TH H(11-48) 為極小的充分條件是J2220TJH Hx即 為正定矩陣。TH H(11-49)當 時， 為 階方陣，且 存在時，則mnHn1H1xH z在一般情況下， ，因此 要用式(11-48)來求。mn x(11-50)例11-3 觀測值 和觀測時刻如例11-2，試用式(11-48)求 的系數 和 。iz z tf t1x2x解: 121 122121

27、f txx tx h tx hth thtt，則 12121122121222124414551588189919hhhhxhhxhhhh xH，123452.012.98111112.50245895.025.47TzzzzzzH，15281.144570.16867281900.168670.0301TTH HH H122.012.981.144570.16867111111.0162.500.168670.0301245890.4965.025.47xx 因此 = =1.016+0.496 z tf tt二、加權最小二乘法二、加權最小二乘法 當 的測量精度高時， 大；反之， 小。這樣可使

28、擬合曲線接近于測量精度高的點，從而保證擬合曲線有較高的準確度。iziwiw 在式(11-45)中，每個誤差值 的系數都為1，即在 中每個誤差值都是“等權”的。事實上，在 值的不同測量范圍內，測量精度往往是不同的，因而測量誤差也不相同。合理的辦法是對不同的誤差項 加不同的權，即把 寫成2ieJz2ieJ2211mmiiiiiiiJwew zh x(11-51) 式中 為 對稱矩陣，稱為加權矩陣。求 對 的偏導數，并令其等于零，可得WmmJ x220TTTTJH H WzH W xxH WH xH Wz則可得 1TTxH WHH Wz(11-53)把式(11-51)寫成TJ zHxW zHx(11

29、-52)為極小的充分條件為J2220TJH WHx即 為正定矩陣。TH WH(11-54)由于測量 時，存在測量誤差 ，故觀測方程為zvzHxv 和 都為 維向量。假定 ， 不一定是正態分布。xvm 0E vTEvvRv，(11-55)由式(11-53)得加權最小二乘法的估計誤差11TTTTxxxxH WHH WzH WHH W Hxz考慮到式(11-55)，得1TT xH WHH Wv(11-56)由上式可得 10TTE xH WHH W v所以，在上述條件下的加權最小二乘法估計為無偏估計。(11-57)可以證明， 使估計誤差的方差為最小，因而 為最優的加權矩陣。有時，人們把式(11-59)

30、稱為馬爾柯夫估計。1WR1WR考慮式(11-56)，估計誤差的方差程為11VarTTTTERWxxxxxH WHH WH H WH(11-58)如果選取 ，則估計 和 分別為1WR xVarx 111TTxH R HH R z11VarTxH R H(11-59)(11-60)第四節 遞推最小二乘法估計 在前面所討論的最小二乘法和加權最小二乘法，需要同時用到所有的測量數據，在計算時不考慮測量數據的時間順序。當測量數據很多時，要求計算機具有很大的存儲量。在實際處理過程中，測量數據往往是按時間順序逐步給出的，我們可先處理已經得到的一批數據，得到的近似估值，來了新的數據后，再對原估值進行修正，這樣可

31、以減少計算機的存儲量。 下面先討論一維遞推最小二乘法。設觀測值 是一維的，假定已進行了 次觀測，得到 的 個觀測值。用 和 表示相應的向量和矩陣，即zkzkkkzHx、kv 1111211221222212nnkkkkkknkzh ththtzh ththtzh thththhzHh，1122kknrxrxknrx vx，先用加權最小二乘法處理這一批觀測值，加權矩陣 ，一般 為對角線矩陣。1kkWRkR21222000000kkR把觀測方程寫成矩陣形式kkkzH xv(11-61)則 21222100100100kkW處理k個觀測值，可得 的估值 xkx1TTkkkkkkkxH W HH W

32、z(11-62)設 1TkkkkPH W H(11-63)則 TkkkkkxP H W z(11-64)現在又得到了第 次觀測值1k 1kz11 1122111kkknnkkzx h tx htx ht hx式中， 。111211kkknkh ththt h， ，(11-65) 下面通過 個觀測值，求得 的估值 ，然后進一步推出 與 的遞推關系。 1k x1kx1kxkx將式(11-65)與式(11-61)合并，可得111kkkzvHx(11-66)式中 111111kkkkkkkkkvzHvzHvzh，(11-67)由于加權矩陣 一般為對角線矩陣，所以選取W1100kkkwWW(11-68)

33、的估值x1kx11111111TTkkkkkkkxHWHHWz(11-69)利用矩陣求逆引理（附錄二），把求 逆矩陣問題轉變為求階數低的逆矩陣。nn11111111111100kkTTTkkkkkkkTTkkkkkkhwhhwhHWHWHHH W H(11-70)再設 1TkkkkPH W H或 1TkkkkPH W H設 11111TkkkkPHWH(11-71)或 11111TTkkkkkkkhwhPH W H(11-72)根據矩陣求逆引理可得11111111TTTkkkkkkkkkkwPPP hhP hhP(11-73) 這樣，把求 矩陣的逆陣轉變為求標量 的倒數，使計算在為簡化，把式(

34、11-69)寫成nn1111TTkkkkwhP h111111111111111111111111kTTTTTkkkkkkkkkkkkkkkkTTTTkkkkkkkkkkkkkkkTTkkkkkkhhhwhhwhwzP hhhwhhwxPPPPH W zzP H W zPPPH W zz111111111111111111111TTTTkkkkkkkkkkkkkkkkkTTTTkkkkkkkkkkkwwhhhwhhwxxP hzP hhP hP H W zPPPz(11-74)注意到11111111111111111111111111111111TTTTkkkkkkkkkkkkkkkkTTTk

35、kkkkkkkkkkTTkkkkkkkwwwwwwwP hzP hhP hhP hzP hhP hhP hzP hhP hz將式說明(11-75)代入式(11-74)，并考慮到式(11-74)，可得111111111TTkkkkkkkkkkkwzxxP hhP hhx該式說明： 的新估值 是原估值 加上與新的觀測值 和 之差的線性修正項。 按式(11-73)遞推公式計算。x1kxkx1kz1kz1kkhx1kP(11-75)(11-76) 在式(11-73)和式(11-76)中的 最優值為觀測誤差 的方差，可用下式表示：11kw1kv122111kkkwE則式(11-76)代入和式(11-73)變成121111111TTkkkkkkkkkkkzxxP hhP hhx12111111TTkkkkkkkkkkPPP hhP hhP(11-77)(11-78) 按式(11-77)和式(11-78)進行遞推計算，必須知道 和 的初值 和 。初值 和 稱為驗前估計。如果沒有給出 和 ，則可以先解第一批的 個方程，求 逆矩陣，這可大大減少計算工作量。kxkP0P0 x0P0 x0