企业官网建设流程全解析

简介：MATLAB仿真技术在半波振子天线子阵方向图研究中具有重要应用。此份zip资源包面向天线设计与电磁仿真初学者、通信工程专业学生及相关工程人员，提供完整的半波振子子阵方向图仿真实现，帮助理解振子天线辐射特性和子阵方向图生成原理。压缩包共2个文件，包含1个MATLAB脚本文件和1个仿真结果图片，整体大小仅45KB，脚本涵盖天线参数设置、子阵构建、仿真计算与方向图绘制等关键环节，图片直观展示了半波振子子阵的辐射强度分布。目前已有613人学习下载。通过运行该脚本，读者可快速获得可视化方向图，并能通过调整振子长度、馈电点位置、子阵排列等参数深入探究天线性能，为实际阵列设计与教学科研提供可复用的仿真参考，也可为赋形波束、宽角扫描等复杂辐射模式研究提供基础。

1. 半波振子子阵方向图仿真：单元方向图和阵因子到底谁说了算

平时画半波振子的方向图，读到的几乎都是那个“8”字形截面，一旦把它拼成子阵，整个方向图就会多出若干瓣和零陷，主瓣也被压窄。真正决定主瓣方向和零点位置的往往不是单元方向图，而是子阵阵因子；单根振子的单元方向图只相当于给阵因子加了一个“罩子”。banbozhenzi.zip 里的 banbozhenzi.m，做的就是把这套方向图乘积定理落成 MATLAB 代码，最后生成类似 3.2figure.jpg 的极坐标子阵方向图。它适合三类人：正在看阵列天线教材的学生、要做波束扫描的射频工程师、以及想快速对照全波仿真结果的系统设计人员。

这个压缩包里的程序不长，核心逻辑就是“单元方向图 × 阵因子”，但每一行参数设置都直接反映天线工程习惯。下面从单元建模开始，一步步拆开这套 matlab 仿真。

2. 半波振子的单元方向图建模：从解析式到 MATLAB 代码

2.1 为什么半波振子单元方向图可以直接公式化

半波振子是最典型的线天线，振子长度等于半个工作波长，电流沿振子呈近似正弦分布、末端电流趋于零。它的远场方向性有解析表达式，不需要像微带贴片那样非得靠全波求解器才能得到稳定结果。若振子沿 z 轴放置，球坐标系下的归一化电场方向图写作：

[ E_{unit}(\theta) = \frac{\cos\left(\frac{\pi}{2}\cos\theta\right)}{\sin\theta} ]

这个式子的特点是：水平面（θ = 90°）方向值最大，归一化后为 1；沿振子轴线方向（θ = 0° 和 180°）辐射为零，所以方向图是典型的“8”字形。做子阵仿真时，单元方向图通常提前算好并归一化，后续直接和阵因子相乘。

工程里常有人直接把这段公式抄进 MATLAB，结果画出来的方向图在轴向出现一个很大的跳变，原因是分母 sinθ 在 θ 接近 0 或 π 时趋近于零，分子也趋近于零，代码里直接相除会得到 NaN。这不是物理问题，而是数值处理问题，需要在计算网格上特殊处理。

2.2 参数初始化：把频率、波长、阵元数写清楚

打开 banbozhenzi.m，你会发现第一段几乎都在做参数声明。这类脚本我习惯按下面的顺序组织：

c = 2.99792458e8; % 光速，单位 m/s freq = 2.4e9; % 工作频率，单位 Hz，这里以 2.4 GHz 为例 lambda = c / freq; % 波长，约 0.125 m N = 8; % 子阵单元个数 d = lambda / 2; % 单元间距，取半波长 theta = linspace(0, pi, 721); % 极角网格，从 0 到 pi，分辨率 0.25°

参数说明：N 是参与合成的振子数量，它直接决定阵因子波束宽度；d 是相邻振子相位中心的距离，通常取半波长，取太大会在可见区出现栅瓣，这个后面第五章专门讲。theta 网格建议取 721 或 3601 点，点数太少会导致 -3 dB 波束宽度读数粗糙，点数太多则没必要，因为解析公式算起来很快，瓶颈只在绘图。

2.3 单元方向图计算的数值处理

单元方向图的计算代码可以这样写：

theta_safe = theta; % 将轴向（sin theta 接近0）置为 NaN，先避开除零 theta_safe(abs(sin(theta)) < 1e-6) = NaN; E_unit = cos(pi/2 * cos(theta_safe)) ./ sin(theta_safe); % 轴向辐射物理上为 0，这里把之前置 NaN 的位置补回 0 E_unit(abs(sin(theta)) < 1e-6) = 0; % 归一化，让最大值等于 1 E_unit = E_unit / max(abs(E_unit));

