CN103473465B - 基于多目标人工免疫***的土地资源空间布局优化方法 - Google Patents

Abstract

本发明提供了基于多目标人工免疫***的土地资源空间布局优化方法，包括将规划区在基期年的土地适宜性评价结果图采用栅格文件进行存储；编码，将一种土地资源空间布局方案通过整数编码，映射为人工免疫***的抗体，包括用规划区内的一个栅格单元代表实际的一个地块，对应人工抗体中的一个基因位；目标函数设置、约束条件设置、抗体种群初始化，抗体目标值向量计算，然后迭代进行克隆、变异、抗体目标值向量计算、种群更新、终止条件判断，直到终止迭代，并解码得到Pareto优化方案。

Description

基于多目标人工免疫***的土地资源空间布局优化方法

技术领域

本发明属于土地资源空间布局优化配置技术领域，特别是涉及一种基于多目标人工免疫***的土地资源空间布局优化方法。

背景技术

土地资源的空间布局优化是土地利用规划工作的重要组成部分，其本质是在区域规划期土地资源数量结构的约束下，根据土地自身的自然（PH值、土层厚度等）、社会（权属、区位条件）、经济（价格、经济效益等）属性特征，确定每一块土地的用途，以提高区域土地利用的适宜性，使得区域土地资源的空间布局更有利于生产经营管理和生态环境保护。因此，土地资源的空间布局优化问题从本质上看是一类带约束条件的多目标组合优化问题。

相对于传统的多目标组合优化问题，土地资源的空间布局优化配置问题更为复杂，求解难度更大：不仅需要处理传统的数学约束关系，还必须处理复杂的空间约束关系，且决策变量众多（数以万计）。因而，由于技术条件的限制，在早期的土地利用规划中，土地资源的空间布局主要依靠规划专家人为主观的空间上划定各类土地利用分区得以实现。依靠专家对土地利用进行分区的主要问题在于分区结果过于主观、分区方案不精确，使得传统的土地利用规划方案最终无法落实到实际的土地利用中。另一方面，随着我国人口的增长、城市化进程的快速推进，土地资源的不合理利用日益明显，重点表现为城镇建设用地规模扩张过快，农用地（尤其是耕地）流失严重，人地矛盾尖锐，并造成了严重的生态和环境灾难。传统的土地利用规划被视为是一种“强政治，弱技术”的政府行为，已不能满足新时期土地可持续利用和生态环境保护的基本需求。因此，在我国新一轮土地利用总体规划工作的开展和实施过程中，对于村镇级的土地利用规划，必须明确每一块土地的具体用途。由此可见，传统的依靠专家主观确定区域土地用途的方法已不能满足新时期土地利用规划工作的实际需求。

近年来，随着智能优化技术和空间信息技术的快速发展，以及我国土地利用规划工作的强烈需求，一些新理论和新方法开始被广泛应用于土地利用规划工作。在土地资源空间布局优化配置领域，国内外有关研究者开始基于不同的智能优化算法对土地资源的空间布局优化技术进行了探索，并取得了较好的效果，例如基于遗传算法的优化模型有关文献:[1]席一凡，杨茂盛，尚耀华.遗传算法在城市土地功能配置规划中的应用[J].西北建筑工程学院学报(自然科学版)，2001，18(04):190-194.有关文献:[2]张鸿辉，曾永年，刘慧敏.多目标土地利用空间优化配置模型及其应用[J].中南大学学报(自然科学版)，2011，42(04):1056-1065.有关文献:[3]CAO K，YE X.Coarse-grained parallel genetic algorithm applied to a vector based landuse allocation optimization problem:the case study of Tongzhou Newtown，Beijing，China[J].StochEnv Res Risk A，2012，1-10.、基于模拟退火算法的优化模型有关文献:[4]SANTE-RIVEIRA I，BOULLON-MAGAN M，CRECENTE-MASEDA R，et al.Algorithm based on simulated annealing for land-use allocation[J].ComputGeosci，2008，34(3):259-268.有关文献:[5]AERTS J，VAN HERWIJNEN M，STEWART T J.Using simulatedannealing and spatial goal programming for solving a multi site land use allocation problem[M].Berlin;Springer-Verlag Berlin.2003:448-463.、基于微粒群算法的优化模型有关文献:[6]LIU YL，LIU D F，LIU Y F，et al.Rural land use spatial allocation in the semiarid loess hilly area inChina:Using a Particle Swarm Optimization model equipped with multi-objective optimizationtechniques[J].Sci China-Earth Sci，2012，55(7):1166-1177.等。

以上智能优化模型的提出，从技术上极大的推进了土地资源空间布局优化配置相关的工作和研究，是的土地资源空间布局朝着自动化、智能化方向发展。然而，在已有的优化技术中，还存在一定不足，主要表现为：

（1）模型的空间搜索能力较弱

