基于多时相多光谱遥感影像的珊瑚礁面积估算方法研究−以西沙群岛羚羊礁为例

熊媛，黄荣永, 4*，余克服, 4

( 1.广西大学 海洋学院，广西 南宁 530004；2.广西南海珊瑚礁研究重点实验室，广西 南宁 530004；3.广西大学 珊瑚礁研究中心，广西 南宁 430004；4.南方海洋科学与工程广东省实验室（珠海），广东 珠海 519080)

1 引言

珊瑚礁具有极高的初级生产力和丰富的生物多样性，是南海最重要、最具特色的生态系统，在南海生物多样性维护、生态资源供给等方面发挥着极其重要的作用；珊瑚礁也是南海最主要的陆地国土类型，是我国渔业安全生产活动的重要基地[1]。然而，受人类活动和气候变化的影响，南海珊瑚礁退化严重，从活珊瑚覆盖度来看，退化幅度达80%[2]。2015年，西沙群岛的活珊瑚覆盖度仅约16.3%[3]，黄岩岛礁坪的面积在1997−2019年间则以11 289 m2/a的平均速率减少[4]。珊瑚礁的退化严重影响到了南海的生态安全。

分布和面积是珊瑚礁的最基本属性。珊瑚礁面积是了解珊瑚礁资源特征[5]、量化和评价其经济价值[6]、计算碳收支和钙化率[7]的重要指标；是探讨珊瑚礁的动态变化过程和预测珊瑚礁未来发展趋势等的基础[8]；是开展珊瑚礁资源规划、管理和保护的前提[9]。对于我国的珊瑚礁面积来说，Hughes等[10]根据其他学者的研究成果汇总得出我国珊瑚礁的面积保守估计有30 000 km2，这与Spalding等[11]估计的2 450 km2相差甚远；赵焕庭[12]汇总得到南海珊瑚礁总面积则约有37 200 km2；Wang和Li[13]估算南海浅水区（水深小于50 m）现代珊瑚礁面积约为8 000 km2，大亚湾、福建东山的造礁石珊瑚群落和珠江口的珊瑚礁群落，因其不能形成珊瑚礁，只能称之为造礁石珊瑚群落[14]，故未被计算在其中。东沙群岛、西沙群岛、中沙群岛和南沙群岛的珊瑚礁面积数据文献[15-17]均有涉及。Dai[18]根据其之前的研究数据得出东沙群岛的面积约为600 km2，文章未明确说明数据来源。对于南沙群岛珊瑚礁的面积，Hughes等[10]认为为26 059 km2，而Spalding等[11]估算为1 150 km2。由此可见，我国珊瑚礁的面积至今尚未形成共识。目前有关南海珊瑚礁的面积资料来源仍然以20世纪报道为主[15,19-22]，但当时的珊瑚礁面积估算基本上是基于海图，并没有相对全面的现场调查或比较清晰的图片，也没有比较好的珊瑚礁面积估算方法，因此迄今对南海珊瑚礁的面积仍然没有形成统一的认识。

遥感是迄今为数不多的能够非接触地对南海珊瑚礁进行大范围同步测量，且时间和空间分辨率均较高的低成本探测手段[9,23]。现有的遥感估算珊瑚礁面积的方法主要是图像分类法（基于训练样本的机器学习和基于先验知识的监督分类方法），即通过统计图像分类像元个数来进行珊瑚礁面积的估算。例如，美国千年珊瑚礁研究计划[16,24]使用半自动化基于图像分析和结合历史海图等方法对珊瑚礁进行分类与面积计算；美国国家海岸带海洋科学中心（National Centers for Coastal Ocean Science, NCCOS：https://coastalscience.noaa.gov/products-explorer/）利用美国国家海洋和大气管理局（National Oceanic and Atmospheric Administration, NOAA）栖息地数字化仪扩展对IKONOS影像进行目视解译从而创建珊瑚礁地图集；哈立德·本·苏丹海洋生物基金会的全球珊瑚礁探测项目（KSLOF-GRE）[25]基于米级高分辨率遥感图像的评估，采用面向对象和随机森林的方法对珊瑚礁进行识别分类；艾伦珊瑚礁地图集[17]采用基于对象的图像分析和生态建模结合的方式进行珊瑚礁的面积估算等。首先，这些研究对于珊瑚礁面积测算方法和准则比较模糊且不统一，导致不同珊瑚礁面积估算结果会存在比较大的差异；其次，NCCOS数据集和KSLOF-GRE数据集未包含南海珊瑚礁的信息；第三，美国千年珊瑚礁研究计划和艾伦珊瑚礁地图集数据集基本上忽略了礁前斜坡、潟湖和暗礁的存在，进而导致珊瑚礁面积被低估，因此其精度和准确性仍存在较大的问题。此外，KSLOF-GRE采用高分辨率遥感影像的方式虽然有利于提高珊瑚礁识别精度，但其成本较高，且解译时间较长，使得优于10 m的更高分辨率的遥感影像通常难以用于估算较大范围珊瑚礁的面积。

