靳慧斌，馮朝輝，張召悅，王志森

（中國民航大學(xué)a.安全科學(xué)與工程學(xué)院；b.空中交通管理學(xué)院，天津 300300）

廣播式自動相關(guān)監(jiān)視（ADS-B，automatic dependent surveillance broadcast）是航空器周期性廣播自身位置及相關(guān)信息，供空管對其進(jìn)行監(jiān)視的一種監(jiān)視技術(shù)[1]。ADS-B 數(shù)據(jù)傳輸過程中可能會因為定位設(shè)備故障、惡劣天氣、飛行姿態(tài)變化等因素導(dǎo)致部分報文數(shù)據(jù)缺失或數(shù)據(jù)跳點，使地面站接收的ADS-B 軌跡數(shù)據(jù)出現(xiàn)缺失點、離群點等異常點[2]。異常點的產(chǎn)生不利于系統(tǒng)和相關(guān)用戶對軌跡位置數(shù)據(jù)進(jìn)行應(yīng)用分析，并有可能產(chǎn)生軌跡點間隙過大、錯誤位置告警、圖形軌跡誤判等影響航空安全水平的事件。

為解決數(shù)據(jù)缺失問題，國內(nèi)外學(xué)者從不同方向展開了研究。鞠時光等[3]從連續(xù)性及保形性兩方面提出了用λ 參數(shù)三次樣條函數(shù)對等高線軌跡數(shù)據(jù)進(jìn)行修補(bǔ)論證，有效地解決了繪制地形圖過程中出現(xiàn)的等高線相交或缺失問題，但該方法只適用于線性函數(shù)。針對非線性數(shù)據(jù)插補(bǔ)修復(fù)，Zamani[4]采用二次樣條函數(shù)處理函數(shù)邊界非線性問題；Farhan 等[5]建立了一種隨機(jī)多重插補(bǔ)（RMI，random multiple interpolation）方法來克服非線性軌跡不確定性問題。但兩種方法對于處理大量數(shù)據(jù)具有局限性。為解決大數(shù)據(jù)插補(bǔ)局限性問題，Peeters 等[6]評估了處理丟失數(shù)據(jù)的不同技術(shù)，得到了對于大數(shù)據(jù)插補(bǔ)問題使用預(yù)測均值匹配的多重插補(bǔ)（MI，multiple interpolation）是最佳方法的結(jié)論。為提升ADS-B 數(shù)據(jù)完整性，鄒文華等[7]建立了一種引入時間序列的基于貝葉斯常均值模型的數(shù)據(jù)增廣算法，解決ADS-B 軌跡點回跳現(xiàn)象。此外為確保插補(bǔ)數(shù)據(jù)精度，Sun 等[8]構(gòu)建了一種使用最小二乘樣條逼近（LSB，least square b-spline）的軌跡數(shù)據(jù)重建算法，能夠避免因數(shù)據(jù)差值而導(dǎo)致插補(bǔ)誤差過大問題。

上述研究多側(cè)重從單一角度對數(shù)據(jù)缺失問題進(jìn)行插補(bǔ)分析，或以局部數(shù)據(jù)是否間斷為依據(jù)進(jìn)行插補(bǔ)驗證，未考慮數(shù)據(jù)可能存在異常點的問題。針對這一情況，本文采用多約束條件篩選出航空器軌跡缺失點并將其清洗插入，然后利用改進(jìn)后的K 最近鄰（KNN，K-nearest neighbor）方法檢測ADS-B 異常點，解決了ADS-B 樣本數(shù)據(jù)量過大帶來的聚類點權(quán)重分配不平衡而無法有效檢測異常點的問題，并對軌跡進(jìn)行多重回歸插補(bǔ)修復(fù)，進(jìn)一步提高軌跡數(shù)據(jù)完整性，從而保障航空安全監(jiān)管水平。

1 基于KNN 的軌跡數(shù)據(jù)異常點檢測

1.1 缺失點篩選

為確保經(jīng)度、緯度和高度使用同一量綱（m），需要利用米勒投影完成經(jīng)緯度和平面坐標(biāo)的相互轉(zhuǎn)換[9]。轉(zhuǎn)換過程如圖1 所示。

圖1 經(jīng)緯度米勒投影轉(zhuǎn)換Fig.1 Miller projection conversion of longitude and latitude

圖1 中：ARP（aerodrome reference point）為機(jī)場基準(zhǔn)點；L 為地球周長；W 為地球周長的平面展開，并將周長視為X 軸；H 為Y 軸，約等于地球周長一半（L/2）；mill 為米勒投影常數(shù)。

