面向電子地圖更新的水體自動提取方法研究

何正國，毛海亞，鐘家暉，樊惠萍，苗立新，周淑芳，劉 烽，陳 婷

(1. 廣州市城市規劃自動化中心，廣東 廣州 510030；2. 二十一世紀空間技術應用股份有限公司，北京 100096)

水體提取的難點主要由光譜多樣性、形態多樣性和季節變化因素引起。水體包括面狀、線狀等不同的形態，同時隨年份、季節波動，造成水體在不同區域和不同時相影像上形成迥異的特征，這都對特征選取、信息提取和質量控制等提出了較高要求，使得陸表水體信息提取技術更加復雜[1-4]。本文選取廣州市作為研究區，基于高分辨率的北京二號遙感數據和面向對象的技術流程，提出一套基于全局最優閾值法的水體自動提取方法，旨在提高電子地圖更新的自動化程度，縮短更新周期。

1 研究區與數據

1.1 研究區概況

廣州市是我國的國家中心城市、廣東省省會城市和珠三角城市群中心城市。國家“一帶一路”建設已經進入縱深發展階段，粵港澳大灣區城市群發展亦寫入了2017年國務院《政府工作報告》，廣州作為國家中心城市以及粵港澳大灣區核心門戶城市之一，城市化進程越來越快。其空間數據的更新應與之相符，廣州市電子地圖數據作為典型的空間數據之一，應與城市發展同步，實現高頻度的地圖更新。

廣州市目前正在開展“四標四實”工作，每天有大量的網格員采集數據，因此有必要對傳統人工更新電子地圖的方式進行革新，研究電子地圖要素信息的自動提取方法。選取廣州市南部作為研究區，見圖1。該區域內水體類型豐富，既包括了天然的大型河流、湖泊，又包括了人工修筑的水庫、坑塘，水體周邊的地物環境也包含了建設用地、植被、裸地等不同類型，可以充分檢驗水體自動提取方法在不同地物構成下的適用性和有效性。

1.2 數據處理

由于城市建設步伐的加快，城市基礎地理信息的變化也越來越快，公眾對電子地圖的現勢性要求越來越高，常規的電子地圖更新時效已經無法滿足用戶的需求。遙感技術對地物信息的快速、大面積獲取能力，使其在電子地圖的數據獲取領域受到越來越多的重視。研究中采用了亞米級的北京二號衛星影像數據，用于水體的自動提取實驗。

圖1 廣州研究區位置示意圖

北京二號遙感衛星星座于2015年7月11日在印度的斯里赫里戈達島(Sriharikota)衛星發射場成功發射。星座由3顆光學遙感衛星組成，均搭載有多光譜和全色成像儀，對全球任一地點的重訪周期為1～2 d，在時間尺度上為水體的高密度動態監測提供了可能。每顆衛星在軌提供幅寬約24 km、0.8 m分辨率全色(450～650 nm)和3.2 m分辨率藍(440～510nm)、綠(510～590 nm)、紅(600～670 nm)、近紅外(760～910 nm)多光譜圖像，波段設置滿足遙感水體信息提取的一般要求。

使用1景北京二號影像，成像時間為2018年1月12日，影像預處理主要包括系統幾何校正、正射校正、波段融合等。結合地形圖，通過選取地面同名控制點的方法對影像進行幾何校正，糾正后各點誤差在一個像元內。通過控制點采集及其優化、區域網平差模型構建及其解算，完成全色和多光譜影像的正射校正處理。采用Gram-Schmidt變換方法對校正后的全色、多光譜數據做融合處理，融合后的影像既保留了多光譜信息，又在原有基礎上提高了其空間分辨率，得到亞米級的影像數據。

2 技術方法

2.1 面向對象分割

研究采用面向對象的技術流程，通過多尺度分割(multi-resolution segmentation)算法，使影像對象成為特征信息的載體[5]。該算法根據像元間的異質性，通過定義一定的閾值將影像分割成相對同質的對象，并以對象作為基本的分類單元。其綜合利用了影像的光譜特征、幾何特征，避免了對象破碎，提高了影像的分類精度[6-7]。

