基于最小累积阻力模型的洱海流域上游面源污染“源—汇”风险格局研究

董亚坤王 钰曾维军

(1.云南农业大学 水利学院,云南 昆明650201;2.自然资源部 云南山间盆地土地利用野外科学观测研究站,云南 昆明650201;3.云南省土地资源利用与保护工程实验室,云南 昆明650201)

洱海是云南9大高原湖泊之一[1]。洱海流域上游为高原山间盆地,使得该地区集中了整个流域约57%的耕地[2],同时洱海流域上游也是洱海最为主要的水源补给区[3]。面源污染具有分布广,污染源多,治理难,成本高等特点,成为区域生态环境污染的主要贡献者[4]。因而控制洱海流域上游面源污染是改善洱海水质的决定因素。

在水质受到污染前对面源污染进行识别和防控具有更好的治理效果[5]。目前涉及洱海流域面源污染的相关研究多围绕水质的理化性质及变化特点进行探讨[6-8],而近年来利用“源—汇”景观理论,以构建阻力面来识别和防控面源污染成为热门方向[9-10]。“源—汇”风险是“源—汇”景观理论的延伸,最初是用来研究大气污染中CO2的来源与吸收,后来被引入到景观生态学当中,发展为“源—汇”景观理论[11]。“源—汇”方法将流域景观赋予“源”“汇”属性,并依据源汇景观指数评价流域内景观对面源污染的贡献,识别流域内面源污染的风险等级[12]。面源污染风险等级识别有多种方法。Jiang,Zhou和Xin[10,13-14]分别采用景观空间负荷比指数、网格景观空间负荷比指数来划分面源污染风险等级,而Zhu、王金亮和刘帆[5,15-16]借助最小累积阻力模型构建阻力面方式分别划分了三峡库区、王家沟、重庆北碚区的面源污染“源—汇”风险等级。相较于SWAT,Ann AGNPS等模型在大中尺度流域的应用,最小累积阻力模型偏定性研究,对数据要求较低,精度较高[17],更适合洱海流域上游这种面积较小(1 230.00 km2),数据易获取的小尺度流域。面源污染研究的方法很多,但有的对研究区域的基础资料要求高[18-19],有的缺乏预测的精确性或适用性不强[20-21]。同时水文水质气象等数据具有较高的保密性,对于非相关研究领域人员获取较难。遥感技术具有全天时、全天候、易获取等优点[22],方便了土地利用数据的获取。通过构建氮磷等面源污染物的评价体系,以阻力成本的方式来识别面源污染的“源—汇”风险等级,是研究洱海流域上游面源污染的有益尝试。因此,本文以洱海流域上游为例,选取影响面源污染转移的阻力因子构建评价体系,使用最小累积阻力模型形成阻力面,识别2005,2010,2015和2020年洱海流域上游面源污染“源—汇”风险等级,分析面源污染风险等级及其转移变化。通过面源污染“源—汇”风险等级的划分,快速识别面源污染发生的高风险区域并开展重点治理,为保护洱海流域生态环境提供理论与决策支持。

1 材料与方法

1.1 研究区概况

洱海流域位于云南省大理州境内,处于三江交汇地带[1]。洱海流域上游属中亚热带西南季风气候带,一年当中只有干季和湿季,无春夏秋冬四季之分,气候宜人,降雨量丰沛,但空间分布不均,95%的雨量集中在雨季的5—10月[7]。研究区独特的地形地貌,优越的环境与气候,丰富的水资源,使得该区域人口密集,耕地面积大而集中,经济活跃,土地利用率高,受面源污染严重[8]。

1.2 基础数据及处理

本文基础数据包括2005,2010,2015和2020年的土地利用类型数据、归一化植被指数、人口密度以及DEM数据、土壤数据、地形湿润指数。其中归一化植被指数和2005,2010,2015和2020年土地利用类型数据来自中国科学院资源环境科学与数据中心(http:∥www.resdc.cn/data.aspx?DATAID=184),分辨率均为30 m(见图1);DEM数据来自地理空间数据云(http:∥www.gscloud.cn),分辨率为30 m;土壤数据来自中科院南京土壤所(http:∥www.issas.ac.cn/),分辨率为30 m;人口密度来自WorldPop(https:∥www.worldpop.org/project/categories?id=18),分辨率为1 km。

图1 洱海流域上游2005—2020年土地利用类型分布

以上数据经处理全部采用CGCS 2000坐标系,分辨率统一为30 m的栅格数据,数据处理以及后续的空间分析使用软件ArcGIS 10.5。

1.3 研究方法

1.3.1 最小累积阻力模型 最小累积阻力模型(the minimal cumulative resistance model,MCRM)是耗费距离模型的扩展应用,其核心观点是景观要素通过克服阻力来实现景观类型的相互转化[23]。最初用来反映物种从“源”到目的地所耗费的最小代价[24],后被广泛应用于生态安全[17]、生态廊道或风险路径[25]、景观格局[26]等。该模型考虑源、空间距离和阻力基面3方面因素[23],计算公式为:

