基于InVEST模型的大理州永平县水生态系统服务功能时空变化分析*

杜佳衡，王锦

(西南林业大学，云南 昆明 650224)

水生态系统作为自然生态系统中的一个重要类型，为人类的生存繁衍提供了多种物质产品及资源，具有多种服务功能，如集水、产水、水质净化，同时兼具一定的美学和娱乐游憩价值[1-3]。其中，产水和水质净化作为水生态系统主要的服务功能，具有提高民生福祉及改善城镇居民饮用水安全的重要作用[4-5]。然而，在当前快速城市化进程中，人类社会生产活动对水生态系统的干扰逐渐增强，导致其服务功能日益衰退[6]。为明确区域水生态系统的健康状况，开展产水和水质净化服务功能的空间可视化定量分析，揭示两者的时空变化分布特征，对于实现区域水资源管理和可持续发展意义深远。

生态系统服务自20世纪70年代提出以来，成为生态学、生态经济学等学科领域中的重要研究方向。在SCEP(study of critical environmental problems)[7]、Costanza等[8]、Daily[9]的研究基础上，生态系统服务的概念、分类和评估方法逐渐被明确，水生态系统服务功能的相关研究也越来越广泛[10]。从研究对象和尺度上看，以往学者主要以流域[11-12]、湿地[13]、河流[14]等类型来进行水生态系统服务功能评估研究。从研究方法来看，目前针对水生态系统产水服务的研究主要基于水量平衡法、蓄水能力法等，运用实测数据开展产水量的评估[15-16]。水质净化服务往往通过单因子法、灰色综合评价法等方法分析评估区域水质状况[17]。以上方法受水文站点实测数据影响较大，仅适用于小尺度范围内的研究，且多采用价值量评价法进行估算，结果未能实现空间化。近年来，在生态系统服务物质量评价法的基础上，随着ArcGIS、RS等技术的迅速发展，对于生态系统服务“模型化、空间化”的研究成为新趋势，各类水文研究模型随之得到广泛应用，包括InVEST模型、SWAT模型等[18]。其中，InVEST模型以简便及快速的空间表达手段成为当下最热门的研究工具，越来越多的学者借助该模型的产水模块(water yield model)和水质净化模块(nutrient delivery ratio model)进行区域水生态系统服务功能的研究。如白杨等[19]、陈龙等[20]利用InVEST模型分别对美国肯塔基州和澜沧江流域的水源涵养、水土保持、水质净化功能进行了评估；吴瑞等[21]、韩会庆等[17]分别对官厅水库流域和珠江流域的产水和水质净化服务进行了模拟，并对其时空分布差异性进行驱动因素的探究。综上所述，InVEST模型因易操作、数据易获取等特点已被国内外学者广泛应用于不同尺度的研究当中，但目前针对国内西南部山区县域尺度的相关研究仍然较少。

永平县地处云南大理州西部山区，近年来，当地发展面临着诸多压力和挑战，特别是开发活动的不断扩张造成了区域内资源的快速消耗甚至是生态破坏。为缓解当地工农业生产与自然过程之间的矛盾，有必要对其水生态系统产水服务和水质净化服务功能展开分析与评估，而目前，该区域内对于这两种服务功能的研究尚未形成一套较为完整的研究思路及方法。因此，本文基于InVEST模型产水量和水质净化量两个子模块模拟2000年、2010年、2020年大理州永平县的年均产水量及年均氮磷输出量，分析这两个水生态系统主要服务功能的时空变化，并划分重要性等级，以期为当地经济健康发展、水源水质保护提供科学依据。

1 研究地概况

永平县(99°17′～99°56′E、25°03′～25°45′N)位于云南省西部，总面积约为2 790 km2，属大理州下辖县。县内与云台、博南两山相连，地势狭长且高差明显，自西北向东南降低，海拔最高2 974 m，最低1 105 m。银江河分布于永平县中部成为主要干流，澜沧江、顺濞河经县域东西两侧穿越而过，形成“两山三河”的天然山水格局和山、谷、坝错综复杂的地形地貌[22]。水资源丰富，平均水资源总量超过10×108m3，人均水资源量约7 000 m3，有开发利用价值的河流27条(图1)。