多尺度分割算法是一種區域生長算法。它基于影像對象的特征，針對不同土地覆蓋類型的尺度效應，解決產生斑塊與真實地物的邊界擬合問題。分割對象的平均異質性最小化與像素的平均異質性最小化是影像分割的標準。當異質性大于等于給定的尺度時，區域就會停止生長，區域內像元將合并形成對象。異質性的計算公式為：

f=wcolor·[∑cwc·δc]+(1-wcolor)·

(1)

式中：wcolor為光譜異質性權重；wc為數據層權重；δc為數據層上光譜標準差；wcpt為緊致度異質性；l為對象長度(像素個數表示)；b為目標多邊形最短邊(像素個數表示)；n為像元個數。

對于多尺度分割，分割尺度的確定至關重要，是后續分類精度的重要保證。分割尺度越大，誤分割圖斑的數量會減少，但誤分割的總面積會增大。為綜合考慮兩種因素，對誤分割圖斑個數和誤分割圖斑面積進行歸一化處理，建立誤分割指數模型，以期得到最佳分割尺度使得誤分割指數達到最小值。誤分割指數計算公式如式(2)所示。

ESI=n2+r2

(2)

式中：n為歸一化后誤分割圖斑數；r為歸一化后誤分割圖斑面積；ESI為誤分割指數。

2.2 特征提取規則集

基于表觀要素中的光譜特征和紋理特征等，構建水體異于其他地物類型的特征集。相較于傳統的分類技術而言，面向對象的信息提取方法是基于對象的操作，能夠充分利用提取對象的空間、質地和光譜特征，集合臨近像元進行分割和分類，最終實現分類結果的矢量輸出。

(1)波段均值(Mean_X)。是指每個對象內第X波段所有像元灰度值的平均。其中，紅波段(Mean_R)與藍波段(Mean_B)的灰度值均值，可以分別剔除大部分建設用地與藍色屋頂的工廠板房。此外，部分水體表面被水生植被所覆蓋，歸一化植被指數(NDVI)取值范圍與地表植被很接近，基于NDVI的提取混分嚴重，采用近紅波段均值(Mean_NIR)可以擴大這部分水體與植被的差異，減少錯提。亮度(Brightness)也屬于波段均值的一種，是通過所有影像層計算出來的，該特征可以區分水體與周邊地物，特別是與建筑物陰影的差異。

(2)歸一化水體指數。水體因對入射能量(太陽光)具有強吸收性，所以在大部分遙感傳感器的波長范圍內，總體上呈現較弱的反射。Mcfeeters提出了歸一化差異水體指數(NDWI)，旨在利用比值運算快速提取水體信息[8-9]。由于水體的反射從可見光到中紅外波段逐漸減弱，在近紅外和中紅外波長范圍內吸收性最強，幾乎無反射，因此基于可見光波段和近紅外波段的反差構成的NDWI可以突出影像中的水體信息。此外，由于植被在近紅外波段的反射率一般最強，因此采用綠光波段與近紅外波段的比值可以最大程度地抑制植被的信息，從而達到突出水體信息的目的。針對廣州城市水體進行分析，多數水體對象的NDWI取值在0.5以上，可以比較明顯地區別于植被與建設用地。

(3)歸一化植被指數。比值型指數創建的基本原理就是在多光譜波段內，尋找出所要研究地類的最強反射波段和最弱反射波段，將強者置于分子，弱者置于分母。通過比值運算，進一步擴大二者的差距，使感興趣的地物在所生成的指數影像上得到最大的亮度增強，而其他背景地物則受到普遍的抑制，從而達到突出感興趣地物的目的。基于此原理，歸一化植被指數是反映植物生長狀況及植被覆蓋度的最佳指示因子。研究區內水體對象的NDVI指數普遍處于-0.3以下，區別于大多數植被。

