一種結合多重插補的ADS-B數據篩選算法*

鄒文華,左 諦,金開研

(民航數據通信有限責任公司 空管數據部，北京 100191)

0 引 言

廣播式自動相關監視技術(Automatic Dependent Surveillance-Broadcast，ADS-B)，是利用機載星基導航和定位系統獲得定位信息，通過空-空、地-空數據鏈通信以廣播形式發送數據的航空器運行監視技術[1]。ADS-B地面站接收其覆蓋范圍內的航空器下發的空地DF17格式的ADS-B報文[2]，地面站對此類報文進行整合、組裝后送出CAT21格式的ADS-B報文，報文里包含位置、時間、速度、航空器身份等信息[3]。雖然原則上地面站的處理過程遵循ED129系列協議[4]，但是各家廠商對協議內容有不同解讀，因此處理過程存在差異，并且不同地面站在硬件配置、運行狀態和輸送數據的網絡延遲等方面也均存在不同。上述原因最終會導致不同地面站送出的報文在數據項完整性、數據質量和數據實時性等方面存在差異，如果不做處理便直接送出，則會出現大量的重復報文、少量的回跳情況，以及各種數據項不全或低質量報文的數據，因此必須對數據進行處理。

隨著國內ADS-B的不斷發展[5]，在西部的無雷達地區，ADS-B逐漸作為主用且唯一的監視手段給管制提供服務，地面系統需對新監視技術的監視數據進行融合處理。對于傳統雷達數據的處理，地面處理系統對雷達目標的融合一般在時間和空間統一后進行加權融合算法，同時結合各種濾波方式，如卡爾曼濾波，消除雷達系統的隨機誤差，平滑航跡[6]。在有雷達數據參與融合的算法一般均采用與多雷達融合類似的思想，即加權和濾波結合的方法[7-9]。而在無雷達監視且新監視技術作為主用且唯一監視手段的地區，使用雷達區域的處理方法處理新監視數據并不合適，特別是對于ADS-B技術。ADS-B與雷達相比，有以下兩個本質的不同：一是ADS-B報文中最重要的位置信息由機載GPS獲得后下發，任何一個地面站獲得報文后解出的位置誤差理論上應該都一樣，與航空器距離地面站的距離無關，不需要對位置進行加權計算；二是機載下發的GPS位置信息本身已經在機載端進行了濾波處理[10-11]，ADS-B地面系統解析輸出的位置信息平滑度較高。因此，多地面站覆蓋時對多ADS-B數據的融合不用套用傳統的多雷達融合算法。

為此，本文利用ADS-B報文中的位置信息具有唯一確定且相對較平滑的特點，采用選優的方式從多條報文中挑選質量最優的ADS-B報文作為篩選后的數據，對于篩選后不合格的數據，根據實際情況，在數據嚴重缺失或條件極度不理想的條件下，可采用多重插補法提升其完整性。通過對實際多個航班ADS-B數據的算法驗證，證實了算法能有效去除重復和回跳報文，并優化提高數據整體質量。

1 航跡建立

本算法為周期型輸出的篩選算法，間隔以1 s為例，在篩選程序啟動后，收集每一秒內的所有報文，查看是否已經建立航跡，如果未建立航跡，則先建立航跡；如果已經完成航跡的建立，則進入篩選程序。對于尚未建立航跡的報文，應用以下方法確立可信賴的航跡。

