亚热带典型农林流域面源磷负荷估算与源解析

李情，唐代生，孟岑，李裕元，吴金水

（1.中南林业科技大学林学院，长沙 410004；2.中国科学院亚热带农业生态研究所亚热带农业生态过程重点实验室，长沙410125；3.中国科学院亚热带农业生态研究所长沙农业环境观测研究站，长沙 410125）

磷作为水生生态系统富营养化的限制性营养元素，其主要来源包括点源和面源，相较于点源，面源磷污染的控制和管理仍存在巨大的挑战[1]，主要因为其具有空间分布的广泛性和不确定性等特点[2]。面源磷负荷的精准核算，关键源区和空间分布规律的明晰，是实现基于流域面源污染高效控制的重要前提。

面源磷污染的随机性、迁移和转化的复杂性和滞后性决定了其负荷估算和控制难度大，因此通过模型模拟是量化污染负荷、识别关键源区的重要手段[3]。近年来我国面源磷负荷估算模型在构建与应用方面取得了长足的进展。根据参数特征，模型可分为集总式和分布式。分布式模型如AGNPS、APEX、HSPF、SWAT 等应用广泛，分布式模型考虑了污染物迁移过程，涉及多介质、多尺度的污染物产汇模拟和农业管理等措施评估[4]。这些模型多建立在大量基础资料和物理机制的基础上，对于特定区域环境能够实现精确模拟，但却限制了其在区域外的应用和推广[5]。通用土壤流失模型作为国内外广泛应用的土壤侵蚀模型，具有结构简洁、参数物理意义明确、计算简单等优势[6]，但其仅能用于吸附态磷的模拟估算。近年来有研究将其与输出系数等集总式模型进行耦合模拟流域总磷负荷量及空间分布[7-8]，但多数研究未对模型模拟效率进行验证，且模型缺乏参数本地化校准。

我国长江中下游亚热带丘陵区是我国重要的农业粮食产区，农业集约化程度高，水系分布密集，面源磷污染给当地生态环境带来了严重的负面影响[9]。同时该区域具有相似的景观格局[10]，相比中大尺度区域，选取典型农林流域为对象构建磷负荷估算模型具有较好的可操作性和推广性。本研究选取长沙县金井河流域为研究区域，基于改进的输出系数和RULSE 模型构建流域磷负荷模型，对面源磷负荷进行模拟估算，解析不同形态磷负荷的空间分布特征，以期为亚热带农业流域面源磷负荷的核算及控制管理提供有效手段和科学依据。

1 材料与方法

1.1 研究区域概况

金井河流域位于湖南省长沙县金井镇（27°55′~28°40′N、112°56′~113°30′E），湘江一级支流捞刀河上游区域，流域总面积为134.92 km2。流域内以低山丘陵为主，中北部地势高，南部地势低，海拔在56~440 m 之间，平均坡度为18°，坡度大于25°的面积约占流域面积的1/3。该地为典型南方红壤丘陵地貌，土壤主要类型为红壤和水稻土，其气候属于典型亚热带湿润季风气候，年平均降水量为1 200~1 500 mm。流域水利设施以水库为主，其中有1 座Ⅰ中型水库（金井水库，位于流域东南部），2 座小Ⅰ型水库，15 座小Ⅱ型水库，2 422 口山塘，总蓄水容量为2 150.38 万m3。土地利用方式以林地、农田、居民地为主，其面积占比分别为65.1%、30.2%和1.4%。林地主要分布于丘陵岗地等海拔坡度较高的区域，主要树种有马尾松（Pinus massoniana）、杉木（Cunninghamia lanceolata）和油茶（Camellia oleifera）等。农田分布在南部河流冲击平原及丘陵沟谷之间，以水稻田为主，平均施磷量为 6 500 kg·km-2·a-1，并辅以有机肥。流域内无中大型城镇和工厂等大型点源污染。金井河流域研究区域如图1所示。

1.2 模型框架构建

1.2.1 模型数据库

流域内设有1个标准气象观测站和3个小型气象观测站。流域出口设有长期定位水文监测点（2014—2020年），选取2014—2019年观测的溶解态磷负荷用于模型参数优化。利用ArcGIS 软件的水文分析模块和D8 流向法将金井河流域划分为20 个集水区，选取其中部分集水区建立水文水质监测点（图1），并将2020 年观测得到的溶解态磷负荷与总磷负荷用于模型效果验证。本研究所需其他基础资料包括遥感影像数据、DEM 数据（分辨率5 m×5 m）、土地利用数据、土壤数据、降雨数据等的来源见表1。

