基于高分辨率影像的道路中心線提取技術(shù)研究

周紹光，向 晶，邱 偉，孫金彥，凡 莉

(河海大學(xué)地球科學(xué)與工程學(xué)院測繪科學(xué)與工程系，南京 210098)

0 引言

道路的檢測與提取為地理信息系統(tǒng)(GIS)數(shù)據(jù)庫的更新以及 GIS應(yīng)用［1－2］提供了基礎(chǔ)數(shù)據(jù)，具有十分重要的現(xiàn)實意義。針對遙感影像道路提取問題，國內(nèi)外學(xué)者提出了很多算法和模型，朱長青等［3］提出了基于形態(tài)分割技術(shù)的道路網(wǎng)絡(luò)提取方法;李卉［4］結(jié)合LiDAR數(shù)據(jù)和遙感影像信息提取城市道路并進行建模;Zhang等［5］使用拉東變換(radon transform)探測道路中心線和估計道路寬度;周紹光等［6］結(jié)合圖割(graph cuts)理論和形狀先驗的思想，從高分辨率遙感影像中檢測道路段。此外，曾發(fā)明等［7］、陳卓等［8］、蔡紅玥等［9］等也給出了很好的道路提取策略。

部分研究著力于從獲取的道路條帶中提取道路中心線［10－11］，其中細(xì)化算法［12］是最普遍的，但該算法會產(chǎn)生很多“毛刺”，降低了道路中心線的光滑性和準(zhǔn)確性。近期，一種叫做子空間約束的均值飄移算法［13］被應(yīng)用到由道路條帶獲取精確中心線的研究中，該方法可以去除噪聲點，獲取的道路中心線比較光滑、沒有“毛刺”;但該方法是有偏估計［14］，嚴(yán)重影響了道路中心線提取結(jié)果的精確性。

本文提出了一種基于道路條帶自動檢測道路中心線的方法。首先基于概率增進樹算法獲得道路候選點，并結(jié)合形態(tài)學(xué)運算獲得道路條帶;然后對獲取的道路條帶進行細(xì)化、獲得道路骨架數(shù)組，通過迭代提取最長測地線，獲得初始道路中心線;再根據(jù)Dijkstra最短路徑算法優(yōu)化每條道路中心線;最后斷開道路中心線中角度變化大的部位，并根據(jù)方向一致性和道路連續(xù)性重新連接道路中心線段，獲得最終結(jié)果。

1 道路中心線提取算法

圖1示出道路中心線提取的技術(shù)流程。

圖1 道路中心線提取技術(shù)流程Fig.1 Technical flowchart of road centerline extraction

1.1 道路條帶的獲取

利用概率增進樹算法［15］從遙感影像中提取感興趣目標(biāo)的輪廓線是一種獨具特色的方法，相關(guān)的研究成果已證明了這一策略的效果。本文借鑒該思想，探索了從遙感影像中提取道路的方法。該方法分為學(xué)習(xí)和測試2個階段:①在學(xué)習(xí)階段，選擇并組合大量不同的特征，將概率增進樹與決策樹相結(jié)合，通過訓(xùn)練得到分類模型;②在測試階段，利用訓(xùn)練好的模型計算出影像中各點屬于道路的概率，即P(S/IN(c))，其中:IN(c)為實驗影像I中以c為圓心、具有一定半徑的圖斑(本文中半徑取25個像素的長度);S為判斷道路點函數(shù)，若S=1則屬于道路點，S=0則為非道路點。通過以上所建模型對測試影像進行實驗，獲得了道路候選點。

概率增進樹算法幾乎能檢測出所有的道路點。但鑒于概率點不能作為分割結(jié)果，并且不可避免地會有些灰度、紋理與道路相似的房屋點被誤判為道路點，而且路面上存在“孔洞”;故本文結(jié)合路徑形態(tài)學(xué)［16－17］以及數(shù)學(xué)形態(tài)學(xué)中的一些其他運算對道路候選點進行處理，獲得了令人滿意的道路條帶。

1.2 初始道路中心線的獲取

