結合路徑形態(tài)學的高分遙感影像道路提取方法

戴激光，苗志鵬，葛連茂，王曉桐，朱婷婷

(1.東華理工大學 江西省數字國土重點實驗室，南昌 330013；2.遼寧工程技術大學 測繪與地理科學學院，遼寧 阜新 123000)

0 引言

利用遙感影像進行道路信息采集在交通管理、城市規(guī)劃、自動車輛導航、地理信息系統(tǒng)數據庫更新、制作電子地圖等方面有重要應用價值，針對該領域的研究一直是國內外學者研究的熱點。隨著空間分辨率的提高，高分辨率影像可對道路幾何光譜特征進行更加精細地刻畫；但另一方面，冗余的細節(jié)信息諸如房屋、隔離帶、陰影、車輛等因素均會對道路提取構成巨大的干擾，進而大幅增加了道路提取的難度。

Vosselman等[1]對道路幾何與紋理信息進行了描述，道路路面存在一定的紋理相似性，并且道路具有長度較長的幾何特性。作為一種形態(tài)學方法，路徑形態(tài)學可提取圖像中長而窄的結構，這一特點亦符合道路的特征，故很多學者將路徑形態(tài)學方法應用于道路提取領域[2]。例如Talboth等人[3]提出不完整的路徑開運算，能越過一些“缺失的像素”搜索路徑，可應用于陰影遮擋、紋理變化較大或破損道路的提取，但該方法的缺點是效率低并且誤提取率較高。劉小丹等人[4]提出一種將Hough變換與路徑形態(tài)學結合的道路提取方法，通過Hough變換構造合適的結構元，再利用路徑形態(tài)學提取道路。但是Hough變換產生的結構元方向單一，難以提取曲率較大的道路。Appleton等[5]在Heijmans等[2]方法基礎上設計了一種有效的分解方法，將空間復雜度降低到對數級。同時，設定了4種三鄰接結構元，可提取曲率較大的道路。Schubert等[6]改進了Appleton等[5]提出的方法結構，增加4種對角線方向的結構元，可提取更多方向的道路。王雙等[7]設計了2種四鄰接結構元，可提取曲率更大的道路，但同時也會連接較多的非道路區(qū)域。

路徑開運算應用于道路提取工作中，雖然取得了較好的實驗成果，但由于該方法僅考慮道路紋理和長度特征，并未顧及道路寬度一致性和連通性，因而提取結果中存在道路與非道路區(qū)域粘連和斷裂問題。通過深入研究后發(fā)現引發(fā)這些問題的主要因素是：(1)道路兩側界限不清晰。當道路與其周圍地物間存在“異物同譜”現象時，路徑開運算便將其視為一個整體，導致道路與非道路區(qū)域粘連問題。(2)路面紋理信息不均一。路面上的車輛、人行橫道、建筑物或樹木陰影會造成局部紋理信息發(fā)生突變，這使得路面紋理一致性的特點不再突出，直接影響影像分割結果，從而易引發(fā)道路斷裂。有鑒于此，針對上述問題，在傳統(tǒng)路徑形態(tài)學的基礎上，本文提出了一種改進道路提取方法。

1 原理

1.1 數學形態(tài)學

數學形態(tài)學(mathematical morphology)是建立在集合論、拓撲論、隨機集等眾多學科的理論基礎之上描述形狀結構的科學。其基本思想是將攜帶形狀、大小、灰度等信息的結構元作為“探針”遍歷圖像，從圖像中“探測”出相應信息。根據輸入圖像的不同，數學形態(tài)學可以分為二值形態(tài)學和灰度形態(tài)學，兩者都包含膨脹、腐蝕、開運算、閉運算、頂帽運算、黑帽運算、擊中擊不中變換等運算。

由于路面上存在路面油漬、油返砂等問題，這使得道路內部光譜相似性降低，并將直接影響道路提取結果，所以本文使用灰度形態(tài)學對影像進行預處理，以減弱或消除這些干擾。根據上述對路面問題的分析，對于深色路面上的淺色干擾信息，本文采用腐蝕或開運算處理；對于淺色路面上的深色干擾信息，本文利用膨脹或閉運算處理。

