简介这是一份面向雷达数据处理与可视化开发者的C工程资源包聚焦基数据解析与雷达显示两大环节。压缩包内共80个文件以头文件和源文件为主涵盖VC工程配置、可执行文件及图标位图等资源整体仅7.83MB。该资源已有344人浏览学习。其中包含PPI显示模块、CRadarBase基类封装、图像类型转换与JPEG处理等核心代码覆盖从原始二进制数据读取、解码到图像显示的实现路径可帮助理解脉冲压缩、频率分集等雷达信号处理技术在工程中的实际落地方式。PPI显示模块与图像格式转换部分对二维雷达成像和GUI更新有直接参考价值。同时完整的MFC工程架构也为读者提供了雷达显示界面与数据管理模块的参考范式适合具备一定C与C编程和数字信号处理基础、希望快速上手雷达数据解析或借鉴完整项目框架的开发者。1. RadarYHM拿到一个雷达基数据压缩包之后你要做的第一件事做过气象或者空管雷达数据处理的人大概率都有过这种经历甲方丢来一个名为“RadarYHM.rar”的压缩包里面藏着一坨看起来像乱码的二进制文件文件名写着“基数据”。你说它是雷达基数据没错但你不知道它的仰角是怎么排的、方位分辨率是多少、强度场和速度场的存储顺序是什么更别说怎么把它变成一张能看的雷达图。RadarYHM.rar 代表的就是这个场景雷达基数据解析程序包负责把厂家原始格式的基数据读出来、换算成物理量再交给显示模块做 PPI 或 CAPPI 渲染。这篇笔记面向的是要用现成程序或自己写脚本把基数据跑通的工作人员目标是让你从拿到一个只有文件名的基数据包到看到第一张完整雷达图全程不靠猜。2. 基数据到底长什么样先搞懂极坐标存储再聊解析很多人第一次打开基数据文件就懵了因为里面全是十六进制字节根本看不到“距离”“方位”“反射率”这些字眼。其实雷达基数据的存储思路几十年没变过雷达绕垂直轴一圈一圈地扫描每扫一圈沿一条径向线按距离库range bin采样把回波强度记录下来。所以整包数据本质上是一个极坐标三维数组——方位角、距离库、仰角三个维度定一个采样点。2.1 径向数据块结构方位头与距离库的读取方式不同厂家的基数据格式差异很大但只要拿到一份格式文档解析逻辑几乎都是同一个套路文件先是一个文件头描述雷达站点经纬度、雷达型号、体扫模式接着就是很多个“径向数据块”每一块以一个径向头开头径向头里写着本径仰角、方位角、径向序号、距离库数量、库长等元信息径向头后面跟着三个等长的数据段分别对应反射率、速度、谱宽。常见做法是在解析前先读一段诊断信息确认为大端还是小端字节序。国内气象雷达常见的基数据文件是 SA/SB 波段标准格式还有不少厂家在标准格式基础上做了“私有扩展”所以程序里通常要留一个“格式探测开关”。用 Python 写最直接的读取逻辑时我一般先用struct.unpack按文档把径向头拆出来。import struct # 假设已按文档确认径向头 144 字节强度段在前速度段在中 def read_radial_head(buf, offset): # 这里以某常见格式为例各厂家的字段偏移以实际文档为准 fields struct.unpack_from(HHI, buf, offset) # 仰角*64、方位角*64、径向号 elev_raw, azim_raw, radial_no [x for x in fields] # 第 3 个位置是库长单位米第 4 个位置是距离库数 range_len, bin_count struct.unpack_from(HH, buf, offset 6) return { elevation: elev_raw / 64.0, azimuth: azim_raw / 64.0, range_len_m: range_len, bin_count: bin_count } # 调用时开一个二进制文件按步长跳着读 with open(radar_base_data.bin, rb) as f: raw f.read() offset 0 radials [] while offset len(raw): head read_radial_head(raw, offset) # 径向头长度不同数据段长度 库数 * 2short 存储 data_len head[bin_count] * 2 offset 144 # 跳到强度段 # 强度、速度、谱宽依次排列这里一次性把三段读走 ref struct.unpack_from(%dH % head[bin_count], raw, offset) vel struct.unpack_from(%dH % head[bin_count], raw, offset data_len) wid struct.unpack_from(%dH % head[bin_count], raw, offset data_len * 2) offset 144 data_len * 3先说明一个关键点elevation和azimuth经常不是直接存浮点数而是存一个整数读出来除以一个比例因子上面的 64 就是最常见的档位有些雷达用 256。如果你读了数据之后发现角度全部偏大或偏小几十倍先怀疑这个比例因子而不是怀疑字节序。另一个常见情况是私有格式里每个距离库不是一个 short 而是两个 short前一个存原始计数后一个存质量标记这种情况可以把上面代码的144改成64或其他实际头长再把data_len乘以 2、4 或者 6。2.2 字节序、仰角扫描顺序与雷达常数的影响字节序判断不靠猜。一个很土但有效的办法是读一个小文件找到文件头里写的站点经纬度如果解析出的纬度在 30 到 55 之间、经度在 70 到 140 之间那字节序基本就对了如果读了能算出纬度 11000 多这种数字那就是断反了。注意一下这里说的基数据解析核心是“把字节翻译成场”而不是“计算雷达气象方程”。雷达常数比如发射功率、天线增益影响的是从原始计数值换算到 dBZ 的斜率反映在标定参数里解析程序通常只负责把已标定或者半标定的计数读出来真正换算成物理量要再看格式文档给的公式。仰角的组织方式也需要注意。雷达一次体扫一般包含 9 到 14 个仰角从低到高依次扫描。文件里的径向数据不一定按仰角分组排列有的厂家按时间顺序写也就是先写第 1 层仰角扫完的几千个径向再写第 2 层也有的厂家先把第 1 层和第 2 层各自扫的中间结果交叉存放最后靠径向头里的“仰角编号”来做重排。我一般取出所有径向头后做一次排序排序的键不是一个 angle 值而是(elevation_index, radial_no)这两个整数靠浮点数来排序容易因为 14.99 和 15.01 这种精度问题把顺序打乱。提示不要在拿到新文件时直接把旧配置硬套进去。每个基数据的“数据段布局”都可能不同尤其是外厂雷达同一型号不同固件版本都能改偏移量。程序里留一个“radar_config”字典里面放头长、库长、字节序、比例因子、是否包含质量标记位换文件先改配置后解析能省掉大量反复打开 hexdump 的时间。3. 基数据解析与物理量换算把字节变成反射率和速度场解析字节只是第一步要把“原始计数”变成“雷达产品”还得做两步换算一是把计数值转换成物理量反射率因子、径向速度、谱宽二是把极坐标网格贴到经纬度或平面投影上。这两个环节做不好出来的图就是一块花斑雷达显示基本没法用。3.1 反射率、速度与谱宽的换算公式反射率的下限和精度一般格式文档里都会给一组阈值比如“小于 2 的计数值视为无效回波”“大于等于 10 的计数值剔除旁瓣干扰”。常见做法是程序启动时读一个标定文件里面包含cal_ref (count - threshold) * 0.5这种线性映射的斜率和偏置或者更常见的查表方式直接将计数值作为索引查一张 256 长度的查询表表中存储对应的 dBZ 值。查询表的好处是雷达厂家可以针对硬件非线性单独修正不至于动解析代码。速度场换算更依赖硬件参数。多普勒速度的换算一般是velocity (count - velocity_nyquist_offset) * velocity_resolution但要注意速度的模糊Nyquist 速度问题如果解析出的速度场在某一径向范围内出现整条从正到负的突变带那基本就是速度模糊需要在显示端做退模糊而数据解析端要保持原始值不被改动。# 反射率和速度换算的参考实现 def convert_reflectivity(counts, lookup_table): # 256 长查询表索引 0..255值单位 dBZ import numpy as np arr np.asarray(counts, dtypenp.uint16) # 无效阈值由文档给出一般是 0 或 1 或 2 invalid_mask arr 255 dBZ lookup_table[np.clip(arr, 0, 255)] dBZ[invalid_mask] np.nan # 无效回波设为 NaN显示时渲染色板跳过 return dBZ def convert_velocity(counts, invalid_threshold, nyquist_velocity): import numpy as np counts np.asarray(counts, dtypenp.float32) counts[counts invalid_threshold] np.nan # 典型公式速度 (count - 阈值) * 速度分辨率 - 奈奎斯特速度 vel (counts - invalid_threshold) * (2.0 * nyquist_velocity / 255.0) - nyquist_velocity return vel这段代码里的lookup_table通常不是解析程序写死的而是由雷达厂家提供的标定文件生成。搜 RadarYHM 相关资源时你会发现不同版本的基数据解析程序带的标定文件有时是不一样的老版本用ref.tbl固定文件新版本改成calibration.xml里面除了反射率还写了仰角补偿和方位增益修正系数。换标定文件后反射率整体偏移 1 到 2 个 dB 是常有的事所以雷达显示程序在左下角会显示标定文件名方便现场比对。3.2 体扫重组按仰角生成扫描层再做坐标变化径向数据处理完之后下一步是“体扫重组”。你需要把整包径向按仰角分组每个仰角组生成一个极坐标矩阵矩阵的形状是(方位径向数, 距离库数)。方位径向数取决于雷达扫描模式如果雷达按 0.5 度方位步进扫描一圈那一圈就是 720 条径向如果是连续扫描并回传角度头径向数可能不是整数最后一条径向和第一条径向的夹角覆盖不到 360 度这时一般要做方位插值。体扫重组的代码如下import numpy as np # radial_list: 每个元素是 {elevation: float, azimuth: float, ref: np.ndarray} def build_elevation_layer(radials, num_azimuths720, max_range_bins460): # 极坐标矩阵第 0 维度是方位第 1 维度是距离库 layer np.full((num_azimuths, max_range_bins), np.nan, dtypenp.float32) for r in radials: azimuth_idx int(round(r[azimuth] / 0.5)) % num_azimuths # 防止某一径向的库数比矩阵宽 bins r[ref][:max_range_bins] layer[azimuth_idx, :len(bins)] bins # 对没有径向覆盖的方位角做最近邻填补 for az in range(num_azimuths): if np.isnan(layer[az]).all(): # 找左右两边最近的完整径向拷贝过来 left az - 1 while left 0 and np.isnan(layer[left]).all(): left - 1 right az 1 while right num_azimuths and np.isnan(layer[right]).all(): right 1 if left 0: layer[az] layer[left] elif right num_azimuths: layer[az] layer[right] return layer这段代码解决的是“扫描径向数不够一圈”的坑。注意这里int(round(r[azimuth] / 0.5))要求方位角头是按 0.5 度为单位存储的如果格式文档里给的是 0.1 度步进那么 azimuth_idx 的计算要把除数改成 0.1。填补最近邻这个操作看着粗暴但在雷达显示里很实用因为不填补会在扇区边界画出一条黑色的“断带”影响对回波连续性的判断。想要更平滑的效果可以改成线性插值但要注意方位角是环形的不能直接线性插值 359 度和 1 度之间的补数要么拆开角度要么把数据复制一份搬到 [0, 720) 再插。4. 从极坐标到屏幕坐标雷达显示的基础映射解析出极坐标场之后真正难的反而是显示这步。雷达显示不能把极坐标数据直接画到屏幕上你必须为每个距离库的经纬度做坐标变换。这里有两个选择一是直接把数据网格化成经纬度网格再做标准的等经纬度渲染二是用 GPU 或图形库直接做极坐标扇区绘制。对从业人员来说前者实现门槛低更贴近大多数气象雷达显示的尺寸和精度。4.1 雷达站点为中心的极坐标到经纬度网格插值把极坐标数据网格化到经纬度核心是正变换已知雷达站点经纬度、某距离库中心的距离、某个径向的方位角求该点的经纬度。然后建立输出网格比如 0.01 度分辨率在每个网格点上判断它落在哪个距离库和哪个方位角之间取对应值。逆插值比正变换更常用因为这相当于在查询经纬度点时直接算出源的极坐标位置。import numpy as np def polar_to_latlon(radar_lon, radar_lat, azimuth_deg, range_km): # 简化的正变换假设地球为球体单位度、千米 R 6371.0 delta range_km / R * 180.0 / np.pi az_rad np.deg2rad(azimuth_deg) lat radar_lat delta * np.cos(az_rad) lon radar_lon delta * np.sin(az_rad) / np.cos(np.deg2rad(radar_lat)) return lon, lat # 逆插值从经纬度网格点反查极坐标方位与距离 def latlon_to_polar(radar_lon, radar_lat, lon, lat): dx (lon - radar_lon) * np.cos(np.deg2rad(radar_lat)) dy lat - radar_lat az (np.degrees(np.arctan2(dx, dy)) 360.0) % 360.0 dist np.sqrt(dx * dx dy * dy) * 111.195 return az, dist这里的正变换公式没有考虑地球曲率和雷达射线的弯曲。说实话对几百公里内的单雷达 PPI 显示来说不做折射修正不会造成肉眼可辨的偏差但对远距离超过 200 公里的强度值按地理位置填到图上会有接近一个网格的偏移。真要做精确的话需要用标准大气折射模型把雷达射线的高度算出来但这只影响高度不影响方位和水平距离。坐标转换之后下一步是插值。距离库方向用“就近取库”方位方向用“线性插值”基本够用。我自己一般先给每个经纬度网格点算出方位角和距离然后直接算它在极坐标矩阵中的浮点坐标再做一次双线性插值这比先正变换到经纬度网格再重采样少一个步骤速度也更快。注意横跨 0 度方位线的处理方位角差要用(az1 - az2 180) % 360 - 180来算不然扇区边界会出裂缝这个细节经常让图上一圈一圈出现“拉丝”的翻车现场。4.2 显示色标与底图叠加的约定雷达显示程序里反射率和速度的颜色是行业规范不是随便挑的。反射率常用“蓝绿黄红紫”档位0 dBZ 以下用蓝黑渐变0 到 5 dBZ 用浅蓝10 到 20 dBZ 用绿色到黄色40 以上直接用红到紫速度场则用正负颜色区分朝向和远离。解析程序只输出物理量数组显示模块负责把物理量映射到颜色表。这里有一个容易做错的地方不要让显示模块把 NaN 值映射成“0”而要映射为透明色或背景色否则远处没回波的地区会显示成一片绿色或紫色这会让新手误判为大范围降水。底图叠加建议使用标准矢量底图分块加载而不是现场去拼在线地图服务。雷达显示页面一般只需要国界、省界、河流和城市点加载完经纬度网格插值数据后直接用图形库按经纬度叠加就行。最多见的定位错误是雷达站点经纬度和底图坐标系不一致基数据文件里如果是 CGCS2000 坐标系的经纬度而底图用的是 WGS84两者在局部地区会差上百米雷达回波会被平移。解析层面要多一步坐标基准统一别把雷达中心点直接当作底图原点来平移。5. 基数据解析和雷达显示避坑五条血泪经验做雷达基数据解析的从业者十有八九都踩过下面这些坑。有些问题看着像数据坏了其实是解析配置错了有些问题看似是显示程序 bug实际上是插值算法的问题。5.1 反射率场全部偏大 10 dBZ标定文件没加载对现象是雷达显示的降水范围看着特别夸张毛毛雨都显示成中等强度降水强回波区几乎全是紫色。你先怀疑标定文件而不是怀疑程序算法。原因很多基数据解析程序源码里有标定文件的默认路径换机器部署时路径不对程序悄悄退回默认值导致查表用的还是出厂标定而雷达硬件近期已经做过标定维护真实偏置变了。有些老格式里反射率计数和 dBZ 的换算斜率是 0.5新格式变成 1.0斜率不对直接把数值翻倍。解决解析程序加一行日志启动时打印实际加载的标定文件路径和查询表第 128 个值对应的 dBZ。如果路径不对去雷达控制软件里读当前标定参数重生成查询表再启动。线上最稳妥的做法是每次雷达定标后导出的calibration.xml都执行一次备份并且把当前正在用的路径写进配置文件的注释里。5.2 雷达图出现整圈稀疏的“齿梳”状纹理径向坐标映射错误现象是强回波区边缘一圈一圈的锯齿尤其是远距离处一帧图里同一个回波团的边缘像被咬了一口每个齿间距完全一致。原因方位径向数按 720 生成但实际基数据雷达是按 0.1 度步进扫描一圈有 3600 条径向。做坐标映射时把 3600 条径向强塞进 720 个方位槽里相当于每 5 条径向数据只保留了 1 条漏掉了 80% 的方位精度。解决先扫描径向头统计方位角的最小差值和最大差值程序用统计结果自动决定方位槽数。不要写死“720”。如果你拿到的是一个多普勒雷达同时做体扫和扇区扫描的数据还会出现某个仰角只有 120 度扇区覆盖这时按 360 度槽位生成全图扇区外全是补出来的假数据显示为一条直线状的假回波边界。对策是要在坐标映射里加一个 coverage_mask记录真实径向覆盖范围显示端只渲染有覆盖的扇区其他置为背景。5.3 速度场出现“黑白棋盘”交错折射修正和速度解模糊冲突现象是速度场图上出现大面积的 0 值附近正负交错看着像棋盘格。原因速度场比较敏感极坐标插值时若对速度场用了反射率场那套“最近邻填补”会把雷达测量不到的远距离库填出一个接近 0 的速度值而这些地方实际没有信号。特别是薄云或层状云降水的边缘插值出的 0 速度与真实的速度场混在一起看起来就是随机正负交替。解决速度场插值时对无覆盖区域一律保持 NaN不做最近邻填充如果有业务需要风场连续性只做范围不超过 2 个距离库的小窗口平滑不要全图插值。另一个坑是速度场先做了退模糊再插值退模糊后的场与原始场跳变特性不同插值引入的伪像更明显。我现在的习惯是原始速度值只做掩码在极坐标域完成退模糊最后才插值到经纬度网格顺序不能反。5.4 体扫数据各个仰角错位仰角分组用错字段现象是 PPI 高仰角图上的回波和一个低仰角图上的回波看起来毫无关联甚至高低仰角的回波位置完全对不上。原因正常体扫中高层仰角看到的回波会相对低层仰角向雷达中心收缩因为雷达射线随距离抬升远距离处低层仰角能看到降水高层仰角已经超出降水顶。如果仰角分组字段选错比如把“俯仰角原始计数值”当作“仰角索引”使用会把不同扫描层的径向混在一起显示端画出的是几层数据的杂合体位置自然对不上。解决读取径向头之后画一张仰角索引-径向序号的关系表比如打印(radial_no, elevation_index)的前 100 条正常情况应该是仰角索引不变、径向序号递增扫完一层后仰角索引跳变。如果前 100 条里仰角索引来回跳检查一下解析的字段偏移量是不是偏了一位。经验是绝大多数解析错误都能通过打印前 10 条径向头元数据排查出来先看元数据再看数据。5.5 远程显示卡顿经纬度重采样太慢现象是雷达站本机显示流畅但远程查看端一加载体扫就卡死浏览器内存直接拉满。原因气象雷达基数据一张体扫如果是 14 层、每层 3600 条径向、每条 800 个库那每个仰角矩阵是 3600 × 800 个 float14 层接近 1.6 亿个数值。如果远程显示端在浏览器里逐个采样点做经纬度反算插值即使每层重采样只算一次也要循环几千万次浏览器直接崩。解决常见的做法是服务器端先把体扫数据在极坐标域插值成经纬度网格比如每层 600×700 网格然后对每个网格编码成 PNG 瓦片或者 GeoTIFF 切片显示端只拼图不重算插值。这样远程显示端只接收图片内存占用和 CPU 消耗都压到很低。雷达显示速度优化越靠近数据源头做的计算越多越靠近用户做的计算越少这条原则不会错。6. 进阶用法用解析出的基数据做 CAPPI 和剖面图验证解析正确性当 PPI 能正常显示之后你其实还可以用同一套解析结果继续做 CAPPI等高平面位置显示和垂直剖面。这两个产品对解析正确性的验证远远强于肉眼看 PPI因为 CAPPI 会把各个仰角的数据融合到同一个高度面上如果仰角元数据解析有误CAPPI 会出现明显的“接缝”或“断崖”。CAPPI 的简单做法是对每个经纬度网格点计算该点在每个仰角上的高度用标准大气折射近似公式h sqrt(R^2 r^2 2Rr) - R其中 R 是等效地球半径 8500 公里r 是斜距然后找出距离目标高度最近上下两个仰角的反射率值做高度线性插值。其实不需要对每个点都算三角函数可以先把每个仰角的距离库先转成高度层做成“高度-库数”映射表经纬度网格点只做一次径向距离查找再查表取高度这样性能能快五倍左右。剖面图的实现更简单但也更考验解析质量选一条穿过雷达站的线段沿线段每个点用 4.1 里的逆变换找到它在每个仰角层中的位置把不同仰角的数值拼到一个二维剖面上。做剖面图之前务必先验证速度场和反射率场是否对齐——有的基数据格式里速度场距离库的起点比反射率场晚一个库不做对齐直接画剖面会在近地面出现一条假的水平亮带。最后想分享一个习惯不管从哪个渠道拿到 RadarYHM 这类基数据解析程序第一件事不是跑代码而是先做一个最小验证——生成一帧 PPI存为 GeoTIFF 带坐标信息叠加在卫星云图或者上一帧雷达图上对比一分钟。这样做的目的是确认坐标转换和仰角分组没问题而不是只是“看到图了”。坐标差半度、仰角错一层这类问题肉眼看实时动画很难发现但叠加底图或者叠加前一时次数据一眼就能看穿。希望这篇笔记帮你少走这一段弯路。本文还有配套的精品资源点击获取
