汽车雷达多径环境下的幽灵目标检测技术

在对接收测量的快时间执行 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) + ∑{k=1}^{K1} β_{k,1} e^{j2πfd(l−1)} a_R(ϕ_k)a_T^T(ϑ_k)x(l) + ∑{k=1}^{K1} β{k,2} e^{j2πfd(l−1)} a_R(ϑ_k)a_T^T(ϕ_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 是归一化多普勒频率
• a_T(·) ∈ C^{MT × 1} 和 a_R(·) ∈ C^{MR × 1} 是导向矢量

导向矢量具体定义为：

a_T(θ) = \frac{1}{\sqrt{M_T}}[e^{j2πd_{T,1}sin(θ)/λ}, e^{j2πd_{T,2}sin(θ)/λ}, …, e^{j2πd_{T,M_T}sin(θ)/λ}]^T

a_R(ϕ) = \frac{1}{\sqrt{M_R}}[e^{j2πd_{R,1}sin(ϕ)/λ}, e^{j2πd_{R,2}sin(ϕ)/λ}, …, e^{j2πd_{R,M_R}sin(ϕ)/λ}]^T

其中 θ 和 ϕ 分别表示 a_T(·) 和 a_R(·) 的角度，λ 表示波长，d_{T,m} 和 d_{R,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)X P(fd) + ∑{k=1}^{K1} β_{k,1} a_R(ϕ_k)a_T^T(ϑ_k)X P(fd) + ∑{k=1}^{K1} β{k,2} a_R(ϑ_k)a_T^T(ϕ_k)X P(fd) + W

经过匹配滤波 Z = Y(XP(fd))^H 并向量化后，虚拟 MIMO 阵列信号的一般模型为：

z = (Rx ⊗ I_{MR})A(Θ,Φ)β + r

其中 Rx = X*X^T，A(Θ,Φ) = A_T(Θ) ∘ A_R(Φ) 表示响应矩阵。

3. 多径检测

3.1 GLRT 检测器

在前一节概述的一般设置中，幽灵检测相当于解决一个耦合的检测 - 估计问题，其中我们必须区分复合假设 H0（观测仅包含来自未知不同方向的未知数量 K0 个直接路径）与复合替代假设 H1（观测还包含未知数量 K1 个一阶路径，每个由未知角度对表征）。

假设首先矩阵已知，我们需要解决复合二元假设检验：

{H0: z = (Rx ⊗ I_{MR})A(Θ0)α + r H1: z = (Rx ⊗ I_{MR})A(Θ,Φ)β + r}

其中 α ∈ C^{K0 × 1} 和 β ∈ C^{(K0+2K1) × 1} 是未知参数。

由于 E(rr^H) = σ^2Rx ⊗ I_{MR}，我们有 r ~ CN(0, σ^2Σ_x)，其中 Σ_x = Rx ⊗ I_{MR}。通过噪声白化变换，测试变为：

{H0: z_bar ~ CN(Σ_x^{1/2}A(Θ0)α, σ^2I_{MTMR}) H1: z_bar ~ CN(Σ_x^{1/2}A(Θ,Φ)β, σ^2I_{MTMR})}

其中 z_bar = Σ_x^{-1/2}z。GLRT 为：

T_GLRT = ||P(Θ0)z_bar||^2 / ||P(Θ,Φ)z_bar||^2 ≷_{H0}^{H1} λ_G

其中 P(Θ0) 和 P(Θ,Φ) 是相应的正交投影矩阵。

3.2 性能界限和波形优化

图 2 描述：图 2 显示了虚警概率 P_fa 与检测阈值 λ_G 的关系，针对 MTMR = 48 的不同 (K0, K1) 值。可以观察到，随着 K1 增加，给定阈值下的虚警概率降低，这是因为假设之间的可区分性增加。

在 H0 下，测试统计量比率 X 具有 Fisher-Snedecor 分布，密度为：

f_{X|H0}(x) = \frac{1}{B(2K1; m)} x^{2K1−1}(1 + x)^{-(m+2K1)}

其中 m = MTMR − K0 − 2K1，B(a;b) 表示参数为 a 和 b 的贝塔函数。

虚警概率和检测概率的闭式表达式为：

P_fa = 1 - \frac{1}{B(2K1; m)} ∑_{i=0}^{m−1} (-1)^i \binom{m−1}{i} \frac{2K1 + i}{(1 - 1/λ_G)^{2K1+i}}

P_d = 1 - \frac{1}{B(2K1; m)} ∑_{i=0}^{m−1} (-1)^i \binom{m−1}{i} \frac{2K1 + i}{((λ_G - 1)/(λ_G + ρ_1))^{2K1+i}}

其中 ρ_1 是一个适合的品质因数，定义为：

ρ_1 = \frac{σ_β^2}{2K1σ^2} Tr(E^HΣ_x^{1/2}P_0Σ_x^{1/2}E)

图 3 描述：图 3 展示了检测概率 P_d 与 σ_β^2/σ^2 的关系，比较了正交波形和优化波形的性能。可以看到，通过波形优化可以显著提高检测性能，特别是在低 SNR 条件下。

波形优化问题可以形式化为以下凸优化问题：

max_{Rx, Π} Tr(E^H(Rx ⊗ I_{MR})E - Π) s.t. [Rx]{m,m} = 1, m = 1, 2, ..., MT Λ ⪰ 0 ||Rx - I{MT}||^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_bar。直接路径角度集合初始化为空集，即 Θ0^(0) = ∅ 且 K_hat_0^(0) = 0。

在第 t 次迭代中，我们插入一条路径到集合中，K_hat_0^(t) = K_hat_0^(t-1) + 1。通过评估以下式子最小化残差的 2-范数：

θ_hat^(t) = arg max_{θ^(t) ∈ {θ~_1, θ~_2, ..., θ~_G}} |(r^(t-1))^H a_bar(θ^(t))|

并更新角度矩阵为 Θ0^(t,0) = [(Θ0^(t-1))^T, θ_hat^(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 中给出。

4.2 H1 假设下的估计器

在 H1 下，算法是前一个算法的扩展，现在必须估计直接路径和一阶路径的角度。为了减少直接路径和一阶路径之间的干扰，我们分别在直接路径和一阶路径上实现估计过程。

搜索额外一阶路径对的粗略估计通过在两个均匀 G 维网格上搜索获得：

(ϑ_hat^(t), ϕ_hat^(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_hat^(t)})^{-1}g^(t,i)

其中 μ^(t,i) 是阻尼参数，通过增益比控制：

ϱ^(t,i) = (F_bar(Θ_bar^(t,i)) - F_bar(Θ_bar^(t,i) + h^(t,i))) / (1/2(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, σ_β^2 I_{2K1})，α ~ CN(0, σ_α^2 I_{K0}) 生成
• 通过以 2°步长离散化角度空间 [-90°, 90°] 获得网格

5.2 估计性能

图 6 描述：图 6 显示了不同算法在 ULA 和 SLA 配置下的 RMSE 性能。图 6(a) 和 6(b) 显示了 ULA 中直接路径（RMSE_0）和一阶路径（RMSE_1）的估计误差。图 6(c) 和 6(d) 显示了 SLA 的相应结果。在所有情况下，所提出的 CSCD 方法都优于基于网格的 OMP、IAA 和 LASSO 方法，特别是在高 SNR 条件下，这归因于连续域优化避免了网格失配问题。

我们评估所提出算法的角度估计均方根误差（RMSE）。注意算法返回一组估计，对应于真实路径或错误路径，如果没有接近其方向的估计，则无法检测到路径。因此，我们参考基于正确路径估计的 RMSE。具体地，一阶路径和直接路径的 RMSE 计算为：

RMSE_1 = sqrt(1/M_C ∑{m=1}^{M_C} 1/(2|Ω_1^m|) ∑{j ∈ Ω_1^m} [(ϑ_j^(m) - ϑ_hat_j^(m))^2 + (ϕ_j^(m) - ϕ_hat_j^(m))^2])

RMSE_0 = sqrt(1/M_C ∑{m=1}^{M_C} 1/|Ω_0^m| ∑{j ∈ Ω_0^m} (θ_j^(m) - θ_hat_j^(m))^2)

其中 M_C 是运行次数，Ω_1^m 和 Ω_0^m 是第 m 次仿真中识别的一阶路径和直接路径的索引集。

5.3 检测性能

图 7 描述：图 7 比较了 H0 场景中 a(θ), a(ψ) 的相关性（图 7a）和 H1 场景中 a_T(ϑ) ∘ a_R(ϕ), a(ψ) 的相关性（图 7b）。ULA 显示了更清晰的主瓣和较低的旁瓣，而 SLA 虽然具有更窄的波束宽度但旁瓣更高，这解释了 SLA 配置中性能下降的原因。

图 8 描述：图 8 展示了 ULA 和 SLA 配置下不同 (K1, K0) 组合的检测概率 P_d 与 σ_β^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)-(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_bar_0 = A_bar(Θ0^(t,i))。

定义 F = f^H f，其中 f = z_bar - A_bar A_bar^† z_bar，F 相对于 Θ0 ∈ R^{K0 × 1} 的梯度可以计算为：

g_0 = [∂F/∂θ_1, ∂F/∂θ_2, …, ∂F/∂θ_{K0}]^T

其中第 q 个元素 [g_0]_q 给出为 ∂F/∂θ_q = 2Re((∂f/∂θ_q)^H f)。

根据 [49] 中的推导，我们得到：

[g_0]_q = -2Re{Tr{A_bar_0^† z_bar z_bar^H P_0 A_bar_q}}

其中 A_bar_q = ∂A_bar_0/∂θ_q = [0, 0, …, ∂a_bar/∂θ_q, …, 0]，其中 ∂a_bar/∂θ_q = ∂a_bar(θ_q)/∂θ_q。

Hessian 矩阵 H_0 表示 F 相对于 Θ0 的近似二阶偏导数。在此矩阵中，第 (q,p) 元素表示为 [H_0]_{q,p} = 2Re{(∂f/∂θ_q)^H ∂f/∂θ_p} 并且可以计算为：

[H_0]_{q,p} = 2Re{Tr{A_bar_p A_bar_0^† z_bar z_bar^H (A_bar_0^†)^H A_bar_q^H P_0}} + 2Re{Tr{A_bar_p^H P_0 z_bar z_bar^H P_0 A_bar_q A_bar_0^† (A_bar_0^†)^H}}

定义偏导矩阵 D_0 = [∂a_bar/∂θ_1, ∂a_bar/∂θ_2, …, ∂a_bar/∂θ_{K0}]，则 g_0 和 H_0 的矩阵形式可以给出为式 (28) 和 (29)。

附录 B：H1 下梯度和 Hessian 矩阵的推导

为了清晰起见，我们在以下推导中省略上标和函数的输入变量，即 F_bar = F_bar(Θ_bar^(t,i)) 和 A_bar = A_bar(Θ_bar^(t,i), Φ_bar^(t,i))。在下文中，我们推导 g_T 和 H_TT 的矩阵表达式，g_R、g_0、H_TR、H_RR、H_RT、H_0T、H_T0、H_R0、H_0R、H_00 的推导遵循类似的论证，为简洁起见省略。

类似于式 (40)，我们知道 g_T 的第 q 个元素可以给出为：

[g_T]_q = -2Re{Tr{A_bar^† z_bar z_bar^H P_1 A_bar_q'}}

= -2Re{Tr{Γ A_bar_q'}}

其中 Γ = A_bar^† z_bar z_bar^H P_1 ∈ C^{(2K1+K0) × MTMR}，A_bar_q' = ∂A_bar/∂ϑ_q = [0, 0, …, ∂a_1/∂ϑ_q, …, 0, …, ∂a_2/∂ϑ_q, …, 0]，其中：

∂a_1/∂ϑ_q = Σ_x^{1/2} ∂(a_T(ϑ_q) ⊗ a_R(ϕ_q))/∂ϑ_q

∂a_2/∂ϑ_q = Σ_x^{1/2} ∂(a_T(ϕ_q) ⊗ a_R(ϑ_q))/∂ϑ_q

我们将矩阵 Γ 划分为三个子矩阵，表示为 Γ = [Γ_1, Γ_2, Γ_0]，其中 Γ_1, Γ_2 ∈ C^{K1 × MTMR}，Γ_0 ∈ C^{K0 × MTMR}。则式 (42) 可以重写为：

[g_T]_q = -2Re{Γ_1^(q) (∂a_1/∂ϑ_q)^T + Γ_2^(q) (∂a_2/∂ϑ_q)^T}

其中 Γ_1^(q) 和 Γ_2^(q) 分别表示 Γ_1 和 Γ_2 的第 q 行。定义两个偏导矩阵：

D_T1 = [∂a_1/∂ϑ_1, ∂a_1/∂ϑ_2, …, ∂a_1/∂ϑ_{K1}]

D_T2 = [∂a_2/∂ϑ_1, ∂a_2/∂ϑ_2, …, ∂a_2/∂ϑ_{K1}]

我们可以得到 g_T 的矩阵形式：

g_T = -2Re{diag{Γ_1 D_T1 + Γ_2 D_T2}}

类似地，定义 D_R1、D_R2 和 D_0，得到：

g_R = -2Re{diag{Γ_1 D_R1 + Γ_2 D_R2}}

g_0 = -2Re{diag{Γ_0 D_0}}

Hessian 矩阵 H_TT 表示相对于 Θ_1 的二阶偏导数，其中第 (q,p) 元素为：

[H_TT]_{q,p} = 2Re{Tr{A_bar_p' A_bar^† z_bar z_bar^H (A_bar^†)^H (A_bar_q')^H P_1}} + 2Re{Tr{(A_bar_p')^H P_1 z_bar z_bar^H P_1 A_bar_q' A_bar^† (A_bar^†)^H}}

经过详细推导，我们可以得到各个 Hessian 块的矩阵形式（式 50-58）。这些表达式涉及矩阵 S 和 C 的分块形式以及偏导矩阵 D_T、D_R 和 D_0 的各种组合。