基于地统计学和GIS技术的北京市大兴区礼贤镇土壤养分空间变异性研究

安永龙，杜子图，黄 勇

(1.北京市地质勘察技术院，北京 102218；2.中国地质调查局 发展研究中心，北京 100037)

0 引 言

土壤是由母质经过长时间的成土作用而形成的三维自然体，成土过程中由于受到物理、化学、生物等外界因素的不断影响，加之人类活动对其造成的不确定性干扰，使得土壤在时间和空间上具有高度的异质性[1-2]。单凭传统的统计学方法无法解决空间数据与属性数据关联的问题，目前许多国内外学者将地统计学与GIS技术相结合研究土壤中微量元素的空间变异性取得了丰富的成果。如，武婕等[3]研究了南四湖区农田土壤有机质和微量元素的空间分布特征及其影响因素；朱菊兰等[4]分析了浑河太子河流域表层土壤养分空间分布以及地形因子、土地利用对土壤养分空间变异的影响；Kay等[5]通过不同土地利用方式及地形因素，对稻田土壤养分的空间分布进行了预测；Roger等[6]以部分环境因子作为辅助数据，对土壤养分指数进行回归克里格插值，分析出线性回归和克里格相结合可以提高土壤养分预测精度；Cordova等[7]分析了英国洛桑实验区土壤中氮元素的空间变异性。

虽然目前已经有不少关于北京地区土壤养分指标空间变异性的研究，但是这些研究的采样密度较稀，赵汝东[8]在研究北京市平原地区耕地的土壤养分空间变异性时，共采集299件土壤样品；董士伟等[9]在揭示北京市大兴区土壤养分空间分布特征时，共采集1 240件土壤样品；由于土壤的空间变异性较强，采样精度的高低和不同插值方法的选择会直接影响研究成果的准确性。因此本文采用地统计学和GIS方法，对大兴区礼贤镇进行高精度高密度采样，共采集土壤样品1 789件，对其中全氮、全磷、全钾、全硼、碱解氮等10项养分指标的空间变异性、丰缺格局及其变化特征进行了详细的论证，旨在为当地土壤中养分管理、平衡施肥和提高土地利用率提供重要依据。如此对乡镇级进行高密度采样来研究其空间变异性在全国可谓首例，同一地形单元内的克里格插值精度随采样密度增加而提高[10]，对于提高乡镇级养分指标空间克里格插值精度具有先行性作用。

1 材料与方法

1.1 研究区概况

礼贤镇地处北京市大兴区最南部，北邻庞各庄镇、魏善庄镇、安定镇，西靠榆垡镇，东接河北省廊坊市，位于永定河洪积冲积平原，地势平坦。研究区地处中纬度，属暖带半湿润大陆季风气候，气候特点是四季分明，春旱多风，夏季炎热多雨，冬季寒冷干燥。年平均气温11.6 ℃，全年无霜期平均210天左右，年平均降水量556.4 mm。研究区土壤类型主要有砂质潮土、壤质潮土、硫酸盐盐潮土、草甸沼泽土、风沙土等。沉积物主要由永定河冲洪积而成，岩性主要为卵石、砂、黏性土及粉土。研究区总面积为92.06 km2，土地利用类型多样，其中基本农田保护区约占总面积的73%，建设用地区约占14%，一般农地区约占10%，其他土地利用地区约占3%，全区农业用地类型占比较大。

1.2 样品采集与分析方法

本次取样采用GPS定位，充分考虑土壤类型及土地利用方式，按照64个样/km2的采样密度进行采集，共采集1 789件表层土壤样品，如图1所示。土壤采样深度为0～20 cm，要求定位点四周辐射20 m，采集3～5个分样点，等组分混合为一个组合样，并基本保证了每个采样单元内采集的土样类型的一致性。采样工作中选择在农用大田、菜地、林带地、山坡下侧土层较厚处等地采样，尽量避开水土流失严重或表土被破坏、明显点状污染、新近搬运的堆积土、垃圾土、田埂等地段，以保证样品具有代表性。

