1、管理系統分析專題報告1/30最佳狀態估計理論卡爾曼濾波方法研究摘要：在任何系統中，為了對系統形成有效的控制，系統狀態的準確把握顯得尤為重要。我們可以通過一系列的手段對特定系統進行觀測，以估計出系統的過去、現在、未來的狀態，具體應用可分別表現為對過去狀況的評估、當前狀態的實時控制、趨勢的準確預測等等。最優估計出系統狀態過程中，實際的量測往往是存在諸如來自系統自身、測量工具等所帶來的干擾，控制論中將這種干擾定義為噪聲。如何去尋求濾除這種噪聲干擾，便成為最佳系統狀態估計首先必須解決的問題。Kalman濾波等一些濾波算法便因此應運而生，其作為一種最優估計理論與方法，由于它的實時遞推、存儲量小和簡單易行

2、的特點，在工程應用中受到了重視，廣泛應用于信號處理、控制、通信、航天、制導、目標跟蹤、石油勘探、故障診斷、衛星測控、GPS定位、檢測與估計及機器人等等領域。卡爾曼濾波隨時間及研究的發展，已形成了多種多樣的理論和應用的形式。本次的學習帶著了解認識該濾波算法思想和數學思維的目的，只對一般卡爾曼濾波問題（最優預測與最優濾波）進行了基本的研究。據此本文主要內容安排包括：1、卡爾曼濾波的理論背景；2、理論基礎（選取離散系統的最優預測問題作為代表）闡述；3、工程擴展應用狀況；4、理論局限性等方面。關鍵字：卡爾曼濾波；最優估計；狀態估計；控制2/30目錄1緒論 . 41.1卡爾曼濾波理論研究及應用概述 .

3、41.2維納濾波簡述 . 51.卡爾曼濾波理論概述 . 72卡爾曼濾波理論基礎 . 92.1卡爾曼濾波問題的提法 . 92.2離散系統卡爾曼最優預測基本方程的推導（舉例） . 122.2.1求解基本過程推演. 132.2.2卡爾曼預測估計遞推方程的計算步驟 . 192.2.3應用算例 . 203工程擴展應用舉例 . 223.1卡爾曼濾波在飛機控制中的應用 . 223.2基于卡爾曼濾波方法的時用水量預測（定量描述略） . 264卡爾曼濾波局限性分析 . 284.1穩定性定理 . 284.2濾波的發散問題 . 284.3卡爾曼濾波理論的進一步發展 . 29參考文獻. 303/301緒論1.1卡爾曼

4、濾波理論研究及應用概述在自動控制、航空與航天、通訊、導航和工業生產等領域中，越來越多地遇到“估計”問題。所謂“估計”，簡單地說，就是從觀測數據中提取信息。在飛行器導航中，要從帶有隨機干擾的觀測數據中，估計出飛行器的位置，速度和加速度等運動狀態變量，這就遇到狀態變量的估計問題，這些變量都是隨機過程。因此，“估計”的任務就是從帶有隨機誤差的觀測數據中估計出某些參數或某些狀態變量，這些被估計參數或狀態變量可統稱為被稱為被估量。我們希望估計出來的參數或狀態變量越接近實際值越好。我們希望估計是最好的，因此提出最優估計問題。所謂最優估計，是指在某一確定的準則條件下，按照某種統計意義上來說，估計達到最優。為

5、了正確地解決參數估計和狀態估計問題，首先要研究估計方法。最早的估計方法是高斯（Gauss,K.F.）于1795年在他的關于天體運動理論一書中提出的最小二乘法。最小二乘法沒有考慮被估參數和觀測數據的統計特性，因此這種方法不是最優估計方法。在192年費歇（Fisher,R.A.）提出了極大似然估計方法：從概率密度出發來考慮估計問題，對估計理論做出了重大貢獻。對于隨機過程的估計問題，到本世紀三十年代才積極開展起來。主要成果為1940年美國學者維納（Wiener,N）所提出的在頻域中設計統計最優濾波器的方法，這一方法稱為維納濾波。同時期，蘇聯學者哥爾莫郭洛夫（A.H.）提出并初次解決了離散平穩隨機序列

