极干旱区稀疏荒漠植被地上生物量遥感估算

叶静芸， 吴 波,2， 贾晓红,2， 费兵强,2， 高君亮，成 龙,2， 庞营军,2， 姚 斌,2， 孔德庸

（1.中国林业科学研究院荒漠化研究所，北京 100091；2.国家林业和草原局荒漠生态系统与全球变化重点实验室，北京 100091；3.中国林业科学研究院沙漠林业实验中心，内蒙古 磴口 015200；4.韶关学院，广东 韶关 512005）

我国干旱区总面积约为2.15×106km2，约占国土面积的22%。包括新疆全境、甘肃河西走廊地区和内蒙古阿拉善盟及宁夏贺兰山以西地区［1］。荒漠植被主要分布于干旱区，以旱生、超旱生的小乔木、灌木、半灌木和小半灌木为建群种，是衡量荒漠生态系统健康状况的重要指征［2-3］。

植被地上生物量是指单位面积内所有地上植物的总重量，通常以鲜重或干重表示，其估算方法一直是生态环境监测领域的研究热点［4-5］。传统的植被地上生物量测算主要依靠地面样方调查，虽然精度高，但费时费力、实施困难，对于复杂生态系统存在一定局限，只适用于小尺度研究［6］。与传统方法相比，遥感技术在植被空间异质性识别和不同时相植被信息获取等方面的优势，使其成为大尺度植被地上生物量估算的主要手段［6-9］。

近年来，光学遥感植被信息提取技术被广泛应用于森林和草原植被，但是在荒漠植被中应用较少［10］。这一方面是由于荒漠植被结构简单、低矮稀疏、单位面积地上生物量非常低，传感器对稀疏荒漠植被的敏感性弱，通用遥感模型的估算精度低；另一方面在提取荒漠植被信息过程中容易出现背景（主要是土壤）信息对目标信号的“污染”，进一步降低了遥感模型的估算精度。为了解决植被信息提取精度较低的问题，Cunliffe等［11］在小尺度范围内，利用无人机遥感影像构建改进的多端元光谱分解模型来提取干旱区植被覆盖度，取得了较好的效果。Yan等［12］在毛乌素沙地采用Landsat TM 和MODIS 遥感数据，综合分析了归一化植被指数（NDVI）、比值植被指数（RVI）、差值植被指数（DVI）、土壤调节植被指数（SAVI）、修正的土壤调节植被指数（MSAVI）5 种植被指数与植被覆盖度的相关性，结果显示MSAVI 与稀疏的沙地植被覆盖度的相关性更高。叶静芸等［13］采用QuickBird遥感数据，基于地面调查数据，对乌兰布和沙漠东北缘荒漠-绿洲过渡带的植被地上生物量进行了估算，发现采用RVI线性模型估算荒漠植被地上生物量效果较好。在植被地上生物量模型方面，Zandler等［14］利用多种方法对塔吉克斯坦帕米尔高原东部稀疏荒漠植被的地上生物量估算结果进行比较后发现，套索回归（LASSO）的估算效果最好。综上所述，空间分辨率更高的遥感数据和不断优化的数据分析方法可以有效提高稀疏荒漠植被遥感信息提取的精度。

为了提高稀疏荒漠植被地上生物量的估算精度，本研究以库姆塔格沙漠南部阿尔金山北麓山前戈壁区为研究区，从遥感数据源和模型方法2个方面进行探索和尝试，以建立更加准确、实用的稀疏荒漠植被地上生物量遥感估算方法，以期为干旱区荒漠生态系统碳储量估算和荒漠化监测提供科学依据。

1 研究区概况

研究区位于阿尔金山北麓洪积扇平原，北接库姆塔格沙漠和库鲁克塔格山，总面积为9397.25 km2，地理位置介于39°09′～39°50′N与90°55′～94°28′E之间，海拔1400～2900 m，地势南高北低，东高西低（图1）。根据2017—2020年的观测，研究区年降水量大约在50～180 mm 之间，在空间上由东向西递减，由南向北递减。降水多集中在5—7 月，常常形成洪水，使洪积扇上沟壑纵横。在地势较低的沟床内主要分布有霸王（Zygophyllum xanthoxylon）、裸果木（Gymnocarpos przewalskii）和膜果麻黄（Ephedra przewalskii）等灌木，局部地段分布有零星的盐生草（Halogeton glomeratus）、刺沙蓬（Salsola ruthenica）和雾冰藜（Bassia dasyphylla）等一年生草本植物，群落盖度大多在10%～20%之间，灌木植株高度一般小于1 m，草本植物高度一般不超过10 cm；在地势较高的沟岸地段主要分布着合头藜（Sympegma regelii）、红砂（Reaumuria soongarica）等灌木，群落盖度在10%以下，许多地段甚至在1%以下，植株高度大多为10～50 cm；在洪积扇北部沙漠边缘地带分布有梭梭（Haloxylon ammodendron）、沙拐枣（Calligonum mongolicum）、柽柳（Tamarix ramosissima）等灌木，群落盖度大约10%～20%，植株高度有时可达2.0～2.5 m［15-16］。该地区荒漠植被的植物群落组成、物种多样性和生物量的空间异质性均处于较高水平，对采用遥感方法进行植被生物量的准确估算提出了挑战。

