基于多模型的黑河中游甘临高地区土地利用情景模拟

蒋小芳 段翰晨 廖 杰 宋 翔 薛 娴

(1.中国科学院西北生态环境资源研究院沙漠与沙漠化重点实验室， 兰州 730000； 2.中国科学院西北生态环境资源研究院干旱区盐渍化研究站， 兰州 730000)

0 引言

土地利用变化关系到生态系统的稳定和持续发展，是土地利用/覆被变化(Land use and cover change，LUCC)研究的热点问题[1]。土地利用模型为这一研究提供了强有力的技术支撑。CA模型出现于20世纪40年代，TOBLER[2]于1979年首次将CA模型应用于地理建模。CA模型涉及人类系统与城市土地系统之间的相互作用，具有同质性、空间离散性和时空局部性等特点，空间模拟能力优越[3-5]。然而传统CA模型中每个像元内的土地利用类型单一，而实际的土地利用类型分布情况复杂。MCCA模型利用的混合元胞包括了不同土地利用类型，符合现实中复杂的土地利用结构，在元胞空间、元胞状态和转化规则的元胞邻域方面均与传统的CA模型存在差异，实现了从定性、静态模拟到定量、动态模拟的跨越[6]。

然而单一模型难以兼具数量模拟和空间格局预测的功能[7]。随着研究的深入，综合数量模型与空间模型的耦合模型出现[8-14]。目前已有的数量模拟模型主要有马尔可夫链(Markov chain，MC)、灰色预测模型(Grey model，GM)和MOP模型等[5,15-20]。MOP模型可用于计算土地系统中多维目标函数，具有强大的计算功能，是一种自上而下的方法，但在空间优化方面存在缺陷[21]。CA和MCCA模型均为自下而上的空间格局模拟方法，能够生成优化的土地利用结构[22]。因此有必要同时采用自顶向下和自底向上的耦合模型来优化土地利用数量-空间格局[23]。

部分学者采用包括小尺度土地利用变化及其空间效益模型(Conversion of land use and its effects at small regional extent，CLUE-S)在内的纯净元胞研究模型对黑河中游甘州区或张掖市等区域的土地利用结构进行模拟[24-25]，但是缺乏对混合元胞的探究。鉴于此，本文以黑河中游甘临高地区为例，确定MCCA模型在研究区的适用性后，分别将MOP和普通线性回归这两个数量预测模型与MCCA空间模拟模型耦合，模拟得到2035年可持续发展(Sustainable development，SUD)情景和基本发展(Basic development，BAD)情景下土地利用空间结构，采用效益定量评价方法对比分析，优选最佳的土地配置方案。

1 研究区和研究方法

1.1 研究区域概况

黑河流域地处中国西北干旱气候区。本文研究范围为黑河流域中游的张掖市甘州区、临泽县和高台县(甘临高地区，图1)。该区域年均气温约为8.2℃，年均降水量在80～150 mm之间，生态环境较为脆弱。土壤类型以水稻土、风沙土、草甸土、石质土和棕钙土为主。土地利用类型较多，主要包括耕地、草地、建设用地和荒漠等。黑河流域中游经济发展速度较快，人口较为稠密，人地矛盾突出。研究该区土地利用变化对于西北地区的经济发展和生态环境保护具有重要意义。

图1 研究区简图Fig.1 Chart of study area

1.2 土地利用研究模型

1.2.1MCCA模型

(1)MCCA模型原理

MCCA模型主要包括3部分[6](图2)。第1部分是混合元胞内多维土地利用变量连续变化的机理挖掘。MCCA模型基于起始年份和目标年份两期土地利用数据得出不同地类的连续变化，采用随机森林回归(Random forest regression，RFR)算法挖掘不同地类与驱动力因素之间的映射关系。RFR有利于克服不同空间变量的多重共线性问题。RFR首先训练各地类的变化量(即增加量)与驱动因素之间的关系，因变量为每个地类的变化量，自变量为各驱动因素；然后预测得到各地类中每个元胞的变化潜力。与以前的模型不同，RFR并未同时拟合多个地类，而是分别训练每种地类的变化潜力；而某一地类的增加表示其他地类的减少，因此模型只关注增加的地类，减少表示不增加，避免重复计算。

