土地的功能性体现了土地资源通过提供多样化的商品和服务, 满足人类的需求的内涵^[1]。2001年, 经济合作与发展组织(OECD)进一步发展了这一概念, 并将其分析框架扩展到各种土地利用的多功能性。中国在工业和经济发展方面处于过渡期, 土地功能相应地从生产主导条件转变为生产, 生活和生态等多种功能相互协调的协调状态^[2]。但是, 随着社会经济的快速发展, 有限的土地资源与人类对各种土地功能的需求日益增加之间的矛盾不断加深。这种矛盾具体表现在城镇化发展、农田保护、生态保护空间的竞争更加突出^[3]。因此, 为了实现生态-生产-生活功能的协同优化和统筹发展, 政府部门已出台了一系列规划政策。从《全国主体功能区规划》对国土空间开发格局的构建^[4]到《全国国土规划纲要(2016—2030年)》对“三区三线”的划定要求, “三生功能”已经成为中国土地利用功能规划和布局的基本框架, 以实现可持续的社会、经济和生态发展。

近年来, 学者们对不同尺度上的土地利用功能, 特别是“三生功能”的最佳分配模式进行了大量研究, 而城市群作为中国主体功能区战略中的重点, 已有专家学者对长江中游城市群、珠三角城市群等经济发展速度较快的东部城市群进行了深入的研究^[5], 但西部地区城市群的研究仍旧较少。兰西城市群作为维护国家生态安全的重要支撑、西北地区发展的重要增长极, 在土地利用功能布局优化研究上具有一定典型意义。在研究方法上数学模型作为土地功能优化的热点和前沿, 可分为数量结构优化和空间分布优化两类。首先, 定量结构优化模型主要解决了土地利用结构的比例关系和承载上限的问题, 如：多目标规划、Markov链、灰色预测、多准则规划^[6—7]等。但由于土地类型之间存在复杂的非线性关系, 因此很难用传统的数学方法构建优化模型, 因此本文选择灰色多目标规划模型(GMOP), 更好的反映土地各功能的动态过程演变特征。其次, 空间分布优化模型主要解决土地利用分类的空间配置问题, 以实现最优适宜性, 出现了如：FLUS模型、CA-Markov模型、CLUE-S模型^[8—9]等地理信息软件进行土地空间优化。但是, 这些自下而上的方法将各类土地类型进行合理分配, 依赖于自身的转换规则, 缺乏土地利用结构的定量优化。因此, 为了解决这个问题, 本文使用了PLUS模型, 不仅保留了自下而上的土地利用模拟的优点^[10], 而且适应惯性竞争和轮盘竞争机制, 进一步完善了挖掘空间转换规则。GMOP-PLUS模型作为当前较新的耦合模型, 较少应用于国内研究中^[11—13], 模型的有效性也需要进一步探索。此外, 目前模型主要侧重于模拟引起土地利用功能变化的经济和社会指标, 而往往忽视了土地利用变化后造成的生态问题。土地利用功能优化布局中考虑环境生态安全状况, 有利于及时调整不合理的、超出土地承载力土地利用方式^[14], 是不可忽略的重要环节。生态安全的核心在于源地识别, 基于“源-汇”理论的生态安全是目前应用最广泛的方法之一^[15]。其中, 形态学空间格局(morphological spatial pattern analysis, MSPA)由于可以对区域土地利用数据的空间格局进行识别和分割, 更客观的识别生态源地, 近年来被广泛应用于生态空间结构研究中^[16]。在此基础上, 本文构建了一个新的模型——MSPA-GMOP-PLUS耦合模型, 将通过新的耦合模型, 系统有机地将区域的生态环境与城市群经济发展结合起来, 为兰西城市群土地利用功能布局优化提供更可持续的规划方案。

基于此, 本文以兰西城市群作为典型地区, 展开土地利用功能布局优化研究。创新性的将MSPA融入到GMOP-PLUS模型中, 把维护区域生态安全贯穿空间布局优化过程和结果中, 以保证兰西城市群国土空间资源可持续利用。并利用不同模拟情景对比, 进一步探讨了城市群土地利用功能规划布局模式, 为优化城市群可持续发展提供基础。此外, 构建MSPA- GMOP-PLUS模型可以最大限度地规避生态系统脆弱地区, 探寻“人-地”和谐发展的新思路, 刻画土地功能协调下的城市群空间布局优化新模式。研究结果以期为兰西城市群生态环境保护、资源有效利用, 国土空间规划编制提供理论参考。

