1. 阵列的方向图

由相同阵元构成的天线阵列，其方向图由两部分相乘得到，第一部分是阵元的方向图，只与阵元本身有关；第二部分取决于阵元间的电流比及相位差，与阵元本身无关，称为阵因子。不妨令阵列的方向图为f(θ,ϕ)f(\theta,\phi)f(θ,ϕ)，则有：
f(θ,ϕ)=f0(θ,ϕ)farr(θ,ϕ)f(\theta,\phi) = f_0(\theta,\phi) f_{arr}(\theta, \phi) f(θ,ϕ)=f0​(θ,ϕ)farr​(θ,ϕ)
其中，f0(θ,ϕ)f_0(\theta,\phi)f0​(θ,ϕ)为阵元的方向性函数（方向图）；farr(θ,ϕ)f_{arr}(\theta, \phi)farr​(θ,ϕ)为阵因子；θ\thetaθ与ϕ\phiϕ的定义如下图所示。

在后文的讨论中，将分析Farr(θ,ϕ)F_{arr}(\theta, \phi)Farr​(θ,ϕ)为归一化的阵因子，其满足Farr(θ,ϕ)=farr(θ,ϕ)max⁡[farr(θ,ϕ)]F_{arr}(\theta, \phi) = \frac{f_{arr}(\theta, \phi) }{ \max{ \left[ f_{arr}(\theta, \phi) \right] } }Farr​(θ,ϕ)=max[farr​(θ,ϕ)]farr​(θ,ϕ)​

2. 理论模型

假设任意阵列中共有MMM个阵元，第mmm个阵元在三维空间中的坐标为pm=[xm,ym,zm]T\boldsymbol{p}_m=\left[ x_m,y_m,z_m \right]^Tpm​=[xm​,ym​,zm​]T。此外，假设在距离原点rrr，方位角为θ\thetaθ，与zzz轴夹角为ϕ\phiϕ的位置存在信号源；信号源向三维空间中辐射的均匀球面波，频率为fff、波长为λ\lambdaλ。显然，可以得到信号源在三维空间中的坐标为ps=[rsin⁡ϕcos⁡θ,rsin⁡ϕsin⁡θ,rcos⁡ϕ]T\boldsymbol{p}_s = \left[ r \sin\phi \cos \theta, r \sin\phi \sin \theta, r \cos \phi \right]^Tps​=[rsinϕcosθ,rsinϕsinθ,rcosϕ]T；则在第mmm个阵元处接收到的信号可以表示为，
s(t,pm)=ej(ωt−km⋅pm);m=0,1,⋯,M−1s(t,\boldsymbol{p}_m) = e^{j \left( \omega t - \boldsymbol{k}_m \cdot \boldsymbol{p}_m \right) };m=0,1,\cdots,M-1 s(t,pm​)=ej(ωt−km​⋅pm​);m=0,1,⋯,M−1
其中，km=2πλpm−ps∥pm−ps∥2\boldsymbol{k}_m=\frac{2 \pi}{\lambda}\frac{\boldsymbol{p}_m - \boldsymbol{p}_s}{\left\Vert \boldsymbol{p}_m - \boldsymbol{p}_s \right\Vert_2}km​=λ2π​∥pm​−ps​∥2​pm​−ps​​为波数，其方向表征了电磁波传输的方向。
在绝大多数应用场景中，为便于分析，往往进行“远场”假设，即，
∀i,j=0,1,⋯,M−1,∥pi−pj∥2≪r\forall i,j=0,1,\cdots,M-1, \left\Vert \boldsymbol{p}_i - \boldsymbol{p}_j \right\Vert_2 \ll r ∀i,j=0,1,⋯,M−1,​pi​−pj​​2​≪r
在这种情况下，从信号源辐射的均匀球面波可以近似地视为均匀平面波，进而在第mmm个阵元处接收到的信号可以简化为：
s(t,pm)=ej(ωt−k0⋅pm);m=0,1,⋯,M−1s(t,\boldsymbol{p}_m) = e^{j \left( \omega t - \boldsymbol{k}_0 \cdot \boldsymbol{p}_m \right) };m=0,1,\cdots,M-1 s(t,pm​)=ej(ωt−k0​⋅pm​);m=0,1,⋯,M−1
其中，若第000个阵元放置在原点，则有，
k0=−2πλps∥ps∥2=−2πλ[sin⁡ϕcos⁡θ,sin⁡ϕsin⁡θ,cos⁡ϕ]T\boldsymbol{k}_0 = - \frac{2 \pi}{\lambda}\frac{\boldsymbol{p}_s}{\left\Vert \boldsymbol{p}_s \right\Vert_2} = - \frac{2 \pi}{\lambda} \left[ \sin\phi \cos \theta, \sin\phi \sin \theta, \cos \phi \right]^T k0​=−λ2π​∥ps​∥2​ps​​=−λ2π​[sinϕcosθ,sinϕsinθ,cosϕ]T

