基于降维方法、聚类算法与众数滤波窗口尺寸的精准土壤采样管理分区质量研究

《Smart Agricultural Technology》：Effects of Dimensionality-Reduction Approach, Clustering Algorithm, and Modal Filter Window Size on Management Zone Quality for Precision Soil Sampling

　　

由环境协变量（environmental covariates）划分的管理分区（management zones，MZs）可支撑特定场地管理（site-specific management）与精准土壤采样，但椒盐效应（salt-and-pepper effect）会降低其空间连续性与操作可行性。众数滤波（modal filtering）可减轻这种破碎化，但目前尚不清楚更平滑的图件是否同时改善区内同质性，以及该响应是否取决于窗口尺寸、聚类算法（clustering algorithm）与降维方法（dimensionality-reduction approach）。本研究评估了众数滤波窗口尺寸是否影响采用不同聚类算法与多元降维方法所生成MZs的质量。在两块农田中，由环境协变量采用常规主成分分析（PCA）或MULTISPATI-PCA结合K-means或Ward聚类生成MZs，并在五个逐步增大的众数滤波窗口（W3–W11）处理前后进行评估。分区质量采用土壤肥力属性的区内方差降低率（VR%）与平滑度指数（Smoothness Index，SI）评估。众数滤波一致提高了空间连续性，但其对区内同质性的影响因配置而异。K-means在较大窗口（W9–W11）下趋于呈现较高的VR%值，而Ward在较小或中等窗口（W3–W7）下趋于呈现较高值。在Field 2中，PCA–K-means配合W11对六个所评估土壤属性中的五个产生最高VR%。多数情景下，常规PCA呈现的实测VR%值高于MULTISPATI-PCA，而MULTISPATI-PCA生成的图件更为连续，但并不总是更同质。上述结果表明，空间连续性的提高并不必然意味着土壤属性的更大同质性，且在两块所评估农田中，滤波表现取决于降维方法、聚类算法与窗口尺寸的组合。

