常用合成孔径雷达（SAR）成像算法

核心任务：把原始回波数据“聚焦”成一幅清晰图像

算法	主要优点	主要缺点	更适合的场景
RD（Range-Doppler）	结构最经典，流程清晰；基于 FFT，计算效率较高；工程实现成熟，教材和传统系统里最常见	需要 RCMC 插值，实现和硬件流水不如 CS 规整；在高斜视、宽孔径、强距离-方位耦合时精度和适应性会下降	常规条带式 SAR、几何条件较温和、需要成熟稳定实现的系统
BP（Back Projection）	几何适应性最强；对非理想轨迹、复杂平台运动、近场/大视角、更复杂成像几何更鲁棒；成像物理意义直接，精度潜力高	计算量通常最大，往往远高于频域算法；大幅面、高分辨率、实时成像时硬件压力明显	UAV/机载、轨迹抖动明显、复杂地形、非直线轨迹、需要更强通用性的场景
CS（Chirp Scaling）	用相位补偿代替 RD 的 RCMC 插值；结构通常由 FFT + 相位补偿组成，硬件实现友好；效率高，适合流水化和 FPGA/ASIC	本质仍是频域近似法；在高斜视、超宽孔径、强高阶耦合条件下，基础 CS 也会遇到补偿精度不足的问题	想保留频域高效率、又希望避免大量 RCMC 插值的工程实现
ω-k（Omega-K）	在二维频域统一处理耦合；通常借助 Stolt 映射处理距离徙动和聚焦问题；对宽孔径、高斜视、大带宽场景通常比 RD/CS 更稳健	实现理解门槛更高；需要 Stolt 重采样/插值，参数和边界处理更敏感；工程调试复杂度通常高于 RD/CS	高斜视、高分辨率、宽孔径、距离-方位耦合显著的场景

工程选用原则：

先做系统原型、常规条带式场景：优先看 RD。
想保留频域高效率，同时减少 RCMC 插值负担：选 CS。
高斜视、宽孔径、耦合明显，追求更强精度：优先看 ω-k。
平台轨迹不规则、地形复杂、近场/机载/UAV、需要最高几何通用性：选 BP。

Range-Doppler算法

思路：

先在 距离向 做匹配滤波，得到距离压缩。
再到 方位多普勒域 去处理方位向聚焦。
因为同一个目标在合成孔径期间会跨越多个距离单元，所以要做 RCMC（range cell migration correction，距离单元徙动校正）。
最后做方位压缩，得到图像。

1. 成像几何与点目标原始回波

设平台沿方位向以恒速 v 飞行，目标位于参考斜距 R_0 处；令慢时间（方位时间）为 \eta，快时间（距离时间）为 \tau。在零斜视 stripmap 条件下，目标的瞬时斜距为

R(\eta)=\sqrt{R_0^2+(v\eta)^2}

这就是 SAR 点目标的“距离历程”。在目标通过波束中心附近，若 R_0 \gg v|\eta|，可作二阶展开

R(\eta)\approx R_0+\frac{v^2\eta^2}{2R_0}

这个几何模型以及二阶近似，是后面把方位回波看成线性调频信号的基础。Moreira 等的教程明确给出了这一斜距表达式与二阶展开，并指出更精确处理应使用完整相位历程。

发射线性调频脉冲后，一个点目标的接收基带回波可写成

s_{\text{rb}}(\tau,\eta) = A_0\,w_a(\eta)\, w_r\!\left(\tau-\frac{2R(\eta)}{c}\right)\, \exp\!\left[-j\frac{4\pi f_c}{c}R(\eta)\right]\, \exp\!\left\{j\pi K_r\left[\tau-\frac{2R(\eta)}{c}\right]^2\right\}

其中 f_c 为载频，K_r 为距离向调频率，w_r 是距离包络，w_a 是天线方位包络，A_0 含目标散射系数和常数幅度项。这个形式正是 stripmap SAR 点目标原始回波的标准基带模型。

2. 第一步：距离压缩

对每条方位线，沿快时间 \tau 用发射 LFM 的共轭做匹配滤波：

h_r(\tau)=\exp(-j\pi K_r\tau^2)

频域等价为乘以其频响共轭。距离压缩后的结果，在忽略窗函数与常数因子时，可写为

s_{\text{rc}}(\tau,\eta) \approx A_0\,w_a(\eta)\, \mathrm{sinc}\!\left[B_r\left(\tau-\frac{2R(\eta)}{c}\right)\right] \exp\!\left[-j\frac{4\pi}{\lambda}R(\eta)\right]

这里 \lambda=c/f_c，B_r 是发射带宽。

这一步的物理意义是：把快时间上的 LFM 脉冲压缩成一条窄脉冲，但脉冲中心仍随 \eta 按 2R(\eta)/c 变化。也就是说，距离压缩后，目标能量沿着一条弯曲轨迹分布，这就是后面要校正的 RCM（range cell migration）。教程对 SAR 处理的第一步“距离压缩”和 RCM 的来源有直接说明。

3. 距离压缩后为什么方位向还是“啁啾”

如果固定看目标主瓣附近，即取 \tau \approx 2R(\eta)/c，距离压缩后的相位主项只剩

\phi_a(\eta)=-\frac{4\pi}{\lambda}R(\eta)

代入二阶展开

R(\eta)\approx R_0+\frac{v^2\eta^2}{2R_0}

得到

\phi_a(\eta)\approx -\frac{4\pi R_0}{\lambda} -\frac{2\pi v^2}{\lambda R_0}\eta^2

于是

s_a(\eta)\propto w_a(\eta)\exp\!\left(j\pi K_a\eta^2\right), \qquad K_a=-\frac{2v^2}{\lambda R_0}

所以，方位回波本质上也是 LFM（azimuth chirp）。Moreira 的教程从相位导数出发，给出了方位频率（即 Doppler 频率）随慢时间近似线性变化的结果，说明了为什么方位压缩本质上也是匹配滤波问题。

进一步，瞬时多普勒频率为

f_\eta(\eta)=\frac{1}{2\pi}\frac{d\phi_a}{d\eta} \approx -\frac{2v^2}{\lambda R_0}\eta =K_a\eta

这说明慢时间 \eta 与方位频率 f_\eta 一一对应；做方位 FFT 后，目标会在 Doppler 域被“展开”，而这一域正是 RD 算法工作的主舞台。

4. 从二维点目标谱看 RD 的本质

现在把距离压缩后的信号再对方位时间 \eta 做 Fourier 变换，进入 (\tau,f_\eta) 的 Range–Doppler 域。严格推导这一步通常用 驻相法；对零斜视单基地 stripmap 点目标，其结果可写成“一个随 f_\eta 变化的距离延迟 + 一个纯方位相位项”的结构。工程上常写成与 RDA 滤波器一致的形式：

S_{\text{RD}}(f_\tau,f_\eta;R_0) \propto W_r(f_\tau)\,W_a(f_\eta)\, \exp\!\left[-j\frac{4\pi R_0}{cD(f_\eta)}\,f_\tau\right]\, \exp\!\left[-j\frac{4\pi R_0}{\lambda}D(f_\eta)\right]

其中

D(f_\eta)=\sqrt{1-\left(\frac{\lambda f_\eta}{2v}\right)^2}

这正是传统 RDA 里方位压缩滤波器和 RCMC 滤波器出现的原因：

第二个指数项只含 f_\eta，它是 方位匹配滤波 对应的相位。
第一个指数项对 f_\tau 线性，意味着 目标在 RD 域中的等效距离延迟依赖于 f_\eta，这就是 距离徙动在 Doppler 域的表现。

Chen 与 Kiang 给出的传统 RDA 滤波器形式，恰好对应这个二维点目标谱的分解：

H_{\text{ac}}(f_\tau,f_\eta)=\exp\!\left(j\frac{4\pi R_0}{\lambda}D(f_\eta)\right)

以及零斜视时

H_{\text{rcmc}}(f_\tau,f_\eta)= \exp\!\left[ j\frac{4\pi R_0 f_\tau}{c}\left(\frac{1}{D(f_\eta)}-1\right) \right]

5. RCMC 是怎么从二维谱里“长出来”的

由上式中关于 f_\tau 的线性相位

\exp\!\left[-j\frac{4\pi R_0}{cD(f_\eta)}f_\tau\right]

可知，在不同 Doppler 频率 f_\eta 上，目标的等效距离时延为

\tau_d(f_\eta)=\frac{2R_0}{cD(f_\eta)}

而若没有距离徙动，理想参考时延应是

\tau_0=\frac{2R_0}{c}

因此需要补偿的时延差为

\Delta\tau(f_\eta;R_0) = \frac{2R_0}{c}\left(\frac{1}{D(f_\eta)}-1\right)

对应的斜距迁移量

\Delta R(f_\eta;R_0) = R_0\left(\frac{1}{D(f_\eta)}-1\right)

这就是 RD 算法里 RCMC 的核心公式。Moreira 的教程指出：RCM 来自目标斜距 R(\eta) 随合成孔径时间变化；若不校正，目标能量会跨多个距离单元，导致方位失焦。Chen 与 Kiang 则给出了与上式等价的 RCMC 相位滤波器形式。

