全球黑土区耕地动态变化及其驱动力Logistic 回归分析

杨佳佳，吕骏超，金珊合，刘 洋，李文博

中国地质调查局沈阳地质调查中心，辽宁沈阳 110034

0 引言

土地利用/土地覆盖变化（land use and land cover change，LUCC）是全球气候变化、环境变化研究内容的重点.地球系统科学的原理告诉我们，土地资源与生态环境是紧密相连的，土地资源数量与质量的变化是生态环境变化的重要特征，生态环境信息也反映出土地资源的开发、利用状况.只有掌握土地资源数量、质量和可持续利用的综合信息，才能为区域的持续发展提供有力保障［1］.另外，土地利用变化会对区域经济增长、收入水平以及自然资源生态环境产生多方面影响，同时，这些自然、社会、经济因素变化又会反馈于土地利用系统，造成土地利用变化.所以，对土地利用变化机制的深入研究有助于决策者在不同尺度水平上实施政策干预，保证土地利用的可持续［2］.

黑土土壤肥沃，生产力高，已成为一种宝贵的农业经济资源.世界四大黑土区均为重要的农牧业生产基地，保护黑土已经是相关国家土地管理者的共识.目前，国内外学者在黑土区开展了水土流失［3-5］、土地流转［6］、土壤属性［7-8］、土壤光谱［9］、有机碳变化［10-12］等诸多研究，但针对全球黑土区耕地利用变化的驱动力研究尚无.耕地变化模拟与预测对黑土耕地可持续发展及更好地理解黑土区土地利用系统机制有重要研究意义.

土地利用空间格局变化模拟研究一是预测土地利用变化将在什么地方发生，二是研究导致土地利用变化的驱动因素.土地利用这种空间变量的增减变化无法定量化描述，一般的线性回归模型因变量必须为连续变量［13］.Logistic 回归模型作为一种非概率性线性回归，是针对二分类或多分类因变量进行回归分析时采用的分类统计方法，其中的自变量既可以是连续变量，也可以是分类变量.通过Logistic 回归分析，可以计算自变量的相关权重，确定相关驱动力的影响作用，同时还可以解释因变量的发生概率，能够有效地解决以往土地利用变化驱动力研究中存在的局限性问题.该方法应用于土地利用变化驱动力研究中在国内外已有较多资料可寻［14-18］.部分学者将元胞自动机、马尔科夫模型、神经网络等方法与Logistic 回归模型结合来预测未来土地利用的空间格局［19-23］.本研究基于全球黑土区不同时期的土地利用数据、自然环境数据和社会经济数据，系统分析2005—2019 年间耕地变化的过程和特征，构建全球黑土区耕地变化的Logistic 回归模型，定量分析耕地变化的空间驱动因素，并对模型预测能力进行检验，以期为理解全球黑土区耕地变化机制提供科学参考.

1 研究区概况及数据来源

1.1 全球黑土区耕地分布范围

美国土壤分类系统（Soil Taxonomy，ST，1999 年）目前被联合国粮农组织（FAO）、联合国教科文组织（UNESCO）以及全球45 个国家普遍采用.该分类系统将所有土壤分为12 个土纲，其中暗沃土纲（Mollisols）定义为草原生态土壤，以深厚、暗色表层为显著特征，因其英文名称来自拉丁语，这一肥沃的表层也称为软土层，是源于植物根系的有机质长期加入的结果.暗沃土表层厚度不小于25 cm，该层有机质含量大于或等于1%（C有机≥0.6%），盐基饱和度（NH4OAc 法）大于或等于50%，润态的明度和彩度小于或等于3［24］.其定义与中国土壤系统分类（Chinese Soil Taxonomy，CST，1999 年）中均腐土纲（Isohumosols）、世界土壤资源参比基础（World Reference Base for Soil Resources，WRB，1998 年）中黑土（Phacozems）、黑钙土（Chemozems）、栗钙土（Kastanozems）的定义比较接近但又有区别.暗沃土被视为广义的黑土，按照其定义，全球共分布四大片黑土区，包括北美洲密西西比河流域地区、中国松辽流域地区、乌克兰与俄罗斯大平原以及阿根廷和乌拉圭的潘帕斯大草原地区［25］（图1），其分布区内耕地将作为本次全球黑土研究的明确范围.

