植被是全球环境变化研究的重要环节，是链接大气与土壤的重要物质载体^[1-2]。植被分布与变化过程与全球气候格局关系密切，植被既受全球大气候背景的影响，也会对区域小气候产生作用^[3-4]。全球气候变暖大背景下，探讨不同区域植被对气候变化响应过程，对于科学认识陆地植被变化特征，辨识不同区域植被对气候变化的响应差异，以及指导生态环境政策科学制定具有重要意义。黄河是中华民族的母亲河，黄河流域是中华文明的发源地，随着黄河流域中部地区黄土高原退耕还林还草工程实施，黄河流域植被及河流径流量、输沙量发生了显著变化，对整个黄河流域产生影响^[5-8]。开展黄河流域植被对气候变化响应及其与水沙关系研究，对于科学认识黄河流域生态环境变化过程及政策调控具有重要理论与现实意义。

黄河流域植被时空变化及其水热响应、黄河水沙关系变化是多学科学者关注的热点问题^[9-10]。2019年9月18日，习近平总书记在“黄河流域生态保护和高质量发展”座谈会上强调“共同抓好大保护、协同推进大治理、让黄河成为造福人民的幸福河”^[11]，为黄河流域生态文明建设提出了要求，指明了方向。黄河流域开展的相关研究较多，主要涉及黄河流域总体、黄河源区、黄河中游黄土高原等区域，其中，已有研究认为黄河流域总体植被呈改善趋势，尤其是退耕还林还草工程实施以来，植被显著改善^[12-13]，但流域内部存在空间差异，如受气候变化影响，黄河源区气温升高，蒸发增强，对植被生长产生不利影响^[14-15]，而黄河中游流域黄土高原和下游植被变化更多受到人类活动影响，趋于改善^[16-17]。黄河水沙时间序列变化及其与植被变化关系等方面的研究开展较多，其中黄河径流量、输沙量锐减尤为受到关注^[18-19]，但水沙锐减与流域内植被变化间的物理机制尚未厘清^{[5, 8-9]}。

MODIS (moderate resolution imaging spectroradiometer) NDVI(normalized difference vegetation index)具有时间尺度短但空间分辨率高的特点，在三江源地区、黄土高原地区的研究中使用频率较高^[20-22]，而GIMMS (global inventory modelling and mapping studies) NDVI具有空间分辨率低但时间尺度长的特点，是开展大尺度、大范围研究的首选数据来源，应用也较为广泛^[23-25]。本研究以1982—2015年黄河流域GIMMS NDVI，气温，降水和2000—2015年黄河上中下游径流量、输沙量数据为基础，利用线性回归方法、相关分析方法等，开展黄河流域NDVI时空变化和对气温、降水的响应特征，及与径流量和输沙量关系研究，以期为该流域生态环境政策调整提供参考。

1 研究方法 1.1 研究区概况

黄河流域位于95°50'29″—119°6'26″E和32°6'53″—41°48'18″N之间，涉及青海、四川、甘肃省、宁夏回族自治区、内蒙古自治区、陕西、山西、河南和山东省，面积总计8.09×10⁵ km²。海拔高度0~6 255 m，其中流域西部为青藏高原，中部主要为黄土高原，东部为华北平原，地形自西向东由高原转为平原。黄河流域气候类型多样，黄河源区以高原气候为主，上游地区主要为温带大陆性干旱—半干旱气候，中下游地区为温带大陆性半干旱—半湿润气候。流域内部季节特征差异明显，温度东西向梯度特征高于南北向；降水集中在夏季，但分布不均匀，年际变化较大；总体上流域湿度小，蒸发大，无霜期短。

利用ArcGIS中水文分析模块，基于90 m分辨率STRM (shuttle radar topography mission) DEM (digital elevation model)提取黄河流域各小流域范围和河网，并以内蒙古自治区托克托县、河南省孟津县作为黄河流域上中下游的分界点，将各小流域进行归并，获得黄河流域上中下游范围，其中上游范围面积41.21×10⁴ km²，占流域总面积的50.94%；中游范围面积3.29×10⁵ km²，占40.66%；下游范围面积6.79×10⁴ km²，占8.40%。黄河流域上游地区尤其是青藏高原、高山地区，植被变化与气候因素关系密切；中游黄土高原腹地除受气候变化影响外，人类活动尤其是退耕还林还草工程的实施对植被恢复与重建效果明显；下游地区河南、山东省境内植被变化主要受到人类活动的干扰^[26-28]。

1.2 数据来源与处理

