典型白浆土区耕层土壤速效钾空间异质性研究

侯永华，王 霞，赵映慧，罗 冲，刘焕军，马士耐，王 烁

（1. 东北农业大学 公共管理与法学院，黑龙江 哈尔滨 150030；2. 中国科学院 东北地理与农业生态研究所，吉林 长春 130012）

白浆土主要分布在黑龙江省和吉林省的东北部[1]，是两省的主要耕作土壤类型之一。其中，黑龙江省白浆土面积约为331.2 hm2，主要分布于该省的三江平原地区。白浆土常与黑土、草甸土等混杂分布，三者中黑土层在相邻土壤中最浅薄，厚度介于6～20 cm，平均为12.8 cm[2]。白浆土成土母质主要是第四纪河湖沉积物，有一层瘠薄的白浆层，质地黏重，且透水性差[3]。因此，该类型土壤表层旱涝频发，且在长期耕作后易造成土壤偏酸性，对作物生长发育产生不良影响。土壤养分是土壤可持续利用与管理的重要基础，为作物生长发育提供了各种营养元素，是不可再生资源[4-5]。土壤速效钾（Available kalium，AK）作为作物生长过程中可获取的主要养分之一，直接影响着农产品的质量与品质[6]。预测白浆土区AK 含量及其空间分布，对于白浆土改良、粮食产量提升等具有重要意义。

国内外学者主要利用克里格空间插值法对土壤属性的异质性进行研究[7-8]。伴随精准农业和测土配方技术的应用与发展，学者们对土壤全氮（Total nitrogen，TN）[9-10]、有 效 磷（Available phosphorus，AP）[11-12]、有 机 质（Organic matter，OM）[13-14]等养分及微量元素[15]开展了定量研究，以分析土壤的空间变异特性。克里格插值可实现未采样区域土壤养分变化量的无偏最优估计，常被应用于土壤养分空间异质性研究。普通克里格（Ordinary kriging，OK）是GIS（地理信息系统）地统计分析中重要的分析方法之一。与OK 相比，协同克里格（Cokriging，COK）可添加与目标变量相关的其他环境辅助变量，根据与土壤养分相关性强弱实现诸多变量的叠加分析，削弱样本采集过程中如环境、人工等其他因素的影响，从而更加客观和全面地反映土壤养分的空间异质性。因此，COK 成为当前土壤养分研究常用的分析方法之一。石淑芹等[16]的研究表明，通过添加区域协同变量的COK 插值，可以提高土壤养分空间预测准确度；黄安等[17]的研究表明，控制辅助变量进入COK 插值顺序有助于提升模型预测精度；章清等[18]基于主成分分析筛选了辅助变量，并以此提升了插值模型的拟合精度；冯家东等[19]认为，筛选相关性高的辅助变量更有利于反映土壤养分的空间异质性。所研究区域不同，往往选取的辅助变量会有不同，但总体而言，COK 插值的拟合效果明显好于OK 插值法。目前尚未有白浆土区域速效钾空间异质性相关研究，因此，拟以黑龙江省桦南县曙光农场（46°13′～48°23′N，130°17′～130°39′E；典型白浆土区域）为研究区，通过分析和量化AK 空间影响因素，并将与AK 相关性较强的辅助因子作为辅助环境变量进行COK 插值，研究典型白浆土区耕层土壤AK 的空间变异规律，为典型白浆土区土壤改良提供理论依据，进一步提升典型白浆土区作物产量，保障粮食安全。

1 材料和方法

1.1 研究区概况

曙光农场位于佳木斯市桦南县，属于三江平原地区。场区东西长35 km，土地总面积177 km2，其中，耕地面积超1万hm2。地势北高南低，东高西低，呈现为由东向西倾斜的鱼脊形。土壤类型以白浆土为主，并有少部分黑土、草甸土分布，成土母质主要是第四纪河湖沉积物。该农场主要种植作物有玉米、大豆、水稻等。

1.2 土壤采样与土壤数据处理

