基于多光谱数据与小型无人机的甘南草地非生长季植被覆盖度

冷若琳，张瑶瑶，谢建全，李芙凝，胥 刚，崔 霞

（1. 兰州大学资源环境学院西部环境教育部重点实验室，甘肃 兰州 730000；2. 草地农业生态系统国家重点实验室 /兰州大学草地农业科技学院，甘肃 兰州 730020）

草地生态系统是除森林、苔原和农田生态系统之外的最主要陆地生态系统，是重要的畜牧业生产基地和生态屏障，为家养牲畜和野生动物提供了栖息地[1]。草地植被覆盖度(fractional vegetation cover，FVC)能够定量描述地表草类植物覆盖情况及生长状况，在植被蒸腾、光合作用、土壤水分蒸发、水土流失等变化中都是一个重要的影响因子[2]。

估测草地植被覆盖度的传统方法是地面测量法，但是由于不同区域空间异质性的影响，很难对大范围的草地植被覆盖度进行精准估测[3]。遥感技术的发展为草地植被覆盖度的估测提供了新的技术手段，尤其为大范围快速和准确估测提供了可能。近几十年来，开展了大量的基于遥感技术估算草地植被覆盖度的研究，但这些研究大多是关注生长季绿色植被覆盖度，对草地非生长季植被覆盖度的遥感估算研究甚少。

草地非生长季植被是指在非生长季(11月到次年4月之间)出现的不能进行光合作用的草地植被。甘南州草地超载过牧现象较为严重，尤其是在冬春季草畜矛盾尤为突出[4]，在整个非生长季，牲畜都要以冬季牧场饲草作为主要的饲料来源，很多地区家畜出现“夏壮、秋肥、冬弱、春死亡”的现象，冬春季节严重缺草造成牲畜掉膘和死亡是制约牧区草地畜牧业发展的瓶颈，在确定牧区载畜量时，不仅要依据饲草总产量，同时也要考虑季节变化的影响。监测草地非生长季植被的空间分布状况与时间动态对冬季饲草贮备、计算牧区季节载畜量等问题都有重要意义，可以为当地畜牧业发展提供科学依据，有利于甘南州畜牧业可持续发展。

由于草地非生长季植被中所含纤维素、半纤维素以及木质素等非结构组分在1.7、2.1以及2.35 μm波段附近的光谱吸收相对深度与裸土不同[5]，可为草地非生长季植被覆盖度遥感估算提供理论基础。由于高光谱数据含有丰富的地物光谱信息，Daughtry和Craig[6]基于高光谱数据建立的纤维素吸收指数(cellulose absorption index，CAI)有效估算了农作物秸秆残茬覆盖度，目前CAI是估算非光合植被覆盖度的最佳指标。但是计算CAI指数所需要的高光谱数据目前只来源于少数卫星(如EO-1 Hyperion)、机载成像光谱仪(airborne visible infrared imaging spectrometer，AVIRIS)以及地物光谱仪，可获得性差。Serbin等[7]利用ASTER数据建立了短波红外归一化残留指数(the shortwave infrared normalized difference residue index，SINDIR)，可 以代替CAI进行草地非生长季植被覆盖度监测，但是由于ASTER数据时间分辨率较低，且每景数据覆盖面积小，并不适合区域范围的植被监测。

中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer，MODIS)具有大空间性、高时效性以及低成本等特点[8]，可用于草地非生长季植被覆盖度研究。长久以来学者们一直试图挖掘多光谱遥感数据(如TM、MODIS等)在草地非生长季植被覆盖度监测方面的潜力。Van Deventer等[9]首次提出利用TM数据5和7波段组成土壤耕作指数(soil tillage index，STI)和 归 一 化 差 异 耕 作 指 数(normalized difference tillage index，NDTI)来对土壤特性进行分类，十几年后Kergoat等[8]将STI指数与草地非生长季植被覆盖度以及生物量联系起来，并证明TM短波红外(shortwave infrared，SWIR) 1.6与2.1 μm波段在监测非光合植被方面有很大优势。任鸿瑞等[10]基于MODIS的SWIR波段建立了与CAI高度相关的NDTI指数，证明SWIR两个波段组合成的指数可以替代CAI进行草地非生长季植被生物量监测。曹鑫等[11]在STI的基础上，基于MODIS数据中绿色植被-非光合植被-土壤3种组分混合物的光谱特性建立了干枯燃料指数(dead fuel index，DFI)，并证明其对与绿色植被混合的草地非生长季植被覆盖度具有较好的估测能力。