在国内，朱海天等[26]、李成鹏等[27]和董娟等[28]结合珊瑚礁不同地貌带的异质性，利用面向对象法分割影像获得了南海珊瑚礁地貌的整体信息；霍艳辉[29]基于固定阈值分割方法提取三亚珊瑚礁的分布信息。这些方法虽然原则上能通过统计其分割的地貌像素个数获得珊瑚礁面积，但是，这些基于像素的地貌边缘提取比较粗糙，影像分辨率的变化仍然会对面积估算产生较大的影响。Liu等[8]在关于典型灰沙岛的研究中曾涉及南海珊瑚礁部分地貌单元面积变化的估计问题，但亦尚未涉及珊瑚礁整体面积的概念界定与科学计算。即使Dong等[30]曾提出“淹没频次”的原理，并据此利用Landsat 8影像通过模糊C均值聚类法（Fuzzy C-means Algorithm, FCM）计算了南沙群岛的137个珊瑚岛礁/沙洲/暗沙/浅滩面积。该方法仍较简单地基于像素统计方式来估算珊瑚礁的面积，亦尚未能够理清其所估算的面积与珊瑚礁地貌领域的对应关系。

因此，就目前的文献来看，国内尚未发现有关南海珊瑚礁面积遥感估算的直接探讨，仅出现较为分散的部分相关技术的探讨。

研究区的环境条件是遥感影像解译的一个重要限制因素[31-32]。利用遥感影像估算珊瑚礁面积的主要依据是珊瑚礁与周围海水之间存在的亮度差异，但诸如潮汐引起的水深变化和海水组分引起的透明度变化之类的环境因素都会导致遥感影像的这一亮度差异发生变化[33-34]，进而导致遥感影像中珊瑚礁与周围海水之间分界线的位置往往存在较大的不确定性。因此，珊瑚礁遥感影像识别的结果通常会存在较大的误差[35-36]，最终造成珊瑚礁面积遥感估算的不确定性。

本文选择南海西沙群岛的羚羊礁为研究对象，基于多时相多光谱卫星遥感影像，结合现场测量与调查资料，提出了一种半自动且低成本的珊瑚礁面积估算方法。具体而言，本文首先利用数字图像边缘检测算法对珊瑚礁地貌带分界线位置进行自动精细化提取，然后再对精细化提取得到的边界线进行多时相融合以便提高面积测算结果的精度和可靠性。由于引入数字图像处理理论中的图像边缘检测算法，所提出的方法能够比较好地降低遥感影像分辨率所带来的影响；而由于采用多时相遥感影像融合的方式，所提出的方法能够较好地减小潮位等水体性质变化所导致的面积估算偏差。该方法的提出将为珊瑚礁面积的计算提供技术参考。

2 研究区域和数据

2.1 研究区域

羚羊礁是西沙群岛的一个环礁，位于永乐环礁西南部。礁体近似水滴状，无口门，礁盘东南边缘有沙洲发育，高出水面约2 m。其礁体处于16°26′～16°30′N，111°34′～111°37′E，南北宽约 6 km，东西长约4 km，如图1所示。

图1 研究区域① 底图来自自然资源部标准地图服务系统（http://hism.mnr.gov.cn/sjkf/bzdt/201902/t20190214_3124659.html），审图号为琼S（2020）038号。Fig.1 The study area

2.2 研究数据

2.2.1 遥感影像

本文采用Sentinel-2卫星2019年前的L1C级和2019−2021年的L2A级多光谱影像进行羚羊礁面积估算，而Landsat 8和WorldView-2多光谱影像则主要用于所提方案性能的比较分析。通过目视方法筛选出研究区域上空无云且海表无强烈太阳耀光影响的影像数据。其中，Sentinel-2 多光谱成像仪（Multispectral Scan Imaging, MSI）影像从欧空局−哥白尼项目（European Space Agency-Copernicus， ESA-Copernicus）的开源门户网站（https://scihub.copernicus.eu/）下载，分辨率为10 m，影像处理基线编号为0204，相对轨道编号为032，拼接域编号为T49QEU；Landsat 8影像则从美国地质调查局（U.S.Geological Survey,USGS）网站（https:// Earthexplorer.usgs.gov/）下载，分辨率为30 m，条带号为122，行编号为049。同时，购买了2014年10月9日永乐环礁的WorldView-2 高分辨率卫星遥感影像，其分辨率为1.8 m，影像编号为1030010036B4B000。所用遥感影像详细列表见表A1。

2.2.2 测深数据

