基于最大值合成和最大類間方差法的東洞庭湖灘涂提取研究

鄧家航 ，姜鎵偉 ，莫登奎，嚴恩萍

(中南林業(yè)科技大學林學院，湖南 長沙 410004)

湖泊與河流在自然條件的控制下，年際內(nèi)存在汛期和旱季變化，水陸交替演變，灘涂作為陸地和水體之間的過渡地帶，主要由淤泥質(zhì)或砂質(zhì)的無植被區(qū)域組成[1]。灘涂作為濕地的重要組成部分，提供著多種重要的生態(tài)服務功能，如促進碳匯、保護生物多樣性、維持區(qū)域生態(tài)平衡等[2]。此外，從古至今人類就通過圍墾灘涂發(fā)展農(nóng)業(yè)和水產(chǎn)養(yǎng)殖業(yè)取得了良好的經(jīng)濟效益[3]。但過度的開墾灘涂，容易導致湖泊生態(tài)環(huán)境退化，灘涂大面積露出，會引發(fā)一系列生態(tài)問題[4]。因此，在湖泊生態(tài)面臨嚴峻威脅的當下，及時、準確地掌握灘涂資源現(xiàn)狀，是實現(xiàn)灘涂可持續(xù)管理的前提[5]。

灘涂與水體密切相關(guān)，受水體周期性漲落的影響，灘涂在水位漲高時被淹沒，在水位下降時暴露在水面之外，與水體處于動態(tài)變化之中。受水位變化的影響，衛(wèi)星過境時灘涂裸露面積也有著極大的隨機性，因此僅靠單一時相的遙感影像不能確保灘涂信息提取工作的順利進行[6]。針對水位實時變化帶來的限制，國內(nèi)外學者進行了相關(guān)研究。如Murray等[7]使用區(qū)域潮汐模型，確定極端潮位圖像的瞬時水邊線，以低潮瞬時水邊線減去高潮瞬時水邊線，最終完成了東亞濱海灘涂的繪制；Wang等[8]使用GEE和決策樹算法，分析了所有Landsat時序影像，確定了最高潮位線和低水線的水體頻率閾值范圍，然后識別并劃分潮間帶，從而生成1986—2016年30 m分辨率的中國東部沿海灘涂地圖；李振等[9]通過建立決策樹法對兩期Landsat影像進行分類，對這兩期土地利用類型進行疊加分析，得到了膠州灣海岸帶土地利用變化矩陣。然而，目前國內(nèi)外學者多聚焦于濱海灘涂提取，尚缺乏對內(nèi)陸灘涂的相關(guān)研究。

隨著遙感技術(shù)的快速發(fā)展，各類衛(wèi)星傳感器在空間分辨率、時間分辨率和光譜分辨率等方面不斷取得突破，為灘涂資源的動態(tài)監(jiān)測和可持續(xù)管理提供了更多機遇[10-13]。在灘涂資源遙感監(jiān)測工作中，空間分辨率和時間分辨率是選擇遙感數(shù)據(jù)的重要考量[14]。在眾多遙感數(shù)據(jù)中，灘涂制圖研究常采用Sentinel系列和Landsat系列。而GEE(Google Earth Engine)云平臺存儲了豐富的歷史遙感影像數(shù)據(jù)(如Landsat、Sentinel等)，支持云端并行計算，讓GEE用戶可以免費且方便地調(diào)用和分析海量的遙感數(shù)據(jù)資源，不用耗費大量時間下載數(shù)據(jù)，并且不依賴硬件條件，是大規(guī)模的長時間序列遙感影像處理的絕佳技術(shù)平臺[15-16]。例如，張康永[17]基于GEE平臺，結(jié)合隨機森林分類法和面向?qū)ο蠹夹g(shù)對時間序列影像數(shù)據(jù)進行制圖，最終得到了杭州灣以北濱海灘涂30年來的變化情況；Jia等[18]開發(fā)了一種結(jié)合最大光譜指數(shù)合成法(maximum spectral index composite，MSIC)和最大類間算法(Otsu algorithm，OTSU)名為MSIC-OA的方法，以Sentinel-2為數(shù)據(jù)源繪制了10 m空間分辨率的中國潮灘圖(China_Tidal Flat，CTF)；程麗娜等[19]構(gòu)建了高質(zhì)量密集時序Sentinel-2影像堆棧，通過MSIC-OA建立分類模型，完成了福建漳江口紅樹林自然保護區(qū)的潮間帶濕地分類。

