水资源是人类社会经济发展和日常生活不可或缺的部分。水资源短缺已成为制约许多国家和地区经济发展的主要瓶颈^[1]。水源涵养功能是生态系统服务功能中不可或缺的一部分，对河川径流量的变化和植被涵养水源的能力评估具有重要意义^[2]。随着3S技术飞速发展，空间可视化和定量化的评估已成为现阶段水源涵养服务评价的主要发展趋势^[3]。目前，常用的评估模型有SWAT模型^[4]、TerrainLab模型^[5]和InVEST模型^[6]等。InVEST模型因数据易获取、相关参数易计算，大小尺度皆有良好的模拟效果，且输出数据可视化强，具有灵活的情景模拟功能等优点，被全世界学者广泛运用。Hame等^[7]在美国加州运用InVEST模型对当地水资源进行保护规划和决策方案筛选，Leh M.D.K等^[8]运用InVEST模型评价西非由土地利用变化引起的生态系统服务功能变化；Trusurat等^[9]评估泰国Thadee流域水源涵养能力，白杨等^[10]发现白洋淀流域西、北部和流域中部山区是水源涵养的重要区域；潘韬等^[11]发现三江源生态保护区生态系统退化严重，水源供给能力减弱。

小滦河是滦河最大的一级支流，是京津冀水源的重要供给区，属于国家重点生态功能区，对京津冀地区的防风固沙、涵养水源等方面发挥着不可替代的作用。然而近20 a来根据水文站点的记录数据发现，自1999年以来小滦河流域持续干旱，2001—2009年地表水资源量相较1956—2000年减少30.9亿m^3[12]，水资源供需矛盾日益突出，严重制约京津冀地区人民的用水安全。目前小滦河流域水源涵养功能研究较为缺乏，亟待进行深入研究。InVEST模型中的Water Yield(产水量)模块可通过分析不同土地利用类型下土壤渗透性差异，结合地形及地表粗糙程度对地表径流的影响，定量评估流域的产水和水源涵养能力。

笔者以小滦河流域为研究对象，采用InVEST模型对流域2000—2020年的产水和水源涵养功能进行全面评估，系统分析不同土地利用类型下产水和水源涵养量的时空分布格局及空间自相关特征，以期为流域水资源优化管理、保障京津冀用水安全提供科学指导。

1 研究区概况

小滦河流域(E 116°32′~118°14′，N 41°35′~42°40′)位于河北省承德市围场满族自治县境内(图 1)，流域总面积约1 646.42 km²，海拔909 m~1 871 m。小滦河是滦河上游的一级支流，北起围场县塞罕坝机械林场，南至围场县石桌子乡，东北部为以山地森林为主，向西南高原地势平坦过渡。研究区为半干旱半湿润的大陆性季风气候，冬季寒冷干燥，夏季多雨。降水时空分布不均，年际变化大，年均降水量439 mm，多集中在6—9月，蒸发量991.8 mm。土壤类型主要有灰色森林土、风沙土、沼泽土、草甸土等。围场县地处内蒙古高原与冀北山地过渡边缘，是典型的农林牧交错区。流域内土地利用方式交错复杂，主要类型为林地、草地和耕地。粮食作物以玉米、马铃薯为主，森林多为天然次生林和部分人工林，主要植被包括樟子松(Pinus sylvestris)、落叶松(Larix gmelinii)、白桦(Betula platyphylla)、山杨(Populus davidiana)等。

2 数据来源与研究方法 2.1 参数来源与预处理

Water Yield模块输入参数包括：降水量、潜在蒸散发、土地利用、土壤深度、植被可利用水含率、研究区流域、子流域、生物物理参数表以及Zhang系数。所有栅格数据分辨率为30 m，空间参考坐标系为Albers_Conical_Equal_Area，椭球体参考坐标为D_Krasovsky_1940。以2020年所需数据为例：

