基于ADS-B數據的一次雷達系統誤差配準方法

張召悅, 黃誠昊

(中國民航大學空中交通管理學院,天津,300300)

一次雷達和廣播式自動相關監視系統(automatic dependent surveillance broadcast, ADS-B)是空中監視的主要手段,兩者對目標監視原理不同,ADS-B只能監視合作目標,而雷達可實現合作與非合作目標的全空域監視。雖然在民航領域一次雷達已逐漸被二次雷達和ADS-B所替代,但在軍用領域及航空進近區域依舊是主要的監視手段。大部分一次雷達在出廠前會進行系統誤差的配準,但雷達的系統誤差會隨著時間的推移而累積[1],長時間使用后需再次進行系統誤差配準。因此研發一套高效的一次雷達系統誤差配準方法具有重要的現實意義。

當前主要的標校方案包括靜態有源標校方案與動態標校方案,前者采用信標機進行標校,其操作復雜且成本較高。后者以靶標(航空器)為合作目標,在靶標(航空器)上加裝差分全球定位系統(differential global position system , DGPS)以提升航空器的定位精度,并記錄DGPS的定位數據作為標校真值。該方法存在諸多不利因素,如調度靶標受航線、空域限制,導致校準周期長、花費成本高、實施難度大等。

在國際民航組織(international civil aviation organization , ICAO)的推動下,多數民航飛機均已裝備ADS-B IN/OUT設備并實時廣播本機的位置等參數,這些信息來自機載全球衛星導航系統(global navigation satellite system , GNSS),可靠穩定、易于選取[2],精度一般是普通雷達的3倍以上,可間接或直接作為真值來配準雷達系統誤差[3-4]。

當前雷達系統誤差的配準研究主要集中在2個方向:①基于雷達系統誤差在空間中分布均勻假設下的常值函數配準方法,如最大似然配準[5-6]、最小平方準則[7]、基于傅里葉變換的方法[8]、實時質量控制法[9]等;②基于雷達系統誤差在空間中分布不均勻假設下的配準方法,如曲面網格逼近方法[10]、擬合分區配準[11]、基于數據驅動的配準方法[12]等。上述方法大多忽略了隨機誤差對雷達系統誤差配準的影響[13],將整體誤差直接進行誤差分析與配準工作,對A/C模式的二次雷達配準十分有效。但是一次雷達配準時,由于存在較大的隨機誤差及軌跡剛性偏轉現象,傳統方法無法有效提升雷達監視精度。對此研究人員引入了卡爾曼濾波算法來減小隨機誤差并進行系統誤差配準,但傳統卡爾曼濾波算法通常需要加入偏置函數來抵消雷達系統誤差,偏置函數的精確性就決定了雷達系統誤差的配準精度。

針對上述問題,本文基于雷達系統誤差在空間中分布均勻的假設下提出了一種基于點云的跟蹤濾波配準方法。

1 航跡點時空對準方法

1.1 監視數據坐標轉化

卡爾曼濾波作為線性濾波器,需要在笛卡爾坐標系下進行計算。因此需將極坐標系(azimuth elevation range, AER)中的雷達航跡數據與空間大地坐標系(geographic coordinate system, GEO)中的ADS-B航跡數據轉化到同一空間直角坐標系(earth centered earth fixed , ECEF)中。

1.2 航跡點時空對準

因一次雷達數據缺少關鍵的航班號信息,無法直接進行系統誤差配準研究,需將雷達航跡與對應的ADS-B航跡相關聯。如果兩條航跡在一定時間內有足夠多的航跡點可被中心點聚類算法劃分為同一簇,則這兩條航跡表征的是同一個目標。基于此原理,提出了一種基于中心點聚類的航跡關聯方法,該方法以航空器航跡點的時間戳和航空器飛行高度為粗篩選條件:

步驟 1對將要進行中心點聚類的航跡點進行粗篩選,以雷達設備的掃描周期tperiod為篩選范圍,雷達航跡點高度與ADS-B航跡點高度的差值閾值為Hlim。假設選取的雷達航跡點時間戳為Tr,其對應的高度為Hr,ADS-B航跡點的時間戳為Ta,i,其對應的高度為Ha,i。篩選ADS-B航跡點集中時間戳滿足|Tr-Ta,i|

步驟 2采用中心點聚類算法,將步驟1通過粗篩選的所有ADS-B和雷達目標點進行聚類。在K-means算法中,異常航跡點將對算法產生較大影響,因而選取K-Medoids方法作為ADS-B和雷達航跡點的航跡關聯算法。

將雷達與ADS-B航跡點的集合作為中心點聚類算法的輸入,以不同ADS-B航班號的數量K作為輸出,使每個航跡點與其所屬簇的中心航跡點距離總和最小。隨機選擇K個航跡點作為初始中心點,反復用其他航跡點(非中心點)代替中心航跡點,直至迭代至最優的中心航跡點。