智能优化算法大多通过随机搜索、反复迭代的方式获得最优解，其计算效率随着决策变量的增多而急剧下降，这一现象也被称为“维度灾难”。随着地块数量的增加，现有的土地利用空间布局优化模型的空间搜索能力不足的问题也随之凸显，甚至出现算法不收敛的情况。智能算法存在的上述问题表明，以随机改变地块类型的搜索方式已不能满足大规模土地利用空间布局优化搜索的需求。因此，必须从土地资源空间布局优化这一空间优化问题的基本特性出发，改进模型的空间搜索策略，增强算法的搜索效率，提高算法的收敛速度和计算效率。

（2）模型的多目标处理方法不能满足土地利用多目标决策支持的需要

土地利用优化配置问题是一类带约束条件的多目标优化问题。在土地资源空间优化配置模型研究领域，目前大多数研究都采用了基于线性加权和法的对多目标进行处理。该方法的优势在于计算简单，然而缺陷也是显而易见的，如目标权重难以确定、对Pareto前端的形状非常敏感，不能处理其前端的凹部等。因此，需要引入更先进的多目标处理方法对土地资源空间布局优化模型进行改进。

人工免疫***（Artificial Immune System，AIS）是一类模拟自然免疫***的各种原理和机制而开发的智能优化算法。相对于遗传算法，模拟退火算法、微粒群算法和蚁群算法等智能优化算法，人工免疫***在寻优能力、收敛速度和保持中秋多样性等方面具有其独特的优势，因而被广泛应用于信息安全、机器学习、数据挖掘和智能优化等多个领域。然而，人工免疫***目前在土地资源空间布局优化配置领域还没有得到足够重视，国内外还未见相关的报导。因此，针对已有优化模型的上述不足，本发明将以多目标人工免疫***模型中的NICA算法为基础，针对土地资源空间布局问题的基本特点和求解需求，提出一种空间搜索能力强、基于Pareto占优策略的新型面向土地资源空间布局优化技术，实现土地资源空间布局的自动化和智能化配置。

发明内容

针对现有土地资源空间布局优化技术的局限性，发明一种智能化、自动化的基于多目标人工免疫***的土地资源空间布局优化技术。为土地利用规划提供科学、合理、精确、切实可行的土地资源空间布局优化方案，进而为实现土地资源的合理、高效、可持续利用提供重要保障。

本发明提供的技术方案是一种基于多目标人工免疫***的土地资源空间布局优化方法，包括以下步骤：

步骤一、将规划区土地适宜性评价结果图采用栅格文件进行存储；

步骤二、编码，将一种土地资源空间布局方案通过二维整数编码，映射为人工免疫***的抗体；包括用规划区内的一个栅格单元代表实际的一个地块，对应人工抗体中的一个基因位，每个基因位包含了对应地块的空间位置行列号和地类的特征信息，所述特征信息包括对所有地类的适宜性分值和当前地块配置的用地类型信息，适宜性分值根据土地适宜性评价结果图的栅格文件取得；

步骤三、目标函数设置，包含两个优化目标，分别为最大适宜度S和最大紧凑度Comp，

$S=\mathrm{Max}\left({\mathrm{\Σ}}_{i=1}^N{\mathrm{Suit}}_i\right)---\left(1\right)$

式中，N表示规划区内栅格单元的总数，Suit_i为第i个栅格单元对其当前所配置的地类的适宜性分值；

$\mathrm{Comp}=\mathrm{Min}\left({\mathrm{\Σ}}_{i=1}^n{\mathrm{LSI}}_j\right)---\left(2\right)$

式中，Comp为最大紧凑度目标的取值，n为规划区内斑块的个数，空间上相邻且用地类型相同的栅格单元构成斑块，LSI_j为第j个斑块的形状指数；

步骤四、约束条件设置，包括设定土地资源空间布局方案对应的土地资源数量结构必须等于区域规划设定的土地资源数量结构；

步骤五、抗体种群初始化，设染色体种群规模为AN，结合规划期土地资源数量结构，采取随机的方式产生AN个初始抗体；

步骤六、抗体目标值向量计算，包括采用式（1）和式（2）计算各抗体的目标值向量，并进行Pareto排序，获得Pareto优化解集，作为初始的记忆抗体种群；

步骤七、克隆，包括将当前的记忆抗体根据预设的克隆系数C复制C份，形成新的抗体种群；

步骤八、变异，包括遍历新的抗体种群，对每个抗体实施变异操作；

步骤九、抗体目标值向量计算，采用式（1）和式（2）计算各抗体的目标值向量。

步骤十、种群更新，将抗体新种群和原有的记忆抗体种群进行合并，并对合并后的种群重新进行Pareto排序，获得新的Pareto解集，得到新的记忆抗体种群；

