莫惠萍，黄宇斌，胡欣雨，林美玲，陆东芳

(福建农林大学 园林学院，福州 350002)

随着城市化进程的加快，人类对于自然资源的开发和利用不断加大，同时也带来了气候、环境、土地利用变化等各方面的问题.

景观格局变化是人为和自然因素共同作用下生态环境变化的综合反映[1]，它能够反映区域内的生态状况及空间变化[2].土地是人们生存和发展的基本场所[3]，它的变化一直受到人类活动的影响，但人类对于土地的不断开发和对土地覆被的改造都是人类发展过程中不可避免的.对土地利用景观格局变化与预测进行研究，能够帮助人们有效了解其成因与机制，对合理的土地利用计划和管理具有重大意义.

目前，国内外对土地利用景观格局的研究主要集中在时间和空间的动态格局演变、各种类型的数量转移、时间和空间的模拟和预测等方面.所运用到的模型主要有马尔可夫链(Markov)模型[4-5]、元胞自动机(CA)模型[6-8]、CA-Markov 模型[9-11]、Ann-CA-Markov 模型[12]、FLUS 模型[13-14]等.其中CA-Markov 模型主要用于模拟各种用地类型间的相互转化的趋势[15]，它很好地将CA模型和Markov 模型两者的优点有效结合，在模拟土地变化的时空格局和避免受到人为因素的干扰方面具有较好的表现.本研究所采用的MCE-CA-Markov 模型由多准则评价(MCE)方法与CA-Markov 模型结合，它量化了自然和社会经济等因素，补充了CA 模型的转换规则，使模型更加优化[16-17]，具有科学性和实用性.

珠三角地区是我国改革开放的前沿阵地.21世纪以来，随着工业化和城市化进程的大大加快，区域和生态系统土地利用结构频繁调整[18]，城镇建设与生态环境的协调成为珠江三角洲城市群需迫切解决的难题[19].佛山市地处珠江三角洲腹地，是我国工业化、城镇化快速发展的典型代表[20].改革开放以来，佛山市始终以工业产业为主导，导致区域内生态环境遭到破坏，严重制约了佛山市的可持续发展.如今随着广佛都市圈和粤港澳大湾区战略的大力推进，城市扩张加剧、人地矛盾剧烈等问题仍然存在，生态保护与经济社会发展之间的矛盾愈发突出，佛山市未来可持续发展面临严峻挑战.

因此，本文选取佛山市为研究区，以佛山市2005，2010 及2015 年的Landsat 遥感影像作为数据源，对佛山市2005—2015 年土地利用景观类型的分布格局和空间演变进行分析，并考虑自然和社会2 大影响因素，基于MCE-CA-Markov 模型对佛山市未来的土地利用景观格局进行模拟和预测，揭示研究区域的景观格局时空变化特征及规律；结合佛山市发展特点提出区域内可持续发展的对策，以期为同类型区域的景观格局演变和模拟预测提供参考，为研究区内未来土地利用计划和生态保护提供理论支持.

1 研究区概况

佛山市位于广东省的中南部(112°2 2′～113°23′E、22°28′～23°34′N)，东接广州，西靠中山，目前是粤港澳大湾区的极点城市以及广佛都市圈的重要组成部分；全市土地总面积3 797.72 km2，下辖禅城、南海、顺德、高明和三水5 区(见图1)；境内地势平缓，坡度多＜15°.改革开放以来，佛山市一直是广东省经济最发达的城市之一，其经济总量在中国大陆城市中居于前列.然而，其经济社会的快速发展也导致了区域内城市的快速蔓延、农业用地不断减少、生态环境恶化等问题的不断涌现.

图1 研究区概况

1.1 数据来源及预处理

