CN104502329B - 一种libs成分分析中元素最优工作曲线的建立方法 - Google Patents

Abstract

一种LIBS成分分析中元素最优工作曲线的建立方法，属于材料成分分析领域。本方法的基本思想是在获得材料每个元素的多条工作曲线的基础上，按照曲线确定度值高和曲线各点浓度偏差加权和小的约束条件，优选出各元素的最佳工作曲线。本方法由两个约束条件和六个操作步骤组成。运用该方法，以铅黄铜标样GSB04‑2146‑2008为实施对象，得到了样本中各元素的最优工作曲线。

Description

一种LIBS成分分析中元素最优工作曲线的建立方法

技术领域

本发明属于材料成分分析技术领域，特别涉及一种LIBS成分分析中元素最优工作曲线的建立方法。

背景技术

激光诱导击穿光谱技术(Laser Introduced breakdown spectroscopy，简称LIBS)是一种新型光谱分析技术，其原理是：当一束高能脉冲激光聚焦到样本上时，激光会将焦点附近的微量样本烧蚀成为高温等离子体，随后等离子体迅速冷却并发射出反映元素特征的光谱，通过测量特征光谱的波长和强度就能得到元素的类型和含量。该技术具有实时检测多种元素、无需样本预处理、样本损伤小等特点(王海舟.冶金分析前沿[M].北京：科学出版社，2004:255-289)。

在LIBS成分分析中，常常使用标准样本来建立元素的工作曲线。现有文献(Loree,T.R,Radziemski,L.J.Laser-induced breakdown spectroscopy:a technique foratomic detection and molecular identification[J].Los Alamos Conference onOptics,1981,288,232-241.)、(Radziemski L.,Cremers D.,Benelli K.,Khoo C.,andHarris R.D.,Use of the vacuum ultraviolet spectral region for LIBS-basedMartian geology and exploration[J],Spectrochimica Acta B,2005,60,237-248.)、(Sallé,B.Cremers D.A.,Maurice S.,and Wiens R.C.,Laser-induced breakdownspectroscopy for space xploration applications:Influence of ambient pressureon the calibration curves prepared from soil and clay samples[J],Spectrochimica Acta B,2005,60,479-490.)一般讨论的是元素的某条曲线的建立和修正问题，对于元素的多条工作曲线优选，还没有系统化论述。

不同的工作曲线通常对于成分分析的精度存在不同的影响，所以如何在多条曲线间进行优选就成了一个很重要的问题。

发明内容

本发明提供一种LIBS成分分析中元素最优工作曲线的建立方法，实现了多条曲线间进行优选。

本发明建立曲线取舍的优化模型，按此模型从元素的多条曲线中选择最优的工作曲线。其内容如下：

(1)对于每个给定元素，工作曲线最优化目标函数为：

1)对于该元素的所有工作曲线，曲线拟合的确定度值R最高；

2)对于该元素的所有工作曲线，曲线上各样本浓度的相对偏差加权和sum_dif最小。

其中R和sum_dif的定义如下：

dif(i)＝|(f(i)-y(i))/y(i)| i＝1,2…n (1.2)

i：样本号；

n：样本总数；

R：工作曲线的确定度值；

f(i)：元素第i个样本的浓度的拟合值；

：元素n个样本的浓度的平均值；

y(i)：元素第i个样本的浓度值；

dif(i)：元素第i个样本的浓度的相对偏差；

W(i)：元素第i个样本的相对偏差的加权值，满足W(i)>＝0,求和；

sum_dif：元素的相对偏差加权和。

W(i)值的选择由分析要求来决定。如果对某个浓度点的精度要求高，就选择较大的权重值，反之，选择较小的权重值。

(2)各元素工作曲线最优化计算方法为(如图1所示)：

1)测量一组标准样本(P个)的实际光谱数据，分别存放在不同路径下的PracSpectraData文件中。每个文件中第1列是光谱波长，第2列到第N+1列是光谱强度。

2)选择第1个样本，进行谱线自动识别，获得各个元素的实际谱线波长InstrMatPracWaveLen4(i,j)。其中，i为元素号，j为谱线号。

其中自动识别方法使用专利CN201310598593(一种LIBS成分分析中元素谱线的自动识别方法)的方法：

(a)定义限制参数

