联合Sentinel-1 与Sentinel-2 数据的青藏高原东缘草地地上生物量反演

孙 剑，杜 忠，林用智，王 杰

（西华师范大学地理科学学院, 四川 南充 637009）

草地资源具有重要的经济价值和生态功能，在畜牧业发展、生物多样性保护、缓解全球气候变化等方面起着至关重要的作用[1-4]。青藏高原的高寒天然草地受人类干扰和气候变化的影响较大。近年来，供水不足、过度放牧、草地利用率低下等原因导致高寒草地出现了不同程度退化现象[5]。灌丛化的加重，导致牧草生产减少，草地生态系统逐步遭到破坏。因此，利用遥感技术进行大面积、及时有效地评估草地地上生物量(aboveground biomass, AGB)，对高寒草地及青藏高原生境质量的保护具有重要现实意义。

目前，光学遥感数据已被广泛用于基于统计模型的生物量估算[6-7]。但光学影像易受云雾影响，使得传感器无法准确获取地物信息，且在中度或高度的植被密度中存在过饱和问题[8]。相较于光学传感器，合成孔径雷达(synthetic aperture radar, SAR)波长较长，不受云雾影响，对地表植被有一定的穿透性，能较为准确地获取植被的结构信息。因此，SAR 数据逐渐成为草地生物量研究中的重要数据源[9]。Moreau 和Toan[10]发现C 波段雷达数据 ERSSAR 的后向散射系数对植被生物量很敏感，并与草地生物量之间建立对数统计回归模型，反演了安第斯湿地牧草生物量。Santos 等[11]使用P 波段多极化AirSAR 数据，采用指数模型和多项式模型成功反演了巴西亚马逊河流域热带雨林的主生林和次生林植被生物量。随着遥感技术和研究方法的发展，有效地结合光学数据和微波数据的优势，利用经验模型、多元线性回归等方法反演植被生物量已成为一种研究趋势。同时，高空间分辨率的SAR 和多光谱数据为提升AGB 反演精度提供了更多可能性，如欧空局(https://www.esa.int/)哨兵系列卫星数据，其中Sentinel-1 为SAR 数据，Sentinel-2 为多光谱数据，两者分辨率可增强到10 m，适用于AGB 的估算[12]。潘磊等[13]结合Sentinel-1 与Sentinel-2 数据，以福建省三明市将乐国有林场为研究区，利用多元线性逐步回归方法，探索联合Sentinel-1 与Sentinel-2 在估算森林AGB 方面的优势，其调整决定系数(R2adj)达到了0.575，相较于仅使用Sentinel-2数据提升了0.064。而Chang 和Shoshany[14]结合Sentinel-2 NDVI及Sentinel-1 VH 极化的后向散射系数建立的半经验物理模型，联合估算了地中海灌木林AGB，融合模型的决定系数(R2)达到了0.866，也取得了很好的反演精度。综上，基于仅使用光学数据或微波数据估算草地AGB 已有大量研究，但多模态数据和相应特征参数的选择对草地AGB 反演建模仍有发展空间，另外关于草地AGB 估算的研究方法多种多样，但旨在对比分析不同模型反演精度和区别的研究较少。

因此，本研究以红原县为研究区，联合Sentinel-2数据、Sentinel-1 数据，利用多元线性回归、半经验物理模型等方式建立草地AGB 遥感反演模型，用留一交叉验证法对反演精度进行评价，最终得到了红原县草地AGB，实现了红原县草地AGB 的数字制图，以期为优化AGB 遥感反演模型、青藏高原地区的草地保护与生态建设等提供科学依据和技术支持。

1 材料与方法

1.1 研究区概况

红原县地处四川省西北部、阿坝藏族羌族自治州的中部，位于101°51′～103°22′ E，31°51′～33°33′ N，属青藏高原的东南缘。年平均气温为1.1 ℃，无明显季节变化。境内平均海拔3 500 m 左右，红原县地势由东南向西北倾斜，地貌由山原向丘状高原过渡，小盆地分布广泛。红原以寒草地和亚高山草地为主，境内天然草地面积达77.6 × 104hm2，占县境总面积的92.97%，属纯牧业县[15]。

1.2 数据来源及预处理

1.2.1 样本采集与处理

野外数据采集时间为2021 年8 月。共设置了10 个样区(图1)，样区的大小为1 000 m × 1 000 m，每个样区内设置3～5 个样点。由于样区内草种类型较多，且空间分布不均匀，在一个样点(10 m × 10 m)大小范围内按等边三角形布设3 个0.5 m × 0.5 m的样方。实测时使用GPS 记录3 个样方中心(等边三角形中心)坐标，并记录实测点位海拔、土壤类型和质地、草样本的株高等参数。草地AGB 的获取方法：在设定样方内齐地刈割，并现场记录其鲜重，返回实验室后用烘箱烘干至恒重，并记录干草重，最后求取3 个样方干草重量的平均值作为该样点的生物量。本次在10 个样区内共设定132 个样方，采集草地AGB 样本共44 个，剔除数据处理过程中造成AGB 损失的数据，确定样本数据40 个。从每个样区均选取2～4 个样点(共30 个样点)用于构建反演模型，其余10 个样点用于模型精度检验。

图1 研究区位置及样区示意图Figure 1 Location of the study area and schematic diagram of the sample area

1.2.2 遥感数据预处理

Sentinel-2 数据获取及预处理。Sentinel-2 影像为光学影像，受云层干扰较大。本研究以红原试验区域同一时期4 景云量较少的Sentinel-2 影像(云量 < 10%)为研究数据，数据来源于美国USGS 网站(https://earthexplorer.usgs.gov/)。使 用 欧 洲 空 间 局SNAP 8.0 对L2A 级数据(已经过大气校正)进行超级分辨率增强，将影像分辨率增强至10 m，然后导入ENVI 5.3 进行镶嵌、裁剪等操作。处理完成后，使用B12、B8A、B2 波段合成显示，结果如图2(a)所示，从波段合成图可看出研究区内无大片云雾遮挡，色彩对比明显，数据质量较高。

图2 Sentinel 影像预处理结果Figure 2 Results of Sentinel image pre-processing

Sentinel-1 数据获取及预处理。Sentinel-1 卫星由两颗极轨卫星A 星和B 星组成，搭载有C 波段的合成孔径雷达(SAR)，其空间分辨率为10 m × 10 m，包含VH 和VV 两种极化方式。为了减少由天气因素、Sentinel-2 与Sentinel-1 过境时间不同而引起的误差，选阿拉斯加数据中心(https://search.asf.alaska.edu/)提供的与Sentinel-2 同日过境的Sentinel-1 GRD 产品为研究数据，该影像覆盖了整个研究区，并用SNAP 8.0 对GRD 数据进行预处理和信息提取。将提取的VV、VH 极化的C 波段数据导入ENVI 5.3 进行裁剪，并交叉合成显示，结果如图2(b)所示。

植被指数和叶面积指数。植被指数(vegetation index，VI)能反映植被生长和覆盖状况，可作为植被覆盖及变化监测的标志[8]，能从宏观上快速准确地分析草地生长状况。本研究从Sentinel-2A 影像中获取了5 个植被指数(表1)。草地AGB 与植被指数之间有着良好的正相关关系[16]。因此，可用这5 个植被指数与实测样本相结合，构建AGB 估算模型。

表1 Sentinel-2 植被指数及计算公式Table 1 Sentinel-2 vegetation index and calculation formula

叶面积指数(leaf area index，LAI)是监测植被生长状况和预测产量的重要变量，也常被用于全球空间尺度的作物生长监测、作物病虫害监测和产量估算等[17]。针对Sentinel-2 影像，SNAP 8.0 提供了特定的模块来获取叶面积指数，步骤如下：先将影像重采样到10 m；再使用SNAP 8.0 内置的生物物理处理器导出LAI。这个生物物理过程使用了Sentinel-2影像的8 个波段(B3、B4、B5、B6、B7、B8A、B11、B12)并考虑了太阳天顶角、方位角等参数。SNAP 反演的LAI 和实测的LAI 总体一致，表明可用于大规模的植被监测[18-19]。

NIR代表Sentinel-2影像的B8 (近红外)波段，RED代表B4 (红)波段，BLUE代表B2 (蓝)波段，GREEN代表B3 (绿)波段。

NIR represents the B8 (Near Infrared) band of Sentinel-2 imagery; RED represents the B4 (Red) band; BLUE represents the B2 (Blue) band; GREEN represents the B3 (Green) band.

1.3 研究方法

1.3.1 回归分析

回归分析(regression analysis)是生物量遥感估测中使用较广泛的分析方法，能够表现两个(多个)相关变量间的定量关系。其中多元线性回归模型(复线性回归模型)通常用于研究一个因变量与多个自变量的函数关系。其数学表达式为：

式中：Y是因变量(生物量的值)，a0是常数项，X是自 变 量(遥 感 影 像 因 子)，X1,…,Xn是 随 机 的 变 量，a1,…,an为回归系数，C为误差项。

逐步回归的基本方法是将变量逐一引入模型，并对引入变量计算偏回归平方和，经检验后剔除不显著的自变量，从而建立最优的多元线性回归模型[20]。

1.3.2 半经验物理模型构建

半经验物理模型综合了经验统计模型和物理模型的优点。模型所用的参数既有经验参数，又具有一定的物理意义，且模型表达也较为简洁。在宏观上，草地的分布较为密集，植被覆盖度高。相较于Chang 和Shoshany[14]研究的灌木，红原县草地植被更浅，8 月份草层高度大部分介于0.15～0.45 m。因此，利用Sentinel-1 影像C 波段能穿透草层到达地表的特点，便能获取草层的结构信息。但草地上植被覆盖度(fractional vegetation cover, FVC)的空间差异不大，故在灌木AGB 估算模型中调整FVC 参数为LAI，构建了如下模型：

式中：AGB为预测草地AGB，f(S2)为根据Sentinel-2影像构建的草地AGB 模型，f(S1)为根据Sentinel-1影像构建的体积回归函数，LAI 为叶面积指数，e为系数。

1.3.3 精度评价

用均方根误差(root mean squared error，RMSE)、平均相对误差(mean relative error，MRE)和决定系数(R2)来评价反演预测值与实测值之间的精度和可信度。

式中：Yi、Yi′分别为实测值、反演值，n为样本数量，N为预留样本数。R2反映反演值与实测AGB 之间的拟合程度，当R2越趋近1，表示反演AGB 与实测AGB 相关性越好。

2 结果与分析

2.1 草地地上生物量的模型构建

线性回归模型。经试验发现，研究区内，Sentinel-1 VH dB 值与AGB 之间的相关性优于VV dB 值与AGB 之间的相关性，于是采用VH dB 值构建模型(表2)；将Sentinel-2 影像与采样点位置相结合，提 取 出 对 应 位 置 的NDVI、DVI、RVI、EVI、GNDVI，并结合实测生物量，采用最小二乘回归的方法分别构建红原县草地AGB 与上述各个植被指数之间的线性模型(表2)。

表2 回归模型建模精度Table 2 Modeling accuracy of the regression model

结果显示，基于单一指标的线性回归模型中，NDVI、RVI、GNDVI 以及LAI 与生物量之间有较高的拟合性，其中，RVI 模型的拟合性最高，决定系数R2达到0.61。多元回归模型中，用Sentinel-2 影像提取的6 个因子建立的多元线性回归六参数模型以及Sentinel-2 影像提取的6 个因子结合VH dB值建立的多元线性回归七参数模型，都有着较好的拟合精度，在加入VH dB 值后，模型的R2提升了0.09，RMSE 减小了9.37 g·m-2，MRE 减小了21.82%。采用逐步回归保留了RVI 和VH dB 值两个参数，建立的模型相较于AGBRVI模型，决定系数R2提升了0.17，RMSE 减 小 了14.12 g·m-2，MRE 减 小 了26.56%。

半经验物理模型。通过比较Sentinel-2 影像的5 个AGB 估算模型，选取最佳的AGBRVI模型，带入公式(2)，替换f(S2)参数，再将公式(5)带入公式(2)，替换参数f(S1)，便得到了一个包含5 个系数[其中f(S2)含两个系数，f(S1)含两个系数，再加上系数e]的半经验物理模型，通过IDL 8.5 程序设计平台，运用曲线拟合方式，以30 个样本作为训练数据，解算出5 个系数，分别为110.945、7.432 55、-21.071 6、-1.361 20、0.145 245,得 到 红 原 县 草 地AGB 估算模型[公式(6)]。结果显示，基于Sentinel-2 和Sentinel-1 的半经验物理模型，R2达到0.77，RMSE 为43.88 g·m-2，MRE 为11.00%。

2.2 模型反演生物量精度验证

为验证模型精度，本研究选用剩余的10 个样本作为验证数据，并以Sentinel-2 多元回归六参数模型、Sentinel-2 结合Sentinel-1 多元回归七参数模型、Sentinel-2 结合Sentinel-1 逐步回归模型以及半经验物理模型4 个较优草地AGB 反演模型为代表，获取10 个样本对应的草地AGB 反演值，将其与实测数据进行对比分析，使用误差值与实测值的百分比大小来评价反演的精度。

在Sentinel-2 建立的多元回归六参数模型中 (图3)，反演值与实测值之间的误差有7 个点在实测值的10%以内，有两个点在10%～20%，有1 个点在20%～30%，实测值与估测值对比分析结果的R2为0.644 8，RMSE 为37.32 g·m-2，MRE 为8.20%；以Sentinel-2 结合Sentinel-1 建立的多元线性回归七参数模型，有4 个点在实测值的10%以内，有5 个点在10%～20%，有1 个点在20%～30%，实测值与估测值对比分析结果的R2为0.714 5，RMSE 为38.99 g·m-2，MRE 为9.93%；以Sentinel-2 结合Sentinel-1 逐步回归模型，模型的反演值与实测值之间的误差有7 个点在实测值的10%以内，有3 个点在10%～20%，实测值与估测值对比分析结果的R2为0.789 1，RMSE 为33.49 g·m-2，MRE 为8.11%；用Sentinel-2结合Sentinel-1 建立的半经验物理模型，模型的反演值与实测值之间的误差有5 个点在实测值的10%以内，有5 个点在10%～20%，实测值与估测值对比分析结果的R2为0.821 4，RMSE 为36.89 g·m-2，MRE为8.37%。

图3 模型反演生物量结果精度验证Figure 3 Accuracy validation of model retrieval biomass results

总体上，上述的4 个模型总体精度均较好，在加入Sentinel-1 的影响因子后，与只用Sentinel-2 建立的多元回归六参数模型相比，Sentinel-1 结合Sentinel-2 多元回归七参数模型精度更高，可用于红原县地区草地AGB 的估算，为草地生物量提供精确的评估和研究。

2.3 试验区域草地反演地上生物量制图

根据上述研究结果，采用模型反演精度较高的Sentinel-2 多元回归六参数模型、Sentinel-2 结合Sentinel-1 多元回归七参数模型、Sentinel-2 结合Sentinel-1 逐步回归模型以及半经验物理模型反演研究区8 月份草地AGB (图4)。制图结果与草地实际分布状况一致，在红原县境内大部分区域均有草地的分布，尤其以中北部丘状高原最为典型。同时表现出明显的空间分布差异：沙化地、裸地、水体内草地AGB 接近于0；东南山原草地AGB 大部分低于100 g·m-2；4 种模型反演的生物量制图结果整体一致，均能准确展现红原县草地AGB 的空间分布及其差异。

3 讨论

从红原县草地AGB 反演结果可以看出，草地AGB 较高值分布在中部和北部丘状高原地区，其原因在于该区域有优质的牧草基地，为保护草地，合理利用草地资源，当地实行草地围栏，将草地划分为冬夏两季牧场，极大程度上减小了放牧压力，有利于推进牧场畜牧业可持续发展。东南山原草地AGB 较低，可能与东南山地海拔高差较大、地势较陡有关，不利于牧草的生长。另外，在适合草地生长但草地AGB 较低的区域，建议合理确定载畜量，实行合理放牧，继续加强防沙治沙、灭鼠种草等草地保护措施，研究不同栽种模式的固沙效果，改善土壤有机质及土壤水分状况。通过人力保护措施，使草地能够继续发挥防风固沙、涵养水源、净化空气的作用，提高草地植被盖度和草地生产力，实现红原县草地生态修复和畜牧业可持续发展。

本研究与前人的研究都证明了多光谱数据在估算AGB 方面有较强的应用潜力[6-7]，但受遥感技术、数据质量和建模方法等的限制，AGB 反演效果仍有提升空间。本研究中基于Sentinel-2 多光谱数据单一指标的线性回归模型能较好地反演AGB，结合Sentinel-1 与Sentinel-2 特征参数的反演模型能进一步提升反演精度，但对于生物量较大的样点，仍存在预测草地AGB 值小于实测草地AGB 值的问题(图3)，可能是Sentinel-1 影像的C 波段对植被生长旺盛、植被覆盖度较大区域的穿透能力有限，不能较好地解决Sentinel-2 多光谱影像反演生物量存在的过饱和问题。后续研究可以考虑使用波长更长的L 波段或P 波段来反演草地AGB。

经验模型不涉及机理过程[21]，与简单经验统计模型相比，半经验物理模型所需参数较多，但其以物理机理为依据，并结合区域特征性参数，能较好地反映区域实际情况。模型精度主要受自变量与因变量间相互关系的显著程度的影响。如本研究中仅根据实测的草地AGB 参数与遥感影像特征参数构建的经验模型，按照误差最小的原则筛选出构建模型所需的显著相关性因子，从而得到较好的模型反演效果，这在前人的研究中有所体现[22-23]。而半经验模型既具有一定的物理学意义，又通过引入经验系数改进模型假设的不足，进而提升模型效果[24]。如Chang 和Shoshany[14]提出的基于植被的生长形态形成的融合生物量模型(半经验物理模型)的准确率较传统的融合模型高14%左右。本研究也根据红原县的草地植被形态结构调整了Chang 和Shoshany[14]构建的反演模型中的FVC 参数，再结合草地植被的叶面积指数构建了半经验物理模型，从而达到和经验模型等同的效果，证明了半经验物理模型在草地AGB 反演中的实用性。

目前采用的反演模型相对传统，建模的数据来源较单一，后续可以基于Lidar、SAR、多光谱影像等遥感数据，使用机器学习方法以多模态遥感信息融合方式进行建模，挖掘多模态数据在估测草地地上生物量的潜力。另外，青藏高原东缘多为丘状高原，区域内湿地较多，地形起伏、土壤水分对雷达影像后向散射系数的影响较大，今后可根据研究区的实际情况，将地形、土壤水分的影响因素加入模型，进一步细化草地地上及地下生物量的反演研究，这对红原县乃至整个青藏高原的草地生态质量监测及其固碳效应研究等都具有重要意义。

4 结论

本研究以青藏高原东缘阿坝藏族自治州红原县为研究区，结合 Sentinel-1、Sentinel-2 主被动遥感影像和实测样地数据，采用单一指标线性回归、多元线性回归、逐步回归及半经验物理模型方式建模，对红原县草地AGB 的反演模型效果和空间分布进行研究。结果表明：Sentinel-2 多光谱数据单一指标线性回归模型、Sentinel-2 多光谱数据多元线性回归模型与草地AGB 之间有较高的拟合性；结合Sentinel-1影响因子VH dB 值建立的多元回归模型能进一步提高建模精度，多元线性回归模型的模型精度R2从0.74 增加到了0.83；逐步回归模型R2达到了0.78，RMSE 为42.21 g·m-2；半经验物理模型R2达到了0.77，RMSE 为43.88 g·m-2，MRE 为11.00%，总体效果较好，精度可靠，可用于红原县草地AGB 的估算；草地AGB 制图结果显示，红原县草地AGB 在中部和北部丘状高原地区较高(300 g·m-2以上区域较多)，东南山原草地AGB 较低(100 g·m-2以下区域较多)，沙化地、裸地、水体内生物量接近于0。