岸外沙洲潮间带地形的增强型遥感构建方法

周永，张东,2*，邓慧丽，徐南，张慧铭，郝昕，沈永明

( 1. 南京师范大学 海洋科学与工程学院，江苏 南京 210023；2. 江苏省地理信息资源开发与利用协同创新中心，江苏 南京210023)

1 引言

潮滩位于海陆交互区域，具有独特的地形[1]、沉积物[2]和水文特征[3]。作为众多生物的栖息场所，潮滩在保持生物多样性、维持生态系统稳定等方面具有极高的生态价值[4]。随着沿海经济的不断发展，潮间带的经济价值日益凸显[5]。但人类开发活动如滩涂围垦、围海养殖、港口建设等在带来大量土地和经济效应的同时，也给沿海地区尤其是潮滩地区造成了巨大影响[6]。近年来，随着自然影响和人类开发活动力度的加大，潮滩环境变化加快，冲淤变化显著，潮滩资源正在面临极大的危机。针对潮滩的动态变化监测与利用保护刻不容缓[7]。

江苏岸外辐射沙洲作为江苏淤泥质潮滩的重要部分，在资源开发利用和海岸演变研究等方面具有重要价值[8]。但受较远的离岸距离、恶劣的勘测环境和周期性淹没下滩面持续变动等因素影响，岸外辐射沙洲潮间带区域长期缺少全面、准确的地形监测数据。基于LiDAR[9]、InSAR[10]技术的潮滩地形获取方法因其在数据获取性和经济性上的不足，制约着在该区域的大面积监测应用。遥感水边线法利用多时相卫星遥感影像数据构建潮滩数字高程模型（DEM），具有数据获取方便、简单实用、精度较高等优势，成为潮间带地形获取的重要手段[8,11]。目前该方法已运用于德国瓦登海[12]、韩国戈姆索和咸平郡湾[13]等多个区域的潮滩地形反演，国内Wang等[14]、Kang等[15]和郑宗生[16]也利用遥感水边线方法完成了长江口、江苏沿海等淤泥质潮滩区域的地形构建，并取得较好结果。

辐射沙洲区域水动力条件复杂，日益频繁的人类活动加剧了水动力环境的调整。虽然在较短周期内，相似潮位下的沙体轮廓范围大体一致，但沙体表面形态变化（如潮沟摆动）会引起水边线的位置发生变化，从而导致在较窄的平面条带范围内集中的大量水边线存在相互交叉现象[17]，造成生成的DEM产生地形破碎化和异常起伏。因此，选择合适的水边线是构建高精度潮间带DEM的基础。Sagar等[18]提出了以长时间序列影像为基础的潮滩高程信息提取方法，通过将序列影像的时间域转化为潮高域来构建合成图层，在获取平均地形变化特征的同时减少影像异常值的影响。该方法在一定程度上可以借鉴用于解决交叉水边线的问题，但其应用于澳大利亚大陆尺度的潮滩地形构建中，只关注了单一潮位站点的潮位序时变化以及长周期下大量像元灰度序列数据的统计特征，而忽略了区域内的潮差空间分布差异以及长时间周期内可能的潮滩冲淤变化，因而难以直接应用到辐射沙洲等潮滩地形变化较快的区域。此外，潮间带沙体的坡度相对平缓，高差不大，滩面分布有大量的树枝状和贯通状潮沟，然而现有的DEM构建方法只关注沙体整体的形态，却忽略了潮沟的微地貌表达[19]，这也直接影响了潮间带DEM的构建精度和呈现效果。

因此，本文拟针对常规水边线方法在应用于地形多变的岸外沙洲潮间带DEM遥感构建中出现的水边线交叉误差和缺乏考虑潮沟微地形的问题，研究增强型潮滩DEM构建方法，实现岸外沙体DEM的快速构建，为利用卫星遥感数据监控潮间带地形变化提供技术支撑。

2 研究区与数据

2.1 研究区概况

