Sentinel-2密集时间序列数据和Google Earth Engine的潮间带湿地快速自动分类

程丽娜，钟才荣，李晓燕，贾明明，王宗明，4，毛德华

1.吉林大学地球科学学院,长春130061;

2.海南省林业科学研究院（海南省红树林研究院）,海口571100;

3.中国科学院东北地理与农业生态研究所湿地生态与环境重点实验室,长春130102;

4.国家地球系统科学数据中心,北京100101

1 引言

潮间带是平均最高潮位和最低潮位间的海陆过渡地带，是海平面上升等全球气候变化的敏感区和影响显著区（鲍毅新等，2007；胡嘉东等，2009；史本伟，2012）。潮间带湿地作为滨海湿地的重要组成部分，在维持海岸线稳定（万峻等，2010）、提供生物栖息地（胡嘉东等，2009）、加速物质交换（孟伟等，2009）和促进碳汇（仲启铖等，2015）等方面都发挥着重要的作用。此外，通过围垦潮间带湿地发展农业和水产养殖业可创造极高的经济效益（张文开，2001）。因此，受海岸侵蚀、海平面上升及人类活动的影响，潮间带湿地生境面临着严峻的威胁。及时准确地掌握潮间带湿地的现状和演变规律，对实现潮间带湿地的可持续管理目标具有重要意义。

遥感技术的推陈出新为潮间带湿地信息的分类提供了更多的机遇。局地尺度高精度潮间带湿地制图通常采用高空间分辨率遥感影像为数据源，如QuickBird、IKONOS、WorldView 以及无人机数据等（Neukermans 等，2008；Wang 等，2004；徐逸等，2021）。然而，由于上述数据的价格和可获取性等原因，在大范围潮间带制图的研究中常采用易获取的Landsat系列或Sentinel系列卫星遥感影像（宫鹏等，2016；宫鹏，2021；姚云长等，2016）。

周期性的潮水淹没是潮间带湿地制图的最大挑战，快速、准确地获得最高潮和最低潮时期的遥感影像是准确提取潮间带湿地信息的关键。近年来，国内外学者在潮间带湿地制图方面开展了一系列研究。如，Murray 等（2012）依据区域潮汐模型，分别获取高潮和低潮影像以确定瞬时水边线，通过将二者叠加，进行东亚滨海地区潮间带滩涂的绘制。胥为（2017）融合Landsat 8 和Sentinel-1 数据，构建基于专家知识的决策树分类方法，对研究区潮间带盐沼湿地进行分类制图。Wang 等（2020）利用时间序列Landsat 影像，在GEE 平台提出了一种基于像素和物候的算法，并制作了30 m 空间分辨率的中国滨海湿地地图。Murray 等（2019）利用时间序列的Landsat 数据，在GEE 云平台上开发了全球范围内包括滩涂、永久水体和其他3种土地覆盖类型的训练数据集，应用随机森林算法获取了具有滩涂时间序列的全球分布图。董迪等（2020）结合分层分类与随机森林算法，采用Sentinel-1/2 卫星数据，实现了红树林与互花米草的自动提取。

GEE（Google Earth Engine）云平台将Landsat/Sentinel 等可以公开获取的遥感图像数据存储在谷歌的磁盘阵列中，使得GEE 用户可以方便的提取、调用和分析海量的遥感大数据资源，为大尺度长时间序列遥感分析提供了巨大的潜力（Dong 等，2016）。例如，吴亚茜等（2018）基于GEE 云平台提取MOD09A1 与Landsat 各数据源的时间序列数据，对盐城潮间带湿地盐沼植被的多年遥感物候特征及变化进行了分析；Li 等（2019）通过构建Sentinel-2 时间序列的归一化植被指数表征红树林物种物候轨迹特征，依托GEE 云平台验证随机森林算法对红树林物种分类的可行性；Tian等（2020）在GEE 云平台开发了一种基于像素的物候特征合成方法（Ppf-CM），对比深度学习方法和支持向量机的性能，对1995年—2017年中国北部湾互花米草的空间分布变化进行了研究。Xu 等（2021）通过GEE 平台，提出了基于Sentinel-2 时间序列影像的像素合成物候植被指数法（PVI），并结合Otsu阈值，实现互花米草与其他土地覆盖类型的区分。Jia 等（2021）基于GEE 平台，选用Sentinel-2 数据，结合最大光谱指数合成算法（MSIC）和Otsu算法（OA），绘制了10 m 空间分辨率的中国滩涂分布图。然而，目前尚缺乏基于Sentinel-2 和GEE云平台的潮间带湿地自动分类研究。

