我第一次打开Nimbus-2 HRIR的L1条带数据时,愣了一下。这个HDF5文件干净得像刚生成的现代产品——根属性里写着卫星名、仪器名、版本号,ScienceData组里躺着辐射率,Geolocation组里躺着经纬度和扫描时间。可它记录的却是上世纪六十年代的一天夜里,这台老红外辐射计扫过地球表面的原始痕迹。五十年过去,卫星早就葬身大气层,磁带和胶片也换了三代存档介质,这批数据却被重新整理成了HDF5 V001版本,还在公开分发。这篇文章我想从数据处理的角度,把这个产品从头到尾拆一遍,从仪器背景、条带几何、文件结构,到Python读取、定标换算和实际踩坑,一次说清楚。
如果你正在处理Nimbus系列的HRIR、MRIR、THIR老数据,或者手头有类似“上个世纪气象卫星数据被重新发布为HDF5”的产品,这篇应该能帮你省下不少查资料的时间。就算你只是对HDF5格式怎么组织遥感数据感兴趣,跟着后面的代码走一遍,也能理解为什么这种格式能在一个文件里装下几十年前的“扫描痕迹”。
1. Nimbus与HRIR:为什么1960年代的红外数据今天还是气候研究的硬通货
1.1 从Nimbus计划到高分辨率红外辐射计
Nimbus系列是1964年到1978年间美国发射的一组极轨气象卫星研究计划,前后一共飞了7颗。它的名字取的是拉丁语“雨云”,定位是试验新一代对地观测技术。现在大家熟悉的极轨气象卫星,很多设计思路都能追溯到Nimbus——三轴稳定姿态控制、太阳同步轨道、星上数据存储回放、多波段扫描辐射计,这些概念在那个年代都是第一次被系统性地搬到轨道上验证。
HRIR全称High Resolution Infrared Radiometer,高分辨率红外辐射计,主要搭载在Nimbus-1和Nimbus-2上。在当时的语境下,“高分辨率”是相对此前的低分辨率红外探测而言的。它用的波段是3.4到4.2微米的大气窗口区,星下点分辨率大约8公里,用旋转扫描镜在垂直于卫星飞行方向上来回扫,每扫一行就是一个条带单元。Nimbus-1在1964年8月28日入轨后虽然因为姿态问题只工作了一个月左右,但已经拿到了人类历史上第一批真正可用的夜间红外云图。
Nimbus-2在1966年5月15日发射,在轨正常工作了大约八个月。这段时间里HRIR持续开机,积累了覆盖全球的昼夜红外扫描数据。现在公开分发的L1条带数据,大部分来自这个时期。这批数据的价值在于它是全球范围内最早的系统性高分辨率夜间红外遥感观测。1960年代是地球系统科学从“定性描述”转向“定量观测”的起点,HRIR就是那个转折点上的第一批数据。
1.2 L1条带数据到底能拿来做什么研究
先明确一下L1在这里的含义。遥感卫星数据产品通常分成几个层级:0级是原始码流,L1是经过辐射定标和地理定位后的工程物理量数据,L2及以上是在L1基础上反演出的云参数、海温、地表温度等地球物理产品。这个Nimbus HRIR产品定位在L1,说明它不直接把“温度”给你,而是给出辐射率或者定标后的计数值,同时告诉你每一像元对应的经纬度、扫描时间和质量标记。
研究上这意味着一件事:你拿到的是半成品,但正因为是半成品,你可以按自己的需求做处理,不用被下游产品的反演假设绑住。比如做气候长序列分析的人,最怕不同卫星时期产品之间的系统偏差。如果你想对比1966年HRIR和现在的VIIRS在夜间云量上的差异,用L1自己控制定标流程,可能比直接拿官方L2云产品更灵活。
实际用途大致有三类。第一类是气候基准研究,把1960年代的数据作为卫星观测的最早端点,用于研究云量、海面温度在半个多世纪里的变化趋势,这需要非常小心地做交叉定标,但也正是老数据最有分量的应用。第二类是历史天气过程复盘,比如风暴系统、平流层增温事件,在天气图上结合当时的卫星云图,可以还原天气系统的三维结构。第三类是和再分析资料做对比验证,用HRIR观测来检验ERA系列或NCEP再分析产品在1960年代的重现能力。
2. “条带”是怎么来的:HRIR扫描成像机制与L1产品结构
2.1 一次扫描一行像元:条带数据的几何来源
很多人第一次看到“条带数据”这几个字会有点懵,到文件里一看数据结构更懵:辐射率是一个二维数组,行方向对应扫描线序号,列方向对应一条扫描线内的像元序号,旁边的经纬度也是同样形状的二维数组。为什么不是均匀网格?因为HRIR本来就是沿轨逐条扫描的。
把卫星想象成在轨道上往前飞,辐射计里面的旋转扫描镜垂直轨道方向从左扫到右,每扫一次形成一条垂直于飞行方向的地面条带。探测器在这条扫过的带上连续采样,取样得到一个一个像元,所以一条扫描线就是一行像元。卫星继续往前飞,下一条扫描线紧挨着上一条,逐行累积就得到了一条“带子”——这就是Swath,中文常说的条带或扫描带。由于地球曲率、卫星姿态和扫描角度的综合影响,每条扫描线对应的地面轨迹并不完全平行,远离星下点的位置像元会拉宽变形。所以每个像元都必须单独配一组经纬度,这就是为什么L1数据里经纬度数组和辐射率数组形状完全相同。
HRIR的轨道属于近极地太阳同步轨道,卫星高度大约700到900公里的区间(不同轨道段有差异),星下点分辨率约8公里,扫描幅宽在千公里量级。这意味着卫星每绕地球一圈,在地面形成一条从北到南或从南到北的宽条带,昼夜交替地掠过不同区域。白天波段里3.4到4.2微米会混入反射太阳光,所以HRIR的夜间数据在热辐射分析上更“干净”,这也是很多研究只用夜轨数据的原因。
2.2 L1产品里装进了一个轨道的信息
HDF5版本的L1产品在组织方式上遵循了很典型的“科学数据集+辅助定位数据”思路。你打开文件,通常能在根目录下看到几大块内容。第一块是元数据属性,以Attribute的形式挂在根组或某个Metadata组上,里面记录卫星名、仪器名、数据级别、产品版本、轨道号、数据起止时间等。第二块是辐射测量主体,可能是原始DN(数字计数值),也可能是经过定标后的辐射率,单位一般是瓦每平方米每球面度每微米。第三块是地理定位数据集,包括每个像元的经纬度、扫描时间、太阳天顶角、卫星天顶角这些。第四块是质量标记,有些产品会把每条扫描线或每个像元的质量状态单独放一个数据集,比如正常、缺行、标定异常、边缘几何畸变等。
这种“主体数据+每条像元定位”的结构,是条带数据的标准格式。和常见的L2格点化产品不同,L1条带数据里的行列号不代表经纬规则网格,必须通过经纬度数组才能把数据投到地图上。后面做可视化的时候,这一步很容易出错——直接拿数组行列当坐标画图,画出来的是一张扭曲变形的图,贴在球面上毫无意义。
3. 用Python解剖一个Nimbus HRIR L1 HDF5文件
3.1 先学会把文件结构“看”出来
不管数据产品文档写得多详细,拿到文件第一件事永远是把它当“洋葱”一层层剥开看。HDF5的优势就在于自描述性,即使没有任何PDF文档,只靠HDF5文件本身,你也能搞清楚里面有哪些组、哪些数据集、每个数据集的形状和类型、单位属性写的是什么。
推荐用h5py作为入门工具。环境准备很简单,用conda或者pip安装:
conda install h5py numpy matplotlib cartopy我习惯先把整个文件结构打印出来。用一段很短的脚本:
import h5py with h5py.File("NIMBUS2_HRIR_L1_1966_XXX_v001.h5", "r") as f: def show(name, obj): if isinstance(obj, h5py.Dataset): print(f"[D] {name} shape={obj.shape} dtype={obj.dtype}") for key, val in obj.attrs.items(): print(f" attr: {key} = {val}") else: print(f"[G] {name}/") f.visititems(show)这段代码做的事情很朴素:DFS遍历文件里所有对象,Group打上[G]标记,Dataset打上[D]标记,并把每个数据集的shape、dtype和属性打印出来。实际打开这批数据的时候,你会看到类似这些名字的路径:
[G] / [G] /Metadata [G] /ScienceData [D] /ScienceData/Radiance shape=(4056, 421) dtype=int16 [D] /ScienceData/Radiance_Scale shape=(1,) dtype=float64 [D] /ScienceData/Radiance_Offset shape=(1,) dtype=float64 [G] /Geolocation [D] /Geolocation/Latitude shape=(4056, 421) dtype=float32 [D] /Geolocation/Longitude shape=(4056, 421) dtype=float32 [D] /Geolocation/ScanTime shape=(4056,) dtype=float64 [D] /Geolocation/SolarZenithAngle shape=(4056, 421) dtype=float32 [D] /Geolocation/SensorZenithAngle shape=(4056, 421) dtype=float32这里的(4056, 421)意味着文件里有4056条扫描线,每条扫描线421个像元,对应的是一次相对完整的轨道片段。Radiance用的是int16,说明数据为了节省空间做了定点缩放,真实辐射率要靠Radiance_Scale和Radiance_Offset还原。这种存储策略在HDF5遥感数据里极其常见,务必养成“先看Scale和Offset再算物理量”的习惯。
3.2 读取辐射、经纬度与质量标记
确认结构以后,下一步就是把需要的数组载入内存。这里有几个值得注意的技巧。HDF5文件往往很大,如果一次把整个Radiance读进内存,4056×421的int16数组其实还好,但如果以后遇到更长的轨道片段或者更宽的扫描线,就要考虑按块读取。h5py的Dataset本身支持类似numpy的切片操作,f["/ScienceData/Radiance"][1000:2000, :]这样就能读到指定块,不会把整个文件全载入。
我的推荐流程是,先读全局属性确认轨道号和时间范围,再读定位数据,最后再读主体辐射数据。这样如果经纬度范围跟你的研究区对不上,可以直接换文件,省得重复读大数组。
import h5py import numpy as np with h5py.File("NIMBUS2_HRIR_L1_1966_XXX_v001.h5", "r") as f: meta = {key: val for key, val in f.attrs.items()} print(meta) lat = f["/Geolocation/Latitude"][:] # (scan_line, pixel) lon = f["/Geolocation/Longitude"][:] scan_time = f["/Geolocation/ScanTime"][:] # (scan_line,) solzen = f["/Geolocation/SolarZenithAngle"][:] senzen = f["/Geolocation/SensorZenithAngle"][:] rad_raw = f["/ScienceData/Radiance"][:] # (scan_line, pixel) scale = f["/ScienceData/Radiance_Scale"][0] offset = f["/ScienceData/Radiance_Offset"][0] radiance = rad_raw.astype(np.float64) * scale + offset有一点要特别提醒:Radiance_Scale在HDF5里一般是一个一维长度为1的数据集。读取时加[0]把它取成标量,后面才不会出现数组形状问题。定标公式本身很简单——物理量等于原始计数值乘以比例因子再加上偏移量——但在老数据里,你很可能会遇到Radiance字段里混着填充值(FillValue),这些填充值必须提前掩膜,否则乘以scale之后会得到极其离谱的正负大数。
关于FillValue,HDF5的数据集属性里一般会写清楚。如果你在attrs里看到:
f["/ScienceData/Radiance"].attrs["_FillValue"]读出来的是比如-9999或-32768,那么一个比较安全的做法是把radiance数组里小于某个合理下限的值全部设成NaN。别嫌这一步麻烦,老数据的FillValue类型五花八门,有的填0,有的填-9999,有的甚至填的是32767(int16最大值)。多花三十秒做掩膜,后面能少掉一堆莫名其妙的色斑。
3.3 把整轨数据画成一张夜间云图
拿到经纬度和辐射率之后,最直观的质量检查就是画图。条带数据的正确画法是用pcolormesh把经纬度直接作为坐标轴,而不是用行列号当坐标。比例尺大约对应几千公里宽,用普通的PlateCarree投影就能看全貌。
import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature mask = radiance < 0 radiance_ma = np.ma.masked_where(mask | ~np.isfinite(radiance), radiance) fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=ccrs.PlateCarree()) ax.set_global() ax.coastlines(resolution="50m", linewidth=0.5) ax.add_feature(cfeature.LAND, facecolor="lightgray") mesh = ax.pcolormesh(lon, lat, radiance_ma, cmap="magma", shading="auto", vmin=20, vmax=70) plt.colorbar(mesh, ax=ax, label="Radiance (W m-2 sr-1 um-1)") plt.title(f"Nimbus-2 HRIR L1 Swath - {meta.get('StartTime', '')}") plt.show()如果pcolormesh画得太慢,尤其是整轨几百万个像元时,可以用ax.scatter配ravel()后的经纬度和辐射率来降采样绘制。但pcolormesh一般够用。真正容易出现的问题是投影范围——如果lon和lat数组里包含了南极和北极的扫描段,面的拼接可能会有畸变,观察上不影响判断。画出来的图上应该能看到一条条淡亮的云带,云顶在3.7微米窗口的辐射率通常低于地表,所以云区颜色偏暗;地表尤其是夜间陆地或海面温度较高,亮一些。如果整幅图出现条纹状异常,或者某几行全黑全白,基本可以判定是定标或填值处理出了问题,可以直接定位到ScanTime对应的行,用质量标记排查。
4. 从DN到亮温:L1数据里的定标与物理量换算
4.1 辐射定标与普朗克反函数
不少做天气气候分析的人拿到L1辐射率之后,还是更习惯转成亮温来用。亮温的本质是:如果一个黑体发出的辐射率等于你观测到的辐射率,那么这个黑体应该是什么样的物理温度。它只是一个辐射率到温度的数学映射,不代表目标的真实热力学温度。
HRIR的通道在3.4到4.2微米,通常取中心波长大约3.8微米来做单色近似。已知观测辐射率L,用普朗克函数反解出来的亮温公式是:
h = 6.626e-34 k = 1.381e-23 c = 2.998e8 lambda_c = 3.8e-6 def radiance_to_brightness_temperature(L): # 处理无效值 L = np.maximum(L, 1e-20) # 普朗克反函数 with np.errstate(over="ignore", invalid="ignore"): Tb = (h * c / (lambda_c * k)) / np.log(2 * h * c**2 / (lambda_c**5 * L) + 1.0) return Tb bt = radiance_to_brightness_temperature(radiance_ma)在这个公式里,注意单位全部要用国际单位。辐射率的单位如果写的是W m-2 sr-1 um-1,需要先换算成W m-2 sr-1 m-1,也就是乘以10的6次方;如果文档里本身就是W m-2 sr-1 m-1,那就直接算。方向搞反的人不在少数,我有一回就是漏了换算,反演出来的亮温整体偏了几十K,白花了一个下午排查。
转换后的亮温在夜间地表场景一般落在250K到300K之间,云顶亮温可能降到220K甚至更低。如果出现200K以下的大片数值,要么是数据处理有问题,要么是卷云或深对流云顶——后者是有物理意义的,不能直接当坏值剔除。这块要结合扫描角度和太阳天顶角一起判断。
4.2 老仪器的质量标记与“看起来正常”的隐患
老数据的定标和当代卫星有很大区别。Nimbus时代的HRIR在轨定标依赖星上的基准源,稳定的黑体和空间冷空观测,但那个年代的星上定标设备精度有限,加上仪器在轨长期工作后响应会发生漂移,所以同一个仪器在不同阶段获取的数据,辐射率可能存在百分之几到百分之十几的偏差。L1产品通常会把这些偏差部分考虑进去,在质量标记上反映出来。
质量标记一般有两种粒度。一种是整轨或整段时间的状态标记,比如“仪器温度正常”“定标源正常”“回放数据无丢失”;另一种是每条扫描线的质量状态,老数据更常见的是扫描线级异常——某一行扫描数据丢失,那一行像元会被填充成FillValue。当你在图像上看到一条或多条黑色横线时,大概率不是仪器抽风,而是该行数据缺失。
老数据的隐蔽问题在于:有些行数据不是完全缺失,而是部分像元被污染。比如扫描镜在特定角度上的信号抖动,或者姿态控制误差导致的地理定位偏差。单看Radiance数组,这些行和正常数据差别不大,但做定量反演时会引入噪声。所以我的习惯是拿到文件后先做一个“按行统计”:计算每条扫描线的均值、标准差和有效像元比例,把哪些明显偏离的扫描线编号记下来,在后续分析中单独处理或排除。这种简单的一维统计图能暴露很多肉眼看不出来的问题。
4.3 存档机构为什么选HDF5,而不是GeoTIFF或NetCDF
接触这批数据的时候,我也想过一个问题:为什么存档机构最终选择用HDF5来重发布这批60年代的模拟转数字数据,而不是用GeoTIFF或者NetCDF?其实答案藏在数据本身的特征里。
L1条带数据不是一个单波段的平面图。它一个文件里同时容纳了二维辐射数组、每像元经纬度、扫描时间、观测几何角、质量标记、全局元数据。GeoTIFF最适合的是规则的二维栅格,要硬塞条带数据也能塞,但时间、角度、质量标记这些辅助信息就无处安放,只能塞进一堆自定义的TIFF Tag里,可读性和通用性都很差。NetCDF在气象界用得非常多,本质上和HDF5同属自描述科学数据格式,但在上世纪九十年代到本世纪初,美国NASA地球观测系统那一套数据分发体系已经深度绑定HDF/HDF-EOS。大量EOS历史数据都存放在HDF5这个技术栈上。Nimbus数据作为早期地球观测数据的一部分,重发布为HDF5,既方便统一工具链,也和后续的MODIS、VIIRS等数据格式保持一致。
这批数据在数据工程上和现在的机器人数据集有很奇妙的共通点。Hugging Face的lerobot项目把机器人示教数据封装成HDF5,原因和遥感如出一辙——多维数组、异构元数据、随机访问切片,一套容器全部解决。自动驾驶领域的MCAP格式走的是偏流式记录的技术路线,而HDF5则更适合“一次写入、长期存档、按需分析”的科学场景。这说明HDF5作为“同一语言”已经跨出了传统的科学计算范围,遥感、机器人、自动驾驶都在用它组织复杂的传感器数据。所以学会读HDF5,收益不止在遥感这一个领域。
5. 复盘:处理这批数据时我踩过的坑和反复确认的资料
5.1 时间是这批数据里最容易理解错的东西
做时间处理之前,先确认时间基准。不同数据产品的ScanTime字段含义差别巨大,有的存的是自某一天的0点以来的秒数,有的存的是远地点的相对时间,有的存的是真实的UTC字符串。Nimbus时代的数据还有一个很麻烦的地方:当年轨道的传递时间精度不高,和现在GPS授时的卫星数据完全不同。ScanTime的绝对精度可能只有分钟级或者几十秒级,这在做时间匹配时会造成几十公里的定位误差。
我踩过的坑是:想直接把ScanTime和现代的再分析数据按小时做匹配,结果发现同一条扫描线上的所有像元都被“平均”到了一个时间点上,导致扫描线内部几十秒的时间梯度被抹平了。后来我改用ScanTime插值到每个像元,虽然绝对误差依然存在,但至少保留了条带内部的时间梯度。建议你在处理前先打印几条扫描线的ScanTime,看看相邻行的间隔是多少秒、有没有跳变。如果间隔忽大忽小,说明当年数据回放过程有丢帧,后面做叠加分析时要注意。
5.2 填充值、比例因子和HDF5的切片读取
填充值和比例因子这两个坑我已经在前面反复提过,这里再强调一遍:老数据的FillValue没有任何国际标准,不同文件、不同数据集可能各不相同。有的存成-9999,有的存成-32768,还有的直接把缺行填成0。如果你直接用0当正常值参与计算,定标之后辐射率为负几百甚至负几千,图上一片“黑洞”。正确做法是每次打开文件都先把每个数据集的_FillValue属性打印出来,建立一个掩膜数组,后续所有物理量转换都带着这个掩膜走。
HDF5本身的存储方式也值得一说。L1条带数据虽然不大,但同系列的HDF5重发布产品里也有几十GB的版本,那种就不能无脑加载。HDF5支持chunked存储和压缩,合理的做法是先看数据集的chunks属性,然后结合你的研究区切片读取。比如你只要南纬60度以南的数据,可以先根据经纬度数组判断扫描线范围,再按行切片读取Radiance。这样文件I/O性能会高很多。还有一点,如果读取非常频繁,可以考虑用h5py的高层API设置读取缓存,或者用dask.array包装HDF5数据集做惰性加载,处理起来体验会好很多。
5.3 和现代卫星数据交叉对比前,先想清楚波段差异
最后说一个最容易被忽视的问题:波段差异。HRIR的中心波长是3.8微米附近,而现在的极轨气象卫星常见的红外窗口通道是10到12微米。同一个云顶、同一个地表,在3.7微米测出来的亮温和在11微米测出来的亮温根本不相等,因为普朗克函数的斜率在这个区间不是线性的,而且云滴和冰晶在这两个波段的吸收特性和发射率也有差异。
所以如果你想拿Nimbus HRIR和MODIS/VIIRS做历史对比,绝对不能直接“亮温对亮温”画散点图。正确做法是,要么用辐射传输模式把两个波段的观测统一到同一个物理量上,要么只用其中一个波段作为“晴空、云量、相对温度梯度”的定性参考,而不是绝对温度。我在和ERA5再分析对比时,采取的办法是把HRIR的观测区域缩小到海面上,用3.8微米的亮温结合海面发射率模型反演等效海表温度,再和再分析的海温做相对变化比较。这样做的结果趋势是能对上的,但绝对偏差仍然可能达到1到2K,不能当作仪器级精度的验证数据。
老数据就是这样——你得先接受它那个年代的不完美,再去挖掘它不可替代的那部分价值。HRIR是现存最早的全球夜间高分辨率红外遥感档案之一,用它做现代仪器的定标基准那是强人所难,但用它做跨六十年的空间格局长序列分析,它几乎是唯一的起点。我自己的习惯是每次打开这类数据都先做一遍完整的结构体检和定标链路验证,把每一条扫描线的统计特征记录下来,再进入正式分析。这个过程听起来琐碎,但在老数据上节省的时间绝对值得。