基于密度擬合的區域航道挖掘方法

聶熠文，劉軍偉

（1.中國電子科技集團公司第三十八研究所，安徽合肥 230088；2.孔徑陣列與空間探測安徽省實驗室，安徽合肥 230088）

0 引 言

雷達探測空中目標，不同于地面目標跟蹤，可以利用交通運輸網等先驗信息輔助目標跟蹤，雷達通常只能利用探測到目標的空間相對位置、運動特征等特征對目標進行跟蹤，難以獲得觀測區域內的飛行器航道信息進行輔助。實際上，由于地形、流量控制、洋流、氣象條件、燃料等因素，空中、海上等雷達目標也遵循特定的航道運動，以便快速到達目標位置，而這些航道信息蘊含在目標的運動軌跡中。因此，可以對一定區域內一定時間段內積累的目標運動軌跡進行分析處理，從中獲得隱含的區域內航道信息，作為先驗輔助信息，在氣象、電磁環境復雜的情況下幫助雷達進行目標跟蹤，提高探測效率和情報質量。

目前大部分工作是對地面道路［1-2］、海（河）面艦船［3］進行航道提取。Li 等人［4］提出改進的動態時間規整算法（DTW）實現時間序列的分類與聚類，該方法能夠較為精確的度量序列差異，但DTW時間復雜度為O（m·n），即與對比的兩條序列的采樣點數之積正相關，因此當軌跡數目多，長度大時，該算法效率較低。Guo等人利用高斯分布對道路上的GPS點進行建模，從中抽取路網信息［5］。由于GPS 精度高，道路上的運動點偏移概率小，因此采用高斯分布對GPS 點進行建模行之有效；對于雷達探測的空中目標而言，航道并不可見且探測獲得點跡的位置精度和穩定性遠低于GPS，導致高斯分布難以在大量航跡點上對其準確建模，提取航道。另一類方法將軌跡點視為圖像中的像素點，提取航道［6-7］。但在由不同采樣率獲得的航跡集合上，這種方法難以進行合適的平面網格化，獲取“像素”。對于水面目標，大部分航道提取工作基于位置精度較高的AIS 信息提取水面航道，如Lee 等人［8-12］利用AIS 建立核密度函數或密度聚類等方法，獲取航道輪廓，再提取輪廓中心線作為航道。

目前大部分提取航道的方法是基于合作式目標，可獲取其主動提供的準確定位信息，然而在雷達探測中，有時難以獲取由目標提供的準確定位，只有被動探測數據。被動探測由于系統誤差、隨機誤差、環境噪聲等原因，精度有限，即使目標是沿著某航道運動，其結果也可能存在偏差。因此我們提出一種基于密度擬合的航道挖掘方法，在無先驗知識且不引入利用航跡點時間關聯性的情況下，將航跡作為最小處理單元，實現精細化航道提取的目標。

1 問題背景與基本定義

本文我們將給出在無先驗知識的條件下，從積累的雷達數據中提取熱門航路的方法，基于密度擬合區域航道挖掘方法（Density based Fitting Method for Route Extraction，DFRE）。具體來說，積累的雷達數據由一個區域內的目標航跡點組成，包含經度、緯度、時間等信息。我們的方法目前只考慮二維平面內的航道挖掘，并未考慮空間高度；忽略航跡中各點的時間關系，將每條航跡視為一個處理單元，找到區域內的熱點航路。熱點航路由一系列斜率參數及起始點經緯坐標組成的線段集合表示。

1.1 方法框架

本文提出的方法包含航跡參數擬合及擬合參數聚類兩個部分。根據目標航跡特征可知，大部分航道呈直線狀，少量航道可呈多條直線組成的折線/曲線狀。因此本文方法首先利用最小二乘方法對航跡數據集進行擬合，獲得各航跡在二維平面上斜率和截距兩個擬合參數，再使用基于密度的空間聚類方法DBSCAN 對擬合參數進行聚類。DBSCAN方法是一種無監督聚類方法，無需提前進行模型訓練，不受聚類形狀限制（只要空間分布連續），可應用于未知聚類數目的聚類。因此對于本方法中未知擬合參數空間分布的情況，DBSCAN方法適用。

1.2 參數聚類

