随着城镇化进程的推进, 我国土地利用结构发生了极大变化, 许多自然土地, 尤其是林地和草地被进一步开垦成耕地、建设用地。随着研究领域不断发展成熟, 出现了针对海岸带^[1]、河流流域^[2]、喀斯特地区^[3]等极具地方特色的特殊地域、特殊对象的土地利用模拟研究, 但较为缺乏针对不同分区之间的过渡地区的土地利用模拟研究。现阶段, 随着国土空间规划与主体功能区规划的制定与实施, 各大分区之间定位合理、分工明确, 但对于地处不同功能分区之间的过渡与交汇区域而言, 各类用地空间之间的权衡与竞争关系得以加强, 如何合理配置区域土地资源、兼顾“三生”空间需求已成为亟待解决的问题, 因此本文着眼于地处由生态示范区向经济发展区过渡、高原特色牧业区与粮油主产区交汇的龙门山过渡带区域, 以期合理配置土地资源、提出更好协调农业生产、城镇化发展与生态保护之间的关系的土地利用方案, 为类似的山地-平原过渡区域、山周盆地地区社会经济发展与生态环境保护协调发展、生态保护与经济发展双线并行提供思路与参考。

土地利用模拟当今国内外学术界对于土地利用变化预测模拟的相关研究领域发展得相当成熟, 方法理论上有系统动力学、多元统计学、指数分析法、地理加权回归、遗传算法、蚁群算法等^[4—6]; 模型工具上有元胞自动机模型(CA)、Markov模型、地理模拟优化系统(GeoSOS)、小尺度土地利用变化及其空间效应模型(CLUE-S)、FLUS模型等^[7—9]; 已取得的重大成果上有黎夏团队的高精度全球未来100年模拟产品、宫鹏团队的全球未来土地利用预测数据集、George C. Hurtt团队的Land Use Harmonization数据集等^[10]; 土地利用模拟与研究区实际情况的适配程度不断提高, 模拟对象不再局限于城市、大城市群或平原地区, 研究视角不再局限于单纯的用地模拟, 土地利用模拟与国土空间、生态安全、“三生”空间等进行交汇^[11], 形成多角度、多层次、多维度的土地利用模拟研究, 近年来随着研究的深入, 针对不同情景、不同用途、不同研究对象的多模型耦合的土地利用模拟模型得以研发并应用^[12], 进一步提高了模拟精度。

“过渡带”这一概念被广泛应用于地质学、生物学、气象学等许多方面的研究中^[13], 如生态过渡带、降水量过渡带、农业牧业过渡带等, 主要表现为多种因素在一定区域内具有一定空间规律的交错分布。龙门山地处四川盆地与青藏高原东缘交界处, 跨越了我国第一、第二级阶梯^[14], 无论是地形地势、气候条件还是社会经济条件都十分特殊且非常具有代表性, 龙门山及其周围地区呈现出典型的山地-平原过渡、都市-乡村过渡、人口密集区-人口稀疏区过渡的情形, 因此, 本研究致力于研究龙门山过渡带的土地利用变化规律、预测土地利用演变趋势, 不仅有利于合理开发利用土地资源, 促进土地资源可持续发展, 也能为情况类似的山地-平原过渡区、都市-乡村过渡区的土地利用发展与规划提供参考与思路。

随着城镇化与工业化进程的发展, 粗放式的土地利用方式导致的生态环境恶化、土地利用矛盾加剧, 如何统筹“三生”空间, 如何优化国土空间结构已成为需要回答的重大问题, 未来土地利用的权衡与配置优化研究已成为研究热点。从土地利用优化研究方法上看, 现有优化方法主要分为数量结构优化与空间结构优化, 数量结构优化研究主要方法主要有线性规划、多目标规划、系统动力学模型等, 空间结构优化主要通过CLUE模型、CA模型、多智能体系统(MAS)、仿生智能算法(如遗传算法、蚁群算法等)实现^[11], 总体看来目前的土地利用权衡优化研究计算量较大, 计算耗时较长且权衡判定的自动化程度较低。本研究根据研究区实际情况与数据, 通过设定多重判定条件, 对待权衡斑块进行多层次逻辑判定, 提出了权衡优化的程序化逻辑。

1 研究区概况及研究方法 1.1 研究区的界定与概况

