基于GIS的烏梁素海東岸上游地區(qū)水土流失動態(tài)變化研究

譚 梟， 王 希， 王秀茹， 劉蘭妹， 王紅雷

(1.北京林業(yè)大學 水土保持學院 教育部水土保持與荒漠化防治重點實驗室，北京 100083； 2.北京大學 環(huán)境科學與工程學院， 北京 100871)

水土流失是指在水力、風力、重力及凍融等自然營力和人類活動作用下，水土資源和土地生產(chǎn)能力的破壞和損失，包括土地表層侵蝕及水的損失[1]。在中國干旱山區(qū)及平原蓄水保土區(qū)，由于自然因素及人為不合理開發(fā)利用，水土流失日趨嚴重。烏梁素海是黃河中上游重要的保水、蓄水場地，是世界半荒漠地區(qū)極少見的高生態(tài)效益的多功能湖泊[2]。目前烏梁素海東岸每年流失表土近0.5 cm，草場沙化嚴重，溝頭不斷延伸，嚴重的水土流失導致湖體萎縮[3]，形成惡性循環(huán)。近年來，隨著3S技術(shù)在水土保持中的廣泛應用，利用遙感影像和數(shù)字高程模型(digital elevation model，DEM)作為基本信息源[4-5]，利用通用土壤流失方程(USLE)開展土壤侵蝕定量評價和研究[6-8]。同時，基于RS和GIS對水土流失的動態(tài)變化展開研究和預測[9-11]，均取得了一定成果，這也是未來區(qū)域水土流失評價和監(jiān)測的主要方法。

本研究以1985，2000和2011年遙感影像和DEM等數(shù)據(jù)，結(jié)合調(diào)查資料，從水土流失強度面積的時空變化和各類型水土流失強度的轉(zhuǎn)換兩方面研究了烏梁素海東岸上游地區(qū)水土流失的動態(tài)變化，對水土流失嚴重區(qū)域提出了水土保持措施，并運用馬爾可夫模型預測水土流失動態(tài)變化。研究有利于掌握區(qū)域水土流失變化規(guī)律，為挽救烏梁素海采取生態(tài)環(huán)境保護措施提供科學依據(jù)。

1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古自治區(qū)巴彥淖爾市烏拉特前旗的烏梁素海東岸上游地區(qū)，地處108°41′—109°40′E，40°38′—41°18′N。主要涉及烏拉特前旗的大佘太鎮(zhèn)、額爾登布拉格蘇木及明安鎮(zhèn)3個鄉(xiāng)鎮(zhèn)，總面積2 451.95 km2。該區(qū)位于陰山山地西段，地貌類型多樣，其間分布有佘太平原、明安洼地及半固定沙丘等，總地勢東北高、西南低。屬中溫帶干旱大陸性氣候，多年平均降雨量200～250 mm，時空分布不均，集中在夏季，蒸發(fā)量2 167～2 500 mm，遠遠大于降雨量。土壤類型以灰褐土和栗鈣土為主，草本植物主要為針茅(StipacapillataL.)、狗尾草〔Setariaviridis(L.)Beauv〕和駱駝蓬(PeganumharmalaL.)。研究區(qū)大范圍屬陰山山地丘陵蓄水保土區(qū)，以溝蝕、面蝕等水力侵蝕為主，中部沖積平原為風、水復合侵蝕，生態(tài)環(huán)境脆弱。

2 數(shù)據(jù)與方法

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

