简介:压缩包内是一套围绕iVISSA光谱特征波段筛选方法的MATLAB实现与配套示例数据,面向遥感、环境科学、农业检测等方向的研究人员与工程师,适用于从高维光谱中挑选有效波段、降低数据冗余并提升定量预测模型精度的场景。包内共12个文件,其中9个.m脚本构成核心算法链,覆盖光谱预处理、PLS建模、iVISSA特征筛选、交叉验证与结果预测等环节;2个.mat文件提供大豆水分光谱原始数据及筛选结果;1份txt文档说明软件许可范围。整个压缩包仅157KB,部署和使用都非常轻量。目前已有1682人学习下载,适合具备一定光谱分析或建模基础、希望快速搭建特征筛选流程的读者。两份可直接运行的示例脚本串联了从数据预处理、特征提取到波段选择、模型评估的完整链路,既能用于理解算法内部原理,也可作为科研实验对照、复现已有结果或在此基础上做二次开发的基础。
1. 光谱特征波段筛选的可选解法中,iVISSA 的核心优势藏在哪
一条近红外光谱通常包含几百甚至上千个波长点,但真正与目标成分挂钩的波段往往只有十几个。直接拿全谱去建 PLS 模型,最典型的后果是训练集表现尚可,一换仪器或者换批次就失控,因为大量不相关变量会把主成分方向带偏。光谱特征波段筛选就是为了解决这类“高维小样本”问题:在尽量不损失预测精度的前提下,把参与建模的变量压到可解释、可转移、可维护的规模。
iVISSA(Iterative Variable Space Shrinkage Analysis,迭代变量空间收缩分析)正是面向这个需求的一种变量选择算法。它不靠穷举组合,而是通过反复随机抽样、模型误差反馈和变量空间收缩,把候选波长逐步推向可信区域。相比 CARS、SPA 这类常见方法,它在强相关的近红外连续谱区里往往能得到更连贯、更容易从分子振动归属上解释的波段,这也是近几年很多定标流程里把它作为默认筛选方案的原因。
这篇文章会把 VISSA 到 iVISSA 的演进思路、一段能直接跑的 Python 实现、四个必调参数的真实边界,以及实践中容易翻车的五个细节讲清楚。适合刚开始接触光谱变量选择的研究生,也适合从 Matlab 转向 Python 的检测一线工程师。
2. 从 VISSA 到 iVISSA:变量空间收缩的数学直觉和前向选择做不到的事
2.1 全谱建模的瓶颈与变量选择的价值换算
光谱矩阵的形状通常是样本数 N×变量数 P,当 P 到几千、N 只有几十时,数据天然病态。偏最小二乘回归虽然通过潜变量做了降维,但潜变量是所有原始变量的加权组合,不相关变量会稀释有效信息。换句话说,每增加一个无关波长,都在给主成分方向拖后腿——不是“不加分”,而是直接“扣分”。
变量选择的价值怎么衡量?工程上我习惯看三个数:筛后变量占比、RMSECV 和 RMSEP 的变化、模型转移成本。从两千个波长筛到三十个,RMSECV 不升、RMSEP 降,这是最理想的剧本。更现实的价值是模型转移:特征波段被压缩到特定化学键振动区间后,两台仪器在同一批窄波段上的响应一致性通常远好于全谱,模型直接共用或只做一次斜率截距校正就够,这比任何精细预处理都好使。
这个方向上有几个成熟算法,CARS 用自适应重加权采样淘汰低权重变量,SPA 用投影方法找近似线性独立的波长,遗传算法在变量组合上做全局搜索。但它们各自有软肋:结果对随机种子敏感,停机条件要靠经验,更重要的是选出的变量往往零散分布,化学意义解释困难。iVISSA 的思路不同,它是从完整变量空间出发,靠迭代收缩逼近一个可信的变量子集,保留的波段通常是完整谱区连段,归因远比其他方法省事。
2.2 VISSA 的两步迭代:蒙特卡洛抽样、加权投票、空间收缩
VISSA 的核心逻辑是空间收缩。初始变量空间就是全谱,从它开始,每轮在当前空间内多次随机抽样,对每个子集建立 PLS 模型,用交叉验证误差评估变量贡献,然后根据贡献权重裁掉低分变量,进入下一轮。收缩到最后留下的变量集合,就是特征波段。
关键动作有两个,缺一不可。第一是蒙特卡洛抽样。每轮在当前空间里按固定比例随机抽取变量,重复几十次。只有抽样次数够多,变量间的共线性和冗余才能在统计意义上体现出来。一个完全无关的变量被抽进好子集是偶尔事件,而一个贡献稳定的变量会在大量随机组合中反复出现、持续拉低误差,这种差异会被逐步放大。
第二是加权投票机制。我常用的实现是:记录每个变量在所有包含它的子集里的平均交叉验证误差,和它在未出现时的平均误差做差。差值越明显,说明这个变量的加入对模型误差有实质改善,分数就越高。也有人用 PLS 回归系数加权、或结合互信息做评分,但本质都是给变量一个可信度排序,然后截断。
每轮收缩之后,需要判断是否收敛。工程里的常见做法有两种:一是设定最大迭代次数,跑完即停;二是看当前候选集合是否连续两轮不再变化,同时最小 RMSECV 的下降幅度小于 0.1%,就提前终止。前者省心,后者省时间,两类判据我后面代码里都会体现。
2.3 iVISSA 相比 VISSA 的改进逻辑:收敛效率与稳定性
VISSA 的基本框架没有问题,但落地时会暴露两个痛点。一是抽样次数固定不变时,早期变量空间巨大,每个变量被抽到的频率低,得分波动很大,收缩方向不稳定。二是收缩率固定,收敛速度偏慢,迭代十几轮后还可能留下几百个变量,后续建模依然复杂。
iVISSA 的改进方向集中在动态性和历史累积。动态收缩是指收缩比例随迭代进度变化,前期少收缩、中期稳步收缩、后期精调,避免早期误杀重要变量。自适应抽样是把每轮抽样次数和当前变量空间大小挂钩,变量多时多抽,样本量少时适当降低,省下的算力留给更需要的轮次。我见过的一种实现还引入了加权历史累积,当前得分会叠加前一轮分数的衰减系数,相当于给投票做了平滑,能明显降低随机波动导致的误选。
这些改进的实际收益很直接。同一个玉米近红外数据集,我自己跑 VISSA 的时候需要重复三次才能拿到可靠结果,iVISSA 跑一轮就稳定了,而且筛出来的波段在 1200nm 到 1800nm 区间保持连续,和淀粉、水分的分子振动归属能一一对上。这种稳定性对工业定标来说比精度本身更有价值,因为筛选结果不可复现,模型验收就无从谈起。
3. 用 Python 从零实现 iVISSA:核心函数、完整代码与筛选结果解读
3.1 算法数据结构:候选变量集、权重向量与历史记录
写 iVISSA 之前,先把数据结构和轨迹想清楚。核心对象有三个:当前候选变量索引集合 active、每个变量的得分向量 score、以及每一轮的历史快照 history。
候选变量集合用 numpy 数组存索引,好处是切片、随机抽样都方便。得分向量在每一轮里重建,因为本轮只对当前候选变量打分,变量被裁掉后不再参与后续轮次。历史快照则是关键设计——很多人跑完算法只看最终结果,但 iVISSA 的价值恰恰在中间过程。如果你能画出每个变量得分随迭代次数的变化曲线,就能直观看到哪些变量是早期高分、后期被淘汰的“伪重要变量”,哪些是从头稳到尾的真特征。这个观察对后续调参和波段解释极其重要。
我还会额外用一个字典来记录每一轮的全局均方根误差和当前变量数,用于判断收敛。代码里我倾向于用 random.default_rng 生成随机数,保证可复现,而不是用 np.random 全局状态。
3.2 可直接运行的 iVISSA 核心代码(含 PLS 交叉验证)
下面这段代码是我在工程里经常改用的一个最小版本。为了便于读者直接复现,先构造一个简单的模拟光谱数据集作为演示对象。
import numpy as np from sklearn.cross_decomposition import PLSRegression from sklearn.model_selection import KFold from sklearn.metrics import mean_squared_error def make_demo_spectra(n_samples=80, n_vars=1500, seed=42): """构造一个带基线漂移和三个高斯峰的模拟光谱""" rng = np.random.default_rng(seed) X = np.zeros((n_samples, n_vars)) wave = np.arange(n_vars) for i in range(n_samples): baseline = 0.01 * i + 0.5 * np.sin(wave / 100) + rng.normal(0, 0.05, n_vars) peaks = np.zeros(n_vars) # 三个峰区域,中心波长分别为 500、700、900,宽度不同 for center, width, coef in [(500, 20, 1.2), (700, 35, 0.8), (900, 50, 1.0)]: peaks += coef * np.exp(-0.5 * ((wave - center) / width) ** 2) X[i] = baseline + peaks * rng.uniform(0.8, 1.5) # 目标变量只与 700 和 900 附近的两个峰相关,其余是干扰 y = X[:, 700] * 0.15 + X[:, 900] * 0.06 + rng.normal(0, 0.1, n_samples) return X, y, wave def pls_cv_rmse(X, y, ncomp=5, cv_fold=5, seed=42): """K 折交叉验证的 PLS 均方根误差""" ncomp = max(1, min(ncomp, X.shape[0] - 1, X.shape[1] - 1)) kf = KFold(n_splits=cv_fold, shuffle=True, random_state=seed) y_pred = np.zeros_like(y, dtype=float) for tr_idx, te_idx in kf.split(X): pls = PLSRegression(n_components=ncomp) pls.fit(X[tr_idx], y[tr_idx]) y_pred[te_idx] = pls.predict(X[te_idx]).ravel() return float(np.sqrt(mean_squared_error(y, y_pred)))make_demo_spectra构造了一个带有基线漂移和三个高斯峰的模拟光谱矩阵,其中只有两个峰与目标变量相关,第三个峰用于干扰测试。pls_cv_rmse封装了 K 折交叉验证的 PLS 评估,内部对 ncomp 做了下限和上限保护,避免子集变量数太少时 PLS 崩溃。
现在进入 iVISSA 主函数:
def ivissa(X, y, max_iter=15, n_subset=40, subset_ratio=0.7, shrink_ratio=0.8, ncomp=5, seed=42): rng = np.random.default_rng(seed) active = np.arange(X.shape[1]) # 初始候选变量 = 全谱 history = [] for it in range(max_iter): if len(active) <= 3: break # 每轮抽样子集大小与当前候选变量数挂钩 subset_size = max(2, int(len(active) * subset_ratio)) scores_dict = {v: [] for v in active} subset_rmses = [] for _ in range(n_subset): idx = rng.choice(active, size=subset_size, replace=False) rmse = pls_cv_rmse(X[:, idx], y, ncomp=ncomp, seed=seed) subset_rmses.append(rmse) for v in idx: scores_dict[v].append(rmse) # 平均误差越低,说明该变量在子集里越可能带来正向贡献 mean_scores = np.array([np.mean(scores_dict[v]) for v in active]) order = np.argsort(mean_scores) keep_num = max(3, int(len(active) * shrink_ratio)) new_active = active[order[:keep_num]] history.append({ "iteration": it, "active_size": len(new_active), "min_rmses": float(min(subset_rmses)), "mean_rmse": float(np.mean(subset_rmses)), "score": {int(v): float(np.mean(scores_dict[v])) for v in active}, }) active = new_active # 若候选集不再变化,提前收敛 if it > 0 and np.array_equal(active, history[-2]["active_vars"]): break return active, history主函数里最需要注意的点是order = np.argsort(mean_scores):因为平均交叉验证误差越低,说明变量越重要,所以是升序排列,然后保留前shrink_ratio比例的变量。每次迭代结束后,用np.array_equal对比当前候选集和上一轮候选集,如果完全一致就提前收敛,避免在局部稳定点空转。
代码里的历史记录history会在每一次迭代中保存当前候选变量数、本轮最小误差、平均误差和每个变量的得分。这个数据结构在后面分析和画变量得分轨迹时直接够用,不需要重新跑算法。
3.3 输出什么、怎么判定收敛、怎么恢复最终波段
运行模拟数据,输出三段关键信息,分别是最终保留的变量索引、收敛情况和交叉验证误差趋势。
X, y, wave = make_demo_spectra() selected, history = ivissa(X, y, max_iter=12, n_subset=40, subset_ratio=0.7, shrink_ratio=0.8, ncomp=5, seed=42) print("最终变量数:", len(selected)) print("选中的变量索引:", np.sort(selected)) print("迭代历史:") for h in history: print(f"{h['iteration']:>2d}轮 变量数 {h['active_size']:>5d} " f"最小RMSE {h['min_rmses']:.5f} 平均RMSE {h['mean_rmse']:.5f}")运行后你会看到变量数从 1500 逐轮收缩到几十个,最终停在一个既有数量优势又保持误差最低的集合上。如果两个峰区域的波长点保留比例明显高于第三峰,说明筛选方向正确;反之就得回头查抽样次数和收缩率。
判断收敛的标准不只看变量数是否稳定,还要看最小 RMSE 是否还在下降。很多时候变量数还在缓慢减少,但误差已经进入平台期,这时候继续迭代只会带来过拟合风险。我会建议把history里的min_rmses单独画出来,观察是否两轮之间下降幅度小于千分之一,达到就可以手动停掉。收敛后的候选变量索引就是最终筛选结果,直接用于后续 PLS 建模或导出为定标波段表。
4. iVISSA 的 4 个必调参数:经验边界与默认值速查
4.1 迭代次数:12 到 15 是常态,20 以上要警惕过拟合
迭代次数决定收缩的轮数。太少,变量空间还没收缩到稳定区域,结果接近随机抽样;太多,后期每个变量都被反复验证,选择慢慢偏向训练集噪声,放到外部验证一测就现原形。
我通常从 12 轮起步。变量数不多于 1000 时,12 轮足够把 0.8 收缩率的累积效果推到只剩 100 个以内。变量数超过 3000,比如某些高分辨光谱仪的数据,我会把迭代加到 20 轮。判断标准很简单:观察最后三轮的变量集合是否完全一致,如果一致但 RMSE 还在明显变化,加迭代次数;如果变量集合不再变化且 RMSE 稳定,就别再加了。
过拟合的信号也很明显。训练集 RMSECV 持续下降但外部验证 RMSEP 回弹,这个现象在 iVISSA 长迭代里特别常见。我的习惯是记录每一轮留下的变量数,如果第 15 轮和第 20 轮选出的变量重叠度超过 90%,把迭代数设在 15 而不是继续加码,稳定性和精度可以同时保住。
4.2 蒙特卡洛抽样次数和子集比例:全局探索与局部收缩的平衡
抽样次数 n_subset 控制每轮生成多少个子集。抽样次数太小,每个变量的得分方差大,筛选结果随机性高;太大,计算时间成倍增长。近红外全谱 1500 个变量时,我建议 30 到 50 次,变量降到 100 个以内后可以减少到 15 次,因为变量之间的共线性结构已经在上轮被破坏,继续大量抽样边际收益很低。
子集比例 subset_ratio 则决定每个抽样子集包含当前候选变量的比例。0.7 的默认值意味着每轮抽样时随机选择 70% 的当前变量。这个比例高,变量间协同作用看得全,但每个变量单独贡献的区分度会下降;比例低,子集探索更激进,但可能漏掉强相关变量组合。
共线性强的谱区,比如水分吸收峰附近的连续谱带,子集比例建议提到 0.8 甚至 0.85。因为分子振动峰本身带宽几十纳米,变量之间高度冗余,抽太少会误判“这整段都没用”,实际是该段任何一个变量都很有代表性。如果谱区分离度高、峰窄,0.6 到 0.7 就够用。
4.3 收缩率:0.7 到 0.85 之间的取舍依据
收缩率 shrink_ratio 是每轮保留变量的比例,0.8 意味着每轮裁掉 20%。这个参数直接决定收敛速度和误杀概率,也是我最常被问到、最像“玄学”的一个参数。
收缩率高于 0.85,收敛慢,计算开销大,但每轮只裁掉少量变量,安全性最高,适合对化学归属要求严格的项目。收缩率低于 0.6,收敛快,但非常危险。第一轮全谱 1500 个变量,随机组合里多数子集的主导因素其实是噪声,得分排序的尾部很可能包含真正有效的变量,一旦按比例裁掉就没有后悔药可吃,后续轮次里它永远不会再出现。
我的经验是首次运行用 0.8,看到历史记录里前两三轮裁掉的变量中有任何位于已知特征峰的变量,立刻把收缩率回调到 0.85。如果前几轮没有误杀,最后筛出的变量过多,再逐步降低收缩率。记住一个原则:收缩率宁可偏高,后期手动剔除多余变量,也不要偏低导致早期把关键谱区整段丢掉。
4.4 PLS 潜变量数:决定变量评分的公正性
PLS 潜变量数 ncomp 是 iVISSA 里最容易被忽略的一个参数,但它的作用比表面看上去大得多。ncomp 设置过小,模型欠拟合,所有变量子集的误差都偏高,变量间差异不明显,筛选退化成随机;ncomp 设置过大,模型过拟合,训练误差极低,误差的平台效应会掩盖单个变量的贡献。
我会先用全谱数据做一次 5 折交叉验证,扫一遍 ncomp 从 2 到 15 的 RMSECV 曲线,选 RMSECV 最低且潜变量数偏小的那个点作为固定值。后续整个 iVISSA 过程都用这个值,不要中途改动。因为变量数在收缩,同样 ncomp 下模型复杂度不同,中途改动会让前几轮得分和后几轮得分不具可比性。
补充一个细节:当候选变量减少到几十个时,ncomp 必须不能超过当前子集变量数和样本数的最小值。我的代码里用 min(ncomp, X.shape[0]-1, X.shape[1]-1) 做了下限保护,如果你换用自己的实现,这行避免会直接报错。
下表是参数边界速查:
| 参数 | 建议默认值 | 调参方向 | 常见错误 |
|---|---|---|---|
| max_iter | 12~15 | 变量多加到 20,稳定后别再涨 | 设置过大导致过拟合 |
| n_subset | 30~50 | 候选变量小于 100 时可降到 15 | 太少时选择结果波动大 |
| subset_ratio | 0.7 | 强共线性谱区提到 0.8+ | 过低漏掉相关变量组合 |
| shrink_ratio | 0.8 | 重跑时按误杀情况微调 0.7~0.85 | 低于 0.6 极其危险 |
| ncomp | 交叉验证决定 | 固定使用,不随迭代改动 | 中途变更导致得分不可比 |
5. iVISSA 避坑笔记:从变量消失到精度倒退的 5 个实战翻车现场
5.1 现象:同一数据集运行两次,选出来的变量差异很大
原因:随机种子未固定。iVISSA 依赖蒙特卡洛抽样,没设置随机种子时,每次运行产生的子集组合不同,得分排序在变量间差异不明显的区域会出现翻转。差异太大时,说明变量本身在模型里的可分性低,或者抽样次数不足以稳定排序。
解决:代码里用 np.random.default_rng(seed) 而不是 np.random 全局函数。固定 seed 后,同一数据同一参数必然得到同一结果。如果想要更稳健的结论,用多个种子各跑一遍,取交集作为最终波段,通常比单次结果更抗噪。
5.2 现象:经 iVISSA 筛选后的模型精度反而低于全谱
原因:变量收缩太狠,丢掉的变量里包含与目标相关的弱信号,而这些弱信号在全谱 PLS 中通过潜变量组合被间接利用了。另一个常见原因是 PLS 潜变量数没有随变量集调整——全谱 5 个主成分和筛选后 30 个变量的 5 个主成分含义完全不同,后者已经过拟合。
解决:筛完后单独做一次 PLS 潜变量数重扫,不要沿用全谱的 ncomp。如果重扫后 RMSEP 仍比全谱差,回到上一轮的历史记录,把被裁掉的变量按得分排序重新放回,逐步放回直到精度不低于全谱。这一步是 iVISSA 的兜底手段,很多人都忽略。
5.3 现象:迭代进行到一半,候选变量缩减为零
原因:代码里收缩率设置过低或子集比例过高,导致某轮抽样后所有变量的平均误差完全一致,排序后存活数量被截断成了一个很小的集合,下一轮子集大小不足、PLS 交叉验证报错直接中断。更隐蔽的原因是光谱矩阵里有 NaN 或零方差列,这些变量在 PLS 拟合时产生奇异矩阵。
解决:进入循环前先检查 X 是否有零方差列并剔除;在子集大小和收缩下限上加保护,比如 max(3, int(...)),避免候选变量数在早期降到 10 以下。如果必须处理缺失光谱,先做插值或剔除样本,不要带着 NaN 进 iVISSA。
5.4 现象:在非线性或强散射影响的谱区,筛选结果不如 CARS
原因:iVISSA 的核心评价指标是线性 PLS 的交叉验证误差,对非线性光谱响应天生不敏感。比如粉末样品的颗粒散射导致光谱基线漂移,或者浓度与吸光度不成线性关系时,误差排序会被散射主导,选出的波段集中在散射影响小的区域,而不是真正的化学信息区。
解决:先做光谱预处理(SNV、MSC、一阶导数),把散射和非线性趋势压掉再做筛选。如果预处理后仍不理想,可以对比 CARS 或 SPA 的结果,结合化学知识判断哪个更可信。iVISSA 不是万能方法,它在线性条件好的近红外透射数据上表现佳,在非线性强烈的拉曼定量场景里不如专门的非线性变量选择策略。
5.5 现象:光谱预处理在筛选前做和筛选后做,结果完全不同
原因:这不是 iVISSA 的问题,是变量选择对数据分布的敏感性差异。预处理改变变量间的相关结构和噪声幅度,先预处理再筛选,筛选结果对应的是处理后的光谱语义;先筛选再预处理,筛选过程面对的是原始噪声,很可能把本来干净的波段判为无用。
解决:统一流水线顺序:原始光谱 → 预处理 → iVISSA 筛选 → 筛选后数据重新预处理(如再次归一化)→ PLS 建模。不要在实验中途调整预处理方式,否则所有历史得分都失去参考价值。这个血泪经验是从一个药片近红外定标项目里得到的,从头到尾保持管线顺序一致,模型转移才真正稳得住。
6. iVISSA 在新数据集上的落地验证清单:一张参数表和三轮检验
6.1 快速开始的最小参数建议
新数据集落地时,不要上来就调参。用一组保守的默认参数跑通流程:max_iter=12,n_subset=40,subset_ratio=0.7,shrink_ratio=0.8,ncomp 用全谱 5 折交叉验证选定。先看筛选结果是否包含已知化学归属区间,再决定是否动参数。
最小验证流程是三段式:先看稳定性,再看对比,最后看外部验证。三轮都通过,这个筛选结果才敢放进定标模型里生产。
6.2 第一轮验证:稳定性检验
把 seed 从 1 换到 10,跑十次 iVISSA,统计两两之间变量集的 Jaccard 相似系数。相似系数平均值高于 0.8,说明筛选结果稳定;低于 0.5,说明当前参数下变量可分性太差不适合本数据。这时候优先提高 n_subset 到 60,再不行就降低收缩率到 0.85,而不是直接换算法。
6.3 第二轮验证:与全谱、CARS、SPA 的对比基准
在同一个交叉验证框架下,比较全谱 PLS、iVISSA 筛选、CARS 筛选和 SPA 筛选的 RMSECV 与 RMSEP。重点关注筛选后变量数相近时的精度差距。如果 iVISSA 精度略低但波段连续性明显更好,我通常会选 iVISSA——因为连续波段在模型转移和化学解释上的价值远大于零点几个 RMSE 百分比。
6.4 第三轮验证:外部预测与模型可转移性
用独立于训练集的验证样本做外部预测,这是所有筛选算法的终审。筛选结果稳定、内测精度良好,但外部 RMSEP 偏高,要回去检查是否过拟合或预处理泄漏。如果项目涉及两台仪器,把筛选后的波段分别从两仪器光谱里提取出来,做一次斜率截距校正,看预测偏差是否在可接受范围内。
我现在的习惯是每到一个新数据集,第一件事不是跑算法,而是先画光谱、看已知特征峰、确定预处理策略,再进 iVISSA。这个顺序帮我避开了大半的翻车场景。希望帮到你。
本文还有配套的精品资源,点击获取