承德市武烈河流域水源涵养功能的时空特征

王 冶，薛忠财，王 琳，常佳宁，姜白杨

（1. 河北民族师范学院，河北 承德 067000；2. 河南省地质环境监测院，河南 郑州 450016）

水源涵养是生态系统在一定时空范围和条件下将水分保持在系统内的过程和能力，是陆地生态系统中的重要生态服务功能之一，包含着大气、水分、植被和土壤等自然过程，其变化将直接影响区域气候水文、植被和土壤等状况[1-3]，在社会经济的可持续发展中发挥着重要作用[4]，常被视为区域生态系统状况的重要指示器。因此，科学地评价水源涵养功能一直以来是生态学和水文学等相关学科的研究热点。随着学科交叉研究的不断发展，水源涵养功能评价已从单一的径流调节功能评价提升为综合考虑各种表现形式的全面水源涵养服务评价，人们逐步认识到水源涵养服务在区域生态安全和人类福祉等方面的重要性[3,5]，并开展了大量研究。早期水源涵养评估方法多是基于大量实测数据进行估算的，主要包括水量平衡法[5]、土壤蓄水法[6]、分布式水文法[7]、降水存储量法[8]等，但上述方法所需数据量过于庞大，难以获取较高精度数据，从而进行大区域估算时会产生较大误差。近年来，随着数学方法的更新、多学科理论的交叉、计算机科学的不断进步，诸如水文评估模型Swat[9]、水文生态模型InVEST、生态系统评价模型GUMBO[10]等已被应用于水源涵养量分布及动态特征研究。其中，由斯坦福大学开发的InVEST 模型由于其数据易获取、评价结果实现空间化、动态量化、灵活性强等优点，在国内外得到广泛应用。Leh 等[11]成功利用InVEST模型对西非国家加纳和科特迪瓦的产水量进行了评估；我国学者利用InVEST 模型对黑河流域[12]、三江源地区[13]、石羊河流域[14]、黄土高原[15]等地区的水源涵养时空分布特征进行了评价，并取得了较好的评估效果。

承德市位于河北省东北部，是连接京津冀辽蒙的区域性中心城市，为京津冀地区“阻沙源、涵水源”发挥着重要的生态作用。承德市境内有滦河、潮河、辽河、大凌河、武烈河等多条水系，年均自产水量37.6 × 108m3，年均向京津等下游地区供水25.2 × 108m3。其中，武烈河流域是承德地区主要的供水水源地，也是京津冀水源涵养功能区的重要组成部分，具有重要的生态价值与研究意义。但是，对该区域长时间序列的水源涵养功能时空分布特征研究仍未见报道。因此，本研究以武烈河流域为研究区域，基于InVEST 模型产水模块，通过参数率定与模型校准，精确地估算1990–2018 年武烈河流域的水源涵养量时空分布特征，分析气候变化与土地利用变化对流域水源涵养功能的影响，并采用情景模拟的方法进一步分析导致武烈河流域水源涵养功能发生变化的主导因素，以期为京津冀水源涵养功能区的建设与合理规划提供科学支撑。

1 研究区概况与研究方法

1.1 研究区概况

武烈河流域位于燕山山脉以东，属于暖温带和寒温带过渡地带，界于117°45′–118°26′ E，40°54′–41°42′ N，干流全长110 km，流域面积2 589 km2，主要包括武烈河、玉带河、鹦鹉河、兴隆河4 条水系(图1)。气候为典型的大陆燕山山地气候，年均气温8.9 ℃，年均降水量为537.2 mm，年均径流量2.60 ×108m3，多年平均流量6.93 m3·s-1，降水时空分布与径流的年内分配十分不均匀，降水期与汛期主要集中在6 月 – 9 月，占全年降水量的77.9%和多年平均径流量的73.7%。流域内主要土地覆被类型为林地，占总面积的50%，其次为草地和耕地，主要土壤类型为棕壤、褐土和草甸土。流域水源主要补给方式为大气降水、地表径流输入，主要涵养形式为地表径流、土壤贮水、枯落物持水等。

图1 研究区示意图Figure 1 Location of Wulie River Basin

1.2 数据获取与研究方法

1.2.1 评估方法

本研究基于InVEST 模型产水模块进行分析与评价。依据水量平衡原理，通过各栅格单元的降水量、蒸发量、土壤深度、土壤质地和植被根系深度等参数估算产水量，再利用地形指数、土壤基础数据计算得到各单元水源涵养量及深度，该模型综合考虑不同土地利用类型的土壤渗透性、地形差异、地表粗糙程度等对地表径流的影响，能够较好地表达流域水源涵养量的空间分布状况及其影响水源涵养量的主要因素[16]。

