基于最大值合成和最大类间方差法的东洞庭湖滩涂提取研究

邓家航 ，姜镓伟 ，莫登奎，严恩萍

(中南林业科技大学林学院，湖南 长沙 410004)

湖泊与河流在自然条件的控制下，年际内存在汛期和旱季变化，水陆交替演变，滩涂作为陆地和水体之间的过渡地带，主要由淤泥质或砂质的无植被区域组成[1]。滩涂作为湿地的重要组成部分，提供着多种重要的生态服务功能，如促进碳汇、保护生物多样性、维持区域生态平衡等[2]。此外，从古至今人类就通过围垦滩涂发展农业和水产养殖业取得了良好的经济效益[3]。但过度的开垦滩涂，容易导致湖泊生态环境退化，滩涂大面积露出，会引发一系列生态问题[4]。因此，在湖泊生态面临严峻威胁的当下，及时、准确地掌握滩涂资源现状，是实现滩涂可持续管理的前提[5]。

滩涂与水体密切相关，受水体周期性涨落的影响，滩涂在水位涨高时被淹没，在水位下降时暴露在水面之外，与水体处于动态变化之中。受水位变化的影响，卫星过境时滩涂裸露面积也有着极大的随机性，因此仅靠单一时相的遥感影像不能确保滩涂信息提取工作的顺利进行[6]。针对水位实时变化带来的限制，国内外学者进行了相关研究。如Murray等[7]使用区域潮汐模型，确定极端潮位图像的瞬时水边线，以低潮瞬时水边线减去高潮瞬时水边线，最终完成了东亚滨海滩涂的绘制；Wang等[8]使用GEE和决策树算法，分析了所有Landsat时序影像，确定了最高潮位线和低水线的水体频率阈值范围，然后识别并划分潮间带，从而生成1986—2016年30 m分辨率的中国东部沿海滩涂地图；李振等[9]通过建立决策树法对两期Landsat影像进行分类，对这两期土地利用类型进行叠加分析，得到了胶州湾海岸带土地利用变化矩阵。然而，目前国内外学者多聚焦于滨海滩涂提取，尚缺乏对内陆滩涂的相关研究。

随着遥感技术的快速发展，各类卫星传感器在空间分辨率、时间分辨率和光谱分辨率等方面不断取得突破，为滩涂资源的动态监测和可持续管理提供了更多机遇[10-13]。在滩涂资源遥感监测工作中，空间分辨率和时间分辨率是选择遥感数据的重要考量[14]。在众多遥感数据中，滩涂制图研究常采用Sentinel系列和Landsat系列。而GEE(Google Earth Engine)云平台存储了丰富的历史遥感影像数据(如Landsat、Sentinel等)，支持云端并行计算，让GEE用户可以免费且方便地调用和分析海量的遥感数据资源，不用耗费大量时间下载数据，并且不依赖硬件条件，是大规模的长时间序列遥感影像处理的绝佳技术平台[15-16]。例如，张康永[17]基于GEE平台，结合随机森林分类法和面向对象技术对时间序列影像数据进行制图，最终得到了杭州湾以北滨海滩涂30年来的变化情况；Jia等[18]开发了一种结合最大光谱指数合成法(maximum spectral index composite，MSIC)和最大类间算法(Otsu algorithm，OTSU)名为MSIC-OA的方法，以Sentinel-2为数据源绘制了10 m空间分辨率的中国潮滩图(China_Tidal Flat，CTF)；程丽娜等[19]构建了高质量密集时序Sentinel-2影像堆栈，通过MSIC-OA建立分类模型，完成了福建漳江口红树林自然保护区的潮间带湿地分类。

目前，基于GEE和Sentinel-2，通过MSIC-OA对滨海滩涂和湿地进行提取都取得了精准的分类效果，但尚未有研究将此方法用于提取内陆滩涂。且此前的研究都聚焦于滩涂年际最大面积提取，而滩涂完全暴露于水面之外的时间是极为短暂的。因此，有必要量化年际内滩涂在不同水位条件下的分布情况，分析滩涂裸露的时间，实现根据水位高度即能预测滩涂资源情况。综上，本研究旨在实现基于Sentinel-2和GEE云平台，通过MSIC-OA算法自动、高精度地提取内陆滩涂。

1 研究区概况及数据来源

1.1 研究区概况

东洞庭湖位于长江中游荆江河段南侧，地理坐标范围为112°19′—113°05′E，28°56′—29°35′N。研究区的总体面积约为14万hm2，地势平坦，平均海拔50 m，洞庭湖滩涂主要分布在此处。该区域地处亚热带湿润气候区，日照充足，雨量充沛，年均气温17 ℃，年均降水量1 200～1 300 mm[20](见图1)。

1.2 数据来源

