1982—2015年永定河山区植被变化及对天然径流的影响

杜 勇，李建柱，牛凯杰，冯 平，胡庆芳，郑彦辰

（1. 天津大学 水利工程仿真与安全国家重点实验室，天津 300072；2. 永定河流域投资有限公司，北京 100193；3. 南京水利科学研究院 水文水资源与水利工程科学国家重点实验室，江苏 南京 210029）

1 研究背景

植被作为陆地生态系统的重要组成部分，具有连接大气、土壤、水文等生态要素的功能，在自然生态系统和人类生产生活中发挥着不可替代的作用［1-3］。归一化植被指数（Normalized Difference Vegeta⁃tion Index，NDVI）是公认的反映大尺度地表植被覆盖和生长状况的有效指标，与光合有效辐射、绿叶密度、植被生产力以及累积生物量等均具有密切的关系，通常用来表征植被活动的强弱［4-6］。近年来，由于气候变化和人类活动的影响，区域乃至全球尺度的植被均呈现显著增加趋势［7-10］，植被变化不仅能影响气候，而且可以通过改变植被蒸腾、冠层截流和土壤蒸发等水文过程而改变河川径流，故能显著影响流域水文循环过程［11-14］。因此，分析植被的时空演替特征，揭示在不同时空尺度下导致植被变化的驱动因子的空间异质规律，并解析以植被变化为重要特征的下垫面变化对径流变化的贡献及影响，可为应对气候变化和人类活动下生态系统的适应能力提供理论依据，同时有助于深刻理解植被活动对水循环过程的影响，并对完善流域水资源管理措施具有重要理论和实践意义。

作为植物生长发育所必需的外界环境因子，气候要素（水分和温度等）对植被活动起重要作用，同时人类活动对植被活动也具有不同程度的促进或抑制作用［15-18］。以黄土高原为例，已有研究表明过去30年气候变化对黄土高原的植被恢复有一定促进，同时植被恢复工程等人类活动也有利于NDVI的增加，因此气候变化和人类活动共同驱动黄土高原的植被变化，但两者的驱动作用呈现出显著的空间异质性［2，19-20］。同时大量研究证实，植被覆盖作为下垫面的重要性质，其变化将导致水循环过程变化，影响河川径流过程［21-23］。目前有关植被变化等研究已在黄土高原等生态脆弱区和人类活动密集区全面开展，但这些研究多采用相关性分析等方法对植被变化的动因进行定性研究，关于不同时间尺度下各种驱动因素影响的空间异质性的定量解析还不够深入［24-25］，且当前关于NDVI 变化对天然径流量影响的文献较少，不能清晰地从源头理解植被变化对河川径流过程的影响。

本文以黄土高原上的典型流域永定河山区为代表，基于气温、降水、径流以及长时序NDVI 数据集，通过构建地理加权回归模型，研究流域内NDVI 和降水、平均气温的回归参数、残差均值及其趋势的空间异质性，并进一步探讨以植被变化为突出特征的下垫面变化对流域典型断面天然径流的影响。

2 研究方法

2.1 NDVI 与影响因素的空间趋势分析方法NDVI 及其影响因素的空间趋势分析可采用基于栅格尺度的普通最小二乘法（OLS），通过模拟每个栅格NDVI 及其影响因素在相应时段的变化趋势，从而反映研究区在不同研究时段内的NDVI 及其影响因素空间特征。空间趋势分析的计算公式为：

式中：θslope为线性回归的斜率，表示不同时间尺度下NDVI、回归参数、残差等研究对象在流域或像元内的年际变化趋势；n 为研究时段年份数；Yi为第i年研究对象数值。

2.2 NDVI 与影响因素的空间异质性表征方法地理加权回归（Geographically Weighted Regression，GWR）是由Brunsdon 等［26］提出的一种空间分析方法，有助于揭示一定的研究区域内自变量与因变量之间的非平稳空间相关性。GWR 模型是对普通线性回归的重要拓展，该方法的参数是空间位置的函数，该模型的表达形式为：