研究背景：农业土壤具有强烈的时空变异性，特定场地管理（site-specific management）旨在通过合理配置投入品以最大化作物生产力。管理分区（management zones，MZs）是具有相对均一的生物、土壤、气候与地形因素、并导致共同作物生产力限制的子区域，其划分可将田块细分为区内方差低于全区方差的区组方案，因此MZs内部降低目标属性变异的能力是其用于土壤采样与特定场地管理的重要方面，也是精准农业中经济可行的途径。MZ分割常用无监督机器学习技术，主要基于环境协变量（environmental covariates）、作物产量数据、土壤理化属性或上述数据组合，数据源选择取决于目标变量。MZ划分中需审慎选择输入变量，因为分区有效性取决于变量表征目标属性变异的能力；且诸多协变量高度相关，会引入冗余信息而影响聚类分析，故需采用变量选择或变换技术，其中主成分分析（principal component analysis，PCA）最为常用。PCA通过将相关的原始变量变换为互不相关的综合变量即主成分（principal components，PCs）实现降维，可按Kaiser与Jolliffe提出的准则依据主成分对数据方差的解释贡献确定保留的PCs数量，随后对所选PCs采用不同算法进行聚类以划分管理分区。可用于MZ划分的无监督算法中，K-means与Ward层次聚类较为突出，此外还包括模糊聚类及显式纳入空间约束或关系的算法。尽管K-means与Ward常用于分割多变量数据，但二者在属性空间中划分观测的方式不同：K-means基于到质心的距离直接划分数据，Ward采用凝聚层次方法以最小化组内方差。然而二者在常规形式下均未在聚类过程中显式纳入空间自相关或像元间邻域关系，因此分类主要由属性间的多变量相似性决定，而不必然由类别在田间的空间连续性决定。由此，分割后出现的主要问题之一即为椒盐噪声（salt-and-pepper noise），表现为一类小簇或孤立像元插入另一类之中，该格局增加空间破碎化并可能限制田间操作可行性，因为小而孤立的斑块可能需频繁调整肥料、改良剂或其他投入品的施用。促进分区空间连续性并减轻椒盐效应的替代方案是采用MULTISPATI-PCA，该多元降维技术由Dray等首次实现并由Córdoba等适配于管理分区概念。MULTISPATI-PCA是PCA的扩展，考虑原始变量中的空间自相关以生成综合变量，属于空间约束型PCA，其成分构建同时兼顾多元方差与空间自相关；常规PCA识别数据中的主要变异轴，而MULTISPATI-PCA优先考虑同时具有空间结构的变异轴，有利于生成受局部噪声影响较小、更为连续的管理分区。Córdoba等与Gavioli等均在MZ划分中以MULTISPATI-PCA替代PCA作为聚类前预处理而获得空间上更连续的分区。然而，更大的空间连续性并不必然意味着分区内目标属性的更大同质性；尽管MULTISPATI-PCA可能有利于更连续的分区，但其是否同时降低分区内土壤属性变异仍不清楚，且由MULTISPATI-PCA生成的更连续分区仍可能需要某种后处理以应对椒盐噪声效应并使MZ在操作上更易管理。降低椒盐效应的一种策略是对最终管理分区图应用众数滤波（modal filter），该后处理已用于修正分区边界并减少可能妨碍MZ操作可行性的小孤立斑块。例如Schepers等应用9×9像元移动窗口的众数滤波去除分区内小于1.25 ha区域所代表的噪声；Oldoni等采用15×15像元滤波以平衡空间噪声削减与原始分区边界保持，Melo等应用9×9像元窗口；Córdoba等亦强调采用3×3、5×5或7×7像元窗口的众数滤波平滑分区图。尽管该方法被反复使用，滤波窗口尺寸仍常被主观确定。该选择可能直接影响所得分区，因为小窗口对空间破碎化的作用可能有限，而较大窗口可能过度平滑类别边界、消除小斑块并合并本应保持区别的异质区域。例如Betzek等指出5×5像元滤波在分区操作性与维持空间异质性之间提供了较好平衡，但该研究仅比较了3×3与5×5像元两种移动窗口。因此，仍需理解增强滤波强度如何同时影响分区的空间连续性及其降低目标属性内部变异的能力；此外，由于不同聚类算法在滤波前产生不同的空间构型，对窗口尺寸增大的响应可能因算法而异。鉴于此，仍需以整合方式理解MZ划分中所采用的方法学决策如何影响其对众数滤波后处理的响应。

研究内容与结论：本研究从两个互补视角界定管理分区质量，即由区内变异降低所表达的区内土壤属性同质性，以及所得分区的空间连续性；二者并不必然等价，因为图件可能变得更平滑而区内土壤属性同质性并未相应提高。本研究假设众数滤波窗口尺寸对管理分区质量的影响取决于所采用的划分方法，因而评估了作为后处理步骤应用于不同聚类算法与多元降维技术的众数滤波窗口尺寸变化是否影响用于土壤采样的管理分区质量。为此，研究评估了众数滤波对由PCA、MULTISPATI-PCA、K-means与Ward层次聚类等多元降维技术与聚类算法不同组合所生成分区的影响。研究框架包括三主要步骤：环境协变量输入数据的预处理、采用两种替代方法的降维、以及采用两种聚类算法的分区分割，合计每块田产生四种不同分区划分方法；随后对每种方法生成的分区施以3、5、7、9、11像元参数的众数滤波，并在原始未滤波图件与不同滤波窗口所得图件上评估分区质量。研究选取巴西圣保罗州37 ha的Field 1（以极黏质 dystrophic Red Latosol 为主，按USDA土壤系统分类为Oxisol，海拔533–566 m）与美国佐治亚州30 ha的Field 2（以Tifton系列土壤为主，按USDA土壤系统分类为细壤质、高岭石型、热性 Plinthic Kandiudults（Ultisols），海拔79.2–91.1 m），以代表热带与温带气候条件下具对比性地形条件的 prevalent 农业系统。两田块分别按40×40 m与30×30 m规则网格进行密集土壤采样，对应密度分别约6与9.5 samples ha?1；Field 1每采样点由距中心5 m半径内六个子样组成混合样、以配备自动土壤探头的全地形车采集，Field 2由距中心6 m半径内约十二个子样组成、采用手动土壤采样器采集；采样深度分别为0–0.20 m与0–0.15 m，为各自区域谷物生产区推荐施肥深度。样品送土壤实验室常规分析，评价pH、阳离子交换量（cation exchange capacity，CEC）、Ca、Mg、K与P有效性。研究结论为：众数滤波通过减少椒盐噪声改善管理分区的空间连续性，但其对区内同质性的影响取决于所评估的配置，且该响应随所用算法而变化。对K-means，最佳结果趋于出现在较大的移动窗口，尤其W9与W11（约0.636与0.950 ha），表明本研究评估的K-means配置受益于更强的后处理以降低空间破碎化；对在所评估田块中产生破碎化较少分区的Ward方法，最大增益出现在较小或中等窗口，尤其W3、W5与W7，分别对应约0.071、0.196与0.385 ha；进一步增加滤波强度时，平滑可能促进异质区域合并并降低潜在增益。因此，众数滤波是优化管理分区的重要步骤，但不存在同时适用于两种算法的单一窗口尺寸，其效益取决于降维方法、聚类算法与窗口尺寸的组合，故其应用不应任意确定，而应以所采用配置保持或提高分区内目标属性同质性的能力为依据。另一点是，不存在在所有田块中均优越的单一算法，因为K-means与Ward的表现随所评估区域而变化；但对两种算法而言，当目标是划分采样分区时，采用常规PCA进行降维倾向于产生更好的区内同质性。换言之，在MULTISPATI-PCA中显式纳入空间成分倾向于平滑分区边界，这可能降低其表征化学属性中更频繁出现的土壤突变的?力。因此，在采样管理分区时应对空间平滑审慎解读，因为更连续的图件并不总是意味着分区内更大的属性同质性。

