王詩(shī)蕾,羅 晉,陳澤強(qiáng)

1(武漢大學(xué) 遙感信息工程學(xué)院,武漢 430079)

2(武漢大學(xué) 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢 430079)

東洞庭湖國(guó)家級(jí)自然保護(hù)區(qū)位于長(zhǎng)江中游荊江江段南側(cè),是典型的濕地生態(tài)系統(tǒng),也是我國(guó)首批進(jìn)入“國(guó)際重要濕地公約”的濕地之一.作為洞庭湖湖系中最大的湖泊,東洞庭湖氣候溫和濕潤(rùn),孕育了鳥類、魚類、植物等豐富的物種資源,也是重要的漁業(yè)生產(chǎn)基地,對(duì)長(zhǎng)江流域的水位調(diào)控起重要作用.近年來,由于人類經(jīng)濟(jì)活動(dòng)和長(zhǎng)江中游地區(qū)洪旱災(zāi)害的影響,東洞庭湖的水文和水質(zhì)特征發(fā)生一定變化,生物多樣性受到威脅,調(diào)洪蓄水功能下降.水位-水體面積關(guān)系是描述湖區(qū)水情變化的指標(biāo)之一,對(duì)湖區(qū)生態(tài)環(huán)境的保護(hù)和政策制定有重要意義.

易波琳等[1]通過解譯岳陽和常德的TM 影像和計(jì)算機(jī)面積量算,再應(yīng)用最小二乘法原理,求得洞庭湖內(nèi)各湖面積與水位的回歸方程.彭定志等[2]利用NDVI對(duì)MODIS 影像進(jìn)行水體面積的計(jì)算,得到洞庭湖水位-面積曲線.杜濤等[3]利用2010年3 種質(zhì)量較好的MODIS 影像數(shù)據(jù),分別運(yùn)用NDVI和RVI 進(jìn)行水體提取,擬合出面積-水位曲線.黃一凡等[4]利用ENVI+IDL 4.7,對(duì)332 個(gè)時(shí)相的MOD13Q1 影像數(shù)據(jù)進(jìn)行批量分類處理,得到東洞庭湖水體濕地面積統(tǒng)計(jì)數(shù)據(jù),從而建立岳陽水文站和水體面積的回歸方程.柯文莉等[5]利用MRT (MODIS Reprojection Tool) 軟件,對(duì)MODIS L1 B數(shù)據(jù)進(jìn)行波段提取、幾何糾正、投影轉(zhuǎn)換,用ENVI進(jìn)行重采樣、影像裁剪等,得到研究區(qū)影像數(shù)據(jù),再利用NDWI 提取水體.眾多學(xué)者不斷改進(jìn)東洞庭湖水位-面積關(guān)系的研究方法,但是已有的研究在進(jìn)行水體提取時(shí)下載與處理數(shù)據(jù)繁瑣,需要使用許多專業(yè)軟件,且數(shù)據(jù)存儲(chǔ)量大,人工干預(yù)較多,耗時(shí)長(zhǎng),處理效率低,不適合長(zhǎng)時(shí)間尺度的地理分析和應(yīng)用[6].因此,需要提高水體提取的工作效率.

Google Earth Engine (GEE)是由Google 建立的處理遙感數(shù)據(jù)的云端計(jì)算平臺(tái),提供大量地球觀測(cè)數(shù)據(jù),具有極強(qiáng)的數(shù)據(jù)存儲(chǔ)和處理能力[7].用戶通過應(yīng)用編程接口API 可以進(jìn)行算法編寫、數(shù)據(jù)集生產(chǎn)和地理空間分析等.GEE 可以有效解決數(shù)據(jù)獲取困難、數(shù)據(jù)處理復(fù)雜耗時(shí)等問題,提高用戶的工作效率.

鑒于此,本文基于Google Earth Engine,選擇2001年至2016年的Landsat 遙感影像,利用NDWI 模型提取東洞庭湖水體面積,根據(jù)城陵磯水文站實(shí)測(cè)水位信息,建立東洞庭湖水位-水體面積的關(guān)系曲線.

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域簡(jiǎn)介

