FS-MD-PGA_PGA_SAR_FS成像_自聚焦

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%产生Stripmap SAR的回波 clear all clc %close all c=3e8;%光速 fc=94e9;%载频 lambda=c/fc;%波长 %%测绘带区域 X0=30;%方位向[-X0,X0] Rtc=0; Rrc=0; Ref=0; H=100; Rc=200; %R0=3;%距离向[Rc-R0,Rc+R0] %%距离向(Range),r/t domain Tr=1.6e-3*3;%LFM信号脉宽1.33us (200m) Br=1.4e9; %LFM信号带宽 150MHz Kr=Br/Tr; %调频斜率 fs=2*940e4;%快时间采样频率 aRou=c/2/Br;%距离分辨率 D=2*aRou; Lsar=lambda/D; Nr=1000; Na=5000; t=(0:Nr-1)/fs;%t域序列 % t=linspace(2*(Ref-R0)/c,2*(Ref+R0)/c+Tr,Nr); %f=linspace(Kr*2*(Rc-R0-Ref)/c,Kr*2*(Rc+R0-Ref)/c,Nr); f = linspace(-1/2*fs, 1/2*fs, length(t)); %%方位向(Azimuth,Cross-Range),x/u domain x=linspace(-X0,X0,Na);%u域序列 v=100; u=x/v;%慢时间 du=2*X0/v/Na; fu=linspace(-1/2/du,1/2/du,Na);%fu域序列 %dx=x(2) - x(1); %ku=2*pi*linspace(-1/dx/2, 1/dx/2, Na);%fu域序列 ku=fu*2*pi./v;%fu域序列 %％目标位置 Ptar=[500,0,1;%斜距,方位向坐标,sigma ]; Ntar=size(Ptar,1);%目标个数 %%产生回波 s_ut=zeros(Nr,Na); X=ones(Nr,1)*x;%扩充为矩阵 U=ones(Nr,1)*u; T=t'*ones(1,Na); Y=1*sin(70*U);%外加误差 Z=1*sin(70*U);%外加误差; for i=1:1:Ntar rn=Ptar(i,1);xn=Ptar(i,2);sigma=Ptar(i,3); rtn=rn+Rtc-Rrc; RT=sqrt((Y-sqrt(rtn^2-H.^2)).^2+(xn-X).^2+(Z+H).^2); R=2*RT; DT=T-R/c; phase=-j*4*pi/c*Kr*(T-2*Ref/c).*(RT-Ref) -j*4*pi/c*fc*(RT-Ref)+ j*4*pi*Kr/(c^2)*(RT-Ref).^2; s_ut=s_ut+sigma*exp(phase).*(abs(DT - Tr/2)<Tr/2).*(abs(X-xn)<(Lsar.*rn)/2); end %% KR=4*pi*Kr/c*(t-2*Ref/c); KRC=4*pi*fc/c; kr=KR+KRC; dkr = kr(2) - kr(1); r = linspace(-1/dkr/2, 1/dkr/2, length(kr))*2*pi; RB=r.'*ones(1,Na); b=8*pi*Kr/(c^2); KR=KR.'*ones(1,Na); Kx=(ones(Nr,1)*ku); Ax=sqrt(1-(Kx./(KRC)).^2); theta_r0=abs(sqrt((RB-Ref).^2-H^2)./(RB-Ref)); theta_r00=abs(H./(RB-Ref)); error_r0=-Y.*theta_r0+Z.*theta_r00; theta_rc=sqrt(Rc^2-H^2)/Rc; theta_rcc=H/Rc; error_rc=-Y*theta_rc+Z*theta_rcc; error_rv=error_r0-error_rc; C1=1; s_ut=s_ut.*C1;clear C1 s_ut=ifftshift(ifft(s_ut, [], 1), 1);%距离向IFFT C2=1; s_ut=s_ut.*C2;clear C2;clear error_r0;clear error_rc;clear error_rv; s_ut=fft(fftshift(s_ut,1),[], 1);%距离向FFT %*************FS************** HFS=exp(j*(KR).^2./(2*b).*(1-Ax)); HFS(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; s_if=fftshift(fft(s_ut, [], 2), 2);%方位向FFT S1=HFS.*s_if;clear HFS;clear s_if%去除弯曲空变性 S11=ifftshift(ifft(S1, [], 1), 1);clear S1;%距离向IFFT %*********RVP****************** H1=exp(-j*b/2*(1./Ax).*((RB-Ref).^2)); H1(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S111=S11.*H1;clear H1;clear S11;%剩余视频相位矫正 S2=fft(fftshift(S111, 1), [], 1);clear S111;%距离向FFT %*********IFS****************** H2=exp(1i*((Ax).^2 - Ax).*(KR).^2/2/b); H2(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S3=H2.*S2;clear H2;clear S2;%逆频率变标 error1=0; error2=0; error3=0; for i=20 Ax1=sqrt(1-((Kx+i*pi/0.11)./(KRC)).^2); error1=1;%(error1+exp(-j*(1-Ax)./2/b./Ax.*(Ax1.^2-Ax.^2).*(KR.^2))); error2=1;%(error2+exp(j*Rc/2/KRC.^3.*((Kx+i*2*pi*5/2).^2./Ax1-Kx.^2./Ax).*(KR+KRC).^1));%.*exp(-j*Rc/2/KRC.^4.*(-(Kx-i*2*pi*5/2).^2./Ax1.^2+Kx.^2./Ax.^2).*(KR+KRC).^3)); error3=1;%(error3+exp(j*1*(Ax1-Ax).*Ref.*KR)); end S3=S3.*error1./abs(error1);clear Ax1 %*************RCM*********************** HRMC=exp(-j*(Ax.*Ref-Ref).*(KR)).*error3; HRMC(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S4=S3.*HRMC;clear HRMC;clear S3;%距离徙动矫正 HSRC=exp(-j*(Rc.*Kx.^2)./(2*KRC^3.*Ax).*((KR).^2)).*exp(j*(Rc.*Kx.^2)./(2*KRC^4.*Ax.^2).*((KR).^3)); S5=S4.*HSRC;clear HSRC;clear S44;clear S4;%距离二次压缩 S55=ifftshift(ifft(S5, [], 1), 1);clear S5%距离向IFFT S55=ifft(ifftshift(S55, 2),[], 2);%方位向IFFT S55=S55;%.*C2; S55=fftshift(fft(S55, [], 2), 2);%方位向FFT HAREF=exp(j*Ax.*KRC.*RB); HAREF(Ax ~= conj(Ax), :) = 0; S6=S55.*HAREF;clear HAREF;%方位压缩 s=ifft(ifftshift(S6, 2),[], 2);clear S6;%方位向IFFT s=s.';%clear S; %% t1=clock; KR=4*pi*Kr/c*(t-2*Ref/c); KRC=4*pi*fc/c; kr=KR+KRC; dkr = kr(2) - kr(1); r = linspace(-1/dkr/2, 1/dkr/2, length(kr))*2*pi; RB=r.'*ones(1,Na); b=8*pi*Kr/(c^2); KR=KR.'*ones(1,Na); Kx=(ones(Nr,1)*ku); Ax=sqrt(1-(Kx./(KRC)).^2); theta_r0=abs(sqrt((RB-Ref).^2-H^2)./(RB-Ref)); theta_r00=abs(H./(RB-Ref)); error_r0=-Y.*theta_r0+Z.*theta_r00; theta_rc=sqrt(Rc^2-H^2)/Rc; theta_rcc=H/Rc; error_rc=-Y*theta_rc+Z*theta_rcc; error_rv=error_r0-error_rc; C1=exp(j*(KR+KRC).*error_rc); s_ut=s_ut.*C1;clear C1 s_ut=ifftshift(ifft(s_ut, [], 1), 1);%距离向IFFT C2=exp(j*KRC*error_rv); s_ut=s_ut.*C2;clear C2;clear error_r0;clear error_rc;clear error_rv; s_ut=fft(fftshift(s_ut,1),[], 1);%距离向FFT %*************FS************** HFS=exp(j*(KR).^2./(2*b).*(1-Ax)); HFS(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; s_if=fftshift(fft(s_ut, [], 2), 2);%方位向FFT S1=HFS.*s_if;clear HFS;clear s_if%去除弯曲空变性 S11=ifftshift(ifft(S1, [], 1), 1);clear S1;%距离向IFFT %*********RVP****************** H1=exp(-j*b/2*(1./Ax).*((RB-Ref).^2)); H1(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S111=S11.*H1;clear H1;clear S11;%剩余视频相位矫正 S2=fft(fftshift(S111, 1), [], 1);clear S111;%距离向FFT %*********IFS****************** H2=exp(1i*((Ax).^2 - Ax).*(KR).^2/2/b); H2(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S3=H2.*S2;clear H2;clear S2;%逆频率变标 error1=0; error2=0; error3=0; for i=20 Ax1=sqrt(1-((Kx+i*pi/0.11)./(KRC)).^2); error1=1;%(error1+exp(-j*(1-Ax)./2/b./Ax.*(Ax1.^2-Ax.^2).*(KR.^2))); error2=1;%(error2+exp(j*Rc/2/KRC.^3.*((Kx+i*2*pi*5/2).^2./Ax1-Kx.^2./Ax).*(KR+KRC).^1));%.*exp(-j*Rc/2/KRC.^4.*(-(Kx-i*2*pi*5/2).^2./Ax1.^2+Kx.^2./Ax.^2).*(KR+KRC).^3)); error3=1;%(error3+exp(j*1*(Ax1-Ax).*Ref.*KR)); end S3=S3.*error1./abs(error1);clear Ax1 %*************RCM*********************** HRMC=exp(-j*(Ax.*Ref-Ref).*(KR)).*error3; HRMC(:, Ax(1,:) ~= conj(Ax(1,:))) = 0; S4=S3.*HRMC;clear HRMC;clear S3;%距离徙动矫正 HSRC=exp(-j*(Rc.*Kx.^2)./(2*KRC^3.*Ax).*((KR).^2)).*exp(j*(Rc.*Kx.^2)./(2*KRC^4.*Ax.^2).*((KR).^3)); S5=S4.*HSRC;clear HSRC;clear S44;clear S4;%距离二次压缩 S55=ifftshift(ifft(S5, [], 1), 1);clear S5%距离向IFFT S55=ifft(ifftshift(S55, 2),[], 2);%方位向IFFT S55=S55;%.*C2; S55=fftshift(fft(S55, [], 2), 2);%方位向FFT HAREF=exp(j*Ax.*KRC.*RB); HAREF(Ax ~= conj(Ax), :) = 0; S6=S55.*HAREF;clear HAREF;%方位压缩 S=ifft(ifftshift(S6, 2),[], 2);clear S6;%方位向IFFT S=S.';%clear S; % figure % % s1=20*log10(abs(S)); % % s2=max(max(s1)); % % s3=20*(s1+20-s2).*(s1>(s2-20)); % % imagesc(r(660:750),x(1500:3500),20*log10(abs(S(1500:3500,660:750)))),axis xy % imagesc(linspace(r(660),r(750),5000),linspace(x(1500),x(3500),5000),interpft(interpft(20*log10(abs(S(1500:3500,660:750))),5000,1),5000,2)),axis xy % title('原始图像','fontsize',15), xlabel('Range[m]','fontsize',15), ylabel('Cross Range[m]','fontsize',15), % set(gca,'Clim',[-45,20]) % % % % % set(gcf,'unit','normalized','position',[0.2,0.2,0.5,0.32]); % figure,plot(x(2460:2540),20*log10(abs(S(2460:2540,708))),'LineWidth',2),hold on ,... % plot(x(2460:2540),20*log10(abs(I(2460:2540,708))),'r','LineWidth',2),... % hold on,plot(x(2460:2540),circshift(20*log10(abs(S1(2460:2540,708))),35),'r','LineWidth',2),xlabel('方位位置[m]','fontsize',15)... % ,ylabel('幅度[dB]','fontsize',15),legend('一阶补偿','二阶补偿','三阶补偿'),title('方位向剖物面对比图','fontsize',15) %% MDA S555=ifft(ifftshift(S55, 2),[], 2);%方位向IFFT %S555=S555.*exp(-j*sin(600*(U))); %S555=S55.*exp(j*pi*1.25e4.*U.^2).*(abs(X-xn)<(Lsar.*rn)/2); AA=708; NUM=length(AA); for ww=1:1:3 for i=1:199 ss_ut(:,1:50,i)=S555(:,1+25*(i-1):50+25*(i-1)); end Tkr=-2*v^2./lambda./r;%多普勒调频率参考值 Tkr=Tkr.'*ones(1,Na); for i=1:199 for uu=1:NUM ss_ut(:,1:25,i)=ss_ut(:,1:25,i);%.*exp(-j*pi*(9.98e3.*U(:,i*50-49:i*50-25).^2));%前一孔径 ss_ut(:,26:50,i)=ss_ut(:,26:50,i);%.*exp(-j*pi*(9.98e3.*U(:,i*50-25+1:i*50).^2));%后一孔径 f_ut=fftshift(fft(ss_ut(:,1:25,i), 5000, 2), 2);%前孔径方位向FFT b_ut=fftshift(fft(ss_ut(:,26:50,i),5000, 2), 2);%前孔径方位向FFT Cc=xcorr(abs(f_ut(AA(uu),:)),abs(b_ut(AA(uu),:))); [y_max