維納濾波與卡爾曼濾波

第四章維納濾波和卡爾曼濾波

引言維納濾波器的離散形式——時(shí)域解離散維納濾波器的z域解維納預(yù)測(cè)卡爾曼(Kalman)濾波

引言信道s(n)x(n)s(n)：原始輸入（發(fā)射）信號(hào)，隨機(jī)平穩(wěn)信道噪聲（測(cè)量噪聲）x(n)：接收（測(cè)量）信號(hào)，隨機(jī)平穩(wěn)問題提出：準(zhǔn)則：最大后驗(yàn)準(zhǔn)則。均方準(zhǔn)則，最大似然準(zhǔn)則，濾波器

h(n)x(n)s(n)FIR,IIR逼近（準(zhǔn)則）已知：s(n),(n)的統(tǒng)計(jì)特性，要求：設(shè)計(jì)線性移不變?yōu)V波器h(n),從x(n)中恢復(fù)s(n)線性均方準(zhǔn)則（最小二乘濾波）M個(gè)權(quán)系數(shù)（抽頭）的橫向?yàn)V波器

定義：：輸入信號(hào)：輸入向量2.橫向?yàn)V波器結(jié)構(gòu)

:濾波器的權(quán)系數(shù)

:濾波器權(quán)向量

:期望響應(yīng):對(duì)期望響應(yīng)的估計(jì):估計(jì)誤差假設(shè)由信號(hào)與噪聲組成⑴如果，圖1的系統(tǒng)稱為濾波（filtering）；⑵如果，圖1的系統(tǒng)稱為預(yù)測(cè)（prediction）；⑶如果，圖1的系統(tǒng)稱為平滑（smoothing）。圖期望信號(hào)、估計(jì)值與誤差信號(hào)的幾何關(guān)系圖表明在濾波器處于最佳工作狀態(tài)時(shí)，估計(jì)值加上估計(jì)偏差等于期望信號(hào)，即注意我們所研究的是隨機(jī)信號(hào)，圖中各矢量的幾何表示應(yīng)理解為相應(yīng)量的統(tǒng)計(jì)平均或者是數(shù)學(xué)期望。再?gòu)哪芰康慕嵌葋砜?，假定輸入信?hào)和期望信號(hào)都是零均值,應(yīng)用正交性原理，則 ,因此在濾波器處于最佳狀態(tài)時(shí)，估計(jì)值的能量總是小于等于期望信號(hào)的能量。4.2維納濾波器的離散形式—時(shí)域解

最小均方誤差準(zhǔn)則：（加性干擾）為此令e(n)和x(n-j)正交（且無關(guān)）是h最優(yōu)的必要充分條件，即任何時(shí)刻估計(jì)誤差都與輸入數(shù)據(jù)正交。由正交方程可得：一、維納—霍夫方程定義可得維納—霍夫（Wiener-Holf)方程或標(biāo)準(zhǔn)方程求和范圍（i）隨濾波器的不同取不同區(qū)間一、維納—霍夫方程定義式可以寫成矩陣的形式，即對(duì)上式求逆，得到

上式表明已知期望信號(hào)與觀測(cè)數(shù)據(jù)的互相關(guān)函數(shù)及觀測(cè)數(shù)據(jù)的自相關(guān)函數(shù)時(shí)，可以通過矩陣求逆運(yùn)算，得到維納濾波器的最佳解。同時(shí)可以看到，直接從時(shí)域求解因果的維納濾波器，當(dāng)選擇的濾波器的長(zhǎng)度M較大時(shí)，計(jì)算工作量很大，并且需要計(jì)算Rxx

的逆矩陣，從而要求的存貯量也很大。此外，在具體實(shí)現(xiàn)時(shí)，濾波器的長(zhǎng)度是由實(shí)驗(yàn)來確定的，如果想通過增加長(zhǎng)度提高逼近的精度，就需要在新M基礎(chǔ)上重新進(jìn)行計(jì)算。因此，從時(shí)域求解維納濾波器，并不是一個(gè)有效的方法。假定所研究的信號(hào)都是零均值的，濾波器為FIR型，長(zhǎng)度等于M，可以得到二、估計(jì)誤差的均方值若h(n)是維納霍夫方程解上式可以進(jìn)一步化簡(jiǎn)得到可以看出，均方誤差與濾波器的單位脈沖響應(yīng)是一個(gè)二次函數(shù)關(guān)系。由于單位脈沖響應(yīng)h(n)為M維向量，因此均方誤差是一個(gè)超橢圓拋物形曲面，該曲面有極小點(diǎn)存在。當(dāng)濾波器工作于最佳狀態(tài)時(shí)，均方誤差取得最小值。

例

設(shè)y(n)=x(n)+v2(n)，v2(n)是一白噪聲，方差σ22=0.1。期望信號(hào)x1(n)的信號(hào)模型如圖(a)所示，其中白噪聲v1(n)的方差σ21=0.27，且b0=0.8458。x(n)的信號(hào)模型如圖（b）所示，b1=0.9458。假定v1(n)與v2(n)、x1(n)與y(n)不相關(guān)，并都是實(shí)信號(hào)。設(shè)計(jì)一個(gè)維納濾波器，得到該信號(hào)的最佳估計(jì)，要求濾波器是一長(zhǎng)度為2的FIR濾波器。圖輸入信號(hào)與觀測(cè)數(shù)據(jù)的模型

