基于Logistic-CA-Markov耦合模型的彬州市LUCC多情景模拟

李世锋, 洪增林,, 薛旭平, 张锋军, 石 卫

(1.长安大学 土地工程学院, 西安710054; 2.陕西省地质调查院, 西安710054; 3.陕西省水工环地质调查中心, 西安710068)

土地利用/覆盖变化(land use/cover change，LUCC)影响区域社会经济、生态环境，进而引发区域生态系统、景观格局及功能的演变[1-3]。由于人类对土地利用的方式具有能动性和调控力，使得土地利用/覆盖变化导致许多自然现象和生态过程的变化，因此对土地利用变化的驱动力分析及模拟预测研究一直是LUCC研究热点领域之一。基于对区域土地利用过程变化及其驱动因子分析，并且进行多情景动态模拟，预测未来土地利用空间分布情况，可为区域土地资源合理开发利用、规划以及生态保护环境提供支持和借鉴。

目前，国内外对于LUCC的分析和模拟有很多种模型，有数量预测模型，如马尔科夫模型(Markov)[4]、Logistic回归模型[5]、灰色预测模型(Gary Forecast Model)等[6]；有空间预测模型，如元胞自动机模型(CA)[7]、CLUE-S模型[8-10]、FLUS模型等[11]；有耦合模型，如CA-Markov模型[12-14]、CA-MAS模型[15]、Logistic-CA-Markov模型等[16]。其中，Logistic回归模型常用于定量分析驱动因子与LUCC之间的相关性，并且能预测在当前的情况下出现某地类的概率，但模拟结果只能是数量上的变化，在空间分布上体现不出来[17]。马尔科夫(Markov)模型预测未来较长时间内土地数量有绝对的优势，但不能预测土地利用空间分布格局的变化[18]。元胞自动机(CA)模型可以有效模拟元胞之间的相互作用，在空间分析和模拟运算有独特且强大的能力，但存在局限性[19]。将Logistic回归模型、CA模型和Markov模型进行耦合，形成Logistic-CA-Markov模型，既能识别预判驱动因子与LUCC的显著关系并预测各地类分布的概率，为后续制作土地转变适宜性图集，校正元胞自动机的转换规则提供基础支撑，又能在数量和空间模拟土地利用变化，这样可以大大提高预测精度。

本文以陕西省彬州市为例，基于2009年和2019年土地利用数据、DEM数据、坡度数据、交通道路数据、规划数据、经济社会数据、水域数据等，通过GIS软件空间分析功能对土地利用结构和利用动态度进行时空演变分析，对影响彬州市土地利用变化的驱动因子进行Logistic回归分析，运用Logistic-CA-Markov模型，以研究区2009年为基期，模拟预测2019年土地利用变化，将预测的结果与2019年实际土地利用数据进行对比，计算其Kappa系数验证模拟精度，最后分别在自然发展、生态保护和限制城市过度开发3种情景下对研究区2029年土地利用变化进行模拟预测，分析得出结果。

1 研究区概况及数据来源

1.1 研究区概况

彬州市位于陕西省渭北高原西部、咸阳市西北部，陕甘边界地区，属陇东黄土高原塬梁丘陵沟壑区，地势西南高东北低，介于东经107°49′—108°22′，北纬34°51′—35°17′，全市土地总面积1 183.59 km2，约占陕西省总面积的6%，区域内海拔为715～1 500 m。泾河自西北向东斜贯中部，将彬州市分割成东北、西南两塬夹川道的地貌格局。研究区属于暖温带半干旱大陆性季风气候，具有雨热同季，寒暑极端，四季分明，光能资源丰富，降水时空分布不均的特点，年平均日照数2 210.8 h，年日照分率为52%，年均气温9.7℃，极端最高温40℃，极端最低温-22.5℃，全年无霜期172～177 d，年均降水量561.4 mm。彬州市共有5个土类，分别为黑垆土、黄绵土、红土、淤土和潮土，其中黑垆土和黄绵土分布最广，也是主要的农业土壤。

