SAR成像-BP后向投影算法及GPU加速

最新推荐文章于 2024-06-24 16:01:03 发布

想到好名再改

最新推荐文章于 2024-06-24 16:01:03 发布

阅读量2.2k

点赞数 3

分类专栏： SAR 文章标签：算法 matlab SAR back projection

本文链接：https://blog.csdn.net/qq_44697805/article/details/129756958

版权

SAR 专栏收录该内容

1 篇文章

订阅专栏

本文介绍了使用MATLAB实现的SAR成像过程，特别是应用了后向投影算法（BP算法），并添加了详细注释和可视化图表。作者还优化了代码，利用GPU加速，测试显示速度提高了约50倍。内容涵盖了关键参数设置、信号处理和成像步骤。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

此文为SAR BP成像matlab实现，主要参考matlab实现-合成孔径雷达(SAR)-后向投影算法(BP算法)公式分析-完整代码-详解，并在其基础上增加更多注释和可视化图，并改写了GPU加速代码，经测试提速约为50倍。
此文旨在记录学习过程，有错误的地方欢迎指出，谢谢。

clear;
close all;
clc;

% 参数
use_gpu = false;% 是否使用GPU加速
Fc = 25.5e9;% 中心频率
Band = 4e9;% 带宽
TimeUp = 640e-6;% 升频率时间
TimeTotal = 1000e-6;% 一个chirp周期总时间
Fs = 25e5;% 采样频率  
% 一个chirp周期升频率采样点数
SamplePerChripUp = double(int32(TimeUp*Fs));% 1600   
C = physconst('lightspeed');% 光速

% 
Nr = SamplePerChripUp;

Tfast = (0:SamplePerChripUp-1)/SamplePerChripUp*TimeUp;
% size(Tfast) = 1 * 1600
% 0 到 0.6396e-3的1600个数，表示chirp升频过程中的时间下标

在这里插入图片描述
plot(Tfast)

% Kr表示调频率，即频率变化速度
Kr = Band/TimeUp;%  6.2500e+12                                  
dr = C/2/Band;% 距离分辨率
lambda = C/Fc;% 信号波长
theta = 120;% 雷达有效探照角度
RadarVelocity = 5;% 雷达运动速度
Timer = 4;% 整个过程持续时间
Na = Timer/TimeTotal;% 方位向位置总点数， 表示有多少个chirp周期
% Na = 4000

% 雷达位置，表示雷达的三维坐标，xyz，其中x为变化的，yz均为0
RadarPos = [(-Timer/2:TimeTotal:Timer/2-1/Fs)*RadarVelocity; zeros(1,Na); zeros(1,Na)];
% size(RadarPos) = 3 * 4000
% 目标点xz=0，y=10
TargetPos = [0,10,0];% 仿真点位置

% 雷达与仿真点距离，表示雷达和目标点的三维空间距离
Range = sqrt( ...
    (RadarPos(1,:)-TargetPos(1)).^2 ...
    +(RadarPos(2,:)-TargetPos(2)).^2 ...
    +(RadarPos(3,:)-TargetPos(3)).^2);

% acosd表示反余弦，./表示数组除法
% TargetPos(2) = 10
% 判断目标是否在探照面内，判断是否在120°范围内
isFull = acosd(TargetPos(2)./Range)<theta/2;                    

N_up = 10;                                                      % 升采样倍数
% tau表示波的双程延时，2倍距离/光速
tau = 2*Range/C;

在这里插入图片描述

雷达延X轴运动，目标为Y轴上的一个点。
在这里插入图片描述
plot(tau)

% 回波
% size(Tfast) = 1*1600
% size(tau) = 1*4000
td = Tfast-tau.';
% td = Tfast-tau.';
% size(td) = 4000*1600
% 4000表示共雷达每个位置的延时，1600表示该时刻的回波，每次的回波区别在于延时不同
% Kr表示调频率
Srnm = exp(1j*2*pi*(Fc*td+1/2*Kr*td.^2)).*isFull.';
% size(Srnm) = 4000 * 1600

SrnmOri = exp(1j*2*pi*(Fc*Tfast+1/2*Kr*Tfast.^2));              % 参考信号
% conj表示复共轭
% size(conj(Srnm)) = 4000 * 1600
% size(SrnmOri) = 1 * 1600
SrnmFs = conj(Srnm).*SrnmOri;                                   % 得到中频信号
% size(SrnmFs) = 4000 * 1600

% size(hamming(SamplePerChripUp)) = 1600 * 1
% Na = 4000
rwin = hamming(SamplePerChripUp)*ones(1,Na);                    % 汉明窗
% size(rwin) = 1600 * 4000

SrnmFs = SrnmFs.*rwin.';

在这里插入图片描述
加汉明窗之前

加汉明窗之后

% size(SrnmFs) = 4000 * 1600
% fft(X, n, dim)表示对X的第dim维进行傅里叶变换，n表示返回的点的数量
% SamplePerChripUp*N_up = 16000
SrnmFt = fft(SrnmFs,SamplePerChripUp*N_up,2);                   % 距离维FFT
% size(SrnmFt) = 4000 * 16000

dr = dr/N_up;                                                   % 更新距离分辨率
Nr = Nr*N_up;                                                   % 更新Nr
figure,plot((0:SamplePerChripUp*N_up-1)*dr,abs(SrnmFt(1000,:)));% 随便找个点看看

在这里插入图片描述
表示其中某一个时刻的中频

imagesc(abs(SrnmFt))
表示全部时刻的中频

% 开始构建成像区域
L = 0.1;                                                        % 基本距离单元（随便设得）
pix = 0.01;                                                     % 像素距离
x = -50*L:pix:50*L;                                             % 构建x轴，我这里是方位维
y = (-50*L:pix:50*L)+TargetPos(2);                              % 以目标为中心构建y轴，对应距离维
[X,Y] = meshgrid(x,y);                                          % 构建
img = zeros(size(x,2),size(y,2));                               % 预构图像
h = waitbar(0,'正在BP计算');                                     % 加个进度条

if use_gpu
    RadarPos = gpuArray(RadarPos);
    X = gpuArray(X);
    Y = gpuArray(Y);
    SrnmFt = gpuArray(SrnmFt);
end
t_start = cputime;
for ix = 1:Na% Na = 4000
    RadarPosNow = RadarPos(:,ix);                               % 找出现在的雷达位置
    Ran = sqrt((RadarPosNow(1)-X).^2+(RadarPosNow(2)-Y).^2 ...  % 计算当前时刻雷达和成像区域每一个点的距离
        +RadarPosNow(3).^2);
    Ran = Ran.*(acosd(Y./Ran)<=60);                             % 只考虑雷达照射范围内的点
    deltaN = round(Ran/dr) + 1;                                 % 计算距离对应的点数
    deltaN_ij = (deltaN > 0 & deltaN <= Nr);                    % 考虑合理区间
    deltaN = deltaN.*deltaN_ij + Nr*(1 - deltaN_ij);            % 考虑合理区间
    sig_rdta=SrnmFt(ix,:);                                      % 找出当前时刻雷达采集到的数据
    sig_rdta(Nr)=0;                                             % 找出合理区间
    img=img+sig_rdta(deltaN).*exp(1i*4*pi/lambda.*Ran);         % 将成像区域每个点在该时刻雷达数据中的点进行对应
    waitbar(ix/Na);                                             % 更新进度条
end
close(h);                                                       % 关闭进度条
figure,imagesc(abs(img));                                       % 成像看看

t_end = cputime;
t_total = t_end - t_start;
t_total