利用羚羊礁的6条实测水深剖面数据（图2）与羚羊礁地貌带识别结果进行比较，以便评估本文地貌带分界线识别的位置精度。这6条实测水深剖面数据来源于国家科技基础条件平台—国家地球系统科学数据中心（http://www.geodata.cn）的南海中北部珊瑚礁本底调查项目（科技基础性工作专项，项目编号：2012FY112400）。由该项目说明文档可知，数据测量时间为2015年5月，采用断面测线连续定位测深，相邻两次测量的水平距离约为5 m；测深与定位仪器分别采用中海达HD370高频测深仪和K3DGPS信标机，测深精度优于0.3 m，定位精度优于10 m。这些实测水深剖面基本涵盖了羚羊礁的各种地貌带类型，如图2所示，W01−W03表示羚羊礁西部的3个断面，E01−E03表示羚羊礁东部的3个断面。

图2 实测水深点位置Fig.2 The position of actually measured water depth

3 研究方法

珊瑚礁面积估算方案如图3所示，由遥感影像预处理、地貌类型划分和珊瑚礁地貌制图与面积估算3个模块构成。

图3 用于珊瑚礁面积估算的工作流程Fig.3 Workflow for area estimation of coral reefs

3.1 遥感影像的预处理

Sentinel-2卫星的L1C级数据是经过正射校正和几何精度校正的大气表观反射率产品。为消除大气吸收和散射等作用，本文使用Sen2cor插件（http://step.esa.int）对2019年之前的Sentinel-2 MSI L1C级数据进行大气校正，以便获得L2A级产品，即地表反射率数据。Landsat 8和WorldView-2影像则通过ENVI软件开展辐射定标和大气校正。影像因海表波浪不规则起伏与太阳辐射的相互作用而产生太阳耀光。本文利用Hedley等[37]所提供的方法进行太阳耀光的消除。随后，以2021年4月24日的Sentinel-2 MSI影像为参考，对其他影像进行几何配准，配准误差控制在1个像元以内，以便保证不同时相影像之间具有共同的参考坐标系。

3.2 珊瑚礁地貌划分

本文将珊瑚礁面积测算的范围界定为理想条件下在光学遥感影像中可见的浅水珊瑚礁（30 m以浅）。30 m是清澈海水中可见光波段的较大透射深度[38]，而珊瑚礁也是在30 m以浅生长最佳[1]。相比于只将出水部分算作珊瑚礁面积来说，上述界定的范围包含了未出露水面的覆水部分，该部分在珊瑚礁生态系统中也起着关键作用[39]。南海珊瑚礁以环礁为主，根据余克服等[40]绘制的南沙群岛环礁的地貌沉积带图，珊瑚礁地貌包括礁前斜坡、礁坪、潟湖坡和潟湖，如图4所示。

图4 珊瑚礁的地貌分区Fig.4 Geographic zonation of coral reefs

（1）礁前斜坡

礁前斜坡位于理论最低潮面以下，由礁坪外缘向海倾斜，宽度为20～30 m，受破浪带作用，是珊瑚礁主要的地貌单元之一。礁前斜坡通常被珊瑚所覆盖，其上发育有珊瑚礁脊槽地貌。脊槽地貌的平行脊由规则间隔的沟槽分隔，形成一种特有的梳齿状图案。礁前斜坡反射率较低，遥感影像上有微弱信息，在图像上表现为较暗区域，呈暗蓝色。

（2）礁坪

礁坪位于理论最低潮面以上，是从礁坪外缘至潟湖坡之间的一个宽阔的地貌带，坡度较缓，是珊瑚生长和生物堆积的主要部位。礁坪基底布满珊瑚礁石和松散砂，广泛发育礁坑。礁坪反射率较高，在遥感图像上表现为较亮区域，呈亮褐色。

（3）潟湖坡

潟湖坡在理论最低潮面以下，但水深不大，是礁坪到潟湖之间的一个缓坡，分布有大量的珊瑚碎屑沉积物和砂砾。潟湖坡反射率较礁坪反射率高，在遥感图像上表现为较亮区域，呈亮蓝色。

（4）潟湖

位于潟湖盆沉积带，在理论最低潮面以下，是环礁的特有地貌。潟湖是一圈珊瑚礁所围成的水深不大的浅湖。潟湖中堆积生物碎屑，有不少点状珊瑚礁，即点礁。潟湖反射率最低，接近海水，在图像上表现为最暗区域。

这些区域之间的边界通过坡度变换线，即坡折线划分[41]，如图4所示的砂坡、礁坑、礁顶和前斜坡处。潟湖与潟湖坡的分界处是坡度较缓的砂坡或小礁崖；礁坪与礁前斜坡的分界处在礁顶；礁坪与潟湖坡的分界线处是小礁坑；同样的，礁前斜坡与外海的分界处在礁缘陡坡带。在遥感影像中，灰度梯度极大值往往出现在地貌分界处[42]，亦是坡折线处。事实上，潟湖坡靠近礁坪处生长有成片珊瑚[1]，礁坪与潟湖坡之间的坡折线易受其影响而难以区分。由于珊瑚礁面积估算不受礁坪与潟湖坡地貌分带的影响，故本文将礁坪与潟湖坡合并为一个地貌带，即“礁坪−潟湖坡”，从而进行地貌划分与面积计算。