步骤十一、终止条件判断，若算当前迭代次数达到预设的最大迭代次数G，则终止迭代，并解码得到Pareto优化方案。

而且，斑块的形状指数LSI按以下公式得到，

$\mathrm{LSI}=\frac{P}{4\sqrt{A}}---\left(3\right)$

式中，P表示斑块的周长，A表示斑块的面积。

而且，步骤五的实现步骤如下，

步骤5.1，构建列表Y，Y的长度为地类的个数，Y中各元素表示对应地类所需要的栅格单元个数；

步骤5.2，产生空的抗体，抗体各基因位的取值为空；

步骤5.3，遍历各基因位，对于当前基因位，基于地块的土地适宜性分值构建轮盘，采用轮盘赌策略随机确定其地类，轮盘构建的原理如下式，

$p_l=\frac{f_l}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{f}_l}---\left(5\right)$

式中，p_l为地类l被选中的概率，f_l为当前地块对地类l的适宜性分值，L为列表Y中对应地类的总数；随机产生一个取值为0-1的随机数rd，判断rd在轮盘中的位置，进而确定当前地块的地类L_i，并将地类L_i在列表Y中对应的取值减去1，若减1后的元素值为0，则将该元素从列表Y中删除；

步骤5.4，重复步骤5.3，直至列表Y为空，完成初始抗体的生成。

而且，遍历抗体的基因位i，以栅格单元为象元，设基因位总数为N，i的初始值为0；步骤八的实现步骤如下，

步骤8.1，i=i+1，当前基因位的地类为L_i，生成一个取值在0-1之间的随机数rnd，若rnd小于预设的抗体变异率P_m，则进入下一步；否则继续重复本步骤跳转到下一个基因位；

步骤8.2，获得当前基因位的邻域象元集P，并统计P中包含的地类及对应的象元个数；

步骤8.3，若邻域P内只包含一种地类，则返回步骤8.1，否则进入下一步；

步骤8.4，统计P中各地类L_j的出现概率，记为概率P_j；

步骤8.5，获得当前基因位对应地块对所有类型用地的适宜性分值S_j；

步骤8.6，计算各地类的综合选择概率Ｃ_j＝P_j×S_j；

步骤8.7，根据各地类的综合选择概率，采用轮盘赌策略随机选择一个地类L_ｃ作为当前基因位实施变异后的地类；

步骤8.8，在抗体中寻找到一个地类编码为L_c的基因位Ｃ，与当前基因位交换地类编码；

步骤8.9，若i<N，则返回步骤8.1继续执行，否则退出循环，抗体变异操作完成。

而且，步骤8.8包括以下步骤，

步骤8.1.1，在抗体内除当前基因位以外的区域随机生成一个位置y；

步骤8.1.2，若位置y的地类L_y为L_c，则进入下一步，否则返回步骤8.1.1；

步骤8.1.3，获得位置y的邻域象元集P，并统计P内地类的数量，若P内只有一种地类，则返回步骤8.1.1，否则进入下一步；

步骤8.1.4，若y的邻域象元中包含L_i则循环终止，以位置y作为抗体基因交换位，否则返回步骤8.1.1继续查找交换位。

而且，步骤五中，对于新的Pareto解集，统计Pareto解的数量TP，若TP大于预设的最大Pareto解数量MP，则根据各解的拥挤距离大小，删除拥挤距离较小的TP-MP个Pareto解，剩下的Pareto解构成新的记忆抗体种群。

而且，各Pareto解相应抗体的拥挤距离计算方法如下：

1）按照第k个优化目标的目标值对抗体种群进行排序；

2）将排序最前和最末的抗体在目标k方向的距离分量设置为一个正无穷大值；

3）根据抗体在种群中的排序位置，计算各抗体在目标k方向的距离分量d_k，得到第i个抗体在第k个目标上的亲和度距离d_ik；

4）对各抗体在所有目标方向的距离分量值求和，得到各抗体的拥挤距离；

相关计算公式如下，

$d_i=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{d_{\mathrm{ik}}}{D_k}$ （式6）

d_ik＝f_(i+1，k)-f_(i-1，k)

D_k＝δ+max(f_k)-min(f_k)

式中，d_i是第i个抗体的拥挤距离，M是优化目标的数量，d_ik取值范围为0-1；f_(i+1，k)是第i+1个抗体的第k个优化目标的目标值，max(f_k)为抗体种群所有抗体的第k个优化目标中最大值，min(f_k)为抗体种群所有抗体的第k个优化目标中最小值，δ为一个正常数。

本发明的特点：总体上具有精确、科学、自动化、智能化特点，相对于目前在实际工程中广泛应用的各种土地资源空间布局优化方法，本发明主要解决的问题有：（1）充分利用多目标人工免疫算法在优化问题求解方面的优势，将多目标人工免疫算法引入土地资源空间布局优化问题的求解；（2）根据土地资源空间布局优化问题的特点，设计了面向土地资源空间布局优化问题求解的抗体编码方案、优化目标函数、约束条件体系、初始化算子和抗体变异算子。