a1)谱线波长匹配阈值D₁

D₁为待匹配实测谱线波长与前一个实测谱线波长的间隔的一半。

a2)谱线波长匹配阈值D₂

D₂为待匹配实测谱线波长与后一个实测谱线波长的间隔的一半。D₁和D₂能不遗漏地和少冗余地覆盖和待匹配实测谱线临近的理论谱线。

a3)光谱强度下限值ILL

ILL可综合光谱仪的光谱强度测量范围和研究者的经验制定。ILL不小于零，一般不小于光谱仪空白噪声均值的2～3倍。

a4)光谱强度上限值IUL

IUL也是综合光谱仪的光谱强度测量范围和研究者的经验来制定的。IUL不大于光谱仪强度测量上限。

(b)制定识别规则

b1)-D₁<＝λ_{元素i，理论谱线j}-λ_{实测谱线k}<＝D₂

其中λ_{元素i，理论谱线j}为元素i的第j条理论谱线的波长，λ_{实测谱线k}为第k条实测谱线的波长。这里的理论谱线即是已知的原子光谱数据库的谱线。D₁＝(λ_{实测谱线k}‐λ_{实测谱线k-1})/2，D₂＝(λ_{实测谱线k+1}‐λ_{实测谱线k})/2，这里实测谱线波长序列{λ_{实测谱线k}，k＝1,2…r}是按从小到大排列的。该规则用于划分实测谱线k为元素i的谱线，这时实测谱线k的波长λ_{实测谱线k}就是元素i的实测谱线波长，记为λ_元素i,实_测谱线k。需要指出的是，一方面，由于单个元素自身的理论谱线波长的间隔常常小于实测光谱波长间隔，所以可能存在一条实测谱线对应于某个元素的多条理论谱线的情形；另一方面，由于实测光谱波长间隔随波长的不同有微小浮动，所以还可能存在同一条理论谱线对应于多条不同实测谱线的情形；另外，由于属于不同元素的谱线波长的间隔也常常小于实测光谱波长间隔，所以可能存在一条实测谱线对应于多个元素所属的谱线的情形。对于这三方面，要分别使用规则b2)、b3)和b4)来剔除重复的谱线。

b2)min{|λ_{元素i，实测谱线k}-λ_{元素i，理论谱线j}|，j＝1,2…m}其中λ_{元素i，实测谱线k}是元素i的第k条实测谱线波长，{λ_{元素i，理论谱线j}，j＝1,2…m}是元素i的与λ_{元素i，实测谱线k}对应的m条理论谱线波长的集合。

该规则用于在λ_{元素i，实测谱线k}对应的所有理论谱线波长中，只选择离此实测谱线波长最近的理论谱线波长。

b3)min{|λ_{元素i，理论谱线j}-λ_{元素i，实测谱线k}|，k＝1,2…n}其中λ_{元素i，理论谱线j}是元素i的第j条理论谱线波长，{λ_{元素i，实测谱线k}，k＝1,2…n}是元素i的与λ_{元素i，理论谱线j}对应的n条实测谱线波长的集合。

该规则用于在λ_{元素i，理论谱线j}对应的所有实测谱线波长中，只选择离此理论谱线波长最近的实测谱线波长。

b4)min{|λ_{元素i，实测谱线j}-λ_{元素i，理论谱线j}|，i＝1,2…p}，这里p条λ_{元素i，实测谱线j}相等其中{λ_{元素i，理论谱线j}，i＝1,2…p}是p个元素的理论波长，这些理论波长对应的实测谱线波长λ_{元素i，实测谱线j}都相等。

该规则用于在同一实测谱线属于不同元素的情形下，选择和此实测谱线波长最接近的理论谱线波长所对应的元素为该实测谱线所属的元素。

b5)ILL≤I_{元素i,实测谱线j}≤IUL

其中，I_{元素i,实测谱线j}是元素i的实测谱线j的强度。该规则用于剔除受到仪器、外部环境或光谱信号自身的影响而导致的谱线光谱强度过大或过小的谱线。

(c)建立识别算法

c1)按光谱仪波长测量范围读取理论谱线波长数据并存储

从各个元素对应的理论波长数据文件TheoWavLen1_元素名中读取波长在光谱仪测量范围[SL,SH]内的数据，并按从小到大顺序存储在各元素的理论谱线波长数组TheoWavLen2_元素名中。这里SL和SH分别为光谱仪的波长测量下限值和上限值。

