汽车雷达在多径存在下的幽灵目标检测
摘要
共置多输入多输出（MIMO）技术已被广泛应用于汽车雷达系统，因为它能够以相对较少的发射和接收天线数量提供精确的角度估计。由于视距目标的发射方向（DOD）和到达方向（DOA）重合，MIMO 信号处理允许形成更大的虚拟阵列用于角度查找。然而，多径反射是一个主要的限制因素，雷达信号可能从障碍物反弹，创建 DOD 不等于 DOA 的回波。因此，在具有多个散射体的复杂场景中，目标的直接路径可能被其他物体的间接路径破坏，导致不准确的角度估计或产生幽灵目标。
本文将多径反射引起的幽灵存在检测作为一个复合假设决策问题：H0 假设观测仅包含共享相同（未知）DOD 和 DOA 的未知数量直接路径；H1 假设观测还包含 DOD 和 DOA 不重合的未知数量间接路径。我们采用广义似然比检验（GLRT）哲学来确定检测器结构，提供理论检测性能的闭式表达式，以及改善检测性能的凸波形优化方法。在实际场景中，GLRT 哲学的未知参数被精心设计的估计器替代。直接路径和多径的角度通过稀疏增强的压缩感知（CS）方法与 Levenberg-Marquardt（LM）优化进行估计，以在连续域中估计角度参数。
1. 引言
近年来，对更安全驾驶的需求导致了汽车雷达的显著需求增长。共置 MIMO 技术已被证明能够有效地以相对较少的天线数量提供精确的物体角度估计，使其在汽车行业中广受欢迎。
共置 MIMO 系统的一个主要挑战是多径反射，其中目标回波通过多条路径到达接收器，包括直接和间接路径。直接路径涉及信号从雷达传输到目标，然后直接反射回雷达，而间接路径可能在反射器之间多次反弹。通常，由于不同的传播延迟，距离门控可以消除我们试图检测的目标的间接路径。然而，某些间接路径的信号 DOD 不等于 DOA，因此共置 MIMO 的假设不成立。结果是，在多目标场景中，意图目标的直接路径可能被其他物体的间接路径破坏，应用经典角度查找算法可能导致角度估计精度下降和检测到幽灵目标。
为了检测幽灵目标，一些研究者利用延迟 - 多普勒域中检测的几何关系。具体来说，R. Feng 等人采用霍夫变换来探索多径检测的线性关系。F. Ross 等人通过分析移动目标的多普勒分布来检测幽灵目标。这些方法在幽灵目标速度显著时可能有效，有效利用多普勒信息可以帮助从多径中提取几何信息进行识别。值得注意的是，文献 [15] 的作者提出了一种通过波形设计抑制幽灵的新方法，该方法有效地以高精度控制不同延迟 - 多普勒单元的响应。
2. 信号模型和问题形式化
最先进的汽车雷达通常采用调频连续波（FMCW）序列来实现高分辨率的目标距离和径向速度估计，并采用共置 MIMO 技术来合成大的虚拟阵列以进行精确的角度估计。我们考虑一个共置 MIMO 雷达系统，具有 MT 个发射天线发射同样多的编码序列和 MR 个接收天线。在接收端，每个天线的信号经过通常的处理以提取每个发射天线的贡献，并合成具有 MTMR 个元素的 MIMO 通道。然后通过沿快时间和慢时间的快速傅里叶变换（FFT）处理该信号，以获得回波的延迟 - 多普勒轮廓。最后，可以构建检测目标的虚拟阵列响应以估计目标的方向。
2.1 多径场景可视化
图 1 描述：图 1 展示了多径场景的可视化。图 1(a) 显示了直接路径，其中雷达信号直接传输到位于位置 A 的目标并直接反射回来，发射和到达角度相等。图 1(b) 显示了一对一阶路径，其中信号在到达目标之前或从目标返回时在位于点 B 的反射器处反弹一次，导致 DOD 和 DOA 不相等。
多径场景可以可视化为雷达发射信号从目标和反射器反弹。如图 1 所示，目标位于位置 A，反射器位于点 B，雷达接收的信号可以采取不同的路径：
- 直接路径：雷达和目标之间的最短路径，直接路径的发射和到达角度等于目标角度，如图 1(a) 所示。
- 一阶路径：间接路径涉及在发射或到达途中在反射器处的单次反弹，导致比直接路径更长的延迟。如图 1(b) 所示，一阶路径的 DOD 不等于相应的 DOA。
- 高阶路径：间接路径在回波到达接收器之前涉及更多反弹。然而，由于目标和反射器处散射引起的衰减，高阶路径通常很弱，可以忽略。
2.2 信号模型
考虑一个 FMCW MIMO 雷达，从每个发射天线传输 L 个脉冲，并采用慢时编码作为复用方法。令 x(l)=[x1(l),x2(l),⋯ ,xMT(l)]T 为第 l 个时期 MT 个发射天线传输的码矢量，传输码矩阵可以表示为 X=[x(1),x(2),⋯ ,x(L)]∈CMT×L。
在对接收测量的快时间执行 FFT 后，我们在给定的被测延迟单元中考虑 K0 个直接路径和 K1 对一阶路径，将观测 y(l)∈CMR×1 建模为：
y(l)=∑k=1K0αkej2πfd(l−1)aR(θk)aTT(θk)x(l) +∑k=1K1βk,1ej2πfd(l−1)aR(ϕk)aTT(ϑk)x(l) +∑k=1K1βk,2ej2πfd(l−1)aR(ϑk)aTT(ϕk)x(l)+w(l)
其中：
- αk、βk,1 和 βk,2 分别表示第 k 个直接路径（k=1,2,…,K0）和第 k 对一阶路径（k=1,2,…,K1）的复振幅
- θk 表示第 k 个直接路径的 DOD，等于 DOA
- ϑk 和 ϕk 表示第 k 对一阶路径的 DOD 和 DOA，其中 ϑk≠ϕk
- fd 是归一化多普勒频率
- aT(⋅)∈CMT×1 和 aR(⋅)∈CMR×1 是导向矢量
导向矢量具体定义为：
aT(θ)=1√MT[ej2πdT,1sin(θ)/λ,ej2πdT,2sin(θ)/λ,…,ej2πdT,MTsin(θ)/λ]T
aR(ϕ)=1√MR[ej2πdR,1sin(ϕ)/λ,ej2πdR,2sin(ϕ)/λ,…,ej2πdR,MRsin(ϕ)/λ]T
其中 θ 和 ϕ 分别表示 aT(⋅) 和 aR(⋅) 的角度，λ 表示波长，dT,m 和 dR,n 表示第 m 个 TX 元素和第 n 个 RX 元素相对于参考阵列元素的相对距离。
定义 P(fd)=diag([1,ej2πfd,⋯ ,ej2πfd(L−1)])，接收数据矩阵为：
Y=∑k=1K0αkaR(θk)aTT(θk)XP(fd)+∑k=1K1βk,1aR(ϕk)aTT(ϑk)XP(fd) +∑k=1K1βk,2aR(ϑk)aTT(ϕk)XP(fd)+W
经过匹配滤波 Z=Y(XP(fd))H 并向量化后，虚拟 MIMO 阵列信号的一般模型为：
z=(Rx⊗IMR)A(Θ,Φ)β+r
其中 Rx=X*X^T，A(Θ,Φ)=AT(Θ)∘AR(Φ) 表示响应矩阵。
3. 多径检测
3.1 GLRT 检测器
在前一节概述的一般设置中，幽灵检测相当于解决一个耦合的检测 - 估计问题，其中我们必须区分复合假设 H0（观测仅包含来自未知不同方向的未知数量 K0 个直接路径）与复合替代假设 H1（观测还包含未知数量 K1 个一阶路径，每个由未知角度对表征）。
假设首先矩阵已知，我们需要解决复合二元假设检验：
{H0:z=(Rx⊗IMR)A(Θ0)α+r H1:z=(Rx⊗IMR)A(Θ,Φ)β+r}
其中 α∈CK0×1 和 β∈C(K0+2K1)×1 是未知参数。
由于 E(rrH)=σ2Rx⊗IMR，我们有 r∼CN(0,σ2Σx)，其中 Σx=Rx⊗IMR。通过噪声白化变换，测试变为：
{H0:zˉ∼CN(Σx1/2A(Θ0)α,σ2IMTMR) H1:zˉ∼CN(Σx1/2A(Θ,Φ)β,σ2IMTMR)}
其中 zˉ=Σx−1/2z。GLRT 为：
TGLRT=∥P(Θ0)zˉ∥2∥P(Θ,Φ)zˉ∥2≷H1H0λG
其中 P(Θ0) 和 P(Θ,Φ) 是相应的正交投影矩阵。
3.2 性能界限和波形优化
图 2 描述：图 2 显示了虚警概率 Pfa 与检测阈值 λG 的关系，针对 MTMR=48 的不同 (K0,K1) 值。可以观察到，随着 K1 增加，给定阈值下的虚警概率降低，这是因为假设之间的可区分性增加。
在 H0 下，测试统计量比率 X 具有 Fisher-Snedecor 分布，密度为：
fX|H0(x)=1B(2K1;m)x2K1−1(1+x)−(m+2K1)
其中 m=MTMR−K0−2K1，B(a;b) 表示参数为 a 和 b 的贝塔函数。
虚警概率和检测概率的闭式表达式为：
Pfa=1−1B(2K1;m)∑i=0m−1(−1)i(m−1i)2K1+i(1−1/λG)2K1+i
Pd=1−1B(2K1;m)∑i=0m−1(−1)i(m−1i)2K1+i(λG−1λG+ρ1)2K1+i
其中 ρ1 是一个适合的品质因数，定义为：
ρ1=σβ22K1σ2Tr(EHΣx1/2P0Σx1/2E)
图 3 描述：图 3 展示了检测概率 Pd 与 σβ2/σ2 的关系，比较了正交波形和优化波形的性能。可以看到，通过波形优化可以显著提高检测性能，特别是在低 SNR 条件下。
波形优化问题可以形式化为以下凸优化问题：
maxRx,ΠTr(EH(Rx⊗IMR)E−Π) s.t. [Rx]m,m=1,m=1,2,⋯ ,MT Λ⪰0 ∥Rx−IMT∥2≤μ Rx⪰0
这是一个半定规划（SDP）问题，可以通过凸优化方法有效求解。
4. 多径角度估计
由于测试 (13) 不可实现（矩阵 A(Θ0) 和 A(Θ,Φ) 未知），我们需要开发估计这些矩阵的方法。
4.1 H0 假设下的估计器
图 4 描述：图 4 比较了 Gauss-Newton (GN) 和 Levenberg-Marquardt (LM) 方法在优化过程中的收敛行为。图 4(a) 显示了 DOD 和 DOA 角度差异较大时（(-1.9°,-13.2°)），两种方法都表现出相似的收敛行为。图 4(b) 显示了角度差异较小时（(-1.9°,-3.2°)），GN 方法面临收敛挑战，而 LM 方法通过正则化项解决了这个问题并展现出更强的鲁棒性。
我们提出一种迭代过程来解决角度估计问题。定义 r(t) 为第 t 次迭代中的残差，初始化为 r(0)=zˉ。直接路径角度集合初始化为空集，即 Θ0(0)=∅ 且 K^0(0)=0。
在第 t 次迭代中，我们插入一条路径到集合中，K^0(t)=K^0(t−1)+1。通过评估以下式子最小化残差的 2-范数：
θ^(t)=argmaxθ(t)∈{θ1,θ2,⋯ ,θ~G}∣(r(t−1))Haˉ(θ(t))∣
并更新角度矩阵为 Θ0(t,0)=[(Θ0(t−1))T,θ^(t)]T。
随后通过 Gauss-Newton (GN) 迭代来增强此估计的准确性：
Θ0(t,i+1)=Θ0(t,i)−(H0(t,i))−1g0(t,i)
其中梯度 g0(t,i) 和 Hessian 矩阵 H0(t,i) 的表达式在附录 A 中给出。
4.2 H1 假设下的估计器
在 H1 下，算法是前一个算法的扩展，现在必须估计直接路径和一阶路径的角度。为了减少直接路径和一阶路径之间的干扰，我们分别在直接路径和一阶路径上实现估计过程。
搜索额外一阶路径对的粗略估计通过在两个均匀 G 维网格上搜索获得：
(ϑ^(t),ϕ^(t))=argmaxϑ(t)∈Ξtϕ(t)∈Ξrϑ(t)<ϕ(t)[∣(r(t−1))H(aT(ϑ(t))∘aR(ϕ(t)))∣+∣(r(t−1))H(aT(ϕ(t))∘aR(ϑ(t)))∣]
由于在 H1 假设下直接和一阶路径的混合，当 DOD 和 DOA 之间的差异不大时，GN 方法可能由于 Hessian 中的秩缺陷而导致不稳定的估计。因此，我们采用 LM 方法来更新角度估计：
h(t,i)=−(H(t,i)+μ(t,i)IK^(t))−1g(t,i)
其中 μ(t,i) 是阻尼参数，通过增益比控制：
ϱ(t,i)=Fˉ(Θˉ(t,i))−Fˉ(Θˉ(t,i)+h(t,i))12(h(t,i))H(μ(t,i)h(t,i)−g(t,i))
5. 仿真和实验结果
5.1 仿真设置
图 5 描述：图 5 展示了 MIMO 雷达天线的实际和虚拟布局。图 5(a) 显示了均匀线性阵列（ULA）配置，其中发射和接收天线均匀间隔半波长。图 5(b) 显示了稀疏线性阵列（SLA）配置，其中天线不均匀间隔以增加孔径但会引入栅瓣。虚拟阵列（蓝色圆圈）由发射和接收阵列的卷积形成。
仿真参数设置如下：
- 雷达工作频率为 79 GHz，载波波长 λ=3.8 mm
- 发射元素数量 MT=6，接收元素 MR=8
- 首先使用如图 5(a) 所示的均匀线性阵列（ULA）进行仿真
- 通过保持 MT 和 MR 恒定，我们确保检测性能的一致上界，并随后验证图 5(b) 所示 SLA 的性能
- 噪声根据方差 σ2=1 的高斯分布随机生成
- 路径幅度根据 β∼CN(0,σβ2I2K1)，α∼CN(0,σα2IK0) 生成
- 通过以 2°步长离散化角度空间 [-90°,90°] 获得网格
5.2 估计性能
图 6 描述：图 6 显示了不同算法在 ULA 和 SLA 配置下的 RMSE 性能。图 6(a) 和 6(b) 显示了 ULA 中直接路径（RMSE₀）和一阶路径（RMSE₁）的估计误差。图 6(c) 和 6(d) 显示了 SLA 的相应结果。在所有情况下，所提出的 CSCD 方法都优于基于网格的 OMP、IAA 和 LASSO 方法，特别是在高 SNR 条件下，这归因于连续域优化避免了网格失配问题。
我们评估所提出算法的角度估计均方根误差（RMSE）。注意算法返回一组估计，对应于真实路径或错误路径，如果没有接近其方向的估计，则无法检测到路径。因此，我们参考基于正确路径估计的 RMSE。具体地，一阶路径和直接路径的 RMSE 计算为：
RMSE1=√1MC∑m=1MC12∣Ω1m∣∑j∈Ω1m[(ϑj(m)−ϑ^j(m))2+(ϕj(m)−ϕ^j(m))2]
RMSE0=√1MC∑m=1MC1∣Ω0m∣∑j∈Ω0m(θj(m)−θ^j(m))2
其中 MC 是运行次数，Ω1m 和 Ω0m 是第 m 次仿真中识别的一阶路径和直接路径的索引集。
5.3 检测性能
图 7 描述：图 7 比较了 H0 场景中 a(θ),a(ψ) 的相关性（图 7a）和 H1 场景中 aT(ϑ)∘aR(ϕ),a(ψ) 的相关性（图 7b）。ULA 显示了更清晰的主瓣和较低的旁瓣，而 SLA 虽然具有更窄的波束宽度但旁瓣更高，这解释了 SLA 配置中性能下降的原因。
图 8 描述：图 8 展示了 ULA 和 SLA 配置下不同 (K1,K0) 组合的检测概率 Pd 与 σβ2/σ2 的关系。图 8(a) 和 8(b) 显示了 ULA 的结果，其中 GLRT-CSCD 的性能接近理论上界。图 8(c) 和 8(d) 显示了 SLA 的结果，其中由于较高的旁瓣，性能差距更大。在所有情况下，所提出的 GLRT-CSCD 都优于 GLRT-OMP、GLRT-LASSO 和 GLRT-IAA。
图 9 描述：图 9 比较了不同 MTMR 值下的检测性能。随着系统自由度的增加（从 MTMR=12 到 48），检测性能提高，但收益逐渐递减。这表明对于给定的 (K0,K1) 值，存在一个最优的天线配置超过该配置后性能改善有限。
5.4 实验结果
图 10 描述：图 10 展示了实验场景。图 10(a) 是实验环境的照片，显示了被混凝土墙包围的道路，这创建了多径传播的理想条件。图 10(b) 显示了雷达点云，蓝色椭圆标记了由一阶路径引起的幽灵目标。
图 11 描述：图 11 展示了使用不同方法检测和消除幽灵目标的结果。图 11(a)-11(c) 分别显示了 GLRT-OMP、GLRT-LASSO 和 GLRT-IAA 的结果，这些方法未能成功移除所有幽灵目标，并且错误地移除了一些静止目标的直接路径。图 11(d) 显示了所提出的 GLRT-CSCD 方法有效地消除了所有幽灵目标，同时保留了静止目标的直接路径。
我们使用实验数据评估所提出检测器的目标检测性能。数据由毫米波 f0=77 GHz MIMO 雷达获得，其中 MT=8 个发射天线和 MR=16 个接收天线，均匀间隔。发射侧间距为 4.5λ，接收侧间距为 4λ。
6. 结论
本文研究了汽车雷达在多径存在下的幽灵目标检测。间接路径的存在被建模为二元复合假设检验，提出了 GLRT 检测器来确定延迟 - 多普勒单元中是否存在间接路径。如果单元包含间接路径，可以移除幽灵目标并保留所需的直接路径。基于完美角度估计下 GLRT 检测性能的理论分析，我们推导了一种凸波形优化方法，可以增强检测性能。考虑到直接和间接路径角度均未知的实际场景，我们提出了一种稀疏增强的 CS 方法来估计连续域中的角度参数。仿真结果表明，所提出的算法优于现有的基于网格的估计器，从而导致更好的检测性能。提出的检测器的虚警率可以控制，ULA 情况下的检测性能接近理论界限。最后，实验结果证明了所提出方法的有效性。
附录
附录 A：H0 下梯度和 Hessian 矩阵的推导
为了清晰起见，我们在以下一些推导中省略函数的上标 (t,i) 和输入变量，即 F=F(Θ0(t,i)) 和 Aˉ0=Aˉ(Θ0(t,i))。
定义 F=fHf，其中 f=zˉ−AˉAˉ†zˉ，FFF 相对于 Θ0∈RK0×1 的梯度可以计算为：
g0=[∂F∂θ1,∂F∂θ2,…,∂F∂θK0]T
其中第 q 个元素 [g0]q 给出为 ∂F∂θq=2Re((∂f∂θq)Hf)。
根据 [49] 中的推导，我们得到：
[g0]q=−2Re{Tr{Aˉ0†zˉzˉHP0Aˉq}}
其中 Aˉq=∂Aˉ0∂θq=[0,0,…,∂aˉ∂θq,…,0]，其中 ∂aˉ∂θq=∂aˉ(θq)∂θq。
Hessian 矩阵 H0 表示 FFF 相对于 Θ0 的近似二阶偏导数。在此矩阵中，第 (q,p) 元素表示为 [H0]q,p=2Re{(∂f∂θq)H∂f∂θp} 并且可以计算为：
[H0]q,p=2Re{Tr{AˉpAˉ0†zˉzˉH(Aˉ0†)HAˉqHP0}} +2Re{Tr{AˉpHP0zˉzˉHP0AˉqAˉ0†(Aˉ0†)H}}
定义偏导矩阵 D0=[∂aˉ∂θ1,∂aˉ∂θ2,…,∂aˉ∂θK0]，则 g0 和 H0 的矩阵形式可以给出为式 (28) 和 (29)。
附录 B：H1 下梯度和 Hessian 矩阵的推导
为了清晰起见，我们在以下推导中省略上标和函数的输入变量，即 Fˉ=Fˉ(Θˉ(t,i)) 和 Aˉ=Aˉ(Θˉ(t,i),Φˉ(t,i))。在下文中，我们推导 gT 和 HTT 的矩阵表达式，gR、g0、HTR、HRR、HRT、H0T、HT0、HR0、H0R、H00 的推导遵循类似的论证，为简洁起见省略。
类似于式 (40)，我们知道 gT 的第 q 个元素可以给出为：
[gT]q=−2Re{Tr{Aˉ†zˉzˉHP1Aˉq′}} =−2Re{Tr{ΓAˉq′}}
其中 Γ=Aˉ†zˉzˉHP1∈C(2K1+K0)×MTMR，Aˉq′=∂Aˉ∂ϑq=[0,0,…,∂a1∂ϑq,…,0,…,∂a2∂ϑq,…,0]，其中：
∂a1∂ϑq=Σx1/2∂aT(ϑq)⊗aR(ϕq)∂ϑq
∂a2∂ϑq=Σx1/2∂aT(ϕq)⊗aR(ϑq)∂ϑq
我们将矩阵 Γ 划分为三个子矩阵，表示为 Γ=[Γ1,Γ2,Γ0]，其中 Γ1,Γ2∈CK1×MTMR，Γ0∈CK0×MTMR。则式 (42) 可以重写为：
[gT]q=−2Re{Γ1(q)(∂a1∂ϑq)T+Γ2(q)(∂a2∂ϑq)T}
其中 Γ1(q) 和 Γ2(q) 分别表示 Γ1 和 Γ2 的第 q 行。定义两个偏导矩阵：
DT1=[∂a1∂ϑ1,∂a1∂ϑ2,…,∂a1∂ϑK1]
DT2=[∂a2∂ϑ1,∂a2∂ϑ2,…,∂a2∂ϑK1]
我们可以得到 gT 的矩阵形式：
gT=−2Re{diag{Γ1DT1+Γ2DT2}}
类似地，定义 DR1、DR2 和 D0，得到：
gR=−2Re{diag{Γ1DR1+Γ2DR2}}
g0=−2Re{diag{Γ0D0}}
Hessian 矩阵 HTT 表示相对于 Θ1 的二阶偏导数，其中第 (q,p) 元素为：
[HTT]q,p=2Re{Tr{Aˉp′Aˉ†zˉzˉH(Aˉ†)H(Aˉq′)HP1}} +2Re{Tr{(Aˉp′)HP1zˉzˉHP1Aˉq′Aˉ†(Aˉ†)H}}
经过详细推导，我们可以得到各个 Hessian 块的矩阵形式（式 50-58）。这些表达式涉及矩阵 S 和 C 的分块形式以及偏导矩阵 DT、DR 和 D0 的各种组合。