关键技术方法：研究在两块田布设密集规则网格土壤采样作为参考数据集（Field 1数据可见于Melo & Amaral，Field 2数据可见于Beasley并可向作者索取），以pH、CEC、Ca、Mg、K与P为评价土壤属性。采用三个环境协变量，即表观电导率（apparent electrical conductivity，ECa）（Field 1用EM38-MK2电磁感应传感器测至37.5 cm深、测线间距30 m、约150 readings ha?1；Field 2用Veris MSP 3100电阻率传感器测至30 cm深、测线间距12 m、约400 readings ha?1）、来自Copernicus DEM的数字高程模型（DEM，原始分辨率30 m）与三季作物产量均值图；产量图在AgDataBox软件中过滤去除异常值后以普通克里金逐幅插值至5 m并逐像元算术平均，ECa亦以普通克里金插值并转为5 m栅格。协变量在QGIS 3.34 Prizren的Precision Zones插件2.2.0中重采样至5 m双线性插值并共定位、行z-score标准化，随后并行采用常规PCA与MULTISPATI-PCA降维；后者在R中以ade4包的dudi.pca()与adespatial包的multispati()实现，邻域矩阵基于像元欧氏距离、以距离≤70 m为邻居，阈值依据Córdoba等与Oldoni等已发表应用确定，经spdep包dnearneigh()与nb2listw()转换为行标准化空间权重矩阵。两方法均仅保留特征值大于1的第一成分（PC1或CS1）。聚类采用K-means（插件中以K-means++实现，n_init=10、random_state=0、max_iter=300）与Ward层次聚类（R中hclust()，欧氏距离与ward.D2凝聚法），每田固定三个采样分区。众数滤波在Precision Zones插件中采用SAGA GIS实现的众数滤波，W3–W11分别对应7×7至23×23像元窗口、圆形邻域，边界外或NoData像元从滤波邻域中排除；5 m分辨率下对应邻域半径约15、25、35、45与55 m，面积约0.071、0.196、0.385、0.636与0.950 ha。分区质量以区内方差降低率VR%与平滑度指数SI评估，并以ΔVR、ΔSI（百分点）表示相对未滤波图件的最大变化；并以Spearman秩相关分析土壤属性与协变量关系、以Queen邻接空间权重矩阵计算Moran指数表征空间自相关。

研究结果：

