彭燕，何國金，王桂周，尹然宇

1. 中國科學院空天信息創(chuàng)新研究院，北京 100094

2. 中國科學院大學，北京 100049

3. 海南省地球觀測重點實驗室，海南三亞 572029

引 言

2016 年3 月23 日，瀾滄江-湄公河區(qū)域六國（中國、緬甸、老撾、柬埔寨、越南、泰國）領(lǐng)導人在三亞舉行“六水合一”儀式，正式啟動瀾湄合作，該合作是因水而生，水資源對于瀾湄流域的重要性不言而喻。瀾湄流域?qū)儆谖髂霞撅L氣候區(qū)，6 月至11 月為濕季，12 月至次年5 月為干季，年降水量較高，但時空分布 不均，上游與下游間的氣候及水資源分布特征存在較大差異，瀾滄江流域降水集中在6 月至8 月，而湄公河流域降水峰值集中在9 月和10 月[1]。近期以來，受強厄爾尼諾現(xiàn)象的影響，越南南部、柬埔寨等湄公河流域遭受了嚴重的干旱襲擊。因此，瀾湄流域30 m 分辨率的陸表水體產(chǎn)品能反映瀾湄流域陸表水體空間分布情況以及變化軌跡，在瀾湄流域水資源安全和管理、氣候變化、生態(tài)環(huán)境等方面具有重要的應(yīng)用價值。清華大學、國家基礎(chǔ)地理信息中心、美國馬里蘭大學、歐盟聯(lián)合研究中心等均發(fā)布了全球陸表水體產(chǎn)品，時間從1985 年到2015 年，空間分辨率從25 km 到30 m 不等，所采用的方法大多為決策樹、指數(shù)閾值法等方法。而目前尚缺少2018 年瀾湄流域的陸表水體產(chǎn)品。

本文提供2018 年30 m 分辨率的瀾湄流域年度陸表水體產(chǎn)品，給出了瀾湄流域30 m 陸表水體產(chǎn)品生產(chǎn)的方法與技術(shù)流程。采用基于多指數(shù)和先驗知識驅(qū)動的全球尺度遙感水體智能提取方法，首先對所需要的Landsat 8 地表反射率數(shù)據(jù)集、GDEM（ASTER Global Digital Elevation）數(shù)據(jù)集以及GLCF-GIW（2000）（the Global Land Cover Facility Global Inland Water）水體產(chǎn)品數(shù)據(jù)集進行數(shù)據(jù)預處理，隨后根據(jù)GLCF-GIW（2000）水體產(chǎn)品進行水體樣本與非水體樣本選擇與優(yōu)化，最后采用決策樹分類法進行自適應(yīng)閾值確定，從而得到初步的水體提取結(jié)果，最后再進行空間過濾等后處理得到最終的陸表水體產(chǎn)品。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)采集方法

本文所生產(chǎn)的瀾湄流域季度陸表水體產(chǎn)品主要是在Landsat 8 地表反射率的基礎(chǔ)上進行生產(chǎn)的，地表反射率數(shù)據(jù)主要來自于中國科學院空天信息創(chuàng)新研究院何國金研究員團隊所生產(chǎn)的中國Landsat 系列衛(wèi)星遙感數(shù)據(jù)地表反射率產(chǎn)品數(shù)據(jù)集[2]。所需要的DEM 數(shù)據(jù)來自于GDEM Version 2.0 dataset[3]；并采用 2000 年 GLCF-GIW version 1.0 水體產(chǎn)品[4]作為樣本庫，該數(shù)據(jù)集可從http://www.landcover.org/data/watercover/免費獲取。

1.2 數(shù)據(jù)處理方法

采用基于多指數(shù)和先驗知識驅(qū)動的全球尺度遙感水體智能提取方法來提取瀾湄流域季度陸表水體專題信息，其技術(shù)路線如圖1 所示，主要分為數(shù)據(jù)預處理、樣本選擇與優(yōu)化、自適應(yīng)確定閾值、連通像元合并及空間過濾后處理4 個步驟。

1.2.1 數(shù)據(jù)預處理

