合成孔径雷达点目标仿真 SAR CS

%%================================= clear;clc;close all; %%================================= % 参数设置 C=3e8; Fc=1e9; lambda=C/Fc; Xmin=-50; %成像范围 Xmax=50; Yc=10000; Y0=500; V=100; %等效雷达速度 H=5000; R0=sqrt(Yc^2+H^2); %场景中心最短斜距 % 天线参数 D=4; Lsar=lambda*R0/D; Tsar=Lsar/V; %%================================= % 方位维参数 Ka=-2*V^2/lambda/R0; %方位向调频率 Ba=abs(Ka*Tsar); %方位向处理带宽 PRF=Ba; %方位向采样频率 PRT=1/PRF; ds=PRT; %时间间隔 Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %慢时间采样个数 Nslow=2^nextpow2(Nslow); %适应FFT后的慢时间采样个数 sn=linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow); %方位向距离向量 PRT=(Xmax-Xmin+Lsar)/V/Nslow; %修正后采样间隔 PRF=1/PRT; %修正后采样频率 ds=PRT; %%================================= % 距离维参数 Tr=5e-6; %脉冲时间 Br=30e6; %距离向脉冲带宽 Kr=Br/Tr; %距离向调频率 Fsr=3*Br; %方位向采样频率 dt=1/Fsr; %采样间隔 Rmin=sqrt((Yc-Y0)^2+H^2); %最短斜距 Rmax=sqrt((Yc+Y0)^2+H^2+(Lsar/2)^2); %最长斜距 Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt); %快时间采样个数 Nfast=2^nextpow2(Nfast); %FFT修正后的快时间采样个数 tm=linspace(2*Rmin/C-Tr/2,2*Rmax/C+Tr/2,Nfast); %时间向量 dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; %时间间隔 Fsr=1/dt; %修正后采样频率 RNear=Rmin-Tr/2*C/2; RFar=Rmax+Tr/2*C/2; %%================================= % 分辨率 DY=C/2/Br; DX=D/2; %%================================= % 目标参数 Ntarget=3; %format [x, y, reflectivity] Ptarget=[(Xmin+Xmax)/2,Yc,1 (Xmin+Xmax)/2,Yc-Y0/2,1 (Xmin+Xmax)/2,Yc+Y0/2,1]; %%================================= % 产生回波数据 T=Ptarget; Srnm=zeros(Nslow,Nfast); for k=1:1:Ntarget sigma=T(k,3); Dslow=sn*V-T(k,1); R=sqrt(Dslow.^2+T(k,2)^2+H^2); tau=2*R/C; Dfast=ones(Nslow,1)*tm-tau'*ones(1,Nfast); phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,Nfast)); Srnm=Srnm+sigma*exp(1i*phase).*(-Tr/2<Dfast&Dfast<Tr/2).*((abs(Dslow)<Lsar/2)'*ones(1,Nfast)); end figure; mesh(real(Srnm)); view(2); xlabel('距离向采样点'); ylabel('方位向采样点'); title('回波数据'); axis( [0 Nfast 0 Nslow] ) ; %%================================= % 傅立叶变换，转换到距离多普勒域 Srnm_rA=ftx(Srnm); %第一步骤 figure; mesh(real(Srnm_rA)); view(2); xlabel('距离向采样点'); ylabel('方位向采样点'); title('距离多普勒域数据'); axis( [0 Nfast 0 Nslow] ) ; DR = linspace(RNear,RFar,Nfast); dx = ds*V; fa = linspace(-1/2,1/2,Nslow)*PRF; %方位向频率 fr = linspace(-1/2,1/2,Nfast)*Fsr; %距离向频率 D_fa_V = sqrt(1-(C*fa).^2/(4*V^2*Fc^2)); %距离徙动因子D 《参考书》P187 7.17 Km = Kr./(1-(Kr*C*(ones(Nslow,1)*DR).*((fa.^2).'*ones(1,Nfast)))./(2*V^2*Fc^3*(D_fa_V.^3).'*ones(1,Nfast))); % 新的调频率km P187 7.18 tmN = ones(Nslow,1)*tm - (2*R0./(C*D_fa_V)).'*ones(1,Nfast);%相对参考时间 《参考书》P189 7.27 Kxref = 0; D_faref_V = sqrt(1-(C*Kxref)^2/(4*V^2*Fc^2)); %参考频率徙动因子 Srnm_rA_nocs=Srnm_rA; % 与变标方程进行相乘 Srnm_rA = Srnm_rA.*exp(1i*pi*Km.*(D_faref_V./(D_fa_V.'*ones(1,Nfast))-1).*tmN.^2); %第二步骤，第一次相位相乘 figure; % colormap('gray'); % imagesc(255-abs(Srnm_rA)); mesh(real(Srnm_rA)); view(2); xlabel('距离向采样点'); ylabel('方位向频率'); title('第一次相位相乘'); axis( [0 Nfast 0 Nslow] ) ; % 距离维匹配滤波 Srnm_RA=fty(Srnm_rA); %变换到二维频域，第三步骤 figure; % colormap('gray'); % imagesc(255-abs(Srnm_RA)); mesh(real(Srnm_RA)); view(2); xlabel('距离向频率采样点'); ylabel('方位向频率采样点'); title('二维频域距离压缩和一致RCM校正前数据'); axis( [0 Nfast 0 Nslow] ) ; Refr = exp(1i*pi*(ones(Nslow,1)*fr).^2.*(D_fa_V.'*ones(1,Nfast)./Km ) );%Refr是《参考书》P191 7.32 中二维频域信号表达式的第二个指数项 %是变标后的距离调制 Srnm_RA = Srnm_RA.*Refr; %相位相乘实现距离压缩和二次压缩，第四步 % 一致RCM的校正 RCMbulk = exp( 1i*4*pi/C*R0*ones(Nslow,1)*fr.*( (1./D_fa_V-1./D_faref_V).'*ones(1,Nfast)) );%精确方式的一致RCM《参考书》P191 7.32 第四个指数项 Srnm_RA = Srnm_RA.*RCMbulk; %第二次相位相乘实现一致RCM校正，第四步 figure; % colormap('gray'); % imagesc(255-abs(Srnm_RA)); mesh(real(Srnm_RA)); view(2); xlabel('距离向频率采样点'); ylabel('方位向频率采样点'); title('二维频域第二次相位相乘后数据'); axis( [0 Nfast 0 Nslow] ) ; Srnm_rA=ifty(Srnm_RA); %变换到距离多普勒域，第五步 figure; colormap('gray'); imagesc(255-abs(Srnm_rA)); % mesh(real(Srnm_rA)); % view(2); xlabel('距离向采样点'); ylabel('方位向频率采样点'); title('距离向处理后距离多普勒域数据'); axis( [0 Nfast 0 Nslow] ) ; figure; subplot(211); plot(abs(Srnm_rA(212,(300:700))));axis tight xlabel('距离向'),ylabel('方位向频率'), title('CS处理后距离多普勒域方位向频率第212个采样点切片'), subplot(212); plot(abs(Srnm_rA(70,(300:700))));axis tight xlabel('距离向'),ylabel('方位向频率'), title('CS处理后距离多普勒域方位向频率第70个采样点切片'), % % %无Chirp Scaling距离压缩%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Srnm_RA_nocs=fty(Srnm_rA_nocs); % Refr = exp(1i*pi*(ones(Nslow,1)*fr).^2.*(D_fa_V.'*ones(1,Nfast)./Km ) );%Refr是《参考书》P191 7.32 中二维频域信号表达式的第二个指数项 % %是变标后的距离调制 % Srnm_RA_nocs = Srnm_RA_nocs.*Refr; %相位相乘实现距离压缩和二次压缩，第四步 % % % 一致RCM的校正 % RCMbulk = exp( 1i*4*pi/C*R0*ones(Nslow,1)*fr.*( (1./D_fa_V-1./D_faref_V).'*ones(1,Nfast)) );%精确方式的一致RCM《参考书》P191 7.32 第四个指数项 % Srnm_RA_nocs = Srnm_RA_nocs.*RCMbulk; %第二次相位相乘实现一致RCM校正，第四步 % % Srnm_rA_nocs=ifty(Srnm_RA_nocs); % figure; % mesh(real(Srnm_rA_nocs)); % view(2); % xlabel('距离向频率采样点'); % ylabel('方位向频率采样点'); % title('距离向处理后距离多普勒域数据'); % axis( [0 Nfast 0 Nslow] ) ; % % figure; % subplot(211); % plot(abs(Srnm_rA_nocs(212,(300:700))));axis tight % xlabel('Range'),ylabel('Azimuth'), % title('SAR after range compression'), % subplot(212); % plot(abs(Srnm_rA_nocs(70,(300:700))));axis tight % xlabel('Range'),ylabel('Azimuth'), % title('SAR after range compression'), % Dphase1 = 4*pi*Fc*(ones(Nslow,1)*DR).*(D_fa_V.'*ones(1,Nfast))./C; %方位向匹配滤波补偿相位 % %附加相位补偿 % Dphase2 = -4*pi/C^2*Km.*(1-(D_fa_V./D_faref_V).'*ones(1,Nfast)).*( (ones(Nslow,1)*DR)./(D_fa_V.'*ones(1,Nfast)) - R0./(D_fa_V.'*ones(1,Nfast))).^2; % Srnm_rA_nocs = Srnm_rA_nocs.*exp(1i*Dphase1).*exp(1i*Dphase2);%第三