1 研究区概况

兰西城市群作为支持西部地区经济发展的重要城市群之一, 依据国土空间均衡发展战略, 横跨甘肃省和青海省两省。同时, 兰西城市群处于黄河上游流域地区, 毗邻青藏高原生态屏障区, 相较于其他城市群, 维护国家生态安全的责任也较大。2020年9月甘肃省和青海省共同签订了《深化甘青合作共同推动兰西城市群高质量协同发展框架协议》, 协议也指出, 生态是兰西城市群建设和发展的最大潜力。此外, 兰西地区气候干燥, 蒸发量远高于降水量, 夏秋两个季节降水密集, 日照充足, 昼夜温差大。自然本地条件较差导致区域内植被覆盖程度低下, 水土流失严重, 生态环境脆弱。

2 数据与研究方法 2.1 数据来源及预处理

本文采用的2010年、2015年、2020年土地利用类型数据为100 m×100 m的栅格数据, 兰西城市群矢量边界是由中国地级市行政区划边界图提取而来, 生态系统服务价值、植被覆盖度(NDVI)、土壤类型、人均GDP、人口密度、年均气温数据为精度1 km×1 km栅格数据, 为保证数据精度一致性, 6类数据均经过重采样后得到100 m×100 m。以上数据均来自于中国科学院资源环境数据中心(http://www.resdc.cn)。文中涉及的地形数据有高程、坡度, 高程采用分辨率100 m的DEM数据, 坡度由DEM数据派生得到, DEM数据来源自地理空间数据云平台(http://www.gscloud.cn)。交通、水系数据为矢量数据, 距离数据通过欧氏距离法计算得到, 数据来源于北京大学地理数据平台(https://geodata.pku.edu.cn)。

2.2 研究方法 2.2.1 扩张“源”选取

结合兰西城市群的实际情况, 已有土地利用类型在本研究中被分为6类, 分别是林地、草地、耕地、水域、建筑用地、未利用地。在生态源地识别中将林地、水域作为MSPA剖析的前景, 其余地类为后景；在城镇建设源中将城镇用地设置为前景, 其余为后景。连通性连接指数(PC)被广泛应用于景观连通性的评价, 该指数可有效判断各个景观斑块之间的连通性强弱^[17]。本文基于Conefor软件将生态、建设斑块连通距离阈值设置为900、1000, 连通概率设为0.5, 对核心区进行重要源地识别, 具体计算公式如下：

(1)

式中：n为斑块数；a_i、a_j分别为斑块i、j的面积；A_L为景观总面积；p_ij^*是物种在斑块i、j中分散的最大概率。PC值在[0, 1], 值越大, 区域景观的连通性越高。

2.2.2 空间阻力面构建

最小累积阻力模型(MCR模型)是指物种从某个“源”点到目标所在地的过程中所需克服的总阻力的模型^[18]。具体公式如下：

(2)

式中, f表示MCR和迁移过程的正相关关系。D_ij是源地j到i类景观的距离；R_i是景观单元i在某种运动过程中的阻力系数。

对MCR模型构建起到的重要影响是其阻力因子的选择。人口发展状况和社会经济状况的改善导致研究区间邻里效应和空间互动的增强, 有利于土地空间变化的空间溢出, 其扩张特征不同于单一城市。这意味着, 优化和规范城市群空间, 不仅要考虑资源环境、空间区位等直接因素, 还要考虑城市群对土地空间扩张的人为作用。此外, 由于本文源地不仅包含生态源地, 还考虑了生活源地, 因此本文参考前人的经验^[19]及研究区实际情况, 选取11个指标作为生态源、生活源阻力因子(表 1)。同时对各阻力因子等级进行划分并采用层次分析法确定权重, 每个因子划分为5个等级：1、2、3、4、5, 等级越低, 阻力越小, 表示生态源地越容易得到保护和维持, 反之则表示生活空间扩张能力越强(EE为生态扩张, LE为生活扩张)。各类阻力因子的分级标准参考《资源环境承载能力和国土空间开发适宜性评价指南(试行)》, 指南中未提及的分类方法在参文献^[20]的基础上, 结合《兰州-西宁城市群发展规划》进行调整。

在上述分析的基础上, 以生态用地扩张和生活用地扩张的最小累积阻力差值为基础, 进而评价国土空间中适宜生态保护及城市发展。具体公式为：

(3)

当MCR_差值＜0时, 适宜作为生态保护区；当MCR_差值>0时, 适宜作为城镇建设发展区；MCR_差值=0时, 为两者的分界线。

2.2.3 优化配置情景模拟

GMOP模型能够充分考虑未来土地利用中目标函数的不确定性, 并利用多重约束解决土地利用类型优化的多目标问题^[20]。本文决策变量的设置在考虑实际情况以及数据的可操作性后, 选取耕地、林地、草地、建设用地、水域和未利用地6个决策变量。本文以生态系统服务价值来表征土地生态效益, 通过参考文献^[21]及中国陆地生态系统生态服务价值当量因子表, 得到适用兰西城市群的生态系统服务价值当量因子表, 再结合GM(1, 1)模型预测的2030年兰西城市群面积粮食经济产值, 计算得到生态系统服务价值系数表(表 2), 单位为万元/km²。

基于此, 研究区各决策变量的生态效益目标函数为：

(4)

式中, X₁为耕地, X₂为林地, X₃为草地, X₄为水域, X₅为未利用地, X₆为建设用地。此外, 约束条件设置为：

(1) 总面积约束：

根据土地面积总量平衡原理, 优化后流域内各土地利用类型总面积应前后保持不变, 即S=X₁+X₂+X₃+X₄+X₅+X₆。

(2) 耕地需求约束：

根据《兰西城市群建设农业农村领域合作发展框架协议》, 到2030年基本农田面积不得低于15000 km², 基本农田保护率不得低于81.37%。即: 81.37%×X₁≥15000。

(3) 林地、草地需求约束：

由于城市群自然本底脆弱性,地表植被受人类活动破坏后极易出现水土流失、荒漠化等现象, 因此, 从林地面积变化趋势上来看, 应以目前2020年为底线, 即X₂≥8910.16 km², X₃≥58624.02 km²。

(4) 水域、未利用地需求约束：

随着城镇化水平的不断提高, 水域、未利用地的开发越来越强, 故2030年未利用地面积会小于现状面积, 而水域面积在政府约束下有所好转。同时参考《兰州-西宁城市群发展规划》得到可开发的水域与未利用地面积。则可得：1523.80 km²≤X₄, 7484.77 km²≥X₆≥6364.91 km²。

(5) 建设用地需求约束：

在满足生态效益最大化的前提下, 为保障研究区经济社会发展, 到2030年城市群建设用地面积不得低于2020年水平, 以满足城市化需要, 即X₅≥2065.94 km²。

PLUS模型利用随机森林分类探索土地利用的未来变化趋势, 在考虑研究区自然、社会两方面后, 结合土地利用数据及驱动因子(表 1), 计算了各类型LULC的适宜性概率。并为保证结果科学性, 利用2000年和2010年模拟2020年, 得到Kappa系数为0.79, 具有较高的一致性且满足研究需要。

2.2.4 土地利用功能重构原则

基于土地利用的多功能理论^[22], 可以将土地利用类型分为农业, 城镇或生态用地的单一或复合功能, 然后将土地利用地块或网格单元合并形成不同功能导向的空间是目前应用最广泛的空间重构方式^[23]。基于此, 本文利用已有文献对“三生”功能划分标准^[24], 根据兰西城市群的农业生产、经济发展、人口增长和生态服务等需求, 基于土地利用功能将土地空间分为三组。分别是生活空间, 生产空间和生态空间^[25], 生活空间是指人类生活生存和发展过程中提供各种空间承载的功能, 如居住、交通、公共服务等功能, 是维持城市群运行的基地, 主要包括城市、农村居住用地。生产空间是以土地为载体进行社会生产活动来保证人类生存及发展的基础性功能, 包含了农业、工业生产的耕地和工矿用地。生态空间主要是包含气候调节、水文调节等维持人类生活的自然条件, 是生态系统中的基本功能, 包括林地, 草地和水域。为了更准确的描述土地利用功能, 本文通过构建3 km×3 km的蜂窝六边形统计空间, 统计每个栅格单元各空间的比例^[26]。以三元图各坐标的0.25为参考线, 将三元图划分为7个区, 区分不同功能的组合类型(图 2)。

基于三元图将土地利用功能结构重构为生产空间、生活空间、生态空间、生产-生活空间、生态-生产空间、生态-生活空间和生态-生活-生产均衡空间七种类型。各类空间的比例分布如表 3所示, 表中ABC代表生产空间、生活空间、生态空间中的任意一个类型, 通过不同组合和土地利用功能占比得到重构结构。

3 结果与分析 3.1 土地利用功能分区 3.1.1 生态安全格局分析

由图 3所示, 基于MSPA分析的兰西城市群建设景观要素面积为1689.06 km², 其中建设核心区面积为547.17 km², 占研究区景观要素面积的32.39%。桥接区是连接核心区与核心区重要区域, 其面积402.12 km², 占研究区景观要素面积的23.80%。兰西城市群的生态景观分类(图 3)的总面积为111834.1 km²。其中, 核心区面积为43355.5 km², 占区域内景观要素的38.74%, 面积较大, 稳定性较强。边缘区面积为31345.23 km², 占景观要素的28.03%, 边缘区是核心区和非绿色景观区域之间的过渡区域, 边缘区面积越多, 适宜物种分散越严重。

对所有核心区斑块作为潜在源地计算各斑块的连通性重要程度并对其进行排序, 最后选取dPC>1的斑块分别作为生活源地和生态源地。兰西城市群生活源地总面积为457.44 km², 共10个源地, 占城市群总面积的0.05%。生态源地总面积为13751.92 km², 占区域面积的1.41%, 共20个生态源地。最终生态、建设源地如图 4所示, 生态源地主要沿祁连山、甘南高原、三江源生态安全屏障带分布, 生活源地多分布于研究区中部河谷地区, 符合兰西城市群实际发展情况。

3.1.2 土地功能区划分结果

根据表 1生态源、生活源的阻力因子分级情况, 重分类后得到各阻力因子空间分布结果, 再对该结果进行加权求和后根据值域的分为5个等级, 得到不同源地的综合阻力成本, 以此构建综合阻力成本。

将生态源和建设源2个源地的最小累积阻力面之差进行分区。分区阈值由MCR差值与面积曲线的突变点确定。MCR差值分别在-256076.973, -148874.097, 2004.025, 45679.271时发生突变, 基于此划定5种国土空间功能。

其中, 生态保护区的面积为29209.24 km², 占兰西城市群面积的30.7%, 分布城市群西的北部, 呈“C”形, 构成了城市群主要的生态屏障区(图 5)。这些区域地势较高, 林地、河流较密集, 作为研究区各类自然保护区的核心地带, 中部及东部分布相对较少。生态边缘区的面积为19016.12 km², 占研究区面积的19.2%(图 5), 以西北部、西南部分布较多农业开发区21884.32 km², 主要以耕地、基本农田、园地和高原农业组成, 占研究区总面积的22.7%, 为城市群农业经济发展的重点区域。农业边缘区面积最小为4240.2 km², 占比为4.4%, 可以作为生活空间扩张用地或是工矿生产空间开发的后备用地。城镇开发区作为今后生活空间集中扩张区域, 面积为21205.92 km², 占兰西城市群总面积的22.2%。主要沿城市群经济发展较快的兰州、西宁、白银等地势较平坦的地区扩张发展, 具有一定开发基础。

3.2 土地利用功能模拟与优化 3.2.1 远景年模拟预测

为探究2030年国土空间发展状况, 以2020年为基期年, 分别利用MSPA-GMOP-PLUS模型和GMOP-PLUS模型模拟2030年兰西城市群不同约束下土地利用状况, 两种模型下土地利用转移成本矩阵值(土地利用转移概率)设定相同, 两种2030年兰西城市群“生活-生产-生态空间”分布模式(图 6)。GMOP-PLUS模型在遵循土地转换规则之上, 重点考虑各阻力因子对城市群发展的限制, 符合城市群自然发展状态。MSPA-GMOP-PLUS模型是在自然发展的前提下, 在模型中融入了MSPA的空间功能分区结果, 将生态保护区、农业开发区为禁止开发区, 生态边缘区、农业边缘区为限制开发区, 最大程度的维护区域生态安全和保护区域耕地资源。

GMOP-PLUS模型下2030年城市群生活空间、生态空间的面积都明显增长(图 6)。其中, 生活空间的面积达到2749.17 km², 是2020年的0.33倍, 主要为农村居民点的蔓延增长。生态空间中的草地面积增长素的较慢, 增长面积为0.3 km²。生产空间中的耕地面积和生态空间中的未利用地面积均有所减少, 未利用地面积幅度最大, 超过1155.55 km²。总的来说, GMOP-PLUS模型下城市群生活空间扩张对生产、生态用地存在一定影响, 其中, 生态用地受生活空间的挤压现象更为严重。MSPA-GMOP-PLUS模型下的土地利用变化是综合生态保护和生活约束的结果(图 6), 是城市群可持续发展的基础。该模式下, 生产空间特别是耕地面积的下降幅度明显缩小, 较2020年减少1.18%, 比未考虑MSPA的模式提高0.12%。生活空间中的林地、草地面积相较于2020年涨幅较大, 增至59826.90 km²、952534 km², 涨幅均超过2%, 水域面积仅增长了154.37 km², 涨幅较小。整体而言, 生活空间的面积为2173.16 km², 生态空间的面积为169231.92 km², 分别比GMOP-PLUS模式下少576.01 km², 高1225.11 km²。土地利用结构变化的实质是自然生态、生产和生活对土地资源的竞争, 而可持续发展实际上是将三者进行平衡, 不仅要对生活空间的发展留有余量, 还要考虑粮食质量安全, 维持生态安全。MSPA-GMOP-PLUS模式在控制城镇空间的同时兼顾了生态保护发展, 模拟结果与《兰州-西宁城市群发展规划》的发展目标一致。

3.2.2 土地利用功能布局优化

为了将城市群内部的土地结构、功能进行进一步细化, 基于三原图得到最终土地利用功能优化布局图(图 7)。其中, P区为生产空间, L区为生活空间, E区为生态空间, PE区为生产-生态, PL区为生活-生产, LE区为生活-生态, PLE为三生空间均衡型。生产空间、生态空间、生活空间、生态-生产空间、生活-生产空间、生产-生态空间和生产-生活-生态空间7大类, 其面积分别为4134.10、61948.63、467.98、26431.77、1094.71、1005.41、972.34 km²。

2030年兰西城市群生态空间面积为最大, 且主要分布于城市群西部地区, 呈连片发展。生态空间的范围与生态保护区的范围几乎一致, 不仅包括有国家级和省级重点风景名胜区、自然保护区、国家森林公园等生态保护区, 还包含龙羊峡、刘家峡、黄河、湟水河等重要水资源, 生态敏感性较高, 是维持兰西城市群整体生态安全的重要区域。生态-生产面积排第二位, 与生态、生产空间紧密相连, 横跨生态边缘区和农业开发区, 主要在城市群东部地区。这类地区不仅承担着黄河流域生态保护、国家级自然保护区生态管理等多个职能, 还是农业生产发展的重要力量。生产空间广泛分布于研究区内各地区, 面积较大但分布零散。主要原因是城市群单一的生产空间主要包括的是农业生产用地, 研究区内海拔落差大, 农业用地难以集中发展。生态空间、生态-生产-生活空间发展空间较小, 且均位于地势较平缓的地区。2030年兰西城市群的生活空间逐渐形成节约集约用地的趋势, 并且依旧为城市群的经济核心地区——兰州市、西宁市主城区及白银区、临夏市等发展基础较好的城市。生态-生活-生产空间呈零星发展, 且多处于几类空间交接的位置, 经济基础较好, 可以在保障生态功能后依托地方特色资源, 发展农畜产品精深加工、商贸物流, 扩大居民生活居住区或活动范围。

4 讨论与结论 4.1 讨论

土地利用功能布局优化是一个复杂的空间和社会经济决策过程^[27]。MSPA-GMOP-PLUS模型基于MSPA先对“三生功能”进行划分, 再对土地利用功能进行优化布局。通过“点-面”的模式对生态、生活空间结构要素的优化重组, 有效的防止城市用地过度增长、自然生态空间占用、生境质量下降等现象。

MSPA-GMOP-PLUS模型在土地利用功能优化布局中考虑了一个关键问题：如何建立最优空间结构^[28]。已有研究证明, 土地利用空间中存在一些重要的斑块, 这些斑块及其内部结构不仅对区域整体发展具有重要意义^[29], 还对土地利用功能合理配置产生间接影响^[30]。由于源地作为区域内有重要意义的斑块, 源地内部结构往往通过连通度表现, 基于此, 本文利用连通度指数对源地进行检验, 以保证土地利用空间结构最优。随斑块粒度的增加生态源地、建设源地连接度指数虽然呈现持续下降趋势, 但是分别在900、1000 m之后迅速下降, 根据质变点的原则^[31], 本文将1000、900 m分别作为建设源地、生态源地最佳粒度范围。基于此得到的土地利用空间结构与传统模拟(GMOP-PLUS)结果相差较大, 特别是对生活空间中农村居民点扩散的抑制作用格外明显(图 6)。农村居民点普遍分布不均, 斑块小且连通程度低, 在MSPA-GMOP-PLUS模拟结果中, 2020—2030年城市群内的农村居民点仅由500.72 km²上升到676.19 km², 面积增加了35.04%, 而传统模拟结果中农村居民点面积增加了40.11%。由此可见, 新模型不仅满足了城市群生活空间集约发展的需求, 还确保了生态空间优先发展的前提。

由于研究区内主要以传统种植农业为主, 对自然条件依赖度高导致生产空间较为分散, MSPA未能为直接对生产空间发展提供有效策略。但是, 考虑到MSPA-GMOP-PLUS模型能通过对生活、生态空间优化分配的角度提升生产空间, 所以认为模型结果依旧有效, 且能为兰西城市群土地利用功能布局优化提供新的思路。

4.2 结论

为了探索城市群最优的土地利用功能布局模式, 本文基于生态安全的角度提出了MSPA-GMOP-PLUS耦合模型, 模拟了2030年城市群布局优化情况。为了优化土地功能以实现可持续发展, 新的模型虽然着重考虑了生态核心区, 但对潜在的城市开发空间也进行了整合。同时, 为了验证模型生成的空间模式的科学性和实用性, 将其与传统的GMOP-PLUS模型进行了比较。对比分析发现, 加入MSPA的土地利用情景更符合促进城市群未来发展。总而言之, 将生态安全纳入土地利用功能优化中对于促进兰西城市群均衡发展具有一定意义, 是国土空间规划的重要规划工具。研究结果表明：

(1) 兰西城市群生态源地总规模较大, 总面积为13751.92 km², 景观组成以林地、水域为主, 但空间分布不均, 东部地区斑块分布较少。建设源地总面积为457.44 km², 以经济发展速度较快的兰州、西宁、白银为主。在此基础上得到2020年兰西城市群国土空间功能区。其中生态保护区面积最大, 为29209.24 km²；农业开发区面积为21884.32 km², 位居第二；城镇开发区面积为21205.92 km²位列第三。

(3) 基于MSPA-GMOP-PLUS耦合模型得到兰西城市群2030年模拟结果与政府规划相一致, 且更适用于研究区未来发展。以此划分出7种功能类型：生活空间、生产空间、生态空间、生活-生产空间、生活-生态空间、生产-生态空间、生态-生活-生产空间, 各类空间面积存在较大差异, 其中生态-生产空间面积最大为35413.91 km², 生活空间面积最小为3088.19 km²。