[阵列信号处理]近场DOA估计算法-2DMUSIC方法（第一集）

  本文的主要参考资料是张小飞老师的《阵列信号处理及MATLAB实现（第3版）》（电子工业出版社出版）。

背景

  近场方向到达（Direction of Arrival, DOA）估计是阵列信号处理领域的一个重要研究课题。DOA估计可以分为远场DOA估计和近场DOA估计。在近场DOA估计中，由于信号源距离阵列较近，因此信号的波前会呈现出球面波特性，这与远场情况下的平面波假设不同。近场DOA估计在很多实际应用中都非常重要，例如，无线通信、雷达、声纳、地震学以及生物医学等领域。
  近场信号和远场信号有以下特点：

1. 近场与远场的区别：
   • 近场区域：信号源到阵列的距离远小于阵列的孔径大小，波前是球面波。
   • 远场区域：信号源到阵列的距离远大于阵列的孔径大小，波前可以近似为平面波。
2. 近场DOA估计的挑战：
   • 近场效应导致信号波前的复杂性增加。
   • 需要精确估计信号源的距离和角度信息。
   • 阵列校准和误差的影响更加敏感。
3. 近场DOA估计算法：
   • 基于子空间的方法：如Music算法、ESPRIT算法等，这些算法需要修改以适应近场条件。
   • 基于最大似然的方法：考虑信号源距离的影响，通过最大化似然函数来估计DOA。
   • 基于稀疏表示的方法：利用信号的稀疏性质，通过求解稀疏重构问题来估计DOA。
   • 基于深度学习的方法：近年来，深度学习技术也被应用于DOA估计，通过训练神经网络来直接从数据中估计DOA。

近场信号模型

  接下来介绍近场DOA估计中的二维MUSIC算法（2D-MUSIC）,在信号模型中，采用均匀线阵ULA(Uniform Linear Array)，所谓二维就是指2个待估计的参数，距离$r$和到达角$\theta$。信号模型图如下所示：

  阵列天线由$M=2N+1$个传感器阵元组成，相邻阵元间的间距为$d$,以中心标号为0的传感器阵元为中心阵元，左右两边各分布$N$个单边阵元，左右两边阵元共享中心标号为0的传感器阵元。因此，阵元$m\in [-N,\cdots ,-2,-1,0,1,2,\cdots ,N]$,并且假设有$K$个窄带近场非相干信号入射，在上图信号模型中画出了第$k（1,2,...,K）$个近场信号入射，距离中心参考阵元距离$r_k$,到达角为${\theta }_k$。那么其接收信号矩阵$\mathbf{x}\left(t\right)$可以写成如下形式
$$\begin{gathered} \mathbf{x}\left(t\right)={\left[x_{-N}\left(t\right),\cdots ,x_{-1}\left(t\right),x_0\left(t\right),x_1\left(t\right),\cdots ,x_N\left(t\right)\right]}^T \\ =\left[\mathbf{a}\left({\theta }_1,r_1\right),\mathbf{a}\left({\theta }_2,r_2\right),\cdots ,\mathbf{a}\left({\theta }_K,r_K\right)\right]\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right) \\ ={\mathbf{A}}_x\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right) \end{gathered}$$
其中，${\mathbf{A}}_x$是阵列流形,$\mathbf{s}\left(t\right)$是信号矩阵，$\mathbf{n}\left(t\right)$是加性高斯白噪声矩阵。
根据近场信号的菲涅尔近似，第$k$个信号的导向矢量$\mathbf{a}\left({\theta }_k,r_k\right)$可以被写成$a\left({\theta }_k,r_k\right)={\left[{\text{e}}^{\text{j}\left({\gamma }_k\left(-N+1\right)+{\varphi }_k{\left(-N+1\right)}^2\right)},...,{\text{e}}^{\text{j}\left({\gamma }_k\left(-1\right)+{\varphi }_k{\left(-1\right)}^2\right)},1,{\text{e}}^{\text{j}\left({\gamma }_{k}1+{\varphi }_k{1}^2\right)},...,{\text{e}}^{\text{j}\left({\gamma }_{k}N+{\varphi }_k{N}^2\right)}\right]}^{\text{T}}$并且，${\gamma }_k=-2\pi \frac{d}{\lambda }\sin {\theta }_k,{\varphi }_k=\pi \frac{d^2}{\lambda r_k}{\cos }^2{\theta }_k$

2D-MUSIC

  假设采样时间大于相干时间，有限次采样得到的接收信号的协方差矩阵可表示为如下所示，并对该协方差矩阵进行特征分解：
$${\hat{R}}_x=x(t)x^{\mathrm{H}}(t)/K={\mathbf{U}}_{\text{S}}{\mathbf{\Sigma }}_{\text{S}}{\mathbf{U}}_{\text{S}}^{\text{H}}+{\mathbf{U}}_{\text{N}}{\mathbf{\Sigma }}_{\text{N}}{\mathbf{U}}_{\text{N}}^{\text{H}}$$
其中，$\mathbf{x}\left(t\right)$是接收信号矩阵，$K$是快拍数。$(\bullet {)}^{\mathrm{H}}$是共轭转置，${\mathbf{U}}_S$是由大特征值对应的特征向量组成的信号子空间，${\mathbf{U}}_N$是由小特征值对应的特征向量组成的噪声子空间。
  根据子空间理论，在近场时，噪声子空间和信号子空间的正交性仍然成立，因此，基于远场的MUSIC方法在近场时仍然适用。
  在搜索区间 $r\in [0.62(D^3/\lambda {)}^{1/2},(2D^2/\lambda )],\theta \in [-\pi /2,\pi /2]$内构造谱峰搜器，$D=Md$。