需要對Landsat 8 地表反射率、DEM 數(shù)據(jù)集以及已有的30 m 水體產(chǎn)品進行數(shù)據(jù)預處理。首先，對于2018 年的所有地表反射率產(chǎn)品根據(jù)質(zhì)量評估文件（QA）進行去云處理，將所有Landsat 8 地表反射進行中值合成，得到2018 年瀾湄流域的中值合成影像，最終為了便于計算，將瀾湄流域的年度中值合成影像按5120×5120 pixels 的大小分塊，共計約179 塊。圖2 為瀾湄流域Landsat 8 地表反射率中值合成影像的分塊示意圖。然后，將全球范圍內(nèi)的DEM 數(shù)據(jù)集鑲嵌成VRT 格式。由于2000 年GLCF-GIW 的水體產(chǎn)品是以景為單位存儲的，坐標系投影為WGS84 UTM（Universal Transverse Mercator Projection，通用橫軸墨卡托投影），為了能適應(yīng)于瀾湄流域甚至是全球區(qū)域的陸表水體信息提取，需要將該水體產(chǎn)品處理成VRT 格式的鑲嵌結(jié)果。

圖1 瀾湄流域陸表水體生產(chǎn)技術(shù)路線圖

圖2 瀾湄流域Landsat 8 地表反射率中值合成影像分塊示意圖

1.2.2 樣本選擇與優(yōu)化

根據(jù)2000 年的GLCF-GIW水體產(chǎn)品對每一分塊影像分別分層隨機選取水體與非水體樣本各150個。由于所利用的訓練樣本數(shù)據(jù)為2000 年的水體產(chǎn)品，然而水體會隨著時間和季節(jié)的變化而變化，難免會出現(xiàn)樣本數(shù)據(jù)為水體，而待提取的數(shù)據(jù)為非水體（如植被）的現(xiàn)象。因此為了避免出現(xiàn)這種因樣本不準確而導致誤分的問題出現(xiàn)，需要對樣本進行進一步優(yōu)化。利用坡度（Slope）和山體陰影（Hillshade）兩種地形指數(shù)以及NDVI（Normalized Difference Vegetation Index，歸一化差值植被指數(shù)）、NDWI（Normalized Difference Water Index，歸一化差值水體指數(shù)）對所選取的樣本進行基于先驗知識驅(qū)動的樣本優(yōu)化。NDVI 指數(shù)用于建立植被掩膜，為避免因時相變化引起的訓練樣本不完全正確時，將植被誤認為水體的問題。NDWI 主要用于優(yōu)化非水體樣本中存在的水體樣本的問題。并將Hillshade＜150 且Slope＞20 的區(qū)域認為是山體陰影[2]，建立山體陰影掩膜，優(yōu)化水體樣本中存在山體陰影的問題。

1.2.3 自適應(yīng)確定閾值

采用了馬里蘭大學發(fā)布的2000 年的水體產(chǎn)品（GLCF-GIW）作為先驗知識選取水體與非水體樣本，并將slope 和hillshade 等地形指數(shù)、NDVI、NDWI、NDBI（Normalized Difference Build Index，歸一化差值建筑指數(shù)）作為輸入，采用決策樹方法對優(yōu)化后的樣本進行了樣本訓練，自適應(yīng)確定分類規(guī)則，該分類規(guī)則便是所選擇指數(shù)的一個閾值組合，從而根據(jù)該分類規(guī)則得到一個初步的陸表水體專題信息提取結(jié)果。

1.2.4 后處理

由于基于像元的分類方法所得到的結(jié)果往往會出現(xiàn)獨立像元或者一兩個像元的孔洞現(xiàn)象，因此需要將得到的初步的陸表水體專題信息提取結(jié)果進行聯(lián)通像元合并及空間過濾。初始結(jié)果是分塊的，因此需要將分塊結(jié)果進行鑲嵌，從而得到最終的瀾湄流域30 m 空間分辨率陸表水體專題信息產(chǎn)品。

2 數(shù)據(jù)樣本描述

2018 年瀾湄流域30 m 陸表水體產(chǎn)品數(shù)據(jù)集包括1 個結(jié)果文件和1 個相應(yīng)的影像成像時間文件。陸表水體產(chǎn)品文件的命名規(guī)則為：water-lanmei-yyyy.TIF（如：water-lanmei-2018.TIF）。影像的空間分辨率為0.00025°（約30 m），投影坐標系是WGS84 經(jīng)緯度。陸表水體產(chǎn)品結(jié)果為灰度二值圖，其中1 表示水體，0 表示非水體。為了降低存儲空間，對柵格結(jié)果進行了“LZW”的無損壓縮。圖3 為2018 年瀾湄流域陸表水體產(chǎn)品示意圖，底圖采用的是2018 年Landsat 地表反射率合成圖，波段組合為R(6)G(5)B(4)。對應(yīng)的影像成像時間文件的命名規(guī)則為waterdate-lanmei-yyyy.TIF（如waterdatelanmei-2018.TIF），為16 位整型灰度圖像，對應(yīng)的灰度值表示該年的儒略日，如18 表示為該像元采用2018 年1 月18 日的Landsat 8 影像進行制圖而成。

