大面積高寒山區(qū)土壤養(yǎng)分空間預(yù)測與管理分區(qū)①

杜龍全，劉 峰，史 舟，趙 霞，李德成，張甘霖,6，董晉鵬，陳東升

大面積高寒山區(qū)土壤養(yǎng)分空間預(yù)測與管理分區(qū)①

杜龍全1,2,3，劉 峰3,4，史 舟5，趙 霞1,2*，李德成3，張甘霖3,4,6，董晉鵬1,2，陳東升1,2

(1 青海師范大學(xué)，青海土壤數(shù)字服務(wù)中心，青海省自然地理與環(huán)境過程重點實驗室，西寧 810008；2 青海省人民政府–北京師范大學(xué)高原科學(xué)與可持續(xù)發(fā)展研究院–高原土壤信息科學(xué)研究團(tuán)隊，西寧 810008；3土壤與農(nóng)業(yè)可持續(xù)發(fā)展國家重點實驗室(中國科學(xué)院南京土壤研究所)，南京 210008；4 中國科學(xué)院大學(xué)，北京 100049；5 浙江大學(xué)環(huán)境與資源學(xué)院，杭州 310058；6 中國科學(xué)院南京地理與湖泊研究所流域地理學(xué)重點實驗室，南京 210008)

本文基于青海省土系調(diào)查的205個土壤剖面數(shù)據(jù)，利用隨機森林模型，分別建立了表層(0 ～ 20 cm)土壤全氮(TN)、全鉀(TK)、全磷(TP)與環(huán)境因素變量(地形、氣候、植被、遙感數(shù)據(jù))之間的定量關(guān)系，對青海省土壤養(yǎng)分含量的空間分布進(jìn)行了預(yù)測；在此基礎(chǔ)上，結(jié)合全國土壤養(yǎng)分的分級標(biāo)準(zhǔn)，利用投影尋蹤法，生成了土壤養(yǎng)分的管理分區(qū)。留一交叉驗證結(jié)果顯示，TN、TK、TP空間預(yù)測的2分別是0.89、0.85、0.82，可解釋土壤養(yǎng)分空間變異的80% 以上，表明隨機森林模型和環(huán)境因素變量結(jié)合可以在稀疏樣本條件下有效地預(yù)測大面積高寒山區(qū)土壤養(yǎng)分空間變異。青海省土壤養(yǎng)分呈現(xiàn)東高西低的分布模式，土壤綜合養(yǎng)分高等級出現(xiàn)在南部的玉樹、果洛、黃南和東部的湟水谷地地區(qū)；低等級主要分布在柴達(dá)木盆地、可可西里、海南州中南部；全省土壤綜合養(yǎng)分分級均在中上等級以上，占全省面積的81%。分析發(fā)現(xiàn)，植被是影響青海省表層土壤養(yǎng)分TN、TP、TK空間分布的主要環(huán)境因素，其中年均降水量、地表溫度是影響青海省表層土壤TN空間模式的重要因素；地表覆被、海拔和地表溫度等環(huán)境因子對表層土壤TP的空間變異占主導(dǎo)作用；年均降水量、晝夜溫差是影響表層土壤TK的空間模式的重要因素。

土壤養(yǎng)分；隨機森林；投影尋蹤法；管理分區(qū)；數(shù)字土壤制圖

土壤養(yǎng)分是土壤提供的植物生長所必需的營養(yǎng)元素[1]，是土壤生產(chǎn)力的基礎(chǔ)，了解土壤養(yǎng)分的空間變異，對于研究土壤退化、生物多樣性缺失、氣候變化和制定土地管理計劃具有重要意義[2-3]。我國西部省區(qū)，尤其是高寒山區(qū)，土壤景觀關(guān)系復(fù)雜，土壤空間變異程度高[4]，但缺乏準(zhǔn)確的土壤空間信息，制約著農(nóng)牧業(yè)的科學(xué)生產(chǎn)與管理。

數(shù)字土壤制圖是將土壤–景觀模型通過計算機輔助生成的以柵格數(shù)據(jù)為表達(dá)形式的制圖方法，可反映較為準(zhǔn)確的土壤空間信息。制圖模型主要包括地統(tǒng)計、專家知識和機器學(xué)習(xí)方法，研究內(nèi)容集中于土壤的多種屬性。在區(qū)域從小到大、景觀從簡單到復(fù)雜、制圖結(jié)果從二維到三維(多維)的制圖研究中，環(huán)境變量選取的是相對較少且易獲取的協(xié)變量集，主要集中于地形因子與氣候因子以及用遙感影像來表征的生物因子。涉及大面積區(qū)域時，成土母質(zhì)這一因子的獲取較難，往往被忽略；同時，表征土壤和環(huán)境之間關(guān)系的一些成土因子定量化很難表示(如時間因子)[5-12]。因此，要實現(xiàn)大尺度、高精度的制圖，開發(fā)挖掘豐富的環(huán)境協(xié)變量數(shù)據(jù)庫是一個重要的突破口。在現(xiàn)有的機器學(xué)習(xí)(集成)模型甚至是人工智能模型中，只能反映小區(qū)域、局部景觀環(huán)境與土壤的關(guān)系，模型對于與研究區(qū)環(huán)境條件相似的地區(qū)可能適用性較好，對于復(fù)雜的大面積區(qū)域適用性較低[12-14]。要精確地模擬土壤性質(zhì)和過程，有必要開發(fā)具有處理不同區(qū)域土壤與環(huán)境之間非線性關(guān)系的針對性的模型。目前國內(nèi)數(shù)字土壤制圖工作主要集中在我國中部和東部地區(qū)，地形梯度相對較小，可達(dá)性較好，土壤調(diào)查樣本通常數(shù)量較多且分布相對均衡。例如，楊煜岑等[5]在陜西省藍(lán)田縣用多元線性回歸模型預(yù)測了陜西省藍(lán)田縣土壤養(yǎng)分的空間分布；Dong等[6]將機器學(xué)習(xí)算法與基于地塊的DSM框架相結(jié)合，對寧夏沖積平原農(nóng)業(yè)區(qū)表層土壤養(yǎng)分(土壤有機質(zhì)、全氮等)特征進(jìn)行了預(yù)測；Liu等[7]用地表動態(tài)反饋方法預(yù)測了江蘇省中部平原地區(qū)土壤有機質(zhì)含量等屬性的空間分布。然而，在我國西部大面積高寒山區(qū)，路網(wǎng)不發(fā)達(dá)，可達(dá)性較差，土壤調(diào)查樣本較為稀疏，已有工作主要集中在景觀或流域尺度以及單一土壤養(yǎng)分屬性(如有機碳或全氮)，如代子俊等[8]針對湟水流域分析了近30 a的土壤全氮的時空變異影響因素；張歡和高小紅[9]基于湟水流域土壤有機質(zhì)屬性對比分析了地統(tǒng)計學(xué)進(jìn)行預(yù)測制圖的精度；周偉等[10]在三江源地區(qū)將原始光譜反射率數(shù)據(jù)及其不同數(shù)據(jù)變換形式下的光譜分別與土壤有機質(zhì)含量進(jìn)行了相關(guān)分析；龐龍輝等[11]對青海省土壤全氮、有機碳、粉粒和pH進(jìn)行了預(yù)測制圖。如何對大面積高寒山區(qū)土壤養(yǎng)分空間分布準(zhǔn)確地預(yù)測，并在此基礎(chǔ)上面向農(nóng)牧業(yè)進(jìn)行合理的生產(chǎn)管理分區(qū)，是一個值得研究的問題。

