基于Landsat 8的若尔盖县沙化监测区NDVI植被覆盖变化特征分析

刘朔, 杨建勇, 蔡凡隆

四川省林业和草原调查规划院，四川 成都 610081

随着遥感技术的发展，植被覆盖研究成了分析生态环境、气候演变和水文过程等众多领域的基础[1]。归一化差值植被指数（NDVI）是反映植被覆盖的一个重要参数，具有较强的植被检测能力，正值表示有植被覆盖，且数值随覆盖度增大而增大，NDVI不仅能更好地检测植被的生长状况，并且能消除部分辐射误差[2-5]。像元二分模型是一种实用的遥感估算模型，是采用NDVI归一化植被指数定量估算植被覆盖度的模型，在植被覆盖度研究和生态研究及评价等方面应用广泛，许多学者基于 NDVI 数据和像元二分模型做了多方面研究[6-10]，一般通过NDVI、植被覆盖度与植被相关的指标分析植被覆盖变化和趋势。国内学者基于NDVI对毛乌素沙地、新疆、西藏高原等沙化广布区进行了分析研究[11-14]，基础数据主要为MODIS_NDVI 数据，该数据空间分辨率为250 m，适合大尺度、广区域的研究区。基于NDVI对川西北沙区的研究鲜有报道，重要原因之一是川西北沙区沙化土地分布分散，主要呈小块状，点状分布，MODIS_NDVI 数据精度难以满足对该区域沙化土地的分析需求。

若尔盖县是川西北沙化最严峻的地区，掌握若尔盖沙地动态变化趋势有利于该区域生态环境保护与恢复，为该区域防沙治沙提供理论依据与数据支持。本文以若尔盖县第五次沙化监测图斑区为研究区，以2013年、2015年、2017年3个年度10月初的 Landsat 8影像作为年际对比基础数据，主要原因是10月初植被已处于生长末期，更有利于提取低植被区沙化土地，减少将沙化土地植被过高判读的可能，加之研究区以草地植被为主，植被单一，季节间生长规律同步性好。通过Landsat 8影像中的近红外波段和红光波段数据提取空间分辨率为30 m的NDVI，再利用像元二分模型反演植被覆盖度，分析研究区植被变化特征，评估研究时段期间植被变化趋势，为若尔盖县沙化治理和生态环境建设对策制定提供科学依据，应用的技术方法可为川西北沙化监测、沙区植被动态变化对比和评估提供科技支撑。

1 研究区概况

研究区为若尔盖县全国第五次沙化监测的图斑区，总面积79 851.78 hm2，主要分布在牧区，以草地沙化为主,平均海拔3 500 m，地质母岩多为三叠系砂岩、板岩、灰岩和第四纪的松散堆积物，其形成的土壤含沙量高。研究区属大陆性季风高原气候区，气候寒冷干燥，长冬无夏，日照强烈，温差大，年均气温0.7℃ ，无绝对无霜期，年均日照2 389 h；年均降雨量656.8 mm，年均蒸发量1232 mm，年均湿度69%；平均风速2.4 m·s-1，最大风速40 m·s-1，年大风日数50 d，年沙尘日数6 d。

2 研究方法

2.1 数据来源及处理

（1）研究采用Landsat 8遥感影像进行NDVI提取。Landsat 8遥感影像来自“地理空间数据云”网站，分别是2013年10月、2015年10月、2017年10月的Landsat 8遥感影像各1景（见表1），空间分辨率为30 m（仅全色波段为15 m）；（2）研究区范围数据来源若尔盖县第五次全国沙化监测数据。

表 1 Landsat 8影像主要信息表Tab. 1 Main information table of Landsat 8 image

2.2 数据处理及研究方法

数据总体处理流程为：（1）在ENVI5.3中主要完成NDVI和FC（植被覆盖度）提取，主要步骤分别为：①辐射定标；②大气校正；③NDVI计算；④FC（采用像元二分模型反演植被覆盖度）。（2）在Arcgis10.2中主要进行NDVI和FC的分级、植被变化的趋势分析、转移矩阵计算。主要步骤分别为：①通过无偏移天地图卫星影像将研究区矢量图层、NDVI栅格、FC栅格校正为大地2000投影坐标系；②按掩膜提取研究区的NDVI和FC栅格；③进行NDVI和FC的分级；④植被变化转移矩阵计算；⑤在研究区创建生成1 000个随机点矢量数据，用来提取各年度NDVI和FC的栅格数据值，并导出。（3）在Spss24中对导出的研究区1 000个矢量点导出数据进行数据类型分布探索的基础上，进行各年度间NDVI和FC差异分析，以及与沙化程度的相关分析。经Spss24数据探索，整体数据呈偏态分布，非正态分布（这是研究区域数据本身决定的，因为植被占绝大部分，而非植被的裸地等占比很小），因此选用Spss24 中“非参数检验”中“k个相关样本”进行差异分析。

