杞麓湖盆地土壤有机质多时相空间分布与演变*

秦发侣，郭春平，陈 红，柳德江†，田育天，谢新乔，李湘伟，王世航

（1. 玉溪师范学院地理与国土工程学院，云南玉溪 653100；2. 土壤与农业可持续发展国家重点实验室（中国科学院南京土壤研究所），南京 210008；3. 玉溪市土壤与肥料工作站，云南玉溪 653100；4. 红塔烟草（集团）有限责任公司，云南玉溪 653100；5. 安徽理工大学测绘学院，安徽淮南 232001）

土壤有机质（Soil organic matter，SOM）含量直接影响耕地质量和功能，同时影响陆地碳库的丰减[1-2]。准确量化区域SOM的空间分布以及时间演变状况不仅能为农田土壤质量演变监测提供翔实数据，也能为识别土壤固碳的有效作物种植模式（如作物的种类、不同品种的种植分布等）和耕作管理措施（如施肥、灌溉等）提供基础数据，从而为区域土壤资源可持续利用、农产品增收和减缓全球气候变化做出积极贡献[3-7]。

众多学者已针对国家、省级、县级和流域尺度的SOM的时空演变与分布做了研究。如杨帆等[8]对比2005—2014年的测土配方施肥项目数据和全国第二次土壤普查数据，量化了近30年中国农田耕层SOM的变化情况。Dai等[3]根据2000年代和20世纪80年代的采样数据，对比了不同的土壤类型、土地利用方式以及自然地理单元的表层SOM的变化情况。赵明松等[9]基于全国第二次土壤普查资料和2006年采样数据，使用地统计方法进行SOM的空间插值，并通过差减2006年和20世纪80年代的SOM空间分布数据，分析了江苏省多类型土壤表层（0～20 cm）SOM的时空变异特征 。张忠启等[10]结合中国科学院耕地土壤质量分等定级与生产潜力评估项目中2007年数据和全国第二次土壤普查数据，明晰了江西省余江县1982—2007年SOC含量的空间分布和时间演变特征。

然而，以上研究均利用时间跨度大于10年的两期数据进行差减，仅以较粗时间分辨率反映SOM的变化率，难以细致且真实地反映SOM年际变化。这主要归因于区域土壤调查成本较高而使历史数据积累不足。因此，有必要在数据较翔实区域基于多期细粒度土壤调查数据量化区域SOM年际变化。

通海县是首批入选“云南省高原特色农业示范县”的县城之一。其农业耕作模式对于全省乃至我国西南高原区现代农业的发展具有典型示范作用。而杞麓湖盆地是通海县的主要农业耕作区，详尽表征该区域内SOM的空间分布与时间演变状况能为区域农业的信息化、精细化管理提供翔实的土壤基础数据。

本研究基于杞麓湖盆地2008年、2011年、2013年和2015年的密集土壤采样数据，应用普通克里格（Ordinary kriging，OK）预测该区域四个时相耕作表层的SOM空间分布，量化区域的SOM年际空间变化情况；并结合统计年鉴中各自然村的蔬菜种植总面积、粮食种植总面积、化肥（氮、磷、钾、复合肥）施用量数据，使用多元回归（Multiple regression）和探索性回归（Exploratory regression）方法分析了耕作表层SOM变化的影响因素。研究结果可为监测区域土壤质量演变提供翔实基础底图，也有助于明晰区域SOM变化影响因素，最终为区域土壤资源可持续利用决策提供支持。

1 材料与方法

1.1 研究区概况

杞麓湖盆地位于云南省通海县（102°34′～102°52′E，24°6′～24°12′ N），属亚热带半湿润高原季风气候区。区内四季温差小、早晚温差大、干湿季分明、雨热同季，多年平均降雨量900 mm，多年平均气温16℃，年日照时间2 274 h。地势由西南向东北逐渐降低，坡度3°～10°，海拔1 796～1 820 m。盆地内主要作物为蔬菜、粮食、油料、烤烟和花卉。全国第二次土壤普查划分的土壤类型为水稻土，水稻土亚类为潴育型、潜育型、淹育型以及沼泽型[11]。

1.2 样品采集与分析

