基于Google Earth Engine的湖泊水位與水體面積關系研究①

王詩蕾,羅 晉,陳澤強

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

2(武漢大學 測繪遙感信息工程國家重點實驗室,武漢 430079)

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

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

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

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

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

1.1 研究區(qū)域簡介

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

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

1.2 Google Earth Engine 平臺簡介

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

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

1.3 實驗數據及處理

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

2 研究方法

2.1 基于NDWI 的水體提取

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

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

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

式中,GREEN為Landsat 數據在綠波段的反射率,NIR為Landsat 數據在近紅外波段的地表反射率.

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

表1 NDWI 閾值

2.2 水位與水體面積關系建立方法

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

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

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

式中,lj為水體影像對應的水位值,aj為該水位值日期的水體面積,a為水體總面積,L為水體影像的加權平均水位值.

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

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

2.3 精度評價

對于水位-水體面積函數關系,使用確定性系數R2(見式(3))進行擬合方程的精度評定.

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

3 結果與討論

3.1 東洞庭湖水域水體提取結果

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

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

表2 水體提取結果與精度評價

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

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

圖2 東洞庭湖水域水體提取結果

3.2 東洞庭湖水位與水體面積關系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結語

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

(1)Google Earth Engine 平臺可以高效下載并處理大量的遙感數據,具有很強的計算能力.

(2)根據東洞庭湖的水文特征,利用Otsu 方法確定利用12 個月閾值分別為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值進行水體提取和面積計算.東洞庭湖水域面積隨季節(jié)變化大,冬夏差異明顯,7～9月水域面積最大(平均面積為1078.1585 km2),1～3月水域面積最小(平均面積為332.6463 km2).研究時期內最大面積為1204.752 km2,最小面積為240.722 km2.

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