珠江三角洲土地利用变化模拟与景观格局脆弱性分析

张经度

（广东省城乡规划设计研究院有限责任公司，广东 广州 510220）

0 引言

景观格局指景观的空间结构特征，是景观组成单元的类型、数目及空间分布与配置，是景观异质性在空间上的综合表现［1］。一般利用有生态学意义的景观格局指数描述景观格局的空间配置特征，景观格局的时空动态变化研究能够揭示景观格局与生态过程的相互作用关系［2］。因此，基于景观格局的时空演化过程研究，不仅能反映区域生态环境变化状况，也能有效预测生态过程的特征，进而有效指引区域景观格局优化和管理［3］。

随着社会经济的迅速发展，人类对土地资源的开发不断加重，生态环境受到的压力和胁迫不断增加。这加剧了景观结构和组成的变化，由此产生的生态问题屡见不鲜［4］。在这一背景下开展景观格局脆弱性研究尤其必要。景观格局脆弱性是指景观格局在受到外界扰动时其结构、功能和特性容易发生改变的属性［5］，体现了景观系统对外界干扰的敏感程度及适应能力［6］。任志远等［7］和梁佳欣等［8］分别利用景观脆弱度指数对银川盆地和山东省南四湖湿地景观格局脆弱度进行时空分异特征研究，认为景观脆弱度指数能够直观反映研究区脆弱度状况，揭示人为活动与景观脆弱性的关系。

当前，学者对景观格局脆弱性的研究主要包括景观格局脆弱性时空分异研究［9-10］、景观格局脆弱性的时空分异成因研究［11-12］等方面，针对景观格局脆弱性的预测性研究还有所欠缺。曾永年等［13］和杨阳等［14］分别利用CLUE-S模型，以青海高原东部和福建省平潭县海坛岛为例进行土地利用变化多情景模拟，并利用景观生态风险指数进行区域生态风险评价，研究结果有助于决策者判断不同情境的优劣程度，从而支撑土地利用科学规划。然而，有研究指出CLUE-S模型存在一定的不足［15］。而FLUS模型作为一种新型土地利用变化模拟模型能够克服CLUE-S模型存在的缺陷，获得更好的土地利用模拟效果［16-17］。因此，有必要在利用FLUS模型模拟土地利用多情景变化的基础上，对模拟结果进行景观格局脆弱性评价。所以，笔者以珠江三角洲地区为例，使用系统动力学（System Dynamics，SD）模型，依据对研究区未来发展的不同预计构建3种不同的土地利用变化路径，在此基础上利用FLUS模型模拟2035年土地利用变化状况（以2020年为基期），再使用景观格局脆弱性指数分析研究区不同情景下面临的景观风险，从而为研究区土地利用规划提供决策依据。

1 研究区概况与数据来源

1.1 研究区概况

珠江三角洲地区位于广东省中南部（21°40′～111°20′～115°30′E），是我国重要的经济中心，当地海运便利、交通通达，总面积5.54万km2（见图1）。2020年，珠江三角洲地区人口为7 801.43万人，生产总值为89 523.93亿元，城市化率为86.28%。随着城镇化的快速发展，珠江三角洲地区土地利用出现了剧烈变化，耕地、水体面积明显减小，建设用地面积显著增加，区域生态环境受到的胁迫明显增加，可持续发展进程受到严峻挑战。

注：该图基于国家基础地理信息中心标准地图服务网站下载的审图号为GS（2019）4342号的标准地图制作，底图无修改，图3和图4同。

1.2 数据来源及处理过程

笔者所使用数据包括土地利用数据和土地利用驱动因子数据，其中驱动因子数据包括基础地理数据、自然环境数据、社会经济数据。土地利用数据来自全球30 m地表覆盖数据集（GlobeLand30），从中提取研究区2020年土地利用数据，并根据研究区实际情况，将土地利用数据划分为耕地、林地、草地、水体、建设用地和未利用地六类。参考以往研究［18］，选取13个影响土地利用变化的驱动因子（见表1），再根据先前研究结论将所有数据重采样为100 m［17］，同时将所有数据统一为UTM投影坐标系。表1中，共有13个土地利用变化驱动因子，余下为指标因子，在系统动力学模型中使用。

表1 研究所使用的数据

2 研究方法

2.1 FLUS模型