細(xì)化算法是一種從道路條帶中提取中心線的有效方法，但提取結(jié)果會伴隨著許多“毛刺”，并且局部方向不一致。本文在保證道路網(wǎng)拓?fù)潢P(guān)系正確的前提下，采取求各對端點間測地距離的策略來消除“毛刺”。初始道路中心線獲取的流程如圖2所示。

圖2 初始道路中心線獲取流程Fig.2 Flow chart of initial road centerline acquisition

圖2中提到的閾值是根據(jù)骨架圖中最短道路中心線的長度確定的，如果閾值過大，大于該閾值的“毛刺”將被錯分為道路中心線。經(jīng)過該流程，起始道路中心線被存入一個新的數(shù)組，而骨架圖中只剩下“毛刺”。

1.3 道路中心線的優(yōu)化

上文去除了骨架圖中的“毛刺”并獲得了起始道路中心線，但局部道路中心線不夠平滑。本文期望獲得的道路中心線走向與參考方向一致，位置必須盡量接近道路條帶中心。為此，引入能量最?。?8］框架來優(yōu)化道路中心線。

1.3.1 直線段獲取

道路大多數(shù)呈直線狀，或者稍微彎曲，對于后者也可將其視為由若干連續(xù)的直線段組成。這里要用到直線匹配法［19］獲取一條初始道路中心線上所有像素的局部方向。直線匹配法的核心思想是將矩形道路視為由一系列的匹配直線構(gòu)成，提取道路段就是檢測出具有一定方向和長度的匹配直線。該方法最終獲得每個像素的距離和方向矩陣，此處只用到了方向數(shù)據(jù)。直線段獲取的具體步驟為:

1)用一個5像素×5像素的移動窗口，對所有位置的值為1的結(jié)構(gòu)元素膨脹起始道路中心線(直線匹配法不適用于單像素寬曲線);

2)利用步驟1)得到的數(shù)組獲取方向數(shù)據(jù)，取出起始道路中心線上每個點的局部方向數(shù)據(jù);

3)將具有相同局部方向(差異在5°范圍內(nèi))的連續(xù)點歸為同一條直線段;

4)一條起始道路中心線被分為若干連續(xù)的直線段，設(shè)鄰近直線段的交點為關(guān)鍵點(見圖3中紫色的點)。值得指出的是，一條道路中心線的2個端點必須被設(shè)為關(guān)鍵點;

圖3 沿道路中心線的關(guān)鍵點和道路截面Fig.3 Crucial points and road sections determ ined along road centerlines

5)對中心線上相鄰關(guān)鍵點間的像素進行直線擬合，并將擬合直線的方向作為相應(yīng)中心線段的參考方向。

1.3.2 道路截面確定

沿著關(guān)鍵點處曲線的法線方向取截面(如圖3中的藍色線段)，截面的長度需包含最優(yōu)中心點。為便于理解，本文用Vi={v1，…，vj，…，vm}代表第i個截面的點集，一個截面中的所有點都分別與其鄰近截面中的任意點相連，這些連接的短線段組成了所有的候選線。優(yōu)化的過程實際上是從每個截面點集中獲取最優(yōu)點，并順次連接成線。該思想充分利用了道路的形狀特征，降低了初始中心線上局部方向波動的影響。

1.3.3 路徑代價設(shè)計

考慮到中心線的走向應(yīng)與參考方向一致、位置應(yīng)盡量接近道路條帶的中心這2條優(yōu)化規(guī)則，路徑代價將受角度約束和影像約束這2個因素制約，故將能量函數(shù)定義為

式中:Eangle為測定候選線方向與相應(yīng)參考方向的差異;Eimage為影像約束能量;α為用于權(quán)衡2個約束能量相對重要性的參數(shù)。如果α值很大，相鄰兩截面間最優(yōu)候選線的方向和相應(yīng)的參考方向可能存在很大的差異;如果α值過小，則優(yōu)化的道路中心線可能會偏離條帶中心很多。

圖4示出角度約束的含義。

圖4 角度能量的含義Fig.4 M eaning of angle energy

