欺骗干扰SAR成像——MATLAB实现.rar

clear all clc close all %% %雷达相关参数 j=sqrt(-1); c=3e8; %光速 fc=10e9; %载频 lambda=c/fc; %波长 PRF=500; %脉冲重复频率 Tp=3e-6; %LFM信号脉宽 Br=50e6; %LFM信号带宽 Fs=2*Br; %距离向采样频率 Kr=Br/Tp; %调频斜率 betal=70*pi/180; %波束中心下视角 %斜视角 theta=0*pi/180; H0=20000; %雷达高度 v_x=1000; %导弹沿x轴的速度 匀速直线飞行 D=10; %天线真实孔径 rou_a=D/2; %方位分辨率 rou_r=c/2/Br/(cos(theta)*sqrt(tan(theta).^2+sin(betal).^2)); %距离分辨率 %% %目标场景 %X0=0; %场景中心的方位向坐标 X0=H0/cos(betal)*tan(theta) ; %场景中心的方位向坐标 20110803修改 Y0=H0*tan(betal); %场景中心的距离向坐标 R0=H0/cos(betal)/cos(theta); %场景中心线到雷达的距离 azirange=100; %方位向成像范围 ranrange=100; %距离向成像范围 Soder_X=50; %方位向点目标间距 Soder_Y=50; %距离向点目标间距 x_range=X0+(-azirange/2:Soder_X:azirange/2); %场景方位向范围 y_range=Y0+(-ranrange/2:Soder_Y:ranrange/2); %场景距离向范围 Ntar=1; %目标个数 target=[0,0,1 %目标相对于场景的坐标(方位坐标 距离坐标 灰度值) 0,1,1 -1,0,1 1,0,1 5,5,1 0,-1,1 ]; target(:,1)=target(:,1)*Soder_X;target(:,2)=target(:,2)*Soder_Y+Y0; %目标在绝对坐标系中的坐标 order=target; %% %方位向采样， fdr=-2*v_x^2*fc/c/R0; Ba=v_x/rou_a ; %满足要求分辨率所需要的多普勒带宽 TT=Ba/abs(fdr); %成像积累时间 Lsar=v_x*TT %合成孔径长度 Tc1=TT+azirange/v_x; %一次成像时间;包括成像区域 Na1=ceil(Tc1*PRF); %一次成像时间内需要的采样点 Na=2^nextpow2(Na1); %为了后面FFT，将Na设置为2的整数次方个 PRF=Na/Tc1; %刷新脉冲重复周期 fangwei=linspace(-azirange/2-Lsar/2-X0,azirange/2+Lsar/2-X0,Na); %20110803 修改 Tm=linspace(-Tc1/2,Tc1/2,Na)-X0/v_x; %方位向慢时间 %20110803 修改 x_rader=Tm*v_x; %雷达方位向历程 fu=linspace(-PRF/2,PRF/2,Na); %多普勒域 %% %距离向采样 %Rmin=R0-(ranrange/2)*sin(beta) %最小照射距离 Rmin=R0-(ranrange/2)*sin(betal) %20110803改 %Rmax=sqrt((R0+(ranrange/2)*sin(beta))^2+(Lsar/2)^2) %最大照射距离 Rmax=sqrt((R0*sec(theta)+(ranrange/2)*sin(betal))^2+(azirange+Lsar/2)^2) %20110803改 Nr=ceil((2*(Rmax-Rmin)/c+Tp)*Fs); %距离向采样点数 Nr=2^nextpow2(Nr); %为FFT，2的整数次方个 % Rrmx=zeros(1,Na); % Rrmn=zeros(1,Na); % for i=1:Na % Rrmx(i)=max(sqrt((xn-x_rader(i)).^2+(Y0+Y0ch-y_rader(i)).^2+(z_rader(i))^2)); % Rrmn(i)=min(sqrt((xn-x_rader(i)).^2+(Y0-Y0ch-y_rader(i)).^2+(z_rader(i))^2)); % end % Rmax=max(Rrmx) %最大照射距离 % Rmin=min(Rrmn) %最小照射距离 % Nr=ceil((2*(Rmax-Rmin)/c+Tp)*Fs); %距离向采样点数 dt=(2*Rmax/c-2*Rmin/c+Tp)/Nr; %重新采样 Fs=1/dt; %刷新采样频率 fr=linspace(-Fs/2,Fs/2,Nr); %距离频域 t=(0:Nr-1)*1/Fs; %距离时域 juli=linspace(-ranrange/2,ranrange/2,Nr); %% %%下面是回波模拟 %stt=rectpuls(t-Tp/2,Tp).*exp(j*2*pi*(0.5*Kr*(t-Tp/2).^2)); %%定义发射波 t=2*Rmin/c+t; %wa=zeros(1,Na); sr=zeros(Na,Nr); for k=1:Na sr(k,:)=zeros(1,Nr); for i=1:Ntar if abs(v_x*Tm(k)+X0-order(i,1))<Lsar/2 r(k)=sqrt((order(i,1)-x_rader(k))^2+(order(i,2))^2+H0^2); tr=2*r(k)/c; %wa(k)=rectpuls(Tm(i),TT); %-order(i,2)*Soder*tan(pha)??? wr=rectpuls(t-Tp/2-tr,Tp); s=order(i,3)*wr.*exp(j*2*pi*(0.5*Kr*(t-Tp/2-tr).^2-fc*tr)); sr(k,:)=sr(k,:)+s; end end end s=sr.'; %Jr=qipian(Na,Nr,Lsar,Tm,order,x_rader,t,Ntar,v_x,X0,H0,c,lambda,Tp,Kr,fc); %s=s+Jr; SAREcho=strcat('huiboshuju.mat'); save(SAREcho,'s'); %sat=rectpuls(Tm,Tc1).*exp(j*pi*4*v_x^2*Tm^2/lambda/R0); % figure(1) % plot_img(fangwei,juli,s,30),grid on % title('原始回波数据'); % xlabel('方位维范围');ylabel('距离维范围'); %% %CS成像 pi2=2*pi; wc =2*pi*fc; DeltaR = c/(2*Fs); %Xc = Xmin+DeltaR*ran_num/2; % Range distance to center of target area X0 = DeltaR*Nr/2; % Target area in range is within [Xc-X0,Xc+X0] Y0 = Na*v_x/PRF/2; % Target area in cross-range is within [Yc-Y0,Yc+Y0] % Rmin=Xc-0.5*ran_num*DeltaR; % Rmax=Xc+0.5*ran_num*DeltaR; Xc = R0; Rc=Xc; Ts=(2/c)*Rmin; % start time of sampling % Tf=(2/C)*Rmax; Rref=Xc; tm = ((0:Na-1) -Na/2)/PRF; t=Ts+(0:Nr-1)*dt; %回波时间加上集体的延迟时间 x = Rmin+(0:Nr-1)*DeltaR; %x=linspace(Rmin,Rmax,Nr) %%cross-range:y domain %y = linspace(-ranrange/2,ranrange/2,Nr); %y域序列:-Y0~Y0 y=linspace(-azirange/2-Lsar/2,azirange/2+Lsar/2,Na); %20110803修改 Deltay = v_x/PRF; %kyc = 0; %20110803 kyc = 4*pi/lambda*sin(theta); dky = pi2/(Na*Deltay); ky = dky*(-Na/2:Na/2-1)+kyc; %%range:x domain dw = pi2/(Nr*dt); % Frequency domain sampling w = dw*(-Nr/2:Nr/2-1); %(魏改：去掉载频wc) fr = w./pi2; %fac = 0; %20110803 fac = 2*v_x*sin(theta)/lambda; fa = (-Na/2:Na/2-1)*(PRF/Na)+fac; kxc = 4*pi/lambda; phi0 = sqrt(kxc^2-ky.^2); %Cs = ones(ran_num,1)*(kxc./phi0-1); Cs = (kxc./phi0-1);%CS因子??? phi1 = phi0./kxc;%cos(theta) phi2 = (ky./kxc).^2;%sin(theta)^2 temp1 = 2*lambda*Xc./c^2; %gamae = ones(ran_num,1)*(1./(1/gama-temp1.*phi2./(phi1.^3))); gamae = (1./(1/Kr-temp1.*phi2./(phi1.^3)));%等效调频斜率 %fac = 0; %20110803 fac = 2*v_x*sin(theta)/lambda; s=s.*exp(-j*2*pi*fac*ones(Nr,1)*tm); %%CSA:7 steps %s=fftshift(fft(fftshift(s).')).'; s = fty(s);%方位向FFT t = t.'; %Chirp Scaling 线频调变标处理 for idx = 1:Nr H1 = exp(j*pi*Cs.*gamae.*(t(idx)-2*(Xc*(1+Cs))./c).^2); s(idx,:) = s(idx,:).*H1; end clear H1; s = ftx(s);%距离向FFT for idx = 1:Na H2 = exp(j*pi*(fr.^2')./(1+Cs(idx))./gamae(idx)... +j*4*pi/c*Rc*Cs(idx).*(fr')); %距离向匹配滤波&距离迁移校正 s(:,idx) = s(:,idx).*H2; end clear H2; %距离向IFFT s = iftx(s); for idx = 1:Nr H3 = exp(j*4*pi/c^2*gamae.*Cs.*(1+Cs).*(((x(idx)-Xc).^2))... +j*4*pi/lambda*x(idx)*phi1);%消除误差函数，方位向匹配滤波 s(idx,:) = s(idx,:).*H3; end clear H3; clear gamae Cs; %方位向IFFT s = ifty(s); figure(3) imagesc(x,y,abs(s)) colormap(gray) Rb = x; figure; colormap(gray(256)); G=abs(s).'; xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng); %%%%%%%%%%%%%%%%%%%%%% 注意下面成像用的坐标是y，不是Yc+y Brightness = 2; axis xy; %yidongliang=round((x-Rc)*0*PRF/v_x);%%%%用这一句比下一句校正的好 yidongliang=round((x-Rc)*sin(theta)*PRF/v_x);%20110803修改 % yidongliang=round((X-Rc)*tan(theta_c)*PRF/v) ; for i=1:Nr G(:,i)=circshift(G(:,i),yidongliang(i)); end % image(x,f