样品由北京一零一生态地质检测有限公司参照国家环境质量相关标准进行测定。元素全量测定方法：全氮使用氧化燃烧-气相色谱法测定，全钾、全磷使用电感耦合等离子体光谱法测定，全硼使用发射光谱法(ES)测定，有机质使用管式炉燃烧红外吸收法测定，pH值使用去二氧化碳水将样品水溶，剧烈搅拌后静止半小时，pH计直接测定。

图1 研究区土壤样点分布Fig.1 Distribution of soil sampling sites in the study area

元素有效态测定方法：使用ICP-OES法测定有效硼，称风干土样(通过2 mm筛)20.0 g于250 mL锥形瓶中，加40 mL蒸馏水，连接冷凝管，文火煮沸5 min，加4滴0.5 mol/L CaCl2溶液，离心分离得清液；使用ICP-OES法测定有效磷，称风干土样(通过2 mm筛)2.5 g于浸提瓶中，加50 mL 0.5 mol/L碳酸氢钠浸提剂，在25 ℃±1 ℃温度中振荡30 min后过滤；使用ICP-OES法测定速效钾，称风干土样(通过2 mm筛)5.0 g于浸提瓶中，加50 mL乙酸铵溶液，加塞振荡30 min后过滤；称风干土样(通过2 mm筛)2.0 g平铺于螺纹丝扣密封类型的塑料扩散皿外室，吸取3 mL 20 g/L硼酸-指示剂溶液(pH=4.5)于扩散皿内室，加10.0 mL 1.8 mol/L氢氧化钠溶液于扩散皿外室，轻摇晃动后迅速旋紧盖子，于40 ℃保温箱内恒温24 h，0.01 mol/L HCl标准溶液滴定内室硼酸中的氨，计算碱解氮含量。

为了保证数据的准确性，选取GSS1、GSS2、GSS3和GSS4共4个国家一级标准物质以及GBW07413a、GBW07414a和GBW07415a共3个有效态标准物质，用选定分析方法对每一个标准物质进行8次分析检验，并参照《多目标区域地球化学调查规范(1∶250 000)》(DZ/T 0258—2014)和《生态地球化学评价样品分析技术要求(试行)》(DD2005—03)中关于合格率标准，分别计算被测项目每个样品平均值与标准值之间的相对误差(RE%)和相对标准偏差(RSD%)，测试结果中相对误差的合格率全部达到100%，相对标准偏差的合格率全部达到90%以上，标样回收率均在93%～107%之间，全部满足分析测试质量的相关要求。

1.3 数据处理

正态分布性检验、相关系数计算、因子分析计算均由SPSS19.0[11]完成，因子分析采用最大四次方值法，最大收敛性迭代次数为25次。只有当数据全部服从正态分布时，才能进行因子分析，克里格插值方法才有效，否则可能存在比例效应[12]。判断数据是否服从正态分布方法较多，如探索性检验中柱状图或Q-Q图[13]、单样本柯尔莫哥洛夫—斯米诺夫(K-S)正态分布检验[14-15]等，本文采用偏度峰度检验法，结果显示有效磷、速效钾、碱解氮的偏度值和峰度值远大于1，而其余各指标值基本在1附近，符合正态分布的要求，随后对有效磷、速效钾、碱解氮值进行对数变换，满足地统计学分析的假设条件，进行半变异函数的计算。

半变异函数的计算和理论模型的拟合均采用地统计学软件GS+9.0[16]，变异函数作为研究空间变异的重要函数，是分析土壤养分指标空间变异随机性或结构性的基础性工具。该函数如下：

(1)

式中：γ(h)表示变异函数；N表示以h为步长的观测数据成对数目；[Z(xi)-Z(xi+h)]表示间隔为h的2个观测点的实测值，也称步长。

在本次研究中为了更加准确地反映土壤养分指标的空间变异性，根据半变异函数的决定系数R2和残差RSS进行拟合，构建各自的理论模型[17]。本次研究所建立的半方差函数模型有指数(Exponential)和线性(Linear)两种[18]。

指数模型：

(2)

线性模型：

γ(h)=C0+Ch/αh>0

(3)

式中：C0，块金方差，即间距为0时的半方差，由随机因素引起的变异；C，结构方差，由系统因素引起的变异；C0+C为基台值，即半方差函数随间距递增到一定程度后出现的平稳值，表示系统内总的变异；h，步长；α，变程，即半方差达到基台值的样本间距，表示随机变量在空间上的自相关尺度。

