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

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

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

信号模型

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

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

其中， A \mathbf{A} A 为阵列流型矩阵，矩阵规模为 N × M N\times M N×M，具体可表示为 M M M 个不同方向对应的阵列导向矢量：
A = [ a ( θ 1 ) , a ( θ 2 ) , … , a ( θ M ) ] \mathbf{A} = [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \dots, \mathbf{a}(\theta_M) ] A=[a(θ1),a(θ2),…,a(θM)]

信号统计特征

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

在上述统计假设下，噪声 n ( t ) \mathbf{n}(t) 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 L L 次快拍数据的接收信号的似然函数为

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

根据文献 [4] 中的推导，可以得到未知量 σ 2 \sigma^2 σ2 和 s ( t ) \mathbf{s}(t) 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=N1tr{PA⊥R^}

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

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

R ^ = 1 L ∑ l = 1 L x ( t l ) x H ( t l ) \mathbf{\hat R} = \frac{1}{L} \sum_{l=1}^{L} \mathbf{x}(t_l)\mathbf{x}^H(t_l) R^=L1l=1∑Lx(tl)xH(tl)
A + = ( A H A ) − 1 A H \mathbf{A}^+ = (\mathbf{A}^H\mathbf{A})^{-1} \mathbf{A}^H A+=(AHA)−1AH
P A ⊥ = I − P A = I − A A + \mathbf{P}_A^\perp = \mathbf{I} - \mathbf{P}_A = \mathbf{I} - \mathbf{A} \mathbf{A}^+ PA⊥=I−PA=I−AA+

因此，最大似然可进一步简化为对 θ \theta θ 的估计：
θ ^ = arg ⁡ min ⁡ Θ t r { P A ⊥ R ^ } \hat \theta = \arg \min_\Theta \mathrm{tr} \{ { \mathbf{P}_A^\perp \mathbf{\hat R} } \} θ^=argΘmintr{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 矩阵为
H ( ψ ) = ∇ 2 f ( ψ ) = 2 ℜ { ( D H P A ⊥ D ) ⊗ ( A + R ^ A + H ) T } \begin{aligned} H(\psi) &= \nabla^2 f(\psi) \\ \\ &= 2 \Re \{ (\mathbf{D}^H \mathbf{P}_A^\perp \mathbf{D} ) \otimes (\mathbf{A}^+ \mathbf{\hat R} \mathbf{A}^{+H})^T \} \end{aligned} H(ψ)=∇2f(ψ)=2ℜ{(DHPA⊥D)⊗(A+R^A+H)T}

上式中， D \mathbf{D} D 为流行矩阵 A \mathbf{A} 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

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

【机器人栅格地图】基于蚁群优化遗传算法求解机器人栅格地图最短路径规划问题附Matlab源码
1 简介通过栅格法建立栅格地图作为机器人路径规划的工作环境,采用蚁群算法结合遗传算法作为机器人路径搜索的规则.将所有机器人放置于初始位置,经过NC次无碰撞迭代运动找到最优路径,到达目标位置.为防止机 ...
附Matlab源码：滚动轴承早期故障优化自适应蝙蝠算法优化的随机共振诊断法
附Matlab源码:滚动轴承早期故障优化自适应蝙蝠算法优化的随机共振诊断法文章目录附Matlab源码:滚动轴承早期故障优化自适应蝙蝠算法优化的随机共振诊断法前言一.自适应蝙蝠算法优化的随机共振 ...
【多目标优化求解】基于matlab粘菌算法MOSMA求解多目标优化问题【含Matlab源码 2279期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[多目标优化求解]基于matlab粘菌算法MOSMA求解多目标优化问题[含Matlab源码 2279期] 点击上面蓝色字体,直接付费下载, ...
【风电功率预测】基于matlab帝国殖民竞争算法优化BP神经网络风电功率预测【含Matlab源码 1314期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源: [风电功率预测]基于matlab帝国殖民竞争算法优化BP神经网络风电功率预测[含Matlab源码 1314期] ⛄二.帝国殖民竞争算法简 ...
【路径规划】基于matlab DWA算法机器人局部避障路径规划【含Matlab源码 890期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[路径规划]基于matlab DWA算法机器人局部避障路径规划[含Matlab源码 890期] 获取代码方式2: 通过订阅紫极神光博客付费 ...
【多目标优化求解】基于matlab飞蛾扑火算法 (NSMFO)求解多目标优化问题【含Matlab源码 2312期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[多目标优化求解]基于matlab飞蛾扑火算法 (NSMFO)求解多目标优化问题 [含Matlab源码 2312期] 点击上面蓝色字体,直 ...
初勘数据处理之预测算法（随机森林random forest）附matlab源码
初勘数据处理之预测算法(随机森林random forest)附matlab源码集成学习方法集成学习方法,是整合多个分类方法的预测结果,以提高分类准确率.集成学习方法先用训练数据构建一组基分类器,再 ...
【数字信号去噪】基于matlab稀疏性BEADS色谱基线估计和去噪【含Matlab源码 1887期】
一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[数字信号去噪]基于matlab稀疏性 BEADS色谱基线估计和去噪[含Matlab源码 1887期] 点击上面蓝色字体,直接付费下载,即可 ...
基于ITD实现的轴承故障信号分解并附Matlab源码
基于ITD实现的轴承故障信号分解并附Matlab源码轴承是旋转机械中重要的支撑部件,其状态的健康程度对设备运行的性能和寿命有着至关重要的影响.因此,轴承状态监测和故障诊断成为了研究热点.本文基于IT ...

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