一種衛星遙感影像的城市水體提取方法與流程

本發明屬于水體遙感影像技術領域，主要涉及一種衛星遙感影像的城市水體提取方法。

背景技術：

城市是人類社會高度發展的體現，城市水體作為城市生態系統中重要的因素，在維持城市生態系統穩定性上具有十分重要的作用。城市水體的改變將會對生活產生巨大的變化，可能會引發一些災害，如城市水源滯留、水資源短缺，甚至引發一些和人類健康生活相關的疾病。因此，了解和掌握城市水體分布以及水域面積的變化已經成為了人們日益關注的焦點。

近年來，隨著遙感技術的應用與發展，遙感影像在自然資源調查、動態監測、自然地表水源規劃等方面發揮著著越來越重要的角色，利用遙感技術進行地表的監測也得到了越來越多的科研工作者的關注；遙感影像可以以一個不同的視角去觀察地球表面的地物，實時的監測地表的變化。在水體提取技術中，通過遙感數據及時準確地獲取城市水體信息成為目前主流的水體提取方式。到目前為止，諸多學者提出了大量遙感影像的水體提取方法，但大部分都是基于中低分辨率遙感影像，由于影像分辨率較低，面積較小的水體未能有效的提取；尤其是對于城市地區水體的提取，城市地區水體面積大小參差不齊，存在諸多小面積人工湖泊和細小的河流。因此，在城市水體提取中應更多的采用高分辨率遙感影像，以資源3號衛星為例，資源3號衛星遙感影像有著5.8m分辨率和5.2km的寬幅，它為城市水體提取提供了理想的多光譜影像數據，資源3號具體參數如下表1：

表1資源三號衛星影像參數

水體是遙感影像中很常見的地物類別，也是重要的基礎地理信息之一，其動態信息的快速獲取，對水資源調查、水利規劃、環境監測與保護等事業都有著十分明顯的實用價值和科學意義。對此，諸多學者很早就對此展開了研究，提出了許多有效的水體信息自動提取模型。大致可以分為4類：(a)單波段和多波段閾值分割法(single-band or multiple-band threshold method)，(b)水體指數法(Water indices)，(c)線性分解模型(linear un-mixing model)，(d)監督與非監督分類方法(supervised or unsupervised classification method)。除此之外，還有一些其他的方法，如：基于數字高程模型的水體提取技術，基于微波遙感(Microwave Remote Sensing)影像的水體提取技術，面向對象(Object Oriented)技術的水體提取技術等。但這些方法并不常用，總得來說，水體指數法由于其模型簡單、方便，且精度較高，在實際中最為常用。

然而隨著遙感影像分辨率的提高，大部分高分辨率遙感影像(如；WorldView-2,IKONOS,RapidEye,and資源3號)沒有像Landsat TM/ETM+/OLI這么多可以利用的波段用于水體的提取，因此，MNDWI(改進的歸一化差異水體指數法)和AWEI(自動水體指數)將無法使用，因為大部分的高分辨率遙感影像只有4個波段(藍色、綠色、紅色、近紅外波段)，缺少MNDWI/AWEI計算所需的短紅外波段(SWIR)。因此，在使用NDWI(歸一化差異水體指數)對高分辨率影像進行水體提取的時候就將會產生一些問題，如陰影無法去除的問題，尤其是城市地區高層建筑物的陰影，在高分辨率影像上表現尤為突出。在高分辨率遙感影像上，城市區域的高層建筑物陰影與水體很難區分，盡管目前有相關學者對這一方面進行了研究，如基于面向對象的技術，通過計算高層建筑物陰影區域的紋理特性來對高層建筑物陰影進行檢測；雖然可以達到預期效果，但由于紋理的描述與計算相對復雜且耗時較長，所以從計算時間上考慮該方法并不是一個理想的陰影檢測方法。也有基于SVM特征訓練進行的陰影檢測，以期達到去除水體檢測中高層建筑物陰影的影響。SVM是一種分類精度較高的方法，但SVM訓練需要花去較多的時間，尤其是當訓練樣本數目較多且樣本特征向量維度較高時。若通過采用高分辨率遙感影像陰影檢測方法(morphological shadow index,MSI)與NDWI相結合的方式對WorldView-2高分辨率影像進行水體提取，以期提高水體檢測精度；盡管該方法原理簡單，但由于該方法以NDWI方法作為水體提取的基礎，其檢測精度并不會很高，尤其是周圍植被茂密的細小面積水域，水體的光譜特性將受到嚴重的污染，水體光譜特征表現極其不穩定，同時城區水體具有懸浮泥沙含量高、水體富營養化嚴重、受各種污染物污染較大等特點，使得城區水體與自然界中未受污染的水體表現不同的光學特性。