图1 基于暗沃土纲（Mollisols）的全球黑土区耕地分布范围Fig.1 Global cropland distribution in black soil（Mollisols）regions

1.2 全球黑土区耕地分布面积

以WGS84 为坐标系，采用正弦投影（sinusoidal projection），依据2019 年MODIS 土地覆盖类型三级数据500 m 产品（MCD12Q1）统计面积如下：全球黑土区耕地分布总面积为375.66×104km2，在六大洲皆有分布，但不均衡，其中欧洲分布147.32×104km2，占比39.22%；北美洲分布119.43×104km2，占比31.79%；亚洲分布74.83×104km2，占比19.92%；其他地区包括南美洲分布25.24×104km2，非洲分布7.44×104km2，大洋洲分布1.40×104km2，总计占比9.07%（图2）.

图2 全球各大洲黑土区分布占比Fig.2 Distribution percentage of black soil areas by continents

1.3 项目数据来源

研究采用的数据说明：MODIS 三级数据描述土地覆盖类型全球500 m 产品（MCD12Q1），根据一年的Terra 和Aqua 卫星观测所得的数据进行处理，信息提取采用监督决策树分类方法.选取国际地圈生物圈计划（IGBP）作为全球植被分类方案，该土地覆被数据集中包含了17 个主要土地覆盖类型，其中包括11 个自然植被类型，3 个土地开发和镶嵌的地类和3 个非草木土地类型定义类.

驱动力因子来源：高程（DEM）数据来源于美国航空航天局（NASA）提供的SRTM3 数据，精度为90 m；世界人口密度来源于WorldPop 发布的地图，是目前比较可靠的长时间序列数据；气象因子数据来源于美国国家海洋和大气管理局（NOAA）官方网站；道路、河流数据来源于Google 全球地图.

2 研究方法

2.1 土地利用数据预处理

根据LUCC 分类体系（参考中科院分类标准），将研究区土地类型合并为6 个一级地类：耕地、草地、林地、水域、建设用地和未利用地.其中将IGBP 中常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林、混交林、郁闭灌丛、稀疏灌丛7 个地类合并为林地类；将多树草原、稀树草原、典型草原3 个地类合并为草地类；将永久湿地、农田/植被（由农田、森林、灌木丛和草原组成的土地，其中没有一个组成部分占景观的60%以上）、冰雪、裸地或低植被覆盖地4 个地类合并为未利用地类；剩下3 种地类分别为水域（水体）、耕地、建设用地（城市和建筑）.

2.2 Logistic 回归模型

Logistic 回归模型是一种因变量与自变量之间为非线性关系的分类统计方法［13，26］.其中，因变量类型为二分类、多分类变量，自变量类型为连续和分类变量，最终以事件发生概率的形式提供结果.

以耕地为例，运用Logistic 逐步回归方法对耕地分布的驱动因子进行分析.此方法解决了一般线性回归方法无法计算自变量为分类变量的问题，对于每个解释变量生成单独系数，并通过一定的权重运算法则被解释为耕地类型的变化概率［13，27］.

式中：pi表示空间上每个图斑可能转为耕地的概率；α 为截距；xi表示相应各候选解释变量；β 为回归系数；k 为解释变量的种类.

发生比（odds）是发生频数与不发生频数之间的比，用来对各种自变量的Logistic 回归系数进行解释，可以较直观看出自变量对事件概率的作用，运用在本研究中可反映出各解释变量对耕地影响程度的差异.

式中：odds 为发生比，即每个图斑成为耕地的概率和成为其他土地类型的概率之比.由式（2）可知，每当驱动因子增加一个单位时，发生比相应变化exp（β）倍.exp（β）用来解释驱动因子的作用，即为发生比率.

3 耕地变化驱动力研究

3.1 基于Logistic 模型的耕地变化驱动力分析

3.1.1 耕地数量动态变化

