利用Sentinel-2影像超分辨率重建的红树林冠层氮含量反演

甄佳宁，蒋侠朋，赵德梅，王俊杰，苗菁，邬国锋

1.深圳大学生命与海洋科学学院,深圳 518060;

2.自然资源部大湾区地理环境监测重点实验室,深圳 518060;

3.深圳大学建筑与城市规划学院,深圳 518060

1 引言

氮素是限制红树林生长最关键的营养元素，是参与其光合作用的重要成分（Camino 等，2018；Lemaire等，2008），也是评价红树林净初级生产力和健康水平的重要指示参数之一（Loozen 等，2020）。因此，红树林冠层氮素含量CNC（Canopy Nitrogen Content）的准确估算对红树林健康及其生态系统碳循环研究具有重要意义。当前迅猛发展的遥感技术为快速、无损、大范围地获取植被冠层信息，进而反演植被氮素（Wang 等，2019）、叶绿素（Ali 等，2020）以及叶面积指数（Estévez等，2020）等理化参数奠定了数据基础。

植被氮素含量遥感反演的原理是植被叶片中化学组分分子结构中的氮化学键在一定辐射水平的照射下发生振动或电子跃迁，引起相关波长的光谱发射和吸收差异，从而产生不同的光谱反射率，且光谱吸收强度与叶片氮含量密切相关（刘良云，2014）。因此植被氮素含量的变化一定程度上影响叶片对特定波长的光谱反射特性（陈永喆等，2017；王莉雯和卫亚星，2013）。目前，国内外研究学者已对植物叶片尺度、冠层尺度氮素含量遥感反演进行了大量研究。按照模型的理论基础，植被氮素含量遥感反演的算法主要可以归纳为两类：（1）统计方法，包括光谱指数法和回归分析法。例如，Inoue 等（2012）在中国和日本不同环境条件下获得水稻地面高光谱数据集和CASI-3机载高光谱影像，分别分析了由原始反射率值和一阶导数值构建的归一化植被指数NDSI（Normalized Difference Spectral Index）和比值植被指数RSI（Ratio Spectral Index）与水稻冠层氮素含量的相关性，研究发现，基于一阶导数的RSI（D740，D522）是反演CNC 最准确的变量。Lee 等（2020）对比了线性回归、随机森林和支持向量机回归方法估算玉米冠层氮素浓度结果的精度，结果表明两种机器学习模型均优于线性模型。（2）辐射传输模型方法，包括PROSPECT 模型和4-scale 模型等。例如，Yang等（2015）将PROSPECT模型中的叶绿素吸收系数替换为等效氮吸收系数，建立了N-PROSPECT模型，实现了冬小麦叶片氮素浓度的反演；Li 等（2018）将N-PROSPECT 与SAIL 模型结合，探究了N-PROSAIL 模型在叶片和冠层尺度上小麦CNC估算的潜力。

目前，国内植被氮素含量遥感反演总体上处于起步阶段（陈永喆等，2017），主要利用地面高光谱数据进行叶片或小尺度范围内的农作物（小麦、水稻及玉米等）氮含量反演（李哲等，2016）。但利用机载或星载高光谱数据反演区域尺度森林冠层氮素含量的研究较少，尤其是针对红树林的氮素含量反演还未见报道。另外，绝大多数星载高光谱影像空间分辨率较低（30 m），且单一像元平均叶片氮素含量实测困难，限制了反演精度。与传统多光谱数据相比（Landsat 和MODIS 等），Sentinel-2 数据进一步提高了遥感数据的时间、空间、光谱分辨率，是目前唯一具有3个红边波段的多光谱传数据，在植被理化参数的精细估算上具有巨大优势（Xie 等，2019；Clevers 和Gitelson，2013）。然而，Sentinel-2 影像中对植被健康和变化非常敏感的“红边”波段的空间分辨率只有20 m，极大限制了它在植被监测中的潜力（Zhang等，2019）。

图像超分辨率重建技术使Sentinel-2 数据部分波段的较低空间分辨率（20 m）提高至10 m 提供了可能。本文以广东省湛江市高桥红树林保护区为例，利用三次卷积重采样、Sen2Res 和SupReMe算法重建10 m空间分辨率的Sentinel-2影像，构建红树林CNC反演模型，进而实现CNC反演和制图。研究成果可为区域尺度红树林冠层健康状况的精细监测提供有效的方法借鉴和数据支撑。