東洞庭湖(29°00′ N 至29°38′ N,112°43′ E 至113°14′ E)位于湖南省北部岳陽市,地處長(zhǎng)江中游荊江南側(cè),屬于亞熱帶濕潤(rùn)季風(fēng)氣候,擁有充足的日照,降水量年際變化大.東洞庭湖是洞庭湖湖區(qū)中保存最為完好的季節(jié)性湖泊,也是長(zhǎng)江上重要的調(diào)蓄湖泊[2].湖區(qū)設(shè)有湖南東洞庭湖國(guó)家級(jí)自然保護(hù)區(qū),1992年加入《國(guó)際重要濕地公約》,1994年升級(jí)為國(guó)家級(jí)自然保護(hù)區(qū),是我國(guó)典型的濕地生態(tài)環(huán)境.城陵磯水文站是長(zhǎng)江流域135 個(gè)水文站之一,位于洞庭湖與長(zhǎng)江的交匯處,對(duì)長(zhǎng)江流域防汛測(cè)報(bào)起重要作用.圖1為東洞庭湖區(qū)域Landsat 影像合成圖.

圖1 東洞庭湖區(qū)域

1.2 Google Earth Engine 平臺(tái)簡(jiǎn)介

Google Earth Engine (GEE)又稱谷歌地球引擎,是一個(gè)基于云操作的平臺(tái),支持快速、交互式探索和分析空間數(shù)據(jù).GEE 包含大量公開的地理空間數(shù)據(jù)集,除了Landsat、MODIS、Sentinel 等遙感影像數(shù)據(jù)外,還提供地形數(shù)據(jù)、土地覆蓋數(shù)據(jù)、氣象數(shù)據(jù)、地球物理數(shù)據(jù)和社會(huì)經(jīng)濟(jì)數(shù)據(jù)等.用戶通過應(yīng)用程序編程接口(API)和基于Web 的交互式開發(fā)環(huán)境(IDE)進(jìn)行訪問控制,可使用Python和JavaScript 語言編寫和調(diào)用算法,利用該平臺(tái)的高性能計(jì)算資源處理大型空間數(shù)據(jù)集[7].

Google Earth Engine 主要從以下方面提高用戶工作效率:(1)數(shù)據(jù)易得:用戶不需要自行采集、存儲(chǔ)和管理遙感數(shù)據(jù),可直接調(diào)用平臺(tái)提供的地理數(shù)據(jù)集;(2)數(shù)據(jù)處理方便:平臺(tái)提供圖像處理、機(jī)器學(xué)習(xí)和地理統(tǒng)計(jì)等操作,用戶通過函數(shù)式編程環(huán)境對(duì)地理信息進(jìn)行分析;平臺(tái)提供大量CPU 的編組和管理,用戶無需提供硬盤和服務(wù)器;平臺(tái)利用Java Just-In-Time (JIT)編譯器優(yōu)化逐像素操作鏈的執(zhí)行;(3)結(jié)果共享:平臺(tái)可方便用戶將研究結(jié)果傳播給其他研究人員、政府工作者和普通公眾等.

1.3 實(shí)驗(yàn)數(shù)據(jù)及處理

本文數(shù)據(jù)有兩類:(1) 影像數(shù)據(jù):選用2002年至2016年東洞庭湖區(qū)域的Landsat 遙感影像作為原始影像數(shù)據(jù),包括Landsat 5 影像63 景,Landsat 8 影像34景,一共97 景影像.利用Google Earth Engine 平臺(tái)進(jìn)行數(shù)據(jù)下載,編程完成輻射定標(biāo)、幾何校正、大氣校正、去云處理、影像拼接、掩模去除建筑區(qū)等預(yù)處理.由于東洞庭湖常年處于多云天氣,很難獲得連續(xù)的無云影像,因此需要將不同時(shí)期的少云影像進(jìn)行拼接,最后得到處理后的2002年至2016年東洞庭湖影像數(shù)據(jù)36 景,每一景為3 個(gè)月的拼接合成影像.(2)水位數(shù)據(jù):城陵磯水文站水位數(shù)據(jù)來源于湖南省水利廳水文信息網(wǎng),包括水位采集時(shí)間、逐日水位(以米為單位).

2 研究方法

2.1 基于NDWI 的水體提取

