1、信息融合大作業(yè)維納最速下降法濾波器，卡爾曼濾波器設(shè)計及 Matlab 仿真 時間： 2010-12-6 專業(yè)：信息工程 班級： 09030702 姓名：馬志強1. 濾波問題淺談估計器或濾波器這一術(shù)語通常用來稱呼一個系統(tǒng)， 設(shè)計這樣的系統(tǒng)是為了從 含有噪聲的數(shù)據(jù)中提取人們感興趣的，接近規(guī)定質(zhì)量的信息。由于這樣一個寬 目標，估計理論應(yīng)用于諸如通信、 雷達、聲納、導(dǎo)航、地震學(xué)、生物醫(yī)學(xué)工程、 金融工程等眾多不同的領(lǐng)域。例如，考慮一個數(shù)字通信系統(tǒng)，其基本形式由發(fā) 射機、信道和接收機連接組成。 發(fā)射機的作用是把數(shù)字源 （例如計算機 ）產(chǎn)生的 0、 1符號序列組成的消息信號變換成為適合于信道上傳送的波形。

2、而由于符號間 干擾和噪聲的存在，信道輸出端收到的信號是含有噪聲的或失真的發(fā)送信號。 接收機的作用是，操作接收信號并把原消息信號的一個可靠估值傳遞給系統(tǒng)輸 出端的某個用戶。隨著通信系統(tǒng)復(fù)雜度的提高，對原消息信號的還原成為通信 系統(tǒng)中最為重要的環(huán)節(jié)，而噪聲是接收端需要排除的最主要的干擾，人們也設(shè) 計出了針對各種不同條件應(yīng)用的濾波器，其中最速下降算法是一種古老的最優(yōu) 化技術(shù)，而卡爾曼濾波器隨著應(yīng)用條件的精簡成為了普適性的高效濾波器。 2維納最速下降算法濾波器2.1 最速下降算法的基本思想考慮一個代價函數(shù) ?（?，） 它是某個未知向量 ?的連續(xù)可微分函數(shù)。 函數(shù) ?（?） 將 ?的元素映射為實數(shù)。這里

