高寒草地生长季/非生长季植被盖度遥感反演

2021-09-22刘佳丽范建容张茜彧杨超徐富宝张晓雪梁博

草业学报 2021年9期

刘佳丽，范建容，张茜彧，杨超，徐富宝，张晓雪，梁博

（1.中国科学院水利部成都山地灾害与环境研究所，四川成都610041；2.中国科学院大学，北京100049；3.西藏自治区水土保持局，西藏拉萨850000）

草地生态系统作为陆地生态系统的重要组成部分，是人类生存重要的自然资源和环境基础要素。青藏高原高寒草地是中国乃至亚洲的重要牧区之一，占青藏高原总面积的59.28%，其生态系统的变化对我国生态安全屏障功能具有重要意义［1］。植被盖度（fractional vegetation cover，FVC）指观测区域内植被投影面积占地表的百分比，是定量描述植被生长状况的重要参数［2］。高寒草地植被盖度变化在截流降雨、防沙治沙、保水固土、调节气候等生态应用方面起着重要作用。因此，高寒草地植被盖度研究对评价区域生态环境具有重要意义。

根据植被的功能属性可以将其分为光合植被（photosynthetic vegetation，PV）和非光合植被（nonphotosynthetic vegetation，NPV）［3］。近年来，开展了大量估算植被盖度的研究［4］，但这些研究多关注的是生长季绿色植被盖度，对非生长季植被盖度的研究很少。进行光合作用的绿色植被无疑是植被重要的组成部分，但并非唯一要素。草地的非生长季通常占全年3/4以上，枯草、枯叶、枯枝、枯干和凋落物等在草地生态系统中发挥着重要作用［5］。NPV作为冬春季牧场主要的牧场饲料，对畜牧业的发展具有重要意义［6］。除此之外，NPV不仅影响着植被生态系统的CO2交换量和碳存储，还能够通过影响土壤（bare soil，BS）与大气间的水热传输条件，进而影响净初级生产力［7］。加之，单一的绿色植被盖度指标无法反映真实的群落结构对植被防蚀功效的内在影响，考虑NPV的影响对减缓径流、减少土壤侵蚀都具有积极的作用［8］。因此，精确估算生长季/非生长季草地植被盖度具有重要意义。

利用像元二分模型提取高寒草原植被盖度具有可行性［9］。大量研究表明，归一化植被指数（normalized difference vegetation index，NDVI）作为像元二分模型的输入变量时，能反映植被的生长分布特征及变化情况，与光合植被盖度有显著的相关性［10－11］。但NDVI难以反映非生长季枯黄植被的盖度信息［12］。在植被光谱响应中，与PV明显的红边位置相比，NPV、BS在可见光－近红外波段（400～1100 nm）区间内具有相似的反射率曲线，因此，无法利用可见光－近红外波段对PV和NPV两者进行区分［13］。但在短波红外（2100 nm）附近波段，NPV具有明显的吸收特征，而BS和PV的光谱曲线中则没有这种吸收特征［14］，这为草地非生长季植被盖度遥感估算提供理论基础。基于高光谱数据的纤维素吸收指数（cellulose absorption index，CAI）是估算NPV的有效指标［15］，但是高光谱图像难以获取，不能满足大面积内估算盖度的需求。相比之下，整个地球表面长期以来已收集了丰富的多光谱数据，如何利用这些数据的时间和空间优势来估算大范围内的PV和NPV变得尤为重要。一些基于多光谱的植被指数被用来估算NPV，例如基于Landsat TM的归一化差异耕作指数（normalized difference tillage index，NDTI）［16］，改进土壤调节作物残茬指数（modified soil-adjusted corn residue index，MSACRI）［17］，基于先进星载热发射和反射辐射计（advanced spaceborne thermal emission and reflection radiometer，ASTER）的短波红外归一化差异残渣指数（the shortwave infrared normalized difference residue index，SINDRI）［18］和基于MODIS的干枯燃料指数（dead fuel index，DFI）［19］。尽管在不同环境中使用这些植被指数估算NPV具有一定可行性，但是这种方法是针对特定区域的，在其他环境中并未得到很好的验证。光谱混合分析（spectral mixture analysis，SMA）提供了另一种有效方法来估算PV、NPV［20］，但该方法的准确性在很大程度上取决于端元成员的选择。Daughtry等［21］研究PV、NPV和BS的NDVI和CAI时发现，NDVI－CAI特征空间近似三角形，在端元已知的情况下可以通过线性模型来分解混合像元，Guerschman等［3］据此首次提出运用NDVI－CAI像元三分模型估算植被盖度的新方法。事实上，像元三分模型中的NDVI和CAI指数可被其他植被指数替代。

