评估使用显式空间模型和非空间模型，结合Sentinel 1卫星数据及相关协变量来估算土壤颗粒大小分布和土壤质地的方法

《Soil Security》：Evaluating the use of explicitly spatial models and non-spatial models to estimate soil particle-size fractions and soil texture using Sentinel 1 and covariates

　　桑德拉·克里斯蒂娜·德奥多罗 | 拉斐尔·德安德拉德·莫拉尔 | 里蒙·费利 | 蒂姆·麦卡锡 | 罗文·费利
爱尔兰气候分析与研究单位（ICARUS），梅努斯大学地理系，基尔代尔郡，爱尔兰

**摘要**
地统计和非地统计方法已被应用于土壤科学中，处理从野外采集的数据到近感和遥感数据的各种类型。然而，关于使用微波传感器估算土壤质地的研究通常仅限于裸土条件。在这里，我们测试了雷达遥感数据和技术（后向散射系数、雷达植被指数和极化分解），并采用四种模型（非空间方法：线性回归和基于树的回归；空间方法：经验贝叶斯克里金回归-EBKR和共克里金-CK）来估算表土层（10厘米深度）的土壤颗粒大小分数和土壤质地类别。研究区域覆盖了爱尔兰的西海岸到东海岸。模型的性能通过留一法交叉验证方法进行了评估。结果显示，CK方法对沙子、粉砂和粘土的预测效果良好，Nash Sutcliffe效率值分别为0.9、0.6和0.5。然而，EBKR方法在捕捉插值表面的趋势和空间模式方面表现更好。重要的是，这项工作提供了一个方法论框架，为设计密集的实地土壤调查提供了依据（模型中不确定性较大的区域可能需要详细的土壤调查）。它还提供了一种在大范围内估算土壤属性的方法，并生成新的数据产品，用于目前使用粗略全球尺度土壤质地产品的基于模型的方法（如土壤安全、气候、水文、地表）。这里介绍的方法也适用于光学遥感数据有限和/或全年都有植被覆盖的地区。

**1. 引言**
土壤质地（即沙子、粉砂和粘土的比例）对于理解其他土壤属性非常重要，因为它影响土壤的持水能力。沙质土壤通常比粘土土壤含有更多的自由水（Srivastava等人，2006；Das & Paul，2015），而粘土土壤含有更多的结合水。然而，使用传统方法（如野外采样、实验室分析）获取土壤质地和其他土壤属性仍然成本高昂且耗时，尤其是在小规模（大面积）地理分析中。作为一种潜在且低成本的替代方案，合成孔径雷达（SAR）已被用于估算土壤属性（Petropoulos等人，2015；Babaeian等人，2019；Domenech等人，2020；Pradipta等人，2022）。这些研究通常集中在利用SAR估算裸地或植被稀疏地面的土壤湿度或碳含量，并将土壤质地作为次要参数来获取这些参数。很少有研究使用SAR估算永久被植被覆盖的土壤的颗粒大小分数（PSF，即沙子、粉砂和粘土的相对比例）（Jagdhuber，2012，2013；Harfenmeister等人，2021；Salma等人，2022；Dave等人，2023）。这可能部分是由于这种土壤属性具有组成性质（即相对比例之和等于100%），以及使用和解释雷达数据时遇到的挑战。

与微波遥感器中的能量强度测量相关的主要参数之一是归一化雷达截面（σ0），也称为后向散射系数，它是极化强度的度量（即垂直发射和接收-VV，垂直发射水平接收-VH）。信号返回或后向散射系数（σ0）在很大程度上取决于土壤质地（Dobson等人，1981）和土壤含水量，因为土壤的介电常数（ε）与土壤孔隙中的水分存在直接相关（Das & Paul，2015）。可以通过建立后向散射系数与感兴趣的土壤属性之间的关系来使用SAR估算土壤属性。

地统计方法主要通过距离和邻域分析来建模、估算、预测和模拟地理现象，基于地理学的三个基本原理或法则：空间自相关（Tobler，1970）、空间异质性（Goodchild，2004）和空间相似性（Zhu等人，2018）。通常有两种不同的方法来预测未采样位置的值（空间预测、空间插值或空间建模）——确定性方法（例如，逆距离加权-IDW、全局多项式插值-GPI、径向基函数-RBF、局部多项式插值-LPI）和非确定性方法（例如，克里金、共克里金、面积插值、经验贝叶斯克里金-EBK）（Singh & Sarma，2023）。空间机器学习是另一种方法，但它与基于地统计的模型不同，因为它需要一些严格的假设，如目标变量的分布和稳定性（Hengl等人，2018）。也可以用这种方法测量空间预测的不确定性，有关方法的详细信息可以在Hengl等人（2018）中找到。

然而，在处理像土壤质地这样的组成数据时，无论是空间建模还是非空间建模都会遇到挑战，因为这种变量的本质是组合性的（即组合响应或组合预测因子），预测/估算需要反映这一限制。Wang & Shi（2017）得出结论，为了提高空间插值的性能，需要结合地理空间方法和辅助变量的使用，因为许多传统方法可能无法满足组成数据插值的特定要求。关于遥感，挑战可能包括影响光学卫星的天气条件（如云层）和植被覆盖的土壤（例如，雷达信号穿透地形表面的能力）。

在这项研究中，测试了两种基于地统计方法的空间模型（经验贝叶斯克里金回归和共克里金）和两种非空间模型（线性回归和基于树的回归），使用微波遥感数据（Sentinel-1 C波段SAR）来估算表土层中沙子、粉砂和粘土的百分比。这里，我们根据模型中使用的输入和输出格式来区分空间模型和非空间建模方法。例如，在非空间方法中，输入和输出都是表格形式，而在空间方法中我们使用地理栅格和形状文件。由于数据格式的不同，这里使用了“明确”一词来区分这些统计建模方法。因此，非空间和空间方法的输出分别是表格和地图。所采用的建模方法明确考虑了土壤质地的组成性质。这项工作旨在解决两个关键研究问题：（i）如何利用基于微波的遥感数据和技术来估算土壤属性以进行土壤安全评估？（ii）空间模型在估算土壤颗粒大小分数方面是否优于非空间模型？本研究重点评估了在模型中使用数值或连续但非地理数据集（例如，栅格、形状文件）与空间数据集（例如，栅格、形状文件）进行土壤成分空间预测的影响。