解這個(gè)問題屬于直接應(yīng)用維納-霍夫方程的典型問題，其關(guān)鍵在于求出觀測(cè)信號(hào)的自相關(guān)函數(shù)和觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)函數(shù)。圖維納濾波器的框圖

根據(jù)題意，畫出這個(gè)維納濾波器的框圖，如圖所示。用H1(z)和H2(z)分別表示x1(n)和x(n)的信號(hào)模型，那么濾波器的輸入信號(hào)x(n)可以看作是v1(n)通過H1(z)和H２(z)級(jí)聯(lián)后的輸出，H1(z)和H２(z)級(jí)聯(lián)后的等效系統(tǒng)用H(z)表示，輸出信號(hào)y(n)就等于x(n)和v2(n)之和。因此求出輸出信號(hào)的自相關(guān)函數(shù)矩陣Ryy和輸出信號(hào)與期望信號(hào)的互相關(guān)矩陣Ryd是解決問題的關(guān)鍵。相關(guān)函數(shù)矩陣由相關(guān)函數(shù)值組成，已知x(n)與v2(n)不相關(guān)，那么

(1)求出期望信號(hào)的方差。根據(jù)圖

(a)，期望信號(hào)的時(shí)間序列模型所對(duì)應(yīng)的差分方程為x1(n)=v1(n)-b0x1(n-1)這里，b0=0.8458,由于x1(n)的均值為零，其方差與自相關(guān)函數(shù)在零點(diǎn)的值相等。

(2)計(jì)算輸入信號(hào)和輸出信號(hào)的自相關(guān)函數(shù)矩陣。根據(jù)自相關(guān)函數(shù)、功率譜密度和時(shí)間序列信號(hào)模型的等價(jià)關(guān)系，已知時(shí)間序列信號(hào)模型，就可以求出自相關(guān)函數(shù)。這里，信號(hào)的模型H(z)可以通過計(jì)算得到。這是一個(gè)二階系統(tǒng)，所對(duì)應(yīng)的差分方程為x(n)+a1x(n-1)+a2x(n-2)=v1(n)式中，a1=-0.1，a2=-0.8。由于v1(n)、v2(n)的均值為零，因此，x(n)的均值為0。將方程兩邊同乘以x*(n-m)，并取數(shù)學(xué)期望，得rxx(m)+a1rxx(m-1)+a2rxx(m-2)=0m＞0

(1)rxx(0)+a1rxx(1)+a2rxx(2)=σ21

m=0

(2)對(duì)方程（1）取m=1,2，得到rxx(1)+a1rxx(0)+a2rxx(1)=0 (3)rxx(2)+a1rxx(1)+a2rxx(0)=0 (4)方程(2)、(3)、(4)聯(lián)立求解，得至此，輸入信號(hào)的自相關(guān)矩陣Rxx可以寫出：

v2(n)是一個(gè)零均值的白噪聲，它的自相關(guān)函數(shù)矩陣呈對(duì)角形，且 ，因此，輸出信號(hào)的自相關(guān)Ryy為

(3)計(jì)算輸出信號(hào)與期望信號(hào)的互相關(guān)函數(shù)矩陣。由于兩個(gè)信號(hào)都是實(shí)信號(hào)，故ryd(m)=E［y(n)d(n-m)］=E［y(n)x1(n-m)］

=E［(x(n)+v2(n))x1(n-m)］

=E［x(n)x1(n-m)］

m=0,1根據(jù)系統(tǒng)H2(z)的輸入與輸出的關(guān)系，有x1(n)-b1x(n-1)=x(n)推出x1(n)=x(n)+b1x(n-1)這樣ryd(m)=E［x(n)x1(n-m)］

=E［x(n)(x(n-m)+b1x(n-1-m))］

=rxx(m)+b1rxx(m-1)將m=0,m=1代入上式,得ryd(0)=rxx(0)+b1rxx(-1)=1-0.9458×0.5=0.5272ryd(1)=rxx(1)+b1rxx(0)=0.5-0.9458×1=-0.4458因此，輸出信號(hào)與期望信號(hào)的互相關(guān)Ryd為求出輸出信號(hào)自相關(guān)的逆矩陣,并乘以Ryd,就可以得到維納濾波器的最佳解Wopt:可以計(jì)算出該維納濾波達(dá)到最佳狀態(tài)最小值E［|e(n)|2］min：