图1 研究区位置以及样方分布Fig.1 Location of the study area and distribution of plots

2 数据与方法

2.1 地面样方调查

2017 年9 月，采用样带与样方相结合的方法进行野外调查。根据阿尔金山北麓戈壁区山前荒漠植被分布的垂直地带性特征，在研究区的东西方向每隔20～50 km 设置一条南北方向样带。沿着样带方向每隔3 km 设置一个30 m×30 m 的样方。调查的12条样带和100个样方的具体分布（图1）。调查内容包括每个样方的地貌特征、土壤特征、海拔和中心点GPS 坐标与植被生长状况，测量并记录样方内每株灌木植物的基径（D，cm）、冠幅长轴（C1，cm）和冠幅短轴（C2，cm）、株高（H，cm），以及每个植株的相对位置。参考《全国生态状况调查评估技术规范——荒漠生态系统野外观测》中的荒漠植被采样方法［17］。在研究区内分不同植物种选取若干株不同大小的健康植株，分别测量其基径、株高、冠幅长轴和冠幅短轴，然后整株贴地剪下，称其鲜重后取样带回实验室，在80 ℃通风干燥箱内烘干48 h称重，经计算得到植株地上生物量干重。少量样方有草本植物分布，在有草本植物分布的样方内，在样方任意1 条对角线的两端和中间各取1 个1 m×1 m的小样方，采用收割法测定小样方内所有草本植物地上生物量鲜重，并采样带回实验室，在80 ℃通风干燥箱内烘干48 h称重，经计算得到每个小样方草本植物地上生物量干重。

2.2 样方地上生物量估算

异速生长模型法是估算灌木地上生物量的有效方法，可以减少对植被的破坏［18］。基于地面样方调查数据，分别以每种灌木的基径（D，cm）、冠面积（A=πC1C2/4，cm2）、株高（H，cm）和植株体积（V=AH，cm3）为预测变量，构建异速生长方程，并选取决定系数（R2）和平均绝对误差（MAE）作为比较和评判异速生长方程拟合效果的指标。计算公式如下：

式中：n为样本个数；yi为单株灌木的生物量实测值（g）；ŷi为单株灌木的生物量估测值（g）；yˉi为每种灌木的生物量实测平均值（g）。

比较估算结果，每种灌木的异速生长方程拟合度均达到了极显著水平。其中泡泡刺（Nitraria sphaerocarpa）、旱蒿（Artemisia xerophytica）、梭梭和柽柳的最佳预测变量为V；沙拐枣、膜果麻黄、裸果木、红砂、合头藜和霸王的最佳预测变量为A，其最优异速生长方程见表1。

表1 灌木植物的最优异速生长方程Tab.1 Best allometric model of shrub plants

根据每个样方调查数据和灌木异速生长方程以及小样方草本植物地上生物量调查数据，对每个样方的地上生物量进行计算，计算结果用于荒漠植被地上生物量遥感估算。

2.3 遥感数据处理

2.3.1Landsat-OLI 遥感数据处理 Landsat-OLI 数据由美国地质勘探局网站下载（https://earthexplorer.usgs.gov/），获取时间是2017 年9 月26 日，覆盖于阿尔金山北麓山前戈壁区，空间分辨率30 m。Landsat-OLI数据包含可见光（VIS）、近红外（NIR）和短红外（SWIR）等9 个光谱波段。Landsat-OLI 数据是经过陆地卫星生态系统干扰自适应处理系统处理后的地表反射率产品，可直接用于遥感特征参数的计算。