还值得注意的是，在这项研究中，土壤质地类别是根据从模型方法获得的沙子、粉砂和粘土的定量估计（%）得出的，而不是通过分类方法预测的。换句话说，土壤质地类别是根据沙子、粉砂和粘土含量的组成预测推断出来的，而不是基于分类建模方法。这种方法类似于在实验室测量PSF时的方法，其中估计值被绘制在三元图上以确定最终的土壤质地类别。

本文的其余部分组织如下：第2节简要回顾了理论和相关文献。第3节概述了进行研究的方法和数据，包括研究区域。第4节和第5节分别展示了建模结果和讨论；第6节描述了工作的优点和局限性，并提出了研究结论。

**2. 背景**
**2.1. 土壤的微波遥感**
自20世纪70年代以来，主动和被动微波遥感技术已被用于土壤属性的估算，采用了空间和非空间方法（Babeian等人，2019；Domenech等人，2020）。SAR后向散射系数（σ0）常用于检索或估算土壤属性。它是与观测几何和地表相关的变量的归一化函数（Beale等人，2019）。然而，从植被覆盖的土壤中通过SAR数据估算土壤属性是一项具有挑战性的任务（Wang等人，2023），因为叶子（树木、作物）或草可能会覆盖整个土壤表面，导致传输的雷达信号和来自下层土壤的后向散射减弱（Beale等人，2019）。因此，在SAR土壤属性估算中，植被通常是必须考虑的最重要因素。一些SAR系统规格，如入射角、极化方式、波长/频率，如果在图像预处理过程中适当应用，可以减少土壤粗糙度和植被对雷达信号的影响（Ulaby等人，1978，1979；Mattia等人，2006；Moran等人，2012；Barrett等人，2009）。例如，较浅的入射角（>40度）可以最大化对作物条件的响应（Moran等人，2012），而入射角低于20度时植被的影响最小（Ulaby等人，1979）。对于裸露田地的土壤粗糙度，入射角在7到17度之间时表面粗糙度的影响最小（Ulaby等人，1979）。

解释土壤在微波光谱中行为的一个理论是土壤中自由水和结合水之间的关系（Dobson等人，1981），这受到土壤颗粒大小分数分布的控制。当土壤含水量增加时，水可以在土壤颗粒周围更自由地移动，重要的是要强调自由水对介电常数有主导影响（Barret等人，2007）。Hallikainen等人（1985）评估了在1.4 GHz至18 GHz频率范围内，土壤-水混合物的微波介电行为作为水分含量、温度和土壤质地组成（粘土、粉砂和沙子）的函数，发现土壤质地在整个频率范围内都有影响，尤其是在5 GHz以下的频率更为明显。因此，土壤质地影响后向散射系数，考虑到沙质土壤含有比粘土土壤更多的自由水（Srivastava等人，2006；Das & Paul，2015），而粘土土壤含有更多的结合水。Jackson和Schmugge（1989）证明了水分子吸附在土壤颗粒上，有效地固定了它们的偶极子，限制了结合水与雷达信号的相互作用（Das & Paul，2015）。使用SAR数据检索土壤属性的全面方法见补充表1——SAR数据中土壤属性检索技术的总结。

**2.2. 空间显式模型和地统计**
空间显式模型至少满足以下要求/测试之一及其组合（Goodchild，2001；Janowicz等人，2020）：不变性测试、表示测试、公式测试和结果测试。Nikparvar和Thill（2021）讨论了影响机器学习性能的空间显式数据的属性——空间依赖性、空间异质性、尺度、数据的特定表示（多边形、线、点、规则网格、文本）和测量过程。有关这些要求/测试和讨论的详细信息可以在Janowicz等人（2020）和Nikparvar和Thill（2021）中找到。

地统计方法基于空间自相关（附近的事物相关）、空间异质性（区域间的方差）和空间相似性（地理配置的相似性）这三个地理学法则（Tobler定律（Tobler，1970）、Goodchild第二定律（Goodchild，2004）和Zhu定律（Zhu等人，2018）。许多地理空间技术基于这些法则，包括Voronoi映射（第一定律）和克里金/共克里金（第二定律）。第二定律意味着Tobler定律的要求在大范围内可能无法满足，因为空间自相关和/或相关性可能会随空间变化。Delaunay三角剖分会产生Voronoi图，其中共享边的特征将被视为邻居——邻近图。这是一种有用的方法，用于获取特征的有影响力区域（Evans & Jones, 1987）（例如，代表土壤颗粒大小分数的点）。Voronoi图，也称为Dirichlet镶嵌或Thiessen多边形，表达了自然界中普遍存在的模式（Dierking et al., 2021）。因此，Voronoi图及其相关图或对偶图Delaunay三角剖分（补充图1）是分析空间数据的潜在有价值工具（Okabe et al., 2000）。空间相关性通常使用半变异函数来测量，这种相关性随距离增加而减小。原点附近的曲线越陡峭（即“块金效应”），表示最近邻居对预测/感兴趣点的影响越大。然而，由于空间自相关性，空间模型容易产生归纳性空间偏差（Smith, 2009; Mialon, 2023）。估计未采样位置数据的最常见地质统计方法是Kriging，它依赖于数学和统计模型以及空间异质性原理（Goodchild的第二定律）。这些特点使Kriging方法区别于更确定性的方法，如IDW。确定性方法由于不涉及随机性而具有可重复性和稳定性（Singh & Sarma, 2023），在需要重复性、可靠性和精确性的情况下非常有用（Singh & Sarma, 2023）。另一方面，它不能提供预测的不确定性估计。此外，非确定性方法试图量化预测中的不确定性。例如，Cokriging（CK）通过使用两个或多个预测因子来结合多个变量之间的空间相关性，从而可以改进模型并提高估计的准确性（Singh & Sarma, 2023）。当模型中考虑协变量时，会使用Cokriging，该方法利用响应变量的自相关性和预测因子与协变量之间的互相关信息来改进预测。Cokriging是Kriging的一种变体。除了创建空间地图外，非确定性方法还试图量化预测中的不确定性。此外，像Cokriging这样的方法可以通过使用两个或多个预测因子来结合多个变量之间的空间相关性，从而改进模型并提高估计的准确性（Singh & Sarma, 2023）。在预测建模之前，对响应变量进行了探索性的地质统计分析，该分析基于模式和邻域方法，以了解数据的空间结构并相应地调整模型中的参数。这种程序基于这样一个事实：空间插值方法依赖于周围测量值的权重来预测未测量位置的估计值。此外，在评估空间影响时，考虑采样单元在空间中的分布是否规则（例如，规则的或网格化的、随机的、分散的、聚集的）是值得的。关于模式分析，分别使用了全局Moran指数（GMI）和平均最近邻（ANN）（(2)，(3)。

