耕地质量由耕地所处的自然与环境条件的状况构成，耕地质量直接影响作物产量和农产品质量^[1]，其对保障粮食安全、生态安全和人类生存安全具有重要的意义。而随着我国经济的快速发展和人口的日益增长，耕地质量面临降低的风险：一是因城市扩张而实施的耕地占补平衡，在实践过程中出现了异化现象，如占优补劣、占水田补旱地等，虽然耕地数量平衡，但耕地质量不平衡，导致大量优质耕地流失；二是因人为活动而导致耕地环境变化，出现水土流失、环境污染等问题，从而导致耕地环境质量下降^[2-4]。因此，保护高质量的耕地并治理退化严重的耕地势在必行。了解并掌握耕地质量状况是耕地质量保护与提升工作的基础，而高效、精确地监测耕地质量变化是其突破点，也日益受到研究学者的关注。耕地质量监测点布设是监测工作开展的前提和基础^[5]，通常采取抽样调查的方法，通过设置一定数量的监测点来监测耕地质量等别变化，以点带面获取区域耕地质量的信息^[6]。耕地质量具有变异性，需要大量的监测点来表征，然而由于人力物力、时间成本等因素的限制，布设大规模的监测点并不现实^[7]，因此，合理布设样点对节省外业调查、内业分析成本和提高监测的准确性具有重要的现实意义。

目前，监测点布设常用方法是经典抽样的概率模型，包括简单随机抽样、分层随机抽样、系统抽样等，此类方法假设样本间相互独立，忽略空间相关性。简单随机抽样布设样点呈不均衡分布，系统和分层抽样效率相对较高，但预测有偏差^[8]。这些方法只考虑了样点被选中的概率，而没有考虑样点的空间相关性，并不能很好地适用于具有空间相关性的耕地质量属性值监测^[6]。为此，祝锦霞等^[9]考虑了空间相关性，采用地统计学方法建立变异函数模型，依据变异函数的变程确定网格大小，基于等别组合优化样点布设。王倩等^[7]对比分析平原地区和丘陵地区的耕地质量的空间变异情况，依据两者不同变异情况，以变程的一半为间隔进行监测点布设。胡晓涛等^[5]依据变异函数的变程网格布点，考虑土地利用规划、自然等的面积比例、土壤类型等修正样点。杨建宇等^[6]对比不同方形网格布设的精度和样点数量，得到恰当网格尺寸的基础监测网，随后采用泰森多边形迭代加入待添加点，直至克里格误差不再减小。

已有的考虑空间相关性的耕地质量监测点布设研究常用规则网格采样法，该方法网格的间距由耕地质量的空间变异情况决定，网格内耕地质量性质相对均一^[6]，在每个网格中的耕地地块上布设样点，根据样点信息利用克里格插值方法估计区域的耕地质量状况。因此，采样网格的设计影响样点的分布及数量，根据耕地质量的空间变异情况来合理优化网格设计，可以提高监测效率、节省监测成本。但现有的网格法只考虑了空间距离，而耕地质量不仅在空间上有自相关性，还与其周边的环境因子表现出相关性^[10]。在土壤属性监测领域已有研究利用协同克里格来考虑与目标变量相关性较好的环境因子信息作为辅助变量，在同等预测精度下可以减少采样点数量。已有学者利用辅助因子，采用协同克里格对土壤铜^[11]、全氮^[12]、锌^[13]、剖面电导率等^[14]土壤属性进行采样数量的优化研究，在满足采样精度的条件下，得到的合理样点数相比普通克里格法减少了10%~40%。例如，李艳等^[15]的研究表明，不论目标变量的样品数目如何减少，利用辅助变量的协同克里格法和回归克里格法的预测精度较普通克里格法均有较大提升。李俊颖等^[16]的研究表明相比于普通克里格法，协同克里格法利用与目标变量相关性较好的辅助变量，可以更合理地确定采样数量。上述学者均利用协同克里格法对土壤养分、盐分等土壤属性进行采样数量优化研究，但使用协同克里格法对耕地质量进行监测网格优化的研究鲜有报道。

