乌江流域是我国长江中上游重要的生态屏障和淡水资源补给地, 同时也是西南地区的重要经济区域, 拥有丰富的水资源、生物资源和生态系统, 对维持区域水源涵养、水质净化和水土保持服务等方面具有重要的作用^[1—2]。乌江作为长江上游南岸最大支流, 优良的流域生态环境对保障和贯彻落实国家“长江经济带战略”生态优先、绿色发展理念发挥着重要的支撑作用, 对保障长江经济带的生态环境安全和可持续发展具有十分重要的意义^[3]。

近年来, 人类已逐步认识到水源涵养与水质净化服务在区域生态安全、环境净化方面的重要性^[4—5]。城市化进程、资源消耗、矿山开采等导致植被资源减少, 土壤质地改变, 使得原本由植被、土壤等保持的水分和吸收转化的氮(N)、磷(P)等营养物质量减少, 水环境中营养物质浓度上升, 区域水源涵养与水质净化服务的加速恶化^[6], 因此, 如何科学评价水源涵养与水质净化服务时空差异及驱动因素一直是生态学研究的难点。越来越多的学者对生态系统的水源涵养与水质净化功能开展了分析, 具备较好研究基础。国外学者用InVEST模型分析了英国小流域^[7]、非洲西部流域^[8]、威拉米特河流^[9]、西班牙弗兰科利河流域^[10]等区域水源涵养功能, 国内学者分别对三江源国家公园^[11]、祁连山区^[12]、闽三角城市群^[13]、官厅水库^[14]、汉江流域^[15]等区域水源涵养与水质净化进行了评估, 国外学者多关注区域水源涵养量的评估及其气象、水文因素分析^[16], 国内学者多数研究集中于国家公园、城市群、小流域等此尺度以及单一生态系统水源涵养与水质净化时空分布特征分析较多^[17—18]。缺乏对气候、地形地貌复杂的较大尺度流域以及不同土地利用类型之间的时空响应关系分析, 对水源涵养与水质净化的时空差异表征缺少定量化工具进行统一描述, 特别是乌江流域这类山地河谷地貌发育强烈, 地表破碎的流域。因此, 开展乌江流域水源涵养和水质净化服务的时空变化及土地利用变化响应研究十分必要。

InVEST模型广泛应用于生态系统服务评估集成, 其机理明晰, 相较于其它方法, InVEST模型优势在于该模型对评价区域内海拔、坡度等环境因子敏感性较高, 且充分考虑了气象和土地覆盖等因素, 模拟结果数值曲线较为平缓, 在模拟结果的空间化上具有优势, 适宜用于大尺度以及地形地貌空间差异较大的流域开展生态系统服务模拟研究^[19]。因此, 本文以乌江流域为研究区, 探讨水源涵养与水质净化时空差异以及影响因素, 科学认识乌江流域土地利用变化过程中对水资源生态服务的影响。以期为典型南方山区流域土地资源优化配置、水污染控制、水生态修复措施提供科学依据。

1 研究区概况

乌江发源于贵州省西北部威宁县乌蒙山东麓, 在重庆市涪陵区注入长江, 是贵州省纵跨经度最大的河流, 也是长江上游南岸最大的一级支流和重要的生态屏障, 本文研究的乌江流域范围为贵州省境内的范围, 流域总面积66303.06 km², 如图 1所示。乌江流域(贵州段)位于省内喀斯特地貌核心区域, 地处乌蒙山区、武陵山区及滇桂黔石漠化区三大集中连片贫困区, 大部分区域属于石漠化敏感地区, 是我国西部四大生态脆弱带之一。地理位置介于东经104°19′—108°47′E, 26°11′—29°13′N之间, 西南高, 东北低, 地表崎岖破碎, 地貌以内陆喀斯特山地丘陵为主, 落差在2600 m以上, 气候属亚热带湿润季风气候类型, 自然景观垂直变化明显, 年均温度11℃—18℃之间, 多年平均降水量在700—1800 mm之间。乌江流域共涉及贵州省7个地州市42个区县, 2015年流域总常驻人口2134.29万人, 其中61.3%为农业人口。矿产资源丰富, 目前已发现的矿产资源有110种以上, 主要以煤矿、磷矿加工(磷化、西洋肥业、中化开磷)等粗放式采掘业和制造业为主, 流域水环境和生态压力大。

2 研究方法与数据来源