工程实现时，RCMC 常在 (\tau,f_\eta) 域通过 沿距离向的插值重采样 完成：

把每个 Doppler 频率通道上的距离剖面按 \Delta\tau(f_\eta;R_0) 平移，使所有 f_\eta 上的目标都对齐到同一个参考距离单元。对小斜视 broadside 场景，这一步就是 RD 的关键近似与关键代价所在。

6. RCMC 后为什么可以做一维方位匹配滤波

做完 RCMC 后，二维谱近似变为

S_{\text{RCMC}}(f_\tau,f_\eta;R_0) \propto W_r(f_\tau)\,W_a(f_\eta)\, \exp\!\left(-j\frac{4\pi R_0}{c}f_\tau\right)\, \exp\!\left[-j\frac{4\pi R_0}{\lambda}D(f_\eta)\right]

现在距离向和方位向基本分离了：

第一项只表示目标在正确距离位置 R_0 上；
第二项是纯方位相位项。

于是方位匹配滤波器就是它的共轭：

H_{\text{az}}(f_\eta;R_0) = \exp\!\left[j\frac{4\pi R_0}{\lambda}D(f_\eta)\right]

乘上该滤波器后，对 f_\eta 做逆 FFT，目标在方位向被压缩成主瓣，得到聚焦图像。Chen 与 Kiang 给出的传统 RDA azimuth-compression filter 正是这个形式；Moreira 的教程则说明了方位压缩与方位参考函数的关系。

若再对 D(f_\eta) 做二阶展开，

D(f_\eta)\approx 1-\frac{\lambda^2 f_\eta^2}{8v^2}

则

-\frac{4\pi R_0}{\lambda}D(f_\eta) \approx -\frac{4\pi R_0}{\lambda} +\frac{\pi \lambda R_0}{2v^2}f_\eta^2

这说明方位参考函数在频域近似也是一个二次相位函数，与时域的 azimuth chirp 完全对应。也就是说，RD 的方位压缩本质就是“方位 LFM 的匹配滤波”。

7. 传统 RD 算法严格步骤

对原始数据 s_{\text{rb}}(\tau,\eta)，RDA 的标准流程就是：

s_{\text{rb}}(\tau,\eta) \;\xrightarrow{\text{range MF}}\; s_{\text{rc}}(\tau,\eta) \;\xrightarrow{\mathcal{F}_\eta}\; S_{\text{RD}}(\tau,f_\eta) \;\xrightarrow{\text{RCMC}}\; S_{\text{RCMC}}(\tau,f_\eta) \;\xrightarrow{\text{azimuth MF}}\; \widetilde S(\tau,f_\eta) \;\xrightarrow{\mathcal{F}_\eta^{-1}}\; I(\tau,\eta)

如果把距离压缩和 RCMC/方位压缩都写到频域里，也可表示为：

距离向 FFT：\mathcal{F}_\tau
乘距离匹配滤波器 H_r(f_\tau)
距离向 IFFT，得到 s_{\text{rc}}(\tau,\eta)
方位向 FFT，进入 RD 域
进行 RCMC（通常是距离插值；也可用相位法写出等效形式）
乘方位匹配滤波器 H_{\text{az}}(f_\eta;R_0)
方位向 IFFT，得到聚焦图像。

8. 推导使用近似

这部分很关键。RD 不是“完全精确算法”；它是从精确点目标模型出发，在以下条件下得到的高效近似算法：

8.1 几何近似

平台直线、匀速、单基地、stop-and-go。斜距模型先是精确的

R(\eta)=\sqrt{R_0^2+(v\eta)^2}

再常在方位参考函数构造时用到二阶展开。

8.2 小斜视/零斜视或弱耦合近似

传统 RDA 最适合小斜视。高斜视时，距离–方位耦合会更强，RCM 更严重，需额外补偿项或改进算法。Chen 与 Kiang专门讨论了高 squint 情况下传统 RDA 的不足并引入额外补偿滤波器。

8.3 窄带近似

在构造 D(f_\eta) 时，通常用载频 f_c 而不是 f_c+f_\tau，这等于忽略了部分距离频率与方位频率的高阶耦合；这也是为什么 RD 在宽带、高分辨率、宽孔径条件下不如 \omega-k 更稳健。综述文献也把 RD 归入“频域高效但带近似”的方法。

9. “核心数学骨架”

把上面的推导压缩成最本质的四条式子：

原始点目标回波

s_{\text{rb}}(\tau,\eta) = A_0\,w_a(\eta)\, w_r\!\left(\tau-\frac{2R(\eta)}{c}\right) e^{-j\frac{4\pi f_c}{c}R(\eta)} e^{j\pi K_r(\tau-\frac{2R(\eta)}{c})^2}

斜距历程

R(\eta)=\sqrt{R_0^2+(v\eta)^2} \approx R_0+\frac{v^2\eta^2}{2R_0}

方位调频率

K_a=-\frac{2v^2}{\lambda R_0}

这来自相位 -4\pi R(\eta)/\lambda 的二阶展开。

RCMC 量

\Delta R(f_\eta;R_0)=R_0\left(\frac{1}{D(f_\eta)}-1\right), \qquad D(f_\eta)=\sqrt{1-\left(\frac{\lambda f_\eta}{2v}\right)^2}

RCMC 做完后，再用

H_{\text{az}}(f_\eta;R_0)=\exp\!\left[j\frac{4\pi R_0}{\lambda}D(f_\eta)\right]

做方位压缩。

MATLAB代码

clc; clear; close all;
rng(0);

%% =========================================================
%  1) 基本参数（Broadside Stripmap SAR）
% ==========================================================
c = 3e8;                      % 光速
fc = 5.3e9;                   % 载频 (Hz)
lambda = c / fc;

v = 150;                      % 平台速度 (m/s)
PRF = 500;                    % 脉冲重复频率 (Hz)
Naz = 1024;                   % 方位向脉冲数（增大以覆盖离轴目标）

fs = 20e6;                    % 距离向采样率 (Hz)
Nr = 2048;                    % 每脉冲距离采样点数
Tp = 10e-6;                   % 发射 LFM 脉宽
B = 15e6;                     % LFM 带宽
Kr = B / Tp;                  % 距离向调频率

R0 = 10e3;                    % 场景中心斜距 (m)

% 慢时间轴（方位时间）
eta = ((0:Naz-1) - Naz/2) / PRF;     % s

% 快时间轴（距离时间）
tau0 = 2 * R0 / c;                    % 场景中心双程延迟
tr = tau0 + ((0:Nr-1) - Nr/2) / fs;   % s

%% =========================================================
%  2) 构造点目标场景
%     target = [x_azimuth(m), y_range_offset(m), amplitude]
%     x: 方位向位置（沿航迹）
%     y: 相对 R0 的距离偏移
% ==========================================================
targets = [
      0,     0,    1.0;
     80,    60,    0.8;
   -120,   -80,    0.7
];

%% =========================================================
%  3) 生成原始 SAR 回波 raw(Nr, Naz)
%     简化模型：Broadside + Baseband LFM
% ==========================================================
raw = complex(zeros(Nr, Naz));

for ia = 1:Naz
    eta_i = eta(ia);
    x_platform = v * eta_i;   % 平台当前位置（沿航迹）

    sig = complex(zeros(Nr, 1));

    for k = 1:size(targets, 1)
        xk = targets(k, 1);
        yk = targets(k, 2);
        ak = targets(k, 3);

        Rk0 = R0 + yk;  % 目标参考斜距

        % 当前脉冲时刻平台到点目标瞬时斜距
        Rk = sqrt(Rk0^2 + (x_platform - xk)^2);

        % 回波双程延迟
        tau_k = 2 * Rk / c;

        % 距离向 LFM 相位 + 传播相位
        t_shift = tr - tau_k;
        gate = abs(t_shift) <= Tp/2;

        echo = ak .* gate .* ...
               exp(1j*pi*Kr.*t_shift.^2) .* ...
               exp(-1j*4*pi*fc*Rk/c);

        sig = sig + echo(:);
    end

    raw(:, ia) = sig;
end

% 加一点复高斯噪声
raw = raw + 0.01 * (randn(size(raw)) + 1j*randn(size(raw)));

%% =========================================================
%  4) 距离压缩（Range Compression）
%     用居中 FFT 写法，避免零频点在数组中间带来的相位斜坡
% ==========================================================
t_ref = ((0:Nr-1) - Nr/2) / fs;
ref_chirp = ((abs(t_ref) <= Tp/2) .* exp(1j*pi*Kr*t_ref.^2)).';   % 列向量

% 匹配滤波器：h(t) = conj(s(-t))
H_range = fft(ifftshift(conj(ref_chirp)), Nr);   % 列向量频响

Raw_f = fft(ifftshift(raw, 1), Nr, 1);
rc = fftshift(ifft(Raw_f .* repmat(H_range, 1, Naz), Nr, 1), 1);