1.2 数据来源

本研究所需的土地利用数据来源于土地调查数据，依托项目为《彬州市国土空间规划》，这与常规的从Landsat-TM遥感数据中识别、解译和提取土地利用/覆被信息相比，避免了解译正确率这一误差，精度更高。研究区高程(DEM)数据来源于中国科学院地理空间数据云，空间分辨率为30 m×30 m。坡度数据由DEM数据在GIS空间分析工具得出，经济社会数据来源于中科院资源环境科学数据中心以及彬州市2009年和2019年统计年鉴。道路交通数据、规划数据及水域数据来源于《彬州市总体规划(2015—2030)》。

2 研究方法

2.1 数据处理

将收集的研究区2009年和2019年土地利用数据在GIS提取土地利用类型矢量数据，并分别将这两年的矢量数据投影成相同的坐标系统，利用GIS数据管理工具中Dissolve对其矢量数据进行融合，使其数据研究范围完全重合。根据《土地利用现状分类》(GB/T21010—2017)以及本次研究的重点，将土地利用类型划分为7类：耕地、园地、林地、草地、建设用地、水域和其他土地。DEM数据通过GIS软件进行镶嵌、整合、裁剪为范围一致的栅格数据，再根据GIS空间分析工具中Slope将DEM数据转换为坡度数据。经济社会数据主要包括GDP和人口密度，以栅格为单元进行处理，利用GIS软件进行裁剪，使其和研究区范围一致。

2.2 土地利用动态变化的空间分析模型

某一区域的某单一土地利用类型i在某段时期的转移速率和新增速率计算公式如下[20]：

(1)

(2)

某一区域综合土地利用动态度的计算公式如下：

(3)

式中:TRLi为转移速率;LA(i,t1)为研究初期t1第i种土地利用类型的面积；ULAi为研究期间第i种土地利用类型没有发生变化的面积；IRLi为新增速率；LA(i,t2)为研究末期t2第i种土地利用类型的面积；S为区域土地利用综合变化率；t1和t2分别为研究初期和末期的时间。面积单位均为km2。

2.3 Logistic回归模型及模拟有效性检验

Logistic回归模型可用于估算某个事件发生的可能性，并且能分析一个因变量与多个自变量之间的多元回归关系。本研究中，在考虑数据的可获取性、数据在时间空间上的一致性以及因子具有空间差异性的基础上，选取高程、坡度、GDP、人口密度、距道路、行政中心、河流的距离7种作为驱动因子进行回归分析。公式如下[21]：

(4)

式中:pi为出现某地类i的概率;β0为常数量;β1,β2,…,βm为回归系数;X1,X2,…,Xm为影响因子。pi取值范围在0～1,数值越大说明出现某地类的概率就越大,反之亦然。

将上式进行对数变换再推导得出土地利用的空间分布概率[22]：

(5)

式中:exp(β0+β1X1+β2X2+…+βmXm)叫做事件发生比,也称似然比,用exp(B)表示,在本研究中,exp(B)表示驱动因子每增加一个单位时土地利用类型发生比的变化情况:exp(B)>1,发生比增加;exp(B)<1,发生比减少;exp(B)=1,发生比不变。

对于Logistic回归模型的拟合效果常用ROC曲线[23]进行验证，一般认为当ROC曲线面积值大于0.7时，模拟效果较好，满足研究需求。

2.4 Logistic-CA-Markov模型构建

元胞自动机(CA)实质上是一种网格动力学模型[24]，其特点是时间、空间、状态都离散，根据空间相互作用和时间因果关系来模拟复杂系统的时空演变过程。CA模型表达式如下[25]：

St+1=f(St，N)

(6)

式中:S为元胞有限、离散的状态集合;N为元胞的邻域;t和t+1为不同时刻;f为元胞的转换规则。

