基于RUSLE的土壤侵蝕建模分析

1、空間信息應用實踐（中級）實驗指導書空間建模基于RUSLE的土壤侵蝕建模分析一實驗背景Soil erosion and gullying in the upper Panuco basin, Sierra Madre Oriental, eastern Mexico土壤侵蝕是地球表面物質運動的一種自然現象，全球除永凍地區(qū)外，均發(fā)生不同程度的土壤侵蝕。人類社會出現后，土壤侵蝕成為自然和人為活動共同作用下的一種動態(tài)過程，構成了特殊的侵蝕環(huán)境背景，并伴隨著人類對自然改造能力的增強，逐漸成為當今世界資源和環(huán)境可持續(xù)發(fā)展所面臨的重要問題之一。土壤侵蝕被稱為“蠕動的災難”，每年因土壤侵蝕造成的經濟損失較諸如

2、滑坡、泥石流和地震等地質災害更大, 土壤侵蝕已成為我國乃至全球的重大環(huán)境問題之一。土壤侵蝕及其產生的泥沙使土壤養(yǎng)分流失、土地生產力下降、湖泊淤積、江河堵塞，并造成諸如洪水等自然災害，泥沙攜帶的大量營養(yǎng)物和污染物質加劇了水體富營養(yǎng)化，水質惡化，不斷嚴重威脅到人類的生存。據估計全球每年因土壤侵蝕損失300萬公頃土地的生產力，造成的損失以百億美元計。我國人口眾多、農耕歷史悠久，加之歷史上戰(zhàn)亂頻仍，以黃土高原為代表的華夏文明發(fā)源地是世界上土壤侵蝕最嚴重的區(qū)域之一， 1990年遙感普查結果，全國水土流失面積達367萬km2，占國土總面積的38.2，其中50%為水蝕地區(qū)，土壤侵蝕以黃土高原、四川紫色土地區(qū)

3、和華南紅壤地區(qū)尤為突出，僅黃土高原地區(qū)一處，平均每年流失泥沙就達到16.3 億t。水土流失已成為中國重要的環(huán)境問題，土壤侵蝕研究已成為目前環(huán)境保護中的一個重要課題。土壤侵蝕預報是有效監(jiān)測水土流失和評價水保措施效益的手段，侵蝕模型則是進行土壤流失監(jiān)測和預報的重要工具。然而傳統(tǒng)預測方法需要在量經費、時間和人力的投入，因此，在一定精度范圍內通過有限的數據輸入，得到滿足要求的土壤侵蝕預測結果成為趨勢。80年代以來，隨著地理信息系統(tǒng) (Geographical Information System, GIS)的成熟，它開始與土壤侵蝕模型通用土壤流失方程 (Universal Soil Loss Equa

4、tion, USLE) 相結合進行流域土壤侵蝕量的預測和估算，業(yè)已成為土壤侵蝕動態(tài)研究的有力工具。GIS與USLE 相結合的分布式方法運用GIS的柵格數據分析功能，可預測出每個柵格的土壤侵蝕量，便于管理者識別關鍵源區(qū)，并通過確定引起水土流失的關鍵因子，針對性地提出最佳管理措施 (Best Management Practices，BMPs)，為流域內土地資源的質量評價、利用規(guī)劃和經營管理等提供科學依據與決策手段。二、實驗目的模型生成器 (ModelBuilder) 為設計和實現空間處理模型提供了一個圖形化的建模環(huán)境。模型是以流程圖的形式表示，它通過工具將數據串起來以創(chuàng)建高級的功能和流程。你可以

5、將工具和數據集拖動到一個模型中，然后按照有序的步驟把它們連接起來以實現復雜的 GIS 任務。通過對本次練習達到以下目的： 掌握如何在ModelBuilder環(huán)境下通過繪制數據處理流程圖的方式實現空間分析過程的自動化； 掌握土壤侵蝕理論的基本知識； 掌握利用腳本文件實現空間建模，加深對地理建模過程的認識，對各種GIS分析工具的用途有深入的理解； 在ModelBuilder環(huán)境下如何計算RUSLE模型的中各個因子，實現RUSLE模型自動化；三、實驗準備實驗環(huán)境：ArcGIS Desktop 9.3實驗數據：矢量和柵格數據矢量數據：研究區(qū)界線(bj.shp)、氣象數據（Climate.shp），土地