本文以地處青藏高原東北部的青海全省為研究區(qū)，將機器學(xué)習(xí)方法與成土環(huán)境協(xié)同變量相結(jié)合，探討土壤養(yǎng)分(全氮(TN)、全鉀(TK)、全磷(TP))含量的空間預(yù)測制圖及生產(chǎn)管理分區(qū)，以期為青海省土壤養(yǎng)分空間管理及農(nóng)牧業(yè)可持續(xù)發(fā)展提供必要參考。

1 材料與方法

1.1 研究區(qū)概況

青海省地處我國西北部內(nèi)陸腹地，地理位置為 31°39′ ～ 39°19′N、89°35′ ～ 103°04′E，總面積72.23萬km2，大部分區(qū)域為高寒山區(qū)，地勢總趨勢為南高北低[15]，海拔介于1 686 ～ 6 477 m，大部分地區(qū)海拔在3 000 m以上。地貌上復(fù)雜多樣兼有高原、山地、草原、湖泊；氣候上，雖然總體上屬于高原大陸性氣候，但屬暖溫帶向亞熱帶過渡型，不同海拔高度和不同地區(qū)氣候有明顯差異，年均溫度 –5.1 ～ 9.0 ℃，年均降水量15 ～ 750 mm，降水量由東向西逐漸遞減。按大地構(gòu)造全省大體可分為三類不同地形區(qū)：南部青南高原、西部柴達(dá)木盆地和北部與東部平行嶺谷區(qū)[4]。區(qū)域內(nèi)成土母質(zhì)類型較多，成因各異。受地形、氣候、土壤環(huán)境的影響，草地面積占全省的3/5，畜牧業(yè)得以發(fā)展，全省的森林、草原、濕地等在應(yīng)對氣候變化上有顯著的作用，是我國重要的生態(tài)屏障。研究區(qū)采樣點空間分布見圖1。

圖1 研究區(qū)采樣點空間分布

1.2 數(shù)據(jù)來源與處理

土壤剖面數(shù)據(jù)來源于科技部基礎(chǔ)專項“我國土系調(diào)查與《中國土系志·青海卷》編制”項目中的205個樣點，主要包括采樣點位置和景觀信息以及土壤理化性質(zhì)(TN，TP，TK，SOC，pH，黏粒、粉粒、砂粒含量等)。土壤TN采用重鉻酸鉀、硫酸消化–蒸餾法(凱氏蒸餾法)測定；TK使用堿熔–火焰光度法測定；TP 使用堿熔–鉬銻抗比色法測定[16]。對剖面各發(fā)生層測定的數(shù)據(jù)進(jìn)行深度加權(quán)，生成表層(0 ～ 20 cm)的TN、TK、TP數(shù)據(jù)。

