結合DEM與淹沒頻率的水庫水體動態遙感提取優化方法

劉詩燕， 蔡曉斌*

(1.中國科學院大學， 北京 100049；2.中國科學院精密測量科學與技術創新研究院環境與災害監測評估湖北省重點實驗室， 武漢 430071)

地表水資源是制約區域經濟發展、保障區域社會穩定、維護區域生態環境的關鍵[1].水庫是地表水資源的重要組成部分，具有供水、防洪、抗旱、發電、灌溉、水產養殖、娛樂等多種生態系統服務價值[2].Chao等通過統計全球水庫建設的歷史資料重構了人工水庫蓄水過程，計算發現由于水庫蓄水使得海平面上升趨勢減少了30 mm[3]，此外全球范圍內大量的水庫建設還會對水資源的時空分布格局產生直接的影響[4].因此，動態監測水庫水資源量的時空變化對于維護生態環境穩定和保障人類用水安全具有十分重要的意義.

水庫水體的水面范圍提取是水庫相關研究的基礎.水體范圍的變化能直觀反映水量的增減[5]，同時對水循環領域相關的其他應用研究也至關重要.例如，水庫蒸發的季節性/年際變化主要受水庫表面積的變化驅動[6]，水環境承載力與生態系統的穩定性也與水面積的變化密切相關[7].隨著對地觀測技術的迅猛發展，遙感技術以其大面積同步觀測、時效性強等優點，目前已經在水體范圍提取方面有了廣泛的應用[8].遙感技術根據傳感器的不同通常分為光學、熱紅外和微波遙感三個大類，其中光學遙感因其歷史數據資源豐富、提取精度高的特點，在長時序地表水識別中應用最為廣泛[9].

目前常用的光學遙感影像提取水體信息的方法包括監督與非監督分類[10]、水體指數閾值分割法[11]、以及決策樹法[12]等，其中閾值分割法應用范圍最為廣泛.張超等利用人工目視確定閾值對柴達木盆地湖泊水面變化進行了遙感提取[13]，吳川等、孫愛民等利用OTSU最大類間方差法確定閾值，分別對丹江口水庫和博斯騰湖水域變化進行了監測研究[14-15].

然而，現有研究中所采用的水體信息提取方法往往運行效率較低，大范圍、長時序的水體監測實現難度較大，此外，常規的光學遙感監測閾值方法在長時序、高精度監測上還存在一些明顯不足.首先，云及云陰影的影響難以消除[16].一些遙感影像上被云遮擋的地物信息無法獲取，在長時序水體提取中通常被舍棄，而未被遮擋區域所包含的有用信息沒有得到充分利用，因此水體監測結果也會存在大量的數據缺失.其次，大量水庫位于山區，山體陰影造成的水體錯分問題難以避免[17]，山體陰影會導致提取結果存在與水庫水體的混分現象，這些明顯的錯分增加了水體動態監測的誤差.此外，在水體的長時序監測中，海量的遙感影像下載及處理工作耗時費力，傳統的提取方法常常難以滿足監測快速響應的需求.

近年來遙感云計算技術的發展為海量遙感數據處理和分析提供了前所未有的機遇[18].谷歌地球引擎(google earth engine，GEE)是一個專門用于衛星遙感數據運算處理的云平臺，通過網頁瀏覽器編程運算，大大加快了算法的測試分析與大尺度的快速分析和應用的進程[19].針對全球尺度的水庫信息遙感提取，已有很多學者基于該平臺進行了相關研究工作.在全球地表水體制圖(joint research centre monthly water history， JRC MWH)數據集[20]的基礎上，Zhao等提出了針對水庫污染像元信息自動糾正的算法，生成了全球6 000多個水庫的面積時間序列[21]，該算法能充分提高遙感影像在水體提取中的利用率而被廣泛使用[22].然而，全球水體數據集往往更新周期長，且提取結果中存在著一些明顯的錯誤[23]，以JRC全球水體數據集為例，直至2021年初才更新前五年的數據，因此要想實現水庫水體范圍的動態監測，仍需探索一種簡單高效且準確的遙感監測方法.

鑒于此，本研究基于GEE平臺，結合Landsat影像和NASA DEM數據，提出了一種結合DEM和淹沒頻率的水庫水體遙感提取DF優化方法.針對全球九個不同的水庫實現了長時序水體動態提取，并結合精度驗證及測高水位數據集，與Zhao等基于JRC MWH數據集的水庫水體動態監測結果進行了對比，分析了新方法的精度與有效性.

