40年来宁乡市稻田土壤肥力时空变异分析

孙 耿，孙 梅，周峻宇，黄 晶，曾跃辉，罗尊长*

（1.湖南土壤肥料研究所，湖南 长沙 410125；2.中国农业科学院农业资源与农业区划研究所/ 耕地培育技术国家工程实验室，北京 100081；3.湖南省农业科学院，湖南 长沙 410125）

作为全国闻名的“鱼米之乡”和优质米生产基地，宁乡市粮食产量保持湖南省第一，连续七年获评全国粮食生产先进县。宁乡市主要耕作土壤是水稻土，以双季稻种植为主，是典型的双季稻产区。2017年全年粮食种植面积13.05万hm2，其中水稻播种面积12.09万hm2，粮食总产量81.8万t［1］。持续耕作及水稻产量的不断增加，对稻田土壤肥力产生深远影响，因此，明确40年来稻田土壤肥力时空变异特征，对宁乡市稻田土壤培肥管理和水稻生产等具有重要的现实意义。

土壤养分具有较强的时空变异特性［2-5］。针对土壤养分时空变异的研究，时间上，对不同时间点土壤养分进行统计分析并比较其变化趋势，一般以1980年左右的全国第二次土壤普查和2005～2014年的测土配方施肥数据作为参考［6-7］，或对长期定位监测点进行跟踪比较［8-9］，也可以采用大面积土壤养分定位监测和肥料定位试验相结合的办法［10］。空间上，主要应用地统计学与GIS相结合的方法，采用泛克里格、简单克里格、普通克里格和张力样条函数等方法进行插值［11-12］，并结合普通克里格法、回归克里格法、随机森林模型等进行空间预 测［13-14］，研究相应区域土壤养分的空间变异特征，并进行影响因素分析。一般认为，结构性因素（地形地貌、母质、气候等）和人为因素（土地利用方式、耕作、施肥、管理水平等）是土壤养分空间分布的影响因子［15-19］，但也有研究将养分空间格局的演变归因于作物残茬、根系生物量的增加和秸秆还田［20］。

目前虽然有针对宁乡市耕地质量的研究［21］，但关于宁乡市稻田土壤养分的时空变异及其影响因素未见报道，本研究基于6项常规指标（pH、有机质、全氮、碱解氮、有效磷、速效钾）和土壤肥力综合肥力指标，通过克里格法进行空间插值与模型预测，系统分析了40年来宁乡市稻田土壤肥力时空变异特征，并对其成因进行深入探讨。

1 材料与方法

1.1 研究区概况

宁乡市地处湘东偏北的洞庭湖南缘地区（27°55′～28°29′N，111°53′～112°46′E），属中亚热带向北亚热带过渡的大陆性季风湿润气候，四季分明，寒冷期短，炎热期长。年日平均气温16.8℃，年平均无霜期274 d，年平均日照1 737.6 h，境内雨水充足，年均降水量1 358.3 mm，相对湿度高达81%。境内多为丘陵地带，西部的沩山区域是雪峰山庞大东部地带的南侧主干区，往东则是雪峰山余脉向东北滨湖平原过渡地带，境内地貌有山地、丘岗、平原。地表轮廓大体是北、西、南缘山地环绕，自西向东呈阶梯状逐级倾斜，东南丘陵起伏，北部岗地平缓，东北低平开阔。水资源丰富，境内有沩水、乌江、楚江、靳江四条主要河流。成土母质主要是花岗岩、板页岩、砂岩、石灰岩、紫色岩、第四纪红土及河流沉积物［21］。

1.2 数据来源