式中:MCR为最小累积阻力值;Dij为某物种从“源”景观j至“汇”景观i的空间距离;Ri为修正后的阻力系数;f为正相关关系函数,表征源景观至空间中某一点路径的相对易达性。

1.3.2 “源”地识别 面源污染中的“源”景观是指能够促进面源污染发生、发展的景观类型,污染风险高,而“汇”景观则与之相反[27]。产生面源污染物的土地利用类型加快面源污染的发展,被认为是“源”景观[28]。依据中国土地分类体系将洱海流域上游的土地划分为耕地、林地、草地、水域、建设用地和其他用地,其中耕地施用农药化肥产成氮磷等面源污染物,为“源”景观;建设用地通常被认定为点源污染,但研究区以农村为主体,生活污水的排放设施在近几年由政府投资修建才得以缓解,因此建设用地被认为是“源”景观(图1)。

1.3.3 阻力基面评价指标体系构建 选取影响面源污染的土地利用类型、相对高程、坡度、地形湿润指数、人口密度、土壤可蚀性和植被指数7个因子,构建阻力基面评价指标体系。

不同因子在面源污染的转移过程中所起的阻力作用不同,权重也不同,参考相关文献和Delphi专家打分法为各因子权重赋值[10,15,29]。而同一因子的不同等级也有不同的阻碍作用,依据文献或自然断点法利用重分类工具划分5个等级,分别赋予对应因子阻力系数值为1,3,5,7,9[30](表1)。其中,相对高程和归一化植被指数,值越高,面源污染发生的风险越小,面源污染风险等级越低;坡度、地形湿润指数[31]、土壤可蚀性[32](K值)和人口密度,值越大,以及土地利用类型中的耕地和建设用地越多,面源风险等级越高。相对高程、坡度和归一化植被指数与面源污染风险等级呈负相关性,主要原因是相对高程越高、坡度越大的区域受到人类活动影响越小,氮磷等污染物排放越低;而归一化植被指数值越大,植被覆盖度越高,生态环境越好,其能有效保持水土,减少氮磷等污染物的流失。地形湿润指数、土壤可蚀性、人口密度和土地利用类型与面源污染风险等级呈正相关性,地形湿润指数可对潜在土壤水分含量和径流产生潜在能力进行量化,值越高,污染物越容易随土壤侵蚀而流失,土壤可蚀性反映不同土壤受到侵蚀的速度,值越高,土壤受到的侵蚀越高,污染物流失也越多;人口密度和土地利用类型是人们生产、生活与经济的反映,也是污染物排放最为集中之地,人口越多,密度越大,耕地和建设用地就越多,面源污染风险等级也随之越高。

表1 洱海流域上游面源污染阻力基面评价指标体系

1.3.4 阻力面构建 阻力面反映了各种“流”(物质和能量等)从“源”克服各种阻力到达目的地的难易程度,也客观表现了事物空间运动的趋势和潜在可能性[24]。基于阻力基面评价体系,运用ArcGIS 10.5中的栅格计算器进行叠加计算,形成2005,2010,2015和2020年的综合阻力基面,再将“源”景观和综合阻力基面运用cost-distance工具进行叠加,生成对应阻力面,即最小累积阻力面。

1.3.5 面源污染风险等级划分 本文将影响流域面源污染的阻力面作为反映面源污染“源—汇”风险等级的指标,该指标反映面源污染物从“源”景观到最终汇集处之间的可达性[23]。阻力面值越大,面源污染“源—汇”风险等级越小,景观单元所起到的汇作用越强,面源污染越不容易发生,否则反之。众多文献利用数据的突变点采用自然断点法进行分类。研究区属于盆地,采用自然断点法会把山地错分到极高风险区当中,综合各种数据分析后采用几何间断方法划分面源污染风险等级,将研究区分为极低风险区、低风险区、中风险区、高风险区和极高风险区。

2 结果与分析

2.1 各因子的阻力系数分布特征

阻力因子的阻力系数分布从1—9等级,表示阻碍污染物的阻力在逐渐增大。由图2分析可知,相对高程、土壤可蚀性和人口密度空间分布规律相近,均是从中间向边界逐渐增大,各等级层次分明,高度集中,呈现圈层分布特征,而坡度、地形湿润指数和归一化植被指数三者空间分布杂乱无序,各阻力等级相互交错,呈小区域聚集,整体分散的特征。各阻力因子的阻力系数面积也存在着巨大差异。相对高程和坡度的阻力等级分布以1,3等级为主,面积是其余阻力因子同等级的数倍。土壤可蚀性的阻力等级分布5等级面积(604.13 km2)最多,远远超过剩余阻力因子最大值(382.92 km2)。人口密度、地形湿润指数和归一化植被指数的阻力等级分布以7,9等级最多,二者面积超过800.00 km2,远超其余阻力因子的两者面积(最高约400.00 km2)之和,与其他阻力因子形成鲜明对比。除了坡度外,各因子阻力等级分布均为中间低边界高,表明平缓区是面源污染物最集中的区域,易受到面源污染。

