CN102841063A - 一种基于光谱技术的生物炭溯源鉴别方法 - Google Patents

Abstract

本发明公开了一种基于光谱技术的生物炭快速溯源鉴别方法：（1）建立生物炭建模光谱库；（2）对生物炭建模光谱库进行主成分分析，按主成分方差贡献率大小，取前3个主成分；（3）根据前3个主成分得分值，计算生物炭品种三维中心坐标；（4）待鉴别样本光谱扫描；（5）对待鉴别样本光谱作主成分分析，按主成分方差贡献率大小，提取前3个主成分，并计算与每一品种中心坐标之间的马氏距离；（6）根据待鉴别样本到每一品种中心坐标之间的最小马氏距离，将待鉴别样本归类到相应的生物炭品种，从而完成生物炭溯源鉴别。本发明属于一种有监督机器学习方法，通过利用光学间接测试和数值处理手段，达到生物炭快速溯源鉴别目的，可在实验室或工作现场应用。

Description

一种基于光谱技术的生物炭溯源鉴别方法

技术领域

本发明涉及一种基于光谱技术的生物炭溯源鉴别方法。

背景技术

生物炭是生物质在缺氧条件下不完全燃烧热裂解后形成的产物，可溶性极低，孔隙度和比表面积大，吸附力、抗氧化力和抗生物分解能力强。生物炭在土壤肥力改良、土壤固碳增汇减排以及受污染环境修复等方面具有显著作用。生物炭理化性能不仅受炭化过程中温度、时间、强度等工艺参数的影响，还与生物质材料有关。据Xu等报道（Xu R,Xiao S,Yuan J,et al.Adsorption of methylviolet from aqueous solutions by the biochars derived from crop residues.Bioresource Technology,2011,102(22):10293-10298.），以蓖麻秸秆、花生秸秆、大豆秸秆和稻壳为原料的生物炭对水溶液中甲基紫（俗称龙胆紫，一种染色剂）的吸附能力存在差别。生物炭使用部门往往也需要对生物炭进行溯源鉴别，从中筛选出具有特定材料性能的生物炭。各种来源的生物质经过炭化工艺加工制备成生物炭后，表观上都呈深黑色，肉眼不易区分甄别。而传统理化测试步骤繁琐，测试时间长，不适合用于生物炭的快速鉴别。因此有必要开展生物炭快速溯源鉴别方法研究。

近年来，红外光谱技术因其测试快速、重复性好、测试精度较高、易于现场使用、免试剂等优点而得到广泛关注。国内外学者对橘色合欢树（Chia C H,Gong B,Joseph S D,et al.Imaging of Mineral-enriched biochar by FTIR,Raman andSEM-EDX.Vibrational Spectroscopy.2012.In Press.）、山核桃壳(Novak J M,Busscher W J,Watts D W，et al.Short-term CO₂ mineralization after additions ofbiochar and switchgrass to a Typic Kandiudult.Geoderma.2010,154(3-4):281-288.)、柳枝稷(Kumar S,Loganathan V A,Gupta R B,et al.An Assessment ofU(VI)removal from groundwater using biochar produced from hydrothermalcarbonization.Journal of Environmental Management.2011,92(10):2504-2512.)、甜菜根(Dong X,Ma L Q,Li Y.Characteristics and mechanisms of hexavalentchromium removal by biochar from sugar beet tailing.Journal of HazardousMaterials.2011,190(1-3):909-915.)等生物炭进行过红外光谱特征分析。但在利用光谱技术对生物炭快速溯源鉴别研究上，还未见文献报道。

发明内容

为了达到简便快速鉴别生物炭品种的目的，本发明提供了一种基于光谱技术的生物炭溯源鉴别方法。

（1）建立生物炭品种光谱库：收集已知不同源质品种的生物炭材料，每种材料5-10个样本，利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，采样间隔1nm,共751个波长点，记录各波长点的光强值，并作基线偏置校正，组建生物炭建模光谱库。

基线偏置校正公式如下：

F(x)＝f(x)-minf(X)

式中，f(x)是波长x的光强，X表示所有波长变量，F(x)表示波长x校正后的光强值，min(X)表示所有波长点的光强值的最小值；。

生物炭建模光谱库由所有样本的光谱记录组成。每一样本的光谱记录包括生物炭品种名称（1字段）、不同波长点校正后的光强值（751字段），通常还可包括归类代码（1字段）。