常用的水體提取方法有單波段閾值法、譜間關(guān)系法和指數(shù)法等.單波段閾值法在TM5 波段影像上選擇一個(gè)閾值進(jìn)行水體提取,原理簡(jiǎn)單,但水體與非水體的交界區(qū)域提取效果差[8].譜間關(guān)系法利用多個(gè)波段影像信息,能有效區(qū)分水體和陰影,適合地勢(shì)起伏較大的地區(qū),而在平原地區(qū)優(yōu)勢(shì)不明顯[9].指數(shù)法根據(jù)研究地類在各波段反射率的強(qiáng)弱,通過比值運(yùn)算突出感興趣地物.

NDWI(Normalized Difference Water Index)指數(shù)對(duì)湖泊等大塊水體的提取效果好,可以突出水體特征,抑制非水體特征[10].EI-Asmar 等[11]利用NDWI研究尼羅河三角洲水體面積.李文波等[12]利用NDWI從ETM+遙感數(shù)據(jù)中提取水體.根據(jù)徐蓉等[13]采用3 種方法對(duì)我國(guó)不同地區(qū)的10 個(gè)湖泊進(jìn)行水體提取的結(jié)果,NDWI比波段插值模型適合東部地區(qū)湖泊提取.根據(jù)張飛等[14]對(duì)新疆艾比湖水體提取的結(jié)果,NDWI比MNDWI精度更高.對(duì)于內(nèi)陸湖泊,利用NDWI指數(shù)水體提取效果好,因此選用NDWI進(jìn)行水體提取.

NDWI利用水體在綠光反射率最大、水體在近紅外波段反射率小、植被土壤在近紅外波段反射率大的特征,根據(jù)式(1)在包含綠波段和近紅外波段的影像中提取水體[15].一般而言,當(dāng)NDWI指數(shù)為正值時(shí),判斷地物為水體;當(dāng)NDWI指數(shù)為零值或負(fù)值時(shí),判斷為植被或土壤.實(shí)際情況中水體含有葉綠素、懸浮物等物質(zhì),因此要根據(jù)研究區(qū)域的具體情況設(shè)定合適的閾值.

式中,GREEN為L(zhǎng)andsat 數(shù)據(jù)在綠波段的反射率,NIR為L(zhǎng)andsat 數(shù)據(jù)在近紅外波段的地表反射率.

確定閾值是利用NDWI區(qū)分水體和非水體的重要步驟,但是水體和其他物體之間的閾值不穩(wěn)定,且隨著場(chǎng)景和位置的變化而變化,因此為每個(gè)圖像手動(dòng)尋找閾值是不切實(shí)際的.Otsu 閾值被廣泛用于確定NDWI在水體檢測(cè)的閾值[16].Otsu 方法從圖像直方圖中確定理想閾值,使類間方差最大,類內(nèi)方差最小.選擇Otsu方法自動(dòng)確定每個(gè)月的閾值,共得到12 個(gè)閾值(見表1),根據(jù)閾值判斷地面像元是否為水體,減少因不同時(shí)期水文特征不同造成水體信息提取遺漏現(xiàn)象.當(dāng)像元NDWI大于閾值時(shí),判斷該像元為水體,當(dāng)像元NDWI小于閾值時(shí),判斷該像元為非水體.由于預(yù)處理后的Landsat影像為拼接影像,包含兩個(gè)或兩個(gè)以上日期的水位信息,因此一幅影像上使用的NDWI閾值同樣多于一個(gè).利用Google Earth Engine 平臺(tái),完成不同閾值的水體提取,得到東洞庭湖地區(qū)水體提取結(jié)果的二值化圖像36 幅.

表1 NDWI 閾值

2.2 水位與水體面積關(guān)系建立方法

本文使用Landsat 數(shù)據(jù)作為影像數(shù)據(jù),由于Landsat影像空間分辨率為30 m,因此利用單個(gè)像元的面積(30 m×30 m=900 m2=0.0009 km2)乘水體像元個(gè)數(shù)可計(jì)算出東洞庭湖水體的總面積.根據(jù)預(yù)處理后得到的Landsat 影像日期信息,可以找到每幅影像對(duì)應(yīng)的城陵磯水文站水位信息.

