基于确定性最大似然算法 DML 的 DoA 估计，用牛顿法实现（附 MATLAB 源码）

本文首次在公众号【零妖阁】上发表，为了方便阅读和分享，我们将在其他平台进行自动同步。由于不同平台的排版格式可能存在差异，为了避免影响阅读体验，建议如有排版问题，可前往公众号查看原文。感谢您的阅读和支持！

在 DoA 估计中，最大似然方法主要分为确定性最大似然（DML）和随机性最大似然（SML）。当源信号是确定性信号时，为确定性最大似然法；当源信号为已知分布的随机信号时，为随机性最大似然法。

下面，我们要用确定性最大似然算法来估计目标的方位。

信号模型

假设空间中存在 $M$ 个不同方向的信号，入射到由 $N$ 个天线单元构成的均匀直线阵上。

令第 $m$ 个信号源的方向为 ${\theta }_m$，对应的信号波形为 $s_m(t)$。令第 $n$ 个天线单元的噪声为 $n_n(t)$。那么，在窄带远场条件下，天线阵的接收信号为
$$\mathbf{x}(t)=\mathbf{As}(t)+\mathbf{n}(t)$$

其中，$\mathbf{A}$ 为阵列流型矩阵，矩阵规模为 $N\times M$，具体可表示为 $M$ 个不同方向对应的阵列导向矢量：
$$\mathbf{A}=[\mathbf{a}({\theta }_1),\mathbf{a}({\theta }_2),\dots ,\mathbf{a}({\theta }_M)]$$

信号统计特征

假设信号模型中的噪声 $\mathbf{n}(t)$ 为圆对称高斯白噪声随机过程，不同阵元的噪声相互独立，信号波形 $\mathbf{s}(t)$ 为确定性信号。

在上述统计假设下，噪声 $\mathbf{n}(t)$ 的一阶矩和二阶矩满足

E { n ( t ) } = 0 \mathrm{E}\{\mathbf{n}(t)\} = 0 E{n(t)}=0

E { n ( t ) n ( t ) H } = σ 2 I \mathrm{E}\{ \mathbf{n}(t) \mathbf{n}(t)^H \} = \sigma^2 \mathbf I E{n(t)n(t)H}=σ2I

由于源信号是确定性信号，接收信号也服从高斯分布，其一阶矩和二阶矩满足

E { x ( t ) } = A s ( t ) \mathrm{E}\{\mathbf{x}(t)\} = \mathbf{A} \mathbf{s}(t) E{x(t)}=As(t)

E { x ( t ) x ( t ) H } = A s ( t ) + σ 2 I \mathrm{E}\{ \mathbf{x}(t) \mathbf{x}(t)^H \} = \mathbf{A} \mathbf{s}(t) + \sigma^2 \mathbf I E{x(t)x(t)H}=As(t)+σ2I

确定性最大似然算法 DML

根据文献 [3] 中圆对称高斯随机向量的概率密度公式，得到包含 $L$ 次快拍数据的接收信号的似然函数为

对数似然函数是关于 $\theta$ 、 ${\sigma }^2$ 和 $\mathbf{s}(t)$ 三个变量的函数。为了得到参数的最大似然估计，就需要求解对数似然函数取最小值时对应的未知量。

根据文献 [4] 中的推导，可以得到未知量 ${\sigma }^2$ 和 $\mathbf{s}(t)$ 的最大似然估计为

σ ^ 2 = 1 N t r { P A ⊥ R ^ } \hat \sigma^2 = \frac{1}{N} \mathrm{tr} \{ { \mathbf{P}_A^\perp \mathbf{\hat R} } \} σ^2=N1​tr{PA⊥​R^}

$$\hat{\mathbf{s}}(t)={\mathbf{A}}^{+}\mathbf{x}(t)$$

上式中，$\hat{\mathbf{R}}$ 为样本协方差矩阵， ${\mathbf{A}}^{+}$ 为 $\mathbf{A}$ 的伪逆矩阵，${\mathbf{P}}_A^{\perp }$ 为 $\mathbf{A}$ 的正交投影矩阵

$$\hat{\mathbf{R}}=\frac{1}{L}\sum _{l=1}^L\mathbf{x}(t_l){\mathbf{x}}^H(t_l)$$
$${\mathbf{A}}^{+}=({\mathbf{A}}^H\mathbf{A}{)}^{-1}{\mathbf{A}}^H$$
$${\mathbf{P}}_A^{\perp }=\mathbf{I}-{\mathbf{P}}_A=\mathbf{I}-\mathbf{A}{\mathbf{A}}^{+}$$