克里格插值由ArcGIS 10.2中Geostatistics analysis模块的普通克里格空间插值法完成[19]。通过变量在点x处影响范围内的n个有效观测值Z(xi)的线性组合得到在该点x处的估计值Zx，如下：

(4)

式中：Zx为采样点之间距离的估计值；Z(xi)为采样点之间距离的观测值；λi是赋予观测值Z(xi)的权重，表示各个观测值Z(xi)对Zx的贡献。

2 结果分析

2.1 土壤养分的统计分析

研究区表层土壤养分指标统计分析是以统计学为基础，通过计算各样点的最小值、最大值、中位数、平均值、标准差等统计特征来描述土壤营养元素的含量和分布特征(表1)。表层土壤中全量养分、pH值、有机碳样品的平均值与中位数相差不大，加之偏度值分布在1附近，表明表层土壤中这些指标的中心趋向分布呈标准正态分布，而有效养分的平均值与中位数相差较大，加之偏度值远超过1，表明有效养分的中心趋向分布可能被异常值影响而使其呈非标准正态分布。由于非正态分布的数据会增加变异函数的波动性，导致误差增大，因此对有效磷、碱解氮和速效钾进行对数转换，转换后的数据服从对数正态分布，因而在计算有效磷、碱解氮和速效钾的变异函数值时用转换后的数据。研究区表层土壤整体为碱性环境，pH值为7.55～9.36。

大量研究表明变异系数(Cv)可以反映土壤元素的离散程度，进而表示土壤元素分布的均匀程度。变异系数越大，元素的分布越不均匀，离散程度越大。本次研究参照测试数据范围将变异程度分为4种类型：Cv<30%为均匀分布，30%≤Cv<60%为弱分异(中等起伏)，60%≤Cv<100%为较强分异(较大起伏)，Cv≥100%为强分异(很大起伏)。参照这一分类标准显示有效磷、速效钾变异系数分别为105.99%、100.43%，表明表层土壤中有效磷和速效钾属于强分异，以有效磷的分异程度最强。相关统计发现，土壤养分的来源、施用率及其在土壤当中的流动性是影响土壤养分空间变异程度的主要因素。说明局部地区在施用磷肥和钾肥之后，短期内植物对土壤中的有效磷和速效钾吸收不充分，导致土壤中有一部分残留，但最终这些成分彼此会达到动态平衡状态[20]；碱解氮的变异系数为80.29%，属于较强分异，由此可见土壤中养分有效态形式的分异性相对较强；有效硼、有机质、全氮的变异系数为32.46%～38.69%，说明三者的分异程度较弱；全磷、全硼、全钾、pH值的变异系数为2.64%～21.28%，四者为均匀分布，说明这几种养分指标在土壤中较为稳定，其中以pH值的分异程度最低，仅为2.64%，此项研究与国内外部分学者的研究结果相似[21-22]，说明表层土壤酸碱性是一个较为稳定的指标，受施肥、灌溉等随机性因素影响较小。

2.2 土壤养分含量之间的相关性分析

表层土壤10个养分指标之间的相关分析(Spearman相关系数)表明，大部分的指标之间存在显著或极显著的相关关系(表2)。土壤中有机质与全磷、全氮、有效硼之间呈极显著的正相关，与有效磷、碱解氮之间呈显著的正相关，与速效钾、全钾、全硼之间相关性较弱。速效钾与全钾呈显著的正相关关系，而两者与全硼之间的相关性较弱。土壤pH值与其余指标间存在负相关关系，这与部分学者研究结果相似[23-25]。

2.3 土壤养分主成分分析

主成分分析法就是要在数据信息丢失最少的前提下，对高维度变量进行降维处理[26-28]。对研究区表层土壤中9种养分指标进行主成分分析，前3个主成分特征参见表3，前3个成分的累计方差贡献率为85.45%，集中了原始数据的大部分信息。第一个主成分的方差贡献率为61.39%，占3种主成分中最大比例，因此可以作为衡量该研究区养分状况的一个综合性指标。为了进一步对主成分所包含信息进行研究，对9种养分指标的3种主成分上的因子载荷进行统计(表4)。结果显示，全磷、有效磷、全氮、碱解氮、有效硼、有机质在第一主成分中普遍具有高载荷，而在第二主成分和第三主成分中均为低载荷；全钾和速效钾在第二主成分中普遍具有高载荷，而在第一主成分和第三主成分中均为低载荷；全硼在第三主成分中普遍具有高载荷，而在第二主成分和第三主成分中均为低载荷。