图2 研究区面源污染各阻力因子等级空间分布

2.2 综合阻力基面和阻力面的空间分布

阻力值的高低表明景观单元受到阻力因子的影响程度,阻力值越大,受到面源污染的概率越小,否则反之。由图3可以看出,2005,2010,2015和2020年的综合阻力基面的最高值均是8.48,而最低值在1.60。结合图1可知,阻力值的分布和土地利用类型的分布较为一致,但局部空间分布也有变化的趋势。主要在于2005—2020年,土地利用类型分布总体保持稳定,但景观内部耕地、草地面积呈现减少,林地、水域面积呈现增加的趋势。同时,茈碧湖镇、三营镇交界处以及凤羽镇的低阻力值空间分布明显减少,尤其在2005—2010年显著减少。从图4可知,4期的阻力面整体分布是基本一致的,都是低阻力值占据绝对主体地位。低阻力值占据平缓区域及附近范围,高阻力值分布流域边缘的极高山地,阻力值在波动中总体呈上升态势,由2005年的55 183.90增加到2020年的62 669.30。16 a来,阻力基面与阻力面的空间分布存在着明显的区域差异,但均呈现中间低边界高的趋势。而阻力面的阻力值增加,表明从2005—2020年,各种保护洱海的政策和工程起到良好作用,产生面源污染物的区域进一步缩小,阻碍范围进一步扩大。

图3 洱海流域上游2005—2020年面源污染综合阻力基面空间分布

图4 洱海流域上游2005—2020年面源污染阻力面空间分布

2.3 面源污染风险等级分析

由图5,表2可知,2005,2010,2015和2020年洱海流域上游的面源污染风险等级分布由里向外依次是极高风险区、高风险区、中风险区、低风险区和极低风险区,中风险区、高风险区和极高风险区占总面积绝大部分。同时受地形影响,研究区南部面源污染风险要高于北部,耕地及其邻近范围内的区域风险等级最高,表明极高风险区即为面源污染的关键源区。耕地“源”作用强,林地、草地和水域“汇”作用弱,因此,“源”景观的面积和分布是影响风险等级的关键因素。16 a间(图6),极高风险区从2005—2015年呈现逐步下降的趋势,到2020年又有所上升;高风险区在2005—2020年持续上升;中风险区在2005—2010年略微下降后,在2010—2020年快速上升;低风险区在2005—2015年快速上升后,到2020年又急剧下降;极低风险区在2005—2010年上升,之后持续下降。

图5 洱海流域上游2005—2020年面源污染风险等级分布

图6 研究区2005—2020年面源污染各风险等级面积趋势变化

表2 洱海流域上游各乡镇2005—2020年面源污染风险等级面积 km2

总之,极低风险区面积减少16.16 km2,低风险区面积减少7.01 km2,极高风险区面积减少35.74 km2,中风险区增加了38.69 km2,高风险区面积增加了20.23 km2。16 a间,茈碧湖镇极高风险区虽有波动,但面积基本上未发生变化,保持在65.00 km2左右;邓川镇和上关镇全镇面积的80%以上被中等风险及以上所占据,而极低风险区几乎为零,受到的面源污染威胁很高;三营镇和右所镇均处在平缓区域,且极高风险区面积最多,都大于70.00 km2,2020年相比2005年极高风险区面积也均有所下降;凤羽镇和牛街乡的各面源污染风险等级较为相似,但牛街乡各面源污染风险等级分布较为分散,而凤羽镇较为连续。总之,各乡镇都受到较高的面源威胁。

2.4 面源污染风险等级转移分析

2005—2020年,洱海流域上游共发生了15种面源污染风险等级转移变化(图7,表3),主要是在极高风险区、高风险区、中风险区之间转移,以转出为主,转出面积分别为80.93,72.51,75.63 km2。其中极低风险区和低风险区的转入转出主要在盆地边缘,而其他面源污染风险等级的转入转出主要在盆地中间。低海拔向高海拔过渡的区域面源污染风险等级变化最为剧烈。

表3 研究区2005—2020年面源污染风险等级转移矩阵 km2

图7 研究区2005—2020年面源污染风险等级转移分布

