CN107992965A - 基于粒子群算法的小流域水土资源协调优化方法 - Google Patents

Abstract

本发明属于水利工程领域，尤其涉及基于粒子群算法的小流域水土资源协调优化方法，包括步骤一：构建水土资源优化配置线性规划数学模型；步骤二：确定水土资源优化方案的空间配置；步骤三：基于粒子群的智能优化求解。步骤一包括确定决策变量、目标函数、约束条件、构建完整的数学模型。本发明的有益效果：优化的目标是求目标函数的最大值即最大的经济效益，通过调整决策变量的价值系数而实现优化目的，以水土资源定田，使规划结果向所希望的方向倾斜，基于粒子群的智能优化求解采用单纯形微粒群算法，继承了微粒群优化算法原有优点的同时，不用进行复杂的矩阵计算，计算工作量小，优化速度快。

Description

基于粒子群算法的小流域水土资源协调优化方法

技术领域

本发明属于水利工程领域，尤其涉及基于粒子群算法的小流域水土资源协调优化方法。

背景技术

资源配置的首要问题就是要看配置客体是否具有稀缺性和多用途性,显然土地和水都具有这两种特性。在以往多是将水资源和土地资源分别进行优化配置。土地资源的配置,是为了达到一定的生态经济最优目标,依据土地特性和土地系统原理，依靠一定的技术和管理手段，对区域有限的土地资源的利用结构、方向,在时空尺度上,分层次进行安排、设计、组合和布局，以提高土地利用效率和效益，维持土地生态系统的相对平衡，实现土地资源的可持续利用。水资源优化配置是一定时期内，在有效、公平和可持续的原则基础上,对特定区域内有限的、不同形式和质量的水资源,通过工程措施与非工程措施在各用水对象之间进行科学分配。比较而言,土地资源优化配置研究较多,水资源优化配置随着水资源供需矛盾日益突出，逐渐引起人们的重视。

水土资源优化配置是利用当地小流域有限的水(坡面径流)、土(坡面土地)资源条件，合理配置蓄水工程，以获取最大的经济、社会效益；水土资源优化配置包括两个层次：第一个层次是水土资源的空间优化配置，第二个层次是蓄水系统的空间优化配置。在第一个层次水土资源的空间优化配置方面，尽管人们已经在水土资源调查、评价、规划，以及以水土资源综合利用为核心的区域经济、社会、生态环境协调发展等方面开展了大量工作，但总体而言,目前将两类资源结合起来所做的研究还比较少，已有研究也多偏重于水资源的配置，与土地资源的结合不够。

发明内容

为解决上述问题，本发明提供基于粒子群算法的小流域水土资源协调优化方法。

本发明的技术方案：基于粒子群算法的小流域水土资源协调优化方法，其特征在于包括以下步骤：

步骤一：构建水土资源优化配置线性规划数学模型；

步骤二：确定水土资源优化方案的空间配置；

步骤三：基于粒子群的智能优化求解。

进一步地，所述步骤一中构建水土资源优化配置线性规划数学模型，包括以下步骤:

步骤A:确定决策变量，包括x₁:一类水平梯田；x₂:二类水平梯田；x₃:栽种林木；x₄:汇水的天然林草地；

步骤B:确定目标函数，表达式为maxZ＝c₁x₁+c₂x₂+c₃x₃+c₄x₄；其中c_i为利润净增值/(hm²·a)(i＝l,2,3,4)；

步骤C:确定约束条件，约束条件方程组为：

式中，e为水平梯田每亩每年的灌溉需水量(m/a)×10⁴m²；Yuxi为鱼鳞坑坡地每亩每年的产流量(m/a)×10⁴m²；caodi为天然林草地每亩每年的产流量(m/a)×10⁴m²；

步骤D:构建完整的数学模型，数学模型表示

进一步地，c_i＝收入/(hm²·a)-支出/(hm²·a)。

进一步地，所述步骤D中构建完整的数学模型是以人定田方案的水土资源优化配置数学模型；

以人定田方案数学模型的初始约束条件方程组表示如下:

进一步地，所述步骤二中确定水土资源优化方案的空间配置，包括以下步骤：

步骤E:确定空间优化配置水土资源的原则，包括就低原则、就平原则、就水原则、就现状原则；

步骤F:细化具体配置原则，包括配置顺序、坡度要求、地形要求和地貌要求。

进一步地，所述步骤三中基于粒子群的智能优化求解，采用标准微粒群算法或单纯形微粒群算法。

进一步地，所述单纯形微粒群算法具体采用单纯形微粒群优化算法。

进一步地，所述单纯形微粒群优化算法的流程包括以下步骤：

1)确定目标函数；

2)初始化设置SMPSO的微粒群个数p、寻优代数、加速因子、惯性权重系数等参数值，随机初始化微粒群的位置和速度；

3)由具体问题确定微粒群的空间维数，用PSO的随机初始微粒作为单纯形搜索法的初始点进行单纯形搜索，按目标函数评价各微粒的初始适应值，将此时微粒群中各微粒的位置作为P_best，从P_best中找到的具有最好适应值的微粒的位置作为g_best；

4)按照公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，对微粒群p和p_best中各微粒的速度和位置进行更新；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)；

5)将更新后的微粒群p和P_best分别作为单纯形搜索法的初始点进行单纯形搜索，把搜索后的微粒群p和P_best中的每个微粒的适应值分别进行比较，取适应值较优的微粒的位置为新的p_best；

6)比较微粒群P_best和g_best的每个微粒的适应值，取适应值较优的微粒的位置为新的g_best；

7)若满足停止条件，搜索停止，输出全局历史最优位置和全局历史最优适应值为所求结果；否则，返回步骤4继续搜索。

进一步地，所述标准微粒群算法的流程包括以下步骤：

a)确定目标函数；

b)设定微粒群微粒个数、寻优代数、加速因子、惯性权重系数参数值，随机初始化各微粒的位置和速度；

