大汶河流域土地利用变化特征及预测

沈振涛，赵 强，高 倩，武 玮，刘玉玉

(济南大学 水利与环境学院，山东 济南 250022)

土地利用/土地覆被变化(LUCC)作为社会、经济与自然环境相互作用的表现形式[1]，在全球环境变化、可持续发展研究中占有重要地位[2]。LUCC模型是探索未来土地利用演变过程的特征规律，揭示人类活动影响下的生态环境变化的重要方法[3]。

目前，用于土地利用变化模拟的方法有很多，常见的有利用变化及效应(CLUE-S)模型[4]、马尔可夫(Markov)预测模型[5]、元胞自动机(CA)-Markov模型[6]、系统动力学模型[7]、逻辑回归(LR)模型[8]、人工神经网络(ANN)[9]等，其中的CA-Markov模型被公认为对模拟和预测土地利用变化具有重要价值[10]。CA-Markov模型可以整合社会、经济、地形和道路数据并进行模拟，从而反映社会经济发展的特征[11]，因此，许多学者使用该模型来预测和优化未来的土地利用，例如：周浩等[12]基于1990、2013年挠力河流域的土地利用数据，对2025年土地利用格局进行模拟预测；侯西勇等[13]基于1990、2000年河西走廊土地利用数据，模拟了20世纪90年代该地区土地利用变化；靳含等[14]运用多时间跨度模拟四道河子镇未来的土地利用变化。目前CA-Markov模型的应用较为广泛，且在土地利用演化模拟中发挥了重要作用。

本文中以大汶河流域为研究对象，选取1985、1995、2005、2017年遥感影像，利用完整的遥感图像处理平台(ENVI)进行解译，结合土地利用变化的面积和土地利用动态度指标分析研究区内1985—2017年土地利用的动态变化特征，最后运用CA-Markov模型对研究区未来的土地利用格局进行模拟预测，为大汶河流域土地利用资源的合理运用及生态环境保护提供参考。

1 研究区概况

大汶河又被称为古汶水，位于山东中部，是黄河下游水系中最大且最后的一条支流。大汶河发源地为山东旋崮山北麓，自东向西流经山东省境内7个县(市)[15]。大汶河流域位于东经116°11′15″—118°0′0″、北纬35°37′30″—36°32′30″之间，东部以鲁山为界，西部以黄河流域为界，南部以蒙山及余脉和淮河流域为界，北部主要以泰山山脉和小清河为界，流域全长209 km，面积8 634 km2。平均海拔为213 m，地势北高南低，东高西低，南北两面环山[16]。大汶河流域地理位置如图1所示。

2 数据来源与研究方法

2.1 数据来源

采用1985、1995、2005年Landsat 4-5 TM数据和2017年Landsat 8 OLI_TIRS数据对研究区进行土地分类。为了便于进行土地分类，选取每年6月份且云量较少时的数据作为基础数据，运用ENVI软件对遥感数据进行辐射纠正、几何校正等处理，参照GB/T 21010—2007《土地利用现状分类》，以及遥感影像分辨率、土地利用方式、土地利用实际情况和后期模型分析的需要，将研究区土地利用类型分为6种，分别为建筑用地(URHD)、林地(FRST)、未利用地(BARR)、水体(WATR)、草地(PAST)、耕地(AGRL)。

2.2 研究方法

2.2.1 土地利用变化特征指标

净转移量是土地转入与转出面积之差，总转移量是指各地类在某时间段的转入和转出面积之和，相关计算公式分别为

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式中：Nc为转移净变化率，%；ΔUin为其他地类转化为该地类的总面积，km2；ΔUout为该地类转化为其他地类的面积，km2；Ua、Ub分别为该地类在时间段内初始、末期面积，km2；Tc为转移总变化率，%；K为单一土地利用类型的动态变化率，%；T为时间，a；Ps为转移状态；Sa为总土地利用净转移变化率，%；ΔUin,i为其他土地利用类型转化为类型i的总面积，km2；ΔUout,i为土地利用类型i转化为其他类型的面积，km2；St为转移总变化率，%；Lc为综合土地利用动态度；ΔUai,j为某时间段内土地利用类型i转移为类型j的面积，km2；ΔUai为土地利用类型i的总面积，km2；n为土地利用类型总数；Pt为某时间段内的转移状态。

