MDL与MEVARC结合的源数目估计：低信噪比下的稳健实现

简介：面向信号处理、阵列信号处理及无线通信领域的MATLAB源数目估计代码包，专门解决在噪声背景下推断观测数据中信源个数的问题。压缩包内共有八个文件，包括七个m格式的源码脚本与一个asv备份文件，整体大小仅五KB，涵盖AIC、IAIC、MDL、IMDL与MEVARC五种经典信息准则算法，并附带主程序与复合矩阵生成脚本，可直接运行对比不同准则的估计输出。通过这些代码，用户能够改变信噪比参数来评估算法在各种噪声环境下的性能，从而获得源数目估计准确率随信噪比变化的曲线，深入理解模型选择原理及各算法的适用条件。包体小巧但功能集中，目前已有183人学习下载，适合信号处理课程设计、科研验证与算法入门参考，可为相关课题提供便捷的算法对比平台。

1. 源数目估计：为什么说 MDL 的低信噪比翻车不是玄学

做阵列信号处理的人都知道，DOA 估计里最容易被低估的一步是源数目估计——跑 MUSIC 或 ESPRIT 之前，得先告诉算法有几个信号源。MDL 是我最早用的准则，可它在 0dB 以下的表现经常让人想摔键盘：明明 3 个源，结果直接给了 1 个。问题往往不在 MDL 本身，而在于它背后的噪声方差参数在实测里根本拿不准。把 MDL 和 MEVARC 这类特征值迭代型信噪比估计绑在一起，先让噪声功率稳定下来，MDL 的判决才站得住。下面把这条链路拆开讲透：从特征值谱到惩罚项，再到能直接复现的 MATLAB 代码、参数表和五个踩坑记录。

2. MDL 准则与 MEVARC：特征值谱、惩罚项和噪声方差估计到底在做什么

2.1 信号模型与特征值谱：K 个源为什么对应 K 个大特征值

假设 M 元均匀线阵，阵元间距为半波长。接收矩阵 X 是 M×N，每一列是一次快拍。如果有 K 个不相关的窄带信号从不同方向到达，接收数据可以写成一个经典的叠加模型：X = A·S + N。A 是由 K 个方向导向矢量拼成的方向矩阵，S 是信号的基带复包络，N 是噪声。源数目估计的本质，就是把 S 和 N 的贡献在统计上分开。

对协方差矩阵 R = E{XX^H} 做特征分解后，理想情况下的特征值谱形态非常清晰：前 K 个特征值对应信号，在信噪比大于 0dB 时明显大于噪声特征值；后 M-K 个特征值全部等于噪声方差 σ²。也就是说，特征值谱上存在一个“台阶”，台阶的位置就是 K。MDL 做的事情，本质上就是用一个统计准则去定位这个台阶，而不是靠肉眼去数。

但实际 R 只能由有限快拍估计出来，特征值谱不会那么干净。信号特征值会被噪声拉平，噪声特征值也会在 σ² 附近上下抖动。当信噪比很低或者快拍数不够时，谱上最末端的“大”特征值到底是弱信号，还是强噪声的抖动，从数据本身很难区分清楚。这就是源数目估计（source-number-estimation）最麻烦的地方，也是后面 MDL 和 MEVARC 合作的空间。

2.2 MDL 准则：从似然比到惩罚项，为什么比 AIC 更稳

MDL 的全称是最小描述长度，最早是信息论里的模型选择思想。把它用在源数目估计上时，思路很直接：对每个候选源数 k，先算一个似然项，描述“把后 M-k 个特征值当成噪声”这个假设和数据有多吻合；再加上一个惩罚项，防止 k 一味往大里取。最后取使总代价最小的 k 作为估计结果。

具体实现上，先把特征值按从大到小排序，记为 λ1 ≥ λ2 ≥ … ≥ λM。对候选源数 k，剩下的 M-k 个特征值是 λk+1 到 λM。分别计算它们的几何均值 g(k) 和算术均值 a(k)，定义似然项：

L(k) = N·(M-k)·log( a(k) / g(k) )

