青藏高原草地非光合植被覆盖度时空动态

牛慧兰，任鸿瑞

（太原理工大学测绘科学与技术系, 山西 太原 030024）

植被是陆地生态系统最基本的组成部分，根据植被的功能可以将植被分为光合植被和非光合植被[1]。非光合植被由地表植被的凋落物和植被衰败后的枯枝、枯茎组成[2]，存储了生态系统中大量有机碳和养分[3]，在平衡地表能量、调节土壤养分、调节净初级生产力和固碳等方面发挥了重要作用[4]，在控制风蚀、水蚀以及监测火灾发生的频率和强度等方面非光合植被发挥的作用也至关重要[1]。因此，估算非光合植被覆盖度对于估测生态系统的碳储量、监测和合理利用自然资源具有重要意义。

目前，利用遥感技术估算大范围非光合植被覆盖度在国内外的研究较少，现有的方法主要是回归模型法和混合像元分解法[5-6]。回归模型法即建立采样点的地表实测覆盖度与植被指数的经验关系，从而估算研究区内非光合植被覆盖度。如Daughtry等[7]利用纤维素吸收指数(cellulose absorption index,CAI)建立回归模型估算了爱荷华州中部的农作物残茬覆盖度，估算结果和实测数据决定系数R2达到0.85；柴国奇等[8]利用干枯燃料指数(dead fuel index, DFI)建立回归模型估算了锡林郭勒草原fNPV，估算结果和实测数据的决定系数R2达到0.63。可以看出，建立回归模型法估算精度较高，但该方法建立经验模型需要大量的实测数据，而且只适用于特定的区域或植被类型。

混合像元分解法是根据某一类像元反映的遥感信息将其从混合像元中分解出来，从而获取植被覆盖度，分为非线性光谱分解法和线性光谱分解法[9-10]。非线性光谱分解法相对于线性光谱分解法模型更加复杂，并且没有明显的精度提升[6,11]，因此大多选用线性光谱分解法估算植被覆盖度。线性光谱分解法目前主要有像元三分模型和像元二分模型，像元三分模型已被国内外学者用于估算fNPV，该方法比像元二分模型更加复杂，且纯净端元的的选取容易受到植被类型的影响[12]。像元二分模型已经广泛应用于光合植被覆盖度的估算，利用归一化植被指数(normalized difference vegetation index, NDVI)构建的像元二分模型估算光合植被覆盖度并且估算精度达到85%[13]，但利用像元二分模型估算fNPV的研究还很少。因此，像元二分模型对fNPV的估算潜力还有待进一步挖掘。

Van Deventer 等[14]在1997 年首次提出利用TM影像的B5 和B7 波段构建NDTI 指数，并用于区分土壤类型。Serbin 等[15]研究了CAI、LCA、NDTI、NDSVI、NDI5、NDI7 指 数 与ƒNPV的 相 关 性，发现NDTI 指数与ƒNPV的相关性仅次于CAI 指数与LCA指数。并且此研究结果还表明NDTI 指数能够分离土壤与非光合植被，并且能够估算fNPV。常伟[10]分析了NDTI 指数与fNPV的相关性，得到NDTI 指数与fNPV的相关系数达到了0.60。冷若琳等[16]用STI、NDTI、DFI 等9 种植被指数建立了草地非生长季植被覆盖度回归模型，对比得到最优植被指数为NDTI，且R2达到0.407。国内外学者的研究表明了利用NDTI 指数估算fNPV的可行性，因此本文利用NDTI 指数构建像元二分模型，并对青藏高原草地fNPV进行估算，从而进一步研究像元二分模型估算fNPV的准确性。

草地是青藏高原分布最为广泛的植被类型，以高寒草甸、高寒草原为主要草地类型。青藏高原草地有丰富的碳储量，对于维持青藏高原生态系统的多功能性和碳循环具有重要意义。本研究以青藏高原草地为研究区，以MODIS 500 m 空间分辨率地表反射率产品 MOD09A1 为数据源，利用NDTI 指数构建像元二分模型估算研究区fNPV，并分析fNPV的时空变化特征，以及与气候因素的关系，为青藏高原的非光合植被动态监测提供依据。

1 研究区概况与数据

1.1 研究区概况

青藏高原面积约占中国陆地面积的26.8%[17]，地处26°00′12″～39°46′50″N，73°18′52″～104°46′59″E。青藏高原呈现西北高东南低的地势特征，平均海拔高度在4 000 m 以上，并且地貌类型复杂多样，有“世界屋脊”之称。青藏高原气候区域差异显著，研究区内气温呈现出从南向西北逐渐递减的变化特征，降水量整体表现为东南部降水多、西北部降水少的分布情况。

