秦發(fā)侶，郭春平，陳 紅，柳德江?，田育天，謝新喬，李湘?zhèn)?，王世?/p>

（1. 玉溪師范學(xué)院地理與國(guó)土工程學(xué)院，云南玉溪 653100；2. 土壤與農(nóng)業(yè)可持續(xù)發(fā)展國(guó)家重點(diǎn)實(shí)驗(yàn)室（中國(guó)科學(xué)院南京土壤研究所），南京 210008；3. 玉溪市土壤與肥料工作站，云南玉溪 653100；4. 紅塔煙草（集團(tuán)）有限責(zé)任公司，云南玉溪 653100；5. 安徽理工大學(xué)測(cè)繪學(xué)院，安徽淮南 232001）

土壤有機(jī)質(zhì)（Soil organic matter，SOM）含量直接影響耕地質(zhì)量和功能，同時(shí)影響陸地碳庫(kù)的豐減[1-2]。準(zhǔn)確量化區(qū)域SOM的空間分布以及時(shí)間演變狀況不僅能為農(nóng)田土壤質(zhì)量演變監(jiān)測(cè)提供翔實(shí)數(shù)據(jù)，也能為識(shí)別土壤固碳的有效作物種植模式（如作物的種類(lèi)、不同品種的種植分布等）和耕作管理措施（如施肥、灌溉等）提供基礎(chǔ)數(shù)據(jù)，從而為區(qū)域土壤資源可持續(xù)利用、農(nóng)產(chǎn)品增收和減緩全球氣候變化做出積極貢獻(xiàn)[3-7]。

眾多學(xué)者已針對(duì)國(guó)家、省級(jí)、縣級(jí)和流域尺度的SOM的時(shí)空演變與分布做了研究。如楊帆等[8]對(duì)比2005—2014年的測(cè)土配方施肥項(xiàng)目數(shù)據(jù)和全國(guó)第二次土壤普查數(shù)據(jù)，量化了近30年中國(guó)農(nóng)田耕層SOM的變化情況。Dai等[3]根據(jù)2000年代和20世紀(jì)80年代的采樣數(shù)據(jù)，對(duì)比了不同的土壤類(lèi)型、土地利用方式以及自然地理單元的表層SOM的變化情況。趙明松等[9]基于全國(guó)第二次土壤普查資料和2006年采樣數(shù)據(jù)，使用地統(tǒng)計(jì)方法進(jìn)行SOM的空間插值，并通過(guò)差減2006年和20世紀(jì)80年代的SOM空間分布數(shù)據(jù)，分析了江蘇省多類(lèi)型土壤表層（0～20 cm）SOM的時(shí)空變異特征 。張忠啟等[10]結(jié)合中國(guó)科學(xué)院耕地土壤質(zhì)量分等定級(jí)與生產(chǎn)潛力評(píng)估項(xiàng)目中2007年數(shù)據(jù)和全國(guó)第二次土壤普查數(shù)據(jù)，明晰了江西省余江縣1982—2007年SOC含量的空間分布和時(shí)間演變特征。

然而，以上研究均利用時(shí)間跨度大于10年的兩期數(shù)據(jù)進(jìn)行差減，僅以較粗時(shí)間分辨率反映SOM的變化率，難以細(xì)致且真實(shí)地反映SOM年際變化。這主要?dú)w因于區(qū)域土壤調(diào)查成本較高而使歷史數(shù)據(jù)積累不足。因此，有必要在數(shù)據(jù)較翔實(shí)區(qū)域基于多期細(xì)粒度土壤調(diào)查數(shù)據(jù)量化區(qū)域SOM年際變化。

通?？h是首批入選“云南省高原特色農(nóng)業(yè)示范縣”的縣城之一。其農(nóng)業(yè)耕作模式對(duì)于全省乃至我國(guó)西南高原區(qū)現(xiàn)代農(nóng)業(yè)的發(fā)展具有典型示范作用。而杞麓湖盆地是通?？h的主要農(nóng)業(yè)耕作區(qū)，詳盡表征該區(qū)域內(nèi)SOM的空間分布與時(shí)間演變狀況能為區(qū)域農(nóng)業(yè)的信息化、精細(xì)化管理提供翔實(shí)的土壤基礎(chǔ)數(shù)據(jù)。