c)由所要优化的目标函数的设计变量数确定粒子的空间维数，由目标函数值计算每个粒子的适应值，并计算P_best和g_best；

式中，P_best为微粒目前为止由自己发现的最优位置；g_best为整个群体中所有微粒发现的最优位置；

d)比较每个微粒的适应值与个体极值P_best，如果f(P_i)＞f(P_best)，则P_best＝P_i；

e)比较每个微粒的适应值与全局极值g_best，如果f(P_i)＞f(g_best)，则g_best＝P_i；

f)根据公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，更新每个微粒的飞行速度和空间位置；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)。

g)判断是否达到设定的收敛控制准则。如未达到返回步骤c继续寻优；若达到，则停止寻优，并输出计算结果。

本发明有益效果是：优化的目标是求目标函数的最大值即最大的经济效益，通过调整决策变量的价值系数而实现优化目的，以水土资源定田，使规划结果向所希望的方向倾斜，基于粒子群的智能优化求解采用标准微粒群算法或单纯形微粒群算法，继承了微粒群优化算法原有优点的同时，不用进行复杂的矩阵计算，计算工作量小，优化速度快。

附图说明

图1是本发明的流程框图。

具体实施方式

下面结合附图对本发明的一种具体实施方式做出说明。

图1是本发明的流程框图。基于粒子群算法的小流域水土资源协调优化方法，包括以下步骤：

步骤一：构建水土资源优化配置线性规划数学模型，包括以下步骤:

步骤A:确定决策变量，包括x₁:一类水平梯田；x₂:二类水平梯田；x₃:栽种林木；x₄:汇水的天然林草地；

步骤B:确定目标函数，表达式为maxZ＝c₁x₁+c₂x₂+c₃x₃+c₄x₄；其中，c_i为利润净增值/(hm²·a)(i＝l,2,3,4)；c_i＝收入/(hm²·a)-支出/(hm²·a)。

步骤C:确定约束条件，约束条件方程组为：

式中，e为水平梯田每亩每年的灌溉需水量(m/a)×10⁴m²；Yuxi为鱼鳞坑坡地每亩每年的产流量(m/a)×10⁴m²；caodi为天然林草地每亩每年的产流量(m/a)×10⁴m²；

步骤D:构建完整的数学模型，数学模型表示

其中，步骤D中构建完整的数学模型是以人定田方案的水土资源优化配置数学模型；

以人定田方案数学模型的初始约束条件方程组表示如下:

步骤二：确定水土资源优化方案的空间配置，包括以下步骤：

步骤E:确定空间优化配置水土资源的原则，包括就低原则、就平原则、就水原则、就现状原则；

步骤F:细化具体配置原则，包括配置顺序、坡度要求、地形要求和地貌要求。

步骤三：基于粒子群的智能优化求解，可采用标准微粒群算法或单纯形微粒群算法。

其中，单纯形微粒群算法具体采用单纯形微粒群优化算法(SMPSO算法)，其流程包括以下步骤：

1)确定目标函数；

2)初始化设置SMPSO的微粒群个数p、寻优代数、加速因子、惯性权重系数等参数值，随机初始化微粒群的位置和速度；

3)由具体问题确定微粒群的空间维数，用PSO的随机初始微粒作为单纯形搜索法的初始点进行单纯形搜索，按目标函数评价各微粒的初始适应值，将此时微粒群中各微粒的位置作为P_best，从p_best中找到的具有最好适应值的微粒的位置作为g_best；

4)按照公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，对微粒群p和P_best中各微粒的速度和位置进行更新；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)；

5)将更新后的微粒群p和P_best分别作为单纯形搜索法的初始点进行单纯形搜索，把搜索后的微粒群p和P_best中的每个微粒的适应值分别进行比较，取适应值较优的微粒的位置为新的P_best；

6)比较微粒群P_best和g_best的每个微粒的适应值，取适应值较优的微粒的位置为新的g_best；

7)若满足停止条件，搜索停止，输出全局历史最优位置和全局历史最优适应值为所求结果；否则，返回步骤4继续搜索。

其中，标准微粒群算法的流程包括以下步骤：

a)确定目标函数；

b)设定微粒群微粒个数、寻优代数、加速因子、惯性权重系数参数值，随机初始化各微粒的位置和速度；

c)由所要优化的目标函数的设计变量数确定粒子的空间维数，由目标函数值计算每个粒子的适应值，并计算P_best和g_best；

式中，P_best为微粒目前为止由自己发现的最优位置；g_best为整个群体中所有微粒发现的最优位置；

d)比较每个微粒的适应值与个体极值p_best，如果f(P_i)＞f(P_best)，则P_best＝P_i；

e)比较每个微粒的适应值与全局极值g_best，如果f(P_i)＞f(g_best)，则g_best＝P_i；

f)根据公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，更新每个微粒的飞行速度和空间位置；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)。

g)判断是否达到设定的收敛控制准则。如未达到返回步骤c继续寻优；若达到，则停止寻优，并输出计算结果。

实施例1以承德市的实际情况，基于粒子群的智能优化求解采用标准微粒群算法的实施例

基于粒子群算法的小流域水土资源协调优化方法，包括以下步骤：

步骤一：构建水土资源优化配置线性规划数学模型，首先分别对构建水土资源优化配置线性规划数学模型的三大部分进行分析，然后，再由这三大部分构建数学模型，具体包括以下步骤:

步骤A:确定决策变量，坡面径流聚散利用工程的配置，首先取决于“土”、“水”两大资源，所谓“土”资源是指土地资源，不同坡度的土地，既是产生、汇集坡面径流、承载坡面径流工程的载体，同时也是坡面径流工程储存、输送的水资源所要灌溉的对象，是决定决策变量个数多少、取值大小的重要的基础数据。土地资源的数量取决于规划区的国土面积。