2 研究区与数据

2.1 研究区概况

高桥红树林保护区（21°31′N—21°35′N，109°45′E—109°48′E，图1）是广东省湛江市红树林国家级自然保护区的核心区之一，位于广东省最西端，北部湾东北角，拥有中国大陆现存连片面积最大的红树林带（毛丽君，2011）。该地气候类型为南亚热带海洋性季风气候，年平均水温25 ℃—27 ℃，且不随气温的急剧升降而变化，这构成了红树林生长的先决条件之一。4—9月为多雨季节，8月雨量最多，10—3月雨量较少。

图1 研究区地理位置与采样点空间分布Fig.1 Location of study area and sampling points

研究区内红树林覆盖度高，树种呈现明显的分带特征，主要树种包括桐花树（Aegiceras corniculatum）、木榄（Bruguiera gymnorrhiza）、秋茄（Kandelia obvolata）、红海榄（Rhizophora stylosa）、白骨壤（Avicennia marina）、 海漆（Excoecaria agallocha）和无瓣海桑（Sonneratia apetala）等，岸带分布有黄槿（Hibiscus tiliaceus）等半红树植物（朱耀军等，2016），其中有中国连片面积最大的木榄群落，达300 hm2。

2.2 实验数据

2.2.1 野外实验设计及地面数据收集

野外工作于2020年8月6日至12日开展，实验共设立了64 个样方，采集143 个样本点。在研究区内选择成熟的桐花树、白骨壤、木榄、红海榄和秋茄5种代表性群落类型设立固定样地。考虑到Sentinel-2 遥感影像空间分辨率最高为10 m，每种群落设置15 m×15 m 的样地。在每个样地内，选择2—3 株、间距5—10 m 的红树林样木，用手持GPS记录其点位信息；在每株样木处，采集每株样木北、南、西、东4 个方向半球影像图片，采用HemiView 数字植物冠层分析系统提取冠层叶面积指数。利用高枝剪采集每株样木冠层上方10 个成熟叶片，并迅速封装在含有冰块的保鲜袋内带回实验室。在实验室内，采用凯氏定氮法（王学奎，2006）测定叶片全氮浓度。

2.2.2 Sentinel-2影像

Sentinel-2 卫星是欧洲空间局发射的新一代高时空分辨率卫星，包括Sentinel-2 A/B 两颗卫星，均搭载多光谱传感器，提供了由可见光、近红外到短波红外的13 个多尺度、中高空间分辨率的遥感影像（表1）。本研究采用从ESA scihub 网站（https://scihub.copernicus.eu/［2021-07-01］）免费下载的Senitnel-2 卫星影像（2020年9月4日获取）L1C级大气顶层反射率产品，该产品已经过辐射定标和几何校正。利用SNAP 软件的Sen2Cor 插件进行大气校正，将影像大气顶层反射率转换成冠层反射率。为了探讨影像分辨率提升对红树林CNC反演的应用潜力，大气校正时，分别将空间分辨率设置为10 m 和20 m。根据Sen2Cor 算法，在进行10 m 大气校正时，B10 波段被剔除，其他波段保持不变，保留12 个波段；在进行20 m 大气校正时，3个60 m波段（B1，B9和B10）和1个10 m波段（B8）被省略，保留了9个波段。考虑到各波段对植被氮素反演的意义，为了使影像的每个波段具有相同数量的像素，本研究移除B1、B9、B10和B8 共4 个波段，并采用三次卷积法将大气校正后的20 m分辨率的波段重采样到10 m。

表1 Sentinel-2卫星波段信息Table 1 Bands information of Sentinel-2 satellite

3 研究方法

3.1 Sentinle-2影像超分辨率重建

本研究Sen2Res和SupReMe的影像超分辨率重建算法分别由SNAP 平台的第三方插件Sen2Res 和MATLAB 平台实现。利用信息熵和平均梯度定量分析影像重建后的质量情况，由于没有参考高分辨率影像，而重采样后的影像与原始影像保持了良好的光谱一致性（Zhang 等，2019），因此，本研究将Sen2Res和SupReMe两种算法重建后的影像与重采样后的影像进行光谱相关性分析来评估影像光谱一致性。