本研究基于杞麓湖盆地2008年、2011年、2013年和2015年的密集土壤采樣數(shù)據(jù)，應(yīng)用普通克里格（Ordinary kriging，OK）預(yù)測(cè)該區(qū)域四個(gè)時(shí)相耕作表層的SOM空間分布，量化區(qū)域的SOM年際空間變化情況；并結(jié)合統(tǒng)計(jì)年鑒中各自然村的蔬菜種植總面積、糧食種植總面積、化肥（氮、磷、鉀、復(fù)合肥）施用量數(shù)據(jù)，使用多元回歸（Multiple regression）和探索性回歸（Exploratory regression）方法分析了耕作表層SOM變化的影響因素。研究結(jié)果可為監(jiān)測(cè)區(qū)域土壤質(zhì)量演變提供翔實(shí)基礎(chǔ)底圖，也有助于明晰區(qū)域SOM變化影響因素，最終為區(qū)域土壤資源可持續(xù)利用決策提供支持。

1 材料與方法

1.1 研究區(qū)概況

杞麓湖盆地位于云南省通?？h（102°34′～102°52′E，24°6′～24°12′ N），屬亞熱帶半濕潤(rùn)高原季風(fēng)氣候區(qū)。區(qū)內(nèi)四季溫差小、早晚溫差大、干濕季分明、雨熱同季，多年平均降雨量900 mm，多年平均氣溫16℃，年日照時(shí)間2 274 h。地勢(shì)由西南向東北逐漸降低，坡度3°～10°，海拔1 796～1 820 m。盆地內(nèi)主要作物為蔬菜、糧食、油料、烤煙和花卉。全國(guó)第二次土壤普查劃分的土壤類(lèi)型為水稻土，水稻土亞類(lèi)為潴育型、潛育型、淹育型以及沼澤型[11]。

1.2 樣品采集與分析

為調(diào)查土壤屬性的空間分布，玉溪市土壤肥料工作站分別于2008年、2011年、2013年和2015年隨機(jī)采集了杞麓湖盆地耕作層（0～20 cm）679、804、484和362個(gè)土壤樣品（圖1）。并采用重鉻酸鉀（K2Cr2O7）氧化—滴定法測(cè)定了SOM含量[12-13]，提取了各樣點(diǎn)的土壤亞類(lèi)（潴育型、潛育型、淹育型、沼澤型）、土壤質(zhì)地（砂壤、壤土、黏壤、黏土）、農(nóng)業(yè)設(shè)施等級(jí)（不配套、基本配套、完全配套）以及灌溉程度（無(wú)灌、有灌溉潛力、能灌、保灌）信息。

1.3 離群值處理

離群值的存在會(huì)影響樣本的概率分布，從而影響根據(jù)樣本數(shù)據(jù)估計(jì)的半方差圖[14]。且回歸分析和普通克里格插值均在數(shù)據(jù)服從近似正態(tài)分布時(shí)能獲得更準(zhǔn)確的結(jié)果。因此采用Hample方法[15]根據(jù)樣本的中值（Median）和中值絕對(duì)偏差（Median Absolute Deviation，MAD）來(lái)識(shí)別并剔除四個(gè)時(shí)相SOM的離群值（outliers）?；驹瓌t是當(dāng)樣本中的數(shù)值大于（或小于）中值加上（或減去）3倍中值絕對(duì)偏差時(shí)，則視為離群值。剔除離群值后分別剩余667、788、469和359個(gè)土壤樣品。

圖1 杞麓湖盆地四個(gè)時(shí)相的土壤采樣點(diǎn)分布情況Fig. 1 Distribution of sampling sites in the Qiluhu basin relative to time phase

1.4 半方差模型擬合與克里格插值

為了將盡可能多的土壤樣點(diǎn)用于半方差模擬的擬合，從而獲得穩(wěn)定的半方差模型，本研究通過(guò)留一交叉驗(yàn)證（Leave One Out Cross Validation，LOOCV）比較土壤研究中常用的指數(shù)和橢球模型。選取的量化指標(biāo)包括誤差均值（Mean errors，ME）、誤差平方的均值（Mean squared errors，MSE）、克里格插值方差均值（Mean kriging variances，以及預(yù)測(cè)值與觀(guān)測(cè)值離差率平方的中值（Median squared deviation ratio，MSDR）。

1.5 統(tǒng)計(jì)年鑒數(shù)據(jù)收集與整理

從《通?？h統(tǒng)計(jì)年鑒》中收集并整理了杞麓湖盆地內(nèi)共47個(gè)自然村的2008年、2011年、2013年和2015年的蔬菜種植總面積、糧食種植總面積、氮肥施用量、磷肥施用量、鉀肥施用量、復(fù)合肥施用量數(shù)據(jù)，并將這些數(shù)據(jù)與村級(jí)行政區(qū)劃圖連接，形成空間屬性一體的矢量文件。統(tǒng)計(jì)每個(gè)自然村2008—2015年的蔬菜種植總面積平均變化量、糧食種植總面積平均變化量、化肥（氮、磷、鉀和復(fù)合肥）施用量平均變化量。

