孫 耿，孫 梅，周峻宇，黃 晶，曾躍輝，羅尊長*

（1.湖南土壤肥料研究所，湖南 長沙 410125；2.中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所/ 耕地培育技術(shù)國家工程實驗室，北京 100081；3.湖南省農(nóng)業(yè)科學(xué)院，湖南 長沙 410125）

作為全國聞名的“魚米之鄉(xiāng)”和優(yōu)質(zhì)米生產(chǎn)基地，寧鄉(xiāng)市糧食產(chǎn)量保持湖南省第一，連續(xù)七年獲評全國糧食生產(chǎn)先進(jìn)縣。寧鄉(xiāng)市主要耕作土壤是水稻土，以雙季稻種植為主，是典型的雙季稻產(chǎn)區(qū)。2017年全年糧食種植面積13.05萬hm2，其中水稻播種面積12.09萬hm2，糧食總產(chǎn)量81.8萬t［1］。持續(xù)耕作及水稻產(chǎn)量的不斷增加，對稻田土壤肥力產(chǎn)生深遠(yuǎn)影響，因此，明確40年來稻田土壤肥力時空變異特征，對寧鄉(xiāng)市稻田土壤培肥管理和水稻生產(chǎn)等具有重要的現(xiàn)實意義。

土壤養(yǎng)分具有較強(qiáng)的時空變異特性［2-5］。針對土壤養(yǎng)分時空變異的研究，時間上，對不同時間點土壤養(yǎng)分進(jìn)行統(tǒng)計分析并比較其變化趨勢，一般以1980年左右的全國第二次土壤普查和2005～2014年的測土配方施肥數(shù)據(jù)作為參考［6-7］，或?qū)﹂L期定位監(jiān)測點進(jìn)行跟蹤比較［8-9］，也可以采用大面積土壤養(yǎng)分定位監(jiān)測和肥料定位試驗相結(jié)合的辦法［10］?？臻g上，主要應(yīng)用地統(tǒng)計學(xué)與GIS相結(jié)合的方法，采用泛克里格、簡單克里格、普通克里格和張力樣條函數(shù)等方法進(jìn)行插值［11-12］，并結(jié)合普通克里格法、回歸克里格法、隨機(jī)森林模型等進(jìn)行空間預(yù) 測［13-14］，研究相應(yīng)區(qū)域土壤養(yǎng)分的空間變異特征，并進(jìn)行影響因素分析。一般認(rèn)為，結(jié)構(gòu)性因素（地形地貌、母質(zhì)、氣候等）和人為因素（土地利用方式、耕作、施肥、管理水平等）是土壤養(yǎng)分空間分布的影響因子［15-19］，但也有研究將養(yǎng)分空間格局的演變歸因于作物殘茬、根系生物量的增加和秸稈還田［20］。

目前雖然有針對寧鄉(xiāng)市耕地質(zhì)量的研究［21］，但關(guān)于寧鄉(xiāng)市稻田土壤養(yǎng)分的時空變異及其影響因素未見報道，本研究基于6項常規(guī)指標(biāo)（pH、有機(jī)質(zhì)、全氮、堿解氮、有效磷、速效鉀）和土壤肥力綜合肥力指標(biāo)，通過克里格法進(jìn)行空間插值與模型預(yù)測，系統(tǒng)分析了40年來寧鄉(xiāng)市稻田土壤肥力時空變異特征，并對其成因進(jìn)行深入探討。

1 材料與方法

1.1 研究區(qū)概況

寧鄉(xiāng)市地處湘東偏北的洞庭湖南緣地區(qū)（27°55′～28°29′N，111°53′～112°46′E），屬中亞熱帶向北亞熱帶過渡的大陸性季風(fēng)濕潤氣候，四季分明，寒冷期短，炎熱期長。年日平均氣溫16.8℃，年平均無霜期274 d，年平均日照1 737.6 h，境內(nèi)雨水充足，年均降水量1 358.3 mm，相對濕度高達(dá)81%。境內(nèi)多為丘陵地帶，西部的溈山區(qū)域是雪峰山龐大東部地帶的南側(cè)主干區(qū)，往東則是雪峰山余脈向東北濱湖平原過渡地帶，境內(nèi)地貌有山地、丘崗、平原。地表輪廓大體是北、西、南緣山地環(huán)繞，自西向東呈階梯狀逐級傾斜，東南丘陵起伏，北部崗地平緩，東北低平開闊。水資源豐富，境內(nèi)有溈水、烏江、楚江、靳江四條主要河流。成土母質(zhì)主要是花崗巖、板頁巖、砂巖、石灰?guī)r、紫色巖、第四紀(jì)紅土及河流沉積物［21］。

