一种光谱归一化结合多元统计模型的玉石原产地鉴定方法与流程

本发明涉及一种对玉石进行原产地鉴定的方法，具体来说，其基本原理是利用激光诱导击穿光谱技术(LIBS)，采用光谱归一化进行预处理，再结合主成分分析(PCA)、支持向量机(SVM)这两种多元统计模型分析光谱数据，进行玉石的原产地鉴定。

背景技术：

玉石是一类珍贵的矿物，通常会被打造为各类首饰、装饰品和艺术品。高品质的玉石具有很高的经济和工艺价值。在中国，玉石文化已有上千年的历史，蕴含了丰富的内涵。而玉石的品质和价格，主要由玉石原料的产地决定。与鉴别玉石真假相比，玉石原产地的鉴定难度更大。通过传统的色泽、经验判断和常规光学检测的原产地鉴别方法正确率只有60％左右，且需要相关专家才能进行，受人为因素影响很大。目前虽然也有采用拉曼光谱技术、电感耦合等离子体光谱技术对不同产地的玉石样品进行比较，却没有形成较为系统的鉴定方法。

激光诱导击穿光谱技术(简称LIBS)是一种新兴的元素分析技术，具有无需样品预处理、对样品损伤小、分析速度快、能实现多元素测量等多种优点，在玉石的原产地鉴定上具有很大的潜力。但LIBS光谱数据通常具有数据量大、维度多、易受实验条件波动影响的特点，难以直接用于原产地的鉴别工作，故本发明采用光谱归一化进行预处理，消除实验条件波动的影响，并结合多元统计学模型，先进行采用主成分分析(PCA)，提取主成分，去掉冗余数据和噪声，降低数据维度；再对提取的主成分采用支持向量机(SVM)进行分类，确定玉石的原产地。

技术实现要素：

本发明的目的是针对传统玉石原产地鉴定方法正确率低、受人为因素影响大的缺陷，提出通过LIBS技术，结合光谱归一化、多元统计模型PCA和SVM来进行鉴定，将定性的人为经验分析提升为更加科学的定量化鉴别，从而提高鉴定结果的正确率。

本发明的技术方案是：

一种光谱归一化结合多元统计模型的玉石原产地鉴定方法，其特征是该方法包括了如下步骤：

1)使用一组原产地已知的玉石样品，作为定标样品进行建模，将来自同一原产地的样品归为同一类，不同产地的样品归为不同类；

2)使用激光诱导击穿光谱实验系统对定标样品进行检测，得到这组定标样品的光谱谱线，里面包含了每个定标样品各种元素的激光诱导击穿光谱的特征谱线以及这些特征谱线的强度；

3)对所有定标样品的激光诱导击穿光谱进行归一化处理：即选择一条强度较高的特征谱线作为标准特征谱线，对于每个定标样品的光谱，各特征谱线的强度同时除以标准特征谱线的强度，作为归一化后的强度保留下来，形成一个光谱归一化后的特征谱线强度矩阵X，

其中，n表示用于定标的玉石样品的数量，p表示特征谱线的数量，x_i1,x_i2,…,x_ip表示第i个玉石样品光谱归一化后各特征谱线的强度；

4)对矩阵X进行主成分分析，提取出主成分：将矩阵X进行对角化，即找到一个正交矩阵A，使得

其中，A^T与X^T分别表示矩阵A和X的转置，λ₁，λ₂，…，λ_p是对角线上的特征值，并满足对角线上的特征值从大到小排序，即λ₁≥λ₂≥…≥λ_p；

选择前m个特征值λ₁，λ₂，…，λ_m，使得这m个特征值的和大于等于特征值总和的95％，即

这前m个特征值对应的维度，就是矩阵X的前m个主成分，分别记为S₁，S₂，…，S_m，满足：

S₁＝XA₁，S₂＝XA₂，…，S_m＝XA_m (4)

其中，A₁，A₂，…，A_m分别为正交矩阵A的第1，第2，…，第m列元素；

将主成分分析得到的主成分矩阵记为S，就有

其中，n表示用于定标的玉石样品的数量，m表示选取的特征值的个数，也就是主成分的个数，s_i1,s_i2,…,s_im表示第i个用于定标的玉石样品的前m个主成分；

5)设这n个用于定标的玉石样品来自k种原产地，每种产地包含C_i个样品，即