1.6 土壤有機(jī)質(zhì)影響因素分析

基于2008年、2011年、2013年和2015年四期土壤樣點(diǎn)數(shù)據(jù)，以SOM為因變量，將土壤亞類(lèi)、土壤質(zhì)地、農(nóng)業(yè)基礎(chǔ)設(shè)施和灌溉程度作為自變量并按各因子的等級(jí)賦予啞變量（dummy variables），進(jìn)行多元回歸分析。此外，基于四期SOM插值數(shù)據(jù)，以自然村為單元統(tǒng)計(jì)四個(gè)年份的SOM均值，并計(jì)算平均變化量。以四期SOM均值平均變化量結(jié)合各村的蔬菜總產(chǎn)量平均變化量、糧食總產(chǎn)量平均變化量和化肥（氮、磷、鉀和復(fù)合肥）施用平均變化量，進(jìn)行探索性回歸。該方法通過(guò)普通最小二乘法（Ordinary Least Squares，OLS）評(píng)估所有待選因子組合，并依據(jù)調(diào)整R2（AdjustR2）、系數(shù)顯著性P（CoefficientPvalue）、方差膨脹因子（Variance Inflation Factor，VIF）和雅克貝拉顯著性P（Jarque-BeraPvalue，JB）來(lái)選擇較為合適的回歸模型。

1.7 數(shù)據(jù)處理與分析

數(shù)據(jù)的探索性分析和多元回歸分析使用R軟件[16]的基礎(chǔ)包完成，土壤數(shù)據(jù)的半方差模型擬合及普通克里格插值使用R軟件的gstat軟件包[17]完成，探索性回歸和部分地圖的制作在ArcMap 10.5中完成。

2 結(jié) 果

2.1 SOM基本統(tǒng)計(jì)特征

2008年、2011年、2013年和2015年的原始樣本的SOM偏度和峰度均較大，概率分布為右向長(zhǎng)尾。雖然剔除離群值后剩余的土壤樣本（667、788、469和359）的概率分布并不完全服從正態(tài)分布，但四期樣本概率分布的偏度均小于0.5，而峰度均明顯降低，較處理前更接近于正態(tài)分布。用剔除離群值后的樣本的變異系數(shù)CV來(lái)表征數(shù)據(jù)的變異性，可見(jiàn)四期SOM的變異性均大于35%。根據(jù)Wilding[18]推薦的分類(lèi)指標(biāo)，當(dāng)CV＞35%時(shí)為強(qiáng)度變異性，CV＞15%為中等變異性，CV＜15%為輕微變異。

統(tǒng)計(jì)對(duì)比剔除離群值后四個(gè)時(shí)相的樣本，結(jié)果表明相對(duì)2008年，2011年的SOM均值顯著降低；而相對(duì)2011年，2013年和2015年SOM均值顯著增加。不同年份SOM均值差異的95%置信區(qū)間最大范圍為均值±3.13 g·kg-1。此值將作為下文中統(tǒng)計(jì)SOM變化面積的臨界參考值。

2.2 SOM的空間結(jié)構(gòu)特征

數(shù)量不小于100且在研究區(qū)內(nèi)分布均勻的采樣點(diǎn)可以較為穩(wěn)定且真實(shí)地估算區(qū)域變量的半方差函數(shù)從而反映其空間相關(guān)關(guān)系[19-20]。根據(jù)剔除離群值后樣本擬合四個(gè)時(shí)相SOM的半方差的指數(shù)型和橢球模型[21]的塊金值（Nugget）、偏基臺(tái)（Partial sill）、變程值（Range）以及塊金-基臺(tái)比（表1）。用ME、MSE、MSDR和R2比較指數(shù)和橢球模型的適用性。其中ME應(yīng)接近于0，而MSE則應(yīng)越小越好；衡量模型好壞的最主要指標(biāo)為MSDR和R2，其值越接近于1，對(duì)應(yīng)模型擬合的越好。結(jié)果表明指數(shù)模型較橢球模型能更好地?cái)M合2008年（MSDR和R2分別為0.988和0.986）和2011年（MSDR和R2分別為0.956和0.947）SOM的半方差函數(shù)，而橢球模型在擬合2013年（MSDR和R2分別為1.065和0.973）和2015年（MSDR和R2分別為1.032和0.988）SOM時(shí)略好于指數(shù)模型。

