基于大尺度低密度样点的东北土壤全氮空间插值方法比较

陈昌华，陈锡云，徐 英

（北京师范大学 地理学与遥感科学学院，北京100875）

土壤属性的空间变异性研究对于实现土壤资源的科学管理与高效利用至关重要。土壤属性空间分布图是土壤属性空间变异最直观的表达形式，然而结合土壤剖面数据与土壤类型图制作土壤属性分布图的传统方法带有较强主观性，难以定量反映土壤属性的空间变异。空间统计是基于数据所代表的属性具有空间自相关特性，根据已知空间数据来估计未知空间数据的数学方法［1］。20世纪70年代后期美国率先利用地理信息系统（GIS）的空间统计分析功能，开展了土壤普查、土壤属性数字制图及土壤变异性方面的研究，推动了利用地统计学方法研究土壤属性空间变异的快速发展，使得反距离加权、普通克里金等方法在该领域的应用趋于成熟。近年来，关于地统计插值方法比较及混合插值（hybrid interpolation）新模型［2］的发展和验证等方面的研究表明，以目标变量的显著相关性环境因子作为辅助信息（auxiliary information）的混合插值对提高不同区域的土壤属性插值精度具有良好效果［3］。常用的辅助信息有地形因子［4］、气候因子［5］及基于遥感卫星影像的土地利用类型［6］、土壤类型等多分类变量［7］。因辅助信息相对容易获得，在一定程度上降低了野外工作量、提高了空间插值的预测能力，使得利用有限野外采样点数据获得更高精度的土壤属性空间分布图成为可能。

我国在20世纪80年代以来也有不少土壤学学者开始这方面的研究，相关文献报道主要集中在小尺度下确定性内插和地统计插值方法比较［8－9］、土壤属性插值效果评价及未采样区土壤属性预测等方面。对于大尺度低密度采样的土壤普查数据，由于采样不均匀且样点分布较为稀疏，土壤属性空间变异的随机性因素增加，使类似方法的应用受限。但从理论上讲，进行空间插值的样点数据满足以下两点［10］即可：一是空间样点之间具有一定的相关性，二是空间变异具有连续性。因此，为充分发挥大尺度低密度采样下土壤普查数据的潜在价值，利用全国第二次土壤普查东北地区的土壤全氮含量（STNC）数据，比较不同插值方法下土壤属性的空间插值效果，同时利用土壤普查数据体系中显著相关的环境因子数据作为辅助信息，探索提高土壤普查数据空间插值精度的有效途径。为尽量避免土壤普查数据的采样设计不合理及数据内在变异性，比较了不同样本容量下插值结果以充分认识不同插值方法的空间趋势预测稳定性，希望能为今后土壤普查工作的土壤属性制图提供指导。

1 材料与方法

1.1 研究区域和数据来源

研究区域包括黑龙江、吉林、辽宁和内蒙古东部的赤峰市、通辽市、兴安盟和呼伦贝尔市（图1），位于115°52′—135°09′E，38°72′—53°55′N 之间，总面积约为1.24×106km2，是世界三大黑土带之一和我国广义的黑土分布区。土壤数据来源于中国第二次土壤普查资料（1979—1992年），排除典型剖面中缺少经纬度坐标和STNC值的记录，共417个有效样点（图1）。每个样点还整理了详细的地理位置、海拔高度、年均温、年降水量等环境特征数据和土种类型、所属土类、剖面厚度、pH、土壤阳离子交换量等土壤理化数据，形成完整的辅助信息因子数据库。

图1 研究区域和有效样点分布

1.2 数据处理

1.2.1 特异值处理与统计分析 由于特异值的存在会造成插值变量空间连续表面中断，采用“样本平均值3倍的标准差”来界定特异值，在该区间以外的数据均定为特异值，并分别用正常的最大值和最小值代替特异值［11］。利用SPSS 18.0软件对特异值处理后的STNC数据进行单样本K—S（Kolmogorov—Smirnov）检验和描述统计分析，以检验数据分布情况。最后基于ArcGIS 10.0软件平台形成土壤氮素空间属性数据库，并对STNC数据作对数转换计算，以满足空间插值需要。