1.2.1 Sentinel-2 数据 本研究采用的遥感数据源为Sentinel-2。Sentinel-2是高分辨率多光谱成像卫星，包含了2颗同时运作的相同卫星，分别是Sentinel-2A和Sentinel-2B，用于对土壤、植被和水覆盖的监测，以及内陆水道和沿海地区的观测等，还可用于紧急救援服务。单颗卫星的重访周期为10d，两颗互补。同时进入运行状态后，双星每5 d更新1次成像数据。Sentinel-2卫星高度为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 目视解译数据 由于水位变化的不确定性，通常在水位最低时刻缺乏卫星影像，因此本研究选取无云影像中水位最低的影像作为替代，分别选取日期为2018-02-12和2019-12-04的影像数据作为2018年和2019年最低水位影像，选取2021-01-02的影像数据作为2020年和2021年最低水位影像。同时，依据水位高差选取了日期为2019-11-19、2019-11-14、2019-11-09、2020-03-18和2020-11-08的影像数据。对所选日期的Sentinel-2影像进行目视解译，区分为滩涂、水体和植被三种地类。所选日期影像水位高度见表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云平台，构建Sentinel-2高质量密集时序影像堆栈；(2)选取合适的光谱指数，基于MSIC算法获取研究区最大水面和最小水面影像；(3)通过Otsu算法进行图像二值化分类。

图2 滩涂提取流程图Fig.2 Workflow of inland beaches extracion

最大光谱指数合成法(MSIC)是由GEE API提供的一种算法,即“imageCollection.qualityMosaic()”。该算法根据时间序列影像集合中具有所选光谱指数最大值的像素来设置合成影像中的每个像素，即在合成影像中，每个像素都是从时间序列影像集合的不同影像中选择的，并且每个像素都代表其自身位置的最大光谱指数值。Otsu算法即最大类间方差法，能根据图像的直方图自动选择最佳阈值，将图像分成背景和目标两类，自动选择使类间方差最大、类内方差最小的最优分割阈值，自动实现图像的二值分割。

2.2 年际最大面积滩涂提取

本文将滩涂定义为水位涨高时被水淹没，水位降低时露出水面的砂质、淤泥质的无植被区域。在年际最大面积滩涂提取工作中，首要任务是提取最大水体边界。为此，需要一个对水体敏感的光谱指数用于合成最高水位影像，归一化差分水体指数NDWI(normalized difference water index)可以很好地突出水体与非水体(滩涂和植被)的特征差异。此外，滩涂和植被的归一化差分植被指数NDVI(normalized difference vegetable index)值高于水体像素的NDVI值，可以突出植被与非植被(水体和滩涂)的特征差异，因此选择NDVI作为指定光谱指数合成最低水位影像。

首先，利用Otsu算法计算最高水位影像最佳分割阈值并进行分割。由于使用的光谱指数为NDWI，因此分类结果为水体和非水体，在水体分类中保留最大面积斑块，即最大水体边界。其次，对最低水位影像进行掩膜提取并再次使用Otsu算法进行二值分类，得到植被和非植被(滩涂和水体)。然后，去除植被，保留滩涂和水体部分，利用该区域对最高水位影像进行掩膜提取并使用Otsu算法进行二值分类，得到滩涂和水体的二值图，去除水体，最终得到滩涂。

2.3 滩涂面积与水位高度定量分析

由于滩涂与水体呈此消彼长的关系，不同高度水位保持时间也不同，因此有必要定量分析滩涂面积与水位高度的关系。针对东洞庭湖滩涂裸露的情况，本研究选取不同水位高度的影像，分别进行滩涂提取。使用GEE获得每景影像的NDWI和NDVI影像，利用最大水体边界对NDVI影像进行掩膜提取并使用Otsu算法进行二值分类，得到植被和非植被(滩涂和水体)。去除植被，利用该区域对最高水位影像进行掩膜提取，再次使用Otsu算法进行二值分类，得到滩涂和水体的二值图，去除水体，最终得到不同水位高度下的滩涂面积。

2.4 精度验证

本研究采用建立混淆矩阵来进行分类结果精度评价。混淆矩阵是精度评价的一种标准格式，根据像元的地表真实情况与该像元分类情况进行比较计算得出，主要指标有总分类精度、Kappa系数、制图精度和用户精度等。总分类精度指被正确分类的样本数占验证样本数的百分比；Kappa系数用于一致性检验，是一种衡量分类精度的指标；用户精度指该类别被正确分类的占比。

3 结果与分析

3.1 东洞庭湖年际滩涂最大面积分类结果及精度分析

东洞庭湖2018—2021年年际滩涂最大面积分类结果如图3所示，滩涂面积分别是20608.30hm2、22956.56hm2、22715.10hm2和22867.30hm2。

图3 2018—2021年东洞庭湖滩涂分类结果Fig.3 Classification results of East Dongting Lake from 2018 to 2021

