C++实现SAR原始数据成像：从RAW到SLC的聚焦处理

简介：面向合成孔径雷达图像处理的C++源码工程，从原始回波数据开始，覆盖数据预处理、聚焦成像、去噪、特征提取、图像增强以及格式转换等完整处理环节。适合遥感科学与技术专业的学生、雷达信号处理方向的研究人员，以及需要借助高效语言处理大规模图像数据的开发者，有助于理解从原始数据到最终图像的算法工程化落地。压缩包内共三十四个文件，以头文件和源文件为主体，同时包含工程配置文件、界面资源、说明文档与示例位图，整体仅一百零八KB，结构紧凑，可快速查看和编译学习。当前已有四百五十一人学习浏览。通过逐段研读代码，能够掌握数据读取、几何校正、辐射校正、频域聚焦等关键模块的实现思路，也能积累在相关框架下进行图像显示、文件交互和工程管理的实战经验，是完整且便于对照学习的样例。

1. RAW 格式的 SAR 数据不是图像，C++ 离 SLC 还差一次聚焦

用二进制查看器打开一包 SAR 原始回波，看到的不是图像，而是大量看似随机的 16 位整数。合成孔径雷达和光学传感器完全不同：回波必须依次经过距离向脉冲压缩、方位向合成孔径聚焦，才能得到可读的 SLC 复数图像；想要输出人眼能看的影像，还要再取模、多视、量化。标题里这套“从 RAW 格式起”的 C++ 处理代码，本质就是把这条链路用编译型语言完整串起来，替代 MATLAB 原型，获得更高的吞吐和更好的集成能力。适合两类人：一类是准备把算法从 MATLAB 迁移到 C++ 的雷达工程师，另一类是要把 SAR 成像作为模块嵌进业务系统的软件工程师。后面的内容就按数据读取、聚焦成像、后处理、验证四段推进，每一段都给出参数边界和可执行的代码。

2. RAW 数据格式解析：先从二进制流里读出复数回波

2.1 元数据决定怎么读，别上来就写解析器

SAR 原始回波以复数采样形式保存，每个采样点由一个 I 路和一个实部、一个 Q 路虚部组成。常见组织方式有两种：I/Q 交织存储，每个采样点的 I、Q 连续存放，16 位量化时 4 个字节构成一个复数；少数格式采用 I/Q 分块，先存所有 I 再存所有 Q。读取代码怎么写的决定因素不是文件后缀，而是旁边的 .hdr、.par 元数据文件。这里还要先澄清一点：检索“RAW 格式”时经常把 U 盘变成 RAW 分区的话题一起带进来，雷达领域说的 RAW 是原始回波数据，和文件系统损坏没有关系。

元数据字段决定了后续每一行代码，处理前必须逐项核对。

PRF 和平台速度必须对得上：如果 PRF 填错，方位向会被整体拉伸或压缩，图像里的点目标会变成一条斜线。星载 SAR 的元数据里还有轨道和时间参数，机载或无人机采集的数据则会有速度、高度、斜距起点。这些值任何一个错了，后面的聚焦都白做。

2.2 用 C++ 读取 16 位 I/Q 交织的最小实现

下面这段代码是完整链路的第一步：把二进制文件逐脉冲读入内存，转成std::complex<float>数组。编译调试阶段，在 VS Code 里配好 C/C++ 环境后，从这个函数开始断点排查比较顺手。

#include <cstdint> #include <fstream> #include <vector> #include <complex> #include <stdexcept> struct SarRawHeader { uint32_t samples_per_pulse; // 距离向每脉冲采样点数 uint32_t num_pulses; // 方位向脉冲总数 bool iq_interleaved; // I/Q 是否交织，预期为 true }; std::vector<std::complex<float>> readRawIq16( const std::string& path, const SarRawHeader& hdr) { std::ifstream ifs(path, std::ios::binary); if (!ifs) throw std::runtime_error("can't open raw file"); std::vector<std::complex<float>> data( hdr.samples_per_pulse * hdr.num_pulses); std::vector<int16_t> buf(2 * hdr.samples_per_pulse); for (uint32_t i = 0; i < hdr.num_pulses; ++i) { ifs.read(reinterpret_cast<char*>(buf.data()), buf.size() * sizeof(int16_t)); if (ifs.gcount() != static_cast<std::streamsize>(buf.size() * sizeof(int16_t))) { throw std::runtime_error("file truncated at pulse " + std::to_string(i)); } for (uint32_t j = 0; j < hdr.samples_per_pulse; ++j) { data[i * hdr.samples_per_pulse + j] = std::complex<float>(buf[2 * j], buf[2 * j + 1]); } } return data; }

