(完整word版)基于RUSLE的土壤侵蚀建模分析.docx

空间信息应用实践（中级）实验指导书

空间建模——基于 RUSLE 的土壤侵蚀建模分析

一．实验背景

Soil erosion and gullying in the upper Panuco basin, Sierra Madre Oriental, eastern Mexico 土壤侵蚀是地球表面物质运动的一种自然现象，全球除永冻地区外，均发生不同程度的

土壤侵蚀。人类社会出现后，土壤侵蚀成为自然和人为活动共同作用下的一种动态过程，构成了特殊的侵蚀环境背景，并伴随着人类对自然改造能力的增强，逐渐成为当今世界资源和

环境可持续发展所面临的重要问题之一。

土壤侵蚀被称为“蠕动的灾难”，每年因土壤侵蚀造成的经济损失较诸如滑坡、泥石流和地震等地质灾害更大, 土壤侵蚀已成为我国乃至全球的重大环境问题之一。

土壤侵蚀及其产生的泥沙使土壤养分流失、土地生产力下降、湖泊淤积、江河堵塞，并造成诸如洪水等自然灾害，泥沙携带的大量营养物和污染物质加剧了水体富营养化，水质恶化，不断严重威胁到人类的生存。

据估计全球每年因土壤侵蚀损失300 万公顷土地的生产力，造成的损失以百亿美元计。

我国人口众多、农耕历史悠久，加之历史上战乱频仍，以黄土高原为代表的华夏文明发源地

是世界上土壤侵蚀最严重的区域之一，1990 年遥感普查结果，全国水土流失面积达367 万km 2，占国土总面积的38.2％，其中 50%为水蚀地区，土壤侵蚀以黄土高原、四川紫色土地

亿 t 。水区和华南红壤地区尤为突出，仅黄土高原地区一处，平均每年流失泥沙就达到16.3 土流

失已成为中国重要的环境问题，土壤侵蚀研究已成为目前环境保护中的一个重要课题。

土壤侵蚀预报是有效监测水土流失和评价水保措施效益的手段，侵蚀模型则是进行土壤流失监测和预报的重要工具。然而传统预测方法需要在量经费、时间和人力的投入，因此，

在一定精度范围内通过有限的数据输入，得到满足要求的土壤侵蚀预测结果成为趋势。80年代以来，随着地理信息系统(Geographical Information System, GIS) 的成熟，它开始与土壤侵蚀模型—通用土壤流失方程(Universal Soil Loss Equation, USLE) 相结合进行流域土壤侵蚀量的预测和估算，业已成为土壤侵蚀动态研究的有力工具。GIS 与 USLE 相结合的分布式方法运用 GIS 的栅格数据分析功能，可预测出每个栅格的土壤侵蚀量，便于管理者识别

关键源区，并通过确定引起水土流失的关键因子，针对性地提出最佳管理措施(Best Management Practices，BMPs) ，为流域内土地资源的质量评价、利用规划和经营管理等提供科学依据与决策手段。

二、实验目的

模型生成器(ModelBuilder)为设计和实现空间处理模型提供了一个图形化的建模环境。

模型是以流程图的形式表示，它通过工具将数据串起来以创建高级的功能和流程。你可以将工具和数据集拖动到一个模型中，然后按照有序的步骤把它们连接起来以实现复杂的GIS 任务。通过对本次练习达到以下目的：

掌握如何在ModelBuilder 环境下通过绘制数据处理流程图的方式实现空间分析过程的

自动化；

掌握土壤侵蚀理论的基本知识；

掌握利用脚本文件实现空间建模，加深对地理建模过程的认识，对各种GIS 分析工具的用途有深入的理解；

在 ModelBuilder 环境下如何计算RUSLE 模型的中各个因子，实现 RUSLE 模型自动化；

三、实验准备

实验环境： ArcGIS Desktop 9.3

实验数据：矢量和栅格数据

矢量数据：研究区界线 (bj.shp)、气象数据（Climate.shp ），土地利用数据 (landuse_Clip.shp,)和土壤数据 (soil_clip) ；

栅格数据：地形数据（DEM ）；

四、实验内容与步骤

(1) 实验准备

ArcGIS进行建模。本次试验需要使用ArcGIS 的建模功能，在实验之前需要掌握如何利用

首先，打开ArcMap ，激活工具箱

”，即可新建一个工

在工具箱中右键单击，选择“ New Toolbox

具箱。可以在此工具箱上右击，通过“ Rename”对工具箱重命名。

在新建的工具箱上右击，按照“ New ”—— >“ Model ”新建一个 Model ，可以按照同样的方法给这个 Model 命名。

然后在此Model 上右击，通过“Edit ”进入模型的编辑模式。

到此，模型准备已经结束，接下来开始逐个建立模型的各个因子。 （ 2）地形因子（ L ， S 因子）

算法：坡长因子采用公式计算 , L

( l )m ，式中： L 为坡长因子， l 为像元坡长， m 为

