城市化进程在我国过去的20年里是最大的社会经济空间运动，在未来的20年里也将受到重点关注^[1]。城市化推动社会经济的发展，与此同时带来的生态问题也不容忽视，平衡经济发展和生态保护这两者的关系，成为我国全面发展的重大挑战^[2]。随着国家政策的实施，大中城市城市化发展迅速，受到较多专家学者关注，但对于县域尺度城镇化问题的研究并不多见^[3]。近些年来，县域城市化发展也是突飞猛进，城镇用地扩张带动当地经济发展，对当地的生态环境也产生了一定影响，故探索县域城镇化的形成机制，对于确立县域城镇化的发展战略、因地制宜地引导城镇化的地域模式和用地政策提供参考依据^[4]。

近年来，元胞自动机(cellular automata，CA)被广泛应用于地理空间模拟研究，与Markov、多智能体、Logistic回归等多种模型进行结合用以提高模拟精度^[5]，其中以CA-Markov模型研究较为热点。CA-Markov模型既具有CA模型模拟复杂系统空间变化的能力，也具有Markov模型长期预测的优势^[6]，且仍然不断被改进。李俊等^[7]通过Logistic改进CA-Markov模型对宁蒙沿黄地带城镇用地进行扩展驱动力分析和情景模拟，张晓娟等^[8]基于多标准评价(multi-criterium evaluation, MCE)改进CA-Markov模型对三峡库区土地利用演变进行模拟预测，均取得较好结果。作者综合前人研究，引入地理探测器(GeoDetector，GD)和支持向量机(Support Vector Machines，SVM)2种方法，从原理上对CA模型进行改进，并结合Markov模型实现模拟预测。地理探测器在分析土地利用变化与驱动因子关系方面具有较大优势^[9]，用以CA模型土地利用适宜性参数的设置；支持向量机在特征提取和目标识别方面具有广泛的应用^[10]，用来确定CA模型非线性转换规则；由此建立GD_SVM_CA-Markov模型。

研究区内蒙古磴口县是典型的干旱半干旱区域，沙漠化问题显著，生态环境保护与治理亟需解决^[11]，且处在城镇化的进程当中，平衡好当地建设用地、沙漠化土地以及绿色生态用地这3类用地景观显得尤为重要。笔者基于GD_SVM_CA-Markov模型对磴口县从时间和空间2个维度对其用地景观动态分布变化进行分析，探究其变化规律，并进行模拟预测，从而希望能为当地城镇发展规划、沙漠化治理以及生态环境保护提供一定的决策参考。

磴口县地处中国西北(E107°05′，N40°13′)，位于内蒙古河套平原源头，黄河中上游，靠狼山山脉，西邻乌兰布和沙漠，如图 1所示。磴口县属中温带大陆性季风气候，气候干旱少雨，区域蒸发量大，水资源较为短缺，土地沙漠化严重。多年平均降水143.9 mm，多年平均蒸发量2 327 mm，多年平均气温7.6℃。黄河流经磴口县52 km，年径流量310亿m³，共有水域面积2 406.67 hm²。境内海拔1 030~2 046 m，整个地形除山区外，呈东南高西北低，逐渐倾斜^[12]。

图 1(Figure 1)

图 1 磴口县区位图 Figure 1 Location map of Dengkou county

本研究所用原始数据主要包括磴口县2006年、2011和2016年Landsat TM/ETM+/OLI分辨率为30 m的遥感影像、磴口县分辨率为30 m的DEM数据以及自然、社会经济统计数据等。对原始数据进行处理后可分为土地利用数据和驱动因子数据2类。利用ENVI和ArcGIS软件对3期的遥感影像进行监督分类和人工目视解译，根据研究区特点，分为建设用地、水域、耕地、林草地、沙地、未利用地6类，分类精度达到了85.6%。驱动因子数据一部分来自于磴口县相关水利部门的测点数据以及磴口县各年社会经济统计年鉴，一部分是通过在ArcGIS、ENVI软件中分析处理得到，可分为自然驱动因子数据、社会经济驱动因子数据和邻域驱动因子数据3类，共10种。