综上，本研究的目的是实现基于Sentinel-2 和GEE 云平台的潮间带湿地自动、高精度提取。为实现这一目标，本文将从以下3点展开：（1）构建高质量密集时序Sentinel-2 影像堆栈；（2）基于GEE 云平台开发一种自动、快速、高精度的潮间带湿地信息分类方案；（3）以福建漳江口潮间带为研究区，验证该方案的可行性。

2 研究区及数据

2.1 研究区概况

本研究的研究区是漳江口红树林国家级自然保护区，位于福建省漳州市云霄县（图1），23°54′17″N—23°57′29″N，117°23′53″E—117°27′33″E，海拔-6—8 m。研究区为亚热带海洋季风气候，光、热、水资源丰富。受不正规半日潮影响，平均潮差为2.33 m，最大潮差为4.66 m。保护区于1992年建立，1997年由县级自然保护区升级为省级自然保护区，2003年晋升为国家级自然保护区，2008年被列入《国际重要湿地目录》，总面积2360 hm2，其中核心区面积700 hm2，缓冲区460 hm2，实验区1200 hm2。区内潮间带湿地类型有滩涂、红树林和互花米草，另外区内还有水体和养殖池两种土地覆盖类型（非潮间带湿地）。红树林作为当地重点保护的资源，集中分布在核心区南侧，主要树种有秋茄、桐花树、白骨壤等。近年来，该区互花米草入侵形势严峻，扩张速度和面积均居国内前列。

图1 研究区与验证点分布Fig.1 Study area and field samples

2.2 数据获取与处理

2.2.1 Sentinel-2影像

Sentinel-2包括两颗极地轨道卫星（2015-06发射的Sentinel-2A和2017-03发射的Sentinel-2B），单颗卫星的重访周期为10 d，组合的双星重访周期为3—5 d。Sentinel-2 卫星搭载的多光谱成像仪（MSI），包含可见光、近红外至短波红外的13个光谱波段，波段参数信息见表1。快速重访和丰富的光谱信息为连续、快速、大范围地监测地球表面资源环境变化提供了新的契机。GEE 平台可直接访问 Sentinel-2 level-1C（L1C）和 level-2A（L2A）两个处理级别的数据产品。本研究所采用的Sentinel-2 数据都是L2A 产品，该产品主要包含经过大气校正的地表反射率数据。

表1 Sentinel-2 MSI波段参数信息Table 1 Sentinel-2 band parameter information

2.2.2 实地调查数据

实地调查工作开展于2020-11-12—25。研究人员使用RTK-GPS 深入研究区进行土地覆盖类型采样，记录相应位置的地物信息，包括土地覆盖类型、植被类型、树高、潮汐条件等。对于人力难以抵达的地区，采用无人机航拍进行调查。精灵PHANTOM 4 RTK 飞行拍摄的航高200 m，航向重叠90%，旁向重叠70%，采集的图像空间分辨率5 cm。研究区北岸的验证点采集于高空间分辨率谷歌地球影像。最终，采集湿地调查样点200 个，其中红树林、互花米草、滩涂各60 个，水体20个，空间分布见图1（b）。

3 研究方法

3.1 总体思路

图2是本文构建的自动、快速、高精度提取滩涂、互花米草和红树林的方法流程，主要包括以下几个步骤：（1）基于GEE云平台，构建高质量密集时序Sentinel-2 影像堆栈（详见3.3 节）；（2）深入分析潮间带湿地的时间序列变化特征，选取适合的光谱指数，基于MSIC 算法，从高质量时序数据集中分别合成最大水面、最小水面影像以及植被差异增强影像（详见3.4 节）；（3）基于Otsu 算法的多层决策树自动分类模型，由简单到复杂分层提取不同潮间带湿地（详见3.5节）。

图2 潮间带湿地分类流程图Fig.2 Workflow of intertidal wetland classification

3.2 MSIC算法和Otsu算法

