万字长文详解：毫米波雷达物体检测技术【硬件&软件】
本文来源：无人车情报局
/ 导读 /
本文由汽车电子与软件授权发布，作者为William。
自动驾驶的技术路线上，对激光雷达的选用各有权衡。但无论是“强感知+强智能”路线还是“弱感知+超强智能”路线，毫米波雷达可以说是一个必选项，这个传感器有什么吸引人的优势和不得不面对的弱势呢？本文从硬件到软件层面，对毫米波雷达的检测技术进行详解。
硬件
Why Radars?
自动驾驶的感知系统是由多个传感器的输入组成, 每个传感器都有其自身的优缺点, 我们使用雷达来补偿其他传感器的缺点. 雷达最大的优点就是它可以直接测量径向速度, 尽管激光, 摄像头等其他传感器需要通过求导来计算速度, 雷达却可以直接确定速度, 这也有助于区分静止物体和运动物体. 另一个优点是雷达很少受到环境因素的影响, 可以在复杂光照, 雨雪等条件下全天候工作.
FMCW automotive radars
雷达早在20世纪30年代就已经投入使用, 当时它们主要被军方用于探测飞机, 但自那以后, 雷达技术已经取得了长足的进步. 如今, 它们越来越多地被用作先进驾驶辅助系统(ADAS)中的汽车雷达传感器. 下面的图片显示了雷达传感器和其他一整套传感器在整车中的使用情况:
整车传感器配置(source : cdn.rohde-schwarz.com)
射频技术和数字信号处理技术的进步使得低成本、小尺寸、高效率的雷达成为可能.
雷达具备精确测量远距离目标准确速度和空间信息的能力, 同时无惧恶劣光照和天气的影响, 因此成为自动驾驶应用的重要传感器. 目前大多数汽车都配备了用于自适应巡航控制ACC的雷达传感器, 用来检测车辆过近,保持距离和跟随前车速度. 随着自动驾驶级别的提高, 越来越多的雷达传感器被使用.
雷达传感器发展趋势
如下图所示, 雷达利用电磁波的传输和探测进行工作:
电磁波遇到障碍物时会被反射. 如果这些反射波在它们的原点再次被接收, 那么这意味着一个障碍物在传播方向上. 用于雷达的电磁能量的频率不受黑暗的影响, 也能穿透雾和云层. 这使得雷达系统能够确定由于距离、黑暗或天气等原因肉眼看不见的道路目标的位置. 现代雷达可以从目标回波信号中提取比距离信号更多的信息.
雷达原理(source: electricalelibrary)
汽车使用的雷达不像军事雷达体积那么庞大, 是小型传感器, 可以很容易地安装在前格栅或保险杠下面. 如下图所示, 一个雷达组件由不同部分组成.
博世汽车雷达
雷达天线罩: 雷达天线罩是一种结构性的、防风雨的外壳, 用来保护雷达天线. 天线罩的材料可以避免过多削弱发送的电磁波或天线接收到的电磁波, 对无线电波来说实际上是透明的.
雷达印刷电路板: 这是模拟硬件, 包括雷达收发信机和无线电波产生所需的天线.
印刷电路板和处理: 这包括数字信号处理(DSP)单元
汽车领域的雷达应用在自动驾驶领域的毫米波雷达主要有3个频段, 分别是24GHz, 77GHz和79GHz, 普遍分为长距雷达和短距雷达, 长距雷达可以测量远达250m的物体, 但视野很小, 而短距雷达被设计为具有更大的视野, 但看不到远处.
Radar vs Lidar
雷达与激光雷达的比较[source: cleantechnica ]
Radar和Lidar的之争是话题热点, 但实际上, 这些传感器是互补的. 激光雷达能够利用目标对激光的反射产生高分辨率成像. 能够提供精确的形状和位置信息, 但是激光雷达在恶劣的天气条件下无法工作, 因为很小的波长不允许在雾或雨中很好地工作. 此外, 激光雷达是一种昂贵的传感器, 成本从35,000美元到100,000美元不等, 截至2019年. 激光雷达是 Waymo 的首选技术. 雷达不具备产生高分辨率图像的能力, 但是它具有基于多普勒现象的高精度速度估计. 此外, 雷达波长允许它在恶劣的天气条件下感知目标. 最重要的是低制造成本的雷达. 一个雷达单元的成本可以低至几百美元, 这使得汽车制造商可以部署多个360度感知的雷达传感器. 特斯拉依赖于雷达作为其主要传感器, 并没有包括激光雷达在其传感器融合系统. 也就是说有钱多多益善, 没钱精挑细选.
Radar Properties
接下来你将看到看到信号特性的一般概述, 包括信号波长的定义和信号方程的一般形式。
单波参数
波长( lambda)是从一个波的一点到下一个波的同一点的物理长度, 计算方法如下:
频率(frequency)是每秒经过的波的数量, 以赫兹(GHz)为单位. 频率越高, 波长越短. 汽车雷达一般工作在 w 波段(76GHz - 81GHz) . 这个频率的信号称为毫米波, 因为波长是毫米. 一个信号的带宽(Bandwidth)是连续频带中最高和最低频率分量之间的差值.
振幅(Amplitute)是信号的强度. 它通常定义为射频信号 / 电磁场的功率( dB/dBm), 与配置雷达输出功率和感知接收信号相关. 雷达信号的幅度越大, 雷达的可见度越高. 汽车雷达可以工作在最大55dBm 输出功率(316w).
信号的频率
信号的振幅
信号的相位
相位(Phase)是波形周期中的一个特定时间点, 以角度为单位测量. 一个完整的循环周期是360°. 频率也可以定义为相位相对于时间、频率和信号相位的一阶导数.
这一特性将用于测量运动目标的多普勒频移. 下图为相位显示图.
正弦波相位
两个周期性信号的相位差称为相位差. 当差值为零时, 两个信号称处于同相位, 否则它们彼此处于同相位.
波的一般方程
A 是信号的振幅
f_c 是信号频率
ϕ 是信号的相位
FMCW
FMCW 调频连续波雷达(Frequency-Modulated Continuous Wave radar)是一种发射连续传输能量的雷达传感器. FMCW 雷达具有测量非常近目标的能力以及同时测量目标距离和相对速度的能力, 这使其成为汽车用雷达的首选类型.
FMCW波形
调频连续波(FMCW)是一种频率随时间增减的信号. 如上图所示, 它们也被称为upramps和downramps. 调频连续波雷达最常用的两种波形方式是锯齿波和三角波. 锯齿波形一般只使用上行坡道, 而三角形波形使用上行坡道和下行坡道.
我们经常需要讨论一个叫做啁啾(Chirp)的雷达信号, 它频率随时间增加或者减少, 频率增加的叫做up-chrip,频率下降的叫做down-chrip. 很多雷达使用重复的啁啾信号来计算距离和速度.啁啾的斜率有频率带宽B_sweep(y轴)和线性调频时间T_s(x轴)决定.
雷达的距离分辨率取决于B, 而最大的速度能力取决于T_s.
一个啁啾序列包括很多个chrip, 为了精确测量多普勒速度, 对每个啁啾进行多次采样以进行多次距离测量, 并由雷达进行传输.
FMCW硬件概述
FMCW 雷达的硬件实现
频率合成器: 在汽车雷达中, 频率合成器是产生频率使线性调频频率一直达到77ghz 的元件.
天线: 天线将电能转换成电磁波, 电磁波在空中辐射, 击中目标, 然后反射回雷达接收天线. 天线通过将能量聚焦在所需方向来增加信号的强度. 此外, 天线方向图决定了雷达的视场.
混频器: 在 FMCW 雷达中, 混频器混频器将回波信号与频率合成器产生的扫频信号相减, 获得频率Delta-也称为频移或中频信号(IF).
IF = Synthesizer Frequency - Return Signal Frequency
处理器: 处理器是所有数字信号处理、检测、跟踪、聚类和其他算法发生的处理单元. 该单元可以是一个微控制器, 甚至是一个 FPGA.
天线方向图
天线方向图
天线方向图是天线发射的相对场强的几何方向图. 天线的波束宽度决定了雷达传感器的视场. 如果对雷达的要求仅仅是探测自己车道上的目标, 那么波束宽度必须足够小, 以覆盖整个车道达到预期的范围. 如果波束宽度大于车道宽度, 它也能感知其他车道的目标. 天线辐射不仅包括主波束, 还包括副瓣. 天线旁瓣非常关键, 因为它们可能会产生错误警报, 并从不希望的方向挑选干扰. 如图所示, 天线的旁瓣指向不同的方向, 可以感知不在主波束内的目标. 为了避免旁瓣检测, 必须将旁瓣电平从主波束的峰值压制到30dB 以上.
天线类型
在77GHz上最常用的汽车雷达天线类型是贴片天线. 贴片阵列天线具有成本低、制造简单、低剖面等优点, 是汽车雷达应用的理想选择.
贴片阵列天线
Radar Cross Section（RCS） Overview
一个目标反射雷达能量的大小和能力由一个单项定义, 称为 sigma, 它的单位是 m ^ 2平方米, 被称为雷达截面积. 这个单位显示了雷达散射截面是一个区域. 目标雷达截面面积取决于:
目标的物理几何形状和外部特征
光滑的边缘或表面会使波向各个方向散射, 从而降低雷达散射截面. 然而, 锐角会使回波信号聚焦到源头的方向上, 从而导致较高的 RCS.
雷达照射的方向
雷达发射器的频率
用于汽车、卡车、自行车的材料, 有时甚至用于行人的服装材料
目标车辆的RCS
RCS单位转换
RCS 也可以使用对数值(dB)定义, 以下为转换方式.
不同目标具有不同的RCS, 我们可以对其进行转换.
雷达作用距离方程
利用雷达作用距离方程可以设计出满足作用距离要求的雷达发射器、接收器和天线, 使其具有理想的功率、增益和噪声性能.
R - 雷达探测目标的最大距离
P_s - 雷达发射功率(dBm)
G - 发射 / 接收天线(dBi)的增益
λ - 信号的波长(m)
σ - 雷达截面(m^2)
P_E - 雷达可以探测到的最小功率
软件
距离估算（Range Estimation）
R_max为最大探测距离.
delta_r为距离分辨率.
多普勒距离估算（Doppler Estimation）
多普勒效应（Doppler effect）
fD: 多普勒效应导致的传输频率频移
νr: 与目标的相对速度
λ : 信号的波长
多普勒相移（Doppler phase shift）
% Doppler Velocity Calculationc = 3*10^8; %speed of lightfrequency = 77e9; %frequency in Hz% TODO : Calculate the wavelengthwavelength = c/frequency;% TODO : Define the doppler shifts in Hz using the information from abovedoppler_shift= [3e3 -4.5e3 11e3 -3e3];% TODO : Calculate the velocity of the targets fd = 2*vr/lambdavr=doppler_shift*wavelength/2;% TODO: Display resultsdisp(vr)
快速傅里叶变换 (FFT)
1D FFT
%% RANGE MEASUREMENT% *%TODO* :%reshape the vector into Nr*Nd array. Nr and Nd here would also define the size of%Range and Doppler FFT respectively.Mix_reshape = reshape(Mix,[Nr,Nd]);% *%TODO* :%run the FFT on the beat signal along the range bins dimension (Nr) and%normalize.sig_fft1=fft(Mix_reshape,Nr)./Nr;% *%TODO* :% Take the absolute value of FFT outputsig_fft1=abs(sig_fft1);% *%TODO* :% Output of FFT is double sided signal, but we are interested in only one side of the spectrum.% Hence we throw out half of the samples.sig_fft1=sig_fft1(1:Nr/2);%plotting the rangefigure ('Name','Range from First FFT')subplot(2,1,1)% *%TODO* :% plot FFT outputplot(sig_fft1);axis ([0 200 0 1]);
2D FFT
%% RANGE DOPPLER RESPONSE% The 2D FFT implementation is already provided here. This will run a 2DFFT% on the mixed signal (beat signal) output and generate a range doppler% map.You will implement CFAR on the generated RDM% Range Doppler Map Generation.% The output of the 2D FFT is an image that has reponse in the range and% doppler FFT bins. So, it is important to convert the axis from bin sizes% to range and doppler based on their Max values.Mix=reshape(Mix,[Nr,Nd]);% 2D FFT using the FFT size for both dimensions.sig_fft2 = fft2(Mix,Nr,Nd);% Taking just one side of signal from Range dimension.sig_fft2 = sig_fft2(1:Nr/2,1:Nd);sig_fft2 = fftshift (sig_fft2);RDM = abs(sig_fft2);RDM = 10*log10(RDM) ;%use the surf function to plot the output of 2DFFT and to show axis in both%dimensionsdoppler_axis = linspace(-100,100,Nd);range_axis = linspace(-200,200,Nr/2)*((Nr/2)/400);figure,surf(doppler_axis,range_axis,RDM);
杂波（Clutter）
恒虚警率（CFAR）
Cell Averaging CFAR (CA-CFAR)
Ordered Statistics CFAR (OS CFAR)
Maximum Minimum Statistic (MAMIS CFAR)
And, multiple variants of CA-CFAR
CA-CFAR
% Implement 1D CFAR using lagging cells on the given noise and target scenario.% Close and delete all currently open figuresclose all;% Data_pointsNs = 1000;% Generate random noises=abs(randn(Ns,1));%Targets location. Assigning bin 100, 200, 300 and 700 as Targets with the amplitudes of 8, 9, 4, 11.s([100 ,200, 350, 700])=[8 15 7 13];%plot the outputplot(s);% TODO: Apply CFAR to detect the targets by filtering the noise.% 1. Define the following:% 1a. Training CellsT = 12;% 1b. Guard CellsG = 4;% Offset : Adding room above noise threshold for desired SNRoffset=5;% Vector to hold threshold valuesthreshold_cfar = [];%Vector to hold final signal after thresholdingsignal_cfar = [];% 2. Slide window across the signal lengthfor i = 1:(Ns-(G+T+1))% 2. - 5. Determine the noise threshold by measuring it within the training cellsnoise_level =sum(s(i:i+T-1));% 6. Measuring the signal within the CUTthreshold = (noise_level/T)*offset;threshold_cfar=[threshold_cfar,{threshold}];signal=s(i+T+G);% 8. Filter the signal above the thresholdif (signal<threshold)signal=0;endsignal_cfar = [signal_cfar, {signal}];end% plot the filtered signalplot (cell2mat(signal_cfar),'g--');% plot original sig, threshold and filtered signal within the same figure.figure,plot(s);hold on,plot(cell2mat(circshift(threshold_cfar,G)),'r--','LineWidth',2)hold on, plot (cell2mat(circshift(signal_cfar,(T+G))),'g--','LineWidth',4);legend('Signal','CFAR Threshold','detection')
2D CA-CFAR
% *%TODO* :%Select the number of Training Cells in both the dimensions.Tr=10;Td=8;% *%TODO* :%Select the number of Guard Cells in both dimensions around the Cell under%test (CUT) for accurate estimationGr=4;Gd=4;% *%TODO* :% offset the threshold by SNR value in dBoff_set=1.4;% *%TODO* :%design a loop such that it slides the CUT across range doppler map by%giving margins at the edges for Training and Guard Cells.%For every iteration sum the signal level within all the training%cells. To sum convert the value from logarithmic to linear using db2pow%function. Average the summed values for all of the training%cells used. After averaging convert it back to logarithimic using pow2db.%Further add the offset to it to determine the threshold. Next, compare the%signal under CUT with this threshold. If the CUT level > threshold assign%it a value of 1, else equate it to 0.% Use RDM[x,y] as the matrix from the output of 2D FFT for implementing% CFARRDM = RDM/max(max(RDM)); % Normalizing% *%TODO* :% The process above will generate a thresholded block, which is smaller%than the Range Doppler Map as the CUT cannot be located at the edges of%matrix. Hence,few cells will not be thresholded. To keep the map size same% set those values to 0.%Slide the cell under test across the complete martix,to note: start point%Tr+Td+1 and Td+Gd+1for i = Tr+Gr+1:(Nr/2)-(Tr+Gr)for j = Td+Gd+1:(Nd)-(Td+Gd)%Create a vector to store noise_level for each iteration on training cellsnoise_level = zeros(1,1);%Step through each of bins and the surroundings of the CUTfor p = i-(Tr+Gr) : i+(Tr+Gr)for q = j-(Td+Gd) : j+(Td+Gd)%Exclude the Guard cells and CUT cellsif (abs(i-p) > Gr || abs(j-q) > Gd)%Convert db to powernoise_level = noise_level + db2pow(RDM(p,q));endendend%Calculate threshould from noise average then add the offsetthreshold = pow2db(noise_level/(2*(Td+Gd+1)*2*(Tr+Gr+1)-(Gr*Gd)-1));%Add the SNR to the thresholdthreshold = threshold + off_set;%Measure the signal in Cell Under Test(CUT) and compare againstCUT = RDM(i,j);if (CUT < threshold)RDM(i,j) = 0;elseRDM(i,j) = 1;endendendRDM(RDM~=0 & RDM~=1) = 0;
Angle of Arrival(AOA)
Φ= incremental phase shift
d= spacing between antenna elements
θ= steering direction from the normal of the antenna surface
λ= wavelength of the signal
Clustering
%%%% *|clusterDetections|*%% This function merges multiple detections suspected to be of the same% vehicle to a single detection. The function looks for detections that are% closer than the size of a vehicle. Detections that fit this criterion are% considered a cluster and are merged to a single detection at the centroid% of the cluster. The measurement noises are modified to represent the% possibility that each detection can be anywhere on the vehicle.% Therefore, the noise should have the same size as the vehicle size.%% In addition, this function removes the third dimension of the measurement% (the height) and reduces the measurement vector to [x;y;vx;vy].functiondetectionClusters =clusterDetections(detections, vehicleSize)N = numel(detections);distances = zeros(N);for i = 1:Nfor j = i+1:Nif detections{i}.SensorIndex == detections{j}.SensorIndexdistances(i,j) = norm(detections{i}.Measurement(1:2) - detections{j}.Measurement(1:2));elsedistances(i,j) = inf;endendendleftToCheck = 1:N;i = 0;detectionClusters = cell(N,1);while ~isempty(leftToCheck)% Remove the detections that are in the same cluster as the one under% consideration%TODO : Complete the clustering loop based on the implementation%discussed in the lessonunderConsideration = leftToCheck(1);clusterInds = (distances(underConsideration,leftToCheck) < vehicleSize);detInds = leftToCheck(clusterInds);clusterDets = [detections{detInds}];clusterMeas = [clusterDets.Measurement];meas = mean(clusterMeas,2);meas2D = [meas(1:2);meas(4:5)];i = i+1;detectionClusters{i} = detections{detInds(1)};detectionClusters{i}.Measurement = meas2D;leftToCheck(clusterInds) = [];enddetectionClusters(i+1:end) = [];
卡尔曼滤波器（Kalman Filter）
% This function initializes a constant velocity filter based on a detection.functionfilter =initSimDemoFilter(detection)% Use a 2-D constant velocity model to initialize a trackingKF filter.% The state vector is [x;vx;y;vy]% The detection measurement vector is [x;y;vx;vy]% As a result, the measurement model is H = [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
%TODO: Implement the Kalman filter using trackingKF function. If stuck%review the implementation discussed in the project walkthroughH = [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1];filter = trackingKF('MotionModel', '2D Constant Velocity', ... 'State', H' * detection.Measurement, ... 'MeasurementModel', H, ... 'StateCovariance', H' * detection.MeasurementNoise * H, ... 'MeasurementNoise', detection.MeasurementNoise);end