1研究區及數據來源

1.1研究區域

全球水庫和大壩數據集(GRanD)[24]包含了全球范圍內7 250座水庫的范圍矢量數據及建成時間、大壩經緯度等基本信息.研究從中選取了九個不同地形和水文地質條件的水庫作為研究案例，用于測試本論文提出的研究方法.九個水庫分布在五大洲不同緯度區，高程分布從230 m到1 740 m.水庫具體分布如圖1所示，基本信息如表1所示.

圖1 九個研究水庫分布圖Fig.1 Distribution maps of nine reservoir

表1 水庫基本信息

1.2遙感數據

1.2.1Landsat影像 所采用的Landsat影像為USGS Landsat 5/7/8 Level 2，Collection 2，Tier 1影像數據集.該數據集包含了經過大氣校正的正射地表反射率(Landsat5/7包含三個可見光波段，一個近紅外波段以及兩個短波紅外波段；Landsat8包含四個可見光波段，一個近紅外波段以及兩個短波紅外波段)以及通過CFMASK算法[25]生成的云、雪卷標信息.數據空間分辨率均為30 m，重訪周期為16 d，獲取時間為1984年1月1日至2020年12月31日.

1.2.2DEM數據 選用的高程數據為NASA-DEM全球數據集，其空間分辨率也為30 m.該數據通過合并ASTER GDEM、ICESat GLAS和PRISM數據集的輔助對NASA SRTM DEM數據進行了再處理，顯著提高了其精度[26].其中最重要的改進處理在于使用了改進的相位解纏技術，通過ICESat GLAS數據進行控制來減少原有SRTM DEM數據中的數據空值區.

1.2.3JRC Monthly Water History地表水月歷史數據 Pekel等利用1984年3月16日至2020年12月31日期間采集的Landsat 5、7和8的400多萬幅影像，通過專家系統對每個像素點分別進行了水/非水分類，并將結果整理成月度歷史，用于動態監測全球水體變化[20].該數據集包含1984年3月至2020年12月間442個月的全球地表水的位置和時間分布的地圖，并提供了有關水體的范圍及其變化的統計數據.

2研究方法

利用Landsat Level2，Collection2，Tier1產品中的云、雪卷標去除影像中的受云、雪影響的無效地表像元，通過OTSU大津法對單一時相影像進行水庫水體信息的初步提取，在時序提取結果的基礎上統計每個像素的淹沒頻率.考慮到光學影像中無法回避的無效地表像元問題，根據單景影像提取結果中云、云影、雪等無效像元的長時序淹沒頻率以及同期正常像元的水體提取結果，修正異常像元的水體覆蓋信息.在此基礎上，結合修正后的水體淹沒頻率和DEM數據對每一景影像的提取結果進行優化.通過DEM和淹沒頻率約束剔除山體陰影及其他異常值的影響，從而最終獲得水庫水體的時序提取結果.工作流程如圖2所示.

圖2 處理流程圖Fig.2 Processing flow chart

2.1基于MNDWI的大津法水體初步提取

為了突出水體特征，便于大津法提取.首先通過波段運算，計算改進的歸一化差分水體指數(MNDWI)[27]，該指數計算公式為：

其中，XMNDWI為XNDWI計算值，Rgreen為綠光波段反射率(TM/ETM+的第2波段，OLI的第3波段)，RSWIR為短波紅外波段反射率(TM/ETM+的第5波段，OLI的第6波段).

改進的歸一化差分水體指數是在歸一化差分水體指數(NDWI)的基礎上利用短波紅外代替近紅外波段，該指數能更好的抑制建筑物和土壤的背景信息從而突出水體.

隨后基于MNDWI數據運用大津法自動化確定分割閾值從而完成水體和非水體的區分.大津法即最大類間方差法，是由Otsu于1979年提出的一種無參數無監督的閾值分割算法[28].其實現算法如下：

假設使用閾值T將灰度圖像分割為前景和背景；Size表示圖像總像素個數，u為圖像的平均灰度，W0為前景像素點占整幅圖像大小的比例，U0為前景像素點的平均值，W1為背景像素點占整幅圖像大小的比例，U1為背景像素點的平均值， 為類間方差，計算公式為：

