简介压缩感知稀疏贝叶斯算法实现包包含SBL、TSBL、TMSBL三种典型方法面向信号处理、压缩感知与稀疏重构方向的研究者、工程师及高年级学生解决从有限线性测量中恢复稀疏或动态变化信号的问题。资源共15个文件压缩包大小479KB以MATLAB脚本11个m文件为主体辅以PDF/docx使用指南和txt说明文件涵盖三种算法的核心实现、对比演示与快速上手文档。已有929人浏览学习代码经过亲测可用对复现实验结果与二次开发具有直接参考价值。压缩包内除算法主程序外还提供多个demo演示脚本如时变信号、相同向量、不同信噪比场景并配备“3分钟/1分钟快速使用”说明与ReadMe可帮助使用者快速理解调用方式、参数设置与效果对比节省算法调试时间适合用于课程设计、科研实验或工程预研。1. 压缩感知稀疏贝叶斯算法SBL/TSBL/TMSBL到底帮你解决什么问题当手里只有 M 个观测、要还原 N 个点上的稀疏系数时压缩感知稀疏贝叶斯算法是少数能把“不确定度”一起算出来的估计器在雷达测向、信道估计、频谱感知和稀疏成像里反复被用。最常见的三个版本是 SBL、TSBL 和 TMSBLSBL 对付单个快照/单个任务TSBL 利用同一段连续时间内的快照相关性TMSBL 面向多个彼此独立但共享支撑集的任务。这套算法最吸引人的地方是不需要人为指定稀疏度 K而是靠迭代把多余原子自动关掉低信噪比下的支撑集往往比 OMP 稳。适合你如果你手里字典已经确定、用 OMP/LASSO 试过但支撑集抖动或者想把估计的不确定度也带出来。2. 用 SBL 跑通第一个稀疏重构高斯先验、ARD 机制与循环更新2.1 ARD 机制SBL 为什么能自动确定稀疏度SBL 和 OMP、LASSO 的思考方式不一样。OMP 是“选一个原子抠掉贡献再选下一个”LASSO 是调一个正则系数 λ 逼出一串零SBL 则是把所有原子都先放进模型里给每个原子配一个未知方差 γ_i然后反复迭代让大多数 γ_i 自己衰减到接近 0这个过程叫自动相关决定ARD。真正留系数 σ² 还是先验方差 γ_i 决定支撑集不需要交叉验证去扫 K。这个机制的核心来自高斯先验的层次化建模观测 y Φx e假设 e 是零均值高斯x 的先验是零均值高斯方差为 γ_i。迭代后后验均值 μ 就是我们要的稀疏系数估计后验方差 Σ 的对角线给出了每一个系数的不确定度。γ_i 的更新式是 |μ_i|² Σ_{i,i}意思是“这个原子既要有比较大的幅度又要是确定的”两者都大的原子才值得留在支撑集里。大多数不相关的原子幅度小、后验方差大γ_i 会一口气掉到 1e-6 以下自然不会被选上。这个机制带来两个实际好处。一是稀疏度 K 是隐式确定的不做模型选择二是结果不是一个点估计而是一个完整后验你能画出每个系数上的置信区间。代价是迭代非线性对初值和噪声方差敏感后面章节我会重点说这两个坑。2.2 最小可跑的 SBL 实现与参数说明单快照 SBL 的迭代可以压缩成下面几步算后验协方差算后验均值更新 γ更新噪声方差 σ²。后验协方差那一步如果直接对 N×N 矩阵求逆字典一大就爆内存常见做法是用 Woodbury 公式把求逆对象换成 M×M 的矩阵。下面的实现可以直接存成mySBL.m跑起来function [x_hat, gamma, sigma2] mySBL(y, Phi, gamma0, sigma2_0, tol, maxIter) % 输入: % y : M x 1 观测向量 % Phi : M x N 字典矩阵建议先做列归一化 % gamma0 : N x 1 先验方差初值没有先验信息就取 ones(N,1) % sigma2_0: 噪声方差初值标量通常取 var(y)/10 量级 % tol : gamma 相对变化阈值例如 1e-6 % maxIter : 最大迭代次数例如 500 % 输出: % x_hat : N x 1 稀疏系数的后验均值 % gamma : N x 1 学习到的先验方差接近0的位置就是非支撑集 % sigma2 : 估计出的噪声方差 [M, N] size(Phi); gamma gamma0(:); sigma2 sigma2_0; Im eye(M); for it 1:maxIter gamma_old gamma; sigma2_old sigma2; D diag(gamma); PhiD Phi * D; % 后验协方差: Sigma D - D*Phi*inv(sigma2*I Phi*D*Phi)*Phi*D A sigma2 * Im PhiD * Phi; Sigma D - PhiD * (A \ PhiD); % 后验均值 mu Sigma * Phi * y / sigma2; % ARDgamma 更新 后验均值平方 后验方差 gamma_new abs(mu).^2 real(diag(Sigma)); % 噪声方差更新分母是有效自由度要保护一下 noise_num norm(y - Phi * mu)^2; noise_den M - sum(real(diag(Sigma)) ./ gamma eps); sigma2_new noise_num / noise_den; gamma max(gamma_new, 1e-12); sigma2 max(sigma2_new, 1e-12); % 用相对变化作为收敛判据避免量级差异 if norm(gamma - gamma_old) / norm(gamma_old) tol break; end end x_hat mu; end这段代码里最关键的是Sigma D - PhiD * (A \ PhiD)。如果不这么做直接对sigma2*eye(N) Phi*Phi*D求逆在 N 是几千上万时内存直接爆炸。Woodbury 版本只需要维护 M×M 的矩阵 AM 是观测数量在压缩感知场景里 M 通常只有几百速度会快很多。配合一个生成随机字典、已知稀疏系数的测试脚本马上能验证支撑集找回得对不对% 测试 mySBL 的最小脚本 N 200; M 60; K 5; Phi randn(M, N); Phi Phi ./ vecnorm(Phi); % 列归一化是 SBL 的默认前提 x_true zeros(N, 1); idx randperm(N, K); x_true(idx) randn(K, 1); y Phi * x_true 0.05 * randn(M, 1); [x_hat, gamma, sigma2] mySBL(y, Phi, ones(N,1), 0.05*var(y), 1e-6, 500); % 对比支撑集 stem(abs(x_hat), b); hold on; stem(idx, 0.2*ones(K,1), r); legend(重构系数, 真实支撑);这里Phi ./ vecnorm(Phi)是我每次拿到新字典第一件事没有列归一化就去跑 SBL结果几乎不可复现。现实的字典例如傅里叶字典、随机卷积字典列范数不均匀很常见而 γ_i 的更新又和列能量耦合列能量大的原子天然容易被选中。2.3 收敛判据和三个参数陷阱第一个陷阱是噪声方差初值。σ² 设得太大会让所有 γ_i 一起往小走迭代很久都选不出支撑σ² 设成 0 会让 A 变成奇异矩阵并且后验均值公式失效。我一般先用var(y)/M做一次粗略估计再乘 0.1 到 0.2 作为初值。第二个陷阱是 γ 的收敛判据用绝对变化。γ 收敛后量级可能从 1e-3 掉到 1e-9绝对变化看起来早就“稳定”了但支撑集还没定改成相对变化norm(gamma-gamma_old)/norm(gamma_old)会把尾部的小量级变化也照顾到。第三个陷阱是 maxIter 太少。SBL 前几十轮变化很大后一百轮是在微调支撑边界如果只给 200 轮支撑集边界常常还会带着一个伪原子。如果发现迭代结束但支撑集还在漂不要急着加迭代次数先看 σ² 估计值是不是比真实噪声大很多。常见做法是先固定 σ² 跑 50 轮让 γ 先把支撑集轮廓逼出来再放开 σ² 一起联合更新。这比盲调迭代次数管用得多。3. 从 SBL 到 TSBL时间相关性怎么建模、B0 矩阵怎么更新3.1 为什么多快照数据不能逐列独立处理单快照 SBL 只处理一列观测。实际系统里很少只有一组观测阵列测向会连续采很多快拍信道估计会有多个导频符号雷达回波本身也是一段时间序列。最直接的做法是每个快照独立跑一次 SBL然后把支撑集做并集或投票。问题在于快照之间的噪声独立但支撑集会抖动低信噪比时同一个目标在前后两帧里可能被选到相邻两个格点上投票后反而把两个格点都留下分辨率被撑宽。TSBL 做的事情是把 L 个快照组成 Y(M×L)让所有快照共享同一组 γ_i但不再假设每个快照之间的系数完全独立。建模上引入一个 L×L 的时间相关矩阵 B0用来描述同一个源在相邻快照间的幅度变化模式。直观地说如果一个目标在第 t 帧出现它在第 t1 帧也不该凭空消失B0 就是用来编码这种连续性的。对时间序列数据B0 往往接近 Toeplitz 结构对角线和次对角线都很强。3.2 TSBL 的模型结构行共享 γ、列相关由 B0 承担TSBL 把 X 看成一个 N×L 的系数矩阵。X 的第 i 行是第 i 个潜在源在 L 个快照里的幅度变化先验是零均值高斯协方差为 γ_i B0。这里 γ_i 控制这个源在整体上是否活跃B0 控制活跃源的时间波形相关结构。观测模型还是 Y ΦX E每个时刻共用一个测量矩阵 Φ。后验计算比 SBL 繁琐一些。按行拆开后设第 i 行的后验均值为 m_iL×1 列向量后验协方差对应第 i 行的块为 Σ_iL×L。超参数更新式变成γ_i 更新γ_i (1/L) × (m_i B0^{-1} m_i tr(B0^{-1} Σ_i))B0 更新B0 (1/N) Σ_i (m_i m_i Σ_i)B0 更新完之后一般还会做一次结构约束常见的是把 B0 强制成 Toeplitz 矩阵或者只保留对角线附近几条带。这样做的理由是时间相关性不应该破坏“均匀时间轴”的结构自由更新的 B0 在 L 较大时会被少数几个快照带偏。3.3 TSBL 的核心循环在 SBL 基础上只多一个 B0下面的代码是 TSBL 的核心循环骨架重点看和 SBL 不同的三行计算 m_i、更新 γ_i、更新 B0。为了不把篇幅拉得过长后验协方差块的计算做了近似简化实际使用请以你拿到的源码包内实现为准function [X_hat, gamma, B0] myTSBL(Y, Phi, gamma0, B0_0, sigma2, tol, maxIter) % Y : M x L 多快照观测L 是快照数 % Phi : M x N 字典 % B0_0 : L x L 时间相关矩阵初值通常取单位阵 % sigma2 : 噪声方差可以先固定 % 输出 X_hat: N x L 重构的系数矩阵 [M, N] size(Phi); [~, L] size(Y); gamma gamma0(:); B0 B0_0; B0inv inv(B0); Im eye(M); % 先按多快照 SBL 的思路算出整体后验均值和协方差对角线块 PhiD Phi * diag(gamma); A sigma2 * Im PhiD * Phi; Sigma_full D - PhiD * (A \ PhiD); % N x N这里没展开块结构 Mu PhiD * (A \ Y); % N x L 的后验均值 gamma_new zeros(N, 1); B0_num zeros(L, L); for i 1:N m_i Mu(i,:).; % L x 1 % 实际实现里 Sigma_i 要从完整后验协方差中抽取这里是占位写法 Sigma_i Sigma_full(i,i) * eye(L); % 简化的对角近似 gamma_new(i) (m_i * B0inv * m_i trace(B0inv * Sigma_i)) / L; B0_num B0_num m_i * m_i Sigma_i; end B0 B0_num / N; % 时间序列场景一般再加 Toeplitz 约束去掉自由 B0 的噪声 gamma max(gamma_new, 1e-12); X_hat Mu; end注意Sigma_full(i,i) * eye(L)只是为了把循环结构讲清楚严格 TSBL 的 Σ_i 是后验协方差的第 i 行对应的 L×L 块需要在计算后验协方差时保留 Kronecker 结构或分块索引。实际成熟实现里会用kron或稀疏分块避免显式构造 N×N 矩阵这也是为什么拿到源码包后先不要改 B0 更新那一段先用默认参数跑通再改。TSBL 的参数设置和 SBL 有明显差别。B0 初值用单位阵最稳如果直接用toeplitz(0.9.^(0:L-1))这种强相关初值头几轮迭代会把所有行都暂时拉成“活跃”收敛变慢。γ 初值还是全 1。噪声方差如果信噪比大致已知直接固定不要更新TSBL 在 σ² 未知时联合估计容易把时间相关性吸收到噪声里导致 B0 被估计成近对角阵。4. TMSBL 多任务联合重构共享支撑集的规则与别忘了先归一度4.1 TMSBL 和 TSBL 的区别任务独立、支撑共享TMSBL 里的“多任务”和 TSBL 的“多快照”在数据形式上很像都是一组观测 Y(M×L)但语义不同。TMSBL 假设每个任务是一组独立的压缩感知观测比如多个用户上报同一频段的压缩采样或者多个极化通道对同一目标成像。任务之间没有时间上的连续关系不能假设幅度波形相关但可以假设它们共享同一个支撑集——目标位置一样、频点一样、原子索引一样。建模上TMSBL 给每个任务单独保留自己的稀疏系数 x_l 和噪声方差 σ_l²所有任务共享同一组 γ_i。γ_i 的更新变成对所有任务取平均γ_i (1/L) Σ_l (μ_{l,i}² Σ_{l,ii})这里 μ_{l,i} 是第 l 个任务里第 i 个原子的后验均值Σ_{l,ii} 是第 l 个任务后验协方差对角线。该式的含义是一个原子只有“在所有任务里都稳定地非零”才被保留。如果某个原子只在两个任务里大、其他任务里小平均之后 γ_i 会被压下去这就是 TMSBL 抑制单任务假目标的机制。4.2 可运行的 TMSBL 最小实现思路TMSBL 的实现可以复用 SBL 里的 Woodbury 求逆区别只是对每个任务分别做一次后验计算然后合并 γ 更新。下面的代码是按“共享 γ、独立 σ²”的常见形式写的function [X_hat, gamma, sigma2_list] myTMSBL(Y, Phi, gamma0, sigma2_init, tol, maxIter) % Y : M x L每个任务一列观测 % Phi : M x N所有任务共享同一个字典 % sigma2_init : 1 x L 或标量每个任务可单独给噪声初值 [M, N] size(Phi); L size(Y, 2); gamma gamma0(:); sigma2_list ones(L, 1) * sigma2_init(:); % 每个任务独立噪声方差 Im eye(M); for it 1:maxIter gamma_old gamma; gamma_sum zeros(N, 1); for l 1:L D diag(gamma); PhiD Phi * D; A_l sigma2_list(l) * Im PhiD * Phi; % 第 l 个任务的后验协方差和均值 Sigma_l D - PhiD * (A_l \ PhiD); Mu_l Sigma_l * Phi * Y(:, l) / sigma2_list(l); % 累加本任务对 gamma 的贡献 gamma_sum gamma_sum abs(Mu_l).^2 real(diag(Sigma_l)); % 第 l 个任务噪声方差更新 noise_num norm(Y(:, l) - Phi * Mu_l)^2; noise_den M - sum(real(diag(Sigma_l)) ./ gamma eps); sigma2_list(l) max(noise_num / noise_den, 1e-12); end gamma max(gamma_sum / L, 1e-12); if norm(gamma - gamma_old) / norm(gamma_old) tol break; end end X_hat zeros(N, L); % 最后再对每个任务算一次后验均值作为输出这里省略重复计算 end这段代码把 TMSBL 和 SBL 的区别拆得很清楚γ 更新变成了跨任务平均σ² 每个任务单独维护。实现上不需要额外引入 B0这是 TMSBL 的最简形式。如果采用完整 TMSBL还要在任务之间估计相关矩阵 B但很多时候共享 γ 的普通多任务 SBL 已经能稳定支撑集B 带来的增益只有任务间统计相关强时才明显。4.3 TMSBL 必调的三个参数第一个是任务间幅值归一化。各个任务如果来自不同通道、不同阵元能量可能差 10 倍。直接丢进 myTMSBL能量大的任务主导 γ 平均小能量任务的支撑信息被淹没。处理方式是每个任务先除以自身的 Frobenius 范数或者给每个任务在观测端做一个标定系数这个步骤在压缩感知多任务里比算法本身更重要。第二个是 γ 的剪枝阈值。迭代结束后 γ 接近 0 的原子不会完全是 0通常留着一批 1e-6 量级的尾巴。取支撑集时不要用固定绝对阈值建议先看 γ 的直方图找那个“突然掉下去”的拐点或者按能量占比保留前若干个原子。第三种做法是按后验均值 μ 的模做阈值但 γ 剪枝更直接。第三个是任务数 L 不能太小。TMSBL 的收敛质量依赖“跨任务投票”的统计稳定性L 只有 2 或 3 时一个任务的异常噪声点就可能把某几个 γ 拉高产生伪原子。我一般建议至少 8 个任务再上 TMSBL任务数不够时退回去用 TSBL 或单任务 SBL 反而更可控。5. SBL/TSBL/TMSBL 避坑实录支撑集错位、收敛黑匣子与多快照对齐排查5.1 支撑集错位重构出的峰不在真实原子附近现象真实稀疏系数在索引 10、50、70 三个位置重构后峰在 11、49、71每个峰都偏了一点相邻原子也带小尾巴。用 OMP 时这个现象偶尔出现换成 SBL 后以为能消失结果仍然有。原因字典原子之间有相关性尤其网格越细相邻原子波形相似度越高SBL 的 γ 更新会把一部分能量摊到相邻原子对应的 γ_i 上。另一个常见原因是我反复提的列能量未归一化列范数大的原子天然更容易被选中。解决先做字典列归一化再看字典的 Gram 矩阵相关性如果相邻原子的相关系数超过 0.9说明是字典设计问题不是算法问题。网格间距尽量按瑞利分辨极限走过密网格对 SBL 的 ARD 机制不友好。还可以在 γ 更新之后加一步最小原子间距约束把距离太近且 γ 都不大的原子合并。5.2 收敛黑匣子迭代结束但 γ 还在缓慢漂移现象maxIter给到 1000tol 设到 1e-7程序正常结束但把前后两次的支撑集打出来会发现偶尔一个原子被换掉。换掉的原子 γ 值都在 1e-4 到 1e-5 之间说大不大说小不小。原因收敛判据用的是绝对变化γ 从 1e-3 变到 1e-4 的变化量和从 1e-7 变到 1e-8 的变化量在绝对值上差很多绝对判据在尾部失效。另一个原因是 σ² 和 γ 在联合迭代时互相拖拽边界原子一直处在“留下还是删掉”的临界状态。解决把判据换成相对变化并且额外加一条支撑集状态判据连续 20 轮支撑集索引集合不变才算收敛。还有一种工程做法是固定噪声方差 σ²只迭代 γ收敛速度会快很多支撑集也更稳。等 γ 基本稳定后再放开 σ² 做一两轮联合精调。5.3 多快照对齐TSBL 估计出的 B0 奇异或近对角现象TSBL 跑完B0 不是预想的时间相关矩阵对角线上是 1其他位置几乎全是 0仿佛算法完全没利用时间相关性。强制加 Toeplitz 约束后性能反而下降。原因多快照数据没有对齐。时间偏移、采样时钟不同步、或者各快照之间幅度差异太大都会让“同一个源的时间波形”在观测里错位统计相关被噪声抹平。B0 更新公式对错位非常敏感两帧波形错开 3 个采样点相关系数直接掉一半以上。解决TSBL 之前先做粗对齐用互相关找延迟再补偿然后对每一帧做能量归一化去掉幅度突变帧最后再看 B0 的次对角线如果仍然接近 0可以尝试把数据降采样或分帧再跑。时间相关性建模不是“多快照自动有效”的魔法数据本身要保证连续、等间隔、稳定。5.4 任务增益差太远TMSBL 的 γ 被单一任务主导现象TMSBL 恢复出的支撑集几乎等于把能量最强那个任务单独跑 SBL 的结果弱势任务里独有的真实原子没被选出来强任务里的伪原子还留着。原因前面 4.3 提到的任务幅值差异。TMSBL 的 γ 更新是算术平均一个任务能量大 100 倍它的 |μ|² 就比别人大 10 倍平均出来的 γ 基本被它一个人决定。解决对每个任务观测 Y(:,l) 除以各自的norm(Y(:,l))让任务能量一致。如果任务之间的信噪比差异仍然很大就把低信噪比任务的噪声方差初值单独设大避免它在 γ 平均里投出假票。这两个动作都要在算法外面做不要依赖 TMSBL 内部自动抵抗。6. 验证 TMSBL 结果三板斧误差、支撑集重合度和参数敏感性拿到结果不要只看重构误差一个数。误差小只能说明系数拟合得好不能说明支撑集对。做验证时我一般固定随机种子用真实支撑集和估计支撑集的重合度作为主指标est_idx find(gamma 1e-3); precision length(intersect(est_idx, idx)) / length(est_idx); recall length(intersect(est_idx, idx)) / K;precion 和 recall 都比重构误差更能反映测向、稀疏成像这类任务的实际收益。还有一个容易被忽略的验证维度是参数敏感性把 σ² 初值从 0.01 扫到 1把 maxIter 从 200 扫到 1000看支撑集重合度会不会抖动。抖得厉害说明你的数据场景在 ARD 机制的临界区需要多试几组初值或者改用固定噪声方差的模式。迁移到自己的数据时第一步永远是拿合成数据做冒烟测试用你真实字典 Φ随机生成 K 个系数加白噪声验证 SBL/TMSBL 能把支撑集找回来。第二步才是换成真实测量数据。换真实数据后最先出问题的通常是字典建模和噪声非高斯不要第一反应去改算法结构。我自己在阵列测向场景里翻过最狠的一次车就是把真实导向矢量字典的相位基准没对齐SBL 迭代了五百轮支撑集稳定在错误格点上换了列归一化都没用最后发现是信号模型里少写了一项。那之后我再也不跳过合成数据验证。这三兄弟本身是成熟的工具箱亲自测试往往只是第一步真正决定成败的是字典、初值和数据格式这三个外围环节。希望帮到你。本文还有配套的精品资源点击获取