3. 阵列流形与阵因子的计算

值得说明的是，阵列信号处理考察的是在特定的采样时间t0t_0t0​，不同阵元接收到的数据，称之为“快拍”。 假设阵列中各个阵元的激励幅度相同，当空间中存在iii个信号源，且其幅度分别为{Ai∣i=1,2,⋯}\left\{ A_i \left| i=1,2,\cdots \right. \right\}{Ai​∣i=1,2,⋯}时，则在t0t_0t0​时刻，阵列接收的快拍数据可以表示为：
x=ejωt0∑iAia(θi,ϕi)\boldsymbol{x} = e^{j \omega t_0} \sum_i A_i \boldsymbol{a}(\theta_i,\phi_i) x=ejωt0​i∑​Ai​a(θi​,ϕi​)
其中，a(θ,ϕ)=[e−jk0⋅p0,e−jk0⋅p1,⋯,e−jk0⋅pM−1]T\boldsymbol{a}(\theta,\phi) = [e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_0},e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_1},\cdots,e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_{M-1}}]^Ta(θ,ϕ)=[e−jk0​⋅p0​,e−jk0​⋅p1​,⋯,e−jk0​⋅pM−1​]T即为阵列流形。同时，阵因子可以通过阵列流形计算获得，即
farr(θ,ϕ)=∥a(θ,ϕ)∥1=∑m=0M−1e−jk0⋅pm\begin{aligned} f_{arr}(\theta, \phi) &= \Vert \boldsymbol{a}(\theta,\phi) \Vert_1 \\ &=\sum_{m=0}^{M-1} e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_m} \end{aligned} farr​(θ,ϕ)​=∥a(θ,ϕ)∥1​=m=0∑M−1​e−jk0​⋅pm​​
值得说明的是，本文假设阵列中各个阵元的激励幅度相同；在激励幅度不同的情况下，可以在a(θ,ϕ)\boldsymbol{a}(\theta, \phi)a(θ,ϕ)中的各个元素前增加权重进行计算。

4. 典型阵列的数字波束形成与MATLAB仿真

数字波束形成的关键在于正确地设置阵列流形，前文我们阐述了其一般形式，即
a(θ,ϕ)=[e−jk0⋅p0,e−jk0⋅p1,⋯,e−jk0⋅pM−1]T\boldsymbol{a}(\theta,\phi) = [e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_0},e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_1},\cdots,e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_{M-1}}]^T a(θ,ϕ)=[e−jk0​⋅p0​,e−jk0​⋅p1​,⋯,e−jk0​⋅pM−1​]T
换言之，阵列流形中的第mmm个元素可以表示为：
[a(θ,ϕ)]m=e−jk0⋅pm=ej2π(xmsin⁡ϕcos⁡θ+ymsin⁡ϕsin⁡θ+zmcos⁡ϕ)λ\left[ \boldsymbol{a}(\theta,\phi) \right]_m = e^{-j\boldsymbol{k}_0 \cdot \boldsymbol{p}_m} = e^{\frac{j2 \pi \left( x_m \sin\phi \cos \theta + y_m \sin\phi \sin \theta + z_m \cos \phi \right) }{\lambda}} [a(θ,ϕ)]m​=e−jk0​⋅pm​=eλj2π(xm​sinϕcosθ+ym​sinϕsinθ+zm​cosϕ)​
依据上式，可以对任意阵列进行波束形成，相关的性能参数由阵列的排布决定。

下面，本节将结合具体的阵列形式，给出仿真MATLAB仿真结果。在仿真中采用相同的信号源，该信号源向空间中辐射高频段的电磁波，其频率f=15MHzf=15MHzf=15MHz，波长λ=20m\lambda = 20mλ=20m；在三维空间中，其矢径为r=1000kmr=1000kmr=1000km，方位角为ϕ=45∘\phi=45^{\circ}ϕ=45∘，仰角为90∘−θ=45∘90^{\circ}-\theta=45^{\circ}90∘−θ=45∘。

4.1 均匀直线阵列