圖4中vj和vl分別為截面點集Vi和Vi+1中的任意一點，紅色實線代表第i個截面和第i+1個截面之間候選線的參考方向 φi，i+1，φvjvl代表候選線段vj vl的方向。則第一個約束條件可定義為

角度差異越大，代價越大。

除了角度約束，影像約束也至關(guān)重要，它決定道路中心線能否盡量位于道路條帶的中心。本文首先計算影像中任一像素最鄰近非零點到該點的歐式距離，并將該距離值作為該像素的灰度值;則道路條帶部分的點像素值大于0，并且越靠近條帶中心、點像素的值越大;非道路條帶部分點的像素值都為0，因為其到自身的距離最近。為計算方便，對距離轉(zhuǎn)換結(jié)果取反，作為最終的參考影像(圖5)。

圖5 影像約束示意圖Fig.5 Sketch map of image constraints

假設(shè)圖4中候選線段vj vl由m個點組成，則局部影像代價定義為

式中I為參考影像。影像灰度值越大，影像代價消耗越大。

1.3.4 基于Dijkstra的最短路徑搜索

Dijkstra 算法是典型的最短路徑算法［20－21］，用于計算一個節(jié)點到其他所有節(jié)點的最小代價。定義C為帶權(quán)有向圖的點集，設(shè)定源點并初始化S=0，S為用于存放已求出最短路徑的頂點集合，則C－S代表尚未求出最短路徑的頂點集合。

本文中點集C由所有的截面點集{V1，…，Vi，…，VN}組成。如圖6所示，假設(shè)在中心線的一端有一點設(shè)為源點，該點與V1點集中的所有點相連;此外，第二個截面中的一點vl分別與V1中的所有點相連。事實上，截面點集V2中的任意一點(如點vl)都與第一個截面中所有點相連，而其它任意2個截面之間點的連接方式也與此類似。

圖6 連接圖Fig.6 Connection image

候選線的傳播代價按式(2)計算。非鄰近的截面點不相連，權(quán)重設(shè)為∞。本文的目的在于搜索出源點到各個節(jié)點的最短路徑，當(dāng)?shù)竭_截面Vi時，便得到了從源點到當(dāng)前截面各點的最小代價，將該代價存于{D1，…，Di，…，DN}，其中 Di={di1，…，dij，…，dim}。

初始化中，將連接源點與V1中各點的候選線的權(quán)重設(shè)為0;添加V1至S。當(dāng)S=C，則執(zhí)行以下步驟:

1)對于Vi+1中的任一點vl，計算該點與Vi中各點相連所得候選線的傳播代價，記為El={e1l，…，ejl，…，em l};

2)計算Di+El，如果所得向量的最小值為dij+ejl，則代表從源點到Vi+1中點vl的最短路徑經(jīng)過Vi中的vj點，更新di+1，l=dij+ejl;

3)當(dāng)順序遍歷完Vi+1中的所有點，則添加Vi+1至S，并更新i=i+1。

由以上步驟獲取了從源點至所有截面點的最短路徑，則DN中的最小值即為優(yōu)化該條中心線的最小代價。找出該代價所對應(yīng)的各截面點并順次連接，即為優(yōu)化結(jié)果。

1.3.5 道路中心線斷開與重新連接

1.3.1 節(jié)中獲得的起始中心線確定了道路骨架的連續(xù)性和延展性，后續(xù)優(yōu)化過程是順序針對每條中心線的。以上操作并沒有考慮道路的方向一致性，因此在交叉口處存在不合理連接(圖7(a)中紅色橢圓內(nèi))。本文的策略是，先找到各個交叉口的近似中心點，然后以一定的半徑去除錯誤的交叉口連接(本文半徑取5個像素，圖7(b)給出了去除交叉口后的道路點方向圖)，最后根據(jù)方向的一致性和道路的連續(xù)性將這些斷開的道路段重新連接，得到最終的道路中心線(如圖7(c)所示，不同的顏色表示不同方向的道路)。

圖7 道路中心線斷開與重新連接Fig.7 Disconnection and reconnection of road central lines

2 實驗結(jié)果與分析