環(huán)境變量的選取是以土壤發(fā)生學(xué)理論為依據(jù)，主要依照土壤形成的五大要素(母質(zhì)、地形、生物、氣候和時間)進(jìn)行選取代表性環(huán)境變量。由于目前尚無統(tǒng)一的定量指標(biāo)來描述時間信息，因此本研究暫不考慮時間類因子。用土壤類型(soil types)近似表征母質(zhì)的差異；地形因子在土壤發(fā)生與演變過程中影響著地表物質(zhì)能量的循環(huán)轉(zhuǎn)換；遙感數(shù)據(jù)具有實時、易獲取的特點，其特征波段及不同的組合可以用來表征生物景觀信息的綜合環(huán)境(地表景觀信息)。氣候因子(降水量、氣溫)是直接影響成土過程中的水熱條件。本研究通過表層土壤屬性及其與環(huán)境變量之間的相關(guān)性確定了比較關(guān)鍵的21個環(huán)境變量。其中，地形因子：DEM數(shù)據(jù)來自于資源環(huán)境科學(xué)與數(shù)據(jù)中心(http://www.resdc.cn/)，基于DEM數(shù)據(jù)利用SAGA軟件(https://sourceforge.net/)計算提取海拔(Elev)、坡度(Slp)、坡向(Asp)、剖面曲率(Curprf)、地形濕度指數(shù)(Twi)、地形開放度(Openness)、多尺度山谷指數(shù)(Mrvbf)；氣候數(shù)據(jù)：年均溫(MAT)、年均降水量(MAP)、氣溫季節(jié)性(Ts)、晝夜溫差(Mdr)、潛在太陽輻射(Radition)數(shù)據(jù)來自于世界氣象數(shù)據(jù)庫網(wǎng)[17](http://worldclim.org/version2)；土壤類型(Soiltype)和土地覆蓋類型(Lucc)均來源于資源環(huán)境科學(xué)與數(shù)據(jù)中心(http://www.resdc.cn/)；植被指數(shù)均值(Evm)，植被指數(shù)標(biāo)準(zhǔn)差(Evs)，白天地表溫度(LSTday)，夜間地表溫度(LSTnight)，7、8月地表溫度(LST78)和地表反射率(band4/7)的處理是在Google Earth Engine (GEE)云平臺通過JavaScript編程完成，選取時間為2015年7—8月的MODIS影像數(shù)據(jù)集，處理包括影像去云、掩膜裁剪、投影變換。所有柵格環(huán)境變量均重采樣(最近鄰法)為1 km。

1.3 土壤養(yǎng)分空間預(yù)測方法

隨機森林[18]是通過建立多個決策樹并將其融合成一個更加準(zhǔn)確和穩(wěn)定的模型，能夠根據(jù)樣本之間的相似性來衡量特征交互，廣泛應(yīng)用于分類和回歸。隨機森林建模在R.4.1中Random Forest包進(jìn)行[19]，該模型主要有3個參數(shù)需要定義，即樹的數(shù)量(ntree)、葉節(jié)點包含最少樣本的個數(shù) (nodesize)、隨機選取變量的個數(shù) (mtry)。使用caret包對nodesize和mtry兩個參數(shù)進(jìn)行優(yōu)化，優(yōu)化值分別是5和3，ntree采用默認(rèn)值500，足以產(chǎn)生穩(wěn)定的預(yù)測。%IncMSE作為判斷預(yù)測變量重要性的指標(biāo)，該值越大表示該變量的重要性越大，即對土壤養(yǎng)分含量的影響度就越大。

本文中制圖預(yù)測精度采用留一交叉驗證法LOOCV(leave one out cross validation)。其原理是：在個數(shù)為數(shù)據(jù)集中隨機抽取1個樣本去驗證-1個數(shù)據(jù)所訓(xùn)練的模型，重復(fù)次取平均值[20]。評價指標(biāo)包括：決定系數(shù)2、一致性指數(shù)(CC)、均方根誤差(RMSE)、預(yù)測偏差值(ME)。其中，2表示模型對自變量解釋的比例，取值范圍0～1；CC表示預(yù)測誤差的方差值占實測值的比，取值范圍0～1，1表示預(yù)測結(jié)果最佳，0表示結(jié)果不太符合；RMSE是預(yù)測值與實測值之間的均方根誤差；ME是預(yù)測值與實測值之間誤差值的平均，反映預(yù)測結(jié)果情況是否存在高估或低估。公式如下：

1.4 土壤養(yǎng)分分區(qū)評價方法

投影尋蹤是將高維數(shù)據(jù)向低維空間投影，通過分析低維空間的投影特性來研究高維數(shù)據(jù)的特征，找到反映數(shù)據(jù)結(jié)構(gòu)特征的最優(yōu)投影，是處理多指標(biāo)問題的統(tǒng)計方法[21]。其具體思路是將影響問題的多因素指標(biāo)通過投影尋蹤分析得到反映其綜合指標(biāo)特性的最優(yōu)投影特征值，然后建立投影特征值與因變量的一一對應(yīng)關(guān)系，通過分析一維的投影值來對樣本做出更加合理的分級和評價[21]。

然后，使用數(shù)量為的多維單位向量(1,2, …, a)表示某一投影方向，用公式(5)可計算樣本在該方向上的投影值，用公式(6)計算投影值的標(biāo)準(zhǔn)差()，用公式(7)計算投影值的局部密度()，用投影值的標(biāo)準(zhǔn)差() 和局部密度() 的乘積表示為投影指標(biāo)函數(shù)()，即公式(8)，目的是當(dāng)把維數(shù)據(jù)綜合成某一方向上的一維投影值時，投影值散布特征為局部投影點盡可能密集，而整體上各個點團(tuán)之間又盡可能分散[21]。

()()() (8)

考慮到()值與投影方向有關(guān)，當(dāng)()取最大值時的方向是最能反映數(shù)據(jù)結(jié)構(gòu)特征的方向，因而利用遺傳算法，通過求解投影指標(biāo)函數(shù)優(yōu)化問題識別最佳投影方向，公式如下

最后，基于最佳投影方向采用公式(5)計算全國土壤養(yǎng)分含量分級標(biāo)準(zhǔn)值的投影值，建立基于投影值的等級評價和分區(qū)方法。

依據(jù)全國第二次土壤普查的分級標(biāo)準(zhǔn)(表1)[22]，將土壤TN、TP、TK的預(yù)測結(jié)果進(jìn)行疊加分析、歸一化處理，基于DPS9.5軟件(http://www.dpsw.cn/)計算最優(yōu)投影方向以及最優(yōu)投影方向的投影值，較為客觀地對土壤養(yǎng)分分區(qū)進(jìn)行劃分。

表1 全國土壤養(yǎng)分含量分級

2 結(jié)果與分析

2.1 表層土壤養(yǎng)分的統(tǒng)計特征

表2為表層樣本的TN、TP、TK的描述性統(tǒng)計情況。Kolmogorov-Smirnov 檢驗結(jié)果顯示，TN、TP經(jīng)對數(shù)轉(zhuǎn)換后均符合正態(tài)分布，TK含量符合正態(tài)分布。樣點表層土壤TN含量范圍0.04 ～ 13.50 g/kg，平均值2.56 g/kg；TP含量范圍0.49 ～ 3.39 g/kg，平均值1.52 g/kg；TK范圍4.67 ～ 26.95 g/kg，平均值17.23 g/kg。從變異系數(shù)來看，表層土壤養(yǎng)分的變異系數(shù)均處于0.1 ～ 1，表明土壤養(yǎng)分均呈中等程度變異，其排序為TN>TP>TK，其中TN的變異程度最大，達(dá)到0.94。

2.2 土壤養(yǎng)分預(yù)測制圖與結(jié)果

利用隨機森林模型預(yù)測青海省表層土壤養(yǎng)分空間分布(圖 2) 表明，土壤TN預(yù)測值的范圍為0.20 ～ 8.67 g/kg，均值2.48 g/kg，在空間上整體呈現(xiàn)東南向西北遞減的趨勢，其中青海省東南、中部和東北部土壤TN含量較高，西北部及湟水谷地的含量較低。土壤TP的預(yù)測結(jié)果范圍0.75 ～ 2.41 g/kg，均值1.33 g/kg，空間分布上的趨勢是東高西低，土壤TP含量高的區(qū)域是在玉樹地區(qū)的東南部、東北部的祁連山地區(qū)和全省東部，而西南部可可西里地區(qū)、西北部的茫崖地區(qū)和中部昆侖山區(qū)含量較低。土壤TK預(yù)測的結(jié)果范圍12.70 ～ 22.25 g/kg，均值16.94 g/kg，整體呈現(xiàn)分布趨勢為東南高西北低，空間分布上較為分散，土壤TK含量高值主要分布在青南高原和祁連山山區(qū)，在唐古拉山中段地區(qū)也呈現(xiàn)高值，低值主要分布柴達(dá)木地區(qū)、昆侖山脈南部、黃南州的西北部和青海省西南部的格爾木市地區(qū)。

表2 表層土壤養(yǎng)分含量描述性統(tǒng)計

圖 2 土壤養(yǎng)分空間預(yù)測圖(A. TN；B. TP；C. TK)

交叉驗證結(jié)果(表3)表明，土壤TN、TP、TK模型擬合度2比較高，分別為0.89、0.85、0.82，即可以解釋土壤養(yǎng)分含量的80% 以上；從一致性系數(shù)(CC)來看，三者均在0.87以上；預(yù)測偏差(ME)三者接近0值，說明預(yù)測的值有被低估。由此，機器學(xué)習(xí)方法與成土環(huán)境協(xié)同變量相結(jié)合對青海省表層土壤養(yǎng)分的預(yù)測，可以很好地解釋土壤養(yǎng)分的空間變異。

表3 土壤養(yǎng)分預(yù)測的交叉驗證結(jié)果

2.3 土壤養(yǎng)分空間變異的影響因素

圖3顯示了土壤養(yǎng)分空間預(yù)測中環(huán)境協(xié)同變量的相對重要性排序。對于青海省表層土壤TN的空間分布格局，植被指數(shù)均值(Evm)為主導(dǎo)因子，年均降水量(MAP)次之。通常在空間范圍較大研究區(qū)域內(nèi)，氣候因素的地帶性可以直觀反映作物的生長狀況，間接影響著土壤養(yǎng)分元素的遷移轉(zhuǎn)化過程。在青海省東北部的祁連山地、湟水谷地及果洛玉樹高原區(qū)域的降水量較高，土壤水分條件好，植被覆蓋度高，且隨海拔增加，覆蓋率增加，地表的枯枝落葉增加，腐殖質(zhì)層較厚，土壤微生物活性較強，土壤養(yǎng)分來源更加豐富，積累的TN含量普遍較高。在青海省西部地區(qū)，氣候干旱炎熱，降水量稀少且地表溫度高，水熱條件的不均衡導(dǎo)致植被覆蓋度小，影響著土壤中微生物的活動，TN含量低。其次，西北區(qū)域分布著廣泛大面積的鹽漬土、風(fēng)沙土、砂質(zhì)土等，且其土壤保水保肥能力較差，土壤養(yǎng)分分解速度快，植物受高溫干旱和鹽分制約生存，氮素的來源和累積弱，土壤TN較低。

從TP來看，植被指數(shù)均值(Evm)、地表覆蓋(Lucc)和海拔(Elev)對其含量重要性排名前三。氣候、地形地貌等環(huán)境因子對土壤中磷元素的循環(huán)產(chǎn)生重要影響，玉樹地區(qū)的東南部、東北部的祁連山地區(qū)為青海省的植被覆蓋度高的區(qū)域，土壤中有機質(zhì)高，可減少磷的吸附位點，增加磷的移動性，其次土壤磷形態(tài)和含量受土壤微生物活性的直接影響，土壤中的微生物，能夠通過其生理活動促進(jìn)土壤中含磷化合物中的磷釋放，從而改變土壤中磷元素形態(tài)和含量。中北部柴達(dá)木盆地，存在綠洲所以出現(xiàn)高值；西南部可可西里地區(qū)、西北部的茫崖地區(qū)和中部昆侖山區(qū)，氣候干旱，同時受海拔、風(fēng)力的作用增強，土壤中砂粒含量高，溫度流失較快，土壤磷素的轉(zhuǎn)換率低。青南高原河谷地區(qū)降水量相對充足，此外海拔(Elev)對降水(MAP)分布模式產(chǎn)生影響，進(jìn)而土壤 pH 空間分布也受之影響，同時氣溫受海拔的影響，在此區(qū)域的蒸發(fā)量遠(yuǎn)低于降水量，使得土壤淋溶過程加速，促進(jìn)土壤酸性升高，在中、酸性的土壤中，磷素形態(tài)轉(zhuǎn)換速率快，土壤TP含量較高；西部柴達(dá)木盆地氣候干燥，蒸發(fā)量大，容易形成鹽沼，導(dǎo)致碳酸鹽集聚，土壤TP含量偏低；在東北部的河湟谷地區(qū)域土壤TP含量高，局部變異也大，該區(qū)域是青海省主要農(nóng)業(yè)區(qū)，長時間的人類農(nóng)耕活動作用，導(dǎo)致土壤堿性增強，土壤TP含量高。

(Soil type：土壤類型；Lucc：地表覆蓋類型；Evm：植被指數(shù)均值；Evs：植被指數(shù)標(biāo)準(zhǔn)差；Elev：海拔；Slp：坡度；Asp：坡向；Curprf：剖面曲率；Radiation：潛在太陽輻射；Twi：地形濕度指數(shù)；Openness：地形開放度；Mrvbf：多尺度山谷指數(shù)；MAP：年均降水量；MAT：年均溫；Mdr：晝夜溫差；Ts：氣溫季節(jié)性；LSTday：白天地表溫度；LSTnight：夜間地表溫度；LST78：7、8月地表溫度； band4、band7：地表反射第4和第7波段)

Fig. 3 Relative importance of environmental covariance in spatial prediction of soil nutrients

對于TK而言，植被指數(shù)均值(Evm)、年均降水量(MAP)、晝夜溫差(Mdr)的影響作用很大。青海省南部的玉樹、東南部的果洛地區(qū)和東北部的祁連山山區(qū)植被覆蓋度高，植被覆蓋較好，土壤受淋溶作用相對較弱，加之林地土壤中的充足草木枯萎被微生物分解，有利于土壤速效鉀的涵養(yǎng)與保持。中部的昆侖山區(qū)和唐古拉山中段地區(qū)的晝夜溫差(Mdr)較大，風(fēng)雪氣候常見，存在嚴(yán)重的土壤侵蝕，由于土壤全鉀主要是由礦物鉀組成，土壤鉀素會隨土壤侵蝕流失，土壤中TK含量較低；在西南部的格爾木地區(qū)，海拔較高，受風(fēng)力與水流的風(fēng)化作用較強，風(fēng)化物中又富含較多的云母、長石等含鉀結(jié)晶體，富鉀風(fēng)化殘積物混入土壤后使得含量變高。潛在太陽輻射(Radiation)、年均降水量(MAP)等協(xié)變量因子對土壤的濕潤和淹水時間持續(xù)的長短起到關(guān)鍵的作用，進(jìn)而影響鉀素固定和釋放過程，土壤失水變干，會產(chǎn)生釋鉀現(xiàn)象，青海省地形起伏度大，同時受氣候和人類活動影響，由此全省的TK含量局部上有所變異。

2.4 土壤養(yǎng)分分區(qū)結(jié)果

從表4可知，土壤養(yǎng)分分級標(biāo)準(zhǔn)值對應(yīng)的投影值分別為0.141、0.234、0.411、0.682和0.814，據(jù)此對土壤養(yǎng)分投影值劃分等級，其結(jié)果為：1級(<0.141)、2級(0.141 ～ 0.234)、3級(0.234 ～ 0.411)、4級(0.411 ～ 0.682)、5級(0.682 ～ 0.814)、6級(>0.814)。

表4 土壤養(yǎng)分的投影特征值及土壤養(yǎng)分分級

圖4為青海省土壤養(yǎng)分分區(qū)圖，表5是每個等級和分區(qū)所占面積的比重。土壤養(yǎng)分等級劃分為極低、低、中、中上、高、極高6個等級，管理分區(qū)分為6個分區(qū)(Q1 ～ Q6)。全省土壤養(yǎng)分等級劃分的大致空間格局為：整體上呈現(xiàn)東南-西北均勻遞減趨勢，極高的主要分布在玉樹的南部地區(qū)、果洛地區(qū)和東北部的祁連山，其分布區(qū)域植被覆蓋度高；昆侖山段也存在極高的分布，面積比例占全省的4.18%，這表明該地區(qū)的土壤人為活動較少，主要是自然發(fā)育的土壤，主要受氣候因子和地形因子影響作用較大；養(yǎng)分等級高的區(qū)域面積比例占全省的面積的13%，其空間分布不均，主要分布在海北州、格木市飛地、玉樹州的西南部和黃南州的西北部等地區(qū)，這些地區(qū)有大面積的林區(qū)、草場和耕地，同時受人類活動作用影響比較大；養(yǎng)分等級中至中上的區(qū)域主要分布西北柴達(dá)木盆地周邊、東部湟水谷地、青海湖周圍以及玉樹州的西南區(qū)域，占面積比例65%，這與該區(qū)域內(nèi)的地形坡度變化及人工施肥量有關(guān)；土壤養(yǎng)分低等級的區(qū)域主要分布在柴達(dá)木盆地、西南的可可西里區(qū)域、扎陵湖和鄂陵湖周邊和海南州的中北部，占面積比16%。土壤養(yǎng)分等級從高到低的面積比例為：4︰13︰40︰24︰16︰3，養(yǎng)分等級低和極低的總面積占全省的19%，表明青海省土壤養(yǎng)分屬于中等級及偏上，主要分布在全省的中南部和東北部。因此，需要因地制宜，科學(xué)地劃分土壤利用區(qū)域，調(diào)整放牧和草場施肥結(jié)構(gòu)，加強對土壤養(yǎng)分含量高的區(qū)域的管理與調(diào)整，恢復(fù)土壤養(yǎng)分含量低的區(qū)域，增加生態(tài)保護(hù)措施。

圖4 土壤養(yǎng)分分區(qū)圖

表5 土壤養(yǎng)分分區(qū)及面積比例

3 討論

通常在研究區(qū)面積相對較小、景觀相對單一、樣點較多且分布相對均衡的條件下，只利用樣點而不用環(huán)境變量進(jìn)行簡單插值預(yù)測，預(yù)測效果比較好。如張鐵嬋等[23]和石小華等[24]采用空間插值法分別對榆林市榆陽區(qū)和周至縣北部獼猴桃適生區(qū)的土壤養(yǎng)分元素進(jìn)行空間預(yù)測，在樣點布設(shè)大的地區(qū)，插值精度很好，反之插值精度不是很理想。插值結(jié)果受樣點布設(shè)密度的影響較大，而預(yù)測精度受不同的插值方法的影響較小。對于大面積高寒地區(qū)，樣點稀少且分布不均，顯然單純基于樣點的插值不適合，需要發(fā)揮環(huán)境因素輔助變量的作用。本文研究發(fā)現(xiàn)，通過把豐富的成土環(huán)境因素變量與復(fù)雜的機器學(xué)習(xí)方法結(jié)合，能夠在大面積高寒山區(qū)實現(xiàn)較為準(zhǔn)確的土壤養(yǎng)分預(yù)測制圖，捕捉土壤空間變異。植被能有效地指示土壤養(yǎng)分在大面積區(qū)域的空間變異。植被指數(shù)均值(Evm)對青海省表層土壤養(yǎng)分的空間變異具有主導(dǎo)作用。Pouladi等[25]在丹麥日德蘭半島，利用NDVI、SAVI、DVI、RVI和地形相關(guān)變量預(yù)測表層SOC含量，發(fā)現(xiàn)這些遙感變量解釋了89% 的SOC空間變化。Wang等[26]在遼寧省表層土壤含量制圖中發(fā)現(xiàn)結(jié)合僅利用遙感數(shù)據(jù)的隨機森林模型預(yù)測效果明顯好于僅使用地形和氣候變量的隨機森林模型。在后續(xù)工作中，在景觀復(fù)雜且草地大面積發(fā)育的青海省，遙感數(shù)據(jù)作為主要環(huán)境變量的土壤制圖值得探究。地表溫度是地表對于太陽輻射的響應(yīng)或反饋[27]，可以作為環(huán)境變量有效地對土壤環(huán)境的空間差異進(jìn)行揭示，對于本研究的高寒山區(qū)，地形地貌復(fù)雜多變，降水少，植被以牧草為主，這種反饋將主要取決于太陽輻射、地形、降水、植被和土壤的綜合作用。這也說明該區(qū)域地表溫度(LST)和年均降水量(MAP)是影響表層土壤養(yǎng)分空間變異的主要環(huán)境因子。地形因子的貢獻(xiàn)率在青海省表層土壤養(yǎng)分的空間變異中相對較低，主要的原因可能是在景觀和小流域尺度上，母質(zhì)、氣候和植被相對一致，地形通過徑流、土壤侵蝕間接影響土壤屬性的空間分布[28]，地形是影響小面積區(qū)域土壤屬性的主要因素，如孫孝林等[29]和楊琳等[30]在小區(qū)域土壤養(yǎng)分預(yù)測中只考慮地形屬性變量就可以進(jìn)行有效的土壤制圖。在大面積高寒山區(qū)，土壤屬性受到多種環(huán)境因素的協(xié)同影響，如母質(zhì)、氣候、地形和植被等。各成土因素變量對土壤養(yǎng)分都有貢獻(xiàn)，有較強的成土因素綜合作用。土壤屬性與環(huán)境變量之間存在非常復(fù)雜的關(guān)系，這與小面積區(qū)域土壤制圖不同。在后續(xù)的工作中，需要發(fā)掘更為有效的成土環(huán)境因素協(xié)同變量，來提高高寒山區(qū)數(shù)字土壤制圖的準(zhǔn)確性。

另一方面，預(yù)測模型選取應(yīng)考慮具體的土壤樣點數(shù)據(jù)和土壤景觀情況。大面積高寒山區(qū)地貌景觀復(fù)雜多樣，土壤空間變異程度高，樣點比較稀少且空間不均衡，導(dǎo)致這類地區(qū)數(shù)字土壤制圖結(jié)果難度往往比較大，簡單景觀是指某個成土因素主導(dǎo)土壤空間分布，土壤–環(huán)境關(guān)系近于線性。復(fù)雜景觀是指多個成土環(huán)境因素，土壤–環(huán)境關(guān)系是非線性的、非平穩(wěn)的[31]。對于前者，線性回歸可以較好地解釋協(xié)變量(如地形屬性)與土壤屬性之間的關(guān)系。對于后者，可以使用機器學(xué)習(xí)算法來模擬土壤變化。青海省作為大面積高寒地區(qū)，土壤-環(huán)境關(guān)系復(fù)雜，因此預(yù)測模型的選取至關(guān)重要，張歡和高小紅[9]比較了地統(tǒng)計中協(xié)同克里格、地理加權(quán)回歸克里格以及結(jié)合土地利用類型的均值、中值修正協(xié)同克里格等方法對湟水流域土壤有機質(zhì)屬性的預(yù)測精度，忽略了復(fù)雜的土壤環(huán)境關(guān)系，在空間分布上能體現(xiàn)其變化規(guī)律，但預(yù)測值結(jié)果存在預(yù)高估和低估現(xiàn)象甚至負(fù)值。地統(tǒng)計學(xué)模型在土壤樣本均勻分布且滿足空間自相關(guān)的條件下，模型的預(yù)測精度較好，這在地形復(fù)雜、土壤調(diào)查樣本較為稀疏的青海省并不適用。

土壤的利用和管理需要對土壤養(yǎng)分進(jìn)行生產(chǎn)管理分區(qū)。通常土壤屬性指標(biāo)和分區(qū)方法的選取不同，會忽略土壤綜合養(yǎng)分的空間差異、土壤質(zhì)量因子的變異性以及土壤環(huán)境的特殊性，進(jìn)而管理分區(qū)結(jié)果也因此產(chǎn)生不同，大部分學(xué)者采用了主成分分析、地統(tǒng)計、神經(jīng)網(wǎng)絡(luò)和模糊聚類模型等方法結(jié)合土壤屬性的單因子來進(jìn)行分區(qū)管理，黃婷等[32]通過主成分分析來計算土壤綜合質(zhì)量指數(shù)進(jìn)行了線性回歸分析，忽略了土壤-景觀關(guān)系的復(fù)雜多變性；朱昌達(dá)等[33]運用地統(tǒng)計方法分析土壤屬性空間變異的同時采用模糊聚類來進(jìn)行分區(qū)，未能考慮土壤屬性的綜合作用，計算分區(qū)指標(biāo)權(quán)重和綜合值。由于土壤綜合養(yǎng)分具有多維性和綜合性，而綜合定量化研究在較大程度上可以避免研究者主觀因素的影響，已然成為土壤綜合養(yǎng)分分區(qū)與管理的一個趨勢[34-36]。對于大面積高寒山區(qū)，土壤養(yǎng)分與環(huán)境變量之間的關(guān)系較為復(fù)雜，需要建立土壤養(yǎng)分與評價值之間較為準(zhǔn)確的關(guān)系，重視土壤養(yǎng)分的綜合指標(biāo)的選取信息，來提高管理分區(qū)結(jié)果的準(zhǔn)確性，本研究基于數(shù)字土壤制圖得到的土壤養(yǎng)分分布圖，采用投影尋蹤法對土壤養(yǎng)分進(jìn)行了管理分區(qū)。盡管投影尋蹤法減少了人為因素干擾，提高了養(yǎng)分管理分區(qū)的客觀和準(zhǔn)確性，但是不同投影方向優(yōu)選算法可能會影響評價結(jié)果，有待進(jìn)一步探討。

4 結(jié)論

1)在土壤調(diào)查樣本較為稀疏的大面積高寒山區(qū)，通過把豐富的成土環(huán)境因素變量與復(fù)雜的機器學(xué)習(xí)方法相結(jié)合，模型解釋度達(dá)83% 以上，能夠?qū)崿F(xiàn)較為準(zhǔn)確的土壤養(yǎng)分預(yù)測制圖，捕捉土壤空間變異。

