读书笔记 | 自动驾驶中的雷达信号处理（第6章 到达方向(DOA)估计算法 ）

本文编辑：调皮哥的小助理

6.1介绍

DOA 估计算法在汽车应用雷达处理中非常重要，它构成了雷达数据立方体（距离、速度和角度）的第三个部分——角度。实际上，会有多个未知数量的源信号同时被接收阵列接收，且幅度未知，因此DOA 估计算法十分复杂。此外，接收到的源信号几乎总是被加性噪声破坏并且存在杂波。

除了这些挑战之外，我们还存在多径问题需要处理。尽管 DOA 估计算法并不简单，但经过50多年的发展，有许多算法可用于估计源信号的数量及其方向。

本文目的：提供一些可用且持续被人们研究的 DOA 估计算法，这些算法是汽车雷达处理算法的关键。但需要注意的是，本文并没有详细分析每种算法的原理，只是简要论述了最流行算法的优势和局限性。若要对DOA估计算法有深入的研究，我推荐各位读者阅读下面这两本书：

两本书都在【雷达开源资料库】中：

或者下面这几本书：

《现代数字信号处理》——王展，李双勋等
《最优阵列处理技术》——Harry L. Van Trees…
《The Difference Between Capon and MUSIC Algorithms》 Huiping Huang

6.2 DOA估计算法分类

DOA 方法可以大致分为二次预测（传统法）、线性预测和子空间法 ，如数字波束形成可分为二次法、前向后线性预测法和MUSIC子空间法。

6.3 DOA估计算法

已被广泛使用的几种 DOA 估计算法如下表中所示，但目前还在进行改进研究，一直在不断地对这些方法在各种期刊、会议论文集等中从性能和应用的角度进行比较。 即使在技术进步的情况下，似乎也没有什么神奇的方法可以解决角度分辨率、算法复杂性和性能鲁棒性等问题。 在汽车应用中，由于天线尺寸和位置的限制，我们对这些算法的要求变得更加严格。

6.3.1信号模型

每种算法都需要定义一个信号模型，故考虑一个由 $M$个传感器（天线阵元）阵列组成的雷达系统，接收来自$K$个源（目标）的信号，则接收到的信号可以表示为：

$X(t)=A(\theta )s(t)+N(t)$（6.1）

其中 $X(t)={\left[x_1(t),\dots ,x_M(t)\right]}^T$ 是 $M\times 1$ 接收到的传感器数据向量，$A(\theta )=\left[a\left({\theta }_1\right),\dots ,a\left({\theta }_K\right)\right]$ 是 $M\times K$流形矩阵（也称为方向矩阵）$s(t)={\left[s_1(t),\dots ,s_K(t)\right]}^T$ 是 $K\times 1$源信号向量， $N(t)={\left[n_1(t),\dots ,n_M(t)\right]}^T$ 表示方差 ${\sigma }^2$ 的 $M\times 1$ 传感器噪声向量。

方向矩阵的元素由定义的转向向量组成：

$\mathbf{a}\left({\theta }_i\right)={\left[1,{\mathrm{e}}^{-\frac{j2\pi d\sin \left({\theta }_i\right)}{\lambda }},\dots ,{\mathrm{e}}^{-\frac{j2\pi d(N-1)\sin \left({\theta }_i\right)}{\lambda }}\right]}^{\mathrm{T}}$（6.2）

其中， $d$ 是天线阵元间距，采用均匀线阵（ ULA），$\lambda$ 是传播信号的波长，${\theta }_i$ 是来自第 ${\theta }_i$ 个源的信号的到达方向（角度）， $"T"$ 表示转置操作， 信号模型如图 1 所示。

（图1 信号模型的示意图。到达方向θ对应于任何源方向θi d是天线单元间距，天线阵元的最小数目为2个。）

在一些需要优化权重的算法中，如 DBF 和 Capon，需要考虑传感器的加权输出， 在这种情况下，传感器输出的加权线性组合可以表示如下：

$y(t)={\sum }_{m=1}^M{w}_i^{\ast }{x}_i(t)=w^{H}X$（6.3）

其中，$w_i^{\ast }$ 是第$i$个传感器的权重， 星号 (*) 表示复共轭，而 $H$ 表示共轭或 Hermitian 转置，则传感器阵列的输出功率 $P(w)$ 可以表示如下：

$P(\mathbf{w})=E\left[|y(t){|}^2\right]=w^{H}E\left[X{X}^H\right]w=w^H\mathbf{R}w$(6.4)

其中 $E[\cdot ]$ 表示期望，其中 $R$ 表示输入信号的协方差矩阵。

6.3.2 DOA估计方法

在下面的几节中，我们将详细介绍各种常见DOA算法的原理。

（1）数字波束形成(DBF)

DBF算法通过优化权值向量，使特定方向的输出功率最大化，最优权向量由下式给出：

${\mathbf{w}}_{\mathrm{opt}}=\frac{\mathbf{a}(\theta )}{\sqrt{{\mathbf{a}}^H(\theta )\mathbf{a}(\theta )}}(6.5)$（6.5）

将最优权向量代入式(6.4)，得到DBF功率谱如下：

