黃窈蕙，范文濤，劉柳楊

(交通運輸部科學研究院， 北京 100029)

道路交通與經(jīng)濟發(fā)展息息相關。當前國內(nèi)經(jīng)濟發(fā)展迅捷，城市化建設速度大大加快，道路網(wǎng)絡更新迅速，現(xiàn)勢性強的道路網(wǎng)數(shù)據(jù)能為經(jīng)濟發(fā)展提供助力，同時也是路徑規(guī)劃、城市建設及位置服務的基石，更為數(shù)字城市、智慧城市的發(fā)展奠定了基礎。傳統(tǒng)的通過測繪人員的道路采集方法，成本高且現(xiàn)勢性低。通過衛(wèi)星遙感及無人機遙感手段獲取遙感影像，然后分析遙感影像從中獲取道路網(wǎng)信息也是一個發(fā)展熱點，但其成本相對較高，且道路網(wǎng)容易被其他地物所遮掩，從遙感影像中獲取準確的道路網(wǎng)數(shù)據(jù)還存在諸多技術難點。近年來，隨著GPS的普及，眾源軌跡數(shù)據(jù)量成爆發(fā)式增長，其數(shù)據(jù)成本相對較低且時效性強，從眾源軌跡數(shù)據(jù)中提取新增道路是一種很有前景的方法。

基于GPS軌跡的交通道路提取算法，從原理上主要可分為聚類、軌跡合成、密度估算等[1-2]。聚類方法不考慮軌跡點間的連接關系，而是基于大規(guī)模軌跡點，利用諸如K-means等算法，先進行軌跡點的空間聚類分析，然后從聚類中心線中提取道路要素[3-4]。該方法受聚類算法的輸入?yún)?shù)影響較大，且軌跡點的空間分布具有空間不均勻性，難以找到適合全局的輸入?yún)?shù)。軌跡合成算法考慮了單個軌跡點數(shù)據(jù)，以及軌跡點所連成的軌跡線，通過合成后者生成道路[5-6]。該方法在GPS軌跡數(shù)據(jù)精確度較高的情況下效果較好，但在數(shù)據(jù)精確度較低時容易生成冗余道路[7]。核密度估算法是另一類從GPS軌跡點提取道路信息的方法。該方法首先將GPS軌跡點轉(zhuǎn)換成柵格圖像，圖像上每一個像素的像素值代表了該像素上經(jīng)過軌跡點的密度，在給定的特定閾值下，對柵格圖像進行篩選將圖像轉(zhuǎn)換為二值柵格圖像，提取柵格骨架線作為道路中心線[8-9]。從軌跡點中提取道路網(wǎng)的另一難點是，提取的道路需要融合到現(xiàn)有的道路網(wǎng)中，而基于軌跡點提取的道路網(wǎng)存在“毛刺”、延伸不到位或過于延伸的問題，因而也存在如何去“毛刺”、如何有效融合現(xiàn)有道路網(wǎng)的技術難點[10-12]。

本文提出一種道路網(wǎng)自動更新和道路融合方法。該方法基于歷史道路網(wǎng)和GPS軌跡數(shù)據(jù)集，通過失配軌跡提取、軌跡融合、新增道路生成3個步驟，完成道路網(wǎng)更新。首先利用最短距離分析和隱馬爾科夫模型分析失配軌跡段，隨后采用基于Delaunay三角網(wǎng)和基于山脊線的兩種軌跡融合方法，對失配軌跡進行融合，得到新增道路骨架線，最后對得到的骨架線進行簡化、延長和節(jié)點匹配操作，得到新增道路，更新歷史道路網(wǎng)。

1 失配軌跡提取

從海量的軌跡數(shù)據(jù)中提取出矢配軌跡是利用GPS軌跡增量更新道路網(wǎng)算法的核心。本文提出的失配軌跡提取方法，結合了幾何最短距離分析和隱馬爾可夫模型分析，在完成失配軌跡提取前，首先給出相關的基本概念。

定義1：道路網(wǎng)眼(road network mesh，RNM)，即道路網(wǎng)中由若干條道路圍成的最小面狀封閉區(qū)域。一個城市可以看成具有相鄰(拓撲)關系的一系列RNM的空間劃分。

定義2：軌跡線(trajectory)，一條軌跡線T由存在偏序關系的一系列連續(xù)軌跡點(p1，p2，…，pn)構成，其中軌跡點pi=(si，ti)由空間位置si和時間ti描述，TM1≤t1<…

1.1 幾何最短距離分析

幾何最短距離分析通過預設兩個距離閾值對道路網(wǎng)眼內(nèi)側進行雙緩沖實現(xiàn)，即最遠候選點距r和最鄰近距γ(r>γ)。計算軌跡點與道路網(wǎng)中道路的幾何最短距離(d)，對給定一個距離閾值r，若d>r，則該軌跡點為失配軌跡點。

