基于CA-Markov模型的哈尔滨市土地利用变化预测研究*

李志明，宋 戈.2※，鲁 帅，王 蓓，徐四桂

(1.东北农业大学资源与环境学院，黑龙江哈尔滨 150030；2.东北大学土地管理研究所，辽宁沈阳 110819)

0 引言

土地利用变化是人类经济社会发展对土地资源利用的直接反映，是全球环境变化和可持续发展研究的主要内容[1-2]。加速的土地利用变化产生更加复杂的土地利用格局，直接或间接引发生态环境变化[3]，如土壤质量[4-5]、气候条件[6-7]、水文状况[8]等。因此，科学合理地分析土地利用演变特征，探究未来土地利用变化趋势，能够为区域土地资源的合理开发利用、生态环境的改善和保护、经济的可持续发展提供科学的参考依据[9-10]。目前，国内外土地利用变化预测模型主要包括侧重时间维度分析的系统动力学模型、Markov模型、Logistic模型、神经网络模型等和侧重空间维度分析的CA模型、CLUE模型、CLUE-S模型等。CA-Markov模型由于结合了Markov模型的时间维度分析优势和CA的空间维度分析能力[11]，得到了广泛的应用和探讨。研究尺度侧重于中尺度的流域[11-13]和城市[14-15]，研究内容集中于模型指标选取对预测精度的影响[16]和未来土地利用趋势对区域生态环境的影响[13]等方面，但多数研究中将Markov模型预测的条件概率图像、较单一的影响因子或单生态方面的影响因子作为元胞自动机的转换规则，转换规则过于简单，缺乏社会经济因子对土地利用变化影响的考虑，使结果不能更好地贴合实际发展情况。因此，文章以哈尔滨市建成区为例，选取经济和生态两方面的GDP、坡度、水域、与道路距离、与水域距离、与建设用地距离共6种影响因子，采用MCE模型中的加权线性合并法单独制作耕地和建设用地适宜性图集参与到元胞自动机转换规则中。应用GIS技术、Markov模型和土地利用重心迁移模型探究1992～2014年土地利用时空演变规律，最后基于CA-Markov模型对2025年研究区土地利用变化进行更加细致的预测。

1 研究区概况和数据来源

1.1 研究区概况

哈尔滨市位于黑龙江省南部，地跨125°42′E～130°10′E，44°04′N～46°40′N，地势东南高，西北低，东南地区多山地丘陵，西北地区位于松嫩高平原上，平坦低洼，松花江贯穿其中。年平均降水570mm，四季分明，属于温带大陆性季风气候。全市总面积5.3万km2， 2014年总人口达987万，是重要的工业城市，城镇化发展迅速。该文选取哈尔滨市的道外区、道里区、南岗区、阿城区、松北区、呼兰区、香坊区和平房区8个主要市辖区为研究区。该区域位于哈尔滨市西北部，总面积约7 083km2，占哈尔滨市总面积的13.36%，土地利用类型丰富，以耕地和林地为主， 2014年土地利用情况与1992相比，耕地面积减少了10.91%，建设用地增加了74.37%，是研究土地利用变化的热点区域。

1.2 数据来源

