[阵列信号处理]近场DOA估计算法-降维MUSIC方法（第二集）

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

背景

  近场方向到达（Direction of Arrival, DOA）估计是阵列信号处理领域的一个重要研究课题。DOA估计可以分为远场DOA估计和近场DOA估计。在近场DOA估计中，由于信号源距离阵列较近，因此信号的波前会呈现出球面波特性，这与远场情况下的平面波假设不同。近场DOA估计在很多实际应用中都非常重要，例如，无线通信、雷达、声纳、地震学以及生物医学等领域。
  前文中已经提过近场源估计的2D-MUSIC算法，其需要复杂度极高的二维联合搜索，接下来将提供另一种子空间方法，降维MUSIC,它将二维联合搜索转变为一维搜索，降低了一定的计算复杂度，在时效性上有更好的效果，并且降维MUSIC方法与2D-MUSIC方法的参数估计性能非常相近。
  接下来介绍近场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$

降维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+1\right)+\varphi {\left(-N+1\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的方法的流程，对于降维MUSIC来说，需要对其第$k$个入射的近场信号的导向矢量$a\left({\theta }_k,r_k\right)$进行进一步处理。
  根据阵列结构的对称性，导向矢量可以被分解为：
$$a\left({\theta }_k,r_k\right)=\left[\begin{array}{cccc}{\text{e}}^{\text{j}\left(-N\right){\gamma }_k}& & & \\ & {\text{e}}^{\text{j}\left(-N+1\right){\gamma }_k}& & \\ & & \ddots & \\ & & & 1\\ & & & \\ & {\text{e}}^{-\text{j}\left(-N+1\right){\gamma }_k}& & \\ {\text{e}}^{-\text{j}\left(-N\right){\gamma }_k}& & & \end{array}\right]\cdot \left[\begin{array}{c}{\text{e}}^{\text{j}{\left(-N\right)}^2{\varphi }_k}\\ {\text{e}}^{\text{j}{\left(-N+1\right)}^2{\varphi }_k}\\ ⋮\\ {\text{e}}^{\text{j}{\left(-1\right)}^2{\varphi }_k}\\ 1\end{array}\right]=a_1\left({\gamma }_k\right)a_2\left({\varphi }_k\right)$$
  其中，${\gamma }_k=-2\pi \frac{d}{\lambda }\sin {\theta }_k,{\varphi }_k=\pi \frac{d^2}{\lambda r_k}{\cos }^2{\theta }_k$,并且$a_1\left({\gamma }_k\right)$仅包含角度信息，$a_2\left({\varphi }_k\right)$同时包含角度和距离信息。
将分解后的导向矢量带入原来的2D-MUSIC谱峰搜索器：
$$\begin{array}{rl}P\left(\theta ,r\right)& =\frac{1}{a_2^{\text{H}}\left(\varphi \right)a_1^{\text{H}}\left(\gamma \right){\mathbf{U}}_n{\mathbf{U}}_n^{\text{H}}{a}_1\left(\gamma \right)a_2\left(\varphi \right)}\\ & =\frac{1}{a_2^{\text{H}}\left(\varphi \right)Q\left(\gamma \right)a_2\left(\varphi \right)}\end{array}$$
其中，$Q\left(\gamma \right)=a_1^{\text{H}}\left(\gamma \right){\mathbf{U}}_n{\mathbf{U}}_n^{\text{H}}{a}_1\left(\gamma \right)$,考虑到一个恒等式$e_1^{\text{H}}{a}_2\left(\varphi \right)=1$,而$e_1={\left[\begin{array}{c}0,\cdots ,1\end{array}\right]}^{\mathrm{T}}\in {\mathbb{R}}^{M\times 1}$,可以将对$P\left(\theta ,r\right)$的角度距离二维搜索看成一个二次优化问题，该二次规划问题的约束条件是$e_1^{\text{H}}{a}_2\left(\varphi \right)=1$，目标则是最大化功率$P\left(\theta ,r\right)$，即最小化$a_2^{\text{H}}\left(\varphi \right)a_1^{\text{H}}\left(\gamma \right){\mathbf{U}}_n{\mathbf{U}}_n^{\text{H}}{a}_1\left(\gamma \right)a_2\left(\varphi \right)=a_2^{\text{H}}\left(\varphi \right)Q\left(\gamma \right)a_2\left(\varphi \right)$
  最终，该二次优化问题可以写成如下标准形式
$$\begin{gathered} \underset{\gamma ,\varphi }{\min }{\mathbf{a}}_2^{\mathrm{H}}(\varphi )\mathbf{Q}(\gamma ){\mathbf{a}}_2(\varphi ) \\ s.t.{\mathbf{e}}_1^{\mathrm{H}}{\mathbf{a}}_2(\varphi )=1 \end{gathered}$$
  根据拉格朗日乘数法，可以构造如下代价函数。
$$L(\gamma ,\varphi )=a_2^{\mathrm{H}}(\varphi )Q(\gamma )a_2(\varphi )-\omega [e_1^{\mathrm{H}}{a}_2(\varphi )-1]$$
  对$a_2^{\mathrm{H}}(\varphi )$求偏导，可得：$\frac{\partial L(\gamma ,\varphi )}{\partial {\mathbf{a}}_2(\varphi )}=2\mathbf{Q}(\gamma ){\mathbf{a}}_2(\varphi )+\omega {\mathbf{e}}_1=0$,
  对上式进行化简可得，
$$a_2\left(\varphi \right)=-\frac{\omega }{2}{Q}^{-1}\left(\gamma \right)e_1=\mu Q^{-1}\left(\gamma \right){\mathbf{e}}_1$$其中，$\mu =-\frac{\omega }{2}$。
  而由于${\mathbf{e}}_1^{\mathrm{H}}{\mathbf{a}}_2(\varphi )=1$,带入上式，可得
$$a_2(\varphi )=\frac{Q(\gamma {)}^{-1}{e}_1}{e_1^{\mathrm{H}}Q(\gamma {)}^{-1}{e}_1}$$
  因此可得到${\gamma }_k$的估计值为：
$${\hat{\gamma }}_k=\arg \underset{\gamma }{\min }\frac{1}{e_1^{\mathrm{H}}Q({\gamma }_k{)}^{-1}{e}_1}=\arg \underset{\gamma }{\max }{e}_1^{\mathrm{H}}Q({\gamma }_k{)}^{-1}{e}_1$$
  而在信号模型中，我们知道${\gamma }_k=-2\pi d\sin {\theta }_k/{\lambda }_k,k=1,2,...,K$，将信号波长归一化后可以知道，$-\pi /2⩽{\gamma }_k⩽\pi /2$,因此在上述区间内进行一维谱峰搜索即可知道对应的估计值${\hat{\gamma }}_k$
  将$K$个估计值${\hat{\gamma }}_k$带入$a_2(\varphi )$中，即可得到：
$${\mathbf{a}}_2({\varphi }_k)=\left[\begin{array}{c}{\mathrm{e}}^{\mathrm{j}(-N{)}^2{\varphi }_k}\\ \\ {\mathrm{e}}^{\mathrm{j}(-N+1{)}^2{\varphi }_k}\\ ⋮\\ \\ {\mathrm{e}}^{\mathrm{j}(-1{)}^2{\varphi }_k}\\ \\ 1\end{array}\right]$$
  对上述矩阵顺时针旋转180度可以得到下式：
$$\begin{gathered} {\mathbf{a}}_2^{\prime }({\varphi }_k)=\left[\begin{array}{c}{\mathrm{e}}^{\mathrm{j}(-N{)}^2{\varphi }_k}\\ \\ {\mathrm{e}}^{\mathrm{j}(-N+1{)}^2{\varphi }_k}\\ ⋮\\ \\ {\mathrm{e}}^{\mathrm{j}(-1{)}^2{\varphi }_k}\\ \\ 1\end{array}\right] \\ {\hat{g}}_k=\mathrm{angle}[a_2^{\prime }({\varphi }_k)] \end{gathered}$$
其中，angle表示取相角。
${\hat{g}}_k$可以表示为：
$${\hat{g}}_k=\left[\begin{array}{c}0\\ {\varphi }_k\\ 4{\varphi }_k\\ ⋮\\ (N-1{)}^2{\varphi }_k\\ \\ N^2{\varphi }_k\end{array}\right]=\left[\begin{array}{c}0\\ 1\\ 4\\ ⋮\\ (N-1{)}^2\\ \\ N^2\end{array}\right]{\varphi }_k=\mathbf{q}{\varphi }_k$$
其中，$\mathbf{q}={\left[\begin{array}{ccc}0,& 1,& 4,\cdots ,(N-1{)}^2,N^2\end{array}\right]}^{\mathrm{T}}$,可以由最小二乘法求出${\hat{\varphi }}_k$:
$$\underset{c_k}{\min }\parallel p{c}_k-{\hat{g}}_k{\parallel }_{\mathrm{F}}^2$$其中，${\hat{g}}_k=\text{angle}[{\mathbf{a}}_2^{\prime }({\varphi }_k)],\mathbf{p}=[1_{N+1},\mathbf{q}],{\mathbf{c}}_k=[c_{k0},{\hat{\varphi }}_k{]}^{\mathrm{T}}\in {\mathbb{R}}^{2\times 1},c_{k0}是最小二乘法的损失，{\hat{\varphi }}_k为待估计的参数，{\mathbf{c}}_k的解为：$
$${\mathbf{c}}_{\mathbf{k}}=(p^{\mathrm{T}}p{)}^{-1}{p}^{\mathrm{T}}{\hat{g}}_k$$
  上述过程中，角度和距离参数是自动配对的，最终可以得到角度的估计值${\hat{\theta }}_k$和距离的估计值${\hat{r}}_k$
$$\begin{gathered} {\hat{\theta }}_k=-\arcsin \left(\frac{{\hat{\gamma }}_k{\lambda }_k}{2\pi d}\right) \\ {\hat{r}}_k=\frac{\pi d^2}{{\lambda }_k{\hat{\varphi }}_k}{\cos }^2{\hat{\theta }}_k \end{gathered}$$

  与上文2d-MUSIC算法相同，使用同样参数的信号，利用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
%%%%%%%%%降维MUSIC算法%%%%%%%%%

derad = pi/180; %角度->弧度
M = 2; % 信源数目

lamta=1; %波长1m

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_k1=-2*pi*d.*sin(xita_k)./lamta;
gama_k1=gama_k1';
phi_k=pi*d*d.*cos(xita_k).*cos(xita_k)./(r_k.*lamta);
phi_k=phi_k';

A_x=exp(1i*(gama_k1.*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;

n=-Nx:1:0;
for iang=1:360
    angleall(iang)=(iang-180)/2;
    gama_k=-2*pi*d.*sind(angleall(iang))./lamta;
    gama_k=gama_k';

    A1_x=zeros(2*Nx+1,Nx+1);
    for k=1:1:Nx+1
        A1_x(k,k)=exp(1i*(k-Nx-1).*gama_k);
    end
    for k=1:1:Nx
        A1_x(2*Nx+2-k,k)=exp(1i*(Nx+1-k).*gama_k);
    end
    A2_x=exp(1i*phi_k.*n.*n);
    A2_x=A2_x';
    Q=A1_x'*Un*Un'*A1_x;
    V=A2_x'*Q*A2_x;

    e1=zeros(length(n)-1,1);
    e1=[e1;1];
    gama_hateqe1=e1'*inv(Q)*e1;

    gama_hateqe(iang)=gama_hateqe1;
    gama_hateqe=abs(gama_hateqe);
    %P_search(num_distance,iang)=1/abs(A_x'*Un*Un'*A_x);

end


[pks locs]=findpeaks(gama_hateqe);
[pks_max2 locs_max2i]=maxk(pks,2);
loc_max2=locs(locs_max2i);
gama_hat=-2*pi*d.*sind(angleall(loc_max2))./lamta;



A1_x_pie1=zeros(2*Nx+1,Nx+1);
for k=1:1:Nx+1
    A1_x_pie1(k,k)=exp(1i*(k-Nx-1)*gama_hat(1));
end
for k=1:1:Nx
    A1_x_pie1(2*Nx+2-k,k)=exp(1i*(Nx+1-k)*gama_hat(1));
end
Q_pie1=A1_x_pie1'*Un*Un'*A1_x_pie1;
A2_x_pie1=inv(Q_pie1)*e1/(e1'*inv(Q_pie1)*e1);
A2_x_pie1=A2_x_pie1(end:-1:1);%翻转矩阵180度
g_hat1=angle(A2_x_pie1);
g_hat1=fliplr(n)'.*fliplr(n)'*g_hat1(2);


A1_x_pie2=zeros(2*Nx+1,Nx+1);
for k=1:1:Nx+1
    A1_x_pie2(k,k)=exp(1i*(k-Nx-1)*gama_hat(2));
end
for k=1:1:Nx
    A1_x_pie2(2*Nx+2-k,k)=exp(1i*(Nx+1-k)*gama_hat(2));
end
Q_pie2=A1_x_pie2'*Un*Un'*A1_x_pie2;
A2_x_pie2=(inv(Q_pie2)*e1)/(e1'*inv(Q_pie2)*e1);
A2_x_pie2=A2_x_pie2(end:-1:1);
g_hat2=angle(A2_x_pie2);
g_hat2=fliplr(n)'.*fliplr(n)'*g_hat2(2);

p=ones(Nx+1,1);
q=[];
for i=1:1:Nx+1
    q(i)=(i-1)^2;
end
q=q';
p=[p,q];


ck1=inv(p'*p)*p'*g_hat1;
ck2=inv(p'*p)*p'*g_hat2;
phi1=-ck1(2);
phi2=-ck2(2);
%a=gama_hat(1)*lamta/(2*pi*d);
%角度估计值
theta1=-(-asind(gama_hat(1)*lamta/(2*pi*d)))
theta2=-(-asind(gama_hat(2)*lamta/(2*pi*d)))
%距离估计值
r1=pi*d^2/(lamta*phi1)*(cosd(theta1)^2)
r2=pi*d^2/(lamta*phi2)*(cosd(theta2)^2)

运行结果中，角度估计值如下图所示：

距离估计值如下图所示：

距离估计值还是有一定的误差，该方法采用一维搜索的方式，精度不如二维搜索，但是大大提高了算法的实时性，损失了一定的估计精度。