c2)读取样本的实测光谱数据并在平均后存储

从样本实测光谱数据文件PracSpectraData中读取波长和多次测量的光谱强度数据，计算各波长处光谱强度的平均值，将波长和平均强度存储在二维数组PracAvgSpectraData中，数组的第一列为波长，第二列为平均强度。

c3)对指定的每个元素i，按规则b1)划分出属于此元素的实测谱线和相应的理论谱线。

对每个元素i，按-D₁<＝λ_{元素i，理论谱线j}-λ_{实测谱线k}<＝D₂匹配规则，得到每个元素的实测谱线波长及相应的平均强度。

c4)对指定的每个元素i，按规则b2)，在每条实测谱线对应的所有理论谱线波长中，只保留与实测谱线最接近的理论谱线波长并存储。

对每条实测谱线，求它和与它对应的每条理论谱线的波长差，取波长差最小的理论谱线作为对应于实测谱线波长的元素波长。将得到的元素实测谱线波长数据、其对应的光谱强度数据和理论谱线波长数据，分别存储到对应的二维数组InstrMatTheoWaveLen1(i,j)、InstrMatstrength1(i,j)、InstrMatPracWaveLen1(i,j)。

c5)对指定的每个元素i,按规则b5)去除强度过小或过大的谱线并存储。

对于元素谱线的光谱强度，按照ILL≤I_{元素i,实测谱线}≤IUL规则，去除光谱强度值过低的谱线以减小背景光谱影响，增大信噪比；去除光谱强度值过高的谱线，以去除光谱强度饱和的谱线。存储元素i对应的理论谱线波长、实测谱线波长和实测谱线强度到对应的二维数组InstrMatTheoWaveLen2(i,j)、InstrMatstrength2(i,j)、InstrMatPracWaveLen2(i,j)。

c6)对指定的每个元素i，按规则b3)，在每条理论谱线对应的所有实测谱线波长中，只保留与理论谱线最接近的实测谱线波长并存储。

对每条理论谱线，求它和与它对应的每条实测谱线的波长差，取波长差最小的实测谱线和该理论谱线对应。将得到的元素实测谱线波长、其对应的光谱强度和理论谱线波长，分别存储到对应的二维数组InstrMatTheoWaveLen3(i,j)、InstrMatstrength3(i,j)、InstrMatPracWaveLen3(i,j)。

c7)按照规则b4)，在同一实测谱线属于不同元素的情形下，选择和此实测谱线波长最接近的理论谱线波长所对应的元素为该实测谱线所属的元素并存储。

存储各元素的理论谱线波长、实际谱线波长、实际谱线强度到二维数组InstrMatTheoWaveLen4(i,j)、InstrMatPracWaveLen4(i,j)、InstrMatstrength4(i,j)。

3)对于第1个到第P个样本的每一个样本，对光谱数据PracSpectraData中的N次光谱数据做平均，得到波长SpectraOrgWav(j,k)，及强度SpectraOrgStren(j,k)，k为样本号，j为谱线号。对于每个样本k，判断其光谱数据的波长SpectraOrgWav(i,j,k)和理论波长InstrMatTheoWaveLen4(i,j)，是否相等。如果相等，就将SpectraOrgStren(i,k)存储在SpectraStren(i,j,k)中，其中i为元素号，j为谱线号，k为样本号。

其中SpectraStren(i,j,k)在i＝1和i>1时分别存放的是基体元素的光谱强度和非基体元素的光谱强度。

4)对于每个样本k，计算每个非基体元素的各条谱线的相对强度值并存储。相对强度值得计算方法是：

ele_relative(m,n,k,s)＝SpectraStren(m,k,s)/SpectraStren(n,k,1)

其中m为非基体元素谱线号，n为基体元素谱线号，k为样本号，s>1为非基体元素号。

5)读取每个非基体元素的样本浓度，并拟合曲线，计算R值和sum_dif值。