(2) GMI = nS0 ∑i=1n ∑j=1n wi,j zi/zj ∑i=1n zi2，S0 = ∑i=1n ∑j=1n wi,j，其中 zi 是特征 i 的属性与其均值的偏差；wi,j 是特征 i 和 j 之间的空间权重；n 是特征的数量；S0 是所有空间权重的总和。

(3) ANN = D?o / D?e，D?o = ∑i=1n di，D?e = 0.5n/A，其中 D?o 是每个特征与其最近邻之间的观测平均距离；D?e 是在随机模式下特征的预期平均距离；di 表示特征 i 与其最近邻特征之间的距离；n 是特征的总数；A 对应于指定的面积值。

计算 GMI 是为了评估响应变量和协变量内的空间自相关性。其范围从 -1 到 1，其中正值表示相似的值聚集在一起（正空间自相关性）；负值表示数据之间的差异，即不同的值彼此相邻（负空间自相关性）；零值表示没有空间自相关性。值得注意的是，自相关性违反了经典统计学的核心原则，即观测值是相互独立的。ANN 有助于发现数据的分布模式。它测量每个特征与其最近邻位置之间的距离。如果 ANN 指数小于 1，则表明存在聚集趋势；如果大于 1，则表明趋势是分散的。

对于邻域分析，使用了距离带（Distance Band）和Delaunay三角剖分（Delaunay Triangulation）方法，因为它们更适合处理点特征。前者是指每个特征被包含为邻居的关键距离。在这里，根据数据的分布模式，找到的距离是 20,000 米。我们选择了曼哈顿距离（Manhattan distance）方法来考虑异常值和空间不连续性（例如，地形、地貌）。

3.3.2 预测建模方法和插值

空间模型
协同克里金法（Co-kriging，CK）：
基于数据的探索性分析（已知均值），采用了简单的 CK（公式 4）：
Z1(s) = μ + ?1(s)，Z2(s) = μ + ?2(s)，其中 Z1 和 Z2 是位于 (s) 的随机变量；μ 是已知且恒定的均值；?1 和 ?2 是随机误差项。协同克里金法假设两个模型和数据的平稳性（均值和方差恒定）。
协同克里金法中使用的变异图（基于方差的交叉变异图，公式 5）是从数据中估计的，用于解释土壤粒子尺寸函数（PSF）的空间变异性（Wang & Shi, 2018）：
γ^u,v(h) = 1/2N(h) ∑i=1N(h) {zu(xi) ? zu(xi+h)} / {zv(xi) ? zv(xi+h)}，其中 γ^ 是估计的变异性；h 是滞后或距离；N(h) 是按 h 分组的样本位置对；i th 是位于 xi+h 的观测值；(z(xj), j=1, 2, …, n) 是用于估计交叉变异图的数据；zu 和 zv 是空间位置 u 和 v 的实际方差（Wang & Shi, 2018）。
简单的 CK 也被采用，以与 EBKR 保持一致，后者依赖于简单的克里金法进行插值步骤。对于滞后大小，使用了 3,000 米的值来拟合交叉变异图，这相当于基于探索性地质统计分析和研究区域内原位土壤采样分布的任意两点之间最大距离的一半。根据探索性分析，并遵循找到任意两点之间最大距离的一半并除以滞后大小的规则，使用了 20 个滞后。由于在探索性分析中观察到沙子和粘土的趋势，本研究考虑了这些成分的各向异性（分别为 199° SSW-NNE 和 134° SE-NW）。在测试了不同的经验半变异图模型（稳定的、圆形的、球形的、指数的和高斯的、旋转二次的以及 k-贝塞尔的）后，选择了基于拟合和捕捉数据空间结构能力的稳定模型。该模型使用留一法交叉验证（LOOCV）进行了验证。

EBK 回归（EBKR）：
对于 EBKR 模型（公式 6）（Krivoruchko & Gribov, 2019），选择了最小方差百分比等于 95% 的指数半变异图，子集大小为 36（以考虑局部影响和局部方差），并进行了 100 次模拟：
z1 = y(si) + εi，i=1...K?，其中 z1 是观测位置 si 处的测量值；y(s) 是在位置 si 处研究的高斯过程；εi 是高斯测量误差；K 是测量次数。
EBKR 模型也使用留一法交叉验证（LOOCV）进行了验证。为了处理 CK 和 EBKR 模型结果中观察到的非线性，本研究在所有协变量上应用了多项式函数（三次的和六次的）和样条函数以及交互项。只有统计上显著的交互项被用作协变量来重新运行 EBKR 和 CK 模型（补充图 1）。

