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

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

3、約束最優化的數學表示。特別適合于自適應濾波的一類無約束最優化算法基于局部迭代下降的算法：從某一初始猜想出發，產生一系列權向量，使得代價函數在算法的每一次迭代都是下降的，即其中是權向量的過去值，而是其更新值。我們希望算法最終收斂到最優值。迭代下降的一種簡單形式是最速下降法，該方法是沿最速下降方向連續調整權向量。為方便起見，我們將梯度向量表示為(2.2)因此，最速下降法可以表示為(2.3)其中代表進程，是正常數，稱為步長參數，1/2因子的引入是為了數學上處理方便。在從到的迭代中，權向量的調整量為(2.4)為了證明最速下降算法滿足式(2.1)，在處進行一階泰勒展開，得到(2.5)此式對于較小時是成立

4、的。在式(2.4)中設為負值向量，因而梯度向量也為負值向量，所以使用埃爾米特轉置。將式(2.4)用到式(2.5)中，得到此式表明當為正數時，。因此，隨著的增加，代價函數減小，當時，代價函數趨于最小值。2.2最速下降算法應用于維納濾波器考慮一個橫向濾波器，其抽頭輸入為，對應的抽頭權值為。抽頭輸入是來自零均值、相關矩陣為的廣義平穩隨機過程的抽樣值。除了這些輸入外，濾波器還要一個期望響應，以便為最優濾波提供一個參考。在時刻抽頭輸入向量表示為，濾波器輸出端期望響應的估計值為,其中是由抽頭輸所張成的空間。空過比較期望響應及其估計值，可以得到一個估計誤差，即(2.6)這里是抽頭權向量與抽頭輸入向量的內積。

5、可以進一步表示為同樣，抽頭輸入向量可表示為如果抽頭輸入向量和期望響應是聯合平穩的，此時均方誤差或者在時刻的代價函數是抽頭權向量的二次函數，于是可以得到(2.7)其中，為目標函數的方差，抽頭輸入向量與期望響應的互相關向量，及為抽頭輸入向量的相關矩陣。從而梯度向量可以寫為(2.8)其中在列向量中和分別是代價函數對應第個抽頭權值的實部和虛部的偏導數。對最速下降算法應用而言，假設式(2.8)中相關矩陣和互相關向量已知，則對于給定的抽頭權向量為(2.9)它描述了為那濾波中最速下降法的數學表達式。3.卡爾曼濾波器3.1卡爾曼濾波器的基本思想卡爾曼濾波器是用狀態空間概念描述其數學公式的，另外新穎的特點是，他

6、的解遞歸運算，可以不加修改地應用于平穩和非平穩環境。尤其是，其狀態的每一次更新估計都由前一次估計和新的輸入數據計算得到，因此只需存儲前一次估計。除了不需要存儲過去的所有觀測數據外，卡爾曼濾波計算比直接根據濾波過程中每一步所有過去數據進行估值的方法都更加有效。+ 圖3.1 線性動態離散時間系統的信號流圖表示“狀態”的概念是這種表示的基礎。狀態向量，簡單地說狀態，定義為數據的最小集合，這組數據足以唯一地描述系統的自然動態行為。換句話說，狀態由預測系統未來特性時所素要的，與系統的過去行為有關的最少的數據組成。典型地，比較有代表性的情況是，狀態是未知的。為了估計它，我們使用一組觀測數據，在途中用向量表

7、示。成為觀測向量或者簡稱觀測值，并假設它是維的。在數學上，圖3.1表示的信號流圖隱含著一下兩個方程：(1) 過程方程(3.1)式中，向量表示噪聲過程，可建模為零均值的白噪聲過程，且其相關矩陣定義為(2) 測量方程(3.2)其中是已知的測量矩陣。向量稱為測量噪聲，建模為零均值的白噪聲過程，其相關矩陣為(3.3)測量方程(3.2)確立了可測系統輸出與狀態之間的關系，如圖3.1所示。3.2 新息過程為了求解卡爾曼濾波問題，我們將應用基于新息過程的方法。根據之前所述，用向量表示時刻到時刻所有觀測數據過去值給定的情況下，你時刻觀測數據的最小均方估計。過去的值用觀測值表示，他們張成的向量空間用表示。從而可

8、以定義新息過程如下：(3.4)其中向量表示觀測數據的新息。3.3 應用新息過程進行狀態估計下面，我們根據信息過程導出狀態的最小均方估計。根據推導，這個估計可以表示成為新息過程序列的線性組合，即(3.5)其中是一組待定的矩陣。根據正交性原理，預測狀態誤差向量與新息過程正交，即(3.6)將式(3.5)代入式(3.6)，并利用新息過程的正交性質，即得(3.7)因此，式(3.7)兩邊同時右乘逆矩陣，可得的表達式為(3.8)最后，將式(3.8)帶入式(3.5)，可得最小軍方差估計(3.9)故對于，有(3.10)然而，時刻的狀態與時刻的狀態的關系式由式可以推導出對于，有(3.11)其中只與觀測數據有關。因