不妨令阵元数量M=30M = 30M=30，阵元间距d=8md = 8md=8m，各个阵元沿xxx轴分布，则第mmm个阵元的坐标可以表示为：
pm=[md,0,0]T;m=0,1,⋯,M−1\boldsymbol{p}_m=\left[ md,0,0 \right]^T;m=0,1,\cdots,M-1 pm​=[md,0,0]T;m=0,1,⋯,M−1
因此，均匀直线阵列的阵列流形可以表示为：
a(θ,ϕ)=[1,ej2πdsin⁡θcos⁡ϕλ,⋯,ej2π(M−1)dsin⁡θcos⁡ϕλ]T\boldsymbol{a}(\theta,\phi) = \left[1,e^{ \frac{j 2 \pi d \sin \theta \cos \phi}{\lambda} },\cdots,e^{ \frac{j 2 \pi (M - 1) d \sin \theta \cos \phi}{\lambda} } \right]^T a(θ,ϕ)=[1,eλj2πdsinθcosϕ​,⋯,eλj2π(M−1)dsinθcosϕ​]T

(1) 本文在三维坐标系中推导均匀直线阵列的阵列流形；而现有文献大多在二维坐标系中推导均匀直线阵列的阵列流形，故不需要考虑仰角，因而在大多数资料中没有cos⁡ϕ\cos \phicosϕ这一项；

(2) 均匀直线阵列仅有一维的角度分辨能力，在三维空间中，为了获得正确的方位角，需要利用cos⁡ϕ\cos \phicosϕ将二维平面从xOyxOyxOy平面转化到“信号源与均匀直线阵列各个阵元所在的平面”。反之，若在阵列流形中没有设置正确的仰角，将会导致方位角的估计错误；

(3) 不难发现cos⁡ϕ=cos⁡(−ϕ)\cos \phi = \cos (-\phi)cosϕ=cos(−ϕ)，故均匀直线阵列会出现测角模糊的问题，其不模糊的测角范围是[0∘,180∘]\left[ 0^{\circ}, 180^{\circ} \right][0∘,180∘]。

(a) 从均匀直线阵列的仿真结果可知，当在仰角为45∘45^\circ45∘处进行波束形成时，可以在正确的方位上（45∘45^\circ45∘）获得峰值；当在仰角为0∘0^\circ0∘处进行波束形成时，估计的方位为60∘60^\circ60∘，与结果显然不符；进而说明在三维空间中，即便没有仰角分辨能力，均匀直线阵列的阵列流形仍然需要考虑仰角。
(b) 在均匀直线阵列的仿真结果中，在仰角维度的波瓣极宽，说明均匀直线阵列没有仰角分辨能力。

4.2 均匀矩形阵列

假设均匀矩形阵列共有Mx=30M_x = 30Mx​=30行，My=30M_y =30My​=30列阵元；在每一行，阵元间距dx=8md_x = 8mdx​=8m；在每一列，阵元间距dy=8md_y = 8mdy​=8m，则第iii行、第jjj列的阵元的坐标可以表示为：
pij=[idx,jdy,0]Ti=0,1,⋯,Mx−1;j=0,1,⋯,My−1\begin{aligned} & \boldsymbol{p}_{ij}=\left[ id_x,jd_y,0 \right]^T \\ & i=0,1,\cdots,M_x-1;j=0,1,\cdots,M_y-1 \end{aligned} ​pij​=[idx​,jdy​,0]Ti=0,1,⋯,Mx​−1;j=0,1,⋯,My​−1​
若按j⊗ij \otimes ij⊗i的形式向量化阵列排布，则均匀矩形阵列的阵列流形可以表示为：
a(θ,ϕ)=ay(θ,ϕ)⊗ax(θ,ϕ)\boldsymbol{a}(\theta,\phi) = \boldsymbol{a}_y(\theta,\phi) \otimes \boldsymbol{a}_x(\theta,\phi) a(θ,ϕ)=ay​(θ,ϕ)⊗ax​(θ,ϕ)
其中，⊗\otimes⊗为克罗内克积，
ax(θ,ϕ)=[1,ej2πdxsin⁡θcos⁡ϕλ,⋯,ej2π(Mx−1)dxsin⁡θcos⁡ϕλ]Tay(θ,ϕ)=[1,ej2πdysin⁡θsin⁡ϕλ,⋯,ej2π(My−1)dysin⁡θsin⁡ϕλ]T\begin{aligned} \boldsymbol{a}_x(\theta,\phi) &= \left[1,e^{ \frac{j 2 \pi d_x \sin \theta \cos \phi}{\lambda} },\cdots,e^{ \frac{j 2 \pi (M_x - 1) d_x \sin \theta \cos \phi}{\lambda} } \right]^T \\ \boldsymbol{a}_y(\theta,\phi) &= \left[1,e^{ \frac{j 2 \pi d_y \sin \theta \sin \phi}{\lambda} },\cdots,e^{ \frac{j 2 \pi (M_y - 1) d_y \sin \theta \sin \phi}{\lambda} } \right]^T \end{aligned} ax​(θ,ϕ)ay​(θ,ϕ)​=[1,eλj2πdx​sinθcosϕ​,⋯,eλj2π(Mx​−1)dx​sinθcosϕ​]T=[1,eλj2πdy​sinθsinϕ​,⋯,eλj2π(My​−1)dy​sinθsinϕ​]T​

