转
自：水土环境研习小组
综合应用水文模型和最佳管理实践（BMP）是管理流域非点源（NPS）污染的关键决策工具。优化和选择BMP方案是管理NPS污染的关键挑战，因为这些过程通常计算成本高昂，而且涉及离散-连续混合决策变量。我们的研究将一种新方法--基于多目标自适应代理建模的约束混合问题优化（MOASMOCH）--与分布式水土评估工具（SWAT）模型相结合，对中国四湖流域的BMP部署进行了有效优化。我们将优化结果与传统的非支配排序遗传算法（NSGA-II）进行了比较。结果表明，MO-ASMOCH 的计算效率明显优于 NSGA-II，与 NSGA-II 所需的 10,000 次模型评估相比，MO-ASMOCH 只需1,150次模型评估就能获得相当的帕累托最优解。这表明，对于具有离散和连续决策变量的 BMP 优化问题，MO-ASMOCH 是一种更高效的优化算法。我们选择了具有代表性的方案来计算总磷(TP)和总氮(TN)污染负荷的湖内浓度。最大的减排方案可分别减少18.3%和20.7%的TN和TP负荷，成本为1.54×10
8
元人民币。在此方案下，TN的水质级别由劣 V类提高到IV-V类，TP全年达到III类。本研究的方法可提高我们有效管理流域内非污染源污染的能力，并为环境管理实践提供有针对性的决策启示。
1.1 研究区概况
四湖流域（图1）位于长江中游和江汉平原腹地。该地区地势由西北向西南逐渐倾斜，中部地区广阔平坦。盆地面积约11,547平方公里，以耕地为主，其中 6,180 平方公里主要用于水稻种植。最初因有洪湖、长湖、三湖和白鹭湖四个湖泊而得名，但随着时间的推移，这些湖泊的面积因自然和人为因素而大幅减少，现仅存长湖和洪湖。四湖流域位于北亚热带季风湿润气候区。年日照时数约为1,800至2,000小时，年平均气温为15.9至16.6摄氏度，年降雨量为1,100至1,300毫米。大部分降雨发生在5月至8月，占全年降水量的53%。长湖和洪湖具有防洪、供水、灌溉和生态保护等多种功能。它们拥有相对丰富的水资源，是湖北省重要的农业生产基地。
2.1 SWAT模型校准
SWAT模型是一种分布式水文模型，用于模拟地表水和地下水的水量和水质。SWAT被广泛用于模拟流域内的新动力源污染控制和区域管理。为了全面模拟流域内的水平衡、侵蚀和水质，SWAT模型需要大量输入数据。表1总结了本研究中使用的具体数据类型和来源。首先，我们采用MO-ASMO方法，根据四湖流域福田寺站观测到的日径流量对SWAT模型进行了校核。校核期为2008-2019 年，其中2008年为预热期，2020-2021年为验证期。然后根据2020-2021年的月度TP和TN污染负荷对SWAT模型进行校核，以确保水质模型的准确性。我们使用纳什效率系数（NSE）和R²（判定系数）来评估模型的性能。优化结果的评价指标见表2。
表1 数据类型和来源
在SWAT模型对四湖流域NPS污染模拟的基础上，采用基于氮、磷污染物输出负荷的单位面积负荷指数法确定了关键源区。以单位面积总氮和总磷的损失强度作为NPS损失评价因子，从低到高分为5个等级。TN和TP损失强度分布如图2所示。全氮和总磷污染负荷在4个湖盆间的分布存在显著差异。靠近干渠、环湖、荆州市附近的子流域污染负荷较高。这种模式与当地人口密度、农业生产污染和外部投入污染密切相关。我们的研究重点是NPS污染从主要湖泊进入四大湖的运河。因此，没有考虑湖区和下游子盆地（57-62）。根据总氮和总磷的损失强度，确定了四大湖流域氮磷污染的关键源区为1、13、14、20和31子流域。
图2 四大湖流域总氮和总磷单位损失强度
2.2 多目标优化方法
2.2.1 MO-ASMOCH
MO-ASMOCH方法是一种基于约束混合问题代理模型的多目标优化方法，由MO-ASMO发展而来。这些算法还包括ASMO、ASMO-PODE和MO-ASMOGS。MO-ASMOCH 方法适用于涉及连续、整数和离散决策变量的复杂多目标优化问题。它包括多个进化算子，可处理不同类型的决策变量，并采用分类代理模型来处理约束函数。MO-ASMOCH方法采用高斯过程回归（GPR）来构建优化代用模型，而不是使用原始数值模型。通过在代用模型上执行多目标优化算法，可以获得近似帕累托最优解集，而代用模型的运行时间几乎可以忽略不计。因此，MO-ASMOCH方法在保证优化效果的同时，还能显著降低优化时间成本。主要步骤包括初始采样、构建代用模型和自适应采样。该方法的流程图详见图3。
图3 MO-ASMOCH方法流程图
2.2.2 NSGA-II
NSGA-II算法是NSGA方法的继承者，是一种高效的多目标优化技术，广泛应用于各个学科。它通过引入几项重大改进，解决了原始NSGA的局限性。首先，NSGA-II采用了更快的非支配排序法，可以更有效地识别个体之间的支配关系，从而提高搜索效率。其次，NSGA-II引入了拥挤度算子，用于衡量群体中个体的密度，确保足够的多样性，防止过早趋同于局部最优。最后，NSGA-II采用二元选择的精英策略，每一步都能比较并保留较好的个体，从而进一步促进群体多样性并提高算法性能。通过联合考虑个体排名和群体多样性，NSGA-II能有效探索帕累托最优解集，使其成为工程设计、资源分配和机器学习等领域的重要工具。
2.3 结构性BMP设计
结构性BMP包括众多参数，具有明确的实施成本，因此适合采用多目标优化算法进行优化。根据已发表的研究和四湖流域的实际项目，本研究选择了两种结构性BMP，即植草水道措施（BMP1）和植被过滤带措施（BMP2），本研究仅考虑BMP的建设成本。利用校准的SWAT模型，模拟两种结构性BMP对农业NPS污染的影响。目标函数为TN和TP的削减效果以及措施成本。TN和TP的减排效果与措施成本之间存在负相关关系，表明减排效果越好，措施成本越高。由于目标函数之间存在冲突关系，增强一个目标函数可能会降低另一个目标函数。因此，针对NPS污染管理优化BMP是一个典型的多目标优化问题。我们采用了两种多目标优化算法（MO-ASMOCH和NSGA-II）来优化其决策变量。在关键污染源区域部署BMP1和BMP2，同时考虑到不同子流域的空间变化，我们为五个子流域分别定义了一组决策变量（GWATL、GWATW、VFSR、FILTERW），决策变量的维度为20，目标函数为三个。本研究涉及以下多目标优化问题：
其中，y
1
为TN负荷削减率，y
2
为TP负荷削减率，C为BMP的总成本（单位：人民币）；x
1
、x
2
、......、x
20
均为结构性BMP的决策变量。在这种情况下，更高的氮和磷负荷降低率是首选，部署BMPs的成本更低是可取的。
NSGA-II算法的优化性能受种群大小和代数的影响。更大的种群规模和更高的进化世代通常会导致更好的优化结果。在本研究中，NSGA-II算法的种群大小分别设置为50和100，我们测试了不同代数对优化性能的影响。比较了NSGA-II算法在不同进化代间的收敛性。由于真正的Pareto前沿是未知的，所以最著名的近似是通过运行NSGA-II方法进行大量的目标函数评估（本研究中有40000次）得到的。MO-ASMOCH算法的优化性能依赖于初始样本的数量，由优化问题的复杂程度决定。在本研究中，我们测试了不同大小的初始样本，分别设置为600、800、1000、1200次。
为了充分评价所使用的MOASMOCH方法和NSGA-II的优化结果，我们采用了可视化和性能指标互补的方法。可视化提供了在客观空间中直接直观地显示估计的帕累托前沿和真实的（或参考的）帕累托前沿。仅仅依靠对实际Pareto前沿和两种算法的优化结果的目视检查，并不能直观和严格地评估优化结果。评估多目标优化算法得到的近似pareto最优解需要评估收敛性和多样性。收敛性表示得到的近似Pareto最优点的准确性，多样性表示解在目标空间中的分布。
理想情况下，得到的非支配解应位于真帕累托前沿，并均匀分布于真帕累托前沿。因此，我们选择利用代际距离(GD)和Hypervolume (HV)作为评价指标，定量评价两种算法的优化结果，使其优化结果更直接的比较。
其中，n代表近似帕累托最优点的数量，d
i
代表第i个最优点与真正帕累托前沿之间在目标空间中的最短欧氏距离。
2.4 BMPs的成本效益分析
经济成本是影响BMPs方案部署的一个重要因素。成本效益（CB）分析首先将环境效益计算为基线NPS污染物负荷与BMP情景下NPS污染物负荷之间的差值。然后将这种环境效益与所测量情景的计算成本之比定义为CB值。CB代表单位成本可实现的污染负荷减少量，是衡量BMPs项目有效性和经济成本的重要指标。CB值越高，表明BMP场景的成本效益越好。CB的计算公式如下：
其中，LOAD
base
和LOAD
BMP
分别表示基线和BMP情景下的污染物负荷（单位：kg）；成本为BMPs情景对应的措施成本，单位：万元。
3.1 MO-ASMOCH与NSGA-II优化结果的比较
3.1.1 优化方法设置的影响
对于MO-ASMOCH方法，当总模型评估从750次增加到1,150次时，所获得的非支配解在图4(a)中得到了显著改善。此外，还观察到溶液有更好的扩散。在900次运行时，优化结果似乎还没有达到收敛。然而，当运行次数达到1150次时，优化结果与参考的Pareto front基本一致，表明已经实现了收敛。图4(b)和(c)进一步显示了优化结果的TN还原率-成本和TP还原率-成本的平面投影。可以看到，当运行次数达到1,150次时，投影后的两条帕累托曲线与参考的帕累托锋基本重合，表明已经达到收敛。
在模型评估时，观察到帕累托最优前沿的显著改进，以及解决方案从2000增加到10000，如图5(a)所示，说明了NSGA-II的优化结果。当运行次数达到10000次时，优化结果与参考的帕累托前沿基本一致，表明已经达到收敛。对三维目标函数空间进行平面投影，得到TN还原率-成本与TP还原率-成本的优化关系，如图5(b)和(c)所示。从图中可以看出，当运行次数达到10000次时，投影后的两条帕累托曲线与参考的帕累托锋几乎重合。
节不同主要和次要地形属性的空间不确定性，因此异质站点比同质站点的误差更大。Grohmann利用穿梭雷达地形图任务(SRTM)的数据，对南美洲进行了一项区域尺度研究，报告了类似的发现。这可能是由于特定地点的地形和坡度的对称或不对称分布造成的；在任何均质地点，误差预计会比异质地点更呈正态分布，因为在异质地点，误差是倾斜的，并因分辨率更粗而恶化。此外，所研究的地形属性或指数是地表特征的数学量化受到每个地表不确定性的影响。
除了站点之间的差异外，每个站点内地形属性的显著不确定性也因地形属性的基础空间分辨率而异，因为MAE的不确定性随着空间分辨率的提高而增加。当分辨率大于15米时，原始值在统计上有明显差异，特别是在异质站点。这种情况可能是由于栅格像素数量的减少与之前讨论过的地表特定不确定性相结合产生的累积误差造成的。这是因为在重新采样基点高程并将其插值到各种地形属性时会产生周期性数学误差。理论上，这些误差会累积起来，正误差会被负误差抵消，并呈现在最终输出中。这些发现得到了以往研究的广泛支持。然而，与此形成鲜明对比的结果是，更高的分辨率会产生相当大的误差，尤其是在使用激光雷达创建凹陷识别输入数据的情况下。出现这种情况的原因可能是：1）LiDAR 数据采集过程中产生的伪影；2）研究边界以外的点产生的凹陷或不想要的表面特征。后一种影响可能是由于在特定区域进行研究时原始DEM值发生了改变，尽管与异质区域相比，同质区域的研究地点显示出最小的离散性。有鉴于此，可能值得研究点密度和预处理的不同组合，包括先验设置边界，这项任务超出了当前研究的范围，仍然存在局限性。
不同地形属性的值随地块测量程序的不同而变化很大，而且往往与分辨率有明显的交互作用。一般曲率和地形位置指数仅受地块测量程序的影响，地块平均值的误差低于从地块中心提取的值。地块测量程序和分辨率之间的相互作用影响了坡度和TRI属性的误差。Smith 等人分析了DEM分辨率和邻域大小对数字土壤测量精度的类似影响。从表面上看，任何地貌的对比性质都会造成这种情况。从缩放的角度来看也是如此；通过平均地块的数值，可以在一定程度上缓解对比性。Hodgson也分析了类似的现象，并建议采用更精细的分辨率或地块平均值来确定合理的数值。总体而言，15米之前的分辨率影响很小，这与本研究之前的发现一致，也解释了本分析的累积效应。有趣的是，海拔高度既不受地块测量程序的影响，也不受分辨率的影响。这可能是由于局部最小值和最大值的影响，因为最高测量值会影响地块内的地块值。另一方面，本研究仅使用了静态邻域大小来限制地块，可以使用移动窗口类型分析或增加地块大小来进一步扩展和测试，以检查局部效应。
图4 MO-ASMOCH方法的多目标优化结果：(a)三维目标函数空间；(b) TN的目标函数空间平面投影；(c) TP目标函数空间的平面投影。
图5 NSGA-II方法的多目标优化结果：(a)三维目标函数空间；(b) TN的目标函数空间平面投影；(c) TP目标函数空间的平面投影。
3.1.2 两种方法的比较
图6更直观地展示了两种多目标优化算法的运行次数与评价指标之间的关系。图中分别显示了随着TN和TP运行次数的增加，Hypervolume和GD的变化情况。很明显，MO-ASMOCH的收敛速度比NSGA-II快得多。在评估次数达到 1,150 次后，MO-ASMOCH的Hypervolume和GD指标变化都很小。此外，MO-ASMOCH运行SWAT模型一千次的结果与NSGA-II运行同一模型一万次以上的结果非常接近。在确定MO-ASMOCH和NSGA-II的收敛点后，我们进一步比较了两种算法收敛时的帕累托前沿。从图7可以看出，在计算预算相同的情况下，NSGA-II的结果比MO-ASMOCH差很多。MO-ASMOCH算法在评估1,150个模型后的结果与NSGA-II算法在评估10,000个模型后的结果非常相似，后者的表现稍好一些。
图7 MO-ASMOCH与NSGA-II优化结果的比较
我们对MOASMOCH和NSGA-II的评价指标进行了比较，以探索它们的优化结果，如图8所示。我们观察到，对于TP优化，MO-ASMOCH在运行1,150次后获得的Hypervolume值优于运行10,000次的NSGA-II。同样，运行1,150 次的 MO-ASMOCH的超体积值也优于运行1,150次和10,000次的NSGA-II。不过，在TN优化方面，运行10,000次的NSGA-II比运行1,150次的 MO-ASMOCH获得了更好的生成距离值。与运行1,150次的NSGA-II相比，运行10,000次的NSGA-II的GD值要好得多。总体而言，运行1,150次的 MOASMOCH和运行 10,000 次的NSGA-II的优化结果非常接近，且优于运行 1,150 次的NSGA-II。这表明，与NSGA-II相比，MOASMOCH可以显著提高计算效率并降低时间成本。
图8 MO-ASMOCH与NSGA-II评价指标的比较
3.2 情景分析
BMP优化后，在关键源区TN的最大减少率约为20%。对成本在200~25,000,000元人民币之间的情景进行了CB分析。结果表明，在同等成本下，优化BMPs
1
和BMPs
2
对TN比对TP更具成本效益。因此，我们使用TN的Pareto边界来选择场景设计。在优化结果对比的基础上，通过1150次模型评价，从MO-ASMOCH优化结果中选择流域NPS管理措施方案解集。TN还原率分别为5.54%、9.75%和18.33%。然后从帕累托前线选择经济成本效益最高的点。具体来说，当所需的TN还原率超过5.54%时，图底部对应的解决点的建设成本为4.09 ×10
7
元。当TN还原率超过9.75%时，图中对应的解决点的建设成本为9.07×10
7
元。当TN还原率达到18.33%时，图顶部对应的解决点的建设成本为15.4×107元。这些上、中、下解决方案点被建立为最优管理措施的场景点。方案1 ~方案3的TN减量率及相应成本如表3所示。
表3 从Pareto前端选择的代表性BMP场景
为了更直观地评估各方案对水质的影响，对方案1至3进行了模拟，以评估湖中TN和TP的浓度。图10和图11从上到下显示了水质浓度曲线，代表了不同方案下湖水中的污染物浓度。方案1-3的工程措施成本分别为4.09×10
7
、9.07×10
7
和15.4×10
7
。其中，方案3中TN的水质分级由劣V类降为IV-V类，TP 全年水质分级达到III类。如果进一步考虑减少化肥投入量的可行性和污染物减排的效果，在方案3的基础上减少25%的化肥投入量，水质将得到显著改善。具体而言，这种组合可使TN的全年水质等级达到IV级或以下，TP的全年水质等级达到III级或以下。
在本研究中，我们使用SWAT模型评估了两种结构性BMP（植草水道和植被过滤带）在减轻四湖流域的NPS污染方面的功效。一种基于代用建模的新型多目标优化方法MO-ASMOCH与SWAT模型相结合，以获得具有连续和离散决策变量的BMP的帕累托最优解。将该方法的优化结果与经典的NSGA-II方法进行了比较。主要结论如下：
(1)根据校准的SWAT模型，使用两种多目标优化方法对两种结构性BMP的部署进行了优化。与NSGA-II相比，MO-ASMOCH方法的效率明显更高。具体来说，MOASMOCH方法至少需要1,150次模型评估才能获得与 NSGA-II 方法相当的非优势解，而 NSGA-II 方法需要10,000次模型评估。这表明，MO-ASMOCH方法实现类似优化结果所需的计算时间仅为 NSGA-II的11.5%。MO-ASMOCH是一种高效的多目标优化算法，可用于设计各种流域情况下的BMP。
(2)使用MO-ASMOCH和NSGA-II方法对BMP1和BMP2进行优化，结果表明TN的成本效益高于TP。根据优化结果，确定了三个具有代表性的成本优化方案，以最低成本实现5.54 %、9.75 %和18.33 %的TN减排率，每个方案都符合利益相关者的利益，同时为实际决策者提供支持。
(3)结果表明，在仅实施结构性BMP的情况下，高成本方案可将TN的水质级别从劣V级提升至IV-V级。此外，在该方案下，TP的水质等级也达到了III级。如果同时考虑结构性和非结构性BMP，方案3以及减少25%的化肥施用量可将TN的水质等级全年维持在IV级或以下，同时将TP的等级全年维持在II~III级。
总之，MO-ASMOCH方法为优化BMP提供了一种高效、可扩展的解决方案，因此在大规模流域中尤为重要，因为在这些流域中，传统方法的计算成本可能过高。这种用户友好型方法可随时适用于各种流域，从而在优化资源利用的同时，针对当地特定的污染挑战制定量身定制的BMP战略。环境机构可以利用MOASMOCH框架来制定综合管理计划，优先采用影响较大的BMP。
本期编辑：王若文
指导教师：吴磊
西北农林科技大学
水利与建筑工程学院