鉴于此，本研究以广东省广州市增城区为研究区，以国家自然质量等指数为研究对象，结合耕地质量分等因素，通过相关性分析筛选易于获取的辅助变量，运用地统计学方法分析耕地质量的空间分布特征，结合协同克里格法优化耕地质量监测样点布设，对快速准确掌握耕地质量的变化、指导耕地质量保护和提升具有重要意义。

1 材料与方法 1.1 研究区概况

广州市增城区位于广东省中东部（图 1），地处23°05′~23°37′N、113°32′~114°00′E，总面积1 616.47 km²。全区地势北高南低，境内山地以低山为主，山脉整体上呈东北-西南走向，气候类型是亚热带海洋性季风气候，年平均气温12.1~28.5 ℃，耕地土壤以赤红壤和渗育型水稻土为主。截止到2016年末，增城区耕地面积共229.76 km²，其中水浇地76.56 km²、旱地7.21 km²、水田145.99 km²。耕地呈现均匀分散分布。为了更好地保护增城区的耕地质量，监测是必要手段，而形成高精度、低成本的监测样点方案是难点。因此本研究选取增城区作为研究区进行耕地质量监测点优化布设研究。

1.2 数据来源

本研究采用广州市增城区2016年度耕地质量等别更新评价成果，样点布设主要考虑了耕地自然等别，将国家自然质量等指数作为研究对象，原因如下：①国家自然等指数评定考虑了地形地貌、水文、土壤和基本设施情况，能有效反映耕地自然状况和农业设施综合情况；②耕地利用等指数和经济等指数是在自然质量等指数的基础上，通过计算土地利用系数和经济系数并进行修正得到，两者均与自然等指数有较强的相关性。

1.3 研究方法

首先利用相关性分析筛选辅助变量，运用ArcGIS进行不同规格网格布设，进而获得监测点，采用基于辅助变量信息的协同克里格法和普通克里格法进行插值模拟，比较分析两者不同网格尺度下的监测效果，进而得到优化后的样点布设。

1.3.1 辅助变量选取

以国家自然质量等指数为研究对象，计算其与环境因子的相关系数，筛选出其中相关性较强的因子作为辅助变量，以协助进行协同克里格插值。

耕地质量不仅受到自然要素、工程要素以及社会因素的直接或间接影响，而且还会受到三者两两组合特征以及全组合特征的影响^[17]。结合研究区实际情况，耕地质量变化主要受到自然条件，如有机质、pH值、地下水位等因素的制约，同时根据辅助变量的易获取性对因素进行筛选^[18-19]。因而，为定量分析耕地质量的影响因子，以国家自然质量等指数为对象，选取了耕地质量分等因素：有效土层厚度、有机质、pH值、土壤质地、剖面构型、地下水位、排水条件、灌溉保证率、地形坡度、田面坡度。以上变量数据均采用增城区耕地分等数据库中的标准化分值，该数据库中所采用的指标大部分来源于《农用地质量分等规程》（GB/T 28407—2012）（以下简称规程）中的指标体系。

1.3.2 网格设置

在ArcGIS中设置不同规格的网格，形成不同规格的网格中心点，分别与耕地矢量数据进行空间叠加，取落在网格中心点的耕地图斑作为初步监测点。对于落在耕地图斑以外的中心点，在初步监测点的基础上，通过ArcGIS近邻分析工具选取距离中心点最近的耕地图斑作为监测点，形成监测点预布设。

1.3.3 插值模拟

应用GS+工具模拟半变异函数模型，得到最优模型。根据最优模型的参数在ArcGIS中进行协同克里格和普通克里格插值模拟。

（1）普通克里格

在地统计分析中，区域化变量具有结构性和随机性特点，一般采用半变异函数来反映。同时，其空间变异情况与空间距离间的关系也可形象体现^[20]。在满足二阶平稳性的前提下，半变异函数定义为：

