简介面向阵列天线设计与电磁仿真方向的研究生、工程师及通信专业学习者这份Matlab天线仿真工具包聚焦阵列天线的波束形成与方向图分析适合天线课程作业、毕业设计及工程预研中的方案验证。压缩包内含1个m脚本文件整体仅620B轻量但实用脚本可配合Phased Array System Toolbox等工具箱调用完成阵列模型参数设置、单元相位/幅度计算以及方向图绘制适合用作快速验证与教学演示的起步工具。已有279人浏览学习说明其在相关课题中被持续关注。通过该脚本读者不仅能直观理解阵元数量、间距与激励信号对主瓣宽度、副瓣电平及波束指向的影响还可自行调整参数重跑脚本观察不同配置下的辐射特性变化并在此基础上扩展出更复杂的波束控制算法为雷达、无线通信等场景下的阵列天线系统仿真与性能优化打下可复用的代码基础。1. 天线方向图仿真一份能直接跑出图像的 MATLAB 脚本集天线方向图仿真听着像电磁场硬核算法真跑起来以后才发现大半时间都耗在“图不对”上直角坐标和极坐标混用、dB 归一化出错、扫描角范围漏了背面出来的结果既不像课本也不像实测。这份 gbo.rar 就是一套 matlab天线 方向图与阵列天线仿真脚本集我把里面最常用的几条路径捋了一遍从单元方向图、阵列因子到切比雪夫加权和波束扫描全部能直接在 MATLAB 里跑通出图。适合做电磁场与微波课程设计、天线方向图大作业以及想快速预研阵列方案但没时间从零推导的人。别指望它替你完成理论推导但它能把“出图、调参、验证”这条链路压缩到半天以内。2. 把方向图拆成两项单元因子与阵列因子各管一段2.1 为什么方向图仿真必须先拆成两项而不是直接堆相位很多人拿到阵列天线仿真第一反应是写一坨 for 循环把每个阵元的辐射场叠加起来最后画出来的图连自己都不知道哪儿对哪儿。工程上更稳的做法是把方向图拆成“单元方向图”和“阵列因子”两项相乘即总方向图 F(θ) 单元方向图(θ) × 阵列因子(θ)。单元方向图描述单个天线本身对空间的响应阵列因子描述阵列排布和相位加权对空间的响应两者独立。拆开写的好处是你想换阵元类型只需要换单元项想改阵元间距、加权系数、扫描角只需要改阵列因子公式和代码互不污染。我拆这类资源包时习惯性先找有没有把这两项分开写的脚本。如果没有我会自己建一套文件名和职责如下表所示脚本文件职责核心输入参数unit_pattern.m计算单元方向图默认半波偶极子可近似全向观察角 theta、极化方式array_factor.m计算阵因子均匀或加权激励N、d、theta0、加权向量 wtotal_pattern.m组合单元项与阵列因子并归一化出图上述全部scan_beam.m批量扫描波束指向并叠加显示theta0 扫描范围、步进extract_metrics.m提取主瓣宽度、第一副瓣电平方向图数据、阈值2.2 用一份可改参数的脚本画出第一个均匀直线阵方向图下面这段是均匀直线阵最基础的版本也是我每次新开阵列项目时最先跑的“冒烟脚本”确认环境没问题以后再往里面加功能。% 均匀直线阵方向图仿真边射阵 clear; clc; N 8; % 阵元数 d 0.5; % 阵元间距单位波长 theta0 0; % 主波束指向单位度0度即边射阵 theta -90:0.5:90; % 观察角度范围单位度 AF zeros(size(theta)); for n 1:N AF AF exp(1j * 2 * pi * d * (n - 1) * (sind(theta) - sind(theta0))); end AF AF / N; % 峰值归一化方便后面转 dB pattern_dB 20 * log10(abs(AF) eps); % 加 eps 防止 log10(0) pattern_dB pattern_dB - max(pattern_dB); % 归一到 0 dB figure; plot(theta, pattern_dB, LineWidth, 1.5); grid on; xlabel(观察角度 (°)); ylabel(归一化方向图 (dB)); title([均匀直线阵 N , num2str(N), , d , num2str(d), \lambda]); ylim([-40 5]);这段代码的核心逻辑是每个阵元在远场的贡献是一个相位项 exp(j·2π·d·(n-1)·(sinθ - sinθ0))其中 theta0 控制波束指向d 是阵元间距。循环里把 N 个阵元的贡献累加再取模转 dB。代码里加了 eps 是为了防止 abs(AF) 恰好为 0 时 log10 出现负无穷这个细节能避免后面一大堆排查时间。参数方面N 影响主瓣宽度N 越大主瓣越窄d 在 0.5 波长时是最常用的折中想扫大角度再看第 4 章怎么改。theta 步进取 0.5° 足够画图如果要精确读-3dB 宽度建议改成 0.1°。2.3 从方向图里提取三个工程指标主瓣宽度、第一副瓣电平、零深画出一张漂亮的图只算完成一半真正落地的指标提取才是能写进验收报告的东西。我一般在包里放一个提取脚本自动算出主瓣半功率宽度和第一副瓣电平。半功率宽度的做法是先找主瓣峰值再往左右两侧找第一个降到峰值-3dB 的位置两位置之差就是波束宽度。% 从方向图数据中提取半功率波束宽度 idx_peak find(pattern_dB max(pattern_dB), 1); left idx_peak; while left 1 pattern_dB(left) max(pattern_dB) - 3 left left - 1; end right idx_peak; while right length(pattern_dB) pattern_dB(right) max(pattern_dB) - 3 right right 1; end HPBW theta(right) - theta(left); fprintf(半功率波束宽度 %.2f°\n, HPBW);这段脚本思路很直白先定位主瓣峰值下标然后从峰值往左右两侧“下山”直到越过-3dB 线两侧角度差就是半功率宽度。注意这里要求 theta 是等间隔的所以下标差可以直接换算成角度差如果 theta 用了不等间隔要改用角度查找而不是下标查找。第一副瓣电平更简单主瓣右边第一个局部极大值就是第一副瓣直接读它的 dB 值就行。均匀激励的直线阵第一副瓣理论值约 -13.2dB如果算出来差太多优先检查是不是归一化时把 AF 多除了一个 N。3. 均匀阵与切比雪夫加权阵列天线仿真的两种主流套路3.1 均匀激励为什么旁瓣压不下去切比雪夫加权解决了什么均匀直线阵能给出的第一副瓣电平极限就压在那里约 -13.2dB不管你加到 16 个阵元还是 32 个阵元这个值几乎不变。原因在于旁瓣电平主要由幅度锥削形状决定而不是阵元数量。想压旁瓣就得引入幅度加权也就是给每个阵元分配不一样的激励幅度让阵列中间“强”、两端“弱”。经典做法有泰勒加权、汉明加权、切比雪夫加权。切比雪夫加权最大的优点是在给定副瓣电平设计指标时它能保证所有副瓣都压到该电平以下且主瓣宽度最窄——也就是说它在“旁瓣约束”和“主瓣宽度”之间取到了理论上的最优。这类脚本在 gbo.rar 这样的阵列天线仿真包里出现频率很高因为做天线波束设计时副瓣电平往往被系统指标卡死比如 -25dB 或 -30dB这时候均匀激励完全没法交差切比雪夫是第一步。需要提醒的是切比雪夫加权的代价是主瓣比均匀阵略微展宽且副瓣是等波纹分布看起来整条曲线都在目标电平附近“贴地飞行”不像泰勒加权那样副瓣单调递减。选择哪种取决于你的应用雷达旁瓣干扰场景喜欢切比雪夫的“硬约束”通信场景更常看到泰勒。3.2 用 chebwin 实现切比雪夫加权副瓣电平直接可控MATLAB 自带 chebwin 函数一行就能算出切比雪夫加权系数省去手写 Dolph-Chebyshev 多项式的麻烦。下面我把上一章的均匀直线阵改成加权版本方便你直接对比。% 切比雪夫加权直线阵方向图 clear; clc; N 8; % 阵元数 d 0.5; % 阵元间距单位波长 theta0 0; % 波束指向 R 30; % 目标旁瓣电平单位 dB正数代表 -30dB w chebwin(N, R); % 计算切比雪夫加权系数 theta -90:0.1:90; AF zeros(size(theta)); for n 1:N AF AF w(n) * exp(1j * 2 * pi * d * (n - 1) * (sind(theta) - sind(theta0))); end AF AF / max(abs(AF)); pattern_dB 20 * log10(abs(AF) eps); pattern_dB pattern_dB - max(pattern_dB); figure; plot(theta, pattern_dB, LineWidth, 1.5); grid on; title([切比雪夫加权 N , num2str(N), , 设计旁瓣 -, num2str(R), dB]); ylim([-50 5]);和均匀版相比唯一的结构变化是循环里把 1 换成了 w(n)向量化写法就是 AF w. * exp(...)两者等价。chebwin 的第一个参数是阵元数第二个参数 R 是旁瓣衰减值用正数传入比如 R30 表示目标旁瓣约 -30dB。实际画出来的第一旁瓣可能和 -30dB 有 1~2dB 偏差这是切比雪夫近似计算和离散采样导致的正常误差不是代码写错。要注意 chebwin 对 N 的奇偶性很敏感偶数的加权系数在两端接近最小值但不会归零奇数阵元时中心阵元权重最大。3.3 非均匀间隔怎么改代码稀疏阵列的旁瓣代价均匀直线阵只是阵列天线仿真里最基础的一类实际项目里阵元位置经常不是等间距的。比如稀疏阵列天线雷达里为了降低成本和硬件通道数会刻意抽掉一部分阵元用“伪随机”位置保住主瓣宽度但副瓣会明显抬高。这种情况下只需要把间距变量从标量改成向量。% 非均匀间隔直线阵位置向量写法 pos [0 0.45 0.9 1.55 2.1 2.8 3.4 4.0]; % 单位波长 N length(pos); w ones(1, N); % 加权向量可替换为 chebwin(N, R) theta0 0; theta -90:0.1:90; AF zeros(size(theta)); for n 1:N AF AF w(n) * exp(1j * 2 * pi * pos(n) * (sind(theta) - sind(theta0))); end核心差异在第 2 行位置向量直接替代了原来 (n-1)*d 的计算。这种方法能处理任意一维布阵包括稀疏阵、非规则阵、稀布阵。代码里我故意把位置间距写得有大有小模拟稀疏抽阵后的效果。跑出来你会发现主瓣宽度基本没变但副瓣明显比满阵高这就是稀疏阵列的典型代价。如果你做的是二维面阵只需要把一维位置扩展成二维坐标相位项变成 exp(j·2π·(x·sinθcosφ y·sinθsinφ))代码结构完全复用。4. 波束扫描仿真相移量、栅瓣边界与一次扫参脚本4.1 相移量与波束指向的关系别把符号搞反电扫描的本质是通过阵元间相位差补偿波程差让所有阵元的辐射在目标方向同相叠加。阵元间距为 d 时要让主瓣指向 theta0相邻阵元需要的相移量是 Δφ 2π·d·sin(theta0) / λ。当 d 以波长为单位时公式简化为 2π·d·sin(theta0)。注意相移量的符号约定如果第 n 个阵元相对第 1 个阵元的相位是 -n·Δφ波束朝正方向偏把负号去掉就会偏向负方向。我见过不止一个人在这里把符号搞反出来的方向图和预期左右镜像检查半天才发现相位项少了个负号白白浪费两小时。阵列天线仿真里把波束扫描脚本写成一个函数是最实用的习惯下面这份 function 可以放进包里反复调用。4.2 写一个波束扫描主函数一次画出多根方向图function pattern_dB array_pattern(N, d, theta0, theta, w) % 计算直线阵方向图并返回归一化 dB 值 % 输入 % N 阵元数 % d 阵元间距单位波长标量或长度 N 的位置向量 % theta0 主波束指向单位度 % theta 观察角度向量单位度 % w 阵元幅度加权长度 N默认全 1 if nargin 5 w ones(1, N); end AF zeros(size(theta)); for n 1:N AF AF w(n) * exp(1j * 2 * pi * d * (n - 1) * (sind(theta) - sind(theta0))); end pattern_dB 20 * log10(abs(AF) / max(abs(AF)) eps); end主脚本里调用这个函数循环扫过一组 theta0 就能叠出一张“波束扫描瀑布图”这是阵列天线仿真里最直观的展示方式。% 批量扫描主波束方向并叠加绘图 theta -90:0.2:90; scan_angles -60:15:60; % 扫描角度序列 figure; hold on; for theta0 scan_angles pat array_pattern(8, 0.5, theta0, theta); plot(theta, pat, LineWidth, 1.2); end grid on; xlabel(观察角度 (°)); ylabel(归一化方向图 (dB)); ylim([-40 5]); legend(arrayfun((x) [num2str(x) °], scan_angles, UniformOutput, false), ... Location, southwest);这段代码的意义在于把“改单一参数”升级成“批量扫描”一次性看出主瓣从 -60° 扫到 60° 的过程中旁瓣怎么变化、有没有栅瓣冒出来。注意我把 theta 步进取 0.2°比单根方向图更密因为扫描到边缘角度时曲线变化更剧烈步长太大容易漏掉栅瓣峰值。4.3 栅瓣的边界为什么 d 不能随便加大只要有周期性排布就会出现栅瓣风险相位加权也压不掉它。栅瓣出现条件可以用公式判断sin(theta_g) sin(theta0) ± m·λ/d。也就是当 d 对应到某个扫描角刚好满足“整数倍相位周期”时额外方向会出现和主瓣一样高的波束。边射阵 d0.5λ 时栅瓣正好不发生但一旦扫描角度拉大同样一组阵元间距就可能出事。下表给出了不同最大扫描角下 d 的安全上限最大扫描角°d 安全上限λ0边射1.0300.667450.586600.535这里的上限是根据 d 1 / (1 |sin(theta_max)|) 算出来的。比如你的系统要求扫描到 ±45°那 d 就不能超过 0.586λ如果沿用常见的 0.5λ 当然安全但阵元数量、通道数、成本都会跟着涨。做阵列天线仿真时栅瓣问题是最容易被忽略的有些人把 d 加到 0.8λ 想看窄波束结果方向图里冒出好几个“主瓣”还以为代码写错了实际上就是栅瓣。遇到这种情况先回头算一遍扫描角边界别急着改代码。5. 避坑排查方向图仿真最容易翻车的五个地方5.1 方向图曲线出现断崖式缺口甚至整段消失现象plot 出来的曲线在某几个角度直接“掉到地底”极坐标下更明显一整块区域没有数据。原因abs(AF) 在某些角度精确为 020*log10(0) 产生 -Inf绘图时这些点被当作无效值处理。常见于等间距均匀阵阵因子零点恰好落在观察角度网格上。解决取对数前加一个极小偏移量即 20log10(abs(AF) eps)或者用 20log10(max(abs(AF), 1e-6))。这两种做法本质都是给零值一个下限保证曲线连续。注意 eps 在 MATLAB 里约 2.2e-16对 -300dB 以下的细节有影响但对工程读数完全无感。5.2 方向图只显示半个空间背面直接消失现象极坐标方向图只画出一个“半圆”另一侧全空或者直角坐标里只有 0° 到 180° 的曲线-90° 到 0° 看不到。原因theta 取值范围写成了 0:180 而不是 -90:90。很多教材画方向图时习惯取 0 到 180°但那只是“从法线一侧看过去”的视角实际阵列的背瓣信息全丢了。另外早期 polar 函数对负半径的处理是取绝对值会把负角度数据折叠到正角度。解决一律用 -90:0.5:90 作为一维扫描角范围画极坐标图时用 polarplot(theta*pi/180, max(pattern_dB, -40))把低于 -40dB 的数据截断到 -40dB避免极坐标中心区域被杂乱小旁瓣填满。检查方向图对称性时对比主瓣两侧的第一副瓣高度是否一致不一致优先怀疑角度范围。5.3 主瓣指向和设定值对不上总是差一个角度现象设定 theta0 30结果主瓣出现在 -30° 或 60°甚至随频率变化位置漂移。原因相位项里的 sind(theta) 被写成了 sin(theta)或者 theta0 忘了转弧度又或者 d 没有以波长为单位而是直接用了物理间距。这些都属于单位制混用在便携脚本里特别容易发生因为从别处复制来的代码可能用了不同约定。解决统一用 sind 处理所有角度相位项里的 d 除以波长后再代入公式。我习惯在文件头部写一行注释所有角度单位 度所有距离单位 波长防止自己三天后再看时犯迷糊。符号问题导致的左右镜像就在相位项加个负号对比一下哪边对用哪边。5.4 切比雪夫加权的实际旁瓣电平和设计值差一大截现象设计 R30画出来第一旁瓣只有 -22dB 左右或者出现比设计值更高的旁瓣。原因两种情况最常见。一是 chebwin 的参数正负含义搞混MATLAB 里 R 是正数代表旁瓣衰减量你传负值就相当于设计了一个“负衰减”结果完全不是预期。二是加权系数没有落实到方向图计算里只算了均匀激励或者加权向量长度和阵元数对不上导致自动截断。解决先打印 w 向量确认 chebwin 返回值在 0~1 之间最大值在中心阵元再检查 array_pattern 函数里乘的是 w(n) 而不是 w。最后跑一个对照实验把切比雪夫加权和均匀激励的结果画在同一张图里确认加权版本的旁瓣整体下降、主瓣略微展宽这个趋势对了基本就没错。5.5 扫描大角度时方向图冒出“鬼影波束”现象theta0 扫到 45° 以上时方向图里出现一个和主瓣几乎等高的波束位置在 -30° 或更低角度。换不同的 N这个“鬼影”还在。原因这就是栅瓣不是代码 bug。前面讲过的边界公式 sin(theta_g) sin(theta0) ± m·λ/dd 超过安全上限后必然出现。很多人第一反应是查相位加权公式查了半天没问题最后发现是 d 设成了 0.8λ。解决按 4.3 节的表格反推当前扫描范围下 d 的安全上限改回来即可。如果系统指标不允许缩小 d就考虑非均匀布阵或者把阵列拆成子阵再错开排布用打破周期性的方式抑制栅瓣。这个坑在稀疏阵列天线雷达里面尤其容易踩因为稀疏抽阵会无意中放大周期性。6. 验证与出图两个检查习惯加一次干净的导出第一个检查习惯叫退化验证。把阵列参数退化到最简单的状态看结果是否符合物理直觉。比如把 N 设为 1方向图应该变成一条直线全向如果还有明显的波瓣起伏说明阵因子计算里混入了单元项或者归一化有问题。再把 d 设为 0N 设为 8所有阵元在同一位置等效还是一个全向阵方向图也应该是平的。这两步几秒钟就能跑完能筛掉八成的基础性错误。我每次新写方向图脚本第一件事不是画漂亮曲线而是先做 N1 和 d0 两个退化测试。第二个习惯是用理论公式核对数值。均匀直线阵边射时半功率宽度可近似为 HPBW ≈ 0.886·λ / (N·d)单位是弧度。用这个公式估算你的 N8、d0.5 的阵列HPBW ≈ 0.886 / 4 ≈ 12.7°。用第 2 章的提取脚本读一下实际仿真值如果两者差距超过 10%优先检查 theta 扫描范围和 -3dB 判据是否用对了。注意这个公式只在边射阵、d 接近 0.5λ 时误差较小扫描到 30° 以上时需要用 cos(theta0) 修正公式变为 HPBW ≈ 0.886·λ / (N·d·cos(theta0))扫描角越大主瓣越宽。验证通过以后就是出图。写论文或报告建议导出矢量格式而不是截图MATLAB 里用 exportgraphics 最省事exportgraphics(gcf, pattern_scan.pdf, ContentType, vector);PDF 矢量图插到 Word 或 LaTeX 里都不会糊方便后续标注。如果你要交 EPS就把后缀换成 .eps新版 MATLAB 也支持直接导出。注意 exportgraphics 对极坐标图的支持在个别版本有兼容问题遇到导出空白时先改用直角坐标出图或者用 print(-depsc, pattern_scan.eps) 保全图。我自己的习惯是无论多简单的仿真只要涉及方向图都强制走一遍“退化验证 理论宽度对照”再出图。这个习惯救过我一次那时候做 16 阵元扫描阵列d 设成 1.0λ扫到 40° 时冒出一个和主瓣一样高的栅瓣差点当成“多波束天线”写进报告。后来发现边界条件算错回头把 d 改回 0.55λ方向图才正常。从那以后每次仿真出图前我都强制走一遍这两步检查再导出矢量图。希望帮到你。本文还有配套的精品资源点击获取