DOA估计基于互质阵列的DOA估计

量大龙神魔

已于 2023-11-25 16:43:03 修改

阅读量5.1k

点赞数 21

文章标签：算法 matlab 矩阵

于 2023-04-11 11:31:15 首次发布

本文链接：https://blog.csdn.net/qq_45471796/article/details/130048892

版权

前言

传统阵列的配置方式是均匀线阵，该阵列要求相邻阵元的间距为半波长，易产生耦合效应，影响 DOA 估计精度。而稀疏阵列利用协方差矩阵构建差分共阵方式在虚拟域上生成虚拟阵列，并利用虚拟阵列实现波达方向角的估计。由于虚拟阵列的自由度不在局限于阵列的物理阵元个数，而是和阵列的配置方式有关。其中应用的较为多的稀疏阵列形式有最小冗余阵列、互质阵列、嵌套阵列三种。最小冗余阵结构是1968年由Moffet提出的一种冗余度最小的稀疏阵列结构，该阵列结构将阵列孔径归一化为L，保证阵列中所有阵元位置差必须是从1到L，同时阵元数量M需要满足最小条件。嵌套阵阵列结构可以采用多级嵌套形式，研究范围较广。互质阵由Palghat P. Vaidyanathan与Piya Pal在2011年于文献Sparse Sensing With Co-Prime Samplers and Arrays中提出，典型的互质阵主要由两个阵元间距分别为半波长整数倍的均匀线阵组成，且两个均匀子阵的阵元数M,N是互质关系。后来有学者发现，将互质阵中阵元较少、间距较大的子阵进行二倍阵列扩展可获得更多的自由度，该种阵列结构即为扩展的互质阵。

互质阵的阵列结构与空间谱估计

首先以典型的互质阵为例，如上图所示，两个子阵阵元数分别为M和N，阵元间距分别为Nd和Md， $M < N$ 。将这两个均匀子阵进行匹配即可得到一个标准的互质阵。在图中，参考阵元位于最左端开始的位置，且剩下的阵元彼此间阵元位置不会发生重合，阵元总数为 $M+N-1$ 。在构造的互质阵中，令 $M+N-1$ 个阵元的位置集合P为：

$P=\left [ P_{1}, P_{2}, ..., P_{M+N-1} \right ]$

其中 $i=1,2,...,M+N-1$ ，则第i个阵元位置可被表示为：

$P_{i} = \left \{ Mnd,0\leq n\leq N-1 \right \}\cup \left \{ Nmd,0\leq m\leq M-1 \right \}$

通过对该阵元集中的任意两个阵元位置作差集，即可得到一个新的阵列差集，该互质阵列中阵元位置差集S可表示为：

$S=\left \{ \pm \left ( Mn-Nm \right ) d,0\leq n\leq N-1, 0\leq m\leq M-1 \right \}$

以M=3，N=4为例，该互质阵的阵列差集有如下分布：

接下来便是通过真实阵元位置得到其阵列差集。

以第一个阵元为参考阵元，则有 $P_1 = 0$ 。假设Q个信号源被阵列接收，就可以得到某个信号源的导向矢量 $a(\theta_q)$ ，有：

$a(\theta_q) = \left [ 1,e^{-j2\pi P_2 sin \theta_q/\lambda},...,e^{-j2\pi P_{M+N-1} sin \theta_q/\lambda} \right ]^T$

阵列流形矩阵A可表示为：

$A = \left [ a(\theta_1) ,a(\theta_2),...,a(\theta_{Q})\right ]$

故阵列的接收信号可表示为：

$x(t) = As(t)+n(t) = \sum_{q = 1}^{Q} a(\theta_q)s_q(t)+n(t)$

其中 $s(t) = \left [ s_1(t),s_2(t),...,s_Q(t) \right ]^T$ 为信源矢量， $n(t)$ 为噪声矢量，一般情况下 $n(t) \sim (0,\sigma_{n}^{2})$ ，则 $x(t)$ 的理想协方差阵 $R_{xx}$ 可表示为：

$R_{xx} = E\{x(t)x^H(t)\} = \sum_{q=1}^{Q}\sigma_q^2a(\theta_q)a^H(\theta_q)+\sigma_n^2 I_{M+N-1}$

其有限快拍数下的近似形式为：

$R_{xx} = \frac{1}{T} \sum_{t=1}^{T} x(t)x^H(t)$

由 $R_{xx}$ 的理想形式可知， $R_{xx}$ 的第i行第j列的元素可以表示为：