研究使用数据包括：①90 m分辨率STRM DEM数据，下载自地理空间数据云(http://www.gscloud.cn/)，用于黄河流域小流域提取及参与上、中、下游范围的划分。②1982—2015年15d 8 km分辨率GIMMS NDVI数据，下载自美国国家航空航天局(http://neo.sci.gsfc.nasa.gov/)。对每年24张影像，采用最大值合成(maximum value composition，MVC)方法获得年值序列数据。③1980—2015年1 km分辨率中国逐年平均气温、降水量空间插值数据集，该数据集是利用全国2 400多个气相站点日监测数据，应用澳大利亚ANUSPLIN插值软件生成。下载自资源环境数据云平台(http://www.resdc.cn/)。利用黄河流域边界裁切得到黄河流域年平均气温和年降水量空间插值数据，并重采样为8 km分辨率进行相关分析。④2000—2015年黄河流域《水资源公报》和《泥沙公报》，下载自水利部黄河水利委员会网站(http://www.yrcc.gov.cn/)。用于统计黄河及上中下游(头道拐、三门峡、利津水文站)年径流量、年输沙量数据。

1.3 研究方法 1.3.1 线性回归方法

线性回归方法用于分析1982—2015年NDVI时间序列和空间序列(像元尺度)变化过程及2000—2015年NDVI与径流量、输沙量线性拟合关系，计算公式如下^[29]：

$Y=a X+b $

(1)

在分析NDVI变化趋势时，式中Y为NDVI，X为年份；在分析NDVI与径流量和输沙量关系时，Y为径流量或输沙量，X为NDVI。a为线性回归方程系数，b为常数。a值计算公式如下^[29]：

$a=\frac{\sum\limits_{i=1}^{n} x_{i} y_{i}-\frac{1}{n}\left(\sum\limits_{i=1}^{n} x_{i}\right)\left(\sum\limits_{i=1}^{n} y_{i}\right)}{\sum\limits_{i=1}^{n} x_{i}^{2}-\frac{1}{n}\left(\sum\limits_{i=1}^{n} x_{i}\right)^{2}} $

(2)

在分析NDVI变化趋势时；x为年份(1982—2015年)；y为1982—2015年NDVI影像；i为样本数(i=1, 2, …, 34)。其中x₁=1982，…，x₃₄=2015；y₁为1982年NDVI影像，…，y₃₄为2015年NDVI影像。在分析NDVI与年径流量和输沙量关系时，x为年径流量或年输沙量(2000—2015年)；y为2000—2015年对应年份NDVI影像平均值；i为样本数(i=1, 2, …, 16)。其中x₁为2000年NDVI影像平均值，…，x₁₆为2015年NDVI影像平均值；y₁为2000年径流量或输沙量，…，y₁₆为2015年径流量或输沙量。a值同前。

1982—2015年NDVI变化趋势显著性采用F检验，计算公式如下^[29]：

$F=\frac{U}{Q /(n-2)} $

(3)

$Q=\sum\limits_{i=1}^{n}\left(y_{i}-y_{i}^{\prime}\right)^{2} $

(4)

$U=b^{2} \sum\limits_{i=1}^{n}\left(x_{i}-x^{\prime}\right)^{2} $

(5)

式中：x为年份(1982—2015年); y为1982—2015年NDVI影像; i为样本数(i=1, 2, …, 34)，n=34;y'为y拟合值; x'为x的平均值; x_i和y_i含义同前。

1.3.2 相关系数

用于分析1982—2015年黄河流域NDVI与年平均气温、年降水量的相关关系，计算公式如下^[29]：

$r=\frac{\sum\limits_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sqrt{\sum\limits_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}} \sqrt{\sum\limits_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}}} $

(6)

式中：r为相关系数；x为1982—2015年NDVI影像；y为年平均气温或年降水量，i为样本数(i=1, 2, …, 34)，其中x₁为1982年NDVI影像，…，x₃₄为2015年NDVI影像；y₁为1982年平均气温或降水量插值数据，…，y₃₄为2015年平均气温或降水量插值数据。

2 结果与分析 2.1 NDVI时空变化 2.1.1 NDVI随时间变化特征

