黄淮海平原农区农用地土壤肥力评价及时空变化特征——以山东省博兴县为例

李颖慧，姜小三，2*

（1.南京农业大学资源与环境科学学院，南京 210095；2.南京农业大学泰州研究院，江苏 泰州 225311）

土壤肥力是土壤质量的本质属性，可为植物生长提供必需的养分，土壤养分含量及其空间分布特征是土壤肥力的重要标志，其含量的高低直接影响粮食生产[1-3]。在现代农业生产中，随着土地利用方式、施肥模式和管理措施的改变，土壤肥力出现明显变化，受土壤类型、地形、气候、生物等因素的影响，土壤养分状况也存在较大差异[4-5]。因此，了解土壤肥力变化规律以及揭示土壤养分的时空演变特征对合理利用土壤资源和科学施肥具有重要意义。

我国土壤肥力在第二次全国土壤普查后已发生明显的变化，除东北黑土外，其他区域土壤肥力总体上虽有提升，但存在养分非均衡化、空间变异较大等特点[4，6]。孙晓兵等[7]的研究发现，河北省曲周县耕地土壤养分的最大值和最小值具有显著性差异，其比值范围为7.14～28.89。王远鹏等[8]对1982—2017年江西省进贤县红壤稻区土壤养分含量进行分析，结果表明有机质、有效磷、速效钾和碱解氮在空间上呈现不同程度的上升趋势。然而，也有研究表明长期耕作和较高施肥量条件下土壤有机质含量明显下降[5]。黄淮海平原农区土壤类型为潮土和盐渍土，是我国重要的农业土壤之一，同时也是重要的粮食生产基地，小麦-玉米复种连作是该区域主要的种植制度，复种指数高，在粮食连年高产的背后，却存在生态环境可持续性显著不足、速效养分含量高和有机质含量缺乏等问题[9-11]。山东省博兴县作为黄淮海平原粮食高产的典型县域之一，随着连续多年农业生产要素的大量投入，长期以来人类对土地资源改造规模和强度在不断增大、增强，使得土壤养分关系发生较大变化，生物及人类农业生产活动增强了土壤养分的空间异质性，原有的土壤普查数据已不能代表如今的土壤现状。因此，迫切需要研究近年来的耕地土壤肥力变化情况，摸清农业生产和管理模式对土壤养分的影响，从而为土壤资源的可持续利用提供科学依据。

本研究选取山东省博兴县2007年和2020年土壤pH、有机质（Organic matter，OM）、有效磷（Available phosphorus，AP）、速效钾（Available potassium，AK）、碱解氮（Alkali-hydrolyzable，AN）等5 种受人为因素影响较敏感且对作物生产具有限制作用的养分指标[12-13]，运用地统计学方法与指数和法对土壤肥力进行综合评价，揭示农用地土壤养分的时空变化和分布特征，并探讨土壤养分及土壤肥力在时间和空间上产生变化的原因，以期为土壤资源的管理和利用提供科学依据。

1 材料与方法

1.1 研究材料

1.1.1 研究区概况

研究区位于山东省滨州市博兴县（118°02′19″～118°22′54″E，36°57′42″～37°23′04″N），地处黄淮海平原东部，山东省北部、滨州市东南部，总面积900.7 km2，下辖3个街道、9个镇。博兴县地势平坦，总体趋势呈现南北高、中间低、西高东低，平均海拔为8 m，属于温带大陆性季风气候，四季分明，年平均气温13.6 ℃，无霜期191 d，日照时间2 479.8 h。境内有黄河、小清河、支脉河三大水系，年平均降水量582.8 mm，主要集中在6—8 月。成土母质为黄河冲积物，土壤类型主要为盐土和潮土，土壤质地有砂壤、轻壤、中壤和重壤。博兴县第一产业总产值中，农业产值位居第一，占比为57.3%，主要种植作物为玉米、小麦、棉花，作物多为一年两熟，土地的集约化利用程度较高。自2007—2020年博兴县粮食作物种植面积由6.58×104hm2增加到7.40×104hm2，粮食单产由5 445 kg·hm-2增加到7 185 kg·hm-2[14]。

1.1.2 数据来源

