杞麓湖盆地土壤有機質多時相空間分布與演變*

秦發侶，郭春平，陳 紅，柳德江?，田育天，謝新喬，李湘偉，王世航

（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含量。