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

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

3、)J(w)(2.1)這也是無約束最優(yōu)化的數(shù)學(xué)表示。特別適合于自適應(yīng)濾波的一類無約束最優(yōu)化算法基于局部迭代下降的算法：從某一初始猜想w(0)出發(fā)，產(chǎn)生一系列權(quán)向量w1,w2,，使得代價(jià)函數(shù)J(w)在算法的每一次迭代都是下降的，即Jwn+1<J(w(n)其中w(n)是權(quán)向量的過去值，而wn+1是其更新值。我們希望算法最終收斂到最優(yōu)值w0。迭代下降的一種簡單形式是最速下降法，該方法是沿最速下降方向連續(xù)調(diào)整權(quán)向量。為方便起見，我們將梯度向量表示為g=Jw=J(w)w(2.2)因此，最速下降法可以表示為wn+1=wn-12g(n)(2.3)其中n代表進(jìn)程，是正常數(shù)，稱為步長參數(shù)，1/2因子的引入是

4、為了數(shù)學(xué)上處理方便。在從n到n+1的迭代中，權(quán)向量的調(diào)整量為wn=wn+1-wn=-12g(n)(2.4)為了證明最速下降算法滿足式(2.1)，在wn處進(jìn)行一階泰勒展開，得到Jwn+1Jwn+gHnw(n)(2.5)此式對于較小時是成立的。在式(2.4)中設(shè)w為負(fù)值向量，因而梯度向量g也為負(fù)值向量，所以使用埃爾米特轉(zhuǎn)置。將式(2.4)用到式(2.5)中，得到Jwn+1Jwn-12g(n)2此式表明當(dāng)為正數(shù)時，Jwn+1<J(w(n)。因此，隨著n的增加，代價(jià)函數(shù)Jn減小，當(dāng)n=時，代價(jià)函數(shù)趨于最小值J。2.2最速下降算法應(yīng)用于維納濾波器考慮一個橫向?yàn)V波器，其抽頭輸入為un,un-1,u(

5、n-M+1)，對應(yīng)的抽頭權(quán)值為w0n,w1n,wM-1n。抽頭輸入是來自零均值、相關(guān)矩陣為R的廣義平穩(wěn)隨機(jī)過程的抽樣值。除了這些輸入外，濾波器還要一個期望響應(yīng)d(n)，以便為最優(yōu)濾波提供一個參考。在時刻n抽頭輸入向量表示為u(n)，濾波器輸出端期望響應(yīng)的估計(jì)值為d(n|Un),其中Un是由抽頭輸un,un-1,u(n-M+1)所張成的空間。空過比較期望響應(yīng)d(n)及其估計(jì)值，可以得到一個估計(jì)誤差e(n)，即en=dn-dnUn=dn-wHnu(n)(2.6)這里wHnu(n)是抽頭權(quán)向量w(n)與抽頭輸入向量u(n)的內(nèi)積。w(n)可以進(jìn)一步表示為wn=w0n,w1n,wM-1nT同樣，抽頭輸

6、入向量u(n)可表示為un=un,un-1,un-M+1T如果抽頭輸入向量un和期望響應(yīng)dn是聯(lián)合平穩(wěn)的，此時均方誤差或者在時刻n的代價(jià)函數(shù)J(n)是抽頭權(quán)向量的二次函數(shù)，于是可以得到Jn=d2-wHnp-pHwn+wHnRw(n)(2.7)其中，d2為目標(biāo)函數(shù)dn的方差，p抽頭輸入向量un與期望響應(yīng)dn的互相關(guān)向量，及R為抽頭輸入向量un的相關(guān)矩陣。從而梯度向量可以寫為Jn=J(n)a0(n)+jJ(n)b0(n)J(n)a1(n)J(n)aM-1(n)+jJ(n)b1(n)J(n)bM-1(n)=-2p+2Rw(n)(2.8)其中在列向量中J(n)ak(n)和J(n)bk(n)分別是代價(jià)函

7、數(shù)Jn對應(yīng)第k個抽頭權(quán)值wkn的實(shí)部ak(n)和虛部bk(n)的偏導(dǎo)數(shù)。對最速下降算法應(yīng)用而言，假設(shè)式(2.8)中相關(guān)矩陣R和互相關(guān)向量p已知，則對于給定的抽頭權(quán)向量wn+1為wn+1=wn+p-Rwn(2.9)它描述了為那濾波中最速下降法的數(shù)學(xué)表達(dá)式。3.卡爾曼濾波器3.1卡爾曼濾波器的基本思想卡爾曼濾波器是用狀態(tài)空間概念描述其數(shù)學(xué)公式的，另外新穎的特點(diǎn)是，他的解遞歸運(yùn)算，可以不加修改地應(yīng)用于平穩(wěn)和非平穩(wěn)環(huán)境。尤其是，其狀態(tài)的每一次更新估計(jì)都由前一次估計(jì)和新的輸入數(shù)據(jù)計(jì)算得到，因此只需存儲前一次估計(jì)。除了不需要存儲過去的所有觀測數(shù)據(jù)外，卡爾曼濾波計(jì)算比直接根據(jù)濾波過程中每一步所有過去數(shù)據(jù)進(jìn)行

8、估值的方法都更加有效。z-1I+v1(n) x(n+1) xn y(n)F(n+1,n) C(n) v2(n)圖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)x(n)是未知的。為了估計(jì)它，我們使用一組觀測數(shù)據(jù)，在途中用向量y(n)表示。y(n)成為觀測向量或者簡稱觀測值，并假設(shè)它是N維的。在數(shù)學(xué)上，圖3.1表示的信號流圖隱含著一下兩個方程：(1) 過程方程xn+1=