逻辑说明：第一处替换是为了让 MATLAB 在计算 NaN 时不会报警告；第二处替换是利用半波振子轴方向电场为零的物理结论，把轴向值强制设为 0。这里的阈值 1e-6 不是固定标准，如果 theta 网格点足够密，取 1e-8 也一样，关键是保证 sinθ 不为零。

有人在网上抄的代码喜欢写E_unit(isnan(E_unit)) = 1，对半波振子来说这是错的。轴向零点是半波振子的固有特性，不是数值异常，改成 1 会把后续子阵方向图的零点全部抬高，副瓣电平直接失真。

3. 方向图乘积定理：子阵方向图的合成逻辑

3.1 阵因子和单元方向图之间的边界在哪里

子阵方向图不是把每根振子的方向图简单相加，而是先算阵因子，再和单元方向图相乘。方向图乘积定理的适用条件是：单元之间互耦可忽略、阵列处于远场观察区域、各单元方向图相同。它的表达式很简洁：

[ F_{total}(\theta) = F_{unit}(\theta) \times AF(\theta) ]

AF 就是阵因子，它只跟阵列几何布局、单元间距、激励幅度和相位有关。对于沿 z 轴等间距排布的 N 元线阵，阵因子可写成：

[ AF(\theta) = \sum_{n=0}^{N-1} w_n e^{j n k d \cos\theta} ]

其中 k = 2π/λ 是波数，w_n 是第 n 个单元的复激励。这个公式的物理含义是：每个单元到远场观察点的路径差不同，产生的相位差叠加后形成干涉图样。均匀激励时 w_n = 1，阵因子有解析闭合形式，但用循环求和的方式写更直观，也方便以后加幅度锥削。

3.2 均匀线阵的法向阵因子实现

先看不扫描、无幅度加权的情况，MATLAB 代码可以这样写：

psi = 2 * pi * d / lambda * cos(theta); % 相邻单元相对观察方向的空间相位差 AF = zeros(size(theta)); for n = 0:N-1 AF = AF + exp(1j * n * psi); end AF = abs(AF) / N; % 除以 N，使最大值为 1

代码逻辑：循环变量 n 从 0 开始，对应第 0 个阵元取参考相位；psi 是一个随 theta 变化的一维数组，exp(1jnpsi) 表示第 n 个单元相对参考单元的相位偏移。最后取模并除以 N 完成归一化。

如果想写成更精简的闭合形式，可以用abs(sin(N*psi/2) ./ (N*sin(psi/2)))，但注意 psi 接近 0 和 2π 的整数倍时会出现 0/0，需要额外做数值保护。工程上我更推荐循环写法，因为后续加相位扫描、加幅度锥削时改动更小。

均匀激励的阵列，副瓣电平理论上固定在 -13.26 dB，和单元数无关，但主瓣宽度会随 N 增加而变窄。这是判断代码是否写对的第一条标准。

3.3 合成子阵方向图并转换成分贝值

把单元方向图和阵因子乘起来，就得到完整的子阵方向图：

E_total = E_unit .* AF; E_total_db = 20 * log10(abs(E_total) + eps);

后面加上 eps 是为了避免 0 dB 以下的极小值在取对数时出现 -Inf，影响后续绘图。要注意的是，阵列做波束扫描以后，单元方向图 E_unit 基本不变，变的只有 AF，所以如果看到扫描后方向图整体被压下去，那不是代码错了，而是单元方向图的“罩子”作用。

不同单元数下的均匀线阵方向图特征如下：

这个表在初步验证 matlab 仿真结果时非常有用。看到 -13.26 dB 附近的副瓣，说明阵因子代码基本正确；如果副瓣变成了 -6 dB 左右，大概率是归一化时少除了 N。

4. 复现 banbozhenzi.m 的实际流程：从工作区到 3.2figure.jpg

4.1 脚本结构按四段式组织

这一类半波振子仿真脚本，结构上通常分成四段：第一段清空工作区并声明参数，第二段计算单元方向图，第三段生成阵因子，第四段合成并绘图。banbozhenzi.m 虽然文件名简短，但跑通后输出的 3.2figure.jpg 正是第四段绘图的产物。

我重写这类脚本时，会先把骨架立出来：