3.3 珊瑚礁地貌分界线的提取

根据坡折线附近存在灰度跳跃变化的特点，利用Xu和Prince[43]提出的基于梯度向量场的主动轮廓线（Gradient Vector Flow-Snake, GVF-Snake）模型对珊瑚礁地貌分界线进行提取。GVF-Snake模型是一种基于偏微分方程矢量扩散的方法，其主要通过偏微分方程的矢量扩散来生成新的外力场，从而降低模型对初始化轮廓的敏感性，同时增强曲线收敛到目标边界凹面的能力。GVF-Snake模型已在相关领域中得到一些应用，例如周旻曦等[44]利用GVF-Snake模型对珊瑚礁地貌进行识别，减少了地貌分类中离散斑块的数量；刘嘉鎏等[42]通过GVF-Snake模型和数字化海岸线分析系统（Digital Shoreline Analysis System, DSAS）分析了南海黄岩岛环礁地貌近40年的变化。

本文涉及的Sentinel-2和Landsat 8多光谱影像均包括有蓝波段、绿波段和红波段3个波段。蓝波段的水体衰减系数小、穿透海水能力强；绿波段和红波段的信息量较蓝光波段更丰富，故对于浅水珊瑚礁地貌有较好的区分度，如图5所示。

图5 不同波段中珊瑚礁类型影像特征比较Fig.5 Comparison of image features of coral reef types in different bands

如3.2节所述，潟湖坡与礁坪之间模糊的坡折线会导致利用GVF-Snake模型提取的该处分界线具有较大的不确定性[45]，这是本文不判读礁坪与潟湖坡分界线的原因。因此，根据蓝、绿和红波段的特点，本文进行所需3条地貌带（外海与礁前斜坡、礁前斜坡与礁坪、潟湖坡与潟湖）瞬时分界线的提取，步骤如图6所示。具体地，外海与珊瑚礁分界线处水深较深，而蓝波段穿透能力最强，提供的水下信息最丰富最全面。因此首先对预处理结果中的蓝光波段进行均值滤波运算，以便消除孤立噪声点，增大外海与礁前斜坡的差异，需要说明的是，在滤波窗口尺度的选择上，我们对比了研究区域的 3×3窗口、5×5窗口、7×7窗口、9×9窗口、15×15窗口和 27×27窗口的滤波结果，发现7×7窗口滤波效果较好，即在平滑礁前斜坡边缘的同时可以更多地保留影像上的细节信息，增大了外海与礁前斜坡在影像上的差异。因此我们使用的滤波窗口大小是7×7。然后利用GVF-Snake模型进行外海与礁前斜坡分界线的提取。礁前斜坡位于水下，而礁坪位于潮间带，这些区域能清晰地显示在绿波段内，且相对于蓝波段，礁前斜坡与礁坪的对比更为明显。因此选取绿波段，利用GVF-Snake模型进行礁前斜坡与礁坪分界线的提取。潟湖坡以珊瑚砂为主，潟湖靠近潟湖坡处亦分布有珊瑚砂，使得两者在影像中的对比度会降低，进而导致获得的分界线位置的不确定性增大。相比于蓝波段和绿波段，透水性较差的红波段能有效避免潟湖珊瑚砂的干扰，即潟湖坡与潟湖反射率的对比在红波段更为明显。因此选取红波段，利用GVF-Snake模型进行潟湖坡与潟湖分界线的提取。

图6 不同时相的地貌带分界线提取示意图Fig.6 Diagram of geomorphic zone boundaries extracted in different time

为了验证利用以上方法进行分界线提取的稳定性，本文先将2019年12月26日外海与礁前斜坡的瞬时分界线向外海扩张1个像素以便作为基准线；然后通过美国地质调查局推出的ArcGIS软件插件—DSAS[46]生成围绕羚羊礁的垂直于基线的517条断面线；进而计算不同时相提取获得的地貌带分界线之间沿着断面线的距离偏差绝对值的平均值。结果表明，Sentinel-2各时相MSI影像提取获得的外海与礁前斜坡分界线之间的平均偏差在6.23～9.55 m的范围内，礁前斜坡与礁坪分界线之间的平均偏差在3.47～5.42 m的范围内，潟湖坡与潟湖分界线之间的平均偏差在7.79～9.81 m的范围内。各地貌分界线的平均偏差均小于1个像素，故可以认为利用GVF-Snake模型提取获得的珊瑚礁地貌带分界线具有比较好的稳定性。

