李 劍

(中國民用航空呼倫貝爾空中交通管理站 內(nèi)蒙古 海拉爾 021008)

卡爾曼濾波算法是一種遞歸計(jì)算方法，最初是通過最優(yōu)估計(jì)的方法，來獲得離散數(shù)據(jù)的線性濾波，隨著計(jì)算機(jī)技術(shù)的提升，其應(yīng)用滲透多個(gè)領(lǐng)域。在航空方面，卡爾曼濾波算法由于功能強(qiáng)大，不止能夠估計(jì)信號(hào)之前及當(dāng)前狀態(tài)，甚至可以據(jù)此估計(jì)信號(hào)的未來狀態(tài)。當(dāng)前民航系統(tǒng)使用的空管二次雷達(dá)及ADS-B等定位設(shè)備，ATC終端使用的自動(dòng)化設(shè)備等，也均使用卡爾曼濾波算法解決目標(biāo)的航跡定位問題。

1 卡爾曼濾波算法形式

1.1 標(biāo)準(zhǔn)形式

卡爾曼濾波算法的核心思想是提供一種算法，通過對(duì)目標(biāo)狀態(tài)的估計(jì)來使過程的均方誤差達(dá)到最小，在估計(jì)狀態(tài)時(shí)，使用的測(cè)量數(shù)據(jù)也受到干擾噪聲影響下進(jìn)行。[1]在做推導(dǎo)時(shí)，隨機(jī)信號(hào)的目標(biāo)狀態(tài)方程使用x(n)表示；而測(cè)量數(shù)據(jù)k(n)跟隨目標(biāo)也為線性狀態(tài)，其中觀測(cè)噪聲服從零均值分布，則卡爾曼濾波算法的模型方程為：

(1)

式中，w(n)和v(n)為目標(biāo)及測(cè)量白噪聲,a和c為常數(shù)，模值均小于1。

(2)

以上即卡爾曼濾波算法的標(biāo)準(zhǔn)形式，它把目標(biāo)的狀態(tài)描述為一個(gè)離散的隨機(jī)過程，該過程中系統(tǒng)狀態(tài)及其測(cè)量均不可避免受噪聲影響?？柭鼮V波算法實(shí)際是維納濾波的一種特殊表現(xiàn)形式。雷達(dá)跟蹤觀測(cè)所使用的濾波算法為一種矢量算法，以下介紹其推導(dǎo)過程。

1.2 卡爾曼濾波算法矢量方程

可以考慮兩種情況，一種是根據(jù)觀測(cè)數(shù)據(jù)同時(shí)推算估計(jì)若干個(gè)類似信號(hào)，例如q個(gè)信號(hào)x1,x2…,xq；另一種是得到一個(gè)高階自回歸過程,如一個(gè)q階的自回歸過程：

(3)

以上兩種情況下，如果引入矢量卡爾曼濾波算法，分析計(jì)算可以變得更加便捷。

以上情況中q個(gè)一階自回歸信號(hào)，我們?cè)O(shè)其為同時(shí)相互獨(dú)立的，則在n時(shí)刻，其取樣分別為x1,x2…，xq,它們的狀態(tài)方程表示為：

xi(n)=aixi(n-1)+wi(n),i=1,2,3…,q

(4)

同標(biāo)量卡爾曼相似，使用wi(n)表示零均值白噪聲序列，它們相互之間相關(guān)。q個(gè)信號(hào)xi(n)使用矩陣表示為q維矢量x(n)=[x1(n)x2(n)…xq(n)]T，則我們可以將以上信號(hào)狀態(tài)方程簡(jiǎn)化為以下矢量方程：

x(n)=Ax(n-1)+w(n)

(5)

分量w(n)是由wi(n)第q維構(gòu)成的矢量，A是一個(gè)q階對(duì)角矩陣，分量為ai。

同理在n時(shí)刻測(cè)得k個(gè)數(shù)據(jù)k1(n),k2(n),…,kk(n)，它們與xi(n)的關(guān)系為：

ki(n)=aixi(n-1)+vi(n),i=1,2,3…,k

(6)

其中,k≤q,測(cè)量噪聲使用vi(n)表示。根據(jù)定義數(shù)據(jù)矢量和噪聲矢量分別為以下等式：

k(n)=[k1(n)k2(n)…kk(n)]T和v(n)=

[v1(n)v2(n)…vk(n)]T

(7)

則K個(gè)測(cè)量方程的矢量方程簡(jiǎn)化為：

k(n)=Cx(n)+w(n)

(8)

系數(shù)矩陣C是一個(gè)K×q矩陣。

已知以上標(biāo)量卡爾曼濾波算法器遞推公式，同理可推導(dǎo)出矢量卡爾曼濾波算法器的相應(yīng)公式為：

(9)