1.2.2 空间插值 采用确定性插值和地统计插值方法中较为常用且具有代表性的反距离加权法（inverse distance weighting，IDW）、径向基函数（radial basis functions，RBF）、普 通 克 里 金 （ordinary kriging，OK）和回归克里金（regression kriging，RK）4种方法对STNC进行空间插值，前3种方法的相关报道较多［9－11］，其理论依据这里不再赘述。回归克里金（RK）方法是一种可使用连续变量又可使用类别变量作为辅助信息，同时结合Kriging插值方法来预测土壤属性的混合插值技术。可用方程式（1）—（2）表示［12］：

式中：f（x）——地理要素空间预测趋势项，用一个随辅助信息平稳变化的线性函数表示；a0，a1——系数项，使预测模型适宜任何偏离率；m（x）——影响空间预测趋势的辅助信息因子；Z（x）——土壤属性预测值；ε（x）——f（x）趋势预测土壤属性值的残差，模拟空间分布趋势及随机因素导致的不确定性。

本文采用趋势项f（x）值与残差项的Kriging预测值之和作为RK方法的预测结果，即STNC与环境变量之间的线性回归预测值和回归残差值的OK插值之和作为RK方法的预测值［13］。其中，参与线性回归分析的辅助信息因子共8个：经度（longitude，Lo）、纬度（latitude，La）、海拔（H）、年均温（mean annul temperature，MAT）、年均降水量（mean annul precipitation，MAP）、土层厚度（d）、pH 值和阳离子交换容量（CEC）。

参考目前研究较为成熟的县域尺度空间插值模型及本研究中不同样本容量下样点尺度的变化，在ArcGIS 10.0和 GS＋Version 7.0中，对4种空间插值方法在大尺度低密度采样下的函数模型和参数设定详见表1。OK和RK方法中涉及到Kriging插值均首先在GS＋软件中进行相应变量的半方差分析，再基于以上所得最优半方差模型参数，在ArcGIS软件平台下完成空间插值。相关研究表明：通过软件得到的各向同性半方差模型实际上是对各向异性半方差函数模型均化的结果［11］，其差别不会对Kriging插值产生明显影响，因此在OK和RK半方差分析中均采用各向同性半方差模型。IDW和RBF两种方法由于拟合模型参数调试相对方便，直接在ArcGIS平台上完成。

表1 不同插值方法模型基本参数设定情况

1.2.3 样本容量设定 为充分认识土壤调查数据STNC的空间变异性特征，以110［14］为最小个数设计了7个不同容量的样点系列。具体操作为：利用Arc－GIS中的Create Subsets功能，随机抽取417条样点记录中80%的记录作为样本A1，抽取A1中80%的记录作为样本A2，从A2中抽取80%的记录作为样本A3，…，依次按比例抽取形成6个不同容量的样本数据，再加上全部的417条样点记录，共形成7个样本系列。由于STNC原始数据中部分参与回归分析的环境因子数据缺失，RK方法初始样本容量为337，并单独记为样本B，之后从A3开始计算，尽量保证4种空间插值方法使用统一样点系列的STNC值进行空间趋势预测。

1.2.4 插值结果精度检验与比较 选择交叉验证方法中的平均误差（ME，mean error）、均方根误差（RMSE，root mean square error）和预测值与观测值之间的相关系数R作为不同空间插值模型精度检验指标。为了消除不同样本容量下STNC平均值和标准差的波动性，将不同样本容量下计算的ME和RMSE分别除以各自的平均值（mean）和标准差（SD）作标准化处理。标准化后的ME／mean和RMSE／SD值越小，R值越大则说明空间插值精度越高。同时参照相关文献［15］引入相对提高指数（RI，Relative Improvement）来定量化表示RK方法对RMSE精度的

式中：RMSE′——标准化 RMSE，即对应的 RMSE／SD值；i——IDW，RBF，OK方法。

RI越大则说明RK方法的精度提高程度越高，即辅助信息的添加越有效。提高度，其计算公式为：

2 结果与讨论

2.1 STNC统计特征

417条有效样点数据的统计分析（表2）表明，STNC原始取值范围为0.08～21.48g／kg，剔除特异值后的取值范围为0.08～5.31g／kg。

