如何用Matlab处理近红外光谱数据？5个关键步骤搞定定量分析

关注

原创最新推荐文章于 2026-09-15 16:39:29 发布 403 阅读 AI 写同款· GEO 优化 ›

MATLAB近红外光谱分析技术 1、主成分分析（PCA）、偏最小二乘（PLS）的基本原理（PCA与PLS的区别与联系；1、深度学习与传统机器学习的区别与联系（神经网络的隐含层数越多越好吗？2、BP神经网络参数的优化（隐含层神经元个数、学习率、初始权值和阈值等如何设置？3、MATLAB文件读写（mat、txt、xls、csv、jpg、wav、avi等格式）二。4、案例演示：一维卷积神经网络的MATLAB实现（基于卷积神经网络的近红外光谱建模）3、值得研究的若干问题（欠拟合与过拟合、泛化性能评价指标的设计、样本不平衡问题等）阅读详情

从光谱噪声到精准预测：用Matlab驾驭近红外数据的实战指南

如果你手头有一堆从近红外光谱仪导出的、看起来杂乱无章的数据曲线，心里盘算着如何从中提取出可靠的浓度信息，那么这篇文章就是为你准备的。近红外光谱分析，听起来是化学、制药或农业质检领域的专业课题，但其核心——从复杂信号中挖掘定量关系——却是许多科研和工程实践中都会遇到的通用挑战。Matlab，作为一款强大的数值计算与可视化平台，绝非仅仅是一个“画图工具”，它在处理这类高维、共线性严重的数据时，展现出的建模与优化能力，能让你从“看图说话”的定性分析，真正迈入“用数据决策”的定量世界。无论你是刚开始接触光谱数据分析的研究生，还是需要快速验证工艺模型的工程师，接下来的内容将抛开泛泛而谈的理论，直接切入五个环环相扣的关键操作步骤，并附上我实践中积累的代码片段和避坑心得，助你搭建起一套稳健、可复现的定量分析流程。

1. 数据导入与初诊：别让第一步就埋下隐患

拿到光谱数据，很多人会迫不及待地直接投入算法建模。但我的经验是，花在数据初诊上的每一分钟，都能在后续避免数小时的调试和错误结论。近红外数据通常以文本文件（.txt, .csv）或Excel表格的形式提供，每一列代表一个波长点，每一行代表一个样本，往往还会附带一列参考测量值（如浓度、水分含量等）。

1.1 稳健的数据读取策略

使用 readmatrix 或 readtable 函数通常比老旧的 xlsread 或 importdata 更可靠，尤其是在处理包含文本头或混合数据类型的文件时。我习惯先检查文件结构：

% 查看文件前几行，了解分隔符和标题行
opts = detectImportOptions('my_spectral_data.csv');
preview('my_spectral_data.csv', opts);

% 稳健地读取数据，指定数值数据范围
data = readmatrix('my_spectral_data.csv', 'Range', 'A2:K100'); % 假设数据从A2到K100
wavelength = data(:, 1); % 第一列为波长
spectra = data(:, 2:end-1); % 中间列为光谱强度
reference_value = data(:, end); % 最后一列为参考值（如浓度）

注意：务必确认波长向量是否单调递增。有些光谱仪导出的数据顺序可能是颠倒的，这会影响后续基于波长的处理（如卷积运算）。一个简单的检查：if ~all(diff(wavelength) > 0); wavelength = flip(wavelength); spectra = flip(spectra, 2); end

1.2 可视化初诊：发现隐藏的问题

在预处理之前，绘制原始光谱图至关重要。不要只画平均光谱，要把所有样本的光谱叠在一起看。

figure;
plot(wavelength, spectra);
xlabel('波长 (nm)');
ylabel('吸光度 (AU)');
title('原始近红外光谱（所有样本）');
grid on;

通过这张图，你需要敏锐地观察几个关键点：

基线漂移：曲线整体是否在Y轴方向上下移动？这通常由样品散射或仪器背景引起。
噪声水平：曲线是否毛刺很多？高频噪声会严重干扰弱特征峰的识别。
异常样本：是否有某一条光谱线明显偏离群体？它可能是测量失误或异常样品，需要在后续分析中决定是否剔除。

我曾遇到一个案例，由于比色皿洁净度问题，导致一个样本的基线异常抬高，如果直接参与建模，会成为一个强杠杆点，将整个模型“拉偏”。通过初诊图，我们迅速定位并复查了该样本。

2. 预处理“组合拳”：为模型准备好“干净”的输入

预处理不是简单地套用公式，其目标是增强信号中与分析物相关的信息，同时抑制或消除不相关的干扰。对于近红外光谱，我通常会按顺序考虑以下几步，但具体组合需根据初诊结果调整。

2.1 降噪：平滑与求导

高频随机噪声常用Savitzky-Golay滤波器（sgolayfilt）处理，它能在平滑的同时保留信号的原有形状，甚至直接计算一阶、二阶导数。导数光谱特别有助于消除基线漂移和增强重叠峰的分离。