C₁+C₂+…C_k＝n (6)

采用支持向量机方法建立模型进行定标，对这n个用于定标的玉石样品按原产地两两分类；

先将来自第一个原产地的样品视作一类，剩下k-1个原产地的样品视作另一类；对于这两类样品，采用支持向量机方法，即找到一个m维的向量ω和一个常数b，使得第一个产地的每个玉石样品i的主成分[s_i1 s_i2 … s_im]均满足

ω^T[s_i1 s_i2 … s_im]+b≥+1 (7)

且存在样品使得等号成立，即第一个产地的某个玉石样品i^*的主成分[s_i*1 s_i*2… s_i*m]满足

ω^T[s_i*1 s_i*2 … s_i*m]+b＝+1 (8)

另一类中的每个玉石样品j的主成分[s_j1 s_j2 … s_jm]均满足

ω^T[s_j1 s_j2 … s_jm]+b≤-1 (9)

且存在样品使得等号成立，即另一类中的某个玉石样品j^*的主成分[s_j*1 s_j*2 … s_j*m]满足

ω^T[s_j*1 s_j*2 … s_j*m]+b＝-1 (10)

这样两类样品就被一个线性平面ω^Ts+b＝0分隔开了，且间隔距离为其中‖ω‖表示向量ω的模值；

满足上述条件的向量ω和常数b不止一组，取其中‖ω‖最小，也就是使两类样品之间的间隔距离最大的那组(ω₁^*,b₁^*)，作为划分第一个产地和剩下k-1个产地的玉石样品的最佳方式；

划分出第一个产地的玉石样品后，再将第二个产地的样品视作一类，剩下k-2个产地的样品视作第二类，采用上述方式进行划分，记录对应的(ω₂^*,b₂^*)；以此类推，直至划分出所有产地的样品；

6)对于一组原产地未知的玉石样品，作为测试样品进行预测，具体做法如下：

首先使用激光诱导击穿光谱实验系统对测试样品进行检测，得到光谱谱线；并对光谱进行归一化，得到归一化后的特征谱线强度矩阵X′；提取出归一化后的特征谱线强度矩阵X′的前m个主成分，记为S′，有

其中q表示用于测试的玉石样品的数量，m表示选取的特征值的个数，也就是主成分的个数，s_i1′,s_i2′,…,s_im′表示第i个用于测试的玉石样品的各个主成分；

采用支持向量机方法对这q个用于测试的玉石样品进行两两分类，预测其原产地；先划分出第一个产地的样品；对于第i个玉石样品，若主成分[s_i1′ s_i2′ … s_im′]满足(ω₁^*)^T[s_i1′ s_i2′ … s_im′]+b₁^*≥0，将其视为第一个产地的样品；若满足(ω₁^*)^T[s_i1′ s_i2′ … s_im′]+b₁^*<0，将其视为其他产地的样品；再在其他产地的样品中划分出第二个产地的样品，即对于第i个玉石样品，

若主成分满足(ω₂^*)^T[s_i1′ s_i2′ … s_im′]+b₂^*≥0，将其视为第二个产地的样品；

若满足(ω₂^*)^T[s_i1′ s_i2′ … s_im′]+b₂^*<0，将其视为其他产地的样品；以此类推，所有测试样品的原产地都预测完毕；

7)将测试样品的预测原产地与真实原产地进行对比，校验原产地鉴定的正确性。。

本发明具有以下优点及突出性的技术效果：

LIBS技术对玉石样品的烧蚀质量在纳克级别，损伤非常小，鉴定过程几乎不会对玉石品相产生影响，可以实现无损鉴定；具体操作过程中每个样品只需采集一个光谱，一批样品从采谱到分析鉴定完毕只需几分钟的时间，可实现快速鉴定；采用LIBS谱线数据作为分类指标，将定性的人为经验分析提升为更加科学的定量化鉴别，显著提高了原产地鉴定的正确性。使用光谱归一化进行预处理，避免了实验条件波动对玉石原产地鉴定的影响；采用PCA对初始数据进行预处理，只保留重要的主成分，去除了不必要的维度，大大降低了模型计算时的时间和空间复杂度；对于主成分采用SVM建模，进行分类实现原产地鉴定，将SVM在优越的分类能力与PCA降维的特点相结合，可以得到极高的鉴定正确率。

附图说明

图1是本发明的流程示意图。

