炉火纯青：毫米波雷达开发手册之大话空间谱估计

写在前面

​ 深知新手在接触毫米波雷达板硬件时需要花费的沉没成本，因此在行将告别毫米波雷达之际，总结这两年以来在毫米波雷达上的一些经验和教训。

​ 本文档用于为实现基于AWR1243BOOST等单板毫米波雷达开发提供参考指南与解决方案，主要包括硬件配置、基础参数、信号模型、应用DEMO开发以及可深入研究方向思考等；为更好地匹配后续级联雷达应用的学习路线，在本手册中会尽可能同化单板雷达和级联雷达中的相关表述。

​ 本指南作者信息：Xuliang，联系方式：22134033@zju.edu.cn。未经本人允许，请勿用于商业和学术用途。

​ 希望后者在使用本指南时可以考虑引用作者在毫米波雷达旅途中的相关工作，如本文参考文献[1].
本章节为可深入研究方向思考章节之空间谱估计，主要解读子空间方法和压缩感知方法。
欢迎各位读者通过邮件形式与笔者交流讨论，本章节完整程序请私信笔者，希望使用本代码时能够提供一份引用和Star，以表示对笔者工作的尊重，谢谢！在后续将定时维护更新。

往期内容：
登堂入室：毫米波雷达开发手册之信号模型
眼观四海：自动驾驶&4D成像毫米波雷达 如今几何？
扬帆起航：毫米波雷达开发手册之硬件配置
初出茅庐：毫米波雷达开发手册之基础参数

空间谱估计算法

信号模型

空间谱估计是利用空间阵列实现空间信号的参数估计的技术，空间谱估计系统由空间信号入射、空间阵列接收以及参数估计等三部分组成，相应地可分为三个空间即目标空间、观察空间以及估计空间。

目标空间通常是一个由信号源的参数与复杂环境参数张成的空间。观察空间是利用空间按一定方式排列的阵元来接收目标空间的辐射信号，接收数据中往往包含信号特征（方位、距离、极化等）和空间环境特征（噪声、杂波、干扰等），观察空间是一个多维空间。估计空间是利用空间谱估计技术从复杂的观察数据中提取信号的特征参数。

阵列接收信号为(s(t))，目标空间源信号载波为(exp(jomega t))，信号在空间沿波束向量(k)的方向传播，基准点处信号为(s(t)exp(jomega t))，则距离基准点处的阵元接收信号为(s_{r}(t)=s(t-frac{1}{c}boldsymbol{r}^{T}/t)e^{j(omega t-r^{T}k)}).

对M阵列而言，基于窄带信号的假设，可以认为(frac{1}{c}boldsymbol{r}^{T}boldsymbol{a}llfrac{1}{B},)因此阵列信号以向量表示为(boldsymbol{s(t)}=s(t)[e^{-jr_1^Tk},e^{-jr_2^Tk},....,e^{-jr_M^Tk}])，,若设第一个阵元为基准点且初始相位为0，那么可以得到导向矢量（方向矢量，(M×1)维）为:(boldsymbol a(boldsymboltheta)=left[boldsymbol1,,boldsymbol e^{-jr_2^Tk},....,boldsymbol e^{-jr_M^Tk}right]^T=left[boldsymbol1,,boldsymbol e^{jr_2^Tk},....,boldsymbol e^{jr_M^Tk}right]^H.)

对M阵列-K信号源模型而言，M个具有全向性的阵元按任意排列构成，并设有K个具有相同中心频率(omega_0)、波长为(lambda)的空间窄带平面波分别以角(Phi_1,Phi_2,dotsPhi_K)入射，(Phi_{i}=(theta_{i},phi_{i}))，阵列第(m)个阵元的输出可以表示为：(x_m(t)=Sigma_{i=1}^K s_i(t)e^{jomega_0tau_m(Phi_i)} + n_m(t))

其中，(s_i(t))表示入射到阵列的第(i)个源信号，(n_m(t))表示第个阵元的加性噪声，(tau_m(Phi_i))为来自(Phi_i)方向的源信号投射至第(m)个阵元时相对选定参考点的时延，(theta_i)表示俯仰角，(phi_i)表示方位角。在AWR1243BOOST等毫米波雷达中需要分析方位图、点云图，方位角和俯仰角是实现这些分析的基础。

根据第(m)个阵元的输出，可以推广至空间的表示，记观察空间为(X(t)),噪声空间为(N(t))，流形为(A(Phi)), (S(t))为源信号。那么有如下：