根據(jù)原始Landsat 影像的日期,將東洞庭湖水體面積序列與對(duì)應(yīng)的水位值序列進(jìn)行配對(duì),得到水位-水體面積組.由于一幅影像對(duì)應(yīng)多個(gè)水位值(1～n),因此在進(jìn)行水位-水體面積函數(shù)關(guān)系的確定時(shí),設(shè)計(jì)3 種方法得到水體提取結(jié)果對(duì)應(yīng)的最佳水位值,再通過線性函數(shù)、指數(shù)函數(shù)、對(duì)數(shù)函數(shù)、多項(xiàng)式函數(shù)等方式建立最優(yōu)水位-水體面積關(guān)系曲線.

第1 種方法為水位加權(quán)平均法.水體提取的結(jié)果取決于該區(qū)域不同時(shí)期的拼接影像,拼接區(qū)域的大小直接影響提取結(jié)果,因此不同時(shí)期的水位對(duì)該影像水位信息的貢獻(xiàn)有差異.基于對(duì)各日期水位值貢獻(xiàn)率的綜合考慮,將水體影像中各個(gè)水位值(即各個(gè)日期的影像區(qū)域)所占面積的比重作為權(quán)重,得到加權(quán)平均的水位值(如式(2))作為最佳水位.

式中,lj為水體影像對(duì)應(yīng)的水位值,aj為該水位值日期的水體面積,a為水體總面積,L為水體影像的加權(quán)平均水位值.

第2 種方法為單點(diǎn)水位法.由于城陵磯水文站位于東洞庭湖與長(zhǎng)江的交界處,同時(shí)也處于水體提取的邊界區(qū)域,相對(duì)于東洞庭湖其他區(qū)域,該水文站所在像元的水位值對(duì)水體提取的結(jié)果貢獻(xiàn)更大.因此,選用該像元水位代表東洞庭湖區(qū)域的水位值.首先確定影像中城陵磯水文站所在位置的影像日期,再選取該日期的水位值為最優(yōu)水位.

第3 種方法為閾值法.東洞庭湖水情特征比較復(fù)雜,按照水情的變化可以分為漲水期、退水期、豐水期和枯水期,每個(gè)時(shí)期的平均水位值差距較大[4].對(duì)于不同水文時(shí)期的東洞庭湖區(qū)域,采用不同的標(biāo)準(zhǔn)來確定最優(yōu)水位.為便于分析,將東洞庭湖水文時(shí)期分為兩類:豐水期和枯水期.豐水期水位值相對(duì)更大,枯水期水位值相對(duì)更小.根據(jù)水位-水體面積的統(tǒng)計(jì)信息,確定一個(gè)面積閾值將水位-水體面積序列分為兩類.如果水體面積大于該閾值,視為豐水期,選取影像最大水位值作為最優(yōu)水位;如果面積小于該閾值則視為枯水期,則取影像最小水位值作為最優(yōu)水位.

2.3 精度評(píng)價(jià)

對(duì)于水位-水體面積函數(shù)關(guān)系,使用確定性系數(shù)R2(見式(3))進(jìn)行擬合方程的精度評(píng)定.

式中,為由擬合方程解算的水體面積,ai為水體面積真值(即由水體提取計(jì)算得到的水體面積),為水體面積真值的平均值.

3 結(jié)果與討論

3.1 東洞庭湖水域水體提取結(jié)果

Landsat 影像預(yù)處理后得到36 景東洞庭湖影像數(shù)據(jù),根據(jù)選取的NDWI閾值,利用Google Earth Engine提取水體.

利用定量統(tǒng)計(jì)分析方法,隨機(jī)選取其中4 景影像的水體提取結(jié)果作為樣本進(jìn)行精度評(píng)價(jià).首先,基于Google Earth Engine 在影像上選取300 個(gè)隨機(jī)點(diǎn),其次,將原始影像與NDWI提取結(jié)果進(jìn)行疊加分析,通過混淆矩陣法計(jì)算水體提取精度(如表2).

表2 水體提取結(jié)果與精度評(píng)價(jià)

從表2可看出,水體提取結(jié)果的總體精度都超過95%,Kappa 系數(shù)大于0.9,說明NDWI進(jìn)行水體提取的結(jié)果理想,利用Ostu 方法得到的NDWI閾值可以準(zhǔn)確提取東洞庭湖水體面積.