由于SWIR1、2组成的指数可以替代CAI，有效估算草地非生长季植被覆盖度[8,10,12]，那么众多植被指数中哪一种最适于非光合植被覆盖度的监测？对于常用的两种多光谱数据源(MODIS、Landsat)，在甘南草地非生长季植被覆盖度监测中的适用性如何？

为此，本研究比较了基于不同空间分辨率遥感数据(Landsat8 OLI、MODIS)的9种植被指数估算甘南州草地非生长季植被覆盖度的能力，分析草地实测数据与多种植被指数之间的相关关系，根据精度验证结果筛选出甘南草地非生长季植被覆盖度最优反演模型，并分析甘南州2018年草地4月 -5月植被覆盖度空间分布特征，以期为甘南州草地非生长季载畜量计算和植被动态监测提供依据。

1 材料与方法

1.1 研究区概况

甘南藏族自治州地处甘肃省西南部，位于黄土高原的西部和青藏高原东北边缘的过渡地带，地理坐标100°46′-104°44′ E，33°06′-36°10′ N[3-4,13-14]。州内海拔1 000～4 900 m，多数地区大于3 000 m，大部分区域属于高原大陆性气候，年均降水量400～800 mm，年均温度1～3 ℃，年均日照时长2 200～2 400 h，自西北向东南逐渐降低[4,15]。甘南州南部为岷迭山区，气候温和，是全省重要的林区之一；东部为丘陵山地，高寒阴湿；西北部为广阔的草甸草原，牧草丰茂，草地面积约为2.063 × 106hm2，占甘南州土地总面积的67.64%，其中可以利用的 草 地 面 积 为2.495 × 106hm2，占 草 地 总 面 积 的95.86%[13-16]。高寒草甸、高寒灌丛草甸、温性草甸草原、温性草原、沼泽、暖性草丛、低平地草甸是甘南州的主要草地类型，玛曲、碌曲、夏河、合作4个县(市)为甘南州的主要牧区[15]。

1.2 数据获取

1.2.1 地面实测数据获取

野外实测数据采集于2018年4月25日-29日进行，在甘南州主要牧区(玛曲县、碌曲县、夏河县、合作市)选择地势较为平坦、具有代表性的典型群落地段设置样地，共获得13个样地数据和65个样方数据(图1)，样地大小为60 m × 60 m，在每个样地四角及中心点设置一个样方，利用大疆精灵无人机Phantom 3 Professional设置航线垂直对地进行拍摄，获取每个样方的照片(图2)，每个样地获取5个样方照片，无人机飞行高度为30 m[15]。每个样方的经纬度从无人机拍摄的照片上获取，样地的位置为中心样方经纬度。

图 1 甘南州样地分布图Figure 1 Distribution map of the sampled land in Gannan

1.2.2 遥感数据获取

图 2 样地和样方设置图(黑色正方形为样方)Figure 2 Plot and quadrat design (The black square is the quadrat)

MODIS数据：采用Terra和Aqua卫星结合的中分辨率成像光谱仪MODIS双向反射分布函数和地表反射率产品MCD43A4[11,16]，该产品经过大气分子散射、臭氧、气溶胶及云处理，对MODIS 7个波段的数据均进行了大气双向反射和太阳高度矫正 以 及NBAR(Nadir BRDF Adjusted Reflectance)纠正，时间分辨率为1 d，空间分辨率为500 m[16]。该数据来自美国国家航空航天局(National Aeronautics and Space Administration，NASA)对 地 观 测 系 统 数据共享平台(EOS data gateway)，获取从2018年4月25日-29日共5 d数据。

Landsat8 OLI数据：空间分辨率为30 m，时间分辨率为16 d，本研究中获取与实测数据获取时间相近的覆盖研究区的2018年4月-5月的4期云量较少的影像(表1)。该数据在美国地质调查局(USGS，United States Geological Survey)官方网站下载[5]。

表 1 研究区Landsat遥感数据Table 1 Landsat data list of the research area

2 研究方法

2.1 非生长季植被覆盖度信息提取