$R_{xx}(i,j) = \left\{\begin{matrix} \sum_{q=1}^{Q} \sigma_q^2e^{-j2\pi(P_i-P_j)sin(\theta_q)/\lambda},i\neq j \\ \sum_{q=1}^{Q} \sigma_q^2e^{-j2\pi(P_i-P_j)sin(\theta_q)/\lambda}+\sigma_n^2,i=j \end{matrix}\right. 0\leq~ i,j\leq~ M+~N-~1$

可以发现协方差矩阵中的元素表实现形式与导向矢量中元素的表现形式相似，不同的是阵列的位置变成了任意两个阵元位置的差，信号部分变成信号源的功率。因而，就可以通过阵列中任意两个物理阵元的位置之差构造虚拟阵元，不同的虚拟阵元组合在一起就构成了虚拟阵列。所有的延迟即构成了差分共阵。想要获得阵列中虚拟阵元的信息，需要对 $R_{xx}$ 向量化，定义 $z$ 为 $R_{xx}$ 向量化后的列向量，则 $z$ 可表示为：

$z = vec(R_{xx}) = Mp+\sigma_n^^2vec(I_{M+N-1})$

其中 $M = \left [ A^{*}\circ A\right ]$ 为扩大孔径的虚拟阵列流形矩阵； $\circ$ 为 Khatri-Rao 积；可以发现 $M$ 包含了该互质阵中所有的两两阵元之差的信息。

空间平滑需要取出z中连续的虚拟阵元，然而经过向量化处理后得到的列向量 z 含有许多冗余项，并且每一项对应的虚拟阵元位置是杂乱的，需要进一步处理。

因此，引入如下两个定义：

定义1：虚拟阵元权值函数：令虚拟阵列中第p个虚拟阵元重复的次数为该虚拟阵元的权值函数 $w(p)$ ，则 $w(p)$ 可表示为：

$w(p) = \{ P_i,P_j \in P| P_i-P_j = P \}$

定义2：冗余阵元整合矩阵：令 $J\in C^{(M+N-1)^2\times (2MN+1) }$ 为冗余阵元整合矩阵，则 $J$ 中对应的第p个虚拟阵元的行向量可以表示为：

$J(p,:) = vec(I(p))$

$I(p)(P_i,P_j)=\left\{\begin{matrix} \frac{1}{w(p)},P_i-P_j \\ 0,else \end{matrix}\right.$

则去冗余之后的信号的表达式为：

$v=Jz$

包含了任意的虚拟阵元，但由于该向量 $v$ 的秩为1，而且此时该信号中的虚拟孔径非连续，为一非均匀稀疏阵列。因此需要首先取连续的虚拟孔径部分作为虚拟均匀线性阵列。

在获得虚拟线性阵列后，再进一步对获得的向量进行平滑处理或者托普利茨化成为一个虚拟协方差矩阵。

空间平滑处理：假设z中连续的虚拟孔径为 $U=[-L_U,L_U]$ ，截取 v 中连续孔径部分 $v_U$ ，则 $v_U$ 中含有 $2L_U+1$ 个元素。在进行空间平滑操作时，每次选取 $L_U+1$ 个连续阵元，一共需要选取 $L_U+1$ 次，记第 $i$ 次选取的连续阵元用 $v_{Ui}$ 表示，则经过空间平滑后的协方差矩阵 $R_{UU}$ 可表示为：

$R_{UU} = \frac {1}{L_U+1}\sum_{i=1}^{L_U+1}v_{Ui}v_{Ui}^H$

$R_{UU}$ 就是经过空间平滑获得的满秩协方差矩阵，进而可利用该矩阵 $R_{UU}$ 进行 DOA 估计。

托普利茨构造：简单来说，均匀线阵的协方差阵在理想情况下由于其范德蒙德结构的特性可以等价一个复托普利茨矩阵，即斜对角线上的各个元素都相等，仍然假设z中连续的虚拟孔径为 $U=[-L_U,L_U]$ ，则根据该虚拟孔径中的 $[-L_U,0]$ 部分及 $[0,L_U]$ 部分即可通过对这两个部分的元素进行托普利茨化得到该协方差阵。

代码与仿真

如果诸位大佬发现仿真过程中有错误，请不吝指出谢谢。具体代码如下所示：

dis_sub_array_1 = 3;%%阵元数量多，阵元间距短的子阵的间距
dis_sub_array_2 = 5;%%阵元数量少，阵元间距大的子阵的间距
multi = 2;%%阵列数量少的子阵的扩展倍数
sub_array_1_num = dis_sub_array_2;
sub_array_2_num = dis_sub_array_1;
sub_array_2_num = sub_array_2_num*multi;
array_num = sub_array_1_num+sub_array_2_num-1;
c = 340;%%波速
f = 1000;%%频率
lambda = c/f;%%波长
d = 0.5*lambda;
snr = 10;
array_structure = [0,dis_sub_array_1:dis_sub_array_1:(sub_array_1_num-1)*dis_sub_array_1,dis_sub_array_2:dis_sub_array_2:(sub_array_2_num-1)*dis_sub_array_2];
array_structure = sort(array_structure,"ascend");%%阵列真实结构构造
source_num = 4;
theta = linspace(-60,60,source_num);%%角度
snapshot_num = 100;%%快拍数
A = exp(-1i*array_structure'*d*f/c*2*pi*sind(theta));
S=randn(source_num,snapshot_num);
X = A*S;
X = awgn(X,snr,'measured');
RR = X*X'/snapshot_num;
z = reshape(RR,[],1);%%协方差矩阵向量化
W = zeros(array_num,array_num);
for ik = 1:array_num
    for jk = 1:array_num
        W(ik,jk) = array_structure(ik)-array_structure(jk);
    end
end
W = reshape(W,1,[]);%%虚拟阵列阵元分布
max_flag = round(max(W));
min_flag = round(min(W));
histog = zeros(1,max_flag-min_flag+1);
for ik = 1:length(W)
    histog(round(W(ik))+abs(min_flag)+1) = histog(round(W(ik))+abs(min_flag)+1)+1;
end%%虚拟权值函数构建
J = [];
for ik = min(W):max(W)
    II = zeros(array_num,array_num);
    for jk = 1:array_num
        for kk = 1:array_num
            if array_structure(jk) - array_structure(kk) == ik
                II(jk,kk) = 1/(histog(ik+abs(min(W))+1));
            end
        end
    end
    J = [J;reshape(II,1,[])];%%冗余阵元整合矩阵构建
end
v = J*z;%%去冗余之后的信号的表达式
histog_flag = round(abs(min_flag)+1);
for ik = 0:min(abs(min(W)),abs(max(W)))
    if histog(histog_flag+ik) == 0 || histog(histog_flag-ik) == 0
        flag = ik-1;
        break;
    end
end
vir_array = v(histog_flag-flag:histog_flag+flag);%%连续虚拟阵列孔径划分
R = zeros(flag+1,flag+1);
for ik = 1:flag+1
    R = R+vir_array(ik:ik+flag)*vir_array(ik:ik+flag)';
end%%空间平滑处理
% R = toeplitz(v(histog_flag:histog_flag+flag,1),flipud(v(histog_flag-flag:histog_flag,1)));%%向量托普利茨化
[U,D,V] = svd(R);
Un = U(:,source_num+1:end);
deg = -90:0.01:90;
kkk = 0;
for ik = 1:length(deg)
    kkk = kkk+1;
    A_bar = exp(-1i*(0:size(R,1)-1)'*d*f/c*2*pi*sind(deg(ik)));
    music_p(kkk) = 1/(A_bar'*Un*Un'*A_bar);
end
figure();
semilogy(deg,abs(music_p));grid on;
xlabel('degree');ylabel('amplitude');
title('CPA-MUSIC谱估计结果');
figure();
stem(min_flag:max_flag,histog);grid on;
xlabel("位置");
ylabel("重复次数");
title("虚拟阵元分布及重复次数");

仿真1：单信源条件下的互质阵空间谱估计

部分仿真参数如下所示：

子阵1阵元数	5
子阵2阵元数	3
子阵2扩展倍数	1
阵元间距	半波长
角度	15°
快拍数	100
信噪比	10dB

仿真结果如下所示：

仿真2：多信源下的互质阵空间谱估计

同样以阵元数为6的阵列为例，与同等阵元数目下的均匀阵DOA估计角度数量作对比可发现，均匀阵只能最多估计5个来波方向，而互质阵可以估计6个方向，使得阵列自由度增加。

部分仿真参数如下所示：

子阵1阵元数	5
子阵2阵元数	3
子阵2扩展倍数	1
阵元间距	半波长
角度	-60°，-36°，-12°，12°，36°，60°，
快拍数	100
信噪比	10dB

仿真结果如下所示：

仿真3：多信源下的扩展互质阵空间谱估计

以前两个仿真中的互质阵结构为基础，对该结构中阵元间距较大的子阵进行二倍扩展，即可获得更多的扩展自由度。

部分仿真参数如下所示：

子阵1阵元数	4
子阵2阵元数	3(实际6)
子阵2扩展倍数	2
阵元总数	9
阵元间距	半波长
角度	-60°至60°中均匀取14个角度
快拍数	100
信噪比	10dB