如果后 M-k 个特征值分布集中，g 和 a 接近，L(k) 接近 0；如果里面混着信号特征值，分布散，a 明显大于 g，L(k) 就会偏大。加上惩罚项后，MDL 的判据是：

MDL(k) = L(k) + 0.5·k·(2M-k)·log(N)

对 k 从 0 到 M-1 逐个算，取使 MDL 最小的那个 k。AIC 和 MDL 长得几乎一样，只是把惩罚项里的 0.5·log(N) 换成了常数 2，即 AIC(k) = L(k) + k·(2M-k)。

从渐近性看，N 越大，MDL 的惩罚项最终会超过 AIC，所以 MDL 不容易多估源数，而 AIC 在高信噪比大快拍下倾向于把源数估计得多一点。反过来，在低信噪比和小快拍场景，MDL 又因为惩罚项过大而偏向低估。所以“MDL 比 AIC 稳”是有前提的，它说的是渐近一致性，不是所有信噪比下都好用。

2.3 MEVARC：一类用迭代特征值统计量逼近噪声功率的估计器

MEVARC 这个名字在不同代码库里出现过多次，细节不完全一致，核心思路是一致的：从特征值谱里估计出稳定的噪声方差，作为 MDL 判决时的参考。我习惯把它看成一种“带截断和迭代的噪声功率平均”。为什么用截断平均而不是直接把尾部特征值做平均？因为低信噪比时尾部特征值里可能混入弱信号，或者有一些特别小的离群值，直接平均对这些离群值太敏感。

一个典型实现如下：

function sig2 = mevarc(lam, K0, iters) % lam: 降序排列的特征值列向量 % K0: 当前估计出的源数 % iters: 迭代次数 lam = sort(lam(:), 'descend'); M = numel(lam); sig2 = mean(lam(K0+1:end)); for it = 1:iters mask = (lam >= 0.5*sig2) & (lam <= 1.5*sig2); if ~any(mask) break; end sig2_new = mean(lam(mask)); if abs(sig2_new - sig2) / sig2 < 1e-3 sig2 = sig2_new; break; end sig2 = sig2_new; end end

代码里最关键的是 0.5 到 1.5 倍的截断窗口。第一次循环用全部尾部特征值算一个初值，然后只保留落在窗口内的特征值做下一次平均，相当于把特别大和特别小的离群值都剔掉。窗口参数是经验值，后面第 4 章会专门讲怎么调。

得到 σ² 之后，信噪比估计就顺理成章了。特征值总和与总功率成正比，所以可以估算阵列端的平均信噪比：

SNR_hat = ( mean(lam) - sig2 ) / sig2

这里 mean(lam) 是所有特征值的平均，也就是总功率的估计；减去噪声功率 σ² 再除以 σ²，就是信号分量和噪声分量的功率比。换算成 dB 后，可以作为选择 MDL 还是 AIC 的开关条件。MEVARC 的价值就在这里：它不只给 MDL 提供一个更稳的噪声功率，还给整个判决链提供了一个可信的信噪比估计值。

3. 用 MATLAB 把 MDL+MEVARC 跑通：最小实现、参数表与判读方法

3.1 造一组合适的阵列数据：阵元数、快拍数与信噪比的控制

先说数据怎么来。实测数据当然最好，但在复现代码前，先用仿真数据把逻辑跑通。下面这段代码生成一个 M 元均匀线阵的接收矩阵：

M = 8; % 阵元数 K_true = 3; % 真实源数 N = 200; % 快拍数 SNR_dB = -3; % 单源信噪比 snr = 10^(SNR_dB/10); theta = [-20 10 35] * pi/180; % 三个来波方向 A = exp(1j*pi*(0:M-1)' * sin(theta)); % 导向矢量矩阵 S = sqrt(snr) * randn(K_true, N); % 信号复包络 X = A * S + randn(M, N); % 叠加单位方差噪声 R = (X * X') / N; % 样本协方差矩阵