附图说明

图1为现有技术中多目标人工免疫优化算法（NICA）的基本流程图。

图2为本发明实施例的面向土地资源空间布局优化问题求解的抗体编码原理图。

图3为本发明实施例的土地资源空间布局优化中抗体空间变异的基本原理。

图4为本发明实施例的抗体空间变异算子的基本流程图。

图5为本发明实施例的多目标优化的Pareto前端示意图。

具体实施方式

以下结合附图和实施例说明本发明技术方案。

为便于理解起见，首先解释理论基础：

土地资源的空间布局优化问题从本质上看是一类带约束条件的多目标组合优化问题。多目标人工免疫***作为一种优秀的多目标优化算法，已在许多领域内得到了成功应用。已有研究表明，在多目标优化领域，人工免疫***在算法的收敛速度、保持种群的多样性等多方面较遗传算法、蚁群优化算法等智能优化算法更为出色。在多目标人工免疫***中，NICA算法被证明是一种最为出色的多目标人工免疫优化算法之一。因此，本发明以多目标人工免疫***中的NICA算法为基础，构建面向土地资源空间布局的多目标人工免疫智能优化模型。其中，NICA算法的基本流程见附图1，多目标优化的相关概念如下：

（1）Pareto占优：对于给定的两个决策向量x1和x2，及其对应的目标向量y1和y2，若决策向量y1中的某一目标值得优于y2，且y1的其余目标值不劣于y2的其余目标值，则称x1对x2是Pareto占优的，或者说x1是非支配的。Pareto占优的数学描述可参考相关文献。（式（2）Pareto最优解：对于解集X及其中的一个可行解x1，若不存在任何x2对x1Pareto占优，则称x1为Pareto最优解或非支配解（Non-Dominated）。

（3）Pareto最优解集：所有Pareto最优解所组成的集合称为Pareto最优解集。

附图1中，NICA算法的基本原理如下：

（1）对于任意人工免疫抗体种群，可根据Pareto占优的原则，将其划分为Pareto解集和非Pareto解集。

（2）以随机的方式生成初始种群，并对种群内各抗体的目标值向量进行计算；

（3）依据Pareto占优的原则，划分为Pareto解集和非Pareto解集，即支配解集Ab_{Dom}和非支配解集Ab_{Non}；

（4）克隆（Cloning），将种群的Pareto解集中的各个抗体，复制c倍，形成新的临时种群Ab_{New}；

（5）变异（Mutation），对新种群中的各个抗体进行变异，得到新种群Ab_{New}*；

（6）评价（Evaluation），对变异后的新种群中的各个抗体的目标值向量进行计算；

（7）Pareto排序，将新种群Ab_{New}和原有的Pareto解集合并，并根据各抗体的目标值向量进行Pareto占优比较，得到新的Pareto解集，即非支配解集Ab_{Non}*，并替换原有的Pareto解集。

（8）种群更新，得到非支配解集Ab_{Non}**。在算法迭代早期，由于得到的Pareto解集规模过小，可采取随机的方式生成一定数量的解以期增大Pareto解集的规模；在算法迭代后期，由于Pareto解集规模过大，影响算法执行效率时，则根据抗体之间的拥挤距离，删除过于密集的抗体。有关NICA算法的具体实现，可参阅相关文献，本发明不予赘述。

本发明实施例的基于多目标人工免疫***的土地资源空间布局优化方法，具体实现过程如下：

一、获得规划区在基期年(现状)土地适宜性评价结果图，并采用栅格文件进行存储。其中，土地适宜性评价结果包含了各块土地针对不同适宜性程度，其取值范围是(1-100)，适宜性值越大表明越适宜。例如，某块土地对耕地的适宜性为50，对林地的适宜性为80，则表明该土地更适宜于作为园地使用。栅格大小根据规划区的范围大小进行设置，通常可取10米10米，20米×20米，25米×25米，50米×50米，100米×100米等几种规格。当规划区面积较大（县级）时，可采用50米×50米，100米×100米或更大规格的栅格，当规划区面积较小时（村、镇级）用10米10米，20米×20米，25米×25米大小的栅格。适宜性分值根据土地适宜性评价工作结果得到，可采用现有技术实现或者直接导入预先评价的结果，本发明不做赘述。

二、编码，将一种土地资源空间布局方案通过整数编码，映射为人工免疫***的人工抗体，并在计算机中以一个大小M*N的二维数组进行存储。其中，M表示行数，N表示列数。规划区内的一个栅格代表实际的一个地块，对应人工抗体中的一个基因，对应二维数组中的一个数组元素。其中，一个基因包含了对应地块的空间位置（行列号）和地类的特征信息（包括对所有地类的适宜性分值和当前地块配置的用地类型信息），具体编码原理见附图2；因此，本发明中提到的二维数组中的数组元素、地块、栅格单元、象元和基因所指均为同一概念。

