天线阵列1——均匀线阵和互质阵
1.基本假设 在论文当中，所有的研究主要基于以下几个基本假设：假设入射源位于阵列无限远处，从而阵列接收到的信号可以看作是平面波；假设信号带宽远小于信号的中心频率，即窄带信号，并且满足各态历经性；所有阵元为全向接收阵元，且阵元之间不存在互耦效应。阵元已完全校正，不存在位置误差和幅相不一致的情况，且阵元之间的间隔以半波长为单位；不同时刻和不同阵元上所接收到的噪声相互独立。2.阵列信号模型 按照阵元摆放的不同形式，天线阵列主要可分为面阵和线阵两大类。其中面阵包括均匀矩形阵（Uniform Rectangular Array: URA）、稀疏矩形阵（Sparse Rectangular Array: SRA）、L/T型阵和圆阵（Circular Array: CA）等，线阵包括 ULA 和 SLA。本文将主要介绍均匀线阵和互质阵。2.1 均匀线阵模型 以一个由 M 个阵元构成的一维线性阵列为例，假设有 K 个远场窄带非相关入射信号源，来波方向分别为 θk，k = 1, 2, · · · , K，如图2.1所示。
图2.1 一维线性阵列模型 将阵列首个阵元视为参考阵元，则在 t 时刻，第 k 个信源使得首个阵元的接收信号为：(2-1)Sk=uk(t)exp(j(ω0t+φ(t))),k=1,2,…,K 其中，j 表示虚部单位，Uk(t) 、 ω0 、 φ(t)表示接收信号的幅度、频率和相位。因此，在 t 时刻，第 m 个阵元的接收信号可表示为：(2-2)xm(t)=∑k=1KSk(t−τm,k)+nm(t),m=1,2,…,M 其中，nm(t)表示第 m个阵元接收到的加性高斯白噪声，τmk)表示对应于第 k 个信号在第 M个阵元的参考接收信号时延。<br/><br/>
由远场窄带信号假设可得：(2-3)Sk(t−τ)≈Sk(t)exp(−jω0τ),m=1,2,…,M 因此，式(2-2)可转化为：(2-4)xm(t)=∑k=1KSk(t)exp(−jω0τm,k)+nm(t),m=1,2,…,M 则，t时刻的阵列输出为：(2-5)x(t)=∑k=1Ka(θk)Sk(t)+n(t)=A(θ)s(t)+n(t) 其中x(t)=[x1(t),…,xN(t)]T∈CN,为阵列输出，A(θ)=[a(θ1),…,a(θK)]∈CN×K,为阵列流形矩阵，s(t)=[s1(t),…,sK(t)]T∈CK ,为入射信号，n(t)∈CN 为阵列接收到的高斯白噪声.式（2-5）一般写作如下(2-6)X=AS+N[1] 当θ∈(−π2,π2),即如图2.1所示建模(2-7)a(θk)=[1,e(−j2πdcos(θk)λ),…,e(−j2π(M−1)dcos(θk)λ)]T[2] 当θ∈(−π,π),即如图2.2所示建模
图2.2 一维线性阵列模型2(2-8)a(θk)=[1,e(−j2πdsin(θk)λ),…,e(−j2π(M−1)dsin(θk)λ)]T2.2 互质阵列模型2.2.1 简单互质阵 考虑两个分别有M和N个传感器的均匀线阵，其中M和N是质数。如图2.3所示，由M个传感器组成的阵列的间距为Nd，而由N个传感器组成的阵列的间距为Md。
图2.3 简单互质阵 假设M<N,且 d=λ2 ，由于质数的特性，上图中由（a）中的两个天线阵组合成的阵列将会共享第一个天线，且后续的天线不会发生覆盖。因此，互质线阵的物理传感器（天线）一共M+N−1个，且最后一个阵元位置在(N−1)Md处。差分阵列可以计算为(2-9)x(m,n)=Mn−Nm 其中0<m<N−1,0<n<M−1,因为M、N互为质数，在MN个不（m,n）组合下，由x(m,n)可以获得MN个不同的差值。因此，仅仅使用M+N−1个实际天线就可以获得具有MN个名义天线的虚拟阵列，其位于(Mn−Nm)d,0<m<N−1,0<n<M−1，这意味着阵列自由度（DOF）可以由O(M+N)增加到O(MN)，更进一步，此处只考虑Mn−Nm,因为还包含Nm−Mn,即−x(m,n),因此可以增加到(2-10)−M(N−1)≤x(m,n)≤M(N−1) 然而，所得到的的虚拟阵列并不是连续的，其存在孔洞（即不连续的部分），如下图2.4中的±11、±8位置。然而，一般常规处理方法是需要获得尽可能长的连续虚拟阵列（这里先不考虑基于差值的凸优化算法），因此需要是从虚拟阵列中获得尽可能长的连续阵列，如下图即是-7~7,那么有没有增长虚拟阵列长度的方法呢？从阵列的角度考虑，就引出了接下来的扩展互质阵。互质阵列如下图所示
图2.4 互质阵及其所得虚拟阵2.2.2 扩展互质阵 为了解决上一节的问题，我们将组成互质阵中的阵元个数为M的子阵阵元数由M调整为2M,以获得一个扩展的互质阵，如图2.5 扩展互质阵所示。
| 图2.5 扩展互质阵 因此得到的差分阵列为(2-11)Sd={x^(m,n) |
x^(m,n)=±(Mn−Nm),0≤m≤2M−1,0≤n≤N−1} 因此在我们得到的差分阵列中，可以获得一个最长连续子阵列范围为−MN到MN。2.2.1 举例SS-MUSIC(以扩展互质阵为例，上一节)（1）接收信号模型 假设有K个远场、不相关的、窄带信源信号作用于扩展互质阵接收天线，其角度，θ=[θ1,θ2，…,θK]T,在时间t时，接收信号向量为(2-12)y(t)=∑k=1Ka(θk)sk(t)+n(t)=A(θ)s(t)+n(t) 其中A(θ)=[a(θ1),a(θ2),…,a(θK)]∈C(2M+N−1)×K,为阵列流形矩阵,，s(t)=[s1(t),s2(t)，…,sK(t)]T∈CK表示信号波形矢量n(t)∽CN(0,δn2I),这里δn2表示噪声功率水平。以L=[l1,l2,…,l2M+N−1]T表示阵元位置，则(2-13)a(θk)=[1,e−j2πl2dsin(θk)λ,…,e−j2πl2M+N−1dsin(θk)λ]接收信号y(t)形成的协方差矩阵可以表示为，其中δk2表示第k个信源的功率。P=diag([δ12,δ22,…,δk2])(2-14)R=E[y(t)y(t)H]=∑k=1Kδk2a(θk)a(θk)H+δk2I=APAH+δk2I考虑到上述的R在实际中无法直接获，一般用采样协方差替代(2-15)R^=1T∑t=1Ty(t)y(t)H=1TYYH（2）SS-MUSIC引入如下两个定义：定义一：虚拟阵元权值函数：令虚拟阵列中第p个虚拟阵元重复的次数为该虚拟阵元的权值函数w(p)，则 w(p)可表示为：w(p)={Pi,Pj∈P |
Pi−Pj=P}定义二：冗余阵元整合矩阵：令J∈C(M+N−1)2×(2MN+1)为冗余阵元整合矩阵，则J中对应的第p个虚拟阵元的行向量可以表示为：(2-16)J(p.:)=vec(I(p))(2-17)I(p)(Pi,Pj)={1w(p),Pi−Pj0,else 对接收协方差矩阵矢量化得到z，从其中选取连续虚拟阵元对应的部分得到v(2-18)z=vec(R^)(2-19)v=Jz此时得到的v∈C(2MN+1)×1,因为其秩为1，因而限制了其可探测的目标数，因此这里采取空间平滑的操作进行秩的恢复，如下图2.6 空间平滑。 |
图2.6 空间平滑计算所的子阵平均协方差Rss，假设取每个字阵长度为MN，则有MN+1个子阵(2-20)Rss=1MN+1∑i=1MN+1vUivUiH最后对获得的Rss应用MUSIC算法，即可。仿真设置信源角度为。。[−50。,30。]，仿真结果如下图
图2.7 仿真结果3.参考文献[1] Zhou, C., et al. (2018). "Direction-of-Arrival Estimation for Coprime Array via Virtual Array Interpolation." IEEE Transactions on Signal Processing 66(22): 5956-5971.
[2]Zhou, C., et al. (2018). "Off-Grid Direction-of-Arrival Estimation Using Coprime Array Interpolation." IEEE Signal Processing Letters 25(11): 1710-1714.
[3] C.-L. Liu, P. P. Vaidyanathan, and P. Pal, “Coprime coarray interpolation
for DOA estimation via nuclear norm minimization,” in Proc. IEEE Int.
Symp. Circuits Syst., Montr´eal, QC, Canada, May 2016, pp. 2639–2642.4.附录（程序）4.1 信号发生程序xxxxxxxxxx261Params.c = 3 * 1e8;2Params.f = 75 * 1e6;3Params.lambda = Params.c / Params.f; %载波波长4Params.fs = 3 * Params.f;5Params.L = 500; %快拍数6Params.snr = 10; %信噪比7Params.t = [1:Params.L] / Params.fs; %时域采样8Params.AntennasNum = 10;9Params.Theta = [-50 30];10Params.Grid = -90:1:90;11% Antennas12M = 3;13N = 5;14Params.Array = Array_f(‘coprime’, 2, M, N, (Params.lambda / 2));15%信号发生，非相干16A = exp(-1 * 1j * 2 * pi * Params.Array.element_positions’ * sin(Params.Theta * (pi / 180)) / Params.lambda); %方向向量17K = length(Params.Theta);18S = randn(K, Params.L) .* exp(1j * 2 * pi * ones(K, 1) * Params.f * Params.t);19Y0 = A * S;20Params.Y = awgn(Y0, Params.snr, ‘measured’); %加入高斯白噪声21Params.R = Params.Y * Params.Y’ / Params.L;22%% 23%最简方法——重排（空间平滑、拓普利兹）24
25% figure(1)26[Pmusic_db, Theta] = Coprime_f(Params, Params.R, M, N);4.2 函数xxxxxxxxxx1051%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2%3% Copyright (c), 2022-2025, xxxx. Co., Ltd.4%5% @file :6%7% @brief :paper “Direction-of-Arrival Estimation for Coprime Array8% via Virtual Array Interpolation”9% a novel virtual array interpolation-based algorithm10% for coprime array11% @version :12%13% @author :hjk14%15% @date :2023.8.2316%17% Details :mainly for coprime array18%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%19function [Pmusic_db, Theta] = Coprime_f(Params, R, M, N)20 z = vec(R);21 %构造去除冗余阵列22 AntennasNum = Params.AntennasNum;23 D = zeros(AntennasNum, AntennasNum);24 d = Params.Array.element_positions / ((Params.lambda / 2));25
26 for m = 1:AntennasNum27 D(m, :) = d - d(m);28 end29
30 D1 = D;31 D = vec(D)’;32 [loc, index] = sort(D);33 indD = unique(D); %去除重复项34 dof = length(indD); %阵列自由度35 %36 J = zeros(max(D) - min(D) + 1, length(z));37 array_num = Params.AntennasNum;38
39 for ik = min(D):max(D)40 II = zeros(array_num, array_num);41
42 for jk = 1:array_num43
44 for kk = 1:array_num45
46 if d(jk) - d(kk) == ik47 II(jk, kk) = 1 / length(find(D == ik));48 end49
50 end51
52 end53
54 J(ik - min(D) + 1, :) = vec(II).’;55 end56
57 v = J * z;58
59 %寻找虚拟阵列连续孔径位置60 cnt = M * N + M - 1;61 flag = (length(v) + 1) / 2;62 pos = flag - cnt:1:flag + cnt;63 %64 v_U = v(pos); %虚拟阵列孔径连续65 %空间平滑法66 subarray_num = M * N + M;67 subarray_num_each = length(v_U) - subarray_num + 1;68
69 Rs = zeros(subarray_num_each, subarray_num_each); % reconstructed covariance matrix70
71 for i = 1:subarray_num72 z_i = v_U(i:i + subarray_num_each - 1);73 Rs = Rs + z_i * z_i’;74 end75
76 % Rf = Rf / subarray_num;77 Params.Array = ula_1d_f(subarray_num_each, (Params.lambda / 2));78 [Pmusic_db, Theta] = music_f(Params, Params.Grid, Rs);79 %80
81 h1 = plot(Params.Grid, Pmusic_db, ‘g’);82 set(gca, ‘XTick’, [-90:30:90])83 xlim([-90 90])84 grid on85 hold on86 %托普利茨构造法87
88 cnt = M * N + M;89 Rf = zeros(cnt, cnt);90
91 for i = 1:1:cnt92 Rf(:, i) = v_U(end - i - cnt + 2:1:end - i + 1);93 end94
95 % Rf = toeplitz(v_U(end - cnt+1:1:end));96
97 Params.Array = ula_1d_f(cnt, (Params.lambda / 2));98 [Pmusic_db, ~] = music_f(Params, Params.Grid, Rf);99 %100 h2 = plot(Params.Grid, Pmusic_db, ‘b*’);101 set(gca, ‘XTick’, [-90:30:90])102 xlim([-90 90])103 grid on104 hold on105end
HJK
Engineer, industry executive, research enthusiast.
Comments