本研究所需的基础图件及土壤肥力数据均由博兴县农业农村局提供。利用ArcGIS 10.2 软件以土地利用现状图为底图设置规则网格进行采样点的布设，如图1 所示，每个采样点用GPS 实地坐标定位后，按梅花型布设五点采集0～20 cm 土样，混合均匀后用四分法取约1 kg土样装入聚乙烯密封袋内，带回室内自然风干，用玛瑙研钵进行研磨后过筛备用。

图1 研究区采样点分布图Figure 1 The location of soil sampling sites

土样测定分析所用试剂、标准溶液配制以及仪器工作参数均遵循相关技术要求规范，分析数据由山东省农业农村局提供。其中，有机质含量测定采用重铬酸钾外热源法，碱解氮含量测定采用NaOH 碱解扩散法，有效磷含量测定采用NaHCO3浸提-钼锑抗吸光光度法，速效钾含量测定采用醋酸铵浸提-火焰光度法，pH值测定采用电位法[15-16]。

1.2 研究方法

1.2.1 土壤肥力综合评价

（1）构建隶属度函数。由于选取的土壤肥力指标具有多样性，且各指标间的量纲具有差异，需要采用隶属度函数对各指标进行标准化处理，使得各个指标之间具有可比性和统一性。隶属度函数是能够反映评价指标与作物生长效应之间关系的数学模型，本研究选取的隶属度函数理论模型为S型和抛物线型。

采用S 型隶属度函数的指标包括土壤有机质、有效磷、速效钾、碱解氮，其函数表达式如下：

采用抛物线型隶属度函数的指标为pH，其函数表达式如下：

式中：x1、x2、x3、x4表示隶属度函数的拐点值，参考文献[17]以及农用地等级划分标准[18]，同时结合博兴县土壤的实际情况，确定隶属度函数的拐点值，见表1。

表1 土壤肥力指标隶属度函数拐点值Table 1 The turning points of membership function of soil fertility indexes

（2）确定指标权重。指标权重能够反映出各个评价指标对事物影响程度的大小，权重的计算方法有特尔菲法、层次分析法、主成分分析法、相关系数法等。其中相关系数法能够客观地对评价指标赋予适宜的权重，避免了人为主观因素的干扰，故本研究采用相关系数法计算权重[19]。相关系数法确定权重首先要计算各评价指标之间的相关系数，然后计算某一评价指标与其他评价指标的相关系数的平均值，最终获得该单项指标相关系数平均值与所有指标相关系数平均值的比值，即为该单项指标的权重。本研究将2007年和2020年两期指标的权重平均值作为土壤肥力指标的最终权重[20]，得到pH、有机质、碱解氮、有效磷、速效钾权重值分别为0.075、0.280、0.259、0.165、0.221。

（3）土壤肥力评价。指数和法是划分土壤肥力等级最普遍的方法[21]，土壤肥力评价指数是以模糊数学中的加法和乘法原则为基础，将每个评价单元中各评价指标的隶属度函数值与权重值分别相乘后求和，得到土壤肥力评价指数（IFI），计算公式如下：

式中：IFI表示土壤肥力的综合评价指数；fi表示第i个评价指标的隶属度函数值；ai表示第i个评价指标的权重值。

土壤肥力等级的划分，基于《全国耕地质量等级评价指标体系》中的分级标准，并根据黄淮海平原农区农用地的实际情况将研究区土壤肥力划分成五个等级（表2），Ⅰ～Ⅴ级土壤肥力水平依次降低。

表2 土壤肥力等级划分标准Table 2 The classification standard of soil fertility degree

1.2.2 统计分析与制图

