基于MATLAB的均匀面阵MUSIC算法DOA估计仿真

---

前言

\;\;\;\;\; 在波达角估计算法中，MUSIC 算法与ESPRIT算法属于特征结构子空间算法，是波达角估计算法中的基石。在前面的文章 一文读懂MUSIC算法DOA估计的数学原理并仿真 中详细介绍了一维MUSIC算法即线阵MUSIC算法DOA估计的原理及仿真，本文将介绍二维MUSIC算法即均匀面阵的MUSIC算法DOA估计原理及MATLAB仿真。

---

提示：以下是本篇文章正文内容，尊重版权，引用请附上链接。

一、二维MUSIC算法原理

下图为面阵入射信号模型，
​
\;\;\;\;\; 假设从远场有 $K$ 个互不相关的窄带信号，入射到一个阵元个数为 $M\times N$ 的平面阵列上。记第$i$个入射信号的方位角和俯仰角分别为${\theta }_i$和${\phi }_i$ ，则阵列接收信号可以表示为：
$$\mathbf{z}(t)=\mathbf{A}\mathbf{s}(t)+\mathbf{n}(t)$$其中$\mathbf{A}$是维度为(MN×K)的均匀矩形阵列的阵列流形，可以表示为如下所示的式子：
$$\mathbf{A}={\left[\begin{array}{c}\mathbf{a}({\theta }_k,{\phi }_1),\mathbf{a}({\theta }_2,{\phi }_2),\cdots ,\mathbf{a}({\theta }_K,{\phi }_K)\end{array}\right]}^T$$ $\mathbf{a}({\theta }_k,{\phi }_k)$为第k个入射信号的导向矢量，仅仅由阵列的阵元排布和参考阵元的选择所决定，用公式可以表示为：
$$\mathbf{a}({\theta }_k,{\phi }_k)={\mathbf{a}}_x({\theta }_k,{\phi }_k)\otimes {\mathbf{a}}_y({\theta }_k,{\phi }_k)\in C^{MN\times 1}$$ 其中$\otimes$表示的是克罗内克内积（Kronecker Product），${\mathbf{a}}_x({\theta }_k,{\phi }_k)$表示x轴方向上均匀线阵接收信号的方向矢量，${\mathbf{a}}_y({\theta }_k,{\phi }_k)$表示y轴方向上均匀线阵接收信号的方向矢量，可分别写为如下数学表达式：
$${\mathbf{a}}_x({\theta }_k,{\phi }_k)={\left[\begin{array}{c}{a}_{x,0}({\theta }_k,{\phi }_k),a_{x,1}({\theta }_k,{\phi }_k),\cdots ,a_{x,M-1}({\theta }_k,{\phi }_k)\end{array}\right]}^T$$ $${\mathbf{a}}_y({\theta }_k,{\phi }_k)={\left[\begin{array}{c}{a}_{y,0}({\theta }_k,{\phi }_k),a_{y,1}({\theta }_k,{\phi }_k),\cdots ,a_{y,N-1}({\theta }_k,{\phi }_k)\end{array}\right]}^T$$ 式中的$\mathbf{s}(t)$是信号源矢量，$\mathbf{n}(t)$为高斯白噪声矢量，服从$N(0,{\sigma }^2)$分布，可以分别表示如下式子：
$$\mathbf{s}(t)={\left[{\mathbf{s}}_0(t),{\mathbf{s}}_1(t),\cdots ,{\mathbf{s}}_{K-1}(t)\right]}^T$$ $$\mathbf{n}(t)={\left[{\mathbf{n}}_0(t),{\mathbf{n}}_1(t),\cdots ,{\mathbf{n}}_{MN}(t)\right]}^T$$ \;\;\;\;\; 阵列接收信号的协方差矩阵可以表示为：$$\mathbf{R}=\mathbb{E}[\mathbf{z}{\mathbf{z}}^H]$$ $$=\mathbf{A}\mathbb{E}[\mathbf{s}{\mathbf{s}}^H]{\mathbf{A}}^H+{\sigma }^2\mathbf{I}$$ $$=\mathbf{A}{\mathbf{R}}_S{\mathbf{A}}^H+{\sigma }^2\mathbf{I}$$ 其中${\mathbf{R}}_S$表示入射信号的协方差矩阵，${\sigma }^2\mathbf{I}$表示功率为${\sigma }^2$的高斯白噪声的协方差矩阵。
\;\;\;\;\; 实际应用中天线阵列获取的信息是有限次的快拍，因此只能得到协方差矩阵的估计值$\hat{\mathbf{R}}$，其计算公式如下：
$$\hat{\mathbf{R}}=\frac{1}{J}\sum _{j=1}^J\mathbf{z}(j){\mathbf{z}}^H(j)$$ \;\;\;\;\; 由于接收信号的协方差矩阵$\mathbf{R}$是对称矩阵，因此可以对其进行特征值分解，可以得到：
$$\mathbf{R}=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^T$$ 其中$\mathbf{U}$为$\mathbf{R}$的特征向量构成的矩阵，$\mathbf{\Lambda }$是一个由特征值构成的对角矩阵。
Λ = d i a g { λ 1 , λ 2 , . . . , λ M N } \boldsymbol{\Lambda} = diag\{ \lambda_1,\lambda_2,...,\lambda_{MN} \} Λ=diag{λ1​,λ2​,...,λMN​} \;\;\;\;\; 假设对角矩阵中的特征值降序排列，满足如下关系：
$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_K>{\lambda }_K+1=\cdots ={\lambda }_{MN}={\sigma }^2$$ 由前$K$个较大的特征值构成的对角矩阵${\mathbf{\Lambda }}_S$，其对应的特征向量构成的矩阵${\mathbf{U}}_S$为信号子空间。由后$M-K$个较小的特征值构成的对角矩阵${\mathbf{A}}_N$，其对应的特征向量构成的矩阵${\mathbf{U}}_N$为噪声子空间。

