SAR成像系列：【3】合成孔径雷达（SAR）的二维回波信号与简单距离多普勒（RD）算法 （附matlab代码）

合成孔径雷达发射信号以线性调频信号（LFM）为基础，目前大部分合成孔径雷达都是LFM体制，为了减轻雷达重量也采用线性调频连续波（FMCW）体制；为了获得大带宽亦采用线性调频步进频（FMSF）体制。

（1）LFM信号

LFM的主要特点在于可以使载波的瞬时频率随调制信号的变化而变化，当其频率线性增加时，称为正调频；当其频率线性减少时，称为负调频。LFM信号的幅度频谱存在部分起伏现象，这是由菲涅尔积分造成的；信号的频谱并不完全限制在-B/2~B/2之内，随着时宽带宽积的增大，信号的幅频特性越接近矩形，顶部起伏也会减小。LFM解决了探测距离和分辨率之间的矛盾，在雷达和制导武器上得到广泛应用。LFM的时域表示为：

其中，T为时宽，K为调频率，rect为矩形窗。

求LFM信号的频谱需要对其作FFT，得到：

存在指数二次项积分，需采用驻定相位原理（POSP）求解，驻定相位点为：

因此求得频谱为

LFM的时域图和频谱图如下（附matlab代码）。

 

 %%%%%%%%% 线性调频信号LFM%%%%%%%%
%% LFM时域波形
clc; clear; close all;
f0 = 0;    %雷达中心频率
T = 3e-7;  %脉宽
B = 3e8;    %带宽
fs = 2*B; %采样率
Ts = 1/fs; %采样时间
N = T/Ts; %采样数
k = B/T;   %调频率
t = linspace(-T/2,T/2,N);%时间
y = exp(1j*(2*pi*f0*t + pi*k*t.^2));%LFM信号
figure;
plot(t*1e6,real(y));xlabel('时间（us）');ylabel('幅度');
title('LFM信号时域波形（实部）');
figure;
plot(t*1e6,imag(y));xlabel('时间（us）');ylabel('幅度');
title('LFM信号时域波形（虚部）');
grid on; axis tight;
%% LFM频谱图
S = fftshift(fft(y));     %频谱
f = linspace(-fs/2,fs/2,N);%频率轴
figure;
plot(f*1e-6,abs(S)./max(max(abs(S))));
xlabel('频率（MHz）')
ylabel('归一化频谱幅度');
title('LFM信号频谱');

 （2）快时间与慢时间

在合成孔径雷达成像中，快时间和慢时间是一个相对的概念。在工程上，快时间指的是脉内时间变化，慢时间是脉间时间变化。

雷达发射信号是以脉冲的形式发射的，发射频率称为脉冲重复频率（PRF），PRF的设定是根据雷达的功能和性能确定的。从几K到几十K，甚至几百K。一般星载SAR的PRF为几K，每一个脉冲都有一个时间戳，这个连续的时间戳合起来叫做慢时间。如下图t1到t6为慢时间标记。

 雷达发射脉冲信号的周期为T，T=1/PRF。雷达有效信号能量时间占一个雷达信号周期的比例称为占空比，也就是说发射的有效信号不总是充满整个雷达信号周期。周期内的有效能量部分的时间称为快时间，一般用来表示。

这种慢时间和快时间在主观上是交替进行的，慢时间沿方位向播放，快时间沿距离向播放，这就形成了SAR的“停-走-停”成像模式。

（3）SAR的二维回波信号

在低轨星载和机载SAR使用中，回波是以“停-走-停”的方式录取的，它是以快时间轴和慢时间轴构成的二维平面。如下图所示。

 合成孔径雷达一般发射的是LFM信号,快时间的数学表达式为

fc为载频，Kr为调频率，rect（.）为脉冲包络。假设，一个散射系数为A0的点目标在雷达发射波束内，距离雷达的距离为R，则雷达接受到的回波是发射信号与目标散射系数的卷积，因此接收信号为(忽略后向散射引起的相位变化):

上式看起来是一个距离向一维回波，但是它的慢时间包含在R中。由于录取平台是运动的，因此R会随着平台位置的变化而改变，如下图所示。

 慢时间在t1-t11时刻距离目标的R1-R11一直发生变化。R关于慢时间的表达式为(正侧式SAR为例)：

R0为雷达到目标的最短参考距离。将R(t)带入回波中，得到SAR的二维回波信号：

为方位信号包络。对Rt进行泰勒展开取近似得到

带入到二维回波信号中，整理得到：

其中，Ka为方位向多普勒调频率：

由回波信号可以看到，在方位向和距离向分别存在两个线性调频信号。

（4）正侧视距离多普勒（RD）算法

一个简单的距离多普勒算法包括距离压缩、距离徙动矫正、方位压缩。其中距离压缩和方位压缩通过LFM匹配滤波或去斜处理实现。距离徙动矫正是补偿掉不同方位位置回波引起的Rt的变化，若距离徙动远小于距离分辨率，则不需要进行距离徙动矫正，RD成像仅剩下距离压缩和方位压缩，且这两步操作没有先后之分。

把二维回波信号分解，仅考虑相位项，距离向的对应LFM信号为

方位向的对应LFM信号为

因此，简单RD算法的算法流程如下图所示：

 算法中的距离参考信号为

其中为距离包络，为距离向参考时间（快时间）。

方位参考信号为

经过距离压缩和方位压缩，得到简单RD算法的成像效果。5个点目标未进行RCM的matlab代码如下：

%% RD算法   含距离徙动矫正（最近邻插值和sinc插值）
%%%
%%%Authed  by Piaobo 氵茶花彡
clear;close all;clc;
SNR = -15;                          % 信噪比