图2中，土地利用类型采取整数编码的方式，将耕地、园地、林地、草地、城镇建设用地和农村居民点用地分别对应为整数1，2，3，4，5，6。初始状态的用地类型信息为空，后续可经过抗体种群初始化进行配置。空间上的二维网格划分及其相应的地类编码构成了一个土地利用空间布局方案，其中的编码是单元格的某一个状态，有可能是土地利用现状，有可能是随机确定的，也有可能是根据适宜性分值的大小确定的。；任意土地利用空间布局方案均可采用上述方式进行描述，并以二维数组的方式进行存储。二维数组中数组元素的行号（Row）和列号（Col）确定了该数组元素对应在空间上的位置；A1，A2等属性字段用于记录该地块对应每一种用地类型的适宜性分值，从土地适宜性评价结果图相应栅格中提取。例如，A1可用于存储该地块对应耕地的适宜性分值，A2用于存储该地块对应园地的适宜性分值，依次类推。

三、目标函数设置，土地资源空间布局优化的目标函数主要根据景观生态学的相关原理进行设计，主要包含两个优化目标:

（1）最大适宜度（Maxmazation of the overall land suitability）

土地资源空间布局的首要原则是要因地制宜的安排各类用地的空间分布。在进行土地资源优化空间布局优化时，必须尽可能提高区域内土地利用的适宜性，做到“地尽其用”.其数学描述为：

$S=\mathrm{Max}\left({\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{Suit}}_i\right)---\left(1\right)$

式中，N表示规划区内栅格单元的总数，Suit_i为第i个栅格单元对其当前所配置的地类的适宜性分值。

（2）最大紧凑度（Maxmazation of the spatial compatness）

根据景观生态学原理，土地资源空间布局方案中，空间上相邻、且用地类型相同的栅格象元构成了斑块。斑块的大小、形状对区域生态***功能保护、物质和能量流动产生重要影响。例如，对于耕地、园地、林地等农业或生态用地的安排应当尽可能形成集中连片的大斑块，空间上实现一定程度的集聚以方便管理、节约成本，同时也有利于农用地保护；而对于农村居民点用地则要体现“集中与分散相结合”的原则，一方面一定程度适当的集中由于农村基础设施的建设，而另一方面又需要一定程度的分散，以减少耕作距离和生产经营管理成本。

土地资源空间布局优化中的“最大紧凑度”目标的设立就是为了对空间上各类土地的图斑大小和形状进行优化，使之能更有利于农用地保护、降低生产经营成本，提高生产便利度等，其数学描述为：

$\mathrm{Comp}=\mathrm{Min}\left({\mathrm{\&Sigma;}}_{i=1}^n{\mathrm{LSI}}_j\right)---\left(2\right)$

式中，Comp为最大紧凑度目标的取值，n为规划区内土地利用簇（斑块）的个数。Comp的理论最小值为n。Comp的取值越小，则表明空间布局方案中各斑块的形状越优。LSI_j为第j个斑块的形状指数，其数学描述为：

$\mathrm{LSI}=\frac{P}{4\sqrt{A}}---\left(3\right)$

式中，P表示斑块的周长，A表示斑块的面积。

四、约束条件设置，土地资源空间布局优化方案对应的土地资源数量结构必须等于区域规划设定的土地资源数量结构，其数学描述为：

${\mathrm{\&Sigma;}}_{t=1}^N{a}_t\&times;x_{\mathrm{tl}}=A_l,x_{\mathrm{tl}}\&Element;\left\{\mathrm{0,1}\right\}---\left(4\right)$

式中，A_l为区域内规划年土地资源数量结构中规定的第l类用地的面积，N表示规划区内栅格单元的个数，a_t为第t个栅格单元所代表地块的真实面积。x_tl为一个决策变量，当栅格单元的地类为l时取值为1，否则取值为0。

五、抗体种群初始化，设染色体种群规模为AN，结合规划期土地资源数量结构，采取随机的方式产生AN个初始抗体，具体步骤为：

（1）构建列表Y，Y的长度为区域土地利用类型的个数，Y中各元素表示对应地类所需要的栅格个数。

（2）产生空的抗体，抗体各基因的取值为空。

（3）第i个栅格单元可标记为地块i，即基因位i。遍历各基因位，对于当前基因位i，基于象元的土地适宜性分值构建轮盘，采用轮盘赌策略随机确定其地类。轮盘构建的原理如下

$p_l=\frac{f_l}{\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{f}_l}---\left(5\right)$

