企业官网建设流程全解析

简介：连续投影算法（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(axis=1) 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_idx=0): """ 连续投影算法最小实现 X : 校正集光谱矩阵 (n_samples, n_bands)，行样本，列波长 n_selected : 要选的波长数 init_idx : 初始波长下标 """ Xc = X - X.mean(axis=0) # 均值中心化 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_selected=8, init_idx=0) 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_splits=5, shuffle=True, random_state=42) def cv_rmse(X_block, y, ncomp): rmse_list = [] for tr, te in kf.split(X_block): pls = PLSRegression(n_components=ncomp) pls.fit(X_block[tr], y[tr]) pred = pls.predict(X_block[te]) rmse_list.append(mean_squared_error(y[te], pred, squared=False)) 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_selected=k, init_idx=0) 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_idx=ini) err = cv_rmse(X[:, sel], y, ncomp=min(3, k)) if err < best_rmse: best_rmse = err best_ini = ini return best_ini, best_rmse # 固定k=8，扫描全部初始波长 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值只看交叉验证最低点，独立测试直接过拟合

现象：交叉验证曲线在k=23处RMSE最低，选23个波长之后独立测试集表现反而比k=12差。

原因：SPA只保证选出的变量彼此低相关，不能保证后几个变量不与噪声相关。k足够大时，后几个变量开始充当残差拟合的补丁，交叉验证曲线在尾部仍然下降，但这是过拟合信号。

解决：把k上限先按样本量的三分之一收紧，再把数据切成训练、验证、测试三份。训练加验证集负责选k，测试集只允许用一次。如果曲线在达到上限前一直持续下降，说明信号强度不够或预处理不合理，需要检查而不是加k。

5.5 高光谱转反射率不干净，选出的波段在ENVI里对不上矿化异常

现象：实验室漫反射光谱上SPA选出2200nm附近的波长，但遥感影像对应位置的反射率曲线有明显异常，后续在ENVI里做矿化蚀变信息提取时没有突出目标异常。

原因：实验室光谱经过白板校正，遥感数据没有经过完整的大气校正，光谱基准不一致。SPA在实验室尺度上选出的波长直接套到遥感尺度，波段中心对不上，吸收特征被大气窗口残留噪声盖住。

解决：遥感流程里先做辐射定标和大气校正，把数据转成地表反射率，再做光谱重采样，让实验室光谱和影像光谱的波段中心、带宽一致。SPA选出的特征波长先映射到距离最近的影像波段，再进入波段运算或光谱角制图。

6. 进阶验证：SPA选出的波长，怎样才算真正落地可用

6.1 选前选后模型对比：值不值得投入，一表看清

我会把这份对比表当成项目决策依据。如果两条曲线的独立测试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再用全波段模型对照，而不是直接信任筛选结果。这个习惯帮我避开了很多重复劳动，希望帮到你。

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