1) 平均降水量。从中国气象数据网(http://data.cma.cn)获得小滦河流域及周边6个气象站点(表 1)2000、2010和2020 3年日降雨数据，用反距离加权法进行空间插值，获得年均降水量的栅格图层。

2) 潜在蒸散发(potential evapotranspiration, ET₀)。FAO推荐改进的Penman-Monteith公式是计算精度最高的方法之一，因此本研究使用Penman-Monteith公式^[13]计算研究区的潜在蒸散发。

3) 土地利用类型。来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn/)，数据解译精度为30 m。

4) 土壤深度、植被可利用含水率。土壤数据由世界土壤数据库(HWSD，http://data.tpdc.ac.cn/zh-hans/)提供，比例尺为1 ∶100万。相关土壤数据经投影、提取和裁剪等步骤，获得2个栅格数据集。其中，植被可利用含水率采用周文佐^[14]提出的公式，基于土壤质地和土壤有机质数据的非线性拟合土壤AWC估算模型计算。

5) 流域及子流域划分。基于DEM数据(http://www.gscloud.cn)，结合Arc Hydro工具，通过填洼、计算流向、汇流计算等水文分析过程得到流域边界。

6) 生物物理系数表。该参数表征区域不同土地利用类型下的生物物理属性，包括土地利用类型及代码、蒸散系数、植被根系深度等(表 2)。植被根系深度取值结合相关地区的研究成果^[15]及用户指南手册进行调整，K值与不同土地利用的自身属性密切相关，计算公式如下:

$ K=\frac{\sum\limits_{m=1}^{12} K_m \mathrm{ET}_m}{\sum\limits_{m=1}^{12} \mathrm{ET}_m} 。$

(1)

式中：K为年均植被蒸散发系数K_m为m月平均植被蒸散发系数，ET_m为m月的潜在蒸散发量。

7) Zhang系数。代表区域的降水分布，为1~30之间的常数。Budyko干燥度指数理论^[16]证明Z值越高，模型受季节常数限制越小。

2.2 产水量的计算

Water Yield模块根据Budyko水量平衡原理，结合气候、地形和森林景观类型等通过降水量减去实际蒸散发量计算每个栅格的水源供给量，即产水量^[17]。产水量与降雨、蒸发两者间的平衡关系受气象要素、土壤特征和地表覆盖等的影响。模型计算时简化汇流过程，忽略地表水和地下水的交互作用，从而得出研究区的产水量。模型计算原理如下：

$ Y_{j x}=\left(1-\frac{\mathrm{AET}_{x j}}{P_x}\right) P_x ; $

(2)

$ \frac{\mathrm{AET}_{x j}}{P_x}=\frac{1+\omega_x R_{x j}}{1+\omega_x R_{x j}+\frac{1}{R_{x j}}} ; $

(3)

$ R_{x j}=\frac{k \times \mathrm{ET}_0}{P_x} ; $

(4)

$ \omega_x=Z \frac{\mathrm{AWC}_x}{P_x}+1.25 。$

(5)

式中: Y_jx为第j类土地利用栅格x中的年产水量, mm；AET_xj为第j类土地利用栅格x中的实际年蒸散量, mm；P_x为栅格单元的年降水量, mm；$\frac{\mathrm{AET}_{x j}}{P_x} $的计算采用Zhang等^[18]提出的水热耦合平衡原理计算; R_xj为第j类土地利用栅格x中干燥度指数; ω_x为土壤性质非物理参数；k为植被蒸散系数; ET₀为参考作物潜在蒸散量, mm；AWC_x为栅格 x的植被可利用水含量, mm；Z为Zhang系数。

2.3 水源涵养量的计算

InVEST模型认为植被截留的降雨最终会通过蒸散作用返回大气，对水源涵养贡献甚微，因此主要计算的是生态系统中土壤水源涵养量。基于InVEST模型计算出产水量后，需通过降水量减去蒸散发和地表径流后，结合土壤饱和导水率、流速系数、地形指数^[19]对水源涵养量进行修正，以获得最终结果。

$ \mathrm{TI}=\log (D / S \times P) ; $

(6)

$ \begin{gathered} \mathrm{WR}=\min \left(\frac{1, 249}{V}\right) \times \min (1, 0.3 \times \mathrm{TI}) \times \\ \min \left(1, \frac{K_{\mathrm{s}}}{300}\right) \times Y。\end{gathered} $