2)青海省表層土壤養(yǎng)分在整體上呈現(xiàn)東南–西北均勻遞減趨勢。植被能有效地指示青海省表層土壤養(yǎng)分全氮(TN)、全磷(TP)、全鉀(TK)空間分布模式，其中年均降水量、地表溫度是影響青海省表層土壤TN空間模式的重要因素；地表覆被、海拔和地表溫度等環(huán)境因子對表層土壤TP的空間變異占主導(dǎo)作用；年均降水量、晝夜溫差對表層土壤TK的空間變異占主導(dǎo)作用。

3)土壤養(yǎng)分預(yù)測結(jié)果高等級出現(xiàn)在南部的玉樹、果洛、黃南和東部的湟水谷地地區(qū)，低等級主要分布在柴達(dá)木盆地、可可西里、海南州中南部，全省養(yǎng)分等級呈中上等級以上，占全省面積的81%。由此建立的土壤養(yǎng)分管理分區(qū)，可為全省土壤養(yǎng)分的監(jiān)測與農(nóng)牧業(yè)區(qū)劃的管理提供參考依據(jù)。

[1] 賈三滿, 杜麗娜, 馮輝, 等. 北京大清河流域生態(tài)涵養(yǎng)區(qū)生態(tài)地質(zhì)環(huán)境質(zhì)量綜合評價[J]. 北京地質(zhì), 2020, 15(2): 125–136.

