基于空间马尔可夫模型的流域非点源优先管理区识别

Identification of Non-Point Source Priority Management Areas in Watershed Based on Spatial Markov Model

摘要: 非点源污染是导致地表和地下水环境恶化的重要原因之一。优先管理区(Priority Management Areas, PMAs)识别对非点源污染精准高效管控至关重要。然而，在已有PMAs识别研究中很少考虑河流系统中的传输损失，且缺乏对识别结果准确性和合理性的验证。基于此，本文通过耦合SWAT模型与空间Markov算法，构建了一个考虑河流滞留效应和削减负荷的多级PMAs识别框架，应用于海河水系的王快水库流域，并验证PMAs识别结果在不同管理措施下的准确性。主要结果如下：(1) 在校准和验证期SWAT模型性能良好，纳什效率系数NSE > 0.60、R² > 0.79；(2) 河道氮磷滞留系数与径流量呈非线性关系，存在显著的径流阈值。当径流量介于1.14~8.23 m³/s，总氮(TN)和总磷(TP)的滞留系数趋于稳定。(3) TN的PMAs均为#29和#30号子流域，负荷贡献为74.35% (枯水年)、58.45% (平水年)和76.46% (丰水年)，面积占流域总面积的14.22%。(4) TP的PMAs为#33 (枯水年)、#28和#33 (平水年)、#29和#30 (丰水年)号子流域，负荷贡献分别达到53.19%、78.83%和69.69%，面积占1.20%、2.95%和14.22%。(5) 仿真结果显示，PMAs削减效率约为非PMAs的4.25倍(TN)和5.37倍(TP)，证实了PMAs识别结果的合理性。该研究结果可为相关流域的非点源污染控制提供有效方法支撑。

Abstract: Non-point source (NPS) pollution is one of the important contributors to the deterioration of the surface and groundwater environment. The identification of priority management areas (PMAs) is crucial for effective management of NPS pollution. However, transport losses within river ecosystems are seldom considered when identifying PMAs, and there is a lack of validation for the accuracy of the results. Therefore, a framework for multi-level PMAs identification was conducted, considering the effect of river retention and retention load by integrating the SWAT model with the Markov algorithm. This framework is applied to the Wangkuai Reservoir Watershed in the Haihe River system, and the accuracy of the PMAs identification results under different best management measures (BMPs) is verified. The main results are illustrated as follows: (1) The SWAT model showed strong performance with high statistical indicators (NSE > 0.60, R² > 0.79) for the calibration and validation periods. (2) A non-linear relationship between the retention coefficient and runoff with a significant runoff threshold was detected. For both total nitrogen (TN) and total phosphorus (TP), the retention coefficients held steady when runoff surpassed 1.14~8.23 m³/s. (3) The PMAs for TN are sub-watersheds #29 and #30, with load contributions of 74.35% (dry year), 58.45% (normal year), and 76.46% (wet year), and the area accounts for 14.22% of the total basin area. (4) The PMAs for TP are sub-watershed #33 (dry year), sub-watershed #28 and sub-watershed #33 (normal year), and sub-watershed #29 and #30 (wet year), with load contributions reaching 53.19%, 78.83%, and 69.69% respectively, and the areas account for 1.20%, 2.95%, and 14.22% respectively. (5) The reduction efficiency of PMAs is approximately 4.25 times (TN) and 5.37 times (TP) that of non-PMAs, confirming the accuracy of the PMAs identification results. The results will provide a reference for NPS pollution control in similar watersheds.

1. 引言

随着人口增长和经济发展，水污染已被认为是全球范围内的一个重要环境问题。通过有效的立法和技术解决方案，如限制工业污水流出和污水处理厂，点源污染得到了明显的抑制[1]。因此，非点源(Non-point Pollution Sources, NPS)污染成为了地表和地下水环境的主要污染来源。为了减轻水污染对水生态系统的不利影响，应制定有效的缓解方法来评估NPS污染的风险。NPS污染由降雨产生的地表径流或地下径流将地表污染物迁移至河道造成，存在明显的不确定性和时空异质性[2]。因此，识别流域优先管理区(Priority Management Areas, PMAs)是高效污染管理的关键。