读取每个非基体元素的所有样本的浓度，并存储在concentration(k,s)中，其中k为样本号，s>1为非基体元素号。对相对光谱强度ele_relative(m,n,k,s)和元素浓度concentration(k,s)进行线性拟合，得到工作曲线curve(m,n,s)。根据公式(1.1)和(1.3)分别计算出R值和sum_dif值，并分别存储在R(m,n,s)和sum_dif(m,n,s)中，其中m为非基体元素的有效谱线号，n为基体元素的有效谱线号，s>1为非基体元素号。

6)对于每个非基体元素，将其各条曲线的R值按从大到小的顺序排序，得到R值较高的前N条曲线。这里排序过的R值和相应的曲线分别记为order_R(m,n,s)和order_curve(m,n,s)。找到这些曲线对应的sum_dif值，存储在order_sum_dif(m,n,s)中。寻找这些sum_dif值的最小值，其相应的工作曲线即为该元素的最优工作曲线OptCurve(s)。

本发明的基本思想是在获得材料每个元素的多条工作曲线的基础上，按照曲线确定度值高和曲线各点浓度偏差加权和小的约束条件，优选出各元素的最佳工作曲线。本方法由两个约束条件和六个操作步骤组成。运用该方法，以铅黄铜标样GSB04-2146-2008为实施对象，得到了样本中各元素的最优工作曲线。

本发明的优点是能从元素的多条工作曲线中获得最优工作曲线，有提高LIBS成分分析的检测精度的潜力。

附图说明

图1为建立最优工作曲线流程图。

图2为LIBS实验装置。

图3为Pb元素最优工作曲线。

图4为Fe元素最优工作曲线。

图5为Ni元素最优工作曲线。

图6为Mn元素最优工作曲线。

图7为Al元素最优工作曲线。

图8为Si元素最优工作曲线。

图9为Sn元素最优工作曲线。

图10为Bi元素最优工作曲线。

图11为P元素最优工作曲线。

图12为S元素最优工作曲线。

具体实施方式

1.实施本方法的实验准备工作

LIBS实验系统如图2所示(吴少波,叶连慧,孙彦广,等.国产Nd:Yag美容用激光器用于LIBS成分分析的可行性研究[J].仪器仪表学报，2013，34(6):181-186)。使用MK-101Nd:Yag脉冲激光器，激光波1024nm、能量106mJ、脉宽14.2ns。使用Avaspace-2048-4-USB2光谱仪，波长范围200nm到850nm、实测光谱波长间隔0.02～0.13nm。信号发生器采用DG535脉冲延时发生器，它能通过激光器触发电路触发激光器以指定频率重复发出激光脉冲。所用计算机为本单位云计算平台中的一台虚拟机：Windows Server2003操作系统，3.75G内存，8核CPU@1800MHz/核。实验采用铅黄铜标样GSB04-2146-2008，其样本成分如表1所示。

表1铅黄铜标样GSB04-2146-2008中各成分浓度

	Cu	Pb	Fe	Ni	Mn	Al	Si	Sn	Bi	Sb	P	S	Zn
														1	59.98	0.414	0.205	1.474	0.119	0.0749	0.0049	0.0099	0.011	0.0056	0.004	0.00077	37.69693
2	57.77	0.76	0.0426	0.795	0.184	0.0116	0.0093	0.0979	0.005	0.0021	0.0075	0.0015	40.3135
														3	57.09	1.421	0.0104	0.347	--*	0.177	0.052	0.192	0.003	0.0116	0.0141	0.00099	40.68091
4	58.64	1.81	0.703	0.104	0.0188	0.452	0.007	0.293	0.0049	0.0102	0.0111	0.0006	37.9454
														5	58.76	2.405	0.107	0.0286	0.0517	0.761	0.209	0.689	0.001	0.0048	0.1	0.001	36.8819
6	59.6	1.393	0.613	0.386	0.0345	0.262	0.0425	0.206	0.003	0.0291	0.0274	0.00094	37.40256

*说明：标准样本没有给出该元素的含量，在分析时设此含量为零。

启动试验装置自动重复进行下述过程250次：信号发生器高能脉冲激光以1Hz重复频率通过聚焦镜汇聚到样本表面，烧蚀表面产生高温等离子体。同时，光电转换电路接收样本表面处激光的漫反射光并转换为电信号，此信号经过放大整形处理后形成TTL信号来触发光谱仪。光谱仪按照1.28us延时时间和1.05ms积分时间通过光纤采集到光谱数据存储在计算机中。

