燕山北部农牧交错区土地利用变化分析与模拟

刘雪莹， 郝弯弯，赵鹏，马佳明，董大鹏，温哲华，谷建才

(1河北农业大学 林学院，河北 保定 071000；2河北省林业和草原信息中心，河北 石家庄 050081)

土地利用变化是人类经济社会发展对土地资源利用的直接反映[1]，土地利用方式受自然和人文因素的综合影响[2]。通过LULC土地利用类型的相关变化情况，可将它与人类的关联、自然的发展以及形态等清晰地体现出来[3]。实际上，该问题对于环境、资源以及人口等的可持续性发展来说是非常重要的，并涉及到土地的利用强度以及覆被等方面的变化情况[2]。虽然它在各个地区中所产生的具体作用及影响机制还有待深入认识，但已逐渐演变成影响全球土地变化的主要因素之一[4]。人类(文化，社会经济和政治)与生物物理环境在不同空间尺度中相互作用，是研究区域环境变化的重要途径[5]，目前受到各国学者的广泛关注。自21世纪以来，在全球范围内由于城市化和智能化战略所带来的人地关系的发展，俨然已经成为一种新趋势，人类活动主导了全球环境的演变[5]。该研究区属于干旱、半干旱气候的牧区和半湿润气候的农业区之间的过渡区域[6]。实际上，从整个陆地生态系统的角度审视可知，农牧交错区是一类具有比较高敏感性的地区，更容易受到全球变化和人为因素的共同干扰，从而导致了一系列的生态问题，例如环境污染、生物多样性丧失、不合理的土地利用方式以及生态环境随着各种影响因素逐渐恶化[4]。耕地—林地、耕地—草地、林地—草地的交替变化对区域生态环境具有重要影响，人地关系矛盾更为突出。人类的高强度活动，如过度开垦、放牧、挖掘和遗弃等，都会给地区带来了巨大的负面效应[7]。气候变暖、人口和粮食需求增加等破坏了天然植被，导致草地的退化，不断地耗尽了地区生态系统服务功能[8]。因此，有效解决区域发展带来的土地问题，正确处理人地关系是保证城市经济健康稳定发展的关键，也是城市建设与发展的重要保证。因此，需要借鉴较为成熟的内陆土地利用研究成果，进一步探索研究区的土地利用强度。

基于此，以燕山北部农牧交错区——围场县为实例研究对象，借助遥感和GIS 等空间分析手段，对2009—2014—2019年土地利用方面的具体变化情况展开理论模拟和深入的分析。

1 研究区概况

研究区地处在燕山余脉与蒙古高原之间的交接区域，位于大兴安岭的南部地区。所处地理位置为北纬39°40′～42°10′，东经115°45′～119°50′[9]。山势陡峭，地势西北高，东南低。山体呈东西走向，海拔高度为500～1 500 m。该研究区有比较低的水土容量，其土壤中含有相对比较低的有机质，在其东、西部的相关区域中，分别是黑色森林土以及栗钙土的组成[10]。研究区处于暖温带大陆性季风气候区。年均温6 ℃～10 ℃，1月均温-12 ℃～6 ℃，7月为20 ℃～25 ℃。10 ℃以上持续期195～205 d，活动积温2 600 ℃～3 800 ℃，年降水量700 mm左右，具有较为明显的时空分布方面的差异性。

植被茂盛，灌木、杂草丛生，森林面积广阔。地带性植被为落叶阔叶林(以栎类为主)，并混生暖性针叶油松林，垂直带谱[11]。700 m以下为落叶阔叶林，树种有蒙古栎、槲栎、栓皮栎、槲树等。700～1 500 m为针阔叶混交林，树种以臭冷杉为主。1 500～2 000 m为针叶林，树种有华北落叶松等，但以次生林为主。

2 数据及方法

2.1 数据来源及信息提取

研究区主要采用2009年、2014年、2019年3期 Landsat 系列影像作为遥感数据源。通过地理空间数据云(www.gscloud.cn/)来进行卫星遥感相关图像信息的完整下载。尽量选择植被生长较为旺盛的 7—9 月份，并且云层覆盖度不高于5%[12]或研究区范围内无云的数据，确保相关的影像数据资料中无条带以及噪声等现象[13]。除了从各种来源获得的辅助信息外，还借助全国森林资源二类调查数据，评估从卫星数据中提取的土地覆盖图的准确性。

2.2 土地利用分类