(1)

式中：r (h)为半方差；N（h）表示两两样点间隔为h的样本量有N (h) 对，h为样点间的距离；Z (Xi) 和Z (X_i + h)分别表示变量Z（X）在X_i和X_i+h处的实测值。

（2）协同克里格

在满足本征平稳假设，同时存在2个协同区域化变量的情况下，交互半方差函数定义为：

(2)

式中：r_ij (h)为2个协同变量的交互半方差；N (h)为样点间隔为h的成对量，h为样点间的距离；Z_i (X_a)和Z_i (X_a + h)分别为变量Z_i (X)在位置X_a和X_a + h处的实测值；Z_j (X_a)和Z_j (X_a + h)同理。

若某一属性经过拟合变异函数和相关性分析，得出该属性存在空间相关性，便可运用普通克里格进行插值^[21]。计算公式为：

(3)

式中：z′(x₀)代表预测点x₀的预测值；z (x_i)则代表实测值；λ_i代表各实测值的权重；n代表参与x₀点预测的z的实测点有n个。

协同克里格是普通克里格的一种扩展，从单一利用待预测点的空间自相关性，转变到既利用前者又利用多个变量间的相关性^[14]。此方法需要有2种变量，一是主变量，二是辅助变量（可为多个）。计算公式为：

(4)

式中：z₁是主变量；z₂是辅助变量；z₁ (x_i)和z₂ (x_j)分别为各变量的实测值；λ_1i与λ_2j分别是主变量和辅助变量的各实测值的权重，且∑λ_1i=∑λ_2j=1；n和p分别是参与预测的z₁和z₂的实测点数量。

1.3.4 方法比较

以自然质量等指数为对象，选择均方根误差及相关系数2个精度比较指标，采用交叉验证的方式，分析比较不同网格尺度下结合辅助变量的协同克里格法和普通克里格法的监测效果，进而得到优化的监测点布设。

1.4 数据处理

自然因素量化分级的标准是其对作物生长发育的适宜程度。对于数据库中的标准化分值数据，若是数值型和定序型，将继续沿用；若存储形式为文本的定性型，如成土母质和剖面构型，将根据两者种植作物的适宜性在SPSS 22.0实现量化分级；正态分布性检验、相关性分析通过SPSS 22.0实现；半方差函数分析、交互半方差分析在GS+9.0中实现；网格布设样点、插值绘制分布图采用ArcGIS 10.2实现。

2 结果与分析 2.1 辅助变量的选取与处理

利用相关性分析初步分析国家自然质量等指数与环境因子间的关联关系，根据相关系数来判断关系强弱。由于排水条件、灌溉保证率等环境因子属于定序数据，本研究采用Spearman等级相关系数进行相关性分析。定序数据是一种可排序但不可做四则运算的数据，仅反映对象间的等级、顺序关系。根据规程，国家自然质量等指数值域是1~6 000，值越大，耕地质量越好；排水条件和灌溉保证率值域划分为4级，级越低，排水条件和灌溉条件越佳；地形坡度和田面坡度分别划分为6级和5级，级越高，坡度越大。结果见表 1。

从表 1可知国家自然质量等指数与各环境因子的关联程度：地下水位最强，其次为灌溉保证率和pH值，之后按排水条件、有效土层厚度、土壤质地、地形坡度、田面坡度、有机质及剖面构型依次排列，且都达到了极显著水平，说明国家自然质量等指数与上述环境因子在研究区受同样的区域化现象或空间过程的影响，属于协同区域化变量^[22]。其中国家自然质量等指数与地下水位相关性最强，且呈负相关，说明增城区地下水位越高，耕地质量越差。综上所述，将地下水位确定为辅助变量。

2.2 空间变异结构分析