FLUS模型由基于SD模型的土地利用数量预测模块和基于元胞自动机（Cellular Automaton，CA）的地类分配模块两部分组成。首先，利用SD模型根据研究区社会经济发展、开发政策导向的不同预期推演研究区未来土地利用数量变化。其次，使用人工神经网络（Artificial Neural Network，ANN）基于土地利用变化驱动因子估计不同地类转化概率，接着将地类转化概率、元胞间相互作用、地类变化趋势等因素结合起来计算元胞转化总概率。最后，利用轮盘赌来确定元胞是否发生类型转化，从而实现土地利用变化模拟。FLUS模型总体结构详见张经度等和Liu等研究［17-18］。

2.2 SD模型

笔者使用Vensim软件构建研究区未来土地利用需求的SD模型（见图2）。模型包括社会经济和土地利用两个子系统，前者以人口增长率和地区生产总值增长率为主导因素，主要模拟人口增长对粮食需求的影响及地区生产总值增长对固定资产投资的影响；后者包括6种土地利用类型，重点模拟社会经济因素及地类间相互作用对土地利用的影响，另选用城镇化率、技术进步（以粮食单位面积产量衡量）及退耕还林政策作为辅助因素。

2.3 景观脆弱度指数

景观脆弱度指数（Landscape Vulnerability Index，LVI）是生态系统对外界干扰的敏感性及受到干扰后自我恢复能力的表征［19］，由景观敏感度指数（Landscape Sensitivity Index，LSI）和景观适应度指数（Landscape Adaptation Index，LAI）组成。其计算公式为

景观敏感度指数体现了景观受到外界干扰时自身的反映程度，取决于外界干扰因素的强弱及景观变化的方向。其计算公式为

式（2）中：n为景观类型数量；i为用地类型；Vi为景观易损度指数，指景观在外界干扰下的损害程度［6］，参考前人研究和研究区实际，设定耕地、林地、草地、水体、建设用地和未利用土地的景观类型易损度分别为0.14、0.23、0.23、0.04、0.04和0.32，实际计算时若某一单元格网内不包含某种地类，则其余地类的易损度值按比例扩大；Ui为景观干扰度指数，指景观受到外界干扰时，景观内部由单一、规则、均值和连续的整体向破碎、零散、异质和不连续的破碎斑块变化［5］，其计算公式为

式（3）中：FNi为景观破碎度指数，值越低表示景观破碎程度越低，反之则越高；FDi为分形维数，值越接近1说明景观形状越简单，反之则越复杂，自然度越强；DOi为景观类型优势度，该指数表示景观内部某种景观的优势程度，值越大则景观优势程度越大，值越小则景观优势程度越小。

景观适应度指数是景观在外界干扰下的适应和恢复能力，与景观的结构、功能、多样性和分布均匀程度密切相关，其计算公式为

式（4）中：PRD为斑块丰度密度指数，表示斑块内物种的丰度；SHDI为香农多样性指数，表示斑块内物种多样性；SHEI为香农均匀度指数，表示不同生态系统的分配均匀度，值越高表明景观格局越均匀，系统越稳定。

2.4 空间变异理论

空间变异理论研究对象为区域化变量，即在空间上分布的变量。通常，空间上分布的变量具有2个性质：①随机性，即在局部某一点，空间变量的取值是随机的；②结构性，即对整个区域而言，存在一个总体或平均的结构，相邻变量的取值具有相关关系［20］。由于土地利用数据为空间变量，针对土地利用数据中存在的随机性和结构性特征，可以利用半变异函数描述其性质［21］。笔者使用GS+7.0软件对研究区景观脆弱性分布进行半变异函数拟合。半变异函数公式为

式（5）中：γ（h）为间距为h的半方差，N（h）是以h为间距的所有观测点的成对数目，Z（xi）和Z（xi+h）为不同距离的观测值。

3 实证研究

3.1 2020—2035年土地利用情景模拟与分析

根据对研究区未来发展的不同预计，通过设置不同模型参数及约束条件，构建基准情景、经济扩张情景、协调发展情景，利用FLUS模型分别模拟各情景下2035年研究区土地利用分布格局。