%% =========================================================
%  5) 方位 FFT，进入 Range-Doppler 域
% ==========================================================
RD = fftshift(fft(ifftshift(rc, 2), Naz, 2), 2);

fa = ((-Naz/2):(Naz/2-1)) * (PRF / Naz);  % 方位频率轴
fa = fa(:).';                             % 行向量

%% =========================================================
%  6) RCMC（Range Cell Migration Correction）
%     使用按距离单元变化的二次近似：
%     DeltaR(fa, R) ≈ lambda^2 * R * fa^2 / (8*v^2)
% ==========================================================
range_axis = tr * c / 2;      % 斜距轴
r_vec = range_axis(:);        % 列向量

RD_rcmc = complex(zeros(size(RD)));

for ia = 1:Naz
    fa_i = fa(ia);

    deltaR = (lambda^2 * r_vec * fa_i^2) / (8 * v^2);
    src_r = r_vec + deltaR;   % 从更大的斜距位置反取样

    RD_rcmc(:, ia) = interp1(r_vec, RD(:, ia), src_r, 'linear', 0);
end

%% =========================================================
%  7) 方位压缩（Azimuth Compression）
%     修正后的匹配滤波器符号：
%     Ka_mag(R) = 2*v^2 / (lambda*R)
%     H_az(fa, R) = exp(-j*pi*fa^2 / Ka_mag(R))
% ==========================================================
Ka_mag = 2 * v^2 ./ (lambda * r_vec);    % 列向量，取正值幅度

H_az = exp(-1j*pi * bsxfun(@rdivide, fa.^2, Ka_mag));
RD_focused = RD_rcmc .* H_az;

%% =========================================================
%  8) 方位 IFFT，得到聚焦图像
% ==========================================================
img = fftshift(ifft(ifftshift(RD_focused, 2), Naz, 2), 2);

img_abs = abs(img);
img_dB = 20 * log10(img_abs / (max(img_abs(:)) + eps) + 1e-12);

%% =========================================================
%  9) 坐标轴
% ==========================================================
az_axis = eta * v;            % 方位位置轴（平台投影近似）

%% =========================================================
%  10) 显示结果
% ==========================================================
figure;
imagesc(az_axis, range_axis, img_dB);
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('SAR Range-Doppler Focused Image');
colorbar;
caxis([-40 0]);
hold on;
plot(targets(:,1), R0 + targets(:,2), 'wx', 'LineWidth', 1.5, 'MarkerSize', 8);
hold off;

%% =========================================================
%  11) 显示原始数据与距离压缩结果
% ==========================================================
figure;
imagesc(az_axis, range_axis, 20*log10(abs(raw)/(max(abs(raw(:))) + eps) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('Raw SAR Data (dB)');
colorbar;
caxis([-40 0]);

figure;
imagesc(az_axis, range_axis, 20*log10(abs(rc)/(max(abs(rc(:))) + eps) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('After Range Compression (dB)');
colorbar;
caxis([-40 0]);

Back Projection算法

思路：对图像中的每个像素，计算该像素到平台在每个脉冲时刻的距离；然后从对应回波里取出那个距离上的复数样本，做相位校正并累加。所有脉冲都加完后，这个像素就聚焦了。相关资料直接指出，BP 需要 显式知道每个脉冲对应的目标距离；其基本处理就是对每个脉冲、每个像素做这种投影。

1. 建模目标

BP 的本质不是先把数据搬到某个统一频域里处理，而是：

对每个候选像素位置 \mathbf{x}，针对每个脉冲 n，计算该像素的理论双程距离 d_n(\mathbf{x})，然后从该脉冲的距离压缩数据中取出对应延迟处的复值样本，再补偿传播相位并做相干累加。这个“每个像素单独定制核函数”的观点，正是 BP 与 RD / CS / \omega-k 的根本区别。

2. 发射信号模型

先取最常见的单基地、LFM 脉冲 SAR。发射信号可写成

s_{\text{tx}}(t)=w(t)\exp\Big(j2\pi f_0 t+j\pi K t^2+j\phi_0\Big)

其中：

t：快时间
f_0：载频
K：调频率
w(t)：发射脉冲包络
\phi_0：初相

这正是 Duersch & Long 在其推导开头使用的 LFM 发射模型。

3. 单点散射体回波模型

设场景中有一个点散射体，位置为 \mathbf{x}。对第 n 个脉冲，平台位置为 \mathbf{r}_n。在单基地、stop-and-go 近似下，双程传播延迟为

\tau_n(\mathbf{x})=\frac{2R_n(\mathbf{x})}{c}, \qquad R_n(\mathbf{x})=\|\mathbf{r}_n-\mathbf{x}\|

对应接收信号是发射信号的延时缩放版本。Duersch & Long 写成“幅度调整、时间延迟的发射副本”，即

s^{\text{rx}}_n(t) = A_n(t-\tau_n)\,G_n(t-\tau_n)\,w(t-\tau_n) \exp\!\Big(j2\pi f_0(t-\tau_n)+j\pi K(t-\tau_n)^2+j\phi_0\Big)+\eta

其中 A_n,G_n 吸收了散射系数、天线增益、传播损耗、系统响应等项。

把信号下变频到基带后，载频的 t 项被去掉，保留传播相位：

s_n(t) = a_n(\mathbf{x})\,w(t-\tau_n) \exp\!\Big(-j2\pi f_0\tau_n+j\pi K(t-\tau_n)^2\Big)

这里我把所有慢变幅度合并为 a_n(\mathbf{x})。

4. 场景扩展到连续反射率

设地面复反射率为 \sigma(\mathbf{x})。在线性叠加、单次散射近似下，第 n 个脉冲的总回波为

s_n(t) = \int \sigma(\mathbf{x})\,a_n(\mathbf{x})\,w\!\big(t-\tau_n(\mathbf{x})\big) \exp\!\Big(-j2\pi f_0\tau_n(\mathbf{x})+j\pi K(t-\tau_n(\mathbf{x}))^2\Big)\,d\mathbf{x}

这一步是把单点目标模型推广成积分模型。BP 的任务就是由 \{s_n(t)\} 反演出 \sigma(\mathbf{x})。

5. 从“最优检测”角度得到匹配滤波形式

对于一个给定候选像素 \mathbf{x}_0，如果它真是一个点散射体，那么第 n 个脉冲下的期望模板应为

h_n(t;\mathbf{x}_0) = w\!\big(t-\tilde\tau_n(\mathbf{x}_0)\big) \exp\!\Big(-j2\pi f_0\tilde\tau_n(\mathbf{x}_0)+j\pi K(t-\tilde\tau_n(\mathbf{x}_0))^2\Big)

其中 \tilde\tau_n(\mathbf{x}_0) 是像素 \mathbf{x}_0 对应的预测延迟。

在加性白噪声下，对该像素做匹配滤波的输出就是

I(\mathbf{x}_0) = \sum_{n\in\mathcal N} \int s_n(\xi)\,h_n^*(\xi;\mathbf{x}_0)\,d\xi

Duersch & Long 直接给出了这个“单像素匹配滤波”的表达，并强调：因为核函数随像素变化，所以 BP 不是普通卷积，而是 空间变核匹配滤波。

把上式展开：

I(\mathbf{x}_0) = \sum_n \int s_n(\xi)\, w(\xi-\tilde\tau_n) \exp\!\Big(+j2\pi f_0\tilde\tau_n-j\pi K(\xi-\tilde\tau_n)^2\Big)\,d\xi

这就是 严格意义上的 BP 连续时间定义。到这里为止，还没有“先距离压缩、后方位叠加”的近似拆分；这是完整二维匹配滤波。Duersch & Long 也明确指出，这样做同时处理了传播相位和 LFM 啁啾，也就是同时做了方位与距离压缩。

6. 像素正确对焦时会发生什么

若候选像素正好等于真实散射体位置，即

\tilde\tau_n(\mathbf{x}_0)=\tau_n(\mathbf{x})

则匹配滤波器的相位与信号相位完全对消，输出在该像素处最大。Duersch & Long 在文中指出：如果每个脉冲的时延都已精确知道，那么滤波器相位和信号相位可以“完美匹配”；此时 BP 就是理想匹配滤波。

若存在延迟失配 \Delta\tau_n=\tilde\tau_n-\tau_n，就会产生残余传播相位

\exp\!\big(j2\pi f_0\Delta\tau_n\big) = \exp\!\big(jk\,\Delta d_n\big), \qquad k=\frac{2\pi}{\lambda}

其中 \Delta d_n=c\Delta\tau_n 为双程残余距离。文中把它称作 residual phase，并强调这是 BP 聚焦误差分析的关键量。

7. 推到工程常用形式：先距离压缩，再方位回投

上面是“严格二维匹配滤波”形式。工程实现通常进一步利用一个常见假设：

脉冲内平台运动可忽略，且距离压缩与方位相干积累可以分离。

在这个假设下，可先用只针对 LFM 啁啾的匹配滤波器做距离压缩。Duersch & Long 给出的距离压缩参考函数为

h_n^R(t)=w(t-\tau_n)\exp\!\big(j\pi K (t-\tau_n)^2\big)

对应距离压缩信号 g_n(t)。当包络是矩形窗时，文中给出近似结果：距离压缩后的回波等于一个 sinc 形主瓣乘以传播相位项。

把时间变量映射到双程距离，定义

d_n(\mathbf{x}) = 2R_n(\mathbf{x}), \qquad \tilde d_n(\mathbf{x}_0)=2R_n(\mathbf{x}_0), \qquad \Delta d_n=\tilde d_n-d_n.

文中随后给出方位压缩核

h_n^A = w_n \exp(-jk\tilde d_n)

并得到回投影结果

I(\mathbf{x}_0) = \sum_{n\in\mathcal N} g_n(\tilde d_n)\,w_n\,\exp(jk\tilde d_n)

这就是 BP 最常见的实用表达式。Duersch & Long 的式(22)正是这个形式。

对单基地 SAR，d_n=2R_n，因此

I(\mathbf{x}_0) = \sum_n w_n\, g_n\!\left(\frac{2R_n(\mathbf{x}_0)}{c}\right) \exp\!\left(j\frac{4\pi}{\lambda}R_n(\mathbf{x}_0)\right)

这条式子就是你在绝大多数代码里看到的 BP 核心。

8. 为什么叫“Back Projection”

现在可以严格解释“后向投影”：

第 n 个脉冲的距离压缩数据 g_n(t) 本质上只说明“能量落在哪个双曲等距面上”；
对某个像素 \mathbf{x}_0，我们计算它相对于该脉冲的理论双程延迟 \tau_n(\mathbf{x}_0)；
从 g_n(t) 中取出 t=\tau_n(\mathbf{x}_0) 处的值，相当于把这一脉冲在该等距面上的能量“投回”到像素 \mathbf{x}_0；
再补偿该像素的传播相位并对所有脉冲相干求和。

所以 BP 不是在二维频域里做统一重采样，而是 把每个脉冲对每个像素的贡献沿真实传播几何直接回投到图像域。

9. 离散实现公式

设距离压缩后的离散数据为 g[n,m]，其中

n：脉冲索引
m：距离采样索引
采样间隔为 \Delta t_r

对像素 \mathbf{x}_p，先算

\tau_n(\mathbf{x}_p)=\frac{2R_n(\mathbf{x}_p)}{c}, \qquad m_n(\mathbf{x}_p)=\frac{\tau_n(\mathbf{x}_p)-t_0}{\Delta t_r}

因为 m_n 往往不是整数，需要插值：

\hat g_n(\mathbf{x}_p)=\operatorname{Interp}\big(g[n,:],\,m_n(\mathbf{x}_p)\big)

于是离散 BP 写成

I[\mathbf{x}_p] = \sum_{n\in\mathcal N_p} w_n(\mathbf{x}_p)\, \hat g_n(\mathbf{x}_p)\, \exp\!\left(j\frac{4\pi}{\lambda}R_n(\mathbf{x}_p)\right)

其中 \mathcal N_p 可加上天线波束支持约束，只累加真正照射到该像素的脉冲。Duersch & Long 还指出，实际中像天线增益这类幅度项也可并入 w_n。

MATLAB代码

clc; clear; close all;

%% =========================================================
%  1) 系统参数
% ==========================================================
c  = 3e8;                 % 光速
fc = 5.3e9;               % 载频 (Hz)
lambda = c / fc; %#ok<NASGU>

v   = 120;                % 平台速度 (m/s)
PRF = 400;                % 脉冲重复频率 (Hz)
Na  = 256;                % 方位脉冲数

fs = 20e6;                % 距离向采样率 (Hz)
Nr = 2048;                % 每脉冲距离采样点数
Tp = 8e-6;                % LFM 脉宽 (s)
B  = 10e6;                % LFM 带宽 (Hz)
Kr = B / Tp;              % 距离调频率

R0 = 8000;                % 场景中心斜距 (m)

%% =========================================================
%  2) 平台轨迹（直线 broadside 示例）
%     这里采用 2D 平面模型：平台沿 x 方向运动，目标位于 z=0 平面
% ==========================================================
eta = ((0:Na-1) - (Na-1)/2) / PRF;  % 用 (N-1)/2 避免半个脉冲偏置
xPlat = v * eta;
yPlat = zeros(size(xPlat));
zPlat = zeros(size(xPlat));

posPlat = [xPlat(:), yPlat(:), zPlat(:)];

%% =========================================================
%  3) 距离快时间轴
%     tFast: 相对场景中心双程时延 tau0 的快时间
%     tr   : 绝对快时间
% ==========================================================
tau0  = 2 * R0 / c;
tFast = ((0:Nr-1) - (Nr-1)/2) / fs;
tr    = tau0 + tFast;

%% =========================================================
%  4) 发射 LFM 参考信号（基带）
% ==========================================================
tref = (-Tp/2 : 1/fs : Tp/2 - 1/fs);
tx = exp(1j * pi * Kr * tref.^2); %#ok<NASGU>