DBSCAN方法將聚類對象抽象為空間中的點，某點的密度由其半徑δ（δ-鄰域）內的點跡數目定義。根據點的密度，空間中的點可分為核心點、邊緣點和噪聲點。其中，點密度不低于閾值γ的為核心點，位于核心點密度范圍內但自身點密度小于γ的為邊緣點，其他則為噪聲點。DBSCAN通過循環迭代完成聚類。

所有的聚類點初始標記為未訪問。隨機選取一個未訪問的p點，若p的δ鄰域內點數不少于γ，則標記p為類別Li且已訪問。將其δ-鄰域中不帶類別標簽（噪聲不是簇標簽）的點加入Li候選集；若點數小于γ，那么p被標記為噪聲且已訪問。對于Li候選集合中的某個點p*，聚類時可能有兩種情況。

情況1:p*沒有被訪問。在這種情況下，p*被標記為已訪問且類別Li。若p*δ-鄰域中點數不小于γ，則將其鄰域中所有未被標記類別的點放入Li候選集。

情況2:p*已被訪問。p*是Li的邊緣點，標記為Li。

上述過程在所有聚類點上迭代進行，直到它們都已訪問且被標記類別（或噪聲）。

2 基于密度擬合的區域航道挖掘流程

2.1 航跡參數擬合

具體航道提取流程如算法1所示，首先對每條航跡進行航跡擬合（算法2），輸出每條航跡的參數集合H={(kh,bh,starth,endh)}Hh=1。集合中由擬合參數kh、bh及起始點位置starth、endh表示相應航跡，利用聚類算法將相近的航跡參數聚為一類（算法3），此處選擇數值較小的聚類鄰近閾值，保證聚類的精確性。之后對屬于同一類的航跡根據起始終結位置進行接續處理（算法4），保證初步提取航道的連續性，避免由于航跡點跡缺失或航跡結束導致的航道不完整。較小的聚類鄰近閾值雖能保證聚類的精確性，但易導致一些屬于同一類的航跡被分為幾類。因此在完成航道接續后，根據各條航道可能存在或交叉或鄰近的位置關系，利用算法再進行一次航道合并（算法5）和接續（算法4），得到最終的區域航道信息。

算法1:基于密度擬合區域航道挖掘（DFRE）

輸入:航跡集合T={tm|tm=(xm1,ym1),(xm2,ym2),…,(xmN,ymN)}Mm=1，擬合精度acc_threshold，劃分單元length_unit，距離閾值δk,δb，連接閾值σc，合并閾值σk,σlat,σlon輸出:航道參數集合H={(kh,bh,starth,endh)}Hh=1（starth,endh為航道起始、終結坐標）1Γ=FFT(T，acc_threshold，length_unit)；2Λ=FPC(Γ,δk,δb)；3Θ=CFT(Λ,δc)；4H=MNR(Θ,σk,σlat,σlon)；

對航跡的參數擬合過程由算法2所示，假設某目標航跡tm由點跡序列組成{(xm1,ym1),(xm2,ym2),(xm3,ym3),…,(xmN,ymN)}，其中x為經度，y為緯度，各航跡點跡數目N并不相同。由該點跡擬合獲得的直線段斜率和截距如下式所示，對于沿某經線運動的航跡，其k值無窮大（即x=c）。對于此類航跡，采用以緯度為x，經度為y進行擬合。將式中x與y互換，即可得到以沿經線運動的航跡的擬合參數。用(km,bm)參數對表示每條航跡（第2行）。

隨后判斷當前擬合是否滿足擬合精度閾值，本方法以航跡中各點的緯度與由各點經度經擬合參數獲得的擬合緯度間的絕對差值的均值作為衡量標準（第3行）。若滿足精度閾值，則該擬合參數及對應的航跡起始終結點坐標組成向量(km,bm,startl,endl)放入集合Γ（第5 行），若不滿足閾值，則將其放入分段航跡集合（第7行）。對該集合內的各條航跡進行分段，每p個航跡點作為一獨立航跡（p取值視數據情況調整，本文實驗中p取20），進行與上述過程類似的參數擬合及精度判斷，滿足閾值要求的分段航跡補充進集合Γ，反之則舍棄不用（第9～13 行）。注:本文中所指航跡起始為當前航跡中經度最小點和終結為當前航跡中經度最大點，與航跡實際產生過程無關。此外，當前算法基于k值存在的情況，若趨向于無窮大（航跡平行于經線），則將緯度視為x、經度視為y，同時將算法中涉及到的經緯度互換即可。