3 數(shù)據(jù)質(zhì)量控制和評估

圖4 給出了本數(shù)據(jù)集水體提取精度驗證的抽樣分布情況。為了對本方法的水體提取結(jié)果進行驗證，在瀾湄流域進行分層隨機選取樣本點，水體和非水體各約2000 個。以原始影像結(jié)合GoogleEarth、GF-1/2 等高空間分辨率遙感影像作為參考影像，進行精度驗證。表1 為水體提取結(jié)果精度驗證混淆矩陣，水體的制圖精度達99.34%，用戶精度達98.25%，總體精度達98%以上。

陸表水體信息提取的難點在于高建筑物以及山體陰影的影響，因此為了驗證本數(shù)據(jù)集在此方面的表現(xiàn)，將本數(shù)據(jù)集與國際上已有的陸表水體產(chǎn)品進行交叉驗證。由于目前國際上尚無2018 年瀾湄流域的陸表水體產(chǎn)品，因此將本數(shù)據(jù)集分別與清華大學發(fā)布的2010 年30 m 全球土地覆蓋產(chǎn)品中的水體類別(FROM-GLC water mask)[5]以及美國馬里蘭大學發(fā)布的2000 年30 m 全球陸表水體產(chǎn)品（GLCF-GIW）[6]進行對比分析，如圖5 所示。分別選取了城區(qū)和山區(qū)兩處的陸表水體結(jié)果進行對比，具體位置如圖5 中瀾湄流域示意圖上標識的1 處和2 處。位置1 處各產(chǎn)品的局部放大圖如圖5(af)所示，位置2 處各產(chǎn)品的局部放大圖如圖5(g-h)所示。為了避免由于時相引起的水體變化，圖5 分別將各產(chǎn)品所對應(yīng)的Landsat 數(shù)據(jù)也展示出來。從圖5 可以看出，位置1 處，F(xiàn)ROM-GLC water mask產(chǎn)品和GLCF-GIW 產(chǎn)品均存在將城區(qū)誤分成水體的現(xiàn)象，本數(shù)據(jù)集較好地區(qū)分了建筑物陰影與水體（如圖5(f)所示）；位置2 處，F(xiàn)ROM-GLC water mask 產(chǎn)品和GLCF-GIW 產(chǎn)品均存在將山體陰影誤分成水體的現(xiàn)象，本數(shù)據(jù)集較好地區(qū)分了山體陰影與水體（如圖5(l)所示）。以上結(jié)果均表明本數(shù)據(jù)集水體提取精度較高。

圖3 2018 年瀾湄流域陸表水體產(chǎn)品示意圖

表1 水體提取結(jié)果精度驗證混淆矩陣

圖4 瀾湄流域季度陸表水體產(chǎn)品精度驗證樣本分布情況

圖5 本數(shù)據(jù)集與國際上已有陸表水體產(chǎn)品的對比驗證圖

4 數(shù)據(jù)價值

本文推出2018 年瀾湄流域30 m 分辨率陸表水體產(chǎn)品，采用基于多指數(shù)和先驗知識驅(qū)動的全球尺度遙感水體智能提取方法，精度較高，后續(xù)本數(shù)據(jù)集將會補充其他年份的季度產(chǎn)品。本數(shù)據(jù)集在瀾湄流域水資源管理與災(zāi)害應(yīng)急、生態(tài)環(huán)境監(jiān)測等方面具有重大的應(yīng)用價值。

致 謝

衷心感謝劉慧嬋和江威在產(chǎn)品質(zhì)量檢驗時給予的建設(shè)性意見。

數(shù)據(jù)作者分工職責

彭燕（1988—），女，湖南郴州市人，在讀博士，工程師，研究方向為遙感圖像智能處理。主要承擔工作：算法集成程序編寫，數(shù)據(jù)生產(chǎn)流程設(shè)計，論文撰寫。

何國金（1968—），男，福建龍巖人，博士，研究員，研究方向為遙感數(shù)據(jù)智能處理與信息挖掘。主要承擔工作：總體思路與方案設(shè)計，論文修改。

王桂周（1984—），男，山東省濟寧市人，博士，高級工程師，研究方向為遙感圖像智能處理。主要承擔工作：技術(shù)指導。

尹然宇（1996—），男，山東省臨沂人，在讀博士，研究方向為遙感圖像智能處理。主要承擔工作：數(shù)據(jù)挑選、整合與預處理。