若不考慮濾波器的因果性，維納濾波式可以寫為設(shè)d(n)=s(n)，對(duì)上式兩邊做Z變換，得到Sxs(z)=Hopt(z)Sxx(z)2.3離散維納濾波器的z域解假設(shè)信號(hào)和噪聲不相關(guān)，即rsv(m)=0，則Sxs(z)=Sss(z)Sxx(z)=Sss(z)+Svv(z)則有當(dāng)噪聲為0時(shí)，信號(hào)全部通過；當(dāng)信號(hào)為0時(shí)，噪聲全部被抑制掉，因此維納濾波確有濾除噪聲的能力。把信號(hào)的頻譜用Pss(ejω)表示，噪聲的頻譜用Pvv(ejω)表示，那么非因果的維納濾波器的傳輸函數(shù)Hopt(ejω)的幅頻特性如圖所示。Pss(ejω)≠0,Pvv(ejω)=0Pss(ejω)≠0,Pvv(ejω)≠0Pss(ejω)=0,Pvv(ejω)≠0然而實(shí)際的系統(tǒng)都是因果的。對(duì)于一個(gè)因果系統(tǒng)，不能直接轉(zhuǎn)入頻域求解的原因是由于輸入信號(hào)與期望信號(hào)的互相關(guān)序列是一個(gè)因果序列，如果能夠把因果維納濾波器的求解問題轉(zhuǎn)化為非因果問題，求解方法將大大簡(jiǎn)化。那么怎樣把一個(gè)因果序列轉(zhuǎn)化為一個(gè)非因果序列呢？圖非因果維納濾波器的傳輸函數(shù)的幅頻特性

回顧前面講到的時(shí)間序列信號(hào)模型，假設(shè)x(n)的信號(hào)模型B(z)已知（如圖(a)所示），求出信號(hào)模型的逆系統(tǒng)B-1(z)，并將x(n)作為輸入，那么逆系統(tǒng)B-1(z)的輸出ω(n)為白噪聲。一般把信號(hào)轉(zhuǎn)化為白噪聲的過程稱為白化，對(duì)應(yīng)的濾波器稱為白化濾波器（如圖(b)所示）。圖x(n)的時(shí)間序列信號(hào)模型及其白化濾波器具體思路如圖所示。用白噪聲作為待求的維納濾波器的輸入，設(shè)定1/B(z)為信號(hào)x(n)的白化濾波器的傳輸函數(shù)，那么維納濾波器的傳輸函數(shù)G(z)的關(guān)系為因此，維納濾波器的傳輸函數(shù)H(z)的求解轉(zhuǎn)化為G(z)

的求解。圖維納濾波解題思路

假設(shè)待求維納濾波器的單位脈沖響應(yīng)為ω(n)，期望信號(hào)d(n)=s(n)，系統(tǒng)的輸出信號(hào)y(n)=s(n)，g(n)是G(z)的逆Z變換，如圖所示。一非因果維納濾波器的求解可以看出，均方誤差的第一項(xiàng)和第三項(xiàng)都是非負(fù)數(shù),要使均方誤差為最小，當(dāng)且僅當(dāng)

-∞＜k＜∞因此g(n)的最佳值為

-∞＜k＜∞對(duì)上式兩邊同時(shí)做Z變換，得到這樣，非因果維納濾波器的最佳解為因?yàn)閟(n)=s(n)*δ(n)，且x(n)=ω(n)*b(n)，根據(jù)相關(guān)卷積定理，得到rxs(m)=rωs(m)*b(-m)對(duì)上式兩邊做Z變換，得到Sxs(z)=Sωs(z)B(z-1)因此根據(jù)x(n)的信號(hào)模型，得到非因果的維納濾波器的復(fù)頻域最佳解的一般表達(dá)式假定信號(hào)與噪聲不相關(guān)，即當(dāng)E［s(n)v(n)］=0時(shí)，有rxs(m)=E［(s(n)+v(n))*s(n+m)］=rss(m)rxx(m)=E［(s(n)+v(n))*(s(n+m)+v(n+m))］

=rss(m)+rvv(m)對(duì)上邊兩式做Z變換，得到Sxs(z)=Sss(z)Sxx(z)=Sss(z)+Svv(z)信號(hào)和噪聲不相關(guān)時(shí)，非因果維納濾波器的復(fù)頻域最佳解和頻率響應(yīng)分別為濾波器的最小均方誤差E［|e(n)|2］min的計(jì)算，根據(jù)圍線積分法求逆Z變換的公式，rss(m)用下式表示:得出

由復(fù)卷積定理

取y(n)=x(n)，有

因此

因?yàn)閷?shí)信號(hào)的自相關(guān)函數(shù)是偶函數(shù)，即rss(m)=rss(-m)，因此

Sss(z)=Sss(z-1)假定信號(hào)與噪聲不相關(guān)，E［s(n)v(n)］=0，

則

若維納濾波器是一個(gè)因果濾波器，

要求

g(n)=0n＜0則濾波器的輸出信號(hào)

估計(jì)誤差的均方值

E［|e(n)|2］=E［|s(n)-y(n)|2］

得到

二因果維納濾波器的求解要使均方誤差取得最小值，

當(dāng)且僅當(dāng)

令

所以因果維納濾波器的復(fù)頻域最佳解為

維納濾波的最小均方誤差為