對迭代完成的K個航跡點簇進行分析,判斷Tr時刻該雷達點歸屬的航跡簇,并檢驗該簇中ADS-B航跡點是否都屬于同一個航班號,若檢驗結果為真,則將該雷達航跡點與對應的ADS-B航跡形成關聯關系。

步驟 3引入航跡關聯百分比(track match percentage,TMP)概念,如果某一時刻Tr,第i條雷達航跡的航跡點和第j條ADS-B航跡的航跡點可被中心點聚類算法分為同一簇,則第i條雷達航跡中存入1個第j條ADS-B航跡的航班號,一段時間內求出第i條雷達航跡中存入的各個航班號的占比即為TMP,TMP(i,j)=max{TMP(i,n),n=1,2,…,n}。當TMP(i,j)≥85%時,則認為第j條ADS-B航跡與第i條雷達航跡為關聯狀態。

步驟 4由于雷達數據和ADS-B數據更新周期不同,對雷達系統誤差分析時,應將兩者的航跡數據進行時間戳同步,獲得相同時刻標度下的航跡數據。由于ADS-B數據的更新頻率約為雷達數據更新頻率的2～4倍,因此以雷達數據的時間戳為基準,對ADS-B的航跡數據進行插值,完成航跡數據的時間戳同步。3次樣條插值法的插值次數較低,待定系數易求解,插值速度較快,且在局部擬合中插值函數的光滑性較好,故選擇該方法[14]。

1.3 雷達誤差分布模型

雷達數據和ADS-B數據進行時空對準后,將ADS-B航跡數據轉化至以雷達為中心的極坐標系中,計算雷達與ADS-B數據的測角誤差(誤差統計為直方圖)。對測角誤差直方圖進行分析,結果見圖1。測角的整體誤差(隨機誤差+系統誤差)并不呈現為正態分布,文獻[15]發現該誤差服從某一自由度的t分布。但整體貼合度依舊不高,因此結合雷達系統誤差和隨機誤差的混合影響因素,考慮采用混合高斯分布模型反映直方圖的特征:

(1)

式中:N(μk,Σk)為正態分布;αk為混合高斯分布的權值;K為單高斯模型數量。在圖1中可以看出,混合高斯分布模型對雷達測角誤差直方圖有較好的擬合效果。整體誤差呈現混合高斯分布,反映出雷達誤差存在多個誤差源,傳統的配準方法忽略了隨機誤差而重點考慮系統誤差的方法,在大誤差一次雷達的配準研究中并不適用。

圖1 測角誤差分布擬合圖

將誤差源分離,并對不同誤差源采用不同的配準方法較為合理,因此本文提出了一種改進的軌跡跟蹤濾波配準方法。軌跡跟蹤濾波算法可減小大誤差一次雷達測量過程中隨機誤差對修正雷達系統誤差的影響,減小雷達航跡中因隨機誤差而產生的航跡波動,而點云配準算法則可預先計算軌跡跟蹤濾波算法中的偏置函數或補償函數,進而提升跟蹤濾波算法在進行雷達配準時系統誤差的配準精度。

2 改進的雷達軌跡跟蹤配準算法

2.1 偏置函數的確定

為確定卡爾曼濾波算法進行誤差配準時最優偏置函數,本文將點云配準算法中的迭代最近點法(iterative closest point,ICP)引入卡爾曼濾波算法中,ICP算法通過不斷迭代尋找雷達觀測數據與ADS-B參考數據之間最小平均誤差的剛性變換矩陣,從而獲得最優的偏置函數。

為避免因位姿初始值選擇不當導致ICP算法陷入局部最優的情況,引入隨機抽樣一致性算法(random sample consensus,RANSAC)進行粗匹配。RANSAC算法通過不斷選擇隨機點對,并記錄最優的匹配點對集合,獲得可靠的變換矩陣作為ICP算法中的位姿初始值。

基于點云配準算法獲得偏置函數主要分為5步:①確定一段歷史數據中的雷達航跡點(待配準點云)與ADS-B航跡點(基準點云);②確定RANSAC算法的迭代次數、兩點的閾值距離、內點的占比,采用RANSAC算法將雷達航跡與ADS-B航跡進行粗配準,輸出初始變換矩陣;③將RANSAC粗配準后輸出的變換矩陣作為ICP算法的初始值,ICP配準中通過K-D樹查找雷達航跡與ADS-B航跡之間的對應點,并迭代求解目標函數;④最終輸出ICP精配準后的位姿變換矩陣,將得到的位姿變換矩陣作為粗偏置函數;⑤將一段歷史數據中所有的ADS-B航跡與雷達航跡按上述方法求解粗偏置函數,并存入集合,利用遺傳算法求解偏置函數集合中最優的偏置函數。

