Matlab预处理近红外光谱

Matlab预处理近红外光谱

1. 微分处理光谱 一阶微分公式：

𝑥𝑖,1𝑠𝑡

二阶微分公式：

𝑥𝑖,2𝑠𝑡

𝑥𝑖+𝑥𝑖+𝑔−2𝑥𝑖=

𝑔2𝑥𝑖−𝑥𝑖+𝑔= 𝑔式中，𝑥𝑖是谱图数据中i波数下的透过率，g为窗口宽度 一阶导数MATLAB代码：

X1st=diff(X,1);%X 为输入光谱矩阵，X1st 为输出一阶导数光谱矩阵 二阶导数MATLAB代码：

X2st=diff(X,2); %X 为输入光谱矩阵，X2st 为输出二阶导数光谱矩阵 2. 多元散射校正处理谱图 （1）计算平均光谱：

𝑥̅𝑖,𝑗

（2）一元线性回归：

̅+𝑏𝑖 𝑋𝑖=𝑚𝑖𝑋

（3）多元散射校正：

𝑋𝑖(𝑀𝑆𝐶)

(𝑋𝑖−𝑏𝑖)

=

𝑚𝑖∑𝑛𝑖=1𝑋𝑖,𝑗=

𝑛公式中X表示n×p维定标光谱数据矩阵，n为样品数，p为波点数。̅表示所有样品的原始近红外光谱在各个波长点处求平均值所得到𝑋

1

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.

的平均光谱矢量，𝑋𝑖是1×p维矩阵，表示单个样品光谱矢量，𝑚𝑖和𝑏𝑖分别表示各样品近红外光谱𝑋𝑖与平均光谱X进行一元线性回归后得到的相对偏移系数和平移变量。I代表第i个样本，j代表第j个波数。 MATLAB代码如下： me=mean(X); [m,~]=size(X); for i=1:m,

p=polyfit(me,X(i,:),1);

Xmsc(i,:)=(X(i,:)- p(2)*ones(1,n))./(p(1)*ones(1,n)); End

代码中输入光谱矩阵为 X，输出光谱矩阵为 Xmsc。 3. 标准正态变量变换（SNV）

𝑋𝑖,𝑆𝑁𝑉=

𝑥𝑖,𝑘−𝑥𝑖∑𝑚(𝑥𝑖,𝑘−𝑥𝑖)2𝑘=𝑖√

(𝑚−1) 式中𝑥𝑖是第i样品光谱的平均值，k=1,2,…,m。m为波长点数；i=1,2,..,n。n为校正样品数；𝑋𝑖,𝑆𝑁𝑉是变换后的光谱。 SNV的MATLAB代码如下： [~,n]=size(X); rmean=mean(X,2); dr=X-repmat(rmean,1,n);

Xsnv=dr./repmat(sqrt(sum(dr.^2,2)/(n-1)),1,n);

代码中 X 为光谱输入矩阵，Xsnv 为 SNV 处理后光谱输出矩阵。

2

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.

4. 小波变换处理谱图

傅里叶变换是小波变换的基础，傅里叶变换是一种分析信号的方法。傅里叶变换的基本公式为：

+∞

F(ω)=∫𝑓(𝑡)𝑒−𝑖𝜔𝑡𝑑𝑡

−∞

公式中 f (t )是输入信号,在这里可看作光谱信号，t是时域信号，可以看作是波数，F (ω)是在信号频谱。

傅里叶变换原理表明：任何连续测量的时序或信号，都可以表示为不同频率正弦信号的无限叠加。傅里叶变换将原来难以处理的时域信号转换成了易于分析的频域信号，在频域表示的信号图称为信号的频谱。在频谱中可以看到时域谱图中看不到的特征，处理这些特征后利用傅里叶逆变换可以将这些频域信号转换成时域信号。

小波变换是傅里叶变换的延伸。傅里叶变换有很多缺点：不能求解变系数线性微分方程，只适用于平移不变的线性系统，不能反映信号在局部时间范围内的信息。小波变换继承了傅里叶变换的长处，同时解决了傅里叶变换的一些缺点。利用小波变换能够“放大”信号局部时间、频率范围内的信息，消除信号噪声。连续小波变换基本公式为：

