SAR成像系列：【5】合成孔径雷达（SAR）成像算法-距离多普勒（RD）算法（附Matlab代码）

HIT夜枭

已于 2022-08-25 15:27:39 修改

阅读量1w

点赞数 13

分类专栏：合成孔径雷达（SAR）成像与目标识别文章标签：算法

于 2022-08-24 17:16:11 首次发布

本文链接：https://blog.csdn.net/zr2006_7/article/details/126485437

版权

合成孔径雷达（SAR）成像与目标识别专栏收录该内容

19 篇文章 402 订阅

订阅专栏

该博客详细介绍了合成孔径雷达（SAR）成像中的距离多普勒算法，包括距离压缩、距离走动矫正和方位压缩等步骤。针对大斜视角情况，提出了双曲线模型和二次距离压缩（SRC）以解决距离徙动问题。通过Matlab代码展示了算法的仿真过程，并给出了不同阶段的图像结果，如距离压缩、距离走动矫正和最终的方位压缩成像结果。

摘要由CSDN通过智能技术生成

完整的距离多普勒算法主要包括距离压缩、距离徙动矫正（矫正距离走动和距离弯曲）、方位压缩等步骤。其中距离走动矫正即可在时域进行也可在频域进行，而距离弯曲矫正一般在多普勒域进行。在距离多普勒域叫作RCMC是算法的主要特定，因此被称为“距离多普勒（RD）”算法。具体算法流程如下图。

（1）小斜视情况

前面已经给出里SAR原始信号模型及距离历程的数学模型，分别为

$s_{r,a}( \tau ,t)=A_{0}\omega _{r}(\tau -\frac{2R(t)}{c})\omega _{a}(t-t_{c})e^{j2\pi f_{c}(\tau -\frac{2R(t)}{c})+\frac{1}{2}K_{r}(\tau -\frac{2R(t)}{c})^{2}}$

$R_{t}=\sqrt{R_{0}^{2}+(V_{a}t)^{2}}\approx R_{0}+\frac{(V_{a}t)^{2}}{2R_{0}}$

根据驻定相位原理，距离FFT结果为：

$S_{r,a}( f_{\tau} ,t)=A_{0}W_{r}(f_{\tau})\omega_{a} (t-t_{c})exp(-j4\pi \frac{(f_{0}+f_{\tau})R(t) }{c})exp(-j\frac{\pi f_{\tau }^{2}}{K_{r}})$

距离频域匹配滤波函数:

$H_{1}(f_{\tau})=rect(\frac{f_{\tau }}{\left | K_{r}\right |T})exp(j\frac{\pi f_{\tau }^{2}}{K_{r}})$

距离压缩的输出为：

$s_{rc}( \tau ,t)=IFFT_{\tau}(S_{r,a}( f_{\tau} ,t)H_{1}(f_{\tau}))=A_{0}p_{r}[\tau -\frac{2R(t)}{c}]\omega_{a} (t-t_{c})exp(-j4\pi f_{0}\frac{R(t)}{c})$

其中，pr(.)为sinc函数。

方位频率为 $f_{\tau}=-K_{a}\tau$ .

进行方位FFT：

$S1_{r,a}(\tau ,f_{t})=FFT_{t}(s_{rc}( \tau ,t))=A_{0}p_{r}[\tau -\frac{2R_{rd}(f_{t})}{c}]W_{a} (f_{t}-f_{t_{c}})exp(-j4\pi f_{0}\frac{R_{0}}{c})exp(j\pi \frac{f_{t}^{2}}{K_{a}})$

第一个相位项为固定信息，在成像中一般没用，在干涉成像中非常重要；第二项为方位调制。

$R_{rd}(f_{t})$ 为距离多普勒域内的距离等式，根据距离的泰勒展开约等式和 $f_{\tau}=-K_{a}\tau$ ，得到

$R_{rd}(f_{t})=R_{0}+\frac{V_{a}^{2}}{2R_{0}}(\frac{f_{t}}{K_{a}})^{2}=R_{0}+\frac{\lambda ^{2}R_{0}f_{t}^{2}}{8V_{a}^{2}}$

上式中的第一项与R0有关，表明不同距离单元的距离走动不同，它随距离线性变化；第二项与方位频率有关，是频率的二次函数。对距离徙动矫正可以在距离多普勒域进行sinc插值完成。

完成sinc插值距离徙动矫正后的信号为：