（2）对生物炭建模光谱库进行主成分分析，按主成分方差贡献率大小，取所述的贡献率最大的前3个主成分；

所述主成分分析的计算步骤如下：

对于n个样本，每个样本用用p个变量来描述研究对象，分别用X₁、X₂…X_p表示，X_p表示第p个波长点所有样本的光强值按顺序组成的向量，这p个变量构成的p维随机向量为X=(X₁，X₂…X_p)^t。设随机向量X的均值为μ，相关系数矩阵为R，对X进行线性变换，考虑原始变量的线性组合：

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=u_{21}{X}_1+u_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ \·\\ \·\\ \·\\ Z_p=u_{p1}{X}_1+u_{p2}{X}_2+\·\·\·+u_{\mathrm{pp}}{X}_p\end{array}\right.$

主成分是不相关的线性组合Z₁，Z₂…Z_p，并且Z₁是X₁，X₂…X_p的线性组合中方差最大者，Z₂是与Z₁不相关的线性组合中方差最大者，Z_p是与Z₁，Z₂…Z_p-1都不相关的线性组合中方差最大者。

计算相关系数矩阵R，R可以反映数据之间相关程度，

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1p}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

其中R_ij（i，j＝1,2…p）为原始变量X_i与X_j的相关系数。由于R为实对称矩阵（即R_ij=R_ji），只需计算其上三角元素或下三角元素即可，R_ij的计算公式为：

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}^2}}$

n表示所有样本的样本数；

根据相关系数矩阵R解特征方程

|λE-R|=0

求出特征值λ_i（i=1,2…p）。因为R是正定矩阵，所以其特征值λ_i都为正数，将λ_i按大小顺序排序，即

λ₁≥λ₂≥…≥λ_i≥0

分别求出对应于特征值λ_i的特征向量e_i(i=1,2,…,p)，要求||e_i||=1，即

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

其中e_ij表示向量e_i的第j个分量;

特征值是各主成分的方差，它的大小反映了各个主成分的影响力。主成分Z_i的贡献率W_i的计算公式为

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

根据主成分方差贡献率大小，取所述的贡献率最大的前3个主成分Z₁，Z₂，Z₃。

根据大小顺序排前三的λ₁、λ₂、λ₃，按照下式计算三个主成分的主成分载荷

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}(i,j=\mathrm{1,2},\·\·\·,p)$

得到各主成分的载荷以后，可以按照下式计算各主成分的得分

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=l_{21}{X}_1+l_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

则Z₁、Z₂、Z₃表示n个样本的3个新变量指标；[Z₁ Z₂ Z₃]是一个n×3的新变量矩阵；

（3）根据前3个主成分Z₁，Z₂，Z₃，将生物炭同一品种不同样本的数据取均值，计算生物炭每一品种的三维中心坐标。

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m},(k=\mathrm{1,2,3}),$

k表示第k种主成分，也就是第k维；

生物炭建模光谱库第i品种样本光谱第k维主成分平均值；

x_i，j,k:生物炭建模光谱库第i品种第j样本光谱第k维主成分得分值；

m:生物炭建模光谱库第i品种的样本数；

则生物碳第i品种的三维中心坐标即为

（4）待鉴别样本光谱扫描：利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，记录与步骤（1）相同的751个波长点的光强值，并作基线偏置校正，方法同步骤（1），组建待鉴别样本光谱库。

（5）对待鉴别样本光谱库作主成分分析，按主成分贡献率大小，提取贡献率最大的前3个主成分，方法同步骤（2），计算各主成分的得分值，从而得到每一待鉴别样本的三维中心坐标；。并计算待鉴别样本与每一品种的三维中心坐标之间的马氏距离。对任待鉴别一样本x=(x₁,x₂,x₃)^t，第i品种的中心坐标

马氏距离公式为：

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}$

式中，S^-1是S的逆矩阵，S是生物炭建模光谱库中第i品种所有样本的主成分得分矩阵的协方差矩阵,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

其中S_kl（k，l=1,2,3）为第i品种中所有样本第k维的主成分得分值向量Z_k与所有样本第l维的主成分得分值向量Z_l之间的协方差，且S_kl=S_lk，其计算公式为：

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

m:生物炭建模光谱库第i品种的样本数；