黄河流域东西长南北窄，自西向东自然地理条件与社会经济条件差异极大，也造成流域内NDVI具有一定的差异。黄河流域1982—2015年NDVI平均值为0.520 1，上、中、下游流域NDVI平均值分别为0.470 1, 0.538 6和0.725 0，即NDVI平均值：下游流域>中游>流域总体>上游流域，表现出与降水自东向西逐渐减少的一致性，其变化趋势见图 1。黄河流域总体、上游、中游和下游流域在1982—2015年NDVI均表现为增加趋势，线性递增速率分别为0.001 6/a(p < 0.001)，0.001 0/a(p < 0.001)，0.002 6/a(p < 0.001)和0.000 7/a(p>0.05)，显著性检验结果表明，仅有下游流域增加趋势不显著。

2.1.2 NDVI随空间变化特征

在像元尺度上，将计算出的1982—2015年黄河流域NDVI变化速率划分为增加和减少2类，并将计算出的F检验结果根据阈值(4.15)划分为显著和不显著2类，最后将变化速率和F检验结果进行叠加，则可将NDVI变化趋势划分为4类，分别为显著增加、(不显著)增加、(不显著)减少、显著减少(见图 2)。

由图 2可知，1982—2015年黄河流域NDVI以增加趋势为主，占流域总面积的91.66%，其中显著增加面积占75.18%，集中分布在中游及其附近区域；显著减少面积仅占1.87%，集中分布在中游关中平原、上游青海省境内，下游河南省、山东省黄河河滩、上游内蒙古自治区境内也有分布。上游流域NDVI呈增加趋势面积总流占域面积的90.05%，其中显著增加趋势占70.37%；中游NDVI增加趋势占95.04%，其中显著增加趋势占85.33%；下游NDVI增加趋势占84.97%，其中显著增加趋势占55.09%。

2.2 NDVI与气温及降水的关系

气温、降水是植被生长发育的必须条件，充足的热量和水分条件有利于植被的生长。黄河流域自南而北热量递减，水分自东向西递减。黄河源区地处青藏高原区东北部，相对流域其他区域较湿润，黄河流域中北部和西北部(陕西北部、内蒙古南部、甘肃中北部、宁夏回族自治区等)由于降水相对较少，属于半干旱和干旱地区。以相关系数0值为界，将1982—2015年黄河流域NDVI与年平均气温、年降水量的相关关系划分为负相关和正相关，再以相关系数临界值(0.349 4)将负相关和正相关划分为显著负相关、不显著负相关、不显著正相关、显著正相关4类。1982—2015年黄河流域NDVI与年平均气温相关性分布(见图 3a)表明，流域总体上NDVI与年平均气温呈显著正相关的面积占22.39%，集中分布在陕西省北部、内蒙古自治区南部和宁夏回族自治区东部，其他地区有零散分布；显著负相关比例较低(占3.41%)。黄河流域上游、中游和下游地区NDVI与年平均气温呈显著正相关的面积分别为占21.23%, 25.97%, 11.88%。黄河流域NDVI与年降水量相关性分布(见图 3b)表明，总体上NDVI与年降水量呈显著正相关面积占21.99%，集中分布在流域中北部(陕西省北部、内蒙古自治区境内)；显著负相关比例不足1%，集中分布在流域的西部(青海省境内)和中南部(关中平原等)。黄河流域上游、中游和下游地区NDVI与年降水量呈显著正相关面积分别占22.09%, 25.72%, 3.25%。

2.3 NDVI变化与年径流量及年输沙量的关系

黄河流域水沙变化是多种因素共同作用的结果。本文仅从统计学角度探讨NDVI变化与年径流量、年输沙量关系，是对黄河流域植被变化与水沙关系的粗浅认识，无法揭示其物理过程。采用2000—2015年利津水文站数据表征黄河流域总体情况，头道拐水文站数据表征黄河流域上游情况，三门峡水文站与头道拐水文站差值表征黄河流域中游情况，利津水文站与三门峡水文站差值表征黄河流域下游情况。

利用黄河流域及上中下游NDVI年值和年径流量、年输沙量制作散点图，并进行线性拟合，结果见图 4。黄河流域及上游、中游和下游流域NDVI与年径流量之间均正相关关系，即随着NDVI增加，年径流量呈增加趋势。线性拟合表明，黄河流域NDVI增加0.001，则年径流量增加1.39×10⁸ m³(p < 0.1)(图 4a)，上游、中游和下游流域年径流量分别增加1.85×10⁸ m³(p < 0.001)、2.79×10⁷ m³和5.28×10⁷ m³(图 4b，4c，4d)，其中仅黄河流域总体和上游流域分别通过0.1和0.001水平显著性检验。