（1）地统计分析。地统计分析是以区域化变量为基础，借助变异函数来研究空间数据的结构性与随机性、空间相关性与依赖性、空间格局与变异，并对这些数据进行最优无偏内插估计。半方差函数能够反映区域内各距离观测值间土壤性质的变化，是研究土壤养分空间变异规律的有效工具[22]。克里金插值作为一种最优无偏估计的方法，基于采样数据计算半方差函数来反映区域化变量结构信息，在每种土壤养分含量符合正态分布的情况下，可以模拟每种土壤养分的空间结构[23]。半方差函数能够说明局部地区变量的空间变异及其相关程度，在半方差函数理论模型中，块金值表示区域变量随机性部分。当采样点的距离h增大时，半方差函数γ（h）从初始的块金值达到一个相对稳定的常数，该常数称为基台值，此时的h称为变程。块金效应指块金值与基台值的比值，表示随机性因素引起的空间异质性占系统总体变异的比例，当比值<25%时，变量具有强烈的空间自相关性，空间差异性主要受结构性因素影响；比值在25%～75%之间时，变量具有中度的空间自相关性，空间差异性受结构性和随机性因素共同影响；当比值>75%时，变量空间自相关性较弱[24]。

区域变化量Z（x）在空间位置xi和xi+h的值为Z（xi）和Z（xi+h），两者差值的方差的一半被称为区域变化量Z（x）的半方差函数，计算公式如下：

式中：γ（h）为半方差函数；Z（x）为区域变化量，满足二阶平稳假设；h为两个样本点的分隔距离；N（h）为以h为间距的所有观测点的成对数目；Z（xi）和Z（xi+h）分别为Z（x）在位置xi和xi+h的实测值。

在计算半方差函数前需要对样本数据进行正态检验，并对不符合正态分布的数据进行转换处理。

（2）数据统计与制图。用Excel记录数据并统计，通过GS+软件进行函数计算和模型拟合，运用SPSS对原始数据和土壤肥力指数进行描述性统计分析，描述性统计只能揭示数据的点位特征，需要进一步利用空间插值的方法对博兴县土壤养分和土壤肥力的空间分布情况进行直观表达。

2 结果与分析

2.1 描述性统计分析

2.1.1 土壤养分描述性统计分析

2007年和2020年博兴县土壤养分的统计特征值见表3，由表3 可知，2007 年土壤有机质和碱解氮、有效磷、速效钾的极差分别为28 g·kg-1和138.50、97.00、460.00 mg·kg-1，平均值分别为15.13 g·kg-1和103.20、39.07、164.73 mg·kg-1。2020年土壤有机质和碱解氮、有效磷、速效钾的极差分别为18.94 g·kg-1和132.40、64.00、486.00 mg·kg-1，平均值分别为17.48 g·kg-1和84.73、24.99、242.80 mg·kg-1。除土壤pH 值外，同种土壤养分的含量在相同年份差异化明显。对比2007年发现，2020 年土壤有机质、碱解氮和有效磷的极差呈降低的趋势，有机质和速效钾含量的平均值升高，增幅分别为15.53%、47.39%，碱解氮和有效磷含量的平均值降低，降幅分别为17.90%、36.04%。参照全国第二次土壤普查土壤养分的分级标准可知，土壤有机质含量在2007 年和2020 年处于中下等水平（10～20 g·kg-1）的样点数最多，占总样点比例分别为44.56%和60.00%；有效磷含量在2007 年和2020 年处于较丰富水平（20～40 mg·kg-1）的样点数最多，占总样点比例分别为53.89%和48.00%；速效钾含量在2007 年处于中上等水平（100～150 mg·kg-1）的样点数最多，占总样点比例为51.22%，在2020年处于丰富水平（>200 mg·kg-1）的样点数最多，占总样点比例为54.00%；碱解氮含量在2007 年处于中上等水平（90～120 mg·kg-1）的样点数最多，占总样点比例为48.22%，在2020年处于中下等水平（60～90 mg·kg-1）的样点数最多，占总样点比例为54.00%；土壤pH 值在2007 年和2020 年均处于碱性水平（7.5～8.7）。CV通常用于表示变量的空间变异程度，CV<10%表示弱变异性；CV处于10%～90%之间表示中等变异性；CV>90%表示强变异性[24]。结果表明，仅土壤pH 值为弱变异性，其余土壤养分指标都表现为中等变异性。K-S（P>0.05）检验表明土壤pH 值和速效钾呈现对数正态分布，其余指标均为正态分布。

表3 土壤养分指标描述性统计特征Table 3 The descriptive statistics of soil nutrient index

