基于CA－Markov模型的海湾土地利用动态预测

吴季秋，俞花美，肖明，葛成军，唐文浩*

1.热带作物种质资源保护与开发利用教育部重点实验室，海南大学环境与植物保护学院，海南 海口 570228

2.中国科学院生态系统网络观测与模拟重点实验室，北京 100101

3.海南大学经济与管理学院，海南 海口 570228

土地作为人类生存和发展的基本物质条件，其资源配置的合理性已成为地区经济的主要制约因素，关乎着土地资源的可持续发展与国家和民族的命运前途。2009年全国土地利用变更调查结果显示，我国人均耕地面积只有1.38亩，还不到世界人均水平的1/3，其中九个省区人均耕地面积低于1亩，三个省区人均耕地面积低于0.5亩［1］。我国的土地资源本身内部禀赋不足，而土地退化、耕地锐减、生态环境质量恶化、用地结构不够合理等现象进一步加剧了经济发展诉求与土地资源稀缺间的矛盾。面临严峻的现状，开展土地利用趋势变化分析，探求与生态环境协调发展相适宜的土地利用模式，一方面可为区域土地利用的可持续发展提供理论依据，另一方面为土地管理和资源优化配置提供支持和借鉴。运用GIS等先进手段进行分析，既是GIS与土地管理决策结合的理论尝试，又是落实《国土资源部关于加快推进土地利用规划数据库建设的通知》(国土资发〔2011〕3号)的指示精神的具体体现。

元胞自动机-马尔科夫(CA-Markov)模型综合了元胞自动机(cellular automata，CA)和马尔科夫模型各自的优点，即元胞自动机模型强大的模拟复杂系统空间演变功能和马尔科夫模型长期预测的优势，将二者有机结合在一起，有利于精确模拟土地利用变化的时空演变过程，是目前区域土地利用空间格局动态模拟预测研究中应用最多且较为有效的研究方法［2-3］。笔者以海南省文昌市八门湾地区为研究对象，采用马尔科夫链分析研究区域1998年、2009年土地利用面积转移矩阵及土地利用转移概率矩阵，并根据马尔科夫的结论，依据土地适宜性原则和环境友好型土地利用模式［4］，运用CA-Markov模型动态模拟预测该区域2020年的土地利用变化趋势，进一步揭示不同的土地利用约束条件对于土地利用/土地覆盖变化的影响，由此分析和探讨土地利用与生态环境协调发展的理论途径。

1 研究区域概况

研究区域八门湾位于海南省文昌市东南部(19°28'15″N ～ 19°43'44″N，110°42'35″E ～110°57'14″E)。八门湾因东边的文教河与西边文昌河在海湾出海口交汇，形成“八”字而得名，所辖湾内海域面积约680多hm2。流经研究区域的河流众多，主要以文教河、文昌河及其各支流所形成的水系为主。研究范围包括清澜镇、文城镇、头苑镇、东阁镇、文教镇和东郊镇，总面积约35377.73 hm2。

2 数据及预处理

2.1 数据来源

选取1988年、1998年、2009年三期TM影像，目视解译标志采用海南省生态监测2009年、2010年两期全国生态环境监测野外定点数据(海南站)建立目视解译数据库。数据处理按照《海洋生态调查指南》、《土地利用现状调查技术规范》要求制定。相关参考资料有《文昌市土地利用总体规划(2006—2020年)》、《文昌市产业布局图》、《海南省环境公报(1998—2009)》、《海南统计年鉴》。

2.2 数据预处理

采用ERDAS IMAGINE 2011的raster模块，在其中完成影像的波段组合、影像镶嵌、大气校正、去云等操作。利用数字化海南省1∶10万地形图对影像进行几何校正，校正后的像元误差控制在0.5个像元以内进行重采样。土地覆盖类型分类参照《土地利用现状调查技术规范》和全国生态环境监测评价体系要求分为19个类型(表1)。

表1 研究区域土地覆盖分类体系Table 1 Land covers classification system in study region

3 分析方法

