张承水源涵养区土地利用演变及产水量的响应

张 薇，王凤春，3，4，贾 悦，万虹麟，5，牛天浩云，杜雨森

（1.河北水利电力学院水利工程系，河北沧州 061001；2.沧州市遥感与智慧水利技术创新中心，河北沧州 061001；3.河北省高校水利自动化与信息化应用技术研发中心，河北沧州 061001；4.河北省数据中心相变热管理技术创新中心，河北沧州 061001；5.河北省岩土工程安全与变形控制重点实验室，河北沧州 061001）

0 引 言

近年来随着水资源短缺问题的日益突出和生态系统服务功能研究的不断深入，水源供给作为生态系统服务的重要功能之一，被广泛关注［1］。生态系统的水源供给功能概念较广，其内涵随着研究的不断深入而不断完善，主要表现为生态系统拦蓄降水、调节径流、影响降雨、净化水质等，研究表明不同的土地利用/覆被（LULC）类型通过影响流域水循环过程中的蒸散、下渗和持水，进而改变流域水源供给服务功能［2］。目前，众多学者在水源供给服务功能方面做了大量研究，取得了系列成果，为相关研究提供了有益参考［3-6］。特别是InVEST 模型，相较于其他水资源评估模型以其原理简单，易于被不同学科背景的使用者理解，且需要的数据量相对较少、易于获取，近年来被广泛应用于水源供给服务功能评价相关研究中［7-9］。

密云水库上游流域张承地区作为密云水库的水源涵养区，担负着向北京市提供饮用水源的重要任务，其水源供给服务功能质量将直接关系首都的水生态安全。近年来，随着经济的发展和人类活动干扰的加剧，流域生态环境脆弱、库区水生态质量下降。其中涉及水质评价、污染源及污染评价、流域水生态改善机制等方面的研究较多，但分析产水量时空变化及其与土地利用变化的响应关系方面的研究相对较少，且不同学者研究的流域范围也不尽相同［10-12］。从土地利用结构驱动水源涵养生态服务功能的角度出发，定量评估区域产水量并分析其时空变化特征，揭示其变化与土地利用类型变化的内在关系，深入分析哪种土地利用方式有利于提高区域水源供给能力，有助于实施精准靶向生态补偿，提高补偿效率和效果。

基于此，本文借助InVEST 模型产水模块，定量评估张承水源涵养区2000-2019年的产水量并分析其时空变化特征，进一步探讨不同土地利用类型变化与产水量变化的响应关系，以期为张承水源涵养区的国土规划和京冀生态协同发展提供科学依据。

1 数据来源与研究方法

1.1 研究区概况

密云水库有白河和潮河两大入库河流。其中，白河起源于张家口沽源县，经赤城、延庆、怀柔入库；潮河源于承德丰宁县，经滦平、古北入库。本文研究区为密云水库上游流域的河北省赤城、丰宁和滦平三县，位于40°30′～41°37′N，115°30′～117°32′E，总面积约111.30 万hm2（图1）。该区域地处山地与平原的过渡地区，地势北高南低。该区域属于中温带向暖温带过渡、半干旱向半湿润过渡的中纬度大陆性季风气候，年平均气温在7 ℃左右，且随海拔的上升有所下降。全年降雨量约为490 mm，且大多数集中在7-9月份，降雨所形成的地表径流是河流的主要补给形式。人工林和天然次生林是区域内的主要植被类型，同时，灌丛有较广的分布。

图1 研究区范围Fig.1 Scope of study

1.2 研究方法

（1）土地利用转移矩阵。引进土地利用转移矩阵来定量描述土地利用类型转移情况。转移矩阵是将各地类变化的转移面积按矩阵形式排列，定量分析不同时期土地利用的变化结构和方向，直观反映一定时期内各地类间的转移方向和数量［13］，其数学模型为：

式中：U为土地利用转移矩阵；Uij表示转移前的i地类转换成转以后的j地类的面积。

（2）InVEST 模型产水模块。InVEST 模型产水模块基于Budyko 水热耦合平衡假设（1974）和年平均降水量数据，年产水量等于年降雨总量与年蒸散总量的差值［14］。该模型运行基于栅格地图，确定研究区每个栅格单元的年产水量［15］，公式如下：