研究区地处102°—107°E, 30°—34°N的四川省北部偏东地区(图 1), 共约九万平方千米。地形单元由西至东大致为川西高原山地区-成都平原区-川东北低山丘陵区, 从主体功能上看, 研究区自西向东呈现明显的生态保护区-经济发展区的过渡趋势; 在气候条件方面, 龙门山以东地区受季风影响强烈, 降水集中在夏季, 年降水量在800—1200mm之间, 年平均气温约15—18℃, 昼夜温差相对较小, 而研究区西部气候呈现较明显的垂直地带性, 汶川、平武等县区海拔相对较低, 年均温约为10—15℃, 而地处更西部的松潘、小金、黑水等海拔高的区县年均温仅为5—10℃, 同时降水受地形影响较大, 研究区西部汶川、平武以西区县年平均降水量不足800mm, 日间太阳辐射强度大, 昼夜温差大。总体看来, 研究区气候条件自西向东呈现明显的高原高山气候-亚热带季风气候的过渡趋势。

研究区土地利用地区发展较不均衡, 如西部耕地资源在数量上与质量上均显不足, 农业生产空间受到制约, 东部地区城镇空间与建设用地扩展迅速, 挤压其他用地空间等。除此之外, 研究区内还存在数量众多、占地面积广阔的自然保护区、森林公园、生态优先保护区等各级生态保护地^[15—16], 随着研究区内各类用地的发展与竞争, 研究区生态压力增大, 如何较好地协调经济发展与生态保护之间的矛盾, 权衡农业生产空间、城镇化发展空间与生态保护空间三者之间的关系已成为研究区未来国土空间规划面临的重要问题。

1.2 主要数据来源及预处理

土地利用变化是在不同时空范围内自然条件和人类活动相互作用产生的结果^[17], 参照FLUS模型模拟土地利用变化的相关研究^[18—20], 选取自然因子、社会经济因子、可达性因子一共9个因子。土地利用数据来自中科院资源环境与数据中心(www.resdc.cn), 包括2010年、2015年与2020年三期四川省土地利用栅格数据, 初始分辨率为30m×30m(其他主要数据及其来源详见表 1)。考虑到研究区实际尺度同时为方便FLUS模型计算及模拟, 将各类数据分辨率统一为100m×100m, 行列数统一为4853×4218。

1.3 研究方法 1.3.1 SRP模型

生态敏感性-压力度-恢复力模型(SRP)主要分为敏感性、压力度、恢复力三个部分, 本次生态敏感性计算指标选取高程、地形起伏度、坡度、坡向、气温、干旱度、降水、土壤侵蚀度, 生态恢复力选取植被覆盖度与植被净初级生产力(NPP)两个指标, 生态压力度选取国内生产总值(GDP)、耕地面积、人口密度三个指标进行评价^[21]。根据各个指标的具体性质对各个指标分别进行标准化处理, 并进行主成份分析, 依据主成分分析结果, 生态脆弱性指数计算公式如下：

(1)

式中, E表示生态脆弱性指数; m_i是第i个主成分所对应的贡献率; P_i为第i个主成分; n是累计贡献率为85%以上的前n个主成分。

1.3.2 Markov模型

采用马尔可夫链来预测未来像元总量参数, 即未来的各土地利用类型所占像元的总量。马尔科夫链可以通过两期土地利用数据生成土地利用转移矩阵, 计算出研究区域内土地利用变化转移概率矩阵, 在此基础上, 预测出研究区未来各土地利用类型的数量。其公式如下^[22]：

(2)

式中, S表示t时刻的土地利用状态, S_t+1表示t+1时刻的土地利用状态, P_ab表示土地利用类型a转变为土地利用类型b的概率。

1.3.3 FLUS模型

FLUS模型以历史土地利用变化数据为基础, 基于神经网络元胞自动机模型(ANN-CA)提取土地利用变化规律, 结合土地利用变化驱动因子数据模拟未来各类用地的数量与空间分布^[23]。通过神经网络训练可以产生较高概率的分布适宜性, 建立各种土地类型的发生概率与驱动因素之间的关系, 神经网络的模拟过程是建立起初始土地利用类型与各个驱动因子在空间上相互作用的关系^[24], 公式如下：

(3)