比較可以看出因果維納濾波器的最小均方誤差與非因果維納濾波器的最小均方誤差的形式相同，但公式中的Hopt(z)的表達(dá)式不同。前面已經(jīng)導(dǎo)出，對(duì)于非因果情況，對(duì)于因果情況，比較兩式，它們的第二項(xiàng)求和域不同，因?yàn)橐蚬闆r下,k=0~+∞，因此可以說明非因果情況的E［|e(n)|2］min一定小于等于因果情況E［|e(n)|2］min。在具體計(jì)算時(shí),可以選擇單位圓作為積分曲線，應(yīng)用留數(shù)定理，計(jì)算積分函數(shù)在單位圓內(nèi)的極點(diǎn)的留數(shù)來得到。因果維納濾波器設(shè)計(jì)的一般步驟：

(1)根據(jù)觀測(cè)信號(hào)x(n)的功率譜求出它所對(duì)應(yīng)的信號(hào)模型的傳輸函數(shù)，即采用譜分解的方法得到B(z)。具體方法為Sxx(z)=σ2ωB(z)B(z-1)，把單位圓內(nèi)的零極點(diǎn)分配給B(z)，單位圓外的零極點(diǎn)分配給B(z-1)，系數(shù)分配給σ2ω。

(2)求 的Z反變換，取其因果部分再做Z變換，即舍掉單位圓外的極點(diǎn)，得

(3)計(jì)算Hopt(z)，E［|e(n)|2］min。例

已知信號(hào)和噪聲不相關(guān)，即rsv(m)=0，噪聲v(n)是零均值、單位功率的白噪聲（σ2v=1，mv=0)，求Hopt(z)和E［|e(n)|2］min。

解根據(jù)白噪聲的特點(diǎn)得出Svv(z)=1，由噪聲和信號(hào)不相關(guān)，得到rxx(m)=rss(m)+rvv(m)。對(duì)上式兩邊做Z變換，并代入已知條件，對(duì)x(n)進(jìn)行功率譜分解:考慮到B(z)必須是因果穩(wěn)定的系統(tǒng)，得到連

(1)首先分析物理可實(shí)現(xiàn)情況，應(yīng)用公式（2.3.38）:令

F(z)的極點(diǎn)為0.8和2，考慮到因果性、穩(wěn)定性，僅取單位圓內(nèi)的極點(diǎn)zi=0.8，f(n)為F(z)的Z反變換。用Res表示留數(shù)，應(yīng)用留數(shù)定理，有取因果部分，

f+(n)=0.6×0.8n×u(n)令

單位圓內(nèi)只有極點(diǎn)zi=0.5，未經(jīng)濾波器的均方誤差

(2)對(duì)于非物理可實(shí)現(xiàn)情況，有

令

單位圓內(nèi)有兩個(gè)極點(diǎn)0.8和0.5，應(yīng)用留數(shù)定理，有

比較兩種情況下的最小均方誤差,可以看出非物理可實(shí)現(xiàn)情況的最小均方誤差小于物理可實(shí)現(xiàn)情況的均方誤差。

3.4.1維納預(yù)測(cè)的計(jì)算在維納濾波中，期望的輸出信號(hào)yd(n)=s(n)，實(shí)際的輸出為y(n)=s(n)。在維納預(yù)測(cè)中，期望的輸出信號(hào)yd(n)=s(n+N)，實(shí)際的輸出y(n)=s(n+N)。前面已經(jīng)推導(dǎo)得到維納濾波的最佳解為^^其中,Sxx(z)是觀測(cè)數(shù)據(jù)的功率譜;Sxyd(z)是觀測(cè)數(shù)據(jù)與期望信號(hào)的互功率譜，即互相關(guān)函數(shù)rxyd(k)的傅里葉變換4.4維納預(yù)測(cè)對(duì)應(yīng)于維納預(yù)測(cè)器，其輸出信號(hào)y(n)和預(yù)測(cè)誤差信號(hào)e(n+N)分別為

同理，要使預(yù)測(cè)誤差的均方值為最小，須滿足

其中，hk表示h(k)。

觀測(cè)數(shù)據(jù)與期望的輸出的互相關(guān)函數(shù)rxyd(k)和互譜密度Sxyd(z)分別為這樣，非因果維納預(yù)測(cè)器的最佳解為

因果維納預(yù)測(cè)器的最佳解為

維納預(yù)測(cè)的最小均方誤差為

從上面分析可以看出，

維納預(yù)測(cè)的求解和維納濾波器的求解方法是一致的。

假設(shè)x(n)=s(n)+v(n)，式中v(n)是噪聲，且v(n)=0，期望信號(hào)為s(n+N)，N＞0，此種情況稱為純預(yù)測(cè)。假定維納預(yù)測(cè)器是因果的，仍設(shè)s(n)與v(n)不相關(guān)，純預(yù)測(cè)情況下的輸入信號(hào)的功率譜及維納預(yù)測(cè)器的最佳解分別為二、純預(yù)測(cè)純預(yù)測(cè)器的最小均方誤差為

應(yīng)用復(fù)卷積定理

取y(n)=x(n)

并考慮到b(n)是因果系統(tǒng)，得到

