当前, 黄土高原地区水土流失依然严重^[1], 生物多样性遭到破坏, 由于气候变化加剧^[2-3], 该地区面临极端水文事件频发带来的一系列环境问题^[4-7], 生态环境治理遭受严重挑战。随着我国社会经济的发展, 人民对美好生活更加向往, 生态环境保护刻不容缓^[8]。黄河流域是中华文明的重要诞生地, 在气候变化和人类活动的共同作用下, 流域生态环境与自然条件显著变化, 对当地人民的生活和生产带来巨大的影响, 黄河流域的健康和高质量成为一项重大国家战略^[9-10]。黄河中游地区主要位于中国黄土高原地区, 该区域水土流失极为严重, 生态脆弱, 环境恶劣, 是入黄泥沙的主要来源地^[11-14], 因此黄土高原地区的水土流失治理工作关系到黄河流域的健康发展。

黄河流域治理的关键问题是水沙变化^[15]。从20世纪50年代开始, 大规模的水土保持工程在黄土高原地区进行了实施, 包括林草措施等^[16]。特别是在1999年, 退耕还林政策在黄土高原地区进行了大规模实施^[17-20], 区域植被显著增加, 有效的减轻流域土壤侵蚀, 黄河流域泥沙显著减少^{[16, 21]}。植被在陆-气系统之间发挥着能量传输、维持和优化生态系统服务的重要作用, 是全球生态环境变化的重要“指示器”^[22-23]。随着全球气候变化与人类活动影响对区域发展影响的课题受到学术界的高度关注, 陆地生态系统与气候变化的关系研究成为当前国际研究的热点科学问题之一^[24-25]。由于归一化植被指数(NDVI)能够指示大尺度的区域植被覆盖状况, 表征植被活动的强弱, 因此常被广泛应用于研究生态变化的研究中^[26-27]。许多研究成果表明植被变化对气候因子变化具有明显的响应, 在区域尺度上气候条件和背景环境的空间异质性对植被变化具有决定作用^[28-29], 干旱区的植被净初级生产力(NPP)更容易在降水的影响下增加^[30]。虽然过去的研究已经取得了丰富的成果, 但是以往气候因子与植被变化的研究主要关注的是气候变化的趋势或者线性变化对植被的影响, 缺乏关于气候因子和植被之间二元关系变化的理论和模型等研究。有一个至关重要事实是, 由于过去对于植被变化的研究大多数是基于单变量统计学手段对其进行分析, 例如变异系数和Mann-Kendall趋势检验等, 对于植被与气候因子的响应关系多基于线性方法进行简单的研究, 例如线性相关检验、线性回归模型等等, 所以这容易造成变量之间的变化存在巨大的差异。自然界中的地表水和地下水, 植被, 大气, 土壤等环境要素之间互相作用, 非线性关系极为复杂^[31], 使用单变量的分析方法有利于简化模型, 减轻工作量。然而, 单变量的分析方法只关注了单个变量的变化(例如降水或植被), 这不符合现实的具体情况, 很难完全揭示对气候变化和植被变化之间的响应关系。因此, 有必要构建植被与降水之间的双变量耦合关系理论, 诊断植被与降水之间变异特征。

本研究的主要目的是, 通过Mann-Kendall趋势检验和Hurst指数分析黄土高原地区1998—2017年NDVI的变化趋势与持续性, 利用RDR指数与重心转移模型分析黄土高原地区NDVI的时空变化差异, 基于耦合协调度理论和Pettitt检验方法构建区域NDVI和降水关系的突变点的识别方法并进行诊断, 探讨了造成黄土高原地区NDVI和降水关系变化的可能驱动因素。研究成果有助于进一步深化对黄土高原地区植被变化和降水关系变化的认识, 有利于为黄土高原水土保持工作开展提供一定的科学依据。

1 材料与方法 1.1 研究区概况

黄土高原地区作为世界上黄土覆盖面积最大的高原(34°—42°N, 101°—114°E), 黄土平均厚度在50—80 m, 部分地区在100 m以上, 总面积大约为64万km²。黄土高原地区位于属于我国地势上的第二阶梯, 海拔大约81—4918 m, 平均海拔为1404 m。区域内主要包括甘肃省、宁夏回族自治区、青海省、内蒙古自治区、山西省及河南省等省区的46个地区(盟、州、市), 282个县。黄土高原地区河水的主要来补给源是降水, 该区域干旱缺水。黄土高原地区特有的气候特征和湿陷性黄土以及历史上的植被破坏, 使得该地区水土流失异常严重, 沟壑交错, 流域支离破碎, 生态环境脆弱。黄土高原地区的地貌分区大致可分为4种, 依次为北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区, 如图 1所示。