技術實現要素：

針對如何從高分辨率衛星遙感影像上，進行城市區域的水體提取問題，特別是有效的區分高層建筑物陰影與水體和提高水體提取的精度問題，本發明提出了一種衛星遙感影像的城市水體提取方法(Automatic urban water extraction method,AUWEM)。

本發明的具體技術方案如下：

一種衛星遙感影像的城市水體提取方法，包括以下步驟：

步驟1：遙感影像數據的預處理，即對遙感影像數據進行正射校正和大氣校正；

步驟2：預處理后的遙感影像數據，包括藍色band1數據、綠色band2數據、紅色band3數據和近紅外band4數據，選取預處理后的遙感影像數據中的藍色band1數據代替歸一化差異水體指數NDWI的計算公式中的綠色band2數據，獲得新的歸一化差異水體指數NNDWI1，新的歸一化差異水體指數NNDWI1的計算公式為：

此計算公式即為NNDWI1指數模型，利用此模型通過閾值分割即得到NNDWI1的閾值分割結果，也即為NNDWI1水體提取結果；

步驟3：對預處理后的遙感影像數據中包括的四個波段數據，即藍色band1數據、綠色band2數據、紅色band3數據和近紅外band4數據進行PCA變換，并將PCA變換后的第一主成分分量Component1替代歸一化差異水體指數NDWI的計算公式中的綠色band2數據，獲得另一個新的歸一化差異水體指數NNDWI2，即：

其中，Component1表示PCA變換的第一主成分分量，此計算公式即為NNDWI2指數模型，利用此模型通過閾值分割即得到NNDWI2的閾值分割結果，也即為NNDWI2的水體提取結果；

步驟4：將步驟2得到的NNDWI1的閾值分割結果和步驟3中得到的NNDWI2的閾值分割結果進行疊加，將得到的結果定義為新的歸一化差異水體指數NNDWI的閾值分割結果，即NNDWI1的閾值分割結果與NNDWI2的閾值分割結果進行疊加，其計算公式為：

NNDWI＝(segmentation_NNDWI1)∪(segmentation_NNDWI2)

式中segmentation_NNDWI1表示NNDWI1的閾值分割結果，segmentation_NNDWI2表示NNDWI2的閾值分割結果，此計算公式即為NNDWI指數模型，利用此指數模型得到NNDWI的水體提取結果；

步驟5：對預處理后的遙感影像數據中的近紅外band4數據進行閾值分割，得到近紅外band4數據的閾值分割結果；

步驟6：對NNDWI的水體提取結果中的大面積水體對象和小面積對象進行分割，NNDWI的水體提取結果中，像素個數大于設定閾值的為大面積水體對象，像素個數小于等于設定閾值的為小面積對象；

步驟7：對步驟6中得到的小面積對象進行數學形態學膨脹處理，得到膨脹后的小面積對象，將步驟5得到的近紅外band4數據的閾值分割結果作為約束條件，即采用膨脹后的小面積對象和近紅外band4數據的閾值分割結果求交集的方式對膨脹后的小面積對象進行約束，約束的數學表達式為：

component2＝(dilate_component)∩(segmentation_band4)

式中，dilate_component表示膨脹后的小面積對象，segmentation_band4表示近紅外band4數據的閾值分割結果，component2表示約束后的小面積對象；

步驟8：對步驟7得到的約束后的小面積對象進行陰影檢測與去除，得到小面積水體對象；

陰影檢測與去除，是指對每個小面積對象中的每個像元進行波譜關系的描述，并判斷該像元是否滿足陰影像元的條件，記錄并統計每個小面積對象中陰影像元的個數，當一個小面積對象中陰影像元所占比例大于閾值T時，把該小面積對象判定為建筑物陰影對象，小于等于閾值T時小面積對象則判定為小面積水體對象，陰影像元所占比例即為小面積對象中陰影像元的個數與該小面積對象中總像元個數的比值，區分小面積對象中小面積水體對象和陰影對象的函數表達式為：