图2 技术流程图Fig.2 Flow chart of this research

第2部分是耦合亚元胞尺度土地竞争与定量转化机制的CA模型。传统的CA模型基于纯净元胞进行模拟，每个元胞代表一个地类，这脱离了实际情况，现实中每个元胞内部含有多个地类；而MCCA模型基于混合元胞内部每种地类的占比进行模拟，即亚元胞。该模型通过驱动系数调节不同地类的变化方向，首先需要计算土地利用成分的总变化概率，计算公式为

(1)

式中t——迭代次数i——混合元胞

k——地类组分O——总变化率

P——变化潜力Ω——领域效应

D——需求反馈

第3部分是MCCA模型模拟结果评价体系，主要包括亚像元混淆矩阵(Sub-pixel confusion matrix，SCM)、平均相对熵(Relative entropy，RE)和混合元胞质量系数(Mixed-cell figure of merit，mcFoM)，分别表征模拟结果的整体精度、土地利用结构相似度和混合元胞的变化精度。SCM越低说明模拟结果的整体精度越高，RE越接近0说明模拟过程中信息丢失越少，mcFoM越接近1说明模拟结果越准确。除此之外还有总体分类精度(Overall accuracy，OA)、生产者精度(Producer’s accuracy，PA)和用户精度(User’s accuracy，UA)，均是越接近1模拟精度越高。

MCCA模型采样率设置为10%，回归树数目为100，领域为3；转换矩阵中0表示两种土地利用类型之间无法转换，1表示可以转换，本研究在模拟过程中未设置限制转换区域；步长表示土地利用类型的转化率，经多次测试发现步长均为1时各项精度预测指标值最大。

(2)MCCA模型数据源及预处理

2000年和2015年的土地利用数据来源于中国科学院地理科学与资源研究所的资源环境科学与数据中心(http:∥www.resdc.cn)，为纯像元数据。该数据的每个像元表征一种地类，空间分辨率为30 m(表1)。通过ArcGIS 10.2软件的分区统计制表工具将这一数据聚合成空间分辨率为250 m的混合像元数据，并计算得到每个像元内部耕地、林地、草地、水域、建设用地和未利用地的覆盖比例数据。

表1 MCCA模型数据源信息Tab.1 Research data information of MCCA model

其他表征土地利用变化的驱动力数据主要有自然环境数据(高程、坡度、坡向、河流矢量数据、2015年的年均气温及年均降水量数据)和社会经济驱动力数据(道路矢量数据、居民点矢量数据、2015年的人口和国内生产总值(Gross domestic product，GDP))，与土地利用数据的来源相同。ArcGIS 10.2软件中的欧氏距离计算工具可处理得到每个像元与道路、河流、居民点的距离栅格数据，该数据的空间分辨率为30 m。MCCA模型提供类似于ArcGIS软件重采样的坐标对齐机制，不要求驱动力栅格数据具有统一的像元分辨率(见MCCA模型操作手册，下载网址：https: ∥github.com/HPSCIL/Mixed_Cell_Cellullar_Automata)。考虑到大部分驱动力数据的空间分辨率为30 m，为提高MCCA模型运行效率，通过ArcGIS 10.2软件中的重采样工具将驱动力数据的空间分辨率统一为30 m。

1.2.2普通线性回归模型

MCCA模型主要用于土地利用空间结构可视化模拟，该模型开发者研发的MCCA软件同时具有普通线性回归功能(图2中简称为MCCA的线性预测模块)，该功能基于一次函数进行预测[6]。除MCCA模型空间模拟需要的2000年和2015年的土地利用数据外，本文引入2005年和2010年各土地利用类型的面积数据对未来各地类面积进行线性回归预测，其来源和处理方法与2000年和2015年土地利用数据一致。

