Off-Grid Direction of Arrival Estimation Using Sparse Bayesian Inference (II)

IV. NUMERICAL SIMULATIONS

在本节中，我们展示了DOA估计的数值模拟结果。我们考虑一个由 $M=8$ 个传感器组成的标准均匀线性阵列（ULA）。原点设置在ULA的中点，以减小(2)中的近似误差。因此我们有 $A_{mn}=\exp \left\{j\pi (m-\frac{M+1}{2})\cos {\tilde{\theta }}_n\right\}$ 和 $B_{mn}=-j\pi (m-\frac{M+1}{2})\sin {\tilde{\theta }}_n\cdot A_{mn}$，其中 $m=1,\cdots ,M$，$n=1,\cdots ,N$，且 $j=\sqrt{-1}$。我们考虑一个均匀采样网格 {0∘,r,2r,⋯ ,180∘}\{0^\circ, r, 2r, \cdots, 180^\circ\}{0∘,r,2r,⋯,180∘}，其中 $r$ 是网格间隔。在MMV情况下，快照数设置为 $T=200$。在MMV情况下，我们只考虑OGSBI-SVD，因为根据经验观察，与OGSBI相比，它收敛更快、更准确。

在OGSBI-SVD中，我们设置 $\rho =0.01$ 和 $c=d=1\times 10^{-4}$。我们初始化 α0=100K∑t=1KVar{YSV}t\alpha_0 = \frac{100K}{\sum_{t=1}^K \text{Var}\{\mathbf{Y}_{SV}\}_t}α0​=∑t=1K​Var{YSV​}t​100K​，$\mathbf{\alpha }=\frac{1}{MK}{\sum }_{t=1}^K|{\mathbf{A}}^H({\mathbf{Y}}_{SV}{)}_t|$ 和 $\mathbf{\beta }=0$，其中 $|\cdot |$ 逐元素作用。我们设置 $\tau =10^{-3}$ 并且最大迭代次数为1000。

我们注意到，所提出的算法对 ${\alpha }_0$、$\mathbf{\alpha }$ 和 $\mathbf{\beta }$ 的初始化不敏感，如果 $\rho$ 不是太大，对 $\rho$ 也不敏感。正如[12]中所报道的，在某些情况下，对 ${\alpha }_0$ 的估计可能不准确。但我们观察到这对DOA估计的结果影响极小。详细的讨论请读者参考[18]。所有实验都在一台配有Windows XP系统和3 GHz CPU的个人电脑上，使用Matlab v.7.7.0进行。Matlab代码已在线提供，网址为 https://sites.google.com/site/zaiyang0248/publication 。

A. Comparison With $l_1$-SVD

我们将 [5] 中的 $l_1$-SVD 作为基于在网格模型的代表性方法，并就均方误差（MSE）和计算时间方面，在不同网格间隔 $r$ 和信噪比 SNR 下将其与 OGSBI-SVD 进行比较。

在我们的实验中，我们考虑 SNR = 10 和 0 dB，以及 $r$ = 0.5°、1°、2° 和 4°。在每次试验中，我们在 [58°, 62°] 和 [86°, 90°] 的区间内分别均匀地生成 $K=2$ 个信源 ${\theta }_1,{\theta }_2$。在展示我们的比较结果之前，我们使用柯尔莫哥洛夫-斯米尔诺夫检验（Kolmogorov-Smirnov test）表明，在所有场景中，(4)中的总噪声（测量噪声加近似误差）以至少94%的概率服从高斯分布。这从经验上验证了离网格模型。对于每种（SNR, r）组合，MSE 在 $R=200$ 次试验中取平均：$MSE=\frac{1}{RK}{\sum }_{i=1}^R{\sum }_{k=1}^K({\theta }_k^i-{\hat{\theta }}_k^i{)}^2$，其中上标 $i$ 指的是第 $i$ 次试验。

应当注意，无论信噪比如何，$l_1$-SVD 的 MSE 存在一个下界，因为 $l_1$-SVD 能得到的最佳DOA估计是距离真实DOA最近的网格点。事实上，这个下界是所有基于在网格模型的方法所共有的，包括CS-MUSIC、SA-MUSIC和[13]中的算法。
通过假设真实DOA是均匀分布的，该下界可以很容易地计算出为 $LB=\frac{r^2}{12}$。³