所谓“水”资源是指坡面径流资源，它是坡面径流工程所要储存、输送、利用的对象，是决定决策变量个数多少、取值大小的又一类重要的基础数据。坡面径流资源的数量既取决于规划区的降水量，又取决于规划区可用于汇集坡面径流的国土面积。

决策变量是数学模型待求的未知数，决策变量的解也就是水土资源优化配置的结果，它取决于“水”、“土”两大资源的数量与两大资源的合理搭配，因此，决策变量的个数、每个决策变量的含义及命名应包含“水”、“土”两大资源要素。

决策变量既然是数学模型待求的未知数，那么，在构建数学模型时，是很难准确判断各种决策变量的取舍的,为了防止遗漏，应当把“水”、“土”资源两大要素的各种组合都考虑、包括进来。下面列出数学模型应当包括的四个决策变量。

包括x₁:一类水平梯田；x₂:二类水平梯田；x₃:栽种林木；x₄:汇水的天然林草地；

其中，(l)x₁:一类水平梯田(能灌溉的水平梯田)

一类水平梯田可以充分利用高、中、低位蓄水灌溉，包括自流灌溉、提水灌溉；这类水平梯田的分布高度可以差别较大，只要有对应的汇水区能满足其灌溉的水平梯田均可归入此类；同级位置高度中可以包含多个地块；与各级别位置高度相对应的特定的高、中、低汇水区的供水水源则相应编号为s_i(i＝l,...,n)，同一个s_i可以分别储存在很多个蓄水工程(比如很多个蓄水池)中。

一类水平梯田的灌溉需水总量与汇流区的汇流总量必须均衡。

一类水平梯田的储水、供水原则:同属于x₁的不同图斑地块的水平梯田靠相对应的特定的高、中、低汇水区的蓄水工程供水，而其特定的汇水区的蓄水工程在满足了其相对应的图斑地块梯田自流灌溉用水之后仍有剩余时，可以给高度更低的另外图斑地块的一类水平梯田自流供水。不能自流灌溉的一类水平梯田则需提水灌溉，其提水灌溉的费用在总收益中单独扣除。

(2)x₂:二类水平梯田(水平旱梯田)

二类水平梯田靠自然降水生长，无灌溉之水；因二类水平梯田不灌溉，因此不存在供水汇水区、供水蓄水工程。

(3)x₃:栽种林木(经济林木、果树林木、用材林木等)

栽种林木的鱼鳞坑、坡地可拦截部分坡面径流供栽种的林木用，另外的部分坡面径流则流出鱼鳞坑区汇入蓄水工程供水平梯田灌溉之用，鱼鳞坑内栽种的林木则靠天生长，无灌溉。因x₃不需专门灌溉，但产生部分坡面径流，因此属不完全的汇水区，计算汇水区面积时，按占“完全汇水区”每亩的产流比例取值。

(4)x₄:汇水的天然林草地

汇水的天然林草地、坡地，靠天生长，无蓄水、无灌溉；坡面产流全部汇入蓄水工程。因x₄不灌溉，但产生坡面径流，因此属完全的汇水区。

上述四个决策变量必须全部满足非负要求，即:x_i≥0(i＝l,2,3,4)。

上述四个决策变量中，只有获得正值的决策变量才是被优化选中的工程类型，但“水”、“土”资源条件不同的计算区域，数学模型的求解结果为不同决策变量的组合及决策变量的不同取值，其中的某些决策变量可能会等于零，说明该类型的工程没有被选中。例如，降水量足够大、坡面径流资源丰富的地域，能灌溉的一类水平梯田己能满足人均基本农田的亩数，二类水平梯田已没有存在的必要，其求解结果必然等于零。

步骤B:确定目标函数，目标函数使优化配置后的“收益”最大，用可以量化为货币的“经济效益”为衡量的尺度；目标函数中的价值系数应当是优化配置后单位面积经济效益的增加值——增收的货币值，单位面积每年货币的增收值为价值系数。

水、土资源配置、优化的总目标是争取经济效益最大化。经济效益用目标函数计算。目标函数是由各决策变量构成的一个线性关系式，价值系数作为各决策变量的系数。

价值系数是目标函数最主要的参数，各决策变量的价值系数可以统一按两种取值方法中的一种取值:

(1)按每公顷每年的纯利润取值，即：c_i＝收入/(hm²·a)-支出/(hm²·a)

(2)按改造前、后，单位面积每年的纯利润净增值取值，即c_i＝利润净增值/(hm²·a)(i＝l,2,3,4)；

目标函数表达式为：maxZ＝c₁x₁+c₂x₂+c₃x₃+c₄x₄；其中，c_i为利润净增值/(hm²·a)

确定价值系数时，应全面反映如下因素，并数值化、定量化。

具体地，价值系数按每公顷每年的纯利润取值

应先分别计算其收入项、支出项，计算出c_i值：c_i＝收入/(hm²·a)-支出/(hm²·a)

收入项

收入项:按单位面积(hm²)计算收入项，收入项主要包括作物产值/(hm²·年):按地域的不同分别按每年一季(或一季半或两季)收成计算，每季的收成按同一种作物折算成单位面积产值(单价×产量/(hm²×每季))，然后计算每年每hm²土地的总产值。林木、草，折算成每公顷每年的产值。因为作物的年产值是农户直接受益的，所以可称之为“局部利益”。

支出项

支出项:按单位面积(hm²)计算支出项。

支出项分多年收益的投资、当年收益的投资两种:

一种是为多年服务、由多年分担的工程成本费，包括高位蓄水池、低位蓄水工程(大蓄水池或小水库)、与之配套的扬水机械设备、相对固定的扬水管道、引水(输水)渠、修筑水平梯田、鱼鳞坑等的费用，按使用20年分摊到每公顷每年。

其中，低位蓄水工程的工程费用、与之配套的扬水机械设备、相对固定的扬水管道、引水(输水)渠费用、提水电力费用等，也可以按每亩地每年分摊水费的方式计算。