目前，基于GEE和Sentinel-2，通過MSIC-OA對濱海灘涂和濕地進行提取都取得了精準的分類效果，但尚未有研究將此方法用于提取內(nèi)陸灘涂。且此前的研究都聚焦于灘涂年際最大面積提取，而灘涂完全暴露于水面之外的時間是極為短暫的。因此，有必要量化年際內(nèi)灘涂在不同水位條件下的分布情況，分析灘涂裸露的時間，實現(xiàn)根據(jù)水位高度即能預測灘涂資源情況。綜上，本研究旨在實現(xiàn)基于Sentinel-2和GEE云平臺，通過MSIC-OA算法自動、高精度地提取內(nèi)陸灘涂。

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

1.1 研究區(qū)概況

東洞庭湖位于長江中游荊江河段南側(cè)，地理坐標范圍為112°19′—113°05′E，28°56′—29°35′N。研究區(qū)的總體面積約為14萬hm2，地勢平坦，平均海拔50 m，洞庭湖灘涂主要分布在此處。該區(qū)域地處亞熱帶濕潤氣候區(qū)，日照充足，雨量充沛，年均氣溫17 ℃，年均降水量1 200～1 300 mm[20](見圖1)。

1.2 數(shù)據(jù)來源

1.2.1 Sentinel-2 數(shù)據(jù) 本研究采用的遙感數(shù)據(jù)源為Sentinel-2。Sentinel-2是高分辨率多光譜成像衛(wèi)星，包含了2顆同時運作的相同衛(wèi)星，分別是Sentinel-2A和Sentinel-2B，用于對土壤、植被和水覆蓋的監(jiān)測，以及內(nèi)陸水道和沿海地區(qū)的觀測等，還可用于緊急救援服務。單顆衛(wèi)星的重訪周期為10d，兩顆互補。同時進入運行狀態(tài)后，雙星每5 d更新1次成像數(shù)據(jù)。Sentinel-2衛(wèi)星高度為786 km，可覆蓋13個光譜波段，影像幅寬290 km，最高空間分辨率可達10 m[21]。Sentinel-2各波段信息見表1。

表1 Sentinel-2波段介紹Tab. 1 Sentinel-2 band information波段中心波長/nm空間分辨率/m描述B144360氣溶膠B249010藍B356010綠B466510紅B570520紅邊1B674020紅邊2B778320紅邊3B884210近紅外B8A86520紅邊4B994060水汽波段B101 37560卷云B111 61020短波紅外1(SWIR1)B122 19020短波紅外2(SWIR2)

1.2.2 目視解譯數(shù)據(jù) 由于水位變化的不確定性，通常在水位最低時刻缺乏衛(wèi)星影像，因此本研究選取無云影像中水位最低的影像作為替代，分別選取日期為2018-02-12和2019-12-04的影像數(shù)據(jù)作為2018年和2019年最低水位影像，選取2021-01-02的影像數(shù)據(jù)作為2020年和2021年最低水位影像。同時，依據(jù)水位高差選取了日期為2019-11-19、2019-11-14、2019-11-09、2020-03-18和2020-11-08的影像數(shù)據(jù)。對所選日期的Sentinel-2影像進行目視解譯，區(qū)分為灘涂、水體和植被三種地類。所選日期影像水位高度見表2。

表2 影像日期和水位影像高度Tab.2 Dates and water level of Images影像日期/(年-月-日)水位/m2018-02-1221.542019-12-0420.692021-01-0221.322019-11-1921.672019-11-1422.782019-11-0923.432020-03-1824.702020-11-0825.69

2 研究方法

