基于像元尺度耕地质量局部空间自相关的基本农田划定

刘彦文 刘成武 何宗宜 周 霞

(1.武汉大学资源与环境科学学院， 武汉 430079; 2.湖北科技学院资源环境科学与工程学院， 咸宁 437100; 3.中南民族大学公共管理学院， 武汉 430074)

0 引言

农田为人类粮食安全提供了保证，划定永久性基本农田将保障中国农业发展的生产基线[1]，我国农田资源的大量流失引起了各级政府部门的关注[2]，近年来划定了基本农田保护区、粮食生产功能区和重要农产品生产保护区，对优质农田进行永久性保护。基本农田划定是一个复杂的过程[3]，虽2017年完成了全国1.03亿hm2永久基本农田的划定工作，且较以往划定成果质量有明显提升[4]，但由于优质耕地选择标准较难统一，划定成果质量仍有改进空间。为此，需探讨相关划定方法，为实现数量、质量、生态“三位一体”耕地保护新格局，制定后续两区科学合理的划定标准，具有重要意义。

目前研究对象以基本农田、高标准基本农田、永久基本农田等为主；研究内容主要集中在从耕地质量、土地利用、空间布局等出发的农田保护、划定/分区、结果评价等方面；GIS是主要研究方法，也有层次分析法、局部空间自相关、LESA、TOPSIS、熵权法、理想点法等。在制定优质耕地入选为基本农田的标准时，基于耕地质量的评价体系占有较大比重。在评价指标体系中，耕地自然质量、区位条件最为常用[5-11]。有关学者考虑耕地空间布局与社会发展因素，引入空间形态、社会经济等相应指标对指标体系进行完善[12-14]。耕地兼有生态服务功能，生态环境与耕地稳定性指标也被相关学者关注[15-17]。在研究尺度方面，多以市县级为对象，从耕地图斑角度出发进行分析[18]，从栅格网与像元粒度的分析则更加细致[19-20]。鉴于耕地在空间上分布存在一定的规律性，文献[18,21-23]对基于局部空间自相关的方法进行了研究。

基本农田必定是耕地，但耕地不一定是基本农田，选择优质耕地的标准除了常用的耕地自然条件、区位条件等之外，可加入反映耕地利用情况、作物长势等情况的内因指标。本文以关联反映耕地利用状态与作物自身信息的NDVImax-min、NDVIsd作为完善耕地质量的综合评价指标。以湖北省嘉鱼县为研究对象，将所有资料归一化至30 m×30 m栅格像元尺度，在采用局部自相关分析像元耕地质量综合指标值的基础上，对基本农田进行划定，以期为进一步完善耕地优选为基本农田的指标体系、提升县级基本农田划定综合质量评价至像元级提供相关参考。

1 研究区与数据来源

1.1 研究区概况

嘉鱼县位于湖北省东南部、长江中游南岸(图1)，地理坐标为113°39′～114°22′E，29°48′～30°19′N。县域地处江汉冲积平原，属副热带湿润季风气候，日照充足，雨量充沛，土地肥沃，全县辖8个乡镇1个农场，人均土地0.28 hm2，人均耕地约0.09 hm2。地势自西南丘岗向东北平原交迭过渡，总体海拔较低，大部分在19～50 m之间，全县大体形成“一山三水四分田、两分道路与庄园”的地貌格局。

图1 研究区示意图Fig.1 Location of study area

1.2 数据来源与处理

数据主要包括国土部门提供的土地利用现状、土地利用规划、基本农田和土地利用分等定级等数据库资料；研究区域土地利用现状、市县乡村四级土地利用意见等调研资料；从地理空间数据云(http:∥www.gscloud.cn/)下载的研究区自然年内15期Landsat遥感影像，购买研究区GF1影像，研究区社会经济数据以及Google Earth影像数据；通过实地考察与调研进行资料分析判断与综合取舍，与已有相关数据、互联网收集资料等进行交叉验证，以检验数据来源的真实可靠性。将所有数据按照评价指标体系归类，在ArcGIS中将每个指标数据处理为一层栅格数据[24]。

