基于GIS-CA-Markov的白马河流域土地利用变化模拟分析

荆延德1,李婷婷

(1.山东省高校南四湖湿地生态与环境保护重点实验室,山东 济宁 273165;2.曲阜师范大学 地理与旅游学院,山东 日照 276826)

1 引言

土地利用/植被变化(LUCC)对多种生态环境因素产生直接或间接的影响,进一步影响到区域和全球的可持续发展[1,2]。随着社会经济的快速发展,土地利用的变化速度和规模发生了巨大转变,探讨土地利用变化的时空变化规律,预测其未来变化速度与方向成为研究的热点问题[3]。

目前进行土地利用变化研究的模型主要有遗传神经网络模型(GNN)[4]、CLUE模型[5]、CLUES-S模型[6]、系统动力学模型(SD)[7]、马尔科夫链(Markov Process)[8]、元胞自动机模型(CA)[9,10]、CA-Markov模型[11,12]等。其中,元胞自动机—马尔科夫(CA-Markov)模型利用Markov模型对土地利用变化进行预测并得出状态转移图集,并利用元胞自动机模型对土地利用变化进行空间上的处理,这样既结合了马尔科夫链长期预测模拟土地类型数量结构的优势,又具有元胞自动机模型对土地利用景观格局的时空动态模拟功能,能很好地预测土地利用变化的发展趋势[13]。目前学者们主要研究模型转换规则的获得[14,15],对时空尺度敏感性的研究很少考虑,这对土地利用模拟预测结果的准确性造成了很大的影响。研究表明,CA模型在预测过程中具有时空尺度敏感性,而只有在一定的时间尺度和空间尺度范围内进行预测模拟才能得到与实际情况较符合的结果[16-20]。本文以已有的预测模拟研究中常用的像元大小和时间步长为候选条件,对其进行模拟精度的对比,选择较合适的条件进行模拟预测,以增加结果的准确性。

白马河发源于山东省邹城市北部,是“南水北调”东线工程的重要水资源承载地之一,也是作为南四湖流域最重要的子流域之一,因此该区域的土地利用变化对南四湖的水资源调控和生态管理起着举足轻重的作用。作为山东省区域经济发展中的重要部分,该区域景观格局变化对流经地经济发展具有很大的影响[21,22]。本文以济宁市白马河流域为研究区域,在CA-Markov模型运行过程中通过改变研究区域的元胞空间尺度和时间尺度进行模拟预测,对比模拟精度,试图探究元胞尺度、时间间隔等参数对CA-Markov模型模拟结果的影响机制,为白马河流域的土地利用规划和相关政策的制定提供理论依据。

2 研究区域与数据

2.1 研究区域

白马河流域位于山东省西南部,是华北地区较大的淡水湖。白马河全长60km,流经曲阜、兖州、邹城、微山等县(市)地区,于微山县鲁桥镇九孔桥村入独山湖。为了更好地研究土地利用变化机理,本文以白马河及其所流经县市区为研究区域进行研究，见图1。

图1 白马河流域的地理位置

本研究区域地理位置位于116°41′—117°11′E、35°07′—35°33′N之间,气候属于暖温带大陆性季风气候,四季分明、雨热同期,降水多集中在夏秋两季,年平均降水量约691.19mm,年平均气温约14.7℃;地势较低平,大致由东向西倾斜,大支流有大沙河、望云河、石强河、十里沟和七里沟等。白马河流域作为南四湖最重要的子流域之一,其土地利用/覆被变化研究可为流域周围土地合理利用提供理论基础,同时也对南四湖流域的水资源保护和生态管理具有重要意义。

2.2 数据来源

本研究需要的数据为分辨率30m的Landsat8 OLI-TIRS TM遥感影像和行政区划图。运用ENVI5.1和ArcGIS10.2软件对图像数据进行处理,再采用ENVI5.1软件对遥感影像进行解译,从而得到分辨率为30m×30m、100m×100m的1990年、2000年、2005年、2010年、2015年的五期土地利用类型图。根据我国土地利用分类系统和白马河流域实际土地利用情况对土地利用类型进行划分,共分为耕地、林地、草地、水域、建设用地、未利用土地六大类。

3 研究方法

3.1 土地利用动态度

土地利用类型变化速率是反映区域土地利用变化剧烈程度的重要指标,一般可分为单一土地利用动态度和综合土地利用动态度。动态度表示该土地利用类型的稳定情况[23],单一土地利用动态度模型的公式为:

式中,k表示研究时段内某一土地利用类型的动态度;At1、At2分别表示评测初期和末期各种土地利用类型的面积;t1、t2分别为变化初期和末期。

3.2 CA-Markov模型

