BiSARBP_sarbp_SARBP成像_SAR_BP成像_双基地SAR成像

clc; close all; % 恒定参数 C = 2.9979e8; % 电磁波速度 j=sqrt(-1); % 雷达参数 Fc = 10e9; % 载频 Lambda = C/Fc; % 波长 Tr = 5e-6; % 脉冲宽度 Br = 300e6; % 脉冲信号带宽 Kr = Br / Tr; % 脉冲调频斜率 Fs = 400e6; % 脉冲采样频率 PRF = 0.5e3; % 脉冲重复频率 Na = 1024; % 方位向采样点数 Nr = 4096; % 距离向采样点数 Beamangle_T = 6 * pi / 180; % 发射机波束宽度(弧度) Beamangle_R = 8 * pi / 180; % 接收机波束宽度(弧度) % 场景信息 Pcenter = [0,0,0,0]; % 场景中心点坐标[x,y,z,掩码] Ntarget = 1; % 散射点数目 Dx = 0.5; % 像素间隔 X 方向 （米） Dy = 0.5; % 像素间隔 Y 方向 （米） ArrayposTar = [Pcenter+[0*Dx,0*Dy,0,1]; % 散射点位置及掩码，当ArrayposTar(i, 4) == 0 时，该点不产生回波 Pcenter+[-10*Dx,0*Dy,0,1]; Pcenter+[10*Dx,0*Dy,0,1]; Pcenter+[0*Dx,-10*Dy,0,1]; Pcenter+[0*Dx,10*Dy,0,0]; Pcenter+[10*Dx,-10*Dy,0,0]; Pcenter+[-10*Dx,10*Dy,0,0]; Pcenter+[-10*Dx,-10*Dy,0,0]; Pcenter+[10*Dx,10*Dy,0,0];]; figure('Name','原始场景'); hold on for iii = 1 : Ntarget if (ArrayposTar(iii, 4) == 0) continue; end plot3(ArrayposTar(iii, 1), ArrayposTar(iii, 2), ArrayposTar(iii, 3), '*'); end title('原始场景');xlabel('距离向（米）');ylabel('方位向（米）'); ArrayRCSTar = ones(1, Ntarget); % 散射点散射系数，与ArrayposTar一一对应； % 发射机平台参数 Ht = 8000; % 发射机高度 Pt0 = [-10000, 0, Ht]; % 发射机初始位置 LOS_T0 = (Pt0 - Pcenter(1:3)) / norm(Pt0 - Pcenter(1:3));% 发射机视线方向，需归一化 Vt = [0, 170, 0]; % 发射机速度 % 接收机平台参数 Hr = 6000; % 接收机高度 Pr0 = [0, -24000, Hr]; % 接收机初始位置 LOS_R0 = (Pr0 - Pcenter(1:3)) / norm(Pr0 - Pcenter(1:3));% 发射机视线方向，需归一化 Vr = [0, 140, 0]; % %For Debug % Pr0 = Pt0; % LOS_R0 = LOS_T0; % Vr = Vt; TRDelay = (norm(Pt0 - Pcenter(1:3)) + norm(Pr0 - Pcenter(1:3))) / C - 0.3E-5; % 接收波门延时，单位：秒 % 初始化中间变量 pulse_len = round(Tr * Fs); if mod(pulse_len, 2) == 1 pulse_len = pulse_len + 1; end tfast = 1 : pulse_len; tfast = tfast - round(pulse_len / 2); tfast = tfast / Fs; Sr = zeros(Na,Nr); for ID_n = 1 : Na IDID = ID_n; Pt = Pt0 + Vt * (ID_n - 1 ) / PRF; % 发射第n个脉冲时发射机的位置 Pr = Pr0 + Vr * (ID_n - 1) / PRF; % 发射第n个脉冲时接收机的位置 LOS_T = LOS_T0; % 设置当前时刻发射机视线方向 LOS_R = LOS_R0; % 设置当前时刻接收机视线方向 for ID_Tar = 1 : Ntarget if ArrayposTar(ID_Tar, 4) == 0 % 如果散射点掩码为0，则不产生回波 continue; end Ptar = ArrayposTar(ID_Tar,1:3); % 将Ptarget中的位置信息赋值给Ptar TMP1 = Pt - Ptar; % 发射机与散射点连线的矢量 Rt = norm(TMP1); % 发射机与散射点的斜距 omiga = acos(abs(LOS_T * TMP1.') / Rt); % 计算本散射点与发射机视线(波束中心)方向夹角（单位为弧度），LOS必须事先归一化 TMP_T = omiga / (Beamangle_T / 2); if TMP_T > 4 % 如果本散射点与发射机视线方向夹角大于2倍波束角，则认为散射点在孔径以外 continue; end Gain = sinc(TMP_T * 0.4430); % 计算发射天线增益 sinc(0.4430) = sqrt(0.5) TMP1 = Pr - Ptar; % 接收机与散射点连线的矢量 Rr = norm(TMP1); % 接收机与散射点的斜距 omiga = acos(abs(LOS_R * TMP1.') / Rr); % 计算本散射点与接收机视线(波束中心)方向夹角（单位为弧度），LOS必须事先归一化 TMP_R = omiga / (Beamangle_R / 2); if TMP_R > 4 % 如果本散射点与接收视线方向夹角大于2倍波束角，则认为散射点在孔径以外 continue; end Gain = Gain * sinc(0.4430 * TMP_R); % 计算发射、接收天线增益 R = Rt + Rr; % 发射第n个脉冲时，第i个散射点到发射机和接收机的距离和(米) tao = R / C; % 发射第n个脉冲时，第i个散射点的回波延时(秒) TMP3 = round((tao - TRDelay) * Fs); % 计算散射点回波到接收波门距离，如果TMP3小于0 或 TMP3 + pulse_len 大于 Nr，则存在部分接收数据丢失； Residu = (tao - TRDelay) - TMP3 / Fs; TMP6 = tfast - Residu; TMP5 = exp(j * pi * Kr * TMP6 .^ 2) * exp(-j * 2 * pi * Fc * tao); TMP4 = zeros(1, Nr); if TMP3 >= Nr % 当回波完全在接收波门之后，则未接到回波 continue; elseif (pulse_len + TMP3) <= 1 % 当回波完全在接收波门之前，则未接到回波 continue; end if TMP3 < 0 % 当回波在接收波门之前，则接到部分回波 TMP4(1 : pulse_len + TMP3 + 1) = TMP5(-1 * TMP3 : pulse_len); elseif (TMP3 + pulse_len) > Nr % 当回波在接收波门之后，则接到部分回波 TMP4(TMP3 + 1: Nr) = TMP5(1 : Nr - TMP3); else % 当回波在接收波门之内，则接到全部回波 TMP4(TMP3 + 1 : TMP3 + pulse_len) = TMP5; end Sr(ID_n,:) = Sr(ID_n,:) + Gain * ArrayRCSTar(ID_Tar) * TMP4; % 产生回波 end end figure('Name','回波信号');imagesc(abs(real(Sr))); title('回波信号');xlabel('距离向（像素）');ylabel('方位向（像素）'); % ------------------------------------- BP成像 ------------------------------------------------ % Range_ref = exp(-j * pi * Kr * tfast .^ 2); % 距离压缩 % S1 = zeros(Na, Nr + pulse_len - 1); % for ID_n = 1 : Na % S1(ID_n,:) = conv(Sr(ID_n, :), Range_ref); % end % figure; % plot(S1); % figure('Name','距离压缩结果'); imagesc(abs(S1)); % title('距离压缩结果', 'fontsize',16);xlabel('距离向（像素）');ylabel('方位向（像素）'); t1=(linspace(0,Nr-1,Nr)-Nr/2)*Fs/Nr;% fr=circshift(t1,[0,Nr/2]); for ii=1:Na Sr(ii,:)=fft(Sr(ii,:)).*exp((j*pi*fr.^2)/Kr)/Nr; %加窗函数 Sr(ii,:)=circshift(ifft(circshift(Sr(ii,:).*4*pi.*(fr+Fc)/C,[0,-0/2]),Nr),[0,-0/2]); end figure('Name','距离压缩结果'); imagesc(abs(Sr)); title('距离压缩结果', 'fontsize',16);xlabel('距离向（像素）');ylabel('方位向（像素）'); Na2=1024; Nr2=4096; oc=10000; t=(linspace(0,Na-1,Na)-Na/2)/PRF; %从0到Na-1均匀分布的Na个数 gama=zeros(3,Na); %3*Na的0矩阵 gama(1,:)=(Vr(2)+Vt(2))*t; gama(2,:)=oc; gama(3,:)=0; z1=(linspace(0,Na2-1,Na2)-Na2/2)*4;%4; z2=(linspace(0,Nr2-1,Nr2)-Nr2/2)*400;%0.1; rindex=C/2*(linspace(0,Nr-1,Nr)-Nr/2)/Fs+oc; data=zeros(Na2,Nr2); for i=1:Na2 temp=zeros(1,Nr2); for j=1:Na qtr=zeros(1,Nr2); rtr=sqrt((z1(i)-gama(1,j))^2+(z2-gama(2,j)).^2+(gama(3,j))^2); temp=temp+interp1(rindex,Sr(j,:),rtr,'pchip').*exp(j*4*pi/Lambda*rtr); end data(i,:)=temp; end figure('Name','BP算法后得到的图像');imagesc(abs(data)); title('BP算法后得到的图像');xlabel('距离向');ylabel('方位向');