由表4 可知，2020 年土壤pH 与土壤碱解氮呈极显著负相关（P<0.01），与有效磷呈显著负相关（P<0.05）；2007 年和2020 年土壤有机质与碱解氮、速效钾均呈极显著正相关（P<0.01），土壤碱解氮与速效钾呈极显著正相关（P<0.01）。2007 年土壤有效磷与速效钾呈极显著正相关（P<0.01）。

表4 2007年和2020年土壤养分指标之间的相关系数Table 4 The correlation coefficient of soil nutrient indexes in 2007 and 2020

2.1.2 土壤肥力指数描述性统计分析

表5 表明，2007 年和2020 年土壤肥力指数平均值分别为0.56、0.54，由表2中的土壤肥力等级划分标准可知，该地区的土壤肥力水平总体表现为中等（Ⅲ级）水平。2007 年和2020 年土壤肥力指数分别为0.17～0.99、0.10～0.94，极差分别为0.82、0.84，同年份的土壤肥力指数变化范围较大。具体来看，2007 年土壤肥力指数的最小值比2020 年高0.07，最大值比2020 年高0.05，但2007 年土壤肥力平均值比2020 年高0.02，初步判定2007—2020 年期间不同级别的土壤肥力呈现出不同程度的变化。土壤肥力指数变异系数由39.16%下降到38.77%，但下降幅度小且均为中等变异。

表5 土壤肥力指数（IFI）描述性统计特征Table 5 The descriptive statistics of soil IFI

2.2 地统计分析

2.2.1 土壤养分的变异性分析

对2007年和2020年土壤养分各指标进行半方差函数理论模型拟合，结果见表6。2007年和2020年土壤有机质和速效钾的块金效应在28.64%～49.71%之间，表明其均具有中度的空间自相关性。2007 年土壤pH 的块金效应为11.45%，呈现强烈的空间自相关性，2020 年升高至50.00%，呈现中度的空间自相关性。虽然在2007年和2020年土壤碱解氮和有效磷均呈现出空间自相关性，但对比2007年发现，2020年土壤碱解氮的块金效应由10.82%上升至18.13%，土壤有效磷的块金效应由13.77%上升至19.61%，块金效应均呈现出上升趋势。这可能是由人为因素造成的，说明其随机性变异逐渐增强，区域内空间变异一般受到施肥、管理水平、作物种植类型、灌溉条件等随机因素的影响。由此可知，土壤碱解氮、有效磷、有机质、速效钾和pH的空间结构变异减弱、随机变异增强，这是由结构性因素和随机性因素共同决定的。

表6 土壤养分变异函数理论模型及相关参数Table 6 Variogram theory models and corresponding parameters of soil nutrient

2.2.2 土壤养分的空间分布特征

基于地统计学的半方差函数理论模型及其相关参数，采用普通Kriging 插值法绘制土壤养分（pH、OM、AN、AP、AK）和土壤肥力指数空间分布图（图2），图2表明，2007年和2020年土壤pH 值为7.75～7.83和8.22～8.30，面积占比分别为73.02%、51.97%，2007 年土 壤pH 最大值为8.00，2020 年土壤pH 最小值为7.94，表明土壤pH 总体略微增长。2007 年土壤有机质≤16.5 g·kg-1的面积占比最大，为62.58%，2020年土壤有机质>16.5 g·kg-1的面积占比最大，为78.55%，有机质含量总体呈现增长的趋势，但增幅较低。2007年乔庄镇、庞家镇、陈户镇和纯化镇土壤有机质≤14.5 g·kg-1，2020 年，在以上四个镇中，除乔庄镇部分区域土壤有机质含量≤14.5 g·kg-1外，其余区域的土壤有机质含量>14.5 g·kg-1。2007 年土壤碱解氮含量>92.5 mg·kg-1的面积最大，有效磷含量>34.5 mg·kg-1的面积最大，占比分别为78.85%、68.16%；2020 年土壤碱解氮和有效磷含量分别为75.5～92.5 mg·kg-1和25.5～34.5 mg·kg-1的面积最大，占比分别为81.35%、68.16%，表明2020 年土壤碱解氮和有效磷均大幅度降低。2007年和2020年土壤速效钾含量均主要分布在150～250 mg·kg-1之间，面积占比分别为50.34%、53.65%；与2007年相比，2020年速效钾含量>350 mg·kg-1的面积占比增加了20.82%，高值区（>350 mg·kg-1）基本分布在城东街道、博昌街道和锦秋街道。综上所述，2007 年和2020 年研究区土壤养分在空间分布上存在一定的差异性，土壤碱解氮和有效磷含量均明显降低，土壤pH和有机质含量有所升高，但增幅较小，土壤速效钾含量在不同范围内呈现不同程度的增长。