2.3.2WorldView-3 遥感数据处理 WorldView-3数据获取时间为2016 年6 月28 日，覆盖面积381 km2，空间分辨率2.5 m，包含蓝光（Blue）、绿光（Green）、红光（RED）和NIR波段。根据课题组在地面样方附近山麓地带设置的自动气象站（距World-View-3 遥感影像覆盖区域西部边缘约10 km）的观测数据，2016年降水量75.5 mm，2017年降水量56.8 mm；该地区气候极端干旱，荒漠植物生长缓慢，1 a内植被生物量变化很小。遥感影像数据获取与地面样方调查之间的时间差异不会给植被地上生物量估算结果带来显著影响。

在ENVI 5.3环境中，利用数据元文件对原始的WorldView-3 数据进行辐射定标；通过高光谱遥感数据的大气纠正模块（FLAASH）对WorldView-3 数据进行大气校正；利用影像的正射校正文件（PRC）和全球定位系统（GPS）差分定位对影像进行几何精校正，误差可以控制在2 m 以内；提取遥感特征参数［19］。

2.3.3遥感特征参数计算 选取单反射率波段和植被指数作为遥感特征参数，用于荒漠植被覆盖区的提取和植被地上生物量回归模型的构建。单反射率波段中RED 和NIR 波段对荒漠植被的敏感性较强［20-21］；植被指数中RVI、NDVI 和MSAVI 对荒漠植被的敏感性更强［21-25］。因此，本研究采用的遥感特征参数有5个，即RED、NIR、RVI、NDVI和MSAVI。

公式（3）、公式（4）和公式（5）为植被指数的计算公式，式中RED和NIR分别表示红光和近红光反射率波段。

2.3.4无人机遥感数据处理 2020 年9 月26—30日，采用DJI 精灵4（RTK）无人机在WorldView-3 数据覆盖区内设置了一条样带进行详细调查，覆盖区域内的荒漠植被具有代表性。该样带南北向长11 km，东西向宽4 km，由44 个1 km×1 km 的飞行单元组成（图1）。为避免太阳入射角度太小而导致的地物阴影拖曳现象，拍摄时间选取在北京时间09:00—16:00 之间。飞行作业期间天气晴好，云量极低，风速很小，满足无人机摄影测量平台工作的要求；飞行高度为200 m，飞行速度9 m·s-1，航向重叠率80%，旁向重叠率75%；空间分辨率0.05 m。数据的平面坐标系统为WGS84坐标系，通用墨卡尔投影。无人机影像覆盖区域的荒漠植被具有代表性，灌木冠幅一般大于20 cm，在无人机影像上可被精准识别。

2.4 植被覆盖区的提取与精度评价

2.4.1植被覆盖区提取原理 利用像元二分法提取植被覆盖区，以提高植被地上生物量模型估算精度［26］。其原理是假设地表信息由纯植被和纯土壤信息构成，则与之相对应的像元光谱信息也由一定覆盖比例的纯植被像元与纯土壤像元构成。表示为：

对公式（6）进行变换，得到植被覆盖度公式：

式中：FVC 为像元中纯植被覆盖面积的比例（植被覆盖度）；S为像元的遥感特征参数值；Sveg和Ssoil分别为像元中纯植被和纯土壤像元的遥感特征参数值。

2.4.2植被和土壤分类 利用目视解译从无人机遥感影像中提取纯土壤像元分布区，其他像元对应的区域则为植被覆盖区［26-31］。通过非监督分类方法对所有无人机影像进行植被和土壤分类，拼接分类结果得到植被和土壤分类图。图2为一个地面调查样方的无人机影像及其分类结果，样方中心地理坐标为39°20′34″N、91°59′4″E，大小为30 m×30 m。从图2 可以看出，无人机影像能够清晰地反映样方内低矮稀疏荒漠植被的冠层信息与纯土壤像元信息。

图2 样方无人机影像及其分类结果Fig.2 Plot UAV image and its classification results

Landsat-OLI和WorldView-3影像的植被覆盖区的提取方法一致。采用非监督分类方法分别对遥感影像进行植被和土壤分类；掩膜提取土壤覆盖区；统计土壤覆盖区遥感特征参数的累计频率。对照无人机影像数据的植被和土壤分类结果并参考相关文献，经过反复调试，选取频率累计表中累计频率为5%的遥感特征参数作为纯土壤像元Ssoil值，遥感特征参数累计频率>5%的像元为植被覆盖区域［26］。

2.4.3植被覆盖区提取与验证 以无人机遥感数据的植被分类结果作为实测值，以遥感影像的分类结果作为预测值，利用Kappa 系数对分类结果进行评价，计算公式如下：

式中：P0为实测值与预测值相同的像元数量；Pe为实测值与预测值不同的像元数量。