图1 研究地概况

2 研究方法

2.1 数据来源及预处理

InVEST产水(water yield model)和水质净化(nutrient delivery ratio model)模型运行参数主要包括以下内容。

(1)土地利用数据 数据来自资源环境数据云，选取2000年、2010年、2020年3期栅格数据，统一分辨率为30 m。根据刘纪远等[23]提出的体系，按照二级分类标准将研究区土地利用划分为5大类12小类(图2)。

图2 大理州永平县土地利用类型

(2)气象数据 包括降水、气温等日值数据，均来自中国气象数据网，包括研究区及其周边28个气象站点的数据，为使数据具备代表性，分别选取1996—2000年、2006—2010年、2016—2020年3个时段的数据滑动平均值，表征2000年、2010年、2020年的气象数据。其中，年均降水量在ArcGIS中通过反距离权重法(IDW)处理插值后获得栅格数据；潜在蒸散发量在参考FAO推荐的Penman-Monteith公式[24-25]计算后，通过IDW插值为栅格数据。

(3)DEM 数据来自地理空间数据云，经投影、填洼处理后借助水文分析工具提取流域及子流域。

(4)土壤数据 数据来自中科院南京土壤所，土壤深度数据在ArcGIS中经投影变换、裁剪处理获取栅格数据，植被可利用含水率数据参考周文佐等[26]拟合的经验公式计算提取。

(5)生物物理系数表 参考前人研究成果[27-28]及模型用户使用手册[29]赋值设定(表1)。

表1 InVEST模型运行参数生物物理系数表

(6)经验参数 产水模块Z系数通过《大理白族自治州水资源公报》调整设定。由于数据获取限制，仅收集到研究区2010年水资源统计成果。基于此对产水模块输出结果进行校验，发现当Z系数赋值为30时，模型的产水量更接近实测结果。其他经验参数如水质净化模块水文连通性K系数等使用模型默认值设定。

2.2 产水服务模拟

InVEST 产水模型(water yield model)是基于Budyko水热耦合平衡原理[30]提出的一种水量平衡估算方法，即以各栅格单元降水量与实际蒸散量之差作为该栅格单元产水量，模型不考虑地下水的补给。产水量计算公式如下：

式中：Y(xj)为地类j在栅格单元上x的年均产水量(mm)，AET(xj)为栅格单元x上地类j的年均实际蒸散量(mm)，P(x)为栅格单元x上的年均降水量(mm)。其中，AET(xj)/P(x)为Zhang等[31]在参考Budyko曲线后提出的方程，计算公式如下。

PET(xj)=Kc(lc)×ET0(x)

AWC(c)=Min(Rest.layar.depth,root.depth)×PAWC

AET(x)=Min〔Kc(lc)×ET0(x),P(x)〕

式中：PET(xj)为地类j在栅格单元x上的年均潜在蒸散发量(mm)，ω(x)为自然气候—土壤属性的非物理参数；ET0(x)为参考植被在栅格单元x上的蒸散量(mm)，Kc(lx)为地类x的植被蒸散系数；AWC(x)为栅格单元x的土壤有效含水率(mm)，PAWC为植被可利用水；Z即Zhang系数，季节常数，取值在1～30之间。

2.3 水质净化服务模拟

水质净化服务通过InVEST NDR Model进行计算，模型基于土地利用/覆被类型在年尺度上计算栅格单元氮、磷营养盐的保持量和输出量，利用氮、磷含量表征水质状况，TN、TP输出量越大表明栅格单元水质净化能力越弱。模型仅考虑植被和土壤对氮、磷营养盐非点源污染的截留及过滤功能，不涉及化学和生物之间的相互作用。水质净化服务计算公式如下：

Xexpj=loadsurf,j×NDRsurf,j+loadsubs,j×NDRsubs,j

式中：xexptot为流域年均氮磷营养盐输出总量(kg/a)；xexpj为栅格单元j的氮磷营养盐输出量(kg/a)，loadsurf,j和loadsubs,j为栅格单元j地表和地下氮磷营养盐负荷量(kg/a)，NDRsurf,j和NDRsubs,j为栅格单元j地表和地下氮磷营养盐输送效率。其中，氮磷营养盐负荷量计算公式如下：