图2a为用第一和第二主成分为横纵坐标绘制的各产地样品的数据分布，图2b为前七个主成分的特征值占所有特征值总和的百分比。

图3是产地鉴定结果图。

具体实施方式

下面结合附图和实施例对本发明做进一步的说明。

本发明提供的一种光谱归一化结合多元统计模型的玉石原产地鉴定方法，其具体包括如下步骤：

1)使用一组原产地已知的玉石样品，作为定标样品进行建模，将来自同一原产地的样品归为同一类，不同产地的样品归为不同类；

2)使用激光诱导击穿光谱实验系统对定标样品进行检测，得到这组定标样品的光谱谱线，里面包含了每个定标样品各种元素的激光诱导击穿光谱的特征谱线以及这些特征谱线的强度，主要是Ca、Mg、Si等元素的各类原子谱线和离子谱线及其强度；

3)对所有定标样品的激光诱导击穿光谱进行归一化处理：即选择一条强度较高的特征谱线作为标准特征谱线，比如Ca元素在波长616.129nm处的原子谱线，谱线强度高且线型较好，适合作为标准特征谱线；对于每个定标样品的光谱，各特征谱线的强度同时除以标准特征谱线的强度，作为归一化后的强度保留下来，形成一个光谱归一化后的特征谱线强度矩阵X，

其中，n表示用于定标的玉石样品的数量，p表示特征谱线的数量，x_i1,x_i2,…,x_ip表示第i个玉石样品光谱归一化后各特征谱线的强度；

4)对矩阵X进行主成分分析，提取出主成分：将矩阵X进行对角化，即找到一个正交矩阵A，使得：

其中，A^T与X^T分别表示矩阵A和X的转置，λ₁，λ₂，…，λ_p是对角线上的特征值，并满足对角线上的特征值从大到小排序，即λ₁≥λ₂≥…≥λ_p；

选择前m个特征值λ₁，λ₂，…，λ_m，使得这m个特征值的和大于等于特征值总和的95％(舍去后5％数值较小的特征值，这些特征值对应的维度对分类作用很小，可以直接剔除；

舍去的百分比可以调整，通常5％-10％均可，既能保证提取的主成分信息充足，也可避免掺杂无用成分)，即

这前m个特征值对应的维度，就是矩阵X的前m个主成分，分别记为S₁，S₂，…，S_m，满足：

S₁＝XA₁，S₂＝XA₂，…，S_m＝XA_m (4)

其中，A₁，A₂，…，A_m分别为正交矩阵A的第1，第2，…，第m列元素；

将主成分分析得到的主成分矩阵记为S，就有

其中，n表示用于定标的玉石样品的数量，m表示选取的特征值的个数，也就是主成分的个数，s_i1,s_i2,…,s_im表示第i个用于定标的玉石样品的前m个主成分，这m个主成分，实际上是原光谱强度数据的线性组合，但比起原光谱进一步地提取了对于玉石样品进行原产地鉴定的有用信息，剔除了噪声等无用信息；

5)设这n个用于定标的玉石样品来自k种原产地，每种产地包含C_i个样品，即

C₁+C₂+…C_k＝n (6)

采用支持向量机方法建立模型进行定标，对这n个用于定标的玉石样品按原产地两两分类；

先将来自第一个原产地的样品(比如原产地为罗甸的玉石)视作一类，剩下k-1个原产地的样品(新疆、青海、俄罗斯、韩国等其他原产地的玉石)视作另一类；对于这两类样品，采用支持向量机方法，即找到一个m维的向量ω和一个常数b，使得第一个产地的每个玉石样品i的主成分[s_i1 s_i2 … s_im]均满足

ω^T[s_i1 s_i2 … s_im]+b≥+1 (7)

且存在样品使得等号成立，即第一个产地的某个玉石样品i^*的主成分[s_i*1 s_i*2… s_i*m]满足

ω^T[s_i*1 s_i*2 … s_i*m]+b＝+1 (8)

另一类中的每个玉石样品j的主成分[s_j1 s_j2 … s_jm]均满足

ω^T[s_j1 s_j2 … s_jm]+b≤-1 (9)

且存在样品使得等号成立，即另一类中的某个玉石样品j^*的主成分[s_j*1 s_j*2 … s_j*m]满足

ω^T[s_j*1 s_j*2 … s_j*m]+b＝-1 (10)