参数方面需要重点说明四个。M 是阵元数，直接决定特征值谱的长度，M 太小比如只有 4 个，最多只能区分 3 到 4 个源。K_true 是真实源数，只用于事后评判，实际处理时是未知的。N 是快拍数，也就是处理用到的采样点数，N 越大，协方差矩阵 R 越接近统计意义上的真实 R。SNR_dB 是每个信号源的单源信噪比，定义是信号功率与噪声功率之比，代码里噪声功率归一化为 1。

这里每个源的功率设成一样，是为了先看清楚 MDL 在理想均匀情况下的行为。实际信号强弱往往差别很大，这一点在避坑章节会专门展开。快拍数 N 和阵元数 M 的关系也要注意，一般至少要有 N 大于 2M，否则 R 的秩不够，特征值谱分布会很糟糕。我在实测数据上只处理 N 大于等于 10M 的情况，低于这个值通常要加平滑或对角加载。

3.2 经典 MDL 实现：从公式到代码的逐行映射

拿到 R 后，第一步是对 R 做特征分解，然后把 MDL 和 AIC 的代价全部算出来。完整代码如下：

lam = sort(eig(R), 'descend'); mdl = zeros(1, M); aic = zeros(1, M); for k = 0 : M-1 lam_n = lam(k+1:end); % 取后 M-k 个特征值作为噪声候选 r = M - k; g = geomean(lam_n); a = mean(lam_n); like = N * r * log(a / g); mdl(k+1) = like + 0.5 * k * (2*M - k) * log(N); aic(k+1) = like + k * (2*M - k); end [~, idxM] = min(mdl); Khat_mdl = idxM - 1; [~, idxA] = min(aic); Khat_aic = idxA - 1; fprintf('MDL估计源数: %d, AIC估计源数: %d\n', Khat_mdl, Khat_aic);

这段代码有一个地方特别容易搞错：候选 k 的取值范围。某些实现会从 1 循环到 M，但那样 k=M 时噪声特征值集合为空，没有意义。所以这里从 0 开始，最多到 M-1。对应到 MATLAB 索引上，k=0 时取全部特征值，k=M-1 时只取最后一个特征值。因为 MATLAB 数组下标从 1 开始，所以数组里存的是 k+1 位置。

似然项 like 等于 N·r·log(a/g)。当候选 k 等于真实源数时，后面的特征值都来自噪声，它们围绕 σ² 分布，a 和 g 接近，like 接近 0。如果 k 小于真实源数，尾部特征值里混进了信号，a 明显大于 g，like 变大。惩罚项则随 k 增大而增大，两者相抵后的最小点就是判决结果。

注意：log(a/g) 在 a 等于 g 时为 0，浮点误差可能让它变成微小的负数，后续如果要做阈值判断，建议先对 like 做 max(like, 0) 处理。

跑完这段，大概率会看到在 -3dB 这个信噪比下 MDL 给出的是 2 或者 1。这不是代码错了，而是经典 MDL 在低信噪比下的固有毛病，下面直接上 MEVARC。

3.3 用 MEVARC 替换解析噪声方差：迭代判决与信噪比估计的联动

MEVARC 不是替代 MDL，而是给 MDL 提供一个可靠的噪声方差参考，再用它估计信噪比来决定要不要切换到 AIC。完整用法如下：

sig2 = mevarc(lam, Khat_mdl, 5); snr_hat = (mean(lam) - sig2) / max(sig2, 1e-12); snr_hat_dB = 10 * log10(max(snr_hat, 1e-6)); if snr_hat_dB < 0 Khat = Khat_aic; else Khat = Khat_mdl; end fprintf('MEVARC噪声方差: %.3f, 估计信噪比: %.1f dB, 最终源数: %d\n', ... sig2, snr_hat_dB, Khat);

这里的逻辑是，先拿经典 MDL 的结果作为 MEVARC 的初值，再用 MEVARC 得到更稳的噪声方差，进而估计信噪比。如果估计出的信噪比低于 0dB，说明当前处于 MDL 最容易低估的区间，改用 AIC 的判决结果；如果信噪比高于 0dB，MDL 渐近一致性的优势能发挥出来，继续用 MDL。这个切换阈值不是拍脑袋定的，我在不同阵元数和快拍数下试过，0dB 附近是一条比较合理的分界线。

