中國水蝕區土壤可蝕性因子更新方法與應用

田芷源, 梁 音?, 趙 院, 曹龍熹, 趙 艷,4, 武逸杭,4

(1.土壤與農業可持續發展國家重點實驗室,中國科學院南京土壤研究所,210008,南京;2.水利部水土保持監測中心,100053,北京;3.成都理工大學生態環境學院,610059,成都;4.中國科學院大學南京學院,211135,南京)

土壤可蝕性因子K是土壤流失方程的關鍵參數之一。K因子是侵蝕預報模型計算和水土流失動態監測的基礎數據。K值的定義是在標準小區上單位降雨侵蝕力所引起土壤流失量的多少[1],反映土壤對侵蝕的敏感程度,其大小與土壤本身性質有關。

我國第1次水利普查時獲得了全國范圍首張K值分布圖[2]。由于實測K值使用的徑流小區建造昂貴,花費時間長且缺少統一的測量規范,能收集到的K值數據十分有限,因此利用經驗公式對K值進行估算[3]。然而,使用的土種志數據距今已有近40 a,難以反映目前的土壤情況。為了獲得更準確的K值結果,亟需在大范圍內對計算K值的基礎數據進行更新。

依靠有限樣點很難滿足越來越高的制圖精度要求,因此還需探索一種能夠實現高空間分辨率的K值點面擴展方法。以往研究在利用樣點制作K值分布圖時,常使用土壤圖斑鏈接法及地統計插值[4]。然而,圖斑鏈接法導致K值在土壤類型邊界發生突變,且不能反映同一土壤圖斑內存在的K值變異性[5]。使用地統計學方法對K值進行制圖時,由于采樣密度對插值精度影響較大,需空間結構因素大于隨機部分才能得到較準確的結果,因此研究范圍大多局限在流域尺度[4]。隨著數字土壤制圖的發展,遙感數據已被應用于K值空間預測,例如高程作為協變量被引入克里格插值進行K值制圖[6],高程及植被作為空間變量可以提高K值的預測精度[7]。機器學習算法中的隨機森林模型,是一種通過自助重采樣技術對多個訓練集隨機選取特征變量,分別進行回歸建模并平均多個預測值得到最終結果的方法[8]。該方法已被成功運用于預測土壤有機質含量、質地與土層厚度等關鍵土壤屬性[9]。

為了實現全國土壤可蝕性因子更新目標,筆者提出一套基于近期匯編的中國土系志數據,利用隨機森林模型進行點面擴展獲得全國水蝕區K值分布圖的方法,解決更新的2個關鍵問題:第一是如何獲取全國地區最新的樣點K值,第二是如何將樣點的更新結果擴展至區域尺度。本次更新結果將有助于改善當前全國K值圖數據陳舊及制圖方法落后的問題。該成果可為土壤侵蝕預報模型提供基礎數據,為年度水土流失動態監測工作提供服務。

1 材料與方法

1.1 土壤理化數據來源與處理

中國土系志數據是我國自第2次土壤普查以來最新且系統性的土壤調查數據[10]。全國土系調查項目于2008—2018年基于定量標準和統一分類原則開展系統性調查研究[11]。筆者基于各省份土系調查資料,收集并整理4 327個采樣點的土壤剖面數據,查閱其調查地點(經緯度)和理化性質。其中有機碳質量分數通過除以系數0.58(采用Van Bemmelen因子,假定土壤有機質的含碳質量分數為58%)統一轉換為有機質質量分數[12]。將數據進行預處理,主要采用加權平均的方法提取各樣點0～30 cm耕層土壤的有機質和機械組成信息,當土體厚度<30 cm時則以實際厚度進行加權平均。

1.2 土壤可蝕性K值計算

采用全國第1次水利普查時K值的計算方法,使用通用土壤流失方程中的諾模圖(Nomo)計算K值(Nomo-K)[2]。

K=(2.1×10-4M1.14(12-O)+
3.25(S-2)+2.5(P-3))/(100×7.593);

(1)

M=N1(100-N2)

(2)

式中:K為土壤可蝕性,t·hm2·h/(MJ·mm·hm2);O為土壤有機質質量分數,%;S為土壤結構系數(計算土壤粒徑幾何平均直徑后查表所得,量綱為1);P為土壤滲透性等級(根據質地分類查表所得,量綱為1);M為質地指數,由N1和N2計算所得,%;N1為粒徑在0.002～0.100 mm之間的土壤顆粒質量分數比例,%;N2為粒徑< 0.002 mm的土壤顆粒質量分數,%。