本研究所用数据主要包括佛山市2005 和2010 年Landsat5 TM 及2015 年Landsat8 OLI 遥感影像数据、ASTER 30 m 数字高程数据和2019年全国路网数据.其中，Landsat 遥感数据主要来源于USGS 平台(https://earthexplorer.usgs.gov/)；ASTER 30 m 数字高程数据来自地理空间数据云平台(http://www.gscloud.cn/)；路网数据源自OpenStreetMap.遥感数据的处理主要通过ENVI完成，以2005，2010 和2015 年的遥感影像为数据源，分别对其进行几何校正、辐射定标、大气校正和监督分类，参照《土地利用现状分类标准》[21](GB/T21010—2017)和研究区域的实际情况，将各年份的景观类型划分为耕地、林地、草地、水体、建设用地和未利用地6 种用地类型，并通过精度验证Kappa 系数使其达到0.8 以上，最终获得2005，2010 和2015 年3 期景观格局现状数据，如图2 所示.

图2 佛山市2005，2010 和2015 年景观格局现状

1.2 研究方法

1.2.1 土地利用动态度

土地利用动态度可以表现研究区的土地利用动态变化情况.本文主要运用了单一土地利用动态度和综合土地利用动态度2 个指标来反映研究区土地利用面积变化的剧烈程度.其中，单一土地利用动态度一般用以表达某一地类在一定时间范围内的数量变化情况[22]；综合土地利用动态度则用来描述区域内总的土地利用变化速度[23].其具体公式如下：

其中，k为单一土地利用动态度；Ua和Ub分别为研究初期和末期某一地类的面积；T为研究初期与末期的间隔时长，单位一般为a；本研究k为某一地类的年变化率，年变化量大小与k成正比，k值越大，则表示该地类的年变化量越大；R表示综合土地利用动态度；Ui为第i类土地利用类型研究初期的面积.

1.2.2 土地利用转移矩阵

土地利用转移矩阵一般用来表示某时期中，期初与期末各土地利用类型之间的转换关系[16]，它能够通过定量描述的方式，更为直观地将各种土地利用类型的变化方向全面、具体地表现出来[24].其表达式如下：

其中，S为土地利用类型的面积；i为期初的土地利用类型；j为期末的土地利用类型；n为研究区域各土地利用类型.

1.2.3 景观格局指数

景观格局分析是进一步研究景观功能及其变化的基础[20]，而景观格局指数能够定量描述景观的差异及其变化.本研究通过对同类型研究中与本研究地块相似的景观格局指数的选取进行综合分析，并根据各指标的生态学意义分别从斑块类型水平和景观水平2 方面选取更为符合本研究区域的指标共11 个.其中，斑块类型水平方面主要选取了斑块类型面积、斑块数量、景观聚集度指数、最大斑块面积指数、景观破碎化指数、散布与并列指数；景观水平方面选择了斑块密度、面积加权平均形状指数、蔓延度指数、香农多样性指数和香农均匀度指数.各指标生态学意义见表1.

表1 景观格局特征指标及其生态学意义

1.2.4 CA-Markov 模型

CA-Markov 模型是CA 模型和Markov 模型的组合模型[25]，该模型可以定量并有效地模拟土地利用格局在时间上的时空动态变化[26].其模拟主要通过以下步骤实现.

1)设置转化规则.本研究区域面积较小，为了使模拟结果更符合实际，将佛山市景观类型栅格图像的像元大小设置为30 m×30 m，然后以2005 和2010 年分别作为初期和末期，通过IDRISI软件的Markov 模块，计算佛山市2005—2010 年的土地利用转移面积矩阵和转移概率矩阵.

2)利用多标准评价模块(MCE)制作适宜性图集.首先，根据研究区已出台政策中对生态的保护以及未来各用地类型的管控方向，设置空间分布概率图中的限制因素和影响因素；其次，通过层次分析法(AHP)确定各因子的权重，结合限制因素和影响因素对景观类型转化的影响，根据MCE 和COLLETION EDIT 模块将各景观类型的空间分布概率图集合，生成适宜性图集.

3)确定滤波器、起始年份和CA 循环次数.使用5×5 滤波器，以佛山市2010 年为起始年，将CA 循环次数设为5，结合2005—2010 年的土地利用转移概率矩阵和适宜性图集，对研究区2015 年的景观类型进行模拟.