基准情景表示研究区人口、经济发展与技术进步速度延续现有趋势，土地利用需求由2010—2020年土地利用变化量线性拟合得到。该情景体现了研究区土地利用无约束发展情况。系统动力学模型预测得到该情景下研究区2035年耕地、林地、草地、水体、建设用地和未利用土地的面积分别为8 026.41、24 524.79、2 470.72、4 341.90、13 031.51 km2和26.07 km2。

经济扩张情景表示研究区实行经济主导的发展政策，当地政府为经济发展提供充足的土地供应。该情景展现了在较快经济增速的背景下，满足建设用地扩张需求的研究区土地利用变化情况。研究预测得到该情景下研究区2035年耕地、林地、草地、水体、建设用地和未利用土地的面积分别为7 529.50、24 432.77、1 576.30、3 960.79、14 879.45 km2和42.60 km2。

《广东省国土空间规划（2020—2035年）》要求严格落实耕地和永久基本农田保护任务［22］，据此在协调发展情景中引入耕地面积反馈调节机制，使研究区耕地面积维持在基准水平之上，同时以《广东省国土空间规划（2020—2035年）》中珠江三角洲陆域生态屏障为土地利用变化限制区域。该情景体现了在保证研究区生态屏障完整及控制耕地面积前提下的研究区土地利用变化情况。研究预测得到该情景下研究区2035年耕地、林地、草地、水体、建设用地和未利用土地的面积分别为8 778.65、24 783.62、2 671.74、3 803.93、12 357.27 km2和26.19 km2。

基于不同情景预测得到的地类数量，模拟得到研究区3种情景下土地利用分布格局（见图3）。在3种情景下，地类转化的主导过程均为耕地向建设用地转移，其中经济扩张情景下建设用地扩张最为剧烈，侵占了城郊大片耕地。从各土地利用类型变化情况来看，3种情景中耕地均呈现减少趋势且减少的空间分布接近。在基准情景中，佛山市耕地面积减少近600 km2，是研究区耕地减少幅度最大的地级市，东莞市、广州市、中山市耕地也出现明显减少；在经济扩张情景中，佛山市耕地面积流失近700 km2，而广州市、东莞市、深圳市耕地面积大幅减少，减少幅度超过基准情景；在协调发展情景中，佛山市、广州市、东莞市依然是耕地面积减少幅度最大的城市，但减少幅度小于前两种情景。3种情景的林地变化幅度均较小，但变化的空间差异较大。3种情景建设用地均出现增长且空间分布类似，都呈现中心剧烈、外围缓和的增长趋势，但经济扩张情景下建设用地增长幅度明显高于其他两种情景，其中佛山市建设用地增长量达762 km2，其他两种情景建设用地增幅较接近，但协调发展情景下建设用地增幅略小。

3.2 研究区2020—2035年景观格局变化分析

3.2.1 空间结构变异性分析。为了更好地分析研究区景观格局脆弱性分布特征，对研究区基期及3种情景模拟结果进行空间结构变异性分析，得到最佳拟合模型及参数如表2所示。基期及3种情景模拟结果的最佳拟合模型均为指数模型，基期的变程值远高于情景模拟结果的变程值，表明基期数据中空间自相关效应的作用范围最大。3种情景中，经济扩张情景的块金值（C0）最大，块金效应最弱，表明该情景随机部分造成的空间异质性较多；基准情景的块金值最小，块金效应最强，说明该情景下结构性因素对空间异质性影响最大。基期及2035年模拟结果的决定系数及残差平方和均表现较好，表明所选拟合模型能正确描述研究区景观格局脆弱性的空间结构特征。

表2 研究区基期及2035年3种情景景观脆弱度分布的半变异函数的拟合参数

3.2.2 景观格局脆弱性的时空演化分析。基于3.2.1得到的最佳拟合模型和参数，分析得到研究区各情景下景观格局脆弱性分布图（见图4）。为了更好地比较2020年基期与2035年各情景模拟结果的差异，使用等间距法将研究区景观格局脆弱性分布图分为低脆弱度（0.013 5～0.235 6）、较低脆弱度（0.235 7～0.321 9）、中脆弱度（0.322 0～0.389 2）、较高脆弱度（0.389 3～0.452 6）和高脆弱度（0.452 7～0.502 9）5个等级。

