简介：星载激光雷达CALIPSO-VFM产品数据读取与显示Matlab代码，适用于研究地球遥感、气象卫星探测数据的科研人员和学生，重点解决垂直特性产品从HDF文件读取到图形化展示的实际问题。压缩包共8个文件，由6个Matlab脚本和2个HDF示例数据组成，总体仅1.23MB，体积轻巧；脚本覆盖主程序、图形绘制、动态颜色表创建、类型判断、行块转换等核心环节，两个示例数据来自2019年3月不同时刻的标准产品，可直接运行查看真实读取与显示结果，便于理解数据结构和处理流程。目前已有1501人学习或下载，对刚接触星载激光雷达数据处理的人来说是一份容易上手的入门范例。这套代码模块划分清晰，读者既能快速复现VFM产品的显示效果，也能掌握分块读取、分类着色等关键思路，进而扩展针对不同区域或时段的星载激光雷达分析，适合在教学、科研或项目预研中使用。轻量体积和模块化设计还可减少调试成本，方便在具体课题中直接套用与二次开发。

1. 星载激光雷达Calipso-VFM产品最常被问的一句话：44字节到底是什么

第一次打开 Calipso-VFM 产品的人，十有八九会愣一下：说好的激光雷达图像，读进来不是一张剖面强度图，而是一堆 44 字节的掩膜码。这不是数据错了，而是 VFM（Vertical Feature Mask）本身就是分类产品——CALIOP 激光雷达把每条视场柱里的目标识别成云、气溶胶、晴空，再按类型给出一串编码。做气溶胶和云研究的人用它做类型判据，做数据质量控制的人用它剔除无效剖面，做算法验证的人拿它当标签真值。这篇文章解决的就是从 HDF4 文件里把这些掩膜读出来、解码、画成能放进报告里的垂直剖面图，并把最容易翻车的五个细节一次讲透。

2. 先说清楚VFM是什么：从L1后向散射到特征掩膜的三个关键概念

2.1 VFM在Calipso产品链里的定位：不是强度图，是分类图

Calipso 卫星搭载的 CALIOP 激光雷达，向下发射 532nm 和 1064nm 两束激光，接收大气后向散射信号。一级产品（L1）给的是沿轨-高度剖面上的衰减后向散射强度，能看到云层和气溶胶层的轮廓，但你无法直接从强度值判断“这层是沙尘还是烟”。VFM 属于二级产品，它把每一条视场柱、每一个高度 bin 做特征分类，输出一个掩膜：这个点是不是特征，特征类型是云还是气溶胶，子类型是什么，分类置信度如何。

所以读 VFM 时心态要先切换：你手上不是一张遥感影像，而是一张标注图。每个格点的数值不是物理量，而是分类代码。这个区别几乎决定了后面所有处理流程的写法——不能用连续色带渲染，不能做平滑插值，不能直接把数值当强度用。

2.2 文件命名规则与需要读取的HDF4变量

VFM 产品以 HDF4 格式分发，文件名形如CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf。文件名里值得留意两处：一是ValStage1，它表示处理阶段，不是数据质量等级，不能看见 Stage1 就以为是中间产品；二是日期时间字段，它表示轨段起始的 UTC 时间，处理多轨数据时排序别按本地时间排。

打开文件后，真正要读的变量不多，核心是这些：

Feature_Flight_2048_Column_44_Byte这个名字很长，拆开看：2048 表示沿轨最大剖面数，44 表示每个剖面由 44 字节的特征标志组成。有些版本还会出 5km 水平均一的版本，变量名里的数字会变成 2060 之类，读取逻辑相同。拿到文件后先用下面的脚本把变量名和维度全部打印出来，再决定怎么解码。

2.3 44字节掩膜的解码思路：先拆字节，再谈显示

44 字节不是 44 个普通整数值，而是按位打包的特征标志。每个位或几个连续位对应一个分类结果：特征类型、云层类型、气溶胶类型、CAD 分数相关标志、晴空置信度等。直接把这 44 个字节当成 44 个灰度值画图，能得到云层轮廓但拿不到类型信息，而且不同字节的数值范围差异很大，画出来一片噪声。

解码的正确思路是：先把 HDF4 里的数组按uint8读进来，再按位运算提取目标字段。严格来说，不同版本（V3-01、V3-40、V4-10）的位号存在细微调整，所以我在实际项目里不会硬编码一套位表，而是先把掩膜数组里前几个字节的最小值、最大值打印出来，再对照所下载版本的文件说明确认位偏移。这一点后面避坑章节还会展开讲。

