GrIMP DEM全解:基于立体摄影测量的格陵兰冰盖数字高程模型
2026/9/15 2:10:58 网站建设 项目流程

做格陵兰冰盖研究的人,手里多少都存过这样一套数据:文件名前缀是GrIMP_DEM,标签上挂着MEaSUREs,格式是GeoTIFF,分辨率30米,覆盖范围从格陵兰最南端的过渡带一路铺到北部皮里地。我第一次认真用它是在处理格陵兰东北部一条潮水冰川的末端变化,当时ArcticDEM还没有铺满整个冰盖,手里的Landsat立体像对又做不了亚像素级的地形改正,我需要一套能稳定覆盖全岛、并且有统一高程基准的面状DEM作为底图。翻来翻去,最后落在MEaSUREs格陵兰冰盖测绘项目(GrIMP)这套基于GeoEye和WorldView卫星影像生产的数字高程模型上,也就是大家常说的GrIMP DEM。

这套数据适合谁?简单说:做冰川变化、冰面高程、冰流速、末端崩解研究的遥感工作者和研究生,以及在极地地区搞地形分析但不想从原始卫星影像重新跑一遍摄影测量流程的人。它解决的核心问题是——如何在没有地面控制点的冰盖上,获得一套空间连续、精度可控、时间上可对比的全岛高程模型。今天这篇就把项目背景、生产原理、版本差异、实操流程和踩坑经验一次性讲透。

1. 先搞清楚这套数据解决什么问题

1.1 格陵兰冰盖高程观测的三个难点

格陵兰冰盖约170万平方公里,差不多是墨西哥面积的两倍。要在这么大面积上持续观测高程变化,地面测量和站点观测全都覆盖不过来。这是第一个难点:观测范围与空间密度之间的矛盾。第二个难点是地形——冰盖边缘不是平缓的斜坡,而是布满陡峭的冰裂隙、狭窄的峡湾和快速流动的冰川。地形起伏大,对遥感立体测高的要求就高。第三个难点是环境限制:高纬度地区极夜漫长,云层频繁,太阳高度角常年偏低,光学遥感影像的获取窗口非常有限。

这三个难点叠加在一起,决定了格陵兰高程观测不能靠单一手段。卫星测高(比如ICESat和ICESat-2)能给出精确的沿轨剖面,但轨道之间的间隔往往几十公里,无法构成连续面。机载激光雷达精度高,但覆盖范围有限,做全岛测量成本极高。剩下可行的大范围连续面状测高手段,就是高分辨率光学立体像对摄影测量,这正是GrIMP DEM的诞生背景。

1.2 GrIMP 项目在整个数据链里的位置

MEaSUREs 全称 Making Earth System data records for Use in Research Environments,是NASA为地球系统科学研究建设长期数据记录的计划。GrIMP 是其中专门针对格陵兰的子项目,全称 Greenland Ice Mapping Project。它产出的不只是DEM,还包括冰流速场、冰川末端位置、冰面高程变化时间序列等一系列产品。

在GrIMP的整个产品体系里,DEM扮演的是“地基”角色。冰流速图反演时需要DEM做透视改正和地形去斜,末端位置分析需要DEM提取测线高程,质量平衡计算需要DEM差分得到体积变化。可以这么说:DEM的精度直接影响下游所有派生产品的可靠性。所以GrIMP团队才会投入那么多精力,专门用高分辨率商业卫星影像来做一套全岛DEM,而不是直接拿全球粗分辨率DEM凑合。

1.3 为什么选 GeoEye 和 WorldView 立体像对

这个问题我在使用过程中想过很多次:为什么不直接用光学卫星图像加雷达干涉测高,或者干脆等ICESat-2的轨道数据密起来?原因其实很实际。

GeoEye-1的全色分辨率约0.41米,WorldView系列(WorldView-1/2/3)的全色分辨率在0.31到0.5米之间。这么高的空间分辨率可以捕捉冰裂隙、雪丘、岩石纹路等地表细节,给立体匹配提供充足的特征点。相比之下,Landsat的15米全色波段虽然也能做立体像对,但匹配精度完全不在一个量级上。