在流域系统中，PMAs被定义为相对较小的部区域，这部分区域贡献了绝大部分的污染，并将大多数污染物输送到附近水体中[3] [4]。与传统关键源区相比，PMAs中包含了负荷强度、水质浓度和污染物的潜在削减程度，为改善水环境提供了更有效的途径[5]。污染物进入河道后，从上游向下游传输污染物的效率会减弱。这主要是由于污染物降解、吸附和沉积的转化过程的影响。因此，河流生态系统中的运输损失评估对于识别PMAs极为重要。流域内的污染物负荷很容易监测，但对下游的贡献无法直接获取。因此，将污染贡献作为PMAs识别的一个标准更能反映实际情况。研究表明，空间马尔可夫模型是模拟上游污染物向评估点迁移的有效工具，将流域上下游河网关系转化为矩阵，量化子流域污染贡献[6]。

PMAs的识别是开展最佳管理措施(Best Management Practices, BMPs)和应对生态环境问题的基础[7]。BMPs是为了减少流域内的NPS污染、改善水质和保护生态环境等采取的一系列措施。BMPs分为源控措施(如化肥减施和耕作管理)和拦截措施(例如植被缓冲带和植草河道)。由于源控措施主要是通过改变农业管理方式达到从源头控制NPS污染的目的，具有成本低，环境友好等特点，被广泛应用[8]。拦截措施属于过程中的措施，通过采用一系列物理、生物或工程手段，直接拦截、过滤或吸附污染物，极大降低污染物向水体、土壤等环境介质的迁移扩散。以往的研究通过对不同BMPs进行生态效益评价，筛选出最适合治理流域NPS污染的BMPs，实现NPS污染的高效治理。然而，当前PMAs的识别研究主要为空间定位，缺乏对布设不同BMPs下流域PMAs的识别有效性进行验证，无法从BMPs削减负荷的角度对流域PMAs进行生态效益评价[9]。

在NPS污染的研究中，通过河道监测收集数据是最直接的方法。这可能成本高昂且耗时，并且无法有效地进行情景分析，而建模为这个问题提供了一个可行的解决方案。以往的研究通常采用经验模型(如输出系数模型、磷(P)指数模型等)或机制模型(如SWAT模型、AnnAGNPS模型、GWLF模型等) [10]。近年来，考虑空间异质性的分布式水文模型已成为水文研究的热点之一，并且更广泛地用于分析流域的水文过程[11]。SWAT模型是生态、水文和环境研究中使用最广泛的水文模型之一。与其他模型相比，SWAT模型需要手动调节参数，耗时较长，但广泛应用于模拟水文和物质的迁移转化。

本文以中国河北省保定市王快水库流域为研究区域。研究目标为：(1) 构建SWAT模型并验证其适用性；(2) 估算TN、TP营养负荷；(3) 探讨河道对TN和TP污染物的削减能力；(4) 识别PMAs，并讨论污染物削减量与PMAs面积之间的关系。(5) 验证不同BMPs下PMAs识别结果的准确性。

2. 研究区域

王快水库位于中国河北省保定市曲阳县和阜平县交界处，东经114.15˚，北纬38.45˚，如图1。王快水库流域面积约3770 km²，全年气候变化明显。年均降水量为626.4 mm，年平均相对湿度为52%，年平均蒸发量1200~1300 mm，年平均径流量为7.53亿m³，年平均输沙量约295万t。流域地势西北高东南低，海拔181~2781 m，主要土壤类型为褐土、棕壤和潮土，土地利用类型以草地、林地和耕地为主。

3. 数据与方法

3.1. SWAT模型

SWAT模型是由美国农业部-农业研究局(USDA-ARS)开发的半分布式物理模型，包含水文、泥沙和营养盐等多个模块[12]。模型主要用于预测流域尺度上不同土壤类型和土地利用方式对水、沉积物和营养盐的长期影响[13]。基于地形、土地利用和土壤等要素，利用ArcGIS软件将流域细分为37个子流域，并基于土壤类型、坡度坡向和土地利用类型等条件将子流域进一步划分为400个水文响应单元(Hydrological Response Units, HRUs)。HRUs是分布式水文模型中表征流域空间异质性的基本单元，精准量化降水截流、径流生成等过程，广泛应用于水土资源管理与气候变化水文效应评估。

