基于协同克里格的耕层土壤速效钾空间异质性研究

江叶枫，郭 熙

(江西农业大学 江西省鄱阳湖流域农业资源与生态重点实验室，江西 南昌 330045)

土壤是有限且不可再生的自然资源，在提供生态系统服务方面发挥着至关重要的作用[1-2]。为维护自然资源以支持未来的可持续发展，亟须改善和维护土壤的生态系统服务功能[3]。土壤属性的变化与土壤生态系统服务功能的可持续性息息相关，土壤侵蚀、土壤酸化、土壤重金属污染与土壤养分供需不平衡等，都不可避免地会影响土壤质量和生态系统服务功能，降低土壤对人类的支持能力[3-5]。土壤养分是土壤可持续利用和管理的重要基础[3-4]。土壤速效钾(AK)能够直观反映土壤钾素水平和供应能力，作为作物生长周期内可以获取的主要养分，直接影响农产品产量和品质[6]。近年来，由于各地区生产水平、耕作技术措施、种植制度和耕地资源利用强度不同，土壤速效钾的空间异质性普遍存在[7]。随着精准农业和配方施肥技术的发展，土壤速效钾的空间异质性研究不断深入[8]。准确掌握土壤速效钾的空间异质性对科学制定配方施肥方案、保护土壤、实现土壤的可持续利用等具有重要意义[6-8]。

土壤是连续的，科学家难以检测到每一处土壤的属性值，不得不依赖于空间插值技术[9]。国内外学者主要运用克里格空间插值法对土壤属性进行空间插值并研究土壤属性的异质性[10-11]。该方法基于半变异函数，是对土壤属性进行的最优无偏估测，能有效揭示土壤属性空间分布和变异特征。其中，普通克里格法(ordinary Kriging，OK)只需要考虑目标变量，简单易行，应用广泛，但插值结果在土壤属性异质性描述方面受采样数量和密度的影响[7，12]。协同克里格法(co-Kriging，COK)是一种基于目标变量与辅助变量的半变异函数，由于辅助变量与土壤属性相关性较高，因而以这种方式生成的可视化图通常能够更加高效、客观、真实地描述土壤属性空间变异特征，在土壤属性空间变异研究中越发受到研究人员的青睐[13-14]。

COK较OK能提高空间插值精度，但在土壤属性的异质性描述方面是否也较OK存在一定优势尚不清楚[13-14]。当下，大部分研究对土壤属性的异质性描述仍处于定性分析层面，缺乏对影响因子作用的量化分析[7]。从相关文献来看，关于土壤性质变异的影响因素定量分析报道仍存在明显不足：1)对土壤属性的辅助变量考虑不充分，地形因子因能较好地刻画土壤性质的地带性分异规律而受到较大关注[12，15]，但仅少量学者定量分析了成土母质和土壤类型等对土壤属性空间异质性的影响[11，16]；2)相关的定量研究较多针对土壤多种性质或针对有机碳[17]、氮素[4]和磷素[18]等单一性质，鲜见关于土壤钾素的报道，特别是在我国耕地复种指数增高和土壤钾素供应严峻的背景下，对速效钾变异的影响因素缺乏量化认识[7]；3)不同区域COK辅助变量的选取各有不同，前人曾对川中丘陵[19]、黄土高原[20]、西北山区[21]和城乡交错带[22]等进行研究，而鄱阳湖平原区则鲜见报道。

为此，本研究以鄱阳湖平原区典型县——万年县为案例区，运用测土配方施肥采集的耕层(0～20 cm)土壤样品数据，以土壤速效钾为研究对象，量化土壤速效钾空间变异的影响因素，运用OK和COK对土壤速效钾的空间异质性进行研究。结果可为类似地区土壤钾素调控、土地可持续管理和利用，以及土壤属性空间分布预测模型构建提供参考。

1 材料与方法

1.1 研究区概况

