合成孔径雷达回波生成，距离多普勒算法前提

@Air-

已于 2024-08-19 20:17:34 修改

阅读量469

点赞数 11

文章标签： matlab

于 2024-08-19 20:15:56 首次发布

本文链接：https://blog.csdn.net/m0_56333457/article/details/141333733

版权

“革命的道路，同世界上一切事物活动的道路一样，总是曲折的，不是笔直的。”

1.1 回波生成

1.1.1 工作模式

正侧式条带SAR

1.1.2 参数设置

相关参数	大小
带宽	150MHz
中心频率	3.0GHz
载机速度	200m/s
载机高度	10km
下视角	45°
脉冲重复频率（PRF）	300Hz
脉冲宽度	4.0 $\mu$ s
方位向分辨率	1.0m
测绘带宽度( $T_p$ )	80m×80m

这部分matlab代码如下所示：

%% 已知参数设置
c = 3e8;        %光速
fc = 3e9;       %中心频率
Br = 150e6;     %带宽
rho_a = 1;      %方位向分辨率
v = 200;        %载机速度
H = 10e3;       %载机高度
Tp = 4e-6;      %脉冲宽度
phi = deg2rad(45);          %中心下视角
SquintAngle = deg2rad(0);   %斜视角
MapBandr = 80;  %距离向测绘带宽度
MapBanda = 80;  %方位向测绘带宽度

1.1.3 求解未知参数

分辨率公式如下：
$\begin{equation} \begin{cases} \rho_r = \frac{c}{2B}\\ \rho_a = \frac{\lambda}{2\theta_{BW}}\\ \end{cases} \end{equation}$
式中， $\theta_{BW}$ 称之为波束宽度。

二维平面几何关系示意图如下图所示，设场景中仅存在有一个散射点P，雷达自左向右飞行，发射波束照射散射点P，形成合成孔径 $L_s$ ,且近似有 $L_s = \theta_{BW}·R_0$ ，其中 $R_0$ 为零多普勒斜距。如果横向有多个散射点，合成孔径保持不变，但雷达运动长度变为 $L_s+$ 散射点横向场景宽度。

二维平面几何关系示意图

对多普勒带宽做如下推导：

众所周知，多普勒公式为： $f_d = \frac{2v}{\lambda}$ ，雷达位于A处时，雷达靠近目标运动，此时多普勒为正，水平方向的速度 $v$ 在雷达视线上的分量（LOS）为 $v·\rm{sin}(\theta_{BW}/2)\approx v·\theta_{BW}/2$ ，故此时多普勒为 $f_d = \frac{2v·\theta_{BW}/2}{\lambda} =\frac{v·\theta_{BW}}{\lambda}$ ；而雷达运动到B处时，多普勒分量为负，大小与A处一致，故多普勒带宽为 $B_d =\frac{2v·\theta_{BW}}{\lambda}$ 。
实际雷达工作中，发射机发出信号，过一段时间后，雷达接收机才能接收到信号，这段时间，称之为时延，通过这个时延便可计算出雷达到目标的粗略距离。我们已经知道，雷达接收回来的信号，要通过A/D采样，将模拟信号转换为数字信号，便于处理。如果雷达在在发射完信号就立即开始采样，势必会有很大程度的噪声影响，而且过多的数据量会使后端处理压力增大。所以，应该要在回波回来的同时，或者放宽一点时间去开启采样（即开启波门），这样便接收到的数据便包含有用信号，且同时减少了不必要的数据量。通过计算雷达到目标的最近、最远距离，便可以得到对应的时延，从而求解波门宽度。

我们仍然以几何图形来说明。
三维平面几何关系示意图

三维平面几何关系示意图

最近距离 $R_{near}$ 、最远距离 $R_{far}$ 如图所示，由此便可计算相关时延，进而求解波门。详细计算见后续代码。

此时，剩余参数计算完成，还有些许没有说明的参数计算，见代码说明。

此部分MATLAB代码如下所示。

%% 剩余参数确定
lambda = c/fc;              %中心波长
G = H*tan(phi);             %地距
R0 = sqrt(G.^2+H.^2);       %零多普勒平面斜距
Rsq = R0/cos(SquintAngle);  %斜距
rho_r = c/(2*Br);           %距离向分辨率
Kr = Br/Tp;                 %调频率
fs  = 1.2*Br;               %采样频率
theta_BW = lambda/(2*rho_a);%阵元波束宽度
theta_BW_deg = rad2deg(theta_BW);%弧度转度
Bd = 2*v/lambda*theta_BW;   %多普勒带宽
Ls = theta_BW * R0;         %有效孔径长度
PRF = 300;                  %脉冲重复频率1.2*Bd
PRI = 1/PRF;                %脉冲重复周期
Rnear = sqrt(H.^2 + (G - MapBandr/2).^2);
Rfar = sqrt(H.^2 + (G + MapBandr/2).^2)/cos(theta_BW/2);
Tstart = 2*Rnear/c;           %波门开启时间
Tw = (2*(Rfar - Rnear)/c+Tp); %距离向接收时宽
Tm = (Ls+MapBanda)/v;       %方位向接收时宽
Nr = round(fs*Tw);          %距离向采样点数
if mod(Nr,2) == 1 
    Nr = Nr + 1;
end
Na = round(PRF*Tm);         %方位向采样点数
if mod(Na,2) == 1 
    Na = Na + 1;
