基于GIS的乌梁素海东岸上游地区水土流失动态变化研究

谭 枭， 王 希， 王秀茹， 刘兰妹， 王红雷

(1.北京林业大学 水土保持学院 教育部水土保持与荒漠化防治重点实验室，北京 100083； 2.北京大学 环境科学与工程学院， 北京 100871)

水土流失是指在水力、风力、重力及冻融等自然营力和人类活动作用下，水土资源和土地生产能力的破坏和损失，包括土地表层侵蚀及水的损失[1]。在中国干旱山区及平原蓄水保土区，由于自然因素及人为不合理开发利用，水土流失日趋严重。乌梁素海是黄河中上游重要的保水、蓄水场地，是世界半荒漠地区极少见的高生态效益的多功能湖泊[2]。目前乌梁素海东岸每年流失表土近0.5 cm，草场沙化严重，沟头不断延伸，严重的水土流失导致湖体萎缩[3]，形成恶性循环。近年来，随着3S技术在水土保持中的广泛应用，利用遥感影像和数字高程模型(digital elevation model，DEM)作为基本信息源[4-5]，利用通用土壤流失方程(USLE)开展土壤侵蚀定量评价和研究[6-8]。同时，基于RS和GIS对水土流失的动态变化展开研究和预测[9-11]，均取得了一定成果，这也是未来区域水土流失评价和监测的主要方法。

本研究以1985，2000和2011年遥感影像和DEM等数据，结合调查资料，从水土流失强度面积的时空变化和各类型水土流失强度的转换两方面研究了乌梁素海东岸上游地区水土流失的动态变化，对水土流失严重区域提出了水土保持措施，并运用马尔可夫模型预测水土流失动态变化。研究有利于掌握区域水土流失变化规律，为挽救乌梁素海采取生态环境保护措施提供科学依据。

1 研究区概况

研究区位于内蒙古自治区巴彦淖尔市乌拉特前旗的乌梁素海东岸上游地区，地处108°41′—109°40′E，40°38′—41°18′N。主要涉及乌拉特前旗的大佘太镇、额尔登布拉格苏木及明安镇3个乡镇，总面积2 451.95 km2。该区位于阴山山地西段，地貌类型多样，其间分布有佘太平原、明安洼地及半固定沙丘等，总地势东北高、西南低。属中温带干旱大陆性气候，多年平均降雨量200～250 mm，时空分布不均，集中在夏季，蒸发量2 167～2 500 mm，远远大于降雨量。土壤类型以灰褐土和栗钙土为主，草本植物主要为针茅(StipacapillataL.)、狗尾草〔Setariaviridis(L.)Beauv〕和骆驼蓬(PeganumharmalaL.)。研究区大范围属阴山山地丘陵蓄水保土区，以沟蚀、面蚀等水力侵蚀为主，中部冲积平原为风、水复合侵蚀，生态环境脆弱。

2 数据与方法

2.1 数据来源与处理