[2] 李瑩瑩, 趙正勇, 楊旗. 數(shù)字土壤制圖在土壤養(yǎng)分方面的研究綜述[J], 江西農(nóng)業(yè)學(xué)報, 2021, 33(7): 61–67.

[3] Dai F Q, Zhou Q G, Lv Z Q, et al. Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau[J]. Ecological Indicators, 2014, 45: 184–194.

[4] 李德成, 趙霞, 張甘霖. 中國土系志·青海卷[M], 北京: 科學(xué)出版社, 2018.

[5] 楊煜岑, 楊聯(lián)安, 王晶, 等. 基于多元線性回歸模型的土壤養(yǎng)分空間預(yù)測——以陜西省藍(lán)田縣農(nóng)耕區(qū)為例[J]. 土壤通報, 2017, 48(5): 1102–1113.

[6] Dong W, Wu T J, Luo J C, et al. Land parcel-based digital soil mapping of soil nutrient properties in an alluvial-diluvia plain agricultural area in China[J]. Geoderma, 2019, 340: 234–248.

[7] Liu F, Rossiter D G, Song X D, et al. An approach for broad-scale predictive soil properties mapping in low-relief areas based on responses to solar radiation[J]. Soil Science Society of America Journal, 2020, 84(1): 144–162.

[8] 代子俊, 趙霞, 李德成, 等. 近30年湟水流域土壤全氮時空變異及影響因素[J]. 土壤學(xué)報, 2018, 55(2): 338–350.