1.2 数据源与预处理

1.2.1 MODIS 数据产品

本研究中使用的地表反射率产品MOD09A1 和植被指数产品MOD13A1 数据来自Terra 卫星的MODIS 传感器，由美国国家航空航天局(https://nasa.gov)提供，空间分辨率均为500 m，时间分辨率分别为8 d 和16 d。青藏高原植被枯黄期主要集中在每年的第260-300 天[10]，由于随着时间推移积雪覆盖会对估算fNPV产生影响，故选择第273-280 天内的MOD09A1 影像数据对研究区fNPV进行估算。选择2020 年第209-224 天内的MOD13A1 数据用于估算研究区绿色植被覆盖度。由于MOD09A1 产品和MOD13A1 产品均已经过辐射校正、云体掩膜、大气校正处理，仅需对产品进行拼接、裁剪、重投影处理即可。

1.2.2 土地利用数据

土地利用数据引自全国地理信息资源目录服务系统提供的2010 期全球30 m 地表覆盖率数据集(https://www.webmap.cn/mapDataAction.do?method=globalLandCover)，该数据集提供了10 个主要的地表覆盖类型，分别为耕地、林地、草地、灌木地、湿地、水体、苔原、人造地表、裸地、冰川和永久积雪。本研究将30 m 空间分辨率青藏高原的土地利用数据重采样至500 m，提取出青藏高原草地部分进行研究。

1.2.3 气象数据

气温和降水数据是由国家地球系统科学数据共享服务平台(http://www.geodata.cn/) 提供的1 km分辨率逐月气温数据和逐月降水数据。该数据集根据CRU 发布的全球0.5°气候数据集和WorldClim发布的全球高分辨率气候数据集在中国地区降尺度生成[18]。将对应月份的数据计算平均值得到平均气温，将对应月份的月降水量数据相加得到累积降水量。

2 研究方法

2.1 归一化差异耕作指数

Van Deventer 等[14]将TM 影 像 的1-7 波 段 以加、减、乘、除以及归一化的形式组合，并进行回归分析，结果表明B5 和B7 波段对于识别农耕地更有效，故利用TM 影像B5 和B7 波段构建了归一化差异耕作指数(NDTI)。Jacques 等[19]研究表明基于TM数据B5、B7 波段构建的NDTI 指数同样适用于MODIS 数据对应的B6、B7 波段。计算公式如式(1)所示：

2.2 非光合植被像元二分模型

像元二分模型假设像元只由植被和土壤两部分构成，混合像元S由植被覆盖部分Sv和土壤覆盖部分Ss组成：

设只有植被覆盖的纯像元为Sveg，植被在混合像元中所占的面积比例为fveg，土壤在混合像元中所占的面积比例为(1 -fveg)，混合像元中植被覆盖部分Sv可以表示为：

同理，设只有土壤覆盖的纯像元为Ssoil，混合像元中土壤覆盖部分Ss可以表示为：

混合像元可以看作植被和土壤两个组分的线性组合，其各自在混合像元中所占的面积比例分别为两因子的系数，其中植被所占的面积比例为该像元的植被覆盖度。即：

本研究构建的非光合植被的像元二分模型，以NDTI 指数作为遥感信息的观测值，将混合像元看作只由非光合植被和裸土两部分组成。混合像元的NDTI 值为非光合植被和土壤两个组分的线性组合，即：

式中：fNPV为混合像元中非光合植被所占的面积比例；NDTINPV为纯非光合植被像元的NDTI值；NDTIsoil为纯土壤像元的NDTI值。本研究设置NDTI累积频率为5% 的像元值作为NDTIsoil，NDTI累积频率为95%的像元值作为NDTINPV。

将式(6) 进行变形，非光合植被覆盖度fNPV可以表示为式(7)的形式：

2.3 一元线性回归模型

一元线性回归模型能够分析变量长时间序列变化趋势和演变规律，本研究采用一元线性回归模型逐像元对青藏高原草地非光合植被覆盖度变化趋势(slope)进行分析，计算方法如式(8)所示：

式中：n为监测年数，i= 1, 2, …, 21 (从2000 年开始i= 1 递增至2020 年)，fNPV为青藏高原草地非光合植被覆盖度。当slope ＞ 0 时表现为增长趋势，当slope ＜0 时表现为减少趋势。

显著性检验采用F检验：

2.4 相关性分析

两个要素之间的相关系数rxy体现了要素之间相互关系的密切程度。rxy取值范围在[-1, 1]区间内，当rxy＜ 0 时表示x和y负相关，当rxy＞ 0 时表示x和y正相关，rxy绝对值越接近于1 表示两要素的相关关系越密切，rxy越接近于0 表示两要素的相关关系越不密切。计算方法如式(11)所示：

式中：rxy、rxz和ryz分别表示要素x和y、x和z、y和z之间的简单相关系数。

本研究以青藏高原草地fNPV、降水量和气温作为要素，两两组合，分别进行相关性分析和偏相关性分析。为研究青藏高原草地fNPV受气温和降水影响的差异程度，分别对1 月-12 月，5 月-9 月和10 月-翌年9 月3 个时间段内fNPV与气温和降水的相关性、偏相关性进行对比分析。

3 结果和分析

3.1 青藏高原草地非光合植被空间分布特征

通过逐像元计算21 年内研究区fNPV的平均值，得到2000-2020 年青藏高原草地fNPV空间分布特征(图1)，青藏高原草地西北、西南地区fNPV较低，东北、东南地区fNPV较高，且呈现出由西北至东南逐渐增加的分布情况。

图1 2000-2020 年研究区平均fNPV 分布情况Figure 1 Average value of fNPV in the study area from the year 2000 to 2020

3.2 青藏高原草地非光合植被年际变化特征

根据2000-2020 年青藏高原草地fNPV年际变化趋势(图2a)，从2000 年至2020 年青藏高原草地年 平 均fNPV在0.341 (2015 年) 至0.444 (2018 年)之间波动，上升趋势不显著(P＞ 0.05)，上升速率为每年 0.001 2。

以2000-2020 年青藏高原草地fNPV空间变化趋势图揭示了青藏高原草地fNPV空间变化趋势(图2b)，从2000-2020 年青藏高原草地fNPV趋于上升的面积占比(13.39%)较小，主要分布于西北地区部分草地，fNPV趋于下降的面积占比(86.61%) 较大，主要分布在东南地区。青藏高原草地fNPV变化显著(0.01 ≤P＜ 0.05)的区域占比较少，其中，显著增加的草地区域占0.001 5%，显著减少的草地区域占0.000 4%。从2000-2020 年青藏高原草地fNPV整体变化趋势处于稳定状态。

图2 2000-2020 年研究区fNPV 年际变化趋势(a)和空间变化趋势(b)Figure 2 The interannual (a) and spatial (b) trend of fNPV in the study area from the year 2000 to 2020

3.3 青藏高原草地非光合植被覆盖度与气象因素的关系

从图3 中可以看出fNPV变化与累积降水量变化呈正相关关系的区域主要分布在青藏高原的中部、东北部和东南部，约占青藏高原草地面积的75%以上。fNPV变化与温度变化呈正相关关系的区域主要分布在青藏高原的中部和东北部，约占青藏高原草地面积的50%，fNPV与温度呈负相关关系的地区主要分布在青藏高原草地的南部、西南部和东南部。青藏高原草地fNPV与5 月-9 月均温呈正相关关系所占的面积为54.6%，高于1 月-12 月和10 月-翌年9 月所占面积。fNPV与全年降水量呈正相关关系的面积为75.6%，高于5 月-9 月和10 月-翌年9 月所占面积。

图3 研究区2000-2020 年fNPV 分别于1 月-12 月、5 月-9 月、10 月-翌年9 月累积降水量、气温的相关关系Figure 3 Correlation between the fNPV of the study area and the cumulative precipitation and between the fNPV of the study area and the temperature for the months of January-December, May-September, and October-September of the year 2000 to 2020

从图4 中得到剔除降水影响后fNPV与温度呈正相关关系所占的面积平均增加了2.7%，剔除温度的影响后fNPV与降水呈正相关关系的所占面积无明显变化。fNPV与年均温的相关性受降水因素的影响呈负相关关系，相关系数为 -0.008 (表1)，剔除降水因素的影响后fNPV与年均温的偏相关性为正相关性。fNPV与降水的相关性几乎不受到温度的影响。

图4 研究区2000-2020 年fNPV 分别于1 月-12 月、5 月-9 月、10 月-翌年9 月累积降水量、气温的偏相关关系Figure 4 Partial correlation between the fNPV of the study area and the cumulative precipitation and between the fNPV of the study area and the temperature for the months of January-December, May-September,and October-September of the year 2000 to 2020

表1 研究区2000-2020 年fNPV 与累积降水量、气温的相关关系Table 1 Correlation coefficient between fNPV and cumulative precipitation and between fNPV and temperature in the study area from the year 2000 to 2020

4 讨论

本研究利用NDVI 指数构建像元二分模型估算了研究区2020 年生长季绿色植被覆盖度(fPV) (图5)，并于2020 年非光合植被覆盖度进行对比。结果表明，研究区非光合植被覆盖度分布情况(图5a)与生长季绿色植被分布情况(图5b)相同，呈现从西北至东南方向逐渐升高的特征，从而间接证明了非光合植被像元二分模型估算青藏高原草地非光合植被覆盖度的准确性。NDTI 指数也会受到土壤类型、植被类型、土壤水分等多种因素影响[20]，NDTI指数更适用于土壤(BS) 和非光合植被(NPV) 两种组分存在的情况[4]。青藏高原地形复杂，植被枯黄期不同，因此当研究区内存在还未完全枯黄的绿色植被时，可能会对估算结果产生影响。

图5 2020 年青藏高原草地光合植被覆盖度(a)和光合植被覆盖度(b)Figure 5 Grassland fNPV (a) and fPV (b) in the Qinghai-Tibet Plateau in the year 2020

利用MOD09A1 产品估算fNPV，发挥了MODIS产品低成本、长时间、大范围观测的优势，但MODIS产品空间分辨率(500 m × 500 m)较低，难以开展大范围的地面调查与估算结果对应。未能利用实测数据对估算的非光合植被覆盖度进行验证是本研究的不足之处。

青藏高原草地fNPV呈现西北低东南高的分布特征，这与青藏高原复杂的地形、气候条件等关系密切。fNPV低的地区主要分布在高海拔的藏北高原地区、祁连山地西部和青海高原，这些地区由于受到高大山脉的阻挡使得暖湿空气很难到达，降水稀少，导致植被生长稀疏。fNPV高的地区主要分布在海拔较低的藏南谷地、川藏高山峡谷区和祁连山地东部，这些地区属于湿润、半湿润气候区，大多为森林草原地带，降水量多，植被生长茂盛。2000-2020年青藏高原草地平均fNPV呈现不显著上升趋势，上升速率为每年0.001 2。在空间分布上，研究区西北部fNPV没有明显变化，在川藏高山峡谷区的部分地区fNPV有不显著下降趋势。可能是相对于fNPV显著下降地区，fNPV显著上升地区面积占比更大。青藏高原草地的中部、东北部和东南部的草地类型主要为高寒草原和高寒草甸，这些地区降水年际变化率高，降水量多的地方土壤含水量增加促进植被生长，由于非光合植被由绿色植被转化而来，形成的非光合植被覆盖度增加，因此非光合植被覆盖度与降水量相关系数较高。在青藏高原的西南地区降水量年际变化率低，降水量对植被覆盖度的影响差异较小。在青藏高原的东南和西南地区在研究时间段内正处于由光合植被转变为非光合植被过程，温度越低光合植被转变为非光合植被的速度越快，因此fNPV与温度呈负相关关系；而在青藏高原草地的中部和东北部植被已经枯黄，由于温度的升高能促进绿色植被的光合作用使植被生长更加茂盛，因此绿色植被转变为非光合植被的数量更多。

5 结论

本文利用MOD09A1 产品计算了青藏高原草地的NDTI 指数，并根据NDTI 指数构建非光合植被像元二分模型，对2000-2020 年青藏高原草地fNPV进行估算，并将估算结果与气候因素进行分析，得到以下结论： 1) 2020 年青藏高原草地fNPV分布情况与生长季绿色植被分布情况相同，证明非光合植被像元二分模型可以有效估算青藏高原草地fNPV。2)青藏高原草地fNPV呈现西北低东南高的分布情况，主要受到青藏高原复杂的地形与气候的影响。从2000-2020 年青藏高原草地平均fNPV呈现不明显上升趋势，在高原东南地区，部分草地的fNPV存在减少情况，但变化不明显。3)青藏高原草地fNPV与降水主要呈正相关关系，在西南部分地区呈负相关关系，剔除温度影响后fNPV与降水的相关性几乎没有变化。在研究区中部和东部地区fNPV与温度呈正相关关系，东南、南部和西南地区与温度呈负相关关系，研究区内约2.7% 的草地在剔除降水因素的影响后，fNPV与温度由负相关关系转变为正相关关系。4)青藏高原草地fNPV与5 月-9 月均温呈正相关关系所占的面积，高于1 月-12 月和10-翌年9 月所占面积。fNPV与年降水量呈正相关关系所占面积，高于5 月-9 月和10 月-翌年9 月所占面积。