阵列信号处理仿真指南：导向矢量、MUSIC与参数排查

简介：面向阵列信号处理学习者与研究者的MATLAB仿真方法参考文献，以PDF格式收录了重庆大学曾浩等发表于《计算机工程与应用》的期刊论文，重点讲解如何用MATLAB构建阵列信号处理系统模型并完成仿真。包内仅1个PDF文件，压缩包大小232KB，短小精悍但信息密度高，适合需要快速了解阵列信号处理建模仿真框架的读者。目前已有510人学习下载。内容系统覆盖阵列接收信号模型、基于前向平滑的协方差矩阵产生方法、子空间类DOA估计、最小功率波束合成器权值求解以及基本系统参数仿真，并给出关键公式与实现思路；按照文中步骤可进一步扩展完成更复杂的阵列信号处理仿真。对于刚接触DOA估计与波束形成的读者，是一份不错的专业指导资料。

1. 为什么阵列信号仿真先要把信号模型坐实

阵列信号处理的仿真项目里，真正消耗时间的往往不是 MUSIC 或 MVDR 那几行算法代码，而是让“接收数据模型”成立的过程。同一个 8 元均匀线阵，有人仿真得出的谱峰干净利落，有人跑出来全是毛刺——差别通常不在编程水平，而在信号模型与真实采集条件是否匹配。这篇内容以 MATLAB 为工具，围绕阵列信号处理中的模型构建与仿真方法，按最小可复现的路径展开：从导向矢量、协方差矩阵这些地基，到波束形成与 DOA 估计，再到蒙特卡洛验证和参数失效排查。适合雷达、声呐、通信和麦克风阵列方向的工程师，也适合想把仿真结果整理成可汇报结论的研究生。读完你大概能判断：一组乱掉的谱峰，到底该去查参数、查快拍数，还是查模型本身。

2. 均匀线阵信号模型与导向矢量的 MATLAB 构建

2.1 窄带远场假设：为什么要先做这个选择题

阵列信号处理的信号模型，最常用的是窄带远场模型。窄带意味着信号带宽远小于载频，因此信号到达不同阵元时，包络形状基本不变，差别的只是相位。远场意味着信号源离阵列足够远，到达波前可以近似为平面波，阵元间只存在由几何位置差引起的传播延迟。这两个条件共同决定了导向矢量只与来向角度和阵元几何有关，与距离无关。

实际工程中最常见的错误，是在条件不满足时硬套这个模型。比如处理超声脉冲或宽带语音信号，信号带宽达到赫兹量级的相对带宽，直接用一个中心频率的导向矢量建模，得到的协方差矩阵里混入了大量带外分量，谱峰自然不稳定。这种情况下一般先把宽带信号做子带分解，在每个频点分别构建导向矢量，再用聚焦矩阵把不同频点的协方差矩阵对齐到参考频率。判断标准很简单：带宽相对载频小于 1% 时窄带模型基本可靠，超过 10% 就得走宽带路线。

2.2 用 MATLAB 生成阵列流型矩阵的最小代码

均匀线阵（ULA）是最容易上手也是理解阵列信号处理的起点。阵元位置等间距排列在一条直线上，第 m 个阵元相对参考阵元的相位差由 m 倍的路径差决定。下面这段代码生成两个信源的阵列流型矩阵 A：

% 均匀线阵（ULA）参数 M = 8; % 阵元数量 d_lambda = 0.5; % 阵元间距，单位：波长 theta = [-20, 10]; % 两个信源的来向，单位：度 % 构建阵列流型矩阵 A，维度 M x length(theta) idx = (0:M-1).'; % 阵元索引列向量，8x1 A = exp(1j * 2 * pi * d_lambda * idx * sind(theta));

逻辑说明：idx是列向量，sind(theta)是 1×2 的行向量，两者做矩阵乘法得到 8×2 的阵列流型矩阵，每一列对应一个来向角度的导向矢量。导向矢量里每一位元素代表该阵元相对参考阵元的相位延迟，实部虚部共同构成复数基带表示。代码里用sind而不是sin，就是因为sind可以直接接受角度制输入，避免deg2rad转换遗漏。

参数调整：d_lambda是阵元间距对波长的归一化值，0.5是标准配置，对应空间采样率刚好满足奈奎斯特条件。theta的取值决定信源在空间中的位置，仿真时如果想验证分辨率，可以试着把两个角度改成[-3, 3]，你会看到常规波束形成的谱峰完全重叠在一起。

2.3 接收数据、协方差矩阵与快拍数

