[0001] 本申请涉及地理信息科学技术领域，尤其涉及一种斑块生成土地利用模拟方法。

[0002] 在土地利用变化模拟中，地理元胞自动机可以模拟土地利用类型的空间分布和演化过程，以便更好地了解土地利用变化的机制和趋势，是研究土地利用变化模拟的主流模型。元胞自动机是一种重要的模拟工具，其模拟了一个由许多小单元(元胞)组成的空间系统，每个元胞的状态会随着时间的推移而发生变化，并受到其周围邻居的影响。

[0003] 元胞自动机的重要特点是通过自下而上的方式操作，即每个元胞只与其周围的有限个元胞进行交互，并通过规则的局部行为来实现全局模式的形成。相比传统的数学模型，元胞自动机能够以更为直观的方式模拟不规则、复杂的现象，并且可以用计算机进行精度无损的建模，以此来模拟各种自然系统的行为。通过元胞自动机，能够更好地了解到复杂系统中个体与局部的行为是如何产生全局有序的模式。

[0004] 地理元胞自动机是元胞自动机的一个特殊形式，它将空间系统分解成规则的、具有一定属性和状态的地理单元，每个单元都可以代表一个特定的土地利用类型。每个元胞周围都有一个邻域，它包括该元胞本身以及与其相邻的其他元胞，邻域的大小和形状可以根据需要进行定义，常见的邻域形状有四邻域、八邻域和十二邻域等。元胞自动机中的元胞状态随时间推移而发生变化，变化的规则由转移规则确定。转移规则是一个函数，它定义了每个元胞的下一个状态如何根据当前状态和邻居状态计算得出。转移规则通常基于元胞的邻域状态，可以是简单的局部函数，也可以是更复杂的全局函数。通过对元胞自动机的状态、邻域和转移规则进行定义和设计，可以对土地利用变化的演化过程进行模拟。

[0005] 目前地理元胞自动机常用两种转换规则挖掘策略：转化分析策略(TAS)和格局分析策略(PAS)。基于TAS的CA模型提取两期土地利用数据之间的所有土地利用的相互转化，但分析所有转化类型会增加CA模型结构的计算复杂度，并且会降低模型的灵活性和普适性。

[0006] 而基于PAS的CA模型只用一期土地利用数据来计算每个单元中土地利用类型出现的概率，并通过土地利用类型之间的竞争来确定未来的元胞的状态，这就避免了在挖掘转化规则的过程中，转化类型组合会随着土地利用类型数量的增加而指数增加的情况。然而基于PAS的CA模型本身就缺乏揭示驱动因素如何引起土地利用变化的能力，因为模型不是基于历史“变化”分析而构建的。因此，基于PAS的CA模型无法获得特定时间间隔的土地利用变化规则，如FLUS模型缺乏揭示驱动因素如何引起土地利用变化的能力。

[0007] 此外，现有模型模拟土地利用斑块演化的能力不足，基于矢量的CA不能用于模拟自然土地利用类型(如森林和草地)的动态变化。而基于栅格的CA模型在模拟斑块演化时难以对多类别和多类型的斑块实现同步演化，如Fore‑SCE模型难以处理存在非常多小斑块的情况，且该模型不是时空动态的模型。

[0031] 为了对本申请的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本申请的具体实施方式。

[0032] 本申请的实施例提供了一种斑块生成土地利用模拟方法。

[0033] 请参考图1，图1是本申请实施例中一种斑块生成土地利用模拟方法的流程图，包括：S1：获取初始土地利用数据；通过初始土地利用数据，得到变化栅格数据；
获取预设个数的影响土地利用和土地覆盖变化的驱动因素，组成驱动因素数据；
S2：在驱动因素数据与变化栅格数据上进行随机点采样，得到采样数据；
S3：对采样数据进行预处理；构建基于随机森林分类算法的随机森林模型；通过预处理后的采样数据，对随机森林模型进行训练；
S4：将所有驱动因素数据输入至训练好的随机森林模型，得到各类用地的发展概率以及每种土地利用类型在区域内的发展潜力图；
具体的，本实例中随机森林的输出为图4所示的7种土地利用类型的发展潜力。

