简介一套覆盖合成孔径雷达三大方向的Matlab源码合集面向遥感和信号处理方向的学生与研究人员聚焦SAR成像、InSAR干涉处理及PolSAR极化定标可满足毕业设计、课程设计和算法复现需求。包内共191个文件以116个m源码为核心辅以35个mat仿真数据、21个pdf原理文档、13个txt运行说明和3个md总结笔记整体约161.84MB目录按SAR、InSAR、PolSAR模块划分并附带理论设计示意图便于对照学习。SAR部分实现RD与CS两种经典成像算法涵盖点目标仿真及Radarsat-1实测数据处理InSAR部分针对平地与圆锥两种人造场景给出从原始回波仿真、成像到干涉处理的完整流程PolSAR部分提供极化定标算法并配有项目文档与算法解析。所有源码均经过测试已有234人学习适合在此基础上扩展也适合作为相关课题的起点。1. 合成孔径雷达源码库里的三套链路到底在解决什么问题拿到一份合成孔径雷达源码包你会发现目录里通常躺着三套互不相同的代码SAR 成像、InSAR 干涉处理、PolSAR 极化分析。很多人以为这是同一个处理流程的三个模块实际上一上来就各自为政。SAR 解决的是从原始回波到一张高质量图像InSAR 解决的是用两张或多张图像反演地表形变与高程PolSAR 则直接把后向散射拆成多种散射机制——植被、裸地、建筑物、水体在不同的极化通道里呈现出完全不同的响应。这三件事共享同一个平台运动模型和信号模型但数据处理几何关系各不相同。本文把这套源码拆开讲先跑通 RD 成像算法拿到聚焦图像再做人造场景的原始回波仿真与干涉处理最后延伸到极化分解。适合要做 SAR 课程设计、星载 SAR 数据预处理、或者想把干涉测量误差源头定位到仿真阶段的工程师。下文给出可复现代码与参数边界。2. SAR成像算法从原始回波到二维图像的RD成像原理与最小实现2.1 为什么成像算法的核心卡在距离徙动校正合成孔径雷达在工作时平台一直向前运动雷达波束以正侧视或斜视方式照射地面目标。目标到天线之间的斜距随时间变化这带来两个独立问题一是方位向的多普勒历史这是 SAR 为什么能有极高方位分辨率的物理基础二是由于斜距变化导致的距离单元走动和弯曲统称为距离徙动。脉冲压缩原理使距离向分辨率做到足够高但若在距离向压缩之后直接把每一行数据按慢时间摆放目标能量会跨距离单元蔓延图像就是模糊。距离徙动校正是 RD 成像算法中决定聚焦质量的核心步骤。RD 算法通过在距离频域乘以一个与方位频率相关的校正因子等效地把弯曲的轨迹拉直这之后方位向匹配滤波才能正确压缩。成像算法的完整流程可以用三个基本步骤概括范围向压缩、距离徙动校正、方位向压缩。把每一列数据沿慢时间做一条一维匹配滤波目标能量集中于一点。常见做法是先对原始回波数据做范围向快速傅里叶变换在频域乘以匹配滤波器后变换回时域然后在距离频域完成距离徙动校正最后沿方位向做匹配滤波。这样做的好处是每一步都可以独立验证中间结果的物理意义出现问题容易定位。2.2 RD算法三步拆解与Python实现2.2.1 距离向压缩的最小实现先构造线性调频信号作为发射波形import numpy as np def range_compress(raw_data, Kr, Tr, fs): 距离向压缩频域匹配滤波 raw_data: (方位向, 距离向) 复数原始数据 Kr: 调频率 (Hz/s) Tr: 脉冲宽度 (s) fs: 距离向采样率 (Hz) nr raw_data.shape[1] # 构造匹配滤波器参考函数 t np.arange(-nr//2, nr//2) / fs ref np.exp(-1j * np.pi * Kr * t**2) ref np.conj(ref) # 匹配滤波 参考信号的共轭反转 ref_fft np.fft.fft(ref, nr) range_fft np.fft.fft(raw_data, axis1) compressed np.fft.ifft(range_fft * ref_fft, axis1) return compressed这段代码里ref的构建是重点发射信号是正调频率的线性调频脉冲匹配滤波器是它的时间反转共轭。在频域里这个操作等价于对所有频率分量施加与信号相位共轭的相位变化使得目标能量集中在零距离单元附近。注意这里的t是从-nr/2到nr/2因为匹配滤波器的群延迟天然会使峰值出现在参考时间位置。如果你发现压缩后的峰值位置整体偏移多半是这里的时间轴定义与发射信号的时间基准不匹配。2.2.2 距离徙动校正的插值实现距离压缩完成后每一个方位慢时间时刻的回波峰值能量实际上落在不同的距离单元里。对这个量做二维数组上的插值搬移就是距离徙动校正的常见做法。def rcm_correction(compressed, kr, ka, fdc, v_r, lambda_c, nr): 距离徙动校正距离频域插值实现 compressed: 距离压缩后的二维数组 v_r: 等效雷达速度 (m/s) fdc: 多普勒中心频率 (Hz) n_az, n_rg compressed.shape f_range np.fft.fftfreq(n_rg, d1/fs) # 距离频率 # 先沿方位向做FFT进入距离多普勒域 rd_domain np.fft.fft(compressed, axis0) f_az np.fft.fftfreq(n_az, d1/prf) # 方位频率 for i, fa in enumerate(f_az): # 计算目标在距离多普勒域中的斜距 R_fa lambda_c * (fdc - fa) / (2 * v_r**2) * (fdc - fa) / (4 * (fdc - fa)**2 / lambda_c**2) # 实际工程中用的是近似式 R_fa lambda_c * (fdc - fa) / (4*ka) # 距离偏移量以距离采样单元为单位 delta_r R_fa - np.mean(R_fa) shift_bins delta_r * 2 * fs / c # 在距离频域乘以线性相位实现时域搬移 phase np.exp(1j * 2 * np.pi * f_range * shift_bins / fs) rd_domain[i, :] np.fft.ifft(np.fft.fft(rd_domain[i, :]) * phase) # 反变换回二维时域 corrected np.fft.ifft(rd_domain, axis0) return corrected这段代码的关键点是shift_bins的计算距离频域乘以线性相位exp(j2π f_range shift_bins / fs)等价于时域里把信号搬移shift_bins个采样点。这样一来不同方位频率分量在各目的距离频域中独立完成偏移不需要做逐像素的插值循环。实际参数下shift_bins通常在零点几到几十个距离单元之间。平台轨道越高、斜距越远徙动量就越大不做校正时目标能量会倾斜分布在多个距离单元中导致方位向压缩后的旁瓣不对称。2.2.3 方位向压缩距离徙动校正完成后信号在距离多普勒域里已经是一条直线。方位向压缩与距离向压缩形式完全相同只是参考函数换成方位向的多普勒调频率def azimuth_compress(corrected, prf, ka, t_az): 方位向压缩 n_az corrected.shape[0] # 方位向参考函数线性调频 ref_az np.exp(-1j * np.pi * ka * np.linspace(-t_az/2, t_az/2, n_az)**2) ref_az np.conj(ref_az) # 逐距离单元做方位向匹配滤波 az_fft np.fft.fft(ref_az, n_az) result np.zeros_like(corrected) for r in range(corrected.shape[1]): row_fft np.fft.fft(corrected[:, r]) result[:, r] np.fft.ifft(row_fft * az_fft) return resultka代表方位向多普勒调频率直接由平台速度与斜距决定。星载 SAR 的ka值远大于机载意味着方位向积累时间更短所需的脉冲数更少但相位误差对运动的敏感度也更高。2.3 成像后的质量检查点扩散函数与分辨率验证用点目标测试合成孔径雷达成像算法是否实现正确需要通过点扩散函数的 3 个指标验证峰值旁瓣比、积分旁瓣比、分辨率。下表是通用检查标准指标理想值实际可接受范围说明方位向峰值旁瓣比-13.26 dB优于 -12 dB过低的旁瓣意味着加窗这会同时展宽主瓣距离向峰值旁瓣比-13.26 dB优于 -12 dB未加窗时与理论值接近方位向分辨率λR/(2L)理论值的 1.0~1.4 倍L 为合成孔径长度加窗后会展宽距离向分辨率c/(2B)理论值的 1.0~1.4 倍B 为发射信号带宽成像完成后找孤立点目标取峰值周围 9×9 像素区域做插值搜索主瓣峰值位置量测 -3 dB 处的宽度。如果测出的旁瓣比明显低于 -13 dB 且主瓣不对称优先怀疑距离徙动校正里shift_bins的计算公式是否有符号错误。3. InSAR人造场景回波仿真与干涉处理从场景建模到干涉图完整闭环3.1 人造场景的几何精度可控性优势干涉合成孔径雷达测量InSAR需要从两幅复图像中提取相位差从而反演高程或微小形变。实测数据面临的难点在于真实地形精确已知度有限、轨道误差与大气延迟会混入干涉相位。人造场景的好处是地面目标的三维坐标、后向散射系数、形变场完全由脚本设置仿真回波经过成像后得到的干涉相位可以与初始模型严格比对。这样就把测量链路的误差源头从环境误差中分离出来。在做源码级实验时我一般会用程序生成一个包含已知高程梯度的数字高程模型DEM如斜坡、台阶、凸起三种形态再将 DEM 叠加到基准平面。这种方式能快速区分成像算法误差与干涉处理误差。3.2 人造场景的原始回波仿真与平台参数表原始回波仿真的基本思路是把场景离散成散射点对每个散射点计算其到雷达的瞬时斜距生成对应的回波并叠加到原始数据矩阵中。接收信号模型可以写成s_echo(t, η) Σ A_i · rect((t - 2R_i(η)/c) / Tr) · exp(-j·4π·R_i(η)/λ) · w_a(η - η_c,i) n(t, η)其中 R_i(η) 是第 i 个散射点在方位慢时间 η 时的斜距A_i 包含后向散射系数与天线方向图加权n(t, η) 是热噪声项。# 点阵场景回波仿真 Python 骨架 def simulate_raw_echo(scene_points, az_time, range_time, params): scene_points: [(x, y, z, sigma), ...] 地面点坐标 az_time: 方位向慢时间向量 (s) range_time: 距离向快时间向量 (s) c 3e8 raw np.zeros((len(az_time), len(range_time)), dtypecomplex) for (x0, y0, z0, sigma) in scene_points: # 正侧视几何雷达位于 (0, v_r * eta, H) # 目标在雷达坐标中的斜距 R0 np.sqrt((params[H] - z0)**2 x0**2 y0**2) for i, eta in enumerate(az_time): # 平台位置 R np.sqrt((params[H] - z0)**2 x0**2 (y0 - v_r*eta)**2) # 距离向回波位置 tau 2*R/c idx np.argmin(np.abs(range_time - tau)) if abs(range_time[idx] - tau) 1.0/fs: # 线性调频信号回波 phase -4*np.pi*R/params[lambda] raw[i, idx] sigma * np.exp(1j*phase) return raw这段仿真代码的关键参数需要显式设置参数符号推荐值作用平台速度v_r150 m/s方位向积累长度与多普勒调频率控制平台高度H3000 m斜距范围控制载频fc9.6 GHz波长 → 干涉相位对高程的敏感度脉冲宽度Tr0.5 μs距离向分辨率的一阶控制距离向带宽B150 MHz距离向分辨率脉冲重复频率PRF300 Hz方位向采样率地面点数N500~2000点数过多仿真会慢一个数量级平台速度与 PRF 的组合关系必须满足方位向采样定理PRF 必须大于方位向多普勒带宽的两倍否则会因欠采样出现方位模糊。仿真完成后回波数据直接送入第 2 章的 RD 成像算法得到两幅单视复图像。人工场景不需要额外做地形相位补偿这一步这正是仿真数据的便利之处。3.3 干涉处理最小流程配准、干涉图、滤波与解缠3.3.1 主辅图像亚像素配准干涉测量要求两幅复图像在像素级对齐否则相干性显著下降。常用做法是先通过星历参数估算一个粗配准偏移量再做基于相关函数的亚像素配准def cross_correlation_register(master, slave, search_range8): 基于幅度互相关的亚像素配准单点计算示例 # 在频域计算移位后的互相关 m_fft np.fft.fft2(master) s_fft np.fft.fft2(slave) # 幅度归一化互相关 cross m_fft * np.conj(s_fft) cross_map np.abs(np.fft.ifft2(cross / np.abs(cross))) max_idx np.unravel_index(np.argmax(cross_map), cross_map.shape) shift (max_idx[0] - cross_map.shape[0]//2, max_idx[1] - cross_map.shape[1]//2) return shift距离向配准误差在 0.1 像素以内方位向在 0.1 像素以内才能保证相干性损失在可接受范围内。实际操作中如果配准偏移量在方位向超过 1 个像素说明两幅数据的零多普勒时刻偏移没算准这时候要先检查轨道插值参数。3.3.2 多视处理与相位解缠配准后按水平和垂直方向的视数对干涉图做多视平滑。对于地形平坦的仿真场景干涉图主要体现目标高程引起的平缓相位变化此时应采用加权圆均值滤波窗口大小 3×3 与 5×5 两种def multilook_filter(interferogram, look_range3, look_az3): 多视滤波 from scipy.ndimage import uniform_filter phase np.angle(interferogram) mag np.abs(interferogram) # 拆成相互独立的实部/虚部分别做窗口均值 re uniform_filter(mag*np.cos(phase), size(look_az, look_range)) im uniform_filter(mag*np.sin(phase), size(look_az, look_range)) return re 1j*im相位解缠是干涉处理中最容易出问题的一环。不同解缠算法的选择直接影响高程反演精度解缠算法适用场景已知局限工程推荐参数枝切法低噪声、相干性高、残差点稀少残差点密集时出现大面积非连通区域枝切长度上限32 像素最小二乘残差点较多、区域连续高程跳变处相位被平滑加权系数设为相干性区域增长法不规则区域边界初始种子点选择敏感种子阈值0.5SNAPHU统计费用流形变场复杂、高噪声耗时显著需调参数默认参数即可解缠后的相位通过 ∮ -(4πB⊥)/(λR sinθ) · h 反演高程。在卫星轨道安全、基线估计精确的仿真条件下如果解缠相位与真值 DEM 的残差仍然呈条纹状分布说明基线参数与成像几何参数不一致需要回到平台参数设置中检查基线相位贡献。4. PolSAR极化定标与特征分解四通道不是四个单通道4.1 合成孔径雷达极化测量为什么单独占一套源码SAR 与 InSAR 处理的对象是标量复数数据PolSAR 处理的对象是每个像素一个 2×2 复数散射矩阵包含 HH、HV、VH、VV 四个极化通道。这四个通道之间不是独立的电磁波在散射过程中发生极化状态变化矩阵各元素之间存在明确的相位关系。直接忽略交叉极化通道、只保留同极化幅度信息会丢失地物分类中的关键物理线索。PolSAR 源码要解决的第一个问题是极化定标。系统硬件在发射和接收路径上会引入通道不平衡和交叉串扰如果不校准后续分解得到的散射功率分布完全不可信。4.2 极化定标矩阵与最小验证集雷达测量得到的极化散射矩阵 M 与真实散射矩阵 S 的关系为M A·R·S·T N其中 T 与 R 分别是发射和接收失真矩阵。做定标时常用一个二面角反射器或三面角反射器其理论散射矩阵已知。定标后计算剩余误差如幅度不平衡优于 0.3 dB串扰优于 -30 dB。def pauli_decompose(hh, hv, vh, vv): Pauli基分解把散射矩阵投影到Pauli基上 Pauli基奇次散射、偶次散射、45度偶次散射、螺旋散射 # 对应 S a[1 0;0 1] b[1 0;0 -1] c[0 1;1 0] d[0 -1j;1j 0] a (hh vv) / np.sqrt(2) # 奇次散射球体、三面角 b (hh - vv) / np.sqrt(2) # 偶次散射二面角 c (hv vh) / np.sqrt(2) # 45°旋转的偶次散射 d (hv - vh) / np.sqrt(2) # 螺旋散射分量 # 散射功率 p_a np.abs(a)**2 p_b np.abs(b)**2 p_c np.abs(c)**2 return p_a, p_b, p_c, np.ones_like(p_a)*np.abs(d)**2对于一个三面角反射器理论上 p_b、p_c 应为零实际中如果 p_b 与 p_a 比值高于 -20 dB需要回去排查接收通道的一致性。d 分量不为零通常说明系统存在极化串扰。另一种常用分解是 Freeman-Wang 三分量分解它把像素分解为表面散射、二面角散射和体散射三种机制# 三分量分解的功率分配过程 def freeman_durand(hh, hv, vv): Freeman-Dur等三分量分解简化版 re_hh np.real(hh) re_vv np.real(vv) # 体散射项由 HV 通道功率决定 p_v 3 * np.abs(hv)**2 # 表面散射与二面角散射的功率分配 # 通过 HH 与 VV 的实部比值判断主导机制 ratio re_hh / re_vv has_surface (ratio 1.0) (ratio 1.5) p_s np.where(has_surface, np.abs(hh)**2, 0) p_d np.where(~has_surface, np.abs(hh)**2, 0) return p_s, p_d, p_v实际使用中发现在有几栋建筑的场景里体散射功率会被高估。原因是植被冠层的去极化效应在 HV 通道的功率贡献远大于建筑墙面。此时需要在像素级检查 HV 通道的信噪比若低于 10 dB则体散射分量不可信。4.3 极化特征用于图像分类的边界条件PolSAR 特征分解最常被用于土地利用分类但有一个容易被忽略的前提入射角变化会显著改变各散射机制的功率分配。同一个地物在大入射角下表现出更强的表面散射特性而在小入射角下二面角散射贡献上升。因此分类模型如果用固定阈值切分 p_s/p_d/p_v入射角差异大的图幅之间必然出现系统性误分类。处理这一问题的常见做法是估计每个像素的本地入射角代入散射模型做归一化后再分类。代码里如果直接拿 PolSAR 分解结果的 RGB 合成图去做深度学习训练务必把本地入射角作为额外特征通道输入。5. 源码调试的核心技巧从点目标到全场景逐步加复杂度的复现路径拿到这套源码后建议按点目标 → 稀疏场景 → 全场景三段路径做调试。这个方法在调试 linux 内核源码时我用过类似思路——先跑最小可执行环境再逐步打开功能模块每一步都能验证中间状态。先用单点目标验证成像算法。把场景缩小到 1 个散射点成像后检查点扩散函数是否符合理论值。然后逐步增加散射点数量每次翻一倍检查新增散射点的旁瓣是否抬升了噪声底。这个阶段最容易暴露距离徙动校正的插值精度问题——稀疏场景下旁瓣可能刚好落在采样点上干涉测量中会表现为虚假的条纹起伏需要做子像素级相位检查才能发现。这一步完成后开始加入 InSAR 的轨道规划。干涉基线 B⊥ 设置从 30 米到 150 米量级递增对应高程模糊度从 80 米变化到 16 米干涉图上能看到明显的条纹密度变化。如果源码里给出的高程反演结果与理论值之间的误差不随基线变化而稳定说明相位解缠模块的参数需要调整。最后加入 PolSAR 处理时用三面角反射器作为定标参考。反射器在 HH 与 VV 通道的理论相位差为零实测若不为零说明收发通道之间存在耦合误差需要在定标矩阵中补一组相位补偿项。这里用到的技巧是把误差从系统层面拆到信号层面用已知目标的测量值反推系统矩阵——与从源码里逆推算法流程是同一套思路。第 4 章里给出的 Freeman 分解是简化版本实际工程中还需要考虑去方向性与地形起伏校正。验证时可以把分解得到的表面散射功率与光学影像中的裸土区域做对比偏差大于 3 dB 时排查入射角归一化与定标余差。整条调试路径跑通后这套源码就不仅是算法学习材料而是一份可以迭代出可用干涉测量结果的工具链。本文还有配套的精品资源点击获取 合成孔径雷达源码解析：SAR成像、InSAR干涉与PolSAR极化分解全链路实战
简介一套覆盖合成孔径雷达三大方向的Matlab源码合集面向遥感和信号处理方向的学生与研究人员聚焦SAR成像、InSAR干涉处理及PolSAR极化定标可满足毕业设计、课程设计和算法复现需求。包内共191个文件以116个m源码为核心辅以35个mat仿真数据、21个pdf原理文档、13个txt运行说明和3个md总结笔记整体约161.84MB目录按SAR、InSAR、PolSAR模块划分并附带理论设计示意图便于对照学习。SAR部分实现RD与CS两种经典成像算法涵盖点目标仿真及Radarsat-1实测数据处理InSAR部分针对平地与圆锥两种人造场景给出从原始回波仿真、成像到干涉处理的完整流程PolSAR部分提供极化定标算法并配有项目文档与算法解析。所有源码均经过测试已有234人学习适合在此基础上扩展也适合作为相关课题的起点。1. 合成孔径雷达源码库里的三套链路到底在解决什么问题拿到一份合成孔径雷达源码包你会发现目录里通常躺着三套互不相同的代码SAR 成像、InSAR 干涉处理、PolSAR 极化分析。很多人以为这是同一个处理流程的三个模块实际上一上来就各自为政。SAR 解决的是从原始回波到一张高质量图像InSAR 解决的是用两张或多张图像反演地表形变与高程PolSAR 则直接把后向散射拆成多种散射机制——植被、裸地、建筑物、水体在不同的极化通道里呈现出完全不同的响应。这三件事共享同一个平台运动模型和信号模型但数据处理几何关系各不相同。本文把这套源码拆开讲先跑通 RD 成像算法拿到聚焦图像再做人造场景的原始回波仿真与干涉处理最后延伸到极化分解。适合要做 SAR 课程设计、星载 SAR 数据预处理、或者想把干涉测量误差源头定位到仿真阶段的工程师。下文给出可复现代码与参数边界。2. SAR成像算法从原始回波到二维图像的RD成像原理与最小实现2.1 为什么成像算法的核心卡在距离徙动校正合成孔径雷达在工作时平台一直向前运动雷达波束以正侧视或斜视方式照射地面目标。目标到天线之间的斜距随时间变化这带来两个独立问题一是方位向的多普勒历史这是 SAR 为什么能有极高方位分辨率的物理基础二是由于斜距变化导致的距离单元走动和弯曲统称为距离徙动。脉冲压缩原理使距离向分辨率做到足够高但若在距离向压缩之后直接把每一行数据按慢时间摆放目标能量会跨距离单元蔓延图像就是模糊。距离徙动校正是 RD 成像算法中决定聚焦质量的核心步骤。RD 算法通过在距离频域乘以一个与方位频率相关的校正因子等效地把弯曲的轨迹拉直这之后方位向匹配滤波才能正确压缩。成像算法的完整流程可以用三个基本步骤概括范围向压缩、距离徙动校正、方位向压缩。把每一列数据沿慢时间做一条一维匹配滤波目标能量集中于一点。常见做法是先对原始回波数据做范围向快速傅里叶变换在频域乘以匹配滤波器后变换回时域然后在距离频域完成距离徙动校正最后沿方位向做匹配滤波。这样做的好处是每一步都可以独立验证中间结果的物理意义出现问题容易定位。2.2 RD算法三步拆解与Python实现2.2.1 距离向压缩的最小实现先构造线性调频信号作为发射波形import numpy as np def range_compress(raw_data, Kr, Tr, fs): 距离向压缩频域匹配滤波 raw_data: (方位向, 距离向) 复数原始数据 Kr: 调频率 (Hz/s) Tr: 脉冲宽度 (s) fs: 距离向采样率 (Hz) nr raw_data.shape[1] # 构造匹配滤波器参考函数 t np.arange(-nr//2, nr//2) / fs ref np.exp(-1j * np.pi * Kr * t**2) ref np.conj(ref) # 匹配滤波 参考信号的共轭反转 ref_fft np.fft.fft(ref, nr) range_fft np.fft.fft(raw_data, axis1) compressed np.fft.ifft(range_fft * ref_fft, axis1) return compressed这段代码里ref的构建是重点发射信号是正调频率的线性调频脉冲匹配滤波器是它的时间反转共轭。在频域里这个操作等价于对所有频率分量施加与信号相位共轭的相位变化使得目标能量集中在零距离单元附近。注意这里的t是从-nr/2到nr/2因为匹配滤波器的群延迟天然会使峰值出现在参考时间位置。如果你发现压缩后的峰值位置整体偏移多半是这里的时间轴定义与发射信号的时间基准不匹配。2.2.2 距离徙动校正的插值实现距离压缩完成后每一个方位慢时间时刻的回波峰值能量实际上落在不同的距离单元里。对这个量做二维数组上的插值搬移就是距离徙动校正的常见做法。def rcm_correction(compressed, kr, ka, fdc, v_r, lambda_c, nr): 距离徙动校正距离频域插值实现 compressed: 距离压缩后的二维数组 v_r: 等效雷达速度 (m/s) fdc: 多普勒中心频率 (Hz) n_az, n_rg compressed.shape f_range np.fft.fftfreq(n_rg, d1/fs) # 距离频率 # 先沿方位向做FFT进入距离多普勒域 rd_domain np.fft.fft(compressed, axis0) f_az np.fft.fftfreq(n_az, d1/prf) # 方位频率 for i, fa in enumerate(f_az): # 计算目标在距离多普勒域中的斜距 R_fa lambda_c * (fdc - fa) / (2 * v_r**2) * (fdc - fa) / (4 * (fdc - fa)**2 / lambda_c**2) # 实际工程中用的是近似式 R_fa lambda_c * (fdc - fa) / (4*ka) # 距离偏移量以距离采样单元为单位 delta_r R_fa - np.mean(R_fa) shift_bins delta_r * 2 * fs / c # 在距离频域乘以线性相位实现时域搬移 phase np.exp(1j * 2 * np.pi * f_range * shift_bins / fs) rd_domain[i, :] np.fft.ifft(np.fft.fft(rd_domain[i, :]) * phase) # 反变换回二维时域 corrected np.fft.ifft(rd_domain, axis0) return corrected这段代码的关键点是shift_bins的计算距离频域乘以线性相位exp(j2π f_range shift_bins / fs)等价于时域里把信号搬移shift_bins个采样点。这样一来不同方位频率分量在各目的距离频域中独立完成偏移不需要做逐像素的插值循环。实际参数下shift_bins通常在零点几到几十个距离单元之间。平台轨道越高、斜距越远徙动量就越大不做校正时目标能量会倾斜分布在多个距离单元中导致方位向压缩后的旁瓣不对称。2.2.3 方位向压缩距离徙动校正完成后信号在距离多普勒域里已经是一条直线。方位向压缩与距离向压缩形式完全相同只是参考函数换成方位向的多普勒调频率def azimuth_compress(corrected, prf, ka, t_az): 方位向压缩 n_az corrected.shape[0] # 方位向参考函数线性调频 ref_az np.exp(-1j * np.pi * ka * np.linspace(-t_az/2, t_az/2, n_az)**2) ref_az np.conj(ref_az) # 逐距离单元做方位向匹配滤波 az_fft np.fft.fft(ref_az, n_az) result np.zeros_like(corrected) for r in range(corrected.shape[1]): row_fft np.fft.fft(corrected[:, r]) result[:, r] np.fft.ifft(row_fft * az_fft) return resultka代表方位向多普勒调频率直接由平台速度与斜距决定。星载 SAR 的ka值远大于机载意味着方位向积累时间更短所需的脉冲数更少但相位误差对运动的敏感度也更高。2.3 成像后的质量检查点扩散函数与分辨率验证用点目标测试合成孔径雷达成像算法是否实现正确需要通过点扩散函数的 3 个指标验证峰值旁瓣比、积分旁瓣比、分辨率。下表是通用检查标准指标理想值实际可接受范围说明方位向峰值旁瓣比-13.26 dB优于 -12 dB过低的旁瓣意味着加窗这会同时展宽主瓣距离向峰值旁瓣比-13.26 dB优于 -12 dB未加窗时与理论值接近方位向分辨率λR/(2L)理论值的 1.0~1.4 倍L 为合成孔径长度加窗后会展宽距离向分辨率c/(2B)理论值的 1.0~1.4 倍B 为发射信号带宽成像完成后找孤立点目标取峰值周围 9×9 像素区域做插值搜索主瓣峰值位置量测 -3 dB 处的宽度。如果测出的旁瓣比明显低于 -13 dB 且主瓣不对称优先怀疑距离徙动校正里shift_bins的计算公式是否有符号错误。3. InSAR人造场景回波仿真与干涉处理从场景建模到干涉图完整闭环3.1 人造场景的几何精度可控性优势干涉合成孔径雷达测量InSAR需要从两幅复图像中提取相位差从而反演高程或微小形变。实测数据面临的难点在于真实地形精确已知度有限、轨道误差与大气延迟会混入干涉相位。人造场景的好处是地面目标的三维坐标、后向散射系数、形变场完全由脚本设置仿真回波经过成像后得到的干涉相位可以与初始模型严格比对。这样就把测量链路的误差源头从环境误差中分离出来。在做源码级实验时我一般会用程序生成一个包含已知高程梯度的数字高程模型DEM如斜坡、台阶、凸起三种形态再将 DEM 叠加到基准平面。这种方式能快速区分成像算法误差与干涉处理误差。3.2 人造场景的原始回波仿真与平台参数表原始回波仿真的基本思路是把场景离散成散射点对每个散射点计算其到雷达的瞬时斜距生成对应的回波并叠加到原始数据矩阵中。接收信号模型可以写成s_echo(t, η) Σ A_i · rect((t - 2R_i(η)/c) / Tr) · exp(-j·4π·R_i(η)/λ) · w_a(η - η_c,i) n(t, η)其中 R_i(η) 是第 i 个散射点在方位慢时间 η 时的斜距A_i 包含后向散射系数与天线方向图加权n(t, η) 是热噪声项。# 点阵场景回波仿真 Python 骨架 def simulate_raw_echo(scene_points, az_time, range_time, params): scene_points: [(x, y, z, sigma), ...] 地面点坐标 az_time: 方位向慢时间向量 (s) range_time: 距离向快时间向量 (s) c 3e8 raw np.zeros((len(az_time), len(range_time)), dtypecomplex) for (x0, y0, z0, sigma) in scene_points: # 正侧视几何雷达位于 (0, v_r * eta, H) # 目标在雷达坐标中的斜距 R0 np.sqrt((params[H] - z0)**2 x0**2 y0**2) for i, eta in enumerate(az_time): # 平台位置 R np.sqrt((params[H] - z0)**2 x0**2 (y0 - v_r*eta)**2) # 距离向回波位置 tau 2*R/c idx np.argmin(np.abs(range_time - tau)) if abs(range_time[idx] - tau) 1.0/fs: # 线性调频信号回波 phase -4*np.pi*R/params[lambda] raw[i, idx] sigma * np.exp(1j*phase) return raw这段仿真代码的关键参数需要显式设置参数符号推荐值作用平台速度v_r150 m/s方位向积累长度与多普勒调频率控制平台高度H3000 m斜距范围控制载频fc9.6 GHz波长 → 干涉相位对高程的敏感度脉冲宽度Tr0.5 μs距离向分辨率的一阶控制距离向带宽B150 MHz距离向分辨率脉冲重复频率PRF300 Hz方位向采样率地面点数N500~2000点数过多仿真会慢一个数量级平台速度与 PRF 的组合关系必须满足方位向采样定理PRF 必须大于方位向多普勒带宽的两倍否则会因欠采样出现方位模糊。仿真完成后回波数据直接送入第 2 章的 RD 成像算法得到两幅单视复图像。人工场景不需要额外做地形相位补偿这一步这正是仿真数据的便利之处。3.3 干涉处理最小流程配准、干涉图、滤波与解缠3.3.1 主辅图像亚像素配准干涉测量要求两幅复图像在像素级对齐否则相干性显著下降。常用做法是先通过星历参数估算一个粗配准偏移量再做基于相关函数的亚像素配准def cross_correlation_register(master, slave, search_range8): 基于幅度互相关的亚像素配准单点计算示例 # 在频域计算移位后的互相关 m_fft np.fft.fft2(master) s_fft np.fft.fft2(slave) # 幅度归一化互相关 cross m_fft * np.conj(s_fft) cross_map np.abs(np.fft.ifft2(cross / np.abs(cross))) max_idx np.unravel_index(np.argmax(cross_map), cross_map.shape) shift (max_idx[0] - cross_map.shape[0]//2, max_idx[1] - cross_map.shape[1]//2) return shift距离向配准误差在 0.1 像素以内方位向在 0.1 像素以内才能保证相干性损失在可接受范围内。实际操作中如果配准偏移量在方位向超过 1 个像素说明两幅数据的零多普勒时刻偏移没算准这时候要先检查轨道插值参数。3.3.2 多视处理与相位解缠配准后按水平和垂直方向的视数对干涉图做多视平滑。对于地形平坦的仿真场景干涉图主要体现目标高程引起的平缓相位变化此时应采用加权圆均值滤波窗口大小 3×3 与 5×5 两种def multilook_filter(interferogram, look_range3, look_az3): 多视滤波 from scipy.ndimage import uniform_filter phase np.angle(interferogram) mag np.abs(interferogram) # 拆成相互独立的实部/虚部分别做窗口均值 re uniform_filter(mag*np.cos(phase), size(look_az, look_range)) im uniform_filter(mag*np.sin(phase), size(look_az, look_range)) return re 1j*im相位解缠是干涉处理中最容易出问题的一环。不同解缠算法的选择直接影响高程反演精度解缠算法适用场景已知局限工程推荐参数枝切法低噪声、相干性高、残差点稀少残差点密集时出现大面积非连通区域枝切长度上限32 像素最小二乘残差点较多、区域连续高程跳变处相位被平滑加权系数设为相干性区域增长法不规则区域边界初始种子点选择敏感种子阈值0.5SNAPHU统计费用流形变场复杂、高噪声耗时显著需调参数默认参数即可解缠后的相位通过 ∮ -(4πB⊥)/(λR sinθ) · h 反演高程。在卫星轨道安全、基线估计精确的仿真条件下如果解缠相位与真值 DEM 的残差仍然呈条纹状分布说明基线参数与成像几何参数不一致需要回到平台参数设置中检查基线相位贡献。4. PolSAR极化定标与特征分解四通道不是四个单通道4.1 合成孔径雷达极化测量为什么单独占一套源码SAR 与 InSAR 处理的对象是标量复数数据PolSAR 处理的对象是每个像素一个 2×2 复数散射矩阵包含 HH、HV、VH、VV 四个极化通道。这四个通道之间不是独立的电磁波在散射过程中发生极化状态变化矩阵各元素之间存在明确的相位关系。直接忽略交叉极化通道、只保留同极化幅度信息会丢失地物分类中的关键物理线索。PolSAR 源码要解决的第一个问题是极化定标。系统硬件在发射和接收路径上会引入通道不平衡和交叉串扰如果不校准后续分解得到的散射功率分布完全不可信。4.2 极化定标矩阵与最小验证集雷达测量得到的极化散射矩阵 M 与真实散射矩阵 S 的关系为M A·R·S·T N其中 T 与 R 分别是发射和接收失真矩阵。做定标时常用一个二面角反射器或三面角反射器其理论散射矩阵已知。定标后计算剩余误差如幅度不平衡优于 0.3 dB串扰优于 -30 dB。def pauli_decompose(hh, hv, vh, vv): Pauli基分解把散射矩阵投影到Pauli基上 Pauli基奇次散射、偶次散射、45度偶次散射、螺旋散射 # 对应 S a[1 0;0 1] b[1 0;0 -1] c[0 1;1 0] d[0 -1j;1j 0] a (hh vv) / np.sqrt(2) # 奇次散射球体、三面角 b (hh - vv) / np.sqrt(2) # 偶次散射二面角 c (hv vh) / np.sqrt(2) # 45°旋转的偶次散射 d (hv - vh) / np.sqrt(2) # 螺旋散射分量 # 散射功率 p_a np.abs(a)**2 p_b np.abs(b)**2 p_c np.abs(c)**2 return p_a, p_b, p_c, np.ones_like(p_a)*np.abs(d)**2对于一个三面角反射器理论上 p_b、p_c 应为零实际中如果 p_b 与 p_a 比值高于 -20 dB需要回去排查接收通道的一致性。d 分量不为零通常说明系统存在极化串扰。另一种常用分解是 Freeman-Wang 三分量分解它把像素分解为表面散射、二面角散射和体散射三种机制# 三分量分解的功率分配过程 def freeman_durand(hh, hv, vv): Freeman-Dur等三分量分解简化版 re_hh np.real(hh) re_vv np.real(vv) # 体散射项由 HV 通道功率决定 p_v 3 * np.abs(hv)**2 # 表面散射与二面角散射的功率分配 # 通过 HH 与 VV 的实部比值判断主导机制 ratio re_hh / re_vv has_surface (ratio 1.0) (ratio 1.5) p_s np.where(has_surface, np.abs(hh)**2, 0) p_d np.where(~has_surface, np.abs(hh)**2, 0) return p_s, p_d, p_v实际使用中发现在有几栋建筑的场景里体散射功率会被高估。原因是植被冠层的去极化效应在 HV 通道的功率贡献远大于建筑墙面。此时需要在像素级检查 HV 通道的信噪比若低于 10 dB则体散射分量不可信。4.3 极化特征用于图像分类的边界条件PolSAR 特征分解最常被用于土地利用分类但有一个容易被忽略的前提入射角变化会显著改变各散射机制的功率分配。同一个地物在大入射角下表现出更强的表面散射特性而在小入射角下二面角散射贡献上升。因此分类模型如果用固定阈值切分 p_s/p_d/p_v入射角差异大的图幅之间必然出现系统性误分类。处理这一问题的常见做法是估计每个像素的本地入射角代入散射模型做归一化后再分类。代码里如果直接拿 PolSAR 分解结果的 RGB 合成图去做深度学习训练务必把本地入射角作为额外特征通道输入。5. 源码调试的核心技巧从点目标到全场景逐步加复杂度的复现路径拿到这套源码后建议按点目标 → 稀疏场景 → 全场景三段路径做调试。这个方法在调试 linux 内核源码时我用过类似思路——先跑最小可执行环境再逐步打开功能模块每一步都能验证中间状态。先用单点目标验证成像算法。把场景缩小到 1 个散射点成像后检查点扩散函数是否符合理论值。然后逐步增加散射点数量每次翻一倍检查新增散射点的旁瓣是否抬升了噪声底。这个阶段最容易暴露距离徙动校正的插值精度问题——稀疏场景下旁瓣可能刚好落在采样点上干涉测量中会表现为虚假的条纹起伏需要做子像素级相位检查才能发现。这一步完成后开始加入 InSAR 的轨道规划。干涉基线 B⊥ 设置从 30 米到 150 米量级递增对应高程模糊度从 80 米变化到 16 米干涉图上能看到明显的条纹密度变化。如果源码里给出的高程反演结果与理论值之间的误差不随基线变化而稳定说明相位解缠模块的参数需要调整。最后加入 PolSAR 处理时用三面角反射器作为定标参考。反射器在 HH 与 VV 通道的理论相位差为零实测若不为零说明收发通道之间存在耦合误差需要在定标矩阵中补一组相位补偿项。这里用到的技巧是把误差从系统层面拆到信号层面用已知目标的测量值反推系统矩阵——与从源码里逆推算法流程是同一套思路。第 4 章里给出的 Freeman 分解是简化版本实际工程中还需要考虑去方向性与地形起伏校正。验证时可以把分解得到的表面散射功率与光学影像中的裸土区域做对比偏差大于 3 dB 时排查入射角归一化与定标余差。整条调试路径跑通后这套源码就不仅是算法学习材料而是一份可以迭代出可用干涉测量结果的工具链。本文还有配套的精品资源点击获取