(X(t)=[x_1(t),x_2(t),...,x_M(t)]^T)

(N(t)=[n_1(t),n_2(t),...,n_M(t)]^T)

(S(t)=[s_1(t),s_2(t),...,s_M(t)]^T)

(A(Phi)=[a(Phi_1),a(Phi_2),...,a(Phi_K)])

观察空间可用噪声空间、目标空间和流形表示，即满足下式：(boldsymbol{X}(boldsymbol{t})=boldsymbol{A}(boldsymbol{Phi})S(boldsymbol{t})+N(boldsymbol{t})overline{})

流形(boldsymbol{A}(boldsymbol{Phi}))与阵列的形状、信号源的来向有关，通常在实际应用中天线阵的形状一旦固定则不会发生改变，因此流形任意一列总是和目标空间源信号的来向有关。阵列形状通常有均匀线阵、均匀圆阵、L型线阵、平面阵列、任意阵列等。

名词解释

波数向量：(|boldsymbol{k}|=frac{omega}{c}=frac{2pi}{lambda})，空间距离变化1m时相位的变化量

信号相对于基准点的延时时间：(frac{1}{c}boldsymbol{r}^Tboldsymbol{a})

电磁波传播至离基准点(r)处的阵元相对于电磁波传播到基准点的滞后相位：(r^Tk)

方向矢量/导向矢量：阵列相对基准阵列的相位向量(a(Phi))

阵列流形：流形表示由K个信号源表示的方向向量组成的矩阵(A(Phi))

快时间维&慢时间维：脉冲雷达往往伴随快时间和慢时间维这两个概念，毫米波FMCW雷达本质上也是一种脉冲雷达。快时间维是针对单个脉冲而言，是对接收回波信号按行存储。快时间维的采样频率为采样率。慢时间维是针对多个脉冲而言，是对接收回波信号按列存储。每M个脉冲参与处理，脉冲间的时间间隔为脉冲重复间隔（Pulse Repetition Interval，PRI），脉冲重复频率（Pulse Repetition Frequency，PRF）是1/PRI。慢时间维的采样频率为PRF。暴力来说，在毫米波雷达里面快时间维可以理解为距离维，慢时间维可以理解为多普勒或速度维。

毫米波雷达与空间信号关系

AWR1243BOOST等雷达采集数据通常由DCA1000回传，回传数据文件需要经过解析得到雷达信号矩阵，该矩阵的维度为4D，分别为距离维ADC采样点数×虚拟天线通道数×Chirp Loop数目（单帧发送的Chirp数目）×帧数，但是我们这里研究的空间信号是2D，这要如何匹配呢？

在标准的毫米波雷达信号处理流程中，4D tensor数据实际在处理的时候需要逐帧处理，每帧的3d tensor数据经过快时间维FFT和慢时间维FFT（或称2D-FFT）得到距离多普勒谱图，再在距离多普勒谱图上做恒虚警率检测得到对应距离和速度的目标信号，此时我们得到的目标信号维度是关于天线维度的向量，这个向量可以理解为单快拍下的空间信号矩阵(X)；实际上，我们也可以通过仅对快时间FFT后的数据作恒虚警率检测得到特定距离的目标，这个时候我们得到的目标信号维度是和天线维度、Chirp数目相关的矩阵，那这个矩阵可以理解为多快拍下的空间信号矩阵(X)，这里谈到的(X)是和信号模型中提到的(X)等价的。

无论是相对运动还是静止，固定帧和距离单元后，由目标空间入射的源信号方向可以认为不发生变化或保持平稳随机，统计特性不随时间变化，因此定义阵列的协方差矩阵(R)为：