CA-Markov模型既包括了元胞自动机模型模拟复杂空间变化的能力,也综合了马尔科夫模型长时间的预测优势,可更精确地对土地利用类型变化进行模拟预测[30]。本研究中,像元可视为元胞自动机中的元胞,土地利用类型可视为元胞状态,像元大小即为元胞的尺度大小,并利用滤波器对元胞邻域进行定义[31]。通过Markov模型得到的两期土地利用状态转移概率矩阵、转移面积矩阵和状态转移图集作为转换规则进行预测下一时刻元胞状态的变化趋势,并对研究区域的土地利用/覆被变化进行模拟预测,具体步骤为:①数据格式转化。本文在利用前期获得的1990年、2000年、2005年、2010年、2015年五期土地利用分类图的基础上,为了方便将其导入IDRISI软件中进行操作分析,转换为软件可用的数据格式(.rst的格式),包括矢量数据转为栅格数据，栅格数据转为ASCII文件，ASCII文件转为IDRISI软件识别的数据格式。②由于CA-Markov模型具有时空尺度敏感性,同时为了兼顾研究意义与效率,本文选取像元为30m和100m的两种规格研究CA-Markov模型的空间尺度效应。CA滤波器是用来创建具有显著空间意义的权重因子,作用于元胞以此来确定元胞下一时期的改变状况,本文选用CA标准的5×5的滤波器。即确定中心元胞周围为5×5的元胞矩形空间,迭代次数分别为5次和10次。③获取转移面积矩阵。将两期土地利用类型图导入IDRISI软件中,利用Markov模块得到该期间内土地利用转移面积矩阵、转移概率矩阵和状态适宜图集,其中转移面积矩阵作为转换规则参与下一步模拟与预测,将模拟得到的土地利用图与实际图进行一致性检验;以模拟精度较高时的元胞规格和迭代次数为条件对2020年的土地利用类型进行模拟预测。④模型预测的精度检验。为比较两年土地利用图的相似程度,利用Kappa系数[32]进行一致性检验,Kappa指数计算公式为:

式中,Pa为正确模拟的比例;Pc为随机情况下期望的正确模拟比例;Pp为理想分类情况下的正确模拟比例,为100%。

本文将模型运行的迭代次数与其两年间的时间步长设为一致。根据Kappa系数检验标准[33,34],当Kappa系数为0—0.20时,表明一致性极低;当Kappa系数为0.21—0.40时,表明一致性一般;当Kappa系数为0.41—0.60时,表明一致性中等;当Kappa系数为0.61—0.80时,表明具有高度的一致性;当Kappa系数为0.80—1.0时,表明两图像几乎完全一样。

4 模型检验与模拟分析

4.1 土地利用现状分析

由ArcGIS10.2和ENVI5.1软件分析所得到的五期土地利用类型图,根据我国土地利用分类系统和白马河流域实际土地利用情况进行土地利用类型的划分,共分为六类土地利用类型,处理后的最终结果见图2。从图2可知,耕地面积广布研究区域,草地和林地主要分布在研究区域的南部和东南部,东北部很少,而西部和北部几乎没有分布;建设用地面积主要占用其周围的耕地而快速增加;未利用地在逐渐减少,水域面积变化较小。

图2 白马河流域五期土地利用类型

4.2 模拟精度检验

为了更加精确地预测研究区白马河流域2020年的土地利用,首先针对该流域进行初步模拟精度检验分析,选择一个较合适的时空尺度条件进行预测。本文对不同的时间步长分类进行了土地利用预测模拟。当时间步长即迭代次数设置为10次时,将1990年与2000年30m×30m的土地利用类型图导入IDRISI软件中,获得1990—2000年期间内土地利用转移面积矩阵、转移概率矩阵和状态适宜图集,将转移面积矩阵作为转换规则参与2010年土地利用的模拟与预测,得到元胞为30m×30m的2010年土地利用预测图。同理,将1990年与2000年100m×100m的土地利用类型图导入IDRISI软件进行相同操作,得到元胞为100m×100m的2010年土地利用预测图。当时间步长即迭代次数设置为5次时,同上述步骤,依次将2005年与2010年元胞为30m×30m、100m×100m的土地利用类型图导入IDRISI软件,利用转移面积矩阵分别得到2015年元胞为30m×30m、100m×100m的土地利用模拟预测图(图3)。首先，将2010年的土地利用模拟结果与实际2010年的土地利用数据在Database Query-CROSSTAB模块中进行一致性检验,其中元胞为30m×30m、100m×100m所得到的Kappa系数值分别0.8876、0.8877;其次将2015年的模拟结果与实际2015年的土地利用数据进行一致性检验,其中元胞为30m×30m、100m×100m所得到的Kappa系数值分别为0.9121、0.9068。由此可以得出在元胞为30m×30m、迭代次数为5次的条件下对2015年进行模拟预测时得到的Kappa系数最高为0.9121。由Kappa系数检验标准可知,在此条件下CA-Markov模型预测模拟的2015年土地利用图和实际几乎完全一致,即此时空尺度为该研究区域预测条件中较适宜的模拟条件。

