耦合GMOP与PLUS模型的干旱区土地利用格局模拟

2023-02-16马瑞范燕敏武红旗朱磊马轩范昕张皓宇

农业资源与环境学报 2023年1期

马瑞，范燕敏*，武红旗，朱磊，马轩，范昕，张皓宇

（1.新疆农业大学资源与环境学院，乌鲁木齐 830046；2.中南财经政法大学法学院，武汉 430073；3.中国地质大学公共管理学院，武汉 430074；4.新疆大学资源与环境科学学院，乌鲁木齐 830046）

土地利用/土地覆被（Land Use/Land Cover，LULC）是反映人类社会和自然环境之间相互依存的基本指标，在全球环境变化过程中发挥了重要作用[1-3]。随着社会经济发展及城市化进程不断加快，人类对自然的改造极大程度上改变了LULC，对生态环境、粮食安全、生物多样性等造成了一定威胁[4-5]。因此，兼顾经济、生态效益的土地利用优化，有利于合理高效利用土地资源，从而实现人与自然和谐发展，保障土地可持续发展[6-7]。

国内外诸多研究将国土空间优化配置作为热点，有研究学者采用线性规划、综合平衡法、遗传算法等对土地利用需求数量进行优化[8-9]。但简单的定量分析技术无法对土地利用进行全面系统分析和动态定量研究[10]，因此有学者将模型构建法与智能算法结合，开辟了土地利用需求数量和空间分布格局优化的新途径[11-12]。模型构建法是通过对土地利用变化影响因素的相关性分析、回归分析等，获得各用地类型在空间上的分布概率，继而实现土地利用空间优化配置，其实质主要是对未来土地利用模拟[13]。目前，基于模型构建的优化配置方法主要有多智能体系统（MAS）、最小累积阻力模型（MCR）、元胞自动机模型（CA）等[14]。MAS 可以将一个庞杂的土地利用系统分割成多个彼此间相互交流、协调的子系统，从而将复杂的问题简单化，但是MAS 的驱动因素层较多，寻优程序繁杂，工作量大[15-17]。MCR 能够较好地反映景观格局变化和生态过程演变的相互作用关系，适用于中尺度的土地利用功能分区研究，但是其阻力值的设定具有较高的主观性[18-19]。CA 模型是一种时间、空间、状态均离散且这三者关系为局部的网格动力学模型，在模拟系统的空间格局变化方面具有强大能力，但CA模型紧密的结构与固定的转换规则使其无法准确解释土地利用与驱动因素之间的相互关系，更无法描述土地利用的成因与过程[16]。本研究采用的PLUS 模型原理源自CA，并在传统的CA模型基础上进行了较大改进，结合了转化分析策略（TAS）和格局分析策略（PAS）优势，既能揭示土地利用变化潜在驱动因素，也能很好地模拟多种土地利用的斑块演变[20-21]。

昌吉回族自治州（简称“昌吉州”）作为建设陆上丝绸之路经济带的重要节点城市之一，经济发展较快，但由于其地处干旱区，存在水土资源配置不平衡、生态较脆弱等问题，主要表现为：南部山区森林系统和部分荒漠地区生态系统功能恶化；草原和湿地面积减少；过度开采地下水资源，导致水位降低；沙漠边缘的生态系统不稳定；南部山区的小煤矿占地面积大，生态破坏严重；林地、草原、水源涵养区等生态环境保护困难[22]。因此，探讨如何科学合理规划土地利用，实现经济社会与资源环境协调发展尤为重要。鉴于此，本研究以昌吉回族自治州为研究区，探究2010—2020年的LULC规律，构建2030年昌吉州土地利用分布格局。为了使模拟结果分布合理并达到研究区经济发展、环境保护、资源有效利用的长期平衡[23]，本研究结合社会经济统计数据，利用PLUS-GMOP 耦合模型，实现昌吉州土地利用在数量和空间上的优化配置情景模拟，指导决策者制定未来土地利用总体规划，合理开发利用土地资源，为干旱地区实现区域社会经济和生态环境协调发展提供参考。

1 材料与方法

1.1 研究区概况

昌吉州位于天山北麓、准噶尔盆地东南缘，总面积7.39×106hm2[24-25]（图1）。该州地处欧亚大陆腹地，地势南高北低，由于受地形地貌、地域降水的影响，南北气候差异较大，南部夏季降水较多，北部沙漠性气候特征显著，是典型的大陆性干旱气候[26]。全州土地利用的总格局地域分明，利用方式比较合理。南部是天山山地，主要以牧业、林业、生态保护区为主；中部是冲积平原，以农业用地、城镇为主；北部是沙漠盆地，为牧业及工矿利用区。截至2020 年末，全州总人口139.94 万人，区域生产总值1 387.25 亿元，其中农林牧渔业生产总值246.2亿元。