(a) 不难发现，均匀矩形阵列具有方位角、仰角的二维分辨能力，且均匀矩形阵列在[0∘,360∘][0^\circ,360^\circ][0∘,360∘]的方位角范围内，不会出现测角模糊。

4.3 均匀圆形阵列

假设均匀圆形阵列共有M=96M=96M=96个阵元，阵列半径为R=120mR = 120mR=120m，阵列圆心位于原点，则第mmm个阵元的坐标可以表示为：
pm=[Rcos⁡2πmM,Rcos⁡2πmM,0]T;m=0,1,⋯,M−1\boldsymbol{p}_m = \left[ R \cos{\frac{2 \pi m}{M}}, R \cos{\frac{2 \pi m}{M}}, 0 \right]^T;m=0,1,\cdots,M-1 pm​=[RcosM2πm​,RcosM2πm​,0]T;m=0,1,⋯,M−1
因此，均匀圆形阵列的阵列流形可以表示为：
a(θ,ϕ)=[1,ej2πdR(cos⁡2πMsin⁡ϕcos⁡θ+cos⁡2πMsin⁡ϕsin⁡θ)λ,⋯,ej2πdR[cos⁡2π(M−1)Msin⁡ϕcos⁡θ+cos⁡2π(M−1)Msin⁡ϕsin⁡θ]λ]T\boldsymbol{a}(\theta,\phi) = \left[1,e^{ \frac{j 2 \pi d R \left( \cos{\frac{2 \pi}{M}} \sin\phi \cos \theta + \cos{\frac{2 \pi}{M}} \sin\phi \sin \theta \right)}{\lambda} },\cdots,e^{ \frac{j 2 \pi d R \left[ \cos{\frac{2 \pi (M-1) }{M}} \sin\phi \cos \theta + \cos{\frac{2 \pi (M-1)}{M}} \sin\phi \sin \theta \right]}{\lambda} } \right]^T a(θ,ϕ)=[1,eλj2πdR(cosM2π​sinϕcosθ+cosM2π​sinϕsinθ)​,⋯,eλj2πdR[cosM2π(M−1)​sinϕcosθ+cosM2π(M−1)​sinϕsinθ]​]T

(a) 不难发现，均匀圆形阵列具有方位角、仰角的二维分辨能力，且均匀矩形阵列在[0∘,360∘][0^\circ,360^\circ][0∘,360∘]的方位角范围内，不会出现测角模糊。

5. MATLAB仿真程序源码

