简介这份资源面向信号处理、无线通信与阵列信号处理方向的学习者和研究者聚焦复杂噪声环境下的源信号数目估计问题。压缩包内提供AIC、IAIC、MDL、IMDL与MEVARC五种经典信息准则与改进算法的MATLAB实现可用于对比不同信噪比条件下各算法的估计准确率与适应性。资源共8个文件以7个m脚本和1个asv备份文件为主脚本分别对应各算法主程序、复合矩阵构造与统一测试入口压缩包约5KB结构紧凑便于直接运行与二次修改。目前已有183人学习下载。读者可借助这些代码复现源数目估计流程通过调整信噪比与数据矩阵观察五种准则的性能差异并可视化结果为算法选型、论文实验与教学演示提供可复用的脚本基础。1. 源数目估计从一段被噪声淹没的阵列数据说起做阵列信号处理的人迟早会撞上同一个问题接收到的多通道数据里到底有几个信号源波达方向估计、盲源分离、频谱感知这些算法的第一步几乎都是把源数目当成已知量喂进去。可现实中没人会提前告诉你答案你手里只有一段采样矩阵信噪比可能低到 0 dB 以下快拍数可能只有几十个。源数目估计算法就是干这个的——从特征值谱里把信号子空间和噪声子空间切开数出信号有几个。这个方向最经典的落地组合是 MDL 和 MEVARC 两条路线配合信噪比估计做前置判断。标题里那串关键词拆开看就是三件事源数目怎么估、信噪比怎么估、MDL 和 MEVARC 这两个准则怎么用。适合正在做阵列处理、频谱监测、雷达或通信信号检测的工程师尤其是那些已经跑通 MUSIC 但卡在“源数怎么自动定”这一步的人。下面按我实际调过的顺序从原理到代码到踩坑一步步拆开讲。2. MDL 与 MEVARC两个准则的数学底子与选型逻辑2.1 特征值分解之后信息论准则在切什么阵列协方差矩阵做特征分解得到一组从大到小排列的特征值。理想情况下前 K 个特征值对应信号子空间后面 M-K 个相等且等于噪声功率。实际数据里噪声特征值不会完全相等会散开源数目估计的本质就是找一个判决边界把“明显偏大”的特征值数出来。MDL 准则来自信息论里的最小描述长度原理它构造一个代价函数把似然函数和模型复杂度惩罚项加在一起。对数似然部分衡量用假设的源数目去拟合数据有多好惩罚项随源数目增加而增大防止过拟合。取使代价函数最小的那个 K 就是估计结果。它的核心假设是噪声为高斯白噪声快拍数足够多。MEVARC 走的是另一条路它不依赖严格的似然模型而是直接对特征值的离散程度做统计。基本思路是如果第 K 个特征值之后都是噪声特征值它们的方差应该很小一旦某个特征值开始明显偏离噪声特征值的分布就说明信号子空间到此为止。MEVARC 对非高斯噪声和低快拍场景的鲁棒性通常比 MDL 好一些代价是门限需要根据场景调。选型上我一般这样分快拍数大于 200、信噪比高于 5 dB优先用 MDL它稳定且不需要调参快拍数低于 100 或者信噪比在 0 dB 附近晃切到 MEVARC但要把门限系数从默认值往下压。如果连信噪比都不知道先跑一遍信噪比估计再决定用哪个准则。2.2 用 Python 把 MDL 准则跑通的最小实现下面这段代码假设你已经有了阵列接收数据矩阵 X形状是 M 个阵元乘 N 个快拍。先算协方差再特征分解然后按 MDL 公式遍历可能的源数目。import numpy as np def estimate_source_number_mdl(X, max_sourcesNone): X: 阵列数据矩阵, 形状 (M, N), M 为阵元数, N 为快拍数 max_sources: 最大可能源数目, 默认取 M-1 返回: 估计的源数目 K M, N X.shape if max_sources is None: max_sources M - 1 # 计算样本协方差矩阵 R X X.conj().T / N # 特征分解, 取特征值并从大到小排序 eigvals np.linalg.eigvalsh(R) eigvals np.sort(eigvals)[::-1] mdl_values [] for k in range(max_sources 1): # 噪声特征值部分 noise_eigs eigvals[k:] if len(noise_eigs) 0: break # 几何平均与算术平均 geo_mean np.prod(noise_eigs) ** (1.0 / len(noise_eigs)) arith_mean np.mean(noise_eigs) if geo_mean 0: break # 似然项 log_likelihood N * (M - k) * np.log(arith_mean / geo_mean) # 惩罚项 penalty 0.5 * k * (2 * M - k) * np.log(N) mdl_values.append(log_likelihood penalty) if not mdl_values: return 0 # 取使 MDL 最小的 k K int(np.argmin(mdl_values)) return K逻辑说明第 8 行算协方差用的是样本协方差的无偏估计形式除以 N 而不是 N-1这是阵列处理里的惯例。第 11 行用eigvalsh是因为协方差矩阵是 Hermitian 的比通用eig快且数值稳定。第 22 行的似然项里(M-k)是噪声子空间维度log(arith_mean/geo_mean)衡量噪声特征值的离散程度。第 24 行的惩罚项是 MDL 的标准形式系数 0.5 和(2M-k)来自参数自由度的计算。参数说明max_sources一般设成 M-1如果你明确知道源数不会超过 3设成 3 能减少计算量也降低虚警。N是快拍数低于 50 时 MDL 的惩罚项会偏大容易把源数估少这时候建议换 MEVARC。2.3 MEVARC 的实现与门限调整MEVARC 的核心是构造相邻特征值比值的统计量然后和门限比较。下面是一个可用的版本。def estimate_source_number_mevarc(X, threshold1.5): X: 阵列数据矩阵, 形状 (M, N) threshold: 判决门限, 默认 1.5, 低信噪比时降到 1.2 左右 返回: 估计的源数目 K M, N X.shape R X X.conj().T / N eigvals np.linalg.eigvalsh(R) eigvals np.sort(eigvals)[::-1] # 计算相邻特征值的比值 ratios eigvals[:-1] / (eigvals[1:] 1e-12) # 从大到小找第一个超过门限的位置 for i, r in enumerate(ratios): if r threshold: return i return M - 1逻辑说明第 10 行加1e-12是防止除零实际数据里噪声特征值可能非常小。第 12 行遍历比值第一个低于门限的位置就是信号子空间的边界。门限 1.5 是我在快拍数 100、信噪比 5 dB 左右常用的值信噪比每降 5 dB门限往下调 0.1 到 0.2。参数说明threshold是 MEVARC 唯一需要调的参数也是最容易翻车的地方。门限太高会把源数估少太低会把噪声起伏当成信号。建议先用一段已知源数的仿真数据扫一遍门限画出估计正确率随门限变化的曲线取正确率平台段的中点。2.4 信噪比估计怎么给源数目估计兜底信噪比估计在这里的作用是给准则选择提供依据。常用的做法是用特征值分解后的噪声特征值均值作为噪声功率信号特征值减去噪声功率再求和作为信号功率两者比值取对数。def estimate_snr_from_eig(X, K): 根据估计的源数目 K 反推信噪比 X: 阵列数据矩阵 (M, N) K: 估计的源数目 返回: 信噪比估计值 (dB) M, N X.shape R X X.conj().T / N eigvals np.linalg.eigvalsh(R) eigvals np.sort(eigvals)[::-1] noise_power np.mean(eigvals[K:]) if K M else eigvals[-1] signal_power np.sum(eigvals[:K] - noise_power) / K if K 0 else 0 if noise_power 0 or signal_power 0: return -np.inf return 10 * np.log10(signal_power / noise_power)逻辑说明第 11 行用噪声特征值均值估计噪声功率这是最直接的做法。第 12 行把每个信号特征值减去噪声功率再平均得到单源平均信号功率。第 15 行转 dB。参数说明这个估计值在低快拍时会有偏因为噪声特征值本身散开均值不能完全代表真实噪声功率。如果估出来的信噪比低于 0 dBMDL 的结果基本不可信直接切 MEVARC 并把门限降到 1.2 以下。3. 从仿真到实测源数目估计的完整落地流程3.1 仿真数据怎么造才接近真实场景很多人写仿真就造一个干净的高斯信号加高斯白噪声跑出来 MDL 准确率 100%一到实测就崩。问题出在仿真太理想。我一般会加三个东西阵元幅相误差、快拍数限制、色噪声。def generate_array_data(M, N, K, snr_db, amp_err_std0.05, phase_err_std0.05): 生成带阵元误差的阵列接收数据 M: 阵元数, N: 快拍数, K: 真实源数 snr_db: 信噪比 (dB) amp_err_std: 幅度误差标准差 phase_err_std: 相位误差标准差 (弧度) # 随机来波方向 doa np.sort(np.random.uniform(-60, 60, K)) * np.pi / 180 # 均匀线阵导向矢量 array_pos np.arange(M) A np.exp(1j * np.pi * array_pos[:, None] * np.sin(doa)[None, :]) # 加入阵元幅相误差 amp_err 1 amp_err_std * np.random.randn(M, 1) phase_err phase_err_std * np.random.randn(M, 1) A A * amp_err * np.exp(1j * phase_err) # 信号源 S (np.random.randn(K, N) 1j * np.random.randn(K, N)) / np.sqrt(2) # 噪声功率根据信噪比计算 signal_power np.mean(np.abs(A S) ** 2) noise_power signal_power / (10 ** (snr_db / 10)) noise np.sqrt(noise_power / 2) * (np.random.randn(M, N) 1j * np.random.randn(M, N)) X A S noise return X, doa逻辑说明第 12 行到第 15 行是关键阵元幅相误差会让特征值谱的噪声部分不再平坦MDL 的似然项会受影响。第 20 行按实际信号功率反推噪声功率保证信噪比准确。参数说明amp_err_std和phase_err_std设成 0.05 对应 5% 的幅度误差和约 3 度的相位误差这是一般射频前端的水平。如果实测数据来自低成本采集板可以加到 0.1。3.2 用蒙特卡洛扫出准则的适用边界单次实验说明不了问题我习惯跑 500 次蒙特卡洛看不同信噪比和快拍数下的估计正确率。def monte_carlo_test(M, N, K, snr_db, n_trials500): correct_mdl 0 correct_mevarc 0 for _ in range(n_trials): X, _ generate_array_data(M, N, K, snr_db) k_mdl estimate_source_number_mdl(X) k_mevarc estimate_source_number_mevarc(X, threshold1.3) if k_mdl K: correct_mdl 1 if k_mevarc K: correct_mevarc 1 return correct_mdl / n_trials, correct_mevarc / n_trials逻辑说明第 5 行和第 6 行分别调两个准则MEVARC 门限固定 1.3 是为了对比。返回的是正确率。参数说明n_trials至少 500低于 200 时正确率的置信区间太宽。跑完把结果画成正确率随信噪比变化的曲线能直观看到 MDL 在哪个信噪比以下开始掉MEVARC 在哪个区间反超。3.3 实测数据的前处理去均值和通道校准实测数据直接丢进准则函数之前有两步不能省。第一步去均值每个通道减去自己的均值消除直流偏置。第二步通道校准如果各通道增益不一致协方差矩阵的对角线会被拉偏特征值谱跟着变形。def preprocess_array_data(X): X: 原始阵列数据 (M, N) 返回: 预处理后的数据 # 去均值 X X - np.mean(X, axis1, keepdimsTrue) # 通道幅度校准: 用每个通道的功率归一化 channel_power np.mean(np.abs(X) ** 2, axis1, keepdimsTrue) X X / np.sqrt(channel_power 1e-12) return X逻辑说明第 7 行去均值第 10 行用通道功率归一化把各通道增益拉平。注意归一化之后协方差矩阵的对角线变成 1噪声功率的绝对值信息丢了但源数目估计只关心特征值的相对分布所以不影响。参数说明1e-12防止除零。如果某个通道完全没信号归一化会把它放大成噪声这种情况要先检查硬件。4. 避坑与排查源数目估计翻车的五个典型场景4.1 现象MDL 总是把源数估少一个原因快拍数不足时MDL 的惩罚项0.5*k*(2M-k)*log(N)随 k 增长较慢但似然项在低快拍下对信号特征值的估计偏低导致两个源的特征值差异被噪声淹没准则把第二个源归到噪声里。解决把快拍数提到至少 5 倍阵元数。如果硬件限制做不到换 MEVARC 并把门限降到 1.2。还有一个办法是对协方差矩阵做前向-后向平均等效增加快拍数。def forward_backward_avg(X): 前向-后向平均, 等效快拍数翻倍 M, N X.shape J np.fliplr(np.eye(M)) R X X.conj().T / N R_fb (R J R.conj() J) / 2 return R_fb4.2 现象MEVARC 在低信噪比下把噪声起伏当成源原因门限设得太低噪声特征值之间的比值偶尔会超过门限尤其是在快拍数少的时候噪声特征值方差大。解决门限不能低于 1.15。更稳的做法是先对特征值做平滑比如用相邻三个特征值的均值代替单个特征值再算比值。另外可以加一个最小特征值间隔约束只有比值连续两个点超过门限才判定为源。4.3 现象估计结果随快拍数增加反而变差原因信号源之间存在相关性比如多径场景协方差矩阵秩亏信号特征值会合并。快拍数增加不会改善秩亏反而让噪声特征值更集中MDL 的似然项对合并的特征值不敏感。解决对协方差矩阵做空间平滑牺牲一点有效阵元数换取去相关。均匀线阵的话取子阵长度约为 M/2做前后向空间平滑。def spatial_smoothing(X, subarray_len): 空间平滑去相关 M, N X.shape n_sub M - subarray_len 1 R np.zeros((subarray_len, subarray_len), dtypecomplex) for i in range(n_sub): Xi X[i:isubarray_len, :] R Xi Xi.conj().T / N return R / n_sub4.4 现象实测数据跑出来源数忽大忽小原因通道间存在相位不同步或者某个通道接触不良。相位误差会让信号子空间泄漏到噪声子空间特征值谱出现多个伪峰。解决先做通道间互相关检查正常情况下来波方向对应的通道间相位差应该随阵元位置线性变化。如果某个通道和相邻通道的互相关异常低先修硬件再跑算法。软件上可以用对角加载给协方差矩阵对角线加一个小的正值。def diagonal_loading(R, loading_factor0.01): 对角加载, 提高数值稳定性 M R.shape[0] return R loading_factor * np.trace(R) / M * np.eye(M)4.5 现象信噪比估计值比实际高很多原因用估计的源数目 K 去算噪声功率时如果 K 估少了噪声特征值里混进了信号特征值噪声功率被高估信噪比反而被低估如果 K 估多了噪声功率被低估信噪比被高估。这是一个耦合问题。解决不要单次估计就信。用两个准则各跑一遍如果结果一致信噪比估计可信如果不一致取源数目较大的那个去算信噪比然后根据信噪比重新选准则再跑一遍迭代两到三次通常收敛。5. 进阶技巧把源数目估计嵌进 DOA 估计的迭代回路单独跑源数目估计只是第一步真正有价值的是把它和 DOA 估计串起来。我常用的做法是先估源数跑 MUSIC 得到角度谱检查谱峰数量是否和估计的源数一致。如果不一致把源数加一或减一重新跑 MUSIC看哪个源数下谱峰最尖锐、背景最平坦。这个“谱峰一致性”判据比单纯的信息论准则更贴近实际需求。具体操作上定义一个谱峰锐度指标每个谱峰的高度除以峰两侧谷值的均值取所有峰的最小值。源数从 1 遍历到 M-1画锐度随源数的变化曲线取锐度开始下降前的最后一个点。这个方法在低信噪比下比 MDL 稳代价是计算量大因为每个候选源数都要跑一次 MUSIC。还有一个技巧是用估计的源数目反过来修正协方差矩阵。如果估计 K 个源把噪声子空间对应的特征值全部置为它们的均值重构协方差矩阵再跑一次源数目估计。这一步相当于去噪通常能让估计结果稳定下来。迭代两三次就够了再多会过拟合。最后说一个我踩过的坑不要在不同快拍段上分别估源数然后投票。信号源可能只在某些时间段出现分段估计会把间歇性源漏掉。正确做法是用全部快拍估一次如果怀疑源有间歇性先做时间维的检测把有信号的段挑出来再估。这套流程我在几个阵列处理项目里反复用过MDL 和 MEVARC 配合信噪比估计加上前处理和对角加载大部分场景下源数目估计的准确率能到 90% 以上。低信噪比和少快拍是永远的敌人没有银弹只能靠多准则交叉验证和迭代。希望帮到你。本文还有配套的精品资源点击获取