基于光学遥感的稀疏乔灌木地上部分生物量反演方法*

石永磊 王志慧 李世明 李春意 肖培青 张 攀 常晓格

(1.河南理工大学测绘与国土信息工程学院 焦作 454003； 2.黄河水利科学研究院水利部黄土高原水土保持重点实验室 郑州 450003； 3.中国林业科学研究院资源信息研究所 北京 100091)

植被生态系统固碳可抵消部分人为的CO2排放，是间接减排的战略选择，也是我国实现“碳中和”目标的重要途径(方精云等， 2011； Xiaoetal.， 2019)。干旱/半干旱面积地区占全球陆地面积的40%以上，碳循环过程对气候变化和人类活动的响应较为敏感，植被碳储量估算不确定性大，已成为碳储量遥感反演研究的热点区域(Smithetal.， 2019)。自1978年以来，“三北防护林、退耕还林草”等生态工程在干旱半干旱地区相继实施，我国人工林面积已跃居世界第一，但该区域降水量少，生态系统脆弱，亟需开展稀疏人工林生态系统结构与功能评估(Zhangetal.， 2015)。地上生物量是表征植被生态系统碳储量的关键指标，更能反映生态系统结构与生态功能优劣(Liu， 2012)。开展稀疏人工林地上生物量遥感反演研究，对干旱地区生态工程效益评估具有重要意义。

植被地上生物量获取的传统方法多以地面实测为主，虽精度较高，但耗时费力且观测站点有限，难以获取植被地上生物量的空间分布信息(冯宗炜等， 1999)。遥感具有宏观、动态、可重复等特点，已成为区域植被地上生物量反演的主要技术手段(李德仁等， 2012； Liu， 2012)。根据遥感平台搭载传感器类型，可分为光学遥感(Liuetal.， 2014； Pengetal.， 2019)、微波遥感(陈尔学， 2007； Nietal.， 2014)、激光雷达遥感(庞勇等， 2012； Latifietal.， 2015)。不同类型遥感数据反演地上生物量的物理机制不同，且各有优缺点(刘茜等， 2015)。由于中低分辨率光学遥感数据的优势，光谱指数-地上生物量统计模型是目前最简单易行的方法。但光谱信号对植被冠层结构变化存在饱和效应，且会受到辐射传输过程中各种环境因素影响，因而光谱指数-地上生物量统计模型的普适性较差。针对这一问题，前人先后将降雨(Shoshanyetal.， 2011)、林龄(Liuetal.， 2014)、纹理(庞勇等， 2017)、地形(Pengetal.， 2019)、树种类型(Jiangetal.， 2020)等信息融入光谱指数-地上生物量模型，应用结果表明可提高茂密森林地上生物量反演精度。但是，干旱/半干旱地区天然林或人工林较为稀疏，遥感观测混合像元问题导致光谱指数-地上生物量模型的反演精度较低(Zandleretal.， 2015； Gaoetal.， 2020)。为解决这一问题，部分学者利用高空间分辨率遥感数据首先获取植被冠层覆盖度，再构建覆盖度-地上生物量模型反演地上生物量(Suganumaetal.， 2006； Ozdemir， 2008； 张华等， 2012； Wangetal.， 2016)，但该类方法大多应用于纯乔木林或纯灌木林区域，由于该区域木本植被类型和草本植被覆盖度背景具有高异质性，乔灌木地上部分生物量遥感估算难度较大。探究基于分层方案的光学遥感的稀疏乔灌木地上部分生物量反演方法，可为干旱区人工林乔灌木地上部分生物量高精度反演提供理论依据。鉴于此，本研究基于无人机数据采用3种分层方案构建冠层盖度-乔灌木地上部分生物量模型以及基于Landsat8 OLI数据采用3种分层方案构建不同光谱指数-乔灌木地上生物量模型，对比分析不同分层方案的乔灌木地上部分生物量模型精度，以期为基于遥感数据的干旱区人工林乔灌木地上部分生物量高精度反演提供科学理论依据。

1 研究区概况

毛乌素沙地位于陕西省榆林市与内蒙古自治区鄂尔多斯市之间(107°20′—111°30′E，37°27′—39°22′N)，面积约4万km2，是中国四大沙地之一。潜在蒸散多年均值为1 048.81 mm，年均气温6.0～8.5 ℃，年均降水量250～400 mm，降水集中在7—9月份。由于历史上剧烈人类活动，该地区植被严重退化(白壮壮， 2019)，风蚀沙化和水土流失严重，生存环境恶劣。经过近30年植树造林，乔灌植被显著增加，实现了“人进沙退”的治沙奇迹(郑玉峰， 2015)。该区域人工乔木主要有新疆杨(Populusalbavar.pyramidalis)、旱柳(Salixmatsudana)和樟子松(Pinussylvestrisvar.mongolica)，灌木主要有沙蒿(Artemisiaordosida)、沙柳(Salixpsammophila)和柠条锦鸡儿(Caraganakorshinskii)等。