1.2 Otsu雙閾值分割

當深色路面存在陰影等深色干擾時，路面灰度值的范圍介于淺色地物和深色干擾(陰影等)的灰度值之間，因此需用區(qū)間閾值才能獲得道路區(qū)域，故本文采用Otsu雙閾值分割方法對影像進行分割。利用Otsu方法得到2個閾值，將灰度值處于兩閾值之間的像素視為道路(灰度值為255)，其余視為非道路(灰度值為0)，從而得到二值圖像。

Otsu雙閾值分割是基于灰度直方圖得到各分割特性值發(fā)生的概率，用兩個閾值變量T1、T2將影像分為三類(C0、C1和C2)，其灰度值范圍分別是[0，T1]、[T1+1，T2]、[T2+1，m]，然后求出每一類的類內方差及類間方差，選取類間方差最大或類內方差最小的T1、T2作為最佳閾值[8]。其公式如下：

(1)

1.3 路徑開運算

路徑開運算(path opening)是路徑形態(tài)學的一種基本算法，屬于形態(tài)學的范疇。形態(tài)學中的開運算利用結構元遍歷圖像，可以剔除小于結構元大小的像素塊。而路徑開運算則利用結構元設定連接規(guī)則，將灰度值相近且鄰接的像素點連接成長度不同的路徑，然后剔除長度小于閾值L路徑點。定義如下：在二值圖像E中定義一種連接關系：ab，表明由a到b存在一條路徑。其中，a是b的后繼節(jié)點，b是a的前驅節(jié)點，符號“”表明路徑方向由a到b。將這種連接關系和圖像E中所有像素點Ω結合在一起(Ω，)，即為鄰接圖(directed acyclic graphs，DAG)。鄰接圖可以清晰地表示像素點間的連接關系。圖1(a)中，黑色部分為一種結構元定義的連接關系(方向為S-N)，灰色部分為該結構元構成的鄰接圖。其中，任一像素點F0最多擁有S1、S2、S3 3個后繼節(jié)點，任一點S0最多擁有F1、F2、F3 3個前驅節(jié)點。

定義1將每個點x關于連接關系“”的后繼節(jié)點的集合，記為δ({x})，即δ({x})={y∈Ω|yx}；x關于連接關系的前驅節(jié)點的集合，記為即y}。設X為前景圖像，對于?X?Ω，有：

δ(X)={y∈Ω|xyfor somex∈X}

(2)

(3)

定義2設Ta=(a1，a2，…，ak)為a的k元組，那么當且僅當aiai+1，且i∈[1，k- 1]時稱Ta為a長度為k的路徑。

定義3設Пk為所有長度為k的路徑集合，設路徑Ta里的所有像素點的集合為σ(Ta)。路徑開運算的定義為：

Ak(X)=∪{σ(Ta)|Ta∈Пkandσ(Ta)?X}

(4)

Ak給出前景圖像X中路徑長度大于或等于k的像素點集合，即為路徑開運算結果。圖1為路徑開運算示意圖，圖1(b)是一些前景黑色像素點在圖1(a)給出連接關系下的所有路徑。其中，每個點旁邊的數字是該點所在路徑的長度值(在一種連接關系下，若一個點同時存在于多條路徑上，規(guī)定其屬于最長的路徑)。圖1(c)為長度閾值L為6時路徑開運算結果，對比圖1(b)可以看出，路徑長度大于或等于6的像素點被保留下來，其余像素點被去除。當L為7時，僅保留了路徑長度大于或等于7的像素點(圖1(d))。

圖1 路徑開運算原理

基于以上分析，結合影像中道路區(qū)域與非道路區(qū)域的特點，發(fā)現道路區(qū)域的路徑長度一般長于非道路區(qū)域。因此將路徑開運算應用于道路提取中，通過人工設定合適的路徑長度閾值L，就可以有效地剔除非道路區(qū)域，保留道路區(qū)域。