研究采用1985，2000年和2011年Landsat TM影像和1∶1萬地形圖，結(jié)合國家自然科學基金委“中國西部環(huán)境與生態(tài)科學數(shù)據(jù)中心”(http://westdc.westgis.ac.cn)1985和2000年的部分數(shù)據(jù)，校正遙感影像，構(gòu)建數(shù)字高程模型(DEM)。在ERDAS Imagine 9.1軟件中通過人機交互解譯得到1985，2000和2011年的矢量化土地利用數(shù)據(jù)，然后在ArcGIS 9.3軟件中將其轉(zhuǎn)化為30 m×30 m的柵格數(shù)據(jù)，再由DEM數(shù)據(jù)提取坡度、坡長等水土流失因子值，同時調(diào)查收集了研究區(qū)植被、土壤和水土流失現(xiàn)狀數(shù)據(jù)，依據(jù)國家土壤侵蝕分類分級標準(SL190—2007)將研究區(qū)水力侵蝕按平均侵蝕模數(shù)〔t/(km2·a)〕分為6個強度：<1 000微度(為了表現(xiàn)微度侵蝕的層次，在水土流失分布圖上又將其分為0～200和200～1 000兩個層級)；1 000～2 500輕度；2 500～5 000中度；5 000～8 000強烈；8 000～15 000極強烈。

2.2 研究方法

2.2.1 水土流失定量研究 水土流失定量研究采用中國土壤流失預報方程[12](China soil loss erosion， CSLE)，表達式如下：

A=R×K×L×S×B×E×T

(1)

式中：A——年平均土壤水力侵蝕流失量〔t/(km2·a)〕；R——降雨侵蝕力因子〔t/(km2·a)〕；K——土壤可蝕性因子；L——坡長因子；S——坡度因子；B——水土保持生物措施因子；E——水土保持工程措施因子；T——水土保持耕作措施因子。在ArcGIS 9.3軟件中實現(xiàn)數(shù)據(jù)的采集和處理，通過對上述各因子進行疊加運算，獲取各時段水土流失強度分布結(jié)果。

(1)降雨侵蝕力因子R。采用Wischmeier提出的基于多年月平均降雨量和多年平均降雨量的經(jīng)驗公式[13]計算R因子：

(2)

式中：R——降雨侵蝕力因子；i,j——月份；Pmj——多年月平均降水量(mm)；Pm——多年平均降水量(mm)。

(2)土壤可蝕性因子K。土壤可蝕性因子K值，通過實地采樣分析土壤機械組成及有機質(zhì)含量，結(jié)合Wischmeier等[13]在侵蝕模型中的算法進行估算，推薦使用的簡易方程表達式為：

K= 7.504{0.0034+0.0405exp

(3)

Dg=exp(0.01∑filnmi)

(4)

式中：K——土壤可蝕性因子；Dg——土壤顆粒平均粒徑(mm)；mi——第i級粒級下組分限值的平均值(mm)；fi——第i級粒級組分的質(zhì)量百分比(%)。

(3)坡長因子L。利用研究區(qū)數(shù)字高程模型(DEM)數(shù)據(jù)，通過ArcGIS 9.3軟件中的水文分析模型，提取了研究區(qū)坡長因子，公式如下：

L=(λ/22.13)m

(5)

式中：L——坡長因子；λ——像元坡長(m)；m——坡長指數(shù)。

像元坡長的計算式如下：

(6)

式中：Di——沿徑流方向每像元坡長的水平投影距(m)，在柵格圖像中為兩相鄰像元中心距，隨方向而異；θi——每個像元的坡度(°)；i——自山脊像元至待求像元個數(shù)。

(4)坡度因子S。

(7)

式中：S——坡度因子，θ——坡度(°)； 運用ArcGIS軟件對研究區(qū)DEM數(shù)據(jù)處理，獲取坡度圖，通過編程語言實現(xiàn)公式(7)程序化，利用柵格計算器計算得出坡度因子S。

(5)水土保持生物措施因子B[14]。利用如下公式計算生物措施因子B與植被覆蓋度C之間的關(guān)系：

B=0.650 8-0.343 6×lgC

(8)

式中：B——水土保持生物措施因子；C——植被覆蓋度(利用NDVI計算)。

(6)水土保持工程措施因子E和水土保持耕作措施因子T。通過研究區(qū)1985，2000和2011年統(tǒng)計的水平溝和魚鱗坑整地面積，淤地壩數(shù)量和谷坊數(shù)量等確定工程措施因子E的取值[15]。由于研究區(qū)耕地面積僅占約15%，T因子對水土流失量計算結(jié)果影響較小。結(jié)合不同土地利用類型和耕地坡度可確定T因子的取值。坡度0°～5.5°，T值為0.15；坡度5.5°～13.5°，T值為0.35；坡度大于13.5°，T值取1。

2.2.2 馬爾可夫模型 馬爾可夫模型是對一種特殊的隨機運動過程預測的模型，這一過程每次狀態(tài)的轉(zhuǎn)移都只與前一時刻的狀態(tài)有關(guān)，而與過去的狀態(tài)無關(guān)[16]。運用馬爾可夫模型可以預測烏梁素海東岸上游地區(qū)的水土流失動態(tài)變化。其數(shù)學表達式如下：

V(t2)=V(t1)×Pk

(9)

式中：k——始末年份t1到t2的整步長，預測值的年份為t2+k；V(t1)——系統(tǒng)的初始狀態(tài)；V(t2)——系統(tǒng)的終止狀態(tài)；P——轉(zhuǎn)移概率矩陣； 轉(zhuǎn)移概率矩陣P的數(shù)學表達式如下：

(10)

3 結(jié)果與分析

3.1 水土流失強度分析

基于ArcGIS 9.3軟件平臺，將中國土壤流失預報方程的7項因子進行柵格計算，得到1985，2000，2011年3期烏梁素海東岸上游地區(qū)水土流失強度空間分布圖(圖1)。水土流失強度顏色加深的區(qū)域面積增加十分明顯。1985—2011年水土流失發(fā)生變化的區(qū)域均主要集中于烏拉山和色爾騰山的山間洼地、山前臺地等，中部平原風沙區(qū)以及烏梁素海東部沿岸農(nóng)牧區(qū)。其中，烏拉山和色爾騰山區(qū)域的水土流失強度由1985和2000年以200～1 000 t/(km2·a)的微度侵蝕為主，兼有輕度和中度侵蝕，變?yōu)?011年以<200 t/(km2·a)的微度侵蝕為主，幾乎沒有輕度和中度侵蝕的面積。這主要得益于近年來實施的流域治理措施，如自然封育、植草造林等，提高了山區(qū)植被覆蓋度，有效地減少了烏梁素海東岸上游山區(qū)的水土流失。平原農(nóng)牧交錯帶和風沙區(qū)的水土流失則逐年加劇，該區(qū)域由1985年以輕度侵蝕為主演變?yōu)?000年以強烈和輕度、中度侵蝕為主，到2011年強烈侵蝕的面積進一步擴大并出現(xiàn)了極強烈的侵蝕區(qū)域。近年來，隨著社會經(jīng)濟的快速發(fā)展，人口增長，畜牧業(yè)的快速發(fā)展以及不合理的農(nóng)業(yè)灌溉帶來的地下水位快速下降，對烏梁素海東岸上游平原區(qū)生態(tài)系統(tǒng)的破壞日趨嚴重，土壤穩(wěn)定性下降，植被覆蓋度降低，是造成該區(qū)域水土流失逐年加劇的主要原因。水土流失逐年加劇的區(qū)域還有烏梁素海東岸的濱海農(nóng)牧區(qū)，由1985年以微度水土流失為主演變?yōu)?011年以強烈水土流失為主。