研究区位于万年县，属鄱阳湖平原区，地理坐标介于116°46′—117°15′E、28°30′—28°54′N之间，总面积1 140.76 km2，东西长47 km，南北宽43 km。研究区属亚热带季风性气候，地貌类型以平原为主，呈“西北部高、东南部低”的总体趋势(图1-a)，海拔12～651 m。土类主要是红壤和水稻土，并伴随有极少量的潮土、石灰土和紫色土，亚类主要为红壤和潴育型水稻土，土属主要包括黄泥红壤和鳝泥红壤，潴育型潮砂泥田、潴育型黄泥田和潴育型鳝泥田。成土母质包括紫色泥页岩类风化物、第四纪红色黏土、河流冲积物、泥质岩类风化物、酸性结晶岩类风化物、碳酸岩类风化物和紫色砂砾岩类风化物。境内水资源比较丰富，河流主要包括乐安河、珠溪河和万年河等182条，总长806 km，河网密度707 m·km-2，年均降水量1 766 mm，年均无霜期263 d，年均相对湿度82%。土地利用分布格局为“六山一水二分田，一分道路和庄园”，耕地面积约为230 km2，耕层土壤pH介于4.5～7.0。万年县是世界稻作文化的发源地，也是人工栽培稻起源地和贡米的原产地，被称为“稻米之乡”，传统的稻米习俗已在该地传承了几千年，万年贡米是晚籼稻“坞源早”品种加工而成的产品，原只产于万年县裴梅镇部分乡村，后经种植推广，截至2016年，万年县贡米种植面积达170 km2。

1.2 土壤采样与数据处理