最小預(yù)測(cè)均方誤差：p(n)=A(n)ξ(n-1)AT(n)+Q(n)

(10)

(11)

(12)

2 雷達(dá)及ADS-B跟蹤中的的卡爾曼濾波算法器

2.1 雷達(dá)跟蹤

二次雷達(dá)系統(tǒng)是民航對(duì)飛行目標(biāo)監(jiān)視使用最廣泛的設(shè)備，是通過地面詢問機(jī)對(duì)目標(biāo)進(jìn)行詢問，目標(biāo)在收到詢問信號(hào)并對(duì)其進(jìn)行應(yīng)答來實(shí)現(xiàn)對(duì)目標(biāo)定位的一種設(shè)備。一般對(duì)目標(biāo)的定位采用極坐標(biāo)的形式，即可用徑向距離及方位角表示目標(biāo)的實(shí)際位置。

在跟蹤系統(tǒng)中，一般利用詢問與應(yīng)答的時(shí)間延遲去計(jì)算目標(biāo)斜距，并近似作為飛機(jī)徑向距離，而目標(biāo)的方位角是在天線波束掃描到目標(biāo)時(shí)雷達(dá)系統(tǒng)碼盤刻度。二次雷達(dá)天線轉(zhuǎn)動(dòng)周期一般設(shè)置為4S，我們用表示其轉(zhuǎn)動(dòng)周期，同時(shí)這也是二次雷達(dá)位置更新的時(shí)間間隔。

(13)

(14)

這在濾波算法中便等效于零均值白噪聲即：

(15)

(16)

得出以下矢量矩陣公式：

x(n+1)=Ax(n)+w(n)

(17)

在雷達(dá)天線半波束輻射范圍內(nèi)，目標(biāo)的方位角由天線碼盤及單脈沖信息獲得，而其距離則借助光在空間勻速直線傳播得出(此時(shí)忽略目標(biāo)高度)。目標(biāo)離開雷達(dá)天線波束后，雷達(dá)輸出點(diǎn)跡信息，此過程中將測(cè)量結(jié)果k1(n)和k2(n)寫成矢量矩陣形式：

k(n)=Cx(n)+v(n)

(18)

矢量w(n)和v(n)各自的方差為自相關(guān)矩陣Q(n)和R(n)：

(19)

(20)

2.2 初始值確定及影響濾波算法的因素

在跟蹤算法中，卡爾曼濾波算法器需要得到增益矩陣Gn的初值。為此，在點(diǎn)跡或航跡形成之初，需確定均方誤差ξ(n)的初始值。一般點(diǎn)跡或航跡的初始化n=3時(shí)形成，及利用前兩次n=1和n=2測(cè)量得到的距離及方位角參數(shù)，得到四個(gè)數(shù)據(jù)k1(2),k2(1),k1(2),k2(2),并對(duì)四個(gè)數(shù)據(jù)k作以下估計(jì)：

(21)

代入以上測(cè)量方程得：

(22)

(23)

將狀態(tài)方程代入得：

(24)

綜合以上一系列公式，可得：

(25)

對(duì)均方誤差做出以下計(jì)算：(2.14)

(26)

即可確定均方誤差ξ(n)的初始值為：

(27)

由于海拉爾二次雷達(dá)建設(shè)受周邊建筑影響，雷達(dá)站周邊環(huán)境急劇惡化，近距離鏡像及假目標(biāo)出現(xiàn)頻繁，反射使得在對(duì)初值、測(cè)距及測(cè)角均方差值賦值時(shí)偏離實(shí)際較大，我們知道,鏡像的出現(xiàn)是存在較為固定的反射體,這時(shí)假目標(biāo)形成的濾波航跡較為穩(wěn)定,雷達(dá)本身很難識(shí)別,技術(shù)人員對(duì)功率可編程、STC曲線進(jìn)行了調(diào)節(jié)，降低了近距離物體反射的影響，情況得以好轉(zhuǎn)，但依舊不能解決零散假目標(biāo)形成的點(diǎn)航跡問題。

對(duì)于鏡像，可以使用抗反射的方法處理，但由于附近居民樓影響，零散雜波較為普遍，就需要?jiǎng)e的措施解決雷達(dá)濾波的方差賦值影響，SELEX二次雷達(dá)在航跡跟蹤門限設(shè)置中，提供了航跡速度初始值及跟蹤維護(hù)的門限設(shè)置，根據(jù)日常目標(biāo)在高空巡航及中低空進(jìn)近的常規(guī)速度，我們對(duì)以上參數(shù)做出門限設(shè)置，雜波假目標(biāo)的狀況得到明顯改善。

圖1 二次雷達(dá)航跡濾波初始值及跟蹤維護(hù)門限設(shè)置

2.3 卡爾曼濾波算法器在ADS-B中的應(yīng)用

