基于Sentinel-2数据的山口红树林生理参数反演及时空特征

王诗文， 何宏昌， 付波霖， 孙宏亮

(桂林理工大学测绘地理信息学院, 桂林 541004)

红树林生长在热带及亚热带的海岸潮间带，属于特殊的湿地沼泽类型之一。它具有很独特的生态特性，在防浪护堤、保持生物多样性和净化环境污染等作用[1]，同时在保障沿海安全和生物多样性等方面有着举足轻重的位置。近几个世纪，湿地受到了不同程度的影响，导致退化速度加快，包括海平面上升、人类对资源的不合理开发以及气候的变化[2-3]。中国红树林减少了50%以上，尽管恢复了部分地区对红树林的生态状况，但是对红树林的保护和修复刻不容缓[4]。目前，有效地开展红树林生态健康评估，可为科学保护与管理红树林湿地资源提供技术支持，对于维持生态系统平衡和可持续发展有着重要的意义。而生理参数则是评价红树林的生态系统状况的重要基础，准确地获取生理参数的信息成为关键[5]。

山口红树林湿地是位于陆地生态系统与海洋生态系统交错的过渡地带，是一个动态且复杂的生态区域，对于传统的野外调查存在一定的局限性，难以满足实时监测需求。遥感技术具有覆盖面积大、数据更新周期短、空间分辨率高等特点，已成为中外红树林监测的主要技术之一[6-7]。遥感为红树林的测绘和评价提供了一种准确、高效、重复性强的方法。由于红树林(周期性海水淹没的潮间滩涂)特殊的生长环境和茂密的森林，人们很难进入红树林进行广泛的野外调查和采样。因此，在过去20年里，遥感数据被广泛用于评估红树林[8]，例如绘制全球红树林分布图[9]，监测国家[10]或区域[11]红树林的范围和动态，确定当地红树林的物种组成[12]，估计生物物理指标[13]。在绘制红树林范围时，中分辨率多光谱图像是最常用的数据。Chen等[14]采用了大量30 m陆地卫星7/8图像与Sentinel-1A图像相结合绘制中国红树林地图。Leon等[15]利用多时相10 m SPOT影像揭示了红树林土地覆盖变化和破碎的时空动态。然而，这些中分辨率多光谱数据所提供的空间和光谱信息可能不足以详细研究红树林及其物种组成。例如，海岸线上总是存在狭窄范围的红树林(<30 m)或孤立的红树林群落(<900 m2)。为了区分红树林物种的类型，最常用的数据是商业高分辨率多光谱图像，其次是高光谱图像。Wang等[16]比较了IKONOS和QuickBird数据在巴拿马加勒比海海岸红树林物种分类方面的表现。虽然高分辨率图像的使用已得到相当大的关注，但其购买成本高、条带面积小的局限性妨碍了连续或大规模的地理覆盖监测。因此，大规模红树林范围和物种组成的测绘可能取决于新兴传感器的能力。“哨兵-2A”是两颗完全相同的卫星中的第一颗，于2015年6月23日发射升空；第二颗“哨兵-2B”于2017年3月7日发射升空。Sentinel-2提供13个多光谱波段，其中4个传统波段(红/绿/蓝/近红外)的空间分辨率为10 m，3个波段是为植被探测设计的红边波段[17]。Sentinel-2数据已用于植被[18]，如估算北方森林冠层覆盖率和叶面积指数[19]，但很少用于红树林。由于WorldView-2(如数据复杂性、高成本和小条带大小)和Landsat 8图像(如缺乏红边带、中等空间和时间分辨率)，Shapiro等[20]和Pastor等[21]建议分别使用Sentinel-2图像进行红树林范围监测和红树林叶绿素浓度估算。Castillo等[13]评估了Sentinel-1和Sentinel-2检索上述红树林类群生物量及其替代土地利用的能力，结果表明，红边波段的表现优于可见光波段和短波红外波段。Valderrama等[22]比较了Landsat 8、SPOT-5、Sentinel-2和WorldView-2对3种红树林物种进行了遥感分类，2018、10、1468-3，共27组，但仅使用归一化差异植被指数(normalized difference vegetation index,NDVI)的阈值来区分物种。研究表明基于含有红边范围内光谱数据对反演生理参数有着较好的相关性。还有学者如谢军飞等[23]利用MODIS产品分析2010—2012年北京植被光合有效辐射吸收率空间分布特征，得出光合有效辐射吸收率(fraction of absorbed photosynthetically active radiation,FAPAR)空间分布呈现东西部高，并向中心逐渐递减的分布特征。吴国训等[24]基于MODIS数据反演江西省2000—2011年植被叶面积指数，江西省植被叶面积指数 (leaf area index,LAI)呈现明显季节变化。可见对红树林长时间序列生理参数的分析是对红树林生态评估提供参考。

