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

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

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 樣品采集與分析

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

1.3 離群值處理

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

圖1 杞麓湖盆地四個時相的土壤采樣點分布情況Fig. 1 Distribution of sampling sites in the Qiluhu basin relative to time phase

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

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

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

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

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

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

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

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

2 結(jié) 果

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

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

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

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

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

表1 根據(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ù)模型預(yù)測2008年和2011年、橢球模型預(yù)測2013年和2015年SOM的空間分布，并輸出粒度為50 m×50 m的SOM分布圖（圖2）?？梢钥闯觯铰春璧豐OM的空間分布特征在不同年份總體格局相似，但局部差異明顯。總體上四個時相的SOM從西南向東北呈現(xiàn)出低-高-低的分布特征。局部來看，盆地西北角的SOM從2008年至2011年有較小程度的降低，而在2013年和2015年均呈現(xiàn)出明顯的增加；中部高值核心區(qū)位置基本一致，但高值區(qū)的輻射范圍從2008年至2011年表現(xiàn)為收縮，而2011年至2013年和2015年則表現(xiàn)為逐漸擴張。

為了明晰不同年份間SOM的空間變化情況，以2008年的SOM數(shù)據(jù)為基底，比較2011年、2013年和2015年的SOM空間變化（圖3）。相比2008年，SOM增加的區(qū)域面積有逐年擴大的趨勢，增加的區(qū)域主要集中在盆地的西北角、西南角、南部以及杞麓湖的西北岸和南岸大部。據(jù)插值結(jié)果統(tǒng)計，杞麓湖盆地2011年、2013年和2015年相對于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五個等級，則各等級面積的變化情況明顯（圖4）。從2008年至2015年的三個時間段內(nèi)SOM連續(xù)增加的區(qū)域面積約占68.7%，為SOM連續(xù)減少的區(qū)域面積的兩倍多。

2.4 耕層土壤有機碳的影響因素

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

圖3 四個時期土壤有機質(zhì)變化分布圖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

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

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

3 討 論

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

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

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

4 結(jié) 論

借助普通克里格插值方法能獲得杞麓湖盆地2008年、2011年、2013年和2015年四個時相SOM的空間分布圖。指數(shù)和橢球半方差模型在空間插值的精度評價方面差別較小，前者更適用于2008年和2011年的SOM數(shù)據(jù)，而后者更適用于2013年和2015年的數(shù)據(jù)。整個盆地的SOM分布從西南到中部再到東北部呈現(xiàn)出低-高-低的分布特征，并且在時間系列上SOM高值核心區(qū)的位置并未呈現(xiàn)顯著變化，但其輻射范圍有逐年擴大的趨勢。從2008年到2015年的三個時間段內(nèi)，SOM呈現(xiàn)出連續(xù)增加的區(qū)域面積為SOM連續(xù)減小的區(qū)域面積的兩倍多。據(jù)土壤樣點調(diào)查數(shù)據(jù)分析，土壤亞類類型對整個區(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含量。