3.1.1 Sen2Res算法

Sen2Res 插件为Brodu（2017）针对Sentinel-2卫星提出的一种利用相邻像素之间“共享”几何信息的超分辨率多光谱多分辨率影像重建方法。该算法基于一个像素区域内具有相同性质地物的比例（几何信息）是该像素的物理属性，与波段无关，而地物的反射率会随波段而变化；利用高分辨率波段拟合整体混合模型，识别相邻像素间的通用共享信息作为独立于反射率的像素属性，进行几何编码。在保持光谱特征不变的同时，将通用几何信息用于解混低分辨率像素，实现低分辨率波段超分辨率重建。此过程在SNAP 平台的Sen2Res插件中完成，算法具体细节见Brodu（2017）。

3.1.2 SupReMe算法

SupReMe （Super-resolution for multispectral multiresolution estimation）模型针对单一传感器在不同光谱范围内获取不同空间分辨率的多光谱影像，实现低分辨率波段超分辨率到最高分辨率的算法（Lanaras等，2017）。该模型依赖影像模糊和下采样过程的观察模型，充分利用高分辨率波段的纹理信息，通过对高分辨率不连续数据进行编码并传递给低分辨率波段制定规则器。在给定二次数据和二次正则化条件下，给出一个凸优化问题来反演观测模型，并提出一个有效的数值求解方法（Lanaras等，2017）。

3.2 红树林冠层氮素含量估算

3.2.1 红树林冠层氮素含量获得

Sentinel-2 影像获取的是红树林冠层尺度的反射率光谱，决定氮素敏感的光谱特征的通常是其冠层氮素含量。因此，本研究参考Wu 等（2010）的方法，将红树林冠层氮素含量近似定义为LAI和LNC的乘积，即：

式中，LNC为叶片氮素含量，LAI为叶面积指数。

3.2.2 植被指数计算

除了探讨前人常用的对氮素敏感的光谱指数，本研究还尝试以红边波段代替可见光红波段构建新的指数，从而充分利用Sentinel-2 影像的红边波段，并探讨其用于红树林冠层氮素含量反演的潜力。研究中所选的植被指数及计算公式见表2。

表2 本研究选用的植被指数Table 2 Vegetation indices used in this study

3.2.3 特征变量选取

SVM-RFE 是一种基于支持向量机最大间隔原理的递归特征消除法。该算法的基本思想是利用SVM 算法对所有特征向量进行排序，并评估每个特征向量的重要程度，然后按照后向迭代策略逐一剔除权重最低的特征，直到数据集中剩下期望的特征子集为止（张睿和马建文，2009）。本研究中的特征变量一共有49 个，其中包括9 个Sentinel-2卫星影像波段和40个相关植被指数波段。

续表

3.2.4 核岭回归KRR（Kernel Ridge Regression）

核岭回归算法由Cristianini 和Shawe-Taylor（2000）等提出，是一种基于核方法的非线性回归算法，结合了核技巧和岭回归，不但具有较好的稳定性和泛化能力，还能够有效地解决高维多重共线性问题。核岭回归的基本原理为：将某一原始空间中的n个样本利用一个核函数，令x→φ( )x，映射到高维特征空间F中，数据在高维空间的特征会更容易划分，然后应用岭回归方法，对映射后的数据进行回归处理。回归方程为

相应的代价函数的矩阵形式为

式中，X=[φ(x1)，φ(x2)，…，φ(xn)]为样本数据，w=(w0，w1，…，wn)为权系数向量，λ为正则化参数。

对式（3）中的参数w求偏导，并令其等于0，可得w=(XTX+λI)-1XTy（I为单位矩阵），令α=(K+λI)-1y，核函数K=K(X，X)=φ(x)φ(x)T，则

w其实就是所有样本的一个加权平均，所以新样本的预测值为

本研究将由野外测得的红树林LNC 和LAI 数据计算的CNC 作为模型因变量，3 种Sentinel-2 重建影像中的9 个波段（B2、B3、B4、B5、B6、B7、B8a、B11 和B12）和40 个植被指数共计49 个特征参数为自变量，利用Python 中的sklearn 包进行SVM-RFE特征选择，并利用KRR模型建立红树林CNC的估算模型。