式中：p_l为用地类型l被选中的概率，f_l为当前地块对用地类型l的适宜性分值，L为列表Y中对应地类的总数。地块i对某类用地的适宜性分值越高，则该用地类型被选中作为当前地块的用地类型的概率越大。随机产生一个（0-1）的随机数rd，判断rd在轮盘中的位置，进而确定当前地块i的地类L_i，并将地类L_i在列表Y中对应的取值减去1，若减1后的数组元素值为0，则将其从列表Y中删除。

（4）重复步骤（3）直至列表Y为空，即可完成初始抗体的生成。

六、抗体目标值向量计算，采用式（1）和式（2）计算各抗体的目标值向量，进行Pareto占优排序，获得Pareto优化解集，作为记忆抗体种群。具体实现为现有技术，本发明不予赘述

七、克隆，将当前的记忆抗体根据克隆系数C，复制C份，形成新的抗体种群；

八、变异，遍历新的抗体种群，对每个抗体实施变异操作，基本原理见附图3，算子流程见附图4。图中，首先以随机的方式确定要进行变异的地块（基因位）（如图中的①），即作为变异位，进而分析当前地块周边地块土地利用情况（即邻域单元状态）及紧凑度、适宜性信息，确定当前地块要转变的目标用地类型；根据目标用地类型以随机的方式在空间上搜索用地类型为当前地块要转变的目标用地类型的地块（图4中的②）；最终两个地块互相交换用地类型（图4中的③，类型6和2交换）。针对土地资源空间布局优化问题规模巨大，传统随机优化算法收敛过慢的问题，本发明对人工免疫***的变异策略进行改进，加入领域知识的指导，算子流程见附图4。遍历抗体基因位i，设基因位总数（象元总数）为N，i的初始值为0。

对照附图4的流程图，抗体基因变异的具体步骤如下：

（1）i=i+1，当前基因位的地类为L_i。生成一个（0-1）之间的随机数rnd，若rnd小于抗体的变异率P_m，则进入下一步。否则跳转到下一个基因位，继续重复本步骤。其中P_m作为算法参数通过用户设置，通常取1/N。

（2）获得当前基因位的邻域象元集P，并统计P中包含的地类及其对应的象元个数。

（3）若邻域P内只包含一种地类，则返回步骤（1），否则进入下一步。

（4）统计P中各地类的出现概率，即各地类L_j在邻域P中所占的百分比，记为概率P_j。

（5）获得当前基因位对应地块对所有类型用地的适宜性分值S_j。

（6）计算各地类的综合选择概率Ｃ_j＝P_j×S_j。

（7）根据各地类的综合选择概率，采用轮盘赌策略（即赌轮法）随机选择一个地类L_ｃ作为当前基因位实施变异后的地类。

（8）在抗体中寻找到一个地类编码为L_c的基因位Ｃ，作为基因交换位，并与当前基因位交换地类编码，即交换位置i和位置Ｃ的地类。

（9）若i<N，则返回（1）继续执行，否则退出循环，抗体变异操作完成。

上述步骤中，步骤（8）在抗体内寻找一个地类编码为L_c并与当前基因位交换地类编码的基因位的算法又可进一步分为以下几步：

（1）在抗体基因内除当前基因位以外的区域随机生成一个位置y。

（2）若位置y的地类编码L_y为L_c，则进入下一步，否则返回步骤（1）。

（3）获得位置y的邻域象元集P，并统计P内地类的数量。若P内只有一种地类，则返回步骤（1），否则进入下一步。

（4）若位置y的邻域象元集P中包含L_i则循环终止，以位置y作为抗体基因交换位。否则返回步骤（1）继续查找交换位。

九、抗体目标值向量计算，采用式（1）和式（2）计算各抗体的目标值向量。

十、种群更新，将抗体新种群和上一次迭代周期内获得的记忆抗体种群进行合并，并对合并后的种群重新进行Pareto排序，获得新的Pareto解集，得到新的记忆抗体种群，下一迭代时采用此次生成的新的记忆抗体种群。

为限制最大Pareto解数量，在进行Pareto排序获得新的Pareto解集后，本发明进一步提出：统计Pareto解集中Pareto解的数量TP，若TP大于预设的最大Pareto解数量MP，则根据各解的拥挤距离大小，删除拥挤距离较小的TP-MP个Pareto解，剩下的Pareto解构成新的记忆抗体种群。其中，抗体拥挤距离的计算方法如下：

1）按照第k个优化目标的目标值对抗体种群进行排序；

2）将排序最前和最末的抗体在目标k方向的距离分量设置为一个正无穷大值；

3）根据抗体在种群中的排序位置，计算各抗体在目标k方向的距离分量d_k，得到第i个抗体在第k个目标上的亲和度距离d_ik；

4）对各抗体在所有目标方向的距离分量值求和，得到各抗体的拥挤距离。

计算公式如下：

$d_i=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\frac{d_{\mathrm{ik}}}{D_k}$ （式6）

