矩形阵列三维波束形成：Python方向图绘制与FFT验证

简介：一套针对矩形阵列的波束形成MATLAB代码包，面向信号处理与通信工程领域的学生和研究人员，核心覆盖三维波束形成与平面波束形成，并包含波束三维图的绘制方法。资源共9个文件，其中7个.m脚本为主程序，实现坐标变换、等高线绘制、多角度方向图分析、三维波束可视化等；2个.asv为MATLAB自动备份文件，便于恢复历史版本。压缩包仅6KB，轻量易用。通过修改脚本中的阵列参数、入射角度等，可直观观察不同条件下的波束指向和旁瓣特性，有助于深入理解矩形阵列的波束形成原理，也便于在此基础上开展扩展实验。在雷达、无线通信、声呐等场景中，波束形成可显著提升方向性增益与抗干扰能力，这套代码为此类研究提供了便捷的仿真起点。已有342人学习下载，适合作为课程设计或科研入门的参考。

1. 矩形阵列三维波束形成：为什么从线阵换到面阵

当阵列从一条线变成平面，波束形成的维度从一维跳到了二维。线性阵列只能控制一个维度的指向角，而矩形平面阵列通过将阵元分布在二维平面上，可以同时调节方位角与俯仰角，形成真正意义上的三维波束。所谓“波束三维图”，就是把方向图幅度响应按球坐标画成曲面，主瓣像一个锥形山丘，栅瓣和旁瓣则是周围的小突起。这个图不是炫技用的，它直接告诉你阵列在哪个角度能收到信号、哪个角度会泄露能量。

做雷达、声呐、无线通信大规模天线阵的人，几乎每天都要和这类图打交道。区别只是频率段不同——从几百兆赫兹到几十吉赫兹，再到超声和水声，几何关系完全相同。本文围绕矩形阵列的平面波束形成，把方向图推导、三维图绘制、权重设置和栅瓣验证串起来，给出可以直接在本地运行的 Python 代码和参数边界。新手能照着画出第一张 3D 方向图，老手可以复用后半部分的 FFT 验证思路来排查波束指向偏移和栅瓣问题。

2. 矩形平面阵列的阵列流形与方向图推导

2.1 从线阵到面阵：相位差如何累加

线性阵列中，相邻阵元的相位差由投影距离决定。设波达方向与阵列轴向夹角为 θ，阵元间距为 d，则相邻阵元接收信号的相位差为 (2\pi d \cos\theta / λ)。矩形阵列把这条线“扩展”成网格：阵元在 x 轴方向间距为 dx，在 y 轴方向间距为 dy，共 Mx 行、My 列。任意阵元的位置是 (m·dx, n·dy)，m = 0,1,…,Mx-1，n = 0,1,…,My-1。

当平面波从方向 (θ, φ) 入射时，这里 θ 是俯仰角（从 z 轴正方向算起，0° 表示阵列正面法向），φ 是方位角（从 x 轴正方向算起），波程差投影到 x 轴的量是 dx·m·sinθ·cosφ，投影到 y 轴的量是 dy·n·sinθ·sinφ。因此第 (m,n) 个阵元相对于原点阵元的相位延迟为：

[ \Delta \psi(m,n) = \frac{2\pi}{\lambda} (m dx \sin\theta \cos\phi + n dy \sin\theta \sin\phi) ]

注意这里的方向约定：如果我们要让波束指向 (θ0, φ0)，就需要在加权时补偿这个相位，也就是取负指数。而这个相位表达式同时出现在阵列流形向量和导向向量中，是后续所有计算的基础。

2.2 矩形阵列方向图的向量化计算

工程实现时，最直接的做法是把每个阵元的坐标铺成网格，然后对一组待观察角度计算“阵列流形向量”，再与权向量求内积。下面这段 Python 代码使用 NumPy 的广播机制，不需要循环遍历每个阵元，速度足够应对 128×128 阵列的常规仿真。