3.3 精度验证

本研究采用留一法交叉验证法LOOCV（Leave-One-Out Cross-Validation）对模型进行验证。评价模型预测能力的指标为：模型决定系数R2（coefficient of determination）、均方根误差RMSE（Root Mean Square Errors）和相对误差RE（Relative Error）。R2反映反演值与实测值之间的线性拟合情况，RMSE 反映反演值与实测值之间的不确定性，RE反映了估算结果的可信度。

4 结果与分析

4.1 红树林叶片和冠层氮素含量分析

受实验条件限制，野外实验数据采集有限，且每个物种采集不均匀。从表3可以看出，不同种类的红树植物叶片氮素含量（LNC）、冠层氮素含量（CNC）和叶面积指数（LAI）存在差异。5 种红树植物LCN 分布在7.2—19.3，平均值为10.6，其中白骨壤LNC最高，平均值达16.4，而桐花树平均LNC 最低，为9.5。5 种红树植物LAI 的范围为0.9—4.2，其中木榄的LAI平均值最大（2.5），其次是桐花树、红海榄和秋茄，最小为白骨壤。而5种红树植物CNC 值分布在8.4—52.6；从各物种CNC平均值看，以榄最大（26.3），其次是桐花树、红海榄、秋茄和白骨壤。5 种红树植物LNC、LAI 和CNC 的变异系数分别为3.2%—17.2%、15.8%—29.5%和16.7%—26.3%，其中均是木榄和红海榄的变异系数较大。总体来说，不同红树物种CNC与LNC有显著差异，且种内CNC变异丰富。

表3 研究区内各红树物种LNC、CNC和LAI值统计分析Table 3 Statistical analysis of LNC，CNC and LAI values of each species of mangroves in the study area

4.2 影像重建质量评价

本研究从定性与定量两个方面分别对3种重建影像方法（Resampled、Sen2Res 和SupReMe）进行评估。对比原始20 m分辨率的影像，目视可见，利用3 种算法重建的10 m 分辨率影像上，红树林、鱼塘和居民地等地物清晰可见，纹理细节明显增强（图2）。表4 反映了将Sentinel-2 影像从20 m 分辨率重建至10 m 分辨率后，影像在信息熵和平均梯度指标值定量变化情况。对比重采样方法可以看出，Sen2Res 算法重建后影像B5、B11 和B12 波段信息熵分别略微下降0.201、0.008 和0.093，其他波段均有所提高；SupReMe 算法重建后影像除B5 波段信息熵略微下降0.083，其他各波段均有所提高。对于平均梯度，Sen2Res 和SupReMe 两种算法重建后影像各波段的平均梯度均有所提升，其中B8a 的增加幅度最大。此外，由Sen2Res 和SupReMe 重建后影像与重采样影像光谱相关性分析（图3 和图4）可看出，两种算法重建后均与原始影像保持了高度相关性（R2均达0.96 以上）。以上分析说明，利用Sen2Res和SupReMe两种方法对Sentinel-2 影像重建，不仅保持了光谱一致性，还提高了影像的清晰度和空间细节。

图3 基于Sen2Res算法的影像重建波段与重采样波段反射率的相关性Fig.3 The correlation of the reconstruction band using Sen2Res algorithm against the resampled band reflectance value

图4 基于利用SupReMe算法的影像重建波段与重采样波段反射率之间的相关性Fig.4 The correlation of the reconstruction band using SupReMe algorithm against the resampled band reflectance value

表4 Sentinel-2影像重建方法评价Table 4 Evaluation results by three methods of Sentinel-2 images reconstruction

图2 影像重建后结果Fig.2 The results of image reconstruction

4.3 植被指数与红树林冠层氮素含量的相关性

图5显示了基于4种影像（原始20 m和3种10 m重建影像）光谱反射率构建的40 个植被指数与对应样本CNC 之间的相关性分析结果。基于4 种影像的植被指数与CNC 的相关性的分布基本相似，30 个植被指数的相关性均达到了显著水平（P<0.05）。3种方法重建影像的相关性均高于原始20 m分辨率影像。具体地，25 个植被指数与CNC 呈正相关关系，4 种影像的MCARI3/MTVI2、MCARI3/MTVI2、SR 及SR 指数分别与CNC 相关性最高，相关系数r分别为0.635、0.733、0.727 和0.732；5 个植被指数与CNC 呈负相关关系，4 种影像的TCARI3/OSAVI 指数与CNC 均具有最高的相关性，相关系数分别为-0.602、-0.725、-0.716和-0.734。由此可知，植被指数对红树林CNC 很敏感，是红树林冠层氮素含量十分有效的预测变量。