[0034] S5：构建基于多类型随机斑块种子机制的CA模型；设定CA模型的转换矩阵以及邻域权重；根据每种土地利用类型的需求面积，设定CA模型的目标像元个数；S6：将发展潜力图与初始土地利用数据输入CA模型，迭代扫描初始土地利用数据的各个像元，得到各个像元的在邻域内包含的土地利用类型和在邻域内所占的比例以及自适应惯性系数；
通过各个像元在邻域内包含的土地利用类型、在邻域内所占的比例、自适应惯性系数、与各类用地的发展概率，组成得到各个像元上的各类土地利用类型的总体概率；
S7：通过总体概率，构建轮盘；通过轮盘赌方法，使像元区域内各种土地利用类型向竞争获胜的土地利用类型进行转化；
当土地利用类型的邻接效应等于0时，通过计算总体概率，产生不同类型的随机种子；通过随机种子生成新的土地利用类型；
S8：基于阈值递减规则，确定竞争胜出的新的土地利用类型的转化是否发生；
S9：重复步骤S6至步骤S8，直至迭代达到目标像元数目后，停止迭代，输出由土地利用模拟图。

[0035] 具体的，图1中的模型包括基于土地扩张分析策略的转化规则挖掘框架（LEAS）模块和基于多类型随机斑块种子机制的CA模型（CARS）模块。

[0036] 步骤S1包括：S11：获取初始土地利用分类影像；对初始土地利用分类影像进行格式转换或重分类，得到初始土地利用数据；初始土地利用数据包括：第一期土地利用数据以及第二期土地利用数据；
S12：提取初始土地利用数据两个周期之间变化的部份，即从第二期土地利用数据中提取状态发生变化的栅格；
S13：根据第一期土地利用数据以及第二期土地利用数据，提取两个周期之间状态发生变化的栅格，即变化栅格数据。

[0037] 具体的，图2为搜集的初始土地利用变化数据，Landsat图像被分为七种土地利用类型，包括建筑区、耕地、落叶阔叶林、常绿阔叶林、草地、水体和裸地。图3为本实例所使用的多种驱动因素的栅格图，包括：武汉市高程，武汉市坡度，武汉市土壤类型，武汉市年平均降水量，武汉市年平均温度，武汉市邻近开放水域，武汉市高速公路线状矢量，武汉市铁路线状矢量，武汉市交通主干道矢量，武汉市交通次要道路矢量，武汉市交通三级道路矢量，武汉市邻近高铁站，武汉市邻近政府，武汉市GDP数据，武汉市人口数据。

[0038] 步骤S2包括：采用均匀采样策略或比例采样策略，对驱动因素数据与变化栅格数据进行随机点采样，得到采样数据，采样后的样本公式表示为：

[0039] 其中 表示第 个采样点抽取的第 个驱动因素数据； 表示第 个采样点所抽取的变化栅格数据的数据；为转置。

[0040] 具体的，均匀采样策略使得各类土地利用类型的采样点一致，适用于不同土地利用类型面积差异较大的区域；比例采样法使得各类土地利用类型的采样点占总采样点的比例与各类初始土地利用类型与区域面积的比例一致，计算量较小，适合较大尺度区域采样。

[0041] 步骤S3包括：对采样数据进行归一化处理，归一化处理公式表示如下：

[0042] 其中 和 分别表示第 个驱动因素的最大值和最小值；。

[0043] 步骤S3还包括：所述随机森林分类算法，表示如下：

[0044] 其中 是常数，是0或1；是由预设个数的驱动因素组成的向量； 是决策树集的指示函数； 是向量 的第n棵决策树的预测类型； 是决策树的总计数。

[0045] 具体的，本实例使用随机森林算法挖掘城市功能用地类型k的发展变化与驱动因子的映射关系 。 是一个二分类器，所以d的值是0或1。

[0046] 步骤S6包括：总体概率合成公式为：