因此，最大似然可进一步简化为对 $\theta$ 的估计：
θ ^ = arg ⁡ min ⁡ Θ t r { P A ⊥ R ^ } \hat \theta = \arg \min_\Theta \mathrm{tr} \{ { \mathbf{P}_A^\perp \mathbf{\hat R} } \} θ^=argΘmin​tr{PA⊥​R^}

为了计算上述的最大似然估计，必须求解非线性多维优化问题。现在有很多算法可以解决这类优化问题，但是计算量都比较大，比如交替投影（AP）算法、MODE、迭代二次型最大似然（IQML）法、牛顿法、遗传算法等。

最大似然算法的优势在于其对相干信号也有效。

牛顿法寻优

本文选择牛顿法来搜索最大似然估计问题的最优解。

牛顿法是求解无约束优化问题最早使用的经典算法之一。其基本思想是用迭代点处的一阶导数（梯度）和二阶导数（Hesse 矩阵）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小点。

这里，将 $\psi = kd\sin \theta $ 作为最大似然估计问题的待优化变量。

目标函数为
f ( ψ ) = t r { P a ⊥ R ^ } = t r { [ I − A ( A H A ) − 1 A H ] R ^ } \begin{aligned} f(\psi) &= \mathrm{tr} \{ { \mathbf{P}_a^\perp \mathbf{\hat R} } \} \\ \\ &= \mathrm{tr} \{ { [\mathbf{I} - \mathbf{A}(\mathbf{A}^H\mathbf{A})^{-1} \mathbf{A}^H ] \mathbf{\hat R} } \} \\ \end{aligned} f(ψ)​=tr{Pa⊥​R^}=tr{[I−A(AHA)−1AH]R^}​

通过推导，得到目标函数的一阶导数为
F ( ψ ) = ∇ f ( ψ ) = − 2 ℜ { d i a g [ A + R ^ P A ⊥ D ] } \begin{aligned} F(\psi) &= \nabla f(\psi) \\ \\ &= -2 \Re \{ \mathrm{diag} [\mathbf{A}^+ \mathbf{\hat R} \mathbf{P}_A^\perp \mathbf{D} ] \} \end{aligned} F(ψ)​=∇f(ψ)=−2ℜ{diag[A+R^PA⊥​D]}​

目标函数的 Hessen 矩阵为
$$\begin{array}{rl}H(\psi )& ={\nabla }^{2}f(\psi )\\ & \\ & =2\Re \{({\mathbf{D}}^H{\mathbf{P}}_A^{\perp }\mathbf{D})\otimes ({\mathbf{A}}^{+}\hat{\mathbf{R}}{\mathbf{A}}^{+H}{)}^T\}\end{array}$$

上式中，$\mathbf{D}$ 为流行矩阵 $\mathbf{A}$ 相对于 $\psi$ 的一阶导矩阵。

牛顿法的性能对初始值的选取比较敏感，如果有一个很好的初始值，用牛顿法就能迅速收敛到极小值；如果初始值比较差，那么可能会收敛到局部极小值。

因此，在牛顿法寻优之前，需要用一些方法来确定迭代初始值，如谱搜索方法、交替投影法等。

MATLAB 仿真实验

为了提高阅读质量，MATLAB 源码在另一篇文章中给出：【代码纯享】DoA 估计 | 基于牛顿法的确定性最大似然算法 DML 的原理与实现（附 MATLAB 源码）

首先，直接根据最大似然谱函数绘制出空间谱分布，如下图

用这种谱搜索方法可以大致确定目标方位角的迭代初始值，但这种方法只适用于单个目标或者多个间隔一定距离的目标。对于多个相隔较近的目标，需要用交替迭代法确定初始值。

后面，会另写文章分享如何用交替投影法、遗传算法等实现最大似然估计，可以更高效更精确地识别多个独立或者相关目标。

参考文献

[1] 王永良. 空间谱估计理论与算法[M]. 清华大学出版社, 2004.
[2] 张小飞, 陈华伟, 仇小锋. 阵列信号处理及MATLAB实现[M]. 电子工业出版社, 2015.
[3] https://www.ee.iitb.ac.in/~sarva/courses/EE703/2013/Slides/CircularlySymmetricGaussian.pdf
[4] Bohme J . Estimation of source parameters by maximum likelihood and nonlinear regression[C]// ICASSP '84. IEEE International Conference on Acoustics, Speech, and Signal Processing. IEEE, 1984.
[5] https://sites.google.com/site/touchesofbreeze/sensor-array-processing-matlab-source-code文章来源地址https://www.toymoban.com/news/detail-583575.html

到了这里，关于基于确定性最大似然算法 DML 的 DoA 估计，用牛顿法实现（附 MATLAB 源码）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！