剔出特异值前后STNC变异系数均大于80%，表现出较强的变异性，这与研究区域自然条件复杂，土壤类型和土地利用类型多样有关。但剔除特异值降低了数据的偏度和峰度，且使偏度和峰度更为接近。剔除特异值后的STNC进行单样本K—S检验表明，对数转换之前STNC属于偏正态分布（p＜0.001）；经过对数转换后STNC的Z值由对数转换之前的4.12减小为1.03，且 P 值变为0.244（＞0.01）。说明经对数转换后STNC数据呈正态分布，空间趋势预测中需要对有效样点的STNC数据进行对数转换，以保证空间插值数据满足正态分布。不同样点系列之间，随样本容量减小，样点最大间距减小，而最小间距的变化幅度不大（表2）。样本系列A5以后最大采样间距变化较大，但空间插值范围仍保持A1时的90%以上，对STNC空间变异性在大尺度上的分辨不会存在影响。不同样本容量下STNC基本描述统计量波动不大，偏度系数在2.57～2.89之间；峰度系数在7.86～9.73之间，说明各样本系列具有相似的统计特征，不同容量的样点系列设计比较合理。空间插值过程中不再针对每个样本进行正态分布检验，空间预测完成后再将对数转换后的数据计算结果进行还原。

表2 STNC数据描述性统计特征及不同样本系列数据特征比较

2.2 辅助信息因子的分析与筛查

根据相关性分析和回归分析结果得到STNC回归预测值及回归残差值，是采用RK方法分离趋势项和残差项的重要计算过程。由于进行逐步回归时因子变量间的共线性会影响结果的真实性，首先对8个环境变量进行共线性判别（表3），以容忍度（tolerance，T）大于0.1和因子膨胀系数（variance inflation factor，VIF）小于10作为因子非共线性的判断标准［16］。

结果表明，8个因子之间均满足非共线性，可以全部作为自变量参与到STNC的逐步线性回归分析中。通过相关性分析得到STNC与Lo，La，H，CEC呈显著正相关，与MAT，d，pH呈显著负相关，而与 MAP相关性不显著。这可能是由于人类耕作过程中灌溉措施在一定程度上降低了降水对土壤养分的影响。

表3 STNC与影响因子的相关性及变量共线性诊断分析

通过逐步线性回归，设定因子变量进入方程的概率水平为0.05，剔除的概率水平为0.10，得到STNC预测的最优因子为CEC，MAT，d和pH，最优线性回归方程式如公式（4）。

考虑到本研究空间尺度较大且自然因素和土地利用变化较为复杂，同时对比其他研究尺度下RK方法回归方程的决定系数［3，15，17］，该回归拟合方 程的决定系数（R2＝0.328）处于可接受范围。以此计算得到的STNC回归预测残差值呈近似正态分布（图2），故用OK方法进行回归预测残差值的空间插值时不再做数据转换处理。

图2 STNC回归预测值的残差分布情况

2.3 不同样本容量下STNC空间变异性分析

利用最优半方差函数模型定量化描述土壤属性空间变异性，是利用地统计方法进行土壤属性空间插值的基础。通过比较不同模型拟合的R2值和RSS值大小，STNC除样本A3满足球状模型（Sh）外（表4），其他样本容量下均满足高斯模型（Gs），说明样本容量变化对拟合空间结构的半方差函数模型影响不大。在实际参数设定过程中，有效滞后距和步长的乘积为最大采样间距的1／2，由于样本容量的减少对样点最大间距的影响，有效滞后距有所下降。块金效应（Nugget／Sill，N／S）是耕作措施、种植制度等随机因素引起的空间异质性占总变异的比例，当样本容量≥171时，STNC的 N／S 值为31.81%～49.93%，表现出中等水平的空间自相关性，且随着样本容量的减少N／S值基本呈现上升趋势，说明随机性因素影响比重上升造成了空间自相关性的减弱。当样本容量＜171时，N／S值均小于25%，此时半方差模型对空间变异性的辨识能力下降使随机性因素被掩盖，空间变异结构发生变化，得到的空间分析结果难以确信。STNC回归预测残差值的N／S取值范围基本落在STNC块金效应波动范围以内，说明回归预测残差值保留了STNC固有空间结构。STNC变程为1 643.72～3 635.57km，当样本容量≥171时，STNC样点最大间距均落在空间自相关范围内，此时各样本容量下的平均变程为2 040.53km。不同样本容量下STNC回归预测残差值的空间自相关范围也均大于2 000km。

表4 STNC半方差模型及拟合参数

