雅魯藏布江貢嘎段流域土壤可蝕性的不同估算方法對比研究

趙萌萌 張 根 古明雙 李云飛

（西藏大學工學院，西藏 拉薩 850000）

0 引言

土壤可蝕性能夠反映土壤對侵蝕的敏感性，是衡量土壤侵蝕程度的重要指標，在通用土壤流失方程（USLE）中用K來表示。對土壤可蝕性進行深入研究是有效治理土壤侵蝕的必要環節。我國的土壤可蝕性研究起步于20世紀50年代。幾十年，來我國學者十分重視土壤可蝕性研究，并將其作為水土保持研究的重要內容之一［1］。

由于經典的土壤可蝕性估算方法——自然徑流小區內直接測定法所耗費的時間較長，數據收集較為困難［2］，因而國內外學者開始嘗試建立更為精確的數學模型用以計算土壤可蝕性K值，如諾謨圖法、修正諾謨圖法、Shirazi公式、EPIC模型和Torri.D模型等。其中，EPIC模型是目前國內外學者估算土壤可蝕性K值普遍采用的模型之一［3］；Torri.D模型、Shirazi公式多應用于不同估算模型的比較研究中，而且將其應用于特定地區進行土壤可蝕性K值估算時，其準確性比EPIC模型高［4］。

我國學者利用上述經驗模型對不同地區的土壤可蝕性K值進行了測算，但大多數學者并未對所得結果進行不確定性評價，導致土壤侵蝕預報的可信區間不夠精確［5］。基于此，筆者以雅魯藏布江貢嘎段流域為研究對象，先測定野外采集土壤樣品的相關參數，再利用諾謨圖法、修正諾謨圖法、Shirazi公式、EPIC模型和Torri.D模型估算雅魯藏布江貢嘎段流域土壤可蝕性K值，而后對結果進行數值對比、正態性檢驗及相關性分析，最后選用平均絕對誤差、均方根誤差、平均相對誤差及精度因子4個因子對5種估算方法所得K值進行不確定性評價，為研究雅魯藏布江貢嘎段流域土壤可蝕性提供方法支撐。

1 試驗材料與方法

1.1 研究區概況

雅魯藏布江貢嘎段流域位于西藏自治區貢嘎縣北部，海拔3 500～3 650 m，全長約81 km，地勢上大體為西高東低。雅魯藏布江貢嘎段兩岸岸坡較陡，風化作用強烈，植被覆蓋率較低。由于自然條件惡劣，當地荒漠化區域逐年擴大。同時，長期以來的過度放牧、大量垃圾傾倒等人為活動，造成該流域草場退化、環境污染嚴重、植被覆蓋率下降，生態環境極其脆弱［6］。另外，當地泥石流、滑坡等自然災害頻發，致使大量泥沙被沖刷進入雅魯藏布江［7］。上述因素導致該流域土壤侵蝕嚴重。近年來，雅魯藏布江中游河段的植被覆蓋率因人工林的栽種而有所提高，建成的生態防護林體系已初具規模［8］。

1.2 土壤樣品處理及K值估算

此研究始于嘎雜村附近（東經91°13′42.79″、北 緯29 ° 16′18.18″，海 拔3 566.40 m），終于 拉 玉 附 近（東 經90 ° 29′35.99″、北 緯29°14′57.93″，海拔3 653.41 m），以8 km為間隔標記10個取樣點，根據現行《水土保持試驗規程》（SL 419—2008）［9］中規定的樣品采集方法及過程，在每個取樣點取深度為20 cm的土壤樣品1 kg，共得10組土壤樣品。該研究采用的K值估算方法為諾謨圖法、修正諾謨圖法、Shirazi公式、EPIC模型、Torri.D模型。利用上述5種方法進行估算時，需要測定的土壤參數有土壤顆粒組成、有機質含量、有機碳含量（有機質含量除以1.724）、土壤結構參數及滲透性級別參數的取值等級［10-11］。各取樣點土壤顆粒組成及有機質含量數據見表1。

表1 各取樣點土壤理化性質

此外，在估算土壤可蝕性K值的過程中，為確保結果的準確性，需要確定樣品的土壤結構參數及土壤滲透性級別參數的取值等級。土壤結構參數（Ss）取值等級標準見表2［12］。表2中，土壤團聚體的結構類型根據其粒徑來確定。在該研究中，利用干篩法對所取樣品的土壤團聚體粒徑進行測定，得到樣品的土壤團聚體粒徑為2～5 mm。根據各個結構類型的粒徑劃分［13］，所取樣品的土壤團聚體結構類型為中等團粒，故Ss取3。

表2 土壤結構參數取值等級

根據現行的《水利水電工程地質勘察規范》（GB 50487—2008）［14］，土壤的滲透性級別參數（Pr）取值等級見表3。

表3 土壤滲透性級別參數取值等級

根據已有的研究成果，西藏自治區的土壤滲透系數為7.82×10-6～0.93×10-3cm/s［15］。結合樣品中各土壤類型的滲透系數經驗值［16］，得到樣品的滲透等級為慢速至中等，故該研究中Pr取4。