不同类型的坡地工程，其成本费分为：

①一类水平梯田的工程成本费：根据一类水平梯田的面积配置自流灌溉蓄水池(水池)的蓄水容积及配套引水、输水渠，这些相当于“水费”；除此之外，还有修梯田的工程成本费，相当于“地费”；一类水平梯田坡面径流聚散利用工程的成本费由这两部分组成。

②二类水平梯田的工程成本费:二类水平梯实际上是无水可灌溉的基本农田，当规划区的坡面径流数量有限，只能满足基本农田(人均0.067～0.013hm²)中的一部分农田的灌溉用水，则剩余的无水可灌溉的基本农田只能是旱梯田。因此，二类水平梯田的工程成本费不包括“水费”，只包括修梯田的工程成本费，相当于“地费”，按使用20年分摊到每公顷每年。

优化结果中二类水平梯田的数量受规划区坡面径流产流量的控制，规划区降水量越大，则坡面径流产流量就越多，那么规划区可灌溉的一类水平梯田就越多，剩余的无水可灌的二类水平梯田就越少，当人均基本农田(0.067～0.013hm²)的水浇地满足之后，如果坡面径流的产流量还有富裕，那么二类水平梯田的数量就是零。

③经果林费用:包括修鱼鳞坑的成本费，及栽种经果林的费用，按使用、收益20年分摊到每公顷每年。

④草地的工程成本费：只包括初建草地的工程成本费，相当于“地费”，按使用20年分摊到每公顷每年。

支出项的另一种是为当年服务的成本投入，包括种籽、肥料、农药、提水电力、投工等支出项。

价值系数按改造前、后，单位面积每年的纯利润净增值取值

c_i＝利润净增值/(hm²·a)(i＝l,2,3,4)；确定价值系数的一种简易方法是确定出不同土地类型改造前后单位面积的收益净增加值，直接作为该决策变量的价值系数使用，则目标函数的计算结果就是规划区因工程改造每年所能增加的最大效益值。

步骤C:确定约束条件，约束条件全面反映当地居民拥有、开发利用“水”、“土”资源的状况。地域不同，人均拥有“水”、“土”资源的差别可能会很大，计算的地域有可能会是“水”、“土”资源均很充足、“水”、“土”资源均不充足、“水”足“土”不足或“土”足“水”不足等情况之一。人均拥有的“水”、“土”资源量不同，决定着开发利用的方案不同，其约束条件的构成也就不同。

约束条件必须在规划区坡地坡度分区、坡面径流汇水分区及面积量测的基础上设置。

坡地应按不同的高度等级统计各坡度区的面积，汇水区也应按不同的高度等级划分。但决策变量及约束条件设置不可过于繁琐，可按规划区水资源(坡面径流)总量、每种坡度土地资源的总量来设置和计算，以避免无最优解的情况出现。

在“决策变量”设置的基础上，约束条件一般包含如下几个:

(l)规划区的总土地约束:x₁+x₂+x₃+x₄≤规划区总面积；

(2)坡面径流总量的约束:一类水平梯田的灌溉需水量≤汇水区的坡面产流总量。

坡面径流总量约束条件的设置最为复杂，需要先计算出几个有关的参数，然后再用这些参数写出坡面径流总量的约束条件。需要先计算出的参数是:

①水平梯田每公顷每年的灌溉需水量[m³/(hm²·a)]:

设水平梯田每亩每年的灌溉需水量为e，则:

e＝频率年的灌溉需水量(m/a)×10⁴m²

＝[作物生长的年需水深(m/a)-频率年的有效降水(m/a)]/0.8×10⁴m²；

②天然林草地每公顷每年的产流量[m³/(hm²·a)]:

设天然林草地每亩每年的产流量为caodi，则:

caodi＝频率年的有效降水量(m/a)×频率年的林草地产流率(％)×10⁴m²＝频率年的降水量×0.8×频率年的林草地产流率(％)×10⁴m²；

③鱼鳞坑坡地每公顷每年的产流量[m³/(hm²·a)]:

设鱼鳞坑坡地每亩每年的产流量为Yuxi，则:

Yuxi＝频率年的有效降水量(m/a)×频率年的林草地产流率(％)×鱼鳞坑坡地产流率占主汇水区(天然林草地)产流率的百分比(％)×10⁴m²

＝频率年的降水量(m/a)×0.8×频率年的林草地产流率(％)×鱼鳞坑坡地产流率占主汇水区(天然林草地)产流率的百分比(％)×10⁴m²；

上述参数在构建约束条件之前应计算出，然后构建约束条件为:

e×x₁≤Yuxi×x₃+caodi×x₄

④草地汇水区的土地约束:x₄≥35～90度的坡地总面积

⑤一、二类水平梯田的土地约束:

x₁+x₂≤适于修水平梯田的坡地总面积；适于修水平梯田的坡地为0～20度坡地。

⑥规划方案特殊要求的约束条件:按照规划方案的不同，可以设置特定的约束条件.例如:

当适于修水平梯田的坡地资源足够时，要求一、二类水平梯田按人口数修建的约束，实现“以人定田”方案，则该约束为:x₁+x₂＝按人定地数；

当适于修水平梯田的坡地资源足够、且坡面径流资源也足够时，要求一类水平梯田按人口数修建的约束，实现“以人定水田”方案，则该约束为:x₁＝按人定地数；

充分发挥土地资源修水平梯田的土地约束:x₁+x₂＝适于修水平梯田的坡地总面积。

其中的第⑥条(规划方案特殊要求的约束条件)，往往通过对第⑤条的改造、替代来实现。

总结以上，约束条件方程组为：

步骤D:构建完整的数学模型，可以把上述三部分按照不同的配置原则分别构建出不同的水土资源优化配置线性规划数学模型，分别为“以人数定田”配置方案、“以水土资源定田”等配置方案的水土资源优化配置线性规划数学模型。针对同一个规划区(例如一个小流域)，可以构建其中之一，也可以同时构建多种配置方案供选用。