按照本发明所述的建立最优工作曲线的方法的各步骤来操作：

2.根据各元素工作曲线最优化计算方法1)和2)，进行谱线的读取和自动识别操作，获得各个元素的实际谱线波长InstrMatPracWaveLen4(i,j)。其中，i为元素号，j为谱线号。其中设置识别参数，考虑到仪器的空白噪声约为200，ILL＝1000使得强度的信噪比不小于5：1，设置IUL＝60000，使得强度远离仪器的饱和强度65535。

3.根据各元素工作曲线最优化计算方法3)读取每个非基体元素的所有有效谱线的每个样本的光谱强度值，SpectraStren(i,j,k)中，其中i为元素号(i＝1～12)，j为谱线号(j＝1,2,...)，k为样本号(k＝1～6)。

其中SpectraStren(i,j,k)在i＝1和i>1时分别存放的是基体元素的光谱强度和非基体元素的光谱强度。

4.根据各元素工作曲线最优化计算方法4)计算每个非基体元素的各条谱线的相对强度值并存储。并存储在ele_relative(m,n,k,s)中。其中m为非基体元素谱线号，n为基体元素谱线号，k为样本号，s>1为非基体元素号。

5.根据各元素工作曲线最优化计算方法5)读取每个非基体元素的样本浓度，并拟合曲线，计算R值和sum_dif值。

6.根据各元素工作曲线最优化计算方法6)，令W(i)＝1/Num，Num为样本数6。根据公式(1.2)计算得到每个非基体元素的所有有效谱线的每个样本的浓度相对偏差。根据公式(1.3)计算出曲线order_curve(m,n,s)的浓度相对偏差总和，并存储在sum_dif(m,n,s)中。对于每个非基体元素，将其各条曲线的R值按从大到小的顺序排序，得到R值较高的前5条曲线。这里排序过的R值和相应的曲线分别记为order_R(m,n,s)和order_curve(m,n,s)。找到这些曲线对应的sum_dif值，存储在order_sum_dif(m,n,s)中。寻找这些sum_dif值的最小值，其相应的工作曲线即为该元素的最优工作曲线OptCurve(s)。其中m为非基体元素序号，n为基体元素谱线号。根据各元素工作曲线最优化计算方法6)，分别得到以Cu为基体元素的Pb、Fe、Ni、Mn、Al、Si、Sn、Bi、Sb、P、S等11个非基体元素的最优工作曲线。

7.Pb元素的最大的5个R值分别为：0.999824、0.999104、0.998510、0.998504、0.998444，他们对应的sum_dif分别为0.012350、0.018478、0.024907、0.027192、0.017388。故选择R值为0.999824，sum_dif为0.01235的工作曲线做为最优工作曲线，如图3所示。

8.Fe元素的前5个最大R值分别为：0.999993、0.999958、0.999951、0.999948、0.999926，其对应的sum_dif分别为：0.031332、0.064399、0.052531、0.084727、0.107730，故选择R值为0.999993、sum_dif为0.031332的工作曲线作为最优工作曲线，如图4所示。

9.Ni元素的前5个最大R值分别为：0.999929、0.999916、0.999897、0.999887、0.999874，其对应的sum_dif分别为：0.088712、0.051129、0.063886、0.058346、0.045780，故选择R值为0.999874、sum_dif为0.045780的工作曲线作为最优工作曲线，如图5所示。

10.Mn元素的前5个最大R值分别为：0.999943、0.999927、0.999821、0.999808、0.999759，其对应的sum_dif分别为：0.012804、0.021452、0.020505、0.008376、0.02835，故选择R值为0.999808、sum_dif为0.008376的工作曲线作为最优工作曲线，如图6所示。

11.Al元素的前5个最大R值分别为：0.999831、0.999573、0.999448、0.999374、0.999354，其对应的sum_dif分别为：0.022892、0.1696、0.107312、0.079634、0.258836，故选择R值为0.999831、sum_dif为0.022892的工作曲线作为最优工作曲线，如图7所示。

12.Si元素的前5个最大R值分别为：0.998747、0.998114、0.997483、0.997461、0.997111，其对应的sum_dif分别为：0.309462、0.482691、0.511293、0.248241、0.596587，故选择R值为0.997461、sum_dif为0.248241的工作曲线作为最优工作曲线，如图8所示。