这段代码的逻辑很清楚：每次读入一个脉冲的 I/Q 原始整数，再按“奇数位置是 Q、偶数位置是 I”的布局转成复数。std::complex<float>的内存布局与两个连续 float 一致，后续传给 FFT 库时可以直接用reinterpret_cast。如果数据超过 2 GB，常见做法是改用mmap按行读取，避免一次性分配巨大 vector；如果量化位数为 8 bit，把int16_t换成int8_t，并在构造复数时除以 128.0f 做归一化。

2.3 行距、辅助字节与尾部检查

不少真实数据文件会在每个脉冲之间夹几个辅助字节，比如 GPS 时间戳、事件计数、RCS 标定值。元数据里如果给了“行距”或“行字节数”，读取时就不能简单地顺序读，而要seekg到每行起始位置。调试时先取前 1000 个脉冲、每脉冲 1024 个采样点跑通流程，比直接灌完整帧数据高效得多。读取结束时再核对实际脉冲数是否等于元数据声明值，文件尾部多出的零填充对聚焦没有影响，可以忽略。

提示：所有读入的数据先打印第一个脉冲前 16 个复数的实部虚部，如果数值范围在几千到几万之间，说明量化正常；如果全是 0 或者全部是同一个值，说明字节序或行距配错了。

3. 距离压缩与方位压缩：把回波聚焦成 SLC 图像

3.1 匹配滤波的频域实现为什么是标准做法

SAR 发射的是线性调频信号，基带形式为s(t) = exp(j * PI * Kr * t²)，其中Kr是距离向调频率，Tp是脉冲宽度。回波是一个延时后的线性调频信号，时域上和匹配滤波器做卷积可以压成窄脉冲，带宽B = |Kr| * Tp决定距离分辨率rho_r = c / (2B)。匹配滤波时域卷积的复杂度是 O(N²)，而 FFT 是 O(N log N)。一条回波长度通常有几千个采样点，整个场景有数万条脉冲，频域实现省下的时间不止一个量级。频域做法是固定的：对回波做 FFT，乘参考信号频谱的共轭，再 IFFT 回来。参考信号直接构造为发射波形s(t)，频域里乘它的共轭即可。

#include <fftw3.h> #include <vector> #include <complex> #include <cmath> // 构造距离向参考信号：直接取发射波形 s(t) std::vector<std::complex<float>> buildRangeRef( int refLen, float fs, float Kr) { std::vector<std::complex<float>> ref(refLen); for (int n = 0; n < refLen; ++n) { float t = (n - refLen / 2) / fs; // 以序列中心为时间零点 ref[n] = std::exp(std::complex<float>(0.0f, PI * Kr * t * t)); } return ref; } // 单条回波做距离向匹配滤波，原位修改 line void rangeCompress(std::vector<std::complex<float>>& line, const std::vector<std::complex<float>>& ref) { const int N = static_cast<int>(line.size()); std::vector<std::complex<float>> refPad(N, {0.0f, 0.0f}); for (int i = 0; i < N && i < static_cast<int>(ref.size()); ++i) { refPad[i] = ref[i]; } fftwf_plan pFwd = fftwf_plan_dft_1d( N, reinterpret_cast<fftwf_complex*>(line.data()), reinterpret_cast<fftwf_complex*>(line.data()), FFTW_FORWARD, FFTW_ESTIMATE); fftwf_plan pRef = fftwf_plan_dft_1d( N, reinterpret_cast<fftwf_complex*>(refPad.data()), reinterpret_cast<fftwf_complex*>(refPad.data()), FFTW_FORWARD, FFTW_ESTIMATE); fftwf_execute(pFwd); fftwf_execute(pRef); for (int i = 0; i < N; ++i) { line[i] *= std::conj(refPad[i]); // 频域匹配滤波 } fftwf_plan pInv = fftwf_plan_dft_1d( N, reinterpret_cast<fftwf_complex*>(line.data()), reinterpret_cast<fftwf_complex*>(line.data()), FFTW_BACKWARD, FFTW_ESTIMATE); fftwf_execute(pInv); const float norm = 1.0f / N; for (auto& v : line) v *= norm; // IFFT 归一化 fftwf_destroy_plan(pFwd); fftwf_destroy_plan(pRef); fftwf_destroy_plan(pInv); }

这段代码里line既是输入也是输出，FFTW 的原位变换允许这样用。refPad补零到和回波等长，是为了让循环长度是 2 的幂或者 FFTW 擅长的长度。距离压缩是以脉冲为独立单位的，各脉冲之间没有依赖，可以自然并行。如果场景很大，调用处用 OpenMP 对脉冲循环做并行，加速比接近核心数。fftwf_plan建议放进类成员或函数外部复用，否则每次创建 plan 都有额外开销。