更重要的是,这些商业卫星具备同轨道立体采集能力——卫星在沿轨道飞行时,可以连续拍下同一地区的两幅影像,前后视角不同,形成天然立体像对。两幅影像获取间隔只有几十秒,地表在这几十秒内的变化可以忽略不计,这在动态的冰川地区是巨大的优势。如果是不同轨道、不同时间的立体像对,冰川本身已经流动了几十米,立体匹配就会彻底失败。

从时间维度看,WorldView-1于2007年发射,GeoEye-1于2008年发射,WorldView-3于2014年发射。这些卫星的影像覆盖从2007年延续到现在,恰好对应格陵兰冰盖变化的剧烈期。GrIMP DEM就把这个时间段的高精度影像压缩成了一套连续的高程产品,这是其他数据源很难做到的。

2. 数据生产链路:从两张卫星图到一张高程图

2.1 光学立体摄影测量的基本原理

理解这套DEM的质量,绕不开立体摄影测量的原理。简单说,当一个目标点被两个不同位置的传感器拍摄时,它在两幅影像中的位置会有差异,这个差异叫视差。视差越大,目标离传感器越近;视差越小,目标越远。人类大脑就是靠两只眼睛的视差感知深度,卫星立体像对也是同一个道理。

具体处理流程大概分四步:第一步,对两幅影像进行几何纠正,确定每幅影像拍摄时的外方位元素;第二步,在左右影像中寻找同名点,这个步骤叫密集影像匹配,全色0.4米分辨率的影像能产生密密麻麻的点云;第三步,利用前方交会原理把这些点的三维坐标计算出来,形成稀疏或密集的三维点云;第四步,把点云网格化,插值成规则格网的DEM。

这里有一个关键点:没有控制点的摄影测量,绝对高程精度会明显漂移。卫星定轨和姿态确定的误差,会直接传导到高程值上。几百公里的影像范围里,哪怕是微小的角误差,也会在末端造成几十米的高程偏差。所以GrIMP生产流程里,光靠立体像对远远不够,还需要用激光高程数据把整个网络“钉”住。

2.2 激光测高数据扮演的“地面控制点”角色

摄影测量里最经典的作业方式,是先布设地面控制点,用RTK或全站仪精确量测,再反求影像外方位元素。但在格陵兰冰盖上布地面控制点,既危险也不现实。GrIMP的替代方案非常聪明——用卫星和机载激光测高中的高精度点云充当虚拟控制点。

ICESat卫星装有GLAS激光测高仪,ICESat-2发射后装备了更强大的ATLAS激光雷达,机载的ATM(Airborne Topographic Mapper)系统也在格陵兰许多地区测量过。这些激光测高数据的特点是:沿轨方向精度极高(厘米级),但覆盖面不连续,像一根根针插在冰盖上。GrIMP团队把这些激光点作为高程控制源,用它们在立体像对网络做区域网平差,把影像解算出的三维点云整体对齐到一个精确的椭球高度上。

这个思路很像装修贴瓷砖时用的水平仪:激光尺提供一个绝对水平基准,工人拿着瓷砖一块块找平,只要每一块都以水平仪为参照,整体地坪就不会歪。没有这道工序,卫星影像本身的姿态误差就会让DEM“歪歪扭扭”,局部精度再高也白搭。这也是为什么我在实际使用中明显感觉,GrIMP DEM的整体高程基准比很多直接用原始影像生成的点云更可靠。

2.3 从2米原始DEM到30米产品,V001到V002做了哪些改进

这套数据最终对外发布的网格分辨率是30米,但它的原始处理并不仅限于30米。基于WorldView和GeoEye影像的立体匹配,原始DEM可以生成到2米甚至更高的分辨率。但2米全岛DEM的存储量和处理代价实在恐怖,而且很多区域因为没有好的控制点、云影干扰或低太阳角,2米网格反而会有大量噪声。所以GrIMP团队采用了分区域处理然后融合降尺度到30米的策略,兼顾了空间细节和稳定性。

V001版本其实存在几处让我在使用时很头疼的问题:一是冰盖陡峭边缘有比较明显的条带条纹,尤其是快速流动区的冰裂隙地带,DEM数值出现锯齿状波动;二是部分区域和激光测高数据相比存在系统性高程偏差;三是相邻影像拼接处偶尔有几十米的错位,看起来像“台阶”。