特征识别被认为是土地利用分类工作中十分重要的一个环节。由于土地利用类型方面的差异性，将会在遥感图像中形成各自的特征光谱。基于不同的光谱信息可以进行不同土地利用类型信息的辨别[14]。充分考虑当前我国在该领域中的工作基础以及研究区的具体情况，对该地区土地利用类型进行分类，具体如下(共计6类)，即耕地、林地、草地、水域、居民建设用地、未利用地6种基本的土地类型[15]。其中，旱地与普通农田是耕地的主要地类，林地通常为疏林地、灌木林、常规林地以及除此之外的其他类型的林地，而草地通常为放牧场以及天然的草场，水域则是湖泊和河流等，居民及生产用地包括各类建筑区域(如居民区以及城乡用地等)，除此之外的其他尚未利用的土地则主要为各类废弃土地、裸地、盐碱地以及沙地等。

2.3 土地利用动态度变化

土地利用动态度是真实反映土地利用类型动态变化程度的指数。它可以表达特定区域一定时间范围内某种土地利用类型变化情况[16]。

(1)

式中：K为单一土地利用动态度；Ua表示初期土地利用面积；Ub分别为末期土地利用面积；T表示变化时间段。

2.4 CA—Markov模型

在Markov 模型中，作为一个具有随机性特征的运动过程，状态与时间均具有离散性。该模型主要侧重于通过时间尺度来开展各土地利用类型的变化情况[17]。在预测土地利用所对应的状态时，采用的公式具体为：

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

(2)

式中：Pij代表了状态转移的相关矩阵，而S(t)与S(t+1)则分别代表对应于t和t+1时的研究土地所对应的具体结构状态。

CA 模型则是空间维度上进行预测[18]，其状态改变的规则在时间和空间是都是局部的，每个变量只取有限状态。其公式表示如下：

S(t+1)=f(S(t),N)

(3)

式中：f代表了局部空间中所对应的元胞转化方面的有关规则；N代表了元胞邻域；而S(t)与S(t+1)则分别代表了对应于t和t+1等2个时刻相关的元胞有限离散所对应的集合状态。CA模型中邻域结构的设置，系统默认设置为5×5，如下图所示:

本研究所采用的CA-Markov 模型则很好地融合了上述模型所具有的一系列优势，从而可从空间及时间等角度来科学预测及模拟相关地区的土地利用情况，从而对未来土地利用格局的发展提供依据[19]。该研究是以2019年为起始年，预测的时间间隔为5 a，并且进行CA循环次数(5)、CA滤波器所对应的大小(此处设为5×5)以及元胞的具体大小(此处设为30 m× 30 m)等参数的合理设置，具体步骤有以下几点：

1)数据格式的转换：再ArcGIS中将栅格数据转化呈IDRISI识别的ASCII格式的数据；

2)重分类：对转化后的数据进行重新分类，数值根据已有地类编码的情况而定；

3)生成转移矩阵：在Markov模块生成转移面积矩阵和概率矩阵；

4)模拟预测：利用Markov预测所得的转移矩阵进行CA-Markov的预测。

2.5 模拟预测精度验证

运用Kappa系数进行精度检验是现在众多学者广泛使用的一种验证方法(式1)。本研究采用IDRISI 18.3软件中CROSSTAB模块计算各地类Kappa系数以验证其空间模拟精度。精度的判别标准为 Kappa 值从0～1 表示2幅图像从完全不吻合到完全吻合的转变[20]。具体标准为：

表1 精度的判别标准Table 1 Discriminant criteria for Kaap coefficient accuracy

公式表达如下 ：

(4)

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

3 结果分析

3.1 土地利用格局及其变化

2009年、2014年及2019年土地利用现状及其变化如表2所示。分析其中的数据可知，对于当前的研究区来说，其占据主要地位的土地利用类型为林地，在整个研究区中，达到了61.70%～63.51%的比例，而且有相对较广的分布情况，占据主导地位，且集中连片，呈现出较小的破碎程度[21-22]；此外，在整个研究区域里，耕地及草地在总体面积中均有超过10%的占比水平，耕地面积广，散乱分布，随着时间的更替变化，耕地面积具有持续降低的趋势；在该区域的西北部，主要分布的土地利用类型是草地，2009年、2014年、2019年草地类型的面积占比依次是17.5%、16.65%、16.74%，呈波动变化的趋势，草地流失现象时有发生；居民生产用地2009年、2014年、2019年面积占比依次是1.50%、1.64%、1.74%，面积由2009年的114.31 hm2，增加到2019年的230.24 hm2，呈现出以居民区为中心的向外扩张趋势，说明在城镇化水平提高、人口压力增大、经济不断发展的背景下，人们对建设用地扩展的需求增加[23]；水域面积在2009年、2014年、2019年间面积所占比例分别为3.20%、2.82%、3.92%，未利用地面积所占比例分别为1.80%、1.47%、2.94%，水域面积和未利用地面积均处于先减少后增加趋势，变化浮动明显。

