多光譜光學遙感影像水體提取模型

鄧開元，任 超

桂林理工大學測繪地理信息學院, 廣西 桂林 541004

地表水變化對生態(tài)、社會和經(jīng)濟的影響是當前研究的熱點[1]。光學遙感水資源研究主要可以分為：氣候和水資源[2]、不同類型水體的提取[3-5]、洪水監(jiān)測[6]、河流形態(tài)[7]、水體長時間序列[8]、海岸帶變化[9]。目前，Landsat影像是用于變化分析，識別和檢測地表水的最廣泛使用的數(shù)據(jù)源[10]。在基于Landsat遙感影像水體提取的研究中，多數(shù)是針對大區(qū)域的河流、湖泊等水體的提取[11-14]。近年來，一些學者采用Sentinel 2影像對大區(qū)域的地表水進行提取研究，取得了良好的效果[15-16]。目視解譯是最初的光學遙感影像提取水體的方法[17]，其缺點是效率低、成本高，沒有利用不同地物具有不同反射率的特性。遙感影像的水體信息提取是基于電磁波譜上的反映特性[18]，文獻[19]通過使用Landsat 4影像的6個波段提出纓帽變換濕度指數(shù)(tasseled cap wetness,TCW)。受到歸一化差異植被指數(shù)[20]的啟發(fā)，文獻[21]基于Landsat影像的綠波段和近紅外波段提出了歸一化差異水體指數(shù)(normalized difference water index,NDWI)，NDWI是簡單的水體指數(shù),該指數(shù)能較好地抑制植被信息，但受薄云和山體陰影影響較大。文獻[22]使用短波近紅外波段代替近紅外波段，改良了NDWI，提出了改進的歸一化差異水體指數(shù)(modified normalized difference water index,MNDWI)。MNDWI指數(shù)能突出水體信息，但也容易把水體和陰影混淆?；诜诸惗鏄涞乃枷?，文獻[23]提出了針對有陰影的場景自動水體提取指數(shù)AWEIsh和針對沒有陰影的場景自動水體提取指數(shù)AWEInsh，二者都是多波段水體指數(shù)。AWEI指數(shù)在提出時，聲明該指數(shù)不適用于高反射率的地物，比如冰雪、白色建筑物等。文獻[24]在全球大范圍多個地區(qū)進行測試了自身的穩(wěn)定性和可行性,提出了水體指數(shù)(water index 2015,WI2015)。針對復雜環(huán)境，文獻[25]成功地構建了一個多波段水體指(multi-band water index,MBWI)，MBWI一共使用了5個波段并且在多個復雜環(huán)境的地區(qū)測試了MBWI的有效性。MBWI能較好地抑制非水體信息的能力，但其缺點之一是不適用提取細小水體。面向?qū)ο筇崛∷w信息的方法[26],在抑制陰影信息的同時，能夠準確提取細小水體。但是面向?qū)ο蟮姆椒ǖ娜秉c之一是：不具有普遍性，難以在大范圍下實現(xiàn)高精度的水體信息提取。

在地形復雜的環(huán)境下，如高山積雪較多的區(qū)域，對水體提取產(chǎn)生干擾的因素有：陰影、冰雪等[27-28]。如何減少遙感影像的噪聲干擾仍是水體信息提取的難點。本文分析多種地物的反射率，綜合考慮多波段間的關系，通過建立六元不等式組，設計一個顧及多個地物噪聲影響的光學遙感影像歸一化多波段水體指數(shù)(normalized difference multi-band water index,NDMBWI)。

1 多光譜影像水體指數(shù)模型構建

1.1 已有的水體指數(shù)

選取7種典型的、常用的水體指數(shù)進行研究，7種常用的水體指數(shù)見表1。這7種水體指數(shù)以綠波段和近紅外波段、短波近紅外波段為主要影響成分。TCW、MNDWI、AWEI、WI2015主要考慮綠波段和短波近紅外波段的反射率，因為水體的綠波段反射率大于短波近紅外波段的反射率，而其他非水體地物的綠波段小于短波近紅外波段的反射率。正常情況下，這些指數(shù)都可以區(qū)分水體與非水體。但是，TCW、MNDWI、AWEIsh、WI2015會把可見光波段反射率大于短波近紅外波段的非水體地物分類為水。文獻[25]指出，雪就是可見光波段的反射率遠遠大于短波近紅外波段的地物。

1.2 模型的機理

