做全国尺度的地形分析项目,我手里压了几个不同来源的DEM数据,实际跑下来以后,长期锁定的数据源还是它:优化后的哥白尼中国30米高精度DEM。这套数据不是官方给你一捆原始瓦片直接砸过来,而是已经按中国区域做过水域调平、空洞填补、坐标系统一和分幅整理。对于做防洪模拟、选址规划、坡度坡向提取的同行来说,省掉的不只是一周的数据清洗时间,还有一堆只有踩过坑才会懂的麻烦。这篇文章就把我选它的理由、数据处理细节、实操流程和避坑经验完整说一遍。
1. 为什么首选锁定哥白尼30米,而不是其它DEM
1.1 哥白尼30米到底什么来头
很多刚接触DEM的朋友会把它和SRTM混在一起。哥白尼30米高程数据(一般叫GLO-30)来自雷达干涉测量,用双星编队飞行,一个发射雷达波、一个接收回波,通过相位差反演出高精度地形。这个生产方式决定了它有几个底层优势:数据一致性很高,不像光学立体像对那样受云和阴影影响明显;全球范围基本覆盖完整,在中国大陆区域几乎没有SRTM那种明显的空白带。
它的网格间隔是按经纬度1弧秒划分的,赤道附近大约对应30米实际距离,所以行业里直接叫30米DEM。官方产品的垂直精度标称在米级以内,很多局部试验和机载雷达测深对比下来,平原地区表现尤其稳。对于需要全国统一标准的项目,这套数据在一致性上是第一梯队。
1.2 和SRTM、ASTER、ALOS放一起比,差距在哪
我整理了一张表格,把所有常见30米级DEM放在同一个尺度下对比,方便你做选型判断。
| 对比项 | 哥白尼30米(GLO-30) | SRTM 1弧秒 | ASTER GDEM v3 | ALOS AW3D30 |
|---|---|---|---|---|
| 原始分辨率 | 1弧秒约30米 | 1弧秒约30米 | 1弧秒约30米 | 1弧秒约30米 |
| 覆盖范围 | 全球,极区也有 | 北纬60°到南纬56° | 全球 | 全球 |
| 山区空洞 | 较少,局部峡谷有空洞 | 高山峡谷空洞明显 | 噪声多、局部有异常 | 多版本拼接后改善 |
| 水体处理 | 带内陆水体识别标志 | 基本未处理水体 | 水面噪声明显 | 有水体修正 |
| 数据来源 | 双星雷达干涉 | 单星雷达干涉 | 光学立体像对 | 雷达干涉+光学修正 |
| 国内使用体验 | 地形连续性好,整体稳定 | 经典稳定,但细节略粗 | 异常尖峰多,需大量清洗 | 细节不错,处理成本高 |
这个表格不是概念对比,是实际项目感受。SRTM当然够经典,但它的原始数据在陡峭山区会留下比较明显的雷达阴影空洞,需要额外插值补救;ASTER GDEM的高程噪声在平原和裸地上非常扎眼,经常出现一米左右的起伏噪声,做坡度分析时会被误判成微地形;ALOS数据细节表现很好,可是要把全国范围的各分幅拼到一起,版本差异带来的接边问题让人头疼,处理成本明显偏高。
哥白尼这套数据最让我放心的点是它的高程连续性。同一座山、同一段河谷,从南到北跑过去,高程不会突然跳变,也不会因为分幅不同出现台阶。做全国尺度分析时,这个特性比某一块区域的极高精度更值钱。
1.3 原始版本的问题,正好是优化版的卖点
如果直接下载官方发布的全球原始分幅,你会发现它并不能拿来做“开箱即用”的国内项目,主要问题集中在四个方面。
第一,内陆水体的高程没有统一压平。湖泊表面本来应该是近乎平面的,但原始数据里湖面会有一层“雷达波纹”,高高低低很影响周边坡度和径流分析。第二,无数据区域和无效值很零碎。有些雷达阴影覆盖的峡谷、冰缘区,高程值缺失,直接参与计算会把结果带偏。第三,全球原始分幅按经纬度1度切分,中国边界附近的瓦片会包含邻国地形,裁剪时很容易带进来多余的地物高程。第四,坐标系是WGS84经纬度网格,拿来和国内常见的CGCS2000成果数据叠加时还要做一轮转换,分幅数量又多,转换过程容易出错。
优化版做的事情,通俗讲就是出厂前帮你做了一次全面质检和修整:水域调平、空洞填补、坐标系统一、边界裁剪干净、无效值重新编码。这些处理不改变DEM的原始地形骨架,只消除掉那些影响分析结果的“毛刺”。这也是为什么我更倾向于推荐优化后的版本,而不是让大家自己去官网一格格拼接。
2. 优化版到底动了哪些“刀”
2.1 负值高程处理:真负值保留,假负值统一
拿到DEM以后,很多人第一反应是看到负值就条件反射地认为是无数据。这个坑我见过不止一次。吐鲁番盆地的高程在海拔以下,最低处大约负154米,这是真实地形。还有不少内陆洼地,也是低于海平面的,比如柴达木盆地内部的局部低点。处理这类负值时,绝对不能直接“负值置零”或者“负值设无数据”,否则整个盆地的地形骨架就废了。
优化版在做负值处理时有两条标准:海洋区域统一处理成0高程,并把海岸线附近的高程过渡处理干净;内陆真实负值区域保留原始负高程,同时把由于雷达噪声造成的零下小波动做平滑修正。这样拿到数据后,你不用再去区分“这个负值是真实洼地还是数据错误”,非常省心。
2.2 水体调平与流域一致性
水面数据处理是整个优化里最见功夫的一步。雷达信号打到平静水面上时,大部分能量会被镜面反射掉,回波很弱,导致水面高程往往是一堆不可靠的噪声。不同时间、不同轨道测得的水面,还可能相差一两米。做洪水淹没分析时,一两米的误差不是小问题,它会直接影响淹没范围边界的判断。
优化版的做法是先把内陆水体矢量边界叠加到DEM上,识别出湖泊、水库、稳定河道等水面区域,然后取水体边界周围可靠的地面高程和水体出口高程做参考,把这些区域统一压平到同一个合理的高程面上。处理好以后,水面周边再按自然地形过渡,这样生成的坡度在水陆交界处是平滑且连续的,不会出现一圈突兀的“堤坝”假象。
我用这套数据做过一次水库下游的淹没模拟,和实测水面线对比,效果比直接用原始数据好了不少。关键就在于水域边界那几行像素的高程靠谱了,后面所有的径流方向计算才靠谱。
2.3 空洞填补:插值不是万能,但局部插值够用
雷达地形数据在极端陡峭的峡谷、高山阴影区还是会有空洞。这些空洞的面积通常不大,可能只是几个像素到几百个像素的破碎区域,但视觉上是一块块黑斑,计算坡度时就是一圈异常值。
优化版采取的填补策略很朴素有效:对每一个空洞区域,用周边可靠高程做局部插值,而不是把整个中国区域丢给一个全局插值算法去拟合。局部插值的地方地形连续性保持得好,不会出现大范围的过平滑。实际操作中,我倾向于先看空洞面积占整个瓦片的比例,如果小于1%,插值处理完全没问题;如果空洞连成片且面积很大,那就要评估这个区域是不是极端地形区,后续使用时要谨慎对待插值结果。
2.4 坐标系与分幅规则:拿到手就能直接对接
官方全球分幅是按WGS84经纬度网格铺开的,而国内日常成果数据多用CGCS2000坐标系,两套坐标之间只做平移转换,效率很低。优化版在做后处理时,已经把数据统一到了CGCS2000经纬度框架下,高程单位统一为米,无效值统一设置为-9999。
分幅规则也做了调整,不再是全球统一的1度网格,而是按国内常用的图幅方式重新切分。单幅大小控制在几百兆以内,方便按区域下载。每幅数据都有对应的索引文件,下载前可以按经纬度范围快速定位。这个整理动作看起来不大,但在大批量项目中,省下来的时间和避免坐标混乱的意义非常明显。
3. 实操:拿到优化版数据后,怎么正确用起来
3.1 先做一次数据体检
不管数据来源多靠谱,解压以后第一步永远是体检,不要直接拿去出图。体检内容包含范围检查、坐标检查、像素类型检查、NoData设置检查。
把分幅文件解压到同一个目录后,先列目录看看文件数量是否和索引一致。然后随便挑覆盖目标区域的几幅,用命令行工具看影像属性,确认坐标系是CGCS2000、像素深度是16位整型或浮点型、无效值标记是-9999。这一步能让你在早期发现文件损坏、下载不完整的问题,而不是等到出图时才发现中间缺了一块。
# 列出目录下所有tif文件,确认数量和命名规则 ls -la ./tiles/*.tif | head # 查看某一幅的详细元数据,重点看坐标、尺寸、波段、NoData gdalinfo ./tiles/N35E105.tif # 查看整幅高程统计,确认无异常跳变 gdalinfo -stats ./tiles/N35E105.tifgdalinfo -stats输出的最小值如果出现-9999或-32768,说明这块区域有无效值或空洞未填充,需要进一步处理;如果最小值和最大值都落在合理范围内,说明数据状态正常。
3.2 分幅瓦片拼接成虚拟栅格
全国范围分析不需要把几百幅数据物理拼接成一个巨大的TIFF,那样既费硬盘又费内存,还特别容易失败。更稳妥的做法是用虚拟栅格文件把所有分幅串起来处理。虚拟栅格不复制像素,相当于一个索引文件,把所有分幅拼在一起,后面的裁剪、重投影、分析都能直接操作。
# 用所有tif构建虚拟栅格 gdalbuildvrt china_dem.vrt ./tiles/*.tif # 查看虚拟栅格整体范围 gdalinfo china_dem.vrt | head -n 30 # 按目标范围裁剪出一块工作区 gdal_translate -projwin 95 42 105 30 -of GTiff china_dem.vrt work_area.tif-projwin后面的四个数字分别是左上角经度、纬度、右下角经度、纬度。裁剪出来的工作区尺寸按需控制,做局部试验时尽量控制在几万乘几万的像素以内,避免后续计算卡顿。
3.3 重投影与插值方法选择
优化版数据一般已经统一到CGCS2000经纬度,但经纬度坐标不是做面积、坡度计算的理想坐标系。经纬度网格在不同纬度上的地面距离不一致,越往北一个经度跨度对应的实际距离越短。直接在这种坐标下计算坡度,相当于用不等比的尺子去量地形,结果一定偏。
因此,在做正式分析之前,建议先转换成适合项目区域的投影坐标系。常见做法是根据项目覆盖范围选择:全国范围做面积统计,可以选用双标准纬线等积投影,变形较小;省市级做坡度坡向分析,可以直接转到高斯-克吕格对应的分带;小范围流域项目,则用UTM分带投影就足够。
重投影时有一个非常关键的参数:重采样方法。默认的近邻法会把地形像素的细节弄出锯齿,坡度结果也会出现台阶状条纹。做地形分析时建议使用bilinear双线性插值或cubic三次卷积插值。双线性插值相对保守,适合地形平缓区域;三次卷积对地形细节保留更好,但偶尔会产生轻微过冲,在陡峭山区要留意。
# 将虚拟栅格重投影到CGCS2000高斯克吕格3度带带号39,输出LZW压缩GeoTIFF gdalwarp -t_srs "EPSG:4549" -r bilinear -co COMPRESS=LZW -co BIGTIFF=YES china_dem.vrt gk39_dem.tif这里EPSG代号只是一个示例,具体分带代码要根据目标区域所在经度选择。转换完成以后,再用gdalinfo -stats看一眼统计值,确认高程范围没有因重采样引入异常。
3.4 直接出地形成果的三件套
有了标准投影的DEM,最常见的三件事就是山体阴影、坡度和坡向。山体阴影用于目视判读,坡度坡向用于定量分析。在桌面GIS软件里操作都一样,只是需要注意参数设置。
生成山体阴影时,方位角和高度角不一定要用默认值。做展示图时,方位角315度、高度角45度的效果比较自然;做特定分析时,要根据太阳轨迹和地形走向调整。坡度计算一定要先确认数据没有负的无效值混入,否则会出现一圈诡异的低坡度异常。坡向是环形变量,导出成果时注意从0到360度的标准化,不要出现负角度。
生成的坡度图如果出现有规律的横向条纹,八成是重采样方法选成了近邻,或者投影转换阶梯效应没有处理好。这时候回去检查投影步骤,改成三次卷积重新来一遍,条纹基本会消失。
4. 常见问题与排查技巧实录
| 问题现象 | 可能原因 | 解决方法 |
|---|---|---|
| 数据里出现-9999或-32768 | 有效无效值未被识别 | 统一设置NoData为-9999后再计算 |
| 湖泊表面高低不平 | 原始DEM水面波纹未完全压平 | 结合水体矢量边界局部重新调平 |
| 坡度图出现条带纹理 | 重采样方法不当或投影转换变形 | 改用双线性或三次卷积重采样 |
| 山体阴影边缘有锯齿 | 像素深度转换导致细节丢失 | 保留原始16位高程,避免中途转8位 |
| 分幅拼接处高程跳变 | 分幅间坐标基准差异 | 用虚拟栅格处理,避免物理拼接 |
| 重投影后出现黑边 | 无效值被当成0参与了插值 | 先设置NoData,再做gdalwarp |
大多数问题并不是数据本身有问题,而是处理顺序不对。处理DEM有一个铁律:所有清洗、投影、裁剪工作做完以后,才去计算坡度坡向。很多人图省事,直接在原始分幅上切一块就开始算坡度,后面再做投影转换,结果就是各种条带、锯齿、假地形。
排查异常时,我习惯先把出问题的区域单独裁剪出来,做一个山体阴影和等高线叠加对比。如果等高线出现明显直角转弯,说明该区域高程有问题,需要拉回来重新处理;如果只是看起来粗糙,大概率是可视化拉伸方式和数据本身无关。
内存不足是做全国处理最常见的崩溃原因。物理拼接超大TIFF会导致内存爆炸,正确做法是用虚拟栅格管理所有瓦片,再按需裁剪。如果裁剪出来的工作区仍然太大,可以用金字塔降低读取压力,或者按图幅分批处理,最后合并结果。
还有一个小经验:在做流域分析时,最好先根据水流方向对DEM做一次填洼处理,而不是直接拿原始DEM计算。优化版虽然做了水域调平,但地形洼地是自然存在的,不填洼直接计算汇水路径会得到零碎的、断断续续的河网。填洼阈值要综合分析区域地形设置,平地区域阈值可以适当调大,山区尽量用小值,避免把真实地形磨平。
5. 一点个人实操体会
用这个数据的过程中,我最常用的一个流程是在重投影之后顺手做一次局部水域边界的抽稀和平滑核对。具体做法是把湖泊水库的矢量边界叠到DEM上,沿着岸边像元采样一圈高程值,主要看岸边有没有明显的坡坎跳跃。如果岸边一圈高程值极差小于1.5米,说明水域调平做得干净,可以直接进入水文分析;如果落差能到好几米,我会再针对这一片做一次岸边高程平均值替换,把这几个异常像元修正掉。这个习惯帮我避免过很多次因水面噪声导致下游坡度计算失真。
对于刚开始接触DEM项目的读者,我的建议是从省域范围的一个小流域开始练手,把这套下载、体检、拼接、投影、出坡度的流程完整跑通两三遍,再上全国级任务。优化后的数据底子已经很好,但工具链和流程的熟练程度,才是最终成果可靠的关键。