CoNMF高光谱解混：协同稀疏约束与sunsal初始化工程实践

简介：面向高光谱遥感解混研究的一份CoNMF算法资源包，适合遥感图像处理、目标检测与环境监测方向的科研人员及相关专业学生。内容包括约束非负矩阵分解的完整实现与演示流程，覆盖端元提取、丰度估计和混合像元分解等关键环节。压缩包共30个文件，以.m脚本为主，配合.mat光谱库数据、.asv/.bak备份及.bib参考文献，整体约20.21MB，便于按需调用。已有198人学习。通过示例可直接运行观察不同光谱纯度与噪声条件下的解混效果，并可结合SUNSAL、VCA等算法理解优化迭代策略，为土地覆盖分类、污染探测等遥感应用提供可复现的实验支撑。

1. CoNMF 高光谱解混：当 sunsal 的稀疏先验还不够时，试着让丰度行一起变稀疏

拿到一景机载高光谱影像，比如 224 个波段的 AVIRIS 数据，最常被问的不是分类精度，而是「这块地到底有哪些矿物、各占多少」。这就是解混：从混合像元里拆出端元光谱和丰度比例。传统做法先提端元再做丰度约束反演，但误差会在两步之间滚雪球。CoNMF 的思路是把端元和丰度放在同一个非负矩阵分解框架里调，同时给丰度矩阵加一个协同稀疏约束，而 sunsal 恰好是在这个框架里做初始化和对比验证的好搭档。这篇就按工程落地的顺序，把 CoNMF 的原理、MATLAB 入手代码、参数调节和常见失败画面串一遍，适合刚转向高光谱定量反演、对稀疏解混有基础但没跑通过完整流程的从业者。

2. 从 NMF 到 CoNMF：协同稀疏约束的逻辑与数学形式

2.1 丰度矩阵的行稀疏与像元稀疏不是一回事

高光谱解混的观测模型可以写成一矩阵分解问题：X ≈ E * A，其中 X 是 L×P 的光谱矩阵，L 是波段数，P 是像元数；E 是 L×K 的端元矩阵，K 是端元数量；A 是 K×P 的丰度矩阵。传统 NMF 只要求 E 和 A 非负，结果不唯一，解出来的端元常常没有物理意义。后来大家往目标函数里加稀疏惩罚，其中最有代表性的是 sunsal 这类基于变量分裂和增广拉格朗日的稀疏回归方法，它默认每个像元只由少数几个端元构成，对丰度矩阵逐列做稀疏约束。但这里有个容易被忽略的点：sunsal 的稀疏是逐像元的，它并不关心某个端元是否在整幅影像里都不出现。

CoNMF 的出发点正好在「行」这一维。它认为在一个场景内，某个端元要么至关重要，要么几乎处处不重要，所以丰度矩阵里不该存在「这行有点值、那行也有点值」的均匀分布。数学上就是把稀疏惩罚作用在丰度矩阵的行上，优先让一整行同时趋近于零，这叫行稀疏或协同稀疏。实现时常用两种正则项：l2,1 范数和 l2,0 范数。l2,1 是先把每行的 l2 范数加总，对行向量做软阈值；l2,0 直接统计非零行数量，更强的选择效果。实际代码里往往用 l2,1 替代 l2,0，因为它连续、可导、好优化，而且解出的行稀疏模式已经足够明显。

2.2 CoNMF 的优化目标与锚点约束

CoNMF 的目标函数一般写成下面这种形式：

min ||X - E A||_F^2 + λ * Σ_k ||A(k,:)||_2 约束条件：A >= 0，E >= 0，1^T A = 1^T

公式里||X - E A||_F^2是重建误差，保证分解后的 E、A 能还原原始观测；Σ_k ||A(k,:)||_2是协同稀疏惩罚，k 遍历所有端元，某个端元在全局不活跃时，对应行向量的 l2 范数会收缩到接近 0；λ 控制稀疏强度。约束A >= 0和1^T A = 1^T分别对应丰度的非负性和归一性，即每个像元的丰度比例加起来等于 1。这两个约束合称锚点约束，因为端元通常也被限制在正数范围，相当于把解固定在一个有物理意义的单纯形上。