目前大多数的植被盖度研究中，相对于光合植被参数的定量反演研究，针对非光合植被参数的定量反演研究较少，仅使用单一的植被盖度和植被指数关系进行估算，这给植被盖度的季节估算带来了严重的不确定性。加之，用于植被盖度季节建模的植被指数种类繁多，不同植被指数对植被盖度的表征能力及敏感性存在一定的差异。面对如此多的植被指数，如何筛选出最优的植被指数，进一步提高植被盖度季节动态监测的精度仍是一项重要研究内容。因此，本研究以西藏当雄县为例，将Landsat-8OLI数据作为遥感信息源，筛选适用于估算高寒草地植被盖度的最优植被指数，并开展了高寒草地生长季/非生长季植被盖度反演研究。

1 材料与方法

1.1 研究区概况

当雄县位于西藏自治区中部，藏南与藏北的交界地带。地理范围为90°45′－91°31′E，29°31′－31°04′N（图1），境内平均海拔4300 m，总地势由西北向东南倾斜，念青唐古拉山脉沿县西北贯穿全境，东北部为高平原，北部高平原分布着纳木错湖，西北和东南皆为山地，中部间夹宽谷盆地。当雄县属高原亚寒带半干旱季风气候，年均温为1.3℃，年均降水量为456.8 mm，降水主要集中在6－8月，年均蒸发量为1725.7 mm。当雄县是寒冷半湿润高寒草甸与寒冷半干旱高寒草原类型之间的过渡地带。植被类型以高寒草甸为主，较高海拔（4700～5200 m）是以高山嵩草（Kobresia pygmaea）和垫状点地梅（Androsace tapete）为优势种的高寒草甸。较低海拔（4300～4650 m）是以丝颖针茅（Stipa capillacea）和矮生嵩草（Kobresia humilis）为优势种的草原化草甸，在海拔4300 m以下地带主要是由河流与低洼地形成的沼泽湿地［22］。草场占当雄县土地总面积的68.86%，全县牧草生长期90～120 d，生长季为5－9月，其余时期为非生长期。当雄县的土壤类型主要为黑褐色酸性至中性的高山草甸土，分布于较高海拔，且含水量较高。而在较低海拔地区，土壤类型则为中性到微碱性的草原化草甸土，砾石含量较高［23］。

图1 当雄县样地分布Fig.1 Distribution of sample plots in Dangxiong County

1.2 数据获取

1.2.1 地面实测数据获取地面实测数据采集时间为2019年4月16－19日（非生长季）以及2019年8月21－27日（生长季）。选择地势平坦、生长状况均一、可达性高的典型草地植物群落布设样地，样地大小为90 m×90 m。4月共布设108个样地，8月共布设109个样地（图1）。采用网格法测定植被盖度，在样地内布设3～5个样方，大小为1 m×1 m，将同一样地内所有样方植被盖度的平均值作为该样地的植被盖度。同时记录样方名称、经纬度、高程、草高、砾石盖度、植被种类、照片编号等信息。为更好地计算植被指数，本研究使用SVCHR-1024地物光谱仪测定土壤高光谱数据。测量时间选择晴朗无云无风、太阳光强度稳定的时段（10：00－14：00）。测量前，先对照白板进行定标，再选择31个地面不同盖度的样方，将25°视场角的传感器探头垂直置于样点上方1.35 m处，探头接收垂直中心点直径0.6 m圆形范围内的光谱。每个样点重复测量5次，后期对光谱曲线进行去重叠和平滑处理，将5次测量的光谱反射率的平均值作为每个样点的反射光谱数据。