表1 金井河流域面源磷污染负荷模型数据来源Table 1 Data sources of nonpoint source phosphorus pollution load model in Jinjing watershed

图1 金井河流域研究范围Figure 1 Study area of Jinjing watershed

1.2.2 模型构建

面源磷负荷可分为溶解态磷负荷和吸附态磷负荷。由于土壤侵蚀是驱动吸附态磷输出的主要因素，因此基于RUSLE 实现模型对吸附态磷的模拟。由于磷在河道迁移过程中存在衰减，因此将传统输出系数模型与一阶河道衰减系数耦合构建溶解态磷模拟模块。模型框架示意见图2。

图2 金井河流域面源磷负荷模型框架Figure 2 Modeling framework of nonpoint source phosphorus load model in Jinjing watershed

1.3 溶解态磷负荷模型构建

本研究基于不同土地利用类型（林地、农田、居民地）输出系数及一阶河道衰减系数构建溶解态磷负荷模型，公式如下：

式中：LSRP为溶解态磷负荷，kg·a-1；ei、mi分别为第i种土地利用溶解态磷输出系数和面积，kg·km-2·a-1和km2；k为溶解态磷河流衰减系数，d−1；t为溶解态磷从集水区河道到流域出水口的输移时间，由河道长度与平均流速计算得到，d；Lj为第j个支流或河道上游入口的溶解态磷负荷，kg·a-1；tj为第j个支流或上游流入的溶解态磷从对应河道起点到集水区出口的输移时间，d。

此外对于活塞流远大于分散流的中小流域，一般假设流域面源污染物沿集水区内的给定河段进入，平均损失超过河流长度的1/2（相当于t/2）[11]。本研究仅对瞬时流速进行监测，因此基于已有研究中瞬时流速与河道平均流速的经验公式计算得到平均流速[12]，计算公式为：

式中：d为流程，m；v为平均流速，m·s-1；Q为集水区出口河流检测瞬时流量，m3·s-1；Q1为河流年均流量，m3·s-1；s为集水区面积，m2；s1为无量纲的积水区面积；Q2为无量纲的流量；g为重力加速度，9.8 m·s-2。

采用贝叶斯-蒙特卡洛方法（Markov Chain Monte Carlo，MCMC）对公式（1）中不同土地利用磷输出系数和河道磷衰减系数后验分布进行校准计算[13]。具体方法如下：假设Y为响应变量，即流域溶解态磷负荷的自然对数。设Yh（h=1，…，l）服从正态分布，其均值为μh（对应每一次采样流域日溶解态磷负荷Yh）、方差为σ2（对应所有Yh）：

溶解态磷瞬时负荷的对数服从正态分布。均值μh是一个空间过程，通过模型得到；方差σ2是一个常数，采用观测数据得到。

通过前期不同土地利用溶解态磷负荷和磷的河道衰减系数获取模型参数先验分布，设模型参数为相互独立的先验分布：

由于模型方差没有先验信息，故对所有模型参数使用模糊先验分布，即对输出系数（ei）先验分布取自然对数，对河道衰减系数（k）先验分布开方。为模型参数指定正态分布，并为参数指定逆gamma（IG）[14]：

采用 Gibbs sampling 算法实现 MCMC 计算[15]，利用已知数据（LSRP、mi、t）计算获得溶解态磷负荷模型中e1（林地）、e2（农田）、e3（居民地）和k的后验分布，并进行校正后输出模型参数结果[13，16]。

1.4 吸附态磷负荷模型构建

吸附态磷负荷主要来源于土壤侵蚀流失，通用土壤流失方程（USLE）和RUSLE 是估算土壤流失的最经典模型，本研究以RUSLE 与污染物负荷方程构建吸附态磷负荷模型[17]，计算公式如下：

式中：LXFP为吸附态磷负荷，kg·a-1；λ为泥沙输移比；XP为土壤侵蚀模数，t·km-2·a-1；CP为磷污染物背景含量，g·kg-1；η为磷污染物富集比；mi为第i种土地利用面积，km2；R为降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1；K为土壤可蚀性因子，t·hm2·h·hm-2·MJ-1·mm-1；LS为坡长坡度因子，无量纲；C为植被覆盖因子，无量纲；P为水土保持措施因子，无量纲。

