基于MCR-ANN-CA模型的包头市生态用地演变模拟

刘建华 张启斌 YANG Di 岳德鹏 于 强 杨 斓

(1.北京林业大学精准林业北京市重点实验室， 北京 100083； 2.弗罗里达大学地理系， 盖恩斯维尔 FL 32611)

0 引言

《全国生态环境保护纲要(国发[2000]38号)》中提到了生态用地一词，明确指出了生态用地的重要性[1]。众多学者基于此不断发展生态用地的概念和内涵。一种观点认为应该从土地类型上来划分生态用地，即提供自然生态系统服务价值的土地都可以被视为生态用地。另一种观点认为应该从土地的主体功能来区分，主张农业用地不应作为生态用地[2-3]。根据前人理论成果，结合包头市实地情况，本文将林地、草地、水体作为生态用地范围。

生态用地作为自然生态系统服务的基本载体，是解决城市建设用地扩张与生态保护矛盾的综合途径[4]。在我国西北干旱地区，生态用地在防止土地沙漠化、水土保持等方面的生态功能更加突出，对于维持区域生态系统健康稳定具有至关重要的作用[5-6]。对生态用地的空间变化进行模拟，可以为干旱地区生态用地科学规划和生态安全格局构建提供重要参考价值。

生态用地的演变模拟是土地利用演变模拟的一部分，当前土地利用格局的演变模拟主要包括两种类型，一种是宏观回顾模型，如逻辑回归(Logistic regression)模型[7]、马尔可夫(Markov)模型[8]、系统动力学(System dynamic)模型等[9-10]，这些模型根据土地利用的时空数据，量化土地利用在宏观时空尺度上的变化规律，对未来情景进行预测。另一种是微观预测模型，包括多智能体模型(Multi-agent)[11]、元胞自动机(Cellular automata, CA)[12]等，这种模型从微观尺度入手，通过研究土地利用与生态、人类活动等因子的相互作用，预测宏观上土地利用的未来演变。将两种模型结合，可得到二者的混合模型，既可以考虑全局宏观效应，又可以考虑多个因子对土地利用的微观作用机制，相比单一模型具有较大优势。将宏观回顾模型与元胞自动机结合，是这种模型的一种常见形式，如CA-Markov模型[13]。在这种模型中，有两个关键点对其模拟精度有较大影响，一是元胞邻域内不同土地利用类型间转变规则，即邻域规则；二是地形、生态环境等因子对土地利用演变影响，即适宜性规则[14]。当前传统元胞自动机邻域规则制定方法较多，如支持向量机、多目标决策、Markov模型等[15]，用户往往难以选择合适的规则定义方法，部分线性模型也难以根据邻域特征精确模拟生态用地演化等非线性过程；基于逻辑回归、多标准评价模型及层次分析法等制定的适宜性规则虽然考虑了多种因素对土地利用变化的影响，但多基于单个像元利用空间叠加手段独立评价其适宜性，不能根据生态用地演化的过程量化研究区内各景观单元演化为生态用地的适宜性。

针对上述问题，选取西北干旱地区典型城市包头市为研究区，基于源汇理论，利用最小累积耗费阻力(Minimum cumulative resistance, MCR)模型量化生态用地演化过程中，从“生态源地”到其他土地利用类型的适宜性，同时利用人工神经网络(Artificial neural network, ANN)提取CA邻域转换规则，构建MCR-ANN-CA模型对包头市生态用地的演变过程进行模拟，以期为区域生态用地规划及生态建设提供理论与方法支持。

1 材料与方法

1.1 研究区概况

包头市位于内蒙古自治区西部(东经109°13′～111°26′，北纬40°13′～42°44′)，面积27 768 km2。包头市深处内陆，气候为典型的温带干旱、半干旱大陆性气候，冬季寒冷干燥、夏季炎热多雨[16]，年平均气温2.0～7.7℃，年均降水量175～400 mm，年均蒸发量为2 100～2 700 mm。包头市可利用地表水总量为9×108m3，地下水补给量为8.6×109m3。黄河流经包头境内214 km，水面宽130～458 m，最大流量6 400 m3/s，年平均径流量为2.6×1011m3，是包头市主要用水来源[17]。包头市海拔976～2 317 m，地势中间高南北低，北部丘陵、中部山地、南部平原分别占总面积的14.49%、75.51%和10%[18]。干旱的气候条件与起伏的地貌特征使得包头市生态环境较为脆弱，面临较高的水土流失与土地沙漠化、荒漠化风险。近年来，包头市城市建设规模扩张迅速，房地产开发项目、工业园区等建设项目不断涌现，导致草地、林地、湿地等生态用地遭到占用与破坏，生态风险有所升高[19]。

1.2 数据来源与处理

本研究主要数据源包括包头市土地利用数据(2006、2011、2016年)、归一化植被指数(Normalized difference vegetation index, NDVI)、数字地面高程模型(Digital elevation model, DEM)、各区县旗人口密度及工业园区分布数据。其中土地利用数据由Landsat-8遥感数据解译得到，工业园区分布数据来自《包头市土地利用总体规划》和《包头市城市总体规划》；归一化植被指数(NDVI)、数字地面高程模型(DEM)来自地理空间数据云(http:∥www.gscloud.cn)；各区人口密度数据来自包头市统计年鉴。