主要研究方法如下：

（1）NDVI提取

NDVI 是遥感监测中反映植被在可见光、近红外波段与土壤背景之间光学差异的指标，是对地表植被活动的简单、有效和经验的度量[15]。计算公式为：

NDVI=（Rnir - Rr ） /（Rnir + Rr）

式中：Rnir为近红外波段反射率； Rr为红波段反射率。

在ENVI5.3中必须先对影像进行辐射定标和大气校正，再通过 “band math” 计算工具中输入公式或者使用软件中专门的NDVI工具进行求算。在Arcgis10.2中根据NDVI数据取值及分布特征并结合经验成果，将其重分类为4个等级，分别为：低（-1～0.2）；较低（0.2～0.4）；中等（0.4～0.5）；较高（0.5～1）。

（2）反演植被盖度

植被覆盖度(vegetation fractional cover)简称FC，在众多遥感测量植被覆盖度的方法中，较为实用和普遍的方法是采用像元二分模型利用植被指数近似估算植被覆盖度，通过一系列的模型运算[16]，可以得到：

式中：NDVIsoil为完全是裸土或无植被覆盖区域的 NDVI 值，NDVIveg代表完全被植被所覆盖像元的NDVI 值，即纯植被像元的 NDVI 值。

受众多因素影响， NDVIsoil并不是理论上的零值，通常在 -1～0.2 之间变化，NDVIveg也随植被类型以及时间和空间而所差异，因此NDVIsoil和NDVIveg没有固定值[17]。实际应用中，在没有实测数据的情况下，可根据经验取一定置信度范围内的NDVI max 和 NDVI min 作为和DVIveg和NDVIsoil则有：

FC在ENVI5.3中可通过 “band math” 计算工具中输入公式进行求算，本研究通过多次测算和对比，将NDVI为-1～0之间的累计频率作为NDVI min，将累计频率为 90% 的 NDVI作为NDVI max。并将FC数据取值及分布特征并结合前人经验，将其重分类为4个等级：低（0～0.2）；较低（0.2～0.4）；中等（0.4～0.5）；较高（0.5～1）

（3）转移矩阵分析

Arcgis10.2中将各期NDVI和FC栅格数据转化为矢量数据，再分别进行叠加相交分析后导出数据，并在EXCEL中通过数据透视表整理提取面积转移矩阵。

3 结果与分析

3.1 NDVI及FC年总体变化特征

由图1和图2和表2可知，研究区NDVI和FC从2013年、2015年、2017年均呈显著增加（差异均达到极显著），NDVI中位数值从2013年的0.477 5，增加到2015年0.537 4，至2017年达0.592 1；FC中位数值从2013年的0.530 5，增加到2015年0.597 1，至2017年达0.657 8，说明研究区植被覆盖度明显增加。主要原因是2007年起若尔盖县逐步开始了持续沙化治理工程，包括：省级沙化治理试点工程以及后续的成果巩固工程、 2013年开始实施的国家川西藏区沙化治理工程、国家重点生态功能区转移支付资金项目等，工程中相关围栏封育、植灌种草、人工种草、生物沙障等措施成效显著[18-20]。同时，近年来若尔盖县持续开展畜牧业提档升级，推行畜种改良，以草定畜等，草地超载情况得到明显改善，2012年末全县牲畜存栏114.1万混合头[21]，到2017年降低为95.96万混合头[22]。

图 1 NDVI各年度中位数变化图Fig. 1 Annual median change chart of NDVI

图 2 FC各年度中位数变化图Fig. 2 Annual median change chart of FC

图 3 NDVI 变化趋势图Fig. 3 NDVI change trend

图 4 FC变化趋势图Fig. 4 FC change trend

3.2 NDVI及FC动态变化转移矩阵

由表3、表4、表5可知，研究区分级后的NDVI从2013年、2015年、2017年均呈显著增加。其中：2013年至2015年，较高等级面积占比从7.26%提高到30.69%；中等级面积占比从75.97%降低到59.31%；较低等级面积占比从13.56%降低到7.32%；低等级面积占比从3.2%降低到2.67%。其具体转化方向主要是：中等级向较高等级转化，转化面积占自身面积31.9%；中等级向较高等级转化，转化面积占自身面积56.9%；低等级向较低等级转化，转化面积占自身面积30.2%。

表 2 研究区各年度NDVI和FC差异显著性系数（Bonferroni调整）Tab. 2 Significant coefficient of NDVI and FC in each year in the study area (Bonferroni adjustment)

表 3 研究区NDVI和FC分级差异显著性系数（Bonferroni调整）Tab. 3 Significant coefficient of NDVI and FC classification grades in the study area (Bonferroni adjustment)

表 4 研究区2013—2015年NDVI动态变化矩阵Tab. 4 Dynamic change matrix of NDVI in the study area from 2013 to 2015