辐射沙洲位于江苏中部沿海，北起射阳河口，南至蒿枝港，由70多条沙脊和潮流通道组成，沙体表面平坦开阔，在潮汐作用下被周期性淹没。中部条子泥-蒋家沙-竹根沙区域作为辐射沙洲的核心区，紫菜养殖[20]、风电开发[21]等人类开发活动频繁，沙体短周期内变化监测需求较大，因而被选作研究区域，进行潮间带DEM构建试验。研究区范围如图1所示，位于32°38′～32°59′ N，120°58′～121°25′ E之间，南北跨度约37 km，东西跨度超过40 km。潮汐类型以正规半日潮为主，浅海分潮显著。潮波系统受东海前进潮波与南黄海旋转潮波共同影响，二者在弶港岸外辐合，波能集中，区域内平均潮差为3.9 m，最大潮差可达9.28 m，且潮差以弶港为中心分别向南、北逐渐减小[22]。

2.2 遥感影像数据

综合可覆盖研究区范围、清晰少云等要求，选择了25景10 m空间分辨率的Sentinel 2 MSI影像和8景15 m分辨率的Landsat 8 OLI影像作为遥感数据源，进行潮间带DEM构建，影像的成像时间分布如图2所示，时间跨度为2017年10月至2019年6月。所有影像经过大气校正、几何精校正、图像增强等预处理，几何精校正精度优于1个像元。由于空间分辨率不同，影像均重采样至30 m。

图2 影像成像时间分布Fig. 2 Image acquisition time

2.3 潮位数据

由于地球曲率变化、浅海潮波变形等影响，辐射沙洲岸段潮位变化呈现出的“跳绳效应”[22]，条子泥-蒋家沙-竹根沙海域的潮差变化较大，无法用单一站点的潮位代表水边线上不同位置的潮位空间差异。因此，在研究区北、西、南3个方向上分别选择了大丰港、弶港和洋口港3个潮位站点，进行潮位数据收集，站点位置如图1所示。根据潮汐调和分析原理，利用站点的短周期（2个月以上）逐时潮位资料，即可计算其潮汐调和常数，实现站点在影像成像时刻的潮位模拟并具有较好精度[23]。因此，收集了上述3个站点2018年6月至2019年1月间逐小时的潮汐表潮位过程数据（http://global-tide.nmdis.org.cn/），利用T_TIDE潮汐调和分析工具进行分潮模拟[24]，根据信噪比大小确定出显著分潮，模拟得到研究区的潮位误差小于20 cm。

2.4 DEM验证数据

收集了2018年蒋家沙-竹根沙局部海域的水下声呐实测地形数据，空间分辨率为250 m，高程基准为1985国家高程基准。在蒋家沙沙体上按照经纬度方向、在竹根沙沙体上按照平行和垂直于沙脊走向方向划定4条验证剖面，分别命名为JJS-H、JJS-Z、ZGS-H和ZGS-Z，位置如图1所示。通过提取断面上的实测点高程，验证遥感方法构建的潮间带DEM精度，评价DEM质量。

图1 研究区概况Fig. 1 Sketch map of the study area

3 研究方法

常规遥感水边线方法是基于潮间带区域在潮汐作用下被周期性淹没的特性，将从遥感影像上提取并赋予潮位信息的瞬时水边线视为“等高线”，利用空间插值处理，得到潮间带DEM。本文在遥感水边线法的基础上，提出面向岸外沙洲潮间带地形的增强型遥感构建方法，实现以下两方面的改进：（1）根据较短周期、较小潮差内淹没概率相似的原理，引入影像序列合成方法，提取合成影像的海陆边界线，解决常规水边线方法中对于较窄区域内多时相水边线交叉重叠造成的DEM高程异常和破碎化难题；（2）提出潮沟DEM构建方法，通过在常规水边线方法得到的DEM上镶嵌潮沟DEM，实现对潮滩表面微地形的表达，提高潮间带DEM精度，改善可视化效果。具体技术流程如图3所示，主要技术方法阐述如下。

