考虑地块拓扑关系的流域最佳管理措施优化方法与流程

本发明属于农业非点源污染模拟技术和最佳管理措施优化技术领域，具体涉及一种考虑地块拓扑关系的流域最佳管理措施优化方法。

背景技术：

最佳管理措施(Beneficial management practices，BMP)是为了控制和削减非点源污染，防止水土流失、保护土壤资源，改善水质及流域生态环境而采取的一系列措施。在实际流域中进行BMP选择和优化配置时，通常上游管理单元措施的选择一定程度上由其下游单元措施的选择决定，例如下游管理单元实施了草地过滤带措施，则上游管理单元就不用实施等高耕作或者免耕措施，因此，BMP优化应该考虑管理单元的上下游关系。

现有的BMP优化方法有两种，一种是根据专家经验和关键源区识别，在一定程度上考虑了管理单元上下游关系来设计BMP方案，然后通过流域模型模拟评价进行择优选择，但是该方案评价时大多采用半分布式流域模型，无法模拟方案中管理单元的空间相互作用(如上游排污对下游的影响)；另一种是基于智能优化算法对流域BMP空间配置方案进行优化，但是由于优化算法以随机法生成初始种群(初始解)，根据模型评价结果来引导产生新种群(新解)进行BMP配置优化，而使用的模型在现有研究中多是采用半分布式流域模型，采用的随机法也忽略了管理单元的上下游关系。综上所述，现有的两种BMP优化方法没有充分考虑管理单元的上下游关系及其空间相互作用。

基于优化算法在流域中进行BMP空间优化配置时，如果考虑管理单元的上下游关系及管理措施空间相互作用的先验知识进行优化，并且采用能够体现空间单元之间物质交互作用的全分布式流域模型作为BMP评价模型，那么不仅可提高优化算法的收敛速度，而且可以提高最优解的质量。

技术实现要素：

鉴于上述，本发明提供了一种考虑地块拓扑关系的流域最佳管理措施优化方法，其将现有的流域最佳管理措施空间相互作用知识应用于优化算法初始，以提高算法搜索的起点。

一种考虑地块拓扑关系的流域最佳管理措施优化方法，包括如下步骤：

(1)将流域整体划分成具有上下游关系的多个独立管理的地块，并构建地块树，所述的地块具有单一的土地利用类型和明确的上下游关系；

(2)初始化生成具有一定数量规模的种群，种群中的每一个体为一套潜在的流域整体管理方案，即对应一棵地块树且其中每一地块通过空间相互作用规则确定了对应的BMP；

(3)根据泥沙削减率和成本利用快速非支配排序法对种群中的个体进行排序，按次序使个体两两配对进行交叉变异，从而生成下一代种群；

(4)根据步骤(3)不断进行遗传进化，直至达到设定的迭代终止条件，最终获得的新一代种群即为精英集，进而根据流域治理的环境目标以及预算从精英集中选取出一个个体作为流域整体的最优管理方案。

所述步骤(1)的具体过程如下：

1.1利用DEM(数字高程模型)建立流域整体的流向图；

1.2根据流向关系将流向图中的所有栅格对应作为节点组建成树形结构；

1.3根据土地利用图使树形结构中具有相同土地利用类型以及直接流向关系的节点进行合并，合并后的树形结构即为所述的地块树且其中的每一节点即对应所述的地块。

所述步骤(2)的具体过程为：

2.1以流域出口所在的地块作为地块树的根节点，根据土地利用类型为根节点选择合适的BMP或随机生成根节点的BMP；

2.2对于地块树中除根节点以外的任一地块，若其直接流向的下游地块BMP为无措施，则从无措施以及其他具体措施中随机选择其一作为该地块的BMP；若其直接流向的下游地块BMP为某一具体措施，则令该地块的BMP为无措施；依此确定地块树中各地块的BMP，从而得到一个个体；