對3期水土流失強度空間分布圖的屬性表進行統(tǒng)計分析可知(表1)，1985年微度水土流失面積占總面積的89.6%，到了2000年為81.0%，2011年為60.4%，呈現(xiàn)明顯的下降趨勢，1985—2000年年均減少0.64%，2000—2011年年均減少1.70%，下降速率逐漸增大。1985，2000和2011年輕度水土流失面積占總面積比例分別為0%，8.2%和21.3%，雖然中度流失面積在研究時段內(nèi)面積逐年減少，但強烈流失面積卻由1985年的0%增長為2011年的15.9%，凈增389.23 km2，年均增長14.97 km2。2011年極強烈流失面積也增加為11.52 km2。究其原因，主要是烏梁素海東岸上游山區(qū)植被覆蓋度逐年降低，土壤松散，降雨徑流攜帶大量泥沙沖刷、淤積溝道，造成了沖積平原等地區(qū)嚴重的水土流失。同時在地形平緩地區(qū)，人為頻繁活動(開墾、伐林、放牧和基礎(chǔ)建設(shè)等)改變了原有土地利用類型從而導致土壤有機質(zhì)流失，理化性質(zhì)發(fā)生變化，加劇了水土流失。

圖1 烏梁素海東岸上游地區(qū)1985，2000，2011年各水土流失強度空間分布

