MATLAB波束形成实现:线阵、面阵与圆阵的导向矢量解析
2026/9/14 12:48:42 网站建设 项目流程

简介:MATLAB波束形成程序源码包面向信号处理与阵列信号处理方向的初学者及有一定经验的开发人员,涵盖线阵、平面阵、圆阵等常见阵列结构下的波束形成实现,可直接用于算法验证、课程设计与项目参考。压缩包共2个文件,包含1个M程序源文件和1个Word文档,M文件提供带详细注释的完整代码,Word文档用于补充算法原理与使用说明,整体大小仅13KB,轻量易用。已有1368人学习浏览,具备实践参考价值。通过源码可理解常规波束形成、阵列响应计算等核心流程,借助注释和文档还能快速梳理不同阵列构型下的参数设置与实现差异,适合在MATLAB环境中边读边练,帮助掌握波束形成的基本编程方法并支持后续二次开发。

1. 波束形成不是调参,是把阵列响应写成方程

接手这个 MATLAB 波束形成源码包之前,我先说一个反直觉的结论:波束形成的核心难点不在“形成”本身,而在“导向矢量”的构造。阵列流型矩阵写对了,后面加窗、加权、扫描都是矩阵乘法和循环;写错了,仿真图再漂亮,拿到实际阵上也白搭。

这个源码包覆盖了线阵、平面阵和圆阵三种几何构型,每一段都带注释,适合两类人:一类是刚接触阵列信号处理、想把手推公式变成可运行代码的初学者;另一类是做雷达或通信系统、需要快速验证波束指向和旁瓣电平的工程师。源码里没有用到工具箱里封装好的phased.ArrayResponse,而是直接用exp(-j*k*... )构造导向矢量,这个选择很关键——它能让你看清波束形成的底层数学,而不是把 API 当黑盒。

下文从线阵的导向矢量开始,逐步扩展到平面阵和圆阵,最后再讲一个我实际项目里踩过的坑:阵元位置误差对波束指向的影响有多大,以及怎么用程序去评估它。

2. 均匀线阵的导向矢量与延迟-相加波束形成器

2.1 窄带假设下,波束形成就是一组复加权

波束形成的物理基础是干涉:同频信号从不同阵元到达时存在波程差,导致相位差。对窄带信号(带宽远小于载频),时延可以等价为相移exp(-j*2*pi*f*τ)。因此,一个 M 元阵列的输出就是各阵元接收信号乘上复权值后求和:

y(t) = w^H * x(t)

其中w是权重向量,x(t)是各阵元接收信号向量。常规波束形成(CBF,Conventional Beamforming)的权重就是导向矢量本身,即对各阵元信号做相位补偿后同相叠加,让某一方向的来波获得最大增益。

2.2 线阵导向矢量:间距、波长和入射角的关系

对于 M 元均匀线阵(ULA,Uniform Linear Array),阵元间距为 d,入射方向与阵列法线夹角为 θ,则相邻阵元间的波程差为d*sin(θ),对应相位差为2*pi*d*sin(θ)/λ。以第一个阵元为参考,第 m 个阵元的导向矢量分量是:

s_m(θ) = exp(-j*2*pi*m*d*sin(θ)/λ)

注意这里用负号还是正号,取决于你定义的入射方向正方向。源码里用的是负号,对应的物理含义是:从 θ 方向来的平面波,到达第 m 个阵元的时间比参考阵元晚,相位滞后。

把上面公式翻译成 MATLAB,就是整个源码包里复用率最高的一段:

function a = ula_steering_vector(M, d_lambda, theta_deg) % 均匀线阵导向矢量 % 输入: % M - 阵元数量 % d_lambda - 阵元间距,以波长为单位(如 d_lambda=0.5 表示半波长间距) % theta_deg - 入射方向,单位度(法线方向为 0 度,正方向为顺时针) % 输出: % a - M x 1 的复导向矢量 theta = theta_deg * pi / 180; % 阵列位置:以第一个阵元为参考零点 pos = (0:M-1) * d_lambda; % 单位:波长 % 空间相位:2*pi * 位置(波长) * sin(入射角) phase = -2 * pi * pos.' * sin(theta); a = exp(1j * phase); end