综合曙光农场作物类型、土壤类型和地形等不同因素，在遵循全面性、均衡性和代表性的基础上布设采样点，利用S 形取样法采集0～20 cm 耕层土壤。采样时间为2019 年10 月，共采集土样186 份（图1）。在采样时记录采样点经纬度、高程（Elevation，E）等信息，带回实验室晾干，研磨，过0.149 mm 筛后测pH 值、OM 含量、AP 含量等。AK含量用乙酸铵浸提-火焰光度法测定；pH 值用电位法测定；OM 含量采用油浴加热重铬酸钾氧化-容量法测定；AP 含量采用Olsen 法测定;AN（碱解氮）含量使用碱解扩散法测定。该农场DEM（数字高程模型，图2）数据分辨率为30 m，数据来源于美国地质勘探局网站（http：//www.usgs.gov/），坡度、坡向等信息利用DEM数据通过ArcGIS 10.2软件提取。

图1 典型白浆土区耕层土壤AK实地采样点分布Fig.1 AK field sampling point distribution of topsoil in typical albic soil area

图2 曙光农场DEMFig.2 Shuguang farm DEM

土壤数据通过SPSS 26.0 进行描述性统计分析和正态分布检验（K-S 检验），利用GS 9.0+软件对AK 进行半变异函数分析和理论模型拟合，利用ArcGIS 10.2软件中地统计模块OK插值和COK插值进行空间分布特征分析。主要使用皮尔森相关系数对AK 的相关因素进行分析，在此基础上选择最佳辅助变量用于COK插值。

1.3 地统计分析

地统计学利用原始数据和半方差函数的结构性对未采样点进行无偏差估计。本研究主要采用地统计学中的半变异和克里格插值法研究AK 的空间变异性，主要公式：

式中，γ(h)为间距的半方差，在一定范围内随x（两样本间的距离）的增加而增大，当测点间距大于最大相关距离时，该值趋于稳定。Q(x)、Q(x+h)分别为随机变量Q在空间位置x、x+h上的取值。

1.4 遥感数据获取与处理

获取裸土时期（2019 年10 月3 日）Sentinel-2A影 像（https：//code.earthengine.google.com/）数 据。Sentinel-2A 影像为几何校正和大气处理后的地表反射率数据，Sentinel-2A 通过SNAP（哨兵应用平台）重采样为10 m 空间分辨率。在此基础上计算该时期归一化差异植被指数（Normalized difference vegetation index，NDVI）、归 一 化 差 分 水 体 指 数（Normalized difference water index，NDWI）、比值植被指数（Ratio vegetation index，RVI）等。将这些遥感指数与AK 含量进行相关性分析后择优选取作为辅助变量进行COK 插值，将所得结果与OK 插值所得结果及其他COK插值所得结果进行对比分析。

1.5 拟合精度评价

通过交叉检验对OK和COK插值所得结果进行验证。主要选取均方根误差（RMSE）、平均相对误差（MRE）和平均绝对误差（MAE）作为检验精度的评价指标，计算公式：

式中，X*i为样点的预测值；Xi为样点实测值；n为样点数。RMSE、MRE和MAE值越小，表示插值拟合效果越好，预测精度越高。

2 结果与分析

2.1 典型白浆土区耕层土壤AK描述性统计分析

研究区AK 含量介于92.00～364.00 mg/kg，平均为171.78 mg/kg（表1）。为使该数据符合正态分布，将其进行对数变换。变换后AK 含量介于1.96～2.56 mg/kg，平均为2.22 mg/kg。AK含量变异系数由27.52%降低为对数变换后的5.32%。根据全国第2次土壤普查标准[AK含量＞100.00 mg/kg处于三级高水平]，本研究区AK含量大部分处于高水平范围内。AK 含量最大值与最小值差异较大，表明AK 含量在研究区内的值域分布较广。

表1 典型白浆土区耕层土壤AK总体描述性统计Tab.1 Overall descriptive statistics of AK in typical albic soil area

2.2 典型白浆土区耕层土壤AK主要相关因素分析

2.2.1 典型白浆土区耕层土壤AK 辅助变量相关性分析 皮尔森相关分析结果表明，AK 含量与AP 含量、OM 含量、AN 含量的r值分别为0.363、0.173、0.194，这3 种土壤养分含量与AK 含量呈极显著正相关，均为区域协同变量，在受外界环境影响方面具有相似性。AK 含量与E、地形起伏度（Topographic relief，TR） 、 地 表 粗 糙 度（Surface roughness，SR）这3 个地形要素的r值分别为-0.371、-0.221、-0.278，均呈极显著负相关，证明地形为AK 含量的重要影响因子。另外，在研究区内AK 含量随海拔升高和地形起伏变大而降低。研究区位于三江平原，属于温带大陆性气候，降水多集中于夏季，对土壤侵蚀效应明显，受雨水及田块微地形影响，AK 流失明显。AK 含量与pH 值、坡度、坡向等无显著相关关系，因此，在后续插值中剔除这些因素。