算法2:航跡參數擬合（FFT）

輸入:航跡集合{tm|tm=(xm1,ym1),(xm2,ym2),…,(xmN,ymN)}Mm=1，擬合精度acc_threshold，航跡劃分單元length_unit輸出:滿足精度的航跡參數集合Γ={(kl,bl,startl,endl)}L l=1（startl,endl分別為航跡起始、終結點坐標）1for m=1 to M 2按照公式（1）計算航跡tm的擬合參數(km,bm)；3計算擬合精度accm=(∑n=1 N|kmxmn+ bm-ymn|)/N；4if accm ≤acc_threshold 5(km,bm,startl,endl)放入Γ；6else 7將tm放入分段航跡集合{}t′s S s=1;8endfor 9for s=1 to S 10將航跡t′s劃分為長度為length_unit的多個子航跡；11計算各個子航跡的擬合參數并計算擬合精度；12滿足擬合精度的子航跡參數(ks,bs,starts,ends)加入集合Γ；13endfor

2.2 擬合參數聚類

在實際情況中，飛行器沿著航道運動，因此在收集到的航跡數據中可以看出在一些區域，某些位置航跡沿某個方向呈現平行聚集情況，如圖1（a）所示。因此在對航跡進行參數擬合后，聚類擬合參數，從平行聚集的航跡集中抽取相應的航道參數。具體聚類提取方法如算法3 所示。以航跡參數集合Γ={(kl,bl,startl,endl)}Ll=1、用于聚類度量的距離閾值δk,δb為輸入，對參數k與參數b分別聚類。對航跡l，將參數k與kl距離不超過δk的航跡編號分別與l構成元素(l,lδik)分別放入集合Ak，對參數kb進行同樣操作（第1～4 行），得到集合Ab。兩個集合交集即為同時滿足k與b聚類要求的航跡對（第5 行）。對集合中出現的每一個航跡編號l，將與l滿足距離關系的編號與l歸為同一類，同時以該類中的編號為參數距離衡量參考，將與其滿足距離關系的編號歸為當前類別，以此類推，當前類別不斷延伸，直至不再有新編號加入該類別，則該類別完成聚類（第7～23行）。

圖1 航跡分布及聚類參數統計

由于聚類過程是基于距離關系類別集合不斷擴充的過程，若距離閾值δk,δb過大，則會導致所有航跡編號聚為同一類；若閾值過小，則會導致聚類失敗。因此，我們隨機選取了訓練集中20%的數據，計算并統計這些航跡參數間的距離，如圖1（b）所示。為保證同一類別內的高內聚性，避免粗粒度聚類導致航道提取不準確，我們以拐點所在區域內較小的距離值作為距離閾值，達到航跡參數聚類效果。

算法3:擬合參數聚類（FPC）

輸入:航跡參數集合Γ={(kl,bl,startl,endl)}Ll=1，距離閾值δk,δb輸出:擬合航跡聚類集合Λ={(kl,bl,startl,endl,C)}L l=1 1for l=1 to L 2將擬合斜率k位于[kl-δk,kl+ δk]區間內的航跡編號lδk i(l,lδk i)I 存入集合Ak={}L i=1l=1；3將擬合截距b位于[bl-δb,bl+ δb]區間內的航跡編號lδb i(l,lδb j)I 存入集合Ab={}L j=1l=1；4endfor 5Akb=Ak ?Ab；6C=1,BC=?；7for l=1 to L 8if l已被訪問9 continue；10else 11將l放入集合BC，每個滿足(l,lδkb) ∈Akb且lδkb ?BC的航跡編號lδkb放入集合BC；12l標記為已訪問；13e指向BC中第一個元素；14while e不為?15每個滿足(e,eδkb)∈Akb且eδkb ?BC的航跡編號eδkb放入集合BC；16e標記為已訪問；17e指向BC中下一個元素；18endwhile 19BC中編號的類別標記為C；20C++,BC=?；21將BC中航跡編號對應的參數向量(kg,bg,startg,endg,C)加入集合Λ；22endif 23endfor

