基于CA-Markov模型的挠力河流域土地利用动态模拟

周　浩，雷国平①，赵宇辉，路　昌，张　博

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



基于CA-Markov模型的挠力河流域土地利用动态模拟

周浩1，雷国平1①，赵宇辉2，路昌1，张博2

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

摘要：以三江平原腹地挠力河流域为研究对象，基于1990和2013年2期遥感数据，结合马尔科夫(Markov)转移矩阵、土地利用动态度模型和土地利用程度模型分析流域1990—2013年土地利用变化情况，运用CA-Markov模型对2025年土地利用格局进行模拟预测。结果表明：流域土地利用结构变化以耕地的内部和外部转换为主，水田化现象是最主要的景观变化特征，未利用地面积由于大量垦荒而大幅度减少；土地利用程度指数由1990年的251.46上升至2013年的270.48，土地开发利用程度不断加大；模型模拟显示流域Kappa系数为0.82，预测结果可信；2013—2025年流域土地利用类型变化趋势和速率与1990—2013年存在差异，林地变化幅度最大，水田面积仍然保持增长，但增长速度有所下降，单一动态度降为2%，旱地面积基本保持不变；未来12 a土地利用程度指数持续上升，流域受人类社会因素扰动作用的强度将会进一步加大。

关键词：CA-Markov模型；挠力河流域；土地利用/覆盖变化(LUCC)；空间分析模型

土地利用/覆盖变化(LUCC)作为区域土地资源利用方式的直接反映[1],其相互作用、相互依存的各要素通过能量、物质和信息的交流而形成的时空系统[2],是全球变化和可持续发展的核心和难点[3]。LUCC是相当复杂的过程,在不同时空尺度上进行模拟研究有助于揭示人类社会影响下区域及生态环境的变化趋势及机理等[4-5],为制定科学、合理的土地利用管理及区域发展策略提供决策支持。包括挠力河流域在内的三江平原是全球变化最为敏感的区域之一,在全球变化研究中意义重大[6]。

目前用来分析和模拟LUCC的模型包括元胞自动机(CA)[7]、小尺度土地利用变化及效应模型(CLUE-S)[8]、马尔科夫(Markov)预测模型和基于智能体[3]等。其中CA模型主要从时空耦合的角度出发,着眼于元胞的局部相互作用,存在局限性;而传统的Markov模型仅从数量上进行长期预测,难以预测土地利用空间格局变化[3]。近年来被学者广泛探讨和应用的CA-Markov模型综合了CA模型模拟复杂系统行为的能力和Markov模型的长期预测优势[9-12],既能提高土地利用类型转化的预测精度,又可以有效模拟土地利用空间格局变化[10],有利于精确挖掘区域土地利用时空演变信息。位于三江平原腹地的挠力河流域作为该地区土地开发历史最早的流域,是我国最重要的垦区之一[13-14],近年来当地农业结构调整促使大量易涝旱地转变为水田,土地利用变化剧烈,开展该地区的土地利用变化及模拟研究具有重要的现实意义。

笔者以挠力河流域为研究对象,在ArcGIS和IDRISI平台支持下,基于遥感影像数据分析流域土地利用变化空间格局,运用CA-Markov模型模拟预测未来土地利用变化情况,以期为土地资源的合理利用、规划与保护提供决策依据。

1研究区概况

挠力河流域位于黑龙江省三江平原境内,地理位置为45°43′～47°45′ N,131°31′～134°10′ E。东南以完达山为界,东与乌苏里江相接,流域面积约为2.49×104km2,其地理位置如图1所示。

图1　挠力河流域地理位置示意

该流域属中温带大陆性季风气候区,夏季高温多雨,冬季寒冷漫长,最热月7月平均气温为21.9 ℃,最冷月1月平均气温为-18.6 ℃。挠力河流域多年平均降水量为518 mm,降水分布不均匀,6—9月降水占全年降水量的72%,春季干旱频繁,秋季多洪涝灾害。流域范围内地形呈现西南高、东北低的态势,水系自西南流向东北。地貌类型主要由山地与平原组成,山地占流域面积的38.3%,主要分布于流域西南部和南部;平原占61.7%,主要分布于流域北部和中部的内、外七星河及挠力河中游地区。挠力河流域农业开发活动非常活跃,尤其在1980年国家进入经济迅速发展时期后。目前挠力河流域已建成6个县,有7个现代化农场,总人口达125万人,其中农业人口占65.4%,该区已成为三江平原主要产粮区和国家重要的商品粮基地。