本研究以1979年宁乡第二次土壤普查为基础，当时取样缺乏GPS定位，仅有文字记载，但在40年的土地流转生产过程中，田块名基本保持不变。而2005～2007年的测土配方施肥有完整的田块信息，因此对1979年和2007年的相同田块进行筛选，共获取139组数据。2018年则在历史点的基础上，随机采集了133个表层土壤（0～20 cm）样品，并调查了取样点的水稻产量。所有样品经风干、剔除杂质、研磨、分别过0.85和0.15 mm筛后，用于土壤理化性质分析。所有测定指标均采用常规方法进行分析［22］。

为使两个时间节点的土壤肥力指标相对应，特选取土壤pH、有机质、全氮、碱解氮、有效磷、速效钾6项指标及其综合肥力指数进行土壤肥力时空变异分析。

图1 宁乡市稻田土壤采样点位分布

1.3 土壤综合肥力指数评价方法

土壤综合肥力指数采用Fuzzy法进行评价。首先建立评价指标的隶属度函数模型，土壤有机质、全氮、碱解氮、有效磷、速效钾等均采用S型隶属函数，pH采用抛物线型隶属函数［23］。

S型隶属函数表达式如方程（1）：

抛物线型隶属函数表达式如方程（2）：

根据前人研究结果［24］以及结合稻田土壤的实际，本研究确定各指标的隶属度函数曲线中转折点的相应取值如表1所示：

表1 S型隶属度函数转折点取值

各个指标的单因子指数相加，即为土壤综合肥力指数（IFI）。

式中，Ni为第i项肥力指标的隶属度值，Wi为第i项指标的权重系数。IFI取值为0～1，本研究将土壤综合肥力指数分为5个等级，分别是高（IFI≥ 0.8）、较高（0.8>IFI≥0.6）、中等（0.6>IFI≥0.4）、较低（0.4>IFI≥0.2）及低（IFI<0.2）。

1.4 数据处理

采用SPSS 22.0软件进行数据统计分析；土壤肥力指标的半方差函数模型通过GS+7.0实现；土壤肥力指标的空间分析则通过ArcGIS 10.2相关分析模块完成。利用R语言（3.4.4）进行数据的随机森林重要性分析。

2 结果与分析

2.1 土壤肥力描述性统计

由表2可知，2018年，土壤pH平均值为5.69，较1979年降低0.58；土壤有机质、全氮、碱解氮、有效磷和速效钾平均值分别为38.71、2.17、187.9、9.20及95.46 mg/kg，分别较1979年增加了13.1%、98.5%、20.1%、98.2%和13.7%，以全氮和有效磷增加最为明显。2018年土壤综合肥力指数为0.78，较1979年增加0.25，增幅达47.9%。由此可见，经过40年的发展，研究区稻田土壤肥力增加明显，但同时也呈现明显酸化趋势。

2018年，土壤pH的变异系数为9.37%，属于弱变异程度。1979年，土壤有效磷变异系数为106.5%，属于强变异程度，其他指标均介于10%～100%之间，为中等变异程度。整体上，相对于1979年，2018年土壤肥力指标的变异系数均有所减小，表明变异程度有所降低。

表2 土壤肥力统计性描述

2.2 土壤肥力指标空间结构分析

利用ArcGIS地统计模块对各肥力指标半方差模型进行拟合，并进行预测误差比较，比较标准是：标准平均值（MSE）最接近于0；标准均方根预测误差（RMSSE）最接近于1［25］。由表3可以看出，就预测误差而言，1979年和2018年有效磷的RMSSE分别为0.747和0.893，模型拟合性较差，其他各肥力指标的半方差函数拟合较好，理论模型能较好反映各指标的空间结构 特征。

块金比值［C0/（C0+C）］可表示空间变异程度（随机因素引起的空间变异占系统总变异的比例），若该比值较高，说明由随机部分引起的空间变异性程度较大；反之，则由空间自相关部分引起的空间变异性程度较大；若该比值接近于1，则说明该变量在整个尺度上具有恒定的变异。从结构性因素角度看，C0/（C0+C）可表示系统变量的空间相关性程度，如果比例<25%，说明变量具有强烈的空间相关性；若比例在25%～75%之间，变量具有中等空间相关性；比例>75%时，变量空间相关性 很弱［26］。