GEE API中提供了“imageCollection.qualityMosaic（）”，即最大光谱指数合成算法（MSIC）。MSIC通过选用一个质量波段作为逐像素排序函数，提取该像元指定波段叠加排序后最大值的数据，质量波段作为最终合成影像中每个像素的像素源。进一步详细解释可见https：//developers.google.com/earth-engine/guides/ic_composite_mosaic［2021-05-21］。

Otsu算法，又名最大类间方差法（汪政辉等，2019），是Otsu（1975）提出的一种自动选择阈值的非参数方法。根据影像的灰度特性，将图像分割成背景和前景两类，选择使其类间方差最大、类内方差最小的分割阈值为最优阈值，进行图像的自动二值分类。

3.3 构建高质量密集时序Sentinel-2影像堆栈

构建高质量密集时序Sentinel-2 影像堆栈，实质是去除时间序列影像集合中的云像素，流程如下：首先，在GEE 中调用2020年全年所有覆盖研究区的Sentinel-2 L2A 产品，建立研究区初始密集时间序列影像集合；然后，利用QA60 波段去识别和掩膜初始密集时间序列集合中的厚云和卷云像素以获得高质量观测。新的时间序列影像集合称为高质量密集时序Sentinel-2 影像堆栈。以上步骤全部在GEE 平台完成，经过计算，2020年漳江口红树林保护区内每个像素的高质量观测次数（即有效观测数量）为53—68（图3，即2020年每个像素所在空间位置有53—68次有效重访数据）。

图3 2020年研究区有效观测数量的空间分布Fig.3 Spatial distribution of the number of effective observations over the study area in 2020

3.4 影像合成

根据Jia等（2021）和Duan等（2020）的研究结果，改进的归一化差异水体指数mNDWI（modified Normalized Difference Water Index）增强了堤坝、道路、池塘外缘等人工岸线与水体的差异，在提取潮间带和内陆地表水覆盖方面优于归一化差异水体指数NDWI（Normalized Difference Water Index）。为了抑制养殖池等其他类型背景水体的影响，本研究选用mNDWI 进行水体信息的提取。另外，非水（植被和滩涂）像素的归一化差分植被指数NDVI（Normalized Difference Vegetation Index）值高于水像素的NDVI 值，创建NDVI 最大光谱指数合成影像，可用于提取最小水面（Jia 等，2021）。该合成影像不仅可以增强植被和水体的差异，还可以去除海水周期性淹没对滩涂提取的影响，有效地将滩涂与水体分离。本研究直接采用mNDWI和NDVI 这两个指数构建MSIC 合成影像，生成最大（图4（a））和最小水面（图4（c））影像。本研究中，最大水面影像获取的是全年中研究区水面的最大覆盖范围，即最高潮；最小水面影像获取的是全年研究区水面的最小覆盖范围，即最低潮。

图4 潮间带合成影像及湿地自动分类各步骤结果Fig.4 Composite image of intertidal zone and results of each step of automatic wetland classification

在区分红树林和互花米草方面，本研究对两种植被的NDVI 时间序列变化特征进行分析。首先，采用时间序列谐波分析法HANTS（Harmonic Analysis of Time Series）对红树林和互花米草的NDVI时序曲线进行滤波和插值，结果如图5所示。互花米草与红树林具有明显不同的NDVI时序特征。红树林是常绿木本植物，全年NDVI 值较高；互花米草的NDVI曲线存在明显的季节变化。6—10月，互花米草处于生长最旺盛状态，红树林和互花米草的NDVI 差异较小。春冬季节，互花米草NDVI值较低，与红树林存在明显差异。

图5 2020年研究区典型像素的NDVI时间剖面及谐波分析（HANTS）Fig.5 NDVI time profile and harmonic analysis of typical pixels in the research area in 2020（HANTS）

MSIC 算法合成的影像，每个像素均是所选光谱指数的最大值像素，因此直接采用NDVI 进行合成，合成影像中的红树林和互花米草像素值相似。值得声明的一点是，本研究进行MSIC 合成时均采用原始时序NDVI。因此，需要构建一个新的指数，该指数的最大值可以区分红树林和互花米草。通过分析这两种湿地类型的NDVI 时序曲线，本研究构建了负NDVI 指数。基于负NDVI 指数进行MSIC 合成时，合成影像中的每个像素是红树林和互花米草NDVI 值最低的像素，如图5所示，红树林的NDVI 最低值和互花米草的NDVI 最低值相差较大。因此，该指数的MSIC 合成影像可以较好的区分这两种植被，该合成影像命名为植被差异增强影像（图4（f））。