目前针对大尺度低密度采样的土壤属性空间插值的报道很少，对比东北地区已有的相关研究，在小流域尺度（992hm2）上292个样点的STNC空间插值的N／S块金效应为89.2%，变程为1 899m［18］。县域尺度上，张素梅等［3］对自采的吉林省农安县471个样点的STNC空间插值的N／S值为36.7%，变程为109.3km。Liu等［19］利用吉林省德惠县的354个土壤肥力普查样点数据所作的SOCC空间插值的N／S值为16.1%，变程达到632km。区域尺度上，Zhang等［20］利用995个自采样点对黑龙江和吉林两省黑土区的STNC空间插值的N／S值为5.97%，变程为1 037.5km。上述结果不难发现，不同空间尺度下不同样本容量的土壤碳或氮的块金效应较为一致地表现出中强度空间相关性，且变程随着各研究尺度的增大而增大。本研究区域为整个东北地区，研究数据的样点间距远大于上述同类研究，得到的空间自相关范围也较大。

2.4 不同方法的空间插值结果

比较4种插值方法的空间趋势预测结果（图3），整体来看各空间插值方法均能较直观的反映出该区STNC空间分布的基本态势，STNC均呈现从东北向西南逐渐递减规律。高值区（＞1.10g／kg）呈倒“V”形条带状贯穿于研究区域的东北部地区，包括黑龙江省大部分地区和内蒙古东北部的呼伦贝尔市及三江平原地区。所对应的基本是暗棕壤、漂灰土、黑钙土林区及沼泽湿地等理论上均有利于有机质积累的区域。西北部中值区（0.63～1.10g／kg）包括松嫩平原在内的吉林省和内蒙古的兴安盟、赤峰市以及通辽市的西北部地区。该区域主要为农业耕作区和草甸草原，其中受人类活动影响的松嫩平原多发育一些草甸白浆土、白浆化暗棕壤等过渡土壤类型，STNC水平适中。低值区（＜0.63g／kg）则呈斑块状分布在通辽市的东南部和辽宁省的西北部沿海地区，是沙漠、旱地和水田区，土壤类型以风沙土、潮土和滨海盐土等有机质积累较弱的土壤类型为主。此插值结果分布趋势与Yang等［21］使用全国第二次土壤普查数据得到的土壤有机碳密度的空间插值结果在该区的分布态势保持很好的一致性。

图3 STNC不同空间插值方法结果比较

对比不同样本容量下插值结果（A1，A4，A73个样本系列，图3）可以发现，随着样本容量降低，4种方法插值结果的中、低值区范围扩大，高值区缩小且有了明显的带状化表现，尤其是采样点比较稀疏的东北部三江平原地区、内蒙古东北和黑龙江西北部，过分的带状化削弱了空间变化特征。在空间分布的细节表达上，相同样本容量下IDW和RBF插值图中出现了特异值斑块或“牛眼”现象，造成空间趋势表达不连续，OK和RK方法由于Kriging插值模型对空间特异性具有一定的平滑作用，STNC空间变化规律明显。同一种空间插值方法下，随着样本容量的减少，整体变化规律一致但趋势变化的连续性变差。对比RK方法下A1和A4两个样点系列的插值含量水平和斑块面积变化，在样本容量较少（A4）情况下RK方法在空间趋势预测上仍较为稳定并在西北部和东北部的样点稀疏区体现出更多的细节变化，说明辅助信息的加入对这些区域的空间预测具有一定帮助。

2.5 不同空间插值方法的精度分析

插值结果的精度比较是评价不同插值方法有效性与适用性的重要手段，不同插值方法对土壤属性插值精度在不同地区相关研究中的表现不一，部分研究表明Kriging的效果优于IDW［22］和RBF中的Spline函数模型［23］，也有研究表明IDW和Kriging插值方法均有较好的表现［8］。仅针对不同Kriging插值方法的研究多体现为RK方法或其他添加辅助信息的Kriging插值方法和模型在提高插值精度方面表现较好［24］。比较本研究中4种空间插值方法的标准化ME，标准化RMSE，R三个精度指标可见（图4），其空间预测效果顺序为：RK＞OK＞RBF＞IDW。不同空间插值方法下3个精度指标值随着样本容量的变化规律不明显，计算不同样本容量下各精度指标平均值，得到RK方法的ME值同OK，RBF和IDW方法相比平均降低0.18，0.17和0.16，RMSE值平均分别降低了0.15，0.18和0.19，R 平均分别提高了0.40，0.42和0.43。