马尔科夫模型(Markov)是离散时间随机模型，在土地利用预测领域运用比较多，它主要是根据过去和当前土地利用信息具体情况，计算出这两个时间段土地利用类型之间相互转移的面积数量和转移概率，进行下一时间点土地利用状态的预测。公式表达如下[26]：

St+1=Pij·St

(7)

式中:St,St+1为t,t+1时刻土地利用系统的状态;Pij为状态转移矩阵。

本研究CA-Markov模型构建是基于IDRISI 17.0软件的支持下完成的[27]，具体过程如下：(1) 数据转换和重分类；(2) 确定研究区2009—2019年土地利用类型的转移面积矩阵和转移概率矩阵；(3) 根据前文回归分析结果，运用地图代数根据公式(5)计算得到各土地利用类型概率图，并用IDRISI软件集合编辑器建立各土地利用类型的适宜性图集；(4) 构造CA滤波器，本研究采用5×5邻近滤波器为邻域定义；(5) 确定起始时刻及迭代次数。

本研究利用Logistic回归分析模型得出土地适宜性图集定义CA模型的转换规则，CA-Markov模型进行预测，三者结合形成Logistic-CA-Markov模型。

2.5 模拟预测精度检验

Kappa系数是一个用于一致性检验的指标，常用来检验模型预测结果和实际分类结果是否一致。其计算公式如下[28]：

(8)

式中:P0为预测一致的栅格占实际总栅格的比例；Pc为随机状态下模拟预测一致的栅格比例。当0

3 结果与分析

3.1 土地利用动态演变

通过ArcGIS 10.2对2009年和2019年研究区土地利用的矢量数据进行空间叠置分析，得到2009—2019年土地利用变化面积转移矩阵(表1)，依据该面积转移矩阵和公式(1)—(3)统计分析并计算各土地利用类型的转移速率，新增速率以及研究区总体土地利用动态度(表2)。

由表1可知，研究区2009—2019年10年间各个土地利用类型的面积发生了不同的变化，其中，林地、建设用地、园地和其他土地的面积分别增加了191.94，10.7，9.38，0.79 km2；草地、耕地和水域面积分别减少了176.47，35.69，0.64 km2。在各土地利用类型转移面积中，草地和耕地转移的面积最多，分别为232.29，190.61 km2，其转移速率为7.54%，5.26%；在各土地利用类型新增面积中，林地、耕地和园地新增面积最大，分别为266.12，154.92，104.83 km2，其中新增速率最快的是林地，草地转变为林地面积最多。研究区在2009—2019年这10年间的土地利用变化显著，再根据表2研究区土地类型在10年间发生转变的面积为653.31 km2，总面积为1 173.36 km2，依据公式(3)计算求得的研究区土地总体利用动态度为5.57%，可以看出研究区10年间土地利用类型的变化比较明显。

表1 彬州市2009-2019年土地利用变化的面积转移矩阵 km2

表2 彬州市2009-2019年土地利用动态变化率

3.2 研究区土地利用变化的驱动因子分析及有效性检验

将土地利用数据和各驱动因子的矢量数据转化ASCII文件，再利用CLUE软件将每种地类和所有的驱动因子ASCII文件转换成7个单一记录文件，每个记录文件中一个地类对应有7个驱动因子数据，最后利用SPSS软件对7个单一记录文件进行logistic回归分析，根据得出的各地类空间分布概率，运用地图代数方法制作7种地类的空间分布概率图(图1)，并利用ROC曲线对7种地类模拟有效性进行检验，结果见表3。