图3 基于遥感水边线的增强型地形构建方法技术流程图Fig. 3 Flowchart of enhanced terrain construction method based on waterlines

3.1 基于空间及潮位关系特征的影像排序

构建高精度的潮间带DEM，关键是用于空间插值的瞬时水边线在平面上具有位置分离性，在高度变化上具有良好趋势性。因此，首先需要对从33景预处理后的遥感影像上提取得到瞬时水边线按照先平面位置、后潮位关系进行排序。

由于潮位滞后效应影响，大丰港、弶港和洋口港3个潮位观测站点的潮位序时变化不同，需要选择合适的站点作为影像序列的潮高排序基准。基于此，分别以其中1个站点的潮高为基础，对全部33景影像按高低排序，然后对3个站点的潮位变化情况进行线性拟合，通过对比趋势线的拟合决定系数R2和斜率，选择潮位排序的基准站点。图4a至图4c分别显示了以大丰港、弶港和洋口港潮位站为基准进行影像成像时刻潮位排序后的潮位分布和线性拟合结果，表1列出了对应的拟合决定系数R2和斜率值。

表1 不同排序基准下各站点的潮位变化拟合决定系数与斜率Table 1 The coefficient of determination and slope for fitting tide level changes at different tide stations under different sequencing conditions

相较而言，如图4c所示，以洋口港站为排序基准下，大丰港站和弶港站潮位变化趋势线性拟合效果相对较好，大丰港站、弶港站、洋口港站的拟合决定系数R2分别为0.74、0.32和0.93，潮位点变化离散起伏较小，表明以洋口港站为排序基准，大丰港与弶港站的潮位变化趋势更具有规律性。进一步对比趋势线斜率，以洋口港站为排序基准时，大丰港站、弶港站、洋口港站潮位趋势斜率分别为-0.11、-0.06和-0.13，3个站点的潮位变化斜率变化差异最小，潮位变化速率最为接近。因此选定洋口港站作为基准站点进行潮位排序，其中，12条具有正确空间与潮位关系特征的水边线，后续直接参与DEM插值构建；其余21条水边线根据潮高大小相似、平面位置相近、线段相交的原则，分成5组，各组潮差范围在0.13～0.4 m之间，平均潮差间隔0.24 m。5组水边线对应的影像组如表2所示，用于下一步的影像序列合成及海陆边界线提取。

表2 各影像组中的影像与成像时刻潮位高度Table 2 Images and their correspondent tide levels for each image group

图4 各潮位站点在不同排序条件下影像成像时刻的潮位分布和线性拟合Fig. 4 Tide level distribution and linear fitting of imaging time at various tide stations under different sequencing conditions

3.2 影像序列合成与海陆边界线提取

影像序列合成方法通过统计一个特定周期内影像各像元的灰度值域分布，来指示该像元位置被潮水淹没的频率状况，然后结合潮高变动范围，将时间序列影像的灰度值域转变为潮高域，在特定潮高范围内获得预设淹没频率下的海陆边界线[18]。以其代替分布于较窄区域内的多时相水边线，可去除水边线交叉所带来的高度异常起伏。

将用于序列合成的遥感影像按组裁剪为相同大小，采用改进的归一化差值水体指数（MNDWI）进行图像增强处理[25]，来加大影像中的海陆差异，并进行影像序列合成。MNDWI公式如下：

式中，Rg为 绿波段反射率，在Landsat 8 OLI和Sentinel 2 MSI中同为B3波段；Rswir为短波红外波段反射率，在Landsat 8 OLI和Sentinel 2 MSI中分别对应B6波段和B11波段。在淤泥质潮滩，由于短波红外波段对浑浊海水的光谱反射率较低，因此使用MNDWI处理后，海陆边界色调对比度更强，有利于海陆边界线的识别。