因为优化的目标是求目标函数的最大值，对我们要解决的问题而言，是求最大的经济效益，所以，各个决策变量价值系数的大小将直接影响着决策变量的求解结果，价值系数大的决策变量会得到模型的“偏爱”而取得较大的求解结果，对规划者来说，可以利用这一特点，通过调整决策变量的价值系数而实现优化目的。其中的“以水土资源定田”的配置方案，可以根据当地梯田与经果林的收益大小、或者规划者的意图来调整价值系数，使规划结果向规划者所希望的方向倾斜。此时的价值系数只能起到按规划者的意图优化水土资源的作用，水土资源优化完成后，目标函数值需代入实际的价值系数来重新计算得到。

数学模型表示

其中，步骤D中构建完整的数学模型是以人定田方案的水土资源优化配置数学模型；

以人定田方案是有前提条件的，即:

(l)当适于修水平梯田的坡地资源足够时，才可以要求一、二类水平梯田按人口数修建x₁+x₂＝按人定地数；

(2)当适于修水平梯田的坡地资源足够、且坡面径流资源也足够时，才可以要求一类水平梯田按人口数修建:x₁＝按人定地数。

以人定田方案数学模型的初始约束条件方程组表示如下:

当适于修水平梯田的坡地资源足够、且坡面径流资源也足够时，约束条件(4)可以直接写成

x₁＝按人定地数

分析数学模型，可以看出它完全具备了线性规划数学模型的三大特点:

①目标函数为线性；

②约束条件全部为线性；

③约束条件线性无关。

步骤二：确定水土资源优化方案的空间配置，以上计算已将一个小流域的土地资源进行了总量上的优化分配，确定了汇水区与灌溉区的面积，从而也分配了小流域的坡面径流资源，但这只是完成了初步的工作，要真正做到空间优化配置水上资源，就必须将规划的结果布置到最优的位置，这就要求首先确定空间优化配置水土资源的原则。

根据坡地坡度不同，将小流域坡地类型划分为三类:<25°的一类坡地，25°～35°的二类坡地，35°～90°的三类坡地。根据地貌位置不同，将其分为坡顶、坡面和沟底。

总的来说，空间优化配置水土资源的原则，要遵循以下基本原则:

(l)就低原则。在其它条件相同情况下，我们将尽可能将梯田、鱼鳞坑修建在低山区。

(2)就平原则。在其它条件相同情况下，我们将尽可能将梯田、鱼鳞坑修建在坡度小的区域。

(3)就水原则。在其它条件相同情况下，我们将尽可能将梯田、鱼鳞坑修建在距汇水区下游出口较近的区域。

(4)就现状原则。经果林的安排尽可能不改变现状，保障经济合理，降低成本。

细化具体配置原则，包括配置顺序、坡度要求、地形要求和地貌要求。

具体地：①配置顺序:优先安排水平梯田，再安排经果林，剩下的地安排天然林草。

②坡度要求。从坡度图上判读:坡改梯田只能在<25°的坡面(一类坡地)进行，经果林种植只能在<35°坡面(一、二类坡地)进行，一类坡地有剩余可并入二类并优先供经果林配置，二类坡地有剩余则并入三类坡地(35°～90°坡面)配置天然林草。

③地形要求。从地形图上判读:当可布置某一使用类型的坡地有富裕时，从低往高布置，剩余部分优先供下一类土地利用形式配置，最后剩下的土地配置天然林草。

④地貌要求。从地貌图上判读:地貌分坡顶，坡面、沟底。坡顶及沟头，只能配置天然林草；坡面考虑坡度限制并按配置顺序安排，一类坡地三种土地利用形式都能配置，二类坡地可用来种经果林或天然林草，三类坡地只能保持种植天然林草；对于沟底，三种土地利用形式都能配置，但依配置顺序优先配置水平梯田。

⑤经济要求。从土地利用现状图上判读:尽量顺应原来土地利用情况进行改造，改变越小花费越少，改造重点在坡面及沟底。

步骤三：基于粒子群的智能优化求解。

步骤三中基于粒子群的智能优化求解采用标准微粒群算法。

微粒群算法(Particle Swarm Optimization，PSO)是一种新兴的演化计算技术，由Kennedy和Eberhart博士受鸟群觅食行为的启发于1995年提出的。PSO研究的最初设想是仿真简单的社会系统，研究并解释复杂的社会行为，但后来发现PSO是解决复杂优化问题的有效技术。PSO的基本思想是通过群体中个体之间的协作和信息共享来寻找最优解。PSO与遗传算法等其它演化算法类似，也是基于群体的，根据对环境的适应度将群体中的个体移动到较好的区域，然而它不像其它演化算法那样对个体使用演化算子，也没有遗传算法复杂的编码、交叉、变异操作和多参数的调整，而是将每个个体看作D维搜索空间中的一个没有体积的微粒，在搜索空间中，微粒以一定的速度和方向飞行，并通过群体间的信息共享和个体自身经验的总结来不断修正个体的行为策略，从而使整体逐渐“飞行”到最佳区域。

在PSO中，D维空间中的每一个“微粒”代表了优化问题的一个可行解，在整个寻优过程中，每个微粒都具有以下几类信息：一个由被优化的函数决定的适应值(fitnessvalue)，一个速度决定其飞行的方向和距离；微粒当前所处位置和到目前为止由自己发现的最优位置(P_best)，这些可视为微粒自身飞行的经验；到目前为止，整个群体中所有微粒发现的最优位置(g_best)，这个可视为微粒群的同伴的经验。这样，各微粒的运动速度受到自身和群体历史运动状态信息的影响，并以自身和群体的历史最优位置来对微粒当前的运动方向和运动速度加以限制和调节，很好地协调了微粒自身运动和群体运动之间的关系。