2 研究方法

2.1 耕地质量综合评价

2.1.1评价单元与指标体系

县级是基本农田划定的落实单位，评价单元是划定成果构成的基本单元，其尺度与划定成果质量直接相关。在市县级尺度上的评价单元包括耕地图斑[18, 22]、不同尺度栅格单元[19,24]等，本文以30 m×30 m栅格为基本评价单元。

《基本农田划定技术规程(TD/T1032—2011)》是基本农田划定的基本依据，结合区域农用地分等成果、区域基本农田划定细则等，现有指标体系基本覆盖耕地自然质量、耕地立地条件、耕地利用、社会经济和生态环境条件等基本准则。以农田作物自身为参照，可以将现有指标体系归为外部因素，其反映了外界环境对划定成果的影响。所有的外部因素最终反映在农田自身方面，植被NDVI时序数据记录了不同时段作物自身的状况，文献[25-26]对NDVI时序数据与耕地利用状况进行了相关研究，同一耕地空间单元如在自然年内NDVI发生了变化，说明对应的土地被实际利用过，NDVI时序数据的统计变化特征可以辅助完善基本农田划定指标体系。本文选择自然年内像元NDVI时序数据中的最大值与最小值之差NDVImax-min间接反映耕地利用与作物自身的差异，用NDVI数据序列的方差NDVIsd间接反映耕地利用频率等的差异，并将其归入耕地景观生态条件，最后所建立的耕地质量综合评价指标体系如表1所示。

表1 嘉鱼县耕地质量综合评价指标体系Tab.1 Comprehensive quality evaluation system of cultivated land in Jiayu County

注：Pi=C表示将P指标对应变量i赋值为C，Ed指欧氏距离。

2.1.2评价指标赋值方法

对于距离连续型变量，首先运用ArcGIS中Euclidean Distance计算县域内每个像元到对应目标的距离，然后视距离指标与优选基本农田准则的正负相关性，将所有值都归一化在[0，100]以内。对于常规离散型指标，Q1、Q4赋值计算公式为

F+=100(Xi-Xmin)/(Xmax-Xmin)

(1)

F-=100(Xmax-Xi)/(Xmax-Xmin)

(2)

Q2、P1、P2参考《农用地质量分等规程》及已有文献如表1赋值，pH值的赋值计算公式为[27]

(3)

其中a1～a4分别为5.0、6.5、7.5和8.5，Abs(·)为取绝对值函数。

像元连片度计算，针对每个像元，运用Qeen空间关系邻接规则，判断其周边8个相邻方向上是否有基本农田像元，每发现一个相邻的基本农田像元，其像元连片度增加1。单个像元的连片度最小为0、最大为8，数值越大，像元连片度越好。优先保留原有基本农田中的高等级耕地、集中连片耕地，要求划定后基本农田集中连片程度有所提高[28]。以耕地图斑为基本单元进行研究划定时着重关注图斑的连片性，由于图斑的平均尺度一般相对像元(约900 m2)而言较大，划定精度不及像元尺度，所以本文以像元为基础，逐像元计算其连片程度，最后再运用式(1)将像元连片程度进行归一化。

2.1.3网络层次分析法确定指标权重

耕地质量综合评价指标部分相互独立、部分存在关联，网络层次分析法(ANP)不但具有层次分析法(AHP)递阶式层次结构，而且还具有内部依赖性和反馈性的层次结构[29]。采用ANP与专家打分法相结合的方式确定指标权重。依据ANP各指标权重求解过程，首先对指标关联情况进行分析，并构建判断矩阵；其次，将专家评分结果输入SuperDecision软件进行一致性检验，对不满足一致性检验的结果进行调整并反馈给专家再次征求意见；最后，输入评分结果、求解超级矩阵计算指标权重(表1)。

2.1.4像元耕地质量综合分值计算

像元是本文的基本评价单元，采用多因素综合评价法对评价单元进行综合质量评分，分值越高耕地质量越好[18]，像元综合分值计算公式为

