基于NDVI的2009—2018年黄河流域林草植被覆盖变化*

郝家田 胡云云 杜一尘 侯晓巍 向安民

(国家林业和草原局西北调查规划设计院 旱区生态水文与灾害防治国家林业和草原局重点实验室 西安 710048)

植被覆盖动态变化是土地资源和环境变化研究的重要内容(王情等， 2013)。地表植被覆盖退化能引发水土流失、沙漠化、湖泊萎缩、碳汇功能丧失等一系列严重问题(郑子彦等， 2020)。目前，遥感技术已成为区域植被覆盖变化监测的主要手段(Martinezetal.， 2009)，能够宏观、快速地获取植被变化。归一化植被指数(NDVI)能很好地反映植被生长状况，是研究植被覆盖变化的最常用遥感监测指标(Baretetal.， 1991； Pouliotetal.， 2009)。Chen等(2019)基于MODIS、NOAA等数据，研究全球2000～2017年的植被变化，发现全球植被覆盖呈增加趋势，其中以中国和印度贡献最大。许翔驰(2019)、李杰等(2019)用MODIS NDVI时间序列数据集，分别研究了我国和中老缅交界区的植被覆盖时空变化特征。

黄河流域不仅是人类文明发展的摇篮，是国家重要的农牧业生产基地和能源基地，也是我国重要的生态安全屏障。黄河流域上游的水能资源、中游的煤炭资源、下游的石油和天然气资源都十分丰富，被誉为中国的“能源流域”(贺振等， 2012)。黄河流域是连接青藏高原、黄土高原和华北平原的生态廊道，拥有三江源、祁连山等国家公园及多个自然保护区等国家重点生态功能区。2019年9月18日，习近平总书记在黄河流域生态保护和高质量发展座谈会上强调，黄河流域是我国重要的生态屏障和重要的经济地带，是打赢脱贫攻坚战的重要区域，在我国经济社会发展和生态安全方面具有十分重要的地位(习近平， 2019)。黄河流域的植被覆盖变化以及引起的生态环境变化，对流域经济社会可持续发展有着重要影响。因此，开展黄河流域植被覆盖时空变化特征研究非常重要。目前，黄河流域植被覆盖变化研究取得了一定成果，刘绿柳等(2006)用1982—1999年的NOAA/AVHRR分析了黄河流域NDVI与气温、降水的年际变化趋势及相互关系的时空变化规律，发现流域NDVI与降水、温度相关显著的植被类型以草地、灌木为主，但相关区域的空间位置随时间发生变化； 杨胜天等(2002)利用1982—1999年AVHRR-NDVI数据分析得出，黄河流域植被覆盖状况总体一直处于上升趋势； 袁丽华等(2013)研究了黄河流域植被2000—2010年NDVI的时空变化特征，发现植被改善区域远大于退化区域； 贺振等(2017)基于GIMMS-NDVI数据，研究了黄河流域1982—2013年植被覆盖演化趋势，表明植被覆盖呈现较缓慢的增长趋势，覆盖度从北至南逐渐提高； 颜明等(2018)基于GIMMS和MODIS的NDVI数据分析了黄河流域1982—2012年不同时间尺度植被覆盖的空间特征，发现总体特征为西北小、东南大； 胡必波等(2018)用MODIS-LAI产品分析黄河流域2000—2014年植被覆盖度变化，呈上升趋势，但不显著。

目前黄河流域植被覆盖变化研究中的植被覆盖包括了林地、草地、耕地、园地和其他土地等，鲜有只针对林地和草地的研究，且近几年黄河流域植被覆盖动态变化分析研究也较少。鉴于此，本研究以黄河流域2009—2018年SPOT/VEGETATION NDVI时间序列数据为数据源，采用Theil-Sen median趋势分析和Mann-Kendall检验方法，仅研究林地和草地植被覆盖的时空分布及变化特征，并用Hurst指数方法分析预测未来发展趋势，以期为流域林草植被科学保护和精准修复、促进生态系统正向演替提供参考依据。

1 研究区概况

黄河发源于青藏高原巴颜喀拉山北麓的约古宗列盆地，流经青海、四川、甘肃、宁夏、内蒙古、陕西、山西、河南和山东9省(区)，在山东省东营市垦利区流入渤海，干流长5 464 km，水面落差4 480 m。黄河流域(96°～119°E、32°～42°N)位于我国干旱、半干旱和半湿润地区，总面积79.5万 km2(含内流区4.2万 km2)，其中上游42.8万 km2(平均海拔2 563 m)，中游34.4万 km2(平均海拔1 199 m)，下游2.3万 km2(平均海拔124 m)。黄河流域属大陆性气候，东南部为半湿润气候，中部为半干旱气候，西北部为干旱气候。土壤类型主要有草甸土、褐土、栗钙土、棕钙土、灰钙土、黑垆土、黄绵土和风沙土等。流域内地貌类型多样，包括高山、中山、低山、丘陵和平原等。植被类型丰富，包括针叶林、阔叶林、灌丛和灌草丛、草原和稀树草原、荒漠等。黄河流域位置见图1。

