基于Grabcut融合多源數(shù)據(jù)提取不透水面

易佳思， 胡翔云

(武漢大學遙感信息工程學院，武漢 430079)

0 引言

近年來城市防災減災任務艱巨，城市暴雨澇災常發(fā)。城鎮(zhèn)建設面積迅速擴張、以植被為主的自然景觀逐漸被人工不透水面所取代，是導致城市發(fā)生暴雨內(nèi)澇的重要影響因素[1]。不透水面指天然或人造的難以被水穿透的地表，城市不透水面主要由人工地物構(gòu)成，如屋頂、停車場、水泥及瀝青路面等。不透水面作為衡量城市化程度和環(huán)境質(zhì)量的重要指標，受到越來越多的關(guān)注。研究表明，遙感數(shù)據(jù)在地物提取方面有著很大的潛力和優(yōu)勢，利用遙感手段提取不透水面已成為研究熱點之一[2]。

目前常用的不透水面遙感提取方法有人工解譯法、指數(shù)模型法、回歸分析法、線性光譜混合模型法、決策樹分類法、人工神經(jīng)網(wǎng)絡和面向?qū)ο蠓ǖ取＿@些研究取得了較好的成果，但是在精確高效提取不透水面方面仍存在不足。如： 基于V-I-S模型利用線性光譜混合分解將光譜信息分解成高反照度、低反照度、植被和土壤端元，但是由于光譜相似造成端元之間存在不確定性[3-5]； 回歸分析法對于不同時相不同區(qū)域難以實現(xiàn)自適應歸一化，且對數(shù)據(jù)噪聲敏感[3,6-7]； 人工智能算法在遙感影像上應用還不夠成熟，網(wǎng)絡連接權(quán)值的物理意義不易推斷，學習收斂速度慢，且對訓練樣本的數(shù)量和質(zhì)量有較高要求[8]； 面向?qū)ο蠓m然包含了形狀、紋理和上下文關(guān)系等空間信息，但是分割尺度的確定對分類精度影響很大[9]。基于上述方法，也開展了大量改進研究，以提高算法的魯棒性和精度。同時，傳統(tǒng)不透水面提取的研究方法只采用單一衛(wèi)星數(shù)據(jù)源，沒有利用多源數(shù)據(jù)的優(yōu)勢。因此引入了夜間數(shù)據(jù)、合成孔徑雷達數(shù)據(jù)(synthetic aperture Radar，SAR)、機載激光點云數(shù)據(jù)(light detection and ranging，LiDAR)等新型數(shù)據(jù)。將這些數(shù)據(jù)與傳統(tǒng)遙感影像相結(jié)合[10-11]，已成為了當前不透水面提取的重要手段。Hodgson等[12]利用航片數(shù)據(jù)和LiDAR數(shù)據(jù)，對城區(qū)不透水面進行了提取，相比使用傳統(tǒng)分類方法精度得到了提高； Im等[13]使用WorldView數(shù)據(jù)和LiDAR數(shù)據(jù)對城市不透水面進行了提取，既豐富了信息源，也提高了不透水面的提取精度。

廣州市是國內(nèi)城市化發(fā)展最為快速的地區(qū)之一，曾多次遭遇連續(xù)強降雨，導致市區(qū)大面積內(nèi)澇，因而實時準確獲取廣州市不透水面分布情況迫在眉睫。以該地區(qū)作為研究區(qū)，把建筑物、水泥瀝青道路和停車場等人工構(gòu)造物歸為不透水面，把植被、裸地和水系歸為透水面。基于圖論的最優(yōu)分割方法(Grabcut)，將不透水面的提取轉(zhuǎn)換為不同數(shù)據(jù)源下的最優(yōu)標記問題，建立能量函數(shù)，融合Landsat8和我國高分一號(GF-1)衛(wèi)星影像的光譜和空間信息，并且加入了LiDAR點云數(shù)據(jù)的高程和強度信息，結(jié)合多源特征的優(yōu)勢，精確提取不透水面。

1 數(shù)據(jù)源及其預處理

1.1 數(shù)據(jù)源

本研究多源數(shù)據(jù)選用30 m和15 m空間分辨率Landsat8影像、8 m和2 m空間分辨率GF-1影像以及0.5 m空間分辨率LiDAR點云數(shù)據(jù)。全部數(shù)據(jù)均在2013年10—12月間采集，獲取時間相近。其中，Landsat8影像波段數(shù)較多，光譜信息豐富； GF-1影像空間分辨率較高，空間信息豐富； LiDAR點云數(shù)據(jù)擁有高程和強度信息。

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

