基于動態(tài)規(guī)劃的城市道路提取改進(jìn)算法

劉曉平 ， 朱曉強， 余 燁， 李 琳

（1. 特種顯示技術(shù)教育部重點實驗室（合肥工業(yè)大學(xué)），安徽 合肥 230009；2. 合肥工業(yè)大學(xué)計算機與信息學(xué)院VCC研究室，安徽 合肥 230009）

從航空或衛(wèi)星影像中提取道路一直都是一個很具有挑戰(zhàn)性的問題。根據(jù)道路特征提取的自動化程度，道路特征提取分為自動特征提取和半自動特征提取。自動提取方法包括基于平行線對[1]、基于二值化和知識[2]、基于窗口模型特征[3]等；半自動提取的常用算法有模板匹配[4]、動態(tài)規(guī)劃[6]、可變模型或Snakes方法[7]等。到目前為止，從航空或衛(wèi)星影像中全自動識別道路仍然處于不成熟的狀態(tài)[8]，而且在未來的短期時間內(nèi)也很難開發(fā)一套完全脫離人工干預(yù)而自動識別道路的系統(tǒng)。在各種半自動道路提取方法中，由于動態(tài)規(guī)劃是一個解決優(yōu)化問題的穩(wěn)定算法，能夠確保結(jié)果的最優(yōu)性，成為近年來道路半自動提取的首選方法之一。動態(tài)規(guī)劃法首先應(yīng)用在低分辨率影像提取道路中（Sakoda等，1993；Gruen和Li，1997；Merlet和 Zerubia，1996）。Dal Poz等[6]修改了“代價”函數(shù)，使之適用于從中高分辨率影像中提取道路。

另一方面，從機載LiDAR（Light Detection and Ranging）數(shù)據(jù)中提取道路在近年來也有了一定的研究，Hatger和 Brenner在文獻(xiàn)[10]中結(jié)合LiDAR數(shù)據(jù)和已經(jīng)存在的數(shù)據(jù)庫信息提取道路的幾何屬性；在Rieger等人的論文[11]中研究了林區(qū)范圍內(nèi)由 LiDAR點云提取道路方法；為充分利用LiDAR數(shù)據(jù)信息，文獻(xiàn)[12]和文獻(xiàn)[13]引入了強度信息輔助道路的提取。

由于從航拍圖像中提取道路時一些經(jīng)驗條件難以滿足，如要求道路的灰度比背景亮，道路的灰度變換不大等，而從 LiDAR數(shù)據(jù)中提取道路對點的要求比較高，如點集密度高、需要強度信息等，因此，本文融合 LiDAR點云和配準(zhǔn)的航空圖像兩種數(shù)據(jù)源來提取道路，降低了對單一數(shù)據(jù)源的依賴性。文章首先介紹了現(xiàn)有的基于動態(tài)規(guī)劃的道路提取方法和代價函數(shù)的意義；然后對 LiDAR點云數(shù)據(jù)進(jìn)行重新采樣，并根據(jù)提取的道路特征對原有代價函數(shù)進(jìn)行擴展；為精確配準(zhǔn)航拍圖像和 LiDAR數(shù)據(jù)，本文利用局部控制點對當(dāng)前位置進(jìn)行插值；最后通過試驗對本文算法進(jìn)行驗證。

1 基于動態(tài)規(guī)劃的道路提取及其改進(jìn)

1.1 現(xiàn)有基于動態(tài)規(guī)劃道路提取方法

1.1.1 在低分辨率航拍圖像中的應(yīng)用

利用動態(tài)規(guī)劃從航拍圖像中提取道路首先從較低分辨率的圖像中[14]，將道路建立成線狀的模型。道路模型中同時考慮道路的幾何屬性和物理屬性（如道路為線狀、灰度和寬度變化慢、方向變化平緩等）。由于道路可以用多邊形來表達(dá){ p1, p2,… ,pn}，其中 pi=( xi,yi)是多邊形在圖像上的第i個定點。則道路模型可以用下面的代價函數(shù)和約束不等式表示

