简介：本资源是一套面向通信工程、电磁场与微波技术方向本科生及初学者的均匀面阵波束合成方向图仿真与分析MATLAB代码，聚焦相控阵天线设计基础能力培养，适用于课程设计、本科毕设选题与自主实验学习。压缩包共6个文件（5个.m函数脚本+1个说明txt），总大小仅5KB，轻量易用；其中uniform_surface_array为核心计算函数，其余主函数分别实现球坐标系、直角坐标系及UV坐标系下的三维/二维方向图绘制，并支持主瓣宽度、副瓣电平的一维切片对比分析。已有183人学习下载，所有代码均含密集中文注释，参数（阵元数MN、间距d、频率f、波束指向角）全部外置可调，绘图横纵坐标物理意义明确，便于直观理解阵列规模与几何参数对辐射特性的影响机制。
1. 从“点”到“面”：为什么均匀面阵波束合成值得深究
在信号处理、雷达、声纳以及无线通信领域，阵列天线是实现空间信号处理的核心。我们常常从最简单的均匀线阵（ULA）开始学习，理解波束形成、方向图、主瓣宽度、旁瓣电平等基本概念。然而，当问题从一维扩展到二维，从“线”变成“面”时，均匀面阵（Uniform Planar Array, UPA）带来的不仅仅是维度的增加，更是波束控制能力和空间分辨率的质变。很多朋友在掌握了线阵的MATLAB仿真后，面对面阵时可能会感到一丝迷茫：代码结构似乎更复杂了，方向图从二维曲线变成了三维曲面，参数也多了起来。
这篇内容，我们就来彻底拆解均匀面阵波束合成的MATLAB实现。我不会仅仅给你一段能“跑通”的代码，那样意义不大。我们将从最根本的阵列流形向量构建开始，一步步推导并实现波束加权、方向图计算与三维可视化，并重点解读方向图中每一个特征——主瓣的指向性与形状、栅瓣的产生条件与规避、旁瓣的分布规律——背后的物理意义和数学原理。无论你是正在完成课程作业的学生，还是需要快速上手阵列仿真的工程师，这篇超详细的指南都将带你越过“照猫画虎”的阶段，真正理解并掌握均匀面阵波束合成的核心。
2. 均匀面阵的数学模型：构建你的“天线面板”
在写代码之前，我们必须夯实理论基础。均匀面阵，顾名思义，是所有阵元在同一个平面上，并且以均匀的间距排列成矩形网格的阵列。这是最经典、最常用的面阵形式。
2.1 阵列几何与坐标定义
假设我们有一个包含M × N个阵元的均匀面阵。通常，我们定义阵元沿 x 轴方向有M个，间距为dx；沿 y 轴方向有N个，间距为dy。为了简化分析，我们通常将阵列平面置于 x-y 平面，并且假设所有阵元都是各向同性的（即其本身的方向图是一个球体）。
那么，第(m, n)个阵元的位置坐标可以表示为：p_mn = [(m-1)*dx, (n-1)*dy, 0]，其中m = 0, 1, ..., M-1，n = 0, 1, ..., N-1。这里我们将参考点（相位中心）设置在(0,0,0)，即第一个阵元的位置。
注意：索引从0开始还是从1开始，会影响相位计算的表达式。在MATLAB中，矩阵索引从1开始，但在理论公式中，从0开始更为简洁。我们会在代码中进行统一和转换。
2.2 阵列流形向量：空间响应的“指纹”
阵列流形向量a(θ, φ)是整个波束合成理论的基石。它描述了来自空间某个方向(θ, φ)的平面波，到达阵列中每个阵元时，相对于参考点的相对相位延迟。
- θ (theta)：俯仰角，从正 z 轴开始测量，范围通常是
[0, 180]度或[0, π]弧度。θ=0代表正上方（法线方向），θ=90度代表阵列平面内。 - φ (phi)：方位角，在 x-y 平面内从正 x 轴逆时针测量，范围是
[0, 360]度或[0, 2π]弧度。
对于一个来自方向(θ, φ)的单位幅度平面波，其在位置p_mn = (x_mn, y_mn, 0)处的波前，相对于原点(0,0,0)的波前，存在一个空间传播延迟。这个延迟表现为相位差。
波的方向矢量（单位矢量）为：k = [sinθ cosφ, sinθ sinφ, cosθ]。 那么，位置p_mn处的阵元相对于原点的波程差为：Δτ_mn = p_mn · k / c（c为波速），对应的相位差为：ψ_mn = 2π * (p_mn · k) / λ，其中 λ 是波长。
因此，阵列流形向量a(θ, φ)的第(m, n)个元素（对应第(m, n)个阵元）为：a_mn(θ, φ) = exp(-j * 2π/λ * (x_mn*sinθ cosφ + y_mn*sinθ sinφ))由于z_mn = 0，cosθ项消失。将x_mn = m*dx,y_mn = n*dy代入，得到：a_mn(θ, φ) = exp(-j * 2π/λ * (m*dx*sinθ cosφ + n*dy*sinθ sinφ))
这个复数向量包含了阵列对来自该方向信号的“固有”相位响应。波束合成的本质，就是通过一个加权向量w（与a维度相同）来与a进行内积，从而对特定方向的信号进行相干增强，对其他方向的信号进行抑制。w通常被设计为与期望方向(θ0, φ0)的流形向量共轭，即w = a*(θ0, φ0)，这就是最简单的延迟求和（DAS）波束形成器。
3. MATLAB代码实现：从理论到三维方向图
理论清晰后，我们开始动手实现。下面的代码将分为几个模块：参数定义、阵列流形计算、波束合成权重计算、方向图计算与归一化、三维可视化。我会对每一行关键代码进行注释。
%% 均匀面阵波束合成与方向图绘制 - 超详细注释版 clear; close all; clc; %% 1. 参数定义 fc = 3e9; % 载波频率 3GHz c = 3e8; % 光速 lambda = c / fc; % 波长 0.1m % 阵列规模 M = 8; % x方向阵元数 N = 8; % y方向阵元数 % 阵元间距 (通常设为半波长以避免栅瓣) dx = lambda / 2; dy = lambda / 2; % 期望波束指向 (角度制) theta0 = 30; % 俯仰角，0度为法线方向 phi0 = 45; % 方位角 % 转换为弧度制，便于计算 theta0_rad = deg2rad(theta0); phi0_rad = deg2rad(phi0); % 方向图计算的角度范围与分辨率 theta_scan = linspace(0, 90, 181); % 俯仰角扫描范围 0-90度，181个点 phi_scan = linspace(-180, 180, 361); % 方位角扫描范围 -180到180度，361个点 [Theta, Phi] = meshgrid(deg2rad(theta_scan), deg2rad(phi_scan)); % 生成网格角度，用于三维绘图 %% 2. 生成阵列流形向量 (针对期望方向) % 初始化权重向量 w，维度为 M*N x 1 w = zeros(M*N, 1); % 遍历所有阵元，计算其在期望方向下的相位，并赋给权重向量 % 权重 w 即 a*(θ0, φ0)，实现相位补偿（电子转向） idx = 0; for m = 0:M-1 for n = 0:N-1 idx = idx + 1; % 计算第(m,n)个阵元相对于原点的相位延迟 phase_shift = 2*pi/lambda * (m*dx*sin(theta0_rad)*cos(phi0_rad) + n*dy*sin(theta0_rad)*sin(phi0_rad)); % 波束合成权重是流形向量的共轭，以补偿该延迟，使该方向信号同相相加 w(idx) = exp(1j * phase_shift); % 注意这里是 +j，因为是共轭（补偿延迟） end end % 通常会对权重进行归一化，使得波束形成器对白噪声的增益为0dB w = w / sqrt(M*N); % 归一化权重 %% 3. 计算阵列方向图 % 方向图是阵列对空间中所有来波方向的响应幅度 % AF(theta, phi) = | w^H * a(theta, phi) |，其中 ^H 表示共轭转置 % 由于 w 已经是 a*(θ0,φ0)，所以 AF 在 (θ0,φ0) 处取得最大值。 % 初始化方向图矩阵 pattern = zeros(size(Theta)); % 遍历所有扫描角度，计算响应 % 警告：此部分使用双重循环，计算量较大。对于大规模阵列或高分辨率扫描，可考虑向量化优化。 for i = 1:size(Theta, 1) for j = 1:size(Theta, 2) theta = Theta(i, j); phi = Phi(i, j); % 计算当前扫描方向 (theta, phi) 的阵列流形向量 a a = zeros(M*N, 1); idx = 0; for m = 0:M-1 for n = 0:N-1 idx = idx + 1; phase = 2*pi/lambda * (m*dx*sin(theta)*cos(phi) + n*dy*sin(theta)*sin(phi)); a(idx) = exp(-1j * phase); % 注意这里是 -j，是信号到达的相位 end end % 计算波束形成器输出功率（或幅度） pattern(i, j) = abs(w' * a); % w' 是 w 的共轭转置 end end % 将方向图归一化，最大值设为 0 dB pattern_dB = 20 * log10(pattern / max(pattern(:))); %% 4. 三维方向图可视化 figure('Position', [100, 100, 1200, 500]); % 子图1：三维曲面方向图 subplot(1, 2, 1); % 将球坐标下的方向图转换为直角坐标用于绘图 [X, Y, Z] = sph2cart(Phi, pi/2 - Theta, pattern); % 注意：sph2cart输入为(方位角，仰角，半径) surf(X, Y, Z, pattern_dB, 'EdgeColor', 'none', 'FaceAlpha', 0.9); hold on; % 标记波束指向 [x0, y0, z0] = sph2cart(phi0_rad, pi/2 - theta0_rad, 1); plot3([0, x0], [0, y0], [0, z0], 'r-', 'LineWidth', 2); plot3(x0, y0, z0, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); hold off; axis equal; axis tight; grid on; view(135, 30); xlabel('X'); ylabel('Y'); zlabel('Z'); title(sprintf('3D波束方向图 (指向: \\theta=%d^\\circ, \\phi=%d^\\circ)', theta0, phi0)); cbar = colorbar; cbar.Label.String = '增益 (dB)'; colormap jet; % 子图2：二维等高线方向图（在方位-俯仰平面上的投影） subplot(1, 2, 2); contourf(rad2deg(Phi), rad2deg(Theta), pattern_dB', -40:5:0); % 转置以匹配维度 hold on; plot(phi0, theta0, 'wx', 'MarkerSize', 15, 'LineWidth', 2); % 标记波束指向 hold off; axis([-180, 180, 0, 90]); xlabel('方位角 \phi (度)'); ylabel('俯仰角 \theta (度)'); title('方向图等高线 (dB)'); colorbar; grid on; %% 5. 绘制特定切面的方向图 figure('Position', [100, 100, 1000, 400]); % 切面1：固定方位角 phi = phi0，观察俯仰角 theta 方向图 subplot(1, 2, 1); phi_fixed = phi0_rad; % 需要从全局方向图数据中提取这条“线”，这里我们重新计算更精确 theta_line = deg2rad(linspace(0, 90, 901)); resp_line = zeros(size(theta_line)); for i = 1:length(theta_line) a_line = zeros(M*N, 1); idx = 0; for m = 0:M-1 for n = 0:N-1 idx = idx + 1; phase = 2*pi/lambda * (m*dx*sin(theta_line(i))*cos(phi_fixed) + n*dy*sin(theta_line(i))*sin(phi_fixed)); a_line(idx) = exp(-1j * phase); end end resp_line(i) = abs(w' * a_line); end resp_line_dB = 20*log10(resp_line / max(resp_line)); plot(rad2deg(theta_line), resp_line_dB, 'b-', 'LineWidth', 1.5); xlabel('俯仰角 \theta (度)'); ylabel('归一化增益 (dB)'); title(sprintf('方位角固定切面 (\\phi = %d^\\circ)', phi0)); grid on; ylim([-50, 0]); xlim([0, 90]); % 切面2：固定俯仰角 theta = theta0，观察方位角 phi 方向图 subplot(1, 2, 2); theta_fixed = theta0_rad; phi_line = deg2rad(linspace(-180, 180, 1801)); resp_line2 = zeros(size(phi_line)); for i = 1:length(phi_line) a_line = zeros(M*N, 1); idx = 0; for m = 0:M-1 for n = 0:N-1 idx = idx + 1; phase = 2*pi/lambda * (m*dx*sin(theta_fixed)*cos(phi_line(i)) + n*dy*sin(theta_fixed)*sin(phi_line(i))); a_line(idx) = exp(-1j * phase); end end resp_line2(i) = abs(w' * a_line); end resp_line2_dB = 20*log10(resp_line2 / max(resp_line2)); plot(rad2deg(phi_line), resp_line2_dB, 'r-', 'LineWidth', 1.5); xlabel('方位角 \phi (度)'); ylabel('归一化增益 (dB)'); title(sprintf('俯仰角固定切面 (\\theta = %d^\\circ)', theta0)); grid on; ylim([-50, 0]); xlim([-180, 180]);运行这段代码，你将得到三张图：一张三维方向图曲面，一张二维等高线图，以及两张分别固定方位角和俯仰角的切面方向图。现在，波束应该清晰地指向你设定的(30°, 45°)方向。
4. 方向图核心特征深度解读：主瓣、栅瓣与旁瓣
生成方向图只是第一步，更重要的是读懂它。均匀面阵的方向图具有一些标志性的特征，理解它们对于阵列设计至关重要。
4.1 主瓣形状与波束宽度
在三维方向图上，主瓣不再是一个简单的扇形，而是一个指向特定方向的“铅笔状”波束。其横截面（在等增益面上）近似一个椭圆。
半功率波束宽度（HPBW）：在主瓣峰值功率下降3dB（幅度下降至约0.707倍）处，波束的角宽度。对于均匀面阵，当波束指向法线方向
(θ0=0°)时，两个主平面（φ=0°和φ=90°）的HPBW可以近似估算：- 方位面HPBW (φ面)：
HPBW_φ ≈ 0.886 * λ / (M * dx * cosθ0)（弧度） - 俯仰面HPBW (θ面)：
HPBW_θ ≈ 0.886 * λ / (N * dy * cosθ0)（弧度） 注意，当波束扫描离轴时，有效孔径投影会减小，导致波束展宽，公式中体现为cosθ0项。这就是为什么我们代码中指向30°时，主瓣比指向0°时要宽一些。你可以尝试修改theta0为0和60，对比切面图中的主瓣宽度变化。
- 方位面HPBW (φ面)：
主瓣指向精度：我们的权重
w是理想相位补偿，因此主瓣峰值能精确指向(θ0, φ0)。在实际系统中，阵元通道间的幅度/相位误差、阵元位置误差都会导致指向偏差和主瓣畸变。
4.2 栅瓣：必须警惕的“幽灵”波束
栅瓣是周期性阵列（如均匀阵列）特有的现象。当阵元间距大于半波长时，在某些空间方向上，相邻阵元接收信号的相位差可能恰好是2π的整数倍。这会导致这些方向上的信号也能同相相加，形成与主瓣增益相同的波束，即栅瓣。
栅瓣出现条件（以x方向为例）：dx/λ * (sinθ cosφ - sinθ0 cosφ0) = ±k, 其中 k 为整数。 当波束指向法线方向(0,0)时，条件简化为：dx/λ * sinθ cosφ = ±k。 为了避免在可见空间（-90°≤θ≤90°）内出现栅瓣，通常要求dx < λ且dy < λ。最严格且常用的选择是dx = dy = λ/2，这能保证在任意波束指向上，可见空间内只有一个主瓣，没有栅瓣。
实操心得：在代码中，你可以尝试将
dx和dy改为0.7*lambda或lambda，然后重新运行。观察三维方向图，特别是在与主瓣对称的位置，是否出现了另一个增益很高的波瓣？那就是栅瓣。栅瓣会带来严重的空间模糊，在雷达中导致虚警，在通信中导致干扰，因此阵列设计时必须避免。
4.3 旁瓣结构：能量的“泄漏”
旁瓣是主瓣之外的其他局部增益峰值。均匀加权（即我们代码中使用的等幅加权）的阵列，其第一旁瓣电平是固定的，约为 -13.2 dB。这是因为均匀分布的孔径，其傅里叶变换（方向图）是 sinc 函数形状。
- 旁瓣分布：在三维方向图中，旁瓣呈现复杂的结构，不是简单的环状。在等高线图上，你可以看到围绕主瓣的多层闭合等高线，它们代表了不同电平的旁瓣。
- 降低旁瓣：-13.2dB的旁瓣在很多应用中仍然太高。为了降低旁瓣，需要对阵列权重
w进行幅度加权（也称为“窗函数”或“锥削”）。常见的窗函数有汉明窗、汉宁窗、切比雪夫窗等。对权重向量w的每个元素乘上一个幅度衰减系数（中心阵元增益高，边缘阵元增益低），可以显著压低旁瓣，但代价是主瓣会略微展宽，并且阵列增益（方向性系数）会有所下降。
%% 示例：应用汉明窗进行幅度加权以降低旁瓣 % 创建二维汉明窗 hamming_x = hamming(M); % Mx1 汉明窗 hamming_y = hamming(N); % Nx1 汉明窗 hamming_2d = hamming_x * hamming_y'; % 外积，得到 MxN 的二维汉明窗 % 将二维窗拉成向量，并与相位权重逐点相乘 w_tapered = w .* reshape(hamming_2d, M*N, 1); % 重新归一化 w_tapered = w_tapered / sqrt(sum(abs(w_tapered).^2)); % 使用 w_tapered 代替原来的 w 重新计算方向图，观察旁瓣电平的变化。5. 性能分析与进阶探索：超越基础DAS
我们目前实现的是最基本的延迟求和波束形成器。在实际应用中，我们关心更多指标。
5.1 方向性系数与阵列增益
- 方向性系数 (Directivity)：描述阵列将能量集中到某个方向的能力，定义为最大辐射强度与平均辐射强度之比。对于均匀加权、无损耗的均匀面阵，其最大方向性系数
D_max可近似为：D_max ≈ 4π * A / λ^2，其中A = M*dx * N*dy是阵列的物理面积。这直观地告诉我们，阵列越大（孔径越大），波束越窄，方向性越好。 - 阵列增益 (Array Gain)：在存在噪声和干扰的场景下，阵列增益比方向性更实用。它考虑了加权向量
w和信号/噪声的统计特性。对于白噪声背景下的点源信号，阵列增益等于|w^H a|^2 / (w^H w)。在我们均匀加权且归一化的情况下，阵列增益为M*N（线性值），即10*log10(M*N)dB。这就是我们常说的“8x8的阵列能提供约18dB的增益”。
5.2 扫描角度的限制与栅瓣再讨论
当波束扫描到较大角度时（例如θ0 > 60°），即使阵元间距为半波长，也可能在可见空间边缘(θ = ±90°)附近出现栅瓣。这是因为波程差公式中的sinθ项在角度大时变化剧烈。因此，对于需要宽角度扫描的相控阵，阵元间距通常需要取得更小（例如λ/3或更小），以在整个扫描空域内保证无栅瓣。这被称为“阵元间距与扫描角的折衷”。
5.3 从DAS到自适应波束形成
我们的代码是数据无关的（Data-Independent Beamforming），权重只取决于期望方向。更高级的方法是数据相关的自适应波束形成（如MVDR, LCMV），它能根据实际接收到的干扰和噪声环境，自适应地调整权重，在期望信号方向形成主瓣的同时，在干扰源方向形成零陷。其核心是求解一个约束优化问题，需要估计接收数据的协方差矩阵Rxx。
% 伪代码示意：MVDR波束形成器 % 假设已有接收数据矩阵 X (阵元数 x 快拍数) Rxx = X * X' / size(X, 2); % 估计协方差矩阵 a0 = a(theta0, phi0); % 期望方向流形向量 % MVDR最优权重： w_mvdr = Rxx^-1 * a0 / (a0^H * Rxx^-1 * a0) w_mvdr = (Rxx \ a0) / (a0' * (Rxx \ a0)); % 然后用 w_mvdr 计算方向图实现自适应波束形成需要模拟或采集包含干扰和噪声的多通道数据，这超出了本篇基础内容的范围，但它是均匀面阵应用的自然延伸。
6. 代码优化与工程实践建议
最初的演示代码为了清晰，使用了多重循环，计算效率较低。在实际工程或大规模仿真中，我们需要进行向量化优化。
6.1 向量化计算阵列流形
我们可以利用MATLAB的矩阵运算能力，避免最内层的双重循环。
%% 优化后的方向图计算（向量化版本） % 生成阵元位置索引矩阵 [m_idx, n_idx] = meshgrid(0:M-1, 0:N-1); % 维度为 N x M m_idx_vec = m_idx(:); % 拉成列向量 n_idx_vec = n_idx(:); % 拉成列向量 % 预计算常数 kx = 2*pi/lambda * dx; ky = 2*pi/lambda * dy; % 计算波束合成权重（向量化） phase_shift_vec = kx * sin(theta0_rad) * cos(phi0_rad) * m_idx_vec + ... ky * sin(theta0_rad) * sin(phi0_rad) * n_idx_vec; w_vec = exp(1j * phase_shift_vec); w_vec = w_vec / sqrt(M*N); % 计算整个方向图（向量化，效率提升巨大） % 思路：对于每个扫描角度(θ, φ)，计算所有阵元的相位，然后与权重向量做点积。 % 我们可以利用矩阵乘法一次性计算多个角度。 % 这里展示一个更高效的方案：重构方向图计算循环 pattern_opt = zeros(length(phi_scan), length(theta_scan)); for i = 1:length(theta_scan) theta = deg2rad(theta_scan(i)); for j = 1:length(phi_scan) phi = deg2rad(phi_scan(j)); % 计算当前角度的流形向量（向量化） phase_vec = kx * sin(theta) * cos(phi) * m_idx_vec + ... ky * sin(theta) * sin(phi) * n_idx_vec; a_vec = exp(-1j * phase_vec); % 到达相位 pattern_opt(j, i) = abs(w_vec' * a_vec); end end pattern_opt_dB = 20*log10(pattern_opt / max(pattern_opt(:)));向量化后，特别是如果进一步将phi循环也向量化（通过更大的矩阵运算），速度可以提升数十甚至上百倍。对于8x8这样的小阵列差别不大，但对于32x32或更大的阵列，优化是必须的。
6.2 工程中的非理想因素考量
仿真代码是理想的，但实际系统充满挑战：
- 阵元互耦：相邻阵元之间的电磁耦合会改变每个阵元的有效方向图和阻抗，导致计算的流形向量
a失真。需要在仿真中引入互耦矩阵C，此时实际流形变为a_actual = C * a_ideal。 - 通道不一致性：每个接收通道的放大器、滤波器、ADC等器件存在增益和相位偏差。这需要在权重
w中引入额外的校准系数，或通过通道校正算法来补偿。 - 宽带信号：上述分析针对窄带信号。对于宽带信号，不同频率分量在固定延迟下的相位差不同，会导致波束色散（主瓣指向随频率变化）。需要使用真时延线（TDL）或频域处理（如子带分解）来形成宽带波束。
- 有限快拍数：在自适应波束形成中，协方差矩阵
Rxx需要用有限长度的数据样本估计，估计误差会导致性能下降，特别是高维情况下（大阵列）。
理解均匀面阵的理想模型是第一步，认识到这些非理想因素并学会在仿真中建模它们，是迈向实际工程应用的关键一步。你可以尝试在代码中引入随机的小幅度/相位误差到权重向量w中，观察方向图（特别是旁瓣和零点深度）如何变得粗糙和不对称，这能让你对系统鲁棒性有更直观的认识。
本文还有配套的精品资源，点击获取
