海岸带淤积与侵蚀变化, 是学者及社会关注的热点问题。红树林湿地, 作为世界上四大高生产力海洋生态系统之一, 分布在热带、亚热带低能海岸潮间带上部, 在全球生态平衡中起着不可替代的重要作用。近几十年来全球气温急剧上升, 导致最明显和最直接的反应就是局部区域海平面上升。生存在平均潮位线与平均大潮高潮位线之间的红树林, 对海平面上升极为敏感^[1]。国内外学者们对未来海平面上升影响下红树林分布变化趋势进行了相关研究^[2-4]。结果表明, 红树林系统稳定性很大程度上取决于相对海平面上升和林区内部沉积速率之间的相对关系^{[3, 5-7]}。但由于人工堤坝以及农业用地和养殖池的建立束缚了红树林的向陆发展。而红树林向海边缘的广阔淤泥质光滩, 是一种有利于红树林和盐沼发展的特殊地貌, 其冲淤状态即高程变化, 同样对红树林的现在与未来的生存环境起着非常重要的作用^[8]。

对于红树林区域外侧的广阔光滩, 由于潮汐周期性的淹没与干出, 及沉积物颗粒极细、含水量较大等特殊性, 使得该区实地地形测量极为困难。遥感技术在海岸带信息提取方面具有宏观、快速、同步、高频度动态观测和节省资金等突出优势。遥感结合GIS技术在用于构建历史时期潮滩地形^[9-12], 以及分析不同时期岸线变迁^[13-14]、潮滩冲淤变化^[15-16]等方面的研究, 已经广泛应用于河口及开阔海域滩涂。而针对于红树林这一特殊生物海岸地貌下的光滩冲淤演变研究涉及较少, 因此本文以具有广阔红树林和潮滩分布的狭长型海湾-英罗湾为研究区, 基于多时相遥感影像提取水边线, 结合FVCOM模型通过迭代模拟得到水边线的潮位高程, 利用GIS方法构建2008年和2015年两个时期的光滩数字高程模型(Digital Elevation Model, DEM), 通过叠置分析对红树林区域外光滩的冲淤变化情况进行研究。

1 研究区概况

我国红树林主要集中在南海西北海域, 包括广西、广东及海南省。本文研究区英罗湾位于北部湾东北部, 雷州半岛西海岸北端的浅海湾, 是两个国家级红树林保护区的结合部, 即省界以西为广西省山口国家级红树林保护区的英罗港林区, 以及省界以东的湛江国家级红树林保护区高桥核心区。湾内分布面积约占80%以上的有砂泥质光滩, 滩上分布纵横交错的小潮沟, 可为红树林和其他动植物的生长发育提供良好的条件, 同时, 湾内向陆侧筑有人工护岸堤, 再向内则为上万亩的人工养殖区, 其地貌如图 1所示^[13]。

2 数据与方法 2.1 数据介绍

