这是 Keystone 变换系列的第四篇。前面几篇我把公式推导和物理图像都讲了一遍重点解释了为什么运动目标的回波会“跑”出距离单元以及 Keystone 变换为什么能通过重采样把慢时间轴“掰弯”来校正距离走动。这篇直接落地用 MATLAB 把整个流程跑通。我会给出一组完整可复现的雷达仿真参数从回波生成、脉冲压缩到 Keystone 变换、相参积累全部用代码说话。内容适合做雷达信号处理、在准备相参积累相关的课程设计或者毕业设计的同学也适合已经从公式层面理解了 Keystone但还不知道“频域重采样”在 MATLAB 里到底怎么写的人。1. 仿真之前必须想明白的三个问题1.1 距离走动在二维图上到底长什么样仿真之前先建立直觉。雷达发射一串脉冲每个脉冲的快时间维对应距离脉冲之间的慢时间维对应观测时刻。目标如果静止不动每个脉冲回波经过匹配滤波之后峰值包络都落在完全相同的距离单元上画成“距离-慢时间”二维图就是一条水平的亮线。目标一旦运动每个脉冲对应的时间延迟 (2R(t_m)/c) 在变包络峰值就会在一个一个距离单元之间移动二维图上看起来就是一条倾斜的线。这条斜线就是距离走动。走动的速度直接和目标径向速度挂钩走得越快积累时间内跨过的距离单元就越多。仿真里最容易犯的错是只盯着波形看而忽略距离单元的定义。距离单元宽度是 (c/(2B))带宽 100 MHz 时也就是 1.5 米而一个高速目标在 0.1 秒内可能飞出几十米。换句话说如果不做任何校正直接做慢时间维傅里叶变换去测速度目标能量会被摊到十几个距离单元上多普勒谱严重展宽积累增益大幅下降。1.2 Keystone变换的输入输出到底是谁很多文章把 Keystone 变换讲得很玄其实落到仿真里就是一句话输入端是“脉冲压缩后的二维矩阵”输出端是“重采样后的二维矩阵”。矩阵的两个维度分别是慢时间和距离频率。注意这里说的距离频率是指对快时间维做傅里叶变换之后的频率轴不是慢时间维的多普勒频率。为什么要到距离频域去做因为距离走动在时域表现为包络位置的平移而在频域表现为相位项 ( -4\pi (f_c f_\tau) v t_m / c )。这个相位项是距离频率 (f_\tau) 和慢时间 (t_m) 的乘积耦合所以没法用简单的时移补偿。Keystone 的思路很直接对每个距离频率点把慢时间轴按比例 (f_c/(f_cf_\tau)) 重新缩放让耦合项中的 (f_\tau) 和 (t_m) 不再相乘距离走动自然就被拉回来了。1.3 为什么不能在时域直接“对齐峰值”容易想到的一个替代方案是既然知道包络是斜线就估计斜率然后把每个脉冲平移回来。这种做法叫包络对齐在 ISAR 成像里很常用但问题在于它需要先估计目标速度或者包络斜率而且对多个不同速度的目标处理起来很麻烦。Keystone 变换的好处是完全不需要速度先验它用重采样操作把“速度-频率”耦合统一消除不管你是一个目标还是多个目标只要信噪比不是太差都能一次性拉直。2. 仿真参数与信号构造2.1 一组留有余量的参数表仿真参数直接决定了距离走动是否明显也决定了 Keystone 变换能不能体现出效果。我用的这组参数比较典型覆盖了常见的 X 波段雷达场景。参数数值说明载频 (f_c)10 GHzX 波段带宽 (B)100 MHz距离分辨率 1.5 m脉宽 (T_p)10 μs匹配滤波增益的基础采样率 (f_s)100 MHz取等于带宽距离窗无过采样脉冲重复频率 PRF2000 Hz无多普勒模糊对应的最大速度约 30 m/s后面细讲积累脉冲数 M256积累时间 128 ms初始距离 R010000 m远场目标速度 v200 m/s高速明显距离走动计算一下积累时间 (T_{obs} M / PRF 128) ms目标在这段时间里移动约 25.6 m。距离单元宽度只有 1.5 m所以包络会跨越约 17 个距离单元。如果这条路走通了Keystone 变换带来的效果在图上会非常直观。2.2 生成回波信号的MATLAB代码下面这段代码生成基带回波。为了教学直观我用了循环逐脉冲生成实际工程中可以改成矩阵化运算提速。%% 参数设置 c 3e8; fc 10e9; B 100e6; Tp 10e-6; K B / Tp; PRF 2000; M 256; dt 1 / PRF; fs 100e6; Ts 1 / fs; Nfast 1024; R0 10000; v 200; tm (0:M-1) * dt; R R0 - v * tm; % 目标距离随慢时间线性变化 %% 快时间轴 t_ref 2 * R0 / c; % 参考时延 tau t_ref - Tp/2 (0:Nfast-1) * Ts; %% 发射基带信号 tx exp(1j * pi * K * ((tau - t_ref).^2)); % 以参考时延为中心的内向线性调频 tx(abs(tau - t_ref) Tp/2) 0; %% 回波生成 s zeros(M, Nfast); for m 1:M td 2 * R(m) / c; t tau - td; s(m, :) exp(1j * pi * K * t.^2) .* (abs(t) Tp/2); end s s .* exp(-1j * 4 * pi * fc * R(:) / c); % 慢时间相位速度信息所在这段代码里最关键的是那个慢时间相位 (\exp(-j4\pi f_c R(t_m)/c))。它包含了目标运动引起的多普勒信息后面测速度就要靠它。2.3 脉冲压缩与运动补偿前的基线结果脉冲压缩用频域匹配滤波实现。构造匹配滤波器的频响时要注意把发射信号补零到和回波矩阵相同的列数再取共轭。Nfft 2^nextpow2(Nfast * 2); S fft(s, Nfft, 2); ref exp(1j * pi * K * ((tau - t_ref).^2)); % 参考信号 H fft(ref, Nfft, 2); Y ifft(S .* conj(H), [], 2); Y Y(:, 1:Nfast);这里参考信号取以参考时延为中心是为了让匹配滤波后的峰值出现在参考距离附近。如果不做任何走动校正直接查看二维图的幅度会看到目标峰值包络随慢时间明显偏移一条斜线贯穿多个距离门。这就是我们后面要校正的对象。3. 核心实现MATLAB里的三种Keystone变换写法3.1 方法一慢时间维Sinc插值最符合原理、也最好理解的实现方式是 Sinc 插值。对每个距离频率点 (f_\tau)计算缩放因子 (\alpha f_c / (f_c f_\tau))然后对慢时间序列在非均匀的 (t_m \alpha t_m) 处重新采样。Yf fftshift(fft(Y, Nfft, 2), 2); freq ((-Nfft/2 : Nfft/2-1) / Nfft) * fs; ValidF abs(freq) B/2; Ykt zeros(size(Yf)); for i 1:Nfft if ValidF(i) alpha fc / (fc freq(i)); tm_new tm * alpha; Ykt(:, i) ks_sinc_interp(Yf(:, i), tm, tm_new, dt, 16); else Ykt(:, i) Yf(:, i); end end y_kt ifft(ifftshift(Ykt, 2), [], 2); y_kt y_kt(:, 1:Nfast);自定义 Sinc 插值函数如下。这里的思路是对每个目标采样点用带限 Sinc 核卷积估计新时刻的值核长截断到 16 个采样点。function y ks_sinc_interp(x, t, t_new, dt, Ntap) % x: Mx1 原始慢时间序列 % t: Mx1 原始慢时间 % t_new: Lx1 新慢时间 % dt: 慢时间采样间隔 % Ntap: 核截断长度 L length(t_new); M length(x); y zeros(L, 1); for n 1:L tau (t_new(n) - t) / dt; idx find(abs(tau) Ntap); if ~isempty(idx) y(n) sum(x(idx) .* sinc(tau(idx))); end end end这种写法最贴近公式适合讲原理和验证算法。缺点是慢时间维每个距离频率都需要做一次插值256 个脉冲、1024 个距离采样点就要做 1024 次循环速度比较慢。实际仿真可以先跑通再考虑优化。3.2 方法二Chirp-Z变换实现工程上更推荐的实现是利用 Chirp-Z 变换。Keystone 重采样本质上是一种尺度变换而尺度变换可以由“频域相位调制 时域卷积”来实现正好是 Chirp-Z 的典型应用。这样做的好处是不需要显式地构造新时刻网格也不会有 Sinc 核截断造成的边缘误差计算速度更快。MATLAB 自带czt函数可以对慢时间序列做任意起止频率和点数的 Z 变换。对每个距离频率点把重采样看成在频域做一次变尺度 Z 变换即可。核心代码框架如下for i 1:Nfft if ValidF(i) alpha fc / (fc freq(i)); % 对慢时间序列做 Chirp-Z 变换等效于重采样到 alpha 倍 temp czt(Yf(:, i), M, exp(-1j * 2 * pi * alpha / M), 1); Ykt(:, i) temp(:); end end这里我只是给出最基本的思路。实际使用时要仔细确认czt的输入参数定义尤其是起点和螺旋因子的设置否则很容易把尺度关系搞反。由于czt在信号处理工具箱中可用大部分机器上都能直接跑速度比逐点 Sinc 插值快一个数量级。3.3 方法三基于FFT的等效实现第三种思路是利用 FFT 在频域完成插值。对慢时间序列先做 FFT 到多普勒域通过补零和相位修正来近似实现分数阶重采样。这个方法实现起来最快但精度受补零方式和尺度因子整数化影响较大适合快速验证或者对精度要求不高的场景。for i 1:Nfft if ValidF(i) alpha fc / (fc freq(i)); Nout round(M * alpha); spec fft(Yf(:, i), M); spec_pad zeros(Nout, 1); if Nout M spec_pad(1:M/2) spec(1:M/2); spec_pad(end-M/21:end) spec(M/21:end); else spec_pad(1:Nout/2) spec(1:Nout/2); spec_pad(end-Nout/21:end) spec(M-Nout/21:end); end Ykt(:, i) ifft(spec_pad, Nout); % 长度会变化 end end注意这个方法会导致每个距离频率上的输出长度不一致后续需要对齐到统一的慢时间网格。实际工程中通常不建议直接这样用但用来理解“重采样即变尺度频谱搬移”的物理意义很不错。三种方法的取舍我整理成了表格。方法精度速度实现复杂度场景建议Sinc插值高慢低原理验证、教学Chirp-Z高快中工程部署、批量处理FFT近似中最快中快速预览、实时性要求高4. 仿真结果校正前与校正后的对比4.1 B-scan图与距离包络对比把脉冲压缩后的幅度沿慢时间画出来校正前的 B-scan 图上目标是一条倾斜的亮线跨过的距离单元数大约为 (2vT_{obs}/c / (1/(2B)) \approx 17) 个。做完 Keystone 变换之后再画斜线变成接近水平的直线目标能量被“收拢”到同一距离单元。这里有个实操细节画图时最好把显示范围固定在同一动态范围内比如统一用分贝值20*log10(abs(Y))否则人眼很容易被颜色映射欺骗觉得“好像没什么变化”。另外校正后如果 Sinc 插值用了 0 外插图像边缘可能会出现暗区这是正常的不代表算法失效。4.2 慢时间FFT与速度测量距离走动校正的最终目的是为了相参积累。校正前对每个距离门做慢时间 FFT目标能量分散在多个距离门多普勒谱被展宽。校正后目标集中在同一个距离门对这个距离门做慢时间 FFT 会得到一个尖锐的峰值峰值频率对应的多普勒频率 (f_d -2vf_c/c)。具体到我这组参数(f_d -2 \times 200 \times 10^{10} / 3\times 10^8 \approx -13333) Hz。PRF 是 2000 Hz所以这个频率已经远远超出了 ([-PRF/2, PRF/2]) 的范围出现多普勒模糊。测出来的模糊频率是 (13333 \mod 2000 1333) Hz对应速度约 (-20) m/s。这说明仿真里 PRF 选择并不适合这么高的速度真实工程中要么提高 PRF要么用多普勒模糊数解算。4.3 信噪比改善和计算代价评价 Keystone 效果不能只看图更要看量化指标。简单做法是在校正前的二维矩阵里找到目标所在距离门和附近距离门的所有慢时间 FFT 幅度计算峰值幅度和噪声底之间的比值然后用同样方法计算校正后同一距离门上的比值。参考仿真条件下校正后峰值幅度相比校正前通常能提升 10 dB 以上具体数字取决于距离走动跨越的单元数。代价方面用 Sinc 插值跑 256×1024 的数据量在普通笔记本上可能要几秒到十几秒而用 Chirp-Z 可以压到 1 秒以内。做科研时建议先用 Sinc 验证正确性再换成 Chirp-Z 跑蒙特卡洛实验。5. 仿真过程中最容易踩的坑5.1 插值核截断导致边缘失真Sinc 插值理论上需要无限长的核但实际只能截断。截断会带来两个问题一是回波边缘产生振铃二是当缩放因子偏离 1 较多时新时刻可能落在原始慢时间范围之外直接外插会把虚假能量带进来。我的处理办法是在插值函数里对超出范围的值强制置 0只保留中间约 80% 的慢时间数据进行后续测速和多普勒分析。5.2 快时间频率轴偏移导致重采样比例错误频域做 Keystone 时距离频率轴必须对应基带频率也就是范围是 ([-f_s/2, f_s/2])。很多人直接用fft后的下标当频率忘记做fftshift结果频率正负颠倒重采样比例全部反了效果自然不对。建议在生成频率轴之后先用一个直流信号自检确保频率为 0 的点位于Nfft/21的位置。5.3 多普勒模糊对测速结果的影响前面已经提到PRF 太低时 Keystone 变换本身仍然可以校正距离走动但测速会模糊因为慢时间采样不满足奈奎斯特条件。遇到这种情况可以先盲速分割或者用多 PRF 解模糊。如果只是在仿真里验证 Keystone 本身建议把目标速度降到 30 m/s 以内这样 PRF 为 2000 Hz 时不会模糊结果看起来更干净。5.4 数据量太大时如何加速Keystone 最耗时的部分在慢时间维插值。除了 Chirp-Z 之外还可以用 MATLAB 的interp1并行计算。如果处理的是采信机采回来的数据可能还涉及 ADC 量化、十六进制转有符号数这些环节建议先把数据切成长度合适的块再逐块做 Keystone。实在需要极致性能就把插值核写成 C 代码用 MEX 编译调进 MATLAB这样比纯脚本快很多。个人体会是Keystone 变换在 MATLAB 里跑通不难难的是把边界条件和频率轴关系理清楚。刚开始做仿真时我第一版代码跑出来的结果是一条完全乱掉的曲线检查了很久才发现是频率轴没有做fftshift。所以建议你拿到代码后先跑一组低速目标确认基线正确再加大速度这样能快速定位问题。这个系列后面的内容我打算继续深入聊 Keystone 在多目标场景下的表现以及和长时间相参积累算法的配合如果大家在复现过程中遇到其他锯齿欢迎一起交流。