非空间模型：
为了比较目的，还进行了非空间线性建模。考虑到土壤质地的组成性质（颗粒大小成分的组成），本研究采用了基于对称性的对数比率变换方法，因为它们保留了距离。对数比率变换（ILR）被应用于响应变量 Y（沙子、粉砂和粘土），以拟合组成线性回归模型（Y-LRM），遵循 CoDaPack 软件中内置的默认分区（Comas-Cufí & Thió-Henestrosa, 2011）。对响应变量应用了中心对数比率变换以训练基于树的模型（Deodoro 等人，2023）。当预测值被反变换回原始单位时，保证 PSF 组件的总和等于 100（或 1）（Comas-Cufí & Thió-Henestrosa, 2011）。对响应变量应用了 CLR 变换以拟合基于树的回归模型（随机森林模型-RF），使用 Ranger 包（Wright & Ziegler, 2017）。对 RF 模型的参数进行了收敛性测试（补充表 2 和 3），并选择了默认值（即 ntree = 500 和 mtry = 3），因为不同的 ntree 和 mtry 值并没有显著提高模型的准确性。关于 CLR 和 ILR 变换的详细信息在补充信息 - 建模中提供。在组成 LRM 中也应用了交互项，以考虑预测变量和响应变量之间的相互作用，就像 RF 和地质统计模型一样，其中影响过程之间的相互作用已经嵌入在算法中。

地质统计分析和建模方法是使用免费和开源软件 GeoDA（Anselin, 2004）以及 ArcGIS Pro? 中可用的扩展 Geostatistical Analyst 和 Spatial Statistics 执行的。非空间 RF 模型是使用 R 语言和“Ranger”包（Wright & Ziegler, 2017）执行的。

3.4 土壤颗粒大小预测的评估指标
计算了均方根误差（RMSE）（公式 7）、平均绝对误差（MAE）（公式 8）和 Nash-Sutcliffe 效率（NSE）（Nash 和 Sutcliffe, 1970）（公式 9），以评估使用 LOOCV 验证的模型的性能。NSE 衡量残差方差与测量数据方差的幅度（Nash 和 Sutcliffe, 1970），并指示观测数据与预测数据（模拟数据）的图表与 1:1 线的吻合程度。值为 1 表示模型与观测数据完美匹配；NSE = 0 表示预测数据和观测数据的准确性相同；-Inf < NSE < 0 表示模型性能较差，观测均值是更好的预测器。
(7) RMSE = 1/n ∑i=1n (yi ? y^i)2，其中 yi 是因变量的实际值，y^i 是因变量的预测值，n 是观测值的数量
(8) MAE = 1/n ∑i=1n |yi ? y^i|，其中 yi、y^i 和 n 表示与 RMSE 相同的描述符。
(9) NSE = 1 ? [∑i=1n (yi ? y^i)2 / ∑i=1n (yi ? yi?)2]，其中 yi 是第 i 个观测值，y^i 是第 i 个预测值，yi? 是观测数据的均值，n 是观测值的总数。

4. 结果
沙子和粉砂的响应变量分布近似为正态分布，而粘土的分布则向右偏斜（补充图 2）。根据均值的分组，沙子和粘土成分（测量的百分比含量）在最高值方面显示出地理模式（图 2a-b）。前者在研究区域的 S-SE 部分显示出较大的值（> 45%），而粘土含量较高的值（> 22%）出现在西部区域（NW-SW）。粉砂含量较高的值（> 39%）倾向于 W-E 方向，尽管其模式不明显——最低和最高的值都分散在整个区域内。趋势分析图（图 2 和补充图 3）确认了土壤 PSF 的分布模式。这些结果表明，局部数据采样代表了一个非平稳过程，这在 Voronoi 图中描绘的局部方差和局部变化中可以观察到（补充图 4）。这些结果与 Singh 和 Sarma (2023) 的观点一致，他们强调土壤性质的空间相关性可能因地区而异，因此平稳性和 Tobler 定律可能不成立。这意味着由于面积较大，实现了地理学的第二定律（空间异质性），并与克里金方法一致。

下载：下载高分辨率图像（1MB）
下载：下载全尺寸图像
图 2. 根据 a) 距离带和 b) Delaunay 三角剖分（Voronoi 镶嵌）方法对响应变量（测量数据，含量百分比）进行的邻域分析。图 “c” 中展示了数据的三维视角，其中样本点的位置绘制在 x, y 平面上。x, z 平面和 y, z 平面是投影到这些平面上的值的散点图，表示趋势方向（蓝色和绿色）。每个样本点上方的值由 z 维中棒的高度表示。
平均最近邻比率大于 1（表 2），因此数据表现出分散趋势，这与全局 Moran 指数-GMI 给出的分布模式一致（表 3）。GMI 还表明响应变量缺乏空间自相关性（另见补充图 17）。关于协变量，VH、表层土壤湿度和坡向显示出聚集模式。只有土壤湿度和坡向协变量是自相关的或空间依赖的。

表 2. 与邻居的距离带。
距离带（至少有一个邻居）
最小（米）
最大（米）
平均（米）
最近邻比率
土壤样本
16
43
2,932
6,312
1.28
曼哈顿距离

