简介连续投影算法SPA是一种光谱分析中常用的特征波长选择方法常与主成分分析结合实现高维光谱数据的降维。这份资料包装载了SPA算法的MATLAB实现、图形界面演示文件、验证与评估脚本以及使用指南面向从事光谱数据建模的科研人员和工程师可解决高维数据带来的过拟合与计算复杂性问题。压缩包共11个文件涵盖.m源码、.p加密程序、.doc/.ppt文档、.fig图形与.mat数据等类型大小约1.8MB各文件分别对应算法实现、交互演示、指标计算和理论讲解等不同用途。已有854人学习下载。通过阅读源码可掌握SPA的正交投影步骤借助GUI可直观观察特征选择过程配合验证脚本能评估模型性能配套读书报告则有助于理解SPA与PCA的联合应用适用于食品安全、环境监测等领域的分类任务。1. 连续投影算法不是黑匣子光谱选波长先搞懂它在选什么做近红外光谱回归时我见过不少同事把整条光谱几千个点直接丢进PLSR模型在训练集上漂亮同一批样品换一台仪器就崩。变量之间高度相关带来的共线性是光谱建模翻车的主要来源之一。连续投影算法Successive Projections Algorithm简称 SPA就是针对这个问题设计的波长选择方法。它不按“重要性”排名而是反复做正交投影保证最后选出来的变量彼此冗余最小、原始波长标签可迁移。想降低过拟合、做仪器选型或把模型搬上在线设备的人值得把它放在工具箱第一层。2. 为什么连续投影算法能选出“最有用”的波段投影逻辑与预处理前提2.1 共线性是光谱数据的老毛病SPA的解法是“剩余信息最大化”光谱数据里相邻波段的吸光度相关性经常超过0.99。原因并不复杂分子吸收带跨度远大于采样间隔同一个吸收带会被几十个采样点重复描述。回归模型想拿到稳定解必须先剥离这种重叠。SPA的思路非常直观把所有波段看成高维空间里的一组列向量从初始波长开始每一步都把尚未选中的波段向量投影到已选波长张成的子空间上。投影之后剩下的残差向量代表“还没被已选变量解释掉的信息”。残差向量越长说明这个波段与已选波段的线性相关性越弱于是选中它。循环k次就得到k个相互低相关的波段。有人会先计算每个波长与响应变量的相关系数取top-k这种做法在光谱选波长上不太够用。相关系数只衡量单个波长与y的线性关系没有考虑波长之间的冗余top10里很可能有9个落在同一个吸收带上。SPA相当于在“与y相关性”和“变量间冗余度”两个维度上做平衡这正是它比单纯排序可靠的地方。这里的计算量并不大。每一步投影的复杂度大约是候选波段数乘已选波长数几百个波段、k不超过30时普通台式机几秒钟就能跑完。它不是一个黑匣子只要把矩阵和投影关系写清楚任何一步都能在中间结果里检查。2.2 数据用什么形态进SPA均值中心化、导数与反射率转换SPA的正交投影基于向量内积内积对列向量的量纲非常敏感。同样一组数据不预处理和做过均值中心化选出的波长经常不一致。我会先把每一列减去列均值这是最基础的均值中心化。如果基线漂移明显先做一阶导或二阶导再中心化导数处理可以移除加性基线突出峰形变化。但要注意导数同时放大高频噪声信噪比较低的高光谱数据用过强预处理等于把噪声喂给SPA。确定用什么形态进SPA我一般会用一小批样本先跑一个SPA加PLSR的小实验比较原始吸光度、SNV、一阶导三种形态下交叉验证RMSE的差别。导数窗口宽度和中心化方式固定之后再正式跑全样本的变量选择。这样做的理由很实际预处理没有绝对正确的答案只有和当前仪器噪声水平匹配的答案。还有一道比预处理更基本的工序确认光谱到底是不是反射率。不少高光谱原始影像存的是DN值或辐射亮度。如果在DN上跑SPA选出的通常是仪器响应很强的波段而不是物质吸收特征。在ENVI工作流里一般先做辐射定标把DN转成辐亮度再用大气校正把辐亮度转成反射率。做完之后要检查曲线形状植被或土壤样本的反射率曲线在红边和水分吸收位置应该符合常识。我常拿ICVL高光谱数据集这类公开数据做快速验证。用loadmat读出来的结构通常包含影像数据块和波长向量先把三维影像reshape成“像素数乘波段数”的光谱矩阵再按反射率阈值去掉背景像素。这一步省略的后果是SPA把权重花在暗背景和目标区域的过渡带上对实际回归毫无贡献。3. 用Python复现连续投影算法从ICVL高光谱.mat到最小可用代码3.1 把高光谱图像拆成光谱矩阵loadmat、reshape和去背景ICVL高光谱数据集常见格式是.mat里面是一个三维影像块和对应的波长向量。第一步先把文件读出来确认字段名再动手。import numpy as np from scipy.io import loadmat # 读取ICVL高光谱数据集mat文件字段名以实际文件为准 mat loadmat(icvl_hsi.mat) print(mat.keys()) # 假设影像块字段是HData形状为(rows, cols, bands) img mat[HData] rows, cols, bands img.shape print(影像形状:, img.shape) # 把三维影像展平成“像素 x 波段”的光谱矩阵 X_all img.reshape(rows * cols, bands) # 用每个像素的平均值做阈值去掉暗背景像素 bright X_all.mean(axis1) mask bright 0.05 X X_all[mask] print(保留像素数:, X.shape[0])reshape成(行*列, 波段)是后面所有列运算的基础。SPA关注的是波段列之间的关系而不是像素的空间位置所以这一步展平没有任何信息损失。用每个像素在全部波段上的均值做背景筛选原理很简单暗背景像素在所有波段上响应都低。阈值0.05不能硬套最好先画一张bright直方图看目标和背景是否形成双峰再取峰谷位置作为阈值。提示如果mat文件是v7.3格式scipy的loadmat会直接报“不支持”这时候改用h5py打开。h5py读出来的数组维度顺序经常是倒着的需要转置后再用。3.2 SPA主循环逐层正交化、选残差最大波长SPA的最小实现并不复杂。核心是维护一组标准正交基让每次候选波长的投影计算都在同一套坐标系下进行。我写过的最小可用版本如下def spa_select(X, n_selected, init_idx0): 连续投影算法最小实现 X : 校正集光谱矩阵 (n_samples, n_bands)行样本列波长 n_selected : 要选的波长数 init_idx : 初始波长下标 Xc X - X.mean(axis0) # 均值中心化 n_bands Xc.shape[1] selected [init_idx] # 已选波长集合 basis [] # 已选波长张成空间的规范正交基 # 先把初始波长向量归一化放入正交基 v0 Xc[:, init_idx] n0 np.linalg.norm(v0) if n0 1e-12: raise ValueError(初始波长向量接近零检查预处理和去背景) basis.append(v0 / n0) while len(selected) n_selected: avail [i for i in range(n_bands) if i not in selected] best_idx None best_norm -1.0 for idx in avail: vec Xc[:, idx].copy() # 把候选向量中与已选波长线性相关的成分全部减掉 for b in basis: vec - np.dot(vec, b) * b norm np.linalg.norm(vec) if norm best_norm: best_norm norm best_idx idx selected.append(best_idx) # 把刚选中的波长也归一化后加入正交基 vec_new Xc[:, best_idx].copy() for b in basis: vec_new - np.dot(vec_new, b) * b n_new np.linalg.norm(vec_new) basis.append(vec_new / (n_new 1e-12)) return np.array(selected) # 在X上试选8个波长初始点用第0号波段 wavelength_idx spa_select(X, n_selected8, init_idx0) print(选中波长下标:, wavelength_idx)这个版本的关键在于维护basis列表。如果把每一步才做一次Gram-Schmidt、而不维护全局正交基后面选出的向量会和前面产生二次相关。每次迭代都让已选波长集合张成一个正交基候选波长的残差范数才是真正“没被解释掉”的信息量。vec - np.dot(vec, b) * b这一行把候选向量在基向量方向的分量减去循环完所有基向量后剩下的vec就是在正交补空间里的投影。n_selected不要一次给太大。SPA是贪婪算法越到后面加入的变量其投影残差越小边际价值越低。判断要不要继续加变量看的应该是下一章讲的交叉验证曲线。3.3 用PLSR交叉验证确定k肘部曲线怎么读SPA只负责选出波长位置选多少个波长由后续模型决定。常见做法是接PLSR因为选出的波段少PLSR可以处理变量间残留的轻微相关性。from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import KFold from sklearn.metrics import mean_squared_error kf KFold(n_splits5, shuffleTrue, random_state42) def cv_rmse(X_block, y, ncomp): rmse_list [] for tr, te in kf.split(X_block): pls PLSRegression(n_componentsncomp) pls.fit(X_block[tr], y[tr]) pred pls.predict(X_block[te]) rmse_list.append(mean_squared_error(y[te], pred, squaredFalse)) return np.mean(rmse_list) # 演示用y正式建模必须换成实测理化值 y X[:, 0] # k上限不超过样本量/3也不宜超过15~20 max_k min(15, X.shape[0] - 1) curve [] for k in range(1, max_k 1): sel spa_select(X, n_selectedk, init_idx0) ncomp min(3, k) # PLSR潜变量数不能超过波长数 rmse cv_rmse(X[:, sel], y, ncomp) curve.append(rmse) print(k , k, RMSE , round(rmse, 6))交叉验证误差曲线通常会先快速下降然后进入平台期偶尔在尾部略微反弹。正确的读取方式是找“肘部”也就是曲线斜率第一次明显变缓的位置。如果只看最低点很容易把噪声也选进模型。我一般会在肘部对应的k值附近多试两三个点再用独立测试集定最终k。4. 连续投影算法参数怎么设初始波长、k值上限与搜索策略4.1 初始波长不要盲选全扫描与载荷辅助的两种做法SPA的初始波长直接影响整个选择路径不同初始点可能收敛到完全不同的波长组合。通常做法有两个全扫描或者用PCA载荷辅助定位。全扫描就是把每个波段都作为一次初始点跑一遍SPA用交叉验证RMSE做比较。def scan_initial_wavelength(X, y, k): best_rmse np.inf best_ini 0 for ini in range(X.shape[1]): sel spa_select(X, k, init_idxini) err cv_rmse(X[:, sel], y, ncompmin(3, k)) if err best_rmse: best_rmse err best_ini ini return best_ini, best_rmse # 固定k8扫描全部初始波长 best_ini, best_rmse scan_initial_wavelength(X, y, 8) print(最佳初始波长下标:, best_ini, RMSE:, best_rmse)全扫描在数百个波段的高光谱数据上还能接受但到了上千波段每轮都跑交叉验证会非常慢。常见做法是先用PCA或PLSR的载荷图选出载荷绝对值最高的几个波段作为候选初始点只扫描这几十个位置再把最优位置附近几个点做细扫。这样计算量缩小一个数量级结果通常接近全扫描。初始波长的选择没有绝对标准经验上优先找吸收峰所在波段或载荷极值波段避免把初始点落在纯噪声区。4.2 k值上限与样本数的关系工程经验值k值上限我在实际项目中基本遵循三个约束不超过30不超过样本量的三分之一不超过波段数的十分之一。三个值取最小的一个。为什么是样本量的三分之一波长数一旦超过样本量的一半后面回归模型的自由度迅速消耗交叉验证曲线往往在尾部仍然下降但那是模型在拟合噪声不是真信号。ICVL这类高光谱数据展开后像素数很多样本量约束不明显但近红外建模经常只有几十个样本k取到20以上就要警惕。把k上限设小还有一个好处SPA是贪婪算法前面几步选出的往往是真正的强信号后面几步进入弱信号和噪声的模糊地带。与其让算法自己决定不如用交叉验证曲线在肘部附近截断。4.3 后悔药SPA粗选之后再用逐步回归收窄变量SPA不保证全局最优它是每步取局部最优的贪婪过程。如果发现选出的波长组合里有个别波长回归系数不稳定一个常见做法是在SPA选出的候选集上再做一轮逐步回归。把SPA选出的15个波长作为初始变量集向前引入或向后剔除按p值或AIC准则收窄到8~10个。这个过程不是SPA的替代而是把SPA当作粗筛器让逐步回归的搜索空间小到可控。我在实际工作中会把SPA和逐步回归组合使用先跑SPA得到15个波长再对波长做相关性矩阵检查把相关系数高于0.9的两个波长中回归系数更弱的那个去掉最后用PLSR验证。这样既保留SPA的原始物理解释性又规避局部最优问题。5. 连续投影算法避坑指南光谱选波长的五个真实翻车记录5.1 原始吸光度直接跑SPA选出的“重要波段”全是噪声峰现象SPA选出的波长集中在某个高频噪声区查看原始光谱曲线时发现这些位置有明显毛刺建模后测试集RMSE反而比全波段模型差。原因噪声列向量的方差并不小它们的投影残差范数看起来很大SPA把这些伪方差当成有效信息。没有做平滑和归一化时算法分不清吸收峰和仪器噪声。解决先做Savitzky-Golay平滑或小波去噪再做均值中心化或SNV。选完变量后把选定波长位置标记到原始光谱曲线上如果发现连续两个选点落在噪声毛刺区返回去调整平滑参数。这个检查步骤看起来普通但能拦住大部分无效选择。5.2 样本量太少SPA结果三次跑三次不一样现象一份只有40个样本的数据每次重新划分校正集和验证集SPA选出的波长组合都不同有时连初始波长都不一致。原因样本量小且光谱相似度高时候选波长的投影残差范数相差非常小随机抽样带来的微小波动足以改变下一轮最优路径。SPA本身对数据分布敏感样本量不足时它的搜索路径并不稳定。解决不要指望一次性得到稳定结果用重复采样的策略。把校正集随机划分30次每次跑一遍SPA加PLSR统计每个波长被选中的频率只保留出现频率高于0.7的波长。这个做法的代价是计算量多三十倍但对小样本光谱项目来说非常值。5.3 波长列单位混用投影计算完全失真现象数据合并时一部分文件用纳米记录波长另一部分用微米记录SPA选出的下标看似连续实际对应物理波长位置错乱模型迁移到新仪器时完全失效。原因索引基于列位置如果单位不统一相同下标在不同批次里代表的光谱位置完全不同。SPA对列位置运算而下游应用依赖物理波长二者一旦脱节就会出问题。解决加载数据后立刻把波长统一到同一单位建议统一为纳米。保存结果时同时输出“下标、波长、单位”三列遇到新数据先验证波长向量单调递增再做任何运算。5.4 k值只看交叉验证最低点独立测试直接过拟合现象交叉验证曲线在k23处RMSE最低选23个波长之后独立测试集表现反而比k12差。原因SPA只保证选出的变量彼此低相关不能保证后几个变量不与噪声相关。k足够大时后几个变量开始充当残差拟合的补丁交叉验证曲线在尾部仍然下降但这是过拟合信号。解决把k上限先按样本量的三分之一收紧再把数据切成训练、验证、测试三份。训练加验证集负责选k测试集只允许用一次。如果曲线在达到上限前一直持续下降说明信号强度不够或预处理不合理需要检查而不是加k。5.5 高光谱转反射率不干净选出的波段在ENVI里对不上矿化异常现象实验室漫反射光谱上SPA选出2200nm附近的波长但遥感影像对应位置的反射率曲线有明显异常后续在ENVI里做矿化蚀变信息提取时没有突出目标异常。原因实验室光谱经过白板校正遥感数据没有经过完整的大气校正光谱基准不一致。SPA在实验室尺度上选出的波长直接套到遥感尺度波段中心对不上吸收特征被大气窗口残留噪声盖住。解决遥感流程里先做辐射定标和大气校正把数据转成地表反射率再做光谱重采样让实验室光谱和影像光谱的波段中心、带宽一致。SPA选出的特征波长先映射到距离最近的影像波段再进入波段运算或光谱角制图。6. 进阶验证SPA选出的波长怎样才算真正落地可用6.1 选前选后模型对比值不值得投入一表看清对比项全波段模型SPA选波段模型输入变量数全波段可能数百到上千通常5到15个训练RMSE通常更低略高但不代表泛化差交叉验证RMSE视共线性程度而定与全波段接近方差更小独立测试RMSE变量多时容易过拟合更稳定硬迁移成本需要完整光谱采集模块可对应窄带滤光片或特定通道我会把这份对比表当成项目决策依据。如果两条曲线的独立测试RMSE差别在5%以内SPA的主要价值不是提高精度而是降低成本、提高稳定性。如果全波段模型测试集明显变差而SPA模型保持稳定那就是选择了正确的方向。6.2 把光谱波长映射到硬件通道带宽与中心波长检查SPA给出的中心波长只是一个数值落地到硬件前还要检查带宽。窄带滤光片的半高宽如果远大于光谱仪采样间隔会把SPA选出的窄吸收峰抹平。我一般会把选中波长相邻±半高宽区间的光谱取平均再做一次模型验证。如果RMSE比直接用单点波长上升超过10%说明中心波长选择没问题但硬件通道的带宽不匹配需要换滤光片或调整k值。6.3 从实验室走向GF-5高光谱蚀变信息提取与ENVI工作流的衔接GF-5高光谱影像覆盖可见光到短波红外常规蚀变提取流程是在ENVI里做辐射定标和大气校正得到地表反射率再用光谱角制图或波段运算圈定蚀变异常。SPA在这个流程里负责压缩冗余波段减少噪声传播。把实验室光谱的SPA选择结果映射到GF-5波段时重点检查2200nm附近的吸收带是否落在有效波段上选出的波长必须与目标矿物的诊断性吸收特征匹配。我做这类工作前总会先跑一次SPA再用全波段模型对照而不是直接信任筛选结果。这个习惯帮我避开了很多重复劳动希望帮到你。本文还有配套的精品资源点击获取