黄河流域及中游流域NDVI与年输沙量之间呈负相关关系，而上游和下游流域呈正相关关系。线性拟合表明，NDVI每增加0.001，黄河流域及中游流域年输沙量下降3.80×10⁵ t和1.94×10⁶ t(p < 0.1)(图 4e，4g)，上游和下游流域分别增加4.90×10⁵ t(p < 0.1)和1.65×10⁶ t(图 4f，4h)，其中仅有上游和中游流域通过0.1水平显著性检验。

3 讨论

植被变化主要受制于气候变化和人类活动。尽管总体上黄河流域及上、中、下游流域NDVI均呈增加趋势，但在空间上存在明显差异。NDVI表现为下降的区域主要集中在黄河源区、流域北部内蒙古境内、关中平原和下游河滩区，而中部大面积区域呈增加过程。黄河源区1982—2015年NDVI下降明显的区域(图 2)，其年平均气温多呈升高趋势，而年降水量多呈下降趋势，呈暖干化特征。气候的暖干化抑制了黄河源区植被的生长发育。近些年黄河源区暖干化趋势有所下降^{[15, 21]}，抑制作用有所减弱并表现出海拔梯度差异^[30]。气候变化对上游流域影响高于中游和下游流域。人类活动在中游和下游流域产生了截然不同的结果，其中黄河流域中部大面积区域的NDVI增加，与人类活动密不可分，退耕还林还草工程实施，促使大面积坡耕地转为林草地，极大的改善了下垫面状况^{[16-17, 22]}，人类活动对植被变化起到了积极作用，但随之引起的其他自然要素连锁反应需要深入开展相关研究，如下垫面改善引起的蒸散发变化、黄河水沙变化等问题；而NDVI下降的区域集中在流域北部内蒙古境内的呼和浩特市附近、关中平原的西安市附近、及下游河滩区的郑州市和济南市附近区域，这些区域NDVI的下降与城市的快速扩张具有密切关系^{[20, 31]}。

黄河流域水沙变化尤其是水沙锐减是学者们一直关注的问题，但其物理机制尚不清晰^[9]。本研究涉及的NDVI与径流量、输沙量关系仅从统计学角度开展了分析，结果反映出黄河流域在总体植被改善情况下，黄河流域和上游流域年径流量呈增加趋势，而输沙量在上游流域表现为增加趋势，在中游流域表现为下降趋势。黄河流域内部植被变化与径流量、输沙量变化关系的空间异质性反映出黄河流域水沙变化驱动机制的复杂性。在植被与水沙变化方面开展的相关研究表明植被改善对产沙影响高于产流^[32]，退耕还林还草工程引起的植被覆盖变化和水土流失治理工程实施，削弱了流域水沙动力关系^[33]，黄河中游植被增加引起径流系数下降^[17]、河龙区间输沙锐减半数由植被恢复所引起^[34]，而黄河流域植被覆盖达到25%时具有最高的水土保持效益^[19]，这些研究反映出植被变化与径流量、输沙量变化之间具有一定的关系，但仅从统计学角度无法认识这种复杂关系，需要利用水文模型开展黄河流域植被变化与水沙变化关系研究^[35-36]。

4 结论

(1) 1982—2015年黄河流域NDVI平均值表现为：下游>中游>流域总体>上游流域。黄河流域及上、中游流域NDVI均呈显著线性增加趋势，其中中游流域增速最快，高于流域总体增速，达到0.002 6/a。黄河流域NDVI显著增加区域面积占75.18%，集中分布在中游及其附近区域；上游、中游和下游流域NDVI显著增加区域面积分别占流域总面积的70.37%，85.33%和55.09%。

(2) 1982—2015年黄河流域NDVI与年平均气温、年降水量均以正相关为主，其中显著正相关分别占22.39%和21.99%，集中分布在中北部区域(陕西省北部、内蒙古自治区南部和宁夏回族自治区东部)。上、中和下游流域中，以中游流域NDVI与年平均气温和年降水量的显著正相关比例最高，分别达到25.97%和25.72%，均高于黄河流域总体。

(3) 随着黄河流域植被总体改善，2000—2015年黄河流域NDVI与年径流量、年输沙量关系表现出明显区域差异。黄河流域及上游流域年径流量表现出随NDVI增加而增加的趋势，而上游和中游流域年输沙量随NDVI增加的表现趋势不一致，上游流域表现为输沙量增加趋势，而中游流域表现为下降趋势。黄河流域水沙变化是多种因素共同作用的结果，单纯进行植被指数与水沙关系的统计学分析，不能揭示其物理过程机制，需要深入开展相关研究工作。