李 劍
(中國民用航空呼倫貝爾空中交通管理站 內(nèi)蒙古 海拉爾 021008)
卡爾曼濾波算法是一種遞歸計算方法,最初是通過最優(yōu)估計的方法,來獲得離散數(shù)據(jù)的線性濾波,隨著計算機技術的提升,其應用滲透多個領域。在航空方面,卡爾曼濾波算法由于功能強大,不止能夠估計信號之前及當前狀態(tài),甚至可以據(jù)此估計信號的未來狀態(tài)。當前民航系統(tǒng)使用的空管二次雷達及ADS-B等定位設備,ATC終端使用的自動化設備等,也均使用卡爾曼濾波算法解決目標的航跡定位問題。
卡爾曼濾波算法的核心思想是提供一種算法,通過對目標狀態(tài)的估計來使過程的均方誤差達到最小,在估計狀態(tài)時,使用的測量數(shù)據(jù)也受到干擾噪聲影響下進行。[1]在做推導時,隨機信號的目標狀態(tài)方程使用x(n)表示;而測量數(shù)據(jù)k(n)跟隨目標也為線性狀態(tài),其中觀測噪聲服從零均值分布,則卡爾曼濾波算法的模型方程為:
(1)
式中,w(n)和v(n)為目標及測量白噪聲,a和c為常數(shù),模值均小于1。
(2)
以上即卡爾曼濾波算法的標準形式,它把目標的狀態(tài)描述為一個離散的隨機過程,該過程中系統(tǒng)狀態(tài)及其測量均不可避免受噪聲影響??柭鼮V波算法實際是維納濾波的一種特殊表現(xiàn)形式。雷達跟蹤觀測所使用的濾波算法為一種矢量算法,以下介紹其推導過程。
1.2 卡爾曼濾波算法矢量方程
可以考慮兩種情況,一種是根據(jù)觀測數(shù)據(jù)同時推算估計若干個類似信號,例如q個信號x1,x2…,xq;另一種是得到一個高階自回歸過程,如一個q階的自回歸過程:
(3)
以上兩種情況下,如果引入矢量卡爾曼濾波算法,分析計算可以變得更加便捷。
以上情況中q個一階自回歸信號,我們設其為同時相互獨立的,則在n時刻,其取樣分別為x1,x2…,xq,它們的狀態(tài)方程表示為:
xi(n)=aixi(n-1)+wi(n),i=1,2,3…,q
(4)
同標量卡爾曼相似,使用wi(n)表示零均值白噪聲序列,它們相互之間相關。q個信號xi(n)使用矩陣表示為q維矢量x(n)=[x1(n)x2(n)…xq(n)]T,則我們可以將以上信號狀態(tài)方程簡化為以下矢量方程:
x(n)=Ax(n-1)+w(n)
(5)
分量w(n)是由wi(n)第q維構成的矢量,A是一個q階對角矩陣,分量為ai。
同理在n時刻測得k個數(shù)據(jù)k1(n),k2(n),…,kk(n),它們與xi(n)的關系為:
ki(n)=aixi(n-1)+vi(n),i=1,2,3…,k
(6)
其中,k≤q,測量噪聲使用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個測量方程的矢量方程簡化為:
k(n)=Cx(n)+w(n)
(8)
系數(shù)矩陣C是一個K×q矩陣。
已知以上標量卡爾曼濾波算法器遞推公式,同理可推導出矢量卡爾曼濾波算法器的相應公式為:
(9)
最小預測均方誤差:p(n)=A(n)ξ(n-1)AT(n)+Q(n)
(10)
(11)
(12)
2.1 雷達跟蹤
二次雷達系統(tǒng)是民航對飛行目標監(jiān)視使用最廣泛的設備,是通過地面詢問機對目標進行詢問,目標在收到詢問信號并對其進行應答來實現(xiàn)對目標定位的一種設備。一般對目標的定位采用極坐標的形式,即可用徑向距離及方位角表示目標的實際位置。
在跟蹤系統(tǒng)中,一般利用詢問與應答的時間延遲去計算目標斜距,并近似作為飛機徑向距離,而目標的方位角是在天線波束掃描到目標時雷達系統(tǒng)碼盤刻度。二次雷達天線轉(zhuǎn)動周期一般設置為4S,我們用表示其轉(zhuǎn)動周期,同時這也是二次雷達位置更新的時間間隔。
(13)
(14)
這在濾波算法中便等效于零均值白噪聲即:
(15)
(16)
得出以下矢量矩陣公式:
x(n+1)=Ax(n)+w(n)
(17)
在雷達天線半波束輻射范圍內(nèi),目標的方位角由天線碼盤及單脈沖信息獲得,而其距離則借助光在空間勻速直線傳播得出(此時忽略目標高度)。目標離開雷達天線波束后,雷達輸出點跡信息,此過程中將測量結果k1(n)和k2(n)寫成矢量矩陣形式:
k(n)=Cx(n)+v(n)
(18)
矢量w(n)和v(n)各自的方差為自相關矩陣Q(n)和R(n):
(19)
(20)
2.2 初始值確定及影響濾波算法的因素
在跟蹤算法中,卡爾曼濾波算法器需要得到增益矩陣Gn的初值。為此,在點跡或航跡形成之初,需確定均方誤差ξ(n)的初始值。一般點跡或航跡的初始化n=3時形成,及利用前兩次n=1和n=2測量得到的距離及方位角參數(shù),得到四個數(shù)據(jù)k1(2),k2(1),k1(2),k2(2),并對四個數(shù)據(jù)k作以下估計:
(21)
代入以上測量方程得:
(22)
(23)
將狀態(tài)方程代入得:
(24)
綜合以上一系列公式,可得:
(25)
對均方誤差做出以下計算:(2.14)
(26)
即可確定均方誤差ξ(n)的初始值為:
(27)
由于海拉爾二次雷達建設受周邊建筑影響,雷達站周邊環(huán)境急劇惡化,近距離鏡像及假目標出現(xiàn)頻繁,反射使得在對初值、測距及測角均方差值賦值時偏離實際較大,我們知道,鏡像的出現(xiàn)是存在較為固定的反射體,這時假目標形成的濾波航跡較為穩(wěn)定,雷達本身很難識別,技術人員對功率可編程、STC曲線進行了調(diào)節(jié),降低了近距離物體反射的影響,情況得以好轉(zhuǎn),但依舊不能解決零散假目標形成的點航跡問題。
對于鏡像,可以使用抗反射的方法處理,但由于附近居民樓影響,零散雜波較為普遍,就需要別的措施解決雷達濾波的方差賦值影響,SELEX二次雷達在航跡跟蹤門限設置中,提供了航跡速度初始值及跟蹤維護的門限設置,根據(jù)日常目標在高空巡航及中低空進近的常規(guī)速度,我們對以上參數(shù)做出門限設置,雜波假目標的狀況得到明顯改善。
圖1 二次雷達航跡濾波初始值及跟蹤維護門限設置
2.3 卡爾曼濾波算法器在ADS-B中的應用
廣播式自動相關(ADS-B)的基本原理是飛機(航空器)通過空地、空空數(shù)據(jù)鏈,采用周期無方向性廣播方式,自動發(fā)送和接收記載設備所提取的監(jiān)控信息。[2]ADS-B數(shù)據(jù)更新率及準確率均優(yōu)于二次雷達。近年來,海拉爾處于支線機場航空運輸快速增長,通航機場不斷增多的現(xiàn)狀,中低空空域機構日漸復雜,需要大量的監(jiān)視源設備填補管制監(jiān)控盲區(qū)。目前,中國民航正在全力推行民航客機二次雷達應答機加裝ADS-B應答功能,而且ADS-B建站成本低廉,可很好提升海拉爾管制區(qū)監(jiān)視冗余備份能力。
與雷達跟蹤類似,根據(jù)上述推算,ADS-B的卡爾曼增益方程可以表示為[3]:
Gk=P(k)HT(k)[H(k)P(k)·HT(k)+R(k)]-1
(28)
航空目標的狀態(tài)方程即為:
(29)
對應的協(xié)方差更新公式為:
ξ(k)=[I-H(k)GK]p(k)
(30)
其中HT(k)為航空目標的觀測矩陣。
不同于二次雷達,ADS-B應用中需要跟蹤的目標狀態(tài)變量表示為T=(x,vx,y,vy,w,z);分別代表目標在X和Y方向的距離、速度值、轉(zhuǎn)彎速率及高度,其中表示實時估計值;模型建立后,使用同一目標不同周期的位置信息,來確定目標目前狀態(tài)、預計狀態(tài)以及相應的協(xié)方差矩陣。
2.4 空管自動化中的濾波算法應用
呼倫貝爾ATC系統(tǒng)目前使用民航二所的AirNet自動化系統(tǒng),接入數(shù)據(jù)將不同監(jiān)視源的同一目標的數(shù)據(jù)進行融合,系統(tǒng)將各監(jiān)視源的航跡融合成系統(tǒng)綜合航跡。并對每個周期的綜合航跡跟蹤濾波,使處理后的系統(tǒng)航跡更加接近真實航跡。其跟蹤使用IMM算法(InteractingMultipleModel,交互式多模型算法)。將目標在水平面內(nèi)的運動分三個分量。
IMM算法的核心思想是同時維護三個濾波[4],兩個標準CT模型和一個CV模型,分別對應左轉(zhuǎn)彎、右轉(zhuǎn)彎及直線運動,運算流程如圖2所示,F(xiàn)ilter0,Filter1,Filter2分別代表三種不同的運動模型,每個模型單獨濾波,初始值賦值后,所有狀態(tài)值的更新依賴于上次濾波運算結果,不斷更新概率模型,同時將更新作為下一次估計的初始輸入,從而實現(xiàn)對目標狀態(tài)的持續(xù)跟蹤。
圖2 IMM基本流程圖
監(jiān)視源權重取值與監(jiān)視源目標信息可靠性相關。目前二次雷達仍然為區(qū)域管制的主要監(jiān)視手段,但ADS-B數(shù)據(jù)處理中心匯接了多個地面站信息,在中低空覆蓋方面與二次雷達形成顯著優(yōu)勢,且數(shù)據(jù)更新率快,精度質(zhì)量高,有利于目標的平滑處理。根據(jù)實際情況,AirNet自動化系統(tǒng)對雷達及ADS-B靜態(tài)權重進行了以下分配:
圖3 AirNet自動化系統(tǒng)監(jiān)視源權重分配情況
目前海拉爾管制區(qū)內(nèi)新建多臺ADS-B地面站,地面站數(shù)據(jù)經(jīng)數(shù)據(jù)處理中心匯總后送入自動化系統(tǒng),考慮數(shù)據(jù)完整性,自動化系統(tǒng)接收所有NUCP取值范圍內(nèi)的數(shù)據(jù)進行數(shù)據(jù)融合,初期出現(xiàn)了較多假目標,經(jīng)分析及確認,先后有GPS干擾、RAIM失效(NUCP取值低)、與雷達信號目標分裂等情況,影響濾波算法準確性。由于GNSS不確定因素會造成NUCP值對應的位置精度及完好性與實際情況不符[5],經(jīng)與廠家溝通,增加了對ADS-B數(shù)據(jù)質(zhì)量驗證的設置選項,改善了航跡濾波效果。
實際中,濾波器的性能在一定條件下會受到Q(n)和R(n)陣的取值影響。無論單雷達航跡跟蹤、靜態(tài)融合或IMM算法,Q(n)取值變化與狀態(tài)濾波估計值依賴外推預測值的權重呈反比狀態(tài);對于R(n),濾波器在性能上主要跟隨觀測方程權重亦呈反比狀態(tài)。在大多數(shù)目標跟蹤中,目標狀態(tài)呈現(xiàn)多種變化趨勢,觀測值也受設備性能及現(xiàn)場環(huán)境影響產(chǎn)生較多變數(shù),Q(n)和R(n)的取值應該呈現(xiàn)動態(tài)變化的趨勢。
實際上,像自動化系統(tǒng)中對多監(jiān)視源的數(shù)據(jù)進行靜態(tài)及動態(tài)權重的分配,也可看作自適應卡爾曼濾波算法的一種表現(xiàn)形式,即便是擴展卡爾曼濾波算法的應用,也是在將非線性目標通過泰勒級數(shù)展開、樣本點獲取權重、平方根取樣等手段,達到目標跟蹤濾波算法線性化的效果。另外一種優(yōu)化方法,是同時建立多個卡爾曼濾波算法模型,對應多個系統(tǒng)噪聲協(xié)方差陣Q(n)和觀測噪聲協(xié)方差陣R(n),在建立初始值及推算過程中,最能貼合目標實際狀態(tài)的濾波權值逐漸變大,最終接近或達到1,其它濾波模型逐漸弱化。這種濾波模型始終處于動態(tài)推算過程,并按照目標的運動規(guī)律周期性改善Q(n)及R(n),以使狀態(tài)及測量方程能動態(tài)反映目標的實際狀態(tài),第i個卡爾曼濾波算法器被分配的概率為:
(31)
(32)
這種算法稱做多模型自適應估計(MultipleModelAdaptiveEstimation,MMAE)。其局限性是計算量較大,但Q(n)和R(n)在濾波器組中并行運算、相互獨立,因此較傳統(tǒng)濾波器更趨于穩(wěn)定,在雷達跟蹤中,鑒于目標機動飛行或完成特殊動作等情景下,更能反映運動軌跡多樣化的實際情況。
面對監(jiān)視源信息的多種不確定性,自適應濾波算法是未來發(fā)展的趨勢。自動化系統(tǒng)中的MME算法也可看作多模自適應應用的個例。ADS-B數(shù)據(jù)中心往往需要處理多個地面站信息,在航跡信息被多個地面站接收時,可以將此濾波算法作為研究方向,此方法具有收斂速度快、濾波精度高,可以提高跟蹤系統(tǒng)的精度和可靠性,得到的系統(tǒng)航跡同時可送至各運輸及通航機場中小顯示設備,以提高其空管運行效率。