6、的觀測和外推問題。維納濾波和哥爾莫郭洛夫濾波方法，局限于處理平穩隨機過程，并只能提供穩態的最優估值。這些濾波方法在工程實踐上也遇到許多困難，在實際應用上受到一定的限制。1960年卡爾曼（Kalman,R.E.）和布西（Bucy,R.S）提出了最優遞推濾波方法，可用數字計算機來實現。卡爾曼濾波既適用于平穩隨機過程，又適用于非平穩隨機過程，因此卡爾曼濾波方法得到廣泛的應用。本文要學習的卡爾曼濾波理論，正是源于他的博士論文和1960年發表的論文A New Approach to Linear Filtering and Prediction Problems（線性濾波與預測問題的新方法）。卡爾曼濾波

7、在工程實踐中，特別在航空空間技術中迅速得到應用。例如在測4/30軌問題和慣性導航等方面都應用卡爾曼濾波。所謂測軌問題，就是在不同時刻，對飛行器進行觀測，根據觀測數據，應用卡爾曼濾波方法，估計出這個飛行器每時每刻的狀態變量，如飛行器的位置、速度、加速度及阻力系數等物理量。以便對飛行器進行導航、制導和攔截。斯坦利.施密特(Stanley Schmidt)首次實現了卡爾曼濾波器。卡爾曼在NASA埃姆斯研究中心訪問時，發現他的方法對于解決阿波羅計劃的軌道預測很有用，后來阿波羅飛船的導航電腦便使用了這種濾波器。關于這種濾波器的論文由Swerling(1958)、Kalman(1960)與Kalman a

8、nd Bucy(1961)發表。在其它方面卡爾曼濾波也得到應用。目前,卡爾曼濾波已經有很多不同的實現。卡爾曼最初提出的形式現在一般稱為簡單卡爾曼濾波器。除此以外，還有施密特擴展濾波器、信息濾波器以及很多Bierman，Thornton開發的平方根濾波器的變種。也許最常見的卡爾曼濾波器是鎖相環，它在收音機、計算機和幾乎任何視頻或通訊設備中廣泛存在。當前對于卡爾曼濾波較前沿的研究與應用有很多，諸如機器人導航、控制；傳感器數據融合；雷達系統以及導彈追蹤；計算機圖像處理，包括頭臉識別、圖像分割、圖像邊緣檢測等。1.2維納濾波簡述為了實現最優估計，最佳線性濾波理論起源于40年代美國科學家Wiener和前

9、蘇聯科學家等人的研究工作，后人統稱為維納濾波理論。這是在信號和干擾都表示成有理譜密度的情況下，找出最優濾波器，使得實際輸出與希望輸出之間的均方誤差最小。維納濾波問題的關鍵是導出維納霍夫積分方程，解這一積分方程可得最優濾波器的脈沖過渡函數，從脈沖過渡函數可得濾波器的傳遞函數。通常，解維納霍夫積分方程是很困難的，即使對少數情況，能得解析解，但在工程上也往往難以實現。特別對于有限記憶或非常平穩過程，維納濾波問題變得更為復雜了。維納濾波問題的提法簡單羅列如下：設系統的觀測方程為z(t)=x(t)+v(t)；式中，x(t)為有用信號z(t)為觀測信號5/30v(t)為觀測誤差。設x(t)、z(t)和v(

10、t)都是均值為零并具有各態歷經性的平穩隨機過程。(t)。維納濾波的任務就是根據觀測值z(t)估計x(t)，使x(t)估值接近于x設計出一個線性定常系統L，如圖2-1所示，使得系統的輸出y(t)與x(t)具有最小方差，即J=Ex(t)-y(t)=min 2 (1-3)(t)。 這樣y(t)就作為x(t)的估值x如果系統的脈沖過渡函數為h(l)，則y(t)=0h(l)z(t-l)dl (1-4)y(t)是系統L根據輸入信號z(t)在(-，t)上的全部過去值所給出的實際輸出，如圖2-2所示。y(t)是z(t-l)的線性函數(l0)。根據問題的性質，維納濾波有下列三個條件：信號與噪聲都是均值為零并具有