2005—2019 年间，全球黑土区耕地面积发生很大变化.耕地由2005 年的434.42×104km2锐减到2019年的375.66×104km2，共计减少58.77×104km2.其中2010—2015 年间耕地面积消减最大，共计减少58.19×104km2.从全球分布情况来看，亚洲黑土区耕地减少占比31.70%，北美黑土区耕地减少占比29.46%，南美黑土区耕地减少占比27.79%，其他黑土区耕地共计减少11.06%.从地类转变情况来看，其中有76.04%的耕地转变为未利用地（包括永久湿地、农田/植被、冰雪、裸地或低植被覆盖地），有19.95%的耕地转变为草地，1.94%转变为建设用地，1.90%转变为林地，0.16%转变为水体.造成这种情况的原因有很多，比如耕地退化、土地灾害等自然条件变化因素，生态退耕、土地结构调整、政策干预等人为干涉因素.

3.1.2 模型建立与抽样过程

模型建立所需的因变量（耕地变化状态）与自变量（影响因子）如表1 所示.将2005—2010、2010—2015、2015—2019 年3 个时期的耕地分布数据及各驱动因子用ArcGIS 软件组织成匹配的空间数据，针对各时期内的耕地变化情况分别建立Logistic 模型.为减弱数据空间自相关性，依据两种因变量结果（耕地增加赋值为1、耕地减少赋值为0）在设定20 km 距离的条件下随机抽取30%的点，再用这些点去提取相应自变量的值，通过异常值去除、多重共线性诊断后导入SPSS 进行分析.SPSS 提供了卡方值（wald）统计量作为回归系数的估计，表示在模型中每个解释变量的相对权重，用来评价每个解释变量对事件预测的贡献力［28］.需要说明的是：在Logistic 模型中设置坡度Ⅱ（2～6°）、坡度Ⅲ（6～15°）、坡度Ⅳ（15～25°）、坡度Ⅴ（＞25°）为哑变量，坡度Ⅰ（≤2°）作为它们的参照水平；设置坡向Ⅱ（阳坡）、坡向Ⅲ（半阳坡）、坡向Ⅳ（半阴坡）、坡向Ⅴ（阴坡）为哑变量，坡向Ⅰ（平坡）作为它们的参照水平；在对DEM、温度、降水提取过程中采用了双线性插值方法.

表1 耕地利用变化及影响因子变量Table 1 Cultivated land use changes and impact factor variables

3.1.3 结果分析

模型似然比检验对整体模型进行有效性分析发现（表2）：各模型p 值均小于0.01，模型有效；对比各阶段模型AIC 值（衡量统计模型拟合优良性的一种标准），2015—2019 年阶段模型构建最好（AIC 值越低越好）.

表2 Logistic 回归模型似然比检验结果Table 2 Results of logistic regression model likelihood ratio test

3 个阶段模型中影响耕地变化的驱动因子类型变动较大，其中到最近道路距离、DEM、温度分别出现2次，降水、坡度Ⅲ、到最近河流距离分别出现1 次.另外，在3 个阶段模型中，人口密度变化这一变量都没有进入模型，可能的原因是该人口密度数据未能反映从事农业人口变化的真实情况（表3）.

表3 耕地变化的logistic 回归模型相关系数Table 3 Parameters of logistic regression model for cultivated land change

第1 阶段（2005—2010 年）：本阶段耕地变化重要的解释变量按照重要程度（Wald 统计量）排列顺序是到最近道路距离、降水、DEM、温度、坡度Ⅲ.其中，到最近道路距离、DEM、温度、降水会对耕地变化产生显著的正向影响关系，坡度Ⅲ会对耕地变化产生显著的负向影响关系.其中，到最近道路距离每增加1 km，其他地类转为耕地的概率将增加0.99 倍（即1 km 的发生率是e0.011=1.011），分析原因可能是在这一阶段，随着农业经济发展，道路建设投入增加成为耕地变化的主导原因.