2研究方法

2.1数据来源与处理

土地利用数据原始信息源来自美国陆地资源卫星1990和2013年Landsat TM/OLI多光谱遥感影像,精度分别为30和15 m,经几何纠正及RGB假彩色合成。参照全国土地利用分类体系并结合流域土地利用现状和研究目标,确定研究区土地利用类型为耕地(包括旱地和水田)、林地、草地、水域、建设用地和未利用地6大类。根据不同土地利用类型的影像色调和纹理等特征建立解译标志,参照流域地形图和各种历史资料配合进行人工目视解译,对所得数据进行查错、修改和拼接,并结合Google Earth软件进行精度验证,最终得到1990和2013年挠力河流域土地利用现状图。将矢量数据进行数据格式的转换、数据的提取和统计。

2.2CA-Markov模型

2.2.1Markov模型

土地利用变化没有固定而精确的变化模式。为便于研究,常采用Markov箱式模糊模型进行拟合预测,它是基于Markov过程而形成的系统预测随机事件发生概率和优化控制的理论方法,常用于具有无后续性特征的地理事件预测,Markov过程下的土地利用类型变化只与前一时刻土地利用类型相关,它们之间相互转换的面积数量或者比例即为状态转移概率,依据贝叶斯条件概论公式,Markov过程下的土地利用变化如下所示[11]:

S(t+1)=Pij×S(t)。

(1)

式(1)中,S(t+1)和S(t)分别为t+1和t时刻的系统状态;Pij为状态转移概率矩阵。

2.2.2CA模型

CA是具有时空计算特征的动力学模型,其特点表现为时间、空间和状态是离散的,模型变量都只有有限个状态,且状态的改变均表现出局部特征,特别适用于复杂的地理空间系统动态模拟研究[6],该模型如下所示:

S(t,t+1)=f[S(t),N]。

(2)

式(2)中,S为元胞有限、离散的状态集合;f为转换规则;N为元胞的邻域。

2.2.3CA-Markov预测模型

CA-Markov模型综合了Markov模型长期序列预测和CA模型空间变化模拟的优势,可以较好地在数量和空间上对土地利用变化时空格局进行模拟[15]。在土地利用栅格图中,每个像元相当于1个元胞,其土地利用类型为元胞状态。模型在IDRISI软件的支持下利用面积转移矩阵和条件概率图像实现模型的运算,模拟未来土地利用变化格局,具体操作如下:

(1)基于Markov模型对1990和2013年2期土地利用现状图进行叠加,得到1990—2013年研究区土地利用类型转移概率矩阵和转移面积矩阵,其中转移概率矩阵作为转换规则参与模拟运算。

(2)建立转换适宜性图像集(MCE)。为了能够较好地保持1990—2013年的转移趋势,使用Markov模型输出的条件概率图像作为转换适宜性图像集。

(3)确定CA滤波器及循环次数。CA滤波器用于创建相邻元胞单元具有显著空间意义的权重因子,使其作用于元胞而确定其状态的改变,确定CA标准5×5的滤波器,即认为一个元胞周围150m×150m的矩形空间对该元胞状态影响显著。由于预测基期土地利用图为2013年的,模拟目标为2025年,因此将CA循环次数设为12。

2.3指数模型

2.3.1土地利用动态度模型

在自然和人为因素影响下,区域时间段、空间位置存在差异,不同土地利用类型的变化幅度和速度是不同的[5],变化幅度为研究时段初、末的面积差值,变化速度可以用动态度模型进行定量描述,计算公式如下:

(3)

式(3)中,K为土地利用动态度;Ua和Ub分别为研究时段初、末的面积,km2;T为研究时段长度,a。

2.3.2土地利用综合程度模型