4)检验模拟精度.第3)步的模拟完成后，将2015 年的模拟结果与2015 年景观格局现状进行对比，检验该模型的模拟精度，精度≥85%则表示模型的模拟精度较高.

5)预测研究区未来的景观类型.以2015 年为预测起始年，CA 循环次数设定为15，结合2010—2015 年的土地利用转移概率矩阵和适宜性图集，预测佛山市2030 年的景观类型分布.

2 结果与分析

2.1 土地利用结构变化及面积流转

佛山市2005—2015 年土地利用结构变化及动态度变化情况分别如表2、表3 和表4 所示.从土地利用结构变化进行分析(见表2)：2005—2015 年，耕地、林地和建设用地始终处于主导地位，其中林地面积所占比重逐渐下降，耕地和建设用地的面积则呈上升趋势，若进一步分析2015年土地利用情况，发现该态势依旧持续；水域的面积则呈逐年减少的趋势，其中2005—2015 年其占比由25.96%持续减少到了13.26%.

表2 佛山市2005，2010 和2015 年土地利用结构变化

从土地利用动态度分析(见表3 和表4)可知，2005—2015 年间各土地利用类型的动态度变化大小依次为：水域＞耕地＞建设用地＞未利用地＞草地＞林地.其中，在2005—2010 年间，佛山市土地利用各类型中，水域的变化最明显，其面积减少了470.542 5 km2，动态度为-7.95%；然后依次是未利用地、耕地和建设用地；草地、林地的动态度较小，基本处于稳定状态.2010—2015 年间，耕地较2005—2010 年动态度明显减小，这说明耕地流失现象依然严峻；水域在该研究时段由-7.95%上升到了-0.39%，这说明水域面积大量减少的情况得到了有效控制；建设用地在该研究时段内较2005—2010 年动态度也从3.35%下降到了1.03%，这说明建设用地的扩展速度在一定程度上得到有效控制；林地由于基数大，动态度变化虽然减弱，面积却也减少了16.403 4 km2，这说明林地的保护措施尚待加强；草地和未利用地的动态度虽然较高，但由于本身初期面积较小，数量上仍然没有大的变化，这说明草地和未利用地都相对稳定.

表3 佛山市2005—2010 年土地利用动态变化情况

表4 佛山市2010—2015 年土地利用动态变化情况

佛山市2005—2015 年的土地利用流转趋势如表5 所示.不考虑自身转化即未变化的部分，可以发现，在该时段内，耕地转化为建设用地量较大，转移了159.719 4 km2；林地同样是转化为建设用地的概率较大，达65.9%；草地较为稳定，转出量较少，但也有51.9%转为了建设用地；就水域而言，大多数转化为了耕地，占转出量的72.8%；未利用地大部分转成了建设用地；建设用地中有69.0%转为了耕地，但建设用地总体是在不断扩张的，并且占用了大量的耕地以及部分林地和草地，这也是城镇化的必然趋势.

表5 佛山市2005—2015 年土地利用类型转移矩阵 km2

2.2 景观格局变化分析

2.2.1 基于类型水平的景观格局指数分析

佛山市2005—2015 年基于类型水平的景观格局指数如表6 所示.从CA、NP和AI这3 个指数分析，研究区2005—2015 年间，耕地的斑块面积最大，且呈现持续变大的趋势，斑块数量则呈先增后减的变化规律，AI指数变化规律相同，以上变化说明该景观类型受到人类活动的干扰，开始朝破碎化方向发展，在后期趋于稳定；林地的CA、NP指数逐渐减少，AI指数变化幅度较小，说明该景观类型破碎化程度较低，但各斑块间的连通性不强；水域在该时段内的变化也十分明显，CA、NP和AI这3 个指数均减少，表明水域的破碎化程度较高；建设用地则表现为斑块数量减少，斑块面积和聚集度指数逐渐增加的变化规律，说明该时段内建设用地逐渐增多且连片发展；草地的NP、CA指数在研究时段内呈现前5 年不断减少，后5 年开始增加的变化趋势；未利用地整体变化幅度不大.

表6 佛山市2005—2015 年类型水平的景观格局指数