参数方面，MEVARC 的迭代次数一般取 3 到 5 次就够。第 1 次迭代主要用来剔除明显离群的噪声特征值，第 2 到第 3 次收敛，超过 5 次基本没有变化，反而可能因为把真实信号特征值也排除在窗口外而产生偏差。窗口上下界 0.5 和 1.5 是输入参数，如果已知接收机底噪非常稳定，可以收窄到 0.7 到 1.3；如果现场干扰比较大，放宽到 0.3 到 2.0 更稳。

跑通这套流程后，手上就有了一个“MDL+MEVARC”的源数目估计器。下一步不是直接上实测，而是用蒙特卡洛把正确率摸清楚，否则你不知道它在什么条件下会失效。

4. 正确率怎么验证：蒙特卡洛脚本与三个必调参数

4.1 蒙特卡洛循环：500 次试验画出准确率曲线

源数目估计是一个随机过程，单次跑得对说明不了问题。我习惯的做法是固定一组参数，重复几百次随机试验，统计 Khat 等于 K_true 的比例，作为该条件下的正确率。下面这个脚本把 MDL、AIC、MDL+MEVARC 三种方案放在同一个循环里比较：

rng(2024); trials = 500; snr_list = -10:5:10; acc_mdl = zeros(size(snr_list)); acc_aic = zeros(size(snr_list)); acc_hyb = zeros(size(snr_list)); for si = 1:numel(snr_list) cnt_mdl = 0; cnt_aic = 0; cnt_hyb = 0; for t = 1:trials % 生成数据 snr = 10^(snr_list(si)/10); A = exp(1j*pi*(0:M-1)' * sin(theta)); S = sqrt(snr) * randn(K_true, N); X = A * S + randn(M, N); R = (X * X') / N; % 特征值分解 lam = sort(eig(R), 'descend'); % 封装函数：返回 MDL 和 AIC 的估计源数 [Khat_mdl, Khat_aic] = mdl_aic_from_lam(lam, N); % 混合方案 sig2 = mevarc(lam, Khat_mdl, 5); snr_hat_dB = 10*log10((mean(lam)-sig2)/sig2); if snr_hat_dB < 0 Khat_hyb = Khat_aic; else Khat_hyb = Khat_mdl; end cnt_mdl = cnt_mdl + (Khat_mdl == K_true); cnt_aic = cnt_aic + (Khat_aic == K_true); cnt_hyb = cnt_hyb + (Khat_hyb == K_true); end acc_mdl(si) = cnt_mdl / trials; acc_aic(si) = cnt_aic / trials; acc_hyb(si) = cnt_hyb / trials; end

判定标准是严格相等，估计值差一个源都算错。这在工程上有点苛刻，因为当源数较多、信噪比不平衡时，漏掉一个弱源在 DOA 后处理里往往还有补救空间，但作为方法对比，严格相等最公平。

我跑过的典型结果大致如下，数值会随随机种子和方向设置浮动，看趋势即可：

趋势是稳定的：信噪比越低，MEVARC 介入带来的提升越明显。到 5dB 以上三者基本没有差别，这时候拼的就是谁在高信噪比下更不容易多估。

4.2 快拍数与平滑窗口：两个直接影响 MDL 硬指标的参数

快拍数 N 是除信噪比之外影响最大的参数。N 太小，协方差矩阵的统计波动太强，特征值谱完全变形，MDL 的似然项对候选 k 的区分度下降。工程上的经验线是 N 至少大于 2M，一般建议 10M 以上。如果数据长度不够，有三条路可以走。

第一个是对角加载，给协方差矩阵对角线加一个很小的量，抑制噪声特征值过分离散。第二个是时间平滑，把连续多个快拍的数据做滑动平均后再算协方差。第三个是空间平滑，用子阵滑动来恢复协方差矩阵的秩，这个对相干源场景特别有效。

常见做法是先看快拍数够不够，不够就加对角加载。加载量取 R 对角线均值的千分之一到百分之一，加得太多会把小信号特征值抹平，导致低估源数。空间平滑放在相干场景里再讲，那里才是它的主场。

