ArcGIS实验Ex18利用水文分析方法提取山脊山谷线.docx

1、ArcGIS实验Ex18利用水文分析方法提取山脊山谷线第十一章 水文分析练习 1：利用水文分析方法提取山脊、山谷线山脊线、 山谷线是地形特征线， 它们对地形、 地貌具有一定的控制作用。 它们与山顶点、 谷底点以及鞍部点等一起构成了地形及其起伏变化的骨架结构。 因此在数字地形分析中， 山 脊线和山谷线以及地形特征点等的提取和分析是很有必要的。二、目的理解基于 DEM结合水文分析的方法提取出研究区域的山脊线和山谷线的原理；掌握水流 方向、汇流累积量的提取方法以及它们的提取原理； 能将水文分析的方法和其它的空间分析 方法相结合以解决应用问题。三、要求1、利用水文分析思想和工具提取研究区域的山脊线；2

2、、利用水文分析思想和工具提取研究区域的山谷线。四、数据一幅25m分辨率的黄土地貌 DEM数据，数据的区域大概有 140 km2。数据存于 /ChP11/Ex1 中，请将其拷贝到 E： /ChP11/Ex1 。结果数据保存在 /ChP11/Ex1/Result 中。五、算法思想对于水文物理过程研究而言， 由于山脊、 山谷分别表示分水性与汇水性， 山脊线和山谷 线的提取实质上也是分水线与汇水线的提取。 因此， 对于山脊线和山谷线就可以利用水文分 析的方法进行提取。基于 DEM的这种地形表面流水物理模拟分析的原理是：对于山脊线而言，由于它同时也 是分水线， 那么对于分水线上的那些栅格， 由于分水线的

3、性质是水流的起源点， 通过地表径 流模拟计算之后这些栅格的水流方向都应该只具有流出方向而不存在流入方向， 也就是其栅 格的汇流累积量为零。 通过对零值的汇流累积值的栅格的提取， 就可以得到分水线， 也就得 到了山脊线；对于山谷线而言， 由于其具有汇水的性质， 那么对于山谷线的提取， 可以利用 反地形的特点，即是利用一个较大的数值减去原始的 DEM数据，而得到了与原始地形完全相反的地形数据， 也就是原始的 DEM中的山脊变成负地形的山谷， 而原始 DEM中的山谷在负地形 中就变成了山脊，那么，山谷线的提取就可以在负地形中利用提取山脊线的方法进行提取。六、操作步骤1、正负地形的提取(1)启动 Ar

4、cToolbox ，展开 Analysis Tools 工具箱，打开 hydrology 工具集。在图层管 理器中加载研究区域的原始 DEM数据。(2)加载 Spatial Analyst 模块，点击 Spatial Analyst 模块的下拉箭头，点击 neighborhood statistics 菜单工具，利用邻域分析的方法以 11 11的窗口计算平均值，如 图 1。分析结果命名为 meandem，如图 2所示。图1 用邻域分析的方法以 1111的窗口计算平均值图2 DEM数据的平均值(3) 点击 spatial analyst 中的 raster calculator 菜单工具，对原始

5、 DEM数据与邻域分析之后的数据 meandem做减法运算，并将运算结果重分为两级，分级界线为 0，那么大于 0的区域在原始 DEM上就是正地形区域，小于 0的区域在原始 DEM上就是负地形区域。(图 3)图 3 Raster Calculator菜单工具图4 正负地形区域(4) 对上一步得到的二值化数据进行两次重分类，如图 5。一次将正地形区域属性值赋值为 1，负地形区域属性赋值为 0，命名为 zhengdixing ；另一次将正地形区域属性值赋值为 0，负地形区域属性赋值为 1，命名为 fudixing 。分别如图 6、图 7所示。图5 重分类图 6 正地形区域（图中深色区域） 图 7 正

6、地形区域（图中深色区域）2、山脊线的提取：（1） 在ArcMap中加载研究区域的原始 DEM数据。（2） 洼地填充：双击 hydrology 工具集中的 fill 工具，进行原始 DEM的洼地点填充。在 Input surface raster 文本框中选择原始 DEM数据 dem，将输出数据命名为 filldem ，因为选图9洼地点填充的结果(3)基于无洼地的水流方向的计算：双击 hydrology 工具集中的 flow direction 工具， 在 Input surface raster 文本框中选择填充过的无洼地 DEM数据 filldem ，将输出的水流方向(4)数据命名为 flo

7、wdirfill ，如图 10、 11。图 11 水流方向的计算结果(5)汇流累积量的计算：双击 hydrology 工具集中的 flow accumulation 工具。选择flowdirfill 作为输入的水流方向数据；输出数据命名为 flowacc1图 13 汇流累积量的计算结果(6)汇流累积量零值的提取：加载 Spatial Analyst 模块，点击 Spatial Analyst 模块的 下拉箭头，然后单击 raster calculator 菜单，打开栅格计算对话框，在文本框中填写汇流 累积量零值的提取公式： facc0 = (flowacc = 0) ，然后点击 evaluat