11、各態歷經性的平穩隨機過程；濾波器是一個物理可實現的線性定常系統。當lt稱為預測（或外推）問題第二類：t1=t，稱為濾波問題；第三類：t1k，稱為預測（或外推）問題；第二類：j=k，稱為濾波問題；第三類：jk稱為平滑（或內插）問題。本文這里只討論連續系統和離散系統的最優預測和最優濾波問題。并通過簡述離散系統卡爾曼濾波最優預測基本方程的推導來說明卡爾曼濾波的具體過程，其余類型過程基本相似，在此只描述它們的推演處理方法及思想，詳細過程略去。2.2離散系統卡爾曼最優預測基本方程的推導（舉例）在推導卡爾曼預測基本方程時，為了簡便起見，先不考慮控制信號的作用，這樣，離散系統的差分方程(2-11)變成x(k

12、+1)=F(k+1,k)x(k)+G(k+1,k)w(k)觀測方程仍為式(2-2)： (2-16)z(k)=H(k)x(k)+v(k)式中，w(k)和v(k)都均值為零的白噪聲序列，w(k)和v(k)相互獨立。 在采樣間隔內w(k)和v(k)為常值，其統計特性如式(2-13)所示，即Ew(k)=Ev(k)=0Ew(k)wT(j)=QdkkjEv(k)vT(j)=RdkkjEw(k)vT(j)=0 狀態向量的初值x(0)，其統計特性是給定的，即Ex(0)=m0TEx(0)-mxt-m=P0 000()給出觀測序列z(0)，z(1)，L，z(k)，要求找出x(k+1)的線性最優估計(k+1/k)，

13、使得估計誤差x%(k+1/k)=x(k+1)-x(k+1/k)的方差為最小，x即12/30T(k+!/k)Exk+1-xk+1/k()()x(k+1)-x=min (k+1/k)是z(0)，z(1)，要求估值xL，z(k)的線性函數，并且要求估計是無偏的，即(k+1/k)=EExx(k+1)2.2.1求解基本過程推演下面推導卡爾曼預測基本公式。推導的方法有幾種，比較簡易的方法是利用正交定理，用數學歸納法推導卡爾曼估計的基本遞推估計公式。當獲得觀測值z(0)，z(1)，L，z(k-1)之后，假定已經得到狀態向量x(k)的(k/k-1)。當還未獲得k時刻的新觀測值z(k)時，根據一個最優線性預測估

14、計x已有的觀測值，可得k+1時刻的系統狀態向量x(k+1)的兩步預測估值(k+1/k-1)： x(k+1/k-1)=F(k+1,k)x(k/k-1) x (2-17)(k/k-1)是x(k)的一步最優線性估計，x(k+1/k-1)也是x(k+1)的由于x最優線性預測估計，這可用正交定理來證明。由式(2-16)減式(2-17)，可得%(k+1/k-1)=x(k+1)-x(k+1/k-1)x(k/k-1)=F(k+1,k)x(k)-x+G(k+1,k)w(k)%(k+1/k-1)=F(k+1,k)x%(k/k-1)+G(k+1,k)w(k) x%(k/k-1)=x(k)-x(k/k-1) 式中，x

15、 (2-18)(k/k-1)是x(k)的最優線性預測估值，%(k/k-1)因為x根據正交定理，估計誤差x%(k/k-1)的線性變換必須正交于z(0)，z(1)，L，z(k-1)，所以x%(k/k-1)也必須正交于z(0)，F(k+1,k)xL，z(k-1)。式(2-18)中的w(k)是均值為零的白噪聲序列，與z(0)，L，z(k-1)相互獨立，因此w(k)正交于(k+1/k-1)是x(k+1)的最優兩z(0)，L，z(k-1)。所以在沒有獲得之前z(k)，x步線性預測。13/30(k+1/k-1)來得到x(k+1)的在新觀測值z(k)獲取之后，可通過修正兩步估值x一步預測估值x(k+1/k)。

16、設z(k)的預測估值(k/k-1) z(k/k-1)=H(k)x由式(2-2) (2-19)z(k)=H(k)x(k)+v(k)減式(2-19)，得z(k)的預測估計誤差為%(k/k-1)=z(k)-z(k/k-1)z(k/k-1)=H(k)x(k)-x+v(k)%(k/k-1)+v(k)=H(k)x%(k/k-1)的原因有兩個：造成z對k時刻狀態向量x(k)的預測估計有誤差；附加了v(k)時刻的觀測噪聲干擾。%(k/k-1)包含修正x(k+1/k-1)的新信息。這樣，在獲得之后，在兩步估顯然z(k+1/k-1)的基礎上，用z%(k/k-1)去修正，便可得到k+1時刻狀態向量值x(k+1/k)