6、利用數據(landuse_Clip.shp,)和土壤數據(soil_clip)；柵格數據：地形數據（DEM）； 四、實驗內容與步驟(1)實驗準備本次試驗需要使用ArcGIS的建模功能，在實驗之前需要掌握如何利用ArcGIS進行建模。首先，打開ArcMap，激活工具箱在工具箱中右鍵單擊，選擇“New Toolbox ”，即可新建一個工具箱?？梢栽诖斯ぞ呦渖嫌覔?，通過“Rename”對工具箱重命名。 在新建的工具箱上右擊，按照“New”“Model”新建一個Model，可以按照同樣的方法給這個Model命名。然后在此Model上右擊，通過“Edit”進入模型的編輯模式。到此，模型準備已經結束，接下

7、來開始逐個建立模型的各個因子。（2）地形因子（L，S因子）算法：坡長因子采用公式計算, ，式中：L為坡長因子，l為像元坡長，m為坡長指數，像元坡長的計算式如下：，m取值如下式： 式中，為像元坡度 (%)式中，li為像元坡長，Di為沿徑流方向每像元坡長的水平投影距 (在柵格圖像中為兩相鄰像元中心距，隨方向而異)，i為每個像元的坡度 ()，i為自山脊像元至待求像元個數。坡度因子S分段計算：L和S因子的模型建立：首先在工具箱中找到Resample工具，可以使用工具箱自帶的搜索功能快速定位到。在工具箱的下方有一行標簽，選擇Search標簽，在搜索框中輸入要查找的工具名，如Resample，點擊Sear

8、ch進行查詢，查詢結束后選中查詢結果，點擊下方的Locate可以快速定位需要查找的工具。可以將這個工具直接拖到Model的編輯窗口中，如圖：現在需要給這個工具添加一個參數，在編輯窗口的Resample上右擊，通過“Make Variable”“From Parameter”“Input Raster”添加。注意：這里不建議使用右鍵菜單的“Create Variable”來添加輸入輸出參數，因為很多工具拖入到編輯窗口后會自帶一個輸出參數，而且它們也有自己的默認輸入參數。如果另外新建一個參數，可能會因為這個新建參數類型不與工具要求的輸入參數類型對應而出現錯誤。按照同樣的方法拖入Slope工具，Si

9、ngle Output Map Algebra工具。通過編輯窗口上的工具將這些工具首尾連接起來。雙擊Input Raster，輸入dem數據輸入數據之后，編輯窗口中的工具顏色會相應的變化，說明這些工具已經相互連接起來，還是白色的工具代表它還沒有和前面的工具聯(lián)系起來構成“流水線”，同時，這也是判斷Single Output Map Algebra工具中的腳本語言是否和前面的輸出文件關聯(lián)起來的依據。在相互連接的工具中，只要有一個工具是白色的，就說明這條“流水線”不能正常運行?？梢园l(fā)現Single Output Map Algebra工具還是白的，這是因為我們沒有添加算法，下面添加用于計算S因子的算

