基于地统计学和GIS的通州区于家务乡土壤肥力综合评价

杨全合，安永龙

(北京市地质勘察技术院，北京 102218)

【研究意义】土壤肥力是衡量土壤养分针对特定植物的供应的能力，也是土壤区别于成土母质和其他自然体的最本质的特征。土壤肥力综合评价是依据对土壤资源调查的数据，针对土壤资源的利用价值和可利用潜力进行的评估与预算[1]。以往由于受到技术和方法条件的限制，对土壤肥力的评价一般采用定性的方法进行描述，评价的准确度偏低。随着地统计学和GIS技术在土壤领域研究中的不断发展和应用，土壤评价的定量研究技术脱颖而出，大幅度提升了土壤评价工作的精准性和确定性。【前人研究进展】目前国内已有不少的学者成功应用定量化技术对土壤肥力进行综合评价，郑城等[2]利用统计学方法分析了小流域内7种地类的土壤有机质、全氮、速效磷和速效钾的特征，并对土壤养分进行了综合评价；叶回春等[3]采用传统统计学、隶属度函数和地统计方法，对北京市延庆盆地表层土壤肥力进行综合评价，并研究其空间变异特征；王妙星等[4]运用统计学方法分析第十三师耕层土壤养分现状和演化趋势，结合层次分析法和特尔斐法开展了地力评价，国外也有许多学者对此开展过研究[5-8]。【本研究切入点】通州区位于北京市东南部，京杭大运河北端，自2012年北京市委明确提出要将通州区建立成为北京唯一的行政副中心[9]，其影响力迅速上升，成为各界人士关注的焦点。通州是北京市重要的农副产品生产、加工基地，适时调整农业结构，及时提高作物肥力供给效率显得尤为重要，由此可见及时掌握通州区土壤肥力做好评价工作意义重大。【拟解决的关键问题】本次研究对通州区于家务乡进行高精度高密度采样，共采集土壤样品1399件，综合运用地统计学方法和GIS技术对其中全氮、全磷、全钾、碱解氮等8项肥力指标进行了综合评价，为于家务乡土壤规划管理和改良提供基础资料。

1 材料与方法

1.1 研究区概况

研究区位于通州区南部的于家务回族乡，其东邻永乐店镇，西南部与大兴区采育镇接壤，是通州区唯一的少数民族乡，总面积约为6570 hm2，其中耕地面积约为4640.8 hm2，国家级民族农产品加工基地、航天育种基地都位于此。通州属温带大陆性半湿润气候区，夏季炎热多雨，冬季寒冷干燥，年日照2730 h，年降水量620 mm左右。地处永定河、潮白河冲积平原地带，从属华北大平原[10]。研究区地势呈现由东北至西南向上的波状起伏之势。通州区土壤质地主要以轻壤质和砂壤质为主，研究区内主要以轻壤质为主；通州区土壤主要包括褐土、潮土、风沙土、沼泽土4种类型，本次研究区中以壤质潮土、粘质潮土、硫酸盐盐潮土为主[11]。

1.2 土壤样品采集与测定方法

布置采样点时充分考虑到了土地利用类型、土壤类型、流域分布等因素，为了尽量降低空间变异性所带来的误差，本次研究加大了采样密度、尽量缩小采样间距，最终按照32～64/km2的采样密度进行布设，共布设1399件表层土壤样品(图1)。野外用校正准确的GPS寻找采样点位，一件完整的土壤表层土壤样品是由3～5个子样点混合均匀后组合成，并去除杂草、草根、砾石、砖块、肥料团块等杂物，样品原始重量不低于1 kg，取样位置0～20 cm。采样时要尽量选择在农用大田、林带地、菜地、山坡下侧土层较厚处等地，避开明显污染地块、新近搬运的堆积土、垃圾土、田埂等地段，以保证样品的代表性。