图1 研究区地形示意图Figure 1 Topographic map of the study area

1.2 数据来源及处理

研究数据集主要包括：①昌吉州2000、2010、2020 年土地利用数据：依据国家土地利用分类标准和本研究目的，将土地利用类型重分类为耕地、林地、草地、水域、建设用地和未利用地六大类；②土地利用限制因素：将新疆西天山和天池区域作为土地利用转换禁区进行二值化处理；③气候和环境数据：驱动土地利用变化的自然因素主要有地形、气候、地下水，其中河流水系数据进行欧氏距离分析处理，坡度数据利用DEM 提取得到；④社会经济数据：分为驱动土地利用变化的栅格数据和计算多目标土地利用数量结构的统计数据，栅格数据主要包括人口空间分布、国内生产总值（GDP）、夜间灯光数据等，统计数据主要涉及乡村城镇人口、GDP、粮食播种面积和产量等。数据来源见表1。

表1 数据类型及来源Table 1 Data types and sources

1.3 研究方法

本研究基于2010、2020 年昌吉州土地利用数据和社会经济、自然数据，探究2010—2020 年的土地利用变化规律，通过马尔可夫链（Markov chain）和GMOP 模型对土地利用需求进行拟合和优化，并利用PLUS模型模拟预测自然发展情景和多目标优化情景下昌吉州土地利用分布格局（图2）。

图2 基于GMOP-PLUS耦合模型的情景模拟研究框架Figure 2 Scenario simulation research framework based on GMOP-PLUS model

1.3.1 Markov chain

Markov chain 是一种随机级数，通过一段时期内土地利用类型之间的转移概率来预测未来土地利用类型需求总量[27]。若一个Markov chain 的状态空间是有限的，将所有状态的转移概率按矩阵排列，可得到下列转移矩阵：

式中：n为土地利用类型数量；Pij为i类土地利用类型转化为j类的概率，由概率分布的性质可得0≤Pi，j≤1，且每行元素之和等于1。Markov chain 完全决定转移矩阵，转移矩阵也完全决定Markov chain。

1.3.2 土地利用变化动态度

土地利用动态度表示某一时期研究区不同土地利用类型的变化速率和变化幅度[28]。

式中：K为研究区某一土地利用类型的动态度；Ua、Ub为研究初期和末期某种土地利用类型的面积；T为研究时长。

1.3.3 灰色多目标线性规划模型

灰色多目标动态规划模型（GMOP）是由灰色预测理论与多目标线性规划交叉而衍生的优化体系，其主体思想在于将不确定的目标函数和约束条件定义在某一确定范围的满意区域内[29]。GMOP 模型充分表达了决策者的意愿并提供决策者一种最佳的土地利用优化配置方案[30]。

1.3.3.1 构建目标体系

决策者在优化土地资源配置时，是以因地制宜可持续发展的思想角度来构建目标函数。目标体系包括经济发展目标、生态保护目标和社会效益目标[31-33]。考虑到社会效益难以量化，本研究从生态价值和经济效益两个方面建立目标函数。

最大生态价值目标函数为：

式中：Z1为昌吉州内生态系统服务功能价值总量，万元；gi表示某一土地利用类型的生态系统服务价值系数，万元·hm-2；xi为土地利用类型i的规划面积，hm2。

本研究采用谢高地等[34]构建的“中国陆地生态系统单位面积生态服务价值系数表”对相关统计年鉴数据进行分析，得到2010—2020年单位面积粮食经济产值，并通过灰色模型预测，得到适用于昌吉州目标年份且与土地利用类型相对应的生态系统服务价值系数（表2）。将生态系统服务价值系数代入式（3）得到昌吉州生态价值目标函数：

表2 2030年昌吉州不同土地利用类型单位面积生态价值系数Table 2 Ecological value coefficient table of different land use types per unit area in Changji in 2030

最大经济效益目标函数为：

式中：Z2为研究区内经济效益总值，万元；Ci为经济效益系数，万元·hm-2；xi为土地利用类型i的规划面积，hm2。