使用ENVI 5.3中监督分类的方法，对无人机所拍摄到的每个样方照片进行分类并提取非生长季植被覆盖度数据。对于Landsat影像来说，由于其空间分辨率为30 m，因而每个样方的大小都能与之匹配，共计获取65个样方的非生长季植被覆盖度数据。对于MODIS数据需要将每个样地内5个样方的覆盖度取平均值作为该样地非生长季植被覆盖度值与MODIS影像相匹配[5,15]，共计获取13个样地的非生长季植被覆盖度数据。

2.2 遥感影像处理

MODIS影像数据用MRT(modis projection tool)软件进行读取[3]，并提取出其中7个波段的数据，同时转换遥感影像的数据格式和地图投影，最后利用ArcGIS提取样地的波段数据以及质量数据，结合质量数据，对样地所对应的像元反射率数据质量进行评价，剔除受云和气溶胶影响较大的数据。

Landsat8 OLI影像数据用ENVI 5.3的Radiometric Calibration模块和FLAASH Atmosphere Correction模块分别对影像数据进行辐射定标与大气校正处理，并提取出波段1～7的TIFF格式文件，在ENVI中进行波段计算，最终提取每个样地所对应像元的相应反射率数据[15]。

2.3 植被指数计算

基于非光合植被在SWIR1、SWIR2波段反射率差异构建的植被指数在非生长季植被覆盖度估算中得到了较为广泛的应用[10-12]。同样，也有研究表明，非光合植被在可见光波段的光谱差异也为可见光植被指数的构建创造了可能[17-19]。本研究总结了目前所提出的可用于草地非生长季植被覆盖度估算的植被指数(表2)。根据提取的Landsat8 OLI、MODIS数据的对应波段数据，根据表2中的公式计算得到基于不同数据源的9种植被指数。

2.4 非生长季植被覆盖度反演模型及精度评价

将9种植被指数作为自变量，基于无人机图像的非生长季植被覆盖度数据作为因变量，利用SPSS软件分别计算各指数与非生长季植被覆盖度之间线性、对数、乘幂以及指数回归模型的参数。由于本研究中获取的样地数量较少，因而采用适用于样地数量少的留一法交叉验证法(leaveone-out cross-validation，LOOCV)，对预测误差进行无偏差估计[25]。这种验证方法每次从总样本里面挑出一个样本作为验证集，余下的n-1个样本建立回归模型，用挑出的样本验证模型精度，进行n次运算后，得到全部样本的RMSEP(root-mean-square error of prediction) (式1)和观测值与模型预测结果之间的交叉系数(cross-validation r)。

表 2 用于非生长季植被覆盖度估算的几种主要植被指数Table 2 Several major vegetation indexes (VIs) used in estimating vegetation coverages in non-growing seasons

2.5 草地非生长季植被覆盖度空间分布

利用最终确定的遥感监测模型以及遥感数据，反演2018年4月-5月甘南州草地非生长季植被覆盖度，并分析其空间分布特征。

3 结果与分析

3.1 不同植被指数与草地非生长季植被覆盖度相关性分析

基于Landsat8 OLI与MODIS两种数据源计算不同植被指数，对比9个指数与草地非生长季植被覆盖度之间的相关系数r (表3)。总体来说基于Landsat8的植被指数与草地非生长季覆盖度相关性较高，当显著水平为0.01时，与覆盖度达到极显著的指数分别为NDTI、STI、Ratio、GRCI、BRCI，r介于0.569～0.638；当显著水平为0.05时，与覆盖度相关系数达到显著的指数为NDVI、NDSVI、GNDVI；基于MODIS的植被指数与草地非生长季植被覆盖度相关系数较低，NDTI、Ratio、STI、DFI这4个指数与覆盖度显著相关(P ＜ 0.05)，相关系数介于0.512～0.565，其中DFI与覆盖度的相关性最佳，相关系数r达到0.565，其余指数与覆盖度均无明显关系。对比基于不同数据源的9种植被指数与草地非生长季植被覆盖度之间的相关系数，在构建反演模型时选取了相关性均较为显著的两个指数：NDTI和STI (因Ratio与STI互为倒数，故本研究仅计算了STI)，同时由于DFI是专为MODIS数据设计的，且本研究基于MODIS的DFI与植被覆盖度有显著相关性，因此选取NDTI、STI和DFI这3种植被指数分别建立草地非生长季植被覆盖度反演模型。