1.4.1 泥沙输移比（λ）

泥沙输移比反映泥沙从土壤侵蚀到进入河道某拦截面处的减少程度，定义为流域出口某断面输沙量与断面以上土壤侵蚀产沙量之比，与流域面积、地形、土壤性质及降雨径流密切相关。计算方法包括直接计算和模型估算，参照长江流域和洞庭湖流域已有的研究结果确定本研究流域内的泥沙输移比为0.4[18-19]。

1.4.2 磷污染物背景含量（CP）

磷污染物背景含量指土壤磷污染物含量，即土壤中总磷的含量。根据流域内土壤采样分析，得到流域内土壤磷含量均值为0.51 g·kg-1。

1.4.3 磷污染物富集比（η）

磷污染物富集比是土壤侵蚀产生的泥沙磷含量与雨前土壤表层磷含量之比。以往研究发现泥沙中磷含量常高于土壤中的磷含量，富集比在1~4 之间，且多在2 附近变动，故本研究的磷污染物富集比定为2[20-22]。

1.4.4 降雨侵蚀力因子（R）

降雨侵蚀力因子是土壤侵蚀的动力因子，表征降雨引起的潜在动能，与雨量、雨强、降雨历时有关，本研究选用WISCHEMIER等[23]的研究中的公式计算：

式中：P为多年平均降雨量，mm；Pi为第i月多年平均降雨量，mm。通过流域内4个气象观测站2014—2020 年的降雨数据得到平均月降雨量和年降雨量，利用Arc⁃GIS软件采用泰森多边形进行空间插值绘制R因子空间分布（图3）。

1.4.5 土壤可蚀性因子（K）

土壤可蚀性因子表示降雨侵蚀力作用下土壤被冲离的难易程度，用土壤侵蚀速率来衡量各种土壤的侵蚀敏感度，该因子与土壤性质密切相关。本研究选用WILLIAMS等[24]提出的EPIC模型计算：

式中：SAN为土壤砂粒含量，%；SIL为土壤粉粒含量，%；CLA为土壤黏粒含量，%；C为土壤有机碳含量，%；SOM为土壤有机质含量，%；SN为土壤粉粒和黏粒含量总和，%。通过采样得到本研究区域土壤砂粒、粉粒、黏粒、有机质含量，并利用ArcGIS 软件绘制K因子空间分布（图3）。

1.4.6 坡长坡度因子（LS）

坡长坡度因子可根据地形地貌直接反映坡长因子与坡度因子对土壤侵蚀的影响。本研究利用Arc⁃GIS水文分析功能，根据坡长坡度公式计算得到[25]：

式中：λ是由 DEM 提取出的水平投影坡长，m；θ是由DEM 提取出的坡度，（°）。本研究区域LS因子空间分布如图3所示。

1.4.7 植被覆盖因子（C）

植被覆盖因子用于反映植被覆盖度和管理水平对土壤流失量的作用关系。本研究采用蔡崇法等[26]提出的模型方法计算：

式中：NDVI为归一化的植被指数，通过遥感影像波谱计算得到。利用ArcGIS软件绘制的C因子空间分布见图3。

图3 金井河流域R、K、LS和C因子空间分布Figure 3 Spatial distribution of R，K，LS and C factor in Jinjing watershed

1.4.8 水土保持措施因子（P）

水土保持措施因子是采用水土保持措施下与直接顺坡耕种下土壤流失量之比。本研究采用赋值法，综合参照已有研究成果[27-28]，将林地、农田、居民地、水域的P因子分别赋值为1、0.35、0.01、0。

2 结果与讨论

2.1 流域磷负荷验证

基于改进的输出系数模型和RUSLE 方程分别模拟得到流域不同土地利用溶解态磷输出负荷以及空间尺度上的吸附态磷负荷，并进一步拟合得到流域面源磷输出负荷。为了对模型的模拟效果进行验证，分别将总磷、溶解态磷、吸附态磷负荷模型模拟值与2020 年观测值进行对比，结果表明总磷模拟值与观测值的误差绝对值为3.3%~27.1%（表2）。溶解态磷负荷决定系数（R2）和纳什效率系数（NSE）分别为0.84和0.82，误差绝对值为3.9%~24.7%。颗粒态磷负荷R2和NSE分别为 0.68 和 0.65，误差绝对 值为 8.9%~30.0%（表2，图4）。模型模拟结果具有一定精度，可用于亚热带典型农林流域面源磷污染的空间估算，且溶解态磷模拟效果优于吸附态磷。