通過對所取的10組土壤樣品進行初步觀察及送檢分析，參照已有的西藏自治區各地土壤分析成果［17-19］，結合《西藏土種志——基于全國第二次土壤普查的數據集》［20］，確定所取土壤樣品的土壤類型及土地利用方式。綜合上述土壤參數，分別采用5種模型估算得到的10組K值（每組包含K諾謨、K修正、KEPIC、KShirazi、KTorri.D5個數據）。由于只有Torri.D模型的K值單位是國際制，其余4種模型的單位均為美國制，因此，進行結果分析之前將K值的國際制單位［t·hm2·h/（MJ·mm·hm2）］除以0.131 7，即變為美國制K值單位。

2 試驗結果與分析

2.1 不同估算方法下K值特征對比

采用SPSS 17.0軟件對土壤參數及5種方法所得的K值進行正態性檢驗，結果如表4所示。

表4 土壤參數及K值的統計分析

表4數據表明，各土壤粒徑等級的含量多少為砂粒＞粉粒＞黏粒＞極細砂粒。各模型所得K值平均值為KEPIC＞KShirazi＞KTorri.D＞K諾謨＞K修正。標準差是方差的算術平方根，其能反映數據的離散程度。5種模型所得K值的標準差中，KTorri.D最大，K諾謨最小，可見KTorri.D的數據離散性最強，K諾謨最弱。變異系數為樣本數據的標準差除以平均值所得的數值，其劃分等級為弱變異性（＜10%）、中等變異性（10%～100%）、強變異性（＞100%）。表4中，各土壤粒徑等級及不同模型所得K值的變異系數均處于10%～100%，因此，各粒級及不同模型估算的K值屬中等變異性。此外，各個土壤樣品的有機質含量的變異性也較強，且具有較強的正態性。由非參數K-S檢驗的結果可見，K值符合正態分布，且KShirazi的正態性顯著較強，KTorri.D的正態性相對較弱。

利用Excel 2019軟件繪制出在各個取樣點利用不同估算方法所得的K值簇狀柱形圖，如圖1所示。

圖1 各模型所得K值

由圖1可知，各模型所得的K值差異較大，最小值為0.12，最大值為0.77。取樣點4為農用地，但其K值明顯高于同種土壤類型的非農用地取樣點1。其原因可能是不合理的耕種造成土壤養分結構發生變化，使土壤質量下降。此外，利用修正諾謨圖法所得的K值顯著低于利用另外4種模型所得K值，可能是因為所取土壤的砂粒和黏粒含量較高，致使該方法的結果不準確。

2.2 各模型所得K值與土壤相關參數的關系分析

利用SPSS 17.0對各個土壤參數與不同模型所得的K值進行Pearson雙變量相關分析，結果如表5所示。

表5 各模型所得K值與土壤參數的相關性分析

由表5可知，K修正、KEPIC、KShirazi均與砂粒和有機質含量呈負相關；K諾謨與黏粒含量呈顯著正相關，但與有機質含量的相關性不顯著。總體來看，各模型所得的K值均與黏粒、砂粒及極細砂粒含量的相關性較強，與有機質及粉粒含量的相關性較弱。

2.3 K值的不確定性分析

利用數學模型估算K值后，對所得結果進行不確定分析是必要的。該研究選用平均絕對誤差（MAE）、均方根誤差（RMSE）、平均相對誤差（MRE）及精度因子（Af）4個因子對所得K值進行不確定性評價。4個因子的計算式為

式（1）至式（4）中：NKi為已有研究所得的第i種土壤類型的理論K值［21］；SKi為不同模型估算所得的第i種土壤類型的實際K值；n為取樣點數量，該研究中n=10。根據統計學原理，MAE、MRE和RMSE的數值越趨近于0，Af的數值越接近于1，則K值預測的不確定性越小。對5種模型估算結果的不確定性評價結果見表6。

表6 各模型K值的不確定性評價

由表6可知，修正諾謨圖法的4個因子數值均為最大，而諾謨圖法的MAE、MRE和RMSE因子值相對較小，其Af值相比另外4個模型的Af值更趨近于1。對4種因子的數值進行綜合排序，得出5種估算模型對雅魯藏布江貢嘎段流域K值估算的不確定性從小到大依次為諾謨圖法、Torri.D模型、EPIC模型、Shirazi公式、修正諾謨圖法，即諾謨圖法對雅魯藏布江貢嘎段流域K值估算的適用性最好，修正諾謨圖法適用性最差。

3 結論

①雅魯藏布江貢嘎段流域的土壤粒級含量從大到小為砂粒＞粉粒＞黏粒＞極細砂粒，各個取樣點所測得的土壤有機質含量差異較大。

②各粒級及不同模型估算的K值屬中等變異性，極細砂粒的變異系數最高，砂粒最低；Torri.D模型變異系數最高，EPIC模型最低。

③各模型所得K值平均值為KEPIC＞KShirazi＞KTorri.D＞K諾謨＞K修正，KShirazi的數據離散性最強，K諾謨最弱。

④各模型所得K值均與黏粒、砂粒和極細砂粒的相關性較強，與有機質及粉粒的相關性較弱。

⑤K值不確定性分析表明，5種估算模型對雅魯藏布江貢嘎段流域K值估算的不確定性從小到大為諾謨圖法、Torri.D模型、EPIC模型、Shirazi公式、修正諾謨圖法。在該流域的K值估算中，諾謨圖法適用性優于另外4種模型，即可使用該模型在雅魯藏布江貢嘎段流域進行K值預測和水土流失敏感度評價。