9、Fn+1,nxn+v1(n)(3.1)式中，M×1向量v1(n)表示噪聲過程，可建模為零均值的白噪聲過程，且其相關(guān)矩陣定義為Ev1nv1Hk=Q1n n=kO nk(2) 測量方程yn=Cnxn+v2(n)(3.2)其中Cn是已知的N×M測量矩陣。N×1向量v2(n)稱為測量噪聲，建模為零均值的白噪聲過程，其相關(guān)矩陣為Ev2nv2Hk=Q2n n=kO nk(3.3)測量方程(3.2)確立了可測系統(tǒng)輸出yn與狀態(tài)xn之間的關(guān)系，如圖3.1所示。3.2 新息過程為了求解卡爾曼濾波問題，我們將應(yīng)用基于新息過程的方法。根據(jù)之前所述，用向量y(n|yn-1)表示n=1時刻

10、到n-1時刻所有觀測數(shù)據(jù)過去值給定的情況下，你時刻觀測數(shù)據(jù)y(n)的最小均方估計(jì)。過去的值用觀測值y1,y2,y(n-1)表示，他們張成的向量空間用yn-1表示。從而可以定義新息過程如下：n=yn-y(n|yn-1)(3.4)其中M×1向量n表示觀測數(shù)據(jù)yn的新息。3.3 應(yīng)用新息過程進(jìn)行狀態(tài)估計(jì)下面，我們根據(jù)信息過程導(dǎo)出狀態(tài)x(i)的最小均方估計(jì)。根據(jù)推導(dǎo)，這個估計(jì)可以表示成為新息過程1,2,n序列的線性組合，即xiyn=k=1nBik(k)(3.5)其中Bi(k)k=1n是一組待定的M×N矩陣。根據(jù)正交性原理，預(yù)測狀態(tài)誤差向量與新息過程正交，即Ei,nHm=Exi-xi

