李爱娟, 徐光来,3, 杨强强, 池建宇
(1.安徽师范大学地理与旅游学院, 安徽 芜湖 241003;2.安徽省江淮流域地表过程与区域响应重点实验室, 安徽 芜湖 241003;3.安徽师范大学 皖江流域退化生态系统的恢复与重建省部共建协同创新中心, 安徽 芜湖 241003)
生态系统服务是对人类的生存和生活质量做出贡献的生态系统产品和生态系统功能,人类从生态系统中获得的所有利益,包括供给服务、文化服务、调整服务和支持服务[1-2].随着经济发展人类对资源和环境的需求增加,生态系统的服务价值与人类的福祉有着更加密切的联系.生态系统服务已成为生态学、地理学研究的热点问题[3-4].
产水量是生态系统供给服务中重要的供给功能之一,受到气候、植被、土壤和地形等多种因素的影响,这些因素使得产水功能在空间上具有异质性[5].我国是水资源相对匮乏的国家,产水量对流域内及下游地区的生产、生活、生态都起着至关重要的作用.产水功能作为生态系统供给服务中的重要指标,受多种因素的影响,在时空上出现不平衡的情况.侯文娟等[6]构建乌江三岔河流域SWAT模型,结合地理加权回归模型分析了喀斯特流域产流服务的空间分异规律及与景观破碎化指数、NDVI等因子的关系.戴尔阜等[7]利用InVEST模型模拟横断山区的产水量空间分布,研究结果表明气候因子是影响产水量的空间异质性的主导因素,不同地貌及气候分区的因子解释力不同.
SWAT模型的水文模块可以综合流域多种因素影响得出空间上的产水量分布,将其用于产水服务模型构建可以科学评估流域产水量的时空分布规律.利用地理探测器对模拟产水模拟结果进行归因分析,可以从空间上探测变量可能存在的因果关系,科学定量研究产水量的影响因素,可以对流域的生态系统产水服务进行科学评价,为水资源的可持续利用和水资源调配管理方法提供参考,同时对流域生态补偿研究奠定基础,促进流域健康可持续发展.
选取淮河干流淮滨水文站以上流域为研究区,该区域总面积为141 12 km2,流经河南、湖北两省,流域范围主要位于信阳市,部分位于南阳市、驻马店市、随州市、孝感市及黄冈市.该流域内地形单元丰富,有低山、丘陵和相对平坦的平原地区,属于亚热带与温暖带过渡的气候,自然水系发育较好.南部及西部地区受人类活动影响较小,平原地区主要为信阳市辖区范围,城镇化的建设使得地面硬化率提高会对产水量产生一定的影响.流域内土地利用类型多样,以耕地和林地为主,流域内地理环境具有较大的空间异质性,如图1所示.
图1 研究区概况
本文中DEM(30 m×30 m),土地利用(30 m×30 m)及影响因子中人口密度数据(1 km×1 km)、NDVI数据均来自中国科学院资源环境科学与数据中心(http://www.resdc.cn/);土壤数据(1 km×1 km)来自世界和谐土壤数据库HWSD;降水、气温等构建模型的气象数据来自中国气象数据网(https://data.cma.cn/),分辨率为日.用于校准模型的实测径流数据来自河南省信阳市水文局.
借助ArcGIS将DEM、土地利用等空间数据库定义统一的投影坐标系.原始的18种土地利用类型进行适用于模型内部数据库的重分类,将研究区分成6种土地利用类型,分别为耕地、林地、草地、水域、居民区、荒地、裸地.利用HWSD土壤数据库,通过SPAW软件及土壤水文学分组进行一定的计算转换得到研究区的土壤数据库.将每个站点的降雨量、最高最低气温、风速、相对湿度及太阳辐射数据整理成SWAT模型所识别的数据格式,时间尺度为1981-2018年,步长为日,并通过站点经纬度与空间数据库连接.
SWAT模型对流域水文过程模拟主要包括陆地、汇流两个阶段的水循环.在陆地阶段的水循环中可以获得流域总径流量与产沙量,前者是通过演算划分的每个HRU的径流量得到,后者则由修正的通用土壤流失方程计算[8-9].模型的水量平衡公式如下[10]:
(1)
式中:SWt代表了土壤最后保持水量(mm);t代表了时间步长(d);SWo和Rday分别代表了第i天的土壤最开始保持水量和降水量(mm);Qsurt和Ea分别代表了第i天的地表产流量和蒸散发量(mm);Wseep和Qgw分别代表了第i天地下和包气带的保持水量(mm).
SWAT模型会在流域内根据土地利用、土壤、坡度等因素划分的具有相同水文特性的水文相应单元(HRU),再根据水量平衡法计算每个HRU内的产水量,即在给定期间内单位流域面积产生的注入水道的总水量,最终读取得出HRU尺度的产水量.本文利用SWAT模型水文模块对HRU尺度的产水量输出数据进行流域水供给量的计算.产水量计算公式如下[11]:
WYLD=SURQ+LATQ+GWQ-TLOSS-pond
(2)
式中,WYLD为总产水量(mm),SURQ为主水道中的地表径流量(mm),LATQ为主水道中的侧向流量(mm),GWQ为主水道中的地下径流量(mm),TLOSS为河床传输损失量(mm),pond为池塘截留量(mm).
产水的空间分异是由多种因素导致的,假设某因素对产水产生影响,则其空间分布就会和产水的空间分布有一定的相似性.本文运用地理探测器对产水服务进行单因子及交互因子探测.单因子探测因变量的空间分异性,明确不同影响因子对产水量空间分异的解释程度,而交互作用探测器则可以评估两个不同自变量同时作用时对产水量空间分异的解释能力.
某影响因子X对产水量的解释力为q,公式如下[12]:
(3)
(4)
土地利用类型改变了地表状况,地表状况影响了地面的蒸发量,从而对区域的产水量产生影响[13].植物的蒸腾作用是水循环的环节之一,其影响着水量的产生.地形和坡度是重要的下垫面因子,是气候形成的重要影响因素,直接影响着降水量及产水量[14].受人类活动的影响,土地利用景观格局发生显著变化,流域径流也因此受到影响,水文过程发生改变.因此综合考虑了自然和人为因子共9个因子进行分析,自然因子中选取降水、气温、DEM、坡度、NDVI等5个因子,人为因子中选取了人口密度、香农多样性指数(SHDI)、蔓延度指数(CONTAG)和有效粒度尺寸(MESH)4个因子分别表示景观的多样性、聚集度以及破碎度.对数据进行ArcGIS渔网重采样,按照1 km间隔生成,使各影响因子在空间单位上统一,便于进一步分析.
SWAT模型内参数众多,同一个参数在不同的研究区对模型的贡献率也不同.但是受限于实测数据的获取难度,很难对每一个参数每一个子流域进行分别校准.本文利用大坡岭长台关2001-2015年实测月径流数据对模型进行率定和验证,其中2001-2010年为率定期,2011-2015年为验证期.参考以往对淮河流域的研究以及自我调参的方式最终选取了8个参数来进行模型的校准.
率定方法采用SWAP-CUP软件中的反演建模算法SUFI2进行参数率定[15],本文修改后的方法用赋值法,率定时采用先上游后下游的顺序,经过200次迭代后得出最佳的取值,进而代入SWAT中进行重新模拟取值验证.大坡岭站和长台关站实测值和模拟值对比如图2和图3所示.经过计算得到大坡岭率定期确定性系数(R2)和纳什效率系数(NSE)分别为0.864、0.818;长台关率定期R2和NSE分别为0.873、0.829;大坡岭验证期R2和NSE分别为0.920、0.747;长台关验证期R2和NSE分别为0.862、0.719.总体来看,模拟效果较理想,率定后的模型在淮河上游具有较强的适用性.
图2 大坡岭月平均径流实测值与模拟值比较
图3 长台关月平均径流实测值与模拟值比较
淮河上游流域产水量的变化范围为8~531 mm,平均产水量为345 mm(图4).产水量整体呈现中间低四周高的分布格局,较高的地区集中在南部及东北部,其中产水量在500 mm以上的地区主要集中在南部的高程较高的地区,中部为低值区,产水量范围主要在300 mm以下.
图4 产水服务空间分布图
3.3.1 单因子归因分析
利用地理探测器对淮河上游流域的产水量空间分异影响因素进行探究.在选取的5个自然因子和4个人为因子中,降水和DEM对产水量空间分异解释力最高,分别达到了61.4%和44.2%.气温和坡度对产水量的空间差异的解释力差异不显著,坡度略高为28.5%,气温为27.7%.人为因子中解释力最高的为表征景观破碎度的有效粒度尺寸,解释力为40.2%,其次为人口密度,解释力为29.0%.所有因子中解释力排名靠后的分别为香农多样性指数的解释力22.8%、蔓延度指数的解释力16.0%以及NDVI的解释力10.0%.总体来看,自然因子的总体解释力大于人为因子,降水量及DEM对产水量的空间分异影响较大.
为探究不同因子在不同海拔高度上对产水量的影响程度,将淮河上游流域不同DEM分为7个不同分区并进行影响因子分析,不同DEM分区各因子对产水服务空间分异的影响程度(图5).由图5可以看出,同一因子在不同海拔分区的解释力存在较大差异.以降水量为例,降水量在大于400 m的地区中对产水量的解释力达到了88%,而在小于50 m的地区,解释力为31.1%.随着海拔的降低,降水量对产水量的解释力逐渐降低.坡度因子在50 m以下的地区对产水空间分布解释力不显著,而在200~300 m却可达26.7%,总体来说随着高程的上升,自然因子的解释力增强显著.在人为因子中,有效力度尺寸在100~200 m的分区中对产水量的解释力最小,为27.4%,在400 m以上的分区中解释力最高,为51.7%.总体来看,降水量对产水量的解释力在不同的DEM分区中的差异最大,即降水对产水量的影响受海拔的影响最大.气温及人口密度对产水量的影响受不同海拔的影响最低.自然因子在不同海拔上对产水量影响差异较大,人为因子在不同海拔上对产水量影响差异较小.
图5 不同DEM分区各因子对产水服务空间分异的影响程度
3.3.2 自然因子与人为因子的综合分析
产水量的空间结构是由多因素间相互作用共同决定的[5],因子之间的交互作用对产水量空间分布格局的影响高于单个因子的影响程度.在单因子的影响中,最高的为降水量的61.4%,在多因子交互作用中,流域内400 m以上的地区,降水和坡度的交互作用解释力可达96.4%,所有因子的交互作用解释力靠前且具有代表性的因子数据如图6所示.
图6 不同DEM分区多因子交互作用产水服务空间分异的影响程度
总体来看,降水因子和其他因子的交互作用对产水量的解释力最高,自然因子内部降水∩坡度的交互作用解释能力最强.人为因子中,人口密度∩景观破碎化程度的解释力较强.
在不同的地形分区中,主导的因子交互作用不同.在海拔相对较低的地区,主导的交互作用主要为降水∩人口密度.在50 m以下的地区解释力为93.7%,随着海拔的上升,解释力有所下降,在400 m以上的地区解释力为89.1%.降水及坡度的交互作用在海拔相对较高的地区对产水量的解释力最高,在400 m以上的地区降水∩坡度的解释力为96.4%,而在50 m以下的地区降水∩坡度的解释力为87.3%.次要的主导交互作用在海拔较低的地区为降水∩坡度,在海拔大于150 m的分区后,降水∩坡度成为主导交互作用,降水∩DEM为次要的主导交互作用.自然因子内部的交互作用解释力都呈现随着海拔升高而增加的趋势,人口密度∩景观破碎化程度的解释力大小受海拔的影响较小.
解释力较低的交互作用因子主要为温度与其他因子的交互及NDVI与其他因子的交互.温度与其他因子的交互解释力在10~50%之间,NDVI与其他因子的交互解释力20~65%之间.人为因子之间的交互作用的解释力也相对较低.
利用SWAT模型模拟的产水量为自变量,以2010年为代表年份,利用地理探测器,对淮河上游流域产水量的空间分异影响因素进行分析,主要结论如下:
1) 淮河上游流域产水量的变化范围为8~531 mm,平均产水量为345 mm.产水量整体呈现中间低四周高的分布格局,较高的地区集中在南部及东北部,中部为低值区.
2) 对产水量影响最大的为降水和DEM,解释力为61.4%和44.2%.人为因子解释力最高的为有效粒度尺寸,解释力为40.2%.自然因子的解释力大于人为因子.对于不同海拔的分区,同一个因子的解释力有显著差异.降水量在大于400 m的地区中对产水量的解释力达到了88%,在小于50 m的地区,解释力为31.1%.随着高程的上升,自然因子的解释力增强显著.
3) 因子之间的交互作用对产水量空间分布格局的影响高于单个因子的影响程度.在海拔相对较低的地区,主导的交互作用主要为降水∩人口密度,在50 m以下的地区解释力为93.7%.在海拔大于150 m的分区后,降水∩坡度为主导交互作用,降水∩DEM为次要主导交互作用.