1、 第 1 页 (共 13 页 ) 基于 SRTM-DEM 水系提取及验证方法研究 以汉江为例 摘要 ： 介绍基于地表流漫流模型从 数字高程模型 (DEM)提取水系的 提取方法，以汉江为例提取汉江流域的水系，并 用从汉江 TM影像中提取 的 水系作为验证 依据验证 上述提取水系方法的精确性。 水系提取方法 的主要流程是洼地的识别和填充、水流方向的确定、 汇流累积量的确定和河流栅格网络的生成以及 矢量化 水系的生成； 从 TM影像 提取 的水系为验证依据验证上述提取水系方法的精确性。 研究结果显示，从 DEM中提取水系基本与真实水系吻合，但是从某些局部看提取的水系还是存在偏差 。 关键字 ： DE
2、M数据水系提取； TM数据水系提取；验证 引言 以往的水文水资源的研究和实践工作中，人们 通常 从 地形图或 其它图中 手动 获取所需水系的分布 ，工作量巨大， 耗时费力，加上 人为的水系分级或者忽略低级水系 使 水系提取 标准不 统一 ，提取的水系网又很难做进一步利用 。 随着GIS在 当前 各方面的应用越来越广， 高精度的 DEM 获取渠道逐步拓宽， 从 DEM 中提取 的 水系 满足了人们某些领域的需要 ， GIS中 水文分析模型 (Hydrology model)提供 的 水系提取方法 使水系网能自动生成， 并且生成的水系 趋近真实的水系网 。水系提取精 度在水系提取中是必须考虑的因素
3、，科学合理的验证方法可以促进提取方法的成熟。过去的验证方法中很少使用 遥感 影像来验证水系提取的 质量 ，本研究尝试 使用从 Landsat TM 影像中提取的水系情况作为验证资料 ，这 种验证方法是否可行可靠十分值得探索和研究。 GIS 中 水文分析模型的主要功能是从数字高程模型 （ DEM） 中提取地表水文要素特征信息 ,进行 地表水文的空间 分析 ，并将分析结果 可视化显示 ，能基于一定分辨率 DEM 分析显示出 地形 的 局部特征， 按某种算法 可以自动 的将 一定地理空间范围内的自然水系 提取出来 。 在使用 DEM 数据提取水系以 来，学者们对提取方法的研究层出不穷，其中 一种是
4、O Callaghan 和 Mark 提出的坡面流模拟方法 ， 通过 模拟地表径流在地表的流动来产生水系 3。该算法依据水总是沿斜坡第 2 页 (共 13 页 ) 最陡方向流动的原理 ， 确定 DEM 中每一个栅格单元的水流方向；然后根据栅格单元的水流方向计算每一个栅格单元的上游 汇流量 ， 再选择合适 汇流量 阈值来确定河网 ，这种 方法 可直接生成相互连接的河网 ， 操作简单，过程明确，结果能较好的反映真实水文要素，常常作为 DEM 水系提取时最好的方法 。 近几年里很多学者实践了这种方法，不仅提取出研究区水系网，而且在影响提取的水系网的多个因素和应用提取的水系网求得其他水文特征方面也做出
5、详细的论述 ,他们验证出这种方法虽然仍存在一些缺陷，但能取得其它方法无法比拟的效果 1 5 。 1 研究区概况 研究区 汉江又称汉水，古时曾叫沔水，与长江、黄河、淮河并称“江河淮汉”。汉江 就 长度而言为长江第一大支流 ，它 发源 于 陕西省西南部秦岭与米仓山之间的宁强县（隶属陕西省汉中市，旧称宁羌）冢山，而后向东南穿越秦巴山地的陕南汉中、安康等市，进入鄂西后北过十堰流入丹江水库，出水库后继续向东南流，过襄阳、荆门等市，在武汉市汇入长江。汉江流域涉及鄂、陕、 豫、川、渝、甘 6 省市的 多个地区 。 流域地势西北高，东南低， 西为褶皱隆起中低山区 ， 东以平原丘陵为主。汉江流域属亚热带季风区，
6、气候温和湿润 ，降 水量较丰沛；但年内分配不均， 5 10 月径流量占全年 75%左右，年际变化较大，是长江各大支流中变化最大的河流。汉江干流中下游区，是长江中游重点保护区之一。由于河槽泄洪能力与洪水来量严重不平衡，历史上洪水灾害严重，也是长江支流中洪水灾害最严重的一条。 2 研究方法 2.1 从 DEM数据中提取水系 2.1.1 DEM 数据预处理 SRTM（ Shuttle Radar Topography Mission，即航天飞机雷达地形测绘使命 ）数据是美国 太空总署（ NASA）和国防部国家测绘局（ NIMA） 联测，由 2000年 2 月 11 日美国发射的 奋进 号航天飞机上搭
7、载 的 SRTM 系统 ，花费 222 小时23 分钟 ，采集到 北纬 60度至南纬 56度之间 的 地表面数据 。 雷达影像数据总面第 3 页 (共 13 页 ) 积超过 1.19亿平方公里 ，采集到全球超过 4/5的陆地表面数据。 SRTM 数据每个文件里包含一个经纬方格，精度 1 arc-second 和 3 arc-seconds，称作 SRTM1和 SRTM3，或者称作 30M 和 90M 数据 ， SRTM1文件里 包含 3601*3601 个采样点的高程数据， SRTM3文件里包含 1201*1201 个采样点的高程数据。由于受下载限制权限的影响，高分辨率的数据较少用于普通科研项
8、目中。 目 前中国境内可免费下载的是 90m 分辨率高程数据， 30m 数据 点算术平均得来的 。采集数据 经过两年多时间 进行了处理，生成我们现在可用的 SRTM 地形产品数据 数字地形高程模型。 原始 DEM 数据在存在洼地、坐标系等方面存在问题， 为使数据正确使用，提高计算速度等目的，本研究中对原数据进行了下面 一 系列的处理 。 1.从 STRM Data Search 下载网页中下载 分辨率为 90m 的 STRM DEM 数据，汉江流域一共占四个分区， 为了得到完整的水系网，需将分隔开的四个区域拼接在一起。 将下载数据添加到 ArcMap 中，点开工具箱里 Data Managem
9、ent Tools-Raster-Raster Dataset-Mosaic To New Raster，将四 幅 DEM作为 输入，设置拼接后的文件名及存储位置， 其他参数均为默认，将四幅数据合成一幅 DEM。 2.由于下载数据只有地理坐标系 GCS_WGS_1984，要想与 TM 数据提取水系进行比较就必须在同一坐标系中即 WGS_1984_UTM_Zone_49N，所以要将下载数据进行投影变换。 点开工具箱里 Data Management Tools-Projections and Transformation-Raster-Project Raster，以合成 DEM 为输入数据，设
10、置输出文件名，输出坐标系选中 Import 输入跟原 TM 数据一样的投影坐标系， 坐标系参数参见表 1。 表 1 投影坐标系参数 投影坐标参数名称 参数取值 Projection Transverse Mercator False Easting 500000.000000 False Northing 0.000000 Central Meridian 111.000000 Scale Factor 0.999600 Latitude Of Origin 0.000000 Linear Unit Meter (1.000000) Geographic Coordinate System G
11、CS_WGS_1984 Angular Unit Degree (0.017453292519943299) Prime Meridian Greenwich(0.000000000000000000) 第 4 页 (共 13 页 ) 续表 投影坐标参数名称 参数取值 Datum D_WGS_1984 Spheroid WGS_1984 Semimajor Axis 6378137.000000000000000000 Semiminor Axis 6356752.314245179300000000 Inverse Flattening 298.257223563000030000 3.为降
12、低数据计算量 ,将上述 DEM重采样，降低栅格数据的数量。点开 Data Management Tools-Raster-Raster Processing-Resample，以上述投影变换后的合成 DEM 为输入，设置输出名，输出栅格单元大小设为 500m，以默认重采样方式进行重采样。 4.为进一步降低计算量，缩短计算时间，根据所了解的汉江流域的 大致 范围，可将计算范围缩小，以达到提高计算效率的目的。在 ArcCatalog 里建立一Shapefile 文件，文件名设为掩膜，要素类型为面 Polygon。单击 Edit 在坐标系选项卡里选中 Import 输入跟投影变换后的合成 DEM 数
13、据一样的投影坐标系，坐标系参数 参见 （表 1）。将掩膜文件，投影变换后的合成 DEM 添加到 ArcMap 中，单击 Editor工具条，在其下拉菜单中点击 Start Editing 开始编辑该面文件，根据以往所学知识在投影变换后的合成 DEM 数据上画出汉江流域的范围，以减少计算量，画出 的 结果用 Save Edits 保存编辑 ，点 Stop Editing 停止编辑。点开工具箱里空间分析中 Extraction-Extract by Mask，以投影变换后的合成DEM数据作为输入，掩膜文件为掩膜数据，设置输出结果名，提取出汉江流域范围内 的 DEM。如 图 1所示 。 第 5 页
14、(共 13 页 ) 图 1 汉江流域 大致范围 的 DEM 5 由于存在“数据洼 地” （即数据误差造成的假洼地） ， 水流方向计算时 会出现逆流现象，水系生成时出现水流断流，使水系网不连续 。 水系提取时，首先要对原始 DEM 进行获取和填充，修正 DEM 数据。原理是搜索 DEM 中的每一个栅格点，找出比其周围点的值都小的栅格点，将其标为洼地，然后将其高程值抬升至周围点的最低值将洼地填充。很多研究者认为 DEM 数据的洼地填充与否，填充方法的优劣对水系提取结果有着莫大的影响，为此提出了各种洼地获取与填充的方法，对洼地进行了正确而合理的填充 6 10 。在对汉江的研究中， 应用水文分析模块中
15、 Fill 工具 ， 以原 始汉江 DEM数据作为输入进行洼地填充，其中，洼地填充阀值设置为系统默认即将所有洼地进行填充，并进行 两 次 洼地 填充，得到汉江无洼地 DEM。 用 Spatial Analyst-Raster Calculator 工具，求出原DEM与第一次洼地填充的差值，可以得到原 DEM数据中洼地的分布如图 2 所示，第二次填洼后，求得第一次填洼与第二次填洼的差值为 0，表示研究区洼地已全部被填充，生成无洼地 DEM。 第 6 页 (共 13 页 ) 图 2 洼地分布及高程 2.1.2 水系的提取 水系提取大致可以分为四步： 首先确定汉江水流流向，栅格单元的水流方向是指 水
16、流流出该单元格的方向。水流方向的确定采用 D8 算法 ， 计算出每一中心栅格与其周围栅格的坡度，找出坡度最大的那个栅格，该栅格就是水流方向。 两个相邻单元格 i 和 j 之间的坡度计算公式为 ： j = arctg hi hj D （ 式 1） 式中 ， hi和 hj为两个单元格的高程值 ， D为两单元格中心之间的距离 6 。 再计算汇流量，上游 汇流量 指所有 到达某一栅格的水流面积 。原理是假想在 流域内 的每一网格上降下一单位的水量 ， 按水流方向每移动一个 栅格，汇流量将增加一个单位， 因此 可以计算出每一栅格的上游累计汇流量，因为得出结果为累计流经的栅格单元数量，需将其乘以栅格单元面
17、积才是最后所求汇流量 6 。 再设置汇流量阀值获取水系栅格网络， 阀值设置应以当地的实际情况为基础，综合各方面影响因素，多次尝试设置合适的汇流量值。 最后转换为矢量水系网。 利用 空间分析工具中水文分析模块的 Stream to Feature工具可将水系矢量化。 在研究区汉江的应用中过程： 应用 水文分析模块中 的 水流方向工具 Flow Direction，以无洼地 DEM 为输入得修正后的水流方向，然后应用水文分析模块第 7 页 (共 13 页 ) 中的 流向累积量 Flow Accumulation 工具，以水流方向和无洼地 DEM 作为输入计算出汇流累积量，其他参数设置为默认值即可。
18、结果如图 3所示 ： 图 3 汇流累积量 以 国家基础地理信息系统数据 中的各级水系为参照，分别将汇流量阀值设置为1000、 800、 500、 400、 550 提取水系网，其中阀值为 550 时，提取的水系网与参照水系网最相近，所以设定汇流量的阀值 为 550，在 Spatial Analyst 下拉菜单中选择 Raster Calculator，在对话框内输入 “Flow Acc Flow 550“，得出汇流量值大于 该阀值的栅格水系网，打开空间分析工具中水文分析模块的Stream to Feature 工具 ，以栅格水系网 、 水流方向为输入， 设置文件名， 得到提取的矢量水系网，其中
19、汉江流域 水系网结果如图 4所示 ： 第 8 页 (共 13 页 ) 图 4 汉江 流域 水系网 2 2 从 TM数据水系的提取 LANDSAT是美国陆地探测卫星系统，从 1972年开始发射第一颗卫星 LANDSAT 1，到目前为止已经发射的 LANDSAT 7 的 7 颗卫星中，只有 LANDSAT 5 从 1984年 3 月 1日 发射至今仍在超期运行。 LANDSAT 5 卫星主要轨道参数： 近极近环形太阳同步轨道轨道高度 为 705km， 倾角 为 98.2， 24 小时绕地球 15 圈 ，辐射宽带为 185km， 重访周期 为 16 天 。研究中用于验证的 TM 数据正是 LANDS
20、AT 5 遥感数据生成的数据产品，该数据信息参见 （表 2）， TM 数据空间参考参数参见 （表1）。经过校正过的 TM 数据可以真实的反应地表真实水系网，由于目前还从 TM中直接提取水系的 方法还不成熟 ，所以 研究中 人工手动矢量化影像 来 提取 其 水系网。 表 2 TM 数据相关的信息 信息名称 信息取值 Acquisition Date 2007-9-15 WOID L1200911 第 9 页 (共 13 页 ) 续表 信息名称 信息取值 Path/Row 125 / 037 Lor Reference Image L51IKR1007258020100_HDF.102331334
21、 Origin Image courtesy of the U.S. Geological Survey Product Creation Time 2010-08-21T13:53:58Z Raster Information Columns and Rows 8241,7101 Cell size(X,Y) 30, 30 Format TIFF Pixel Type unsigned integer Extent Top 3776715 Left 428385 Right 675615 Bottom 3563685 Statistics Build Parameters skipped c
22、olumns:1, rows:1, ignored value: Min Max 0 255 Mean 48.70397175183501 Std dev 37.13696697000832 Classes 0 2.2.1创建线文件 在 ArcCatalog 里建立一 Shapefile 文件， 设置文件名 ，要素类型为线Polyline。单击 Edit 在坐标系选项卡里选中 Import 输入与 无洼地 DEM 数据一样的投影坐标系 （参见表 1） 。 2.2.2提取水系 将 创建的线文件 ， TM 数据添加到 ArcMap 中，单击 Editor工具条，在其下拉菜单中点击 Start Ed
23、iting 开始编辑该线文件，对影像进行手动 画出水系网 ，矢量结果用 Save Edits 保存编辑 ，点击 Stop Editing 停止编辑，得出结果即为验证依据。如图 5所示 ： 第 10 页 (共 13 页 ) 图 5 从 TM 中提取 水系网 3 研究结果的 验证 及分析 3.1 用 TM 提取的水系验证 将 从 DEM提取的汉江 矢量 化 水系网， 从 TM 矢量 中提取 水系网加载到 ArcMap中，比较两矢量 化 水系可以看出水系网的基本轮廓大致相同，但在有些部分重合度较好，有些部分却有很大差别。当水系网曲折度较小或水系级别较高时重合的较好，反映了当地真实水系网；在一些拐弯较多的地方 DEM 矢量化水系网较平直，就很难反应真实水系，提取结果的精度也降低，验证结果 为图 6 所示。 3.2 用 TM 影像验证 将 从 DEM提取的 矢量水系网， TM 影像加载到 ArcMap中，可以观察到 DEM提取 的水系网 作为理论化提取，提取结果会出现很多因地形，应用模拟地表径流在地表的流动来产生水系的坡面流模拟方法等原因导致的非真实水系，空间位置大致与影像相吻合，但也有些支流偏离了真实水系位置。通过实验， DEM 数据的水系提取方法受 提取 方法参数的影响，尤其是 DEM 数据填充和汇流量阀值的
