我第一次打开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) dtypeint16 [D] /ScienceData/Radiance_Scale shape(1,) dtypefloat64 [D] /ScienceData/Radiance_Offset shape(1,) dtypefloat64 [G] /Geolocation [D] /Geolocation/Latitude shape(4056, 421) dtypefloat32 [D] /Geolocation/Longitude shape(4056, 421) dtypefloat32 [D] /Geolocation/ScanTime shape(4056,) dtypefloat64 [D] /Geolocation/SolarZenithAngle shape(4056, 421) dtypefloat32 [D] /Geolocation/SensorZenithAngle shape(4056, 421) dtypefloat32这里的(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之后会得到极其离谱的正负大数。关于FillValueHDF5的数据集属性里一般会写清楚。如果你在attrs里看到f[/ScienceData/Radiance].attrs[_FillValue]读出来的是比如-9999或-32768那么一个比较安全的做法是把radiance数组里小于某个合理下限的值全部设成NaN。别嫌这一步麻烦老数据的FillValue类型五花八门有的填0有的填-9999有的甚至填的是32767int16最大值。多花三十秒做掩膜后面能少掉一堆莫名其妙的色斑。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(projectionccrs.PlateCarree()) ax.set_global() ax.coastlines(resolution50m, linewidth0.5) ax.add_feature(cfeature.LAND, facecolorlightgray) mesh ax.pcolormesh(lon, lat, radiance_ma, cmapmagma, shadingauto, vmin20, vmax70) plt.colorbar(mesh, axax, labelRadiance (W m-2 sr-1 um-1)) plt.title(fNimbus-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(overignore, invalidignore): 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是现存最早的全球夜间高分辨率红外遥感档案之一用它做现代仪器的定标基准那是强人所难但用它做跨六十年的空间格局长序列分析它几乎是唯一的起点。我自己的习惯是每次打开这类数据都先做一遍完整的结构体检和定标链路验证把每一条扫描线的统计特征记录下来再进入正式分析。这个过程听起来琐碎但在老数据上节省的时间绝对值得。
我第一次打开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) dtypeint16 [D] /ScienceData/Radiance_Scale shape(1,) dtypefloat64 [D] /ScienceData/Radiance_Offset shape(1,) dtypefloat64 [G] /Geolocation [D] /Geolocation/Latitude shape(4056, 421) dtypefloat32 [D] /Geolocation/Longitude shape(4056, 421) dtypefloat32 [D] /Geolocation/ScanTime shape(4056,) dtypefloat64 [D] /Geolocation/SolarZenithAngle shape(4056, 421) dtypefloat32 [D] /Geolocation/SensorZenithAngle shape(4056, 421) dtypefloat32这里的(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之后会得到极其离谱的正负大数。关于FillValueHDF5的数据集属性里一般会写清楚。如果你在attrs里看到f[/ScienceData/Radiance].attrs[_FillValue]读出来的是比如-9999或-32768那么一个比较安全的做法是把radiance数组里小于某个合理下限的值全部设成NaN。别嫌这一步麻烦老数据的FillValue类型五花八门有的填0有的填-9999有的甚至填的是32767int16最大值。多花三十秒做掩膜后面能少掉一堆莫名其妙的色斑。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(projectionccrs.PlateCarree()) ax.set_global() ax.coastlines(resolution50m, linewidth0.5) ax.add_feature(cfeature.LAND, facecolorlightgray) mesh ax.pcolormesh(lon, lat, radiance_ma, cmapmagma, shadingauto, vmin20, vmax70) plt.colorbar(mesh, axax, labelRadiance (W m-2 sr-1 um-1)) plt.title(fNimbus-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(overignore, invalidignore): 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是现存最早的全球夜间高分辨率红外遥感档案之一用它做现代仪器的定标基准那是强人所难但用它做跨六十年的空间格局长序列分析它几乎是唯一的起点。我自己的习惯是每次打开这类数据都先做一遍完整的结构体检和定标链路验证把每一条扫描线的统计特征记录下来再进入正式分析。这个过程听起来琐碎但在老数据上节省的时间绝对值得。
