王 冶，薛忠财*，武智勇，靳 松，姜白杨

(1.河北民族师范学院资源与环境科学学院，河北承德 067000；2.河北省山区地质环境重点实验室，河北承德 067000；3.承德市古代建筑保护研究所，河北承德 067000)

产水服务与氮磷净化能力对流域的水质及水量的变化有直接的影响，是流域生态系统研究中非常重要的两项服务[1]。其中，产水服务常被作为区域生态系统状况的重要指示器，而氮磷净化能力则是水质净化功能的重要体现[2]。产水服务研究大多是基于大量实测数据进行估算，主要包括水量平衡法、降水存储量法等[3]，但上述方法难以获取较高精度数据，从而进行大区域估算时会产生较大误差。水质评价方法有单因子法、水质指数法、多变量分析法等[4]，这些方法多用来评价水质状况，无法针对水质净化功能进行评估[5]。近年来，随着计算机科学的不断进步，诸如水文评估模型Swat[6]、半分布式水文水质模型HSPF[7]、水文生态模型InVEST[8]等已被应用于产水服务与氮磷净化能力分布及动态特征研究。其中，由斯坦福大学开发的InVEST模型在国内外得到广泛应用[9]。Leh等[10]成功利用InVEST模型对西非国家加纳和科特迪瓦的产水量进行了评估，Goldstein等[11]利用InVEST模型分析了不同情境下的夏威夷岛屿水质净化功能空间分异情况；我国学者利用InVEST模型的产水模块分析了黄河流域[12]、太湖流域[13]、黑河流域[14]、三江源地区[15]、石羊河流域[16]等地区的产水量及水源涵养时空分布特征，并且利用InVEST模型的NDR模块(水质净化模块)对珠江流域[5]、海河流域[17]、官厅水库流域[18]、海南岛[19]等流域的氮磷净化功能开展研究，并取得较好的评估效果。

武烈河流域是承德地区主要供水水源地，也是京津冀水源涵养功能区重要组成部分，其水量与水质有重要的生态价值与研究意义，但是截至目前仍缺乏针对该区域产水服务和水质净化服务的研究。当前针对武烈河流域的研究主要集中于径流模拟[20]、水体污染特征[21]、水质污染源解析[22]等，但从生态系统服务角度对水量与水质的研究鲜见报道。笔者以武烈河流域为研究区域，基于InVEST模型产水模块与NDR模块，分析了1990—2018年产水服务及氮磷净化能力的时空分布特征，并利用情景模拟法量化了气候要素与土地利用/覆被变化对其影响的程度，以期为京津冀水源涵养功能区的合理规划提供科学支撑。

1 资料与方法

1.1 研究区概况武烈河流域位于河北省承德市境内，属滦河一级支流，古有热河之称，是承德市最重要的供水水源地之一，干流全长110 km，流域面积2 589 km2，主要由武烈河干流、玉带河、鹦鹉河、兴隆河这4条水系构成(图1a)，界于117°45′～118°26′E、40°54′～41°42′N；属于暖温带和寒温带过渡地带，四季分明，气候为典型大陆燕山山地气候，年均气温8.9 ℃，年降雨量为537.2 mm，年均径流量2.60×108m3，多年平均流量6.93 m3/s，降水的年内分配十分不均匀，降水期与汛期主要集中在6—9月，占全年降雨量的77.9%。流域内主要土地覆被类型为林地、草地，占总面积的 75.1%，其次为耕地与建筑用地，分别占总面积的19.3%、4.7%(图1b)，主要土壤类型为棕壤、褐土和草甸土。

图1 研究区地理位置(a)和土地利用类型(b)Fig.1 Geographical location(a)and land use type(b)of the study area

1.2 研究方法

1.2.1产水服务计算原理。InVEST模型产水量计算是根据水循环的原理，通过降水、植物蒸腾、地表蒸散发、根系深度、土壤深度和植被可利用水等参数经过模型计算获得。计算公式如下：