由于全國土系調查中的質地分布僅包含砂粒(≥0.050～2.000 mm)、粉粒(≥0.002～0.050 mm)和黏粒(<0.002 mm)分級,Nomo公式計算時需要把已有土壤顆粒分析結果轉化為極細砂質量分數(≥0.050～0.100 mm)。本研究比較了3種插值方法,分別是自然對數線性插值法、三次樣條函數法、三次樣條函數結合自然對數。計算過程如下:首先將各級別土壤顆粒質量分數累加計算得到小于某一粒徑的質量分數(<0.002 mm、<0.050 mm、<2.000 mm),并計算該粒徑對應的自然對數(ln(0.002)、ln(0.050)和ln(2.000))。自然對數線性插值法可建立顆粒累加質量分數與粒級自然對數之間的線性關系;三次樣條函數法可建立顆粒累加質量分數與粒徑之間的三次樣條函數;三次樣條函數結合自然對數法則可建立顆粒累加質量分數與粒徑自然對數之間的三次樣條函數?；谝陨戏椒ǖ玫搅?0.100 mm粒徑所對應的質量分數,并將<0.050 mm和<0.100 mm粒徑質量分數相減得到極細砂質量分數,通過分析極細砂質量分數的合理性和比較插值方程R2確定最優的插值方法。

當有機質質量分數>12%時,使用Nomo公式可能會導致負值出現,因此對有機質>12%的樣點采用EPIC模型計算K值[3]

(3)

Sn1=1-Sa/100。

(4)

式中:Sa為粒徑≥0.050～2.000 mm的砂粒質量分數,%;Si為≥0.002～0.050 mm 的土壤粉粒質量分數,%;Cl為粒徑< 0.002 mm的土壤黏粒質量分數,%;C為土壤有機碳質量分數,%,通過有機質乘以轉換系數58%所得;Sn1為土壤碳酸鈣質量分數,%。

由于Nomo和EPIC公式的計算過程不同,采用EPIC公式得到的K值(EPIC-K)不能直接替代 Nomo-K值。因此,筆者使用線性、指數、冪函數和3階多項式函數分別建立2種公式計算K值的轉換關系。比較4種擬合方程的R2,除線性方程擬合精度略低以外(0.795 5),其他3種非線性方程的擬合精度十分接近(0.806 4～0.807 2)。從趨勢來看,Nomo-K和EPIC-K值之間并非簡單線性關系,隨著EPIC-K值的增大,Nomo-K值出現先緩慢增加后快速增加的趨勢。這與指數方程所顯示的規律是一致的?？傊?當土壤有機質≤12%時,使用Nomo公式計算K值;當有機質>12%時,采用EPIC公式計算K值并利用其與Nomo-K值的指數關系式進行修正。

1.3 環境因素分布圖收集與處理

沿用在太湖流域片預測K值時使用的遙感方法獲取的環境因子指標作為預測變量,包括氣候、地表溫度、植被、光波、地形和母質[13]。由于隨機森林模型具有避免過擬合的優勢,且模型訓練集樣本量足夠大,沒有對環境因子進行剔除。

各項環境指標的空間分辨率從30 m到1 km不等,通過最鄰近法將所有指標統一重采樣為30 m分辨率后進行疊加運算。由于氣候和地表溫度在較小尺度下不存在大的差異,重采樣后在1 km范圍內氣候及地表溫度值不變符合實際情況。各樣點對應的環境變量值通過坐標位置提取。

1.4 土壤可蝕性K值制圖

采用隨機森林回歸方法建立K值預測模型。首先對樣點的環境要素及K值進行訓練,然后利用多源遙感影像將建立起K值預測模型推廣至全國地區,完成K值由點及面的反演計算。其中模型的2個關鍵參數通過調參選擇:模型誤差隨著決策樹數量(ntree)的增加而降低,當決策樹數量超過400以后,模型誤差趨于穩定,因此選擇使用默認參數ntree=500;隨機森林節點變量數(mtry)一般選擇環境變量數的1/3(本研究對應mtry=20),通過對節點變量數在mtry=20左右逐個比較模型精度,確定最佳的mtry值為17。

K值預測模型的建模精度通過外部驗證的方法進行評估,采取常用的5折交叉驗證法,重復100次取平均值作為評價結果。評價指標包括決定系數 (R2),均方根誤差 (E1)以及平均絕對誤差 (E2),計算公式如下:

(5)

(6)

(7)

式中:n為訓練樣本個數;Yi為觀測K值;i為模型預測K值;為平均觀測K值。

基于K值預測模型,以61項環境因子空間分布圖作為自變量,使用隨機森林方法進行預測,輸出全國尺度K因子柵格圖,K值更新圖的空間分辨率為30 m×30 m。

1.5 數據處理與統計分析