import numpy as np def rect_array_factor(theta, phi, Mx, My, dx, dy, freq, weights=None): """ 计算矩形平面阵列的三维方向图因子 theta: 俯仰角数组，[deg]，从z轴算起 phi: 方位角数组，[deg]，从x轴算起 Mx: x方向阵元数 My: y方向阵元数 dx, dy: 阵元间距，单位m freq: 工作频率，Hz weights: 权向量，shape (Mx, My)，默认全1 """ c = 299792458.0 lam = c / freq k = 2 * np.pi / lam # 阵元坐标网格 m = np.arange(Mx) - (Mx - 1) / 2.0 n = np.arange(My) - (My - 1) / 2.0 m_grid, n_grid = np.meshgrid(m, n, indexing='ij') # 角度网格 TH, PH = np.meshgrid(theta, phi, indexing='ij') # 球坐标到方向余弦 u = np.sin(np.radians(TH)) * np.cos(np.radians(PH)) v = np.sin(np.radians(TH)) * np.sin(np.radians(PH)) # 相位项 exp(j*k*(x*u + y*v))，注意是正向传播相位 # 方向图因子 = sum(a(m,n) * exp(j*k*(m*dx*u + n*dy*v))) phase = k * (m_grid[:, :, np.newaxis, np.newaxis] * dx * u + n_grid[:, :, np.newaxis, np.newaxis] * dy * v) # 广播后形状: (Mx, My, Ntheta, Nphi) AF = np.sum(np.exp(1j * phase), axis=(0, 1)) if weights is not None: # 权向量逐阵元相乘后累加 AF = np.sum(weights[:, :, np.newaxis, np.newaxis] * np.exp(1j * phase), axis=(0, 1)) return np.abs(AF) / np.max(np.abs(AF))

这段代码的关键在于把相位拆成 x、y 两项，再用np.newaxis升维做广播。AF的形状是(Ntheta, Nphi)，对应每个俯仰角和方位角组合下的幅度响应。默认全 1 权重对应均匀加权，此时方向图就是阵列因子本身。如果需要后续做切比雪夫窗或泰勒窗，则把窗函数矩阵传给weights。使用前要确认角度定义与目标一致：本函数采用“俯仰角从 z 轴算起”的物理学约定，很多天线教材用“仰角从 xy 平面算起”，差 90°，换算时不注意会在主瓣位置上产生明显的偏移。

3. 用 Python 绘制波束三维图：从方向图到可视化

3.1 准备方向图数据并映射到三维坐标

有了方向图因子之后，画三维图前还需要做一步映射：把幅度响应从 (θ, φ) 网格转换成笛卡尔坐标。常见做法是用球坐标转直角坐标，半径方向表示归一化幅度。这样画出来的曲面，每个点离原点的距离就是这个角度下阵列响应的幅度。主瓣会鼓出来一个大包，零点则是向原点凹进去的深坑。

下面的代码承接上一节的rect_array_factor，生成 121×241 的角度网格，并计算 10GHz 下 16×16 阵列的方向图。然后使用matplotlib的plot_surface绘制三维曲面。

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 参数设置 freq = 10e9 # 10 GHz c = 299792458.0 lam = c / freq Mx, My = 16, 16 dx = dy = 0.5 * lam # 半波长间距 # 角度扫描范围：俯仰0~90度，方位0~360度 theta = np.linspace(0, 90, 121) phi = np.linspace(0, 360, 241) # 计算方向图 AF = rect_array_factor(theta, phi, Mx, My, dx, dy, freq) # 网格化 TH, PH = np.meshgrid(theta, phi, indexing='ij') th_r = np.radians(TH) ph_r = np.radians(PH) # 球坐标转直角坐标：r = 归一化幅度，角度按标准球坐标 r = AF x = r * np.sin(th_r) * np.cos(ph_r) y = r * np.sin(th_r) * np.sin(ph_r) z = r * np.cos(th_r) fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') surf = ax.plot_surface(x, y, z, cmap='viridis', linewidth=0, antialiased=True) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_title('16x16 Rectangular Array 3D Beam Pattern') plt.tight_layout() plt.show()