研究基于InVEST模型的水源涵养和水质净化模块分析流域水源涵养和水质净化服务功能, 模型由美国斯坦福大学环境森林研究所、世界自然基金和大自然保护协会共同研发。水源涵养模块基于水量平衡的原理, 在估算流域产水量的基础上考虑气候、土壤、地形、土地利用状况因素, 计算水源涵养量^[20]。水质净化模块是评估分析流域范围内的N、P营养物质从外界进入地表后, 模拟经汇流流向河流湖库的空间运输过程。

2.1 Aunsplin气象插值

乌江流域峡谷分布广泛, 区域气候空间变化幅度大, 降水时空分布差异明显, 传统的气象插值方法误差较大^[21—22]。相关研究结果表明, Anusplin气象插值方法在复杂地形区域的插值精度优于反距离权重、克里金、协同克里金等传统气象插值方法^[22—24]。因此, 本研究通过Aunsplin基于普通薄盘和局部薄盘样条函数插值方法, 引入地形因子作为协变量进行降水插值获得研究区降水空间分布, 如图 2所示。

2.2 水源涵养评估 2.2.1 产水量计算

InVEST模型基于水热耦合平衡模型(Budyko)来计算流域产水量, 定义区域产水量由每个栅格单元的降水量减去蒸散量(包括地表蒸发量、土壤含水量、冠层截留量以及枯落物持水量)后剩下的水量, 未严格区分地表水、地下水和基流, 每个栅格单元产水经汇流后达到流域出口, 公式如下：

(1)

式中, Y_xj为栅格单元年产水量(mm), 包含地表径流、土壤含水、冠层截留水量; P_x为栅格单元x的年降水量(mm); AET_xj为土地利用类型j中栅格单元x的实际蒸散量(mm); AET_xj/P_x为实际蒸散量与降水量比值, 由公式(2)计算。

(2)

(3)

(4)

式中, R_xj为土地利用类型j中栅格单元x上的Budyko干燥系数, 定义为潜在蒸发量与降水量的比值, 无量纲, 由公式(3)计算; ω_x为修正植被可利用水量与降水量的比值, 与区域自然气候和土壤背景性质有关, 无量纲; Z是Zhang系数, 是表征季节性降水的特征常数^[25]; AWC_x为可利用水量(mm), 由土壤的深度和理化性质决定, 由公式(5) 计算, k为植被蒸散系数, ET₀为年潜在蒸发量(mm); 由式中(7)计算。

(5)

(6)

(7)

式中, SoilDepth_x为最大土壤深度; RootDepth_x为植物根系深度, PAWC_x为植被可利用有效含水率; sand为土壤砂粒含量(%); silt为土壤粉粒含量(%); clay为土壤粘粒含量(%); OM为有机质含量(%)。; R_n为净辐射MJ/(m²·d); G为土壤热通量; γ为干湿表常数; T为日均温度℃; μ₂为2 m高处风速(m/s); e_s和e_a分别为饱和水汽压和实际水汽压(kPa); Δ为饱和水汽压曲线斜率(kPa/℃)。

2.2.2 水源涵养计算

(8)

(9)

式中, WR_ij为j类土地利用上栅格单元i的年持水量, Y_ij为j类土地利用上栅格单元i的年产水量, Runoff_ij为j类土地利用上栅格单元i的地表净流量, P_ij为j类土地利用上栅格单元i的降水量, C_j为j类土地利用上栅格单元i的地表径流系数。

2.3 水质净化评估

模型水质净化营养物输出量的表达公式如下：

(10)

(11)

式中, X_{exp, i}为流域范围中每个栅格单元i的营养物输出量, load_{surf, i}为地表营养物的负荷量, NDR_{surf, i}为地表营养物的传输率, load_{subs, i}为地下营养物的负荷量, NDR_{subs, i}为地下营养物的传输率, X_exptot为区域划分子流域的营养物总输出量。

(1) 地表营养物计算公式

(12)

(13)

(14)

(15)

(16)

(17)

(18)

式中, NDR_{0, i}是指下游栅格单元未保留的营养物传输率, 与栅格单元在地表上的空间位置无关, IC_t是区域地形指数, IC₀和k是矫正系数; eff′_i为营养物在地表栅格单元i传输到河流之间的最大截留效率; eff_LULCi是土地利用类型i到达流的最大截留效率, eff′_downi是下游栅格单元i上的有效截留效率, s_i为步长因子; l_idown是栅格单元i到流域下游地区相邻栅格的路径距离, l_LULCi是土地利用类型栅格单元i的有效截留距离; D_up是上坡区域的平均斜率梯度(m/m), A是上坡贡献面积(m²), d_i是区域坡度最大的栅格下坡方向沿第i个单元的流动路径距离, S_i是第i个单元的斜率梯度。