研究利用Excel 2019計算樣點Nomo-K值與EPIC-K值,并建立兩者的回歸關系;利用ArcGIS 10.8處理空間數據,包括投影轉換、空間重采樣,提取樣點對應的環境變量,統計柵格極值、平均值及標準差等;利用R語言編程實現土壤粒徑的批量插值轉換,同時構建隨機森林回歸模型,進行全國范圍K值的空間預測。在分析K值分布規律時,土壤類型圖使用1995年編制的《1∶100萬中華人民共和國土壤圖》。

2 結果與分析

2.1 土壤質地插值轉換

插值獲得諾模圖計算所需極細砂的質量分數,以江西省關山系樣點為例,使用3種插值方法得到粒徑質量分數結果如表1所示。其中,自然對數線性插值法擬合精度較差,獲得的<0.100 mm粒徑質量分數低于<0.050 mm粒徑,導致極細砂質量分數為負值;三次樣條函數法受樣本數量過少(n=3)的影響,獲得的累加極細砂粒徑(<0.100 mm)質量分數高于100%,明顯不合理;三次樣條函數結合自然對數法克服了函數分布不合理的缺點,以高精度方程擬合了不同土壤粒徑的質量分數(R2≈1),同時避免了質量分數高于100%的情況(圖1)。

圓點為已知粒徑質量分數,三角形為預測粒徑質量分數。下同。The dot represents the known particle size mass fraction, and the triangle represents the predicted particle size mass fraction. The same below. 圖1 土壤質地插值轉換方法對比Fig.1 Comparison of soil texture interpolation conversion methods

表1 土壤質地插值轉換結果對比(江西關山系)Tab.1 Comparison of soil texture interpolation conversion results(Guanshan series in Jiangxi)

在少數樣點中,土壤<0.050 mm和<2.000 mm的粒徑質量分數十分接近,使用三次樣條函數結合自然對數法計算所得<0.100 mm的粒徑結果可能略高于100%。針對這種情況,以山東省小島河系樣點為例(圖2),利用線性函數在<0.050 mm和<2.000 mm分段之間插值得到<0.100 mm的粒徑質量分數,線性插值使<0.100 mm粒徑質量分數既小于<2.000 mm粒徑質量分數,也大于<0.050 mm粒徑質量分數,能確保結果的合理性。

實線為三次樣條函數結合自然對數法擬合,虛線為自然對數線性插值法。The solid line is cubic spline function combined with natural logarithm, and the dotted line is linear interpolation with natural logarithm. 圖2 土壤粒徑插值結果過高處理Fig.2 Processing of excessive soil particle size interpolation results

2.2 土壤可蝕性K值的計算與預測模型

采用Nomo和EPIC公式分別計算4 327個土系樣點K值,統計結果如表2所示。EPIC計算K因子的最小值比Nomo-K值大,最大值比Nomo-K值小,且標準差更小,EPIC計算K值的平均值及中位數比Nomo-K值大16%。因此,不能直接使用EPIC-K值替換高土壤有機質條件下的Nomo-K值。

表2 全國土壤可蝕性K值計算結果統計Tab.2 Statistics of calculated results of soil erodibility K value nationwide

建立2種公式的轉換關系時,R2均>0.79(圖3),但散點圖顯示隨著EPIC-K值的增大,Nomo-K值出現先緩慢增加后快速增加的趨勢,這與指數方程的規律相符。最終選擇該方程對全國有機質>12%的55個樣點進行K值修正。這些樣點主要分布在東北和青藏高原的濕潤或者滯水條件下,采樣部位在上、中、下坡和臺地均有分布,表明并非沉積導致K出現負值。根據定義,K因子為在標準小區上由單位降雨侵蝕力所引起土壤流失量的多少[1]。土壤流失量不可能為負數,所以K值不應存在負值的情況。校正前55個樣點的EPIC-K均值為0.034 4,校正后均值為0.026 6,彌補了EPIC-K值大于Nomo-K值的差異。

EPIC-K is soil erodibility calculated by EPIC model. Nomo-K is soil erodibility calculated by Nomo formula.圖3 土壤可蝕性K值計算結果擬合方程對比Fig.3 Comparison of fitting equations for the calculated results of soil erodibility K

以土系調查項目所建立的更新K值作為因變量,以包含氣候、地形、植被、母質類型、地表溫度、可見光及近紅外波段在內的61項環境因子指標作為自變量,在mtry=17,ntree=500的模型參數下,建立的全國K因子隨機森林預測模型,運用5折交叉驗證方法重復100次的評估結果為:R2=0.381,E1=0.012 47,E2=0.009 63。這與太湖流域片建立K值模型的精度類似[13],但低于使用隨機森林預測土壤有機質和全氮的空間分布[14]。這與本研究采樣密度稀疏、空間建模尺度大以及缺乏與K值強相關性的預測變量等原因有關。