2.2.2 CA-Markov模型

Markov模型的原理是基于栅格的空间概率，根据现有状况预测未来各个时刻的变化[17]，已被广泛应用于土地利用变化模拟[18]、景观格局变化模拟[19]和森林和植被的可持续性模拟[20]。利用遥感与地理信息结合应用系统(IDRISI Selva)中的CA-Markov模块，以1995、2005年土地利用状况为初始矩阵，得出1995—2005年间大汶河流域土地利用面积转移矩阵并将其作为模型转移规则，验证2017年土地利用状况，然后以2005、2017年土地利用状况为初始矩阵，得出2005—2017年间土地面积转移矩阵，并预测研究区2030年土地利用状况。

3 大汶河流域土地利用变化特征

3.1 土地利用结构变化特征

根据大汶河流域1985、1995、2005、2017年共4期土地利用状况，基于ArcGIS的空间分析功能，对研究区4期的土地利用类型图进行处理，同时提取大汶河流域1985、1995、2005、2017年的4期土地利用数据，如图2、表1所示。

表1 大汶河流域不同时期土地利用类型面积对比

从相关数据可以看出，林地、耕地和草地是研究区主要的3种土地利用类型，覆盖了研究区的绝大部分地区。1985年，耕地面积所占比例最大，达到35.59%，面积为3 072.88 km2；之后是草地、林地，所占比例分别为28.12%、24.05%，面积分别为2 427.72、2 076.73 km2；建筑用地占比为7.73%，面积为667.78 km2；水体、未利用地所占比例相对较小，分别为1.21%、3.30%，面积分别为104.20、284.74 km2。在空间上，草地、耕地、林地交叉分布，林地遍布研究区各地，林地主要位于丘陵山区，建筑用地则主要分布于河流附近及平原地区，而耕地位于中间位置，主要依附水体且与建筑用地同存，未利用地主要零散分布于丘陵山区未开垦处。

1995年，草地面积所占比例最大，达到37.01%，面积为3 195.31 km2；之后是耕地，所占比例为26.57%，面积为2 427.72 km2；林地、建筑用地占比分别为19.96%、10.69%，面积分别为1 723.63、923.33 km2；水体、未利用地所占比例分别为1.36%、4.40%，面积分别为117.81、379.71 km2。与1985年相比，草地面积增加了8.89%，耕地和林地面积分别减少了9.02%、4.09%，建筑用地面积增加了2.96%，未利用地和水体面积变化相对较小。

2005年，各土地利用类型中所占比例最大的依然是草地，达到43.71%，面积为3 773.67 km2；林地和耕地次之，所占比例分别为20.95%、19.64%，面积分别为1 808.53、1 695.54 km2；建筑用地占比分别为13.91%，面积为1 201.01 km2；水体、未利用地所占比例分别为1.32%、0.48%，面积分别为114.27、41.04 km2。与1995年相比，草地面积持续增加，增加率为6.7%；耕地面积持续减少，减少率为6.93%；林地、建筑用地面积分别增加了0.99%、3.22%；未利用地面积减少了3.92%，水体面积变化相对较小。

2017年，耕地面积所占比例最大，达到38.42%，面积为3 317.32 km2；草地次之，所占比例为32.16%，面积为2 776.88 km2；之后是林地、建筑用地，所占比例分别为11.84%、16.20%，面积分别为1 022.22、1 398.63 km2；水体、未利用地所占比例较小，分别为1.34%、0.04%，面积分别为115.9、3.09 km2。与2005年相比，草地面积减少了11.55%，耕地面积增加了18.78%，林地面积减少了9.11%，建筑用地面积增加了2.29%，未利用地面积减少了0.44%，水体面积变化依然不大。

为了直观分析大汶河流域近30年来的土地利用变化，从不同土地利用类型在各年中所占比例(表1)可知：1985—2017年间，研究区建筑用地面积逐年增加，从667.78 km2扩大到1 398.63 km2；林地面积只在2005年有微弱增加，总体呈减少趋势，从2 076.73 km2减少到1 022.22 km2；耕地面积先减少后增加，与1985年相比，2017年耕地面积增加了244.44 km2；草地面积变化与耕地相反，呈先增加后减少趋势，与1985年相比，2017年草地面积仍增加了349.16 km2；未利用地面积逐渐减少，而水体面积无明显变化。