根据研究需要，从遥感影像中提取耕地、林地、草地、建设用地、水体、未利用地并通过实地验证确保其精度，其中林地、草地、水体为本文中生态用地范围；利用DEM数据提取高程、坡度数据。

1.3 最小累积耗费阻力模型

MCR模型由KANPPEN提出，最早应用于物种迁徙过程研究，之后在物种保护、景观格局分析等方面取得了广泛应用[20]。该模型主要考虑3个因素，即“源”、阻力和累积代价，通过对3个因素的分析，对“源”克服阻力向外传播所耗费的代价或者所做的功进行描述[21]。MCR模型的一般形式为

R=f∑DijRi

(1)

式中R——最小累积阻力

f——未知负函数，表示最小累积阻力与生态适宜性的负相关关系

Dij——从源j到景观单元i的空间距离

Ri——景观单元i处的阻力

本文将生态用地的演化过程看作生态用地对其他景观的竞争性控制过程，且这种演化必须通过克服阻力实现，这样生态用地的演化过程就可以抽象为从源(现有生态用地斑块)到汇(其他景观单元)克服阻力做功的水平过程[22]。由于区域下垫面差异，不同空间位置的土地演化为生态用地的阻力是不同的，通过由“源”到当前像元的累积阻力可量化当前像元演化为生态用地的概率，即阻力越大，该像元演化为生态用地的概率越小，演化为其他用地类型的概率越大。本文利用MCR模型，综合考虑土地利用、NDVI、坡度、政府规划工业园区、人口密度、水体距离、高程因子构建累积耗费阻力面，利用该累积耗费阻力面构建CA模型的适宜性规则，对上述过程进行模拟，提高模型预测精度。

1.4 元胞自动机模型

CA模型是一种时间、空间、状态都离散的动力学模型，具有明显的时间与空间特征，CA模型的一般形式为[23]

Si,t+1=f(Si,t,SN,t)

(2)

式中Si,t+1——元胞i在t+1时刻的状态

Si,t——元胞i在t时刻的状态

SN,t——元胞i的邻域集合在t时刻的状态

f——转换规则

CA模型中，中心元胞在下一时刻的状态是其在当前时刻状态及其邻域集合状态的函数，准确定义该函数对于CA模型的模拟精度具有关键作用[24-27]。由于生态用地的演化过程是一种非线性的复杂动力学过程，当前元胞及其邻域状态影响中心元胞过程很难用简单的规则定义，因此本文采用人工神经网络模型提取CA模型转换规则，对土地利用模拟过程中的现有规则进行改进[28-29]，以提高CA模型的模拟精度。

本研究中，元胞形状为正方形的栅格像元，尺寸为30 m×30 m，采用3×3经典摩尔邻域定义邻域空间，模型初始栅格数据中的元胞状态定义为研究区景观格局类型，依次为耕地、林地、草地、建设用地、水体、未利用地，根据中心元胞及邻域元胞状态，利用ANN模型提取元胞转换规则，利用该规则计算元胞演化为生态用地的概率。

1.5 人工神经网络模型

ANN模型是20世纪80年代以来人工智能领域兴起的研究热点[30]。ANN模型在复杂的非线性系统的模拟方面具有明显优势，其自组织、自学习、联想以及记忆的优势能够有效简化CA模型，从原始训练数据中提取CA转换规则，避免主观因素影响，提高模拟精度[31]。

由于ANN模型在解决此类非线性问题中的明显优势，国内外学者很早就开始了相关研究，尝试利用ANN模型提取CA模型的转换规则[32]。此类研究往往把邻域规则和适宜性规则统一放入ANN模型进行规则提取，然而在模拟生态用地演化的CA模型中，相比适宜性规则，邻域规则部分更难定义也更难被人类理解，因此本文利用ANN模型提取CA模型邻域部分转换规则，适宜性规则采用MCR模型构建。

采用经典BP神经网络算法，网络结构共3层，分别是输入层、隐含层、输出层。模型的输入层共7个节点，隐含层共5个节点，采用tansing激励函数，输出层共1个节点，采用sigmoid激励函数，其输出为中心元胞转换为生态用地的概率，值域为0～1。输入变量及其取值范围如表1所示。

表1 ANN模型输入变量及取值范围Tab.1 Input variables and their ranges of ANN model

1.6 MCR-ANN-CA模型

以CA模型为基础框架，耦合MCR模型与ANN模型，构建MCR-ANN-CA模型用于包头市生态用地的演化模拟，技术路线如图1所示。模型模拟的基本步骤为：

图1 技术路线Fig.1 Technical roadmap

(1)利用MCR模型，综合考虑多种土地适宜性因子，构建研究区范围内各像元演化为生态用地的累积阻力面。

