基于GEE的杭州湾慈溪段潮滩提取及时空变化分析

郑修诚， 周斌， 雷惠， 黄祺宇， 叶浩林

(1.杭州师范大学遥感与地球科学研究院，杭州 311121； 2.浙江省城市湿地与区域变化研究重点实验室，杭州 311121； 3.杭州师范大学信息科学与技术学院，杭州 311121)

0 引言

潮滩是指淤泥质海岸潮间带浅滩[1]，处在海陆交汇的敏感地带，是沿海城市发展的重要物质基础和后备资源，也是沿海地区生态安全体系不可或缺的一部分，因此研究潮滩的发展和变化具有重要的意义[2]。由于潮滩范围广，变化快，并且在潮滩区域行动极为不便，部分区域难以到达，使用传统的实地勘测法调查潮滩具有明显的复杂性，勘测效率低下，难以实现动态监测，而遥感技术则因其具有大范围、高分辨率、多光谱和多时序等优势，非常适用于潮滩的动态监测。

目前基于遥感技术的潮滩提取方法主要有“八分算潮法”[3]和“相似三角形原理”[4]等。张春桂[3]采用“八分算潮法”来确定高、低潮时，选用过境时间与这些时刻最为接近的MODIS数据，实现了福建省海岸带潮滩的提取，但MODIS数据空间分辨率偏低且该方法只是一种潮时的近似估算方法，存在一定的误差； 王小龙等[4]采用“相似三角形原理”，结合海岛多年的潮汐数据，利用高分辨率遥感数据计算了东沙岛潮滩和湿地的范围，该方法将潮滩地形剖面视作一个直角三角形，是一种理想的估算方法，王小丹等[5]研究曹妃甸地区潮滩演变时使用的“相同潮位对比法”以及王靖雯等[6]和韩倩倩等[7]研究中使用的“潮位校正法”也都以此方法为基础。总体来看，当前常见的潮滩遥感提取方法往往通过估算来确定潮滩的边界，难以保证较高的提取精度。而张媛媛等[8]采用穷举法，应用多时相环境卫星数据作为数据源提取序列水边线，通过面向对象分类，实现了江苏省如东县潮滩的提取，该方法摆脱了对高、低潮线的“估算”，使用若干条水边线拟合高、低潮线，具备较高的精度。

近年来，随着Google Earth Engine(GEE)平台的诞生和广泛应用，穷举法所面临的数据量大、依赖人工等问题迎刃而解，本研究基于GEE平台，选用时序覆盖较广的Landsat系列影像数据，目视解译人工海岸线作为平均高潮线，利用水边线拟合平均低潮线，实现了对杭州湾潮滩变化最为显著的慈溪段潮滩区域的提取，并分析了其变化情况，为该地区潮滩湿地的保护提供技术参考。

1 研究区概况与数据源

1.1 研究区概况

以杭州湾南岸的慈溪市为研究区(图1)，其隶属于浙江省宁波市，东接镇海区，西连余姚市，地势南高北低，面向杭州湾呈现丘陵、平原、滩涂3级阶梯状。由于杭州湾独特的空间形态，大量的泥沙在杭州湾南岸北凸弧段淤积，形成了典型的平原型潮滩，为慈溪市填海造陆提供了得天独厚的自然条件，使其成为浙江省土地后备资源最富足的地区之一。

图1 研究区地理位置

慈溪市也具有悠久的填海造陆历史，从明末至20世纪末，慈溪市海岸线外移的速度为28 m·a-1[9]。20世纪90年代以后，围涂筑塘技术趋于成熟，施工设备和施工工艺逐步先进，围塘砌筑标准级别提高，单块围涂规模不断扩大[10]。

1.2 数据源

GEE是由Google提供的基于云计算的全球尺度地理空间分析平台，其存储的大量公开的地理空间数据集能够省去影像下载的时间，而其强大的计算力则能批量处理较大规模的数据，为较大时空尺度的研究提供便捷。出于对数据可获得性和逐年监测的需求，本研究选用30 m空间分辨率的Landsat系列数据作为遥感影像数据源，所有数据均来自于Earth Engine Data Catalog。本研究对数据源进行了2轮筛选。首先筛选覆盖完整研究区(条带号118/39)的1990—2021年间云量低于20%的影像，累计得到142景影像； 再对这些影像逐景目视筛选，确保潮滩区域无明显厚云覆盖，并尽量保证1 a内有2景间隔90 d以上的影像，最终得到77景符合条件的影像(图2)。