与两步法最大的区别在于，CoNMF 对 E 的更新也要考虑 A 的状态。我在实际项目中通常把 E 的更新写成乘性规则：E <- E ⊙ (X A^T) / (E A A^T)，分母加一个极小量防止除零。代数上这一步是在最小二乘方向上做乘法步长，能自动保持非负性，实现成本比投影梯度低很多。A 的更新则是乘性更新之后，加一次行方向软阈值，再投影回非负和归一约束。整个迭代过程交替进行，端元和丰度在每一轮相互纠偏，这是 CoNMF 比「先 VCA 提端元、再 FCLS 求丰度」不易漂移的原因。

2.3 与 sunsal 的关系：先稀疏后协同

搞清楚 CoNMF 和 sunsal 的关系，比会调参更重要。sunsal 解决的是固定端元 E、未知丰度 A 的稀疏反演问题，目标函数是min ||E A - X||_F^2 + λ||A||_1，约束 A 非负，它的稀疏是面向单个像元的。CoNMF 则把 E 与 A 联合优化，稀疏约束放在 A 的行方向。两者不是替代关系，而是递进关系：先用 sunsal 算出 A 的稀疏初值，再进入 CoNMF 主循环迭代细化，是这套组合最常见的用法。如果代码库里只给了 CoNMF 的参数入口而没有配套初始化，用 sunsal 的结果替换随机初始化，收敛速度和端元准确性通常都会明显改善。

sunsal 的调用形式简单，但有几个参数值得理解。lambda控制稀疏强度，值越大丰度越稀疏；POSITIVITY开启非负约束；ADDONE控制是否施加和为 1 的约束。在给 CoNMF 提供初始 A 时，我一般会把ADDONE设为 no，同时在 CoNMF 主循环里统一加锚点约束，避免双重归一化造成丰度整体偏小。下面是三种方法的对比。

3. 用 sunsal 初始化 CoNMF：MATLAB 最小可复现脚本

3.1 数据准备：把高光谱影像整理成 L×P 矩阵

在 MATLAB 里跑这套流程，第一步是把影像从三维立方体转成二维矩阵。ENVI 格式的遥感数据通常有.hdr头文件和.dat数据文件，社区里常见的做法是用 enviread 这类函数读取，但若只是自己调试，可以先读成三维数组再 reshape。假设img是 R×C×L 的三维数组，执行下面代码就能得到 L×P 光谱矩阵：

[R, C, L] = size(img); X = reshape(img, R * C, L)'; % L×P，P = R*C X = double(X); % 防止 uint16 参与矩阵乘法截断 % 可选：去除水汽吸收波段，例如 1350~1450nm 和 1800~2000nm bad = [find(wl >= 1350 & wl <= 1450), find(wl >= 1800 & wl <= 2000)]; X(bad, :) = [];

这里 reshape 后转置得到 L×P 的原因是，解混算法统一以波段为行、像元为列。X的每一列是一个像元光谱，每一行是一个波段在所有像元上的灰度。水汽波段噪声大且吸收强，会干扰端元提取，预处理阶段删掉比在算法里加权重更直接。

3.2 CoNMF 主迭代的 MATLAB 代码

以下代码是实现 CoNMF 思想的一个工程简化版，亮点在于用 sunsal 初始化 A、用乘性更新保持非负性、再用行软阈值完成协同稀疏。把它保存成conmf_sunsal_demo.m可以直接跑。