2.1 總體思路

圖2是本文提取東洞庭湖灘涂的方法流程，主要步驟如下：(1)基于GEE云平臺，構(gòu)建Sentinel-2高質(zhì)量密集時序影像堆棧；(2)選取合適的光譜指數(shù)，基于MSIC算法獲取研究區(qū)最大水面和最小水面影像；(3)通過Otsu算法進行圖像二值化分類。

圖2 灘涂提取流程圖Fig.2 Workflow of inland beaches extracion

最大光譜指數(shù)合成法(MSIC)是由GEE API提供的一種算法,即“imageCollection.qualityMosaic()”。該算法根據(jù)時間序列影像集合中具有所選光譜指數(shù)最大值的像素來設置合成影像中的每個像素，即在合成影像中，每個像素都是從時間序列影像集合的不同影像中選擇的，并且每個像素都代表其自身位置的最大光譜指數(shù)值。Otsu算法即最大類間方差法，能根據(jù)圖像的直方圖自動選擇最佳閾值，將圖像分成背景和目標兩類，自動選擇使類間方差最大、類內(nèi)方差最小的最優(yōu)分割閾值，自動實現(xiàn)圖像的二值分割。

2.2 年際最大面積灘涂提取

本文將灘涂定義為水位漲高時被水淹沒，水位降低時露出水面的砂質(zhì)、淤泥質(zhì)的無植被區(qū)域。在年際最大面積灘涂提取工作中，首要任務是提取最大水體邊界。為此，需要一個對水體敏感的光譜指數(shù)用于合成最高水位影像，歸一化差分水體指數(shù)NDWI(normalized difference water index)可以很好地突出水體與非水體(灘涂和植被)的特征差異。此外，灘涂和植被的歸一化差分植被指數(shù)NDVI(normalized difference vegetable index)值高于水體像素的NDVI值，可以突出植被與非植被(水體和灘涂)的特征差異，因此選擇NDVI作為指定光譜指數(shù)合成最低水位影像。

首先，利用Otsu算法計算最高水位影像最佳分割閾值并進行分割。由于使用的光譜指數(shù)為NDWI，因此分類結(jié)果為水體和非水體，在水體分類中保留最大面積斑塊，即最大水體邊界。其次，對最低水位影像進行掩膜提取并再次使用Otsu算法進行二值分類，得到植被和非植被(灘涂和水體)。然后，去除植被，保留灘涂和水體部分，利用該區(qū)域?qū)ψ罡咚挥跋襁M行掩膜提取并使用Otsu算法進行二值分類，得到灘涂和水體的二值圖，去除水體，最終得到灘涂。

2.3 灘涂面積與水位高度定量分析

由于灘涂與水體呈此消彼長的關(guān)系，不同高度水位保持時間也不同，因此有必要定量分析灘涂面積與水位高度的關(guān)系。針對東洞庭湖灘涂裸露的情況，本研究選取不同水位高度的影像，分別進行灘涂提取。使用GEE獲得每景影像的NDWI和NDVI影像，利用最大水體邊界對NDVI影像進行掩膜提取并使用Otsu算法進行二值分類，得到植被和非植被(灘涂和水體)。去除植被，利用該區(qū)域?qū)ψ罡咚挥跋襁M行掩膜提取，再次使用Otsu算法進行二值分類，得到灘涂和水體的二值圖，去除水體，最終得到不同水位高度下的灘涂面積。

2.4 精度驗證

本研究采用建立混淆矩陣來進行分類結(jié)果精度評價?；煜仃囀蔷仍u價的一種標準格式，根據(jù)像元的地表真實情況與該像元分類情況進行比較計算得出，主要指標有總分類精度、Kappa系數(shù)、制圖精度和用戶精度等?？偡诸惥戎副徽_分類的樣本數(shù)占驗證樣本數(shù)的百分比；Kappa系數(shù)用于一致性檢驗，是一種衡量分類精度的指標；用戶精度指該類別被正確分類的占比。

3 結(jié)果與分析