图3 2015年、2010年白马河流域土地利用预测

4.3 土地利用特征分析

由精度检验分析可知,分辨率为30m×30m、迭代次数为5次时预测精度达到最大,因此本文选择在此条件下进行后续模拟分析,只对与该条件相关的几期(2005年、2010年、2015年三期)土地利用进行特征分析。

土地利用面积变化:2005年、2010年、2015年三期的土地利用/覆被变化图经过ArcGIS10.2统计分析,得到流域三个时期的土地利用类型面积与变化动态度,变化动态度即可计算各土地利用类型面积在某段时期内的变化速度(表1)。由表1可知,在2005—2015年研究期间,研究区域2005年、2010年、2015年不同土地利用类型的面积表现为耕地>建设用地>草地>水域>林地>未利用地,耕地是白马河流域的主要土地利用类型(2005年、2010年、2015年分别占66.82%、66.87%、65.63%)。土地利用格局在该期间发生了巨大变化。整体看,耕地、林地、草地和未利用地均减少。其中,草地、林地减少较明显。草地从2005年的446.57km2减少到2015年的328.38km2,减少幅度达26.48%;林地从2005年的114.83km2减少到2015年的92.31km2,减少幅度达19.67%。其次为未利用地和耕地,减少幅度分别为9.95%和1.93%。建设用地和水域逐渐扩张,建设用地呈现出十分强劲的上升趋势,变化幅度在2005—2010年达28.2%,2010—2015年增加了6.91%;水域扩张不明显,十年间仅增加了0.24%。

表1 三期土地利用面积

从整体来看(表1),各地类的土地变化速度相差较大。耕地在研究期间虽有增有减,但总体变化动态度较小;林地在整个研究期间呈减少趋势,其中2005—2010年减少了19.52%,2010—2015年减少较缓慢;草地主要在2005—2010年减少且变化动态度较大,2010—2015年减少相对缓慢;水域整体增长较缓慢;建设用地急剧增长,在2005—2010年增加最快,动态度高达28.20%,2010—2015年增加相对缓慢;未利用地一直在减少,但变化动态度相对较慢。

土地利用结构的空间转换分析:本文主要利用IDRISI软件中的Markov模块分别对2005—2010年和2010—2015年的土地利用类型图进行转移矩阵分析,得到转移面积矩阵见表2、表3。由表2、表3可知,研究区域总体呈现出“两增四减”的变化规律。即建设用地和水域呈现增加趋势,耕地、林地、草地和未利用地呈现减少趋势。水域变化不明显,主要转入地类为耕地、草地和建设用地,耕地转入占优,主要转出地类为耕地和建设用地。2005—2015年建设用地一直在增加,总体净转入576.34km2,转入类型主要为耕地,占转入的98%以上;转出地类同样主要为耕地,但转出的耕地面积远远小于耕地的转入面积。2005—2015年耕地一直在减少,转出地类主要为建设用地和水域。其中,2005—2010年和2010—2015年耕地转为建设用地,分别占总面积的25.23%和17.80%;转入地类主要为草地和建设用地。林地转换不平衡,2005—2010年净转出27.50km2,2010—2015年净转入0.07km2,转出地类主要为耕地、草地和建设用地,其中主要为耕地和林地;转入地类同样主要为耕地和草地,尤其在2010—2015年草地转林地面积为11.04km2。草地转出地类为主要耕地,2005—2010年草地转为耕地的面积占其总转出面积的90%以上;草地转入地类同样主要为耕地。未利用地总体变化趋势并不太明显,但转入转出地类仍主要为耕地。

表2 2005—2010年土地利用转移面积矩阵(km2)

表3 2010—2015年土地利用转移面积矩阵(km2)

4.4 2020年土地利用模拟结果分析

本文选择在元胞为30m×30m、迭代次数为5次的条件下，以2015年的土地利用数据为基础对2020年的土地利用变化进行预测模拟。利用IDRISI软件中的Markov模块得到2010年和2015年两期土地利用类型图的转移概率矩阵、转移面积矩阵和状态适宜性图集。以2015年为基础,以状态适宜性图集为转换规则,迭代次数设置为5次(利用2010年和2015年预测2020年,时间间隔为五年),利用其转移概率矩阵对2020年土地利用进行预测模拟。2015—2020年的土地利用转化与2020年的土地利用类型见表4、表5和图4。