本文的路徑開運算方法來自Schubert等方法，該方法包含8種三鄰接結構元(圖2)。其中每個點有3個連接方向，表示在搜索路徑時有角度方向上的容差，可提取具有一定曲率的道路。每種結構元可以探測一定方向的道路，不同的探測結果進行組合即可獲得完整道路網(圖3(c))。但由于前文所述的原因，有些道路提取結果不盡人意。如圖3(d)所示，影像中部存在的樹木陰影導致了路徑開運算結果中的道路斷裂問題；影像中的房屋與道路灰度近似并且鄰接，導致結果中存在道路與非道路區(qū)域的粘連問題(圖3(f))。

圖2 鄰接圖

圖3 路徑開運算提取道路

1.4 非道路粘連削減方法

針對上述路徑開運算中存在的問題，本文提出一種非道路粘連削減方法。道路與非道路粘連處(以下簡稱粘連處)可分為3個區(qū)域：道路、非道路、連接區(qū)域。連接區(qū)域即道路與非道路銜接區(qū)域。如圖4(a)所示，可以看到粘連處存在以下2個特征：(1) 連接區(qū)域寬度較小。(2) 非道路區(qū)域的路徑長度小于道路區(qū)域的路徑長度。

基于以上特征，本文首先設定大小為3×3的矩形結構元，通過形態(tài)學開運算將連接區(qū)域斷開，得到孤立的非道路區(qū)域；其次設定路徑長度閾值L進行路徑開運算。該閾值可依據實際影像區(qū)域道路目視分析，由人工進行設定；最后使用閉運算將道路網復原。

該方法有一個重要的前提：連接區(qū)域寬度應當很小，這樣形態(tài)學開運算才能在不影響道路區(qū)域的前提下斷開連接區(qū)域。若連接區(qū)域很寬，那么此處粘連無法去除。圖4展示了一條道路線上的兩處粘連問題的處理結果。圖4(a)中矩形框內的粘連，其連接區(qū)域寬度較小，本文方法可以將其去除。而圓圈內的粘連由于其連接區(qū)域寬度較大，因而無法進行刪除。

圖4 削減粘連示意圖

1.5 斷裂端點判讀

由于路徑開運算只能探測相鄰區(qū)域灰度值近似的像素，無法跨越灰度值相差很大的像素，故僅依靠形態(tài)學方法難以解決道路斷裂問題。由于模板跟蹤在道路提取過程中對道路光譜變化有一定的魯棒性[10]。因此本文在路徑開運算的結果中，取其斷裂處的端點作為種子點，并由此提出改進的圓形模板跟蹤方法，以此完成斷裂處道路中心線的跟蹤，但在開展該工作前首先需要確定斷裂處端點。

考慮到路徑開運算提取的道路結果中存在很多端點，因此在道路跟蹤前首先需要判斷哪2個端點屬于同一斷裂區(qū)域。仔細觀察道路斷裂處存在下述2個特征：(1) 道路斷裂是由陰影、人行橫道等造成的，因而道路斷裂處長度較短。因此可依據這一特征，比較所有端點的間距，找出間距較小的端點對。(2) 斷裂處曲率變化較小。即斷裂處兩側線段的角度差異應當較小。依據上述特征，可采用最小二乘法將兩端點所在區(qū)域進行擬合，提取2條最佳擬合線段，同時將2個端點連接起來形成連接線段，分析3條線段的角度差異是否低于π/4。若2個端點同時滿足上述2個特征，則認為這2個端點應屬于同一斷裂區(qū)域。

1.6 改進的圓形模板匹配跟蹤方法

本文提出的跟蹤方法是對文獻[10]方法的改進。文獻[10]方法的思想是：人工選取若干種子點，在兩兩種子點間以一定規(guī)則內插新種子點，不斷重復內插過程直至達到閾值。該方法用直線段逼近曲線，對曲率較小的道路提取效果較好。其優(yōu)點在于，若前幾次內插的種子點精度較高，隨著種子點的不斷增多，其形狀會不斷逼近道路形狀，其計算過程會越來越快，精度會越來越高。其缺點在于：(1) 若初始種子點間道路曲率過大，第一次內插過程的難度就會比較高。(2) 已有內插點的精度會影響后續(xù)內插點精度，進而影響道路提取質量。

