做阵列信号处理的同学多半都绕不开均匀直线阵Uniform Linear ArrayULA和常规波束形成Conventional Beamforming这两座入门的大山。原因很简单几乎所有更高级的算法——MVDR、MUSIC、ESPRIT甚至后来深度学习的DOA估计网络——底层都跑在同一个基础之上阵列模型和导向矢量。MATLAB作为最常用的仿真工具恰好能把这套理论从公式变成看得见摸得着的方向图。这篇文章我会从数学模型推导开始给出可直接复现的MATLAB代码再结合我调仿真时踩过的坑把几个关键参数的影响讲透目标是让看完的朋友能在自己的机器上跑出一张正确的波束方向图并且能理解每个参数为什么这么取。1. 均匀直线阵与常规波束形成先搞清楚在做什么1.1 什么是均匀直线阵它为什么重要均匀直线阵就是M个特性完全相同的阵元等间距地排成一条直线。这里“均匀”指的是阵元间距d固定不变每个阵元的响应也一致。这种阵列结构是阵列信号处理里最简单、最经典的一种也是理解其他阵列构型面阵、圆阵、共形阵的基础。它的重要性可以从实际物理场景去理解。假设你有一个单天线接收机它只能收到一个点上的电磁波场强没有空间分辨能力。但当你把多个天线按一定几何位置摆放并同时采样它们收到的信号就能利用波到达不同阵元的时间差或者说相位差来反推波来自哪个方向甚至能从强干扰中把弱信号提出来。这就是阵列带来的“空间增益”。均匀直线阵在雷达、声纳、无线通信、麦克风阵列、射电天文等领域都有广泛应用。比如5G基站的大规模天线阵列本质上就是很多个阵元按某种排布方式工作的阵列系统相控阵雷达则通过调整各阵元信号的相位让整个阵列的波束在空间快速扫描而无需转动天线。理解了均匀直线阵的波束形成再去看这些复杂系统就不会觉得隔着一层纸。1.2 常规波束形成的一句话原理延迟补偿再加权求和常规波束形成也叫延迟求和Delay-and-Sum波束形成、Bartlett波束形成。它的核心思想非常朴素既然同一个平面波前到达不同阵元的时间有先有后那我就把每个阵元的输出信号“补上”各自的延迟或者等效地补偿相位使所有阵元在某个期望方向上对齐再加权求和。用生活里的话说就像一排合唱团成员站在不同位置为了让某个方向的观众听到整齐划一的声音每个人需要根据自己的位置调整发声时间。调整好了从这个方向传来的声音会因相干叠加而增强而其他方向的信号因为相位没有对齐叠加时部分抵消形成相对抑制。这个操作在窄带信号信号带宽远小于载频条件下可以简化为“相位补偿”。因为窄带信号的包络变化很慢时间延迟主要体现为一个相位旋转所以不需要真的去估计时延直接用复数加权就能完成波束形成。这也是为什么接下来几乎所有公式里都围绕复数指数展开本质上都是在描述相位关系。1.3 角度、相位差和导向矢量几个绕不开的基础概念在推导公式之前需要固定坐标约定。通常把均匀直线阵放在x轴上第一个阵元放到原点法线方向定义为与阵列轴线垂直的方向也就是z轴或y轴方向。信号来波方向θ定义为与法线的夹角范围取-90°到90°正半平面通常是阵元编号增大的那一侧。当平面波以角度θ入射时相邻阵元之间的波程差为d·sinθ。如果阵元间距d和波长λ的比值固定那么相邻阵元的相位差就是k·d·sinθ其中k2π/λ是波数。也就是说第m个阵元相对第0个阵元的相位延迟为m·k·d·sinθ。导向矢量Steering Vector就是把这种相位关系写成向量的形式a(θ) [1, e^(j2πd·sinθ/λ), e^(j2π·2d·sinθ/λ), …, e^(j2π·(M-1)·d·sinθ/λ)]^T这个向量描述的是如果有一个方向为θ的单位幅度信号到达阵列各阵元采样到的复数值之比。它是整个波束形成算法的基石。指向一个方向的波束形成权向量就是该方向导向矢量的共轭窄带远场情况下。提示有些教材把角度定义为与阵列轴线的夹角那相位差就变成d·cosθ。两种定义都能用但仿真时混用会导致波束方向完全对不上。建议从头到尾统一用“与法线夹角”这个定义。2. MATLAB实现前的建模与核心公式推导2.1 阵列信号模型与窄带假设有了前面的物理图像现在把阵列接收信号写成矩阵形式。设有K个远场窄带信号从角度θ1, θ2, ..., θK入射阵列在某个快拍时刻的接收数据为x(t) A·s(t) n(t)其中x(t)是M×1的复数向量A [a(θ1), a(θ2), …, a(θK)]是M×K的阵列流形矩阵每一列是对应方向的导向矢量s(t)是K×1信号复幅度向量n(t)是M×1噪声向量通常假设为高斯白噪声。这里“窄带”假设非常关键。所谓窄带是说信号的带宽B远远小于载频也就是信号包络在跨越整个阵列所需的时间最大约M·d/c内几乎不变。这样波束形成只需要补偿相位而不必考虑包络延迟的影响。如果信号是宽带信号比如语音、宽带通信信号简单的相位补偿就不够需要做子带分解或者真时延补偿我在后面章节会提一下。在MATLAB里做这类仿真一般不用真的生成一个高频载波而是直接在复基带上做。因为线性系统的分析在中频/射频与基带是等效的直接建复基带模型可以省掉大量不必要的计算只看待关心的相位和幅度关系。2.2 导向矢量与阵列流形矩阵导向矢量的本质是一个空间采样的离散傅里叶变换核。你可以把“阵列对某个方向信号的响应”看成对空间中连续场在离散阵元位置上取样再按相位关系排列。导向矢量也就是波数域上的“频点”。在MATLAB里生成导向矢量非常简单关键是处理好角度到弧度的转换。以下是一个可以直接调用的函数function a steering_vector(M, d_lambda, theta_deg) % M: 阵元数量 % d_lambda: 阵元间距与波长之比 d/lambda % theta_deg: 来波方向与法线夹角单位度 q (0:M-1); % 阵元索引列向量 phase 2 * pi * d_lambda * q * sind(theta_deg); a exp(1j * phase); end注意这里用了MATLAB自带的sind()它直接把角度转成弧度再计算正弦比手动写sin(theta_deg*pi/180)更直观也不容易漏掉转换。我个人在早期写仿真时因为角度单位混淆吃过不少亏后来就养成习惯所有函数接口里角度一律用度内部只用弧度计算两者严格区分。如果需要一次生成多个方向的导向矢量矩阵可以这样写function A array_manifold(M, d_lambda, theta_list) % theta_list: 方向列表单位度行向量 q (0:M-1); A exp(1j * 2 * pi * d_lambda * q * sind(theta_list)); end输出的A是一个M×length(theta_list)的复矩阵每一列对应一个方向的导向矢量。这个向量化写法比for循环逐个角度生成要快一个数量级而且在后面做全角度波束扫描时非常方便。2.3 常规波束形成器与方向图的表达式设波束指向为θ0则权向量为w a(θ0)这里省略了共轭符号因为导向矢量本身就是共轭对称的直接用a(θ0)做权向量与实际相符。有的书会写成w (1/M)·a(θ0)多了一个1/M归一化因子这个因子不影响波束形状只影响输出幅度。阵列对某个方向θ的响应幅度就是 (|w^H a(θ)|)把θ在整个扫描范围内代入就得到方向图Beam PatternF(θ) |a^H(θ0)·a(θ)| |Σ_{m0}^{M-1} e^(j2πm·d·(sinθ - sinθ0)/λ)|这是一个等比数列求和可以写成闭式解F(θ) |sin(M·x/2) / (M·sin(x/2))|其中x 2πd·(sinθ - sinθ0)/λ这其实就是Dirichlet核的形式它的性质决定了常规波束形成方向图的主要特征主瓣出现在θθ0第一旁瓣约-13.3dB主瓣宽度随M增大而变窄。这些特性在后续仿真中都能直接看到。方向图通常有两种画法线性刻度幅度或对数刻度dB。实际工程中基本都看dB值因为旁瓣-13dB这种相对幅度差距在线性刻度上完全显示不出来。dB转换时用20·log10因为这里是电压/幅度响应不是功率响应。3. MATLAB代码实现从零开始写一个波束形成仿真3.1 最精简的版本方向图扫描先写一个最基础的版本目的是快速看到波束方向图。这个版本不需要真实信号只需要计算权向量和扫描范围内各角度的导向矢量然后做内积即可。%% 参数设置 clear; clc; close all; c 3e8; % 光速 (m/s) f 1e9; % 载频 1 GHz lambda c / f; % 波长 M 8; % 阵元数 d lambda / 2; % 阵元间距通常取半波长 theta0 30; % 波束指向角度度 %% 生成权向量 w steering_vector(M, d/lambda, theta0); % 常规波束形成权向量就是指向方向的导向矢量 %% 全角度扫描 theta_scan -90:0.1:90; % 扫描角度范围 A_scan array_manifold(M, d/lambda, theta_scan); % 所有角度的导向矢量矩阵M × N F w * A_scan; % 阵列响应复数1 × N %% 归一化并转dB F_dB 20 * log10(abs(F) / max(abs(F)) eps); %% 绘图 figure; plot(theta_scan, F_dB, LineWidth, 1.5); xlabel(扫描角度 (deg)); ylabel(归一化方向图 (dB)); title([ULA波束方向图, M, num2str(M), , dλ/2, 波束指向 , num2str(theta0), °]); grid on; ylim([-40 5]);运行这段代码你会看到一个在30°处达到峰值、峰值两侧依次下降的方向图。主瓣宽度大概十几度第一旁瓣大约-13dB附近方向图整体形状和前面公式预期一致。这里插一个细节abs(F)/max(abs(F))是归一化到最大值的电压响应加eps是为了防止在接近零点时log10(0)产生-inf导致绘图时曲线断裂到底部不好看。很多新手没加这个直接画结果方向图像锯齿一样断断续续实际上不是程序算错了只是log10的数值边界问题。3.2 带信号仿真的版本多个来波方向只有方向图还不够我们常常还想验证当多个方向的信号同时入射时波束形成器是不是真的能把指向方向的信号保留下来同时抑制其他方向的信号。下面这段代码生成了一个简单的信号环境并用波束形成器提取期望信号。%% 信号环境设置 clear; clc; close all; c 3e8; f 1e9; lambda c / f; M 8; d lambda / 2; theta_desired 40; % 期望信号方向度 theta_interf -20; % 干扰信号方向度 N 2048; % 快拍数 fs 10e6; % 基带采样率不必真实对应只是仿真用 t (0:N-1) / fs; % 两个窄带信号用复正弦表示频率不同便于区分 s1 exp(1j * 2 * pi * 100e3 * t); % 期望信号 s2 0.7 * exp(1j * 2 * pi * 150e3 * t); % 干扰信号幅度较弱 % 阵列流形矩阵 A [steering_vector(M, d/lambda, theta_desired), ... steering_vector(M, d/lambda, theta_interf)]; S [s1; s2]; % 合成阵列接收数据不考虑噪声 X A * S; % M × N % 加入复高斯白噪声信噪比大约 10dB snr 10; noise_power mean(abs(X(:)).^2) / (10^(snr/10)); noise sqrt(noise_power/2) * (randn(M, N) 1j*randn(M, N)); X X noise; %% 常规波束形成指向40° w steering_vector(M, d/lambda, theta_desired); y w * X; % 1 × N %% 对比恢复信号与原始期望信号 % 计算两个已知频率处的幅度 s1_est_amp abs(mean(y .* conj(s1))) * 2 / N; % 近似提取幅度 s2_est_amp abs(mean(y .* conj(s2))) * 2 / N; % 近似提取幅度 fprintf(期望信号幅度理论 1.00波束输出 %.3f\n, s1_est_amp); fprintf(干扰信号幅度理论 0.70波束输出 %.3f\n, s2_est_amp);这个例子直接给出了“波束形成后哪个信号被保留、哪个被抑制”的量化结果。由于方向图在40°增益最大干扰方向落在旁瓣区被压掉不少所以y中主要是期望信号成分。实际跑一下你会发现干扰抑制的幅度大概率在10dB上下和方向图上-20°相对40°的增益差基本吻合。这里面有个工程上的小习惯验证波束形成效果时不要只看时域波形直接计算输出中各个频率分量的幅度一眼就能看出增益关系比肉眼比较波形高效得多。3.3 参数封装成函数方便反复调用写仿真时最忌讳每改一次参数就复制粘贴一大段代码。我习惯把常规波束形成封装成函数参数变化只动函数输入这样换场景、做参数扫描都非常快。function [P_dB, theta_scan] cb_ula_direction_pattern(M, d_lambda, theta0) % 常规波束形成方向图计算 % 输入 % M: 阵元数 % d_lambda: 阵元间距与波长之比 % theta0: 波束指向度 % 输出 % P_dB: 归一化方向图dB-90°到90°步进0.1° % theta_scan: 对应的角度轴 theta_scan -90:0.1:90; w steering_vector(M, d_lambda, theta0); A_scan array_manifold(M, d_lambda, theta_scan); F w * A_scan; P_dB 20 * log10(abs(F) / max(abs(F)) eps); end然后再写一个测试脚本比如对比不同阵元数下的方向图figure; hold on; for M_i [4 8 16 32] [P_dB, theta_scan] cb_ula_direction_pattern(M_i, 0.5, 0); plot(theta_scan, P_dB, DisplayName, [M, num2str(M_i)]); end xlabel(扫描角度 (deg)); ylabel(归一化方向图 (dB)); legend; grid on; ylim([-40 5]);参数扫描只需要改一个数字整个对比图就出来了。这个封装思路也方便后续扩展到MVDR、MUSIC等算法因为它们的输入结构几乎一样换的只是权向量的计算方式。提示封装函数时尽量让角度处理统一比如内部固定用sind()、固定用“与法线夹角”定义。一旦同一个工程里出现两种角度定义调试到怀疑人生的概率会翻好几倍。4. 仿真结果怎么看主瓣、旁瓣、栅瓣与关键参数的影响4.1 阵元数M如何决定波束宽度阵元数M是均匀直线阵最直接的“资源”它直接决定孔径长度。对于半波长间距的均匀直线阵阵列物理孔径为(M-1)·λ/2近似等于M·λ/2。孔径越大波束越窄空间分辨能力越强。3dB波束宽度半功率波束宽度的近似公式为BW_3dB ≈ 0.886λ/(M·d·cosθ0)弧度代入dλ/2、θ00°可以得到BW_3dB ≈ 101.5°/M。也就是说M8时波束宽度约12.7°M16时约6.3°M32时约3.2°。这个是“近似值”实际从方向图上量测会略有偏差但趋势非常稳。需要注意波束指向偏离法线后等效投影孔径变小波束会展宽。指向60°时cos60°0.5波束宽度翻倍。这就是为什么相控阵雷达在超大扫描角度下性能下降通常只扫描±60°以内的区域超出后主瓣太宽、增益损失太大。你在MATLAB里把theta0分别设为0°、30°、60°跑一遍方向图能直观看到主瓣逐渐变宽的过程。旁瓣高度几乎不随M变化第一旁瓣始终在-13.3dB左右。这是常规波束形成的固有缺点旁瓣电平较高意味着强干扰或杂波可能从旁瓣进入。如果想压低旁瓣需要在权向量上加窗函数切比雪夫窗、汉明窗等但代价是主瓣变宽。这也是经典波束形成器的性能天花板理解这一点对后面学自适应波束形成很有帮助。4.2 阵元间距d与栅瓣为什么默认取λ/2阵元间距d是另一个关键参数。当d过大时方向图会在多个方向出现与主瓣等高的峰这些多余的峰称为栅瓣Grating Lobe。栅瓣的出现是因为阵列采样的空间相位出现了混叠——不同方向角可能产生相同的阵元间相位差波束形成器无法区分它们。空间采样和时域采样完全可以类比时域采样率不够会导致频谱混叠空间采样间距太大同样会导致“空间谱混叠”。时域采样的奈奎斯特条件是采样间隔小于半个周期空间采样的对应条件就是阵元间距d ≤ λ/2。在扫描范围±90°、波束指向0°的情况下只要d/λ ≤ 0.5就不会出现栅瓣。如果波束指向大角度条件会更苛刻。更精确的判断方法是看是否存在除θ0之外的另一个角度θg使得sinθg - sinθ0 ±λ/d成立。满足这个关系的方向就会出现栅瓣。MATLAB里验证栅瓣特别直观把d改成1.2λ再画方向图你会看到除了主瓣外-90°和90°附近都出现了明显的等幅峰。这时候如果真实信号从这些栅瓣方向进来会被当成主瓣方向的信号产生严重的方向模糊。在工程中比如稀疏阵列设计的初衷之一就是为了在保证不出现栅瓣的前提下尽量分散阵元位置从而获得更大的等效孔径。4.3 波束指向与扫描范围从-90°到90°的注意事项常规波束形成的一个天然优势是可以通过改变theta0实现电扫描不需要任何额外运算。MATLAB里做波束扫描本质上就是循环更新权向量计算每个指向角的方向图或者输出功率。但有几个细节容易踩坑一是扫描范围。对于一维均匀直线阵我们通常只扫描-90°到90°的前向半空间。理论上因为阵列在正反两个方向性能对称来自后方的信号也可能被波束捕获所以有的应用需要扫描到-180°到180°。这时候就需要区分“法线方向”和“端射方向”的定义以及是否需要使用双面阵元。在我的经验里一上来就扫±180°的仿真十有八九画出来是左右对称的两个主瓣容易被误判成程序错误其实那是阵列前后镜像响应的正常表现。二是扫描步进。步进太大比如5°会导致方向图峰值不够尖锐波束指向附近出现锯齿步进太小比如0.01°在参数扫描频繁运行时会让计算量成倍增加。一般选0.1°到1°就足够。三是波束指向30°时权向量本身就是在该方向归一化的导向矢量方向图主瓣最大值理应为10dB。如果发现最大值不在theta0处优先检查角度是否转成了弧度、是否用了sind和sin混用、以及d/lambda是否传错。5. 常见问题与调试经验我踩过的坑5.1 方向图两端开裂或出现“伪峰”怎么办方向图在扫描范围两端-90°或90°异常爬升通常有两个原因。第一是阵元间距超过λ/2导致栅瓣靠近边缘第二是本身没有异常但两端因为投影孔径减小波束变宽归一化后看上去会有抬升趋势。实际上在90°附近即使没有栅瓣方向图也会比中心区域更“钝”。排查时先把d改成λ/2、theta0改成0°如果曲线恢复正常的尖峰旁瓣结构说明参数设置没问题如果依然有额外的与主瓣等高的峰说明栅瓣因素没有消除需要检查d/lambda是否确实传入了0.5。另一个常见做法是把扫描范围暂时缩到-60°到60°重画排除边缘效应干扰。5.2 dB归一化后曲线断裂这个问题太常见了几乎是新手必踩。方向图在很多角度上响应接近零取log10之后会得到-Infplot时曲线直接掉出画图窗口看上去像断裂。解决方案就是我在代码里演示过的那一句20*log10(abs(F)/max(abs(F)) eps)。加eps本质上就是给归一化响应一个极小的地板不让log10的参数取到0。如果你希望动态范围更合适也可以用max(dir, 1e-6)之类的操作。还有一个小技巧画图时把ylim设置在[-40 5]或者[-60 5]低于这个电平的部分直接截断可以让曲线更干净。5.3 多快拍数据的处理不要一次只循环一个快拍我在早期写多快拍仿真时犯过一个低级错误用for循环逐快拍做波束形成再把所有快拍的输出存下来代码又慢又冗长。后来才想明白MATLAB的矩阵运算天然适合批量处理快拍。正确的做法是把接收数据X组织成M×N矩阵N为快拍数波束形成的权向量w是M×1那么波束输出就是y w * X一次性得到1×N的输出序列。如果要对多个波束指向同时处理可以组织成W为M×L的权矩阵输出Y W * X得到L×N的输出矩阵。整套运算完全向量化不仅代码简洁而且效率提升几个量级。类似的计算波束输出功率时可以先算样本协方差矩阵R (X*X)/N再算P w*R*w。这种方法在做波束扫描时特别高效因为R只用算一次后面每个角度只需要一次矩阵乘法。5.4 关于角度的两种写法sin、sind与deg2rad这是一个很小但影响很大的编码习惯问题。MATLAB同时提供了sin、sind、deg2rad这些函数用错一个方向图就直接错乱。比如sin(90)会算的是90弧度的正弦不是90度。我的习惯是所有自定义函数的输入输出角度一律用“度”在函数内部计算相位时统一用sind()或deg2rad()转成弧度。这样做的好处是脚本里到处传的都是30、40、-20这样的整数角度好读也好改。如果你喜欢用sin(theta*pi/180)也没有问题关键是全工程保持一致不要这一行用sind下一行用sin还传了度数这种错误用了好多年编程的人也偶尔会犯。这里再提供一个检查技巧如果画出来的方向图主瓣位置和设定值差了很大一截先别怀疑算法先检查所有角度的转换是否一致。我见过太多因为这个问题而反复调试半天的案例了。5.5 版本与工具箱哪些函数需要额外安装本文里的代码只用了MATLAB基础功能矩阵运算、复数运算、绘图不依赖任何特定工具箱所以不管是R2018b还是更新的版本都能直接跑。唯一建议是确保绘图时用了legend、grid这类基础函数它们在所有版本里都稳定。如果你后续想做更复杂的扩展比如用优化工具箱设计最优权向量或者用信号处理工具箱做滤波器设计就要检查当前环境是否安装了对应工具箱。判断方法很简单在命令行输入ver查看已安装的工具箱列表即可。这里顺便提一句MATLAB官方对教育和科研用户有home版和校园授权正版软件对学生来说其实有比较实惠的渠道不太建议去网上找来路不明的破解版本既不稳定也容易遇到版权问题。6. 从常规波束形成往后走值得继续深入的方向写到这里均匀直线阵的常规波束形成已经从原理到MATLAB仿真都覆盖到了。但我个人觉得入门归入门真正把这块吃透之后有几个方向非常值得沿着往下走。第一个方向是自适应波束形成典型算法就是MVDRMinimum Variance Distortionless Response。它和常规波束形成的核心区别在于权向量不再单纯取导向矢量而是在保持期望方向增益为1的约束下最小化输出功率从而自动在干扰方向形成零陷。效果很惊艳但需要注意信号模型失配时性能会急剧下降这也是学界和工程界研究多年的问题。第二个方向是波达方向估计也就是DOA估计算法代表算法包括MUSIC和ESPRIT。它们的思路和波束形成完全不同波束形成是在“猜测方向”上扫描看哪个方向的能量强MUSIC则是利用信号子空间与噪声子空间的正交性在角度谱上做超分辨搜索。这部分内容对矩阵特征分解的要求更高但MATLAB实现起来依然很顺手。第三个方向是把模型推广到更复杂的阵列构型。比如矩形面阵可以同时估计方位角和俯仰角均匀圆阵能实现360°全方位覆盖共形阵贴近载体表面但流形计算更复杂。这些场景下虽然导向矢量的表达式变了但常规波束形成的基本逻辑——补偿相位、加权求和——是完全一致的。还有一点如果信号带宽比较大窄带相位补偿就不再准确需要考虑宽带波束形成。简单来说就是把宽带信号分到若干子带每个子带单独做窄带波束形成最后再把输出合起来。MATLAB里可以做这种多子带处理也可以用滤波器组实现真时延补偿思路都是相通的。我自己的切身体会是均匀直线阵和常规波束形成这套基础内容看起来公式不多、代码也不长但它把阵列信号处理里最核心的几个物理概念——空间采样、相位差、方向图、主瓣旁瓣、栅瓣——全都串联了起来。后面学再多高级算法回归到本质还是在和这些概念打交道。这个仿真做熟练之后再去看相控阵雷达的扫描控制、5G天线阵列的用户对齐、甚至声学麦克风阵列的语音增强都会有一种“原来都是同一套底层逻辑”的感觉。最后再分享一个小技巧做参数扫描时比如同时扫描M和d的取值建议在脚本开头用clear; clc; close all三件套并在代码里用DisplayName给每条曲线命名配合legend导出对比图这会大幅提升你分析参数影响的效率。