(boldsymbol{R}=boldsymbol{E}(bigl(X(boldsymbol{t})-boldsymbol{m}_x(boldsymbol{t})bigr)bigl(X(boldsymbol{t})-boldsymbol{m}_x(boldsymbol{t})bigr)^H)

({boldsymbol m}_{{boldsymbol x}}({boldsymbol t})={boldsymbol E}[{boldsymbol X}({boldsymbol t})],{boldsymbol m}_{{boldsymbol x}}({boldsymbol t})={boldsymbol0})

(textbf{R}=E{boldsymbol{X}(boldsymbol{t})boldsymbol{X}(boldsymbol{t})^{boldsymbol{H}}}=Eleft{big(boldsymbol{A}(boldsymbol{Phi})boldsymbol{S}(boldsymbol{t})+boldsymbol{N}(boldsymbol{t})big)big(boldsymbol{A}(boldsymbol{Phi})boldsymbol{S}(boldsymbol{t})+boldsymbol{N}(boldsymbol{t})big)^{boldsymbol{H}}right}=A(Phi)S(boldsymbol{t})boldsymbol{S}(boldsymbol{t})^H A(boldsymbol{Phi})^H+boldsymbol{sigma}^2I=A(boldsymbol{Phi})boldsymbol{R}_s A(boldsymbol{Phi})^H+boldsymbol{R}_N)

阵列信号(X(t))的协方差可以表示为目标子空间与噪声子空间的加形式，也可以通过特征分解表示为(M)个特征值与特征矢量的积和形式。

(boldsymbolSigma=mathbf{diag}(boldsymbollambda_1,boldsymbollambda_2,...,boldsymbollambda_M),)通过特征值排序，(boldsymbol{lambda}_{1}geqboldsymbol{lambda}_{2}geqcdotsgeqboldsymbol{lambda}_{K}succboldsymbol{lambda}_{K+1}=cdots=boldsymbol{lambda}_{M}=boldsymbol{sigma}^{2},)前(K)个个特征值可以认为是与目标子空间信号相关的，由其对应的特征向量可以表示目标信号子空间(U_s)，后(M-K)个特征值完全取决于噪声，数值等于(sigma^2)，，由其对应的特征向量构成噪声子空间(U_N)。因此，阵列信号(X(t))的协方差可以进一步表示为：(boldsymbol{R}=boldsymbol{USigma}U^H=boldsymbol{U_S}boldsymbol{Sigma_S}U_S^H+boldsymbol{U_N}boldsymbol{Sigma_N}U_N^H)

子空间方法

Vanilla-MUSIC

这种将阵列信号的协方差矩阵进行特征分解，得到与信号分量对应的信号子空间和信号分量正交的噪声子空间，利用两个子空间的正交性来估计信号参数。根据信号子空间与噪声子空间的正交性，可以得到以下表述：

(boldsymbol{R}U_N=[U_S,U_N]boldsymbol{Sigma}begin{bmatrix}U_S^H\U_N^Hend{bmatrix}boldsymbol{U}_N=[U_S,U_N]boldsymbol{Sigma}begin{bmatrix}boldsymbol{O}\Iend{bmatrix}=[U_S,U_N]begin{bmatrix}boldsymbol{sigma}_S^2&cdot\.&boldsymbol{sigma}_N^2end{bmatrix}begin{bmatrix}boldsymbol{O}\Iend{bmatrix}=sigma_N^2U_N)

(A(Phi)R_{s}A(Phi)^{H}U_{N}=0)

(U_N^H A(Phi)R_s A(Phi)^H U_N=(A(Phi)^H U_N)^H R_s A(Phi)^H U_N=0)

因为(boldsymbol{Aneq0},boldsymbol{R_{s}neq0},)故(boldsymbol{A(Phi)}^{H}boldsymbol{U}_{N}=boldsymbol{0})，这表明在无噪声情况下的信号(x(t)=Sigma_{i=1}^M s_i(t)a_i(theta)to x(t)in span{a_1,a_2,...,a_K}.)协方差矩阵大特征值对应的特征向量张成的空间与入射信号的导向矢量张成的空间是同一个空间。即(begin{matrix}U_{S}=[e_{1},e_{2},...,e_{K}],U_{N}=[e_{K+1},e_{K+2},...,e_{N}]\ textit{span}{e_{1},e_{2},...,e_{K}}=textit{span}{a_{1},a_{2},...,a_{K}}\ end{matrix})

因信号子空间和噪声子空间相互正交，可知信号子空间的导向矢量(A(Phi))与噪声子空间正交，但实际由于干扰(boldsymbol{a}^{H}(boldsymbol{Phi})U_{N}=boldsymbol{0})不完全成立，即不完全满足正交性。故采用最小优化搜索（零谱搜索）来估计波达方向：(Phi_{text{music}}=boldsymbol{arg_Phi min}a^H(boldsymbolPhi)U_N U_N^H a(boldsymbolPhi))

根据信号子空间的导向矢量与噪声子空间正交的原理，信号入射方向上会出现极小值，因此空间谱函数可以表示为：(mathbf{P}_{text{music}}(theta)=frac{1}{boldsymbol{a}^H(Phi)U_N U_N^Hboldsymbol{a}(Phi)})

那么通过遍历(theta)，当(mathbf{P}_{text{music}}(theta))由极大值时，对应角为估计的角度。

下面给出Vanilla-MUSIC的单独求解方位角和联合求解角案例。

function [PoutMusic] = DOA_MUSIC(X, P, searchGrids) 	% By Xuliang     % X: 输入信号 Channel * ChirpNum     % P: 目标数目     % PoutMusic: 输出功率谱          M = size(X, 1); % 阵元数     snap = size(X, 2); % 快拍数     RX = X * X' / snap; % 协方差矩阵          [V, D] = eig(RX); % 特征值分解     eig_value = real(diag(D)); % 提取特征值     [B, I] = sort(eig_value, 'descend'); % 排序特征值     EN = V(:, I(P+1:end)); % 提取噪声子空间          PoutMusic = zeros(1, length(searchGrids));          for id = 1 : length(searchGrids)         atheta_vec = exp(-1j * 2 * pi * [0:M-1]' * 1 / 2 * sind(searchGrids(id))); % 导向矢量         PoutMusic(id) = (abs(1 / (atheta_vec' * EN * EN' * atheta_vec))) ; % 功率谱计算     end end  function [index1,index2,Pmusic] = MUSIC_2D(X,P) %% MUSIC算法：用于方位角和俯仰角的联合估计 % By Xuliang % X: 输入信号 Channel * ChirpNum % P: 目标数目  	M = size(X, 1); % 阵元数     snap = size(X, 2); % 快拍数          d2rad = pi / 180; % pi rad = 180 ° 1°= pi/180 rad     lambda = physconst('lightspeed') / 77e9; % 波长          D = 0.5 * lambda; % 阵元间距     D1 = 0:D:(M-1)*D; % TX0和TX2的线阵间距     D2 = 2 * D : D : 5 * D; % TX1的线阵间距     D_AZI = [D1,D2]; % TX0 TX2 和 TX1的水平间距     D_ELE = [zeros(1,8) D*ones(1,4)]; % TX0/TX2和TX1的纵向间距      %% 计算协方差矩阵     RX = X1 * X1' / snap;     [EV,D] = eig(RX); % 特征值分解     EVA = diag(D);% 特征值对角线提取并转为一行     [~,I] = sort(EVA);% 特征值排序从小到大，eig默认有排序     EV = fliplr(EV(:,I));% 特征矢量排序     EN = EV(:,P+1:M); % 噪声子空间      %% 遍历每个方位角和俯仰角 计算空间谱     for ele = 1:181 % 俯仰角遍历         phim(ele) = ele - 91;          phim1 = d2rad * phim(ele);         for azi = 1:181 % 方位角遍历             theta(azi) = azi - 91;             theta1 = d2rad * theta(azi);             a = exp(-1j * 2 * pi * (D_AZI * cos(phim1) * sin(theta1) + D_ELE * sin(phim1)) / lambda).'; % 构建导向矢量             Pmusic(azi,ele) = 1 / (a' * EN *EN' * a); % 空间谱函数         end     end     Pmusic = abs(Pmusic);     [index1,index2] = find(Pmusic == max(max(Pmusic))); % 找到空间谱谱峰并返回俯仰角和方位角 end 

% 如何对特征值可视化判断信源数？ [V, D] = eig(RX); % 特征值分解 SP = V(:, M-P+1:M); EN = V(:, 1:M-P);  figure(1); % 特征值分布的可视化 plot(diag(D),'kd-'); xlabel('Number of Eigenvalues'); ylabel('Eigenvalues');  % 如何验证信号和噪声空间的正交性？ s_eigen_fft = fft(SP); n_eigen_fft = fft(EN); plot(abs(s_eigen_fft(:,1:2)), 'ks-','LineWidth',1.5); hold on; plot(abs(n_eigen_fft(:,1:2)), 'rd-','LineWidth',1.5); % legend('Signal Space', 'Noise Space'); 

ESPRIT

Waiting ..

DML

压缩感知

网格模型

无网格模型

参考文献

[1] X. Yu, Z. Cao, Z. Wu, C. Song, J. Zhu and Z. Xu, "A Novel Potential Drowning Detection System Based on Millimeter-Wave Radar," 2022 17th International Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, Singapore, 2022, pp. 659-664, doi: 10.1109/ICARCV57592.2022.10004245.