選取1景500像素×500像素的高分辨率航空影像作為實驗圖像(圖8)，該圖像覆蓋的區(qū)域為英國倫敦某一城區(qū)。

圖8 原始圖像Fig.8 Original image

首先采用概率增進樹算法獲得初始道路點。由于得到的初始道路點是屬于道路的概率，所以非道路點也會有值，但都比較小(圖9(a)中綠色橢圓內(nèi));然后結(jié)合路徑形態(tài)學(xué)和其他形態(tài)學(xué)運算進行處理，得到光滑完整的道路條帶(圖9(b))。

圖9 道路條帶獲取Fig.9 Road stripe acquision

圖10(a)給出對圖9(b)細(xì)化運算得到的骨架圖;圖10(b)為去除“毛刺”后的結(jié)果，但道路局部方向不一致，線條不夠平滑。根據(jù)方向信息構(gòu)建截面，并擬合得到參考方向;基于Dijkstra算法解算最小代價路徑，優(yōu)化后的道路中心線與原始圖像疊置(圖10(c))的效果很好。

圖10 道路中心線獲取Fig.10 Road centralline acquision

利用下列公式進行道路條帶和中心線提取精度的定量評價，即

式(4)—(6)中:Rm為人工解譯的道路;Rc，Rl和Re分別為正確提取的道路、遺漏的道路和錯誤提取的道路。

根據(jù)上述評價指標(biāo)，對本文方法的道路提取精度做出統(tǒng)計(表1)。

表1 本文方法道路提取結(jié)果統(tǒng)計Tab.1 Statistic results of road extracted with method proposed in this paper (%)

從表1可看出，提取出的道路條帶和中心線的準(zhǔn)確率均大于90%，提取結(jié)果較為理想。

3 結(jié)論

本文給出的從高分辨率遙感影像中提取道路中心線的方法，其創(chuàng)新點在于借鑒概率增進樹算法，并將其應(yīng)用于道路提取領(lǐng)域，獲得了完整的道路段;基于迭代求取最長測地距離的思想，獲取初始道路中心線，構(gòu)建能量函數(shù);利用Dijkstra算法優(yōu)化道路中心線，最終得到了精確的道路中心線。試驗結(jié)果證明了該策略的有效性和可行性。

本文尚未驗證該方法在大幅影像中的性能;另外，將本文策略用于道路邊緣提取領(lǐng)域，還沒有考慮立交路、環(huán)形路等形狀特別復(fù)雜的道路。這些將是今后進一步研究的內(nèi)容。

［1］ 史文中，朱長青，王 昱.從遙感影像提取道路特征的方法綜述與展望［J］.測繪學(xué)報，2001，30(3):257－262.ShiW Z，Zhu CQ，Wang Y.Road feature extraction from remotely sensed image:Review and prospects［J］.Acta Geodaetica et Cartographica Sinica，2001，30(3):257－262.

［2］ 趙曉峰.高分辨率遙感影像城區(qū)道路提取方法研究［D］.武漢:華中科技大學(xué)，2010:1－2.Zhao X F.A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of Master of Engineering［D］.Wuhan:Huazhong University of Science and Technology，2010:1－2.

［3］ 朱長青，王耀革，馬秋禾，等.基于形態(tài)分割的高分辨率遙感影像道路提?。跩］.測繪學(xué)報，2004，33(4):347－351.Zhu C Q，Wang Y G，Ma Q H，et al.Road extraction from high－resolution remotely sensed image based on morphological segmentation［J］.Acta Geodaetica et Cartographic a Sinica，2004，33(4):347－351.

［4］ 李 卉.集成LiDAR和遙感影像城市道路提取與三維建模［J］.測繪學(xué)報，2011，40(1):133－133.Li H.Road extraction and modeling with LiDAR and RS image in urban area［J］.Acta Geodaetica et Cartographic a Sinica，2011，40(1):133－133.

［5］ Zhang Q P，Couloigner I.Accurate centerline detection and line width estimation of thick lines using the radon transform［J］.IEEE Transactions on Image Processing，2007，16(2):310－316.