從圖2提取結(jié)果可知,隨著季節(jié)的變化,東洞庭湖水體面積變化較大,最小面積常出現(xiàn)在春季(1～3月),夏季(4～6月)開始上升,到秋季(7～9月)時(shí)湖面面積達(dá)到一年的最大值,在冬季(10～12月)水面積逐漸減小.研究時(shí)間內(nèi)的最大面積為1204.7517 km2(2002年7～9月),最小面積為240.7221 km2(2008年1～3月).春季平均面積為332.6463 km2,夏季平均面積為715.1662 km2,秋季平均面積為1078.1585 km2,冬季平均面積為696.5282 km2.由提取影像面積變化可以看出,東洞庭湖是一個(gè)季節(jié)性的湖泊.

圖2 東洞庭湖水域水體提取結(jié)果

3.2 東洞庭湖水位與水體面積關(guān)系

使用水位加權(quán)平均法、單點(diǎn)水位法和閾值法(根據(jù)統(tǒng)計(jì)值確定閾值為750 km2)這3 種方法對(duì)每幅水體提取影像進(jìn)行計(jì)算,得到最佳水位-水體面積序列,如表3所示.由計(jì)算可知,加權(quán)平均水位法與閾值法的相關(guān)系數(shù)為88.67%,單點(diǎn)水位法與閾值法相關(guān)系數(shù)為90.80%,加權(quán)平均水位法與單點(diǎn)水位法的相關(guān)系數(shù)為86.11%.3 種水位序列相關(guān)性不大于95.00%,具有一定差異性.

表3 最佳水位-水體面積序列

利用線性函數(shù)、指數(shù)函數(shù)、對(duì)數(shù)函數(shù)、三次多項(xiàng)式、五次多項(xiàng)式這5 種函數(shù)模型,以最優(yōu)水位值為橫坐標(biāo),東洞庭湖水體面積為縱坐標(biāo),建立最優(yōu)水位-水體面積關(guān)系曲線.由關(guān)系曲線(如圖3)可知,水體面積隨水位增長(zhǎng)整體呈上升趨勢(shì).

圖3 東洞庭湖水位-水體面積線性模型對(duì)比圖

使用3 種方法求解得到5 種模型的確定性系數(shù)R2如表4所示.由表4可看出,閾值法擬合的水位-水體面積函數(shù)模型精度最好,大部分模型優(yōu)于其他方法的模型,且5 種模型的確定性系數(shù)均大于90%.分析原因,由于水體影像由同一季節(jié)若干幅不同日期的影像拼接而成,在3 個(gè)月內(nèi)水位變化可能較大(如2016年9月25日與2016年7月23日水位差為10.07 m),如使用加權(quán)平均水位值或單點(diǎn)水位法則不能準(zhǔn)確表達(dá)該季節(jié)水位情況.使用閾值法選擇最優(yōu)水位可在一定程度上避免此類情況.

表4 3 種方法的確定性系數(shù)(R2)

對(duì)比閾值法的5 種模型,可知五次多項(xiàng)式模型(如式(4))確定性系數(shù)最高,其次是三次多項(xiàng)式模型、線性函數(shù)模型、對(duì)數(shù)函數(shù)模型,指數(shù)函數(shù)模型確定性系數(shù)最低.

其中,l為東洞庭湖水位值,a為東洞庭湖水域面積.

由五次多項(xiàng)式模型擬合的水位-水體面積函數(shù)關(guān)系曲線(如圖4)可知,88.89%的水體面積擬合誤差小于50 km2,94.44%的水體面積擬合誤差小于100 km2,五次多項(xiàng)式模型擬合結(jié)果較好.水位值在23 m 至29 m之間時(shí)擬合效果最好,擬合誤差較大的時(shí)期集中在秋季(7～9月)和春季(1～3月).該模型對(duì)水體面積在250 km2至1200 km2范圍外的水位-面積序列擬合誤差較大,不適用于水體面積過大或過小的湖區(qū)情況.

圖4 閾值法五次多項(xiàng)式模型擬合方程(R2=0.9628)