廣播式自動(dòng)相關(guān)(ADS-B)的基本原理是飛機(jī)(航空器)通過空地、空空數(shù)據(jù)鏈，采用周期無方向性廣播方式，自動(dòng)發(fā)送和接收記載設(shè)備所提取的監(jiān)控信息。[2]ADS-B數(shù)據(jù)更新率及準(zhǔn)確率均優(yōu)于二次雷達(dá)。近年來，海拉爾處于支線機(jī)場(chǎng)航空運(yùn)輸快速增長(zhǎng),通航機(jī)場(chǎng)不斷增多的現(xiàn)狀，中低空空域機(jī)構(gòu)日漸復(fù)雜，需要大量的監(jiān)視源設(shè)備填補(bǔ)管制監(jiān)控盲區(qū)。目前，中國民航正在全力推行民航客機(jī)二次雷達(dá)應(yīng)答機(jī)加裝ADS-B應(yīng)答功能，而且ADS-B建站成本低廉，可很好提升海拉爾管制區(qū)監(jiān)視冗余備份能力。

與雷達(dá)跟蹤類似，根據(jù)上述推算，ADS-B的卡爾曼增益方程可以表示為[3]：

Gk=P(k)HT(k)[H(k)P(k)·HT(k)+R(k)]-1

(28)

航空目標(biāo)的狀態(tài)方程即為：

(29)

對(duì)應(yīng)的協(xié)方差更新公式為：

ξ(k)=[I-H(k)GK]p(k)

(30)

其中HT(k)為航空目標(biāo)的觀測(cè)矩陣。

不同于二次雷達(dá)，ADS-B應(yīng)用中需要跟蹤的目標(biāo)狀態(tài)變量表示為T=(x,vx,y,vy,w,z)；分別代表目標(biāo)在X和Y方向的距離、速度值、轉(zhuǎn)彎速率及高度，其中表示實(shí)時(shí)估計(jì)值；模型建立后，使用同一目標(biāo)不同周期的位置信息，來確定目標(biāo)目前狀態(tài)、預(yù)計(jì)狀態(tài)以及相應(yīng)的協(xié)方差矩陣。

2.4 空管自動(dòng)化中的濾波算法應(yīng)用

呼倫貝爾ATC系統(tǒng)目前使用民航二所的AirNet自動(dòng)化系統(tǒng)，接入數(shù)據(jù)將不同監(jiān)視源的同一目標(biāo)的數(shù)據(jù)進(jìn)行融合，系統(tǒng)將各監(jiān)視源的航跡融合成系統(tǒng)綜合航跡。并對(duì)每個(gè)周期的綜合航跡跟蹤濾波，使處理后的系統(tǒng)航跡更加接近真實(shí)航跡。其跟蹤使用IMM算法(InteractingMultipleModel，交互式多模型算法)。將目標(biāo)在水平面內(nèi)的運(yùn)動(dòng)分三個(gè)分量。

IMM算法的核心思想是同時(shí)維護(hù)三個(gè)濾波[4]，兩個(gè)標(biāo)準(zhǔn)CT模型和一個(gè)CV模型，分別對(duì)應(yīng)左轉(zhuǎn)彎、右轉(zhuǎn)彎及直線運(yùn)動(dòng)，運(yùn)算流程如圖2所示，F(xiàn)ilter0,Filter1,Filter2分別代表三種不同的運(yùn)動(dòng)模型，每個(gè)模型單獨(dú)濾波，初始值賦值后，所有狀態(tài)值的更新依賴于上次濾波運(yùn)算結(jié)果，不斷更新概率模型，同時(shí)將更新作為下一次估計(jì)的初始輸入，從而實(shí)現(xiàn)對(duì)目標(biāo)狀態(tài)的持續(xù)跟蹤。

圖2 IMM基本流程圖

監(jiān)視源權(quán)重取值與監(jiān)視源目標(biāo)信息可靠性相關(guān)。目前二次雷達(dá)仍然為區(qū)域管制的主要監(jiān)視手段，但ADS-B數(shù)據(jù)處理中心匯接了多個(gè)地面站信息，在中低空覆蓋方面與二次雷達(dá)形成顯著優(yōu)勢(shì)，且數(shù)據(jù)更新率快，精度質(zhì)量高，有利于目標(biāo)的平滑處理。根據(jù)實(shí)際情況，AirNet自動(dòng)化系統(tǒng)對(duì)雷達(dá)及ADS-B靜態(tài)權(quán)重進(jìn)行了以下分配：

圖3 AirNet自動(dòng)化系統(tǒng)監(jiān)視源權(quán)重分配情況