1.2.2 遥感数据获取以30 m空间分辨率的Landsat-8OLI地表反射率产品为数据源，影像时间范围选取为2019年4月1－30日以及2019年8月1－31日。对Landsat-8OLI地表反射率产品进行质量筛选、辐射校正、大气校正、几何校正、去除云污染、影像拼接等预处理，获得当雄县2019年4、8月影像。中国省级行政区划矢量数据来源于中国科学院资源环境科学数据中心（http：//www.resdc.cn）。

1.3 研究方法

1.3.1 植被指数计算结合前人的研究结果及本研究实际情况，共选取NDVI［24］、比值植被指数（ratio vegetation index，RVI）［25］、土壤调整植被指数（soiladjusted vegetation index，SAVI）［26］、转换型土壤调整植被指数（transformed soil-adjusted vegetation index，TSAVI）［27］、修改型土壤调整植被指数（modified soiladjusted vegetation index，MSAVI）［28］、NDTI［16］、MSACRI［17］、土壤调整型归一化差异耕作指数（soiladjusted normalized difference tillage index，S-NDTI）［29］、DFI［19］和土壤耕作指数（soil tillage index，STI）［16］共10种植被指数（表1），以比较各种不同植被指数对研究区植被盖度的估测能力。为避免云和水体的影响造成的异常值，对植被指数均做掩膜处理。其中，TSAVI、MSACRI的计算需要土壤线（土壤线是可见光红、近红外波段发射率或高度值之间的线性关系，是对土壤反射特征的综合描述。）斜率、截距因子，利用野外实测的31组土壤光谱曲线建立土壤线方程获得图2。

图2 土壤线Fig.2 Soil line

表1 用于植被盖度估算的几种主要植被指数Table 1 Several main vegetation indices for grassland coverage estimation

1.3.2 相关性分析相关性分析可以研究现象之间是否存在某种依存关系。采用Pearson相关系数计算10种植被指数与生长季/非生长季植被盖度的相关性，公式如下：

式中：rxy为相关系数；n表示样本数；xi表示样本中第i个因素的值；xˉ表示样本中某个因素的平均值；yi表示样本中第i个样本植被盖度的值；yˉ表示样本植被盖度的平均值。相关系数介于－1～1，绝对值越接近于1，则相关性越强，反之，相关性越弱。

1.3.3 植被盖度遥感反演模型像元二分模型是被广泛使用的植被盖度遥感反演模型之一，该模型要求所选植被指数与植被盖度呈线性相关，因此需要筛选出适宜表达生长季/非生长季光谱信息的植被指数构建光谱特征空间。而端元值的选择直接影响着模型的精度。仅利用植被指数构建的一维特征空间难以确定端元位置，加之影像噪声对端元值的选取有一定影响，使得像元二分模型的应用有着一定局限性。因此，本研究引入像元三分法提取端元特征值［3］，基于像元二分模型反演研究区生长季/非生长季植被盖度。

像元二分模型的原理是假设遥感数据中单个像元由土壤和植被两部分组成。生长季时，地表混合像元仅由光合作用植被和裸土组成；而非生长季，地表混合像元则由非光合作用植被和裸土组成。像元二分模型可如下式表达：

式中：S、SPV、SNPV和SBS分别代表混合像元、光合植被端元、非光合植被端元和裸土端元的遥感信息，fNPV为非光合植被盖度，fPV为光合植被盖度。

