系列文章目录

近红外光谱分析技术属于交叉领域，需要化学、计算机科学、生物科学等多领域的合作。为此，在（北京邮电大学杨辉华老师团队）指导下，近期准备开源传统的PLS，SVM，ANN，RF等经典算和SG，MSC，一阶导，二阶导等预处理以及GA等波长选择算法以及CNN、AE等最新深度学习算法，以帮助其他专业的更容易建立具有良好预测能力和鲁棒性的近红外光谱模型。

---

---

前言

NIRS是介于可见光和中红外光之间的电磁波，其波长范围为（1100∼2526 nm。 由于近红外光谱区与有机分子中含氢基团（OH、NH、CH、SH）振动的合频和 各级倍频的吸收区一致，通过扫描样品的近红外光谱，可以得到样品中有机分子含氢 基团的特征信息，常被作为获取样本信息的一种有效的载体。 基于NIRS的检测方法具有方便、高效、准确、成本低、可现场检测、不 破坏样品等优势，被广泛应用于各类检测领域。但 近红外光谱存在谱带宽、重叠较严重、吸收信号弱、信息解析复杂等问题，与常用的 化学分析方法不同，仅能作为一种间接测量方法，无法直接分析出被测样本的含量或 类别，它依赖于化学计量学方法，在样品待测属性值与近红外光谱数据之间建立一个 关联模型(或称校正模型，Calibration Model) ，再通过模型对未知样品的近红外光谱 进行预测来得到各性质成分的预测值。现有近红外建模方法主要为经典建模 （预处理+波长筛选进行特征降维和突出，再通过pls、svm算法进行建模）以及深度学习方法（端到端的建模，对预处理、波长选择等依赖性很低）

---

本篇主要讲述基于matlab语言版本的光谱预处理方法,

一、预处理算法

%% 标准化Autoscales
function [ax,mx,stdx] = auto(x)
% Autoscales matrix to zero mean and unit variance
%
% [ax,mx,stdx] = auto(x)
%
% input:
% x 	data to autoscale
%
% output:
% ax	autoscaled data
% mx	means of data
% stdx	stantard deviations of data
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil



[m,n] = size(x);
mx    = mean(x);
stdx  = std(x);
ax    = (x-mx(ones(m,1),:))./stdx(ones(m,1),:);

%# 数据中心化(Mean centering)
function [mcx,mx] = center(x)
% Mean center scales matrix to zero mean
%
% [mcx,mx] = center(x)
%
% input:
% x 	data to mean center
%
% output:
% ax	mean center data
% mx	means of data
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil

[m,n] = size(x);
mx    = mean(x);
mcx   = (x-mx(ones(m,1),:));

%#归一化Normalize
function [nx] = normaliz(x)
% Normalize matrix rows dividing by its norm
% 
% [nx] = normaliz(x)
%
% input:
% x (samples x variables)   data to normalize
%
% output:
% nx (samples x variables)  normalized data
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil

[m,n]=size(x);
nx=x;
nm=zeros(m,1);
for i = 1:m
nm(i)=norm(nx(i,:));
nx(i,:)=nx(i,:)/nm(i);
end
% z=[1,2,3;3,4,5];
% [m,n]=size(z);
% nx=z;
% nm=zeros(m,1);
% for i = 1:m
% nm(i)=norm(nx(i,:));
% nx(i,:)=nx(i,:)/nm(i);

%% 多元散射校正MSC
function [x_msc]=msc(x,xref)
% Multiplicative Scatter Correction
%
% [x_msc]=msc(x,xref)
%
% input 
% x (samples x variables)      spectra to correct
% xref (1 x variables)         reference spectra (in general mean(x) is used)
%
% Output
% x_msc (samples x variables)  corrected spectra  
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil

[m n]=size(x);
rs=xref;cw=ones(1,n);
mz=[];mz=[mz ones(1,n)'];mz=[mz rs'];
[mm,nm]=size(mz);
wmz=mz.*(cw'*ones(1,nm));
wz=x.*(ones(m,1)*cw);
z=wmz'*wmz;
[u,s,v]=svd(z);sd=diag(s)'; 
cn=10^12;
ms=sd(1)/sqrt(cn);
cs=max(sd,ms ); 
cz=u*(diag(cs))*v';  
zi=inv(cz);
b=zi*wmz'*wz';B=b';
x_msc=x; 
p=B(:,1);x_msc=x_msc-(p*ones(1,mm));
p=B(:,2);x_msc=x_msc./(p*ones(mm,1)');

%# 标准正态变换
function [x_snv] = snv(x)
% Standard Normal Variate
%
% [x_snv] = snv(x) 
%
% input:
% x (samples x variables) data to preprocess
%
% output:
% x_snv (samples x variables) preprocessed data
%
% By Cleiton A. Nunes
% UFLA,MG,Brazil

[m,n]=size(x);
rmean=mean(x,2);
dr=x-repmat(rmean,1,n);
x_snv=dr./repmat(sqrt(sum(dr.^2,2)/(n-1)),1,n);

%# 移动平均平滑
%# function [mdata]=nirmaf(data,window)
%#
%#  AIM:   移动窗口平滑光谱矩阵，用于光谱预处理
%#
%#  INPUT:  data        m×n的矩阵，m个光谱，n个变量
%#          window      光谱窗口宽度，必须是奇数(默认: 3)
%# 
%#  OUTPUT: mdata       平滑过的光谱，m个光谱，n个变量
%#
%#  AUTHOR: 王毅 
%#  EMAIL:  wang727yi@hotmail.com
%#  VERSION:1.0 (15/03/2009)

function [mdata]=nirmaf(data,window)

[m,n] = size(data);	
mdata = zeros(m,n);	

if nargin == 1
	window = 3;	    % 默认窗口宽度
elseif round(window/2) == window/2
	error('光谱窗口宽度必须是奇数')
end

wcenter = floor(window/2);
extdata = [zeros(m,wcenter) data zeros(m,wcenter)];	% 延展原始光谱矩阵


%# 延展后的光谱矩阵曲线拟合
for k = 1:m
	bstart = polyfit(wcenter+1:wcenter+4,extdata(k,wcenter+1:wcenter+4),2);		% 左端曲线拟合
	bend = polyfit(n-4+wcenter:n+wcenter,extdata(k,n-4+wcenter:n+wcenter),2);	% 右端曲线拟合
	extdata(k,1:wcenter) = polyval(bstart,1:wcenter);				% 返回左端曲线拟合的值
	extdata(k,n+wcenter+1:n+window-1) = polyval(bend,n+wcenter+1:n+window-1);	% 返回右端曲线拟合的值
end



%# 均值平滑
for i = 1:n
	mdata(:,i) = mean(extdata(:,i:i+window-1)')';	
end

%# Savitzky-Golay平滑滤波

%%%%问题：设置N为3，F为5，A97中的每一行（一个样本）的到的平滑数据、一阶导数和二阶导数均是是一行104个数据；有97个样本，得到的数据应该是97*104，但按以下程序运行得到的平滑数据和一阶导数，二阶导数为一行数据9897个，图也不对？？？

%光谱预处理，包括Savitzky-Golay平滑和Savitzky-Golay求导
clear
clc
load A97 %导入光谱,第一列为编号，最后一列为浓度值
dx=1;
x=215:321; %波长范围是215nm-321nm
a=length(x); %x的数据个数
k=0:96;
z=A97(:,109); %97个样本对应的浓度值 （这一句对光谱平滑、求导没有多大关系）
N=input('请输入拟合次数');
F=input('设置窗口参数(奇数)');

[b,g]=sgolay(N,F);
Halfwin=((F+1)/2)-1;
hold on            %这是为了保存每次计算的图
for k=k+1
hold on 
for n=(F+1)/2:a-(F+1)/2
y=A97(k,2:108); %第k个样本 215nm-321nm应的吸光度对
%SG平滑
SG0(n,k)=dot(g(:,1),y(n-Halfwin:n+Halfwin));   %红色是之前没有添加的
%一阶差分求导
SG1(n,k)=dot(g(:,2),y(n-Halfwin:n+Halfwin));
%二阶差分求导
SG2(n,k)=2*dot(g(:,3)',y(n-Halfwin:n+Halfwin)); 

end
SG1=SG1/dx;
SG2=SG2/(dx*dx);

subplot(3,1,1);
plot(SG0)
legend('S-G Smooth');

subplot(3,1,2);
plot(SG1)
legend('S-G Smooth 1st derivative');

subplot(3,1,3);
plot(SG2)
legend('S-G Smooth 2st derivativa');
end

代码如下（示例）：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
import  ssl
ssl._create_default_https_context = ssl._create_unverified_context

二.使用教程

1.安装matlab

受matlab版权保护，matlab安装教程自行查找

2.读入数据、预处理以及展示

clc;
clear;
% close all;
%%  打开文件
%%%%%%%%%
CA1 = csvread('D:\ProgramData\Inception\A1.csv');
Data = CA1;
%%%%%%%%%

 %% 作图，原始图
figure(1);
title('光谱图');
xlabel('Wavenumber(cm-1)');
ylabel('Absorbance');
hold on;
plot(Data');

%% 数据处理

Absorbance = Data;
[Absorbance_m,Absorbance_n]=size(Absorbance);
Absorbance_mean = mean(Absorbance);
% %% 吸光度波数
% Absorbance=data; %得到吸光度
% % Absorbance=Absorbance;
% [Absorbance_m,Absorbance_n]=size(Absorbance);
% Wavenumber=data(:,1:2:end); %得到波数
% % Wavenumber=1:141;
% Absorbance_mean=mean(Absorbance);%每个样本吸光度平均

% %% 多元散射校正MSC
% Absorbance_msc=msc(Absorbance,Absorbance_mean);
% figure(2);
% plot(Wavenumber,Absorbance_msc);
% %set(gca,'XDir','reverse'); % 横坐标从大到小
% title('多元散射校正MSC');
% xlabel('Wavenumber(cm-1)');
% ylabel('Absorbance');
%% 多元散射校正MSC
msc_file_name =  'D:\DsekTop\MSCA1.csv';
Absorbance_msc=msc(Absorbance,Absorbance_mean);
csvwrite(msc_file_name,Absorbance_msc,0,0);
figure(2);
title('多元散射校正MSC');
xlabel('Wavenumber(cm-1)');
ylabel('Absorbance');
hold on;
plot(Absorbance_msc');

3.结果(以msc为例)

原始光谱

MSC预处理后光谱

总结

完整代码可从获得GitHub仓库
代码仅供学术使用，如有问题，联系方式：QQ：1427950662，微信：Fu_siry