ADS-B 地面接收機(jī)在同一時間接收不同航空器機(jī)載ADS-B OUT 發(fā)出的報文信息是有限的[10]，中國民用航空局規(guī)定接收機(jī)標(biāo)準(zhǔn)是600 架飛機(jī)的解碼[11]，一臺ADS-B OUT 發(fā)射全類型報文為4 次/s，則ADS-B地面接收機(jī)每秒最少接收2 400 條全類型報文。

假設(shè)一臺ADS-B OUT 單位時間T 內(nèi)發(fā)射報文次數(shù)為C，則每條報文時間間隔

從航路進(jìn)入終端區(qū)的進(jìn)近航路相對地面高度一般為6 000 m 左右，無線電傳播速度為299 792 458 m/s，因此在終端區(qū)范圍內(nèi)，因距離而導(dǎo)致的接收延遲問題可以忽略。

一條全類型報文解碼平均時間長度為Ta。地面接收機(jī)接收報文并完成解碼，傳入終端顯示界面的時間間隔為Tu，則同一航班的兩條報文時間間隔應(yīng)為

一條軌跡點的報文信息包括經(jīng)度、緯度、高度、地速、航向角、時間，所以設(shè)軌跡點P={Plo，Pla，Pal，Psp，Pha，Pt}，軌跡點速度均為瞬時速度，因此將軌跡點之間的時間段位移作為平均速度位移，通過距離和時間差判斷是否有缺失點，其判斷公式如下

式中：時間差Δt=Pt，i-Pt，i-1；速度差Δv=Psp，i-Psp，i-1；Px為兩點之間距離；i 為軌跡點序號。其缺失點判斷如圖2 所示，其中Pvt為同一架飛機(jī)t 時間段內(nèi)飛行的距離。

圖2 缺失點判斷示意圖Fig.2 Diagram of missing point judgment

當(dāng)兩點之間出現(xiàn)缺失數(shù)據(jù)，應(yīng)插入空白軌跡點，其插入軌跡點數(shù)目滿足

1.2 基于標(biāo)準(zhǔn)歐氏度量的離群點篩選

KNN 算法具有最近鄰特性，對于非線性數(shù)據(jù)篩選可避免偏離數(shù)據(jù)與正常數(shù)據(jù)出現(xiàn)較高擬合[12-13]，因此適用于終端區(qū)ADS-B 數(shù)據(jù)檢測。

1.2.1 數(shù)據(jù)處理

在1.1 節(jié)中已經(jīng)完成缺失點的篩選和空白軌跡點插入。但在計算軌跡點相似度時，需要賦予空白軌跡點標(biāo)簽，避免出現(xiàn)因無標(biāo)簽而導(dǎo)致分類錯誤的問題，通過下三角距離矩陣對空白軌跡點進(jìn)行賦值，其矩陣如下

式中t 為軌跡點時間，當(dāng)矩陣中元素Wi，j＜0 時，對應(yīng)Pt，i為原數(shù)據(jù)；當(dāng)Wi，j＞0 時，對應(yīng)Pt，i=0。因ADS-B 軌跡在界面只顯示二維坐標(biāo)軌跡，所以只需要插補(bǔ)經(jīng)度和緯度。

1.2.2 距離計算

KNN 距離度量有兩種，一種是歐氏度量（Euclidean metric），另一種是曼哈頓度量。歐氏度量用于衡量各點之間的絕對距離。如三維坐標(biāo)點M1（X1，Y1，Z1）和M2（X2，Y2，Z2），則M1和M2兩點之間的距離為

在n 維空間中，兩個n 維向量A=（X1，X2，…，Xn）、B=（Y1，Y2，…，Yn）間的歐式度量為

曼哈頓度量是向量之間距離的范數(shù)。如向量A=（X1，X2，…，Xn）和B=（Y1，Y2，…，Yn）之間的曼哈頓距離為

選擇歐氏度量或曼哈頓度量主要取決數(shù)據(jù)維度特征[13]。對于高維向量，曼哈頓度量比歐氏度量計算距離的效果更好；對于低維向量，歐氏度量要優(yōu)于曼哈頓度量。因軌跡點P 經(jīng)過篩選，只考慮經(jīng)度、緯度和時間，所以本文選擇歐式度量。