(7)

式中: WR为水源涵养量，mm；V为流速系数; Y为产水量，mm；K_s为土壤饱和导水率，mm/d；TI为地形指数；D为集水区栅格数；S为土壤深度；mm；P为坡度比。

2.4 空间自相关分析

全局自相关可以揭示地理现象中的空间集聚性，局部空间自相关则能反映全局中一个局部区域单元与相邻单元同一属性的相关度，二者分别用Moran's I和LISA指数表示。Moran's I指数是运用最广泛的全局指数，反映研究区域中邻近单元是相异、相似还是非相关，判断属性值在空间分布上是否具备集聚特征。计算公式如下：

1) 全局自相关：

$ \text { Moran's } I=\frac{n \sum\limits_{i=1}^n \sum\limits_{j=1}^n\left(x_i-\bar{x}\right)\left(x_i-\bar{x}\right)}{\sum\limits_{i=1}^n \sum\limits_{j=1}^n W_{i j} \sum\limits_{i=1}^n\left(x_i-\bar{x}\right)^2} ; $

(8)

2) 局部自相关：

$ \mathrm{LISA}_i=\frac{\left(x_i-\bar{x}\right)}{\sum\limits_i \frac{\left(x_i-\bar{x}\right)^2}{n}} \sum\limits_j W_{i j}\left(x_j-\bar{x}\right)(i \neq j) 。$

(9)

式中: n为空间单元总数; W_ij为单元i与j的连接矩阵; x_i、x_j分别为单元i与j的属性值; x为单元属性值的平均值。LISA_i无固定取值范围，LISA_i值＞0表示相似单元集聚，表现为高—高值或低—低值集聚；LISA_i值＜0表示相异单元集聚，表现为高—低值或低—高值集聚。

Moran's I的取值范围为[-1, 1]，当I＞0时，正相关; 当I＜0时，负相关; 当I=0时，随机不相关。检验Z值越高，聚类程度越高，Z值越接近0，则不存在明显集聚。本研究数据通过在ArcGIS10.3中建立500 m×500 m的渔网，经过裁剪、分区统计、连接、筛选并剔除空白值后，运用SpatialStatistics工具对小滦河流域水源涵养功能的高低聚类区进行识别。

3 结果分析 3.1 参数率定与模型验证

由于不同地区地形和自然条件有较大差异，所以运用该模型的首要任务是对模型进行参数率定和验证。

3.1.1 Zhang系数率定

Z系数是InVEST模型产水模块重要参数，主要通过经验和公式法确定。本研究选取小滦河流域沟台子水文观测站(小滦河干流控制站)1980—2017年平均天然径流量(图 2)对Z系数进行率定。将不同的Z系数依次输入Water Yield模块，对模拟产水量和水文站37年平均天然径流量0.74亿m³进行对比分析，结果表明当Z系数取值为19.8时，与实测数值最为接近。

3.1.2 模型验证

根据水量平衡原理，以水量的进入和输出为两端点, 降水量与蒸散量以及其他消耗的差即为此循环中所持有的水资源总量。模拟产水量与天然径流量具有显著相关性(P＜0.05)，从产水量模拟结果(表 3)来看，模型模拟出的产水量与小滦河流域沟台子水文站统计的30年实测径流值十分接近。

3.2 产水量和水源涵养量的时空变化

从时间上看，小滦河流域2000、2010和2020年产水量和水源涵养量输出差异显著。总产水量分别为4 250万、6 690万和1.147亿m³，流域平均产水量为147万、505和626万m³，基于栅格像元的总体产水深度为291.28、415.90、404.26 mm，整个流域在2000年产水量最低，2020年是小滦河流域产水水平最高的年份。从流域尺度上看，小滦河流域产水量的变化呈现出较强的空间异质性(图 3、图 4)。2000年时，产水量的高值聚集在研究区中下部的御道口牧场、御道口乡及老窝铺乡，低值处则出现在北部林牧场交界及南部的西龙头乡、南山嘴乡。2010年相较于2000年，产水量增长29.99%，高值区从中下部往北部御道口牧场和塞罕坝林场转移，南部的老窝铺乡、西龙头乡等产水量明显降低。