低风险区、中风险区、高风险区具有相似的转移特征,与其相邻的等级转移最多,最高分别为41.12,33.88和43.29 km2,相距越远彼此转移越少。极高风险区和极低风险区均只相邻一个面源污染风险等级,故只沿一侧单向递减转移。极低风险区几乎只和低风险区相互转移,而极高风险区因距离递减作用,向高风险区和中风险区转移最多,分别为51.10和29.47 km2,向其余面源污染风险等级几乎不转移。在空间分布上,茈碧湖镇、凤羽镇和右所镇的交界处以及牛街乡的北部是主要转移的区域,这些区域是“源—汇”景观交界处及“汇”景观所在地,也是面源污染治理重点关注的区域。风险等级转移的方向和程度受到多种因素的影响,而人类活动和洱海流域生态保护政策的实施是最主要的原因。

3 讨论

本文采用四期数据研究面源污染风险等级变化,相较于一期数据而言,能够体现时间和空间变化特征,表明面源污染“源—汇”风险识别方法具有可行性[5,9]。虽然划分出5个级别,但数据段是不同的。这种划分方法是否可以进行不同年份之间的比较有待探讨。但可以确定的是基于最小累积阻力模型,可以初步识别面源污染的风险等级,找到面源污染的关键区域[12,23]。同时阻力因子的数量以及权重会受人为因素影响,不同的学者、不同的研究区域会选择不同的阻力因子,具有一定主观性。此外,耕地分类只考虑到一级地类,未进一步细分旱地和水田,而旱地的氮磷流失要高于水田[27],也未分析氮磷污染物从“源”到“汇”过程中的机理。上述局限性需要进行更深入研究。

洱海流域上游的高风险区和极高风险区主要分布在耕地及其附近,这些区域与经济、人口、城镇建设等具有正相关性。限制城镇规模,管控耕地,进行景观优化[33]、农业产业结构调整[4]是减少氮磷等污染物的有效方式,可从“源”头上有效降低面源污染的产生和转移。对耕地质量较低、易造成水土流失的坡耕地要逐步退出[34],对茈碧湖、西湖和海西湖所在的区域重点管控,实施缓湖滨带生态修复。此外,极低风险区、低风险区具有波动性,和中高风险区呈现此消彼长,在风险等级转变的重点区域加大植树造林力度,增加林地和草地的面积,提高“汇”景观的优势度及主导作用。

农药、化肥的禁用,大蒜等高水高肥的经济作物禁止种植,大面积的退耕还林还湿等众多保护洱海的政策实施,改变了洱海流域上游的面源污染风险等级,使得面源污染物大幅度减少,洱海水质得到较大提升。上述研究结果与《中国环境状况公报》《洱海保护治理与流域生态建设“十三五”规划纲要》《云南省环境状况公报》等文件总体保持一致,如《中国环境状况公报》《云南省环境状况公报》指出:洱海氮磷含量显著降低,水质由Ⅲ类接近Ⅱ类。而王琦[35]、徐建锋[36]、孔佩儒[37]等人基于最小累积阻力模型分别识别了查干湖、丹江口水库和海河流域的面源污染风险格局,指出耕地是高风险区域,与本文研究结果相吻合。同时农户对于环境保护的积极性较高,甚至可以损失自己的部分利益,而一系列保护洱海的政策的实施使得当地农民经济损失严重,导致大量人员外出务工,众多农田被荒。因此,以绿色高效低成本方式,识别和控制洱海流域上游面源污染,降低因洱海保护造成的对农户及农业的损失,实现洱海保护和农民利益的双赢,会对洱海的保护起到更加非凡的意义。

4 结论

(1)阻力因子和源的分布使阻力基面和阻力面的空间分布、阻力值的大小存在区域性差异。随着海拔的升高,阻力值的分布呈现中间低边界高,“源”景观的作用被“汇”景观逐渐取代。2020年比2005年阻力面值提高了7 485.40,表明洱海流域上游的生态环境有所改善。

(2)研究区4期极高风险区面积全都大于360.00 km2,所占比例达30%,而中风险区以上占比超60%,整体面源污染风险等级偏高;16 a间,极低风险区、低风险区和极高风险区面积分别减少了16.16,7.01和35.74 km2,中风险区和高风险区面积分别增加了38.69和20.23 km2。

(3)面源污染风险等级由“源”向外依次分布:极高风险区、高风险区、中风险区、低风险区和极低风险区,即距离“源”越近,风险等级越高;受空间距离和“源”的影响,风险等级表现为中部高于边界,南部高于北部。

(4)16 a间的中、高、极高风险区面积转出最多,分别为80.93,72.51和75.63 km2,分别占其总面积的36.44%,53.68,19.89%。茈碧湖镇、凤羽镇和右所镇三镇交界处以及牛街乡的北部是面源污染风险等级转移的主要区域,主要是中高风险等级所在地;极高风险区是洱海流域上游面源污染的关键区域,而减少氮磷的流失需要围绕关键区域展开。