(4)灰度共生矩陣(GLCM)。灰度共生矩陣原理是：如果有n×n個灰度大小的影像，則該影像能夠轉化為n×n階的灰度共生矩陣，在影像平面(x,y)中的某處(i,j)對應一個向量D=(Δx,Δy)。對于隨意一個這樣的向量和影像，可以計算由D向量偏移的成對出現的亮度值之間的聯合概率密度P(i,j)。灰度共生矩陣的紋理信息都包含在移動方向為θ=0°、θ=45°、θ=90°和θ=135°的亮度值空間依賴矩陣中。能夠用于提取紋理特征的紋理測度有很多統計值，而其中最常用的幾個二階概率紋理特征統計值主要是：均值(Mean)、方差(Variance)、同質性(Homogeneity)、對比度(Contrast)、非相似度(Dissimilarity)、角二階矩(Angular Second Moment, ASM),以及熵(Entropy)和相關性(Correlation)，其中的每一個GLCM紋理測度都可以單獨使用。

值得注意的是，在人工確定閾值的過程中，各單波段均值、亮度值及NDWI與NDVI的設定要依區域位置和時相的情況而定，不同區域、時相的影像差別較大。特別是在設定NDWI取值范圍時，要根據研究區實際情況結合參考數據反復調試提取閾值的大小。這些雖然會達到理想的提取效果，但耗費較多人力與時間，在業務化實施中并非是最佳選擇。

2.3 全局最優閾值算法

水體包括面狀、線狀等不同的形態，其提取難點主要由光譜多樣性、形態多樣性和季節變化因素引起。在采用水體指數、植被指數等進行水體信息提取過程中，指標閾值的確定是非常關鍵的。為了提高生產效率，同時針對本項目開展的業務化應用需求，采取了自動計算閾值的方法進行水體提取。

水體通常具有較大面積，呈現面狀分布的特點，且內部特征相對均質，適用于根據水體指數、植被指數等設定閾值對影像進行前景與背景的劃分。采用全局最優閾值的方法，其中心思想是基于迭代的方式將灰度圖像用兩個或兩個以上正態分布的概率密度函數做近似表示，每次都取其中最顯著的波峰來劃分區域，然后依據各個區域的平均值選擇合適閾值，重復過程直到該值收斂。方法實現如下。

遍歷圖像中各點讀取像素，獲得分割圖像的最小灰度值T1以及最大灰度值T2，設定初始閾值T(0)，T(0)=(T1+T2)/2；根據閾值T(0)進行圖像分割，將圖像分為前景和背景兩個區域。再次對分割后的圖像進行遍歷，計算出此時的前景T1和背景的平均灰度值T2，再次平均求出閾值T(1)。用gray(i,j)表示待處理圖像像素(i,j)點的灰度值，N1表示圖像前景像素的個數，N2表示圖像背景像素的個數，計算如式(3)～式(5)。

(3)

(4)

T(1)=(T1+T2)/2

(5)

重復上述運算，直到T(k)=T(k+1)，迭代完成，此時用T(k)作為最終閾值對圖像進行分割。

全局最優閾值法的水體提取技術首先根據影像的NDWI指數將其分為前景和背景值，然后利用全局最優閾值法自動獲取NDWI、NDVI、Brightness指數的閾值，分別提取水體和剔除混分的植被、建設用地等要素。同時，考慮到在高分辨率遙感影像中，建筑物陰影是較為明顯的干擾因素，且往往被錯提為水體，因而在規則集中增加了紋理特征的灰度共生矩陣(GLCM)，用于進一步剔除建筑物陰影等干擾。經過對比實驗，特別采用了紋理同質性指標(GLCM Homogeneity)作為對影像分布平滑性的測度，其計算公式如式(6)所示，勻質區域具有較高的特征值，反之則較低。

hom=∑∑p(i,j)/[1+(i-j)2]

(6)