分别统计Landsat-OLI 和WorldView-3 遥感数据的5个特征参数影像的土壤覆盖区像元信息。土壤覆盖区内累计频率为5%的特征参数值代表纯土壤像元信息，Landsat-OLI 遥感数据的5 个特征参数对应的Ssoil值分别为0.27、0.31、1.13、0.06 和0.05，WorldView-3遥感数据的5个特征参数对应的Ssoil值分别为0.26、0.22、1.11、0.04 和0.03。累计频率>5%的像元覆盖区域为植被覆盖区。利用无人机遥感数据分别对Landsat-OLI 和WorldView-3 遥感数据的植被和土壤分类结果进行精度评价。由表2 可知，Landsat-OLI的RVI分类效果最优，其Kappa系数为0.77；WorldView-3 的NDVI 分类效果最优，其Kappa 系数为0.95。分别采用最优植被指数RVI 和NDVI 对研究区进行植被覆盖区的提取，World-View-3 影像覆盖区域的植被覆盖区面积为223 km2，Landsat-OLI 影像覆盖区域的植被覆盖区面积为3592 km2。在无人机影像覆盖区随机选取200个像元对植被覆盖区提取结果进行验证，Landsat-OLI影像的植被覆盖区提取精度为84.54%，WorldView-3影像的植被覆盖区提取精度为91.26%。

表2 植被和土壤分类植被指数Tab.2 Vegetation indices for vegetation and soil classification

2.5 植被地上生物量遥感模型的构建与验证

将研究区内100个样方的中心点GPS坐标叠加于Landsat-OLI影像上，从中随机选取70%（n=70）的样方用于模型构建，其余30%（n=30）的样方用于模型验证。建立70 个样方的植被地上生物量数据与Landsat-OLI 遥感影像上对应像元的光谱特征值之间的回归模型。100 个样方中有21 个样方分布于WorldView-3 影像覆盖区，将这21 个样方均分为84个15 m×15 m的小样方。样方内所有植株的相对位置已知，可以准确计算每个小样方的植被地上生物量，将其中心点的GPS 坐标叠加于WorldView-3 影像上，随机选取70%（n=58）的小样方，建立小样方地上生物量数据与WorldView-3遥感特征参数之间的回归模型，剩余30%的小样方（n=26）用于回归模型验证。Townshend 等［30］研究表明，遥感数据验证的最佳样方大小为E=F（1+2L），其中E为样方大小（m2），F为像元大小（m），L为像元几何精度（m）。WorldView-3 遥感数据的像元大小为2.5 m，几何精度2 m 以内，则最佳样方大小为12.5 m，其与15 m×15 m的小样方大小相近。因此，样方数据与遥感像元之间存在匹配关系。

SLR、LASSO、RR 模型常被应用于植被地上生物量的预测［32-33］。在R 语言环境中构建回归模型，利用“glmnet”包构建遥感特征参数与植被地上生物量之间的回归模型，并利用“sperrorest”包进行模型的交叉验证［34］。交叉验证数据集被随机分配至G（G=10）组，交叉验证重复10 次，选取均方根误差（RMSE）、标准偏差（SD）和均方根误差百分比（RMSErel）作为模型精度的指标，用于对比分析模型的回归效果。RMSE用来衡量预测值与实测值之间的偏差；SD为各数据偏离平均数的距离；RMSErel为均方根的偏移。3项指标越小表示模型的误差越小。

式中：yi为第i个样方植被地上生物量实测值（g·m-2）；ŷi为第i个样方植被地上生物量预测值（g·m-2）；yˉ为实测植被地上生物量的平均值（g·m-2）；n为模型中使用的样方个数。

2.6 植被地上生物量遥感估算结果的修正

构建WorldView-3 和Landsat-OLI 影像植被地上生物量分布数据之间的线性关系，通过该线性关系对基于Landsat-OLI 影像得到的阿尔金山北麓山前戈壁区的荒漠植被地上生物量分布图进行修正。

3 结果与分析

3.1 估算结果验证

采用交叉验证方法对所有样方的估算结果进行对比分析（表3）。WorldView-3 数据的反演效果明显优于Landsat-OLI，RMSE平均降低58.65 g·m-2，RMSErel平均降低了72.27%。LASSO、RR和SLR模型中，RR 模型的估算效果最好，其中WorldView-3 RR 模型的RMSErel 为44.68%，略优于WorldView-3 SLR 模型的58.33%和WorldView-3 LASSO 模型的45.26%；Landsat-OLI RR模型的RMSErel为114.34%，略优于Landsat-OLI SLR 模型的134.94%和Landsat-OLI LASSO模型的115.81%。