1.2 数据来源

本文所使用的降水数据是来源于国家地球系统科学数据中心的栅格降水量数据(http://loess.geodata.cn), 时间序列为1901—2017年, 空间分辨率为1 km×1 km, 时间分辨率为月, 该数据使用全国496个独立气象观测点数据进行了验证, 验证结果可信。本研究通过ArcGIS软件裁剪了黄土高原地区1998—2017年逐月的降水量数据, 并通过对逐月数据进行合成, 获得了黄土高原地区1998—2017年逐年年降水量栅格数据。植被归一化指数(NDVI)是一个长期的SPOT/VEGETATION NDVI卫星遥感数据集, 数据集由中国科学院地理科学与资源研究所网站提供(http://www.resdc.cn), 数据时间为1998—2017年, 空间分辨率为1 km×1 km, 时间分辨率为年。通过使用ArcGIS软件将NDVI数据与栅格降水量数据进行投影转换, 获得逐年相互匹配的NDVI数据和栅格降水量数据。

1.3.1 Mann-Kendall趋势检验与Hurst指数

Mann-Kendall(M-K)趋势检验^[32]作为一种非参数的时间序列趋势检验方法, 由于简单实用常被广泛使用在气象和水文序列的趋势检验中, 本文中使用Mann-Kendall趋势检验识别降水量、NDVI以及其他变量的变化趋势。

Hurst指数(H)常被用在定量表征时间序列的持续性或长期相关性^[33], 本研究中使用Hurst指数来描述NDVI的未来变化趋势。其中H代表Hurst指数。H的值在[0, 1]的范围内。当H=0.5时, 表明NDVI的时间序列是随机序列, 是不可持续的。当H>0.5时, NDVI的变化与目前的趋势基本一致, 表明植被的可持续性是正向的。H < 0.5表示负可持续性, 未来植被变化将与当前趋势相反。

1.3.2 RDR指数

相对发展率(RDR)被用来分析研究区内^[34], 不同栅格NDVI在研究时段初期与末期的变化量与研究区整体NDVI在研究时段初期与末期变化量的变化差异, 计算公式如下:

(1)

式中: Y_2i为2017年黄土高原地区第i个栅格NDVI；Y_1i为1998年黄土高原地区第i个栅格NDVI；Y₂为2017年黄土高原地区NDVI平均值；Y₁为1998年黄土高原地区NDVI平均值。

1.3.3 重心转移模型

重心转移模型^[34]常被用来描述地理变量分布重心在空间上的变化特征, 可以被用来研究NDVI空间重心随时间的变化, 重心转移模型在本研究中被用来研究黄土高原地区1998—2017年NDVI的重心变化, 计算公式如下:

(2)

(3)

式中: X_m和Y_m为NDVI空间分布重心的经度和纬度；c_mi为第i个栅格的NDVI；x_i为第i个栅格的经度；y_i为第i个栅格的纬度。

1.3.4 耦合协调度理论

两个或两个以上的系统或运动形式通过各种相互作用而彼此影响的现象被称为耦合^[35-37], 随着科学理论的深入研究与应用, 这一理论逐渐被广泛使用在生态学领域。NDVI和降水两个子系统之间的耦合度模型, 表示成以下公式^[38]:

(4)

式中: C为NDVI和降水二元系统的耦合度, u₁, u₂分别为NDVI子系统和降水系统对总系统有序度的贡献。在本研究中NDVI和降水子系统的指标为黄土高原地区所有栅格1998—2017的NDVI和降水量, 假定NDVI和降水量越大对系统越好。

(5)

(6)

式中u_i为第i个子系统对总系统有序度的贡献；u_ij为第i个子系统中第j个指标的归一化值；w_ij为第i个子系统中第j个指标的权重, 每个子系统中指标的权重计算使用熵权法进行计算。

在计算每个子系统的熵权时, 必须通过进行归一化先对数据处理, 常用的方法是极值法:

当u_ij越大对系统越有利(正向归一化):

(7)

当u_ij越小对系统越有害(负向归一化):

(8)

式中u_ij为子系统第i个指标的第j个时序的归一化值；x_ij为子系统第i个指标的第j个时序的值。

为了更好的反映出子系统整体“功效”与“协同”效应, 耦合协调度模型用来评判NDVI和降水两系统的交互耦合的协调程度, 计算公式如下:

(9)

(10)

式中, D为耦合协调度；C为耦合度；T为NDVI子系统与降水子系统的综合调和指数, 它反映NDVI子系统与降水子系统的整体协同效应或贡献；a、b为待定系数, 实际中常常认为两个子系统的重要性相同, 所以a=b=0.5。

1.3.5 Pettittt检验

Pettitt检验法^[39]采用Mann-Whitney中U_{t, n}值检验同一总体中两个样本X₁, …, X_t和X_t+1, …, X_N, Pettitt检验的零假设为没有变化点, 当|U_{t, n}|取最大值时对应的X_t被认为是可能的突变点。其显著性水平可由下式计算:

(11)

当P≤0.05时认为数据中存在均值变异点。

2 结果与分析 2.1 NDVI变化特征

黄土高原地区NDVI在1998—2017年73.49%面积呈现显著增加趋势(P < 0.05), 变化趋势时空分布如图 2所示。北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区的NDVI均有大面积呈现显著增加趋势(P < 0.05), 4个分区NDVI显著增加面积占各分区面积的比例为54.32%、72.91%、93.97%和74.98%(P < 0.05)。北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区的多年平均NDVI分别为0.34、0.53、0.55和0.58, 东部河谷及土石山区的植被覆盖状况最好, 北部风沙区和农灌区的植被覆盖状况最差(图 2)。1998—2017年黄土高原NDVI的Hurst指数空间分布最大值为0.75, 最小值为0.25, 黄土高原只有5.84%面积NDVI的Hurst指数低于0.5, 主要集中在北部风沙区和农灌区。黄土高原大部分区域NDVI在未来依旧呈现增加趋势。

2.2 NDVI时空变化差异

黄土高原地区NDVI在1998—2017年间变化存在时空差异, 如图 3所示。RDR指数最大值为8.14, 最小值为-6.42, RDR指数大于1的区域主要集中在丘陵沟壑区, RDR指数小于1的区域主要集中在北部风沙区和农灌区以及黄土高原边界区域。在高原沟壑区也有较大面积NDVI的RDR指数大于1, 这表明在1998—2017年间黄土高原地区NDVI在丘陵沟壑区与高原沟壑区的增加幅度远远大于黄土高原整体的增加幅度, 北部风沙区和农灌区以及黄土高原边界区域的NDVI增加相对总体来说更为滞后。

图 3的结果显示, 1998—2017年间黄土高原NDVI重心在不断发生转移。从1998—1999年、1999—2000年、2000—2001年、2001—2002以及2015—2016年转移距离比较长, 这表明这期间植被在黄土高原地区大范围增加。从总体看, NDVI重心由1998年的109.36°E 35.79°N变化为2017年的109.36°E 35.81°N, 呈现从波动变化的趋势, 表明黄土高原地区植被在年际间存在明显的变化差异。

2.3 NDVI子系统与的降水子系统权重分布变化

根据耦合协调度理论, 构建了黄土高原地区1998—2017年间NDVI子系统与降水子系统之间的二元耦合系统, NDVI子系统与降水子系统的构建指标分别为黄土高原地区1998—2017年每个栅格的NDVI时间序列和降水量时间序列, 使用熵权法计算了每个指标的权重如图 4所示。NDVI的权重分布最大值为90.6×10^-7, 最小值为4.7×10^-7, 北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区的NDVI子系统权重均值分别为17.9×10^-7、16.66×10^-7、14.29×10^-7和15.05×10^-7, 高原沟壑区和北部风沙区和农灌区的NDVI变化对黄土高原地区整体的重要性较大(图 4)。图 4表明北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区的降水子系统权重均值分别为13.88×10^-7、19.50×10^-7、13.95×10^-7和15.21×10^-7, 高原沟壑区和东部河谷及土石山区的降水变化对黄土高原地区整体的重要性较大。

2.4 NDVI与降水的耦合协调度变化及其突变点

根据耦合协调理论计算, 1998—2017年间黄土高原地区NDVI与降水的耦合协调度变化如图 5所示。由M-K趋势检验分析可知, NDVI与降水之间的耦合协调度呈现显著增加趋势, 1998年耦合协调度为0.61, 2017年耦合协调度为0.81, 增加幅度为32.79%。由Hurst指数分析可知, 1998—2017年间黄土高原地区NDVI与降水的耦合协调度Hurst指数为0.65, 未来依然呈现增加趋势, 这表明黄土高原地区NDVI与降水之间的耦合协调程度在持续增强。由于耦合协调度可以反映黄土高原地区NDVI与降水之间的相互作用关系, 因此通过Pettitt方法识别了NDVI与降水的耦合协调度变化的突变点。1998—2017年间黄土高原地区NDVI与降水的耦合协调度在2006年发生了显著突变(P < 0.05), 这表明黄土高原地区NDVI与降水之间的关系可能在2006年存在变异。

3 讨论 3.1 降雨对NDVI的影响

降水作为重要的气候因子, 对植被的生长和变化具有不可忽视的作用^{[34, 40-41]}。许多研究成果表明, 降水可以促进植被更好的生长和存活^[28]。图 6所示为1998—2017年黄土高原地区NDVI显著增加区域(P < 0.05)和不显著变化区域(P>0.05)逐年平均年降水量的时间变化过程。NDVI显著增加区域的降水在1998—2017年间均高于NDVI变化不显著的区域, 这说明降水确实在一定程度上影响着黄土高原NDVI的变化。但根据1998—2017年黄土高原地区NDVI和降水量的相关性分布图可知, 黄土高原大部分地区NDVI和降水之间不存在显著的相关关系, 部分面积NDVI和降水存在显著正相关关系(P < 0.05), 主要集中在丘陵沟壑区、高原沟壑区北部和东部河谷及土石山区北部。这说明, 降水确实在黄土高原地区的局部地区发挥了促进植被增加的作用。

3.2 造成NDVI与降水之间突变点的原因

植被的生长确实需要充足的水分供应, 但是许多研究表明黄土高原地区大部分区域年降水并未发生明显的增加趋势^[42], 造成黄土高原地区植被在1998—2017年间增加的主要驱动力应该是退耕还林、植树造林等水土保持工程^[21]。黄土高原地区1998—2017年多年平均年降水量为445 mm, 空间变化范围为109—878 mm, 北部风沙区和农灌区、高原沟壑区、丘陵沟壑区和东部河谷及土石山区四大区域的多年平均年降水量分别为271、483、459、546mm, 北部风沙区和农灌区的多年平均年降水量最小, 东部河谷及土石山区的多年平均年降水量最大(图 7)。根据M-K趋势检验分析, 1998—2017年黄土高原地区年降水的变化趋势如图 7所示。黄土高原地区年降水显著增加的区域主要集中在丘陵沟壑区和东部河谷及土石山区中部(P < 0.05), 这可能是该区域NDVI显著增加的一个重要原因, 但是黄土高原地区大部分区域的降水并不存在显著变化趋势。图 7显示了1998—2017年间黄土高原地区NDVI均值的变化过程。根据M-K趋势检验, 黄土高原地区NDVI均值在1998—2017年呈现显著增加趋势。由Pettitt检验可知, 黄土高原NDVI均值在2006年出现显著突变点, 这可能是造成NDVI和降水关系在2006年发生显著突变点的主要原因(7), 而1999年中国政府开始在黄土高原地区实施了退耕还林政策, 极大的改善了黄土高原地区的生态环境, 这可能是造成NDVI在2006年发生突变的可能原因, 这同时表明黄土高原地区整体NDVI发生显著突变相对于黄土高原地区实施退耕还林政策开始实施的时间滞后了7年, 植被恢复需要一个较长的时间周期。

4 结论

(1) 在1998—2017年, 黄土高原地区73.49%面积NDVI存在显著变化趋势, 大面积区域植被在未来依旧呈现增加的趋势, 东部河谷及土石山区的植被覆盖状况最好, 北部风沙区和农灌区的植被覆盖状况最差；

(2) 1998—2017年间黄土高原地区丘陵沟壑区与高原沟壑区的NDVI增加幅度大于黄土高原地区整体的增加幅度, 北部风沙区和农灌区以及黄土高原地区边界区域的NDVI增加滞后于整体, NDVI逐渐从东南向西北增加；

(3) 本研究构建了基于耦合协调度的NDVI降水关系突变点诊断方法, 识别出黄土高原地区NDVI降水关系在2006年存在显著突变点, 高原沟壑区和北部风沙区和农灌区的NDVI变化对黄土高原地区整体的重要性较大, 同时高原沟壑区和东部河谷及土石山区的降水变化在黄土高原地区占有重要地位。

(4) 降水对黄土高原地区NDVI的增加具有积极的促进作用, 在丘陵沟壑区、高原沟壑区北部和东部河谷及土石山区北部NDVI和降水存在显著正相关关系, 但是黄土高原地区大部分区域的降水并不存在显著变化趋势, 因此造成黄土高原地区NDVI与降水关系在2006年发生显著突变的主要原因应该是退耕还林等水土保持人类工程措施, 黄土高原地区整体NDVI发生显著突变相对于黄土高原地区实施退耕还林政策开始实施的时间滞后了7年, 植被恢复需要一个较长的时间周期。