3.2 土地利用动态转移过程

在近35 a间，大汶河流域存在几十种土地利用转移关系，其中绝大部分来自建筑用地、林地、耕地以及草地的转移(见表2)。由表可见：1985—2017年间，建筑用地主要转化为耕地，转化面积为186.43 km2，转化比例为28%；林地主要转化为耕地，转化面积为670.14 km2，转化比例达到32%；林地向草地转化的面积为354.04 km2，转化比例为17%；未利用地主要转化为草地、耕地，转化面积分别为161.41、87.34 km2，转化比例分别为58%、31%；耕地主要转化为草地，其次转化为建筑用地,转化面积分别为1 102.33、536.76 km2，转化比例分别为36%、17%；草地则主要转化为耕地，转化面积为1 038.96 km2，转化比例为43%。水体面积的转化量相对较小。

表2 1985—2017年大汶河流域土地利用变化转移矩阵 km2

根据式(1)—(8)计算可得到大汶河流域1985—2017年土地利用动态度和土地利用转移态势，如表3所示。由表中数据可知，大汶河流域在1985—2017年间，土地利用动态度Lc为0.91，土地利用转移状态Pt为0.27，流域整体土地利用属于弱平衡型转换；水体、耕地以及草地Ps分别为0.13、0.07和0.11，说明水体、耕地和草地呈缓慢扩张趋势；建筑用地Ps为0.52，说明建筑用地呈扩张趋势，且变化趋势较为明显。林地Ps为-0.73，说明大汶河流域林地呈显著减少趋势；未利用地Ps为-0.98，说明未利用地显著减少，呈逆向转换。

表3 1985—2017年大汶河流域土地利用动态度和转移状态

3.3 土地利用驱动力分析

人口与经济都是城市化发展进程的主要影响因素。在人类发展过程中，为了满足生产、生活的需要，人们对土地的开发和利用改变了当地土地利用类型和分布。大汶河自东向西流经莱芜、泰安市，最后注入东平湖。2005年泰安市总人口为549.95万，到2017年全市总人口达到570.70万，13年间人口增加20.75万；2005年莱芜市总人口为124.29万，2017年全市总人口达到129.68万，13年间人口增加5.39万。流域内人口数量的增加使大汶河流域土地利用格局发生变化，建筑用地面积不断扩大。经济发展过程实际上是土地利用结构的优化调整与土地利用集约度提高的过程。2017年，泰安、莱芜市固定资产分别达到2 992.5、568.14亿元，经济的快速发展，使得市区建筑面积不断扩大，计算结果显示，2005—2017年研究区建筑用地增加16.45%，草地和未利用地转化为建筑用地的比例分别为11.98%和7.58%。政策因素是短期内推动土地利用类型变化的关键因素之一。1987年的《泰安市林业区划》中规定，山区建设要以水源涵养林为主，使林地的种植面积得到了增加。2006年的《莱芜市土地利用总体规划》中规定，要优化土地利用空间格局，加强中心城区土地利用控制，增强农用地保护和基本农田建设，推进建筑用地节约和集约利用等，对莱芜市土地利用变化起到了控制作用。

4 大汶河土地利用变化预测

4.1 驱动力因子选取

驱动力因子是影响土地利用模拟的重要因素，因此选取适宜的驱动力因子对模型的模拟非常重要[21]。结合研究区经济、生态条件，参考相关研究结果[22-24]，本文中选取了自然及社会因素共10种驱动力因子(见表4)。A1、A2分别为高程和坡度，代表自然因素；A3—A10代表社会因素，其中A4—A9

表4 元胞自动机-马尔科夫(CA-Markov)模拟采用的驱动力因子

表示道路到达的可能性，A10代表与河流的距离，表示植被获取水资源的可能性。利用ArcGIS软件将各类影响因子与土地利用进行统计分析，对土地利用的影响进行归一化处理，结果见图3。

4.2 适宜性图集

土地适宜性是土地在一定条件下对不同用途的适应程度[25]。在模拟土地利用转化过程的适宜性时，选取了以上10个驱动力因子。利用逻辑回归计算土地利用状况跟驱动力因子之间的关系，然后依据回归关系得到全部6种土地利用类型的适宜性图集(见图4)。该过程利用IDRISI Selva软件进行模拟。