CA-Markov的模拟过程大致分为两个阶段:先用Markov模型对于土地利用的时空推移概率进行计算，然后建立土地利用适宜性规则集，综合前一阶段的概率集在CA-Markov模型中推导出最终结果(图1)。

图1 研究路线Fig.1 Research route

3.1 Markov模型描述

目前，在众多的土地利用分析模型、土地利用/土地覆盖变化模拟和城市景观格局变化预测研究中，Markov模型因其长期预测的明显优势而成为应用较为广泛的一种土地利用变化建模方法［5-6］。Markov过程是基于马尔科夫链的具有无后效性特征的一种状态和时间都离散的随机运动过程，即一个运动系统在未来某个t+1时刻的状态仅与当前已知的t时刻所处的状态有关，而与t时刻之前的状态无关。马尔科夫链分析的原理是通过一个运动系统的不同初始状态及各种状态的转移概率来研究各种状态之间的转换情况，进而实现对未来状态变化趋势的预测。在马尔科夫链中，系统状态的转移可用概率矩阵来描述，其数学表达式为［7］:

式中，P为某一土地利用类型从初始状态到下一状态的转移概率矩阵;m和n为土地利用类型个数;Pij为土地利用类型i转变为土地类型j的转移概率。

自任何一个状态开始，经过一次转移，即会出现该系统所有状态中的一种，该过程可表示为:

式中，Rt为初始状态;Rt+1为下一状态。因此，一次马尔科夫模型可由已知t时刻的状态推知t+1时刻的状态。

在Idrisi的 Modeling菜单下点击 Ecological environmental/Simulation Models，选择 Markov模块，选择1998年、2009年八门湾区域的土地利用/土地覆盖图像，将这两期图像的间隔时间和预测时间周期均设置为11年，比例误差设置为0.15，运算得出该区域1998年、2009年的土地利用面积转移矩阵及土地利用转移概率矩阵。

3.2 CA-Markov模型原理

元胞自动机(CA)，是时间、空间、状态都离散，每个变量都只有有限多个状态，且状态改变的规则在空间的相互作用及时间上的关系皆表现为局部特征的网格动力学模型［8-12］。元胞自动机实际上是一个包含元胞(cell)及状态信息(state)、元胞空间(lattice)、邻居(neighbor)、规则(rule)及与时间(time)相关的转换机制的函数［13］。

CA 模型可表示为［14］:

式中，S为元胞有限、离散的状态集合;N为元胞的邻域;f为局部空间的元胞状态的转化规则。

在Idrisi的Modeling菜单中调用CA-Markov模型。以2009年土地利用类型图为基础，预测时间间隔为11年，模拟2020年八门湾的土地利用类型转化状况。上述元胞自动机构成的六部分参数设置如下。

元胞:是元胞自动机的最小运算单元，按照一定的空间大小划分元胞单位。文中元胞自动机的元胞即为图像栅格，元胞大小设置为30 m×30 m。

状态:每个元胞的状态可以定义为多个变量。在空间分析中，其状态是指元胞所具有的空间位置、大小、土地利用类型、周长等空间属性。文中的元胞状态即为其土地利用类型数字化属性。

元胞空间:是指元胞在空间位置分布上网点的集合。在地理空间模拟的应用实践中，元胞空间是具有空间位置属性的概念，与空间栅格数紧密相关。

邻居:元胞下一时刻的状态是由元胞及其邻居这一时刻的状态所决定的。滤波器的工作原理是依据邻居距离元胞的远近位置关系来创建具有明显空间意义的权重因子，并作用于元胞，进而促使元胞状态的转变，文中选用的是5×5滤波器。

规则:元胞自动机的函数运算准则是依据其设定的规则，将形成的约束函数用以指导元胞转变的方向。在地理空间分析中，规则可由一系列的土地利用控制条件构成，按照经济适用性、距离、高程和环境保护等适宜性实际需要，制定运行规则。Idrisi中依据多目标决策模块(multi-criteria evaluation，MCE)确定的适宜性图像集来确定［15］。以距离、坡度、降雨量为约束条件，将道路、水库、居民点、坡度、降雨七个因子二值标准化为布尔图像，即将适宜开发的区域赋值为1，反之赋值为0。