V002版本推出的核心改进,在我看来可以归纳为三点。第一,引入了更多新获取的影像,把时间覆盖向后延伸了好几年,直接可用于近十年的高程变化分析。第二,大规模使用了ICESat-2的ATL06高程产品作为控制源,相比V001时代的ICESat和ATM数据,控制点的密度和精度都有了明显提升,整体高程偏差被压下来了。第三,优化了区域网平差和镶嵌策略,之前那种接缝处“台阶”的问题大幅度减少。一句话:能用V002就尽量不要用V001

2.4 投影、基准与单位:最容易搞混的细节

极地数据处理里最阴险的坑往往不在算法,而在坐标参考系。这套DEM产品的常见投影是北极极射赤面投影,在EPSG代码里对应的是3413。你需要记住:EPSG:3413以WGS84椭球为基准,使用极射赤面投影方式,标准纬线设为北纬70°。所以在ArcticDEM、ITS_LIVE等极地产品里常看到“3413”这个数字,因为大家约定俗成统一在这个坐标系下工作。

更关键的是高程基准。GrIMP DEM的高程值是以WGS84椭球面为参考的椭球高,不是日常地图上的海拔高度。椭球高和海拔正高的差值,在格陵兰地区可以达到几十米。如果你把DEM高程和大地水准面模型、海面高度数据或者GPS大地高混用,必然出现系统性偏移。我自己见过不少新手把Ellipsoidal Height当成Orthometric Height去和海平面变化对比,结果得出了虚假的“冰面抬升”结论,回头还得重新处理。

单位方面,高程值单位是米,平面坐标单位也是米(在极射赤面投影下)。这个看似简单的问题,在实际处理中却经常被忽略——尤其是当DEM被转成经纬度坐标后,再被某些工具强行当“米”来处理时,计算出来的坡度和面积全都失真。

3. 产品档案与选型对比

3.1 数据档案速览

这里整理一份我常用的数据档案卡,方便查阅:

项目具体参数
数据集全称MEaSUREs Greenland Ice Mapping Project (GrIMP) Digital Elevation Model from GeoEye and WorldView Imagery
数据版本Version 2(V002)
数据中心NSIDC DAAC,数据集编号可检索NSIDC-0715
覆盖区域格陵兰岛冰盖及周边冰川,基本覆盖全岛冰面
时间跨度影像获取时间约2007年至2020年前后
网格分辨率30米
高程参考WGS84椭球高,单位为米
常见投影北极极射赤面投影(EPSG:3413),分块文件以GeoTIFF内嵌投影信息为准
原始影像GeoEye-1、WorldView-1/2/3全色波段立体像对
辅助控制数据ICESat、机载ATM、ICESat-2 ATL06激光测高数据

这套数据以分块(tile)形式发布,每个文件对应一个地理范围。从数据目录的命名可以看出区域范围,比如文件名里带经纬度或行列号信息。实际操作中我一般直接把整个目录下载下来,再用VRT(虚拟栅格)把它们拼成一张全岛图,这样后续处理效率最高。

3.2 V001 与 V002 版本选择建议

如果现在去NSIDC下载,默认拿到的就是V002。但有些同学早年间存过V001的旧文件,或者在某些第三方服务器上找到了旧版数据,这里我强烈建议不要为了省流量和省事继续用V001。

为什么?因为我在两个版本上都跑过同样的高程差分实验。在相对平坦的冰盖内部,两个版本差异不大,差值基本在5米以内。但在格陵兰东南部那些陡峭的出口冰川区域,V001和V002之间的差距经常达到10到20米,局部甚至更大。这个量级的误差会直接淹没掉真实的年际高程变化信号——格陵兰冰盖边缘区平均每年变化也就是几米以内。

V002对V001的关键修正包括:区域性高程漂移修正、更严格的影像控制网平差、新数据补充了旧版本的空洞区域。如果你要做时间序列分析,务必统一版本,绝对不要把V001和V002的tile混用在同一套镶嵌结果里,否则那些本以为是物理信号的异常差值,很可能只是版本差异在作怪。

3.3 和ArcticDEM、其他高程产品怎么选

