示例 2：四湖流域的最佳管理实践（BMP）优化

背景介绍

当流域受非点源污染影响时，在SWAT模型的基础上应用最佳管理实践（Best Management Practices，BMPs）进行流域污染管理，已被证明是一种可靠且有效的方案。本例参考了 Long 等（2025）对四湖流域最佳管理实践的研究成果。

四湖流域概况

四湖流域位于长江中游、江汉平原腹地。为研究该流域的水质输移过程，本例构建了关于四湖流域的SWAT模型。构建过程使用的数据包括：

DEM（数字高程模型） - ASTER GDEM，分辨率为 30 米
土地利用 - 中国土地利用/覆盖变化（CNLUCC）数据集
土壤数据 - 南京土壤研究所的第二次全国土壤普查（1:100 万）
气象数据 - 中国区域地面气象因子驱动数据集
径流观测 - 水文年鉴数据（2008.01.01 至 2021.12.31）
水质观测 - 中国环境监测总站数据（2020.11 至 2021.12）

本例省略了径流和水质的校准过程，重点介绍最佳管理实践（BMP）的优化过程。在中国，总氮（TN）和总磷（TP）是评估湖泊水质的重要指标。以下是四湖流域在2021年12月31日的TN及TP分布图：

四湖流域内TN及TP分布情况

在本例的SWAT模型，四湖流域的主湖区位于子流域 32，而其主要入湖流量来自子流域 51。此外，基于上图，确定污染重点管理区域主要包括子流域1、13、14、20及31。

工程管理优化

总览

SWAT模型内置了多种工程管理措施（BMP），例如梯田操作（BMP1）、瓦管排水（BMP2）、植被过滤带（BMP4）、植草水道（BMP7）等。

为了减少TN和TP排放，植被过滤带（BMP4）和植草水道（BMP7）最为常用。考虑布置成本，本例计划仅在关键子流域（1、13、14、20及31）布置植被过滤带和植草水道。

SWAT项目中的 .ops文件用于控制BMPs的基础设置，其中，植被过滤带的相关参数包括：

FILTER_I：是否启用植被过滤带（1 启用，0 不启用）
FILTER_RATIO：农田面积与植被过滤带面积的比例（单位：ha/ha，范围：0–300）
FILTER_CON：植被过滤带中高密植区域占比（未启用）
FILTER_CH：通道化流量比例（未启用）

植草水道的相关参数包括：

GWATI：是否启用植草水道（1 启用，0 不启用）
GWATN、GWATSPCON：水力特性（略）
GWATD：沟渠深度（默认 3/64 * GWATW）
GWATW：植草水道宽度（单位：米）
GWATL：植草水道长度（单位：千米）
GWATS：坡度（单位：m/m）

为简化问题，本例每个子流域仅优化5个关键参数（即FILTER_I、FILTER_RATIO、GWATI、GWATW、GWATL），因此共计25 个变量。而优化目标则包括：目标1：TN减少率、目标2：TP减少率、目标3：实施成本。总体而言，此问题是一个典型的 混合型多目标优化问题。

参数信息

本例所优化的参数信息包括：

名称	类型	范围	单位
FILTER_I	整数	0–1	无
FILTER_RATIO	浮点数	1–300	无
GWATI	整数	0–1	无
GWATW	离散	1, 5, 10, 15, 20, 25, 30	米
GWATL	浮点数	10–1000	千米

由于每个子流域的BMPs设置不同，建立参数文件时，需要为每个子流域单独定义这5个参数。

同时注意，离散型参数GWATW的所有可能值需通过下划线分隔列在 Min_Max 字段中：

GWATW v d 1_5_10_15_20_25_30 1

完整的参数文件如下：

文件名：para_bmp.par

Name Mode Type Min_Max Scope
FILTER_I v i 0_1 1
FILTER_RATIO v f 1_300 1
GWATI v i 0_1 1
GWATW v d 1_5_10_15_20_25_30 1
GWATL v f 10_1000 1
FILTER_I v i 0_1 13
FILTER_RATIO v f 1_300 13
GWATI v i 0_1 13
GWATW v d 1_5_10_15_20_25_30 13
GWATL v f 10_1000 13
FILTER_I v i 0_1 14
FILTER_RATIO v f 1_300 14
GWATI v i 0_1 14
GWATW v d 1_5_10_15_20_25_30 14
GWATL v f 10_1000 14
FILTER_I v i 0_1 20
FILTER_RATIO v f 1_300 20
GWATI v i 0_1 20
GWATW v d 1_5_10_15_20_25_30 20
GWATL v f 10_1000 20
FILTER_I v i 0_1 31
FILTER_RATIO v f 1_300 31
GWATI v i 0_1 31
GWATW v d 1_5_10_15_20_25_30 31
GWATL v f 10_1000 31

💡 提示： 参数文件支持重复名称参数，因为SWAT-UQ将通过文件中参数出现的顺序进行唯一ID编号。

目标函数

本例包含三个目标函数：

1.TN减排率： $$ Obj_1 = \frac{TN_{base} - TN_{now}}{TN_{base}} $$ 表示应用BMPs前后，子流域51对主湖区的TN流出总量的减少率。

2.TP减排率： $$ Obj_2 = \frac{TP_{base} - TP_{now}}{TP_{base}} $$ 表示应用BMsP前后，子流域51对主湖区的TP流出总量的减少率。

3.BMPs布置成本：

植被过滤带单价：420元/ha
植草水道单价：200元/ha