［6］ 周紹光，陳 超，岳建平.形狀先驗和圖割的高分辨率遙感影像道路段提?。跩］.測繪學(xué)報，2014，43(1):60－65.Zhou SG，Chen C，Yue JP.Extracting roads from high－resolution RS images based on shape priors and graph cuts［J］.Acta Geodaetica et Cartographic a Sinica，2014，43(1):60－65.

［7］ 曾發(fā)明，楊 波，吳德文，等.基于Canny邊緣檢測算子的礦區(qū)道路提取［J］.國土資源遙感，2013，25(4):72－78.doi:10.6046/gtzyyg.2013.04.12.Zeng FM，Yang B，Wu DW，et al.Extraction of roads in mining area based on Canny edge detection operator［J］.Remote Sensing for Land and Resources，2013，25(4):72－78.doi:10.6046/gtzyyg.2013.04.12.

［8］ 陳 卓，馬洪超，李云帆.結(jié)合角度紋理信息和Snake方法從Li-DAR點云數(shù)據(jù)中提取道路交叉口［J］.國土資源遙感，2013，25(4):79－84.doi:10.6046/gtzyyg.2013.04.13.Chen Z，Ma H C，Li Y F.Extraction of road intersection from Li-DAR point cloud data based on ATSand Snake［J］.Remote Sensing for Land and Resources，2013，25(4):79－84.doi:10.6046/gtzyyg.2013.04.13.

［9］ 蔡紅玥，姚國清.基于分水嶺算法的高分遙感圖像道路提取優(yōu)化方法［J］.國土資源遙感，2013，25(3):25－29.doi:10.6046/gtzyyg.2013.03.05.CaiH Y，Yao GQ.Optimized method for road extraction from high resolution remote sensing image based on watershed algorithm［J］.Remote Sensing for Land and Resources，2013，25(3):25－29.doi:10.6046/gtzyyg.2013.03.05.

［10］ Zhang Q P，Couloigner I.Accurate centerline detection and line width estimation of thick lines using the Radon transform［J］.IEEE Transactions on Image Processing，2007，16(2):310－316.

［11］ Miao Z L，Shi W Z.Road centreline extraction from classified images by using the geodesic method［J］.Remote Sensing Letters，2014，5(4):367－376.

［12］ Soille P.Morphological Image Analysis:Principles and Applications［M］.New York:Springer Verlag，1998.

［13］ Ozertem U，Erdogmus D.Locally defined principal curves and surfaces［J］.Journal of Machine Learning Research，2011，12:1249－1286.

［14］ Genovese C R，Perone－Pacifico M，Verdinelli I，et al.Nonparametric ridge estimation［J］.Institute of Mathematical Statistics，2014，42(4):1511－1545.

［15］ Tu ZW.Probabilistic boosting－tree:Learning discriminative models for classification，recognition，and clustering［C］//Proceedings of the Tenth IEEE International Conference on Compute Vision.Beijing:IEEE，2005:1589－1596.

［16］ Randen T，Husoy JH.Filtering for texture classification:A comparative study［J］.IEEE Transactions on Pattern Analysis ＆ Machine Intelligence，1999，21(4):291－310.

［17］ Meli L，Santiago JM，Lodge TP.Path－dependent morphology and relaxation kinetics of highly amphiphilic diblock copolymer micelles in lonic liquids［J］.Macromolecules，2010，43(4):2018－2027.

［18］ Deo N，Pang C Y.Shortest－path algorithms:Taxonomy and annotation［J］.Networks，1984，14(2):275－323.

［19］ Sperling G.Binocular vision:A physical and a neural theory［J］.The American Journal of Psychology，1970，83(4):461－534.

［20］ ShiW Z，Zhu CQ.The line segment match method for extracting road network from high－resolution satellite images［J］.IEEE Transaction on Geoscience and Remote Sensing，2002，40(2):511－514.

［21］ Ahuja R K，Magnanti T L，Orlin JB.Network Flows:Theory，Algorithms and Applications［M］.Englewood Cliffs，NJ:Prentice－Hall，1993.