选取研究区1992年、2003年和2014年3期遥感影像(1992年和2003年为Landsat 5 TM影像， 2014年为Landsat 8 OLI影像，分辨率均为30m)为基础数据(数据下载于http：//www.gscloud.cn/)，利用ENVI 5.1对遥感影像进行波段融合、图像校正等预处理后，应用ArcGIS10.2对数据进行解译，解译过程借助地形图等相关材料进行校准，参考《第二次全国土地调查土地分类》标准和研究区实际情况，将土地利用类型分为耕地、林地、草地、水域、建设用地和未利用地6个一级分类(图1)。其他数据包括研究区30m分辨率的DEM数据，哈尔滨市社会经济统计数据等。

2 研究方法

基于研究区3期土地利用解译数据，运用GIS空间分析技术，利用Markov模型和土地利用重心迁移模型分析1992～2014年土地利用动态变化特征，选取一定约束条件和适宜性因子，结合MCE模型，建立适宜性图集，采用CA-Markov模型对土地利用变化进行预测研究，并应用土地利用动态度对预测结果进行评价。

2.1 土地利用重心迁移模型

土地利用重心是在人口分布重心理论的背景下提出，用于表征土地开发利用过程的重心空间位置以及变化情况。计算方法如下：

(1)

式(1)中，Xk、Yk分别表示k种土地利用类型重心X、Y坐标值；Aki表示k种土地利用类型i图斑的面积；Xki、Yki分别表示k种土地利用类型i图斑的几何中心X、Y坐标值；Ak表示k种土地利用类型总面积。

2.2 MCE模型

MCE是为实现一个具体目标，而设定一系列的评价准则的过程[17]，其最基本的问题是如何整合不同的决策准则信息，即如何整合约束条件和适宜性因子两类决策准则信息。该文采用MCE中的加权线性合并法，能够较好的表达连续性的适宜性因子，其运算过程如下：

S=∑wixi∏cj

(2)

式(2)中，S为适宜性指数；wi为第i个适宜性因子的权重指数；xi是第i个适宜性因子的得分。该文将适宜性因子统一模糊处理至0～255范围内，数值越大，适宜度越高；cj是约束条件j的数值，其值为0或1，表示不流转或可流转。

结合研究区经济、生态条件，选取GDP、坡度、水域、与道路距离、与水域距离、与建设用地距离6种影响因子，针对建设用地和耕地，单独设置约束条件和适宜性因子用于制作适宜性图集，其余土地利用类型由于限制因素较少，故选择Markov模型生成的条件概率矩阵作为适宜性图集。具体如下：

(1)建设用地。设定约束条件为： ①水域，考虑到水资源短缺、水域污染等生态环境问题，出于对研究区的环境保护，将水域设定为不发生用地类型转换，赋值为0，其余赋值为1； ②坡度，根据《水土保持工作条例》和研究区自然条件，设定坡度大于25°为禁止开发区域，赋值为0，其余赋值为1。设定适宜性因子为： ①与建设用地距离，建设用地对周边土地的转入有较高的吸引，且随着距离的增大，吸引力不断弱化，直至形成无差别影响。故设定500m内适宜流转，适宜度随距离的增加单调递[12]； ②与水域距离，出于对环境的保护，不考虑距水域100m范围内土地的流转， 100m外随距离增加，转变为建设用地适宜度增高[13]； ③与道路距离，道路对经济发展的影响显著，设定研究区主要道路100m内，适宜度不断减弱， 100m外适宜度无差别[10]； ④GDP，由于建设用地对GDP贡献最大，故设定GDP越高，建设用地转入适宜度越高。

(2)耕地。设定的约束条件与建设用地相同。设定适宜性因子为： ①与建设用地距离， 500m内适宜性最弱， 500m外适宜度不断增加； ②与水域距离，考虑到水源对耕地的重要性，设置100m内，适宜度最高， 100m外适宜度逐渐减弱； ③与道路距离，与主要道路100m内适宜度最弱， 100m外适宜度逐渐增强； ④GDP，耕地对GDP的贡献相对较小，故设定适宜度随GDP的降低而降低。

2.3 CA-Markov模型及预测方法

Markov模型是用不同时间的系统状态给出一个从一种土地利用类型到另一种类型的转换概率矩阵，预测未来某一时刻的系统状态，是一种空间概率模型，强调的是从时间维度上分析土地利用的变化趋势。CA模型是基于元胞的前一个状态及周边邻域状态按一定转换规则预测未来元胞状态的模型，是局部网络动力学模型，其主要强调在空间维度上进行预测。CA-Markov模型结合Markov模型时间维度分析的优势和CA模型在空间维度分析的优势，降低了制定转换规则的难度，减少了人为因素的干扰。具体预测过程如下：

(1)计算土地转移矩阵。将1992年和2003年土地利用解译数据进行叠置分析，利用Markov模型计算土地转移概率矩阵和面积矩阵，并作为转换规则参与后面预测。

(2)适宜性图集的建立。应用MCE整合不同准则标准，得到不同土地类型的转移适宜性图，合并成适宜性图集，参与土地利用变化的预测。

(3)参数设定。该文以2003年为预测时间起点，设定5×5的元胞滤波器，对2014年土地利用空间分布进行预测。在预测精度较好的情况下，以2014年为起始年，进行2025年的土地利用空间分布预测。

2.4 土地利用动态度

土地利用动态度是反映土地利用变化程度的重要指标，其动态模型可分为单一和综合土地利用动态度，该文选择单一土地利用动态度指标进行预测结果的评价，可以反映某一土地利用类型变化的剧烈程度。

(3)

式(3)中，i表示土地利用类型；Ki表示一定时间段内i土地利用类型的动态度；At2和At1分别表示变化目标年和变化基准年土地利用面积；t2-t1表示研究期时间间隔。

3 结果分析

3.1 土地利用结构变化分析

基于ArcGIS10.2分析了1992、2003和2014年3期解译的土地利用数据，得到土地利用情况(表1)。

由表1可见，研究区内耕地占主体地位，所占比例分别为72.35%、70.45%和64.45%，在此期间，耕地以每年0.35个百分点的速率不断减少，共减少5.592 7万hm2，耕地压力不断增加。草地和未利用地也呈现减少态势，减少速率较为缓和。林地和水域先减少后增加， 2014年相较于1992年所占比例有所提升，体现出近年来哈尔滨对生态环境保护工作的重视及取得的成果。建设用地不断增加， 22年期间增加5.64个百分点，增加面积为3.996 8万hm2，反映出在城镇化水平不断提高、人口压力增大、经济不断发展的背景下，对建设用地扩展的需求。

表1 1992年、2003年和2014年研究区各类土地利用面积及所占比例

土地利用类型199220032014面积(万hm2)百分比(%)面积(万hm2)百分比(%)面积(万hm2)百分比(%)耕地51.248072.3549.900670.4545.655364.45林地11.377616.0611.178915.7811.415016.12草地0.48340.680.44290.630.33650.48水域2.14213.022.01392.844.01495.67建设用地5.36977.587.179310.149.363513.22未利用地0.21280.300.11790.170.04850.07

3.2 土地利用变化的空间转换分析

单纯的土地利用类型面积的增减不能较好地反映各类型间的转换情况，该文通过IDRISI中的空间分析功能，得到1992～2003年和2003～2014年土地利用面积转移矩阵(表2)，从而充分了解土地利用空间变化过程。

表2 1992～2014年研究区土地利用面积转移矩阵 hm 2

1992～2003年间，耕地净转出1.350 4万hm2，主要转出为建设用地2.354 7万hm2和林地1.294 0万hm2； 林地净转出1 436hm2，转入转出基本平衡； 建设用地净转入1.810 7万hm2，耕地转入的2.354 7万hm2是建设用地增加的主要原因； 水域少量减少，净转出1 811hm2，水域和耕地之间的相互转换明显； 草地和未利用地由于面积较小，总体转入转出不明显。但空间位置变化明显。2003～2014年间，耕地持续减少，净转出4.195 2万hm2，主要转出为建设用地2.932 2万hm2，水域2.142 7万hm2和林地1.298 3万hm2，主要转入未建设用地8 797hm2和林地7 399hm2； 林地由净转出变为净转入4 297hm2，主要由于转出为耕地的面积减少和耕地的转入增多； 建设用地不断增加，增加速度高于上时段，净增加2.184 6万hm2，耕地仍然是建设用地重要来源，转入2.932 2万hm2，转出主要为耕地8 797hm2； 水域较上时段明显不同，不降反增加1.756 5万hm2，转入主要为耕地2.142 7万hm2，总转出由上时段的1.197 8万hm2减少到4 942hm2，可见环境保护效果明显。

图1 1992、2003和2014年研究区土地利用情况

图2 土地利用重心点空间位置

3.3 土地利用重心变化

从上述分析可知， 1992～2014年各土地利用类型中，耕地和建设用地变化最为显著，同时也是区域经济发展和粮食安全的重要影响因子，故该文以耕地和建设用地两种地类为代表，分析其重心变化特征，借助GIS空间分析技术和式(1)得到图2，并将2025年预测结果展布到图2，用于后期预测分析。1992～2003年，耕地重心向东南方向移动356m，移动速度较慢； 建设用地向西移动2 307m，表明建设用地向四周扩张的同时，更偏向于哈尔滨建成区西南方向发展。2003～2014年，耕地重心向北移动886m，较上时间段速度明显加快，未来几年有加快移动的趋势； 建设用地向南移动1 791m，移动速度减缓。整体来看， 1992～2014年间，受研究区北部粮食主产区的影响，耕地重心整体向北方向移动，向呼兰区、巴彦县方向靠拢，建设用地由于松花江流域的阻隔和气候的限制，整体向有利于经济发展的西南方向移动。

3.4 CA-Markov模型预测精度检验

模型的一个重要性步骤是检验。一般情况下，利用某一土地覆盖状况已知时段的预测，估算对一个过程的认识和模仿的功效，这就是用于验证的测试。Kappa统计是较为常用的精度检验模型(式3)，其判别标准为Kappa值从0～1表示两幅图像从完全不吻合到完全吻合的转变，其中Kappa≥0.75表明模型预测精度较高，Kappa≤0.4预测精度较差。该文应用CA-Markov模型预测的2014年土地利用类型预测结果(图3b)与2014年已知土地利用数据(图3a)进行Kappa验证，结果为Kappa=0.878 0，检验精度较高，能够进行后续预测研究，公式表达如下[18]。

(3)

式(3)中，p0表示土地利用分布预测正确比例；pc表示随机情形下预测正确比例；pp表示理想状态下预测完全正确的比例。

图3 2014年研究区土地利用现状图(a)和预测图(b)

图4 2014年研究区土地利用现状图和2025年研究区土地利用预测

表3 2014～2025年研究区土地利用转移概率矩阵

土地类型耕地林地草地水域建设用地未利用地耕地0.7410.0520.0020.0870.1180.000林地0.1710.7760.0000.0200.0330.000草地0.3390.0240.5360.0070.0950.000水域0.3280.0150.0000.6410.0150.000建设用地0.2410.0170.0020.0030.7360.001未利用地0.3820.0810.0030.0760.2100.249 注:横向代表变化基准年2014年,纵向代表变化目标年2025年

表4 2003～2025年研究区土地利用变化情况

土地利用类型20252003～20142014～2025预测面积(万hm2)面积变化(万hm2)动态度(%)面积变化(万hm2)动态度(%)耕地39.4527-4.2453-8.51-6.2026-13.59林地11.63790.23602.110.22301.95草地0.2895-0.1065-24.03-0.0470-13.96水域6.75432.001099.362.739568.23建设用地12.77642.184230.423.412936.45未利用地0.0398-0.0694-58.87-0.0087-17.96

3.5 CA-Markov模型预测结果分析

在模型预测精度检验检验合格的基础上，以2014年土地利用数据为基础，制作用于2025年预测的适宜性图集，并与Markov模型预测的2014～2025年土地面积转移矩阵(表3)同时作为转换规则进行2025年研究区土地利用变化的预测，结果见图3、图4和表4。

预测结果表明， 2025年研究区耕地面积39.452 7万hm2，延续1992年以来减少的趋势，较2014年共减少6.202 6万hm2，且减少速率与2003～2014年时段相比不降反增，耕地保护压力加大，耕地重心加速向北移动，向粮食主产区靠近， 2014～2025年间预计移动1 780m； 林地面积11.637 9万hm2，较2014年稳中有增，变化不显著，空间分布多集中于东南部； 草地面积减少到2 895hm2，与2014年相比减少41hm2，相比于上一时间段减少速度降低，但动态度反映草地变化仍然活跃； 水域面积6.254 3万hm2，增加2.739 5万hm2，空间分布变化体现为在2014年原有水域的基础上扩宽、延长，水域的大量增加将使生态环境得以改善； 建设用地面积12.776 4万hm2，自1992年以来不断增加， 2014～2025年间增加3.412 9万hm2，增速高于2003～2014时段。土地利用变化活跃，空间分布结果表明，建设用地增加过程集聚效应显著，增加的面积多分布于原建设用地周边，且占用耕地现象明显，建设用地重心继续向西南方向移动， 2014～2025年预计移动1 114m，移动速度减缓； 未利用地面积398hm2，比2014年减少87hm2。

4 结论与讨论

(1)该文基于Markov模型和GIS、RS技术，分析哈尔滨市建成区1992～2014年间土地利用类型演变特征。选取经济生态方面的6个影响因子，水域和坡度两个约束条件和GDP、与水域距离、与建设用地距离、与道路距离4个适宜性因子参与到转换规则的制定中，假设预测时间段内未发生突发性事件，应用CA-Markov模型预测了2025年土地利用情况。试验得到的2014年土地利用结果和实际解译数据Kappa系数达0.878 0，表明转换规则的制定可行，模型可信度较高。转换规则的制定提高了预测结果的准确合理性，使预测结果能够为土地利用规划的修编和生态环境的保护提供决策支持。

(2)土地利用变化演变特征主要集中于： 1992～2014年间，耕地面积持续减少，耕地重心向北加速移动，向粮食主产区呼兰、巴彦县方向靠近； 建设用地由于耕地的大量转入使其在原有空间位置上不断扩张，建设用地重心向西南方向移动； 林地在数量和空间分布上相对稳定； 草地和未利用地面积变化微弱，但空间分布变化明显； 水域面积先减后增，反映出环境保护政策中水域保护效果明显，耕地和水域之间大量转换。

(3)2025年预测表明，耕地继续减少，且减少速率不见缓和，耕地保护压力较大，耕地重心持续向北移动，且速度加快； 林地延续前期变化趋势，稳中有升； 水域面积明显增加，土地利用变化活跃，对未来环境改善有积极促进作用； 建设用地不断增加，占用大量耕地，扩展明显，建设用地重心持续向西南移动，移动速度减慢； 草地和未利用地总面积较小，呈零星分布。

研究中，对模型参数选择采用统一的30m×30m元胞大小， 5×5滤波器，未考虑元胞大小以及滤波器对预测结果的影响； 土地利用变化受到社会、经济、生态等多方面影响，该文由于数据的不足和部分因子定量化困难的原因，只考虑了经济生态方面的6个影响因子，对影响土地利用的社会因子如永久基本农田、人口密度等分析不足。在未来的研究中，加大对上述问题的探索力度，实现更好的空间预测效果。

[1] 后立胜，蔡运龙.土地利用/覆被变化研究的实质分析与进展评述.地理科学进展， 2004， 23(6)： 96～104

[2] 王兵， 臧玲.我国土地利用/土地覆被变化研究近期进展.地域研究与开发， 2006, 25(2): 86～91

[3] Weber C，Puissant A.Urbanization pressure and modeling of urban growth：Example of the tunis metropolitan area.Remote Sensing of Environment， 2003， 86(1/2)： 341～352

[4] 李志， 刘文兆，郑粉莉.基于CA-Markov模型的黄土塬区黑河流域土地利用变化.农业工程学报， 2010， 26(1)： 346～352， 391

[5] Bormann H，Breuer L，Graff T，et al.Analysing the effects of soil properties changes associated with land use changes on the simulated water balance：a comparison of three hydrological catchment models foe scenario analysis.Ecol Model， 2007， 209(1)； 29～40

[6] 刘纪远， 邵全琴，延晓冬，等.土地利用变化对全球气候影响的研究进展与方法初探.地球科学进展， 2011, 26(10): 1015～1022

[7] 赵丽， 张蓬涛，周智.环首都贫困区土地利用变化对生态系统的影响及驱动因素分析——以河北涞水县为例.中国农业资源与区划， 2013, 34(5): 74～81

[8] Brath A，Montanan A，Moretti G.Assessing the effect on flood frequency of land ues change via hydrological simulation(with uncertainty).Journal of Hydrology， 2006， 324(1/2/3/4)： 141～153

[9] 汪佳莉， 吴国平，范庆亚，等.基于CA-Markov模型的山东省临沂市土地利用格局变化研究及预测.水土保持研究， 2015, 22(1): 212～216

[10]赵冬玲， 杜萌，杨建宇，等.基于CA-Markov模型的土地利用演化模拟预测研究.农业机械学报， 2016, 47(3): 278～285

[11]王友生， 余新晓，贺康宁，等.基于CA-Markov模型的藉河流域土地利用变化动态模拟.农业工程报， 2011， 27(12)： 330～336， 442

[12]肖明， 吴季秋，陈秋波，等.基于CA-Markov模型的昌化江流域土地利用动态变化.农业工程学报， 2012, 28(10): 231～238

[13]何丹， 周璟，高伟，等.基于CA-Markov模型的滇池流域土地利用变化动态模拟研究.北京大学学报(自然科学版)， 2014, 50(6): 1095～1105

[14]龚文峰， 袁力，范文义.基于CA-Markov的哈尔滨市土地利用变化及预测.农业工程学报， 2012, 28(14): 216～222

[15]许月卿， 崔丽，孟繁盈.大城市边缘区土地利用变化与社会经济发展关系分析——以北京市平谷区为例.中国农业资源与区划， 2008, 29(4): 16～21

[16]杜萌， 赵冬玲，杨建宇，等.基于元胞自动机复合模型的土地利用演化模拟——以北京市海淀区为例.测绘学报， 2015， 44(S1)： 68～74

[17]Carver，S.J.Integrating multi-criteria evaluation with geographical information systems，International Journal of Geographical Information Systems， 1991, 5(3): 321～339

[18]刘洁， 李宏，马勇刚.基于CA-Markov模型的中亚典型城市土地利用变化预测分析.水土保持研究， 2014, 21(3): 51～56