有了阵列流型矩阵，下一步是把接收数据 X 按照 X = A·S + N 的形式生成，其中 S 是信源复包络，N 是加性噪声。快拍数 N 在这里不只是采样点数，它直接决定协方差矩阵 R 的估计质量：

N = 1000; % 快拍数 snr = 10; % 信噪比，单位：dB sigma = 10^(-snr/20); % 噪声幅度折算 S = exp(1j * 2 * pi * rand(length(theta), N)); % 随机复基带信源 Noise = (randn(M, N) + 1j * randn(M, N)) / sqrt(2); % 复高斯白噪声 % 接收数据模型：X = A * S + sigma * Noise X = A * S + sigma * Noise; % 样本协方差矩阵 R = X * X' / N;

逻辑说明：S的每个元素是单位幅度的随机复包络，相位在 0 到 2π 之间均匀分布，这保证了不同快拍之间、不同信源之间均不相关。Noise的实部和虚部各是标准正态分布，除以 sqrt(2) 使总功率保持为 1。sigma按信号功率为 1 反推，所以信噪比 10 dB 时噪声幅度约 0.316。协方差矩阵R是 M×M 的厄密矩阵，X * X' / N在数学上是对真实协方差矩阵的极大似然估计。

快拍数的选择是模型构建里的隐性参数。理论上只要 N ≥ M 就能让 R 满秩，但实际仿真里 N 太小，协方差矩阵的特征值分布会和真实值偏差很大，导致后面 MUSIC 的噪声子空间估计不准。常见做法是让 N 在 M 的 10 到 20 倍以上，雷达仿真里一个相干处理间隔内拿到几千个快拍很常见。

提示：如果你不想每次都手写上述数据生成代码，MATLAB 的 Phased Array System Toolbox 提供了phased.ULA和phased.Collector等现成对象。但建议至少手写一遍导向矢量构建过程，因为后面的波束形成、MUSIC 谱峰搜索、CRB 计算都要反复用到导向矢量，自己掌握公式才能在报错时知道该查哪里。

3. 波束形成与 MUSIC 空间谱的仿真参数

3.1 CBF 常规波束形成：输出功率空间谱的起点

常规波束形成（CBF）的思路最直接：用一个匹配期望来向的权向量 w 对接收数据做加权合并，然后求输出功率。权向量取 w = a / M，其中 a 是扫描方向对应的导向矢量，M 是阵元数。扫描整个角度范围，把每个方向上的输出功率画出来，就得到空间谱。

