1998-2018年山东省植被覆盖变化及其水文响应规律研究

庄会波，季 妤，高振勇，谭秀翠

（1.山东省水文局，济南250002；2.山东农业大学水利土木工程学院，山东泰安271018）

植被作为陆地生态环境中的重要组成部分，具有明显的时空变化特点，其覆盖变化可以反映区域内自然因素与人为因素对植被的影响效果［1，2］。植被变化改变流域下垫面条件，并通过影响降雨截留、蒸散发、降水入渗、地下水补给等水文过程，进而直接影响径流的形成与演变［3，4］，因此，分析植被覆盖的时空变化规律，及其对径流过程的影响，有利于阐述流域下垫面变化对水文循环过程的影响机制。毕早莹等［5］基于Budyko 理论定量分析窟野河流域植被变化对径流的影响，研究结果表明植被变化对径流的贡献率达到84.9%，而气候变化和人类活动影响比较小。李慧赟等［6］借助改进的水文模型模拟植被变化后的径流过程，在Crawford 流域有25%的面积植树造林情况下，流域年径流量减少32.4 mm，由植被变化引起的径流变化量占63%，由气候变化引起的径流变化量占37%。一般认为植物覆盖通过增加渗透到土壤基质水分，进而减少地表径流［7］，但这并不适应所有地区，所以，植被覆盖对径流产生怎样的水文效应还需要进一步研究。丁婧祎等［8］研究显示NDVI与径流的相关系数在0.04～0.59 之间，13 个子流域NDVI与径流间的正负相关性并不一致。在不同的时空尺度下，植被覆盖对径流的影响规律也存在差异。李春晖等［9］研究黄河流域降水、植被对径流影响，结果表明，在年内，NDVI变化与降水、径流变化呈明显的正相关，在年际，NDVI波动趋势与降水具有一定的相关性，但与径流、径流系数变化的关系相对复杂，不具有明显的相关性。董喆等［10］研究表明，径流量受降水与NDVI的共同影响，其中7-10月的降水与NDVI对径流的影响显著。史晓亮等［11］研究发现流域尺度上植被NDVI与径流不具有明显的相关性，但从像元尺度来看，植被NDVI和径流的正相关和负相关共存，流域不同空间位置的植被变化与径流的关系并不一致。

归一化植被指数NDVI（Normalized Difference Vegetation Index）作为植被生长状况及植被覆盖度的最佳指示因子，被认为是监测植被和生态环境变化的最有效指标［12］。文中采用归一化植被指数NDVI为植被覆盖表征指标，分析山东省植被覆盖的时空变化规律，与水文要素进行相关分析，并采用Copula函数构建径流与降水量、NDVI的二维联合分布，系统分析径流与影响因素之间的概率分布关系，研究结果可为水资源合理开发利用及下垫面因素影响分析提供技术指导。

1 研究区域

山东省位于中国东部沿海、黄河下游，境内中部山地突起，西南、西北低洼平坦，东部缓丘起伏，山东省地形分布见图1。山东水系比较发达，自然河流的平均密度每平方公里在0.7 km以上。湖泊集中分布在鲁中南山丘区与鲁西南平原之间的鲁西湖带。气候属暖温带季风气候类型，年平均气温11～14 ℃，降水季节分布很不均衡，全年降水量有60%～70%集中于夏季。由2000-2018年《山东省水资源公报》数据，年均降水量675.99 mm，径流深118.53 mm，各年降水量与径流深分布过程线见图2。

图1 山东省地形图Fig.1 Topographic aap of Shandong Province

图2 2000-2018年山东省降水量与径流深分布Fig.2 The distribution of precipitation and runoff depth in Shandong Province from 2000 to 2018

山东省土地利用现状数据由资源环境科学与数据中心（http：//www.resdc.cn/）获取，数据格式为1 km 的栅格数据。2015年，山东省土地利用类型以耕地为主，占总面积的65%，其次为建设用地，占15%，草地、林地、水域分别占比8%、6%和4%，其空间分布见图3。

图3 2015年山东省土地利用数据Fig.3 Land use data of Shandong Province in 2015

归一化植被指数NDVI能反映地表植被覆盖变化状况，是评价植被生长状况的重要指标［13］。文中采用的NDVI数据来源于中国年度植被指数（NDVI）空间分布数据集［14］。1998-2018年山东省平均NDVI为0.71，植被指数NDVI空间分布见图4。

图4 1998-2018年山东省植被指数NDVI空间分布Fig.4 Spatial distribution of vegetation index NDVI in Shandong Province from 1998 to 2018

2 研究方法

2.1 NDVI趋势分析

一元线性回归趋势分析法用于计算每个栅格的变化趋势，反映不同时期NDVI变化趋势的空间特征［15］，计算公式如下所示：