基于上述問題，本文提出改進的圓形模板提取方法。與原方法主要不同之處是：(1) 改進了內插種子點的方式，將其對精度影響均勻分布到每個點上。(2) 加入步長增加措施，可越過“障礙”繼續(xù)搜索道路。具體流程如下：

(1)采用L0濾波方法[11]處理輸入影像，利用3×3矩形結構元計算每個點的形態(tài)學梯度。

(2)將上文得到的斷裂處兩端點作為初始種子點。以一個初始種子點及其八鄰域內的點為圓心，分別建立半徑為1像素的圓形模板，計算模板內點的梯度之和Gradsum。將同一半徑下Gradsum最小的模板視為最優(yōu)模板，并判斷其Gradsum是否超過閾值δ(經驗閾值，本文設為200)，若不超過δ，則以1個像素的步長逐步擴大模板半徑。繼續(xù)該搜索過程直至最優(yōu)模板Gradsum超過δ，將前一半徑的最優(yōu)模板設為標準模板，其圓心p1為道路中心點，其半徑r為道路寬度的一半(圖5(a)過程①)。按照上述流程，在第二個初始種子點處尋找半徑為r的標準模板，獲得道路中心點p2(如圖5(a)過程④)。將p1到p2的方向設為初始跟蹤方向。

(3)若p1到p2的距離小于r，連接上述兩點并停止跟蹤。否則從初始端點p1開始，在跟蹤方向上長度為r的點處，垂直于跟蹤方向上建立一系列標準模板(圖5(a)中過程②的Ti)。定義模板角度θmod為當前模板圓心P與前一個模板圓心M連線的角度(圖5(a)中的θ)。

定義灰度均值差為模板P中每個點q的灰度值與前一個模板M灰度均值之差的絕對值之和，公式如下所示：

(5)

定義梯度方差為模板P中每個點q的梯度Gq與模板P梯度均值之差的平方之和的平方根，公式如下所示：

(6)

定義模板角度差為模板P的θmod與其前6個模板θmod均值(不足6個按實際數量計算)之差的絕對值，公式如下所示：

(7)

分別計算Ti的灰度均值差、梯度方差、以及模板角度差，將3個量歸一化后加權求和(公式(8))，求得權重函數取得最大值時的圓心坐標，將其作為待定道路中心點。若2個模板的權重函數值相同，取模板角度差較小的一個；

(4)為驗證待定道路中心點是否是道路中心，本文采用邊緣判斷法對其進行判斷。由于路面與道路兩側存在灰度差異，所以在待定道路中心點模板A兩側，垂直于跟蹤方向上建立標準模板B、C(圖5(b))，分別計算B與A、C與A的灰度均值之差的絕對值。若任一差值大于閾值Dgary(本文設為30)，說明該點為道路中心點，并將其作為新的初始端點，其模板角度為跟蹤方向，重復步驟(3)；若差值均小于閾值Dgary，說明該步長下求得的待定點存在異常(圖5(a)過程③中障礙物)，需要進一步判斷，進入步驟(5)；

(5)將步長變?yōu)槎?，重復步驟(3)、步驟(4)。若二倍步長無法滿足要求，則變?yōu)槿恫介L(圖5(a)過程③)。若三倍步長也無法滿足要求，停止跟蹤。

圖5 道路邊緣判斷、斷裂處跟蹤示意圖

(8)

1.7 方法步驟

通過上述對不同方法原理的分析，本文給出了具體的方法處理流程(圖6)。

圖6 技術流程圖

(1) 采用灰度形態(tài)學對輸入高分辨率遙感影像進行處理，以消除路面噪聲；

(2) 利用Otsu方法分割影像，得到包含道路區(qū)域和非道路區(qū)域的二值圖；