% Savitzky-Golay平滑（窗口宽度11，2阶多项式）
order = 2;
framelen = 11;
spectra_smooth = sgolayfilt(spectra, order, framelen, [], 2); % 沿列（波长维度）滤波

% 计算一阶导数光谱（同时完成平滑和求导）
spectra_1der = sgolayfilt(spectra, order, framelen, 1, 2); % 第4个参数为求导阶数

参数选择是关键：窗口宽度（framelen）应大于噪声的周期但小于特征峰的宽度。一个实用的技巧是观察一段相对平坦区域的光谱，调整参数直到噪声被抑制而特征峰形无明显展宽。

2.2 基线校正与散射校正

对于明显的基线漂移，多元散射校正（MSC）或标准正态变量变换（SNV）比简单的去趋势（detrend）更有效，因为它们能同时校正乘性和加性散射效应。

% 标准正态变量变换 (SNV)
% 对每个样本的光谱进行中心化并除以自身标准差
spectra_snv = zeros(size(spectra));
for i = 1:size(spectra, 1)
    spec_mean = mean(spectra(i, :));
    spec_std = std(spectra(i, :));
    spectra_snv(i, :) = (spectra(i, :) - spec_mean) / spec_std;
end

2.3 数据标准化与划分

预处理最后一步，通常是将所有光谱进行均值中心化，有时还需标准化（除以标准差）。这能确保不同波长变量处于可比的数量级，对许多建模算法至关重要。

% 均值中心化
spectra_centered = spectra_snv - mean(spectra_snv, 1);

% 划分训练集与测试集（70%-30%是一种常见选择）
rng('default'); % 设置随机种子保证可重复性
cv = cvpartition(size(spectra_centered, 1), 'HoldOut', 0.3);
idx_train = training(cv);
idx_test = test(cv);

X_train = spectra_centered(idx_train, :);
y_train = reference_value(idx_train);
X_test = spectra_centered(idx_test, :);
y_test = reference_value(idx_test);

提示：务必在划分数据集之后再进行任何涉及全体数据统计量（如均值、标准差）的预处理操作（如SNV、标准化），然后用训练集的统计量去变换测试集，避免数据泄露（Data Leakage），这是保证模型评估结果无偏的关键。

3. 特征筛选与降维：从千维波长中抓住核心

近红外光谱数据点密集，相邻波长高度相关（共线性），直接使用全谱建模不仅计算量大，还容易引入噪声导致模型过拟合。因此，我们需要进行特征筛选或降维。

3.1 基于知识的波长区间选择

如果你对分析物的特征吸收峰位置有先验知识（例如，水分子的O-H键在特定波段有强吸收），可以直接选取这些感兴趣的波段（ROI）进行分析。这能极大提升模型的可解释性。

% 假设我们关注1450-1600 nm和1900-2100 nm两个特征区间
idx_band1 = find(wavelength >= 1450 & wavelength <= 1600);
idx_band2 = find(wavelength >= 1900 & wavelength <= 2100);
idx_combined = [idx_band1, idx_band2];

X_train_reduced = X_train(:, idx_combined);
X_test_reduced = X_test(:, idx_combined);

3.2 基于统计的降维：主成分分析（PCA）

PCA是无监督降维的利器，它能将高维数据投影到几个互不相关的主成分（PC）上，这些主成分包含了原始数据的大部分方差。

[coeff, score_train, latent, ~, explained] = pca(X_train);

% 查看前N个主成分的累计贡献率
cum_explained = cumsum(explained);
num_components = find(cum_explained >= 95, 1); % 找到累计贡献率达95%所需的主成分数
fprintf('保留前%d个主成分，可解释%.1f%%的方差。\n', num_components, cum_explained(num_components));

% 使用前num_components个主成分作为新特征
X_train_pca = score_train(:, 1:num_components);

% 将测试集投影到训练集得到的PCA空间上
score_test = (X_test - mean(X_train, 1)) * coeff(:, 1:num_components);

PCA得分图是异常样本检测的又一利器。绘制前两个主成分的散点图，远离样本群的点很可能就是异常值。

figure;
scatter(score_train(:,1), score_train(:,2), 'filled');
xlabel('PC1');
ylabel('PC2');
title('训练集PCA得分图（用于异常值检测）');

4. 建立定量模型：偏最小二乘回归的核心实践

对于光谱定量分析，偏最小二乘回归（PLSR）是当之无愧的“标准答案”。它巧妙地同时分解光谱矩阵X和浓度向量Y，寻找能最大程度解释Y变异的潜在变量（LV），天生适合处理高共线性数据。

4.1 模型训练与主成分数选择

在Matlab中，使用 plsregress 函数。最关键也是最容易被忽视的一步，是确定最佳潜在变量数。数量太少，模型欠拟合；太多，则过拟合。

% 设置最大潜在变量数进行交叉验证
max_lv = 10;
[Xloadings, Yloadings, Xscores, Yscores, beta, PLSstats] = plsregress(X_train, y_train, max_lv);

% 提取预测残差平方和（PRESS），用于交叉验证
PRESS = PLSstats.PRESS;