式中：Yx，j为j类土地利用类型栅格x的产水量；Px为栅格x的年降水量；AETx，j为j类土地利用类型x的年实际蒸散量［16］；Rx，j是土地利用类型j栅格x的布德科干燥度指数，是潜在蒸散与降水的比值；Wx是一个非物理参数；AWCx是植被可利用的体积含水量，mm；Z为Zhang 系数，为一个经验值，其值在1～10，根据实际情况调整［17］。

式中：ET0为栅格x内的潜在蒸散量；RA为太阳大气顶层辐射；Tagv是日最高温平均值与日最低温平均值的均值；TD是最高温平均值与日最低温平均值之差；P为栅格x的年降水量［18］；AWC植被可利用的体积含水量；SAN%土壤沙粒含量百分比；SIL%土壤粉粒含量百分比；CLA%土壤黏粒含量百分比；C%土壤有机碳含量百分比。

（3）情景分析法。InVEST模型基于水量平衡原理模拟产水量，降水和实际蒸散量是影响模型模拟结果的主要因素，而蒸散量主要受植被和土地利用类型影响。为进一步探究产水量与土地利用/覆被变化的关系，这里设计了模拟情景：在2000年模型数据的基础上，只改变土地利用覆被数据，分别输入2010年和2019年土地利用数据和对应的作物生物物理参数，可研究2000-2010 和2000-2019 土地利用覆被变化对产水量的影响［19］。

式中：Rl土地利用覆被变化对研究区产水量变化的贡献率；l0基准产水量；li模拟产水量。

1.3 数据来源及预处理

选取2000年、2010年、2019年3 期数据为研究时序。基于上述方法，收集了研究区的流域界限、行政区划、土地利用、气象、土壤属性、DEM 等数据用于模型计算，同时还整理了研究区历年水资源量数据对模型结果进行验证。其中，土地利用数据主要通过遥感解译获得，遥感数据由地理空间数据云（http：//www.gscloud.cn/sources/）下载；气象数据来源于国家气象信息中心中国气象数据网（http：//data.cma.cn/），包括降雨、气温、日照等数据；土壤数据来源于世界土壤属性数据库（HWSD）；DEM数据也来源于地理空间数据云；水资源量来源于《河北省水资源公报》（http：//slt.hebei.gov.cn/）。

土地利用数据由Landsat OLI/ TM 遥感影像，利用ENVI5.3通过监督分类结合人工目视解译的方法获取，参照中国科学院资源环境数据中心土地分类系统，将土地利用类型分为6 个一级类：耕地、林地、草地、建设用地、水域、裸地，其综合精度达87%以上；降雨量空间数据由河北省13 个站点实测数据，经ArcGIS10.2空间插值、裁剪所得；小流域划分空间数据理论上应该由DEM 经ArcGIS10.2 投影变换、裁剪、水文分析流域提取所得，但这里由于研究区的不是一个完整的流域，且为了后期分析统计方便，小流域划分空间数据由行政区划矢量数据代替；潜在蒸发量空间数据是由气象站点实测数据，经公式计算后进行空间插值、裁剪得到；植物有效含水率基于土壤属性数据（沙粒、粉粒、黏粒、有机碳含量）经ArcGIS10.2 属性字段计算得到；土壤深度数据由全国土壤深度数据经ArcGIS10.2 裁剪得到（见图2）。

图2 基础处理数据（2010年6月）Fig.2 Processed basic data（2010，Jun）

2 结果与讨论

2.1 土地利用结构时空变化分析

由遥感解译得到2000年、2010年和2019年3 期的土地利用分布图（见图3）。分析可知：林地的分布最为广泛，其次是草地、耕地和建设用地，耕地分布在河流沿岸地势平坦地区，草地则沿河流走势分布在农业生产用地的外围，林地则主要位于离河较远的丘陵、山地区域。

图3 2000-2019年研究区土地利用空间分布Fig.3 Spatial distribution of Land use of the study area in 2000-2019

