【SAR成像】SAR成像BP算法【附全部MATLAB代码】

微信公众号：EW Frontier

QQ交流群：949444104

主要内容

MATLAB代码

%BPA 多点目标 单站 SAR
%parameters from Table 6.1
%date: 2010/10/14
clc;
clear all;
close all;
% (1) parameters' definition
%================================================
c=3e8;
j=sqrt(-1);
pi=3.1416;
fc=5.3e9;
lamda=c/fc;
D=4; %antenna size in azimuth direction
Va=150;
Kr=20e12;
Tr=2.5e-6;
sq_ang=3.5/180*pi;
Br=Kr*Tr;
Frfactor=1.2;
Fr=Br*Frfactor;
Ba=0.886*2*Va*cos(sq_ang)/D;
Fafactor=1.2;
Fa=Ba*Fafactor;
R_near=2e4; %near range of the scene
R_far=R_near+1000; %far range of the scene
% R=(R_near+R_far)/2; %assumed a fixed R for implement
% Y_min=-50; %both R_near and R_far are closest approach range
% Y_max=Y_min+100;
La_near=0.886*R_near*lamda/cos(sq_ang)^2/D; %synthetic aperture length of near range
La_far=0.886*R_far*lamda/cos(sq_ang)^2/D; %synthetic aperture length of far range
Tc_near=-R_near*tan(sq_ang)/Va; %beam center cross time of near range
Tc_far=-R_far*tan(sq_ang)/Va; %beam center cross time of far range
fdc=2*Va*sin(sq_ang)/lamda; %doppler centroid
Y_min=Va*Tc_far;
Y_max=Y_min+100;
Rmin=sqrt(R_near^2+(Tc_near*Va+La_near/2)^2); %shortest range between radar and scene
Rmax=sqrt(R_far^2+(Tc_far*Va-La_far/2)^2); %longest range between radar and scene
disp('parameters:');
disp('minimal slant range:');disp(Rmin);
disp('maximal slant range:');disp(Rmax);
disp('range resolution:');disp(0.886*(c/2/Br));
disp('azimuth resolution:');disp(0.886*Va/Ba);
disp('doppler centroid frequency:');disp(fdc);
%================================================
% (2) echo model
%================================================
Nr=(2*Rmax/c-2*Rmin/c+Tr)*Fr;
Nr=2^nextpow2(Nr);
tr=linspace(-Tr/2+2*Rmin/c,Tr/2+2*Rmax/c,Nr);
Fr=(Nr-1)/(Tr/2+2*Rmax/c-(-Tr/2+2*Rmin/c));
Na=((Tc_near+La_near/2/Va)-(Tc_far-La_far/2/Va))*Fa;
Na=2^nextpow2(Na);
ta=linspace(Tc_far-La_far/2/Va,Tc_near+La_near/2/Va,Na);
Fa=(Na-1)/(Tc_near+La_near/2/Va-(Tc_far-La_far/2/Va));
Rpt=[R_near R_near+500 R_near+1000]; %position of point targets
Ypt=[0 0 0]; %Rpt is there closest approach range
La=0.886*Rpt*lamda/(cos(sq_ang)^2)/D; %synthetic aperture length of each target
Tc=-Rpt*tan(sq_ang)/Va; %beam center cross time of each target
Npt=length(Rpt);
Y_high=max(Ypt)+50; %%确定成像网格的范围以使其包括目标点;
Y_low=min(Ypt)-50; %%成像区域是R_right-R_left*(Y_high-Y_low)的一片区域
R_left=R_near-50;
R_right=R_far+50;
disp('number of point targets:');disp(Npt);
disp('range sample number:');disp(Nr);
disp('azimuth sample number:');disp(Na);
disp('range sample rate:');disp(Fr);
disp('azimuth sample rate:');disp(Fa);
sig=zeros(Na,Nr);
for k=1:Npt
delay=2/c*sqrt(Rpt(k)^2+(Ypt(k)-ta*Va).^2);
Dr=ones(Na,1)*tr-delay'*ones(1,Nr);
sig=sig+exp(j*pi*Kr*Dr.^2-j*2*pi*fc*delay'*ones(1,Nr))...
.*(abs((ta-Ypt(k)/Va-Tc(k))'*ones(1,Nr))<=La(k)/2/Va)...
.*(abs(Dr)<=Tr/2);
end
figure('Name','回波信号幅度');imagesc(abs(sig));
title('回波信号幅度');xlabel('距离向(采样点)');ylabel('方位向采样点');
figure('Name','回波信号相位');imagesc(abs(angle(sig)));
title('回波信号相位');xlabel('距离向(采样点)');ylabel('方位向采样点');
%================================================
%%--BP 算法
%%步骤一距离向匹配滤波
sig_rd=fft(sig,[],2);
fr=-1/2:1/Nr:(1/2-1/Nr);
fr=fftshift(fr*Fr);
filter_r=ones(Na,1)*exp(j*pi*fr.^2/Kr);
sig_rd=sig_rd.*filter_r;
nup=100; %距离向升采样系数为 100
Nr_up=Nr*nup;
nz=Nr_up-Nr;
dtr=1/nup/Fr;
sig_rd_up=zeros(Na,Nr_up);
sig_rd_up=[sig_rd(:,1:Nr/2),zeros(Na,nz),sig_rd(:,(Nr/2+1):Nr)];
sig_rdt=ifft(sig_rd_up,[],2);
% figure('Name','距离向匹配滤波后');imagesc(abs(sig_rdt));
% title('距离向匹配滤波后');xlabel('距离向(采样点)');ylabel('方位向采样点');
%%分别压至(2*sqrt(R_near^2+(Va*Tc_near)^2)/c-(2*Rmin/c-Tr/2))*Fr
%%----(2*sqrt((R_near+500)^2+(Va*Tc(2))^2)/c-(2*Rmin/c-Tr/2))*Fr
%%----(2*sqrt((R_near+1000)^2+(Va*Tc(3))^2)/c-(2*Rmin/c-Tr/2))*Fr
%%--步骤二 对成像区域进行网格化Na*Nr
R=zeros(1,Nr);
for ii=1:Nr
% R(1,ii)=R_near+(R_far-R_near)/(Nr-1)*(ii-1);
R(1,ii)=R_left+(R_right-R_left)/(Nr-1)*(ii-1);
end
Y=zeros(1,Na);
for ii=1:Na
Y(1,ii)=Y_low+(Y_high-Y_low)/(Na-1)*(ii-1);
end
R=ones(Na,1)*R;
Y=Y'*ones(1,Nr);
%%--步骤三根据时延在回波域寻找相应位置并进行成像
f_back=zeros(Na,Nr);
for ii=1:Na
% ii=100;
R_ij=sqrt(R.^2+(Y-Va*ta(ii)).^2);
t_ij=2*R_ij/c;
% t_ij=round((t_ij-tr(Nr/2+1))/dtr)+Nr_up/2+1;%%怎么将时域转化为点数
t_ij=round((t_ij-(2*Rmin/c-Tr/2))/dtr);
% t_ij=round(t_ij/dtr);
it_ij=(t_ij>0&t_ij<=Nr_up);
t_ij=t_ij.*it_ij+Nr_up*(1-it_ij); %%将矩阵中为零的元素用 Nr_up 代替
sig_rdta=sig_rdt(ii,:);
sig_rdta(Nr_up)=0;
f_back=f_back+sig_rdta(t_ij).*exp(1i*4*pi*R_ij/lamda);
end
figure('Name','BP 算法处理后');imagesc(abs(f_back));
title('BP 算法处理后');xlabel('距离向');ylabel('方位向');