基于改进CA-Markov模型的滹沱河流域山区段土地利用变化模拟及预测

席梅竹，赵中秋，吴攀升，赵颖丽

(1.中国地质大学(北京) 土地科学技术学院，北京 100083；2.忻州师范学院 地理系，山西 忻州 034000；3.自然资源部 土地整治重点实验室，北京 100035；4.山西省国土资源调查规划院，山西 太原 030024)

土地利用和覆被变化(land use and cover change，LUCC)是全球环境变化和可持续发展的重要组成部分，在全球物质循环和能量流动中起着关键作用，因此通过研究土地利用变化的规律，利用模型拟合并预测区域土地利用变化趋势对于合理开发利用土地资源，促进区域可持续发展具有重要意义[1]。由于土地利用和覆被变化是社会经济复合系统在时间和空间2个维度的反映，所以当前土地利用变化的研究热点已由对土地利用变化的数量预测转向时空模拟预测，同时预测模型也由单一模型向多模型耦合方向发展[2-3]。当前诸多学者已对土地利用变化时空动态预测模型进行了深入研究，主要包括CLUE-S模型(conversion of land use and its effects at small region extent)[4-5]、多智能体系统模型(multi-agent system,简称 MAS)[6-7]、元胞自动机模型(cellular automata,简称CA)[8-9]等。CLUE -S模型是基于统计的数量预测模型与基于经验的空间格局演化分配模型相结合，能够实现多尺度不同情景变化的模拟研究，但是由于其空间分配规则的主观性以及对土地利用数据空间相关性的忽视使模型的不确定性增加。MAS模型的核心思想就是由微观个体的相互作用产生宏观全局的变化[10]，然而由于土地利用决策系统的复杂性以及Agent的决策规则和行为模式的多样性，使大尺度宏观决策的模拟难度增加。CA模型是一种时空状态均离散、相互作用关系皆为局部的网格动力学模型，能够对复杂地理过程进行模拟和预测[11]。Markov 模型是基于马尔科夫过程系统而形成的一个过程理论模型，以达到预测和随机控制的目的[12]。结合两者优势的CA-Markov(时空马尔科夫)模型能够实现对时空演变的预测和模拟[13]。近年来，有不少学者应用CA-Markov模型进行区域土地利用的预测研究[14-16]，该方法虽具有较高的准确度，但是基于土地利用状态转移概率矩阵的土地适应性分析无法准确反映驱动因素与土地利用类型分布之间的数量关系，所以多因素驱动的土地利用类型变化逐渐引起学者的关注[17-18]。一些学者采用Logistic模型与CA-Markov模型相结合的方法，对区域土地利用/覆被类型进行预测研究[19-20]。但是应用CA-Markov模型进行预测研究的区域多集中于城市，对流域的研究较少，并且通常都是以10 a为周期，伴随着城市化进程的加快，利用5 a为研究周期可以对土地利用覆被变化进行更加细致的研究。

滹沱河流域山区段位于京津冀水源地的上游地区，土地利用/覆被影响下的区域生态环境状况会直接影响到下游华北地区的用水安全。同时该区域又属于忻州市的人口密集区和粮食主产区，随着经济社会的发展，区域人地矛盾突出。首先利用Markov模型探讨研究区2000-2015年的土地利用/覆被变化规律，并对影响区域土地利用/覆被类型分布的自然和社会经济因素进行Logistic回归分析，在此基础上借助CA-Markov模型以5 a为预测周期，对区域2020、2025年的土地利用/覆被格局进行预测，以期能够缓解区域人地矛盾，同时为制定更加科学合理的区域土地利用管理策略提供科学依据。

1 研究区概况

滹沱河发源于山西省忻州市繁峙县东北的泰戏山麓，属海河流域的子牙河水系。在忻州市境内向西南流经繁峙县、代县、原平市，至原平市与忻府区交界的界河铺村转向东流，流经定襄县和五台县，横穿太行山支脉系舟山，流入华北平原(图1)。滹沱河流域山区段处于滹沱河上游地区，地理位置介于38°13′-39°27′N，112°17′-113°58′E，研究区涵盖了6个县市区(图1)。境内支流呈羽状分布于山区，主要包括峨河、阳武河、云中河、牧马河、清水河、滹沱河支流，流域面积1 169 319.78 hm2。研究区属温带大陆性季风气候区，年平均气温4.0℃～8.8℃，多年平均降水量400～450 mm，且年内降水时空分布不均，7-8月降水量占全年降水量的53%；山区雨量大，平川、河谷雨量偏少。研究区土壤类型多样，山区以褐土、粗骨土和棕壤为主，平原河谷区以褐土和潮土为主，是忻州市的粮食主产区和人口聚集区。截至2018年底，研究区常住人口达210.46万，区域经济发展中第3产业占优势，3大产业结构比为7.5∶45.2∶47.3。