3.2 距离压缩后的二维布局与方位压缩

距离压缩完成后，数据布局是“行 = 脉冲序号，列 = 距离采样点”，这时多普勒信息还压在列方向里。要把方位向聚焦做出来，先转置，让每一行对应一个固定距离单元、每一列对应方位慢时间，再对每行做方位压缩。这相当于一个可分离的二维匹配滤波：先距离后方位。

// 对某一个距离单元的方位向信号做匹配滤波 void azimuthCompress(std::vector<std::complex<float>>& azLine, float prf, float ka) { const int M = static_cast<int>(azLine.size()); fftwf_plan pFwd = fftwf_plan_dft_1d( M, reinterpret_cast<fftwf_complex*>(azLine.data()), reinterpret_cast<fftwf_complex*>(azLine.data()), FFTW_FORWARD, FFTW_ESTIMATE); fftwf_execute(pFwd); for (int i = 0; i < M; ++i) { float fa = (i <= M / 2 ? i : i - M) * prf / M; // 方位频率轴 azLine[i] *= std::exp(std::complex<float>( 0.0f, PI * fa * fa / ka)); // 匹配滤波相位 } fftwf_plan pInv = fftwf_plan_dft_1d( M, reinterpret_cast<fftwf_complex*>(azLine.data()), reinterpret_cast<fftwf_complex*>(azLine.data()), FFTW_BACKWARD, FFTW_ESTIMATE); fftwf_execute(pInv); const float norm = 1.0f / M; for (auto& v : azLine) v *= norm; fftwf_destroy_plan(pFwd); fftwf_destroy_plan(pInv); }

ka是方位向调频率，按ka = 2 * V² / (lambda * R0)计算，V是平台速度，lambda是波长，R0是场景中心斜距。这里最容易出错的点是符号约定：不同教材里ka正负定义不一致，代码里的相位符号要和回波多普勒变化方向匹配。最稳的验证方式是用一个模拟点目标测试，而不是拿着真实数据看效果。

聚焦参数表中，距离向和方位向的关键量需要一并核对。

真实的全套处理代码不会只有两步压缩。如果目标是星载 SAR，地球自转、轨道弯曲会引入额外的相位误差，必须用轨道状态矢量做精确斜距建模，一阶近似无法聚焦到理论分辨率。机载数据在低空大斜视角下还会有距离单元徙动，最简做法是在距离频域对方位时间做线性相位补偿，也就是常说的 RCMC。很多网上的处理代码省略了这一步，结果是点目标被拉成一条弧线，方位向分辨率明显变差。

4. 幅度图像生成与多视：处理代码里的工程细节

4.1 多视的公式和实现

聚焦完成后的数据是 SLC，每个像素仍是复数。直接取模可以看到图像轮廓，但斑点噪声很强，肉眼难以判读。多视就是牺牲一部分分辨率，把若干相邻像素的功率平均后再开方，换取斑点噪声的抑制。典型配置是距离向 2 至 4 个像素、方位向 2 至 4 个像素合成一个输出像素。注意这里算的是功率平均，不是幅度平均，即先对各像素幅度平方求和，再除以视数后开方，这与“多视”的统计意义一致。

#include <vector> #include <complex> #include <cmath> // 对 SLC 复数数据做多视，输出幅度浮点图像 std::vector<float> multiLook( const std::vector<std::complex<float>>& slc, size_t width, size_t height, size_t rangeLooks, size_t azimuthLooks) { const size_t outW = width / rangeLooks; const size_t outH = height / azimuthLooks; std::vector<float> out(outW * outH, 0.0f); for (size_t az = 0; az < outH; ++az) { for (size_t rg = 0; rg < outW; ++rg) { float acc = 0.0f; for (size_t a = 0; a < azimuthLooks; ++a) { for (size_t r = 0; r < rangeLooks; ++r) { const auto& v = slc[(az * azimuthLooks + a) * width + rg * rangeLooks + r]; acc += v.real() * v.real() + v.imag() * v.imag(); } } out[az * outW + rg] = std::sqrt(acc / (rangeLooks * azimuthLooks)); } } return out; }

rangeLooks和azimuthLooks的取值需要结合任务权衡。

多视之后的分辨率变为原来的rangeLooks倍和azimuthLooks倍，这个损失是显式的。如果后处理还要做边缘检测或目标识别，视数不宜取大，通常 2 或 4 是比较均衡的点。

4.2 动态范围压缩和图像输出