(2) 地下营养物传输率计算公式

(19)

式中, eff_subs是从地表栅格单元能够下渗到地下的营养物最大截留效率, i为栅格单元, l_subs是地下河流的截留长度, 即土壤能够保持营养物最大容量所需的距离, l_i是栅格单元到地下河流的距离。

水质净化模块输入参数包括研究区数字高程模型(DEM)(填洼后)、土地利用类型、子流域范围、流域流量累积阈值、营养物径流代理、N、P营养物输出负荷和最大滞留距离、植被滞留效率以及Boreselli k参数。其中子流域采用ArcGIS水文分析工具对流域范围和子流域边界进行划分, 共分割出52个子流域, 作为InVEST模型在流域尺度水源涵养以及N、P营养物输出量的统计单元; 营养物径流代理采用研究区年总降雨量代替, 流域流量累积阈值与Boreselli k参数经过多次调试试验, 最终确定流域流量累积阈值取值为500, Boreselli k取值2。根据自然环境的相似性, 参考InVEST模型手册和相关学者的研究成果^[26—28], 本研究氮磷营养物输出负荷系数、植被滞留效率取值见表 1。

2.4 Pearson相关分析

本研究借助SPSS统计软件, 使用Pearson相关系数研究不同土地利用类型与水质净化强度之间的相关性, Pearson相关性分析公式如下：

(20)

式中, R为相关系数, x和y分别为变量x和y的均值, x_i和y_i分别量x和y的第i个观测值, n为样本数量, R取值范围为-1—1之间, 绝对值越大, 土地利用类型与水质净化变量之间相关性越强, R>0表明变量之间呈现正相关关系, R < 0表明变量之间呈现负相关关系。

2.5 数据来源