式中： yi、 xik、 εi分别为空间上i 点的因变量、自变量和模型残差；(ui，vi)为i 点的空间位置；m为样本量i 点处自变量的个数，本文取2，分别表示降水量和气温； βk为i 点上的回归系数； β0为截距。在本文中，年均NDVI 为因变量，降水、年平均气温为相应自变量，用以分析NDVI 对气候要素的空间响应关系。

GWR 模型采用局部加权最小二乘法求解，其中的权重是待估点到各观测点的空间距离的函数，具有距离衰减效应。参数公式为：

式中：β( ui，vi)为回归系数的无偏估计；W (ui，vi)为权重矩阵，一般采用高斯模型作为权重函数，具体参见文献［27］；X 和Y 分别为自变量和因变量的矩阵。

本文通过建立GWR 模型得到年降水量和年平均气温对应的回归参数，并计算得到NDVI 的模拟值（NDVIcc），用来表征气候因素对NDVI 的影响；然后计算NDVI 实际观察值（NDVIobs）与模拟值（NDVIcc）之间的差值，表示人类活动对NDVI 的影响：

2.3 植被变化对径流影响的研究方法流域水量平衡方程为：

式中：R 为径流深，mm；P 为降水量，mm；ET 为实际蒸散发量，mm； ΔS 为流域储水量的变化，mm。当研究时间尺度较大时， ΔS ≈0，本文为了减少永定河山区流域ΔS的影响，取各变量时间尺度为10年的滑动平均值。

在流域尺度上，降雨量可直接通过雨量站点观测数据进行空间插值获得，实际蒸散发可基于Bu⁃dyko 假设的水热耦合平衡方程计算：

式中：ET0为年均潜在蒸散发量，mm；n 为反映流域下垫面特征的参数（简称下垫面参数），与地形、土壤和植被等一系列因素有关。对于流域下垫面参数n 的估计，诸多学者提出了不同的估算方法，这些方法均认为植被是影响n 的重要因素。Li 等［27］针对全球典型大流域研究分析发现，参数n 与NDVI存在较好的拟合关系：

式中：a、b 为回归参数，对特定流域为常数，a 的大小表征了n 对NDVI 变化的敏感性程度，b 则反映了NDVI 以外的其它因素对n 的影响。

3 研究区域和数据

3.1 研究区域永定河流域位于海河流域，地理位置介于东经111.95°—116°22′，北纬38°90′—41°16′，流经内蒙古、山西、河北、北京、天津五个省级行政区，集水面积约4.7 万km2，占海河水系的14.7%。2016年流域总人口约813 万，城镇化率为58.3%。永定河流域上游一级支流为桑干河和洋河，两支合并后始称为永定河。永定河流域属半湿润半干旱季风气候，多年平均降水量介于360 ~650 mm，6—8月降水量约占年均降水量的80%。干流官厅水库以上流域属于山区，面积为4.51 万km2，占全流域总面积95.8%，是主要产流区。册田水库、响水堡水库分别为永定河一级支流桑干河、洋河的主要水库（图1），集水面积分别为1.71 万km2和1.45 万km2。

图1 永定河山区气象水文站点和河道断面分布

3.2 数据来源与处理

3.2.1 NDVI 数据 本文收集了1982—2015年的永定河山区逐月NDVI 数据，数据来源于GIMMS ND⁃VI3g 数据集（https：//ecocast.arc.nasa.gov/data/pub/gimms/），空间分辨率为8 km，根据研究区域空间范围，提取了相应网格单元逐月NDVI 数据，经进一步处理后得到季尺度和年尺度的NDVI 数据。

3.2.2 气象数据 本文收集了1982—2015年的MSWEP 逐日降水数据和中国区域地面气象要素驱动数据集中的逐日平均气温数据。两者通过网络共享获取（网址分别为：http：//www.gloh2o.org/和http：//www.tpdc.ac.cn/zh-hans/），空间分辨率均为0.1°×0.1°。同时，在流域内可获取国家气象科学数据中心（http：//cdc.cma.gov.cn/home.do）提供的7 个国家基本气象站点1982—2015年逐日常规地表气象数据（含雨量和气温）。