本研究用单位面积土地所产出的生产总值来表示其经济效益系数Ci[35]。结合《昌吉回族自治州统计年鉴》，将农林牧渔业生产总值中农业、林业、畜牧业、渔业产值作为耕地、林地、草地、水域的产出效益，而建设用地的经济效益则以第二、三产业的GDP 来表示，计算得出2010—2020年内各地类经济效益系数，并利用灰色模型进行预测，得到2030 年耕地、林地、草地、水域、建设用地的经济效益系数分别为1.34 万、0.18万、0.29 万、0.48 万、165.32 万元·hm-2。借鉴相关研究[36]，本研究将未利用地的经济效益设为0.000 1 万元·hm-2，将经济效益系数代入式（4）最终可得昌吉州经济效益目标函数：

1.3.3.2 约束条件设置

约束条件参数的选取准确与否，决定了模拟结果是否可行，也决定了接近离散时动态系统的静态模型能否真实地反映系统自身规律[37]。

（1）土地总面积约束：各决策变量的土地利用面积之和为研究区总面积（7 388 542.17 hm2），即：

（2）人口总量约束：农用地和城镇用地承载的人口总数应控制在预测年份（2030年）人口总量之内。

式中：D1、D2分别为昌吉州目标年份农用地及城镇用地的平均人口密度，人·hm-2；系数0.71为2010—2020年城镇面积占建设用地总面积比例的平均值；P为昌吉州2030年总人口预测值。根据2010—2020年人口及土地利用结构数据并利用灰色预测模型得到：0.25(x1+x2+x3)+7.51×0.71x5≤1 586 593

（3）粮食安全约束：

式中：⊗()a31为灰色粮食单产系数，kg·hm-2，其上限为目标年份粮食单产，下限为现状粮食单产；f0为目标年份规划的复种指数；fr为粮食系数；s0为人均标准用粮（根据国家规定的居民食物消费和营养基本目标，人均粮食消费量为450 kg·a-1）[38]；p0为目标年份预测总人口数。通过对昌吉州国民经济统计资料的分析，得到：⊗（a31）∈（5 750.19，7 880.50），f0=0.55，fr=0.48，s0=450，p0=1 586 593。

（4）农用地面积约束：农业用地变化应保持动态平衡，所以农业用地面积总和应大于或等于现状值，即：x1+x2+x3≥3 409 009.74。

（5）草地面积约束：考虑昌吉地区草地大幅度减少的历史趋势，依据《昌吉回族自治州土地利用总体规划》提出的合理安排牧草地、加强人工草地建设，本研究假设多目标情景下草地面积应低于现状草地面积、高于自然发展情景下草地面积。但在政策影响下，草地面积的减少速率正在下降，Markov chain预测的2030年各地类面积是基于2010—2020年变化趋势得到的，因此多目标情景下草地面积减少速率将低于Markov chain 预测的草地减少速率，因而可得：x3≥2 018 674.71;（2 207 856.69-x3）÷2 207 856.69≤8.57%。

（6）建设用地面积约束：昌吉州作为一带一路核心区的重要组成部分，经济快速发展的同时城镇用地规模不断扩大。因此本研究将Markov chain 预测的2030 年建设用地面积的120%、80%分别作为多目标发展情景优化下的建设用地面积上限和下限[20]，即：116 713.23 ≤x5≤175 069.80。

（7）耕地面积约束：为追求经济效益，昌吉州等干旱区耕地面积长期呈现持续扩张状态，导致灌溉用水量随耕地扩张而迅速增加，但由于受到水资源限制，耕地的扩张保持在一定范围内，即：1 040 624.01≤x1≤1 084 542.21。

（8）林地面积约束：“十三五”期间，昌吉州实施以生态建设为主的林地保护战略，建设“南护天山，北治沙漠，中建绿洲”的生态格局，故目标年份林地面积应高于现状，另外，将Markov chain 预测得到的林地面积的120%设为目标年份的上限，可得：160 529.04≤x2≤185 222.38。

（9）水域面积约束：2011 年，昌吉州各县（市）加大湿地保护力度，在坚决保护现有湿地、努力扩大湿地面积基础上，全面完成州内19 条规模以上河湖水域岸线保护利用规划。因此，将2011 年水域面积设置为优化方案下目标年份水域面积的下限，将Markov chain预测的自然发展情景下水域面积的120%设为上限，即：54 228.69≤x4≤67 434.44。