%% =========================================================
%  5) 构造点目标场景
%     target = [x(m), y(m), amplitude]
%     本例中 y 是场景横向距离坐标；由于 z=0，x=0 处其数值等于斜距
% ==========================================================
targets = [
     0,   R0,      1.0;
    50,   R0+40,   0.8;
   -80,   R0-60,   0.7
];

%% =========================================================
%  6) 生成原始 SAR 回波 raw(Nr, Na)
%     对每个目标，先算真实斜距 R，再转成相对场景中心的时延 dt
% ==========================================================
raw = complex(zeros(Nr, Na));

for ia = 1:Na
    xa = posPlat(ia,1);
    ya = posPlat(ia,2);
    za = posPlat(ia,3);

    sig = complex(zeros(Nr, 1));

    for k = 1:size(targets, 1)
        xt = targets(k,1);
        yt = targets(k,2);
        zt = 0;
        at = targets(k,3);

        R  = sqrt((xa - xt)^2 + (ya - yt)^2 + (za - zt)^2);
        dt = 2 * (R - R0) / c;

        tshift = tFast - dt;
        gate   = abs(tshift) <= Tp/2;

        echo = at .* gate .* ...
               exp(1j * pi * Kr * tshift.^2) .* ...
               exp(-1j * 4 * pi * fc * R / c);

        sig = sig + echo(:);
    end

    raw(:, ia) = sig;
end

% 加噪声
raw = raw + 0.01 * (randn(size(raw)) + 1j * randn(size(raw)));

%% =========================================================
%  7) 距离压缩（matched filter）
%     注意：raw/ref 都是"零时刻在中间"的排列，FFT 前要 ifftshift，
%     IFFT 后再 fftshift 回来，这样压缩峰值仍落在正确的快时间位置
% ==========================================================
ref = (abs(tFast) <= Tp/2) .* exp(1j * pi * Kr * tFast.^2);
Hmf = conj(fft(ifftshift(ref(:)), Nr));

rawFFT = fft(ifftshift(raw, 1), Nr, 1);
rc = fftshift(ifft(rawFFT .* repmat(Hmf, 1, Na), Nr, 1), 1);

%% =========================================================
%  8) 成像网格
% ==========================================================
xImg = -150 : 2 : 150;
yImg = (R0 - 120) : 2 : (R0 + 120);

Nx = length(xImg);
Ny = length(yImg);

img = complex(zeros(Ny, Nx));

%% =========================================================
%  9) Back Projection
%     对每个像素：
%       - 计算每个脉冲对应的斜距 R
%       - 换算成相对 tau0 的时延 dt
%       - 在 rc(:, ia) 上按 tFast 插值
%       - 乘 exp(+j4*pi*fc*R/c) 去传播相位
%       - 相干累加
% ==========================================================
for ix = 1:Nx
    for iy = 1:Ny
        xp = xImg(ix);
        yp = yImg(iy);
        zp = 0;

        s = 0;

        for ia = 1:Na
            xa = posPlat(ia,1);
            ya = posPlat(ia,2);
            za = posPlat(ia,3);

            R  = sqrt((xa - xp)^2 + (ya - yp)^2 + (za - zp)^2);
            dt = 2 * (R - R0) / c;

            val = interp1(tFast, rc(:, ia), dt, 'linear', 0);
            s = s + val * exp(1j * 4 * pi * fc * R / c);
        end

        img(iy, ix) = s;
    end
end

%% =========================================================
%  10) 显示结果
% ==========================================================
imgAbs = abs(img);
imgdB = 20 * log10(imgAbs / max(imgAbs(:)) + 1e-12);

figure;
imagesc(xImg, yImg, imgdB);
axis xy;
xlabel('Azimuth x (m)');
ylabel('Scene y (m)');
title('SAR Back Projection Image');
colorbar;
caxis([-40 0]);

%% =========================================================
%  11) 中间结果查看
% ==========================================================
rangeAxis = tr * c / 2;

figure;
imagesc(xPlat, rangeAxis, 20 * log10(abs(raw) / max(abs(raw(:))) + 1e-12));
axis xy;
xlabel('Azimuth position (m)');
ylabel('Slant range (m)');
title('Raw SAR Data (dB)');
colorbar;
caxis([-40 0]);

