简介面向雷达基数据解析与显示的VC工程源码包适合雷达信号处理、气象监测或交通监控方向的软件开发者用于解决基数据从二进制解码、算法处理到PPI雷达图像可视化的完整流程。压缩包内共80个文件大小约7.83MB其中包含17个头文件与13个C源文件另有工程配置、图标位图资源、静态库和中间编译文件构成一套可直接编译学习的MFC项目。目前已有344人浏览学习。工程提供了RadarYHM主框架、PPI显示视图、CRadarBase基数据解析、图像类型转换、JPEG编码库等模块能够帮助读者理解雷达基数据的存储格式、脉冲压缩与频率分集等信号处理思路掌握在文档视图架构下将原始数据转换为二维/三维雷达图的编码方式。源码目录结构清晰适合作为二次开发的基础便于将相关解析与显示类模块抽取复用搭建自己的雷达数据解析工具。1. Radaryhm 在做什么一套能直接用的雷达基数据解析与显示方案做气象或雷达信号处理的人大概率见过这样的尴尬场面手里拿到一卷基数据文件扩展名五花八门打开全是二进制乱码没有配套软件就完全没法看。RadarYHM 这个名字对应的就是一套「把雷达基数据读出来、画成图」的解析显示工具常见于气象台站、雷达厂家交付目录和高校课题组的内部资料包里通常以一个 rar 压缩包形式分发解压后就是源码加样例基数据。它解决的问题很直接让你不依赖厂家闭源软件也能把反射率、速度、谱宽从原始基数据里解析出来并画成常见的 PPI 显示图甚至输出成标准化数据给后续算法使用。适合做雷达数据解析入门、历史基数据回放、以及需要批处理多文件的从业者。这篇文章就把这套方案的原理、关键代码和踩坑点完整展开照着做就能把属于自己的雷达显示工具跑起来。2. 基数据文件结构读数据前必须先搞懂的三个字段2.1 基数据到底存了什么反射率、速度、谱宽在文件里的组织方式天气雷达基数据不是一张图而是一根根径向扫描的原始观测记录。以国内最常见的 CINRAD SA/SB 雷达为例一个完整的体扫由多个仰角扫描组成每个仰角扫描包含若干根径向每根径向在同一个距离库上记录反射率、径向速度、谱宽等观测值。RadarYHM 这类工具的核心第一步就是把文件里这些径向按顺序解析出来。理解基数据结构先要认识三个关键字段站点与雷达参数区包含雷达站名、经纬度、天线高度、体扫模式、仰角数、每层仰角值、脉冲重复频率等。径向头区每根径向开始处有状态标识、方位角、仰角、径向距离库数、数据格式类型等。数据区按距离库排列的观测值每个库通常对应 1 公里或 0.5 公里距离数值经过缩放和偏移编码需要用公式还原成物理量。很多新手直接跳过文件头去读数据结果画出来的图方位角错乱、距离标尺不对根源就是没搞清这三个区的位置和长度。2.2 二进制读取的骨架通用基数据解析最小实现不同型号雷达的基数据格式略有差异但解析思路完全一致按字节读、按结构拆、按公式还原。下面给一个通用的二进制基数据解析骨架以 SA 格式常用的「磁带头 径向块」布局为例。import struct import numpy as np def parse_base_data(file_path): # 以二进制方式打开基数据文件 with open(file_path, rb) as f: # 1. 读取固定长度的磁带头不同雷达略有差异 header f.read(128) # 站点名通常是 ASCII 字符串前 8 字节 station_name header[0:8].decode(ascii, errorsignore).strip() # 雷达经度、纬度、天线高度常见为 double注意字节序 lon, lat, height struct.unpack(ddd, header[24:48]) # 2. 读取体扫参数仰角层数、每层径向数 scan_count struct.unpack(H, f.read(2))[0] # 3. 循环读取每一根径向 radial_data [] for scan_idx in range(scan_count): # 径向头一般固定 32 字节 radial_header f.read(32) azimuth struct.unpack(H, radial_header[0:2])[0] / 64.0 elevation struct.unpack(H, radial_header[2:4])[0] / 64.0 bin_count struct.unpack(H, radial_header[4:6])[0] # 读取该径向的反射率数据每个库 1 字节或 2 字节 raw_ref np.frombuffer(f.read(bin_count), dtypenp.uint8) # 缩放与偏移实际部分格式用 raw * 0.5 - 32.0 得到 dBZ ref raw_ref.astype(np.float32) * 0.5 - 32.0 radial_data.append({ azimuth: azimuth, elevation: elevation, reflectivity: ref, }) return { station: station_name, lon: lon, lat: lat, height: height, radials: radial_data, }这段代码的逻辑重点在两处。第一处是struct.unpack(ddd, ...)的前缀雷达基数据文件里的大端小端问题极其容易翻车SA 格式不少数据块是大端序但个别字段又按小端处理必须逐字段核对厂家格式说明。第二处是缩放还原公式基数据文件里存的几乎不会直接是 dBZ 浮点数而是用一个字节或两个字节压缩存储的整数RadarYHM 这类工具在还原时一般还有无效值标记比如数值为 0 或一个极值表示缺测解析时要把这些值剔除否则图上会出现一整圈异常的强回波点。2.3 坐标系参数从极坐标到经纬度网格的转换解析出来的反射率数据是极坐标系下的每根径向对应一个固定的方位角和仰角数据点沿径向等距分布。要显示成地图上常见的 PPI 图就要把极坐标投影到笛卡尔坐标或经纬度坐标。这一步的换算公式不复杂但很多人在这个环节把方位角搞错。PPI 显示时常用「雷达极坐标展开」的方式横轴是方位角纵轴是距离。但更常用的叠加地图方式需要把每个库的坐标换算到雷达站点平面甚至经纬度。def polar_to_cartesian(radial_azimuth, gate_range_km, site_lon, site_lat): # 方位角转弧度注意气象雷达方位角通常以正北为 0顺时针增加 az_rad np.deg2rad(radial_azimuth) # 库距离假设每个库 1 km x_km gate_range_km * np.sin(az_rad) y_km gate_range_km * np.cos(az_rad) # 粗略经纬度换算中纬度地区 1 度经度约 111 km * cos(lat) dlon x_km / (111.0 * np.cos(np.deg2rad(site_lat))) dlat y_km / 111.0 return site_lon dlon, site_lat dlat注意这里的sin和cos对应关系。北为 0 度顺时针正东是 90 度那么 x 方向东向用 siny 方向北向用 cos。如果先把方位角当成数学极坐标的逆时针角度画出来的图会左右翻转这在雷达显示里是特别低级的错误。更精确的做法是考虑地球曲率与雷达波束折射做波束高度订正但对 PPI 平面显示来说上述转换已经能定位准确。3. 用 Python 实现解析与显示从文件头到 PPI 图3.1 准备环境与目录结构解析基数据我一般建议用 Python配合 numpy、struct、matplotlib 就能完成全流程不需要引入重型的 GIS 框架。创建一个工作目录结构保持简单radar_workspace/ |-- data/ | -- 20240815_0001.bin |-- parse_radar.py |-- plot_ppi.py |-- output/依赖安装只做三件事numpy负责数组运算struct是 Python 内置库负责二进制读取matplotlib负责画图。如果后续要叠加地理底图再引入cartopy但基础显示不依赖它。RadarYHM 这类工具包解压后一般也是这样的结构核心解析脚本加一两个样例基数据。拿到手先别急着跑用十六进制编辑器打开样例数据看一眼头部确认字段起始偏移和文件字节序这比盲目跑代码省时间。3.2 逐径向读出扫描数据上一章的骨架代码已经完成了逐径向读取。现在要做的是把多个仰角的数据组织成「层」的结构因为一个体扫文件里往往包含多个仰角扫描每个仰角扫描的径向数和起始径向索引都不同。def group_by_elevation(radials): 将径向按仰角分组返回 {仰角: [径向列表]} 部分雷达基数据的仰角值有微小波动要做容差匹配 groups {} for r in radials: # 仰角值四舍五入到 0.1 度避免浮点波动导致的误分组 elev_key round(r[elevation], 1) if elev_key not in groups: groups[elev_key] [] groups[elev_key].append(r) return groups这里的容差处理很关键。基数据文件里记录的仰角不是精确等于名义仰角比如 1.5 度仰角实际可能记录成 1.49 或 1.51。如果直接用原始浮点值做字典 key每一根径向都可能被分到不同组后续生成极坐标网格时会出现径向缺失。另外有些体扫文件是扫描模式交替的比如 VCP21 和 VCP31 的仰角层数不同分组前先打印所有径向的仰角分布确认实际层数。3.3 绘制反射率 PPI一张图看清降水结构极坐标径向数据要先转成扇形网格图。这里可以直接用 matplotlib 的pcolormesh把方位角和距离构造成极坐标网格。import matplotlib.pyplot as plt import numpy as np def plot_reflectivity_ppi(radials, max_range_km150, output_pathoutput/ppi.png): # 取 0.5 度仰角那层径向 azimuths np.array([r[azimuth] for r in radials]) refs np.array([r[reflectivity] for r in radials]) # 距离库轴假设每个库 1 km库数取径向数据长度 gate_count refs.shape[1] ranges_km np.arange(gate_count) # 构造极坐标网格方位角沿 y 方向扩展距离沿 x 方向扩展 az_grid, range_grid np.meshgrid(np.deg2rad(azimuths), ranges_km, indexingij) # 使用极坐标投影画图 fig plt.figure(figsize(10, 10)) ax fig.add_subplot(111, projectionpolar) # 注意方位角顺序基数据里方位角通常升序排列直接画即可 mesh ax.pcolormesh(az_grid, range_grid, refs, cmapjet, shadingauto, vmin0, vmax65) ax.set_theta_zero_location(N) ax.set_theta_direction(-1) ax.set_rlim(0, max_range_km) plt.colorbar(mesh, axax, labeldBZ) plt.title(Reflectivity PPI 0.5 deg) plt.tight_layout() plt.savefig(output_path, dpi200)set_theta_zero_location(N)和set_theta_direction(-1)是画天气雷达 PPI 图最容易忽略的两行。默认的 matplotlib 极坐标图 0 度在正右方且角度逆时针增加但雷达显示习惯是正北为 0、顺时针为正这两行不写整个回波图会旋转 90 度且镜像你看到的降水回波会出现在完全错误的地物方位上。vmin0, vmax65是反射率图常用的色标范围低于 0 dBZ 基本是噪声或弱杂波超过 65 dBZ 属于极端强回波在显示时截断掉可以提升图面对比度。如果做科研写论文建议用pyart_graph或wradlib的标准色标与国内业务显示更接近。3.4 速度图与谱宽显示时要注意的量纲速度数据解析出来后物理量纲是米每秒通常范围是 [-32, 32] 或 [-64, 64] 之间取决于脉冲重复频率。雷达显示速度图时会用红绿两套色标正速度表示远离雷达负速度表示朝向雷达运动。速度显示时最需要注意的是「速度模糊」问题。def plot_velocity_ppi(radials, output_pathoutput/velocity_ppi.png): azimuths np.array([r[azimuth] for r in radials]) vels np.array([r[velocity] for r in radials]) # 速度值范围校验超出正常范围的通常是解模糊失败或无效值 print(velocity range:, np.nanmin(vels), np.nanmax(vels)) fig plt.figure(figsize(10, 10)) ax fig.add_subplot(111, projectionpolar) mesh ax.pcolormesh(np.deg2rad(azimuths), np.arange(vels.shape[1]), vels, cmapBrBG, vmin-32, vmax32) ax.set_theta_zero_location(N) ax.set_theta_direction(-1) plt.colorbar(mesh, axax, labelm/s) plt.savefig(output_path, dpi200)速度图的vmin和vmax要与你解析时使用的最大不模糊速度匹配。如果数据里速度值已经超过设定的显示范围pcolormesh会把超范围的部分截断成色标两端颜色造成视觉上的大范围正负速度区。RadarYHM 这类工具通常会在界面上显示最大不模糊速度参数解析代码里也要把这个参数打印出来供判断。谱宽图的量纲同样是米每秒一般显示范围 0 到 10 m/s超过 10 的谱宽往往是地物杂波或信号质量差导致的。4. 处理多 PRF、径向缺失与数据订正三个必调参数4.1 多脉冲重复频率解模糊现代天气雷达普遍使用双 PRF 或多 PRF 方式扩展最大不模糊速度。基数据文件里会记录每个径向使用的 PRF 组合解析速度数据时如果不做解模糊处理速度图会出现一块块相邻径向速度突变超过最大不模糊速度的「马赛克」。解模糊的常见做法是速度退模糊。径向数据中相邻距离库的速度差超过阈值时加减若干个整数倍的nyquist_v让速度剖面连续。RadarYHM 的处理逻辑一般在这一步给用户提供一个可选开关默认开启因为在业务显示中未解模糊的速度图几乎没有应用价值。def dealias_velocity(velocity, nyquist_v, gate_spacing_m1000): 一维速度退模糊沿距离库方向做连续性订正 dealiased velocity.copy() # 找出有效数据区 valid ~np.isnan(dealiased) if valid.sum() 0: return dealiased diff_threshold nyquist_v * 0.6 # 经验阈值超过即认为发生跳变 for gate in range(1, len(dealiased)): if not valid[gate] or not valid[gate - 1]: continue diff dealiased[gate] - dealiased[gate - 1] n_nyq round(-diff / (2 * nyquist_v)) if abs(diff n_nyq * 2 * nyquist_v) abs(diff): dealiased[gate] n_nyq * 2 * nyquist_v return dealiased这个一维退模糊算法是基础版对大多数体扫数据够用。真正到强对流天气时径向间风速变化剧烈一维连续性假设经常失效退化出来会出现「条状模糊带」。那时候先用一维退模糊再做径向中值滤波把散点状的异常值抹掉可以改善显示效果。这个算法的阈值0.6 * nyquist_v是个经验值PRF 越高阈值越小要自己根据雷达型号调整。4.2 地物杂波过滤解析显示的地物回波通常集中在雷达站周边 20 公里内特征是反射率强、谱宽小、且不随体扫变化。如果不做过滤直接显示PPI 图中心会出现一片红色的强回波把降水信号完全盖住。最常用的过滤思路很简单比较同一位置不同仰角的数据。地物回波在低仰角明显、高仰角迅速消失而降水回波随仰角变化平缓。def filter_clutter(ppi_low, ppi_high, clutter_threshold5.0): ppi_low: 0.5 度仰角反射率网格 ppi_high: 1.5 度仰角反射率网格 当低仰角明显强于高仰角时判为地物并置 NaN filtered ppi_low.copy() diff ppi_low - ppi_high clutter_mask diff clutter_threshold filtered[clutter_mask] np.nan return filtered这个办法在层状云降水时效果好但在强对流时 1.5 度仰角可能已经穿过云体顶回波强度大幅衰减会把真实的强降水误删。所以我在实际使用中会把clutter_threshold调大或者只在显示层做修饰不把过滤结果写回基数据。4.3 显示参数设置RadarYHM 的显示界面里通常有几个高频调整的参数搞明白它们的含义能让出图质量提升一个档次插值方法。极坐标数据转到规则网格时nearest和linear是最常用的两种。nearest速度快但显示有明显锯齿linear平滑但会把小尺度回波抹淡。我建议显示用linear定量分析用nearest因为插值永远会改变数据分布。\n\n最大显示距离。国内多普勒天气雷达体扫最大探测距离一般 230 公里或 460 公里但实际有效反射率距离通常只有 150 到 230 公里。把最大显示距离设置到 230 公里可以在不牺牲中心区域分辨率的情况下显示完整回波。部分基数据文件里距离库数是 460 或 920对应每个库 0.5 公里或 1 公里解析时把gate_length读出来并参与网格构建这是 RadarYHM 这类工具内部会自动处理但用户容易忽视的参数。颜色映射。反射率用NWSReflectivity之类的色标速度用蓝绿红渐变谱宽用灰度。用 matplotlib 默认的jet也能看但在强回波区容易出现色带断层。5. 避坑清单基数据解析显示的常见问题与排查5.1 文件字节序错误读出数据全是噪点现象解析出来的反射率数据在图上显示为均匀的颗粒噪声没有明显的回波区结构数值范围也极其离谱比如反射率全部是 100 以上。原因雷达基数据文件很多按大端序存储而 Python 的struct默认使用小端序。头部字段没指定前缀读取的短整型、浮点数全部错位导致径向头里的角度、库数、数据长度全部变成无意义的大数。解决先用十六进制工具打开文件头部找一个已知字段验证字节序。比如站点名后面的仰角值如果是大端序正常应看到类似0x0001这样的值如果看到0x0100那就是字节序反了。所有struct.unpack调用统一加上前缀并且在解析完成后打印一条径向的头字段人工核对方位角和仰角值是否在合理范围。5.2 雷达径向角度跳变或缺失导致 PPI 图出现缺口现象画出来的 PPI 图上有一个扇形的数据缺口或者回波结构出现沿方位角方向的「撕裂」。原因基数据文件在扫描过程中可能丢失径向尤其是强雷暴导致雷达在某个方位角上数据质量标记异常。如果代码里假设径向数是固定值并按这个值切分文件后面的解析全部错位图上表现为一片混乱的色块。另一种情况是方位角从小到大出现了跳变比如从 350 度直接跳到 2 度而网格构造时没有做环状首尾相接的处理。解决解析时不要按固定径向数读取。先读取每个径向头的径向状态字节跳过标记为无效的径向。构造极坐标网格时改用径向的真实方位角值填充不要假设方位角均匀分布。如果缺口恰好跨过 0 度就把方位角大于 350 度和小于 10 度的数据在网格里拼接起来否则 matplotlib 会在缺口处做白色插值。5.3 反射率出现大片负值和空心回波现象图上中心区域出现一个负 dBZ 的空洞周围回波明显有时候整个径向都显示为负值看起来像回波被挖掉了一圈。原因反射率数据的缩放公式里通常有一个偏移量比如raw * 0.5 - 32如果代码里偏移量写错比如只做了缩放没做偏移数据整体会偏小 32 dBZ大量弱回波变成负值。另一个常见原因是距离库起始位置没对齐雷达探测最近几个距离库没有数据文件里用特殊值标记解析时没把这些特殊值剔除。解决验证缩放公式时直接用样例数据里一个已知强回波点手工计算一遍对比代码输出。对无效值解析时检查原始编码是否是 0 或大于等于某个阈值直接替换成np.nan。绘图时用np.nanmin检查数据分布确认负值比例如果负值超过 30%先怀疑缩放公式而非滤波器。5.4 速度图出现大范围红蓝交替条纹现象速度图上有沿径向排列的、规律性红蓝交替的条纹且条纹间距大致相等看起来像斑马纹。原因这是典型的未解模糊现象。数据里真实速度超过了最大不模糊速度折叠后形成锯齿直接显示就成了条纹。也可能是双 PRF 处理时没匹配每个径向的 PRF 值导致多个模糊速度叠加。解决先打印文件头里的 PRF 参数确认真实的最大不模糊速度。然后用上一章的退模糊算法处理。退模糊后如果条纹还在检查速度数据在距离库方向是否做了差分解算——部分基数据存储的是速度差而不是绝对速度需要累加还原。5.5 叠加地图时回波位置偏移几十公里现象回波中心相对地图上的河流、城市标记有明显偏移所有回波都朝同一个方向平移。原因经纬度网格转换时没有把雷达站点经纬度加进去或者点距换算用了错误的每度公里数。部分雷达基数据站点的经纬度是度分格式比如113.32实际代表 113 度 19 分直接当十进制度数用会偏移约 0.1 度中纬度对应约 10 公里。解决统一把经纬度转成十进制度数用站点经纬度做基准。验证方法很简单把解析出的第一根径向方位角和距离转换为经纬度后在通用地图上打点看是否落在站点周边。如果所有点整体南移或东移检查是不是把纬度的 cos 项丢了。6. 一个可用的快速显示脚本验证你的解析结果6.1 最小可运行代码把前面的解析、分组、显示整合成一个独立脚本作为验证解析结果的标准工具。这个脚本只做一件事输入一个基数据文件输出一张低仰角反射率 PPI 图和一个数据范围摘要。import numpy as np import matplotlib.pyplot as plt import struct def quick_look(file_path, elevation0.5, max_range_km150): # 1. 读取并解析全部径向 parsed parse_base_data(file_path) groups group_by_elevation(parsed[radials]) # 2. 选择目标仰角最近的一组径向 elev_keys sorted(groups.keys()) target_key min(elev_keys, keylambda e: abs(e - elevation)) radials groups[target_key] print(target elevation:, target_key) # 3. 构建极坐标数组 azimuths np.array([r[azimuth] for r in radials]) refs np.array([r[reflectivity] for r in radials]) valid_mask ~np.isnan(refs) print(valid ratio:, valid_mask.sum() / refs.size) # 4. 绘制 PPI az_grid, range_grid np.meshgrid(np.deg2rad(azimuths), np.arange(refs.shape[1]), indexingij) fig plt.figure(figsize(9, 9)) ax fig.add_subplot(111, projectionpolar) ax.pcolormesh(az_grid, range_grid, refs, cmapjet, vmin0, vmax65) ax.set_theta_zero_location(N) ax.set_theta_direction(-1) ax.set_rlim(0, max_range_km) plt.title(f{parsed[station]} {target_key} deg PPI) plt.savefig(output/quick_look.png, dpi200) print(saved to output/quick_look.png)脚本运行后先看数据有效比例。正常体扫数据的有效比例应该在 70% 以上如果低于 50%大概率是无效值没剔除或解析错位。再看回波形态层状云降水应当是片状结构对流降水应当是块状结构边界清晰。如果图上出现完整的圆环状条纹检查是不是把某根径向的数据重复填充到了所有方位角上。6.2 如何对比验证解析结果拿到一个不熟悉的基数据文件时我建议做两步交叉验证。第一步是用通用的第三方库解析同一份文件常见的选择有 wradlib、Py-ART、CINRAD 扩展库读取同一仰角的反射率并导出为一个 npy 文件与自己的解析结果做数值对比。允许的差异在 0.01 dBZ 级别如果偏差明显检查缩放公式和无效值阈值。第二步是姿态验证。找到基数据对应的业务产品图比如雷达网站发布的反射率拼图对比回波中心位置、回波强度分布。这一步虽然不精确但能快速发现坐标转换错误和方位角方向错误。有一次我发现解析结果整体旋转了 90 度就是靠与官方 PPI 图对比发现的。6.3 扩展方向RadarYHM 这个方向做完解析和显示后续可以自然延伸出几条实用链。把基数据解析模块抽出来做批处理用于历史数据的统计和个例库建设把极坐标数据插值成三维网格后可以用 matplotlib 画单仰角 PPI 叠加地形阴影也可以输出成 NetCDF 格式方便与数值模式数据做融合。显示层可以考虑用 Cartopy 做底图叠加这一步能直接提升输出图的可用性尤其适合做汇报材料和论文插图。我自己的习惯是给 quick_look 脚本加一个命令行参数支持一次传入多份基数据文件自动批量出图。但如果你要处理的是实时业务数据流就需要注意内存释放问题因为一个完整体扫的反射率、速度、谱宽三个数据场同时驻留内存时占用可能上百兆连续处理多个文件时要用生成器逐文件解析而不是一次性读入列表。\n\n这套解析显示方案做到这里已经覆盖了基数据从二进制流到可视化 PPI 图的完整链条。每一个环节的代码都不复杂但字节序、角度方向、无效值过滤、退模糊这些细节决定了最终效果是专业雷达图还是「噪声拼贴画」。我早期在这上面翻过最大的车就是拿到一个新格式的基数据不验证字节序就直接跑输出图完全没法看还以为是滤波算法问题排查了两天才发现只是前缀写漏了。希望这些血泪经验能帮你少走这些弯路。这些细节验证之后你手头这套工具就真正可用、可扩展了希望帮到你。本文还有配套的精品资源点击获取