对永定河山区0.1°×0.1°的逐日MSWEP 网格降水数据和中国区域地面气象要素驱动数据集中的平均气温数据，采用IDW 方法，分别采用国家基本气象站的降水与气温数据进行校正，之后对校正后的逐日降水、平均气温网格数据进行重采样处理，得到与NDVI数据空间分辨率一致的逐日降水、平均气温格网数据，最后通过时间尺度上的聚合处理，得到月、季、年等尺度的降水、平均气温数据。

3.2.3 地表径流数据 根据海河流域水文年鉴和各省市水资源公报，收集了册田水库、响水堡水库和官厅水库所在断面1956—2016年的实测径流量，采用分项还原法，对天然年径流量计算有影响的水量进行还原，将被人类开发利用的水量（如农业用水、工业用水、生活用水、环境用水耗损量、跨流域引调水量等）逐项还原到天然年径流量中，得到永定河典型断面的天然径流量。

4 结果分析与讨论

4.1 NDVI 的时空变化特点图2 给出了1980年和2015年研究区域土地利用类型分布，从中可知永定河山区主要土地利用类型是耕地、林地、草地、水域、建设用地及未利用土地，其中又以耕地、林地和草地为主。2015年较1980年林地面积增加了2364.4 km2，相对增幅34.6%；耕地、草地面积分别减少1610.1 和1056.8 km2，相对降幅分别为7.7%和7.5%；建设用地增加了926.9 km2，相对增幅116.2%。

图2 永定河山区1980年和2015年土地利用类型

表1 为永定河山区土地利用面积转移矩阵，反映出1980—2015年土地利用变化的主要特征是耕地、林地和草地之间的相互转移以及部分耕地、草地转化为建设用地；其结果是，流域内耕地占总面积的比重由1980年的47.3%下降至2015年的43.7%，林地和草地面积之和占总面积的比重由47.3%上升至50.2%，而建设用地占总面积的比重由1.8%上升至3.9%，不同类型土地利用变化量占流域总面积的比重较小且变化范围介于0.4% ~ 5.3%之间，因此，总体而言，1980—2015年期间永定河山区土地利用格局变化不大。

表1 1980—2015年永定河山区土地利用变化转移矩阵 （单位：km2）

图3 给出了永定河山区1982—2015年年均NDVI 变化，从3 可知，年均NDVI 呈波动增加趋势。在研究时段内，永定河山区年均NDVI 的波动范围是0.2722 ~ 0.3394，最小值和最大值分别出现在1985年和2013年。总体上，1982—2015年年均NDVI 线性趋势率为0.0013/a（p<0.01），这说明表明永定河山区植被恢复较明显。图4 表明，1982—2015年永定河山区NDVI 的年际变化趋势具有一定的空间异质性；年均NDVI 呈增加和减少趋势的区域面积分别占总面积的95.36%和4.64%。其中，年均NDVI 增长较快（slope≥2×10-3/a）的面积约占研究区域总面积的12.32%，主要分布在桑干河上游南部林地和永定河官厅水库附近的湿地；年均NDVI 呈降低趋势的区域分布比较零散，主要在城市扩张的地区，且降低速度率较小，多在-1×10-3/a 以内。

图3 1982—2015年永定河山区年均NDVI年际变化

图4 1982—2015年永定河山区年均NDVI 趋势空间分布

图5 给出了永定河山区1982—2015年春、夏、秋季的NDVI 变化情况，从中可知春、夏、秋季平均NDVI 均呈波动增加趋势。在研究时段内，永定河山区春季平均NDVI 的变化范围是0.2196 ~0.2746，最小值和最大值分别出现在1985年和2014年，其趋势率为0.0013/a（p<0.01）；夏季平均NDVI的变化范围是0.4317 ~ 0.5512，最小值和最大值分别出现在1985年和2013年，趋势率为0.0018/a（p<0.01）；秋季平均NDVI 的变化范围是0.2626 ~ 0.3669，最小值和最大值分别出现在1984年和2013年趋势率为0.0020/a（p<0.01）。这表明在春、夏、秋季永定河山区植被均有明显恢复。由图6 可知，1982—2015年永定河山区NDVI 变化趋势在春、夏、秋季均具有一定的空间异质性。其中，春季平均NDVI呈增加和减少趋势的区域面积分别占总面积的92.03%和7.97%， NDVI增长较快（slope≥2×10-3/a）的面积约占18.26%，主要分布在桑干河上游南部林地和永定河官厅水库附近的湿地地区；夏季平均NDVI呈增加和减少趋势的区域面积分别站总面积的90.58%和9.42%，NDVI增长较快（slope≥2×10-3/a）的面积所占比例高达46.40%，主要分布在桑干河中上游地区林地、洋河流域北部的林草地和官厅水库附近的湿地地区；秋季平均NDVI 呈增加和减少趋势的区域面积分别站总面积的98.41%和1.59%，NDVI 增长较快（slope≥2×10-3/a）的面积约占51.30%，主要分布情况与夏季类似。