为调查土壤属性的空间分布，玉溪市土壤肥料工作站分别于2008年、2011年、2013年和2015年随机采集了杞麓湖盆地耕作层（0～20 cm）679、804、484和362个土壤样品（图1）。并采用重铬酸钾（K2Cr2O7）氧化—滴定法测定了SOM含量[12-13]，提取了各样点的土壤亚类（潴育型、潜育型、淹育型、沼泽型）、土壤质地（砂壤、壤土、黏壤、黏土）、农业设施等级（不配套、基本配套、完全配套）以及灌溉程度（无灌、有灌溉潜力、能灌、保灌）信息。

1.3 离群值处理

离群值的存在会影响样本的概率分布，从而影响根据样本数据估计的半方差图[14]。且回归分析和普通克里格插值均在数据服从近似正态分布时能获得更准确的结果。因此采用Hample方法[15]根据样本的中值（Median）和中值绝对偏差（Median Absolute Deviation，MAD）来识别并剔除四个时相SOM的离群值（outliers）。基本原则是当样本中的数值大于（或小于）中值加上（或减去）3倍中值绝对偏差时，则视为离群值。剔除离群值后分别剩余667、788、469和359个土壤样品。

图1 杞麓湖盆地四个时相的土壤采样点分布情况Fig. 1 Distribution of sampling sites in the Qiluhu basin relative to time phase

1.4 半方差模型拟合与克里格插值

为了将尽可能多的土壤样点用于半方差模拟的拟合，从而获得稳定的半方差模型，本研究通过留一交叉验证（Leave One Out Cross Validation，LOOCV）比较土壤研究中常用的指数和椭球模型。选取的量化指标包括误差均值（Mean errors，ME）、误差平方的均值（Mean squared errors，MSE）、克里格插值方差均值（Mean kriging variances，以及预测值与观测值离差率平方的中值（Median squared deviation ratio，MSDR）。

1.5 统计年鉴数据收集与整理

从《通海县统计年鉴》中收集并整理了杞麓湖盆地内共47个自然村的2008年、2011年、2013年和2015年的蔬菜种植总面积、粮食种植总面积、氮肥施用量、磷肥施用量、钾肥施用量、复合肥施用量数据，并将这些数据与村级行政区划图连接，形成空间属性一体的矢量文件。统计每个自然村2008—2015年的蔬菜种植总面积平均变化量、粮食种植总面积平均变化量、化肥（氮、磷、钾和复合肥）施用量平均变化量。

1.6 土壤有机质影响因素分析

基于2008年、2011年、2013年和2015年四期土壤样点数据，以SOM为因变量，将土壤亚类、土壤质地、农业基础设施和灌溉程度作为自变量并按各因子的等级赋予哑变量（dummy variables），进行多元回归分析。此外，基于四期SOM插值数据，以自然村为单元统计四个年份的SOM均值，并计算平均变化量。以四期SOM均值平均变化量结合各村的蔬菜总产量平均变化量、粮食总产量平均变化量和化肥（氮、磷、钾和复合肥）施用平均变化量，进行探索性回归。该方法通过普通最小二乘法（Ordinary Least Squares，OLS）评估所有待选因子组合，并依据调整R2（AdjustR2）、系数显著性P（CoefficientPvalue）、方差膨胀因子（Variance Inflation Factor，VIF）和雅克贝拉显著性P（Jarque-BeraPvalue，JB）来选择较为合适的回归模型。

1.7 数据处理与分析

数据的探索性分析和多元回归分析使用R软件[16]的基础包完成，土壤数据的半方差模型拟合及普通克里格插值使用R软件的gstat软件包[17]完成，探索性回归和部分地图的制作在ArcMap 10.5中完成。

2 结 果

2.1 SOM基本统计特征

2008年、2011年、2013年和2015年的原始样本的SOM偏度和峰度均较大，概率分布为右向长尾。虽然剔除离群值后剩余的土壤样本（667、788、469和359）的概率分布并不完全服从正态分布，但四期样本概率分布的偏度均小于0.5，而峰度均明显降低，较处理前更接近于正态分布。用剔除离群值后的样本的变异系数CV来表征数据的变异性，可见四期SOM的变异性均大于35%。根据Wilding[18]推荐的分类指标，当CV＞35%时为强度变异性，CV＞15%为中等变异性，CV＜15%为轻微变异。