从最大斑块面积指数(LPI)分析，2005—2015年间，耕地的LPI指数呈现先减少后增加的变化规律，该景观类型的优势度先降低后凸显；林地和水域的LPI指数逐年降低，说明林地和水域的优势度在逐渐减弱；而建设用地在初始阶段的LPI值虽然偏小，但在2005—2015 年间呈现先减少后增加的态势，优势度已经逐渐高于林地和水域，说明建设用地的优势度也愈来愈明显；草地和未利用地的LPI指数较低且变化幅度较小，这主要是由于二者本身的面积基数小.

从散布与并列指数(IJJ)方面分析可知，在研究时段内，耕地的IJJ指数逐年减少，说明该景观类型聚集程度降低，分离程度增加；林地和草地的IJJ指数均逐年增加，说明2 类景观类型的斑块聚集程度增加；水域的IJJ指数则呈现先减少后增加的趋势，说明该景观类型斑块间近年来受人类活动干扰强度降低，聚集度增加；建设用地的散布与并列指数最小，说明其聚集程度最低，容易受人类活动的影响，且与之相邻的斑块，大多也容易受人类活动的影响，排列比较复杂.

2.2.2 基于景观水平的景观格局指数分析

佛山市2005—2015 年基于景观水平的景观格局指数如表7 所示.由表7 可知，佛山市2005—2015 年整体的景观格局变化相对稳定.其中，研究区的斑块密度出现先减少后增加的变化，说明该时段内景观破碎化程度逐渐增加.SHAPE_AM指数具有相同特征，说明区域内景观的空间复杂程度增加，景观斑块的形状逐渐变得不规则.而SHDI和SHEI指数在2005—2015 年间一直处于减少的状态，说明研究区由于外界的干涉，导致了自然生物多样性减少，景观分布不均匀.CONTAG指数则呈持续上升趋势，说明研究区的各斑块类型的连接性较好.

表7 佛山市2005—2015 年景观水平的景观格局指数

2.3 土地利用景观格局变化模拟及预测

2.3.1 土地利用转移概率矩阵

佛山市2005—2010 年、2010—2015 年土地利用类型单位转移概率矩阵由IDRISI 软件中的Markov 模块计算得到，结果如表8 和表9 所示.

表8 佛山市2005—2010 年土地利用变化转移概率矩阵

表9 佛山市2010—2015 年土地利用变化转移概率矩阵

2.3.2 适宜性图集制作

1)限制因素.限制因素主要选取了水域和建设用地.通过IDRISI 中的reclass 工具对提取的水域和建设用地进行重分类，将水域和建设用地的部分赋值为0，其他背景值赋值为1，生成二值化图像.

2)影响因素.影响因素主要选择了高程、坡度、距道路距离、距城镇用地距离和距水域距离5 个指标.首先，利用IDRISI 中的FUZZY 模块中的“S”型和“J”型函数对各因子进行标准化处理，得到各地类的标准化图层；其次，按照各自的转换规则，分别将制定的限制条件和标准化后的影响因子图层导入MCE 模块，选择加权线性合并法，并赋予各因子权重值，得到各景观类型的空间分布概率图；最后，根据 MCE 和COLLETION EDIT 模块将各景观类型的空间分布概率图集合，生成所有景观类型的适宜性图集.

2.3.3 佛山市2015 年土地利用景观格局演变模拟

基于以上数据的准备，以佛山市2005 年为初始年，利用IDRISI 中的CA-Markov 模块，输入土地利用变化转移概率矩阵和各景观类型的适宜性图集，模拟循环次数设置为10，实现对佛山市2015 年土地利用景观格局的模拟，结果如图3(b)所示.为了验证该模型的模拟精度，利用IDRISI中的CROSSTAB 模块对2015 年模拟预测结果与2015 年现状进行了对比(见图3)，得到的Kappa系数为0.935 2，说明该模型模拟效果较好，与实际较符合，可行性较高.

图3 2015 年佛山市景观类型现状与模拟结果对比

2.3.4 佛山市2030 年土地利用景观格局预测