首先進行影像配準，將不同傳感器獲取的不同空間分辨率衛(wèi)星影像配準到統(tǒng)一的坐標系下，并且對衛(wèi)星影像進行輻射糾正，建立同時相、多尺度、多源數(shù)據(jù)金字塔。其中，由于多源影像間空間分辨率跨度較大，直接配準會導致操作不便且精度較低，故采取根據(jù)空間分辨率大小，逐級配準的方法。以ENVI軟件中同地區(qū)的數(shù)字高程模型(digital elevation model，DEM)為基準，先對全色和多光譜GF-1影像分別進行無控制點的正射校正，再逐次將全色和多光譜Landsat8影像配準到正射校正后的GF-1影像坐標系。

2 提取方法

為了精確提取不透水面，采用Grabcut融合多源衛(wèi)星影像的光譜和空間信息，并輔以LiDAR點云數(shù)據(jù)的高程和強度信息。研究方案如圖1所示。

圖1 提取方法流程Fig.1 Flowchart of extracted method

首先進行影像預處理，統(tǒng)一坐標系； 然后分別從衛(wèi)星影像和點云數(shù)據(jù)中提取光譜特征和道路信息； 再將提取的光譜特征使用模糊C均值聚類(fuzzy C-means clustering，F(xiàn)CM)方法進行聚類，以分類粗結(jié)果作為分類初值； 最后構(gòu)建基于Grabcut的分類模型，建立能量函數(shù)，融合多源特征，迭代優(yōu)化初值。

2.1 多源特征

2.1.1 光譜特征

由于不透水面普遍在熱紅外波段輻射率很高，在近紅外反射率很低[14]，因此衛(wèi)星影像可以在不透水面提取方面發(fā)揮其光譜特征方面的較大優(yōu)勢。本文將歸一化差值植被指數(shù)(normalized differential vegetation index，NDVI)、增強型植被指數(shù)(enhanced vegetation index，EVI)和改進的歸一化差值水體指數(shù)(modified normalized differential water index，MNDWI)選作光譜特征參數(shù)，用于Grabcut的輸入特征和初值獲取。

1)NDVI。通過表現(xiàn)植被強吸收能力的紅光波段和表示植被生物量的近紅外波段的組合計算，可以使感興趣地物信息得到增強[15]，其公式為

(1)

式中：NIR表示近紅外波段的反射率值；RED表示紅光波段的反射率值。

2)EVI。在NDVI的基礎上改進得到EVI，其公式為

(2)

式中：BLUE表示藍光波段的反射率值；L=1，為土壤調(diào)節(jié)參數(shù)； 參數(shù)C1=6.0，C2=7.5，用來描寫通過BLUE修正大氣對RED的影響。

3)MNDWI。由于土壤、水體和沙土、部分水泥瀝青具有類似的光譜特征，單利用NDVI無法準確區(qū)分出不透水面，有必要引入水體指數(shù)。水體在綠光波段中比紅外波段具有更高反射率，MNDVI可以準確提取出水體[16]，其公式為

(3)

式中：GREEN表示綠光波段的反射率值；MIR表示中紅外波段的反射率值。Landsat8影像第6波段SWIR1對應中紅外波段。

2.1.2 點云特征

由于光譜相似性，單從衛(wèi)星影像中難以準確區(qū)分出人工道路和裸地。通過引入LiDAR數(shù)據(jù)，得到較準確的人工道路和停車場等不透水面，彌補衛(wèi)星影像的缺陷。

LiDAR系統(tǒng)可以快速獲取具有地表三維坐標及反射強度等信息的密集點云，對于植被、建筑物和道路等地物的精確提取具有重要意義。然而僅依靠高程信息進行分類的可靠性有限，并且強度信息的噪聲大，影響因素多，還難以真實反映地物反射率信息[17]。因此，可以結(jié)合高程和強度信息，將點云信息作為輔助信息參與到不透水面分類過程中。

1)高程信息。為了精確提取道路，首先需要利用高程信息通過點云濾波區(qū)分地面點和非地面點； 再采用形態(tài)學濾波，依次對移動窗口中的點進行腐蝕和膨脹，搜索窗口內(nèi)的最低點，逐步構(gòu)建DEM。此算法對于城市地區(qū)濾波很有效，且計算簡單。