3. 读取Calipso-VFM：一套可直接跑的Python读取流程

3.1 安装pyhdf并确认HDF4版本：环境准备

VFM 是 HDF4 格式，这一点是第一个坑。很多新同事拿到.hdf后缀就直接h5py.File(path)，立刻报错Unable to open file (file signature not found)。HDF4 和 HDF5 虽然同名，文件签名完全不同，h5py 只认 HDF5。读取 HDF4 的 Python 库是pyhdf。

安装pyhdf在 Linux 下往往要先装 HDF4 运行库。Debian/Ubuntu 系常见做法是先安装系统包，再pip install pyhdf。Windows 环境下 pip 直接装 wheel 经常缺依赖，我一般建议用 conda 从 conda-forge 安装：

conda install -c conda-forge pyhdf

装完做一次快速验证：能成功导入from pyhdf.SD import SD，说明环境就绪。这个验证很重要，因为pyhdf本身是对底层 HDF4 库的封装，库缺失时导入并不会直接报错，而是在打开文件时才崩。

3.2 最小读取脚本：打开文件、列出SD变量、打印shape

拿到一个 VFM 文件后，我建议第一步永远不要急着读数据，而是先列出文件里有哪些科学数据集（SD），以及每个数据集的维度。这样能避免后面维度顺序搞错导致图横竖颠倒。

from pyhdf.SD import SD, SDC path = "CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf" # 以只读方式打开 HDF4 文件 hdf = SD(path, SDC.READ) # 列出所有科学数据集及维度信息 for name, info in hdf.datasets().items(): print(name, info) # 关闭文件，HDF4 在 pyhdf 里需要显式 end hdf.end()

这段代码输出里能看到每个变量的维度元组和数据类型。我见过很多次Feature_Flight_2048_Column_44_Byte的维度在不同版本里发生变化：有的版本存成[2048, 545, 44]，有的存成[545, 2048, 44]，44 这一维是特征标志字节，545 是高度层数，2048 是沿轨剖面数。如果不先打印维度直接按固定下标取数，轻则画图方向错，重则整段数据都是错位值。

参数说明：SD(path, SDC.READ)的SDC.READ是只读模式，避免误写原始文件；hdf.datasets()返回的是变量名到维度信息字典，不是数据本身，这一步不会把大数据量装进内存，适合快速侦察。

3.3 把Feature_Classifier_Flags解码成云/气溶胶/晴空三类掩膜

读取核心掩膜并解码，我用下面这套流程。代码里先动态把 44 字节那一维挪到最后一维，再按位提取特征类型。特征类型的取值按最常见约定：0 表示无效或未分类，1 表示晴空，2 表示云，3 表示气溶胶。不同版本可能把 2 和 3 互换，但三位分类的思路一致。

import numpy as np from pyhdf.SD import SD, SDC path = "CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf" hdf = SD(path, SDC.READ) # 读取掩膜原始数据 sds = hdf.select("Feature_Flight_2048_Column_44_Byte") raw = sds.get() # numpy 数组 sds.end() # 统一维序为 (沿轨, 高度, 44字节) arr = np.array(raw, dtype=np.uint8) if arr.ndim == 3: # 找到长度为44的那个轴，并移到最后一维 for i, dim in enumerate(arr.shape): if dim == 44: arr = np.moveaxis(arr, i, -1) break # 按位提取特征类型：用第0字节的最低2位 feature_type = arr[..., 0] & 0x03 # 构建三分类掩膜：0=无效, 1=晴空, 2=云, 3=气溶胶 cloud_mask = (feature_type == 2) aerosol_mask = (feature_type == 3) clear_mask = (feature_type == 1) print("原始shape:", raw.shape) print("特征类型取值范围:", np.unique(feature_type)) print("云像元数:", cloud_mask.sum(), "气溶胶像元数:", aerosol_mask.sum()) hdf.end()

这里最值得说明的是arr[..., 0] & 0x03这一步。arr[..., 0]取出每个格点的第一个字节，& 0x03是位与操作，只保留最低两个二进制位，高位的其他特征标志全部置零。为什么用位与而不是直接取整数值？因为 44 字节里几乎每个字节都塞了多个标志，直接看整数值会被无关高位干扰，位与之后才能干净地拿到目标字段。

