SAR成像系列：【7】合成孔径雷达（SAR）成像算法-后向投影（Back Projecting）算法（附Matlab代码）

前面介绍了SAR成像的RD算法和CS算法，接下来介绍两种时域成像算法，其一就是后向投影（BP）算法.

BP成像的优点：成像算法简单，鲁棒性好，分辨率高，适用于任何轨道或飞行轨迹模型，不存在斜距近似假设，运动补偿容易。尤其适合双基和多基SAR成像。

BP成像的缺点：使用时域成像方式，计算量大，成像效率低，实时性差。

（1）BP成像算法原理

首先将场景进行网格划分获得各个分辨单元的网格坐标，然后根据时延大小将距离压缩后的回波数据反向投影到成像区域网格的每个像素上。BP成像流程如下所示。

 （2）算法核心

首先，成像场景事先已经划分网格，需要计算雷达发射的每个方位时刻脉冲与网格点的距离，并根据距离计算出双程时延。然后根据时延找到回波中相同时延的数据，即网格点在这一方位时刻收到的数据，这就是“后向投影”。
由于回波信号是数字采样的，采样率有限，这种对应关系存在误差，这就需要插值提高成像精度。

（3）算法说明

由于BP算法是利用时延和相位补偿来进行目标聚焦，在理论上的推导公式较少，所以理解起来比较困难（尽管成像原理非常简单），最好的方式就是结合成像的代码进行理解。

（4）算法仿真

仿真参数如下表所示：

载频	5.3GHz
雷达高度	800Km
斜视角	0°
带宽	30MHz
雷达速度	7000m/s
目标个数	4
脉冲持续时间	40us

仿真结果：

回波信号实部：

回波信号虚部：

 

 距离压缩结果：

 BP成像结果：

 Matlab部分代码：

clc;close all; clear all;
c=3e8;  %光速
%星载SAR
disp('Spaceborne Circle SAR')
%%雷达参数
fc=5.3e9;                                                               %载频 C波段
lamda=c/fc;                                                             %波长 
Br=30e6;                                                                %线性调频信号调频率带宽 30MHz
Tr=40e-6;                                                               %脉冲持续时间
Kr=Br/Tr;                                                               %调频率
Fr=1.2*Br;                                                              %快时间采样率  1.2倍带宽
H=800e3;                                                                %SAR的距离向距离 800km
theta=26.288/180*pi;                                                    %俯仰角θ(下视角)
Yc=H*tan(theta);
R0=H/cos(theta);
rouY = c/(Br*2);                                                        % 距离向分辨率
gama = 0*pi/180;                                                        %方位向斜视角  正侧视 side-looking
D=10;                                                                   %SAR 方位向天线长度
rouX = D/2;                                                             %方位向分辨率
V=7000;                                                                 %SAR速度

Ka=-2*V^2/lamda/R0;                                                    %doppler 调频率
Lsar=0.886*lamda*R0/D;                                                 %SAR 合成孔径长度
Tsar=Lsar/V; 
Ba=abs(Ka*Tsar);                                                        %doppler 调频带宽（方位带宽）
PRF=1.2*Ba;                                                                 %脉冲重复频率

%%目标区域
Xmin=-200;                                                                %目标区的方位向范围[Xmin,Xmax]
Xmax=200; 
%方位时间序列
PRT=1/PRF;                                                              %脉冲重复间隔
ds=PRT;                                                                 %慢时间域采样步长                               
Nslow=2048;                                      %慢时间域采样数
%Nslow=2^nextpow2(Nslow);                                                %写成2的指数 for fft
ta=-PRT*Nslow/2:PRT:PRT*(Nslow/2-1);                     %慢时间域离散时间序列

%距离时间序列
Y0=1000;  %地距向成像范围[Yc-Y0,Yc+Y0]                                                                      
dt=1/Fr;                                                               %快时间域采样步长                                                           %SAR 合成孔径时间
                                            
Rmin=sqrt(H^2+(Yc-Y0)^2);                                                %SAR与一个孔径测绘带的最短距离                                         
Rmax=sqrt(H^2+(Yc+Y0)^2);                                               %SAR与一个孔径测绘带的最长距离
Nfast=ceil(2*(Rmax-Rmin)/c/dt+Tr/dt);                                   %快时间域的采样数
tr=linspace(2*Rmin/c-Tr/2,2*Rmax/c+Tr/2,Nfast);                       %快时间域的离散时间序列
r=tr*c/2;%距离向的距离序列