根据红树林中植被的生物特性，由李珊珊等[25]提出红树林特征指数BSTDEV对估算叶面积指数有着一定的优越性；Zhou等[26]利用SPOT6卫星和UAV数据资料估算三沙湾滨海湿地植被覆盖度且具有较高精度；Francisco Flores-de-Santiago等[27]根据测定叶绿素浓度、叶面积指数以及叶长来对红树林进行生态评估，得出健康的红树林叶绿素含量也没有明显的季节差异，但冠层上部的叶绿素a含量较高。综上所述，构建生理参数反演模型可以较好地反映区域整体特征，但是，目前研究集中于某一领域中某一种或两种生理参数反演鲜有对红树林4种生理结构在长时间序列和不同生长期动态变化的研究。

以广西山口红树林自然保护区为例，以含有红边波段的Sentinel-2影像数据为主要数据源对广西山口红树林生理参数利用SNAP软件对逐月的叶绿素浓度(chlorophyll-A/B，CAB)、叶面积指数[28](LAI)、植被覆盖度[29-30](fractional vegetation cover,FVC)以及光合有效辐射吸收率[31](fraction of absorbed photosynthetically active radiation,FAPAR)4种生理参数进行反演，结合MODIS 影像数据产品采用BP神经网络算法构建回归模型来评价Sentinel-2影像反演生理参数的数据精度；利用4种生理参数分别从年际和不同生长期两个角度分析红树林植被生长变化，对山口红树林植被变化监测提供依据，为湿地长时间序列多生理参数特征分析提供参考。

1 研究区域及数据源

1.1 研究区概况

山口红树林国家级自然保护区地处广西壮族自治区北海市合浦县的沙田半岛，东西两侧分别于丹兜海内和英罗湾西侧左下角，英罗湾东侧为广东省界[32]，地理位置为21°29′00″～21°33′25″N，109°42′03″～109°45′36″ E，如图1所示。研究区总面积8 000 hm2(1 hm2=104m2)，光热条件好，属于南亚热带季风型海洋性气候带，年平均气温22.9 ℃[33]。广西北部湾红树林内有海水呈现片状连续或间断分布，有着很强的分布特征，有着与其他地物明显不同的表现，红树林的植物资源丰富，有维管束植物106科370多种，其中有中国极为罕见连片的红海榄纯林和高大的通直的木榄[34]，具有重要的科考价值且在降低海洋灾难、净化环境等方面发挥了重大作用[35]。随着沿海经济开发，船舶运输业、渔业及农业等人为活动影响，大量污染汇聚于红树林湿地，导致湿地严重退化，对山口红树林湿地研究刻不容缓。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 数据源获取与预处理

选取1—12月影像数据，11个月的影像数据进行研究，其中7月成像均有遮挡无法使用。研究数据包括Sentinel-2A/B影像数据、MODIS MOD13A2和MOD15A2数据产品。Sentinel-2属于欧洲空间局(ESA)卫星[36]，其中包括Sentinel-2A和Sentinel-2B两颗卫星，两颗卫星互补，重访周期为5 d，高度为786 km，可覆盖13个光谱波段，幅宽达290 km。地面分辨率分别为10、20、60 m，从可见光和近红外到短波红外，具有不同的空间分辨率，在光学数据中，Sentinel-2数据是唯一一个在红边范围[37]含有3个波段的数据，这对监测植被健康信息非常有效。

MODIS中分辨率成像光谱仪搭载于Terra、Aqua两颗卫星，是美国对地观测系统(EOS)计划中重要传感器。选择MODIS 全球1 km分辨率植被指数16 d合成产品和全球1 km叶面积指数/FAPAR 8天合成产品(表1)。

表1 Sentinel-2A/B多光谱影像和MODIS数据产品

数据预处理利用ArcGIS软件对Sentinel-2影像几何重采样，将Sentinel-2数据与MODIS数据进行地理配准、辐射定标以及对地形校正等，通过欧洲航天局(ESA)[38]提供的Sen2cor模型[39](Python语言编写)实现，处理后波段数由13个减少为10个，剔除了气溶胶、水汽和卷云波段的影响。MOD13 A2和MOD15 A2产品已经完成了几何纠正。