2)強度信息。LiDAR強度信息是接收到的目標回波信號形成的圖像。瀝青、混凝土的回波強度明顯區(qū)別于植被等自然地物，且在回波特點上，道路點云都為末次回波[18]。對于強度值，有學者采取先對樣本標定再提取的方法，也有采取聚類方法提取點云數(shù)據(jù)中的道路[17]。基于以上理論，本文首先去除多次回波點，排除一部分非道路點； 然后采取FCM方法利用強度值將道路和植被、裸地分開。同時，在實驗過程中，發(fā)現(xiàn)道路和建筑物質(zhì)地均一，回波強度在某一小范圍內(nèi)，而裸地里包含有沙石、土壤和草地等，回波強度起伏較大，因此將回波強度的方差值同強度值一起作為聚類的特征，增強道路信息。

2.2 多源數(shù)據(jù)融合分類模型

2.2.1 圖割與Grabcut

圖割(Graph cut)是一種基于圖論的圖像分割方法，其能量最小化框架由Boykov等提出[19]，近年來在計算機視覺領域中被廣泛應用。該方法可以很好地將圖像的區(qū)域和邊界信息相結(jié)合，建立能量模型，通過最大流最小割算法，從而得到能量最小的分割結(jié)果。

Grabcut是在此研究基礎上進行了改進，將交互方式從種子點選取變成框選目標，并且引入高斯混合模型進行建模，通過不斷分割估計和模型參數(shù)學習進行迭代，實現(xiàn)分割[20]。其能量函數(shù)為

E(α,κ,θ,z)=U(α,κ,θ,z)+V(α,z),

(4)

式中：α為像素標記值；κ為像素的高斯分量；θ為像素屬于前景/背景的概率；z為圖像像素；U(α,κ,θ,z)為數(shù)據(jù)項；V(α,z)為平滑項。

數(shù)據(jù)項表示圖模型上所有節(jié)點的權(quán)值，由高斯混合模型計算得到，即

U(α,κ,θ,z)=∑D(αn,κn,θ,zn),

(5)

式中n表示第n個像素。該項也表示像素被標記為前景或背景的懲罰，即某個像素屬于前景或背景概率的負對數(shù)，代入高斯混合模型可得

D(αn,κn,θ,zn)=-lnp(zn|αn,κn,θ)-ln π(αn,κn) ，

(6)

(7)

式中： det用來計算行列式的值；μ為每個高斯分量的均值向量。

平滑項表示圖模型上所有邊的權(quán)值，由相鄰節(jié)點間的歐式距離計算得到，體現(xiàn)了像素間不連續(xù)的懲罰，即

(8)

式中： ?為常數(shù)，用來調(diào)節(jié)數(shù)據(jù)項和平滑項的比例；i和j表示圖像范圍C內(nèi)的相鄰像素； 參數(shù)β由圖像對比度決定，用來調(diào)節(jié)相鄰像素間的差異。

通過最大流最小割算法，計算能量函數(shù)最小狀態(tài)，更新高斯混合模型的參數(shù)，重新計算新的能量函數(shù)，如此反復迭代得到最優(yōu)分割結(jié)果。

2.2.2 基于改進的Grabcut不透水面提取模型

將Grabcut和多源遙感數(shù)據(jù)結(jié)合，從能量函數(shù)項、初值設置和特征空間3方面進行改進，融合多源特征實現(xiàn)不透水面提取。

1)能量函數(shù)項。由于點云數(shù)據(jù)的不穩(wěn)定性，只能作為輔助條件參與分類，因此在原有函數(shù)基礎上，增加了一項道路信息的約束項L(α)，改進后的能量函數(shù)為

E(α,κ,θ,z)=U(α,κ,θ,z)+V(α,z)+L(α),

(9)

式中： 數(shù)據(jù)項U表示由光譜特征建模得到每個像素隸屬于對應類別的權(quán)值； 平滑項V表示所有相鄰像素間特征距離的加權(quán)累加；L(α)為道路信息的約束項，表示由點云數(shù)據(jù)判斷像素隸屬于對應類別的權(quán)值，即

L(α)=γ∑[ln≠αn],

(10)

式中：γ為新增點云數(shù)據(jù)項的權(quán)值，由于點云數(shù)據(jù)強度信息的不穩(wěn)定性，僅將其作為一種軟約束，可根據(jù)實際情況，通過改變權(quán)值大小來調(diào)節(jié)點云信息的影響力；ln表示點云數(shù)據(jù)中的像素類別標記，當其和當前標記值αn不一致時取1，否則取0。

