简介一套针对矩形阵列的波束形成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-1n 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, weightsNone): 计算矩形平面阵列的三维方向图因子 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, indexingij) # 角度网格 TH, PH np.meshgrid(theta, phi, indexingij) # 球坐标到方向余弦 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, indexingij) 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, projection3d) surf ax.plot_surface(x, y, z, cmapviridis, linewidth0, antialiasedTrue) 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度 切片phi0 idx_phi0 np.argmin(np.abs(phi - 0)) af_phi0 AF[:, idx_phi0] # 俯仰角90度 切片theta90, 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(phi0 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(theta90 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 影响波束三维图的关键参数速查参数典型值对方向图的影响越界后果阵元间距 dx, dy0.5λ间距越大主瓣越窄超过 λ 会出现栅瓣阵元数 Mx, My8~64阵元越多旁瓣越低波束越窄成本与复杂度上升扫描角范围通常 ±60°大扫描角时主瓣展宽接近 90° 时方向图畸变加权方式均匀/切比雪夫/泰勒控制旁瓣电平与主瓣宽度旁瓣抑制过强导致主瓣变宽工作频率由应用决定λ 改变电尺寸变化频率偏移导致波束指向漂移这张表里的 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, indexingij) 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, at30) My_win windows.chebwin(My, at30) win2d np.outer(Mx_win, My_win) # 注意顺序与矩阵尺寸一致 w steering_weights(16, 16, 0.5*lam, 0.5*lam, theta030, phi045, freq10e9) w_total w * win2d # 重新计算方向图 AF_steered rect_array_factor(theta, phi, Mx, My, dx, dy, freq, weightsw_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, indexingij) 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), at30) win_circ interp1d(dist_ref, win_ref, kindlinear)(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_max90): 检查给定的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, theta030, phi045, freq10e9) # 不加窗直接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]], aspectequal, cmapjet, originlower) plt.colorbar(labeldB) 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 以下就说明实现没出错。这个验证技巧既适用于矩形面阵也可以推广到均匀圆阵和稀疏阵的快速预估。本文还有配套的精品资源点击获取