图2 土壤养分指标空间分布Figure 2 The spatial distribution of soil nutrient indexes

续图2 土壤养分指标空间分布Continued figure 2 The spatial distribution of soil nutrient indexes

土壤肥力指数空间分布见图3，2007 年博兴县土壤肥力以Ⅲ级为主，面积占比为31.74%，主要分布在东南和东北部地区；Ⅰ级的面积占比为19.71%，基本分布在研究区的中部和西南地区；Ⅱ级的面积占比较小，为17.86%，主要分布在曹王镇；Ⅴ级的面积占比为12.38%，分布在乔庄镇的北部、纯化镇的东北部和庞家镇的西北部地区。土壤肥力水平在整个研究区的空间分布由低到高依次为北部地区<南部地区<中部地区，各等级地块呈斑块状零星分布在各区域内。2020 年土壤肥力处于Ⅰ级的农用地面积占比为6.09%，以Ⅱ级和Ⅲ级为主，面积占比分别为35.64%、27.60%，其中，除纯化镇和吕艺镇部分区域土壤肥力处于Ⅱ级水平外，Ⅱ级基本分布在研究区的西北和西南地区；Ⅴ级面积占比为11.39%，主要分布在乔庄镇和店子镇。土壤肥力水平在整个研究区的空间分布由低到高表现为东部地区<西部地区。与2007 年对比发现，2020 年土壤肥力Ⅰ级的面积占比减少了13.62 个百分点，店子镇的西北部、湖滨镇的东北部、吕艺镇的西南部和城东街道土壤肥力由Ⅰ级降低为Ⅱ级或Ⅲ级。除纯化镇少数地区土壤肥力为Ⅳ级外，庞家镇、陈户镇和纯化镇土壤肥力均由原来的Ⅲ级及以下水平提升至Ⅲ级及以上水平。店子镇和兴福镇土壤肥力等级呈现出不同程度的退化，基本由Ⅲ级及以上水平退化至Ⅲ级及以下水平。曹王镇的土壤肥力等级基本保持不变。对比2007年和2020年两个时期的土壤肥力水平，总体来看，同一区域内肥力水平差异呈减小的趋势，向均一化方向发展。

图3 土壤肥力指数空间分布Figure 3 The spatial distribution of soil fertility index

3 讨论

3.1 土壤养分变化分析

土壤养分中氮、磷、钾是作物生产的主要限制因素，土壤有机质能够改善土壤的物理性质，同时为植物生长提供大量的养分元素，土壤pH 值会影响植物对其他元素的吸收[6，25]。2007 年和2020 年博兴县土壤pH均呈碱性，这可能与土壤类型主要为潮土有关，潮土的pH值通常表现为中碱性水平。该地区成土母质为黄河冲积物，近年来，受到有机肥来源和价格等因素的影响，农户对有机肥的施用量明显减少[26]，博兴县土壤有机质含量基本维持在16.5～20.5 g·kg-1之间，大部分地区土壤有机质处于中下等水平，同时该地区菜地多采用设施农业栽培方式，且复种指数高，共同导致土壤中有机质的缺乏。碱解氮和有效磷的含量降幅较大，速效钾含量则大幅增加，这与当地的土壤类型、施肥方式、管理模式等自然因素和人为因素有关[27-28]，根据实地调查数据，博兴县自2007 年开始虽采取了测土配方施肥的管理模式，但农户仍然按照经验进行施肥，其实际氮肥施用量由275.68 kg·hm-2降至243.03 kg·hm-2，磷肥施用量由116.63 kg·hm-2降至85.28 kg·hm-2，氮肥和磷肥的投入量降低，长期耕种导致土壤氮和磷逐渐减少[7]；同期钾肥施用量由93.45 kg·hm-2提升至107.48 kg·hm-2，潮土主要由富钾矿物云母等组成[29]，使得土壤中全钾含量较高，钾肥的投入量增加也可能是造成土壤速效钾含量升高的原因。2007 年和2020 年土壤碱解氮、有效磷和速效钾含量均处在中上等水平，这可能是由于土壤pH 值较高，增加了土壤中氮、磷、钾的有效性[29-30]。2016 年我国氮肥、磷肥、钾肥的平均施用量分别为138.64、49.81 kg·hm-2和38.22 kg·hm-2[31]，而博兴县的施肥量约是我国平均施肥量的1.5 倍，除土壤本底影响外，农户的施肥方式也是影响土壤养分含量的主要因素。因此，在以后的农业土壤养分管理过程中应注意增施有机肥料，适当控制化肥的投入。