[9] 張歡, 高小紅. 復(fù)雜地形區(qū)土壤有機質(zhì)空間變異性分析及制圖[J]. 水土保持研究, 2020, 27(5): 93–100.

[10] 周偉, 謝利娟, 楊晗, 等. 基于高光譜的三江源區(qū)土壤有機質(zhì)含量反演[J]. 土壤通報, 2021, 52(3): 564–574.

[11] 龐龍輝, 劉峰, 趙霞, 等. 青海省表層土壤屬性數(shù)字制圖[J]. 土壤通報, 2019, 50(3): 505–513.

[12] 陳榮, 韓浩武, 傅佩紅, 等. 基于多時相遙感影像和隨機森林算法的土壤制圖[J]. 土壤, 2021, 53(5): 1087–1094.

[13] Nussbaum M, Spiess K, Baltensweiler A, et al. Evaluation of digital soil mapping approaches with large sets of environmental covariates[J]. Soil, 2018, 4(1): 1–22.

[14] 朱阿興, 楊琳, 樊乃卿, 等. 數(shù)字土壤制圖研究綜述與展望[J]. 地理科學(xué)進(jìn)展, 2018, 37(1): 66–78.

[15] 青海省統(tǒng)計局. 青海統(tǒng)計年鑒.2002(總第18期)[M]. 北京: 中國統(tǒng)計出版社, 2002.