表 3 草地覆盖度与植被指数之间的相关性比较Table 3 Comparison of correlation coefficients of FVC with VIs

3.2 基于不同数据源的反演模型对比

3.2.1 基于Landsat数据的非生长季植被覆盖度反演模型构建及精度验证

在基于Landsat 8 OLI数据的植被指数与覆盖度进行回归分析时，发现以植被指数为纵坐标、覆盖度为横坐标的散点图中，非生长季植被覆盖度相同或相近的位置，所对应的植被指数却有所差别，采用Xu等[5]的方法将覆盖度相同或相近点(相差5%以内)的植被指数进行平均(图3)，可以将其他因素引起的植被指数差异最小化。

图 3 覆盖度相近的点所对应的3个植被指数平均前后对比Figure 3 Comparison of averaged before and after points with similar coverage

基于Landsat8 OLI的3种植被指数反演草地非生长季植被覆盖度的模型方程及精度评价结果如表4所示。其中NDTI与覆盖度呈极显著相关，其线性模型的R2最高达到0.407，高于对数、乘幂以及指数模型的决定系数；基于NDTI的指数模型具有最小的RMSEP (12.805)及最高的交叉系数(0.672)，但是线性模型的R2远大于指数模型(分别为0.407和0.352)。因此，甘南州基于Landsat8 OLI的草地非生长季植被覆盖度最优模型为基于NDTI的线性模型(R2= 0.407，RMSEP = 12.773，r = 0.638)：

式中：y表示基于Landsat 8 OLI NDTI反演的非生长季植被覆盖度(%)；x表示基于Landsat 8 OLI的植被指数NDTI。

3.2.2 基于MODIS数据的非生长季植被覆盖度反演模型构建及不同数据源结果对比

与Landsat8 OLI数据经过相同方法处理，建立基于MODIS的3种植被指数与草地非生长季植被覆盖度的不同回归模型如表4所示。其中基于DFI的反演模型在R2以及精度评价方面都明显优于STI与NDTI，基于DFI的指数模型R2和交叉系数最高，分别达到0.319和0.544；DFI线性模型RMSEP最小，为4.310。虽然DFI指数模型的RMSEP略高于线性模型，但是R2和交叉系数却比后者大很多。因而，基于MODIS的草地非生长季植被覆盖度最优反演模型为指数模型：

式中：y表示基于MODIS DFI反演的非生长季植被覆盖度(%)；x表示基于MODIS的植被指数DFI。

表 4 基于Landsat 8 OLI与MODIS的植被指数与覆盖度反演模型及精度评价对比Table 4 Comparison of the inversion model and accuracy evaluation between VIs and FVC based on Landsa8 OLI and MODIS

对比基于Landsat 8 OLI与MODIS的植被指数与非生长季植被覆盖度最优反演模型(分别为Landsat NDTI线性模型与MODIS DFI指数模型)，发现Landsat NDTI与覆盖度反演模型的相关系数(R2= 0.407)及交叉系数(r = 0.638)均大于后者(R2=0.319，r = 0.544)。MODIS DFI指数模型的RMSEP较小(RMSEP = 4.321)，但是Landsat NDTI与覆盖度呈极显著相关(MODIS DFI与覆盖度仅为显著相关)，且相关系数远大于MODIS DFI。综上所述，甘南州非生长季植被覆盖度最优反演模型为Landsat NDTI线性模型(式2)。

3.3 甘南州2018年4月- 5月植被覆盖度空间分布

利用甘南州非生长季植被覆盖度最优反演模型，得到甘南州2018年4月-5月草地植被覆盖度空间分布图(图4)。甘南州非生长季草地植被集中分布在西部合作、夏河、碌曲、玛曲4个牧业县中，东部草地植被覆盖较少。总体上甘南州植被覆盖度西部高、东部低，其中玛曲县草地植被覆盖度最高，其西北部覆盖度维持在60%以上；合作、夏河、碌曲大部分地区草地覆盖度在20%～50%，夏河北部、合作中部以及玛曲东南部少部分区域植被覆盖度小于20%。

图 4 甘南州地区2018年4月-5月非生长季的草地植被覆盖度空间分布图Figure 4 The spatial distribution map of grassland FVC in non-growing seasons from April to May, 2018, in Gannan

4 讨论

