调频连续波(FMCW)雷达目标轨迹提取

大白菜菜籽

已于 2022-10-16 11:25:33 修改

阅读量7.4k

点赞数 26

分类专栏：雷达 FMCW 文章标签： matlab 其他

于 2022-04-27 17:32:12 首次发布

本文链接：https://blog.csdn.net/weixin_42686221/article/details/124451239

版权

雷达同时被 2 个专栏收录

10 篇文章

订阅专栏

FMCW

2 篇文章

订阅专栏

调频连续波雷达目标轨迹提取

本篇文章内容主要是基于调频连续波雷达采集的人体运动回波，提取回波中的距离、角度信息，并对目标运动轨迹进行绘制。
主要涉及到的信号处理的相关知识有：脉冲压缩、加窗、MTI/MTD、多天线测角（DOA估计）等

调频连续波雷达简介

调频连续波雷达是一种发射的信号为调频连续波的雷达。调频连续波信号就是频率随时间变化，而且占空比为1的信号。
但是一直发射信号也要考虑一些其他的因素，所以有些调频连续波雷达还是可以设置成脉冲的形式，比如TI的XWR系列。
我们只需要关心雷达采集到的三维回波矩阵（快时间、慢时间、通道）就可以了，所以对于这个问题就不再这里深究了。
调频连续波的两种主要形式：三角波、锯齿波。

下面给出三角波示意图（图片来自百度）

图片来自百度
关于三角波与锯齿波在信号处理中的差距，对本文的目的影响不大，此处便不再赘述，有感兴趣者可自查。
该链接可参考：https://www.zhihu.com/question/401476711

发射信号参数配置

%% 雷达为24GHz的三角波调频连续波雷达，天线为一发两收。
C          = 3e8;                 %光速              m/s
Fc         = 24.125e9;            %中心频率          Hz
B          = 250e6;               %信号带宽          Hz
Tr         = 2e-3;                %重复周期          s
PRF        = 1/Tr;                %重复频率          Hz
Fs         = 500e3;               %快时间维采样频率  Hz
Ts         = 1/Fs;                %快时间维采样周期  s 
D          = 0.0145;              %天线间隔

回波信号形式

回波矩阵

本文的雷达只有两个接收天线，所以通道维的个数=2，分为channelA 和 channelB。
快时间维是每个信号每个周期的采样点数，采样点数=采样时间 X 采样频率=1000。
慢时间就是发射的脉冲个数（可以自己设定，根据每帧多少个脉冲，发射多少帧这样，不过不重要），本文脉冲数为	13838。

回波处理

 在得到回波矩阵之后就可以，对数据进行处理了。首先，为了方便处理，首先将三角波矩阵转变成了锯齿波。
（很明显只需要取快时间前500就可以了，但是这个时候，相当于采样率没变，采样点数少了一半，在计算距离分辨率的时候，我们再进行分析）

距离向脉冲压缩

为了得到目标距离信息，需要进行脉冲压缩。脉冲压缩主要有两种方法：匹配滤波 和 Dechirp+FFT。
前者针对LFM信号，后者针对混频后的中频信号。本文调频连续波雷达采用混频得到中频信号属于第二种形式。
所以只需要对回波快时间维进行FFT即可得到回波的一维距离像。
考虑到有环境背景噪声的影响，本文采取了加窗操作。

下面展示加窗前后距离像对比。为了提高频率分辨率，快时间维进行了2048点的FFT，并对回波位置进行了截取。
加窗前后对比

下面给出距离向处理代码

echoes=load('~');
channelA=echoes.echoes.channelA;
channelB=echoes.echoes.channelB;
channelB=imag(channelB)+1j*real(channelB);   %%由于雷达数据格式问题，B通道实部虚部需要互调
[Nr,Na]=size(channelB);
%三角波转锯齿波
channelAA=channelA(1:500,:);
channelBB=channelB(1:500,:);
clear channelB channelA

%%距离向FFT
w=hamming(500)*ones(1,13838);  %% 窗函数。
figure
channelAA_fw=fftshift(fft(channelAA.*w,2048,1),1);
dataA_FFT=channelAA_fw(minbin_range:maxbin_range,:);
channelBB_fw=fftshift(fft(channelBB.*w,2048,1),1);
dataB_FFT=channelBB_fw(minbin_range:maxbin_range,:);
imagesc(abs(dataA_FFT))
title("通道A加窗后")

得到距离信息之后需要计算实际距离，其实只需要知道距离分辨率就可以了。
在采样率为 ${F_s}$ 的情况下进行 ${N}$ 点的采样，可得频率分辨率为
$Δf=\frac{F_s}{N}$ 根据调频率 ${K=B/\frac{T_r}{2}}$ ，就可以计算时间分辨率 ${Δt=\frac{Δf}{K}}$ ，进而可以计算距离分辨率 ${ΔR=\frac{CΔt}{2}=\frac{C\frac{\frac{F_s}{N}}{B/\frac{T_r}{2}}}{2}=\frac{CT_rF_s}{4BN}}$ 需要注意：由于是三角波，所以调频率的计算只用半个周期。