与以长时间序列影像为基础，取MNDWI中位数作为合成影像像元值的方式不同[18]，由于各影像组内原始影像的时间跨度短、数量少，经过筛选和排序后各像元MNDWI序列中基本不存在离群值，利用中位数去除异常值影响的作用并不显著，因此选取序列的平均值来完成新图像的合成。这样通过计算像元的MNDWI平均值作为合成图像的像元值，将较短周期内各影像中呈现的沙体变化进行均值化处理，可以在弱化单一成像时刻下瞬时沙体形态的同时，获得代表地形平均变化特征的海陆边界位置信息。

图5 显示了以影像组HC-4为例的5景遥感影像的MNDWI序列合成与海陆边界线提取过程。可以看到，取序列MNDWI均值形成的合成图像（图5c）中出露沙体与海水的边界总体清晰，仅在部分区域如蒋家沙东部和西南部呈现海陆灰度值相似、边界不明显状况。通过观察该合成图像的MNDWI分布直方图，取0.42的谷值作为分割阈值，生成海陆二值图像（图5d），进一步利用Sobel边缘检测算子提取出海陆边界线。以边界线为基础、300 m为距离构建缓冲区，叠加合成前5景原始影像提取的瞬时水边线，除部分区域如蒋家沙西北侧与陆岸相连区域以及东南侧潮沟摆动频繁区域以外，各瞬时水边线均位于该缓冲区内，表明序列合成影像提取的海陆边界线结果能够较准确的还原该潮位区间内沙体的平均轮廓信息。图5e显示了5组合成影像提取的海陆边界线结果。

图5 影像组HC-4的MNDWI序列合成与海陆边界线提取过程Fig. 5 MNDWI sequence synthesis and sea-land boundary extraction of image group HC-4

3.3 初始潮间带DEM构建

由于提取的水边线或合成的海陆边界线不能被简单视为等高线，因此将提取出的12条具有正确空间与潮位关系特征的水边线和从合成影像组中提取的5条合成海陆边界线，按照30 m的间隔进行空间离散，针对其环状的形态特征，选用反距离平方加权方法进行空间潮位内插[26]。其中，12条水边线的潮位插值依据为各潮位控制站点在影像成像时刻的潮位，5条合成海陆边界线的潮位插值依据为各潮位控制站点在合成组内各影像成像时刻潮位的平均值。对所有的离散点通过潮位控制站点的潮位值进行反距离平方加权潮位插值及赋值处理后，利用不规则三角网（TIN）方法，构建得到研究区的初始潮间带DEM。

3.4 具有潮沟微地形特征的潮间带DEM构建

潮间带表面潮沟大多宽度窄、深度浅、高差小，通过空间插值方法直接反演难度较大。因此参考平原河道地形嵌入思路，通过单独构建潮沟区域地形DEM，然后将其镶嵌至初始潮间带DEM上，可以实现对潮间带微地形的细致模拟[27]。

选取覆盖研究区的2018年2月23日Sentinel 2 MSI卫星遥感影像进行潮沟边界线和中线提取，该影像潮滩出露基本为最大范围，潮沟形态清晰完整。首先参照水边线提取流程，提取得到矢量格式的潮沟两侧边界线，用于控制潮沟范围，潮沟宽度从0 m至575.92 m；然后将两侧边界线封闭形成面图层，提取出潮沟中线，用于控制潮沟深度。对潮沟边界线与中线按照30 m间隔离散，通过前述3个潮位站点的潮位值利用反距离平方加权进行潮位空间插值，得到离散点的高程。对于提取的潮沟边界线与中线，3个潮位站点的潮位取值方法如下：潮沟边界线直接取所提取的影像的成像时刻潮位；潮沟中线由于并非直接提取自影像，因此假定所提潮沟的底部在全年最低潮位下均被海水淹没，选择排序基准洋口港站在2018年最低潮位时刻对应的3个潮位站点的潮位，作为其潮位插值基准。然后同样利用不规则三角网（TIN）方法，插值构建出潮沟DEM。利用ArcGIS的图像镶嵌工具，将构建后的潮沟DEM以像元为基础嵌入原始潮间带DEM，完成对潮沟位置处的地形替换和边界平滑处理，即可得到最终具有潮沟微地形特征的潮间带DEM。

