汽车雷达多径幽灵目标检测与 GLRT 算法研究

其中：

• α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)=arg⁡max⁡θ(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))=arg⁡max⁡ϑ(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 的各种组合。