1.2.3MOP模型

MOP模型包括目标函数和约束条件两部分，其研究结果对于认识未来的土地利用结构和当前土地规划具有参考价值[26]。研究区面积为1 073 078 hm2，本研究利用MOP算法，基于Lingo 12.0软件，对2035年各地类面积进行了优化。

(1) 目标函数

土地利用结构的优化应该遵循动态性、统筹兼顾和生态优先原则。土地系统内部的各个要素处在动态变化过程中。土地利用结构优化时应基于系统论观点，兼顾整体和局部、协调人地矛盾，保证人与自然和谐发展。

生态效益目标函数的确定能提高研究区的生态价值，增强生态服务功能，维持生态系统的稳定性。当前生态效益主要采用价值量估算，并用货币值量化。本研究采用谢高地等[27]的中国生态系统服务价值当量因子表来估算生态效益，采用黑河中游张掖地区的粮食单位面积产量与全国粮食单位面积产量的比值作为修订系数。参考《中国统计年鉴》和《张掖年鉴》[28-29]，2015年全国粮食单位面积的产量为5 482.900 kg/hm2，全国粮食平均价格为2.370元/kg，张掖市同年的粮食单产为7 198.260 kg/hm2，计算得到修订系数为1.313。单位面积耕地生态系统服务价值当量因子的经济价值为本年度全国粮食单产价值的1/7，其计算结果为1 856元/hm2。校正后耕地、林地、草地、水域和未利用地的生态效益分别为16 840、53 250、17 650、112 040、1 020元/hm2，建设用地的生态效益较低，本研究按零值计算。

黑河中游甘临高地区土地利用结构优化是为获得最大经济效益，经济效益目标的实现有利于提升工农业发展水平。本文采用每种地类单位面积的国内生产总值作为经济效益目标函数的系数。参考2015年相关统计年鉴的经济指标数据[28-31]，耕地、林地、草地、水域和建设用地的经济效益核算结果分别为30 400、400、1 700、200、379 000元/hm2。

社会效益主要包括粮食安全保障和农民生存保障，参考前人的研究，耕地、林地、草地、水域、建设用地和未利用地的社会价值分别设置为23 360、2 520、2 340、5 160、99 040、10元/hm2[32]。

(2) 约束条件

约束条件包括3部分：土地总量约束、社会经济约束和生态环境约束。研究区各种土地利用类型的约束条件主要参考2015年实际土地类型面积比例、各土地面积变化趋势以及普通线性回归模型得到的2035年土地类型面积比例。社会经济约束包括经济总量约束、土地开发强度约束和土地利用率约束。生态环境约束包括生态承载力和生态绿当量两方面。

科研数据管理过程具有明显的周期性和阶段性特征，基于科研活动不同阶段的数据形态和数据处理活动，形成关于科研数据管理的相关生命周期理论。重点介绍了ICPSR社会科学数据存档生命周期管理模型。

生态承载力是产量因子、均衡因子与各土地利用类型面积的乘积，用来衡量生态系统服务能力的上限[33-34]。参考相关研究，耕地、林地、草地、水域和建设用地的均衡因子分别为2.210、1.340、0.490、0.360和2.210[35]，产量因子分别为1.310、0.910、0.390、0.980和1.660[36-37]，研究区的生态承载力应高于2015年(表2)。

生态绿当量主要用于衡量某一区域生态环境的优劣。因地制宜确定本研究区的耕地、林地、草地均为拥有绿当量的土地。生态绿当量计算公式为

(2)

式中F——生态服务功能价值系数

i等于1为耕地，i等于2为林地，i等于3为草地。

2 结果与分析

2.1 MCCA模型适用性评价

2.1.1模型精度评价