表 3. 空间自相关性：全局 Moran 指数。
土壤 PSF*
全局 Moran 指数
分布模式
P 值
显著性水平
空间自相关（或依赖）
沙子 -0.05
0
分散
0.01
0.10
否
粉砂 -0.05
78
分散
0.00
0.01
否
粘土 -0.00
22
随机
0.84
否
Sand.clr -0.04
99
分散
0.01
0.05
否
Silt.clr -0.05
51
分散
0.00
0.01
否
Clay.clr -0.00
22
随机
0.85
否
ilr.1 -0.00
22
随机
0.85
否
ilr.2 -0.06
69
分散
0.00
0.01
否
协变量*
VV -0.00
49
随机
0.98
否
VH
0.07
98
聚集
0.17
0.05
否
RVI
0.02
56
随机
0.39
否
Curvature -0.07
01
分散
0.07
0.10
否
表层土壤湿度（SSM）
0.24
69
聚集
0.00
0.01
是
Alpha -0.03
46
随机
0.42
否
Altitude -0.01
93
随机
0.70
否
Slope
0.03
27
随机
0.28
否
Aspect
0.96
00
聚集
0.00
0.01
是
来源：ArcGIS 输出。距离方法：曼哈顿距离；空间关系方法：逆距离。
* 与土壤样本（观测数据）相关。
增量距离，用于测量每个增加的距离（小规模或大面积）的空间聚类强度，显示只有粘土（测量数据）在两个距离上表现出聚类，大约为 7,500 米和 17,500 米（补充图 5）。这意味着 7.5 公里可以作为粘土的更局部分析尺度，这是第一个统计上显著的峰值（显著性水平 0.10）。对于沙子和粉砂，结果表明在研究区域内可能存在多个不同的空间过程（没有明显的空间模式）。

4.1 空间建模
结果显示 EBKR 和 CK 的行为不同。在估计值方面，根据模型评估指标（表 4）和观测数据与预测数据之间的一致性（图 3），CK 得出了更好的结果。然而，统计诊断显示预测沙子的误差中存在与曲率预测因子相关的模式（补充图 10-12）。这在 EBKR 的误差图中没有观察到（补充图 7 和 8）。CoDa 方法的使用并没有改善模型的预测。

