范亞洲，張 珂,3,4，劉林鑫，晁麗君，姚 成

(1.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098;2.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098; 3.長(zhǎng)江保護(hù)與綠色發(fā)展研究院，江蘇 南京 210098;4.中國(guó)氣象局-河海大學(xué)水文氣象研究聯(lián)合實(shí)驗(yàn)室，江蘇 南京 210098)

水利工程設(shè)施是水資源管理與規(guī)劃的重要設(shè)施[1]。隨著遙感技術(shù)在水利工程領(lǐng)域的廣泛應(yīng)用，通過(guò)遙感影像提取水體信息已經(jīng)成為水資源調(diào)查及宏觀監(jiān)測(cè)[2-6]、復(fù)核庫(kù)容曲線等領(lǐng)域中的關(guān)鍵技術(shù)[7-8]。遙感影像提取水體信息的研究方法包括單波段閾值法[9]、多波段譜間關(guān)系法[10]、水體指數(shù)法以及決策樹法[11]等。王國(guó)興等[12]提取了新安江水庫(kù)不同時(shí)相的水體面積，結(jié)合實(shí)測(cè)水位擬合了水位-面積關(guān)系曲線；曹明亮等[13]通過(guò)遙感影像提取水庫(kù)庫(kù)區(qū)的水體面積，提出了小水庫(kù)攔洪計(jì)算方法；曹榮龍等[14]通過(guò)分析地物光譜特征構(gòu)建了修訂型水體指數(shù)，提取了密云水庫(kù)水面面積的動(dòng)態(tài)變化；張莉芳等[15]基于遙感技術(shù)提取了柘林水庫(kù)庫(kù)區(qū)內(nèi)的水體面積，進(jìn)而對(duì)庫(kù)容曲線進(jìn)行復(fù)核；戚曉明等[16]基于21景ETM+遙感影像提取了洪澤湖的水面變化序列數(shù)據(jù)，并推導(dǎo)了洪澤湖庫(kù)容曲線。上述研究均使用人工判讀或經(jīng)驗(yàn)確定的閾值來(lái)提取庫(kù)區(qū)水體信息，并沒有對(duì)分割閾值做出客觀的判斷。

高分一號(hào)(GF-1)衛(wèi)星于2013年4月26日發(fā)射成功，搭載了兩臺(tái)2 m分辨率全色和8 m分辨率多光譜相機(jī)(PMS)外，還搭載了4臺(tái)16 m分辨率多光譜相機(jī)(WFV)，具有壽命長(zhǎng)、效能高等特點(diǎn)[17]，提高了我國(guó)高分辨率影像數(shù)據(jù)的自給率[18]，為遙感影像信息提取及特征識(shí)別帶來(lái)了新的機(jī)遇[19]。本文開展了以GF-1衛(wèi)星遙感影像為數(shù)據(jù)源的東方紅水庫(kù)庫(kù)區(qū)水體信息優(yōu)化提取研究，提出基于最大類間方差的庫(kù)區(qū)水體迭代優(yōu)化提取方法。

1 方法的提出

1.1 歸一化差分水體指數(shù)法

歸一化差分水體指數(shù)(normalized difference water index， NDWI)法由McFeeter[20]提出，根據(jù)水體在綠光波段上具有很強(qiáng)反射性、在近紅外波段具有很強(qiáng)吸收性的特點(diǎn)，建立波段差異比值，從而增強(qiáng)水體信息，計(jì)算公式為

INDWI=(λ1-λ2)/(λ1+λ2)

(1)

式中：INDWI為歸一化差分水體指數(shù)；λ1為GF-1衛(wèi)星的綠光波段反射率；λ2為GF-1衛(wèi)星近紅外波段反射率。

通過(guò)地物光譜分析及水體提取實(shí)驗(yàn)，在突出庫(kù)區(qū)主要水體信息的同時(shí)，為了較好地抑制植被等其他地物信息，本文將初始閾值設(shè)為T0=-0.2，即INDWI>T0得到的結(jié)果為初步提取的水體信息。