表1 根據(jù)杞麓湖盆地四個(gè)時(shí)相土壤有機(jī)質(zhì)調(diào)查值擬合的指數(shù)和橢球半方差模型的基本參數(shù)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時(shí)空分布特征

本研究采用普通克里格方法，基于指數(shù)模型預(yù)測(cè)2008年和2011年、橢球模型預(yù)測(cè)2013年和2015年SOM的空間分布，并輸出粒度為50 m×50 m的SOM分布圖（圖2）。可以看出，杞麓湖盆地SOM的空間分布特征在不同年份總體格局相似，但局部差異明顯?？傮w上四個(gè)時(shí)相的SOM從西南向東北呈現(xiàn)出低-高-低的分布特征。局部來(lái)看，盆地西北角的SOM從2008年至2011年有較小程度的降低，而在2013年和2015年均呈現(xiàn)出明顯的增加；中部高值核心區(qū)位置基本一致，但高值區(qū)的輻射范圍從2008年至2011年表現(xiàn)為收縮，而2011年至2013年和2015年則表現(xiàn)為逐漸擴(kuò)張。

為了明晰不同年份間SOM的空間變化情況，以2008年的SOM數(shù)據(jù)為基底，比較2011年、2013年和2015年的SOM空間變化（圖3）。相比2008年，SOM增加的區(qū)域面積有逐年擴(kuò)大的趨勢(shì)，增加的區(qū)域主要集中在盆地的西北角、西南角、南部以及杞麓湖的西北岸和南岸大部。據(jù)插值結(jié)果統(tǒng)計(jì)，杞麓湖盆地2011年、2013年和2015年相對(duì)于2008年的SOM的變化值的范圍為-21.87～17.65 g·kg-1。如依據(jù)上文中不同年份SOM均值差異的95%置信區(qū)間將0±3.13 g·kg-1確定為SOM減少或增加的臨界區(qū)，并將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五個(gè)等級(jí)，則各等級(jí)面積的變化情況明顯（圖4）。從2008年至2015年的三個(gè)時(shí)間段內(nèi)SOM連續(xù)增加的區(qū)域面積約占68.7%，為SOM連續(xù)減少的區(qū)域面積的兩倍多。

2.4 耕層土壤有機(jī)碳的影響因素

圖2 土壤有機(jī)質(zhì)的50 m×50 m預(yù)測(cè)值分布圖Fig. 2 Predicted SOM distribution map with 50 m×50 m resolution

圖3 四個(gè)時(shí)期土壤有機(jī)質(zhì)變化分布圖Fig. 3 Soil organic matter distribution map relative to time phase

圖4 杞麓湖盆地五個(gè)土壤有機(jī)碳變化等級(jí)的面積Fig. 4 Areas of the five levels of soil organic matter variation in the Qiluhu Basin

將土壤亞類(lèi)、土壤質(zhì)地、農(nóng)業(yè)基礎(chǔ)設(shè)施和灌溉程度作為自變量并賦予啞變量后進(jìn)行回歸分析。結(jié)果表明，除灌溉程度對(duì)土壤有機(jī)質(zhì)的貢獻(xiàn)不顯著外，土壤亞類(lèi)、土壤質(zhì)地和農(nóng)業(yè)基礎(chǔ)設(shè)施這三個(gè)因子均呈現(xiàn)出顯著的影響（P＜0.001）。其中，對(duì)杞麓湖盆地耕層SOM影響最大的是土壤亞類(lèi)的類(lèi)型，一共能解釋14.3%；而農(nóng)業(yè)設(shè)施和土壤質(zhì)地分別能解釋2.6%和1.3%。在水稻土亞類(lèi)中，潴育型、潛育型和沼澤型分別較淹育型高出10.94±0.74、21.12±1.59和25.36±1.37 g·kg-1。相對(duì)于沒(méi)有配套農(nóng)業(yè)設(shè)施的農(nóng)田，基本配套和完全配套的農(nóng)田耕層SOM均有所提高，但前者提高的程度未達(dá)到統(tǒng)計(jì)顯著性。相對(duì)于砂壤，SOM在壤土和黏壤中的含量顯著更高（P＜0.001），但黏土并未呈現(xiàn)出與砂壤顯著差異。

探索性回歸分析結(jié)果表明，在不同的候選因子組成的模型中，由糧食作物總產(chǎn)量平均變化量、蔬菜總產(chǎn)量平均變化量、氮肥施用量平均變化量和復(fù)合肥施用量平均變化量四個(gè)因子組成的模型的調(diào)整R2最大（0.23），其中前兩個(gè)因子達(dá)到顯著性（P＜0.001），而后兩個(gè)因子未達(dá)到顯著性。四個(gè)因子的VIF（2.95）并未超過(guò)閾值7.5，因子間的共線(xiàn)性不明顯。但JB值（P=0.000）表明該模型擬合后的殘差不符合正態(tài)性分布，說(shuō)明仍有除了以上四個(gè)因子之外的其他重要因子影響村級(jí)的SOM的分布。