目前对绿色植被覆盖度的研究已经较为成熟，NDVI、EVI等植被指数经过多年验证，是草地绿色植被生长状态以及覆盖度的最佳指示因子[27]。但是基于多光谱数据估测草地非生长季植被还处于起步阶段，尚未找出具有普适性的最优监测指数。与绿色植被覆盖度反演模型相比，草地非生长季植被覆盖度反演模型的R2要低一些[3,13,15]。这是由于高寒草地在冬春季节雨雪充沛，降雪过后对土壤湿度有很大影响，土壤湿度是影响非生长季植被监测的一个重要因素[27]。

本研究选取了NDTI、STI、DFI这3个指数分别建立植被覆盖度监测模型。其中STI是由短波红外两波段简单比值组合成的[26]，短波红外波段光谱带位于电磁波谱的吸水区附近，可用于增强枯黄植被的信号，当植被由绿色变为枯黄时，这些区域的光谱响应会由于叶组织中的水分流失而增加[21]，因此STI对草地非生长季植被覆盖度有很好的响应。NDTI在STI的基础上，采用归一化算法能够在一定程度上抵消大气影响[9]，与本研究中基于Landsat8 OLI的NDTI反演模型精度略高于STI的结论一致。Cao等[11]通过分析绿色植被、非光合植被和裸土平均光谱曲线，根据其在MODIS波段范围内的光谱特征提出了DFI指数。基于MODIS的DFI与草地非生长季植被覆盖度具有良好的相关性，但是基于Landsat8 OLI的DFI与覆盖度之间却没有这样的规律，其根本原因是MODIS与Landsat8 OLI数据的波段设置存在差异。

利用Landsat8 OLI以及MODIS数据估算草地非生长季植被覆盖度，充分发挥了多光谱遥感数据低成本、大范围、长时间观测的优势[23,28-29]。本研究表明，SWIR1、SWIR2两波段对草地非生长季植被的监测能力优于可见光三波段(R、G、B)以及近红外波段，反演草地非生长季植被覆盖度的最优指数是基于Landsat8 OLI的NDTI。

此外，由于研究区设置的样方大小(30 m × 30 m)与Landsat8 OLI的空间分辨率更为匹配，MODIS数据空间异质性较大，60 m × 60 m的样方不能很好的与其匹配，导致本研究中基于MODIS的植被指数整体上反演草地非生长季植被覆盖度的能力与Landsat8 OLI相比较差。虽然地面样方法与传统测量方法具有高度相关性[3]，但是如何设置地面样方使其与低分辨率遥感数据(500 m × 500 m)相匹配，提高地面样方的代表性，是今后研究中需要解决的关键问题。目前对非生长季植被监测存在模型适用范围小、精度差等困难，由于不同植被类型、土壤湿度与类型以及非光合植被随时间分解的影响，在短波红外波段会产生光谱差异[28,30]。高海拔地区在非生长季期间的遥感影像会受到大气、云以及雪的影响，多光谱遥感数据与野外实测数据在时间以及空间上不能完全对应，这些都是草地非生长季植被覆盖度反演误差的主要来源[27]。

5 结论

本研究选取了甘南州研究区2018年4月-5月野外调查的草地数据，结合该时段MODIS以及Landsat8 OLI遥感数据，对比了基于不同遥感数据的9种植被指数估算草地非生长季植被覆盖度的能力，并建立了甘南州地区基于MODIS DFI以及Landsat NDTI的覆盖度反演模型。在此基础上，通过精度分析确立了用于估算草地非生长季植被覆盖度的最优反演模型，得出：

1)与MODIS数据相比，Landsat8 OLI空间分辨率更高，与多种植被指数呈显著性相关，更适合用于甘南州草地非生长季植被覆盖度估算。

2)基于Landsat8 OLI的NDTI指数对研究区草地非生长季植被覆盖度反应最为敏感，NDTI与覆盖度之间的线性模型为研究区草地非生长季植被覆盖度的最优反演模型(y = 1 432.074x-166.855)，决定系数R2为0.407，RMSEP为12.773。该模型可用于研究甘南州非生长季植被覆盖度。

3)甘南州2018年4月-5月草地植被覆盖度总体上西部高、东部低，其中玛曲县草地植被覆盖度最高，其西北部覆盖度维持在60%以上；合作、夏河、碌曲大部分地区草地覆盖度在20%～50%，夏河北部、合作中部以及玛曲东南部少部分区域植被覆盖度小于20%。