距离约束条件，从降低开发费用的角度出发，设置距离道路小于500 m区域适宜开发，赋值为1;大于或等于500 m的区域不适宜开发，赋值为0［16］。考虑到生产、生活产生的垃圾会导致地下水污染和周边居民点的非点源污染，将距离水库100 m缓冲区范围设为非适宜开发，其余区域设置为适宜开发。

坡度约束条件，我国《水土保持工作条例》规定，25°以上的陡坡地，禁止开荒种植农作物。据此，采用数字化高程DEM数据，利用ArcGIS中的Slope工具得到坡度图，再重分类为25°以上和以下两个区域划分禁耕区和适宜种植区。

降雨约束条件，采集了分布于区域内文昌河、文教河、清澜红树林保护区的五个气象监测站的降雨数据，统计了2003—2009年的降雨记录，结果显示，年均降雨量为913～1300 mm的西部农业生产总值明显低于年均降雨量为1220～1935 mm的东部。对五个监测站点的降雨量进行克里格插值，并重新分类后得到降雨分布图。年均降雨量在1200 mm以下的划定为不适宜农业种植区域，年均降雨量大于1200 mm的划定为农业生产种植适宜区。

时间:元胞自动机是一个在时间变化上离散的动力系统，即时间t是具有连续等间距特性的整数值。元胞在t+1时刻的状态只取决于t时刻该元胞及其邻近元胞的状态。设置模拟预测图像的时间间隔为11年。

4 结果与分析

4.1 Markov模型模拟结果分析

分析表2、表3可知，在1998年、2009年的土地利用变化中，养殖水面面积增加较大，其主要是由水田、果园和旱作耕地转化而来，其中水田的转移面积最大，占696.60 hm2，转移概率为0.100，果园和旱作耕地转移面积次之，分别为332.91和298.62 hm2，其转移概率分别为0.099和0.093;水田面积的增加主要是由果园和养殖水面改变用途转变而来，其转移概率分别为0.340和0.264，总面积为1721.07 hm2，占该土地利用类型总转化面积的35.54%。分析2009年土地利用类型面积减少的状况，红树林减少的面积主要变成了水田和养殖水面，其转移概率分别为0.262和0.269，共计267.66 hm2，占该土地类型转化总面积的57.88%;减少的旱作耕地主要转变为果园，沿海附近的耕地转变用途变成了养殖水塘，转移概率分别为0.156和0.093，共计799.47 hm2，占该土地利用类型转移面积的28.98%。

通过对 Markov模型模拟结果的分析，对比1998年和2009年土地利用转移矩阵，探讨引起八门湾区域土地利用转移特征的主要原因:1)八门湾海湾地区的养殖水面迅速增长，表明该区域沿海附近为追求经济效益，大力发展高位池塘养殖，从2009年的遥感卫星影像上可以清晰看到沿文教河入海口至八门湾出海口沿岸一线布设建造了大量的围栏养殖塘。这种局面虽然短期内带来了一定的经济回报，但却带来了一系列较为严重的生态问题，未经处理或处理不充分的养殖废水的成倍增加给海洋生态环境造成了较严重的威胁和破坏，危及该地区的生态安全。2)八门湾的红树林地带是海南岛著名的红树林景观之一，但是，近年来由于高位池养殖的方兴未艾，导致大片的红树林生态系统转变为海产品养殖池塘，同时受到经济利益的驱使，不断围海造田，把沿岸红树林滩涂湿地改造成为农田，进而导致红树林面积的减少。

4.2 CA-Markov模型结果分析