17、，即 x(k+1)的一步預測估計x(k+1/k)=x(k+1/k-1)+K(k)z%(k/k-1) x或(k+1/k)=F(k+1,k)x(k/k-1)+K(k)(k/k-1)xz(k)-H(k)x(2-20) 式中K(k)是待定矩陣，稱為最優增益矩陣或加權矩陣。 式(2-20)可改寫成(k+1/k)=F(k+1,k)x(k/k-1)x%(k/k-1)+K(k)v(k)+K(k)H(k)x (2-21)k+1時刻系統狀態方程為x(k+1)=F(k+1,k)x(k)+G(k+1,k)w(k) (2-22)由式(2-22)減式(2-21)得x(k+1)的估計誤差為14/30%(k+1/k)=%(k

18、/k-1)xF(k+1/k)-K(k)H(k)x+G(k+1/k)w(k)-K(k)v(k) (2-23) %(k/k-1)、w(k)、v(k)均分別正交于觀察式(2-23)右邊，x%(k+1/k)()正交于z(0)，z(1)，z(0)，z(1)，L，z(k-1)，因此，xL，z(k-1)。%(k+1/k)與z(k)正交，則x(k+1/k)就是x(k+1)的一步最優線性預測估值。若x%(k+1/k)與z(k)的正交條件： 因此，利用x%(k+1/k)zT(k)=0 Ex (2-24)來確定最優增益矩陣K(k)。把式(2-23)和式(2-2)代入式(2-24)，得%(k/k-1)+G(k+1/k

19、)w(k)-K(k)v(k)E【F(k+1/k)-K(k)H(k)xT(k/k-1)+H(k)x%(k/k-1)+v(k)】H(k)x=0(k/k-1)、w(k)和v(k)相互間都是正交的，因此上式可簡考慮到(k/k-1)、x化為%(k/k-1)x%T(k/k-1)HT(k)-K(k)v(k)vT(k)=0EF(k+1/k)-K(k)H(k)x即Tk/k-1HTk-KkEvkvTk=0x%Fk+1/k-KkHkEk/k-1x()()()()()()()()()%(k/k-1)x%T(k/k-1)=P(k/k-1)，又有Ev(k)vT(k)=R，代入上式后設Exk可得最優增益矩陣為-1TT K(

20、k)=F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)+R(k)(2-25)現在需確定估計誤差方差陣P(k+1/k)的遞推關系式。 根據估計誤差方差陣的定義，有%(k+1/k)x%T(k+1/k) P(k+1/k)=Ex將式(2-23)代入上式得15/30%T(k/k-1)P(k+1/k)=E【F(k+1,k)-K(k)H(k)xT%(k/k-1)+G(k+1,k)w(k)-K(k)v(k)Fk+1,k-KkHk()()()xT+G(k+1/k)w(k)-K(k)v(k)】%(k/k-1)、w(k)和v(k)相互間都正交，可得： 考慮到xP(k+1/k)=F(k+1,k

21、)-K(k)H(k)P(k/k-1)F(k+1,k)-K(k)H(k)T+K(k)RKT(k)+G(k+1,k)QGT(k+1/k)kk式中Q=Ew(k)wT(k) k將上式展開，經整理后得P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-K(k)H(k)P(k/k-1)FT(k+1/k)-F(k+1,k)P(k/k-1)HT(k)KT(k)+K(k)H(k)P(k/k-1)HT(k)KT(k)+K(k)RKT(k)k+G(k+1,k)QGT(k+1/k)k上式右邊第四項與第五項之和為 (2-26)K(k)H(k)P(k/k-1)HT(k)KT(k)+K(k)RKT(k)k

22、=K(k)H(k)P(k/k-1)HT(k)+RKT(k)k-1TT=F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)+RkH(k)P(k/k-1)HT(k)+RKT(k)k=Fk+1,kPk/k-1HTkKTk()()()()顯然，式(2-26)右邊第四項與第五項之和在數值上等于第三項，但符號相反。這 樣，式(2-25)右邊第三、第四和第五項之和為零，所以P(k-1/k)=F(k+1/k)P(k/k-1)FT(k+1,k)-K(k)H(k)P(k/k-1)FT(k+1/k)+G(k+1/k)QGT(k+1/k)k16/30把式(2-25)的K(k)代入上式，可得估計誤