极地领域常用的高分辨率DEM主要是GrIMP DEM和ArcticDEM,两者经常被放在一起比较。ArcticDEM由美国极地地理空间中心生产,同样是基于WorldView/GeoEye影像做立体摄影测量,分辨率达到2米,格陵兰大部分地区都有覆盖。

对比维度GrIMP DEM V002ArcticDEM
网格分辨率30米2米(也可降尺度使用)
基础影像GeoEye-1、WorldView系列WorldView系列为主
时间范围2007-2020前后,按时间版本区分多时相条带产品,可提取相对年代
镶嵌策略区域网平差后融合成稳定基准按时间拼接,局部地区存在云洞和条带
适用场景长时间尺度高程变化、区域物质平衡、冰流速地形改正局地精细地形、冰川形态分析、高分辨率地貌判读

我的选择经验是:核心任务如果是全岛尺度的高程变化或长时间序列分析,优先用GrIMP DEM V002,它的控制网和基准一致性更好;如果关注某条冰川末端的精细形态,比如冰崖高度、冰裂隙细节,那2米分辨率的ArcticDEM更有价值。两个数据并不互斥,经常是组合着用——ArcticDEM负责“看清细节”,GrIMP DEM负责“对齐时间”。

4. 实操全流程:从下载到出图的完整链路

4.1 从NSIDC下载数据:账号、检索、批量获取

格陵兰DEM这类NASA数据产品,托管在美国国家雪冰数据中心(NSIDC)的DAAC上,需要注册一个Earthdata账号才能下载。注册流程不复杂:访问Earthdata Login页面,填邮箱、设密码、勾选服务条款就行。认证方式现在推荐用NASA Earthdata的Token或者.netrc文件,直接在NSIDC的下载页面按提示操作即可。

数据检索可以直接在NSIDC的搜索界面里输入“MEaSUREs Greenland Ice Mapping Project DEM”或者数据集编号NSIDC-0715来定位。产品页会列出所有tile文件。我个人的习惯是先把整个产品目录用HTTrack或者wget脚本同步一遍,因为全岛DEM文件数量不少,一个个点下载太不现实。

这里给一个用wget批量下载的示例,前提是已经在本地配置好了.netrc认证:

wget -r -l1 -np -nH --cut-dirs=3 \ -A "*.tif" \ https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0715.002/

注意把URL替换成你实际看到的目录地址。下载完成后,先看一遍文件列表,确认tile的覆盖范围和命名规律,再进入拼接阶段。

4.2 用GDAL完成拼接、投影和浏览

拿到一堆GeoTIFF之后,第一步不是直接扣进ArcGIS或QGIS,而是先用命令行工具摸清数据的真实状态。GDAL是全套地理信息处理里最稳定的搭档。先检查一个文件的基本信息:

gdalinfo GrIMP_DEM_xx_yy.tif

重点关注这几项:投影信息、行列数、像素尺寸、NoData值、数据范围。如果投影不是EPSG:3413,后续统一投影时就需要指定目标坐标系。

然后做拼接。我通常不用一条gdal_merge到底,而是先构建VRT虚拟栅格,理由是可以先快速浏览所有tile的空间布局,而且VRT不复制像素数据,处理几百个文件几乎瞬间完成:

gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 grimp_dem_all.vrt *.tif

浏览拼接后的成果,可以直接在QGIS里把VRT拖进去看。需要导出为一张完整的大TIFF时,再执行gdal_translate,推荐使用COG(Cloud Optimized GeoTIFF)格式,后续切片发布效率很高:

gdal_translate -of COG grimp_dem_all.vrt grimp_dem_all_cog.tif

如果你需要把数据重投影到常见的经纬度坐标,用gdalwarp:

gdalwarp -t_srs EPSG:4326 -r cubic -dstnodata -9999 \ grimp_dem_all.vrt grimp_dem_wgs84.tif

重投影时建议用三次卷积插值(cubic)而不是最近邻,后者会在陡峭地形上产生明显的锯齿。

4.3 打开文件却一片黑?先检查NoData和拉伸设置