（10）未利用面积约束：昌吉州北部覆盖大面积沙漠，长期以来，干旱少雨的自然气候以及周边草场的过度放牧，使本就脆弱的生态系统出现了较为严重的破坏现象，另外大风也导致昌吉州土地沙漠化不断蔓延。“十三五”以来昌吉州统筹推进山水林田湖草沙冰一体化保护和系统治理，全力推进防沙治沙工程。因此本研究假设未利用地面积低于自然发展情景，高于2010 年未利用地面积，即：3 668 162.76≤x6≤3 928 886.37。

（11）目标约束：假设昌吉州实现土地经济效益、生态价值同时增长，即：

式中：13 693 350.34 为2020 年现状土地经济效益，9 331 105.86为2020年现状土地生态价值总量。

（12）决策变量的非负约束：xi≥0，i=1，2，3…6。

1.3.3.3 GMOP模型求解

通过Lingo 软件对函数进行求解，最终得到6 种土地利用类型的面积，即为2030 年多目标优化情景下的土地需求数量。

2 结果与分析

2.1 土地利用变化特征

在过去10 年间，随着昌吉州社会经济的快速发展，土地利用方式发生了很大改变。根据2010—2020 年昌吉地区土地利用转移矩阵（表3）可知，土地利用发生变化的总面积为3.67×105hm2，占地域总面积的4.96%。各地类之间均存在互相转变，其中草地转出数量最多，草地向其他地类转入面积从高到低依次为未利用地、耕地、建设用地、林地、水域，其中草地退化为未利用地的面积占总退化面积的58.76%；未利用地面积增加最为明显，增加了1.48×105hm2，主要原因是草场退化，其中草地退化成未利用地面积为1.63×105hm2，占草地总面积的6.69%。；耕地和建设用地大幅度增加，但主要都是以草地转入为主，同时耕地和未利用地均不同程度地向建设用地转变；各类型土地净转变量中，水域和林地的变化最不明显，但转出、转入均以草地为主。近10年昌吉州土地结构转变主要是草地退化为未利用地，这是由于不合理的开发利用和大面积的破坏植被使土地受到风沙侵蚀，生产力丧失。

表3 2010—2020年土地利用类型转移矩阵（hm2）Table 3 2010—2020 land use type transfer matrix（hm2）

2.2 模型模拟精度检验

在PLUS 软件中，以昌吉州2000、2010 两年土地利用类型数据为基础，输入自然和社会经济等驱动因素，得到研究区2020 年土地利用模拟图，将其与土地利用现状图（图3）对比发现，南部山区和绿洲区模拟结果较好；虽然中北部荒漠区和东南部山区部分区域有偏差，但各土地利用类型在空间分布上整体一致。模拟精度验证中，Kappa系数为0.810 7，总体精度达到0.882 0，表明PLUS模型能够较好地反映昌吉地区土地利用变化的实际情况，进而也可以认为PLUS模型适用于昌吉州未来年份的土地利用模拟与优化。

图3 2020年昌吉回族自治州土地利用现状图与模拟图Figure 3 Land use status map and simulation map of Changji Hui Autonomous Prefecture in 2020

2.3 2030年土地利用情景模拟结果

2.3.1 基于PLUS单模型的自然发展情景

本研究设定的自然发展情景假设2020—2030 年间土地利用变化不受政策干扰，这使得历史土地利用变化的趋势得以保持。利用Markov chain 的概率转移矩阵计算2030 年各土地利用类型需求面积。为了防止昌吉州自然保护区生态系统被破坏，将新疆西天山和天池设置为限制转化区域，通过设置随机森林参数、模拟参数、各类用地需求并运行PLUS 模型，得到在自然发展情景下2030年昌吉州的土地利用模拟图。

在自然发展情景下，2030 年昌吉地区的土地利用分布格局延续了现有变化趋势（图4，表4）。草地和林地面积持续减少，草地累计净减少189 181.98 hm2，约占草地总面积的8.57%，林地土地利用动态度为-0.384 8%，林草地变化面积主要集中在南部山区（图4A2）。其他3种地类面积均出现增加现象，其中未利用地面积增加最多，较2020 年增加了113 136.75 hm2，主要由草地转入为主，部分沙丘出现活化现象。近年来，水域面积虽呈增加趋势，但其土地利用动态度变化却不明显，约为0.142 2%。昌吉州经济的快速发展使得建设用地变化最为剧烈（图4B2），其土地利用动态度为3.461 7%，增加面积主要来源于草地、未利用地，说明在该情景下，部分未利用地被开发为建设用地，确保了昌吉州土地经济效益从2020 年的1 369.34 亿元增长到2030 年的2 621.26 亿元，增幅高达91.43%。同时大量占用草地、破坏草地生态的人类活动，使昌吉地区土地生态价值出现了负增长，生态效益值相比2020 年减少1.40%（表5），土地荒漠化进一步发展，土地生态安全并不乐观。

