GRAIL-SCML：一种几何正则化自适应集成似然的波达方向估计算法

TL;DR：GRAIL-SCML 是一种融合确定性最大似然、协方差域拟合、弱几何正则化与证据支撑集选择的高精度 DOA 估计算法，在低信噪比、单快拍、强相关性、非均匀噪声等极端场景下全面超越传统方法。
原文链接

---

1. 引言：什么是DOA估计？

波达方向估计（Direction-of-Arrival, DOA）是阵列信号处理中最核心的问题之一，广泛应用于雷达、声呐、无线通信、麦克风阵列等领域。其基本问题可以描述为：

利用一组空间分布的传感器（天线阵列）接收到的多路信号，估计信号的入射方向。

在实际应用中，DOA估计面临诸多挑战：

• 📉 低信噪比：信号功率微弱，淹没在噪声中
• 🔢 单快拍：仅有极少数测量样本可用
• 🔗 源相关性高：多个信号源具有高度相关性甚至相干性
• 📊 非均匀噪声：各阵元噪声功率不一致

传统方法（如 MUSIC、ESPRIT、MVDR）在这些极端条件下性能急剧退化。GRAIL-SCML 正是为了应对这些挑战而提出的。

---

2. 系统模型

2.1 均匀线阵（ULA）窄带模型

考虑一个由 $M$ 个阵元组成的均匀线性阵列（ULA），阵元间距为半波长 $d=\lambda /2$。假设有 $K$ 个远场窄带信号源，其入射角为 ${\theta }_1,{\theta }_2,\dots ,{\theta }_K$，则阵列接收数据可建模为：

$$\mathbf{Y}=\mathbf{A}(\mathbf{\theta })\mathbf{S}+\mathbf{N}$$

其中：

• $\mathbf{Y}\in {\mathbb{C}}^{M\times L}$：$L$ 个快拍的阵列接收数据矩阵
• $\mathbf{A}(\mathbf{\theta })\in {\mathbb{C}}^{M\times K}$：阵列流形矩阵（steering matrix）
• $\mathbf{S}\in {\mathbb{C}}^{K\times L}$：信号波形矩阵
• $\mathbf{N}\in {\mathbb{C}}^{M\times L}$：加性噪声矩阵

ULA的导向矢量定义为：

$$\mathbf{a}(\theta )=\frac{1}{\sqrt{M}}{\left[1,e^{j2\pi \frac{d}{\lambda }\sin \theta },\dots ,e^{j2\pi \frac{d}{\lambda }(M-1)\sin \theta }\right]}^{\top }$$

令 $d/\lambda =0.5$，导向矢量可简写为：

$$\mathbf{a}(\theta )=\frac{1}{\sqrt{M}}{\left[1,e^{j\pi \sin \theta },\dots ,e^{j\pi (M-1)\sin \theta }\right]}^{\top }$$

阵列流形矩阵为各方向导向矢量的拼接：

$$\mathbf{A}(\mathbf{\theta })=[\mathbf{a}({\theta }_1),\mathbf{a}({\theta }_2),\dots ,\mathbf{a}({\theta }_K)]$$

2.2 噪声模型

项目支持两种噪声模型：

• 均匀噪声（uniform）：各阵元噪声方差相同，${\sigma }_1^2={\sigma }_2^2=\cdots ={\sigma }_M^2$
• 非均匀噪声（nonuniform）：各阵元噪声方差不同，具有指数级差异（跨距可达 $12\text{dB}$）

2.3 源相关模型

信号源的相关性由相关系数 $\rho$ 刻画（以 $K=2$ 为例）：

$${\mathbf{R}}_s=\left[\begin{array}{cc}1& \rho \\ {\rho }^{\ast }& 1\end{array}\right]$$

当 $|\rho |\to 1$ 时，两个信号源高度相关（实验中取 $\rho =0.9\times (0.5010+j0.8654)$）。

2.4 宽带模型

对于宽带信号，采用频率聚焦（frequency focusing）框架。在 $\omega$-空间下定义归一化方向参数：

$$\omega =\frac{1}{2}\sin \theta$$

频率因子定义为：

$$f_m=1-\frac{M_f+1-2m}{M_0},m=1,2,\dots ,M_f$$

其中 $M_f$ 是频率子带数量，$M_0$ 是参考频率参数。第 $m$ 个频率子带的导向矢量为：

$${\mathbf{a}}_m(\omega )=\frac{1}{\sqrt{N}}{\left[1,e^{j2\pi f_m\omega },\dots ,e^{j2\pi f_m(N-1)\omega }\right]}^{\top }$$

---

3. GRAIL-SCML 算法详解