13.Sn元素的前5个最大R值分别为：0.999817、0.999707、0.999670、0.999654、0.999645，其对应的sum_dif分别为：0.093692、0.175365、0.232881、0.161750、0.067409，故选择R值为0.999645、sum_dif为0.067409的工作曲线作为最优工作曲线，如图9所示。

14.Bi元素的前5个最大R值分别为：0.991182、0.987128、0.98683、0.981936、0.981269，其对应的sum_dif分别为：0.116655、0.231407、0.199614、0.130868、0.194037，故选择R值为0.991182、sum_dif为0.116655的工作曲线作为最优工作曲线，如图10所示。

15.Sb元素的前5个最大R值分别为：0.974099、0.968224、0.967019、0.955336、0.950818，其对应的sum_dif分别为：0.188088、0.249178、0.392775、0.376703、0.612459，由于所有工作曲线的R值均小于0.98，故认为没有找到Sb元素的最优工作曲线。

16.P元素的前5个最大R值分别为：0.991914、0.981664、0.975140、0.971964、0.957121，其对应的sum_dif分别为：0.606903、0.585476、0.925426、1.099998、0.868065，故选择R值为0.991914、sum_dif为0.606903的工作曲线作为最优工作曲线，如图11所示。

17.S元素的前5个最大R值分别为：0.993783、0.992992、0.991031、0.984171、0.974133，其对应的sum_dif分别为：0.027087、0.030962、0.036993、0.040732、0.073668，故选择R值为0.993783、sum_dif为0.027087的工作曲线作为最优工作曲线，如图12所示。

Claims (2)

1.一种LIBS成分分析中元素最优工作曲线的建立方法，基于内标法建立元素的工作曲线并优化；其特征在于，包括如下工艺步骤：

(1)对于每个给定元素，工作曲线最优化目标函数为：

1)对于该元素的所有工作曲线，曲线拟合的确定度值R最高；

2)对于该元素的所有工作曲线，曲线上各样本浓度的相对偏差加权和sum_dif最小；

其中R和sum_dif的定义如下：

$R=\sqrt{{\mathrm{\&Sigma;}}_1^n|f\left(i\right)-\stackrel{\&OverBar;}{y}{|}^2/{\mathrm{\&Sigma;}}_1^n|y\left(i\right)-\stackrel{\&OverBar;}{y}{|}^2}---\left(1.1\right)$

dif(i)＝|(f(i)-y(i))/y(i)|i＝1,2…n (1.2)

$sum\_dif={\mathrm{\&Sigma;}}_1^{n}W\left(i\right)*dif\left(i\right)---\left(1.3\right)$

i：样本号，

n：样本总数，

R：工作曲线的确定度值，

f(i)：元素第i个样本的浓度的拟合值，

元素n个样本的浓度的平均值，

y(i)：元素第i个样本的浓度值，

dif(i)：元素第i个样本的浓度的相对偏差，

W(i)：元素第i个样本的相对偏差的加权值，满足W(i)>＝0,求和

sum_dif：元素的相对偏差加权和，

(2)各元素工作曲线最优化计算方法为：

1)测量一组P个标准样本的实际光谱数据，分别存放在不同路径下的PracSpectraData文件中；每个文件中第1列是光谱波长，第2列到第N+1列是光谱强度；

2)选择第1个样本，进行谱线自动识别，获得各个元素的实际谱线波长InstrMatPracWaveLen4(i,j)；其中，i为元素号，j为谱线号；

3)对于第1个到第P个样本的每一个样本，对光谱数据PracSpectraData中的N次光谱数据做平均，得到波长SpectraOrgWav(j,k)，及强度SpectraOrgStren(j,k)，k为样本号，j为谱线号；对于每个样本k，判断其光谱数据的波长SpectraOrgWav(j,k)和理论波长InstrMatTheoWaveLen4(i,j)是否相等；如果相等，就将SpectraOrgStren(i,k)存储在SpectraStren(i,j,k)中，其中i为元素号，j为谱线号，k为样本号；

其中SpectraStren(i,j,k)在i＝1和i>1时分别存放的是基体元素的光谱强度和非基体元素的光谱强度；

4)对于每个样本k，计算每个非基体元素的各条谱线的相对强度值并存储；相对强度值的计算方法是：