8、e 进行计算。计算结果为所有的汇流累积量为 0的栅格。（6） 在 ArcMap中打开 facc0 ，发现所提取出的栅格很乱，有很多的地方并不是山脊线的 位置， 因此应对这个数据进行处理。 处理过程可以利用邻域分析的方法， 对提取出的汇流累 积量等于零值的数据进行 33邻域分析进行光滑处理，处理后的数据如图 16所示。图 16 对提取出的数据进行 33邻域分析进行光滑处理(7)单击 spatial analyst 模块中的 surfer analyst 中的 countline 和 hillshade 菜单命令，分别生成原始 DEM的等值线图 ctour （图 17）和晕渲图 hillshade

9、 （图 18）。图 17 生成原始 DEM的等值线图 ctour(8)打开 neiborfacc0 数据属性信息， 进行重新分类， 将分类级别设置为两类， 不断调 整分界数据大小，并以由 DEM生成的等值线图和晕渲图为辅助判断数据。在 neiborfacc0 中，属性值越接近于 1的栅格越有可能是山脊线，这里确定的分界阈值为 0.5541 ，如图 19、 20。(9)将进行过二值化的 neiborfacc0 进行重分类为 reneibor ，将属性值接近 1的那一类的图 19 neiborfacc0 数据进行重新分类图 20 重分类的结果(10)将重分类过后的 neiborfacc0 数据与正

10、地形数据 zhengdixing 利用 spatial analyst 菜单下的 raster calculator 进行相乘运算，这样就消除了那些存在在负地形区域中的错误 的山脊线，结果如图 22所示。3、山谷线的提取(1)在ArcMap中加载原始 DEM数据。(2)加载 Spatial Analyst 模块，点击 Spatial Analyst 模块的下拉箭头，点击 options raster calculator 菜单工具，打开栅格计算对话框；在文本框中填写反地形的计算公式： fandem =Abs (dem-2000) ，点击 evaluate 进行计算。计算结果与原始 DEM地形完

11、全相反的反地形数据，如图 24 所示。（3） 反地形计算完毕之后， 山谷线的提取就和山脊线的提取步骤一样的， 直到最终利用 重分类的方法将重新分级的邻域分析后的结果二值化为止。在这里，是不需要对反地形 DEM 进行洼地填充的。 计算过程中的数据名称分别为： 水流方向数据为 flowdirfan（ 图 25） ，汇流 累积数据为flowacc2（ 图26） ，零值汇流累积量提取数据为flowacc0fan（ 图27），对 flowacc0fan 进行均值 3 3邻域分析后的结果数据为 nbfacc0fan（ 图28） ，并将其分级改为两级，分级阈值为 0.65677 ，结果如图 29所示，初次得

12、到山谷线如图 30。图 25 水流方向数据图 26 汇流累积数据图 28 行均值 3 3邻域分析后的结果数据图29 重分类将结果分为两级calculator 进行相乘运算， 这样就消除了那些存在在负地形区域中的错误的山脊线， 如图 31。然后将计算结果进行重分类，所有属性不为 1的栅格属性值赋为 NOD ATA，如图 32。就得到了新的山谷线，比初次的山谷线更为细化，如图 33、 34所示。图 31 存在在负地形区域中错误的山脊线（绿色部分）图 32 重分类将错误的山脊线消除图33 初次得到的山谷线 （ 左 ）与消除错误后的山谷线 （ 右） 对比图 34 计算出的研究区域的山谷线图中深色区域，

13、背景为该区域的晕渲图）七、遇到问题及解决方法1 、遇到问题：在ArcMap中打开 facc0，发现所提取出的栅格很乱， 有很多的地方并不是山脊线的位置， 因此应对这个数据进行处理。如何处理？解决方法：将重分类过后的数据与负地形数据 fudixing 利用 spatial analyst 菜单下的 raster calculator 进行相乘运算，这样就消除了那些存在在负地形区域中的错误的山脊线。然后 将计算结果进行重分类，所有属性不为 1 的栅格属性值赋为 NOD ATA，如图 32。就得到了新 的山谷线，比初次的山谷线更为细化。2、遇到问题： 打开 neiborfacc0 数据的属性信息，进行重新分类，将分类级别设置为两类，不断的调 整分界数据大小，并以由 DEM 生成的等值线图和晕渲图为辅助判断数据。怎么做？ 解决方法：在 neiborfacc0 中，属性值越接近于 1 的栅格越有可能是山脊线的位置，打开 Classification 对话框，不断的调整分界数据大小（如左图中的红圈部分） ，直到生成的数0.5541 。据清晰度达到最高，这里最终确定的分界阈值为