matlab-SAR回波成像,分别采用CSA处理,RDA处理以及RMA处理-源码

clc; clear; close all; warning off; addpath(genpath(pwd)); %% 雷达发射参数设置 a_resol =5; %方位向分辨率 r_resol =5; %距离向分辨率 %测绘带 W_y = 320; W_x = 520; %景中心斜距 H = 0; %雷达高度 Yc = 10000; %场景中心与雷达的距离 Yw = 200; %target area in range is within[Yc-Yw,Yc+Yw] Rc = sqrt(H^2+Yc^2); %R0 %雷达平台飞行速度 V=200; %LFM参数 fc = 5.0e9; %中心频率 C = 3e8; %光速 lambda = C/fc; %波长 a_os = 1.5; %方位向过采样因子 r_os = 1.5; %距离向过采样因子 %-----------目标参数-----------% Xmin =-150; %target area in azimuth is within[Xmin,Xmax] Xmax=150; %天线孔径设置 D= 2*a_resol; %天线真实孔径长度 theta = 0.886*lambda/D; Q = -1*pi/180; %前斜视角 %Q = 5*pi/180; %前斜视角 La = Rc*theta/cos(Q); %合成孔径长度 %-----------慢时间参数设置-----------% Ts = (La+W_y)/V; %方位向时宽 Ka = 2*V^2*(cos(Q)^3)/lambda/Rc; %doppler frequency modulation rate Ba = Ts*Ka; %多普勒带宽 PRF = a_os*Ba; %脉冲重复频率 PRT = 1/PRF; %脉冲重复时间 dxt = PRT; %sample spacing in slow-time domain Na = ceil(Ts/dxt); %sample number in slow-time domain Na = 2^nextpow2(Na); %for fft x = linspace((-W_y-La)/2/V,(W_y+La)/2/V,Na); %discrete time array in slow-time domain Azimuth = V*x; PRT = Ts/Na; %refresh dxt = PRT; PRF = 1/PRT; %-----------快时间参数设置-----------% Tr = 20e-6; %pulse width 20us Kr = 0.886*C/2/r_resol/Tr; %range FM rate Br = Kr*Tr; %range bandwidth Fs = r_os*Br; %距离采样率 drt = 1/Fs; %sample spacing in fast-time domain Rm = W_x+C*Tr/2; %range swath Nr = ceil(2*Rm/C/drt); Nr = 2^nextpow2(Nr); %块时间域采样点数 t = linspace(2*Rc/C-Tr/2,2*(Rc+W_x)/C+Tr/2,Nr); %discrete time array in fast-time domain Range = C*t/2; drt = (2*(Rc+W_x)/C+Tr-2*Rc/C)/Nr; Fs = 1/drt; %-----------目标参数设置-----------% Ntarget = 21;%目标个数 % dr=50; %load 'Ptarget.mat'; Ptarget=[ -100,Yc-200,1 -100,Yc-133.3,1 -100,Yc-66.7,1 -100,Yc,1 -100,Yc+66.7,1 -100,Yc+133.3,1 -100,Yc+200,1 0,Yc-200,1 0,Yc-133.3,1 0,Yc-66.7,1 0,Yc,1 0,Yc+66.7,1 0,Yc+133.3,1 0,Yc+200,1 100,Yc-200,1 100,Yc-133.3,1 100,Yc-66.7,1 100,Yc,1 100,Yc+66.7,1 100,Yc+133.3,1 100,Yc+200,1 ]; %目标位置 xc = -Ptarget(:,2)*(tan(Q))/V+(Ptarget(:,1))/V; Xc = xc*V;%方位向 fac = 2*V*sin(Q)/lambda;%多普勒频率 fan = ceil(fac/PRF*Na); %% 生成SAR回波 srt = zeros(Na,Nr); h = waitbar(0,'SAR回波生成'); for k = 1:Ntarget R = sqrt(Ptarget(k,2)^2+(Azimuth-Ptarget(k,1)).^2); Delay = 2*R/C; %range cell migration delay = ones(Na,1)*t-Delay.'*ones(1,Nr); phase = 1j*pi*Kr*delay.^2-1j*4*pi*fc*(R.'*ones(1,Nr))/C; % srt=srt+Ptarget(k,3)*rectpuls(delay,Tr).*(((sinc(0.886*atan(v*(x-xc(k))/R0)/theta)).^2).'*ones(1,M)).*exp(phase);%sinc窗 srt=srt+Ptarget(k,2)*rectpuls(delay/Tr).*rectpuls((Azimuth-Ptarget(k,1)).'*ones(1,Nr)/La).*exp(phase);%矩形窗 waitbar(k/Ntarget); end close(h); figure(1); % imagesc(abs(s0)); imagesc(Range*1e-3,Azimuth,abs(srt)); axis tight; xlabel('距离向（km）\rightarrow'); ylabel('\leftarrow方位向（m）'); clear R Dealy delay phase; %% 距离压缩 fr=(-Nr/2:Nr/2-1)*Fs/Nr; %range frequency sequence % fa=(-N/2:N/2-1)*PRF/N; fa=((0:Na-1).'-fix(Na/2))/Na*PRF+fac; fa=circshift(fa,round(+fac*Na/PRF)).'; % 方法一(考虑SRC的距离向匹配滤波) mt=zeros(1,Nr); t1=linspace(-Tr/2,Tr/2,ceil(Tr/drt)); mt1=exp(-1j*pi*Kr*t1.^2); N0=ceil((Nr-length(mt1))/2); mt(N0:N0+length(mt1)-1)=mt1; mf=fty(mt); srt=fty(srt); srt=ifty(srt.*(ones(Na,1)*mf)); % % 方法二 % mf=exp(1j*pi*fr.^2/Kr); %generate range matched filter % s0=ifty(fty(s0).*(ones(N,1)*mf));%range compression % 绘制距离压缩后的信号 figure(2); imagesc(Range*1e-3,Azimuth,abs(srt)); axis tight; xlabel('距离向（km）\rightarrow'); ylabel('\leftarrow方位向（m）'); title('距离压缩后的信号'); %% 方位向fft-距离多普勒RD域 srt=ftx(srt); srt=fty(srt); % 二次距离压缩(二维频谱域) D_fa_Vr=sqrt(1-lambda^2*fa.^2/4/V^2); Ksrc=(2*V^2*fc^3*D_fa_Vr.^3/C./fa.^2).'*(1./Range); h=waitbar(0,'二次距离压缩'); for k=1:Na Hsrc=exp(-1j*pi.*fr.^2./Ksrc(k,:)); srt(k,:)=srt(k,:).*Hsrc; waitbar(k/Na); end close(h); srt=ifty(srt); % 压缩后的距离多普勒域信号 figure(3); imagesc(abs(srt)); axis tight; xlabel('距离向（采样点数）\rightarrow'); ylabel('\leftarrow方位向（采样点数）'); title('方位fft信号-距离多普勒域'); %% 距离迁徙校正RCMC %RCM=lambda^2/8/v^2/(cos(theta_rc)^3)*(fa.^2).'*r; RCM=(1./D_fa_Vr-1).'*Range; % RCM=ones(N,1)*r./(D_fa_Vr'*ones(1,M))-ones(N,1)*r/D_fa_Vr(fan+N/2); RCM_point=RCM*2/C*Fs; % RCM_point=RCM*2*Fs/c; N_interp=8; Src_RCMC=zeros(Na,Nr); h=waitbar(0,'sinc插值，请稍后...'); for k=1:Na n_rcm=round(((1:Nr)+RCM_point(k,:)-1)*N_interp+1); N_add=N_interp*ceil(max(RCM_point(k,:)))+100; Src_interp=zeros(1,Nr*N_interp+N_add); Src_interp(1:Nr*N_interp)=interp(srt(k,:),N_interp); Src_RCMC(k,:)=Src_interp(n_rcm); waitbar(k/Na); end close(h); %----------RCMC后的距离多普勒域信号----------% figure(4); imagesc(abs(Src_RCMC)); axis tight; xlabel('距离向（采样点数）\rightarrow'); ylabel('\leftarrow方位向（采样点数）'); title('RCMC-距离多普勒域'); %% 方位向压缩 mfa=exp((1j*4*pi*fc/C*D_fa_Vr).'*Range); % mfa=exp(-1j*pi*lambda/(cos(theta_rc)^3)/2/v^2*(fa.^2)'*r); % mfa=circshift(mfa,0); sac=iftx(Src_RCMC.*mfa); %----------RDA----------% % figure(5); % imagesc(r*1e-3,X,abs(sac)); Z2=20*log10(abs(sac)+1e-6); Zm2=max(max(Z2)); Zn2=Zm2-40; %显示动态范围40dB Z2=(Zm2-Zn2)*(Z2-Zn2).*(Z2>Zn2); figure(5); imagesc(Range*1e-3,Azimuth,-Z2); colormap(gray); %绘制灰度图 axis tight; xlabel('距离向（km）\rightarrow'); ylabel('\leftarrow方位向（m）'); title('RDA最终信号'); %% 五个点目标分析和处理 %----------A点距离剖面、方位剖面----------% [location_a,location_r]=find(abs(sac)==max(max(abs(sac)))); Ns=8; %8点插值 sac_interp_r=interp(sac(:,location_r).',Ns); sac_interp_r_abs=abs(sac_interp_r); sac_interp_r_absmax=sac_interp_r_abs/max(sac_interp_r_abs); sac_interp_r_log=20*log10(sac_interp_r_absmax); rr=find(sac_interp_r_log==max(sac_interp_r_log)); sac_interp_a=interp(sac(location_a,:),Ns); sac_interp_a_abs=abs(sac_interp_a); sac_interp_a_absmax=sac_interp_a_abs/max(sac_interp_a_abs); sac_interp_a_log=20*log10(sac_interp_a_absmax); aa=find(sac_interp_a_log==max(sac_interp_a_log)); pslr_r1=pslrfunc(sac_interp_r_abs); %计算距离剖面图峰值旁瓣比dB islr_r1=islrfunc(sac_interp_r_abs); %计算距离剖面图积分旁瓣比dB pslr_a1=pslrfunc(sac_interp_a_abs); %计算方位剖面图峰值旁瓣比dB islr_a1=islrfunc(sac_interp_a_abs); %计算方位剖面图积分旁瓣比dB figure(7); subplot(121); plot(sac_interp_r_log); axis([rr-150,rr+150,-50,0]); ylabel('幅度dB'); title('(a)方位剖面图幅度'); subplot(122); plot(sac_interp_a_log); axis([aa-150,aa+150,-30,0]); ylabel('幅度dB'); title('(b)距离剖面图幅度'); point = Ptarget([1 11 17 20 21 ],:); %points to be illustrated a_locat = zeros(5, 1); r_locat = zeros(5, 1); %memorize the new coordinate(peaks) for k1=1:5 % locat = find(abs(Azimuth(:)-point(K1,1))-min(abs(Azimuth(:)-point(K1,1)))==0); %find the nearest azimuth point locat1= find(abs(Azimuth(:)-point(k1,1))-min(abs(Azimuth(:)-point(k1,1)))==0); a_locat(k1) = min(locat1); %record it locat2= find(abs(Range(:)-point(k1,2))-min(abs(Range(:)-point(k1,2)))==0); %find the nearest range point r_locat(k1) = locat2; %record it end oversam = 16; %oversampling ratio = 16 for k1 = 2:2 % t = sac(a_locat(k1),:); sar_interp_r = interp(sac(a_locat(k1),:),oversam); %interpolate 16 point between 2 in range sar_interp_log_r = 20*log10(abs(sar_interp_r)/abs(sar_interp_r(16*(r_locat(k1)-1)+1))); %normalization to the peak of each point figure;