土地系统作为自然和社会因素综合作用下的复杂综合体,其利用程度可以体现土地利用的广度与深度。为进一步充分反映LUCC效应,参照刘纪远等[16]提出的土地利用程度分级标准来定量化分析土地利用程度的高低,土地利用程度综合指数及变化模型如下所示:

(4)

(5)

(6)

式(4)～(6)中,L、ΔLb-a和R分别为土地利用程度综合指数、指数变化量和指数变化率;Ai为第i级土地利用程度分级指数;Ci为第i级土地利用程度分级面积比例;La和Lb分别为研究时段初、末时间土地利用程度综合指数;Ca，i和Cb，i分别为研究时段初、末第i级土地利用程度面积比例。

3结果与分析

3.1土地利用动态演变

3.1.1土地利用结构变化及特征

基于Markov转移矩阵,统计挠力河流域1990—2013年土地利用变化情况(表1)。可以看出,1990—2013年流域土地利用变化幅度较大。42个土地转换类型中共有28个子类型发生转换,其中旱地→水田转换类型面积达3 466km2,占土地转换总面积的56.93%,未利用地→旱地、林地→旱地和水田→旱地转换类型面积依次为1 070、857和695km2,以上4种转换子类型面积之和占土地转换总面积的86.48%,1990—2013年间流域土地利用变化以耕地转换为主(包括内部转换和外部转换)。

具体而言,20多a间水田面积增加3 400km2,其单一动态度最大,达7.52%,水田化现象是该流域最主要的景观变化特征。其次为未利用地,流域内大量垦荒导致其面积持续下降,1990—2013年间下降幅度达1 586km2,动态度为-3.68%。20世纪90年代初期以来,当地农业结构调整促使大量易涝旱地转变为水田,尽管期间大量未利用地转换为旱地,但旱地仍然减少944km2。与1990年相比,林地共减少881km2,其中有857km2转变为旱地。除水田外,建设用地面积也呈增加趋势,20多a间建设用地共增加68km2,通过转移矩阵分析其用地来源可知,旱地、林地、水田和未利用地均有部分转变为建设用地。流域内草地面积相对较小,1990年面积仅为6km2,水域用地多为湖泊和河流水面,转换程度相对有限,20多a间两者面积分别下降51和2km2,其中水域用地单一土地利用动态度为所有用地类型中最低,为-0.04%。

表1挠力河流域1990—2013年土地利用类型动态变化

Table 1Dynamics of land use types in the Naoli River Valley from 1990 to 2013

km2

表中某一行数据为1990年各土地利用类型转变为2013年该行对应的土地利用类型的面积。1)单位为%。

3.1.2土地利用变化综合程度分析

1990和2013年土地利用现状图见图2。为直观反映挠力河流域内部土地利用综合程度差异,消除空间尺度效应,将2期土地利用矢量图进行1 km×1 km的矢量栅格切割,计算单个栅格用地类型面积比例,结合式(6)得到1990和2013年2期土地利用综合程度指数图,并以此为基础计算得到1990—2013年土地利用综合程度指数变化量和变化率(图3)。

20多a间挠力河流域土地利用综合程度指数整体呈上升趋势,土地利用受人类扰动作用大。土地利用综合程度指数由1990年的251.46上升至2013年的270.48,土地利用综合程度变化指数为19.01,该时段流域土地利用整体处于发展时期。从土地利用程度指数图(图3)可以看出,变化主要集中于流域中部的内、外七星河流域和北部平原地带,1990年该地区分布有大量综合指数阈为100～150的土地类型,2013年该阈值范围的土地类型仅见零星分布,这是因为大量未利用地和滩涂用地被开垦为农田,土地综合利用指数不断提高。1990—2013年该地区大片土地利用综合程度指数变化量处于100～200之间。值得注意的是,流域内零星分布有变化量为负值的土地类型,结合土地利用现状图可以发现1990年负值地区多为建设用地,土地利用程度指数相对较高,经合村并点土地整理后多转变为指数较低的旱地。流域平均土地利用程度指数变化率为7.56%,年均变化量为0.33%,土地利用程度指数变化率最高值达200%,多集中于垦荒较多的内、外七星河流域地区。

3.2CA-Markov预测结果分析

3.2.1预测精度检验