9、此可知，與彼此正交(其中)。利用式(3.11)以及當時的計算公式，可將式(3.10)右邊的求和項改寫為(3.12)為了進一步討論，引入如下基本定義。3.4 卡爾曼增益定義矩陣(3.13)其中是狀態向量和新息過程的互相關矩陣。利用這一定義和式(3.12)的結果，可以將式(3.10)簡單重寫為(3.14)式(3.14)具有明確的物理意義。它標明：線性動態系統狀態的最小均方估計可以由前一個估計求得。為了表示對卡爾曼開創性貢獻的認可，將矩陣稱為卡爾曼增益。現在剩下唯一要解決的問題是，怎樣以一種便于計算的形式來表示卡爾曼增益。為此，首先將與乘積的期望表示為(3.15)式中利用了狀態與噪聲向量互不相關這一

10、事實。其次，由于預測狀態誤差向量與估計正交，因此與乘機的期望為零。這樣，用預測狀態誤差向量代替相乘因子，將不會引起式(3.15)變化，故有(3.16)由此，可將上式進一步變化為(3.17)現在我們重新定義卡爾曼增益。為此，將式(3.17)代入式(3.13)得(3.18)現在我們已經了解了卡爾曼濾波的整個過程和相應的參數設置，為了能夠更為方便利用計算機仿真實現，特將其中參數變量進行小結。卡爾曼變量和參數小結變量定義維數時刻狀態時刻狀態值從時刻到時刻的轉移矩陣時刻的測量矩陣過程噪聲的相關矩陣過程噪聲的相關矩陣給定觀測值在時刻狀態的預測估計給定觀測值在時刻狀態的濾波估計時刻卡爾曼增益矩陣時刻新息向量

11、新息向量的相關矩陣中誤差相關矩陣中誤差相關矩陣基于單步預測的卡爾曼濾波器的小結觀測值=轉移矩陣=測量矩陣=過程噪聲的相關矩陣=測量噪聲的相關矩陣= 4 Matlab仿真為了簡化，這里只討論簡單的一維單輸入單輸出線性系統模型，其中加入白噪聲作為系統的擾動，具體仿真結果可以獲得如下4.1 維納最速下降法濾波器仿真結果以上為最速下降法中不同的遞歸步長所導致的跟蹤效果變化，對于最速下降法中的步長是影響其算法穩定的關鍵，最速下降算法穩定的充分必要條件是條件步長因子為小于輸入自相關矩陣的最大特征值倒數的2倍。上面的序列分別從相關矩陣的隨大特征之2倍的0.4倍開始變化至其1倍，最后一幅圖象能夠看出其已經不再

12、收斂，下面是大于輸入相關矩陣的最大特征值2倍步長時所表現的跟蹤結果可以看出其已經明顯發散，不再是我們所期望的濾波算法。因此可以總結出，對于最速下降法來說，步長的選取是很重要的，根據不同條件的需求，選取正確的步長，能為算法的快速高效提供基礎。4.2 卡爾曼濾波器仿真結果從圖中可以發現，卡爾曼濾波器能夠非常有效地在比較大的干擾下比較準確地反映真實值，如果觀測端加入干擾較大時，卡爾曼濾波器能夠較為有效地進行濾波，不過當狀態端的干擾增大時，卡爾曼濾波器的濾波效果也會隨之下降。如下圖，是加大了狀態端的干擾，所呈現的濾波效果。如上圖所示，狀態端的干擾導致狀態不穩定，卡爾曼濾波器的估計值也出現了比較大的波動

13、。如果將狀態端的干擾再增大，則會出現更為嚴峻的濾波考驗，濾波效果如下。這是的狀態已經很勉強了，所以，研究更為有效的多方法卡爾曼濾波器也顯得十分必要了。4.3 一種不需初始化的卡爾曼濾波器仿真這種濾波器只是實現了無需對部分變量進行初始化的設計，沒有特別意義上的改進經典卡爾曼濾波器本身性能的特點。仿真圖如下4.4 后聯平滑濾波的卡爾曼濾波器仿真只是在經典卡爾曼濾波器后端聯接了平滑濾波器，對性能改進的效果并不特別明顯，仿真圖如下如圖中所表示，即使平滑過的估值與觀測值之間的差別也不是特別令人滿意，所以，對于經典卡爾曼濾波的研究還需要更深一步進行，由于時間和能力有限，本次的作業對于卡爾曼及其他濾波器的研

14、究只能達到這種程度，希望在以后的學習中，能發現更好的對經典卡爾曼濾波器的改進方法。5 Matlab源代碼(部分參考自互聯網)5.1經典卡爾曼濾波器clearN=200;w(1)=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

15、;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(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&&am

16、p;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(:,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 后聯平滑濾波器的卡爾曼濾波器clearc

17、lc;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(k); p

18、 = (1 - kg) * p;endFilter_Wid = 10;smooth_res = zeros(1,N);kalman_p = zeros(1,N);for i = Filter_Wid + 1 : N tempsum = 0; kalman_m = 0; for j = i - Filter_Wid : i - 1 tempsum = tempsum + y(j); kalman_m = kalman_m+x(j); end kalman_p(i) = kalman_m/Filter_Wid; smooth_res(i) = tempsum / Filter_Wid;%平滑濾波end% figure(1);% hist(y);t=1:N;figure(1);expValue = zeros(1,N);for i = 1: N expValue(i) = CON;endplot(t,expValue,'r',t,x,'g',t,y,'b