接触这套DEM的新手最容易遇到的问题,就是影像加载后整幅图黑乎乎的,看不出地形起伏。这通常不是数据坏了,而是显示拉伸的问题。DEM的高程值域加上NoData的极值,会让默认的灰度拉伸把有效信息压到极窄的动态范围里。

在QGIS中,打开图层属性,找到“符号系统”的“渲染类型”,选择“单波段假彩色”或“山区阴影”,再设置最小值为-100米左右、最大值为2000米左右,地形立体感就出来了。如果还看不出细节,可以配合“山体阴影”工具生成一个Hillshade图层叠加显示,视觉冲击力立刻不一样。

在Python里也可以用numpy快速检查数据范围:

from osgeo import gdal import numpy as np ds = gdal.Open('GrIMP_DEM_all.vrt') band = ds.GetRasterBand(1) arr = band.ReadAsArray() valid = arr[arr > -9990] # 过滤NoData print('像素统计 -> min:', valid.min(), 'max:', valid.max(), 'mean:', valid.mean(), 'std:', valid.std())

如果统计结果里出现正负几千甚至更大的值,多半是NoData没有被正确识别,需要回到gdalbuildvrt阶段检查srcnodata参数是否设置正确。

4.4 一个完整的多年高程差分析案例

找一条典型的格陵兰出口冰川,比如北部或东南部的冰川,演示一下高程差分析怎么做。假设我们需要2008年前后的DEM和2015年前后的DEM之间的高程变化。

当然,V002的tile可能混有多年的影像,严格的做法是使用逐tile的采集年份信息,把同一年份附近的数据放到一组。这里为了演示,先假设你拿到了三组tile:A组覆盖2008年前后,B组覆盖2015年前后。

第一步,分别对两组tile构建VRT并重投影到统一的EPSG:3413网格,用gdal_calc做差值:

gdal_calc.py \ -A dem_2008.vrt \ -B dem_2015.vrt \ --outfile=dh_2008_2015.tif \ --calc="B-A" \ --NoData=-9999

第二步,用numpy统计高程差在冰盖内部和边缘区的分布。建议先用一个格陵兰冰盖掩膜文件过滤掉基岩区,因为基岩区的高程变化应该接近零,可以用来评估数据误差底噪:

ds = gdal.Open('dh_2008_2015.tif') dh = ds.GetRasterBand(1).ReadAsArray() dh_valid = dh[(dh > -200) & (dh < 200)] print('中位数:', np.median(dh_valid)) print('标准差:', np.std(dh_valid))

如果整个冰盖内部的dh中位数明显偏离0,比如超过5米,说明两组镶嵌数据之间存在系统性高程偏移,这时候不要急着解释成“冰面抬升”,而要回到原始tile检查相对ICESat-2控制点的高程残差。我见过的靠谱研究,一般会先在稳定的基岩区域做验证,再谈论冰盖高程变化。

5. 典型应用场景:这套DEM到底能拿来干什么

5.1 冰面高程变化时间序列

利用不同时间段获取的DEM做逐区域差分,是重建格陵兰冰盖高程变化最直接的路径。GrIMP DEM的tile覆盖有时间标签,可以按时间段重组,生成2008年到2020年之间的多期高程差图。

这个应用的关键在于误差控制。冰盖内部的高程变化本身不大,一年只有半米到一米量级,如果两期DEM之间的配准误差超过两米,信号就全被噪声淹没了。所以实际操作中,我会先用稳定的基岩区和冰盖内部缓慢流动区做配准检验,如果发现系统性偏移,就对其中一期DEM做三维平移改正。GrIMP V002在这方面的底子比较好,因为它在生产时就已经用激光控制数据做了绝对基准统一,很多配准工作被前置了。

5.2 冰川末端测绘与崩解事件监测

冰川末端位置提取通常用Landsat或哨兵影像,但要量化冰崖高度、崩解后的体积损失,就必须有DEM支持。GrIMP DEM的高程数据可以与末端轮廓线叠加,快速提取冰前缘的海拔剖面,估算崩解事件的体积量级。

比如某条冰川在几个月内发生了大范围崩解,我们可以把崩解前后的DEM相减,得到体积变化总量,再除以崩解面积,得到平均厚度损失。这个流程现在几乎成了潮水冰川研究的标配。30米分辨率也许不足以看清单块崩解冰山的轮廓,但对区域尺度的体积核算来说已经非常够用。