3 結果與討論

3.1 水體對象的多尺度分割

主要基于面向對象的方法進行水體信息提取，該方法根據像元間的異質性，通過定義一定的閾值將影像分割成相對同質的對象，并以對象作為基本的分類單元。其綜合利用了影像的光譜特征、幾何特征，避免了對象破碎，提高了影像的分類精度[10]。

先對研究區內的北京二號影像進行多尺度分割，分割后的特征值反映了每個對象范圍內所有像元的平均值。當分割尺度過大時，水體邊界與地面實體相差較大，邊界不準確，對象內其他地物“雜質”過多、同質性差，為后期信息提取帶來困難。而分割尺度過小時，雖然地物邊界清晰，對象同質性較好，能夠準確地反映地物空間分布狀況，但會使分割過程及后期提取效率低下，混分的細碎斑塊不易剔除，也難以形成規整的水體邊界，不利于提取成果通過GIS方法等進行分類后處理(見圖2)。經過反復對比，確定分割尺度為50(形狀參數0.1，緊致度參數0.5)，以盡量將光譜、紋理特征接近的像元分割確定為同一對象，同時避免斑塊過于細碎。經過多尺度分割后的影像對象產生了新的特征信息，相較于單個像元無疑含有更多的信息量。

本研究基于多尺度分割方法充分利用了多個尺度的信息，對分割結果進行約束優化，是目前一種較好的信息提取方法，應用也較廣。多尺度分割是一種自下而上，采用異質性最小的區域合并策略的影像分割方法，分割時基于局部最佳相互適配原則，將每個像元逐步合并成較大的多邊形對象。通過合理設置分割參數，使分割后的對象盡可能體現出研究區范圍內不同構成要素間的差異，以便后續對樣本特征參數的分析與選擇。對水體進行多尺度分割時的參數選擇應結合影像空間分辨率與地表覆蓋的具體情況而定，分割結果既能滿足同類要素的均一性，又能體現不同要素間的差異，并達到分割精度與效率的統一。

3.2 基于規則集的提取結果

獲取了2018年覆蓋廣州研究區的北京二號融合影像，包含可見光、近紅外共4個波段，空間分辨率0.8 m。在該區域內開展了水體(主要為河流、湖泊、坑塘等)的自動提取實驗，確認方法及技術流程的可行性。共采用了NDWI、NDVI、Brightness、GLCM Homogeneity (band 1)與GLCM Homogeneity (band 2)等特征，均基于全局最優閾值法進行計算，確定各提取規則的特征值，見表1。首先充分利用NDWI指數在水體提取中的優勢，基于迭代的方式將NDWI灰度圖像的特征直方圖用正態分布的概率密度函數做近似表示，取其中最顯著的波峰將圖像劃分為前景和背景兩個區域。其中，前景即為初步提取的水體，背景為其他地物類型。隨后，對前景中的水體做進一步優化，不斷剔除錯提的其他地物。需要再次對分割后的影像進行遍歷，分別依據NDWI、NDVI、Brightness、GLCM Homogeneity等每一類特征，計算出前景和背景的平均值，并由此獲取各特征的合適閾值，重復過程直到該值收斂。最后獲取的特征閾值，即是經過不斷迭代優化后的取值。

注：(a)為北京二號融合影像，(b)、(c)、(d)所對應的分割尺度分別為30、50、70圖2 不同分割尺度下的影像對象變化

表1 全局最優閾值法確定的提取規則閾值

水體的初步提取結果如圖3所示，研究區范圍內除了大型河流外，還存在較多的小型水體與坑塘，且零星分布。提取結果中的水體邊界相對規整，后處理工作量小。除主要河流外，坑塘等小型水體往往受到季節、氣象及人為因素的影響，年際變化較大，這也是南方水體普遍存在的明顯特征。因而，使用自動提取方法可以較為準確、快速地獲取其分布情況，有利于減少對此類小型水體的漏提。主要的后處理操作為排查與歸并小圖斑，剔除偽變化，同時也可以應對數據冗余。

