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

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

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

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

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

6、入向量u(n)可表示為un=un,un-1,un-M+1T如果抽頭輸入向量un和期望響應dn是聯合平穩的，此時均方誤差或者在時刻n的代價函數J(n)是抽頭權向量的二次函數，于是可以得到Jn=d2-wHnp-pHwn+wHnRw(n)(2.7)其中，d2為目標函數dn的方差，p抽頭輸入向量un與期望響應dn的互相關向量，及R為抽頭輸入向量un的相關矩陣。從而梯度向量可以寫為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)分別是代價函

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

8、估值的方法都更加有效。z-1I+v1(n) x(n+1) xn y(n)F(n+1,n) C(n) v2(n)圖3.1 線性動態離散時間系統的信號流圖表示“狀態”的概念是這種表示的基礎。狀態向量，簡單地說狀態，定義為數據的最小集合，這組數據足以唯一地描述系統的自然動態行為。換句話說，狀態由預測系統未來特性時所素要的，與系統的過去行為有關的最少的數據組成。典型地，比較有代表性的情況是，狀態x(n)是未知的。為了估計它，我們使用一組觀測數據，在途中用向量y(n)表示。y(n)成為觀測向量或者簡稱觀測值，并假設它是N維的。在數學上，圖3.1表示的信號流圖隱含著一下兩個方程：(1) 過程方程xn+1=

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

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

11、ynH(m)=O(3.6)將式(3.5)代入式(3.6)，并利用新息過程的正交性質，即得ExiHm=BimEmH(m)=BimR(m)(3.7)因此，式(3.7)兩邊同時右乘逆矩陣R-1(m)，可得Bim的表達式為Bim=ExiHmR-1(m)(3.8)最后，將式(3.8)帶入式(3.5)，可得最小軍方差估計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時刻的狀態x(n+1)與n時刻的狀態x(n)的

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

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

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

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

16、陣N×NK(n,n-1)xnyn-1中誤差相關矩陣M×MK(n)xnyn中誤差相關矩陣M×M基于單步預測的卡爾曼濾波器的小結觀測值= y1,y2,y(n-1)轉移矩陣=F(n+1,n)測量矩陣=C(n)過程噪聲的相關矩陣=Q1(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、這里只討論簡單的一維單輸入單輸出線性系統模型，其中加入白噪聲作為系統的擾動，具體仿真結果可以獲得如下4.1 維納最速下降法濾波器仿真結果以上為最速下降法中不同的遞歸步長所導致的跟蹤效果變化，對于最速下降法中的步長是影響其算法穩定的關鍵，最速下降算法穩定的充分必要條件是條件步長因子為小于輸入自相關矩陣的最大特征值倒數的2倍。上面的序列分別從相關矩陣的隨大特征之2倍的0.4倍開始變化至其1倍，最后一幅圖象能夠看出其已經不再收斂，下面是大于輸入相關矩陣的最大特征值2倍步長時所表現的跟蹤結果可以看出其已經明顯發散，不再是我們所期望的濾波算法。因此可以總結出，對于最速下降法來說，步長的選取是很重要的，根

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

19、曼濾波器仿真這種濾波器只是實現了無需對部分變量進行初始化的設計，沒有特別意義上的改進經典卡爾曼濾波器本身性能的特點。仿真圖如下4.4 后聯平滑濾波的卡爾曼濾波器仿真只是在經典卡爾曼濾波器后端聯接了平滑濾波器，對性能改進的效果并不特別明顯，仿真圖如下如圖中所表示，即使平滑過的估值與觀測值之間的差別也不是特別令人滿意，所以，對于經典卡爾曼濾波的研究還需要更深一步進行，由于時間和能力有限，本次的作業對于卡爾曼及其他濾波器的研究只能達到這種程度，希望在以后的學習中，能發現更好的對經典卡爾曼濾波器的改進方法。5 Matlab源代碼(部分參考自互聯網)5.1經典卡爾曼濾波器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%狀態矩陣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');%紅色卡爾曼，綠色觀測值，藍色狀態值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 后聯平滑濾波器的卡爾曼濾波器clearclc;N=300;CON = 5; x = zeros(1,N);x(1) = 1;p = 10;Q = randn(1,N)*0.2;%過程噪聲協方差R = randn(1,N);%觀測噪聲協方差y = R + CON;%加過程噪聲的狀態輸出for k = 2 : N Q1 = cov(Q(1:k-1);%過程噪聲協方差 Q2 = cov(R(1:k-1); x(k) = x(k - 1);%預估計k時刻狀態變量的值 p = p + Q1;%對應于預估值的協方差 kg = p / (p + Q2);%kalman gain x(k) = x(k) + kg * (y(k) - x