以佛山市2015 年的土地利用为基期，利用Markov 模型得到佛山市2010—2015 年的土地利用转移概率矩阵，以及基于MCE 模块构建的适宜性图集，将循环次数设置为 15，重复CA-Markov 操作，实现了对佛山市2030 年景观土地利用景观格局的模拟，结果如图4 所示.

图4 2015 年景观类型现状与2030 年预测结果对比

为了直观分析佛山市2015—2030 年各景观类型的数量变化，本文对佛山市2030 年预测结果与2015 年现状中各景观类型的面积进行了统计对比，结果如表10 所示.由表10 可知，2030 年，佛山市各景观类型变化继续呈现出耕地与林地的适当减少，建设用地的迅速增加，以及水域、草地与未利用地不明显变化等特点，这种变化趋势是经济发展较好的区域普遍存在的.其中，耕地面积较2015 年减少了171.31 km2，占比减少了4.51 %.为防止未来耕地持续减少，建议政府应加大政策管控，严格控制耕地保护红线，限制盲目开垦.林地面积减少了102.45 km2，占比减少了2.70%，其变化强度要大于2005—2015 年该类型的变化，且林地作为研究区域内面积大小仅次于耕地和建设用地的景观类型，其生态地位十分重要，为了实现区域内生态环境和社会经济的可持续发展，政府应采取措施严格控制林地的减少量.草地面积所占百分比由 0.27%增长到了0.35%，面积增加了3.177 km2，占比增加了0.08%，变化幅度较小，但草地具有涵养水源、保持水土等生态功能，其地位也不能忽视.水域面积所占百分比由13.26%上升到了13.59%，占比增加了0.33 %，说明佛山市2015—2030 年间水环境得到改善，这与佛山市人民政府2015 年提出的恢复河流水生态系统等举措有关.建设用地面积由2015 年的31.20%上升到了2030 年的38.00%，面积增加了258.15 km2，虽增加强度较2005—2015年有所减弱，但建设用地的转换具有不可逆性，在规划时，应合理控制建设用地的增加量.未利用地面积整体变化量较小，尚有1.395 9 km2处于闲置，可见，未来佛山市的城市建设应将未利用地充分利用.

表10 佛山市2015 与2030 年土地利用景观格局面积变化对比

3 结论与讨论

3.1 结论

1)从土地利用动态变化来看，佛山市2005—2015 年间，各用地类型的动态度变化依次为：水域＞耕地＞建设用地＞未利用地＞草地＞林地；土地利用流转趋势上，耕地、林地、草地、未利用地转化为建设用地的概率较大，水域大部分转化为耕地，小部分也转为建设用地，建设用地则主要转为耕地.

2)从景观格局变化来看，2005—2015 年佛山市整体景观格局变化相对稳定，景观破碎化程度逐渐升高，空间复杂程度增加，斑块形状趋于不规则，景观分布越发不均匀.耕地斑块由于受到人类活动的干扰，存在一定破碎化的情况但总体趋于稳定；林地和水域的优势度逐渐减弱，破碎化程度较高；建设用地逐渐增多且连片发展，破碎化程度较低，但聚集程度最低.

3)通过MCE-CA-Markov 模型对佛山市2030年景观格局的模拟预测结果看，佛山市林地和耕地的面积持续减少，水域和建设用地面积出现增加，草地和未利用地基数小，总体变化不明显.

3.2 讨论

在利用MCE-CA-Markov 模型模拟佛山市未来的景观格局时，还存在一定的不足.如在选择限制因素制作模拟所需的适宜性图集时考虑的因素不够全面，地理、社会、经济、资源环境等客观因素都应考虑进去，这样模拟结果才会更精准；其次是对尺度的考虑不够，本文主要以30 m×30 m 的栅格大小进行模拟，但栅格大小对于模拟精度会有一定的影响，所以导致模拟结果与实际情况产生一定的差距；另外，影响模拟精度的还有时间间隔长短以及滤波器大小的设置.基于此，今后的研究应进一步探究不同栅格大小、时间间隔长短和滤波器大小对模拟精度的影响，以选择最适合的尺度进行模拟，进而提高模拟精度.