这样两类样品就被一个线性平面ω^Ts+b＝0分隔开了，两类样品分别位于这个线性平面的两侧，并且间隔距离为其中‖ω‖表示向量ω的模值。

满足上述条件的向量ω和常数b不止一组，取其中‖ω‖最小，也就是使两类样品之间的间隔距离最大、区分最明显的那组(ω₁^*,b₁^*)，作为划分第一个产地和剩下k-1个产地的玉石样品的最佳方式；

划分出第一个产地的玉石样品后，再将第二个产地的样品视作一类(比如在划分出罗甸的玉石后，接着划分新疆的玉石样品)，剩下k-2个产地的样品视作第二类，采用上述方式进行划分，记录对应的(ω₂^*,b₂^*)；以此类推，直至划分出所有产地的样品；

6)对于一组原产地未知的玉石样品，作为测试样品进行预测，具体做法如下：

首先使用激光诱导击穿光谱实验系统对测试样品进行检测，得到光谱谱线；并对光谱进行归一化，得到归一化后的特征谱线强度矩阵X′；提取出归一化后的特征谱线强度矩阵X′的前m个主成分，记为S′，有

其中q表示用于测试的玉石样品的数量，m表示选取的特征值的个数，也就是主成分的个数，s_i1′,s_i2′,…,s_im′表示第i个用于测试的玉石样品的各个主成分；

采用支持向量机方法对这q个用于测试的玉石样品进行两两分类，预测其原产地。先划分出第一个产地的样品。对于第i个玉石样品，若主成分[s_i1′ s_i2′ … s_im′]满足(ω₁^*)^T[s_i1′ s_i2′ … s_im′]+b₁^*≥0，将其视为第一个产地的样品；若满足(ω₁^*)^T[s_i1′ s_i2′ … s_im′]+b₁^*<0，将其视为其他产地的样品；再在其他产地的样品中划分出第二个产地的样品，即对于第i个玉石样品，

若主成分满足(ω₂^*)^T[s_i1′ s_i2′ … s_im′]+b₂^*≥0，将其视为第二个产地的样品；

若满足(ω₂^*)^T[s_i1′ s_i2′ … s_im′]+b₂^*<0，将其视为其他产地的样品；以此类推，直至所有测试样品的原产地都预测完毕；

7)将测试样品的预测原产地与真实原产地进行对比，校验原产地鉴定的正确性。

实施例：对来自5种原产地的和田玉样品进行原产地鉴定。

638个和田玉样品的原产地鉴定，其中114个产于罗甸，114个产于新疆，110个产于青海，150个产于俄罗斯，150个产于韩国。

主要步骤为：

1)先使用500个已知原产地(分别产于罗甸、新疆、青海、俄罗斯和韩国，每个产地100个样品)的和田玉样品作为定标样品，建立模型：使用激光诱导击穿光谱实验系统对每个样品进行检测，得到光谱，找到特征谱线，主要包括Mg在285.2nm、382.9nm、383.2nm等处的原子线，Ca在487.812nm、616.129nm、714.814nm等处的原子线，Si在288.157nm、265.977nm等处的原子线和离子线。扣除背景，积分计算各谱线区域的面积作为谱线强度。

2)选择Ca在616.129nm处的原子线作为标准谱线，进行光谱归一化。

3)采用PCA对所有定标样品的谱线强度数据进行预处理，得到前五个主成分的特征值分别占特征值总和的57.5％、25.9％、7.1％、4.0％、1.8％，其占比达到96.3％，超过95％，故提取每个样品前5个主成分的数据。

4)采用SVM方法建立模型，对500个定标样品按原产地进行两两划分。模型准确性达到100％，即建模过程中每个用于定标的玉石样品的原产地均可以被正确鉴定。

5)为了验证鉴定方法的正确性，使用事先未知原产地的138个和田玉样品(14个产于罗甸，14个产于新疆，10个产于青海，50个产于俄罗斯，50个产于韩国)作为测试样品，使用已建立的模型，进行原产地预测。预测结果的正确率也达到了100％，得到了很好地原产地鉴定结果。

图2a为用第一和第二主成分为横纵坐标绘制的各产地样品的数据分布，图2b为前七个主成分的特征值占所有特征值总和的百分比。

图3是产地鉴定结果图，对于5种原产地的玉石样品，定标和校验的鉴定正确率均达到100％。