3.1 東洞庭湖年際灘涂最大面積分類結(jié)果及精度分析

東洞庭湖2018—2021年年際灘涂最大面積分類結(jié)果如圖3所示，灘涂面積分別是20608.30hm2、22956.56hm2、22715.10hm2和22867.30hm2。

圖3 2018—2021年東洞庭湖灘涂分類結(jié)果Fig.3 Classification results of East Dongting Lake from 2018 to 2021

由于灘涂的特殊性質(zhì)，在提取工作中需確保灘涂在指定時間內(nèi)曾被水淹沒，因此本研究使用經(jīng)“最大面積水體斑塊提取”的水體作為灘涂潛在的最大面積區(qū)域，該水體斑塊并不是完整的東洞庭湖水面，但可以用來驗證本研究精度。在得到2018年、2019年、2020年和2021年東洞庭湖灘涂分類后處理的結(jié)果后，采用制圖精度、總分類精度、用戶精度和Kappa系數(shù)等分類精度參數(shù)，利用隨機生成的驗證點對4期分類結(jié)果進行精度分析。

表3為根據(jù)隨機驗證點生成的東洞庭湖2018年地物分類的混淆矩陣，結(jié)果表明：2018年灘涂及水體、植被分類結(jié)果總分類精度和Kappa系數(shù)分別為96.3%和0.95。另經(jīng)計算，2019年、2020年和2021年分類結(jié)果的總分類精度分別為96.0%、95.7%和95.3%，Kappa系數(shù)分別為0.94、0.94和0.93。主要存在的誤分情況在于灘涂與植被之間和灘涂與水體之間，灘涂的NDWI值和NDVI值都介于植被與水體之間，灘涂的高含水量和植被區(qū)域的低覆蓋度都會使得灘涂與水體和植被的光譜特征相似，這是造成錯分的主要原因，而植被與水體則極少存在誤分。

3.2 東洞庭湖灘涂面積與水位高度定量分析及精度分析

東洞庭湖在水位高度分別為20.69 m、21.67 m、22.78 m、23.43 m、24.70 m和25.69 m時，灘涂裸露的面積分別為22589.87hm2、16533.60hm2、14990.73hm2、13072.23hm2、7030.25hm2和0 hm2(見圖4)。此處需要聲明一點：在進行灘涂面積和水位高度定量分析時，NDWI最大值合成影像的時間跨度為2018—2021年，由于不同年份的最高水位不同，湖水的最大淹沒區(qū)域也會有差異，即存在當年未被湖水淹沒的區(qū)域在其他年份被湖水淹沒，因此“最大面積水體斑塊提取”的水體面積相較于前一年合成更大，灘涂面積也會隨之變大。

(a)2019-12-04 / 20.69 m(b)2019-11-19 / 21.67 m(c)2019-11-14 / 22.78 m

表4是根據(jù)隨機驗證點生成的東洞庭湖2019-12-04(水位：20.69 m)地物分類的混淆矩陣，結(jié)果表明：2019-12-04灘涂及水體、植被分類結(jié)果總分類精度和Kappa系數(shù)分別為97.3%和0.96。2019-11-19(水位：21.67 m)、2019-11-14(水位：22.78 m)、2019-11-09(水位：23.43 m)、2020-03-18(水位：24.70 m)和2020-11-08(水位：25.69 m)的總分類精度分別為93.0%、97.7%、97.0%、96.3%和100%，Kappa系數(shù)分別為0.93、0.95、0.96、0.95和1。

表4 2019-12-04東洞庭湖分類結(jié)果精度評價Tab.4 Classification confusion matrix and precision a-nalysis of inland beaches on 2019-12-04實際地物類型植被灘涂水體總計用戶精度/% 植被1000 0100100分類結(jié)果 灘涂297 110097 水體14 9510095總計103101 96300制圖精度/%97.196.399.0總分類精度/%97.3Kappa系數(shù)0.96

3.3 東洞庭湖灘涂面積預測