由表3可知：(1) 彬州市耕地格局主要受坡度、距道路、行政中心、河流的距离以及GDP影响，从回归系数看，坡度影响最大且与耕地呈负相关；从事件发生比exp(B)来看，坡度每增加1个单位，耕地的空间转化概率减少0.990 137倍。(2) 建设用地格局受自然因子和社会经济因子的共同影响，建设用地与人口密度呈正相关，与其他因子呈负相关，说明建设用地分布在海拔低、坡度低、距离道路、行政中心近的区域，从发生比看，人口密度每增加一个单位，建设用地的空间转化概率增加1.004 488倍，说明人口密度越大，建设用地的需求就会越大。(3) 从方程回归系数来看，园地与高程、坡度、GDP和距道路距离因子呈负相关，说明园地分布在高程坡度低的区域，且随着距道路距离和GDP的增大，园地呈减少的趋势；园地与人口密度、距行政中心距离和距河流距离3种因子呈正相关，且人口密度的相关系数最大，从exp(B)来看，人口密度每增加1个单位，园地的空间转化概率增加1.006 747倍。(4) 林地与GDP和人口密度因子呈负相关，与高程、坡度、距道路距离、距河流距离4种因子呈正相关，说明林地分布在人口较少，高程较高，坡度较大，距离河流和道路较远的区域，同时坡度每增加1个单位，林地的空间转化概率增加1.017 682倍。(5) 草地与高程、距河流距离两种因子呈负相关，与坡度、人口密度和距道路距离3种因子呈正相关，说明草地分布在高程低、人口密度大、距离河流近的区域，且随着坡度和距道路距离的增大，草地呈增大的趋势。水域与GDP呈正相关，与其他因子呈负相关。

图1 彬州市7种地类空间分布概率

在本文中，ROC值用于判断计算出的各类土地利用类型在不同因子影响下空间分布概率与该类土地利用类型真实的分布格局之间的拟合度[29]。由表3可知，各地类的ROC值依次是：耕地0.706，园地0.749，林地0.734，草地0.726，建设用地0.774，水域0.968，其他用地0.872，均大于0.7，说明拟合度较好，满足本研究需求。

3.3 彬州市2019年土地利用模拟及精度验证

根据7种地类空间分布概率图，运用IDRISI软件集合编辑器将7种地类概率图按重分类的顺序分别编序1—7，建立各土地利用类型的适宜性图集，在运用IDRISI软件中CA-Markov模型以2009年作为土地利用预测起始时间，设置迭代次数为10，模拟预测2019年土地利用变化情况，将其模拟结果与2019年实际土地利用现状图进行精度验证，见图2。

Kappa系数广泛用于对遥感数据的分类精度和两个图像的相似程度进行评价[30]。本研究通过IDRISI软件中CROSSTAB分析计算2019年彬州市土地利用预测图和土地利用现状图中各土地类型预测准确率(表4)，经统计得到彬州市总栅格数为1 301 880个，2019年预测模拟正确的栅格总数为1 035 403个，占彬州市总栅格数的79.53%，根据公式(8)计算得出2019年土地利用模拟的Kappa系数值为0.761 2，一般认为Kappa系数值大于0.75时，说明模拟效果为优，模拟精度通过检验。

图2 2019年彬州市土地利用现状、预测

3.4 多情景条件下彬州市2029年土地利用预测

本研究在自然发展、生态保护和限制城市过度开发3种情景下对彬州市2029年土地利用进行预测模拟，自然发展情景是按照当前的发展情况，不受任何条件的限制；生态保护情景是把水域和永久基本农田保护区作为约束条件；限制城市过度开发情景是根据彬州市政策把禁止建设区、限制建设区、有条件建设区和允许建设区作为约束条件。在IDRISI软件CA-Markov模块中以研究区2019年土地利用现状为基期，结合Logistic模型分析得出的土地适宜性图集和2009—2019年研究区土地利用面积转移矩阵作为转换规则，设置迭代次数10，分别在以上3种情景下对彬州市2029年土地利用变化情况进行预测，预测结果见图3。

表3 彬州市主要地类的Logistic回归模型

表4 2019年彬州市土地利用类型实际栅格与模拟栅格对比