这段代码里，r = AF等价于把归一化方向图当作球半径。注意当 θ=0° 时，无论 φ 取多少，sinθ=0，所有坐标都汇聚到同一个点，因此目视图中央会出现一个尖点，这是正确的。若想观察旁瓣的细节，把幅度做分贝压缩更合适，比如r_db = 10**(AF_db/20)，但栅瓣和主瓣的动态范围很大，线性图通常会掩盖旁瓣。

3.2 用等高线图和切片图辅助分析

三维曲面图适合看整体形状，却很难准确判断主瓣宽度和零点位置。建议同时绘制两个二维切片：方位角固定时俯仰方向的方向图，俯仰角固定时方位方向的方向图。下面的代码从 AF 数据中切出 φ=0° 和 θ=90° 两条线，注意 θ=90° 对应 xy 平面，此时方向图退化为线阵的方向图。

# 方位角=0度 切片（phi=0） idx_phi0 = np.argmin(np.abs(phi - 0)) af_phi0 = AF[:, idx_phi0] # 俯仰角=90度 切片（theta=90, xy平面） idx_th90 = np.argmin(np.abs(theta - 90)) af_th90 = AF[idx_th90, :] fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(theta, 20*np.log10(af_phi0 + 1e-12)) axes[0].set_title('phi=0 deg, theta scan') axes[0].set_xlabel('theta [deg]') axes[0].set_ylabel('dB') axes[0].grid(True) axes[1].plot(phi, 20*np.log10(af_th90 + 1e-12)) axes[1].set_title('theta=90 deg, phi scan') axes[1].set_xlabel('phi [deg]') axes[1].set_ylabel('dB') axes[1].grid(True) plt.tight_layout() plt.show()

切片图中的分贝转换要注意加一个小量避免log10(0)。从图上能读出主瓣半功率宽度、第一旁瓣高度等指标。均匀加权时矩形阵列的第一旁瓣电平约为 −13.3dB，和线阵一致，这是因为阵列因子在方向余弦域里是矩形窗的二维傅里叶变换，旁瓣特性由窗函数决定，而这里默认一致加权。

3.3 影响波束三维图的关键参数速查

这张表里的 0.5λ 是硬约束但不是绝对边界。在相控阵中，为了避免在可视区内出现栅瓣，最大扫描角 θmax 时需满足 (d / λ \le 1 / (1 + \sin\theta_{max}))。如果只扫描小角度，间距可以适当放宽；如果要求全空间扫描，0.5λ 几乎不能动。这一点在后面的栅瓣验证中会再次用到。

4. 波束形成的权重设计：指向、旁瓣抑制与栅瓣规避

4.1 常规波束形成的权向量计算

要让波束指向 (θ0, φ0)，权向量必须补偿从原点到各阵元的传播延迟，也就是取相位共轭。在上一节的符号约定下，权向量为：

[ w(m,n) = e^{-j k (m dx \sin\theta_0 \cos\phi_0 + n dy \sin\theta_0 \sin\phi_0)} ]

考虑到阵元坐标以阵列中心为零点，这里的 m、n 可能为负值。实现时直接生成与阵元位置对应的权矩阵。下面给出加入切比雪夫窗的完整示例。切比雪夫窗在等旁瓣电平下可使所有旁瓣高度一致，但二维扩展时需要分别沿 x、y 方向生成窗函数，再取外积。