% 绘制PRESS随潜在变量数变化的曲线
figure;
plot(1:max_lv, PRESS, 'bo-', 'LineWidth', 1.5);
xlabel('潜在变量数');
ylabel('预测残差平方和 (PRESS)');
title('交叉验证PRESS图');
grid on;

% 常用规则：选择PRESS首次达到最小值或拐点处的LV数
% 也可以计算每个LV的Q2（预测解释方差）等指标辅助判断

在我的一个药品活性成分含量分析项目中，通过10折交叉验证的PRESS图，我们发现LV数从3增加到4时，PRESS下降已不显著，且第4个LV的回归系数向量开始出现大量无规律的微小波动（这是过拟合的典型信号），因此最终选择3个潜在变量。

4.2 模型系数与解释

得到模型后，不仅要看预测结果，还要理解模型本身。绘制回归系数图能告诉你哪些波长对预测贡献最大。

% beta包含截距项，所以回归系数从第二个元素开始
regression_coeff = beta(2:end);

figure;
plot(wavelength, regression_coeff);
xlabel('波长 (nm)');
ylabel('回归系数');
title('PLSR模型回归系数谱');
grid on;
hold on;
% 可以在图上标出系数绝对值大的波长点
[~, idx_max] = max(abs(regression_coeff));
plot(wavelength(idx_max), regression_coeff(idx_max), 'r*', 'MarkerSize', 10);

系数绝对值大的波长点，意味着该处的光谱变化对预测浓度变化最敏感。你可以将此图与物质的已知特征吸收峰对比，检验模型是否具有物理化学意义上的合理性。

5. 模型评估与部署：跨越从“拟合”到“预测”的鸿沟

一个在训练集上表现完美的模型，可能在新样本面前一败涂地。因此，我们必须用未参与训练的数据来严格评估其泛化能力。

5.1 关键性能指标计算

不要只看R²，尤其是训练集的R²。对于定量模型，我习惯同时计算以下几个指标：

% 对测试集进行预测
y_pred_test = [ones(size(X_test,1),1), X_test] * beta;

% 计算性能指标
n_test = length(y_test);
RMSE_test = sqrt(mean((y_test - y_pred_test).^2)); % 均方根误差
R2_test = 1 - sum((y_test - y_pred_test).^2) / sum((y_test - mean(y_test)).^2); % 决定系数
RPD = std(y_test) / RMSE_test; % 残差预测偏差，>2.5通常认为模型预测能力较好

fprintf('测试集性能:\n');
fprintf('RMSE: %.4f\n', RMSE_test);
fprintf('R²: %.4f\n', R2_test);
fprintf('RPD: %.4f\n', RPD);

RPD是一个非常有用的指标，它衡量了模型预测误差相对于样本本身自然变异的大小。在农产品品质检测中，RPD > 2.0通常被认为模型可用于粗略筛选，而 > 3.0则可用于精确定量。

5.2 可视化诊断：发现系统误差

绘制预测值 vs. 真实值的散点图和残差图，能直观判断模型是否存在系统偏差。

figure;
subplot(1,2,1);
plot(y_test, y_pred_test, 'ko', 'MarkerFaceColor', 'b');
hold on;
plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], 'r--', 'LineWidth', 1); % 绘制y=x参考线
xlabel('参考测量值');
ylabel('模型预测值');
title('预测值 vs. 真实值');
legend('测试样本', '理想拟合线', 'Location', 'best');
grid on;

subplot(1,2,2);
residuals = y_test - y_pred_test;
plot(y_pred_test, residuals, 'ks', 'MarkerFaceColor', 'r');
hold on;
plot([min(y_pred_test), max(y_pred_test)], [0, 0], 'k-'); % 绘制y=0参考线
xlabel('预测值');
ylabel('残差');
title('残差图');
grid on;

理想的残差图应该是点随机、均匀地分布在零线两侧，无明显趋势。如果残差呈现喇叭形（方差齐性问题）或曲线趋势，说明模型可能遗漏了某些非线性信息或关键变量。

5.3 模型保存与应用

最后，将训练好的模型参数（beta、预处理所需的均值、标准差、PCA系数等）保存下来，用于未来新样本的预测。这里建议封装成一个预测函数：

function y_pred = my_pls_predictor(new_spectrum, wavelength, model_params)
    % model_params 是一个结构体，包含beta, train_mean, train_std, pca_coeff等
    % 1. 预处理（使用训练时保存的参数）
    new_spec_processed = (new_spectrum - model_params.train_mean) ./ model_params.train_std;
    % 2. 降维（如果用了PCA）
    new_scores = new_spec_processed * model_params.pca_coeff;
    % 3. 预测
    X_new = [1, new_scores]; % 添加截距项
    y_pred = X_new * model_params.beta;
end

记住，这个模型的有效性边界受限于你训练样本的浓度范围、基质类型和测量条件。当遇到全新的样本类型时，模型可能需要更新或重建。定期用新的标准样本验证模型性能，是保证其长期可靠性的不二法门。