3.5 基于Otsu的潮间带湿地自动分类

3.5.1 滩涂和水面分类

在本研究中，将滩涂定义为高潮与低潮之间的淤泥质或沙质无植被海岸。首先，应用Otsu 算法对“最大水面影像”的mNDWI 波段进行二值分割，自动生成二元分类结果，即水体和非水体，获得了年度最大地表水覆盖。由于mNDWI 有效地抑制了这些水面周围的人工岸线，内陆水面与潮间带区域水面没有直接相连。考虑到滩涂与海水密切相接，因此，通过保留最大面积水体斑块的方法获得“潮间带最大水面”（图4（b）），即滩涂的潜在区域，以除去养殖池、内陆湖泊等水体对潮间带滩涂提取的影响。其次，用“潮间带最大水面”掩膜提取“最小水面影像”（图4（c）），得到“最大潮水淹没区”影像（包括滩涂、海水及淹没植被），通过Otsu算法选取NDVI阈值去除淹没植被覆盖区，获得仅剩“滩涂海水影像”（图4（d）），采用Otsu 算法对NDVI进行二值分割，可实现滩涂和水体的自动分类（图4（e））。

3.5.2 红树林与互花米草的分类

“最小水面影像”（图4（c））中，通过Otsu算法选取NDVI阈值，得到植被和非植被两种类型，然后用海岸线数据裁剪植被区，得到潮间带植被区（红树林和互花米草）。将选取负NDVI 值进行MSIC 合成，生成潮间带“植被差异增强影像”（图4（g））。最后，应用Otsu 算法对“植被差异增强影像”进行二值分割，选取负NDVI 做为阈值特征波段，实现红树林与互花米草的自动分类（图4（h））。

4 结果与讨论

4.1 漳江口红树林保护区潮间带湿地分类精度

本研究利用采集的验证样点对潮间带湿地分类结果进行验证，得到的混淆矩阵如表2所示。滩涂的制图精度最小，水体的用户精度最低，一方面是因为不能确保获取绝对最低潮、最高潮时期的Sentinel-2 数据，致使滩涂的漏分现象多于其他地物，另一方面因为滩涂与水体交界处既有少量滩涂又有部分海水，容易造成二者之间的混分。红树林和互花米草的制图精度和用户精度均达到95%以上，主要是由于NDVI 可以很好的突出植被与滩涂、水体的差异，减少了滩涂与水体对植被信息提取的干扰，且通过负NDVI 最大值合成影像提供的物候信息增强了红树林与互花米草的特征差异，漏分及错分现象明显减少。本研究总体精度为96.5%，Kappa 系数为0.95，表明本文分类结果和验证数据具有很好的一致性。

表2 混淆矩阵及精度分析Table 2 Classification confusion matrix and precision analysis

4.2 漳江口潮间带湿地空间分布特征

2020年漳江口红树林自然保护区内红树林、互花米草、滩涂和水体的面积分别为82.46 hm2、218.26 hm2、496.84 hm2和401.90 hm2。此处声明一点：由于水体经过“最大面积水体斑块提取”，最终获取的不是完整的潮间带水面，但该水体结果可以用来验证本研究精度。由图6可以看出，滩涂作为区域内最丰富的资源，分布于红树林、互花米草的外缘。漳江口红树林分布较为集中，主要分布于研究区西南岸边。而互花米草分布于红树林的外缘与滩涂直接相连，在漳江口的南北两侧均有分布，南侧互花米草斑块面积大、整体性较好，北侧则沿海岸呈条带状分布。

图6 漳江口保护区内潮间带湿地分类结果Fig.6 Classification results of intertidal wetland in the Zhangjiang Estuary

4.3 MSIC和Otsu方法的优势

本研究基于GEE 云平台，构建高质量密集时序Sentinel-2 影像堆栈，结合MSIC 和Otsu 算法（MSIC-Otsu），实现了潮间带湿地信息的自动分类。本研究中潮间带湿地的成功分类主要归因于以下3个因素，即Sentinel-2数据的优势、稳健的MSICOtsu算法以及GEE海量的数据库与云计算能力。