(1)

(2)

(3)

(4)

AWCxj=min(MSDxj,RDxj)×PAWCxj

(5)

式中，Yxj表示土地利用/覆被类型j中栅格单元x上的年产水量；AETxj为土地利用/覆被类型j上栅格单元x的年平均蒸散量；Pxj为栅格单元x土地利用/覆被类型j的年降雨量；Rxj表示土地利用/覆被类型j上栅格单元x的干燥指数；Kcxj为栅格单元x上土地利用/覆被类型j的作物系数；LAI表示植被叶面积指数；ET0表示潜在蒸散发量；Z为Zhang系数，是用于表征研究区降水特征的常数，取值为1～10[23]；ωxj是气候土壤非物理参数；AWCxj表示栅格单元x上土地利用/覆被类型j的有效可利用水；MSDxj为土层厚度最大值(m)；RDxj为根系深度(m)；PAWCxj表示栅格单元x上土地利用/覆被类型j的植被可利用水，计算方法采用周文佐等[24]提出的模型进行估算，公式如下：

PAWCxj=54.509-0.132SSAN-0.003(SSAN)2-0.055SSIL-0.006(SSIL)2-0.738CCLA+0.007(CCLA)2-0.688OM+0.501(OM)2

(6)

式中,SSAN为土壤砂粒含量(%)；SSIL表示土壤粉粒含量(%)；CCLA为土壤黏粒含量(%)；OM为土壤有机质含量(%)。

1.2.2氮磷净化能力计算原理。氮磷净化能力采用InVEST模型中的NDR模块(水质净化模块)进行计算，主要运用植被和土壤可以通过储存和转换等方式移除或减少径流中的营养盐污染物以达到净化水质的作用的机理，该模型只考虑非点源污染中的总氮(TN)、总磷(TP)，二者输出量越高，表明氮磷净化能力越低，用于评估生态系统中植被和土壤的水质净化服务，其主要算法如下：

ALVx=HSSx×polx

(7)

式中，ALVx是栅格x的调整后输出量，HSSx为栅格x的水文敏感性得分，polx为栅格x的输出系数。计算得到TN、TP输出量之后，再结合各土地利用/覆被类型对污染物的移除效率来计算得出营养物保持(持留)量。

1.2.3情景模拟法。基于情景模拟的方法，研究气候要素与土地利用/覆被变化对武烈河流域产水服务功能、氮磷净化能力的贡献程度。为排除异常气候数据的影响，采用1990—1999、2000—2009、2010—2018年3个时期的平均模拟产水量、TN输出量、TP输出量进行研究。设计情景如表1所示。

表1 气候变化和土地利用/覆被变化情景的各时段数据输入

假设各时期的气候要素按前一时期输入，土地利用/覆被按实际情况输入，计算得到气候变化情境下的模拟值，该值与真实情景模拟值的差即为气候变化的影响量；同理，假设各时期的土地利用/覆被按前一时期输入，气候要素按实际情况输入，土地利用/覆被变化情景下的模拟值与真实情景模拟值的差即为土地利用/覆被变化的影响量。基于此，气候要素与土地利用/覆被影响的贡献比值可以用以下公式进行量化：

(8)

(9)

(10)

式中，RWY、RTN、RTP分别表示气候要素与土地利用/覆被对产水量、TN输出量、TP输出量影响的贡献比值，该值越接近1表明气候要素与土地利用/覆被的影响程度相近，越大表明气候要素影响程度大于土地利用/覆被，越小表明土地利用/覆被影响程度大于气候要素；WYO、TNO、TPO分别表示真实情景下的各时期年均产水量(m3)、TN输出量(t)、TP输出量(t)；WYC、TNC、TPC分别表示气候变化情境下的各时期年均模拟产水量(m3)、TN输出量(t)、TP输出量(T)；WYL、TNL、TPL分别表示土地利用/覆被变化情景下的各时期年均模拟产水量(m3)、TN输出量(t)、TP输出量(t)。