(4)

式中Zi——第i个像元的耕地质量综合得分

wj——第j个指标的权重

Zij——第i个像元第j个指标的初始值

gstan(·)——标准化函数，保证像元在每个指标的得分都在[0,100]之间

2.2 局部空间自相关分析

近年来，很多学者基于耕地质量指数采用局部空间自相关方法对耕地合理利用进行了相关研究[21-23,30-31]。局部空间自相关指标(LISA)可以揭示空间参考单元与其邻近单元属性特征值之间的相似性或相关性，运用LISA指标中最常用的局部莫兰指数(Local Moran’sI)和Moran散点图对耕地质量在空间上的聚集、异质或随机的分布特征进行研究，局部莫兰指数计算式为

(5)

(6)

(7)

式中Ii——单元i的局部莫兰指数

wij——单元i与j之间的空间权重

n——与单元i相邻接的空间单元数量

xi、xj——单元i、j的耕地质量指数

σ2——方差

空间权重矩阵是进行局部莫兰指数计算的关键，不同邻接规则得到的空间权重矩阵不同，Rook、Queen是最常用的邻接规则。耕地的连片性在中心像元周边8个方向上均有体现，本文选择一阶Queen邻接规则构建空间权重矩阵。局部莫兰指数高值表明有相似变量值的面积单元在空间集聚，低值表明不相似变量值的面积单元在空间集聚，Moran散点图将空间单元的集聚形式分4个象限进行可视化。第1、3象限分别为HH与LL正相关型，各自表示高值与高值、低值与低值集聚，第2、4象限分别为LH、HL负相关型，空间异质性主要区别是高值包含低值异常和低值包含高值异常。

2.3 耕地利用与像元NDVI时序数据关联分析

NDVI时序数据在反映植被生长季节性和年内特征方面具有良好的可信性[32]，文献[25-26]利用其对耕地、休耕地进行了相关研究，并取得了良好效果。同一像元由于自然年内不同种值季数的差异(如一年一季、一年两季、一年多季等)会在NDVI时序数据中有不同取值，NDVI时序曲线峰值个数及其方差、最大最小值及其极差等都与耕地利用状况相关联。理论上曲线形态与种植收割相对应，裸地从作物种植开始到其生长最旺盛，从成熟到收获，对应曲线的不同取值，曲线方差越小说明耕地种植利用越平稳。极值和极差可反映一个种植收割期内不同作物间的横向生长状况。利用NDVImax-min、NDVIsd指标关联分析耕地利用状况，其计算公式为

NDVImax-min=max(NDVIseries)-min(NDVIseries)

(8)

NDVIsd=SD(NDVIseries)

(9)

式中NDVIseries——NDVI时序数据序列

其中，max(·)、min(·)、SD(·)分别为针对NDVI时序序列取最大值、最小值和方差函数，数据越密集效果越好。

3 结果与分析

3.1 耕地综合质量指标处理结果与分析

本文像元面积约900 m2，与耕地图斑相比其单元面积较小且均一，同时与研究区人均耕地866 m2也最为接近。基本农田的划定要求上图入库、落实到户，像元尺度可在更细粒度上落实基本农田的保护范围与保护责任。将评价耕地质量所需的16个指标在ArcGIS中各自处理成一幅30 m×30 m的栅格数据，统一采用WGS-1984-UTM-Zone-49N投影坐标，并采用Snap Raster保证不同指标层单元匹配对齐。全县范围内所有指标进行归一化处理，范围在[0,100]之间，结果如图2所示，运用式(4)可得各像元耕地质量综合指数。

图2 像元尺度耕地质量评价指标栅格图Fig.2 Raster maps of farmland quality evaluation indexes at pixel scale