%% Author：ZLY
clearvars;
close all;
clc;
%% 信号为单频正弦信号
global c;
global f0;
global lambda;
global t;
global N;
global TxLoc;c = 3e8;
f0 = 15e6;              % 载波频率
lambda = c / f0;    % 波长
fs = 4 * f0;              % 采样频率
N = 1024;               % 采样点数
t = (0:N-1).' / fs;    % 采样时刻x = exp(1j * 2 * pi * f0 * t);  % 单频信号
%% 信号源位于远场
% 1、阵元间距 < lambda / 2
% 2、阵列孔径 v.s. 角度分辨力
% TxLoc = [1000e3, 1000e3, 245.5756e3];                       % 远场源，以正北为参照方位角45°，仰角10°（常见的情况）
TxLoc = [1000e3, 1000e3, sqrt(2) * 1000e3];           % 远场源，以正北为参照方位角45°，仰角45°（夸张的情况）
arrayType = 3;
switch arrayTypecase 1ULA();case 2URA();case 3UCA();
end
%% 均匀线阵的情况
function ULA()global c;global N;global f0;global lambda;global t;global TxLoc;d = 8;M = 30;RxLoc = [ d * (0:M-1).', zeros(M, 2)];% 生成接收信号，N个快拍的数据y = zeros(N, M);for ii = 1:Mtau = norm( RxLoc(ii,:) - TxLoc ) / c;y(:,ii) = exp(1j * 2 * pi * f0 * (t - tau));end% 波束形成k = 2 * pi / lambda;phi = ( 0 : 359 ) / 180 * pi;            % 与x轴夹角theta = ( 90:-1:0 ) / 180 * pi;         % 与z轴夹角DBF = zeros(length(phi), length(theta));for ii = 1:length(phi)for jj = 1:length(theta)a = exp( 1j * k * (0 : M - 1) * d * cos ( phi(ii) ) * sin( theta(jj) ) ).';DBF(ii,jj) = sum( abs( y * conj(a) ), 1 );endendDBF = DBF ./ max(abs(DBF), [], 'all');figure(1);mesh(90 - theta / pi * 180, phi / pi * 180, db(DBF));xlabel('俯仰角'); ylabel('方位角'); colormap jet; colorbar; view(2);axis tight; set(gca, 'CLim', [-40, 0]);title('ULA');end%% 均匀矩形阵列的情况
function URA()global c;global N;global f0;global lambda;global t;global TxLoc;dx = 8;Mx = 30;dy = 8;My = 30;RxLoc = zeros(Mx * My, 3);for ii = 1:Mxfor jj = 1:MyRxLoc( (ii - 1) * My + jj, :) = [(ii - 1) * dx, (jj - 1) * dy, 0];endend% 生成接收信号，N个快拍的数据y = zeros(N, Mx * My);for ii = 1:Mx * Mytau = norm( RxLoc(ii,:) - TxLoc ) / c;y(:,ii) = exp(1j * 2 * pi * f0 * (t - tau));end% 波束形成k = 2 * pi / lambda; phi = ( 0 : 359 ) / 180 * pi;            % 与x轴夹角theta = ( 90:-1:0 ) / 180 * pi;         % 与z轴夹角DBF = zeros(length(phi), length(theta));for ii = 1:length(phi)for jj = 1:length(theta)ax = exp( 1j * k * (0 : Mx - 1) * dx * cos( phi(ii) ) * sin( theta(jj) ) ).';ay = exp( 1j * k * (0 : My - 1) * dy * sin( phi(ii) ) * sin( theta(jj) ) ).';a = kron( ay, ax );DBF(ii,jj) = sum( abs( y * conj(a) ), 1 );endendDBF = DBF ./ max(abs(DBF), [], 'all');figure(2);mesh(90 - theta / pi * 180, phi / pi * 180, db(DBF));xlabel('俯仰角'); ylabel('方位角'); colormap jet; colorbar; view(2);axis tight; set(gca, 'CLim', [-40, 0]);title('URA');end%% 均匀圆形阵列
function UCA()global c;global N;global f0;global lambda;global t;global TxLoc;M = 96;     % 96个阵元R = 120;     % 半径120mRxLoc = zeros(M, 3);for ii = 1:Malpha = (ii - 1) * 2 * pi / M ;RxLoc(ii,:) = [R * cos(alpha), R * sin(alpha), 0];end% 生成接收信号，N个快拍的数据y = zeros(N, M);for ii = 1:Mtau = norm( RxLoc(ii,:) - TxLoc ) / c;y(:,ii) = exp(1j * 2 * pi * f0 * (t - tau));end% 波束形成k = 2 * pi / lambda; phi = ( 0 : 359 ) / 180 * pi;            % 与x轴夹角theta = ( 90:-1:0 ) / 180 * pi;         % 与z轴夹角DBF = zeros(length(phi), length(theta));for ii = 1:length(phi)for jj = 1:length(theta)xm = R * cos( (0: M - 1) * 2 * pi / M ).';ym = R * sin( (0: M - 1) * 2 * pi / M ).';a = exp(1j * k * ( xm * cos(phi(ii)) * sin(theta(jj)) + ym * sin(phi(ii)) * sin(theta(jj)) ) );DBF(ii,jj) = sum( abs( y * conj(a) ), 1 );endendDBF = DBF ./ max(abs(DBF), [], 'all');figure(3);mesh(90 - theta / pi * 180, phi / pi * 180, db(DBF));xlabel('俯仰角'); ylabel('方位角'); colormap jet; colorbar; view(2);axis tight; set(gca, 'CLim', [-40, 0]);title('UCA');end