2.2 跟蹤濾波算法

在進行機動目標跟蹤濾波時,為快速適應機動目標,需實時切換不同的運動模型,然而單模型的運動模式是固定的,不具備實時調整的可能性,因此使用固定結構交互式多模型算法時,建立多個模型并存儲在模型集合中,每個模型代表一種飛機的運動模式。為了保證跟蹤模型可以更加匹配機動目標,需要在模型集合中儲存大量飛機運動模型,然而這勢必會導致跟蹤精度的降低[16]。

因此,需選取變結構交互式多模型算法來對雷達探測目標進行跟蹤濾波,減少隨機誤差帶來的影響。該算法的關鍵在于能夠動態地調整模型的結構和參數,以適應目標在不同模式之間跳變的情況。因此,本文采用基于期望增加模型的變結構交互式多模型算法(expected-mode augmentation for multiple-model estimation,EMA),來抑制雷達航跡中的隨機誤差。

假設飛機運動模型集為M,通過擴展期望模型,找到更好的描述目標運動狀態的運動模型C,使得模型集E=M∪C,模型E稱為期望模型集。一般模型集Mk是隨著時刻的變化而變化的,通過應用模型自適應算法,新的模型集為:Mk=Ek∪(Mk-1-Ek-1),Mk為k時刻的模型集,Ek為k時刻的擴展期望模型集,Mk-1為k-1時刻的模型集,Ek-1為k-1時刻的擴展期望模型集。

模型集自適應分為兩部分,第一部分為模型之間的切換,模型切換過程中假設M中各模型之間的相互轉換服從齊次的馬爾科夫過程;第二部分為期望模型擴展,期望模型擴展將生成一個新的期望模型來代替原有的期望模型。

模型集自適應步驟如下:

步驟1計算有效模型的模型概率和似然函數。

(3)

(6)

步驟3卡爾曼濾波。

將RANSAC-ICP點云配準獲得的最優偏置函數b,代入單一模塊的Kalman濾波器中。

x(k)=Fx(k-1)+b+w(k-1)

(7)

式中:x(k)為系統的狀態向量;F為狀態轉移方程;b為偏置函數;w(k-1)為過程噪聲。

3 算例分析

為驗證上述算法的性能,本文進行了仿真實驗和實測驗證。圖2為雷達誤差配準的流程圖。

圖2 雷達誤差配準流程圖

3.1 仿真驗證

為驗證本算法的可靠性,采用計算機仿真的方法對本文文獻[6]、文獻[7]和Alpha-beta濾波器的算法進行對比。ADS-B數據源于華東某地區的實測數據,雷達仿真數據通過ADS-B航跡數據增加剛性偏移量產生。為保證實驗的可靠性,建立2組仿真場景,分析算法的性能。

3.1.1 場景 1:二次雷達A/C模式的監視場景

X方向添加系統誤差200 m,Y方向添加系統誤差200 m,Z方向添加系統誤差400 m。X、Y、Z方向分別添加均值為50 m的高斯噪聲。圖3為該場景的仿真結果。

(a)配準結果

(b)方位角誤差

(c)俯仰角誤差

(d)測距誤差

3.1.2 場景 2:大誤差一次雷達的監視場景

X方向添加系統誤差200 m,Y方向添加系統誤差200 m,Z方向添加系統誤差400 m。X方向添加均值為50 m的高斯噪聲,Y方向添加均值為100 m的高斯噪聲,Z方向添加均值為100 m的高斯噪聲。添加隨機生成的旋轉誤差0.2°～0.3°。圖4為該場景的仿真結果。

各方法計算效率分別為:文獻[6]為0.20 s,文獻[7]為2.08 s,α-β濾波為0.035 s,本方法1.88 s。此處的計算時間為處理約90條軌跡所花費的時間。

在場景1的仿真結果中,4種配準方法在配準的精度上都具有較好的表現,尤其是文獻[6]、文獻[7]中的方法,在對A/C模式二次雷達的配準中,都具有較高的配準精度。本文提出的方法在測距精度上略低于前2種方法,但在測角精度上具有較好的配準效果。

(a)配準結果

(b)方位角誤差

(c)俯仰角誤差

(d)測距誤差

在場景2的仿真結果中,文獻[6]、文獻[7]和Alpha-beta濾波器在配準的精度上均有降低,尤其是文獻[6]中的方法因其不具備濾波器功能,在出現大誤差場景時精度出現大幅度降低。文獻[7]中的方法在大誤差的場景下,精度也出現了降低,主要歸因于大誤差的場景往往伴隨著雷達軌跡出現剛性偏轉,對剛性偏轉無法較好的配準。本文提出的配準方法能在大誤差的場景下依舊保持較高且穩定的配準精度,主要得益于卡爾曼濾波器中偏置矩陣的優化與卡爾曼濾波器本身對于噪聲軌跡的優化。