function [E_est, A_est] = conmf_sunsal_demo(X, K, lambda, maxIter) % 输入：X 为 L×P 高光谱数据，K 为端元数 % lambda 为协同稀疏系数，maxIter 为最大迭代次数 % 输出：E_est 为 L×K 端元，A_est 为 K×P 丰度 if nargin < 4, maxIter = 200; end if nargin < 3, lambda = 0.01; end % 数据归一化到 [0,1]，让 lambda 有稳定尺度 X = X / max(X(:)); % 1. 用 PCA 均值离差选 K 个像元作为初始端元 [coef, score] = pca(X'); [~, idx] = max(abs(score(:,1:min(K,3))), [], 1); E_est = X(:, idx(1:K)); % 2. 用 sunsal 生成稀疏丰度初值 addpath('sunsal'); A_est = sunsal(E_est, X, 'lambda', lambda * 5, ... 'POSITIVITY', 'yes', 'ADDONE', 'no', 'verbose', 'no'); A_est = max(A_est, 0); prevCost = inf; for iter = 1:maxIter % 更新 E：乘性更新，保持非负 gradNum = X * A_est'; % L×K gradDen = E_est * (A_est * A_est'); % L×K E_est = E_est .* (gradNum ./ max(gradDen, 1e-10)); E_est = max(E_est, 0); % 非负 E_est = min(E_est, 1); % 反射率上界 % 更新 A：乘性更新 + 行软阈值协同稀疏 ANum = E_est' * X; % K×P ADen = E_est' * E_est * A_est; % K×P A_est = A_est .* (ANum ./ max(ADen, 1e-10)); rowNorm = sqrt(sum(A_est.^2, 2)); % K×1 shrink = max(1 - lambda ./ max(rowNorm, 1e-10), 0); A_est = A_est .* shrink; % 锚点约束：非负 + 列和为 1 A_est = max(A_est, 0); A_est = A_est ./ max(sum(A_est, 1), 1e-10); % 重建误差下降小于阈值则提前停止 cost = norm(X - E_est * A_est, 'fro'); if abs(prevCost - cost) / prevCost < 1e-4 break; end prevCost = cost; end end

这套代码把 CoNMF 的核心分解成了三层结构。第一层是端元更新，X * A_est'和E_est * (A_est * A_est')分别相当于最小二乘梯度的分子与分母，逐元素相除再乘上原 E，保证新 E 不会出现负数。第二层是丰度更新，同样采用乘性规则，随后对每一行计算 l2 范数，max(1 - lambda / rowNorm, 0)就是行软阈值操作：某行整体能量小于 lambda 时整行被压成 0，大于 lambda 时按比例缩小，但非零结构得以保留，这就是协同稀疏落地的关键。第三层是锚点投影，用列和归一化满足丰度之和为 1。

3.3 代码里的关键参数说明

sunsal 初始化时把 lambda 乘了 5，是因为它在第一轮要为后续迭代提供一个「偏稀疏、但结构稳定」的起点，比最终 CoNMF 收敛的稀疏度略高一些。PCA 初始化端元的逻辑是取前三个主成分得分绝对值最大的像元，这种方法虽然不如 VCA 严谨，但完全不需要额外工具箱，适合快速验证流程。若项目里已经装了高光谱专用工具箱，可以把这段替换成E_est = vca(X, K)，效果通常更好。

迭代停止条件用的是相对重建误差：当相邻两轮norm(X - E*A, 'fro')的相对变化低于 1e-4 就提前退出，避免固定轮数的盲目性。实际数据里如果 200 轮还没收敛，先检查lambda是否过大，它是导致收敛慢最频繁的原因。运行结束后建议立刻画出 A 的每行能量柱状图，能一眼看出哪些端元在整幅图里冗余。

4. CoNMF 参数调节与高光谱解混翻车现场

4.1 端元数 K：先估后用，别拍脑袋

K 是 CoNMF 里最敏感的参数。K 设小了，两种不同矿物被迫合并成一个端元，丰度图看起来干净但物理意义错位；K 设大了，算法会把噪声拆成一个端元，或者产生高度相似的冗余端元。经验上先用 HySime 或虚拟维度法估一个初始值，再结合场景知识微调。对一幅以矿物为主的 AVIRIS 影像，K 落在 10 到 20 之间比较常见；植被和城镇混合场景则要更大。下面这段用特征值阈值做粗略估计，不依赖额外工具箱：

% 基于噪声累计方差的粗估 [coef, score, latent] = pca(X'); noiseVar = median(latent) * 10; K_rough = sum(latent > noiseVar); K = min(K_rough, 30);

latent 是各主成分的方差，它按从大到小排列。真实端元对应的主成分方差明显大于噪声方差，所以设定一个相对噪声方差的倍数阈值，统计超过阈值的成分个数。这个估计偏保守，配合后续查看丰度行能量，可以快速锁定真实 K 的区间。

4.2 协同稀疏系数 lambda 与迭代轮数的配合

lambda 控制行稀疏的强弱。lambda 太小，约束不起作用，结果退化成普通 NMF；lambda 太大，丰度矩阵几乎所有行都归零，重建误差飙升。实际项目中我会先固定 K，在 1e-3 到 0.1 之间按对数间隔试三到五个值，画出重建误差和非零行数的折线图。一般选重建误差开始明显上升之前、非零行数刚好等于预期端元数的那个点。

迭代轮数方面，200 轮对大多数场景足够，但要注意与 lambda 的联动关系：lambda 越大，行软阈值越强，A 的更新越需要更多轮次重新分配剩余端元的流量。若发现 100 轮内重建误差曲线变成平线但丰度图仍有噪点，可以先提高迭代轮数到 400，若没有改善，再回头调 lambda。不要同时改多个参数，否则无法定位问题。下面是参数速查表：

4.3 高光谱解混常见的 4 个失败画面

第一个画面是丰度图出现整条全零行，原因基本是 lambda 过大，或者该端元在场景里确实不活跃。先调小 lambda 重跑，如果调小后依然全零，说明这个端元是初始化阶段留下来的冗余端元，应该退回上一层减小 K 而不是强行保留。

第二个画面是端元光谱出现负值或超范围值。CoNMF 的乘性更新虽然能维持非负，但若 E 的初始值里含异常像元，或数据没有归一化，端元值会越过物理上界。我习惯在每轮更新后加一行E_est = min(max(E_est, 0), 1)，保证端元对应的反射率或发射率保持在合理区间内。

第三个画面是不同端元光谱相关系数超过 0.99，通常出现在 K 估多了，或者影像中有两类物质光谱非常相似。此时先检查初始端元，若初始端元本身就相似，问题就不在 CoNMF，而在数据预处理：是否遗漏了水汽波段，是否在辐射定标前混入噪声通道。

第四个画面是丰度图出现明显的随机椒盐噪声。这种情况常见于训练时用了固定步长或迭代不充分，丰度没有完全收敛。把 maxIter 调大，同时检查 sunsal 初值是否给出了一组全零行，如果全零行过多，后续乘性更新无法将它重新激活，CoNMF 就会停留在局部最优。

5. 让 CoNMF 在实际高光谱场景里更准的三步收尾

5.1 先把 DN 值转成反射率再进 CoNMF

高光谱解混的物理前提是像元光谱符合线性混合模型，而线性混合严格成立的对象是反射率，不是传感器记录的数字量化值。如果数据提供方没有做大气校正，常见做法是用经验线性法或平场域归一化先把 DN 值转到反射率。经验线性法需要在场景内布置至少两个已知反射率靶标，用回归系数把 DN 映射到反射率；没有靶标时也可以用暗像元减去大气路径辐射项。公式表达的转换关系为：

R = (pi * L) / (Esun * cos(theta) * d^2)

其中 L 是辐射亮度，Esun 是波段太阳辐照度，theta 是太阳天顶角，d 是日地距离。辐射定标和大气校正做完了，CoNMF 提取的端元才能和实验室光谱库直接对比，解混才有跨场景迁移的可能。

5.2 用代表性像元把大规模场景切成小批量

CoNMF 的迭代里包含多次矩阵乘法，像元数 P 达到百万级别时内存和耗时都会失控。我一般先把影像做超像素分割，在每个超像素内部取最接近均值的像元作为代表，把整幅图的解混拆成「代表像元解混」和「全图丰度插值」两步。代表的端元用整幅图的代表像元计算出，保留全部像元，最后用带锚点约束的最小二乘反演把丰度映射回去。这样 K 的规模不变，但每次迭代的 P 降了一到两个数量级，运行时间比值约为初始计算量比值的一次方。

5.3 CoNMF 端元与高光谱 transformer 的协作思路

近几年高光谱影像分类常用到 transformer 结构，但纯数据驱动的 transformer 容易忽视物理约束。可以先用 CoNMF 解出的端元和丰度构造一个光谱先验特征，把它作为额外 token 拼进 transformer 的输入序列，让注意力机制同时看到原始光谱和物质组成信息。反过来，用 transformer 的分类注意力热力图筛出代表性像元再送给 CoNMF，也能显著减少不必要的像元输入。两者互补后，CoNMF 提供可解释的物质含量，transformer 提供全局上下文关系，是当前做高光谱定量反演时比较实用的组合。最后提醒一句：任何算法参数调完，都建议用光谱角匹配计算解混端元与真实光谱库的夹角，角度小于 0.1 弧度才说明端元质量可信。

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