合成孔径雷达(SAR) RD算法 面目标成像Matlab仿真

一、概述

所谓面目标仿真，其实是点目标仿真的升级版，不了解点目标仿真的同学可以先看一下点目标仿真

https://blog.csdn.net/qq_30017409/article/details/133905960

在点目标仿真中，我们主要是通过距离多普勒算法，将回波数据成像为了多个点，但是这个点在内容上是没有任何意义的。而如果有许多点目标，每个点目标都有它的灰度值，组合在一起就能够达到成像的效果了。
例如，一张分辨率为300*300的SAR图像，代表着它有300行×300列的像素，而每个像素都有一个灰度值（注意将SAR图像与RGB三通道图像区分），也就是强度值；那么我们利用点目标成像的方法，将每个点的强度值赋到我们的点目标上，就能够实现面目标成像了。
重新梳理一下面目标成像的流程：首先有一张SAR图像->将SAR图像的每个像素点对应的灰度值，赋到一个点目标上->对许多点目标同时进行成像->得到一副由回波数据，进过距离多普勒算法（其他算法也一样）成像的图片，这张图片对应与我们预设的那张SAR图像
那么肯定有同学要问了，这把一幅图片生成为回波数据，再成像，最后得到的还是原来那一副图片，好像没有什么意义？其实不然，因为各种各样的原因，写论文或是做实验的时候，很多实验并不适合在实测（真实的雷达）数据上进行，转而使用这种仿真的方式，一样能够达到发论文，论证创新方法有效性的目的，所以面目标成像也算是基本功。
本文旨在基于Ian G. Cumming的《合成孔径雷达成像算法与实现》中第六章的距离多普勒算法处，机载雷达的数据进行面目标仿真，辅以讲解，完整代码请见文末

另外：作者也是SAR成像的初学者，代码和讲解当中也有许多瑕疵的存在，仅供读者参考

二、面目标成像的注意点

将一个点目标的成像看成是用打印机打印“一个点”，那么面目标成像就是用打印机打印“很多个点”，而这些点的强度值(灰度值)不相同，从而组成了我们看到的图像。在这个“打印”的过程当中，下面的问题需要注意，其他涉及的部分基本与点目标仿真相同

1.算力的问题

面目标成像十分消耗算力，但不是成像过程消耗算力，而是将一张图片转换成回波，生成回波的计算过程非常消耗算力，这里面涉及到大量的矩阵运算，所以在专业仿真当中一般都会用GPU运算进行优化，接下来的示例使用一张130*150的J20图像用于演示，一般的笔记本也可以进行仿真。

在Matlab中，我们可以在生成回波时使用：

• single函数转换数据为单精度，缩短运行时间
• gpuArray函数将数据转移到GPU上进行计算

%根据Naz*Nrg的点数生成回波
Target_num=Ntarget;
S_echo = zeros( Nrg,Naz);
S_echo=gpuArray(single(S_echo));%转换为单精度并放置到GPU上计算

2.距离向分辨率与方位向分辨率的问题

将面目标成像看做是打印目标点，那么需要非常注意目标点之间的“间距”。

• 在图像中，我们的每个像素点都是紧邻的，所以不存在间距的概念
• 但是在SAR成像中，我们的每个点目标之间的距离分别是“距离分辨率”与“方位分辨率”
  因此产生了不同，我们需要在代码中将像素点与点目标的坐标对应起来
  
  首先我们要计算rho_r和rho_a

%  合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率|La / 2

然后生成点目标的坐标

%% 点目标坐标设置
%  设置目标点相对于景中心之间的距离

img_0   = imread('./800600_SIM/J20.png');
img_0=imrotate(img_0, -6, 'crop');%顺时针旋转5度(可选)
img     = im2double(rgb2gray(img_0));
%img     = im2double(img_0);
w_r =rho_r;           % 图片img中单位像素对应的距离(m)
w_a=rho_a;
NPosition = zeros(numel(img),4);
% 坐标系原点设为参考点零多普勒时刻雷达星下点，x地距向，y方位向，z地表法线向
for m = 1:size(img,1)%行数
    for n = 1:size(img,2)%列数
        img(m,n)
        if(img(m,n)==0)
            continue;%跳过背景点
        end
        NPosition( (m-1)*size(img,2)+n,: ) = [(n-1)*w_a,(m-1)*w_r, 0, img(m,n) ];
    end
