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

空间信息应用实践（中级）实验指导书空间建模——基于RUSLE的土壤侵蚀建模分析一．实验背景SoilerosionandgullyingintheupperPanucobasin,SierraMadreOriental,easternMexico土壤侵蚀是地球表面物质运动的一种自然现象，全球除永冻地区外，均发生不同程度的土壤侵蚀。人类社会出现后，土壤侵蚀成为自然和人为活动共同作用下的一种动态过程，构成了特殊的侵蚀环境背景，并伴随着人类对自然改造能力的增强，逐渐成为当今世界资源和环境可持续发展所面临的重要问题之一。土壤侵蚀被称为“蠕动的灾难”，每年因土壤侵蚀造成的经济损失较诸如滑坡、泥石流和地震等地质灾害更大,土壤侵蚀已成为我国乃至全球的重大环境问题之一。土壤侵蚀及其产生的泥沙使土壤养分流失、土地生产力下降、湖泊淤积、江河堵塞，并造成诸如洪水等自然灾害，泥沙携带的大量营养物和污染物质加剧了水体富营养化，水质恶化，不断严重威胁到人类的生存。据估计全球每年因土壤侵蚀损失300万公顷土地的生产力，造成的损失以百亿美元计。我国人口众多、农耕历史悠久，加之历史上战乱频仍，以黄土高原为代表的华夏文明发源地是世界上土壤侵蚀最严重的区域之一，1990年遥感普查结果，全国水土流失面积达367万km2，占国土总面积的38.2％，其中50%为水蚀地区，土壤侵蚀以黄土高原、四川紫色土地区和华南红壤地区尤为突出，仅黄土高原地区一处，平均每年流失泥沙就达到16.3亿t。水土流失已成为中国重要的环境问题，土壤侵蚀研究已成为目前环境保护中的一个重要课题。土壤侵蚀预报是有效监测水土流失和评价水保措施效益的手段，侵蚀模型则是进行土壤流失监测和预报的重要工具。然而传统预测方法需要在量经费、时间和人力的投入，因此，在一定精度范围内通过有限的数据输入，得到满足要求的土壤侵蚀预测结果成为趋势。80年代以来，随着地理信息系统(GeographicalInformationSystem,GIS)的成熟，它开始与土壤侵蚀模型—通用土壤流失方程(UniversalSoilLossEquation,USLE)相结合进行流域土壤侵蚀量的预测和估算，业已成为土壤侵蚀动态研究的有力工具。GIS与USLE相结合的分布式方法运用GIS的栅格数据分析功能，可预测出每个栅格的土壤侵蚀量，便于管理者识别关键源区，并通过确定引起水土流失的关键因子，针对性地提出最佳管理措施(BestManagementPractices，BMPs)，为流域内土地资源的质量评价、利用规划和经营管理等提供科学依据与决策手段。二、实验目的模型生成器(ModelBuilder)为设计和实现空间处理模型提供了一个图形化的建模环境。模型是以流程图的形式表示，它通过工具将数据串起来以创建高级的功能和流程。你可以将工具和数据集拖动到一个模型中，然后按照有序的步骤把它们连接起来以实现复杂的GIS任务。通过对本次练习达到以下目的：掌握如何在ModelBuilder环境下通过绘制数据处理流程图的方式实现空间分析过程的自动化；掌握土壤侵蚀理论的基本知识；掌握利用脚本文件实现空间建模，加深对地理建模过程的认识，对各种GIS分析工具的用途有深入的理解；在ModelBuilder环境下如何计算RUSLE模型的中各个因子，实现RUSLE模型自动化；三、实验准备实验环境：ArcGISDesktop9.3实验数据：矢量和栅格数据矢量数据：研究区界线(bj.shp)、气象数据（Climate.shp），土地利用数据(landuse_Clip.shp,)和土壤数据(soil_clip)；栅格数据：地形数据（DEM）；四、实验内容与步骤(1)实验准备本次试验需要使用ArcGIS的建模功能，在实验之前需要掌握如何利用ArcGIS进行建模。首先，打开ArcMap，激活工具箱在工具箱中右键单击，选择“NewToolbox”，即可新建一个工具箱。可以在此工具箱上右击，通过“Rename”对工具箱重命名。在新建的工具箱上右击，按照“New”——“Model”新建一个Model，可以按照同样的方法给这个Model命名。然后在此Model上右击，通过“Edit”进入模型的编辑模式。到此，模型准备已经结束，接下来开始逐个建立模型的各个因子。（2）地形因子（L，S因子）算法：坡长因子采用公式计算,mlL)13.22(，式中：L为坡长因子，l为像元坡长，m为坡长指数，像元坡长的计算式如下：iiiiiiiiiDDDlcos/)cos/()cos/(111，m取值如下式：%12.0%31%0.3%53%0.4%50.5m式中，为像元坡度(%)式中，li为像元坡长，Di为沿径流方向每像元坡长的水平投影距(在栅格图像中为两相邻像元中心距，随方向而异)，θi为每个像元的坡度(°)，i为自山脊像元至待求像元个数。坡度因子S分段计算：100.9621.91sin1050.5-16.8sin50.0310.8sinSL和S因子的模型建立：首先在工具箱中找到Resample工具，可以使用工具箱自带的搜索功能快速定位到。在工具箱的下方有一行标签，选择Search标签，在搜索框中输入要查找的工具名，如Resample，点击Search进行查询，查询结束后选中查询结果，点击下方的Locate可以快速定位需要查找的工具。可以将这个工具直接拖到Model的编辑窗口中，如图：现在需要给这个工具添加一个参数，在编辑窗口的Resample上右击，通过“MakeVariable”——“FromParameter”——“InputRaster”添加。注意：这里不建议使用右键菜单的“CreateVariable”来添加输入输出参数，因为很多工具拖入到编辑窗口后会自带一个输出参数，而且它们也有自己的默认输入参数。如果另外新建一个参数，可能会因为这个新建参数类型不与工具要求的输入参数类型对应而出现错误。按照同样的方法拖入Slope工具，SingleOutputMapAlgebra工具。通过编辑窗口上的工具将这些工具首尾连接起来。双击InputRaster，输入dem数据输入数据之后，编辑窗口中的工具颜色会相应的变化，说明这些工具已经相互连接起来，还是白色的工具代表它还没有和前面的工具联系起来构成“流水线”，同时，这也是判断SingleOutputMapAlgebra工具中的脚本语言是否和前面的输出文件关联起来的依据。在相互连接的工具中，只要有一个工具是白色的，就说明这条“流水线”不能正常运行。可以发现SingleOutputMapAlgebra工具还是白的，这是因为我们没有添加算法，下面添加用于计算S因子的算法，依据为：坡度因子S分段计算：100.9621.91sin1050.5-16.8sin50.0310.8sinS双击SingleOutputMapAlgebra工具，添加如下代码：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))这时点击左上方的绿色圆点会发现有错误提示单击会看到SingleOutputMapAlgebra工具所接收到所有参数我们发现，这其中并没有我们计算时所需要的Slope_degree1，所以会提示错误。这时就需要对这个Outputraster参数名进行修改，改为Slope_degree1。修改这个参数其实就是把SingleOutputMapAlgebra工具前的Outputraster输出框进行改名，如图：接下来我们要进行一些参数的设置。（只需要修改红框中的参数，其他的采用默认设置）双击编辑窗口中的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上右击，选择“AddToDisplay”就可以将结果显示出来结果为：接下来对L因子建立模型L因子模型的建立可以在上面的S因子模型基础上进行。需要添加工具Fill和工具FlowDirection以及SingleOutputMapAlgebra工具，然后将他们连接起来SingleOutputMapAlgebra工具中的代码为：Con([FlowDirection]==2|[FlowDirection]==8|[FlowDirection]==32|[FlowDirection]==128,Sqrt(2)*90,1*90)这些步骤和在建立S因子模型的时候是一样的，这里不再赘述。下面解释一下代码：3264128161842对于中间栅格来说，它与邻近的8个栅格中心点之间的距离只有两种（sqrt(2)和1），当流向为32,128,8,2的时候，距离为sqrt(2)，其他情况下距离为1。因此在代码中的表现就是如上的形式。接下来建立mlL)13.22(中计算m的模型这个过程有点和计算S因子类似，不同的是计算S因子使用的是度，而计算m使用的是百分比。依据的公式为：%12.0%31%0.3%53%0.4%50.5m式中，为像元坡度(%)将Slope工具和SingleOutputMapAlgebra工具工具拖入编辑窗口中，然后将他们连接起来。双击Slope(2)，进行如下设置：这里要选择Persent百分比的形式，因为m的判断是利用百分比形式的坡度。然后在SingleOutputMapAlgebra中输入的代码为：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。因为在变量存储值的时候，是不会存储%的，它只能利用浮点型数字来表示百分数。接下来计算mlL)13.22(中的l，依据iiiiiiiiiDDDlcos/)cos/()cos/(111式中，li为像元坡长，Di为沿径流方向每像元坡长的水平投影距(在栅格图像中为两相邻像元中心距，随方向而异)，θi为每个像元的坡度(°)，i为自山脊像元至待求像元个数。从公式中可以看出，我们需要求出D和θ才能计算l，而在前面的部分，我们已经算出了D，θ也已知了，在算S因子的时候就已经算出了θ，因此我们可以再次利用Slope计