1.2 數(shù)據(jù)來源

本研究以1979年寧鄉(xiāng)第二次土壤普查為基礎(chǔ)，當(dāng)時取樣缺乏GPS定位，僅有文字記載，但在40年的土地流轉(zhuǎn)生產(chǎn)過程中，田塊名基本保持不變。而2005～2007年的測土配方施肥有完整的田塊信息，因此對1979年和2007年的相同田塊進(jìn)行篩選，共獲取139組數(shù)據(jù)。2018年則在歷史點的基礎(chǔ)上，隨機(jī)采集了133個表層土壤（0～20 cm）樣品，并調(diào)查了取樣點的水稻產(chǎn)量。所有樣品經(jīng)風(fēng)干、剔除雜質(zhì)、研磨、分別過0.85和0.15 mm篩后，用于土壤理化性質(zhì)分析。所有測定指標(biāo)均采用常規(guī)方法進(jìn)行分析［22］。

為使兩個時間節(jié)點的土壤肥力指標(biāo)相對應(yīng)，特選取土壤pH、有機(jī)質(zhì)、全氮、堿解氮、有效磷、速效鉀6項指標(biāo)及其綜合肥力指數(shù)進(jìn)行土壤肥力時空變異分析。

圖1 寧鄉(xiāng)市稻田土壤采樣點位分布

1.3 土壤綜合肥力指數(shù)評價方法

土壤綜合肥力指數(shù)采用Fuzzy法進(jìn)行評價。首先建立評價指標(biāo)的隸屬度函數(shù)模型，土壤有機(jī)質(zhì)、全氮、堿解氮、有效磷、速效鉀等均采用S型隸屬函數(shù)，pH采用拋物線型隸屬函數(shù)［23］。

S型隸屬函數(shù)表達(dá)式如方程（1）：

拋物線型隸屬函數(shù)表達(dá)式如方程（2）：

根據(jù)前人研究結(jié)果［24］以及結(jié)合稻田土壤的實際，本研究確定各指標(biāo)的隸屬度函數(shù)曲線中轉(zhuǎn)折點的相應(yīng)取值如表1所示：

表1 S型隸屬度函數(shù)轉(zhuǎn)折點取值

各個指標(biāo)的單因子指數(shù)相加，即為土壤綜合肥力指數(shù)（IFI）。

式中，Ni為第i項肥力指標(biāo)的隸屬度值，Wi為第i項指標(biāo)的權(quán)重系數(shù)。IFI取值為0～1，本研究將土壤綜合肥力指數(shù)分為5個等級，分別是高（IFI≥ 0.8）、較高（0.8>IFI≥0.6）、中等（0.6>IFI≥0.4）、較低（0.4>IFI≥0.2）及低（IFI<0.2）。

1.4 數(shù)據(jù)處理

采用SPSS 22.0軟件進(jìn)行數(shù)據(jù)統(tǒng)計分析；土壤肥力指標(biāo)的半方差函數(shù)模型通過GS+7.0實現(xiàn)；土壤肥力指標(biāo)的空間分析則通過ArcGIS 10.2相關(guān)分析模塊完成。利用R語言（3.4.4）進(jìn)行數(shù)據(jù)的隨機(jī)森林重要性分析。

2 結(jié)果與分析

2.1 土壤肥力描述性統(tǒng)計

由表2可知，2018年，土壤pH平均值為5.69，較1979年降低0.58；土壤有機(jī)質(zhì)、全氮、堿解氮、有效磷和速效鉀平均值分別為38.71、2.17、187.9、9.20及95.46 mg/kg，分別較1979年增加了13.1%、98.5%、20.1%、98.2%和13.7%，以全氮和有效磷增加最為明顯。2018年土壤綜合肥力指數(shù)為0.78，較1979年增加0.25，增幅達(dá)47.9%。由此可見，經(jīng)過40年的發(fā)展，研究區(qū)稻田土壤肥力增加明顯，但同時也呈現(xiàn)明顯酸化趨勢。