1.2 隱馬爾可夫模型分析

可以應用隱馬爾可夫模型(hidden Markov model，HMM)分析GPS軌跡點的失配性[2]。在軌跡點的失配判別中，由于軌跡點p的觀測坐標與其真實值存在誤差，可假設該誤差服從正態(tài)分布，假設該軌跡點隸屬于歷史道路網(wǎng)(非新道路)某道路上的cp點，將存在一定的表現(xiàn)概率，其中cp即是HMM的可見狀態(tài)，根據(jù)誤差正態(tài)分布定義，p對應cp的表現(xiàn)概率可表示為

(1)

假設pi和pi+1為軌跡線T上相鄰的兩個軌跡點，這里用下式表示其間的狀態(tài)轉(zhuǎn)換概率Pr。

(2)

2 軌跡融合

將道路網(wǎng)眼內(nèi)的失配軌跡進行融合，從而得到新增道路的骨架線，這是生成新增道路的基礎。本文提出了基于Delaunay三角網(wǎng)的軌跡融合方法和基于山脊線的新增道路骨架線提取方法。在基于Delaunay三角網(wǎng)的軌跡融合方法中，因為約束Delaunay三角剖分算法的時間復雜度為O(nlogn)[13]，融合過程中隨著融合軌跡段幾何復雜度的增加，每次融合需要的時間也相應增加。在基于山脊線的軌跡融合方法中，在數(shù)據(jù)預處理階段需要生成遍布整個網(wǎng)眼的格網(wǎng)，且山脊線提取需要一定的處理時間。因此在道路網(wǎng)眼內(nèi)失配軌跡數(shù)量較少時，使用基于Delaunay三角網(wǎng)的軌跡融合方法效率較高，避免了基于山脊線方法初期處理的操作；當失配軌跡較多時，基于山脊線的方法避免了處理時間過多的增加，從而效率更高。

2.1 基于Delaunay三角網(wǎng)的軌跡融合方法

基于Delaunay三角網(wǎng)的軌跡融合，其基本原理包括：①選取兩條相交或最短距離小于閾值的軌跡段，將其中一條TS1生成半徑為BR(buffer radius)的緩沖區(qū)；②判斷另外一條軌跡段TS2與緩沖區(qū)的空間關系，將落在緩沖區(qū)內(nèi)的軌跡段與TS1相對應的部分提取出來，建立具有權值的Delaunay三角網(wǎng)；③隨后提取三角網(wǎng)的中線，與余下的軌跡段拼接，完成軌跡融合。如圖1所示，在Road1、Road2、Road3和Road4組成的道路網(wǎng)中，有TS1和TS2兩條軌跡段經(jīng)過。首先選定TS1生成半徑為BR的緩沖區(qū)，然后依次判斷TS2中點是否落在緩沖區(qū)內(nèi)。如圖1中虛線矩形框所示，TS2上共有兩個子段落在緩沖區(qū)內(nèi)，這兩段軌跡段子段即為TS1的相似軌跡段。

圖1 相似軌跡段融合示意圖

提取到相似軌跡段后需要將其融合[14]，在得到融合后的軌跡段后，需要處理其與現(xiàn)有軌跡段的邊界處公共軌跡點。在融合軌跡段的兩端，用融合后軌跡段的端點取代原來軌跡段的端點。圖1中虛線圓形內(nèi)即為軌跡段邊界處理示意圖，其中黑色帶箭頭的折線即為處理好邊界公共軌跡點的融合后軌跡段。

2.2 基于山脊線的軌跡融合方法

基于山脊線提取算法的主要步驟為：對每條失配軌跡段進行緩沖區(qū)生成及空間疊加分析，生成失配段緩沖區(qū)的空間密度分布圖，通過提取山脊線即緩沖區(qū)密度骨架線，即可獲得新增道路。設定一個最小密度閾值，去除由于采樣誤差造成的超低密度區(qū)；通過水文(Hydrology)分析中的提取山脊線的計算方法，可獲取密度最高處對應的山脊線，即對應為新增道路。詳細步驟如下：

(1) 在道路網(wǎng)眼內(nèi)生成覆蓋整個網(wǎng)眼的格網(wǎng)，為網(wǎng)眼中每個格子賦權值為0?？紤]試驗中的數(shù)據(jù)質(zhì)量，格子的邊長可以設置在3～10 m內(nèi)。隨后根據(jù)GPS的采集誤差，為每條失配軌跡段生成半徑為10 m的緩沖區(qū)。判斷緩沖區(qū)與格網(wǎng)中格子的關系，若格子在緩沖區(qū)內(nèi)，則格子的權值+1。將所有的失配軌跡線處理完后，即生成了失配段緩沖區(qū)的空間密度分布圖。