clear; close all; clc; % 参数声明 freq = 2.4e9; lambda = 3e8 / freq; N = 8; d = lambda / 2; theta = linspace(0, pi, 721); % 单元方向图 theta_safe = theta; theta_safe(abs(sin(theta)) < 1e-6) = NaN; E_unit = cos(pi/2 * cos(theta_safe)) ./ sin(theta_safe); E_unit(abs(sin(theta)) < 1e-6) = 0; E_unit = E_unit / max(abs(E_unit)); % 阵因子 AF = zeros(size(theta)); for n = 0:N-1 AF = AF + exp(1j * n * 2*pi*d/lambda * cos(theta)); end AF = abs(AF) / N; % 合成方向图 E_total = E_unit .* AF; E_total_db = 20 * log10(abs(E_total) + eps);

这段代码没有任何工具箱依赖，只要有 MATLAB 基础环境就能运行。有些教材会把 theta 范围写成 -90° 到 90°，那是把坐标原点放在阵面法向；这里用 0 到 π，是因为单元方向图的解析式天然以振子轴线为基准，混用两套坐标，容易把方向图画反。

4.2 极坐标绘图与坐标轴参数

生成 3.2figure.jpg 这一步，推荐用 polarplot 而不是 plot，因为天线方向图在极坐标下能直接读出主瓣指向和零陷夹角：

figure('Name', '半波振子子阵方向图'); polarplot(theta, E_total_db, 'LineWidth', 1.5); rlim([-40 0]); % 动态范围 40 dB，低于 -40 dB 的旁瓣压到圆外 thetalim([0 180]); % 只显示 0 到 180 度，对应半空间 grid on;

rlim 参数是这里最关键的一步。方向图全动态范围可能超过 60 dB，直接绘图会让主瓣占满整个极坐标半径，副瓣细节全部被压扁。工程上一般取 30 至 40 dB 动态范围，既能看清主瓣和第一副瓣，又能避免噪声底抬起来干扰判断。

thetalim 是否要做限制，取决于你要看全空间还是只看前半空间。半波振子在轴线方向本来就存在零点，所以 0° 附近自然会凹下去，不用刻意放大。

4.3 从方向图数据里读取副瓣和波束宽度

运行脚本后，除了看图，还应该直接从数组里算指标。手动把鼠标放到图上读数，不同版本的 MATLAB 图形窗口读数习惯不同，效率太低。可以这样自动提取：

E_total_db = E_total_db - max(E_total_db); % 重新归一化到 0 dB % 主瓣位置 [~, idx_max] = max(E_total_db); theta_deg = rad2deg(theta); theta_main = theta_deg(idx_max); % 3 dB 波束宽度 idx_hp = find(E_total_db >= -3); HPBW = theta_deg(max(idx_hp)) - theta_deg(min(idx_hp)); % 第一副瓣电平：把主瓣周围 ±5 度剔除后取最大值 mainlobe_mask = abs(theta_deg - theta_main) <= 5; pattern_for_sll = E_total_db; pattern_for_sll(mainlobe_mask) = -inf; SLL = max(pattern_for_sll);

说明：HPBW 计算使用的是角度网格步长，步长越大结果越粗糙；如果 theta 用 721 点，步长 0.25°，波束宽度读数精度已经足够。SLL 的剔除宽度取 ±5° 是经验值，只适用于窄波束情况，波束扫描到 60° 以后主瓣宽度明显加大，需要按主瓣宽度的一半动态设置剔除区间。

5. 参数如何影响扫描、栅瓣和副瓣：调参实战

5.1 阵元间距 d 与栅瓣边界

把 d 从 λ/2 加大到 λ，甚至 1.5λ，可见区会出现栅瓣，也就是第二个幅度与主瓣相同的主瓣。栅瓣出现的条件是：

[ \sin\theta_{gl} = \frac{m\lambda}{d}, \quad m = \pm 1, \pm 2, \dots ]

当这个方程在 [-1, 1] 范围内有解时，栅瓣就进入方向图。下表总结了不同间距下的表现：

仿真时可以直接改 d 参数重跑脚本，看到第二主瓣出现时不用怀疑代码，那是间距过大导致的物理现象。需要注意的是，均匀线阵在 d = λ 时，栅瓣恰好出现在 θ = 90° 和 θ = 90° 对称位置，d 继续加大，栅瓣会进入前半空间并与主瓣争夺能量。

5.2 用相位梯度实现波束扫描

实际子阵设计里经常要求主瓣指向某个特定方向 θ_scan，这时需要对每个单元施加相位补偿：