23、差方差陣的遞推關系式為P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-1TT(2-F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)-RkH(k)P(k/k-1)FT(k+1,k)+G(k+1,k)QGT(k+1,k)k27)方程(2-20)、(2-25)和(2-27)構成一組完整的最優線性估計方程，現綜合如下： 最優預測估計方程為式(2-20)，即(k+1/k)=F(k+1,k)x(k/k-1)+K(k)(k/k-1)xz(k)-H(k)x最優增益矩陣方程為(2-25)，即-1TT K(k)=F(k+1,k)P(k/k-1)H(k)H(k)P(k/

24、k-1)H(k)+Rk估計誤差方差陣的遞推方程為式(2-27)，即P(k+1/k)=F(k+1,k)P(k/k-1)FT(k+1/k)-1TT -F(k+1,k)P(k/k-1)H(k)H(k)P(k/k-1)H(k)-RkH(k)P(k/k-1)FT(k+1,k)+G(k+1,k)QGT(k+1,k)k從式(2-27)可看出，估計誤差方差陣與Q和R有關，而與觀測值z(k)無關。因kk此，可事先估計誤差方差陣P(k+1/k)，同時也可算出最優增益矩陣K(k)。 按式(2-16)和式(2-2)作出系統模型方塊圖，如圖2-3所示。圖2-4表示由方程(2-21)所描述的卡爾曼最優預測估計方塊圖。17

25、/30從圖2-4可看出，最優預測估計由三部分組成：估值為觀測值的線性函數；最優增益矩K(k)；單位負反饋回路。現在需要驗證第二節時最優估計提出的三項要求：估值為觀測值的線性函數；估計值的方差為最小；估值是無偏的，即(k+1/k)Ex=Ex(k+1)在上面推導預測估計方程時，是按照和兩項要求推導的，現需要說明估值是無偏的問題。對式(2-16)的兩端取數學期望，考慮到Ew(k)=0，可得Ex(k+1)=F(k+1,k)Ex(k) (2-28)再對式(2-21)的兩端取數學期望，考慮到Ev(k)=0，可得(k+1/k)(k/k-1)%(k/k-1)Ex=F(k+1,k)Ex+K(k)H(k)Ex(k

26、/k-1)(k/k-1)=F(k+1,k)Ex+K(k)H(k)Ex(k)-x將式(2-28)減式(2-29)，得到 (2-29)(k+1/k)(k/k-1)Ex(k+1)-x=F(k+1,k)Ex(k)-x(k/k-1)-K(k)H(k)Ex(k)-x如果初始條件為 (2-30)(0/0)Ex=Ex(0)=m0或18/300/0x則 (-)=Ex(0)=m0()%0/0=E0/0=0 Exx(0)-Ex-根據式(2-30)的遞推關系，可得(1/0)Ex(1)-x=0(2/1)Ex(2)-x=0M(k+1/k)Ex(k+1)-x=0()因此(k+1/k)Ex=Ex(k+1)k=0,1,2L0/