(2) 在提取山脊線(緩沖區(qū)的密度骨架線)之前，設置一個最小的密度閾值，格網(wǎng)中格子小于閾值則將其閾值賦值為0，這樣可以去除GPS采樣誤差導致的超低密度區(qū)。此時的空間密度分布圖在數(shù)據(jù)形式上表現(xiàn)為DEM。

(3) 空間密度分布圖在空間上直觀的表現(xiàn)為連綿的山脈，圖中位置越高，意味著密度越大，即意味著新增道路的可能性越大。此時通過提取山脈的山脊線，即提取緩沖區(qū)的密度骨架線，可以得到新增道路。本文采用水文分析來提取山脊線，水文分析的具體步驟為：首先在整個格網(wǎng)區(qū)域求出最高值H，依據(jù)H對地形取反地形DEM，計算方式為(H-DEM)，然后對得到的反地形DEM進行填洼操作，將格網(wǎng)中封閉的洼地區(qū)域填平，最后進行坡度計算得到山脊線，即新增道路。

3 道路網(wǎng)更新

本文提出的基于軌跡數(shù)據(jù)的道路網(wǎng)更新方法主要包括數(shù)據(jù)預處理(包括失配軌跡提取)、軌跡融合、新道路生成3個過程。算法的輸入為歷史道路網(wǎng)及最新的GPS軌跡數(shù)據(jù)，輸出為更新后的道路網(wǎng)。

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

首先對軌跡數(shù)據(jù)進行質(zhì)量控制，依據(jù)相鄰的軌跡點的時序差t及軌跡運行速度，進行幾何距離的合理性分析，主要包括軌跡線去冗余和軌跡線分割。對于軌跡線T=(p1，p2，…，pf，pf+1，…，pn)，給定兩個距離閾值disa、disb，進行如下兩方面的質(zhì)量處理：①軌跡線去冗余，若任意相鄰兩個軌跡點間距dis(pf，pf+1)≤disa，則移除軌跡點pf；②軌跡線分割，若兩個相鄰軌跡點間距dis(pf，pf+1)>disb，則將該軌跡T從pf處截斷，生成兩個新的軌跡T1=(p1，p2，…，pf),T2=(pf+1，pf+2，…，pn)，其中，軌跡線去冗余將間距過于密集的軌跡點(如在同一位置重復采集的GPS軌跡點)刪除，減少了后續(xù)計算量。

隨后根據(jù)現(xiàn)有道路網(wǎng)和處理后的軌跡數(shù)據(jù)生成道路網(wǎng)眼，最后使用幾何最短距離分析和隱馬爾科夫模型分析得到失配軌跡段。

3.2 軌跡道路融合

在得到失配軌跡段后，判斷道路網(wǎng)眼中的失配軌跡段數(shù)量。若數(shù)量小于閾值，則采用基于Delaunay三角網(wǎng)的軌跡融合方法；相反，若數(shù)量大于閾值，則使用基于山脊線的軌跡融合方法。選取耗時最小的提取方法完成新增道路骨架線的提取。

3.3 新增道路生成

算法過程中，由于利用最鄰近距離法排除了道路網(wǎng)眼內(nèi)側所有軌跡點，造成提取的新道路存在懸掛端點而無法結合到道路網(wǎng)中，本文采用延伸端線及節(jié)點匹配的方法，首先將提取的道路進行雙向延長，與道路網(wǎng)產(chǎn)生相交節(jié)點，若相交節(jié)點與道路其他節(jié)點相距較近，則自動將相交節(jié)點移動到該鄰近節(jié)點，否則作為新節(jié)點保留。

4 試驗與分析

4.1 數(shù)據(jù)來源與處理

試驗中軌跡數(shù)據(jù)來自歐洲機器學習和知識發(fā)現(xiàn)聯(lián)合會[15]，道路網(wǎng)數(shù)據(jù)來自OpenStreetMap數(shù)據(jù)，遙感影像來自谷歌影像。

首先，對車輛軌跡數(shù)據(jù)進行質(zhì)量控制，設定閾值disa=15 m，disb=500 m，計算每條軌跡T上相鄰兩個軌跡點pi、pi+1之間的幾何距離dis(pi，pi+1)。若dis(pi，pi+1)≤15 m，則認為出租車沒有移動，刪除pi；若dis(pi，pi+1)>500 m，則認為此處GPS采樣出錯而導致采樣點丟失，將軌跡T從pi、pi+1處斷開，生成兩條新的軌跡。其次，將道路網(wǎng)數(shù)據(jù)和GPS軌跡數(shù)據(jù)生成覆蓋全市的道路網(wǎng)眼。最后，判斷失配軌跡段的數(shù)量分布和網(wǎng)眼的平均大小，當?shù)缆肪W(wǎng)眼內(nèi)失配軌跡段超過8條時，采用基于山脊線的軌跡融合方法，當?shù)缆肪W(wǎng)眼內(nèi)的失配軌跡段小于或等于8條時，采用基于Delaunay三角網(wǎng)的軌跡融合算法來提取新增道路骨架線。為了驗證本文算法在重建道路網(wǎng)的準確性，通過在歷史道路網(wǎng)中隨機刪除一定數(shù)量(30條)的路段，分析通過軌跡線計算而能夠被正確重建的路段比例。試驗在i5-4460的Winodws 10操作系統(tǒng)上完成，圖形顯示用ArcGIS 10.3，數(shù)據(jù)處理和計算利用數(shù)據(jù)庫PgSQL(PostgreSQL9.6)及其空間數(shù)據(jù)操作PostGIS(v=2.3)插件，將軌跡數(shù)據(jù)和OSM數(shù)據(jù)導入到PgSQL數(shù)據(jù)庫。