from scipy.signal import windows def steering_weights(Mx, My, dx, dy, theta0, phi0, freq): c = 299792458.0 lam = c / freq k = 2 * np.pi / lam m = np.arange(Mx) - (Mx - 1) / 2.0 n = np.arange(My) - (My - 1) / 2.0 m_grid, n_grid = np.meshgrid(m, n, indexing='ij') u0 = np.sin(np.radians(theta0)) * np.cos(np.radians(phi0)) v0 = np.sin(np.radians(theta0)) * np.sin(np.radians(phi0)) w = np.exp(-1j * k * (m_grid * dx * u0 + n_grid * dy * v0)) return w # 生成切比雪夫窗，旁瓣电平-30dB Mx_win = windows.chebwin(Mx, at=30) My_win = windows.chebwin(My, at=30) win2d = np.outer(Mx_win, My_win) # 注意顺序与矩阵尺寸一致 w = steering_weights(16, 16, 0.5*lam, 0.5*lam, theta0=30, phi0=45, freq=10e9) w_total = w * win2d # 重新计算方向图 AF_steered = rect_array_factor(theta, phi, Mx, My, dx, dy, freq, weights=w_total)

权向量中的负号与方向图计算中的正号互为共轭，目的是让阵元内积在期望方向同相叠加。窗函数外积生成二维窗时，要注意np.outer的行列顺序——第一个参数对应 x 方向还是 y 方向，取决于你后续如何把矩阵展开成向量。通常建议先把 x 方向权重放在行方向，与坐标网格一致，避免转置错误。实际运算时，weights矩阵通常需要展成一维向量再做内积，但上面的rect_array_factor里直接广播成四维，省去了 reshape 的麻烦。

4.2 窗函数在二维阵面的扩展方式

一维窗函数可以直接乘以线阵激励，但二维矩形阵存在两种常见扩展：可分离窗和圆对称窗。可分离窗就是两个一维窗的外积，优点是计算简单、实现方便，但会在对角线方向产生比主轴更高的旁瓣。圆对称窗则根据阵元到中心的距离计算窗值，如二维泰勒窗或圆孔径切比雪夫窗，能更好地抑制斜向旁瓣，但计算量更大。工程上，如果只关心主平面（φ=0° 和 φ=90°）的性能，可分离窗足够；如果要求全空间低旁瓣，建议用圆对称窗。

# 计算每个阵元到中心的距离（单位：波长） x_pos = (np.arange(Mx) - (Mx-1)/2) * dx y_pos = (np.arange(My) - (My-1)/2) * dy Xp, Yp = np.meshgrid(x_pos, y_pos, indexing='ij') dist_lambda = np.sqrt((Xp/lam)**2 + (Yp/lam)**2) # 使用一维切比雪夫窗按距离插值近似圆对称窗 # 实际实现通常用二维窗函数库，这里示意用法 from scipy.interpolate import interp1d dist_ref = np.linspace(0, dist_lambda.max(), 1024) win_ref = windows.chebwin(len(dist_ref), at=30) win_circ = interp1d(dist_ref, win_ref, kind='linear')(dist_lambda)

这里的圆对称窗用一维窗按半径映射，近似效果取决于阵面形状。对于矩形阵列，严格圆对称窗在边缘处截断会产生新的旁瓣，实际效果需用方向图迭代验证。若你的系统对旁瓣有硬指标，建议直接使用已有天线工具包中的二维窗函数，或采用“采样密度加权”方法。

4.3 栅瓣出现条件与阵元间距选择

栅瓣是阵列方向图中周期性重复的主瓣副本，来源于空间混叠。当阵元间距大于 λ 时，相邻阵元的相位差在某个方向上刚好等于 2π 的整数倍，原本不该出现的方向也会满足同相叠加条件。栅瓣位置由方程 (kd(\sin\theta\cos\phi - \sin\theta_0\cos\phi_0) = 2\pi p) 决定。下面给出一个快速判断脚本，当扫描角为 (θ0, φ0) 时，计算可视区内是否出现栅瓣。