27、0=E所以，只要初始估計選為Exx(0)，所得估計是無偏的。 -2.2.2卡爾曼預測估計遞推方程的計算步驟卡爾曼預測估計遞推方程的計算步驟如下：在t時刻給定初值： 0估值誤差方差陣初值：P0/0()0/0x(-(0)=E)=xx(0)=m0 (T=P(0)=Ex0-x0x0-x0()()()() -根據公式(2-25)計算t時刻最優增益陣K(0)： 0-1TT K(0)=F(1，0)P(0)H(0)H(0)P(0)H(0)+R0(1/0)： 根據公式(2-20)計算x(1)的最優估值x(1/0)=F(1，(0)+K(0)(0)x0)xz(0)-H(0)x根據公式(2-27)計算P(1/0)：-

28、1TTTP(1/0)=F(1，0)P(0)F(1，0)-F(1，0)P(0)H(0)H(0)P(0)H(0)+R0H(0)P(0)FT(1，0)+G(1，0)QG(1，0)019/30根據已知的P(1/0)計算t時刻的K(1)。 1(2/1)。 根據K(1)計算x(2)的估值x重復上述遞推計算步驟，可得(3/2)，(k/k-1)，P(k/k-1),K(k)，x(k+1/k) P(2/1),K(2),xL，x2.2.3應用算例例2-1設系統狀態方程和觀測方程為x(k+1)=0.5x(k)+w(t)z(k)=x(k)+v(k)w(k)和v(k)都是均值為零的白噪聲序列，且不相關，其統計特性如下：E

29、w(k)=0,Ev(k)=0Ew(k)w(j)=1d初值Ex(0)=m0=0,kj,Ev(k)v(j)=2d kjP0/0(-)=1觀測值z(0)=0，z(1)=4，z(2)=2試求x(k)的最優預測估值。解:與前面的有關方程對照，可得F(k+1,k)=0.5，G(k+1,k)=1H(k)=1，Q(k)=1，Rk=2k最優預測估值方程為(k+1/k)=0.5x(k/k-1)+K(k)(k/k-1)xz(k)-x最優增益矩陣為K(k)=0.5P(k/k-1)P(k/k-1)+2估值誤差方差陣遞推方程為 -12-1P(k+1/k)=0.25P(k/k-1)-0.25Pk/k-1Pk/k-1+2()

30、()+1(1/0)，x(2/1)和x(3/2) 下面計算x20/30P0/0(-)=116K(0)=0.51(1+2)-1=1P(1/0)=0.25-0.25+1=1.1673K1=0.51.1671.167+2-1=0.1842 ()()P(2/1)=0.251.167-0.251.167(1.167+2)-1+1=1.1842K(2)=0.51.1842(1.1842+2)-1=0.18590/0取x(-)的初值0/0x(-)=Ex(0)=m0=01(1/0)=0.50+0=0x6(2/1)=0.50+0.1842(4-0)=0.7368x(3/2)=0.50.7368+0.1859(2-

31、0.7368)=0.6032x2.2.4其它系統估計問題基本方程1.連續系統卡爾曼最優濾波2.離散系統卡爾曼最優濾波3.噪聲為有色噪聲情況4.系統噪聲與觀測噪聲相關情況與離散系統卡爾曼最優預測推導過程和基本思想相似，最終得出卡爾曼濾波的三個重要方程。由于時間限制，這里不再詳細闡述。21/303工程擴展應用舉例3.1卡爾曼濾波在飛機控制中的應用卡爾曼濾波在航空和航天中，以及在其他工程領域中都得到應用。對于每一個具體應用問題，都要求深入了解問題的物理實質和工程實際問題，因此比較復雜。這里僅對卡爾曼濾波在飛機自動駕駛中的經典應用進行了探討。圖3-1自動駕駛儀（附注：自動駕駛儀是按技術要求自動控制飛行

32、器軌跡的調節設備，其作用主要是保持飛機姿態和輔助駕駛員操縱飛機。對無人駕駛飛機，它將與其他導航設備配合完成規定的飛行任務。導彈上的自動駕駛儀起穩定導彈姿態的作用，故稱導彈姿態控制系統。自動駕駛儀是模仿駕駛員的動作駕駛飛機的。它由敏感元件、計算機和伺服機構組成。當飛機偏離原有姿態時，敏感元件檢測變化，計算機算出修正舵偏量，伺服機構將舵面操縱到所需位置。）22/30圖3-2飛機縱向運動參量關系圖設飛機在垂直平面內運動（縱向運動）各參量的關系如圖3-2所示。 設O為飛機重心，OX1為飛機縱軸，v為飛機速度向量。OX1軸與地平面的夾角J稱為俯仰角。速度向量v與地平面的夾角q稱為航跡角。a稱為攻角，d為

33、飛機操縱面的偏轉角。J、q和a成下列關系： rrJ=q+a(3-1)飛機的縱向運動方程為&+aJ&+aa=ad(3-2) J123&-a&+a4a=0(3-3) J飛機在飛行過程中往往受到擾動氣流的作用，可把這一擾動作用歸結為在攻角a中出現了附加的隨機干擾項w1，則方程(3-2)和式(3-3)變成&+aJ&J1+a2a=a3d-a2w1(3-4)&-a&-a4a=a1w1(3-5) J把上述二式寫成狀態方程，設狀態變量為&，x=J(3-6) x1=a，x2=J3則狀態方程&1=-a4x1+x2-a1w1(3-7) x&2=-a2x1-a1x2+a3d-a2w1(3-8) x&3=x2(3-9)

34、 x23/30把上述方程組寫成矩陣形式：&=Ax+Bu-Tw(3-10) x式中x=x1，x2，x3，u=d，w=w1T-a1A=-a201-a1100a 0，B=300&，由于有測量誤在飛機上用二自由度測速陀螺儀測量飛機的俯仰角速度J差，測速陀螺的輸出為&+v=x+n z1=J121(3-11)&，由于有測量誤差，三自由度在飛機上用三自由度陀螺測量飛機的俯仰角J陀螺的輸出為z2=J+v2=x3+v2(3-12)把式(3-11)和式(3-12)聯合起來得到矩陣形式z=Hx+v 式中(3-13)z1010v1z=，H=，v=v z00122設w1，v1，v2都是均值為零的白噪聲，且互不相關，即E

35、w1=Ev1=Ev2=02Ewtwt=qd(t-t) ()()111r12Ev(t)v(t)=0T0d(t-t) 2r2TEwtv()(t)1=0&和q的估值。 求a、J解根據前面所述的連續系統卡爾曼濾波基本方程（）三式，可得最優濾波方程為24/30&(t/t)x1-a1&(t/t)=-ax22&0xt/t)3(1-a111(t/t)00xk112(t/t)+a3u+k210x(t/t)00xk3131(t/t)xk12z1010k22-xt/t)z0012(2k32xt/t()3(3-14)最優增益矩陣為k11K=k21k31k12p11pk22=21k32p31p12p22p32p1300

36、-210r1p230p3301p12r1-20-2=pr221r2-2p32r1-2p13r2-2p23r2-2p33r2-2(3-15)&=0時，可得代數黎卡提方程： 當濾波達到穩態時，即當P-a110p11p12p13p11-a-a0p+ppp2121222321100p31p32p33p31-a1p11p12q2-a-a0-p+-a21p2222110p31p32p11010001p21p31p12p22p32p13p23=0p33p13-a1-a201-a1ap22p231p32p33000p1300-r1-20p2310-20-r201p33p12(3-16)P為對稱矩陣，因此p12

37、=p21，p13=p31，p32=p23設a1=1.76秒，a2=4.2秒，a3=7.4秒，a4=0.77秒q1=1.2度，r1=0.85度/秒，r2=0.85度求解式(3-16)得p11=0.86，p22=2.78/秒，p33=0.48 p31=p13=0.41秒，p21=p12=0.58，p23=p32=0.17求解式(3-15)得k1=0.81秒，k12=0.57，k21=3.8625/30k22=0.24秒，k21=0.23秒，k32=0.67&和q的估值。 把求得的K代入式(3-14)，可得a、J&=x&=x1(t/t)，J2(t/t)，J=x3(t/t)a3.2基于卡爾曼濾波方法的

38、時用水量預測（定量描述略）杭州位于錢塘江畔,是一個人口密集城市,其設計供水量能力為140萬t/d。本文以杭州市的時用水量為例,分別采用卡爾曼濾波模型以及卡爾曼濾波季節模型進行預測,對結果進行比較分析(見二圖),研究上述2種方法的利弊。從圖中可以看出,對于1h預測,非季節模型誤差有較大波動,改進后的模型,即26/30卡爾曼濾波季節模型相對較為平穩。總體來說季節模型預測精度相對較高。評價模型預測精度,還需查看其預測誤差分布特征7。圖3為卡爾曼濾波常規模型與季節模型預測誤差分布比較,采用了1個月數據,用2種方法進行了預測,從圖中可以看出,2種方法誤差都主要集中在-0.1至0.1區域,呈對稱分布。從1

39、個月所有數據來看,誤差大致呈正態分布,季節模型較常規模型趨勢圖更為尖銳,表明季節模型比常規模型預測效果更為理想。綜合上圖和下圖,比較得出,卡爾曼濾波季節模型具有更高的預測精度。采用基于卡爾曼濾波方法預測模型,進行時用水量預測。鑒于城市時用水量具有周期性、趨勢性及隨機擾動性,引入季節因子,形成卡爾曼濾波方法季節模型。以杭州市實測數據為例,分別采用卡爾曼濾波模型和卡爾曼濾波季節模型分別對該市時用水量進行預測,預測結果表明,卡爾曼濾波方法季節預測模型精度相對較高。從計算模型看,季節預測模型要求低,但預測精度較高,具有一定實用性,能夠滿足供水系統優化調度的要求。27/304卡爾曼濾波局限性分析4.1穩定性定理前面比較詳細地討論了線性系統卡爾曼濾波基本方程的推導和遞推計算步驟。在估值計算時，需要利用一連串的觀測數據，按照濾波基本方程進步遞推計(t1/t)。這里，需要確切知道x(t)的初值和算，可得狀態向量x(k)的最優估值x估計誤差方差陣的初值P(0/0)。但在許多實際問題中，往往不可能確切知道初值x(0)和P(0/0)，甚至根本不知道這些初值。為了進步濾波計算，只能假定初