图5 1982—2015年永定河山区春、夏、秋季平均NDVI年际变化率年际变化

图6 1982—2015年永定河山区春、夏、秋季平均NDVI 变化率的空间分布

因此，在土地利用格局变化不大的情况下，永定河山区NDVI 在年度和春、夏、秋季均有较大幅度的增加，可能由于受降水、气温的气象因素或流域内植被恢复工程、农林草地管理水平的提高等人类活动的影响，使年平均和各季度的NDVI 均有较大提高。

4.2 NDVI 与相关驱动因子的空间非平稳性

4.2.1 驱动因子回归参数与残差的均值 为验证植被NDVI 与降水、气温的空间相关关系是否具有非平稳性，首先对回归参数均值进行空间自相关分析，得到各要素GWR 回归参数的Moran 指数与Z 值如表2 所示。表2 结果表明，NDVI 与降水、平均气温的GWR 模型回归参数均值的Moran 指数均大于0，说明参数具有正的空间自相关，且Z 值均大于2.58，说明计算参数在0.01 显著水平下具有统计意义。因此采用GWR 方法分析降水、平均气温等要素在不同地理位置对NDVI 的定量影响是合理的。通过GWR 模型残差项可用来近似表征除降水、平均气温以外的其他因素对NDVI 的定量影响，包括地理位置、太阳辐射、人类活动等因素，结合流域的特点和前人研究［28］，人类活动影响最大，本文用残差体现人类活动的影响作用；另外，残差作为NDVI 实测值与模拟值之差具有一定的不确定性，这种不确定性主要来源于输入信息（如遥感数据、实测数据）的不确定性、模型结构的不确定性等方面，对结果影响不大。

表2 各要素GWR 回归参数的Moran 指数与Z 值

永定河山区1982—2015年年均NDVI 与降水、平均气温的GWR 空间回归参数及残差的均值分布如图7 所示。从图7 可知，在1982—2015年期间，总体上NDVI 与降水、平均气温表现为正相关关系，残差也为正值，这说明降水、平均温度、人类活动均对植被活动有一定的促进作用。同时，回归参数及残差的均值表现出显著的空间异质性。NDVI 与年降水量的GWR 空间回归参数在桑干河、洋河中上游地区介于5.9×10-4~ 7×10-4，表明这些地区年均降水量每升高100 mm 会导致NDVI 增加5.9×10-2~ 7×10-2；同时NDVI 对降水变化的敏感性从上游至下游逐渐加大，下游地区的GWR 空间回归参数提升至1.2×10-3~ 1.5×10-3。NDVI 与年均气温的GWR 空间回归参数介于0 ~ 7.2×10-3，其中约45%的区域低于1×10-3，主要分布在桑干河和洋河下游以及官厅水库以下地区。这说明流域内约一半的范围NDVI 对平均温度的变化不甚敏感，温度每升高1℃导致植被NDVI 的增加值低于1×10-3。GWR模型残差值在各像元内均为正，表明永定河山区人类活动对NDVI 为促进作用，且在流域周边地区残差值较高，这可能由于永定河山区周边地区主要为林地、草地，而中心区域主要为耕地，近年来人工植被恢复措施导致林草地NDVI 增长较耕地更快。

图7 1982—2015年永定河山区年均NDVI 与降水、平均温度的GWR 回归参数和残差均值空间分布

1982—2015年永定河山区不同季度的NDVI 与各影响因素的GWR 回归参数均值分布如图8 所示。从图8 可以看出，各季度GWR 回归参数分布与年平均情况下具有相似空间分布特征，且回归参数均为正值。