(a) 年际分布(b) 月际分布

行政区矢量边界数据获取自中国科学院资源环境科学与数据中心(https: //www.resdc.cn/)。

研究区潮汐预报数据获取自国家海洋科学数据中心(http: //mds.nmdis.org.cn/)。

2 研究方法

本文研究的潮滩范围是介于平均大潮高潮线和平均大潮低潮线之间的潮侵地带，结合前人的相关研究[8,11]，平均大潮低潮线可以认为是水边线集中距离陆地最远的外界线，这通常不是集中的某一条线，而是多条线段的组合； 而海岸线则可以被认为是多年平均大潮高潮位形成的痕迹线[12-14]。平均高、低潮线合围形成的区域就是本文的潮滩区域。结合现有数据，本文拟提取1990—2021年间5～6 a时间间隔的6期潮滩区域。

2.1 瞬时水边线提取

瞬时水边线的提取本质上是对数字图像中的陆地区域和水体区域进行分割，2类区域的分界线就是瞬时水边线。本研究主要使用基于水体指数的阈值分割法提取水边线，其中阈值的选取起着至关重要的作用[15]。通常经过水体指数计算得到的灰度图像直方图会出现“双峰”特性[16]，选取谷底所对应的灰度值作为阈值能够有效地对灰度图像进行分割。再对阈值分割后的二值图像矢量化。技术路线如图 3所示。

图3 潮滩提取技术路线

2.1.1 自动水体提取指数AWEI

2014年，Feyisa等[17]提出了一种基于TM多波段的自动水体提取指数(automated water extraction index, AWEI)，并且证实了该指数相较于常用的改进的归一化差异水体指数(modified normalized difference water index，MNDWI)[18]具有更高的精度。AWEI通过给TM影像不同波段赋予相应的系数并进行加减计算实现水体像元和非水体像元之间最大程度的分割，该指数由2个独立的方程组成，公式分别为：

AWEIsh=ρBLUE+2.5ρGREEN-
1.5(ρNIR+ρSWIR1)-0.25ρSWIR2，

(1)

AWEInsh=4(ρGREEN-ρSWIR1)-
(0.25ρNIR+2.75ρSWIR2)，

(2)

式中:ρBLUE，ρGREEN，ρNIR，ρSWIR1和ρSWIR2分别代表Landsat影像蓝光、绿光、近红外、短波红外1和短波红外2波段的反射率值;AWEIsh是为研究场景下有明显阴影存在的情况设计，AWEIsh能够有效去除阴影像元[19]；AWEInsh则是为研究场景下阴影问题并不明显的情况设计[20],AWEInsh能够有效地去除易与水体混淆的黑色建筑地表。本研究聚焦的潮滩区域没有明显的阴影干扰，因此选用AWEInsh进行水边线提取。经过实验，AWEInsh同样适用于Landsat8 OLI数据[19]。

2.1.2 大津算法提取阈值

大津算法(OSTU)[21]是日本学者大津于1979年提出的一种确定图像二值化分割阈值的算法，从原理上讲，用该方法提取的阈值对图像进行二值化分割后，前景与背景图像的类间方差最大，因此该方法又被称为最大类间方差法。通过在GEE平台中编写代码，可以自动获取经过水体指数计算的每景灰度图像的阈值，并对其进行二值化分割。

2.2 平均低潮线提取

受潮汐作用、泥沙淤积等因素的影响，瞬时水边线之间会出现交叉现象，并且随着泥沙的不断淤积、地形不断改变，时间间隔越长、水边线数量越多，交叉现象就越复杂[8]，通过目视解译的方法从水边线集中提取平均低潮线就需要耗费大量的时间，并且受人为干扰较大，因此本研究使用基于ArcGIS软件开发的数字岸线分析系统(digital shoreline analysis system, DSAS)[22]，通过提取离陆地最远的潮位点，拟合一条平均低潮线。

DSAS模块通过对一条Baseline作若干条间隔相同、长度固定并与水边线集充分相交的垂线，记录垂线与水边线集的所有交点数据，这些交点包含了其所在垂线编号、至Baseline距离以及空间坐标等信息。经过多次测试，结合研究区实际情况，选定垂线间隔为500 m，垂线长度为18 000 m，最后通过Python编程对交点数据进行分析，筛选出同一垂线编号中距离Baseline最远的潮位点，将以上潮位点按其属性中的垂线编号依次相连，拟合生成的折线即认为是平均低潮线(图 4)。共生成1995年、2000年、2005年、2010年、2015年和2021年6期平均低潮线。

(a) DSAS模块提取低潮点(b) 低潮线生成

2.3 平均高潮线提取

慈溪市有悠久的围海筑堤历史，域内以人工海岸线为主，因此本研究对平均高潮线的提取实质上是对人工海岸线的提取。人工海岸线在遥感影像中纹理清晰、辨识度高，因此本文采用目视解译法直接提取。结合研究区Landsat影像和Google Earth历史影像，遵循海岸线判别依据[23]，在ArcGIS 10.2软件中绘制与平均低潮线对应的6期海岸线，作为本研究所需的平均高潮线。

2.4 潮滩区域提取及精度评价

潮滩提取实质上是将研究区内地物分为潮滩和其他地物2类，而精度评价是检验分类结果与地表真实信息吻合程度的过程，以此来评价分类结果可信度。误差矩阵又称混淆矩阵，广泛应用于遥感土地利用分类精度评价，能简单地对分类精度信息进行概括[24]。

本研究通过构建缓冲区、创建随机采样点，统计每个采样点在提取结果和影像数据中分别是否属于潮滩，并构建统计结果的混淆矩阵从而对潮滩提取的精度进行评价。为使采样点分布合理、覆盖均匀，以提取的潮滩结果为中心作半径5 km的缓冲区，作为提取结果精度评价的区域(图5)。每期缓冲区内使用Create Random Point工具随机生成240个点作为采样点，每个采样点的间距不小于100 m。借助研

图5 缓冲区和采样点构建

究区低潮期遥感影像以及Google Earth历史影像，目视解译将采样点分为潮滩和其他2类，最后构建混淆矩阵评价潮滩提取精度。

2.5 平均高、低潮线迁移速率计算

对于2.2节中的任一垂线段i，计算其与某2期高(低)潮线的交点间的距离di，作为该时期内在该垂线段上的潮线位移距离。任意2期高(低)潮线的平均迁移速率V的计算公式为：

(3)

式中：n为生成垂线段总数；Δt为2期高(低)潮线间的年份差值。

3 结果与分析

3.1 潮滩提取结果及精度评价

3.1.1 潮滩提取结果

图6反映了1990—2021年间6期潮滩的空间位置，底图为1990年Google Earth历史影像，可见潮滩区域呈现明显向北发展的态势。依次统计各期潮滩的面积(图7)，通过图 7可知，研究区内潮滩的面积大部分时间维持在20 000～24 000 hm2区间，在2000—2005年间陡然减少至10 344 hm2，随后又恢复到20 000 hm2以上。

(a) 1995年(b) 2000年(c) 2005年

(d) 2010年(e) 2015年(f) 2021年

图7 潮滩面积统计

3.1.2 精度评价

根据目视分类结果，构建分类结果混淆矩阵(表1)，并统计各期潮滩精度(表2)。总体精度表示所有正确分类的土地覆盖类别的检验点数所占总抽取的检核点数的百分比； 用户精度表示分类结果中样本点的类别与地面实际类别相同的条件概率﹔生产者精度表示地面获得的实际资料参考点类别与分类结果中对应点类别相一致的条件概率； Kappa系数是1960年由Cohen首先提出的一种应用于遥感影像分类结果评价的一致性检验方法﹐一般而言，Kappa值介于0～1之间，Kappa值越大表示分类精度越高(表2)。

表1 潮滩提取精度评价混淆矩阵(2021年)

表2 潮滩提取精度评价结果汇总

由表2可见本研究提取的6期潮滩总体精度都在90%以上，并且Kappa系数所反映的分类精度都达到了最佳[25]，表明本研究方法能有效地提取潮滩区域。

3.2 潮滩时空变化分析

3.2.1 潮滩地理中心空间变化

逐一计算6期潮滩区域的地理中心，并统计潮滩中心点位移情况，具体如图 8和表3所示。

图8 潮滩地理中心时空变化

表3 潮滩地理中心迁移情况

从中心迁移轨迹(图8)来看，潮滩区域地理中心经历了“西北-东北-西北-东南-西北”的复杂迁移过程，但总体呈现由东南向西北迁移的趋势，累计迁移距离为28 083 m，总迁移距离为8 127 m。从各个迁移方向上看，西移发生在1995—2000年、2005—2010年和2015—2010年，西移距离分别为701 m,8 893 m和171 m，速度分别为140.2 m·a-1,1 778.6 m·a-1和34.2 m·a-1； 东移发生在2000—2005年和2010—2015年，距离分别为1 614 m和4 926 m，速度分别为322.8 m·a-1和821 m·a-1； 南移只发生在2010—2015年，距离为1 150 m，速度为230 m·a-1，这主要是由于2010—2015年间研究区西侧海岸线外移潮滩减少，而东南侧潮滩增加明显，导致潮滩地理中心向东南方向迁移； 北移连续发生在1995—2010年和2015—2021年，距离分别为7 370 m和1 240 m，速度分别为491.3 m·a-1和206.7 m·a-1。

总体上，地理中心向西迁移了3 225 m，速度为124 m·a-1； 向北迁移了7 460 m，速度为286.9 m·a-1，南北方向迁移的跨度距离大于东西方向的迁移跨度表明潮滩南北方向变化大于东西方向变化。从各个迁移阶段看，2005—2010年在东西和南北方向上的迁移距离均为最大，说明这5 a间潮滩发生了最为明显的变化。

3.2.2 高、低潮线空间变化

表4显示，除2000—2005年低潮线向内退缩外，其余时段高、低潮线均以不同速率外扩。高潮线外扩最快的时段是2000—2005年，达到233.87 m·a-1； 最慢的时段是2015—2021年，仅为23.65 m·a-1； 低潮线外扩最快的时段是2005—2010年，达到461.52 m·a-1。

表4 高(低)潮线迁移速率

3.2.3 潮滩面积及空间变化原因分析

研究发现，30余年间，杭州湾南岸慈溪段潮滩基本保持北移的态势，但是潮滩面积基本维持在相对稳定的水平，仅在2000—2005年间出现骤减，对比图6和表4可知，此时段内海岸线发生了明显外移，但潮滩北界相较于2000年并未明显外移，西侧潮滩还出现了消退的情况，以上因素综合导致了2005年的潮滩面积大幅减少。张华国等[23]研究了1986年以来杭州湾围垦淤涨状况，发现2000年以后每年围垦面积增长的区域数量急剧上升，并推测在2001—2003年间有一轮较大规模的围垦工程。2001年，慈溪市作出了开发建设杭州湾新区的战略决策[26]； 2003年，杭州湾跨海大桥开始奠基建设，大桥又恰好位于潮滩消退明显的西侧。这一系列涉及该区域的重大事件，可能是造成慈溪市2000—2005年海岸线剧烈外移、潮滩显著减少的原因。

而此后的5 a间，潮滩则以461.52 m·a-1的速率快速淤涨，因此，海岸线的剧烈外移可能会在未来一定时期内推动潮滩的淤涨。孙超等[27]研究江苏中部沿海盐沼演变与围垦的关系时也认为“围垦活动能够改变原潮滩的沉积环境”，推动盐沼扩张。

另外，2010—2021年间潮滩面积保持在较稳定的水平，2015—2021年潮滩地理中心和高、低潮线迁移也较小，可见近10 a，尤其是2015年之后针对潮滩的开发利用强度大幅减弱。近年来，国家为了应对海洋生态环境日益严峻的形势和海洋开发秩序混乱等问题，相继出台政策管控围填海活动。2018年，国务院发布了《国务院关于加强滨海湿地保护严格管控围填海的通知》[28]，严格限制了“向海索地”。因此，在未来一定时期内，研究区内潮滩在空间形态和空间位置上很可能以稳定的速度缓慢发展。

4 讨论

在本研究中，影像成像时的瞬时潮位会对潮滩提取的精度产生直接影响。潮滩提取精度主要取决于低潮线提取的精度，这就需要研究所使用的影像拍摄时间处于该海域的低潮期，因此需要对影像和其拍摄时间的潮位进行星地数据匹配。潘存鸿等[29-30]在研究杭州湾潮汐特征的研究时介绍了杭州湾沿岸潮位站分布： 杭州湾北岸设有澉浦、乍浦、金山嘴和芦潮港4个长期潮位站，杭州湾南岸仅有镇海口1个长期潮位站。这些潮位站均未在研究区范围内，其提供的潮汐实测数据不能客观地反映研究区潮汐特征。而国家海洋信息中心的全球潮汐预报服务平台(http: //global-tide.nmdis.org.cn/)提供了1980—2022年全球485个主要港口的潮汐预报数据，其中“海黄山”站点(121.50°E, 30.21°N)处于研究区核心地带，可以作为相对客观的潮位参考数据。本研究统计了所使用的每景影像拍摄时间对应的预报潮高，并对多年潮汐数据做了简单分析，发现该站点的低潮潮高通常在80～120 cm之间，经统计，共有14景影像成像时的预报潮位低于120 cm，并在本研究设定的6个时段内都有分布，因此本研究提取的潮滩具有一定的可信度。

然而，本研究主要基于“云量”因素筛选影像，最终完成潮滩的提取，是存在一定随机性的。即便低云量的数据量非常丰富，但影像拍摄时间未在低潮期，提取出的水边线也无法认为是低潮线。因此，可以在筛选影像前进行星地数据的匹配，并以此作为影像筛选的依据，例如Landsat系列卫星在本研究区的过境时间约为北京时间上午10: 00前后，仅需要统计低潮潮时在上午10: 00左右的日期，再与卫星过境日期进行匹配，便可高效筛选目标影像，并大幅降低低潮线提取的不确定性。

当然，这一定程度上还依赖于遥感卫星的重访周期。若非Landsat系列数据在时间尺度上覆盖了此前的几十年，其16 d的重访周期会成为星地数据匹配的消极因素。近年来，随着国产卫星遥感数据(如高分、资源、环境等系列数据)获取愈发便捷和Sentinel系列数据可免费获取，大大弥补了原本匮乏的数据源。尽管不同数据间由于传感器设计不同，几何位置和波段设置上存在差异，多源数据交叉使用时，数据匹配存在一定困难，但遥感技术在潮滩监测领域的应用仍具备非常光明的前景。

5 结论

本研究结合GEE遥感云计算平台和GIS技术，选用1990—2021年间共77景Landsat卫星影像，通过提取水边线和海岸线，实现了对杭州湾南岸慈溪段的潮滩提取和面积估算，并对潮滩区域的变化情况展开了分析，为慈溪潮滩湿地的管理与保护提供了技术参考。取得的主要结论如下：

1)本研究使用Landsat卫星影像作为数据源，在保证影像质量的同时满足了较长时序监测的需要； 使用阈值分割法提取了影像的瞬时水边线； 通过DSAS以及Python语言编程将研究对象从线要素转变为点要素，再由经过筛选的点要素拟合出研究所需的平均低潮线； 基于潮滩的定义，将平均低潮线与海岸线合围，提取出了1995年、2000年、2005年、2010年、2015年和2021年6期潮滩区域，并且达到了较高的精度，说明该方法适用于本研究。

2)本研究大部分影像数据处理工作在GEE遥感云计算平台完成，自动化程度高，对数据量较大的时序监测研究非常友好； DSAS模块通过对研究对象要素的转变，较为科学客观地实现了低潮线的提取； 在海岸线信息比较明晰的(基岩、人工海岸)区域，可以比较便捷地提取出潮滩区域，具备一定的推广性。

3)本研究发现杭州湾南岸慈溪段潮滩面积基本维持在20 000～24 000 hm2区间，2000—2005年间陡然减少至10 344 hm2，2010年又恢复至原先水平； 潮滩空间变化趋势是由南向北迁移，迁移速度为286.9 m·a-1，其中2005—2010年潮滩的空间变化最为明显。

4)本研究认为沿海围垦和潮滩淤涨共同影响了潮滩面积和空间的变化，大范围的沿海围垦可能在未来一定时期内促进潮滩的淤涨，而这一切最主要的驱动力是地方政策，随着针对围填海活动的严控政策相继出台，潮滩将保持稳定的发展态势。

当然，基于“云量”因素筛选影像具有一定的随机性。在后续的研究中将结合潮汐预报数据，进行星地数据匹配，从而实现精准、高效的影像筛选，并尝试引入更高时空分辨率的多源数据，实现更高精度的潮滩提取和监测。