基于GD_SVM_CA-Markov模型对研究区景观格局变化进行模拟，具体的过程如下：

1) 适宜性图集的创建。将土地景观类型发生变化的区域设置为1，土地景观类型未发生变化的区域设置为0，选取DEM、坡度、坡向、NDVI、地下水埋深、蒸散发、人口密度、距水域最近距离、距居民点最近距离、距道路最近距离共10个驱动因子，将土地景观类型变化值(1或0)与10个驱动因子通过地理探测器(GeoDetector，GD)在GeoDetector程序中进行因子探测分析，地理探测器因子探测的原理^[9]为

式中：P_{D, H}为驱动因素解释力；D为地理要素指标变化的驱动因子；n与σ²分别为研究区域整体样本量与方差；H为某种因子的分类数；n_{D, i}为D指标在i类上样本数。P_{D, H}取值区间为[0, 1]，其数值越大则表明该因子对地理要素指标变化解释能力越强^[13]。

得到的结果作为接下来CA转化规则的参数，并结合IDRISI软件提供的MCE模块完成土地利用适宜性图集创建。

2) 定义元胞、状态、邻域和转化规则。笔者利用支持向量机(SVM)制定分类决策^[14]，利用非线性的转换规则来反映土地景观格局系统复杂的特征^[15]，转化规则的非线性不确定性定义主要由支持向量机的设计与开发、训练数据的采样、模型检验等组成，将CA模型与支持向量机(SVM)分类的土地景观格局发展概率结合，通过计算得出各区域空间变量对每个元胞的土地景观格局发展概率的贡献。以CA标准的5×5邻近滤波器为邻域定义，即每个中心元胞周围有5×5元胞组成的矩阵空间对该元胞状态的改变具有显著的影响^[16]。

3) Markov过程。基于2006年和2011年2期的土地景观数据，应用Markov模型生成土地利用转移(概率)矩阵^[17]，作为接下来模拟的各土地景观数量需求。

4) 实现模拟。以2011年作为模拟的起始时刻，加入适宜性图集，利用已训练好的SVM模型，结合经过Markov过程的2006—2011年土地利用转移矩阵以及转移概率矩阵，CA的迭代次数设置为10，对研究区土地景观格局发展趋势进行外推，实现2016年研究区土地景观空间格局预测模拟^[18]。

5) 精度验证。为了检验GD_SVM_CA-Markov模型模拟2016年不同土地景观类型面积和空间分布格局结果的合理性，笔者运用IDRISI软件中的Crosstab模块，引入各类景观的Kappa系数^[19]，将2016年的土地景观模拟结果与实际情况进行比较。

如表 1所示，2006—2016年间，由于城镇扩张，磴口县建设用地景观面积呈现持续增加的趋势，且2011—2016年间增加幅度更大；沙地景观面积呈现持续减少趋势，可见当地沙漠化治理取得一定成效，但在2011—2016年间，沙地面积减少速率放缓，沙漠化治理仍需加大力度；耕地景观面积呈现持续增加的趋势，未利用地景观面积呈现持续减少的趋势；受到城镇扩张的影响，水域、林草地等生态用地景观面积呈现前期增加缓慢、后期减少的趋势，城镇发展的同时，应注意生态环境保护，退耕还林(草)等生态保护工程仍需抓紧。

表 1(Table 1)

表 1 2006—2016年不同土地景观类型面积变化 Table 1 Area changes of different landscape types from 2006—2016 hm2 年份 Year 建设用地 C1 水域 C2 耕地 C3 林草地 C4 未利用地 C5 沙地 C6 2006 5 780.55 23 859.88 60 806.47 117 452.37 82 898.15 76 616.15 2011 6 259.79 26 300.71 70 150.88 129 117.65 81 679.31 57 970.46 2016 8 952.67 24 679.10 81 964.14 128 120.87 73 094.66 56 460.50   Notes: C1:Construction land. C2: Water area. 3: Arable land. 4: Woodland and grassland. C5: Unused land. C6: Sandy land. The same below.