表1 土壤养分的描述性统计

注：全氮、全钾和有机质的单位为%，有效磷、碱解氮、速效钾、全硼和有效硼的单位为10-6。

表2 土壤10个养分指标的相关系数

注：*表示在 0.05 水平(双侧)上显著相关，**表示在0.01水平(双侧)上显著相关，n=1 789。

综合相关性分析和主成分结果进一步分析，可知全磷、有效磷、全氮、碱解氮、有效硼、有机质来源相同，均以结构性因素为主，全钾和速效钾以随机性因素为主，全硼则以两者混合为主。

表3土壤养分的主成分分析结果

Table3Principalcomponentanalysisresultofsoilnutrients

主成分特征值累计特征值方差贡献率/%累计方差贡献率/%主成分18.739.7361.3961.39主成分21.1510.8812.8274.21主成分31.0111.8911.2485.45

表4主成分分析旋转后的成分载荷矩阵

Table4Rotatedcomponentmatrixofprincipalcomponentanalysis

养分指标变量在各主成分上的因子载荷第一主成分第二主成分第三主成分全磷 0.819-0.203-0.048有效磷0.732-0.003-0.396全钾 -0.0120.7350.179速效钾0.4210.614-0.358全氮 0.8810.0110.185碱解氮0.6240.322-0.121全硼 0.2000.0710.759有效硼0.6920.2240.226有机质0.776-0.1430.295

注：旋转方法采用具有Kaiser标准化的四分旋转法，旋转在5次迭代后收敛。

2.4 土壤养分元素的趋势分析

图2 土壤微量养分元素趋势分析图Fig.2 Trend analysis of soil microelements

自然界中由于土壤易受到多种因素的共同影响，因此其空间分布具有较为明显的趋势效应特征，本文运用ArcGIS 10.2中地统计模块的趋势分析功能，探讨了各项指标的空间分布的趋势效应特征。通过异向性轴轴向自动搜索功能，最终确定各向异性的情况如下，选取无(没有趋势效应)、一阶(区域化变量沿一定方向呈直线变化)、二阶(区域化变量沿一定方向呈多项式变化)[29-30]。图2为研究区中选取3类典型指标的空间分布趋势效应分析图，图中x轴表示正东方向，y轴表示正北方向，z轴表示各样点测定值的大小；左后投影面上浅绿色曲线表示东西向的全局性趋势效应变化，右后投影面上浅蓝色曲线表示南北向的全局性趋势效应变化。图2显示，研究区土壤中全氮表现为在东西方向上先减小后增大的趋势，属于二阶类型；有效磷表现为在东西方向和南北方向上呈直线分布，属于一阶类型，该类型还包括全磷、全钾、速效钾、全硼和有效硼；有机质表现为无阶效应类型，该类型还包括碱解氮和pH值。

结合克里格插值法所造成的插值误差进行比较，综合考虑各个参数选择的合理性，评价标准为平均误差(ME)的绝对值最接近0，如果平均标准误差(ASE)大于均方根误差(RMSE)，说明高估了预测值，反之则说明低估了预测值，因此两者尽量相等。标准化平均误差(MSE)的绝对值最接近0，标准化均方根误差(RMSSE)值最接近1，如果标准化均方根误差大于1，说明低估了预测值，反之为高估预测值[31]。由表5可知，全磷、有效磷、全钾、速效钾、全硼、有效硼的一阶趋势预测与无趋势预测和二阶趋势预测相比，ME和MSE的绝对值更加接近0，ASE和RMSE最接近，RMSSE的值最接近1，因此克里格插值时应选择一阶更加精准，同样碱解氮、有机质和pH值应该选择无趋势，而全氮则应该选择二阶。

2.5 土壤养分空间异质性特征