式中：i为年序列号；n为总年数，取值21（1998-2018年）；NDVIi为第i年NDVI值；slope为各栅格NDVI变化趋势的斜率。

2.2 突变检测

曼-肯德尔（Mann-Kendall）法是一种非参数统计检验方法［16］，被广泛应用于气候［17］、水文［18］、植被［19］等变量的突变分析中。

对于有n个样本量的时间序列x，构造一个秩序列：

在时间序列随机独立的假定下，定义统计量：

式中：UF1=0，E（sk）和var（sk）是sk的均值和方差，在x1，x2，…，xn相互独立，且有相同连续分布时。

给定显著性水平α，若|UFi|>Uα，则表明序列存在明显的趋势变化。该方法计算简便，可以明确突变开始的时间，并指出突变区域，因此是常用的突变检测方法。

按时间序列的逆序xn，xn-1，…，x1，再重复上述过程，同时使UBk=-UFk（k=n，n-1，…，1），UB1=0。

2.3 Copula联合分布

Copula 函数是定义域为［0，1］均匀分布的多维联合分布函数，它可以将多个随机变量的边际分布连接起来构造联合分布［20］。二维联合经验分布，计算公式如下：

式中：P为X≤xi，Y≤yi的二维联合概率值；Nmn为X≤xi，Y≤yi的序号；N为总的数据对数目。

对应的二维联合重现期计算公式为：

在水文计算中，二维联合理论分布函数，一般采用最多的是Archimedean Copula 函数，其中Clayton Copula、Gumbel Copu⁃la、Frank Copula最为常用。

Clayton Copula函数：

Gumbel Copula函数：

Frank Copula函数：

为实现对Copula函数的拟合优度检验，采用、AIC信息准则和离差平方和最小准则（OLS）来评价Copula 方法的有效性［21］，选择AIC值越小、OLS值越小的Copula函数构建二维联合分布。

AIC信息准则计算公式为：

式中：Pei、Pi分别为联合经验频率和联合理论频率；i为样本序号；n为样本容量；m为模型参数的个数。

AIC值越小，表明Copula函数的拟合效果越好。

OLS的计算公式为：

3 结果分析

3.1 NDVI特征分析

采用一元线性回归趋势分析法，计算1998-2018年山东省NDVI空间变化趋势，计算结果见图5，可以看出，大部分区域NDVI呈增加趋势，植被覆盖良好，其占总面积的63%（slpoe>0.001），NDVI基本不变的区域占13%（0.001≥slpoe>-0.001），而NDVI呈减小趋势的区域占24%（slpoe≤-0.001）。

图5 1998-2018年山东省NDVI变化趋势Fig.5 Trends of NDVI in Shandong Province from 1998 to 2018

由山东省NDVI突变检测结果可以看出，见图6，2004年之后，NDVI统计值UF>0，NDVI呈增加趋势，但未超过显著水平线|UF|

图6 NDVI突变检验Fig.6 The abruption tests of NDVI

由山东省不同时期土地利用数据可以看出，见图7，2000年至2015年，土地利用结构发生了一定变化，主要是由耕地、草地、未利用土地向建设用地转化，其面积增大2.2%。不同土地利用区域的NDVI值也有较大差异，由表1 可以看出，NDVI值相对较低的区域集中在未利用地、水域和建设用地，其中仅建设用地的NDVI变化趋势是轻度减少，主要受城市化扩张影响。

表1 不同土地利用类型区域NDVI变化趋势Tab.1 Variation trend of NDVI in different land use types

图7 2000-2015年山东省土地利用类型数据分析Fig.7 Data Analysis of Land Use Types in Shandong Province from 2000 to 2015

3.2 水文要素特征分析

根据2000-2018年《山东省水资源公报》中降水、径流深统计资料，进行突变检测。

由径流深突变检测结果可以看出，见图8，径流深|UF|

图8 径流深突变检验Fig.8 The abruption tests of runoff depth

由降水量突变检测结果可以看出，见图9，与径流深突变分析曲线分布比较类似，降水量|UF|

图9 降水量突变检验Fig.9 The abruption tests of precipitation

受土地利用类型与植被覆盖度的双重影响，2000-2015年，耕地、草地面积减少，但区域内NDVI增大，建设用地面积增加，区域内NDVI降低。而2000-2015年，降水量由607.3 mm，减少到575.7 mm，径流深由106.9 mm，减少到53.8。受降水与下垫面因素共同影响，2000年到2015年，山东省径流深减少50%。

3.3 NDVI水文响应分析

为分析NDVI与径流的变化关系，采用相同时间尺度下的统计数据，分析NDVI与径流的变化关系。由图10（a）可以看出，在年尺度下，径流深与NDVI呈正相关，即随着NDVI的增大，径流深增大，但两者的相关系数R2值较低，仅为0.24，而径流深与降水量的相关程度较为密切，两者的相关系数R2值为0.90。表明，在年尺度下，山东省的径流深主要受降水量的影响。