SAR 幅度图动态范围很大，地物后向散射差异可达 30 ~ 60 dB，直接线性映射到 8 位图会让暗区完全看不到。常见做法是先取 dB，再按直方图分位数做裁剪，最后线性映射到 0 ~ 255。下面代码实现的是线性域的 1%/99% 裁剪拉伸，效果与先取 dB 再用分位数类似，区别是线性拉伸对强散射目标更敏感，dB 拉伸对弱散射区域更友好。

#include <algorithm> #include <vector> #include <cstdint> // 用分位数裁剪后线性映射到 uint8 std::vector<uint8_t> stretchTo8bit( const std::vector<float>& amp, float lowPct, float highPct) { std::vector<float> sorted = amp; std::sort(sorted.begin(), sorted.end()); const float amin = sorted[static_cast<size_t>(sorted.size() * lowPct)]; const float amax = sorted[static_cast<size_t>(sorted.size() * highPct)]; if (amax <= amin) return std::vector<uint8_t>(amp.size(), 0); std::vector<uint8_t> out(amp.size()); for (size_t i = 0; i < amp.size(); ++i) { float v = (amp[i] - amin) / (amax - amin); out[i] = static_cast<uint8_t>(v * 255.0f); } return out; }

lowPct=0.01、highPct=0.99适合大多数场景；城区有大量强反射体时，可以把高端改成 0.995 甚至 0.999，避免少数亮目标占用整个灰度范围。如果要输出 16 位 TIFF，把 255 换成 65535，并把中间变量改成uint16_t。如果项目里已经集成了 OpenCV，cv::imwrite可以直接输出 PNG 和 JPG；如果只是验证算法，写 16 位灰度 PGM 最省事，任何看图软件都能打开。FPGA 落地场景里，这一步通常会被重写为流水线式的定点运算，但分位数裁剪的逻辑保持不变。

5. 验证与调试：成像算法的验收从点目标开始

5.1 点目标切片测量，判断聚焦是否到位

不管代码从哪来，成像结果的验收第一件事不是看整幅图顺不顺眼，而是找点目标做切片分析。取幅度图全局峰值附近 64 x 64 的切片，分别沿距离向和方位向画功率剖面。理想情况下，主瓣 -3 dB 宽度应接近理论值：距离向为c / (2B)，方位向为V / PRF。旁瓣峰值比理想 sinc 脉压约为 -13.2 dB，加窗后旁瓣更低但主瓣会展宽。如果实测主瓣宽度是理论的几倍，问题几乎都出在调频率、平台速度或斜距这三个参数上；如果旁瓣明显不对称，先查参考信号的时间零点是否对准。

// 提取点目标附近切片并统计 -3dB 宽度（示意代码） auto it = std::max_element(amp.begin(), amp.end()); int peakIdx = static_cast<int>(it - amp.begin()); int peakRow = peakIdx / width; int peakCol = peakIdx % width; float peakPower = (*it) * (*it); int halfR = 0; for (int r = 0; r < 32; ++r) { float p = amp[(peakRow + r) * width + peakCol]; if (p * p < peakPower / 2.0f) { halfR = r; break; } }

上面的循环只做了一侧主瓣的粗测。完整做法还要把另一侧也量出来，并换算成实际距离或方位长度。多视处理后的幅度图不宜做这类验证，因为视数已经把主瓣展宽了，点目标分析应放在 SLC 幅度或单视幅度图上做。

5.2 拿到一份 SAR C++ 工程后先查这几个点

别人写的处理代码拿到手后，不要先急着换文件路径跑通，先按下面的顺序排查：第一步，确认读入脉冲数和元数据一致；第二步，找一个独立实现的 FFT 结果对照距离压缩输出；第三步，检查转置前后的索引关系有没有把距离向和方位向弄反；第四步，确认 IFFT 后有没有做 1/N 归一化，FFTW 的逆变换默认是不归一化的。

提示：聚焦结果出现“距离向清晰、方位向糊”时，优先级依次检查方位向调频率ka、RCMC 是否缺失，再怀疑方位向加窗过度。如果是“完全没聚焦”，先回到模拟点目标测试，用参数直接生成理想回波，看看处理链本身能不能自洽。

在高分辨率场景下，参数估计误差很难完全消除，这时候需要自动聚焦。最常用的是基于图像熵或图像对比度的最优化方法：构造一个多项式形式的相位误差，搜索多项式系数使得图像熵最小。熵的定义是E = -sum(p_k * log(p_k))，其中p_k是归一化强度。系数搜索用共轭梯度或 Nelder-Mead 都可以，每次迭代都要对整幅图像做一次方位向逆滤波再评价熵，计算量不小，但这是把中等精度代码推向高分辨率成像的必经之路。相位误差的多项式阶数一般取 2 到 3 阶，再多就会把真实地形细节也当成误差消掉。

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