半变异函数(semi-variogram)是地统计学中应用最广泛的空间格局描述的基本工具[32-33]，不但可以利用半变异函数对参数的空间分布进行结构性和变异性分析，而且利用这一结果结合克里格法可以对未知点进行预测和模拟。因此应用地统计学方法对研究区土壤养分指标的空间异质性特征进行半变异分析。其中有3个重要的参数，块金值(C0)代表在最小取样距离内参数的变异，属于由测量误差等随机因素引起的随机变异；变程(A0)也称空间最大相关距离，反映了变量空间自相关范围的大小；块金系数是块金值与基台值的比值，当块金系数小于25%时说明变量具有较强的空间相关性，其空间变异主要由结构性因素如种植制度、施肥措施等控制[34]。当块金系数介于25%～75%，说明变量具有中等强度空间相关性，其空间变异由随机性因素与结构性因素共同控制。若块金系数大于75%，则说明研究变量的空间相关性较弱，其空间变异主要受随机性因子控制[35-37]。随机性因子能够减弱土壤养分的空间相关性，而结构性因素能够增加土壤养分的空间相关性[38]。半变异函数模型的选择是土壤养分空间变异性研究的关键，R2值越接近1.0，RSS值越小，表明拟合模型的效果越好。

表5 不同趋势阶数插值误差比较

由表6可知，10项养分指标的R2值为0.718～0.956之间，说明上述养分指标的实验半变异与相应的理论半变异模型非常吻合，拟合效果好。这10项指标中有5项指标的变程超过了8 000 m，最小的变程也在200 m以上，因此对其进行克里格插值能得到可靠的结果[39]。

全磷的半方差函数属于指数模型，变程为210 m，块金系数为12.4%，有效磷的半方差函数属于指数模型，变程为230 m，块金系数为13.2%，全钾的半方差函数属于指数模型，变程为650 m，块金系数为11.1%，全硼的半方差函数属于指数模型，变程为270 m，块金系数为11.8%，pH值的半方差函数属于指数模型，变程为210 m，块金系数为12.4%，这5项养分指标都具强烈的空间相关性，说明研究区表层土壤中这5项指标主要受到气候、地形、成土母质、土壤类型等结构性因素的影响；全氮的半方差函数属于线性模型，变程为8 409.81 m，块金系数为74.2%，有效硼的半方差函数属于线性模型，变程为8 409.81 m，块金系数为71.9%，有机质的半方差函数属于线性模型，变程为8 409.81 m，块金系数为71.9%，这3项养分指标都具中等强度的空间相关性，说明研究区表层土壤中这3项指标主要受到结构性因素和随机性因素的共同影响；速效钾的半方差函数属于线性模型，变程为8 409.81 m，块金系数为86.4%，碱解氮的半方差函数属于线性模型，变程为8 409.81 m，块金系数为84.5%，这2项养分指标都具较弱的空间相关性，说明研究区表层土壤中这2项指标主要受到施肥、灌溉、耕作措施等随机性因素的影响。

2.6 普通克里格空间插值法

GIS在数据统计模型的建立和制图表达方面具有很强大的功能，它的插值图可视化功能可以直观地反映土壤养分指标的空间变异特征。空间插值的方法较多，如反距离权重法[40]，泛克里格法[41]，规则样条函数法[42]等，本次研究选取地统计学中最常用的普通克里格空间插值法。普通克里格空间插值法也称空间局部估计或空间局部插值法，它是建立在半变异函数理论及结构的分析基础上，在有限区域内对区域化变量未知点进行线性无偏最优估计的一种方法[43]。