由图2可知，数值离散型指标(如Q1、Q2)由于原始评价基于耕地图斑而赋值，其局部区域内取值较为集中，除了边界更替处对像元区分有贡献之外，同质区域内部对像元区分无贡献。对于数值连续型指标(如L1、E4)，以像元为单元计算距离，可以在图斑内部进一步细化距离权重对划定结果的影响。相比以图斑为单元计算图斑质心与目标之间的距离、图斑内部处处同值而言，基于像元的分析可以使划定结果更趋合理，本文有7个指标与距离相关，占总指标的43.75%。

3.2 局部空间自相关结果与分析

采用GeoDa与ArcGIS软件进行局部空间自相关分析与制图，首先将像元转为矢量点，并将像元耕地综合质量增加到点属性，其次将GeoDa计算结果保存到点属性中，最后将矢量点带属性转为与已有指标同样的栅格进行可视化分析。

全县耕地质量局部空间自相关指数为0.864 5，由于耕地像元与非耕地像元的耕地综合质量评价得分呈两级分化趋势，所以Moran散点图(图3a)呈两个团聚状，综合得分对耕地与非耕地具有良好的区分度。就耕地像元而言，其耕地质量局部空间自相关指数为0.991 6，其空间聚集程度较全县更加明显(图3b)。

图3 像元耕地质量空间关联分析Fig.3 Spatial correlation analysis of farmland quality in pixel scale

由图3c、3d可知，像元耕地质量空间关联集聚效果明显，在空间上，全县范围的HH型主要分布在东北平原区，包括牌洲湾镇、潘家湾镇、渡普镇，以及新街镇的大部分范围，该区交通便利、自然条件优越，是典型的水稻蔬菜基地[33]，LL型在长江、斧头湖等水域区最为明显，县域西南部以丘陵岗地为主，长江沿线南门湖村、官洲村等HH相对明显外，主体区内其余类型交替相间，集聚规律不明显。由图3d可知，正相关类型与耕地质量关联情况更为明显，尤其体现在交通便利、地势平坦、土质肥沃的省道S102两侧，但牌洲湾镇北部及其与潘家湾镇连接处相关性不明显，可能受洪涝灾害及近年来新起的畈湖工业园区的影响等所致。中部空白处是县城所在地，其西侧南门湖村是唯一保留的大面积耕地区，其HH型尤为明显。负相关类型集中体现在官桥镇、高铁岭镇以及陆溪镇的大部分耕地中，该区较东北部耕地分布较为零星、地势错落起伏、林地比重较大，像元耕地综合质量较低，LL型聚集与此相对应。

由表2可知，在全县范围内HH与LL正相关类型占72.5%，负相关类型共占1%，且其中90%主要分布在HH型周边、并与其形成了良好的连片分布。在耕地范围内，正相关类型占59.3%，说明耕地质量在空间上集聚特征明显。无显著相关的区域占比40.7%，P1、E3与Q1、Q3等因子在此处高低值交替相间，空间关联错综复杂。对耕地范围内全县集聚情况分析可知，HH型与耕地重合度为81.3%，说明即使在耕地范围未知的情况下，若优先利用HH型进行基本农田划定，也可达到较好的效果，这对于利用遥感反演信息客观进行基本划定提供了良好的依据。

表2 局部空间自相关结果分析Tab.2 Analysis of local spatial autocorrelation results

注：P<0.05，在95%置信度时统计结果；NS型表示非显著型。

3.3 基于像元局部相关性的基本农田划定结果与分析

局部空间自相关值Moran’sI反映了中心像元与周边像元的相关程度，其在[-1，1]之间，越靠近1表示中心像元属性与周边像元同质属性正相关性越强，反之，越靠近-1表示其负相关越强，等于0时表示规律性较弱，呈空间随机分布状态。依据其相关关系，相关学者按HH型、LL型、LH型和HL型从宏观层面定性对耕地利用进行了人为划分区块并分析建设时序的研究，如呈空间扩散效应的HH型为高质量耕地图斑被同样高质量图斑所包围，所以应优先划定、区域内禁止非农建设。本文从像元微观层面出发，利用相关属性值由大到小循环迭代依次选入单个像元，并将划定结果与空间集聚情况进行联合分析，具体实现过程利用Python完成。