表1烏梁素海東岸上游地區(qū)1985，2000和2011年各水土流失強度面積統(tǒng)計

水土流失強 度 1985年面積/km2比例/%2000年面積/km2比例/%2011年面積/km2比例/%1985—2000年變化面積/km22000—2011年變化面積/km2微 度2 197.73 89.601 986.24 81.001 480.3060.40-211.49-505.94輕 度0.00 —201.18 8.20521.5321.30201.18320.35中 度254.22 10.4092.32 3.8049.372.00-161.90-42.95強 烈0.00 —172.21 7.00389.2315.90172.21217.02極強烈0.00 —0.00 —11.520.500.0011.52合 計2 451.95 100.002 451.95 100.00 2 451.95100.00——

3.2 水土流失強度類型動態(tài)變化分析

表2表示1985—2000年研究區(qū)各水土流失強度的面積轉(zhuǎn)換和比例。輕度水土流失面積擴大主要來自微度和中度水土流失面積的轉(zhuǎn)化，分別為167.48和33.70 km2。中度水土流失面積有170.03 km2轉(zhuǎn)化為強烈水土流失的面積。而微度水土流失主要向輕度和中度轉(zhuǎn)化，共轉(zhuǎn)出面積為219.78 km2，占總面積的10.0%。中度水土流失面積轉(zhuǎn)出比例達到了83.4%，為212.02 km2，主要轉(zhuǎn)入強烈水土流失，其次為中度和輕度。由此可見，各種水土流失類型的轉(zhuǎn)換，都主要來自相對更弱的一種或多種水土流失強度類型，這也說明了1985—2000年水土流失加劇。

表2 烏梁素海東岸上游地區(qū)1985-2000年各水土流失強度轉(zhuǎn)移面積和比例

表3表示從2000—2011年各水土流失強度類型的轉(zhuǎn)移面積和比例。2011年微度水土流失面積相對于2000年減少了505.94 km2，轉(zhuǎn)換為其他流失強度的總面積為523.07 km2，主要為輕度和中度水土流失，分別為479.49 km2，24.1%和37.26 km2，1.9%。微度水土流失面積逐年減少的原因主要是人類在山區(qū)的過度放牧和平原區(qū)的過度開墾，破壞了原有地表植被，暴雨沖刷裸露的表土，造成更大的水土流失。

而2000—2011年輕度和中度水土流失面積轉(zhuǎn)化較復雜，其中輕度水土流失主要轉(zhuǎn)換為強烈水土流失，面積為141.69 km2，占其原面積的70.4%，其次還轉(zhuǎn)化為微度和中度水土流失；中度水土流失也主要轉(zhuǎn)換為強烈水土流失，面積為79.89 km2占原面積86.5%。