如图3所示，参照《北京市土壤养分分等定级标准》以及北京市平原区土壤元素背景值[44]，可知全氮的含量为0.05×10-2～0.27×10-2、全钾的含量为1.37×10-2～2.36×10-2，整体较为缺乏；碱解氮的含量为47.9×10-6～241.89×10-6、全磷的含量为652.1×10-6～1 920×10-6、有机质含量为0.55×10-2～2.66×10-2、有效硼含量为0.38×10-6～2.35×10-6，整体处于中等水平；有效磷的含量为6.42×10-6～92.37×10-6、速效钾的含量为79.42×10-6～2 691.17×10-6和全硼含量为45.55×10-6～74.11×10-6，总体较为丰富。研究区土壤pH值范围在8.19～8.96之间。从空间分布来看，10项养分指标具有相对较明显的区域性特征，低值区分布具有以下特点，全氮、碱解氮、全磷、有效硼、有机质主要位于研究区中南部的2条北西向展布的低值带，分别为王化庄—赵家庄—小马坊—紫各庄一带、大马房—辛家安—祁各庄—杨各庄一带，其余低值带零星分布于研究区边缘，有效磷则位于其中一条低值带，即王化庄—赵家庄—小马坊—紫各庄一带。而这6项养分指标的高值区都分布在研究区的西部和东北部边缘地带。全钾和速效钾低值区的分布与之不同，全钾低值区分布在研究区大辛庄—良种场—杨各庄一带、东段家务—于场—内官庄一带，速效钾零星分布在研究区北部和南部边缘。全硼则在李各庄和紫各庄等地含量相对较低。

表6 土壤养分含量变异函数理论模型及其相关参数

图3 土壤养分指标的空间分布Fig.3 Spatial distribution of soil nutrients

研究区中土壤pH值全部为碱性，碱性较高地区分别为王化庄—赵家庄—小马坊—西段家务一带、大马房—辛家安—祁各庄—杨各庄一带，因此可以看出，土壤碱性越高，土壤中全氮、碱解氮、全磷、有效硼、有机质、有效磷含量越低。而全钾、速效钾、全硼分布则与pH值关系较弱，同时三者含量与有机质之间无相关关系，因此得出，全钾、速效钾、全硼分布与土壤结构性因素关系较弱，主要受人为活动影响较大。这与前文中相关分析、主成分分析所得的结果一致。

因此，研究区应按照养分指标空间分布格局以及所种植农作物对养分吸收特点进行科学有效的施肥，按照不同养分指标进行“局部圈定，分类施肥”，尤其要加大氮肥和钾肥的投入量，调整具体肥料的用量比，实现合理施肥，这样既可以避免肥料施用过量后造成的污染和浪费，又可以保障土壤中养分比例得当，增加肥料的利用效率，有效地提高农作物的产量与质量。

3 结 论

农田区分布面积在北京市大兴区礼贤镇占有绝对优势，因此土壤养分指标情况成为当地农业规划管理的关键性因素。通过对礼贤镇表层土壤10项养分指标空间变异及分布格局特征进行研究，得出以下结论：

(1)经典统计学分析中有效磷、速效钾属于强分异，变异系数分别为105.99%、100.43%，碱解氮属于较强分异，变异系数为80.29%；其余指标分异性相对较弱。综合相关性分析和主成分结果进一步分析，可知全磷、有效磷、全氮、碱解氮、有效硼、有机质来源相同，均以结构性因素为主，全钾和速效钾以随机性因素为主，全硼则以两者混合为主。

(2)通过异向性轴轴向自动搜索功能，结合不同趋势阶数插值误差综合对比，确定了全磷、有效磷、全钾、速效钾、全硼、有效硼应该选择一阶，碱解氮、有机质和pH应该选择无趋势，而全氮则应该选择二阶。

(3)地统计学分析结果表明，全磷、有效磷、全钾、全硼、pH值可以用指数模型较好地拟合，速效钾、全氮、碱解氮、有效硼、有机质可以用线性模型较好地拟合。其中全磷、有效磷、全钾、全硼、pH值的块金效应都小于25%，因此具强烈的空间相关性，受到结构性因素的影响；全氮、有效硼、有机质的块金效应为25%～75%，因此具中等强度的空间相关性，主要受到结构性因素和随机性因素的共同影响；速效钾、碱解氮的块金效应大于75%，具较弱的空间相关性，受到随机性因素的影响。

(4)普通克里格插值图直观地反映了礼贤镇土壤养分的分布特征。研究区发现2条全氮、碱解氮、全磷、有效硼、有机质北西向展布的低值区，主要分布于中南部，即王化庄—赵家庄—小马坊—紫各庄一带、大马房—辛家安—祁各庄—杨各庄一带；有效磷的低值区归属于其中1条。全钾和速效钾的含量分布区域相似。同时发现研究区土壤碱性越高，土壤中全氮、碱解氮、全磷、有效硼、有机质、有效磷含量越低，这都与经典统计学结论相似。