end
Ntt =round(fs*Tp);          % 发射信号采样点数
if mod(Ntt,2) == 1 
    Ntt = Ntt + 1;
end
tr = Tstart+(0:Nr-1)/fs;    %快时间域
tm = (-Na/2:Na/2-1)/PRF;    %慢时间域
tt = (0:Ntt-1)/fs;          %发射信号快时域
fr = (-Nr/2:Nr/2-1)*fs/Nr;  %距离频域
fa = (-Na/2:Na/2-1)*PRF/Na; %多普勒域
fdc = 0;
Xc = 0;                     %合成孔径中心
X = tm*v;                   %SAR所处位置
rr = (-Nr/2:Nr/2-1)*c/(2*fs*sin(phi));%距离向距离坐标
ra = tm*v;                  %方位向距离坐标

1.1.4 布置点目标

%% 布目标点阵
points_x = [-20 0 20];    % 散射点横坐标
points_y = [-20 0 20];    % 散射点纵坐标

target = sqrt((G+points_y).^2 + H.^2); %目标到雷达轨迹的最短路径
Pnum = length(points_x); %目标点数

figure;
scatter(points_x,points_y);axis([-MapBanda/2,MapBanda/2,-MapBandr/2,MapBandr/2]);
xlabel("X/m"),ylabel("Y/m");title("点阵布设");

在这里插入图片描述

点阵布设

1.1.5 回波生成

设雷达发射信号为

$\begin{equation} S_t(t)= A·w_r(t)·{\rm exp \{(j2\pi}(f_ct+\frac{1}{2}\gamma t^2)\} \end{equation}$

式中， $w_r(t) = {\rm{rect}}(\frac{t}{T_p})$ ，为距离窗函数。

设点目标P到雷达的距离为 $R_p(\eta)$ ，令 $\tau(\eta) = 2R_p(\eta)/c$ ，则雷达接收到的回波为

$\begin{equation} S_r(t,\eta)= A·w_r(t-\tau(\eta))·w_a(\eta-\eta _c)·{\rm{exp}}\bigg\{ {\rm {j2\pi}}\Big[f_c(t-\tau(\eta))+\frac{1}{2}\gamma (t-\tau(\eta))^2\Big]\bigg\} \end{equation}$

式中， $t、\eta$ 分别为快时间和慢时间， $w_a(\eta) = p_a^2[(\theta(\eta)]$ ，为方位窗函数。其中， $p_a(\theta) = {\rm{sinc}}\Big(\frac{0.886\,\theta}{\theta _{BW}} \Big)$ 。 $\theta$ 为斜距平面内测得的与视线的夹角。由此可见，方位窗与天线波束形状有关。但本人在仿真中，方位窗设置的矩形窗 😐。

正交解调之后回波为：
$\begin{equation} \begin{split} S_r(t,\eta)&= A·w_r(t-\tau(\eta))·w_a(\eta-\eta _c)·{\rm exp\{ j2\pi} [ f_c(-\tau(\eta))+\frac{1}{2}\gamma (t-\tau(\eta))^2]\}\\ &= A·w_r(t-\tau(\eta))·w_a(\eta-\eta _c)·{\rm exp\{-j2\pi} f_c·\tau(\eta)+{\rm j} \pi\gamma (t-\tau(\eta))^2]\} \end{split} \end{equation}$

当然，这只是假设只存在一个点目标的情况，存在多个点目标时，要在回波前多加一个求和符号 $\Sigma$ 。由此可列出此部分代码。

%% 目标回波
sr_tt = zeros(Nr,Na);     %定义存放时域回波数据的数组,距离向Nr,方位向Na
for p =1:Pnum
    for m = 1:Na
        eta = tm(m);            
        Rp = sqrt(target(p).^2 + (v*eta - points_x(p)).^2);   %慢时间为m时，目标距离
        tao = 2*Rp/c;                               %时延
        win_r = (abs(tr - Tp/2 - tao) < Tp/2).';    %距离窗
        win_a = (abs(X - Xc - points_x(p)) < Ls/2); %方位窗
        sr_tt(:,m) =sr_tt(:,m) + win_a(m).* win_r .* exp(-1i*2*pi*fc*tao + 1i*pi*Kr*(tr-tao-Tp/2).^2).';
        %sr_tt为二维时域回波，r表示回波，第一个t表示快时间域，第二个t表示慢时间域，后面同理
    end
end

figure;
imagesc(tm,tr*1e6,abs(sr_tt));axis xy;
xlabel("慢时间/s");ylabel("快时间/\mus");title("二维时域回波");
colorbar;

Sr_ft = (fft(sr_tt)); %距离向傅里叶变换
figure;
imagesc(tm,fr/1e6,fftshift(abs(Sr_ft),1));axis xy;
xlabel("慢时间/s");ylabel("频率/Hz");title("距离频域方位时域回波");
colorbar;

Sr_ff = (fft2(sr_tt));%二维频域
figure;
imagesc(fa,fr/1e6,fftshift(abs(Sr_ff)));axis xy;
xlabel("多普勒/Hz");ylabel("频率/Hz");title("二维频域回波");
colorbar;

Sr_tf = (ifft(Sr_ff));%默认对列向量求ifft
figure;
imagesc(fa,tr*1e6,fftshift(abs(Sr_tf),2));axis xy;
xlabel("多普勒/Hz");ylabel("快时间/\mus");title("距离多普勒回波");
colorbar;