5.3 冰流速反演中的地形改正

冰流速产品的生产通常采用影像特征追踪技术。卫星影像上的特征点在斜坡上移动时,除了冰川本身的运动,还有地形起伏造成的投影变形。如果没有DEM做透视改正,反演出的位移场在陡峭边缘区会出现系统性的“假速度”。

ITS_LIVE这类全球极地冰流速数据产品,在做特征追踪时就会用到DEM做地形扭曲消除。GrIMP DEM在这个环节里充当“基准地形”的角色:把每张卫星影像重投影到DEM对应的几何空间,消除地形起伏引起的像点位移,保留真实的水平位移信号。这也是为什么即使你不直接做高程分析,只要做极地影像相关工作时,手里也需要备一套高质量的DEM。

6. 踩坑实录:高频问题与排查方法

6.1 下载阶段的老大难:连接中断、认证失败、看得见下不动

NSIDC的服务器位于境外,国内用户下载大文件时经常遇到连接不稳定、断点续传支持不好的问题。我的解决办法是用支持断点续传的工具,比如lftp或者wget的-c参数,设置重试次数。如果单文件反复下载失败,按tile逐个下载反而比整体同步更可控。

还有一个容易忽略的细节:NSIDC的HTTPS下载要求客户端支持TLS相关版本,部分老旧wget版本会握手失败。升级到新版wget或者直接用curl --retry 5 --continue-at - 的方式,可以省掉很多麻烦。

6.2 拼接处出现“台阶”错位怎么办

即使V002已经改进了镶嵌质量,在实际拼接时仍有小概率遇到相邻tile之间的高程错位,特别是某些不同年份影像拼接的边界。这种错位的典型表现是:在坡度平缓区域看不出问题,但地形突然变陡的地方出现几米到十几米的“悬崖”伪像。

处理方法分两步。第一步,在QGIS里打开两个相邻tile,用“轮廓线”工具分别提取同一条山脊线或冰脊线的高程,直接看差值。第二步,如果确认错位,可以对其中一个tile做高程平移,用稳定的基岩区域计算中位数偏移量,然后对整个tile减去该偏移值。注意这种平移是全局的、常数偏移,不要对每个像元做匹配变换,否则会引入更复杂的变形。

6.3 冰面空域、NoData和插值

GrIMP DEM存在少量NoData区域,主要集中在极陡峭的峡湾两侧、阴影严重的长坡面以及部分云的残留区。这些空洞如果直接参与差分,会表现为异常大的正值或负值。我的习惯是:先建一个掩膜,把所有NoData区域扩大几圈缓冲带,后续所有分析都避开这些缓冲区域;如果某条测线必须穿过空洞,再用周围的DEM值做空间插值,同时严格记录插值范围,在报告里注明哪些区域是野外实测、哪些是插值估计。

6.4 新手最容易忽略的:高程是椭球高,不是海拔

这个坑我在前面提过一次,但值得再强调一遍,因为它太隐蔽了。GrIMP DEM的高程值本质上是在WGS84椭球面上量取的,不是相对于大地水准面的海拔高度。在格陵兰地区,椭球面和大地水准面之间的差距可能达到二三十米。

如果你把DEM高程和GPS大地高比较,使用的是椭球高,两者可以直接对比;但如果把DEM高程与海面高度或者潮位观测对比,就必须先做大地水准面改正。大家常用的格陵兰大地水准面模型可以从国际大地测量协会下载,转换后在对比。哪怕只是做DEM之间的差分,不同版本的大地水准面处理也可能带来虚假的区域内趋势,务必保持版本一致。

如果让我给刚接触这套数据的人一个建议,我会让他先下载两个相邻tile,叠加到一条典型出口冰川的测线上,看一眼断面形态。现成的山体阴影、冰裂隙纹理和末端崩解崖的立体感,比读十篇文档都直观。这套DEM真正体现了“数据产品”四个字的意义——把一批商业卫星原始影像和激光测高数据,加工成了可以直接上手做科研的标准化成果。以后遇到任何格陵兰高程相关的问题,记得先把它拿出来试试。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询