11、ynH(m)=O(3.6)將式(3.5)代入式(3.6)，并利用新息過程的正交性質(zhì)，即得ExiHm=BimEmH(m)=BimR(m)(3.7)因此，式(3.7)兩邊同時右乘逆矩陣R-1(m)，可得Bim的表達(dá)式為Bim=ExiHmR-1(m)(3.8)最后，將式(3.8)帶入式(3.5)，可得最小軍方差估計(jì)xiyn=k=1nExiHkR-1(k)k=k=1n-1ExiHkR-1(k)k+ExiHnR-1(n)n(3.9)故對于i=n+1，有xn+1yn=k=1n-1Exn+1HkR-1(k)k+Exn+1HnR-1(n)n(3.10)然而，n+1時刻的狀態(tài)x(n+1)與n時刻的狀態(tài)x(n)的

12、關(guān)系式由式可以推導(dǎo)出對于0kn，有Exn+1Hk=EFn+1,nxn+v1nHk=Fn+1,nExnH(k)(3.11)其中k只與觀測數(shù)據(jù)y1,y2,y(k)有關(guān)。因此可知，v1n與k彼此正交(其中0kn)。利用式(3.11)以及當(dāng)i=n時xiyn的計(jì)算公式，可將式(3.10)右邊的求和項(xiàng)改寫為k=1n-1Exn+1HkR-1(k)k=Fn+1,nk=1n-1ExnHkR-1kk=F(n+1,n) xnyn-1(3.12)為了進(jìn)一步討論，引入如下基本定義。3.4 卡爾曼增益定義M×N矩陣Gn=Ex(n+1)H(n) R-1(k)(3.13)其中Ex(n+1)H(n)是狀態(tài)向量x(n+

13、1)和新息過程(n)的互相關(guān)矩陣。利用這一定義和式(3.12)的結(jié)果，可以將式(3.10)簡單重寫為xn+1yn= Fn+1,nxnyn-1+Gn(n)(3.14)式(3.14)具有明確的物理意義。它標(biāo)明：線性動態(tài)系統(tǒng)狀態(tài)的最小均方估計(jì)xn+1yn可以由前一個估計(jì)xnyn-1求得。為了表示對卡爾曼開創(chuàng)性貢獻(xiàn)的認(rèn)可，將矩陣Gn稱為卡爾曼增益。現(xiàn)在剩下唯一要解決的問題是，怎樣以一種便于計(jì)算的形式來表示卡爾曼增益Gn。為此，首先將x(n+1)與H(n)乘積的期望表示為Exn+1Hn=Fn+1,nExnHn,n-1CH(n)(3.15)式中利用了狀態(tài)xn與噪聲向量v2(n)互不相關(guān)這一事實(shí)。其次，由于

14、預(yù)測狀態(tài)誤差向量(n,n-1)與估計(jì)xnyn-1正交，因此xnyn-1與Hn,n-1乘機(jī)的期望為零。這樣，用預(yù)測狀態(tài)誤差向量(n,n-1)代替相乘因子xn，將不會引起式(3.15)變化，故有Exn+1Hn=Fn+1,nE(n,n-1)Hn,n-1CH(n)(3.16)由此，可將上式進(jìn)一步變化為Exn+1Hn=Fn+1,nK(n,n-1)CH(n)(3.17)現(xiàn)在我們重新定義卡爾曼增益。為此，將式(3.17)代入式(3.13)得Gn=Fn+1,nK(n,n-1)CH(n) R-1(k)(3.18)現(xiàn)在我們已經(jīng)了解了卡爾曼濾波的整個過程和相應(yīng)的參數(shù)設(shè)置，為了能夠更為方便利用計(jì)算機(jī)仿真實(shí)現(xiàn)，特將其中

15、參數(shù)變量進(jìn)行小結(jié)。卡爾曼變量和參數(shù)小結(jié)變量定義維數(shù)x(n)n時刻狀態(tài)M×1y(n)n時刻狀態(tài)值N×1F(n+1,n)從n時刻到n+1時刻的轉(zhuǎn)移矩陣M×MC(n)n時刻的測量矩陣N×MQ1(n)過程噪聲v1(n)的相關(guān)矩陣M×MQ2(n)過程噪聲v2(n)的相關(guān)矩陣N×Nxnyn-1給定觀測值y1,y2,y(n-1)在n時刻狀態(tài)的預(yù)測估計(jì)M×1xnyn給定觀測值y1,y2,y(n)在n時刻狀態(tài)的濾波估計(jì)M×1G(n)n時刻卡爾曼增益矩陣M×N(n)n時刻新息向量N×1R(n)新息向量(n)的相關(guān)矩

