精准农业是现代农业发展的方向, 长期以来备受关注。全面分析掌握农田土壤特征, 精准划分管理区, 是采取差异化管理措施, 实现精准农业的基础^[1-2]。近年来, 一些学者针对精准管理分区方法展开研究, 探讨了运用决策树、系统聚类、K均值聚类(灰色聚类)和模糊聚类等方法^[3-6]划分管理分区的可行性, 其中K均值聚类和模糊c-均值聚类(FCM, Fuzzy c-means algorithm)^[7-8]应用较多。由于很多分类对象间没有明确的分类界线、存在亦此亦彼的表现, 因此相比而言, 模糊c-均值聚类作为一种新的空间连续性聚类算法, 克服了K均值聚类只能将每个对象严格地分类到某一聚类中心、与其他聚类中心的空间相关关系得不到体现的缺陷^[9], 对于复杂的土壤管理问题更为适宜。基于FCM划分管理分区的关键在于聚类分区数目c和模糊加权指数φ的确定, 已有研究主要采取经验取值法、c－φ多次组合最优法、指数检验法等^[10-12], 对分区结果进行检验主要分为基于内部变量和外部变量两种方法^{[7, 13]}。

综观现有农田精准管理分区研究, 小尺度的试验田研究较多、大中尺度的农田研究少, 土壤条件较好的内陆田块试验较多、滨海盐渍化区域缺乏涉及；只针对单一土壤属性的研究较多, 综合考虑土壤多个属性的研究较少；从分区结果检验来看, 基于单方面变量的检验都存在较明显的误差^[14]。因此, 本文拟以山东省无棣县农田为例, 针对滨海盐渍化农田盐碱瘠薄、土壤属性空间变异大、粗放管理效益低的现实, 在县域尺度上, 研究基于多个土壤属性变量的农田管理分区精准划分方法和分区精度校验方法, 以提高分区结果的合理性和科学性, 为滨海盐渍化地区采取差异化管理措施, 提升盐渍化土地利用水平提供科学依据^[15]。

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

无棣县(37°41′—38°18′N, 117°31′—118°12′E)地处山东省最北部, 总面积1963 km², 东北部濒临渤海, 属北温带东亚季风区域大陆性气候, 年均降水570.1 mm, 80%集中在夏季, 秋冬春三季干旱, 年均蒸发1285.5 mm, 蒸降比2.3。近年来, 该县作为滨海低平原区新增耕地潜力和增产潜力较大的县区之一, 置身新增千亿斤粮食生产能力规划、盐碱地治理、“渤海粮仓”科技示范工程等国家战略叠加区, 成为农业开发热点区域^[16]。但是, 无棣县淡水匮乏, 地下水位浅、矿化度高, 土壤瘠薄盐碱问题突出^[17], 农业生产效率和效益较低, 如何通过精准分区、差异化管理提高农田利用水平和生产效益成为亟需解决的问题。由于县境东北部主要为盐田、海水养殖区和荒草地, 因此本文研究区仅涉及农田分布区, 总面积约为756 km²。

1.2 数据来源

在土地利用现状图上以2 km×2 km网格覆盖研究区农田(避开建设用地、盐田、滩涂), 面积不足者并入相临网格, 每个网格中心设计1个调查样点(图 1), 野外调查时根据现场情况确定实际采样地点, 使用手持式GPS记录坐标, 并用五点取样法采取对应的土壤样品1 kg左右装袋密封, 带回实验室化验。野外调查在2016年5月24—28日进行, 共采土样193个。土样经自然风干、研磨、过筛, 化验分析^[18]获得193组土壤有机质、有效氮、有效磷、速效钾、含盐量和pH等属性数据。

1.3 研究方法 1.3.1 土壤属性特征与变异性分析

运用SPSS 20.0进行土壤属性描述性统计分析, 并通过Person-双侧检验得到相关系数矩阵, 分析土壤各属性变量的变异性及相关性。在ArcGIS 10.2中运用地统计学方法, 分析土壤属性的空间分布特征。土壤养分等级参照全国第二次土壤普查分级标准划分为极低、低、中、丰、高5级, 含量区间从低至高依次为, 有机质(g/kg)：(0, 6)、[6, 10)、[10, 20)、[20, 30)、[30, ∞)；有效氮(mg/kg)：(0, 30)、[30, 60)、[60, 90)、[90, 120)、[120, ∞)；有效磷(mg/kg)：(0, 5)、[5, 15)、[15, 30)、[30, 40)、[40, ∞)；速效钾(mg/kg)：(0, 40)、[40, 75)、[75, 120)、[120, 150)、[150, ∞)；含盐量(g/kg)结合作物耐盐能力分为脱盐化、轻度、中度、重度和极重度5级：(0, 1)、[1, 2)、[2, 3)、[3, 4)、[4, ∞)；pH分为酸性、适宜、中等碱性、强碱性和极强碱性5级：(0, 6.5)、[6.5, 7.5)、[7.5, 8.5)、[8.5, 9)、[9, ∞)^[17]。