2 研究方法

2.1 数据采集与处理

2.1.1 野外样地设置与调查 2019年9月10—30日，在毛乌素沙地乌审旗境内开展了20天野外调查，共随机调查102块30 m×30 m样地，采用Garmin GPSMAP 669 s测量样地中心点坐标，该设备具有接收CORS(continuous operational reference system)信号的功能，且研究区植被稀疏，地势平坦，可保障GPS测量坐标准确度。测量样地内每株乔木的树高和胸径，每株灌木的高度和冠幅(东西向与南北向)，并记录样地内乔木和灌木的树种信息，所测的102块样地内均至少包含1种乔木或灌木树种。同时使用大疆无人机(Phantom 4 Pro)垂直拍照样地，照片重叠度大于60%，利用无人机后处理软件Photoscan对照片进行拼接和正射校正处理。

2.1.2 Landsat遥感卫星数据 从Google Earth Engine云平台下载2019年9月21日的研究区Landsat8 OLI地表反射率数据，空间分辨率为30 m，且基于F-mask算法对影像中的云和云阴影进行标识。

2.2 基于异速生长方程的人工林地上生物量计算

由于该地区禁止砍伐树木，因此，基于实地调查树种类型及其胸径、冠幅，选用乔灌木树种异速生长方程(Liuetal.， 2014； 李海奎等， 2010； Guoetal.， 2021)计算单株乔灌木地上生物量(kg)。将样地内所有单株木的地上部分生物量求和并除以样地面积，即得样地乔灌木地上部分生物量(t·hm-2)。

2.3 稀疏乔灌木地上部分生物量遥感反演方法

本研究总体技术路线见图1。

图1 基于光学遥感数据的稀疏乔灌木地上部分生物量反演技术路线图

2.3.1 基于冠层覆盖度-乔灌木地上部分生物量模型的分层建模方案 根据前人经验，HSV颜色空间变化可有效提高基于无人机RGB图像的植被冠层识别精度(Yanetal.， 2019)，因此将原始无人机RGB图像进行HSV颜色空间变化。

先根据无人机图像中的样地边线对样地区域图像进行裁剪，再利用人工目视解译从裁剪后的样地图像中对地表覆盖类型样本进行选取，选取的样本分别为乔木、灌木、草地、裸地和阴影，并按6∶4比例将样本集分为训练样本与验证样本。基于eCognition软件平台，首先用多尺度分割算法分割无人机正射图像，提取分割对象的纹理、亮度、饱和度等特征，再用支持向量机(SVM)算法和训练样本对样地图像分类。基于验证样本构建误差矩阵，评价分类精度。冠层覆盖度为植被像素个数占样地总像素个数的比值。由于无人机影像中存在大量阴影，阴影内主要包含草本植被和裸土，本研究采用归一化方法对草本植被覆盖度进行修正，计算公式如下：

(1)

式中： GN为草本植被覆盖度修正值，G、BS和S分别为原始草本、裸土和阴影在样地内的面积占比。

将样地图像的冠层覆盖度作为自变量，异速生长方程计算的乔灌木地上部分生物量作为因变量，两边取对数后建立冠层覆盖度-乔灌木地上部分生物量统计模型：

lnAGB=alnC+b。

(2)

式中：C为乔灌木冠层覆盖度；AGB为乔灌木地上部分生物量(t·hm-2)；a和b均为系数。

为探讨稀疏乔灌混交林木本植被类型异质性对冠层覆盖度-乔灌木地上部分生物量模型精度的影响，采用3种分层方案构建冠层覆盖度-乔灌木地上部分生物量模型： 1) 不分层建模，即用样地乔灌木总覆盖度与总生物量构建乔灌木地上部分生物量模型； 2) 乔木和灌木分别建模； 3) 5个树种(樟子松、新疆杨、旱柳、沙蒿和柠条锦鸡儿)分别建模。对比分析不同分层方案的乔灌木地上部分生物量模型反演精度。

2.3.2 基于光谱指数-乔灌木地上部分生物量模型的分层建模方案 根据前人经验，选取与植被地上生物量相关性较强的6种光谱指数： 归一化植被指数(NDVI)(Rouseetal., 1973)、比值植被指数(RVI)(Pearsonetal., 1972)、改进型调整植被指数(MSAVI)(Qietal., 1993)、植被湿度指数(NDMI)(Gaoetal., 1996)、纯植被近红外反射率(NIRv)(Badgleyetal., 2017)和缨帽变化绿度分量(TCG)(李博伦等， 2016)。

