做冰川研究的人应该都有同感格陵兰冰盖到底在以多快的速度往海里送冰这个问题看起来简单但要给出一个可靠的数字比想象中麻烦得多。MEaSUREs 格陵兰岛年度冰盖速度镶嵌图V005就是目前公认上手门槛低、覆盖面全、可以直接拿来当底图用的数据产品之一。它把 SAR 和 Landsat 两类完全不同的遥感数据整合到统一的网格上给出了整个格陵兰冰盖逐年的二维冰流速度场。这篇文章我尽量把从数据原理、版本差异到下载读取、实际应用里会遇到的问题都讲一遍希望能帮你少走点弯路。这套数据适合谁只要是做冰冻圈、冰川动力学、海平面变化评估或者只是想快速了解格陵兰冰盖整体流动状态的人都值得花点时间把它的脾气摸清楚。它不是单纯的卫星影像而是一套经过配准、差分、镶嵌、误差校正的“研究级产品”所以里面有很多细节直接拿过来用很容易踩坑。1. 这个数据产品解决的是什么问题1.1 MEaSUREs 项目的定位把散装卫星数据变成“研究级”数据记录MEaSUREs 的全称是 Making Earth System Data Records for Use in Research Environments翻译过来就是“为研究环境制作地球系统数据记录”。这一看就知道和普通遥感数据分发不太一样——它的目标不是把原始卫星影像堆给你而是由专业团队把多年的多源卫星数据加工成一套连续、统一、可比较的数据记录。对普通科研人员来说自己从零开始处理 SAR 影像偏移追踪或者光学影像特征匹配单是数据量、配准精度和计算资源就能劝退大部分人了。MEaSUREs 的价值就是把这一层脏活累活打包好让你可以直接从栅格文件里读出某个格陵兰冰盖位置的速度分量。格陵兰冰盖速度产品只是 MEaSUREs 系列中的一个但它很有代表性。冰盖流速是衡量整个冰盖动态变化的最直接指标因为物质损失到底是靠表面融水径流还是靠冰流加速排冰两者对海平面上升的影响路径完全不同。速度场就是用来区分这两种机制的“分层器”。1.2 为什么格陵兰冰盖速度是“硬指标”格陵兰冰盖每年流失的数百亿吨冰里相当一部分是通过冰川快速流动把内陆的冰输送到海岸再以冰山崩解或入海融化结束行程。如果只观测表面高程变化你很难知道底部过程发生了什么。而冰流速数据可以直接告诉你哪条冰川在加速哪个冰流在减速这种变化往往比高程变化出现得更早、更明显。所以冰川动力学里有句行话速度场是冰盖的“心电图”。一条出海的冰流突然从每年 1 公里加速到每年 2 公里就意味着它的阻力机制发生了变化可能是底部融水润滑也可能是末端冰架崩解失去了支撑。没有速度场就只能看到结果看不到过程。1.3 数据产品面向的典型用户和场景这套数据最常用的场景包括冰川流速的长期趋势分析、冰盖数值模型初始化与验证、和 GRACE/ICESat 等其它遥感产品联合做物质平衡归因。哪怕是只做区域水文或者近岸海洋研究的人也经常拿它来界定冰川入海淡水通量的边界条件。我见过不少刚开始接触格陵兰数据的人第一反应是去下载 Sentinel-1 原始影像自己做偏移追踪结果被轨道校正、地形相位、去平、滤波这些步骤折磨得够呛。其实如果研究问题不要求精细到月尺度或者特定冰川的小尺度细节MEaSUREs 的年度镶嵌产品完全可以满足需求而且不同年份之间的一致性好很多。2. 数据源拆解SAR 和 Landsat 是怎么搭配工作的2.1 SAR 成像凭什么能穿透极夜和云层SAR 的全称是 Synthetic Aperture Radar也就是合成孔径雷达。这里必须先提醒一句你在搜索引擎里敲“SAR”经常会蹦出来股票行情软件里的同花顺副图指标那些跟遥感一点关系都没有。我们讨论的 SAR 是主动式微波遥感雷达自己发射微波脉冲并接收地面回波。为什么冰盖测速特别依赖 SAR首先极地地区经常被云覆盖而且极夜期间几乎没有太阳光光学影像直接“瞎掉”。SAR 工作在厘米到分米级波长的微波波段可以穿透云层而且不依赖光照无论白天黑夜都能成像。其次SAR 对地表纹理和后向散射特性非常敏感冰裂隙、雪丘、粒雪层的变化都能在影像上留下可追踪的标记。第三SAR 不仅有强度信息还有相位信息相位差可以用于干涉测量InSAR精度可以达到厘米级甚至毫米级。不过在快速流动的冰川上InSAR 有个致命弱点相位差太大容易“混叠”。流速超过一定程度后干涉条纹变得极其密集解缠困难。所以针对格陵兰这种存在大量急流ice stream的地区数据处理时往往更依赖偏移追踪offset tracking也就是在两期 SAR 影像的强度图上做局部互相关找出地表特征在时间间隔内的位移量。偏移追踪的精度虽然不如 InSAR但胜在鲁棒不会因为速度太快而彻底失效。2.2 Landsat 光学影像在传统 SAR 方案里的角色很多人疑惑既然 SAR 已经可以测冰速了为什么还要掺和 Landsat 光学影像答案很简单空间覆盖和时间配对不够。SAR 卫星的存档数据虽然长但不同轨道、不同入射角的覆盖范围有限单靠 SAR 很难在某个时间窗口内把整个格陵兰冰盖完整拼起来。Landsat 系列特别是 Landsat 8/9的 OLI 传感器提供了 30 米分辨率的多光谱影像重访周期 16 天理论上在一两个季节内可以获取覆盖整个冰盖的多期光学影像。光学影像测速的基本思路叫特征追踪feature tracking找两期影像里冰面共同存在的纹理特征比如裂缝、雪丘、裸露冰面的亮度差异然后计算这些特征在时间间隔里移动了多少像素。根据像素位移和时间间隔就能换算出流速。这里的“特征”其实就是自然冰面自身的纹理不需要人工布设标志物所以是一种完全被动式的遥感测量。不过光学影像对天气要求苛刻只要有云就会挡住地面特征。格陵兰沿海地区夏季多云真正无云且质量好的影像并没那么容易拿到。所以实际产品里Landsat 和 SAR 是互相补位的晴天多的时候用光学不差这一两天的话就用 SAR 补齐缺口。2.3 两种数据源的融合策略与取舍V005 产品的整合逻辑其实并不复杂先把每一对可用的影像对分别解算出各自的位移场然后做时间归一化统一折算成每年的位移量再进行坐标转换和网格重采样最后用加权平均或者中值合成的方式拼成一张覆盖全岛的年度速度图。不同数据源的精度并不一样。SAR 偏移追踪在冰面纹理清晰、时间间隔合适时精度较高Landsat 特征追踪在平坦的粒雪区可能会因为对比度不足而失效。所以在融合阶段处理团队会给不同来源的测量赋予不同的权重还会做异常值剔除。这个“加权”的思路很关键它不是简单把多张图叠加取平均而是对每一条测量记录的质量做出估计然后根据不确定性来分配影响。从我自己的使用体验来看V005 在格陵兰中部的粒雪区冰面比较均匀、缺乏纹理偶尔会出现一些条纹状的噪声伪影这些区域速度本身很慢通常每年只有几米到几十米如果做全区统计时不去掩膜很容易被当成异常加速区。使用前了解这一点很重要。3. V005 版本的关键更新与数据细节3.1 从早期版本到 V005变化在哪里MEaSUREs 格陵兰冰速产品已经迭代了好几个版本。V005 相比之前版本一个明显的改进是整合了更新的 Sentinel-1 和 Landsat 8/9 数据同时引入了更长时间序列的存档数据参与合成。数据越多时间窗口被填补得越完整最终镶嵌图的空洞也就越少。另一个可能是处理流水线的升级更精细的轨道校正、更可靠的误差传播、改进后的加权合成算法。这些内部变化普通人不容易直接感知但体现在产品上就是速度异常值更少、海岸线和冰舌边缘的数值更平滑。V005 里有些冰川前端的异常向量明显减少这是处理链路进步的直接迹象。顺便说一句网上很多教程基于的是旧的 V004 或者更早的版本文件名后缀和变量名可能对不上。如果你照着老教程跑代码发现读不出变量先别急着怀疑自己先确认下载的版本是不是 V005。3.2 投影、分辨率、变量、单位读数据前必须知道的几件事这套数据用的是极地立体投影Polar Stereographic在 NSIDC 的 EASE-Grid 2.0 体系里对应的是 EPSG:3413。中心点大致在格陵兰岛中部纬度 70 度标准纬线。它和常见的地理坐标系经纬度 EPSG:4326完全是两回事。如果你直接把数据当经纬度画图结果会扭曲得非常离谱必须做重投影。分辨率方面标准产品通常输出 500 米网格部分子集或区域产品可能有更高分辨率。这里注意500 米是“镶嵌网格”的分辨率不等于原始遥感的最高分辨率。Landsat 30 米、Sentinel-1 也有十几米的分辨率镶嵌过程相当于把高分辨率信息聚合到了统一的粗网格上。好处是整个岛的数据量可控坏处是小冰川、窄冰舌的空间细节会丢失。变量方面最常用的包括速度在东西方向的分量vx、南北方向的分量vy、速度大小magnitude 或 speed以及对应的误差字段。单位需要特别留意这套产品的速度单位是米每年m/yr不是米每天也不是公里每年。我见过有人把 1000 m/yr 的数据直接当成 1000 m/day 用结果算出来的冰川通量大了 365 倍这个错误在论文里是非常尴尬的。3.3 年度镶嵌图是怎么“镶嵌”出来的“年度”这个词容易让人误解成“整年只有一张图”其实更准确的说法是一个年度的合成结果mosaic是由那个时段内获取的多个影像对共同拼接得到的。由于冰盖边缘的冰川流速快且受季节性融水影响明显一套年度数据实际上混合了不同季节的观测。处理团队通常会优先选择时间间隔在 10 到 100 天左右的影像对间隔太短位移量太小无法从噪声里提取信号间隔太长则可能跨越季节性变化导致速度被平均得不自然。经过筛选的影像对解算出位移场均摊到年尺度后再按质量加权合成。所以你在 V005 里看到的“年度速度”本质上是一个经过时间归一化后的合成的“典型年速度”不一定等于某个具体自然年 365 天的精确平均值。做跨年份趋势分析时这一点必须牢记如果相邻两个年份的合成数据里混入了不同季节的观测比例那么一两年的小幅度变化可能并不代表真实的物理加速或减速至少需要连续 5 年以上的趋势才能下结论。4. 获取数据与本地处理的完整流程4.1 注册 Earthdata 账号并下载数据这套数据由 NSIDC DAAC 分发下载前需要先注册 NASA Earthdata 账号。注册免费过程不复杂但要注意登录 NSIDC 下载页时即使你是通过 NSIDC 的网页跳转过去的也必须使用同一个 Earthdata 账号完成认证。下载方式有几种。最简单的是直接在 NSIDC 的数据检索页面上手动选择时间范围、区域、版本然后下单排队下载。如果你需要下载多年的全部数据手动点下载就很痛苦了我推荐用 Python 的 earthaccess 库几行代码就能批量查找和下载。import earthaccess # 登录 Earthdata auth earthaccess.login(strategynetrc) # 搜索 MEaSUREs 格陵兰冰速数据 results earthaccess.search_data( short_nameNSIDC-0725, # 具体的短名称以 NSIDC 检索结果为准 version5, bounding_box(-73, 59, -12, 84) ) # 批量下载到本地 files earthaccess.download(results, ./data/)注意 short_name 需要以 NSIDC 当前检索页里的实际数据集短名为准版本写 5。如果你不确定先在页面上搜 Greenland Ice Sheet Velocity看到对应条目后点进去看详情页里的 Data Set ID。4.2 用 Python 读取 GeoTIFF 并做基础可视化下载下来的文件通常是 GeoTIFF 格式读取最稳妥的工具是 GDAL或者直接用 rioxarray 和 matplotlib 做可视化。我最常用的是一套组合拳rioxarray 读数据cartopy 处理投影和底图。import rioxarray import matplotlib.pyplot as plt import cartopy.crs as ccrs # 读取速度幅值 ds rioxarray.open_rasterio(Greenland_ice_speed_v005.tif) speed ds.squeeze() # 去掉冗余维度 # 注意这个产品的坐标已经是 EPSG:3413 # 如果需要转成经纬度显示先做重投影 speed_ll speed.rio.reproject(EPSG:4326) fig, ax plt.subplots( figsize(10, 10), subplot_kw{projection: ccrs.PlateCarree()} ) im ax.imshow( speed_ll.values, extent[speed_ll.x.min(), speed_ll.x.max(), speed_ll.y.min(), speed_ll.y.max()], transformccrs.PlateCarree(), cmapviridis, vmin0, vmax3000 ) ax.coastlines() plt.colorbar(im, labelIce speed (m/yr)) plt.show()这里有个小坑GeoTIFF 里可能带有 nodata 值读取后如果不对 nodata 做掩膜速度数值里会出现一个极大或极小的填充值直接影响色标范围和统计结果。用 rioxarray 时最好检查一下ds.rio.nodata属性必要时手动设置speed speed.where(speed ! ds.rio.nodata)另一个容易忽略的是数据文件的坐标系是否被正确识别。如果文件本身没有内嵌坐标系信息GDAL 会把它当成未知坐标rio.reproject就会报错。遇到这种文件需要手工指定坐标系通常就是 EPSG:3413。speed speed.rio.write_crs(EPSG:3413, inplaceTrue)4.3 常见的坐标与单位坑坐标和单位是我见过最多的两个翻车点。先讲坐标如果你用 cartopy 直接叠加 Natural Earth 的海岸线一定要确保海岸线的坐标系和冰速数据的坐标系一致。最省事的方法是先把冰速数据重投影成经纬度再画地图这样底图切割和范围设置都简单。单位的问题前面提过m/yr 和 m/day 差 365 倍计算冰通量时如果从这篇论文拿到的是 m/yr 的栅格从那篇文献里查到的是 m/day 的经验公式不做换算就会出大事。我的习惯是拿到数据后第一件事就是查变量单位并且把单位写进脚本的注释里避免过了一个月再回来看代码时猜谜。还有一个网格对齐的细节不同来源的栅格数据即使分辨率标称都是 500 米网格起点也未必一致。如果你要把冰速数据和其它栅格比如冰厚度、表面高程做逐像元计算一定要先检查网格是否完全对齐。最省心的方式是统一用 xarray 的reproject_match函数把多个栅格重采样到同一个网格上。5. 冰速数据能用在哪些地方5.1 典型应用冰川流速变化与物质平衡估算把逐年速度图叠起来看最直观的应用就是识别哪条冰川在加速。格陵兰最典型的例子是 Jakobshavn Isbræ雅各布港冰川它在 2000 年代经历了一波剧烈的加速和崩解后期又出现过减速。这种变化用速度场的多年序列一眼就能看出来。另一个常见应用是计算跨断面冰通量。你在冰川的某个断面划一条线把线上每个像素的冰厚和速度垂直分量相乘再积分就可以估算这条断面一年向海洋输送了多少体积的冰。配合表面物质平衡积雪积累减消融的数据就能知道某个流域到底是收入大于支出还是支出大于收入。这个计算对网格对齐的要求极高速度场和冰厚数据的网格差一个像素断面积分结果可能差出百分之几。5.2 与其它遥感产品的联合使用思路冰速和表面高程变化来自 ICESat/ICESat-2、ATM 或 ArcticDEM配合使用价值很大。如果某个区域冰速在加快同时表面高程在降低基本可以判断是“动力排冰”在起作用如果冰速没太大变化但高程持续下降那大概率是表面融化主导。这种归因分析在 IPCC 式的评估报告里非常常见。冰速场还可以直接用于冰盖模型的初始化。模型里需要给定初始速度场来反推基底滑动系数如果初始速度本身噪声太大反演出的基底摩擦系数就会出现很多不真实的条纹。MEaSUREs 这类平滑过后的产品反而比原始的高分辨率测量更适合做模型输入。对做海洋研究的人来说冰速变化决定了冰川前缘淡水通量和冰山崩解通量的时空分布。格陵兰周边峡湾的环流、生态系统的营养物质输入都跟这些冰通量变化有关。这种跨领域的需求其实很常见所以哪怕你不是冰冻圈专业出身了解这套数据的结构也有实际意义。5.3 处理技巧从区域裁剪到异常值剔除实际操作里我一般不建议直接对整岛范围做统计边缘区域和峡湾里的噪声会影响结论。如果研究区域集中在某个流域先裁剪流域边界再处理会省很多心。NSIDC 有提供格陵兰流域边界数据或者你也可以用自然地球的水文数据做裁剪。对异常速度的识别我通常用一个简单的经验法则速度方向与局部冰流方向相差超过 60 度且速度又比较低的像元大概率是噪声速度超过 5 km/yr 的像元也可能有问题除非你知道那条冰川真的在快速流动。滤波时可以用中值滤波或高斯滤波降低孤立异常的影响但滤波窗口不要太大否则会把真实的窄冰流信号也抹平。6. 实操中踩过的坑与排查技巧6.1 常见问题速查表问题现象可能原因解决方法速度图整体偏移坐标系没识别或搞错检查 EPSG:3413不要用经纬度直接画图速度数值大得离谱单位混淆 m/yr 和 m/day确认变量说明换算单位后再使用图上有条纹状噪声粒雪区纹理不足导致匹配失败做中值滤波或掩膜低纹理区域某个年份整幅图空洞很多该年度云覆盖严重、SAR 数据不足改用前后年份差值或者使用更高版本与其它栅格叠加后边缘不齐网格原点对不齐用reproject_match统一网格下载需要登录却始终失败Earthdata 账号未与 NSIDC 关联在 NSIDC 页面上使用同一个账号登录6.2 几种容易误判的异常信号有一种情况容易误导新手在冰川的陡峭边缘处速度图会出现局部加速或反向运动的假象。这是因为 SAR 偏移追踪和光学特征追踪在陡峭地形上都有几何误差尤其是雷达侧视造成的阴影和叠掩效应。遇到孤立的高速度像元先别急着激动看一眼它是不是紧贴着高差大的山地边缘如果是大概率是地形伪影。还有一种情况是冰舌前端的季节性崩解造成的“假信号”。某个像元上一年还在冰舌上这一年冰舌崩解后它变成了开阔水面特征追踪会把这种与冰面纹理无关的变化也当成位移计算出来形成看起来像高速运动的错误向量。年度产品对这种瞬时事件的处理能力有限如果你关注的就是冰舌前端变化建议用更短时间尺度的原始影像数据复核。最后提醒一点下载数据后最好记录一下数据的生产时间和版本说明。NSIDC 会不定期更新数据产品版本同一个 short_name 的 V004 和 V005 之间可能有很大的算法差异不同版本的数据不建议直接混用。如果你在某篇论文里看到别人引用了 V005那就尽量用同一版本复现避免因版本差异产生不必要的争议。我自己的体会是MEaSUREs 格陵兰冰速镶嵌图这种“看起来简单、用起来细碎”的遥感产品最考验人的不是下载和画图而是对数据背后观测原理和误差来源的理解。把 SAR 为什么能干这个活、Landsat 为什么来帮忙、V005 到底改了什么这些问题想清楚遇到奇怪的结果时才不会慌。
MEaSUREs格陵兰冰速数据V005：原理、读取与应用全解析
做冰川研究的人应该都有同感格陵兰冰盖到底在以多快的速度往海里送冰这个问题看起来简单但要给出一个可靠的数字比想象中麻烦得多。MEaSUREs 格陵兰岛年度冰盖速度镶嵌图V005就是目前公认上手门槛低、覆盖面全、可以直接拿来当底图用的数据产品之一。它把 SAR 和 Landsat 两类完全不同的遥感数据整合到统一的网格上给出了整个格陵兰冰盖逐年的二维冰流速度场。这篇文章我尽量把从数据原理、版本差异到下载读取、实际应用里会遇到的问题都讲一遍希望能帮你少走点弯路。这套数据适合谁只要是做冰冻圈、冰川动力学、海平面变化评估或者只是想快速了解格陵兰冰盖整体流动状态的人都值得花点时间把它的脾气摸清楚。它不是单纯的卫星影像而是一套经过配准、差分、镶嵌、误差校正的“研究级产品”所以里面有很多细节直接拿过来用很容易踩坑。1. 这个数据产品解决的是什么问题1.1 MEaSUREs 项目的定位把散装卫星数据变成“研究级”数据记录MEaSUREs 的全称是 Making Earth System Data Records for Use in Research Environments翻译过来就是“为研究环境制作地球系统数据记录”。这一看就知道和普通遥感数据分发不太一样——它的目标不是把原始卫星影像堆给你而是由专业团队把多年的多源卫星数据加工成一套连续、统一、可比较的数据记录。对普通科研人员来说自己从零开始处理 SAR 影像偏移追踪或者光学影像特征匹配单是数据量、配准精度和计算资源就能劝退大部分人了。MEaSUREs 的价值就是把这一层脏活累活打包好让你可以直接从栅格文件里读出某个格陵兰冰盖位置的速度分量。格陵兰冰盖速度产品只是 MEaSUREs 系列中的一个但它很有代表性。冰盖流速是衡量整个冰盖动态变化的最直接指标因为物质损失到底是靠表面融水径流还是靠冰流加速排冰两者对海平面上升的影响路径完全不同。速度场就是用来区分这两种机制的“分层器”。1.2 为什么格陵兰冰盖速度是“硬指标”格陵兰冰盖每年流失的数百亿吨冰里相当一部分是通过冰川快速流动把内陆的冰输送到海岸再以冰山崩解或入海融化结束行程。如果只观测表面高程变化你很难知道底部过程发生了什么。而冰流速数据可以直接告诉你哪条冰川在加速哪个冰流在减速这种变化往往比高程变化出现得更早、更明显。所以冰川动力学里有句行话速度场是冰盖的“心电图”。一条出海的冰流突然从每年 1 公里加速到每年 2 公里就意味着它的阻力机制发生了变化可能是底部融水润滑也可能是末端冰架崩解失去了支撑。没有速度场就只能看到结果看不到过程。1.3 数据产品面向的典型用户和场景这套数据最常用的场景包括冰川流速的长期趋势分析、冰盖数值模型初始化与验证、和 GRACE/ICESat 等其它遥感产品联合做物质平衡归因。哪怕是只做区域水文或者近岸海洋研究的人也经常拿它来界定冰川入海淡水通量的边界条件。我见过不少刚开始接触格陵兰数据的人第一反应是去下载 Sentinel-1 原始影像自己做偏移追踪结果被轨道校正、地形相位、去平、滤波这些步骤折磨得够呛。其实如果研究问题不要求精细到月尺度或者特定冰川的小尺度细节MEaSUREs 的年度镶嵌产品完全可以满足需求而且不同年份之间的一致性好很多。2. 数据源拆解SAR 和 Landsat 是怎么搭配工作的2.1 SAR 成像凭什么能穿透极夜和云层SAR 的全称是 Synthetic Aperture Radar也就是合成孔径雷达。这里必须先提醒一句你在搜索引擎里敲“SAR”经常会蹦出来股票行情软件里的同花顺副图指标那些跟遥感一点关系都没有。我们讨论的 SAR 是主动式微波遥感雷达自己发射微波脉冲并接收地面回波。为什么冰盖测速特别依赖 SAR首先极地地区经常被云覆盖而且极夜期间几乎没有太阳光光学影像直接“瞎掉”。SAR 工作在厘米到分米级波长的微波波段可以穿透云层而且不依赖光照无论白天黑夜都能成像。其次SAR 对地表纹理和后向散射特性非常敏感冰裂隙、雪丘、粒雪层的变化都能在影像上留下可追踪的标记。第三SAR 不仅有强度信息还有相位信息相位差可以用于干涉测量InSAR精度可以达到厘米级甚至毫米级。不过在快速流动的冰川上InSAR 有个致命弱点相位差太大容易“混叠”。流速超过一定程度后干涉条纹变得极其密集解缠困难。所以针对格陵兰这种存在大量急流ice stream的地区数据处理时往往更依赖偏移追踪offset tracking也就是在两期 SAR 影像的强度图上做局部互相关找出地表特征在时间间隔内的位移量。偏移追踪的精度虽然不如 InSAR但胜在鲁棒不会因为速度太快而彻底失效。2.2 Landsat 光学影像在传统 SAR 方案里的角色很多人疑惑既然 SAR 已经可以测冰速了为什么还要掺和 Landsat 光学影像答案很简单空间覆盖和时间配对不够。SAR 卫星的存档数据虽然长但不同轨道、不同入射角的覆盖范围有限单靠 SAR 很难在某个时间窗口内把整个格陵兰冰盖完整拼起来。Landsat 系列特别是 Landsat 8/9的 OLI 传感器提供了 30 米分辨率的多光谱影像重访周期 16 天理论上在一两个季节内可以获取覆盖整个冰盖的多期光学影像。光学影像测速的基本思路叫特征追踪feature tracking找两期影像里冰面共同存在的纹理特征比如裂缝、雪丘、裸露冰面的亮度差异然后计算这些特征在时间间隔里移动了多少像素。根据像素位移和时间间隔就能换算出流速。这里的“特征”其实就是自然冰面自身的纹理不需要人工布设标志物所以是一种完全被动式的遥感测量。不过光学影像对天气要求苛刻只要有云就会挡住地面特征。格陵兰沿海地区夏季多云真正无云且质量好的影像并没那么容易拿到。所以实际产品里Landsat 和 SAR 是互相补位的晴天多的时候用光学不差这一两天的话就用 SAR 补齐缺口。2.3 两种数据源的融合策略与取舍V005 产品的整合逻辑其实并不复杂先把每一对可用的影像对分别解算出各自的位移场然后做时间归一化统一折算成每年的位移量再进行坐标转换和网格重采样最后用加权平均或者中值合成的方式拼成一张覆盖全岛的年度速度图。不同数据源的精度并不一样。SAR 偏移追踪在冰面纹理清晰、时间间隔合适时精度较高Landsat 特征追踪在平坦的粒雪区可能会因为对比度不足而失效。所以在融合阶段处理团队会给不同来源的测量赋予不同的权重还会做异常值剔除。这个“加权”的思路很关键它不是简单把多张图叠加取平均而是对每一条测量记录的质量做出估计然后根据不确定性来分配影响。从我自己的使用体验来看V005 在格陵兰中部的粒雪区冰面比较均匀、缺乏纹理偶尔会出现一些条纹状的噪声伪影这些区域速度本身很慢通常每年只有几米到几十米如果做全区统计时不去掩膜很容易被当成异常加速区。使用前了解这一点很重要。3. V005 版本的关键更新与数据细节3.1 从早期版本到 V005变化在哪里MEaSUREs 格陵兰冰速产品已经迭代了好几个版本。V005 相比之前版本一个明显的改进是整合了更新的 Sentinel-1 和 Landsat 8/9 数据同时引入了更长时间序列的存档数据参与合成。数据越多时间窗口被填补得越完整最终镶嵌图的空洞也就越少。另一个可能是处理流水线的升级更精细的轨道校正、更可靠的误差传播、改进后的加权合成算法。这些内部变化普通人不容易直接感知但体现在产品上就是速度异常值更少、海岸线和冰舌边缘的数值更平滑。V005 里有些冰川前端的异常向量明显减少这是处理链路进步的直接迹象。顺便说一句网上很多教程基于的是旧的 V004 或者更早的版本文件名后缀和变量名可能对不上。如果你照着老教程跑代码发现读不出变量先别急着怀疑自己先确认下载的版本是不是 V005。3.2 投影、分辨率、变量、单位读数据前必须知道的几件事这套数据用的是极地立体投影Polar Stereographic在 NSIDC 的 EASE-Grid 2.0 体系里对应的是 EPSG:3413。中心点大致在格陵兰岛中部纬度 70 度标准纬线。它和常见的地理坐标系经纬度 EPSG:4326完全是两回事。如果你直接把数据当经纬度画图结果会扭曲得非常离谱必须做重投影。分辨率方面标准产品通常输出 500 米网格部分子集或区域产品可能有更高分辨率。这里注意500 米是“镶嵌网格”的分辨率不等于原始遥感的最高分辨率。Landsat 30 米、Sentinel-1 也有十几米的分辨率镶嵌过程相当于把高分辨率信息聚合到了统一的粗网格上。好处是整个岛的数据量可控坏处是小冰川、窄冰舌的空间细节会丢失。变量方面最常用的包括速度在东西方向的分量vx、南北方向的分量vy、速度大小magnitude 或 speed以及对应的误差字段。单位需要特别留意这套产品的速度单位是米每年m/yr不是米每天也不是公里每年。我见过有人把 1000 m/yr 的数据直接当成 1000 m/day 用结果算出来的冰川通量大了 365 倍这个错误在论文里是非常尴尬的。3.3 年度镶嵌图是怎么“镶嵌”出来的“年度”这个词容易让人误解成“整年只有一张图”其实更准确的说法是一个年度的合成结果mosaic是由那个时段内获取的多个影像对共同拼接得到的。由于冰盖边缘的冰川流速快且受季节性融水影响明显一套年度数据实际上混合了不同季节的观测。处理团队通常会优先选择时间间隔在 10 到 100 天左右的影像对间隔太短位移量太小无法从噪声里提取信号间隔太长则可能跨越季节性变化导致速度被平均得不自然。经过筛选的影像对解算出位移场均摊到年尺度后再按质量加权合成。所以你在 V005 里看到的“年度速度”本质上是一个经过时间归一化后的合成的“典型年速度”不一定等于某个具体自然年 365 天的精确平均值。做跨年份趋势分析时这一点必须牢记如果相邻两个年份的合成数据里混入了不同季节的观测比例那么一两年的小幅度变化可能并不代表真实的物理加速或减速至少需要连续 5 年以上的趋势才能下结论。4. 获取数据与本地处理的完整流程4.1 注册 Earthdata 账号并下载数据这套数据由 NSIDC DAAC 分发下载前需要先注册 NASA Earthdata 账号。注册免费过程不复杂但要注意登录 NSIDC 下载页时即使你是通过 NSIDC 的网页跳转过去的也必须使用同一个 Earthdata 账号完成认证。下载方式有几种。最简单的是直接在 NSIDC 的数据检索页面上手动选择时间范围、区域、版本然后下单排队下载。如果你需要下载多年的全部数据手动点下载就很痛苦了我推荐用 Python 的 earthaccess 库几行代码就能批量查找和下载。import earthaccess # 登录 Earthdata auth earthaccess.login(strategynetrc) # 搜索 MEaSUREs 格陵兰冰速数据 results earthaccess.search_data( short_nameNSIDC-0725, # 具体的短名称以 NSIDC 检索结果为准 version5, bounding_box(-73, 59, -12, 84) ) # 批量下载到本地 files earthaccess.download(results, ./data/)注意 short_name 需要以 NSIDC 当前检索页里的实际数据集短名为准版本写 5。如果你不确定先在页面上搜 Greenland Ice Sheet Velocity看到对应条目后点进去看详情页里的 Data Set ID。4.2 用 Python 读取 GeoTIFF 并做基础可视化下载下来的文件通常是 GeoTIFF 格式读取最稳妥的工具是 GDAL或者直接用 rioxarray 和 matplotlib 做可视化。我最常用的是一套组合拳rioxarray 读数据cartopy 处理投影和底图。import rioxarray import matplotlib.pyplot as plt import cartopy.crs as ccrs # 读取速度幅值 ds rioxarray.open_rasterio(Greenland_ice_speed_v005.tif) speed ds.squeeze() # 去掉冗余维度 # 注意这个产品的坐标已经是 EPSG:3413 # 如果需要转成经纬度显示先做重投影 speed_ll speed.rio.reproject(EPSG:4326) fig, ax plt.subplots( figsize(10, 10), subplot_kw{projection: ccrs.PlateCarree()} ) im ax.imshow( speed_ll.values, extent[speed_ll.x.min(), speed_ll.x.max(), speed_ll.y.min(), speed_ll.y.max()], transformccrs.PlateCarree(), cmapviridis, vmin0, vmax3000 ) ax.coastlines() plt.colorbar(im, labelIce speed (m/yr)) plt.show()这里有个小坑GeoTIFF 里可能带有 nodata 值读取后如果不对 nodata 做掩膜速度数值里会出现一个极大或极小的填充值直接影响色标范围和统计结果。用 rioxarray 时最好检查一下ds.rio.nodata属性必要时手动设置speed speed.where(speed ! ds.rio.nodata)另一个容易忽略的是数据文件的坐标系是否被正确识别。如果文件本身没有内嵌坐标系信息GDAL 会把它当成未知坐标rio.reproject就会报错。遇到这种文件需要手工指定坐标系通常就是 EPSG:3413。speed speed.rio.write_crs(EPSG:3413, inplaceTrue)4.3 常见的坐标与单位坑坐标和单位是我见过最多的两个翻车点。先讲坐标如果你用 cartopy 直接叠加 Natural Earth 的海岸线一定要确保海岸线的坐标系和冰速数据的坐标系一致。最省事的方法是先把冰速数据重投影成经纬度再画地图这样底图切割和范围设置都简单。单位的问题前面提过m/yr 和 m/day 差 365 倍计算冰通量时如果从这篇论文拿到的是 m/yr 的栅格从那篇文献里查到的是 m/day 的经验公式不做换算就会出大事。我的习惯是拿到数据后第一件事就是查变量单位并且把单位写进脚本的注释里避免过了一个月再回来看代码时猜谜。还有一个网格对齐的细节不同来源的栅格数据即使分辨率标称都是 500 米网格起点也未必一致。如果你要把冰速数据和其它栅格比如冰厚度、表面高程做逐像元计算一定要先检查网格是否完全对齐。最省心的方式是统一用 xarray 的reproject_match函数把多个栅格重采样到同一个网格上。5. 冰速数据能用在哪些地方5.1 典型应用冰川流速变化与物质平衡估算把逐年速度图叠起来看最直观的应用就是识别哪条冰川在加速。格陵兰最典型的例子是 Jakobshavn Isbræ雅各布港冰川它在 2000 年代经历了一波剧烈的加速和崩解后期又出现过减速。这种变化用速度场的多年序列一眼就能看出来。另一个常见应用是计算跨断面冰通量。你在冰川的某个断面划一条线把线上每个像素的冰厚和速度垂直分量相乘再积分就可以估算这条断面一年向海洋输送了多少体积的冰。配合表面物质平衡积雪积累减消融的数据就能知道某个流域到底是收入大于支出还是支出大于收入。这个计算对网格对齐的要求极高速度场和冰厚数据的网格差一个像素断面积分结果可能差出百分之几。5.2 与其它遥感产品的联合使用思路冰速和表面高程变化来自 ICESat/ICESat-2、ATM 或 ArcticDEM配合使用价值很大。如果某个区域冰速在加快同时表面高程在降低基本可以判断是“动力排冰”在起作用如果冰速没太大变化但高程持续下降那大概率是表面融化主导。这种归因分析在 IPCC 式的评估报告里非常常见。冰速场还可以直接用于冰盖模型的初始化。模型里需要给定初始速度场来反推基底滑动系数如果初始速度本身噪声太大反演出的基底摩擦系数就会出现很多不真实的条纹。MEaSUREs 这类平滑过后的产品反而比原始的高分辨率测量更适合做模型输入。对做海洋研究的人来说冰速变化决定了冰川前缘淡水通量和冰山崩解通量的时空分布。格陵兰周边峡湾的环流、生态系统的营养物质输入都跟这些冰通量变化有关。这种跨领域的需求其实很常见所以哪怕你不是冰冻圈专业出身了解这套数据的结构也有实际意义。5.3 处理技巧从区域裁剪到异常值剔除实际操作里我一般不建议直接对整岛范围做统计边缘区域和峡湾里的噪声会影响结论。如果研究区域集中在某个流域先裁剪流域边界再处理会省很多心。NSIDC 有提供格陵兰流域边界数据或者你也可以用自然地球的水文数据做裁剪。对异常速度的识别我通常用一个简单的经验法则速度方向与局部冰流方向相差超过 60 度且速度又比较低的像元大概率是噪声速度超过 5 km/yr 的像元也可能有问题除非你知道那条冰川真的在快速流动。滤波时可以用中值滤波或高斯滤波降低孤立异常的影响但滤波窗口不要太大否则会把真实的窄冰流信号也抹平。6. 实操中踩过的坑与排查技巧6.1 常见问题速查表问题现象可能原因解决方法速度图整体偏移坐标系没识别或搞错检查 EPSG:3413不要用经纬度直接画图速度数值大得离谱单位混淆 m/yr 和 m/day确认变量说明换算单位后再使用图上有条纹状噪声粒雪区纹理不足导致匹配失败做中值滤波或掩膜低纹理区域某个年份整幅图空洞很多该年度云覆盖严重、SAR 数据不足改用前后年份差值或者使用更高版本与其它栅格叠加后边缘不齐网格原点对不齐用reproject_match统一网格下载需要登录却始终失败Earthdata 账号未与 NSIDC 关联在 NSIDC 页面上使用同一个账号登录6.2 几种容易误判的异常信号有一种情况容易误导新手在冰川的陡峭边缘处速度图会出现局部加速或反向运动的假象。这是因为 SAR 偏移追踪和光学特征追踪在陡峭地形上都有几何误差尤其是雷达侧视造成的阴影和叠掩效应。遇到孤立的高速度像元先别急着激动看一眼它是不是紧贴着高差大的山地边缘如果是大概率是地形伪影。还有一种情况是冰舌前端的季节性崩解造成的“假信号”。某个像元上一年还在冰舌上这一年冰舌崩解后它变成了开阔水面特征追踪会把这种与冰面纹理无关的变化也当成位移计算出来形成看起来像高速运动的错误向量。年度产品对这种瞬时事件的处理能力有限如果你关注的就是冰舌前端变化建议用更短时间尺度的原始影像数据复核。最后提醒一点下载数据后最好记录一下数据的生产时间和版本说明。NSIDC 会不定期更新数据产品版本同一个 short_name 的 V004 和 V005 之间可能有很大的算法差异不同版本的数据不建议直接混用。如果你在某篇论文里看到别人引用了 V005那就尽量用同一版本复现避免因版本差异产生不必要的争议。我自己的体会是MEaSUREs 格陵兰冰速镶嵌图这种“看起来简单、用起来细碎”的遥感产品最考验人的不是下载和画图而是对数据背后观测原理和误差来源的理解。把 SAR 为什么能干这个活、Landsat 为什么来帮忙、V005 到底改了什么这些问题想清楚遇到奇怪的结果时才不会慌。