表 4.评估指标是通过空间建模方法得出的，并使用LOOCV进行了验证。建模方法（验证方法：LOOCV）：土壤颗粒大小分布（PSF）模型评估指标包括RMSE和MAE。模型：EBK回归（非CoDa）：沙子14.11、1.10、0.2；粉砂10.08、0.00、0.2；粘土8.06、0.00、0.3；模型：应用二阶交互项的EBK回归（非CoDa）：沙子14.71、1.40、0.1；粉砂10.78、0.40、0.1；粘土8.36、0.10、0.3；模型：EBK回归Y-组成（clr转换）：沙子14.11、1.20、0.2；粉砂10.28、0.10、0.2；粘土8.56、0.20、0.2；模型：Cokriging（非CoDa）：沙子2.01、1.50、0.9；粉砂6.95、0.10、0.6；粘土6.64、0.80、0.5；模型：应用二阶交互项的Cokriging（非CoDa）：沙子3.62、0.70、0.9；粉砂7.86、0.00、0.5；粘土6.64、0.90、0.5；模型：Cokriging Y-组成（clr转换）：沙子3.32、0.40、0.9；粉砂7.05、0.30、0.6；粘土6.95、0.10、0.5。下载：高分辨率图像（1MB）；下载：全尺寸图像。图3：使用LOOCV方法验证的CK方法得到的观测值和预测值的散点图，包括是否应用CLR转换的情况。未对协变量应用交互项的模型。相反，在表面预测（插值表面）方面，EBKR地图似乎遵循了观测（测量）数据中的模式和趋势（图2），尽管观测值和预测值的散点图有所不同（图4、图5、图6）。CoDa方法有助于捕捉空间模型中测量数据的趋势。下载：高分辨率图像（1MB）；下载：全尺寸图像。图4：使用LOOCV方法验证的EBKR方法得到的观测值和预测值的散点图，包括是否应用CLR转换的情况。未对协变量应用交互项的模型。下载：高分辨率图像（1MB）；下载：全尺寸图像。图5：使用LOOCV方法验证的EBKR方法得到的表面预测值，包括是否应用CLR转换的情况（1-3和4-5），以及相应的LOOCV验证图。下载：高分辨率图像（879KB）；下载：全尺寸图像。图6：使用LOOCV方法验证的CK方法得到的表面预测值，包括是否应用CLR转换的情况（1-3和4-5），以及相应的LOOCV验证图。在半变异图中未观察到块金效应（nugget effect），这可能意味着测量误差或微观尺度变化（补充图6和9）。后者可能是由于探索性分析图和半变异图显示短距离内各对数据之间没有自相关性。4.2. 非空间模型对于非空间模型方法（数值模型），未对响应变量应用二阶交互项的组成LRM根据基于误差和基于一致性的指标都取得了更好的结果（表5）。另一方面，这些结果与组成RF的结果相似，特别是对于粘土和沙子组分。总体而言，使用CoDa方法可以改善非空间模型中的预测。表5：使用LOOCV验证的非空间建模方法得到的评估指标。建模方法（验证方法：LOOCV）：SPSF*模型评估指标包括RMSE和MAE。模型：Y-组成LRM（ILR转换；未应用二阶交互项）：沙子4.12、3.37、0.94；粉砂7.08、5.44、0.62；粘土7.00、5.02、0.46；模型：Y-组成LRM（ILR转换；应用二阶交互项）：沙子4.88、3.35、0.91；粉砂9.08、6.45、0.37；粘土8.93、6.16、0.13；模型：Y-组成RF（CLR转换）：沙子4.28、2.91、0.93；粉砂5.55、5.52、0.59；粘土7.13、5.27、0.44；模型：非CoDa方法：简单LRM（单变量响应）：沙子7.10、5.35、0.81；粉砂7.92、6.10、0.52；粘土7.36、5.33、0.41。?SPSF：土壤颗粒大小组分。4.3. 土壤质地类别——所有模型总体而言，从统计模型间接获得的土壤质地类别在空间和非空间模型中是相似的（图7），除了EBKR，它在类别数量和样本的显示模式（非线性）方面有所不同，这是通过三元图评估的。CK方法在三元图（图7f和7g）中显示的土壤属性分布与非空间LRM（未应用二阶交互项）的结果相似（图7b和7d）。下载：高分辨率图像（1MB）；下载：全尺寸图像。图7：以美国农业部（USDA）三元图格式显示的观测土壤质地类别，以及非空间（b-e）和空间模型（f-i）得到的类别。a：观测（测量）数据；b：未应用对数比率转换和交互项的LRM模型；c：应用ILR转换和交互项的LRM；d：应用ILR转换但未应用交互项的LRM；e：应用CLR转换的RF模型；f：应用CLR转换的CK模型；g：未应用CLR转换的CK；h：应用CLR转换的EBKR模型；i：未应用CLR转换的EBKR。5. 讨论探索性分析显示，尽管GMI值为负，但响应变量和协变量的空间结构在数据上存在差异而非相似性。半变异图可以被视为差异函数，因为观测值对之间的差异方差随距离增加而增加。分散的空间模式通常反映了高值和低值之间的某种排斥过程，这与土壤PSF一致。例如，在本研究中，沙子的百分比范围为6-83%。采样大小和密度可能影响了产生零值和负值的GMI值，因为样本配置（即网格化、随机和聚集）会影响数据中的自相关程度。数据的非平稳性也可能与密度和采样配置有关。本研究中土壤PSF之间缺乏空间相关性也表明缺乏空间依赖性，即观测值彼此独立。这并不违反传统统计学的主要假设。然而，这与地统计学的一个核心原则——Tobler的第一定律（空间自相关性）相矛盾。尽管如此，值得注意的是EBKR依赖于经典统计学（即LMM）进行预测。此外，由于空间异质性和研究区域的大小（大面积），地理学的第二定律（Goodchild的第二定律）似乎成立。尽管有这些地统计学结果，但拟合的空间模型对于在未测量站点空间估计土壤PSF以及捕捉较大距离的空间结构（例如，差异性）仍然很有用。尽管土壤PSF没有显示出强烈的局部空间自相关性，空间模型仍然有助于表示更广泛的空间模式（即捕捉整体模式）。此外，CK和EBKR依赖于线性回归进行预测或插值。虽然区域尺度是研究区域的地理分析单位，但增量的空间自相关性表明可以基于聚类的空间过程分析并绘制较小尺度的粘土图（大尺度）。尽管如此，这里进行的地统计学探索性分析对于拟合交叉半变异图和执行CK和EBKR模型非常重要。此外，由于观察到趋势且未检测到自相关性，拟合空间模型的目的是为了模拟响应变量中的空间模式的影响。此外，由于Voronoi图显示非平稳性，因此采用了CK进行方法比较。从这一点来看，本研究的主要发现是CK与观测数据相关，而EBKR与用于拟合局部模型的邻居数据相关，这对表面预测（插值）很有用。再次强调，数据的采样配置（分散）可能会影响结果，这引发了一个问题：更密集的网格是否会导致不同的解释？关于估计值（PSF百分比），根据NSE、RMSE和MAE指标，CK对沙子、粉砂和粘土的预测效果优于EBK回归。因此，CK是一个全局预测的良好模型。尽管有这些结果，但值得注意的是CK的NSE值在PSF之间变化较大（例如0.5-0.9），而EBKR的值相对恒定（例如0.2-0.3）。从误差诊断来看，EBKR在预测值和标准化误差方面显示出更一致的图表，因为CoDa和非CoDa方法在图表中没有观察到模式。相反，在CK模型的沙子图中仅观察到与预测器曲率相关的模式（补充图10-12）。对于粘土的表面预测（插值），CK取得了良好的结果。这表明曲率具有非线性关系，且沙子的预测可能受到CK模型的影响。在尝试处理这种非线性时，对曲率应用了多项式函数（三次和六次），但未成功（补充图14）。还使用了样条函数，但它们也不足以解决这个问题。二阶交互项对于使用CK估计沙子和粉砂特别重要，而在EBKR模型中效果不佳（表4和补充图15及16）。它们对粘土也有好处，无论是CK还是EBKR。尽管交互项导致沙子估计略微过拟合，但改善了模型的诊断。结果表明，在模型中使用交互项时不应考虑各向异性。这一点在粉砂表面预测中很明显，其半变异图是各向同性的（无趋势）。这是空间模型中的一个重要权衡，因为插值的目标就是如此。曲率是沙子交互作用中最常见的协变量；对于粉砂，还有曲率和alpha；对于粘土，则是坡度和alpha。Alpha考虑了来自表面的微波信息和植被。RVI是另一个基于雷达的协变量，它考虑了植被。这些是重要的指标，因为它们与表面湿度有关，并存在于粘土和粉砂模型的交互作用中。此外，它们对湿度敏感。有趣的是，表面土壤湿度（SSM）和方位角协变量显示出空间依赖性（空间自相关性），EBKR和CK模型通过半变异图捕捉到了这一点。从结果来看，雷达变量、土壤湿度和地形（地貌条件）之间的空间关系很清楚，特别是对于粘土，其表面预测图反映了这些关系。关于没有交互作用的表面预测（即PSF插值）地图，EBKR模型在捕捉局部趋势和模式方面表现更好，特别是应用CLR转换时，因为插值表面与Voronoi图和测量数据的邻域分析相似。由于CLR是一种一对一的转换，它保持了样本组分之间的等距（保持距离）对数比率，其中目标组分是分子，几何平均值由所有土壤组成（三个组分）构成分母（补充信息——建模），因此为模型提供了正确的相对PSF信息。Cokriging并不适合映射PSF预测（即插值）。然而，这些地图应谨慎解释，因为它们可能反映了土壤数据采样排列的伪影（即分散模式）。此外，这些地图代表了预测的局部和全局尺度效应，考虑到图例中显示的类别。经验贝叶斯克里金法（EBK）对于处理非平稳性很有用，因为它在地理上进行克里金分析，以捕捉不同空间中的潜在过程，分别处理局部方差。尽管如此，EBKR依赖于经典统计学（即LMM）进行预测。此外，由于空间异质性和研究区域的大小（大面积），地理学的第二定律（Goodchild的第二定律）似乎成立。尽管有这些地统计学结果，拟合的空间模型对于在未测量站点空间估计土壤PSF以及捕捉较大距离的空间结构（例如，差异性）仍然很有用。尽管土壤PSF没有显示出强烈的局部空间自相关性，空间模型仍然有助于表示更广泛的空间模式（即捕捉整体模式）。此外，CK和EBKR依赖于线性回归进行预测或插值。虽然区域尺度是研究区域的地理分析单位，但增量的空间自相关性表明可以基于聚类的空间过程分析和绘制较小尺度的粘土图（大尺度）。尽管如此，地统计学探索性分析对于拟合交叉半变异图和执行CK和EBKR模型非常重要。此外，由于观察到趋势且未检测到自相关性，拟合空间模型的理由是为了模拟响应变量中的空间模式的影响。此外，虽然非平稳性是测量土壤PSF的特征，但CK被用来比较方法。从这一点来看，本研究的主要发现是CK与观测数据相关，而EBKR与用于拟合局部模型的邻居数据相关，这对表面预测（插值）很有用。再次强调，数据的采样配置（分散）可能会影响结果，这引发了另一个问题：更密集的网格是否会导致不同的解释？关于估计值（PSF百分比），根据NSE、RMSE和MAE指标，CK对沙子、粉砂和粘土的预测效果优于EBK回归。因此，CK是一个全局预测的良好模型。尽管有这些结果，但值得注意的是CK的NSE值在PSF之间变化较大（例如0.5-0.9），而EBKR的值相对恒定（例如0.2-0.3）。从误差诊断来看，EBKR在预测值和标准化误差方面显示出更一致的图表，因为CoDa和非CoDa方法在图表中没有观察到模式。相反，在CK模型的沙子图中仅观察到与预测器曲率相关的模式（补充图10-12）。对于粘土的表面预测（插值），CK取得了良好的结果。这表明曲率具有非线性关系，且沙子的预测可能受到CK模型的影响。在尝试处理这种非线性时，对曲率应用了多项式函数（三次和六次），但未成功（补充图13）。还使用了样条函数，但它们也不足以解决这个问题。二阶交互项对于使用CK估计沙子和粉砂特别重要，而在EBKR模型中效果不佳（表4和补充图15及16）。它们对粘土也有好处，无论是CK还是EBKR。尽管交互项导致沙子估计略微过拟合，但改善了模型的诊断。结果表明，在模型中使用交互项时不应考虑各向异性。这一点在粉砂表面预测中很明显，其半变异图是各向同性的（无趋势）。这是空间模型中的一个重要权衡，因为插值的目标就是如此。曲率是沙子交互作用中最常见的协变量；对于粉砂，还有曲率和alpha；对于粘土，则是坡度和alpha。Alpha考虑了来自表面的微波信息和植被。RVI是另一个基于雷达的协变量，它考虑了植被。这些是重要的指标，因为它们与表面湿度有关，并存在于粘土和粉砂模型的交互作用中。此外，它们对湿度敏感。有趣的是，表面土壤湿度（SSM）和方位角协变量显示出空间依赖性（空间自相关性），EBKR和CK模型通过半变异图捕捉到了这一点。从结果来看，雷达变量、土壤湿度和地形（地貌条件）之间的空间关系很清楚，特别是对于粘土，其表面预测图反映了这些关系。关于没有交互作用的表面预测（即PSF插值）地图，EBKR模型在捕捉局部趋势和模式方面表现更好，特别是应用CLR转换时，因为插值表面与Voronoi图和测量数据的邻域分析相似。由于CLR是一种一对一的转换，它保持了样本组分之间的等距（保持距离）对数比率，其中目标组分是分子，几何平均值由所有土壤组成（三个组分）构成分母（补充信息——建模），因此为模型提供了正确的相对PSF信息。Cokriging并不适合映射PSF预测（即插值）。然而，这些地图应谨慎解释，因为它们可能反映了土壤数据采样排列的伪影（即分散模式）。此外，这些地图代表了预测的局部和全局尺度效应，考虑到图例中显示的类别。经验贝叶斯克里金法（EBK）对于处理非平稳性很有用，因为它在地理上进行克里金分析，以捕捉不同空间中的潜在过程，分别处理局部方差和全局方差。因此，EBKR似乎能够捕捉小尺度（较大区域）的响应变量的模式和趋势，因为它通过为每个子集建模每个半变异图来进行局部预测。这一发现与Wadoux等人（2020）的研究结果一致，他们也发现贝叶斯校准能够捕捉瑞士河流流域规模的模型输入、初始状态、参数和结构不确定性。EBKR还处理了非平稳性，这是农业区域的一个内在特征，因为土壤属性的复杂性（Gribov和Krivoruchko，2020）。此外，处理偏差-方差（或准确性-精度）权衡的一种方法是使用混合模型，EBKR模型在其回归阶段就是这种情况。另一方面，根据模型评估指标，CK在统计上取得了更好的结果，因为它利用了响应变量的自相关性和预测变量与协变量之间的交叉相关性。基于方程5中的h和N(h)部分（基于方差的交叉变异图），CK模型可能捕捉了土壤样本（响应变量）的空间变异性和分布，就距离和样本位置对而言——即采样密度或采样配置。在预测（插值）数据和观测（原位）数据之间的比较中，空间模型对于预测粘土很有用，而非空间模型在[数值]估计方面则不然。看来空间模型捕捉到了粘土与所用预测变量之间的一些关系。这可能与空间建模旨在发现和解释地理区域内的空间关系、模式和变化有关。这一点在本研究中得到了证实，因为探索性分析的结果表明响应变量中存在模式和趋势。此外，粘土是唯一在相对较短距离（至少7.5公里）上显示出自相关的组分。模型也似乎捕捉到了这一点。关于土壤质地（分类），CK与美国农业部三元图中的观测类别更为一致，因为CK在PSF百分比方面的估计优于EBKR，根据模型指标。三元图是根据各组分（成分）的数值计算得出的。与土壤PSF不同，CoDa方法对土壤质地没有影响，因为通过数值模型（推断）获得的类别无论是否应用了对土壤组分进行对数比率转换，结果都是相似的。因此，与连续型的PSF不同，在处理分类型的土壤质地时，响应变量的组成性质并不是问题。CoDa方法对于数字土壤制图可能很重要，因为使用该技术时土壤PSF的总和等于100%。尽管如此，结果显示，将土壤质地视为组成数据对于空间模型中的估计并没有带来实质性的改进，除了在插值表面（空间预测）方面略有改善。相反，在非空间模型中，CoDa方法得出了更显著的结果。Amirian-Chakan等人（2019年）在他们的研究中也发现了这一点，他们研究了使用对数比率转换方法（ALR、CLR、ILR）对土壤质地对有效土壤含水量和灌溉总量空间估计的影响，发现不同转换方法之间的估计结果没有显著差异。

