刘 远,周买春
(华南农业大学水利与土木工程学院,广州 510642)
数字高程模型(Digital Elevation Model, DEM)是地表曲面或粗糙度的数字化表达,是计算机地形分析基础数据。像过去地形图一样,世界上许多国家都把DEM数据作为一种国家基础设施来建设,就像道路交通、通讯设施一样,通过政府立项,由具有技术力量的科技机构竞标,按一定标准生成,委派专门机构管理和负责更新,并鼓励广泛地推广应用。随着DEM数据的丰富和社会需求迫切,一种趋势是开放DEM数据供公众免费下载,如几种精度不一、覆盖全球可免费下载的DEM有HYDRO1K[1]、SRTM3[2]和ASTER GDEM[3]。不同的DEM数据产品由不同的机构发布,其数据基础和生成方式的不同隐藏了数据质量的差异。
DEM数据在流域水文研究中有广泛的应用。在流域地表水、地下水运动中,地形代表重力对水的作用,使得流域径流随地形起伏从高处流向低处,许多分布式水文模型基于DEM数据,通过提取的数字河网及相关地形信息,将流域产流汇出流域。除流域汇流外,TOPMODEL模型还通过流域土壤—地形指数的分布来模拟流域土壤蓄满状况及径流形成区域,DEM数据为土壤—地形指数计算提供了极大的方便和效率。此外,在利用站网气象数据(如降雨、大气压强、气温等)通过插值形成流域空间分布数据时,也常用DEM代表的地面高程对插值结果进行修正,如在降雪融雪模型、土壤冻融模型和冰川水文模拟中,近地气温需考虑随高程的直减率。DEM数据质量对流域水文研究成果的可靠性和精度有极大的影响,如笔者根据HYDRO1K、SRTM3和ASTER GDEM 3种DEM数据提取的韩江流域数字河网则存在明显差异[4]。
在DEM易于获得的情况下,有关DEM提取地形信息的研究主要集中在信息提取方法[5-9]和水平分辨率对提取信息精度的影响[10-14]2大方面。然而,地形信息的可靠性和精度归根结底取决于DEM是否精确地反映了地面高程,基于不同DEM产品提取流域水文地形信息进行比较研究,国内外还鲜有报道。为此,本文继续深入比较HYDRO1K、SRTM3(version 4)和ASTER GDEM 提取的流域水文地形信息,分析不同DEM数据源对流域地形特征描述的差异,确定最真实反映流域水文地形状况的DEM数据,为流域水文研究选择恰当的DEM数据提供依据。
流域水文地形信息采用分布式流域水文模型BTOPMC[15-17]的地形子模型来提取。BTOPMC的地形子模型由DEM填洼运算、生成数字河网、划分子流域及其他提取地形相关信息(如水流方向、坡度、地形指数等)的算法等构成。
(1)填洼。DEM实际上是通过某种数学插值方法得到的栅格高程,是连续地面的一种离散表达,生成DEM的高程数据源密度和数学插值方法都会造成DEM中或多或少存在某些虚假洼地,即一些栅格上的水流流不出流域出口(对于一个闭合的外流流域,不存在真实洼地),这些洼地可认为是DEM数据的一种误差。要提取一个完整的流域数字河网,首先必须对原DEM进行洼地处理,即通过修改洼地DEM栅格的高程值,使它高于与之相邻的8个栅格其中之一,使得洼地栅格的水流能够顺畅地流出。由于一个(或一群)洼地栅格可能套合在一个更大的洼地中,填洼通常是一个反复的过程,使全部的洼地得以消除,它是地形子模型中时间消耗最长的步骤。BTOPMC采用逐步填洼,直到全部洼地填平,循环终止,每次填洼采用随地貌倾斜度变化的小高程增量:
(1)
式中:i、j分别是洼地的行坐标和列坐标;hc是用户定义的一个小高程增量阈值(避免造成新的洼地),取hc=0.1 m(考虑到HYDRO1K、SRTM3和ASTER GDEM的高程精度是米);α、β是x、y向地形倾角的权重系数,考虑韩江流域的实际情况取α=β=1;Nr、Nc是填洼前DEM的行、列数;填洼后,洼地高程h(i,j)=hmin(i,j)+dh(i,j);hmin(i,j)是与洼地栅格(i,j)相邻的最低栅格的高程。
(3)上游集水面积和流域边界。上游集水面积是指水流流入当前栅格的所有单元格面积(包括当前栅格)。计算任一个单元格的上游集水面积,只需沿着水流方向反向追踪,直至到达流域边界(上游集水面积为1的栅格),当整个水流方向矩阵搜索完毕后,位于追踪线路上的全部栅格数目即为上游集水面积(以栅格的个数表示)。确定流域边界首先要确定流域的出口,即上游集水面积最大的栅格。水流最终流入到流域出口栅格的所有单元格都处于该流域范围,这些栅格所占区域的边界即为该流域边界。
韩江流域是广东除珠江以外的第2大流域,流域面积30 112 km2,位于粤东、闽西南,经纬度为115°13′~117°09′E、23°17′~26°05′N(见图1)。韩江的主源是梅江,发源于紫金县乌突山七星岽,在大埔县三河坝与发源于宁化县武夷山南段木马山北坡的汀江汇合(三河坝以下始称韩江),最后经韩江三角洲,分北、东、西溪在汕头市出南海,全长470 km。韩江流域以山地为主,约占流域总面积的70%。此外,丘陵和平原大约为25%和5%,丘陵主要分别分布在梅江流域和其他干支流谷地,平原主要分布在韩江三角洲。
图1 韩江流域地理位置和水系Fig.1 The geographic location and drainage network of Hanjiang River basin
2.2.1 3种常用DEM数据的介绍
HYDRO1K[1]是美国地质勘测局地球资源观测系统数据中心(US Geological Survey’s EROS Data Center)与联合国环境计划全球资源信息中心(United Nations Environmental Program/Global Resources Information Database)1996年发布的1 km分辨率的全球DEM(不包括格陵兰岛和极地)。HYDRO1K是在另一个基于多源数据的全球DEM GTOPO30[18](分辨率30″,约1 km)的基础上,结合DCW(Digital Chart of the World)的流域河网和边界,增强了水文特性。
SRTM[2](Shuttle Radar Topography Mission)是在由美国国家航空航天局(National Aeronautics and Space Administration,NASA)和国家空间信息情报局(National Geospatial-Intelligence Agency,NGA)合作完成。2000年2月,通过装载干涉雷达于“奋进”号航天飞机在太空飞行11 d,获得地球北纬60°至南纬56°之间的雷达影像数据,覆盖地球80%以上的陆地表面。SRTM于2003年6月发布,按精度可以分为SRTM1和SRTM3,分辨率分别为1″(约30 m)和3″(约90 m),前者仅有北美地区数据,后者覆盖全球。经过多版的改进,SRTM已更新至第4版。
ASTER GDEM[3](Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model)是由美国国家航空航天局和日本经济产业省于2009年6月发布的30m分辨率的全球DEM。ASTER GDEM覆盖范围占全球陆地面积的99%(北纬83°至南纬83°),是目前覆盖最广的地形数据。该数据是根据NASA的新一代对地观测卫星TERRA的详尽观测结果制作而成。从2000年发射至今,ASTER还在不断更新数据,算法也不断改进,目前已发布了第2版ASTER GDEM数据。
2.2.2 DEM数据处理
用韩江流域所在的矩形边界分别裁剪3种DEM数据,见图2。3种DEM反映的地形变化趋势基本一致,流域地势高低分明,高程较低的山谷、河流清晰可见;具有高分辨率的ASTER GDEM和SRTM3更能反映地形的精细度。统计它们的特征值,结果见表1。同时,统计各个高程值的栅格出现的频率,得到累积分布曲线(见图3)。表1显示,HYDRO1K、SRTM3和ASTER GDEM在矩形边界内的平均值很接近,但DEM分辨率越低,最大高程越小,最小高程越大,即DEM的变化范围越小,反映地形起伏的栅格高程的标准差也越小。由图3可以看到,3条高程累积分布曲线很接近,表明3种DEM反映的研究区域内的地貌变化是一致的,特别是SRTM3和ASTER GDEM;HYDRO1K中,大高程栅格的比例较SRTM3和ASTER GDEM小,说明较低的分辨率使得DEM反映的地形趋于平坦。
图2 研究区域3种DEM的空间分布Fig.2 Spatial distributions of 3 DEMs in the study region
图3 研究区域的高程累积分布曲线Fig.3 Elevation frequency curve of the study region
表1 研究区域3种DEM的基本特征 m
利用BTOPMC地形子模型对裁剪的研究区域内3种DEM的原始数据进行填洼处理,生成各自的无洼地DEM,用于计算流域的水文地形信息。
在填洼后的DEM基础上,采用D8法确定研究区域各栅格的水流方向[见图4(d)],并裁剪韩江流域边界(潮安水文站以上)内的部分,所用的流域边界是由各自DEM提取得到的,结果见图4(a)~(c)。3种DEM提取的水流方向栅格在空间上分布基本一致,分辨率较高的SRTM3和ASTER GDEM吻合程度相对较高;HYDRO1K在汀江中、上游平原地区和ASTER GDEM在梅江的中游出现了较大面积的同流向区域,在这些地方提取的河网出现与实际不相符的平行河道或错误河道[4]。这与采用的D8法有关,根据D8法对流向的定义,其计算得到的流向往往会在 45°或其倍数的方向上产生平行径流。
图4 3种DEM提取的韩江流域水流方向的空间分布Fig.4 Spatial distributions of drainage directions in Hanjiang river basin extracted from 3 DEMs
由图5的对比可以看到,3组数据8个方向的比例基本一致,SRTM3和ASTER GDEM更为接近;以水平和垂直方向居多,HYDRO1K最为明显,它在该4个方向的比例达到66%,又以北向最少;ASTER GDEM虽也以水平和垂直方向居多,但8个方向的比例较其他2组数据均匀。这同样与采用的D8法有关,因为沿水平或垂直方向的路径比沿对角方向短,在相邻8个栅格的高程相差不大的情况下,根据最陡坡度法则确定的流向多为水平或垂直方向;另一方面,DEM分辨率越低,栅格尺寸越大,沿水平或垂直方向的路径与沿对角方向的路径就相差越大,但是低分辨率DEM使地势趋于平坦,相邻栅格间的高程变化变小,所以DEM分辨率越低就越容易产生水平或垂直的流向。
图5 3种DEM提取的韩江流域水流方向对比Fig.5 Comparisons of drainage directions in Hanjiang river basin extracted from 3 DEMs
在填洼后的DEM基础上,沿水流方向计算流域内各栅格的坡度,结果见图6。3种DEM提取的坡度栅格在空间上分布基本一致,高山、平地清晰能辨;分辨率较高的SRTM3和ASTER GDEM的吻合程度相对较高,反映的地面坡度变化更有层次,地形更加精微;分辨率较低的HYDRO1K反映的地形较为
单调,呈现大面积的平坦区域,与实际地形有所不符。3组数据的统计特征见表2,HYDRO1K的平均值和标准差都比SRTM3和ASTER GDEM小很多,说明低分辨率DEM使得韩江流域的地形趋于平坦、均匀。
图7是根据3组坡度数据统计得到的各级坡度栅格出现频率的累积分布曲线。HYDRO1K的曲线明显偏离SRTM3和ASTER GDEM的曲线,其坡角大于3°的栅格百分比还不到50%,严重坦化了流域的地形,将影响由坡度计算的其他地形因子(如地形指数)的精度,同时也影响由HYDRO1K提取的数字河网的精度,其数字河网在坡度较小的成片区域中生成了与实际不相符的平行河道[4]。SRTM3和ASTER GDEM的坡度累积频率曲线吻合较好,说明它们反映的韩江流域地势具有较好一致性;ASTER GDEM的曲线呈跃阶分布,说明ASTER GDEM虽具有较高的水平分辨率,但垂直精度不如SRTM3(前者垂直精度±20 m,后者±16 m),导致坡度的跃阶。
图6 3种DEM提取的韩江流域坡度的空间分布Fig.6 Spatial distributions of surface slope in Hanjiang river basin extracted from 3 DEMs
图7 3种DEM提取的韩江流域坡度的累积分布曲线Fig.7 Surface slope distribution curve of Hanjiang river basin extracted from 3 DEMs
采用单流向算法计算韩江流域的地形指数,结果见图8。3种DEM提取的地形指数栅格在空间上分布基本一致,地形指数较大区域,即土壤缺水容量小或土壤达到了饱和状态,能呈现出流域的河流、沟谷形状;分辨率较高的SRTM3和ASTER GDEM吻合程度相对较高,其地形指数的空间变化更有层次,即更能反映流域土壤水分空间变化,符合山地、丘陵流域的实际情况,这与高分辨率DEM更能反映流域地形变化的结果是一致的。3组数据的统计特征见表2,由于HYDRO1K分辨率低使得地形区域均坦化,流域的坡度变小,所以其地形指数在全流域的平均值和标准差都明显较其他2者大,据此进行流域水文模拟时,流域土壤水分易趋于饱和。
图9是根据3组地形指数数据统计得到的各级地形指数栅格出现频率的累积分布曲线,HYDRO1K的分布曲线与另外2条曲线相差较大,其地形指数大小分布相对均匀,而SRTM3和ASTER GDEM的地形指数集中分布在10以下,分别达到了82.3%和75.0%,说明流域大部分区域的坡度较大,土壤水分不易于饱和,与流域的山地、丘陵地形相一致;SRTM3与ASTER GDEM的分布曲线吻合较好,在地形指数较小处,几乎不存在差别,虽然ASTER GDEM的坡度分布曲线呈明显的阶梯分布,但其地形指数已几乎不存在这一现象,说明由DEM垂直精度引起的误差,在地形因子的计算过程中得到一定程度的衰减。
图8 3种DEM提取的韩江流域地形指数的空间分布Fig.8 Spatial distributions of topographic index in Hanjiang river basin extracted from 3 DEMs
DEM数据坡度平均值最大值最小值标准差地形指数平均值最大值最小值标准差HYDRO1K0.0720.6108.600×10-80.07312.60433.4097.4025.115SRTM30.2422.5692.400×10-80.1838.07037.4883.5473.639ASTERGDEM0.2384.3841.330×10-70.1908.66937.0072.6163.750
图9 3种DEM提取的韩江流域地形指数的累积分布曲线Fig.9 Topographic index distribution curve of Hanjiang river basin extracted from 3 DEMs
笔者[17]对HYDRO1K、SRTM3和ASTER GDEM提取的韩江流域数字河网的对比分析表明,DEM的垂直精度对提取的河网精度起控制作用,所以ASTER GDEM提取的河网精度不如SRTM3;低分辨率的HYDRO1K在大尺度流域上能提取到一定精度的河网,但在流域地势平坦区域的河网精度较低。进一步对比分析这3种DEM数据提取的地形信息,得到与之相辅相成的结论:①在大尺度的韩江流域上,低分辨率的HYDRO1K提取的地形信息在一定程度上能反映流域的地形变化,但坦化了流域地形,使得基于地形指数的流域土壤水分分布趋于均匀,在产流表现上易趋于饱和,从水文意义上达不到刻画像韩江这一尺度流域的山地、丘陵地形。②DEM的水平分辨率和垂直精度决定DEM反映地表变化的精细度,DEM水平分辨率和垂直精度越高,反映的地形越精微。SRTM3和ASTER GDEM提取的韩江流域地形信息能精微地反映流域的地形变化和土壤水分空间格局变化。③当水平分辨率达到一定精度时,与水平分辨率相比,DEM的垂直精度对提取地形信息的真实性起控制作用。ASTER GDEM的垂直精度 不如SRTM3,所以提取的地形信息不如SRTM3精确。④地形信息提取精度与DEM的填洼算法和流向算法有关。BTOPMC采用随地貌倾斜度变化的小高程增量逐步填洼,填洼过程对区域地形信息的扰动较小;在采用D8单流向算法确定栅格水流方向时,沿最陡坡度流向比沿最大落差流向更容易产生水平和垂直流向,这在低分辨率的HYDRO1K中尤为明显;但由于HYDRO1KDEM生成时融合了DCW河网信息,2种流向算法对数字河网和地形指数的影响只限于那些不具有DCW河网信息的源头区域。
不同DEM数据源提取的流域地形信息存在差别,其反映的流域地形变化、土壤水分空间格局也不同,将对分布式水文模型的产、汇流模拟产生影响。所以,流域水文模拟应当选择最能反映流域真实情况的DEM数据。根据HYDRO1K、SRTM3和ASTER GDEM 3种常用DEM提取的流域地形信息和数字河网的对比分析,就韩江流域而言,选择SRTM3进行水文模拟是最合适的。虽HYDRO1K的精度较低,但由于融合了DCW河网信息使得水文特性得到加强,所以提取的韩江流域地形信息和数字河网仍具有相当的精度,可用于流域的水文模拟。可以设想,若在SRTM3中融入河网信息(DCW,或较之精度更高的河网),必将进一步提高其水文应用的精度。
□
[1] U. S. Department of the Interior, U. S. Geological Survey. Hydro1k[EB/OL]. (2015-01-30)[2015-07-16]. https:∥ lta.cr.usgs.gov/HYDRO1K.
[2] Cgiar-Csi. SRTM 90 m digital elevation data[EB/OL]. (2008-08-19)[2015-07-16]. http:∥srtm.csi.cgiar.org.
[3] Systems Japan Space. Aster gdem[EB/OL]. (2012-12-30)[2015-07-16]. http:∥www.jspacesystems.or.jp/ersdac/ GDEM/E/index.html.
[4] 刘 远, 周买春, 陈芷菁, 等. 基于不同DEM数据源的数字河网提取对比分析——以韩江流域为例[J]. 地理科学, 2012,32(9):1 112-1 118.
[5] Quinn P, Bever K, Chevalier V, et al. The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models[J]. Hydrological Processes, 1991,5(1):59-79.
[6] 解河海, 黄国如. 地形指数若干计算方法探讨[J]. 河海大学学报(自然科学版), 2006,34(1):46-50.
[7] 邬 伦, 汪大明, 张 毅. 基于DEM的水流方向算法研究[J]. 中国图象图形学报, 2006,11(7):998-1 003.
[8] 雍 斌, 张万昌, 陈艳华. TOPMODEL中地形指数ln(a/tanβ)的新算法[J]. 地理研究, 2007,26(1):37-45.
[9] O'Callaghan J F, Mark D M. The extraction of drainage networks from digital elevation data[J]. Computer Vision, Graphics, and Image Processing, 1984,28(3):323-344.
[10] 杨 光, 李庆和, 孙保平. 基于GIS的不同精度DEM提取地形指标误差分析——以盐池县南部山区为例[J]. 水土保持研究, 2008,15(5):57-60.
[11] Wolock D M, Price C V. Effects of digital elevation model map scale and data resolution on a topography-based watershed model[J]. Water Resource Research, 1994,30(11):3 041-3 052.
[12] 秦福来, 王晓燕, 窦培谦, 等. 不同精度DEM流域地形特征提取分析及坡度误差控制研究[J]. 国土资源遥感, 2005,(4):56-59.
[13] 吴险峰, 刘昌明, 王中根. 栅格DEM的水平分辨率对流域特征的影响分析[J]. 自然资源学报, 2003,18(2):148-154.
[14] 张微微, 武 伟, 刘洪斌. 不同比例尺DEM提取地形信息的比较分析——以重庆市为例[J]. 西南大学学报(自然科学版), 2007,29(7):153-157.
[15] Takeuchi K, Ao T. Introduction of block-wise use of TOPMODEL and Muskingum-Cunge method for the hydro-environment simulation of a large ungauged basin[J]. Hydrological Sciences Journal, 1999,44(4):633-646.
[16] Ao T. Development of a Distributed Hydrological Model for large river basins and its application to Southeast Asian Rivers[D]. Kofu, Japan: University of Yamanashi, 2001.
[17] Takeuchi K, Hapuarachchi P, Zhou M, et al. A BTOP model to extend TOPMODEL for distributed hydrological simulation of large basins[J]. Hydrological Processes, 2008,22(17):3 236-3 251.
[18] U. S. Department of The Interior, Survey U S G. Global 30 arc-second elevation (GTOPO30)[EB/OL]. (2015-01-01) [2015-07-16]. https:∥lta.cr.usgs.gov/GTOPO30.