式中，n表示某一小面積對象中總像元個數，m為該小面積對象中陰影像元的個數；

陰影像元的條件，是指滿足陰影像元的波譜大小關系，即滿足以下三個不等式條件：

步驟9：將步驟6中得到的大面積水體對象與步驟8中得到的小面積水體對象進行疊加，即將步驟6中得到的大面積水體對象和步驟8中得到的小面積水體對象求并集，得到衛星遙感影像的城市水體提取結果。

本發明的有益效果如下：

(1)新的歸一化差異水體指數NNDWI1和NNDWI2可以有效的對水體進行初始提取，提高后續水體的提取精度。

(2)在陰影對象獲取方面，為了更加準確的獲取陰影對象，對小面積對象進行膨脹處理；同時為了限制在真實地表陰影區域，對膨脹的對象采用近紅外band4數據的閾值分割結果進行約束。

(3)為避免陰影對象檢測過程中時間消耗，對小面積對象采用基于波譜特性的描述，減少傳統中采用紋理特征進行描述的時間消耗，提高方法的計算效率。

(4)AUWEM方法的分類精度要高于NDWI的分類精度和最大似然法的分類精度。在5個實驗區的實驗結果中，AUWEM平均Kappa系數約為93.0％，NDWI的平均Kappa系數約為86.2％，最大似然法分類精度介于兩者之間，平均Kappa系數約為88.6％；同時AUWEM漏檢與虛警總錯誤率也要少于NDWI的分類結果和最大似然法的分類結果，AUWEM平均漏檢與虛警總錯誤率約為11.9％，最大似然法平均漏檢與虛警總錯誤率約為18.2％，NDWI平均漏檢與虛警總錯誤率約為22.1％％。除此之外，AUWEM具有更高的水體邊緣檢測精度。

附圖說明

圖1是本發明具體實施方式中的NNDWI的水體提取結果，其中圖(a)是從資源3號衛星獲得的遙感影像，圖(b)是NNDWI1水體提取結果，圖(c)是NNDWI2水體提取結果，圖(d)NNDWI水體提取結果；

圖2是本發明具體實施方式中的膨脹約束過程示意圖；

圖3(a)～(e)是本發明具體實施方式中的陰影區域不同波譜特性曲線；

圖4是本發明具體實施方式中的衛星遙感影像的城市水體提取方法流程圖；

圖5是本發明具體實施方式中的不同方法在5個實驗區的實驗結果；

圖6(a)～(f)是本發明具體實施方式中的不同方法在5個實驗區的6個指標直方圖；

圖7是本發明具體實施方式中的水體邊緣精度評估中待評估區域獲取示意圖；

圖8(a)～(c)是本發明具體實施方式中的不同方法在5個實驗區水體邊緣檢測精度比較。

具體實施方式

下面結合附圖和具體實施方式對本發明進行詳細說明。

一種衛星遙感影像的城市水體提取方法，以資源3號衛星遙感影像數據為例，包括以下步驟：

步驟1：遙感影像數據的預處理，即對遙感影像數據進行正射校正和大氣校正；對實驗區域的影像利用RPC+30DEM進行無控制點正射校正，采用Feyisa G L et al.的FLAASH(Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes)大氣校正模型進行影像的大氣校正，均在ENVI5.2軟件中完成。其中，資源3號FLAASH大氣校正需要的定標系數可通過網站http://www.cresda.com/CN/Downloads/dbcs/index.shtml進行下載，光譜響應函數可通過網站http://www.cresda.com/CN/Downloads/gpxyhs/index.shtml進行下載。

步驟2：預處理后的遙感影像數據，包括藍色band1數據、綠色band2數據、紅色band3數據和近紅外band4數據，選取預處理后的遙感影像數據中的藍色band1數據代替歸一化差異水體指數NDWI(Normalized Difference Water Index)的計算公式中的綠色band2數據，獲得新的歸一化差異水體指數NNDWI1(New Normalized Difference Water Index 1)，新的歸一化差異水體指數NNDWI1的計算公式為：