通过对于马尔科夫土地利用转移矩阵的特征分析发现，当前的土地利用方式仍然存在许多不合理之处，1998—2009年的11年间，养殖水面和水田面积呈稳定增长的态势，而红树林和旱作耕地面积却明显减少。这对当前八门湾生态环境的保护以及社会经济环境的可持续发展构成了一定的威胁。基于此，在运用CA-Markov模型模拟预测该地区2020年的土地利用状况时，相应地设置了一系列土地开发利用适宜性规则，以距离、坡度和降雨三方面的约束条件为出发点来制定元胞运行规则，其计算结果如图2和表4所示。

图2 CA-Markov模型预测2020年八门湾土地利用类型Fig.2 Land use type in 2020 in Bamen bay

预测结果表明:在土地利用适宜性规则形成的函数约束条件下，2020年防风林的面积增长迅速，较2009年增加了307.91 hm2，土地利用类型变化单一动态度指数达0.1585，天然林和红树林的面积也有所增加，动态度指数分别为0.0243和0.0184;伴随防风林、天然林和红树林面积的增长，养殖水面、果园和水田的面积呈减少趋势，养殖水面减少了1102.75 hm2，动态度指数为0.0322，果园和水田的面积分别减少了1088.18和2057.43 hm2，动态度指数分别达到0.0174和0.0172。由图2可知，沿文昌河、文教河至八门湾海湾沿岸一线，红树林和防风林的覆盖面积有所恢复和增长，相应的养殖水面有所缩减，陆域部分果园周边的防风林和红树林植被覆盖也有所增加。这一预测结果说明人为制定的一系列适宜土地开发利用和保护的经济实用性及生态环境保护规则，不仅有利于保护该区域生态系统的平衡，维护生物多样性及生态安全，而且可以保障该区域经济社会与生态环境的协调可持续发展。

?

?

表4 2009年和2020年八门湾土地利用类型空间变化Table 4 Spatial distribution of land use type change between 2009 and 2020 in Bamen bay hm2

4.3 预测精度检验

CA-Markov模拟结果的精度验证在目前还未形成较为统一的方法，通常使用的主要有整体逐个检测与随机抽样检测两种。随机抽样检测适用于验证对象的分布较为均匀，代表性较强的情况。虽然该方法因为对象数据量少而处理快速，但其主观随意性大，容易因为样本选取的典型性不够而产生偏差。整体逐个检测可以检验到参与模拟的每一个单元，可信度较高，但是计算冗繁，需要较长的处理时间。考虑文中研究区域数据量较小，故采用整体逐个检测法。将模拟出的2009年预测结果〔图3(b)〕与实际解译的2009年土地利用类型结果〔图3(a)〕进行求差栅格运算，前后一致的区域求差后土地利用类型属性编码值为0〔图3(c)〕。最后以提取出土地利用类型编码值为0的区域总栅格数，除以研究区域总栅格数，即得模拟精度。由表5可见，2009年预测的模拟精度为77.30%，反映出该模拟结果可信度较好，CA-Markov预测的土地利用类型具有较好的可靠性和适用性。

图3 八门湾土地利用类型预测精度验证Fig.3 Verification of forecast accuracy for land cover types in Bamen bay

表5 2009年预测精度对比Table 5 Comparison of forecast simulation accuracy in 2009

5 结论

(1)马尔科夫链分析得出，研究区域养殖水面和水田的面积呈增长趋势，相应地使得红树林和旱作耕地等土地利用类型的面积减少。

(2)基于以上的发展变化趋势，从可持续发展的角度出发，根据土地利用适宜性原则，调整驱动因子空间分布规则，由此得到2020年的预测结果，防风林、天然林和红树林的面积将有所增加，动态度指数分别达0.1585、0.0243和0.0184;而养殖水面、果园和水田面积将有所减少，动态度指数分别为0.0322、0.0174 和0.0172。

(3)按照有利于生态环境改善，不影响当前经济发展，降低经济开发成本的原则，将距离、坡度和降雨三个因素对土地利用适宜性约束条件进行了调整，模拟结果显示，研究区域的生态环境在该土地适宜性原则的约束下，天然植被有所恢复，人为的生产性土地扩张得到相应的控制。

6 不足与展望