1.2 最大類間方差迭代法及形態(tài)學(xué)膨脹算法

最大類間方差法是由Otsu[21]1979年提出的一種無(wú)參數(shù)無(wú)監(jiān)督的閾值分割算法[22]。該算法根據(jù)圖像的灰度特性，分割閾值從最小灰度值遍歷到最大灰度值，在某一灰度值下，將圖像分成背景和前景兩部分，當(dāng)兩者類間方差最大時(shí)，背景和前景這兩部分被錯(cuò)分的概率是最小的[23]，對(duì)應(yīng)最大類間方差的閾值即為所求的自適應(yīng)閾值。

本文在最大類間方差法基礎(chǔ)上，提出了最大類間方差迭代法來(lái)提高水體提取的精度。首先通過(guò)邊緣檢測(cè)的方法得到初步提取的庫(kù)區(qū)水體信息范圍，以此基于形態(tài)學(xué)膨脹算法建立庫(kù)區(qū)水體的緩沖區(qū)[24]，再對(duì)緩沖區(qū)包圍的區(qū)域應(yīng)用最大類間方差法，自適應(yīng)地確定庫(kù)區(qū)水體與周邊地物的分割閾值，最后通過(guò)多次迭代后計(jì)算得到優(yōu)化提取的庫(kù)區(qū)水體面積，水體優(yōu)化提取流程見圖1。

圖1 水體優(yōu)化提取流程

本文結(jié)構(gòu)元素B是數(shù)值均為1的行列數(shù)相等且為奇數(shù)的矩陣，當(dāng)其大小發(fā)生變化，經(jīng)膨脹算法計(jì)算得到的緩沖區(qū)范圍也會(huì)不一致，這導(dǎo)致了水體信息不斷變化，尤其當(dāng)結(jié)構(gòu)元素太大時(shí)，會(huì)造成水體信息的誤提。為了提取準(zhǔn)確、穩(wěn)定的庫(kù)區(qū)水體信息，通過(guò)試驗(yàn)，本文對(duì)優(yōu)化提取算法進(jìn)行了6次迭代計(jì)算，迭代過(guò)程中結(jié)構(gòu)元素的大小分別為3×3、5×5、…、13×13。通過(guò)本文提出的最大類間方差迭代法計(jì)算得到的相鄰水體面積，當(dāng)兩者差值的絕對(duì)值最小時(shí)(公式(2))，那么這一水體面積即為最優(yōu)提取的水體面積(公式(3))，相應(yīng)得到的結(jié)構(gòu)元素最優(yōu)，緩沖區(qū)為最優(yōu)緩沖區(qū)。

1.3 精度評(píng)價(jià)方法

采用決定系數(shù)(coefficient of determination，R2)和均方根誤差(root mean square error， RMSE)來(lái)驗(yàn)證遙感提取水體面積與參考水體面積的一致性，同時(shí)為了驗(yàn)證水體信息的分類精度，各景影像均選取了水體和非水體各100個(gè)樣本像元，由混淆矩陣計(jì)算了各次優(yōu)化結(jié)果及最優(yōu)結(jié)果的總體分類精度(overall accuracy， OA)和Kappa系數(shù)(Kappa coefficient， KC)。

為了統(tǒng)一4項(xiàng)基礎(chǔ)評(píng)價(jià)指標(biāo)，本文分別將各項(xiàng)指標(biāo)歸一化后取平均值建立綜合指標(biāo)，即平均歸一化指數(shù)(average normalized index， ANI)，RMSE是水體遙感提取面積與參考水體面積之間誤差的度量，因此最大值歸一化后為0，最小值歸一化為1。

2 實(shí)例驗(yàn)證

2.1 研究區(qū)概況與數(shù)據(jù)處理