3. Results and Discussion（结果与讨论）：滤波窗口尺寸在不同土壤属性上使管理分区表现最优的取值随各田块算法类型而异（Figure 4、Figure 5）。尽管Field 1与Field 2的VR%存在比例差异，但出现收敛格局：每种算法的VR%响应由所用众数滤波窗口尺寸驱动。两田块中Ward层次聚类倾向于在较小或中等窗口（W3–W7）呈现最高VR%值，而K-means倾向于偏好较大窗口，主要为W9与W11。该格局在Field 2采用PCA时尤为明显，K-means评估的六个属性中有五个在W11呈现最高VR%；该配置下相对原始图的增益达CEC 6.45、P 3.50、K 3.65、Mg 2.63与Ca 5.21个百分点（Table 5）。对Field 1（Figure 4），Ward一般对具有更强空间结构的土壤属性如CEC、K、Mg与Ca呈现较高VR%值（Figure 6），表明其形成内部同质分区的能力更强；但窗口尺寸的过度增大降低了表现，尤其在与常规PCA关联时，提示过度平滑可能合并原本不同的分区并消除部分相关空间异质性。另一方面，对K-means尤其Field 2（Figure 5），窗口尺寸的逐步增大对多数属性包括同样具有较强空间结构的CEC、P、K、Mg与Ca促进了VR%增益（Figure 6），表明较大窗口在降低椒盐空间噪声与修正更破碎分区方面有效。因此，众数滤波被证明是管理分区后处理中的相关步骤，但其应用应审慎：在所评估田块中，K-means生成的配置倾向于呈现更大空间破碎化，较大窗口一般改善空间连续性并提高VR%，而基于Ward的配置倾向于更为空间连续并在较小窗口下表现更好，因为过大的窗口可能人为改善分区周界的平滑度但降低分区表征土壤真实空间变异的能力。因此，滤波窗口的选择不应仅基于分区破碎化程度，而应主要考虑聚类方法及其在特定数据集中产生的空间构型，因为这些因素影响MZs对增强滤波强度的响应。

无论采用何种聚类算法，对滤波的响应在所评估组合中并非一致为正。多数情景下至少有一个滤波窗口相对原始图提高VR%，但某些窗口未优于未滤波图（Table 5）。该行为在Field 2部分基于MULTISPATI-PCA的组合尤其与K-means组合时更明显，其中CEC、Mg、Ca与pH在原始图中的VR%值高于任何滤波版本。然而，若干研究不对管理分区应用此类后处理滤波，或仅依据田块大小或固定窗口尺寸确定滤波窗口，而未采用定量准则评估其对区内同质性的影响。结果表明，滤波在采用适当窗口时可提高VR%，但在某些配置下也可能降低表现（Figure 4、Figure 5；Table 5），表明仅基于任意或视觉准则应用滤波可能非但未改善、反而损害区划的最终质量。

另一点是，降维技术影响了分区的初始平滑度、其对众数滤波的响应以及由此产生的VR%值。总体而言，多数所评估情景中PCA呈现的实测VR%值高于MULTISPATI-PCA（Figure 4、Figure 5）。相反，MULTISPATI-PCA初始生成更连续的图件（Figure 7、Figure 8），表现为略高的平滑度指数值（Table 6）。每个田块内PCA与MULTISPATI-PCA第一成分所代表的方差相近，Field 1分别为66.70%与65.96%，Field 2为56.47%与55.59%（Table 3）。然而，这些成分空间结构差异的幅度在田块间变化。Field 1中常规PCA的PC1的Moran's I为0.958，而MULTISPATI-PCA的CS1为0.960（Table 3），两方法对ECa、DEM与平均产量的载荷幅度相近（Table 4）。但在Field 2中，Moran's I由PC1的0.821增至CS1的0.850（Table 3），伴随DEM贡献增大（载荷幅度由0.666增至0.761）与ECa贡献减小（由0.441降至0.328）（Table 4）。该行为与在多元降维步骤中显式纳入空间结构相一致，因为MULTISPATI-PCA平衡方差与空间自相关，有利于空间组织性更强的成分。结果亦表明该纳入在Field 1中作用有限，因该田常规PCA第一成分已呈现高空间自相关，而在Field 2中更为明显。然而，更大的空间平滑并不必然意味着对土壤属性变异的更好表征。若仅依据SI解读，MULTISPATI-PCA倾向于被偏好；但土壤化学属性可能在短距离内突变，故过度平滑可能掩盖真实空间反差。在此意义上，较大比例所评估组合中PCA观察到较高VR%值表明，MULTISPATI-PCA所促进的空间组织性增强并不必然导致分区内土壤属性的更大同质性，尤其在与适当滤波窗口尺寸组合时。

