rd.zip_matlab 图像处理_rd算法_合成孔径雷达

%%%%RD算法 两种方法距离徙动矫正%%%%% clear all;close all;clc; c=3e8; j=sqrt(-1); R0=30e3; %景中心斜距 Vr=150; %速度 Tp=2.5e-6; %脉宽 f0=5.3e9; %载频 Kr=20e12; %距离向调频率 Fr=150e6; %距离向采样率 %Ba=80; %多普勒带宽 theta_R=0; %斜视角 La=2; %天线真实孔径 Ba=0.886*2*Vr*cos(theta_R)/La; %多普勒带宽 %PRF=100; PRF=5*Ba; %PRF lamda=c/f0; %波长 Ls=0.886*lamda*R0/La; %合成孔径长度 x_a=30; %方位向范围,范围估算 x_r=500; %距离向范围 Na=2^nextpow2((2*x_a+Ls)/Vr*PRF); %方位向采样点数 Nr=ceil((2*(x_r)/c+Tp)*Fr); %斜距最远和最近的双程差所用时间再加上脉冲持续时间 Nr=2^nextpow2(Nr); % Na=256; %方位向采样点数 % Nr=320; %距离向采样点数 target=[0,0;0,50;50,50]; %点目标相对中心位置[距离向,方位向] data=zeros(Na,Nr); ta=[-Na/2 :Na/2-1]/PRF; tr=[-Nr/2 :Nr/2-1]/Fr; tr=tr+2*R0/c; %距离向相对时间 for k=1:1:3 Ra=ta*Vr-target(k,2); %方位向距离差 Rmin=R0+target(k,1); %最近斜距 Rn=sqrt(Ra.^2+Rmin^2); %瞬时斜距 tao=2*Rn/c; %时间延迟 t=ones(Na,1)*tr-(tao.')*ones(1,Nr); %方位相对时间 data=data+exp(-j*4*pi*f0*Rn.'/c)*ones(1,Nr).*exp(j*pi*Kr*t.^2).*(abs(t)<Tp/2).*((abs(Vr*t)<Ls/2)); end figure(1); imagesc(abs(data)); title('存储数据') %% 距离向脉压 f_r=(-Nr/2:Nr/2-1)/Nr*Fr; %距离向频率，中心0频 H_r=exp(j*pi/Kr*f_r.^2); for i=1:1:Na data(i,:)=ifftshift(ifft(ifftshift(fftshift(fft(fftshift(data(i,:)))).*H_r))); % data(i,:)=ifft(fftshift(( fftshift(fft(data(i,:))))).*H_r); end figure(2); imagesc(abs(data)); title('距离向脉压'); %% 方位向FFT、RCMC、方位向脉压 Ka=2*Vr^2/lamda/R0; %方位向调频率 % f_a=[1:Na/2,-Na/2:-1]/Na*PRF; %方位向频率 f_a=(-Na/2:Na/2-1)/Na*PRF; H_a=exp(-j*pi/Ka.*f_a.^2); data0=fftshift(fft(fftshift(data))); %用作未矫正的数据 data1=fftshift(fft(fftshift(data))); %方位向fft data=fftshift(fft(fftshift(data))); %方位向fft figure(3); imagesc(abs(data)); title('校正前距离多普勒域图像'); %% GRCM % D=sqrt(1-((lamda^2)*(f_a.^2)/(4*Vr*Vr))); % deltaR=R0*((1-D)./D); %SRC情况下的距离徙动量 % deltaR=deltaR.'*ones(1,Nr); % rf=ones(Na,1)*f_r; % Grcmc=exp(j*4*pi*rf.*deltaR/c); %相位乘法器 % pecho_2f2=Grcmc.*fftshift(fft(fftshift(data.'))).'; % data=fftshift(ifft(fftshift(pecho_2f2.'))).'; % figure(4); % imagesc(abs(data)); % title('GRCM校正后距离多普勒域图像'); %% sinc差值RCMC bar=waitbar(0,'sinc插值'); for k=1:Na %按行扫描数据 r_line=data1(k,:); %取出第K行数据 for n=400:950 D=sqrt(1-lamda^2*f_a(k)^2/4/Vr^2); rma=2*R0*(1/D-1)/c*Fr; %结合上式得出需要校正的距离徙动的采样点数 %data(k,n)=r_line(n+round(rma)); %以最近点代替插值点 h=floor(rma); g=ceil(rma); x=[r_line(n+h-1),r_line(n+h),r_line(n+g),r_line(n+g+1)]; y=[sinc(rma-h+1),sinc(rma-h),sinc(g-rma),sinc(g-rma+1)];% 插值 data1(k,n)=x*y'; end waitbar(k/Na) end close(bar); figure(5); imagesc(abs(data1)); title('sinc插值校正后距离多普勒域图像'); %% 方位向脉压 for n=1:Nr data(:,n)=data(:,n).*H_a.'; end data=fftshift(ifft(fftshift(data))); %RCMC成像结果 for n=1:Nr data0(:,n)=data0(:,n).*H_a.'; end data0=fftshift(ifft(fftshift(data0))); %未经RCMC成像结果 for n=1:Nr data1(:,n)=data1(:,n).*H_a.'; end data1=fftshift(ifft(fftshift(data1))); figure(6); subplot(1,2,1); imagesc(abs(data0)); title('未经RCMC成像结果'); %axis([460,560 1500 2500]); % subplot(1,3,2); % imagesc(abs(data)); % title('GRCM后成像结果'); % %axis([460,560 1500 2500]); subplot(1,2,2); imagesc(abs(data1)); title('sinc插值矫正后成像结果'); %axis([460,560 1500 2500]); %% 质量评估 %选取目标周围33*33大小的区域并进行10倍升采样 sample=data1(Na/2-16:Na/2+16,Nr/2-16:Nr/2+16); sample_freq=fft2(sample); figure(9); imagesc(abs(sample_freq)); sample=zeros(16*10,16*10); sample(1:16,1:16)=sample_freq(1:16,1:16); sample(1:16,9*16:10*16)=sample_freq(1:16,17:33); sample(9*16:10*16,9*16:10*16)=sample_freq(17:33,17:33); sample(9*16:10*16,1:16)=sample_freq(17:33,1:16); sample=ifft2(sample); figure(7); imagesc(abs(sample)); title('目标升采样结果'); figure(71); surf(abs(sample)); %距离向剖面 section_r=sample(80,:); section_r=abs(section_r)/max(abs(section_r)); section_r=20*log10(section_r); %方位向剖面 section_a=sample(:,80); section_a=abs(section_a)/max(abs(section_a)); section_a=20*log10(section_a); figure(8); subplot(2,1,1); plot(section_r); title('距离向剖面图'); subplot(2,1,2); plot(section_a); title('方位向剖面图'); %峰值旁瓣比PSLR peaks_r=findpeaks(section_r); peaks_r=sort(peaks_r,'descend'); PSLR=peaks_r(2)-peaks_r(1) % peaks_a=findpeaks(section_r); % peaks_a=sort(peaks_a,'descend'); % PSLR_a=peaks_a(2)-peaks_a(1) %积分旁瓣比ISLR section_r=abs(sample(80,:)); pos=[]; k=0; for i=2:1:159 if section_r(i)<section_r(i+1)&&section_r(i)<section_r(i-1) k=k+1; pos(k)=i; else k=k; end end Pmain=0; for i=pos(k/2):1:pos(k/2+1) Pmain=Pmain+section_r(i).^2; end Ptotal=0; for i=1:1:160 Ptotal=Ptotal+section_r(i).^2; end ISLR=10*log10((Ptotal-Pmain)/Pmain) %3dB主瓣宽度IRW IRW=0.886/(abs(Kr)*Tp)