此表達式即為NNDWI1指數模型。利用此模型在編程軟件Matlab 2014a中通過閾值分割得到NNDWI1閾值分割結果，也即為NNDWI1水體提取結果。

步驟3：對預處理后的遙感影像數據中包括的四個波段數據，即藍色band1數據、綠色band2數據、紅色band3數據和近紅外band4數據進行PCA變換，并將PCA變換后的第一主成分分量Component1替代歸一化差異水體指數NDWI的計算公式中的綠色band2數據，獲得另一個新的歸一化差異水體指數NNDWI2(New Normalized Difference Water Index 2)，即：

其中，Component1表示PCA變換的第一主成分分量。

此表達式即為NNDWI2指數模型。利用此模型在編程軟件Matlab 2014a中通過閾值分割得到NNDWI2閾值分割結果，也即為NNDWI2水體提取結果。

步驟4：將步驟2得到的NNDWI1的閾值分割結果和步驟3中得到的NNDWI2的閾值分割結果進行疊加，將得到的結果定義為新的歸一化差異水體指數NNDWI(New Normalized Difference Water Index)的閾值分割結果，即NNDWI1的閾值分割結果與NNDWI2的閾值分割結果進行疊加，其計算公式為：

NNDWI＝(segmentation_NNDWI1)∪(segmentation_NNDWI2)

式中segmentation_NNDWI1表示NNDWI1指數影像的閾值分割結果，segmentation_NNDWI2表示NNDWI2指數影像的閾值分割結果。此計算公式即為NNDWI指數模型，利用此模型在編程軟件Matlab 2014a中得到NNDWI水體提取結果。

通過實驗分析，實驗中NNDWI1的閾值分割結果所采用的閾值與NNDWI2的閾值分割結果所采用的閾值都設為0，可獲得較為理想的實驗結果。NNDWI的水體提取結果綜合了NNDWI1的水體提取結果和NNDWI2的水體提取結果，避免了單一指數對水體的漏檢現象發生。如圖1(a)～(d)所示。NNDWI1對含泥沙高的水體檢測效果較為明顯，而NNDWI2對受植被光譜信息干擾嚴重的水體較為敏感，因此，實際水體提取中結合NNDWI1的水體提取結果和NNDWI2的水體提取結果，整合生成新的水體提取結果，提高后續水體提取精度。

步驟5：對預處理后的遙感影像數據中的近紅外band4數據進行閾值分割，得到近紅外band4數據的閾值分割結果。

通過分析大量NNDWI水體提取結果的影像資料發現，除部分城區面積較小的人工池塘與湖泊之外，在一般情況下，城市高層建筑物陰影對象的面積一般要小于水體對象的面積。所以在實際的分析中，只對面積相對較小的地物進行檢測，因為在這部分面積相對較小的對象中包括幾乎所有的陰影對象，同時也包括小面積水體對象。

步驟6：對NNDWI的水體提取結果中的大面積水體對象和小面積對象進行分割，NNDWI的水體提取結果中，像素個數大于設定閾值的為大面積水體對象，像素個數小于等于設定閾值的為小面積對象；

小面積對象獲取模型，可表示為：

其中，t為設定的閾值，其值取最大陰影對象的像素個數，component表示NNDWI水體指數提取結果中的離散對象，area(component)表示NNDWI水體指數提取結果中的離散對象的面積，area(component)>t為大面積水體對象，area(component)≤t為小面積水體或者陰影對象。

步驟7：對步驟6中得到的小面積對象進行數學形態學膨脹處理，得到膨脹后的小面積對象，將步驟5得到的近紅外band4數據的閾值分割結果作為約束條件，即采用膨脹后的小面積對象和近紅外band4數據的閾值分割結果求交集的方式對膨脹后的小面積對象進行約束，約束的數學表達式為：

component2＝(dilate_component)∩(segmentation_band4)

式中，dilate_component表示膨脹后的小面積對象，segmentation_band4表示近紅外band4數據的閾值分割結果，component2表示約束后的小面積對象；

該膨脹后的小面積對象更加完整的包括了小面積水體對象與陰影的陰影像素，對膨脹后的小面積對象進行約束也保持了小面積對象的膨脹操作限制在地表真實陰影區域的范圍內。