3.4 珊瑚礁面积的估算

珊瑚礁若未经历人为施工或大型自然灾害，其地貌带的分布不会在短时间内发生显著变化。但卫星过境时的潮位差异会导致基于单时相影像所提取的地貌带分界线位置存在较大的不确定性。理想情况下，通过单时相影像也仅能获得珊瑚礁地貌带局部区域的高潮/低潮瞬时分界线。理论上，若采用遥感影像的数量足够多，则会使获得珊瑚礁的高潮/低潮瞬时分界线的概率增加。根据图4所示的珊瑚礁各地貌带的定义，外海与礁前斜坡分界线水深大于一般潮时的分界线水深，其特点就是比一般潮时的分界线更深入外海；同理，礁前斜坡与礁坪分界线在理论最低潮位附近，比一般潮时的分界线更深入外海；潟湖坡与潟湖分界线的特点就是比一般潮时的分界线更深入潟湖。以此为依据，为了改善珊瑚礁地貌面积提取的稳定性和可靠性，对各时相遥感影像提取的珊瑚礁地貌带瞬时分界线进行融合（图7），以便进行珊瑚礁面积的估算。具体方法如下：将提取的瞬时分界线转为面，对礁前斜坡和礁坪−潟湖坡的多时相面取并集，对潟湖的多时相面取交集。

图7 珊瑚礁面积估算示意图Fig.7 Schematic diagram of coral reef area estimation

据此得到各类珊瑚礁地貌面积分别为

式中，S1为外海与礁前斜坡分界线内的面积；S2为礁前斜坡与礁坪分界线内的面积；S3为潟湖的面积，S1、S2和S3如图8所示。

图8 S1、S2和 S3示意图Fig.8 Schematic diagrams of S1, S2 and S3

需要明确的是，本文所述“面积”指的是将珊瑚礁垂直投影到WGS84参考椭球相应高斯投影面上所形成的面积，而并未对倾斜的地形进行校正，倾斜地形（例如礁前斜坡、脊槽系统）的实际面积要比垂直投影得到的面积大。

3.5 效果评估

3.5.1 面积与边界线提取的精度验证

本文通过如下3个方面对提取面积的精度进行评估：

（1）通过将提取的地貌分带与实测水深进行比较，来确定本文所提取的珊瑚礁地貌带分界线与实际地貌带分界线之间的一致性。具体地，将实测水深由外海至潟湖绘制成珊瑚礁水深剖面图，进而计算得到相应的坡度剖面。因为不同地貌带之间的坡度变化较大，故珊瑚礁不同地貌带的分界线，即坡折线位于坡度剖面的极值点。同时，30 m水深处往往也是外海与礁前斜坡分界线的坡折线。据此，将提取的各地貌带分界点与剖面中坡度极值点、30 m水深点的位置进行比较，即可验证本文提取获得的珊瑚礁地貌带分界线与实际分界线之间的一致性。

（2）采用人机交互方式绘制出WorldView-2影像外海和礁前斜坡的分界线、礁前斜坡与礁坪的分界线以及潟湖和潟湖坡的分界线。其中外海和礁前斜坡的分界线与实测剖面中30 m水深处平均位置的误差为1.6 m，远小于Sentinel-2和Landsat 8影像1个像素的大小。这能够确保通过目视解译WorldView-2影像估算的珊瑚礁面积具有较高的精度，因而World-View-2影像获取的面积可以作为Sentinel-2影像面积估算结果精度评估的重要参考值。需要说明的是，珊瑚礁最外轮廓线处坡度很大（接近垂直），向海延伸，水深下降很快，往往是礁前斜坡坡折线所在之处，故最外轮廓线与30 m等深线的距离偏差很小（事实上，本文后面的结果也确实验证了这一点，如图9所示）。此外，由于WorldView-2影像分辨率远高于Sentinel-2影像，我们认为WorldView-2影像目视解译的精度相对于Sentinel-2影像和Landsat 8影像可以忽略不计。因此，本文将由WorldView-2影像获得的各地貌类型的面积作为参考面积，用于与多时相Sentinel-2 MSI影像获得的面积进行比较分析。

（3）将岸线精度分析方法（附录）引入到本文珊瑚礁地貌带分界线提取精度的分析中，即通过完整度、正确度和提取质量等3个指标定量地分析本文地貌带分界线提取的精度。

3.5.2 面积提取的稳定性分析

按照如下两种方式分析珊瑚礁面积提取的稳定性：

（1）Landsat 8多时相影像和Sentinel-2多时相影像（表A1）提取获得的珊瑚礁面积分别与由World-View-2影像获得的参考面积进行比较，进而探讨单时相与多时相珊瑚礁面积估算结果之间的差异。