安徽省屯溪流域東方紅水庫(kù)(圖2)地處新安江水系二級(jí)支流東亭河上，水庫(kù)壩址以上的集水面積為60.0 km2，主河道長(zhǎng)11.8 km，平均坡降為21.0‰，總庫(kù)容為2 445萬(wàn)m3，設(shè)計(jì)灌溉面積 30 km2，多年平均發(fā)電量131萬(wàn)kW·h，是以防洪、灌溉為主，結(jié)合發(fā)電、城鎮(zhèn)供水等綜合利用功能的中型水利樞紐工程。

圖2 研究區(qū)示意圖

2020年2月陸面過(guò)程資料分布式數(shù)據(jù)存檔中心(LP DAAC)公布了1弧秒分辨率NASADEM數(shù)據(jù)集產(chǎn)品(https://lpdaac.usgs.gov/news/release-nasadem-data-products/)，本文以此數(shù)據(jù)提取了研究區(qū)流域。通過(guò)現(xiàn)場(chǎng)查勘，收集了東方紅水庫(kù)庫(kù)容曲線及2015年的水庫(kù)水文資料。2015年?yáng)|方紅水庫(kù)集水流域平均降水量為2 411.6 mm，年平均降水頻率為15.1%，其中1—2月降水量較少，5—6月降水量最大。為體現(xiàn)雨期前后庫(kù)區(qū)水體信息的不同特征，選取了2015年10景衛(wèi)星影像(表1)，數(shù)據(jù)源于中國(guó)資源衛(wèi)星應(yīng)用中心(http://www.cresda.com/CN/index.shtml)提供的GF-1 WFV數(shù)據(jù)，影像空間分辨率為16 m。

本研究采用的數(shù)據(jù)為L(zhǎng)evel 1A遙感影像(預(yù)處理級(jí)輻射校正影像產(chǎn)品)，在后續(xù)分析前需對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，包括輻射定標(biāo)、大氣校正[25]、正射校正、配準(zhǔn)及研究區(qū)裁剪。由于數(shù)據(jù)較多，對(duì)遙感影像數(shù)據(jù)分別處理會(huì)受到較大的人為影響，因此對(duì)影像進(jìn)行了批量預(yù)處理。

2.2 水體提取結(jié)果

首先根據(jù)東方紅水庫(kù)的水位-面積關(guān)系(圖3)和每景影像對(duì)應(yīng)日期的實(shí)測(cè)水位(表1)推導(dǎo)相應(yīng)日期水庫(kù)的參考水體面積A(表2)。

圖3 東方紅水庫(kù)水位-面積關(guān)系曲線

表1 影像數(shù)據(jù)集及實(shí)測(cè)庫(kù)水位

從表2看出，初步提取的水體信息與參考水體面積相比，在雨期前4景影像(1月和2月)水庫(kù)缺水期有較大的誤差，平均誤差約為0.80 km2，這首先是因?yàn)殄e(cuò)提了庫(kù)區(qū)淺水區(qū)大部分的裸地，其次本文依據(jù)的水位-面積關(guān)系曲線是在289 m以上，在缺水期通過(guò)擬合得到的參考庫(kù)區(qū)水體面積會(huì)有所偏小，但是在雨期(5月和8月)及雨期后(10月)水庫(kù)蓄水期提取的效果較好，平均誤差約為0.16 km2。在6次優(yōu)化提取的結(jié)果中，僅20150513景影像的提取結(jié)果有所增大，其他影像的水體面積均有減小。根據(jù)本文對(duì)最優(yōu)提取水體的定義，表2列出了各景影像的最優(yōu)水體面積。

圖4為初步提取的庫(kù)區(qū)水體信息及最優(yōu)緩沖區(qū)。從圖4可以看出，NDWI法結(jié)合經(jīng)驗(yàn)閾值初步提取的水體信息在20150207景影像中混淆了庫(kù)區(qū)水體和庫(kù)區(qū)周邊干涸的裸地，在20150212景影像中有大面積的裸地及建筑地物誤提為水體信息。