2.3 擬合航跡連接

由于不同航跡不同的起始位置和終結位置，歸屬同一類的航跡參數會存在位置重疊或者不連續的情況。為了去除重疊，同時保證航道連貫，我們利用擬合航跡連接算法進行航跡去重和接續處理。具體過程如算法4 所示。對于某一航跡參數類別，將該類別中的參數向量按照航跡起始點的經度依升序排序（第2行）。

令s1表示初始連接航跡，s2表示最新連接上的航跡，e表示當前待連接航跡，若e與s2存在重疊或者二者端點的經度差（e的起始點和s2的終結點）不超過連接閾值，則e連接至s2航跡，s2向當前連接航跡中結束點經度最大的航跡，e指向下一待連接航跡（第5～8，14 行）；若e與s2不滿足條件，則中斷當前連接過程，以s1航跡的起始為起始，s2航跡終結為結束，從s1到s2的航跡參數k與b的均值為參數，表示經連接產生的航道，同時下一待連接航跡作為新的初始（最新）連接航跡s1和s2表示，e則指向s1和s2所示航跡的下一條航跡（第10～12，14行）。

上述過程循環直至所有類別內的航跡參數都完成連接評估，完成航道提取。

算法4:擬合航跡連接（CFT）

輸入:擬合航跡聚類集合Λ={(kl,bl,startl,endl,C)}Ll=1，連接閾值σc輸出:航道（跡）鏈集合Θ={(kg,bg,startg,endg)}G g=1 1for each C 2根據start_lon值對同屬C類別的航跡參數向量排序；3s1=1，s2=s1，e=s1+1；4while e ≤||C 5 if starte_lon ≤ends2_lon+ σc //若e航跡開始經度不超過s2航跡結束經度的σc擴展6 if ende_lon ＞ends2_lon //若e航跡結束經度大于s2航跡結束經度7 s2=e；8 endif 9 else 10以s1至e所有航跡的k，b均值為合并航跡k，b均值，以s1航跡的起始為起始，11s2航跡的結束為結束，存入集合Θ；12s1=e，s2=s1；13endif 14e=e+1；15endwhile 16以s1至e所有航跡的k，b均值為連接航跡的k，b均值，以s1航跡的起始為起始，e航跡的結束為結束，存入集合Θ；17endfor

2.4 鄰近航跡合并

雖然在航跡參數聚類后進行航跡連接，初步提取出了連貫的航道，但聚類中使用的距離閾值較小，導致有些鄰近航跡可能并未被歸屬同類。因此我們在初步提取的航道上再進行鄰近航道合并，去除重疊航道。具體過程由算法5所示。對某航道t，根據其參數kt，在集合Θ中找到滿足閾值σk的航道編號，放入同一分組（第2行）。分別計算該分組內的航道i與t的距離，若二者存在交點或二者經度交疊部分的緯度差值不超過閾值σlat（若二者斜率參數絕對值之和大于2，則為緯度交疊部分的經度差值不超過閾值σlon），則記t與i為一組鄰近航道對，放入集合Q（第3～13 行）。對鄰近航道對集合Q，利用FPC算法在航道參數做一次聚類和航跡連接（第15、16行），避免由于鄰近關系引入的參數誤差，保證航道方向、位置的準確性和航道連貫性。

通過本方法提取出的航道可用于后續該區域內的目標跟蹤。當環境復雜，雜波較多時，可將跟蹤航跡與航道匹配，利用航道幫助過濾異常點，防止目標跟偏，航跡斷批。

算法5:鄰近航道合并（MNR）

輸入:航跡（道）集合Θ={(kl,bl,startl,endl)}Ll=1，合并閾值σk,σlat,σlon輸出:航道參數集合H={(kh,bh,starth,endh)}H h=1 1for t=1 to L 2將ki ∈[kt,kt+ σk]的航道編號i放入集合I；3for each i ∈I 4 if航道i與航道t存在交點5(t,i)放入集合Q；6 elseif 航道i與t的經度存在重疊部分且該部分最大距離（緯度差值）不超過σlat 7(t,i)放入集合Q；8 elseif||ki+||kt ＞2 9 if航道i與t的緯度存在重疊部分且該部分最大距離（經度差值）不超過σlon 10(t,i)放入集合Q；11endif 12endif 13endfor 14endfor 15基于集合Q運行FPC算法第6～23行，獲得集合?；16H=CFT(?,δc)；