loadsurf,j=(1-proportion_subsurfacej)×modified.load_nj

loadsubs,j=proportion_subsurfacej×modified.load_nj

modified.loadxj=loadxj×RPIxj

式中：modified.loadxj为修正后栅格单元j的氮磷营养盐负荷量(kg/a)，RPIxj为径流代理指数，RPj为营养盐径流代理，RPav为平均径流代理。

氮磷营养盐输送效率计算公式如下：

式中：NDRsurf,j为地表营养盐输送效率；IC0和k是校准参数；ICj为地形指数；NDR0,j为栅格单元j和径流之间的养分最大持留效率；effLULCj为地类j的最大持留效率；eff’downj为栅格单元j的下游有效持留效率；sj为步长因子；ljdown为栅格单元j到下游相邻栅格单元流动的路径长度；lLULCj为栅格单元土地利用类型j的持留距离；Dup和Ddn为上下坡贡献面积的平均坡度；NDRsubs,j为地下营养盐输送效率；effsubs为通过地下径流可达到的最大营养盐持留效率；lsubs为地下营养盐流动持留效率；li表示从栅格单元到径流的距离。

2.4 冷热点分析

3 结果与分析

3.1 产水服务时空变化模拟

基于InVEST模型产水(Water Yield)模块运行结果，大理州永平县2000年、2010年、2020年产水量空间分布格局如图3所示。

图3 大理州永平县产水服务时空分布特征

20年间，研究区产水量总体趋于下降趋势，不同年份间，产水量的空间分布格局基本保持稳定，大体呈现出西南高东北低的态势。3个时期流域的平均产水深度分别为642.37、348.95、275.59 mm，高值区集中于龙门乡东北部、博南镇中部、北斗乡西南部；较高值区集中于杉阳镇西部、龙街镇西部、厂街乡中部、水泄乡中北部；低值区集中连片分布于北斗乡中部、龙街镇东部、博南镇西北部和南部。

同一年份中，单位面积产水深度在空间上的分布具有异质性。按照一级分类标准，对土地利用重新分类后进行分区统计，结果如图4a所示。从各地类上看，20年来建设用地产水量始终占据高值，且始终高于流域年均产水深度，这可能与快速城市化进程改变的不透水下垫面有关，这一过程改变了水量的产汇流，使得地表产流增加而入渗减少，产水能力相对较高；耕地产水能力次之，这可能与人类引水灌溉等农耕干扰行为存在一定的关联性；林地和草地产水能力中等，这可能是由于植被及其枯落物对降水和地表径流的拦截作用所导致的，同时受到水量入渗土壤的影响，削减了洪峰流量，因此产水能力相对较弱；而水域的产水量始终处于低值，这是由于水域的蒸散发量较高，大于植被的蒸散发量，因此产水能力最弱。

从乡镇尺度上看，3个时期产水深度如图4b所示，产水深度最高的乡镇分别为龙门乡、杉阳镇、水泄乡；产水深度低值区则集中在北斗乡、龙街镇、博南镇。考虑到土地利用类型的影响，这种分布空间差异可能与各时期乡镇间的用地类型比例有关，可以看出杉阳镇、水泄乡、龙门乡等用地类型以耕地为主导的乡镇产水量最高，而北斗乡、龙街镇、博南镇等以林草地为主要用地类型的乡镇产水量则相对较少。同时，各时期气候变化导致的降水量、蒸散发量的大小也会对产水深度产生一定影响，从而影响产水量的时空分布格局特征。

图4 大理州永平县各地类、各乡镇产水量

3.2 水质净化服务时空变化模拟