（2）Sentinel-2影像和 Landsat 8影像（表 A1）所提取的外海与礁前斜坡的瞬时分界线，分别按1条，2条，……，15条共15种组合方式随机选择分界线，每种组合方式15组，各组分别按照本文方法（图7）获取相应的珊瑚礁面积估计值，统计每种组合方式下获取得到的15个珊瑚礁面积估计值的标准差，进而探讨多时相遥感影像珊瑚礁面积估算结果的波动大小，并说明不同分辨率遥感影像获取结果之间的差异。

4 实验结果与讨论

4.1 面积与边界线的提取精度

坡度极值点通常在各地貌带间坡折线的位置，可作为地貌带分界线[33]，本文将6条实测水深剖面、坡度剖面与相应位置的地貌带分界线位置进行比较，结果如图9所示。可以看到，本文提取的外海与礁前斜坡分界线、礁前斜坡与礁坪分界线以及潟湖坡与潟湖分界线的位置均能处在坡度变化极值附近，与坡度剖面所反映的地貌带分界线位置相吻合。同时可见，本文多时相影像提取的羚羊礁最外轮廓线位置，即外海与礁前斜坡分界点位置恰好在30 m水深点附近。因此，所提取的地貌类型与实际地貌类型具有良好的一致性。

图9 水深剖面与地貌分带的对比Fig.9 Comparison of water depth profile and geomorphic zonation

利用多时相提取获得的地貌带分界点与30 m水深点、坡度极值点的位置差异如表1所示。提取的最外轮廓线处水深与30 m水深点的位置差异在5.7～9.5 m范围内，差异不超过1个像素；提取的外海与礁前斜坡、礁前斜坡与礁坪、潟湖坡与潟湖的3类分界点与坡度极值点的距离范围分别为0.1～4.9 m、0.2～4.8 m和0.5～4.8 m，差异最大约为0.5个像素。相比之下，随机挑选10景Sentinel-2影像分别提取3种地貌带分界线各10条，计算最外轮廓线与30 m水深点的位置差异平均为22.9 m，差异超过2个像素；影像提取的地貌带分界点与坡度极值点两者距离平均为18.7 m，差异亦接近2个像素。

表1 提取地貌带分界点与实测剖面30 m水深点、坡度极值点的距离对比（单位：m）Table 1 The distance between extracted boundary points of geomorphic zone with 30 m water depth points and measured slope extreme points is compared (unit: m)

提取得到的地貌带分界线与相应实测数据点位置偏差越小，其面积估算精度越高。总结以上结果可知，提取得到的最外轮廓线的精度由高到低依次为WorldView-2、多时相Sentinel-2和单时相Sentinel-2；提取礁前斜坡与礁坪分界线以及潟湖坡与潟湖分界线的精度由高到低依次为多时相Sentinel-2、单时相Sentinel-2。因此，相比于传统的单时相来说，本文基于多时相影像提取的地貌带分界线与实际地貌带分界线的吻合度更高，而这恰好能够确保珊瑚礁地貌面积地估算能够具有较高的精度。

根据3.4节所述方法得到羚羊礁总面积为17.22 km2，其中礁前斜坡、礁坪−潟湖坡和潟湖的面积分别为1.76 km2、10.29 km2、5.17 km2。计算单时相、多时相的提取面积与参考面积的差异，其结果如图10所示。可以看到，多时相影像的面积提取差异显著低于单时相。由多时相影像得到的珊瑚礁面积与参考面积的差异仅为0.02%，而由单时相得到的珊瑚礁面积与参考面积的差异在0.28%～1.14%之间；由多时相影像得到的礁坪面积与参考面积的差异仅为0.09%，而由单时相得到的礁坪面积与参考面积的差异在0.31%～1.22%之间；由多时相影像得到的潟湖面积与参考面积的差异仅为0.23%，而由单时相得到的潟湖面积与参考面积的差异在0.87%～1.24%之间。换而言之，相对于单时相，本文提出的方法能够显著地提高珊瑚礁面积估算的精度。

图10 不同时相的提取面积与参考面积的差异Fig.10 The difference between extraction area and reference area in different phases

由图11可见，本文多时相融合的边界线的完整度、正确度和提取质量都显著高于单时相。对于外海与礁前斜坡分界线而言，多时相影像的边界线提取完整度为84%、正确度为83%、提取质量精度为72%。单时相提取完整度最高为79%、提取质量精度最高为70%，其中除2018年7月4日外，单时相影像的提取质量精度均低于70%，提取的正确度均小于80%。对于礁前斜坡与礁坪分界线而言，多时相影像的边界线提取完整度为82%、正确度为85%、提取质量精度为81%。单时相提取完整度最高为78%、提取质量精度最高为74%、提取的正确度最高为83%、对于潟湖坡与潟湖分界线而言，多时相影像的边界线提取完整度为80%、正确度为77%、提取质量精度为79%。单时相提取完整度最高为78%，提取质量精度最高为76%、提取的正确度最高为75%，这表明，本文基于多时相遥感影像的珊瑚礁边界线提取效果要优于单时相遥感影像的方式。