3 實驗結果及分析

3.1 實驗設置

本文在實測航跡數據集上進行算法性能測試。航跡數據集包含同一區域不同時間的三次采集結果，共計包含2 501 條航跡。隨機抽取其中1 600條作為訓練集，剩余數據作為測試集，測試提取獲得的航道是否能符合區域內航跡情況。訓練集和測試集中的航跡分布如圖2 所示。測試集和訓練集中的航道分布基本一致，因此可用測試集對從訓練集中抽取出的結果進行航道匹配性測試。

圖2 航跡數據集

3.2 實驗結果分析

我們將基于區域熱點的航道骨架提取算法作為控制對比組，與本文方法做實驗對比。圖3（a）為控制組方法從訓練集中提取的航道，圖3（b）為本方法從訓練集中提取的航道與測試集航跡的對比圖，其中藍色點線為航道，綠色點為測試集航跡。觀察可知，對比組方法雖然航道連通性高，基本勾畫出航跡分布輪廓，但航道走向勾畫精準度不高；本方法航道間連通性較低，但對整體航道走向勾畫較為精準。為衡量提取航道與測試航跡的匹配性，我們以航跡與航道是否處于同一區域、航跡點與航道間的距離，以及航跡匹配上的航道點的整體航向衡量航跡航向的匹配度。

圖3 航道提取對比

具體來說，在航跡與航道存在經度或緯度重疊的部分，統計航跡點與航道間距離不超過設定閾值、航向匹配的航跡點數目，計算匹配的航跡點數目與該航跡總航跡點數目的占比。由于篇幅有限，表1 僅展示了兩種不同距離閾值，不同匹配比例下的兩種航道對測試集航跡的匹配性。從表中可知，在控制變量的條件下，匹配到的航跡數目和距離閾值成正相關。與對比組結果相比，本方法所獲得航道的匹配度普遍大幅度優于對比組方法，要求的匹配度越高，本方法的優勢越明顯。雖然對比組方法將勾畫出連通的整體航道，但受限于航道骨架提取需將區域網格化的要求，難以對航道方向進行精確提取，同時易受少量異常值影響，導致航道精確度差，匹配度較低；而本方法利用了參數擬合的思想，避免了區域網格化引入的不必要誤差，保證了航道的精確度。

表1 航道-航跡匹配統計

圖4 展示了在距離閾值0.1°時，匹配上航道的航跡結果，圖中紅色點跡即為匹配上的航跡點。圖中黑色虛線圓圈標識了航道匹配較差的區域。對比提取航道的訓練航跡集，圓圈1內的訓練集本身的航跡較為稀疏，是否存在航道難確定，為避免提取有誤，本算法忽略相應航跡；圓圈2 內的區域較小，訓練集內航跡雜亂，未能形成明確航道，因此測試集中出現于此處的航跡未能有航道匹配。對比方法中雖然航跡也有匹配，但由于航道不平滑，整體匹配精度較低，此外圓圈3 內的航跡未被匹配，說明對比方法更為偏重對航跡存在區域的輪廓描述，輪廓內航道的精細化提取較弱。

圖4 測試集航跡與提取出的航道對比

通過上述對比分析，可以看出本方法具有較強的精細化航道提取能力，在沒有先驗知識的情況下，能夠從積累的歷史航跡信息抽象區域內目標運動軌跡，精確有效地提取區域內航道，為航跡跟蹤、目標推斷提供信息支撐。

4 結束語

本文研究了在不使用額外輔助信息的條件下，如何從航跡數據集中精細化提取航道，提出了一種基于密度擬合的區域航道挖掘方法。由于燃料、成本、飛行安全等限制，大部分航跡近似直線，因此該方法采用直線擬合的方法獲取航跡表征參數（對于少數非直線航跡進行分段擬合），基于航跡參數實現鄰近航跡的聚類，在同類航跡中進行航道提取，最終根據航道始末位置進行航道連通，完成整個區域的航道挖掘。通過分析和對比試驗，本方法效率高且提取出航道具有較高的精細度，能夠在無先驗知識的條件下有效地挖掘區域內熱點航道信息。