采用OA、Kappa系数、mcFoM、RE、PA和UA等参数评估MCCA模型的预测精度。OA、Kappa系数和mcFoM的计算结果分别为0.934、0.886和0.261。所有地类的PA和UA均高于0.700，且除林地和建设用地外其他土地利用类型的PA和UA均高于0.800。RE越低，模拟过程中损失的信息量越少，说明模型模拟结果与实际结果更相符。模拟结果和实际土地利用数据之间的RE均值为0.508，表明模拟过程中信息损失较少，二者之间的相似度高。各项评估指标值均较理想，这说明模型可用于预测黑河中游甘临高地区未来的土地利用变化情况。

采用MCCA模型分别计算耕地、林地、草地、水域、建设用地和未利用地的2015年模拟结果和实际土地利用分布数据较2000年的变化情况，变化程度为-1～1，值为正表示土地利用类型的面积在扩张，值为负则表示减少。可以发现所有土地利用类型的变化区域基本吻合，变化地类的位置拟合效果较好，变化的斑块数量差异不明显。各土地利用成分的模拟结果与真实情况的变化规律基本一致，仅有少数像元在变化程度上存在差异(图3，图中左图为实际的土地利用成分的变化程度，右图为模拟的土地利用成分的变化程度)。

图3 2000—2015年实际和模拟的土地利用成分变化程度Fig.3 Actual and simulated land use changes during 2000—2015

为更深入探究基于混合元胞的MCCA模型的优越性，使用PLUS模型对空间分辨率250 m的纯像元土地利用结构数据进行模拟，模拟结果的FoM、Kappa系数、OA分别为0.117、0.853、0.916，而MCCA模型的OA、Kappa系数和mcFoM的计算结果分别为0.934、0.886和0.261。据此可以确认MCCA模型的模拟结果与地表实际状况吻合度更高，能够很好地模拟复杂土地利用结构。

2.1.3MCCA模型在研究区的最适模拟尺度

不同模型对于栅格尺度的要求存在差异。LIANG等[6]采用4种空间分辨率(250、500、750、1 250 m)的混合像元验证MCCA模型的模拟精度，发现250 m的模拟精度最高。鉴于此本文将混合像元空间分辨率缩小至250～500 m之间进一步确定最佳的混合像元大小，结果表明在250、300、350、400、450、500 m的混合像元中，依然是250 m混合像元模拟精度最高。不同栅格尺度的模拟准确度评价指标值见图4，6种栅格尺度模拟结果的OA均高于0.920，RE均值均低于0.550，Kappa系数均高于0.870，mcFoM均高于0.210。随着混合像元的增大，RE均值波动上升，说明栅格尺度越大模拟过程中丢失的栅格信息越多；OA、Kappa系数和mcFoM则波动降低，说明MCCA模型的拟合准确度随着栅格尺度的增大而趋于降低。

图4 不同栅格尺度的模拟精度Fig.4 Simulation accuracy of different raster scales

2.2 2035年不同情景预测结果

预知未来的土地利用情况对维持生态系统的稳定性有重要意义，本文基于普通线性回归模型的BAD情景中2015—2035年期间耕地、水域和建设用地逐年递增(图5)。说明随着城市化进程的加快，二三产业的投资明显增多，建设用地面积迅速扩张。尤其是人口老龄化导致二孩政策开放，人口数量的增加对粮食的需求量加大，耕地面积也在逐年上升。BAD情景中草地、林地和未利用地类型逐年递减，表明生态环境在逐渐恶化。

图5 研究区2000—2015年实际的和2015—2035年预测的土地利用类型覆盖度Fig.5 Actual (2000—2015) and predicted (2015—2035) land use proportions in study area

BAD和SUD情景的耕地、水域和建设用地均较2015年有所增长，而草地则低于2015年。SUD情景中2035年的林地面积占比明显高于BAD情景，草地略高于BAD情景，而耕地、建设用地和未利用地则略低于BAD情景。SUD情景占用的未利用土地较少，说明SUD情景的经济发展状况有利于保护脆弱的生态环境(表3)。