式中, α(m, a, t)表示土地利用类型a于时间t出现在栅格m上的概率, w_{j, k}与sigm(net_i(m, t))分别表示隐藏层与输出层间的权重系数与关联函数, net_i(p, t)是在i隐藏层像神经元i在栅格m的t时间发送的信号, w_{i, j}表示输入层与隐藏层之间的信号。

2 基于生态综合评价的限制转换区划定

本研究在进行生态服务功能重要性评价的同时结合了生态脆弱性的评价, 利用生态综合评价分级的方法, 将研究区划分为不同的生态保护需求度等级, 并以其等级为根据对不同模拟情景设置不同的限制转换区。

2.1 生态服务功能重要性评价

生态服务功能重要性评价的方法与指标参考了《生态保护红线划定技术指南》(环办生态〔2017〕48号)文件, 基于2020年现状数据进行计算, 考虑到研究区实际并不存在荒漠化、石漠化问题, 故从水源涵养、水土保持、生物多样性保护三个方面来进行评价, 限于篇幅, 具体计算方法见蔡德峰^[25]与李月臣^[26]等的研究, 将计算值以分位数断点法进行分级, 结果如图 2所示。

2.2 生态脆弱性评价

生态脆弱性评价采用敏感性-压力度-恢复力模型(SRP模型)进行计算, 限于篇幅, 具体计算方法见贾晶晶^[27]的研究, 将计算结果按照分位数断点法进行分级, 从低到高分为一般脆弱、比较脆弱、中等脆弱、脆弱、极度脆弱五个等级, 计算结果如图 3。

2.3 限制转换区的提取与处理

将生态系统服务功能重要性与生态系统脆弱性计算结果进行空间叠加(一般脆弱/重要赋值为1, 比较脆弱/重要赋值为2, 以此类推), 计算并提取每个像元的最大值(值为1的像元对应Ⅰ级, 值为2的像元对应为Ⅱ级, 以此类推), 并在此基础上叠加研究区现阶段已划定的生态保护区数据(现阶段已划定保护区直接赋值为Ⅴ级), 得出生态综合评价结果(如图 4所示), 即等级越高的区域生态保护需求度越高。

最终计算结果显示, 研究区总体生态保护需求度较高, 其中处于Ⅲ级以上的区域约占研究区总面积的79.5%, 其中等级为Ⅴ级的区域占比为42.61%, 主要分布于西部高原高山区; 等级为Ⅰ级与Ⅱ级区域共占比为20.5%, 主要分布在东部平原地区。

根据计算结果, 将研究区生态综合评价等级为Ⅴ级的区域提取, 作为后文规划情景下的土地利用转移限制区; 将研究区生态综合评价等级为Ⅳ级与Ⅴ级的区域进行提取作为后文生态保护情景下的土地利用转移限制区(图 5)。

3 多情景土地利用模拟 3.1 土地利用模拟参数设置

本研究以2020年为现状基期年份对未来年份进行模拟, 在模拟训练时分别以2010年、2015年为基期, 5年为跨度对2015年与2020年土地利用情况进行模拟并以对应年份实际土地利用数据验证, 2010—2015年模拟结果Kappa系数为0.8891, 2015—2020年模拟结果Kappa系数为0.8359, 故本次模拟选择2010—2015年土地利用变化规律进行后续模拟。

3.1.1 转移概率矩阵与需求预测

根据2010—2015年土地利用转移矩阵, 计算研究区2010年至2015年土地利用转移概率矩阵, 再利用马尔科夫链以2015年为基期, 计算出2020年研究区各类土地利用需求数量(表 2)。

3.1.2 基于神经网络的适宜性概率计算

选取研究区2015年土地利用数据作为基期数据, 将归一化处理后的相匹配年份的9个土地利用驱动因子栅格数据导入模型, 设置神经网络训练采样方法为随机采样, 采样参数设置为20, 隐藏层数量设置为13, 计算得到土地利用适宜性概率图集, 其均方根误差值(RMSE)为0.29163, 表示训练结果可信。

3.1.3 转移成本矩阵及邻域权重参数设置

转移成本矩阵设置参考王芳莉的研究^[28]与2015—2020年实际土地利用转移量进行设置, 同时将建设用地与水域间设置为相互不可转换, 转移成本矩阵在此不作赘述。