表2 庫(kù)區(qū)參考水體面積及遙感提取水體面積

(a) GF-1_20150102

圖5為最優(yōu)提取的水體信息，提取結(jié)果有效消除了影像中建筑等非水體信息的干擾，較準(zhǔn)確地提取了東方紅水庫(kù)在雨期前(1月、2月)、雨期(5月、8月)及雨期后(10月)的水體信息。1月和2月降水較少，庫(kù)區(qū)水體也隨著補(bǔ)給下游而逐漸減少；5月和8月降水量增大，庫(kù)區(qū)的來(lái)水量大于出水量，水庫(kù)開始蓄水，水體面積也隨之增加，而且提取結(jié)果較好地保留了8月兩景影像中細(xì)小河道的水體；10月水庫(kù)水體的面積趨于穩(wěn)定。

2.3 精度分析

基于上述水體初步提取、6次迭代提取以及最優(yōu)提取的3類結(jié)果，表3為各類提取結(jié)果4項(xiàng)基礎(chǔ)評(píng)價(jià)指標(biāo)的計(jì)算結(jié)果，其中R2和RMSE分別為十景影像水體提取面積與參考水體面積的決定系數(shù)和均方根誤差，OA和KC分別是十景影像總體分類精度和Kappa系數(shù)的均值。從表3可以看到，本文定義的最優(yōu)水體提取結(jié)果一致性方面相較初步提取結(jié)果顯著提高，R2從0.013 9提高到了0.972 0，RMSE從0.282 2 km2減小到了0.059 2 km2；總體分類精度與Kappa系數(shù)較初步提取分別提高了9.36%和24.09%。此外，第二次和第三次迭代中水體提取的R2分別為0.976 0和0.972 1，而最優(yōu)提取結(jié)果為0.972 0，且RMSE也比最優(yōu)提取的結(jié)果小，可見這兩次計(jì)算的提取結(jié)果一致性要高于本文最優(yōu)提取結(jié)果，但最優(yōu)提取的總體分類精度為0.894 0，Kappa系數(shù)為0.788 0，分類精度比第二次和第三次都更加理想；同理，第六次迭代中的總體分類精度為 0.899 5，Kappa系數(shù)為0.799 0，雖然分類精度高于最優(yōu)提取結(jié)果，但是一致性方面并不滿意，最優(yōu)水體提取結(jié)果的R2和RMSE分別為0.972 0和0.059 2 km2。

(a) GF-1_20150102

表3 各類水體提取結(jié)果精度評(píng)價(jià)

對(duì)比基礎(chǔ)指標(biāo)即可發(fā)現(xiàn)初步提取的結(jié)果較差，因此不參與ANI的計(jì)算，通過(guò)計(jì)算ANI發(fā)現(xiàn)，最優(yōu)提取水體的整體精度最高，ANI為0.7378，相比第二次、第三次和第六次的迭代計(jì)算結(jié)果，精度分別提高了2.67%、13.02%和12.85%，平均提高了9.51%。

3 結(jié) 語(yǔ)

本文在最大類間方差法的基礎(chǔ)上提出了水庫(kù)水體的最大類間方差迭代遙感提取方法，并以東方紅水庫(kù)為研究區(qū)，對(duì)該方法進(jìn)行了驗(yàn)證分析，結(jié)果表明，相較于通過(guò)經(jīng)驗(yàn)閾值初步提取的水體面積信息，該方法較完整地提取了在不同時(shí)期庫(kù)區(qū)水體的動(dòng)態(tài)變化特征，在一致性方面水體提取精度顯著提高；相較于各次迭代中的水體提取結(jié)果，本文提出的庫(kù)區(qū)水體提取方法整體效果良好。但本文僅應(yīng)用廣泛的NDWI方法進(jìn)行了水體提取試驗(yàn)，對(duì)于不同的水體提取方法，提取結(jié)果會(huì)有不同的精度，對(duì)此在未來(lái)的研究中仍需進(jìn)一步深入探討。