常规波束形成MATLAB实现：原理、参数与避坑指南

简介常规波束形成常规聚焦波束形成MATLAB 程序资源面向阵列信号处理与水声/声呐方向的学习者用于完成声源定位仿真、波束扫描成像与声图可视化。资源包内共 1 个 Word 文档约 29KB文档内含完整的 MATLAB 代码实现与步骤说明目前已吸引 1678 人学习浏览。代码覆盖了参数设定、阵元坐标生成、声程差与相位延迟计算、加噪信号建模、自协方差矩阵构建、扫描聚焦矢量计算、归一化处理、mesh/pcolor 声图绘制以及峰值与主瓣尺寸提取等环节结构清晰且参数可调。使用者可修改声速、阵元数、阵元间距、声源位置、信噪比和扫描步长直观观察不同条件对波束形成定位结果的影响既能帮助初学者理解常规聚焦波束形成的基本原理也能为声呐、超声成像或麦克风阵列定位实验提供可直接运行的参考模板。1. 常规波束形成16 阵元 MATLAB 程序到底在算什么常规波束形成Conventional Beamforming也叫延迟求和波束形成是阵列信号处理里最朴素也最不容易出错的一种方案。做过麦克风阵列或水声阵列的工程师应该都有印象当系统里只有 8 个阵元、目标大致在 30° 方向时用一组导向矢量把各路信号在相位上对齐再叠加输出信噪比能直接提升约 9 dB加到 16 阵元就有 12 dB 的阵列增益。它不需要协方差矩阵求逆不依赖迭代收敛唯一要算准的是阵元间距、波长和角度定义。这篇笔记从均匀线阵出发给出一套不依赖工具箱、可直接运行的 MATLAB 程序并把参数设计逻辑和常见翻车点写清楚适合刚接触阵列处理、想快速跑通验证的从业者。2. 延迟求和与导向矢量常规波束形成的原理里最该吃透的 3 个公式2.1 为什么时延补偿能在噪声里把信号捞出来考虑一个 N 元均匀线阵阵元间距为 d一个平面波从 θ 方向入射。相邻阵元之间多走的距离是 d·sinθ对应的时间差是 d·sinθ / c。如果我把每个阵元的信号按照这个时间差向前补偿再全部加起来目标方向的信号会变成同相叠加幅度变成原来的 N 倍功率变成原来的 N² 倍而噪声在阵元之间是独立的叠加后功率只变成原来的 N 倍。于是输出信噪比相比单个阵元提升了 N 倍换算成 dB 就是 10·log10(N)。这里有一个容易被新入门的开发者忽略的点阵列增益提升的是输出 SNR不是方向分辨能力。16 个阵元能把输入 SNR 抬高 12 dB但角度分辨率只和阵列有效孔径有关。在实际工程里我通常把这个 N 倍增益当成整个波束形成链路的第一项验收指标——仿真里如果目标方向输出功率没有比单阵元噪声底高出大约 10log10(N) dB那多半是噪声建模或者导向矢量符号出了问题先别急着往下调参数。2.2 窄带信号如何用相位代替时延对窄带信号时延可以用相位旋转来近似。假设信号是中心频率为 f0 的窄带信号 s(t)那么 s(t - τ) 在频域近似等于 s(t)·exp(-j·2π·f0·τ)。也就是说补偿时延变成了补偿一个与频率和时延有关的相位。对均匀线阵从 θ 方向入射的信号在第 n 个阵元上的相位延迟是φ_n -2π · d · n · sin(θ) / λ (n 0, 1, ..., N-1)把每个阵元的相位补偿因子写成一个列向量就是导向矢量steering vectora(θ) [1, exp(-j·2π·d·sin(θ)/λ), ..., exp(-j·2π·d·(N-1)·sin(θ)/λ)]^T波束形成的权向量直接取 w a(θ)输出 y w^H · x。这里的 H 表示共轭转置。符号约定必须前后一致我在程序里统一采用「信号沿阵列传播方向滞后为负相位」的约定补偿时就取共轭。如果你在别的教程里看到符号相反不要慌它们只是把入射方向定义成相反波束图会关于 0° 翻转但峰值的绝对值不变。2.3 波束图和主瓣宽度阵列孔径说了算波束图的定义很简单把权向量固定指向目标方向 θ_s然后扫描所有可能的 θ看输出幅度随 θ 的变化。对均匀线阵解析表达式可以写成B(θ) |Σ_{n0}^{N-1} exp(-j·2π·d·n·(sinθ - sinθ_s)/λ)| / N这个式子展开后形式上很像离散傅里叶变换所以主瓣宽度和阵元数 N、阵元间距 d 的关系一目了然。半功率波束宽度大约为θ_3dB ≈ 0.886 · λ / (N·d·cosθ_s) 弧度当 N16、dλ/2、θ_s0° 时算出来约 6.3°。注意两个容易被忽略的结论第一主瓣宽度由阵列孔径N·d决定而不是由阵元数单独决定第二波束指向偏离法线后cosθ_s 变小主瓣会变宽。所以我在做水声目标探测时经常提醒自己不要只看 0° 法线方向的波束图30° 以上扫描角的主瓣宽度要按实际指向重新估算。另一个重要参数是第一旁瓣水平。矩形加权也就是不加窗的均匀线阵第一旁瓣固定约 -13.3 dB这个值不随阵元数增加而改善。如果你想压低旁瓣只能靠幅度加权加窗代价是主瓣变宽。这个取舍在第四章会展开。2.4 常规波束形成与自适应波束形成的定位差异常规波束形成在空间谱估计里对应的是 Bartlett 谱P(θ) a^H(θ)·R·a(θ)其中 R 是接收数据的协方差矩阵。它的特点是完全不需要对 R 求逆即使快拍数很少、R 秩亏也能算出结果只是曲线毛刺多一点。与之相对的 MVDR最小方差无失真响应波束形成也叫 Capon 法空间谱里带 R^(-1)能实现超分辨把两个靠得很近的目标分开。但自适应方法有一个绕不开的前提导向矢量必须足够精确协方差矩阵估计要稳定否则性能会急剧下降。实际阵列里阵元位置误差、幅相不一致、互耦都会破坏这个前提。我自己的经验是第一版算法永远先用常规波束形成建立基线确认信号模型和数据通路没问题再考虑上自适应。这个习惯帮我避开了很多「仿真很漂亮、实测翻车」的情况。3. 常规波束形成 MATLAB 最小程序均匀线阵从参数到波束图3.1 参数初始化声速、频率、阵元数与阵元间距先上参数初始化代码这段是整个程序的地基。我用声呐场景常用的参数来写但改成雷达或语音阵列只需换 c、f、lambda 三行。%% 参数初始化 c 1500; % 介质声速单位 m/s水声场景取 1500 f 10000; % 窄带信号中心频率单位 Hz lambda c / f; % 波长 0.15 m N 16; % 阵元数 d lambda / 2; % 阵元间距取半波长 theta_s 30; % 目标真实来波方向单位度 theta_scan -90:0.5:90; % 角度扫描范围单位度 snap 2000; % 快拍数 SNR_dB 10; % 单阵元输入信噪比单位 dBlambda 由 c/f 计算而不是直接写 0.15是为了让你改频率时波长自动更新避免手算出错。d 取 lambda/2 是常规波束形成的标准选择原因在第四章细说。theta_scan 用 0.5° 步长对这个程序对应的 6.3° 主瓣宽度来说主瓣内至少能采到 12 个点画出来足够平滑。3.2 接收信号建模导向矢量、目标信号与复高斯噪声接收信号仿真用窄带模型。这里最关键的是噪声要生成复高斯白噪声实部虚部方差各 0.5总噪声功率才是 1。很多人直接写 randn(N, snap) 当成噪声功率其实是 1但后面信噪比缩放时容易把信号和噪声的功率尺度搞反。%% 导向矢量生成函数内联形式 n_idx (0:N-1).; % 列向量阵元索引 a_s exp(-1j * 2 * pi * d * n_idx * sind(theta_s) / lambda); %% 接收信号仿真 s (randn(1, snap) 1j * randn(1, snap)) / sqrt(2); % 复基带信号功率为 1 noise (randn(N, snap) 1j * randn(N, snap)) / sqrt(2); % 复噪声功率为 1 X sqrt(10^(SNR_dB / 10)) * a_s * s noise; % N x snap 接收矩阵这段的数学含义信号是秩为 1 的矩阵 a_s * s每个阵元收到同一信号的不同相位功率为 SNR_dB 对应的线性值噪声每个阵元独立同分布功率为 1。窄带假设下信号在各阵元上只有相位差、没有幅度差这是常规波束形成建模的核心简化。如果实际信号是宽带这个假设不成立第五章会讲对应的问题。信噪比缩放我更喜欢放在信号上而不是噪声上这样噪声功率固定为 1所有后续统计对比都有统一基准。改动 SNR_dB 只需动一行波束图输出功率的绝对值变化能直接反映输入 SNR 的影响。3.3 Bartlett 波束扫描for 循环里的共轭转置空间谱扫描的核心代码只有几行。用循环而不是矩阵运算一次性算完是为了让每一位读代码的人都能看清每次迭代在做什么生成当前扫描角度的导向矢量、做两次内积、得到该角度的输出功率。%% Bartlett 波束扫描 R X * X / snap; % 样本协方差矩阵N x N P zeros(size(theta_scan)); for i 1:length(theta_scan) a_scan exp(-1j * 2 * pi * d * n_idx * sind(theta_scan(i)) / lambda); P(i) real(a_scan * R * a_scan); % 该角度的输出功率 end %% 归一化并绘图 P_dB 10 * log10(P / max(P)); figure; plot(theta_scan, P_dB, LineWidth, 1.5); xlabel(扫描角度 (deg)); ylabel(归一化输出功率 (dB)); title(常规波束形成波束图Bartlett 谱); grid on;a_scan * R * a_scan 的逻辑是先用共轭转置把导向矢量取共轭对应时延补偿然后与协方差矩阵相乘最后再右乘导向矢量对应加权求和。real() 是为了去掉浮点运算带来的极小虚部不影响结果。R 是 N×N 的样本协方差矩阵用 X * X / snap 得到snap 远大于 N 时估计比较稳定。如果快拍数小于阵元数常规波束形成仍然能算只是谱线毛刺增多这正是它相比自适应方法的稳健之处。3.4 完整脚本运行后你该看到什么样的波束图把上面三段拼起来就是一个完整可运行的脚本MATLAB 里不需要额外安装任何工具箱Phased Array System Toolbox 在这里是可选的不是必需的。运行后你应当看到峰值出现在 30°也就是 theta_s 对应的位置主瓣宽度大约 6°7°第一旁瓣在 -13 dB 附近远离主瓣的底噪整体低于 -20 dB 并伴随随机起伏。如果峰值不在 30°优先检查 sind 是不是写成了 sin以及 a_s 和 a_scan 的符号约定是否一致。这是全部程序里最容易出问题的地方也是下一章要展开的第一个避坑点。我来补一个运行性能的参考snap2000、扫描 361 个角度整个循环在普通笔记本上跑完不到 1 秒。这意味着你完全可以把扫描步长加密到 0.1° 来观察主瓣细节算力不是瓶颈瓶颈反而是角度扫描范围过大时栅瓣会不会混进来。4. 常规波束形成参数怎么定N、d、快拍数与加窗的取舍4.1 阵元数 N主瓣宽度与 10log10(N) 阵列增益阵元数 N 的选取直接决定两个量主瓣宽度和阵列增益。主瓣宽度大约正比于 1/(N·d)阵列增益是 10·log10(N) dB。把 N 从 8 增加到 16主瓣宽度减半增益增加 3 dB从 16 增加到 32再减半、再加 3 dB。看起来是线性收益但实际系统里 N 对应着 AD 通道数、阵元位置校准成本、数据吞吐量每翻一倍都意味着硬件成本近似翻倍。N 的典型值是 8、16、32。做声呐或超声成像时我一般从 16 起步因为单阵元 SNR 通常不高12 dB 增益能把目标从噪声底里抬出来做语音麦克风阵列48 阵元也够用因为语音 SNR 本身不差更在意的是指向性和体积。程序里把 N 单独放在一行就是方便你在不同项目间快速对比——改一个数重新跑一遍看波束图主瓣变化这是最直观的选型依据。4.2 阵元间距 d半波长是上限还是推荐值均匀线阵的阵元间距标准是 λ/2。这个值有明确的数学背景d λ/2 时相邻阵元的相位差为 π·sinθsinθ 的取值范围 [-1,1] 对应相位差范围 [-π, π]整个可见区域没有相位模糊。这个条件下常规波束形成在 -90° 到 90° 全范围扫描不会出现栅瓣。把 d 加大到 λ相位差范围变成 [-2π, 2π]当 sinθ ±0.5即 ±30°时相位差恰好 ±π与 sinθ ∓0.5 的相位差无法区分波束图出现一对假峰。这个现象称为栅瓣是阵列信号处理里最容易踩的坑之一。反过来d 缩小到 λ/4栅瓣不会出现但主瓣宽度按比例变宽阵列孔径小了分辨率下降。所以 d 的最优值不是越小越好λ/2 是分辨率和无模糊条件之间的自然平衡点。4.3 快拍数与信噪比让仿真相差不太远的两种设定快拍数 snap 影响的是样本协方差矩阵 R 的估计质量。snap 越大R 中信号子空间和噪声子空间的分离越清晰波束图背景越平滑。经验上 snap 至少取 10 倍阵元数我一般会取到 50 倍以上。程序里的 2000 对 16 阵元来说冗余度足够当你把 N 加到 64 时snap 建议同步上调到 5000 以上否则波束图旁瓣会出现明显的随机毛刺。仿真里还有一个值得养成的习惯用 rng(0) 固定随机数种子。没有固定种子时每次运行噪声不同主瓣附近会出现 12 dB 的随机起伏你可能误以为是参数改动引起的花大量时间排查。固定种子后所有对比都在同一噪声现实下进行可复现性大大提高。等仿真逻辑完全稳定了再去掉种子、用蒙特卡洛多次运行统计平均性能那是更严谨的做法。4.4 hamming 窗压旁瓣主瓣拓宽一倍的代价前面说过矩形加权的第一旁瓣是 -13.3 dB。如果你的系统有强干扰比如 30° 目标旁边 40° 有一个功率高 20 dB 的干扰源-13dB 的旁瓣根本压不住干扰这时需要给阵列加幅度窗。加窗的 MATLAB 做法是在导向矢量上乘一个窗函数而不是对接收数据直接做乘法——对窄带波束形成来说权向量是窗函数和导向矢量的逐元素乘积win hamming(N); % 或 hanning(N)、blackman(N) w_scan win .* a_scan; % 加窗后的权向量 P_win(i) real(w_scan * R * w_scan);加窗的代价是主瓣变宽。hamming 窗的主瓣宽度约为矩形窗的 1.3 倍第一旁瓣降到约 -43 dBhanning 窗主瓣更宽一些旁瓣约 -31 dBblackman 窗旁瓣能到 -58 dB但主瓣宽度接近矩形窗的两倍。如果你要同时获得低旁瓣和窄主瓣正确的方向是增加阵元数 N而不是纯靠换更陡的窗函数。这个取舍用一句话记窗函数是用主瓣宽度换旁瓣高度N 才是决定绝对分辨率的东西。5. 常规波束形成 MATLAB 避坑5 个把结果带偏的细节5.1 栅瓣d 超过半波长时出现的假目标现象目标设置在 30°波束图在 -30° 附近出现一个幅度几乎相同的峰值怎么调 SNR 都消不掉。原因阵元间距 d 大于 λ/2sinθ 相位超出 [-π, π]不同角度对应相同相位差出现方向模糊。解决检查 d 的赋值确认它是按当前频率的波长计算的。如果必须用固定物理间距覆盖多频段要在程序里按最高工作频率约束 d ≤ λ_min/2并在扫描前加一个合法性判断if d lambda / 2 warning(阵元间距超过半波长可能出现栅瓣); end5.2 sind 与 sin 混用角度一超过 45° 就乱套现象峰值位置在 0° 附近时还算准目标改到 60° 时峰值偏了 2°3°改到 90° 时完全错位。原因sind(theta) 接收的是角度值sin(theta) 接收的是弧度值。程序里 theta_scan 是角度向量写成 sin(theta_scan) 后90° 被当成 90 弧度处理sin(90) ≈ 0.894对应真实方向约 63°当然找不到目标。解决全程序统一用 sind尤其是导向矢量和解析波束图里不要在一个文件里混用两种如果从外部配置文件读入的是弧度用 rad2deg 明确转换。5.3 宽带信号套窄带模型阵列增益直接消失现象用一段语音或超声回波直接代入程序发现输出信号像被低通滤波了一样目标方向增强不明显波束图在扫频后完全变形。原因窄带模型的相位补偿 exp(-j·2π·f·τ) 依赖频率宽带信号里每个频率分量对应的导向矢量不同单一频率的补偿只能对齐一个频点。解决先做短时傅里叶变换把每个频点单独按窄带波束形成处理再合回时域。具体来说对每个频带 bin 用该 bin 中心频率对应的波长重新生成导向矢量扫描同一方向最后把所有频带的输出求和。这也是实战中处理语音阵列的标准套路窄带程序是它的地基。5.4 SNR 定义倒挂改参数后波束图纹丝不动现象SNR_dB 从 0 改到 20波束图峰值位置不变但旁瓣和底噪完全没有变化甚至峰值幅度也没明显增加。原因信噪比缩放被写在了噪声上或者信号幅度用了线性值但没换算 dB。常见写法是把 noize 乘以 10^(SNR/20)方向弄反了SNR 越高噪声越大当然看不出规律。解决统一噪声功率为 1信号幅度按 sqrt(10^(SNR_dB/10)) 缩放输出功率用 10*log10(P) 显示。每次调完 SNR 后检查主瓣峰值处的绝对功率是否近似等于线性 SNR 乘以 N²验证信号链路本身是对的。5.5 0° 定义不一致峰值跑到阵列端射方向现象目标明明放在阵列法线方向程序却在 ±90° 找到峰值。原因坐标系定义不一致。我常用的约定是 0° 为阵列法线方向阵列面朝的方向±90° 为端射方向信号沿阵列轴向入射。但有些算法包默认 0° 是端射导致导向矢量里 sinθ 的映射关系整体偏移了 90°。解决在程序头部用注释写明坐标系% 坐标约定均匀线阵沿 x 轴排列0° 为法线方向y 轴±90° 为端射方向然后让目标角度、扫描范围和实际物理摆放保持一致。排查时先用 0° 目标做自检峰值如果不在 0°一定是坐标约定或者符号反了。6. 验证与扩展用解析波束图给数值程序做体检6.1 解析波束图与数值扫描叠加对比数值程序有一个不变的好习惯用解析解交叉验证。常规波束形成的解析波束图就是导向矢量加权求和MATLAB 里可以直接画出理论曲线theta_th -90:0.01:90; B_theory abs(sum(exp(-1j * 2 * pi * d * n_idx * (sind(theta_th) - sind(theta_s)) / lambda), 1)) / N; B_theory_dB 10 * log10(B_theory.^2 / max(B_theory.^2)); hold on; plot(theta_th, B_theory_dB, r--, LineWidth, 1.2); legend(数值扫描, 解析波束图);数值曲线和解析曲线应当几乎重合。如果主瓣宽度对不上先查阵元数和 d 是否一致如果旁瓣电平对不上查是不是忘了除以 N如果整体偏移查 sind/sin 和符号约定。这一步通过后这个程序才有资格作为后续所有算法的参照基线。6.2 从均匀线阵到平面阵的扩展套路均匀线阵只在一个维度上分辨角度实际系统中更常见的是平面阵或圆阵。平面阵的扩展思路并不复杂阵元位置从一维坐标变成二维坐标 p_n [x_n, y_n]导向矢量变成 exp(-j·k^T·p_n)其中 k 是波数矢量由方位角和俯仰角共同决定。波束形成的主循环结构完全不变变的只是导向矢量生成函数。我做这类扩展时习惯先把线阵程序验证干净再改导向矢量函数这样如果出了问题问题一定只出现在几何映射部分排查范围一下子缩小很多。6.3 常规波束形成的稳健边界为什么它仍是基线最后说一个我自己的教训不要因为常规波束形成「简单」就低估它。有一段时间我为了把两个相距 5° 的目标分开直接换了 MVDR仿真里效果惊艳拿到实测数据后波束图被导向矢量误差打得面目全非。回到常规波束形成做基线才把阵列幅相校准的问题暴露出来。自适应算法上限高但下限也低常规算法上限中规中矩下限却非常稳。任何一个阵列信号处理流程都应该先把常规波束形成的结果跑通、和解析解对齐再谈超分辨。希望这个习惯和这份程序能帮到你。本文还有配套的精品资源点击获取