%%采样分辨率
DR=c*dt/2;                                                              %距离向分辨率
DX=V*ds;                                                                 %方位向分辨率
%% 目标参数   极坐标
% Rpt=[0 20 20 20 20];                                                     %目标点的位置,径向位置
% thetaT=[0,0,0,0,0];                                                        %仰角
% alpha=[0,pi/4,3*pi/4,5*pi/4,7*pi/4];                                       %方位角

Rpt=[R0-5 R0+5 R0-5 R0+5];                                                     %目标点的位置,径向位置
Xpt=[-5 -5 5 5]; 
Spt=[10 20 10 20];

Npt=length(Rpt);                                                           %目标个数
%% 基线参数

K=Npt;                                                                   %目标个数
N=Nslow;                                                                %慢时间域向量长度（即采样个数）
M=Nfast;                                                                %快时间域向量长度

Srnm_sN=zeros(N,M,sN);  
    for k=1:1:K
        sigma=1;
        Dslow=ta*V-Xpt(k);                                                  %SAR与目标方位向距离形成的序列
        R=sqrt(Dslow.^2+(Rpt(k)-By(i))^2+(Bs(i)-Spt(k))^2);                         %生成径向距离序列    加入距离向基线
        tau=2*R/c;                                                          %生成径向距离延迟序列
        Dfast=ones(N,1)*tr-tau'*ones(1,M);                                  %生成t-tao矩阵
        phase=pi*Kr*Dfast.^2-(4*pi/lamda)*(R'*ones(1,M));                  %回波相位矩阵
        Srnm_sN(:,:,i)=Srnm_sN(:,:,i)+sigma*exp(j*phase).*(abs(Dfast)<=Tr/2).*((abs(Dslow)<=Lsar/2)'*ones(1,M));
    end

figure;
imagesc(imag(Srnm(:,:,1)));colormap(gray)
title('点目标原始回波信号虚部');
xlabel('距离向（采样点）');
ylabel('方位向（采样点）');
% colorbar;
figure;
imagesc(real(Srnm_sN(:,:,1)));
colormap(gray);
title('点目标原始回波信号实部');
xlabel('距离向（采样点）');
ylabel('方位向（采样点）');
colorbar;
for i=1:1:1
    Srnm=zeros(N,M);
    Srnm=Srnm_sN(:,:,i);
    %% BP成像算法 %======================================================================
    %BP成像算法
%========================================================================
    % 步骤一：距离向匹配滤波
    t_ref=tr-2*R0/c;                                             %距离向信号发射时刻生成的时间序列
    Refr=exp(j*pi*Kr*t_ref.^2).*(abs(t_ref)<=Tr/2);                 %参考信号 未取负共轭 在下面fft中进行   参考fft的性质
    Sr_temp=fftshift(fft((fftshift(Srnm.')))).'.*(ones(N,1)*conj(fftshift(fft((fftshift(Refr)))))); %距离向频域压缩
    nup=1;                                                                    %距离向升采样系数为50(平滑的作用，提高精度)
    M_up=M*nup;
    nz=M_up-M;                                                                 %补零个数
    delta_tr=1/((M_up-1)*Fr/(M-1));                                            %采样间隔（做了优化修正！）
%     sig_fr_up=[sig_fr(:,1:Nr/2),zeros(Na,nz),sig_fr(:,(Nr/2+1):Nr)];         %不用fftshift时的补零操作
    Srnm_up=[zeros(N,nz/2),Sr_temp(:,1:M/2),Sr_temp(:,(M/2+1):M),zeros(N,nz/2)];%两端补零
    Sr=fftshift(ifft(fftshift(Srnm_up.'))).';
    %显示匹配滤波结果
    figure;
    colormap(gray);colormap(gray)
    imagesc(abs(Sr));
    title('距离向匹配滤波后');
    xlabel('距离向（采样点）');
    ylabel('方位向（采样点）');
    colorbar;
  %=======================================================================
    % 步骤二：对成像区域进行网格化（N，M）
  %=======================================================================
    R=zeros(1,M);
    for ii=1:M
        R(1,ii)=Rmin+(Rmax-Rmin)/(M-1)*(ii-1);
    end
    X=zeros(1,N);
    for ii=1:N
        X(1,ii)=Xmin+(Xmax-Xmin)/(N-1)*(ii-1);
    end
    R=ones(N,1)*R;
    X=X'*ones(1,M);    %========================================================================
    % 步骤三：根据时延在回波域寻找相应的位置并进行成像    %========================================================================
    tic;
    f_back=zeros(N,M);
    pixel=zeros(N,M);
    str=['第',num2str(i),'个雷达方位后向投影计算中....'];
    h1=waitbar(0,str); 

由于代码是从SAR三维成像代码中修改，具体代码参考：(43条消息) SAR成像中后向投影（BP）算法Matlab代码；-行业报告文档类资源-CSDN文库