尽管在空间模型中使用CoDa的方法没有效果，但由于这种方法可以确保预测样本中各组分单位的总和为100%，因此对于土壤PSF来说仍然很重要。

**本研究的优势与局限性**
本研究展示了比较两种对比鲜明的建模方法（基于统计学的方法（非空间模型）和基于地理学的方法（空间模型）在估算土壤PSF方面的相关性，为土壤PSF和土壤质地的制图提供了通用基础，特别是在像爱尔兰这样的地区，由于云层覆盖频率高，光学遥感数据的应用受到限制。
这项研究的另一个重要意义在于，模型表明土壤粒径的估计更多是基于空间模式过程，而不是基于空间自相关过程，这与所使用测量数据的空间配置有关。

然而，这项研究也存在一些局限性。首先，进行地统计分析时，样本大小与研究区域大小之间的相互作用是一个需要注意的问题，因为这种分析依赖于网格配置。其次，模型中缺乏土壤的介电常数作为协变量，而这一参数有助于通过后向散射系数（σVV0和σVH0）理解沙子-粉砂-粘土之间的相互作用以及土壤中的孔隙水含量，从而改进建模结果。

**结论**
研究结果表明，在两种统计建模方法中考虑地理位置（空间条件）和土壤质地的组成性质是重要的，因为这有助于进行土壤质地的数字制图。由于空间信息的复杂性，空间模型中会存在一定的归纳偏差，这是由于数据中存在的某种空间结构（如空间异质性和空间自相关（共线性）所导致的。因此，了解这些偏差背后的原理和公式对于选择合适的方法至关重要。在估算土壤性质的统计建模方法方面，研究人员可以根据研究目的和问题以及地理位置是否重要来决定采用空间/地统计方法还是非空间方法。
本研究提供了一个方法论框架，有助于设计成本高昂且劳动密集型的实地土壤调查（模型中不确定性较大的区域可能需要详细的土壤调查）。它还提供了一种在大范围内估算土壤性质的方法，并生成新的数据产品，用于土壤模型（例如，土壤安全、健康、侵蚀研究、水文、地表模型（如Noah LSM模型）、气候模型等目前使用粗粒度全球尺度土壤质地产品的模型）。作为未来的工作，我们计划利用Copernicus卫星Biomass的数据，该卫星在L波段工作，允许雷达信号穿透土壤约50厘米深度。这将通过扩展SAR数据的应用范围来促进土壤分析，超越表层土壤，并进一步推进SAR技术在土壤分析中的应用。Biomass SAR数据预计将于2026年1月可用。此外，我们还将探索遥感技术在改进土壤安全研究方面的潜力。