由表3可以看出，1979年和2018年土壤有效磷的C0/（C0+C）均为100%，空间结构表现为纯块金效应，在研究区域的取样尺寸下，有效磷变异恒定，土壤有效磷不能建立有效半方差函数模型。两个时间点的土壤pH、有机质、全氮、碱解氮、速效钾和综合肥力的C0/（C0+C）则均<25%，表现较强的空间相关性。

2.3 稻田土壤肥力时空分布特征

2.3.1 稻田土壤肥力时空变异

1979年，有机质西北部地区最高为41.8 g/kg，东北部地区最低为28.1 g/kg；全氮中南部地区最高1.42 g/kg，东南部地区最低为0.81 g/kg；碱解氮东北部地区最高为194.7 mg/kg，南部地区最低为105.2 mg/kg；有效磷区域差异较小，位于1.17～8.64 mg/kg之间；速效钾中南部地区最大为119.0 mg/kg，西南地区最小为51.6 mg/kg；pH东南部地区最高达7.09，西北部地区最低为5.59。

表3 土壤肥力指标半方差函数模型及其参数

2018年，有机质中部地区最高达44.8 g/kg，中南部地区最低为32.6 g/kg；全氮中部地区最高达2.52 g/kg，东北部地区最低为1.87 g/kg；碱解氮东南部地区最高为239.7 mg/kg，东北部地区最低为150.2 mg/kg；有效磷中北部地区最高为14.42 mg/kg，中南部地区最低为4.05 mg/kg；速效钾中北部地区最高为150.1 mg/kg，西部地区最低为62.8 mg/kg；pH北部地区最高为6.81，中部地区最低为5.09。

40年来，有机质北部地区最高增加11.3 g/kg，南部地区降低3.5 g/kg；全氮均有增加，其中东南部地区增加最多达1.58 g/kg，东北部地区增加最小为0.67 g/kg；碱解氮除东北部地区有所降低，其他地区最高增加达81.5 mg/kg；有效磷中北部地区最高增加13.13 mg/kg，东北部地区降低1.94 mg/kg；速效钾中北部地区增加最高达72.9 mg/kg，中南部地区则有一定幅度的降低；pH北部地区略有增加，东南部地区降低最多达1.18。

2.3.2 稻田土壤综合肥力指数时空变异

从图3可以看出，1979年研究区大部分地区稻田土壤综合肥力处于中等水平（0.41

2.4 环境因素与稻田土壤养分空间分布相关性分析

通过下载宁乡市30 m×30 m精度的DEM数据（地理空间数据云http：//www.gscloud.cn/），采用ArcGIS软件提取各个样本的地形因子数据，通过SPSS软件分析，得到研究区内主要稻田土壤养分与海拔、坡度、坡向、曲率等多种因子的相关性。从表4可以看出，1979年，研究区稻田土壤速效钾和坡度与海拔呈显著负相关，综合肥力指数与有机质、全氮和碱解氮显著正相关，2018年，研究区稻田土壤速效钾和坡度呈显著负相关，综合肥力指数与曲率、有机质、全氮、碱解氮和速效钾呈显著正相关，其他指标与地形因子无显著相关性。

2.5 稻田土壤肥力指标重要性分析

由于两个时间点的稻田土壤综合肥力差异较大，因此，利用随机森林对土壤综合肥力的影响因素进行重要性分析（图4）。平均下降精度是指预测误差准确性降低的程度，该值越大表示该变量的重要性越大［27］。1979年，碱解氮对土壤综合肥力指数的重要性得分达52.7%，其次为有机质，重要性得分为30.3%，全氮仅为16.5%，故碱解氮和有机质是土壤综合肥力的主要驱动因素。2018年，全氮和有机质的重要性得分分别为33.1%和32.4%，土壤碱解氮、有效磷、速效钾的重要性得分在10.5%～12.3%之间，故全氮和有机质是土壤综合肥力的主要驱动因素。