[0047] 其中 是一个从0到1的随机值； 是为土地利用类型 生成新的土地利用斑块的阈值，它是由用户挖掘的； 为邻域函数；通过邻域函数构建一个扫描窗口，统计扫描窗口内的各类像元的数量来衡量各种土地利用类型在空间上的相互影响；
邻域函数 的定义如下：

[0048] 其中 表示在n×n窗口内最后一次迭代时土地利用类型 占用的网格单元总数， 是不同土地利用类型之间的权重；
自适应函数 的自适应方法定义如下：

[0049] 其中 和 是第 次和第 次迭代时土地利用类型 的当前数量和未来需求之间的差异。

[0050] 阈值递减规则，公式如下所示：

[0051]

[0052] 其中，为递减阈值 的衰减系数，服从正态分布的随机值；为递减阈值； 为衰减次数； 为递减步长； 是定义土地利用类型k是否允许转换为类型 的转化矩阵；表示随机森林概率，即像元 发展为类型 的概率。

[0053] 本申请提供一个实施例如下：设置研究对象为湖北省省会武汉市，其位于湖北省中部的武汉市，市区域总面积约 8494平方公里，是中国中部的一个特大城市。本研究区域中所采用的测试数据为：2003年及2013的武汉市的土地利用数据作为模拟的初始数据，通过决策树对2003年和2013年的相应Landsat影像进行分类。采用武汉市2013年的Landsat影像解译的土地利用变化数据作为验证数据，2003年武汉土地利用变化数据与2013年的数据具有相同的土地利用分类系统。所有的区域数据文件都要被统一成同尺度的栅格影像以方便进行像元尺度上的同质分析。根据土地利用变化模拟的研究经验结合武汉市的历史数据和现有数据，本次实例共选取十五种土地利用/土地覆盖变化的驱动因素。驱动因素分别为：武汉市高程，武汉市坡度，武汉市土壤类型，武汉市年平均降水量，武汉市年平均温度，武汉市邻近开放水域，武汉市高速公路线状矢量，武汉市铁路线状矢量，武汉市交通主干道矢量，武汉市交通次要道路矢量，武汉市交通三级道路矢量，武汉市邻近高铁站，武汉市邻近政府，武汉市GDP数据，武汉市人口数据。

[0054] 本发明中研究对象选择武汉市的原因主要在于：武汉位于长江中游地区，国内生产总值(GDP)在2011年达到913亿美元，在中国所有城市中排名第十三，是一个重要的交通枢纽和经济中心。武汉市地理条件优越，位于长江和汉江的交汇处，拥有丰富的自然资源，如森林、矿产、水资源等。水体占了武汉的很大面积，森林主要分布在丘陵地区，自然植被主要由落叶阔叶树组成。武汉市在过去的十年中，武汉经历了快速的城市发展，导致建成区土地的广泛扩张(从1988年的4.19×104公顷到2011年的49.39×104公顷)，侵占了周围有生态价值的土地(即耕地、林地、草地和水面)。研究武汉市土地利用变化的潜在原因，确定武汉市可持续土地利用方案对其城市发展有重要的意义。研究区域包括武汉市13个区中的11个区，涵盖了中央都市区。

[0055] 第一步：对初始土地利用分类影像进行格式转换或重分类从而获取初始土地利用数据；选取若干影响土地利用/土地覆盖变化的驱动因素，组成驱动力数据；将两期的土地利用数据进行叠加，提取两期数据之间变化的部份，即从第二期土地利用数据中提取状态发生变化的栅格。

[0056] 第二步：在驱动力数据与S1所输出的变化栅格数据上进行随机点采样，根据初始土地利用数据的特点选择采样方式。模型提供两种采样方式：1、均匀采样策略；2、比例采样策略。本实例中使用均匀采样策略使得各类土地利用类型的采样点一致，用采样点提取图2中的土地扩张结果数据和图3中的驱动力数据，采样率设置为0.01。