Kappa系数(K)常用来评价遥感分类精度和图件间相似程度,能够从数量和空间角度上定量反映土地利用变化模拟过程中丢失的空间信息量。为了检验CA-Markov模型预测不同土地利用类型的合理性,运用IDRISI软件的Crosstab模块对基于1990—2013年Markov转移概率矩阵生成的2013年土地利用预测值进行Kappa系数检验,Kappa系数计算公式为

K=(Po-Pc)/(Pp-Pc)。

(7)

式(7)中,Po为正确模拟的比例;Pc为随机情况下期望的正确模拟比例;Pp为理想分类情况下的正确模拟比例,其值为100%。当K>0.75时,图件一致性较高,变化相对较小;K处于0.4～0.75之间时,图件一致性一般,变化较为明显;K<0.4时,模拟效果差[17]。

将2013年土地利用现状图与模拟图进行求差栅格运算,模拟正确区域的栅格属性值为0,通过对属性值为0的栅格数进行提取,可得到模拟精度。由检验结果可知,共有20 241 082个栅格得到正确模拟,占栅格总数的76.93%,模拟结果精度较高。2013年模拟值与真实值的总体Kappa系数达0.82,大于0.75,反映两者一致性较高,变化相对较小,模拟结果可信。不同土地利用类型的Kappa系数存在差异,建设用地、林地、旱地和水域用地的模拟精度较高,而未利用地和草地的Kappa系数分别仅为0.27和0.10(表2)。

图2　挠力河流域1990、2013年土地利用现状和2025年土地利用格局预测

图3　挠力河流域土地利用综合程度评价

表22013年挠力河流域不同土地利用类型实际面积与模拟面积

Table 2Measured and simulated areas of land use by type in 2013 in the Naoli River Valley

土地利用类型实际面积/km2模拟面积/km2Kappa系数草地 9 580.10旱地949098640.78建设用地4703680.91林地616053760.90水田536533470.59水域用地22813980.76未利用地29015980.27

3.2.2模拟结果分析

根据预测结果(表3、图2),2025年流域草地、旱地、建设用地、林地、水田、水域用地和未利用地面积分别为14、9 463、425、4 580、6 649、769和109 km2。与2013年相比,除林地、水田和水域用地外,其他土地利用类型面积变化幅度较小。林地面积变化量最大,由2013年的6 160 km2减少到2025年的4 580 km2。模拟结果表明,今后随着挠力河流域未利用地和滩涂用地的开发殆尽,耕地后备资源逐渐减少,未来林地更可能受人类扰动影响,逐渐转变为耕地等土地利用程度指数较高的用地类型。

表32013—2025年挠力河流域土地利用变化

Table 3Predicted changes in land use during the period from 2013 to 2025 in the Naoli River Valley

土地利用 类型 2025年预测值2013—2025年变幅面积/km2比例/%面积/km2比例/%动态度/%草地 140.06 50.024.71旱地946343.00-27-0.12-0.02建设用地4251.93-45-0.20-0.80林地458020.81-1580-7.18-2.14水田664930.2112845.832.00水域用地7693.495412.4619.76未利用地1090.50-181-0.82-5.19

2025年水田面积将增加1 284 km2,相对于1990—2013年的水田变化动态度(7.52%),2013—2025年其动态度继续保持为正值(2.00%),可见在转移趋势不变的情景下,未来水田面积仍然保持增长,但增长速度有所下降,该变化趋势与实际情况较为符合,在三江平原旱改水农业结构调整的大背景下,未来仍会有部分易涝旱地转换为水田,但水田面积过度扩张导致的农业灌溉用水供应不足、土壤盐渍化严重等负面因素会抑制水田的扩张。

2025年旱地面积较2013年基本保持不变,面积仅下降27 km2,而相对面积比例也仅变化0.12%。与1990—2013年建设用地面积缓慢增加趋势不同的是,未来农村居民点整理、合村并点等政策的实施将会导致建设用地面积逐渐减少。模拟结果显示,与2013年相比,2025年建设用地面积将会下降至425 km2,变化动态度为-0.80%。草地和未利用地受自然因素影响较大,在保持转移趋势不变的前提下,CA-Markov模型难以考虑气候变化等自然因素的影响,因此模型模拟效果较差,可信度不高。