土地利用邻域权重参数范围为0—1之间, 领域权重参数值越接近0, 表示该类用地越容易转化为其他地类^[29], 本次研究邻域参数采用3×3的邻域窗口并参照研究区2015—2020年实际土地利用转移量进行计算与设置(表 3)。

3.2 土地利用训练模拟结果及验证 3.2.1 ROC曲线验证

由受试者特征曲线(ROC曲线)验证结果(图 6)可知, 耕地、林地、草地、水域、建设用地、未利用地模拟曲线下方面积(AUC值)分别为0.873、0.725、0.790、0.790、0.906、0.881, 各类用地AUC值高于0.7, 说明此次模拟精度较高。

3.2.2 空间精度验证

计算得到Kappa系数值为0.8282, 此外, 计算出2020年研究区土地利用模拟结果整体精度为0.8834, 表示2020年土地利用模拟结果与2020年研究区实际土地利用数据一致性较高(图 7), 模拟结果可信。

3.3 土地利用多情景模拟

自然发展情景  是指在2020—2030年研究区各类土地利用变化不受政策等影响, 主要受自然地理条件以及社会经济因素影响, 其变化发展保持2010—2015年的发展变化规律不变;

规划约束情景  结合研究区实际, 参考研究区所属地的规划规程等相关政策文件, 根据相关指标严格控制耕地、林地、建设用地数量, 在自然发展情景的基础之上, 增设前文计算出的生态综合评价等级为Ⅴ级的区域作为限制转换区, 将耕地向建设用地转出概率降低30%, 减少转化的这部分耕地保持原地类不变, 林地、水域向建设用地转化概率分别降低20%, 减少转化的部分保持原地类不变, 草地向建设用地转化的概率降低10%, 其中减少转化的部分保持原用地类型不变;

生态保护情景  以生态保护优先, 严格限制具有生态功能的林草地向其他用地类型转出的情景, 其限制转换区在规划约束发展情景基础上进一步扩大, 将生态综合评价等级为Ⅳ级的区域也纳入限制转换区, 同时参考相关研究^[30—31], 在自然发展情景基础之上, 将耕地向建设用地转移概率降低30%, 并将此部分增加到耕地向林地的转出之上, 草地、林地向建设用地转变概率分别降低50%, 另50%保持原地类不变。

根据Markov链, 计算不同情景下各用地类型需求如表 4所示。

3.4 多情景模拟结果

根据以上参数设置, 分别计算模拟出在自然发展情景、规划约束情景、生态保护情景三种情景之下2030年研究区土地利用状况(表 5)与空间分布情况(图 8)。

在自然发展情景下, 无其他政策等因素干扰, 研究区2020—2030年各类用地变化趋势与2010—2015年保持一致, 经过单位换算, 得出自然发展情景下, 2020—2030年间, 研究区耕地共减少381.04km², 林地共减少197.54km², 草地共增加168.7km², 水域共减少33.33km², 建设用地共增加449.66km², 未利用地减少4.29km²。

规划约束情景下, 耕地、林地得到一定程度保护, 建设用地总规模得到限制, 在该情景下, 2030年模拟结果耕地减少了320.93km², 相对自然发展情景耕地减少的情况得以好转; 该情景下林地较2020年减少了194.72km², 与自然发展情景相比, 林地缩减数量更少; 规划约束情景下草地增加169.13km², 水域减少33.06km², 建设用地增加386.58km², 未利用地减少4.84km²。

在生态保护情景下, 林地、草地等具有生态价值的用地类型得到相应保护, 在此情景下, 耕地、林地、未利用地与2020年实际土地利用数据相比, 分别减少了381.04km²、190.48km²、5.46m², 其中林地栅格数量虽处于减少趋势, 但与自然发展情景及规划约束情景相比, 生态保护情景下的林地缩减量得到了控制; 生态保护情景下, 2030年研究区草地、建设用地均呈增加趋势, 分别增加的面积为231.02km²、381.43km², 其中草地增加栅格数明显大于其他两种情景, 建设用地扩张程度明显小于自然发展情景且小于规划约束情景。