1.3 数据来源与处理

1.3.1数据来源。该研究收集了承德市武烈河流域气象、土地利用/覆被、土壤属性、流域矢量边界、数字高程模型(DEM)、生物物理属性以及氮磷输出负荷数据等用于模型计算，同时还整理了武烈河水文站径流监测数据与河流断面总氮、总磷监测数据对模型结果进行验证。其中，DEM 数据来源于地理空间数据云平台(http：//www.gscloud.cn)，7期土地利用/覆被数据(分别为1990、1995、2000、2005、2010、2015、 2018年)来源于资源环境数据云平台(http：//www.resdc.cn)；土壤属性数据来源于世界土壤属性数据库(HWSD)；气象数据(包括降水、气温、气压、风速等)来源于国家气象信息中心(http：//data.cma.cn)与国家青藏高原科学数据中心(http：//data.tpdc.ac.cn)；生物物理属性以及氮磷输出负荷数据参考InVEST 用户指南和查阅研究区附近区域的相关文献[25-26]后估算获得，相关参数见表2。

表2 各土地利用/覆被类型生物物理属性及氮磷输出主要参数

1.3.2数据处理。降水数据通过使用ArcGIS对国家青藏高原科学数据中心下载的全国月均降水量栅格数据进行叠加、裁剪获取。土地利用与土壤质地经ArcGIS对世界土壤属性数据库进行投影变换、裁剪获取。集水区利用ArcGIS水文分析模块的工具，根据流域与子流域的边界将研究区 DEM数据分割成小流域，共产生4个子流域。植被可利用水是基于世界土壤属性数据库(HWSD)中研究区内土壤黏粒含量、土壤粉粒含量、土壤砂粒含量、土壤有机质含量由经验公式(6)计算得到。根据史晓亮等[20]对武烈河径流的模拟和产水量的估算，通过公式(2)、(3)、(4)对Zhang系数反复验证校准，最终确定Zhang系数为3.7。潜在蒸散发量采用经过系数校正后的Hargreaves公式进行计算获取[27]。

2 结果与分析

2.1 武烈河流域产水服务功能与氮磷净化能力的时空分布特征经分析，武烈河流域多年平均产水深度为104.22 mm，平均产水总量为2.69×108m3；多年平均TN、TP输出总量分别为291.52、28.22 t，单位面积多年平均TN、TP输出总量分别为1.17、0.11 kg/(hm2·a)。

为精确量化武烈河流域的产水服务功能与氮磷净化能力，将数据分为1990—1999、2000—2009、2010—2018年3个时期进行分析。从时间上看，产水服务功能方面(图2a)，1990—1999年产水量呈减小趋势，2000—2018年产水量呈增加趋势，总体呈微弱下降趋势，多年平均降幅为0.04×108m3/a；最大值出现在2016年，为6.67×108m3，最小值出现在2009年，为0.58×108m3。氮磷净化能力方面(图2b)，TN、TP输出量在1990—2009年变化趋势不明显，2010—2018年显著增加，总体呈增加趋势，增幅分别为0.02、0.06 t/a；TN输出量在2017年最大，为293.86 t，在2009年最小，为290.04 t；TP输出量在2017年最大，为29.66 t，在1999年最小，为27.74 t。

图2 1990—2018年武烈河流域年均产水量(a)与氮、磷输出量(b)变化Fig.2 Variation in annual average water yield(a)and TN and TP output(b)in Wulie River Basin in 1990-2008

从空间上看，产水服务功能方面(图3a)，产水深度呈东高西低、南高北低的空间分布特征，这主要是由降水、蒸散发、土地利用/覆被等因素的空间分异导致的；平均产水深度1990—1999年为139.65 mm，2000—2009年为58.19 mm，2010—2018年为104.38 mm。氮磷净化能力方面(图3b、c)，TN、TP输出量高的地区均分布在地势平坦和人口分布密集区，武烈河流域以种植农作物和果园为主，农药、化肥的使用是氮、磷的主要来源，有部分农业区距离河道较近，因此濒临河道处应扩大林地面积，进一步增强该区水质净化能力；年均TN、TP输出量1990—1999年分别为1.16、0.10 kg/(hm2·a)，2000—2009分别为1.15、0.09 kg/(hm2·a)，2010—2018年分别为1.18、0.12 kg/(hm2·a)。

图3 1990—2018年武烈河流域各时段年均产水量(a)与TN输出量(b)、TP输出量(c)时空分布变化Fig.3 Temporal and spatial distribution changes of the average annual water production(a)，nitrogen output(b)，and phosphorus output(c)in the Wulie River Basin from 1990 to 2018

2.2 气候因素与武烈河流域产水功能、氮磷净化能力的关系基于线性回归方法，分析1990—2018年武烈河流域各年均产水量以及TN、TP输出量与对应年份的降水量、潜在蒸散发量、温度之间的关系。结果表明(图4)，武烈河流域降水量与产水量呈显著正相关(P<0.05)，与TN、TP输出量均呈正相关，但不显著(P>0.05)；潜在蒸散发量与产水量呈显著负相关(P<0.05)，与TN、TP输出量均呈负相关，但不显著(P>0.05)；气温与产水量呈负相关，与TN、TP输出量均呈正相关，但相关性均不显著(P>0.05)。由此可见，武烈河流域的产水量与降水量(R2=0.987 8)、潜在蒸散发量(R2=0.539 5)具有较强的相关性，与气温存在一定相关性但不显著；流域的TN、TP输出量与各气候因素不显著。因此，各类气候因素中，降水量、潜在蒸散发量是影响武烈河流域产水量的主要因素，而对流域的TN、TP输出量影响相对较小。

图4 武烈河流域各气象因素与产水深度、TN输出量、TP输出量的关系Fig.4 Correlation analysis of each climatic factor and water yield depth，TN output and TP output in Wulie River Basin

2.3 土地利用/覆被变化对武烈河流域产水功能与氮磷净化能力的影响1990—2018年武烈河流域各类土地利用/覆被对产水量、TN输出量、TP输出量的多年平均贡献率与单位面积平均值表明(表3)，林地与草地的产水量贡献率最高，占比达73.33%，其次为耕地、建筑用地、水域；建筑用地的单位面积产水能力最强，为1 471.99 m3/(hm2·a)，耕地次之，其次为林地、草地、水域。耕地的TN、TP输出量贡献率最高，占比分别为83.15%、63.38%，其次为草地、林地、建筑用地、水域；耕地的单位面积TN、TP输出量也是最高的，分别为4.65、0.34 kg/(hm2·a)，其次为建筑用地、草地、水域、林地。由此可知，林地、草地虽然产水能力低于耕地和建筑用地，但因其面积占比高，故武烈河流域产水量主要来源于林地、草地；TN、TP输出能力最强的为耕地、建筑用地，原因是耕地化肥农药的使用增加了氮磷的输出量，同时建筑用地城镇化建设明显，人口和人类活动的增加也导致污染物排放递增，从而导致氮磷净化能力较弱，而林地、草地对污染物的截留效用较强，导致氮磷净化能力强。土地利用/覆被的变化会导致局部水汽循环、地表径流及污染物截留效用的变化，进而导致产水服务功能、氮磷净化能力的变化。随着武烈河流域城市化的加剧，大面积的林地、草地和耕地转化为建筑用地(表4)，导致其拦截降水与污染物截留效用的减弱，同时，这些区域更加频繁的人类活动也导致污染物排放的增加，从而导致产水服务功能增强及氮磷净化能力的减弱。

表3 武烈河流域各类土地利用/覆被对产水量、TN输出量与TP输出量的贡献率及单位面积平均值

2.4 气候和土地利用/覆被变化的影响程度量化分析气候要素与土地利用/覆被变化是导致武烈河流域产水服务功能、氮磷净化能力变化的主要原因，两者的影响程度存在差异，基于情景模拟法计算结果表明(表5)，1990—1999年，气候要素在产水服务功能、氮磷净化能力方面分别是土地利用/覆被影响程度的58.17、0.78、0.64倍；2000—2009年，气候要素在产水服务功能、氮磷净化能力方面分别是土地利用/覆被影响程度的26.17、0.49、0.11倍；2010—2018年，气候要素在产水服务功能、氮磷净化能力方面分别是土地利用/覆被影响程度的38.68、0.62、0.32倍。由此可见，1990—2018年在产水服务功能、氮磷净化能力方面，气候要素与土地利用/覆被影响程度比的平均值分别为41.01、0.63、0.36倍。气候要素变化对产水服务功能的影响程度高于土地利用/覆被变化，而土地利用/覆被变化对氮磷净化能力的影响程度高于气候要素变化。该结论与前人研究结果相似，如王亚慧等[28]研究表明横断山区的产水量与降水量存在极显著正相关；赵亚茹等[23]研究表明气候要素是石羊河上游产水量的主要影响因素；窦攀烽等[29]通过量化气候与土地利用/覆被对宁波地区产水量的影响，确定气候变化对产水量的贡献率在97%以上；张建等[2]研究表明丹江口市的氮磷输出量与土地利用/覆被的景观形状指数呈显著正相关；吴瑞等[18]研究表明官厅水库区域的氮磷净化能力的强弱与土地利用/覆被类型呈显著相关。

表4 武烈河流域各土地利用/覆被类型面积占比转移矩阵Table 4 Area proportion transfer matrix of various land use/land cover types in Wulie River Basin %

表5 气候要素与土地利用/覆被变化对武烈河流域产水服务功能、氮磷净化能力的影响程度

3 结论

该研究以武烈河流域为研究对象，运用InVEST模型对武烈河流域1990—2018年产水服务功能、氮磷净化能力的时空变化进行计算，分析了与不同气候因素的相关性，并量化了各类土地利用/覆被对产水量、TN输出量与TP输出量的贡献率，采用情景模拟法对影响武烈河流域产水功能与氮磷净化能力的因素进行量化。具体结果如下：

(1)1990—2018年，武烈河流域产水量呈现先减小后增加，总体呈微弱减小的变化趋势，多年平均产水量为2.69×108m3；TN、TP输出量呈增加趋势，多年平均TN、TP输出量分别为291.52、28.22 t。近29年来，武烈河流域产水服务功能与氮磷净化能力均减弱。

(2)武烈河流域产水量与降水量呈显著正相关，与潜在蒸散发量呈显著负相关，与气温呈负相关但不显著；TN、TP输出量与各气候要素呈一定相关关系，但均不显著。气候要素与产水量的相关性较强，与TN、TP输出量的相关性较弱。

(3)武烈河流域产水量主要来源于林地和草地，占比达73.33%，其次为耕地、建筑用地、水域；TN、TP输出量主要来源于耕地，占比分别达83.15%、63.38%，其次为草地、林地、建筑用地、水域。建筑用地产水能力最强，单位面积产水量为1 471.99 m3/(hm2·a)；耕地TN、TP输出能力最强，单位面积输出量分别为4.65、0.34 kg/(hm2·a)。

(4)气候要素变化对产水服务功能的影响更为显著，而土地利用/覆被变化对氮磷净化能力的影响更为显著。1990—2018年产水服务功能方面，气候要素是土地利用/覆被影响程度的41.01倍；氮磷净化能力方面，气候要素分别是土地利用/覆被影响程度的0.63、0.36倍。