本文的技術路線如圖1所示。在不同的時段、不同的地區(qū)，從多幅Landsat 8、Sentinel 2、高分一號(GF-1)影像里選擇不同地物不同波段的真實地表反射率光譜值進行擬合，具體采樣點的信息見表2。內(nèi)陸水體、海洋水體、陰影、雪、植被、旱地、建筑物、云的平均反射率如圖2所示。圖2曲線是基于表2中的采樣點求出來的平均值。水體的反射主要在藍綠光波段，其他波段的吸收率很強，尤其是在近紅外波段、短波近紅外波段有強烈的吸收帶。研究已有的水體指數(shù)可以從各個光譜間的關系作為切入點，本文以不同波段與綠波段的比值，闡述已有水體指數(shù)模型的機理及其存在的缺陷，并建立本文水體指數(shù)NDMBWI。表3為各個光譜與綠波段的關系，表3中的數(shù)值可看成圖2曲線量化后的結果。圖1中，可以把TCW、NDWI、MNDWI、AWEIsh、AWEInsh、WI2015近似看成雙波段水體指數(shù)。水體指數(shù)的模型機理在于不同波段的對不同地物具有不一樣的敏感性和有效性。雪、云的反射特性與水體相似，表3中，部分地物(內(nèi)陸水體、海洋水體、雪、云)的Q4、Q5、Q6的值都小于1，從機理上說明：主要依賴綠波段、近紅外波段、短波近紅外波段的雙波段水體指數(shù)在有雪環(huán)境下不能有效地提取水體。

表1 已有的水體提取指數(shù)

圖1 本文技術路線Fig.1 Technical route of this paper

表3中的數(shù)值是來源于Landsat 8、Sentinel 2、GF-1在各個地物、各個波段上的平均值，然后各自相除得到的比重。本文試驗數(shù)據(jù)源為Landsat 8、Sentinel 2、GF-1，具體信息見表4。

表4中，Landsat 8、Sentinel 2、GF-1三者的中心波長不完全一致，這也導致了已有的水體指數(shù)模型在用于不同的數(shù)據(jù)源時，表現(xiàn)出閾值跳動的問題。經(jīng)過研究，已有的水體指數(shù)的模型可用式(1)表示

(1)

式中，Y表示水體指數(shù)的值；x1、x2、x3、x4、x5、x6是6個波段前的系數(shù)；ρ1至ρ6是藍波段、綠波段、紅波段、近紅外波段、短波近紅外1波段、短波近紅外2波段的反射率。

圖2 擬合多個地物的反射率曲線Fig.2 Reflectance curves of multiple features

式(1)整體除以綠波段，可變形為

(2)

表2 數(shù)據(jù)采樣點涉及的數(shù)據(jù)

表3 各個光譜關系

式中，Y表示水體指數(shù)的值；x1、x2、x3、x4、x5、x6是6個波段前的系數(shù)；ρ1至ρ6是藍波段、綠波段、紅波段、近紅外波段、短波近紅外1波段、短波近紅外2波段的反射率。

根據(jù)水體的水體指數(shù)的值理論上要大于α，非水體的水體指數(shù)的值理論上要小于α。式(2)中，ρ1、ρ2、ρ3、ρ4、ρ5、ρ6的值用比值系數(shù)代替,得到以下不等式組

(3)

式中，a、b、c、d、e、f、g、h分別代表求得的比值系數(shù)。

表4 數(shù)據(jù)源信息

在本文把α設置為0，目的是希望本文指數(shù)能夠以穩(wěn)定的閾值提取水體信息，到達自動化提取水體的目的。結合式(3)與表3的參數(shù)，組成一個六元不等式組

(4)

式中，x1、x2、x3、x4、x5、x6分別是藍波段、綠波段、紅波段、近紅外波段、短波近紅外1波段、短波近紅外2波段的系數(shù)。

式(4)的系數(shù)源于Landsat 8、Sentinel 2、GF-1在各個地物、各個波段上的平均值，若是直接賦予式(4)每項等于定值，會因為推導不嚴謹，導致水體指數(shù)模型出現(xiàn)錯漏現(xiàn)象，本文通過最小二乘法的方式進行消元、擬合，求出不同定值下的解。

由表3中可以發(fā)現(xiàn)，已有水體指數(shù)Q1、Q2、Q3、Q4、Q5、Q6與1的大小關系并不能完全地實現(xiàn)水體的識別。例如，NDWI由綠波段與近紅外波段組成，即是從Q4的值可以解析NDWI對水體信息識別的能力，大于1的數(shù)值的地物包括陰影、植被、旱地、建筑物、云，小于1的值包括內(nèi)陸水體、海洋水體、雪，此時導致NDWI無法以固定的閾值區(qū)分雪、水體。MNDWI、AWEIsh、AWEInsh、WI2015近似看成由綠波段與短波近紅外波段組成的模型，可從Q5的值可以解析此類型模型對水體信息識別的能力，大于1的數(shù)值的地物包括陰影、植被、旱地、建筑物，小于1的值包括內(nèi)陸水體、海洋水體、雪、云。此時，以上指數(shù)只考慮Q4或Q5(雙波段)，并沒有綜合考慮Q1、Q2、Q3、Q4、Q5、Q6的內(nèi)在聯(lián)系，所以本文建立六元不等式組，綜合考慮多波段間的關系，設計一個顧及多個地物噪聲影響的水體指數(shù)模型。