u=W0×U0+W1×U1,

(1)

g=W0×(u-U0)2+W1×(u-U1)2.

(2)

將(1)式帶入(2)式得到類間方差計算式：

g=W0×W1(U0-U1)2.

(3)

該算法根據圖像的灰度特性，從圖像的最小灰度值遍歷到最大灰度值，找到使得分割的兩類類間方差最大，即g最大時的灰度值T即為最佳分割閾值[29].

2.2高質量數據篩選及淹沒頻率計算

由于光學遙感影像存在云雪覆蓋對地表地物的遮擋，以及Landsat7部分影像因傳感器故障造成了影像中大量的無效像元，研究通過云、雪卷標剔除被影響的無效像元，按有效地表覆蓋像元統計分析后將影像數據集分為三類：一類是缺失像元占比小于5%的遙感影像，將其視為高質量影像；第二類是缺失像元占比大于等于5%小于80%的遙感影像，這類影像存在有效信息可以利用，但需要對污染缺失像元進行信息的回填；第三類是缺失像元占比大于等于80%的影像，該類影像因有效信息參考太少，大津法分類誤差過大而被舍棄.在三類影像數據分類的基礎上，選擇高質量影像通過大津法的水庫水體提取結果計算水庫淹沒頻率分布，并用于無效像元的信息回填以及提取結果的優化.

2.3污染像元信息回填

Zhao等提出了基于淹沒頻率的水體無效像元回填算法[30]，該算法的原理是：通過單景影像中有效水體像元淹沒頻率的頻次直方圖，計算不同淹沒頻率在該景影像中有效水體像元中的出現頻次，當第一個頻次大于當期所有有效水體像元淹沒頻率下的頻次均值的0.17倍時，以該值所對應的淹沒頻率作為閾值.污染像元對應淹沒頻率高于該閾值的像元劃分為水體，反之則劃分為非水體.本研究基于此算法和計算的淹沒頻率，對缺失像元信息進行回填，以獲得更完整、連續的水體信息.

2.4水體DF優化提取結果

水庫水體提取修正基于以下假設：水庫水體應為連續水體，遠離主要水體的破碎斑塊應視為誤提.鑒于水庫水體的特殊性，提取的水體高程應小于壩頂高程，因此DEM超出壩頂高程的水像元應視為誤提；修正操作步驟為：

1) 以20%作為淹沒頻率閾值，得到淹沒頻率二值圖；

2) 以壩頂高程作為DEM閾值，得到DEM二值圖；

3) 提取及回填后同時滿足淹沒頻率大于20%和DEM小于壩頂高程的水像元被分類為水，否則被分為非水；篩選后水體柵格轉為矢量；

4) 選擇水體矢量的最大面要素作為水庫最終提取結果.

3結果與分析

3.1方法精度驗證

為驗證DF優化提取方法的精度，從九個水庫的多時相提取結果中各自隨機選擇了一景影像的提取結果，通過隨機生成的900個驗證樣點，驗證了分類結果的準確性，并與Zhao等基于JRC MWH數據集的直接修正結果[30]進行了對比分析.對比結果如表2所示，采用DF方法優化提取的水庫水面分類結果總體精度為95.67%，Kappa系數為0.912 6，而基于JRC MWH數據集的直接修正結果的總體精度為82.67%，Kappa系數為0.658 6.此外，圖3還分別展示了兩個數據的分類結果，可以看出基于JRC MWH數據集的直接修正結果中水庫水體存在明顯錯分現象.

表2 兩種方法水庫水體分類精度對比表Tab.2 Comparison of classification accuracy of reservoir water body by two methods

圖3 分類結果對比圖Fig.3 Comparison of classification results

在此基礎上，研究還計算了九個水庫建成以來至2020年12月的水體面積時間序列結果，并按月合成得到各月份的水面積均值.基于美國農業部官方網站提供的水庫水位衛星監測結果[31]，建立水面積與其同期水位監測結果的回歸關系，如圖4所示.從圖中可以看出，優化后水面積和水位之間有較好的一致性，決定系數(R2)均在0.6以上.

通過水面積與水位的相關關系，計算均方根誤差，結果如表3所示.實驗的九個水庫通過DF優化方法提取的水體面積，RMSE值均優于基于JRC MWH數據集的直接修正結果.總體RMSE均值從16.984 3減小到了11.796，其中五個水庫減小幅度超過了30%.由此也可證明本論文提出的水庫水體范圍優化提取方法具有更高的準確性.

