AIC、MDL、MEVARC源数目估计对比：选型、实现与实测避坑指南

1. 源数目估计这件事，到底在解决什么问题

阵列信号处理里有个绕不开的前置环节：来波方向估计算法跑之前，你得先知道空间里有几个信号源。MUSIC、ESPRIT这些子空间类算法，本质上都是把协方差矩阵的特征空间劈成信号子空间和噪声子空间两半，劈的位置就是源数目。数目估错了，噪声特征向量混进信号子空间，谱峰直接糊掉；数目估多了，信号子空间里掺进噪声，分辨率断崖式下跌。所以源数目估计不是可选项，是整套流程的地基。

工程上常见的做法是看特征值曲线的拐点，肉眼看拐点这事在实验室里凑合能用，上了实测数据就完全靠不住——低信噪比、小快拍、相干源这几个条件随便来一个，特征值曲线就变得又平又糊，拐点位置漂移得厉害。这时候就需要信息论准则类的算法来兜底，AIC、MDL、MEVARC这三个是文献里出现频率最高的，也是我在实际项目里反复对比过的三套方案。

这篇内容面向的是已经写过阵列处理代码、但被源数目估计坑过的朋友，也适合刚接触子空间算法、想搞清楚“为什么我的MUSIC谱峰数量不对”的读者。我会把三个算法的推导逻辑、参数选择、代码实现、实测踩坑全部摊开讲，重点放在什么条件下该选哪个算法以及为什么你的估计结果总是偏大或偏小这两个实际问题上。所有代码基于MATLAB，但思路换成Python、C++一样成立。

先说结论性的判断，方便你带着预期往下读：AIC在大快拍下过估计概率高，MDL在低信噪比下欠估计概率高，MEVARC在相干源场景下表现相对稳，但它的惩罚项设计决定了它对快拍数更敏感。这三个没有绝对优劣，选型取决于你的信噪比区间、快拍数和源之间的相关性。

2. 三个算法的核心逻辑拆解与选型依据

2.1 AIC和MDL的共同骨架：特征值似然加惩罚项

AIC（赤池信息准则）和MDL（最小描述长度）在源数目估计里的用法，本质上是同一套框架的两个变体。先把阵列协方差矩阵做特征分解，得到M个特征值按降序排列，然后对每个候选源数目k从0到M-1，计算一个准则函数，取最小值对应的k作为估计结果。

准则函数的一般形式是：

J(k) = -N * (M-k) * ln( 几何平均(后M-k个特征值) / 算术平均(后M-k个特征值) ) + 惩罚项(k)

前半部分是对数似然项，衡量的是“假设只有k个源时，剩下的M-k个特征值有多接近相等”。如果剩下的特征值确实都是噪声，它们的几何平均和算术平均会很接近，这个比值趋近1，对数项趋近0；如果剩下的特征值里还混着信号，几何平均会被小特征值拉低，比值变小，对数项变成负数，绝对值变大。所以这一项是在奖励“把信号都归到信号子空间”。

惩罚项是防止过拟合的刹车。AIC的惩罚项是2 * k * (2M - k)，MDL的惩罚项是0.5 * k * (2M - k) * ln(N)，N是快拍数。差别就在这个ln(N)上——MDL的惩罚随快拍数对数增长，快拍越多惩罚越重，所以MDL在大快拍下更保守，不容易过估计；AIC的惩罚与快拍数无关，快拍一多，似然项主导，就容易把噪声也当成信号，过估计概率上升。

这个推导我建议你自己在纸上推一遍，因为很多实现里的系数写错了你根本看不出来。我见过不止一份开源代码把AIC惩罚项的系数写成k*(2M-k)漏了2，结果估计结果系统性偏小。推导的关键在于对数似然函数的二阶展开，AIC取的是渐近无偏的惩罚，MDL取的是贝叶斯准则下的描述长度，两者在N趋于无穷时的收敛性不同——AIC不是一致估计，MDL是一致估计。所谓一致估计，就是快拍数趋于无穷时估计误差趋于零的概率为1，AIC做不到这一点，它有过估计的正概率下界。

2.2 MEVARC的差异化设计：不用全部特征值

MEVARC（Modified EigenVector-based Akaike Information Criterion的变体，不同文献里缩写略有出入，我这里指的是基于修正特征向量的准则）的思路和前两个不太一样。它不直接对全部特征值做似然比，而是构造一个基于特征向量的统计量，核心思想是利用信号子空间和噪声子空间的正交性来构造检验量。