圖3 廣州研究區水體自動提取結果

3.3 結果后處理與精度評價

在完成初步的信息提取后，還需要進一步開展結果的后處理工作，并主要針對于水體中面積較小的圖斑對象，以有效減少錯提和漏提這兩類現象，見圖4。首先將水體圖斑進行合并，使空間位置上接壤的對象合并為一個較大的對象。這樣一方面減小了后續操作的圖斑數據量，提高運行效率；另一方面也便于通過設定面積閾值，剔除掉面積較小的孤立圖斑。基于北京二號影像空間分辨率及研究區內水體分布的實際情況，確定面積閾值為4000像素，將小于該值的水體對象賦類為其他，可以對錯提的陰影等要素進行剔除，見圖4中(a)、(d)、(g)。同時，由于研究區內水運發達，水面船只數量繁多，見圖4中(b)、(e)、(h)，農田水利設施廣泛分布,見圖4中(c)、(f)、(i)，考慮到水面船只、水利設施等對水體提取的干擾，將此類面積小于17 000像素且完全被水體圖斑所包含的對象，也重新歸類于水體，保證了水體圖斑的完整，減少了漏提。

對水體自動提取的精度與效率進行分析，在研究區內生成300個隨機點，通過人工判定方式對自動提取結果開展精度評價。經過混淆矩陣的計算，水體的生產者精度為64%，用戶精度則達到了96%，Kappa系數為0.73。整個研究區面積達到了近300 km2，水體總面積約為98 km2，水體圖斑總數超過了2800個，全流程的對象分割、信息提取可以在1 h之內完成。其中,除多尺度分割與GLCM紋理特征計算較為耗時以外，基于全局最優閾值法的各特征值獲取均在3 min內完成，具有較高的計算效率。從最終提取水體對象的數量、面積及準確度而言，自動提取方法可以極大地提高工作效率，有利于快速獲取水體要素分布情況，而人工參與只需要安排在后續的成果精細化修訂環節即可。

注：(a)、(b)、(c)為北京二號融合影像，(d)、(e)、(f)為水體初步提取結果，(g)、(h)、(i)為后處理結果圖4 針對錯堤和漏堤的水體提取結果后處理

值得注意的是，基于高空間分辨率影像開展水體提取，在城市范圍內往往存在建筑物等地物陰影的干擾問題。在本研究區內，部分高層建筑物、植被及山體的陰影會對相鄰地物造成明顯遮擋。如圖4(a)所示，陰影遮擋下的地物光譜信息受到抑制，特別是植被、裸地等亮度、紋理等特征更接近于水體，導致在多尺度分割時就難以識別其準確邊界，因而產生錯提。如果在此基礎上開展進一步的全覆蓋地物類型提取，則會使誤差繼續累加。針對這一問題，可以考慮在研究區內提前圈定盡可能大的水體矢量范圍，作為水體提取的目標區域，以便排除建成區內易混淆的高大建筑物陰影等。規則集運行中只提取目標區域內的水體，這樣可以比較有效地規避建筑物及植被陰影的影響，減少水體錯提。

4 結 論

本文以廣州市部分水體為研究對象，基于北京二號高分辨率遙感影像，采用面向對象的技術流程，通過全局最優閾值算法實現了水體的自動化信息提取。研究表明，基于全自動的閾值確定方法，在無人工參與情況下，可以在較大面積的區域內實現水體要素的快速提取，應用效果較好。研究區內水體類型復雜，表觀特征差異較大，針對電子地圖更新的數據需求，通過面向對象的方法，提取的水體圖斑邊界規整，能夠滿足圖4的精度要求。水體的提取作為電子地圖數據更新的第一步，也為后續植被、建筑物、道路等信息提取與更新提供了重要參考，和傳統的人工矢量化相比，有效提高了工作效率。提取結果的總體精度可以達到96%，Kappa系數為0.73。