图8 1982—2015年永定河山区各季度NDVI 与降水、平均温度的GWR 回归参数和残差均值空间分布

NDVI 与降水GWR 的空间回归参数值在春、秋季大致介于2.5×10-3~ 7.5×10-3之间，NDVI 变化对降水敏感性最强的区域分布在永定河下游地区；而在夏季，GWR 回归参数值介于7.6×10-4~ 3.68×10-3，表明这一季节NDVI 对降水的敏感性程度相对较弱，洋河流域北部地区敏感性程度相对较高。春、夏、秋季NDVI 与年平均气温的GWR 回归参数大于2×10-3的面积占总面积的比重分别为27.4%、51.2%和61.2%，可以看出春季NDVI受温度影响的作用较夏、秋季明显偏弱；在空间上，各季度NDVI对温度变化较敏感的区域主要分布在桑干河、洋河流域上游。不同季节GWR 模型的残差项均为正，其中春、夏、秋季残差项大于0.05 的范围分别占流域面积的10.4%、22.5%、13.0%，主要分布在流域周边林草地和湿地，这说明人类活动在春、夏、秋季对NDVI 均为促进作用，且对夏季NDVI 受人类活动影响作用最大，秋季次之、春季最少。

4.2.2 驱动因子回归参数与残差的趋势 永定河山区1982—2015年年均NDVI 与降水、平均气温的GWR 回归参数及残差的时间变化趋势如图9 所示。从图9 可知，永定河山区NDVI 与年降水的GWR空间回归参数呈下降趋势的面积占研究区域总面积的比例接近17%，主要分布在洋河流域，而其余83%的面积对应的像元回归参数均呈上升趋势，这表明1982—2015年流域内多数地区NDVI 受降水的影响程度呈增加趋势；NDVI 与年平均气温的GWR 回归参数呈增加和下降趋势的像元比例为62.9%和37.1%，其中洋河流域和官厅水库以下区域主要呈下降趋势，而桑干河流域主要呈上升趋势，这说明洋河流域和官厅水库以下流域NDVI 受气温的影响程度在减弱，而桑干河流域NDVI 受气温的影响程度增强；残差呈增加和下降趋势的像元比例分别为39.1%和60.9%，表明人类活动对大多数区域NDVI的促进有所减弱。

图9 1982—2015年永定河山区年均NDVI 与降水、平均温度的GWR 回归参数和残差趋势空间分布

1982—2015年永定河山区不同季度下NDVI 与各影响因素的GWR 空间回归参数趋势分布如图10所示。降水方面，在春、夏季流域大部分地区（除东部外）NDVI 与降水的GWR 的回归参数斜率为正，而秋季NDVI 与降水的GWR 回归参数斜率在全流域内为负，表明在永定河山区大部分区域春、夏季NDVI 受降水的促进作用不断增加，秋季NDVI 受降水的促进作用在全流域内呈减小趋势；平均气温方面，春季流域内大多数区域平均气温与NDVI 的GWR 空间回归参数斜率为负，其中洋河流域回归参数斜率小于-20×10-5/a，表明春季流域总体上平均温度对NDVI 的促进作用呈减缓趋势，且以洋河流域降低程度最为显著；夏季平均气温与NDVI 的GWR 空间回归参数斜率为负的区域主要分布在洋河流域中下游地区，占总面积比例为31.3%，其余占流域68.7%的地区NDVI 受平均气温的促进作用呈增加趋势；秋季NDVI 受温度的促进作用呈减缓趋势的面积占44.2%，分布区域较夏季有所扩大。残差方面，春、夏、秋季残差趋势呈增加的面积占流域面积的比重分别为45.9%、40.7%和45.5%，且主要分布在永定河山区周边的林草地地区，而中部大部分区域NDVI 受人类活动的促进作用呈减弱趋势。

图10 1982—2015年永定河山区春、夏、秋季NDVI 与降水、平均温度的GWR 回归参数及残差的趋势空间分布