2018—2021年灘涂出現(xiàn)時間及面積如表5所示。年際內(nèi)灘涂完全裸露的時間是極短暫的，完全裸露時間以2019年最長，為33 d；2020年最短，僅9 d。而灘涂被完全淹沒的時間以2021年最長，為184 d；其后分別是2020年、2019年和2018年，灘涂被完全淹沒的時間分別為165 d、131 d和127d。在2018—2021年中，東洞庭湖灘涂在水位為23.43～24.70 m的維持時間最長，為218 d，此時灘涂面積為7 030.25～13 072.23 hm2；其次是21.67～22.78 m，維持時間為183 d，此時灘涂面積為14 990.73～16 533.60 hm2；水位為20.69～21.67 m時維持170 d，此時灘涂面積為16 533.60～22 589.87 hm2。

表5 2018—2021年東洞庭湖灘涂出現(xiàn)時間及面積Tab.5 Exposure time and area of East Dongting Lake inland beaches in 2018-2021水位高度/m灘涂面積/hm2灘涂出現(xiàn)天數(shù)/d2018年2019年2020年2021年總計<20.69>22 589.8720339238520.69～21.6716 533.60～22 589.874430395717021.67～22.7814 990.73～16 533.607224493818322.78～23.4313 072.23～14 990.73182731189423.43～24.707 030.25～13 072.234092563021824.70～25.69<7 030.2544281615103>25.690127131165184607

4 結(jié)論與討論

4.1 結(jié)論

本研究基于GEE云平臺，以Sentinel-2為數(shù)據(jù)源，結(jié)合MSIC和Otsu算法(MSIC-Otsu)對東洞庭湖灘涂進行提取，完成了內(nèi)陸灘涂的快速、自動和高精度提取，并對東洞庭湖灘涂和水位高度進行了定量分析，實現(xiàn)了根據(jù)水位即可預測灘涂面積的功能。研究結(jié)論如下：

(1)MSIC-Otsu算法能夠高度完成內(nèi)陸灘涂提取任務。2018—2021年，東洞庭湖灘涂年際最大裸露面積為20608.30～22956.56hm2，總分類精度達96.3%，Kappa系數(shù)達0.95；在灘涂面積與水位高度定量分析中，總分類精度達97.3%，Kappa系數(shù)達0.96。

(2)Sentinel-2較短的重訪周期，極大地增加了獲取最高和最低水位時期影像數(shù)據(jù)的機會，有效解決了水體動態(tài)變化造成的困難，為內(nèi)陸灘涂資源監(jiān)測提供了更多的機遇。而GEE云平臺海量的數(shù)據(jù)資源、支持云端并行計算的能力以及支持用戶快速訪問、調(diào)用和分析Sentinel-2數(shù)據(jù)產(chǎn)品的服務同樣為監(jiān)測內(nèi)陸灘涂資源提供了巨大助力。

(3)東洞庭湖灘涂在水位高度達25.69 m時被完全淹沒，年際內(nèi)完全淹沒時間約占35%～50%。水位低于25.69 m時，東洞庭湖灘涂會部分或者全部裸露，對灘涂面積和水位高度進行定量分析，實現(xiàn)了根據(jù)水位即可預測灘涂面積的功能，并預估其維持時間。

4.2 討論

首先，本研究使用NDWI作為指定波段用于提取最大水面，但NDWI并不能完全抑制含水量過高的灘涂，因此導致部分水體存在錯誤分類；其次，由于GEE云平臺僅提供2015年3月之后的Sentinel-2數(shù)據(jù)，因此無法獲得此前覆蓋研究區(qū)的歷史觀測數(shù)據(jù)，為長時間序列的研究造成數(shù)據(jù)源限制。此外，盡管本研究通過構(gòu)建密集時序影像堆棧，省去了人工篩選極端水位影像的步驟，但無法保證衛(wèi)星會在極端水位時拍攝影像，這也是本研究能否精確提取灘涂的重要影響因素之一。最后，受制于Sentinel-2數(shù)據(jù)自身的空間分辨率，研究中出現(xiàn)低于10 m的斑塊時，也會對灘涂提取工作的精度造成一定的影響。