1.2.2 数据来源

本研究收集了承德市武烈河流域气象、土地利用、土壤属性、流域矢量边界、数字高程模型(digital elevation model, DEM)等数据用于模型计算，同时还整理了武烈河水文站径流监测数据对模型结果进行验证。其中DEM 数据来源于地理空间数据云平台(http://www.gscloud.cn)，7 期土地利用数据(分别为1980、1990、2000、2005、2010、2015、2018 年)来源于资源环境数据云平台(http://www.resdc.cn)，土壤数据来源于世界土壤属性数据库(Harmonized World Soil Database)，气象数据(包括降水、气温、气压、风速等)来源于国家气象信息中心 (http://data.cma.cn)与国家青藏高原科学数据中心(http://data.tpdc.ac.cn)。

1.2.3 数据处理

降水数据通过使用ArcGIS 10.6 对在国家气象信息中心与国家青藏高原科学数据中心下载的全国月均降水量栅格数据进行叠加、裁剪获取。土地利用与土壤质地数据经ArcGIS 10.6 对在资源环境数据云平台和世界土壤属性数据库下载的数据进行投影变换、裁剪获取。集水区利用ArcGIS 10.6 水文分析模块的工具，根据流域与子流域的边界将研究区DEM 数据分割成小流域，本研究区共产生4 个子流域(图2)。

图2 1990 年武烈河流域基础数据Figure 2 Basic data processing in Wulie River Basin in 1990

2 结果与分析

2.1 参数率定及模型验证

2.1.1 潜在蒸散发量计算与校准

潜在蒸散发量(ET0)是InVEST 模型计算产水量的重要基础数据，表示充分供水条件下的区域蒸散发能力。前人研究表明，根据FAO-56 Penman-Monteith (PM)公式(9)计算的ET0精度较高，但覆盖全流域的数据资料难以获得[18]。因此，本研究基于研究区的现有承德气象站数据资料，采用线性回归的方法，通过PM 公式计算结果对Hargreaves (Har)公式(10)进行校准，计算结果如公式(11)所示。

图3 承德气象站潜在蒸散发量校准前后对比Figure 3 Comparison of potential evapotranspiration before and after calibration at Chengde Weather Station

2.1.2 Zhang 系数率定

Zhang 系数是InVEST 模型产水模块的必要输入参数，因此该值的精确程度对水源涵养模型的计算结果存在较大的影响，取值范围在1～10[17]。传统确定Zhang 系数的方法主要为经验法[15]和公式法[19-20]，与实际情况存在一定偏差。根据水量平衡原理，当人为因素干扰相对较小时，采用InVEST模型产水模块计算的产水量与自然径流量近似相等，因此本研究选取1970-1989 年平均自然径流量对Zhang 系数进行率定，该时期受人为因素干扰较小，且包含流域丰水年与枯水年的径流量和降水量，能够较为客观地反映武烈河流域的自然径流特征。将不同的Zhang 系数依次输入InVEST 模型产水模块，对模拟产水量与1970-1989 年平均自然径流量进行对比分析(图4)，结果表明模拟产水量随Zhang 系数增加而减小，当Zhang 系数取值为3.7 时，模拟效果最佳。

图4 不同Zhang 系数对应的模拟产水量与年均径流量Figure 4 Comparison of simulated water yield and natural runoff corresponding to different Zhang values

2.1.3 其他参数率定

水源涵养计算中涉及到其他参数包括土壤饱和导水率、流速系数、地形指数、植被可利用含水量，均可通过相关数据计算获得。其中土壤饱和导水率的计算是基于世界土壤属性数据库(HWSD)中研究区内土壤粘粒、粉粒和粗砂百分比含量，通过Neuro Theta 软件计算得出的；地形指数的计算是基于土壤深度、百分比坡度和汇水面积由公式(2)计算得到的，其中土壤深度数据来源于世界土壤属性数据，百分比坡度与汇水面积根据DEM 分别基于ArcGIS 10.6 空间分析模块与水文分析模块获得；流速系数根据各期土地利用情况通过模型参数表获取[21]；植被可利用含水量是基于世界土壤属性数据库中研究区内土壤粘粒、粉粒和粗砂百分比含量由经验公式(8)计算得到的。相关参数空间分布如图5 所示。

图5 水源涵养模型所需的其他空间参数Figure 5 Other spatial parameters of the water conservation model

2.1.4 模型验证

根据水量平衡原理，流域总产水量与自然径流量近似相等，因此模型计算结果通过武烈河水文站径流量监测数据进行验证。人为水资源开采会对自然径流量的监测产生干扰，因此本研究采用人为因素干扰较小的1970-1979 年与1980-1989年分别作为模型的校正期与验证期，结果表明(图6)，模拟产水量与自然径流量具有显著相关性(P＜ 0.05)，模型校正期和验证期的决定系数(R2)均大于0.8，纳什效率系数(Nash-Sutcliffe efficiency coefficient)均大于0.7 (表1)，说明该模型模拟效果较好，具有可行性。

表1 模拟产水量误差分析Table 1 Error analysis of simulated runoff

图6 1970−1989 年武烈河流域自然径流量与模拟产水量的相关性Figure 6 Correlation between natural runoff and simulated water yield in Wulie River Basin from 1970 to 1989

2.2 武烈河流域水源涵养功能时空分布特征

利用1990-2018 年的气象数据和7 期土地利用数据对承德武烈河流域29 年来水源涵养量进行估算，结果表明，武烈河流域多年平均水源涵养深度为83.85 mm，平均水源涵养总量为2.18 × 108m3。从年际变化趋势来看(图7)，1990–2018 年武烈河流域水源涵养深度呈现先减小后增大的变化趋势，总体呈下降趋势但趋势不明显。为准确量化武烈河流域多年平均水源涵养深度变化，将数据按照每10 年一段进行分析。1990–1999 年水源涵养深度呈减小趋势；2000–2009 年和2010–2018 年水源涵养深度呈增大趋势但趋势不明显。水源涵养深度最大值出现在2016 年，为187.4 mm，最小值出现在2009 年，为23.27 mm。

图7 武烈河流域多年平均水源涵养深度Figure 7 Annual average water conservation depth in Wulie River Basin

各水系年均水源涵养深度与流域总体变化趋势相同，呈现先减小后增大的变化趋势(表2)。玉带河水系多年平均水源涵养深度最大；武烈河水系多年平均水源涵养总量最大；兴隆河水系多年平均水源涵养深度与水源涵养总量均最小。(表2)。

表2 武烈河流域不同时期各水系水源涵养深度和水源涵养量Table 2 Water conservation and water conservation depth of each water system in different periods of Wulie River Basin

1990-1999 年全流域水源涵养深度减小，东北山地区域变化量最大，平均变化量为-115.99 mm；2000-2009 年河道及河岸地区水源涵养深度减小，山地区域水源涵养深度增大，该时期流域整体水源涵养深度增大，平均变化量为33.97 mm；2010-2018年流域整体水源涵养深度增大，东北山地区域变化量最大，平均变化量为17.14 mm；1990-2018 年武烈河流域水源涵养深度变化量呈现自北向南递增的分布规律，整体水源涵养深度减小，平均变化量为-30.11 mm (图8)。

图8 武烈河流域各时期水源涵养深度变化空间分布Figure 8 Spatial distribution of water conservation depth variation in the Wulie River Basin in different periods

2.3 武烈河流域水源涵养功能影响因素的分析

2.3.1 气候因素

在各类气候因素中，降水量、潜在蒸散发量与气温是影响流域水源涵养功能的主要因素。从年际变化趋势来看(图9)，1990-2018 年流域整体降水量呈减小趋势，潜在蒸散发量与气温均呈增大趋势，但趋势均不显著(P＞ 0.05)。在每10 年变化中，1990-1999 年流域降水量呈减小趋势，潜在蒸散发量显著增大(P＜ 0.05)，气温呈增大趋势；2000-2009年流域降水量呈增大趋势，潜在蒸散发量与气温均呈减小趋势；2010-2018 年流域降水量、潜在蒸散发量与气温均呈增大趋势。相关性分析表明(图9)，武烈河流域降水量与水源涵养深度呈显著正相关关系(P＜ 0.05)；潜在蒸散发量与水源涵养深度呈显著负相关关系(P＜ 0.05)；气温与水源涵养深度呈负相关关系，但不显著(P＞ 0.05)。

图9 武烈河流域年均降水量、年均潜在蒸散发量、年均气温及其与年均水源涵养深度的相关性Figure 9 Variation of annual average rainfall, annual potential evapotranspiration, and annual average temperature and correlation analysis with annual average water conservation depth in Wulie River Basin

2.3.2 土地利用类型

土地利用是水源涵养计算的关键影响因子之一。由于研究区内水域、未利用地面积占比较小，不予以统计，主要统计的土地利用类型为耕地、林地、草地及建筑用地。根据武烈河流域6 期(1990、2000、2005、2010、2015、2018 年)的土地利用类型面积统计表明(图10)，多年平均面积占比最大的为林地，其次分别为草地、耕地、建筑用地。1990-2018 年耕地面积占比、林地面积占比、草地面积占比均显著减小(P＜ 0.05)；建筑用地面积占比显著增大(P＜ 0.05)；1990-2005 年各土地利用类型面积变化不明显，2005-2018 年各土地利用类型面积均显著变化。

图10 武烈河流域各土地利用类型面积占比变化Figure 10 Variation in the proportion of land use types in Wulie River Basin

由于土壤性质、植被类型与覆盖程度的不同，各土地利用类型的水源涵养功能也存在明显的差异。从1990-2018 年武烈河流域各类土地利用引起的水源涵养量变化(表3)来看，流域多年平均水源涵养贡献率与水源涵养深度最大的土地利用类型均为林地，其次为草地、耕地与建筑用地。从用地类型的角度讲，流域内耕地、林地、草地面积减小，水源涵养量也相应减少，建筑用地面积增加，水源涵养量对应增加，其中面积下降幅度最大的耕地水源涵养量下降1.25 × 106m3，面积增加幅度最大的建筑用地水源涵养量对应增加2.86 × 106m3。从水源涵养量的增减幅度来讲，建筑用地的水源涵养增量并不能补偿耕地、林地、草地的水源涵养量的损失。从单位面积变化对水源涵养量影响的绝对量来看，其强度由大到小依次为林地、草地、耕地、建筑用地，对水源涵养影响强度最大的林地用地类型面积每减少1 km2，水源涵养量损失6.83 × 105m3，对水源涵养影响强度最小的建筑用地类型面积每增加1 km2，水源涵养量相应只增加0.87 × 105m3。对比各用地类型面积增减幅度及水源涵养量变化关系发现，虽然面积与水源涵养量呈正相关关系，但面积增减仅意味着水源涵养量的增减，其增减幅度并不能表征水源涵养量的增减幅度。

表3 1990−2018 年武烈河流域各类土地利用引起的水源涵养量变化Table 3 Variation in water conservation caused by land uses in the Wulie River Basin from 1990 to 2018

3 讨论

在全球变暖的大背景下，气候变化是影响武烈河流域水源涵养功能的主导因素，而土地利用变化也是武烈河流域水源涵养功能改变的不可忽视的重要影响因素。

3.1 武烈河流域水源涵养功能变化的主导因素

武烈河流域降水量与水源涵养深度呈显著正相关关系；潜在蒸散发量与水源涵养深度呈显著负相关关系；气温与水源涵养深度呈负相关关系，但不显著，各气候因素的变化趋势与流域水源涵养量的变化趋势基本相同，因此气候因素是武烈河流域水源涵养功能的主要影响因素。该结果与黑河流域[12]、白龙江流域[13]、京津冀地区[22]等研究区的结论一致。土地利用类型不同会导致土壤物理化学性质、地表粗糙程度、植被类型、植被覆盖度、根系深度等不同，从而影响地表蒸散发量与土壤储水、阻水的能力[23]，因此不同土地利用类型的水源涵养功能存在明显差异。本研究表明，对武烈河流域水源涵养影响强度最大的土地利用类型为林地，其次为草地、耕地与建筑用地。随着城市化与工业化的加剧，林地与草地面积占比均呈减小趋势，建筑用地逐年增加，进而影响流域整体的水源涵养功能。

气候与土地利用的变化是导致武烈河流域水源涵养量发生变化的主要因素，但两者的影响程度存在差异，为进一步分析其最主要因素，本研究采用情景模拟的方法，对两者的影响程度进行研究。在1990-1999 年、2000-2009 年、2010-2018 年3 个时期中，假设武烈河流域的土地利用方式不发生变化，气候要素按实际情况输入计算，从而得到土地利用方式不变时的模拟水源涵养深度WRL，将该值与实际模拟量WRO做差，得到土地利用变化影响量ΔWRL；同理，假设武烈河流域的气候不发生变化，土地利用要素按实际情况输入计算，从而得到气候不变时的模拟水源涵养深度WRC，将该值与实际模拟量WRO做差，得到气候变化影响量ΔWRC；ΔWRO表示各时期实际水源涵养深度的变化量。

情景模拟计算结果(表4)显示，1990-1999 年，实际水源涵养深度减小115.99 mm，其中气候变化与土地利用变化分别导致水源涵养深度减小114.78、1.21 mm，两者对水源涵养量均起到抑制作用。该时期降水量大幅度减小，潜在蒸散发量与气温增大，因此导致水源涵养量大幅度减小(图9)；该时期土地利用变化不明显，气候变化影响是土地利用变化影响的94.86 倍，因此气候变化对水源涵养量的减小起主导作用(图10)。

表4 气候与土地利用变化对武烈河流域水源涵养深度的影响Table 4 Impact of climate and land use changes on water conservation depth in the Wulie River Basin

2000-2009 年，实际水源涵养深度增加33.96 mm，其中气候变化导致水源涵养深度增大36.78 mm，对水源涵养量起到正向促进作用，而土地利用变化导致水源涵养深度减小2.82 mm，对水源涵养功能起到反向抑制作用。该时期降水量增加，潜在蒸散发量与气温减小，有利于水源含量的提升(图9)；林地、草地、耕地面积占比减小，建筑用地面积占比增加不利于水源涵养量的提升(图10)。气候变化对水源涵养深度的影响是土地利用变化对水源涵养深度影响的13.04 倍，因此该时期气候变化是水源涵养量增加的主导因素。

2010-2018 年，实际水源涵养深度增加17.14 mm，其中气候变化导致水源涵养深度增大21.46 mm，对水源涵养量起到正向促进作用，而土地利用变化导致水源涵养深度减小4.32 mm，对水源涵养功能起到反向抑制作用。该时期降水量增加，潜在蒸散发量与气温减小，有利于水源含量的提升(图9)；林地、草地、耕地面积占比持续减小，建筑用地面积占比大幅增加，抑制了水源涵养量的提升(图10)。气候变化影响是土地利用变化影响的4.97 倍，因此该时期气候变化是水源涵养量增加的主导因素。

综合表明，1990-2018 年气候变化是武烈河流域水源涵养量变化的主导因素，但是随着城市化的加剧，建筑用地面积的大幅增加和林地、草地、耕地面积的减少，导致土地利用变化逐渐成为了抑制水源涵养量提升的主要抑制因素。

3.2 模型精度和不确定性

InVEST 模型目前已被大量用于不同地区生态系统服务功能的评价与决策，该模型的应用较成熟，但由于模型自身的设定和简化算法等原因，仍存在一定的局限性。首先InVEST 模型在产水量计算中并未考虑土壤深层含水量[24]、人为地下水开采以及地貌等因素的影响，并且由于其算法的简化，导致模拟结果是以年份为单位的产水量平均值，难以在较小的时间尺度上精准计算产水量；其次，InVEST 模型计算所需数据的精度在很大程度上影响了模型计算结果，本研究基础数据来自地理空间数据云平台、资源环境数据云平台、世界土壤属性数据库(HWSD)、国家青藏高原科学数据中心等，然而部分基础图层分辨率较低，以及土壤成分、地表环境随时间的变化导致参数的精度很难保证；最后，InVEST 模型中Zhang 系数的选取在很大程度上会影响模型计算结果，本研究通过反复校验取Zhang系数为3.7，通过模型验证，与实际情况较为符合，但该值精度进一步的提升还需结合大量的实测数据进行修正和验证。

4 结论

本研究基于InVEST 模型产水模块，分析了承德武烈河流域水源涵养功能的时空分布特征，主要得到以下结论：

1)当Zhang 系数为3.7 时，模拟效果最佳；1990-2018 年，武烈河流域水源涵养量呈先减小后增加，总体呈微弱减小的变化趋势，多年平均水源涵养深度为83.85 mm，平均水源涵养总量为2.18 ×108m3。

2)武烈河流域降水量与水源涵养深度呈显著正相关关系；潜在蒸散发量与水源涵养深度呈显著负相关关系；气温与水源涵养深度呈负相关关系，但不显著。

3)流域各土地利用类型多年平均水源涵养贡献率由大到小依次为林地(66.09%)、草地(24.23%)、耕地(8.77%)、建筑用地(0.9%)；单位面积变化对水源涵养量影响强度由大到小依次为林地(1.21 × 105m3·km-2)、草地(6.83 × 105m3·km-2)、耕地(1.84 × 105m3·km-2)、建筑用地(0.87 × 105m3·km-2)。

4)根据情景模拟表明，1990-1999 年气候变化对武烈河流域水源涵养量提升起到反向抑制作用；2000-2018 年气候变化对水源涵养量提升起到正向促进作用；土地利用变化始终对水源涵养量提升起到反向抑制作用。气候变化是武烈河流域水源涵养量变化的主导因素，但是随着城市化的加剧，建筑用地面积的大幅增加和林地、草地、耕地面积的减少，导致土地利用变化逐渐成为了抑制水源涵养量提升的主要抑制因素。