圖4 水庫水面積與水位的回歸關系Fig.4 Regression relationship between water area and observed elevation of reservoirs

表3 各水庫提取結果精度評價

3.2水庫水體面積時間序列變化

利用上述算法，可以從遙感影像中準確地計算出水庫的水體面積，從而構建水庫水體面積的時序產品.然而，當水庫范圍內80%以上像元為無效像元時，DF優化方法無法使用.因此，使用該方法在部分月份無法獲取到水庫水體范圍信息.利用其他年份多年該月面積均值替代的方法對缺失的月份面積估計值進行補全.研究所提取的水庫面積時間序列結果如圖5所示.由圖5可以看出，從總體變化幅度上看，采用DF優化方法提取的水面積在1984年到2020年間波動幅度相較于直接修正的JRC MWH數據面積波動幅度更小，水庫水面積極高或極低的異常值出現的頻率也更低.

圖5 九個水庫月面積動態時間序列Fig.5 Dynamic time series of nine reservoirs per month

4討論

水庫水體范圍的遙感提取受到山體陰影干擾[32]、云的遮擋、雪覆蓋以及狹窄河道造成的水體不連續[33]等多種因素的制約，本文結合DEM和淹沒頻率，建立了DF優化方法消除部分山體陰影造成的水體錯分誤差，通過建立水體緩沖區消除斷裂水體，選擇最大水面要素作為水庫范圍，從而剔除遠離水庫的細小水體.DF方法在傳統大津法水體提取的結果上進行了優化，通過對比JRC MWH數據集的直接修正結果，驗證了DF優化提取方法的準確性.

當水陸像元占比接近時，大津法分割的結果最為理想[34]，但原始的影像邊界來自于GRanD數據庫shapefile矢量緩沖區，水陸占比與水庫的形狀及大小相關，難以控制確切的比例.因此，提取的最終結果也與大津法的分割效果有一定的關系.基于MNDWI大津法分割的水體提取對遙感影像的質量也有要求，當無效像元過多時，提取結果的可靠性降低，因此存在部分時段的影像被棄用的情況.這種影像棄用情況對于基于JRC MWH數據集的直接修正結果而言并不存在，因此其水庫水面積監測數據的時相理論上應更為完整.當然，也必須認識到針對無效值過多的影像基于JRC MWH數據集的直接修正結果雖然不會缺失，但其結果的可靠性也存在較大問題.

表4顯示了各水庫時序覆蓋數據以及月份缺失情況，其中污染像元占比超過80%的影像因無法準確進行大津法提取而被舍棄，水庫建成后至2020年間，有效月份數據占比均值為61.98%.基于JRC MWH數據集的直接修正結果中，無效值占比超過90%的月份被舍棄，有效覆蓋月份占比均值為69.48%，略高于本研究方法的提取結果.但是，由于JRC MWH數據集提取過程中算法的局部適宜性問題，也會出現部分水庫有效覆蓋數量明顯少于DF優化方法(如Tahtakopru水庫).

5結論

本論文提出了一種結合DEM和淹沒頻率的水庫水面積遙感DF優化提取方法，并以全球九個不同類型的水庫為研究區，對該方法進行了驗證，計算其1984年至2020年間的水面積動態變化過程,得出以下主要結論.

表4 水庫覆蓋數據統計Tab.4 Reservoir coverage data statistics

1) 相較于基于JRC MWH數據集的水庫地表水修正結果，利用DF優化提取方法得到的水體范圍具有更高的提取精度，且水庫水面積與水庫水位之間有更強的一致性，均方根誤差更小.

2) 受限于大津法分類原理，DF優化提取方法無法適用于污染像元占比過大的影像水體提取，但總體有效監測月數與已有基于JRC MWH數據集的水庫地表水修正結果中的有效數據占比相差不大.而且，缺失時相的影像數據質量本身無效值較多，其結果的有效性也值得商榷.同時，該提取方法不用依托于JRC MWH數據集的更新，動態監測能力更強.

3) 與傳統遙感監測方法相比，利用GEE云計算平臺能更高效的進行數據處理和統計分析，為后續大流域甚至全球尺度的水庫水體動態監測提供了有力支撐.