利用BP算法进行SAR成像_时域回波数据成像_matlab源码

clc close all clear all c=3e8; j=sqrt(-1); pi=3.1416; fc=5.3e9; lamda=c/fc; D=4; 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; R_far=R_near+1000; La_near=0.886*R_near*lamda/cos(sq_ang)^2/D; La_far=0.886*R_far*lamda/cos(sq_ang)^2/D; Tc_near=-R_near*tan(sq_ang)/Va; Tc_far=-R_far*tan(sq_ang)/Va; fdc=2*Va*sin(sq_ang)/lamda; Y_min=Va*Tc_far; Y_max=Y_min+100; Rmin=sqrt(R_near^2+(Tc_near*Va+La_near/2)^2); Rmax=sqrt(R_far^2+(Tc_far*Va-La_far/2)^2); disp('parameters:'); disp('minimal slant range:');disp(Rmin); disp('maximal slant range:');disp(Rmax); disp('range resolution:');disp(c/2/Br); disp('azimuth resolution:');disp(0.886*Va/Ba); disp('doppler centroid frequency:');disp(fdc); %%%%%%%%%echo 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]; Ypt=[0 0 0]; La=0.886*Rpt*lamda/(cos(sq_ang)^2)/D; Tc=-Rpt*tan(sq_ang)/Va; Npt=length(Rpt); Y_high=max(Ypt)+50; Y_low=min(Ypt)-50; R_left=R_near-50; R_right=R_far+50; disp('number of point target:');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/Na)... .*(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; Nr_up=Nr*nup; nz=Nr_up-Nr; dtr=1/nup/Fr; % sig_rd_up=zero(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('方位向（采样点）'); %%%%%%对成像区域进行网格化 R=zeros(1,Nr); for ii=1:Nr 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 R_ij=sqrt(R.^2+(Y-Va*ta(ii)).^2); t_ij=2*R_ij/c; t_ij=round((t_ij-(2*Rmin/c-Tr/2))/dtr); it_ij=(t_ij>0&t_ij<=Nr_up); t_jj=t_ij.*it_ij+Nr_up*(1-it_ij); 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('方位向（采样点）');