图1 黄河流域位置示意图Fig. 1 Location of the Yellow River Basin

2 研究方法

2.1 数据来源

NDVI数据集采用中国科学院资源环境科学与数据中心的中国年度植被指数空间分布数据集(http://www.resdc.cn/)，该数据集基于连续时间序列的SPOT/VEGETATION NDVI卫星遥感数据，在每年1～12月份的月数据基础上采用最大值合成法生成，选取2009—2018年数据，空间分辨率为1 000 m×1 000 m，共计10幅影像。黄河流域范围矢量数据由水利部黄河水利委员会提供。林地、草地矢量数据来自国家林业和草原局2018年全国森林资源管理“一张图”成果，按土地利用类型提取黄河流域范围内的林地、草地图斑，融合生成林草范围矢量数据，林地和草地面积共507 206.74 km2，占黄河流域总面积的63.83%。

2.2 年际变化与空间格局分析

逐年计算2009—2018年黄河流域林草地NDVI年均值，得到10年数据，分析其时间变化特征； 逐像元计算2009—2018年林草地NDVI平均值，得到10年NDVI平均值空间分布，分析空间格局特征。

2.3 趋势变化分析

采取Sen趋势度结合Mann-Kendall趋势显著性检验的方法开展林草植被覆盖变化趋势分析，该方法已成功用于气象因子(Chenetal.， 2007； Sunetal.， 2014)及植被变化趋势分析(Xuetal.， 2003； 蔡博峰等， 2009)。

Theil-Sen Median趋势分析是一种稳健的非参数统计趋势计算方法，可以减少数据异常值的影响(Theil，1950； Sen，1968)。计算公式如下：

对于时间序列{xt}，t=1，2，3，…，n，定义斜率Q：

式中：xa和xb是时间序列数据； Median代表取中位数； 当Q>0时，时间序列{xt}呈现增长趋势，反之时间序列{xt}呈现下降趋势。

雾化法制取的粉末已占当今世界金属 3D 打印粉末的80%以上，其原理是以快速运动的流体(雾化介质)冲击或以其他方式将金属或合金液流破碎为细小液滴，随之冷凝为固体粉末的粉末制取方法，其原理结构图如图1所示[5]，根据雾化介质不同，雾化法主要分为水雾化和气雾化。

Mann-Kendall非参数检验方法用于判断趋势的显著性，其优点在于不需要样本服从正态分布，也不受缺失值和异常值的干扰(Mann，1945)。计算公式如下：

对于时间序列{xt}，t=1，2，3，…，n，定义统计量S：

S服从正态分布，其方差VS为：

VS=[n(n-1)(2n+5)]/18。

定义统计量Z：

式中：xi和xj为时间序列数据； sgn为符号函数；统计量Z的取值范围为(-∞，+∞)。在给定显著性水平α下，在正态分布表中查得临界值u1-2/2当|Z|>u1-α/2时，表示时间序列{xt}在α水平上存在显著的变化。通过0.05置信水平的显著性检验判断时间序列{xt}变化趋势，当Z>1.96时表明时间序列{xt}显著上升； 当0

2.4 未来演化趋势分析

Hurst指数是定量描述时间序列长期依赖性的有效方法，由水文学家Hurst(1951)提出，在水文学、经济学、地质学等领域有着广泛应用(王孝礼等， 2002； 彭丽娜等， 2006； 贺可强等， 2009)，近期逐渐应用于植被覆盖变化研究中(李杰等， 2019； 严恩萍， 2014； 袁丽华等， 2013)。计算方法如下。

(τ=1 ， 2 ，…，n)；

累积离差序列Xt,τ：

(1≤t≤τ)；

极差序列Rτ：

( 1≤t≤τ；τ=1 ， 2 ，…，n)；

标准差序列Sτ：

(τ=1 ， 2 ，…，n)；

若Rτ/Sτ∝τH成立，则说明时间序列{xt}存在Hurst现象，Hurst指数(H)值可在双对数坐标系[inτ,in(Rτ/Sτ)]中通过最小二乘法回归得到。

根据H的值判断时间序列{xt}是完全随机还是存在持续性。Hurst指数取值包括3种形式： 如果0.5

3 结果与分析

3.1 NDVI空间格局分析

基于2009—2018年NDVI时间序列数据，逐像元计算10年平均值，得到NDVI平均值空间分布(图2)，将其分成6个等级并按上、中、下游进行面积统计(表1)。结果显示，流域高植被覆盖度(NDVI≥0.7)林草地面积为214 728.12 km2，占林草地总面积的42.34%，多为乔木林与灌木林，呈聚集状分布于流域上游的祁连山、三江源东部、若尔盖草原，中游的秦岭北麓、黄土高原胡焕庸线以南地区，以及下游的济南市区周边区域； 中植被覆盖度(0.3

图2 2009—2018年黄河流域林草地NDVI平均值Fig. 2 NDVI mean value in forestland and grassland in the Yellow River Basin in 2009-2018

表1 2009—2018年黄河流域林草地NDVI均值空间分布Tab.1 Spatial distribution of NDVI mean value in forestland and grassland in the Yellow River Basin in 2009-2018

3.2 NDVI年际变化分析

2009—2018年，黄河流域林草地年均NDVI有一定波动，但总体上呈增加趋势，增速约为0.007 4 a-1； 最小值与最大值分别出现在2011和2018年，分别为0.552 9和0.651 6(图3)。分析发现，2011—2013年呈显著增加趋势， 2015—2018期呈稳定增加趋势，总体上这10年内黄河流域林草植被覆盖状况逐步改善。

图3 2009—2018年黄河流域林草地NDVI年际变化Fig. 3 Inter-annual variation of NDVI in forestland and grassland in the Yellow River Basin in 2009-2018

3.3 NDVI空间变化分析

根据Theil-Sen median趋势分析结果，按Q>0、Q=0和Q<0，将林草地植被变化分别划分为改善、稳定不变和退化3种情况。将置信水平0.05上的Mann-Kendall检验结果按∣Z∣>1.96和∣Z∣≤1.96划分为显著变化和不显著变化。将Theil-Sen median趋势分级结果和Mann-Kendall检验分级结果进行空间叠加，根据耦合结果将林草植被覆盖变化趋划分为显著改善、轻微改善、稳定不变、轻微退化和显著退化5类(图4)，并分流域上、中、下游统计各类面积(表2)。

图4和表2表明， 2009—2018年黄河流域林草植被覆盖变化趋势以改善为主，改善区域面积413 972.93 km2，占林草地总面积的81.62%，其中显著改善面积169 177.13 km2，占33.36%，主要集中在上游的祁连山林区、若尔盖草原，中游的六盘山林区、子午岭林区、黄龙山林区、吕梁林区和秦岭北麓。退化区域面积70 317.69 km2，占林草地总面积的13.86%，以轻微退化为主，其中流域上游的青海、宁夏、内蒙3省(区)占轻微退化区域面积的85.00%，主要集中在青海的玛多县和内蒙的乌拉特前旗、乌拉特中旗； 显著退化面积3 974.46 km2，仅占林草地总面积的0.78%。林草植被覆盖变化趋势表现为稳定的区域面积22 916.12 km2，占林草地总面积的4.52%。综上， 2009—2018年间林草高覆盖度区域植被覆盖变化趋势表现为持续改善，低覆盖度区域持续退化。

图4 2009—2018年黄河流域林草植被覆盖变化趋势Fig. 4 Change trend of forestry and grass coverage in the Yellow River Basin during 2009-2018

表2 2009—2018年黄河流域林草植被覆盖变化趋势空间分布Tab.2 Spatial distribution of change trend of forestry and grass coverage in the Yellow River Basin in 2009-2018

图5 黄河流域林草植被覆盖未来发展趋势Fig. 5 Future development trend of forestry and grass coverage in the Yellow River Basin aa: 持续性&显著改善 Sustainability & obvious improvement； ab: 持续性&轻微改善 Sustainability & slight improvement； be: 反持续性&显著退化Reverse sustainability & obvious degradation； bd: 反持续性&轻微退化 Reverse sustainability & slight degradation； ac: 持续性&稳定不变 Sustainability & unchanged； bb: 反持续性&轻微改善 Reverse sustainability & slight improvement； ba: 反持续性&显著改善 Reverse sustainability & obvious improvement； ad: 持续性&轻微退化 Sustainability & slight degradation； ae: 持续性&显著退化 Sustainability & obvious degradation； bc: 反持续性&稳定不变 Reverse sustainability & unchanged.

3.4 NDVI未来演化趋势

3.4.1 NDVI空间可持续性 经R/S分析，研究区Hurst指数均值为0.57，说明黄河流域林草植被覆盖存在时间序列持续性。分类结果显示，0.5

3.4.2 NDVI未来发展趋势 将NDVI空间变化趋势和空间可持续性特征进行空间叠加分析(图5)，将研究区内林草植被覆盖未来发展趋势分为良性方向、恶性方向、稳定不变和不确定4个方向(表3)。林草植被覆盖呈良性发展的区域面积占林草地总面积的61.89%，其空间分布与林草植被覆盖改善区域基本一致，其中持续性轻微改善区域占35.01%、持续性显著改善区域占23.36%； 反持续性退化区域占3.52%，该区域未来植被停止退化，转向良性发展。林草植被覆盖呈恶性发展的区域面积占林草地总面积的33.60%，其空间分布与林草植被覆盖退化区域基本一致，其中持续性轻微退化、显著退化区域分别占9.77%、0.57%； 反持续性改善区域占23.26%，主要集中在三江源地区和毛乌素沙地，该区林草植被覆盖虽有轻微改善，但持续性较差，未来林草植被覆盖将由改善转为退化。林草植被覆盖未来发展方向稳定不变和不确定的区域占林草地总面积的4.51%。

表3 黄河流域林草植被覆盖未来发展方向空间分布Tab.3 Spatial distribution of future development direction of forestry and grass coverage in the Yellow River Basin

4 讨论

本研究分析了黄河流域林草植被覆盖在2009—2018年的时空变化特征，与已有研究结果相比，林草植被覆盖改善区域面积比例明显提高(贺振等， 2012； 2017)，预测未来植被持续改善区域面积继续增加(袁丽华等， 2013)，持续退化面积逐渐减少，说明流域整体生态环境在不断改善。这主要得益于林业重点工程的实施， 2006—2015年黄河流域人工造林68 013.00 km2(国家林业局， 2006—2015)，其中天然林保护工程、三北防护林工程和退耕还林工程分别占25.52%、23.22%和15.74%。根据历年全国人工造林综合核查结果(国家林业局， 2006—2015)，黄河流域各省(区)人工造林平均保存率87.40%，平均成林率41.00%，测算2006—2015年人工造林成林面积24 371.90 km2，可为流域森林覆盖率增长贡献3.07%。然而，黄河流域尚有33.60%的区域面临植被退化风险，主要为三江源地区和毛乌素沙地荒漠地区，应作为重点关注对象。

森林资源管理“一张图”是森林资源保护、经营管理的基础，为森林资源保护管理及生态文明建设重大决策提供科学依据。由于数据的可获得性，本研究仅采用2018年森林资源管理“一张图”的林地和草地范围，分析黄河流域2009—2018年植被覆盖变化，忽视了期间林地与草地范围的增加、减少以及相互转化，今后应进一步收集相关数据，补充林草地范围变化的研究，以使结论更为科学客观。植被覆盖变化受气候、人为等多重因子影响，本研究仅结合人为因素讨论了黄河流域林草植被的改善情况，未能分离出气候变化的影响。此外，对于植被退化区域的驱动机制也未能做出解释，下一步应基于森林资源连续清查、森林资源管理“一张图”成果，结合遥感、大数据、模型等技术，回溯退化区域植被覆盖变化过程，分析植被覆盖变化与气候条件、土地利用以及其他人类活动等因素之间的相关性，探究植被覆盖变化原因，为后期流域生态系统保护修复和综合治理提供科学依据，从而提升生态环境可持续发展能力，更好地满足流域高质量发展要求。

5 结论

本文基于2009—2018年SPOT/VEGETATION NDVI时间序列数据，采用Theil-Sen Median趋势分析和Mann-Kendall检验以及Hurst指数方法，对2009—2018年黄河流域林草植被覆盖变化进行了时空特征分析和未来发展趋势预测。主要结论如下：

1)黄河流域林草植被覆盖整体呈现东南部和西北高山地区高、北部低的特征。高覆盖区聚集在上游的祁连山、三江源东部、若尔盖草原，中游的秦岭北麓和黄土高原胡焕庸线以南地区； 低覆盖区则分布在毛乌素沙地、陇中黄土高原、阴山南麓及三江源北部。

2)2009—2018年间黄河流域林草植被覆盖显著改善，年均NDVI增速为0.007 4 a-1，植被覆盖改善区域面积占81.62%，远大于退化区域面积，改善区域主要为上游的祁连山林区、若尔盖草原，中游的六盘山林区、子午岭林区、黄龙山林区、吕梁林区及秦岭北麓。

3)黄河流域林草植被覆盖未来发展前景整体良好。林草植被覆盖良性发展的区域占林草地总面积的61.89%； 恶性发展区域面积占33.60%，主要集中在三江源地区和毛乌素沙地，该区植被覆盖在2009—2018年间虽轻微改善，但持续性较差，未来可能呈退化趋势，在流域生态保护治理中应重点关注。