通过将2013年土地利用类型分布图和2025年模拟图叠加发现,流域以河流水面、滩涂为主的水域用地未来地理位置基本保持不变,但面积出现较大变化,较2013年增加541 km2。挠力河流域土地利用程度指数将由2013年的270.48上升至2025年的276.57,预示未来流域土地受人类社会因素扰动作用的强度将进一步加大。

4讨论

自建国以来挠力河流域农业开发活动非常活跃,经历了大规模的土地开发,尤其是1980年国家进入经济迅速发展时期后,农业化的发展、人口的增加和市场规律作用以及经济利益的驱使使得该地区土地利用结构变化剧烈。研究结果表明,20世纪90年代初期以来当地“以稻治涝”政策的推行，导致水田化现象成为该流域最主要的景观变化特征,土地利用结构变化以耕地内、外部转换为主。

深入开展土地利用/覆盖过程及模拟研究,有助于为土地利用决策提供更加科学的依据。但由于遥感影像解译结果与真实情况存在一定误差,且CA-Markov模型侧重于本身机制的模拟,其转换规则的设定包括各用地类型的内在联系,难以充分考虑地理要素、经济社会因素等影响。元胞大小、计算时间等内在参数存在不确定性,模型模拟难以规则化。作为国家重要的商品粮基地,挠力河流域土地利用结构变化较其他地区更易受国家政策的影响,模拟预测中难以考虑该因素,如国家近年来围绕该地区采取的保护湿地、有序化开发耕地后备资源等一系列措施。在保持原有转移趋势不变的模型背景下,未利用地的模拟结果可能与实际情况存在较大差异,由此导致其Kappa系数偏低。后续研究可对以上几个方面加以改进,并提高模型的适用性。

5结论

挠力河流域作为我国重要的商品粮基地之一,近几十年来人类的开发活动对流域土地利用结构变化产生强烈影响,该研究以20世纪90年代初当地农业结构调整作为时间切入点,研究流域近20多a间土地利用变化情况,运用CA-Markov模型模拟预测并分析未来土地利用格局变化情况,主要结论如下:

(1)1990—2013年挠力河流域土地利用结构变化明显,其变化以耕地内部转换和外部转换为主;20多a间水田面积增加3 400 km2,而大量垦荒促使未利用地面积下降1 586 km2。林地面积下降881 km2,林地绝大部分转变为旱地。建设用地共增加68 km2,供地来源广;流域各用地类型间面积变化差异大,单一土地利用变化动态度介于-0.04%～7.52%之间。

(2)20多a间流域土地利用受人类扰动作用逐渐加大,综合指数整体呈现明显上升趋势,土地利用程度指数由1990年的251.46上升至2013年的270.48,变化主要集中于流域中部的内、外七星河流域和北部平原地带。

(3)CA-Markov模型模拟的Kappa系数为0.82,整体模拟结果可信度高。2025年的模拟结果表明,与2013年相比,林地、水田和水域用地变化幅度大,其中林地变化幅度达最大,减少1 580 km2,未来林地倾向于开发为耕地等土地利用程度指数较高的用地类型;水田面积仍然保持增长,但增长速度有所下降,旱地面积基本保持不变;未来流域受人类社会因素扰动作用的强度进一步加大。应针对土地利用类型变化特征的差异,制定差别化的土地利用开发与管理政策。

参考文献：

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

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

[3]赵冠伟,陈颖彪,陈健飞,等.CA-Markov模型的空间尺度敏感性研究[J].地理科学,2011,31(8):897-902.

[4]侯西勇,常斌,于信芳.基于CA-Markov的河西走廊土地利用变化研究[J].农业工程学报,2004,20(5):286-291.

[5]宋开山,刘殿伟,王宗明,等.三江平原过去50年耕地动态变化及其驱动力分析[J].水土保持学报,2008,22(4):75-81.

[6]王宗明,宋开山,刘殿伟,等.1954—2005年三江平原沼泽湿地农田化过程研究[J].湿地科学,2009,7(3):208-217.