\;\;\;\;\; 根据前文假设，信号与噪声相互独立，因此信号子空间与噪声子空间是相互正交的，故信号阵列流矢量与噪声子空间也具有正交性。同一维MUSIC算法一样，可构造二维空间谱函数：
$$P_{2D-MUSIC}(\theta ,\varphi )=\frac{1}{{\mathbf{a}}^H(\theta ,\varphi ){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{a}(\theta ,\varphi )}$$ \;\;\;\;\; 当天线阵列的方向矢量与噪声子空间近似正交时，上式分母部分取极小值，空间谱函数在此时取得极大值，得到空间谱的谱峰。对空间谱进行谱峰搜索，就能够得到入射信号的方位角与俯仰角的角度，至此完成了对于信源的二维 DOA估计。

二、二维MUSIC算法MATLAB仿真

\;\;\;\;\; 参数设置如下：改变任何一个参数，仿真结果都会跟着改变，可以通过修改参数观察不同条件对估计结果的影响。

M=3;           % x轴阵元个数
N=2;           % y轴阵元个数
K=1024;        % 快拍数
fc=100e+6;     % 载波
fs=300e+6;     % 采样频率
Pn=1;          % 噪声功率fines=[45 180 250 300]; % 信号入射方位角
thetas=[5 30 55 75];    % 信号入射俯仰角
signal_f=[15e6 30e6 45e6 60e6]; % 信号频率
signal_SNR=[30 30 30 30];       % 信噪比m=(0:M-1)';    % x轴坐标
n=(0:N-1)';    % y轴坐标
c=3e+8;        % 光速
lamda=c/fc;    % 波长
dx=1/2*lamda;  % x轴阵元间距
dy=1/2*lamda;  % y轴阵元间距

\;\;\;\;\; 通过观察参数，可以得出以下结论，可以自己通过改变参数来验证，这里就不贴图了。
1、随着阵元数目的增大，MUSIC 算法的分辨率逐渐增强。
2、随着信号信噪比的增大，MUSIC 算法的分辨率逐渐增强。
3、当阵元间距与波长的比值为二分之一时，MUSIC算法能够有效进行 DOA 估计；当阵元间距小于波长的二分之一时，MUSIC 算法的分辨率会降低；当阵元间距大于波长的二分之一时，由于采样严重不足，MUSIC算法可能会丧失分辨能力。

三、MATLAB源代码

均匀面阵MUSIC算法DOA估计MATLAB仿真源代码

---

总结

\;\;\;\;\; 以上就是今天记录的所有内容，分享了均匀面阵MUSIC算法DOA估计的原理及其在MATLAB软件上仿真的结果。