(3) 采用路徑開運算消除大部分的非道路區(qū)域；

(4) 設定道路長度最低閾值L，利用形態(tài)學和路徑開運算削減粘連的非道路區(qū)域；

(5) 對斷裂端點進行判斷，若存在則轉入步驟(6)，否則轉入步驟(7)；

(6) 在L0濾波影像基礎上，以斷裂端點為種子點，運用改進圓形模板匹配方法跟蹤道路區(qū)域，進入步驟(7)；

(7) 輸出道路結果。

2 實驗結果分析和算法評價

2.1 實驗影像及實驗方法介紹

本文基于VS2013平臺，采用C++編程語言對3幅遙感影像進行實驗。其中，Pleiades影像大小為632像素×679像素，空間分辨率為0.5 m；2幅高分二號影像大小均為1 000像素×1 000像素，空間分辨率為0.8 m。本文首先分別利用TALBOT 等方法中“不完整路徑開運算”、Appleton等方法中“有效的路徑開運算和閉運算”以及Schubert等方法中“有效的灰度路徑開運算”進行路徑開運算實驗，闡述斷裂問題和粘連問題。然后在Schubert等方法的結果中，加入本文的后處理方法，以驗證本文方法的有效性。其中，TALBOT 等方法提供了該方法的C++代碼，Schubert等方法提供了Schubert等方法和Appleton等方法的C++代碼，因而對比實驗結果是可信的。

2.2 實驗一

實驗一選取Pleiades拍攝的農村遙感影像。該影像中道路與居民區(qū)有一定的粘連，路面有樹木陰影遮擋；一些居民區(qū)呈帶狀分布，灰度值與道路相近，形成了類似道路的結構，對于道路識別構成很大的干擾。

基于如圖7(b)所示二值圖，3種方法開展的道路提取實驗，結果分別如圖7(c)、圖7(d)、圖7(e)所示?？梢钥吹剑捎诋斍坝跋窀采w區(qū)域比較復雜，TALBOT 等方法雖然將道路轉彎處提取出來，也解決了樹木陰影遮擋對道路遮擋的問題，但仍然存在很多粘連，右側道路誤提取率較大；而Appleton等方法在道路轉彎區(qū)域提取效果不佳，并且將成排的房屋視為道路；而Schubert等方法提取道路結果精度較高，不僅能夠提取出彎曲的道路并且誤提取率低，但依然存在非道路與道路粘連問題和道路斷裂問題。因此本文以Schubert等方法為基礎，實驗證明是可行的。如圖7(f)所示，依據本文提出的削減粘連方法，經過一系列形態(tài)學與路徑開運算相結合的操作(圖8)，削減了大部分的粘連區(qū)域，仍有一些粘連由于連接處較寬而無法削減，并且產生了一個新的道路斷裂(圖8(c)左下)。如圖7(g)斷裂處放大顯示，改進圓形模板跟蹤的道路中心線與實際道路中心線基本吻合。最后用道路中心線將斷裂處連接起來即可得到該影像的道路網(圖7(h))。

圖7 Pleiades影像實驗

圖8 削減粘連

2.3 實驗二

實驗二選取高分二號城區(qū)遙感影像。高大建筑物產生的陰影將影像上的路面區(qū)域截為兩段(圖9(a)下部)，因而道路內部光譜信息呈現出不一致的現象。同時由于道路上存在大量的車輛，以及較多的人行橫道信息，這均會導致二值圖中出現孔洞。并且右下角橫向道路兩側存在與道路灰度相近的植物以及部分低矮建筑的陰影，這將會造成道路與非道路的粘連。

圖9 高分二號影像實驗