表 1 2006—2016年不同土地景观类型面积变化 Table 1 Area changes of different landscape types from 2006—2016

如表 2、表 3所示，2006—2011年间，土地利用转移过程主要发生在建设用地—耕地、水域—耕地、沙地—林草地之间的转化。就建设用地而言，有15.64%的概率转化为耕地；就水域而言，有11.56%的概率转化为耕地；就沙地而言，有18.37%的概率转化为林草地。建设用地―沙地、水域—沙地、未利用地—沙地以及沙地—建设用地之间基本上未发生转化。

表 2(Table 2)

表 2 2006—2011年土地利用转移矩阵 Table 2 Land use transfer matrix from 2006 to 2011 hm2 土地利用 Land use C1 C2 C3 C4 C5 C6 C1 3 497.75 18.34 707.55 278.93 19.20 0.91 C2 151.29 16 419.75 2 446.92 1 868.54 287.63 0.00 C3 1 297.49 3 670.96 50 072.51 1 842.30 3 897.83 248.48 C4 136.58 2 190.03 9 418.54 103 763.59 1 899.66 155.48 C5 1 170.27 1 906.73 3 338.37 6 324.26 75 022.93 0.00 C6 6.43 94.89 4 166.98 14 040.03 552.06 57 565.58

表 2 2006—2011年土地利用转移矩阵 Table 2 Land use transfer matrix from 2006 to 2011

表 3(Table 3)

表 3 2006—2011年土地利用转移概率矩阵 Table 3 Land use transfer probability matrix from 2006 to 2011 土地利用 Land use C1 C2 C3 C4 C5 C6 C1 0.773 378 0.004 056 0.156 445 0.061 674 0.004 245 0.000 202 C2 0.007 145 0.775 463 0.115 562 0.088 246 0.013 584 0.000 000 C3 0.021 260 0.060 151 0.820 463 0.030 187 0.063 868 0.004 072 C4 0.001 162 0.018 628 0.080 114 0.882 615 0.016 159 0.001 323 C5 0.013 334 0.021 726 0.038 039 0.072 061 0.854 840 0.000 000 C6 0.000 084 0.001 242 0.054 523 0.183 708 0.007 223 0.753 220

表 3 2006—2011年土地利用转移概率矩阵 Table 3 Land use transfer probability matrix from 2006 to 2011

将研究区土地景观类型变化值及10种驱动因子带入GeoDetector程序进行分析计算，得到如表 4所示的结果，可知，P值为驱动因子对土地景观类型变化的解释程度，P值大小与解释程度强弱成正比。10种驱动因子对磴口县土地景观类型变化影响程度的大小依次是：地下水埋深、NDVI、蒸散发、人口密度、距居民点最近距离、距水域最近距离、DEM、距道路最近距离、坡向、坡度，其中地下水埋深、NDVI以及蒸散发对磴口县土地景观类型变化的解释力极强。由此可见，水文条件、植被覆盖以及人口要素是决定当地土地景观类型时空变化的决定性因素。

表 4(Table 4)

表 4 驱动因子地理探测结果 Table 4 GeoDetector results of driving factors DEM 坡度 Slope/(°) 坡向 Aspect NDVI 地下水埋深 Groundwater depth/m 蒸散发 Evapotranspiration/(mm·d-1) 人口密度 Population densityPeople/km2 距水域最近距离 The nearest distance to water area/km 距居民点最近距离 The nearest distance to settlement/km 距道路最近距离 The nearest distance to road/km P 0.248 816 0.048 784 0.134 342 0.951 212 0.975 924 0.873 667 0.520 317 0.256 626 0.413 550 0.178 658

表 4 驱动因子地理探测结果 Table 4 GeoDetector results of driving factors