2018年，土壤pH的變異系數(shù)為9.37%，屬于弱變異程度。1979年，土壤有效磷變異系數(shù)為106.5%，屬于強(qiáng)變異程度，其他指標(biāo)均介于10%～100%之間，為中等變異程度。整體上，相對于1979年，2018年土壤肥力指標(biāo)的變異系數(shù)均有所減小，表明變異程度有所降低。

表2 土壤肥力統(tǒng)計性描述

2.2 土壤肥力指標(biāo)空間結(jié)構(gòu)分析

利用ArcGIS地統(tǒng)計模塊對各肥力指標(biāo)半方差模型進(jìn)行擬合，并進(jìn)行預(yù)測誤差比較，比較標(biāo)準(zhǔn)是：標(biāo)準(zhǔn)平均值（MSE）最接近于0；標(biāo)準(zhǔn)均方根預(yù)測誤差（RMSSE）最接近于1［25］。由表3可以看出，就預(yù)測誤差而言，1979年和2018年有效磷的RMSSE分別為0.747和0.893，模型擬合性較差，其他各肥力指標(biāo)的半方差函數(shù)擬合較好，理論模型能較好反映各指標(biāo)的空間結(jié)構(gòu) 特征。

塊金比值［C0/（C0+C）］可表示空間變異程度（隨機(jī)因素引起的空間變異占系統(tǒng)總變異的比例），若該比值較高，說明由隨機(jī)部分引起的空間變異性程度較大；反之，則由空間自相關(guān)部分引起的空間變異性程度較大；若該比值接近于1，則說明該變量在整個尺度上具有恒定的變異。從結(jié)構(gòu)性因素角度看，C0/（C0+C）可表示系統(tǒng)變量的空間相關(guān)性程度，如果比例<25%，說明變量具有強(qiáng)烈的空間相關(guān)性；若比例在25%～75%之間，變量具有中等空間相關(guān)性；比例>75%時，變量空間相關(guān)性 很弱［26］。

由表3可以看出，1979年和2018年土壤有效磷的C0/（C0+C）均為100%，空間結(jié)構(gòu)表現(xiàn)為純塊金效應(yīng)，在研究區(qū)域的取樣尺寸下，有效磷變異恒定，土壤有效磷不能建立有效半方差函數(shù)模型。兩個時間點的土壤pH、有機(jī)質(zhì)、全氮、堿解氮、速效鉀和綜合肥力的C0/（C0+C）則均<25%，表現(xiàn)較強(qiáng)的空間相關(guān)性。

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

2.3.1 稻田土壤肥力時空變異

1979年，有機(jī)質(zhì)西北部地區(qū)最高為41.8 g/kg，東北部地區(qū)最低為28.1 g/kg；全氮中南部地區(qū)最高1.42 g/kg，東南部地區(qū)最低為0.81 g/kg；堿解氮東北部地區(qū)最高為194.7 mg/kg，南部地區(qū)最低為105.2 mg/kg；有效磷區(qū)域差異較小，位于1.17～8.64 mg/kg之間；速效鉀中南部地區(qū)最大為119.0 mg/kg，西南地區(qū)最小為51.6 mg/kg；pH東南部地區(qū)最高達(dá)7.09，西北部地區(qū)最低為5.59。

表3 土壤肥力指標(biāo)半方差函數(shù)模型及其參數(shù)

2018年，有機(jī)質(zhì)中部地區(qū)最高達(dá)44.8 g/kg，中南部地區(qū)最低為32.6 g/kg；全氮中部地區(qū)最高達(dá)2.52 g/kg，東北部地區(qū)最低為1.87 g/kg；堿解氮東南部地區(qū)最高為239.7 mg/kg，東北部地區(qū)最低為150.2 mg/kg；有效磷中北部地區(qū)最高為14.42 mg/kg，中南部地區(qū)最低為4.05 mg/kg；速效鉀中北部地區(qū)最高為150.1 mg/kg，西部地區(qū)最低為62.8 mg/kg；pH北部地區(qū)最高為6.81，中部地區(qū)最低為5.09。

