简介空间平滑MUSIC算法是一种结合子空间类高分辨方法与空间平滑预处理的技术主要用于均匀线阵下相关或相干信号源的方向估计。该MATLAB程序包给出了完整可运行代码包含主程序、空间平滑函数与前后向平滑函数三个源文件压缩包仅3KB轻量便于阅读。实现中详细展示了构造协方差矩阵、特征分解、划分信号与噪声子空间、修正空间平滑矩阵以及谱峰搜索定位角度的流程能够帮助读者厘清MUSIC算法在相关信号场景下失效的原因以及平滑操作如何恢复秩从而提升估计精度。程序附有必要的注释和结果图形便于直接运行观察谱峰位置。目前已有935人学习适合信息与通信工程、雷达信号处理方向的学生在课题仿真或课程设计中借鉴使用。1. 不要让相干信号源毁掉你的MUSIC测向在阵列信号处理中MUSIC算法凭借特征分解把接收数据分成信号子空间与噪声子空间再通过谱峰搜索完成波达方向DOA估计是经典高分辨测向手段之一。但当你把MUSIC从理想独立信源场景搬到真实环境——比如雷达多径反射、5G室内强反射、或通信中的同频干扰叠加——会发现一个尴尬事实当两个入射信号完全相干相关系数为1时信号协方差矩阵的秩发生亏缺MUSIC的谱峰会消失或严重偏移测向彻底失效。空间平滑MUSIC算法正是为这种场景诞生的通过将均匀线阵ULA划分成若干个相互重叠的子阵把各子阵的协方差矩阵取平均用子阵间的平移不变性人为恢复信号协方差矩阵的满秩性从而让MUSIC在相干源下重新工作。本文从秩亏缺的数学本质出发给出空间平滑的推导与MATLAB完整实现并讨论子阵数选择、双向平滑增益、网格搜索步长这些真实工程里绕不开的参数。适合刚接触阵列信号处理、或已经在用MUSIC但被相干源困扰的工程师与研究人员。2. MUSIC为什么怕相干源从秩亏缺说起2.1 信号模型与协方差矩阵的秩设均匀线阵有 (M) 个阵元阵元间距 (d)波长 (\lambda)。假设有 (K) 个远场窄带信号从方向 (\theta_1,\theta_2,\dots,\theta_K) 入射第 (t) 次快拍接收数据为[ \mathbf{x}(t) \mathbf{A}\mathbf{s}(t) \mathbf{n}(t) ]其中 (\mathbf{A} [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \dots, \mathbf{a}(\theta_K)]) 是 (M\times K) 的导向矢量矩阵第 (k) 个导向矢量[ \mathbf{a}(\theta_k) [1, e^{j2\pi d\sin\theta_k/\lambda}, \dots, e^{j2\pi(M-1)d\sin\theta_k/\lambda}]^T ]理想情况下各信号相互独立信号协方差矩阵 (\mathbf{R}_s E[\mathbf{s}(t)\mathbf{s}^H(t)]) 是对角阵秩为 (K)。此时接收数据协方差矩阵[ \mathbf{R} \mathbf{A}\mathbf{R}_s\mathbf{A}^H \sigma^2\mathbf{I} ](\mathbf{R}) 的特征分解中最大的 (K) 个特征值对应信号子空间其余 (M-K) 个小特征值对应噪声子空间。MUSIC通过搜索使导向矢量与噪声子空间正交的角度完成估计谱函数为[ P_{MUSIC}(\theta) \frac{1}{\mathbf{a}^H(\theta)\mathbf{E}_n\mathbf{E}_n^H\mathbf{a}(\theta)} ]当 (\theta) 等于真实入射角时分母趋近于零谱峰出现。2.2 相干源如何让秩“蒸发”当第 (i) 个和第 (j) 个信号完全相干时存在复常数 (\alpha) 使得 (s_j(t) \alpha s_i(t))。此时信号协方差矩阵不再是对角阵例如两个相干源时[ \mathbf{R}_s E\left{\begin{bmatrix} s_1(t) \ \alpha s_1(t) \end{bmatrix} \begin{bmatrix} s_1^(t) \alpha^s_1^(t) \end{bmatrix}\right} \begin{bmatrix} E[|s_1|^2] \alpha^E[|s_1|^2] \ \alpha E[|s_1|^2] |\alpha|^2 E[|s_1|^2] \end{bmatrix} ]这个矩阵的秩只有 1。更一般地若有 (K) 个信号其中 (L) 个相干(\mathbf{R}_s) 的秩降为 (K - L 1)。MUSIC的谱峰数量直接依赖信号子空间的维数秩亏缺导致信号子空间维度不足部分真实方向在谱函数中对应的分母不再趋近于零谱峰丢失。同一个载波频率上若信号之间存在固定相位关系比如多径中直达波与反射波来自同一发射源它们就是相干信号。此时直接运行MUSIC要么漏峰要么把两个相干源识别成单个虚假峰测向完全不可用。2.3 空间平滑的核心思想用平移平均“伪造”秩空间平滑的做法是把 (M) 元均匀线阵划分成 (L) 个相互重叠的子阵每个子阵有 (P M - L 1) 个阵元。第 (l) 个子阵接收数据为[ \mathbf{x}l(t) [x_l(t), x{l1}(t), \dots, x_{lP-1}(t)]^T ]子阵之间的平移量为一个阵元间距。对每个子阵分别计算 (P\times P) 的协方差矩阵 (\mathbf{R}_l E[\mathbf{x}_l(t)\mathbf{x}_l^H(t)])然后取平均[ \mathbf{R}^{SS} \frac{1}{L}\sum_{l1}^{L} \mathbf{R}_l ]关键数学事实当子阵数 (L) 不小于相干源个数时(K) 个相干源在 (\mathbf{R}^{SS}) 中恢复为秩 (K)。背后的原因是每个子阵看到的相干源的相位关系不同平均后等效于给每个源引入了独立的“空间相位扰动”打破了原有的线性相关。需要限定空间平滑只适用于均匀线阵ULA或可等效为ULA拓扑的阵列。L型阵、圆阵要做平滑需要先通过阵列流形变换映射到虚拟均匀线阵这是后话本文不展开。3. MATLAB实现空间平滑MUSIC最小可复现代码3.1 仿真场景设定我们用 MATLAB 模拟一个 8 阵元均匀线阵阵元间距为半波长两个相干信号分别从 (-10^\circ) 和 (20^\circ) 入射。第一个信号是第二个信号的“母本”第二个信号是其幅度衰减、相位偏移的副本。快拍数 500信噪比 15dB。在这个场景下先跑一次原始MUSIC确认失效再跑空间平滑MUSIC验证恢复。3.2 空间平滑函数前向平滑的实现function R_ss forward_smooth(R, L) % 前向空间平滑 % 输入: % R : M x M 接收数据协方差矩阵 % L : 子阵个数 % 输出: % R_ss: 平滑后协方差矩阵 (M-L1) x (M-L1) M size(R, 1); P M - L 1; % 子阵阵元数 R_ss zeros(P, P); for l 1:L idx l : l P - 1; R_ss R_ss R(idx, idx); end R_ss R_ss / L; end这段代码的核心是从原协方差矩阵 (\mathbf{R}) 中取主对角线上连续的 (P\times P) 子块。idx l : l P - 1表示第 (l) 个子阵对应的阵元下标范围例如 (l1) 时取第 1 到第 (P) 个阵元(l2) 时取第 2 到第 (P1) 个阵元。最后除以 (L) 是取平均。这里有一个容易踩的坑如果你直接对快拍数据矩阵 (\mathbf{X})(M\times N)做子阵划分再估计协方差得到的平滑结果与对全局协方差 (\mathbf{R}) 做上述对角子块平均是等价的。因此直接从 (\mathbf{R}) 操作更简单避免重复计算子阵协方差。3.3 双向平滑用共轭对称性把子阵数翻倍前向平滑只利用了阵列一侧的子阵。由于均匀线阵的导向矢量具有共轭对称结构可以构造“后向”子阵前后向结合能在同样阵元数下获得更好的秩恢复能力。function R_fb forward_backward_smooth(R, L) % 前后向空间平滑 M size(R, 1); P M - L 1; % 前向平滑 R_f zeros(P, P); for l 1:L idx l : l P - 1; R_f R_f R(idx, idx); end R_f R_f / L; % 后向平滑: 利用交换矩阵 J 构造反向子阵协方差 J fliplr(eye(M)); R_b conj(J * R * J); R_fb (R_f R_b) / 2; end后向平滑的数学依据是对于均匀线阵接收数据的协方差矩阵满足共轭对称关系 (\mathbf{R}_b \mathbf{J}\mathbf{R}^*\mathbf{J})。把前向平滑结果和后向平滑结果再平均等效于把子阵数量从 (L) 变成 (2L)因此能处理的相干源数量翻倍。当阵元数紧张时双向平滑几乎是必选方案代价是计算量稍增。3.4 完整仿真脚本% 空间平滑MUSIC完整仿真 % 8阵元ULA, 两个相干信号 clear; clc; close all; M 8; % 阵元数 K 2; % 信号数 theta [-10 20]; % 真实方向(度) SNR 15; % 信噪比(dB) N 500; % 快拍数 d_lambda 0.5; % 阵元间距与波长比 % 生成导向矢量矩阵 A zeros(M, K); for k 1:K A(:, k) exp(1j * 2 * pi * d_lambda * (0:M-1) * sind(theta(k))); end % 生成相干信号: s2 是 s1 的 0.8 倍幅度、45度相移 s1 exp(1j * 2 * pi * 0.1 * (1:N)); s2 0.8 * s1 .* exp(1j * pi/4); S [s1; s2]; % 加噪声 X A * S (randn(M, N) 1j*randn(M, N)) / sqrt(2) * 10^(-SNR/20); % 协方差矩阵 R X * X / N; % 原始MUSIC [EigVec, ~] eig(R); [~, idx] sort(diag(EigVec), ascend); En EigVec(:, idx(1:M-K)); % 噪声子空间 theta_scan -90:0.1:90; P_music zeros(1, length(theta_scan)); for i 1:length(theta_scan) a_theta exp(1j * 2 * pi * d_lambda * (0:M-1) * sind(theta_scan(i))); P_music(i) 1 / (a_theta * (En * En) * a_theta); end P_music 10 * log10(P_music / max(P_music)); % 空间平滑MUSIC (前向平滑) L_sub 4; % 子阵数 R_ss forward_smooth(R, L_sub); P_ss M - L_sub 1; % 子阵阵元数 [EigVec_ss, ~] eig(R_ss); [~, idx_ss] sort(diag(EigVec_ss), ascend); En_ss EigVec_ss(:, idx_ss(1:P_ss-K)); P_ss_music zeros(1, length(theta_scan)); for i 1:length(theta_scan) a_ss exp(1j * 2 * pi * d_lambda * (0:P_ss-1) * sind(theta_scan(i))); P_ss_music(i) 1 / (a_ss * (En_ss * En_ss) * a_ss); end P_ss_music 10 * log10(P_ss_music / max(P_ss_music)); % 绘图 figure; plot(theta_scan, P_music, b--, LineWidth, 1.5); hold on; plot(theta_scan, P_ss_music, r-, LineWidth, 1.5); xline(theta(1), k:, LineWidth, 0.8); xline(theta(2), k:, LineWidth, 0.8); xlabel(角度 (deg)); ylabel(归一化空间谱 (dB)); legend(原始MUSIC, 空间平滑MUSIC, 真实方向); grid on; xlim([-60 60]);运行这段代码会看到原始MUSIC谱在 (-10^\circ) 方位完全没有峰值而空间平滑MUSIC在两个真实方向均有清晰谱峰。需要注意sort(diag(EigVec))只对特征值排序如果特征值有重根[~, idx] sort(...)的用法没有问题但更稳妥的做法是用eig同时获取特征向量按特征值大小手动重排。3.5 关键参数说明参数作用选择建议L_sub子阵数量需大于相干源个数本例至少取 2实际取 4~6P_ss子阵阵元数由M - L_sub 1决定(P) 越大分辨率越高但秩恢复能力越差theta_scan谱峰搜索网格网格越细越精确但计算量线性增加普通场景 0.1° 足够L_sub与P_ss是一对矛盾量(P) 大意味着有效孔径大波束更窄、MUSIC分辨率更高(L) 大意味着秩恢复能力强能处理更多相干源。Matlab 里一般从 (L M/2) 开始调观察谱峰是否分裂或偏移再做增减。4. 空间平滑MUSIC实战多径场景与性能验证4.1 构造真实的多径相干场景前面仿真中的相干源是人为构造的。实际工程中常见的相干来源是镜面反射一个基站信号打到建筑物表面反射路径与直达路径同时到达接收阵列反射路径的幅度、相位与直达路径相关但不完全相同且反射路径的DOA通常远离直达路径。这时候直接用MUSIC会漏掉反射方向尤其当反射路径信号更强时这种情况下漏掉的可能反而是直达方向。构造一个三径场景一个直达波从 (5^\circ) 入射两个反射波分别从 (-30^\circ) 和 (40^\circ) 入射三个信号完全相干。阵元数改为 12保证足够的子阵划分空间。M 12; K 3; theta_true [5 -30 40]; L_sub_vals [2 4 6]; figure; for idx_L 1:3 L_sub L_sub_vals(idx_L); R_fb forward_backward_smooth(R, L_sub); P_fb M - L_sub 1; [EigVec_fb, ~] eig(R_fb); [~, idx_fb] sort(diag(EigVec_fb), ascend); En_fb EigVec_fb(:, 1:P_fb-K); for i 1:length(theta_scan) a_fb exp(1j * 2 * pi * d_lambda * (0:P_fb-1) * sind(theta_scan(i))); P_fb_music(i) 1 / abs(a_fb * (En_fb * En_fb) * a_fb); end P_fb_music 10 * log10(P_fb_music / max(P_fb_music)); subplot(1, 3, idx_L); plot(theta_scan, P_fb_music, LineWidth, 1.2); xline(theta_true, k--); grid on; title(sprintf(L_sub %d, L_sub)); xlabel(角度); ylabel(dB); end这段代码演示了子阵数从 2 增加到 6 的过程。运行后会看到L_sub 2时三个峰中有一个位置明显偏移L_sub 4时三峰位置正确但有伪峰L_sub 6时谱峰干净清晰。这里印证了上一节的矛盾关系——在阵元总数充足时优先保证 (L) 大于相干源数量且留有余量。4.2 快拍数与信噪比对平滑效果的影响空间平滑的本质是对多个子阵的协方差做平均因此它对协方差矩阵估计质量的依赖比普通MUSIC更重。将快拍数从 200 降到 50平滑后的协方差矩阵方差显著增大谱峰会出现抖动伪峰概率升高。信噪比低于 5dB 时即使子阵划分正确秩恢复后的小特征值也可能无法与噪声特征值形成明显间隔MUSIC谱峰的分辨力急剧下降。实际工程中如果快拍数受限比如雷达单次扫描只有几十个脉冲可以叠加时间平滑——把多个相邻距离单元的协方差平均后再做空间平滑代码上只需在构造 (\mathbf{R}) 时改成所有距离单元快拍的累加平均R zeros(M, M); for n 1:N_snapshots R R X(:, :, n) * X(:, :, n); end R R / N_snapshots;4.3 如何验证你的空间平滑代码没有写错一个简单的验证方法把两个相干源改为独立源各自的随机相位序列不相关对比原始MUSIC和空间平滑MUSIC的输出两个谱图都应该在真实方向出现谱峰且峰位置误差在 0.2° 以内。如果空间平滑后谱峰分裂大概率是子阵数选得过大导致有效孔径过小而非算法错误。另一个快速检查是计算平滑后协方差矩阵的条件数cond(R_ss)原始协方差矩阵在相干源下条件数通常大于 (10^6)平滑后应降到 (10^3 \sim 10^4) 量级。条件数仍然巨大说明平滑没有真正恢复秩检查子阵数是否大于相干源数。当接收阵列模型存在幅相误差时矩阵不满足理想共轭对称结构这时forward_backward_smooth中的 (\mathbf{R}_b \mathbf{J}\mathbf{R}^*\mathbf{J}) 会引入系统性偏差最终导致谱峰偏移。建议先用前向平滑做基线确认模型误差不显著后再启用双向平滑。5. 提高空间平滑MUSIC精度的两个实用技巧5.1 利用归一化空间谱压制近旁瓣伪峰MUSIC谱函数的尖峰动态范围很大旁瓣在强信号附近可能被抬高到接近主瓣的水平。对谱函数做归一化处理时使用中值滤波后的底噪作为参考而不是直接除以最大值可以有效凸显真实峰。具体做法是先取全角度谱值的median作为噪声底参考再计算显性峰占比P_db 10 * log10(P_ss_music / max(P_ss_music)); noise_floor median(P_db); peak_mask P_db noise_floor 10;然后对peak_mask做findpeaks提取峰位置。这个方法在多径场景中能过滤掉大量不稳定的旁瓣峰比固定阈值更鲁棒。5.2 搜索网格自适应加密的MATLAB实现固定步长 0.1° 在 180° 范围需要 1800 次导向矢量计算。如果先以 1° 步长粗搜定位峰值大致位置再在峰会附近以 0.01° 步长加密搜索计算量可降约一个数量级同时不损失峰值精度。MATLAB实现[~, idx_peak] findpeaks(P_db, MinPeakHeight, noise_floor 10); coarse_angles theta_scan(idx_peak); fine_angles []; for k 1:length(coarse_angles) local_grid (coarse_angles(k)-1):0.01:(coarse_angles(k)1); P_local zeros(size(local_grid)); for i 1:length(local_grid) a_fine exp(1j * 2 * pi * d_lambda * (0:P_ss-1) * sind(local_grid(i))); P_local(i) 1 / (a_fine * (En_ss * En_ss) * a_fine); end [~, idx_max] max(P_local); fine_angles(end1) local_grid(idx_max); endfindpeaks读入的是谱峰矩阵里面包含了主峰和伪峰因此先用噪声底过滤再对每个局部区域做加密搜索。这样得到的fine_angles就是高精度DOA估计值。5.3 双向平滑之后别忘了重新估计信号源数量平滑操作改变了协方差矩阵的维度信号源个数 (K) 不能沿用原始估计。常用的做法是对平滑后的矩阵使用信息论准则AIC或MDLMATLAB中可借用aic函数手动实现evals eig(R_ss); evals_sorted sort(evals, descend); for k 0:P_ss-2 aic_val(k1) -2 * (P_ss-k) * N * log(sum(evals_sorted(k1:end)) / (P_ss-k) / prod(evals_sorted(k1:end))^(1/(P_ss-k))) 2 * k * (2*P_ss-k); end [~, K_est] min(aic_val);注意这里的N是快拍数如果实际中做了时间平均用等效快拍数代入。(K) 估计偏大会导致把噪声子空间的一部分划到信号子空间谱峰数量虚高偏小则会漏掉真实源。在空间平滑场景下(K) 的估计要基于平滑后的 [P\times P] 协方差矩阵进行而不是原 (M\times M) 矩阵。本文还有配套的精品资源点击获取 空间平滑MUSIC算法破解相干源DOA估计秩亏缺及MATLAB实现
简介空间平滑MUSIC算法是一种结合子空间类高分辨方法与空间平滑预处理的技术主要用于均匀线阵下相关或相干信号源的方向估计。该MATLAB程序包给出了完整可运行代码包含主程序、空间平滑函数与前后向平滑函数三个源文件压缩包仅3KB轻量便于阅读。实现中详细展示了构造协方差矩阵、特征分解、划分信号与噪声子空间、修正空间平滑矩阵以及谱峰搜索定位角度的流程能够帮助读者厘清MUSIC算法在相关信号场景下失效的原因以及平滑操作如何恢复秩从而提升估计精度。程序附有必要的注释和结果图形便于直接运行观察谱峰位置。目前已有935人学习适合信息与通信工程、雷达信号处理方向的学生在课题仿真或课程设计中借鉴使用。1. 不要让相干信号源毁掉你的MUSIC测向在阵列信号处理中MUSIC算法凭借特征分解把接收数据分成信号子空间与噪声子空间再通过谱峰搜索完成波达方向DOA估计是经典高分辨测向手段之一。但当你把MUSIC从理想独立信源场景搬到真实环境——比如雷达多径反射、5G室内强反射、或通信中的同频干扰叠加——会发现一个尴尬事实当两个入射信号完全相干相关系数为1时信号协方差矩阵的秩发生亏缺MUSIC的谱峰会消失或严重偏移测向彻底失效。空间平滑MUSIC算法正是为这种场景诞生的通过将均匀线阵ULA划分成若干个相互重叠的子阵把各子阵的协方差矩阵取平均用子阵间的平移不变性人为恢复信号协方差矩阵的满秩性从而让MUSIC在相干源下重新工作。本文从秩亏缺的数学本质出发给出空间平滑的推导与MATLAB完整实现并讨论子阵数选择、双向平滑增益、网格搜索步长这些真实工程里绕不开的参数。适合刚接触阵列信号处理、或已经在用MUSIC但被相干源困扰的工程师与研究人员。2. MUSIC为什么怕相干源从秩亏缺说起2.1 信号模型与协方差矩阵的秩设均匀线阵有 (M) 个阵元阵元间距 (d)波长 (\lambda)。假设有 (K) 个远场窄带信号从方向 (\theta_1,\theta_2,\dots,\theta_K) 入射第 (t) 次快拍接收数据为[ \mathbf{x}(t) \mathbf{A}\mathbf{s}(t) \mathbf{n}(t) ]其中 (\mathbf{A} [\mathbf{a}(\theta_1), \mathbf{a}(\theta_2), \dots, \mathbf{a}(\theta_K)]) 是 (M\times K) 的导向矢量矩阵第 (k) 个导向矢量[ \mathbf{a}(\theta_k) [1, e^{j2\pi d\sin\theta_k/\lambda}, \dots, e^{j2\pi(M-1)d\sin\theta_k/\lambda}]^T ]理想情况下各信号相互独立信号协方差矩阵 (\mathbf{R}_s E[\mathbf{s}(t)\mathbf{s}^H(t)]) 是对角阵秩为 (K)。此时接收数据协方差矩阵[ \mathbf{R} \mathbf{A}\mathbf{R}_s\mathbf{A}^H \sigma^2\mathbf{I} ](\mathbf{R}) 的特征分解中最大的 (K) 个特征值对应信号子空间其余 (M-K) 个小特征值对应噪声子空间。MUSIC通过搜索使导向矢量与噪声子空间正交的角度完成估计谱函数为[ P_{MUSIC}(\theta) \frac{1}{\mathbf{a}^H(\theta)\mathbf{E}_n\mathbf{E}_n^H\mathbf{a}(\theta)} ]当 (\theta) 等于真实入射角时分母趋近于零谱峰出现。2.2 相干源如何让秩“蒸发”当第 (i) 个和第 (j) 个信号完全相干时存在复常数 (\alpha) 使得 (s_j(t) \alpha s_i(t))。此时信号协方差矩阵不再是对角阵例如两个相干源时[ \mathbf{R}_s E\left{\begin{bmatrix} s_1(t) \ \alpha s_1(t) \end{bmatrix} \begin{bmatrix} s_1^(t) \alpha^s_1^(t) \end{bmatrix}\right} \begin{bmatrix} E[|s_1|^2] \alpha^E[|s_1|^2] \ \alpha E[|s_1|^2] |\alpha|^2 E[|s_1|^2] \end{bmatrix} ]这个矩阵的秩只有 1。更一般地若有 (K) 个信号其中 (L) 个相干(\mathbf{R}_s) 的秩降为 (K - L 1)。MUSIC的谱峰数量直接依赖信号子空间的维数秩亏缺导致信号子空间维度不足部分真实方向在谱函数中对应的分母不再趋近于零谱峰丢失。同一个载波频率上若信号之间存在固定相位关系比如多径中直达波与反射波来自同一发射源它们就是相干信号。此时直接运行MUSIC要么漏峰要么把两个相干源识别成单个虚假峰测向完全不可用。2.3 空间平滑的核心思想用平移平均“伪造”秩空间平滑的做法是把 (M) 元均匀线阵划分成 (L) 个相互重叠的子阵每个子阵有 (P M - L 1) 个阵元。第 (l) 个子阵接收数据为[ \mathbf{x}l(t) [x_l(t), x{l1}(t), \dots, x_{lP-1}(t)]^T ]子阵之间的平移量为一个阵元间距。对每个子阵分别计算 (P\times P) 的协方差矩阵 (\mathbf{R}_l E[\mathbf{x}_l(t)\mathbf{x}_l^H(t)])然后取平均[ \mathbf{R}^{SS} \frac{1}{L}\sum_{l1}^{L} \mathbf{R}_l ]关键数学事实当子阵数 (L) 不小于相干源个数时(K) 个相干源在 (\mathbf{R}^{SS}) 中恢复为秩 (K)。背后的原因是每个子阵看到的相干源的相位关系不同平均后等效于给每个源引入了独立的“空间相位扰动”打破了原有的线性相关。需要限定空间平滑只适用于均匀线阵ULA或可等效为ULA拓扑的阵列。L型阵、圆阵要做平滑需要先通过阵列流形变换映射到虚拟均匀线阵这是后话本文不展开。3. MATLAB实现空间平滑MUSIC最小可复现代码3.1 仿真场景设定我们用 MATLAB 模拟一个 8 阵元均匀线阵阵元间距为半波长两个相干信号分别从 (-10^\circ) 和 (20^\circ) 入射。第一个信号是第二个信号的“母本”第二个信号是其幅度衰减、相位偏移的副本。快拍数 500信噪比 15dB。在这个场景下先跑一次原始MUSIC确认失效再跑空间平滑MUSIC验证恢复。3.2 空间平滑函数前向平滑的实现function R_ss forward_smooth(R, L) % 前向空间平滑 % 输入: % R : M x M 接收数据协方差矩阵 % L : 子阵个数 % 输出: % R_ss: 平滑后协方差矩阵 (M-L1) x (M-L1) M size(R, 1); P M - L 1; % 子阵阵元数 R_ss zeros(P, P); for l 1:L idx l : l P - 1; R_ss R_ss R(idx, idx); end R_ss R_ss / L; end这段代码的核心是从原协方差矩阵 (\mathbf{R}) 中取主对角线上连续的 (P\times P) 子块。idx l : l P - 1表示第 (l) 个子阵对应的阵元下标范围例如 (l1) 时取第 1 到第 (P) 个阵元(l2) 时取第 2 到第 (P1) 个阵元。最后除以 (L) 是取平均。这里有一个容易踩的坑如果你直接对快拍数据矩阵 (\mathbf{X})(M\times N)做子阵划分再估计协方差得到的平滑结果与对全局协方差 (\mathbf{R}) 做上述对角子块平均是等价的。因此直接从 (\mathbf{R}) 操作更简单避免重复计算子阵协方差。3.3 双向平滑用共轭对称性把子阵数翻倍前向平滑只利用了阵列一侧的子阵。由于均匀线阵的导向矢量具有共轭对称结构可以构造“后向”子阵前后向结合能在同样阵元数下获得更好的秩恢复能力。function R_fb forward_backward_smooth(R, L) % 前后向空间平滑 M size(R, 1); P M - L 1; % 前向平滑 R_f zeros(P, P); for l 1:L idx l : l P - 1; R_f R_f R(idx, idx); end R_f R_f / L; % 后向平滑: 利用交换矩阵 J 构造反向子阵协方差 J fliplr(eye(M)); R_b conj(J * R * J); R_fb (R_f R_b) / 2; end后向平滑的数学依据是对于均匀线阵接收数据的协方差矩阵满足共轭对称关系 (\mathbf{R}_b \mathbf{J}\mathbf{R}^*\mathbf{J})。把前向平滑结果和后向平滑结果再平均等效于把子阵数量从 (L) 变成 (2L)因此能处理的相干源数量翻倍。当阵元数紧张时双向平滑几乎是必选方案代价是计算量稍增。3.4 完整仿真脚本% 空间平滑MUSIC完整仿真 % 8阵元ULA, 两个相干信号 clear; clc; close all; M 8; % 阵元数 K 2; % 信号数 theta [-10 20]; % 真实方向(度) SNR 15; % 信噪比(dB) N 500; % 快拍数 d_lambda 0.5; % 阵元间距与波长比 % 生成导向矢量矩阵 A zeros(M, K); for k 1:K A(:, k) exp(1j * 2 * pi * d_lambda * (0:M-1) * sind(theta(k))); end % 生成相干信号: s2 是 s1 的 0.8 倍幅度、45度相移 s1 exp(1j * 2 * pi * 0.1 * (1:N)); s2 0.8 * s1 .* exp(1j * pi/4); S [s1; s2]; % 加噪声 X A * S (randn(M, N) 1j*randn(M, N)) / sqrt(2) * 10^(-SNR/20); % 协方差矩阵 R X * X / N; % 原始MUSIC [EigVec, ~] eig(R); [~, idx] sort(diag(EigVec), ascend); En EigVec(:, idx(1:M-K)); % 噪声子空间 theta_scan -90:0.1:90; P_music zeros(1, length(theta_scan)); for i 1:length(theta_scan) a_theta exp(1j * 2 * pi * d_lambda * (0:M-1) * sind(theta_scan(i))); P_music(i) 1 / (a_theta * (En * En) * a_theta); end P_music 10 * log10(P_music / max(P_music)); % 空间平滑MUSIC (前向平滑) L_sub 4; % 子阵数 R_ss forward_smooth(R, L_sub); P_ss M - L_sub 1; % 子阵阵元数 [EigVec_ss, ~] eig(R_ss); [~, idx_ss] sort(diag(EigVec_ss), ascend); En_ss EigVec_ss(:, idx_ss(1:P_ss-K)); P_ss_music zeros(1, length(theta_scan)); for i 1:length(theta_scan) a_ss exp(1j * 2 * pi * d_lambda * (0:P_ss-1) * sind(theta_scan(i))); P_ss_music(i) 1 / (a_ss * (En_ss * En_ss) * a_ss); end P_ss_music 10 * log10(P_ss_music / max(P_ss_music)); % 绘图 figure; plot(theta_scan, P_music, b--, LineWidth, 1.5); hold on; plot(theta_scan, P_ss_music, r-, LineWidth, 1.5); xline(theta(1), k:, LineWidth, 0.8); xline(theta(2), k:, LineWidth, 0.8); xlabel(角度 (deg)); ylabel(归一化空间谱 (dB)); legend(原始MUSIC, 空间平滑MUSIC, 真实方向); grid on; xlim([-60 60]);运行这段代码会看到原始MUSIC谱在 (-10^\circ) 方位完全没有峰值而空间平滑MUSIC在两个真实方向均有清晰谱峰。需要注意sort(diag(EigVec))只对特征值排序如果特征值有重根[~, idx] sort(...)的用法没有问题但更稳妥的做法是用eig同时获取特征向量按特征值大小手动重排。3.5 关键参数说明参数作用选择建议L_sub子阵数量需大于相干源个数本例至少取 2实际取 4~6P_ss子阵阵元数由M - L_sub 1决定(P) 越大分辨率越高但秩恢复能力越差theta_scan谱峰搜索网格网格越细越精确但计算量线性增加普通场景 0.1° 足够L_sub与P_ss是一对矛盾量(P) 大意味着有效孔径大波束更窄、MUSIC分辨率更高(L) 大意味着秩恢复能力强能处理更多相干源。Matlab 里一般从 (L M/2) 开始调观察谱峰是否分裂或偏移再做增减。4. 空间平滑MUSIC实战多径场景与性能验证4.1 构造真实的多径相干场景前面仿真中的相干源是人为构造的。实际工程中常见的相干来源是镜面反射一个基站信号打到建筑物表面反射路径与直达路径同时到达接收阵列反射路径的幅度、相位与直达路径相关但不完全相同且反射路径的DOA通常远离直达路径。这时候直接用MUSIC会漏掉反射方向尤其当反射路径信号更强时这种情况下漏掉的可能反而是直达方向。构造一个三径场景一个直达波从 (5^\circ) 入射两个反射波分别从 (-30^\circ) 和 (40^\circ) 入射三个信号完全相干。阵元数改为 12保证足够的子阵划分空间。M 12; K 3; theta_true [5 -30 40]; L_sub_vals [2 4 6]; figure; for idx_L 1:3 L_sub L_sub_vals(idx_L); R_fb forward_backward_smooth(R, L_sub); P_fb M - L_sub 1; [EigVec_fb, ~] eig(R_fb); [~, idx_fb] sort(diag(EigVec_fb), ascend); En_fb EigVec_fb(:, 1:P_fb-K); for i 1:length(theta_scan) a_fb exp(1j * 2 * pi * d_lambda * (0:P_fb-1) * sind(theta_scan(i))); P_fb_music(i) 1 / abs(a_fb * (En_fb * En_fb) * a_fb); end P_fb_music 10 * log10(P_fb_music / max(P_fb_music)); subplot(1, 3, idx_L); plot(theta_scan, P_fb_music, LineWidth, 1.2); xline(theta_true, k--); grid on; title(sprintf(L_sub %d, L_sub)); xlabel(角度); ylabel(dB); end这段代码演示了子阵数从 2 增加到 6 的过程。运行后会看到L_sub 2时三个峰中有一个位置明显偏移L_sub 4时三峰位置正确但有伪峰L_sub 6时谱峰干净清晰。这里印证了上一节的矛盾关系——在阵元总数充足时优先保证 (L) 大于相干源数量且留有余量。4.2 快拍数与信噪比对平滑效果的影响空间平滑的本质是对多个子阵的协方差做平均因此它对协方差矩阵估计质量的依赖比普通MUSIC更重。将快拍数从 200 降到 50平滑后的协方差矩阵方差显著增大谱峰会出现抖动伪峰概率升高。信噪比低于 5dB 时即使子阵划分正确秩恢复后的小特征值也可能无法与噪声特征值形成明显间隔MUSIC谱峰的分辨力急剧下降。实际工程中如果快拍数受限比如雷达单次扫描只有几十个脉冲可以叠加时间平滑——把多个相邻距离单元的协方差平均后再做空间平滑代码上只需在构造 (\mathbf{R}) 时改成所有距离单元快拍的累加平均R zeros(M, M); for n 1:N_snapshots R R X(:, :, n) * X(:, :, n); end R R / N_snapshots;4.3 如何验证你的空间平滑代码没有写错一个简单的验证方法把两个相干源改为独立源各自的随机相位序列不相关对比原始MUSIC和空间平滑MUSIC的输出两个谱图都应该在真实方向出现谱峰且峰位置误差在 0.2° 以内。如果空间平滑后谱峰分裂大概率是子阵数选得过大导致有效孔径过小而非算法错误。另一个快速检查是计算平滑后协方差矩阵的条件数cond(R_ss)原始协方差矩阵在相干源下条件数通常大于 (10^6)平滑后应降到 (10^3 \sim 10^4) 量级。条件数仍然巨大说明平滑没有真正恢复秩检查子阵数是否大于相干源数。当接收阵列模型存在幅相误差时矩阵不满足理想共轭对称结构这时forward_backward_smooth中的 (\mathbf{R}_b \mathbf{J}\mathbf{R}^*\mathbf{J}) 会引入系统性偏差最终导致谱峰偏移。建议先用前向平滑做基线确认模型误差不显著后再启用双向平滑。5. 提高空间平滑MUSIC精度的两个实用技巧5.1 利用归一化空间谱压制近旁瓣伪峰MUSIC谱函数的尖峰动态范围很大旁瓣在强信号附近可能被抬高到接近主瓣的水平。对谱函数做归一化处理时使用中值滤波后的底噪作为参考而不是直接除以最大值可以有效凸显真实峰。具体做法是先取全角度谱值的median作为噪声底参考再计算显性峰占比P_db 10 * log10(P_ss_music / max(P_ss_music)); noise_floor median(P_db); peak_mask P_db noise_floor 10;然后对peak_mask做findpeaks提取峰位置。这个方法在多径场景中能过滤掉大量不稳定的旁瓣峰比固定阈值更鲁棒。5.2 搜索网格自适应加密的MATLAB实现固定步长 0.1° 在 180° 范围需要 1800 次导向矢量计算。如果先以 1° 步长粗搜定位峰值大致位置再在峰会附近以 0.01° 步长加密搜索计算量可降约一个数量级同时不损失峰值精度。MATLAB实现[~, idx_peak] findpeaks(P_db, MinPeakHeight, noise_floor 10); coarse_angles theta_scan(idx_peak); fine_angles []; for k 1:length(coarse_angles) local_grid (coarse_angles(k)-1):0.01:(coarse_angles(k)1); P_local zeros(size(local_grid)); for i 1:length(local_grid) a_fine exp(1j * 2 * pi * d_lambda * (0:P_ss-1) * sind(local_grid(i))); P_local(i) 1 / (a_fine * (En_ss * En_ss) * a_fine); end [~, idx_max] max(P_local); fine_angles(end1) local_grid(idx_max); endfindpeaks读入的是谱峰矩阵里面包含了主峰和伪峰因此先用噪声底过滤再对每个局部区域做加密搜索。这样得到的fine_angles就是高精度DOA估计值。5.3 双向平滑之后别忘了重新估计信号源数量平滑操作改变了协方差矩阵的维度信号源个数 (K) 不能沿用原始估计。常用的做法是对平滑后的矩阵使用信息论准则AIC或MDLMATLAB中可借用aic函数手动实现evals eig(R_ss); evals_sorted sort(evals, descend); for k 0:P_ss-2 aic_val(k1) -2 * (P_ss-k) * N * log(sum(evals_sorted(k1:end)) / (P_ss-k) / prod(evals_sorted(k1:end))^(1/(P_ss-k))) 2 * k * (2*P_ss-k); end [~, K_est] min(aic_val);注意这里的N是快拍数如果实际中做了时间平均用等效快拍数代入。(K) 估计偏大会导致把噪声子空间的一部分划到信号子空间谱峰数量虚高偏小则会漏掉真实源。在空间平滑场景下(K) 的估计要基于平滑后的 [P\times P] 协方差矩阵进行而不是原 (M\times M) 矩阵。本文还有配套的精品资源点击获取