2.3根据步骤2.1和2.2执行多次，生成多个个体，从而组成种群。

所述步骤(3)中个体泥沙削减率的计算表达式如下：

其中：F_sed(X)表示按个体X中各地块BMP实施后流域整体的泥沙减少率，V(0)为各地块不采取任何措施下流域整体受侵蚀所产生的泥沙模拟量，V(X)为按个体X中各地块BMP实施情况下流域整体受侵蚀所产生的泥沙模拟量。

所述步骤(3)中个体成本的计算表达式如下：

其中：F_cost(X)表示按个体X中各地块BMP实施的总费用，C(x_i)表示个体X中第i个地块的单位面积BMP实施费用，A_i表示第i个地块的面积，n为地块总数。

所述步骤(3)中使个体两两配对进行交叉变异的具体过程为：对于任一对个体，从其中一个个体A中除根节点外任意选定一个地块a，同时从另一个个体B中选定一个地块b且地块b与地块a在各自个体中的位置相同；使个体A中以地块a为根节点的子树与个体B中以地块b为根节点的子树进行相互交换，得到个体A'和个体B'，从而完成交叉过程；对于交叉后得到的新个体，随机从中选取若干个体，从被选取的新个体中随机选取一地块，并随机选择一种具体措施或无措施替换该地块原有的BMP，从而完成变异过程。

本发明充分利用流域管理单元(即地块)的上下游关系及BMP的空间相互作用的知识进行种群初始化，与随机法生成初始解的传统方法相比，本发明可提高解的质量，有助于算法收敛，快速找到较优的解。

附图说明

图1为本发明方法的优化流程示意图。

图2为本发明中地块划分的流程示意图。

图3为本发明中种群初始化的流程示意图。

图4为本发明中遗传算法的交叉策略示意图。

图5为本发明方法与传统方法的实验步骤对比示意图。

图6(a)为种群规模为60情况下本发明方法与传统方法的实验结果对比图。

图6(b)为种群规模为100情况下本发明方法与传统方法的实验结果对比图。

图6(c)为种群规模为200情况下本发明方法与传统方法的实验结果对比图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

本发明方法以遗传算法为基础，以流域模拟技术为手段，结合现有的流域最佳管理措施BMP空间相互作用知识，充分考虑管理单元的空间拓扑关系，即上下游关系，进行BMP的优化，其方法具体步骤如图1所示：

(1)将流域整体划分成具有上下游关系的流域管理单元，以地块为管理单元，提出一种新的地块划分方法：根据土地利用和DEM提取的流向图，重新划分地块，划分的地块具有单一的土地利用类型和明确的上下游关系；该步骤具体操作如图2所示：

1.1根据流向图构建流向关系树形结构，同时记录在树形结构中每个栅格的ID号(ID与空间位置相对应)；

1.2将土地利用图每个栅格的ID号与树形结构中相同ID的栅格叠加，形成具有土地利用信息和流向关系的树形结构；

1.3将相同土地利用的栅格合并，最终形成具有上下游关系、单一土地利用的地块。

(2)选择一个多目标遗传算法(如ε-NSGA-II)，并构建一个分布式流域模型(如Landscape SWAT模型)和BMP成本计算组件。

(3)提取流域最佳管理措施BMP之间的空间相互作用规则，并将其应用于遗传算法的种群初始化过程，即在地块上按照这些规则生成一定数量的BMP空间配置情景；该步骤具体操作如图3所示：

3.1读取具有上下游关系的地块信息，构建地块树；

3.2找到流域出口所在的地块(树形结构的根节点)，该地块是其他地块的下游地块，根据该地块的土地利用类型，随机生成或者指定某种合适的BMP；

3.3从根节点的孩子节点开始，按照下游到上游的地块关系遍历(树的先序遍历)每个地块，先判断其下游地块是否实施了BMP，如果实施，则该地块不实施BMP，如果没有实施，则根据该地块的土地利用类型，随机生成或者指定某种合适的BMP。假设有3种可选的BMP，分别为：无措施、免耕、退耕还草，则我们有如表1所示的BMP相互作用规则：