3 討 論

若按全國(guó)第二次土壤普查時(shí)設(shè)定的標(biāo)準(zhǔn)，將耕層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標(biāo)準(zhǔn)劃分為六級(jí)[8]。那么，杞麓湖盆地耕層SOM均值除2011年處于第二級(jí)外，其余三個(gè)時(shí)相均為第一級(jí)，說(shuō)明整個(gè)盆地的耕層SOM均值處于較高的水平。與楊帆等[8]研究結(jié)果比較，杞麓湖盆地的SOM顯著高于全國(guó)耕層SOM平均水平（24.65 g·kg-1），接近我國(guó)耕層SOM最高的黑龍江省的平均水平（40.43 g·kg-1），顯著高于云南省耕層SOM平均水平（35.37 g·kg-1）。

從影響機(jī)制來(lái)看，糧食作物總產(chǎn)量平均變化量、蔬菜總產(chǎn)量平均變化量和氮肥施用量平均變化量對(duì)杞麓湖盆地耕層SOM的影響為負(fù)作用，而復(fù)合肥施用量平均變化量的影響為正作用。在杞麓湖盆地內(nèi)糧食和蔬菜產(chǎn)量持續(xù)增加的自然村，由于蔬菜的高復(fù)種指數(shù)（多數(shù)為四茬）、糧食作物的高秸稈移除率（除根系外基本全部移除）以及高度依賴(lài)化肥投入而輕視有機(jī)肥的施用導(dǎo)致了SOM的降低；而調(diào)整化肥施用中氮肥和復(fù)合肥的占比情況，適當(dāng)減少氮肥占比并增加復(fù)合肥的占比在一定程度上有利于維持甚至提升SOM含量。

根據(jù)已有研究，為了進(jìn)一步提高杞麓湖盆地的耕層SOM，可以嘗試有機(jī)肥與化肥的科學(xué)配比施用[22]；在地理環(huán)境類(lèi)似的洱海流域的研究結(jié)果表明，降低化肥施用量，同時(shí)增加有機(jī)肥的施用量能確實(shí)提高蔬菜地的耕層SOM[23]。此外，冬季休耕或少耕也被證明能有效地提高滇中地區(qū)耕層SOM總量[24]。當(dāng)然，科學(xué)地選擇耕作管理措施需要綜合考慮生態(tài)效益、農(nóng)民的經(jīng)濟(jì)效益及由此帶來(lái)的社會(huì)效益。因此，在未來(lái)的研究工作中應(yīng)加強(qiáng)長(zhǎng)期定位對(duì)比實(shí)驗(yàn)的研究，以此來(lái)共同促進(jìn)杞麓湖盆地耕層SOM的科學(xué)管理。

4 結(jié) 論

借助普通克里格插值方法能獲得杞麓湖盆地2008年、2011年、2013年和2015年四個(gè)時(shí)相SOM的空間分布圖。指數(shù)和橢球半方差模型在空間插值的精度評(píng)價(jià)方面差別較小，前者更適用于2008年和2011年的SOM數(shù)據(jù)，而后者更適用于2013年和2015年的數(shù)據(jù)。整個(gè)盆地的SOM分布從西南到中部再到東北部呈現(xiàn)出低-高-低的分布特征，并且在時(shí)間系列上SOM高值核心區(qū)的位置并未呈現(xiàn)顯著變化，但其輻射范圍有逐年擴(kuò)大的趨勢(shì)。從2008年到2015年的三個(gè)時(shí)間段內(nèi)，SOM呈現(xiàn)出連續(xù)增加的區(qū)域面積為SOM連續(xù)減小的區(qū)域面積的兩倍多。據(jù)土壤樣點(diǎn)調(diào)查數(shù)據(jù)分析，土壤亞類(lèi)類(lèi)型對(duì)整個(gè)區(qū)域內(nèi)SOM的變異起主要影響作用，可解釋變異的14.3%，而農(nóng)業(yè)設(shè)施和土壤質(zhì)地分別能解釋2.6%和1.3%。以自然村為基礎(chǔ)單元，SOM的降低主要由糧食和蔬菜的高產(chǎn)出所導(dǎo)致，而適當(dāng)?shù)臏p少氮肥施用量并增加復(fù)合肥施用量能在一定程度上維持甚至提升SOM含量。