像元三分模型将混合像元进一步细分为PV、NPV和BS的线性组合，通过选择NDVI表征fPV、DFI表征fNPV来构建NDVI－DFI像元三分模型（图3）［30］。在理想情况下，NDVI－DFI二维特征空间表现为三角形。PV端元的NDVI值高、DFI值较低，位于三角形的右侧中部；NPV端元的NDVI值较低、DFI值最高，位于三角形的左上角；而BS端元的NDVI值和DFI值都很低，位于三角形的左下角；三角形的内部则是混合像元。在实际应用中，对于生长季而言，右侧中部顶点附近NDVI值高，DFI值较低；三角形左下顶点的DFI、NDVI值都较小，二者分别符合光合植被和裸土在特征空间的分布，因此，选择这两个顶点周围一定区域内的散点，取其平均值分别作为光合植被和裸土的端元值。对于非生长季而言，上部顶点附近DFI值高，NDVI值较低；三角形左下顶点的DFI、NDVI值都较小，两者分别符合非光合植被和裸土在特征空间的分布，因此，选择这两个顶点周围一定区域内的散点，取其平均值分别作为非光合植被和裸土的端元值。

图3 NDVI-DFI像元三分模型Fig.3 Tri-endmember linear mixture model with the NDVI and DFI

1.3.4 精度评价为了综合衡量不同植被指数基于像元二分法反演生长季/非生长季草地植被盖度的效果，本研究利用非生长季108个样点和生长季113个样点的实测盖度作为验证数据，将均方根误差（root mean square error，RMSE）、相对误差δ、决定系数（coefficient of determination，R2）作为精度评价指标。其中，相对误差能体现植被盖度反演值与实测值的偏离度，其值越大，模型精度越低。RMSE能够反映整个样本的总体精度，其值越小，估算精度越高。R2越大，表示拟合优度越大，即自变量对因变量的解释程度越高。公式如下：

2 结果与分析

2.1 不同植被指数与草地植被盖度相关性分析

为选出适宜生长季/非生长季植被盖度反演的植被指数，统计10个植被指数与生长季/非生长季草地植被盖度之间的相关系数（correlation coefficient，r）（表2）。可见，NDTI、MSACRI、S-NDTI、DFI和STI这5种植被指数与非生长季植被盖度的相关系数显著高于其他指数，适用于非生长季植被盖度的反演。NDVI、TSAVI、SAVI、RVI和MSAVI这5种植被指数与生长季植被盖度之间的相关系数高于其他指数，适用于生长季植被盖度的反演。适用于非生长季的指数中，NDTI、MSACRI、S-NDTI和STI的相关系数均为0.92（P＜0.01），DFI的相关系数与最大值相差较小。适用于生长季的指数中，NDVI、SAVI、MSAVI和TSAVI最大，均达到0.88（P＜0.01），RVI的相关系数接近于最大值。因此，仅靠相关系数无法筛选出适用于生长季/非生长季的最优植被指数，需要进一步比较。

表2 植被盖度与植被指数的相关性Table 2 Correlation between grassland coverage and vegetation index

2.2 植被盖度像元二分模型及精度评价

本节分别选取适用于生长季/非生长季的5种植被指数，引入像元三分法提取端元值，基于不同植被指数像元二分模型构建植被盖度反演模型，并根据精度评价结果，选出最优生长季/非生长季植被指数。从10种植被指数反演的高寒植被盖度结果（图4）来看，研究区植被盖度整体偏低，生长季植被盖度高的地方对应非生长季植被盖度也高，生长季植被盖度低的地方对应非生长季植被盖度也低。由像元二分模型植被盖度精度评价（表3）可知，10种植被指数的相对误差均小于45%，生长季最小相对误差可达9.23%，非生长季最小相对误差为10.13%，说明像元二分模型在高寒草地具有较好的适用性。以相对误差较小，RMSE较小，R2较大为优的原则选择最优反演模型。MSACRI的相对误差（10.18%）虽大于S-NDTI（10.13%），但其R2和RMSE均优于S-NDTI，具有最小的RMSE（7.78），其R2为0.81。因此，MSACRI最适用于反演非生长季高寒草地植被盖度，其余优劣顺序依次为S-NDTI、NDTI、STI、DFI。生长季时，NDVI具有最小的相对误差（9.23%）、最小的RMSE（8.43），其R2为0.77。因此，NDVI最适用于生长季高寒草地植被盖度反演建模，其余优劣顺序依次为TSAVI、MSAVI、SAVI、RVI。

