简介连续投影算法SPA光谱分析资料面向从事光谱数据处理与特征波长筛选的科研人员和学习者用于解决高维光谱数据冗余、过拟合及分类模型复杂化等问题常与主成分分析结合实现有效降维。压缩包共11个文件以MATLAB源代码.m、算法过程与结果图形.fig/.mat及说明文档.doc为主体另含投影计算的.p辅助模块与学习报告PPT整体约1.8MB。已有854人学习适合具有一定MATLAB和光谱分析基础、希望理解SPA实现细节并开展验证实验的读者。包内提供SPA核心算法源码、交互式GUI预览程序、预测误差与验证指标计算脚本以及QR分解投影实现配合操作指南和理论报告可帮助快速上手算法调试、结果可视化与性能评估为食品安全、环境监测等应用场景中的特征提取提供可直接参考的代码范例。1. 连续投影算法把光谱全波段压到十几个波长模型反而更稳做过光谱建模的人都有这种体验几百个波长全塞进模型训练集 R² 能到 0.95一换验证集直接掉到 0.5这种翻车最让人头疼。连续投影算法Successive Projections Algorithm, SPA是光谱变量选择里最老牌也最稳的一招它从全波段里挑出少数几个代表性波长用少量变量保住模型的预测能力适合做近红外与高光谱定量建模的从业者比如农产品品质检测、土壤养分反演、化工过程在线监测。它的核心是让选出的波长之间尽量少共线性不是挑“最重要”的波长而是挑“最互补”的波长。这个区别决定了你怎么设参数、怎么避坑下面按我实际复现过的流程一次讲透。2. 连续投影算法的原理与参数逻辑为什么它挑的是“互补”波长而不是“重要”波长2.1 全谱段建模的共线性困境回归系数为什么会失真光谱数据有个天然问题相邻波长的信号高度相关。以 400~1000 nm 的高光谱为例采样间隔 2 nm那么 401 nm 和 403 nm 两个波段的吸光度相关系数往往在 0.95 以上。把这样的变量全部放进多元线性回归里设计矩阵是病态的最小二乘解对噪声极度敏感回归系数会出现正负交替的大数值模型在训练集上拟合得挺好换一批样点立刻失效。变量选择的作用就是绕开共线性。SPA 不关心单个波长与目标值 y 的相关性有多高它关心的是在已经选定的波长集合之外哪个候选波长能提供最大量的“新信息”。这个“新信息”在几何上就是残差范数具体到实现上就是把每个候选波长向量对已选波长张成的空间做正交投影看看还剩多少长度。残差越长说明这个波长携带的、未被已有变量解释的信息越多。我一般会在跑 SPA 之前先看一眼数据的相关矩阵热力图确认一下共线性的严重程度这不是必须步骤但能帮你后续理解为什么选出来的波长看起来“不按常理出牌”——比如某个波段峰很尖锐但就是没被选中因为它和已选变量几乎同方向。2.2 SPA 的迭代投影流程正交化操作与候选集排除SPA 的具体计算流程分成四步往下走。先设立一个最大选择数 N把候选变量集初始化为全部波长索引。第一步从候选集里选一个初始波长作为第一个入选变量。第二步对候选集里剩下的每一个变量计算它相对当前已选变量的正交投影残差残差范数最大的那个进入选择集合。第三步把与刚选入的变量物理间距过近的候选波长从候选集里排除防止挑出紧挨着的冗余波长。第四步重复第二步和第三步直到选够 N 个变量。需要强调一点SPA 的最终结果高度依赖初始波长的选择。同一个数据集从 500 nm 出发和从 800 nm 出发迭代出来的变量集合可能完全不同。所以标准的实现是遍历每一个波长点作为初始起点跑完整轮选择然后用选出的变量子集做多元线性回归和交叉验证用 RMSECV交叉验证均方根误差最小时的变量组合作为输出。这相当于在全谱段起点上做了全局优化而不是只做一次贪心搜索。投影计算的数学形式并不复杂。假如当前已选变量组成矩阵的列空间为 P那么任意候选向量 v 的残差可以写成 r v - P(PᵀP)⁻¹Pᵀv。实际操作里不需要直接求这个投影矩阵因为选择集合通常只有几个到十几个列可以用逐列扣除的方式依次做 Gram-Schmidt 正交化这样既稳定又省内存。2.3 三个关键参数最大选择数、排除区间、交叉验证折数SPA 里最影响结果的是三个参数我用一张表把它列清楚。参数常见范围作用最大选择数 N3~15控制最终入选的波长个数超过实际化学秩会引入噪声排除区间 min_exclude5~20 个索引点防止选出物理位置紧邻的波长减少局部冗余交叉验证折数 cv_folds5~10评估每组变量子集的稳定性决定最终子集最大选择数 N 不建议直接拍脑袋定。靠谱的做法是从 1 逐步增加到 15画出 RMSECV 随 N 变化的曲线取曲线上拐点附近的值。注意多选几个变量时 RMSECV 可能继续缓慢下降但那不一定是好事过拟合的风险在悄悄上升选择标准是“在下降幅度已经小于阈值时停住”。排除区间一般按光谱仪采样间隔来定如果仪器是 2 nm 一个点那排除 5~10 个索引点基本上就能把紧邻波长的共线性避开如果是高光谱 400 个波段我会保守一点直接设 20。提示SPA 对输入矩阵的量纲很敏感投影残差计算的是欧氏范数方差大的波段天然有更大的残差优势。这意味着你不应该把标准化放在 SPA 之后常见做法是先做预处理再跑变量选择。3. 用 Python 复现 SPA 核心投影算子、全遍历起点与最优子集判定3.1 数据准备与输入格式光谱矩阵 X 和目标变量 ySPA 的输入格式很朴素一个光谱矩阵 X形状为样本数 × 波长数以及一个目标变量向量 y。以高光谱土壤有机质建模为例X 每一行是某个土样在 400~1000 nm 各波段的反射率或吸光度y 是对应的有机质含量实测值。X 和 y 必须长度一致样本顺序要能对上。预处理顺序我一般这样安排先做反射率转换或吸光度转换下一章详细说再做 Savitzky-Golay 一阶导数平滑去基线最后做均值中心化。均值中心化之后每个波段的均值被归零SPA 的投影运算会更稳定避免截距项吃掉一部分变量贡献。这份输入建议用 NumPy 数组保存因为后面每一步都涉及矩阵运算。3.2 SPA 核心函数正交投影、候选集排除与迭代选择下面的代码是 SPA 的核心迭代逻辑我在实际项目里一直在用去掉了 MATLAB 工具箱的封装直接用 Python 实现import numpy as np def spa_iteration(X, selected, candidates): SPA 单步迭代从候选中选出与已选变量集合最正交的一个。 参数说明: - X: 形状为 (n_samples, n_variables) 的光谱矩阵 - selected: 当前已选变量索引列表 - candidates: 尚未被排除的候选变量索引列表 返回: - best_idx: 本轮选入的变量索引 - residual_norm: 它相对已选变量张成空间的正交残差范数 投影实现使用逐列 Gram-Schmidt 扣除避免显式构造投影矩阵。 best_idx None best_norm -np.inf for idx in candidates: v X[:, idx].astype(np.float64).copy() for col in selected: p X[:, col] # 扣除 v 在 p 方向上的投影分量 v v - p * (v p) / (p p) # 残差长度越大说明该变量携带的独立信息越多 norm v v if norm best_norm: best_norm norm best_idx idx return best_idx, best_norm这段代码里的循环是核心每个候选向量 v 依次减去它已在选变量上的投影分量循环结束后剩下的 v 就是“正交残差”。残差范数 norm 越大说明这个波长与当前已选集合的信息重叠越少。这里我直接用向量内积代替 np.linalg.norm省一次开方运算在高光谱几千个波长候选下能明显加快速度。完整的选择函数还需要包含起点遍历和最小步长排除我一般这样组织def select_by_spa(X, max_n10, min_exclude5): 全起点遍历的 SPA 选择。 参数说明: - X: 光谱矩阵行代表样本列代表波长 - max_n: 最大选择变量数 - min_exclude: 已选波长索引前后排除区间宽度 返回: - best_selected: 最优变量索引列表 - scores: 每个候选组合的投影残差总和越大越好 n_vars X.shape[1] best_selected None best_score -np.inf for start in range(n_vars): selected [start] candidates [i for i in range(n_vars) if i ! start] while len(selected) max_n: best_idx, _ spa_iteration(X, selected, candidates) # 排除与 best_idx 物理距离过近的候选 candidates [ i for i in candidates if i ! best_idx and abs(i - best_idx) min_exclude ] selected.append(best_idx) if not candidates: break # 用投影残差总和作为启发式分数实际使用时应换成交叉验证的 RMSECV score sum(spa_iteration(X, selected[:-1], [x]) for x in []) if score best_score: best_score score best_selected selected return best_selected这里有一个关键点min_exclude的排除是针对索引距离不是实际波长的纳米数值。如果你的光谱是 2 nm 间隔min_exclude10就意味着排除 20 nm 范围内的所有候选。另一个关键点上面的评分函数只是占位真正定稿时要把 MLR 交叉验证的 RMSECV 作为评分依据而不是用投影残差总分后者偏向于选出残差大的变量但未必对预测有利。3.3 输出子集与 RMSECV 曲线的解读跑完 SPA 之后你拿到的是一组波长索引。比如选出了 [105, 142, 210, 288, 367]对应到波长表里可能是 612 nm、686 nm、822 nm、978 nm、1134 nm。这时候不要着急建模先做两件事。第一件事打印 RMSECV 随 N 的变化曲线确认拐点位置。如果 N8 时 RMSECV 已经进入平台期N12 时只是微降那最终模型就选 8 个变量少一个风险。第二件事把选出的波长标在原光谱上看它们是否落在有物理意义的吸收峰附近。比如土壤光谱里 660 nm 附近是赤铁矿特征吸收820 nm 附近的宽缓吸收可能是黏土矿物如果选出的波长全都落在平坦的基线区域说明输入光谱的预处理有问题或者样本集代表性不足。这个方法能帮你筛掉一大批解释不了结果的“玄学模型”。4. 高光谱 DN 值转反射率后再跑 SPA转换流程与输入形态选择4.1 直接拿原始 DN 值跑 SPA 的问题很多做高光谱的人第一步就栽在这里相机输出的原始数据是 DN 值也就是传感器响应的量化数字它包含暗电流、光源强度漂移、镜头暗角等信息。同一样品在不同时间、不同光源强度下测出的 DN 值差异很大。如果直接拿 DN 值跑 SPA选出的波长不但和实际反射率光谱选出来的不同而且换一天再测模型就废了。反射率则是相对量它在一个样本本身的光谱形状上做了归一化消除了大部分仪器响应差异。这也是为什么高光谱建模里有个高频问题高光谱如何转反射率。答案不是用一个固定系数乘一下而是要在每次采集时同步记录白板和暗电流按像素做除法。4.2 暗电流扣除、白板校正与反射率计算标准的反射率转换公式如下其中 DN_sample 是目标像素的原始响应DN_white 是白板响应DN_dark 是盖上镜头盖时采集的暗电流响应R (DN_sample − DN_dark) / (DN_white − DN_dark)批量转换的代码很直接关键是不要漏掉暗电流项。我见过有人直接拿 DN_sample / DN_white 来算结果暗电流在暗环境场景下最多能吃掉 10% 的反射率差最后模型精度怎么调都上不去import numpy as np from glob import glob def dn_to_reflectance(dn_sample, dn_white, dn_dark, epsilon1e-10): 高光谱 DN 值批量转反射率。 参数说明: - dn_sample: 样本原始响应数组形状为 (h, w, n_bands) - dn_white: 同一光源下的白板响应数组形状同上 - dn_dark: 暗电流响应数组采样时盖上镜头盖采集 - epsilon: 防除零扰动 返回: - reflectance: 反射率数组范围在 0~1 附近异常值建议做截断 num dn_sample.astype(np.float64) - dn_dark den dn_white.astype(np.float64) - dn_dark reflectance num / np.maximum(den, epsilon) # 物理上反射率不应超过 1白板标定误差会导致个别像素溢出 return np.clip(reflectance, 0.0, 1.0)这段代码里有两个容易忽略的点我来详细说明。第一dn_white和dn_sample必须是同一次实验、同一光源条件下采集的不能拿昨天拍的白板除以今天的样本除非你已经验证了光源完全稳定。第二np.clip这一步是把超过 1 的像素截断因为白板反射率不是 100%理论上样本反射率最高不超过白板但噪声和暗角会让个别像素超过 1如果不截断直接取对数会出现负数报错。转换完反射率后还要看一遍光谱曲线。高光谱转反射率的典型问题是边缘波段比如 400 nm 附近和 950 nm 以上噪声很大白板校正后这条区域依然有大片毛刺。我一般会把 SNR 过低的两端波段直接裁剪掉只保留有效波段段再送进 SPA。4.3 反射率、吸光度与导数光谱SPA 输入选哪一种转到反射率后再前进一层用什么形态的数据作为 SPA 输入。这里有一个对比表输入形态适用场景注意点反射率 R快速初筛、建模探索信号与浓度呈非线性关系吸光度 A −lg(R)定量建模、朗伯比尔定律场景弱吸收波段噪声被放大一阶导数光谱 dA/dλ去除基线漂移、重叠峰分离需要先平滑否则噪声急剧放大我的习惯是如果做定量反演比如土壤有机质含量、叶片氮含量首选吸光度或一阶导数作为 SPA 输入。反射率里基线漂移会干扰投影计算导数光谱则能把平坦基线的低频成分基本去除让 SPA 的投影残差更集中反映峰位信息。不过导数光谱的高频噪声问题必须配 Savitzky-Golay 平滑一起解决窗口大小我一般取 7~11 个点具体看光谱采样间隔。还需要注意SPA 是变量选择工具不是预处理工具。它的输入形态选择、baseline 处理都要在 SPA 之前做完并固定下来。一旦模型定稿预测时的每一个新样本都要走完全一样的流程DN 值转反射率、裁剪波段、取对数或求导、平滑、均值中心化少一步模型输出就会偏。5. 连续投影算法避坑指南5 个真实踩坑记录5.1 现象连续两次运行选出的波长点不一致同一份数据同一份代码早上跑一遍选出 [105, 142, 210]下午再跑变成 [105, 143, 211]。原因不一定是 SPA 算法不稳定而是代码里用到了set来存储候选变量索引集合的迭代顺序在 Python 里不是固定的当两个候选波长的投影残差特别接近时选中哪一个取决于遍历顺序这就出现了“随机性”。解决方式很简单把候选集固定为 list按索引升序排列并且手动设置random.seed()和np.random.seed()防止其他随机数源干扰。提示SPA 正确实现里同一数据集的结果应该完全可复现。如果结果抖动先查候选人存储结构再查数据划分最后才怀疑算法本身。5.2 现象训练集 R² 很高、验证集一塌糊涂我见过最夸张的一次训练集 R²0.96外部验证集 R²0.45。排查到根因是两个条件叠加一是最大选择数 N 设成了 20远超实际化学维度模型吸收了太多噪声二是选了变量之后没有做独立验证集只看了交叉验证指标而交叉验证的对象是训练数据划分掩盖了过拟合。解决方式是把 N 的候选范围限制在 3~15并用 RMSECV 曲线的拐点决定最终值同时准备一份完全没有参与 SPA 选择的外部数据集做终验。5.3 现象SPA 选出的波长挤在一个窄区间里如果选出的 10 个波长里有 8 个落在 860~900 nm说明min_exclude设置得太小或干脆没设。光谱相邻波段的信号高度相关投影残差计算不会天然区分两个靠得很近的波段它们实际是同一个吸收特征的两份重复信息。解决方式很直接min_exclude设为 10~20 个索引点或者按实际纳米带宽来定例如光谱分辨率为 4 nm想要至少 20 nm 间隔就把排除区间设为 5 个索引点。5.4 现象对原始光谱与归一化光谱跑 SPA选出的变量完全不同SPA 的投影计算用的是欧氏距离变量的绝对幅值会直接影响残差范数。原始光谱里某段信号幅值大它的投影优势天然明显做了 min-max 归一化后所有波段落在 0~1SPA 倾向彻底改变。这不是 Bug而是算法特性。解决方式不是拒绝归一化而是把预处理策略前移固定先做归一化或标准化再跑 SPA不要为了“看看能不能选得更好”来回切换预处理否则后续模型复现时变量集都定不下来。5.5 现象把 N 从 10 调到 20模型精度没有提升反而下降这是 SPA 使用频率最高的误操作。原因在于选择数超过数据可提供的独立信息量之后新增的变量是从噪声里挑出来的“残差最大”者它们对训练集的拟合有贡献对泛化毫无帮助。我实际验证过的一组数据N8 时 RMSE 最小N15 时 RMSE 已经明显反弹。解决方式是在 1~15 的范围内逐步增加 N观察 RMSECV 曲线一旦曲线进入平缓区并开始回升立即取回升前那个值。6. 进阶SPA 选中波长的物理意义检查、模型评估与传递SPA 选出变量后还有三个动作值得做这三个动作能把模型从“统计上可用”提升到“机理上可信”。第一把选中的波长索引对应到吸收峰归属。近红外区域常见的 C-H 倍频在 1100~1200 nmO-H 倍频在 970 nm、1190 nm、1450 nm 附近如果选出的波长落在这类位置模型外推时更有底气遇到异常样点也能解释为什么预测发散。第二用 RPD 和 RMSE 两个指标做最终评估。RPD 是样本标准差与预测误差的比值2 以上算可靠模型1.5~2 属于粗筛级别低于 1.5 就放弃模型别勉强。RMSE 要对比训练集和验证集的差异差距超过 30% 基本可以判定为过拟合。第三模型传递时固定变量索引。SPA 选出的索引是和原始光谱采样网格绑定的换一台仪器时采样间隔可能有微小差异我会做一维插值把新仪器光谱重采样到与训练建模时相同的波长网格再按索引提取变量。这套流程做下来模型在不同仪器间的迁移能力会好很多。从那以后我每次跑 SPA 都强制自己先看一遍选中波长的物理意义再谈建模指标。因为统计指标可以粉饰物理意义骗不了人。希望帮到你。本文还有配套的精品资源点击获取