**作者贡献**
Sandra Cristina Deodoro：概念化；数据整理；形式分析；方法论；可视化；初稿撰写；审阅和编辑。
Rafael de Andrade Moral：形式分析；方法论；审阅和编辑；监督。
Reamonn Fealy：调查；资源获取；审阅和编辑；监督。
Tim McCarthy：资源获取；审阅和编辑；监督。
Rowan Fealy：形式分析；审阅和编辑；资源获取；监督。

**数据可用性**
土壤数据集可在以下链接获取：
https://esdac.jrc.ec.europa.eu/content/lucas-2009-topsoil-data
https://esdac.jrc.ec.europa.eu/content/lucas2015-topsoil-data
https://www.isric.org/explore/wosis
Sentinel 1 C波段SAR数据可在以下链接获取：https://dataspace.copernicus.eu/
Copernicus DEM ESA EEA-10数据可在以下链接获取：https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model
来自Sentinel 1的地表土壤湿度（SSM）数据可在以下链接获取：https://land.copernicus.eu/global/products/ssm
本研究中提到的其他土壤数据集可根据需求提供。

**未引用的参考文献**
Creamer和O’Sullivan, 2018; Dobson和Ulaby, 1981; Gholizadeha等人, 2018; Met éireann（爱尔兰气象服务），2021; Wang和Shi, 2017; Wang和Shi, 2018; Barrett等人, 2007

**CRediT作者贡献声明**
Sandra Cristina Deodoro：初稿撰写；可视化；方法论；形式分析；数据整理；概念化。
Rafael de Andrade Moral：审阅和编辑；监督；方法论；形式分析。
Reamonn Fealy：审阅和编辑；监督；资源获取；调查。
Tim McCarthy：审阅和编辑；监督；资源获取。
Rowan Fealy：审阅和编辑；监督；资源获取；形式分析。