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

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

0 引言

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

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

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

1 道路中心線提取算法

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

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

1.1 道路條帶的獲取

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

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

1.2 初始道路中心線的獲取

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

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

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

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

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

1.3.1 直線段獲取

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

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

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

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

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

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

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

1.3.2 道路截面確定

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

1.3.3 路徑代價(jià)設(shè)計(jì)

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

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

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

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

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

角度差異越大，代價(jià)越大。

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

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

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

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

1.3.4 基于Dijkstra的最短路徑搜索

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

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

圖6 連接圖Fig.6 Connection image

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

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

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

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

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

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

1.3.5 道路中心線斷開(kāi)與重新連接

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

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

2 實(shí)驗(yàn)結(jié)果與分析

選取1景500像素×500像素的高分辨率航空影像作為實(shí)驗(yàn)圖像(圖8)，該圖像覆蓋的區(qū)域?yàn)橛?guó)倫敦某一城區(qū)。

圖8 原始圖像Fig.8 Original image

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

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

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

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

利用下列公式進(jìn)行道路條帶和中心線提取精度的定量評(píng)價(jià)，即

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

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

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

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

3 結(jié)論

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

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

［1］ 史文中，朱長(zhǎng)青，王 昱.從遙感影像提取道路特征的方法綜述與展望［J］.測(cè)繪學(xué)報(bào)，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］ 朱長(zhǎng)青，王耀革，馬秋禾，等.基于形態(tài)分割的高分辨率遙感影像道路提?。跩］.測(cè)繪學(xué)報(bào)，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］.測(cè)繪學(xué)報(bào)，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］ 周紹光，陳 超，岳建平.形狀先驗(yàn)和圖割的高分辨率遙感影像道路段提?。跩］.測(cè)繪學(xué)報(bào)，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邊緣檢測(cè)算子的礦區(qū)道路提?。跩］.國(guó)土資源遙感，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點(diǎn)云數(shù)據(jù)中提取道路交叉口［J］.國(guó)土資源遙感，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］ 蔡紅玥，姚國(guó)清.基于分水嶺算法的高分遙感圖像道路提取優(yōu)化方法［J］.國(guó)土資源遙感，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.