可以看到，隨著N增加，E［|e(n+N)|2］min也增加。這一點(diǎn)也容易理解，當(dāng)預(yù)測(cè)的距離越遠(yuǎn)，預(yù)測(cè)的效果越差，偏差越大，因而E［|e(n+N)|2］min越大。例

已知其中-1＜a＜1，求:(1)最小均方誤差下的s(n+N)；

(2)E［|e(n+N)|2］min。^解首先對(duì)Sxx(z)進(jìn)行功率譜分解。因?yàn)樗云浯?，求出B(z)的Z反變換

然后，應(yīng)用Z變換的性質(zhì)，得到

圖4.4.1純預(yù)測(cè)維納濾波器由Hopt(z)=aN，此時(shí)可以把純預(yù)測(cè)的維納濾波器看作是一個(gè)線性比例放大器（如圖4.4.1所示）。根據(jù)x(n)的信號(hào)模型可以寫出x(n)的時(shí)間序列模型所對(duì)應(yīng)的輸入輸出方程x(n)=ω(n)+ax(n-1)將信號(hào)x(n)通過純預(yù)測(cè)維納濾波器，隨著時(shí)間的遞增，可以得到當(dāng)N=1時(shí)，x(n+1)=ax(n)=as(n)當(dāng)N=2時(shí)，x(n+2)=ax(n+1)=a2s(n)當(dāng)N=N時(shí)，x(n+N)=ax(N+n-1)=aNs(n)…以上推導(dǎo)結(jié)果相當(dāng)于在n+N時(shí)刻，ω(n+N)=0，即去掉噪聲時(shí)的結(jié)果。設(shè)N＞0時(shí)，ω(n+N)=0，則x(n+N)=ax(n+N-1)此時(shí)，從統(tǒng)計(jì)意義上講，當(dāng)N＞0時(shí)，白噪聲信號(hào)ω(n+N)對(duì)x(n)無影響。這一結(jié)論還可以推廣，對(duì)于任何均值為零的x(n)，要估計(jì)s(n+N)時(shí)，只需要考慮B(z)的慣性，即可認(rèn)為ω(n+N)=0，N＞0，這樣估計(jì)出來的結(jié)果將有最小均方誤差。^表明一個(gè)信號(hào)的功率譜在單位圓上沒有極點(diǎn)與信號(hào)均值等于0等價(jià)，因此對(duì)于功率譜在單位圓上沒有極點(diǎn)的信號(hào)，要估計(jì)s(n+N)時(shí)，可認(rèn)為ω(n+N)=0,N＞0，即僅需要考慮B(z)的慣性，這樣估計(jì)出來的結(jié)果將有最小均方誤差。^終值定理

已知x(n-1),x(n-2)，…，x(n-p)，預(yù)測(cè)x(n)，假設(shè)噪聲v(n)=0，這樣的預(yù)測(cè)稱為一步線性預(yù)測(cè)。設(shè)定系統(tǒng)的單位脈沖響應(yīng)為h(n)，根據(jù)線性系統(tǒng)的基本理論，輸出信號(hào)令apk=-h(k)，則

預(yù)測(cè)誤差

三、一步線性預(yù)測(cè)的時(shí)域解其中，ap0=1,要使均方誤差為最小值，要求

同維納濾波的推導(dǎo)過程一樣,可以得到

E［e*(n)x(n-l)］=0l=1,2,…,p

l=1,2,…,p

由于預(yù)測(cè)器的輸出是輸入信號(hào)的線性組合，

得到

說明誤差信號(hào)與輸入信號(hào)滿足正交性原理預(yù)測(cè)誤差與預(yù)測(cè)的信號(hào)值同樣滿足正交性原理。預(yù)測(cè)誤差的最小均方值得到下面的方程組：

將方程組寫成矩陣形式

這就是有名的Yule-Walker方程，可以看出Yule-Walker方程具有以下特點(diǎn)：

(1)除了第一個(gè)方程外，其余都是齊次方程;(2)與維納-霍夫方程相比，不需要知道觀測(cè)數(shù)據(jù)x(n)與期望信號(hào)s(n)的互相關(guān)函數(shù)。該方程組有p+1個(gè)方程，對(duì)應(yīng)地，可以確定apk，k=1,2,…，p和E［e2(n)］min，共計(jì)p+1個(gè)未知數(shù)，因此可用來求解AR模型參數(shù)。這就是后面要介紹的AR模型法進(jìn)行功率譜估計(jì)的原理，它再一次揭示了時(shí)間序列信號(hào)模型、功率譜和自相關(guān)函數(shù)描述一個(gè)隨機(jī)信號(hào)的等價(jià)性。例

已知x(n)為AR模型，求AR模型參數(shù)。

解求解AR模型參數(shù)包括確定AR模型的階數(shù)p及系數(shù)ap1，ap2，…，app。首先對(duì)Sxx(z)做傅里葉反變換，得到x(n)的自相關(guān)函數(shù)rxx(m),rxx(m)=0.8|m|

采用試驗(yàn)的方法確定模型階數(shù)p。首先取p=2，各相關(guān)函數(shù)值由上式計(jì)算，并代入(3.4.29)式，可得計(jì)算得到a1=-0.8，a2=0,σ2ω=0.36

