数字土壤制图（Digital soil mapping，DSM）的应用越来越广泛。例如，Hengl等^[1]建立了SoilGrids全球数字土壤图系统；Poggio和Gimona^[2]应用DSM估算土壤有机碳（SOC）储量；Dai等^[3]将DSM产品应用于区域陆面过程模拟。

DSM的模型模拟方法主要有：基于土壤发生学理论的环境相关模型；基于地理学第一定律的空间自相关模型；前两者的混合模型^[4]。环境相关模型包括多元线性回归（Multivariable linear regression，MLR）、广义可加模型（Generalized additive model，GAM）、判别分析、分类回归树（Classification and regression tree，CART）、随机森林（Random forest，RF）、人工神经网络（Artificial neural network，ANN）等；空间自相关模型如反距离加权（Inverse distance weighted，IDW）和地统计学等；混合模型主要是回归克里格（Regression kriging，RK），即MLR与普通克里格（Ordinary kriging，OK）的结合。

理论上，混合模型由于同时将环境变量与土壤属性的空间自相关纳入考虑，因此相对其他单一模型而言，具有明显的优势。例如，混合模型适用于复杂地区制图、可以降低制图成本等^[5]。很多研究表明混合模型的准确度较好。例如，李启权等^[6]研究发现径向基函数神经网络（Radial basis function neural network model，RBFNN）-OK混合模型和RK在预测精度上优于OK；Ma等^[7]研究发现决策树（Cubist）-简单克里格（Simple kriging，SK）混合模型在较大尺度上准确度高于单一的决策树模型。然而，也有很多研究表明，混合模型的准确度不如单一模型^[5，8]。这些研究进一步分析后认为，混合模型的准确度与环境变量和土壤属性之间的关系、土壤属性空间自相关程度和样本间距等有关^[5]。但是，这些因素对混合模型准确度的影响难以判断^[5，8]。

广义线性地统计模型（Generalized linear geostatistical model，GLGM）也是一种混合模型。它是广义线性混合模型（Generalized linear mixed model，GLMM）与地统计学的混合，由Diggle等^[9]提出。GLGM与RK类似，均对回归模型的残差项进行空间分析，但与RK不同的是，它包含了随机项，并引入连接函数从而适应响应变量的不同随机分布形式。因此，GLGM相对于广泛使用的RK而言具有一定的优势。最明显的优势是它通过随机项解决了土壤变异的非平稳性。尤其是当前，土壤变异的多尺度性逐渐在DSM研究中受到重视^[10]。例如Sun等^[11]探索了多尺度土壤制图方法。另一个优势是GLGM可以适应响应变量的不同分布形式，尽管这个优势在DSM中并不特别明显，因为土壤属性数据一般符合或经过转化后符合正态分布。综上，GLGM在DSM中的应用具有一定的优势，在一些研究中得到了应用，例如Kempen等^[12]利用GLGM更新荷兰泥炭地区的土壤图；Poggio和Gimona ^[2]利用与GLGM类似的GAM-地统计混合三维模型（Hybrid GAM-geostatistical 3D model）模拟苏格兰的土壤有机碳垂直与横向分布。

另一方面，GLGM也具有明显的劣势，主要是形式复杂多样导致的计算繁琐问题，以及参数和变异函数的最优选择等问题^[12]。一般研究人员要依赖较强的专业知识或者大量的计算，从众多的环境变量和多样的模型组合（即固定变量、随机变量、空间变异函数）中筛选出最优的模型及其参数。因此，GLGM较难使用，在国内外的DSM应用中较少。截至目前，还没有研究表明其准确度是否优于常用的DSM方法。

本文以广西南宁高峰林场林业示范基地内一小面积（3.03km²）丘陵为研究区，利用地形因子，建立SOC、pH、黏粒和阳离子交换量（CEC）的GLMM及其与OK相结合的模型（GLGM），并与MLR、GWR、回归森林（Regression Forest，RF）、OK、RK和GAM等常用的DSM方法进行准确度比较，探讨产生准确度差异的原因，评价GLMM和GLGM在DSM中的使用效果，从而分析GLMM和GLGM在DSM中的应用可行性，为提高DSM的准确度提供参考。

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