2 研究方法

2.1 红树林生理参数遥感定量反演

Sentinel-2 L2A级影像数据反演生理参数方法是利用SNAP 5.0软件[40]下的Biophysical Processor(生物物理处理模块)，采用人工神经网络算法(artifificial neural network，ANN)。ANN算法主要包括生成一个全面的植被特征和相关的天顶冠(TOC)反射率的合成数据集。然后训练神经网络从TOC反射率和定义观测结构的相应角度估计冠层特征。对于每个生物物理变量，都要对一个特定的神经网络进行校准，其中PROSAIL模拟光谱反射率按Sentinel-2波段结构(B2、B3、B4、B5、B6、B7、B8A)，然后用来训练网络。使用反向传播ANN，它由一个输入层、一个或多个隐藏层和一个输出层组成。其中输入层由11个标准化输入数据组成，包括B3、B4、B5、B6、B7、B8a、B11、B12、cos(viewing_zenith)、cos(sun_zenith)、cos(relative_azimuth_angle)；隐层中有5个具有传递功能的神经元；带有线性传递函数为输出层。通过这种方法分别反演出LAI(叶面积指数)、CAB(叶绿素含量)、FVC(植被覆盖度)以及FAPAR(光合有效辐射吸收率),再调用mosaicing(镶嵌)和subset(取子区)模块将6景影像的生物物理指标进行拼接、裁剪、异常值处理，获取空间分辨率为 10 m 的山口红树林自然保护区的生理参数。最后将LAI、CAB、FVC以及FAPAR导入Arcgis10.2软件中进行分析。将反演后的分辨率为10 m的4种结构参数分别重采样，与MODIS影像数据分辨率一致。

MODIS植被数据产品中MOD15A2产品通过Arcgis10.2软件提取出LAI和FAPAR,MOD13A2产品计算的FVC公式如式(1)所示；NDVI(植被指数)与CAB呈现正相关[41],将反演出来的CAB进行归一化处理，再与NDVI进行相关性分析。

(1)

式(1)中：NDVIsoil=NDVImin,NDVIveg=NDVImax；NDVI 为混合像元的植被指数；NDVIveg为纯植被像元最大值；NDVIsoil为纯土壤像元最小值。

2.2 基于BP神经网络的Sentinel-2生理参数结果和MODIS数据产品回归模型的构建

BP(back propagation)神经网络[42]是1986年Rinehart和McClelland为首的科学家小组提出，按误差逆向传播算法训练的多层前馈网络。它的学习规律则是使用最速下降法，通过反向传播来调整权值和阈值，使误差平方和最小。同样BP神经网络存在一个输入层、隐含层和输出层，每层之间互相连接, BP神经网络结构如图2所示。

图2 BP神经网络结构Fig.2 Structure of BP neural network

采用BP神经网络由 MATLAB实现Sentinel-2A/B影像数据反演的生理参数与相对应的MODIS数据产品做一元线性回归和相关性分析。网络共有3层结构，包括输入层、输出层和隐含层。将获取的1—12月的数据分成11组数据；按照不同生长期划分为4组数据。每组数据按照2∶1的比例分成测试集和训练集，训练集的数据输入到模型中用于学习，测试集的数据对于构建模型输出的结果进行验证。通过选取网络的输入层、输出层及隐含节点，利用训练数据进行训练来建立适合山口红树林的验证模型，并利用测试集对构建好的模型进行验证，来评价验证模型的精确度。输入层分别为逐月的Sentinel-2A/B影像数据，输出层分别为相对应逐月的MODIS数据产品。由于输入层神经元为1，因此隐含层个数为3。构成了3层的BP神经网络模型结构。经过反复训练已确定隐含层神经元个数为3，训练设置迭代次数为10 000,训练目标为0.000 1，模型模拟得较好。

3 结果与讨论

3.1 基于Sentinel-2A/B生理结构反演结果与MODIS估算的数据产品对比