图5 基于3种影像重建方法和原始影像的植被指数与红树林冠层氮素含量的相关关系Fig.5 Correlation between vegetation index and canopy nitrogen content of mangroves based on three image reconstruction methods and original image

4.4 红树林冠层氮素含量反演结果

利用SVM-RFE 法分别对4 种影像进行变量选择，基于KRR 方法构建CNC 反演模型，在保留1—49 个特征变量情况下的模型拟合精度和验证精度如图6所示。对于原始20 m 影像和Resampled 方法重建的影像，在保留特征变量数量分别小于13和10 时，R2val 随着特征变量的增加呈波动性变化的趋势，在13 和10 以后呈缓慢增加的趋势，直到34 和27 之后基本不变；对于Sen2Res 和SupReMe方法重建的影像，在保留特征变量数量分别小于20 和10 时，R2val 随着变量的增加呈增加的趋势，在20 和10 之后基本不变。从模型验证精度来看，在保留特征变量数量大于19、22、15 和9 时，4 种影像R2cv基本呈不断下降趋势。

图6 SVM-RFE方法特征变量选择数量与模型精度、验证精度的关系Fig.6 Relationship between the number of selected variables using SVM-RFE method and model accuracy and validation accuracy

综合以上分析，本研究分别选择特征变量数量为34、27、20 和10 为4 种影像数据构建CNC 反演模型，对应的SVM-RFE 方法特征变量选择的结果如表5所示。其中，基于Resampled 影像的模型拟合精度最高（R2val=0.640，RMSE_val=5.064，RE_val=0.179），其次是Sen2Res 方法（R2val=0.630，RMSE_val=5.133，RE_val=0.179）和SupReMe方法（R2val=0.579，RMSE_val=5.473，RE_val=0.194），最低的是原始20 m影像（R2val=0.504，RMSE_val=5.943，RE_val=0.223）。模型交叉验证精度最高的是SupReMe 方法（R2cv=0.501，RMSE_cv=5.958，RE_cv=0.216），其次是Sen2Res 方法（R2cv=0.497，RMSE_cv=5.985，RE_cv=0.214）和Resampled 方法（R2cv=0.493，RMSE_cv=6.008，RE_cv=0.221），最低的也是原始20 m 影像（R2cv=0.333，RMSE_cv=6.893，RE_cv=0.260）。

由表5 可看出，在对4 种影像构建的模型中，都选择了红边波段（B5）、短波红外波段（B11）及植被指数RSSI 和TCARIre1/OSAVI 作为特征变量。此外，红光波段（B4）、近红外—窄波段（B8a）、短波红外波段（B12）以及植被指数（SR和CIG）也多次被选择。这表明，与红树林CNC相关的波段主要集中在红边波段、近红外波段以及短波红外红外波段（Homolová 等，2013）。同时，与“红边波段”相关的植被指数也是构建CNC 反演模型的敏感变量。以Sentinel-2 数据的近红外、绿色和红边波段构建的RSSI 在红树林冠层叶绿素含量估算中也具有良好的应用潜力（Zhen等，2021）；TCARIre1/OSAVI不仅对叶绿素含量变化非常敏感，也能减小植被下垫面、叶面积指数和太阳天顶角变化的影响（Haboudane 等，2002）。由此说明，除了红树林特殊的下垫面外，冠层结构也是红树林CNC 遥感反演的一个重要影响因素（Yu等，2017）。

表5 基于3种重建影像数据SVM_RFE方法特征变量选择策略Table 5 Feature variables selection strategies using SVM RFE based on three image reconstruction methods