如果取p=3，可計(jì)算出a1=-0.8,a2=a3=0,σ2ω=0.36，說明AR模型的階數(shù)只能是一階的。采用譜分解的方法，即對(duì)Sxx(z)進(jìn)行譜分解，得到的模型也是一階的，其時(shí)間序列模型和差分方程為卡爾曼濾波是用狀態(tài)空間法描述系統(tǒng)的，由狀態(tài)方程和量測(cè)方程所組成?？柭鼮V波用前一個(gè)狀態(tài)的估計(jì)值和最近一個(gè)觀測(cè)數(shù)據(jù)來估計(jì)狀態(tài)變量的當(dāng)前值，并以狀態(tài)變量的估計(jì)值的形式給出?？柭鼮V波具有以下的特點(diǎn)：

(1)算法是遞推的，且狀態(tài)空間法采用在時(shí)域內(nèi)設(shè)計(jì)濾波器的方法，因而適用于多維隨機(jī)過程的估計(jì)；離散型卡爾曼算法適用于計(jì)算機(jī)處理。

(2)用遞推法計(jì)算，不需要知道全部過去的值，用狀態(tài)方程描述狀態(tài)變量的動(dòng)態(tài)變化規(guī)律，因此信號(hào)可以是平穩(wěn)的，也可以是非平穩(wěn)的，即卡爾曼濾波適用于非平穩(wěn)過程。

(3)卡爾曼濾波采取的誤差準(zhǔn)則仍為均方誤差最小準(zhǔn)則。4.5卡爾曼(Kalman)濾波假設(shè)某系統(tǒng)k時(shí)刻的狀態(tài)變量為xk，狀態(tài)方程和量測(cè)方程（也稱為輸出方程）表示為(4.5.1a)(4.5.1b)其中，k表示時(shí)間，這里指第k步迭代時(shí)，相應(yīng)信號(hào)的取值；輸入信號(hào)ωk是一白噪聲，輸出信號(hào)的觀測(cè)噪聲vk

也是一個(gè)白噪聲，輸入信號(hào)到狀態(tài)變量的支路增益等于1，即B=1；一、卡爾曼濾波的狀態(tài)方程和量測(cè)方程

A表示狀態(tài)變量之間的增益矩陣，可以隨時(shí)間發(fā)生變化，用Ak表示第k步迭代時(shí)，增益矩陣A的取值；C表示狀態(tài)變量與輸出信號(hào)之間的增益矩陣，可以隨時(shí)間變化，第k步迭代時(shí)，取值用Ck表示，其信號(hào)模型如圖4.5.1所示。將狀態(tài)方程中時(shí)間變量k用k-1代替，得到的狀態(tài)方程和量測(cè)方程如下所示：xk=Ak-1xk-1+ωk-1

yk=Ckxk+vk

其中，xk是狀態(tài)變量;ωk-1表示輸入信號(hào)是白噪聲;vk是觀測(cè)噪聲;yk是觀測(cè)數(shù)據(jù)。圖

4.5.1卡爾曼濾波器的信號(hào)模型

其中

為了后面的推導(dǎo)簡(jiǎn)單起見，假設(shè)狀態(tài)變量的增益矩陣A不隨時(shí)間發(fā)生變化，ωk，vk都是均值為零的正態(tài)白噪聲，方差分別是Qk和Rk，并且初始狀態(tài)與ωk，vk都不相關(guān)，表示相關(guān)系數(shù)。即

卡爾曼濾波是采用遞推的算法實(shí)現(xiàn)的，其基本思想是先不考慮輸入信號(hào)ωk和觀測(cè)噪聲vk的影響，得到狀態(tài)變量和輸出信號(hào)（即觀測(cè)數(shù)據(jù)）的估計(jì)值，再用輸出信號(hào)的估計(jì)誤差加權(quán)后校正狀態(tài)變量的估計(jì)值，使?fàn)顟B(tài)變量估計(jì)誤差的均方值最小。因此,卡爾曼濾波的關(guān)鍵是計(jì)算出加權(quán)矩陣的最佳值。當(dāng)不考慮觀測(cè)噪聲和輸入信號(hào)時(shí)，狀態(tài)方程和量測(cè)方程為(4.5.4)(4.5.5)二卡爾曼濾波的遞推算法顯然，由于不考慮觀測(cè)噪聲的影響，輸出信號(hào)的估計(jì)值與實(shí)際值是有誤差的，用表示

(4.5.6)為了提高狀態(tài)估計(jì)的質(zhì)量，用輸出信號(hào)的估計(jì)誤差來校正狀態(tài)變量

(4.5.7)其中，Hk為增益矩陣，實(shí)質(zhì)是一加權(quán)矩陣。經(jīng)過校正后的狀態(tài)變量的估計(jì)誤差及其均方值分別用和Pk表示，把未經(jīng)校正的狀態(tài)變量的估計(jì)誤差的均方值用表示(4.5.8)(4.5.9)(4.5.10)