具体做法是：对每个候选k，取前k个特征向量构成信号子空间估计，然后计算剩余特征向量在这个子空间上的投影能量，构造一个类似F统计量的比值，再套上AIC形式的惩罚。它的优势在于对特征值分布的假设更弱，不要求噪声特征值严格相等，所以在色噪声或者特征值扩散比较严重的场景下，比纯特征值准则更稳。

但MEVARC的代价是计算量更大，而且它对特征向量的估计精度敏感。快拍数少的时候，特征向量本身估计就不准，投影能量的统计特性会偏离理论值，准则函数的最小值位置会漂。我在一个8阵元、快拍数200的仿真里对比过，MEVARC的估计方差比MDL大，但在相干源场景下它的均值更接近真值。这个取舍很典型：你要稳定性还是要无偏性，取决于你的下游算法对过估计和欠估计哪个更敏感。

2.3 选型的三个判断维度

我把选型逻辑整理成三个维度，你在实际项目里按这个顺序过一遍基本不会选错：

这个表不是绝对的，但方向是对的。快拍少的时候MDL的ln(N)惩罚不够重，反而AIC的固定惩罚相对合适；信噪比低的时候似然项本身不可靠，MDL的重惩罚能压住虚警；相干源会让协方差矩阵秩亏，特征值准则类方法全部失效，这时候只能靠MEVARC这类基于特征向量结构的准则来救。

注意：如果你的源之间存在完全相干（比如多径导致的复制信号），AIC和MDL的估计结果会系统性偏小，因为它们看到的是一个“合并”后的特征值。这时候要么先做解相干处理（空间平滑），要么直接上MEVARC，别在AIC/MDL上浪费时间调参。

3. 实操过程：从协方差矩阵到准则曲线

3.1 数据预处理与协方差矩阵估计

整个流程的第一步是拿到协方差矩阵。假设你的接收数据矩阵是X，维度是M x N，M是阵元数，N是快拍数，那么样本协方差矩阵就是R = X * X' / N。这一步看起来简单，但有两个坑。

第一个坑是是否去均值。阵列处理里通常假设信号零均值，但实测数据往往有直流偏置，不去均值的话第一个特征值会异常大，直接导致估计结果偏大1。我的习惯是无条件先减掉每行的均值，代价几乎为零，收益是避免一类系统性偏差。

第二个坑是前向-后向平均。如果阵列是均匀线阵且信号是实信号或者圆对称复信号，用R = (X*X' + J*conj(X)*conj(X)'*J) / (2N)做前向-后向平均，等效快拍数翻倍，特征值估计的方差能降一半。J是反对角置换矩阵。这个操作在低快拍场景下效果非常明显，我实测过N=50的情况，做了前后向平均后MDL的欠估计概率从30%降到12%左右。