figure;
imagesc(xPlat, rangeAxis, 20 * log10(abs(rc) / max(abs(rc(:))) + 1e-12));
axis xy;
xlabel('Azimuth position (m)');
ylabel('Slant range (m)');
title('Range-Compressed Data (dB)');
colorbar;
caxis([-40 0]);

Chirp Scaling算法

思路：RD 的痛点在 RCMC 插值。于是 CSA 通过一个“啁啾尺度变换”的相位操作，把原本要靠插值完成的距离徙动校正，改成了 频域相位补偿。多篇资料都明确说，CSA 是通过 以相位补偿替代 RDA 的插值 来实现的，算法结构主要由 FFT + 相位补偿 组成。

1. 几何、记号与起点

设平台沿方位向 x 以恒速 v 飞行，点目标位于最近斜距 R_0 处。快时间记为 \tau，慢时间记为 \eta。载频 f_0，波长 \lambda=c/f_0，发射 LFM 脉冲调频率 K_r。

点目标瞬时斜距

R(\eta)=\sqrt{R_0^2+v^2\eta^2}

基带回波可写成

s(\tau,\eta) = A\,w_r\!\left(\tau-\frac{2R(\eta)}{c}\right) w_a(\eta) \exp\!\left[-j\frac{4\pi f_0}{c}R(\eta)\right] \exp\!\left[j\pi K_r\left(\tau-\frac{2R(\eta)}{c}\right)^2\right]

这就是 CSA、RDA、\omega-k 等算法共同的点目标模型起点；CSA 适用于二维频域处理，并用相位补偿完成 differential RCMC 与 bulk RCMC。

2. 先做方位 FFT：进入距离-多普勒域

对 \eta 做 FFT，令方位频率为 f_a。在驻相近似下，点目标信号可写成

S_1(\tau,f_a;R_0) = A_1\, W_a(f_a)\, w_r\!\left(\tau-\frac{2R_0}{cD(f_a)}\right) \exp\!\left[-j\frac{4\pi f_0R_0}{c}D(f_a)\right] \exp\!\left[ j\pi K_m(f_a) \left(\tau-\frac{2R_0}{cD(f_a)}\right)^2 \right]

（1）

其中

D(f_a)=\sqrt{1-\left(\frac{\lambda f_a}{2v}\right)^2}

称为 range migration factor；而 K_m(f_a) 是进入距离-多普勒域后的等效距离调频率。经典 CSA 文献与后续实现文献都采用这种表述：在方位 FFT 后，同一最近斜距上的目标在 (\tau,f_a) 域里“塌缩”为一条轨迹，RCM 体现在轨迹中心

\tau_c(f_a;R_0)=\frac{2R_0}{cD(f_a)}

也就是说，目标的距离徙动不是常数平移，而是随 f_a 变化。

3. 为什么需要 chirp scaling

RDA 的做法是：对每个 f_a 切片，把轨迹中心 \tau_c(f_a;R_0) 插值搬回参考位置。

CSA 的做法是：不用插值，而是先对快时间作一个依赖 f_a 的“尺度变换”，把所有目标的 RCM 变成“和距离无关的公共形状 + 一个常量平移”，这样后面就能用统一的相位函数补掉。CSA 的提出本意正是避免 RCMC 插值，同时保持相位保真。

4. 参考距离与尺度因子

取参考斜距 R_{\rm ref}。在这个参考斜距上定义参考 migration factor

D_{\rm ref}(f_a)=D(f_a)

这里对零斜视直线匀速模型，D 本身不显含 R；真正起作用的是参考轨迹中心

\tau_{\rm ref}(f_a)=\frac{2R_{\rm ref}}{cD(f_a)}

CSA 构造的目标是：

把任意 R_0 的轨迹中心

\frac{2R_0}{cD(f_a)}

经过一个依赖 f_a 的 chirp scaling 变换后，改写成

\frac{2R_{\rm ref}}{cD(f_a)}+\frac{2(R_0-R_{\rm ref})}{c}

即把原本随 f_a 变化的距离徙动，变成“参考轨迹 + 与 f_a 无关的常量距离偏移”。这一步就是 differential RCMC。

为此定义 chirp scaling 系数

C_s(f_a)=\frac{D(f_{a,\rm ref})}{D(f_a)}-1

在零多普勒参考下常取 f_{a,\rm ref}=0，于是 D(0)=1，得到

C_s(f_a)=\frac{1}{D(f_a)}-1 \tag{2}

这正是很多工程实现里出现的 \beta(f_a)=1/D(f_a)-1 形式。

5. 第一相位函数：实现 differential RCMC

CSA 的第一步相位函数写成

H_1(\tau,f_a) = \exp\!\left\{ j\pi K_m(f_a)\,C_s(f_a)\, \left(\tau-\tau_{\rm ref}(f_a)\right)^2 \right\}

把它乘到式 (1) 上：

S_2(\tau,f_a;R_0)=S_1(\tau,f_a;R_0)\,H_1(\tau,f_a)

把二次项展开。设

u=\tau-\frac{2R_0}{cD(f_a)},\qquad u_{\rm ref}=\tau-\frac{2R_{\rm ref}}{cD(f_a)}

则

u^2+C_s u_{\rm ref}^2 = (1+C_s)\tau^2 -\frac{4}{cD}(R_0+C_sR_{\rm ref})\tau +\frac{4}{c^2D^2}(R_0^2+C_sR_{\rm ref}^2)

若取 1+C_s=1/D，即式 (2)，则线性项对应的“新中心”满足

\tau_c'(f_a;R_0) = \frac{2}{c} \frac{R_0+C_sR_{\rm ref}}{1+C_s}\frac{1}{D} = \frac{2}{c}\left(R_{\rm ref}+\frac{R_0-R_{\rm ref}}{1}\,D(0)\right)

在 D(0)=1 的参考取法下化为

\tau_c'(f_a;R_0) = \frac{2R_{\rm ref}}{cD(f_a)}+\frac{2(R_0-R_{\rm ref})}{c} \tag{4}

这就是 CSA 的关键：

原先 RCM 是 \propto R_0/D(f_a)；

经过 chirp scaling 后，变成“参考轨迹 \propto R_{\rm ref}/D(f_a)”加上与 f_a 无关的平移 (R_0-R_{\rm ref})。

因此，所有点目标的徙动形状被“等化”了。这就是文献所说的用第一相位函数完成 differential RCMC / chirp scaling。

6. 再做距离 FFT：进入二维频域

对 \tau 做 FFT，记距离频率为 f_\tau。得到二维频域表达式