(2)利用历史数据的转换情况及邻域特征，构建训练数据，训练ANN模型，提取CA模型邻域规则。

(3)以当前土地利用为CA模型的输入数据，针对任一元胞，利用训练好的ANN模型，根据其邻域状态判断其演化为生态用地的概率P。

图2 阻力因子的阻力分布Fig.2 Spatial distributions of resistance values of resistance factors

(4)将P与累积阻力面中像元最大值相乘，其结果为S(取值范围为0到累积阻力面最大值)，若S大于当前元胞所对应的累积阻力，则该元胞演化为生态用地，否则，该元胞演化为非生态用地。

(5)将模拟结果转换为栅格图像并输出结果。

2 结果与分析

2.1 基于MCR模型的适宜性规则构建

对土地利用、NDVI、坡度、政府规划工业园区、人口密度、水体距离、高程等多种因子进行空间化与标准化处理，结果如图2所示。其中将NDVI、坡度、人口密度、高程因子数据原始值进行归一化处理作为其阻力。通过对不同用地类型阻力赋值，并对赋值结果归一化处理得到土地利用数据。政府规划工业园区阻力赋值为1，其他区域赋值为0，将数据栅格化后作为一个阻力因子。

根据包头市2016年土地利用数据，提取包头市生态用地范围，如图3a所示，由提取结果可知，包头市生态源地主要分布在北部草原带、大青山一带及黄河沿岸，其他区域也有零星生态用地分布，大致形成了三屏多点的格局。

利用GIS叠加分析，将图2中各因子进行空间叠加，得到包头市各个空间位置演化为生态用地的阻力面，结果如图3b所示。基于该阻力面，采用MCR模型，对源地进行累积耗费阻力计算，量化由近及远的空间范围内各像元演化为生态用地面临的阻力，得到如图3c所示的累积生态阻力面。由图3可知，包头市南北部生态阻力较低，中部累积生态阻力较高，尤其图3c中A区域及B区域，累积生态阻力达到包头市最高水平。

图3 包头市生态用地演化累积生态阻力面构建Fig.3 Construction of ecological land transition minimum accumulative resistance surface of Baotou City

2.2 基于ANN模型的CA邻域规则提取

基于研究区2006、2011年景观格局分布数据，通过Matlab随机抽取3 000个景观像元，统计2006年景观格局分布数据中，以3 000个像元为中心的3×3邻域内的各输入变量大小，形成ANN模型训练数据库。选其中70%进行模型训练，15%作为验证数据集，15%作为测试数据集。

采用量化共轭梯度法对ANN进行训练，模型经57次迭代后均方根误差达到最小，为0.12，如图4所示，图中圆圈处为最优拟合点，此时验证集的决定系数R2为0.90，模型拟合精度较好地满足了本研究要求。

图5 MCR-ANN-CA与CA-Markov模拟结果对比Fig.5 Comparison of simulation results of MCR-ANN-CA model and CA-Markov model

图4 ANN模型均方根误差Fig.4 Root mean square error of ANN model

2.3 MCR-ANN-CA生态用地模拟

将2011年景观格局数据输入训练后的ANN模型，判断各像元在邻域规则影响下是否会演变为生态用地，进而根据MCR模型生成的累积耗费阻力面，判断其最终演变方向,将模型模拟结果输出为栅格图像。包头市2016年实际生态用地8 409.32 km2，MCR-ANN-CA模型的模拟结果中包头市生态用地面积为8 670.01 km2，如图5a、5b所示。

2.4 模型模拟精度对比分析

为验证模型模拟精度，利用CA-Markov模型对包头市2016年生态用地进行模拟，结果如图5c所示。以2016年生态用地实际分布为参照，利用Idrisi Selva软件中的Cross Tab模块对2个模型的模拟结果进行分析，定量分析模型的模拟精度，结果如表2所示。表中KIA指数为卡帕一致性指数。

表2 模型模拟精度评价Tab.2 Accuracy assessment of simulation result

由表2可知，MCR-ANN-CA模型与CA-Markov模型模拟结果中，生态用地的面积均大于实际值，但总体上与2016年生态用地的实际分布保持了较高的一致性，其中MCR-ANN-CA模型的模拟精度略高于CA-Markov模型，二者的相对误差分别为3.10%与5.31%，KIA指数分别为0.89和0.87。

相比CA-Markov模型，MCR-ANN-CA模型通过ANN模型提取了元胞自动机邻域内的转换规则，同时利用MCR模型构建累积耗费阻力面，对不同空间位置元胞演化为生态用地的阻力进行了量化，因此模拟精度得到了进一步提高。

3 结论

(1)利用ANN模型提取了元胞自动机的邻域规则，同时利用MCR模型构建累积耗费阻力面，基于MCR-ANN-CA模型对包头市生态用地演化情况进行模拟，模拟精度较高。

(2)将MCR-ANN-CA模型与CA-Markov模型模拟结果进行对比，KIA指数分别为0.89和0.87，相对误差分别为3.10%和5.31%，MCR-ANN-CA模型对包头市生态用地的演化过程具有更高的模拟精度。