统计对比剔除离群值后四个时相的样本，结果表明相对2008年，2011年的SOM均值显著降低；而相对2011年，2013年和2015年SOM均值显著增加。不同年份SOM均值差异的95%置信区间最大范围为均值±3.13 g·kg-1。此值将作为下文中统计SOM变化面积的临界参考值。

2.2 SOM的空间结构特征

数量不小于100且在研究区内分布均匀的采样点可以较为稳定且真实地估算区域变量的半方差函数从而反映其空间相关关系[19-20]。根据剔除离群值后样本拟合四个时相SOM的半方差的指数型和椭球模型[21]的块金值（Nugget）、偏基台（Partial sill）、变程值（Range）以及块金-基台比（表1）。用ME、MSE、MSDR和R2比较指数和椭球模型的适用性。其中ME应接近于0，而MSE则应越小越好；衡量模型好坏的最主要指标为MSDR和R2，其值越接近于1，对应模型拟合的越好。结果表明指数模型较椭球模型能更好地拟合2008年（MSDR和R2分别为0.988和0.986）和2011年（MSDR和R2分别为0.956和0.947）SOM的半方差函数，而椭球模型在拟合2013年（MSDR和R2分别为1.065和0.973）和2015年（MSDR和R2分别为1.032和0.988）SOM时略好于指数模型。

表1 根据杞麓湖盆地四个时相土壤有机质调查值拟合的指数和椭球半方差模型的基本参数Table 1 Basic parameters of the exponential and spherical semi-variance models inferred from the soil survey on soil organic matter in the Qiluhu Basin during the four time phases

2.3 SOM时空分布特征

本研究采用普通克里格方法，基于指数模型预测2008年和2011年、椭球模型预测2013年和2015年SOM的空间分布，并输出粒度为50 m×50 m的SOM分布图（图2）。可以看出，杞麓湖盆地SOM的空间分布特征在不同年份总体格局相似，但局部差异明显。总体上四个时相的SOM从西南向东北呈现出低-高-低的分布特征。局部来看，盆地西北角的SOM从2008年至2011年有较小程度的降低，而在2013年和2015年均呈现出明显的增加；中部高值核心区位置基本一致，但高值区的辐射范围从2008年至2011年表现为收缩，而2011年至2013年和2015年则表现为逐渐扩张。

为了明晰不同年份间SOM的空间变化情况，以2008年的SOM数据为基底，比较2011年、2013年和2015年的SOM空间变化（图3）。相比2008年，SOM增加的区域面积有逐年扩大的趋势，增加的区域主要集中在盆地的西北角、西南角、南部以及杞麓湖的西北岸和南岸大部。据插值结果统计，杞麓湖盆地2011年、2013年和2015年相对于2008年的SOM的变化值的范围为-21.87～17.65 g·kg-1。如依据上文中不同年份SOM均值差异的95%置信区间将0±3.13 g·kg-1确定为SOM减少或增加的临界区，并将SOM划分为＜-10 g·kg-1、-10～-3 g·kg-1、-3～3 g·kg-1、3～10 g·kg-1、＞10 g·kg-1的L1～L5五个等级，则各等级面积的变化情况明显（图4）。从2008年至2015年的三个时间段内SOM连续增加的区域面积约占68.7%，为SOM连续减少的区域面积的两倍多。

2.4 耕层土壤有机碳的影响因素

图2 土壤有机质的50 m×50 m预测值分布图Fig. 2 Predicted SOM distribution map with 50 m×50 m resolution

图3 四个时期土壤有机质变化分布图Fig. 3 Soil organic matter distribution map relative to time phase

图4 杞麓湖盆地五个土壤有机碳变化等级的面积Fig. 4 Areas of the five levels of soil organic matter variation in the Qiluhu Basin