碱解氮的测定：称风干土样(通过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标准溶液滴定内室硼酸中的氨，计算碱解氮含量；有效磷的测定：称风干土样(通过2 mm筛)2.5 g于浸提瓶中，加50 mL 0.5 mol/L碳酸氢钠浸提剂，在(25±1)℃温度中振荡30 min后过滤；速效钾的测定：称风干土样(通过2 mm筛)5.0 g于浸提瓶中，加50 mL乙酸铵溶液，加塞振荡30 min后过滤。

各项元素指标测定方法见表1。

1.3 土壤肥力综合评价方法

目前针对土壤肥力的综合评价方法有很多，如指数和法[12]、灰度关联分析法[13]、投影寻踪模型法[14]、聚类分析法[15]、BP人工神经网络法[16]，本研究采用Fuzzy综合评价法对土壤肥力进行综合[17-18]，这种方法是要以隶属度函数中所求得隶属度来刻画肥力质量指标的状态界线。

1.3.1 土壤肥力指标的隶属函数选择 隶属函数模型主要包括线性模型和非线性模型，确定隶属度函数的方法有很多，如模糊统计法、例证法、分段函数法等[19]。考虑评价工作的实用性和简洁性，本次研究以线性模型为基础，采用峰值型和戒上型。

(1)峰值型：

(1)

表1 土壤样品测定方法

(2)戒上型：

(2)

综合考虑土壤中肥力指标的地球化学特点和专家经验[20]，认为土壤肥力指标中全氮、全磷、全钾、碱解氮、速效磷、有效钾、有机质隶属度函数采用戒上型，pH值的隶属度函数采用峰值型。各评价指标的隶属度函数及阈值见表1。

1.3.2 土壤肥力综合评价指标计算 IFI(Integrated Fertility Index)表示土壤肥力状况的综合评价值，是根据加法法则在相互交叉的同类指标之间由加法而合成[21]，计算公式：

IFI=∑Wi×Ni

(3)

式中，Wi——第i种评价指标的隶属度值，Ni——第i种评价指标的权重系数。

1.3.3 地统计处理方法 采用SPSS19.0对土壤肥力指标值进行正态分布性检验，相关系数的计算，通过对数转化将不服从正态分布的数据转化为正态分布，满足半方差分析的要求。

采用地统计学软件GS+9.0进行理论模型的拟合和半变异函数的计算，变异函数是研究土壤肥力指标空间变异随机性或结构性的重要方法，能够反映不同距离观测值的空间自相关程度。该函数如下：

(4)

式中：γ(h)——半变异函数；h——步长；N(h)——

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

表2 土壤肥力评价指标隶属度函数类型及其转折点

观测样点对数；Z(xi)和Z(xi+h)]——区域化变量Z(x)在空间观测点xi和xi+h处的实测值。

在本次研究中为了更加准确地反映土壤综合指标的空间变异性，根据半变异函数的决定系数R2和残差RSS进行拟合，构建理论模型[22]，其中R2是回归平方和与总平方和的比值，RSS为各采样点实际值与预测值的差的平方和，因此R2值越大，RSS值越小，模型的拟合精度越高[23]。本次研究所建立的半方差函数模型为指数(Exponential)。

γ(h)=0h=0

(5)

式中：C0——块金方差，即间距为0时的半方差，由随机因素引起的变异；C——结构方差，由系统因素引起的变异；C0/(C0+C)——块金效应，反映了系统变量的空间相关性程度，可表明土壤肥力的空间依赖性；a——变程，即半方差达到基台值所对应的距离。

由ArcGIS10.0中地统计模块的普通克里格空间插值法完成Kriging插值[24]。通过变量在点x处影响范围内的n个有效样本值Z(xi)的线性组合得到在该点x处的估计值Zx，如下：

(6)

式中：λi——赋予样本值Z(xi)的权重，表示各个样本值Z(xi)对Zx的贡献。

2 结果与分析

2.1 土壤肥力指标描述性统计分析

