SAR成像程序，包括补偿。

%beginning clear; clc; close all; %% Parameter--constant C=3e8; % speed of light光速 Fc = 1e9; % carrier frequency载频 lambda = C/Fc; % wavelength波长 Rfactor = 2; % 距离维采样倍数因子 Afactor = 1; % 方位维采样倍数因子 %% Parameter--target area Xmin = 0; % Targets area in azimuth is within [Xmin Xmax] Xmax = 200; Yc = 10000; %Center of imaged area 场景中心线 即RB：点目标到航线的最短距离 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); % 载机的中心垂直距离H是常数 Yc是场景中心线距离 theta=90*pi/180; % squint angle 斜视角rad %% Parameter--antenna D=2; % antenna length in azimuth direction px = D/2 Lsar = lambda*R0/D; % SAR integration length 合成孔径长度 Tsar = Lsar/V; % SAR integration time 合成孔径时间 %% Parameter--fast time domain Tr = 10e-6; %pulse duration 10us 发射脉冲宽度 Br = 3e7; %chirp frequency modulation bandwidth 30MHz LFM信号的带宽 由距离维分辨率决定 Kr = -Br/Tr; %chirp slope LFM信号的调频率 这里为了看回波采用负调频 Fsr = Rfactor*Br; %sample frequency in fast-time domain 距离维采样频率 由奈奎斯特原理知要大于等于带宽 Tsr = 1/Fsr; %sample spacing in fast-time domain 距离维的采样间隔 Rmin = sqrt((Yc-Yw)^2+H^2); %最小斜距在距离维 %Rmax = sqrt((Yc+Yw)^2+H^2+(Lsar/2)^2); %最大斜距在距离维 结果差不多 Rmax = sqrt((Yc+Yw)^2+H^2); %最大斜距在距离维 Rwid = Rmax-Rmin; %也就是波门的范围 Twid = 2*Rwid/C; %相干积累时间？ Nfast = ceil(Twid/Tsr+Tr/Tsr); %sample number in fast-time domain 雷达成像技术书P23页说明 采样个数 ceil（1023.9）=1024 Nfast = 2^nextpow2(Nfast); %Ready for fft 此操作后，必须更新采样频率和采样间隔 tn=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast); %discrete time array in fast-time domain 离散化 Tsr=(2*Rmax/C+Tr-2*Rmin/C)/Nfast; % 更新采样间隔 Fsr=1/Tsr; %更新采样频率 %% Parameter--resolution DY =C/(2*Br); %range resolution DX = D/2; %cross-range resolution %% Parameter--slow-time domain Ka = -2*V^2/lambda/R0; % doppler frequency modulation rate 多普勒频率调制率 公式5.19 Ba = abs(Tsar*Ka); % doppler frequency modulation bandwidth 多普勒频率调制带宽 PRF = Afactor*Ba; % pulse repetition frequency 以Afactor倍的带宽去采样 也就是采样频率???? PRT = 1/PRF; % pulse repetition interval 也就是采样间隔 Tsa = PRT; % sample spacing in slow-time domain 采样间隔 % 因为我们知道慢时间域采样频率为PRF (因为。。) ？ Rawid = Xmax-Xmin+Lsar; %慢时间域的距离窗 要根据成像区域来确定慢时间的范围 Tawid = Rawid/V; %慢时间域的时间窗 因为慢时间域距离是 单程的 Nslow = ceil(Tawid/Tsa); %sample number in slow-time domain 方位采样数 Nslow = 2^nextpow2(Nslow); %ready for fft， 此操作后，必须更新采样频率和采样间隔 2048 tm = linspace((Xmin-Lsar/2)/V,(Xmax+Lsar/2)/V,Nslow); %discrete time array in slow-time domain 根据成像区域确定的 Tsa= Tawid/Nslow; % 更新方位采样间隔 也就是上面做了近似 Fsa = 1/Tsa; % 更新方位采样频率 PRF = Fsa; PRT = Tsa; %% Parameter--point targets 点目标 Ntarget=3; %number of targets %format [x, y, reflectivity] Ptarget = [Xmin+20,Yc, 1 Xmin+20,Yc+100,1 Xmin+40,Yc+200,1 ]; % corordinates(position)坐标 %% Display some SAR Parameters 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,'SAR回波生成'); % 生成进度条，h为句柄 for k = 1:1:K sigma = T(k,3); Dslow = tm*V-T(k,1); %计算 V*tm-x 方位向长度 tm 1xN =2048 Dslow 1xN=2048 Ra = sqrt(H^2+Dslow.^2+T(k,2)^2); %斜距计算公式 R（tm;RB） tau=2*Ra/C; %计算得到每个目标的延迟 Dfast=ones(N,1)*tn-tau'*ones(1,M); %计算得到快时间变量 tn是快时间离散变量 tm是慢时间离散变量（tn-2*Ra/C） phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(Ra'*ones(1,M)); %pi*(tn-2*Ra/C)^2-4*pi*R(tm;RB)/lambda Srnm=Srnm+sigma*exp(1i*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M)); %累加 waitbar(k/K); %更新进度条状态,括号里k/K必须是0-1的数 end close(h) %关闭进度条指示栏 %% 距离向匹配滤波 range compression %tr=tn-2*Rmin/C; %为什么？ tr1 = linspace(0,Tr,ceil(Tr/Tsr)); %[0 Tr]脉宽里取Tr/Tsr个点 1x763个 tr=zeros(1,Nfast); %1xN=2048 tr(1:length(tr1))=tr1; %加零补长到tn的长度 也算是一种补零加长数据的一种方式？ Refr=exp(1i*pi*Kr*tr.^2).*(0<tr&tr<Tr); %距离匹配滤波函数 Sr=ifty(fty(Srnm).*(ones(N,1)*conj(fty(Refr)))); %匹配滤波（脉冲压缩）行fft Gr = abs(Sr); %%开始进行距离弯曲补偿，正侧视没有距离走动项 主要是因为斜距的变化引起回波包络的徙动 %%补偿方法：1:时域插值法 2：频域相移法 ：方位向做fft处理 再在频域做距离弯曲补偿 Sa_RD = ftx(Sr); % 方位向FFT 变为距离-多普域（tn-fa域）进行距离弯曲校正 %% 第一种插值方法-- 最近邻插值 % 距离徙动运算,由于是正侧视 ，fdc=2vsinβ/λ=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 %小数部分大于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); %取模 tm=tm-Xmin/V; % 减去 近距时间 Refa=exp(1i*pi*Ka*tm.^2).*(abs(tm)<Tsar/2); %方位匹配滤波参考函数 Tsar/2是为了严谨么？还是固定的，Tsar可以么？ Sa=iftx(Sa_RD.*(conj(ftx(Refa)).'*ones(1,M))); %方位脉压 距离方位域 输出信号 Ga1=abs(Sa); %% %第二种方法进行 截断sinc插值 进行距离徒动校正 % h = waitbar(0,'Sinc插值'); % P=4; %4点sinc插值 % RMCmaxtix = zeros(N,M); %RMCmaxtix这个