表3 不同情景下的2035年土地利用面积占比预测结果Tab.3 Land use percentage under different scenarios in 2035

MCCA模型对2035年地类结构进行空间可视化的结果显示(图6，水域变化程度较小，故未列出)，BAD情景中研究区的耕地不断向外扩张，侵占林地和草地；SUD情景的耕地和草地变化趋势与BAD情景一致，但是占用林地的现象较少。IPCC报告指出2016—2035年全球气温和降水会逐年增加[38]，研究区的水域范围的增加验证了这一结论。随着城市化进程的加快，两种情景中城镇用地和农村用地均在原有基础向外扩张，甘州区的市区范围扩展尤为明显；高台县和临泽县的城市发展依托于交通网络，未来在城市周边的交通集散地将形成较大的卫星城。两种情景中的建设用地和耕地在未来20年将向未利用地延伸，现有的沙漠和耕地交界处在未来有极大可能转变为农耕区，这可能是经济技术发展以及政策引导的结果。但是并不是所有的建设用地均会向外扩展，MCCA模型的模拟结果可能与现实情况存在少数不一致之处，比如临泽北部的工矿用地未来可能会萎缩，而模拟的结果却显示该区域将会扩张，这可能是因为MCCA模型模拟过程中受到驱动因素时间滞后性的影响。

图6 2015年实际土地利用结构、2035年模拟土地结构和2015—2035年土地利用成分变化程度Fig.6 Actual land use structure in 2015, simulated land use structure in 2035, and change of land use percentage during 2015—2035

2.3 2015年与2035年的土地资源价值核算

通过量化2035年BAD情景下的生态效益、经济效益、社会效益(表4)可以发现，2035年综合效益达3 504 943万元，而2015年仅为2 930 312万元。BAD情景下耕地、水域和建设用地的3种效益均增加，且耕地及建设用地的效益增速较水域快。西部大开发战略的实施带动经济发展，企业投资增加，建设用地随之扩张。国家人口政策的实施带动了未来总人数的上涨，粮食需求量大增，农用地面积增加。BAD情景虽然经济效益快速增加，但是生态效益增速较慢，社会效益增速居中。BAD情景下林地生态效益、经济效益和社会效益均明显降低。

表4 研究区不同情景效益对比Tab.4 Benefit comparison of different scenario in study area 万元

2015年研究区建设用地的经济效益、社会效益分别为854 144、223 213万元，而2035年SUD情景下则分别增至1 077 746、281 647万元。黑河中游甘临高地区的经济总量持续增长，建设用地的范围大幅扩张。2035年耕地面积得到了优化，耕地的各项效益高于2015年。但是耕地并未盲目扩张，占用林地的现象较少，增加的耕地主要是未利用地和草地演变而来，这有利于保障资源的可持续利用。2035年SUD情景下草地和未利用地的效益略低于2015年，而林地与2015年持平。这说明政府在发展经济的同时注重保护环境。

2035年SUD情景和BAD情景中林地的生态效益分别为36 057、22 752万元，SUD情景的生态效益较BAD情景增速更快。SUD情景较BAD情景更好地协调了生态、经济和社会三者之间的关系。SUD情景建设用地和耕地各效益的增速低于BAD情景，而林地的三大效益高于BAD情景，工农业用地的适度增加有助于维持社会稳定、保障粮食安全和保护生态环境。尤其是西北干旱区生态脆弱，环境问题突出，保护环境势在必行。

3 讨论

3.1 基于混合元胞的MCCA模型模拟精度和时间跨度分析