经过InVEST模型(NDR)模块计算，大理州永平县2000年、2010年、2020年水质净化能力如图5所示，TN、TP输出量越大则水质净化能力越差。结果显示，20年来研究区TN、TP总输出量先增大后减小，总体处于增加态势，表明研究区水质净化能力总体趋于减弱趋势。其中，3个时期流域尺度TN年均输出量最大值分别为5.90、6.10、6.40 kg/(hm2·a)，总输出量分别为840.75、880.49、877.66 t；TP年均输出量最大值分别为0.14、0.15、0.16 kg/(hm2·a)，总输出量分别为78.46、79.46、79.34 t。从空间上看，3个时期氮磷营养盐输出空间分布格局基本保持一致，高值区主要集中于龙门乡中南部、博南镇中部和西南部、杉阳镇中东部、厂街乡西部且多呈面状分布，其余高值区以点状分布于水泄乡东北部、龙街镇中南部、北斗乡西北部。

图5 大理州永平县水质净化服务时空分布特征

从土地利用类型上看，3个时期TN年均输出量高低值依次为耕地>建设用地>草地>林地>水域；TP输出量依次为耕地>建设用地>林地>草地>水域。从各时期来看，耕地始终为大理州永平县氮磷营养盐输出总量的主要贡献角色，且输出量呈现出逐年增加的趋势，这是由于农耕地所种植的作物生长长期依赖于农药和化肥的使用，而农耕地多呈面状分布于河谷地带，因此成为氮、磷营养盐污染物的主要来源；建设用地TN、TP输出量次之，这可能是由于研究区城镇建设及农村养殖畜牧业发展所导致的；林地、草地氮磷营养盐输出量相对较少，这是因为植被具有过滤和截留营养盐的功能；而水域作为人类生存发展所需的重要资源，则不具备氮磷营养盐输出的能力(图6)。

图6 大理州永平县各地类TN、TP输出量

从不同乡镇间来看，杉阳镇、厂街乡、博南镇氮磷营养盐输出量始终高于其余4个乡镇，从土地利用类型的分布格局进行分析，杉阳镇、厂街乡、博南镇的地类以耕地和建设用地为主导，因此三者的营养盐输出量较高；北斗乡、龙门乡、龙街镇的氮磷输出量相对较少，这是因为这些区域的植被覆盖较为完好，从而对营养盐的输出起到一定的过滤及净化效用(图7)。

图7 大理州永平县各乡镇TN、TP输出量

3.3 产水和水质净化服务冷热点分析

在进行冷热点分析之前，首先运用ArcGIS空间自相关分析判断大理州永平县产水量和水质净化量在空间上是否存在聚集特征，计算结果如表2所示。其中，2000年、2010年、2020年研究区的产水量和氮磷输出量莫兰指数均大于0，表明两项服务功能均具有空间相关性，同时，Z值均大于2.58，概率P值为0，均小于0.01，表明研究区产水量和氮磷输出量的空间分布为非随机状态，具有聚集特征。基于此，开展大理州永平县产水和水质净化服务冷热点分析，结果如图8、图9所示。

表2 大理州永平县产水和水质净化服务功能空间自相关指数

图8 大理州永平县产水服务冷热点空间分布

图9 大理州永平县水质净化服务冷热点空间分布

3个时期，研究区两项服务的冷热点空间分布格局特征具有相似性。其中产水服务冷热点分析结果显示，热点区域主要分布于杉阳镇、龙门乡、博南镇、水泄乡，这些区域地类多以耕地、建设用地为主，冷点区域主要集中在北斗乡，该区以林草地为主要用地类型，植被覆盖度高，人为干扰活动较少，有利于抑制地表径流；水质净化服务热点区域则主要分布于杉阳镇、博南镇、厂街乡等乡镇间，这些区域用地多以耕地为主，氮磷输出量较高，北斗乡则为研究区的主要冷点区域。总体上看，20年来大理州永平县产水服务的热点区域从空间上看主要集中分布于西部地区，并逐渐从北部地区向南部地区发生转移，冷点区域主要集中分布于中部和东北部地区；水质净化服务热点区域分布较为稳定，主要集中在西北部和东南部地区，呈现出与耕地和建设用地一致性的空间分布格局特征。

3.4 大理州永平县水生态系统服务功能重要性等级划分