将土壤亚类、土壤质地、农业基础设施和灌溉程度作为自变量并赋予哑变量后进行回归分析。结果表明，除灌溉程度对土壤有机质的贡献不显著外，土壤亚类、土壤质地和农业基础设施这三个因子均呈现出显著的影响（P＜0.001）。其中，对杞麓湖盆地耕层SOM影响最大的是土壤亚类的类型，一共能解释14.3%；而农业设施和土壤质地分别能解释2.6%和1.3%。在水稻土亚类中，潴育型、潜育型和沼泽型分别较淹育型高出10.94±0.74、21.12±1.59和25.36±1.37 g·kg-1。相对于没有配套农业设施的农田，基本配套和完全配套的农田耕层SOM均有所提高，但前者提高的程度未达到统计显著性。相对于砂壤，SOM在壤土和黏壤中的含量显著更高（P＜0.001），但黏土并未呈现出与砂壤显著差异。

探索性回归分析结果表明，在不同的候选因子组成的模型中，由粮食作物总产量平均变化量、蔬菜总产量平均变化量、氮肥施用量平均变化量和复合肥施用量平均变化量四个因子组成的模型的调整R2最大（0.23），其中前两个因子达到显著性（P＜0.001），而后两个因子未达到显著性。四个因子的VIF（2.95）并未超过阈值7.5，因子间的共线性不明显。但JB值（P=0.000）表明该模型拟合后的残差不符合正态性分布，说明仍有除了以上四个因子之外的其他重要因子影响村级的SOM的分布。

3 讨 论

若按全国第二次土壤普查时设定的标准，将耕层SOM按大于40 g·kg-1、30～40 g·kg-1、20～30 g·kg-1、10～20 g·kg-1、6～10 g·kg-1、小于等于6 g·kg-1标准划分为六级[8]。那么，杞麓湖盆地耕层SOM均值除2011年处于第二级外，其余三个时相均为第一级，说明整个盆地的耕层SOM均值处于较高的水平。与杨帆等[8]研究结果比较，杞麓湖盆地的SOM显著高于全国耕层SOM平均水平（24.65 g·kg-1），接近我国耕层SOM最高的黑龙江省的平均水平（40.43 g·kg-1），显著高于云南省耕层SOM平均水平（35.37 g·kg-1）。

从影响机制来看，粮食作物总产量平均变化量、蔬菜总产量平均变化量和氮肥施用量平均变化量对杞麓湖盆地耕层SOM的影响为负作用，而复合肥施用量平均变化量的影响为正作用。在杞麓湖盆地内粮食和蔬菜产量持续增加的自然村，由于蔬菜的高复种指数（多数为四茬）、粮食作物的高秸秆移除率（除根系外基本全部移除）以及高度依赖化肥投入而轻视有机肥的施用导致了SOM的降低；而调整化肥施用中氮肥和复合肥的占比情况，适当减少氮肥占比并增加复合肥的占比在一定程度上有利于维持甚至提升SOM含量。

根据已有研究，为了进一步提高杞麓湖盆地的耕层SOM，可以尝试有机肥与化肥的科学配比施用[22]；在地理环境类似的洱海流域的研究结果表明，降低化肥施用量，同时增加有机肥的施用量能确实提高蔬菜地的耕层SOM[23]。此外，冬季休耕或少耕也被证明能有效地提高滇中地区耕层SOM总量[24]。当然，科学地选择耕作管理措施需要综合考虑生态效益、农民的经济效益及由此带来的社会效益。因此，在未来的研究工作中应加强长期定位对比实验的研究，以此来共同促进杞麓湖盆地耕层SOM的科学管理。

4 结 论

借助普通克里格插值方法能获得杞麓湖盆地2008年、2011年、2013年和2015年四个时相SOM的空间分布图。指数和椭球半方差模型在空间插值的精度评价方面差别较小，前者更适用于2008年和2011年的SOM数据，而后者更适用于2013年和2015年的数据。整个盆地的SOM分布从西南到中部再到东北部呈现出低-高-低的分布特征，并且在时间系列上SOM高值核心区的位置并未呈现显著变化，但其辐射范围有逐年扩大的趋势。从2008年到2015年的三个时间段内，SOM呈现出连续增加的区域面积为SOM连续减小的区域面积的两倍多。据土壤样点调查数据分析，土壤亚类类型对整个区域内SOM的变异起主要影响作用，可解释变异的14.3%，而农业设施和土壤质地分别能解释2.6%和1.3%。以自然村为基础单元，SOM的降低主要由粮食和蔬菜的高产出所导致，而适当的减少氮肥施用量并增加复合肥施用量能在一定程度上维持甚至提升SOM含量。