Matlab实现雷达成像RMA算法

clear; close all; clc ; thetaT=0; thetaT=thetaT*pi/180; thetaR=0; thetaR=thetaR*pi/180; c=3e8; fc=1.5e9; lambda=c/fc; X0=200; Rtc=3000; Rrc=3000; Rc=(Rtc+Rrc)/2; R0=150; Tr=1.33e-6; Br=150e6; Kr=Br/Tr; Nr=1024; r=Rc+linspace(-R0,R0,Nr); t=2*r/c; dt=R0*4/c/Nr; f=linspace(-1/2/dt,1/2/dt,Nr); v=100; Lsar=300; Na=512; x=linspace(-X0,X0,Na); u=x/v; du=2*X0/v/Na; fu=linspace(-1/2/du,1/2/du,Na); ftdc=v*sin(thetaT); ftdr=-(v*cos(thetaT))^2/lambda/Rtc; frdc=v*sin(thetaR); frdr=-(v*cos(thetaR))^2/lambda/Rrc; fdc=ftdc+frdc; fdr=ftdr+frdr; Ntar=3; Ptar=[Rrc,0,1 Rrc+50,-50,1 Rrc+50,50,1]; s_ut=zeros(Nr,Na); U=ones(Nr,1)*u; T=t'*ones(1,Na); for i=1:1:Ntar rn=Ptar(i,1);xn=Ptar(i,2);sigma=Ptar(i,3); rtn=rn+Rtc-Rrc; RT=sqrt(rtn^2+(rtn*tan(thetaT)+xn-v*U).^2); RR=sqrt(rn^2+(rn*tan(thetaT)+xn-v*U).^2); R=RT+RR; DT=T-R/c; phase=-pi*Kr*DT.^2-2*pi/lambda*R; s_ut=s_ut+sigma*exp(j*phase).*(abs(DT)<Tr/2).*(abs(v*U-xn)<Lsar/2); end; p0_t=exp(-j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2); p0_f=fftshift(fft(fftshift(p0_t))); s_uf=fftshift(fft(fftshift(s_ut))); src_uf=s_uf.*(conj(p0_f).'*ones(1,Na)); src_ut=fftshift(ifft(fftshift(src_uf))); src_kxf=fftshift(fft(fftshift(src_uf.'))).'; w=2*pi*f; w_c=2*pi*fc; kc=2*w_c/c; kx=fu/v; w_c=w_c*ones(1,Nr); W_Wc=(w+w_c)'*ones(1,Na); kr=sqrt(4*W_Wc.^2/c^2-(ones(Nr,1)*kx).^2)-kc; W=w'*ones(1,Na); W=c*sqrt((ones(Nr,1)*kx).^2+(kr+kc).^2)/2-2*pi*fc; src_kxkr=src_kxf; src_fut=fftshift(ifft(fftshift(src_kxkr))); FU=ones(Nr,1)*fu; H_kx=exp(j*pi/fdr*(FU-fdc).^2); I_ut=src_fut.*H_kx; I_ut=fftshift(ifft(fftshift(I_ut.'))).'; subplot(221) G=20*log10(abs(s_ut)+1e-6); gm=max(max(G)); gn=gm-40; G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(x,r-Rc,-G),colormap(gray) grid on,axis tight, xlabel('Azimuth') ylabel('Range') title('(a)成像结果') subplot(222) G=20*log10(abs(src_ut)+1e-6); gm=max(max(G)); gn=gm-40; G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(fu,r-Rc,-G),colormap(gray) grid on,axis tight, xlabel('Azimuth') ylabel('Range') title('(b)成像结果') subplot(223) G=20*log10(abs(src_kxkr)+1e-6); gm=max(max(G)); gn=gm-40; G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(fu,r-Rc,-G),colormap(gray) grid on,axis tight, xlabel('Azimuth') ylabel('Range') title('(c)成像结果') subplot(224) G=20*log10(abs(I_ut)+1e-6); gm=max(max(G)); gn=gm-60; G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(x,r-Rc,G),colormap(gray) grid on,axis tight, xlabel('Azimuth') ylabel('Range') title('(d)成像结果')