基于BP神经网络对MODIS产品影像数据与Sentinel-2A/B影像数据线性回归分析，选取年内11个月4种生理参数得出相关系数R在95%置信区间内均高于0.86，效果较好。由表2中4种生理参数的年平均均方根误差(RMSE)、相关系数R和决定系数R2精度评价指标可以看出，LAI 3种参数值较高，R为0.92，R2为0.85，RMSE为0.19，年内LAI反演效果好，由图3看出，样点向趋势线靠近，并且样点分布紧凑，再由图4(a)可以看出，逐月的LAI的R、R2及RMSE个别月份较低，逐渐趋于稳定，从而证明用Sentinel-2反演LAI可靠，效果好。由表2看出，CAB年际生理参数系数R为0.86，R2较低为0.66，RMSE为0.18，由图3可以看出CAB样点分布较散，多数样点向趋势线靠近，由图4(b)可以看出，逐月的CAB参数系数在5月浮动最大下降到最低，然后逐渐升高至平稳。表2中FAPAR年际生理参数系数平稳，效果较好，R、R2及RMSE分别为0.87、0.76及0.14，由图3看出，年内少量样点较分散，大部分样点向趋势线靠近，呈现较好的相关性，由图4(c)可以看出，逐月的参数系数2月最低逐渐上升趋于平稳；年内FVC由表2看出，R、R2及RMSE分别为0.93、0.84、0.09，由图3看出，FVC年内样点分布靠近于趋势线，样点较分散不够紧凑，由图4(d)看出，逐月的R与R2后两个月系数接近吻合，说明反演效果很好，在4月呈现较低值后逐渐上升趋于平稳。综合来看，逐月和年内的4种生理参数用Sentinel-2数据反演较好，比较可靠。

图3 年际生理参数回归模型Fig.3 Regression model of interannual physiological parameters

表2 年际生理参数模型系数Table 2 Model coefficients for interannual physiological parameters

图4 年内红树林4种生理参数精度指标变化Fig.4 Change of precision index of four physiological parameters of mangrove during the year

3.2 年内红树林生理参数动态变化及空间分布结构特征

山口红树林年内生理参数的动态变化反映了湿地质量的水平对湿地的保护有着重要意义。红树林3个区域分别为核心区、缓冲区和实验区，如图5所示。核心区是重点保护区域，区域内人为活动较少；实验区以从事参考观察、科研试验和积极活动等；缓冲区以保护为目的科研、生态修复工程以及生态旅游等[43]。根据Sentinel-2影像数据利用SNAP软件反演了逐月的生理参数。将反演后的3个区域的生理参数分别导入ArcGIS软件进行区域统计分析。图6和表3展现了研究区内3个区域年内动态变化以及具体数值，整体来看，3个区域1—12月中生理参数发生明显变化的月份在1—2月和9—12月。图5选取了4种生理参数值最低月份和最高月份，由图5可见3个区域中核心区反演效果较好。由图6看出，4种生理参数相对稳定，核心区CAB在2月呈现最低值为16.71，最高值出现在10月为65.08，最高值与最低值之间相差48.37；FVC在2月呈现最低值为0.23，最高值在9月为0.38，最高值与最低值之间相差0.15；FAPAR年内核心区整体很稳定，最低值在2月为0.25，最高值在9月和12月均为0.38，其中最高值与最低值相差0.13；LAI最低值出现在8月为0.67，最高值出现在9月为1.50，最高值和最低值之间相差0.83。实验区CAB最低值在2月为10.58，最高值则在9月为43.67，最高值与最低值之间相差33.09；FVC最低值也出现在2月为0.18，最高值在6—9月平稳于0.29，最高值与最低值相差0.11；FAPAR在12月出现最低值0.04，最高值10月为0.27，最高值与最低值相差0.23；LAI最低值在1月为0.60，最高值在9月为1.04，最高值与最低值之间相差0.44。缓冲区，CAB最低值出现在2月为1.07，最高值出现在10月为17.29，最高值与最低值之间相差16.22；FVC最低值也出现在2月为0.04，最高值出现在6—8月均为0.14，最高值与最低值之间相差0.10；FAPAR在2月出现最低值为-0.03，最高值在12月为0.31，最高值与最低值之间相差0.34；LAI在1月出现最低值为0.24，最高值则在2月为0.65，最高值与最低值之间相差0.41。

图5 红树林自然保护区生理参数空间分布Fig.5 Spatial distribution of physiological parameters in Mangrove Nature Reserves

表3 研究区三个区域植被生理参数Table 3 Physiological parameters of vegetation in three regions of the study

图6 红树林3个区域生理参数年内变化Fig.6 Changes in physiological parameters in three regions of Mangrove forests