2)初值設置。原始Grabcut中，將人工框選區(qū)外的像素標記為絕對背景初值，框內(nèi)標記為可能前景初值，通過迭代一步步在可能的前景中篩選出絕對前景。本文將不透水面和透水面分別作為分類的前景和背景，需要處理的是覆蓋面積大、連通區(qū)域細碎的衛(wèi)星影像，人工標記種子點工作量大，需盡可能減少人工干預。所以首先對多源衛(wèi)星影像進行了FCM聚類，將得到的粗分類結(jié)果作為Grabcut的初值，即將聚類得到的不透水面和透水面分別標記為可能前景和背景初值。

3)特征空間。一般的圖割基于顏色空間(灰度或RGB)對目標和背景建模，經(jīng)典的模型有灰度直方圖、高斯混合模型等。但是對于多源多光譜影像而言，每個波段的重要性不一樣，顏色信息冗余，所以需要先進行篩選組合，將計算得到的一系列突出表征特定地物的指數(shù)作為輸入特征，進行建模。本文采用NDVI，EVI和MNDWI作為輸入特征，采用高斯混合模型建模，并結(jié)合由鄰接關(guān)系得到的空間特征和點云數(shù)據(jù)中提取的道路信息，共同構(gòu)建能量函數(shù)。

通過以上3方面的改進，本文將改進的Grabcut引入不透水面提取中，建立了新的能量函數(shù)，融合植被特征、水系特征、道路特征以及各相鄰點間的光譜和空間特征，使分類過程更有針對性，并且很大程度上減少了人工干預，有利于分類精度的改善和不透水面的準確提取。

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

在Visual Studio2013軟件平臺上，通過編程實現(xiàn)研究區(qū)不透水面的提取。所需的多源特征如圖2所示。

(a) NDVI (b) EVI (c) MNDWI (d) LiDAR道路信息

圖2多源特征

Fig.2Multi-sourcefeatures

首先，對多源衛(wèi)星影像計算指數(shù)值，GF-1具有紅光、綠光、藍光和紅外4個波段，且空間分辨率較高，故采用GF-1計算NDVI和EVI值，如圖2(a)和(b)所示； 由于Landsat8具有11個波段，且具有對應于中紅外的SWIR1波段，故采用Landsat8計算MNDWI，如圖2(c)所示。其次，結(jié)合高程和強度信息從LiDAR數(shù)據(jù)中提取出人工道路，如圖2(d)所示。然后，根據(jù)提取的光譜特征，采用FCM聚類得到初值。最后，構(gòu)建圖割能量函數(shù)，迭代得到能量最小的分類結(jié)果。

3.1 實驗結(jié)果

實驗結(jié)果如圖3所示。

(a) GF-1 B3(R)，B2(G)，B1(B)假彩色合成影像 (b) 聚類初值(c) 最終分類結(jié)果

圖3分類結(jié)果

Fig.3Classificationresults

從圖3(c)中可以看出，經(jīng)過Grabcut分割后，由于空間特征的影響，離散的噪聲點都被消除了； 道路也由于點云信息的加入被較好地提取出； 水系、植被和建筑物也均被較好提取出； 人造田徑場通常是在水泥地面上鋪蓋草皮，應該被分為不透水面，分類正確； 然而裸地和人造道路間存在部分誤分，原因是部分裸露地表含水量低、沙石含量較高、有人工挖掘碾壓痕跡，光譜表現(xiàn)和不透水面相似，且表面平滑無法根據(jù)高程信息計算粗糙度來剔除。局部細節(jié)分類結(jié)果如圖4所示。

(a) 局部影像1(b) 局部影像2 (c) 局部影像3(d) 局部影像4 (e) 局部影像5

(f) 局部結(jié)果1(g) 局部結(jié)果2 (h) 局部結(jié)果3(i) 局部結(jié)果4 (j) 局部結(jié)果5

圖4局部細節(jié)分類結(jié)果

Fig.4Detailsofclassificationresults

從圖4(f)和(g)可以看出植被和建筑群間區(qū)分清晰； (h)中道路和植被分類正確； (i)中左下角為水系，右上角為植被，中間為獨立建筑物，均能被較好地區(qū)分； (j)右上角的植被被正確提取出，但左半部分和左下角的裸地均沒有提取出，被誤分為不透水面。說明本文方法在裸地和不透水面的區(qū)分上還有待改進，對于其他地物均能正確區(qū)分。

3.2 精度評定

為了定量檢驗分類精度，隨機選取了均勻分布在影像中的5 136個測試樣本，每個樣本均為5像素×5像素大小，將各個樣本中的每個點，通過人工解譯賦予屬性作為真值，其中不透水面點有71 550個像素，透水面點有56 850個像素。在ENVI軟件上完成真值勾畫。