式(4)沒有精確解，其原因是：式(4)是超定方程組的類型，即方程個數(shù)大于未知量個數(shù)。式(4)的解集是一個區(qū)域，解多元不等式(組)就是求它的各個不等式的解集的交集。為了使本文水體指數(shù)適用于國產(chǎn)遙感影像(沒有短波近紅外1波段、短波近紅外2波段)，需要提前把x5、x6設置為0。本文通過逐步消元、擬合，求出其中的一種解：x1、x2、x3、x4、x5、x6分別是-1、3、2、-5、0、0，即初步構建的水體指數(shù)為

NDMBWI=3ρ2-ρ1+2ρ3-5ρ4

(5)

式中，ρ1、ρ2、ρ3、ρ4代表藍波段、綠波段、紅波段、近紅外波段的反射率。

式(5)的系數(shù)是可行域中的簡化解，目的是為了更好推廣NDMBWI以及方便記憶。另外，可以通過使用式(5)的參數(shù)，反解式(4)進行驗證。因為雙波段水體指數(shù)難以抑制雪、云的影響，本文指數(shù)沒有將綠波段作為唯一的正參數(shù)，額外引入了紅波段作為正參數(shù)。為了使本文指數(shù)歸一化，最終擬定的水體指數(shù)模型為

(6)

式中，ρ1、ρ2、ρ3、ρ4代表藍波段、綠波段、紅波段、近紅外波段的反射率。

本文推算出來的系數(shù)，實則為式(6)的分子。在提取水體時，起決定性作用的是分子而不是分母。這樣做的好處是：預處理時，在大氣校正環(huán)節(jié)，某些算法為了減少存儲空間會把結果擴大10 000倍。此時，本文使用了歸一化的形式，則不會受以上的影響，無論大氣校正有無擴大10 000倍，NDMBWI的范圍始終在-1至1的范圍。使用本文水體指數(shù)和表1中列舉的已有的水體指數(shù)進行仿真模擬，在不同的地區(qū)、不同的數(shù)據(jù)源選擇水體、陰影、雪、植被、旱地、建筑物、云等7種地物作為對象，每種地物選擇100個樣本點，不同水體指數(shù)的穩(wěn)定性測試結果如圖3所示。圖3表明，0可以作為NDMBWI的閾值，能夠有效地區(qū)分水體與多種非水體；以上的水體指數(shù)均可區(qū)分水體與植被、旱地、建筑物；除了NDMBWI、MBWI，其他的水體指數(shù)均不能區(qū)分水體與冰雪；MBWI會對部分水體漏分；NDWI的閾值提高到0.2左右可以分區(qū)水體與冰雪。

2 水體提取實例分析

2.1 不同衛(wèi)星影像在有雪環(huán)境下的水體提取

試驗1選擇的地區(qū)位于西藏林芝地區(qū)。數(shù)據(jù)源是Landsat 8影像和Sentinel 2影像，影像拍攝時間均為2020年2月13日。圖3中得知：8種指數(shù)以閾值為0均能將水體與部分非水體(植被、旱地、建筑物)區(qū)分。所以在本次試驗中，NDMBWI、MBWI、WI2015、AWEInsh、AWEIsh、MNDWI、TCW統(tǒng)一以0作為閾值進行水體提取，而NDWI則是使用0.2作為閾值。提取雅魯藏布江水體的結果如圖4所示。圖4(a)為Landsat 8影像的原始圖。圖4(b)為Sentinel 2影像的原始圖。對比圖4(c)、(d)可以發(fā)現(xiàn)，NDMBWI應用于Landsat 8影像、Sentinel 2影像時，均可以抑制了雪的影響，完成了提取水體信息。

在試驗1中，選擇100個水體、非水體樣本組成混淆矩陣，計算出的Kappa系數(shù)、總體精度、錯分誤差、漏分誤差、用戶精度、生產(chǎn)者精度見表5，使用Sentinel 2影像的效果均優(yōu)于Landsat 8。為了進一步研究水指數(shù)的閾值，研究NDMBWI、MBWI、WI2015應用于Sentinel 2影像的數(shù)值分布，三者的直方圖如圖5所示。在直方圖中，最優(yōu)閾值一般選定在“雙峰”之間的低谷區(qū)間。圖5(a)驗證了在試驗1中，0可以作為NDMBWI的最優(yōu)閾值。圖5(c)并無明顯的低谷，說明WI2015本身難以抑制冰雪的干擾。

表5 分類精度混淆矩陣分析