本文以2008年、2015年为代表, 共收集不同潮位共计34景多传感器光学遥感影像, 其中包括6景Landsat卫星系列影像数据(http://glovis.usgs.gov/), 26景HJ-1A/1B卫星CCD影像数据(http://www.cresda.com/CN/), 1景GF-1号卫星WFV影像数据, 1景CBRES-04卫星CCD影像数据, 如表 1所示。为了防止波浪对模拟水位的影响, 数据选取避开波浪影响较大的月份(6、7、8月)。由于潮水淹没红树林区域时, 从遥感影像上提取的是红树林的边界, 并非水边线, 所以只挑选含有光滩出露的中低水位时刻影像。

为了三维漫滩模型FVCOM模拟水位的需要, 收集了研究区中国人民解放军海军司令部航海保证部2015年发布的1︰40 000安铺港及附近电子海图、1︰32 000铁山港电子海图以及2013年发布的江洪港至北海港1︰120 000纸质海图^[18], 作为模型中的初始地形。模型开边界处的水位驱动数据来源OTIS (OSU Tidal Data Inversion)(http://volkov.oce.orst.edu/tides)中China Sea范围的水位预测值。

在研究区内, 为了验证水位模拟结果, 收集了研究区附近的铁山港潮汐表数据。在2015年3月利用中海达厘米级RTK-GPS, 在研究区潮滩进行考察, 对部分有植被覆盖和陆地边界进行高程测量, 并且测量了6个断面用于反演地形的验证实验。其中, 研究区内潮汐数据和海图深度基准一致, 而实测高程为85国家高程基准。根据广东省航道局所提供的英罗湾内龙头沙港建港所测量的高程基准换算关系, 所有高程数据均统一为当地平均海平面。

2.2 研究方法

本文以上述遥感影像数据为基础, 时间上分为2008年和2015年两个时间段。通过遥感影像解译技术, 提取水边线并建立水边点集。利用FVCOM模型对相应时刻进行水位模拟, 并将最邻近网格水位值赋给水边点。为了避免未知光滩地形对水位模拟的影响, 通过迭代的方法, 利用空间插值生成的DEM更新FVCOM模型中光滩地形后, 重新驱动模型。比较迭代前后DEM之间的差值, 当差值在阈值范围内, 输出最后生成的DEM为最终光滩的反演地形^[19]。在此基础上, 利用空间叠置的方法, 分析2008年和2015年两个时期光滩的冲淤变化, 其流程如图 2。

2.2.1 水动力模型

本文中模拟水位的水动力模型采用由美国马萨诸塞大学海洋科技研究院和伍兹霍尔海洋研究所联合开发的FVCOM(An Unstructured Grid, Finite-Volume Coastal Ocean Model)无结构网络有限体积近海数值模式。模式采用有限体积积分法, 与有限差分法和有限元模式中使用的微分形式不同, FVCOM将控制方程的积分形式离散化。这种积分方法结合了有限元方法和有限差分方法的优势, 一方面计算效率得到提高, 另一方面其在水平空间上具有很好的几何灵活性。该方法体现了有限元方法能够较好地拟合自然岸线的优点, 在模拟河口、海岸等固体边界复杂区域的应用上效果较好, 且能更好的保证动量、质量、热量、盐度守恒^[20]。模型中底部拖曳系数(C_d)通过拟合底边界层内流速的对数分布来获得, 利用距底部z_ab处高度算的:

$ {C_{\rm{d}}} = \max \left( {\frac{{{\kappa ^2}}}{{{{\left( {\ln \left( {{z_{ab}}/{z_0}} \right)} \right)}^2}}}, 0.002\;5} \right), $

(1)

其中, κ=0.4, 为冯卡曼系数, z₀为海底粗糙度, z_ab是距水底的水柱高度^[16]。底部粗糙度在模型中利用水深进行计算。为了能更好地模拟英罗湾内的潮汐, 模拟边界设在英罗湾外, 利用SMS(Surface Water Modeling System)软件设置地形网格分辨率由研究区内的30 m到海洋开边界的2 000 m。网络边界在研究区以及地形变化剧烈的区域加密, 计算区域共分成15 240个节点, 28 974个三角网, 地形格网如图 3所示。模式的水深结合海岸地形内插求得, 考虑到研究区大面积潮滩, 开启干湿网格处理模块。本次计算的开边界采用了O₁、K₁、M₂、P₁、Q₁、S₂、N₂, 8个分潮的调和常数, 利用TMD工具箱导出开边界的水位值驱动模型。模型中水中的拖曳系数按照式(1)计算。红树林和其根部对水流的拖曳系数较大, C_d在红树林边缘设置为0.047, 在林区内部设置为0.033^[21]。模型采用冷启动, 假设初始时刻海水静止, 即初始水位和初始流速均为零, 温度和盐度初值取自研究区平均值, 未考虑开边界温盐通量对模型的影响。注入英罗湾内河流径流流域面积非常小, 具有溪流性质, 在枯水期多数河流接近断流。由于模拟期避开受波浪影响季节以及台风、热带气旋天气, 未考虑河流和风的作用。外模步长设置为0.1 s, 内模时间步长为1 s。模拟时间覆盖影像过境时间, 计算过程中每10分钟输出一次结果。

2.2.2 提取水边线

如图 2蓝色虚线框内步骤所示, 在ENVI软件下, 对所有影像进行辐射定标和FLAASH大气校正等操作。为了解决不同时相异源的影像几何配准问题, 结合实际地形图提取环绕潮滩分布的22个人工池塘顶点作为地面控制点, 保证控制点均匀分布, 且对所有影像进行几何精校正, 总体均方根误差不超过0.5个像素。

首先文中所处理的影像, 以908专项岸线数据作为处理的边界制作陆地掩膜, 排除陆地水体区域(如鱼塘、居民地)等干扰信号。通过利用归一化水体指数NDWI^[22], 利用大津法(OTSU)^[23-24]选择最佳阈值, 进行二值化处理, 提取光滩与水体的边界。由于研究区存在很多狭长潮沟, 涨落潮时潮沟内的水流和漫滩的水流不同步, 并且会在短期内摆动, 在影像上容易形成混合像元很难自动提取, 所以小于0.5个像元的狭长小潮沟不计入模型中。最后通过目视解译对数字化岸线进行部分调整, 确保误差不超过一个像素。

3 光滩DEM构建与验证 3.1 光滩DEM构建

利用ArcGIS对提取的水边线按照30米的间隔离散化, 并用临近网格的模拟水位值对水边点进行高程赋值。在坡度比较大的地方(例如红树林带前缘陡坎), 受影像分辨率影响多个相近水位的水边线可能由同一个像元提取。它的特点是某些特定区域, 同一像元内提取的水边点具有两个以上不同的高程值。针对这种异常情况, 只保留具有最高高程值的水边点。将优化好的高程点集, 建立不规则三角网TIN, 再通过双线性插值方法转为规则格网, 插值生成研究区光滩DEM。首次生成的光滩DEM用来修正FVCOM地形网格中的光滩部分, 作为用来修正潮汐模型的光滩地形, 重新驱动模型来降低沿岸水位模拟误差, 生成精度更高的DEM。当第N次DEM与第N-1次DEM之间的差值的绝对值满足精度要求停止迭代, 并将第N次DEM作为最终输出结果(流程见图 2)。

通过分析2008年和2015迭代模拟前后DEM差值绝对值概率分布, 得出2008年模拟的第一次和第二次DEM差值绝对值在0~0.15 m区间内为85%, 第二次和第三次为96%; 2015年模拟的第一次和第二次DEM差值绝对值在0~0.15 m内为86%, 第二次和第三次为99%;两个年份均显示三次迭代后模拟结果趋向稳定, 退出迭代, 将第三次模拟生成的两个时期DEM作为最终输出(图 4)。结果表明, 2015年反演的光滩地形更具有斑块化, 即潮滩上的大中型潮沟进一步发育。通过对比2008年和2015年红树林空间分布, 2015年向海方向存在三个比较明显的新扩张林区, 即图中A、B、C区域, 都是2008年在高程0.25 m以上光滩区域发育而成的。

3.2 精度验证

本文利用研究区最邻近铁山港(位置见图 1)潮汐表对进行模拟潮位验证, 利用2015年实测断面数据进行DEM结果验证。通过比较2015年和2008年的水动力模拟中前18天(432小时)的连续水位, 其结果如图 5所示。误差最大的基本上分布在高潮附近(误差 > 0.2 m), 而所选影像光滩可提取水边线位于潮间带中下部, 即水位不超过1 m。综上, 除去高水位数据(> 1 m), 90%以上水位误差在0.15 m以内, 平均误差为0.08 m, 均方根误差为0.11 m。该数值模式能较好的反映研究区真实的潮位变化。

根据研究区龙头沙码头建港时所用的各高程垂直转换关系, 将实测数据(85高程基准)转换成和模拟水位一致的当地平均海平面基准。通过ArcGIS中“值提取至点”方法, 提取实测高程点所在位置的模拟高程值, 由于所实测高程点不能保证处于像元中心位置, 所提取的像元值是根据相邻像元的有效值进行双线性插值所得到。通过6个剖面实测高程点(位置见图 1)和所对应的反演高程值进行差值运算, 来评价DEM的反演精度。如图 6所示, 随着迭代次数的增加, 误差统计频率直方图的范围逐渐向中央值0 m聚拢, 其均方根误差分别为0.19 m、0.13 m、0.09 m, 即其垂直精度随着迭代运算精度提高。该垂直精度虽仅代表局部断面精度(即实地测量的空间布设主要是依赖于人在潮滩上的可通行性), 并且所验证剖面地形在空间上的连续变化使得整体误差比较集中(即所测的地形由岸向海变化较为平缓, 未涉及细小潮沟的测量), 但其误差分析中数值的变化也反映了该迭代的方法提高反演地形的精度。

通过自然资源部连续对研究区附近海域的海平面变化监测(海平面公报: http://www.mnr.gov.cn/sj/sjfw/hy/gbgg/zghpmgb/), 广西和广东沿海在2008年和2015年两个年份的相对海平面分别相差8 mm和0 mm。所以本文中2008年的光滩模拟中的所涉及的平均海平面高度和基准同2015年一致。

4 冲淤变化结果与讨论

以英罗湾口门为界, 根据生成的最终DEM进行空间叠置, 得出2008~2015时间段内光滩的冲淤变化, 按照以下五类进行分层设色:较强侵蚀区(< –0.4 m)、侵蚀区(–0.4~–0.15 m)、稳定区(–0.15~ 0.15 m)、淤积区(0.15~0.4 m)、较强淤积区(> 0.4 m), 如图 7所示。根据栅格统计分析可得, 由ArcGIS中的填挖方功能进行统计可知研究区光滩的冲淤平均值为–0.03 m, 最大值为0.8 m, 最小值为–1.04 m。红树林岸外光滩处于冲刷的面积约为9.51 km², 侵蚀量约为186.36 × 10⁴m³; 淤积的面积约为7.81 km², 淤积量约为169.85×10⁴ m³。本研究区内2008~2015年光滩稳定性在面积上总体表现为:稳定 > 侵蚀 > 淤积。空间分布整体上表现为:光滩上部淤积, 中部稳定, 下部侵蚀; 西岸淤积量大于东岸。为了进一步阐明其冲淤特征, 根据地貌将研究区光滩划分为三个区域: Ⅰ区是从潮汐通道弯曲后突变狭窄的位置至湾顶; Ⅲ区是从位于西南角的小海湾光滩至口门; Ⅱ区则是连接湾顶Ⅰ区和口门Ⅲ区的中部区域。

在英罗湾湾顶Ⅰ区, 大部分处于稳定-淤积和少量侵蚀的状态。这主要是由于潮汐通道东向偏转使得湾顶处于波浪遮蔽区, 并且波浪传入过程中受到大面积光滩造成的能量耗散。淤积区主要集中在红树林边缘, 其中在西岸有成片的较高淤积区(对应图 4中虚线框A), 这和新红树林的定植有很大相关性。

英罗湾的中部区域Ⅱ和口门附近区域Ⅲ的光滩多数处于稳定状态, 并且侵蚀面积大于淤积面积。侵蚀主要发生在光滩下部, 即潮汐通道的两侧, 其可能受潮汐不对称性, 以及夏季西南方向传入海岸的波浪引起的泥沙搬运的影响。光滩中部较为稳定, 其高程变化在0.15~–0.15 m, 伴随一些沟壑状侵蚀, 这与大中型潮沟的发育和侧向摆动有关。光滩的淤积主要分布于光滩的上部, 由图 7所示, 湾内西岸呈带状, 东岸呈斑块, 且淤积量西岸高于东岸。这一差异主要有以下两方面原因: (1)自然原因:在红树林湿地系统达到平衡前, 光滩上原始海拔决定红树林先锋植物是否可以定植生长, 新红树林的定植也可以促进沉积物在其边缘区域堆积。通过红树林空间分布变化分析, 2008~2015年期间红树林总体稳定且缓慢增长, 林带宽度大的东岸以向林内小潮沟扩张为主, 而西岸以向海潮滩扩张为主。红树林向海扩张的西岸在光滩上部淤积量要高于东岸。而东岸红树林稳定区域的向海光滩基本呈现绝大部分稳定、少量淤积的状态, 并伴随中型潮沟发育。(2)人为原因:通过实地考察发现, 广东省管辖的东岸, 其水产养殖主要以在红树林靠陆一侧的人工养殖池为主, 未发现占用滩涂进行水产养殖的活动。而位于广西省管辖的西岸, 在光滩中下部分布有成片插有竹竿的水产养殖区。网帘、竹竿的架设对水动力的阻滞会大大加速该区域光滩的沉积。位于口门区域Ⅲ东岸龙门沙港口内侧的一片淤积区, 是由于在2009年该港口新建一东西向延伸向海的码头拦截湾内泥沙, 并且附近有一些建筑垃圾。

5 结论

本文利用Landsat系列、HJ-1A/1B、GF-1、CBRES-04等多源多时相的遥感影像数据, 采用水边线与迭代的水动力模型FVCOM相结合的方法, 生成了2008年和2015年两个代表年份的DEM, 提取南海西北红树林保护区岸外光滩的冲淤空间分布, 并结合红树林林区空间分布的变化, 得出了以下结论:

(1) 随着我国地球观测卫星的快速发展(如中巴资源卫星、环境系列卫星HJ-1A/B、高分系列卫星等), 多卫星和短重访周期的影像采集增加了短时间内(小于1 a)收集一定数量高质量遥感影像的可能性, 为利用水边线方法反演潮滩地形提供了数据基础。通过FVCOM水动力模型迭代, 可以削弱海湾内潮汐通道两侧宽广的潮滩以及红树林对水位模拟的影响, 并通过实测数据验证该方法可以提高反演的精度。在缺乏历史地形实测资料的情况下, 可以为红树林保护区外光滩地形演变的定量化研究提供一定的技术支持。

(2) 通过对近期(2008—2015年)南海西北英罗湾红树林保护区外光滩冲淤分析, 可以看出:湾内除了湾顶光滩处于稳定-淤积的状态, 海湾中部至口门的光滩从上部到下部在空间上呈现淤积-稳定-侵蚀的沉积格局。其中湾内西岸光滩上部出现的成片高淤积区域与红树林的岸外扩张有极大的空间相关性, 红树林林区边界稳定的区域其地形冲淤状况也相对稳定; 中部沟壑状侵蚀为潮沟的发育和侧向摆动所致。而人类活动, 如滩涂水产养殖、港口建设等行为, 也使得海岸演变趋势变得更加复杂。