表4 2015—2020年土地利用转移面积矩阵(km2)

表5 2015—2020年土地利用类型变化与面积比

图4 2020年土地利用变化状况

由表4和表5可知,2015—2020年土地利用变化呈现出“三增三减”的变化规律,即耕地、林地、草地减少,水域、建设用地、未利用地面积增加。耕地变化幅度高达6.675%,面积比由65.63%减少到58.95%,主要转为建设用地和水域,其中转为建设用地的面积为488.95km2,占耕地总转出的95.18%;转入地类主要为草地和耕地周围的林地。水域变化幅度较小,转出地类主要为耕地、林地、建设用地和未利用地,其中林地占优;转入地类主要为耕地和建设用地。草地面积比由9.18%减少到8.60%,转出主要为耕地、林地和建设用地,其中转为耕地的面积为48.58km2,占草地面积总转出的77.59%;转入地类主要为水域,变化幅度较小。建设用地变化幅度较大,面积比由18.26%增加到25.18%,转出地类主要为水域;转入地类主要为耕地、草地和水域,其中耕地转入的面积占建设用地总转出的97.96%,变化幅度和动态度都较大。林地和未利用地没有明显变化,变化幅度较小。

5 结论与讨论

本文利用白马河流域不同时期的遥感影像,以ArcGIS10.2、ENVI5.1和IDRISI17.00为计算平台,结合CA-Markov模型,对研究区域2020年的地类进行了模拟预测。由于CA-Markov模型具有时空尺度敏感性,因此首先将该模型进行模拟所需的元胞规格和时间步长设为不同,分别对2010年和2015年地类进行预测;其次，将预测结果分别与实际年份的土地利用类型进行对比,得出模拟精度较高时的元胞规格与时间间隔;第三,在次较适宜的时空尺度基础上,预测研究区域2020年的地类。研究分析表明:①以Markov模型中获得的两期土地利用状态适宜图集作为模型转换规则具有可行性,CA-Markov模型在预测白马河研究区域的土地利用类型时具有很好的适用性。②本文对比分析了像元大小为30m×30m,时间尺度5年;100m×100m,时间尺度5年;30m×30m,时间尺度10年;100m×100m,时间尺度10年共四种情况,分别在CA-Markov模型中进行模拟预测,并将预测结果与实际年份土地利用图进行一致性检验。对2015年的土地利用/覆被变化进行模拟预测时得到的Kappa系数最高为91.21%,由此得出在本文研究中像元大小为30m×30m、时间尺度为5年为较适宜的研究模拟条件,模拟精度最高。③在研究期间内,耕地、建设用地和草地为研究区域的主要用地类型,三者合计占白马河流域面积的90%以上。总体看,耕地、林地、草地和未利用地均呈减少趋势,其中草地减少幅度最大,减少了31.06%;建设用地和水域呈增加趋势,建设用地增加最多为47.10%。2005—2020年,白马河流域建设用地面积呈现出急剧增加的趋势,由2005年的476.66km2增加到2020年的901.09km2,主要是随着经济的发展城市化速度加快,城市人口增加,导致建设用地大幅度扩张,建设用地主要是占用周围耕地面积而增加,耕地转为建设用地占建设用地总转入的94.30%。耕地面积虽然在2005—2010年有轻微增加,但总体减少了11.78%,除了被建设用地大量占用外,耕地同样转出为水域,占耕地总转出的7.57%,耕地为了弥补巨大的减少面积而占用建设用地和周围的草地。草地有减少趋势,减少了31.06%,由2005年的446.57km2减少到2020年的307.88km2,草地主要被耕地占用。林地、水域、未利用地并无大幅度变化。

土地利用/覆被变化是一个极其复杂的变化过程,为了使预测结果更加准确,应综合考虑自然因素和人为因素对地理空间系统的影响。在今后的研究中可尝试建立不同的影响因子和不同决策者共同作用的影响因素,加入权重要素建立不同类型的土地利用/覆被变化模型。CA-Markov模型具有时空尺度敏感性,本文仅对比了四种不同时空尺度,选择较适宜的条件对2020年的土地利用进行预测分析,对比分析的时空尺度条件较少。今后研究中可尝试对时空尺度条件进行更细致的对比分析,以得出可靠性、准确性更高的研究结果。CA-Markov模型综合了元胞自动机和马尔科夫链各自的优势,提高了模拟的精度,该模型的预测结果较好地反映了白马河流域的土地利用变化趋势,可在一定程度上对该流域的土地利用/覆被变化动态监测和可持续发展提供参考。