图11 不同时相的边界线提取精度Fig.11 Different phase boundary extraction accuracy

4.2 面积提取的稳定性和应用潜力分析

如图12所示，通过多时相影像提取地貌的分界线与从高分辨率影像中观测到的实际地貌分界线位置吻合程度高于单时相，且多时相不同分辨率影像提取的结果比单时相影像提取的结果具有更好的一致性。

图12 单时相和多时相不同影像提取结果对比Fig.12 Comparison of different extraction results of single phase and multi-temporal phases

基于Landsat 8影像（分辨率为30 m）提取获得的羚羊礁面积为17.29 km2，与从Sentinel-2影像得到的数值（17.22 km2）差异不大，说明不同分辨率影像所提取结果之间的良好一致性。需要指出的是，卫星传感器的观测角对面积估算结果会产生一定的影响，但本文采用的Sentinel-2 MSI和Landsat 8 陆地成像仪（Operational Land Imager, OLI）传感器的观测天顶距都接近于0°，影像幅宽不太大，所以可以认为Sentinel-2和Landsat 8影像都接近于垂直观测影像，能够满足垂直投影面积计算的需要。再结合两种影像所提取结果的一致性，我们认为本文中传感器观测角对珊瑚礁面积计算结果的影响可以忽略不计。不同影像数量和不同分辨率影像提取面积的标准差如图13所示。就影像数量（时相）而言，Sentinel-2和Landsat 9多景影像（即多时相，影像数量大于1）所得的面积标准差均小于单景影像（即单时相，影像数量为1）所得的面积标准差，说明本文多时相影像所提取的面积具有良好的稳定性。而且，随着所包含的时相数量越多，提取得到的面积标准差越小—利用6景以上影像所得的珊瑚礁面积标准差均稳定为较小的数值：Sentinel-2和Landsat 8影像提取的羚羊礁面积标准差分别不超过0.01 km2和0.05 km2，与高分辨率WorldView-2影像目视解译得到的面积差异分别不超过0.2%和0.5%。而影像数量增加至10景以上后，Sentinel-2和Landsat 8影像提取的羚羊礁面积标准差趋于一致，且更加稳定。这表明面积估算所采用的时相影像数量越多，提取得到的面积波动就越小。

图13 Sentinel-2和Landsat 8不同数量影像提取面积的标准差Fig.13 Standard deviation of extraction area of Sentinel-2 and Landsat 8 images in different quantities

前面的结果与分析已表明，本文所提的方案能够减小因潮位变化等因素而导致的珊瑚礁估算的不确定性，亦能够减小不同遥感数据源对估算珊瑚礁面积的影响。这意味着，参照本文所提方案，不同时相不同分辨率的多源遥感影像在珊瑚地貌面积估算方面具有很大的应用潜力：（1）能够在较短的时间跨度内以较高的精度实现珊瑚礁面积的估算；（2）能够用于估算不同时期珊瑚礁地貌的面积，实现长时间珊瑚礁地貌的精密监测。

4.3 与同类数据的比较

本文利用所提方法基于多时相Sentinel-2影像获得羚羊礁面积为17.22 km2，与高分辨率WorldView-2影像目视解译获得的17.22 km2相同。将提取得到的羚羊礁地貌分布与联合国环境规划署的世界保护监测中心[16]和艾伦珊瑚礁[17]地图数据集进行比较，结果如图14所示。

由图14可以看到：就珊瑚礁地貌分类来说，本文将珊瑚礁地貌分为3个类别，包含遥感影像中可见的羚羊礁的所有地貌。艾伦珊瑚礁地图数据集划分了9个珊瑚礁地貌类别，但其大多数地貌带由数个碎小图斑构成，分散且不连续；地貌存在误分现象（如礁顶部分）；潟湖东北部和礁体东南出露沙洲处存在孔洞，这一部分未纳入珊瑚礁面积测算范围，其给出的羚羊礁面积为16.60 km2，与高分辨率WorldView-2影像和本文多时相Sentinel-2影像的结果相差0.62 km2。世界保护监测中心地图数据集只包含低潮出露的礁体这一地貌带，潟湖和礁前斜坡部分亦未纳入珊瑚礁面积的估算范围，其给出的羚羊礁面积为11.48 km2，与高分辨率WorldView-2影像和本文多时相Sentinel-2影像的结果相差5.74 km2。就珊瑚礁地貌带分界线来说，利用本文方法提取的珊瑚礁位置与影像之间能够较好地吻合，线条较平滑。而艾伦珊瑚礁地图数据集边界线不够平滑，有很明显的“马赛克”形态，且珊瑚礁最外轮廓线和礁坪的分界线均向礁体一侧过度收敛。世界保护监测中心地图数据集则明显存在珊瑚礁轮廓与影像不相符的问题，其对于潟湖坡与潟湖的分界线拟合较粗糙。