图1展示了我们的实验结果。在所有考虑的场景中，OGSBI-SVD 都比 $l_1$-SVD 具有更准确的DOA估计。此外，在大多数场景中，OGSBI-SVD 可以超越 $l_1$-SVD 的下界。在较高信噪比或较粗采样网格的情况下，这种现象尤为显著，因为在这些情况下，在网格模型描述真实观测模型的性能较差，而本文中使用的离网格模型在很大程度上克服了这个问题。

表 I 展示了 OGSBI-SVD 和 $l_1$-SVD 在不同信噪比 SNR 和网格间隔 $r$ 下的平均 CPU 时间（在我们的案例中，不包括耗时约 0.003s 的 SVD 过程）。⁴ 对于 OGSBI-SVD 和 $l_1$-SVD 两种算法，当网格变粗时，它们的 CPU 时间都会减少。在 $r=2^{\circ }$ 和 $4^{\circ }$ 时，OGSBI-SVD 比 $l_1$-SVD 更快。所提出方法的一个缺点是在密集采样网格的情况下速度较慢。在实践中，我们建议为所提出的算法使用 $r=2^{\circ }$ 的较粗网格，因为它能给出既准确又快速的DOA估计。

注 2：

1. 我们选择 $l_1$ 优化进行比较，因为它通常被认为在稀疏恢复方面有更好的理论保证，尽管在我们的设置中（两个信源被很好地分离开）一些更简单的CS求解器，例如OMP [22]，也可能成功，并且在这类低维问题中计算成本更低。对于间隔很近的信源的情况，读者可以参考 [18] 获取更多关于 $l_1$ 优化和所提方法的仿真结果。
2. 所提出算法的另一个优点是，与 $l_1$-SVD 相比，它具有更小的DOA估计偏差 [18]。

B. Comparison With STLS

离网格模型最近已在 [17] 中被用于DOA估计。在 [17] 中，提出了一种稀疏总最小二乘（sparse total least-squares，STLS）方法。在SMV情况下，STLS旨在解决以下非凸优化问题：

min⁡x,β{∣∣β∣∣22+∣∣y−(A+Bdiag(β))x∣∣22+λ∣∣x∣∣1}(27) \min_{\mathbf{x},\boldsymbol{\beta}} \{ ||\boldsymbol{\beta}||_2^2 + ||\mathbf{y} - (\mathbf{A}+\mathbf{B}\text{diag}(\boldsymbol{\beta}))\mathbf{x}||_2^2 + \lambda ||\mathbf{x}||_1 \} \tag{27} x,βmin​{∣∣β∣∣22​+∣∣y−(A+Bdiag(β))x∣∣22​+λ∣∣x∣∣1​}(27)

其中 $\mathbf{x}$ 是感兴趣的稀疏源信号，$\mathbf{y}$ 是带噪声的测量值，$\mathbf{A}$、$\mathbf{B}$ 和 $\mathbf{\beta }$ 与离网格模型中的定义相同，并且 $\lambda >0$ 是一个正则化参数。

从贝叶斯角度来看，这等同于在假设测量噪声为白高斯噪声、$\mathbf{x}$ 服从拉普拉斯分布、$\mathbf{\beta }$ 服从高斯分布的条件下，寻求 ($\mathbf{x}$, $\mathbf{\beta }$) 的MAP解。值得注意的是，最后这个关于 $\mathbf{\beta }$ 的高斯假设无法恰当地捕捉 $\mathbf{\beta }$ 的属性。

在 [17] 中，(27) 中的问题是通过一种交替方法来求局部最小值的，即，在固定 $\mathbf{\beta }$ 的情况下交替求解 $\mathbf{x}$（这需要解一个 $l_1$-正则化最小二乘问题），以及在固定 $\mathbf{x}$ 的情况下求解 $\mathbf{\beta }$（这需要解一个 $N$ 维线性系统）。

如 [5] 中所论证的，OGSBI-SVD中使用的SVD可以减轻MMV情况下对测量噪声的敏感性，而 [17] 中并未使用SVD。为了进行公平比较，在将我们的方法与STLS比较时，我们只考虑SMV情况，尽管对于MMV情况的STLS也可以构建类似的问题。在我们的OGSBI实现中，我们初始化 α0=100Var{y}\alpha_0 = \frac{100}{Var\{y\}}α0​=Var{y}100​ 和 $\mathbf{\alpha }=\frac{1}{MK}|{\mathbf{A}}^{H}y|$。其余设置与MMV情况下的OGSBI-SVD相同。

