MATLAB / Octave算法
汽车雷达多径环境下幽灵目标检测算法研究
对汽车雷达多径反射导致的幽灵目标问题，提出基于广义似然比检验（GLRT）的检测方法。通过稀疏增强压缩感知与 Levenberg-Marquardt 优化估计角度参数，推导了理论性能界限并设计了凸波形优化方案。仿真与实验表明，该方法优于传统网格基估计器，能有效消除幽灵目标并保留真实路径。
对汽车雷达多径反射导致的幽灵目标问题，提出基于广义似然比检验（GLRT）的检测方法。通过稀疏增强压缩感知与 Levenberg-Marquardt 优化估计角度参数，推导了理论性能界限并设计了凸波形优化方案。仿真与实验表明，该方法优于传统网格基估计器，能有效消除幽灵目标并保留真实路径。
D. Sharif, S. Murtala and G. S. Choi, "A Survey of Automotive Radar Misalignment Detection Techniques," in IEEE Access, vol. 13, pp. 123314-123324, 2025.
共置多输入多输出（MIMO）技术已被广泛应用于汽车雷达系统，因为它能够以相对较少的发射和接收天线数量提供精确的角度估计。由于视距目标的发射方向（DOD）和到达方向（DOA）重合，MIMO 信号处理允许形成更大的虚拟阵列用于角度查找。然而，多径反射是一个主要的限制因素，雷达信号可能从障碍物反弹，创建 DOD 不等于 DOA 的回波。因此，在具有多个散射体的复杂场景中，目标的直接路径可能被其他物体的间接路径破坏，导致不准确的角度估计或产生幽灵目标。
本文将多径反射引起的幽灵存在检测作为一个复合假设决策问题：H0 假设观测仅包含共享相同（未知）DOD 和 DOA 的未知数量直接路径；H1 假设观测还包含 DOD 和 DOA 不重合的未知数量间接路径。我们采用广义似然比检验（GLRT）哲学来确定检测器结构，提供理论检测性能的闭式表达式，以及改善检测性能的凸波形优化方法。在实际场景中，GLRT 哲学的未知参数被精心设计的估计器替代。直接路径和多径的角度通过稀疏增强的压缩感知（CS）方法与 Levenberg-Marquardt（LM）优化进行估计，以在连续域中估计角度参数。
近年来，对更安全驾驶的需求导致了汽车雷达的显著需求增长。共置 MIMO 技术已被证明能够有效地以相对较少的天线数量提供精确的物体角度估计，使其在汽车行业中广受欢迎。
共置 MIMO 系统的一个主要挑战是多径反射，其中目标回波通过多条路径到达接收器，包括直接和间接路径。直接路径涉及信号从雷达传输到目标，然后直接反射回雷达，而间接路径可能在反射器之间多次反弹。通常，由于不同的传播延迟，距离门控可以消除我们试图检测的目标的间接路径。然而，某些间接路径的信号 DOD 不等于 DOA，因此共置 MIMO 的假设不成立。结果是，在多目标场景中，意图目标的直接路径可能被其他物体的间接路径破坏，应用经典角度查找算法可能导致角度估计精度下降和检测到幽灵目标。
为了检测幽灵目标，一些研究者利用延迟 - 多普勒域中检测的几何关系。具体来说，R. Feng 等人采用霍夫变换来探索多径检测的线性关系。F. Ross 等人通过分析移动目标的多普勒分布来检测幽灵目标。这些方法在幽灵目标速度显著时可能有效，有效利用多普勒信息可以帮助从多径中提取几何信息进行识别。
最先进的汽车雷达通常采用调频连续波（FMCW）序列来实现高分辨率的目标距离和径向速度估计，并采用共置 MIMO 技术来合成大的虚拟阵列以进行精确的角度估计。我们考虑一个共置 MIMO 雷达系统，具有 MT 个发射天线发射同样多的编码序列和 MR 个接收天线。在接收端，每个天线的信号经过通常的处理以提取每个发射天线的贡献，并合成具有 MTMR 个元素的 MIMO 通道。然后通过沿快时间和慢时间的快速傅里叶变换（FFT）处理该信号，以获得回波的延迟 - 多普勒轮廓。最后，可以构建检测目标的虚拟阵列响应以估计目标的方向。
![图 1 多径场景可视化]
图 1 描述：图 1 展示了多径场景的可视化。图 1(a) 显示了直接路径，其中雷达信号直接传输到位于位置 A 的目标并直接反射回来，发射和到达角度相等。图 1(b) 显示了一对一阶路径，其中信号在到达目标之前或从目标返回时在位于点 B 的反射器处反弹一次，导致 DOD 和 DOA 不相等。
多径场景可以可视化为雷达发射信号从目标和反射器反弹。如图 1 所示，目标位于位置 A，反射器位于点 B，雷达接收的信号可以采取不同的路径：
考虑一个 FMCW MIMO 雷达，从每个发射天线传输 L 个脉冲，并采用慢时编码作为复用方法。令 x(l) = [x_1(l), x_2(l), ..., x_MT(l)]^T 为第 l 个时期 MT 个发射天线传输的码矢量，传输码矩阵可以表示为 X = [x(1), x(2), ..., x(L)] ∈ C^{MT×L}。
在对接收测量的快时间执行 FFT 后，我们在给定的被测延迟单元中考虑 K0 个直接路径和 K1 对一阶路径，将观测 y(l) ∈ C^{MR×1} 建模为：
y(l) = Σ_{k=1}^{K0} α_k e^{j2πfd(l-1)} a_R(θ_k)a_T^T(θ_k)x(l)
其中：
导向矢量具体定义为：
a_T(θ) = (1/√MT)[e^{j2πdT,1sin(θ)/λ}, ..., e^{j2πdT,MTsin(θ)/λ}]^T a_R(ϕ) = (1/√MR)[e^{j2πdR,1sin(ϕ)/λ}, ..., e^{j2πdR,MRsin(ϕ)/λ}]^T
其中 θ 和 ϕ 分别表示 a_T(·) 和 a_R(·) 的角度，λ 表示波长，dT,m 和 dR,n 表示第 m 个 TX 元素和第 n 个 RX 元素相对于参考阵列元素的相对距离。
定义 P(fd) = diag([1, e^{j2πfd}, ..., e^{j2πfd(L-1)})，接收数据矩阵为：
Y = Σ_{k=1}^{K0} α_k a_R(θ_k)a_T^T(θ_k)XP(fd) + Σ_{k=1}^{K1} β_{k,1} a_R(ϕ_k)a_T^T(ϑ_k)XP(fd)
经过匹配滤波 Z = Y(XP(fd))^H 并向量化后，虚拟 MIMO 阵列信号的一般模型为：
z = (Rx ⊗ IMR)A(Θ,Φ)β + r
其中 Rx = X*X^T，A(Θ,Φ) = AT(Θ) ∘ AR(Φ) 表示响应矩阵。
在前一节概述的一般设置中，幽灵检测相当于解决一个耦合的检测 - 估计问题，其中我们必须区分复合假设 H0（观测仅包含来自未知不同方向的未知数量 K0 个直接路径）与复合替代假设 H1（观测还包含未知数量 K1 个一阶路径，每个由未知角度对表征）。
假设首先矩阵已知，我们需要解决复合二元假设检验：
{H0: z = (Rx ⊗ IMR)A(Θ0)α + r H1: z = (Rx ⊗ IMR)A(Θ,Φ)β + r}
其中 α ∈ C^{K0×1} 和 β ∈ C^{(K0+2K1)×1} 是未知参数。
由于 E(rr^H) = σ^2Rx ⊗ IMR，我们有 r ~ CN(0, σ^2Σx)，其中 Σx = Rx ⊗ IMR。通过噪声白化变换，测试变为：
{H0: z_bar ~ CN(Σx^{1/2}A(Θ0)α, σ^2IMTMR) H1: z_bar ~ CN(Σx^{1/2}A(Θ,Φ)β, σ^2IMTMR)}
其中 z_bar = Σx^{-1/2}z。GLRT 为：
T_GLRT = ||P(Θ0)z_bar||^2 / ||P(Θ,Φ)z_bar||^2 ≷_{H0}^{H1} λ_G
其中 P(Θ0) 和 P(Θ,Φ) 是相应的正交投影矩阵。
![图 2 虚警概率与阈值关系]
图 2 描述：图 2 显示了虚警概率 Pfa 与检测阈值 λG 的关系，针对 MTMR=48 的不同 (K0, K1) 值。可以观察到，随着 K1 增加，给定阈值下的虚警概率降低，这是因为假设之间的可区分性增加。
在 H0 下，测试统计量比率 X 具有 Fisher-Snedecor 分布，密度为：
f_{X|H0}(x) = (1/B(2K1; m)) x^{2K1-1}(1 + x)^{-(m+2K1)}
其中 m = MTMR - K0 - 2K1，B(a;b) 表示参数为 a 和 b 的贝塔函数。
虚警概率和检测概率的闭式表达式为：
Pfa = 1 - (1/B(2K1; m)) Σ_{i=0}^{m-1} (-1)^i binom(m-1, i) (2K1+i)/(1-1/λG)^{2K1+i} Pd = 1 - (1/B(2K1; m)) Σ_{i=0}^{m-1} (-1)^i binom(m-1, i) (2K1+i)/((λG-1)/(λG+ρ1))^{2K1+i}
其中 ρ1 是一个适合的品质因数。
![图 3 检测概率与信噪比关系]
图 3 描述：图 3 展示了检测概率 Pd 与 σβ^2/σ^2 的关系，比较了正交波形和优化波形的性能。可以看到，通过波形优化可以显著提高检测性能，特别是在低 SNR 条件下。
波形优化问题可以形式化为以下凸优化问题：
max_{Rx, Π} Tr(E^H(Rx ⊗ IMR)E - Π) s.t. [Rx]_{m,m} = 1, m = 1, 2, ..., MT Λ ⪰ 0 ||Rx - IMT||^2 ≤ μ Rx ⪰ 0
这是一个半定规划（SDP）问题，可以通过凸优化方法有效求解。
由于测试不可实现（矩阵 A(Θ0) 和 A(Θ,Φ) 未知），我们需要开发估计这些矩阵的方法。
![图 4 GN 与 LM 收敛行为对比]
图 4 描述：图 4 比较了 Gauss-Newton (GN) 和 Levenberg-Marquardt (LM) 方法在优化过程中的收敛行为。图 4(a) 显示了 DOD 和 DOA 角度差异较大时，两种方法都表现出相似的收敛行为。图 4(b) 显示了角度差异较小时，GN 方法面临收敛挑战，而 LM 方法通过正则化项解决了这个问题并展现出更强的鲁棒性。
我们提出一种迭代过程来解决角度估计问题。定义 r^(t) 为第 t 次迭代中的残差，初始化为 r^(0) = z_bar。直接路径角度集合初始化为空集，即 Θ0^(0) = ∅ 且 K0^(0) = 0。
在第 t 次迭代中，我们插入一条路径到集合中，K0^(t) = K0^(t-1) + 1。通过评估以下式子最小化残差的 2-范数：
θ^(t) = arg max_{θ^(t) ∈ {θ1, ..., θG}} |(r^(t-1))^H a_bar(θ^(t))|
并更新角度矩阵为 Θ0^(t,0) = [(Θ0^(t-1))^T, θ^(t)]^T。
随后通过 Gauss-Newton (GN) 迭代来增强此估计的准确性：
Θ0^(t,i+1) = Θ0^(t,i) - (H0^(t,i))^{-1}g0^(t,i)
其中梯度 g0^(t,i) 和 Hessian 矩阵 H0^(t,i) 的表达式在附录 A 中给出。
在 H1 下，算法是前一个算法的扩展，现在必须估计直接路径和一阶路径的角度。为了减少直接路径和一阶路径之间的干扰，我们分别在直接路径和一阶路径上实现估计过程。
搜索额外一阶路径对的粗略估计通过在两个均匀 G 维网格上搜索获得：
(ϑ^(t), ϕ^(t)) = arg max_{ϑ^(t)∈Ξt, ϕ^(t)∈Ξr, ϑ^(t)<ϕ^(t)} [| (r^(t-1))^H(a_T(ϑ^(t))∘a_R(ϕ^(t))) | + | (r^(t-1))^H(a_T(ϕ^(t))∘a_R(ϑ^(t))) |]
由于在 H1 假设下直接和一阶路径的混合，当 DOD 和 DOA 之间的差异不大时，GN 方法可能由于 Hessian 中的秩缺陷而导致不稳定的估计。因此，我们采用 LM 方法来更新角度估计：
h^(t,i) = -(H^(t,i) + μ^(t,i)I_K^(t))^{-1}g^(t,i)
其中 μ^(t,i) 是阻尼参数，通过增益比控制。
![图 5 MIMO 雷达天线布局]
图 5 描述：图 5 展示了 MIMO 雷达天线的实际和虚拟布局。图 5(a) 显示了均匀线性阵列（ULA）配置，其中发射和接收天线均匀间隔半波长。图 5(b) 显示了稀疏线性阵列（SLA）配置，其中天线不均匀间隔以增加孔径但会引入栅瓣。虚拟阵列由发射和接收阵列的卷积形成。
仿真参数设置如下：
![图 6 不同算法 RMSE 性能]
图 6 描述：图 6 显示了不同算法在 ULA 和 SLA 配置下的 RMSE 性能。图 6(a) 和 6(b) 显示了 ULA 中直接路径（RMSE0）和一阶路径（RMSE1）的估计误差。图 6(c) 和 6(d) 显示了 SLA 的相应结果。在所有情况下，所提出的 CSCD 方法都优于基于网格的 OMP、IAA 和 LASSO 方法，特别是在高 SNR 条件下，这归因于连续域优化避免了网格失配问题。
我们评估所提出算法的角度估计均方根误差（RMSE）。注意算法返回一组估计，对应于真实路径或错误路径，如果没有接近其方向的估计，则无法检测到路径。因此，我们参考基于正确路径估计的 RMSE。具体地，一阶路径和直接路径的 RMSE 计算为：
RMSE1 = sqrt( (1/MC) Σ_{m=1}^{MC} (1/(2|Ω1^m|)) Σ_{j∈Ω1^m} [(ϑ_j^(m) - ϑ̂_j^(m))^2 + (ϕ_j^(m) - ϕ̂_j^(m))^2] ) RMSE0 = sqrt( (1/MC) Σ_{m=1}^{MC} (1/|Ω0^m|) Σ_{j∈Ω0^m} (θ_j^(m) - θ̂_j^(m))^2 )
其中 MC 是运行次数，Ω1^m 和 Ω0^m 是第 m 次仿真中识别的一阶路径和直接路径的索引集。
![图 7 相关性对比]
图 7 描述：图 7 比较了 H0 场景中 a(θ), a(ψ) 的相关性和 H1 场景中 a_T(ϑ)∘a_R(ϕ), a(ψ) 的相关性。ULA 显示了更清晰的主瓣和较低的旁瓣，而 SLA 虽然具有更窄的波束宽度但旁瓣更高，这解释了 SLA 配置中性能下降的原因。
![图 8 检测概率与信噪比]
图 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 描述：图 9 比较了不同 MTMR 值下的检测性能。随着系统自由度的增加（从 MTMR=12 到 48），检测性能提高，但收益逐渐递减。这表明对于给定的 (K0, K1) 值，存在一个最优的天线配置超过该配置后性能改善有限。
![图 10 实验场景]
图 10 描述：图 10 展示了实验场景。图 10(a) 是实验环境的照片，显示了被混凝土墙包围的道路，这创建了多径传播的理想条件。图 10(b) 显示了雷达点云，蓝色椭圆标记了由一阶路径引起的幽灵目标。
![图 11 幽灵目标检测结果]
图 11 描述：图 11 展示了使用不同方法检测和消除幽灵目标的结果。图 11(a)-(c) 分别显示了 GLRT-OMP、GLRT-LASSO 和 GLRT-IAA 的结果，这些方法未能成功移除所有幽灵目标，并且错误地移除了一些静止目标的直接路径。图 11(d) 显示了所提出的 GLRT-CSCD 方法有效地消除了所有幽灵目标，同时保留了静止目标的直接路径。
我们使用实验数据评估所提出检测器的目标检测性能。数据由毫米波 f0 = 77 GHz MIMO 雷达获得，其中 MT = 8 个发射天线和 MR = 16 个接收天线，均匀间隔。发射侧间距为 4.5λ，接收侧间距为 4λ。
本文研究了汽车雷达在多径存在下的幽灵目标检测。间接路径的存在被建模为二元复合假设检验，提出了 GLRT 检测器来确定延迟 - 多普勒单元中是否存在间接路径。如果单元包含间接路径，可以移除幽灵目标并保留所需的直接路径。基于完美角度估计下 GLRT 检测性能的理论分析，我们推导了一种凸波形优化方法，可以增强检测性能。考虑到直接和间接路径角度均未知的实际场景，我们提出了一种稀疏增强的 CS 方法来估计连续域中的角度参数。仿真结果表明，所提出的算法优于现有的基于网格的估计器，从而导致更好的检测性能。提出的检测器的虚警率可以控制，ULA 情况下的检测性能接近理论界限。最后，实验结果证明了所提出方法的有效性。
为了清晰起见，我们在以下一些推导中省略函数的上标 (t,i) 和输入变量，即 F = F(Θ0^(t,i)) 和 Ā0 = Ā(Θ0^(t,i))。
定义 F = f^Hf，其中 f = z_bar - ĀĀ†z_bar，FFF 相对于 Θ0 ∈ R^{K0×1} 的梯度可以计算为：
g0 = [∂F/∂θ1, ∂F/∂θ2, ..., ∂F/∂θK0]^T
其中第 q 个元素 [g0]_q 给出为 ∂F/∂θq = 2Re((∂f/∂θq)^H f)。
根据 [49] 中的推导，我们得到：
[g0]_q = -2Re{Tr{Ā0† z_bar z_bar^H P0 Āq}}
其中 Āq = ∂Ā0/∂θq = [0, 0, ..., ∂ā/∂θq, ..., 0]，其中 ∂ā/∂θq = ∂ā(θq)/∂θq。
Hessian 矩阵 H0 表示 FFF 相对于 Θ0 的近似二阶偏导数。在此矩阵中，第 (q,p) 元素表示为 [H0]_{q,p} = 2Re{(∂f/∂θq)^H ∂f/∂θp} 并且可以计算为：
[H0]_{q,p} = 2Re{Tr{Āp Ā0† z_bar z_bar^H (Ā0†)^H Āq^H P0}}
定义偏导矩阵 D0 = [∂ā/∂θ1, ∂ā/∂θ2, ..., ∂ā/∂θK0]，则 g0 和 H0 的矩阵形式可以给出为式 (28) 和 (29)。
为了清晰起见，我们在以下推导中省略上标和函数的输入变量，即 F_bar = F_bar(Θ_bar^(t,i)) 和 Ā = Ā(Θ_bar^(t,i), Φ_bar^(t,i))。在下文中，我们推导 gT 和 HTT 的矩阵表达式，gR、g0、HTR、HRR、HRT、H0T、HT0、HR0、H0R、H00 的推导遵循类似的论证，为简洁起见省略。
类似于式 (40)，我们知道 gT 的第 q 个元素可以给出为：
[gT]_q = -2Re{Tr{Ā† z_bar z_bar^H P1 Āq'}} = -2Re{Tr{Γ Āq'}}
其中 Γ = Ā† z_bar z_bar^H P1 ∈ C^{(2K1+K0)×MTMR}，Āq' = ∂Ā/∂ϑq = [0, 0, ..., ∂a1/∂ϑq, ..., 0, ..., ∂a2/∂ϑq, ..., 0]，其中：
∂a1/∂ϑq = Σx^{1/2} ∂(aT(ϑq)⊗aR(ϕq))/∂ϑq ∂a2/∂ϑq = Σx^{1/2} ∂(aT(ϕq)⊗aR(ϑq))/∂ϑq
我们将矩阵 Γ 划分为三个子矩阵，表示为 Γ = [Γ1, Γ2, Γ0]，其中 Γ1, Γ2 ∈ C^{K1×MTMR}，Γ0 ∈ C^{K0×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{Γ1 DT1 + Γ2 DT2}}
类似地，定义 DR1、DR2 和 D0，得到：
gR = -2Re{diag{Γ1 DR1 + Γ2 DR2}} g0 = -2Re{diag{Γ0 D0}}
Hessian 矩阵 HTT 表示相对于 Θ1 的二阶偏导数，其中第 (q,p) 元素为：
[HTT]_{q,p} = 2Re{Tr{Āp' Ā† z_bar z_bar^H (Ā†)^H (Āq')^H P1}}
经过详细推导，我们可以得到各个 Hessian 块的矩阵形式（式 50-58）。这些表达式涉及矩阵 S 和 C 的分块形式以及偏导矩阵 DT、DR 和 D0 的各种组合。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML 转 Markdown 互为补充。 在线工具，Markdown 转 HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML 转 Markdown在线工具，online
通过删除不必要的空白来缩小和压缩JSON。 在线工具，JSON 压缩在线工具，online