% 协方差矩阵估计，含去均值和前后向平均 X = X - mean(X, 2); % 按行去均值 N = size(X, 2); M = size(X, 1); R = (X * X') / N; J = fliplr(eye(M)); R_fb = (R + J * conj(R) * J) / 2; % 前后向平均

这里conj(R)而不是R'，因为前后向平均的公式里是共轭不是共轭转置，很多实现这里写错，导致结果和理论对不上。你可以用ishermitian(R_fb)验证一下，结果应该是1。

3.2 特征分解与特征值排序

[V, D] = eig(R_fb)拿到特征向量矩阵V和特征值对角阵D，然后lambda = sort(diag(D), 'descend')降序排列。这里有个细节：eig返回的特征值不保证顺序，必须显式排序，而且特征向量矩阵V的列顺序要和特征值对应，如果你后面要用MEVARC，V的列顺序必须跟着lambda一起排。

[V, D] = eig(R_fb); [lambda, idx] = sort(diag(D), 'descend'); V = V(:, idx); lambda = real(lambda); % 数值误差可能引入极小虚部，取实部

real(lambda)这一步别省。复数协方差矩阵的特征值理论上是实数，但数值计算会带出1e-16量级的虚部，后面取对数的时候如果虚部不为零，log会返回复数，准则函数比较大小就乱套了。我踩过这个坑，排查了半天才发现是虚部惹的祸。

3.3 三个准则函数的实现

AIC和MDL的实现放在一个函数里，因为骨架一样只差惩罚项：

function [k_aic, k_mdl] = aic_mdl_estimate(lambda, N) M = length(lambda); k_aic = 0; k_mdl = 0; J_aic = inf; J_mdl = inf; for k = 0:M-1 noise_eig = lambda(k+1:M); if any(noise_eig <= 0) continue; % 非正特征值跳过，防止log报错 end geo_mean = exp(mean(log(noise_eig))); ari_mean = mean(noise_eig); L = -N * (M-k) * log(geo_mean / ari_mean); pen_aic = 2 * k * (2*M - k); pen_mdl = 0.5 * k * (2*M - k) * log(N); if L + pen_aic < J_aic J_aic = L + pen_aic; k_aic = k; end if L + pen_mdl < J_mdl J_mdl = L + pen_mdl; k_mdl = k; end end end

几个实现要点。noise_eig <= 0的判断必须有，样本协方差矩阵在快拍数小于阵元数时会出现零特征值甚至负特征值（数值误差），log直接报错或者返回-Inf，准则函数就废了。geo_mean用exp(mean(log(...)))而不是prod(...)^(1/n)，前者数值稳定性好得多，后者在特征值动态范围大时容易溢出或下溢。

MEVARC的实现稍微复杂一点，需要用到特征向量：

function k_mevarc = mevarc_estimate(V, lambda, N) M = length(lambda); k_mevarc = 0; J_min = inf; for k = 0:M-1 if k == 0 proj_energy = 0; else Us = V(:, 1:k); Un = V(:, k+1:M); proj_energy = norm(Us' * Un, 'fro')^2; end stat = proj_energy / (M - k + eps); pen = 2 * k * (2*M - k); J = N * stat + pen; if J < J_min J_min = J; k_mevarc = k; end end end

norm(Us' * Un, 'fro')^2衡量的是信号子空间和噪声子空间的正交性偏离程度，理想情况下应该接近0。eps防止除零。这个统计量的理论分布我没有严格推导过，但实测下来它的准则曲线比AIC/MDL更平滑，最小值位置更明确，代价是计算量随M增长得快——Us' * Un是k x (M-k)的矩阵乘法，k从0遍历到M-1，总复杂度是O(M^3)量级，M=16的时候还能接受，M=64就明显慢了。

3.4 准则曲线的可视化与人工复核

我强烈建议你把准则曲线画出来，不要只看最小值。三个算法的曲线叠在一张图上，你能直观看到最小值的位置是否一致、曲线是否平坦、有没有多个局部极小。实测数据里经常出现AIC和MDL估计差2的情况，看曲线就知道是似然项在某个k处有个浅坑，还是惩罚项把最小值推偏了。

k_range = 0:M-1; figure; plot(k_range, J_aic_vec, '-o', k_range, J_mdl_vec, '-s', k_range, J_mevarc_vec, '-^'); legend('AIC', 'MDL', 'MEVARC'); xlabel('候选源数目 k'); ylabel('准则函数值'); grid on;

曲线平坦的区域意味着几个k的准则值差别在数值误差量级，这时候最小值的位置不可信，得结合其他信息判断。我一般的做法是：如果最小值附近两个点的准则值差小于总动态范围的5%，就认为估计不可靠，标记出来人工复核或者换算法。

4. 实测踩坑记录与排查速查表

4.1 估计结果系统性偏大1的三种原因

这是最常见的症状，AIC尤其容易犯。排查顺序按概率从高到低：

第一种，没去均值。直流偏置贡献了一个额外的“源”，特征值谱上第一个值明显突出，所有准则都会把它算进去。验证方法：看第一个特征值和第二个的比值，如果大于10且你的信号没有这么强的直达波，基本就是直流。

第二种，前后向平均用错。前面说过conj(R)和R'的区别，用错了相当于对协方差矩阵做了转置共轭，特征值不变但特征向量全乱，MEVARC会直接崩，AIC/MDL可能表现为偏大。验证方法：R_fb应该满足Hermitian性，norm(R_fb - R_fb', 'fro')应该接近0。

第三种，快拍数N传错了。AIC的惩罚项不含N，但似然项含N，N传大了似然项权重过高，过估计；MDL的惩罚项含ln(N)，N传大了惩罚过重，欠估计。我见过把size(X,1)当成N传进去的，阵元数当快拍数，结果完全不可信。验证方法：打印一下N，和你的数据维度对一遍。

4.2 低信噪比下MDL欠估计的补救

信噪比低于0dB的时候，MDL的ln(N)惩罚会把真实源也压掉，估计结果偏小。这时候有几个补救方向：

• 增加等效快拍数：前后向平均、或者对数据做分段平均（如果信号是分段平稳的）。
• 降维：如果阵元数远大于预期源数，先做波束空间变换把维度降到2到3倍预期源数，特征值谱的动态范围会拉开，MDL的判别力提升。
• 换准则：这个信噪比区间MEVARC通常比MDL稳，代价是计算量。

我做过一组对比仿真，8阵元、3个源、信噪比从-5dB扫到20dB、快拍数固定200，每个点跑500次蒙特卡洛，统计正确估计概率：

这张表的信息量很大。低信噪比段MEVARC明显占优，10dB以上MDL反超，20dB时AIC因为过估计正确率反而掉下来。所以如果你的应用场景信噪比波动大，单一算法不够，得做自适应切换——用粗略的信噪比估计（比如最大特征值和最小特征值的比值）来决定用哪个准则。

4.3 相干源场景的识别与处理

相干源的判断方法：看特征值谱的“台阶”结构。独立源的特征值谱是明显的“大值群加小值群”，相干源会让某个大特征值异常突出，后面的信号特征值被压到接近噪声水平，看起来像源数目少了。这时候AIC和MDL都会欠估计。

处理路径有两条。第一条是空间平滑，把均匀线阵分成若干重叠子阵，子阵协方差矩阵平均，恢复秩。子阵长度的选择有讲究：子阵数要大于等于源数，子阵长度要大于等于源数的2倍，这两个约束经常冲突，需要折中。第二条是直接用MEVARC，它对特征值分布的假设弱，相干源下虽然也会受影响但比前两个轻。

提示：空间平滑会损失有效孔径，DOA估计的分辨率会下降。如果分辨率是硬指标，优先考虑MEVARC或者别的解相干手段，别轻易上平滑。

4.4 常见问题速查表

5. 参数选择背后的计算逻辑与经验值

5.1 快拍数下限的估算

快拍数不是随便定的，它和阵元数、信噪比、期望的估计精度都有关系。一个经验公式是N >= 2 * M / SNR_linear，SNR_linear是线性信噪比。比如8阵元、10dB信噪比（线性值10），N至少16，但这是理论下限，实际要留3到5倍余量，我一般取N >= 10 * M / SNR_linear，8阵元10dB就是80，实测下来这个量级MDL和MEVARC都能稳定工作。

这个公式的推导来自特征值估计的方差分析：样本协方差矩阵的特征值估计误差方差正比于lambda^2 / N，信噪比越低，信号特征值和噪声特征值的差距越小，需要更多快拍才能把两者分开。系数10是经验值，来自我做的蒙特卡洛扫描，保证正确估计概率在0.9以上。

5.2 阵元数与源数的比例约束

子空间类算法要求阵元数大于源数，这是硬约束。但源数目估计本身对比例有更细的要求：M / k >= 2是底线，M / k >= 3比较稳。原因在于噪声子空间的维数是M - k，维数太低的时候噪声特征值的统计涨落会污染似然比统计量，准则函数的最小值位置漂移。如果你的应用里源数可能接近阵元数的一半，考虑增加阵元或者用稀疏阵列扩展自由度。

5.3 准则函数比较时的数值细节

三个准则函数的量纲不同，AIC和MDL的似然项是N量级，MEVARC的统计量也是N量级但系数不同，直接比较绝对值没有意义，只能各自取最小值。如果你想把三个结果融合（比如投票），不要比较准则值，比较估计出的k，多数投票或者取中位数。

另外，log(N)在N很大时增长缓慢，MDL的惩罚项在N从1000增到10000时只增加0.5*k*(2M-k)*2.3，而似然项随N线性增长，所以超大快拍下MDL也会趋向过估计，只是比AIC晚得多。这个 crossover 点大概在N > exp(2)量级之外，实际数据很难达到，可以忽略，但理论上要知道。

6. 把三个算法串成一套自适应流程

单算法跑通之后，我在实际项目里更常用的是一个自适应流程：先算特征值谱，用最大最小特征值比值粗估信噪比，高信噪比走MDL，低信噪比走MEVARC，中间地带两个都跑，结果一致就采纳，不一致就看准则曲线的最小值锐度，选锐度高的那个。相干源检测用特征值谱的台阶结构做前置判断，检测到相干就直接走MEVARC加空间平滑的组合。

这套流程的代码量不大，核心是一个switch加几个if，但省去了大量人工调参。我把它封装成一个estimate_source_number(X, options)函数，options里可以指定强制算法、是否做前后向平均、是否画曲线，返回估计值和诊断信息。诊断信息里包含三个算法的原始估计、准则曲线的最小值锐度、粗估信噪比，出问题的时候看诊断信息比看结果有用得多。

最后分享一个我调试时常用的小技巧：构造一个已知源数目的仿真数据，把信噪比和快拍数设成你实测数据的估计值，跑一遍三个算法，看它们的表现是否符合预期。如果仿真里表现正常但实测偏了，问题一定在数据预处理或者参数传递上，不在算法本身。这个对照实验能把排查范围缩小一半，比盲目调参高效得多。