卡爾曼濾波要求狀態(tài)變量的估計(jì)誤差的均方值Pk為最小，因此卡爾曼濾波的關(guān)鍵就是要得到Pk與Hk的關(guān)系式，即通過選擇合適的Hk,使Pk取得最小值。首先推導(dǎo)狀態(tài)變量的估計(jì)值和狀態(tài)變量的估計(jì)誤差，然后計(jì)算的均方值Pk

，并通過化簡(jiǎn)Pk，得到一組卡爾曼濾波的遞推公式。將（4.5.3）、

(4.5.5)式代入（4.5.7）式

(4.5.11)同理，狀態(tài)變量的估計(jì)誤差為

(4.5.12)由上式可以看出，狀態(tài)變量的估計(jì)誤差由三部分組成，可記為其中(4.5.13b)(4.5.13c)(4.5.13d)那么，狀態(tài)變量的估計(jì)誤差的均方值Pk就由9項(xiàng)組成:(4.5.14a)其中

(4.5.14b)(4.5.14d)(4.5.14c)

下面化簡(jiǎn)Pk的表達(dá)式，根據(jù)假設(shè)的條件，狀態(tài)變量的增益矩陣A不隨時(shí)間發(fā)生變化，起始時(shí)刻為k0，則（4.5.2）式經(jīng)過迭代，得到令l=k-k0-j，得到取k0=0，k=k-1，得到(4.5.15)所以xk-1僅依賴于x0,ω0,ω1,…,ωk-2,與ωk-1不相關(guān)，即(4.5.16)又據(jù)(4.5.7)式和(4.5.3)式，得(4.5.17)所以僅依賴于xk-1,vk-1，而與vk不相關(guān)，即(4.5.18)(4.5.19)

把（4.5.15）～（4.5.19）式代入(4.5.14)式,Pk中的9項(xiàng)可以分別化簡(jiǎn)為(4.5.20a)(4.5.20b)(4.5.20c)(4.5.20d)(4.5.20e)(4.5.20f)(4.5.20g)(4.5.20h)(4.5.20j)也就是說，Pk僅有其中的三項(xiàng)不為零，化簡(jiǎn)成(4.5.21)為了進(jìn)一步化簡(jiǎn)Pk，推導(dǎo)未經(jīng)誤差校正的狀態(tài)估計(jì)誤差的均方值Pk′，由下面推導(dǎo)結(jié)果可以看出，Pk′是一對(duì)稱矩陣，滿足Pk′=(Pk′)T。(4.5.22)將(4.5.22)式代入(4.5.21)式，即把Pk′代入Pk,(4.5.23)其中, 是正定陣，記

(4.5.24)令

(4.5.25)將上式代入（4.5.23）式，得(4.5.26)將(4.5.26)式后三項(xiàng)配對(duì)(4.5.27)第二項(xiàng)和第三項(xiàng)均與Hk無關(guān)，第一項(xiàng)為一半正定陣，因此使Pk最小的Hk應(yīng)滿足(4.5.28)(4.5.29)將Hopt代入Pk，得到最小均方誤差陣將(4.5.7)、(4.5.22)、(4.5.29)式和(4.5.30)式聯(lián)立，得到一組卡爾曼遞推公式(4.5.30)(4.5.31a)(4.5.31b)(4.5.31c)(4.5.31d)假設(shè)初始條件Ak,Ck,Qk,Rk,yk,xk-1,Pk-1已知，其中x0=E［x0］,P0=var［x0］,那么，遞推流程見圖4.5.2。^^圖

4.5.2卡爾曼濾波遞推流程

例

已知信號(hào)與噪聲不相關(guān)，yk=xk+vk,求卡爾曼信號(hào)模型中的Ak和Ck。

解由yk=xk+vk知道，Ck=1。

對(duì)Sxx(z)進(jìn)行譜分解，確定x(n)的信號(hào)模型B(z)，從而確定Ak。根據(jù)Sxx(z)=σ2ωB(z)B(z-1)，得出

上式與卡爾曼狀態(tài)方程相比，不同之處在于輸入信號(hào)ω(n)的時(shí)間不同，因此將Sxx(z)改寫為再對(duì)Sxx(z)進(jìn)行譜分解，得到

(解畢)

卡爾曼濾波和維納濾波都是采用均方誤差最小的準(zhǔn)則來實(shí)現(xiàn)信號(hào)濾波的，但維納濾波是在信號(hào)進(jìn)入了穩(wěn)態(tài)后的分析，卡爾曼濾波是從初始狀態(tài)采用遞推的方法進(jìn)行濾波。對(duì)于平穩(wěn)隨機(jī)信號(hào)，當(dāng)過渡過程結(jié)束以后，卡爾曼濾波與維納濾波的結(jié)果間存在什么關(guān)系呢?下面舉一例說明。

例

已知在k=0時(shí)開始觀察yk,yk=xk+vk，用卡爾曼過濾的計(jì)算公式求xk，并與維納過濾的方法進(jìn)行比較。