水源涵养功能表征生态系统保持水分的能力，是流域蓄水保水的体现，小滦河流域的产水量分布直接影响着水源涵养的分布(图 5和图 6)。从时间上看，2000、2010和2020年，小滦河流域水源涵养总量分别为350万、610万和730万m³，流域平均水源涵养量为279万、371万和415万m³，单位面积水源涵养深度分别为41.44、71.92和90.55 mm，呈现不断上升趋势。2000年整个流域的水源涵养能力低下，2020和2010年的水源涵养量相较2000年，分别提升52.05%和42.62%。从空间上看，2000年时的水源涵养量高值区主要集中在中下部御道口牧场和老窝铺乡，低值区集中在北部的塞罕坝林场和南部的南山嘴乡，高低值分布较分散。2020年的高值区主要聚集在流域上游靠近塞罕坝林场的人工林地和下游大面积天然林地，御道口牧场的水源涵养水平下降明显，牧场地区从高值转变为中低值。总体来看，小滦河流域的整体产水服务功能呈现稳中向好的趋势。

3.3 不同土地利用类型下产水量和水源涵养量的变化

土地利用是评价水源涵养能力的关键影响因子土地利用是评价水源涵养能力之一。不同土地利用类型下垫面粗糙度的改变直接影响植被蒸腾及地表蒸发，进而影响地表径流和下渗，最终影响产水量及水源涵养量的变化。

从不同年份来看，不同土地利用类型下产水量和水源涵养功能时空分布复杂且不均(图 7)。2000年沙地面积高达51.48 km²，为3个时期中最大，土地环境沙化明显；中覆盖度草地的产水量最大，也仅有1 100万m³，2000年流域水源涵养总量仅为350万m³。而旱地和低覆盖度草地由于地表的覆盖保护层被破坏，导致降水蒸发加剧，土壤水分被迅速消耗，蒸散发量远大于产水量，水源涵养量计算结果约为0。2010年，沙地面积从2000年的51.48 km²降至2.62 km²，沙化环境得到明显改善；全域产水量、水源涵养量增幅分别为36.47%、42.63%。灌丛产水量相比2000年增加59.15%，成为该流域对产水量贡献最大的植被覆盖类型之一。相较于2010年，2020年流域内沙地面积持续缩减，总产水量和水源涵养量稳步提升，分别增长41.67%和16.44%。2020年，林地的产水量和水源涵养量分别为2 020万和180万m³，中覆盖度草地、灌丛产水量分别占整个流域的35.92%和16.74%，与2000年相比，灌丛和中覆盖度草地水源涵养量的增幅达到56.25%和66.67%。

3.4 水源涵养量空间自相关分析

对小滦河流域水源涵养量进行全局空间自相关分析，结果见图 8，2000、2010、2020年全域Moran's I指数分别为0.64、0.75、0.81，Z值和P值(P＜0.001)均通过显著性检验(表 5)，说明小滦河流域水源涵养量空间分布呈非随机状态，集聚特征极为明显，单元高低聚类分布差异显著，进一步证明其水源涵养量空间异质性高。高—高和低—低型为正相关空间格局聚集，表明该单元及周边单元水源涵养量同步高低变化，高—低和低—高型为负相关空间格局聚集，当流域水源涵养量高时，周边则低，流域低时，周边则高。从结果来看，2000年以高—高值集聚为主，空间关联正集聚效应明显。高—高值集聚区多为零散分布，遍布流域的上中下部，低集聚区几乎不存在。2010年高—高值集聚区发生明显转移，从零散分布转变为自御道口牧场往北集聚，高—高值集聚区出现在流域北部的草地和林地附近。2010年出现低—低值集聚区域，多分散在御道口乡以南的耕地、建设用地附近，此区域实际蒸散发高于降水量，水源涵养力弱，呈现低值单元集聚分布。2020年时高—高值集聚区转变的更加集中，水源涵养功能发展稳中向好，高—高值集聚主要出现北部林场、灌丛和南部山区天然林附近，中部以草地为主的牧场低—低区集聚面积相比2010年有较大扩张，说明中部草地受外界影响较大有退化趋势，生态逐渐脆弱，导致水源涵养能力下降。

