简介面向高光谱遥感解混研究的一份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 0E 01^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 主循环里统一加锚点约束避免双重归一化造成丰度整体偏小。下面是三种方法的对比。方法稀疏作用的维度端元是否参与迭代典型输出sunsal逐像元列方向否E 固定单应稀疏丰度NMF 类算法无显式稀疏是端元丰度CoNMF协同行稀疏行方向是全局稀疏丰度3. 用 sunsal 初始化 CoNMFMATLAB 最小可复现脚本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×PP 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 太小约束不起作用结果退化成普通 NMFlambda 太大丰度矩阵几乎所有行都归零重建误差飙升。实际项目中我会先固定 K在 1e-3 到 0.1 之间按对数间隔试三到五个值画出重建误差和非零行数的折线图。一般选重建误差开始明显上升之前、非零行数刚好等于预期端元数的那个点。迭代轮数方面200 轮对大多数场景足够但要注意与 lambda 的联动关系lambda 越大行软阈值越强A 的更新越需要更多轮次重新分配剩余端元的流量。若发现 100 轮内重建误差曲线变成平线但丰度图仍有噪点可以先提高迭代轮数到 400若没有改善再回头调 lambda。不要同时改多个参数否则无法定位问题。下面是参数速查表参数常见范围过小的表现过大的表现建议调试顺序K 端元数830端元合并、丰度图缺类冗余端元、光谱高度相似先估 HySime再看行能量lambda1e-31e-1丰度行不稀疏大量丰度行全零对数网格搜索maxIter200400未收敛、端元偏模糊计算时间成倍增加看重建误差曲线是否平缓列归一化约束建议开启丰度比例没有物理意义部分像元丰度过小做结果验证时观察和是否为 14.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 弧度才说明端元质量可信。本文还有配套的精品资源点击获取