图14 不同数据集羚羊礁地貌带提取结果对比Fig.14 Comparison of Lingyang Reef geomorphic zone extracted by different dataset

因此，目前艾伦珊瑚礁地图数据集和世界保护监测中心地图数据集均低估了珊瑚礁的面积。故结合4.1节所述的结果，可以认为本文基于多时相遥感影像的珊瑚礁面积估算的精度要高于艾伦珊瑚礁地图数据集和世界保护监测中心地图数据集。

5 结论

本文针对南海珊瑚礁面积至今仍然缺乏共识和可靠估算方法的问题，兼顾珊瑚礁和遥感影像两者的特点，从珊瑚礁地貌单元的定义开始，提出了一种低成本利用多时相多光谱遥感影像进行珊瑚礁面积的半自动化估算方法。具体地，首先利用GVF-Snake模型自动精细化提取珊瑚礁地貌带分界线，然后再对精细化提取得到的边界线转化为面要素进行多时相融合，从而实现珊瑚礁面积的估算。利用自2015年12月至2021年4月共计53景Sentinel-2 MSI遥感影像按本文方法得到的羚羊礁的总面积为17.22 km2，其中礁前斜坡、礁坪−潟湖坡和潟湖的面积分别为1.76 km2、10.29 km2和 5.17 km2。实验与分析表明，本文在珊瑚礁面积估算方面，能够用低成本的10 m分辨率Sentinel-2 MSI影像和30 m分辨率Landsat 8 OLI影像获得接近1.8 m分辨率WorldView-2影像的面积估算精度，且估算结果稳定可靠，具体如下：

（1）该方法具有较高的珊瑚礁面积估算精度：获得的珊瑚礁最外轮廓线与30 m等深线的位置差异在5.7～9.5 m范围内（不超过1个像素），与实测水深所指示分界点的位置偏差能控制在0.2～4.9 m的范围内（不超过0.5个像素），以上数值均小于基于单时相结果相应差异。同时，估算的面积与高分辨率World-View-2影像目视解译得到的面积的差异仅为0.02%，优于由单时相得到的面积差异（0.28%～1.14%），亦优于同类的艾伦珊瑚礁地图数据集（3.63%）和世界保护监测中心地图数据集（33.32%）。且该方法获得的珊瑚礁边界线的完整度、正确度、提取质量能够由单时相平均的60%、64%和54%分别提高至84%、83%和72%。

（2）该方法能够降低潮位的影响，减小基于不同遥感数据源的珊瑚礁面积估算结果之间的差异，具有良好的可靠性：基于Landsat 8影像（分辨率为30 m）提取获得的羚羊礁面积（17.29 km2）与从Sentinel-2影像（分辨率为10 m）得到的数值（17.22 km2）差异不大。另外，与高分辨率WorldView-2影像（分辨率为1.8 m）目视解译得到的面积相比，6景以上多时相Sentinel-2 MSI和Landsat 8 OLI影像提取的珊瑚礁面积的最大差异从常规单时相的1.3%和1.9%分别降至0.2%和0.5%；而最大标准差从单时相的0.08 km2和0.20 km2分别降至0.01 km2和0.05 km2。

这些结果表明，本文提出的方法在珊瑚礁面积估算方面具有可行性和有效性。随着我国卫星技术的高速发展，发射的高分辨率卫星将会越来越多，这为大范围珊瑚礁面积估算提供了良好的条件。未来有望结合本文方法，利用国内外免费的高分辨率多光谱卫星影像，实现低成本、稳定、精确地估算南海乃至全球珊瑚礁面积。

附录 地貌边界线提取精度的分析方法

参照周亚男等[47]和乔学瑾等[48]所述的海岸线精度分析方法，以参考线为基线，分别向线内、外建立10 m的缓冲区，记落入缓冲区中的自动提取的分界线总长度为TP1（True-Positive），反之记其总长度为FP（False-Positive）；以自动提取的分界线为基线，分别向线内、外建立10 m的缓冲区，记落入缓冲区中的参考线总长度为TP2（True-Positive），反之记其总长度为FN（False-Negative），如附图 A1 所示。

图 A1 精度评价示意图[2]Fig.A1 Diagram of accuracy evaluation

定义以下参数：

式中，参数Complete表示提取结果完整度（值越接近1表示结果越完整）；Correct表示提取结果正确度（值越接近1表示结果越正确）；Quality通过正确度与完整度综合评价分界线提取质量（值越接近1表示结果质量越好）。

表 A1 本文所用遥感影像信息列表Table A1 List of satellite images used in this paper

续表 A1

续表 A1