theta_scan = deg2rad(30); % 期望主瓣指向 30° % 计算相邻单元之间的扫描相位差 phase_step = 2 * pi * d / lambda * cos(theta_scan); w = exp(-1j * phase_step * (0:N-1)); % 共轭相位，使主瓣向指定方向移动 AF = zeros(size(theta)); for n = 1:N AF = AF + w(n) * exp(1j * 2*pi*d/lambda * cos(theta) * (n-1)); end AF = abs(AF) / max(abs(AF));

相位为什么要取共轭？因为阵因子求和里，第 n 个单元在 θ_scan 方向本来就自带一个正相位exp(j n k d cosθ_scan)，想让它们在该方向同相叠加，就给每个单元乘一个大小相等、符号相反的相位。这样在 θ_scan 方向所有单元贡献同相相加，而其他方向相位对消。

扫描以后有两个边界条件要同时检查：一是主瓣扫描到 60° 以上时，单元方向图增益明显下降，这是阵元的“罩子效应”，不是阵列故障；二是 d 必须满足：

[ \frac{d}{\lambda} \le \frac{1}{1 + |\cos\theta_{scan}|} ]

例如扫描到 60°，cosθ = 0.5，此时 d/λ 必须小于 0.667，超过这个值就可能看到栅瓣进入可见区。这是调 d 参数时必须对照的判据。

5.3 幅度锥削对副瓣的抑制

均匀分布的副瓣电平固定 -13.26 dB，想进一步压副瓣就得做幅度锥削。最常用的是道尔夫-切比雪夫加权，它能在给定副瓣电平下让主瓣宽度最窄：

N = 8; SLL_desired = -30; % 目标副瓣电平 w_amp = chebwin(N, -SLL_desired); % 输入正值，例如 30 % 把幅度加权乘进阵因子 AF = zeros(size(theta)); for n = 1:N AF = AF + w_amp(n) * exp(1j * 2*pi*d/lambda * cos(theta) * (n-1)); end AF = abs(AF) / max(abs(AF));

chebwin 需要 Signal Processing Toolbox，如果没有该工具箱，可以用汉明窗做近似，代价是副瓣电平会均匀下降到约 -40 dB，但主瓣会展宽 1.3 倍左右。做 matlab 仿真时，副瓣和波束宽度之间的取舍关系比具体窗函数更重要，代码写对的前提下，副瓣压得越低，主瓣必然越宽，这不是数组问题。

6. 子阵方向图的快速校验与全波对拍方法

6.1 用 3 dB 波束宽度反向验证方向性

仿真结果对不对，不能只看主瓣是不是在 0°。均匀线阵的法向波束宽度有理论近似公式：

[ HPBW \approx 0.886 \frac{\lambda}{N d} ]

以 N = 8、d = λ/2 为例，HPBW ≈ 0.886 / 4 rad ≈ 12.7°。代码里算出的 HPBW 如果明显偏离这个值，优先检查 theta 网格是否包含了完整的可见区，其次检查阵因子归一化是否正确。扫描到 θ_scan 后，波束宽度按 1/cosθ_scan 展宽，如果扫描到 60° 还不变宽，说明相位加权代码写错了。

6.2 与 HFSS、CST 全波结果对拍

方向图乘积定理忽略了单元互耦，而全波仿真会包含互耦、边缘截断效应和实际馈电结构的影响。工程上最常见的做法是先跑 HFSS 得到一个单元的辐射方向图和 S 参数，再用 MATLAB 做阵因子合成。如果两者在副瓣区域差异超过 2 dB，通常是互耦导致的单元方向图畸变在起作用。修正办法不是把全波结果直接当作单元方向图，而是提取阵列中每个单元的有源方向图，再重新合成，这是阵列天线仿真中更贴近真实的一步。

6.3 把脚本结果导出成对比图

多次调参后，最好把不同扫描角度的方向图叠在一张图里保存。用 R2020a 之后的版本，可以这样导出：

figure('Name', '多波束对比'); hold on; for scan_deg = [0 15 30 45] % 此处重新计算 AF 并合成 E_total_db polarplot(theta, E_total_db, 'DisplayName', [num2str(scan_deg) '°']); end rlim([-40 0]); legend('show', 'Location', 'southoutside'); exportgraphics(gcf, 'subarray_scan_compare.png', 'Resolution', 300);

exportgraphics 比 saveas 更稳定，不会出现白色边距裁不掉的问题。这样保存出来的图可以直接贴进仿真报告，审阅人一眼就能看出波束扫描过程中主瓣宽度和副瓣电平的变化趋势。

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