4.3 三个必调参数：MEVARC 窗口、迭代次数与判决切换阈值

把参数说透是复现的关键。MEVARC 窗口就是截断平均的上下边界，默认 0.5 到 1.5。窗口收窄到 0.7 到 1.3，能更严格地剔除噪声离群值，但风险是真实的弱信号特征值也被剔除，导致噪声方差偏低、信噪比估计偏高；窗口放宽到 0.3 到 2.0，则更保守，适合现场干扰明显、特征值谱重尾很长的场景。

迭代次数一般设 3 到 5 次。你可以打印每次迭代的 sig2 变化，如果第二次和第三次的相对变化小于千分之一，就说明已经收敛。如果始终不收敛，往往是初值 K0 给得太离谱，需要回到 MDL/AIC 的结果检查。

判决切换阈值默认 0dB，但也要看任务。如果你宁可多估一个源也不愿意漏源，比如后面还有 MUSIC 谱峰复核步骤兜底，可以把阈值提高到 3dB 或 5dB，让 AIC 介入更频繁；反过来，如果对虚警敏感，阈值可以下调到 -3dB。

提示：参数微调请固定其他变量，一次只动一个。只靠感觉同时调三个，你根本分不清性能变化是谁引起的。我自己吃过这个亏，后来参数微调一律用控制变量法。

5. 源数目估计避坑手册：五个实测翻车场景与对应修复手段

5.1 低信噪比低估源数：MDL 把 3 个源估成 1 个

现象：真实源数 3，信噪比 -5dB，MDL 输出 1 或 2，多次试验结果稳定偏低。

原因：低信噪比时弱信号特征值与噪声特征值在幅度上无法区分，MDL 的似然项 L(k) 对“少一个源”和“多一个噪声特征值”的敏感度下降。与此同时惩罚项仍然随 k 增长，于是模型倾向于选择更小的 k。

解决：第一步用 MEVARC 先估计噪声方差和信噪比；如果信噪比估计低于 0dB，直接切换到 AIC 结果，并保留 MEVARC 估计的信噪比作为可信度标记。低信噪比下的信噪比估计是整个源数目估计成败的分水岭，越早介入越好。

5.2 特征值弥散严重：MEVARC 算出的噪声功率比接收机底噪高一个量级

现象：MEVARC 输出的 sig2 明显大于接收机标定的底噪功率，MDL 随之把源数估大。

原因：有限快拍下，噪声特征值本身在 σ² 附近随机分布，其中最大的几个噪声特征值可能比 σ² 大出数倍。直接做算术平均时，这几个偏大的值把均值拉高了，导致噪声功率高估。

解决：把截断窗口从 0.5~1.5 收窄到 0.5~1.2，或者改用中位数代替均值，只保留特征值分布的主体。MEVARC 的价值就在于这个截断，而不是单纯的平均。中位数对重尾更鲁棒，但在小快拍下会偏低，需要结合迭代次数来平衡。

5.3 通道幅相不一致引入伪源

现象：实测 8 通道接收机存在 1 到 2dB 的增益差异，MDL 在无源测试时输出 1 个“源”，但暗室其实是空的。

原因：通道增益不一致导致噪声功率在通道间不同，协方差矩阵不再满足 σ²I 的结构，特征值谱出现接近小信号特征值的凸起，MDL 把硬件失配当成信号。

解决：在跑 MDL 之前先做通道校准，要么用校准数据估计增益并归一化，要么对协方差矩阵做对角归一化，把对角元素都标准化为 1。校准后特征值谱的噪声部分会平很多。这个坑最容易在从仿真转实测时出现：仿真里永远不会遇到通道失配，实测里却是家常便饭。

5.4 快拍太少：MDL 直接估出 0 个源

现象：N=50，M=8，真实源数 3，MDL 输出 0，AIC 输出 4。

原因：快拍太少时，协方差矩阵的秩不够，特征值谱被统计波动抹平，信号与噪声特征值混在一起，a/g 对所有 k 都接近 1，似然项失去区分度。此时惩罚项主导判决，k=0 的惩罚项为 0，自然成为最优选择。问题不单纯是惩罚项过大，而是似然项已经分不出信号和噪声。

