1. 数据集是什么:先看 MEaSUREs 和 ITS_LIVE 这两个词
经常有同行问我:“冰川流速数据你们都用哪套?”我基本上第一反应就是 ITS_LIVE。ITS_LIVE 全称是 Inter-mission Time Series of Land Ice Velocity and Elevation,属于 NASA MEaSUREs 计划产出的一个长期数据集。MEaSUREs 的全称是 Making Earth System Data Records for Use in Research Environments,翻译成人话就是“让地球系统数据变成能直接用于科研的标准产品”。MEaSUREs 资助了很多这类“数据记录”项目,ITS_LIVE 是里面专门做陆地冰速度和高程时间序列的一个。
这个数据集最核心的特点是:它把 1984 年以来的 Landsat 4/5/7/8 光学影像,两两配对,用影像匹配技术算出冰川和冰盖的表面速度,做成了一套覆盖全球、时间跨度超过三十年的速度场产品。我手里拿到的这一版是 V001,也就是第一个正式发布版本。很多研究里引用的 ITSLIVE 数据集,指的就是这个。
它能解决什么问题呢?很直接:过去你想研究某条冰川在某个时间段流得有多快,你得自己去下载几十景 Landsat 影像,自己跑特征匹配程序,还要处理配准误差、投影转换、云遮挡,一套流程下来几个月就过去了。ITS_LIVE 把这些活全替你干了,而且处理得很规范,你拿到手就是一个带坐标、带误差、带日期的速度栅格。对于搞冰川动力学、冰湖溃决风险、冰川物质平衡、冰盖接地线变化的团队来说,这几乎是“开箱即用”的基础数据。
适合谁用?三类人最受益:一是做区域冰川变化监测的研究生和科研人员,他们需要快速得到大范围速度场;二是做冰盖或冰川建模的人,需要时空连续的流速输入;三是做工程或灾害评估的从业者,比如评估冰湖溃决、冰川跃动风险时,速度序列是最直观的动态指标。这篇内容我尽量按照“先讲原理、再拆数据、然后给实操、最后讲坑”的顺序写,你照着试一遍基本就能上手。
2. 影像对测速的基本逻辑:两张图怎么就变成了速度场
2.1 核心算法:不是“像素位移”,而是“特征追踪”
很多第一次接触 ITS_LIVE 的人会以为:两张影像的时间间隔已知,把同一块冰在前后两张影像里的位置找出来,位移除以时间就是速度。大方向没错,但实际操作远不是“找同一块冰”这么简单。
冰川表面不是均匀的白板,它上面有裂缝、冰碛、岩石碎屑、融水通道等纹理。Landsat 影像虽然只有 30 米分辨率,但在大冰川上这些纹理足够作为“特征点”。ITS_LIVE 使用的核心算法是 auto-RIFT(Autonomous Repeat Image Feature Tracking),它会在两张影像上做金字塔式的多尺度匹配:先在低分辨率上粗略找到位移量,再逐步加密到原始分辨率上精确定位。匹配的相似度指标常用的是归一化互相关(NCC),也就是在目标影像上取一个窗口,在参考影像上滑窗寻找与之最相似的窗口。
这里有一个很关键的细节:速度场里的每个格网值,不是单个像素的位移,而是一个窗口内所有特征的统计结果。所以 ITS_LIVE 输出的速度网格分辨率(120 米、240 米、480 米、960 米)远大于 Landsat 原始像元 30 米。你用 30 米影像,却得到 120 米甚至更粗的网格,这是正常的,不是数据缩水了。窗口匹配天然需要空间范围,窗口越大,对纹理的要求越低,但空间分辨率和细节捕捉能力会下降。ITSLIVE 同时输出多个尺度,就是为了让用户在“覆盖率”和“分辨率”之间做选择。
2.2 误差和不确定性是从哪里来的
速度产品的误差项,不是事后加的,而是算法每一步累积下来的。主要来源有三个。
第一是配准误差。Landsat 影像本身做了地理定位,但不同轨道、不同年份的影像之间总会存在系统性偏移,这个偏移可能来自地形、轨道姿态和 DEM 误差。auto-RIFT 在匹配前会先做“稳定地面”的配准估计,把非冰川区的位移中位数当作整体配准偏差扣除掉。
第二是特征失相干。如果在两景影像之间,冰雪表面发生了强烈变化——比如新雪覆盖、消融季大量融水、跃动冰川发生剧烈形变——那么窗口里的纹理在两景之间已经对不上了,相关性很低,匹配到的位移就是噪声。ITSLIVE 会通过信噪比或相关系数阈值来标记这种低质量结果,你在数据里看到的 NaN 或大误差值,很多就是这么来的。
第三是时间间隔的不合理。影像对间隔太短,比如只有 16 天,而冰川在这个周期内可能只移动了几米甚至不到一米,这时位移和配准误差在同一个量级上,速度值的相对误差就非常大。反过来,如果跨越好几年,特征很可能完全失相干。所以用这套数据时,你首先要养成一个习惯:看每个像元的速度时,同时看它对应的影像对时间间隔和误差值,不能只看速度大小。
2.3 为什么偏偏是 Landsat
可能有人会问:现在有 Sentinel-2 甚至更高分辨率的商业影像,为什么不都用这些?原因很简单:时间连续性。ITSLIVE 的核心目标是提供 1985 到现在的长序列,Sentinel-2 2015 年才发射,之前的历史只能靠 Landsat。Landsat 4/5 从 1982 年就开始提供 TM 影像,Landsat 7 在 1999 年上天,Landsat 8 在 2013 年接棒,加上 2021 年发射的 Landsat 9,光谱能力和成像质量虽然代际不同,但经过处理后能形成一条相对一致的长序列。
这是 ITS_LIVE 最大的价值所在:它不是一套只覆盖几年的“快照”,而是跨越了至少三四个传感器、三十多年的“时间档案”。你可以用它做长时间尺度的趋势分析,比如对比 1990 年代和 2010 年代同一区域冰川流速的整体变化,这是没有长序列数据时完全做不到的。
3. V001 数据文件里到底有什么:拿到一个 NetCDF 先看这些
3.1 V001 版本与早期预发布版本的区别
V001 是正式发布版本。在它之前,ITSLIVE 网站上有过若干预发布数据,主要是给内部合作者用,算法还在迭代。V001 最大的变化是处理流程固定下来了,输入数据源统一到了 Landsat Collection 2 级别,并且对全球所有区域采用了同一套处理参数。这意味着不同区域、不同时间的数据之间具备可比性,这对做全球尺度的统计非常关键。
我实际用下来还有两个感受:第一,V001 对南极高纬度区域做了更细致的投影处理;第二,文件里自带的误差字段比早期版本更靠谱。早期版本有些地方的误差值是同一个常数,V001 你会看到误差场随地形、纹理和季节变化,这才是真实情况。
3.2 关键变量:vx、vy、v 和它们的误差
打开一个 ITS_LIVE NetCDF 文件,核心变量名不算多,但每个都要理解清楚。常见的是这样几组:
| 变量 | 含义 | 单位 |
|---|---|---|
| v | 速度大小,即合成速度 | m/s |
| vx | 东向速度分量 | m/s |
| vy | 北向速度分量 | m/s |
| vx_err | 东向速度分量的误差 | m/s |
| vy_err | 北向速度分量的误差 | m/s |
| v_err | 速度误差(有些版本中给出) | m/s |
| date1 | 影像对中较早一景的日期 | 天数,相对某个历元 |
| date2 | 影像对中较晚一景的日期 | 天数,相对某个历元 |
这里有个容易踩坑的点:ITSLIVE 变量的时间基准不是标准 Unix 时间,很多文件里 date1 和 date2 是相对 2000-01-01 之类的参考日期计算的天数偏移。你用 xarray 读出来后,必须先搞清楚文件的 time_units 属性,再转成真正的日期,否则画时间轴时就会差一大截。
坐标方面,全球陆地冰区域的 V001 产品通常采用极方位立体投影或 EASE 网格,不同区域对应的投影不同。文件里的 x 和 y 一般已经是投影坐标,单位是米。所以你在做研究时,如果要和其他来源的矢量数据叠加,不要忘了先统一坐标系。
3.3 网格、文件命名和每个文件的时间含义
V001 的全球陆地冰产品按区域分块发布,常见网格间距包括 120 米、240 米、480 米和 960 米。你下载时看到的文件名里通常会包含传感器类型、时间信息和版本号。以典型的 ITS_LIVE 文件名为例,里面会有一段日期,说明这一对影像跨越的时间范围。一个文件就是一个影像对的速度场,不是一个月的平均速度场。
这条很重要,我多说两句。很多人拿到的第一个 ITS_LIVE 文件,看到“日期 2015”,以为这是 2015 年的平均速度图。其实不是,这是“2015 年某两个日期之间”的瞬时平均速度。如果这两个日期分别跨了冬夏,速度会表现出很强的季节性差异。所以你要做“年度平均速度”或“月平均速度”,绝不能直接用一个文件代表一段时间,而是要把多个影像对的结果聚合起来。
4. 获取数据与预处理实操:从下载到能画图
4.1 下载途径与筛选参数
现在获取 ITS_LIVE 数据主要有两个渠道:一是 NSIDC DAAC 的 Earthdata Search 界面搜索 “ITS_LIVE”,二是 ITS_LIVE 项目官网的交互式地图浏览工具。前者适合批量下载,后者适合你先快速看一眼区域覆盖情况。
如果你用 NASA Earthdata Search,搜索到数据集后,需要先用注册的 Earthdata 账号登录。建议按区域和日期范围筛选,一次不要下载整个全球产品。文件体积不一定会很小,120 米网格的单景场景可能是几十到上百 MB,960 米网格会小很多。以我自己的习惯,先下载 960 米或 480 米网格做探索性分析,确认研究区域和时间段没问题后再去拿 120 米精细数据。
提示:如果连续下载大量文件,建议用 Python 脚本配合 Earthdata 的 token 做批量下载,而不是在网页上一个一个点。文件数量多时,手动下载很容易漏。
4.2 用 xarray 打开一个 NetCDF 并快速出图
拿到文件后,推荐直接用 Python 的 xarray 打开,配合 rioxarray 处理空间参考。以下是最基本的读取和可视化流程。
import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs # 打开文件 ds = xr.open_dataset("ITS_LIVE_...nc") # 查看全局属性,确认单位、投影和时间基准 print(ds.attrs) # 看一眼速度变量统计 print(ds["v"].values.min(), ds["v"].values.max()) # 简单出图 fig, ax = plt.subplots(figsize=(8, 8)) im = ax.imshow(ds["v"].values, cmap="viridis", vmin=0, vmax=100) plt.colorbar(im, label="velocity (m/s)") ax.set_title("ITS_LIVE velocity field") plt.show()注意,大多数 ITS_LIVE NetCDF 里的空间维是垂直的 y 和水平的 x,而 imshow 默认把第一维当作行,第二维当作列,所以直接用 ds["v"] 画图时方向一般没问题。只是要记住,因为投影不是经纬度,所以你直接打印坐标的时候,x 和 y 单位是米。
4.3 多文件合成:如何得到区域年度合成速度场
单个影像对的速度场会有缺口:可能是云遮挡,可能是冬季纹理太弱,也可能是时间间隔太短导致误差过大。实际使用中,我们要把一段时间内的多个影像对合成一个“代表速度场”。
常用的做法是按年份分组,对一年的所有影像对像元取中位数。不要取平均值,因为个别影像对的错误匹配会产生极端速度值,中位数能有效抑制这些尖峰。具体可以这样做:先把所有文件的时间字段转成真实年份,然后用 xarray 的 concat 拼接,再 groupby 年份后取中位数。我强烈建议在合成前先把误差大的像元清除,常见阈值是去掉 v_err 大于某个值(比如 20 m/yr)的像元。这个阈值没有绝对标准,要根据区域流速和影像对密度来定。
5. 实际案例:提取一条冰川的速度时间序列
5.1 从速度场到冰川中轴速度序列
假设你研究域内有一条冰川,想看看它过去 20 年的流速变化。ITSLIVE 的做法很直接:对每一年合成速度场均提取冰川中轴线上的速度值,然后画时间序列。
具体步骤我总结为四步。
第一步,准备冰川中轴线或中心线矢量,最好基于最新的冰川边界数据绘制。第二步,对每个年份合成速度场,在用 GDAL 或 rioxarray 按矢量裁剪后,提取线上像元。第三步,把线上所有像元的速度取中位数,代表该年该冰川的特征速度。第四步,把所有年份的结果串起来,画散点或折线图。
这一步看起来简单,但有一个大坑是冰川面积不大时,速度场上的像元很少。如果你的冰川宽度不足 1 公里,120 米网格上可能只有 5 到 8 个像元,再抠掉云遮挡和低质量像元,剩下的可能只有两三个。这时中位数会很抖。更稳的方法是:沿着中轴线做 500 米缓冲区,把缓冲区内的速度全部纳入统计。这样能显著提高样本量,代价是小冰川两侧速度差异会被略微抹平。
5.2 跃动冰川场景:为什么单一年份可能误导
ITSLIVE 数据在研究跃动冰川时特别好用,但也要特别小心。跃动冰川在活跃期速度可以在几个月内从每年几十米暴涨到每年几公里。如果你手里只有一两个影像对,恰好落在跃动爆发前后,你的“年度速度”可能会严重偏离真实过程。
我的建议是:一旦发现速度序列里有异常陡升陡降,不要直接当传感器噪声,先回去把当年的原始影像对文件全部找出来逐一看。ITSLIVE 网站上的交互工具支持按时间滑块浏览单个影像对,这是排查异常值最快的方式。如果几个影像对在相近时间段内都显示出高速,那很可能是真的跃动;如果只有一个点突出而周围时间点正常,那大概率是匹配错误。
6. 常见问题与排查技巧实录
| 问题 | 可能原因 | 排查与解决 |
|---|---|---|
| 速度场大片 NaN | 云覆盖、冬季新雪覆盖导致纹理失相干 | 改用其他月份影像对,或用低分辨率网格尝试 |
| 速度值明显偏大(如几千 m/s) | 影像对匹配错误或配准未完全消除 | 检查相关系数,删除误差过大的像元再合成 |
| 时间轴显示异常 | date 字段基准时间未转换 | 查看 attrs 中 time_units,用 cftime 正确解码 |
| 文件打开后投影单位是米无法叠加经纬度矢量 | 投影坐标系不同 | 用 rasterio 读取 CRS,用to_crs统一矢量坐标系 |
| 遥感和实地速度对不上 | 影像对瞬时速度与野外年均值定义不同 | 明确数据时间口径,尽量用同一季节或作年合成 |
6.1 影像对时间间隔的判断经验
做时间序列分析时,我自己的经验是优先选择时间间隔在 32 到 96 天之间的影像对。为什么?间隔太短(16 天)时,如果冰川流速较慢(比如每年只有几十米),位移就很小,误差占比太高;间隔超过 96 天往往进入不同季节,雪盖和消融状态差异大,失相干风险升高。当然这只是一般原则,对于速度极快的冰流或跃动期冰川,16 天反而更有优势。你不该用一个固定规则死套,而是要结合研究区域的量级去筛选。
6.2 低质量像元过滤阈值怎么定
很多新手喜欢照搬文献里的固定阈值,例如把误差大于 30 m/yr 的点全部过滤。这有时管用,有时会误伤。更好的做法是打印出研究区域速度场的误差分布直方图,看看哪里有明显分界。如果大多数像元误差在 5 m/yr 以内,个别像元误差到了 50 m/yr,那显然可以去掉;但如果整个区域误差都在 20 到 30 m/yr 之间,你硬要卡 20,可能大半个研究区就没了。误差阈值应该由数据本身决定,而不是随手填一个数。
6.3 与其他传感器数据拼接时注意的系统差
现在 ITLIVE 还发布了基于 Sentinel-1 雷达和 Sentinel-2 光学影像的速度产品。很多研究想把 Landsat 和 Sentinel-2 序列拼接起来延长或加密时间线。这里要提醒:不同传感器的速度场之间存在系统性差异,尤其是在光学影像上,传感器特性、处理流程和匹配算法参数都可能略有不同。拼接前最好在重叠时间段内做逐像元比较,对差值进行线性拟合,再决定是否要加偏差项。不要想当然认为“都是 ITS_LIVE 产品就完全一致”,跨传感器序列的一致性检验在任何时间序列研究里都是必须具备的一步。
7. 一些实操中的心得体会
最后聊几点我在用这套数据时积累的个人感受。
第一,下载数据前先花半小时看网站上的区域覆盖图。全球陆地冰产品虽然覆盖广,但某些山区小冰川因为目标太小或纹理太弱,实际有效数据可能稀少。提早看覆盖情况能避免你下载一堆“看上去有冰、统计出来全 NaN”的数据。
第二,ITSLIVE 的 NetCDF 文件属性信息非常值得仔细读。里面记录了影像对来源、处理软件版本、投影参数、时间基准等细节。我在写方法部分时,很多重要信息直接从 attrs 里复制即可,不需要再去翻处理文档。
第三,做长时间序列时,宁可牺牲分辨率也别牺牲重叠覆盖度。如果你用了 120 米网格,但每年只有一两个有效影像对,时间序列会很破碎;改用 480 米网格后,虽然空间细节少了,但每年可能能得到 5 到 10 对有效数据,序列的连续性和统计稳健性会好得多。对于大多数区域尺度的趋势研究,480 米网格完全够用。
第四,别怕处理速度场上的“异常箭头”。所有基于影像匹配的速度产品都会有少量错误匹配,这是无法彻底消除的。关键是你的后处理流程里有没有“误差过滤 + 中位数合成 + 关键时段人工核查”这三道工序。把这套流程走顺了,不管面对哪套影像测速产品,你都能做出一份可信的结果。
ITSLIVE 这套数据贵在“长”和“全”,但再好的数据也得靠合理的处理才能变成科学结论。希望上面的内容能让你少走几步弯路,如果之后再碰到 Landsat 测速的问题,随时可以循着这些思路去处理。