对于CA-Markov的适宜性原则主要来自大量文献的比较研究和土地利用政策综合得出，由于缺少八门湾其他生态影响因子的实际研究成果和依据，因此只从距离、坡度和降雨三个方面的自然条件，同时仅结合政策或经济约束条件进行模拟，未考虑更多其他可能的潜在影响因素，如人口数量的快速增长与产业布局的不断调整都会使得区域的土地利用发生变化，从而对整个区域的生态格局产生各种影响。与此同时，随着八门湾地区部分吸引旅游开发投资的产业政策的出台，也会产生与之对应的土地利用空间格局的变化。因此，在今后的研究中，不仅要丰富其他潜在影响因子，如放牧、捕鱼等人为影响因子的影响半径测定，而且需要进一步探讨按照不同的自然资源条件和经济发展水平、政策导向影响下的土地利用空间变化及其与土地承载力和适宜性的评价机制，以完善土地利用变化趋势研究的理论和手段。

［1］颜玉华:对我国18亿亩耕地“红线”的忧思［EB/OL］.中国改革论坛.(2012-01-21)［2012-03-01］.http://www.chinareform.org.cn/Economy/Agriculture/Forward/201201/t20120122_132973.htm.

［2］赵建军，张洪岩，乔志和.基于CA-Markov模型的向海湿地土地覆被变化动态模拟研究［J］.自然资源学报，2009，24(12):2178-2184.

［3］赵冠伟，陈颖彪，陈健飞，等.CA-Markov模型的空间尺度敏感性研究［J］.地理科学，2011，31(8):897.

［4］张秋玲，李东敏，邵亚杰.基于环境友好型模式的土地利用与生态环境协调研究:以河南省信阳市为例［J］.安徽农业科学，2009，37(23):11130-11132.

［5］ALEJANDRO F S，MIGUEL M R，OMAR R M.Assessing implications of land-use and land-cover change dynamics for conservation of a highly diverse tropicalrain forest［J］.Biological Conservation，2007，138(1/2):131-145.

［6］杨国清，刘耀林，吴志峰.基于CA-Markov模型的土地利用格局变化研究［J］.武汉大学学报:信息科学版，2007，32(5):414-418.

［7］唐启义，冯明光.实用统计分析及其DPS数据处理系统［M］.北京:科学出版社，2000:95-105.

［8］朱海涌，李新琪，仲嘉亮.基于CA-Markov模型的艾比湖流域平原区景观格局动态模拟预测［J］.干旱环境监测，2008，22(3):134-139.

［9］李书娟，曾辉，夏洁，等.景观空间动态模型研究现状和应重点解决的问题［J］.应用生态学报，2004，15(4):701-706.

［10］LI X.Neural-network-based cellular automata for simulating multiple land use changes using GIS［J］.International Journal of Geographical Information Science，2002，16(4):323-343.

［11］LOPEZ E，BOCCO G，MENDOZA M，et al.Predicting land-cover and land-use change in the urban fringe:a case in Morelia City，Mexico［J］.Landscape and Urban Planning，2001，55(4):271-285.

［12］侯西勇，常斌，于信芳.基于CA-Markov的河西走廊土地利用变化研究［J］.农业工程学报，2004，20(5):286-291.

［13］HAN J，HAYASHIA Y，CAO X，etal.Application of an integrated system dynamics and cellular automata model for urban growth assessment:a case study of Shanghai，China［J］.Landscape and Urban Planning，2009，91(3):133-141.

［14］孙贤斌，刘红玉，李玉凤，等.基于CA-Markov模型土地利用对景观格局影响辨识［J］.生态与农村环境学报，2009，25(1):1-7.

［15］何丹，金凤君，周璟.基于Logistic-CA-Markov的土地利用景观格局变化:以京津冀都市圈为例［J］.地理科学，2011，31(8):903-909.

［16］毛蒋兴，闫小培.城市交通干道对土地利用的廊道效应研究:以广州大道为例［J］.地理与地理信息科学，2004，20(5):58-61.▷