1.3.2 模糊均值聚类

模糊c-均值聚类(FCM)是一种非监督的连续性动态聚类算法, 定义一个n×p的土壤属性多源数据集(n=193, p= 6), 最常用的目标函数为：

(1)

式中, μ_{i j}为i个样本的土壤属性对于第j个聚类类别的隶属度值, φ(φ≥1)是模糊加权指数(Fuzziness exponent), d_ij²是第i个样点的土壤属性与第j个聚类中心之间距离的平方。

聚类有效性检验通过模糊性能指数(FPI, Fuzziness performance index)和归一化分类熵(NCE, Normalized classification entropy)评价。FPI和NCE的取值范围在0—1之间, 取值越接近0, 分类界线越明显^[19], 反之分类界线越模糊。公式为：

(2)

(3)

式中, H(U; c)是分类熵的函数表达式, a可以取任意正整数。

1.3.3 管理分区划分及可视化

运用MAPGIS 6.7进行研究区底图及采样点绘制, 并将得到的隶属度赋属性到各采样点中, 实现管理分区的可视化表达。由于在MAPGIS中只能依据最大隶属度原则进行分类, 不能体现模糊连续分类在空间预测中的重要作用, 因此运用ArcGIS地统计分析模块对土壤属性模糊隶属度进行插值分析, 实现其空间连续性表达。

1.3.4 分区精度检验

为了评价基于模糊聚类分析进行管理分区的精度, 确定其变异程度是否满足管理分区的要求, 对不同管理分区内的土壤属性进行变异性分析, 并用最小极差法(LSR, Least significant range method)对各分区结果进行差异显著性检验, 校验管理分区结果的精度。

2 结果与分析 2.1 土壤属性特征及其变异性

在SPSS 20.0支持下, 采用邻域法对获取的193组土壤属性数据进行特异值处理^[20], 并进行描述性统计分析(表 1)。由表 1可见, 研究区农田土壤有机质含量均值19.253 g/kg, 处于中等水平；有效氮均值59.349 mg/kg, 处于较低水平；有效磷均值29.527 mg/kg, 处于中等水平；速效钾均值179.676 mg/kg, 含量较高；含盐量均值1.204 g/kg, 总体为轻度, 这与采样时已进入夏季、土壤水分含量较高有关；pH均值8.184, 为中度碱性, 土壤盐渍化程度总体呈轻至中度水平。有机质、有效氮、有效磷、速效钾和含盐量均呈中等变异性(C_V < 10%为弱变异性, 10% ≤ C_V < 100%为中等变异性, 100% ≤ C_V为强变异性), pH值总体较高但变异性较弱, 是否作为管理分区的依据有待于分析与其他土壤属性的相关性。综合来看, 研究区土壤属性总体呈中等变异, 有必要进行分区管理。

2.2 土壤属性相关性分析

在SPSS 20.0中对土壤属性变量进行相关性分析, 得到相关系数矩阵(表 2)。由表 2可得, 有机质、有效氮、有效磷、速效钾和pH值等5个指标的相关性较高。含盐量虽相关性较低, 但由于研究区位于环渤海低平原区, 地下水埋深浅且矿化度高, 土壤盐渍化风险大, 有必要将含盐量作为分区评价指标之一。因此, 确定将上述6个土壤属性作为管理分区的评价指标。

2.3 土壤属性空间变异特征

在SPSS 20.0中用单样本K-S(Kolmogorov-Smironov)工具进行正态分布检验, 选择默认显著性为0.05、置信区间为95%, 若数据的渐进显著性P值>0.05, 则假设成立, 符合正态分布(表 3)。由表 3可见, 有机质、有效氮、pH的K-S值均大于0.05, 为正态分布；有效磷、速效钾和含盐量为偏态分布, 经对数变换后K-S值大于0.05, 达到显著水平, 呈正态分布, 可以进行克里格插值分析。

运用ArcGIS 10.2地统计工具进行克里格插值, 采用C-V交叉检验法, 通过多次试验选择最优拟合模型^[21]。最终确定有机质、速效钾的拟合模型均为指数模型, 有效磷、有效氮、含盐量和pH的拟合模型为球面模型。此时, 土壤属性各拟合模型的标准平均值最接近于0, 标准均方根最接近于1(表 4)。