[0057] 第三步：使用随机森林算法训练采样数据对采样数据进行训练，设定好训练参数与采样率，然后对回归树进行训练。本实例中设置决策树的数量为20；设置随机森林模型的最大特征数为16，等于驱动因素的数量。

[0058] 第4步：将全部驱动力数据输入训练好的随机森林，输出每种土地利用类型在区域内的发展潜力。

[0059] 第5步：将发展潜力图与初始土地利用数据输入多类型随机斑块种子机制的CA模型模块，设定好转换矩阵和邻域权重，以及根据每种土地利用类型的需求面积设定目标像元个数，即未来各类用地需要达到的像元数；本实例中的武汉地区影像像元为100m，因而邻域范围选择3×3摩尔邻域。本实例模拟的时间为：从2003年模拟到2013年，因而目标像元个数与2013年的各类别土地利用分类影像一致，当所有土地利用类型的面积达到其对应的需求面积时模拟自动停止。转换矩阵限制为落叶林、城镇不能转换为其他用地，常绿林除耕地外不能转换为其他用地。每个土地利用类型的邻域权重值可以根据专家知识和一系列模型测试来确定，范围从0‑1。

[0060] 第6步：迭代扫描初始土地利用数据的像元，通过计算该像元在邻域内包含的土地利用类型和在邻域内所占的比例、自适应惯性系数与各类用地发展概率共同合成该像元上各类土地利用类型的总体概率；另外，当像元邻域概率为0时，该类土地利用类型的总体概率由发展概率、自适应惯性系数、随机数、该类土地利用类型的生长阈值组成。

[0061] 本实例中，通过邻域函数构建一个扫描窗口，统计窗口内的各类像元的数量并衡量各个不同土地利用类型在空间上的相互影响，邻域规则的定义如下：

[0062] 其中 表示在n×n窗口内最后一次迭代时土地利用类型k占用的网格单元总数， 是不同土地利用类型之间的权重，因为不同土地利用类型有不同的邻里效应。实例中取n=3，构建摩尔邻域。

[0063] 根据本发明的实施例，在构成轮盘赌的总体概率合成公式为：

[0064] 其中 是一个从0到1的随机值； 是为土地利用类型k生成新的土地利用斑块的阈值，它是由用户挖掘的, 为领域函数， 为自适应函数。本实例设定的生成新种子数量的最大阈值为0.0001，意味着只有0.01%的元胞可以生成新的斑块。 定义如下：

[0065] 其中 和 是第 次和第 次迭代时土地利用类型 的当前数量和未来需求之间的差异。

[0066] 第7步：该像元上的各类土地利用类型的总体概率构成轮盘，通过采用轮盘赌的方法，使区域内各种土地利用类型在像元上竞争，本方法提出了一种阈值递减机制：当一种新的土地利用类型在一轮竞争中胜出，则使用一个递减的阈值来判断转化是否发生。本实例中设定的生成新斑块的递减阈值 为0.5，设定的衰减系数为0.1。

[0067] 第8步:迭代完一幅影像的全部有效像元，返回第六步刷新土地利用变化影像进入下一次迭代，计算到目标像元数目的差值。达到目标像元数目，停止迭代输出结果。

[0068] 图5为最终输出的土地利用模拟图，图6为新的模型和之前旧模型的模拟结果对比，参与对比的模型有：退化的PLUS模型(I)(采用PAS)、退化的PLUS模型(II)(不采用多类型随机斑块种子机制)和FLUS模型(基于人工神经网络的CA模型)。

[0069] 本申请还公开了计算机可读存储介质，该计算机可读存储介质存储有多条指令，指令适于处理器进行加载，以执行上述的一种斑块生成土地利用模拟方法。

[0070] 以上者，仅为本公开的示例性实施例，不能以此限定本公开的范围。即但凡依本公开教导所作的等效变化与修饰，皆仍属本公开涵盖的范围内。本领域技术人员在考虑说明书及实践真理的公开后，将容易想到本公开的其他实施方案。

[0071] 本申请旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未记载的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的范围和精神由权利要求限定。