于2014年10—11月在作物收割后，根据《测土配方施肥技术模式》，在充分考虑土壤类型、成土母质和地形条件的基础上，遵循均匀性、代表性和连续性的原则，运用“S”形采样策略和多点混合的方法采集耕层(0～20 cm)土壤样品，每个样点采集1 kg土样，4个重复样本在该点周围5 m范围内采集，并运用GPS详细记录该点的经纬度、成土母质、土壤类型和海拔等信息，共获取100个耕层土壤样品。土壤样品经自然风干后，带回实验室磨碎过筛，采用乙酸铵浸提—原子吸收分光光度法测定土壤速效钾含量[7]。数字高程模型(digital elevation model， DEM)如图1-a所示，分辨率30 m，通过地理空间数据云(http://www.gscloud.cn/)获取。

受采样及化学分析误差的影响，土壤速效钾含量测量结果存在离群值。采用拉依达准则法[4]对测定数据进行检验，剔除离群值，得到剔除后有效样点96个(图1-b)。高程、坡度、坡向、曲率、坡度变率和坡向变率均由数字高程模型通过ArcGIS 10.2表面分析模块进行分析和提取，各地形因子的计算公式详见文献[23]。土壤AK根据第二次土壤普查土壤养分分析标准分为6级，1～6级分别对应极高(>200 mg·kg-1)、丰富(>150～200 mg·kg-1)、中等(>100～150 mg·kg-1)、较低(>50～100 mg·kg-1)、低(>30～50 mg·kg-1)和极低(≤30 mg·kg-1)。成土母质和土壤类型为定性变量，采用哑变量[4，16]进行赋值。

1.3 模拟精度评价

采用交叉验证评价OK和COK的预测精度。交叉验证利用每个实测点周围点对该点进行预测，将该点预测值与实测值进行比较。以均方根误差(RMSE)、平均绝对误差(MAE)、平均相对误差(MRE)对预测值与实际采样值进行对比分析，得出精度评价结果，计算公式如下：

(1)

(2)

(3)

图1 研究区数字高程模型和样点分布图Fig.1 Maps of DEM and sample points in study area

描述性统计分析、相关性分析和回归分析在IBM SPSS statistics 22.0中实现，半变异函数分析在GS+9.0中实现，OK和COK插值通过ArcGIS 10.2实现。

2 结果与分析

2.1 描述性统计分析

2.1.1 土壤速效钾

研究区土壤速效钾值域范围为33.46～164.84 mg·kg-1，均值为82.04 mg·kg-1。根据第二次土壤普查土壤养分分级标准，研究区速效钾含量为4级，属于较低水平，勉强能满足当地耕地生产需要。因部分地区土壤钾素盈余或不足，以及耕作技术措施差异，同时考虑耕地复种指数和化肥施用比例不协调等原因，建议重视土壤钾素保育[7]。研究区土壤速效钾的最大值与最小值之比为4.93，表明研究区土壤钾素分布极不平衡，变异系数为30.49%，属于中等变异。Kolmogorov-Smirnov(K-S)检验的结果表明，土壤速效钾含量符合正态分布。

2.1.2 土壤类型与成土母质

土壤类型反映区域条件的综合变化，因成土过程、矿物组成和发育程度不同，土壤速效钾存在空间异质性[7，16]。由表1可知，不同类型土壤速效钾均值差异显著(P<0.05)，表现出水稻土>红壤的总体趋势。这主要是由于水稻土受长期人为活动施钾肥的影响[4]，导致其速效钾均值要高于红壤。各土属中，以潴育型潮砂泥田均值最高，达到100.88 mg·kg-1，鳝泥红壤最低，为73.95 mg·kg-1。

成土母质的矿物分解是土壤钾素的基本来源，不同成土母质因土壤团聚体数量及其稳定性、物理化学组成和风化淋溶进程等引起土壤速效钾变异[24]。从表2可以看出，成土母质对土壤速效钾有显著(P<0.05)影响。不同母质间土壤速效钾含量为61.10～92.04 mg·kg-1，其中，以白垩纪紫色泥页岩类风化物发育而来的土壤速效钾均值含量最高，而紫色砂砾岩类风化物最低。白垩纪紫色泥页岩一般含碳酸钙，呈中性或微碱性反应，有机质含量低，但磷、钾丰富；酸性结晶岩类风化物以残积物和坡积物为主，质地主要为黏壤土，土壤保水保肥能力强，土壤速效钾含量也相对较高；紫色砂砾岩类风化物易风化，水土流失严重，因此土壤速效钾较低[7，24]。其他研究也得到类似结论[16]。

表1不同土壤类型速效钾描述性统计特征

Table1Descriptive statistic characteristics of soil available potassium relative to soil types

土类Soil group亚类Subgroup土属Soil family样点数Sample No.均值Mean/(mg·kg-1)SD/(mg·kg-1)CV/%红壤Red soil红壤Red soil黄泥红壤Yellow mud red soil584.84 b27.0131.84鳝泥红壤Muddy red soil3673.95 c18.3624.83水稻土潴育型水稻土潴育型潮砂泥田Breeding type tidal mud field8100.88 a29.7629.50Paddy soilWaterloggogenic潴育型黄泥田Breeding yellow mud field1489.61 b20.8023.21type paddy soil潴育型鳝泥田Breeding muddy field3382.66 b27.5233.29

SD，标准差；CV，变异系数。同列数据后无相同字母的表示差异显著(P<0.05)。下同。

SD， Standard deviation; CV， Coefficient variation. Data marked without the same letters indicated significant difference atP<0.05. The same as below.

表2不同成土母质土壤速效钾描述性统计特征

Table2Descriptive statistic characteristics of soil available potassium relative to parent materials

成土母质Parent material样点数Sample No.均值Mean/(mg·kg-1)SD/(mg·kg-1)CV/%紫色泥页岩类风化物Purple mud shale weathering1492.04 a29.6232.18第四纪红色黏土Quaternary red clay2989.72 a26.2029.20河流冲积物River alluvial1081.78 b14.6417.90泥质岩类风化物Mudstone weathering1269.57 b13.6919.68酸性结晶岩类风化物Acidic crystalline rock weathering589.74 a19.9022.18碳酸岩类风化物Carbonate weathering2175.15 b24.3232.36紫色砂砾岩类风化物Purple glutenite weathering561.10 c10.8917.82

2.1.3 土壤速效钾与影响因素的相关性分析

从表3可以看出，土壤速效钾与pH的相关性达极显著水平(P<0.01)，相关系数为-0.258，说明在一定范围内土壤速效钾随pH降低而逐渐增加，这与其他区域研究结果一致[25]。土壤阳离子交换量(CEC)与速效钾相关性同样达极显著水平(P<0.01)，且相关系数达到0.602，表明CEC与AK呈较强的正相关关系，CEC越高，土壤速效钾越容易累积。高程与土壤速效钾的相关性达极显著水平(P<0.01)，相关系数为-0.414，高程越高处土壤受到的冲刷侵蚀越严重，土壤速效钾越易流失，而地势低洼处速效钾易随地表物质积聚[7]。土壤速效钾与其他地形因子，如坡度、坡向、曲率、坡度变率、坡向变率的相关性不显著，可能与鄱阳湖平原区有关，地形因子本身的变异较小。同时，本研究获取的DEM数据精度较低，其派生的地形变量精度也相对较低，这也会在一定程度上影响地形因子与土壤速效钾的相关性[16]。

表3土壤速效钾与pH、CEC和地形因子的Pearson相关系数

Table3Pearson correlation coefficients of soil available potassium with pH， CEC， and terrain factors

指标IndexrpH-0.258**CEC0.602**高程Elevation-0.414**坡度Slope-0.098坡向Aspect-0.072曲率Curvature0.18坡度变率Slope variability-0.073坡向变率Aspect variability-0.06

*，P<0.05; **，P<0.01.

2.1.4 土壤速效钾与影响因素的回归分析

为定量揭示成土母质、土壤类型(亚类和土属)、高程、pH和CEC对研究区耕层土壤速效钾空间变异的影响，对影响因素进行回归分析(表4)。结果表明：成土母质、土壤类型(亚类和土属)、高程、pH和CEC对土壤速效钾空间变异影响均达极显著水平(P<0.01)。成土母质对土壤速效钾空间变异的影响程度较低，为8.3%，这与相关研究结果[16]一致。江叶枫等[12，16]研究表明，该区主要成土母质包括白垩纪紫色泥页岩、第四纪红色黏土和几种风化物，多由花岗岩、石灰岩和紫色砂岩经风化而成，岩性总体上较为一致，导致土壤机械组成相近[16]，成土母质主要通过影响土壤机械组成来影响土壤速效钾分异，因此，成土母质对本研究区土壤速效钾影响程度较低。土属的影响程度要高于土类和亚类。江叶枫等[26]认为，由于土属较土类和亚类而言反映的成土过程、化学成分和生物活性等信息更多，因此影响程度更高。研究区pH介于4.5～7.0，当pH值降低时，土壤胶体微粒表面所负电荷也减少，导致K+吸附量随之增加。陈洋等[7]研究表明，在酸性环境中K+运移量较少，固定量增多。CEC基本上代表了土壤可能保持的养分数量，其影响因素包括土壤胶体类型、土壤质地和pH值，这些都是影响土壤速效钾变异的因素。因此，CEC对土壤速效钾变异影响程度较高。

2.2 土壤速效钾的半变异函数分析

表5给出了土壤速效钾半变异函数拟合结果。

表4不同因素对土壤速效钾的回归分析

Table4Regression analysis of soil available potassium with affecting factors

影响因素FactorF决定系数R2调整决定系数Adjusted R2P成土母质 Parent material1.9130.0920.083<0.01土壤类型 Soil group土类 Soil group5.4180.0710.066<0.01亚类 Subgroup5.4180.0710.066<0.01土属 Soil family6.2980.0850.080<0.01高程 Elevation12.6900.1690.165<0.01pH4.9060.0650.059<0.01CEC30.3760.3560.354<0.01

表5土壤速效钾的半变异函数值

Table5Semivariance parameters for soil available potassium

土壤属性Soil properties模型Model块金值Nugget基台值Partial sill块金效应Ratio of nuggest to sill/%变程Range/km决定系数Determination coefficient残差Residual速效钾AK指数Exponential260.54560.0046.532.640.6581.37×10-4

基于最大的拟合系数和最小的残差得到土壤速效钾的最佳拟合模型为指数模型，块金值为260.54，说明存在由实验误差和田间采样等人为因素造成的空间变异，块金效应值为46.53%，表明速效钾的空间变异受自然特征和人为活动的共同影响，为中等程度的空间相关性[16-20]。

2.3 土壤速效钾空间分布模拟

为直观反映研究区土壤速效钾的空间异质性，在半变异函数的基础上运用OK、COK1、COK2和COK3对研究区土壤速效钾进行空间插值，用可视化图作为手段来评估4种方法描述土壤速效钾空间异质性的能力。其中，COK1代表以CEC作为辅助变量进行插值，COK2代表以CEC和pH作为辅助变量进行插值，COK3代表以CEC、pH和高程同时作为辅助变量进行插值。

交叉验证的精度评价结果(表6)表明，基于辅助变量的3种COK明显优于仅基于目标变量的OK，COK1、COK2和COK3较OK的RMSE分别降低了1.03、1.92、4.86 mg·kg-1，MAE分别降低了1.10、2.21、5.37 mg·kg-1，MRE分别降低了1.41、2.74、5.50个百分点。基于辅助变量的COK预测精度得到了较为明显的提高。从表6还可以看出，随着辅助变量增加，COK的精度也在相应提升。这一方面表明增加与目标变量相关的辅助变量协助空间插值可以提高精度，另一方面也说明运用多个辅助变量较单个辅助变量精度更高，后期研究中可以适当增加辅助变量个数以提高预测精度。

从图2可以看出，OK、COK1、COK2和COK3对研究区土壤速效钾的模拟均表现出“西北部高、东南部低”的总体趋势，这与数字高程模型的空间变化较吻合，比较符合地学分布规律。从空间分布模拟效果看，4种方法模拟的局部特征差异明显。OK得到的空间分布模拟结果较平滑，高低值界限较清晰，难以准确地表达土壤速效钾的空间异质性，预测值域范围58.83～121.83 mg·kg-1，与统计分析值有较大差距。COK1运用相关性最强的CEC作为辅助变量进行空间插值，预测精度较OK有所提升，但空间分布模拟效果提升不明显，从预测范围(63.59～109.48 mg·kg-1)和均值(78.59 mg·kg-1)来看，COK1预测更趋向于均值，平滑效应明显，难以较清晰地刻画土壤速效钾的空间异质性信息。COK2在COK1的基础上引入pH作为辅助变量进行协同插值，不仅预测精度有所提升，在空间异质性信息描述方面也有很大提升。COK2预测土壤速效钾值域介于58.19～122.08 mg·kg-1，比较接近统计分析值，预测的空间分布模拟图高低值呈块状分布，同时出现了较多高值区域包含低值或低值区域包含高值部分，能更详细地描述土壤速效钾的空间异质性。COK3预测结果高低值呈块状分布，空间连续性较其他3种方法有所增强，插值精度明显提高，能刻画更多速效钾空间异质性的细节信息，土壤速效钾为突变而非渐变，且预测范围为36.78～124.31 mg·kg-1，与描述性数据分析结果最为接近。

表6不同方法精度对比

Table6Precision comparison of different methods

方法MethodsRMSE/(mg·kg-1)MAE/(mg·kg-1)MRE/%OK25.1119.2625.23COK124.0818.1623.82COK223.1917.0522.49COK320.2513.8919.73

3 讨论

本研究表明，研究区土壤速效钾普遍较低，空间分布受自然特征和人为活动及其协同作用的共同影响，但空间变异主要是由人为活动导致的。Pearson相关性分析结果表明pH、CEC和高程与土壤速效钾相关性达极显著水平。单因素方差分析结果表明，不同成土母质和土壤类型的土壤速效钾差异显著，不同因素的影响程度由大到小依次为CEC>高程>成土母质>土壤类型(亚类和土属)>pH。以3个辅助变量协同进行空间插值的COK3模拟精度最高，空间分布表现出“西北部高、东南部低”的总体趋势，空间异质性描述方面更加符合研究区实际情况和地学分布规律。

图2 不同方法的土壤速效钾空间分布预测结果Fig.2 Maps of soil available potassium by different methods

研究区耕层土壤速效钾均值为82.04 mg kg-1，含量处于较低水平，与孙凯等[27]、江叶枫等[28]和刘雪梅等[29]对鄱阳湖平原区的研究结果较为接近。从不同地貌类型看，均值要低于川西山区耕地[21]、秦岭山地区[30]和黄土高原区[20]。从不同土地利用类型来看，也要低于川中丘陵植烟区[25]、江南茶区[31]和冬小麦-夏玉米轮作区[32]，同时也低于福建省耕地[33]和江苏省农田[34]。这表明研究区土壤速效钾的分异可能是由于地貌类型和土地利用方式引起的。鄱阳湖平原区水热资源丰富，基岩风化、分解速率快，淋溶作用强烈，同时作为中国的商品粮基地，其土地利用强度也必然较大，土壤中钾素损失严重。

由于土壤属性获取的昂贵性与费时性，为了得到精确的土壤属性空间分布，应充分考虑影响土壤属性空间分布的因素及土壤属性的空间自相关性[28]。目前，地统计学中的克里格插值方法基于半变异函数的结构性，在考虑空间自相关的基础上对未知采样区域的土壤属性进行无偏最优估值，广泛应用于土壤属性的空间异质性研究[10-16]。OK只简单地考虑被预测土壤属性的空间自相关，仅关注于目标变量，而没有充分利用影响土壤属性空间分布的各种因素，相比之下，COK可以利用与目标变量相关性较好的辅助变量来提高目标变量的估测精度，能充分考虑影响目标变量的其他因子，将目标变量的空间自相关性和辅助变量间的交互相关性结合起来用于无偏最优估计，因而更有利于提高估测精度，这与前人研究结果一致[13-14]。本研究还发现，对于平原区，单个的辅助变量可能会在土壤属性异质性描述方面存在一定的偏见，而增加辅助变量的个数可以更加详细地刻画其异质性，在后期研究中应注意这一现象。

一般认为，耕地土壤速效养分含量受农业生产活动影响较大[7]。半变异函数分析结果表明，研究区土壤速效钾空间变异受自然特征和人为活动及其协同效应的共同作用[4]。成土母质和土壤类型是土壤速效钾空间分异的自然特征，对土壤速效钾空间变异影响程度显著。黏土矿物是土壤钾素供应的初始来源。形成各种土壤类型的基岩、黏土矿物、成土母质等类型不尽相同，母质遗传了基岩矿物特性，并构成了钾素有效性的载体[35]。尽管土壤钾素的供应潜力受制于自然特征，但自然特征在一定时期内不会发生显著变化[4]，因此研究区土壤速效钾空间变异主要是由于人为活动导致的。施肥、灌溉、翻耕，以及复种指数等均会引起土壤速效钾空间变异，然而输入与输出并非均衡，因此不同地区含量也有所差异。

本研究基于COK3统计表明，AK含量处于4级的面积为87.41%，而第二次土壤普查结果显示3级及以上的面积为55.97%，表明研究区绝大部分耕地土壤速效钾含量下降程度较大。考虑到鄱阳湖平原区是商品粮基地，为维持土壤和农业可持续发展，建议调控复种指数，增加有效钾施用。