采用留—法交叉验证法对以上4种影像数据建立的CNC 反演模型进行验证，实测CNC 值与预测CNC 值散点图如图7所示。总体来说，3 种重建影像的模型预测结果与反演结果均优于原始20 m 影像结果，预测值和真实值较为均匀地分布在1∶1线附近。基于SupReMe 重建影像构建的反演模型取得了最好的预测效果，模型预测精度R2最高为0.501，RMSE和RE均最小，分别为5.958和0.216。Sen2Res 重建影像预测精度与SupReMe 基本相当，而重采样后影像预测精度略差。由图7 可以看出，3 种模型对于CNC 在10—40 范围内的预测结果较好，而当CNC>40时，模型预测结果较低，也就是模型对于高CNC值有饱和现象。

图7 针对3种影像重建数据KRR模型实测CNC值与预测CNC值的拟合结果Fig.7 The fitting results of measured CNC values and estimated CNC values using KRR model for three image reconstruction data

通过以上结果和分析，基于Sen2Res方法重建影像构建的反演模型拟合精度与重采样基本相当，模型预测精度也较高且模型变量选择数量更为合理。因此，本研究基于Sen2Res 方法重建的Sentinel-2 影像，采用SVM-RFE 法选择的20 个特征变量的反演模型对红树林CNC 进行空间制图（图8）。由图8 可知，研究区内红树林CNC 的范围为0—61.8，主要集中在18.8—42.5。从空间分布规律来看，研究区红树林CNC 值分布整体上与红树物种分布相似（高常军等，2022）。研究区西北部中间区域红树林CNC 值较高，东南角区域的CNC 值较低，主要原因是西北部主要分布着树龄大且覆盖度高的木榄、桐花树、红海榄群落；而东南部主要分布着覆盖度较低的白骨壤群落，以及人工种植的较为低矮的木榄和红海榄群落。

图8 基于Sen2Res重建影像的红树林CNC制图Fig.8 The map of mangroves CNC derived from Sen2Res Sentinel-2 data

5 结论

（1）本研究利用Sen2Res和SupReMe两种算法实现了Sentinel-2 影像波段20 m 的分辨率提升至10 m，在保持光谱一致性的同时提高了影像的空间细节。对比原始20 m分辨率影像，2种方法重建后影像大部分波段信息熵均有所提高（除个别波段），平均梯度均有所提升，明显提高了地物的清晰度和空间细节。重建后影像与重采样影像反射率具有高度相关性，所有波段之间的决定系数（R2）均大于0.96。

（2）相关性分析和SVM-RFF 特征选择结果表明，与红树林CNC 相关的波段主要集中在红色波段（B4）、红边波段（B5）、近红外波段（B8a）以及短波红外红外波段（B11 和B12）。植被指数（SR 和CIG），以及与“红边波段”的植被指数（RSSI 和TCARIre1/OSAVI）也是构建CNC 反演的敏感变量。

（3）综合考虑特征变量选择、模型拟合和验证精度，本研究选用Sen2Res算法重建的Sentinel-2影像数据进行研究区红树林CNC 空间制图。基于3种方法重建的10 m影像构建的模型反演精度均优于原始20 m 的影像。3 种影像重建方法（Resampled、Sen2Res 和SupReMe）的特征变量选择数量分别为27、20 和10，重采样后影像的模型拟合精度最高（R2val=0.640，RMSE_val=5.064，RE_val=0.179），其次是Sen2Res （R2val=0.630，RMSE_val=5.133，RE_val=0.179），SupReMe 最低。SupReMe 方法重建后影像模型验证精度最高（R2cv=0.501，RMSE_cv=5.958，RE_cv=0.216），其次是Sen2Res（R2cv=0.497，RMSE_cv=5.985，RE_cv=0.214）。

综合考虑重建影像光谱细节及模型精度，基于Sen2Res 算法重建的Sentinel-2 影像在红树林CNC 反演中具有良好的应用潜力。本研究基于重建的Sentinel-2B 影像首次在国内成功构建了红树林CNC 反演模型，并完成了红树林冠层氮素含量空间制图。Sentinel-2B 卫星影像，尤其是它的“红边波段”和短波红外波段对于红树林氮含量反演具有较强的应用潜力。相比农作物和草地等植被，红树林CNC 遥感反演影响因素更为复杂。尽管本研究考虑了主要冠层结构因素（LAI）的影响，但红树物种、群落结构、叶倾角以及其他生化组分的协同变化等因子对反演结果的影响还有待进一步研究。