航跡建立和處理是將同一目標的點跡連成航跡的處理過程，高效且可靠的航跡建立是保障航跡穩定更新的前提。一般直接起始法的基本原理是，對于連續收到的目標位置信息Pi(i=1,2，…，N)，如果其中有M(M

圖1 航跡建立流程圖

Step1 收集N條ADS-B報文，記為{REPORTi}(i=1,2，…，N),N為報文數量，為自定義參數，一般取大于3的數，本文中取5。每條報文中含位置(即經緯度)、時間和速度。

Step2 采用倒敘遍歷方式，即時間值從大到小順序，從{REPORTi}中挑選時間差大于0.4 s且小于4 s的3條報文，記為{RPi}(i=1,2，…，N)。0.4 s和4 s這個兩個時間差選擇的標準主要基于業務考慮，ADS-B位置更新的頻率為0.4～0.6 s[2]，因此0.4 s ADS-B位置更新的最小時間間隔，如果前后報文時間小于0.4 s，則有理由相信這兩條報文中的位置可能是相同的，不適宜用來做航跡建立判斷。雷達更新的頻率為4～12 s。如果前后兩個位置的時間差達到了ADS-B自身更新頻率的10倍，接近雷達的更新頻率，則認為此報文串連續性有欠缺，不適宜用來做航跡建立判斷。

Step3 對{RPi}(i=1,2,3)中的3條報文進行互相驗證。建立階段要滿足3條連續報文都是可靠的，方法是位置和速度相匹配，且不能超出最大位置差，即判定3條報文是否滿足下述公式：

(1)

(P1與P2點的計算方向差-P2與P3點的計算方向差)<120°。

(2)

式中:Vmax為航空器最大可能速度。

如果3條報文滿足式(1)和式(2)，則航跡建立成功，輸出RP1作為航跡起始點；如果3條報文不滿足式(1)和式(2)，則航跡建立不成功，重新遍歷尋找滿足Step 2的3條報文，找到后執行Step 3，如果遍歷完畢沒有找到滿足Step 2的3條報文，則說明當前{REPORTi}隊列中N個報文無法建立航跡，需要釋放清空此隊列，重新收集報文進行航跡建立過程，即回到Step 1。

2 ADS-B數據篩選算法

2.1 多約束條件下的數據篩選

航跡建立成功后，按第1節中的兩個篩選步驟進行融合。首先挑選進入融合的報文。滿足公式(3)和(4)的報文可進入2.2節進一步篩選：

Tn+1<([Tn]+1),

(3)

Tn+1>(Tn+0.4)。

(4)

式(3)和(4)中的Tn表示第n個融合周期優選輸出的點的報文時間(Time of Message Reception for position，TMR)，[Tn]為第n個融合周期融合后的點的報文中TMR值的整數，Tn+1表示第n+1個融合周期收到的報文中的TMR時間。

公式(3)中的1 s，是基于本文融合算法的周期為1 s；公式(4)中的0.4 s，即進入融合的報文時間必須比上一條報文時間大0.4 s，原因見1.1節中關于0.4 s的說明。

按照上面挑選標準，如果只挑選出了一條報文，且質量達到要求，則作為第n+1個融合周期融合點輸出；如果質量未達到要求可采用多重插補法(Multiple Imputation,MI)對缺失數據進行補全達到質量要求后輸出。如果有多條報文，則進入2.2節，進行優選融合。

用圖2舉例說明。U為報文接收時間軸，T為報文TMR時間軸，假設第一個周期輸出的融合點為U1/T1點，下個融合周期的時間范圍為[[U1]+1,[U1]+2)，[U1]為U1值的整數，那么在這個融合周期時間段內，這架飛機共收到4條報文，這4個點的接收時間和TMR時間(格式記為U2/T2)相對上個融合輸出點分別是(U1+0.8)/(T1+0.1)、(U1+1.0)/(T1+1.5)、(U1+1.2)/(T1+1.6)和(U1+1.4)/(T1+0.6)，再根據式(3)和式(4)進行挑選。

圖2 進入融合程序的報文挑選示例圖

假設上一個點的融合點U1/T1的U1值為0.4，T1值為0.1，那么可以得出(U1+0.8)/(T1+0.1)點不滿足式(3)和式(4)，(U1+1.4)/(T1+0.6)不滿足式(3)，(U1+1.0)/(T1+1.5)和(U1+1.2)/(T1+1.6)滿足條件，可按照2.2節對這兩個點進行進一步質量評估，挑選出質量高點，作為此融合周期的融合點輸出。

2.2 根據數據質量指標的篩選

2.1節中，第n+1個融合周期挑選出來多條報文，如何從這些報文中挑選出最優的報文，則需要從多方面進行評估。本文選取兩個方面評估，分別是報文中數據項的完整性和報文中位置的質量。數據項的完整性是報文中用戶關心的數據項的存在情況，用戶關心的數據項越豐富，說明報文的數據項完整性越好。報文位置的質量由航空器下發報文中攜帶的位置數據項不確定性(Navigation Uncertain Category，NUC)進行評價。NUC表征的是ADS-B報文中下發的位置的精度和完好性，取值范圍為0～9，數值越大，位置的精度和完好性即數據質量越高。通過式(5)計算參與融合的點的質量水平參數值Qi，挑選Qi值最大的點作為融合后的點:

Qi=Witm×Qitm+Wnuc×Qnuc。

(5)

式中：Qitm為數據項完整性標識參數；Qnuc為報文位置質量標識參數；Witm和Wnuc分別是數據項完整性標識參數和報文位置質量標識參數的權重，兩個權重值的和為1。

通過計算所有報文的Qi，選取Qi值最大的報文作為此融合周期的優選融合點輸出。

2.2.1Qitm的計算公式

在[Ti,Ti+1)時間段內，參與優選融合的某個報文k的數據項完整性參數Qitm(i,k)由報文中重要數據項的存在性確定。根據民航局空管局在2018年頒布的《民航空管ADS-B及S模式雷達數據應用實施計劃(修訂版)》[12]中確定的第一階段中對數據的要求，并結合航空公司在應用過程中的實際需求，本文選擇氣壓高度、航班號、二次代碼、地速和時間5個數據項的有無表征完整性分項參數，各個參數根據專家評級法取合適的權重，但各分項權重之和需等于Witm，即

θ1+θ2+θ3+θ4+θ5=Witm。

(6)

式(6)中的θ1、θ2、θ3、θ4和θ5分別是氣壓高度、航班號、二次代碼、地速和時間5個數據項的完整性分項標識參數。

2.2.2Qnuc的計算公式

在[Ti,Ti+1)時間段內，參與融合的某個報文k的NUC參數Qnuc(i,k)的計算公式如下：

Qnuc(i,k)=NUCk/9。

(7)

由于NUC是一個0～9的整數，需要做歸一化處理后用于表征ADS-B報文質量參數。

3 數據極度不理想情況下的補全

對于ADS-B報文中的各項指標，根據用途要求綜合專家評級法會有不同的權重需求，根據不同的權重要求，可以算出不同要求下的Qi，根據應用場景以及引用要求的差異，結合歷史數據確定閾值M，閾值M的確定需結合航空器飛行空域內的環境及流量共同確定。在篩選后，若所篩選Qi>M，則判定報文數據合格；若所篩選Qi

3.1 基于MI算法對缺失數據進行補全

多重填補方法分為三個步驟:

Step1 為每個空值產生一套可能的填補值，這些值反映了無響應模型的不確定性；每個值都被用來填補數據集中的缺失值，產生若干個完整數據集合。

Step2 每個填補數據集合都用針對完整數據集的統計方法進行統計分析。

Step3 對來自各個填補數據集的結果進行綜合，產生最終的統計推斷，這一推斷考慮到了由于數據填補而產生的不確定性。該方法將空缺值視為隨機樣本。

3.2 基于DA方法的MI插補方法

對于上述篩選過程，如果數據質量達不到要求，則需要對數據進行插補處理。考慮到ADS-B數據類型特點，即報文之間彼此獨立，每條報文之間的特征值也是彼此獨立，所以采用基于DA的MI算法對數據進行插補。數據擴張(Data Augmentation，DA)算法是馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)算法之一，特點在于可以處理任意缺失模式[14]。