一般来说, 块金效应(Nugget / Sill)的取值小于25%时, 具有强空间相关性；在25%—75%之间时, 具有中等空间相关性；大于75%时, 具有弱空间相关性。由表 4可见, 土壤有机质、有效氮和有效磷的块金效应值分别为51.9%、67.4%、67.0%, 说明这三种土壤属性的空间相关性为中等水平, 空间变异特征受随机因素和结构因素的共同影响。速效钾、含盐量的块金效应值分别为20.8%、22.4%, 都小于25%, 说明含盐量和速效钾具有强空间相关性, 其变异特征受随机因素影响较小, 主要由结构性因素引起。pH的块金效应值为0, 说明pH具有完全的空间自相关性, 其空间变异特征完全由结构性因素引起。

运用克里格法确定土壤属性最佳拟合模型后, 生成土壤属性的空间插值图(图 2), 以对土壤属性空间变异规律进行最优预测, 实现土壤属性空间变异的数量化及可视化。由图 2可见, 有机质高值区主要分布在研究区西南和西北部, 呈块状分布。有效氮和有效磷高值区的分布趋势较为一致, 中部和南部分散分布, 含量最高或最低处呈小块状或点状分布。速效钾的分布具有明显的地带性特征, 高值区主要分布在中东部地区, 呈东高西低现象。含盐量在中部地区较低, 且分布范围广, 含量约为1.0—2.0 g/kg, 属于轻度盐渍化, 少数点状分布地区含盐量在2.0—4.0 g/kg之间, 属于中重度盐渍化。pH值呈块状分布, 值高(8.5—9.0)的地方主要分布在中东部地区。

2.4 管理分区划分及其可视化

在MATLAB R2016a中通过“z-score”数据标准化处理, 得到的n×p(n=193, p=6)的土壤属性模糊数据集。调用FCM函数, 取幂指数为1.6, 最大迭代次数为200, 收敛阈值0.000001, 进行模糊聚类分析。为了确定最佳分类数c, 选取分类数2—8分别进行模糊聚类分析, 并根据模糊性能指数(FPI)和归一化分类熵(NCE)进行评价(图 3), 两者的值越小, 聚类效果越明显。由图 3可以看出, 当分类数c取值为3时, FPI和NCE的值最小, 说明此时同一分类中的数据差异性最小, 不同分区间的差异性最大, 确定本研究区的最佳分类数为3个。

由模糊聚类分析得出的隶属度矩阵, 依据最大隶属度原则将采样点分为3类, 在MAPGIS 6.7中根据各采样点的最大隶属度划分管理分区, 实现管理分区的空间可视化(图 4)。由于最大隶属度原则是以网格为单位进行分类, 不能完全体现模糊分类的空间预测性。本文根据分类后各采样点的模糊隶属度矩阵, 在ArcGIS 10.2中利用地统计模块进行插值分析, 实现土壤属性模糊隶属度的空间连续性表达(图 5)。

对比图 5和图 2可得, 分区Ⅲ主要分布于有机质含量丰等级、有效氮含量中和丰等级、有效磷含量高等级区域；分区Ⅰ的中西部主要分布于有机质含量中等级区域, 东南部地区主要分布于有效氮含量丰等级区域；分区Ⅱ主要分布于速效钾含量高等级区域, 其他养分含量较低, 可以考虑多施钾肥除外的其他肥料。且盐渍化严重的北部地区和pH强碱性区域主要为分区Ⅰ和分区Ⅱ, 与分区Ⅲ土壤的拟合度较低。因此, 管理分区和土壤各属性的空间分布具有很强的空间相关性和空间拟合度, 各分区内土壤属性向均一化方向发展。对比图 4和图 5可得, 各分区土壤模糊隶属度空间预测图与最大隶属度原则划分的管理分区图在空间分布上基本一致, 各分区间的隶属度较为明显, 交叉重叠程度低, 说明以最大隶属度原则划分的管理分区具有空间预测的有效性。

统计各分区土壤属性平均值(表 5)可见, 分区Ⅲ中的有机质、有效氮、有效磷含量均最高, 速效钾含量较高, 土壤养分水平总体最高, 且含盐量、pH值较低；分区Ⅱ中的速效钾含量最高, 有机质、有效氮的含量在分区Ⅲ和分区Ⅰ之间, 有效磷比分区Ⅰ略低, 土壤养分水平总体居中；分区Ⅰ的土壤养分水平总体最低, 含盐量和pH值也最高。