表 5 研究区2015—2017年NDVI动态变化矩阵Tab. 5 Dynamic change matrix of NDVI in the study area from 2015 to 2017

2015年至2017年，较高等级面积占比从30.69%提高到46.22%；中等级面积占比从59.31%降低到45.87%；较低等级面积占比从7.32%降低到4.94%；低等级面积占比从2.67%增加到2.96%。其具体转化方向主要是：中等级向较高等级转化，转化面积占自身面积34.05%；中等级向较高等级转化，转化面积占自身面积49.39%；低等级向较低等级转化，转化面积占自身面积27.88%。

由表3、表6、表7可知，研究区分级后的FC从2013年、2015年、2017年均呈显著增加。其中：2013年至2015年，较高等级面积占比从64.92%提高到82.14%；中等级面积占比从25.45%降低到10.71%；较低等级面积占比从6.88%降低到4.85%；低等级面积占比从2.74%降低到2.29%。其具体转化方向主要是：中等级向较高等级转化，转化面积占自身面积71.41%；中等级向较高等级转化，转化面积占自身面积30.08%；低等级向较低等级转化，转化面积占自身面积29.42%。

表 6 研究区2013—2015年FC动态变化矩阵Tab. 6 Dynamic change matrix of FC in the study area from 2013 to 2015

表 7 研究区2015—2017年FC动态变化矩阵Tab. 7 Dynamic change matrix of FC in the study area from 2015 to 2017

2015年至2017年，较高等级面积占比从82.14%提高到89.41%；中等级面积占比从10.71%降低到4.07%；较低等级面积占比从4.85%降低到3.92%；低等级面积占比从2.29%增加到2. 60%。其具体转化方向主要是：中等级向较高等级转化，转化面积占自身面积77.41%；中等级向较高等级和高等级转化，转化面积分别占自身面积23.24%和25.03%；低等级向较低等级转化，转化面积占自身面积28.31%。

3.3 NDVI及FC与沙化程度相关性分析

将研究区NDVI 及FC等级与研究区2015年监测的沙化程度（轻度、中度、重度、极重度）进行相关分析，发现它们均呈显著正相关（见表8）。主要原因是沙化程度的划分主要依据之一就是植被覆盖度，而NDVI和FC均是反映植被生长情况的指标。另一方面也说明，本次NDVI 及FC重分类等级的划分基本符合研究区沙化植被覆盖度现状。由于全国沙化监测工作是每5年开展一次，时间跨度长，利用NDVI及其反演的植被覆盖度的方法理论上可在不同时段（如年际、月际），对研究区沙化植被动态变化进行分析评估，是一种非常灵活高效的方法。

表 8 研究区NDVI 及FC等级与沙化程度相关性Tab. 8 Correlation between NDVI, FC grades and the desertification degree in the study area

4 讨论

（1）Landsat 8作为免费的影像数据，自2013年来已被广泛应用于自然资源调查、生态环境监测、地质灾害响应、教育科研和政府管理等领域，其空间分辨率为30 m，较MODIS影像（空间分辨率为250 m）高很多，可进行更加细致和精确的研究。但是，其拍摄周期间隔16 d，相对较长，因此影像可选择的数量有限，加之数据若受到气候、天气等因素的影响，云量较大的现象也较普遍。笔者搜索了研究区2013—2020年4—10月的Landsat 8影像，再排除受云量影响较大影像后，结合研究内容特点，最终筛选确定了2013年、2015年、2017年共3个年度10月的影像数据进行研究。不同于MODIS可直接下载合成好的NDVI数据，Landsat 8需要根据红光波段和近红外波段求算NDVI值，在此过程中，辐射定标和大气校正是必须进行预处理工作，这是基于Landsat 8准确求算NDVI的前提，未进行大气校正测算的NDVI值要普遍小于大气校正后测算的NDVI值。

（2）研究区2013年、2015年、2017年的NDVI和FC均呈显著增加，由于NDVI和FC指标变化趋势一致，都可反映植被盖度变化趋势，实际生产过程中可根据研究需要灵活选择指标进行研究分析。

（3）在实际生产工作中往往需要对植被盖度进行分级，确定各等级间的区划阈值最为关键，但区划等级间的阈值没有统一标准，要综合影像、研究目的、研究方法和研究区本底情况或结合已有研究而确定。本研究中，研究区NDVI 及FC等级与研究区2015年监测的沙化程度（轻度、中度、重度、极重度）均呈显著正相关，说明本次NDVI 及FC重分类等级的划分基本符合研究区沙化植被覆盖度现状。在沙化监测及成效对比评估工作中，充分利用遥感数据的同步性、时效性、周期性和成本优势，深入研究NDVI与工作区植被盖度关联特征，将有效提升监测、评估工作的时效，降低因调查时间差异或调查者个体差异产生的判读误差，为若尔盖县沙化治理和生态环境建设对策制定提供科学依据。