简介：这份资源面向信号处理、无线通信与阵列信号处理方向的学习者和研究者，聚焦复杂噪声环境下的源信号数目估计问题。压缩包内提供AIC、IAIC、MDL、IMDL与MEVARC五种经典信息准则的MATLAB实现，可用于对比不同信噪比条件下各算法的估计准确率，理解模型复杂度与数据拟合度之间的权衡关系。资源共8个文件，以m脚本为主要类型，另有1个asv自动备份文件，整体约5KB，体积轻量便于快速部署与二次修改。目前已有183人学习下载。读者可借助完整代码复现五种算法的性能曲线，通过调整信噪比与数据集观察各准则的适应性差异，并可视化对比结果，为论文实验、课程设计或工程选型提供可运行的参考脚本与排错思路。
1. 源数目估计为什么总在低信噪比下翻车
阵列信号处理里有一类问题看着不起眼，却决定了后续所有环节能不能跑通：源数目估计。你拿到一段多通道采样数据，第一件事不是画 MUSIC 谱，而是先回答"到底有几个源"。这个数字错了，波束形成零陷打偏、MUSIC 谱峰多一个假峰、DOA 跟踪直接发散。标题里的 MDL、MEVARC、信噪比估计，本质上都是围绕这个数字在做判决。
我见过太多人把源数目估计当成"调个阈值"的小事，结果在实验室信噪比 20 dB 时一切正常，一到外场 0 dB 附近就集体翻车。原因很直接：信息论类准则（MDL、AIC）的判决量是特征值函数的差，低信噪比下信号特征值和噪声特征值挤在一起，判决量曲线变得又平又抖，最小值位置随机漂移。这篇笔记就按"这是什么 → 怎么算 → 参数怎么设 → 坑在哪"的顺序，把 MDL 和 MEVARC 这两条主流路线讲透，让你能自己复现并判断值不值得上。
2. MDL 与 MEVARC 的判决逻辑：从特征值到源数目
2.1 协方差矩阵特征分解是共同起点
不管走哪条路，第一步都一样：对阵列接收数据估计协方差矩阵，然后做特征分解。设阵元数 M，快拍数 N，接收数据矩阵 X 是 M×N，协方差估计为 R = (1/N)·X·Xᴴ。对 R 做特征分解得到 M 个特征值 λ₁ ≥ λ₂ ≥ … ≥ λ_M。理论上，前 K 个特征值对应信号子空间，后 M−K 个对应噪声子空间且相等（都等于噪声功率 σ²）。源数目估计的全部难点，就是在这条"台阶"被噪声抹平之后，还能不能找到那个拐点。
import numpy as np def estimate_covariance(X): """ X: 接收数据矩阵, shape=(M, N), M 为阵元数, N 为快拍数 返回: 特征值降序数组, 特征向量矩阵 """ M, N = X.shape # 用样本协方差估计, 注意除以 N 而非 N-1, 与多数文献一致 R = (X @ X.conj().T) / N eigvals, eigvecs = np.linalg.eigh(R) # eigh 针对 Hermitian 矩阵, 更快更稳 # eigh 返回升序, 翻转成降序 idx = np.argsort(eigvals)[::-1] return eigvals[idx], eigvecs[:, idx]逻辑说明：用eigh而不是eig，因为协方差矩阵是 Hermitian 的，eigh保证实数特征值且数值更稳。参数上，除以 N 还是 N−1 在 N 较大时差别可忽略，但 N 小于 50 时建议统一口径，否则和文献对比会有系统性偏差。特征值降序排列是为了后面判决量计算时下标一致。
2.2 MDL 判决量：把信息论准则落到公式
MDL（Minimum Description Length）准则的判决量由两部分组成：似然项加惩罚项。对候选源数目 k（从 0 到 M−1），定义：
- 几何平均与算术平均之比：似然项 = N·(M−k)·ln[ (∏_{i=k+1}^{M} λ_i)^(1/(M−k)) / ( (1/(M−k))·∑_{i=k+1}^{M} λ_i ) ]
- 惩罚项 = 0.5·k·(2M−k)·ln(N)
MDL(k) = 似然项 + 惩罚项，取使 MDL(k) 最小的 k 作为估计。AIC 的区别只在惩罚项是 k·(2M−k)，惩罚更轻，所以 AIC 在高信噪比下容易过估计，MDL 更保守、一致性更好，但低信噪比下容易欠估计。
def mdl_criterion(eigvals, N): """ eigvals: 降序特征值, 长度 M N: 快拍数 返回: (估计源数目, 各候选 k 的 MDL 值数组) """ M = len(eigvals) mdl_vals = [] for k in range(M): noise_eigs = eigvals[k:] # 假设后 M-k 个是噪声特征值 m_k = M - k if m_k <= 0: mdl_vals.append(np.inf) continue # 几何平均 / 算术平均 geo_mean = np.prod(noise_eigs) ** (1.0 / m_k) arith_mean = np.mean(noise_eigs) # 防止除零和 log(0) ratio = geo_mean / max(arith_mean, 1e-12) likelihood = N * m_k * np.log(max(ratio, 1e-12)) penalty = 0.5 * k * (2 * M - k) * np.log(N) mdl_vals.append(likelihood + penalty) mdl_vals = np.array(mdl_vals) k_hat = int(np.argmin(mdl_vals)) return k_hat, mdl_vals逻辑说明：noise_eigs取后 M−k 个特征值，这是 MDL 的核心假设——候选 k 之外的全是噪声。ratio恒小于等于 1，所以likelihood是负数，k 越大噪声子空间越小、比值越接近 1、似然项越接近 0，而惩罚项随 k 增大，两者博弈出最小值。参数上，np.log(N)的底数不影响 argmin，但和文献对比时注意别人用的是 ln 还是 log10。
2.3 MEVARC：不靠信息论，靠特征值方差比
MEVARC（Maximum EigenValue to Arithmetic mean Ratio Criterion 类方法，工程上也常指基于特征值方差比的判决）走的是另一条路：它不构造似然函数，而是直接看相邻特征值的"跳变程度"。常见形式是构造比值序列 r_i = λ_i / λ_{i+1}，或者用每个特征值与剩余特征值均值的比，找最大跳变点。它的优势是不依赖 N 的大样本假设，在快拍数少、信噪比低时往往比 MDL 更稳；代价是需要设一个门限，门限选不好会过估计。
def mevarc_criterion(eigvals, threshold=2.0): """ eigvals: 降序特征值 threshold: 跳变门限, 工程经验值 1.5~3.0 返回: 估计源数目 """ M = len(eigvals) ratios = [] for i in range(M - 1): # 当前特征值与剩余特征值均值的比 rest_mean = np.mean(eigvals[i+1:]) ratios.append(eigvals[i] / max(rest_mean, 1e-12)) ratios = np.array(ratios) # 找第一个超过门限的位置, 其下标+1 即源数目 candidates = np.where(ratios > threshold)[0] if len(candidates) == 0: return 0 return int(candidates[-1] + 1)逻辑说明：ratios[i]衡量第 i 个特征值相对"后面所有特征值平均"的突出程度。信号特征值会明显高于噪声均值，噪声特征值之间比值接近 1。取最后一个超门限的位置而不是第一个，是为了对抗噪声导致的局部抖动。参数threshold是最关键的旋钮，后面避坑章节会专门讲怎么标定。
2.4 信噪比估计为什么绕不开
标题里"信噪比估计"和源数目估计是绑在一起的。原因有两个：一是 MDL 的判决性能对信噪比极度敏感，知道信噪比才能判断当前结果可不可信；二是很多改进算法用信噪比做自适应加权，比如低信噪比时切到 MEVARC、高信噪比时用 MDL。一个实用的信噪比估计用信号子空间和噪声子空间的功率比：
def estimate_snr(eigvals, k_hat): """ 基于特征值分解的信噪比估计 (dB) eigvals: 降序特征值 k_hat: 估计的源数目 """ if k_hat <= 0: return -np.inf signal_power = np.sum(eigvals[:k_hat]) - k_hat * np.mean(eigvals[k_hat:]) noise_power = np.mean(eigvals[k_hat:]) snr_db = 10 * np.log10(max(signal_power / max(noise_power, 1e-12), 1e-12)) return snr_db逻辑说明：信号功率要减掉噪声底，因为信号特征值里也混了噪声功率。k_hat来自前面的估计，所以这是个"鸡生蛋"问题——实践中先用 MDL 粗估，再用估出的信噪比反过来校验。参数上，k_hat取 0 时直接返回负无穷，调用方要处理这个分支。
3. 从仿真到落地：一套可复现的验证流程
3.1 造数据：均匀线阵仿真脚本
要验证算法，先得有可控的数据。下面这段造一个 M 元均匀线阵、K 个窄带远场源、加复高斯白噪声的场景，信噪比可调。
def simulate_ula(M=8, K=3, N=500, snr_db=0.0, seed=42): """ M: 阵元数, K: 真实源数, N: 快拍数, snr_db: 信噪比 返回: 接收数据矩阵 X (M, N), 真实角度 """ rng = np.random.default_rng(seed) # 随机来波方向, 避开端射 angles = np.sort(rng.uniform(-60, 60, K)) # 半波长间距均匀线阵导向矢量 d_lambda = 0.5 A = np.zeros((M, K), dtype=complex) for i, theta in enumerate(angles): phase = 2 * np.pi * d_lambda * np.arange(M) * np.sin(np.deg2rad(theta)) A[:, i] = np.exp(1j * phase) # 源信号, 单位功率复高斯 S = (rng.standard_normal((K, N)) + 1j * rng.standard_normal((K, N))) / np.sqrt(2) # 噪声功率由信噪比反推 noise_power = 10 ** (-snr_db / 10.0) Noise = np.sqrt(noise_power / 2) * ( rng.standard_normal((M, N)) + 1j * rng.standard_normal((M, N))) X = A @ S + Noise return X, angles逻辑说明：导向矢量用半波长间距，sin里是角度转弧度。源信号和噪声都除以 √2 保证复高斯总功率为 1。noise_power由信噪比反推，这样snr_db就是真实的单阵元信噪比。参数seed固定是为了结果可复现，做蒙特卡洛时循环改 seed 即可。
3.2 跑通 MDL 与 MEVARC 的对比
把前面三个函数串起来，做一次完整的估计并对比。
if __name__ == "__main__": for snr in [10, 5, 0, -5]: X, true_angles = simulate_ula(M=8, K=3, N=500, snr_db=snr) eigvals, _ = estimate_covariance(X) k_mdl, _ = mdl_criterion(eigvals, N=500) k_mevarc = mevarc_criterion(eigvals, threshold=2.0) snr_est = estimate_snr(eigvals, k_mdl) print(f"SNR={snr:>3}dB | 真实K=3 | MDL={k_mdl} | MEVARC={k_mevarc} | 估计SNR={snr_est:.1f}dB")典型输出会是：10 dB 时两者都对；0 dB 时 MDL 偶尔欠估计成 2，MEVARC 靠门限还能撑住；−5 dB 时两者都开始不稳。这个对比就是判断"值不值得上"的依据——如果你的场景信噪比长期低于 0 dB，单靠 MDL 不够，必须上 MEVARC 或做多帧累积。
3.3 参数标定：门限和快拍数怎么定
MEVARC 的门限不是拍脑袋。做法是：固定真实源数，扫一遍信噪比和快拍数，画出正确检测概率曲线，选在目标工作点上正确率最高的门限。
| 快拍数 N | 信噪比 | 推荐门限 | 说明 |
|---|---|---|---|
| 200 | 5 dB | 2.5 | 快拍少, 门限抬高防过估计 |
| 500 | 0 dB | 2.0 | 常用工作点 |
| 1000 | −5 dB | 1.6 | 低信噪比, 门限压低防欠估计 |
| 2000 | −5 dB | 1.8 | 快拍多, 统计更稳可略抬 |
这张表是经验区间，不是标准答案。核心逻辑：快拍数少时特征值估计方差大，门限要高；信噪比低时信号特征值被压低，门限要低。两者冲突时优先保信噪比维度，因为低信噪比是主要矛盾。
4. 避坑与排查：源数目估计的五个血泪教训
4.1 现象：高信噪比下 MDL 反而过估计
原因：快拍数 N 远大于阵元数 M 时，惩罚项 0.5·k·(2M−k)·ln(N) 增长过快被误判为"惩罚不够"，实际是似然项在数值上被浮点误差污染。当特征值动态范围超过 1e12，np.prod会下溢成 0，ratio变 0，log出 −inf。
解决：计算几何平均时改用对数域求和，log_geo = np.mean(np.log(noise_eigs))，再likelihood = N * m_k * (log_geo - np.log(arith_mean))。这一步能消掉绝大多数数值翻车。
4.2 现象：MEVARC 在相干源下直接失效
原因：相干源会让协方差矩阵秩亏，两个信号特征值合并成一个，比值序列里少一个跳变，估计值偏小。
解决：先做空间平滑（前向或前后向）恢复秩，再送进 MEVARC。平滑会牺牲有效阵元数，M=8 做 3 阶平滑后有效阵元只剩 6，源数上限跟着降，这个代价要提前算清楚。
4.3 现象：估计结果随快拍数剧烈跳变
原因：快拍数不足时样本协方差估计偏差大，特征值本身就不稳，判决量自然抖。
解决：要么增加快拍，要么做多帧累积——把多段数据的协方差矩阵平均后再分解。累积帧数建议不少于 5，且各帧之间源方位不能变，否则平均会把信号抹掉。
4.4 现象：噪声不是白噪声时全盘失效
原因：MDL 和 MEVARC 都假设噪声特征值相等。实际外场有色噪声，噪声特征值有斜率，判决量最小值位置偏移。
解决：先做白化，用估计的噪声协方差矩阵对数据预白化。没有先验噪声协方差时，可以用对角加载（R + δI，δ 取噪声功率的 1%~5%）压制色噪声影响，代价是低信噪比下性能进一步下降。
4.5 现象：阵元数 M 和源数 K 接近时估计崩溃
原因：可分辨的源数上限是 M−1，当 K 接近 M 时信号子空间和噪声子空间维度接近，任何准则都难分。
解决：这是物理限制，算法救不了。要么加阵元，要么降源数（比如先做波束域降维只保留感兴趣区域）。提前跟需求方说清楚这个上限，比事后解释强。
5. 把估计结果用起来：自适应切换与在线校验
单跑一个准则不够稳，工程上我一般做双准则交叉校验加自适应切换。思路是：同时算 MDL 和 MEVARC，两者一致就直接采信；不一致时看估计信噪比——信噪比高于 5 dB 信 MDL，低于 0 dB 信 MEVARC，中间区间取两者较小值（保守策略，宁可少估一个也不要多估）。
def robust_source_number(eigvals, N, snr_switch_high=5.0, snr_switch_low=0.0): """ 双准则交叉校验 + 信噪比自适应切换 """ k_mdl, _ = mdl_criterion(eigvals, N) k_mevarc = mevarc_criterion(eigvals, threshold=2.0) if k_mdl == k_mevarc: return k_mdl, "consistent" snr_est = estimate_snr(eigvals, max(k_mdl, k_mevarc)) if snr_est >= snr_switch_high: return k_mdl, "high_snr_trust_mdl" elif snr_est <= snr_switch_low: return k_mevarc, "low_snr_trust_mevarc" else: return min(k_mdl, k_mevarc), "conservative"逻辑说明：snr_switch_high和snr_switch_low是两个切换门限，按你的场景标定。中间区间取min是保守策略，因为过估计会让 MUSIC 多出假峰，比欠估计更难排查。返回第二个值是为了记录决策来源，方便事后复盘。
在线校验还有一个技巧：估计出源数目后，用 MUSIC 谱峰个数反查。如果 MUSIC 谱上明显有 4 个峰但估计值是 3，说明欠估计，可以临时把 MEVARC 门限下调 0.2 重算。这个反馈环能把低信噪比下的正确率再拉几个百分点。
最后说个我自己的习惯：任何源数目估计结果，我都会同时打印特征值序列和判决量曲线，而不是只看那个数字。数字会骗人，曲线不会——特征值台阶在哪、判决量最小值平不平，一眼就能看出这次估计可不可信。这个习惯帮我省掉了无数次"结果看着对、实际错"的后悔药。希望帮到你。
本文还有配套的精品资源，点击获取