以步驟6得到的小面積對象為建筑物陰影對象為例的膨脹約束過程示意圖如圖2所示，從資源3號假彩色影像中，分別獲取建筑物陰影對象和Band4影像數據閾值分割結果，將得到的建筑物陰影對象進行數學形態學膨脹處理，獲得建筑物陰影對象的膨脹結果，將建筑物陰影對象的膨脹結果與Band4影像數據閾值分割結果進行求交集處理，即得到在Band4影像數據閾值分割結果約束下的膨脹后的建筑物陰影對象。

步驟8：對步驟7得到的約束后的小面積對象進行陰影檢測與去除，得到小面積水體對象；

在NNDWI水體提取結果中，基本上只包含小面積水體對象與陰影，所以只需要對小面積水體對象與陰影的特征進行研究與分析，找到適合區分小面積水體對象與陰影的特征。在實驗中發現，盡管紋理特征可以很好的描述小面積水體對象與陰影，但由于地物的紋理特征，如：灰度共生矩陣，其計算相對復雜、耗時較長，不太適于用于小面積水體對象與陰影的區分，所以在實驗中采用地物的光譜特征描述小面積水體對象和陰影的像素，以此作為依據區分小面積水體對象與陰影。

通過分析大量的水體與陰影的波譜特性曲線，得到水體像元的波譜關系滿足不等式：

band2>band4

而陰影的像元波譜曲線較為復雜，實驗中分析總結出了以下5種波譜特性曲線，如圖3(a)～(e)所示。

根據上述波譜曲線各波段對應的大小關系，實驗中總結出陰影像元的波譜關系滿足以下三個不等式條件：

由此可知，陰影像元的條件，是指陰影像元的波譜關系滿足以上三個不等式條件。

陰影檢測與去除，是指對每個小面積對象中的每個像元進行波譜關系的描述，并判斷該像元是否滿足陰影像元的條件，記錄并統計每個小面積對象中陰影像元的個數，當一個小面積對象中陰影像元所占比例大于閾值T時，把該小面積對象判定為建筑物陰影對象，小于等于閾值T時小面積對象則判定為小面積水體對象，陰影像元所占比例即為小面積對象中陰影像元的個數與該小面積對象中總像元個數的比值，區分小面積對象中小面積水體對象和陰影對象的函數表達式為：

式中，n表示某一小面積對象中總像元個數，m為該小面積對象中陰影像元的個數。閾值T是通過實驗統計得來的，對資源3號遙感影像數據的陰影像元進行統計發現，當T取0.5時可以很好的區分水體與陰影對象。

步驟9：將步驟6中得到的大面積水體對象與步驟8中得到的小面積水體對象進行疊加，即將步驟6中得到的大面積水體對象和步驟8中得到的小面積水體對象求并集，得到衛星遙感影像的城市水體提取結果。

一種衛星遙感影像的城市水體提取方法的總體流程圖如圖4所示。

為驗證方法的有效性，分別采用NDWI方法水體提取結果和最大似然法(MaxLike)水體提取結果進行對比實驗。在中國境內選取了5處不同地區且具有不同周邊環境的影像用于實驗，它們包括湖泊和河流，分別位于中國地區的北京市、武漢市、蘇州市、廣州市，其中武漢市選用了兩個不同覆蓋區域的影像。資源3號影像詳細信息描述如表2所示，實驗中詳細的參數的設置情況如表3所示，其中band4的數值首先按下式歸化到0-255取值范圍，然后再選取閾值進行分割的。

表2實驗數據詳細信息

表3不同實驗地區閾值設定(其中T1，T2，T3，分別為NNDWI1、NNDWI2、band4的分割閾值)

為便于不同方法分類結果的目視判讀與分析，對正確分類的水體像元賦予灰色，正確分類的非水體像元賦予黑色，錯誤分類的像元賦予白色，實驗結果如圖5所示。從圖5的分類結果目視判讀可以發現，本發明提出的AUWEM的水體提取分類精度要好于NDWI的水體提取分類精度和最大似然法的水體提取分類精度。AUWEM對水體邊緣混合像元可以很好的分類(結合Beijing、Wuhan_1和Wuhan_2地區的水體分類結果)、對細小的池塘水體檢測性能好于NDWI和最大似然法(結合Suzhou地區水體分類結果)、對房屋建筑物的陰影可以很好的去除(結合Suzhou和Wuhan_2地區水體分類結果)。