因此，子流域 $i$ 的实施成本表示为：

$$ cost_{filter}^i = Area_{AGRI}^i \times FILTER_RATIO \times FILTER_I \times 420 $$

$$ cost_{gwat}^i = GWATW \times GWATL / 10 \times GWATI \times 200 $$

总成本为：

$$ Obj_3 = \sum_{i \in {1,13,14,20,31}} (cost_{filter}^i + cost_{gwat}^i) $$

本例三个目标函数的定义无法通过评估文件直接完成。但先通过评估文件提取模拟结果中的关键数据，进而为用户自定义的目标函数userObjFunc或约束函数userConFunc提供支持。

评估文件内容如下：

文件名：obj_bmp.evl

SER_1 : ID of series data
OBJ_1 : ID of objective function
WGT_1.0 : Weight of series combination
RCH_51 : ID of RCH, or SUB, or HRU
COL_42 : Extract Variable. The 'NUM' is differences with *.rch, *.sub, *.hru.
FUNC_7 : Func Type ( 1 - NSE, 2 - RMSE, 3 - PCC, 4 - Pbias, 5 - KGE, 6 - Mean, 7 - Sum, 8 - Max, 9 - Min )
2021/01/01 to 2021/12/31 : 

SER_2 : ID of series data
OBJ_2 : ID of objective function
WGT_1.0 : Weight of series combination
RCH_51 : ID of RCH, or SUB, or HRU
COL_43 : Extract Variable. The 'NUM' is differences with *.rch, *.sub, *.hru.
FUNC_7 : Func Type ( 1 - NSE, 2 - RMSE, 3 - PCC, 4 - Pbias, 5 - KGE, 6 - Mean, 7 - Sum, 8 - Max, 9 - Min )
2021/01/01 to 2021/12/31 : Period for data extraction

通过Python代码自定义userObjFunc函数：

# 定义基准 TN / TP 值（用于归一化计算）
TN_Base = 3.314e7  # 总氮基准值（单位视情况而定）
TP_Base = 3.717e6  # 总磷基准值

# 定义应用 BMP（最佳管理措施）措施的子流域编号列表
Basins = [1, 13, 14, 20, 31]

def userObjFunc(attr):
    """
    用户自定义目标函数。

    参数：
    - attr: dict
        包含决策变量、目标值、约束值、时间序列、HRU 信息等。

    返回：
    - objs: np.ndarray
        包含 3 个目标函数值的数组 [obj_1, obj_2, obj_3]
    """

    objs = np.zeros(3)  # 初始化三个目标值
    x = attr["x"]       # 获取决策变量

    # 第一个目标函数：TN 相对减少率
    objs[0] = (TN_Base - attr['objs'][1]) / TN_Base

    # 第二个目标函数：TP 相对减少率
    objs[1] = (TP_Base - attr['objs'][2]) / TP_Base

    # 第三个目标函数：BMP 实施的总成本
    HRUInfosTable = attr["HRUInfos"]
    cost = 0

    for i, ID in enumerate(Basins):
        # 计算当前子流域的总面积
        areas = np.sum(
            HRUInfosTable.loc[
                (HRUInfosTable.SUB_ID == ID),
                "Area"
            ].tolist()
        )

        # 从决策变量中提取 BMP 设计参数
        filter_I = x[5 * i]          # 滤带开关
        filter_ratio = x[5 * i + 1]  # 滤带处理比例
        graw_I = x[5 * i + 2]        # Graw BMP 开关
        graw_W = x[5 * i + 3]        # Graw BMP 宽度
        graw_L = x[5 * i + 4]        # Graw BMP 长度

        # 滤带措施成本计算（单位成本：420 元/公顷）
        cost_filter = areas * filter_ratio * filter_I * 420

        # Graw 措施成本计算（单位成本：200 元/公顷）
        cost_graw = graw_W * graw_L * graw_I / 10 * 200

        cost += cost_filter + cost_graw  # 累计总成本

    objs[2] = cost  # 第三个目标函数为总成本
    return objs

至此，完成所有准备工作。

执行优化

基于Python环境的执行优化示例如下：

import numpy as np
from swatuq import SWAT_UQ
from UQPyL.optimization.multi_objective import NSGAII

nInput = 25
nOutput = 3

projectPath = "E:\\BMPs\\TxtInOut"        # SWAT 项目路径
exeName = "swat.exe"                      # 执行程序名
workPath = "E:\\DJ_FSB"                   # 工作路径
paraFileName = "para_bmp.par"            # 参数文件名
evalFileName = "obj_bmp.evl"             # 评价文件名
specialFileName = "special_paras1.txt"   # 特殊参数文件名

problem = SWAT_UQ(
    projectPath = projectPath,
    swatExeName = exeName,
    specialFileName = specialFileName,
    workPath = workPath,
    paraFileName = paraFileName,
    evalFileName = evalFileName,
    verboseFlag = True,
    numParallel = 10,
    userObjFunc = userObjFunc,
    nInput = 25,
    nOutput = 3,
    optType = ["max", "max", "min"]  # 目标优化类型：最大化、最大化、最小化
)

nsgaii = NSGAII(
    nPop = 100,
    maxFEs = 20000,
    saveFlag = True,
    verboseFlag = True,
    verboseFreq = 5
)

nsgaii.run(problem = problem)

# 最终结果将保存在目录：Result\Data\NSGAII_SWAT-UQ_D25_M3.hdf

BMPs优化的结果如下：