4 结果

4.1 DEM构建结果

对比常规遥感水边线方法与增强型地形构建方法构建的沙体地形，常规遥感水边线方法以全部33景影像所提取的瞬时水边线为基础进行空间插值，构建的潮间带地形如图6a所示，高度范围为-2.30～1.69 m；增强型潮滩构建方法构建的DEM如图6b所示，高程范围为-2.68～1.69 m。将二者所构建的DEM作差值运算，得到的差值分布呈现明显正态分布特征，平均差值为-2.93 cm，差值范围在±0.3 m以内的区域占64.64%，可见两者构建的DEM结果整体上保持了相似的地形信息。但是由于增强型地形构建方法进行了潮沟镶嵌处理，将原始潮滩上的部分平缓地形替换为V字型潮沟地形，造成滩面高程的最低值下探，因此高程变化范围略大于前者。

图6 地形构建结果对比Fig. 6 Comparison of terrain construction results

在地形呈现方面，常规遥感水边线方法构建的DEM表面多呈点状或线状分布，形态较为破碎，起伏变化杂乱，且无明显的潮沟形态。高程较高的沙脊区域分布于蒋家沙西部和高泥的西北部，虽然能看出呈块状趋势，但总体较为破碎。而增强型潮滩构建方法减少了因为水边线交叉冗余所造成的高度异常起伏与地形破碎，增强了潮间带滩面的潮沟微地形特征，因此可以看到构建的潮间带DEM表面起伏过渡平滑自然，地貌单元多呈块状分布，潮沟形态清晰，多以线状或树枝状分布于沙体边缘，尤其在蒋家沙南部较为集中。沙脊区域位置与前者一致，但空间形态更加完整。

4.2 与DEM实测地形对比分析

4.2.1地形精度分析

增强型地形构建方法获得的潮间带DEM空间分辨率为30 m，能够较为精细的描述地形起伏细节，但实测剖面是由250 m间隔下的剖面点组成，不同分辨率基准无法直接进行精度对比，因此首先采用升尺度方法，将模拟DEM的空间分辨率插值为250 m；然后采用平均绝对误差（Mean Absolute Error, MAE）、均方根误差（Root Mean Square Error, RMSE）和相关系数（r）3个指标进行DEM高程精度验证。

根据4条剖面的同名点高程数据对比，利用增强型地形构建方法得到的条子泥-蒋家沙-竹根沙DEM平均的MAE为0.43 m，RMSE为0.54 m，地形起伏变化的相关系数r为0.75。从各断面来看，如图7所示，误差最大剖面为JJS-Z，其MAE与RMSE分别为0.53 m和0.65 m；误差最小的为ZGS-Z，MAE为0.37 m，RMSE为0.34 m。地形趋势变化中相关系数最高的断面同样为ZGS-Z，r为0.85，最低为JJS-H，r为0.65。4条剖面中，竹根沙沙体的误差结果较蒋家沙沙体偏小，主要是由于蒋家沙、条子泥区域与陆岸的距离较近，受陆岸围垦等人类开发活动影响较大，水动力环境复杂，沙体东侧的东大港大型潮沟系统向东南方向摆动，潮沟拓宽浚深；西侧的西大港潮沟向南贯通，造成蒋家沙沙体冲淤频繁，滩面小型潮沟变化活跃，对水边线的提取与DEM插值精度产生影响。而竹根沙区域距离陆岸较远，水动力环境相对蒋家沙而言较为稳定，沙体在研究周期内变化较小，因此DEM精度相对较高。此外，增强型地形构建方法得到的地形DEM反映的是模拟时段内的平均地形信息，具有一定的地形概化与平均效果，而实测数据显示的是地表测量时点的瞬时形态，因而DEM获取方法的差异及数据所代表的时效性也是造成两者地形之间存在差异的一个重要原因。

图7 各剖面实测地形与模拟地形的精度和趋势相关性对比Fig. 7 Comparison of the accuracy and trend correlation between the measured and simulated terrain at each section