DA算法可以通過以下步驟實現:

從Pr(Ymis|Yobs,θ(t))中抽取Ymis(t+1), 從Pr(θ|Yobs,Ymis(t+1))中抽取θ(t+1)，重復該過程多次，這樣就建立了一條馬爾科夫鏈：Ymis(1)，θ(1)，Ymis(2)，θ(2)……而且收斂于P(Ymis,θ|Yobs)。其中，θ表示條件概率的參數，Yobs代表觀測值，Ymis代表缺失值，P(Ymis,θ|Yobs)則表示在觀測值已知的情況下對缺失值進行預測的條件分布。DA算法估計的目的是從收斂的分布中隨機抽取Ymis值，替代缺失數據。

3.3 基于DA方法的MI插補方法有效性論證

在通過文獻查詢和理論分析論證了插補算法在ADS-B數據補全方面的適用性之后，本節通過使用人工生成模擬數據的方法驗證插補方法的有效性。選取地面站系統中2019年6月20日航班號為CES5810的1 000個完整ADS-B數據，人工產生缺失率分別為5%、10%、20%、30%、40%、50%，缺失類型為完全隨機缺失的缺失數據，再使用基于DA方法的MI插補方法對缺失數據進行插補。而在上述文本中選擇的氣壓高度、航班號、二次代碼、地速和時間5項參數中，航班號、二次代碼屬于定性數據，由飛行員輸入和設置，無法通過已有數據建立預測分布，而其中的時間項由于對數據精度要求過高，不適用于用插補法進行插補，故僅對氣壓高度、地速進行插補。