其中 Ep1(pi,pi+1), Ep2(pi,pi+1)，和 Ep3(pi,pi+1)描述道路的各種幾何和物理屬性，每個表達(dá)式只依賴于連續(xù)鄰接的兩個點pi和pi+1；αi表示由點pi-1和pi構(gòu)成的向量方向；β和γ為正常數(shù)；|ΔSi|為pi-1和pi之間的距離；T為用戶定義的由控制點構(gòu)成的臨近兩個向量的方向改變量的上限。

分析式(1)可以知道，代價函數(shù)為一系列Ei的加權(quán)和，每一個Ei依賴于多邊形P的 3個連續(xù)的控制點(pi-1, pi,pi+1)，即每個點pi只和它的前一個pi-1、后一個點pi+1直接相關(guān)。這使得可以利用動態(tài)規(guī)劃來該問題，通過一系列的決策過程來達(dá)到優(yōu)化的目的。根據(jù)精度要求的不同，迭代的過程可以在原有的控制點之間添加一些額外的輔助控制點。

1.1.2 在中、高分辨率航拍圖像中的應(yīng)用

為使式(1)也能適合從高分辨率航空影像中提取道路，文獻(xiàn)[6]研究了中等和較高分辨率航空影像中道路中心線的提取方法，將道路的邊緣信息和寬度信息融入到式(1)得到了新的代價函數(shù)如下

1.2 代價函數(shù)的改進(jìn)

在式(3)中雖然融入了航空圖像中道路的各種屬性，但是很多條件是基于經(jīng)驗，當(dāng)條件改變的時候，該方法很難保證道路提取的正確性，如Ep中假定路面的灰度高于兩側(cè)的灰度值[14]，路

1面的灰度值變化緩慢，這些條件在大多數(shù)情況下是成立的，但是不能絕對保證（如斑馬線處灰度值變化劇烈等）。為提高算法的魯棒性，本文中引入了另外一個數(shù)據(jù)源——LiDAR點云，從而可以為模型式(3)添加更多的路面幾何信息，如路面的邊緣、高度等。

1.2.1 LiDAR數(shù)據(jù)的重采樣和特征提取

為使LiDAR數(shù)據(jù)能夠與航空圖像結(jié)合使用，并減少對數(shù)據(jù)密度的依賴性，本文對 LiDAR點云進(jìn)行重新采樣。重采樣的方法就是，按LiDAR點云的X、Y范圍進(jìn)行等網(wǎng)格劃分，每個網(wǎng)格對應(yīng)像素的灰度值正比于網(wǎng)格內(nèi)的點的平均高度。由于 LiDAR點云分布是雜亂無章、密度不均的（如圖1所示），所以生成的高度圖包含空洞[12]，空洞處像素用離最近的像素灰度值代替。至于網(wǎng)格的大小可以根據(jù)實際的點云密度來設(shè)定。該階段生成的高度圖由于灰度對應(yīng) LiDAR點集的高度信息，為后續(xù)處理和敘述的方便，記該高度圖為Image1（如圖2所示），并用邊緣檢測方法（如Canny方法）生成邊緣圖像Image2。因為路面是平滑的，檢測出的邊緣不會出現(xiàn)在路面上，而道路與兩旁的建筑物以及立交橋與兩旁的地面有明顯的高度差，因此在道路邊緣處會檢測到邊緣點的存在（如圖3所示）。

圖1 城市LiDAR點云數(shù)據(jù)

圖2 高度圖像

圖3 邊緣圖像

圖4 邊緣距離圖像

為確定道路中心線的精確位置，根據(jù)邊緣檢測圖像Image2生成邊緣距離圖像Image3，其不僅在求取道路中心線中用到，在最后求取道路寬度時也要使用該圖像中的信息。Image3生成方法為，對于 Image2中任一非邊緣點P，求取其最近的邊緣點E1，若E1有關(guān)于P對稱的邊緣點E2，則將P和E1的距離作為像素P點的灰度值，如圖4所示。

1.2.2 LiDAR數(shù)據(jù)道路特征