在我们的实验中，我们考虑两个DOA，分别为 $63.2^{\circ }$ 和 $90.3^{\circ }$，信噪比 SNR = 20 dB。我们为OGSBI和STLS都考虑了 $r=2^{\circ }$ 和 $4^{\circ }$ 的情况。(27)中的参数 $\lambda$ 被我们调整到最佳，以使STLS达到最小的误差。表II展示了STLS和OGSBI在 $R=200$ 次试验中的平均MSE和CPU时间。在两种场景中，OGSBI都获得了比STLS更准确的DOA估计，并且计算时间显著减少。我们还注意到，使用最先进的CS算法来加速STLS是可能的。

1. 两种等价视角的定义
   问题可通过两种核心视角描述，二者存在内在等价性：
   + 优化视角：核心为最小化成本函数（公式27），表达式为：
   min⁡x,β{∣∣β∣∣22+∣∣y−(A+Bdiag(β))x∣∣22+λ∣∣x∣∣1} \min_{\mathbf{x},\boldsymbol{\beta}} \{ ||\boldsymbol{\beta}||_2^2 + ||\mathbf{y} - (\mathbf{A}+\mathbf{B}\text{diag}(\boldsymbol{\beta}))\mathbf{x}||_2^2 + \lambda ||\mathbf{x}||_1 \} x,βmin​{∣∣β∣∣22​+∣∣y−(A+Bdiag(β))x∣∣22​+λ∣∣x∣∣1​}
   - 贝叶斯视角：核心为最大后验概率（MAP）估计，即寻找给定数据$\mathbf{y}$时，参数$(\mathbf{x},\mathbf{\beta })$最可能的取值。
2. 核心关联：从后验概率到负对数形式的转化
   根据贝叶斯定理，后验概率$p(\mathbf{x},\mathbf{\beta }|\mathbf{y})\propto p(\mathbf{y}|\mathbf{x},\mathbf{\beta })\cdot p(\mathbf{x})\cdot p(\mathbf{\beta })$。
   最大化后验概率等价于最大化其对数形式，而对数的单调性与加法转化特性，进一步等价于最小化其负对数形式：
   $$\underset{\mathbf{x},\mathbf{\beta }}{\min }\left[-\log p(\mathbf{y}|\mathbf{x},\mathbf{\beta })-\log p(\mathbf{x})-\log p(\mathbf{\beta })\right]$$
3. 概率假设与成本函数项的一一对应
   通过三个关键概率假设，可建立贝叶斯视角与优化视角的具体对应关系：
   • 高斯噪声假设 $⟹$ 数据拟合项
     假设测量噪声服从高斯分布，则负对数似然$-\log p(\mathbf{y}|\mathbf{x},\mathbf{\beta })$正比于误差的二范数平方，对应成本函数中的数据拟合项：
     $$-\log p(\mathbf{y}|\mathbf{x},\mathbf{\beta })\propto ||\mathbf{y}-(\mathbf{A}+\mathbf{B}\text{diag}(\mathbf{\beta }))\mathbf{x}|{|}_2^2$$
   • $\mathbf{x}$ 的拉普拉斯先验假设 $⟹$ 稀疏正则项
     假设稀疏信号$\mathbf{x}$服从拉普拉斯先验，则负对数先验概率$-\log p(\mathbf{x})$正比于$\mathbf{x}$的一范数，对应成本函数中的稀疏正则项：
     $$-\log p(\mathbf{x})\propto \lambda ||\mathbf{x}|{|}_1$$
   • $\mathbf{\beta }$ 的高斯先验假设 $⟹$ 偏移量正则项
     假设偏移量 $\mathbf{\beta }$ 服从高斯先验，则负对数先验概率$-\log p(\mathbf{\beta })$正比于$\mathbf{\beta }$的二范数平方，对应成本函数中的偏移量正则项：
     $$-\log p(\mathbf{\beta })\propto ||\mathbf{\beta }|{|}_2^2$$
