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-RIFTAutonomous 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/svx东向速度分量m/svy北向速度分量m/svx_err东向速度分量的误差m/svy_err北向速度分量的误差m/sv_err速度误差有些版本中给出m/sdate1影像对中较早一景的日期天数相对某个历元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 米网格的单景场景可能是几十到上百 MB960 米网格会小很多。以我自己的习惯先下载 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, cmapviridis, vmin0, vmax100) plt.colorbar(im, labelvelocity (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 测速的问题随时可以循着这些思路去处理。