結(jié)合DEM與淹沒頻率的水庫水體動態(tài)遙感提取優(yōu)化方法

劉詩燕， 蔡曉斌*

(1.中國科學(xué)院大學(xué)， 北京 100049；2.中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院環(huán)境與災(zāi)害監(jiān)測評估湖北省重點實驗室， 武漢 430071)

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

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

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

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

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

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

1研究區(qū)及數(shù)據(jù)來源

1.1研究區(qū)域

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

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

表1 水庫基本信息

1.2遙感數(shù)據(jù)

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

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

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

2研究方法

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

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

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

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

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

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

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

假設(shè)使用閾值T將灰度圖像分割為前景和背景；Size表示圖像總像素個數(shù)，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)

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

2.2高質(zhì)量數(shù)據(jù)篩選及淹沒頻率計算

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

2.3污染像元信息回填

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

2.4水體DF優(yōu)化提取結(jié)果

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

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

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

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

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

3結(jié)果與分析

3.1方法精度驗證

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

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

圖3 分類結(jié)果對比圖Fig.3 Comparison of classification results

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

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

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

表3 各水庫提取結(jié)果精度評價

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

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

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

4討論

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

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

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

5結(jié)論

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

表4 水庫覆蓋數(shù)據(jù)統(tǒng)計Tab.4 Reservoir coverage data statistics

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

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

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