2.5 管理分区精度校验

为了对基于模糊聚类分析进行管理分区的结果进行精度验证, 对不同分区内土壤属性进行统计分析, 得到基于变异系数和最小极差法(LSR)的各分区差异显著性检验结果(表 5)。

对比表 1和表 5可见, 各分区的土壤属性变异系数较全研究区均有所下降。其中, 有机质变异系数由25.0%下降为14.7%—23.9%, 有效氮变异系数由30.0%下降为17.9%—30.0%, 有效磷变异系数由52.3%下降为33.0%—51.5%, 速效钾变异系数由32.8%下降为18.2%—29.3%。pH值变异系数由3.0%下降为2.7%—3.0%, 变化较小, 均属于弱变异性, 分区结果对pH影响较小；含盐量变异系数由48.4%变为41.0%—54.1%, 分区Ⅲ变异系数下降, 分区Ⅰ和分区Ⅱ均略有升高, 说明含盐量分布不均匀。

LSR检验结果表明(表 5), 土壤有机质、有效氮、速效钾在3个分区均达到极显著差异(P < 0.01)或显著差异(P < 0.05), 有效磷含量在分区Ⅲ和分区Ⅱ之间具有显著差异(P < 0.01), 含盐量在分区Ⅱ和分区Ⅰ、分区Ⅲ间具有极显著差异(P < 0.05), pH值在分区Ⅲ和分区Ⅰ、分区Ⅱ之间具有极显著差异(P < 0.01)。

总体而言, 分区内的土壤属性变异程度较分区前全区的变异程度有所下降, 分区间的差异显著性明显, 分区结果具有可行性。其中分区Ⅲ的变异程度下降最明显, 各采样点土壤属性隶属度最明确, 分区Ⅰ下降程度最低, 分区Ⅱ居中。基于变异系数和LSR的检验结果均表明, 管理分区具有较高的精度。各分区内部的土壤养分变异系数由25%—52.3%减小为14.7%—51.5%, 各分区内部的空间变异性较全研究区减小, 且各分区间具有显著性差异, 分区结果可以作为农田分区调控的统一作业单元。根据图 5统计各分区面积, Ⅰ、Ⅱ、Ⅲ管理区面积分别为2.56万hm²、1.76万hm²、3.24万hm²。

3 结论

该文以无棣县农田为研究区, 采用网格法结合土地利用现状定点野外采样、室内化验分析获取土壤属性数据, 运用ArcGIS 10.2地统计方法分析土壤属性的空间变异特征；在MATLAB R2016a中采用模糊c-均值聚类法(FCM)计算各样点的模糊隶属度, 通过插值预测模糊隶属度的空间分布, 基于最大隶属度原则进行分区, 并根据模糊隶属度分布图分析分区结果的有效性；通过变异性分析和最小极差法(LSR)差异显著性检验, 对分区结果进行精度验证。

研究结果表明：无棣县农田土壤总体呈轻至中度盐渍化, 有效氮含量偏低, 有机质、有效磷含量中等, 速效钾含量较高；有机质、有效氮、有效磷、速效钾和含盐量呈中等变异性(变异系数25%—52.3%), 空间变异性较大, 应分区调控；空间变异分析显示, 速效钾、含盐量和pH的块金效应值都小于25%, 受人为因素影响相对较小, 主要受土壤质地、地下水矿化度等结构因素影响, 有机质、有效氮和有效磷的块金效应值在50%—75%之间, 受耕作方式、施肥等人为因素影响较大；将全县农田划分为Ⅰ、Ⅱ、Ⅲ共3类管理区, 估算面积分别为2.56万hm²、1.76万hm²、3.24万hm²；各分区土壤养分变异系数分别为23.9%—51.5%、15.9%—50.3%、14.7%—33.0%, 检验结果表明各分区间差异显著, 而各分区内部变异性明显低于全研究区。管理分区与土壤属性空间分布特征具有较高的拟合度, 表明分区结果可以作为差异化管理的作业单元。研究结果为各分区内部统一、不同分区间差异化管理提供了依据, 研究有助于推进滨海盐渍化农田精准化管理水平的提高。

由于分区后pH和含盐量的变异系数变化不大, 少数区域甚至出现含盐量变异性增大的现象, 考虑到当地受滨海盐渍化影响严重, 土壤水盐状况短时间难以改变, 可在上述管理分区的基础上, 对变化不明显或者变异性增大的区域进行细化管理, 以更好地改善其水盐状况。