3、，我們要尋找一個最優(yōu)解 ?。使它滿足如下條件?(? < ?(?)(2.1)這也是無約束最優(yōu)化的數(shù)學(xué)表示特別適合于自適應(yīng)濾波的一類無約束最優(yōu)化算法基于局部迭代下降的算法： 從某一初始猜想？?(0)出發(fā)，產(chǎn)生一系列權(quán)向量？21),?22),使得代價函數(shù)?(?在算法的每一次迭代都是下降的，即?(? 1)< ?(?(?)其中?(?是權(quán)向量的過去值，而？2?打1)是其更新值。我們希望算法最終收斂到最優(yōu)值？。迭代下降的一種簡單形式是最速下降 法，該方法是沿最速下降方向連續(xù)調(diào)整權(quán)向量。為方便起見，我們將梯度向量 表示為?= ?凸=竺?)?(2.2)因此，最速下降法可以表示為1?s(？s+ 1)

4、= ?_?(?)2(2.3)其中？代表進程，？是正常數(shù)，稱為步長參數(shù)，1/2因子的引入是為了數(shù)學(xué)上處 理方便。在從？到？?+ 1的迭代中，權(quán)向量的調(diào)整量為1?=?s(？s+ 1) - ?7?= -?(?)2(2.4)為了證明最速下降算法滿足式(2.1)，在？?處進行一階泰勒展開，得到?(?？+ 1)?(?) + ?(?(?)(2.5)此式對于？較小時是成立的。在式(2.4)中設(shè)？?為負值向量，因而梯度向量?也為負值向量，所以使用埃爾米特轉(zhuǎn)置。將式(2.4)用到式(2.5)中，得到1?(?+ 1)?(?e?)-?(?)i2此式表明當？為正數(shù)時，？?(?+ 1) < ?(?(?)因此，隨著？

5、的增加,代價函數(shù)?減小，當？?= X時，代價函數(shù)趨于最小值 也。2.2最速下降算法應(yīng)用于維納濾波器考慮一個橫向濾波器，其抽頭輸入為？?,？?？ 1),?(? ?+1),對應(yīng)的 抽頭權(quán)值為？(？?,？(？?,？?-1 (? 0抽頭輸入是來自零均值、相關(guān)矩陣為？?勺 廣義平穩(wěn)隨機過程的抽樣值。除了這些輸入外，濾波器還要一個期望響應(yīng)？?(?) 以便為最優(yōu)濾波提供一個參考。在時刻 ？抽頭輸入向量表示為？?(?)濾波器輸出端期望響應(yīng)的估計值為其中？?是由抽頭輸？?,？? 1),?(?+ 1)所張成的空間。空過比較期望響應(yīng)?(?及其估計值，可以得到一個估計 誤差?(?)即?= ?- ? = ?- ?(?(

6、?)(2.6)這里？?(？?(?是抽頭權(quán)向量?(?與抽頭輸入向量?(?的內(nèi)積。？?(?可以進- 步表示為?= ?(?,?(?,?-1 (?同樣，抽頭輸入向量?(?可表示為? = ?,? 1),? ?+ 1) ?如果抽頭輸入向量？?和期望響應(yīng)？?是聯(lián)合平穩(wěn)的，此時均方誤差或者 在時刻？的代價函數(shù)？?(?是抽頭權(quán)向量的二次函數(shù)，于是可以得到? = ?- ?b?(?2 ?+ ?b彳？(？)(2.7)其中，？2為目標函數(shù)?的方差，？抽頭輸入向量？?與期望響應(yīng)?的互相關(guān)向量，及？為抽頭輸入向量？?？的相關(guān)矩陣。從而梯度向量可以寫為?(?)? ?(?)?(?)?(?)?(?) ?(?) ? = ?(?)?

7、(?)?-1 (?)?(?)?-1 (?) ?(?) = -2? + 2?(?)(2.8)其中在列向量中籌和鴛分別是代價函數(shù)？?對應(yīng)第？個抽頭權(quán)值？(？?的實?)部？?V?和虛部？?(?的偏導(dǎo)數(shù)。對最速下降算法應(yīng)用而言，假設(shè)式(2.8)中相關(guān) 矩陣？和互相關(guān)向量？已知，則對于給定的抽頭權(quán)向量？2?+ 1)為?+ 1) = ? ? + ? ?1 / /L/(2.9)它描述了為那濾波中最速下降法的數(shù)學(xué)表達式。3.卡爾曼濾波器3.1卡爾曼濾波器的基本思想卡爾曼濾波器是用狀態(tài)空間概念描述其數(shù)學(xué)公式的，另外新穎的特點是， 他的解遞歸運算，可以不加修改地應(yīng)用于平穩(wěn)和非平穩(wěn)環(huán)境。尤其是，其狀態(tài) 的每一次更新

8、估計都由前一次估計和新的輸入數(shù)據(jù)計算得到，因此只需存儲前 一次估計。除了不需要存儲過去的所有觀測數(shù)據(jù)外，卡爾曼濾波計算比直接根 據(jù)濾波過程中每一步所有過去數(shù)據(jù)進行估值的方法都更加有效。?(?)?(? 1)?y(?)O二二三?二0=® 0O才圖3.1線性動態(tài)離散時間系統(tǒng)的信號流圖表示“狀態(tài)”的概念是這種表示的基礎(chǔ)：。狀態(tài)向量，簡單地說狀態(tài)，定義為數(shù)據(jù)的最小集合，這組數(shù)據(jù)足以唯一地描述系統(tǒng)的自然動態(tài)行為。 換句話說，狀態(tài)由預(yù)測系統(tǒng)未來特性時所素要的，與系統(tǒng)的過去行為有關(guān)的最少的數(shù)據(jù)組成。典型地，比較有代表性的情況是，狀態(tài)？?(?是未知的。為了估計它，我們使用一組觀測數(shù)據(jù)，在途中用向量？?

9、(?表示。y(?成為觀測向量或者簡稱觀測值，并假設(shè)它 是？?維的。在數(shù)學(xué)上，圖3.1表示的信號流圖隱含著一下兩個方程:(1)過程方程?+ 1) = ?+ 1, ? + ?(?)(3.1)式中，M X1向量？?(?康示噪聲過程，可建模為零均值的白噪聲過程，且其 相關(guān)矩陣定義為?(? = ?；? ? ?；(2)測量方程? = ? + ?<(?)(3.2)其中？?是已知的n XM測量矩陣。N X1向量？?稱為測量噪聲，建模 為零均值的白噪聲過程，其相關(guān)矩陣為?(? = ? ? ?為？(3.3) 測量方程 (3.2) 確立了可測系統(tǒng)輸出 ?(?)與狀態(tài)?(?)之間的關(guān)系，如圖 3.1 所示。3.

10、2 新息過程為了求解卡爾曼濾波問題， 我們將應(yīng)用基于新息過程的方法。 根據(jù)之前所述， 用向量?(?|?-1)表示?= 1時刻到?- 1時刻所有觀測數(shù)據(jù)過去值給定的情況下， 你時刻觀測數(shù)據(jù) ? (?的)最小均方估計。 過去的值用觀測值 ?(1),?(2),?(?-? 1)表 示，他們張成的向量空間用 ?-1表示。從而可以定義新息過程如下：?(?) = ?(?)- ?(?|?-1)(3.4) 其中M X1向量？?表示觀測數(shù)據(jù)?的新息。3.3 應(yīng)用新息過程進行狀態(tài)估計下面，我們根據(jù)信息過程導(dǎo)出狀態(tài) ?(?的?)最小均方估計。根據(jù)推導(dǎo)，這個估 計可以表示成為新息過程 ?(1),?(2),?(?)序列的

11、線性組合，即?)=刀？?？?(?)?=1(3.5) 其中 ?(?)?=1是一組待定的 ?X ?矩陣。根據(jù)正交性原理，預(yù)測狀態(tài)誤差向量 與新息過程正交，即?(?,?)?(?) = ?(?)?- ?|?(?) = ?(3.6) 將式(3.5)代入式 (3.6)，并利用新息過程的正交性質(zhì)，即得?(?)?(?) = ?(?)?(?)?(?) = ?(?)?(?)(3.7) 因此，式 (3.7)兩邊同時右乘逆矩陣 ?-1 (?)，可得 ?(?)的表達式為? ?(?) = ? ?( ?)?(?) ?-1 (?)(3.8) 最后，將式 (3.8)帶入式 (3.5)，可得最小軍方差估計?l?)=刀?/p?)?

12、(?7 ?(?)?=1?-1=刀？?(? ?9(?)?+ ?(? ?9(?)? ?=1(3.9) 故對于 ?= ?+ 1，有?-1?+ 1|?)=刀？?+ 1)?/?(?) ?(?)?=1+ ?(?+ 1)?(?)?-1 (?)?(?)(3.10) 然而， ?+ 1 時刻的狀態(tài) ?(?+? 1) 與?時刻的狀態(tài) ?(?的) 關(guān)系式由式可以推導(dǎo)出對 于0 < k < n,有?(?+ 1)?(?) = ?(?+ 1,?)?(?) + ?1(?)?(?)= ?(?+ 1,?)?(?)?(?)(3.11) 其中?(?)只與觀測數(shù)據(jù) ?(1),?(2),?(?有) 關(guān)。因此可知, ?1(?)

13、與?(?)彼此正交 (其中0 < k < n)。利用式(3.11)以及當？= ?時？?？?)的計算公式，可將式(3.10) 右邊的求和項改寫為?-1?-1刀 ?+ 1)?/?(?| ? (?)?(? = ?+ 1, ?刀 ？?(？?| 7?1 (?=1?=1= ?(?+ 1 , ?)?(?| ?-1 )(3.12) 為了進一步討論，引入如下基本定義。3.4 卡爾曼增益定義M XN矩陣? = ?(?1)?(?) ? (?)(3.13) 其中？?(?71)?(?)是狀態(tài)向量?(?+ 1)和新息過程?(?的互相關(guān)矩陣。利用 這一定義和式(3.12)的結(jié)果，可以將式(3.10)簡單重寫為?

14、+ 1|?) = ?+ 1, ?-1) + ?(?)(3.14)式(3.14)具有明確的物理意義。它標明：線性動態(tài)系統(tǒng)狀態(tài)的最小均方估計?+ 1|?詢可以由前一個估計？?-1)求得。為了表示對卡爾曼開創(chuàng)性貢獻的 認可，將矩陣？?稱為卡爾曼增益。現(xiàn)在剩下唯一要解決的問題是，怎樣以一種便于計算的形式來表示卡爾曼增 益？?。為此，首先將?(?+ 1)與?:?乘積的期望表示為?+ 1)? = ?+ 1, ?軟？?？ 1) ?(?)(3.15)式中利用了狀態(tài)？?與噪聲向量？?(n)互不相關(guān)這一事實。其次，由于預(yù)測狀態(tài) 誤差向量?(? 1)與估計？?-1)正交，因此？?-1)與？(？? 1)乘機的 期望為

15、零。這樣，用預(yù)測狀態(tài)誤差向量 ?(?- 1)代替相乘因子？?？，將不會引 起式(3.15)變化，故有?+ 1) ?(? = ?+ 1, ?(? 1)?(? 1) ?)(3.16)由此，可將上式進一步變化為?- 1)?(?刃 =?+ 1, ?(? 1)?(?)(3.17)現(xiàn)在我們重新定義卡爾曼增益。為此，將式(3.17)代入式(3.13)得? = ?» 1 ?(? 1)?夕7?) ? (?)1 ，,、， 丿、/ /(3.14) 現(xiàn)在我們已經(jīng)了解了卡爾曼濾波的整個過程和相應(yīng)的參數(shù)設(shè)置，為了能夠更為方 便利用計算機仿真實現(xiàn)，特將其中參數(shù)變量進行小結(jié)。卡爾曼變量和參數(shù)小結(jié)?(?)?時刻狀態(tài)?

16、x 1?(?)?時刻狀態(tài)值?X1?(?+ 1, ?)從？時刻到？?+ 1時刻的轉(zhuǎn)移矩陣?x ?(?)?時刻的測量矩陣?x?欣?)過程噪聲？?4?的相關(guān)矩陣?x ?><?)過程噪聲？賢?的相關(guān)矩陣?x ?-1)給定觀測值?1),?2),?(? 1)?x 1在？時刻狀態(tài)的預(yù)測估計?)給定觀測值?1),?2),?(?在?x 1時刻狀態(tài)的濾波估計?(?)?時刻卡爾曼增益矩陣? x ?(?)?時刻新息向量?x1?(?)新息向量?(?的相關(guān)矩陣?x ?(? 1)?,)中誤差相關(guān)矩陣?x ?(?)?)中誤差相關(guān)矩陣?x ?基于單步預(yù)測的卡爾曼濾波器的小結(jié)觀測值=?1),?2),?(? 1)轉(zhuǎn)移矩

17、陣=?(?+ 1,?)測量矩陣=?(?)過程噪聲的相關(guān)矩陣=?)測量噪聲的相關(guān)矩陣=?)? ? = ?+ 1, ? ? 1) ?( ?2 ? ? 1) ?( ? + ?>?)-1?= ?- ?(?/)?+ 1| ?) = ?+ 1, ?-1) + ?(?)? = ?；? 1) - ?+ 1)?(?7?2 1)?+ 1, ? = ?+ 1, ?(?+ 1,?+ ?(?)4 Matlab 仿真為了簡化，這里只討論簡單的一維單輸入一單輸出線性系統(tǒng)模型，其中加入白噪聲作為系統(tǒng)的擾動，具體仿真結(jié)果可以獲得如下4.1 維納最速下降法濾波器仿真結(jié)果 以上為最速下降法中不同的遞歸步長所導(dǎo)致的跟蹤效果變化

18、， 對于最速下降法中 的步長是影響其算法穩(wěn)定的關(guān)鍵， 最速下降算法穩(wěn)定的充分必要條件是條件步長 因子為小于輸入自相關(guān)矩陣的最大特征值倒數(shù)的 2 倍。上面的序列分別從相關(guān)矩 陣的隨大特征之 2 倍的 0.4 倍開始變化至其 1 倍，最后一幅圖象能夠看出其已經(jīng) 不再收斂，下面是大于輸入相關(guān)矩陣的最大特征值 2 倍步長時所表現(xiàn)的跟蹤結(jié)果 可以看出其已經(jīng)明顯發(fā)散， 不再是我們所期望的濾波算法。 因此可以總結(jié)出， 對 于最速下降法來說， 步長的選取是很重要的， 根據(jù)不同條件的需求， 選取正確的 步長，能為算法的快速高效提供基礎(chǔ)。4.2 卡爾曼濾波器仿真結(jié)果 從圖中可以發(fā)現(xiàn)，卡爾曼濾波器能夠非常有效地在比

19、較大的干擾下比較準確 地反映真實值， 如果觀測端加入干擾較大時， 卡爾曼濾波器能夠較為有效地進行 濾波，不過當狀態(tài)端的干擾增大時， 卡爾曼濾波器的濾波效果也會隨之下降。 如 下圖，是加大了狀態(tài)端的干擾，所呈現(xiàn)的濾波效果。如上圖所示， 狀態(tài)端的干擾導(dǎo)致狀態(tài)不穩(wěn)定， 卡爾曼濾波器的估計值也出現(xiàn) 了比較大的波動。如果將狀態(tài)端的干擾再增大，則會出現(xiàn)更為嚴峻的濾波考驗， 濾波效果如下。這是的狀態(tài)已經(jīng)很勉強了， 所以，研究更為有效的多方法卡爾曼濾波器也顯 得十分必要了。4.3 一種不需初始化的卡爾曼濾波器仿真 這種濾波器只是實現(xiàn)了無需對部分變量進行初始化的設(shè)計， 沒有特別意義上 的改進經(jīng)典卡爾曼濾波器本身

20、性能的特點。仿真圖如下4.4 后聯(lián)平滑濾波的卡爾曼濾波器仿真 只是在經(jīng)典卡爾曼濾波器后端聯(lián)接了平滑濾波器， 對性能改進的效果并不特 別明顯，仿真圖如下如圖中所表示，即使平滑過的估值與觀測值之間的差別也不是特別令人滿意， 所以，對于經(jīng)典卡爾曼濾波的研究還需要更深一步進行，由于時間和能力有限，希望在以后的學(xué)本次的作業(yè)對于卡爾曼及其他濾波器的研究只能達到這種程度， 習(xí)中，能發(fā)現(xiàn)更好的對經(jīng)典卡爾曼濾波器的改進方法。5 Matlab 源代碼 ( 部分參考自互聯(lián)網(wǎng) )5.1 經(jīng)典卡爾曼濾波器clearN=200;w(1)=0;x(1)=5;a=1;c=1;Q1 = randn (1,N)*1;%過程噪聲Q

21、2 = randn (1,N);%測量噪聲for k=2:N;x(k)=a*x(k-1)+Q1(k-1); en d%狀 態(tài)矩陣for k=1:N;Y(k)=c*x(k)+Q2(k);endp(1)=10;s(1)=1;for t=2:N;Rww=cov(Q1(1:t);Rvv=cov(Q2(1:t);p1(t)=a42*p(t-1)+Rww;b(t)=c*p1(t)/(c42*p1(t)+Rvv);%kalman 增益s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1);p(t)=p1(t)-c*b(t)*p1(t);endt=1:N;plot(t,s,'r'

22、;,t,Y,'g',t,x,'b');%紅色卡爾曼，綠色觀測值，藍色狀態(tài)值legend('kalman estimate','ovservations','truth');5.2 最速下降法clcclear allN=30;q=2.1;%q>1 &&q<2/Ryx 最大特征值 hn=zeros(1,N);hn(:)=5;vg=0;Rxx=xcorr(1);Ryx=min(min(corrcoef(1, 1+randn); echo offfor i=1:N-1; %vg=2*Rxx*hn

23、(:,i)-2*Ryx;%hn(:,i+1)=hn(:,i)-1/2*q*vg; vg=2*Rxx*hn(i)-2*Ryx;hn(i+1)=hn(i)-1/2*q*vg; m(i)=1;endt=1:N-1;plot(t,hn(t),'r-',t,m(t),'b-');5.3 后聯(lián)平滑濾波器的卡爾曼濾波器clearclc;N=300;CON = 5;x = zeros(1,N);x(1) = 1;p = 10;Q = randn( 1,N)*0.2;%過程噪聲協(xié)方差R = randn (1,N);%E測噪聲協(xié)方差 y = R + CON;加過程噪聲的狀態(tài)輸出for k = 2 : NQ1 = cov(Q(1:k-1);%過程噪聲協(xié)方差Q2 = cov(R(1:k-1);x(k) = x(k - 1);%預(yù)估計 k 時刻狀態(tài)變量的值 p = p + Q1;%寸應(yīng)于預(yù)估值的協(xié)方差 kg = p / (p + Q2);%kalman gain x(