end
Range_Size=size(img,1);%求像素点大小数量
Azimuth_Size=size(img,2);
nonZeroRows = any(NPosition, 2);%去除全为0的行，也就是刚才的背景点
NPosition = NPosition(nonZeroRows, :);
Ntarget     = size( NPosition,1 );%像素点的个数
 Position_x_r = (NPosition(:,2)-max(NPosition(:,2))/2).';
 Position_y_a = (NPosition(:,1)-max(NPosition(:,1))/2).'; %点目标的坐标矩阵表示
 Position_Abs= NPosition(:,4).';%目标点的幅度

在这个过程中还需要注意跳过背景点，也就是说灰度值为0的点我们不需要进行仿真。 if(img(m,n)==0)

3.随机相位的问题

在SAR成像中，随机相位通常指的是由于各种因素（如运动误差、载荷、地形等）引起的回波信号中的不可预测的相位变化。这些随机相位误差会影响图像的质量，因为理论上，只有当回波信号的相位被准确补偿时，才能得到聚焦的图像1。
添加随机相位的目的是为了对这些误差进行补偿。在SAR/ISAR相位补偿技术中，通过基于图像评价函数的逐级逼近算法，可以对随机干扰相位进行补偿，以确保对未知信号中干扰相位的正确估计2。这样的补偿技术有助于提高图像的质量和准确性。
简而言之，随机相位是SAR成像过程中不可避免的，但通过适当的处理和补偿技术，可以最大限度地减少其对成像质量的负面影响。

在我们的仿真过程中也要对每个点目标添加一个随机生成的相位

%% 生成随机相位(0~2*pi)
% 生成符合正态分布的随机数
mu = pi; % 均值
sigma = pi; % 标准差
num_samples = Naz*Nrg; % 生成的随机数个数
random_numbers = mu + sigma * randn(num_samples, 1);
% 将随机数映射到0到2π之间
random_numbers_mapped = mod(random_numbers, 2*mu);

对每个点进行仿真的时候：

   %后向散射系数
    A0=Position_Abs(i).*exp(1j*random_numbers_mapped(i));

将随机相位组合进后向散射系数，乘到回波上

三、成像步骤与图解

1.回波生成

2.距离向频谱

3.距离压缩结果

4.方位向傅里叶变换结果(距离多普勒域)

5.徙动校正结果

6.方位压缩(成像结果)

左侧为imagesc图片，右侧为灰度图。可以直观地展示成像结果

四、完整代码

注意：由于代码使用的是书上的机载参数，所以仅进行了小斜视角的RD算法仿真

%{
%      机载SAR参数RDA面目标仿真                          
%       作用：读取图片，生成回波，进行RD成像                               
%            作者：CYAN                                  
%}
clc;clear;close all;
%代码格式化命令：MBeautify.formatCurrentEditorPage()

%% 机载轨道参数

% 姿态参数

theta_r_c = 1*pi/180; %斜视角，单位为弧度;
%计算等效雷达速度(卫星做圆周运动的线速度)
Vr = 150; %运行速度150m每秒

%景中心斜距R_eta_c和最近斜距R0，斜视角theta_rc由轨道高度、俯仰角、入射角计算得出
R_eta_c = 20e3; %景中心斜距

R0 = R_eta_c/ cos(theta_r_c);


%% 信号参数设置

%   电磁波参数
c = 3e+8; % 光速
Vs = Vr; % 卫星平台速度
Vg = Vr; % 波束扫描速度

Lr=1.5;%距离向天线尺寸——>椭圆的短轴
f0 = 5.3e+9; % 雷达工作频率
lambda = c / f0; %电磁波波长


%  距离向信号参数
Tr = 2.5e-6; % 发射脉冲时宽
Kr = 20e12; % 距离向调频率
Br = Tr*Kr; % 距离向信号带宽
alpha_os_r = 1.2; % 距离过采样率