40年來，有機(jī)質(zhì)北部地區(qū)最高增加11.3 g/kg，南部地區(qū)降低3.5 g/kg；全氮均有增加，其中東南部地區(qū)增加最多達(dá)1.58 g/kg，東北部地區(qū)增加最小為0.67 g/kg；堿解氮除東北部地區(qū)有所降低，其他地區(qū)最高增加達(dá)81.5 mg/kg；有效磷中北部地區(qū)最高增加13.13 mg/kg，東北部地區(qū)降低1.94 mg/kg；速效鉀中北部地區(qū)增加最高達(dá)72.9 mg/kg，中南部地區(qū)則有一定幅度的降低；pH北部地區(qū)略有增加，東南部地區(qū)降低最多達(dá)1.18。

2.3.2 稻田土壤綜合肥力指數(shù)時空變異

從圖3可以看出，1979年研究區(qū)大部分地區(qū)稻田土壤綜合肥力處于中等水平（0.41

2.4 環(huán)境因素與稻田土壤養(yǎng)分空間分布相關(guān)性分析

通過下載寧鄉(xiāng)市30 m×30 m精度的DEM數(shù)據(jù)（地理空間數(shù)據(jù)云http：//www.gscloud.cn/），采用ArcGIS軟件提取各個樣本的地形因子數(shù)據(jù)，通過SPSS軟件分析，得到研究區(qū)內(nèi)主要稻田土壤養(yǎng)分與海拔、坡度、坡向、曲率等多種因子的相關(guān)性。從表4可以看出，1979年，研究區(qū)稻田土壤速效鉀和坡度與海拔呈顯著負(fù)相關(guān)，綜合肥力指數(shù)與有機(jī)質(zhì)、全氮和堿解氮顯著正相關(guān)，2018年，研究區(qū)稻田土壤速效鉀和坡度呈顯著負(fù)相關(guān)，綜合肥力指數(shù)與曲率、有機(jī)質(zhì)、全氮、堿解氮和速效鉀呈顯著正相關(guān)，其他指標(biāo)與地形因子無顯著相關(guān)性。

2.5 稻田土壤肥力指標(biāo)重要性分析

由于兩個時間點的稻田土壤綜合肥力差異較大，因此，利用隨機(jī)森林對土壤綜合肥力的影響因素進(jìn)行重要性分析（圖4）。平均下降精度是指預(yù)測誤差準(zhǔn)確性降低的程度，該值越大表示該變量的重要性越大［27］。1979年，堿解氮對土壤綜合肥力指數(shù)的重要性得分達(dá)52.7%，其次為有機(jī)質(zhì)，重要性得分為30.3%，全氮僅為16.5%，故堿解氮和有機(jī)質(zhì)是土壤綜合肥力的主要驅(qū)動因素。2018年，全氮和有機(jī)質(zhì)的重要性得分分別為33.1%和32.4%，土壤堿解氮、有效磷、速效鉀的重要性得分在10.5%～12.3%之間，故全氮和有機(jī)質(zhì)是土壤綜合肥力的主要驅(qū)動因素。

圖2 土壤肥力指標(biāo)時空變化

圖3 土壤綜合肥力指數(shù)時空變化

表4 水稻土肥力與地形因子相關(guān)性分析

圖4 土壤肥力指標(biāo)重要性

2.6 水稻相對產(chǎn)量與土壤綜合肥力相關(guān)性分析

由圖5可以看出，當(dāng)前肥力水平下，研究區(qū)水稻相對產(chǎn)量隨著土壤綜合肥力指數(shù)的升高，呈增加趨勢，通過線性擬合，發(fā)現(xiàn)水稻相對產(chǎn)量與土壤綜合肥力指數(shù)呈極顯著正相關(guān)（P<0.01）。

圖5 土壤綜合肥力指數(shù)與相對產(chǎn)量相關(guān)性

3 討論