$S2_{r,a}(\tau ,f_{t})=FFT_{t}(s_{rc}( \tau ,t))=A_{0}p_{r}[\tau -\frac{2R_{0}}{c}]W_{a} (f_{t}-f_{t_{c}})exp(-j4\pi f_{0}\frac{R_{0}}{c})exp(j\pi \frac{f_{t}^{2}}{K_{a}})$

方位压缩的频域参考函数为：

$H_{2} = exp(-j\pi \frac{f_{t}^{2}}{K_{a}})$

方位压缩后的信号变化到时域为：

$s2_{r,a}( \tau ,t)=IFFT_{t}(S2_{r,a}(\tau ,f_{t})H_{2})=A_{0}p_{r}[\tau -\frac{2R_{0}}{c}]p_{a}(t)exp(-j4\pi f_{0}\frac{R_{0}}{c})exp(j2\pi f_{t_{c}}t)$

（2）大斜视情况

随着SAR成像斜视角的增大，距离等式应采用更精确的双曲线模型 $R_{t}=\sqrt{R_{0}^{2}+(V_{a}t)^{2}}$ 。此时的距离多普勒域为：

$S1_{r,a}(\tau ,f_{t})=A\omega_{r}(\frac{1}{1-K_{r}Z}[\tau -\frac{2R_{0}}{cD(f_{t},V_{a})}])W_{a}(f_{t}-f_{t_{c}})exp(-j4\pi R_{0}f_{0}\frac{D(f_{t},V_{a})}{c})exp(j\pi K_{m}[\tau-\frac{2R_{0}}{cD(f_{t},V_{a})}]^{2})$

其中， $D(f_{t},V_{a})=\sqrt{1-\frac{c^{2}f_{t}^{2}}{4V_{a}^{2}f_{0}^{2}}}$ , $Z(R_{0},f_{t})=\frac{cR_{0}f_{t}^{2}}{2V_{a}^{2}f_{0}^{3}D^{3}(f_{t},V_{a})}$ 新的调频率为 $K_{m}=\frac{K_{r}}{1-K_{r}Z}$ .

距离多普勒域中的距离等式为：

$R_{rd}(f_{t})=\frac{R_{0}}{D(f_{t},V_{a})}$

将 $\frac{1}{1-K_{r}Z}$ 近似为1，得到

$S1_{r,a}(\tau ,f_{t})=A\omega_{r}([\tau -\frac{2R_{rd}(f_{t})}{c}])W_{a}(f_{t}-f_{t_{c}})exp(-j4\pi R_{0}f_{0}\frac{D(f_{t},V_{a})}{c})exp(j\pi K_{m}[\tau-\frac{2R_{rd}(f_{t})}{c}]^{2})$
$K_{src}=\frac{1}{Z}$ 为二次距离压缩（SRC）滤波器的调频率。

当斜视角较大时，距离多普勒域中的Km与Kr存在较大差异，此时使用Kr进行距离压缩会造成散焦。对距离多普勒域的距离徙动矫正进行改进，需要矫正的距离徙动量为：

$R_{rd}(f_{t})-R_{0}=R_{0}[\frac{1-D(f_{t},V_{a})}{D(f_{t},V_{a})}]$

大斜视的方位滤波器变为：

$H_{2} = exp(j4\pi R_{0}f_{0}\frac{D(f_{t},V_{a})}{c})$

如果在距离压缩时使用的调频率为Kr，则需按照真实调频率对匹配滤波器进行补偿，补偿的调频率为 $K_{src}$ 。使用 $K_{src}$ 进行第二次距离压缩称为“二次距离压缩（SRC）”,SRC的匹配函数为

$H_{src}=exp(-j\pi \frac{f_{\tau}^{2}}{K_{src}(R_{0},f_{t})})$

(3)算法仿真

下表为仿真参数：

载频	9.875GHz
雷达高度	500Km
斜视角	30°
带宽	30MHz
雷达速度	7100m/s
目标个数	5
脉冲持续时间	20us

仿真结果：

①回波信号实部

②距离压缩结果

③距离走动矫正结果

④距离弯曲矫正和二次距离压缩结果

⑤方位压缩成像结果

Matlab 部分代码：

%% 前斜视RD算法含距离走动矫正和距离徙动矫正（距离多普勒域RCMC与二次距离压缩SRC）
%%% 时域矫正距离走动频域矫正距离弯曲；算法适用于斜视角较大的情况
%%%Authed by Piaobo
clear;close all;clc;
SNR = -15; % 信噪比