在IDRISI软件中运用CROSSTAB计算3种情景下彬州市2029年土地利用预测的7种地类栅格数，以此为基础计算得出各土地类型占彬州市总面积的比例，再与2019年彬州市各土地利用类型的实际情况相对比，结果表明在自然发展情景下，2029年彬州市林地面积基本保持不变，耕地、园地和水域面积将分别减少27.31%，10.96%，3.84%，耕地减少主要分布于新民镇、城关镇、太峪镇和水口镇，园地减少主要集中在新民镇、北极镇和义门镇，水域面积减少集中在义门镇和龙高镇；草地、建设用地和其他土地面积分别增加24.35%，40.54%，21.09%，草地增加的面积主要集中在新民镇、义门镇南部、城关镇西部，建设用地增加的面积分布较均匀，增加的其他土地面积主要集中在龙高镇东部。

图3 3种情景下彬州市2029年土地利用预测

在生态保护情景下，2029年研究区园地面积和水域基本保持不变，草地、建设用地和其他土地面积分别增加41.33%，2.82%，10.62%，草地增加的面积集中在义门镇、北极镇、新民镇和龙高镇，建设用地增加的面积主要集中于城关镇，其他土地面积增加主要分布在新民镇；耕地和林地面积分别减少6.24%，12.84%，耕地面积减少主要集中在城关镇和新民镇两个镇，林地面积减少主要集中在永乐镇、北极镇和义门镇。

在限制城市过度开发情景下，2029年彬州市水域面积基本保持不变，耕地、园地、林地和建设用地面积分别减少13.1%，2.45%，12.05%，2.89%，耕地面积减少主要分布在城关镇、新民镇和北极镇，园地减少面积分布较均匀，林地面积减少主要集中在永乐镇、北极镇和义门镇，建设用地减少的面积主要集中于城关镇；草地和其他土地面积分别增加47.19%，11.38%，草地增加的面积主要分布于韩家镇、义门镇、北极镇和龙高镇，其他土地面积增加主要分布在新民镇。

4 讨论与结论

(1) 土地利用面积转移矩阵表明，彬州市在2009—2019年期间土地利用发生了很大的变化，7种地类面积的转化程度很显著，草地和耕地转移面积最多，具体有155.73 km2的草地转为林地，82.64 km2耕地转为林地，58.54 km2耕地转为园地以及25.39 km2耕地转为建设用地等，其中草地转化最明显；根据计算得出土地利用综合动态度为5.57%也可以看出研究区10年间土地利用有着显著性的变化。

(2) 用Logistic模型对选取的7种自然因子和社会经济因子进行驱动力分析，将得到的各土地类型的空间分布概率图组成土地适宜性图集，用以校正CA-Markov转化规则，并利用ROC曲线对7种地类模拟有效性进行检验，使其在数量和空间上达到较为精确的预测模拟，让预测结果更可信。

(3) 在自然发展情景下预测的2029年彬州市土地利用变化，建设用地面积增长过快，侵占了水域等重要生态用地，耕地面积也大量减少，长期发展下去可能会破坏研究区的生态平衡，这也不符合我国可持续发展观的要求；在生态保护和限制城市过度开发两种情景下，建设用地面积增加得到控制甚至还减小，水域生态用地也得到了保护，在一定程度上缓解了耕地保护的压力，有利于区域可持续发展。对比3种预测结果，生态保护和限制城市过度开发这两种情景下具有明显优势，可为未来彬州市土地规划、生态城市建设和区域生态环境保护提供科学的参考依据。

本文根据研究区2009年和2019年土地利用调查数据做多情景预测，虽然避免了常规的遥感解译产生的误差，精度更高，但是由于在数据获取上有局限，对影响研究区土地利用变化的驱动因子只从自然因素和社会经济因素选取了7个，考虑不够全面，使模拟的土地利用类型与该地类真实的分布格局之间的拟合度不够完全匹配，可能对最后预测精度产生一定的影响，因此，在下一步的研究中，我将对影响土地利用变化的政治制度因素、技术因素等进行全面考虑，选取因子。