歐式度量適用于樣本向量的各個分量度量標(biāo)準(zhǔn)統(tǒng)一的情形。ADS-B 軌跡樣本各分量權(quán)重相同，當(dāng)時間與經(jīng)緯度波動范圍量綱值差距較大時，可能會引起各個分量分類不平衡問題[14]。因此可以通過標(biāo)準(zhǔn)化歐氏度量（standardized Euclidean metric）解決Pt和Plo、Pla分量差值問題。標(biāo)準(zhǔn)化歐氏度量是針對歐氏度量缺點的一種改進(jìn)。標(biāo)準(zhǔn)化歐氏度量的思路為：既然數(shù)據(jù)各維分量的分布不一樣，那先將各個分量都“標(biāo)準(zhǔn)化”到均值、方差相等[15]。假設(shè)樣本集X 的均值為m，標(biāo)準(zhǔn)差為s，X 的“標(biāo)準(zhǔn)化變量”表示為

標(biāo)準(zhǔn)化歐氏度量公式為

由式（10）可以得到不同軌跡點之間標(biāo)準(zhǔn)化歐氏度量為

式中z 為分量度量，因KNN 遵循最近鄰原則，為避免同一航班軌跡點之間不同分量差值較大，導(dǎo)致線性方差誤差過大，因此需要限制軌跡點方差樣本數(shù)目i -n1≤i≤i+n2，n1為第i 點前一個不相鄰異常點位置，n2為第i 點后一個不相鄰異常點位置，將式（11）的變量進(jìn)行如下改進(jìn)，即

通過改進(jìn)后的KNN，檢測出缺失點和其他異常點的總集合{O1，O2，…，On}，然后通過程序遍歷總集合，將異常點清洗為空白軌跡點，便于后續(xù)插補(bǔ)。

2 基于多重插補(bǔ)的ADS-B 軌跡補(bǔ)全

多重回歸插補(bǔ)理論定義[16]為：設(shè)β 是從所要修復(fù)樣本數(shù)據(jù)X 中估計得到的參數(shù)向量，β 的先驗分布為f（β），使用貝葉斯準(zhǔn)則可得到其后驗概率分布為

式中Q 是分析樣本數(shù)據(jù)X 得到似然函數(shù)。

設(shè)Xp為X 中缺失信息的變量。pobs和pmis分別為p中的已存元素和缺失元素，并假設(shè)p 中的信息缺失機(jī)制為隨機(jī)丟失（MAR，missing at random）。

缺失數(shù)據(jù)樣本的參數(shù)向量期望值β|Xp，pobs為

Bernstein-von Mises 定理[17]將貝葉斯估計與經(jīng)典的最大似然估計聯(lián)系起來，總結(jié)這種關(guān)系如下

根據(jù)式（18），在二維多重插補(bǔ)隨機(jī)修復(fù)下（G=1，2），軌跡樣本數(shù)據(jù)的β|Xp，pobs計算公式如下

式（19）可以用于計算軌跡樣本的最終方差。通過從軌跡數(shù)據(jù)集中抽取多個隨機(jī)樣本，可以推斷出總體參數(shù)的均值和方差的無條件估計[18]。均值和方差-協(xié)方差矩陣計算如下

式中：Ex（β|X）為樣本數(shù)據(jù)X 中參數(shù)向量的期望；Varx（E（β|X））為Ex（β|X）期望的方差。軌跡缺失數(shù)據(jù)多重插補(bǔ)的步驟如下：

步驟1從要修復(fù)的軌跡數(shù)據(jù)隨機(jī)抽取部分樣本數(shù)據(jù)xg，g=1，2，…，G；

步驟2引入誤差項e 的方差，計算為

步驟5計算插入缺失值向量，如下

式中Z 為數(shù)據(jù)插補(bǔ)間可變性系數(shù)，在多次插補(bǔ)步驟結(jié)束時，可得到G 個完整的數(shù)據(jù)集。因航空器性能的影響，軌跡數(shù)據(jù)插補(bǔ)不應(yīng)以整個軌跡數(shù)據(jù)方差為系數(shù)，為得到更加準(zhǔn)確的插補(bǔ)值，采用軌跡分隔和軌跡排序判斷對多重插補(bǔ)方法進(jìn)行改進(jìn)，其過程如圖3 所示。

圖3 軌跡修復(fù)流程圖Fig.3 Flow chart of trajectory repair

因軟件界面只顯示航空器ADS-B 二維坐標(biāo)軌跡，所以只需要插補(bǔ)經(jīng)度和緯度即可。

3 算例分析

3.1 數(shù)據(jù)處理

本文選擇天氣狀況良好的時間段，排除天氣干擾因素，篩選出某機(jī)場終端區(qū)4 天內(nèi)軌跡坐標(biāo)數(shù)據(jù)，如圖4 所示。

圖4 終端區(qū)不同維度軌跡顯示Fig.4 Trajectory display of different dimensions in terminal area