表2 2009—2019年研究区土地利用面积及所占比例Table 2 Area and percentage of land use types in 2009, 2014 and 2019

3.2 土地利用类型面积转移情况分析

从土地利用类型的面积转移情况可以清楚地看到每种土地利用类型转移为其他土地利用类型的面积概率。2009—2014及2014—2019年土地利用面积转移概率矩阵分别见表3和表4：

1)林地在2009—2019年呈现持续增加的态势，2009—2014年林地面积增加速度高于2014—2019年。2009—2014年林地净新增量来源依次为草地(61.76%)、耕地(32.63%)、水域(39.77%)、居民生产用地(30.41%)、未利用地(4.34%)，草地是最大的转入者。2014—2019年林地净新增量来源依次为耕地(22.05%)、草地(21.33%)、居民生产用地(7.05%)、水域(6.74%)和未利用地(5.26%)，说明退耕还林项目取得效果，并持续进行中。

2)耕地呈现出先增加后减少的态势，2009—2014年净增加面积增加了0.07%，在从2014年到2019年的时间区间里，出现了大规模的面积减少，林地为其主要的转化方向，而且有一部分的尚未利用开发的土地以及草地的转入。

3)草地随着时间的演替呈现了一种逐渐退化的趋势[24-25]，2009—2014年净减少面积为50.37 hm2，2014—2019年略有增加，具有35.18 hm2的净增量。最为重要的转出类型为林地，相比之下，草地则具有更小的面积。

4)居民建筑用地面积由2009年的114.31 hm2，增加到2019年的230.24 hm2，在2009—2019年期间呈现只增不减的趋势，而未利用地且大部分是在原有城区和居民点的基础上进行扩张。

5)2009—2019年水域面积呈现先减少后增加的波动状态，新增土地利用类型分别对应于生产以及生活所用土地、各类耕地、草地、尚未得到利用开发的土地以及林地。水域的动态变化主要受自然环境和人为因素的干扰[16]。

6)未利用地也呈现出少量增加趋势，通过林地以及生产生活所用土地的有效转入，形成了新增利用土地。

表3 2009—2014年土地利用转移概率矩阵Table 3 Land use transfer probability matrix during 2009—2014

表4 2014—2019年土地利用转移概率矩阵Table 4 Land use transfer probability matrix during 2014—2019

3.3 土地利用动态度变化

土地利用动态度模型可以定量描述研究区不同时间段不同土地利用类型数量变化的幅度和速度。研究区土地利用用动态度变化见表5。

表5 研究区土地利用动态度Table 5 Land use dynamic attitude in the study area %

由表5可知，(1)耕地在2009—2014年间的动态度为0.60%，增长速度缓慢，2014—2019年间的动态度为-0.84%，这段时间经济快速发展，城镇用地呈现向外扩张的趋势，其次退耕还林相关生态项目的开展，占用了大量的耕地。

(2)林地在2009—2014年间的动态度为1.56%，增长速度较快，2014—2019年间的动态度为-0.53%，增长速度较缓慢；草地在2009—2014年间动态度为-0.80%，变化速度缓慢，2014—2019年间动态度为0.54%，呈现缓慢增加的趋势，说明近些年来政府对环境保护的重视，生态项目在有效进行。

(3)居民生产用地2009—2014年间的动态度为2.91%，是6种土地利用类型中速度最快的，2014—2019年间的动态度为0.76% ，呈缓慢变化的趋势，畜牧业带动当地经济的发展，建筑用地不断向外扩张。

(4)水域在2009—2014年间的动态度为-2.05%，2014—2019年间动态度为7.26%。水域的动态度变化受自然环境和人类生产生活影响较大，同时水利设施的变化也对水域动态度变化有一定影响。

(5)未利用地在2009—2014年间动态度为-3.51%，2014—2019年间动态度为19.04%，动态波动明显，政府应加强土地的有效利用。

4 CA—Markov模型预测结果

4.1 预测精度验证

采用CA-Markov模型，以2014年为起始年，时间间隔为5 a，预测2019年的土地利用情况，并将预测结果与实际结果进对比，结果如图1。

图1 CA-Markov模型精度预测对比图Figure 1 Comparison chart of CA-Markov model accuracy prediction

在预测图和实际图对比的基础上，通过IDRISI 18.3软件中CROSSTAB模块将实际数据(见图a)和模拟数据(见图b)展开Kappa 验证，如表6所示。发现Kappa的值为0.912 2，这说明它具有比较好的精度水平，所以可继续开展后面的预测工作。