研究采用1985，2000年和2011年Landsat TM影像和1∶1万地形图，结合国家自然科学基金委“中国西部环境与生态科学数据中心”(http://westdc.westgis.ac.cn)1985和2000年的部分数据，校正遥感影像，构建数字高程模型(DEM)。在ERDAS Imagine 9.1软件中通过人机交互解译得到1985，2000和2011年的矢量化土地利用数据，然后在ArcGIS 9.3软件中将其转化为30 m×30 m的栅格数据，再由DEM数据提取坡度、坡长等水土流失因子值，同时调查收集了研究区植被、土壤和水土流失现状数据，依据国家土壤侵蚀分类分级标准(SL190—2007)将研究区水力侵蚀按平均侵蚀模数〔t/(km2·a)〕分为6个强度：<1 000微度(为了表现微度侵蚀的层次，在水土流失分布图上又将其分为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软件中实现数据的采集和处理，通过对上述各因子进行叠加运算，获取各时段水土流失强度分布结果。

(1)降雨侵蚀力因子R。采用Wischmeier提出的基于多年月平均降雨量和多年平均降雨量的经验公式[13]计算R因子：

(2)

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

(2)土壤可蚀性因子K。土壤可蚀性因子K值，通过实地采样分析土壤机械组成及有机质含量，结合Wischmeier等[13]在侵蚀模型中的算法进行估算，推荐使用的简易方程表达式为：

K= 7.504{0.0034+0.0405exp

(3)

Dg=exp(0.01∑filnmi)

(4)

式中：K——土壤可蚀性因子；Dg——土壤颗粒平均粒径(mm)；mi——第i级粒级下组分限值的平均值(mm)；fi——第i级粒级组分的质量百分比(%)。

(3)坡长因子L。利用研究区数字高程模型(DEM)数据，通过ArcGIS 9.3软件中的水文分析模型，提取了研究区坡长因子，公式如下：

L=(λ/22.13)m

(5)

式中：L——坡长因子；λ——像元坡长(m)；m——坡长指数。

像元坡长的计算式如下：

(6)

式中：Di——沿径流方向每像元坡长的水平投影距(m)，在栅格图像中为两相邻像元中心距，随方向而异；θi——每个像元的坡度(°)；i——自山脊像元至待求像元个数。

(4)坡度因子S。

(7)

式中：S——坡度因子，θ——坡度(°)； 运用ArcGIS软件对研究区DEM数据处理，获取坡度图，通过编程语言实现公式(7)程序化，利用栅格计算器计算得出坡度因子S。

(5)水土保持生物措施因子B[14]。利用如下公式计算生物措施因子B与植被覆盖度C之间的关系：

B=0.650 8-0.343 6×lgC

(8)

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

(6)水土保持工程措施因子E和水土保持耕作措施因子T。通过研究区1985，2000和2011年统计的水平沟和鱼鳞坑整地面积，淤地坝数量和谷坊数量等确定工程措施因子E的取值[15]。由于研究区耕地面积仅占约15%，T因子对水土流失量计算结果影响较小。结合不同土地利用类型和耕地坡度可确定T因子的取值。坡度0°～5.5°，T值为0.15；坡度5.5°～13.5°，T值为0.35；坡度大于13.5°，T值取1。

2.2.2 马尔可夫模型 马尔可夫模型是对一种特殊的随机运动过程预测的模型，这一过程每次状态的转移都只与前一时刻的状态有关，而与过去的状态无关[16]。运用马尔可夫模型可以预测乌梁素海东岸上游地区的水土流失动态变化。其数学表达式如下：

V(t2)=V(t1)×Pk

(9)

式中：k——始末年份t1到t2的整步长，预测值的年份为t2+k；V(t1)——系统的初始状态；V(t2)——系统的终止状态；P——转移概率矩阵； 转移概率矩阵P的数学表达式如下：

(10)

3 结果与分析

3.1 水土流失强度分析

基于ArcGIS 9.3软件平台，将中国土壤流失预报方程的7项因子进行栅格计算，得到1985，2000，2011年3期乌梁素海东岸上游地区水土流失强度空间分布图(图1)。水土流失强度颜色加深的区域面积增加十分明显。1985—2011年水土流失发生变化的区域均主要集中于乌拉山和色尔腾山的山间洼地、山前台地等，中部平原风沙区以及乌梁素海东部沿岸农牧区。其中，乌拉山和色尔腾山区域的水土流失强度由1985和2000年以200～1 000 t/(km2·a)的微度侵蚀为主，兼有轻度和中度侵蚀，变为2011年以<200 t/(km2·a)的微度侵蚀为主，几乎没有轻度和中度侵蚀的面积。这主要得益于近年来实施的流域治理措施，如自然封育、植草造林等，提高了山区植被覆盖度，有效地减少了乌梁素海东岸上游山区的水土流失。平原农牧交错带和风沙区的水土流失则逐年加剧，该区域由1985年以轻度侵蚀为主演变为2000年以强烈和轻度、中度侵蚀为主，到2011年强烈侵蚀的面积进一步扩大并出现了极强烈的侵蚀区域。近年来，随着社会经济的快速发展，人口增长，畜牧业的快速发展以及不合理的农业灌溉带来的地下水位快速下降，对乌梁素海东岸上游平原区生态系统的破坏日趋严重，土壤稳定性下降，植被覆盖度降低，是造成该区域水土流失逐年加剧的主要原因。水土流失逐年加剧的区域还有乌梁素海东岸的滨海农牧区，由1985年以微度水土流失为主演变为2011年以强烈水土流失为主。

对3期水土流失强度空间分布图的属性表进行统计分析可知(表1)，1985年微度水土流失面积占总面积的89.6%，到了2000年为81.0%，2011年为60.4%，呈现明显的下降趋势，1985—2000年年均减少0.64%，2000—2011年年均减少1.70%，下降速率逐渐增大。1985，2000和2011年轻度水土流失面积占总面积比例分别为0%，8.2%和21.3%，虽然中度流失面积在研究时段内面积逐年减少，但强烈流失面积却由1985年的0%增长为2011年的15.9%，净增389.23 km2，年均增长14.97 km2。2011年极强烈流失面积也增加为11.52 km2。究其原因，主要是乌梁素海东岸上游山区植被覆盖度逐年降低，土壤松散，降雨径流携带大量泥沙冲刷、淤积沟道，造成了冲积平原等地区严重的水土流失。同时在地形平缓地区，人为频繁活动(开垦、伐林、放牧和基础建设等)改变了原有土地利用类型从而导致土壤有机质流失，理化性质发生变化，加剧了水土流失。

图1 乌梁素海东岸上游地区1985，2000，2011年各水土流失强度空间分布

表1乌梁素海东岸上游地区1985，2000和2011年各水土流失强度面积统计

水土流失强 度 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 水土流失强度类型动态变化分析

表2表示1985—2000年研究区各水土流失强度的面积转换和比例。轻度水土流失面积扩大主要来自微度和中度水土流失面积的转化，分别为167.48和33.70 km2。中度水土流失面积有170.03 km2转化为强烈水土流失的面积。而微度水土流失主要向轻度和中度转化，共转出面积为219.78 km2，占总面积的10.0%。中度水土流失面积转出比例达到了83.4%，为212.02 km2，主要转入强烈水土流失，其次为中度和轻度。由此可见，各种水土流失类型的转换，都主要来自相对更弱的一种或多种水土流失强度类型，这也说明了1985—2000年水土流失加剧。

表2 乌梁素海东岸上游地区1985-2000年各水土流失强度转移面积和比例

表3表示从2000—2011年各水土流失强度类型的转移面积和比例。2011年微度水土流失面积相对于2000年减少了505.94 km2，转换为其他流失强度的总面积为523.07 km2，主要为轻度和中度水土流失，分别为479.49 km2，24.1%和37.26 km2，1.9%。微度水土流失面积逐年减少的原因主要是人类在山区的过度放牧和平原区的过度开垦，破坏了原有地表植被，暴雨冲刷裸露的表土，造成更大的水土流失。

而2000—2011年轻度和中度水土流失面积转化较复杂，其中轻度水土流失主要转换为强烈水土流失，面积为141.69 km2，占其原面积的70.4%，其次还转化为微度和中度水土流失；中度水土流失也主要转换为强烈水土流失，面积为79.89 km2占原面积86.5%。

该时段强烈水土流失面积除了未转换的面积外，主要来源于轻度和中度水土流失的转化，而极强烈水土流失面积主要来源于强烈水土流失的转化。由此可见，各种水土流失强度类型的转换，主要来自相对更轻的另一种或多种水土流失强度类型，特别是强烈水土流失面积主要来自低两个级别的轻度水土流失，这也说明了研究区1985—2000年水土流失状况持续加重。

表3 乌梁素海东岸上游地区2000—2011年各水土流失强度转移面积和比例

3.3 基于马尔可夫模型的水土流失预测

水土流失的动态变化具有马尔可夫过程的性质。马尔可夫模型是根据变量目前的状态去预测未来可能的变化，只需要变量最近的动态资料，而不需要连续的历史资料。因此，水土流失的动态变化满足马尔可夫模型的特点，此模型可以用于水土流失动态变化预测。运用马尔可夫模型预测水土流失，需要确定初始状态矩阵。因研究区1985年水土流失强度类型缺失了轻度、强烈和极强烈3种，故选择2000年不同水土流失强度的面积占全部流失面积的百分比作为各强度类型的初始概率，组成初始状态矩阵。

而对水土流失的动态变化预测关键在于转换概率矩阵P的确定。利用表3中2000—2011年水土流失面积的变化计算出各流失强度的平均转化率，从而得到转换概率矩阵(表4)。

表4 乌梁素海东岸上游地区2000-2011年水土流失强度类型转换概率矩阵

以2011年各水土流失强度面积所占比例为初始向量，2000—2011年水土流失强度类型转换概率矩阵为预测时段的转换概率矩阵，12 a为一步长。预测结果，2022年微度、轻度、中度、强烈和极强烈水土流失的面积分别为1 130.15，453.16，60.59，779.38和28.67 km2(表5)。2011—2022年微度水土流失面积依然减少，为350.15 km2，相比于2000—2011年的衰减程度变缓。轻度水土流失面积相比过去的2000—2011年呈现了一定的起伏，减少68.37 km2，但依然可以预测乌梁素海东岸上游地区水土流失微度和轻度的区域面积仍会进一步减少，水土流失的治理形势还很严峻。2011—2022年强烈水土流失面积大幅增加390.15 km2，表明研究区水土流失严重的区域亟需治理。同时，到2022年中度和极强烈水土流失面积分别增加11.22和17.15 km2。

表5 乌梁素海东岸上游地区2022年水土流失变化预测

4 结 论

(1)通过对乌梁素海东岸上游地区1985—2000年，2000—2011年不同强度水土流失面积的空间变化以及各强度类型转换情况的分析发现，研究区水土流失呈现逐年加剧的趋势。微度水土流失面积逐年减少，由1985年的2 197.73 km2减少至2011年的1 480.30 km2，主要转换为轻度和中度水土流失，而强烈水土流失面积大幅增加达389.23 km2。各类型水土流失强度面积的转换中，更剧烈的一种类型均是主要来自低一级甚至二级的转化，这也说明了研究区水土流失加重的严峻形势。乌梁素海东岸上游地区水土流失严重的区域主要集中于中部农牧交错区、风沙区以及沿岸农牧区，而乌拉山和色尔腾山的山区水土流失已经逐步减弱，应尽快对以上水土流失持续加剧地区采取综合治理措施，保护生态环境。

(2)基于马尔可夫模型对水土流失强度类型的动态变化预测分析，2022年乌梁素海东岸上游地区的水土流失将出现两极分化的形势，微度水土流失面积相比2011年减少350.15 km2，强烈水土流失面积相比2011年增加390.15 km2，微度和轻度水土流失面积进一步减少，强烈和极强烈水土流失面积仍在急剧扩大。应以乌梁素海东岸人口密集的农牧区为治理重点，发展节水灌溉及牧草改良，营建水土保持林网，配合山区封育、植树种草和修建谷坊、淤地坝等流域治理措施，生物措施与工程措施结合，开展综合治理，方能有效地遏制水土流失。

[参考文献]

[1]余新晓,毕华兴.水土保持学[M].3版.北京:中国林业出版社,2013.

[2]李兴.内蒙古乌梁素海水质动态数值模拟研究[D].呼和浩特：内蒙古农业大学,2009.

[3]莫日根,童伟,段瑞琴,等.乌梁素海生态环境存在的问题和治理措施[J].北方环境,2012,24(4):18-22.

[4]张明阳,王克林,陈洪松.基于RS和GIS的喀斯特区域水土流失动态监测与分析:以广西环江县为例[J].资源科学,2007,29(3):124-131.

[5]魏兴萍.基于RS 和GIS 的重庆南川区水土流失变化研究[J].水土保持研究,2009,16(5):60-65.

[6]张建香,张勃,张华,等.黄土高原的景观格局变化与水土流失研究:以黄土高原马莲河流域为例[J].自然资源学报,2011,26(9):1513-1525.

[7]王文娟,张树文,李颖,等.基于GIS和USLE的三江平原土壤侵蚀定量评价[J].干旱区资源与环境,2008,22(9):112-117.

[8]邸利,张仁陟,张富,等.基于RS与GIS的定西市安定区土地利用变化与土壤侵蚀研究[J].干旱区资源与环境,2011,25(2):40-45.

[9]冯晓刚,李锐.西安市水土流失动态变化及对策研究[J].水土保持通报,2010,30(6):107-111.

[10]石香琼,査轩,陈世发,等.基于马尔柯夫模型的红壤退化地水土流失动态变化预测研究[J].水土保持研究,2009,16(4):19-23.

[11]马义娟,苏志珠.晋西沿黄地区水土流失动态变化及成因分析[J].干旱区资源与环境,2004,18(1):122-128.

[12]刘宝元,谢云,张科利.土壤侵蚀预报模型[M].北京:中国科学技术出版社,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模型与地理信息系统IDRISI预测小流域土壤侵蚀量的研究[J].水土保持学报,2000,24(2):20-25.

[15]杨洁,汪邦稳.赣南地区水土流失时空变化和评价研究[J].中国水土保持,2011(12):10-12.

[16]胡希军.城市化主导的景观结构演变机制研究：以义乌市为例[D].长沙：中南林业科技大学,2006.