MCCA模型能够自下而上生成优化的土地利用空间格局，改进了传统的纯像元CA模型，混合元胞的转化规则是定量的，其邻域和元胞空间的状态具有混合、连续和多维的特点，与传统元胞土地利用成分离散、一维和定性相反，是从定性的静态模拟到定量的动态模拟的跨越[6]。因此，MCCA模型较基于纯净元胞的PLUS模型的模拟准确度更高，本文的研究结果证实了这一点。而LIANG等[6]将MCCA模型与基于纯净像元的未来土地利用(Future land use simulation model，FLUS)模型进行对比分析同样得出基于混合元胞的MCCA模型模拟准确度更高。

受到现实情况的制约，MCCA模型未统一两期现有数据的时间跨度与现有数据-预测年份数据的时间跨度，PLUS、FLUS和CLUE-S等空间模拟模型也存在类似问题。这类模型通常是基于土地利用栅格数据以及驱动力栅格数据采用回归算法(如随机森林回归、人工神经网络和逻辑斯蒂回归等算法)挖掘不同地类与驱动力因素之间的关系[6,39-40]；然后基于现有的土地利用数据中各地类的面积，采用数量预测模型(如马尔科夫链、灰色预测模型、多目标线性优化模型和系统动力学模型等)预测未来指定年份的地类面积，因此时间跨度的不统一在地类面积数值预测方面存在合理性[41]；最后基于基期土地利用栅格数据、地类变化可能性和未来指定年份的地类面积进行预测，获得未来的土地利用栅格数据。当研究中可获得栅格数据之间时间跨度过短时，考虑到社会发展的惯性及预测年份数据的参考作用，本研究未统一时间跨度。

3.2 不同情境效益分析和MOP模型优势

多目标空间优化问题较为复杂，需要同时兼顾多方面因素[42]。MOP模型采用自上而下的方式解决复杂土地系统问题，具备强大的数据处理能力，能够在综合考虑生态效益、经济效益和社会效益的基础上，从多个方面优化土地利用结构[43]。而普通线性回归模型仅从各地类面积的数量变化方面进行分析，难以实现土地优化配置的目的。将基于普通线性回归的BAD情景与基于MOP算法的SUD情景进行对比，计算不同情景下各土地利用类型的生态效益、经济效益和社会效益，发现2035年土地优化配置的SUD情景综合效益远高于2015年，虽然其经济效益和社会效益略低于BAD情景，但生态效益明显高于2015年，这有利于缓解当地脆弱的生态环境。而将MOP和MCCA进行耦合模拟2035年研究区的土地利用空间结构，能够充分利用数量模拟模型和空间模拟模型的优势，实现未来土地利用结构的可视化，为相关部门进行土地规划和生态修复提供参考。

4 结论

(1)采用Kappa系数、OA、PA、UA、RE和mcFoM等指标进行精度验证，各项指标的计算结果均说明MCCA模型的模拟准确度高，且MCCA模型的模拟精度明显优于基于纯净元胞的PLUS模型，这验证了MCCA模型在黑河中游甘临高地区具有高度适用性。为进一步确定最佳像元尺度，选择分辨率250、300、350、400、450、500 m的混合元胞数据进行模拟，结果表明分辨率为250 m时模型模拟精度最高，Kappa系数、mcFoM和平均RE分别为0.886、0.261和0.508，因此采用250 m的混合像元数据预测2035年的土地利用结构。

(2)在基于普通线性回归模型的BAD情景中，2000—2035年期间的耕地、水域和建设用地范围逐年增加，林地、草地和未利用地则相反，存在工农业用地占用林草地的现象。SUD情景的建设用地、草地和未利用地的变化趋势与BAD情景基本一致。但SUD情景耕地的增加速率较BAD情景偏低，其2035年林地的范围大于BAD情景，该情景中工农业用地占用林草地的情况较少。SUD情景的经济效益和社会效益略低于BAD情景，而生态效益明显高于BAD情景。这说明基于多目标线性规划(MOP)模型与MCCA模型的SUD情景能够实现土地优化配置，在发展工农业的同时保护生态环境。