利用SPSS19.0对所采集的1399件土壤样品进行了描述性统计分析，结果见表3，从表中可以看出有机质含量分布在0.23～3.94 g/kg，平均值为1.3 g/kg。碱解氮、有效磷、速效钾的平均含量分别为85.62、27.20、161.81 mg/kg，变化幅度分别为6.96～441.04、4.48～371.60、27.97～2699.56 mg/kg，有效磷的最大值是最小值的82.95倍，变幅最大。从变异系数分析，一般按照变异系数的大小一般分为四级[25]，<30 %为均匀分布，30 %≤c.v.<60 %为弱变异，60 %≤c.v.<100 %为中等分异，≥100 %为强分异。有效磷变异系数为107.04 %，属于强程度变异，这可能与土壤中相关物理化学反应及磷肥的使用有关，施用的磷肥易残留在土壤中，因其利用率低、迁移速度慢，因而有效磷分布不均，变异程度大。速效钾和碱解氮的变异系数分别为76.32 %和58.40 %，属于中等程度变异；全氮、全磷、全钾、有机质、pH值的变异系数全部低于40 %，属于弱程度变异，其中pH值最低，为2.80 %，体现出在空间中分布较为稳定。

表3 研究区土壤肥力质量的评价定量指标测定结果

表4 土壤肥力的评价标准

表4为第二次土壤普查土壤养分分级标准，由此来看，本次研究中全氮、全磷、碱解氮、有效磷、速效钾、有机质含量分布在Ⅰ～Ⅴ，全钾含量分布在Ⅰ～Ⅳ，表明研究区肥力指标分布较为均匀，土壤中养分分布较为平衡。由于北京地区的土壤大部分都为碱性，因此本次研究中pH值全部为Ⅰ。

2.2 土壤肥力综合评价指标分析

2.2.1 土壤肥力指标权重值的计算 在多项肥力指标综合评价过程中，每一项肥力指标由于对结果的影响程度和重要程度不同，因而贡献率即权重也随之不同。如何合理地确定各单项肥力指标的权重成为完成土壤肥力综合评价的关键因素。目前，国内外确定土壤肥力各指标权重值的方法主要有主观经验法、专家征询法、专家打分法、指标法、主因子分析法、层次分析法等[26-28]。由于土壤是一个有机的整体，土壤肥力指标之间是相互作用并相互影响的，之间存在一定的相关关系，因此本次研究选用指标相关系数法对各指标进行权重赋值[29-30]。通过某项肥力指标与其他肥力指标之间的相关系数的平均值占全部肥力指标的相关系数平均值的综合比例作为该肥力指标的权重系数。由表5可知，全氮、全磷、碱解氮、有效磷、有机质的权重值相当，并高于全钾、速效钾、pH值的权重值，其中全钾的权重值最低，仅为2.58 %。

2.2.2 土壤肥力综合评价指标的计算 由土壤中每一项肥力指标的隶属度函数和权重值进行计算，可得土壤肥力状况的综合评价值(IFI)。IFI的取值范围在0～1，如果值越小，则表明土壤肥力水平越低，相反，如果值越大，表明土壤肥力水平越高。经过统计，发现研究区IFI在0.2～0.95，平均值为0.49，标准差为0.23，变异系数为45.98 %，由此可知肥力分布相对均匀。整体表明研究区内土壤整体肥力状况良好。

2.3 土壤肥力指标的趋势分析

经过统计分析发现全磷、有效磷、速效钾、碱解氮和IFI数据呈非正态分布，这样的数据会增大变异函数的波动性，使误差提高，因此对这四项指标的数据进行对数转换，转换后的数据服从对数正态分布，因而在计算这4项指标的变异函数值时用转换后的数据，同时在运用ArcGIS 10.2中地统计模块的趋势分析功能和克里金插值时，变换类型应当选择“Log”选项。

由于土壤与外界物质(大气、水、生物等)不断接处，自然会受这些因素的内在影响，故而土壤内相关元素含量或按一定规律计算所得的指定系数的空间分布具有较为明显的趋势效应特征，本文采用ArcGIS 10.2中趋势分析模块，通过异向性轴轴向自动搜索功能，探讨了土壤肥力各项单指标和综合指标的空间分布趋势效应特征。在讨论各项指标的各向异性特征时，依次选取常数(区域化变量沿一定方向呈常量增加或减少)、一次(区域化变量沿一定方向呈直线变化)、二次(区域化变量沿一定方向呈多项式变化)[31]。图2为研究区全氮、全钾、速效钾和IFI的空间分布趋势效应分析图，图中X轴表示正东方向，Y轴表示正北方向，Z轴表示各样点测定值的大小；右后投影面上浅蓝色线表示南北向的全局性趋势效应变化曲线，左后投影面上浅绿色线表示东西向的全局性趋势效应变化曲线。研究区全氮在南北方向上有一定程度的降低，属于常数类型，碱解氮和pH值也属于该类型；速效钾表现为在东西方向和南北方向上呈直线分布，属于一次类型；全钾表现为在南北方向上先增加后降低的趋势，属于二次类型，全磷、有效磷、有机质也属于该类型；此外土壤IFI表现为在东西方向和南北方向上呈直线分布，也属于一次类型。