本文的研究区位于广西南宁高峰林场林业示范基地内，经纬度范围为22°57′8″—22°58′41″N，108°20′57″—108°21′54″E，南北距离约为2 820 m，东西距离约为1 710 m，面积约为3.03km²。本区高程约为126~301 m，属于典型的丘陵地形（图 1）。本区气候为亚热带湿润季风气候，年平均气温约为21.6℃，年平均降水量约为1 300 mm。土壤母质以古生代的泥岩、泥质页岩、砂页岩等沉积岩系的坡积物为主；主要土壤类型为赤红壤，相当于中国土壤系统分类中的强育湿润富铁土^[13]。本区自有记录以来一直用于林业，近20年以种植桉树人工林为主。

1.2 土壤采样与测定

本文的土壤采样设计共使用了4种方法，包括条件拉丁超立方抽样、网格采样、地形序列和随机采样。条件拉丁超立方抽样以六个地形因子：高程、坡度、坡向、平面曲率、剖面曲率、地形湿度指数为条件变量，共采集45个样点；网格采样法以250 m为间距，共采集45个样点；地形序列采样是在研究区的一个山谷中，沿与等高线相平、近似“Z”字形的道路上，选择等高线凸和凹的顶点处采样，共获得满足条件的20个地形序列样点；随机采样法即从研究区所有位置中，随机抽取45个样点。地形序列样点是在杨谦等^[14]的基础上增加的，其余样点均相同。这些地形序列样点间隔较近，有助于反映土壤属性在小范围内的变异，因而有助于本文研究GLGM。本文的土壤采样是在以样点为中心的1 m× 1 m正方形内，取中心点和四个顶角上0~20 cm的土样，充分混合后形成的土壤样本。这些土壤样本经过风干、研磨过筛后，用重铬酸钾氧化法测定SOC含量；用pH计（1︰2.5土水比）测定pH；用比重计方法测定黏粒含量；用交换法测定CEC。本文以条件拉丁超立方抽样、网格采样和地形序列三种采样方法所得的110个样点作为训练集，以随机采样所得的45个样点作为验证集，用于验证DSM准确度。

1.3 环境变量

本文的研究区面积较小，植被均一，气候和成土母质等环境要素也基本一致。因此，本文仅使用地形因子作为环境变量进行制图。本文在Arcgis10.2和SAGA软件中，使用10 m分辨率的数字高程模型，提取出14个地形因子，包括：高程（m）、坡度（%）、坡向（°）、偏北度（°）、偏东度（°）、偏北指数、偏东指数、太阳辐射（10³ KW·m^–2）、平面曲率、剖面曲率、比汇水面积（m²）、地形湿度指数、径流强度指数和地形位置指数。

1.4 制图模型

本文共使用8种制图模型：MLR、GWR、RF、OK、RK、GAM、GLMM和GLGM。其中，MLR、GWR、RF、OK、RK是DSM常用的模型，已有很多研究对它们进行了详细介绍，例如文献^[15]。因此，本文仅对其他三种模型在下文中进行详细介绍。在MLR模型中，本文首先依据方差膨胀因子（Variance inflation factor，VIF）剔除具有共线性的环境变量，再用逐步回归方法进行变量筛选和最优模型拟合，即当模型中所有环境变量达到显著（P < 0.05），且赤池信息准则（Akaike information criterion，AIC）最小时，为最优变量组合模型。在OK模型中，本文采用残差最大似然法（Residual maximum likelihood，REML）进行半方差函数模拟，从球状、高斯和Matérn函数中选择AIC值最小的作为最优半方差函数。其中Matérn函数的kappa系数以0.1为间隔，由0.1递增至3。在RK模型中，线性回归部分由MLR模型组成，残差的半方差函数选择方法与OK相同，整个模型使用REML方法模拟^[14]。

广义可加模型（Generalized additive model，GAM）由Hastie和Tibshirani提出^[16]。它是在广义线性模型（Generalized linear model，GLM）的基础上扩展而来，用每个环境变量的非线性光滑函数取代GLM中具体的参数来表达相应响应变量的函数，且光滑函数具有可加性。GAM形式为：