这段代码的核心就一行:exp(1j * (-2*pi*pos*sin(θ)))pos是阵元位置向量,单位是波长,这样波数 * 位置 = 2π/λ * (m*d) = 2*pi*m*d/λ,直接得到相位。如果阵元间距不是均匀的,只需要把pos改成实际位置向量,这个函数依然成立——这是后面处理稀布阵或非均匀阵的基础。

2.3 扫描空间谱:把导向矢量当“模板”去匹配来波方向

有了导向矢量,波束形成器的输出功率就可以写成:

P(θ) = a^H(θ) * R * a(θ)

其中R = (1/N) * X * X^H是接收数据的协方差矩阵。对常规波束形成,这就是对各方向做匹配滤波。扫描整个 θ 空间,就能得到空间谱,谱峰位置对应来波方向估计。

MATLAB 源码里的扫描循环通常长这样:

% 假设 X 是 M x N 的接收数据矩阵(M个阵元,N个快拍) M = size(X, 1); R = (X * X') / size(X, 2); % 样本协方差矩阵 theta_scan = -90:0.5:90; % 扫描角度范围,步进 0.5 度 P = zeros(size(theta_scan)); for idx = 1:length(theta_scan) a = ula_steering_vector(M, 0.5, theta_scan(idx)); P(idx) = a' * R * a; % 输出功率(未归一化) end % 归一化为 dB 并绘图 P_dB = 10*log10(P / max(P)); plot(theta_scan, P_dB); xlabel('角度 (deg)'); ylabel('归一化功率 (dB)'); grid on;

这里有几个参数值得说明:theta_scan的步进决定了谱的角分辨率,步进设为 0.5 度在大多数线阵仿真里够用,但如果阵元数很多、波束很窄,建议改到 0.1 度。M是程序开头从数据矩阵行数自动读取的,不需要硬编码。0.5是半波长间距,后面会专门讨论这个值为什么是标配。

3. 平面阵实现:从一维指向到二维方位角和俯仰角估计

3.1 矩形面阵的阵元坐标与二维导向矢量

线阵只能分辨一个维度上的角度。当需求变成“同时给出目标相对于阵列的方位角和俯仰角”时,就需要面阵。源码包里的平面阵程序用的是矩形栅格布局,阵元在 x-y 平面内均匀排列。假设 x 方向有 Mx 个阵元、间距 dx,y 方向有 My 个阵元、间距 dy,那么每个阵元的坐标就是:

[xx, yy] = meshgrid((0:Mx-1)*dx_lambda, (0:My-1)*dy_lambda); xx = xx(:); % 拉成列向量 yy = yy(:); % 拉成列向量

meshgrid生成的坐标矩阵按行展开后,xxyy分别对应所有阵元的 x、y 坐标。总阵元数M = Mx * My

二维导向矢量由方位角 φ 和俯仰角 θ(有些文献里 θ 是阵列法线方向夹角,这里用俯仰角)共同决定。波数向量在直角坐标系下的分量为:

kx = (2π/λ) * sin(θ) * cos(φ)ky = (2π/λ) * sin(θ) * sin(φ)

导向矢量是对每个阵元位置(x, y)计算相位:

function a = planar_steering_vector(xx, yy, theta_deg, phi_deg) % 平面阵二维导向矢量 % 输入: % xx, yy - 阵元坐标,单位波长 % theta_deg - 俯仰角(与 z 轴夹角),单位度 % phi_deg - 方位角(与 x 轴夹角),单位度 % 输出: % a - 总阵元数 x 1 的复导向矢量 theta = theta_deg * pi / 180; phi = phi_deg * pi / 180; % 波数分量(以 2*pi/lambda 为基准,位置已归一化到波长) kx = sin(theta) * cos(phi); ky = sin(theta) * sin(phi); % 每个阵元的空间相位 phase = -2 * pi * (xx .* kx + yy .* ky); a = exp(1j * phase); end

逻辑跟线阵完全一致:xx .* kx + yy .* ky计算的是阵元在波传播方向上的投影距离,乘上 2π 就是相位。程序里如果阵元坐标不是以波长为单位,记得先除以波长,否则相位会差一个比例系数。

3.2 二维扫描的网格化处理与绘图

二维波束扫描就是遍历方位角和俯仰角的组合。源码里做的是双重循环,外层俯仰角、内层方位角:

theta_scan = 0:1:90; % 俯仰角扫描范围 phi_scan = -180:1:180; % 方位角扫描范围 P = zeros(length(theta_scan), length(phi_scan)); for ti = 1:length(theta_scan) for pi_idx = 1:length(phi_scan) a = planar_steering_vector(xx, yy, theta_scan(ti), phi_scan(pi_idx)); P(ti, pi_idx) = a' * R * a; end end % 取最大值做归一化 P_dB = 10*log10(P / max(P(:))); % 画三维图 surf(phi_scan, theta_scan, P_dB, 'EdgeColor', 'none'); xlabel('方位角 (deg)'); ylabel('俯仰角 (deg)'); zlabel('功率 (dB)');

注意俯仰角的扫描范围。如果你的坐标系里 90 度是阵列平面方向(即掠射方向),那么在接近 90 度处会看到严重的栅瓣和波束畸变,这不是程序 bug,而是阵列投影孔径趋于零的物理现象。实际使用中,俯仰角扫描范围不要超过 ±60 度。

surf画出来的三维谱图对论文插图够用,但想快速找峰时效率不高。更实用的做法是输出P矩阵后,用[max_val, idx] = max(P(:))找到峰值位置,再换算成对应角度:

[max_val, max_idx] = max(P(:)); [theta_peak_idx, phi_peak_idx] = ind2sub(size(P), max_idx); theta_est = theta_scan(theta_peak_idx); phi_est = phi_scan(phi_peak_idx); fprintf('估计角度: 俯仰角=%.2f°, 方位角=%.2f°, 峰值=%.2f dB\n', ... theta_est, phi_est, max_val);

4. 圆阵实现与周期阵列的栅瓣抑制策略

4.1 均匀圆阵的导向矢量:没有所谓“法线方向”

圆阵的几何特点是所有阵元均匀分布在半径为 R 的圆周上,没有优先方向,波束可以 360 度旋转。第 m 个阵元的角度位置为:

ψ_m = (m-1) * 2π / M

阵元直角坐标为(R*cos(ψ_m), R*sin(ψ_m))。圆阵和面阵的本质区别在于:它的孔径是圆对称的,这意味着波束形状在扫描到不同方位角时基本保持不变——这对全向监视应用很重要。

圆阵导向矢量的 MATLAB 实现:

function a = circular_steering_vector(M, radius_lambda, theta_deg, phi_deg) % 均匀圆阵导向矢量 % 输入: % M - 阵元数量 % radius_lambda - 阵列半径,单位波长 % theta_deg - 俯仰角(与 z 轴夹角),单位度 % phi_deg - 方位角,单位度 % 输出: % a - M x 1 的复导向矢量 theta = theta_deg * pi / 180; phi = phi_deg * pi / 180; % 阵元位置角度 psi = (0:M-1) * 2*pi / M; % 阵元直角坐标 x_pos = radius_lambda * cos(psi); y_pos = radius_lambda * sin(psi); % 入射方向单位矢量在 x-y 平面的投影 kx = sin(theta) * cos(phi); ky = sin(theta) * sin(phi); % 相位计算 phase = -2*pi * (x_pos * kx + y_pos * ky); a = exp(1j * phase); end

4.2 圆阵的栅瓣:阵元间距不是唯一变量

圆阵设计中一个容易踩的坑:阵元沿圆周均匀放置时,相邻阵元的直线距离不是半径,而是弦长:

d_chord = 2*R*sin(π/M)

如果 M=32、R=1.5λ,那么d_chord = 2*1.5*sin(π/32) ≈ 0.294λ,远小于半波长。这说明圆阵的阵元数可以做得比线阵多,而不必担心栅瓣。但前提是半径 R 的选择要匹配已验证的阵元数:

% 给定阵元数 M 和期望的相邻阵元间距 d_target,反推半径 d_target = 0.5; % 期望间距(波长) M = 32; R = d_target / (2 * sin(pi/M)); fprintf('建议半径: %.3f 波长\n', R);

这个反推公式在源码里没有直接给出,但做工程评估时很有用。如果把目标放宽容一些,允许间距到 0.8λ,那 R 还能更大,但代价是波束扫描到低俯仰角时会出栅瓣。判断依据很简单:对圆阵,阵元间最大相位差不仅取决于间距,还取决于 R 的大小,所以不要只看弦长。

5. 权值优化:从均匀加权到切比雪夫加窗的旁瓣控制

5.1 均匀加权的波束只有一个优点

均匀加权(所有权值幅度相同)的主瓣宽度最窄,在相同阵元数下具有最优的角分辨率。代价是第一旁瓣电平高达约 -13.3 dB。在很多雷达和通信场景里,这会引入来自旁瓣方向的干扰。加窗的目的就是降低旁瓣——权值幅度按一定规律从中心到边缘递减。

MATLAB 里的实现是在导向矢量上做元素级乘法:

% 30 元线阵,半波长间距,目标方向 20 度 M = 30; d_lambda = 0.5; theta_target = 20; % 计算导向矢量 a_target = ula_steering_vector(M, d_lambda, theta_target); % 加窗函数(列向量) wnd = chebwin(M, 30); % 切比雪夫窗,旁瓣电平 -30 dB % 加权后的波束权重 w = a_target .* wnd; % 画方向图(扫描各方向计算增益) theta_scan = -90:0.1:90; pattern = zeros(size(theta_scan)); for idx = 1:length(theta_scan) a = ula_steering_vector(M, d_lambda, theta_scan(idx)); pattern(idx) = abs(w' * a); % 注意这里取绝对值,画幅度方向图 end pattern_dB = 20*log10(pattern / max(pattern)); plot(theta_scan, pattern_dB);

w = a_target .* wnd这行代表了加窗的本质:导向矢量只负责相位对齐,决定波束指向;窗函数包络负责幅度锥削,决定旁瓣电平。两者是乘的关系,不是加的关系。chebwin(M, 30)里的 30 是旁瓣电平参数,单位是 dB。想要更低的旁瓣就把 30 改成 40,代价是主瓣会展宽——波束宽度和旁瓣电平存在固有的折中关系。

5.2 加窗对波束宽度和增益的影响量级

具体量化关系可以用一个简单实验说明。以 16 元线阵为例,工作频率对应半波长间距,目标方向 0 度。测得的参数对比如下:

窗函数第一旁瓣电平 (dB)3dB 波束宽度 (度)阵列增益损失 (dB)
均匀(矩形窗)-13.36.40
Hanning-31.58.91.2
Chebyshev (-30dB)-30.17.80.9
Chebyshev (-40dB)-40.09.11.8

从表格里能直接看到,当旁瓣要求从 -13 dB 压到 -40 dB 时,3dB 主瓣宽度从 6.4 度展宽到 9.1 度,增益损失 1.8 dB。如果设计指标同时要求窄波束和低旁瓣,靠加窗是做不到的,只能增加阵元数——这是物理孔径决定的。

Windows 里的chebwin函数对应的是 Dolph-Chebyshev 加权,它能在所有旁瓣电平均匀的条件下的给出最窄主瓣。在 MATLAB 环境里,直接调用即可。

5.3 一个常见误用:把窗函数加到数据上而不是权值上

我见过不少初学者把窗写成了:

x_win = x .* wnd; % 错:对数据加窗

这在频谱分析里是对的(对时域数据加窗抑制频谱泄漏),但在波束形成里不对。波束形成的加窗对象是权值,因为你是在控制阵列的幅度响应,而不是在改接收信号本身。如果对每个阵元的接收数据单独加窗,等效于引入了与方向无关的幅度扰动,会破坏导向矢量的相位结构。

6. 阵元位置误差对波束指向的影响及快速验证方法

6.1 位置误差导致的不只是增益下降

前面所有推导都假设阵元在精确位置上。实际阵列中,阵元安装位置会存在偏差——对一个半波长间距的 C 波段阵列来说,1% 的波长误差大约是 0.5 mm,别觉得这个量级微小。它会带来两部分影响:一是导向矢量失配,导致目标方向的增益下降;二是波束指向出现随机偏移,该偏移量在波束指向接近阵列法线附近时更为显著。

量化方式很简单:给阵元位置加一个零均值高斯扰动,然后分析波束方向图的统计特性:

% 均匀线阵,阵元位置扰动的影响分析 rng(42); % 可复现的随机种子 M = 32; d_lambda = 0.5; std_dev = 0.01; % 位置误差标准差(波长) % 理想阵元位置 pos_ideal = (0:M-1) * d_lambda; % 加扰动后的位置 pos_error = std_dev * randn(M, 1); pos_real = pos_ideal + pos_error; % 计算实际导向矢量(假设我们希望指向 30 度) theta_target = 30; a_ideal = exp(1j * (-2*pi*pos_ideal * sin(theta_target*pi/180))); a_real = exp(1j * (-2*pi*pos_real * sin(theta_target*pi/180))); % 增益损失(dB) gain_loss = 20*log10(abs(a_ideal' * a_real) / M); fprintf('位置误差 std=%d λ 时的增益损失: %.2f dB\n', std_dev, gain_loss);

std_dev = 0.01λ时,这个阵列的指向精度大约会退化到 0.3 度左右,增益损失只有零点几个 dB,看上去不严重。但把误差加大到0.05λ(四分之一阵元间距的公差内),增益损失会迅速恶化到 2 dB 以上,波束指向抖动超过 1 度。

6.2 校正技巧:在有源校正源下估计位置误差

一个实用的校正做法是:在远场放置一个已知方向的校正源(比如把发射天线架在楼顶,对阵列来说就是远场条件),记录各阵元的接收相位与理想相位的偏差,反推实际位置:

% 校正源方向已知为 theta_cal % 接收数据为 X_cal (M x N 快拍) % 从数据中估计各阵元相对参考阵元的相位差 X_centered = X_cal - mean(X_cal, 2); % 去均值 phase_measured = angle(mean(X_cal, 2) ./ X_cal(1,:)'); % 理论上接收相位应为 -2*pi * pos * sin(theta_cal) pos_estimated = -phase_measured / (2*pi*sin(theta_cal*pi/180)); % 位置误差 = 估计位置 - 设计位置 pos_error = pos_estimated - (0:M-1)'*d_lambda;

校正后,用估计出的实际位置替换理论位置去构造导向矢量,波束指向精度可以恢复到接近理想水平。这里的原理是:相位测量值中同时包含了位置误差和通道幅相误差,如果通道误差在之前已经校正过,那么剩下的相位偏差基本就是位置引起的。

这也是为什么我会建议在使用这个源码包的仿真结果投论文或造假数据之前,先跑一遍这段误差分析——它能告诉你一个方向图的旁瓣水平到底是真实的阵列性能,还是仅仅存在于理想仿真中的理论值。工程上更残酷一些:实测方向图与理论方向图的偏差,经常会重新定义你在论文里敢写的指标。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询