ele_relative(m,n,k,s)＝SpectraStren(m,k,s)/SpectraStren(n,k,1)

其中m为非基体元素谱线号，n为基体元素谱线号，k为样本号，s>1为非基体元素号；

5)读取每个非基体元素的样本浓度，并拟合曲线，计算R值和sum_dif值；

读取每个非基体元素的所有样本的浓度，并存储在concentration(k,s)中，其中k为样本号，s>1为非基体元素号；对相对光谱强度ele_relative(m,n,k,s)和元素浓度concentration(k,s)进行线性拟合，得到工作曲线curve(m,n,s)；根据公式(1.1)和(1.3)分别计算出R值和sum_dif值，并分别存储在R(m,n,s)和sum_dif(m,n,s)中，其中m为非基体元素谱线号，n为基体元素谱线号，s>1为非基体元素号；

6)对于每个非基体元素，将其各条曲线的R值按从大到小的顺序排序，得到R值较高的前N条曲线；这里排序过的R值和相应的曲线分别记为order_R(m,n,s)和order_curve(m,n,s)；找到这些曲线对应的sum_dif值，存储在order_sum_dif(m,n,s)中；寻找这些sum_dif值的最小值，其相应的工作曲线即为该元素的最优工作曲线OptCurve(s)。

2.根据权利要求1所述的方法，其特征在于，步骤(2)中进行谱线自动识别的步骤如下：

(a)定义限制参数

a1)谱线波长匹配阈值D₁

D₁为待匹配实测谱线波长与前一个实测谱线波长的间隔的一半；

a2)谱线波长匹配阈值D₂

D₂为待匹配实测谱线波长与后一个实测谱线波长的间隔的一半；D₁和D₂能不遗漏地和少冗余地覆盖和待匹配实测谱线临近的理论谱线；

a3)光谱强度下限值ILL

ILL可综合光谱仪的光谱强度测量范围和研究者的经验制定；ILL不小于零，一般不小于光谱仪空白噪声均值的2～3倍；

a4)光谱强度上限值IUL

IUL也是综合光谱仪的光谱强度测量范围和研究者的经验来制定的；IUL不大于光谱仪强度测量上限；

(b)制定识别规则

b1)-D₁<＝λ_{元素i，理论谱线j}-λ_{实测谱线k}<＝D₂

其中λ_{元素i，理论谱线j}为元素i的第j条理论谱线的波长，λ_{实测谱线k}为第k条实测谱线的波长；这里的理论谱线即是已知的原子光谱数据库的谱线；D₁＝(λ_{实测谱线k}‐λ_{实测谱线k-1})/2，D₂＝(λ_{实测谱线k+1}‐λ_{实测谱线k})/2，这里实测谱线波长序列{λ_{实测谱线k}，k＝1,2…r}是按从小到大排列的；该规则用于划分实测谱线k为元素i的谱线，这时实测谱线k的波长λ_{实测谱线k}就是元素i的实测谱线波长，记为λ_{元素i,实测谱线k}；一方面，由于单个元素自身的理论谱线波长的间隔常常小于实测光谱波长间隔，所以可能存在一条实测谱线对应于某个元素的多条理论谱线的情形；另一方面，由于实测光谱波长间隔随波长的不同有微小浮动，所以还可能存在同一条理论谱线对应于多条不同实测谱线的情形；另外，由于属于不同元素的谱线波长的间隔也常常小于实测光谱波长间隔，所以可能存在一条实测谱线对应于多个元素所属的谱线的情形；对于这三方面，要分别使用规则b2)、b3)和b4)来剔除重复的谱线；

b2)min{|λ_{元素i，实测谱线k}-λ_{元素i，理论谱线j}|，j＝1,2…m}其中λ_{元素i，实测谱线k}是元素i的第k条实测谱线波长，{λ_{元素i，理论谱线j}，j＝1,2…m}是元素i的与λ_{元素i，实测谱线k}对应的m条理论谱线波长的集合；

该规则用于在λ_{元素i，实测谱线k}对应的所有理论谱线波长中，只选择离此实测谱线波长最近的理论谱线波长；

b3)min{|λ_{元素i，理论谱线j}-λ_{元素i，实测谱线k}|，k＝1,2…n}其中λ_{元素i，理论谱线j}是元素i的第j条理论谱线波长，{λ_{元素i，实测谱线k}，k＝1,2…n}是元素i的与λ_{元素i，理论谱线j}对应的n条实测谱线波长的集合；