2.3 全國土壤可蝕性K值分布規律

全國K因子分布范圍在0.005 1～0.074 5之間,這與前人在樣點尺度建立的我國主要土壤K因子分布范圍0.000 8～0.070 5相似[15]。分布趨勢呈現出華北平原和黃土高原地區高,天山地區及東北平原次之,南方山地丘陵、青藏高原及內蒙古高原較低的規律。這與其他學者總結的我國土壤可蝕性分布規律相似[16],但該研究發現青藏高原的K值最高。也有研究發現青藏高原K值比黃土高原、川渝地區和滇東北地區的K值略低[17],這可能是由于青藏高原K值存在較大差異所造成的。

進一步利用水土保持一級區劃分析全國更新K值的分布規律,發現更新K值在西北黃土高原區和北方土石山區較大,在青藏高原區最小(表3),這與自然及人為活動等因素的綜合作用有關。其中,西北黃土高原區土壤有機質含量普遍很低因此K值偏高[18];北方土石山區長期以來受人類耕作活動的影響導致K值較高;青藏高原區由于地形落差大,K值具有垂直地帶性,在少量低海拔河谷地帶K值較高,而在大面積的高海拔地區由于成土作用弱K值明顯降低[17]。此外,東北黑土區在風砂土分布地區K值最大,在黑土集中的漫川漫崗地帶由于農耕活動導致K值較大,暗棕壤分布的大興安嶺及長白山等地K值較小。

表3 土壤可蝕性更新K值在水土保持一級區劃的統計Tab.3 Statistics of soil erodibility update K value in the first level zoning of soil and water conservation

K值的空間分異與我國主要土壤類型的分布情況有關,前人曾報道黃土高原的黃綿土和東北的黑土K值較大,南方紅壤K值較小[15,19]。本研究從主要土壤類型統計來看,黃綿土和潮土的更新K值較大,平均值均超過0.04,而高山土的平均K值最小,其次為栗鈣土、暗棕壤和紅壤,均不超過0.03(表4)。分析原因認為黃綿土和潮土的粉粒含量較高,土壤更容易遭受侵蝕。此外,潮土K值的標準差最大。前人研究發現,同一種土壤下不同母質和植被覆蓋可以通過影響土壤的水穩性從而改變土壤的抗蝕能力[20],因此相同土壤類型下K值的差別也可能很大。

表4 主要土壤類型下土壤可蝕性更新K值統計Tab.4 Statistics of updated K values of soil erodibility under main soil types

值得注意的是,由于Nomo和EPIC模型是基于美國土壤實驗所建立的,在中國地區存在適用性問題。張科利等[19]利用徑流小區觀測資料在我國東部地區建立了K值的校正公式,分別為K=-0.033 36+0.744 88KNomo(R=0.721),K=-0.013 83+0.515 75KEPIC(R=0.613);更進一步說,Nomo模型在黃土高原和黑土區的應用效果(R2=0.55～0.60)要優于紅壤和紫色土(R2=0.11～0.37),EPIC模型除黑土區以外(R2=0.05)在其他地區的應用效果較好(R2>0.58)[21]。直接應用本研究結果可能會對K值產生不同程度的高估,因此需要在全國范圍內進一步收集實測資料將K值校正以后用于土壤侵蝕模數計算。

3 結論

1)本次更新進一步完善了全國土壤可蝕性因子計算方法。本研究在美國制砂粒、粉粒、黏粒機械組成實測值的基礎上插值得到極細砂質量分數,基礎數據更為準確;通過計算土壤粒徑幾何平均直徑來獲取結構系數,更具有客觀性;建立于全國統一的土系基礎數據庫,消除了不同省份之間的系統誤差;以30 m分辨率柵格為制圖單位,更為詳細地表達K值的空間變異情況。

2)全國土壤可蝕性更新K值的空間分布呈現出以下規律:黃土高原區和北方土石山區居高,西南紫色土區、南方紅壤區、西南巖溶區、東北黑土區和北方風沙區次之,青藏高原區最小。這主要是由于不同土壤類型在各地區的分布所造成的。統計可知黃綿土和潮土的K值較大,紅壤、暗棕壤、栗鈣土的K值較小,高山土的K值最小。

3)本研究采用的Nomo和EPIC公式都是基于美國土壤建立的K值計算經驗公式,可能與實際監測的K真實值有差距。因此需要積累我國不同地區標準徑流小區的長期觀測資料,使用年平均侵蝕量和降雨侵蝕力因子進一步率定K值,以期建立起實測值與經驗公式計算值之間的關系式或直接建立適用于我國的土壤可蝕性計算公式。