由径流深与NDVI、降水量三者的分布关系可以看出，见图10（b），随着NDVI的增大、降水量的增大，径流深呈增加趋势，径流深最小值为33.9 mm，对应的降水与NDVI均为最小值，分别为420.2 mm 和0.65，当径流深达到最大值222.9 mm 时，对应的降水量为936.3 mm，NDVI为0.71。

图10 2000-2018年山东省径流深与降水、NDVI相关关系Fig.10 The correlation between runoff depth，precipitation and NDVI in Shandong Province from 2000 to 2018

3.4 二维Copula联合分布

多变量联合分析比单变量分析能更好地描述水文事件的内在规律和分析各个特征属性之间的相互关系［22］。因此，本文利用Copula 函数分别构建径流深与NDVI、降水量的二维联合分布，该函数不需要假设变量是独立的或者正态分布的或者它们具有相同的边缘分布，可以用来描述水文变量之间的相关性结构，能够灵活地构造边缘分布为任意分布的水文变量联合分布［20］。

通过拟合优度检验，计算结果如表2所示，通过比较可以看出，Gumbel Copula 函数的AIC与OLS值均最小，拟合精度较高。由Gumbel Copula 函数的理论联合分布与经验联合分布图可以看出，见图11、图12，Gumbel Copula 联合分布的相关系数R2分别是0.86和0.97，表明选取Gumbel Copula函数构建二维联合分布是合理的。

表2 Copula函数的拟合优度检验Tab.2 The goodness-of-fit test for Copula function

图11 NDVI与径流深理论联合分布与概率联合分布Fig.11 The theory joint distribution and probability joint distribution of NDVI and runoff depth

图12 降水量与径流深理论联合分布与概率联合分布Fig.12 The theory joint distribution and probability joint distribution of precipitation and runoff depth

由径流深与NDVI的二维联合分布图可以看出［见图13（a）］，随着径流深与NDVI的增大，两者的联合概率值增大。由二维联合分布等值线图［见图13（b）］，可以确定不同径流深与NDVI组合的联合概率值，当径流深一定时，随着NDVI的增大，联合概率值增大；当NDVI一定时，随着径流深的增大，联合概率值增大。

Fig.13 径流深与NDVI二维联合分布与等值线图Fig.13 Two-dimensional joint distribution and contour map of runoff depth and NDVI

径流深与降水的二维联合分布反映了基本相同的变化规律，由图14可以看出，随着径流深与降水量的增大，两者的联合概率值增大。根据等值线图，可以确定不同径流深与降水量组合的概率值。

Fig.14 径流深与降水量二维联合分布与等值线图Fig.14 Two-dimensional joint distribution and contour map of runoff depth and precipitation

由径流深与NDVI、降水量的联合重现期及等值线图，可以看出（见图15、图16），随着各要素的增大，联合重现期增大，但最大不超过20年，且随着重现期的增大，径流深与NDVI、降水量的取值范围变小。联合重现期是表示为“或”的联合分布重现期，即径流深大于等于特定值或是NDVI、降水量大于等于特定值时的重现期情况。

Fig.15 径流深与NDVI联合重现期与等值线图Fig.15 Combined return period and contour map of runoff depth and NDVI

Fig.16 径流深与降水量二维联合分布与等值线图Fig.16 Combined return period and contour map of runoff depth and precipitation

4 结 论

本文分析了1998-2018年山东省植被覆盖时空变化规律及其对径流的影响，基于Copula 函数构建了径流与NDVI、降水的二维联合分布函数，通过以上研究得到如下结论。

（1）在空间上，山东省63%的面积NDVI增大，13%面积NDVI不变，24%的面积NDVI减小。

（2）在时间上，山东省NDVI呈增加趋势，但增加趋势不显著。

（3）径流深与NDVI、降水量呈正相关，在NDVI与降水量的共同影响下，2011年前后，径流深减少30%。

（4）采用Gumbel Copula 函数构建径流深与NDVI、降水量的二维联合分布，可以确定不同要素组合下的联合概率值及重现期。

植被覆盖变化受自然因素与人为因素共同影响，在空间上呈现差异性，其产生的水文效应就更为复杂。如何在空间上，有效的分析植被覆盖的水文响应，遥感技术提供了良好的技术手段，其可以快速获取大范围、高精度的水文资料，如降水［23，24］、径流［25］等，有效弥补了水文气象站点布设不足的缺陷。因此，可在后续的研究中，将综合运用实测资料与遥感卫星资料，提高数据分析的精度与可靠性，为水文响应研究提供有效手段。 □