(1) 土地利用数据, 空间分辨率30 m, 解译于各期Landsat TM5、ETM+、OLI和环境减灾卫星HJ—1影像; (2)DEM数据, 空间分辨率为30 m, 来源于地理空间数据云(http://www.gscloud.cn/); (3)气象数据, 来源于中国气象科学数据共享网; (4)土壤数据, 来源于联合国粮农组织和维也纳国际应用系统研究所构建的世界土壤数据库, 研究区为南京土壤所第二次全国土地调查数据(1 ∶ 100万); (5) 基础地理数据, 包含研究区喀斯特分布数据、行政区划、河流水系等。

各类数据利用研究区矢量边界对其进行裁剪, 空间分辨率均采样为30 m, 均统一转换为WGS84/Albers Equal Area Conic投影参与空间计算。

3 结果分析 3.1 土地利用时空变化特征

乌江流域2000年、2010年和2020年土地利用空间分布及数据统计如图 3、表 2所示。从整体上看, 2000—2020年, 流域土地利用类型呈现以“耕地-林地”为主导的地域特征, 其中耕地占据绝对主导地位, 林地相对次之, 未利用地的面积最少。以2020年为例, 耕地和林地的面积分别达到28817.39 km²和24235.68 km², 占流域总面积的80.02%。耕地集中分布在流域中下游及上游高原台地：①遵义市播州区、绥阳县、湄潭县中北部、凤冈县西部等区域, 此区域中小型坝区分布相对较广, 地势较为平坦, 水田面积占比较大; ②正安县、凤冈县东部等中低山和盆坝地貌密集, 坡耕地、旱地占比较大; ③思南县西南部、石阡县北部、印江县西部等山地丘陵区域, 其余耕地大多分在上游高原台地、峰丛中低山区的山间洼地和峡谷河流阶地区。林地主要分布在流域上游毕节市纳雍县、织金县等乌蒙山区高中山、中山地貌区域以及下游沿河、德江、印江县等大娄山脉和武陵山脉交错区域。建设用地集中分布在流域上游毕节市、中游遵义市、省会贵阳市主城区及其周边, 总面积为1508.42 km², 占流域总面积的2.28%。

从变化趋势上看, 流域20年来土地利用格局发生了明显的变化, 除林地、建设用地以及水域面积呈现增加趋势外, 其余各地类面积均呈现下降趋势。以耕地面积下降最为明显, 在2000—2010年间, 面积减少了225.17 km², 而在2010—2020年间则快速下降, 总面积减少了995.36 km²; 林地2000—2020年面积增加328.47 km²; 城镇建设用地面积增加最为明显, 2020年面积为1508.42 km², 相较于2000年面积增加了1132.92 km², 主要是为满足对居住、工业和基础设施等用地的需求, 城镇建设用地增加的面积是2000年面积的4.02倍, 20年年均增长56.65 km², 建设用地增加的区域主要分布于贵阳市主城区周边、贵安新区、平坝区城区向安顺市区方向延伸以及毕节市和遵义市城区, 主要受城市化进程不断加快以及社会经济水平不断提升的影响。水域的面积20年间增加了224.68 km², 未利用地则减少了1.92 km²。

对乌江流域2000年、2010和2020年3期土地利用数据进行转移矩阵分析, 各地类之间的转换如表 3所示, 流域内由人类活动所引起土地利用类型转移方向主要表现出以下特点：

(1) 2000年—2020年来流域内各土地利用类型之间的相互转化主要在耕地、林地和建设用地之间, 以耕地的转换最为剧烈, 呈现面积在减少, 质量在上升特征。其中有2335.65 km²和1390.25 km²的耕地向林地及草地转移, 主要原因是2000年以来生态脆弱区退耕还林还草、农业结构调整、国土绿化以及水源地保护政策的实施, 使得不适宜和不能耕种的耕地流转变化, 进而导致耕地面积的减少; 其次是部分坡度较大以及土壤质量贫瘠的旱地由于常年撂荒导致杂草丛生从而发生变化; 此外, 城乡建设用地的快速扩张致使周边一部分耕地发生转移变化。

(2) 建设用地增加速度最快, 主要集中在地形平坦, 人口密集的贵阳市主城区、贵安新区、安顺市、毕节市和遵义市等城市周边, 特别是2010年以后城乡建设用地加速扩张, 总共增加1093.76 km², 是2000年到2010年增加的27.93倍, 20年来流域建设用地扩张呈现以省会贵阳市、毕节市和遵义市为节点的“V型”空间拓展结构特征, 主要原因是受城市用地扩张和社会经济发展的影响, 使得原本城市周边的耕地和草地被建设用地所取代; 而在建设用地转出过程中, 总计有68.63 km²的建设用地向耕地、林地和草地转移, 原因是部分地区由于易地扶贫搬迁、生态移民以及占补平衡等政策的影响, 搬迁后拆除旧房, 实施复垦或复绿, 整体上建设用地的变化主要以其它地类转入为主。

(3) 林地主要向耕地、草地和建设用地转移, 有2069.82 km²向耕地转移, 2160.63 km²转移为草地, 172.57 km²转移为建设用地, 主要以低覆盖的灌木林地、疏林地等低效林地转化为主。而草地和水域的变化幅度也相对较为明显, 其中有大量的草地向林地、耕地和建设用地转移; 水域则主要向草地、林地和耕地转移。而未利用地变化则相对稳定。

3.2 水源涵养与水质净化服务时空变化特征 3.2.1 水源涵养量时空变化

基于模型产水量模拟结果, 结合公式得到研究区2000—2020年水源涵养量, 如图 4所示, 从时间上看, 乌江流域水源涵养量数值差距较大, 平均水源涵养量从2000年的709.41 mm减少到2010年的532.80 mm, 减少了24.88%, 再增加到2020年781.32 mm, 从表 4中可以看出喀斯特分布区域水源涵养要明显低于非喀斯特分布区域, 2020年平均水源涵养量要少46.74 mm。2000—2020年水源涵养总量分别为472.54亿m³、354.70亿m³、520.14亿m³, 20年年均增长2.38亿m³, 总体上呈现先剧降、后增加的整体变化趋势。与降水的空间分布对比来看, 降水为流域生态水源涵养空间分布的主要影响因素, 2010年前后(2009—2011)由于贵州省发生百年一遇的特大旱灾, 降水量大幅度降低, 在此影响下, 流域产水量及水源涵养量也大幅度降低。

水源涵养量在空间上分布上具有明显的地区差异性, 高低值区域分布存在一定的规律性, 空间格局基本稳定, 三个年份均表现为西部和中部较低, 南部和东北部较高, 由西到东逐渐递增特征, 且存在共同的高低值区域。其中水源涵养量较高的区域主要分布在流域南部六盘水市普定县、水城县、六枝特区以及流域北部遵义市境内大娄山区域, 包括播州区、湄潭县、绥阳县、务川县一带; 水源涵养量中等的区域主要分布在流域中游黔西县到金沙县一带以及下游思南县等区域。而流域上游毕节市境内乌蒙山区域包括威宁县、赫章县、以及人口密集区周边水源涵养量较低, 主要原因一是该区域降水相对较少, 二是该区域是我国西南地区喀斯特地貌分布最为广泛的地区之一, 也是经济最为落后的区域, 人地矛盾突出, 早期森林砍伐、植被破坏等造成水土流失, 土壤退化等现象较为严重, 虽然近年来政府水资源建设、植树造林、退耕还林还草等工程的大范围实施使得植被覆盖状况改善, 但局部石漠化现象仍然存在, 水土流失还较突出, 区域大面积岩石裸露造成地表水源涵养能力差。

流域内水源涵养量随着高程的增加呈现先增加后减少的趋势, 量低值主要为海拔2100 m以上的区域, 高值区域主要分布在海拔1300 m到1500 m之间。坡度上, 流域内水源涵养量整体上随坡度的增加呈现逐渐增加的趋势, 增加趋势较为平缓。

3.2.2 水质净化时空变化

乌江流域2000年、2010年和2020年TN、TP输出强度如图 5、图 6所示。从空间上看, 流域内TN、TP高强度输出区均主要分布在流域中下游, 其中TN高强度输出区域主要分布在贵阳市南部乌江一级支流鱼梁河、翁安河; 遵义市境内湘江、湄江以及铜仁市境内六池河、石阡河下游所在子流域, TN输出强度最小值均在5.60 kg/hm²以上。TP高强度输出区域主要分布在贵阳市东部清水河、瓮安县瓮安河; 东北部湄潭县、凤冈县, 铜仁市思南县乌江干流以及芙蓉江、清溪河、三江交汇处等子流域区域。所在子流域TN输出强度最小值均在为0.98 kg/hm²以上。TN、TP低强度输出区域均主要分布在流域上游乌江北源毕节市境内赫章县、黔西县, 贵阳市南部独木河等区域。就整个流域TN、TP输出来看, 上游地区相比中下游区域水质净化服务功能较好。主要原因是上游耕地以旱地占绝对优势, 以旱作植物为主, 施肥少, 所以TN、TP较少, 且该区域坡度大、土壤贫瘠等因素导致N、P营养物质加速汇流河湖流向中下游, 而中下游水田在耕地中占比比上游大, 单位面积施肥多; 同时中游瓮安、开阳、息烽还是全省主要的磷矿开采区, 水体中的TP浓度比上游高。

从时间上看, 乌江流域2000、2010和2020年在TN、TP输出总量上存在差异, TN输出总量从2000年输出32883.97 t减少到2020年输出32258.08 t, 呈现逐渐减少的趋势, 表明流域TN净化量持续增强, 水质净化服务效果在持续变好。而TP净化服务效果上20年来有浮动趋势, 2000、2010和2020年TP输出总量分别为5376.29 t、5386.42 t、5323.35 t, 呈现先上升后下降的趋势, 总共减少了52.94 t。总体上看, 乌江流域20年来水质净化功能出现稳中变好趋势。

乌江流域TN、TP输出强度随着区域高程的升高呈现逐渐下降的趋势, 700 m高程以下区域TN、TP输出强度最大, 平均值分别达到602.17 kg/km², 108.15 kg/km², 低值区域为高程2100 m以上区域, TN、TP平均输出强度分别为280.59 kg/km², 43.20 kg/km²。小于5°的区域, TN、TP平均输出强度为616.78 kg/km², 100.92 kg/km², 大于40°以上区域输出强度呈现断崖式下降, TN、TP平均输出强度分别为2.81 kg/km², 0.41 kg/km²。

从乌江流域喀斯特与非喀斯区域TN、TP平均输出强度(表 5)来看, 非喀斯特区域的TN、TP平均输出量高于喀斯特区域。其中TN输出在2010年输出差距最大, 达到48.21 kg/km², 高出10.62%。TP输出在2010年输出差距最大, 达到12.31 kg/km², 高出16.71%。

3.3 土地利用变化对水源涵养及水质净化的影响 3.3.1 不同土地利用类型水源涵养变化

乌江流域2000年、2010年和2020年不同土地利用类型水源涵养量如图 7所示, 可以看出, 不同土地利用类型水源涵养能力存在明显差异, 其中草地和林地水源涵养能力最强, 平均水源涵养深度为809.53 mm、733.41 mm, 其次是耕地, 平均水源涵养深度为602.4 mm, 主要原因是草地和林地在根系深度上要高于耕地, 且同时植被冠层、枯落物部分能够有效截留降水, 且喀斯特发育区林下多裸露, 因而草地的滞留拦截坡面地表径流能力较刚开始恢复的次生林地\\人工林地更强, 而耕地由于长期耕作, 土壤容易板结, 减少了水分的入渗能力, 导致地表径流增加, 水源涵养能力降低。因此水源涵养服务能力草地>林地>耕地, 而水域缺少植被覆盖, 建设用地受人类活动影响强烈, 水源涵养能力相对较差。水源涵养总量上, 流域林地和耕地水源涵养总量最大, 平均值分别达到177.77×10⁹ m³和176.92×10⁹ m³; 其次是草地, 水源涵养总量为92.15×10⁹ m³; 水域和建设用地最少, 平均值不到2×10⁹ m³。整个流域地类主要以耕地、林地和草地为主, 是流域水源涵养服务功能的主要参与者和贡献者, 以2020年为例, 三者总面积占流域面积的96.87%, 贡献流域水源涵养总量的99.77%。

从时间上看, 不同土地利用类型水源涵养量变化呈现波动趋势。耕地、林地和草地是影响水源涵养能力变化的主要地类, 流域耕地转换为其它地类中, 转换为林地和草地的面积最大。相关研究表明, 退耕还林还草有利于减缓气候变暖背景下水源涵养量下降的速度, 增加区域水源涵养能力^[29]。而林地转换为耕地和草地过程中大多数发生在稀疏林地, 间接增加了水源涵养能力。

3.3.2 土地利用类型与TP、TN输出量的响应

将2000、2010和2020年流域不同地类面积与TN、TP输出总量进行相关性分析(表 6), 结果表明, 不同土地利用对TN和TP的输出量存在较大差异, 即对水质净化服务功能存在差异。从TN、TP输出总量与各时期地类面积的相关系数来看, 与耕地、林地和草地均具有较强的相关性, 而与水域、建设用地和未利用地的相关性较弱, 其中TP、TN输出总量与耕地的相关性最强, 2020年相关性系数均达到0.98, 与林地和草地的相关性次之, 2020年TP、TN输出总量与草地相关性系数均在0.84以上, 与未利用地较弱, 说明生态系统类型对流域水质净化的影响与不同土地利用类型的面积有关, 耕地是流域非点源污染负荷量最为密切的土地利用。

在子流域尺度上, 以相关性较强的耕地、草地和林地为例, 分析上述三种地类与TN、TP输出总量的关系, TN、TP输出总量与耕地的面积存在明显的相关性, 当耕地达到峰值时, 两者的输出量同样也存在峰值, 即子流域内耕地面积越大, 则TN、TP输出量越大; 当子流域内林地的面积相对大于耕地时, TN、TP输出总量也相应降低; 而草地则和耕地一样, 也存在较为明显的峰值变化, 表明不同土地利用类型下的TN、TP输出量存在耕地>林地>草地的排序, 水质净化服务功能则相反。

通过对各子流域内TN、TP输出总量与各地类的回归分析来看, 两者的输出总量与地类面积间存在着线性关系, 其回归方程如表 7所示。TN、TP输出总量与耕地的线性拟合的程度最高, R²分别为0.97和0.96;其次为林地, R²分别为0.82和0.80;与草地的R²分别为0.70和0.71。由此可见流域内耕地面积与TN、TP输出总量间呈现高度线性相关; 林地和草地面积次之; 与建设用地和水域面积相关性小; 未利用地面积则与TN、TP输出总量无明显相关。

综上来看, TP、TN输出总量受耕地面积的影响, 耕地面积越大, 则输出总量越多, 进一步说明耕地是导致水质净化服务退化的主要贡献源; 林地和草地具有一定的水质净化服务功能; 建设用地和未利用地则对水质净化无明显作用。

2000—2020年来流域内耕地对TN、TP输出贡献率有减少的趋势, 但仍然贡献了流域内TN、TP输出量的85%以上, 20年来耕地TN、TP输出贡献率分别下降了1.85%、2.63%。林地和草地TN贡献率分别在5%和7%以内, TP贡献率均在3%以内, 建设用地在2000—2010年TN、TP输出贡献率变化幅度不大, 2010—2020年变化幅度较大, 其中TN贡献率增加了2.82倍, TP增加了2.22倍, 主要原因是2010—2020年间建设用地面积增加了1093.76 km², 大量的其它类型地类转换为建设用地增加了地表径流系数, 导致水质净化服务变差, 因此TN、TP输出在2010—2020年间贡献率上升较快。未利用地、水域对整个流域的TN、TP贡献率基本上可忽略不计。综上所述, 由于林地、草地和流域水质密切相关, 因此林地和地草地面积越多, 流域整体的TN、TP贡献率就会越小, 从而水质净化能力越强。因此, 乌江流域水质净化服务功能提升主要基于耕地、林地和草地来考虑。

4 讨论 4.1 乌江流域水源涵养与水质净化影响因素分析

流域水源涵养与水质净化受到地形地貌、地表覆盖、区域气候以及人类活动的共同控制^[30]。研究基于InVEST模型以量化的方式空间化展现流域水源涵养与水质净化服务的本质。

水源涵养服务受到区域降水、蒸散和径流变化的影响, 是一个多途径复合的过程^[31]。降水丰富与否很大程度上影响区域水源涵养功能强弱, 而蒸散受到区域纬度、地表下垫面覆盖状况以及高程决定, 两者与区域水源涵养相关性较为明显。干旱少雨和地形起伏较大的区域水源涵养能力较差, 纬度高, 植被覆盖低的区域是整个流域水源涵养能力最差的区域。

土地利用通过改变下垫面的覆盖结构从而影响水源涵养过程, 是影响水源涵养的因素之一^[32]。流域上游2000—2400 m的高原区域主要以旱地为主。80年代不合理的土地利用方式、过度开发自然资源和主要以木材能源的生活生产方式, 迫使原生态结构遭到破坏, 水土流失、石漠化现象严重, 后期经修复后的植被主要以受人为干扰后封育形成的次生常绿落叶阔叶混交林、喀斯特灌丛、灌草丛、草坡等为主, 截留能力差, 土壤无法恢复到早期原生林时候状态, 甚至石漠化分布区域无土壤供植被生长, 大量水分流失, 水源涵养能力差; 中下游区域以500—1400 m的丘原、低山丘陵为主, 水田占比较多, 植被茂密, 土壤水分保持能力强, 径流流速低, 水源涵养能力强。

水质净化服务方面, 林地和草地的水质净化服务明显好于其它地类, 主要是氮磷营养物输出负荷系数、植被滞留效率较高, N、P吸纳、截留、分解及转化作用强, 输出强度小。耕地作为流域内TN、TP高强度输出土地利用类型, 水质净化服务功能最差, 耕地N、P大量施放是引起流域水质变坏的关键因素, 其中N肥能够迅速转换为硝态氮、P易与土壤颗粒以及其他化合物结合, 未被吸收的部分随径流汇入湖库, 因此退耕政策以及耕地转化为城镇建设用地的过程中能够改善水质净化效率。此外, 流域中游息烽河、洋水河、瓮安河等子流域为省内磷矿资源、磷化工企业分布密集地区, 磷矿开采废水、企业生产废水及渣场渗滤液对流域内地表水、地下水造成一定污染, 造成下游净化难度大, 也是影响流域水质净化差异性的重要因素。此外区域喀斯特地貌发育强烈, 地表切割破碎, 区域小气候和生态水文响应复杂, 难以准确研究界定和量化其变化过程和环境影响因子, 地表N、P等营养物在喀斯特地貌区由于土壤贫瘠, 常不经土壤层的过滤作用直接进入地下水系统, 经地下河汇入河流湖库, 一定程度上对流域水质净化服务产生影响。

4.2 基于水源涵养与水质净化的流域水环境保护

乌江流域不同河流、湖泊水质长期处于中营养状态, 近年来富营养等级有降低的趋势, 但局部区域仍居高不下, 如乌江干流沿江渡、大乌江镇、乌杨树等断面常年监测均总磷浓度超标。因此, 就如何改善流域水环境生态服务功能, 是当前管理部门的重要难点之一。

植被覆盖减少, 林冠层截留以及土壤持水量降低, 应由植被和土壤截留保持的水分和营养物质汇入河流造成水源涵养能力和水质变差, 因此, 改善水源涵养服务总体上要坚持流域内山水林田湖草生态要素的协同治理和一体化保护。目前主要的问题在于上游生态脆弱区水源涵养能力有待加强, 上游高地海拔区土层较薄, 主要以稀疏林地和旱地为主, 可以在此区域继续加强退耕还林还草措施, 重建恢复森林、草原、湿地等生态系统, 并在旱地区域开展水和土壤污染协同防治, 综合防治农业面源污染和生产生活用水污染等。此外, 应当重点保护和修复上游威宁国家草海湿地公园生态系统; 持续推进三岔河、六冲河以及白莆河等重要支流以及黔中水利枢纽工程的水土流失综合治理等生态工程建设, 减少水土流失现象, 防止泥沙和有机污染物汇入乌江主干流。

前文分析表明, 林地大面积分布有助于流域水质净化, 因此从防治和改善流域河流湖库水质的角度, 首先应该增加重点区域林地比例, 特别是流域上游重要河流湖泊周边实施林地种植, 提高林分质量, 提高植被覆盖; 其次可在河流湖库100—200 m缓冲范围内合理配置土地利用格局, 尽可能增加林地、草地覆盖, 湖库水源地消落带实施种植生物修复植物, 以达到增大植物的生物量, 最大程度地去除土壤中N、P污染物目的; 耕地作为流域内面积占比最大的土地利用类型, 占流域面积的44.90%, 说明农业生产活动仍然占主导地位, 是非点源污染的防控重点, 尤其是乌江上游六冲河、三岔河以及支流流域, 应尽量控制农田比例, 科学施肥, 减少其对水质净化服务的影响, 提高流域整体生态对水环境水质净化作用, 在坡耕地、低效耕地以及靠近饮用水源地周边等区域实施退耕还林还草工程, 加强低海拔平坝丘陵区的生态农业建设, 促进流域土地资源空间优化; 对当前农业产业结构进行调整, 实施国土绿化减少耕地对河流水质的污染, 如大力推进现代化生态高效农业, 提高耕地单位面积产出率, 降低对大面积耕地的依赖等。此外, 息烽河、洋水河、瓮安河等重点采矿区域要加强子流域整体管护, 实施主干和重要支流河长制, 及时发现河道水质异常、入河排污口排放异常和侵占河道等问题。虽然建设用地与整个流域水质净化服务相关性不明显, 但不可忽视其对污染物的径流加速汇流作用, 城市化目前仍然是流域主要人类活动, 应当严格控制建设用地平面空间扩展, 以建设用地的垂向扩展代替平面蔓延, 增加建设用地内部绿化、公园以及以乔、灌木为主的绿地网络, 增加内部可透水性面积, 缓解流域的水质净化压力。

5 结论

基于InVEST模型, 本文分析了乌江流域20年来水源涵养与水质净化服务时空格局特征, 并探究了流域土地利用变化对水源涵养与水质净化服务的影响机理, 得出以下结论：

(1) 乌江流域各年份水源涵养服务在空间上呈现出一致性, 均表现为由西到东呈现逐渐递增特征, 高值区域集中分布在流域中下游低中山丘陵和宽谷盆地区, 低值区域主要分布在上游高原山地生态脆弱区。各土地利用类型水源涵养方面, 草地、林地是流域水源涵养的主要贡献者, 两者水源涵养量占流域总量的78.97%, 其次为耕地、水域、建设用地、未利用地。2000—2020年间水源涵养量呈现先减少后增加的变化趋势, 平均水源涵养量在532.80—781.32 mm之间。

(2) 流域水质净化服务表现为TN输出整体上呈现减少的趋势, TP输出总量变化甚微。表明流域N的净化能力呈现逐渐增强趋势。TN、TP高强度输出区域均主要分布在流域中下游地区, 各土地利用类型中, 林地和草地的分布对水质净化起正向促进作用, 而耕地对TN、TP输出占绝对主导地位, 贡献了流域内85%以上的营养物质。

(3) 土地利用结构多元化、破碎化一定程度上影响流域水源涵养与水质净化服务功能, 流域内草地、林地与水源涵养和水质净化呈明显的正相关关系, 耕地与水质净化呈负相关关系, 退耕还林还草有利于改善流域水源涵养与水质净化服务。