d_ik＝f_(i+1，k)-f_(i-1，k)

D_k＝δ+max(f_k)-min(f_k)

式中，d_i是第i个抗体的拥挤距离，M是优化目标的数量，d_ik的取值范围为（0-1）；f_(i-1，k)是第i+1个抗体的第k个优化目标的目标值，max(f_k)为抗体种群所有抗体的第k个优化目标中最大值，min(f_k)为抗体种群所有抗体的第k个优化目标中最小值。为了防止第k个目标的距离区间D_k（第为0，设置δ为一个非常小的正常数。本发明中有两个优化目标，因此M=2，k的取值为1或2。

十一、终止条件判断，若算法当前迭代次数达到用户预设的最大迭代次数G，则程序终止，并解码得到Pareto优化方案，否则返回步骤七。由规划决策者根据实际需要，从中选择一个或多个作为最终的规划备选方案。

参数设置建议：算法的关键参数设置是否合理也决定了算法的性能和运行结果的准确度，涉及参数及建议的参数设置如下：

参数	参数描述	建议取值
			AN	初始种群规模	建议取值50-200
C	克隆倍数	建议取值3-6
			MP	最大pareto解集规模	建议取值，50-200
Pm	变异概率	建议取值1/N，N为象元总数
			G	迭代次数	建议取值范围(1～10)*N

为说明效果起见，选取某县土地资源空间布局优化问题作为案例，需要对其耕地、园地、林地、草地、城镇建设用地和农村居民点用地的空间布局进行智能优化。其中，栅格象元大小为50×50米，优化配置单元为86万个，在计算机中编写程序实现本发明设计的优化方法，设置初始种群规模为100，最大非支配解集规模为64，克隆倍数为4，变异率为90万分之一，进化代数为90万代。执行优化后，得到Pareto前端见图5。其中，图5中，三角形的点为该区域土地资源空间布局现状对应的适宜性值和紧凑度值。由图5可知，优化后的Pareto解集的中各Pareto解的适宜性目标和紧凑度目标均优于现状。在实际工作中，可综合考虑优化方案的适宜性目标和紧凑度目标，选择方案A作为最终的规划方案

以上仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围。熟悉本领域的技术人员在不违背本发明所指示的精神与原理下所完成的一切等效变形、替换或修饰，仍包含在本发明权利要求所限定的范围内。

Claims (6)

1.一种基于多目标人工免疫***的土地资源空间布局优化方法，其特征是，包括以下步骤：

步骤一、将规划区土地适宜性评价结果图采用栅格文件进行存储；

步骤二、编码，将一种土地资源空间布局方案通过二维整数编码，映射为人工免疫***的抗体；包括用规划区内的一个栅格单元代表实际的一个地块，对应人工抗体中的一个基因位，每个基因位包含了对应地块的空间位置行列号和地类的特征信息，所述特征信息包括对所有地类的适宜性分值和当前地块配置的用地类型信息，适宜性分值根据土地适宜性评价结果图的栅格文件取得；

步骤三、目标函数设置，包含两个优化目标，分别为最大适宜度S和最大紧凑度Comp，

$S=Max\left({\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{Suit}}_i\right)---\left(1\right)$

式中，N表示规划区内栅格单元的总数，Suit_i为第i个栅格单元对其当前所配置的地类的适宜性分值；

$Comp=Min\left({\mathrm{\&Sigma;}}_{j=1}^n{\mathrm{LSI}}_j\right)---\left(2\right)$

式中，Comp为最大紧凑度目标的取值，n为规划区内斑块的个数，空间上相邻且用地类型相同的栅格单元构成斑块，LSI_j为第j个斑块的形状指数；

步骤四、约束条件设置，包括设定土地资源空间布局方案对应的土地资源数量结构必须等于区域规划设定的土地资源数量结构；

步骤五、抗体种群初始化，设染色体种群规模为AN，结合规划期土地资源数量结构，采取随机的方式产生AN个初始抗体；步骤五的实现步骤如下，

步骤5.1，构建列表Y，Y的长度为地类的个数，Y中各元素表示对应地类所需要的栅格单元个数；

步骤5.2，产生空的抗体，抗体各基因位的取值为空；

步骤5.3，遍历各基因位，对于当前基因位，基于地块的土地适宜性分值构建轮盘，采用轮盘赌策略随机确定其地类，轮盘构建的原理如下式，

$p_l=\frac{f_l}{\underset{l=1}{\overset{L}{\&Sigma;}}{f}_l}---\left(5\right)$

式中，p_l为地类l被选中的概率，f_l为当前地块对地类l的适宜性分值，L为列表Y中对应地类的总数；随机产生一个取值为0-1的随机数rd，判断rd在轮盘中的位置，进而确定当前地块的地类L_i，并将地类L_i在列表Y中对应的取值减去1，若减1后的元素值为0，则将该元素从列表Y中删除；