S_3(f_\tau,f_a;R_0) = A_3\,W_a(f_a)\,W_r\!\left(\frac{f_\tau}{K_m'(f_a)}\right) \exp\{j\Phi_3(f_\tau,f_a;R_0)\} \tag{5}

其中相位可以整理成五类项：

\Phi_3 = \Phi_{\rm azi}(f_a;R_0) +\Phi_{\rm rc/src}(f_\tau,f_a) +\Phi_{\rm lin}(f_\tau;R_0) +\Phi_{\rm bulk}(f_\tau,f_a) +\Phi_{\rm const}(f_a;R_0)

在经典 CSA 的解释里，这些项中最重要的是：

一个 新的距离调频项，代表缩放后的 range chirp，并含有 SRC 所需的距离-方位耦合补偿；
一个 bulk RCMC 线性项，对应参考距离 R_{\rm ref} 的公共徙动；
一个只与目标真实距离偏移 R_0-R_{\rm ref} 有关的线性项。

后续第二相位函数就是专门针对这几项而设计。实现文献对这一点说得很明确：二维频域里第二个相位函数同时完成 range compression、SRC、bulk RCMC。

7. 第二相位函数：同时做 RC、SRC、bulk RCMC

在二维频域中乘第二相位函数

H_2(f_\tau,f_a) = \exp\!\left[ -j\pi \frac{f_\tau^2}{K_{\rm eq}(f_a)} \right] \, \exp\!\left[ -j\frac{4\pi}{c}R_{\rm ref}\Big(\frac{1}{D(f_a)}-1\Big)f_\tau \right] \tag{6}

这里：

第一因子对消缩放后的距离 chirp，完成 range compression；
同时因为 K_{\rm eq}(f_a) 依赖 f_a，它也完成 secondary range compression (SRC)；
第二因子对消参考轨迹的线性频率项，完成 bulk RCMC。

这样得到

S_4(f_\tau,f_a;R_0)=S_3(f_\tau,f_a;R_0)\,H_2(f_\tau,f_a) \tag{7}

这一点与后续工程论文的概括完全一致：第一相位函数做 differential RCMC，第二相位函数做 RC + SRC + bulk RCMC，第三相位函数做 azimuth compression 和 residual phase compensation。

8. 距离 IFFT 之后会变成什么

对 f_\tau 做 IFFT 回到距离-多普勒域，得到

S_5(\tau,f_a;R_0) = A_5\, W_a(f_a)\, p_r\!\left(\tau-\frac{2R_0}{c}\right) \exp\!\left[j\Phi_5(f_a;R_0)\right] \tag{8}

这里 p_r(\cdot) 是距离压缩后的主瓣。最重要的是：

\tau=\frac{2R_0}{c}

已经 不再依赖 f_a。

也就是说，RCMC 已经全部完成了，而且是通过相位补偿完成的，不需要 RDA 那种按 f_a 插值搬移。这个性质正是 CSA 的根本优势。

9. 第三相位函数：方位压缩与残余相位补偿

此时剩余相位 \Phi_5(f_a;R_0) 中仍有两部分：

目标的方位调频项（azimuth chirp）；
经过前两步变换留下的 residual phase。

因此设计第三相位函数

H_3(\tau,f_a) = \exp\!\left[-j\Phi_{\rm az\_mf}(f_a;\tau)\right] \exp\!\left[-j\Phi_{\rm res}(f_a;\tau)\right] \tag{9}

其中 \Phi_{\rm az\_mf} 是方位匹配滤波器，相当于 RDA 中的 azimuth compression；\Phi_{\rm res} 是 residual phase correction。乘上它以后得到

S_6(\tau,f_a;R_0)=S_5(\tau,f_a;R_0)\,H_3(\tau,f_a)

最后做方位 IFFT，得到聚焦图像

s_{\rm img}(\tau,\eta) = A_{\rm img}\, p_r\!\left(\tau-\frac{2R_0}{c}\right) p_a(\eta-\eta_0) e^{j\phi_{\rm targ}} \tag{10}

这就是一个在距离和方位两个维度都被压缩好的点目标响应。后续实现文献也把第三相位函数明确描述为：补偿 azimuth modulation 与 residual phase，然后 azimuth IFFT 得到图像。

10. 整体“严格的代数逻辑”

CSA 的数学本质可以概括为下面四步：

\underbrace{\frac{R_0}{D(f_a)}}_{\text{原始 RCM}} \;\xrightarrow{\;H_1\;} \; \underbrace{\frac{R_{\rm ref}}{D(f_a)} + (R_0-R_{\rm ref})}_{\text{等化后的 RCM}} \;\xrightarrow{\;H_2\;} \; \underbrace{R_0}_{\text{RC + SRC + bulk RCMC 后}} \;\xrightarrow{\;H_3\;} \; \underbrace{\delta_{\rm range}\delta_{\rm azimuth}}_{\text{聚焦点目标}}

所以它不是“神奇地不做 RCMC”，而是：

先用 chirp scaling 把不同距离处目标的徙动轨迹变成同一个参考轨迹族；
再在二维频域用统一的相位函数消掉这个公共轨迹和缩放后 chirp；
最后做方位匹配滤波。

这也正是原始论文和后续综述/实现文献对 CSA 的共同描述。

11. 这套推导的适用前提

上面的“严格推导”实际上依赖以下条件：

直线匀速平台；
stripmap 几何；
stop-go 近似；
点目标二维频谱可用驻相近似表达；
原始 CSA 主要针对零/小斜视与标准机/星载几何。

当斜视很大、轨迹误差显著、ScanSAR/Spotlight/TOPS、或高阶耦合很强时，原始 CSA 需要扩展为 ECS、NLCS 等变体。经典论文明确说明了原始 CSA 的优势是避免插值且适合 wide-swath / large-beamwidth / large-squint；而后续 ECS、squint CSA、NLCS 论文则是在更复杂几何下继续修正其相位模型。

12. 易混淆点：CSA 为什么“看起来像 RDA，但又不是 RDA”

它和 RDA 的共同点：

都先做 azimuth FFT；
都要做 RC、SRC、RCMC、azimuth compression；
都在距离-多普勒/二维频域里工作。

不同点在于：

RDA：RCMC 主要表现为插值重采样；
CSA：RCMC 被拆成
- differential RCMC（靠 H_1 的 chirp scaling），
- bulk RCMC（靠 H_2 的线性相位），
  于是整个流程变成 FFT/IFFT + phase multiply。

这正是 CSA 被认为特别适合硬件实现的原因。

MATLAB代码

clc; clear; close all;

%% =========================================================
%  1) Basic parameters
% ==========================================================
c  = 3e8;                  % Speed of light
fc = 5.3e9;                % Carrier frequency (Hz)
lambda = c / fc;

v   = 150;                 % Platform speed (m/s)
PRF = 600;                 % Pulse repetition frequency (Hz)
Naz = 512;                 % Number of azimuth pulses

fs = 20e6;                 % Range sampling rate (Hz)
Nr = 2048;                 % Number of fast-time samples per pulse
Tp = 10e-6;                % LFM pulse width (s)
B  = 15e6;                 % LFM bandwidth (Hz)
Kr = B / Tp;               % Range chirp rate

Rref = 10e3;               % Reference slant range (m)

%% =========================================================
%  2) Centered axes
%     Keep time/frequency axes consistent with fftshift/ifftshift usage.
% ==========================================================
eta   = ((0:Naz-1) - Naz/2) / PRF;
xPlat = v * eta;

tauRef = 2 * Rref / c;
tFast  = ((0:Nr-1) - Nr/2) / fs;
tr     = tauRef + tFast;

fa = ((0:Naz-1) - Naz/2) * (PRF / Naz);
fr = ((0:Nr-1) - Nr/2) * (fs  / Nr);

%% =========================================================
%  3) Point-target scene
%     [azimuth_x(m), range_offset_from_Rref(m), amplitude]
% ==========================================================
targets = [
   -20,   0,   1.0;
     0,  30,   0.8;
    20, -30,   0.7
];

%% =========================================================
%  4) Raw SAR echo simulation raw(Nr, Naz)
%     Fast time is referenced to tauRef, so delay is dt = 2*(R-Rref)/c.
% ==========================================================
raw = complex(zeros(Nr, Naz));

for ia = 1:Naz
    xp = xPlat(ia);
    sig = complex(zeros(Nr, 1));

    for k = 1:size(targets, 1)
        xt = targets(k, 1);
        dR = targets(k, 2);
        A  = targets(k, 3);

        R0 = Rref + dR;
        R  = sqrt(R0^2 + (xp - xt)^2);

        dt     = 2 * (R - Rref) / c;
        tShift = tFast - dt;
        gate   = abs(tShift) <= Tp / 2;

        echo = A .* gate .* ...
               exp(1j * pi * Kr * tShift.^2) .* ...
               exp(-1j * 4 * pi * fc * R / c);

        sig = sig + echo(:);
    end

    raw(:, ia) = sig;
end

raw = raw + 0.01 * (randn(size(raw)) + 1j * randn(size(raw)));

%% =========================================================
%  5) CSA Step-1: Azimuth FFT -> Range-Doppler domain
%     Because eta is centered, use ifftshift before FFT.
% ==========================================================
S_rd = fftshift(fft(ifftshift(raw, 2), Naz, 2), 2);

%% =========================================================
%  6) CSA Step-2: Chirp scaling phase
%     Teaching-version broadside approximation.
% ==========================================================
argD = 1 - (lambda * fa / (2 * v)).^2;
validFa = argD > 0;

Dfa = ones(1, Naz);
Dfa(validFa) = sqrt(argD(validFa));

S_rd(:, ~validFa) = 0;

[TFast, DFa] = ndgrid(tFast, Dfa);
dTauCS = tauRef * (1 ./ DFa - 1);

H_cs = exp(1j * pi * Kr .* (1 ./ DFa - 1) .* (TFast - dTauCS).^2);
S_cs = S_rd .* H_cs;

%% =========================================================
%  7) CSA Step-3: Range FFT -> 2D frequency domain
% ==========================================================
S_2df = fftshift(fft(ifftshift(S_cs, 1), Nr, 1), 1);

[FR, DFa2] = ndgrid(fr, Dfa);

%% =========================================================
%  8) CSA Step-4: Range compression + bulk RCMC
%
%  Fixes:
%    1) Range matched filter sign must be positive.
%    2) Bulk RCMC delay uses (1/D(fa) - 1), not (D(fa) - 1).
% ==========================================================
Km = Kr ./ DFa2;   % Teaching approximation: chirp rate scales with D(fa)

H_r = exp(1j * pi * (FR.^2) ./ Km);

dTauBulk = tauRef * (1 ./ DFa2 - 1);
H_bulk = exp(1j * 2 * pi * FR .* dTauBulk);

S_filt = S_2df .* H_r .* H_bulk;

%% =========================================================
%  9) CSA Step-5: Range IFFT -> back to Range-Doppler domain
% ==========================================================
S_rd2 = fftshift(ifft(ifftshift(S_filt, 1), Nr, 1), 1);

%% =========================================================
% 10) CSA Step-6: Azimuth compression
%
%     The original script used exp(-j*pi*fa^2/Ka), but Ka < 0 here,
%     so the matched-filter sign must be flipped.
% ==========================================================
rangeAxis = Rref + tFast * c / 2;
[rangeMat, faMat] = ndgrid(rangeAxis, fa);

Ka = -2 * v^2 ./ (lambda * rangeMat);
H_az = exp(1j * pi * (faMat.^2) ./ Ka);