图2 土壤肥力指标时空变化

图3 土壤综合肥力指数时空变化

表4 水稻土肥力与地形因子相关性分析

图4 土壤肥力指标重要性

2.6 水稻相对产量与土壤综合肥力相关性分析

由图5可以看出，当前肥力水平下，研究区水稻相对产量随着土壤综合肥力指数的升高，呈增加趋势，通过线性拟合，发现水稻相对产量与土壤综合肥力指数呈极显著正相关（P<0.01）。

图5 土壤综合肥力指数与相对产量相关性

3 讨论

研究区40年间土壤主要肥力指标表现出不同程度的增加，以全氮和有效磷增加最为明显，增幅分别达98.5%和98.2%，这与前期研究结果一致［22］。产生这一变化趋势的原因，一方面是肥料的长期投入，增加了土壤养分含量［28］，据估算，从1983年到2010年，湖南省化肥年均施用量从1.02×106t 增加到2.36×106t，年均增长率为3.12%，单位面积年均增长率为4.20%［29］。另一方面随着生产水平 的不断提高，水稻单产不断增加［30］，作物残茬和 根系生物量的增加，并且秸秆还田，对提高土壤养分库容具有更为重要的影响［31-32］。但在土壤养分提高的同时，土壤pH降低0.52，有酸化趋势，这 一变化趋势可能成为稻田土壤肥力和生产力合理演变的制约因素［33］。土壤酸化程度取决于作物产量、施氮量与降水的相互作用及土壤类型和母质的缓冲 能力［34］，因此，要合理调控氮肥的施用方法和施用量，并注重平衡施肥。

基于主要肥力指标的增加，研究区稻田土壤综合肥力指数提高达47.9%。从因子重要性分析，40年来，有机质、全氮和碱解氮一直都是土壤综合肥力指数的重要因素，但碱解氮对土壤综合肥力指数的重要性得分从52.7%降至12.3%，而全氮则从16.5%提高至33.1%（图4），有机质则保持在30%左右，可见，全氮的增加是稻田土壤综合肥力指数提升的关键驱动因素。进一步相关性分析表明，水稻相对产量与土壤综合肥力指数呈极显著正相关。因此，土壤氮素和有机质的维持是稻田土壤肥力和水稻生产的重要措施［35］，稻田生产在合理调控氮素的前提下，更要提倡有机无机配合施用［36］。

水稻土隶属于人为土纲，在长期淹水耕种条件下，受到人为活动和自然成土的双重因素影 响［37］。本研究中，土壤有效磷40年来块金比值均为100%，空间变异恒定，空间相关性较弱，表明其受随机因素影响较大，土壤pH、有机质、全氮、碱解氮、速效钾和综合肥力的块金比值则<25%，表现较强的空间相关性。环境因素空间相关性分析发现，1979年稻田土壤速效钾和坡度与海拔呈显著负相关，2018年稻田土壤速效钾和坡向呈显著负相关，肥力指数与曲率呈显著正相关，表明本区域土壤肥力空间异质性受结构性因素（地形、土壤类型等）和人为因素（施肥、耕作措施、种植制度等）共同影响［19，38］。

4 结论

40年间，研究区稻田土壤有机质、养分及综合肥力指数均有不同程度的增加，但pH降低明显，呈逐渐酸化趋势。

土壤有效磷空间相关性较弱，受随机因素影响较大；有机质、全氮、碱解氮、速效钾、pH和综合肥力指数空间相关性较强，受海拔等地形因子和施肥等人为因素的共同作用。

有机质、全氮和碱解氮是土壤综合肥力指数的重要因素，全氮的大幅增加是土壤综合肥力指数提升的关键驱动原因，当前水稻产量与土壤综合肥力指数呈极显著正相关。