由于研究区在地理位置以及社会经济结构上的特殊性, 在三种情景下, 草地面积均被压缩; 在规划约束情景下耕地减少得到较为明显的改善, 但由于西部耕地资源短缺, 人地矛盾加剧, 耕地保护力度还需加大; 在自然发展情景下, 建设用地增长未被管控, 建设用地扩张大大压缩了其他地类空间, 而在规划约束情景特别是生态保护情景下研究区建设用地增涨幅度较小, 较难满足研究区经济发展的需要。因此需要提出一种综合集成各种情景优点, 既能更好地保护耕地与具有生态功能的林草地, 又能更好地满足研究区经济发展需求的土地利用方案。

4 斑块权衡与多情景综合集成 4.1 待权衡斑块的定义与提取

在未来国土空间开发与土地利用规划过程中, 一定存在重点关注区域与可优化空间。将三种情景模拟结果进行叠加分析, 将某像元的用地类型设置为A, 其中自然情景下该像元用地类型为A1, 规划约束情景下该像元用地类型为A2, 生态保护情景下该像元用地类型为A3, 若自然情景下该像元用地类型与规划约束情景下用地类型相同, 则有A1=A2, 若自然情景下该像元用地类型与规划约束情景不同, 则有A1≠A2, 在此基础上将三种情景不同结果进行排列组合, 共五种可能情况(图 9), 根据不同的情况提取斑块, 其中A1=A2=A3时, 该斑块三种情景下用地类型相同, 不存在用地类型冲突, 将此类对象进行提取, 当A1、A2、A3中有任意二者不等时, 表示该斑块不同情景间用地类型存在冲突, 故判定其为待权衡斑块, 将此类斑块进行提取。在未来的国土空间规划与开发中, 应重点关注与慎重考虑这些存在差异的区域。

4.2 待权衡斑块的多层级权衡判定

参考耕地质量等别、城建用地标准等相关政策与规程文件, 结合研究区实际情况, 设置生态综合评价等级是否为Ⅴ级为第一级判定依据, 坡度是否大于25°为第二级判定依据, 研究区2020年土地利用现状作为第三级判定依据, 林线分布海拔^[32—33]作为第四级判定依据, 主体功能区划分作为第五级判定依据, 整理出待权衡斑块多层级权衡判定的基本逻辑如下：

生态综合评价为Ⅴ级的区域应以生态效益优先, 保护林草地, 限制林草地转出; 坡度＞25°的区域更低于建设用地适宜性过低, 应考虑还林还草; 当待权衡图斑的生态保护需求度、坡度等条件较为相似时, 通过现状用地类型、林线分布、主体功能区划分进行判定; 判定应综合考虑经济效益与生态效益, 其中生态效益以林草地优先, 经济效益以建设用地优先; 根据研究区实际情况, 研究区绝大部分水域坡度不超过25°。

根据以上判定逻辑, 将待权衡斑块转为矢量数据并新建生态保护需求度是否为Ⅴ级(1是, 0否)、坡度是否大于25°(1是, 0否)、2020年土地利用类型(1耕地, 2林地, 3草地, 4水域, 5建设用地, 6未利用地)、海拔是否大于3500m(1是, 0否)、主体功能区划分(1农产品主产区, 2重点开发区)等多个字段, 通过Python等编程手段实现自动权衡判定。

4.3 基于权衡判定的综合集成结果

通过权衡判定将待权衡斑块赋予相对应的土地利用类型后转为100m×100m的栅格, 与三种情景下均为耕地、林草地、水域及其他用地、建设用地的图层进行合并集成, 得到基于权衡判定结果的土地利用方案(图 10)。

为方便计算统计, 也为了更直观地对比三种情景与权衡判定集成方案的差异, 将单位换算为平方千米进行比较, 计算结果如表 6所示。

由上表可知权衡集成方案中, 耕地面积较自然发展情景模拟结果增加了133.16km², 较规划约束情景模拟结果增加了73.05km², 较生态保护情景模拟结果增加了133.16km², 由此可见在权衡集成土地利用方案之下, 研究区耕地得到了更好的保护; 林地较自然发展情景模拟结果增加了209.48km², 较规划约束情景增加了206.66km², 较生态保护情景增加了202.42km², 从林地面积上看, 权衡集成的土地利用方案林地面积明显大于其他三种情景, 林地得到了较好的保护; 权衡集成方案中建设用地面积为1944.83km², 较自然发展情景减少了46.99km², 较规划约束情景增加了16.09km², 较生态保护情景增加了21.24km², 从建设用地数量上看, 权衡判定方案的建设用地面积仅次于自然发展情景, 高于规划约束情景与生态保护情景。权衡集成土地利用方案下, 耕地、林地都得到了更好的保护, 同时建设用地也能较好的发展, 因此, 合理的权衡判定对未来土地利用模拟研究区起到了优化的作用, 能在未来的土地利用规划中更好地兼顾生态保护与经济发展的需要。