坡 22.13

i

i 1

长指数，像元坡长的计算式如下：

l i

( D i / cos i )(D i / cos i ) D i / cos i ，

m 取

1

1

值如下式：

0.5 5%

0.4 3% 5% 为像元坡度 (%)

m

1% 式中， 0.3 3%

0.2 1%

式中， l i 为像元坡长， D i 为沿径流方向每像元坡长的水平投影距 (在栅格图像中为两相邻像

元中心距，随方向而异

)， θi 为每个像元的坡度 ( °)， i 为自山脊像元至待求像元个数。

10.8sin

0.03 5

坡度因子 S 分段计算： S 16.8sin

- 0.55 10 21.91sin

0.96 10

L 和 S 因子的模型建立：

首先在工具箱中找到

Resample 工具，可以使用工具箱自带的搜索功能快速定位到。在工具

箱的下方有一行标签 ，选择 Search 标签，在搜索框中输入要查找的工具名，如 Resample ，点击

Search 进行查询，查询结束后选中查询结果，点击

下方的 Locate 可以快速定位需要查找的工具。

可以将这个工具直接拖到Model 的编辑窗口中，如图：

现在需要给这个工具添加一个参数，在编辑窗口的Resample 上右击，通过“ Make Variable ”

—— >“ From Parameter”—— >“ Input Raster”添加。

注意：这里不建议使用右键菜单的“Create Variable”来添加输入输出参数，因为很多工具

拖入到编辑窗口后会自带一个输出参数，而且它们也有自己的默认输入参数。如果另外新建一个参数，可能会因为这个新建参数类型不与工具要求的输入参数类型对应而出现错误。

按照同样的方法拖入Slope 工具， Single Output Map Algebra工具。

通过编辑窗口上的工具将这些工具首尾连接起来。

双击 Input Raster ，输入 dem 数据

输入数据之后，编辑窗口中的工具颜色会相应的变化，说明这些工具已经相互连接起来，还是白色的工具代表它还没有和前面的工具联系起来构成“流水线”，同时，这也是判断Single

Output Map Algebra 工具中的脚本语言是否和前面的输出文件关联起来的依据。在相互连接

的工具中，只要有一个工具是白色的，就说明这条“流水线”不能正常运行。

可以发现 Single Output Map Algebra 工具还是白的，这是因为我们没有添加算法，下面添加

用于计算 S 因子的算法，依据为：

10.8sin0.035

坡度因子 S分段计算：S16.8sin- 0.5510

21.91sin0.9610

双击 Single Output Map Algebra 工具，添加如下代码：

Con([Slope_degree1]< 5 ,10.8 * Sin([Slope_degree1]* 3.14/ 180) + 0.03,Con([Slope_degree1]>=

5 & [Slope_degree1]< 10, 16.8 * Sin([Slope_degree1]* 3.14 /180) - 0.5, 21.91 * Sin([Slope_degree1]

* 3.14 / 180) - 0.96))

这时点击左上方的绿色圆点会发现有错误提示

单击

会看到 Single Output Map Algebra 工具所接收到所有参数

我，其中并没有我算所需要的个 Output raster 参数名行修改，改Slope_degree1 ，所以会提示。就需要Slope_degree1 。修改个参数其就是把

Single Output Map Algebra 工具前的Output raster 出框行改名，如：

接下来我要行一些参数的置。（只需要修改框中的参数，其他的采用默置）

双窗口中的Resample 工具，行如下置，注意，里置格大小90，并且本中个置都一采用90。 .

Slope 置

注意里的是DEGREE （度数）。其坡度有两种表示方式，一种是用我平常用的

度数；有一种是用百分数表示。DEGREE表示是用度数来表示坡度，即的。

（里要解一下代：在里我使用的是条件句CON ，似于 C 言中的if 句，其写法CON（条件 1,如果条件真行，如果条件假行）。如果多个条件行

嵌套，就要写成CON(W1,T1,CON(W2,T2,CON(W3,T3,CON(W4,T4,??))))W代表条件，T

代表条件真行的句）

（关于变量：这里的变量要用[] 括起来，如 [Slope_degree1] > 5 ）

注意：之所以使用 [Slope_degree1] * 3.14 / 180 这是将原来的角度制转化为弧度制，

能识别角度制。

注意：在运算符（如+， -， * ，/）的左右要有空格，如[Slope_degree1] * 3.14 / 180计算机不不要写成

[Slope_degree1]*3.14/180 。前者的运算符左右有空格，后者的运算符左右没有空格。这一点必须严格遵守，否则相同的代码会出现不同的错误。

这样一来， S 因子的模型就建立好了。点击运行

运行成功口在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 因子模型的时候是一样的，这里不再赘述。下面解释一下代码：

32 64

128

16 1

8

4 2

对于中间栅格来说， 它与邻近的

8 个栅格中心点之间的距离只有两种 （ sqrt(2) 和 1），当流向