不同實驗區域三種方法水體提取分類精度比較統計結果如表4所示，從表4的統計結果中發現，AUWEM水體提取分類精度要高于NDWI和最大似然法。AUWEM在這5個實驗區域的分類精度是最高，平均Kappa系數達93.0％；而NDWI的分類精度最低，平均Kappa系數約為86.2％；最大似然法分類精度介于兩者之間，平均Kappa系數約為88.6％。

表4不同實驗區域三種方法的精度統計

為了更加詳細的評估AUWEM對水體的提取精度，采用生產者精度、用戶精度、Kappa系數、漏檢率、虛警率、總錯誤率6個指標來描述方法的水體提取精度。三種不同方法在5個實驗地區的6個指標直方圖如圖6(a)～(f)所示，從各個指標的直方圖可以發現，AUWEM的水體提取分類精度要高于NDWI的水體提取分類精度和最大似然法的水體提取分類精度。AUWEM在水體提取的虛警率方面，除在Suzhou地區達到9.1％左右，其他的實驗區域都是低于5％；在水體漏檢率方面，5個地區漏檢率都要明顯少于NDWI和最大似然法。當水體提取的虛警率與水體漏檢率都低時，總錯誤率自然也就是最低。從圖6(a)～(f)中可以看出AUWEM總錯誤率最低，其次是最大似然法，總錯誤率最高的是NDWI，其中AUWEM平均總錯誤率約為11.9％，最大似然法平均總錯誤率約為18.2％，NDWI平均總錯誤率約為22.1％％。

對于水體提取分類生產者精度，AUWEM水體提取分類生產者精度最高，平均約為91.6％；最大似然法水體提取分類生產者精度次之，平均約為84.8％；NDWI的水體提取分類生產者精度最低，平均約為81.6％。在水體提取用戶精度方面，最大似然法的水體提取用戶精度要大于的AUWEM和NDWI，平均精度用戶高達96.6％；AUWEM方法次之，平均精度用戶約為96.4％；最差的仍然屬于NDWI方法，平均精度用戶約為95.7％。

為了更加客觀評價三種方法提取的水體的邊緣檢測精度，設計以下方法來進行精度評價，方法具體實施描述如下：

首先利用參考影像采用Canny算子獲取三種方法提取的水體的邊界線；

對獲取的邊界線進行數學形態學的膨脹處理，建立以邊界線為中心半徑為4個像元的緩沖區域；

然后對緩沖區域的像元進行判斷，假設緩沖區域總像元值為N，正確分類的像元數目的N_R，漏檢像元數目為No，虛警數目為Nc，則：

其中，A+Eo+Ec＝100％。A表示邊緣像元正確劃分類別的比例，在這里稱它為邊緣像元正確分類精度；Eo表示邊緣像元漏檢的比例，稱它為邊緣像元漏檢率；Ec表示邊緣像元虛警的比例，稱它為邊緣像元虛警率。將經過數學形態學膨脹處理的水體的邊界線與從資源3號衛星獲取的遙感影像進行疊加，得到邊緣精度評估中待評估區域。邊緣精度評估中待評估區域獲取示意圖如圖7所示。

根據上述方法統計了實驗區域的水體邊緣檢測精度，分別統計了邊緣像元虛警率(Commission Error)、邊緣像元漏檢率(Omission Error)和邊緣像元正確分類精度(Accuracy of edge detection)，并對實驗結果進行比較，如圖8(a)～(c)所示。從圖8的比較結果可以很清晰的發現：AUWEM方法邊緣像元正確分類精度要好于NDWI和最大似然法。其中AUWEM方法水體邊緣像元正確分類精度最高達93.7691％(Guangzhou地區)，最小精度為79.5798％(Wuhan_2地區)；NDWI水體邊緣像元正確分類精度最高達84.0917％(Suzhou地區)，最小精度為69.8310％(Beijing地區)；最大似然法水體邊緣像元正確分類精度最高達85.8149％(Guangzhou地區)，最小精度為69.7974％(Wuhan_2地區)。