Figure 1. The DEM of the Wangkuai Reservoir Watershed

图1. 王快水库流域数字高程图

3.2. 数据来源

SWAT模型的输入数据包括DEM、地形、土壤属性、土地利用和气象数据。数据来源见表1。

Table 1. Data sources of the Wangkuai reservoir watershed

表1. 王快水库流域数据来源

数据	内容	来源
数字高程	海拔、坡度、坡向等	中国地理空间数据云网站 (http://www.gscloud.cn)
土壤数据	土壤类型、物化属性	世界土壤数据库 (https://www.ncdc.ac.cn)
土地利用	土地利用类型	国家地球系统科学数据共享平台 (http://www.geodata.cn)
气象数据	降水、气温、风速、相对湿度、太阳辐射等	资源环境科学与数据中心 (https://www.resdc.cn)
径流和水质	径流量(逐日)、总氮浓度(逐月)、总磷浓度(逐月)	中国河北省保定市阜平水文站

3.3. 校准与验证

SWAT模型通常需要引入多个参数来描述流域的特征。选择合适的参数对于SWAT模型的成功应用至关重要，并且可以通过校准和验证来增强模型性能。本文结合全局敏感性分析和one-at-a-time (OAT)分析方法来识别更敏感的参数。首先，选取显著等级较高的参数，设置阈值和模拟次数。其次，将模拟数据与观测数据进行对比。最后，利用SWAT-CUP的顺序不确定性拟合算法(SUFI2)对月径流量、TN和TP数据进行校准和验证。SWAT模型的性能通过百分比偏差(PBIAS)、决定系数(R²)和纳什效率(NSE)系数进行评估[14]。PBIAS表示观测值和模拟值之间的差异。R²表示模拟数据和观测数据之间的拟合程度。NSE系数用于验证水文模型模拟结果的准确性。

$PBIAS=\frac{{\displaystyle \sum _{i=1}^n\left(P_i-O_i\right)\ast 100}}{{\displaystyle \sum _{i=1}^n{O}_i}}$ (1)

$R^2=\frac{{\left[{\displaystyle \sum _{i=1}^n\left(O_i-\overline{O}\right)\left(P_i-\overline{P}\right)}\right]}^2}{{\displaystyle \sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}{\displaystyle \sum _{i=1}^n{\left(P_i-\overline{P}\right)}^2}}$ (2)

$NSE=1-\frac{{\displaystyle \sum _{i=1}^n{\left(P_i-O_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}}$ (3)

式中，$P_i$ 是被评估的第i个模拟值，$O_i$ 是第i个观测值，n是观测总数。根据现有文献，当NSE > 0.5且R² > 0.6时，模拟结果良好[15]。当径流和污染物的PBIAS分别在±25%和±70%内，模拟结果可信[16]。

3.4. 方法

马尔可夫转移矩阵是用来描述系统在不同状态之间转移概率的矩阵。对于一个具有n个状态的马尔可夫链，转移矩阵为$n\times n$ 的方阵，矩阵中元素表示为在当前处于状态i的情况下，下一步转移到状态j的概率。而空间马尔可夫模型是在传统马尔可夫模型的基础上考虑空间位置信息，转化为描述空间上的结构状态，将矩阵中的元素与子流域的空间结构关系一一对应，表征为流域空间单元的上下游关系，状态i和j分别被描述为子流域的上下游之间的关系，即污染物在流域i状态下，下一步流向流域j的概率。Grimvall和Stalnacke基于马尔可夫链的理论将子流域的空间结构定义为上下游之间有效关系的转移矩阵，通过流域上下游之间的滞留系数描述任意物质在子流域系统中的输移，使数据与流域空间结构有效结合，对污染负荷进行空间分配[17]。Chen等人基于流域模型、流模型和马尔可夫链模型，描述流域过程的动态变化以及各种水质状况。结果表明河道迁移过程对于河网中的水质恶化至关重要[6]。因此，本文基于SWAT模型和空间马尔可夫算法，结合河网矩阵和滞留系数矩阵，以流域河道出口处为评估点，计算子流域对污染物在污染迁移过程的削减量并进行分配，确定PMAs。最后，以不同BMPs下的削减率验证PMAs识别结果。

3.4.1. 确定河网结构

根据SWAT模型中王快水库流域的子流域划分情况，按照河网水系的分布情况和矢量关系，将河网结构转化为矩阵。对于n个子流域$S_1,S_2,\cdots ,S_n$ ，河网结构用$n\times n$ 的矩阵H表示，如图2。

$H\left(i,j\right)=\{\begin{array}{l}1,如果S_j是S_i下游最近的点\hfill \\ 0,否则\hfill \end{array}$ (4)

式中，矩阵H为马尔可夫链的转移矩阵，表示n个不相交的子流域上下游之间的水流连通状况，所有概率均为0或1。其中，1表示水作为污染物的载体存在直接的水流连接，从通过上游$S_i$ 河道流向下游$S_j$ 。0表示$S_i$ 与$S_j$ 之间没有直接的联系，不存在水流连通，阻断了污染物的直接传输路径。

Figure 2. River network and matrix formulation

图2. 河网与矩阵计算

3.4.2. 污染物总负荷

SWAT模型将子流域划分为若干个由性质均一的土地利用类型、土壤和坡度构成的HRUs，并对每个HRUs的污染物负荷量进行计算。每个子流域的负荷就是HRUs的污染物负荷量累加的结果。为了计算各子流域的污染物负荷，还需要根据流域内污染物的不同赋存形态，明确污染物负荷的构成。因此，本文的氮、磷负荷构成分别如下：

$C_{TN}=O_{ORGN}+N_{NSURQ}+N_{NLATQ}+N_{N{O}_{3}L}+N_{N{O}_{3}GW}$ (5)

$C_{TP}=O_{ORGP}+S_{SEDP}+S_{SOLP}+P_{PGW}$ (6)

式中，$C_{TN}$ 和$C_{TP}$ 分别为氮元素和磷元素的组成。$O_{ORGN}$ 代表有机氮负荷。$N_{NSURQ}$ 表示地表径流硝态氮负荷。$N_{NLATQ}$ 表示横向流动中的硝酸盐负荷。$N_{N{O}_{3}L}$ 和$N_{N{O}_{3}GW}$ 分别是土壤表面的硝酸盐负荷和地下水中的硝酸盐负荷。$O_{ORGP}$ 表示沉积物有机磷负荷，$S_{SEDP}$ 表示泥沙产磷量。$S_{SOLP}$ 和$P_{PGW}$ 分别表示可溶性磷负荷和地下水中磷负荷。

$L_i={\displaystyle \sum _{j=1}^n\left(A_{HRU}\left(i,j\right),C_{HRU}\left(i,j\right)\right)}$ (7)

式中，$L_i$ 表示子流域i排放的污染物总负荷量(kg)，A和C分别为第i个子流域，第j个HRUs的面积(ha)和负荷强度(kg/ha)。矩阵L为各子流域释放的污染物负荷矩阵。

$L={\left(L_1,L_2,\cdots ,L_i,\cdots ,L_n\right)}^{\text{T}}$ (8)

3.4.3. 污染物滞留系数

污染物在随水迁移的过程中会受到污染物降解、吸附、沉积、和阻滞剂等因素的影响，使上游向下游传输污染物的效率降低。因此，河道对污染物的滞留在污染物迁移过程中具有重要影响[6]。滞留系数可视为对污染物的削减能力，作为衡量河流运输效率的关键指标。对于流域内的每条河流，污染物的输入和输出的差异可以用滞留系数来表示。

$R_i=\frac{L_{in}-L_{out}}{L_{in}}$ (9)

式中，$R_i$ 为子流域i的滞留系数。$L_{in}$ 表示上游负荷与子流域i的负荷总和。$L_{out}$ 表示子流域i出口处的负荷。$R_i$ 取值范围为[0, 1]，解释了滞留在河道中的污染物负荷的转移概率，即$R_i$ 越大，说明该河道段对污染物的截留、吸附或沉淀作用越强，从而降低了污染物向下游的传输量。在流域尺度上，

$L_{inTN}=O_{inORGN}+N_{inNO}+N_{inNH4}+N_{inNO2}$ (10)

$L_{outTN}=O_{outORGN}+N_{outN{O}_3}+N_{outN{H}_4}+N_{outN{O}_2}$ (11)

$L_{inTP}=O_{inORGP}+M_{inMINP}$ (12)

$L_{outTP}=O_{outORGP}+M_{outMINP}$ (13)

式中，$L_{inTN}$ 和$L_{outTN}$ 分别为河流入口处和出口处氮的组成。$L_{inTP}$ 和$L_{outTP}$ 分别表示河流入口处和出口处氮的生成量。$O_{ORGN}$ 代表有机氮负荷。$N_{N{O}_3}$ 表示硝酸盐，$N_{N{H}_4}$ 表示铵态氮。$N_{N{O}_2}$ 代表亚硝酸盐氮。$O_{ORGP}$ 表示有机磷负荷。$M_{MINP}$ 表示矿物磷。空间单元的滞留系数用$n\times n$ 的矩阵R来表示，

$R=\left[\begin{array}{cccc}{R}_1& & & \\ & R_2& & \\ & & \ddots & \\ & & & R_n\end{array}\right]$ (14)

3.4.4. 污染物负荷贡献

步骤1：确定评估点。评估点被定义为河网中的典型位置，例如流域出口、水质监测点等。本文以河道出口32号子流域和37号子流域作为评估点，评估上游子流域的氮磷在河道迁移过程的转移贡献。

步骤2：计算污染贡献率。流域河网空间关系矩阵的连通性和河道滞留系数矩阵的削减损失性，能明确污染传输的路径和负荷量。在污染负荷沿着河道向下游传输的过程中，依据流经河道段对应的滞留系数，对分配的污染负荷进行相应调整。因此，引入马尔可夫矩阵计算，结合河网空间关系矩阵和滞留系数矩阵，量化上游对下游的污染贡献[6] [17]。

$L^{\prime }={\left({\tilde{H}}_k\right)}^{n_k}{V}_k\ast L$ (15)

${\tilde{H}}_k\left(i,j\right)=\{\begin{array}{l}\tilde{H}\left(i,j\right),i\ne k\\ 1,i=j=k\\ 0,\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}i=k,j\ne k\end{array}$ (16)

$V_k\left(i\right)=\{\begin{array}{l}1,i=k\hfill \\ 0,否则\hfill \end{array}$ (17)

式中，$L^{\prime }$ 表示污染物负荷贡献矩阵。*表示元素乘法。k表示河道评估点的子流域。$\tilde{H}$ 由河网空间关系矩阵H和贡献矩阵1-R的乘积定义，表示子流域产生的污染负荷在河网中上游子流域i向下游子流域j的迁移传输。$\tilde{H}$ 到${\tilde{H}}_k$ 的转换表示子流域$S_1,S_2,\cdots ,S_n$ 传输到评估点k的有效迁移传输。$n_k$ 表示迁移传输过程中河网最长河流的长度。$V_k$ 是一个$n\times 1$ 矩阵，用于提取评估点所在列的值。

3.4.5. 评估点削减量

当流域的地表水达到生态环境质量的水质要求，从而有效控制水污染。根据中国的《地表水环境质量标准》(GB3838-2002)中II类水质标准，TN的标准浓度小于0.5 mg/L，TP的标准浓度不大于0.1 mg/L。将评估点的超标负荷作为各子流域的削减量。

${\alpha }_k=\{\begin{array}{l}\frac{C_k-C_{limit}}{C_k},C_k>C_{limit}\hfill \\ 0,C_k\le C_{limit}\hfill \end{array}$ (18)

$C_k=\frac{L_k}{F_k}$ (19)

式中，${\alpha }_k$ 表示评估点k处的污染物超标比例。$C_k$ 表示评估点k处的实际污染物浓度(mg/L)。$C_{limit}$ 表示水质标准浓度限值。根据$L_k$ 表示评估点k上游的污染物负荷。$F_k$ 表示评估点k处的流量(m³/s)。

$E_{LA}{}_k=L_{oadk}\ast {\alpha }_k$ (20)

式中，$E_{LA}{}_k$ 表示污染物的超标负荷量(kg)。$L_{oadk}$ 表示评估点k处的实际污染物负荷量(kg)。

3.4.6. PMAs识别

PMAs的识别标准是将评估点超出水质浓度标准限值的部分作为削减量，按照各个子流域贡献的比例将其分配到子流域，计算每个子流域的拟削减量。

$D=E_{LA}{}_k\ast \frac{L^{\prime }}{{\displaystyle \sum L^{\prime }}}$ (21)

式中，D表示子流域的削减量。ArcGIS软件中自然断点法是通过统计优化自动划分数据区间的分类方法，最大化类间差异并最小化类内差异。因此，利用自然断点法对流域污染削减量进行分类，更能揭示数据自然聚类特征的空间分级。

3.4.7. BMPs下的削减率

BMPs的评估验证结果通过措施实施前后的削减率来表示，

$N=\frac{L_{pre}-L_{now}}{L_{pre}}$ (22)

$N^{\prime }=\frac{{\displaystyle \sum N}}{a}\times 100\%$ (23)

式中，N表示实施措施的污染负荷削减比例，$L_{pre}$ 表示措施实施前的负荷量，$L_{pre}$ 表示措施实施后的负荷量，$N^{\prime }$ 表示分级PMAs的削减比率，a表示各分类级别的子流域数量。

BMPs的实施设置集中在源控措施和拦截措施。其中，源控措施主要集中在化肥减施情景，施肥量的减少与污染负荷的削减率成正相关，布设在流域的农业用地[18]。拦截措施主要集中在植被缓冲带和植草河道情景下，布设在所有土地利用类型。当田间与缓冲带的面积比超过一定范围后，大部分污染物在植被过滤带的前端已经被过滤截留，之后增加面积比并不能明显提高污染物的削减率[19]。随着河道植草河道长度的增加，对氮磷负荷的削减率不断增加，但河道植草到达一定长度后，削减率不再显著增加[8]。因此，结合王快水库流域的农业及自然地理环境等因素，本文共设置四大类BMPs，包括10种情景设置，布设在PMAs和非PMAs。(1) 基准情景：不添加任何BMPs。(2) 3个化肥减施情景：化肥减施10% (FR10)、20% (FR20)和30% (FR30)。(3) 田间和缓冲带面积比为1 (FS1)、5 (FS5)和10 (FS10)的3个植被缓冲带情景。(4) 长度为100 m (GW1)、200 m (GW2)和500 m (GW5)的3个植草河道情景。具体情景设置与参数设置见表2。

Table 2. BMPs scenarios setting

表2. 最佳管理措施情景设置

序号	BMPs	措施设置	代码	情景设置	输入文件	参数设置
1	基础情景	/	BAS	基准情景	/	/
2	源控措施	化肥减施	FR10	化肥减施10%	.mgt	FRT_KG减少10%
3			FR20	化肥减施20%	.mgt	FRT_KG减少20%
4			FR30	化肥减施30%	.mgt	FRT_KG减少30%
5	拦截措施	植被缓冲带	FS1	田间与缓冲带面积比为1	.ops	FSRATIO = 1
6			FS5	田间与缓冲带面积比为5	.ops	FSRATIO = 5
7			FS10	田间与缓冲带面积比为10	.ops	FSRATIO = 10
8		植草河道	GW1	植草河道长度为100 m	.ops	GWATL = 0.1
9			GW2	植草河道长度为200 m	.ops	GWATL = 0.2
10			GW5	植草河道长度为500 m	.ops	GWATL = 0.5

4. 结果与讨论

4.1. 模型校准与验证

本文SWAT模型的运行时间为2019~2023年。径流量、TN和TP的校准期为2019-2021年，验证期为2022~2023年。

Figure 3. The calibration and validation results of monthly

图3. 逐月校准与验证结果

在SWAT模型的校准与验证期间，NSE > 0.60、R² > 0.79，表明SWAT模型均处于可接受水平，如图3。校准与验证结果与观测数据吻合较好。其中，径流量、TN和TP的PBIAS最大分别为32%、24%和33%。总体而言，尽管模拟结果存在一些差异，但统计评估指标均在可接受范围内。因此，根据评估标准可以认为模拟结果是可信的，SWAT模型可以成功应用于王快水库流域[15] [16]。

4.2. 子流域氮磷负荷时空分布

根据年径流量大小，将2020年、2022年和2023年分别确定为典型平水年、枯水年和丰水年。王快水库流域TN、TP的时空分布，如图4所示。在不同水文年情景模拟下，子流域TN和TP的负荷输出量分别为0.008~3393.932 t和0.000~29.165 t。丰水年受降雨、径流量影响，TN和TP负荷输出偏高，而枯水年和平水年较低，表明降雨和径流在TN和TP负荷流失和迁移过程中发挥重要的作用。降雨量和径流量的大小成正比，流域污染负荷总体呈现南多北少的分布格局，原因可能是南部地区存在大量耕地，人类活动密集，农业集约化程度高。流域土地利用类型以耕地、林地为主，其中耕地主要集中在东南部，约23052.60 ha。耕地区集中于流域河道附近，分布在土壤表层的营养盐物质进入收纳水体，在降雨、径流的驱动下更容易造成TN和TP的迁移，造成流域污染输出量增加。同时，受流域水文条件和自然环境的影响，TN和TP负荷在空间呈现TN多，TP少的分布格局[20]。在不同水文情景下，TN负荷输出量约是TP负荷的117~13925倍，重度污染的TN负荷输出量约是TP的300~8513倍，TN重度污染面积约是TP的2~5倍，如表3所示。这表明TN污染具有显著的空间覆盖和较高的贡献率。

Table 3. Statistical table of TN and TP output

表3. TN和TP输出统计表

指标	水文年	重度污染子流域	负荷量 (kg)	负荷占比 (%)	面积(ha)	面积占比(%)
TN	枯水年	12、13、14、18、19、24、29、35	1344728.17	53.42	132226.99	40.11
	平水年	24、28、29、35、37	1524983.03	37.23	104661.96	31.75
	丰水年	2、6、12、13、14、18、19、24、29、30、35、37	22382646.90	68.87	202704.79	61.49
TP	枯水年	28、37	157.95	87.38	26993.03	8.19
	平水年	28、37	200.48	66.04	26993.03	8.19
	丰水年	24、29、37	74412.18	26.76	81081.63	24.60

Figure 4. Spatial distribution of pollution load

图4. 污染时空分布

4.3. 河道污染物削减能力

各子流域枯水年TN和TP的削减能力分别高达0.51和0.54。在平水年，各子流域TN和TP的最大削减能力分别为0.52和0.54。丰水年削减能力分别为0.51和0.52，相对于枯水年和平水年，丰水年削减能力相对较弱。滞留系数和径流量之间存在非线性关系，具有显著的径流量阈值。当径流滞留系数介于1.14~8.23 m³/s时，TN和TP的削减能力趋于稳定，如图5。由于水环境中转运过程和赋存形式的差异，TP滞留系数的径流阈值一般高于TN。这些具有较强污染物削减能力的子流域一般位于流域河道附近，在控制河道内的污染物方面发挥着至关重要的作用。与靠近汇水区的子流域相比，距离汇水区较远的子流域的负荷转移效率较低[7]。

Figure 5. The relationship between retention capacity and runoff in sub-watershed

图5. 子流域削减能力与径流量关系

4.4. PMAs识别

利用ArcGIS中的自然断点法对各子流域不同水文年情景下的TN和TP削减量进行了划分。得到了三个等级的分类结果。将一级作为PMAs的识别结果，如图6。在三个水文年情景下，TN的PMAs均为#29和#30号子流域，负荷贡献为74.35%、58.45%和76.46%，面积占流域总面积的14.22%。在枯水年，TP的PMAs为#33号子流域，平水年为#28和#33号子流域，丰水年为#29和#30号子流域，负荷贡献分别达到53.19%、78.83%和69.69%，面积分别占流域总面积的1.20%、2.95%和14.22%。

在三个水文情景模拟下，TN的PMAs无显著变化，这是因为TN的原有负荷量较大，并显著超出该流域水环境容量，但TN削减量与降雨和径流作用呈正比关系，削减量介于0.07~8933976.10 kg，在丰水年达到最大。TP的PMAs呈现显著的变化趋势，削减量介于0.00~101666.96 kg之间，在丰水年达到最大值。这表明，丰水年的负荷潜在损失量高于平水年和枯水年，降雨径流对PMAs的识别产生直接影响。

PMAs主要分布在流域东南区域。此区域为流域出口，污染负荷贡献随河流自上而下累积的污染物负荷量增大而增大，但流域的自净能力有限，无法及时、有效地去除污染物的影响，导致下游处的污染负荷量多于上游处。相反，王快水库流域的中部和北部区域没有明显的PMAs变化，因为PMAs主要关注的是实际到达评估点的负荷，而不是污染物的源输出[6] [10]。

为了更直观地解释PMAs来自流域的较小部分，根据污染物削减量对子流域进行降序排序，计算累积削减量和面积。TN和TP的削减率临界区间58.45%~76.46%和53.19%~78.83%之间，而面积为1.20%~14.22%区间。如图7。因此，随后的控制措施应集中在PMAs地区，才能有效控制污染源的扩散[3]。

4.5. 基于BMPs效率的PMAs合理性验证

在源控措施情景(FR10、FR20和FR30) BMPs下，王快水库流域PMAs的TN削减率分别达到44.50%~48.80% (枯水年)，33.06%~33.26% (平水年)和51.59%~51.80% (丰水年)，对TP分别削减了

Figure 6. The identification results of PMAs

图6. PMAs的识别结果

Figure 7. The relationship between cumulative reductions and cumulative areas

图7. 累积削减量和累积面积的关系

52.08%~57.24% (枯水年)，32.10%~38.08% (平水年)和54.96%~55.86% (丰水年)。在BMPs下，减少化肥施用量对TP的削减效率要优于对TN的削减效率，且丰水年削减效率最高。与先前的研究相比，本文布控化肥减施在PMAs处，能更有效削减流域TN和TP负荷[21]，如图8。

在拦截情景措施(FS1、FS5、FS10和GW1、GW2、GW5)的BMPs布控下，王快水库流域PMAs的TN削减率为68.64%~80.95% (枯水年)，64.85%~76.98% (平水年)和48.05%~52.78% (丰水年)，对TP分别削减了44.74%~74.73% (枯水年)，34.80%~49.78% (平水年)和54.20%~57.59% (丰水年)。如图8。与源控措施比，拦截措施的削减效果更优，但波动幅度较大[9]。这可能是因为植被缓冲带和植草河道受降雨和径流影响较大，更宽的植被缓冲带可以更有效的减少污染物流入水体。植草河道通过植被减缓径流，但过多又会破坏河道结构，导致泥沙堆积。

Figure 8. Reduction rates of PMAs and non-PMAs under different BMPs

图8. 不同BMPs下PMAs与非PMAs的削减率

因此，PMAs的定位能有效削减流域TN和TP负荷的输出。与TN的削减相比，流域BMPs对TP的削减率更高。与非PMAs相比，BMPs对PMAs的削减率更高，最高是非PMAs的4.25倍(TN)和5.37倍(TP)，验证了PMAs识别结果的准确性。

5. 结论

本文基于SWAT模型与马尔可夫算法的结合，建立了考虑河流滞留效应和削减负荷的王快水库流域多级PMAs识别框架，并验证了识别PMAs在不同BMPs下的准确性。主要结论如下：

(1) SWAT模型在王快水库流域的校准期和验证期的性能均被接受。

(2) 滞留系数与径流量呈非线性关系，径流量阈值显著。TN和TP在径流量介于1.14~8.23 m³/s时系数趋于稳定。

(3) TN的PMAs均为#29和#30号子流域，负荷贡献为74.35% (枯水年)、58.45% (平水年)和76.46% (丰水年)，面积占流域总面积的14.22%。

(4) TP的PMAs为#33 (枯水年)、#28和#33 (平水年)、#29和#30 (丰水年)，负荷贡献分别达到53.19%、78.83%和69.69%，面积占1.20%、2.95%和14.22%。

(5) 王快水库流域BMPs对TP的削减率更高。与非PMAs相比，BMPs对PMAs的削减效率更高，约为非PMAs的4.25倍(TN)和5.37倍(TP)，验证了PMAs识别结果的准确性。

基金项目

本研究由国家重点研发计划项目资助(2024YFD700803)。

NOTES

^*通讯作者。