广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用.pdf

第 60 卷 第 4 期 土 壤 学 报 Vol.60，No.4 2023 年 7 月 ACTA PEDOLOGICA SINICA Jul.，2023 *国家自然科学基金项目(42071062,41771246)资助 Supported by the National Natural Science Foundation of China(Nos.42071062,41771246)通讯作者 Corresponding author，E-mail：sun_ 作者简介：郝辰恺(1997)，男，广西桂林人，硕士研究生，主要从事数字土壤制图研究。E-mail： 收稿日期：20210729；收到修改稿日期：20211108；网络首发日期（）：20220321 http:/ DOI：10.11766/trxb202107290386 郝辰恺，孙孝林，王会利.广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用J.土壤学报，2023，60（4）：9931006.HAO Chenkai，SUN Xiaolin，WANG Huili.Application of Generalized Linear Geostatistical Model for Digital Soil Mapping in a Typical Subtropical Hilly AreaJ.Acta Pedologica Sinica，2023，60（4）：9931006.广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用*郝辰恺1，孙孝林1，王会利2（1.中山大学地理科学与规划学院，广州 510275；2.广西壮族自治区林业科学研究院，南宁 530002）摘 要：数字土壤制图在当前的应用越来越广泛，其方法主要包括环境相关模型、空间自相关模型，以及这些模型的混合模型。理论上，混合模型相对单一模型具有明显的优势。广义线性地统计模型（GLGM）也是一种混合模型，相对于最常用的混合模型回归克里格（RK），又具有能加入随机效应来解决土壤变异的非平稳性等优势。然而，GLGM 因计算繁琐等缺点，在国内外应用较少。本文以广西南宁高峰林场内一小面积（3.03 km2）丘陵为研究区，以 14 个地形因子为预测变量，使用广义线性混合模型（GLMM）及其与普通克里格（OK）相结合的模型（即 GLGM），对土壤有机碳（SOC）、pH、黏粒和阳离子交换量（CEC）的空间分布进行预测，并与常用的多元线性回归（MLR）、地理加权回归（GWR）、回归森林（RF）、OK、RK 和广义可加模型（GAM）进行比较。结果表明：GLMM 在预测黏粒上准确度较高；GLMM 和 GLGM 在预测 CEC上准确度中等，但在预测 SOC 和 pH 上准确度较低。综合线性回归模型的调整决定系数、块金效应和全局莫兰指数，本文认为，当土壤属性与环境变量具有较低的线性回归调整决定系数（即小于 5%）、土壤属性具有较弱的空间自相关性（即块金效应大于 71%）和较强的局部空间变异（即全局莫兰指数小于 0.09）时，GLMM 和 GLGM 具有较高的适用性，例如本文中的黏粒。反之，GLMM 和 GLGM 的适用性不好，例如 SOC 和 pH。鉴于土壤空间变异的高度异质性和多尺度性，GLMM和 GLGM 具有较好的应用前景。但是，今后研究还需进一步探讨如何提高 GLMM 和 GLGM 的模拟效率。关键词：数字土壤制图；土壤计量学；广义线性地统计模型 中图分类号：K903 文献标志码：A Application of Generalized Linear Geostatistical Model for Digital Soil Mapping in a Typical Subtropical Hilly Area HAO Chenkai1,SUN Xiaolin1,WANG Huili2(1.School of Geography and Planning,Sun Yat-sen University,Guangzhou 510275,China;2.Guangxi Forestry Research Institute,Nanning 530002,China)Abstract:【Objective】Digital Soil Mapping is receiving more attention and becoming widely used.Its methods mainly include 994 土 壤 学 报 60 卷 http:/ environmental correlation-based models,spatial auto-correlation based models,and a mixture of these two kinds of models.The mixed model is expected to be advantageous over the single models.A generalized linear geostatistical model(GLGM)is a kind of mixed model.Compared with the commonly used mixed model,i.e.,regression kriging(RK),GLGM has advantages such as having random effects to account for the non-stationarity of soil variability.However,GLGM is seldomly used due to its major disadvantages,i.e.,complicated computation.【Method】In this study,within a small hilly area(3.03 km2)in Gaofeng Forest of Nanning,Guangxi,generalized linear mixed model(GLMM)and its combination with ordinary kriging(OK),i.e.,GLGM,were used to predict the spatial distribution of soil organic carbon(SOC),pH,clay and cation exchange capacity(CEC).Performances of the two models were then compared with commonly used models,including multivariable linear regression(MLR),geographically weighted regression(GWR),regression forest(RF),OK,RK and generalized additive model(GAM).【Result】The results showed that GLMM had higher accuracy in predicting clay,while GLMM and GLGM had medium accuracy in predicting CEC.Further,GLMM and GLGM had lower accuracy in predicting SOC and pH.【Conclusion】Based on the adjusted R2 of the linear regression model,nugget effect and global Morans I,it is concluded that GLMM and GLGM are appropriate when there is a low adjusted R2 of linear soil-landscape regression(less than 5%),weak spatial auto-correlation of soil(nugget-to-sill ratio large than 71%),and strong local spatial variability of soil(Morans I less than 0.09),e.g.,clay in this paper.Otherwise,GLMM and GLGM are not appropriate,e.g.,for SOC and pH in this paper.For the high spatial heterogeneity and multi-scale variability of soil,we think that GLMM and GLGM are promising for DSM,although more studies are needed to improve the efficiency of GLMM and GLGM modelling.Key words:Digital soil mapping;Pedometrics;Generalized linear geostatistical model 数字土壤制图（Digital soil mapping，DSM）的应用越来越广泛。例如，Hengl 等1建立了 SoilGrids全球数字土壤图系统；Poggio 和 Gimona2应用 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 中并不特别明显，因为土壤属性4 期 郝辰恺等：广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用 995 http:/ 数据一般符合或经过转化后符合正态分布。综上，GLGM 在 DSM 中的应用具有一定的优势，在一些研究中得到了应用，例如 Kempen 等12利用 GLGM更新荷兰泥炭地区的土壤图；Poggio 和 Gimona 2利用与 GLGM 类似的 GAM-地统计混合三维模型（Hybrid GAM-geostatistical 3D model）模拟苏格兰的土壤有机碳垂直与横向分布。另一方面，GLGM 也具有明显的劣势，主要是形式复杂多样导致的计算繁琐问题，以及参数和变异函数的最优选择等问题12。一般研究人员要依赖较强的专业知识或者大量的计算，从众多的环境变量和多样的模型组合（即固定变量、随机变量、空间变异函数）中筛选出最优的模型及其参数。因此，GLGM 较难使用，在国内外的 DSM 应用中较少。截至目前，还没有研究表明其准确度是否优于常用的 DSM 方法。本文以广西南宁高峰林场林业示范基地内一小面积（3.03 km2）丘陵为研究区，利用地形因子，建立 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 研究区概况 本文的研究区位于广西南宁高峰林场林业示范基地内，经纬度范围为 22578225841N，10820571082154E，南北距离约为 2 820 m，东西距离约为 1 710 m，面积约为 3.03km2。本区高程约为 126301 m，属于典型的丘陵地形（图 1）。本区气候为亚热带湿润季风气候，年平均气温约为21.6，年平均降水量约为 1 300 mm。土壤母质以古生代的泥岩、泥质页岩、砂页岩等沉积岩系的坡积物为主；主要土壤类型为赤红壤，相当于中国土壤系统分类中的强育湿润富铁土13。本区自有记录以来一直用于林业，近 20 年以种植桉树人工林为主。图 1 研究区的数字高程图与土壤样点分布 Fig.1 Digital elevation map of the study area and distributions of sampling sites 1.2 土壤采样与测定 本文的土壤采样设计共使用了 4 种方法，包括条件拉丁超立方抽样、网格采样、地形序列和随机采样。条件拉丁超立方抽样以六个地形因子：高程、坡度、坡向、平面曲率、剖面曲率、地形湿度指数为条件变量，共采集 45 个样点；网格采样法以 250 m为间距，共采集 45 个样点；地形序列采样是在研究区的一个山谷中，沿与等高线相平、近似“Z”字形的道路上，选择等高线凸和凹的顶点处采样，共获得满足条件的 20 个地形序列样点；随机采样法即从研究区所有位置中，随机抽取 45 个样点。地形序列样点是在杨谦等14的基础上增加的，其余样点均相同。这些地形序列样点间隔较近，有助于反映土壤属性在小范围内的变异，因而有助于本文研究GLGM。本文的土壤采样是在以样点为中心的 1 m 1 m 正方形内，取中心点和四个顶角上 020 cm 的土样，充分混合后形成的土壤样本。这些土壤样本经过风干、研磨过筛后，用重铬酸钾氧化法测定 SOC含量；用 pH 计（12.5 土水比）测定 pH；用比重计方法测定黏粒含量；用交换法测定 CEC。本文以996 土 壤 学 报 60 卷 http:/ 条件拉丁超立方抽样、网格采样和地形序列三种采样方法所得的 110 个样点作为训练集，以随机采样所得的 45 个样点作为验证集，用于验证 DSM 准确度。1.3 环境变量 本文的研究区面积较小，植被均一，气候和成土母质等环境要素也基本一致。因此，本文仅使用地形因子作为环境变量进行制图。本文在 ArcGIS10.2和 SAGA 软件中，使用 10 m 分辨率的数字高程模型，提取出 14 个地形因子，包括：高程（m）、坡度（%）、坡向（）、偏北度（）、偏东度（）、偏北指数、偏东指数、太阳辐射（103 KWm2）、平面曲率、剖面曲率、比汇水面积（m2）、地形湿度指数、径流强度指数和地形位置指数。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）进行半方差函数模拟，从球状、高斯和 Matrn 函数中选择 AIC 值最小的作为最优半方差函数。其中 Matrn 函数的 kappa系数以 0.1 为间隔，由 0.1 递增至 3。在 RK 模型中，线性回归部分由 MLR 模型组成，残差的半方差函数选择方法与 OK 相同，整个模型使用 REML 方法模拟14。广义可加模型（Generalized additive model，GAM）由 Hastie 和 Tibshirani 提出16。它是在广义线性模型（Generalized linear model，GLM）的基础上扩展而来，用每个环境变量的非线性光滑函数取代 GLM 中具体的参数来表达相应响应变量的函数，且光滑函数具有可加性。GAM 形式为：1()()pjjjg E yfx （1）式中，E(y)为响应变量y的期望值，g()为连接函数，为截距，fj(xj)表示第 j 个预测变量 xj的非线性光滑函数，为随机误差。光滑函数可以是样条函数、loess函数和核函数等。在实际研究中，光滑函数多使用样条函数进行拟合。本文使用样条函数中广泛使用的三次样条函数进行拟合，因为三次样条函数的光滑性和拟合效应均较好，也能良好适应数据的变化趋势。三次样条函数是由一系列节点分隔开的分段连续三次多项式组成：32111111232222222332(1)(1)(1)(1)(1)1(),(),()(),jjjjjjjjjjjjjjjjjjjjjjj njj njj njj njj nnjnSxa xb xc xdtxtSxaxbxcxdtxtfxSxaxbxcxdtxt （2）式中，()jjfx为第j个环境变量jx在一系列递增序列12,nt tt所划分的区间12231,nnt ttttt上的三次样条函数，12(1),jjj naaa、12,jjbb(1)j nb、12(1),jjj nccc、12(1),jjj nddd为这些多项式中的待定系数，12,nt tt称为节点，1(),jjSx 2(1)(),()jjj njSxSx为这些区间中的三次多项式。在模型拟合时，三次样条函数要求这些位于不同区 间的分段多项式在节点处保持连续，同时也保证两个相邻的多项式在节点处的一阶导数和二阶导数连 续16。本文将 MLR 模型中 VIF 筛选出的环境变量进行所有可能组合后，用于 GAM 建模。为了选择合适的 n 值，本文在模型拟合时，对每个环境变量光滑项的 n 值由小到大进行了试验，最后选取所有光滑项显著（P0.05）、光滑度指数（k-index）大于1 并小于 1.1，且 AIC 值最小的模型为最优模型。GLMM 在形式上是 MLR 与 GAM 结合而成，其中 MLR 作为固定项表达全局变异，而 GAM 作为随机项表达局部变异。GLGM 则是 GLMM 与 OK的混合。两者的形式分别为：4 期 郝辰恺等：广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用 997 http:/ 01()()pTjjjg E y uXfx （3）01()()()pTjjiijig E y uXfxZ u（4）式中，()g为连接函数，0()y u为待估点0u上响应变量的值，0()E y u为0()y u的期望值，为截距，X表示作为固定项的环境变量，为X的参数，()jjfx表示作为随机项的环境变量，()iZ u为采样点iu上的 GLMM 残差，i为采样点iu上的权重（由半方差函数估计后计算得到），为随机误差。与GAM 模型模拟一样，本文采用 VIF 筛选后的环境变量进行所有可能的固定项与随机项的组合后，进行 GLMM 模拟，并选择 AIC 值最小的模型为最优模型。GLMM 残差项的 OK 建模与上述 OK 建立过程相同。1.5 模型准确度评价 在验证集的 45 个土壤样点上，计算土壤制图预测值的平均误差（Mean error，ME）、平均绝对误差（Mean absolute error，MAE）、均方根误差（Root mean square error，RMSE）、一致性相关系数（Lins concordance correlation coefficient，CCC）和拟合优度（Goodness-of-fit，R2），用于评价土壤制图的准 确度。ME 越接近 0，MAE 和 RMSE 越小，以及 CCC和 R2越接近 1，表明模型预测的准确度越高，反之则越低。2 结果与讨论 2.1 土壤属性描述性统计特征 表 1 列出了土壤属性测定结果的统计特征。通过比较训练集、验证集的统计特征可以看到，各土壤属性在训练集和验证集上的统计特征相近，表明本文所采用的训练集与验证集对本研究区的各土壤属性均具有较好的代表性。训练集与验证集的变异系数值表明，除 pH 属于弱变异外，SOC、黏粒和CEC 均属中等变异。八分位偏度（Octile skews）结果表明训练集各土壤属性均符合正态分布。2.2 模型模拟结果 2.2.1 MLR、GWR 和 RF 建模 SOC、pH、黏粒和 CEC 的 MLR 模型及其决定系数（R2）和调整决定系数（Adjusted R2）如表 2 所示。调整决定系数结果表明，MLR 模型分别能解释 SOC、pH、黏粒和 CEC 空间变异的 7.5%、19%、4.9%和 18%，说明各土壤属性与地形因子之间的线性关系不强，以及 MLR 模型对各土壤属性空间变异的解释程度不高。在杨谦等14的研究中，SOC、pH 和黏粒的 表 1 研究区土壤属性的描述性统计特征 Table 1 Statistics of soil properties in the study area 最小值 Min 最大值 Max 平均值Mean 中值 Median标准差 Standarddeviation偏度 Skewness峰度 Kurtosis 变异系数/%Coefficient of variation 八分位偏度Octile skew训练集（样本数 n=110）Calibration dataset（n=110）SOC/（gkg1）11.58 40.11 21.62 21.20 5.34 0.67 3.58 24.72 0.07 pH 3.95 5.10 4.41 4.41 0.20 0.49 3.80 4.55 0.06 黏粒 Clay/%28.00 52.00 38.74 39.00 4.56 0.08 2.70 11.77 0.01 CEC/（cmolkg1）11.40 26.80 16.15 16.00 2.44 1.00 5.80 15.13 0.01 验证集（样本数 n=45）Validation dataset（n=45）SOC/（gkg1）10.72 31.64 20.35 19.63 4.18 0.55 3.56 20.54 0.38 pH 4.01 5.04 4.42 4.41 0.22 0.63 3.69 4.97 0.09 黏粒 Clay/%25.80 52.00 37.26 36.80 6.39 0.38 2.53 17.16 0.05 CEC/（cmolkg1）10.20 20.40 15.52 15.40 2.54 0.13 2.26 16.34 0.02 998 土 壤 学 报 60 卷 http:/ 表 2 各土壤属性的 MLR 模型 Table 2 MLR models for soil properties 土壤属性 Soil property MLR 方程 MLR formulae 决定系数R2 调整决定系数 Adjusted R2 SOC/（gkg1）SOC=31.41 0.04 高程+83.88剖面曲率 0.092*0.075 pH pH=4.79 4.55 103 坡度 3.35 剖面曲率 4.30 平面曲率+5.81 106 径流强度指数 0.05 地形湿度指数 0.227*0.189 黏粒 Clay/%Clay=38.55 1.79 偏北指数 0.058*0.049 CEC/（cmolkg1）CEC=16.44+0.89 偏东指数+49.29 平面曲率 0.193*0.178 注：*表示显著性水平为 0.05，*表示显著性水平为 0.01，*表示显著性水平为 0.001。下同。Note：*represents significant at the 0.05 level，*represents significant at the 0.01 level，*represents significant at the 0.001 level.The same as below.MLR 模型的调整决定系数分别为 2.0%、20%和9.8%。可见，由于本文多使用了 20 个地形序列样点进行模型模拟，SOC、pH 和黏粒的 MLR 模型的调整决定系数分别增强、略降低和降低。Lai 等17的研究也表明，同一个研究区上不同数量的土壤采样点会影响模型的模拟结果及土壤制图的准确性。使用各土壤属性 MLR 模型中的环境变量建立各自的 GWR 模型，经最小 AIC 选择，得到各自的GWR 模型。模型的最优带宽分别处于 266880 m、5381 143 m、4151 012 m、7871 799 m 区间范围内。所有带宽均处于训练集样点最小间距（4 m）与最大间距（2 652 m）之间，模型的 R2分别为 0.293（P0.001）、0.343（P0.001）、0.126（P0.001）和0.222（P0.001），调整决定系数分别为 0.280、0.311、0.118 和 0.208。使用所有环境变量进行各土壤属性的 RF 建模，结果表明，SOC、pH、黏粒和 CEC 的 RF 模型的 R2分别为 0.564（P0.001）、0.540（P0.001）、0.536（P0.001）和 0.542（PGWR 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 m200 m 处，SOC 与高程呈 4 期 郝辰恺等：广义线性地统计模型在典型亚热带丘陵区数字土壤制图中的应用 999 http:/ 表 3 各土壤属性的 GAM 模型 Table 3 GAM models for soil properties 土壤属性 Soil property GAM 方程 GAM formulae 决定系数 R2 调整决定系数 Adjusted R2 SOC/（gkg1）21.62+s（高程，6.07）+s（坡度，11.34）+s（偏东度，12.73）+s（剖面曲率，6.32）0.625*0.610 pH 4.41+s（高程，3.04）+s（偏北度，4.97）+s（剖面曲率，1.64）+s（地形湿度指数，4.73）0.422*0.400 黏粒 Clay/%38.74+s（偏北指数，13.01）0.273*0.266 CEC/（cmolkg1）16.15+s（坡度，12.94）+s（偏东指数，1.00）+s（太阳辐射，1.00）+s（剖面曲率，6.24）+s（平面曲率，1.00）+s（径流强度指数，12.27）0.631*0.610 注：s（）表示 GAM 模型中的光滑项。括号内数字表示该光滑项的估计自由度，接近 1 时表明曲线形状接近直线，值越大表明曲线形状的波动越大。Note：s（）represents smooth terms within GAM formulae.Numbers in the bracket represent the estimated degree of freedom（EDF）of the smooth term.When EDF is close to 1，it indicates that the curve is close to a straight line.The larger the EDF，the more flexible the curve.图 2 SOC 的 GAM 模型中部分光滑项.图 a-c 的估计自由度分别为 6.07、11.34、12.73，见表 3；阴影部分表示光滑项估计的2 倍标准差范围.Fig.2 Some smooth terms of GAM for SOC.The degree of freedom in a-c are 6.07，11.34 and 12.73，respectively，as shown in Table 3.The shaded area represents two times of standard deviations.正相关关系，可能是因为：森林经营活动逐渐减少，使土壤遭受的侵蚀降低，因而 SOC 含量有所提升；随着海拔继续升高至 200 m260 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 低，在训练集上的调整决定系数较低。1000 土 壤 学 报 60 卷 http:/ 表 4 各土壤属性的 GLMM 模型 Table 4 GLMM models for soil properties 土壤属性 Soil property GLMM 方程 GLMM formulae 决定系数 R2 调整决定系数 Adjusted R2 SOC/（gkg1）19.99+s（高程，5.88）+s（坡度，7.62）+s（剖面曲率，5.13）+0.02 偏东度 0.422*0.399 pH 4.59+s（高程，3.64）+s（偏北度，4.87）+s（剖面曲率，6.27）+s（地形湿度指数，6.74）3.48 103 坡度 0.532*0.509 黏粒 Clay/%32.20+s（高程，3.68）+s（坡向，9.42）58.39 偏北指数 66.60 剖面曲率 0.328*0.302 CEC/（cmolkg1）16.24+s（坡度，13.14）+s（坡向，4.56）+s（剖面曲率，6.45）+s（径流强度指数，12.52）+41.68 平面曲率 0.682*0.667 注：s（）表示 GLMM 模型中的光滑项。括号内数字表示该光滑项的估计自由度，接近 1 时表明曲线形状接近直线，值越大表明曲线形状的波动越大。Note：s（）represents smooth terms within GLMM formulae.Numbers in the bracket represent the estimated degree of freedom（EDF）of the smooth term.When EDF is close to 1，it indicates that the curve is close to a straight line.The larger the EDF，the more flexible the curve.综上，对于所有土壤