SI在两田块中均随众数滤波的应用而增大（Table 6），证实滤波在减少椒盐效应与提高管理分区空间连续性方面的效率。然而，SI增益应结合图件的初始平滑水平来解读。Field 1中初始SI值已较高，表明原始图件破碎化低（Figure 7），因此滤波后的增量在数值上较小并主要反映平滑度的精细化。该田中各组合ΔSI仅为0.1至0.6个百分点，最大增幅见于PCA–K-means。但在Field 2中初始SI值较低，尤其基于PCA的组合，表明图件破碎化更高且后处理后改善空间更大（Figure 8）。因此，增幅更为明显，PCA–K-means与PCA–Ward分别达3.3与2.6个百分点，而基于MULTISPATI-PCA的组合ΔSI为1.6个百分点（Table 6）。因此，尽管某些SI绝对差异看似较小，其可能代表空间上相关的改变，尤其当与去除孤立像元、减少不规则边界与增加分区邻接性相关联时。该行为与关于管理分区图件修正的研究一致，其中邻域滤波即使在对已具高初始连续性的图件数值SI增益有限时也提高平滑度并降低破碎化。因此，SI应作为VR%的补充指标解读：VR%反映分区相对目标属性的内部同质性，而SI描述最终图件的空间组织与操作可行性。因此，当SI增大未伴随VR%增大时，仅凭更大的空间连续性不应被视为土壤采样表现更好的充分指标，从而强化VR%作为主要评价准则。

总体而言，Field 2观察到的VR%值低于Field 1，无论对K-means聚类还是Ward层次聚类均如此。无论采用何种算法，良好的区划结果取决于用于分区划分的协变量能否表征目标属性的变异。Field 1中观察到ECa与CEC（r = 0.74）、Ca（r = 0.71）、Mg（r = 0.69）与K（r = 0.59）之间以及DEM与CEC（r = ?0.56）、Ca（r = ?0.54）、Mg（r = ?0.62）与K（r = ?0.50）之间特别强的关系（Figure 9）。但在Field 2中，ECa与这些相同属性之间的关系较弱，为0.12至0.20，尽管产量与Ca（r = 0.59）、CEC（r = 0.48）与K（r = 0.45）呈中等关联。该反差与Field 1中观察到的较高VR%值相一致，提示协变量表征目标属性的能力可能因田块而异。本研究中由于土壤数据已可获得，故能测定该关系；但在实际使用条件下，农民与顾问在确定区划前并无此信息。在这些情况下，对属性与协变量之间是否存在关系的不确定性可能是该方法的重要局限，因为其倾向于降低区划增益，无论采用何种算法。一种替代方案是采用两步法，首先进行初步低密度采样以表征土壤属性的变异并评估其与可用协变量的关系。另一种可能性适用于已有历史土壤采样数据的田块，尤其此前采用网格采样生成土壤肥力图的田块。在此类情况下，现有数据可用于评估哪些协变量最能表征目标属性的变异，且若确定一致的关系，未来采样活动可从基于网格的策略转向基于分区的策略。基于该评估，将可能选择更适于区划的变量，从而降低初始不确定性并提高通过管理分区划分获得有效增益的可能性。

本研究所得结果亦提示未来研究进一步评估区划方法的机会。未来研究可考察较本研究所评估更强的滤波强度，采用大于W11的窗口，以确定增强空间平滑在损害区内同质性前仍在多大程度上具有益处。MULTISPATI-PCA中不同的邻域配置亦可评估，以评价分区对该方法学选择的敏感性。此外，可采用其他聚类算法包括模糊方法与显式纳入空间约束或关系的方法考察对众数滤波的响应，以评估K-means与Ward所观察到的格局是否亦出现于具不同分割特征的方法。将该方法应用于更多田块，涵盖不同空间变异水平、土壤类型、田块大小与种植制度，将可评估本研究所观察格局的一致性。其他环境协变量亦可纳入分区划分。例如Almeida等将卫星衍生植被指数与土壤电导率结合使用，而Gon?alves等展示了Sentinel-1时间序列数据用于管理分区划分。其他来源如裸土影像及与土壤属性与作物响应相关的其他指标亦可研究。未来研究还可通过考虑例如所需样本数与分区内土壤属性估计相关误差，直接评估这些配置对采样效率的影响。