Fr = alpha_os_r * Br; % 距离向采样率

%  方位向信号参数
alpha_os_a =1.25; % 方位过采样率(高过采样率避免鬼影目标)

delta_f_dop=80;%多普勒带宽80Hz
La = 0.886*2*Vs*cos(theta_r_c)/delta_f_dop;

Fa = alpha_os_a * delta_f_dop; % 方位向采样率
Ta = 0.886 * lambda * R_eta_c / (La * Vg * cos(theta_r_c)); %目标照射时间
L_La=R_eta_c*lambda*delta_f_dop/(2*Vr*cos(theta_r_c));
%  景中心点(原点)的参数
time_eta_c = -R_eta_c * sin(theta_r_c) / Vr; % 波束中心穿越时刻
f_eta_c = 2 * Vr * sin(theta_r_c) / lambda; % 多普勒中心频率

La = 0.886*2*Vs*cos(theta_r_c)/delta_f_dop;              % 实际天线长度       
                          

%  合成孔径参数
rho_r = c / (2 * Fr); % 距离向分辨率
rho_a = Vr/Fa; % 距离向分辨率|La / 2
theta_bw = 0.886 * lambda / Lr; % 方位向3dB波束宽度
theta_syn = Vs / Vg * theta_bw; % 合成角宽度(斜面上的合成角)
Ls = R_eta_c * theta_syn; % 合成孔径长度



%% 点目标坐标设置
%  设置目标点相对于景中心之间的距离

img_0   = imread('./800600_SIM/J20.png');
img_0=imrotate(img_0, -6, 'crop');%顺时针旋转5度(可选)
img     = im2double(rgb2gray(img_0));
%img     = im2double(img_0);
w_r =rho_r;           % 图片img中单位像素对应的距离(m)
w_a=rho_a;
NPosition = zeros(numel(img),4);
% 坐标系原点设为参考点零多普勒时刻雷达星下点，x地距向，y方位向，z地表法线向
for m = 1:size(img,1)%行数
    for n = 1:size(img,2)%列数
        img(m,n)
        if(img(m,n)==0)
            continue;%跳过背景点
        end
        NPosition( (m-1)*size(img,2)+n,: ) = [(n-1)*w_a,(m-1)*w_r, 0, img(m,n) ];
    end
end
Range_Size=size(img,1);%求像素点大小数量
Azimuth_Size=size(img,2);
nonZeroRows = any(NPosition, 2);%去除全为0的行，也就是刚才的背景点
NPosition = NPosition(nonZeroRows, :);
Ntarget     = size( NPosition,1 );%像素点的个数
 Position_x_r = (NPosition(:,2)-max(NPosition(:,2))/2).';
 Position_y_a = (NPosition(:,1)-max(NPosition(:,1))/2).'; %点目标的坐标矩阵表示
 Position_Abs= NPosition(:,4).';%目标点的幅度