区域化变量是指在空间上分布的变量，具有2个基本性质，一是空间相关性，二是空间依赖性。在此基础上，地统计学以半方差函数为主要工具，可定量地描述这两大性质，并最大限度地解释空间变异信息。而耕地质量正是符合这2大性质的区域化变量^[23-24]。因此，有必要对其进行空间变异结构分析，模拟出一个最优半方差函数的空间理论模型，定量地描述耕地质量的空间相关性和依赖性，并从变量背景和半方差函数理论模型的角度进行科学分析和解释。一般情况下，半方差函数的计算要求数据符合正态分布，否则影响结果的有效性^[25]。经S-W检验可知，P= 0.865 06>0.05，服从正态分布。

表 2给出了国家自然质量等指数拟合半变异函数结果，运用GS+，以不同半方差模型对国家自然质量等指数进行拟合，以最大的拟合系数和最小残差为评判标准，得到国家自然质量等指数最佳拟合模型是球面模型，块金值是3 500，说明存在实验误差和田间采样等人为因素造成的空间变异。根据CAMBARDELLA等^[26]提出的空间相关度划分原则，分为三个层次：第一层次，块金系数小于25%，表示变量呈现强空间相关性；第二层次，块金系数处于25%~75%之间，表示变量呈现中等空间相关性；第三层次，块金系数大于75%，表示变量呈现较弱空间相关性。因此，本研究中块金系数等于5.24% 表明国家自然质量等指数的空间变异特性受到了自然特征和人为活动的共同作用，表现为强空间相关性。

2.3 合理样点布设分析

耕地质量变化监测通常采取抽样调查的方式。为了获取详细的耕地质量变化信息，需要充足的采样数据，但其具有高成本性和耗时性，进行大规模的采样并不现实。因此，样点布设的合理最小化成为一个日益重要的研究内容，这对耕地质量动态变化的长期监测研究具有重要意义。许多研究表明，协同克里格在充分利用前期数据的基础上，可减少样本数，优化样点布设，主要是利用更易获取的变量作为辅助变量协助插值以提高主变量的插值精度。本研究采用交叉验证的方式，以均方根误差和相关系数为判断标准，分别代表预测值与实测值间的误差和两者的相关性，以此来验证插值精度。RMSE值越小，表明预测值越接近实测值。r越大，则表明预测误差越小，精度就越高。

2.3.1 网格布设样点

本研究以国家自然质量等指数为目标变量，采用网格法布设样点，通过协同克里格来优化样点布设。在研究区按规则选取样点，样点的位置处于规则网格的中心或是最邻近中心点的耕地图斑。分别以2 km×2 km至6 km×6 km 5种网格为研究对象，以1 km为梯度分别对增城区耕地布设监测样点，布设的样点数分别为390、176、104、66个和47个（图 2）。基于原始数据中各等别样点数占总样点数量的比例考虑以上样点是否合理分配到各个等别。分析在不同网格布设下利用地下水位为辅助变量的协同克里格插值精度，进行耕地质量监测样点布设优化研究。

2.3.2 半方差分析

由表 3分析可知，不同网格布设下国家自然质量等指数的半方差和交互半方差函数模型拟合的块金系数具有强或中等空间相关性。单变量表现出强烈的空间自相关，主要受土壤母质、地形、气候等结构因素影响；交互变量表现出中等自相关，不仅受到结构因素影响，还受到施肥、管理、耕作等随机因素作用。两者的差异可能是因为：相比单变量，交互变量加入了辅助变量——地下水位的信息构建模型，而地下水位随着城市化进程加快，人类活动越加频繁，受到随机因素影响程度加深，所以影响了构建模型的块金系数，进而得到交互变量表现出中等自相关的结果。从变程这一角度分析，较单变量而言，交互变量的变程更长，表明在单变量的基础上，交互变量由于融合了辅助变量提供的更丰富更多元空间信息，扩展了变量具有空间相关性和依赖性的范围。因此，在相同网格布设下，与普通克里格法相比，协同克里格法拥有更广的适用范围。

2.3.3 克里格插值精度比较