NDVI=(NIR-R)/(NIR+R)；

(3)

RVI=NIR/R；

(4)

MSAVI=0.5×

(5)

NDMI=(NIR-SWIR1)/(NIR+SWIR1)；

(6)

NIRv=NDVI×NIR。

(7)

式中： NIR、R和SWIR1分别为Landsat8 OLI数据的近红外波段(0.845～0.885 μm)、红光波段(0.63～0.68 μm)和短波红外波段(1.56～1.66 μm)的地表反射率。

提取样地中心点GPS坐标所对应的Landsat像元地表反射率并计算相应光谱指数，将光谱指数作为自变量，异速生长方程计算的生物量作为因变量，构建光谱指数-乔灌木地上部分生物量统计模型：

AGB=a×VI+b。

(8)

式中：VI为光谱指数。

为探讨稀疏乔灌混交林草本覆盖度背景异质性对光谱指数-乔灌木地上部分生物量模型精度的影响，采用3种分层方案构建光谱指数-乔灌木地上部分生物量模型： 1)不分层建模，即所有样地的光谱指数与乔灌木地上部分生物量构建模型； 2)有草本植被与无草本植被的样地分别建模； 3)3种草本覆盖度等级(0%、0～30%、>30%)的样地分别建模。对比分析同一分层方案下，不同光谱指数反演乔灌木地上部分生物量的精度差异，以及同一光谱指数在不同分层方案下反演乔灌木地上部分生物量的精度差异。

3 结果与分析

3.1 样地无人机影像分类

无人机影响与基于面向对象的支持向量机分类结果如图2。利用误差矩阵法评价无人机影像的总体分类精度、各覆盖类型的生产者精度和用户精度。所有样地无人机影像分类精度评价指标的分位数如图3所示，表明基于面向对象的支持向量机分类可获取高空间分辨率无人机影像的高精度地表覆盖分类结果，从而计算得到高精度植被冠层覆盖度。

图2 无人机影像(a)与基于面向对象的支持向量机分类结果(b)

图3 所有样地无人机影像分类精度指标分位数统计图

3.2 冠层覆盖度-乔灌木地上部分生物量模型分层建模反演结果

3种分层方案所构建的样地冠层覆盖度-乔灌木地上部分生物量模型见表1。图4为乔灌木地上部分生物量模型模拟结果与实测值对比结果。不同分层方案构建的模型均通过了P<0.01的显著性检验，但R2存在差异。不分层方案所构建的模型鲁棒性最差(R2= 0.22)，且反演结果精度最低(RMSE = 14.98 t·hm-2)。与不分层的冠层覆盖度-乔灌木地上部分生物量模型相比，基于乔木和灌木2种植被类型分层建模(RMSE = 7.44 t·hm-2)和基于5个树种分层建模(RMSE = 5.82 t·hm-2)的反演误差分别减少了50.32%和61.1%。结果表明，冠层覆盖度-乔灌木地上部分生物量模型方法对木本植被类型较为敏感，模型反演精度会随木本植被类型分类精细化程度增加而显著提高。

3.3 光谱指数-乔灌木地上部分生物量模型分层建模反演结果

在3种分层方案下，用6种光谱指数分别构建光谱指数-乔灌木地上部分生物量模型(表2)，图5为不同模型的精度验证结果。由图5可以看出，NIRv-乔灌木地上部分生物量模型的精度最高(3种分层方案平均RMSE = 7.25 t·hm-2)，这表明NIRv对稀疏乔灌木地上部分生物量变异的解释能力优于其他光谱指数； 常用的NDVI-乔灌木地上部分生物量模型反演精度最低(3种分层方案平均RMSE = 9.43 t·hm-2)； 不同光谱指数对稀疏乔灌木地上部分生物量变异的解释能力排序为NIRv>NDMI>TCG>MSAVI>RVI>NDVI。光谱指数-乔灌木地上部分生物量模型反演精度会随草本覆盖度等级分类精细化程度的增加而增加。其中，NIRv-乔灌木地上部分生物量模型对草本覆盖度异质性最敏感，分层方案使RMSE减少了16.62%，而NDMI-乔灌木地上部分生物量模型对草本覆盖度异质性最不敏感，分层方案仅使RMSE减少了8.13%。不同光谱指数-乔灌木地上部分生物量模型对草本覆盖度背景的敏感性排序为NIRv>TCG>NDVI>MSAVI>RVI>NDMI。

表1 3种分层方案所构建的冠层覆盖度-乔灌木地上部分生物量模型

图4 3种分层方案下冠层覆盖度-乔灌木地上部分生物量模型反演精度

4 讨论