def check_grating_lobe(dx, dy, theta0, phi0, theta_max=90): """ 检查给定的dx/dy和扫描角是否产生栅瓣 返回可视区内是否存在栅瓣（布尔值） """ lam = 1.0 # 归一化，使用波长单位 k = 2 * np.pi / lam u0 = np.sin(np.radians(theta0)) * np.cos(np.radians(phi0)) v0 = np.sin(np.radians(theta0)) * np.sin(np.radians(phi0)) # 遍历所有可能的整数对 (px, py) for px in range(-5, 6): for py in range(-5, 6): if px == 0 and py == 0: continue # 栅瓣条件，解出u和v delta_u = px * lam / dx delta_v = py * lam / dy u_gl = u0 + delta_u v_gl = v0 + delta_v if u_gl**2 + v_gl**2 <= 1: # 在可视圆锥内 return True return False # 示例：0.6λ 间距，扫描到60度 print(check_grating_lobe(0.6, 0.6, 60, 0))

上面脚本用归一化波长，省去频率换算的麻烦。实际上，只要可视区不包含任何非零整数对对应的投影坐标，就不会有栅瓣。对于矩形阵列，这个条件可以简化为 (d_x / \lambda \le 1/(1+\sin\theta_{max})) 同时满足 x 和 y 两个方向。上面代码中的u_gl**2 + v_gl**2 <= 1就是判断该方向是否落在实空间内（方向余弦平方和小于等于1）。若返回 True，就说明该间距在对应扫描角下会看到栅瓣。

5. 用二维 FFT 快速验证波束三维图的正确性

前面写的方向图函数逐角度计算相位累加，准确但速度慢。调试阶段更推荐用二维 FFT 来验证：把阵面激励当作一个二维数组，对激励做二维傅里叶变换，其结果就是方向图在方向余弦平面上的采样。这个方法的输出可以直接和rect_array_factor的结果对比，确认是否出现权重矩阵转置、相位符号搞反等问题。

作法很简单：创建一个与阵面大小相同的复数激励矩阵，阵元值等于权向量乘上加窗系数，然后对矩阵做np.fft.fft2，再用fftshift把零频移到中心。横纵轴分别对应方向余弦 u 和 v，每个像素点的物理角度可以通过 (u = \lambda / (Mx dx) * ix) 换算。注意 FFT 默认索引顺序为 [行][列]，与我们的 x、y 网格存在对应关系，需要确认矩阵第一维对应 y 还是 x。下面的代码演示了 32×32 阵列在 10GHz 下的 FFT 方向图切片。

Mx = 32 My = 32 dx = dy = 0.5 * lam w = steering_weights(Mx, My, dx, dy, theta0=30, phi0=45, freq=10e9) # 不加窗，直接FFT AF_fft = np.fft.fftshift(np.fft.fft2(w)) AF_db = 20 * np.log10(np.abs(AF_fft) / np.max(np.abs(AF_fft)) + 1e-12) # 建立方向余弦坐标 ix = np.arange(Mx) - Mx // 2 iy = np.arange(My) - My // 2 u_axis = ix / (Mx * dx / lam) # u 坐标，单位为1 v_axis = iy / (My * dy / lam) plt.imshow(AF_db, extent=[u_axis[0], u_axis[-1], v_axis[0], v_axis[-1]], aspect='equal', cmap='jet', origin='lower') plt.colorbar(label='dB') plt.xlabel('u = sin(theta)cos(phi)') plt.ylabel('v = sin(theta)sin(phi)') plt.title('2D FFT of 2D Array Weights') plt.show()

运行后，主瓣中心应当出现在 (u0, v0) = (sin30°cos45°, sin30°sin45°) ≈ (0.3536, 0.3536) 的位置。如果你发现主瓣在镜像位置，说明权向量取相位时符号搞反了；如果 FFT 图比逐项计算的方向图多出一圈重复图案，那就是 0.5λ 间距在 FFT 中周期性延拓的正常现象，并不是栅瓣——FFT 的输出本身是周期的，只有落在 unit circle 内的部分才是真实可视区。对比时建议把 FFT 的结果按方向余弦映射到球坐标，再和方向图函数的输出做归一化均方误差，差值在 −80dB 以下就说明实现没出错。这个验证技巧既适用于矩形面阵，也可以推广到均匀圆阵和稀疏阵的快速预估。

本文还有配套的精品资源，点击获取