经统计得2000-2019年各地类面积及变化情况，如表1 所示。张承水源涵养区土地利用结构中林地和草地是主要的用地类型，2019年两种用地类型面积分别为793 461.70 和138 779.25 hm2，各占总面积的71.3%和12.5%，作为密云水库上游流域涵养区，林地、草地有利于水量的储存和净化。2000-2019年的土地利用结构转型中，耕地、裸地和林地总体减少；草地、建设用地和水域总体增加，其中建设用地面积增长较快，19年间增加了4 268.32 hm2，年均涨幅2.1%，表明随着城镇化水平的不断提高，研究区建设用地需求量增加。

表1 2000-2019年各地类面积及变化 hm2Tab.1 2000-2019 Area and change of each category

进一步分析各土地利用类型相互转化的方向和数量，在图3 的3 期土地利用分类图的基础上，利用ArcGIS 的空间分析功能对不同时期的土地利用图进行空间分析，获得研究区2000-2010年和2010-2019年的土地利用功能类型的转移矩阵（表2、3）。结果总体表现为建设用地面、草地面积增加，林地、耕地面积减少。2000-2010年，耕地、林地、建设用地面积均有所增加，增加的面积大多都是由水域、草地转化而来，其中草地分别向耕地、林地、建设用地转化了1 115.59、5 773.18 和92.08 hm2，裸地有86.3%开垦成了耕地。2010-2019年，林草地面积之和增加，耕地和裸地面积均有减少，且减少的耕地主要转化为林草地、水域和建设用地，分别为34 655.8、384.68 和4 231.44 hm2；建设用地持续增加，林草地、耕地、裸地、水域均不同程度向其转化，转化率分别为：0.8%、2.4%、2.2%和1.5%；水域面积有所增加，增加的面积由林地、耕地、建设用地转化而来。研究区作为密云水库的水源涵养区，长期以来始终以保护生态环境作为地区发展的首要条件，充分利用地区林地、草地面积较多的优势，优化区域水源涵养生态服务功能，按照《京津冀协同规划纲要》的相关规定，多措并举推进生态清洁小流域建设。特别是随着《河北省张承地区生态保护和修复实施方案》等一系列政策的落地实施，对研究区土地利用结构变化起到一定的促进作用。

表2 2000-2010土地利用变化矩阵 hm2Tab.2 2000-2010 Land use change matrix

表3 2010-2019土地利用变化矩阵 hm2Tab.3 2010-2019 Land use change matrix

2.2 产水量时空分布及变化特征

InVEST 模型在估算年产水量时除了数据预处理中所列的数据外，还需要输入Zhang系数和生物物理参数表。Zhang系数是InVEST 模型需要输入的重要参数，取值范围1～10，可以由研究区年降雨发生次数进行初步估算，再结合模型检校数据进行反复调试。生物物理参数表至少包含各土地利用/覆被类型的植被蒸散系数Kc和最大根系深度两项参数，其中植被蒸散系数Kc是作物实际蒸散量与参考蒸散量的比值，经多次试验，发现研究区模型产水量的预测结果对该指标敏感性较高。而以往的研究中Kc的确定多是直接采用联合国粮农组织作物参考值［20］，由于植被蒸散系数受植被的生长状况、气候及水分条件等诸多因素的影响，最好结合当地实际情况进行校正。目前基于遥感估算作物蒸散系数Kc的研究相对较多［21，22］，在韩文霆、李霞［23，24］等人研究成果基础上，结合杨洁、吴瑞［25，26］等人的研究确定本文3 期生物物理参数如表4 所示。以河北省水资源公报数据进行模型检校，经过反复调试当3 期Zhang 系数分别为1.1、3.5、3.0 时，模型模拟产水量与实际产水量最接近，相对误差均控在20%以内（表5）。

电装公司是汽车零部件及系统的全球知名供应商，于1949年在日本成立，目前在全球30多个国家和地区,设有180多家关联公司，拥有员工约14万名。电装公司在环保、安全、舒适和便利等广泛领域中，推进技术开发，并提供众多产品。2003年，电装（中国）投资有限公司正式成立，并在长春、天津、上海、广州、武汉、济南和重庆设有分公司。

表4 InVEST模型的各土地利用类型生物物理参数Tab.4 Biophysical parameters of land use description in InVEST model

表5 InVEST年产水模型结果校验Tab.5 Verification of result about InVEST water yield model