表1總結比較了原數據與插補后數據的均值、標準差以及皮爾森相似度[15]，用以衡量插補前后數據的相似度。在R語言中完成上述過程[16-17]。由表1可知，插補后的數據的均值與方差與原數據相比均沒有太多差異，而插補后的數據與原數據的皮爾森相似度隨著缺失率的增長而緩慢下降，但仍然保持著較高相似度，證明了基于DA方法的MI插補方法的有效性。

表1 插補前后對照表

4 篩選算法的仿真過程

4.1 仿真流程

本文取本公司地面站系統中2019年6月19日3個航班的所有地面站數據進行融合，驗證融合及選優效果。所選航班信息、跨時、接收到報文的地面站數量見表2。

表2 模擬數據信息表

取氣壓高度、航班號、二次代碼、地速和報文時間5個數據項的有無表征報文完整性的分項參數，取NUC值的歸一化后表征報文質量參數，各參數的權重設定如表3所示。

表3 權重參數表

表3中權重值在滿足局方文件中對數據項的存在性要求的條件下，兼顧了用戶的實際使用需求。NUC、氣壓高度、航班號、二次代碼的權重較高，NUC和氣壓高度權重高是對位置精度、可靠性和高度存在性的要求，體現了監視技術的核心是四維航跡信息。航班號和二次代碼是表征航空器身份的信息，在管制系統中是ADS-B航跡與飛行計劃進行相關的關鍵信息，重要性不言而喻。報文時間項和地速項存在性相對較弱，但也是體現了ADS-B數據自身的優勢。ADS-B數據更新率高、及時性強并自身攜帶高穩定性的速度信息，無需地面系統計算，優勢明顯。

4.2 仿真結果

根據第2節的參數和權重對這3架航空器進行組合優選融合處理，對融合前后報文數量、報文NUC質量平均值、數據完整性指標、整體質量進行評估，結果如表4所示，可以確定算法去除了大量重復的報文，并且整體質量參數有提高。

表4 仿真結果

結合該運行場景對數據的要求，選取閾值M為0.85，由表4可知，經過篩選后，除地址碼為780E81的航班外，另外兩個航班的Qi均大于M，符合質量指標要求，而對于地址碼為780E81的航班，則需使用上述MI插補方法進行補全。使用上述DA方法的MI插補方法對極端不理想數據地址碼為780E81的航班進行插補。

4.3 插補過程模擬

基于DA算法的算法原理，需要基于已有數據生成預測分布，如3.3節中所述，在上述文本中選擇的氣壓高度、航班號、二次代碼、地速和時間5項參數中，航班號、二次代碼屬于定性數據，由飛行員輸入和設置，無法通過已有數據建立預測分布，可由人工填寫，故依舊僅對氣壓高度、地速進行插補。地址碼為780E81的航班的ADS-B數據(氣壓高度、地速)的數據缺失情況如圖3所示。

圖3 航班 “780E81”ADS-B數據缺失情況

由圖3可知，地址碼為780E81的航班的ADS-B數據的缺失類型為完全隨機缺失，且由飛機航行以及氣壓高度測量原理可知，氣壓高度與地速之間存在線性關系，所以在DA算法中使用線性回歸預測模型，基于已有值，建立預測分布進行插補，插補結果如圖4所示，可見經過插補氣壓高度及地速的缺失率下降為0，Qi值顯著增加至0.890 9，且大于閾值M，符合數據質量要求。

圖4 插補前后缺失情況對比圖

5 結束語

本文所提算法可有效解決不理想條件下ADS-B數據中的重復及回跳現象，確定數據質量指標以快速對報文數據進行篩選，在惡劣環境下亦有良好應用。對于缺失數據，使用多重插補法進行補全，極大地提高了數據完整性。

但同時，在上述流程中，無論是篩選還是插補依舊存在信息的浪費以及誤差，想要從根本上解決數據的質量問題，仍應該不斷發展航空通信導航技術，提高地面站及航空器下發數據質量。