解(1)由x(n)功率譜及量測(cè)方程，確定卡爾曼遞推算法。首先對(duì)Sxx(z)進(jìn)行功率譜分解，由例4.5.1的結(jié)果，得到卡爾曼濾波的狀態(tài)方程為xk=0.8xk-1+ωk-1,確定Ak=0.8由量測(cè)方程yk=xk+vk,確定Ck=1，將參數(shù)矩陣Ak,Ck,Rk代入卡爾曼遞推公式（4.5.30），得到(4.5.32a)(4.5.32b)(4.5.32c)(4.5.32d)

(2)求出卡爾曼濾波的輸出。由卡爾曼遞推公式，以及 ，P0=var［x0］=1,可得到Pk′,Hk,Pk及xk（k表示迭代次數(shù)），迭代流程為： 由具體迭代結(jié)果可以看出，原先的增益矩陣Ak,由于只選擇了一個(gè)狀態(tài)變量，變成了加權(quán)系數(shù)。見表4.5.1。表4.5.1Kalman濾波迭代結(jié)果

(3)求出卡爾曼濾波的穩(wěn)態(tài)解。將(4.5.32b)式代入方程(4.5.32d),得到第5個(gè)方程(4.5.32e)將方程(4.5.32c)、(4.5.32e)代入方程(4.5.32d),消去Pk′，可以得到Pk的遞推關(guān)系:Pk=(1-Pk)［0.64Pk-1+0.36］

=0.64Pk-1-0.64Pk

-1

Pk+0.36-0.36Pk

化簡(jiǎn)上式，得到

1.36Pk+0.64Pk-1Pk=0.64Pk-1+0.36要求的是穩(wěn)態(tài)解，因此將Pk,Pk-1都用P∞代替,得到

根據(jù)P∞，可以確定達(dá)到穩(wěn)態(tài)后的卡爾曼濾波的狀態(tài)方程:(4.5.33)

(4)用維納濾波的方法分析。采用功率譜分解的方法，得到x(n)的時(shí)間序列信號(hào)模型的傳輸函數(shù)H(z):上式說明x是一階AR模型，對(duì)H(z)做Z反變換得到(4.5.34)比較（4.5.33）式和（4.5.34）式，可以看出卡爾曼濾波的穩(wěn)態(tài)解與維納解是相等的。(解畢)通過上面的例題，可以看出維納濾波是已知前p個(gè)觀測(cè)數(shù)據(jù)及信號(hào)與噪聲的相關(guān)函數(shù)，通過建立模型的方法分析的?？柭鼮V波要求已知前一個(gè)時(shí)刻的狀態(tài)估計(jì)值x(k-1)和當(dāng)前的觀測(cè)值yk，由狀態(tài)方程和量測(cè)方程遞推得到結(jié)果。維納濾波的解以H(z)的形式給出，卡爾曼濾波是以狀態(tài)變量的估計(jì)值給出解的形式。它們都采用均方誤差最小的準(zhǔn)則,但卡爾曼濾波有一個(gè)過渡過程，其結(jié)果與維納濾波不完全相同，但到達(dá)穩(wěn)態(tài)后，結(jié)果相同。^下面舉一個(gè)雷達(dá)跟蹤目標(biāo)物的例子說明卡爾曼濾波的應(yīng)用。雷達(dá)跟蹤目標(biāo)的基本原理是通過發(fā)射脈沖，根據(jù)接收到的脈沖與發(fā)射脈沖的時(shí)間間隔，來確定目標(biāo)物的距離和速度。由于干擾的影響，接收到的脈沖波形變化很大，那么一次的測(cè)量結(jié)果可能存在很大的誤差。為了減小誤差，往往采取發(fā)射一串脈沖的方法進(jìn)行測(cè)量。

三、應(yīng)用舉例我們知道，空間中的一點(diǎn)需要由徑向距離和方位角來確定。假設(shè)雷達(dá)跟蹤的目標(biāo)為飛行器，發(fā)射的脈沖時(shí)間間隔為T，在時(shí)間k，徑向距離為R+ρ(k)，在時(shí)間k+1，距離為R+ρ(k+1)，兩者之間有秒的延時(shí)，因此T表示空間一次掃描的時(shí)間間隔。平均距離用R表示，ρ(k)和ρ(k+1)表示對(duì)平均值的偏差。假定偏差是統(tǒng)計(jì)隨機(jī)的，均值為零,那么可以寫出距離方程式中，表示速度。用u表示加速度，則可以得到加速度方程假定加速度u(k)是零均值的平穩(wěn)白噪聲，即滿足為了后面敘述方便，定義x(k)表示第k個(gè)雷達(dá)回波脈沖獲得的目標(biāo)距離，z(k)表示在第k個(gè)雷達(dá)脈沖進(jìn)行數(shù)據(jù)處理之后的目標(biāo)距離估計(jì)，z(k)表示在第k個(gè)雷達(dá)脈沖進(jìn)行數(shù)據(jù)處理之后的目標(biāo)速度估計(jì)。設(shè)定狀態(tài)變量為x(k)，選擇的狀態(tài)變量有4個(gè)，分別表示徑向距離、徑向速度、方位角和方位角速度，即

.（4.5.35）

根據(jù)狀態(tài)變量的物理含義，得到以下方程：

式中，u1