[7]周成虎,孙战利,谢一春.地理元胞自动机研究[M].北京:科学出版社,2001：79-82.

[8]唐华俊,吴文斌,杨鹏,等.土地利用/土地覆被变化(LUCC)模型研究进展[J].地理学报,2009,64(4):456-468.

[9]何丹,金凤君,周璟.基于Logistic-CA-Markov的土地利用景观格局变化:以京津冀都市圈为例[J].地理科学,2011,31(8):903-910.

[10]郑青华,罗格平,朱磊,等.基于CA-Markov模型的伊犁河三角洲景观格局预测[J].应用生态学报,2010,21(4):873-882.

[11]杨国清,刘耀林,吴志峰.基于CA-Markov模型的土地利用格局变化研究[J].武汉大学学报(信息科学版),2007,32(5):424-418.

[12]姜广辉,张凤荣,孔祥斌.北京山区农村居民点整理用地转换方向模拟[J].农业工程学报,2009,25(2):214-221.

[13]栾兆擎,胡金明,邓伟,等.人类活动对挠力河流域径流情势的影响[J].资源科学,2007,29(2):46-51.

[14]孙贤斌,刘红玉,李玉凤,等.基于CA-Markov模型土地利用对景观格局影响辨识[J].生态与农村环境学报,2009,25(1):1-7,31.

[15]程刚,张祖陆,吕建树.基于CA-Markov模型的三川流域景观格局分析及动态预测[J].生态学杂志,2013,32(4):999-1005.

[16]刘纪远,布和敖斯尔.中国土地利用变化现代过程时空特征的研究:基于卫星遥感数据[J].第四纪研究,2002,20(3):229-239.

[17]布仁仓,常禹,胡远满,等.基于Kappa系数的景观变化测度:以辽宁省中部城市群为例[J].生态学报,2005,25(4)：778-785.

(责任编辑： 许素)

Simulation of Dynamics of Land Use in Naoli River Valley Based on CA-Markov Model.

ZHOUHao1,LEIGuo-ping1,ZHAOYu-hui2,LUChang1，ZHANGBo2

(1.Land Management Institute, Northeastern University, Shenyang 110004, China;2.College of Resource and Environment, Northeast Agricultural University, Harbin 150030, China)

Abstract:Based on the Landsat TM/ETM remote sensing image data of 1990 and 2013 of the Naoli River Valley in the heart of the Sanjiang Plain, the land use dynamics model and land use intensity model were used in combination with the Markov transfer matrix to analyze changes in land use during the period of 1990-2013 and the CA-Markov model was to simulate and predict land use pattern in the Naoli River Valley in 2025. Results show that the changes in land use structure was mainly dominated by conversion of farmland from paddy to upland or vice versa and alienation of farmlands and reclamation of wastelands. The conversion of upland to paddy was the main feature of the changes in landscape. Reclamation of unused land led to substantial reduction of wasteland in area. The land use intensity index of the valley rose from 251.46 in 1990 to 270.48 in 2013, which indicates that land use increased steadily in exploitation degree in the Naoli River Valley. The simulation using the models demonstrates that the Kappa coefficient of the valley was 0.82, indicating that the prediction was reliable. The variation of land use types in the valley in 2013-2025 may differ from that in 1990-2013 in trend and rate. The change in woodlands will be the widest in range and the area of paddy fields will keep on rising, but at a declining rate with the single dynamic degree falling down to 2%; The area of upland will stay almost unchanged. The land use intensity index will continue rising and the valley will continue to be disturbed by human activities with rising intensity during the period from 2013 to 2025.

Key words:CA-Markov model;Naoli River Valley;LUCC;spatial analysis model

作者简介：周浩(1990—),男,安徽安庆人,博士生,主要从事土地利用与规划方面的研究。E-mail: zhouhao7404@163.com

DOI：10.11934/j.issn.1673-4831.2016.02.013

中图分类号：X87；F301.24

文献标志码：A

文章编号：1673-4831(2016)02-0252-07

通信作者①E-mail: guopinglei@126.com

基金项目：教育部博士学科点基金博导类项目(20112325110007);黑龙江省国土资源科研项目(黑国土科研201411)

收稿日期：2015-04-24