从插值精度的变化幅度来看，随着样本容量的变化，RK和OK方法对不同样本系列表达效果较为稳定。由上述分析可知，当样本容量＜171时，STNC空间变异结构发生变化，此时A6，A7两个样本容量下反映出的精度变化水平也不可信。

图4 不同方法STNC空间插值精度随样本容量的变化

通过计算各样本容量的平均RI值（图4），RK方法相对于其他空间插值方法将预测精度分别提高了16.15%（OK），18.50%（RBF）和19.40%（IDW），表明利用环境因子作为辅助信息使该区空间插值精度得到了有效提高。相关研究中，Simbahan［4］等利用地表反射率、相对高程及土壤电导率3个环境因子作为辅助性的RK方法提高了CS元素制图精度的16%。Zhang［7］等利用水流功率指数、湿润指数及坡度3个环境变量作为辅助信息的RK方法对土壤有机质空间插值提高了13.93%的精度。Sumfleth［15］等利用相对高程和河网高程2个环境变量作为辅助信息的RK方法对STNC空间预测精度提高了13%。本研究中使用的辅助信息是土壤普查工作中必须测定的项目或者是相应土壤点位上配套的气象站点数据，同上述研究相比，在辅助信息易获得性以及与土壤属性的时空匹配性上具有明显的优势，这对于今后提高土壤普查数据的空间插值精度价值巨大。

3 结论

（1）我国第二次土壤普查期间东北地区STNC在0.08～21.48g／kg之间，剔除特异值前后STNC数值水平具有较大变异性。大尺度低密度采样下STNC表现出中等水平空间自相关性，但变程大于同区域的小尺度研究结果。目前针对大尺度（＞1∶100 000）低密度采样的土壤属性空间插值方法的研究较少，本研究得到的STNC空间变异性特征只能与同区域小尺度上的相关研究进行类比分析，得到的变程值对于STNC空间自相关范围的指示意义难以进行有效评价，但空间插值方法是利用土壤普查数据进行土壤属性制图的有效途径，因此对于大尺度低密度采样的土壤属性空间插值研究需要进一步的加强。

（2）添加了辅助信息的回归Kriging（RK）是对经度（Lo）、纬度（La）、海拔（H）、年均温（MAT）、年降水量（MAP）、土层厚度（d）、pH 值、土壤阳离子交换量（CEC）8个环境变量与STNC值进行线性回归拟合，将得到的回归拟合残差值运用普通Kriging（OK）进行OK空间插值模拟，由OK模拟后的残差值与回归拟合的预测值之和得到RK方法的空间预测值。共线性判别表明，8个因子之间均满足非共线性，可以全部作为自变量参与到STNC的逐步线性回归分析中；通过逐步线性回归得到本研究RK方法的辅助信息为CEC，MAT，d和pH值，其中，STNC与CEC呈显著正相关，与 MAT，d，pH 值呈显著负相关关系。

（3）样本容量变化对不同插值方法的STNC空间趋势预测影响不明显。但在细节表达上，相同样本容量下OK和RK方法的STNC空间插值结果空间变化特征更鲜明。随着样本容量的减少，同一种空间插值方法的空间表达更简略。RK方法对无样点区的空间趋势预测体现出更多的细节变化。此外，在以后的土壤普查工作中，应该着重加强本研究中氮含量插值结果的高、中、低值区的过渡地带以及内蒙古东北部的风沙地、黑龙江东部与东北部的多年冻土区和吉林东南沿海的盐渍地等生态敏感地区的布点设计。

（4）不同样本容量下4种方法空间插值精度总体表现为：RK＞OK＞RBF＞IDW，大尺度低密度采样下地统计插值方法要优于确定性内插。RK方法相对于其他空间插值方法将预测精度提高程度（RI）分别为 16.15%（OK），18.50%（RBF）和 19.40%（IDW）。土壤普查过程中所收集的环境因子作为辅助信息即可使该区土壤氮含量的空间插值精度得到有效地提高，而且在低样本容量下RK方法在空间趋势预测上仍较为稳定，预示着辅助信息的添加对大尺度低密度采样下STNC的空间表达具有较好的效果。

［1］ Lam N S.Spatial Interpolation［M］∥International Encyclopedia of Human Geography.Elsevier Science Ltd.，2009.