在具体优化过程中，PSO算法首先初始化生成一个随机微粒群(随机解)，每个微粒在其多维解空间中飞行(寻优)，飞行速度根据其自身的和其它微粒的飞行状况及经验进行动态调整。设在D维空间中有n个微粒，第i个微粒在空间中的位置X_i、速度V_i以及经历的历史最好位置P_i分别表示为：

式中：P_i也称为P_best。每一个微粒都有与优化目标函数f(X)相对应的适应值，一般以优化目标函数作为适应值函数。整个种群中微粒所经历过的具有最好适应值的位置为g_best＝(g₁,g₂,...,g_D)。函数

f(P_best)＝min{f(P₁(t)),f(P₂(t)),...,f(P_n(t))} (4-2)

f(P_best)表示优化目标函数f(X)的全局最优解。对第t代的第i个微粒，PSO算法根据下列进化方程更新第t+1代的第j维的速度和位置：

v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t)) (4-3)

x_ij(t+1)＝x_ij(t)+v_ij(t+1) (4-4)

式中：C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)，j＝1,2,...D。

此外，微粒的速度V_i被一个最大速度V_max所限制。如果当前对微粒的加速导致它在某维的速度v_id超过该维的最大速度v_max,d，则该维的速度被限制为该维最大速度v_max,d。

从运动方程(4-3)中可以发现，微粒飞行的速度由三项组成：第1部分表示微粒维持先前速度的程度，它维持算法拓展搜索空间的能力；第2部分为“认知(cognition)”部分，表示微粒本身的思考，微粒对自身成功经验的肯定和倾向，同时存在适当的随机变动，反映学习的不确定因素；第3部分“社会(social)”部分，表示微粒间的信息共享与合作。

PSO算法迭代终止的条件一般为最大迭代次数或微粒群迄今为止搜索到的最优位置的适应值满足预定的最小适应度阈值。由于所有微粒都根据自身经验和群体经验不断向最优解的方向靠近，当所有微粒趋向同一点时，即认为达到了最优位置。

适应度函数的确定

PSO算法在搜索过程中一般不需要其它外部信息，仅用适应度函数值评估个体的优劣。适应度函数值又称适应度(Fitness)，根据PSO算法原理，微粒的适应值是衡量微粒飞行位置优劣的指标，也是实现微粒由初始位置不断向最佳位置靠近的动力源。在确定适应度函数时，不同问题的适应度函数有不同的意义，但是无外乎以下三种：适应度函数值最小、最大、接近某个值。在基于寿命周期成本的结构优化中，适应度函数在优化中希望LCC值越小越好。

标准微粒群算法(标准PSO算法)的参数分析与设置

微粒群算法自出现至今，已经经历了很多的调整与修正，很多研究人员对参数的选择及其对算法性能的影响进行了大量的分析与实验，为微粒群算法理论和应用的研究奠定了坚实的基础。目前公认的对PSO性能影响较大的参数为：惯性权重w，最大速度V_max，加速常数C₁和C₂。

①惯性权重w

w对于微粒群算法的收敛性起到了很大的作用，w值越大，则全局寻优能力越强，局部寻优能力越弱，反之，则局部寻优能力增强，全局寻优能力减弱。通过调整w的大小来控制以前速度对当前速度的影响，使其成为兼顾全局搜索与局部搜索的一个折中。因为w大，则速度v就大，有利于微粒搜索更大的空间，可能发现新的解域；而w小，则速度v就小，有利于在当前解空间里挖掘更好的解。

在公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))中，如果没有第1部分，即w＝0，则速度只取决于粒子当前位置和其历史最好位置p_best和g_best，速度本身没有记忆性。假设一个粒子位于全局最好位置，它将保持静止。而其他粒子则飞向它本身最好位置p_best和全局最好位置g_best的加权中心。在这种条件下，微粒群将收缩到当前的全局最好位置，更像一个局部算法。在加上第1部分后，微粒有扩展搜索空间的趋势，即第1部分有全局搜索能力。这也使得w的作用为针对不同的搜索问题，调整算法，使全局和局部搜索能力达到平衡。早期的试验将w固定为1.0，但对全局搜索来说，通常的好方法是在前期有较高的探索能力以得到合适的种子,而在后期有较高的开发能力以加快收敛速度。为此可将w设为随时间线性减小,例如由1.4到0，由0.9到0.4,由0.95到0.2等。w按下式进行调节：

式中:MAITER为最大迭代次数。

②最大速度V_max

微粒群算法是通过调整每一次迭代时每一个粒子在每一维上移动的距离来进行的。速度的改变是随机的，但是并不希望不受控制的微粒轨道扩展到越来越广阔的空间，并最终达到无穷。如果微粒要有效地搜索，必须采取某些措施使振幅衰减。传统的方法是使用一个常数V_max，当该维的速度超过该维的最大速度的时候，就限制该维的速度为V_max。参数V_max有利于防止搜索发散(爆炸)，实现人工学习和态度转变。V_max值的选择需对问题有一定的先验知识。为了跳出局部最优需要较大的步长，而在接近最优值时，采用更小的步长会更好。但由于引入了惯性权重w，可消除对V_max的需要,因为它们的作用都是维护全局和局部搜索能力的平衡。这样,当V_max增加时,可通过减小w来达到平衡搜索。而w的减小可使得所需的迭代次数变小。从这个意义上看,可以将V_max固定为每维变量的变化范围,只对w进行调节。

③加速常数C₁和C₂

加速常数C₁、C₂是用来调整微粒的自身经验与社会(群体)经验在其运动中所起作用的权重。在公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))中，如果没有第2部分，即C₁＝0，则微粒没有认知能力，也就是“只有社会(social-only)”的模型。在微粒的相互作用下，有能力到达新的搜索空间。它的收敛速度比标准版本更快，但对复杂问题，则比标准版本更容易陷入局部优值点。如果没有第3部分，即C₂＝0，则微粒之间没有社会信息共享，也就是“只有认知(cognition-only)”的模型。因为个体间没有交互，一个规模为M的群体等价于运行了M个单个粒子的运行，因而得到解的几率非常小。如果C₁＝C₂＝0，则粒子没有任何经验的信息，粒子的运动则显得杂乱无章。Suganthan的实验表明,C₁和C₂为常数时可以得到较好的解。如今的PSO应用中一般都取C₁＝C₂＝2。