2.2 有云環(huán)境下的高分一號影像水體提取

云一直是光學遙感影像難以處理的問題之一，本文提出的NDMBWI適用于存在少量云的情景。云一般都有陰影，傳統(tǒng)的水體指數(shù)無法處理陰影。由仿真模擬的結果可知：NDMBWI能夠以0為閾值區(qū)分水體與云、陰影。GF-1是中國自主研制的第一顆民用高分辨率光學衛(wèi)星，為我國水環(huán)境的遙感監(jiān)測提供更高空間分辨率的光學遙感數(shù)據(jù)源。但是，GF-1只有可見光波段與近紅外波段，使用了短波近紅外波段的MBWI、WI2015、AWEI、MNDWI、TCW等指數(shù)不能應用在GF-1。試驗2以中國香港的迪士尼公園為試驗區(qū)域，影像的拍攝日期為2018年10月2日。使用NDMBWI提取水體的對比結果如圖6所示。

圖3 不同水體指數(shù)的穩(wěn)定性測試結果Fig.3 Stability test results of different water indexes

圖4 雅魯藏布江水體提取結果Fig.4 Water extraction results of Yarlung Zangbo river

圖5 NDMBWI、MBWI、WI2015的直方圖Fig.5 Histogram of NDMBWI, MBWI and WI2015

圖6 香港迪士尼地區(qū)水體提取結果Fig.6 Water extraction results in Hong Kong Disneyland

圖7作為NDMBWI在試驗2的直方圖，驗證了0可以作為NDMBWI的最優(yōu)閾值。云是光學影像無法處理的痛點，去云步驟一般分為兩步：①識別云像素；②把云像素設置為背景值(無效值)。本文在設計指數(shù)模型時把云考慮到非水體地物的因子，所以通常情況下NDMBWI在會把云識別為非水體地物，云通過NDMBWI得到的值為負值。即NDMBWI無須做去云步驟。但是，值得說明的是NDMBWI適用于云量為20%以下的光學遙感影像。也就是說NDMBWI是無法適用于云量過多的場景。

圖7 NDMBWI的直方圖Fig.7 histogram of NDMBWI

2.3 多地區(qū)水體提取

選取表2中的齊齊哈爾、巴彥淖爾、梧州、三亞作為試驗3，從多個不同的地區(qū)，復雜的地形環(huán)境下提取不同類型的水體，驗證本文水體指數(shù)模型的穩(wěn)定性，圖8中的NDMBWI均以0作為閾值。

3 結論和討論

如何減少噪聲干擾對遙感影像的不利影響是水體提取的難點，本文設計一種歸一化多波段水體指數(shù)NDMBWI，避免積雪、陰影等的影響。試驗1使用相同時相、同一地區(qū)的Landsat 8影像、Sentinel 2影像，以西藏林芝地區(qū)為試驗區(qū)域，試驗結果說明NDMBWI能抑制雪的影響。試驗2使用GF-1影像，在影像含有少量云的情況下完成水體提取，驗證NDMBWI可以抑制云、陰影的影響。閾值，對于水體指數(shù)是至關重要的。從水體指數(shù)的構建開始，要注意閾值的選取問題。試驗3選擇多個不同地區(qū)進行試驗，以說明本文水體指數(shù)模型的穩(wěn)健性。今后，構建新的光學水體指數(shù)時，應該考慮多個波段間的敏感性和有效性。本文的創(chuàng)新點在于：充分利用可見光波段、近紅外波段的信息，不需額外借助輔助數(shù)據(jù)，即可排除雪、云、陰影的影響，完成水體信息的提取。

圖8 多地區(qū)水體提取結果Fig.8 Water extraction results in multiple regions

從理論上，任何的水體指數(shù)都無法保證提取所有的水體，本文提出水體指數(shù)雖然在某些情況下能得到較高的精度，并在冰雪、陰影、建筑、云等區(qū)分具有一定的優(yōu)勢，然而因為自然水體的復雜性，并不是任意一種水體指數(shù)能區(qū)分所有水體。也就是說，NDMBWI的最優(yōu)閾值不能一直保持在0，例如在錢塘江懸浮泥沙較多的出?？诘貐^(qū)提取水體時，NDMBWI的閾值應該適當降低。目前主流的水體指數(shù)的研究大多數(shù)是固定閾值。本文的下一步工作是：基于原有的水體指數(shù)模型，構建動態(tài)閾值的模型。例如，以NDVI作為NDMBWI的閾值，使用NDMBWI與NDVI模型進行雙通道并行處理，在NDMBWI的矩陣中，大于NDVI的像素點可以歸為水體，反之為非水體。本文方法穩(wěn)健，可以作為基于深度學習的水體檢測分類的基礎，解決難以區(qū)分雪、水的問題。