结合影像与实地调查发现，研究区部分耕层土壤长期深翻后白浆层被打破，在遥感影像上出现真彩色状态下发白发亮的现象。通过NDVI、NDWI、RVI 等遥感指数与裸土时期影像特征对比分析发现，研究区内裸土期NDWI、NDVI 与表层白浆化土壤在影像上的特征具有高度一致性。另将遥感指数与AK 含量进行相关性分析后发现，裸土期NDWI、NDVI 与AK 含量的r值分别为0.164、-0.194，存在显著相关关系，择优选取NDVI进行COK插值。

2.2.2 典型白浆土区耕层土壤AK 辅助变量回归分析 通过SPSS 26.0 进行线性回归分析，选取了AP、E、TR、NDVI等多种辅助变量进行回归分析，部分结果如表2 所示。结果显示，上述辅助变量和研究区AK 含量存在显著或极显著相关关系。其中，AN 和AP的标准化系数分别为0.194、0.363，E、TR和SR标准化系数分别为-0.371、-0.221、-0.278，均与AK 含量呈极显著相关；NDWI、NDVI 的标准化系数分别为0.164、-0.194，均与AK 含量呈显著相关。该结果与皮尔森相关系数分析结果一致。

表2 典型白浆土区耕层土壤AK不同辅助变量回归分析Tab.2 Regression analysis of different auxiliary variables of topsoil AK in typical albic soil area

2.2.3 典型白浆土区耕层土壤AK 半变异函数分析 在GS 9.0+中对AK 进行半变异函数拟合分析发现，最优拟合模型为球状模型。半变异函数块金值为0.04，变程为8.26 km。按照半方差分析法和空间相关性程度的分级标准，研究区较大范围内AK存在强烈空间变异。R2为0.886，表明模型具有良好的拟合优度。

2.3 典型白浆土区耕层土壤AK模拟空间分析

为分析研究区内AK的空间异质性，将2种土壤养分因子（AP、AN）和2 种地形因子（E、TR）以及遥感指数（NDVI）这5 种因子作为辅助变量分别与AK进行COK 插值。为了比较OK 插值与COK 插值的拟合精度，将AK 进行OK 和COK 插值，得到6 幅AK空间分布图，以可视化图表征AK 空间异质性，如图3 所示。以AP、AN 和TR 为辅助变量的COK1、COK2和COK4插值结果均为斑块状分布，以E 和NDVI 为辅助变量的COK3和COK5插值结果均为条带状分布，且多种插值所得结果具有总体趋势一致性，AK含量表现为南高北低。

图3 典型白浆土区耕层土壤AK空间分布Fig.3 Spatial distribution of AK in typical albic soil area

交叉验证结果见表3。由多种不同验证数据来看，增加了辅助变量的克里格（协同克里格，COK）插值法优于单一目标变量的克里格（普通克里格，OK）插值。COK1—COK5的RMSE分别较OK 下降了2.096、0.992、0.300、1.459、0.173 mg/kg；MRE分别较OK 降 低 了1.51、0.25、0.22、0.65、0.08 个 百 分 点，MAE分别较OK 降低了0.50、0.16、0.07、0.35、0.02 mg/kg。从结果来看，以土壤养分因子（AP、AN）、地形因子（E、TR）和遥感指数（NDVI）作为辅助变量的COK 插值法精度更高，即COK＞OK，且以AP 作为辅助变量的效果最好。当加入2种及以上因子作为辅助变量时COK 精度会有所降低，如将AP+TR 作为辅助变量时，COK 交叉验证RMSE为38.188 mg/kg，该结果精度低于基于AP 的COK（COK1，RMSE为36.225）mg/kg。以上表明，COK 插值并非辅助变量越多越好，应该充分考虑与目标信息的相关程度，择优选取协同变量，这样才能获得更高的拟合精度。

表3 典型白浆土区耕层土壤AK不同辅助变量插值精度系数Tab.3 Interpolation precision coefficient of different auxiliary variables of topsoil AK in typical albic soil area