目前比较常用的种群初始化方法是随机生成微粒的初始速度和位置。在种群规模选择方面，因为其对具体问题模型的依赖性比较明显，所以还没有任何确切的原则。

标准微粒群算法(标准PSO算法)流程包括以下步骤：

a)确定目标函数；

b)设定微粒群微粒个数、寻优代数、加速因子、惯性权重系数等参数值，随机初始化各微粒的位置和速度；

c)由所要优化的目标函数的设计变量数确定粒子的空间维数，由目标函数值计算每个粒子的适应值，并计算p_best和g_best；

d)比较每个微粒的适应值与个体极值p_best，如果f(P_i)＞f(P_best)，则P_best＝P_i；

e)比较每个微粒的适应值与全局极值g_best，如果f(P_i)＞f(g_best)，则g_best＝P_i；

f)根据公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，更新每个微粒的飞行速度和空间位置；

g)判断是否达到设定的收敛控制准则。如未达到返回步骤c继续寻优；若达到，则停止寻优，并输出计算结果。

标准PSO算法的特点

虽然标准PSO算法的功能与遗传算法非常相似，但是其实现技术却有如下几个显著的优点：

①算法通用,不依赖于问题信息；

②无交叉和变异运算,依靠微粒速度完成搜索；

③有记忆性，微粒和群体的历史最好位置可以记忆并传递给其他微粒；

④需调整的参数较少,结构简单,易于实现；

⑤采用实数编码，直接由问题的解决定，问题解的变量数直接作为微粒的维数；

⑥收敛速度快，在迭代进化中只有最优的微粒把信息传递给其他微粒，属于单向信息流动，在大多数情况下，所有微粒将更快的收敛到最优解；

⑦PSO采用基于邻域的搜索技术，能够利用较小的种群规模保持足够的多样性，从而降低了种群规模。

虽然PSO算法有着明显的优越性，并在求解单目标优化问题中已经取得了一定的成功；但在标准PSO算法中，参数的选择对问题的收敛精度和速度起着重要的作用，往往随着迭代次数的增多，算法的收敛速度降低，全局搜索能力下降，有时无法保证解的全局最优性。

实施例2以承德市的实际情况，基于粒子群的智能优化求解采用单纯形微粒群算法的实施例

实施例2与实施例1的步骤一和步骤二相同，具体参见实施例1，步骤三不同，具体如下：

步骤三：基于粒子群的智能优化求解采用单纯形微粒群算法,具体采用单纯形微粒群算法具体采用单纯形微粒群优化算法。

单纯形微粒群算法

基本的单纯形法是由Spendley等人于1962年提出的一种求解函数最小化的连续优化方法。随后，Nelder和Mead在此基础上发展出非线性单纯形搜索。

非线性单纯形搜索首先在D维欧氏空间中构造一个包含1+D个顶点的凸多面体，由一个随机生成的初始多面体开始，求出各顶点的函数值，并确定其中的最大值、次最大值、最小值，然后通过反射、收缩、延伸等策略求出一个较好解，取代最大值，从而构成新的多面体，这样重复迭代可以逼近一个局部最小值。

该方法直接根据函数值就可以完成优化过程，不需要目标函数的导数信息，是一种广泛应用于非线性无约束优化问题的直接局部搜索技术。其优点是不用求目标函数的一次导数矩阵和赫森矩阵，不用进行复杂的矩阵计算，计算工作量小，优化速度快。缺点是对初始解依赖性较强，容易陷入局部极小点，但是如果给定一个良好的起始点，单纯形法能很快地找到全局极小点。

单纯形微粒群优化算法

PSO是具有随机特征的单目标优化算法，微粒通过跟踪两个最优解来更新自己，与其它随机优化算法相比，不仅能较好地保持多样性，而且收敛速度快，从而给单纯形法(SM)提供了良好的起始点。SM是确定性的单目标优化方法，从一个优良的初始点出发，可以迅速地得到单目标的局部最优解。PSO与SM的结合,从某种意义上是将强大的全局搜索与精确的局部搜索相结合。从PSO所提供的位置出发，在对粒子实施SM搜索的过程中，若SM找到全局最优解，则整个算法成功结束；若SM没有找到全局最优解，它经过一定次数的迭代(预先设定)后，一般来说也会找到更好的位置，然后把由SM产生的新位置放回到PSO算法的微粒群中，形成新一代的微粒群。如此交替使用PSO和SM，直到找到全局最优解，算法成功收敛或达到最大迭代次数。

PSO和SM方法都不需要求导，而只需要计算目标函数和约束函数。将SM嵌入PSO中，不增加PSO的使用限制，使单纯形微粒群优化算法(SMPSO算法)有着广泛的适用性。

单纯形微粒群算法具体采用单纯形微粒群优化算法(SMPSO算法)，其流程包括以下步骤：

1)确定目标函数；

2)初始化设置SMPSO的微粒群个数p、寻优代数、加速因子、惯性权重系数等参数值，随机初始化微粒群的位置和速度；

3)由具体问题确定微粒群的空间维数，用PSO的随机初始微粒作为单纯形搜索法的初始点进行单纯形搜索，按目标函数评价各微粒的初始适应值，将此时微粒群中各微粒的位置作为p_best，从p_best中找到的具有最好适应值的微粒的位置作为g_best；

4)按照公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，更新每个微粒的飞行速度和空间位置；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)；

5)将更新后的微粒群p和p_best分别作为单纯形搜索法的初始点进行单纯形搜索，把搜索后的微粒群p和P_best中的每个微粒的适应值分别进行比较，取适应值较优的微粒的位置为新的P_best；