$$P\left(\theta ,r\right)=\frac{1}{a^H\left(\theta ,r\right)U_n{U}_n^{H}a\left(\theta ,r\right)}$$
其中，$a\left(\theta ,r\right)={\left[{\text{e}}^{\text{j}\left(\gamma \left(-N\right)+\varphi {\left(-N\right)}^2\right)},...,{\text{e}}^{\text{j}\left(\gamma \left(-1\right)+\varphi {\left(-1\right)}^2\right)},1,{\text{e}}^{\text{j}\left(\gamma 1+\varphi 1^2\right)},...,{\text{e}}^{\text{j}\left(\gamma N+\varphi N^2\right)}\right]}^{\text{T}}$并且，$\gamma =-2\pi \frac{d}{\lambda }\sin \theta ,\varphi =\pi \frac{d^2}{\lambda r}{\cos }^2\theta$。
得到对应峰值的下标，即目标对应的距离与角度:
$$\begin{array}{rl}& \hat{\theta }=[{\hat{\theta }}_1,{\hat{\theta }}_2,\cdots ,{\hat{\theta }}_K]\\ & \hat{r}=[{\hat{r}}_1,{\hat{r}}_2,\cdots ,{\hat{r}}_K]\end{array}$$
其中$\hat{\theta }$,$\hat{r}$分别是对应信号相对于参考阵元的距离和到达角度。
  2D-MUSIC的算法流程可以表示如下：
  Step1:根据采样得到的信号矩阵构造接收信号的协方差矩阵的估计值。
  Step2:对接收信号的协方差矩阵进行特征值分解，然后按照特征值的大小顺序，把与信号个数相等的较大的特征值和对应的特征向量看成信号子空间，把剩下的对应部分看成噪声子空间$U_n$。
  Step3:根据噪声特征向量和信号向量的正交性，构造空间谱峰搜索函数，对其进行谱峰搜索可以得到$K$组配对的信源目标角度与距离参数的估计结果。
接下来利用2个近场非相干信号源进行仿真，其距离参考阵元的距离$r_k$和到达角${\theta }_k$分别为$\left(20°,1.85\lambda \right)$、$\left(40°,4.8\lambda \right)$,阵元间距$d=\lambda /4$,信噪比$SNR=20dB$,单边阵元数$N=4$,快拍数$K=2000$。
matlab代码如下：

clc
clear all
close all
%%%%%%近场2D-MUSIC算法
derad = pi/180; %角度->弧度
M = 2; % 信源数目

lamta=1;% 信号波长
Nx =4; % 单边阵元个数
Mx=2*Nx+1; %总阵元个数
d =lamta/4; % 阵元间距

D = (Mx-1)*d; %阵列孔径
distance_left = 0.62*sqrt((D.^3)/lamta);%距离搜索的下界
distance_right = 2*(D.^2)/lamta;%距离搜索的上界

xita_k = [20 40].*derad; % x轴待估计角度

snr =20; % 信噪比
K =2000; % 快拍数
r_k=[1.85*lamta,4.8*lamta]; %信号源距离
m=-Nx:1:Nx;%x轴入射点位

gama_k=-2*pi*d.*sin(xita_k)./lamta;
gama_k=gama_k';
phi_k=pi*d*d.*cos(xita_k).*cos(xita_k)./(r_k.*lamta);
phi_k=phi_k';
A_x=exp(1i*(gama_k.*m+phi_k.*m.*m));%构造阵列流形矩阵

% %%%%%%构建信号模型%%%%%
S=randn(M,K);
X1=A_x'*S;
X1=awgn(X1,snr,'measured'); %将白色高斯噪声添加到信号中

Rxx= X1*X1'/K;%获取接收信号的协方差矩阵
% 特征值分解
[EV,D]=eig(Rxx); %特征值分解
EVA=diag(D)'; %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA); %将特征值排序 从小到大
EV=fliplr(EV(:,I)); % 对应特征矢量排序
Un=EV(:,M+1:end);%获取x轴接收数据的噪声子空间
%进行角度和距离的二维搜索
num_distance=0;
for distance=0.01:0.01:distance_right
num_distance=num_distance+1;
for iang=1:360
angle(iang)=(iang-180)/2;
gama_k=-2*pi*d.*sind(angle(iang))./lamta;
gama_k=gama_k';
phi_k=pi*d*d.*cosd(angle(iang)).*cosd(angle(iang))./(distance.*lamta);
phi_k=phi_k';

A_x=exp(1i*(gama_k.*m+phi_k.*m.*m));
A_x=A_x';
P_search(num_distance,iang)=1/abs(A_x'*Un*Un'*A_x);

end
end

distance=0.01:0.01:distance_right;
angle(iang)=(iang-180)/2;

[distance,angle]=meshgrid(distance,angle);
%绘制三维图
mesh(distance',angle',P_search);
yticks([-90:10:90]);
ylabel('角度/°');
xticks([0:0.5:distance_right]);
xlabel('r/\lambda');
figure
contour(distance',angle',P_search)
yticks([-90:10:90]);
ylabel('角度/°');
xticks([0:0.5:distance_right]);
xlabel('r/\lambda');

  画出相应的三维谱峰搜索图和二维等高线图，可以看出算法可以有效估计信源的角度和距离参数，且参数自动配对。