表1

3.4根据用户指定的种群规模(个体数目)，按照步骤3.2～3.3，生成种群。

(4)多目标评价(个体适应度计算)，即调用流域模型和成本计算组件对每个BMP情景的多个目标(如减沙率最大、成本最小)进行计算。其中多目标评价的目标包括环境目标和经济目标，优化问题描述为最大化泥沙削减率和最小化BMP成本；进行BMP情景多目标评价时，由优化算法产生的新的情景需要与一个基准情景(baseline scenario)进行比较，基准情景是指当前未实施BMP时的模型模拟结果；这里的模型模拟结果是指流域模型在参数率定后，泥沙流入河道的模拟值。

每个BMP情景的泥沙削减率计算是通过该情景的模型模拟值与基准情景做差值来计算，具体计算公式如下：

其中，F_sed(X)表示实施BMP后的泥沙减少率(％)，V(0)是基准情景下侵蚀所产生的泥沙模拟量(kg)，V(X)是某个BMP情景下侵蚀所产生的泥沙模拟量(kg)。

BMP情景的成本计算采用如下公式：

其中，F_cost(X)是该BMP情景的费用(元)，C(x_i)表示在第i个地块上单位面积BMP实施费用(元/公顷)；如果该地块上没有实施BMP，则费用为0；如果实施了BMP，则按照该BMP单位面积费用计算；A_i表示第i个地块的面积。BMP费用数据是通过野外调查或者查相关文献获得；本实施例中用到的BMP成本如表2所示：

表2

(5)根据个体适应度计算结果，对种群中的所有个体采用快速非支配排序法进行排序，排序得到的非支配解利用ε-支配调整策略进行调整，并保留到一个集合中(称为精英集)。

(6)判断是否达到算法的终止条件(如最大进化代数)，如果没有达到终止条件，则对种群进行交叉、变异操作，产生新的种群，并转到步骤4；如果达到了终止条件，则算法停止执行，最终获得的精英集为Pareto最优集或近似最优集。

其中遗传算法种群的交叉策略采用树形编码的交叉策略，算法在遗传操作时的交叉算子采用随机选择某个树节点的策略，然后将父代对应树节点交叉，并且树结构保持不变，如图4所示。交叉算子在随机选择树节点后，首先判断父代1和父代2以该节点为根的子树BMP编码是否相同，如果不同则交换两棵子树；如果相同，则重新选择树节点，再比较父代以该节点为根的子树BMP编码是否相同，原因是当两棵子树的编码相同时，无需交换。变异算子根据用户给定的变异概率和随机数(范围为0～1)来确定对种群中的个体是否变异，如果随机数小于用户给定变异概率，则该个体不执行变异操作，否则执行变异操作，即采用随机选择法挑选树的某个节点号，并以随机法生成新的BMP类型代替该节点原来的BMP类型。

以下我们在其他条件完全相同的情况下，使本发明方法与传统方法进行对比验证，实验步骤如图5所示。由于对比实验主要涉及ε-NSGA-II算法的种群初始化过程，因此，种群规模参数可能对优化结果有所影响。考虑到上述情况，设计了种群规模分别为60、100、200的三组对比实验，60、100、200均为遗传算法中较为常见的种群规模参数值。在设计的三组实验中，地块数目为79；ε-NSGA-II的主要参数设置如下：最大进化代数为100代，变异概率为0.1，ε(网格大小)为0.075。

当种群规模分别为60、100、200时，本发明方法与传统的方法的BMP优化结果如图6所示。结果显示，本发明方法获得的Pareto最优解均比传统的方法好，此结果表明本发明方法在考虑地块空间拓扑关系进行BMP优化，可提高解的质量，有助于算法收敛，快速找到较优的解。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。