該時段強烈水土流失面積除了未轉(zhuǎn)換的面積外，主要來源于輕度和中度水土流失的轉(zhuǎn)化，而極強烈水土流失面積主要來源于強烈水土流失的轉(zhuǎn)化。由此可見，各種水土流失強度類型的轉(zhuǎn)換，主要來自相對更輕的另一種或多種水土流失強度類型，特別是強烈水土流失面積主要來自低兩個級別的輕度水土流失，這也說明了研究區(qū)1985—2000年水土流失狀況持續(xù)加重。

表3 烏梁素海東岸上游地區(qū)2000—2011年各水土流失強度轉(zhuǎn)移面積和比例

3.3 基于馬爾可夫模型的水土流失預測

水土流失的動態(tài)變化具有馬爾可夫過程的性質(zhì)。馬爾可夫模型是根據(jù)變量目前的狀態(tài)去預測未來可能的變化，只需要變量最近的動態(tài)資料，而不需要連續(xù)的歷史資料。因此，水土流失的動態(tài)變化滿足馬爾可夫模型的特點，此模型可以用于水土流失動態(tài)變化預測。運用馬爾可夫模型預測水土流失，需要確定初始狀態(tài)矩陣。因研究區(qū)1985年水土流失強度類型缺失了輕度、強烈和極強烈3種，故選擇2000年不同水土流失強度的面積占全部流失面積的百分比作為各強度類型的初始概率，組成初始狀態(tài)矩陣。

而對水土流失的動態(tài)變化預測關(guān)鍵在于轉(zhuǎn)換概率矩陣P的確定。利用表3中2000—2011年水土流失面積的變化計算出各流失強度的平均轉(zhuǎn)化率，從而得到轉(zhuǎn)換概率矩陣(表4)。

表4 烏梁素海東岸上游地區(qū)2000-2011年水土流失強度類型轉(zhuǎn)換概率矩陣

以2011年各水土流失強度面積所占比例為初始向量，2000—2011年水土流失強度類型轉(zhuǎn)換概率矩陣為預測時段的轉(zhuǎn)換概率矩陣，12 a為一步長。預測結(jié)果，2022年微度、輕度、中度、強烈和極強烈水土流失的面積分別為1 130.15，453.16，60.59，779.38和28.67 km2(表5)。2011—2022年微度水土流失面積依然減少，為350.15 km2，相比于2000—2011年的衰減程度變緩。輕度水土流失面積相比過去的2000—2011年呈現(xiàn)了一定的起伏，減少68.37 km2，但依然可以預測烏梁素海東岸上游地區(qū)水土流失微度和輕度的區(qū)域面積仍會進一步減少，水土流失的治理形勢還很嚴峻。2011—2022年強烈水土流失面積大幅增加390.15 km2，表明研究區(qū)水土流失嚴重的區(qū)域亟需治理。同時，到2022年中度和極強烈水土流失面積分別增加11.22和17.15 km2。

表5 烏梁素海東岸上游地區(qū)2022年水土流失變化預測

4 結(jié) 論

(1)通過對烏梁素海東岸上游地區(qū)1985—2000年，2000—2011年不同強度水土流失面積的空間變化以及各強度類型轉(zhuǎn)換情況的分析發(fā)現(xiàn)，研究區(qū)水土流失呈現(xiàn)逐年加劇的趨勢。微度水土流失面積逐年減少，由1985年的2 197.73 km2減少至2011年的1 480.30 km2，主要轉(zhuǎn)換為輕度和中度水土流失，而強烈水土流失面積大幅增加達389.23 km2。各類型水土流失強度面積的轉(zhuǎn)換中，更劇烈的一種類型均是主要來自低一級甚至二級的轉(zhuǎn)化，這也說明了研究區(qū)水土流失加重的嚴峻形勢。烏梁素海東岸上游地區(qū)水土流失嚴重的區(qū)域主要集中于中部農(nóng)牧交錯區(qū)、風沙區(qū)以及沿岸農(nóng)牧區(qū)，而烏拉山和色爾騰山的山區(qū)水土流失已經(jīng)逐步減弱，應盡快對以上水土流失持續(xù)加劇地區(qū)采取綜合治理措施，保護生態(tài)環(huán)境。