表3 像元二分模型植被盖度精度评价Table 3 The accuracy evaluation of grassland coverage based on dimidiate pixel model

图4 当雄县生长季/非生长季高寒草地植被盖度分布Fig.4 Distribution of alpine grassland coverage in growing season/non-growing season in Dangxiong County

2.3 当雄县草地覆盖现状

基于以上结果，该部分使用MSACRI和NDVI指数反演的盖度数据分析研究区生长季/非生长季植被盖度的时空变化特征（图5～图6）。当雄县草地总面积为693172.10 hm2，主要分布于海拔4300～5100 m，生长季面积占比为72.45%。生长季平均植被盖度为50.54%，高寒特征明显，植被分布相对稀疏。低盖度（＜20%）草地主要分布于高海拔区域（＞5300 m），其面积占比为5.92%。中低盖度（20%～40%）草地一部分位于海拔4700～5100 m的念青唐古拉山北坡、纳木错湖周边，面积占比为11.15%；一部分位于海拔4300～4700 m的当雄县的西南部及县城周边，面积占比为6.72%。当植被盖度为40%～80%时，多分布于海拔4700～5100 m处，面积占比为27.66%。高盖度（＞80%）草地多分布于低海拔区域（＜4500 m），主要是河流与洼地形成的沼泽草甸。此外，随着海拔的增加，低盖度（＜20%）的草地逐步增加，植被盖度为20%～80%的草地于海拔4700～4900 m附近达到峰值，高盖度草地（＞80%）的比例逐渐减少。与生长季相比，非生长季草地盖度减小，其平均盖度为31.67%。约70.83%的草地盖度小于40%，草地仍集中分布于海拔4300～5100 m处，面积占比为82.17%。低海拔区域（＜4500 m）的沼泽草地在非生长季时盖度略微下降，盖度可达60%以上。

图5 不同海拔生长季植被盖度面积统计Fig.5 Statistics of vegetation coverage area in different altitude growing seasons

图6 不同海拔非生长季植被盖度面积统计Fig.6 Statistics of vegetation coverage area in different altitude non-growing seasons

3 讨论

3.1 不同植被指数与植被盖度间相关性差异

枯黄干草与绿色鲜草的光谱响应存在显著的差异，因此植被指数的拟合效果在生长季和非生长季也存在明显差异。草地在4月仍处于枯黄期，叶片水分和叶绿素下降，而纤维素和木质素增多。多光谱空间中，NPV与土壤可分性最高的是OLI6和OLI7波段，这两个波段位于吸水区附近，当植被变黄时，光谱响应随水分流失而增强［31］，因此由这两个波段构建的DFI、STI、S-NDTI、MSACRI和NDTI这5种植被指数与非生长季植被盖度相关性高。而NDVI、RVI、SAVI、MSAVI和TSAVI这5种植被指数由可见光－近红外波段构建，因此难以用于反演枯黄的盖度信息，这与Dai等［32］证明NDVI、SAVI无法区分NPV和土壤，NDTI可以区分NPV和土壤的结论一致。8月草地正处于生长最旺盛的状态，生物量、盖度等指标均达到较高水平，绿色植被在近红外波段具有独特的强反射的特性，且其反射率明显大于干枯植被［33］，因此NDVI、RVI、SAVI、MSAVI和TSAVI这5种植被指数与生长季植被盖度的相关系数远远高于其他指数。

3.2 生长季/非生长季像元二分模型构建