对比5种网格布设下协同克里格插值的精度，主要选取均方根误差及预测值与实测值间的相关系数作为比较指标。由表 4可知，5种网格布设下，与普通克里格法相比较而言，协同克里格法预测的均方根误差分别减少了8.83%、11.23%、10.35%、0.11% 及9.94%，而预测值与实测值的相关系数分别提高了146.93%、79.37%、80.09%、136.36% 及574.19%。表明相同的网格布设下，协同克里格法由于融合了更多更丰富的空间信息，预测的精度高于普通克里格法，可利用该方法优化监测点布设。当网格尺度加大至5 km×5 km时，协同克里格法的预测精度仍高于2 km×2 km网格布设下的普通克里格法。因此，在保证监测精度的条件下，本研究区耕地质量监测点布设利用协同克里格法优化研究，网格尺度可加大至5 km× 5 km，样点数量可减少至66个。

为了更直观地比较不同规格网格布设下，不同监测点数量对掌握耕地质量空间变异情况的影响，分别用普通克里格法和协同克里格法绘制了不同监测点数量下的耕地质量的空间变异分布图（图 3）。随着网格规格的扩大、监测点数量的减少，插值结果出现平滑效应，说明随着监测点数量的减少，预测精度降低，会丢失一些空间信息。对比普通克里格和协同克里格插值结果，可以看出协同克里格法相比普通克里格法可以获得更丰富的空间信息，预测效果更优，说明协同克里格法结合辅助变量的信息可弥补目标变量在实测点和待测点之间位置关系的缺失。

3 讨论

协同克里格法和普通克里格法均可反映耕地质量的空间变异分布特征和规律。在相同网格布设下，相比普通克里格插值法，以地下水位为辅助变量的耕地质量协同克里格插值法融合了更丰富多元的空间信息，监测效果更优。

耕地质量监测点布设目的是通过监测点推断全区的耕地质量等别变化，难点在于形成高精度、低成本的样点方案。传统的监测点布设较少考虑耕地质量与周边环境因子的协同变化关系，本研究提出的方法利用了单一环境因素作为辅助变量的信息，在5 km×5 km网格尺度下，样点数更少的协同克里格插值精度仍高于2 km×2 km网格布设普通克里格插值。因此，可利用基于辅助变量的协同克里格法优化监测点布设，以获取更少监测点，这样不仅降低成本，而且保证了监测精度。

综上所述，本研究提出的方法对于优化县域尺度耕地质量监测点布设具有一定的适用性，可为其他类似地区提供借鉴。但本方法还存在一定的局限性，例如仅考虑了单一环境因素作为辅助变量，可在今后的工作中进一步开展经济因素等其他因素作为辅助变量或多个因素作为辅助变量的研究。

4 结论

将广州市增城区作为研究区，以国家自然质量等指数为研究对象，通过相关性分析和地统计分析，以及克里格插值和基于地下水位辅助变量的协同克里格插值进行监测点布设优化研究。主要研究结论如下：

（1）在相同网格布设下，结合地下水位辅助变量信息可以显著提高预测结果精度，协同克里格法预测耕地质量等别结果精度明显优于普通克里格法。相比普通克里格法，协同克里格法适用性也更强。

（2）在5 km×5 km网格尺度下，协同克里格法利用少量样点得到耕地质量等别预测精度高的样点布设方案，且优于普通克里格法的布设方案。协同克里格利用辅助变量信息可弥补网格规格变大、样点数量减少带来的信息缺失，从而优化网格布设，实现少量样点的布设方案。

以广州市增城区为实例的研究结果表明，在5 km×5 km网格尺度下，最终监测点数量为66时，本研究提出的基于协同克里格法的县域耕地质量监测点优化布设方案在耕地等别预测方面的精度高于2 km×2 km的普通克里格法，能有效提高监测效率，预测县域耕地等别并监控其变化情况，满足县域耕地质量等别变化监测的需求。