%% 时间轴参数
Nrg = Range_Size*2; % 距离线采样点数
Naz = Azimuth_Size*2; % 距离线数
Trg = Nrg / Fr;Taz = Naz / Fa;
%距离向/方位向采样时间间隔
Gap_t_tau = 1 / Fr;Gap_t_eta = 1 / Fa;
%距离向/方位向采样频率间隔
Gap_f_tau = Fr / Nrg;Gap_f_eta = Fa / Naz;
%  时间轴变量
time_tau_r = 2 * R_eta_c / c + (-Trg / 2:Gap_t_tau:Trg / 2 - Gap_t_tau); % 距离时间变量
time_eta_a = time_eta_c+(-Taz / 2:Gap_t_eta:Taz / 2 - Gap_t_eta); % 方位时间变量
%  随着距离向时间变化的最近斜距；c/2是因为距离向上一个时间包含两次电磁波来回
R0_tau_r = (time_tau_r * c / 2) * cos(theta_r_c);
Ext_R0_tau_r = repmat(R0_tau_r.', 1, Naz); %扩展R0，用于计算变量Ka
%  频率变量
f_tau = (-Fr / 2:Gap_f_tau:Fr / 2 - Gap_f_tau); % 距离频率变量
%f_tau=f_tau-(round(f_tau/Fr))/Fr;%混叠方程
f_eta = f_eta_c + (-Fa / 2:Gap_f_eta:Fa / 2 - Gap_f_eta); % 方位频率变量
%f_eta=f_eta-(round(f_eta/Fa))/Fa;

%  时间轴
[Ext_time_eta_a,Ext_time_tau_r] = meshgrid(time_eta_a,time_tau_r); % 设置距离时域-方位时域二维网络坐标
%  频率轴
[Ext_f_eta,Ext_f_tau] = meshgrid( f_eta,f_tau); % 设置频率时域-方位频域二维网络坐标





%% 生成随机相位(0~2*pi)
% 生成符合正态分布的随机数
mu = pi; % 均值
sigma = pi; % 标准差
num_samples = Naz*Nrg; % 生成的随机数个数
random_numbers = mu + sigma * randn(num_samples, 1);

% 将随机数映射到0到2π之间
random_numbers_mapped = mod(random_numbers, 2*mu);


%% 生成回波S_echo
%根据Naz*Nrg的点数生成回波
Target_num=Ntarget;
S_echo = zeros( Nrg,Naz);
S_echo=gpuArray(single(S_echo));%转换为单精度并放置到GPU上计算
f = waitbar(0,'生成回波数据中...'); % waitbar显示进度条

for i = 1:Ntarget

    R0_Target = R0  + Position_x_r(i); %对每个目标计算瞬时斜距

    time_eta_c_Target =(Position_y_a(i) -R0_Target * tan(theta_r_c)) / Vr; %波束中心穿越时刻

    %  计算目标点的瞬时斜距
    R_eta = sqrt((R0_Target)^2+(Vr^2)*((Ext_time_eta_a - Position_y_a(i) / Vr).^2));
    
    %  距离向包络
    Wr = (abs(Ext_time_tau_r-2*R_eta/c) <= Tr / 2);
    %  方位向包络
    Wa = (abs(Ext_time_eta_a-time_eta_c_Target)) <=Ta/2 ;
    
    %   相位
    Phase = exp(-1j*4*pi*f0*R_eta/c) .* exp(+1j*pi*Kr*(Ext_time_tau_r - 2 * R_eta / c).^2);
    %后向散射系数
    A0=Position_Abs(i).*exp(1j*random_numbers_mapped(i));
    %  接收信号叠加
    S_echo_Target = A0*Wr .* Wa .* Phase;
    S_echo = S_echo + S_echo_Target;
    % 进度条
    str=['回波数据生成：',num2str(100*i/Target_num),'%'];% 百分比形式显示处理进程,不需要删掉这行代码就行

    %fprintf("当前回波点坐标:%.2f,%.2f\n",Position_x_r(i),Position_y_a(i));
    waitbar(i/Target_num,f,str)% 更新进度条，配合句柄使用
end

S_echo = S_echo .* exp(-2j*pi*f_eta_c*Ext_time_eta_a); %校正时间轴



%% 距离压缩

Hf = (abs(Ext_f_tau) <= Br / 2) .* exp(+1j*pi*Ext_f_tau.^2/Kr);%滤波器
%  匹配滤波
S1_ftau_eta = fftshift(fft(fftshift(S_echo, 1), Nrg, 1), 1);
S1_ftau_eta = S1_ftau_eta .* Hf;
S1_tau_eta = fftshift(ifft(fftshift(S1_ftau_eta, 1), Nrg, 1), 1);

%% 方位向傅里叶变换
S2_tau_feta = fftshift(fft(fftshift(S1_tau_eta, 2), Naz, 2), 2);

%% 距离徙动校正RCMC:采用相位补偿法

%虽然Ka是随着R0变化的，但是在相位补偿时需要假设R0是不变的
delta_R = (((lambda * Ext_f_eta).^2) .* R0) ./ (8 * (Vr^2)); %距离徙动表达式
G_rcmc = exp((+4j * pi .* Ext_f_tau .* delta_R)./c); %补偿相位
S3_ftau_feta = fftshift(fft(fftshift(S2_tau_feta, 1), Nrg, 1),1); %在方位向傅里叶变换的基础上进行距离向傅里叶变换

S3_ftau_feta = S3_ftau_feta.*G_rcmc ; %与补偿相位相乘
S3_tau_feta_RCMC = fftshift(ifft(fftshift(S3_ftau_feta, 1), Nrg, 1), 1); %距离向傅里叶逆变换

%距离徙动校正结束

%% 方位压缩
%  根据变化的R0计算出相应的Ka矩阵(距离向变化，方位向不变)
Ka = 2 * Vr^2 * cos(theta_r_c)^2 ./ (lambda * Ext_R0_tau_r);
%  方位向匹配滤波器
Haz = exp(-1j*pi*(Ext_f_eta.^2)./Ka);
Offset = exp(-1j*2*pi*Ext_f_eta.*time_eta_c);%偏移滤波器，将原点搬移到Naz/2的位置，校准坐标
%  匹配滤波
S4_tau_feta = S3_tau_feta_RCMC .* Haz.*Offset;
S4_tau_eta = fftshift(ifft(fftshift(S4_tau_feta, 2), Naz, 2), 2);




%% 回波图可视化

figure('name', "回波可视化")
subplot(2, 2, 1);
imagesc(real(S_echo));
title('原始仿真信号实部');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');

subplot(2, 2, 2);
imagesc(imag(S_echo));
title('原始仿真信号虚部');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');


subplot(2, 2, 3);
imagesc(abs(S_echo));
title('原始仿真信号幅度');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');

subplot(2, 2, 4);
imagesc(angle(S_echo));
title('原始仿真信号相位');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');

%% 距离压缩可视化

%距离频谱可视化
figure('name', "回波频谱可视化")
subplot(2, 2, 1);
imagesc(real(S1_ftau_eta));
title('回波距离向实部');
ylabel('距离向频谱点数f_\tau');
xlabel('方位向时间 \eta');

subplot(2, 2, 2);
imagesc(imag(S1_ftau_eta));
title('回波距离向虚部');
ylabel('距离向频谱点数f_\tau');
xlabel('方位向时间 \eta');

subplot(2, 2, 3);
imagesc(abs(S1_ftau_eta));
title('回波距离向频谱幅度');
ylabel('距离向频谱点数f_\tau');
xlabel('方位向时间 \eta');

subplot(2, 2, 4);
imagesc(angle(S1_ftau_eta));
title('回波距离向相位');
ylabel('距离向频谱点数f_\tau');
xlabel('方位向时间 \eta');

%距离压缩结果
figure('name', "距离压缩时域结果")
subplot(1, 2, 1);
imagesc(real(S1_tau_eta));
title('实部');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');

subplot(1, 2, 2);
imagesc(abs(S1_tau_eta));
title('幅度');
ylabel('距离向时间 \tau');
xlabel('方位向时间 \eta');

%% 方位向傅里叶变换可视化

figure('name', "方位向傅里叶变换结果")
subplot(1, 2, 1);
imagesc(real(S2_tau_feta));
title('实部');
ylabel('距离向时间\tau');
xlabel('方位向频谱点数f_\eta');

subplot(1, 2, 2);
imagesc(abs(S2_tau_feta));
title('幅度');
ylabel('距离向时间\tau');
xlabel('方位向频谱点数f_\eta');

%% 距离徙动校正可视化

figure('name', "距离徙动校正结果")
subplot(1, 2, 1);
imagesc(real(S3_tau_feta_RCMC));
title('实部');
ylabel('距离向时间\tau');
xlabel('方位向频率点数 f_\eta');

subplot(1, 2, 2);
imagesc(abs(S3_tau_feta_RCMC));
title('幅度');
ylabel('距离向时间\tau');
xlabel('方位向频率点数 f_\eta');

%% 回波成像
figure('name', "点目标成像结果")
subplot(1, 2, 1);
imagesc(abs(S4_tau_eta));
title('幅度');
ylabel('距离向时间\tau');
xlabel('方位向时间 \eta');

subplot(1, 2, 2);
imshow(abs(S4_tau_eta)./max(max(abs(S4_tau_eta))));
title('幅度');
ylabel('距离向时间\tau');
xlabel('方位向时间 \eta');