黑龙江省耕地系统安全预警及其驱动因素空间分异

宋 戈，陈藜藜，邹朝晖

0 引 言

耕地系统安全预警是科学把握区域耕地安全状态的有效途径。当前国内外耕地系统安全预警研究多以耕地资源为研究对象[1-2]，主要集中在预警理论[3-5]、预警评价框架体系[6-9]以及预警方法[5,7,10-11]方面，而以耕地系统作为研究对象进行安全预警研究极少。耕地系统是以耕地资源为主体的自然生态系统和以人类活动为主体的社会经济系统在特定区域内通过协同和拮抗作用而形成的复合系统，是耕地资源子系统自身的脆弱性（隐患）以及受人为活动影响（压力）产生波动时，通过自我调解与自我恢复（免疫）以维持自身稳定，同时人类社会也为减少和避免耕地环境风险，通过经济、环境、土地等政策、决策或管理措施做出调控（响应）的自组织复合系统[12]，其影响系统变动因素关系十分复杂，即使是单因子变化也会引起其有关因子的“共振”，从而导致系统整体的重大变动[13]，使耕地系统安全受到威胁，发生警情，更重要的是不同区域耕地系统安全程度又存在较大的差异[14]，而当前对耕地系统安全预警驱动因素空间差异的研究极少，亟需加强。目前研究大多利用预警评定指标作为驱动因素进行分析，忽略了直接影响因素与根本影响因素的区别，导致无法识别警情真正根本影响因素。同时，现有的预警研究多在宏观尺度上开展，且仅对警情进行评价，没有落实到空间上，在一定程度上不能满足因地制宜保护耕地的要求，亟待明确耕地利用系统安全预警的空间分布状态和空间分布特征，为有针对性的进行耕地保护提供科学依据。

黑龙江省是中国东北粮食主产区，也是国家重要的商品粮基地，商品粮率高达 75%，承担着保障国家粮食安全的重任。尽管近年来黑龙江省耕地面积呈增长的趋势，但其耕地土壤污染、耕地生态环境破坏以及农药和化肥施用等问题凸显[15-16]，严重威胁着耕地系统安全，识别和掌握黑龙江省耕地系统安全预警状态，可为制定耕地系统安全保障措施提供科学依据。鉴于此，在课题组前期省域尺度下耕地系统安全预警体系构建的基础上[12]，本文以黑龙江省各县域为研究单元，从系统角度对 2014年研究区耕地系统安全预警进行评定，将预警状态落实到具体空间，采用探索性空间数据分析方法，分析预警状态的空间差异特征；并从区分预警评定指标（直接影响因素）与驱动因素（根本影响因素）的角度，从自然生态和社会经济方面选取耕地系统安全预警的根本影响因素作为驱动因素，引入地理加权回归模型，探寻耕地系统安全警情的驱动因素及其空间差异，弥补当前耕地系统安全预警研究无法准确判断各驱动因素对其空间差异化作用的缺陷，为制定更为精细的耕地系统安全警情缓解措施及农业可持续发展的决策提供科学依据。

1 研究区及数据来源

黑龙江省土地面积 47.3万 km2，占全国土地的4.9%，是世界三大黑土带之一。耕地面积 1 594.09万hm2，占全国耕地面积的十分之一以上，居全国各省耕地资源总量第一，是全国第一产粮大省和重要的商品粮生产基地，承担着保障国家粮食安全的重任。但由于掠夺式生产和经营，毁林开荒，重种轻养，过度施加化肥、农药等，导致耕地养分收支长期处于赤字状态，耕地黑土层逐渐变浅、变薄，土壤肥力退化[17-18]，耕作性能和抗旱、抗涝能力逐渐降低，严重威胁着耕地质量、农民收益和农业可持续。