4.3 NDVI 变化对流域天然径流的影响人类活动对河川径流的影响主要包括两部分内容，其一主要为流域内人类取用水和水利工程调节的影响，其二主要为下垫面变化（土地利用和覆被变化）的影响［29］。对于天然径流，鉴于已对社会经济取用水和水利工程调节进行还原，天然径流变化可不考虑前一种影响，仅考虑下垫面变化的的影响。图11 为永定河山区1956—2016年典型断面（册田水库断面、响水堡水库断面、官厅水库断面）天然径流量，可以看出流域天然径流量的年际衰减趋势十分明显，这是气候变化、土地利用和覆被变化共同作用的结果。

为了进一步探究NDVI 变化对永定河山区天然径流的影响，首先对天然径流变异性进行了诊断，采用滑动T 检验法对各断面1956—2016年天然径流量进行了突变诊断（采用的显著性水平α=0.05，统计量临界值为1.67），诊断结果如图12 所示。由该图可知，册田水库、响水堡水库、官厅水库断面的滑动T 检验统计量分别在1982年、1983年、1983年取得最大值，且均显著高于临界值，这说明桑干河、洋河和永定河上游年际天然径流量在1982年、1983年、1983年发生了突变。为便于分析，本文统一将1983年作为3 个断面天然年径流量的变异点，同时下文将1956—1982年记作基准期，将1983—2016年记作变异期。结合图11分析，官厅水库断面基准期年均天然径流量43.4 mm，变异期为21.2 mm，相对减少51.1%；册田水库断面基准期、变异期年均天然径流量分别为37.8 和22.5 mm，相对减少40.5%；响水堡水库断面基准期、变异期年均天然径流量分别为51.0 和21.9 mm，相对减少57.1%。

图11 1956—2016年永定河山区主要断面地表天然径流量

图12 1956—2016年永定河山区控制断面天然地表径流量滑动T 检验统计量

图13 给出了永定河山区干支流断面在不同阶段的年降水量-天然径流量散点图及线性回归关系。从图13 可知，基准期内降水量与天然径流量有较好相关性，册田水库、响水堡水库和官厅水库断面降水-天然径流量线性回归方程的决定系数分别为0.57、0.63 和0.59，表明基准期降雨量对天然径流有较强的解释能力；而变异期内降水-径流量线性回归方程的决定系数变为0.17、0.07 和0.10，降水量对天然径流的解释能力明显变差。变异期年降水-天然径流量散点基本位于基准期散点下方，同等年降水量条件下永定河山区产流能力大幅度下降。根据两个时期的年降水-天然径流回归关系，计算得到变异期内由降水减少导致各断面天然径流累计缩减量分别为3.64、4.26 和3.34 mm，对天然径流减少的贡献率分别为16.4%、14.6%和21.8%。相关研究［30］表明，海河流域潜在蒸散发变化对径流的变化影响微弱，降水减少和由于人类活动导致的下垫面变化是流域径流量减少的主导因素。由于降水减少对永定河径流减少的贡献率仅为20%左右，因此导致永定河流域天然径流减少的主要原因是土地利用和覆被等下垫面条件的变化，由前述分析可知，永定河山区土地利用格局总体变化不大，而NDVI 呈显著变大的趋势，这进一步表明永定河山区天然径流衰减是由NDVI 的增加所引起的。

图13 永定河山区干支流断面年降雨量-天然径流量关系

通过式（5）与式（6）得到册田水库、响水堡水库、官厅水库控制流域内在不同年份的下垫面参数n，发现永定河山区三个水库控制流域n 的多年平均值分别为3.18、2.72、2.74，均大于傅抱璞模型推荐的n 的缺省值［31］，因此永定河山区的植被水分利用效率要高于平均水平，这是由于永定河山区处于中国半干旱地区，实际蒸散发量占降水量的比重高于90%，属于典型的水分限制区，因此参数n 的均值大于平均值符合水热平衡规律，且在各控制流域中，册田水库控制流域的持水能力最大，而响水堡水库控制流域的持水能力最小。同时发现，与天然径流的变化趋势相反，永定河山区下垫面参数n和年均NDVI 随时间呈增加趋势，其中册田水库、响水堡水库、官厅水库下垫面参数n 的年均增长率分别为0.0359、0.0426、0.0433，对应各流域年均NDVI 的增长率分别为0.0011、0.0013、0.0013，其中，响水堡水库与官厅水库的下垫面参数n 和NDVI 增长情况表现出高度一致性，均大于册田水库控制流域。绘制了不同年份滑动平均的下垫面参数n 与NDVI 点据分布，对下垫面参数n 与NDVI 进行拟合，如图14 和表3 所示，可以看出下垫面参数n 与NDVI 有较好的正相关关系，线性回归系数a 分别为17.90、22.14、20.97，截距b 分别为-2.19、-2.84、-3.45，表明NDVI 每增加1%时，下垫面参数n将增加20%左右，这意味着永定河流域产流性能对植被变化有极强的敏感性。回归方程的决定系数分别为0.44、0.45、0.46，对应相关系数拟合效果均达到显著性水平（P<0.01）。