$P_{\mathrm{DBF}}(\theta )=\frac{a^H(\theta )\mathbf{Ra}(\theta )}{a^H(\theta )a(\theta }$（6.6）

性质

DBF只需要计算接收数据向量的协方差矩阵 $R$ ，计算复杂度较低。通过在期望的间隔内建立一个等间距角度的转向向量，谱搜索后则可以得到与最大功率相对应的角度作为DOA。

但DBF无法分辨密集目标，不过它可以作为第一步来缩小高分辨率方法的目标范围，如先进行初步角度搜索，然后再进行精细搜索。**

（2）Capon算法

Capon 算法旨在为从特定方向到达的信号保持恒定增益，同时赋予噪声较小的权重 。 它是最大似然法的替代算法，用于求解阵列的最小方差无失真响应 (MVDR)，从而使信干比最大化。 优化问题可以表示为：

$\min (P(\mathbf{w}))$ subject to $w^{H}a(\theta )=1$(6.7)

得到的最佳权重向量由下式给出：
$w_{\mathrm{opt}}=\frac{{\mathbf{R}}^{-1}a(\theta )}{a^H(\theta ){\mathbf{R}}^{-1}a(\theta )}$(6.8)

将最优权重向量代入方程。 (6.4) 给出 Capon 功率谱如下：

$P_{\text{Capon }}(\theta )=\frac{1}{a^H(\theta )R^{-1}a(\theta )}$（6.9）

性质

Capon 方法在分辨紧密相邻的目标方面为 DBF 提供了卓越的性能 ，此外，主要计算是从数据中确定逆协方差矩阵，这使得它在许多需要高分辨率的情况下非常有吸引力。还值得注意的是，Capon 和 DBF 都不需要事先了解信号源的数量，可以同时估计源信号的 DOA 和功率。

Capon 分离多个目标的能力仅受 SNR 和接收器阵列大小的限制，这是大多数高分辨率 DOA 方法的常见限制。

Capon 方法的缺点在于需要稳定计算，确定近奇异协方差矩阵的逆矩阵并非易事，尤其是在接收信号因噪声和杂波而劣化的汽车应用中。由于天线的孔径有限，如果没有某种形式的方法作为预处理的一部分，就无法可靠地计算协方差矩阵。然而，计算协方差矩阵的这个障碍并没有阻止 Capon 的使用。当然，如上所述，必须应用一些预处理技术。

（3）Multiple Signal Classifier (MUSIC)

MUSIC 采用子空间法来解决 DOA 估计问题 ，它利用特征分解作为主要工具。 特征分解背后的理论可以在文献[8]中找到，在这里，我们给出要点。 子空间方法背后的主要假设是信号和噪声子空间是正交的，则使用正交性原理，可以计算伪谱， 协方差矩阵可以分为信号和噪声分量。

根据定义，$R=E\left[X{X}^H\right]$ 밈 可以写为:

$\begin{array}{rl}\mathbf{R}& =E\left[(\mathbf{A}(\theta )\mathbf{s}(t)+\mathbf{N}(t))(\mathbf{A}(\theta )\mathbf{s}(t)+\mathbf{N}(t){)}^H\right)]\\ & =\mathbf{A}(\theta )\mathbf{E}\left[s(t){\mathbf{s}}^H(t)\right]\mathbf{A}(\theta {)}^H+E[N(t)(\mathbf{N}(t){)}^H\\ & =\mathbf{A}(\theta ){\mathbf{R}}_s\mathbf{A}(\theta {)}^H+{\sigma }^2\mathbf{I}\end{array}$(6.10)

其中 $R_s$ 是未观察到的源信号协方差矩阵，$I$ 表示单位矩阵。 为简单起见，等式 (6.10)中对时间 $t$已被省略。 可以观察到，协方差矩阵可以分为与 $\theta$ 无关的信号子空间和噪声子空间。 通过特征分解，协方差矩阵 $R$ 可以分解为特征向量和特征值矩阵。

由于信号子空间的大小为 $K$ ，因此 $R$ 的 $M-K$个特征值属于噪声子空间， 如果$u_i$是 $R$ 的特征向量，那么我们可以写成：

$\mathbf{R}{\mathbf{u}}_i=\left(\mathbf{A}(\theta ){\mathbf{R}}_s\mathbf{A}(\theta {)}^H+{\sigma }^2\mathbf{I}\right){\mathbf{u}}_i={\lambda }_i{\mathbf{u}}_i$（6.11）

其中${\lambda }_i,i=1,\dots ,M$ 是 $R$ 的第 $i$ 个特征值。 $R$ 的特征值可以划分为 ${\lambda }_i={\sigma }_i^2+{\sigma }^2$ , $i=1,\dots ,K$ 和 ${\lambda }_i={\sigma }^2,i=K+1,\dots ,M$ ，由于 $R_s$ 是正定的，我们可以写成：

${\mathbf{u}}_i^H\left(\mathbf{A}(\theta ){\mathbf{R}}_s\mathbf{A}(\theta {)}^H\right){\mathbf{u}}_i=0,i=K+1,\dots ,M$（6.12）
$\left(A(\theta {)}^H\right){\mathbf{u}}_i=0,i=K+1,\dots ,M$（6.13）
$\left(a{\left({\theta }_k\right)}^H\right)u_i=0,k=1,\dots ,K,i=K+1,\dots ,M$（6.14）

根据公式(6.14)，$MUSIC$ 算法找到信号子空间与噪声子空间正交的角度 ${\theta }_k$ ，噪声子空间的特征向量定义为：

${\mathbf{U}}_N=\left[{\mathbf{u}}_{K+1},\dots ,{\mathbf{u}}_M\right]$（6.15）

基于公式(6.15)，可以使用以下表达式计算 MUSIC 伪谱：

$P_{\text{MUSIC }}(\theta )=\frac{1}{a^H(\theta ){\mathbf{U}}_N{\mathbf{U}}_N^{H}a(\theta )}$(6.16)

性质

MUSIC算法是DOA估计的高分辨率技术之一。 它能够可靠地从严重劣化的接收信号中分辨出角度间隔低至 1° 的目标，伪谱中的杂散峰值非常少。 由于 MUSIC 计算接收信号的伪频谱，它不能用于估计信号源的功率。 对于功率估计，DBF 或 Capon 可以与 MUSIC 结合使用。 通过增加用于频谱估计的快拍数量，可以提高 MUSIC 分离多个目标的能力，调皮哥提醒各位读者，这句话要深刻理解。
尽管具有高分辨能力，但 MUSIC 需要事先了解信号子空间的大小，另外还取决于需要高计算复杂度的特征分解。

（4）Root-MUSIC

这是一种 $MUSIC$ 算法的变体，其中 DOA 估计通过找到多项式的根来执行 ，等式(6.16) 中 $P_{\text{MUSIC }}(\theta )$的分母可以表示为多项式$P(z)$，使得：

$P(z)=a^H(\theta ){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{a}(\theta )$（6.17）

其中：

$\mathbf{a}(\theta )={\left[\begin{array}{lllll}1& z^1& z^2& \dots & z^{M-1}\end{array}\right]}^{\mathrm{T}}$(6.18)

且有：

$z={\mathrm{e}}^{-\frac{j2\pi d}{\lambda }\sin (\theta )}$(6.19)

因此，我们可以将 $P(z)$ 展开成这种形式的多项式：

$P(z)={\sum }_{k=-M+1}^{M-1}{C}_k{z}^k$（6.20）

在理想条件下，$P(z)$ 的根将在单位圆上，其角度对应于接收信号的 DOA。 在实际情况下，噪声的存在会导致幅度从单位圆发生偏移。 因此，最接近单位圆的 $P(z)$ 的根可以被视为 MUSIC 伪谱的极点， 则根可以写成：

$z_k=\left|z_k\right|{\mathrm{e}}^{\arg \left(z_k\right)}$（6.21）

从幅度最接近 1 的根，DOA 可以估计为：

${\theta }_k=\sin \left(\frac{\lambda }{2\pi d}\arg \left(z_k\right)\right),k=1,\dots ,K$（6.22）

性质

Root-MUSIC 具有与 MUSIC 相似的性能。 在某些需要考虑相位噪声的情况下，Root-MUSIC可以提供更好的性能。 与其他子空间算法一样，仍然需要提前知道目标源的数量。 在噪声和杂波严重的情况下，根的大小可能会远离单位圆，从而难以得到取感兴趣的根。

（5）Estimation of SignalParameters viaRotational InvarianceTechnique (ESPRIT)

ESPRIT 是一种子空间方法，它采用与 MUSIC 不同的方法，因为它将阵列分解为两个子阵列，由一个已知向量置换，无需计算昂贵的峰值搜索过程就可以从中估计 DOA 。

图 2 显示了将单个阵列分解为两个子阵列的示意图。 应该注意，物理上存在一个单独的 ULA 阵列，从该阵列创建子阵列。 ESPRIT 还通过放宽天线几何结构减轻了校准任务的负担，因为子阵列不必具有相同的天线元件，并且它们的物理位置可以任意选择 。

（图2 （a）原始阵列）

（图2 （b）两个子阵列）

对于由 $x$ 和 $y$ 表示的大小为 $M$ 的两个子阵列，接收到的数据向量可以由以下等式表示：

$\mathbf{x}(t)={\sum }_{i=1}^{K}a\left({\theta }_i\right)s_i(t)+n_x(t)=As(t)+n_x(t)$
$\mathbf{y}(t)={\sum }_{i=1}^{K}a\left({\theta }_i\right){\mathrm{e}}^{j{\gamma }_i}{s}_i(t)+n_x(t)=A\Phi s(t)+n_y(t)$（6.23）

其中，${\gamma }_i=\frac{{\omega }_0\Delta }{c}\sin {\theta }_i,\Phi =\mathrm{diag}\left({\mathrm{e}}^{j{\gamma }_1},\dots ,{\mathrm{e}}^{j{\gamma }_K}\right)$ ，c 是光速，${\omega }_0$ 是中心频率， $\Delta$ 是两个阵列之间的平移位移。 $\mathrm{diag}(\cdot )$运算符将序列转换为对角矩阵，流形矩阵 $A$ 的大小为 $M\times K$ 。ESPRIT 旨在从旋转算子估计DOA， 这可以通过首先定义一个 $2M\times 1$ 向量 $z(t)$ 来实现，使得：

$\mathbf{z}(t)=\left[\begin{array}{l}\mathbf{x}(t)\\ \mathbf{y}(t)\end{array}\right]=\left[\begin{array}{c}A\\ \mathbf{A}\Phi \end{array}\right]s(t)+\left[\begin{array}{l}{n}_x(t)\\ n_y(t)\end{array}\right]=A_{z}s(t)+n_z(t)$（6.24）

$z(t)$ 的协方差矩阵由下式给出：

${\mathbf{R}}_z=E\left[\mathbf{z}(t)\mathbf{z}(t{)}^H\right]=A_z{\mathbf{R}}_s{A}_z^H+{\sigma }^2\mathbf{I}$（6.25）

的特征分解使得特征向量具有大小为 $K$ 的信号子空间和大小为 $2N-K$ 的噪声子空间。假设 $K\le M$ ，并将 $2N\times K$ 信号子空间特征向量表示为 $E_s$，从广义特征分解可知$E_s$ 和 $A_z$具有相同的跨度，存在一个非奇异变换矩阵 $T$ 使得：

$E_s=A_{z}T$（6.26）

此外， $E_s$ 可以划分为 $K$ 秩的 $E_x$ 和 $E_y$ ，如下所示：

$E_s=A_{z}T=\left[\begin{array}{c}A\\ \mathbf{A}\Phi \end{array}\right]T=\left[\begin{array}{c}AT\\ \mathbf{A}\Phi \mathrm{T}\end{array}\right]=\left[\begin{array}{l}{E}_x\\ E_y\end{array}\right]$（6.27）

定义一个大小为 $N\times 2K$ 且秩为 $K$ 的 $Exy=\left[E_x{E}_y\right]$ 矩阵，存在一个 $2K\times K$ 矩阵 $F$ 使得:

$0=E_{xy}F=\left[E_x{E}_y\right]\left[\begin{array}{l}{\mathbf{F}}_x\\ {\mathbf{F}}_y\end{array}\right]$
${\mathbf{E}}_x{\mathbf{F}}_x+{\mathbf{E}}_y{\mathbf{F}}_y=AT{\mathbf{F}}_x+\mathbf{A}\mathbf{\Phi }\mathbf{T}{\mathbf{F}}_y=0$
$\mathbf{A}\Phi \mathrm{T}=-AT{\mathbf{F}}_x{\mathbf{F}}_y^{-1}$（6.28）

上面的表达式意味着$F$跨越 $\left[E_x{E}_y\right]$ 的零空间，因为我们正在假设$K\le M$ 。由于 $T$是满列秩，它的逆存在，则有：
$$\begin{array}{rl}& \mathrm{A}\Phi =-AT{F}_x{F}_y^{-1}{\mathrm{T}}^{-1}\\ & \Phi =T{F}_x{F}_y^{-1}{\mathrm{T}}^{-1}=T\psi {\mathrm{T}}^{-1}\end{array}$$（6.29）

其中， $\psi =-F_x{F}_y^{-1}$ 。

$\psi$ 的特征值等于 $\Phi$ 的对角元素，由此可以估计DOA。由于测量中存在噪声和校准误差，通常使用TLS算法计算$\Phi$ 。

（6）TLS ESPRIT 算法

TLS算法可以概括为以下步骤：

1. 根据测量值计算 $R_z$ 的估计值。

2. 对 $R_z$ 进行广义特征分解，得到 $R_z{E}_z={\Sigma }_z{E}_z{\Lambda }_z$ 。

3）估计信号子空间大小K。

4. 估计 $E_s$ 并分解为$E_x$ 和 $E_y$。

5. 定义矩阵 $E_D$ 并进行特征分解如下：

$E_D=\left[\begin{array}{c}{E}_x^H\\ E_y^H\end{array}\right]\left[\begin{array}{ll}{\mathbf{E}}_x& {\mathbf{E}}_y\end{array}\right]=\mathbf{E}\mathbf{\Lambda }{\mathbf{E}}^H$（6.30）

6. 将$E$ 分解为 $4K\times K$矩阵，如下所示：

$E=\left[\begin{array}{ll}{\mathbf{E}}_{11}& {\mathbf{E}}_{12}\\ {\mathbf{E}}_{21}& {\mathbf{E}}_{22}\end{array}\right]$（6.31）

7. 计算 $\psi =-E_{12}{E}_{22}^{-1}$的特征值得到:

${\hat{\varphi }}_k={\lambda }_k\left(-{\mathbf{E}}_{12}{\mathbf{E}}_{22}\right),k=1,\dots ,K$（6.32）

8. DOA 估计是通过下面等式计算得到：

${\hat{\theta }}_k={\sin }^{-1}\left(c\frac{\arg \left({\hat{\varphi }}_k\right)}{{\omega }_0\Delta }\right)$（6.33）

其中，特征值是从 $\psi$ 中提取的，如上所述。文献中存在 ESPRIT 算法的几种变体和改进，例如单一 ESPRIT 和 波束空间 ESPRIT，其详细信息可在文献 [1] 中找到。

性质

如上所述，ESPRIT 的优点是不需要穷举峰值搜索来确定 DOA ，因此，它在计算复杂度方面更具优势。 然而，额外传感器的需求可能会导致噪声 DOA 估计。 与所有算法一样，算法性能的关键在于准确计算协方差矩阵的能力，不过，随着更多传感器的可用性，协方差估计可以得到改进。

6.3.3空间平滑

从上述算法可以看出，协方差矩阵 $R$ 的估计是 DOA 估计过程的关键。 校正矩阵的大小由接收数据向量的大小决定，而向量大小又取决于接收天线阵列的长度。

例如，在汽车雷达应用中，由于对雷达传感器尺寸的限制，通常使用小于10个阵元的接收天线阵列， 通常，天线元件的数量可以少到4个。 这些阵元接收的数据是协方差矩阵估计方法的输入（这是重点内容）， 在大多数情况下，为了降低计算复杂度，接收信号通过离散傅里叶变换变换到频域，几乎总是采用FFT实现。

对于基于特征的算法，DOA 估计只有在协方差矩阵是非奇异的情况下才有可能，这意味着信号是非相干或不相关的。 相干性的来源是由于多径等自然传播特性，也可能是干扰等人工信号的结果。 用于去相关信号的一种方法是空间平滑（SS）。

空间平滑是通过首先将天线阵列从可以计算的平均协方差矩阵划分为子阵列来实现的。 假设数组可以划分为 $M$ 个子数组，并将每个子数组计算的协方差矩阵表示为 $R_m^f,m=1,\dots ,M$ ，上标表示正向计算，即增加元素 如图 3 所示，我们可以使用以下表达式计算协方差矩阵。

$R_{\mathrm{SS}}^f=\frac{1}{M}{\sum }_{m=1}^M{R}_m^f$（6.34）

（图3 空间平滑处理子阵列构造示意图）

空间平滑协方差矩阵 $R_{SS}^f$ 可用于计算 DOA 估计的特征值和特征向量。 如果子阵列中的元素个数是 $p$ ，那么子阵列的数量 $M=N-p+1$ 。为了检测 $K$个源， $p$ 必须大于等于 $K+1$ ，则所需的阵元将是 $2K$ 。

因此，空间平滑有效地将可检测目标的数量减少了一半，对于最多 N-1 个的 N 阵元阵列，与相干情况下的 N/2 相比，在非相干情况下无需空间平滑即可检测到目标。 这意味着对于固定数量的阵列阵元，必须在子阵列数量和可检测目标之间进行权衡。

作为对空间平滑目标检测问题的改进，在文献[12]中提出了前向后空间平滑（FBSS）。 顾名思义，FBSS 方法对前向计算的协方差矩阵进行平均，如 SS 方法中所述，以及从第 N 个阵元计算回第一个阵元的协方差矩阵。 使用交换矩阵 $J$ 可以有效地完成反向计算。使用以下表达式计算整体协方差矩阵 RFBSS：

$R_{\mathrm{FBSS}}=\frac{{\mathbf{R}}_{\mathrm{SS}}^f+\mathbf{J}{\left({\mathbf{R}}_{\mathrm{SS}}^f\right)}^H{}^H}{2}$（6.35）

其中：

$\mathbf{J}=\left[\begin{array}{ccc}0& \cdots & 1\\ ⋮& \ddots & ⋮\\ 1& \cdots & 0\end{array}\right]$（6.36）

对于 N 个阵元的阵列，可检测的相干源数量增加到 2N/3。

作为DOA估计的例子，我们考虑了三种广泛使用的方法：DBF、Capon和MUSIC。

光速，c：3.0e8 [m/s]
中心频率，f0：77.5e9 [Hz]
ULA 天线单元数，N：5
天线间距，d：0.5 * λ，（λ为波长）
目标源数量，N：2
目标源 1， 的到达方向：-10 [deg]，-2 [deg]
目标源 2， 的到达方向：+10 [deg]、+2 [deg]
噪声方差：0.015。

使用四阵元长度的子阵列将空间平滑应用于协方差矩阵。 这总共给出了两个子阵列用于平均，仿真结果如图 4、5 和 6 所示。

(图 4 DBF方法在−10[deg]和+10[deg]两个源下使用五元ULA天线估计DOA)

(图5 使用五元ULA天线，在−10度和+10度两个源处使用Capon方法估计DOA)

(图6 使用五元ULA天线，在−10和+10两个源处使用MUSIC方法估计DOA)

结果表明，与更传统的 DBF 和 Capon 方法相比，子空间方法 MUSIC 提供了更清晰的频谱选择。 如图 7、8 和 9 所示，当源在−2 [deg] 和 +2 [deg] 处更靠近时，MUSIC 比 DBF 和 Capon 更好地解析它们。

(图7 DBF方法在−2[度]和+2[度]处使用五元ULA天线进行DOA估计)

(图8 使用五元ULA天线对源2在−2[度]和2[度]处的Capon算法估计DOA)

(图9 使用五元ULA天线，在−2[度]和2[度]两个源下，MUSIC方法估计DOA)

关于几种DOA估计方法的效果，【雷达技术交流群】里有大佬做出过总结：

通过瑞利准则算出来角度分辨率是4，那么FFT算法的分辨率就是4，MUSIC算法（子空间类）可能是2，压缩感知算法可能是1。如果再增大孔径，瑞利准则算出来2，FFT算法的角度分辨率为2，MUSIC算法（子空间类）可能是1，压缩感知算法0.5，以此类推。
然后再结合上面的三个效果图，自然就能够分析清楚各种算法的在角度分辨率上的优劣，但仅仅是在角度分辨率上做比较，其他方面则需要重新考虑，各种算法计算量不同。

6.3.4其他DOA算法

除上述算法外，还有一些算法如ML、LP和WSF也存在[1]。然而，其稳定性和计算复杂度的问题使其难以在汽车应用中应用。这并不一定意味着这些算法较差，但它们受到当前可用技术的限制，无法用于实时情况。最近的一些方法，如传播子方法，完全避免了特征分解，而采用最小二乘法，这使它们值得探索。下面我们将简要地解释其中的一些算法。

（1）Minimum Norm Method（最小范数法）

最小范数法是最古老的 DOA 估计高分辨率方法之一 [15]。 从估计的协方差矩阵 R 中，执行奇异值分解 (SVD) 以获得 U、S 和 V矩阵。噪声子空间特征向量被提取为 $E_N=U(:,K+1:N)$，这意味着所有列 从第 (K + 1) 列到第 N 列的 U矩阵。 频谱是基于位于第一个元素等于 1 的噪声子空间中的最小范数向量构造的。频谱由以下表达式计算：

$P_{\mathrm{MinNorm}}(\theta )=\frac{1}{{|a(\theta )E_N{E}_N^{H}u|}^2}$(6.37)

其中$u=[1,0,0,\dots 0{]}^T$ 和 $a(\theta )$ 是前面定义的转向向量。 由于最小范数方法属于基于特征的方法，它不能用于可靠地估计源信号的功率。 它主要用于DOA的估计。 必须事先知道有关来源数量的信息。 在退化环境中估计 DOA 是可能的，但在频谱中获得杂散峰的趋势增加。

（2）Maximum Entropy Method (MEM)最大熵法

MEM 是基于自相关函数外推的频谱估计方法，旨在通过使用自回归 (AR) 系数 [16] 来最大化信号熵（“不确定性”）。 AR 系数 a 使预测误差最小化，即受制于 $a^H{e}_1=1$ 的约束，其中 $e_1=[1,0,\dots 0{]}^T$ 。 估计 AR 参数的其他抗噪方法的详细信息可以在 [17] 中找到。 使用拉格朗日方法，AR 系数可以计算为

$\mathbf{a}=\mathrm{argmin}\left\{a^{H}R\right\}$(6.38)
$a=\frac{R^{-1}{e}_1}{e_1^T{R}^{-1}{e}_1}$(6.39)

频谱估计是使用:

$P_{\mathrm{MEM}}(\theta )=\frac{1}{{|a(\theta {)}^H{C}_j|}^2}$(6.40)

其中$C_j$表示协方差矩阵的逆矩阵的第 $j$ 列。j 的选择是任意的，会影响 MEM 的性能。AR 参数和谱估计之间的关系基于自相关函数存在且外推有效的。 这些假设的有效性存在一些问题，这些问题已在 [16] 中进行了概述。 MEM 可能很难扩展到多维 DOA 估计问题。

（3）Linear Prediction 线性预测

线性预测（LP）方法广泛用于音频和语音处理。 LP 背后的想法是最小化平均输出功率，前提是阵列中任意选择的元素的权重是统一的。 数组权重向量由下式给出：

$w=\frac{R^{-1}u}{u^H{R}^{-1}u}$(6.41)

其中 $u$ 是 $N\times N$ 单位矩阵的第 $i$ 个列向量，对应于数组的第 $i$ 个选定元素。 功率谱可以使用表达式计算:

$P_{\mathrm{LP}}(\theta )=\frac{u^H{R}^{-1}u}{{|u^H{R}^{-1}a(\theta )|}^2}$(6.42)

LP 方法的性能会因噪声的存在而降低，并且通常在 SNR 较高时运行良好。

（4）Propagator Method 传播算子法

传播子方法是一种基于线性算子的方法，该算子可以通过对流形矩阵的分块从接收数据矩阵中轻松提取出来。不需要对协方差矩阵进行特征分解或奇异值分解来估计DOA。 到达方向估计可以表述如下。

流形矩阵被划分使得我们有：

$\mathrm{A}={\left[\begin{array}{ll}{A}_1& A_2\end{array}\right]}^{\mathrm{T}}$（6.43）

其中$A_1$ 是 $K\times K$矩阵，$A_2$ 是 $(N-K)\times K$矩阵。 矩阵 A2 是 A1 的线性变换，假设 A1 是非奇异的。 子矩阵之间的关系由下式给出:

$A_2=P^H{A}_1$(6.44)

其中 P 称为投影矩阵。 根据投影矩阵，可以构造矩阵 Q 使得:

$Q=\left[P^H-I_{N-P}\right]$(6.45)

其中 I 是单位矩阵。 由于$Q^{H}A=0$ ，可以从以下关系中提取伪谱：

$P(\theta )=\frac{1}{|a(\theta )Q{Q}^{H}a(\theta )|}$(6.46)

与 MUSIC 一样，伪频谱中的峰值对应于接收信号的到达角。 从协方差矩阵 $R=\left[EX{X}^H\right]$ 估计投影矩阵 P。 然后将协方差矩阵分成两个矩阵，使得:

$\mathbf{R}=\left[\begin{array}{ll}{R}_1& R_2\end{array}\right]$(6.47)

在零噪声接收信号的情况下，投影矩阵计算如下:

${\mathbf{R}}_2={\mathbf{P}}^H{\mathbf{R}}_1$(6.48)

在正常情况下存在噪声的情况下，可以使用最小二乘法通过 Frobenius 范数(F-范数)的最小化来估计投影矩阵$\Vert R_2-P^H{R}_{1F}\Vert$ 。 这导致表达式：

$P^H=R_2{\left(R_1^H{R}_1\right)}^{-1}{R}_1^H$(6.49)

（5）Weighted Subspace Fitting (WSF)

在像 MUSIC 这样的子空间方法中，必须从有限的数据中估计协方差矩阵$R$ 。 这导致违反信号和噪声子空间之间的正交性假设。 加权子空间试图通过使用最小二乘法进行子空间估计来克服这一限制。加权子空间拟合模型可以表示为 :

$[\hat{A},\hat{T}]=\underset{A,T}{\mathrm{argmin}}{\Vert M{W}^{1/2}-A(\theta )T\Vert }_F^2$(6.50)

其中,W是正定加权矩阵，下标“F”表示 Frobenius 范数。 经过一些操作，问题减少到:

$\hat{\theta }=\underset{\mathbf{\theta }}{\mathrm{argmax}}\mathrm{Tr}\left\{{\mathbf{P}}_{\mathbf{A}(\mathbf{\theta })}{\hat{E}}_s\mathbf{W}{\hat{E}}_s^H\right\}$(6.51)

其中${\mathbf{P}}_A=\mathbf{A}{\left({\mathbf{A}}^H\mathbf{A}\right)}^{-1}{A}^H$是投影在 A 的列空间上的投影矩阵， $T_r$ 表示矩阵的迹线，${\hat{E}}_S$ 是信号的估计子空间特征向量。 可以以不同的方式选择数据的矩阵 M 表示，从而导致不同的成本函数。 T 表示 A 和 M 之间的匹配程度的度量。例如，选择 M 可以导致 WSF、多维 MUSIC、加权 ESPRIT 等。对于 WSF，选择 M 为 ${\hat{E}}_s{\mathbf{W}}_{\mathrm{opt}}^{1/2}$ ，其中 通过将协方差矩阵特征分解为信号和噪声子空间来计算最佳加权矩阵，从而得到以下表达式：

$\mathbf{R}=E_s{\Lambda }_s{E}_s^H+E_n{\Lambda }_n{E}_n^H$（6.52）

$W_{opt}$ 的估计值由下式给出:

${\mathbf{W}}_{\mathrm{opt}}={\left({\hat{\mathbf{\Lambda }}}_s-{\hat{\sigma }}^2\mathbf{I}\right)}^2{\hat{\Lambda }}_s^{-1}$(6.53)

其中 ${\sigma }^2$ 是噪声方差的估计值，可以通过平均噪声子空间特征值来计算。

6.3.5多维DOA算法

在本章中，我们主要关注 DOA 估计的一维（1D）算法。在汽车应用中，出于多种原因，通常对执行多维 DOA 估计非常感兴趣。其中包括越来越需要以动态方式测量道路物体和道路基础设施（如行人、车辆和桥梁）的高度。通常，通过前面部分中概述的一维技术来估计方位角 DOA 就足够了，以固定目标位置。由于需要增加仰角，因此有必要采用 2D DOA 估计算法。有两种方法可用。第一种方法是在方位角和仰角方向上分别执行一维 DOA，如图 10 所示。

(图 10二维DOA估计的一维方法示意图)

这种方法的优点是可以方便地使用已经过测试和信任的一维算法，缺点是必须对估计的方位角和仰角进行配对，这不是一项简单的任务，并且可能导致错误配对和错误的位置估计，进而可能导致致命事故。另一方面，第二个可用选项是使用 2D DOA 算法，从该算法可以直接提取方位角和仰角，无需配对，如图 11 所示。

(图11 2D DOA估计)

这些算法已经以 2D DBF、2D MUSIC、2D ESPRIT 等形式存在。缺点是算法所需的计算负载增加。增加的计算也使实时实现和天线设计复杂化。但是，如果有选择，人们会更喜欢 2D 算法并找到优化计算的策略。

6.3.6 最邻近估计

可以可靠分离的目标数量受阵列孔径的限制，或者简单地说，受天线元件数量的限制。然而，增加的孔径会增加天线的尺寸，这从成本和车辆集成的角度来看是不可取的。为了克服孔径问题，最近研究了虚拟阵列处理作为一种可能的解决方案，这个想法是使用多个发射天线来扩大阵列孔径，如图12 所示。

(图12 由两个中转天线单元和三个接收天线单元得到六元虚拟阵列天线构造示意图)

通过基于接收器天线单元间距适当选择发射天线间隔，可以将虚拟阵列天线单元的数量增加到 $Ntx\ast Nrx$ ，其中 $Ntx$ 是发射天线单元的数量，Nrx 是接收阵列天线单元的数目。

对于 ULA，发射天线间距由 $dtx=d\ast Nrx$ 给出，其中 $d$ 是接收天线间距。在图 12 中，我们从2个发射天线和3个接收天线中得到了一个由6个接收天线组成的虚拟阵列。

有时被称为虚拟阵列或 MIMO 雷达的优势在于能够使用小尺寸物理天线扩展孔径。将高分辨率算法应用于生成的虚拟阵列可以显着增加可分离的目标数量。此外，虚拟阵列方法可应用于旨在减少接收器天线元件数量的稀疏阵列，稀疏阵列中缺失的阵元可以作为虚拟阵元内插。

参考文献

1. Van Trees, H.: Optimum Array Processing Part IV. Wiley-Interscience (2002)

2. Karim, H., Viberg, M.: Two decades of array signal processing research. IEEE Sign. Process.
   Mag. pp. 67–94 (July 1996)

3. Gupta, P., Aditya, K., Datta, A.: Comparison of conventional and subspace based algorithms
   to estimate Direction of Arrival (DOA). In: 2016 International Conference on Communication
   and Signal Processing (ICCSP) (6–8 April 2016)

4. Cho S., et al.: A new direction-of-arrival estimation method using automotive radar sensor
   arrays. Int. J. Distrib. Sens. Netw. 13(6), 1–12 (June 2017)

5. Van Veen, B.D., Buckley, K.M.: Beamforming-A versatile approach to spatial filtering. IEEE
   ASSP Mag. pp. 4–24 (April 1988)

6. Capon, J.: High-resolution frequency-wavenumber spectrum analysis. Proc. IEEE 57,
   1408–1418 (1969)

7. Schmidt, R.O.: Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas
   Propag. 34, 276–280 (1986)

8. Golub, G.H., Van Loan, C.F.: Matrix Computations, 3rd edn. Johns Hopkins University Press,
   Baltimore, MD (1996)

9. Forsythe, K.W.: Utilizing waveform features for adaptive beamforming and direction finding
   with narrowband signals. Lincoln Lab. J. 10(2), 99–126 (1997)

10. Roy, R., Kailath, T.: ESPRIT—estimation of signal parameters via rotational invariance techniques.IEEE Trans. Acoust Speech Sign. Process. 37, 984–995 (1989)

11. Shan, T.J., Wax, M., Kailath, T.: On spatial smoothing for direction-of-arrival estimation of
    coherent sources. IEEE Trans. Acoust. Speech Sign. Process. 33(4), 806–811 (August 1985)

12. Pillai, U., Kwon, B.H.: Forward/backward spatial smoothing techniques for coherent signal
    identification. IEEE Trans. Acoust. Speech Sign. Process. 37(1), 8–15 (January 1989)

13. Munier, J., Delisle, G.Y.: Spatial analysis using new properties of the cross-spectral matrix.
    IEEE Trans. Sig. Process 39(3), 746–749 (1991)

14. Marcos, S., Marsal, A., Benidir, M.: The propagatormethod for source bearing estimation. Sig.
    Process. 42(2), 121–138 (1995)

15. Tufts, D.W., Kumaresan, R.: Estimation of frequencies of multiple sinusoids: making linear
    prediction perform like maximum likelihood. Proc. IEEE 70, 975–989 (1982)

16. Kay, S.M.: Modern Spectral Estimation: Theory and Application. Prentice Hall, Englewood
    Cliffs, NJ (1988)

17. Gamba, J.:OnNoise-CompensatedTechniques forTimeDelayEstimation. Ph.D.Thesis (2005)

18. Makhoul, J.: Linear prediction: a tutorial review. Proc. IEEE 63(4), 561–580 (1975)

19. Paulraj, A., Ottersten, B., Roy, R., Swindlehurst, A., Xu, G., Kailath, T.: Subspace methods for
    directions-of-arrival estimation. In: Bose, N.K., Rao, C.R. (eds.) Handbook of Statistics, vol.
    10, pp. 693–739. Elsevier Science Publishers B.V. (1993)

20. Viberg,M., Ottersten, B.: Sensor array processing based on subspace fitting. IEEE Trans. Sign.
    Process. 39, 1110–1121 (May 1991)

21. Chen, C.-Y., Vaidyanathan, P.P.: Minimum redundancy MIMO radars. In: 2008 IEEE International
    Symposium on Circuits and Systems, pp. 45–48 (2008)