4 讨论

本文基于小滦河流域3期土地利用、土壤以及多年气象数据，定量评估该流域的总产水量和水源涵养功能的时空变化。本研究发现InVEST模型对中小面积流域产水、水源涵养量的定量评估具有更高的精确度，说明模拟结果对中小流域具有更高的适用性和参考性。20年间小滦河流域内沙地面积骤减，城镇建设用地不断扩张，不同覆盖度草地面积剧增，流域产水量和水源涵养量大幅提升。这与巩飞等^[20]、黄菁等^[21]在张家口地区的研究成果相似。单个地类中灌丛、林地和中覆盖度草地是研究区产水量和水源涵养量贡献率最高的土地利用类型。这与杨昀则等^[22]在甬江流域得出相似结果，说明这3种生态用地均具有较强的水文功能，应加强保护。然而，包玉斌等^[23]、刘娇等^[24]通过模拟实验得出黄土高原2000—2010年产水量和水源涵养量先增后降，干热河谷2000—2019年产水量则逐年下降的不同结论。可见，不同时间段流域内土地利用类型、地形、气候、植被条件等因素都对产水量变化有影响。2020年小滦河流域产水量高值区逐渐从中部牧场向北部林场转移，南部天然林地随着时间的增加，植物根系逐渐稳固扩张，植被逐渐茂盛，森林郁闭度大幅度提升，水源涵养能力明显提升，但高—高值聚集区有所减少，低—低值聚集区有所增加，说明水源涵养服务是多因素复合驱动的影响过程，小滦河流域作为典型的农牧交错区，中部大面积草、牧场受不科学放牧影响，草地生态系统明显退化，高覆盖草地多退化为旱地，地表下垫面覆盖、蒸散发量及拦截径流能力都随之发生改变，水源涵养能力也随之减弱。虽然整体上水源涵养能力稳步提升，但是从单个地类的低—低值聚集有所增加，当地应重视并采取相关措施，如围栏封育，松土改良，施有机肥，提高土壤肥力或者人工种草等。到2020年为止，小滦河流域产水量和水源涵养量达到1.147亿和0.073亿m³。与其他地区相比^[25]，小滦河流域现阶段的产水和水源涵养功能仍处于较低水平。因此，建议加强对流域内林地和草地资源的保护工作，实行弹性放牧积极恢复草地生态系统的生产力，定期进行育林护林，完善土地利用规划，有效提升小滦河流域水源涵养能力，促进流域内生态用地的良性转换。

尽管InVEST模型已被广泛运用于各种区域和情景下的生态系统服务功能评估，但是该模型的部分输入参数依赖经验取值，未来仍需进一步进行参数优化，以期得到更精密的结果。

5 结论

1) 2000—2020年小滦河流域产水量、水源涵养量均呈现稳步上升趋势，但上升速率有所减缓。空间上水源涵养量总体呈现东北高，中、南部低的分布格局，高值区从中部逐渐向上、下游转移，空间异质性显著。

2) 20年间，小滦河流域建设用地不断扩张、沙地环境明显改善。各地类中，灌丛、林地及中覆盖草地占比较大，单位面积平均产水量高且蒸散发小，在小滦河流域生态系统的水源涵养功能中发挥主要作用，而高覆盖度草地逐渐向旱地退化，水源涵养功能下降明显，需引起一定重视。

3) 小滦河流域水源涵养功能呈现出显著的空间自相关性，单元集聚特征突出。高—高值集聚区主要分布在东北部林地和灌丛周围，低—低值集聚区则分布在中部草牧场及南部耕地。在水源涵养量的高值集聚区应避免优势地类被建设用地侵占扩张，应建立严格的生态保护措施，以草牧场为主的低值区应采取相关措施减缓草地退化，以维持流域生态系统水源涵养功能的长久稳定。