简介:这份资源是面向天线设计与电磁仿真方向学习者的Matlab阵列天线仿真脚本,适合通信工程、电磁场与微波技术等专业学生及工程师用于波束形成与方向图分析实战。脚本可用于交互式设置阵元数量、间距、激励信号等参数,自动计算并绘制天线方向图,输出波束形成性能指标,帮助理解主瓣宽度、副瓣电平及相位幅度控制等核心概念。压缩包内仅包含一个m格式脚本文件,整体大小约620B,结构精简却完整覆盖了阵列天线设计的关键流程。目前已有279人学习下载,适合作为中高级学习者快速验证阵列天线算法的参考工具。通过运行该脚本,可直观对比不同阵列参数对辐射效果的影响,为后续深入掌握Matlab Phased Array System Toolbox等工具箱的应用打下基础,也可在此基础上二次开发扩展功能。
1. 从 gbo.rar 说起:这套阵列天线仿真代码到底能帮你做什么
做阵列天线和天线方向图仿真的人,大概率都在某个资料包里见过“gbo.rar”这类名字的压缩包,解压后里面是一堆 matlab 脚本,跑出来一张极坐标方向图,主瓣、副瓣、零陷都有,但很少有人一次跑明白:波束指向为什么偏了、副瓣为什么压不下去、换个频率怎么方向图就乱了。这背后其实是一套完整的阵列天线仿真链路——从阵元排布到方向图计算,再到波束成形和优化。这个标题指向的,就是用 matlab 做阵列天线方向图仿真、天线波束控制、以及用 GBO(Gradient-based Optimizer,梯度优化器)这类算法做阵列综合的一套可落地脚本方案。适合于刚接触阵列天线仿真、手里有一份代码但不知道怎么调参数的人,也适合想把方向图综合从“玄学”变成“可调参工程”的工程师。
2. 阵列天线方向图仿真的底层逻辑:先搞懂方向图是怎么算出来的
2.1 方向图乘积定理:阵元方向图和阵列因子是两回事
很多人拿到一个阵列天线仿真代码,第一步就去看复杂的加权系数和优化算法,结果方向图算歪了还找不到原因。实际上,阵列天线远场方向图的根基是方向图乘积定理:总方向图 = 阵元方向图 × 阵列因子。这个公式看起来简单,但它是整个阵列天线仿真里最重要的一个概念。
阵元方向图指的是单个天线单元(比如微带贴片、偶极子)自身的方向图,它只取决于阵元结构,和阵列怎么排没关系。阵列因子(Array Factor, AF)则只取决于阵元的空间位置、激励幅度和相位,和阵元长什么样没关系。二者相乘才得到完整的阵列方向图。
我平时在 matlab 里做阵列方向图仿真,第一步永远是先把阵列因子单独算一遍、画出来,确认阵列因子没问题,再乘上阵元方向图。因为阵列因子是后续所有波束成形和优化算法的作用对象,它的计算一旦出错,后面所有结果都是错的。
% 均匀线阵的阵列因子计算 N = 16; % 阵元数量 d = 0.5; % 阵元间距,单位:波长(lambda) theta = -90:0.1:90; % 俯仰角扫描范围,单位:度 theta_rad = deg2rad(theta); % 导向矢量(steering vector)——核心变量 % 对均匀线阵,第n个阵元的空间相位差是 n*d*sin(theta) AF = zeros(size(theta_rad)); for n = 0:N-1 AF = AF + exp(1j * 2 * pi * d * n * sin(theta_rad)); end AF_mag = abs(AF) / N; % 归一化幅度 % 转成dB并绘图 AF_dB = 20 * log10(AF_mag + eps); plot(theta, AF_dB, 'LineWidth', 1.5); xlabel('角度 (deg)'); ylabel('归一化幅度 (dB)'); title('均匀线阵阵列因子 (N=16, d=0.5\lambda)'); grid on;这段代码逻辑上最重要的变量是导向矢量exp(1j * 2 * pi * d * n * sin(theta_rad)),它表示第 n 个阵元相对于第 0 个阵元在 θ 方向上的相位延迟。d以波长为单位,sin(theta_rad)是空间角到相位差的转换关系。用 for 循环累加而不是直接用 matlab 的矩阵运算,是为了让每一步都看得见,出了错能定位。AF_mag = abs(AF) / N做了归一化,因为 N 个同相相加的最大幅度就是 N,除完以后主瓣峰值就是 1,方便后面转 dB。
如果你手里是一份别人给的 gbo.rar 这类代码包,我建议你进代码后的第一件事就是把它的阵列因子部分单独拉出来,用上面的逻辑重写一遍,对比结果是否一致。这是最快确认代码可靠性的方式,也是后面调波束指向和旁瓣电平的基础。
2.2 matlab 仿真方案的选型:相控阵工具箱还是手写脚本
做阵列天线仿真,matlab 里有两条路:一是直接调用 Phased Array System Toolbox 里的phased.ULA、pattern这些封装好的函数,点几下就能出方向图;二是像我上面那样手写导向矢量、手算阵列因子。两条路我都走过,给一个真实的选型判断。
工具箱路线适合验证思路、快速出图、以及做雷达系统级仿真时和波束成形算法联动。比如你要仿真一个相控阵雷达的扫描过程,用工具箱的phased.ArrayResponse配合phased.RadarTarget会很顺手,代码量少,不容易出低级错误。但工具箱的缺点是黑匣子——你很难从函数内部看到方向图是如何一步步算出来的,一旦结果异常,排查手段有限。
手写脚本路线适合做阵列综合算法研究、学习方向图原理、以及调参精细化控制。GBO 这类优化算法通常需要把目标函数(比如峰值旁瓣电平、主瓣宽度)显式写出来,工具箱的封装对象反而不方便做这种自定义目标函数。我在做方向图综合时基本是手写为主,工具箱只用来做交叉验证。
两条路还有一个重要差异是速度。相控阵工具箱的pattern函数在大规模阵列(比如几百个阵元)下经过内部优化,计算速度比新手手写的 for 循环快得多。但如果只是仿真 8 到 32 元的阵列,手写脚本的性能差距可以忽略,而且代码更容易理解、改起来更直接。
3. 用 matlab 跑通阵列天线方向图:阵元布局到方向图输出的完整实现
3.1 阵元坐标生成:线阵、平面阵和稀疏阵列的建模方式
一份阵列天线仿真代码里,阵元坐标是整个仿真空间的骨架。多数 gbo.rar 这类代码包默认给的是均匀线阵(ULA),因为它的阵列因子有解析表达式、实现简单、适合作为优化算法的起点。但实际工程里,平面阵和稀疏阵列的建模需求更常见。
我的做法是先把阵元坐标生成单独写成一段脚本,方便在不同阵列构型之间切换。下面这段代码生成了矩形平面阵的坐标,并计算了对应的导向矢量。
% 矩形平面阵阵元坐标生成 Nx = 8; % x方向阵元数 Ny = 8; % y方向阵元数 dx = 0.5; % x方向阵元间距(波长单位) dy = 0.5; % y方向阵元间距(波长单位) % 生成所有阵元坐标 [x, y] = meshgrid((0:Nx-1)*dx, (0:Ny-1)*dy); x = x(:)'; % 拉直成行向量 y = y(:)'; % 拉直成行向量 z = zeros(size(x)); % 平面阵在z=0平面 % 方位角phi和俯仰角theta的扫描网格 phi = 0:1:360; theta = 0:1:90; [PHI, THETA] = meshgrid(deg2rad(phi), deg2rad(theta)); % 平面阵的导向矢量:exp(j * k * (x*sin(theta)*cos(phi) + y*sin(theta)*sin(phi))) % k = 2*pi,因为坐标以波长为单位时波数归一化为2*pi AF = zeros(size(PHI)); for n = 1:length(x) AF = AF + exp(1j * 2 * pi * (x(n) * sin(THETA) .* cos(PHI) + ... y(n) * sin(THETA) .* sin(PHI))); end平面阵的导向矢量计算比线阵复杂的地方在于相位差是二维的:x 方向贡献x*sin(θ)*cos(φ),y 方向贡献y*sin(θ)*sin(φ)。注意这里的波数k = 2π/λ,因为坐标已经用波长归一化(dx = 0.5表示半波长),所以 k 就简化为2π。这个归一化约定很常见,但也是后面最容易出错的地方——如果坐标单位是米,k 就要写成2*pi*f/c。
生成坐标时用meshgrid然后(:)'拉直,是为了让坐标向量和后面循环里的阵元索引一一对应。用deg2rad把所有角度统一转成弧度,避免 sin/cos 函数里混用单位出问题。如果你要做的是稀疏阵列(比如标题相关热词里提到的“稀疏阵列天线雷达”),坐标就不是网格生成了,而是用一个 0/1 布尔向量表示阵元开或关,或者用优化算法直接优化阵元位置。GBO 优化算法的一个重要应用场景就是稀疏阵列的阵元位置布局优化。
3.2 方向图的可视化输出:极坐标图、直角坐标图和 3D 方向图
方向图算出来之后,怎么画直接影响你对结果的判断。我用三个层次来画方向图,每一步应对不同的排查需求。
第一层是直角坐标图,也就是最常见的角度-幅度曲线,适合看主瓣宽度、旁瓣电平、零陷位置这些具体的数值指标。第二层是极坐标图,适合看方向图的整体覆盖形态,尤其是波束指向的直观判断——极坐标图里主瓣在哪一个角度,一眼就能看出来。第三层是 3D 方向图(平面阵时用 surf 或 mesh 画),适合看整个球面空间的辐射覆盖情况,做波束扫描仿真时比较直观。
% 三种方向图绘制方式 % 假设已经算好了AF_dB:1xN的向量,对应角度theta_deg % 1. 直角坐标图 figure(1); plot(phi, AF_dB, 'b-', 'LineWidth', 1.5); xlabel('方位角 (deg)'); ylabel('幅度 (dB)'); title('直角坐标方向图'); grid on; ylim([-50 5]); % 2. 极坐标图 figure(2); polarplot(deg2rad(phi), max(AF_dB, -40)); % 限制最小显示值 title('极坐标方向图'); % 3. 3D方向图(平面阵用) figure(3); [X, Y] = meshgrid(deg2rad(phi), deg2rad(theta)); surf(X.*cos(Y), X.*sin(Y), AF_db_2d, 'EdgeColor', 'none'); xlabel('x'); ylabel('y'); zlabel('dB'); title('3D方向图');这段代码里有几个细节是我踩过坑之后才加上的。polarplot之前一定要用max(AF_dB, -40)把极小值截断,否则副瓣区域可能出现大量毛刺,极坐标图看起来像刺猬一样,信息全被噪声淹没。ylim([-50 5])同理,方向图的动态范围太大时,如果不限制 y 轴范围,主瓣和远区副瓣之间的差异会让副瓣细节完全看不见。
surf画 3D 方向图时,X.*cos(Y)和X.*sin(Y)是把球坐标投影到二维平面上的常见做法,适合做可视化。如果你需要精确读取某个方向上的增益值,还是得回到数据本身,不要直接在 3D 图上估读。
4. 天线波束成形与方向图综合:从均匀加权到 GBO 优化
4.1 波束指向控制:相位加权向量是怎么算出来的
均匀激励的阵列天线,主瓣永远指在法向(也就是线阵的 0 度方向)。要让波束偏转到某个角度,需要给每个阵元补偿一个额外的相位,让所有阵元在目标方向上的辐射同相叠加。这个相位补偿就是波束指向控制的核心。
以均匀线阵为例,假设阵元间距为 d,波束指向 θ0,则第 n 个阵元的相位补偿是-2π * d * n * sin(θ0)/λ。注意这里是负号,因为导向矢量里是正号,而加权向量要取共轭才能把主瓣搬过去。这是我见过出错率最高的地方——很多人直接拿导向矢量的共轭来当加权,方向却算反了。
% 波束指向控制的相位加权计算 N = 16; d = 0.5; % 半波长间距 theta0 = 30; % 目标波束指向,单位:度 % 计算每个阵元的相位补偿(单位:弧度) n = 0:N-1; phase_shift = -2 * pi * d * n * sin(deg2rad(theta0)); % 生成复数加权向量 w = exp(1j * phase_shift); % 用加权向量计算偏转后的阵列因子 theta = -90:0.1:90; AF_steered = zeros(size(theta)); for t = 1:length(theta) s = exp(1j * 2 * pi * d * n' * sin(deg2rad(theta(t)))); AF_steered(t) = w * s; % 加权向量乘导向矢量 end AF_steered_dB = 20 * log10(abs(AF_steered) / N); plot(theta, AF_steered_dB);代码里w * s是向量内积,等价于把所有阵元的加权限信号求和。phase_shift里用负号,不是数学上的必然,而是我约定“导向矢量取正相位累加、加权取负相位补偿”的惯例。如果你的代码库里约定相反,那波束指向结果会是镜像对称的——这在排查问题时要警惕,不是算法错了,是符号约定变了。
4.2 低旁瓣加权:切比雪夫加权和 Taylor 加权的实现对比
雷达和通信系统里,均匀加权虽然主瓣最窄、增益最高,但第一旁瓣电平高达 -13.26 dB,这个值在很多场景下不满足指标要求。压低旁瓣的经典手段是幅度加权,其中最常用的是 Dolph-Chebyshev 加权和 Taylor 加权。matlab 甚至自带了chebwin函数,可以直接生成切比雪夫加权系数。
% 切比雪夫加权与Taylor加权对比 N = 24; SLL = -30; % 目标旁瓣电平,单位:dB % 方法1:matlab自带chebwin w_cheb = chebwin(N, abs(SLL))'; % 转成行向量 % 方法2:手写Taylor加权(nbar = 4通常够用) nbar = 4; % 泰勒级数截断项数,控制主瓣宽度 w_taylor = taylorwin(N, nbar, abs(SLL))'; % 验证两种加权的方向图 theta = -90:0.1:90; AF_cheb = zeros(size(theta)); AF_taylor = zeros(size(theta)); d = 0.5; for t = 1:length(theta) s = exp(1j * 2 * pi * d * (0:N-1)' * sin(deg2rad(theta(t)))); AF_cheb(t) = w_cheb * s; AF_taylor(t) = w_taylor * s; end AF_cheb_dB = 20 * log10(abs(AF_cheb) / sum(w_cheb)); AF_taylor_dB = 20 * log10(abs(AF_taylor) / sum(w_taylor)); plot(theta, AF_cheb_dB, 'b-', 'LineWidth', 1.5); hold on; plot(theta, AF_taylor_dB, 'r--', 'LineWidth', 1.5); legend('Dolph-Chebyshev', 'Taylor'); ylim([-60 5]); grid on;两种加权的核心思想都是“用主瓣变宽换旁瓣降低”,但实现上有区别。chebwin生成的权重,理论上能把所有旁瓣压缩到指定电平以下,但代价是两端阵元的幅度权重会很小,激励幅度动态范围很大。Taylor 加权给了 nbar 这个参数来控制远旁瓣的衰减速度,实际工程里出现更多,因为它的激励幅度分布更平坦、更容易在馈电网络里实现。
参数说明:SLL = -30表示目标旁瓣 -30 dB,这是雷达系统里比较常见的要求;nbar = 4控制泰勒级数展开的项数,增大 nbar 会让方向图更接近切比雪夫,但归一化后的幅度动态范围也变大。对比图中如果发现切比雪夫和 Taylor 方向图高度重合,是正常的,两者区别主要在远旁瓣区域和激励幅度分布上,近主瓣区域差异不大。
4.3 用 GBO 做阵列方向图综合:目标函数设计与优化循环
切比雪夫和 Taylor 加权是解析方法,只能处理规则阵列。当你面对的是一个稀疏阵列、或者要求同时控制多个零陷方向时,解析方法就失效了,这时候要用优化算法。GBO 是近年在阵列天线综合里用得比较多的算法之一,它在梯度信息的基础上结合了种群迭代,对方向图综合这类非凸优化问题有不错的收敛效果。
GBO 优化方向图综合的核心是设计目标函数(适应度函数)。我最常用的一版目标函数是:主瓣指向误差的惩罚 + 旁瓣峰值电平 + 零陷深度约束。下面给出的是 GBO 方向图综合中最关键的目标函数代码框架。
function fitness = gbo_antenna_fitness(w, params) % 方向图综合适应度函数 % w: 优化变量。如果w是复数,表示幅度和相位联合优化; % 如果w是实数,前一半是幅度、后一半是相位,调用时要自行拆分 % params: 包含阵列几何、目标指向、旁瓣区域等参数的struct theta = params.theta; % 角度扫描向量(弧度) AF_pattern = zeros(size(theta)); for t = 1:length(theta) s = exp(1j * 2 * pi * params.d * params.n' * sin(theta(t))); AF_pattern(t) = w * s; end AF_db = 20 * log10(abs(AF_pattern) / max(abs(AF_pattern)) + eps); % 寻找主瓣区域:以目标指向为中心,第一零点的位置作为主瓣边界 % 简化处理:找到离目标指向最近的角度,取其两侧第一个零点作为主瓣范围 [~, mainlobe_idx] = min(abs(theta - params.theta0)); mainlobe_half_width = params.mainlobe_half_width; % 主瓣半宽,需根据阵列尺寸估算 mainlobe_mask = abs(theta - params.theta0) < mainlobe_half_width; % 旁瓣区域方向图的峰值 sll_region = AF_db(~mainlobe_mask); peak_sll = max(sll_region); % 目标指向误差惩罚 mainlobe_peak_db = max(AF_db(mainlobe_mask)); [~, actual_peak_idx] = max(AF_db(mainlobe_mask)); actual_angle = theta(mainlobe_mask); angle_error = abs(actual_angle(actual_peak_idx) - params.theta0); % 加权求和:旁瓣峰值电平(越小越好)+ 指向误差惩罚 fitness = peak_sll + params.lambda_angle * angle_error; end这段目标函数至少要理解三点。第一,mainlobe_half_width是一个需要根据阵列口径预先估算的参数,比如 16 元半波长间距线阵的主瓣半宽大约在 6 度左右,估算不准会导致旁瓣区域混入主瓣、优化方向错误。第二,+ eps是防止 log10 里出现 0 导致结果为负无穷,方向图计算中这个细节必须注意,否则优化过程中 fitness 可能出现 NaN。第三,params.lambda_angle是角度误差惩罚权重,取值一般 1 到 10,它决定了优化算法更偏向压低旁瓣还是更偏向指向精确。GBO 的种群大小通常设 30 到 50,迭代次数 100 到 300,这几个参数写在 GBO 主循环里,直接影响优化结果的收敛速度和稳定性。
5. 阵列天线仿真避坑指南:五条血泪经验
5.1 栅瓣问题:阵元间距超过半波长,方向图出现“伪主瓣”
现象:方向图在目标指向之外的角度出现了一个和主瓣幅度几乎一样的波峰,而且波峰位置随频率变化。原因:阵元间距 d 超过 λ/2 时,阵列因子的空间采样不满足奈奎斯特条件,出现空间混叠,也就是栅瓣。这是阵列天线仿真里最常见的翻车现场,而且用优化算法也修不掉——因为它在物理层面就发生了。
解决:把阵元间距改到 0.5λ 以下,或者明确你的扫描范围之后,按d ≤ λ / (1 + |sinθ_max|)来约束间距。比如扫描范围 ±60 度,d 要小于 0.536λ,只取 0.5λ 是安全的通用值。仿真时先用这个解析条件验证一下你的阵列参数,能过滤掉一大批方向图异常。
5.2 坐标约定混乱:极坐标图画出来主瓣位置不对
现象:用polarplot画方向图,主瓣应该指向 30 度,画出来却发现指向了 -30 度或 60 度。原因:matlab 的polarplot中 0 度在水平向右方向,角度逆时针增加,而天线方向图通常用的是“法向为 0 度、顺时针或逆时针计角”的球坐标系约定,两者如果没对齐就会出现镜像或旋转偏差。
解决:画图前做一次角度映射,把 θ 转换成polarplot需要的角。具体做法是:先画出不加任何相位加权的法向方向图,确认主瓣落在 90 度(polarplot 的垂直方向),再把你的角度坐标整体偏移修正。如果发现主瓣跑到 270 度了,说明角度方向反了,把 θ 取负或加 180 度即可。
5.3 单位混淆:dB 和线性幅度混用,副瓣电平算错
现象:用优化算法算出来的旁瓣电平显示为 -60 dB,实际转为线性幅度后发现旁瓣比主瓣只低 1 个数量级,完全对不上。原因:目标函数里计算旁瓣峰值时用了线性幅度,却在 fitness 值里直接和 dB 阈值的数值作比较,单位不统一。
解决:目标函数里全部使用 dB 值,或者在全部使用线性值时把阈值 0.1 当成 -20 dB 对应的线性幅度。我的惯例是在适应度函数入口处统一转成 dB,所有比较、惩罚、加权都用 dB 数值,最后输出报告也全部用 dB。这样从fitness到结果图的单位链条是一致的,排查问题才不必反复换算。
5.4 优化算法结果不稳定:每次跑出来的旁瓣电平不一样
现象:GBO 优化算法每次运行得到的方向图旁瓣电平不同,有时差 3 到 5 dB,结果不可复现。原因:GBO 和其他种群优化算法一样,初始化种群是随机的,如果不固定随机种子,每次迭代路径完全不同。另一个原因是迭代次数不够,算法还没收敛就被提前终止了。
解决:在 GBO 主循环前加一句rng(2024)(任意固定整数),保证随机种子固定。同时检查收敛曲线——如果适应度值在最后几十次迭代里还在明显下降,说明迭代次数设少了,加大到 300 或 500 次。如果固定种子后结果还是跳,那就是目标函数里存在数值不稳定,重点检查+ eps是否加对了位置。
5.5 大阵列仿真速度慢:for 循环遍历角度矩阵卡死
现象:200 元以上的阵列、0.1 度间隔扫描角度,for 循环算一次方向图要几十秒,加上优化算法循环几百次,整个仿真要跑一晚上。原因:嵌套 for 循环没有利用 matlab 的矩阵运算能力。这不是功能问题,是代码性能问题。
解决:用矩阵运算重写阵列因子计算,把回答速度提升两个数量级。核心是把阵元索引向量n和角度向量theta组织成一个二维矩阵,一次算出所有阵元在扫描角度上的导向矢量。这在 2.1 和 3.1 的代码里已经体现了雏形——用n' * sin(theta)生成二维矩阵,而不是逐角度循环。对 32 元阵列来说两种写法时间差异不大,但到了 256 元,矩阵运算的优势就是决定性的了。
6. 进阶验证技巧:拿解析公式和第三方电磁仿真给结果兜底
无论 matlab 里跑出来的方向图多漂亮,最后一步都应该做验证。我最常用的验证手段有两个。
第一个是用解析公式交叉验证。均匀线阵在等幅激励下的阵列因子有闭式解:|AF(θ)| = |sin(Nπd sinθ/λ) / sin(πd sinθ/λ)|。用这个公式把方向图画出来,和手写脚本算出来的结果对比,如果两条曲线重合,那说明导向矢量、坐标定义、归一化逻辑都没问题。我建议把这段验证代码写成一个独立脚本,每次调整阵列参数后先跑一遍,确认基础方向图没算错再继续做复杂优化。
第二个是用第三方电磁仿真软件做抽样验证。在 matlab 里做阵列方向图综合,本质上是基于理想点源模型——每个阵元都是全向辐射、阵元间没有互耦。但真实天线的阵元方向图不是全向的,阵元之间有耦合、边缘效应等额外影响。所以在方案定型前,取两三个典型工作状态(比如法向波束、30 度偏转、低旁瓣加权)到 CST 或 HFSS 里建模对照,比对主瓣指向和旁瓣包络。
这个对照过程里,一个很现实的差异是:matlab 结果的主瓣宽度和旁瓣电平和 CST 仿真结果会有偏差,通常来自互耦导致的阵元阻抗失配。如果指标要求不高,偏差在 1 到 2 dB 内可接受;如果要求严格,就需要用实测方向图数据来校准 matlab 模型。我个人习惯是把校准系数记录在脚本注释里,方便后续追溯。
另一个容易被忽略的验证是栅瓣判定。用d ≤ λ/(1+|sinθ_max|)这个条件在脚本里写一个断言,当阵列参数不满足时直接报警。这算是我给自己留的一个“后悔药”——与其在方向图里找诡异的大瓣,不如一开始就让脚本告诉你参数不合理。
做阵列仿真这几年,最大的教训就是:方向图算出来的那一刻永远不要急着下结论。先和解析解对比,再上电磁仿真验证,最后才是调参数做优化。这个过程能帮你筛掉 80% 的低级错误。GBO 这类优化算法确实能压低旁瓣,但它不是魔法,好代码 + 好验证流程才能让仿真结果真正可信。希望帮到你。
本文还有配套的精品资源,点击获取