1. 引言
土壤有机质是土壤重要理化性状，是土壤性能和肥力 [1] 的重要体现。土壤有机质可为植物的生长发育提供营养；改善土壤的物理性质 [2]，增强其通水透气性；为土壤微生物提供能量来源，促进土壤的生物化学过程，继而又加速有机质的分解。研究土壤有机质含量的时空变异特征，掌握其空间分布规律，可为合理种植、科学施肥、实现农业现代化提供有效依据 [3]。
目前对土壤有机质的研究多集中在耕地 [4] [5]、林地 [6] 等。研究区域主要集中在黄土高原区、黄淮海平原区、西北干旱半干旱区等。研究方法主要运用传统统计学和地统计学 [7] 结合GIS技术两种方法 [8]。杨刚等 [9] 运用传统统计学方法对成都平原水田的有机质进行分析，得出其属于中等偏上水平；解文艳等 [10] 运用传统统计学方法分析山西省潇河流域的1982年和2017年的有机质含量，探究35年间整体呈上升趋势。传统统计学方法只能分析土壤有机质静态的变化特征 [11]，对于动态的时空分布特征，则需要运用地统计学结合GIS技术的方法，例如，刘东海等 [12] 对武汉市汉南区土壤有机质的时空变异特征进行研究；张世熔等 [13] 研究了黄淮海冲积平原区的河北省曲周县土壤有机质的时空变异特征，为实现土地资源的持续利用提供参考。目前针对某地区土壤有机质的研究年限多为两个时期的整体变异特征的研究，但对于多个时期的研究较少，不利于系统地掌握研究区有机质变化的动态特征，亟需加强研究。
山东省乐陵市是农业大市、重要粮食产区，非常重视土壤肥力及长期动态监测。2005年以来，每年设置样点、采集土壤样品、分析土壤有机质等理化性状，但未进行长时间序列动态时空特征分析，不能很好地指导农业生产。基于此，本文以乐陵市为研究区，以土壤有机质含量为研究对象，分2005年、2010年、2015年和2020年4个时段，运用经典统计学、地统计学和GIS空间分析方法，研究土壤有机质含量时序变化、空间变异，揭示其时空动态变化特征，以期为乐陵市土壤改良、农业生产和生态保护提供依据 [14]。
2. 材料与方法
2.1. 研究区概况
乐陵市位于山东省西北部，为黄河下游冲积平原，地势平坦。境内有漳卫新河、马颊河、德惠新河等，水资源丰富。属暖温带湿润季风气候区，年平均气温13.6℃，年降水量610.3毫米。土壤有潮土、盐土、风沙土三大类，耕层质地以壤土为主，耕性优良，保肥透气性好，适宜农业生产。土地总面积1173.06平方千米，其中耕地面积为663.34平方千米，占全市土地总面积的56.55%。农作物主要为小麦和玉米。辖4街道9镇3乡，截止2020年11月1日 [15]，常住人口55.92万人。2020年全市实现生产总值(GDP) 242.90亿元，其中第一产业增加值为41.39亿元，占生产总值的17.04%。
2.2. 数据来源
土壤有机质含量数据来源于2005年起的乐陵市耕地地力调查、测土配方施肥和耕地质量等级评价等项目。利用GPS定位，土壤样品采样深度为0~20 cm，土壤有机质含量测定方法为重铬酸钾容量法。
2.3. 研究方法
2.3.1. 经典统计学方法
本文利用经典统计学方法进行不同年份土壤有机质含量的异常值分析、正态分布检验、变异性分析和时序变化分析。
首先对不同年份的土壤有机质含量进行异常值识别。利用SPSS 24.0采用三倍标准差检验法对各年份数据的异常值进行修正，保留平均值 ± 3倍标准差内的值 [16]，在此范围外的值用正常的最大值或最小值代替。识别出异常值并对其按上述方法进行处理，以剔除异常值。
数据符合正态分布才能进行普通Kriging插值分析 [17]。对各年份土壤有机质含量原始数据进行K-S检验，以探究数据的正态分布性，为后续的插值分析做准备。当P > 0.05时，说明数据的分布特征符合正态分布；若原始数据不符合正态分布，需对数据进行对数或平方根转换，使其符合正态分布，才能进行插值。
变异系数(Coefficient of variation, CV)的大小反映数据的离散程度，当CV ≤ 10%时，为弱变异性；当10% < CV < 100%时，为中等变异性；当CV ≥ 100%时，为强变异性。本文采用CV分析不同年份土壤有机质含量的变异性。
2.3.2. 地统计学方法
本文利用地统计学方法分析土壤有机质含量的空间趋势、最佳拟合模型、空间相关性，在此基础上进行空间插值，得到空间分布图。
利用ArcGIS 10.3的地统计分析模块中的趋势分析功能研究各年份土壤有机质含量在空间方向(经度和纬度方向)的总体规律，得到空间趋势图，以便在空间插值时去除此趋势，更准确的描述有机质的空间分布特征。ArcGIS 10.3地统计中趋势分析有5种类型：变量没有趋势效应记做“0”，按常量增加或减少记为常量；采样点的投影线呈直线增加或减少为一阶趋势；投影线呈“U”型增加或减少为二阶或多阶趋势。
半方差函数是地统计学中研究土壤变异性的关键函数。利用GS+7.0软件探究各年份土壤有机质含量半方差函数的最佳拟合模型。通过对比各类型模型的决定系数(R2)和残差(RSS)，按照R2最接近于1和RSS最接近于0的原则选择半方差函数的最佳理论模型，探究出后续空间插值工作中的模型类型。
块金比值Co/(Co+C)反映区域变量的空间自相关性程度，该比值反应随机部分引起的变异占系统总变异的比例。空间相关性受研究区内自然和人为因素的影响。一般认为，自然等结构性因素导致强的空间相关性，人为等随机性因素使得土壤有机质的空间相关性减弱。Co/(Co+C)) ≤ 25%，表示变量具有强烈的空间自相关性，其分布主要受结构性因素如地形、气候、土壤类型、植被等自然因素影响；25% < Co/(Co+C)) < 75%，表示变量具有中等的空间自相关性，受结构性因素和随机性因素共同作用影响；Co/(Co+C)) ≥ 75%，表示变量的空间自相关性弱，主要受随机性因素如灌溉、施肥、土地利用方式、耕作制度等人为因素影响 [18]。变程(Ao)表示区域化变量的空间自相关的范围程度，在变程范围之内的具有空间自相关性，反之则空间自相关性不存在。研究有机质数据的变程可为设置采样点间距大小提供有效依据。
利用ArcGIS 10.3软件，以有机质含量为数据字段，采用普通Kriging方法，选择相应的变换类型以及趋势移除阶数，基于GS+7.0得到的最佳半方差函数理论模型及各种参数值，绘制各年份土壤有机质含量的空间分布图。本文参考《中国土壤》有机质含量等级分类(表1) [19] 及结合当地的实际情况，将有机质含量分为5级(表2)。
Table 1. Classification Standard of soil organic matter content in China
表1. 中国土壤有机质含量分级标准
Table 2. Classification standard of soil organic matter content
表2. 土壤有机质含量分级标准
3. 结果与分析
3.1. 土壤有机质时序含量变化特征
利用SPSS 24.0中Kolmogorov-Smirnow对乐陵市不同年份的土壤有机质含量进行检验，剔除并替代异常值。经处理后，2005年、2010年、2015年和2020年样点分别为1029个、412个、561个和178个。经K-S正态分布检验(P > 0.05)，2005、2010、2015年的数据均符合正态分布，2020年的数据经过对数变换后也符合正态分布(表3)。
Table 3. Descriptive statistical results of soil organic matter content
表3. 土壤有机质含量描述性统计结果
2005年土壤有机质含量为5.4~18.9 g∙kg−1，均值为11.81 g∙kg−1；2010年土壤有机质的范围为5.7~17.4 g∙kg−1，均值为11.88 g∙kg−1，年变化率为0.014 g∙kg−1；2015年土壤有机质的范围为5~20 g∙kg−1，均值为13.11 g∙kg−1，年变化率为0.246 g∙kg−1；2020年土壤有机质的范围为6.3~26 g∙kg−1，均值为14.88 g∙kg−1，年变化率为0.354 g∙kg−1；2005~2020年间，年平均增速为0.20 g∙kg−1，也可以看出，年变化率总体也呈上升趋势。根据中国土壤有机质含量分级标准(表1)，整体上，乐陵市土壤有机质含量处于中等水平。近15年间，乐陵市土壤有机质明显提高，说明随着农业机械化水平提升及政府和农民对土地的保护意识的提高，作物根茬、秸秆还田等方式增加了农业土壤有机物质的投入来源，研究区土壤的固碳作用明显提高。不同年份的有机质含量的变异系数分别为17%，16.73%，30.37%，34.16%，均在10%~100%之间，揭示了土壤有机质的离散程度均属于中等变异强度。
3.2. 土壤有机质含量趋势效应
利用ArcGIS 10.3的地统计分析模块得到土壤有机质含量分布的趋势效应图(图1)，可为变异函数模型和空间分布插值方法的选择提供依据。X轴表示正东方向，Y轴表示正北方向，Z轴表示土壤有机质含量。从图中可以看出，2005年土壤有机质含量在南北方向中部略有下凹，在东西方向呈直线式缓慢增加；2010年土壤有机质含量在南北方向和东西方向均呈倒“U”型，呈先缓慢增加再缓慢减少的趋势；2015年土壤有机质含量在南北方向和东西方向均呈“U”型，先减少再增加；2020年土壤有机质含量在南北方向呈先缓慢增加，后急速减少的趋势，在东西方向呈指数型增加的趋势。
(a) 2005年 
(b) 2010年 
(c) 2015年 
(d) 2020年
Figure 1. Spatial trend map of soil organic matter
图1. 土壤有机质空间趋势图
3.3. 土壤有机质含量空间变异特征
利用GS+7.0探究土壤有机质含量的最佳半方差函数模型及其相关参数(表4)。根据决定系数(R2)最接近于1和残差(RSS)最接近于0的原则，得到2005年、2015年、2020年的土壤有机质含量最佳变异函数模型均为指数模型，其中2015年的R2最高，为0.953，拟合效果最好；2005年和2020的RSS分别为0.883和0.855，拟合效果较好。2010年土壤有机质含量的最佳变异函数模型为高斯模型，R2为0.820，拟合效果较好。
块金比值(Co/(Co+C))反映变量的空间自相关性。2005年、2010年和2015年土壤有机质含量的块金比值分别为12%、14%、12%，均小于25%，说明变量具有强空间自相关性，其分布主要受结构性因素的影响。2020年的土壤有机质块金比值为31%，属于中等强度的空间自相关性，受结构性因素因素和随机性因素的共同影响，说明随着年份的增加，灌溉、施肥等一系列人为因素对土壤有机质含量的变化具有一定的影响。块金比值说明四个时期的土壤有机质含量均适合采用普通Kriging空间插值法。
变程(Ao)表示区域变量空间自相关性的范围。2005年和2010年的土壤有机质含量的变程分别为1980 m和1957 m，均在2000 m以内；2015年的数据变程为2700 m；2020年的数据变程为14970 m。总的来看，乐陵市变程多在2000 m范围内。土壤有机质含量数据在相应的变程范围内才具有空间自相关性，以此可作为采样设计的依据。
Table 4. Characteristic parameters of semivariance function of soil organic matter
表4. 土壤有机质半方差函数特征参数
以下为根据确定的最佳半方差函数模型及相关参数绘制的半方差函数模拟模型(图2)。
3.4. 土壤有机质含量时空特征
利用普通Kriging插值法，结合半方差函数理论模型、趋势效应以及相关参数的设定，得到土壤有机质等级空间分布图(图3)。
(a) 2005年 
(b) 2010年 
(c) 2015年 
(d) 2020年
Figure 2. Semi variance function simulation model of soil organic matter content
图2. 土壤有机质含量的半方差函数模拟模型
(a) 2005年 
(b) 2010年 
(c) 2015年 
(d) 2020年
Figure 3. Spatial distribution map of soil organic matter grade
图3. 土壤有机质等级空间分布图
从图3可以看出，2005年乐陵市土壤有机质含量主要以10~12 g∙kg−1和12~15 g∙kg−1等级为主，在全域内均匀分布，分别占总面积的56.8%和39.7%，有较少 < 10 g∙kg−1等级的部分集中在中部和西北部地区；2010年土壤有机质含量12~15 g∙kg−1等级所占的面积最大，占总面积的57.4%，主要集中分布在中东部和西南地区，其次是10~12 g∙kg−1，占总面积的31.9%，主要分布在中部地区，从2005年的有机质含量以10~12 g∙kg−1为主到2010年以12~15 g∙kg−1为主，土壤有机质整体提高了一个等级，说明这5年间乐陵市科学种田、合理布局、增施肥料、种养结合等使农业土壤条件有了一定的提高；2015年全域土壤有机质以12~15 g∙kg−1等级为主，占总面积的47.1%，15~20 g∙kg−1等级的部分较之前年份占比明显增加，占比14.35%，主要以东部地区最为集中；2020年乐陵市土壤有机质仍以12~15 g∙kg−1等级为主，占总面积的53.22%，15~20 g∙kg−1等级的部分主要分布在东部地区，占总面积的19.04%，并且还有少量 > 20 g∙kg−1等级的部分，<10 g∙kg−1等级的部分主要集中在西北部和北部地区。
可以看出，2005~2020年乐陵市土壤有机质的高等级部分逐渐向东部地区和西南地区迁移，说明东部和西南地区土壤有机质积累较好，与当地的有机肥和化肥配施增多，农田管理技术更加合理有关；低等级部分逐渐向西北地区集中，其地形较高，主要为砂土，保水保肥性差。总的来看，乐陵市全域土壤有机质含量逐年提高，这与近年来山东省积极推行秸秆还田和配方施肥工作，使得有机质积累情况较好有关。
4. 结论与讨论
4.1. 讨论
1) 近15年乐陵市土壤有机质呈增加趋势，说明随着秸秆还田、合理増施化肥和有机肥等，土壤的固碳作用明显增强。随着时间的推移，块金比值逐渐增大，随机性因素如灌溉、施肥、土地利用方式等人为因素对土壤有机质含量变化的影响越来越突出，应注意人类活动对土地的影响，对乐陵市的土壤有机质进行长期监测和培肥管理，注重养地用地相结合，以促进农业的可持续发展。
2) 整体来看，乐陵市土壤有机质含量属于中等肥力水平，还需要积极的育土培肥。像东部和西南地区，有机质含量相对高的区域，土壤较肥沃，有利于农业耕作，但也要注意养地和用地相结合；对于西北地区，有机质含量低，可考虑科学施肥或调整产业结构。
3) 由于4个时段的土壤有机质样点个数不同，是否影响插值结果，应怎样设计采样方案，既保证区域有机质含量特征的提取，又能节省人力、物力，值得进一步探讨。
4.2. 结论
1) 乐陵市土壤有机质含量随年份的增长而逐渐增加，年平均增速0.20 g∙kg−1；从2015年，15~20 g∙kg−1等级部分明显增加，到2020年，占比达到19.04%。各年份的变异系数均在10%~100%，均属于中等强度的变异。
2) 各年份的土壤有机质含量数据均具有一定的趋势效应。2005年、2010年和2015年土壤有机质含量的块金比值均小于25%，具有强空间自相关性，其分布主要受地形、气候、土壤类型、植被等结构性因素的影响。2020年的土壤有机质块金比值为31%，属于中等强度的空间自相关性，受结构性因素和随机性因素的共同影响。
3) 2005年、2015年、2020年的土壤有机质含量最佳变异函数模型均为指数模型，2010年土壤有机质含量的最佳变异函数模型为高斯模型，其中2015年的决定系数(R2)最高，为0.953，拟合效果最好。
4) 利用普通Kriging插值法的结果显示，从空间分布来看，乐陵市东部和西南地区土壤有机质含量较高，西北部地区有机质含量较低。从时间跨度来看，乐陵市全域土壤有机质含量逐年提高，从2005年的以10~12 g∙kg−1为主到2020年10~12 g∙kg−1和12~15 g∙kg−1两等级所占比例基本一致。
基金项目
山东省自然科学基金(ZR2019MD014)；国家农业科技战略智库联盟2021年度重点研究选题(耕地保护重大科技战略问题琛研究)。
NOTES
*通讯作者。