第2 阶段（2010—2015 年）：本阶段耕地变化重要的解释变量按照重要程度排列顺序是到最近河流的距离、DEM，该两项因子均会对耕地变化产生显著的负向影响关系.其中，到最近河流距离每增加1 km，耕地转为其他地类的概率将增加1.01 倍，分析原因可能是在这一阶段，随着区域社会、工业化与城镇化的发展，河流附近土地利用的转变成为耕地变化的主导原因.

第3 阶段（2015—2019 年）：本阶段重要的解释变量按照重要程度排列顺序是温度、到最近道路的距离，该两项因子均会对耕地变化产生显著的正向影响关系.其中，温度每增加1 ℃，其他地类转为耕地的概率将增加0.84 倍，分析原因可能是在这一阶段，区域社会的发展趋于稳定，人为因素的干扰作用在减退，温度的变化反而成为耕地变化的主导原因.

3.2 耕地转入模型空间格局模拟

利用2005—2010 年的Logistic 耕地转入模型，以2005 年为基期更新驱动因子，来检验在2005—2019年阶段可能转为耕地的准确性，并将预测结果与2005—2019 年实际变化的耕地进行空间比对，检测模型的预测能力.根据2005—2010 年耕地转入模型运行结果，耕地转入概率P 的预测模型为：

其中，Y=ln［P/（1-P）］=-1.206＋0.011×距道路距离＋0.848×DEM＋0.024×温度＋1.750×降水-1.311×坡度Ⅲ.

通过研究区耕地转入概率模拟结果与2000—2019 年耕地变化分布的对比发现（图3、4），2005—2019 年间呈片状增加的耕地大部分落入耕地转入预测概率高值区，其间在实际转入的耕地中，有71.64%图斑的预测概率值在0.5 以上，52.07%图斑的预测概率值在0.6 以上.由于土地利用系统的复杂性，土地利用空间分布的预测准确性很难达到85%以上［2，29］.可见，虽然2 个阶段影响耕地变化的解释变量不同，但从模拟的结果来看，该模型对转入的耕地预测效果较好，说明在驱动过程保持不变或者近似的条件下，近期的耕地增加模型可以有效地预测未来一定时期内耕地转入空间分布概率.

图3 2005—2019 年中国黑土区耕地转入预测结果及概率分布Fig.3 Predicted probability distribution of cultivated land in black soil regions of China during 2005-2019

图4 2005—2019 年北美黑土区耕地转入预测结果及概率分布Fig.4 Predicted probability distribution of cultivated land in black soil regions of North America during 2005-2019

4 结论

1）全球黑土区3 个阶段（2005—2010、2010—2015、2015—2019 年）耕地变化显著，耕地变化特征先后出现小幅增加—大幅减少—小幅增加的趋势.14 年间耕地共减少58.77×104km2，其中2010—2015 年间耕地面积消减58.19×104km2.从全球分布情况来看，亚洲黑土区耕地减少幅度最大，占比31.70%.从地类转变情况来看，耕地转变为未利用地占比最大，约占76.04%.

2）全球黑土区耕地变化在各个阶段有着不同的重要驱动因素.在2005—2010 年阶段最主要的驱动因素是到最近道路距离、降水和DEM；2010—2015 年主要的驱动因素是到最近河流距离和DEM；2015—2019年主要的驱动因素是到最近道路距离和温度.通过对耕地变化Logistic 回归模型的检验，得出Logistic 回归模型进行耕地变化驱动力分析具有可构建性、可操作性等优点，能较好地揭示全球耕地变化过程的主要驱动力.同时，本研究选取的自然因子、社会因子指标中还缺乏一些重要变量，比如土壤类型、地貌类型、国内生产总值（GDP）等，在一定程度上影响了模型的精度.

3）全球黑土区耕地变化是一个相当复杂的系统过程，与全球各地区的经济、社会、人口、技术、自然因素、政治结构与价值观念等直接相关，有些因素难以量化纳入模型，因此为了提高耕地变化的预测概率，我们应考虑建立驱动力的动态模型.这些驱动力应包括政策的制定与调控、农产品的价格变化等动态变化因素，这样，代表变化可能性原因的驱动力在空间模型中就可以显现.因此，从这方面选取一些可以量化的指标应是下一步研究的重点内容.