S_az = S_rd2 .* H_az;

%% =========================================================
% 11) Azimuth IFFT -> focused image
% ==========================================================
img = fftshift(ifft(ifftshift(S_az, 2), Naz, 2), 2);

imgAbs = abs(img);
imgdB  = 20 * log10(imgAbs / max(imgAbs(:)) + 1e-12);

azAxis = xPlat;

%% =========================================================
% 12) Display results
% ==========================================================
figure;
imagesc(azAxis, rangeAxis, imgdB);
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('SAR Chirp Scaling Image');
colorbar;
caxis([-40 0]);
hold on;
plot(targets(:,1), Rref + targets(:,2), 'rx', 'LineWidth', 1.5, 'MarkerSize', 10);
hold off;

figure;
imagesc(azAxis, rangeAxis, 20 * log10(abs(raw) / max(abs(raw(:))) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('Raw SAR Data (dB)');
colorbar;
caxis([-40 0]);

figure;
imagesc(azAxis, rangeAxis, 20 * log10(abs(S_rd2) / max(abs(S_rd2(:))) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('After CSA Range Processing (dB)');
colorbar;
caxis([-40 0]);

ω-k算法

思路：ω-k 是更“二维频域化”的方法。它的标志性步骤是 Stolt mapping（Stolt 映射）。经典比较论文指出，ω-k 的这个步骤可以 同时完成残余 RCMC、二次距离压缩（SRC）和方位压缩。

1. 几何与信号模型

设平台沿方位向以等效速度 V_r 匀速运动，慢时间为 \eta，快时间为 \tau。一个地面点目标在最近距时刻的斜距为 R_0。在标准 stripmap 近似下，它的斜距历程是双曲线

R(\eta)=\sqrt{R_0^2+V_r^2\eta^2}

这一“hyperbolic range equation”正是 ω-k 推导的基础；Cumming 文中后续二维频域相位就是建立在这个距离历程上的。

发射线性调频脉冲（中心频率 f_0，调频率 K_r）后，经解调得到点目标回波，可写成

s(\tau,\eta) = \sigma\, w_r\!\left(\tau-\frac{2R(\eta)}{c}\right) w_a(\eta)\, \exp\!\left[ -j\frac{4\pi f_0}{c}R(\eta) \right] \exp\!\left[ j\pi K_r\left(\tau-\frac{2R(\eta)}{c}\right)^2 \right]

这里 \sigma 是散射系数，w_r,w_a 分别是距离向和方位向加窗。这个式子本质上就是“延时的 LFM 回波 + 传播相位”。这是 SAR 标准点目标模型；后面 ω-k 的所有步骤都只是对这个模型做更高效的匹配滤波实现。

2. 距离向傅里叶变换：得到距离频域信号

对快时间 \tau 做傅里叶变换，记距离频率为 f_\tau。LFM 脉冲的频域形式会给出一个二次相位项，于是得到

S(f_\tau,\eta) = A_r(f_\tau)\,w_a(\eta)\, \exp\!\left[ -j\frac{4\pi}{c}(f_0+f_\tau)R(\eta) \right] \exp\!\left[ -j\pi \frac{f_\tau^2}{K_r} \right]

其中 A_r(f_\tau) 吸收了幅度和距离包络。这个形式和后面文献里的二维频域相位是一致的：最后那项 -\pi f_\tau^2/K_r 就是距离 LFM 的频域二次相位项。

3. 方位向傅里叶变换：得到二维频域相位

再对慢时间 \eta 做傅里叶变换，方位频率记为 f_\eta：

S_{2D}(f_\tau,f_\eta) = \int S(f_\tau,\eta)\, e^{-j2\pi f_\eta \eta}\,d\eta.

将上式的总相位记为

\Phi(\eta) = -\frac{4\pi}{c}(f_0+f_\tau)R(\eta)-2\pi f_\eta\eta

用驻相法（POSP）求主贡献点 \eta_s，令

\frac{d\Phi}{d\eta}=0

因为

\frac{dR}{d\eta}=\frac{V_r^2\eta}{R(\eta)}

所以有

-\frac{4\pi}{c}(f_0+f_\tau)\frac{V_r^2\eta_s}{R(\eta_s)} -2\pi f_\eta = 0

整理得

f_\eta = -\frac{2(f_0+f_\tau)}{c}\frac{V_r^2\eta_s}{R(\eta_s)}

把它和

R(\eta_s)^2=R_0^2+V_r^2\eta_s^2

联立消去 \eta_s，可得驻相点对应的斜距

R_s = R(\eta_s) = R_0\, \frac{f_0+f_\tau} {\sqrt{(f_0+f_\tau)^2-\dfrac{c^2f_\eta^2}{4V_r^2}}}

把 \eta_s 代回相位，得到二维频域相位主项

\theta_{2Df}(f_\tau,f_\eta) = -\frac{4\pi R_0}{c} \sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}} -\pi\frac{f_\tau^2}{K_r}

这正是 Cumming 和综述文献列出的 ω-k 关键二维频域相位表达式。

因此二维频域信号可写成

S_{2D}(f_\tau,f_\eta) = A(f_\tau,f_\eta)\, \exp\!\big(j\theta_{2Df}(f_\tau,f_\eta)\big).

4. 参考函数乘法：bulk compression

ω-k 的第一步聚焦不是直接“完全匹配”，而是先选一个参考斜距 R_{\text{ref}} 做参考函数乘法（RFM）：

H_{\text{ref}}(f_\tau,f_\eta) = \exp\!\left[ +j\frac{4\pi R_{\text{ref}}}{c} \sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}} +j\pi\frac{f_\tau^2}{K_r} \right]

把 S_{2D} 乘上它以后，距离 LFM 项被消掉，且参考距离 R_{\text{ref}} 处的目标被完全聚焦。剩余相位变成

\theta_{\text{RFM}}(f_\tau,f_\eta) \approx -\frac{4\pi(R_0-R_{\text{ref}})}{c} \sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}}

这一步就是所谓 bulk compression；文献中给出的形式与上式一致，并明确说“参考距离上的目标被正确聚焦，其他距离只部分聚焦”。

5. 为什么还没聚焦：耦合没有消掉

注意上式中还保留着

\sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}}

它把距离频率 f_\tau 和方位频率 f_\eta 耦合在一起。只要这个根号还在，二维逆 FFT 后的响应就不是一个在 (\tau,\eta) 上可分离的 sinc 型聚焦点。Cumming 明确指出，这个根号项里同时含有残余 RCMC、方位调制和距离-方位耦合，Stolt 映射正是用来一次性把它们校正掉的。

6. Stolt 映射：把耦合变量“拉直”

现在定义新的距离频率变量 f'_\tau，使

f_0+f'_\tau = \sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}}

这就是 Stolt 映射。文献里直接把这一步写成“沿距离频率轴的映射”，用它把上面的根号替换成新的线性变量。

于是剩余相位立刻变成

\theta_{\text{stolt}}(f'_\tau,f_\eta) = -\frac{4\pi(R_0-R_{\text{ref}})}{c}(f_0+f'_\tau)

这已经是对 f'_\tau 线性的相位。综述文献把这个结果直接写成 Stolt 后的相位公式。

更直观地，在波数域写会更漂亮。令

k_r=\frac{4\pi}{c}(f_0+f_\tau),\qquad k_x=\frac{2\pi f_\eta}{V_r}

则二维频域相位主项可写成

\theta_{\text{RFM}} = -(R_0-R_{\text{ref}}) \sqrt{k_r^2-k_x^2}

再定义

k_y=\sqrt{k_r^2-k_x^2}

则相位直接变为

\theta_{\text{stolt}} = -(R_0-R_{\text{ref}})k_y

这说明 ω-k 的本质就是：把原先采样在曲线 k_y=\sqrt{k_r^2-k_x^2} 上的数据，重采样到规则的 (k_x,k_y) 直角网格 上。完成这一步后，信号相位就与目标位置线性对应，可以直接用二维 IFFT 成像。Cumming 也用“scaling + shifting”的解释说明了 Stolt 映射如何同时完成残余 RCMC 和方位压缩。

7. 二维逆 FFT：得到聚焦图像