GRAIL-SCML 的全称揭示了其核心设计思想：

缩写	全称	含义
G	Geometry-Regularized	几何正则化
RA	Regularized Adaptive	正则化自适应
IL	Integrated Likelihood	集成似然
SCML	Structure-Constrained Maximum Likelihood	结构约束最大似然

整个算法由三个核心模块和两种关键技术组成。

3.1 波形域确定性最大似然（DML）拟合

给定角度 $\mathbf{\theta }$，信号波形的最小二乘估计为：

$$\hat{\mathbf{S}}=({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H\mathbf{Y}$$

因此确定性最大似然的代价函数为：

$$J_{\text{DML}}(\mathbf{\theta })=\frac{\Vert \mathbf{Y}-{\mathbf{P}}_A\mathbf{Y}{\Vert }_F^2}{\max (\Vert \mathbf{Y}{\Vert }_F^2,\epsilon )}$$

其中 ${\mathbf{P}}_A=\mathbf{A}({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H$ 是投影矩阵。该代价函数衡量阵列流形张成空间对接收数据的拟合程度——值越小，说明选定的角度越接近真实值。

直觉理解：DML 试图找到一组角度，使得信号在该方向的投影尽可能"解释"所有观测数据。

3.2 协方差域拟合

仅依赖波形域拟合在样本稀少时不够鲁棒。GRAIL-SCML 引入协方差域拟合作为互补信息源。

首先估计样本协方差矩阵：

$$\hat{\mathbf{R}}=\frac{1}{L}\mathbf{Y}{\mathbf{Y}}^H$$

然后应用前向-后向平均：

$${\hat{\mathbf{R}}}_{\text{FB}}=\frac{1}{2}\left(\hat{\mathbf{R}}+\mathbf{J}{\hat{\mathbf{R}}}^{\ast }\mathbf{J}\right)$$

其中 $\mathbf{J}$ 是反恒等矩阵（反对角线上的元素为 $1$）。接着进行收缩正则化：

$$\tilde{\mathbf{R}}=(1-\alpha ){\hat{\mathbf{R}}}_{\text{FB}}+\alpha \cdot \frac{\text{tr}({\hat{\mathbf{R}}}_{\text{FB}})}{M}\mathbf{I}$$

协方差域拟合的思想是：将估计的协方差矩阵分解为信号分量和噪声分量：

$$\mathbf{R}=\sum _{k=1}^K{p}_k\mathbf{a}({\theta }_k){\mathbf{a}}^H({\theta }_k)+p_{K+1}\mathbf{I}$$

通过求解如下非负最小二乘问题获得功率系数：

min ⁡ p ≥ 0 ∥ [ Re { B } Im { B } ] p − [ Re { r } Im { r } ] ∥ 2 2 \min_{\mathbf{p} \geq \mathbf{0}} \left\|\begin{bmatrix}\text{Re}\{\mathbf{B}\} \\ \text{Im}\{\mathbf{B}\}\end{bmatrix} \mathbf{p} - \begin{bmatrix}\text{Re}\{\mathbf{r}\} \\ \text{Im}\{\mathbf{r}\}\end{bmatrix}\right\|_2^2 p≥0min​ ​[Re{B}Im{B}​]p−[Re{r}Im{r}​] ​22​

其中 $\mathbf{r}=\text{vec}(\tilde{\mathbf{R}})$，$\mathbf{B}$ 的第 $k$ 列为 $\text{vec}(\mathbf{a}({\theta }_k){\mathbf{a}}^H({\theta }_k))$。协方差拟合代价为：

$$J_{\text{COV}}(\mathbf{\theta })=\frac{\Vert \tilde{\mathbf{R}}-{\hat{\mathbf{R}}}_{\text{fit}}{\Vert }_F^2}{\max (\Vert \tilde{\mathbf{R}}{\Vert }_F^2,\epsilon )}$$

为什么需要协方差域？ 波形域对信号波形的随机性敏感，而协方差域通过统计平均能有效压制噪声波动，在低快拍场景提供更稳定的梯度方向。

3.3 弱几何正则化

为了防止估计值坍缩（例如多个角度估计聚集到同一点），引入两项弱正则化：

（1）分离度惩罚项

$$J_{\text{sep}}(\mathbf{\theta })=\frac{1}{K-1}\sum _{k=2}^K\frac{1}{\max ({\theta }_k-{\theta }_{k-1}{)}^2,\epsilon }$$

当相邻角度过于接近时，该惩罚项急剧增大，推动角度分离。

（2）边界惩罚项

$$J_{\text{bound}}(\mathbf{\theta })=\sum _k\left[({\theta }_{\min }-{\theta }_k{)}^2\cdot {\mathbb{I}}_{{\theta }_k<{\theta }_{\min }}+({\theta }_k-{\theta }_{\max }{)}^2\cdot {\mathbb{I}}_{{\theta }_k>{\theta }_{\max }}\right]$$

防止搜索超出预设的角度范围。

关键设计：正则化权重极小（分离惩罚权重仅 $0.02$，边界惩罚权重 $10^3$ 但因触发概率低而影响有限），确保正则化在正常区域几乎无影响，仅在估计即将"出错"时提供纠正。这就是"弱正则化"的含义。

3.4 总代价函数

将上述三项组合为总目标函数：

$$\boxed{J(\mathbf{\theta })=w_{\text{DML}}\cdot J_{\text{DML}}(\mathbf{\theta })+w_{\text{COV}}\cdot J_{\text{COV}}(\mathbf{\theta })+w_{\text{sep}}\cdot J_{\text{sep}}(\mathbf{\theta })}$$

默认权重为：

• $w_{\text{DML}}=0.65$
• $w_{\text{COV}}=0.35$
• $w_{\text{sep}}=0.02$

权重设计哲学：DML 为主、COV 为辅（65:35），分离惩罚作为最后的"安全阀"。这一设计平衡了拟合精度和鲁棒性。

3.5 连续支撑集精炼

传统DOA方法通常在离散网格上搜索，精度受限于网格分辨率。GRAIL-SCML 采用无网格连续优化彻底消除这一限制。

使用 Nelder-Mead 单纯形算法（fminsearch）直接在连续空间中最小化 $J(\mathbf{\theta })$：

$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }\in [{\theta }_{\min },{\theta }_{\max }{]}^K}{\min }J(\mathbf{\theta })$$

搜索参数：

• 最大迭代次数：$70$
• 最大函数评估次数：$220$
• 角度容差：$10^{-4}$
• 函数容差：$10^{-5}$

3.6 证据支撑集选择：多初始化器策略

连续非凸优化对初始值高度敏感。为解决这一问题，GRAIL-SCML 采用多初始化器 + 择优策略。

初始化器池包括 $9$ 种算法（最多取前 $7$ 种有效结果）：

类别	初始化器	特点
子空间类	FB-MUSIC	前向-后向MUSIC
子空间类	GSB-DOA	广义子空间DOA
子空间类	Root-MUSIC	求根MUSIC
子空间类	TLS-ESPRIT	总体最小二乘ESPRIT
稀疏类	IAA-APES	迭代自适应方法
稀疏类	SPICE	稀疏迭代协方差估计
束波类	MVDR (Capon)	最小方差无畸变响应
束波类	CBF (Bartlett)	常规波束形成
贪心类	MMV-OMP	多测量向量正交匹配追踪
顺序类	Seq-DML	顺序确定性最大似然初始化

处理流程：

1. 候选生成：依次运行各初始化器，获取 $K$ 个角度估计
2. 候选清洗：剔除重复候选（角度差 $<0.25^{\circ }$ 视为重复）
3. 独立精炼：对每个候选，运行 Nelder-Mead 优化器
4. 证据选择：选择精炼后代价函数最小的结果作为最终输出

$${\hat{\mathbf{\theta }}}_{\text{final}}=\arg \underset{{\hat{\mathbf{\theta }}}_i}{\min }J({\hat{\mathbf{\theta }}}_i^{\text{refined}})$$

为什么叫"证据"选择？ 最终的判决依据是总代价函数 $J(\cdot )$，它集成了波形域和协方差域两方面的"证据"，选出最一致的解。

3.7 算法伪代码

算法: GRAIL-SCML (窄带)

输入: 阵列数据 Y, 源数量 K, 配置参数 cfg
输出: DOA估计 θ̂

1. 初始化候选集 C ← ∅
2. for 每个初始化器:
3.     θ_init ← initializer(Y, K, cfg)
4.     if θ_init 有效: C ← C ∪ {θ_init}
5. end for
6. C ← 去重(C, 容差=0.25°)
7.
8. bestCost ← +∞, bestθ ← C[1]
9. for each θ_cand in C (最多 maxInitializers 个):
10.    θ_refined ← NelderMead(J(·), θ_cand)
11.    cost ← J(θ_refined)
12.    if cost < bestCost:
13.        bestCost ← cost, bestθ ← θ_refined
14. end for
15.
16. 返回 sort(bestθ)

---

4. 宽带扩展：WB-GRAIL-SCML

宽带 GRAIL-SCML 将窄带框架扩展到多频带场景，命名为 WB-GRAIL-SCML。

4.1 核心差异

维度	窄带 GRAIL-SCML	宽带 WB-GRAIL-SCML
数据	单频阵列数据 $\mathbf{Y}$	多频阵列数据 { Y m } m = 1 M f \{\mathbf{Y}_m\}_{m=1}^{M_f} {Ym​}m=1Mf​​
导向矢量	$\mathbf{a}(\theta )$	${\mathbf{a}}_m(\omega )$
搜索空间	$\theta$ (角度/度)	$\omega$ (归一化频率)
DML代价	单频DML	各子带DML平均
COV代价	单频COV	各子带COV平均
初始化器	9种窄带方法	4种宽带方法 + Seq-WB

4.2 宽带代价函数

$${\bar{J}}_{\text{DML}}(\mathbf{\omega })=\frac{1}{M_f}\sum _{m=1}^{M_f}\frac{\Vert {\mathbf{Y}}_m-{\mathbf{P}}_{A_m}{\mathbf{Y}}_m{\Vert }_F^2}{\max (\Vert {\mathbf{Y}}_m{\Vert }_F^2,\epsilon )}$$

$${\bar{J}}_{\text{COV}}(\mathbf{\omega })=\frac{1}{M_f}\sum _{m=1}^{M_f}\frac{\Vert {\tilde{\mathbf{R}}}_m-{\hat{\mathbf{R}}}_{\text{fit},m}{\Vert }_F^2}{\max (\Vert {\tilde{\mathbf{R}}}_m{\Vert }_F^2,\epsilon )}$$

总代价为：

$$J_{\text{WB}}(\mathbf{\omega })=0.65{\bar{J}}_{\text{DML}}+0.35{\bar{J}}_{\text{COV}}+0.02J_{\text{sep}}$$

---

5. 实验设计与结果分析

项目提供了四个精心设计的实验场景，使用 $500$ 次蒙特卡洛实验评估 RMSE 性能。

实验一：非均匀噪声 + 有限快拍

参数	设置
阵元数 $M$	$8$
源数 $K$	$2$
快拍数 $L$	$40$
真实DOA	$[-2^{\circ },7^{\circ }]$
SNR范围	$-20\to 20\text{ dB}$
噪声模型	非均匀（跨距 $12\text{dB}$）
源相关性	不相关

测试目的：验证算法在非均匀噪声环境下的鲁棒性。非均匀噪声严重破坏子空间类方法的假设（噪声子空间结构），而 GRAIL-SCML 的协方差域拟合直接建模噪声分量的对角结构，天然适应此场景。

实验二：单快拍 + 紧密源

参数	设置
阵元数 $M$	$10$
源数 $K$	$2$
快拍数 $L$	$1$
真实DOA	$\theta =\arcsin ([-0.1,0.1])$
SNR范围	$-5\to 20\text{ dB}$
度量空间	$u=\sin \theta$（因角度极小）

测试目的：极限条件下的分辨能力。$L=1$ 意味着大多数统计方法（如MUSIC、MVDR）直接失效——它们依赖协方差矩阵的秩结构，而单快拍下协方差矩阵秩为 $1$。GRAIL-SCML 的 DML 项直接在波形域工作，不受此限制。

实验三：宽带分离度扫描

参数	设置
子带数 $M_f$	$11$
各带阵元 $N$	$11$
源数 $K$	$2$
快拍数 $L$	$100$
SNR	$0\text{ dB}$
$\omega$ 分离度	$0.02\to 0.50$

测试目的：考察宽带 GRAIL-SCML（WB-GRAIL-SCML）在源逐渐分离时的性能收敛特性，验证其"超分辨"能力——即使在 $\Delta \omega =0.02$ 的极小分离度下仍能分辨双源。

实验四：高度相关窄带源

参数	设置
阵元数 $M$	$6$
源数 $K$	$2$
快拍数 $L$	$500$
相关系数	$\rho =0.9(0.5010+j0.8654)$
SNR范围	$-25\to 20\text{ dB}$

测试目的：高相关性破坏信号子空间的秩结构，严重削弱 MUSIC 等子空间方法的性能。GRAIL-SCML 同时在波形域和协方差域评估，两个域各自提供独立证据，协同克服相关性干扰。

对比方法汇总

窄带实验（12种方法）：

• 🟦 经典类：CBF、MVDR
• 🟩 子空间类：FB-MUSIC、Root-MUSIC、TLS-ESPRIT、GSB-DOA
• 🟨 稀疏类：MMV-OMP、IAA-APES、SPICE、SLIM、SCMR-SBL-GR
• 🟥 本文方法：GRAIL-SCML

宽带实验（5种方法）：

• I-MUSIC、CSSM、WB-SBL、COFFEE-CF
• 本文方法：WB-GRAIL-SCML

---