通過匹配機(jī)型，從中隨機(jī)挑選2 個不同機(jī)型的航班，其軌跡信息如表1 所示。

表1 樣本軌跡信息Tab.1 Information of sample trajectory

樣本軌跡數(shù)據(jù)通過米勒投影轉(zhuǎn)換和缺失點篩選得到結(jié)果如表2 所示。

表2 缺失點信息Tab.2 Information of missing points

每條全類型報文平均解碼時間Ta=1.12×10-4s，每條全類型報文發(fā)射間隔To=0.25 s，數(shù)據(jù)經(jīng)過處理傳入終端界面時間Tu＜0.1，從表2 可以得到最小間隔Tmin＜Ta+To+Tu均成立，說明缺失點判別模型算法正確。通過距離矩陣對缺失點賦值，便于改進(jìn)后的KNN算法檢測分類。

3.2 算法對比

在軌跡點分類的具體應(yīng)用中，選擇適用的距離計算方法是決定KNN 分類結(jié)果質(zhì)量的關(guān)鍵性因素。常見的空間距離算法有歐氏度量和曼哈頓度量，如1.2節(jié)所述。本節(jié)將歐氏度量、曼哈頓度量和標(biāo)準(zhǔn)化歐氏度量進(jìn)行實驗對比，通過對比3 種距離算法檢測ADS-B 異常軌跡點的結(jié)果，驗證3 種不同算法檢測ADS-B 異常點的準(zhǔn)確率。所有樣本訓(xùn)練集與測試集比例設(shè)置為6 ∶4，每個樣本設(shè)置不同k 值（樣本包括780310、781181），交叉驗證采用不同距離算法的KNN檢測ADS-B 異常軌跡點，并計算每種算法檢測出的異常點數(shù)量與每個樣本異常點總數(shù)之間占比，其結(jié)果如圖5 所示。

圖5 算法結(jié)果對比Fig.5 Comparison of algorithm results

從圖5 可得標(biāo)準(zhǔn)化歐氏度量的精確度均高于歐氏度量和曼哈頓度量。通過對比不同k 值各算法準(zhǔn)確率，當(dāng)k=2 時為最佳選擇。通過同一終端區(qū)的不同航班軌跡、不同k 值和不同距離算法的交叉驗證，驗證了改進(jìn)后的KNN 算法檢測ADS-B 異常點是完全可行的。

3.3 軌跡修復(fù)

在進(jìn)行插補(bǔ)修復(fù)前，需要將檢測出的異常值的經(jīng)度和緯度歸0，并賦值序號標(biāo)簽，方便插補(bǔ)算法識別需插補(bǔ)位置，以免忽略。完成賦值后，篩選出樣本中緯度和經(jīng)度將其作為修補(bǔ)序列C[X，Y]輸入多重插補(bǔ)算法中，其對比結(jié)果如圖6 所示。

圖6 實驗結(jié)果對比Fig.6 Comparison of experimental results

由圖6 可知，ADS-B 軌跡包含線性軌跡與非線性軌跡，經(jīng)過缺失點篩選和改進(jìn)后的KNN 檢測出異常點，將其清洗后通過多重插補(bǔ)能夠有效完善軌跡圖像，進(jìn)一步提高數(shù)據(jù)的完整性。

4 結(jié)語

航空器ADS-B 技術(shù)的普及應(yīng)用，使得ADS-B數(shù)據(jù)成為當(dāng)前航空器交通流分析、軌跡監(jiān)管和空域運(yùn)行態(tài)勢分析的重要數(shù)據(jù)來源之一。然而ADS-B 數(shù)據(jù)的缺失及異常點的出現(xiàn)給相關(guān)系統(tǒng)及用戶的應(yīng)用分析造成一定困難。本文方法不僅可以有效判別是否存在缺失點，而且可以有效檢測出除缺失點外的其他異常點，能夠有效插補(bǔ)航空器ADS-B 軌跡，提升軌跡圖像連續(xù)性，降低因異常點或圖形問題而造成的不安全事件發(fā)生概率。針對ADS-B 數(shù)據(jù)缺失水平問題，應(yīng)進(jìn)一步研究ADS-B 數(shù)據(jù)質(zhì)量分析指標(biāo)體系，并根據(jù)用戶需求和實際應(yīng)用場景不斷改善和發(fā)展ADS-B 技術(shù)，提高數(shù)據(jù)鏈傳輸質(zhì)量和地面站接收能力，從根本上解決數(shù)據(jù)缺失和異常點出現(xiàn)問題。