c=3e8;
f0 = 9.875e9; % 雷达工作频率Hz
lamda = c/f0; % 雷达工作波长m
H = 500e3; % 高度
Yc=750e3; % 成像区域中线
R0 = sqrt(Yc^2+H^2); % 中心斜距
theta = asind(H/R0); % 下视角
alpha = 30; %波束斜视角
Rc = R0/cosd(alpha); %景中心到天线距离
Xc = R0*tand(alpha); %景中心的方位坐标
Br=30e6; % 带宽
Vr = 7100; % 雷达有效速度m/s
Tr =20e-6; % 脉冲持续时间s
Kr = Br/Tr; % 线性调频率
Fr = 1.2*Br; % 距离采样频率，1.2为过采样率
Ts = 1/Fr; % 距离采样时间间隔s

Nk = ceil((2 * 1000/ c + Tr) / Ts); %距离向前后500m
Nf = 2^nextpow2(Nk); % 距离向的采样点个数
tf_ori = [-Nf/2:1:Nf/2-1]*Ts; % 距离向采样时序
tf = [-Nf/2:1:Nf/2-1]*Ts+2*Rc/c; % 实际快时间采样值

La = 6; % 等效天线尺寸
Ls = R0*lamda/La*cosd(alpha); % 合成孔径时长度m,Ls=（0.886*R0*lamda）/(La*cos(Theta)）
Ta = Ls/Vr; % 目标照射时间s
% Ta = 0.8; % 目标照射时间s
% Ls =Ta*Vr;
Ka = -2 * Vr^2*(cosd(alpha)).^2 / (lamda * R0); % 方位多普勒调频率Hz
Ba = 2*Vr*cosd(alpha)/La;%多普勒带宽
% Ba=abs(Ka*Ta); % 多普勒频率调制带宽
PRF = ceil(1.2*Ba); % 方位采样率Hz
% PRF = 1000; % 方位采样率Hz
PRT = 1/PRF; % 方位向采样时间间隔s
Ns = 2^nextpow2((200/Vr+Ta)*PRF); % 方位向的采样点个数左右各100m
ts = [-Ns/2 : (Ns/2 - 1)] * PRT; % 方位向采样时序
% ts = [-Ns/2 : (Ns/2 - 1)] * PRT+Xc/Vr; % 实际方位向采样时序
fdc = 2*Vr*sind(alpha)/lamda;%多普勒中心频率式4.33
fdc_ref=fdc;%参考目标的多普勒中心频率
% 理论分辨率
rho_r=c/2/Br;
rho_a=Vr*PRT;
% rho_a=La/2;

% 目标参数
X0 = [Xc-100 Xc+100 Xc Xc-100 Xc+100]; % 目标1位置坐标
Z0 = [0 0 0 0 0];
Y0 = [Yc+100 Yc+100 Yc Yc-100 Yc-100];
NT=size(X0,2);

%目标各自的波束中心穿越时刻
eta_c=(X0-sqrt(H^2+Y0.^2)*tand(alpha))/Vr;

%%================================================================
%%生成回波信号
Sb = zeros(Ns,Nf);
sigma = 1; % 回波幅度
for ii=1:NT
R = sqrt((Vr*ts-X0(ii)).^2+Y0(ii).^2+(Z0(ii)-H).^2);
tau = 2*R/c;
Dfast = ones(Ns,1) * tf - tau' * ones(1, Nf);
phase = pi*Kr*Dfast.^2 - (2 * pi *f0 * tau') * ones(1,Nf);
Sb = Sb+sigma * exp(1j*phase) .* (abs(Dfast) <= Tr/2) .* ((abs(ts -eta_c(ii)) <=Ta/2)' * ones(1,Nf));
end
% Sb = awgn(Sb,SNR,0); % 回波加噪

figure
imagesc(real(Sb)),colormap(gray);

%% 距离向压缩
x0 = ones(Ns,1)*(exp(-1j*pi*Kr*(tf_ori).^2).* (abs(tf_ori) <= Tr/2)); % 距离向匹配函数
fftx1 = fftshift(fft(fftshift(x0.'))).'; % 距离向匹配函数FFT
fftSb = fftshift(fft(fftshift(Sb.'))).'; % 原信号FFT
y0 = fftshift(ifft(fftshift((fftSb.*fftx1).'))).'; % 距离向压缩后信号
%显示
ta = ts * Vr; % 方位向的距离序列
tr = tf * c / 2; % 快时间采样对应的距离域(单程距)
figure;imagesc(abs(y0));colormap(gray);
xlabel('距离向');ylabel('方位向');title('距离向压缩');

%% 距离走动矫正与多普勒中心补偿

具体代码参考：(39条消息) 合成孔径雷达(SAR)成像距离多普勒(RD)算法仿真-行业报告文档类资源-CSDN文库