综上所述，由图5总体来看，4种生理参数最高值均出现在核心区，最低值均出现在缓冲区，主要原因可能是实验区人员活动较频繁核心区人员活动较少，对湿地破坏较少。而缓冲区人为活动较多对湿地破坏较为严重。CAB和FVC年内3个区域相对稳定，其中个别月份出现骤降，骤降的原因与自然和人为方面密切相关，例如海平面上升、气候突然变化以及人为乱砍滥伐等因素[44]，这与祝萍等[45]在中国典型自然保护区生境的时空变化中人类扰动对植被覆盖度变化相关吻合。FAPAR年内核心区和缓冲区稳定上升的趋势，而实验区变化较大，极有可能是人为原因导致的。年内LAI也变化浮动较大，3个区域变化浮动都很大，环境因素应该是主要原因，可能受全球变暖的影响，气候的变化导致LAI变化较大，植物的自身机理会为适应环境而产生变化[47-48]。

由表4和图7可知,由于植被还没有生长起来，在萌芽期(1—3月)植被覆盖度还呈现最低值，FVC由0.19开始逐渐升高至0.30再下降至0.20处于稳定；在生长期(7—8月)处于夏季，夏季温度高，植被生长较快，植被覆盖度处于最高值，在成熟期气温逐渐降低降水量也减少，落叶植被会受到影响植被覆盖度降低；FAPAR由萌芽期0.13上升至最高值0.37后逐渐平稳于0.23，FAPAR呈现先增大后减少的趋势，与植物的生长规律基本吻合[49]。LAI在萌芽期处于最低值0.31，逐渐升高到0.47走向稳定；发育期(4—6月)植被开始展叶，LAI逐渐增加，在1—3月萌芽期植被绿色器官刚开始生长，LAI呈现最低值,这与夏传福等[48]在基于MODIS叶面指数的遥感物候分析基本吻合在年季的时间尺度上，基本反映了植被的生长轨迹；成熟期(9—12月)植被开始走向逐渐升高至稳定的过程，由于气候的原因FAPAR和FVC呈现的值较低，而CAB和LAI呈现出的值较高。

图7 不同生长期四种生理参数平均值Fig.7 Average values of four physiological parameters at different growth stages

表4 不同生长期Sentinel-2A/B反演植被结构参数

综上，不同生长期红树林的植被生理参数则变化不同，由此可见不同生长期对4种生理参数影响不同，FVC、LAI和FAPAR在萌芽期和发育期呈上趋势，CAB和LAI在成熟期(9—12月)呈上升趋势。

4 结论

红树林具有很独特的生态特性，在定量研究山口红树林植被生理参数时序变化特征成为研究热点之一。在获取研究区的Sentinel-2和MODIS影像数据基础上，对逐月的Sentinel-2影像数据定量反演了时序的山口红树林4种生理参数CAB、FVC、FAPAR和LAI，并利用BP神经网络算法与MODIS数据进行回归分析，验证了Sentinel-2 数据反演生理参数可靠性较高，并且分析了逐月以及不同生长期红树林植被的动态变化以及生理参数的差异。

以MODIS数据产品为基础，利用BP神经网络模型，基于Sentinel-2影像数据反演4种生理参数包括CAB、FVC、LAI和FAPAR，从时间序列和空间上掌握广西山口红树林年内整体植被动态变化，对红树林4种生理参数反演、对比以及分析得到以下结论：Sentinel-2影像数据反演4种生理参数与MODIS数据产品趋势基本相同，呈现较好的相关性，因此，利用Sentinel-2A/B多光谱影像反演红树林生理参数具有较高的精度；山口红树林年内的生理参数存在波动变化的趋势，在2月LAI、CAB、FAPAR及FVC 4种生理参数均出现最低值，分别为0.30、0.08、0.08、0.13。最高值呈现在不同月份，分别为10月、12月、8月、4月，存在某个月份环境的突然变化或砍伐等造成指标降低等可能，还需进一步研究；从不同生长期来看整体反演效果好，其中CAB最低值在生长期0.20，最高值在成熟期0.25；FAPAR 在萌芽期最低0.13，最高在发育期0.37；LAI最低值同是在萌芽期0.31，最高值在成熟期0.47；FVC最低值在萌芽期0.19，最高值在生长期0.31；因此得出成熟期植物反演效果最好。研究结果对红树林植被变化监测提供依据，为生态评估及保护提供参考。