为 32,128,8,2 的时候，距离为 sqrt(2) ，其他情况下距离为 1。因此在代码中的表现就是如上的形式。

接下来建立 L

( l )m 中计算 m 的模型

22.13

这个过程有点和计算 S 因子类似，不同的是计算 S 因子使用的是度，而计算

m 使用的是百

分比。依据的公式为：

0.5 5%

0.4 3%

5% 为像元坡度 (%)

m

1%

式中，

0.3

3%

0.2 1%

将 Slope 工具和 Single Output Map Algebra 工具工具拖入编辑窗口中， 然后将他们连接起来。双击 Slope(2)，进行如下设置：

这里要选择Persent 百分比的形式，因为m 的判断是利用百分比形式的坡度。

然后在 Single Output Map Algebra 中输入的代码为：

CON([Slope_p] >= 0.05,0.5,CON([Slope_p] < 0.05 & [Slope_p] >= 0.03,0.4,CON([Slope_p] >= 0.01 &

[Slope_p] < 0.03,0.3,0.2)))

注意：这里使用0.05 来表示5%，而不要将代码写成[Slope_p] >= 5 。因为在变量存储值的时候，是不会存储 %的，它只能利用浮点型数字来表示百分数。

接下来计算 L (l)m中的l，依据l i i i 1

( D i / cos i )(D i / cos i )D i / cos i

22.1311

式中， l i为像元坡长， D i为沿径流方向每像元坡长的水平投影距(在栅格图像中为两相邻像

元中心距，随方向而异)，θi为每个像元的坡度( °)， i 为自山脊像元至待求像元个数。

从公式中可以看出，我们需要求出 D 和θ才能计算l，而在前面的部分，我们已经算出了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 ”添加到视图中。

到此， L 和 S 因子完成。

(3)降雨侵蚀力因子（ R 因子）

利用日降雨量估算降雨侵蚀力的多参数模型来计算流域的降雨侵蚀力，公式如下：

R = - 0.0334 P + 0.006661 P2(1)

式中 R表示的侵蚀力值( MJ ·mm·hm-2·h-1 )， P表示年雨量（mm）。

为了方便，我们新建一个Model ，命名为“ K 因子”。

在编辑窗口中先后拖入Create Thiessen Polygons 工具， Feature to Raster 工具，还有Single Output Map Algebra 工具，并把它们连接起来

双击 Input Features 输入实验数据。

双击 Create Thiessen Polygons 进行如下设置

双击 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 值估算采用Williams 等在EP IC 模型中的方法, 利用土壤有机质和颗粒组成因子进行估算，计算式如下：

K ={0.2+0.3exp[-0.0256S d(1 - S i/100)]} *[ S i/(C l+

S i)] 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- S d/100)/{1- S d

/100+exp[-5.51+22.9(1- S d/100)]}}*0.1317

式中： Sd 为砂粒含量，Si 为粉粒含量， Cl 为粘粒含量， C 为有机质含量。

Single Output Map Algebra工具，然将 4 个 Feature to Raster 工具拖入编辑窗口，并拖入一个

后将他们按照如下方式连接和命名

双击 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 - [SAND] / 100) / (1 - [SAND] / 100 + exp(- 5.51 + 22.9 * (1 - [SAND] / 100)))))

模型如下

运行这个模型，并将结果添加到视图中

到此， K 因子完成。

（ 5）植被覆盖度因子（ C 因子）和水保措施因子（P 因子）

植被覆盖度因子, 又称作物经营管理因子。经验指出, 植被覆盖度与土壤侵蚀量关系极大。

在其它地理环境因子值相同的情况下, 植被覆盖度越大, 土壤流失量越小；反之, 则越大。

流域的 C 因子值赋值如表。

表 1不同土地利用 C 因子值

土地利用类型旱地水田交通用地和水体草地居民地林地

C 因子0. 310.1800.060.20.006水土保持措施因子是采取水保措施后, 土壤流失量与顺坡种植时的土壤流失量的比值。通常，包含于这一因子中的控制措施有：等高耕作、等高带状种植和修梯田等。将土地利用图与P 值属性库文件记录建立链接，再分别将P 值赋给土地利用图，得到P 值因子图。以自然植

被 P 因子为1，坡耕地为0.35，水稻是梯田修筑最好的一种土地利用，P 值为0.01。

编号土地利用类型

11水田

12旱地

20~23林地

31~33草地

43水体

52居民地

新建一个模型，命名为“CP 因子”

分别将 Feature to Raster 工具和 Single Output Map Algebra 工具拖入编辑窗口，并按如下方式连接它们和命名

双击 Input features 进行数据的输入

双击 Feature to Raster 工具进行如下设置

在Single Output Map Algebra 中输入如下代码：

Con([Landuse] == 11,0.18,Con([Landuse] == 12,0.31,Con([Landuse] >= 20 && [Landuse] <=

23,0.006,Con([Landuse] >= 31 && [Landuse] <= 33,0.06,Con([Landuse] == 52,0.2,0)))))