图1 研究区地理位置

2 数据来源与研究方法

2.1 数据来源及处理

研究区DEM数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn/)，空间分辨率为30 m。滹沱河流域山区段边界数据由SWAT水文模型计算得到的流域边界与已有行政区划边界相交得到。

土地利用遥感监测数据来源于美国地质勘探局官网(https://www.usgs.gov/)，轨道号为 125/33，选择下载2000年7月1日、2005年7月31日、2010年6月27日、2015年7月3日4期相近时段的Landsat TM/ETM/OLI影像，空间分辨率为30 m，首先经过影像融合、影像校正和图像增强等预处理；然后根据我国《土地利用现状分类》(GB/T 21010-2017)系统结合研究区实际，将土地利用类型分为6类，分别为耕地、林地、草地、建设用地(包括5～10类)、水域(水域及水利设施用地)、未利用地，并在此基础上通过监督分类中最优波段组合和人工目视解译结合，再经实地勘正校准，最终得到研究区4 a的土地利用类型栅格数据，2000、2005、2010、2015年研究区土地利用类型解译数据的Kappa系数分别为0.861 6、0.854 6、0.867 7、0.880 9，均满足遥感影像解译精度>0.85的要求。

研究区气温、降水气象数据来源于忻州市气象局6个气象站点的观测统计数据。社会经济数据来源于相应年份的《忻州市统计年鉴》。研究区的道路、河流数据由精度为15 m的遥感影像解译得到。研究区自然和社会经济空间数据均为30 m×30 m的栅格数据，由相应的属性数据利用ARCGIS10.5和IDRIS软件对其进行空间插值、属性连接、缓冲区分析、栅格重分类等处理得到。研究区所用数据统一采用WGS84坐标系、UTM投影。

2.2 研究方法

2.2.1 土地利用变化动态度 本研究综合采用单一土地利用动态度和综合土地利用动态度对滹沱河流域各种土地利用类型的数量变化及总体土地利用数量变化情况进行动态分析。其数学公式为：

(1)

式中，K—研究区某种土地利用类型的动态度；Ub、Ua—研究期末和期初某种土地利用类型数量；T—研究时长

(2)

式中，Lc—综合土地利用动态度；LUi—研究期初第i种土地利用类型的数量；ΔLUi-j—研究期间第i种土地利用类型转化为非i种土地利用类型面积的绝对值；T—研究时长。

2.2.2 土地利用变化图谱 土地利用图谱能够以空间单元和时间单元合并而成的图谱单元定量表示土地利用/覆被格局变化[21]。利用研究区三期遥感影像数据，在ARCGIS10.5中采用地图代数方法生成研究区2000-2005、2005-2010、2010-2015年土地利用格局变化图谱，能够准确反映研究区土地利用时空变化规律。

2.2.3 CA-Markov模型 Markov模型是对不同时期系统状态转移的定量描述，利用该模型对土地利用过程进行转换研究，可以获得不同时期土地利用状态之间的面积转移矩阵和面积转移概率矩阵[9]。其数学公式可描述为：

(3)

(4)

式中，i，j=1，2，…，n分别代表转移前后的土地利用类型；Sij—土地利用面积转移矩阵；Pij—土地利用面积转移概率矩阵。

CA模型中散布在规则格网中的每个元胞取有限的离散状态，遵循同样的作用规则、依据确定的局部规则作同步更新，其优点在于能够通过定义局部的元胞邻近关系，使用相对简单的局部转换规则，实现复杂系统的时空动态变化模拟，计算模型如下：

St+1=f(St,N)

(5)

式中，St+1—元胞在下一时刻的状态；St—元胞在上一时刻的状态；N—元胞邻域；f—局部空间元胞转换规则。

CA-Markov模型综合了CA 模型模拟复杂空间动态变化的能力和Markov 模型长期预测的优势[17]，在模拟土地利用格局的时空变化时，具有预测精度高、空间分布格局可视化的优点。

本研究利用IDRISI软件中的Markov模块和CA-Markov模块进行土地利用的模拟预测分析。首先用Markov模块计算研究区2000-2005、2005-2010年的土地利用面积转移矩阵，然后在CA-Markov模型中通过设置30 m×30 m元胞、5×5邻域结构、元胞自动机循环次数为10等参数，以2000-2005年的土地利用面积转移矩阵为条件1；以Logistic回归模型中得到的2005年土地利用适宜性图集为条件2；模拟得到2010年研究区土地利用类型分布图，最后运用IDRISI软件中的CROSSTAB模型对预测精度进行评价。预测精度满足要求的条件下，以2015年为预测基期，对2020、2025年的土地利用格局进行预测。

2.2.4 二元Logistic回归模型 Logistic回归模型是响应变量为事件发生概率的二值回归模型[14]，其回归模型为

(6)

(7)

式中，pi—区域内某种土地利用类型出现的概率；x—土地利用类型分布的驱动力因素；β0—回归模型常数项；β1～βm—回归系数。

二值Logistic回归模型的拟合优度通常采用ROC曲线线下面积值来评价，其取值范围介于0.5～1，取值越接近1，诊断效果越好，取值>0.75表示回归拟合结果准确性较高，能够满足要求[22]。

3 结果与分析

3.1 研究区土地利用/覆被变化特征

基于2000、2005、2010年和2015年的研究区土地利用/覆被数据，经过Markov模型分析得到研究区近20 a的土地利用面积转移矩阵，具体见表1～表4。

3.1.1 建设用地 2000-2015年，研究区土地利用动态度变化最大的为建设用地，面积占比呈不断上升的态势，由2000年的2.73%上升至2015年的3.44%，通过分析三个时段的转移矩阵可知，建设用地主要由耕地转入，随着城镇化进程的加快，城市快速扩张，城镇周边的耕地被划入规划范围，从而使耕地转入建设用地。同时在建设用地转为其他用地中，耕地和草地占比最大，区域内煤炭资源、有色金属资源开采之后的工矿用地、五台山核心景区企业村庄拆迁之后的工矿建设用地都是通过工程措施复垦为耕地和草地。

3.1.2 耕地 研究期间耕地面积比例呈缓慢下降趋势，面积比例由2000年的31.46%下降至2015年的31.13%，土地利用动态度由前期的-0.08%，中期-0.07%到后期的-0.06%，前期耕地主要转出为草地、建设用地和林地，转移面积分别为17 220.15、3 610.71、3 008.16 hm2，中期转移面积为17 371.17、3 615.48、2 816.01 hm2，后期转移面积为17094.24、3529.62、2989.17 hm2。根据土地利用转移图谱与坡度图的空间叠加分析可知，研究期间耕地转为草地的区域主要集中于坡度<10°的边缘地区，面积比例前中后期分别为84.56%、85.72%、84.48%，而且斑块空间分布破碎，其余转移区域集中于坡度为15°～25°的区域。耕地转移为林地的区域中，前期67.92%位于坡度<10°地区，中期比例上升到75.77% 后期降为67.87%，主要是由于经济林相较于传统种植业存在比较优势，当地政府增加农民收入、促进农业增长，积极进行农业结构调整由传统种植业向经济林果转变；其余集中于坡度15°区域，主要是为了实现区域生态环境的保护，积极落实国家对易造成水土流失的耕地实施退耕还林还草政策要求。耕地转移为建设用地的区域主要集中分布于河谷地区城镇周边坡度<5°的区域，主要是近20 a随着区域城镇化进程加快，城乡建设占用较多耕地；同期研究区域内其他地类也不断向耕地转移，转入耕地比例最大的是草地，主要集中于坡度<10°的区域。其次为林地和建设用地，分布于代县、原平、忻府区境内河谷地带，主要是由于区域废弃工矿用地等建设用地通过土地整治工程转入。

3.1.3 林地 研究期间林地面积比例呈不断下降趋势，面积比例由2000年的23%下降到2015年面积22.84%，而且土地利用动态度呈现不断增加趋势，由前期的-0.03%增加到中期-0.04%，后期增为-0.07%。林地在15 a间主要转出为草地和耕地，前、中、后3期林地转为草地的面积分别为12 009.60、11 631.87 hm2和11 954.52 hm2，空间集中分布于研究区东西两侧的山地地区，其中坡度为10°～15°的区域占比最大，主要原因是陡坡地区黄土沙砾沉积物极易产生水土流失导致土壤贫瘠化，进而出现林地的退化。前、中、后3期林地转为耕地的面积分别为2 754.99、2 998.71 hm2和2 770.02 hm2，其中坡度为0°～5°的区域占比最大；空间分布于忻府区南部地区及滹沱河河谷原平市、代县一带。

3.1.4 草地 草地是区域土地利用格局的基质景观。研究期间，草地面积比例出现不断下降的趋势，由2000年的41.31%下降到2015年的41.07%，中后期土地利用动态度明显高于前期土地利用动态度，草地转出为耕地、林地的同时，也有大量的耕地和林地转入草地。研究期间，草地主要转出为林地的区域中，主要是由于农户对宜林荒山、荒坡的承包种植以及区域生态环境不断改善在一定程度上促进草地向林地的转移。同时，耕地、林地及建设用地也在向草地转移。

3.1.5 水域 水域面积基本保持稳定略有上升趋势，主要是由于近年来滹沱河两岸绿色廊道建设及城镇周边水上公园、水利风景区的建设增加了河流的积蓄水面积，面积比例由2000年的1.38%上升为后期的1.39%。

3.1.6 未利用地 研究区未利用地面积所占比例最少，面积比例由2000年的0.12%上升为2015年的0.13%。后期增加的未利用地主要由林地和草地转入，空间分布最大的区域位于原平市境内西侧的坡度为10°～15°区域，说明区域内均存在不同程度的林地和草地的退化现象。

表1 2000-2015年研究区各类土地利用动态变化

表2 2000-2005年研究区土地利用变化面积转移矩阵

表3 2005-2010年研究区土地利用变化面积转移矩阵

3.2 影响土地利用格局变化的驱动因素分析

利用IDRISI软件中Logistic回归分析得到气温、降水、海拔、坡度、坡向5项自然因素，人口密度、地均GDP、距居民点距离、距河流距离、距路网距离5项社会因素与6种土地利用类型空间分布之间的回归方程，关系如下：

由表5标准化回归系数可知，耕地分布的关键影响因素包括正向因素——人口密度以及负向因素——坡度、高程、距道路距离，即坡度和高程值越大，耕地分布的可能性越小，距离道路距离越近，耕地分布的可能性越大，由于山区谷地水热资源条件好，人口分布较集中，农业生产耕地面积大；影响林地分布的关键因素包括正向因素——高程以及负向因素——人口密度，山区林地在海拔高的区域分布概率大，人口分布稠密的地区受人类活动的影响，经济收益相对较少的林地分布概率小；影响草地分布的关键因素包括负向因素——人口密度和正向因素——降水，降水量大和人类活动干预少的地区草地分布概率高；影响水域分布的关键因素为负向因素——坡度和距道路距离，地势越低，距离道路越近，水域分布的可能性越大；影响建设用地分布的关键因素是负向因素——距道路距离、距居民点距离，因为人口聚居、交通便捷的区域建设用地分布的可能性越大；影响未利用地分布的关键因素为负向因素——距居民点距离，距离居民点距离越大，未利用地分布概率越小，因为区域未利用地主要分布在人类活动破坏土地之后导致无法自我恢复的区域。

Logistic回归结果拟合精度评价中，耕地、林地、草地、水域、建设用地和未利用地的ROC曲线线下面积值分别为0.935 3、0.887 8、0.860 2、0.877 5、0.942 0、0.765 9，均>0.75，表明各种土地利用类型与驱动因素之间建立的回归模型拟合精度均满足要求。

表4 2010-2015年研究区土地利用变化面积转移矩阵

表5 各因素对土地利用类型分布的标准化回归系数

3.3 土地利用/覆被预测精度检验

目前CA-Markov模型预测结果没有统一的精度检验标准，本研究将预测结果栅格数据与2010年土地利用/覆被现状数据进行比较(表6)，预测精度通过IDRISI软件中的CROSSTAB模型进行评价，依据Cohen提出的模拟预测精度评价标准[23]，Kappa系数介于0.6～0.8，预测结果一致性程度属于显著，而本研究Kappa系数为0.825 1，说明改进的CA-Markov模型模拟滹沱河流域山区段土地利用变化趋势具有较高精度。

3.4 土地利用/覆被预测

利用通过精度验证的Logistic-CA-Markov模型对2020、2025年研究区的土地利用格局进行预测，预测结果显示：在自然环境因素不发生突变，社会发展保持现有速度的条件下，到2020年研究区各种用地面积分别为耕地 36.43万hm2、林地 26.71万hm2、草地 48.0 3万hm2、水域 1.62万hm2、建设用地3.99万hm2、未利用地0.1 5万hm2(表7)。通过与2000-2015年的各类用地的变化趋势对比发现，区域土地利用格局基本与前一时期保持一致，2015-2020年耕地面积缓慢增加，土地利用动态度由2010-2015年的-0.06%升高为0.01%，在转入耕地的所有地类中，草地转移面积最大为2.48万hm2，空间集中分布于北部和西南部河谷与山地的过度地带(图5)；建设用地面积减少，土地利用动态度由2010-2015年的1.64%降低到-0.15%，林地、草地和水域面积基本保持稳定，未利用地面积有小幅下降。2025年研究区各种用地面积分别为耕地36.52万hm2、林地26.18 万hm2、草地47.78万hm2、水域1.95万hm2、建设用地4.32 万hm2、未利用地面积分别为0.18万hm2，林地和草地面积继续呈减少趋势，耕地和建设用地面积继续增加，水域面积有所增加。

表6 2010年研究区土地利用类型实际栅格数与预测栅格数对照

表7 2020、2025研究区土地利用变化预测

图2 2000-2015年研究区土地利用/覆被变化图谱

图3 2020、2025年研究区土地利用/覆被预测

4 讨论

山区土地利用/覆被变化是一个非常复杂的过程，虽然本研究中对CA-MARKOV模型中土地适宜性图集进行了改进，但是在土地利用政策不变的前提下选取了30 m×30 m的尺度进行区域土地利用预测研究，未来土地利用/覆被变化预测研究中应该加入政策因素的驱动以及进行不同时间、空间尺度优选，从而使研究结果具有更高实效性。当前预测结果可以在一定程度上为区域土地利用管理提供参考。随着城镇化进程的加快，政策积极引导的重要性不容忽视。相关部门应加快推进土地流转，在保证粮食安全的前提下，推广种植优质高产经济林果作物从而提高农民收入、落实就地城镇化政策。同时，林业部门应对林地退化区进行深入研究，优选适生森林植被保护区域水土资源，政府应加大对废弃工矿用地以及实施搬迁后废弃的村庄整理复垦的力度，补充耕地数量，统筹推进区域山水田林湖草的综合治理。

5 结论与讨论

2000-2015年研究区草地占主体地位，面积逐年减少，林地和耕地也呈减少趋势，耕地主要转变为草地，在一定程度上说明区域退耕还草生态工程成效显著，同时也表明区域存在弃耕现象。林、草地之间转换频繁，而且林地转草地的总面积大于草地转为林地的面积；由于区域生态环境脆弱，今后应加强对林地的保护，防治林地退化。建设用地和未利用地面积呈不断增加的趋势，水域面积有小幅增加的趋势。

2000-2015年研究区综合土地利用动态度0.69%，整体土地利用变化速率较低，但各种土地利用类型变化速率差异较大，其中建设用地变化速率最大，土地利用变化动态度由1.33%上升为1.64%，林地和草地土地利用动态度呈现不断上升趋势，说明林、草地面积减少的速度在增加，在退耕还林还草的同时，应加强对林草地的管护。同时要严格限制工矿用地等对林、草地的占用。耕地的土地利用动态度逐年降低，说明耕地减少的速度在降低。

2020、2025年的土地利用/覆被格局预测结果显示，2025年研究区土地覆被将会继续保持建设用地增加的趋势，但建设用地增加速率降低，同时林地面积有所减少，会转化为草地，由于当地自然条件及人类活动的影响导致林地植被出现不同程度的退化趋势。草地面积出现下降趋势；耕地出现增加趋势；水域出现增加趋势，反映了滹沱河流域上游湿地系统的保护成效。