Sentinel-2 数据的优势。在开源数据中，Sentinel-2 是时空分辨率最高的卫星数据。其较短的重访周期，为构建高质量密集时序影像堆栈提供了可能，大大增加了获取最高和最低潮时期影像的机会，为处于不断动态变化的潮间带湿地监测提供了巨大的潜力；更精细的空间分辨率不仅可以获取更多的空间细节，同时也减少了混合像元；其丰富的波段信息，也为潮间带湿地植被的监测提供了更广阔的应用潜力。

稳健的MSIC-Otsu 算法。本研究结合MSIC 和Otsu 算法进行潮间带湿地信息自动、快速、高精度的分类，该方法具有以下4 个优势：（1）MSICOtsu为全自动的分类方法，与现有的潮间带湿地分类方法相比，不依赖于训练样本、人工设定阈值等，鲁棒性强，具有广泛应用的潜力；（2）MSIC算法是基于像素合成影像的方法，即使某景影像存在云等噪音像素，也可合成高质量的影像，在一定程度上可以缓解滨海地区影像稀缺的影响；（3）由于滩涂和植被的NDVI 均大于云的NDVI，因此基于MSIC 算法合成的最低潮影像可以避免云噪声像素对其产生影响；（4）本研究选用MSIC 算法代替了传统的目视选取最低潮影像的方法，该算法通过选取指定波段最大值的像素合成一景影像，即每个像素均代表指定波段的极端值，不仅可以解决因潮汐作用引起的潮位不确定性问题，同时也可以提供充足的植被物候信息，为潮间带湿地分类提供有效手段。

GEE 海量的数据库与云计算能力。GEE 服务器支持用户快速访问、调用和分析在线存档的Sentinel-2 数据集产品，并采用实时分布式并行计算模型（Gorelick 等，2017），提高了潮间带湿地分类的时效性。

4.4 研究方法的不足

首先，尽管构建了高质量的密集时序影像堆栈，但是很难确保卫星会在绝对最低潮和最高潮位时拍摄影像，这也是本研究能否对潮间带湿地进行精确分类的重要影响因素之一；其次，由于在GEE 平台仅可获取2017-03 以后的Sentinel-2 L2A 产品，因此无法提供更为广泛的历史观测数据，且受制于Sentinel-2 系列数据的空间分辨率，难以识别低于10 m 的斑块；再者，本研究不能确保经过云掩膜后，每个像元都含有充足的无云观测以获取地物的关键信息；此外，mNDWI 并不能完全抑制过湿甚至常年被海水淹没的养殖池堤坝区域，这也是本研究部分海水养殖池存在错误分类的原因。

5 结论

本文基于GEE 平台构建高质量密集时序Sentinel-2 影像堆栈，深入分析潮间带湿地的时间序列变化特征，选取适合的光谱指数，结合MSIC和Otsu 算法实现了漳江口潮间带湿地的快速自动分类，主要得出以下结论：

（1）本研究构建的MSIC-Otsu算法得到的漳江口红树林自然保护区内的潮间带湿地的总体分类精度为为96.5%，Kappa 系数为0.95。2020年漳江口红树林自然保护区内，红树林、互花米草和滩涂的面积分别为82.46 hm2、218.26 hm2和496.84 hm2。

（2）高质量密集时序Sentinel-2 影像堆栈，一方面增加了获得最高潮、最低潮影像的机会，并可以自动、快速的合成最高潮和最低潮影像；另一方面，为红树林与互花米草分类提供了关键物候信息。

（3）结合mNDWI 最大值合成影像和提取最大面积斑块的方法，在抑制养殖池、堤坝及内陆水体方面取得了较好的效果。NDVI 最大值合成影像很好地突出了滩涂与植被、水体的差异；负NDVI最大值合成影像增强了红树林与互花米草的特征差异，为二者的分类奠定了基础。基于Otsu 算法的多层决策树自动分类模型，实现了研究区潮间带湿地的快速自动分类。

志 谢感谢厦门大学博士生郭先仙，中国科学院东北地理与农业生态研究所博士生张蓉，吉林大学硕士生闻馨，延边大学硕士生陈高在野外调查中的辛苦工作！