與前人的研究對(duì)比,本文水位-水體面積關(guān)系建立方法主要有3 點(diǎn)優(yōu)勢(shì):

(1)探索合成影像上多個(gè)水位值與水體面積的關(guān)系.由于拼接影像由多個(gè)不同時(shí)間采集的影像合成,對(duì)應(yīng)多個(gè)水位值,現(xiàn)研究主要選取無云影像或質(zhì)量較好的影像,很少對(duì)拼接影像進(jìn)行研究,因此實(shí)際可用的影像數(shù)量較少,直接影響水位-水體面積模型精度.如何將拼接影像的多個(gè)水位值與獲取的水域面積建立關(guān)系,本文進(jìn)行了探索,設(shè)計(jì)并討論3 種水位選取方法,提高了研究時(shí)間范圍內(nèi)遙感影像的利用率,同時(shí)在一定程度上避免了拼接影像水位選取的誤差.

(2)影像數(shù)據(jù)量大,時(shí)間尺度大.受到影像質(zhì)量和處理效率的限制,現(xiàn)研究多采用小尺度時(shí)間序列(某年或數(shù)年內(nèi))的若干幅遙感影像進(jìn)行分析,水位-水體關(guān)系模型樣本少,且存在時(shí)間上的局限性,無法滿足對(duì)長(zhǎng)時(shí)間序列關(guān)系研究的要求.本文基于Google Earth Engine高效處理并得到較大數(shù)據(jù)量的影像數(shù)據(jù),所建立的水位-水體關(guān)系模型時(shí)間尺度大,為東洞庭湖區(qū)域長(zhǎng)期水資源變化研究和管理提供依據(jù).

(3)模型精度高.柯文莉等[4]對(duì)2000 至2012年東洞庭湖水面面積與城陵磯水位月平均繩套曲線進(jìn)行研究,得到R2為93.7%.宋求明等[17]利用二次多項(xiàng)式對(duì)2003 至2006年洞庭湖水位-面積關(guān)系進(jìn)行研究,得到R2為91.32%.本文利用閾值法選取影像最優(yōu)水位值,利用五次多項(xiàng)式模型得到R2為96.28%,說明該方法更加準(zhǔn)確地反映實(shí)際水位與湖泊水體面積關(guān)系.

4 結(jié)語

本文基于Google Earth Engine 平臺(tái),利用長(zhǎng)時(shí)間序列的東洞庭湖區(qū)域Landsat 歷史影像,使用NDWI閾值法提取水體得到東洞庭湖水體面積,探索影像所對(duì)應(yīng)的最優(yōu)水位值,確定水位-水體面積關(guān)系曲線模型,得到以下結(jié)論:

(1)Google Earth Engine 平臺(tái)可以高效下載并處理大量的遙感數(shù)據(jù),具有很強(qiáng)的計(jì)算能力.

(2)根據(jù)東洞庭湖的水文特征,利用Otsu 方法確定利用12 個(gè)月閾值分別為0.08、0.06、0.03、0.01、0.00、0.00、0.00、0.00、0.01、0.02、0.05、0.08 的NDWI值進(jìn)行水體提取和面積計(jì)算.東洞庭湖水域面積隨季節(jié)變化大,冬夏差異明顯,7～9月水域面積最大(平均面積為1078.1585 km2),1～3月水域面積最小(平均面積為332.6463 km2).研究時(shí)期內(nèi)最大面積為1204.752 km2,最小面積為240.722 km2.

(3)設(shè)計(jì)3 種方法確定影像水位,經(jīng)實(shí)驗(yàn)得到利用閾值法計(jì)算影像最優(yōu)水位值的效果最好(最大R2為0.9628),優(yōu)于加權(quán)平均水位法(最大R2為0.8322)和單點(diǎn)水位法(最大R2為0.9457).對(duì)于閾值法計(jì)算的最優(yōu)水位值,通過五次多項(xiàng)式模型建立的水位-水體面積函數(shù)關(guān)系精度最高(R2=0.9628),優(yōu)于線性函數(shù)(R2=0.9317)、指數(shù)函數(shù)(R2=0.908)、對(duì)數(shù)函數(shù)(R2=0.9287)和三次多項(xiàng)式(R2=0.9566).