1. 引言
20世纪50年代以来，全球城市空间不断扩大，对可持续发展构成严重威胁。为缓解城镇发展的相关问题，各国政府都提出很多对应的措施，如国外的城市增长边界（Urban Growth Boundary, UGB）、城市增长区域（Urban Growth Area, UGA）、城市服务区域（Urban Service Area, USA）和中国的城镇开发边界。
将主体功能区规划、土地利用规划、城乡规划等空间规划融合为统一的国土空间规划，实现"多规合一"，是党中央、国务院做出的重大部署，"三区三线"（农业空间、生态空间、城镇空间，永久基本农田、生态保护红线、城镇开发边界）划定是国土空间规划的重要基础工作，"三线"是国土空间规划的核心要素。
2. 研究区概况与数据来源
2.1 研究区概况
澄江市
重点生态功能区，面积984km²，2021年人口17.3万人，城镇化率48%，GDP153亿元，域内有著名的抚仙湖，水域占比达到18.61%，形成了"七山、二水、一分坝"自然地理格局。
陆良县
农产品主产区，面积1989km²，2021年人口59.1万人，城镇化率40%，GDP307亿元，域内有云南省第一平坝，粮食年总产量长期稳定在32万t以上，具有"滇东粮仓"的美誉。
昆明市主城区
城市化地区，面积3206km²，2021年人口538.9万人，城镇化率达96%，GDP5339亿元，占全省比重约20%，是云南省经济和政治中心。
图1 滇中城市群及研究区的位置
2.2 数据来源
土地利用数据
- 2009、2019、2020、2021年LULC来源于"二调"和"三调"数据
- 生态保护红线、永久基本农田和城镇开发边界数据来源于国土空间规划中"三区三线"划定成果
驱动因子数据
- 降水、温度、日照时间、人口、GDP和DEM分辨率为1000m的栅格数据
- 河流、道路、政府、铁路、高速公路、医院和学校的矢量数据
3. 研究方法
3.1 技术路线
研究基于2009、2019、2020和2021年4期土地利用数据，选取不同的驱动因子和限制因子，构建耦合马尔可夫链、PLUS模型和层次分析法的预测模型，基于2009、2019和2020年数据，按照主体功能区定位设置不同情景模拟2021年土地利用变化，进行精度检验后，模拟2035年研究区土地利用变化情况，通过形态学腐蚀膨胀算法划定不同情景下2035年城镇开发边界。
图2 城镇开发边界划定技术路线
3.2 PLUS模型
PLUS模型是一种基于栅格数据生成新像元来计算土地类型变化的模型，在FLUS模型的基础上进一步开发完善。可以更好地挖掘土地利用变化的影响因素，从而更好地模拟土地类型的变化，通常认为PLUS模型较FLUS模型具有更好的模拟性能。
模型公式：
Wi = TAi / ∑TAi
式中：Wi代表第i种土地利用类型的邻域权重；TAi代表第i种土地利用类型向其他土地利用类型转换的栅格数量；n代表土地利用类型的种类。
采用Fom（Figure of Merit）值对模拟结果进行精度检验，Fom值为正确预测的栅格数与研究区栅格总数的比值。一般认为Fom值在0.1~0.3之间证明模型模拟精度达到要求。
精度检验公式：
Fom = B / (A + B + C + D)
式中：A代表预测地类发生转变但实际未变化的栅格数量；B代表预测地类转变且变化正确的栅格数量；C代表预测地类转变但变化错误的栅格数量；D代表预测地类未发生转变但实际发生转变的栅格数量。
3.3 匀速递增马尔可夫链
马尔可夫链能够预测事件发生概率的过程，常用于无后效性特征的事件，即下一个步长的状态只由当前状态决定，与过去状态无关。本文将其用于土地利用变化预测。
基本公式：
XY+1(x1,...,xn) = XY × Pij
式中：XY代表某一期的地类状态；XY+1代表下一步长的地类状态；x1,...,xn代表不同土地类型的栅格数量；Pij代表矩阵中Y时刻处于状态i转移到Y+1时刻处于状态j的概率。
普通的马尔可夫链仅能进行等步长预测，若采用2009年和2021年（步长为12年）进行预测，只能得到2033年的数量关系，不能完成2035年国土空间规划末期目标，为此本文提出匀速递增的马尔可夫链。
改进公式：
X'(x1,...,xn) = (z/N)(XY+1-XY) + XY (z=1,2,...,N)
式中：X'代表在一个步长中某一年的地类状态；z代表距该步长初期的年份；N代表步长。
3.4 多情景设置
本文结合不同主体功能区设置了对应情景，进行土地利用变化模拟。不同情景下的转换成本矩阵如表1所示。
生态保护情景
突出生态保护，在该情景下，林地和水域将不转化为其他地类。
农业生产情景
突出耕地保护，在该情景下，耕地仅可以向城镇建设用地转换，其他地类均可向耕地转换。
城镇发展情景
突出城镇建设发展，在该情景下，城镇建设用地不向任何地类转换，同时除水域以外的地类均可以向城镇建设用地转换。
自然发展情景
该情景中土地利用变化继续延用现有趋势，在模拟过程中延续2009—2021年的参数设置。
表1 2035年不同情景下的地类面积(km²)
| 研究区 | 情景 | 耕地 | 园地 | 林地 | 草地 | 城镇建设用地 | 其他建设用地 | 水域 | 其他土地 |
|---|---|---|---|---|---|---|---|---|---|
| 澄江市 | 自然发展 | 153.74 | 21.12 | 353.29 | 29.03 | 12.24 | 38.77 | 146.87 | 1.39 |
| 生态保护 | 153.74 | 15.12 | 368.16 | 25.39 | 11.76 | 34.01 | 146.87 | 1.39 | |
| 陆良县 | 自然发展 | 816.80 | 107.46 | 732.84 | 49.54 | 40.39 | 180.48 | 40.21 | 21.06 |
| 农业生产 | 860.49 | 106.84 | 704.66 | 49.54 | 32.81 | 156.74 | 53.49 | 24.21 | |
| 昆明市主城区 | 自然发展 | 296.84 | 169.97 | 1660.94 | 93.13 | 429.05 | 303.87 | 239.80 | 12.17 |
| 城镇发展 | 303.43 | 114.42 | 1604.53 | 93.13 | 488.65 | 344.16 | 245.26 | 12.17 |
4. 结果及分析
4.1 城市建设用地发展驱动因素贡献分析
澄江市 重点生态功能区
主要驱动因素：社会因素（36.33%）
次要驱动因素：区位因素（14.88%）和经济因素（8.85%）
陆良县 农产品主产区
主要驱动因素：社会因素（52.24%）
次要驱动因素：区位因素（8.17%）
昆明主城区 城市化地区
主要驱动因素：自然因素（35.36%）
次要驱动因素：经济因素（9.21%）和区位因素（8.81%）
4.2 土地利用变化情景分析
澄江市
生态保护情景下：城镇建设用地较自然发展情景减少0.48km²，其他建设用地减少4.75km²，林地增加39.06km²，更符合重点生态功能区定位。
陆良县
农业生产情景下：耕地面积较自然发展增加43.68km²，城镇建设用地减少7.58km²，避免耕地占用，符合农产品主产区定位。
昆明主城区
城镇发展情景下：城镇建设用地增加59.60km²，其他建设用地增加30.66km²，更符合城市化地区发展需求。
4.3 城镇开发边界划定
澄江市
现状规模：7.06km²
优化建议：缩减1.68km²，调出占用连片耕地区域
陆良县
现状规模：26.06km²
优化建议：布局优化，避让连片耕地，缩减0.54km²
昆明主城区
现状规模：494.40km²
优化建议：规模合理，局部布局优化
图3 2035年昆明市主城区国土空间规划与主体功能区下的城镇开发边界
5. 结论与讨论
2035年为国土空间规划的目标年，选择2035年作为研究节点具有重要的现实意义。传统的马尔可夫模型是基于等步长进行预测，预测目标年和现有数据步长不一致时，难以完成预测。在国土空间规划中需按照主体功能区定位，优化城市功能布局和空间结构，以城镇开发边界着手，引导城镇空间集约增长。
主要结论：
- 不同类型主体功能地区城镇建设用地变化驱动因素差异显著，生态功能区和农产品主产区等的城镇建设用地变化驱动因素主要为人口数、到医院的距离等社会因素；城市化地区的主要驱动因素为高程、日照时间等自然因素。
- 提出的匀速递增马尔可夫链对2019年和2020年的各个地类数量关系的实际值与预测值进行整体精度对比，分别为96.3%和97.8%；对2020年和2021年的各个地类数量关系的实际值与预测值进行整体精度对比，分别为98.0%和97.7%，证明该方法具有较好的预测精度和较强的适用性。
- 作为重点生态功能区的澄江市，城镇开发边界规模可在现划定成果上进行适当扣减；作为农产品主产区的陆良县，城镇开发边界规模基本合理，需优化布局；作为城市化地区的昆明市主城区，城镇开发边界规模已充分考虑了城市的发展需求。
本研究对国土空间规划中城镇开发边界的划定提出了一种新思路和方法。在主体功能区战略要求下，每个地区功能定位不同，发展方向不同，应按照主体功能区划分进行城镇开发边界的差异化划定，此方法适用于国内大多数地区。
研究展望
由于基础数据时段的限制，所提出的匀速递增马尔可夫链进行非等步长时间节点的土地类数量预测仅应用在滇中城市群，其可靠性尚需进一步验证；同时引起土地利用变化差异的驱动因素也需进一步梳理与研究。