［2］ Takata Y，Funakawa S，Akshalov K，et al.Spatial prediction of soil organic matter in northern Kazakhstan based on topographic and vegetation information［J］.Soil Science and Plant Nutrition，2007，53（3）：289-299.

［3］ 张素梅，王宗明，张柏，等.利用地形和遥感数据预测土壤养分空间分布［J］.农业工程学报，2010，26（5）：188-194.

［4］ Simbahan G C，Dobermann A，Goovaerts P，et al.Fineresolution mapping of soil organic carbon based on multivariate secondary data［J］.Geoderma，2006，132（3／4）：471-489.

［5］ McKenzie N J，Ryan P J.Spatial prediction of soil properties using environmental correlation［J］.Geoderma，1999，89（1／2）：67-94.

［6］ Bui E N，Moran C J.A strategy to fill gaps in soil survey over large spatial extents an example from the Murray-Darling basin of Australia［J］.Geoderma，2003，111（1／2）：21-44.

［7］ Zhang Shiwen，Huang Yuanfang，Shen Chonggan，et al.Spatial prediction of soil organic matter using terrain indices and categorical variables as auxiliary information［J］.Geoderma，2012，171（S）：35-43.

［8］ Schloeder C A，Zimmerman N E，Jacobs M J.Comparison of methods for interpolating soil properties using limited data［J］.Soil Science Society of America Journal，2001，65（2）：470-479.

［9］ 石小华，杨联安，张蕾.土壤速效钾养分含量空间插值方法比较研究［J］.水土保持学报，2006，20（2）：68-72.

［10］ Webster R.Soil resources and their assessment［J］.Phil.Trans.R.London：B，1997，352（1356）：963-973.

［11］ 刘付程，史学正，于东升，等.基于地统计学和GIS的太湖典型地区土壤属性制图研究：以土壤全氮制图为例［J］.土壤学报，2004，41（1）：20-27.

［12］ Stacey K F，Lark R M，Whitmore A P，et al.Using a process model and regression Kriging to improve predictions of nitrous oxide emissions from soil［J］.Geo-derma，2006，135：107-117.

［13］ Odeh I O A，McBratney A B，Chittleborough D J.Further results on prediction of soil properties from terrain attributes：Heterotopic Cokriging and regression-Kriging［J］.Geoderma，1995，67（3-4）：215-226.

［14］ Webster R，Oliver M A.Sample adequately to estimate variograms of soil properties［J］.Journal of Soil Science，1992，43（1）：177-192.

［15］ Sumfleth K，Duttmann R.Prediction of soil property distribution in paddy soil landscapes using terrain data and satellite information as indicators［J］.Ecological Indicators，2008，8（5）：485-501.

［16］ Menard S W.Applied logistic regression analysis［M］.Thousand Oaks：CA：SAGE Publication，Inc.，1995.

［17］ Ziadat F M.Analyzing digital terrain attributes to predict soil attributes for a relatively area［J］.Soil Science Society of America Journal，2005，69（5）：1590-1599.

［18］ Wei Jianbing，Xiao Duning，Zhang Xingyi，et al.Topography and land use effects on the spatial variation of soil organic carbon：A case study in a tipical small watershed of the black soil region in Northeast China［J］.Euras Soil Sci.，2008，41（1）：39-47.

［19］ Liu Dianwei，Wang Zongming，Zhang Bai，et al.Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region，Northeast China［J］.Agriculture，Ecosystems and Environment，2006，113（1／4）：73-81.

［20］ Zhang Xingyi，Sui Yueyu，Zhang Xudong，et al.Spatial variability of nutrient properties in black soil of Northeast China［J］.Pedosphere，2007，17（1）：19-29.

［21］ Yang Yuanhe，Mohammat A，Feng Jianming，et al.Storage，patterns and environmental controls of soil organic carbon in China［J］.Biogeochemistry，2007，84（2）：131-141.

［22］ Kravchenko A，Bullock D G.A comparative study of interpolation methods for mapping soil properties［J］.Agronomy Journal，1999，91（3）：393-400.

［23］ Panagopoulos T，Jesus J，Antunes M D C，et al.A-nalysis of spatial interpolation for optimising management of a salinized field cultivated with lettuce［J］.European Journal of Agronomy，2006，24（1）：1-10.

［24］ Hengl T，Heuvelink G，Alfred S.A generic framework for spatial prediction of soil variables based on regression-Kriging［J］.Geoderma，2004，120（1／2）：75-93.