3 结论与讨论

一般认为，土壤速效养分与农业生产活动密切相关[20]。曙光农场内AK 含量的空间变异系数为27.52%，与半变异函数的分析结果相一致，均属于中等变异程度，并且与AK 在空间上的分布存在较强相关性。因此，认为该地的AK 同时受自然环境因素和人为因素影响。2 种插值方式描述的AK 含量空间分布规律基本一致，均呈南高北低趋势。插值结果除COK1、COK2和COK4插值均为斑块状分布外，其余均为条带状分布。从研究区整个农场范围来看，AK 的空间分布存在一定规律性，考虑到该地地形及耕地类型，研究区总体北高南低，且南部均为水田，各种土壤养分随水分流失后聚集于南部区域，提高了南部区域多种养分含量。

地统计学中的克里格插值利用原始数据的空间依赖性，能够实现未知空间的无偏最优估计，进而达到土壤养分空间变异规律的最优拟合，此方法多应用于土壤属性的空间异质性研究[21]。与OK 插值法相比，COK 插值法将待测值无偏最优估计方法从单一土壤属性提高到2 个及以上协同变量，COK插值无论从模型拟合程度还是交叉验证精度来说都优于OK插值[22]。

本研究结果显示，无论辅助变量增加土壤养分因子（AP、AN）的COK1、COK2，还是增加地形因子（E、TR）的COK3、COK4，亦 或 是 增 加 遥 感 指 数（NDVI）的COK5，其RMSE、MRE和MAE值都低于OK。综合比较后认为，拟合精度表现为COK＞OK，证明了在同样条件下添加具有较强相关性协同变量的COK 插值比OK 插值精度更高。这与赵业婷等[23]的研究结果相一致，COK 对土壤属性变异特征的拟合优度更高。本研究结果表明，COK 在农场尺度典型白浆土区土壤养分属性的模拟和预测中效果更优。

在众多环境变量中选取相关性强且更能反映目标环境信息的辅助因子可获取更加精确的土壤养分空间分布信息。本研究比较了添加土壤养分因子、地形因子、遥感指数与单一变量的OK 插值交叉精度验证结果，最终模型拟合效果为COK1＞COK4＞COK2＞COK3＞COK5＞OK。本研究选取E 作为辅助因子时插值效果不佳，可能是研究区位于三江平原地区，属于典型白浆土区，土壤表层旱涝频发，受降雨影响土壤被侵蚀后带走大量养分，导致土壤本身原始养分含量降低，最终导致预测结果具有不稳定性。虽然在本研究中坡度与AK 相关性较低，但以坡度为辅助变量的插值精度仍然较高。原因可能是数字高程模型分辨率不足以说明坡度与养分之间的关系。江叶枫等[24]的研究表明，辅助变量数量越多模拟精度越高。但本研究中以AP+TR 作为辅助变量的插值精度低于以AP 为辅助变量时的插值精度，证明了虽然辅助变量有助于模型精度的提高，但并非越多效果越好。与目标变量相关性较小的辅助变量可能存在未知因素带来的影响，甚至可能会对拟合精度产生负影响，这与冯家东等[19]的研究结果一致。另外，选取恰当的遥感指数加入克里格插值也有助于提升插值精度。例如，加入NDVI 的COK5插值精度优于单一的OK 插值精度。随着遥感技术的发展，后续可进一步分析其他遥感指数的适用性，以刻画更加完善的环境变量，获取更加详尽的耕地土壤养分空间分布信息。

本研究主要利用地统计学方法并结合克里格空间插值技术，对典型白浆土区AK 的空间变异特征进行分析。结果表明，曙光农场耕层土壤AK 含量为92.00～364.00 mg/kg，平均为171.78 mg/kg，耕层土壤AK 含量整体处于较高水平。空间插值结果表明，研究区内土壤钾素含量南高北低，部分斑块状区域含量较高，为实现土壤改良，应分区域调整钾肥施用量。本研究中，模型拟合效果依次为COK1＞COK4＞COK2＞COK3＞COK5＞OK，以土壤养分因子AP为辅助变量的COK1较适合于典型白浆土区耕层AK的空间分布研究，而单一变量的OK 插值效果较差。综上，选择最优辅助变量后的COK 插值法获取典型白浆土区耕层土壤养分空间分布信息切实可行。