5 讨论

多情景土地利用模拟的意义在于研究不同情景不同导向下的土地利用结构发展趋势, 为未来的国土空间规划布局提供建议与参考。从空间分布上看, 多情景土地利用模拟结果基本符合新一轮国土空间规划布局：未来耕地主要集中分布于研究区东南部地势平坦的“粮油主产区”, 如中江县、三台县、盐亭县、梓潼县等区县; 未来建设用地及其增长主要集中分布于“成都平原经济区”及其周边地区, 如大邑县、什邡市、广汉市等区县; 林草地主要分布于西部高原山地广布的“川西北生态示范区”, 该地区作为高生态组分与落后发展交织区域, 在发展高原特色牧业的同时也要注重对林草地的保护。

龙门山过渡带作为多种地形单元、多类国土空间交界处, 各类用地空间需求存在多样性与更强的竞争性, 传统的自然发展情景、规划约束情景以及生态保护情景三种情景的用地情况对比之下各有优劣, 均较难满足龙门山过渡带的实际土地利用需求, 如何权衡生态环境保护、社会经济发展与粮食生产安全之间的关系显得尤为重要。故本文在多情景模拟基础上综合考虑了三种不同情景的优点, 并通过斑块权衡与集成的方式提出新的土地利用集成方案, 在此方案之下(表 6), 建设用地扩张规模既能保证经济发展需求, 又能得到合理控制, 耕地、林地、草地面积均大于传统多情景土地利用模拟结果, 表明在该方案下, 龙门山过渡带的耕地、林地、草地均能得到较好的保护。因此, 基于多情景土地利用结果提出的权衡集成方案可以更好地协调生态保护、农业生产、经济发展之间的关系, 更能适应未来土地利用的多方面多层次的需求。

从基于神经网络CA与Markov链的FLUS模型虽然能很好地对未来土地利用情况进行模拟, 但在FLUS模型模拟未来土地利用过程中, 各土地利用类型需求量是以转移概率或历史数据为基础来进行预测; 在对驱动因子的选取方面, 由于FLUS模型自身要求必须使用可量化可栅格化的因子, 限制了因子的选取。二者均可能会在一定程度上影响模拟精度。

6 结论

本研究选取龙门山过渡带这一典型山地-平原过渡区, 在生态服务功能重要性与生态脆弱性综合评价基础上划定了不同情景的限制转换区, 对2030年该区域土地利用情况进行了多情景模拟, 在此基础上进行了用地斑块权衡, 提出了基于多情景模拟结果的权衡集成土地利用方案, 主要结论如下：

FLUS模型在山地-平原过渡区域仍具有较强的适用性, 模拟结果Kappa系数为0.8282, 整体精度为0.8834, 可以较好地模拟该区域未来土地利用变化进行模拟, 未来该区域国土空间规划可以参考此模型模拟结果。

在三种不同情景下, 研究区耕地、林地、水域、未利用地总体处于缩减状态, 其中耕地减少尤为明显, 减少面积均在300km²以上, 耕地仅在规划约束情景下缩减较少, 得到了相对较好的保护; 草地、建设用地总体处于扩张状态, 在自然发展情景下建设用地增长程度较高, 在规划发展情景与生态保护情景增量分别为386.58km²与381.43km²。

在基于用地斑块权衡提出的多情景权衡集成土地利用方案下, 耕地、林地都得到了更好的保护, 同时建设用地也能较好的发展, 综合考虑并了耕地保护、生态保护、经济发展等多种需求。通过生态综合评价等级、坡度、海拔、主体功能区等多条件多层次的逻辑判定方式, 可以实现较好地兼顾耕地保护、生态保护、经济发展多种需求的用地斑块自动化权衡与判定。