目前海拉爾管制區(qū)內(nèi)新建多臺(tái)ADS-B地面站，地面站數(shù)據(jù)經(jīng)數(shù)據(jù)處理中心匯總后送入自動(dòng)化系統(tǒng)，考慮數(shù)據(jù)完整性，自動(dòng)化系統(tǒng)接收所有NUCP取值范圍內(nèi)的數(shù)據(jù)進(jìn)行數(shù)據(jù)融合，初期出現(xiàn)了較多假目標(biāo)，經(jīng)分析及確認(rèn)，先后有GPS干擾、RAIM失效(NUCP取值低)、與雷達(dá)信號(hào)目標(biāo)分裂等情況，影響濾波算法準(zhǔn)確性。由于GNSS不確定因素會(huì)造成NUCP值對(duì)應(yīng)的位置精度及完好性與實(shí)際情況不符[5],經(jīng)與廠家溝通，增加了對(duì)ADS-B數(shù)據(jù)質(zhì)量驗(yàn)證的設(shè)置選項(xiàng)，改善了航跡濾波效果。

3 自適應(yīng)卡爾曼濾波算法器的應(yīng)用

實(shí)際中，濾波器的性能在一定條件下會(huì)受到Q(n)和R(n)陣的取值影響。無論單雷達(dá)航跡跟蹤、靜態(tài)融合或IMM算法，Q(n)取值變化與狀態(tài)濾波估計(jì)值依賴外推預(yù)測(cè)值的權(quán)重呈反比狀態(tài)；對(duì)于R(n)，濾波器在性能上主要跟隨觀測(cè)方程權(quán)重亦呈反比狀態(tài)。在大多數(shù)目標(biāo)跟蹤中，目標(biāo)狀態(tài)呈現(xiàn)多種變化趨勢(shì)，觀測(cè)值也受設(shè)備性能及現(xiàn)場(chǎng)環(huán)境影響產(chǎn)生較多變數(shù)，Q(n)和R(n)的取值應(yīng)該呈現(xiàn)動(dòng)態(tài)變化的趨勢(shì)。

實(shí)際上，像自動(dòng)化系統(tǒng)中對(duì)多監(jiān)視源的數(shù)據(jù)進(jìn)行靜態(tài)及動(dòng)態(tài)權(quán)重的分配，也可看作自適應(yīng)卡爾曼濾波算法的一種表現(xiàn)形式，即便是擴(kuò)展卡爾曼濾波算法的應(yīng)用，也是在將非線性目標(biāo)通過泰勒級(jí)數(shù)展開、樣本點(diǎn)獲取權(quán)重、平方根取樣等手段，達(dá)到目標(biāo)跟蹤濾波算法線性化的效果。另外一種優(yōu)化方法，是同時(shí)建立多個(gè)卡爾曼濾波算法模型，對(duì)應(yīng)多個(gè)系統(tǒng)噪聲協(xié)方差陣Q(n)和觀測(cè)噪聲協(xié)方差陣R(n)，在建立初始值及推算過程中，最能貼合目標(biāo)實(shí)際狀態(tài)的濾波權(quán)值逐漸變大，最終接近或達(dá)到1，其它濾波模型逐漸弱化。這種濾波模型始終處于動(dòng)態(tài)推算過程，并按照目標(biāo)的運(yùn)動(dòng)規(guī)律周期性改善Q(n)及R(n)，以使?fàn)顟B(tài)及測(cè)量方程能動(dòng)態(tài)反映目標(biāo)的實(shí)際狀態(tài)，第i個(gè)卡爾曼濾波算法器被分配的概率為：

(31)

(32)

這種算法稱做多模型自適應(yīng)估計(jì)(MultipleModelAdaptiveEstimation,MMAE)。其局限性是計(jì)算量較大，但Q(n)和R(n)在濾波器組中并行運(yùn)算、相互獨(dú)立，因此較傳統(tǒng)濾波器更趨于穩(wěn)定,在雷達(dá)跟蹤中，鑒于目標(biāo)機(jī)動(dòng)飛行或完成特殊動(dòng)作等情景下，更能反映運(yùn)動(dòng)軌跡多樣化的實(shí)際情況。

4 總結(jié)

面對(duì)監(jiān)視源信息的多種不確定性，自適應(yīng)濾波算法是未來發(fā)展的趨勢(shì)。自動(dòng)化系統(tǒng)中的MME算法也可看作多模自適應(yīng)應(yīng)用的個(gè)例。ADS-B數(shù)據(jù)中心往往需要處理多個(gè)地面站信息，在航跡信息被多個(gè)地面站接收時(shí)，可以將此濾波算法作為研究方向，此方法具有收斂速度快、濾波精度高，可以提高跟蹤系統(tǒng)的精度和可靠性，得到的系統(tǒng)航跡同時(shí)可送至各運(yùn)輸及通航機(jī)場(chǎng)中小顯示設(shè)備，以提高其空管運(yùn)行效率。