根據(jù)圖像Image1, Image2和Image3的生成過程和道路本身固有的屬性，可以得到如下有關(guān)道路特征的描述：

（1）小段范圍內(nèi)道路高度變化小，從而對應(yīng)Image1中像素的灰度值變化小，則有

其中 f(s)為道路曲線的表達(dá)式；s為道路曲線的弧長；Δsi為取定的一段道路；G可分別表示灰度圖像Image1和航拍圖像 Image4的灰度值；Gm( Δ si)為對應(yīng)圖像中在道路段Δsi上的平均灰度值，即

（2）路面是平滑的，在邊緣圖像 Image2中道路曲線一般不會橫跨邊緣點，也即 Image2中的邊緣點對代價函數(shù)E的貢獻(xiàn)應(yīng)該最小，為

其中 G表示圖像Image2對應(yīng)圖像灰度值。

（3）提取的道路應(yīng)離道路的中心線很近，則邊緣距離圖像中 Image3中每點的灰度值對代價函數(shù)的影響為

其中 G表示圖像Image3對應(yīng)圖像灰度值。

（4）道路呈現(xiàn)帶狀形態(tài)，在短距離內(nèi)不會有太大的寬度變化，基于此得出道路寬度變化對代價函數(shù)的影響

其中 G表示圖像Image3對應(yīng)圖像灰度值。

1.2.3 代價函數(shù)的擴充

綜合航空圖像和 LiDAR數(shù)據(jù)的信息，可以從式(3)～式(7)中得出最終的道路模型的代價函數(shù)

其中 Ep5(pi,pi+1)=μ1E1+μ2E2-μ3E3+μ4E4，

式(2)、式(8)共同組成了本文道路識別的動態(tài)規(guī)劃半自動提取模型。

2 航拍圖像與LiDAR數(shù)據(jù)配準(zhǔn)

由于距離較遠(yuǎn)和影像獲取環(huán)境的不穩(wěn)定等因素會使獲取的圖像有一定的畸變和傾斜，影響道路的提取效果和精度。因此，在結(jié)合航空圖像和 LiDAR點云提取道路之前需要對兩者進(jìn)行配準(zhǔn)。配準(zhǔn)過程需要選定一些控制點作為配準(zhǔn)依據(jù)，所提供的控制點個數(shù)和分布情況影響配準(zhǔn)結(jié)果。配準(zhǔn)的每一步都使用全部的控制點不僅計算量很大，而且當(dāng)控制點離當(dāng)前像素點較遠(yuǎn)時反而會降低配準(zhǔn)的精度。

本文利用基于局部插值的思想，首先對控制點進(jìn)行構(gòu)造Delaunay三角網(wǎng)，然后對每個三角形內(nèi)的像素只使用該三角形的頂點進(jìn)行配準(zhǔn)[15]，既減少了計算量，也能保證配準(zhǔn)的精確性。如圖5和圖6分別為配準(zhǔn)前的航拍圖像和由LiDAR數(shù)據(jù)生成的高度圖，以及根據(jù)選取的控制點所構(gòu)建的三角網(wǎng)，圖7為與圖6高度圖配準(zhǔn)的航拍圖像。

圖5 配準(zhǔn)前的航拍圖像

圖6 LiDAR高度圖

圖7 配準(zhǔn)后的航拍圖像

3 試驗結(jié)果與分析

本文算法是半自動提取方法，首先由用戶在道路附近選定一些控制點，之后的工作就由計算機來自動完成。它會自動算出道路代價函數(shù)的最大值，并根據(jù)列表保存的中間步驟信息，逆向求出使代價函數(shù)達(dá)到最大時的道路中心線位置。此外，通過道路中心線的位置和邊緣距離圖像Image3，算出道路中心線上任一位置的道路寬度；根據(jù)道路中心線的位置和高度圖 Image1的生成過程計算道路的高度。