非生长季时，地表混合像元分解为NPV和BS。而DFI常用于区分NPV、PV和BS这3种组分共存的情况［30］，因此其相对误差最大，为44.02%。NDTI在STI的基础上进一步减弱大气噪声影响［34］，这与本研究中NDTI相对误差（20.09%）低于STI（29.96%）的结论一致。S-NDTI适用于无法估算土壤线的情况，其相对误差（10.13%）小于NDTI（20.09%），可见S-NDTI估算NPV的能力优于NDTI［29］。MSACRI中引入土壤线参数，能有效弱化土壤背景对草地光谱的影响，抗土壤噪声能力优于NDTI，适合于低盖度草地的监测［35］。MSACRI具有最小的RMSE（7.78），其R2高达0.81。因此，基于MSACRI的像元二分模型是反演非生长季植被盖度的最优模型。

研究区位于半干旱区，受高寒气候、土壤和水热条件的制约，生长季草地还是相对稀疏。土壤对植被指数的影响较为显著，而RVI对高植被覆盖区域较为敏感［36］，因此RVI并不适用于该地区。SAVI在削弱土壤对植被信息影响的同时也降低了与植被盖度的相关性［37］，其反演精度最低（77.36%）。尽管MSAVI适用于低植被盖度的反演，但野外采集的高光谱曲线构建的TSAVI考虑了土壤斜率和截距［38］，因此其相对误差（11.20%）略小于MSAVI（13.70%）。研究区为典型的高寒草甸与寒冷半干旱高寒草原类型之间的过渡地，生长季植被盖度集中于20%～80%，NDVI对植被盖度的检测幅度较宽，这可能也是TSAVI反演植被盖度的能力不如NDVI的原因。因此基于NDVI的像元二分模型是反演生长季植被盖度的最优模型。

3.3 当雄县植被覆盖现状分析

当雄草地多位于海拔4300～5100 m，高寒特征明显，海拔对植被盖度具有决定性的影响。随海拔的增加，气温显著降低，而降水则呈先增加后降低的趋势，盖度呈现先增加后减少的单峰变化格局［39］。西南部低海拔地区（4300～4700 m）植被受到干旱胁迫，形成以针茅属（Stipa）物种为主的草原化草甸，土壤水分条件较差，沙化严重，因此植被盖度普遍偏低，生长季约占20%～40%；当雄县城附近，受到频繁的人类活动及过度放牧的影响，植被盖度相对较低；而该海拔梯度高盖度（＞80%）的区域多位于河流与洼地形成的以藏北蒿草（Kobresia littledalei）、矮生嵩草为主的沼泽草甸，土壤水分充足，因此植被长势较好。高海拔区域（4500～5100 m）降水充足，形成以高山嵩草为主的高寒草甸，植被长势较好，盖度多为60%～80%，面积占比为35.33%。随着海拔的进一步增加，当海拔超过5300 m时，气温的下降限制了植被对水分的吸收，高覆盖植被（＞80%）几乎消失。

8月伴随着雨季开始，气温升高，草地达到生长高峰期。4月为冬春交接的时间段，温度低，降水少，非生长季叶片枯萎蜷缩导致盖度减小，各海拔区间非生长季植被盖度相较于生长季植被盖度均有所下降，非生长季平均植被盖度较生长季下降18.87%。对于盖度较高的沼泽草甸而言，其盖度随季节的变化小，非生长季植被盖度相较于生长季而言，盖度减少幅度不大。

此外，本研究只考虑了非生长季NPV、BS和生长季PV、BS仅有两种组分存在的情况，缺少对NPV、PV和BS这3种组分共同存在情况的考虑。植被生态系统自身的复杂性会影响植被盖度估算的精度，因此后续将进一步考虑混合时期的植被盖度情况。另外，由于研究区位于高海拔地区，容易受到大气、云和雪的影响，实测数据与遥感数据在时空上常不能完全吻合。且需要考虑不同草地类型、含水量、分解程度及土壤类型、土壤湿度对拟合效果的影响。