全县划定任务为304 843个像元单位，由于在全县范围内土地利用类型复杂、质量评定指标多样，局部空间自相关莫兰指数理论上存在非耕地高于耕地的情况(如部分水域)，所以，选择以下两种方法进行划定，方法1是在全县范围内按前述计算的像元耕地综合质量得分，由高到低依次选入像元确定划定区域；方法2是在耕地范围内按莫兰指数由大到小顺序依次选入像元确定划定区域，结果如图4所示，对比分析如表3所示。由图4a可知，不管是利用全县范围内耕地质量还是耕地范围内莫兰指数，其划定结果都保持了高度的一致性，二者重合像元数为259 040个，占划定任务的85.0%。结合耕地质量评价数据与Google高分影像叠置分析，基于耕地质量的划定，其优势主要在于划入的像元是现有指标下耕地区内质量最优的区域，且与二者共有部分的连片性保持较好，基本在其周围分布，但其对丘陵区陇形地带等效果欠佳。基于莫兰指数的划定结果对县域西南部丘陵岗地区存在的条形、陇形相对碎片耕地的区域划定效果占优，但整体而言连片性不及前者。

由图4b与表3可知，两种方法都优先将空间上呈HH型聚集的99.97%像元划入了保护范围，这是划定成果的重点保护区域，区域内应该禁止非农建设[21]。HH型中未划入的29个像元虽然耕地质量较好，但大部分呈独立像元形式存在，小部分仅一边与HH相连，极大地影响了片块景观布局形态，在相关项目占用耕地时，可以考虑优先利用。两种方法的差异主要在于对LL型和NS型的处理上，基于耕地综合质量的划定优先将质量得分较高但指数得分低、在空间上呈NS型聚集的部分划入，其主要分布在牌洲湾镇北部、畈湖工业园周边、渡普镇与斧头湖连接处等。基于莫兰指数的划定则优先考虑了指数得分，与前者差异在于质量得分相对前者较低，在空间上呈LL型。此部分主要分布在陆溪镇西部、高铁岭镇东侧，后者明显可见丘陵岗地内陇形耕地轮廓特征。由于基本农田的划定是一个复杂多层次的过程[34]，二者差异部分也体现出了各自的优势，前者在片块连接性方面相对较好，后者在丘陵岗地区划定较为接近现实，二者空间上的直观分布，为最终决策提供了良好参考。

图4 划定结果Fig.4 Graph of demarcation results

类型LISA莫兰指数/耕地质量重合部分莫兰指数独立部分耕地质量独立部分像元像元比例/%像元比例/%像元比例/%比例/%HH型11515311512432.90029032.9LL型926454695813.44566113.026026.4LH型7800007800HL型100000010000NS型1426709695827.714204557013.040.7小计35064625904074.04580313.04580313.0100

4 结论

(1)在县级空间尺度内耕地综合质量涉及距离类指标时，一般基于像元的评价较以耕地图斑为单元的评价更具良好的区分性，划入的基本农田理论上更趋合理，对人均、户均基本农田保护责任的落实也更加精确。但基于像元的相关分析涉及数据量较大，对运算资源的要求也相对较高。

(2)湖北省嘉鱼县耕地综合质量在空间上表现出较强的空间相关性，全域和耕地范围的局部空间自相关指数分别为0.864 5和0.991 6。基于基本农田包含于耕地的先验知识，运用后者的划定更趋合理，其HH正相关类型主要集中在潘家湾镇、渡普镇等东北平原主导区，LL负相关类型主要集中在官桥镇、高铁岭镇等西南部丘陵地带，且分布较为零星。

(3)基于耕地质量与莫兰指数的划定结果在空间上保持了高度的一致性，两种方法都优先将空间上呈HH型聚集的99.97%的像元划入了保护区，前者在综合质量方面占优，但对丘陵区陇形地带等划定效果不及后者。基于像元的分析对于嘉鱼县的划定更为合理，NDVI等指标间接关联了土地利用状态与农作物自身信息，对耕地质量综合评价体系的完善有益。