16、陣N×NK(n,n-1)xnyn-1中誤差相關(guān)矩陣M×MK(n)xnyn中誤差相關(guān)矩陣M×M基于單步預(yù)測的卡爾曼濾波器的小結(jié)觀測值= y1,y2,y(n-1)轉(zhuǎn)移矩陣=F(n+1,n)測量矩陣=C(n)過程噪聲的相關(guān)矩陣=Q1(n)測量噪聲的相關(guān)矩陣=Q2(n) Gn=Fn+1,nKn,n-1CHnCnKn,n-1CHn+Q2(n)-1 n=yn-C(n)xnyn-1 xn+1yn= Fn+1,nxnyn-1+Gn(n) Kn=Kn,n-1-Fn,n+1GnCnK(n,n-1) Kn+1,n=Fn+1,nKnFHn+1,n+Q1(n)4 Matlab仿真為了簡化，

17、這里只討論簡單的一維單輸入單輸出線性系統(tǒng)模型，其中加入白噪聲作為系統(tǒng)的擾動，具體仿真結(jié)果可以獲得如下4.1 維納最速下降法濾波器仿真結(jié)果以上為最速下降法中不同的遞歸步長所導(dǎo)致的跟蹤效果變化，對于最速下降法中的步長是影響其算法穩(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é)出，對于最速下降法來說，步長的選取是很重要的，根

18、據(jù)不同條件的需求，選取正確的步長，能為算法的快速高效提供基礎(chǔ)。4.2 卡爾曼濾波器仿真結(jié)果從圖中可以發(fā)現(xiàn)，卡爾曼濾波器能夠非常有效地在比較大的干擾下比較準(zhǔn)確地反映真實(shí)值，如果觀測端加入干擾較大時，卡爾曼濾波器能夠較為有效地進(jìn)行濾波，不過當(dāng)狀態(tài)端的干擾增大時，卡爾曼濾波器的濾波效果也會隨之下降。如下圖，是加大了狀態(tài)端的干擾，所呈現(xiàn)的濾波效果。如上圖所示，狀態(tài)端的干擾導(dǎo)致狀態(tài)不穩(wěn)定，卡爾曼濾波器的估計(jì)值也出現(xiàn)了比較大的波動。如果將狀態(tài)端的干擾再增大，則會出現(xiàn)更為嚴(yán)峻的濾波考驗(yàn)，濾波效果如下。這是的狀態(tài)已經(jīng)很勉強(qiáng)了，所以，研究更為有效的多方法卡爾曼濾波器也顯得十分必要了。4.3 一種不需初始化的卡爾

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

20、0;x(1)=5;a=1;c=1;Q1 = randn(1,N)*1;%過程噪聲Q2 = randn(1,N);%測量噪聲for k=2:N;x(k)=a*x(k-1)+Q1(k-1); end%狀態(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)=a.2*p(t-1)+Rww;b(t)=c*p1(t)/(c.2*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

21、(t)*p1(t);endt=1:N;plot(t,s,'r',t,Y,'g',t,x,'b');%紅色卡爾曼，綠色觀測值，藍(lán)色狀態(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)

22、;echo offfor i=1:N-1; %vg=2*Rxx*hn(:,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);%觀測噪聲協(xié)方差y = R + CON;%加過程噪聲的狀態(tài)輸出for k = 2 : N Q1 = cov(Q(1:k-1);%過程噪聲協(xié)方差 Q2 = cov(R(1:k-1); x(k) = x(k - 1);%預(yù)估計(jì)k時刻狀態(tài)變量的值 p = p + Q1;%對應(yīng)于預(yù)估值的協(xié)方差 kg = p / (p + Q2);%kalman gain x(k) = x(k) + kg * (y(k) - x