讨论与结论总结：讨论部分围绕三条主线展开。其一，滤波窗口的最优取值依赖于聚类算法：K-means因常规形式未显式纳入空间邻域关系而易产生更破碎的分区，较大窗口（W9–W11）一般改善空间连续性并提高VR%，而Ward为基于最小化组内方差的凝聚层次方法，在所评估田块中产生破碎化较少的分区，其最大增益见于较小或中等窗口（W3–W7），进一步增大滤波强度可能合并异质区域而降低潜在增益；因此滤波窗口不应仅依据分区破碎化程度任意确定，而应主要考虑聚类方法及其在特定数据集中的空间构型。其二，滤波并非一致为正，部分MULTISPATI-PCA尤其结合K-means的组合中出现了未优于原始图的窗口，表明仅凭任意或视觉准则应用滤波可能损害区划最终质量；降维技术影响分区的初始平滑度、对滤波的响应及VR%值，常规PCA在多数情景下呈现较高实测VR%，而MULTISPATI-PCA初始生成更连续的图件但并未必然提高区内同质性，Field 1中常规PCA第一成分已具高空间自相关故纳入空间结构作用有限、Field 2中作用更明显，提示空间平滑不应被解读为对土壤属性变异表征的改善。其三，SI随滤波一致增大，但其增益需结合初始平滑水平解读，在初始SI已较高的Field 1中增量有限而在Field 2中更明显；SI仅描述空间组织与操作可行性，当SI增大未伴随VR%增大时不应视为采样表现更好的充分证据。此外，Field 2的VR%整体低于Field 1，与协变量和目标属性关系的强弱差异一致，提示协变量表征目标属性能力的不确定性是实际应用中的重要局限，并提出低密度预采样或利用既有网格采样历史数据的两步式思路，以及在未来研究中考察更大窗口、不同MULTISPATI-PCA邻域配置、其他聚类算法、更多田块与其他协变量并直接评估采样效率等方向。研究结论部分可译述为：本研究将K-means与Ward算法应用于环境协变量以生成用于土壤采样的管理分区，采用常规PCA与MULTISPATI-PCA作为降维方法；分区划分后，评估了移动窗口众数滤波相对原始图在减少椒盐噪声方面的效果，最终图件质量以土壤属性的VR%评估并辅以应用于分区的平滑度指数。总体而言，众数滤波通过减少椒盐噪声改善管理分区的空间连续性，但其对区内同质性的影响取决于所评估配置，且该响应随所用算法而异；对K-means最佳结果趋于出现在较大移动窗口尤其W9与W11（约0.636与0.950 ha），表明本研究评估的K-means配置受益于更强的后处理以降低空间破碎化，而对在所评估田块中产生破碎化较少分区的Ward方法最大增益出现在较小或中等窗口尤其W3、W5与W7（分别约0.071、0.196与0.385 ha），进一步增加滤波强度时平滑可能促进异质区域合并并降低潜在增益；因此众数滤波是优化管理分区的重要步骤，但不存在同时适用于两种算法的单一窗口尺寸，其效益取决于降维方法、聚类算法与窗口尺寸的组合，其应用不应任意确定，而应以所采用配置保持或提高分区内目标属性同质性的能力为依据；另一点是不存在在所有田块中均优越的单一算法，K-means与Ward的表现随所评估区域而变化，但对两种算法而言当目标是划分采样分区时采用常规PCA降维倾向于产生更好的区内同质性，换言之在MULTISPATI-PCA中显式纳入空间成分倾向于平滑分区边界，可能降低其表征化学属性中更频繁出现的土壤突变的能力，因此在采样管理分区时应审慎解读空间平滑，因为更连续的图件并不总是意味着分区内更大的属性同质性。