图4 溶解态磷和吸附态磷负荷观测值与模拟值比较Figure 4 Observed and the modeled of dissolved and adsorbed phosphorus load

表2 金井河流域面源磷负荷估算与验证Table 2 Estimation and validation of nonpoint source phosphorus load in Jinjing watershed

2.2 流域磷负荷空间分布

模拟结果表明，研究区溶解态磷和吸附态磷负荷的空间分布存在较大的差异（图5），其中溶解态磷负荷平均为29.6 kg·km-2·a-1（12.5~37.4 kg·km-2·a-1），高值区多分布在中部和西南部的农田和居民地聚集区，而东南部和北部区域相对较低。吸附态磷负荷平均为34.7 kg·km-2·a-1（20.5~59.4 kg·km-2·a-1），高值区集中在北部多林地区域。研究表明不同形态面源磷负荷污染来源和输移机制存在差异，溶解态磷主要来源于化肥、生活污水和畜禽粪污，而吸附态磷主要源自高流量水文过程下的土壤侵蚀。

图5 金井河流域溶解态磷和吸附态磷负荷集水区空间分布Figure 5 Sub-watershed spatial distribution of dissolved and adsorbed phosphorus load in Jinjing watershed

2.3 流域磷负荷源解析

流域面源磷负荷年输出总量为7 310.3 kg·a-1，其中溶解态磷和吸附态磷负荷分别占39.0%和61.0%（图6）。不同土地利用输出负荷总量为林地（5 004.6 kg·a-1）>农田（2 170.8 kg·a-1）>居民地（134.9 kg·a-1）。林地溶解态磷和吸附态磷负荷总量分别为1 186.3 kg·a-1和 3 818.3 kg·a-1，农田溶解态磷和吸附态磷负荷总量分别为1 544.2、626.6 kg·a-1，居民地溶解态磷和吸附态磷负荷总量分别为157.8、17.2 kg·a-1。

图6 金井河流域不同土地利用类型面源磷输出情况Figure 6 Nonpoint source phosphorus exports of different land use types in Jinjing watershed

不同土地利用面源磷输出负荷分别为林地56.5 kg·km-2·a-1、农 田 56.9 kg·km-2·a-1、居民 地 35.6 kg·km-2·a-1。林地吸附态磷负荷（43.1 kg·km-2·a-1）显著高于溶解态磷负荷（13.4 kg·km-2·a-1，P<0.01）。农田溶解态磷负荷（40.5 kg·km-2·a-1）显著高于吸附态磷负荷（16.5 kg·km-2·a-1）。居民地溶解态磷负荷（31.1 kg·km-2·a-1）显著高于吸附态磷负荷（4.5 kg·km-2·a-1）。研究流域作为典型的亚热带丘陵区的典型农林流域，其农业生产活动强度大、施肥量高是造成农田溶解态磷负荷高的主要原因。此外居民地普遍缺乏污水收集和处理设施，地表硬化率较高，溶解态磷易随地表径流进入河道。林地溶解态磷负荷虽然较低，但其吸附态磷负荷高，主要由于林地多在流域内的山岗高地，坡度大、土壤侵蚀量大。

3 结论

（1）基于相关评价指标所构建的面源磷负荷估算模型模拟具有一定精度。基于改进的输出系数模型构建溶解态磷负荷模型，并采用贝叶斯方法求解模型参数能够有效解决模型参数本地化问题；基于RULSE 模型构建的吸附态磷负荷模型能够较为准确地模拟亚热带丘陵区土壤侵蚀状况，但受降雨监测频率及模型参数限制，时间精度有待进一步提升。

（2）金井河流域面源总磷负荷平均为64.3 kg·km-2·a-1，其中溶解态磷负荷为29.6 kg·km-2·a-1，吸附态磷负荷为34.7 kg·km-2·a-1。流域溶解态磷负荷高值区集中在西南部平原区，吸附态磷负荷高值区主要分布在西北部林坡地区域。

（3）农田、居民地溶解态磷负荷显著高于林地，平均输出负荷分别为 40.5、31.1、13.4 kg·km-2·a-1；林地吸附态磷负荷显著高于农田和居民地，负荷均值分别为 43.1、16.5、4.5 kg·km-2·a-1。亚热带丘陵区面源磷负荷的治理需要根据不同形态磷负荷输出空间分布差异有针对性地开展管控措施。