参数说明：np.moveaxis(arr, i, -1)的作用是把长度为 44 的轴移动到最后一维，这样后面所有下标都用[..., 0]、[..., 1]访问字节位置，代码不依赖具体版本里 44 维在哪个位置，可移植性好。

3.4 数据质量位QA的读取与过滤：避免把无效格点画进图里

VFM 的 44 字节里除了特征类型，还有质量相关标志。我在批量处理时发现，如果不过滤质量标志，有些夜侧数据或信号衰减严重的剖面会出现一整列“假云”。常见做法是同时读取另一个质量字节，把低置信度格点统一改成无效值。

质量位的位置同样随版本变化，我一般先打印前 8 个字节的统计值，观察哪些字节在无效区域恒为 0，再确定质量位偏移。下面的代码展示的是典型的过滤流程，QA_BYTE和QA_BIT需要按你下载的版本说明调整：

# 假设质量标志位于第1字节的最高位，1表示高质量 QA_BYTE = 1 QA_BIT = 7 quality_flag = (arr[..., QA_BYTE] >> QA_BIT) & 0x01 # 低质量格点与无效格点统一置为0 feature_type[quality_flag == 0] = 0 feature_type[(arr[..., 0] & 0x03) == 0] = 0 # 重新统计 print("过滤后云像元数:", (feature_type == 2).sum())

这段代码的核心逻辑是先提质量位，再用布尔索引把低质量格点清零。注意不要直接对raw做过滤，而是生成新的feature_type数组，这样保留原始数据方便回退。参数说明里的两个常量值得单独强调：QA_BYTE和QA_BIT在不同 VFM 版本里有差异，我习惯把它们放在脚本顶部，换版本时只改这两行，不碰主体逻辑。

4. 显示Calipso-VFM垂直剖面：坐标、配色与输出

4.1 用matplotlib的pcolormesh画离散色块

分类掩膜不能像强度图那样用imshow加连续色带渲染，必须用离散颜色映射，每个类别固定一个颜色，否则云、气溶胶、晴空之间的边界会糊成一片渐变。我用matplotlib的pcolormesh和BoundaryNorm来做，因为pcolormesh天然支持非均匀网格，而且对离散分类的色块边界控制比imshow干净。

import matplotlib.pyplot as plt import matplotlib.colors as mcolors import numpy as np # 沿用上一节的 feature_type, 高度网格和沿轨距离 # heights: 每个高度层中心，单位 km # along_dist: 沿轨距离，单位 km # 定义类别与颜色：0-无效, 1-晴空, 2-云, 3-气溶胶 cmap = mcolors.ListedColormap(["#f0f0f0", "#a6d96a", "#1a9641", "#d7191c"]) # 定义离散边界：4个类别需要5个边界值 bounds = [-0.5, 0.5, 1.5, 2.5, 3.5] norm = mcolors.BoundaryNorm(bounds, cmap.N) fig, ax = plt.subplots(figsize=(10, 4)) # X轴用沿轨距离，Y轴用高度 mesh = ax.pcolormesh(along_dist, heights, feature_type.T, cmap=cmap, norm=norm, shading="nearest") ax.set_xlabel("Along-track distance (km)") ax.set_ylabel("Altitude (km)") ax.set_ylim(-0.5, 20) cbar = fig.colorbar(mesh, ax=ax, ticks=[0, 1, 2, 3], extend="neither") cbar.ax.set_yticklabels(["Invalid/Undefined", "Clear air", "Cloud", "Aerosol"]) plt.tight_layout() plt.savefig("vfm_profile.png", dpi=200)

注意feature_type.T这一步：feature_type的 shape 是(沿轨数, 高度数)，而pcolormesh要求传入的二维数组形状是(高度数, 沿轨数)，所以需要转置。这正是前面统一维序的回报——如果开头维度没有理顺，这里转置转错一个轴，整张图就是镜像的。

BoundaryNorm的参数说明：bounds列表把数值区间切成 4 段，[-0.5, 0.5]对应类型 0，[0.5, 1.5]对应类型 1，以此类推。这样即使feature_type里有异常数值，也不会被归一化拉伸成渐变。

4.2 高度轴与沿轨轴的处理：让横轴真正代表地面距离

VFM 的横轴经常被直接画成剖面序号，这在单轨展示时能看，但一旦要对比多轨或叠加地理信息，序号轴完全不够用。我更习惯把横轴换算成沿轨距离，再用副刻度标出经纬度，这样读者一眼能看出“这段云层大概在哪个区域”。