4. 结论：两种视角的数学等价性
   将贝叶斯视角下的三个负对数项相加，可完整重构优化视角下的成本函数：
   min⁡x,β{∣∣β∣∣22+∣∣y−(A+Bdiag(β))x∣∣22+λ∣∣x∣∣1}=min⁡x,β[−log⁡p(y∣x,β)−log⁡p(x)−log⁡p(β)] \min_{\mathbf{x},\boldsymbol{\beta}} \{ ||\boldsymbol{\beta}||_2^2 + ||\mathbf{y} - (\mathbf{A}+\mathbf{B}\text{diag}(\boldsymbol{\beta}))\mathbf{x}||_2^2 + \lambda ||\mathbf{x}||_1 \} = \min_{\mathbf{x},\boldsymbol{\beta}} \left[ -\log p(\mathbf{y} | \mathbf{x}, \boldsymbol{\beta}) - \log p(\mathbf{x}) - \log p(\boldsymbol{\beta}) \right] x,βmin​{∣∣β∣∣22​+∣∣y−(A+Bdiag(β))x∣∣22​+λ∣∣x∣∣1​}=x,βmin​[−logp(y∣x,β)−logp(x)−logp(β)]
   因此，在文中三个概率假设下，求解STLS优化问题与寻求MAP解在数学上完全等价。

V. CONCLUSION

在本文中，我们研究了[17]中首次提出的离网格DOA估计模型，该模型用于减少因连续范围离散化而引起的建模误差。我们从贝叶斯角度提出了一种基于离网格模型的算法，该算法同时适用于单快照和多快照情况。我们利用一种基于子空间的想法来降低信号恢复过程的计算复杂度和对噪声的敏感性。

我们通过仿真说明，所提出的方法优于标准的CS方法，因为后者的性能受到其内在的标准在网格模型（on-grid model）的限制。它也比[17]中基于离网格模型的算法更准确。

所提出算法的一个缺点是在密集采样网格的情况下速度较慢，不过可以采用更粗糙的网格来获得准确而快速的DOA估计。未来的一个工作是开发我们算法的快速版本。在这项工作之后，我们在[24]中表明，$l_1$ 优化也适用于离网格DOA估计问题，并且在某些条件下也提供了性能保证。

VI. 收敛性证明

1. 第一步：建立核心等式
   证明始于一个巧妙的恒等变换。对于任意的参数 $\theta$ 和隐变量 $Z$，我们有：
   $$\log P(X|\theta )=\log P(X,Z|\theta )-\log P(Z|X,\theta )$$
   这个式子只是条件概率公式 $P(X,Z|\theta )=P(Z|X,\theta )P(X|\theta )$ 的简单变形。
   接下来，我们在等式两边同时取关于后验分布 $P(Z|X,{\theta }^{(t)})$ 的期望（这里的 ${\theta }^{(t)}$ 是我们上一步迭代得到的参数）。

• 等式左边: 因为 $\log P(X|\theta )$ 与隐变量 $Z$ 无关，所以对它求期望，结果还是它本身。
• 等式右边: 展开为两项期望的差。
  这样我们就得到了一个关键的分解式：
  $$\log P(X|\theta )=\underset{\text{定义为 }Q(\theta ,{\theta }^{(t)})}{\underbrace{{\mathbb{E}}_{Z|X,{\theta }^{(t)}}[\log P(X,Z|\theta )]}}-\underset{\text{我们暂时称之为 }H(\theta ,{\theta }^{(t)})}{\underbrace{{\mathbb{E}}_{Z|X,{\theta }^{(t)}}[\log P(Z|X,\theta )]}}$$
  其中，$Q(\theta ,{\theta }^{(t)})$ 就是我们在E-步中计算和在M-步中要最大化的Q函数。

2. 第二步：分析两次迭代的似然函数之差
   为了证明似然函数是不减的，我们来考察两次迭代之间的差值：
   $$\log P(X|{\theta }^{(t+1)})-\log P(X|{\theta }^{(t)})$$
   利用上面的核心等式，这个差值可以写为：
   $$[Q({\theta }^{(t+1)},{\theta }^{(t)})-Q({\theta }^{(t)},{\theta }^{(t)})]-[H({\theta }^{(t+1)},{\theta }^{(t)})-H({\theta }^{(t)},{\theta }^{(t)})]$$
   现在，我们只需要分别证明第一部分 [...] 是非负的，第二部分 [...] 是非正的，就能完成证明。

3. 第三步：分析第一部分 (Q函数的变化)
   根据M-步的定义，我们选择的 ${\theta }^{(t+1)}$ 正是使 $Q(\theta ,{\theta }^{(t)})$ 最大的那个 $\theta$：
   $${\theta }^{(t+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(t)})$$
   这意味着，新的 Q 函数值一定大于等于旧的 Q 函数值：
   $$Q({\theta }^{(t+1)},{\theta }^{(t)})\ge Q({\theta }^{(t)},{\theta }^{(t)})$$
   所以，第一部分 $[Q(\dots )-Q(\dots )]$ 必然大于等于0。