(2)基于馬爾可夫模型對水土流失強度類型的動態(tài)變化預測分析，2022年烏梁素海東岸上游地區(qū)的水土流失將出現(xiàn)兩極分化的形勢，微度水土流失面積相比2011年減少350.15 km2，強烈水土流失面積相比2011年增加390.15 km2，微度和輕度水土流失面積進一步減少，強烈和極強烈水土流失面積仍在急劇擴大。應以烏梁素海東岸人口密集的農(nóng)牧區(qū)為治理重點，發(fā)展節(jié)水灌溉及牧草改良，營建水土保持林網(wǎng)，配合山區(qū)封育、植樹種草和修建谷坊、淤地壩等流域治理措施，生物措施與工程措施結(jié)合，開展綜合治理，方能有效地遏制水土流失。

[參考文獻]

[1]余新曉,畢華興.水土保持學[M].3版.北京:中國林業(yè)出版社,2013.

[2]李興.內(nèi)蒙古烏梁素海水質(zhì)動態(tài)數(shù)值模擬研究[D].呼和浩特：內(nèi)蒙古農(nóng)業(yè)大學,2009.

[3]莫日根,童偉,段瑞琴,等.烏梁素海生態(tài)環(huán)境存在的問題和治理措施[J].北方環(huán)境,2012,24(4):18-22.

[4]張明陽,王克林,陳洪松.基于RS和GIS的喀斯特區(qū)域水土流失動態(tài)監(jiān)測與分析:以廣西環(huán)江縣為例[J].資源科學,2007,29(3):124-131.

[5]魏興萍.基于RS 和GIS 的重慶南川區(qū)水土流失變化研究[J].水土保持研究,2009,16(5):60-65.

[6]張建香,張勃,張華,等.黃土高原的景觀格局變化與水土流失研究:以黃土高原馬蓮河流域為例[J].自然資源學報,2011,26(9):1513-1525.

[7]王文娟,張樹文,李穎,等.基于GIS和USLE的三江平原土壤侵蝕定量評價[J].干旱區(qū)資源與環(huán)境,2008,22(9):112-117.

[8]邸利,張仁陟,張富,等.基于RS與GIS的定西市安定區(qū)土地利用變化與土壤侵蝕研究[J].干旱區(qū)資源與環(huán)境,2011,25(2):40-45.

[9]馮曉剛,李銳.西安市水土流失動態(tài)變化及對策研究[J].水土保持通報,2010,30(6):107-111.

[10]石香瓊,査軒,陳世發(fā),等.基于馬爾柯夫模型的紅壤退化地水土流失動態(tài)變化預測研究[J].水土保持研究,2009,16(4):19-23.

[11]馬義娟,蘇志珠.晉西沿黃地區(qū)水土流失動態(tài)變化及成因分析[J].干旱區(qū)資源與環(huán)境,2004,18(1):122-128.

[12]劉寶元,謝云,張科利.土壤侵蝕預報模型[M].北京:中國科學技術(shù)出版社,2001.

[13]Wischmeier W H, Johnson C B, Cross B V. Soil erodibility nomograph for farm land and construction sites[J]. Journal of Soil and Water Conservation, 1971,26(5):189-193.

[14]蔡崇法,丁樹文,史志華,等.應用USLE模型與地理信息系統(tǒng)IDRISI預測小流域土壤侵蝕量的研究[J].水土保持學報,2000,24(2):20-25.

[15]楊潔,汪邦穩(wěn).贛南地區(qū)水土流失時空變化和評價研究[J].中國水土保持,2011(12):10-12.

[16]胡希軍.城市化主導的景觀結(jié)構(gòu)演變機制研究：以義烏市為例[D].長沙：中南林業(yè)科技大學,2006.