4.2.2剖面形态分析

图8 显示了增强型地形构建方法获得的沙体DEM与实测DEM的剖面形态对比。可以看出，JJSH、JJS-Z和ZGS-Z 3条剖面中实测与模拟地形的起伏位置总体基本吻合，但在地形起伏程度上，除了ZGS-H剖面的模拟地形起伏较实测地形稍大以外，其余剖面实测地形的起伏更为显著。如图8a所示，JJSH剖面中实测地形在距离起点3 000 m和9 000 m处有各有一个明显的沙脊隆起，高程分别达到0.3 m和0.8 m；而在模拟地形中，相同位置上同样有小的地形上升，但上升高度并不明显。图8b所示的JJS-Z剖面中，实测地形在距起点2 000 m和7 500 m处都存在明显的潮沟下切，深度可达到约-2.5 m；在模拟剖面中两处位置同样存在潮沟，但下切深度仅达到-2 m左右，较实测地形存在一定差异。这主要是由于序列合成处理将沙体在不同时相的变化平均化，且潮位赋值也同样以潮高变化的平均值为主，因此相比于实测地形，模拟DEM的高程变化幅度呈现被压缩状态。图8d中ZGS-Z剖面的反演结果在地形起伏方面基本与实测剖面一致，其中在距起点约0～2 000 m的范围内呈现出地形上升，而后2 000～6 500 m呈现出缓慢下降趋势，表明竹根沙中沙体北侧较南侧更陡。

ZGS-H剖面（图8c）与前面3个剖面呈现出相反特征，在距离起点2 000～8 000 m的范围中，实测剖面相对较为平缓，但是模拟地形起伏较为明显。对比影像数据可知，起伏位置都表现为潮沟分布。由于模拟DEM的空间分辨率高，结果能呈现更为细致的起伏结果，且升尺度后的同名位置点仍然能保持周期影像内的起伏信息，但实测地形在250 m的分辨率下，部分较小潮沟信息无法呈现，因而该段区域内的地形呈现出轻微起伏下的平缓状况。

图8 剖面形态对比Fig. 8 Morphology comparison of each section

5 结论

本文提出了一种岸外沙洲潮间带增强型地形遥感构建方法，通过将常规遥感水边线法与序列影像合成法相结合，完成小范围和短周期下的潮滩地形反演；进一步引入潮沟DEM镶嵌技术，实现了表现潮滩微地貌特征的岸外沙洲潮间带DEM快速构建。研究结果表明：

（1）根据模拟的潮间带DEM与实测地形之间的4条验证剖面对比，平均的MAE为0.43 m，RMSE为0.54 m，相关系数为 0.75，两者的地形起伏位置较为吻合。但是由于增强型地形构建方法得到的是模拟周期内的平均地形变化状况，因此地形整体起伏偏小。

（2）对比增强型地形遥感构建方法与常规遥感水边线法得到的潮滩DEM结果，两者的平均差值为-2.93 cm，差值呈现明显的正态分布特征且范围在±0.3 m以内的区域占64.64%，可见增强型地形遥感构建方法能够较好地保留潮滩平均高度变化信息。同时该方法减少了因滩面潮沟摆动导致水边线交叉造成的高程误差和地形破碎，构建的潮间带DEM表面起伏过渡自然，潮沟形态清晰，因此适用于淤泥质潮间带地形的快速构建，并保证较高的DEM模拟精度。

该方法在合成影像组的确定过程中，需要首先对所有影像的水边线进行提取，而后在水边线空间排序的基础上再对比潮位高度的变化范围，从而确定影像组对应的影像数据。但是从以上影像组提取的水边线并未参与后续的潮间带DEM构建，所以后续研究需要对影像中沙体出露范围判断和影像成像时刻潮位匹配流程进行改进优化，来提高本潮间带增强型地形遥感构建方法的建模效率。

致谢：感谢欧洲航天局（ESA）和美国地质调查局（USGS）提供的免费卫星影像资料，以及国家海洋信息中心提供的潮位资料信息。