6)比较微粒群P_best和g_best的每个微粒的适应值，取适应值较优的微粒的位置为新的g_best；

7)若满足停止条件(适应值误差小于设定的适应值误差限度或迭代次数超过最大允许迭代次数)，搜索停止，输出全局历史最优位置和全局历史最优适应值为所求结果；否则，返回步骤4继续搜索。

通过对比实施例1和实施例2，实施例2采用单纯形微粒群算法(Hybrid Algorithmof Sim-plex and Particle Swarm Optimization)，算法在继承了微粒群优化算法原有优点的同时，可以有效增强PSO算法的局部搜索能力，提高PSO的稳定鲁棒性，比标准PSO算法更容易找到全局最优解，是求解高维、复杂问题的一种高效算法。

与现有技术相比，优化的目标是求目标函数的最大值即最大的经济效益，通过调整决策变量的价值系数从而实现优化目的，以水土资源定田，使规划结果向规划者所希望的方向倾斜，基于粒子群的智能优化求解采用单纯形微粒群算法，继承了微粒群优化算法原有优点的同时，不用进行复杂的矩阵计算，计算工作量小，优化速度快。

以上对本发明的一个实例进行了详细说明，但内容仅为本发明的较佳实施例，不能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖范围之内。

Claims (9)

1.基于粒子群算法的小流域水土资源协调优化方法，其特征在于包括以下步骤：

步骤一：构建水土资源优化配置线性规划数学模型；

步骤二：确定水土资源优化方案的空间配置；

步骤三：基于粒子群的智能优化求解。

2.根据权利要求1所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述步骤一中构建水土资源优化配置线性规划数学模型，包括以下步骤:

步骤A:确定决策变量，包括x₁:一类水平梯田；x₂:二类水平梯田；x₃:栽种林木；x₄:汇水的天然林草地；

步骤B:确定目标函数，表达式为maxZ＝c₁x₁+c₂x₂+c₃x₃+c₄x₄；其中c_i为利润净增值/(hm²·a)(i＝l,2,3,4)；

步骤C:确定约束条件，约束条件方程组为：

式中，e为水平梯田每亩每年的灌溉需水量(m/a)×10⁴m²；Yuxi为鱼鳞坑坡地每亩每年的产流量(m/a)×10⁴m²；caodi为天然林草地每亩每年的产流量(m/a)×10⁴m²；

步骤D:构建完整的数学模型，数学模型表示

3.根据权利要求2所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于c_i＝收入/(hm²·a)-支出/(hm²·a)。

4.根据权利要求2所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述步骤D中构建完整的数学模型是以人定田方案的水土资源优化配置数学模型；

以人定田方案数学模型的初始约束条件方程组表示如下:

5.根据权利要求1或2所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述步骤二中确定水土资源优化方案的空间配置，包括以下步骤：

步骤E:确定空间优化配置水土资源的原则，包括就低原则、就平原则、就水原则、就现状原则；

步骤F:细化具体配置原则，包括配置顺序、坡度要求、地形要求和地貌要求。

6.根据权利要求1或2所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述步骤三中基于粒子群的智能优化求解，采用标准微粒群算法或单纯形微粒群算法。

7.根据权利要求6所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述单纯形微粒群算法具体采用单纯形微粒群优化算法。

8.根据权利要求7所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述单纯形微粒群优化算法的流程包括以下步骤：

1)确定目标函数；

2)初始化设置SMPSO的微粒群个数p、寻优代数、加速因子、惯性权重系数参数值，随机初始化微粒群的位置和速度；

3)由具体问题确定微粒群的空间维数，用PSO的随机初始微粒作为单纯形搜索法的初始点进行单纯形搜索，按目标函数评价各微粒的初始适应值，将此时微粒群中各微粒的位置作为P_best，从P_best中找到的具有最好适应值的微粒的位置作为g_best；

4)按照公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，对微粒群p和p_best中各微粒的速度和位置进行更新；

式中，C₁和C₂为加速常数(acceleration constants)，rand₁()和rand₂()为两个在[0,1]范围内变化的随机函数，w为惯性权重(inertia weight)；

5)将更新后的微粒群p和P_best分别作为单纯形搜索法的初始点进行单纯形搜索，把搜索后的微粒群p和P_best中的每个微粒的适应值分别进行比较，取适应值较优的微粒的位置为新的p_best；

6)比较微粒群P_best和g_best的每个微粒的适应值，取适应值较优的微粒的位置为新的g_best；

7)若满足停止条件，搜索停止，输出全局历史最优位置和全局历史最优适应值为所求结果；否则，返回步骤4继续搜索。

9.根据权利要求6所述的基于粒子群算法的小流域水土资源协调优化方法，其特征在于所述标准微粒群算法的流程包括以下步骤：

a)确定目标函数；

b)设定微粒群微粒个数、寻优代数、加速因子、惯性权重系数参数值，随机初始化各微粒的位置和速度；

c)由所要优化的目标函数的设计变量数确定粒子的空间维数，由目标函数值计算每个粒子的适应值，并计算P_best和g_best；

式中，P_best为微粒目前为止由自己发现的最优位置；g_best为整个群体中所有微粒发现的最优位置；

d)比较每个微粒的适应值与个体极值P_best，如果f(P_i)＞f(P_best)，则P_best＝P_i；

e)比较每个微粒的适应值与全局极值g_best，如果f(P_i)＞f(g_best)，则g_best＝P_i；

f)根据公式v_ij(t+1)＝wv_ij(t)+C₁rand₁()(p_ij-x_ij(t))+C₂rand₂()(g_ij-x_ij(t))和x_ij(t+1)＝x_ij(t)+v_ij(t+1)，更新每个微粒的飞行速度和空间位置；

g)判断是否达到设定的收敛控制准则。如未达到返回步骤c继续寻优；若达到，则停止寻优，并输出计算结果。