本文研究区主要包括黑龙江省各县域和市辖区共 72个（除大兴安岭地区外，下同）（图1），数据来源主要包括4方面：1）经济社会数据主要来源于2015年黑龙江省12个地级市统计年鉴、统计公报及环境公报；2）DEM数据、降水量、土壤侵蚀和水土流失数据主要来源于中国科学院地理科学与资源研究所(http://www.resdc.cn/Login.aspx)；3）耕地后备资源来源于黑龙江省县级耕地后备资源调查评价数据库；4）各县域（市辖区，下同）土地利用数据来源于2015年土地利用变更数据库。

图1 黑龙江省行政区划图Fig.1 Administrative map of Heilongjiang Province

2 研究方法

2.1 耕地系统安全预警指标体系构建

基于本课题组前期研究[12]，本文构建基于抗原-抗体的县域尺度下黑龙江省耕地系统安全预警评定指标体系，分为目标层、准则层、要素层和指标层（表1）。

2.2 探索性空间数据分析

本文采用探索性空间数据分析（exploratory spatial data analysis，ESDA）的全局和局部空间自相关指数（Global Moran’s I和 Local Moran’s I）分析黑龙江省县域耕地系统安全预警的空间关联性和差异程度。耕地自然地理条件的相似性和差异性，以及人口的区域流动性和农产品的区域流通性，表明一定区域范围内耕地系统安全预警状况可能存在较强的空间关联。本文分别采用Global Moran’s I（I）和 Local Moran’s I（Ii）指数检验区域总体和内部空间之间各自的相似性、差异程度以及不稳定性[19-22]。

式中xi和xj分别表示县域i和j的观测值，x为均值，wij为空间权重矩阵，n为空间单元数；I和Ii取值范围都为[-1,1]，若I (Ii) > 0，则表示观测值在区域空间内存在集聚特征（高高集聚或低低集聚，分别表示为H-H型和L-L型）；若I (Ii) < 0，则表示观测值在区域空间内存在异质集聚特征（高低集聚或低高集聚，表示为H-L型和H-L型）；若I (Ii) = 0，则表示观测值在区域空间不存在空间关联，相互独立。同时，I (Ii)越趋近于1，研究区域空间集聚特征越显著。

2.3 最小二乘法和地理加权回归

本文引入地理加权回归（geographical weighted regression，GWR）[23-24]进行耕地系统安全预警根本驱动因素分析，为判别耕地系统安全预警的根本驱动因素提供更为准确结果及空间信息。耕地系统安全预警的驱动因素众多，自然地理要素和社会经济要素都可能具有空间上的相关性和差异性，准确把握局部空间变异特征尤为重要。但是否应对预警的驱动因素进行 GWR分析有待检验，若各驱动因素的回归参数不受空间位置影响，保持固定不变，则应使用普通最小二乘法（ordinary least squares，OLS）[25]回归模型，计算公式为

式中y为因变量，x为自变量，β0为回归常数，β1~βm为回归系数，ε为残差，m为研究县域的个数。

与传统线性 OLS模型相比，GWR模型具有以下优点：1）克服了OLS模型因无法对局部参数估计而不能反映区域局部情况的缺陷；2）能够识别两个或多个变量的多重共线性问题，使回归结果更为符合客观实际；3）借助 GIS估计模型的参数，以空间表达，有利于进一步构建地理模型，探索空间变异特征和空间规律[26]。GWR模型结构如下

表1 县域尺度下黑龙江省耕地系统安全预警指标体系Table 1 Early-warning index system of cultivated land system security under county level in Heilongjiang province

式中(ui, vi)是第i个样本的地理中心坐标，β0(ui,vi)和βk(ui,vi)分别为第 i个采样点的统计回归常数项和第 k个回归参数；xik是第i个采样点上第k个变量；p为某一采样点上参与回归的变量个数；εi为误差项。对于 βk(ui,vi)来说，越临近的数据会有更大的影响。

GWR模型的确定取决于权重函数的确定和最优带宽的选择。空间权重的确定则决定了空间地理位置对空间因素的影响规则，本文采用高斯函数来确定权重函数

式中b是带宽，dij为样本点(ui, vi) 到样本点(uj, vj)的距离。如果i的数据被观测，则其他观测值的权重会根据高斯曲线，随着距离dij的增加而降低。同时采用更能适应不同模型存在不同自由度特性的赤池信息量准则AICc（akaike information criterion）[24]确定最优带宽。

3 黑龙江省耕地系统安全预警空间分异分析

3.1 各县域耕地系统安全预警空间分异特征

采用改进的级数突变模型[12]对2014年黑龙江省各县域耕地系统安全预警状态进行评定，得到各县域耕地系统安全预警值，预警最低值为0.314 6，其次为0.322 6，分别出现在牡丹江东宁县和伊春市辖区；最高值为0.511 1，出现在大庆市辖区。在预警结果基础上，借助ArcGIS10.2自然断点法和IBM SPSS 23.0层次聚类分析法[27]对预警值（early-warning value, EWV）进行分类，确定预警为4个等级：较安全（0.2227 < EWV < 0.3533）、轻警（0.3534< EWV < 0.3990）、中警（0.3991 < EWV < 0.4502）和重警（0.4503 < EWV < 0.5329）（图2a），2014年，全省的72个县域中，有69个县域处于预警状态，处于轻警、中警和重警状态的县域分别为28、32和9个，仅有3个县域处于较安全状态，整体警情水平较高，耕地系统安全状况不容乐观。

从预警状态空间格局可知，整体表现为西部警情水平最高，东部次之，南北轴带地区（包括中部）最低。西部属于经济较发达区域，人类活动频繁，经济扩张带来了建设占用耕地面积的升高，对耕地系统安全的干预程度较高；同时，西部水土流失较为严重，土壤侵蚀程度较高，导致耕地系统安全警情水平较高。在南北轴带区域中，北部和南部地势较高，森林和草地较多，有助于促进良好的耕地环境，加强耕地生态恢复力，提升耕地系统安全，降低其警情水平。东部区域主要属于三江平原地带，自然地理条件优越，高标准农田建设力度较大，但耕地系统也受人类活动影响较大，且东临黑龙江，由于排水条件限制，当降雨季节，临江县域耕地受洪涝灾害严重，导致耕地易涝面积较多，土壤侵蚀加剧，致使此区域耕地系统安全警情水平偏高。此外，各市辖区及周边县域的警情水平较高于远离市区的县域预警水平（图 2a）。越近市区，区域经济越发达，尽管对耕地系统投入和保护的响应程度可能越高，但耕地系统受到人类社会干预程度也越高，可能导致耕地系统警情水平越高。

图2 2014年黑龙江省各县耕地系统安全预警状态空间格局和局部空间关联聚集Fig.2 Spatial pattern and local spatial agglomeration of early-warning status of cultivated land system security in Heilongjiang Province in 2014

3.2 各县域耕地系统安全预警空间自相关分析

基于Geoda1.6.7和ArcGIS10.2软件，根据各县域之间的邻接关系，分别采用Rook、Queen 和Bishop邻接方式，构建空间邻接矩阵（一阶邻接）。计算结果表明，基本符合正态分布特征的仅有Rook邻接关系的空间邻接性频率直方图，因而以Rook邻接原则确定空间权重，对黑龙江省2014年各县域耕地系统安全预警值进行全局空间自相关分析，得Moran’I为0.377，Z值为3.769，在0.01显著性水平下通过检验，表明黑龙江省各县域耕地系统安全预警值存在高度全局空间正相关性。进一步运用局部空间自相关分析研究区内部的集聚特征（图2b），从聚集结果类型看，正相关类型与负相关类型在空间分布上呈现鲜明的对比特征，正相关类型以“组团”形式凸显，H-H型和L-L型县域分别有11和10个，共占总县域个数29.17%。H-H型主要聚集在西部的齐齐哈尔市和大庆市，L-L型主要聚集在南北轴带上的中部区域，主要为伊春市和哈尔滨的北部和东部，及牡丹江北部；L-H型和H-L型县域较少，分别为4和5个，且分别零星分散在H-H型和L-L型周围。在实际中，研究区西部多年来土壤侵蚀和水土流失较为严重，且大庆县域多盐碱耕地；随着经济增长加速，工业污染和环境恶化问题凸显，在高度集约化的现代农业化生产下，农业污染逐渐影响到耕地土壤生产能力，促进了耕地系统安全警情的升高，与H-H型耕地系统安全预警值聚集区在空间上相吻合。在中北部的伊春市，森林面积高达 86.9%（来源于 2014年土地利用变更数据库），具有较好的水土条件，粮食生产稳定，耕地生态环境良好。南部的牡丹江处于农林交错地带，森林面积也较为丰富，耕地利用强度适中，粮食产量较高，耕地压力指数较低。更重要的是，近年来耕地生态修复和退耕还林举措加速了耕地系统安全状况的提升，这与L-L型在空间上分布较为一致。

4 黑龙江省耕地系统安全预警驱动力空间分异分析

4.1 驱动因素选取

本文从区分耕地系统安全预警的直接表征指标（预警评定指标）与其根本驱动力角度出发，并从系统、宏观角度考虑耕地系统安全预警根本驱动因素，结合根本驱动因素对耕地系统安全预警直接表征指标的影响，将根本因素分为自然因素和社会经济因素。考虑黑龙江省作为粮食生产大省的特点，本文在初期选取高程（elevation，EL）、坡度（slope，SL）、年降水量（annual average precipitation，AVP）、年均气温（annual average temperature，AVT）、森林覆盖率（Forest coverage，FC）5个变量为自然因素，取人口自然增长率（natural population growth rate，NPGR）、人均 GDP（per capita GDP，GDPP）、城镇化水平（urbanization level，UL）、单位面积农业固定投资（agricultural investment in agricultural fixed assets，AIN）、单位面积农业从业人员数（agricultural employees，PE）、单位面积机械总动力（the total power of agricultural machinery，TPAM）、水土协调度（water-soil coordination degree，WSCD）、单位面积农业三废指数（agricultural wastes index，AWI）、路网密度（road network density，RND）9个变量为经济社会因素。考虑到有些因素之间数量级相差较大，可能导致结果偏差，因而对因变量及14个自变量无量纲化。以2014年黑龙江省各县域耕地系统安全预警值作为因变量，初步对各县域预警值与14个变量进行相关性分析，结果显示NPGR和TPAM与各县域预警值相关性不显著，故本文在后续分析中删除这2个变量。

4.2 OLS模型结果

为初步判断各因素对黑龙江省县域耕地系统安全预警的整体影响程度，检验驱动因素的分析是否适应GWR模型，并为避免由于多重共线性而导致GWR模型出现设计错误提供可靠依据，应首先进行OLS模型构建和分析。因此，以2014年黑龙江省各县域耕地系统安全预警值为因变量，以12个自然因素和经济社会因素作为解释变量，构建 OLS模型 1。从回归结果可知，解释变量 AVP和FC的方差膨胀因子（variance inflation factor，VIF）[28]值都远大于7.5，表明存在冗余解释变量，因此逐步移除解释变量AVP和FC，直至冗余消失，最后得到OLS模型2（表2）。以校正赤池信息量准则AICc最小和校正的决定系数（R2）最大为最优的选择原则，则模型2优于模型 1，模型 2整体上通过了 1%水平的显著性检验，且Jarque-Bera检验结果表明模型没有出现偏差，其结果显示，EL、AVT、GDPP、UL、AWI和RND 6个因素对县域整体预警值有显著的正向影响，而SL、AIN 和WSCD 3个因素有着显著的负向影响。

表2 OLS估计结果Table 2 Results of ordinary least squares estimation

随着高程和坡度的增加，耕地耕作难度加大，人们对耕地管理和保护的程度越小，灌溉排水条件越加欠缺，但同时人类活动对耕地干扰的负面影响越小，对耕地系统安全预警的正负影响取决于二者中某一方面的偏向，这可能是EL和SL分别与耕地系统安全预警值成显著正相关和负相关的原因。GDPP和UL是经济发展状况的体现，其与耕地系统安全预警值呈显著的正相关，这可能是经济的发展带来高强度的农业现代化和城市的扩张，高强度的农业现代化所带来的过度干预降低了耕地系统安全状况；且城市的扩张带来耕地的占用，增加耕地压力，也增加了耕地系统安全警情程度。有研究表明，农业化肥、农药和塑料薄膜的使用在短期内可能增加粮食产量，但长期累积将导致耕地土壤板结，水土失衡，破坏土壤微生物的生存环境，最终影响耕地的粮食生产能力[29-30]。近年来，黑龙江省农业化肥、农药和塑料薄膜三废使用量明显增加，对耕地系统安全的负面影响日益凸显，这可能是AWI与耕地系统安全预警值成显著正相关的主要原因。路网密度的增加，增强了人们耕作的便捷度，可能增加对耕地的保护力度，但也在一定程度上加深了人类活动对耕地的影响程度，RND与县域预警值整体上成正相关，即表明人类活动对耕地造成的负面作用大于增加路网密度带来的对耕地保护的正向作用。农业固定资产投资的增加，可增强对耕地的保护力度，提升耕地系统安全程度，因而AIN与县域预警值整体呈现负相关关系。水土协调度是有效灌溉耕地面积与耕地总面积的比值，其比例的增加，有利于维持农作物所需的土壤水分，促进作物生长，也为耕地良好的土壤环境提供了条件，这可能是 WSCD回归估计系数为正的主要原因。OLS模型回归估计结果符合理论与实际，但 Koenker（BP）统计量显著性结果表明解释变量在空间上是非平稳的，应进一步采用 GWR模型进行局部空间预警驱动因素回归分析。

4.3 GWR模型结果

OLS回归分析结果为地理加权分析提供了可行依据，以OLS模型2最终通过检验的10个因素为解释变量，以2014年黑龙江省各县域预警值为因变量，并依据方程（4）构建GWR模型，以自适应型空间核回归为基础，以固定高斯函数（fixed Gaussian）为权属函数，以AICc方法确定最优带宽，运用GRW 4.0软件，进行回归计算并得到估计结果（表 3）。结果表明，EL、SL和 AVT 3个自然因素和GDPP、UL、AIN、WSCD、AWI 和RND 6个社会经济因素对耕地系统安全预警值有显著的影响。与OLS模型相比，GWR模型的R2和Adjusted R2值都有所提高，且模型拟合显著性标识 AICc值从 OLS模型2的-101.776下降至-109.065，数值差大于 3，表明 GWR模型改进回归的拟合结果优于OLS模型。对比表2和表3，可发现GWR模型估计的10个变量局部系数的平均值与OLS模型估计的全局变量系数值相差较小，说明OLS模型仅能反映驱动因素在全局空间上的平均相关联水平，而无法反映驱动因素作用的空间差异特征。GWR模型可估算每个变量对每个空间单元的影响系数，估计结果可以更为详细地反映不同因素对耕地系统安全预警影响的空间差异。

表3 地理加权回归模型估计结果Table 3 Results of geographical weighted regression estimation

4.4 预警驱动因素的空间差异分析

借助ArcGIS，将不同驱动因素作用的空间差异直观表达（图3）。黑龙江省各县域的高程和坡度（图3a和3b）分别与耕地系统安全预警值呈正相关和负相关。高程回归系数的正高值和坡度回归系数的负高值集聚区都出现在研究区西北部，且由西北部向东南部回归系数绝对值逐渐递减，表明高程和坡度对耕地系统安全预警状态的影响呈现由西北地区逐渐向东南部地区减弱的趋势。

年均气温与耕地系统安全预警值呈正相关关系。回归系数空间分布上呈现西部和南部高，东部和北部低的特征（图3c），表明西部和南部年均气温对耕地系统安全预警值的影响较大。回归系数较大值出现在哈尔滨的市辖区、五常市和尚志市等，大庆市的市辖区、肇源县和肇州县等，这些区域人口密度大，经济较其他区域发达，人类活动频繁，年均温度高于其他区域，这可能与温度的升高（尤其在夏季）不利于农作物和土壤微生物的生长有关；回归系数较小值出现在佳木斯市的同江市和抚远市，黑河市辖区和孙吴县，鹤岗市的嘉荫县等，这些地区临近黑龙江流域，受水域比热容大的影响，气温不会太高，对土壤微生物生长，农作物的生长都没有明显抑制作用，因而与耕地系统安全预警值的相关系数较低。

水土协调度与耕地系统安全预警值呈负相关关系。回归系数绝对值由东部向西部逐渐递减（图3g），说明东部耕地系统安全预警受水土协调度影响较大，这主要是因为东部县域临近黑龙江流域，充沛的水资源为耕地土壤水分的获取提供了便利条件，促进了有效灌溉耕地面积的增加，水土匹配较优，提升了耕地系统安全程度，降低了其警情水平。

图3 2014年黑龙江省各驱动因素回归系数空间分布Fig.3 Spatial distribution of regression coefficient of factors in Heilongjiang Province in 2014

人均 GDP、城市化水平、农业三废指数以及路网密度都与耕地系统安全预警值成正相关关系（图3d、图3e、图 3h和图 3i），而单位耕地面积农业机械总动力与预警值成负相关（图3f）。从5个变量回归系数的空间分布来看，其绝对值高值都在西部和西北部聚集，主要分布在齐齐哈尔市和大庆市各县域，路网密度、城市化水平和单位耕地面积农业机械总动力回归系数绝对值高值区还包括黑河市辖区和嫩江县，且由西（西北）向东逐渐递减。原因可能是西部地区经济较为发达，人均 GDP、城市化水平和路网密度都较高，一方面，在促进了城市扩张的同时带来了建设占用耕地面积的增加，即使有“耕地占补平衡”来保障耕地面积不减少，但补充耕地质量上的下降，终使得耕地产能下降；另一方面，路网密度的增加，导致人类活动对耕地利用的干扰程度加剧，增加了耕地系统安全的负面影响。经济的快速发展，也带来了高强度的农业现代化利用，过度的农业三废施用导致土壤板结，破坏土壤的水土保持能力和生产能力，更显著地致使耕地系统安全警情的升高。

5 结 论

1）黑龙江省耕地系统有 69个县域处于预警状态，处于轻警、中警和重警的县域分别有28、32和9个，其中预警最低值出现在牡丹江东宁县，最高值出现在大庆市辖区，仅有 3个县域处于较安全状态，整体警情水平较高，耕地系统安全状况不容乐观。

2）黑龙江省耕地系统安全预警空间差异性和聚集性特征明显。西部预警水平最高，东部次之，南北轴带地区（包括中部）最低。正、负相关类型在空间分布上呈现鲜明的对比特征，正相关类型以“组团”形式凸显，聚集性较强，高-高（H-H）型主要聚集在西部的齐齐哈尔市和大庆市县域，低-低（L-L）型则聚集在南北轴带上的中部区域，主要为伊春，哈尔滨北部和东部，以及牡丹江北部县域；低-高（L-H）型和高-低（H-L）型县域较少，且分别零星分散在高-高（H-H）型和低-低（L-L）型周围。

3）高程、坡度和年均气温3个自然因素和人均GDP、城镇化水平、地均固定资产投资、水土协调度、单位耕地面积农业三废指数和路网密度 6个社会经济因素对耕地系统安全预警值有着显著的影响，且影响的空间差异较为显著。

传统OLS模型通过最小二乘法对参数进行“全局”和“平均”估计，忽视了耕地系统安全预警根本驱动因素作用的空间异质性和不平稳性，导致结果难以准确解释各县域耕地系统安全预警驱动力的空间作用程度；而构建的GWR模型拟合估算结果精度较高，效果更优，充分考虑了各县域空间异质特征，基于GWR模型的耕地系统安全预警状态驱动因素的可视化，能够精准地刻画各因素的空间变化特征，可为研究区耕地系统保护提供新的思路和视角。

县域尺度下粮食主产区黑龙江省耕地系统安全预警状态空间分异特征，结合其驱动因素的影响系数空间分布，可以判断出研究区不同县域内耕地系统安全预警状态变化的主导因素及限制因素，分别识别耕地保护和警情亟待缓解的重点区域，为研究区制定分区缓解耕地系统安全警情的方案和措施提供了参考。

因考虑耕地质量作为耕地系统安全预警影响的直接影响因素，体现为耕地系统安全预警的直接表征（耕地系统安全预警评定指标），而本研究耕地系统安全预警指标体系的构建是基于抗原-抗体两方面考虑，因此将耕地质量以粮食生产稳定性指数、耕地标准系数、耕地农业产值综合性指标体现，且本研究从系统、宏观角度进行考虑耕地系统安全预警根本驱动因素，因而对于涉及耕地质量的微观指标，如有机质含量、地下水水位等没有纳入驱动因素进行分析。此外，本研究从定性角度选取耕地系统安全预警指标（直接影响因子）和耕地系统安全预警驱动因素（根本影响因子），如何进一步结合定量分析厘清二者之间的关系，并探究耕地系统安全预警影响机理，构建精准的耕地系统安全警情空间调控方案是下一步研究的方向和重点。

[1] 宋戈，连臣. 黑龙江省耕地资源安全预警分析及预警系统的构建[J]. 农业工程学报，2012，28(6)：247－252.Song Ge, Lian Chen. Analysis and system construction of safety early warning for cultivated land resources in Heilongjiang province[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(6): 247－252. (in Chinese with English abstract)

[2] 赵婷婷，姜广辉，张凤荣，等. 大城市近郊县域耕地保护预警方法[J]. 农业工程学报，2010，26(1)：335－339.Zhao Tingting, Jiang Guanghui, Zhang Fengrong, et al.Methods for early-warning of cultivated land protection in outskirts of megapolis[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2010, 26(1): 335－340. (in Chinese with English abstract)

[3] 蔡运龙，傅泽强，戴尔阜. 区域最小人均耕地面积与耕地资源调控[J]．地理学报，2002，57(2)：127－134.Cai Yunlong, Fu Zeqiang, Dai Erfu. The minimum area per capita of cultivated land and its implication for the optimization of land resource allocation[J]. Acta Geographica Sinica. 2002, 57(2):127－134. (in Chinese with English abstract)

[4] 黄贤金，曲福田. 耕地生态经济预警的理论与方法[J]. 生态经济，1998 (5)：14－17.Huang Xianjin, Qu Futian The theory and methods for early warning system of cultivated land eco-economy[J].Ecological Economy, 1998, (5):14－17. (in Chinese with English abstract)

[5] 刘友兆，马欣，徐茂. 耕地质量预警[J]. 中国土地科学，2003，73(12)：9－12.Liu Zhaoyou, Ma Xin, Xu Mao. Preliminary study on the early warning of cultivated land quality. China Land Science,2003,73(12): 9－12. (in Chinese with English abstract)

[6] 代路，雷国平，姜博，等. 基于生态经济系统的佳木斯市耕地质量预警分析[J]. 东北农业大学学报：社会科学版，2010，8(2)：6－9.Dai Lu, Lei Guoping, Jiang Bo, et al. The early-warning analysis of cultivated land quality of eco-economy in Jiamusi[J]. Journal of Northeast Agricultural University:Social Science Edition. 2010, 8(2): 6－9. (in Chinese with English abstract)

[7] Ye L, Ranst E V. Production scenarios and the effect of soil degradation on long-term food security in China[J]. Global Environmental Change, 2009, 19(4): 464－481.

[8] 刘艳芳. 多尺度耕地供需动态平衡预警体系研究[J]. 武汉大学学报：信息科学版，2004，29(5)：420－424.Liu Yanfang. Forewarning system of dynamic balance between cultivated land demands and supplies at multi-measures[J]. Geomatics and Information Science of Wuhan University. 2004, 29(5): 420－424. (in Chinese with English abstract)

[9] McBratney A, Field D J, Koch A. The dimensions of soil security[J]. Geoderma, 2014, 213: 203－213.

[10] Han B L, Liu H X, Wang R S. Urban ecological security assessment for cities in the Beijing–Tianjin–Hebei metropolitan region based on fuzzy and entropy methods[J].Ecological Modelling, 2015, 318: 217－225.

[11] 徐美，朱翔，刘春腊. 基于RBF的湖南省土地生态安全动态预警[J]. 地理学报，2012，67(10)：1411－1422.Xu Mei, Zhu Xiang, Liu Chunla. Early-warning of land ecological security in Hunan Province based on RBF. Acta Geographica Sinica. 2012, 67(10): 1411－1422. (in Chinese with English abstract)

[12] 陈藜藜，宋戈，邹朝晖.基于免疫机理的黑龙江省耕地系统安全预警研究[J]. 中国土地科学，2017，31(5)：71－81.Chen Lili, Song Ge, Zou Chaohui. Study on early-warning of cultivated land system security in Heilongjiang Province based on the immune mechanism[J]. China Land Science,2017, 31(5): 71－81. (in Chinese with English abstract)

[13] 宋戈，王越，雷国平. 松嫩高平原黑土区耕地利用系统安全影响因子作用机理研究：以黑龙江省巴彦县为例[J]. 自然资源学报，2014，29(01)：13－26.Song Ge, Wang Yue, Lei Guoping. Effect mechanism research of influential factors of cultivated land use system security of black soil region in Songnen High Plain: A case study of Bayan County in Heilongjiang Province. Journal of Natural Resources, 2014, 29(01): 13－26. (in Chinese with English abstract)

[14] 肖思思，吴春笃，储金宇. 近15年中国耕地驱动因素的空间差异分析[J]. 中国土地科学，2013，27(7)：55－61.Xiao Sisi, Wu Chundu, ChuJinYu. Spatial variation of driving factors on cultivated land change in China for the past 15 years[J]. China Land Sciences, 2013, 27(7): 55－61. (in Chinese with English abstract)

[15] 王粟，史风梅，裴占江，等. 松嫩平原农田土壤污染现状分析与评价：以黑龙江省绥化地区为例[J]. 东北农业大学学报，2015，46(05)：75－83.Wang Su, Shi Fengmei, Pei Zhanjiang, et al. Evaluation and analysis of farmland soils pollution status of Songnen Plain:For example Suihua area of Heilongjiang Province[J].Journal of Northeast Agricultural University, 2015, 46(5): 75－83. (in Chinese with English abstract)

[16] Chen J. Rapid urbanization in China: A real challenge to soil protection and food security[J]. Catena. 2007, 69: 1－15.

[17] 成琨，付强，任永泰，等. 基于熵权与云模型的黑龙江省水资源承载力评价[J]. 东北农业大学学报，2015，46(8)：75－80.Cheng Kun, Fu Qiang, Ren Yongtai, et al. Evaluation of bearing capacity of water resources in Heilongjiang Province based on entropy weight and cloud model[J]. Journal of Northeast Agricultural University, 2015, 46(8): 75－80. (in Chinese with English abstract)

[18] 张德华，娄厦，刘慧萍. 黑龙江省粮食生产地位及潜在问题研究[J]. 农业现代化研究，2014，35(5)：525—530.Zhang Dehua, Lou Sha, Liu Huiping. Analysis on the status and potential problems of grain production in Heilongjiang Province[J]. Research of Agricultural Modernization, 2014,35(5): 525－530. (in Chinese with English abstract)

[19] 聂英. 中国粮食安全的耕地贡献分析[J]. 经济学家，2015，1：83－93.

[20] 陆丽珍，詹远增，叶艳妹，等. 基于土地利用空间格局的区域生态系统健康评价：以舟山岛为例[J]. 生态学报，2010，1：245－252.Lu Lizhen, Zhan Yuanzeng, Ye Yanmei, et al. Regional ecosystem health assessment based on land use pattern: A case of study of Zhoushan island[J]. Acta Ecologica Sinica,2010, 1: 245－252. (in Chinese with English abstract)

[21] 陈藜藜，宋戈，邹朝晖，等. 基于改进局部空间自相关分析的萝北县耕地保护分区[J]. 资源科学, 2016, 38(10)：1871－1882.Chen Lili, Song Ge, Zou Chaohui, et al. Zoning of protected cultivated land based on improved local indicators of spatial association[J]. Resources Science, 2016, 38 (10): 1871－1882. (in Chinese with English abstract)

[22] 刘玉，郭丽英，刘彦随. 1980－2008年环渤海地区县域粮食的时空动态变化及分析[J]. 农业工程学报， 2012，28(12)：230－236.Liu Yu, Guo Liying, Liu Yansui. Spatial-temporal dynamic change and analysis of grain production in county scale in Bohai Rim from 1980 to 2008[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(12): 230－236. (in Chinese with English abstract)

[23] 陈强，朱慧敏，何溶，等. 基于地理加权回归模型评估土地利用对地表水质的影响[J]. 环境科学学报，2015，35(5)：1571－1580.Chen Qiang, Zhu Huimin, He Rong, et al. Evaluating the impacts of land use on surface water quality using geographically weighted regression [J]. Acta Scientiae Circumstantiae, 2015, 35(5): 1571－1580. (in Chinese with English abstract)

[24] Tu J, Xia Z G. Examining spatially varying relationships between land use and water quality using geographically weighted regression I: model design and evaluation[J].Science of the Total Environment, 2008, 407(1): 358.

[25] 吴莉，侯西勇，徐新良. 环渤海沿海区域耕地格局及影响因子分析[J]. 农业工程学报，2014，30(9)：1－10.Wu Li, Hou Xiyong, Xu Xinliang. Analysis of spatial pattern of farmland and its impacting factors in coastal zone of Circum Bohai [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014,30(9): 1－10. (in Chinese with English abstract)

[26] 张耀军，任正委. 基于地理加权回归的山区人口分布影响因素实证研究：以贵州省毕节地区为例[J]. 人口研究，2012，36(4)：53－63.Zhang Yaojun, Ren Zhengwei. Factors affecting population distribution in mountainous areas: geographically weighted regression using data from Bijie[J]. Population Research,2012, 36(4): 53－63. (in Chinese with English abstract)

[27] 陈晓红，王玉娟，万鲁河，等. 基于层次聚类分析东北地区生态农业区划研究[J]. 经济地理，2012，32(01)：137－140.Chen Xiaohong, Wang Yujuan, Wang Luhe. Study on the eco-agriculture regionalization based on hierarchical cluster analysis in Northeast China[J]. Economic geography, 2012,32(1): 137－140. (in Chinese with English abstract)

[28] Hu S, Yang S, Li W, et al. Spatially non-stationary relationships between urban residential land price and impact factors in Wuhan city, China[J]. Applied Geography, 2016,68: 48－56.

[29] German R N, Thompson C E, Benton T G. Relationships among multiple aspects of agriculture's environmental impact and productivity: A meta-analysis to guide sustainable agriculture[J]. Biological Reviews, 2017, 92: 716－738.

[30] Li Y X, Zhang W F, Ma L, et al. An analysis of China’s fertilizer policies: Impacts on the industry, food security, and the environment[J]. Journal of environmental quality, 2013,42, 972－981.