[16] 張甘霖, 龔子同. 土壤調(diào)查實驗室分析方法[M]. 北京: 科學(xué)出版社, 2012.

[17] Fick S E, Hijmans R J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas[J]. International Journal of Climatology, 2017, 37(12): 4302–4315.

[18] Leo Breimal. Random Forests[J]. Machine Learning, 2001, 45(1): 5–32.

[19] R DeveiopmentCore Team R A Language and Environmentfor Statistical Computing[M]. Vienna: R Foundation for Statistical Computing, 2014.

[20] Picard R R, Cook R D. Cross-validation of regression models[J]. Journal of the American Statical Association, 1984, 79(387): 575–583.

[21] 謝賢健, 韓光中. 基于普通克里格和投影尋蹤模型的城市土壤重金屬污染評價[J]. 生態(tài)環(huán)境學(xué)報, 2017, 26(9): 1584–1590.

[22] 全國土壤普查辦公室. 中國土壤[M]. 北京: 中國農(nóng)業(yè)出版社, 1998.

[23] 張鐵嬋, 常慶瑞, 劉京. 土壤養(yǎng)分元素空間分布不同插值方法研究——以榆林市榆陽區(qū)為例[J]. 干旱地區(qū)農(nóng)業(yè)研究, 2010, 28(2): 177–182.

[24] 石小華, 楊聯(lián)安, 張蕾. 土壤速效鉀養(yǎng)分含量空間插值方法比較研究[J]. 水土保持學(xué)報, 2006, 20(2): 68–72.

