RD.zip_RD_RD成像算法_多普勒rd算法_距离多普勒算法_雷达

%% %%使用邻近插值和sinc插值的距离迁徙校正仿真 %%========================================================================= %% Filename: %% Help file: .doc %% Project: %% Author: %% E-mail: %%========================================================================= %%========================================================================= % clear all; clc; close all; %%========================================================================= %% Parameter--constant C = 3e8; % speed of light Fc = 1e9; % carrier frequency lambda = C/Fc; % wavelength Rfactor = 2; % 距离维采样倍数因子 Afactor = 1; % 方位向采样倍数因子 %% Parameter--resolution 距离维和方位向的分辨率要求是多少 DY = 3; % range resolution DY = C/2/Br DX = 1; % azimuth resolution DX = D/2 = Va / Ba 分辨率都是两米 %% Parameter--target area 观察场景相关参数 Xmin = 0; % Targets area in azimuth is within [Xmin Xmax] Xmax = 200; Yc = 10000; % Center of imaged area 场景中心线 Yw = 300; % Target area in range is within[Yc-Yw,Yc+Yw] imaged width is 2*Yw %% Parameter--orbital information 载机平台相关参数 V = 100; % SAR velocity 100m/s 载机速度 H = 6000; % height of airborne 6000m 载机运行高度 R0 = sqrt(H^2+Yc^2); % 载机的中心垂直距离 是常数 Yc是场景中心线距离 theta=90*pi/180; % squint angle 斜视角 %% Parameter--antenna 天线参数 D = 2*DX; % antenna length in azimuth direction 由方位向分辨率决定 Lsar = lambda*R0/D; % SAR integration length 合成孔径长度 Tsar = Lsar/V; % SAR integration time 合成孔径时间 %% Parameter--fast-time domain 快时间域参数设置 Tr = 10e-6; % pulse duration 10us 发射脉冲宽度 Br = C/(2*DY); % chirp frequency modulation bandwidth 30MHz LFM信号的带宽 由距离维分辨率决定 Kr = -Br/Tr; % chirp slope LFM信号的调频率 这里为了看回波采用负调频 Fsr = Rfactor*Br; % sample frequency in fast-time domain 距离维采样频率 由奈奎斯特原理知要大于等于带宽 dt = 1/Fsr; % sample spacing in fast-time domain 距离维的采样间隔 Rmin = sqrt((Yc-Yw)^2+H^2); % 最小斜距在距离维 Rmax = sqrt((Yc+Yw)^2+H^2); % 最大斜距在距离维 Rwid = Rmax-Rmin; % 也就是波门的范围 Twid = 2*Rwid/C; Nfast = ceil((Twid +Tr)/dt); % sample number in fast-time domain 雷达成像技术书P23页说明 Nfast = 2^nextpow2(Nfast); % Ready for fft t_hat = linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);% discrete time array in fast-time domain 离散化 dt = (2*Rmax/C+Tr-2*Rmin/C)/Nfast; % 修正和更新采样间隔 Fsr=1/dt; %% Parameter--slow-time domain 慢时间域参数设置 Ka = -2*V^2/lambda/R0; % doppler frequency modulation rate 多普勒频率调制率 Ba = abs(Tsar*Ka); % doppler frequency modulation bandwidth 多普勒频率调制带宽 PRF = Afactor*Ba; % pulse repetition frequency 以Afactor倍的带宽去采样 也就是采样频率 PRT = 1/PRF; % pulse repetition interval 也就是采样间隔 ds = PRT; % sample spacing in slow-time domain 采样间隔 %因为我们知道慢时间域采样频率为PRF Rawid = Xmax-Xmin+Lsar; % 慢时间域的距离窗 要根据成像区域来确定慢时间的范围 Tawid = Rawid/V; % 慢时间域的时间窗 因为慢时间域距离是单程的 Nslow = ceil(Tawid/ds); % sample number in slow-time domain Nslow = 2^nextpow2(Nslow); % ready for fft; eta = linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow); % discrete time array in slow-time domain 根据成像区域确定的 PRT = Tawid/Nslow; % 修正和更新采样间隔 也就是上面做了近似 PRF = 1/PRT; ds = PRT; %% Parameter--point targets 点目标 %format [x, y, reflectivity] Ptarget = [Xmin+20,Yc-300,1 Xmin+30,Yc+300,1 Xmin+40,Yc+150,1 Xmin+50,Yc+50,1 Xmin+60,Yc+10,1 Xmin+70,Yc+20,1 Xmin+80,Yc+30,1 Xmin+90,Yc+40,1 Xmin+100,Yc+60,1 Xmin+110,Yc+80,1 Xmin+120,Yc+70,1 Xmin+130,Yc+90,1 Xmin+140,Yc+110,1 Xmin+150,Yc+120,1 Xmin+160,Yc+100,1 Xmin+170,Yc+140,1]; %一个点目标 Ntarget = size(Ptarget,1); % number of targets %% ----Display some SAR Parameter---- disp('Parameters:') disp('Sampling Rate in fast-time domain');disp(Fsr/Br) disp('Sampling Number in fast-time domain');disp(Nfast) disp('Sampling Rate in slow-time domain');disp(PRF/Ba) disp('Sampling Number in slow-time domain');disp(Nslow) disp('Range Resolution');disp(DY) disp('Cross-range Resolution');disp(DX) disp('SAR integration length');disp(Lsar) disp('Position of targets');disp(Ptarget) %% Generate the raw signal data 产生需要处理的数据矩阵 K = Ntarget; % number of targets N = Nslow; % number of vector in slow-time domain M = Nfast; % number of vector in fast-time domain T = Ptarget; % position of targets Srnm = zeros(N,M); % 初始化回波数据矩阵 h = waitbar(0,'Generating echo signal'); % 生成进度条 for k = 1:1:K sigma = T(k,3); Dslow = eta*V - T(k,1); % 计算 V*eta - x Ra = sqrt(H^2 + Dslow.^2 + T(k,2)^2); % 斜距计算公式 tau = 2*Ra/C; % 计算得到每个目标的延迟 Dfast=ones(N,1)*t_hat-tau'*ones(1,M); % 计算得到快时间变量 tm是快时间离散变量 sn是慢时间离散变量 phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(Ra'*ones(1,M)); Srnm=Srnm + sigma * exp(1i*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M)); %累加 waitbar(k/K); %更新进度条状态 end close(h) %关闭进度条指示栏 %% 距离向匹配滤波 range compression tr=t_hat -2*Rmin/C; %为什么那？ % tr1 = linspace(0,Tr,ceil(Tr/dt)); % tr=zeros(1,Nfast); % tr(1:length(tr1))=tr1; %加零补长到tm的长度 也算是一种补零加长数据的一种方式 Refr=exp(1i*pi*Kr*tr.^2).*(0<tr&tr<Tr); Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); Gr = abs(Sr); figure waterfall(real(Srnm(500,:))); figure waterfall(Gr(500,:)); %% 开始进行距离弯曲补偿正侧视没有距离走动项 主要是因为斜距的变化引起回波包络的徙动 %%补偿方法：1:时域插值法 2：频域相移法 %%方位向做fft处理 再在频域做距离弯曲补偿 Sa_RD = ftx(Sr); % 方位向FFT 变为距离多普域进行距离弯曲校正 %%第一种插值方法--最近邻插值 %距离徙动运算,由于是正侧视 ，fdc=0,只需要进行距离弯曲补偿。 % Kp=1; %计算或者预设预滤波比 % h = waitbar(0,'最近邻域插值'); % %%首先计算距离迁移量矩阵 % for n=1:N % for m=1:M % delta_R = (1/8)*(lambda/V)^2*(R0+(m-M/2)*C/2/Fsr)*((n-N/2)*PRF/N)^2;%距离迁移量 % RMC=2*delta_R*Fsr/C; % delta_RMC = RMC-round(RMC); % if m+round(RMC)>M %判断是否超出边界 % Sa_RD(n,m)=Sa_RD(n,M/2); % else % if delta_RMC>=0.5 % Sa_RD(n,m)=Sa_RD(n,m+round(RMC)+1); % else % Sa_RD(n,m)=Sa_RD(n,m+round(RMC)); % end % end % end % waitbar(n/N) % end % close(h) % %======================== % Sr_rmc=iftx(Sa_RD); %%距离徙动校正后还原到时域 % Ga = abs(Sr_rmc); %%第二种方法进行截断sinc插值进行距离徒动校正 h = waitbar(0,'Sinc插值'); P=4;%4点sinc插值 RMCmaxtix = zeros(N,M); for n=1:N % N - # of azimuth samples for m=P:M delta_R = (1/8)*(lambda/V)^2*(R0+(m-M/2)*C/2/Fsr)*((n-N/2)*PRF/N)^2;%首先计算距离迁移量 计算方法就是把斜距变换到距离多普勒域就知道了 RMC=2*delta_R*Fsr/C; %距离徒动了几个距离单元 delta_RMC = RMC-round(RMC); %距离徒动量的小数部分 for i = -P/2:P/2-1