采用灰度形態(tài)學對圖像預處理，采用Otsu方法對影像進行二值化處理，采用3種方法進行路徑開運算，其結果分別如圖9(c)、圖9(d)、圖9(e)所示。由于該影像道路結構比較簡單，3種路徑開運算方法提取道路的完整性很高，但是提取質量各不相同。其中，TALBOT 等方法在能夠較好地提取部分道路，但誤提取率很高。Appleton等方法由于其結構元比較適合提取直線，其誤提取率最低，但完整性略低(圖9(d)下部道路缺失)。Schubert等方法誤提取率介于其他2種方法之間，完整性較好。對Schubert等方法的結果進行削減粘連后，發(fā)現路面仍存在空洞，用二值形態(tài)學閉運算對空洞進行填補，并通過圓形模板匹配跟蹤方法對道路斷裂處進行連接，得到道路提取結果(圖9(f))。

2.4 實驗三

實驗三選取高分二號城市郊區(qū)遙感影像。對影像的目視分析可以發(fā)現，場景中存在的建筑物陰影會造成道路線斷裂，一些道路兩側界線不清晰可能會造成道路與非道路的粘連。

采用Otsu方法對影像進行二值化處理后，如圖10(b)可以發(fā)現道路邊緣比較規(guī)整，橫向部分的道路雖然曲率較小，但斷裂情況依然非常明顯。3種路徑開運算提取結果分別如圖10(c)、圖10(d)、圖10(e)所示。其中，TALBOT 等方法雖對斷裂問題有所改善，但誤提取較高并且粘連問題較嚴重；Appleton等方法和Schubert等方法完整性相似，但Schubert等方法的粘連較少。如圖10(f)所示，針對道路三岔口處的斷裂問題，利用本文的方法，不僅能夠道路粘連處進行削減，同時能夠解決由陰影導致的道路斷裂問題，并如圖10(f)所示由此得到完整的道路提取結果。

圖10 高分二號影像實驗

2.5 算法的精度和效率評價

本文利用文獻[12]提出的完整性(complete percent，CP)、正確率(correct percent，CR)、提取質量(quality percentage，QL)3個指標對實驗結果進行量化分析評判，3個指標的計算方法如下所示。

(9)

(10)

(11)

式中：TP表示匹配的道路長度；FP表示未匹配的提取道路的長度；TN表示匹配的參考實際道路長度；FN表示未匹配的實際道路長度。本文通過統(tǒng)計結果中各部分對應的像素的個數，根據公式(9)、公式(10)、公式(11)，得到各類方法的精度統(tǒng)計如表1所示。

表1 不同方法精度統(tǒng)計表

由表1可以看出，Appleton等方法由于其采用的結構元數量較少(4個)，因而運行速度較快，但提取的完整性不高。Schubert等方法和TALBOT 等方法包含較多方向的結構元(8個)，因而提取道路的完整性較高，效率較低。另外，TALBOT 等方法還包含對道路斷裂問題的改善，因此運行速度最低，完整度最高但出現很多的誤提取問題。本文所用的Schubert等方法相較于Appleton等方法，優(yōu)勢在于可提取彎曲的道路并且道路提取完整性較好；相較于TALBOT 等方法，優(yōu)勢在于誤提取率較低、運行速度快。而本文方法在Schubert等方法基礎上進行改進，后處理結果的完整性、正確率和提取質量均表現良好。本文方法主要涉及3個步驟，分別是路徑開運算、削減粘連和連接斷裂。其中路徑開運算方法效率與二值圖像的像素數和結構元的數量有關。削減粘連過程效率取決于形態(tài)學運算和路徑開運算的方法效率。連接斷裂方法效率與斷裂區(qū)域的灰度穩(wěn)定程度有關。從整體來看，本文提出的改進步驟速度均控制在1 s之內，效率較高。

3 結束語

本文首先對基于路徑開運算的高分辨率道路提取結果存在的問題進行分析探討；然后提出削減道路粘連的處理策略以及解決道路斷裂問題的改進圓形模板跟蹤方法；最后選用3幅有典型意義的影像進行實驗。實驗結果表明本文方法對于解決道路粘連問題和道路斷裂問題有較好的效果。但本文方法也有不足之處，即需要人工設定路徑開運算閾值L，這也是路徑開運算難以避免的問題。因此如何改進路徑開運算方法，自適應設定閾值L，這將是未來研究的重點。