该规则用于在λ_{元素i，理论谱线j}对应的所有实测谱线波长中，只选择离此理论谱线波长最近的实测谱线波长；

b4)min{|λ_{元素i，实测谱线j}-λ_{元素i，理论谱线j}|，i＝1,2…p}，这里p条λ_{元素i，实测谱线j}相等，其中{λ_{元素i，理论谱线j}，i＝1,2…p}是p个元素的理论波长，这些理论波长对应的实测谱线波长λ_{元素i，实测谱线j}都相等；

该规则用于在同一实测谱线属于不同元素的情形下，选择和此实测谱线波长最接近的理论谱线波长所对应的元素为该实测谱线所属的元素；

b5)ILL≤I_{元素i,实测谱线j}≤IUL

其中，I_{元素i,实测谱线j}是元素i的实测谱线j的强度；该规则用于剔除受到仪器、外部环境或光谱信号自身的影响而导致的谱线光谱强度过大或过小的谱线；

(c)建立识别算法

c1)按光谱仪波长测量范围读取理论谱线波长数据并存储

从各个元素对应的理论波长数据文件TheoWavLen1_元素名中读取波长在光谱仪测量范围[SL,SH]内的数据，并按从小到大顺序存储在各元素的理论谱线波长数组TheoWavLen2_元素名中；这里SL和SH分别为光谱仪的波长测量下限值和上限值；

c2)读取样本的实测光谱数据并在平均后存储

从样本实测光谱数据文件PracSpectraData中读取波长和多次测量的光谱强度数据，计算各波长处光谱强度的平均值，将波长和平均强度存储在二维数组PracAvgSpectraData中，数组的第一列为波长，第二列为平均强度；

c3)对指定的每个元素i，按规则b1)划分出属于此元素的实测谱线和相应的理论谱线；

对每个元素i，按D₁<＝λ_{元素i，理论谱线j}-λ_{实测谱线k}<＝D₂匹配规则，得到每个元素的实测谱线波长及相应的平均强度；

c4)对指定的每个元素i，按规则b2)，在每条实测谱线对应的所有理论谱线波长中，只保留与实测谱线最接近的理论谱线波长并存储；

对每条实测谱线，求它和与它对应的每条理论谱线的波长差，取波长差最小的理论谱线作为对应于实测谱线波长的元素波长；将得到的元素实测谱线波长数据、其对应的光谱强度数据和理论谱线波长数据，分别存储到对应的二维数组InstrMatTheoWaveLen1(i,j)、InstrMatstrength1(i,j)、InstrMatPracWaveLen1(i,j)；

c5)对指定的每个元素i,按规则b5)去除强度过小或过大的谱线并存储；

对于元素谱线的光谱强度，按照ILL≤I_{元素i,实测谱线}≤IUL规则，去除光谱强度值过低的谱线以减小背景光谱影响，增大信噪比；去除光谱强度值过高的谱线，以去除光谱强度饱和的谱线；存储元素i对应的理论谱线波长、实测谱线波长和实测谱线强度到对应的二维数组InstrMatTheoWaveLen2(i,j)、InstrMatstrength2(i,j)、InstrMatPracWaveLen2(i,j)；

c6)对指定的每个元素i，按规则b3)，在每条理论谱线对应的所有实测谱线波长中，只保留与理论谱线最接近的实测谱线波长并存储；

对每条理论谱线，求它和与它对应的每条实测谱线的波长差，取波长差最小的实测谱线和该理论谱线对应；将得到的元素实测谱线波长、其对应的光谱强度和理论谱线波长，分别存储到对应的二维数组InstrMatTheoWaveLen3(i,j)、InstrMatstrength3(i,j)、InstrMatPracWaveLen3(i,j)；

c7)按照规则b4)，在同一实测谱线属于不同元素的情形下，选择和此实测谱线波长最接近的理论谱线波长所对应的元素为该实测谱线所属的元素并存储；

存储各元素的理论谱线波长、实际谱线波长、实际谱线强度到二维数组InstrMatTheoWaveLen4(i,j)、InstrMatPracWaveLen4(i,j)、InstrMatstrength4(i,j)。