2000-2019年，研究区产水深度在7.2～65.1 mm。2000、2010 和2019年的总产水量分别为2.61、3.46 和2.71 亿m3，说明研究区产水量先增加后减少（图4）。不同时期的产水量表现出空间分异格局，整体变化趋势为产水量由西高东低逐渐变为东高西低。2000年产水主要集中在赤城县和丰宁县西北部；2010年赤城中东部和丰宁西部地区产水量减少，丰宁东部和滦平县地区产水量增加；2019年研究区产水量呈东高西低的趋势更加明显。这种格局与研究区年均降水量和土地利用类型分布有直接关系，降水量多且植被蒸散量低的区域产水能力强，反之产水能力较弱。

图4 2000-2019产水深度空间分布Fig.4 Spatial distribution of the water yield depth in 2000-2019

由图5 可以看出，2000-2019年，研究区产水量先增加后减少，总产水量略有增加。2000-2010年增加区域主要集中在滦平县和丰宁县东北部，以及赤城县的西北部，增加的面积占总面积的76.5%；2010-2019年产水量总体减少，其中赤城县、丰宁县中部和滦平减少最多。2000-2019年总产水量略有增加，其中以巴克什营镇、涝洼乡和两间房乡产水深度增加最多，分别增加了28.5、31.5 和29.4 mm；产水量减少的区域主要集中在赤城县和丰宁县西部地区。2000-2019年研究区的产水量变化与不同时期的降雨和土地利用结构有直接的关系。

图5 2000-2019产水深度变化的空间分布Fig.5 Spatial distribution of the water yield depth change in 2000-2019

为了定量描述3 个区县2000-2019年产水量的变化情况，分别对赤城、丰宁、滦平3个县域三期的总产水情况进行了统计分析（表6）。2000-2019年3 个县域内产水整体呈现先增加后减少趋势，跟整个研究区三期产水趋势吻合。其中，赤城县2000-2019年产水总量总体减少0.177 亿m3，丰宁县和滦平县2000-2019年产水总量总体分别增加了0.005和0.277 亿m3。

表6 2000-2019各县模拟产水量及变化 亿m3Tab.6 Water yield and change in each county of 2000-2019

2.3 不同土地利用类型的产水深度和产水量

InVEST 模型的产水模块是一种基于水量平衡法的估算方法，某栅格单元的降水量减去实际蒸发量即为该栅格的产水量，即总体的产水能力与降雨量成正比，与蒸散发量成反比，与坡度高程等地形因子无关。由于不同土地利用类型的蒸散发能力、土壤含水量、凋落物持水能力及冠层截留量存在差异，导致同一时期不同土地利用类型的产水深度不同。2000-2019年研究区的各地类平均产水深度如图6所示，分析可知：不同土地利用类型平均产水深度从大到小依次为：建设用地、裸地、水域、耕地、林地、草地，各地类单位面积的平均产水深度分别为：30.786、27.638、27.564、27.304、26.572 和26.547 mm。其中建设用地植被覆盖较少，相对蒸散发量也就较少，再加之不透水面的增加，使得降水渗入地下较少，导致该地类产水量最高；草地和林地地表调落物拦截地表径流，延迟降雨汇流时间，增加土壤入渗量，同时植被冠层蒸散量相对较强，导致林地和草地产水量相对较少。

各地类的产水能力和面积占比是影响该地类总产水量的主要因素，通过图6 分析除建设用地外，其他5 种地类平均产水能力相差不是很大，这样研究区各地类总产水量主要取决于各地类面积占比。研究区主要的用地类型是林地，2019年，林地占整个研究区的71.3%，产水量占产水总量的72.0%。

图6 不同土地覆被类型平均产水深度的时间变化Fig.6 Temporal variation of average water yield depth of different land cover types

2.4 产水量与土地利用相关性