MTI/MTD

MTI本质就是个频域滤波器：可以抑制固定目标和慢速杂波。进行双脉冲对消，直接慢时间维差分就可以了。
MTD是一个带通滤波器组，一般使用慢时间FFT来实现，也就意味着在多普勒域对目标进行锁定，提高检测性能

下面给出示例代码

%%MTI/MTD 脉冲对消
%两脉冲对消
figure
dataA_FFT_mti=diff(dataA_FFT,1,2);
dataB_FFT_mti=diff(dataB_FFT,1,2);
%%MTD处理：方位向FFT，然后在距离-多普勒图中搜索最大值
figure
%方位向分帧数，帧数决定方位FFT个数，影响速度分辨率，这里Nfft_A取了64
framenum=floor((Na-1)/Nfft_A);
dataA_MTD=zeros(size(dataA_FFT_mti,1),framenum);
dataB_MTD=zeros(size(dataA_FFT_mti,1),framenum);
for k=1:framenum
    dataA_mtd=fftshift(fft(dataA_FFT_mti(:,(k-1)*Nfft_A+1:k*Nfft_A),Nfft_A,2),2);
    dataB_mtd=fftshift(fft(dataB_FFT_mti(:,(k-1)*Nfft_A+1:k*Nfft_A),Nfft_A,2),2);
    dataA_RV=dataA_mtd(:,minbin_v:maxbin_v);
    dataB_RV=dataB_mtd(:,minbin_v:maxbin_v);
    [maxrow,maxcol]=find(abs(dataA_RV)==max(max(abs(dataA_RV))));
    dataA_MTD(:,k)=dataA_RV(:,maxcol);
    dataB_MTD(:,k)=dataB_RV(:,maxcol);
end
imagesc(abs(dataA_MTD))

下面给出MTI/MTD处理过程图
在这里插入图片描述
结果显示MTI+MTD后得到的距离像已经较为理想，下一步可以进行角度提取。

角度提取

对角度信息的提取主要有天线维比相：长短基线法、天线维FFT，超分辨DOA估计：MUSIC、ESPRIT，等具体理论知识自查哈...
由于本文只针对单目标，且只有两个平面天线，所以可以简单的进行双天线比相，得到平面角。

双通道比相具体步骤：
获得channelA距离向最大值位置，将该位置上对应的channelA和channelB的数据进行共轭相乘后取相角：
$\begin{array}{l} \Delta \varphi {\rm{ = }}angle(channel{A_{\max }}.*channelB_{\max }^*)\end{array}$ 得到的相角即为两个通道的相位角度差，根据几何关系可以得到来波方向： $\begin{array}{l}\theta = \arcsin (\frac{{\lambda \Delta \varphi }}{{2\pi d}}) \end{array}$ 下面给出代码示例

%%比相测角
target_theta=[];
target_range=[];
dataA_MTD_dB=20*log10(abs(dataA_MTD)/max(max(abs(dataA_MTD))));
for k=1:framenum
    [maxrow,maxValue]=find(abs(dataA_MTD(:,k))==max(abs(dataA_MTD(:,k))));
    if dataA_MTD_dB(maxrow,k)>-25
        target_range=[target_range rangRes*maxrow];
        target_theta=[target_theta angleEstimate(dataA_MTD(maxrow,k),dataB_MTD(maxrow,k))];
    end
end
%%%%   angleEstimate函数   %%%%%
%%   phi=angle(a.*conj(b));
%%   theta=asin(lamba* phi/2/pi/D)/pi*180;
%平滑滤波
target_theta=smooth(target_theta);
target_range=smooth(target_range);

轨迹绘制

得到目标距离、角度信息之后就可以绘制目标轨迹了，下面分别给出绘制动态极坐标以及直角坐标图，并保存为了gif。

%图像显示
filename='track.gif';
h=figure;
for i=1:length(target_range)
    polarplot(target_theta(i)*pi/180,target_range(i),'bo')
    thetalim([-90,90]);
    hold on
%     pause(0.001)
    drawnow 
      % Capture the plot as an image 
      frame = getframe(h); 
      im = frame2im(frame); 
      [imind,cm] = rgb2ind(im,256); 
     if i == 1 
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.05); 
      else 
          imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.05); 
     end
end
%直角坐标系
X=target_range.*sind(target_theta);
Y=target_range.*cosd(target_theta);
b=figure;
filename='track_xy.gif';
for i=1:length(target_range)
    plot(X(i),Y(i),'bo')
    xlabel('水平距离（m）')
    ylabel('垂直距离（m）')
    axis([-4 4 0 8])
    hold on
%     pause(0.001)
        drawnow 
      % Capture the plot as an image 
      frame = getframe(b); 
      im = frame2im(frame); 
      [imind,cm] = rgb2ind(im,256); 
     if i == 1 
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.05); 
      else 
          imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.05); 
     end
end