适宜性图集中的因变量是展开的像元，以二进制栅格表示，其中1表示特定像元的变化，说明可转化为该类地；0表示这一段时间内没有变化，说明可转化为其他类地[26]。从图4中可以看出，耕地和草地转化的适宜性变化最大，而建筑用地、林地及未利用地转化的适宜性变化较小，水体以保持原样为主。

4.3 模型精度

Kappa系数可以用来检验模拟结果与实际数据的一致性，因此可用于2个图件一致性评价和遥感解译精度评价[27]。Kappa系数K的计算公式为

(9)

式中：Po为模拟正确的栅格比例；Pc为期望模拟正确的栅格比例；Pp理想状态下模拟正确的栅格比例；N为栅格总数；N1为模拟正确的栅格数。

基于1995、2005年土地利用数据，利用CA-Markov模型预测大汶河流域2017年土地利用状况，并将2017年土地利用类型图与模型的模拟图对比，进行CA-Markov模型精度检验。大汶河流域2017年土地利用类型图共有栅格个数为9 586 156。

利用IDRISI软件中的CROSSTAB模块计算模型精度，得到Kappa系数为0.705，表明模型可以较好地模拟大汶河流域2017年土地利用状况，且精度较高，可以用来模拟研究区未来土地利用状况。

4.4 未来土地利用预测

以2017年土地利用状况为基准，利用2005—2017年土地利用转移矩阵以及完成的适宜性图集，利用CA-Markov模型，循环13次，预测大汶河流域2030年土地利用状况，结果如图5所示。

将1985—2030年间共5期土地利用类型图中各类地所占比例进行对比，结果如图6所示。与2017年相比，大汶河流域草地面积有所增加，耕地、未利用地以及水体面积基本不变，林地面积大幅度减少，而建筑用地有所增加。未来研究区草地面积呈先增加、后减少的趋势，到2030年又有所增加；耕地面积则正好相反，先减少、后增加，但2017—2030年间建筑用地和耕地面积逐渐增加；未利用地面积在1995—2005年间大幅度减少，直到2030年变化都不大；未来大汶河流域水体面积基本保持不变。预测结果显示，未来一段时间内研究区建筑用地增加有所放缓，说明通过政府的政策影响对莱芜市土地利用变化起到了控制作用，因此要严格贯彻执行《莱芜市土地利用总体规划》，积极响应国家号召，加强中心城区土地利用控制，加大农用地保护和基本农田建设、保护生态环境可持续发展的力度，推进建筑用地节约和集约利用。

大汶河流域地图从国家测绘地理信息局标准地图服务网站下载，地图审批编号为GS(2019)1697(http://bzdt.ch.mnr.gov.cn/browse.html?picId=%224o28b0625501ad13015501ad2bfc0258%22)，经过Arcgis 10.7软件数字化处理后得到。图5 大汶河流域2030年土地利用状况预测

图6 1985—2030年大汶河流域土地利用类型变化

5 结论

本文中以大汶河流域为研究对象，通过分析1985—2017年土地利用变化情况、影响因子的驱动作用进行分析，并结合CA-Markov模型对2030年研究区土地利用变化趋势进行了模拟预测，得到主要结论如下：

1)通过对1985年以来大汶河流域土地利用变化进行分析，水体、耕地和草地呈缓慢扩张趋势，建筑用地呈明显扩张趋势，林地呈显著减少趋势，未利用地呈逆向转换。究其原因，人口数量增加和经济发展使得大汶河流域土地利用格局发生变化，城市化进程处于不断扩张阶段。

2)对大汶河流域2017年土地利用进行验证，Kappa系数为0.705，说明利用CA-Markov模型可以较好地模拟大汶河流域2017年土地利用状况，且精度较高，可以用来模拟未来土地利用状况。

3)运用CA-Markov模型对大汶河流域2030年土地利用变化进行模拟，结果显示，草地面积逐渐增加，耕地、未利用地以及水体面积基本没有变化，林地面积大幅度减少，而建筑用地有所增加，表明政府的政策调控措施优化了中心城区土地利用格局，对城市土地利用变化起到了一定的控制作用。