[25] Pouladi N, M?ller A B, Tabatabai S, et al. Mapping soil organic matter contents at field level with Cubist, Random Forest and kriging[J]. Geoderma, 2019, 342: 85–92.

[26] Wang S, Jin X X, Adhikari K, et al. Mapping total soil nitrogen from a site in northeastern China[J]. CATENA, 2018, 166: 134–146.

[27] 王俊雅, 劉峰, 宋效東, 等. 基于地表溫度的干旱平緩區(qū)土壤屬性制圖[J]. 土壤通報, 2018, 49(6): 1270–1278.

[28] 郭澎濤, 李茂芬, 羅微, 等. 基于多源環(huán)境變量和隨機森林的橡膠園土壤全氮含量預(yù)測[J]. 農(nóng)業(yè)工程學(xué)報, 2015, 31(5): 194–200, 202, 201.

[29] 孫孝林, 趙玉國, 趙量, 等. 應(yīng)用土壤-景觀定量模型預(yù)測土壤屬性空間分布及制圖[J]. 土壤, 2008, 40(5): 837–842.

[30] 楊琳, 朱阿興, 秦承志, 等. 運用模糊隸屬度進(jìn)行土壤屬性制圖的研究——以黑龍江鶴山農(nóng)場研究區(qū)為例[J]. 土壤學(xué)報, 2009, 46(1): 9–15.

[31] Brungard C W, Boettinger J L, Duniway M C, et al. Machine learning for predicting soil classes in three semi-arid landscapes[J]. Geoderma, 2015, 239/240: 68–83.

[32] 黃婷, 岳西杰, 葛璽祖, 等. 基于主成分分析的黃土溝壑區(qū)土壤肥力質(zhì)量評價——以長武縣耕地土壤為例[J]. 干旱地區(qū)農(nóng)業(yè)研究, 2010, 28(3): 141–147, 187.

[33] 朱昌達(dá), 高明秀, 王文倩, 等. 基于GIS的濱海鹽漬化農(nóng)田土壤空間變異及其分區(qū)管理[J]. 生態(tài)學(xué)報, 2020, 40(19): 6982–6990.

[34] 張彬, 楊聯(lián)安, 楊粉莉, 等. 基于投影尋蹤的土壤養(yǎng)分綜合評價及影響因素研究[J]. 土壤, 2020, 52(6): 1239–1247.

[35] 梁斌, 齊實. 北京山區(qū)土壤養(yǎng)分空間變化特征研究[J]. 土壤, 2018, 50(4): 769–777.

[36] Li H, Leng W J, Zhou Y B, et al. Evaluation models for soil nutrient based on support vector machine and artificial neural networks[J]. The Scientific World Journal, 2014, 2014: 478569.

Spatial Prediction and Management Zoning of Soil Nutrients in Large-scale Alpine Mountainous Areas

DU Longquan1,2,3, LIU Feng3,4, SHI Zhou5, ZHAO Xia1,2*, LI Decheng3, ZHANG Ganlin3,4,6, DONG Jinpeng1,2, CHEN Dongsheng1,2

(1 Qinghai Normal University, Qinghai Soil Digital Service Center, Physical Geography and Environmental Process Key Laboratory of Qinhai Province, Xining 810008, China; 2 Qinghai Normal University Academy of Plateau Science and Sustainability-Plateau Soil Information Science Research Team, People’s Government of Qinghai Province & Beijing Normal University, Xining 810008, China; 3 State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China; 4 University of Chinese Academy of Sciences, Beijing 100049, China;5 College of Environmental＆Resource Sciences, Zhejiang University, Hangzhou 310058, China; 6 Key Laboratory of Watershed Geography Sciences, Nanjing Institute of Geography & Limnology, Chinese Academy of Sciences, Nanjing 210008, China)

Based on the data of 205 sample points of soil series survey in Qinghai Province in recent years, the quantitative relationships between the contents of total nitrogen (TN), total potassium (TK) and total phosphorus (TP) of topsoils (0–20 cm) and environmental factor variables (terrain, climate, vegetation and remote sensing data) were established respectively by using random forest model, and the spatial distribution of soil nutrient contents in Qinghai Province was predicted, then the management zoning of soil nutrients was generated by using the projection pursuit method and the national soil nutrient classification standard. The cross validation results show that2of spatial prediction of TN, TK and TP are 0.89, 0.85 and 0.82, respectively. The model can explain more than 80% of the spatial variation of soil nutrients, indicating that the combination of random forest model and environmental factor variables can effectively predict the spatial variation of soil nutrients in large-area alpine mountainous areas under the condition of sparse samples. The distribution pattern of soil nutrients in Qinghai Province is high in the east and low in the west. The high levels of soil comprehensive nutrients appear in Yushu, Guoluo, Huangnan in the south and Huangshui Valley in the east; The lower grades are mainly distributed in Qaidam Basin, Hoh Xil and the south central part of Hainan prefecture; The soil nutrient classification of the whole province is above the middle and upper grades, accounting for 81% of the total area of the whole province. It is found that vegetation is the main environmental factor affecting the spatial distribution of soil nutrients in topsoil in Qinghai Province, among which annual precipitation and surface temperature are important factors affecting the spatial model of TN in topsoil in Qinghai Province. The spatial variation of TP in topsoil was dominated by environmental factors such as surface cover, altitude and surface temperature. Annual precipitation and temperature difference between day and night are important factors affecting the spatial model of TK in topsoil.

Soil nutrients; Random forest; Projection pursuit method; Management partition; Digital soil mapping

S159

A

10.13758/j.cnki.tr.2022.06.023

杜龍全, 劉峰, 史舟, 等. 大面積高寒山區(qū)土壤養(yǎng)分空間預(yù)測與管理分區(qū). 土壤, 2022, 54(6): 1273–1282.

國家自然科學(xué)基金重點項目(41930754)、國家自然科學(xué)基金地區(qū)項目(42067001)、青海省科技廳科技合作專項(2018-HZ-804)資助。

通訊作者(zhaoxia@qhnu.edu.cn)

杜龍全(1994—)，男，山西呂梁人，碩士研究生，主要從事資源環(huán)境與生態(tài)評估研究。E-mail: 15135132501@163.com