聶熠文，劉軍偉

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

0 引 言

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

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

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

1 問題背景與基本定義

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

1.1 方法框架

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

1.2 參數(shù)聚類

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

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

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

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

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

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

2.1 航跡參數(shù)擬合

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

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

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

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

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

算法2:航跡參數(shù)擬合（FFT）

輸入:航跡集合{tm|tm=(xm1,ym1),(xm2,ym2),…,(xmN,ymN)}Mm=1，擬合精度acc_threshold，航跡劃分單元length_unit輸出:滿足精度的航跡參數(shù)集合Γ={(kl,bl,startl,endl)}L l=1（startl,endl分別為航跡起始、終結(jié)點坐標(biāo)）1for m=1 to M 2按照公式（1）計算航跡tm的擬合參數(shù)(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計算各個子航跡的擬合參數(shù)并計算擬合精度；12滿足擬合精度的子航跡參數(shù)(ks,bs,starts,ends)加入集合Γ；13endfor

2.2 擬合參數(shù)聚類

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

圖1 航跡分布及聚類參數(shù)統(tǒng)計

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

算法3:擬合參數(shù)聚類（FPC）

輸入:航跡參數(shù)集合Γ={(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]區(qū)間內(nèi)的航跡編號lδk i(l,lδk i)I 存入集合Ak={}L i=1l=1；3將擬合截距b位于[bl-δb,bl+ δb]區(qū)間內(nèi)的航跡編號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標(biāo)記為已訪問；13e指向BC中第一個元素；14while e不為?15每個滿足(e,eδkb)∈Akb且eδkb ?BC的航跡編號eδkb放入集合BC；16e標(biāo)記為已訪問；17e指向BC中下一個元素；18endwhile 19BC中編號的類別標(biāo)記為C；20C++,BC=?；21將BC中航跡編號對應(yīng)的參數(shù)向量(kg,bg,startg,endg,C)加入集合Λ；22endif 23endfor

2.3 擬合航跡連接

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

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

上述過程循環(huán)直至所有類別內(nèi)的航跡參數(shù)都完成連接評估，完成航道提取。

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

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

2.4 鄰近航跡合并

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

通過本方法提取出的航道可用于后續(xù)該區(qū)域內(nèi)的目標(biāo)跟蹤。當(dāng)環(huán)境復(fù)雜，雜波較多時，可將跟蹤航跡與航道匹配，利用航道幫助過濾異常點，防止目標(biāo)跟偏，航跡斷批。

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

輸入:航跡（道）集合Θ={(kl,bl,startl,endl)}Ll=1，合并閾值σk,σlat,σlon輸出:航道參數(shù)集合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的經(jīng)度存在重疊部分且該部分最大距離（緯度差值）不超過σlat 7(t,i)放入集合Q；8 elseif||ki+||kt ＞2 9 if航道i與t的緯度存在重疊部分且該部分最大距離（經(jīng)度差值）不超過σlon 10(t,i)放入集合Q；11endif 12endif 13endfor 14endfor 15基于集合Q運行FPC算法第6～23行，獲得集合?；16H=CFT(?,δc)；

3 實驗結(jié)果及分析

3.1 實驗設(shè)置

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

圖2 航跡數(shù)據(jù)集

3.2 實驗結(jié)果分析

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

圖3 航道提取對比

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

表1 航道-航跡匹配統(tǒng)計

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

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

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

4 結(jié)束語

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