表5 肥力指标的相关系数平均值和权重值

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

结合克里格插值法所造成的插值误差进行比较，综合研究所选参数的合理性，参照标准为平均误差(ME)的绝对值最接近0，若平均标准误差(ASE)大于均方根误差(RMSE)，则表示人为提高了预测值，反之则表示人为降低了预测值，故而两者应尽量相等。标准化平均误差(MSE)的绝对值尽量接近0，标准化均方根误差(RMSSE)值尽量接近1，如果标准化均方根误差大于1，说明人为降低了预测值，反之人为提高了预测值[32]。由表6可知，二次趋势预测中MSE的值为1153.77，远远大于0值，因此不予考虑，IFI的一次趋势预测与常数趋势预测相比，ME和MSE的绝对值更加接近0，ASE和RMSE最接近，的绝对值更加接近0，RMSSE的值最接近1，因此在之后克里格插值时趋势参数选择一次更加精准，同理可知，全氮、碱解氮和pH值选择常数，速效钾选择一次，其余指标均选择二次。通过运用不同趋势次数插值误差分析，佐证了肥力指标趋势分析图的结论。

2.4 土壤肥力指标的空间变异

半变异函数是地统计学中应用最广泛的空间格局描述的基本工具[33]，不但可以利用半变异函数对参数的空间分布进行结构性和变异性分析，而且利用这一结果结合克里格法可以对未知样点进行预测和模拟。块金系数(C0/C+C0)是块金值与基台值的比值，若块金系数大于75 %，则说明研究变量的空间相关性较弱，其空间变异主要受随机性因子控制，如工业污染、耕作、施肥等[34]，它会对变量空间变异的结构性起到消减作用，使土壤朝均一化方向发展[35-36]；当块金系数介于25 %～75 %说明变量具有中等强度空间相关性，其空间变异性是由随机性因素与结构性因素共同决定；当块金系数小于25 %时说明变量具有较强的空间相关性，其空间变异主要由结构性因素如成土母质、地形、地貌、气候等控制[37-39]。本次研究所得全钾的块金系数为29.73 %，具有较强的空间相关性，主要由随机性和结构性因素共同控制，其余指标的块金系数9.63 %～20 %，具有很强的空间相关性，主要由结构性因素控制。由表7可知，研究区各项指标的R2值为0.691～0.976，说明各项指标的实验半变异与理论半变异指数模型拟合效果很好。变程分布范围在468～2550，因此保障了对其进行克里格插值得到结果的可靠性[40-41]。

表6 不同趋势次数插值误差比较

续表6 Continued table 6

指标Index趋势效应Trend effect平均误差ME均方根误差RMSE平均标准误差ASE标准化平均误差MSE标准化均方根误差RMSSE二次-0.000040.12378-0.000450.127520.97442碱解氮常数0.1440245.654-0.0142046.5911.01602一次0.1416845.645-0.0152546.5261.01940二次0.1844045.798-0.0236546.4581.03993有效磷常数-1.311626.9451-0.0684818.4261.44717一次-1.312326.9452-0.0683718.4271.44677二次-1.311326.9414-0.0688218.4301.44726速效钾常数-3.9305117.80-0.0629771.5091.57570一次-4.0402117.24-0.0600070.2301.54445二次-4.0609117.25-0.0614370.2231.54452有机质常数0.001410.410240.002660.399921.02292一次0.001400.410230.002620.399071.02488二次0.000870.410870.001310.395431.03527pH值常数0.0000220.215410.000040.209281.02718一次0.0000360.215400.000110.209391.02661二次0.0000240.215450.000070.209611.02578IFI常数0.010280.209320.015450.658940.31763一次0.009290.210290.001290.295790.77613二次0.006980.26114-0.095121153.770.86881

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