$ g(E(y)) = \alpha + \sum\limits_{j = 1}^p {{f_j}({x_j})} + \varepsilon $

(1)

式中，E(y)为响应变量$ y $的期望值，g()为连接函数，α为截距，f_j(x_j)表示第j个预测变量x_j的非参数光滑函数，ε为随机误差。光滑函数可以是样条函数、loess函数和核函数等。在实际研究中，光滑函数多使用样条函数进行拟合。本文使用样条函数中广泛使用的三次样条函数进行拟合，因为三次样条函数的光滑性和拟合效应均较好，也能良好适应数据的变化趋势。三次样条函数是由一系列节点分隔开的分段连续三次多项式组成：

$ f_j\left(x_j\right)=\left\{\begin{array}{l} S_{j 1}\left(x_j\right)=a_{j 1} x_j{ }^3+b_{j 1} x_j{ }^2+c_{j 1} x_j+d_{j 1}, t_1 \leq x_j \leq t_2 \\ S_{j 2}\left(x_j\right)=a_{j 2} x_j{ }^3+b_{j 2} x_j^2+c_{j 2} x_j+d_{j 2}, t_2 \leq x_j \leq t_3 \\ \cdots \cdots \\ S_{j(n-1)}\left(x_j\right)=a_{j(n-1)} x_j{ }^3+b_{j(n-1)} x_j^2+c_{j(n-1)} x_j+d_{j(n-1)}, t_{n-1} \leq x_j \leq t_n \end{array}\right. $

(2)

式中，$ {f_j}({x_j}) $为第$ j $个环境变量$ {x_j} $在一系列递增序列$ {t_1},{t_2},...,{t_n} $所划分的区间$ [{t_1},{t_2}],[{t_2},{t_3}],...,[{t_{n - 1}},{t_n}] $上的三次样条函数，$ {a_{j1}},{a_{j2}},...,{a_{j(n - 1)}} $、$ {b_{j1}},{b_{j2}},..., $ $ {b_{j(n - 1)}} $、$ {c_{j1}},{c_{j2}},...,{c_{j(n - 1)}} $、$ {d_{j1}},{d_{j2}},...,{d_{j(n - 1)}} $为这些多项式中的待定系数，$ {t_1},{t_2},...,{t_n} $称为节点，$ {S_{j1}}({x_j}), $ $ {S_{j2}}({x_j}),...,{S_{j(n - 1)}}({x_j}) $为这些区间中的三次多项式。在模型拟合时，三次样条函数要求这些位于不同区间的分段多项式在节点处保持连续，同时也保证两个相邻的多项式在节点处的一阶导数和二阶导数连续^[16]。本文将MLR模型中VIF筛选出的环境变量进行所有可能组合后，用于GAM建模。为了选择合适的n值，本文在模型拟合时，对每个环境变量光滑项的n值由小到大进行了试验，最后选取所有光滑项显著（P < 0.05）、光滑度指数（k-index）大于1并小于1.1，且AIC值最小的模型为最优模型。

GLMM在形式上是MLR与GAM结合而成，其中MLR作为固定项表达全局变异，而GAM作为随机项表达局部变异。GLGM则是GLMM与OK的混合。两者的形式分别为：

$ g(E(y({u_0}))) = \alpha + {X^T}\beta + \sum\limits_{j = 1}^p {{f_j}({x_j})} + \varepsilon $

(3)

$ g(E(y({u_0}))) = \alpha + {X^T}\beta + \sum\limits_{j = 1}^p {{f_j}({x_j})} + \sum\limits_i^{} {{\omega _i}} Z({u_i}) + \varepsilon $

(4)

式中，$ g() $为连接函数，$ y({u_0}) $为待估点$ {u_0} $上响应变量的值，$ E(y({u_0})) $为$ y({u_0}) $的期望值，$ \alpha $为截距，$ X $表示作为固定项的环境变量，$ \beta $为$ X $的参数，$ {f_j}({x_j}) $表示作为随机项的环境变量，$ Z({u_i}) $为采样点$ {u_i} $上的GLMM残差，$ {\omega _i} $为采样点$ {u_i} $上的权重（由半方差函数估计后计算得到），$ \varepsilon $为随机误差。与GAM模型模拟一样，本文采用VIF筛选后的环境变量进行所有可能的固定项与随机项的组合后，进行GLMM模拟，并选择AIC值最小的模型为最优模型。GLMM残差项的OK建模与上述OK建立过程相同。

1.5 模型准确度评价

在验证集的45个土壤样点上，计算土壤制图预测值的平均误差（Mean error，ME）、平均绝对误差（Mean absolute error，MAE）、均方根误差（Root mean square error，RMSE）、一致性相关系数（Lin’s concordance correlation coefficient，CCC）和拟合优度（Goodness-of-fit，R²），用于评价土壤制图的准确度。ME越接近0，MAE和RMSE越小，以及CCC和R²越接近1，表明模型预测的准确度越高，反之则越低。

2 结果与讨论 2.1 土壤属性描述性统计特征

表 1列出了土壤属性测定结果的统计特征。通过比较训练集、验证集的统计特征可以看到，各土壤属性在训练集和验证集上的统计特征相近，表明本文所采用的训练集与验证集对本研究区的各土壤属性均具有较好的代表性。训练集与验证集的变异系数值表明，除pH属于弱变异外，SOC、黏粒和CEC均属中等变异。八分位偏度（Octile skews）结果表明训练集各土壤属性均符合正态分布。

2.2 模型模拟结果 2.2.1 MLR、GWR和RF建模

SOC、pH、黏粒和CEC的MLR模型及其决定系数（R²）和调整决定系数（Adjusted R²）如表 2所示。调整决定系数结果表明，MLR模型分别能解释SOC、pH、黏粒和CEC空间变异的7.5%、19%、4.9%和18%，说明各土壤属性与地形因子之间的线性关系不强，以及MLR模型对各土壤属性空间变异的解释程度不高。在杨谦等^[14]的研究中，SOC、pH和黏粒的MLR模型的调整决定系数分别为2.0%、20%和9.8%。可见，由于本文多使用了20个地形序列样点进行模型模拟，SOC、pH和黏粒的MLR模型的调整决定系数分别增强、略降低和降低。Lai等^[17]的研究也表明，同一个研究区上不同数量的土壤采样点会影响模型的模拟结果及土壤制图的准确性。

使用各土壤属性MLR模型中的环境变量建立各自的GWR模型，经最小AIC选择，得到各自的GWR模型。模型的最优带宽分别处于266~880 m、538~1 143 m、415~1 012 m、787~1 799 m区间范围内。所有带宽均处于训练集样点最小间距（4 m）与最大间距（2 652 m）之间，模型的R²分别为0.293（P < 0.001）、0.343（P < 0.001）、0.126（P < 0.001）和0.222（P < 0.001），调整决定系数分别为0.280、0.311、0.118和0.208。

使用所有环境变量进行各土壤属性的RF建模，结果表明，SOC、pH、黏粒和CEC的RF模型的R²分别为0.564（P < 0.001）、0.540（P < 0.001）、0.536（P < 0.001）和0.542（P < 0.001），调整决定系数分别为0.499、0.471、0.466和0.475。

综上，对于所有土壤属性，MLR、GWR和RF模型拟合的调整决定系数的排序均为RF > GWR > MLR，说明非线性模型拟合程度优于线性模型，且局部回归的线性模型拟合程度优于全局回归的线性模型。这与郭澎涛等^[18]研究中RF的相关系数大于逐步线性回归（Stepwise linear regression，SLR）、赵明松等^[19]研究中GWR的调整决定系数大于MLR等的结果一致。

2.2.2 OK和RK建模

SOC、pH、黏粒和CEC的半方差分析结果表明，块金效应分别为41.9%、71.3%、100%和0%，表明除黏粒外，其他三种土壤属性分别具有中等、中等和很强的空间自相关性，而黏粒不具有空间自相关性。本文中SOC、pH和黏粒的半方差模型的块金效应分别大于、小于和大于杨谦等^[14]研究（即分别为19.3%、85.8%和91.1%）。如前所述，这主要是因为本文较杨谦等^[14]多使用了20个地形序列采样点。在Sun等^[20]的研究中，在同一区域内不同采样结果上建立的半方差模型也具有大小不同的块金效应。

RK的半方差分析结果表明，SOC、pH、黏粒和CEC经过MLR模拟之后，残差的块金效应分别为44%、80%、100%和0%，分别表明具有中等、很弱、无和很强的空间自相关性。其中SOC与pH的RK模型块金效应相对于OK均有提高，表明SOC与pH的空间变异适合MLR和地统计的混合模型；CEC的RK块金效应与OK相比无变化，表明CEC的空间变异不适合使用MLR和地统计的混合模型模拟。

2.2.3 GAM建模

GAM的模型模拟结果如表 3所示。其中CEC的GAM模型中存在3个自由度为1的光滑项，即s（偏东指数）、s（太阳辐射）和s（平面曲率）（表 3），表明这3个光滑项近似于线性项。图 2以SOC的GAM模型中部分光滑项为例，反映了土壤属性与环境变量的非线性光滑项趋势。在小于180 m处，SOC与高程呈负相关关系，可能是因为：在山脚海拔较低处，积水较多，植物生长茂盛，归还量高，因而SOC含量一般较高；随着海拔升高，积水减少，森林经营活动较为剧烈，植物生长略差，归还量降低，因而SOC含量降低；随着海拔继续升高至180 m~200 m处，SOC与高程呈正相关关系，可能是因为：森林经营活动逐渐减少，使土壤遭受的侵蚀降低，因而SOC含量有所提升；随着海拔继续升高至200 m~260 m处，SOC与高程呈负相关关系，可能是因为：随着海拔提升，汇水面积逐渐减小，土壤水分较低，植物生长较差，归还量降低，因而SOC含量降低；随着海拔继续升高至大于260 m处，SOC与高程呈正相关关系，可能是因为：海拔再升高时，林下草本植物逐渐增加，因而有利于SOC的累积^[21]。表 3中的调整决定系数表明GAM模型分别能解释SOC、pH、黏粒和CEC空间变异的61%、40%、27%、61%，均大于对应的MLR模型。这表明GAM模型中各土壤属性与地形因子的非线性拟合较好。

2.2.4 GLMM和GLGM建模

GLMM的模型模拟结果如表 4所示。调整决定系数表明GLMM模型能分别能解释SOC、pH、黏粒和CEC空间变异的40%、51%、30%、67%，均大于对应的MLR模型。并且，除SOC之外，GLMM的调整决定系数均大于GAM模型。SOC的GLMM模型与GAM模型具有相同的环境变量，但GAM模型存在2个自由度 > 10的环境变量，而GLMM模型中不存在自由度 > 10的环境变量，且有1个环境变量为固定项。因此，SOC的GLMM模型拟合程度较GAM低，在训练集上的调整决定系数较低。

综上，对于所有土壤属性的MLR、GWR、RF、GAM和GLMM模型，整体而言GLMM的调整决定系数最大，其次是RF、GAM和GWR，MLR最小，表示线性与非线性结合的模型拟合程度最高，其次是纯非线性模型，线性模型拟合程度最低。

如前文所述，GLGM为GLMM趋势项加上GLMM残差的OK插值。由于黏粒在GLMM模拟后残差的空间自相关性很弱，只能建立纯块金效应的半方差函数，因此不能建立黏粒的GLGM模型。其他三种土壤属性的GLGM模型的块金效应分别为84%、96%和4.0%，表明这三种土壤属性在GLMM模拟之后的残差分别具有很弱、很弱和很强的空间自相关性。而且，这三种土壤属性的GLGM模型的块金效应均较RK和OK高，表明GLMM更适合模拟这些土壤属性变异。

2.3 土壤制图与模型验证

上述模型预测SOC、pH、黏粒和CEC的空间分布图如图 3所示。MLR由于使用了土壤属性与地形因子之间的线性关系，因而预测图直接反映了研究区中所使用的地形因子分布。例如，SOC的MLR方程中使用了高程（表 2），因而预测图与DEM相似（图 3a）。由于pH的MLR模型中使用了较多的地形因子（表 2），因而其预测图对地形因子分布的反映较不明显（图 3b）。GWR使用的变量与MLR一致，且建模原理相似，因此各土壤属性的GWR预测图（图 3e~图 3h）与MLR模型的较为相似（图 3a~图 3d），但在高低值区域有较为集中的趋势，且空间邻近值之间的变化较为平滑。RF的预测图也显示出了地形分布，但由于使用了所有地形因子的非线性项，因此RF预测图与地形因子分布的差别较MLR大，并且细节更为丰富，值域范围相对于其他模型较窄（图 3i~图 3l）。OK的预测图中，值的变化较GWR更为平滑，且高、低值分别呈集中环状分布，反映了土壤属性分布的大致趋势（图 3m~图 3o）。RK由于结合了MLR与OK，因此既反映了地形因子分布细节，又反映了土壤属性空间相关性结构的大致分布格局（图 3p~图 3r）。

GAM由于使用了非线性光滑项，因此相较于其他模型细节更丰富，更详细地反映了土壤属性的局部变异（图 3s~图 3v）。由于pH与CEC的GAM模型使用了自由度接近于1的项，因此二者的GAM模型预测图均显示了地形信息。反之，SOC与黏粒的GAM预测图对地形信息的反映较不明显。GLMM与GAM建模方式类似，因此两者的预测图也类似。但由于SOC和pH的GLMM与GAM模型相近，而黏粒和CEC的GLMM与GAM模型差别较大，因此SOC和pH的GLMM与GAM的预测图更相近，而黏粒和CEC的GLMM与GAM预测图差别较大。此外，由于SOC、pH的GLGM中块金效应较高，而CEC的GLGM中块金效应较低，因而SOC、pH的GLMM与GLGM的预测图差别不大，而CEC的GLMM与GLGM预测图则差别较大。

基于45个随机验证样点的预测结果准确性检验结果如表 5所示。各指标均没有一致地表明哪种方法最优。以SOC为例，RMSE表明，最优模型为RF，最差模型为GAM，但CCC表明，最优模型为OK，最差模型为MLR。根据Khaledian和Miller^[22]分析，RMSE与ME和MAE相似，而CCC包含了R²。然而ME中正误差和负误差会相互抵消，使误差变小；MAE中所有残差的权重相等；R²不能捕捉模型预测中的偏差。因此本文综合使用RMSE和CCC作为主要的准确度指标进行分析。这两个指标表明，预测SOC准确度最好的模型是RF和OK，其次是RK，最差是MLR、GWR、GAM、GLMM、GLGM；预测pH准确度最好的模型是MLR和RK，其次是GWR，再次是OK和GAM，最差是RF、GLMM、GLGM；在黏粒上，排序为GLMM > RF > GWR > MLR > GAM；在CEC上，排序为OK > RK、GLMM、GLGM > RF > GAM > GWR > MLR。这些结果表明，GLMM与GLGM在黏粒上表现较好，在CEC上表现中等，在SOC、pH上表现较差，而GAM整体表现均较差。

为了分析预测制图准确度的影响因素，图 4展示了各土壤属性在模型模拟中的调整决定系数和块金效应，以及主要的准确度指标——RMSE和CCC。MLR模型的调整决定系数均在0.2以下，说明这些土壤属性的空间变异仅有少部分能由全局单一的线性模型解释。因此，MLR模型的准确度一般较差，仅在pH上因为块金效应较大（即较低的空间自相关性）且调整决定系数也较大时表现出较好的准确度。OK模型的块金效应表明，SOC和CEC具有较好的空间自相关性，而pH的空间自相关性较差。因此，OK在SOC和CEC上的准确度较好，而在pH上的准确度较差。同理，可以总结出其他模型在不同土壤属性上准确度较好或较差的原因。例如，RF在SOC上的决定系数较高，因而其在SOC上的准确度较高。

然而，上述规律也并不能完全解释所有结果。例如，GAM在四个土壤属性上均具有较高的调整决定系数，但其准确度总是较差。进一步分析认为，这些模型的准确度差异还与模型自身的特点有关。尤其是GWR、GAM、GLMM和GLGM等四个模型并不是全局模型，而是局部模型，即描述局部土壤空间变异的模型。这些模型一般具有较好的决定系数（GLGM具有与GLMM一样的决定系数），但它们并不一定具有较好的准确度。这是因为这些模型容易过度拟合，导致模型模拟结果实际上并不准确。例如，GAM将所有变量视为随机项，因而非常容易过度拟合而使模型结果失真。为此，本文计算了各土壤属性的全局莫兰指数^[23]，用以反映土壤属性的局部变异强弱。SOC的莫兰指数为0.31（极显著）；pH为0.09（不显著）；黏粒为–0.004（不显著）；CEC为0.18（极显著）。这些结果表明：SOC、CEC和pH的局部变异性较弱，而黏粒具有较强的局部变异性。因此，GAM虽然在SOC上具有最高的调整决定系数，但这可能是由于模型过度拟合，因而预测准确度实际较差。反之，GLMM、GLGM不仅用随机项来描述局部变异，同时也用固定项来描述全局变异，因而具有一定的优势，例如，GLMM在黏粒预测上具有最好的准确度。

综合以上分析，本文认为当土壤属性与环境变量具有较低的线性回归调整决定系数（约小于5%），且土壤属性具有较弱的空间自相关性（块金效应约大于71%）和较强的局部空间变异（全局莫兰指数约小于0.09）时，GLMM及其与地统计的混合模型——GLGM具有较高的适用性。对比本文中的多个模型，本文还认为，在选用模型进行DSM时，应根据土壤属性的实际变异情况来选用合适的模型。例如，Sun等^[5]在大量模拟数据的基础上，系统比较了OK和RK，也认为当模型能比较真实地反映土壤空间变异时，才表现出较好的预测准确度。

尽管本文初步分析得到了GLMM和GLGM在何种情况下准确度较好，但是使二者准确度高的线性回归调整决定系数、块金效应和莫兰指数等指标的具体定量范围仍有待进一步研究。而且本文的研究区面积较小，而在较大区域尺度中，由于气候、植被、母质等成土因素的尺度和作用强度不同，土壤属性的变异包含了多个尺度，以及不同的强度。因此，GLMM和GLGM在较大区域尺度中可能具有较好的适用性。值得一提的是，本文研究区的地形为丘陵，土壤属性变异较为明显，也有较为丰富的地形变量用于GLMM和GLGM。在地形平缓区，土壤属性的变异比较平缓，而可用于DSM的环境变量并不多^[24]。因此，GLMM和GLGM在大面积平缓区的应用还有待深入分析。此外，目前还缺乏高效地为GLMM和GLGM筛选合适的环境变量作为固定项和随机项的方法。因此，本文在建立GLMM和GLGM模型时，运用全变量组合方法来海量寻找最优的GLMM和GLGM模型，未来需进一步研究如何高效地为GLMM和GLGM筛选合适的环境变量作为固定项和随机项。

3 结论

本文以广西南宁高峰林场林业示范基地内的小面积典型丘陵为研究区，以14个地形因子为预测变量，比较了GLMM、GLGM与DSM常用模型（即MLR、GWR、RF、OK和RK）及GAM模型在SOC、pH、黏粒和CEC空间分布预测制图上的准确性。结果表明，当土壤属性与环境变量具有较低的线性回归调整决定系数（约小于5%），且土壤属性具有较弱的空间自相关性（块金效应约大于71%）和较强的局部空间变异（全局莫兰指数约小于0.09）时，GLMM和GLGM具有较高的适用性。同时，GLMM和GLGM预测的土壤图较其他常规模型具有更丰富的细节，更详细地反映了土壤属性的局部变异。鉴于土壤空间变异的高度异质性和多尺度性，本文认为GLMM和GLGM具有较好的应用前景，尤其是在包含了不同尺度、不同强度变异的较大区域上。这将为DSM研究中模型方法的选取提供有价值的参考。此外，GLMM和GLGM在模型模拟时的效率很低，未来还需研究如何高效地为GLMM和GLGM筛选合适的环境变量作为固定项和随机项，以提高效率。