theta_scan = -90:0.1:90; % 扫描角度范围 P_cbf = zeros(size(theta_scan)); for k = 1:numel(theta_scan) a_scan = exp(1j * 2 * pi * d_lambda * idx * sind(theta_scan(k))); P_cbf(k) = abs(a_scan' * R * a_scan) / M^2; end % 归一化并转成 dB 显示 P_cbf_dB = 10 * log10(P_cbf / max(P_cbf)); plot(theta_scan, P_cbf_dB); grid on; xlabel('角度 (deg)'); ylabel('归一化功率 (dB)');

逻辑说明：a_scan' * R * a_scan是权向量为a_scan/M时的输出功率，除以M^2是为了让主瓣增益归一化。扫描步长取 0.1 度时，角度分辨率足够覆盖主瓣宽度；如果你只做粗略分析，0.5 度的步长也能用，但谱峰定位精度会直接受影响。

CBF 的优点是稳健，权向量只依赖阵列几何，与接收数据无关，所以低信噪比下也不太会出幺蛾子。缺点是分辨率受瑞利限约束，两个来向间隔小于主瓣宽度时，谱峰无法分辨。比如阵元数为 8、间距半波长时，主瓣宽度约 2/8 弧度，换算到约 14 度——两个相隔 5 度的信号在 CBF 谱里只能看到一个缝都没开的单峰。

3.2 MVDR 自适应波束形成的 3 个必调参数

MVDR（最小方差无失真响应）是 CBF 的改进版：在保证期望方向增益不变的约束下，最小化输出功率。这样做能自适应地在干扰方向形成零陷，但也带来了协方差矩阵求逆的需求。MVDR 的功率谱计算公式从滤波器输出功率推导而来，不需要显式构造权向量，直接写成P = 1 / (a'·R⁻¹·a)的形式。

delta = 1e-3; % 对角加载系数 R_loaded = R + delta * trace(R) / M * eye(M); R_inv = inv(R_loaded); % 求逆 P_mvdr = zeros(size(theta_scan)); for k = 1:numel(theta_scan) a_scan = exp(1j * 2 * pi * d_lambda * idx * sind(theta_scan(k))); P_mvdr(k) = 1 / abs(a_scan' * R_inv * a_scan); end

3 个必调参数按重要程度排序如下：

对角加载是最容易被忽略的参数。R 是由有限快拍估计得到的，小特征值对应的噪声子空间分量很不稳定，直接求逆会把这种不稳定放大，导致谱峰位置乱跳。加载的本质是给矩阵对角线加一个稳定项，相当于给特征值设了下限。delta 太小等于没加，太大则 MVDR 逐渐退化回 CBF——你可以把 delta 改成 1 试试，谱图会明显变钝。

3.3 MUSIC 谱峰搜索：协方差矩阵子空间分解的威力

MUSIC 利用协方差矩阵的特征分解，把特征空间划分为信号子空间和噪声子空间，然后让导向矢量在所有可能方向上与噪声子空间做正交性检验。信号方向上的导向矢量应当正交于噪声子空间，所以1 / (a'·Un·Un'·a)会在信号方向出现尖锐的峰值。

% 特征值分解并按从大到小排序 [U, S_mat] = eig(R); [~, sort_idx] = sort(diag(S_mat), 'descend'); U = U(:, sort_idx); K = size(A, 2); % 信源数，这里为 2 Un = U(:, K+1:end); % 取后 M-K 列为噪声子空间 P_music = zeros(size(theta_scan)); for k = 1:numel(theta_scan) a_scan = exp(1j * 2 * pi * d_lambda * idx * sind(theta_scan(k))); P_music(k) = 1 / abs(a_scan' * Un * Un' * a_scan); end

这段代码里最容易出错的地方是eig的默认排序。MATLAB 的eig返回的特征值矩阵并没有保证有序，必须先排序再取噪声子空间的列索引，否则Un里装的可能全是信号子空间的基。K 的取值直接决定 Un 的构造，如果 K 比真实信源数多，一个真实信号会被漏进噪声子空间，谱峰直接消失；K 少取则会出现虚假峰。信源数估计方法在第 5 章会专门展开。

如果不想做谱峰搜索，可以用 ESPRIT 算法。ESPRIT 利用均匀线阵相邻阵元间的旋转不变性，把角度估计转化为特征值的相位提取，计算量远小于 MUSIC。它的代价是必须要求阵列具备严格的平移不变结构——阵元位置稍有偏差，相位映射到角度时的误差就会被放大。MUSIC 对阵列结构的容错性更好，实测数据里如果阵元位置校准没法保证，优先考虑 MUSIC。

提示：MUSIC 谱峰是无限尖锐的理想极限，实际因为有限快拍和噪声，峰有宽度。要精确定位峰位置，不要直接用扫描网格上的最大值，而是在峰值附近做抛物线插值，至少可以把估计精度从 0.1 度提高到 0.01 度量级。

4. 蒙特卡洛仿真：RMSE 与 CRB 的配合验证

4.1 用蒙特卡洛循环实测算法 RMSE

单次仿真的谱图漂亮不代表算法可靠。随机噪声每次都在变，单次结果可能恰好落在好的一方。工程上一律用蒙特卡洛循环去估计算法的统计性能：固定阵列参数和信噪比，重复生成新数据、重复做估计，最后统计估计值与真值的均方根误差（RMSE）。

rng(42); % 固定随机种子，保证结果可复现 n_trial = 500; % 蒙特卡洛次数 est_angles = zeros(n_trial, 1); theta_true = theta(1); % 取第一个信源作为估计对象 for trial = 1:n_trial % 每次循环重新生成信源和噪声 S_cur = exp(1j * 2 * pi * rand(K, N)); N_cur = (randn(M, N) + 1j * randn(M, N)) / sqrt(2); X_cur = A * S_cur + sigma * N_cur; R_cur = X_cur * X_cur' / N; % 运行 MUSIC，得到谱峰位置 % 这里把 3.3 节的谱计算封装成函数 music_peak() est_angles(trial) = music_peak(R_cur, K, d_lambda, M); end rmse = sqrt(mean((est_angles - theta_true).^2)); fprintf('MUSIC RMSE = %.4f deg\n', rmse);

逻辑说明：循环里每次都要重新生成 S 和 N，这很关键。如果只换噪声不换信源，MUSIC 对同一组信源包络的估计结果会存在系统性偏差，统计出来的 RMSE 不代表真实性能。rng(42)让随机数流固定下来，下次运行同一段代码能得到完全相同的结果，这在调试时特别重要。

蒙特卡洛次数 n_trial 的选择有讲究。500 次是一个起步值，RMSE 的置信区间大约还在 10% 量级；如果文章里的结论要求误差条窄，通常要跑 2000 到 5000 次。代价是计算时间线性上涨，调试阶段先用 100 次，确认算法流程没毛病了再放量跑。

4.2 CRB 克拉美罗界：判断算法还有多少余量

RMSE 反映的是某个算法的实际表现，但没法回答“这个算法距离理论上限还有多远”。克拉美罗界（CRB）给出了任何无偏估计器方差的理论下界，是评估算法性能的基准线。对单信源均匀线阵，CRB 可以由导向矢量对角度的一阶导数解析计算：

theta0 = theta(1); % 目标来向，单位：度 delta_ang = 1e-6; % 数值微分步长 % 导向矢量及其数值导数 a0 = exp(1j * 2 * pi * d_lambda * idx * sind(theta0)); a_p = (exp(1j * 2 * pi * d_lambda * idx * sind(theta0 + delta_ang)) - ... exp(1j * 2 * pi * d_lambda * idx * sind(theta0 - delta_ang))) / (2 * delta_ang); % 正交投影矩阵 P_perp = eye(M) - a0 * inv(a0' * a0) * a0'; % CRB（弧度）并转成角度 crb_rad = sigma^2 / (2 * N * real(a_p' * P_perp * a_p)); crb_deg = rad2deg(sqrt(crb_rad));

代码里a_p是导向矢量对来向角的数值导数，P_perp把导数投影到信号子空间的正交补。crb_deg的物理含义是：在给定信噪比、快拍数和阵列构型下，无偏估计的标准差下界。如果 MUSIC 的 RMSE 已经逼近 CRB，说明算法已没有明显提升空间；如果差了几十倍，优先检查是不是信源数估错或者协方差矩阵构造有问题。

4.3 三种算法的典型性能对比表

在 M=8、d=0.5λ、N=1000、单信源的条件下，三种算法加 CRB 的 RMSE 大致处于以下量级：

这张表的重点是量级关系而不是精确数值，不同随机种子下结果会浮动。CBF 的 RMSE 随信噪比下降得慢，MVDR 在中高信噪比明显优于 CBF，但低信噪比时因为协方差矩阵求逆放大噪声，反而可能不如 CBF。MUSIC 在中高信噪比下性能最好，前提是你已经准确知道了信源数。CRB 作为参考线，能让你一眼看出算法离理论极限还有多远。

5. 阵列信号模型失配与参数失效排查

5.1 阵元间距超过半波长：栅瓣不是算法问题

把第 2 章仿真里的d_lambda从 0.5 改成 0.8，MUSIC 谱里会在真实来向之外多出几个假峰，这些峰称为栅瓣。栅瓣的条件是空间相位差超过 2π，导致多个角度方向的导向矢量产生相位混叠。用公式判断：sin(θ_g) = sin(θ) + m·λ/d，其中 m 是整数。只有当右侧数值落在 [-1, 1] 区间内，才会真的出现栅瓣。

举一个具体例子：d = 0.8λ 时，真实来向 θ = -20°，sin θ = -0.342，m = 1 时 sin θ_g = -0.342 + 1.25 = 0.908，所以 θ_g ≈ 65° 处会出现一个假峰。但如果真实来向是 0°，sin θ = 0，m = 1 时 sin θ_g = 1.25，超出定义域，反而没有栅瓣。排查时先按这个公式算一遍，判断谱峰到底是真的还是几何混叠的产物。分布式阵列经常遇到这个问题，阵元间距按物理空间排布但工作频段变化时，λ 变小导致 d/λ 超过 0.5，栅瓣随之而来。

5.2 相干信号导致协方差矩阵秩亏

两个信号源在物理上完全相干，比如同一个发射信号经过多径到达阵列，S 的第二行和第一行成比例，S 的秩从 2 降到 1。这会导致 R = A·E[SS^H]·A^H + σ²I 中的信号部分秩不足，特征分解后信号子空间少了一个维度，MUSIC 的噪声子空间里混入信号成分，谱峰直接消失或错位。

常见解决手段是空间平滑：把均匀线阵划分成若干重叠子阵，对各个子阵的协方差矩阵取平均，让重排后的矩阵恢复满秩。一段最小实现如下：

L = 4; % 子阵长度，要求 L 大于信源数 sub_M = M - L + 1; % 子阵数量 R_smooth = zeros(sub_M, sub_M); for s = 1:L R_smooth = R_smooth + R(s:s+sub_M-1, s:s+sub_M-1); end R_smooth = R_smooth / L;

逻辑说明：每个子阵对应协方差矩阵的一个对角子块，平均以后相当于人为制造了多个“快照视角”。代价是有效阵元数从 M 降到了 M - L + 1，孔径变小，波束变宽，分辨率下降。L 的选择是权衡：L 越大去相干能力越强，但分辨率损失越明显。前后向平滑可以在同样 L 下提升效果，把R_smooth加上翻转共轭项再取平均即可。

5.3 信源数估计：AIC/MDL 的做法

MUSIC 的性能强依赖信源数 K，但仿真和实际场景里 K 是未知的。一个常规做法是用信息论准则从特征值里自动推断 K：特征值从小到大排序后，前 K 个明显偏大，后面 M-K 个接近噪声功率。AIC 和 MDL 把这个问题建模成模型选择问题，给出可计算的代价函数：

sv = sort(diag(S_mat), 'descend'); % 特征值降序 K_max = M - 1; aic = zeros(K_max, 1); mdl = zeros(K_max, 1); for k = 1:K_max noise_var = mean(sv(k+1:end)); lr = sv(k+1:end) / noise_var; % 特征值与噪声功率之比 % 对数似然项 log_like = -N * (M - k) * log(geomean(lr) / mean(lr)); % AIC / MDL 代价 aic(k) = -2 * log_like + 2 * k * (2 * M - k); mdl(k) = -log_like + 0.5 * k * (2 * M - k) * log(N); end [~, K_aic] = min(aic); [~, K_mdl] = min(mdl);

逻辑说明：geomean(lr) / mean(lr)衡量噪声子空间特征值的分布均匀程度——如果 k 选得合适，剩余特征值应都接近同一个噪声功率，比值接近 1，对数项趋近 0；而 k 选大了之后，模型复杂度惩罚项会迅速增大。AIC 在高信噪比下容易多估信源数，MDL 则倾向于少估，两者结合看：如果 AIC 和 MDL 给出同一个 K，基本可以放心用；不一致时，高信噪比下信 MDL，低信噪比下信 AIC。

6. 模型构建里值得固定下来的 3 个技巧

6.1 用函数句柄封装阵列模型

把“给定角度算导向矢量”这个操作封装成函数句柄，后续 CBF、MVDR、MUSIC、CRB 处处复用，避免每段代码里重复出现那行exp(1j * 2 * pi * d * idx * sind(theta))，也降低改阵元数时漏改某个副本的风险。习惯写法是让句柄只依赖角度参数：

make_steer = @(M, d) @(theta_deg) exp(1j * 2 * pi * d * (0:M-1).' * sind(theta_deg)); steer = make_steer(8, 0.5);

这样写的好处是把阵列构型参数和算法逻辑分开，换 ULA 阵元数或改间距时，一行函数调用即可全局生效。配合matlabFunction还能生成更快的 C 代码版本，但仿真阶段函数句柄的开销完全可以忽略。

6.2 统一随机种子，让仿真结果可复现

蒙特卡洛仿真如果每次运行结果都不一样，排错时很难判断谱峰变化是代码改动引起的还是随机波动引起的。在脚本开头固定三件套：rng(42)、快拍数 N、蒙特卡洛次数，并把 M、d、SNR、theta 打包成一个结构体。改参数时复制这个结构体再修改，留档的仿真结果对应具体的参数快照，汇报时能原样重跑。注意固定种子只对主随机数流有效，如果在代码中间调用过randn或randi，调用顺序变了，结果也会变。

6.3 复基带模型到实测数据的桥接

仿真里从头到尾都在用复基带信号，实测场景中天线直接采样的是实信号，中间隔着下变频和 I/Q 解调。切换之前先在仿真数据上验证一遍完整链路，再引入实际采集数据。实际数据最常见的两个坑：一是阵元位置误差，标称 0.5λ 实际偏差 5% 就会让高信噪比下的 MUSIC 出现系统偏差；二是通道幅相不一致，需要在测向前先做校正，把幅相误差矩阵乘到导向矢量里。仿真里可以在数据生成阶段加入幅相误差项Gamma * A * S，Gamma为对角矩阵，用来模拟通道失配对估计结果的影响，这也是从“仿真好看”走向“实测可用”的必要一步。

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