基于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),?,\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)?{\mathbf{a}}_y({\theta }_k,{\phi }_k)\in C^{MN\times 1}$$ 其中$?$表示的是克羅內克內積（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),?,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),?,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),?,{\mathbf{s}}_{K?1}(t)\right]}^T$$ $$\mathbf{n}(t)={\left[{\mathbf{n}}_0(t),{\mathbf{n}}_1(t),?,{\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 ?\ge {\lambda }_K>{\lambda }_K+1=?={\lambda }_{MN}={\sigma }^2$$ 由前$K$個較大的特征值構成的對角矩陣${\mathbf{\Lambda }}_S$，其對應的特征向量構成的矩陣${\mathbf{U}}_S$為信號子空間。由后$M?K$個較小的特征值構成的對角矩陣${\mathbf{A}}_N$，其對應的特征向量構成的矩陣${\mathbf{U}}_N$為噪聲子空間。

\;\;\;\;\; 根據前文假設，信號與噪聲相互獨立，因此信號子空間與噪聲子空間是相互正交的，故信號陣列流矢量與噪聲子空間也具有正交性。同一維MUSIC算法一樣，可構造二維空間譜函數：
$$P_{2D?MUSIC}(\theta ,?)=\frac{1}{{\mathbf{a}}^H(\theta ,?){\mathbf{U}}_N{\mathbf{U}}_N^H\mathbf{a}(\theta ,?)}$$ \;\;\;\;\; 當天線陣列的方向矢量與噪聲子空間近似正交時，上式分母部分取極小值，空間譜函數在此時取得極大值，得到空間譜的譜峰。對空間譜進行譜峰搜索，就能夠得到入射信號的方位角與俯仰角的角度，至此完成了對于信源的二維 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軟件上仿真的結果。