完成 Stolt 重采样后，数据已在规则的 (f'_\tau,f_\eta) 或 (k_y,k_x) 网格上。此时做二维逆 FFT：

I(r,x) = \mathcal{F}^{-1}_{2D} \left\{ \widetilde{S}(k_x,k_y) \right\}

即可得到图像。文献标准流程就是：

2D FFT
参考函数乘法（bulk compression）
Stolt 映射（differential compression）
2D IFFT。

8. 从“匹配滤波”角度看 ω-k 为什么对

理想点目标的二维匹配滤波核，本质上要共轭补偿传播相位

\exp\!\left[ -j\frac{4\pi}{c}(f_0+f_\tau)R(\eta) \right]

问题在于它在二维频域里不是一个对所有距离都统一可分离的简单乘法核。ω-k 的做法不是直接在原 (f_\tau,f_\eta) 网格上硬乘一个完全匹配滤波器，而是：

先用 R_{\text{ref}} 做一遍统一的“粗聚焦”；
再用 Stolt 映射把所有不同距离的点目标频谱都搬运到同一种线性相位结构；
最后用 IFFT 完成精聚焦。

因此，Stolt 映射不是附属小修正，而是 ω-k 的核心。Cumming 明确指出它在一步里同时完成残余 RCMC、SRC 和方位压缩。

9. 近似到底在哪里

很多人说 ω-k“严格”，但它并不是对所有几何都完全无近似。当前这套推导至少依赖下面几条：

(a) 双曲线距离历程

要求点目标的斜距可由

R(\eta)=\sqrt{R_0^2+V_r^2\eta^2}

很好描述，也就是直线匀速、标准 stripmap 几何。

(b) 驻相法

从 (f_\tau,\eta) 到 (f_\tau,f_\eta) 的解析结果使用了 POSP。

(c) 等效速度对距离不变

综述文献明确指出，ω-k 的经典形式假定 V_r 是距离不变的，所以它不适合特别大的距离向幅宽；Cumming 也指出 ω-k 仍有一个限制大 swath 的近似。

所以更准确地说，ω-k 是 在标准单基地 stripmap 几何下，对距离-方位耦合处理得比 RD/CS 更彻底的二维频域算法，而不是“所有条件下都无近似”的算法。

10. 最终得到的实现公式

工程上，经典 ω-k 可以概括为：

\boxed{ \begin{aligned} &\text{(1) } s(\tau,\eta) \;\xrightarrow{\mathcal F_{\tau,\eta}}\; S_{2D}(f_\tau,f_\eta)\\[2mm] &\text{(2) } S_{2D}\leftarrow S_{2D}\,H_{\text{ref}}(f_\tau,f_\eta)\\[2mm] &\text{(3) } f_0+f'_\tau = \sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}} \qquad\text{(Stolt)}\\[2mm] &\text{(4) } I(r,x) = \mathcal F^{-1}_{2D}\big\{\widetilde S(f'_\tau,f_\eta)\big\}. \end{aligned} }

如果换成波数变量，就是

k_y=\sqrt{k_r^2-k_x^2}

然后把数据重采样到均匀的 (k_x,k_y) 网格，再做二维 IFFT。这个写法最能体现“ω-k”名字的来历：它完全是在二维频率/波数域里工作的。

11. 总结

Omega-K 的严格推导核心只有一句：

先把点目标回波通过二维 FFT 写成带有

\sqrt{(f_0+f_\tau)^2-\frac{c^2f_\eta^2}{4V_r^2}}

这个耦合项的二维频谱；再通过参考函数乘法消掉参考距离处的相位；最后用 Stolt 映射把这个根号项改写成新的线性频率变量，使剩余相位线性化，于是 2D IFFT 就直接给出聚焦图像。

MATLAB代码

clc; clear; close all;

%% =========================================================
% 1) System parameters
% ==========================================================
c  = 3e8;
fc = 5.3e9;
lambda = c / fc; %#ok<NASGU>

v   = 150;
PRF = 600;
Naz = 512;

fs = 20e6;
Nr = 2048;
Tp = 10e-6;
B  = 15e6;
Kr = B / Tp;

Rref = 10e3;

%% =========================================================
% 2) Centered axes
%    Use FFT-consistent centered grids:
%      time  : ((0:N-1) - N/2) / Fs
%      freq  : ((0:N-1) - N/2) * (Fs/N)
% ==========================================================
eta   = ((0:Naz-1) - Naz/2) / PRF;
xPlat = v * eta;

tauRef = 2 * Rref / c;
tFast  = ((0:Nr-1) - Nr/2) / fs;
tr     = tauRef + tFast;

%% =========================================================
% 3) Point-target scene
%    [azimuth_x(m), range_offset_from_Rref(m), amplitude]
% ==========================================================
targets = [
   -20,   0,   1.0;
     0,  30,   0.8;
    20, -30,   0.7
];

%% =========================================================
% 4) Generate raw SAR data raw(Nr, Naz)
%    Fast time is referenced to tauRef, so the delay is dt = 2*(R-Rref)/c
% ==========================================================
raw = complex(zeros(Nr, Naz));

for ia = 1:Naz
    xp = xPlat(ia);
    sig = complex(zeros(Nr, 1));

    for k = 1:size(targets, 1)
        xt = targets(k, 1);
        dR = targets(k, 2);
        A  = targets(k, 3);

        R0 = Rref + dR;
        R  = sqrt(R0^2 + (xp - xt)^2);

        dt     = 2 * (R - Rref) / c;
        tShift = tFast - dt;
        gate   = abs(tShift) <= Tp / 2;

        echo = A .* gate .* ...
               exp(1j * pi * Kr * tShift.^2) .* ...
               exp(-1j * 4 * pi * fc * R / c);

        sig = sig + echo(:);
    end

    raw(:, ia) = sig;
end

raw = raw + 0.01 * (randn(size(raw)) + 1j * randn(size(raw)));

%% =========================================================
% 5) Range compression
%    raw/ref are centered in time, so FFT must be done with ifftshift first
% ==========================================================
ref = (abs(tFast) <= Tp / 2) .* exp(1j * pi * Kr * tFast.^2);
Hrc = conj(fft(ifftshift(ref(:)), Nr));

rawFFT = fft(ifftshift(raw, 1), Nr, 1);
rc = fftshift(ifft(rawFFT .* repmat(Hrc, 1, Naz), Nr, 1), 1);

%% =========================================================
% 6) Wavenumber axes
%    fr : range frequency after FFT
%    fa : azimuth Doppler frequency
%    KR : two-way range wavenumber
%    KX : azimuth spatial wavenumber
% ==========================================================
fr = ((0:Nr-1) - Nr/2) * (fs / Nr);
fa = ((0:Naz-1) - Naz/2) * (PRF / Naz);

kr = 4 * pi * (fc + fr) / c;
kx = 2 * pi * fa / v;

[KR, KX] = ndgrid(kr, kx);

%% =========================================================
% 7) 2D FFT to (KR, KX) domain
% ==========================================================
S2 = fftshift(fft(ifftshift(rc, 1), Nr, 1), 1);
S2 = fftshift(fft(ifftshift(S2, 2), Naz, 2), 2);

%% =========================================================
% 8) Bulk compression
%    KY is the range-direction spatial wavenumber after dispersion relation
%
%    Important:
%      Because fast time is referenced to Rref, the reference function is
%      exp(j * Rref * (KY - KR)), not only exp(j * Rref * KY).
% ==========================================================
valid = KR > abs(KX);

KY = zeros(size(KR));
KY(valid) = sqrt(KR(valid).^2 - KX(valid).^2);

S2(~valid) = 0;
Href = exp(1j * Rref * (KY - KR));
Sref = S2 .* Href;

%% =========================================================
% 9) Stolt interpolation
%    Map uniform KR samples to uniform KY samples
%    A Jacobian factor KY/KR is included for better amplitude fidelity
% ==========================================================
kyGrid = kr(:);
Sstolt = complex(zeros(size(Sref)));

for ia = 1:Naz
    kySrc = KY(:, ia);
    sSrc  = Sref(:, ia) .* (kySrc ./ kr(:));

    validCol = valid(:, ia) & isfinite(kySrc) & isfinite(sSrc) & (kySrc > 0);
    kyv = kySrc(validCol);
    sv  = sSrc(validCol);

    [kyv, uniqIdx] = unique(kyv, 'stable');
    sv = sv(uniqIdx);

    if numel(kyv) >= 2
        Sstolt(:, ia) = interp1(kyv, sv, kyGrid, 'linear', 0);
    end
end

%% =========================================================
% 10) 2D IFFT to image domain
% ==========================================================
img = fftshift(ifft(ifftshift(Sstolt, 1), Nr, 1), 1);
img = fftshift(ifft(ifftshift(img, 2), Naz, 2), 2);

imgAbs = abs(img);
imgdB  = 20 * log10(imgAbs / max(imgAbs(:)) + 1e-12);

%% =========================================================
% 11) Image axes
% ==========================================================
rangeAxis = Rref + tFast * c / 2;
azAxis    = xPlat;

%% =========================================================
% 12) Display results
% ==========================================================
figure;
imagesc(azAxis, rangeAxis, imgdB);
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('\omega-k SAR Image');
colorbar;
caxis([-40 0]);
hold on;
plot(targets(:,1), Rref + targets(:,2), 'rx', 'LineWidth', 1.5, 'MarkerSize', 10);
hold off;

figure;
imagesc(azAxis, rangeAxis, 20 * log10(abs(raw) / max(abs(raw(:))) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('Raw SAR Data (dB)');
colorbar;
caxis([-40 0]);

figure;
imagesc(azAxis, rangeAxis, 20 * log10(abs(rc) / max(abs(rc(:))) + 1e-12));
axis xy;
xlabel('Azimuth (m)');
ylabel('Slant Range (m)');
title('Range Compressed Data (dB)');
colorbar;
caxis([-40 0]);