4. 第四步：分析第二部分 (KL散度的出现)
   现在我们来看第二部分 $[H(\dots )-H(\dots )]$。它的展开形式是：
   $${\mathbb{E}}_{Z|X,{\theta }^{(t)}}[\log P(Z|X,{\theta }^{(t+1)})]-{\mathbb{E}}_{Z|X,{\theta }^{(t)}}[\log P(Z|X,{\theta }^{(t)})]$$
   利用对数的性质 $\log a-\log b=\log (a/b)$，上式可以合并为：
   $${\mathbb{E}}_{Z|X,{\theta }^{(t)}}\left[\log \frac{P(Z|X,{\theta }^{(t+1)})}{P(Z|X,{\theta }^{(t)})}\right]$$
   这正是负的KL散度。KL散度的定义是 $D_{KL}(P||Q)={\mathbb{E}}_P[\log \frac{P}{Q}]$。
   所以，上式等于：
   $$-D_{KL}\left(P(Z|X,{\theta }^{(t)})||P(Z|X,{\theta }^{(t+1)})\right)$$
   根据相对熵（KL散度）的性质，KL散度永远大于等于0 ($D_{KL}\ge 0$)。因此，负的KL散度必然小于等于0。

5. 第五步：得出结论
   我们回到第二步的差值表达式：
   $$\log P(X|{\theta }^{(t+1)})-\log P(X|{\theta }^{(t)})=\underset{\ge 0}{\underbrace{[Q(\dots )-Q(\dots )]}}-\underset{\le 0}{\underbrace{[H(\dots )-H(\dots )]}}$$
   一个非负数减去一个非正数，结果必然是非负数。
   因此，我们证明了：
   $$\log P(X|{\theta }^{(t+1)})-\log P(X|{\theta }^{(t)})\ge 0$$
   即对数似然函数在EM算法的每次迭代中都是单调不减的，从而保证了算法的收敛性。

VII. Code

1. 初始化阶段

• 从输入的 paras 结构体中提取测量数据 Y、字典矩阵 A、导数矩阵 B 及其他算法控制参数
• 将输入的初始噪声方差 sigma2 转换为噪声精度 alpha0
• 根据是否已知噪声方差，设置 alpha0 的伽马先验分布参数 a 和 b
• 确定信源数量 K，并为提高计算效率预先计算 B' * B

2. EM 迭代主循环
   算法进入主循环，交替执行 E-步 和 M-步 直至收敛

   • 2.1 E-步 (Expectation Step): 计算信号后验分布

     • 根据当前估计的离网格参数 beta，构建完整的观测矩阵 Phi
     • 应用 Woodbury矩阵恒等式 来高效地计算信号的后验协方差矩阵 Sigma (Σ)
     • 根据公式(15)，计算信号的后验均值矩阵 mu (μ)

   • 2.2 M-步 (Maximization Step): 更新超参数
     利用 E-步 的计算结果来更新所有模型的超参数

     • 2.2.1 更新 alpha (α)
       根据公式(17)，结合后验均值和方差更新稀疏度参数 alpha

     • 2.2.2 更新 alpha0 (α₀)
       根据公式(18)，结合当前的残差等信息更新噪声精度 alpha0

     • 2.2.3 更新 beta (β)

       1. 确定活动集: 算法从 alpha 向量中找到最大的 K 个值，将这些值对应的索引作为当前迭代的"活动集"
       2. 计算P和v: 仅针对活动集，根据公式(20)和(21) 计算用于 beta 优化的二次型矩阵 P 和向量 v
       3. 求解beta: 首先尝试直接计算全局最优解 beta = P⁻¹v。若该解不满足约束，则切换到坐标下降法，根据公式(25)和(26) 对活动集中的 beta 元素进行逐个迭代更新

3. 收敛判断与终止

• 在每次迭代后，计算 alpha 向量的相对变化量，并与设定的容差 tol 比较
• 如果变化量小于容差，或达到最大迭代次数 maxiter，则算法终止
• 代码还计算了边缘似然函数值 ML，用于监控算法的收敛性（该值应单调递增）

5. 输出结果

• 循环结束后，函数将最终的后验均值 mu、协方差 Sigma、更新后的参数 alpha 和 beta 等所有结果打包并返回
• 最终的DOA估计值可由 alpha 峰值对应的网格点角度加上相应的 beta 偏移量得到

实际上TSP文章没有给出EM算法收敛性的证明过程