∗

Υ(s,τ)=∫𝑓(𝑡)𝜑𝑠,𝜏(𝑡)𝑑𝑡

小波逆变换为：

𝑓(𝑡)=∬Υ(𝑠,𝜏)𝜑𝑠,𝜏(𝑡)𝑑𝑠𝑑𝜏

3

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.

公式中*表示共轭，𝜑𝑠,𝜏表示小波基函数。小波基函数是由一个基本小波φ(𝑡)放缩、平移形成。

𝑡−𝜏

𝜑𝑠,𝜏(𝑡)=𝜑()

𝑠√𝑠本研究使用 MATLAB 的 wavelet 工具箱里的连续小波变换函数，MATLAB 代码如下：

Xcwt(i,:)=CWT(X(i,:),scale,basis);

代码中 X(i,:)为第i个样本原始光谱矩阵，scale是小波分解尺度，取值范围1到正无穷，basis是小波基函数，Xcwt(i,:)是第i个样本光谱小波变换后的小波系数。 5. 提取图谱主成分

主成分分析的实质是对原坐标系进行平移和转换，使得新坐标原点与样本点集合的重心重合，新坐标的第一主轴与数据变异的最大方向对应，新坐标的第二主轴与第一主轴标准正交，并且对应于数据变异的第二大方向。其基本步骤如下：

（1）对数据进行标准化处理，得到矩阵X。也可以不处理，但是标准化处理可以使数据计算方便。

（2）计算矩阵X的协方差矩阵V。这时V又是X的相关系数矩阵。

（3）求V的前n个特征值Υ1>Υ2>…>Υ𝑛，以及特征值对应的特征向量𝛼1>𝛼2>…>𝛼𝑚要求它们相互正交。

（4）求第h主成分𝐹ℎ，有：

1

4

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.

𝑝

𝐹ℎ=𝑋𝛼ℎ=∑𝛼ℎ𝑗𝑥𝑗

𝑗=1

式中：𝛼ℎ𝑗是主轴𝛼ℎ的第j个分量。所以主成分𝐹ℎ是原变量𝑥1,𝑥2, ,…, 𝑥𝑝的线性组合，组合系数恰好为𝛼ℎ𝑗。从这个角度，又可以说𝐹ℎ是一个新的综合变量。

MATLAB代码如下：

[~,XSelected_SCORE,~]=princomp(X); f1=XSelected_SCORE(:,1); f2=XSelected_SCORE(:,2); . .

fn=XSelected_SCORE(:,n); Xprincompn=[f1,f2,…,fn];

X 为输入光谱矩阵，XSelected_SCORE 是提取的全部主成分得分，f1,f2,…,fn 是第 1，2，…，n 个主成分。Xprincompn 是提取前 n 个主成分。

6. 相关系数法度量相似度

相关系数（通常用 R 表示）表明变量之间线性相关程度。在这里计算两个光谱谱图之间的相关系数来评价相似度，相关系数小于等于 1，越接近 1 越相似。公式如下：

R=

∑𝑛̅𝑘)𝑖=1(𝑥𝑗𝑖−𝑥̅𝑗)(𝑦𝑘𝑖−𝑦

𝑛(√∑𝑛∑−𝑥̅̅𝑘)2(𝑥)𝑗𝑖𝑗𝑖=1𝑖=1𝑦𝑘𝑖−𝑦

2 5

文档来源为:从网络收集整理.word版本可编辑.欢迎下载支持.

R为相关系数，𝑥𝑗𝑖，𝑦𝑘𝑖是𝑗，𝑘两个谱图第𝑖点的透射率，𝑥̅𝑗和𝑦̅𝑘是𝑗，𝑘两个谱图在各个波数的均值，n根据谱图选取。MATLAB 实现代码如下：

r= corrcoef(X1, X2); R=r(1,2);

代码中 X1,X2 代表两个光谱谱图数据，R 为要求的相关系数。corrcoef 是 MATLAB 计算相关系数矩阵的内置函数。

6