表6 Kappa系数预测结果Table 6 Kappa coefficient prediction results

4.2 预测结果

以2019年土地利用数据为基础，基于CA-Markov模型对研究区的2024年进行土地利用类型的预测，预测结果见图2。

图2 2024年研究区土地利用预测图Figure 2 Land use prediction map of the study area in 2024

结果图2表明：研究区2024年较2019年变化显著，林地面积仍占主导地位，集中连片；耕地和草地分布散乱，有减少趋势，草地退化现象加剧了未利用地面积的增加，居民生产用地增加显著，以居民点向外扩张，城镇化水平提高，对建筑用地的需求增加。

表7 2009—2024年研究区土地利用面积变化Table 7 2009—2024 change in land use area of the study area

同2009年、2014年及2019年3个时段相比：

1)2024年耕地面积减少，呈现波浪状的变化趋势，2014年耕地面积较其他年份最多，此后相继减少，这与当前农业政策直接相关。

2)林地面积净减153.32 km2，波动不太明显，且相对集中，新增林地面积小于林地转移面积，由于林地面积较大，林地的损失和收益强度基本处于休眠状态。

3)草地面积呈现出逐年下降的趋势，与2019年相比下降了118.49 km2，过牧、滥垦、滥樵和滥挖药材等人为因素的干扰，导致草场退化。

4)居民生产用地持续不断增加，居民安置点和城区的持续扩大就是具体的体现。耕地以及林地是建筑收益的最重要来源。

5)尚未得到利用及开发的土地出现了持续性的面积扩大，它和林地、草地以及耕地间的类型转换相对活跃，水域面积大量增加，在原来水域面积的基础上进行了扩宽和延长。

5 结论

本研究以典型农牧交错区围场县为例，结合3s和遥感技术相关原理，进行2009年、2014年和2019年土地利用类型的解译。并通过CA-Markov模型来预测2024年该区域的土地利用变化情况，由此可促进当地经济和社会的良性发展，尤其是科学有效地利用土地资源提供重要的基础信息。结果显示：

(1)不同时期遥感图像通过解译显示出各土地利用类型的结构分布。研究区的土地利用类型呈现出“三增三减”的变化趋势。其主要的土地利用类型是草地、耕地以及林地，3期林地面积比例均超过60%，耕地面积均在10%，草地、居民生产用地、水域、未利用地所占份额较小。

(2)转移概率矩阵表明：林地收支相对平衡；耕地呈现出先增加后减少的态势，它们的主要转化为草地和林地，其中林地达到32.63 %，而且林地和耕地两者之间存在着一定程度的转化。

(3)根据模拟，2009—2024年该研究区有着较为显著的土地利用类型方面的改变趋势，围场县的林地和耕地面积出现了逐渐下降的趋势，居民生活点和各类建筑用土地将会逐渐以居民区为中心向外持续进行扩展，是最活跃的土地利用类型，它不同于耕地和草地，这直接关系到当地社会和经济发展水平的提高以及城市化进程的加快。

6 讨论

(1)从动态的视角把握物质流动的规律，优化土地资源的利用方式，调控经济系统资源流动的良性路径，已经逐渐成为学术界研究的热点问题[26]。因此，土地覆被变化也成为了全球环境变化的重要内容，但大多数侧重于宏观尺度，对于县域的研究还相对较少，为此本研究区以2009年、2014年、2019年遥感影像和全国森林资源二类调查数据对该研究区土地利用类型的变换规律进行分析。高凤杰[6]等在东北山地丘陵地带农林交错区的研究中表明在这几类土地利用类型中，变化最明显的是建筑用地的持续增加，这与本研究的观点是一致的。

(2)本研究所采用的CA-Markov模型有效地融合了CA模型和Markov模型所具有的优点。Markov模型可有效地完成较长时间范围内的预测，以及CA模型在系统较为复杂时可进行空间变化的模拟，所以备受关注。已经有许多科研工作者通过该模型来对黄土塬区[27]、黄土丘陵[28]及籍河流域[29]等地区进行土地利用变化模拟，并取得良好的效果。研究区域是典型的农牧交错区，土地利用覆被变化剧烈，有研究表明这与产业结构，土地适度开发等政策息息相关[2]。同时，受到技术以及研究时间方面的限制，关于经济以及社会等层面的政策影响的探讨还相对欠缺[12]，如何精确量化及调整是今后研究的重点，在未来的研究中，应加强对这方面的研究，以提高预测结果的准确性[21]。