（6）根据待鉴别样本到每一品种中心坐标之间的最小马氏距离，将待鉴别样本归类到相应的生物炭品种，从而完成生物炭溯源鉴别。

本发明所述不同源质品种的生物炭材料可以来自农林废弃物、生活垃圾、炭基肥、畜禽***物等。

具体地，以下14种生物炭的光谱数据按本发明方法检测建模，其三维主成分中心坐标如下表所示:

总体来说，本发明是根据生物炭光谱特征对未知样本进行鉴别测试，是一种有监督机器学习方法。本方法利用光学间接测试和数值处理手段，达到生物炭快速溯源鉴别目的。

本发明较常规理化分析技术有多项优点：

（1）快速检测。单个生物炭样本测试时间仅为几秒钟。

（2）现场使用。光谱仪携带方便，不受场地限制。

（3）批量测试。能够对大量样本进行分析检测，免试剂，无污染。

附图说明

图1为基于光谱技术的生物炭快速溯源鉴别的方法示意图。

图2实施例中14种生物炭的典型发射光谱图。

具体实施方式

收集到以下14种源质生物炭，其中农业废弃物7种（玉米、砻糠、花生壳、山核桃壳、稻杆、麦杆），林业废弃物4种（毛竹、食用菌菌渣、木屑、柠条），生活垃圾1种，生物炭复合肥1种、养殖***物1种（猪粪）。为方便起见，表1列出了这些生物炭中英文名称及缩略词。

表1

生物炭用研钵磨细，获得颗粒直径小于0.1mm样本。每种生物炭制备10个样本，共140个样本用来分析测试。

生物炭样本的测量装置包括海洋光学公司生产的Maya2000Pro光谱仪、10W卤素灯光源、传输光纤、样本架、计算机、电源等组成。波长范围选用350-1100波段。光纤探头垂直放置于样本上方1cm左右处，光源光线由光纤连接到探头，光线经样本反射后，由探头接收再经光纤传输到光谱仪，光谱数据经USB接口发送到计算机上。使用计算机软件对光谱数据进行处理，并给出计算结果。

（1）建立14种生物炭建模光谱库：对以上14种不同源质品种的生物炭材料，每种材料8个样本，利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，采样间隔1nm，共751个波长点，记录光强值，并作基线偏置校正，组建生物炭建模光谱库。典型生物炭发射光谱如图2所示。

基线偏置校正公式如下：

F(x)＝f(x)-minf(X)

式中，f(x)是波长x的光强值，X表示所有波长变量，F(x)表示波长x校正后的光强值，minf(X)表示所有波长点的光强值的最小值；

生物炭建模光谱库由所有样本的光谱记录组成。每一样本的光谱记录包括生物炭品种名称（1字段）、不同波长点校正后的光强值（751字段），归类代码（1字段）。

（2）对生物炭建模光谱库，利用Unscrambler X10.1软件主成分分析功能模块进行计算，按主成分贡献率大小，取贡献率最大的前3个主成分。

所述主成分分析的计算原理如下：

对于n个样本，每个样本用p个变量来描述研究对象，分别用X₁、X₂…X_p表示，X_p表示第p个波长点所有样本的光强值按顺序组成的向量，这p个变量构成的p维随机向量为X=(X₁，X₂…X_p)^t；

计算相关系数矩阵R，

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1p}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

其中R_ij（i，j＝1,2…p）为原始变量X_i与X_j的相关系数，且R_ij=R_ji，其计算公式为：

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}^2}}$

根据相关系数矩阵R解特征方程，

|λE-R|=0

求出特征值λ_i（i=1,2…p）；将λ_i按大小顺序排序，即

λ₁≥λ₂≥…≥λ_i≥0

分别求出对应于特征值λ_i的特征向量e_i(i=1,2,…,p)，要求||e_i||=1，即

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

其中e_ij表示向量e_i的第j个分量;

计算主成分Z_i的贡献率W_i，计算公式为

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

根据主成分贡献率大小，取所述的贡献率最大的前3个主成分Z₁，Z₂，Z₃；

根据大小顺序排前三的λ₁、λ₂、λ₃，按照下式计算三个主成分的主成分载荷

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}(i=\mathrm{1,2,3},j=\mathrm{1,2},\·\·\·,p)$

得到各主成分的载荷以后，按照下式计算各主成分的得分

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=l_{21}{X}_1+l_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