表4 昌吉州2020—2030年不同情景下土地利用面积及土地利用动态度Table 4 Land use area and land use dynamics degree under different scenarios in Changji from 2020 to 2030

表5 昌吉州2020—2030年不同情景下土地经济、生态效益值与变化率Table 5 Benefit value and change rate of land economic and ecological benefits under different scenario in Changji from 2020 to 2030

图4 2020年和2030年不同情景下土地利用类型分布的局部细节图Figure 4 Map of local details of land use type distribution under different scenarios for 2020 and 2030

2.3.2 基于GMOP-PLUS耦合模型

在经济与生态兼顾情景下，2030 年草地面积依旧呈减少趋势，但速率有所减缓，减少速率为3.11%。未利用地面积较2020 年减少了41 372.62 hm2，主要转化为耕地和建设用地。耕地面积增加了43 917.99 hm2，主要以草地和未利用地转入为主，其中耕地以散点式分布状态占用未利用地（图4C3），决策者可根据土地实际情况，在散点式分布范围内集中连片开发未利用地为耕地。建设用地虽延续历史趋势呈现扩张趋势，但其扩张程度得到控制（图4B3），土地利用动态度从2010—2020 年的5.726 4%降到2.707 7%，并以昌吉市、阜康市为中心，众多小城镇为基点进行分布。水域面积的增加幅度仅次于建设用地，其土地利用动态度为2.170 7%。林地面积增加了24 693.36 hm2，在南部山区增加效果明显（图4A3）。优化后的土地经济效益达2 490.77 亿元，相比优化前增加了81.90%；土地生态价值达985.49 亿元，较优化前增加了5.61%，实现了土地经济效益和生态价值的双提升。

2.3.3 PLUS与GMOP-PLUS模拟结果对比

对比分析自然发展情景下模拟结果与多目标优化情景下模拟结果可知：从土地利用转化来看，相比于PLUS 单模型模拟结果，多目标优化情景下林草地面积表现为增加趋势且主要集中在南部山区（图4A2、图4A3），草地和林地面积分别增加了120 570.29、30 870.42 hm2。其中林地面积增加的主要原因是当地规定严禁林地向其他地类转出，这很大程度上促进了土地生态价值的增加。建设用地面积减少了8 170.83 hm2，有效减缓了对草地的占用。在多目标情景下，未利用地面积减少了154 509.37 hm2，说明其较大程度地被开发利用，提高了土地利用率，有效防止了农用地荒漠化。

从土地效益来看，自然发展情景下土地经济效益发展飞速，增速高达91.43%，这得益于建设用地的迅速扩张，但与此同时，生态价值却出现了负增长（-1.40%），说明在自然发展情景下昌吉州追求土地经济发展的同时破坏了生态环境；而在调整优化土地利用结构情景下，通过对未利用地的开发利用、对建设用地的管控以及林地和草地的稳步增加，昌吉州土地荒漠化控制在一定范围内，土地生态环境明显改善，既保证了土地经济效益的增加，又提升了生态价值。

3 讨论

根据自治区实施“一棱二轴多组群”和围绕乌鲁木齐打造五大经济圈的总体部署，昌吉州凭借其优越的地理区位和资源优势，具有很大的城市化快速发展潜力。这给原本脆弱的生态环境带来了更大挑战，而预测该区域在未来不同发展情景下土地利用变化及其对生态系统可能产生的影响将有助于应对这一挑战。

与自然发展情景相比，2020—2030 年期间，多目标优化情景下未利用地将被大量开发，生态用地转化为城市用地的面积减少。更重要的是，优化后的用地方案更具可持续性和弹性，因为其不仅考虑了研究区未来10 年“养活”和“容纳”不断增长的人口所需的土地，而且林地、草地和水域面积也满足相关森林和水资源保护规划的要求。新疆西天山和天池风景名胜区为保护区，在多目标发展情景下其用途不会改变。优化方案结果表明，更多的未利用地被开发为城市用地，使经济效益得到提高。而防止草地退化和保护水资源、林地资源等措施将有助于建设一个更加生态友好的城市。但本研究也存在一定不足，如在设置PLUS模拟参数时，需要结合专家知识和区域实际情况，通过连续测试最终确定参数值，这具有一定主观性。