先利用格网均值评估研究区景观格局脆弱性整体水平，计算得到2020年基期、2035年基准情景下、2035年经济扩张情景下、2035年协调发展情景下的景观格局脆弱度指数均值分别为0.332 5、0.356 1、0.394 6、0.349 7，可以看到相对2020年基期，3种情景下研究区脆弱性整体水平均有所上升。这说明随着社会经济的发展，研究区面临的生态风险将不断增加，而快速扩张情景面临的景观风险最大，协调发展情景面临的景观风险最小。

从空间分布来看，研究区2020年基期和2035年情景模拟结果的景观脆弱性分布由内向外呈现高、低、高的特征，高脆弱度区域主要分布在研究区东部、西北、南部及研究区中部的城市连绵区，低脆弱度区域则主要分布在研究中部建成区的外围。如图4（a）所示，2020年基期较高和高脆弱度区域主要分布在研究区东北的惠东县、博罗县，西北的广宁县、封开县及中部的广州市。图4（b）为基准情景下研究区景观格局脆弱性分布情况，与基期相比，研究区中部较高和高脆弱度区域面积明显增加，表明建设用地扩张将造成区域景观脆弱度上升，区域面临的景观风险将增加，研究区西北部分高脆弱度区域增幅明显。图4（c）显示了经济扩张情景下研究区的景观格局脆弱性分布情况，与基准情景相比，研究区中部建成区较高和高脆弱度分布范围明显更广，研究区南部的深圳市高脆弱度区域大量增加，研究区内低景观脆弱度区域明显减少。图4（d）展示了协调发展情景下研究区景观格局脆弱性分布状况，与其他两种情景相比，该情景研究区中部建成区较高和高景观脆弱度范围最小，全域较低和低景观脆弱度面积最大。

4 讨论与建议

笔者基于FLUS模型模拟了2035年珠江三角洲地区3种情景下土地利用变化状况，在此基础上利用景观脆弱度指数评价基期及各模拟情景下景观格局脆弱度分布状况，并借助空间统计分析方法揭示其分异规律，展现了在不同社会经济发展模式和政策的影响下研究区未来土地利用变化情况及对应的景观格局脆弱性变化特征。

在3种情景下，2020—2035年研究区耕地均呈减少趋势，建设用地均呈增长趋势，但不同情景间存在差异。在经济扩张情景下，研究区中部尤其是佛山市、广州市建设用地迅速增加，耕地、水体等较高生态价值地类被大量侵占，生态环境急速恶化。在基准情景下，建设用地扩张幅度仍然较大，无约束的建设用地扩张将导致研究区生态屏障被破坏，人居环境质量下降。在协调发展情景下，相应约束条件的设置使得耕地、水体得到较好保护，建设用地布局趋向合理，研究区土地利用的可持续性在3种模拟情景中最好。

此研究的基期及情景模拟结果均显示，研究区景观格局脆弱性从内向外呈现高、低、高的特征，高脆弱度区域主要分布在研究区东北、西北和中部建成区，低脆弱度区域主要分布在研究区中部建成区的外围地区。与基期相比，经济扩张情景下研究区中部高脆弱度区域面积明显增加，表明建设用地的蔓延式扩张将导致区域生态风险加重；基准情景下研究区中部高脆弱度区域面积增幅依然明显；协调发展情景下中部建成区脆弱性与基准情景类似，但研究区西北和东部脆弱性分布范围较小，同时该情景下全域景观格局脆弱性最低，表明该情景下研究区面临的生态风险最小，土地利用布局最合理。

该研究结果表明，3种情景都显示未来研究区中部城市连绵区将随着建设用地增加而扩张，研究区西北林地连绵区将持续存在。这表明这些区域将长期处于高景观脆弱性状态，针对因城市连绵出现的高景观脆弱性状态，未来应加强风险防御体系建设，降低自然灾害、人为活动可能造成的生态风险；针对林地连绵出现的高景观脆弱性状态，未来应增加土地利用的空间异质性，以降低区域景观脆弱性水平。

笔者在土地利用多情景模拟的基础上，使用景观脆弱性指数探究研究区未来生态风险的演变状况，实现了区域生态风险的表征和可视化。然而，尽管景观脆弱性指数能表征景观面临的生态风险，但是无法解释研究区生态系统演化的前因后果，因而对生态风险预测缺乏指向性，所以接下来应针对特定生态灾害过程展开模拟，从而加深对景观脆弱性指数生态内涵的理解。