则Z₁、Z₂、Z₃表示n个样本的3个新变量指标；[Z₁ Z₂ Z₃]是一个n×3的新变量矩阵；

（3）根据前3个主成分，将生物炭同一品种不同样本的数据取均值，计算生物炭每一品种的三维中心坐标。

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m},(k=\mathrm{1,2,3})$

k表示第k种主成分，也就是第k维；

生物炭建模光谱库第i品种样本光谱第k维主成分平均值；

x_i，j,k:生物炭建模光谱库第i品种第j样本光谱第k维主成分值；

m:生物炭建模光谱库第i品种的样本数；本实施例中，m=8。

则生物碳第i品种的三维中心坐标即为

所得结果如表2。

表2

（4）待鉴别样本光谱扫描：利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，采样间隔1nm，共751个波长点，记录光强值，并作基线偏置校正，方法同步骤（1），每个品种2个样本，共28个样本，组建待鉴别样本光谱库。

（5）按步骤（2）同样方法对待鉴别样本光谱库作主成分分析，按主成分方差贡献率大小，提取贡献率最大的前3个主成分，计算各主成分的得分值，从而得到每一待鉴别样本的三维中心坐标；如表3所示。

表3

并计算每个待鉴别样本与每一品种的三维中心坐标之间的马氏距离，结果见表4。

对任一待鉴别样本x=(x₁,x₂,x₃)^t，，第i品种的中心坐标

马氏距离公式为：

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}$

式中，S^-1是S的逆矩阵，S是生物炭建模光谱库中第i品种所有样本的主成分得分矩阵的协方差矩阵,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

其中S_kl（k，l=1,2,3）为第i品种中所有样本第k维的主成分得分值向量Z_k与所有样本第l维的主成分得分值向量Z_l之间的协方差，且S_kl=S_lk，其计算公式为：

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

m:生物炭建模光谱库第i品种的样本数；

表4

（6）根据待鉴别样本到每一品种中心坐标之间的最小马氏距离，将待鉴别样本归类到相应的生物炭品种，如表4所示，从而完成生物炭溯源鉴别。根据本发明所述方法对上述28个待鉴别样本进行归类测试，其中27个样本归类正确，1个F类样本（生活垃圾生物炭）误判为L类（柠条生物炭），鉴别误差为3.6%。

Claims (2)

1.一种基于光谱技术的生物炭溯源鉴别方法，其特征在于该方法的步骤如下：

（1）建立生物炭品种光谱库：收集已知不同源质品种的生物炭材料，每种材料5-10个样本，利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，采样间隔1nm，共记录751个波长点的光强值，并作基线偏置校正，组建生物炭建模光谱库；

基线偏置校正公式如下：

F(x)f(x)-minf(X)

式中，f(x)是波长x的光强值，X表示所有波长变量，F(x)表示波长x校正后的光强值，minf(X)表示所有波长点的光强值的最小值；

生物炭品种光谱库由所有样本的光谱记录组成，每一样本的光谱记录包括生物炭品种名称、不同波长点校正后的光强值；

（2）对生物炭建模光谱库进行主成分分析，按主成分方差贡献率大小，取所述的贡献率最大的前3个主成分；

所述主成分分析的计算步骤如下：

对于n个样本，每个样本用p个变量来描述研究对象，分别用X₁、X₂…X_p表示，X_p表示第p个波长点所有样本的光强值按顺序组成的向量，这p个变量构成的p维随机向量为X=(X₁，X₂…X_p)^t；

计算相关系数矩阵R，

$R=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \·\·\·& R_{1p}\\ R_{21}& R_{22}& \·\·\·& R_{2p}\\ \·& \·& & \·\\ \·& \·& & \·\\ \·& \·& & \·\\ R_{p1}& R_{p2}& \·\·\·& R_{\mathrm{pp}}\end{array}\right]$

其中R_ij（i，j＝1,2…p）为原始变量X_i与X_j的相关系数，且R_ij=R_ji，其计算公式为：

$R_{\mathrm{ij}}=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}{\sqrt{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{ki}}-{\stackrel{\‾}{X}}_i)}^2\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(X_{\mathrm{kj}}-{\stackrel{\‾}{X}}_j)}^2}}$

n表示所有样本的样本数；

根据相关系数矩阵R解特征方程，

|λE-R|=0

求出特征值λ_i（i=1,2…p）；将λ_i按大小顺序排序，即

