技术领域
[0001] 本发明涉及一种土地利用自动分区方法，尤指应用于土地利用空间优化配置中的基于粒子群优化算法的土地利用自动分区方法，属于土地利用规划领域。
相关背景技术
[0002] 土地利用分区是土地资源用途管制的主要技术手段，也是土地利用空间优化配置的关键内容。国际上土地利用分区研究开始于20世纪初美国开创的区域土地利用综合分区，随后在20世纪中期开始了以解决各部门之间土地资源合理配置与布局的土地利用分区研究，如韩国、日本开展的国土规划的分级分区体系。近些年来，土地利用分区研究开始注重分区技术手段和分区理念的研究，分区准确性和实用性也得到了一定提高。我国的土地利用分区研究开始较晚，第一个代表性成果是上世纪60年代完成的全国土地利用现状区划，80年代中国科学院地理综合考察委员会在编制全国1∶100万土地资源图过程中提出的全国土地资源潜力区划思想，将全国土地资源分为9个潜力区；之后土地利用分区技术方法研究得到较广泛的关注，提出了基于改进聚类分析、空间叠加分析、星座图法、模比系数法、元胞自动机等土地利用分区方法，这些方法的提出部分克服了传统的土地利用分区以定性为主，分区边界模糊、分区随意性大、分区工作效率低下、分区结果对土地利用指导意义不强的缺陷，并在前两轮土地利用总体规划编制过程中发挥了重要作用。然而由于土地利用问题的复杂性，土地利用系统是一个涉及自然、社会、经济的土地生态经济复合系统，简单的计量地理学方法和地理计算系统能有效分析和模拟土地利用系统的自然属性，而在土地利用的社会经济和人主观意愿等方面的对土地利用的影响则难以得到较好的分析效果，粒子群优化算法也是起源于对简单社会系统的模拟，最初是模拟鸟群觅食的过程。
[0003] 传统土地利用分区方法基于计量地理和地理计算方法难以模拟土地利用系统中人的复杂行为对土地利用系统的影响，导致土地利用分区结果合理性不高，对土地资源合理利用和用途管制指导性不强的缺陷。
具体实施方式
[0022] 本发明所采用的模型流程如图1示。
[0023] 该土地利用自动分区方法包括如下步骤：
[0024] 步骤1.提取土地利用分区的基础数据并进行整合，整合后的数据为土地用途组合及其中的各个图斑信息，提取基础数据是指从其他数据库或其他系统获取土地利用分区所需的数据，如土地利用现状数据、土地利用适宜性评价数据、统计年鉴数据等，整合数据是指这些数据是多种来源、不同格式的数据，要对它进行统一规范化处理，形成一个统一的基础数据库。
[0025] 步骤2.基于粒子群优化算法构建土地利用自动分区模型，以步骤(1)的数据为模型输入数据，以土地利用图斑为数据处理单元，建立待求问题与粒子群之间的映射关系，最终求解土地利用优化布局结果。土地利用图斑抽象为粒子群优化算法中的粒子，每个图斑的位置和土地利用类型对应为粒子的位置(x，y)和种属(i)，土地利用综合效益函数为粒子群适应度函数。
[0026] 步骤3.设置粒子群规模，粒子的最大飞行速度(υxmax，υymax)以及加速权重(c1，c2)，最大迭代次数，误差值，初始化每个粒子的位置、速度和属性；
[0027]
[0028] 式中：υxij和υyij分别为粒子的x和y轴方向上的速度。
[0029] 步骤4.根据适应度函数，计算每个粒子的适应度。
[0030] 粒子的适应度函数如下：
[0031]
[0032]
[0033]
[0034] F＝(Ck，Sk，Zk)
[0035] 满足： xik∈{0，1}
[0036] 式中：n为土地利用图斑总数；cik为第i个图斑的土地用途转变为第k种土地利用类型所需要的费用；sik为第i个土地利用图斑的土地用途k的适宜性评价指数；Ti为第i个土地利用图斑的邻接图斑所构成的集合；nij为土地利用图斑i与图斑j的公共边数；lijh为图斑i与图斑j的第h条公共边的边长；ai为图斑i的土地面积；A1k和A2k分别为第k种土地利用类型所需求的最小和最大面积；若第i个土地单元为第k中土地利用类型时xik＝1，否则xik＝0，Ck为土地变更费用目标；Sk为第k种土地用途的土地适宜性目标；Zk为图斑紧凑性即形状目标；F为多目标综合评价函数。
[0037] 步骤5.计算个体最优值pb和全局最优值pg，并计算出粒子的个体最优位置和全局最优位置，假设f(X)中确定的适应度函数，那么每次更新后的粒子的最佳位置公式为：
[0038]
[0039] 假设粒子群规模为s，全局最优值pg，则：
[0040] Pg(t)∈{P1b(t)，P1b(t)，........，Psb(t)}＝min{f(P1b(t)，f(P2b(t))，.......，f(Psb(t)))}
[0041] 步骤6.计算惯性权重，更新每个粒子的飞行速度，根据步骤5得出个体最优位置和全局最优位置，更新每个粒子的当前位置。
[0042] 惯性权重的函数为：ω(t)＝ωmax-t·(ωmax-ωmin)/Imax
[0043] 式中：t为迭代次数t＝1，2，…，Imax，其中Imax为最大迭代次数。ωmax为最大惯性权重，ωmin为最小惯性权重。
[0044] 飞行速度更新函数为：
[0045]
[0046] 其中：“i”表示粒子i，“j”表示粒子的第j维，w为惯性权重，t表示迭代次数，c1、c2表示加速系数，c1调节粒子自身飞行速度的步长，c1则调节粒子飞向全局最好位置的步长，r1～U(0，1)，r2～U(0，1)为相互独立的随机函数，表示粒子上一时刻速度对当前速度的影响程度，也即影响系数。
[0047] 为了防止粒子飞行速度的无限增大，呈现爆炸无序状态，需增加一组约束条件：
[0048]
[0049] 粒子最大飞行速度υmax限制了粒子的飞行速度υi，υmax决定了粒子在空间内搜索的精度，当取值太大，则粒子容易越过最优解；当υmax取值太小，则粒子飞行速度慢，容易陷入局部搜索空间而无法得到全局最优解。
[0050] 更新粒子的当前位置最近邻法则，若当前位置有多个时，最近邻法则采用面积占优法。该法则针对如何确定模型中不规则性图斑的最邻近关系而提出的，是更新粒子位置时用来判断粒子飞行的下一个位置的准则。当已知粒子的当前位置和速度时，以当前位置对应的图斑中心为圆心，当前速度矢量为半径做圆，则它的下一时刻可能位置应该是与这个圆相交的某个图斑。当存在多个相交图斑时，首先计算每个相交图斑落入圆中的面积占该图斑面积的比例，取比例最大的图斑所在的中心位置作为其下一个位置，如附图2所示，当已知粒子的当前位置和速度时，它的下一个可能位置应该是在以当前位置为圆心(A)，速度矢量为半径的圆上，附图2中所示的虚线圆圈，则它的下一个位置可能是图斑114、111、204或112，根据基于面积占优法的最近邻法则，取这些图斑中落入圆中的面积占图斑面积比例最大的那个图斑所在的中心位置为其下一个位置，即为图中A1位置。
[0051] 步骤7.循环迭代，当满足判断条件时，即前后两次全局最优值之差的绝对值不大于误差值，并且循环迭代的次数不大于最大迭代次数时，搜索结束，否则将更新的值作为步骤4的输入值，继续搜索得到全局最优值和全局最优位置。
[0052] 步骤8.对步骤2)得到的土地利用优化布局，通过对邻近相同土地利用类型的图斑合并处理，提取合并后的图斑界线即得到土地利用分区界线，并最终生成土地利用分区图。
[0053] 实施例1：
[0054] 1.提取某一乡镇土地利用的基础数据并进行整合，整合后的数据为土地用途组合及其中的各个图斑信息，如附图3所示。
[0055] 2.应用粒子群优化算法对上述的数据优化，进行建模和选择最优的土地利用布局。
[0056] 3.设置粒子群规模20个，粒子的最大飞行速度(3，4)以及加速权重(1.44，-61.44)，最大迭代次数为200次，误差值为10 ，初始粒子对应的位置由图斑的重心坐标确定，属性对应十种土地利用类型，范围在(0，10)之间取值，初始速度为最大飞行速度乘以(0，1)间的随机数随机产生的速度。
[0057] 4.根据适应度函数，计算每个粒子的适应度。
[0058] F＝(Ck，Sk，Zk)
[0059] 其中：
[0060]
[0061]
[0062] 5.计算个体最优值pb和全局最优值pg，并计算出粒子的个体最优位置和全局最优位置。
[0063] 6.计算惯性权重，更新每个粒子的飞行速度，根据步骤5得出个体最优位置和全局最优位置，更新每个粒子的当前位置。
[0064] 惯性权重w看作迭代次数的线性函数，设置其从0.98到0.48线性减小，最大惯性权重为0.98，最小惯性权重为0.48，惯性权重的函数为：
[0065]
[0066] 7.循环迭代。当满足判断条件时，即前后两次全局最优值之差的绝对值不大于-610 ，并且循环迭代的次数不大于最大迭代次数200时，搜索结束，否则将更新的值作为步骤4的输入值，继续搜索得到全局最优值和全局最优位置。
[0067] 8.优化结果的处理、分区界线的提取与确定。对完成循环迭代结果的土地利用优化布局图，进行邻近同类图斑的合并处理，提取不同用地类型分区界线，制作土地利用分区结果图，如附图4所示。