由于滩涂的特殊性质，在提取工作中需确保滩涂在指定时间内曾被水淹没，因此本研究使用经“最大面积水体斑块提取”的水体作为滩涂潜在的最大面积区域，该水体斑块并不是完整的东洞庭湖水面，但可以用来验证本研究精度。在得到2018年、2019年、2020年和2021年东洞庭湖滩涂分类后处理的结果后，采用制图精度、总分类精度、用户精度和Kappa系数等分类精度参数，利用随机生成的验证点对4期分类结果进行精度分析。

表3为根据随机验证点生成的东洞庭湖2018年地物分类的混淆矩阵，结果表明：2018年滩涂及水体、植被分类结果总分类精度和Kappa系数分别为96.3%和0.95。另经计算，2019年、2020年和2021年分类结果的总分类精度分别为96.0%、95.7%和95.3%，Kappa系数分别为0.94、0.94和0.93。主要存在的误分情况在于滩涂与植被之间和滩涂与水体之间，滩涂的NDWI值和NDVI值都介于植被与水体之间，滩涂的高含水量和植被区域的低覆盖度都会使得滩涂与水体和植被的光谱特征相似，这是造成错分的主要原因，而植被与水体则极少存在误分。

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年，由于不同年份的最高水位不同，湖水的最大淹没区域也会有差异，即存在当年未被湖水淹没的区域在其他年份被湖水淹没，因此“最大面积水体斑块提取”的水体面积相较于前一年合成更大，滩涂面积也会随之变大。

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

表4是根据随机验证点生成的东洞庭湖2019-12-04(水位：20.69 m)地物分类的混淆矩阵，结果表明：2019-12-04滩涂及水体、植被分类结果总分类精度和Kappa系数分别为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系数分别为0.93、0.95、0.96、0.95和1。

表4 2019-12-04东洞庭湖分类结果精度评价Tab.4 Classification confusion matrix and precision a-nalysis of inland beaches on 2019-12-04实际地物类型植被滩涂水体总计用户精度/% 植被1000 0100100分类结果 滩涂297 110097 水体14 9510095总计103101 96300制图精度/%97.196.399.0总分类精度/%97.3Kappa系数0.96

3.3 东洞庭湖滩涂面积预测

2018—2021年滩涂出现时间及面积如表5所示。年际内滩涂完全裸露的时间是极短暂的，完全裸露时间以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年东洞庭湖滩涂出现时间及面积Tab.5 Exposure time and area of East Dongting Lake inland beaches in 2018-2021水位高度/m滩涂面积/hm2滩涂出现天数/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 结论与讨论

4.1 结论

本研究基于GEE云平台，以Sentinel-2为数据源，结合MSIC和Otsu算法(MSIC-Otsu)对东洞庭湖滩涂进行提取，完成了内陆滩涂的快速、自动和高精度提取，并对东洞庭湖滩涂和水位高度进行了定量分析，实现了根据水位即可预测滩涂面积的功能。研究结论如下：

(1)MSIC-Otsu算法能够高度完成内陆滩涂提取任务。2018—2021年，东洞庭湖滩涂年际最大裸露面积为20608.30～22956.56hm2，总分类精度达96.3%，Kappa系数达0.95；在滩涂面积与水位高度定量分析中，总分类精度达97.3%，Kappa系数达0.96。

(2)Sentinel-2较短的重访周期，极大地增加了获取最高和最低水位时期影像数据的机会，有效解决了水体动态变化造成的困难，为内陆滩涂资源监测提供了更多的机遇。而GEE云平台海量的数据资源、支持云端并行计算的能力以及支持用户快速访问、调用和分析Sentinel-2数据产品的服务同样为监测内陆滩涂资源提供了巨大助力。

(3)东洞庭湖滩涂在水位高度达25.69 m时被完全淹没，年际内完全淹没时间约占35%～50%。水位低于25.69 m时，东洞庭湖滩涂会部分或者全部裸露，对滩涂面积和水位高度进行定量分析，实现了根据水位即可预测滩涂面积的功能，并预估其维持时间。

4.2 讨论

首先，本研究使用NDWI作为指定波段用于提取最大水面，但NDWI并不能完全抑制含水量过高的滩涂，因此导致部分水体存在错误分类；其次，由于GEE云平台仅提供2015年3月之后的Sentinel-2数据，因此无法获得此前覆盖研究区的历史观测数据，为长时间序列的研究造成数据源限制。此外，尽管本研究通过构建密集时序影像堆栈，省去了人工筛选极端水位影像的步骤，但无法保证卫星会在极端水位时拍摄影像，这也是本研究能否精确提取滩涂的重要影响因素之一。最后，受制于Sentinel-2数据自身的空间分辨率，研究中出现低于10 m的斑块时，也会对滩涂提取工作的精度造成一定的影响。