10、法，依據為：坡度因子S分段計算：雙擊Single Output Map Algebra工具，添加如下代碼：Con(Slope_degree1 = 5 & Slope_degree1 5）注意：之所以使用Slope_degree1 * 3.14 / 180這是將原來的角度制轉化為弧度制，計算機不能識別角度制。注意：在運算符（如+，-，*，/）的左右要有空格，如Slope_degree1 * 3.14 / 180不要寫成Slope_degree1*3.14/180。前者的運算符左右有空格，后者的運算符左右沒有空格。這一點必須嚴格遵守，否則相同的代碼會出現不同的錯誤。這樣一來，S因子的模型就建立好了

11、。點擊運行運行成功口在S上右擊，選擇“Add To Display”就可以將結果顯示出來結果為：接下來對L因子建立模型L因子模型的建立可以在上面的S因子模型基礎上進行。需要添加工具Fill和工具Flow Direction以及Single Output Map Algebra工具，然后將他們連接起來Single Output Map Algebra工具中的代碼為：Con(FlowDirection = 2 | FlowDirection = 8 |FlowDirection = 32 | FlowDirection = 128 , Sqrt(2) * 90 , 1 * 90)這些步驟和在建立S

12、因子模型的時候是一樣的，這里不再贅述。下面解釋一下代碼：3264128161842對于中間柵格來說，它與鄰近的8個柵格中心點之間的距離只有兩種（sqrt(2)和1），當流向為32,128,8,2的時候，距離為sqrt(2)，其他情況下距離為1。因此在代碼中的表現就是如上的形式。接下來建立中計算m的模型這個過程有點和計算S因子類似，不同的是計算S因子使用的是度，而計算m使用的是百分比。依據的公式為： 式中，為像元坡度 (%)將Slope工具和Single Output Map Algebra工具工具拖入編輯窗口中，然后將他們連接起來。雙擊Slope(2)，進行如下設置：這里要選擇Persent百

13、分比的形式，因為m的判斷是利用百分比形式的坡度。然后在Single Output Map Algebra中輸入的代碼為：CON(Slope_p = 0.05,0.5,CON(Slope_p = 0.03,0.4,CON(Slope_p = 0.01 & Slope_p = 5。因為在變量存儲值的時候，是不會存儲%的，它只能利用浮點型數字來表示百分數。接下來計算中的l，依據式中，li為像元坡長，Di為沿徑流方向每像元坡長的水平投影距 (在柵格圖像中為兩相鄰像元中心距，隨方向而異)，i為每個像元的坡度 ()，i為自山脊像元至待求像元個數。從公式中可以看出，我們需要求出D和才能計算l，而在前面的部分

14、，我們已經算出了D，也已知了，在算S因子的時候就已經算出了，因此我們可以再次利用Slope計算后的結果來計算l。添加一個Single Output Map Algebra工具，將其與Slope_degree1和D相連。在Single Output Map Algebra中輸入如下代碼：D / Cos(Slope_degree1 * 3.14 / 180)最后就是L因子的計算了再添加一個Single Output Map Algebra工具，將其與L1和m相連，并輸入代碼：Pow(L1 / 22.13,m)最終結果模型為運行這個模型，并將L通過“Add To Display”添加到視圖中。到此，

15、L和S因子完成。(3)降雨侵蝕力因子（R因子）利用日降雨量估算降雨侵蝕力的多參數模型來計算流域的降雨侵蝕力，公式如下：R = - 0.0334 P + 0.006661 P2 (1)式中R表示的侵蝕力值 ( MJmmhm-2 h-1)，P表示年雨量（mm）。為了方便，我們新建一個Model，命名為“K因子”。在編輯窗口中先后拖入Create Thiessen Polygons工具，Feature to Raster 工具，還有Single Output Map Algebra工具，并把它們連接起來雙擊Input Features輸入實驗數據。雙擊Create Thiessen Polygons

16、進行如下設置雙擊Feature to Raster進行如下設置在Single Output Map Algebra工具中輸入如下代碼：（不要忘記修改參數Output_raster）0.0066611 * Output_raster * Output_raster - 0.0334 * Output_raster最終的模型為：運行這個模型，并將最終的R結果添加到視圖中到此，R因子完成。（4）土壤侵蝕力因子（K因子）新建一個Model，并且命名為“K因子”，如圖K因子反映了土壤對侵蝕的敏感性。影響K因子的因素是多方面，一般說來，質地越粗或越細的土壤有較低K值，而質地適中的反而有較高的K值。K值估算

17、采用Williams等在EP IC模型中的方法, 利用土壤有機質和顆粒組成因子進行估算，計算式如下：K =0.2+0.3exp-0.0256Sd(1 - Si/100) * Si/(Cl+ Si)0.3*1.0-0.25*(C/1.724)*(C/1.724)+exp(3.72-2.95*(C/1.724)*1.0-0.7(1- Sd/100)/1- Sd /100+exp-5.51+22.9(1- Sd/100)*0.1317式中：Sd為砂粒含量，Si 為粉粒含量，Cl為粘粒含量， C為有機質含量。將4個Feature to Raster工具拖入編輯窗口，并拖入一個Single Output

18、 Map Algebra工具，然后將他們按照如下方式連接和命名雙擊Input features導入數據雙擊Feature to Raster分別進行如下設置：在Single Output Map Algebra工具中輸入如下代碼：0.1317 * (0.2 + 0.3 * exp(- 0.0256 * SAND * (1 - SILT / 100) * pow(SILT / (CLAY + SILT),0.3) * (1 - 0.25 * (OM / 1.724) / (OM / 1.724) + exp(3.72 - 2.95 * (OM / 1.724) * (1.0 - 0.7 * (1

19、 - SAND / 100) / (1 - SAND / 100 + exp(- 5.51 + 22.9 * (1 - SAND / 100)模型如下運行這個模型，并將結果添加到視圖中到此，K因子完成。（5）植被覆蓋度因子（C因子）和水保措施因子（P因子）植被覆蓋度因子, 又稱作物經營管理因子。經驗指出, 植被覆蓋度與土壤侵蝕量關系極大。在其它地理環(huán)境因子值相同的情況下, 植被覆蓋度越大, 土壤流失量越小；反之, 則越大。流域的C 因子值賦值如表。表1不同土地利用C因子值土地利用類型旱地水田交通用地和水體草地居民地林地C 因子0. 310.1800.060.20.006水土保持措施因子是采取水

20、保措施后, 土壤流失量與順坡種植時的土壤流失量的比值。通常，包含于這一因子中的控制措施有：等高耕作、等高帶狀種植和修梯田等。將土地利用圖與P值屬性庫文件記錄建立鏈接，再分別將P值賦給土地利用圖，得到P值因子圖。以自然植被P因子為1，坡耕地為0.35，水稻是梯田修筑最好的一種土地利用，P值為0.01。編號土地利用類型11水田12旱地2023林地3133草地43水體52居民地新建一個模型，命名為“CP因子”分別將Feature to Raster工具和Single Output Map Algebra工具拖入編輯窗口，并按如下方式連接它們和命名雙擊Input features進行數據的輸入雙擊Fe

21、ature to Raster工具進行如下設置在Single Output Map Algebra中輸入如下代碼：Con(Landuse = 11,0.18,Con(Landuse = 12,0.31,Con(Landuse = 20 & Landuse = 31 & Landuse = 21 & Landuse = 31 & Landuse = 33,1,0)最終模型為運行這個模型，并將結果添加到視圖中到此，C和P因子結束。（6）基于The Revised Universal Soil Loss Equation (RUSLE)流域土壤侵蝕量預測采用的土壤侵蝕模型修訂的通用土壤流失方程式(R

22、USLE)，其方程式形式簡單，參數容易獲取，且各因子具有明顯的物理解釋意義, 是目前預測土壤侵蝕量較為廣泛使用的方法之一。其公式為: 式中，A 為單位面積上的年均土壤流失量 (thm-2a-1)； R為降雨/徑流侵蝕指數(MJmmhm-2h-1a-1)，用多年平均年降雨侵蝕力指數表示；為土壤可蝕性因子(thm2hMJ-1hm-2mm-1)；為坡長因子；S 為坡度因子；為作物栽培管理因子； 為水土保持工程措施因子。注意: 該方程對于坡長格網中大于800 ft的(244m, 1ft = 0.3048)坡長不適合。分別將“CP因子”“R因子”“K因子”中的模型復制粘貼到“土壤侵蝕模型”中利用將他們排

23、列好然后添加一個Single Output Map Algebra工具，將它與L，S，P，C，K，R連接起來并在Single Output Map Algebra中輸入如下代碼K * R * P * C * L * S * 100最終的結果為到此，A計算完畢。五、實驗結果與分析（1）流域土壤侵蝕等級圖的生成運用GIS的數據庫管理功能和柵格空間分析功能，生成模型的R、L、S、K、C和P共6個因子專題圖，將各因子連乘后，得各像元土壤侵蝕量 (thm2a-1)圖，再乘以系數100，即轉換為tkm2a-1（這個過程直接在模型中的公式中實現了）。對侵蝕柵格圖進行分類，獲得土壤侵蝕等級圖。根據水利部頒布的土壤侵蝕分級標準，將其分為微度、輕度、中度、強度、極強度和劇烈侵蝕6類 (表1)。（在地圖整飾的步驟中實現。）表1流域土壤侵蝕強度分級 侵蝕分級侵蝕模數 (tkm2a-1) 面積 (km2) 百分比 (%) 微度侵蝕 15 000 105.267621.40在上面的模型A之后，再拖入一個Reclassify工具，將其與A相連雙擊Reclassify工具，點擊，在彈出的Classification窗口中將Method設置為Equal Interval，