步骤5.4，重复步骤5.3，直至列表Y为空，完成初始抗体的生成；

步骤六、抗体目标值向量计算，包括采用式(1)和式(2)计算各抗体的目标值向量，并进行Pareto排序，获得Pareto优化解集，作为初始的记忆抗体种群；

步骤七、克隆，包括将当前的记忆抗体根据预设的克隆系数C复制C份，形成新的抗体种群；

步骤八、变异，包括遍历新的抗体种群，对每个抗体实施变异操作；

步骤九、抗体目标值向量计算，采用式(1)和式(2)计算各抗体的目标值向量；

步骤十、种群更新，将抗体新种群和原有的记忆抗体种群进行合并，并对合并后的种群重新进行Pareto排序，获得新的Pareto解集，得到新的记忆抗体种群；

步骤十一、终止条件判断，若当前迭代次数达到预设的最大迭代次数G，则终止迭代，并解码得到Pareto优化方案。

2.根据权利要求1所述基于多目标人工免疫***的土地资源空间布局优化方法，其特征是：斑块的形状指数LSI按以下公式得到，

$LSI=\frac{P}{4\sqrt{A}}---\left(3\right)$

式中，P表示斑块的周长，A表示斑块的面积。

3.根据权利要求1所述基于多目标人工免疫***的土地资源空间布局优化方法，其特征是：设遍历抗体的基因位i，以栅格单元为象元，设基因位总数为Num，i的初始值为0；步骤八的实现步骤如下，

步骤8.1，i＝i+1，当前基因位的地类为L_i，生成一个取值在0-1之间的随机数rnd，若rnd小于预设的抗体变异率P_m，则进入下一步；否则继续重复本步骤跳转到下一个基因位；

步骤8.2，获得当前基因位的邻域象元集P，并统计P中包含的地类及对应的象元个数；

步骤8.3，若邻域P内只包含一种地类，则返回步骤8.1，否则进入下一步；

步骤8.4，统计P中各地类L_j的出现概率，记为概率P_j；

步骤8.5，获得当前基因位对应地块对所有类型用地的适宜性分值S_j；

步骤8.6，计算各地类的综合选择概率C_j＝P_j×S_j；

步骤8.7，根据各地类的综合选择概率，采用轮盘赌策略随机选择一个地类L_c作为当前基因位实施变异后的地类；

步骤8.8，在抗体中寻找到一个地类编码为L_c的基因位C，与当前基因位交换地类编码；

步骤8.9，若i<Num，则返回步骤8.1继续执行，否则退出循环，抗体变异操作完成。

4.根据权利要求3所述基于多目标人工免疫***的土地资源空间布局优化方法，其特征是：

步骤8.8包括以下步骤，

步骤8.8.1，在抗体内除当前基因位以外的区域随机生成一个位置y；

步骤8.8.2，若位置y的地类L_y为L_c，则进入下一步，否则返回步骤8.8.1；

步骤8.8.3，获得位置y的邻域象元集P，并统计P内地类的数量，若P内只有一种地类，则返回步骤8.8.1，否则进入下一步；

步骤8.8.4，若y的邻域象元中包含L_i则循环终止，以位置y作为抗体基因交换位，否则返回步骤8.8.1继续查找交换位。

5.根据权利要求1或2或3或4所述基于多目标人工免疫***的土地资源空间布局优化方法，其特征是：步骤十中，对于新的Pareto解集，统计Pareto解的数量TP，若TP大于预设的最大Pareto解数量MP，则根据各解的拥挤距离大小，删除拥挤距离较小的TP-MP个Pareto解，剩下的Pareto解构成新的记忆抗体种群。

6.根据权利要求5所述基于多目标人工免疫***的土地资源空间布局优化方法，其特征是：各Pareto解相应抗体的拥挤距离计算方法如下：

1)按照第k个优化目标的目标值对抗体种群进行排序；

2)将排序最前和最末的抗体在目标k方向的距离分量设置为一个正无穷大值；

3)根据抗体在种群中的排序位置，计算各抗体在目标k方向的距离分量d_k，得到第i个抗体在第k个目标上的亲和度距离d_ik；

4)对各抗体在所有目标方向的距离分量值求和，得到各抗体的拥挤距离；

相关计算公式如下，

d_ik＝f_(i+1，k)-f_(i-1，k)

D_k＝δ+max(f_k)-min(f_k)

式中，d_i是第i个抗体的拥挤距离，M是优化目标的数量，d_ik取值范围为0-1；f_(i+1，k)是第i+1个抗体的第k个优化目标的目标值，max(f_k)为抗体种群所有抗体的第k个优化目标中最大值，min(f_k)为抗体种群所有抗体的第k个优化目标中最小值，δ为一个正常数。