前人研究已表明基于高空间分辨率遥感的冠层覆盖度-乔灌木地上部分生物量模型方法可有效消除草本植被盖度背景对稀疏乔灌木地上部分生物量反演精度的影响(彭守璋等， 2010； 张华等， 2012； Wangetal., 2015)，但对稀疏乔灌混交林区域，木本植被类型异质性会严重影响冠层覆盖度-乔灌木地上部分生物量模型反演精度。这主要是由于干旱地区乔木和灌木的冠幅与树高之比存在明显差异，导致二者的冠层覆盖度与乔灌木地上部分生物量的响应关系差异显著。若基于乔灌木混合数据(即不分层)构建乔灌木地上部分生物量模型，必然降低模型鲁棒性与反演精度(Caoetal.， 2014)。因此，对稀疏乔灌混交林，必须考虑乔木和灌木2种植被类型才能使冠层覆盖度-乔灌木地上部分生物量模型方法精度达到实用要求。

木本植被类型的冠层覆盖度-乔灌木地上部分生物量模型精度更高，但基于高空间分辨率遥感提取乔木、灌木两种植被类型难度较大，且不适合大范围推广应用。本研究继续探讨了基于Landsat卫星数据的光谱指数-乔灌木地上部分生物量模型在稀疏乔灌混交林区域的适用性。大量研究已证实纯植被冠层近红外反射率NIRv与植被初级生产力(GPP)有较强相关性(Badgleyetal.， 2017)，本研究发现NIRv与植被地上部分生物量的相关性也优于其他光谱指数，这主要是因为NIRv指数能反映更多植被冠层结构信息，同时能消除草本植被覆盖度背景对冠层光谱的部分影响(Badgleyetal.， 2017)。

根据植被辐射传输过程的物理机理，影响卫星光谱信号的关键植被结构参数为叶面积指数，而不是冠层覆盖度(Heetal.， 2019)，因此光谱指数-乔灌木地上生物量模型反映的是叶面积指数与乔灌木地上部分生物量之间的生物物理响应关系。考虑草本覆盖度背景异质性的分层方案可有效提高光谱指数-乔灌木地上部分生物量模型的反演精度，但不同光谱指数-乔灌木地上部分生物量模型对草本覆盖度背景的敏感性存在差异。NIRv是NDVI与近红外反射率的乘积，而NDVI受草本植被绿度影响较大，因此草本植被变化对NDVI和NIRv反演乔灌木地上部分生物量精度影响较大。而干旱地区草本植被叶片含水量较低，NDMI对草本覆盖度变化并不敏感，因此，草本植被覆盖变化对NDMI反演精度影响较小。

表2 3种分层方案所构建的光谱指数-乔灌木地上部分生物量模型

图5 3种分层方案下6种光谱指数-乔灌木地上部分生物量模型反演精度

本研究所用实测乔灌木地上部分生物量数据是由前人构建的乔木、灌木异速生长方程计算得到的，与真实乔灌木地上部分生物量的绝对值存在一定偏差，但该偏差为系统误差，并不影响不同乔灌木地上部分生物量模型的精度评价和相互对比关系。本研究利用面向对象的机器学习算法分类地表覆盖精度较高，但分类误差仍会影响乔灌木地上部分生物量反演精度，需结合深度学习算法，以期降低由图像分类误差引起的乔灌木地上部分生物量反演不确定性。本研究仅利用了统计回归方法对乔灌木地上部分生物量遥感反演进行研究，今后需综合利用混合像元分解与辐射传输模型，研究光学卫星观测信号与植被乔灌木地上生物量的物理响应关系(Jos′eetal.， 2021； 陈晋等， 2016； 柳钦火等， 2016)，有机结合微波和激光雷达数据，提高卫星反演乔灌木地上生物量方法的普适性与精度水平。

5 结论

1) 木本植被类型对冠层覆盖度-乔灌木地上部分生物量模型精度影响较大，考虑木本植被类型的分层方案可显著提高反演精度，且反演精度随分层精细程度增加而增加，至少需区分乔木和灌木两种木本植被类型才能保证冠层覆盖度-乔灌木地上部分生物量模型在乔灌混交林区域的反演精度满足实用需求。

2) 基于草本植被覆盖度的分层方案可显著提高光谱指数-生物量模型反演精度，但每种光谱指数-乔灌木地上部分生物量模型对草本覆盖度背景的敏感性又存在差异。

3) 高空间分辨率遥感可用于获取木本植被类型及其草本覆盖度背景信息等先验知识，基于Landsat卫星数据的考虑草本植被覆盖度背景的分层方案的NIRv-乔灌木地上部分生物量统计模型适于大区域稀疏乔灌木地上部分生物量的遥感反演。