雷达成像的三种不同算法，分别为CS算法、RD算法、RMA算法

%%产生Stripmap SAR的回波 clear all clc thetaT=0;%T平台波束斜视角 thetaT=thetaT*pi/180;%rad thetaR=0;%R平台波束斜视角 thetaR=thetaR*pi/180; c=3e8;%光速 fc=1.5e9;%载频 lambda=c/fc;%波长 %%测绘带区域 X0=200;%方位向[-X0,X0] Rtc=3000; Rrc=3000; Rc=(Rtc+Rrc)/2; R0=150;%距离向[Rc-R0,Rc+R0] %%距离向(Range),r/t domain Tr=1.33e-6;%LFM信号脉宽1.33us (200m) Br=150e6; %LFM信号带宽 150MHz Kr=Br/Tr; %调频斜率 Nr=1024; r=Rc+linspace(-R0,R0,Nr); t=2*r/c;%t域序列 dt=R0*4/c/Nr;%采样周期 f=linspace(-1/2/dt,1/2/dt,Nr);%f域序列 %%方位向(Azimuth,Cross-Range),x/u domain v=100;%SAR 平台速度 Lsar=300;%合成孔径长度 Na=512; x=linspace(-X0,X0,Na);%u域序列 u=x/v; du=2*X0/v/Na; fu=linspace(-1/2/du,1/2/du,Na);%fu域序列 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;%Doppler调频中心频率 fdr=ftdr+frdr;%Doppler调频斜率 %％目标位置 Ntar=3;%目标个数 Ptar=[Rrc,0,1 %距离向坐标,方位向坐标,sigma 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; %方位向fft s_kt=fftshift(fft(fftshift(s_ut).')).'; %CS变换 kc=4*pi/lambda; kc=kc*ones(1,Na); kx=fu/v; p_kx0=-sqrt(kc.^2-kx.^2);%相位项泰勒展开的系数函数 p_kx1=2*kc/c/p_kx0; p_kx2=-2.*kx.^2/c^2./p_kx0.^3; C_kx=-(c*p_kx1/2+1); Ks_r=1-2*Kr*Rc.*p_kx2; Ks_kx_r=Kr/pi./Ks_r; r0=Rc; s2_ut=exp(j*pi*C_kx.*ones(Nr,1)*Ks_kx_r.*(t'*ones(1,Na)-2*r0*(1+C_kx)/c).^2);%设计的线性调频信号 S_cs=s_kt.*s2_ut; %距离向fft S_kw=fftshift(fft(fftshift(S_cs))); %距离向匹配滤波 w=2*pi*f; rmc_r=exp(j.*w*2*C_kx*r0/c).*exp(j.*w.^2/4/pi/Kr/(1+C_kx)); rmc_r=rmc_r'*ones(1,Na); S_rmc=S_kw.*rmc_r; %距离向ifft S_kt=fftshift(ifft(fftshift(S_rmc))); d_kxr=4*pi/c^2*Kr*C_kx*(1+C_kx).*(Rc-r0).^2;%CS变换带来的相位误差 S_kt=S_kt.*exp(-j*d_kxr);%消除相位误差 %方位向匹配滤波 FU=ones(Nr,1)*fu; H_kx=exp(j*pi/fdr*(FU-fdc).^2);%方位向压缩因子 I_ut=S_kt.*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;%显示动态范围40dB 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(S_rmc)+1e-6); gm=max(max(G)); gn=gm-40;%显示动态范围40dB G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(x,r-Rc,-G),colormap(gray) grid on,axis tight, xlabel('Azimuth') ylabel('Range') title('(b)距离向匹配滤波后频谱') subplot(223) G=20*log10(abs(S_kt)+1e-6); gm=max(max(G)); gn=gm-40;%显示动态范围40dB G=255/(gm-gn)*(G-gn).*(G>gn); imagesc(x,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;%显示动态范围40dB 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)目标图象')