为划定大理州永平县水生态系统服务功能重要性等级，本文仅以2020年产水和水质净化服务作为划分对象。研究运用ArcGIS将产水服务和水质净化服务计算结果进行归一化处理，并进行空间叠加分析，将研究区水生态系统服务功能重分类为极重要、较重要、一般重要、不重要四个等级，分别归为极重要区、较重要区、一般重要区和不重要区[33-35]，面积分别为103 089.87、108 449.37、43 441.38、14 644.8 hm2，占总面积的比例为38.24%、40.22%、16.11%、5.43%，结果如图10所示。

图10 大理州永平县2020年水生态系统服务功能重要性等级分布

由图中各区空间分布特征可知，水生态系统服务功能极重要区主要分布在研究区中部和东部区域，用地以林草地为主，该区植被覆盖度高，具备良好的水源涵养功能，人为干扰活动少，是继续推进实施退耕还林还草的重要区域，同时应加强该区内水源地的保护工作，杜绝毁林开荒、毁林采矿等行为；较重要区和一般重要区主要分布在用地以灌草和旱地为主的西北部和南部区域，这些区域应加强退化草地的恢复治理力度，控制好环境污染，切实保护水源水质，减少对该区的干扰；不重要区分布在研究区用地类型以耕地和建设用地为主的区域，该区植被覆盖度低，人为干扰活动强烈，是重点实施退耕还林还草和污染防治的区域，以提升水生态系统的水源涵养能力。

4 讨论与结论

4.1 讨论

本文通过运用InVEST模型产水模块和水质净化模块对大理州永平县产水量和水质净化量时空分布变化特征进行模拟，并利用冷热点分析工具识别两项服务的空间分布聚集情况。研究结论与已有的研究相类似[36-39]，结果均呈现出产水服务以建设用地和耕地占据高值，而林草地占据低值的特征，这是由于快速城市化进程所改变的不透水下垫面导致的[40]。因此，当地应加强林草地种植及恢复力度，提高水生态系统水源涵养能力；水质净化服务中，耕地的氮磷输出量始终最高，该结果强调了农药和化肥的大量使用是造成区域水质净化服务功能衰退的主要因素，实施退耕还林还草工程是提高当地水质净化能力的有效措施[41]。同时，在氮磷输出量高值区，应多种植喜氮、喜磷植物，充分发挥植物的拦截与净化效用，合理使用农药和化肥，制定排污标准，降低氮磷进入水体的概率，确保当地居民饮用水安全。此外，通过模拟大理州永平县两项服务功能的时空变化，本文将水生态系统服务功能划分为4个重要区，确定优先重点保护的区域，以此为当地开展水源水质保护提供科学指导。其中，极重要区植被覆盖度高，具有良好的调节地表径流、水源涵养能力，是当地优先考虑实施保护的区域，应加强区域内退耕还林还草工程的建设工作，以此促进水源地的健康与稳定。

4.2 结论

(1)2000年、2010年、2020年大理州永平县的产水服务年均产水深度分别为642.37、348.95、275.59 mm，呈下降趋势。建设用地产水量始终占据高值，其次为耕地，水域产水量最低；用地类型以耕地为主导的乡镇产水量最高。

(2)2000年、2010年、2020年大理州永平县的水质净化服务TN年均输出量最大值分别为5.90、6.10、6.40 kg/(hm2·a)，TP年均输出量最大值分别为0.14、0.15、0.16 kg/(hm2·a)，均呈上升趋势。耕地氮磷输出量最高，其次为建设用地；以耕地和建设用地为主要用地类型的乡镇氮磷输出量较高。

(3)大理州永平县产水服务热点区域主要集中在西部地区，并逐渐从北部地区向南部地区发生转移，冷点区域主要集中分布在中部和东北部地区；水质净化服务热点区域呈现出与耕地和建设用地一致性的空间分布规律，冷点区域则集中连片分布于植被覆盖度高，较少人为活动干扰的区域。

(4)大理州永平县水生态系统服务功能极重要区占研究区总面积的38.24%，主要分布在用地以林草地为主的中部和东部区域；较重要区和一般重要区主要分布在西北部和南部用地以灌草和旱地为主的区域，占总面积的40.22%、16.11%；不重要区占总面积的5.43%，分布在用地以耕地和建设用地为主的区域。