研究區(qū)40年間土壤主要肥力指標(biāo)表現(xiàn)出不同程度的增加，以全氮和有效磷增加最為明顯，增幅分別達(dá)98.5%和98.2%，這與前期研究結(jié)果一致［22］。產(chǎn)生這一變化趨勢的原因，一方面是肥料的長期投入，增加了土壤養(yǎng)分含量［28］，據(jù)估算，從1983年到2010年，湖南省化肥年均施用量從1.02×106t 增加到2.36×106t，年均增長率為3.12%，單位面積年均增長率為4.20%［29］。另一方面隨著生產(chǎn)水平 的不斷提高，水稻單產(chǎn)不斷增加［30］，作物殘茬和 根系生物量的增加，并且秸稈還田，對提高土壤養(yǎng)分庫容具有更為重要的影響［31-32］。但在土壤養(yǎng)分提高的同時，土壤pH降低0.52，有酸化趨勢，這 一變化趨勢可能成為稻田土壤肥力和生產(chǎn)力合理演變的制約因素［33］。土壤酸化程度取決于作物產(chǎn)量、施氮量與降水的相互作用及土壤類型和母質(zhì)的緩沖 能力［34］，因此，要合理調(diào)控氮肥的施用方法和施用量，并注重平衡施肥。

基于主要肥力指標(biāo)的增加，研究區(qū)稻田土壤綜合肥力指數(shù)提高達(dá)47.9%。從因子重要性分析，40年來，有機(jī)質(zhì)、全氮和堿解氮一直都是土壤綜合肥力指數(shù)的重要因素，但堿解氮對土壤綜合肥力指數(shù)的重要性得分從52.7%降至12.3%，而全氮則從16.5%提高至33.1%（圖4），有機(jī)質(zhì)則保持在30%左右，可見，全氮的增加是稻田土壤綜合肥力指數(shù)提升的關(guān)鍵驅(qū)動因素。進(jìn)一步相關(guān)性分析表明，水稻相對產(chǎn)量與土壤綜合肥力指數(shù)呈極顯著正相關(guān)。因此，土壤氮素和有機(jī)質(zhì)的維持是稻田土壤肥力和水稻生產(chǎn)的重要措施［35］，稻田生產(chǎn)在合理調(diào)控氮素的前提下，更要提倡有機(jī)無機(jī)配合施用［36］。

水稻土隸屬于人為土綱，在長期淹水耕種條件下，受到人為活動和自然成土的雙重因素影 響［37］。本研究中，土壤有效磷40年來塊金比值均為100%，空間變異恒定，空間相關(guān)性較弱，表明其受隨機(jī)因素影響較大，土壤pH、有機(jī)質(zhì)、全氮、堿解氮、速效鉀和綜合肥力的塊金比值則<25%，表現(xiàn)較強(qiáng)的空間相關(guān)性。環(huán)境因素空間相關(guān)性分析發(fā)現(xiàn)，1979年稻田土壤速效鉀和坡度與海拔呈顯著負(fù)相關(guān)，2018年稻田土壤速效鉀和坡向呈顯著負(fù)相關(guān)，肥力指數(shù)與曲率呈顯著正相關(guān)，表明本區(qū)域土壤肥力空間異質(zhì)性受結(jié)構(gòu)性因素（地形、土壤類型等）和人為因素（施肥、耕作措施、種植制度等）共同影響［19，38］。

4 結(jié)論

40年間，研究區(qū)稻田土壤有機(jī)質(zhì)、養(yǎng)分及綜合肥力指數(shù)均有不同程度的增加，但pH降低明顯，呈逐漸酸化趨勢。

土壤有效磷空間相關(guān)性較弱，受隨機(jī)因素影響較大；有機(jī)質(zhì)、全氮、堿解氮、速效鉀、pH和綜合肥力指數(shù)空間相關(guān)性較強(qiáng)，受海拔等地形因子和施肥等人為因素的共同作用。

有機(jī)質(zhì)、全氮和堿解氮是土壤綜合肥力指數(shù)的重要因素，全氮的大幅增加是土壤綜合肥力指數(shù)提升的關(guān)鍵驅(qū)動原因，當(dāng)前水稻產(chǎn)量與土壤綜合肥力指數(shù)呈極顯著正相關(guān)。