表3 永定河山区Budyko 模型中下垫面参数n 与植被NDVI 拟合结果

图14 永定河山区典型流域内年均NDVI 与下垫面参数的关系

永定河山区NDVI 的增加极大促进了流域河川天然径流量的衰减作用，结合前述NDVI 与相关驱动因素的GWR 模型计算结果，NDVI 增加与降水量、平均气温作用下的自然生长有关，也与植被恢复相关的退耕还林、风沙源治理等大规模人类活动有关。根据调查，自21 世纪以来永定河流域内开展实施了多项植被恢复工程，主要有2000年开始实施的京津风沙源治理一期工程、2002年实施的退耕还林工程以及2013年实施的京津唐风沙源治理二期工程等，而2000年以前流域内无相关大型植被恢复工程，在一定程度上说明，2000年前植被的生长主要受气候条件的影响，而2000年后植被生长同时受气候条件和人类活动的共同影响作用。结合图3 永定河山区年均NDVI 变化图可知，1982—2000年期间年均NDVI 增加值为0.012，而2000 至2015年年均NDVI 增加0.043，后者远大于前者，表明2000年后大型植被恢复工程极大促进了永定河山区NDVI 的增加，同时也引起流域生态系统蒸散作用增强，导致永定河山区在同等降水条件下产流能力显著降低。有研究已指出，大型植被措施对湿润地区的水资源影响较小，而在干旱半干旱地区往往会显著减少水资源［32］，显然永定河山区植被与河川天然径流的变化情况与后者是相符的。1982—2015年永定河山区册田水库、响水堡水库、官厅水库控制流域的平均湿润指数（即降水量与潜在蒸散发的比值）分别为0.39、0.40、0.41，相应下垫面参数n 的平均值分别为3.18、2.71、2.74，表明永定河山区属于持水能力高的半干旱地区，结合文献［33］给出的中国人工植被恢复的原则，建议未来永定河流域植被恢复应以自然修复为主，对人工植被恢复工程应进行适当的控制，同时遵守植被分布规律来界定植被分区［34］。

5 结论

通过构建地理加权回归模型探究了1982—2015年永定河山区NDVI 和降水、平均气温的空间非平稳关系，并进一步剖析了NDVI 变化对流域典型断面天然径流的影响，结果如下：（1）1982—2015年永定河山区土地利用格局变化不大，但流域尺度上植被恢复明显，年均和春、夏、秋季NDVI平均趋势率为0.0013/a、0.0013/a、0.0018/a、0.0020/a，且流域内绝大多数空间像元NDVI 的年际增加趋势均显著。（2）1982—2015年，永定河山区年、季尺度下的NDVI 和降水、平均气温的回归参数、残差均值以及两者的年际变化率均具有明显空间异质性。其中回归参数及残差均值在空间像元内均为正，表明流域内降水、平均气温等气候因素和人类活动因素对NDVI均表现为增加作用，但在不同空间位置的差异亦较为显著。（3）1956—2016年永定河山区天然径流量具有显著的衰减趋势，结合水热耦合平衡方程发现流域天然径流量对植被变化具有高度敏感性，以植被变化为突出特征的下垫面变化对永定河天然径流衰减具有主导作用。（4）永定河山区属于持水能力高的半湿润半干旱地区，流域产流能力对植被覆盖状况的变化较为敏感。大规模的恢复工程的建设对河川天然径流衰减的可能影响应给予高度重视，建议未来永定河流域植被恢复应以自然修复为主，对人工植被恢复工程应进行适当控制。