λ₁≥λ₂≥…≥λ_i≥0

分别求出对应于特征值λ_i的特征向量e_i(i=1,2,…,p)，要求||e_i||=1，即

$\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{e}_{\mathrm{ij}}^2=1$

其中e_ij表示向量e_i的第j个分量;

计算主成分Z_i的贡献率W_i，计算公式为

$W_i=\frac{{\mathrm{\λ}}_i}{{\mathrm{\Σ}}_{j=1}^p{\mathrm{\λ}}_i}$

根据主成分方差贡献率大小，取所述的贡献率最大的前3个主成分Z₁，Z₂，Z₃；

根据大小顺序排前三的λ₁、λ₂、λ₃，按照下式计算三个主成分的主成分载荷

$u_{\mathrm{ij}}=p(Z_i,X_j)=\sqrt{{\mathrm{\λ}}_i}{e}_{\mathrm{ij}}(i=\mathrm{1,2,3},j=\mathrm{1,2},\·\·\·,p)$

得到各主成分的载荷以后，按照下式计算各主成分的得分

$\left\{\begin{array}{c}{Z}_1=u_{11}{X}_1+u_{12}{X}_2+\·\·\·+u_{1p}{X}_p\\ Z_2=l_{21}{X}_1+l_{22}{X}_2+\·\·\·+u_{2p}{X}_p\\ Z_3=u_{31}{X}_1+u_{32}{X}_2+\·\·\·+u_{3p}{X}_p\end{array}\right.$

则Z₁、Z₂、Z₃表示n个样本的3个新变量指标；[Z₁ Z₂ Z₃]是一个n×3的新变量矩阵；

（3）根据前3个主成分Z₁，Z₂，Z₃，将生物炭同一品种不同样本的数据取均值，计算生物炭每一品种的三维中心坐标；

${\stackrel{\‾}{x}}_{i,k}=\frac{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{x}_{i,j,k}}{m},(k=\mathrm{1,2,3})$

k表示第k种主成分，也就是第k维；

生物炭建模光谱库第i品种样本光谱第k维主成分平均值；

x_i，j,k:生物炭建模光谱库第i品种第j样本光谱第k维主成分得分值；

m:生物炭建模光谱库第i品种的样本数；

则生物碳第i品种的三维中心坐标即为

（4）待鉴别样本光谱扫描：利用光谱仪采集每个样本的反射光谱数据，记录波长范围350-1100nm的光谱数据，记录与步骤（1）相同的751个波长点的光强值，并作基线偏置校正，方法同步骤（1），组建待鉴别样本光谱库；

（5）对待鉴别样本光谱库作主成分分析，按主成分方差贡献率大小，提取贡献率最大的前3个主成分，方法同步骤（2），计算各主成分的得分值，从而得到每一待鉴别样本的三维中心坐标；计算待鉴别样本与每一品种的三维中心坐标之间的马氏距离，对任一待鉴别样本x=(x₁,x₂,x₃)^t，第i品种的中心坐标马氏距离公式为：

$D(x,\stackrel{\‾}{x})=\sqrt{{(x-\stackrel{\‾}{x})}^t{S}^{-1}(x-\stackrel{\‾}{x})}$

式中，S^-1是S的逆矩阵，S是生物炭建模光谱库中第i品种所有样本的主成分得分矩阵的协方差矩阵,

$S=\left[\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& S_{23}\\ S_{31}& S_{32}& S_{33}\end{array}\right]$

其中S_kl（k，l=1,2,3）为第i品种中所有样本第k维的主成分得分值向量Z_k与所有样本第l维的主成分得分值向量Z_l之间的协方差，且S_kl=S_lk，其计算公式为：

$S_{\mathrm{kl}}=\frac{\underset{t=1}{\overset{m}{\mathrm{\Σ}}}(Z_{\mathrm{tk}}-{\stackrel{\‾}{Z}}_k)(Z_{\mathrm{tl}}-{\stackrel{\‾}{Z}}_l)}{m-1}$

m:生物炭建模光谱库第i品种的样本数；

（6）根据待鉴别样本到每一品种中心坐标之间的最小马氏距离，将待鉴别样本归类到相应的生物炭品种，从而完成生物炭溯源鉴别。

2.如权利要求1所述的方法，其特征在于：以下14种生物炭的光谱三维主成分中心坐标如下表所示：

Images