解决：先检查数据长度能不能延长；不能延长就用空间平滑或前向-后向平均增加有效快拍数；如果还是不行，至少把 AIC 的结果一起打印出来，综合两者判断。快拍太少时任何准则的置信度都不高，必须同时输出信噪比估计值辅助判断。

5.5 相干源场景：MDL 输出比实际来波数少

现象：两个来波是同一个信号经多径到达，MUSIC 谱上能看到两个峰，但 MDL 只输出 1。

原因：相干信号使得信号协方差矩阵秩亏，原本应该突出的第二个信号特征值被合并进噪声区，MDL 在统计上无法把它们区分开。

解决：先做前向-后向空间平滑，再对平滑后的协方差矩阵做 MDL。核心代码如下：

L = M - 2; % 子阵长度，经验值取 M/2 到 2M/3 Rf = zeros(L, L); for i = 1 : M-L+1 Rf = Rf + R(i:i+L-1, i:i+L-1); end % 反向平滑：用反对角交换矩阵 J 对 R 做共轭翻转 J = fliplr(eye(M)); Rb = J * conj(R) * J; Rf2 = zeros(L, L); for i = 1 : M-L+1 Rf2 = Rf2 + Rb(i:i+L-1, i:i+L-1); end Rsm = (Rf + Rf2) / (2 * (M - L + 1)); % 对 Rsm 重新做特征值分解和 MDL 判决

子阵长度 L 决定了空间平滑能恢复的秩。L 越大，平滑用的子阵数越少；L 太小，阵元损失太多。经验上 L 取 M/2 到 2M/3 之间比较合适。平滑之后源数目估计的物理含义也变了：此时估出来的是空间上可分辨的相干组数，而不是严格意义上的来波数，后端的 DOA 算法需要配合起来理解。

6. 把 MDL 结果当先验：用 MUSIC 谱峰梯度做二次验证的下班技巧

6.1 用 MDL 定子空间维度，再看谱峰有没有“多余”的峰

MDL 和 AIC 说到底是在统计意义上做模型选择，它们不感知空间谱的形状。MUSIC 谱的峰是几何意义上的来波方向，两种信息可以互相印证。我现在的习惯是不让 MDL 当唯一裁判，而是把它当成一个粗糙先验，再用 MUSIC 谱的峰结构做二次验证。

做法是先得到 Khat，然后计算 MUSIC 空间谱，取幅度最大的前 Khat+2 个峰，比较第 Khat 个峰和第 Khat+1 个峰之间的 dB 差：

music_spectrum_dB = 10*log10(music_spectrum + 1e-12); [peaks, locs] = findpeaks(music_spectrum_dB, 'SortStr', 'descend'); if numel(peaks) < Khat + 1 warning('谱峰数量不足：Khat=%d，实际峰数=%d', Khat, numel(peaks)); else gap = peaks(Khat) - peaks(Khat+1); fprintf('Khat=%d, 第%d峰与第%d峰差 %.1f dB\n', Khat, Khat, Khat+1, gap); if gap < 3 fprintf('谱峰梯度偏小，MDL可能漏源，建议改用AIC或检查数据质量\n'); end end

3dB 是我常用的经验阈值。如果第 Khat 个峰和第 Khat+1 个峰只差不到 3dB，说明 MUSIC 谱上存在一个接近真实源的峰，而 MDL 把它归成了噪声，源数很可能被低估了。反过来，如果第 Khat+1 个峰比第 Khat 个峰低 10dB 以上，基本可以确认 MDL 的结果是干净的。

我自己用这招最值回票价的一次，是暗室实测里 0dB 信噪比下 MDL 给了 2，MUSIC 谱上第 3 个峰只比第 2 个低了 2.8dB。当时差点当成旁瓣忽略，后来补了一次窄带滤波重测，确认第 3 个峰是真实目标。从那以后，任何源数目估计做完，我都会顺手打印 Khat 和这个 gap 值。两个数字并排放在一起，比任何单一准则都让人放心，希望帮到你。

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