由InVEST 模型产水模块的估算原理可知，某栅格产水量主要受降雨和蒸散发两个因素的影响，由于不同土地利用类型的蒸散发能力不同，且不同地类的凋落物持水能力及冠层截留量也存在一定的差异，导致不同土地利用类型将直接导致产水量的变化。本文采用情景模拟法，在2000年模型数据的基础上，分别输入2010年和2019年土地利用数据和对应的作物生物物理参数，进行情景模拟得2010年和2019年研究区总产水量，并分析两期各地类的平均产水深度，以探究研究区产水量与土地利用的相关性。结果显示2010年和2019年研究区总产水量相较2000年实际产水总量2.61 亿m3分别减少1.42 亿m3和增加0.15 亿m3，地表覆被变化对产水量变化的贡献率分别为-54.4%和5.7%。对比分析对应的生物物理参数，2019年各地类作物系数与2000年的相同，而2010年各地类作物系数较2000年略有增加，作物系数代表各类地物的蒸散发能力，进一步说明该产水模型对各地类蒸散系数较敏感，且产水量与作物系数成反比。通过对比分析2019年各地类模拟产水深度，如图7 所示，发现与2000年相比，2019年各地类平均产水深度均有所增加，其中水域平均产水深度增幅最大由2000年的22.5 mm增加到2019年的27.8 mm，林地平均产水深度增幅最小由2000年的23.5 mm 增加到2019年的24.7 mm。通过情景模拟分析说明土地利用类型之间的相互转移可能会导致产水量的增加或减少，例如在其他条件不变的情况下，建设用地增加会增加产水量，而林地增加会减少产水量。

图7 模拟情景下各地类平均产水深度变化Fig.7 Change of average water yield depth in different regions in simulated situations

图9 2000-2019产水量变化分布Fig.9 Spatial distribution of water yield and change in 2000-2019

为了进一步探讨区域土地利用与产水量之间的时空分布关系，将2019年模拟产水量与2000年实际产水量变化以及两期的土地利用变化做时空对比分析，如图8、9 所示。通过对比分析发现：赤城县的西南部、丰宁县的中部以及滦平县的南部产水量明显增多，该地区主要土地利用类型变化为耕地或草地转为建设用地，林地或草地转为耕地，转换后的土地利用类型相较耕地产水能力强，因此该地区整体产水量增加；产水量减少的4 个乡镇：赤城县的白草镇、丰宁县的胡麻营乡、滦平县的邓厂满族乡和付家店满族，主要土地利用类型变化为耕地、建设用地、林地转为草地，转化后的土地利用类型较转化前的地类产水能力较弱，导致该区域产水量减少。通过土地利用与产水量时空变化关系的分析，再次证明了在6种土地利用类型中，建设用地产水能力最强，其次是耕地，林地和草地产水能力最弱，蓄水能力较强。

图8 2000-2019土地利用变化分布Fig.8 Spatial distribution of the land use change in 2000-2019

3 结 语

通过遥感解译张承水源涵养区2000年、2010年和2019年三期土地利用数据，用InVEST 模型产水模块从空间上量化评估了2000-2019年研究区的产水量，同时模拟不同土地利用变化对研究区产水功能的影响，最终借助GIS 空间分析功能分析了研究区产水功能的时空分布特征及其变化，并进一步探讨了产水功能变化与土地利用类型变化的内在关系，通过研究分析得出以下结论。

（1）2000-2019年，研究区土地利用整体呈现耕地面积减少，林地、草地、建设用地不断增加的特点。将其划分为2000-2010年和2010-2019年两个阶段，前期水域和草地减少，转化为耕地、林地和建设用地；后期耕地和林地减少，转化为草地、水域和建设用地。

（2）2000年、2010年和2019年研究区的总产水量分别为2.61、3.46 和2.71 亿m3。不同时期的产水量表现出空间异质性，整体变化趋势为由西高东低逐渐变为东高西低。2000-2010年产水量增加区域占总面积的76.5%，主要集中在滦平县和丰宁县东北部，以及赤城县的西北部；2010-2019年产水量总体减少，其中赤城县、丰宁县中部和滦平减少最多。

（3）研究区不同土地利用类型平均产水深度从大到小依次为：建设用地、裸地、水域、耕地、林地、草地，通过情景模拟分析发现产水模型对各地类蒸散系数较敏感，在蒸散系数不变的情况下，土地利用由产水能力低的地类转为产水能力高的地类会造成区域产水量的增加，反之减少。

（4）采用了土地利用的6个一级分类进行了产水量分析，但由于水源涵养功能的涵义较广，未来有必要将土地利用分类进一步细化，进一步考虑区域土壤的渗透性、地表径流系数等因素，以便更好地分析土地利用和人类活动对张承地区水源涵养功能的影响。

（5）未来应进一步校验模型所需的各类参数，精确研究区域的本地化参数，提高模型精度。