所使用的試驗數(shù)據(jù)來源于在美國新奧爾良市區(qū)的LiDAR數(shù)據(jù)。圖2～圖4為LiDAR數(shù)據(jù)信息，圖7為該區(qū)域的配準(zhǔn)航拍圖像，圖8中為用戶選取的控制點位置，圖9顯示所提取的道路位置和寬度信息，圖10則為道路重建后的三維模型。從提取過程來看，用戶提供的控制點甚至不需要在道路上，該算法仍然能夠取得很好的效果。

圖8 控制點位置

圖9 提取的道路

圖10 道路重建后

4 結(jié)束語

本文首先總結(jié)了動態(tài)規(guī)劃方法在道路提取過程中的發(fā)展，其最先應(yīng)用于較低分辨率的航拍圖像中，將道路抽象為一條很細(xì)的線狀，為了將該方法應(yīng)用于中、高分辨率的航拍圖像中，研究人員引入了道路的寬度信息，從而能夠提取航拍圖像中的道路中心線。

為了進(jìn)一步提高算法的魯棒性，本文融合航拍圖像和LiDAR點云兩種數(shù)據(jù)源，對LiDAR點云數(shù)據(jù)進(jìn)行重新采樣，分別生成對應(yīng)的高度圖、邊緣檢測圖像和邊緣距離圖像；結(jié)合生成的三幅圖像和對應(yīng)的航拍圖像提取道路的特征并用于改進(jìn)后的動態(tài)規(guī)劃算法；最后通過試驗驗證了本文方法的可行性。

[1]Ton J, Jain A K, Enslin W R, et al. Automatic road identification and labeling in Landsat 4 TM images [J].Photogrammetric Engineering & Remote Sensing,1989, 43(2): 257-276.

[2]Wang F, Newkirk R. A knowledge-based system for highway network extraction [J]. IEEE Transactions on Geoscience and Remote Sensing, 1988, 26(5):525-531.

[3]Barzohar M, Cooper D B. Automatic finding of main roads in aerial images by using geometric-stochastic models and estimation [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1996,18(7): 707-721.

[4]楊 云, 朱長青, 張 德. 高分辨率遙感影像上道路中心線的半自動提取[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2007, 19(6): 781-785.

[5]German D, Jedynak B. An active testing model for tracking roads in satellite images [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1996, 18(1): 1-14.

[6]Pos A, Vale G. Dynamic programming approach for semi-automated road extraction from medium and high-resolution images [J]. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2003, 34(8): 87-91.

[7]Armin Gruen, Li Haihong. Semi-automatic linear feature extraction by dynamic programming and LSB-snakes [J]. Photogrammetric Engineering and Romote Sensing, 1997, 63(8): 985-995.

[8]Shackelford A K, Davis C H. Fully automated road network extraction from high-resolution satellite multispectral imagery [C]//Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, 2003: 461-463.

[9]Zhu C, Shi W, Pesaresi M, et al. The recognition of road network from high-resolution satellite remotely sensed data using image morphological characteristics [J].International Journal of Remote Sensing, 2005, 26(24):5493-5508.

[10]Hatger C, Brenner C. Extraction of road geometry parameters from laser scanning and existing databases[C]//International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol. XXXIV, Part 3/W13,2003.

[11]Rieger W, Kerschner M, Reiter T, et al. Roads and buildings from laser scanner data within a forest enterprise [C]//IAPRS, ISPRS Workshop, Vol. 32,La Jolla, California, 1999: 185-191.

[12]王 濤, 楊建思, 廖明生. 從機載激光掃描數(shù)據(jù)中提取道路[J]. 測繪信息與工程, 2006, 31(5): 10-11.

[13]Simon Clode, Peter Kootsookos, Franz Rottensteiner.The automatic extraction of roads from lidar data [C]//International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXV-B3,2004: 231-236.

[14]Armin Gruen, Li Haihong. Road extraction from aerial and satellite images by dynamic programming [J].JSPRS Journal of Photogrammetry and Remote Sensing, 1995, 50(4): 11-20.

[15]羅月童, 朱曉強, 劉曉平. 基于三點的航拍圖像傾斜校正方法及其應(yīng)用[C]//中國儀器儀表學(xué)會第九屆青年學(xué)術(shù)會議論文集, 2007: 590-593.