在運算時間效能上,Alpha-beta濾波的方法在測距精度上稍弱于其他3種配準方法,但在運算時間效能上,Alpha-beta濾波器作為一種構造簡單的濾波器,在時間性能上大幅度領先于其他方法,更加適合于工程應用。文獻[7]的方法大部分時間用于實時計算變換矩陣H,而本文提出的方法通過歷史軌跡數據率先進行了最優偏置矩陣的計算,因此本文提出的方法比文獻[7]方法花費的時間短。

3.2 實測分析

本次實測使用的雷達精度參考標準如下:方位角精度為0.3°,測高誤差為800 m,測距誤差為60 m,是一種大誤差的一次雷達。該雷達測角誤差較大,且探測的目標距離在100 km左右,因此俯仰角度的輕微誤差都會導致幾百米的測高誤差。

本實驗首先對該雷達的歷史數據進行分析,利用RANSAC-ICP算法獲得卡爾曼濾波算法的最優偏置函數,同時復現文獻[6]、文獻[7]、α-β濾波器的雷達配準算法與本算法作對比。其中文獻[7]采用了一種實時動態的卡爾曼濾波算法,采用變換矩陣H來進行偏置或狀態補償。文獻[6]采用了最大似然估計來進行雷達系統誤差的配準。Alpha-beta濾波器計算效率高,是一種偏向工程的配準方法。圖5為某一時段的實測雷達數據配準結果。

(a)實測配準結果

(b)測距精度

(c)俯仰角精度

(d)方位角誤差修正結果

在俯仰角的配準中,不具備濾波器功能的文獻[6]與濾波性能較低的α-β濾波器配準結果并不理想。這是因為該雷達的大誤差主要反映在俯仰角測量精度較低,導致俯仰角噪聲和隨機誤差較大。而基于卡爾曼濾波算法的2類方法都保持了理想的配準結果,其歸因于卡爾曼濾波算法對于隨機誤差具有較好的抑制作用。

在方位角的配準中,本文方法與文獻[7]方法相比,具有更好的配準精度,這是因為文獻[7]中獲得的變換矩陣H來自同一時段下周圍航空器的航跡偏移數據,該小數據量的估計方法,提升了配準方法的實時性,但犧牲了配準方法的準確性。本文方法一部分基于文獻[7]的思想,利用卡爾曼濾波降低隨機誤差對雷達配準的影響,同時又對此算法提出改進,提前對雷達的歷史數據進行分析,借助點云配準思想來確定最優的偏置矩陣,在配準過程中實時性與準確性都優于傳統卡爾曼濾波算法。

3.3 結果分析

對該型號雷達收到的大量雷達航跡實驗驗證的結果見表1。處理效率為每處理100條軌跡所花費的時間。

實測分析符合雷達精度的參考標準。從表1中可以得出α-β濾波器雖然計算速度最快,但同時也損失了一部分精度作為代價。文獻[6]的優勢在于不需要歷史數據做支撐,可直接投入雷達配準的應用中,但在大誤差一次雷達的配準中修正精度和修正速度均低于本文的方法。文獻[8]的修正方法在修正波動較小的方位角時,修正效果較好,且計算速度較快,但因不具備濾波器功能,在波動較大的俯仰角修正時,得不到較好的修正結果。通過本文提出的方法修正后該雷達俯仰角精度提升至0.04°,方位角精度提升至0.07°,相較于文獻[7]、文獻[6]和α-β濾波器的方法,精度得到較好的提升。同時與同類型的文獻[7]方法相比,具有較高的計算效率。因此本文提出的方法在配準大誤差一次雷達時,具備更好的配準效果。

表1 各方法修正精度結果對比

4 結語

針對一次雷達系統誤差配準問題,本文首先采用中心點聚類對雷達數據和ADS-B數據進行航跡匹配,其次分析雷達誤差的分布特征,發現雷達誤差呈現混合高斯分布的特性,驗證了雷達誤差由多種誤差源構成。同時針對傳統卡爾曼濾波配準方法忽略了偏置函數或補償函數的問題,本文提出了改進EMA-IMM的配準方法。該方法利用RANSAC-ICP算法來確定最佳的卡爾曼濾波偏置函數或補償函數。研究結果表明,該算法可以提高大誤差一次雷達監視定位的精度,具有實際的工程應用價值。后續的研究將集中于ADS-B數據的精度分析,以滿足高精度雷達標定需求。