表3 植被地上生物量预测结果交叉验证Tab.3 Result of cross-validated aboveground biomass in the samples

3.2 估算结果分析

利 用WorldView-3 RR 和Landsat-OLI RR 模型估算植被地上生物量，得到模型预测值与实测值的1:1 线性图（图3）。利用WorldView-3 RR 模型反演得到的植被地上生物量预测值与样方实测值更贴近于1:1线，基于WorldView-3数据的估算效果明显优于Landsat-OLI RR。Landsat-OLI RR 模型植被地上生物量估算值基本上不随实测值的变化而发生变化，表明Landsat-OLI遥感数据对于稀疏荒漠植被的敏感性非常低，在实际应用中会使植被地上生物量估算结果出现较大偏差。

图3 植被地上生物量实测值和预测值的1:1线性图Fig.3 Scatterplots of measured and predicted aboveground biomass

在WorldView-3 数据覆盖区内随机选取500 个像元，每个像元之间的距离不小于100 m，构建WorldView-3 RR 和Landsat-OLI RR 植被地上生物量估算数据之间的线性关系，对Landsat-OLI RR 植被地上生物量估算结果进行修正。WorldView-3 RR与Landsat-OLI RR估算值之间的关系式为（图4）：

图4 植被地上生物量回归分析Fig.4 Regression analysis of aboveground biomass

式中：OLI 与WV 分别为Landsat-OLI RR 与World-View-3 RR模型的植被地上生物量预测值（g·m-2）。

3.3 植被地上生物量遥感估算

利用WorldView-3 RR 估算WorldView-3 影像覆盖区的植被地上生物量，其均值为122.43 g·m-2（RMSErel=44.68%）（图5a）；利用Landsat-OLI RR 估算WorldView-3 影像覆盖区的植被地上生物量，其均值为84.59 g·m-2（RMSErel=114.34%）（图5b）。Landsat-OLI RR的RMSErel比WorldView-3 RR高出69.66%。利用公式（12）对Landsat-OLI RR 模型估算结果修正，修正的植被地上生物量均值为90.35 g·m-2（RMSErel=83.16%）（图5c）。修正方程在一定程度上提高了Landsat-OLI RR 模型的估算精度，RMSErel降低了31.18%。对比图5b和图5c可见，修正后的植被地上生物量分布更加“精细”。洪积扇戈壁表面因洪水而沟壑纵横，受微地形影响植被分布空间异质性较高，图5c对植被地上生物量分布的空间异质性刻画得更加清晰。

图5 WorldView-3遥感影像覆盖区域植被地上生物量估算结果Fig.5 Estimation results of aboveground biomass in the area covered by WorldView-3 remote sensing data

利用Landsat-OLI RR模型估算研究区植被地上生物量，其均值为85.28 g·m-2（RMSErel=84.93%）（图6a）。修正后的植被生物量均值为118.28 g·m-2（RMSErel=71.51%）（图6b）。从植被地上生物量分布来看，图6b的植被地上生物量分布格局与地面调查数据更为符合，南部海拔较高的山麓地区植被地上生物量较高，而北部接近于沙漠边缘的海拔较低区域植被地上生物量较低。对比两幅地上生物量分布图可以发现，在植被地上生物量较低的93°00′E 附近地区和植被地上生物量普遍较高的94°20′E 附近地区，修正前后分别出现了高值变低和低值变高的现象，表明低植被地上生物量区域被高估和高植被地上生物量区域被低估的问题得到了一定程度的修正。

图6 植被地上生物量估算结果Fig.6 Estimation results of aboveground biomass

4 结论

（1）遥感数据空间分辨率是影响稀疏荒漠植被地上生物量回归模型估算精度的主要因素。与中空间分辨率Landsat-OLI遥感数据相比，基于高空间分辨率WorldView-3遥感数据建立的稀疏荒漠植被地上生物量估算模型估算精度更高。

（2）传统的植被地上生物量回归模型多是利用植被指数等遥感特征参数与地面数据的相关关系直接构建。虽然采用该方法可能会得到较高的估算精度，但是遥感特征参数之间的共线问题会使模型的稳定性降低。采用非线性的多元回归方法SLR、LASSO和RR构建遥感特征参数与地面样方数据之间的回归模型可以有效提高模型的稳定性。

（3）提取稀疏荒漠植被信息时，采用高空间分辨率遥感数据的估算结果对中空间分辨率遥感数据的估算结果进行修正，可以有效提高中空间分辨率遥感数据对稀疏荒漠植被信息的提取精度，并可以在具有相同自然条件的较大尺度上应用。