3.2 土壤肥力变化分析

土壤肥力是决定农业生产可持续发展的基础，土壤肥力状况因土壤养分元素的差异性而具有复杂的时空变化特征[13]。博兴县土壤肥力等级总体表现为北部最低、中部最高、南部介于二者之间，这可能与盐碱地基本分布在博兴县的北部地区有关，盐渍化土壤具有肥力低、有障碍层和营养元素缺乏等特点[32]。两时期对比发现，2020 年博兴县偏北部地区的庞家镇、陈户镇和纯化镇的土壤肥力基本上升至Ⅲ级及以上水平，而店子镇和兴福镇的土壤肥力等级呈现不同程度的下降，该变化除了与土壤养分元素的本底属性不同有关之外，还与土地利用和管理等因素密切相关。对比2007 年和2020 年土地利用现状（图4）发现，2007 年城镇用地面积为8 794.13 hm2，零散均匀地分布在研究区内，随着经济的快速发展和城镇化进程的不断加快，2020 年城镇用地面积增加至13 478.65 hm2，增加面积主要分布在店子镇的西北部、湖滨镇的东北部、吕艺镇的西南部和城东街道，该区域农用地面积占比少、人口密集、地块零散、集约化程度不高；而庞家镇、陈户镇和纯化镇主要受农作活动影响，其他人为活动干扰较少，方便对土地进行集约化利用和规模化管理，适宜进行大面积的作物耕种。除此之外，近几年两区（粮食生产功能区和重要农产品生产保护区）地块的建设、管理和保护得到加强也可能是该区域土壤肥力升高的原因之一。康日峰等[33]研究发现，黑土区经过10～26 年农民常规施肥管理的耕作，土壤有机质、全氮、碱解氮、有效磷和速效钾含量呈逐年上升趋势[33]，但华北平原属于典型的黄河下游地区，近10 年来在社会经济和自然因素的共同作用下，耕地土壤养分综合水平退化[26]。博兴县成土母质为黄河冲积物，土壤类型主要为盐土和潮土，土壤质地有砂壤、轻壤、中壤和重壤，不同区域土壤养分存在显著差异，加上土地利用方式和管理模式的变化，共同影响着土壤肥力。

图4 2007年和2020年博兴县土地利用现状Figure 4 Current land use map of Boxing County in 2007 and 2020

4 结论

（1）与2007 年相比，2020年博兴县土壤有机质和速效钾的平均含量增加15.53%和47.39%，碱解氮和有效磷的含量降低17.90%、36.04%，土壤pH值由7.79升高至8.28。土壤肥力各指标的空间结构变异减弱、随机变异增强，这是由结构性因素和随机性因素共同决定的。

（2）与2007 年相比，2020 年博兴县土壤肥力Ⅰ级、Ⅲ级和Ⅴ级地的面积占比分别下降了13.62、4.14、0.99 个百分点，Ⅱ级地的面积占比提高了17.78个百分点，Ⅳ级地的面积基本保持不变。

（3）博兴县2007 年土壤肥力水平在空间分布上整体表现为北部地区<南部地区<中部地区，各等级地块呈斑块状零星分布在各区域内；2020 年土壤肥力水平表现为东部地区<西部地区，同区域内土壤肥力基本处于同一水平。与2007 年相比，2020 年同区域内的土壤肥力水平差异有所减小，表现为均一化发展趋势。