# 读取经纬度 lat = hdf.select("Latitude").get() lon = hdf.select("Longitude").get() # 沿轨距离用球面距离累积计算 from numpy import deg2rad, sin, cos, arcsin, sqrt def haversine(lat1, lon1, lat2, lon2): lat1, lon1, lat2, lon2 = map(deg2rad, [lat1, lon1, lat2, lon2]) dlat = lat2 - lat1 dlon = lon2 - lon1 a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 return 6371.0 * 2 * arcsin(sqrt(a)) # 计算相邻剖面间的距离并累加 deltas = [0.0] for i in range(len(lat)-1): deltas.append(haversine(lat[i], lon[i], lat[i+1], lon[i+1])) along_dist = np.cumsum(deltas)

这段代码用 Haversine 公式逐段计算剖面间的地面距离并累计，横轴单位是千米。为什么不用np.arange(n_profile)？因为卫星在极区附近的星下点轨迹疏密不同，剖面序号和实际地面距离不是线性关系，尤其在轨道最高纬度的拐弯段，序号轴会明显拉伸或压缩距离。用真实距离画出来的云层横向尺度才可信。

高度轴的处理相对简单，直接用文件里的Lidar_Data_Altitudes就行，但要注意这个高度是相对海平面的几何高度，单位是千米，低层可能从负值开始，画图时通常只显示-0.5到20km区间，因为 20km 以上只有极其罕见的极地平流层云，留太多空白反而削弱图面信息。

4.3 配色与图例的工程选择：可读性优先，不要用彩虹色

VFM 显示这块，我见过最典型的翻车是直接用viridis或jet渲染掩膜，结果云、气溶胶、晴空三类颜色在色带上距离很近，打印成黑白稿后完全分不清。离散分类图的核心是类别间的视觉区分度，不是美观。

我的配色习惯是：无效用浅灰，晴空用浅绿，云用深绿，气溶胶用红色。这套方案的逻辑是，云和气溶胶是研究者最关注的两类，深绿和红在视觉上差异最大；晴空作为背景用浅色，不会抢前景注意力。如果你需要打印成灰度图，建议把云的灰度值压暗、气溶胶提亮，确保印刷后仍可分。

图例不要用colorbar的默认数值刻度，直接改成类别名。上面代码里ticks=[0,1,2,3]配set_yticklabels的做法，比在图例右上角放一堆数字直观得多，也方便直接贴进报告。

5. Calipso-VFM读取与显示避坑：5条排查记录

5.1 现象：h5py打开HDF4文件直接报错，文件签名找不到

原因：VFM 是 HDF4 格式，h5py 只支持 HDF5，两者文件头签名完全不同。这不是文件损坏，是库选错了。

解决：换pyhdf.SD打开。如果项目里已经大量使用 h5py 的读取流程，可以把 VFM 先转成 HDF5 或 NetCDF4 再走原流程，但我不推荐，因为转格式会引入属性丢失和维度重排的额外风险。直接在读取函数里加一个 HDF4 分支更干净。

5.2 现象：读出来的掩膜数组维度是(545, 2048, 44)，按(2048, 545, 44)处理导致剖面方向错乱

原因：VFM 的维度排列在不同处理版本里不统一，代码写死了下标位置，换一个文件就出错。

解决：统一维序后再操作。具体做法是扫描数组各维长度，找到等于 44 的那一维挪到末尾，再用[..., 0]方式访问字节。这段逻辑我已经放在第 3 章的代码里，每次处理新版本文件时先运行维度检查脚本，确认后再解码。

5.3 现象：特征类型提取出来全是0或全是255，画出来一片灰

原因：位偏移用错。VFM 的 44 字节在不同版本里，特征类型位所在字节号和位号有调整。V3-01 和 V4-10 的位表就有差异，如果直接照抄网上的旧代码，很可能取到的是填充字节或质量标志字节。

解决：先打印掩膜数组每个字节的统计值，找到在云区有明显分段特征的字节，再用位与提取。我在实际项目里会在脚本顶部维护一个版本到位号的对照字典，换数据版本时改一行配置即可。

5.4 现象：云层剖面画出来是一条紧贴横轴的扁平细线，高度结构完全看不出来

原因：pcolormesh默认的纵横比不是按数据单位等比缩放，而显示窗口宽度远大于高度，导致垂直方向被严重压缩。

解决：画出图后先看一眼云层的厚度占比，再手动设置ax.set_aspect或用figsize控制画布比例。我一般把figsize设为(12, 4)，并将高度轴范围限制在-0.5到20km，这样云层的垂直结构能看到，又不会因为显示到 30km 高空而浪费大量画布。