c=3e8;
f0 = 9.875e9;                     % 雷达工作频率Hz
lamda = c/f0;                 % 雷达工作波长m
H = 1000;                           % 高度
Yc=2000;                         % 成像区域中线
R0 = sqrt(Yc^2+H^2);       % 中心斜距m

theta = asind(H/R0);          % 下视角
Br=50e6;                           % 带宽
Vr = 200;                           % 雷达有效速度m/s
Tr =5e-6;                          % 脉冲持续时间s
Kr = Br/Tr;                           % 线性调频率
Fr = 1.2*Br;                          % 距离采样频率，1.2为过采样率
Ts = 1/Fr;                             % 距离采样时间间隔s

Nk = ceil((2 * 800/ c + Tr) / Ts);  %距离向前后500m
Nf = 2^nextpow2(Nk);                % 距离向的采样点个数
tf_ori = [-Nf/2:1:Nf/2-1]*Ts;                     % 距离向采样时序
tf = [-Nf/2:1:Nf/2-1]*Ts+2*R0/c;                  % 实际快时间采样值

La = 6;                                                % 等效天线尺寸
Ls = R0*lamda/La;                   % 合成孔径时长度m,Ls=（0.886*R0*lamda）/(La*cos(Theta)）
% Ta = Ls/Vr;                                         % 目标照射时间s
Ta = 0.8;                                         % 目标照射时间s
Ls =Ta*Vr;
Ka = -2 * Vr^2 / (lamda * R0);                          % 方位多普勒调频率Hz
Ba=abs(Ka*Ta);                                             % 多普勒频率调制带宽
PRF = ceil(1.8*Ba);                                       % 方位采样率Hz
% PRF = 1000;                                       % 方位采样率Hz
PRT = 1/PRF;                                      % 方位向采样时间间隔s
Ns = 2^nextpow2((80/Vr+Ta)*PRF);             % 方位向的采样点个数 左右各100m
ts = [-Ns/2 : (Ns/2 - 1)] * PRT;                         % 方位向采样时序
% 理论分辨率
rho_r=c/2/Br;
% rho_a=Vr*PRT;
rho_a=La/2;

% 目标参数
X0 = [-20 20 0 -20 20];                         % 目标1位置坐标
Z0 = [0 0 0 0 0];
Y0 = [Yc+50 Yc+50 Yc Yc-50 Yc-50];
NT=size(X0,2);

%%================================================================
%%生成回波信号
Sb = zeros(Ns,Nf);
sigma = 1; % 回波幅度
for ii=1:NT
    R = sqrt((Vr*ts-X0(ii)).^2+Y0(ii).^2+(Z0(ii)-H).^2);
    tau = 2*R/c;
    Dfast = ones(Ns,1) * tf - tau' * ones(1, Nf);
    phase = pi*Kr*Dfast.^2 - (2 * pi *f0 * tau') * ones(1,Nf);                                           
    Sb = Sb+sigma * exp(1j*phase) .* (abs(Dfast) <= Tr/2) .* ((abs(ts * Vr-X0(ii)) <=Ls/2)' * ones(1,Nf));
end
% Sb = awgn(Sb,SNR,0);                                % 回波加噪

figure
imagesc(real(Sb)),colormap(gray);

%% 距离向压缩
x0 = ones(Ns,1)*(exp(-1j*pi*Kr*(tf_ori).^2).* (abs(tf_ori) <= Tr/2)); % 距离向匹配函数
fftx1 = fftshift(fft(fftshift(x0.'))).'; % 距离向匹配函数FFT
fftSb = fftshift(fft(fftshift(Sb.'))).'; % 原信号FFT
y0 = fftshift(ifft(fftshift((fftSb.*fftx1).'))).';  % 距离向压缩后信号
%显示
ta = ts * Vr;                                     % 方位向的距离序列
tr = tf * c / 2;                                  % 快时间采样对应的距离域(单程距)
% figure;imagesc(tr,ta,abs(y0));colormap(gray);
% axis([tr(Nf/2-Nf/2^8),tr(Nf/2+Nf/2^8),ta(1),ta(end)]);
% xlabel('距离单元(m)');ylabel('方位单元(m)');title('距离向压缩');
figure;imagesc(abs(y0));colormap(gray);
xlabel('距离向');ylabel('方位向');title('距离向压缩');
 %% 距离徙动矫正1    SINC插值
 

%% 距离徙动矫正2    最近邻插值
 

%% 方位向压缩
ffty0 = fftshift(fft(fftshift(y0))); % 距离向压缩后信号FFT
% x1 = exp(-1j*4*pi/lamda*sqrt((X0-Vr*ts-Bx(1)).^2+(Y0)^2+(Z0-Bz(1)-H)^2))'*ones(1,Nf); % 方位向匹配函数
x1 = exp(1j*pi*(ts.^2)*Ka)'*ones(1,Nf); % 方位向匹配函数
fftx1 = fftshift(fft(fftshift(x1))); % 方位向匹配函数FFT
y1 = fftshift(ifft(fftshift((ffty0.*fftx1)))); % 方位向压缩后信号

figure;imagesc(abs(y1));colormap(gray);
xlabel('距离向');ylabel('方位向');title('方位向压缩');
% xlabel('Range(m)');ylabel('Azimuth(m)');title('Uncompensated');
成像结果如下：

原始信号实部：

原始信号虚部：

距离压缩结果：

方位压缩成像结果(此图未进行距离徙动补偿，非中心目标出现散焦)：

 中心目标的db图（边上有其他目标的扩散能量进来）：