详解：光谱变量筛选原理与Python实战)
简介本资源是一份面向数据科学初学者与光谱分析从业者的SPA连续投影算法原理精讲与工程实践教程聚焦高维数据降维这一核心问题特别适用于生物信息学、化学计量学及机器学习预处理场景。压缩包共14个文件6.55MB含11个MATLAB源码.m——涵盖spa_train、spa、validation等主流程与子函数2个实测光谱数据集.mat包括jasperRidge2_R198.mat等典型高维光谱样本1份结构清晰的中文教程文档.docx系统梳理算法原理、迭代逻辑、调试要点与光谱应用案例。已有2413人学习下载读者可直接运行调试后的完整代码链复现从数据加载、投影迭代、方差评估到分类预测的全流程并获得作者在内存管理与收敛性优化方面的实操经验显著降低算法落地门槛。 做光谱建模的人十有八九都会碰到同一个问题波长变量太多样本数又不够模型老是过拟合今天调出来的效果明天就翻车。我早年做近红外定量分析的时候一套光谱数据经常是上千个波长点但实际有效信息可能就集中在几十个波段里。那时候最苦恼的就是怎么从一堆冗余变量里挑出真正有用的那部分直到后来接触了连续投影算法Successive Projections Algorithm简称SPA才算是把这条弯路走顺了。SPA的核心工作很简单它能在高维光谱数据中自动挑出一组共线性最小、信息冗余最低的特征波长再用这些筛选出来的变量去建校正模型。它不依赖任何特定模型框架和偏最小二乘PLS、多元线性回归MLR都能搭配使用很适合做光谱定量分析的预先降维步骤。这篇教程我会从算法原理一步步拆开讲配上完整的Python实现和实际案例最后聊聊我在工程里踩过的坑。无论是刚入门化学计量学的学生还是已经在做光谱建模的工程师这篇文章都应该能帮到你。1. 为什么需要变量选择光谱数据里的共线性问题1.1 光谱建模的痛点高维、共线性、过拟合先说说我们面对的原始数据长什么样。以近红外光谱为例一台光谱仪扫描一个样本可能在4000到10000波数范围内采集上千个数据点每个数据点就是一个变量。如果我有200个样本那数据矩阵就是200乘1000的维度。用这样的矩阵去建常规的多元线性回归模型数学上是不可行的因为变量数远远大于样本数即使强行用偏最小二乘去压缩潜变量也会面临信息分散、模型解释困难的问题。更麻烦的是光谱数据里相邻波长的响应不是独立的而是高度相关的。某个官能团的吸收峰往往横跨几十个波数对应几十个变量之间都有极强的关联。这种共线性会让很多回归算法不稳定你今天采样环境稍微变一点模型系数就可能发生剧烈波动。说白了变量之间大量重复的信息不仅没帮上忙反而成了噪声放大器。这个问题的本质是我们需要在保证预测能力的前提下找到一组信息互补、冗余最小的波长变量。这好比你请一个团队干活如果所有人都是同一种能力结构团队整体战斗力反而很弱真正理想的团队是每个人都能填补其他人的盲区。SPA想要做的就是这件事找出那些彼此最“不重复”的波长变量。1.2 变量选择算法选型对比SPA、CARS、UVE、VIP光谱领域里变量筛选算法不少我自己先后试过VIP变量重要性投影、UVE无信息变量消除、CARS竞争性自适应重加权采样还有SPA。每种方法的思路差别挺大适用场景也不一样。先看VIP它是在PLS模型的框架下根据每个变量对模型解释方差的贡献度来打分的好处是计算快、直观缺点是它本质上只衡量了“重要性”并没有考虑变量之间的共线性结构。UVE的做法是把变量和一组随机噪声变量放在一起比较稳定性剔除那些稳定性低于噪声的变量这个方法能有效去噪但有时候也会误删一些有用但贡献不强的小信号。CARS是模拟达尔文进化论的思想通过自适应重加权采样不断淘汰“劣质”变量收敛快但因为有随机抽样环节每次运行结果会有一定波动需要跑多次取共识。SPA和上面几个都不太一样。它不是从统计显著性出发而是从几何角度衡量变量之间的关系。它逐步选出与当前已选变量最少冗余的波长每一步都保证新选的变量和前面选出的变量在向量空间上尽量正交从根源上砍掉共线性。这种思路在光谱数据这种共线性很强的场景里效果特别明显。1.3 连续投影算法的核心思路一句大白话SPA到底在做什么我用一句最直白的话概括它每次从剩下的变量里挑一个让这个变量的光谱响应向量与所有已选变量向量之间的“重叠程度”最低也就是在向量空间中对应的投影距离最大然后把信息最不重复的那个波长收进来。这样一轮一轮挑下去最终得到一组彼此“各自独立”的波长组合。这里的“投影”就是一个向量在另一个向量上的正交投影长度。如果两个变量谱图长得几乎一模一样那么一个变量在另一个变量上的投影就会很长两者冗余程度就很高。反过来如果两个变量谱图差异很大投影长度就越小。SPA的策略就是每一轮都选投影距离最大的那个变量让新选中的变量和已有变量整体冗余度最低。2. SPA的数学原理与步骤拆解2.1 投影计算与共线性度量从向量投影说起要动手实现SPA第一步得懂向量投影运算。假设我们有两个波长变量对应的响应向量分别记作向量x和z每个向量都是样本维度的列向量。我们把x投影到z上的投影向量记作proj_z(x)计算公式是proj_z(x) (x · z / (z · z)) · z这里x · z是内积z · z是z自身的内积。投影向量是z方向上的一个向量它表示x中能被z“解释”掉的部分。如果x和z完全一样完全线性相关那么投影向量就等于x自身也就是说x的信息完全冗余于z如果x和z垂直则投影向量为零向量说明两者完全没有冗余。SPA在筛选变量时对于当前选中向量z它计算每个候选变量x在z上的投影向量长度即P(x) || x - proj_z(x) ||这个值表示x中与z不重叠的信息量。P(x)越大说明x携带的独立信息越多。SPA每一轮就选P(x)最大的那个变量作为下一个入选波长。这也解释了为什么它天然适合抵抗共线性两个相关性极高的变量之间P(x)会非常小几乎没机会同时被选中。2.2 完整算法流程初始化、迭代投影、候选变量集合完整SPA算法的设计比上面描述的要周详一些它并不是只指定一个起始变量然后一路选到底而是从不同的初始变量出发生成多组候选变量子集再挑选整体效果最好的那一组。我将SPA的完整流程拆解成四个阶段阶段一确定初始变量。在全部变量中分别以每一个变量作为初始入选变量然后分别执行后面的投影迭代。也就是说如果原始数据有1000个波长那么会有1000轮独立的搜索每一轮都从不同的初始波长开始。阶段二循环投影选择。在每一轮独立的搜索里设置一个预设的最大选变量数N一般取10到30之间或者根据经验公式定然后从初始变量出发迭代执行N-1次投影计算。每次迭代时计算所有尚未入选的变量与最近一次入选变量之间的投影距离选出投影距离最大的变量作为新的入选变量。阶段三形成候选子集。每一轮搜索结束都会得到一组长度等于N的波长索引序列。把所有初始变量跑完之后我们手上就有了与变量总数相同数量的候选子集每个子集包含N个波长变量。阶段四用校正模型评估。对每一个候选子集使用选出的N个波长变量建立多元线性回归MLR模型计算验证集的均方根误差RMSE选验证误差最小的那组波长作为最终输出。这种“多初始点搜索模型评估”的双层结构保证了SPA不会因为初始变量选得太差而错过全局最优子集。虽然计算量有所增加但对几百个变量的光谱数据来说完全在可接受范围内。2.3 为什么选“投影最大”而不是“相关最小”一个容易混淆的地方很多初学者会在这里卡住。他们觉得SPA要选“共线性最小”的变量那直接找相关系数最小的两个波长不就行了为什么非要绕一圈算投影、选投影距离最大我当年也绕了很久才想明白这中间的差别。两变量之间的相关系数低只说明它们之间没有明显的线性趋势但并不能保证它们在多元回归中不互相干扰。而SPA要的是变量对整体信息空间的“增量贡献”。它关心的是在已经选了某些变量的前提下再增加哪一个变量能带来最多的新信息。投影距离衡量的恰恰是“新增信息量”。比如变量A和变量B的相关系数是0.01看起来很独立但是B其实只是A经过一个微小变换得到的那么B在A方向上的投影可能仍然很大意味着B并没有带来多少新东西。相关系数小并不代表冗余度低这种微妙的差距正是SPA选择投影距离作为准则的原因。2.4 与MLR结合的SPA-MLR变量数怎么定SPA输出波长集合后还需要搭配一个校正模型来形成完整预测方案最常见的搭配是MLR也就是标准的最小二乘多元线性回归。因为经过SPA筛选后的变量数量已经从上千降到十几个MLR在数学上变得可行而且模型解释性比PLS强很多直接能看到每个波长对目标的贡献系数。变量数N怎么定这是一个很关键的问题。定得太少信息不足模型拟合差定得太高筛选变量过多MLR可能重新陷入过拟合。一个实用的经验值是N不超过样本数的1/10比如你有200个样本N最大选到20。实际操作中一般设置一个N的上限然后在1到这个上限之间逐一尝试观察交叉验证RMSECV的变化。通常RMSECV会先随着变量数增加而下降到某个点之后开始回升或趋于平稳选择回升前的那个拐点对应的N即可。3. Python完整实现与实操案例3.1 环境准备与数据说明动手之前我们先准备好环境。这里我直接用最基础的Python科学计算栈需要numpy、scikit-learn、pandas和matplotlib。如果你的环境里还没有安装可以执行这条命令pip install numpy pandas scikit-learn matplotlib为了便于演示我构造一组模拟的近红外光谱数据。在实际项目中你会把数据换成自己仪器导出的光谱矩阵但数据结构都一样X是样本数乘以波长数的二维矩阵y是样本数乘以1的浓度向量。这里我用200个样本、100个波长变量的模拟数据来跑通流程。import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split from sklearn.cross_decomposition import PLSRegression from sklearn.metrics import mean_squared_error, r2_score # 模拟光谱数据200行100个波长变量 np.random.seed(42) n_samples 200 n_wavelengths 100 wavelengths np.linspace(4000, 10000, n_wavelengths) # 真实有效波长索引模拟4个物质成分对应的吸收区域 true_wavelengths [20, 35, 55, 78, 90] X np.zeros((n_samples, n_wavelengths)) y np.zeros(n_samples) for i in range(n_samples): for j in true_wavelengths: # 生成峰宽较大的高斯吸收峰 X[i] (np.random.rand() * 2 - 1) * np.exp(-((wavelengths - wavelengths[j]) / 200) ** 2) noise np.random.normal(0, 0.02, n_wavelengths) X[i] noise y[i] 10 X[i, true_wavelengths].sum() * 2.5 np.random.normal(0, 0.5) X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42)这段代码生成了带高斯峰形状的模拟光谱其中只有5个波长位置对浓度有真实贡献其余变量基本是噪声和冗余信息。用这样的数据来检验SPA我们能直观看到算法能不能把这5个关键变量从噪声里捞出来。3.2 核心代码从零实现SPA算法下面这部分是SPA的核心实现。我尽量写得清楚直接便于理解和二次开发。整个类接收输入矩阵X、候选变量个数上限N以及需要评估的交叉验证折数。class SPA: def __init__(self, max_num_vars20): self.max_num_vars max_num_vars def _projection(self, X): 计算所有变量在指定方向上的投影距离 pass # 在fit里执行 def fit(self, X, y): X: 样本×变量 y: 样本×1 返回最佳波长索引集合 n_samples, n_vars X.shape # 1. 数据标准化避免量纲影响 Xnorm (X - X.mean(axis0)) / X.std(axis0) best_rmse np.inf best_indices None # 2. 分别从每个变量出发构建候选变量子集 for start_idx in range(n_vars): selected [start_idx] # 逐步选出 N-1 个投影距离最大的变量 for _ in range(1, self.max_num_vars): projection_distances np.full(n_vars, np.inf) last_selected selected[-1] z Xnorm[:, last_selected] # 计算所有未选变量的投影距离 for v in range(n_vars): if v in selected: continue x Xnorm[:, v] inner_product np.dot(x, z) z_norm_sq np.dot(z, z) if z_norm_sq 1e-12: continue # 投影向量 proj (inner_product / z_norm_sq) * z proj_dist np.linalg.norm(x - proj) projection_distances[v] proj_dist # 选投影距离最大的 next_var np.argmax(projection_distances) selected.append(next_var) # 3. 用MLR评估当前候选子集 X_selected X[:, selected] X_selected np.column_stack([np.ones(n_samples), X_selected]) # 简单的MLR系数估计beta (XX)^-1 Xy beta np.linalg.pinv(X_selected.T X_selected) X_selected.T y y_pred X_selected beta rmse np.sqrt(mean_squared_error(y, y_pred)) if rmse best_rmse: best_rmse rmse best_indices selected return np.array(best_indices), best_rmse这个实现做了两件事对于从每一个变量出发的搜索路径都完整地跑完投影迭代与MLR评估最后返回RMSE最低的那组波长索引。需要说明的是在真实项目中我们一般用验证集或交叉验证来评估而不是直接用训练集RMSE否则很容易过拟合。后面我会补一个交叉验证的版本。3.3 用交叉验证选择最优N与波长组合直接用训练集RMSE评估候选子集会把那些在训练集上表现好但泛化能力差的组合选出来。更好的做法是在每个候选子集上做K折交叉验证用平均交叉验证均方根误差RMSECV作为评价标准。下面我把评估部分替换成5折交叉验证同时循环尝试不同的变量个数N从1到20逐个跑输出RMSECV随N变化的曲线。from sklearn.model_selection import KFold def spa_cv_selection(X, y, max_n20, cv_folds5): Xnorm (X - X.mean(axis0)) / X.std(axis0) kf KFold(n_splitscv_folds, shuffleTrue, random_state42) results [] best_rmsecv np.inf best_combo None for start_idx in range(X.shape[1]): selected [start_idx] for step in range(1, max_n): last_z Xnorm[:, selected[-1]] proj_dists np.full(X.shape[1], -np.inf) for v in range(X.shape[1]): if v in selected: continue x Xnorm[:, v] z last_z inner np.dot(x, z) znorm_sq np.dot(z, z) if znorm_sq 1e-12: continue proj (inner / znorm_sq) * z proj_dists[v] np.linalg.norm(x - proj) next_v np.argmax(proj_dists) selected.append(next_v) # 对当前N值做交叉验证 rmsecv calculate_rmsecv(X[:, selected], y, kf) results.append((len(selected), rmsecv, selected.copy())) if rmsecv best_rmsecv: best_rmsecv rmsecv best_combo selected.copy() return best_combo, best_rmsecv, results def calculate_rmsecv(X_sub, y, kf): rmses [] for train_idx, val_idx in kf.split(X_sub): X_train, X_val X_sub[train_idx], X_sub[val_idx] y_train, y_val y[train_idx], y[val_idx] model np.linalg.pinv(X_train.T X_train) X_train.T y_train y_pred X_val model rmses.append(np.sqrt(mean_squared_error(y_val, y_pred))) return np.mean(rmses) best_combo, best_rmsecv, results spa_cv_selection(X_train, y_train, max_n20) print(最优波长索引:, best_combo) print(最优RMSECV:, best_rmsecv)在实际运行结果里SPA选出的波长索引通常能落在我们在模拟数据中预设的那5个位置附近哪怕初始噪声很大它也能通过多轮搜索找到接近最优的组合。这个稳定性正是它在真实光谱项目里可靠的原因。3.4 模型效果对比全谱PLS vs SPA-MLR选出了波长变量还要验证效果确实有提升。我习惯把全谱PLS模型和SPA-MLR模型放在同一批测试集上做对比。PLS是近红外建模的常青树和SPA-MLR对比是最有说服力的。# 全谱PLS模型 pls PLSRegression(n_components8) pls.fit(X_train, y_train) y_pred_pls pls.predict(X_test) rmse_pls np.sqrt(mean_squared_error(y_test, y_pred_pls)) r2_pls r2_score(y_test, y_pred_pls) print(fPLS全谱: RMSE{rmse_pls:.4f}, R2{r2_pls:.4f}) # SPA-MLR模型只用选出的波长 X_train_s X_train[:, best_combo] X_test_s X_test[:, best_combo] mlr_coef np.linalg.pinv(X_train_s.T X_train_s) X_train_s.T y_train y_pred_spa X_test_s mlr_coef rmse_spa np.sqrt(mean_squared_error(y_test, y_pred_spa)) r2_spa r2_score(y_test, y_pred_spa) print(fSPA-MLR: RMSE{rmse_spa:.4f}, R2{r2_spa:.4f})运行这个对比你会发现SPA-MLR在测试集上的RMSE通常与全谱PLS相当甚至略好而它使用的变量数量可能只有后者的几十分之一。这意味着更短的建模时间、更易解释的模型以及在实际仪器部署时更少的计算负担。对工业生产线上实时检测这种对速度有要求的场景这个优势是实打实的。3.5 从候选变量中观察吸收峰位置可解释性检查跑完算法之后千万别急着收工。我会习惯性地把选出的波长位置标在原始光谱平均谱图上看看它们是不是落在某些特征吸收峰附近。这一步做的是“物理意义校验”。mean_spectrum X_train.mean(axis0) plt.figure(figsize(12, 5)) plt.plot(wavelengths, mean_spectrum, lw1.2, labelMean spectrum) for idx in best_combo: plt.axvline(xwavelengths[idx], colorred, linestyle--, alpha0.6) plt.xlabel(Wavenumber (cm⁻¹)) plt.ylabel(Absorbance) plt.legend() plt.show()如果选出的波长全部落在平坦无特征的区域而且模型表现又好这通常暗示数据可能有潜在的问题比如谱图没有对齐、样本泄漏或者预处理引入伪信息。反之如果选出的波长恰好对应某些特征基团的吸收区域那模型的物理可解释性就大大增强验证集外推时也更让人放心。4. 常见问题与排查技巧实录4.1 最终变量个数N到底怎么定我在真实项目里反复试过很多定N的策略比较可靠的做法不是机械套用一个固定值而是观察RMSECV随N变化的曲线。当变量个数从1逐步增加到20时RMSECV通常先快速下降然后进入平台期偶尔会有一个小幅度的最低点。我一般选择平台期开始前的那个点而不是绝对最低点因为绝对最低点附近往往已经开始过拟合了。举个例子有次做汽油辛烷值的近红外模型RMSECV在N7时达到最低点但从N5开始下降幅度就已经少于1%了而N7的模型在后续新批次的测试上反而不如N5稳定。从那以后我就坚持一个原则选曲线上“收益递减”的拐点而不是纯数学的最小值。这条经验在多个数据集上都适用。4.2 初始波长怎么选为什么要遍历所有变量SPA的搜索结果依赖初始变量的选择。如果初始波长选在噪声区域且后续搜索路径受局部结构影响最后拿到的候选子集可能不是全局最优。所以实际做SPA的时候我不会指定某一个初始波长而会老老实实把所有变量都作为初始点跑一遍再把每轮得到的子集放在同一评估标准下比较。这种“暴力枚举初始点”的方式虽然会让计算时间乘以变量数的倍数但对几百到几千个变量的光谱数据来说现代计算机跑下来也就几秒钟到几十秒的事完全值回票价。有一个小技巧是如果变量特别多比如高光谱图像上万个波段可以先做个粗筛把明显属于噪声区域的变量剔除再对剩下几百个候选波段跑SPA速度提升非常明显。4.3 数据预处理对SPA的影响标准化和均值中心化SPA是基于向量投影的算法对变量的绝对尺度十分敏感。如果不同波长位置的响应值量级差很大投影距离就会被量级大的变量主导导致“选大不选对”。所以在跑SPA之前我通常会对光谱矩阵做标准化或均值中心化处理。标准化是把每个波长变量的均值变为0、方差变为1这样每个波长在投影计算中是“公平”的均值中心化则只去除基线漂移适合本身量纲一致的光谱数据。在近红外光谱中吸光度本身量纲一致用均值中心化就够。但如果你处理的是其他类型的数据或者光谱数值范围差异大就优先用标准化。我的建议是两种都试一遍看最终模型的交叉验证表现再决定这个对比成本很低。4.4 SPA与CARS、UVE结论不一致怎么办我经常被问到一个问题同一份数据SPA、CARS、UVE选出来的变量不一样而且差异很大到底信谁的很多新手这时候会陷入选择困难。我的回答是变量筛选结果不一致是常态不是异常。不同算法的目标函数不一样SPA追求低共线性CARS追求强相关与稳定性UVE剔除噪声变量它们从不同的角度压缩信息最后留下的变量自然不可能完全相同。正确的处理方式不是争论谁对谁错而是把不同方法的结果放到同一套验证框架里做公平对比让测试集RMSE说话。我之前做过一个项目SPA和CARS选出的变量几乎没有交集但两个模型的测试性能都很接近。最后我选了SPA的结果因为它的变量在谱图上对应已知的吸收特征解释成本更低。4.5 实际化工项目中的踩坑样本划分和模型外推最后聊一个我在真实项目里踩过的大坑。早年间做某个化工在线监测项目实验室里用留出法验证SPA-MLR模型效果特别好R²高达0.98但一到现场就崩。排查了很久才发现问题出在样本划分上。光谱数据天然具有批次效应同一个批次内的样本光谱特征高度相似如果随机划分训练集和测试集同一批次的样本会同时出现在两边模型等于考试前偷看到了答案验证结果自然虚高。正确做法是按批次划分样本集保证训练集和测试集来自不同的批次或者至少使用按时间顺序划分的验证集。这个原则不仅适用于SPA对所有的光谱建模项目都成立。我后来做任何变量筛选评估都会先确认数据划分方式是否真的模拟了未来的使用场景否则算法调得再好也是白搭。在实际操作中我还发现一个很实用的小技巧SPA选出的波长变量可以在不同批次的样本上重新做一次简单的MLR系数稳定性检查。如果系数方向和大小在不同批次间保持稳定说明这些变量携带的是真实化学信息如果系数波动剧烈那大概率选出的只是统计巧合的噪声变量。这个检查几乎不花时间却能在模型部署前帮我们排掉一大半的隐患。本文还有配套的精品资源点击获取 连续投影算法(SPA)详解：光谱变量筛选原理与Python实战
简介本资源是一份面向数据科学初学者与光谱分析从业者的SPA连续投影算法原理精讲与工程实践教程聚焦高维数据降维这一核心问题特别适用于生物信息学、化学计量学及机器学习预处理场景。压缩包共14个文件6.55MB含11个MATLAB源码.m——涵盖spa_train、spa、validation等主流程与子函数2个实测光谱数据集.mat包括jasperRidge2_R198.mat等典型高维光谱样本1份结构清晰的中文教程文档.docx系统梳理算法原理、迭代逻辑、调试要点与光谱应用案例。已有2413人学习下载读者可直接运行调试后的完整代码链复现从数据加载、投影迭代、方差评估到分类预测的全流程并获得作者在内存管理与收敛性优化方面的实操经验显著降低算法落地门槛。 做光谱建模的人十有八九都会碰到同一个问题波长变量太多样本数又不够模型老是过拟合今天调出来的效果明天就翻车。我早年做近红外定量分析的时候一套光谱数据经常是上千个波长点但实际有效信息可能就集中在几十个波段里。那时候最苦恼的就是怎么从一堆冗余变量里挑出真正有用的那部分直到后来接触了连续投影算法Successive Projections Algorithm简称SPA才算是把这条弯路走顺了。SPA的核心工作很简单它能在高维光谱数据中自动挑出一组共线性最小、信息冗余最低的特征波长再用这些筛选出来的变量去建校正模型。它不依赖任何特定模型框架和偏最小二乘PLS、多元线性回归MLR都能搭配使用很适合做光谱定量分析的预先降维步骤。这篇教程我会从算法原理一步步拆开讲配上完整的Python实现和实际案例最后聊聊我在工程里踩过的坑。无论是刚入门化学计量学的学生还是已经在做光谱建模的工程师这篇文章都应该能帮到你。1. 为什么需要变量选择光谱数据里的共线性问题1.1 光谱建模的痛点高维、共线性、过拟合先说说我们面对的原始数据长什么样。以近红外光谱为例一台光谱仪扫描一个样本可能在4000到10000波数范围内采集上千个数据点每个数据点就是一个变量。如果我有200个样本那数据矩阵就是200乘1000的维度。用这样的矩阵去建常规的多元线性回归模型数学上是不可行的因为变量数远远大于样本数即使强行用偏最小二乘去压缩潜变量也会面临信息分散、模型解释困难的问题。更麻烦的是光谱数据里相邻波长的响应不是独立的而是高度相关的。某个官能团的吸收峰往往横跨几十个波数对应几十个变量之间都有极强的关联。这种共线性会让很多回归算法不稳定你今天采样环境稍微变一点模型系数就可能发生剧烈波动。说白了变量之间大量重复的信息不仅没帮上忙反而成了噪声放大器。这个问题的本质是我们需要在保证预测能力的前提下找到一组信息互补、冗余最小的波长变量。这好比你请一个团队干活如果所有人都是同一种能力结构团队整体战斗力反而很弱真正理想的团队是每个人都能填补其他人的盲区。SPA想要做的就是这件事找出那些彼此最“不重复”的波长变量。1.2 变量选择算法选型对比SPA、CARS、UVE、VIP光谱领域里变量筛选算法不少我自己先后试过VIP变量重要性投影、UVE无信息变量消除、CARS竞争性自适应重加权采样还有SPA。每种方法的思路差别挺大适用场景也不一样。先看VIP它是在PLS模型的框架下根据每个变量对模型解释方差的贡献度来打分的好处是计算快、直观缺点是它本质上只衡量了“重要性”并没有考虑变量之间的共线性结构。UVE的做法是把变量和一组随机噪声变量放在一起比较稳定性剔除那些稳定性低于噪声的变量这个方法能有效去噪但有时候也会误删一些有用但贡献不强的小信号。CARS是模拟达尔文进化论的思想通过自适应重加权采样不断淘汰“劣质”变量收敛快但因为有随机抽样环节每次运行结果会有一定波动需要跑多次取共识。SPA和上面几个都不太一样。它不是从统计显著性出发而是从几何角度衡量变量之间的关系。它逐步选出与当前已选变量最少冗余的波长每一步都保证新选的变量和前面选出的变量在向量空间上尽量正交从根源上砍掉共线性。这种思路在光谱数据这种共线性很强的场景里效果特别明显。1.3 连续投影算法的核心思路一句大白话SPA到底在做什么我用一句最直白的话概括它每次从剩下的变量里挑一个让这个变量的光谱响应向量与所有已选变量向量之间的“重叠程度”最低也就是在向量空间中对应的投影距离最大然后把信息最不重复的那个波长收进来。这样一轮一轮挑下去最终得到一组彼此“各自独立”的波长组合。这里的“投影”就是一个向量在另一个向量上的正交投影长度。如果两个变量谱图长得几乎一模一样那么一个变量在另一个变量上的投影就会很长两者冗余程度就很高。反过来如果两个变量谱图差异很大投影长度就越小。SPA的策略就是每一轮都选投影距离最大的那个变量让新选中的变量和已有变量整体冗余度最低。2. SPA的数学原理与步骤拆解2.1 投影计算与共线性度量从向量投影说起要动手实现SPA第一步得懂向量投影运算。假设我们有两个波长变量对应的响应向量分别记作向量x和z每个向量都是样本维度的列向量。我们把x投影到z上的投影向量记作proj_z(x)计算公式是proj_z(x) (x · z / (z · z)) · z这里x · z是内积z · z是z自身的内积。投影向量是z方向上的一个向量它表示x中能被z“解释”掉的部分。如果x和z完全一样完全线性相关那么投影向量就等于x自身也就是说x的信息完全冗余于z如果x和z垂直则投影向量为零向量说明两者完全没有冗余。SPA在筛选变量时对于当前选中向量z它计算每个候选变量x在z上的投影向量长度即P(x) || x - proj_z(x) ||这个值表示x中与z不重叠的信息量。P(x)越大说明x携带的独立信息越多。SPA每一轮就选P(x)最大的那个变量作为下一个入选波长。这也解释了为什么它天然适合抵抗共线性两个相关性极高的变量之间P(x)会非常小几乎没机会同时被选中。2.2 完整算法流程初始化、迭代投影、候选变量集合完整SPA算法的设计比上面描述的要周详一些它并不是只指定一个起始变量然后一路选到底而是从不同的初始变量出发生成多组候选变量子集再挑选整体效果最好的那一组。我将SPA的完整流程拆解成四个阶段阶段一确定初始变量。在全部变量中分别以每一个变量作为初始入选变量然后分别执行后面的投影迭代。也就是说如果原始数据有1000个波长那么会有1000轮独立的搜索每一轮都从不同的初始波长开始。阶段二循环投影选择。在每一轮独立的搜索里设置一个预设的最大选变量数N一般取10到30之间或者根据经验公式定然后从初始变量出发迭代执行N-1次投影计算。每次迭代时计算所有尚未入选的变量与最近一次入选变量之间的投影距离选出投影距离最大的变量作为新的入选变量。阶段三形成候选子集。每一轮搜索结束都会得到一组长度等于N的波长索引序列。把所有初始变量跑完之后我们手上就有了与变量总数相同数量的候选子集每个子集包含N个波长变量。阶段四用校正模型评估。对每一个候选子集使用选出的N个波长变量建立多元线性回归MLR模型计算验证集的均方根误差RMSE选验证误差最小的那组波长作为最终输出。这种“多初始点搜索模型评估”的双层结构保证了SPA不会因为初始变量选得太差而错过全局最优子集。虽然计算量有所增加但对几百个变量的光谱数据来说完全在可接受范围内。2.3 为什么选“投影最大”而不是“相关最小”一个容易混淆的地方很多初学者会在这里卡住。他们觉得SPA要选“共线性最小”的变量那直接找相关系数最小的两个波长不就行了为什么非要绕一圈算投影、选投影距离最大我当年也绕了很久才想明白这中间的差别。两变量之间的相关系数低只说明它们之间没有明显的线性趋势但并不能保证它们在多元回归中不互相干扰。而SPA要的是变量对整体信息空间的“增量贡献”。它关心的是在已经选了某些变量的前提下再增加哪一个变量能带来最多的新信息。投影距离衡量的恰恰是“新增信息量”。比如变量A和变量B的相关系数是0.01看起来很独立但是B其实只是A经过一个微小变换得到的那么B在A方向上的投影可能仍然很大意味着B并没有带来多少新东西。相关系数小并不代表冗余度低这种微妙的差距正是SPA选择投影距离作为准则的原因。2.4 与MLR结合的SPA-MLR变量数怎么定SPA输出波长集合后还需要搭配一个校正模型来形成完整预测方案最常见的搭配是MLR也就是标准的最小二乘多元线性回归。因为经过SPA筛选后的变量数量已经从上千降到十几个MLR在数学上变得可行而且模型解释性比PLS强很多直接能看到每个波长对目标的贡献系数。变量数N怎么定这是一个很关键的问题。定得太少信息不足模型拟合差定得太高筛选变量过多MLR可能重新陷入过拟合。一个实用的经验值是N不超过样本数的1/10比如你有200个样本N最大选到20。实际操作中一般设置一个N的上限然后在1到这个上限之间逐一尝试观察交叉验证RMSECV的变化。通常RMSECV会先随着变量数增加而下降到某个点之后开始回升或趋于平稳选择回升前的那个拐点对应的N即可。3. Python完整实现与实操案例3.1 环境准备与数据说明动手之前我们先准备好环境。这里我直接用最基础的Python科学计算栈需要numpy、scikit-learn、pandas和matplotlib。如果你的环境里还没有安装可以执行这条命令pip install numpy pandas scikit-learn matplotlib为了便于演示我构造一组模拟的近红外光谱数据。在实际项目中你会把数据换成自己仪器导出的光谱矩阵但数据结构都一样X是样本数乘以波长数的二维矩阵y是样本数乘以1的浓度向量。这里我用200个样本、100个波长变量的模拟数据来跑通流程。import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split from sklearn.cross_decomposition import PLSRegression from sklearn.metrics import mean_squared_error, r2_score # 模拟光谱数据200行100个波长变量 np.random.seed(42) n_samples 200 n_wavelengths 100 wavelengths np.linspace(4000, 10000, n_wavelengths) # 真实有效波长索引模拟4个物质成分对应的吸收区域 true_wavelengths [20, 35, 55, 78, 90] X np.zeros((n_samples, n_wavelengths)) y np.zeros(n_samples) for i in range(n_samples): for j in true_wavelengths: # 生成峰宽较大的高斯吸收峰 X[i] (np.random.rand() * 2 - 1) * np.exp(-((wavelengths - wavelengths[j]) / 200) ** 2) noise np.random.normal(0, 0.02, n_wavelengths) X[i] noise y[i] 10 X[i, true_wavelengths].sum() * 2.5 np.random.normal(0, 0.5) X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42)这段代码生成了带高斯峰形状的模拟光谱其中只有5个波长位置对浓度有真实贡献其余变量基本是噪声和冗余信息。用这样的数据来检验SPA我们能直观看到算法能不能把这5个关键变量从噪声里捞出来。3.2 核心代码从零实现SPA算法下面这部分是SPA的核心实现。我尽量写得清楚直接便于理解和二次开发。整个类接收输入矩阵X、候选变量个数上限N以及需要评估的交叉验证折数。class SPA: def __init__(self, max_num_vars20): self.max_num_vars max_num_vars def _projection(self, X): 计算所有变量在指定方向上的投影距离 pass # 在fit里执行 def fit(self, X, y): X: 样本×变量 y: 样本×1 返回最佳波长索引集合 n_samples, n_vars X.shape # 1. 数据标准化避免量纲影响 Xnorm (X - X.mean(axis0)) / X.std(axis0) best_rmse np.inf best_indices None # 2. 分别从每个变量出发构建候选变量子集 for start_idx in range(n_vars): selected [start_idx] # 逐步选出 N-1 个投影距离最大的变量 for _ in range(1, self.max_num_vars): projection_distances np.full(n_vars, np.inf) last_selected selected[-1] z Xnorm[:, last_selected] # 计算所有未选变量的投影距离 for v in range(n_vars): if v in selected: continue x Xnorm[:, v] inner_product np.dot(x, z) z_norm_sq np.dot(z, z) if z_norm_sq 1e-12: continue # 投影向量 proj (inner_product / z_norm_sq) * z proj_dist np.linalg.norm(x - proj) projection_distances[v] proj_dist # 选投影距离最大的 next_var np.argmax(projection_distances) selected.append(next_var) # 3. 用MLR评估当前候选子集 X_selected X[:, selected] X_selected np.column_stack([np.ones(n_samples), X_selected]) # 简单的MLR系数估计beta (XX)^-1 Xy beta np.linalg.pinv(X_selected.T X_selected) X_selected.T y y_pred X_selected beta rmse np.sqrt(mean_squared_error(y, y_pred)) if rmse best_rmse: best_rmse rmse best_indices selected return np.array(best_indices), best_rmse这个实现做了两件事对于从每一个变量出发的搜索路径都完整地跑完投影迭代与MLR评估最后返回RMSE最低的那组波长索引。需要说明的是在真实项目中我们一般用验证集或交叉验证来评估而不是直接用训练集RMSE否则很容易过拟合。后面我会补一个交叉验证的版本。3.3 用交叉验证选择最优N与波长组合直接用训练集RMSE评估候选子集会把那些在训练集上表现好但泛化能力差的组合选出来。更好的做法是在每个候选子集上做K折交叉验证用平均交叉验证均方根误差RMSECV作为评价标准。下面我把评估部分替换成5折交叉验证同时循环尝试不同的变量个数N从1到20逐个跑输出RMSECV随N变化的曲线。from sklearn.model_selection import KFold def spa_cv_selection(X, y, max_n20, cv_folds5): Xnorm (X - X.mean(axis0)) / X.std(axis0) kf KFold(n_splitscv_folds, shuffleTrue, random_state42) results [] best_rmsecv np.inf best_combo None for start_idx in range(X.shape[1]): selected [start_idx] for step in range(1, max_n): last_z Xnorm[:, selected[-1]] proj_dists np.full(X.shape[1], -np.inf) for v in range(X.shape[1]): if v in selected: continue x Xnorm[:, v] z last_z inner np.dot(x, z) znorm_sq np.dot(z, z) if znorm_sq 1e-12: continue proj (inner / znorm_sq) * z proj_dists[v] np.linalg.norm(x - proj) next_v np.argmax(proj_dists) selected.append(next_v) # 对当前N值做交叉验证 rmsecv calculate_rmsecv(X[:, selected], y, kf) results.append((len(selected), rmsecv, selected.copy())) if rmsecv best_rmsecv: best_rmsecv rmsecv best_combo selected.copy() return best_combo, best_rmsecv, results def calculate_rmsecv(X_sub, y, kf): rmses [] for train_idx, val_idx in kf.split(X_sub): X_train, X_val X_sub[train_idx], X_sub[val_idx] y_train, y_val y[train_idx], y[val_idx] model np.linalg.pinv(X_train.T X_train) X_train.T y_train y_pred X_val model rmses.append(np.sqrt(mean_squared_error(y_val, y_pred))) return np.mean(rmses) best_combo, best_rmsecv, results spa_cv_selection(X_train, y_train, max_n20) print(最优波长索引:, best_combo) print(最优RMSECV:, best_rmsecv)在实际运行结果里SPA选出的波长索引通常能落在我们在模拟数据中预设的那5个位置附近哪怕初始噪声很大它也能通过多轮搜索找到接近最优的组合。这个稳定性正是它在真实光谱项目里可靠的原因。3.4 模型效果对比全谱PLS vs SPA-MLR选出了波长变量还要验证效果确实有提升。我习惯把全谱PLS模型和SPA-MLR模型放在同一批测试集上做对比。PLS是近红外建模的常青树和SPA-MLR对比是最有说服力的。# 全谱PLS模型 pls PLSRegression(n_components8) pls.fit(X_train, y_train) y_pred_pls pls.predict(X_test) rmse_pls np.sqrt(mean_squared_error(y_test, y_pred_pls)) r2_pls r2_score(y_test, y_pred_pls) print(fPLS全谱: RMSE{rmse_pls:.4f}, R2{r2_pls:.4f}) # SPA-MLR模型只用选出的波长 X_train_s X_train[:, best_combo] X_test_s X_test[:, best_combo] mlr_coef np.linalg.pinv(X_train_s.T X_train_s) X_train_s.T y_train y_pred_spa X_test_s mlr_coef rmse_spa np.sqrt(mean_squared_error(y_test, y_pred_spa)) r2_spa r2_score(y_test, y_pred_spa) print(fSPA-MLR: RMSE{rmse_spa:.4f}, R2{r2_spa:.4f})运行这个对比你会发现SPA-MLR在测试集上的RMSE通常与全谱PLS相当甚至略好而它使用的变量数量可能只有后者的几十分之一。这意味着更短的建模时间、更易解释的模型以及在实际仪器部署时更少的计算负担。对工业生产线上实时检测这种对速度有要求的场景这个优势是实打实的。3.5 从候选变量中观察吸收峰位置可解释性检查跑完算法之后千万别急着收工。我会习惯性地把选出的波长位置标在原始光谱平均谱图上看看它们是不是落在某些特征吸收峰附近。这一步做的是“物理意义校验”。mean_spectrum X_train.mean(axis0) plt.figure(figsize(12, 5)) plt.plot(wavelengths, mean_spectrum, lw1.2, labelMean spectrum) for idx in best_combo: plt.axvline(xwavelengths[idx], colorred, linestyle--, alpha0.6) plt.xlabel(Wavenumber (cm⁻¹)) plt.ylabel(Absorbance) plt.legend() plt.show()如果选出的波长全部落在平坦无特征的区域而且模型表现又好这通常暗示数据可能有潜在的问题比如谱图没有对齐、样本泄漏或者预处理引入伪信息。反之如果选出的波长恰好对应某些特征基团的吸收区域那模型的物理可解释性就大大增强验证集外推时也更让人放心。4. 常见问题与排查技巧实录4.1 最终变量个数N到底怎么定我在真实项目里反复试过很多定N的策略比较可靠的做法不是机械套用一个固定值而是观察RMSECV随N变化的曲线。当变量个数从1逐步增加到20时RMSECV通常先快速下降然后进入平台期偶尔会有一个小幅度的最低点。我一般选择平台期开始前的那个点而不是绝对最低点因为绝对最低点附近往往已经开始过拟合了。举个例子有次做汽油辛烷值的近红外模型RMSECV在N7时达到最低点但从N5开始下降幅度就已经少于1%了而N7的模型在后续新批次的测试上反而不如N5稳定。从那以后我就坚持一个原则选曲线上“收益递减”的拐点而不是纯数学的最小值。这条经验在多个数据集上都适用。4.2 初始波长怎么选为什么要遍历所有变量SPA的搜索结果依赖初始变量的选择。如果初始波长选在噪声区域且后续搜索路径受局部结构影响最后拿到的候选子集可能不是全局最优。所以实际做SPA的时候我不会指定某一个初始波长而会老老实实把所有变量都作为初始点跑一遍再把每轮得到的子集放在同一评估标准下比较。这种“暴力枚举初始点”的方式虽然会让计算时间乘以变量数的倍数但对几百到几千个变量的光谱数据来说现代计算机跑下来也就几秒钟到几十秒的事完全值回票价。有一个小技巧是如果变量特别多比如高光谱图像上万个波段可以先做个粗筛把明显属于噪声区域的变量剔除再对剩下几百个候选波段跑SPA速度提升非常明显。4.3 数据预处理对SPA的影响标准化和均值中心化SPA是基于向量投影的算法对变量的绝对尺度十分敏感。如果不同波长位置的响应值量级差很大投影距离就会被量级大的变量主导导致“选大不选对”。所以在跑SPA之前我通常会对光谱矩阵做标准化或均值中心化处理。标准化是把每个波长变量的均值变为0、方差变为1这样每个波长在投影计算中是“公平”的均值中心化则只去除基线漂移适合本身量纲一致的光谱数据。在近红外光谱中吸光度本身量纲一致用均值中心化就够。但如果你处理的是其他类型的数据或者光谱数值范围差异大就优先用标准化。我的建议是两种都试一遍看最终模型的交叉验证表现再决定这个对比成本很低。4.4 SPA与CARS、UVE结论不一致怎么办我经常被问到一个问题同一份数据SPA、CARS、UVE选出来的变量不一样而且差异很大到底信谁的很多新手这时候会陷入选择困难。我的回答是变量筛选结果不一致是常态不是异常。不同算法的目标函数不一样SPA追求低共线性CARS追求强相关与稳定性UVE剔除噪声变量它们从不同的角度压缩信息最后留下的变量自然不可能完全相同。正确的处理方式不是争论谁对谁错而是把不同方法的结果放到同一套验证框架里做公平对比让测试集RMSE说话。我之前做过一个项目SPA和CARS选出的变量几乎没有交集但两个模型的测试性能都很接近。最后我选了SPA的结果因为它的变量在谱图上对应已知的吸收特征解释成本更低。4.5 实际化工项目中的踩坑样本划分和模型外推最后聊一个我在真实项目里踩过的大坑。早年间做某个化工在线监测项目实验室里用留出法验证SPA-MLR模型效果特别好R²高达0.98但一到现场就崩。排查了很久才发现问题出在样本划分上。光谱数据天然具有批次效应同一个批次内的样本光谱特征高度相似如果随机划分训练集和测试集同一批次的样本会同时出现在两边模型等于考试前偷看到了答案验证结果自然虚高。正确做法是按批次划分样本集保证训练集和测试集来自不同的批次或者至少使用按时间顺序划分的验证集。这个原则不仅适用于SPA对所有的光谱建模项目都成立。我后来做任何变量筛选评估都会先确认数据划分方式是否真的模拟了未来的使用场景否则算法调得再好也是白搭。在实际操作中我还发现一个很实用的小技巧SPA选出的波长变量可以在不同批次的样本上重新做一次简单的MLR系数稳定性检查。如果系数方向和大小在不同批次间保持稳定说明这些变量携带的是真实化学信息如果系数波动剧烈那大概率选出的只是统计巧合的噪声变量。这个检查几乎不花时间却能在模型部署前帮我们排掉一大半的隐患。本文还有配套的精品资源点击获取