為了證明本文方法能優(yōu)化聚類結(jié)果，且選取的多源數(shù)據(jù)均對分類結(jié)果有利，分別對本文方法結(jié)果、聚類所得初值的分類結(jié)果、不加入LiDAR數(shù)據(jù)的分類結(jié)果、不加入Landsat8數(shù)據(jù)的分類結(jié)果進行精度評定。

精度評價結(jié)果如表1所示。

表1不同數(shù)據(jù)源的分類精度

Tab.1Classificationaccuracyofdifferentsourcedata

數(shù)據(jù)源生產(chǎn)者精度/%用戶精度/%不透水面透水面不透水面透水面總體精度/%Kappa系數(shù)本文方法93.4786.2099.9399.9790.250.820聚類所得初值88.4187.4289.8485.7087.970.783不加入LiDAR數(shù)據(jù)89.3886.2989.1386.5988.010.783不加入Landsat8數(shù)據(jù)90.9381.1885.8887.6786.610.760

由表1中分析，僅采用聚類所得初值、不加入點云數(shù)據(jù)和不加入Landsat8數(shù)據(jù)時，分類精度均比本文方法低，主要體現(xiàn)在不透水面的漏分上。結(jié)合分類結(jié)果影像來分析，原因其一是缺乏鄰接關(guān)系干預導致離散噪聲未被剔除； 其二是缺乏點云信息的干預而導致道路漏分； 其三是Landsat8的中紅外波段對水系敏感，其缺失將使水系被誤分成不透水面。由此證明了本文將多源數(shù)據(jù)在Grabcut框架下融合分類對于不透水面提取有一定意義，精度得到提高。

為對比分析本文和其他不透水面提取方法的精度，本文使用ENVI5.1軟件里的最大似然監(jiān)督分類、線性光譜混合分解和決策樹分類3種方法分別進行分類實驗，精度評價如表2所示。根據(jù)表2結(jié)果可知，本文方法的總精度達到90.25%，Kappa系數(shù)為0.820，均優(yōu)于其他3種傳統(tǒng)遙感方法，在道路和水系提取上更準確。同時，相比于監(jiān)督分類需要采集感興趣區(qū)域作為訓練樣本、線性光譜混合分解需要人工選取端元波譜、決策樹分類需要進行屬性選擇度量，本文方法在人工干預方面大大減少了工作量，具有更高的適應性和穩(wěn)定性。

表2本文方法與傳統(tǒng)遙感方法的分類精度評價

Tab.2Classificationaccuracyoftheproposedmethodandthetraditionalremotesensingmethods

方法生產(chǎn)者精度/%用戶精度/%不透水面透水面不透水面透水面總體精度/%Kappa系數(shù)本文方法93.47 86.20 99.93 99.97 90.25 0.820 最大似然監(jiān)督分類77.47 91.48 91.97 76.34 83.68 0.719 線性光譜混合分解94.20 72.43 81.13 90.84 84.56 0.727 決策樹分類89.22 82.66 86.62 85.90 86.31 0.756

4 結(jié)論

本文結(jié)合多光譜影像Landsat8、高空間分辨率衛(wèi)星影像GF-1和LiDAR點云數(shù)據(jù)的優(yōu)勢，構(gòu)建圖模型，將Grabcut引入不透水面提取中。結(jié)論如下：

1)基于Grabcut的多源數(shù)據(jù)融合提取方法，通過能量優(yōu)化框架，實現(xiàn)了光譜信息、空間信息、高度和強度信息等多源特征融合。根據(jù)多源數(shù)據(jù)特點設計特征，一定程度上解決了裸地和不透水面相互混淆的問題，相比于單一數(shù)據(jù)提取精度得到明顯提高。

2)本文方法能夠?qū)ο袼刂g的空間關(guān)系和局部的相似性信息進行很好的表達，相較于傳統(tǒng)的像素分類方法，有利于獲得噪聲更少的分類結(jié)果。

3)與傳統(tǒng)提取方法相比，本文方法大大減少了工作量，具有更高的適應性和穩(wěn)定性。

基于Grabcut的多源數(shù)據(jù)融合提取方法在以下方面還有待進一步深入研究： 第一，可引入更多突出表現(xiàn)特定地物的波段或指數(shù)，提高分類精度； 第二，采用的LiDAR數(shù)據(jù)空間分辨率僅0.5 m，而若使用空間分辨率更高的密集點云，可以通過計算粗糙度來區(qū)分自然裸地和平坦道路，在一定程度上減少誤分； 第三，需對更多區(qū)域進行實驗，檢驗方法的普適性。