5.5 现象：跨UTC日期的轨段拼接时，相邻两轨之间出现一条斜向的“假云墙”

原因：文件命名里用的是 UTC 日期，当地日期跨 UTC 零点时，同一天的轨道文件可能分布在两个 UTC 日里，按文件名排序会漏掉或错序；另外 VFM 每个文件是半个轨道段，相邻文件之间本身有重叠或间隔，画拼接图时没按时间戳排序就会错位。

解决：用Profile_UTC_Time作为排序键，不依赖文件名日期。拼接前先打印每个文件的起始和结束时间，确认顺序连续后再合并；拼接图的公共横轴用累计的沿轨距离或直接由时间戳换算。

6. 进阶技巧：一天多轨段批量出图，并输出云/气溶胶占比统计

6.1 批处理循环：按日期筛选轨道段

单轨图画得再好，研究里真正需要的是连续一个月甚至一年的统计结果。我处理整日数据时，先按轨道段文件名聚合同一 UTC 日期的文件，遍历每个文件解码，再把同一天的剖面拼到同一个数组里。

import glob file_list = sorted(glob.glob("CAL_LID_L2_VFM-*.hdf")) daily_masks = {} for path in file_list: # 解析UTC日期，从文件名第16到25位 date_str = path.split(".")[2] # 解码掩膜，复用前面封装的函数 feature_type = decode_vfm_feature_type(path) if date_str not in daily_masks: daily_masks[date_str] = [] daily_masks[date_str].append(feature_type)

这段代码的核心是按文件名里的 UTC 日期做分组，把一天内的多个轨段收集到列表里。注意这里排序必须用文件名里的 UTC 字段，不能用本地时间，否则跨 UTC 零点那天会少一段数据。分组之后，每个日期对应的列表里每个元素是一轨的二维掩膜，下一步就可以逐轨画图或做统计。

6.2 统计掩膜占比并输出CSV

批处理不能只出图，还要把每一轨的云占比、气溶胶占比统计成表格，这样后续做纬度带分析时不用重新读数据。

import csv with open("vfm_daily_stats.csv", "w", newline="") as f: writer = csv.writer(f) writer.writerow(["date", "profile_count", "cloud_frac", "aerosol_frac"]) for date_str, masks in daily_masks.items(): # 把同一天所有轨段拼接起来 all_masks = np.concatenate(masks, axis=0) total = all_masks.size cloud_frac = (all_masks == 2).sum() / total aerosol_frac = (all_masks == 3).sum() / total writer.writerow([date_str, all_masks.shape[0], round(cloud_frac, 4), round(aerosol_frac, 4)])

统计时要注意分母是全部格点数还是有效格点数。我一般用全部格点数，同时额外输出一列valid_frac，让下游分析能判断统计稳健性——如果某天有效格点占比只有 30%，那云占比再高也不具代表性。

参数说明：np.concatenate(masks, axis=0)要求所有轨段的feature_type形状一致，包括高度层数。不同版本 VFM 的高度层数可能不同，拼之前先统一，否则会直接报维度不匹配。

6.3 验证方法：把统计结果和原始后向散射剖面对照

批量出完图后，我习惯做一道验证：挑一轨数据，把 VFM 掩膜叠到同轨的 532nm 衰减后向散射剖面上，检查云层和气溶胶层的位置是否与强度高值区吻合。做法是画双面板，上面放 L1 强度，下面放 VFM 掩膜，两图横轴和高度轴完全对齐，肉眼核对边界偏移。这一步能快速暴露位号错误、维度错位这类代码问题，比直接看统计数字可靠得多。

另外一个更廉价的验证是把 VFM 气溶胶出现频率按纬度带统计，画成纬度-高度剖面图。真实数据里，沙尘层通常出现在北非到中亚的纬度带，海洋气溶胶在低纬海洋上空更常见，烟羽则常出现在火点纬度附近。如果你的统计图里气溶胶均匀分布在所有纬度，大概率是解码位号有问题，而不是真实大气如此。

这套批处理流程跑过一轮之后，我自己的习惯是保留每个版本文件的位号配置、解码脚本和输出图放在同一个目录下，换版本时先对照位号再跑批处理，避免辛辛苦苦处理完几万条剖面后发现位号用错，全部返工。也希望这篇文章能帮你少走这几趟弯路。

本文还有配套的精品资源，点击获取