注：*代表已经对数据进行了对数转换。

图3 土壤肥力综合评价得分空间分布Fig.3 Spatial distribution of soil integrated fertility scores

表8 研究区土壤不同等级面积统计

为了更加直观地展现研究区土壤肥力指标的空间分布格局，采用普通克里格插值法(Kriging)对土壤的肥力指标进行插值，其空间分布图见图3。研究区土壤整体pH为碱性，全钾含量较高，碱解氮、有效磷、速效钾和有机质含量中等，全氮和全量含量偏低。整体来看，研究区土壤肥力分布较均匀，土壤肥力高的地区主要分布在南仪垡-大耕垡一带、满庄-果村-西马坊一带、养鸡场和渠头周边。土壤肥力低的地区主要分布在研究区中部，周边也有零星分布。此处土地利用主要以耕地为主，因此，应当合理选肥、合理配肥、合理施肥。切不可盲目滥施肥，不但造成肥料的浪费和污染，而且肥料施用不合适的种类和比例用量都不利于作物生长。相关研究表明，无机肥一般用于增加土壤供肥强度，提高土壤中有效养分比例；有机肥用于增加土壤供肥容量，优化养分库容；两者按照一定比例配置结合使用则会取长补短，在补充土壤肥力和增加作物产量方面优于两者单独使用[42]。

用各指标隶属度函数平均值制作雷达图不但可以反映研究区整体肥力状况，而且能体现出哪项指标对整体肥力贡献率大。从图4可以看出，全氮、全磷、全钾、碱解氮、有效磷、速效钾、有机质隶属度函数平均值均在0.5附近波动，体现了研究区肥力指标整体分布均匀，而pH最小，仅为0.21。总体上各指标隶属度平均值所组成的多边形面积中等，说明该区土壤肥力状况处于中等水平。

本次研究中土壤肥力综合指标取值在0.2～0.95之间，参照以往的研究方法[43]，将土壤肥力综合指标取值划定为5种土壤等级，即：优质(0.74～0.95)、良好(0.57～0.73)、中等(0.45～0.56)、差等(0.3～0.44)、劣等(0.2～0.29)，并统计其分布面积，见表8。研究区主要以良好级、中等级和差等级土壤为主，分别为15.17、26.51和20.32 km2，共占比为94.86 %。优质级和劣等级分布较少，分别为0.84和2.52 km2，共占比为5.14 %。

图4 各指标隶属度函数平均值Fig.4 The membership mean values of indices

3 结论与讨论

土壤是一个非常复杂的生态载体，不但空间上存在着较大的变异性，而且内部某些元素和指标之间在特定的环境下会发生一系列生物、化学反应转化为另外一些产物，直接或者间接地影响着评价结果。因此，对于土壤肥力评价研究最终重要的就是评价指标的筛选。此次评价过程中指标筛选应遵循主导性原则、系统性原则、独立性原则、空间变异性原则，依次降低主观因素或片面原因对评价结果带来的影响。同时，各指标的权重是通过采用指标相关系数法来确定的，再与隶属度函数相结合得到土壤肥力状况系数(IFI)，进而对土壤肥力进行综合评价。

此次对于家务乡土壤肥力进行的综合评价只能代表潜在的生产能力，尤其是对人口密集的北京地区而言，一些自然因素和人为因素的影响不可忽视，下一步研究应当根据地理位置周边的实际情况，将气候、地形、施肥周期、微生物、作物种类等因素与肥力指标本身综合分析，构建一个更加全面的评价体系，这样才能更加精准地得到土壤的实际生产能力。

利用地统计学和GIS技术结合的方法能够更加全面、更加直观地了解于家务乡土壤单项和综合肥力状况。从土壤肥力综合评价得分的空间分布来看，于家务乡土壤肥力整体趋于良好，劣等土壤和优质等级土壤占比最低，仅为3.85 %和1.29 %，差等、中等和良好等级土壤占比为94.86 %；土壤综合肥力主要受到结构性因素影响的结论与变异系数普遍偏小有关。同时，对于综合肥力较低地区要合理施肥尤其是氮肥和磷肥，尽快改善土壤性能，最终达到不断增加农产品附加值的作用。