4.2 結果與討論

圖2顯示了一個典型的以新建區(qū)構成的道路網(wǎng)眼，以此為控制單元展示新建道路的提取過程及結果。該過程的關鍵是提取失配軌跡點和失配軌跡段，通過道路網(wǎng)眼γ內(nèi)側緩沖，排除了網(wǎng)眼內(nèi)共約1/2的非失配軌跡點；同時，將道路網(wǎng)眼進行r內(nèi)側緩沖，通過幾何距離分析，為每個軌跡段尋找r內(nèi)側緩沖區(qū)域外的一個最遠軌跡點(三角點所示)，該步驟過濾出網(wǎng)眼內(nèi)剩余的約90%的失配點，并還原軌跡點所在的軌跡段，獲得所有失配點所在的完整軌跡段；最后通過隱馬爾可夫模型，提取余下的失配軌跡點。通過提取的失配點，還原所有失配軌跡段，則可以獲得的失配軌跡段分布(如圖2(a)所示)，進一步利用ArcGIS的緩沖區(qū)分析和水力模型計算工具，提取出失配軌跡段的骨架線(圖2(b))，利用ArcGIS中ArcScan矢量化工具并利用道格拉斯-普克算法進行壓縮，采用延伸端線及節(jié)點匹配法，為懸掛端點尋找匹配節(jié)點并將新增路段連接到歷史道路網(wǎng)，得到更新道路網(wǎng)(圖2(c))。通過與后期的谷歌地圖的比較發(fā)現(xiàn)，更新后的道路網(wǎng)較好地體現(xiàn)了道路的更新情況(如圖2(d)所示)。此外，在歷史道路網(wǎng)中隨機扣除30條道路段，基于本文方法用GPS軌跡數(shù)據(jù)集進行重建，扣除道路網(wǎng)能夠被正確重建的比例，試驗結果顯示，其平均重建率達87%，而未能被正確重建的，主要是由GPS軌跡點數(shù)據(jù)質(zhì)量(如密度過低)引起。

5 結 論

本文提出了一種基于軌跡數(shù)據(jù)的道路網(wǎng)增量更新算法。該算法利用歷史道路網(wǎng)和現(xiàn)勢GPS軌跡數(shù)據(jù)，高效地實現(xiàn)了道路網(wǎng)的更新，包括失配軌跡提取、軌跡融合和新道路生成3個過程。利用幾何最短距離分析和隱馬爾科夫模型，提取失配軌跡，隨后根據(jù)道路網(wǎng)眼中失配軌跡的數(shù)量，選擇基于Delaunay三角網(wǎng)的軌跡融合算法或基于山脊線的軌跡融合算法，對失配軌跡段進行融合得到新增道路的骨架線，最后對得到的骨架線進行簡化、延長和節(jié)點匹配，完成新道路獲取并以此更新歷史道路網(wǎng)。相比其他用軌跡數(shù)據(jù)更新道路網(wǎng)的方法，本文提出的方法至少具有以下幾個優(yōu)勢：

(1) 使用最短距離和隱馬爾可夫模型來提取失配軌跡，相比現(xiàn)有方法(如文獻[2])，可大量減少軌跡候選點的數(shù)量，降低計算量。

(2) 采用Delaunay三角網(wǎng)和山脊線提取相結合的方式進行失配軌跡融合，可以在保證新增骨架線準確度的同時，進一步提高處理效率。

(3) 由于本算法屬于增量方式更新道路網(wǎng)，從而可保持對歷史道路版本化管理，保證歷史道路數(shù)據(jù)的連續(xù)性和一致性。

本文的主要貢獻在于提出了基于Delaunay三角網(wǎng)和基于山脊線提取相結合的軌跡融合方法，實現(xiàn)了一種高效的新增道路骨架線的提取方法。然而在兩種方法選取的閾值確定上，還需要更多的試驗進行驗證。同時，提取的新增道路骨架線還存在毛刺現(xiàn)象，還需要在后續(xù)研究中進一步改善。