为了验证模型的精度，本研究与传统CA-Markov模型的模拟结果进行精度比较，来说明GD_SVM_CA-Markov模型的优势性。2016年真实解译图、2016年GD_SVM_CA-Markov模拟图、2016年CA-Markov模拟图如图 2所示。经计算(表 5)，2016年CA-Markov模型模拟结果的Kappa系数为0.862 8，2016年GD_SVM_CA-Markov模型模拟结果的Kappa系数为0.925 0，总体模拟精度均较高，但GD_SVM_CA-Markov模型模拟精度更高。因此，利用GD_SVM_CA-Markov模型对磴口县土地景观空间分布格局进行模拟预测，模拟精度较高，结果可靠。

图 2(Figure 2)

图 2 磴口县土地景观解译数据与模拟结果 Figure 2 Landscape interpretation data and simulation results in Dengkou county

表 5(Table 5)

表 5 模型精度对比 Table 5 Model accuracy comparison 土地景观 Landscape 2016实际栅格数 Actual raster number in 2016 2016 GD_SVM_CA-Markov模拟栅格数 Simulated GD_SVM_CA-Markov raster number in 2016 Kappa系数 Kappa coefficient 2016 CA-Markov模拟栅格数 Simulated CA-Markov raster number Kappa系数 Kappa coefficient 总体 Total 20 713 446 21 564 053 0.925 22 252 184 0.862 8

表 5 模型精度对比 Table 5 Model accuracy comparison

基于2016年磴口县土地利用数据以及2011—2016年土地利用转移数据，应用GD_SVM_CA-Markov模型，对2021年土地景观空间分布格局进行模拟预测，模拟结果如图 3所示，可以直观看出，建设用地景观和沙地景观面积增多，水域、林草地等生态用地景观面积减少。在ArcGIS软件中进行面积统计分析后可知，建设用地景观面积由2016年的8 952.67 hm²增加为1万1 619.21 hm²，沙地景观面积由2016年的5万6 460.50 hm²增加为6万7 235.11 hm²，水域、林草地等生态用地景观面积由2016年的15万2 799.97 hm²减少为14万3 670.04 hm²。北部以及东北部的部分村镇出现大量城镇扩张，那仁布鲁格嘎查西南部以及沙拉毛道嘎查西部出现沙漠化扩张现象，对于磴口县而言，在城镇发展的同时，采取沙漠化治理以及生态环境保护措施刻不容缓。

图 3(Figure 3)

图 3 基于GD_SVM_CA-Markov模型的2021年磴口县土地景观格局模拟结果 Figure 3 Simulated landscape pattern in Dengkou county in 2021 based on GD_SVM_CA-Markov model

2006—2016年间，磴口县建设用地景观面积呈现持续增加的趋势，受到城镇扩张的影响，当地水域、林草地等生态用地景观面积呈现前期增加缓慢、后期减少的趋势，沙地景观面积呈现持续减少趋势，可见当地沙漠化治理取得一定成效。水文条件、植被覆盖以及人口要素是决定当地土地景观类型时空变化的决定性因素。故城镇发展的同时，应注意生态环境保护，退耕还林(草)等生态保护工程仍需抓紧。

2016年CA-Markov模型模拟结果的Kappa系数为0.862 8，2016年GD_SVM_CA-Markov模型模拟结果的Kappa系数为0.925 0，总体模拟精度均较高，但总体上GD_SVM_CA-Markov模型模拟精度更高。因此，利用GD_SVM_CA-Markov模型对磴口县土地景观空间分布格局进行模拟预测，模拟精度较高，结果可靠。

本研究中，模拟预测的过程统一采用了30 m×30 m大小的栅格尺寸，滤波器大小设置为5×5，对于尺度的考虑不足，未来会针对不同元胞大小和滤波器大小进行模拟预测的研究，探究最适合的元胞大小和滤波器大小，以期达到更高的模拟精度和结果。社会经济要素对于土地利用变化具有一定影响，本研究对于社会经济要素涉及不多，未来会更多地参考地区相关的社会经济政策。另外，接下来也会就当下热点的多情景模拟问题继续深入研究，通过设定一定情景，来定义元胞的转化规则，实现多情景模拟，提高研究的丰富性、全面性。