CN105510427A - 一种用于多种元素同位素双稀释剂技术的数值处理方法 - Google Patents

Abstract

本发明公开了一种用于多种元素同位素双稀释剂技术中获取最优配比值的数据处理方法。该方法包括以下步骤：根据质谱仪器参数建立误差模型；基于误差模型，利用蒙特-卡罗法模拟得到质谱分析的样品测量值；利用鸟巢迭代算法校正测量值，得到对应的自然分馏因子值及其误差；以自然分馏因子误差和双稀释剂在混合样品中的比例为对象，作二维误差曲线图，最小误差点所对应的值即为最优配比参数值。同时，此算法也可获得自然样品的同位素比值。本发明为同位素双稀释剂的选择、双稀释剂间以及双稀释剂与样品的最佳混合有明显指导作用，能够显著降低该技术的分析成本，提高分析效率。其在非传统稳定同位素地球化学学科的方法论和认识论方面具有重要价值。

Description

一种用于多种元素同位素双稀释剂技术的数值处理方法

技术领域

本发明属于同位素地球化学和分析化学技术领域，尤其涉及一种用于多种元素同位素双稀释剂技术中获取最优配比值的数值优化方法。

背景技术

近年来，非传统稳定同位素地球化学的研究已经成为地球科学、环境科学等领域的重要前沿学科之一。同位素组成的准确分析是开展非传统稳定同位素理论和应用研究的重要技术基础。现有的主流分析技术包括标准-样品交叉法(Standard-SampleBracketing)、元素外标法(ElementDoping)和同位素双稀释剂法(IsotopeDoubleSpike)。由于前两者方法对质谱仪器运行时的稳定性、外标元素加入与所测同位素离子化时行为的一致性要求较高，且严重制约着后续同位素分析过程中的质量歧视校正，同位素双稀释剂分析技术越来越受到诸多研究者的重视。同位素双稀释剂法属于内标校正，能够有效校正元素同位素在化学分离、纯化和质谱仪测定过程中出现的分馏，实现同位素比值分析的长期重现性，显著提高同位素比值分析的准确度和精度，为之后的同位素理论和应用研究奠定基础。

同位素双稀释剂法相比较于标准-样品交叉法和元素外标法具有较为明显的优势，可以获得自然样品的高精度同位素比值和准确浓度值；能够校正样品分离纯化和质谱测试过程中产生的同位素分馏，简化样品分离纯化步骤和缓冲一定的干扰；特别是可用于低浓度自然样品的同位素比值测定和降低所测试元素同位素的回收率等。但其缺点是使用人员如果难以掌握两单稀释剂之间以及稀释剂与样品之间的最佳配比参数值，就难以实现同位素比值的高精度分析。相关参数的最佳配比可以通过实验实测获得，但其时间、经费和昂贵稀释剂的消耗较大，分析成本居高不下。本发明公开的同位素双稀释剂技术的数值优化方法，可以通过数值模拟获得上述参数，能够显著地提高分析效率和降低分析成本。

发明内容

为了解决上述现有技术中存在的问题，本发明提出一种用于多种元素同位素双稀释剂技术中获取最优配比值的数值处理方法，该方法基于蒙特-卡罗和鸟巢迭代算法，通过相关的数值模拟获得所需要的关键参数。

本发明提出的一种用于多种元素同位素双稀释剂中获取最优配比值的数值处理方法，包括以下步骤：

步骤1，根据质谱仪器参数建立误差模型；

步骤2，基于所述步骤1获得的误差模型，利用蒙特-卡罗法模拟得到样品在质谱分析过程中的测量值；

步骤3，利用鸟巢迭代算法校正所述步骤2得到的测量值，并得到对应的自然分馏因子值；

步骤4，计算自然分馏因子的误差；

步骤5，以自然分馏因子误差和双稀释剂在混合样品中的比例为对象，作二维误差曲线图，所述二维误差曲线图中最小误差点所对应的相关配比参数值即为最优配比参数值。

本发明方法不仅简单、快捷，极大地避免同位素稀释剂的浪费和显著降低同位素分析成本，而且也能够给出同位素双稀释剂的最佳选择、稀释剂之间以及稀释剂与样品间的最佳配比等参数，同时也能够利用此算法进行同位素比值的提取，极大推进同位素双稀释剂技术在我国相应科研领域中的应用。另外，基于本发明的数值优化方法是尚未公开的原创性研究成果，国内外尚无本发明内容的相关报道。

附图说明

图1是本发明一种用于多种元素同位素双稀释剂技术中获取最优配比值的数值处理方法的流程图。

图2是蒙特-卡罗法示意图。

图3是⁵⁰Cr-⁵⁴Cr双稀释剂E_α-p二维曲线图，E_α-p二维曲线图即为自然分馏因子的误差与双稀释剂在混合样品中比例的二维曲线图，(⁵⁰Cr-⁵⁴Cr)_稀释剂为两单稀释剂之间的比值。

图4是⁹⁷Mo-¹⁰⁰Mo双稀释剂E_α-p二维曲线图，(⁹⁷Mo-¹⁰⁰Mo)_稀释剂为两单稀释剂之间的比值。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明提出了一种用于多种元素同位素双稀释剂技术中获取最优配比值的数据处理方法，适用于周期表中能够使用同位素双稀释剂技术分析的所有元素同位素(同位素数≥4)，其目的主要为了获取同位素双稀释剂之间和双稀释剂与样品间的最佳混合比例。

图1是本发明一种用于多种元素同位素双稀释剂技术中获取最佳配比值的数值处理方法的流程图，下面基于蒙特-卡罗和鸟巢迭代算法，以⁵⁰Cr-⁵⁴Cr双稀释剂、⁹⁷Mo-¹⁰⁰Mo双稀释剂的相关配比参数优化过程为例，对于本发明方法进行说明，如图1所示，所述方法包括以下步骤：

步骤1，根据质谱仪器参数建立误差模型；

所述步骤1进一步包括以下步骤：

步骤1.1，获取质谱的仪器参数；

所述仪器参数至少包括但不限于仪器的频宽、仪器的阻抗、仪器运行时的温度，仪器的电阻器噪声(Johnsonnoise)和散粒效应(Shotnoise)等数值，其中，所述电阻器噪声和散粒效应的数值由下式计算得到：

a＝4kbTR；d＝ebR，

其中，a为电阻器噪声，d为散粒效应，k为Boltzmann常数，k＝1.3806504×10^-23JK^-1，b为频宽，即每秒可传送的离子数，R为电阻，T为质谱运行时的绝对温度，e为电子电荷，e＝1.602176487×10^-19C(库伦)，b为传送周期t的倒数。在⁵⁰Cr-⁵⁴Cr双稀释剂、⁹⁷Mo-¹⁰⁰Mo双稀释剂以及其它元素同位素双稀释剂相关配比参数的优化过程中，预设质谱仪器的标准配置参数R＝10¹¹Ω，T＝300K，t＝8s。

步骤1.2，根据所述步骤1.1得到的仪器参数建立误差模型。

所述误差模型表示为：

其中，i表示同位素序号，变化范围为1～4，比如在Cr同位素中1、2、3、4分别表示⁵⁰Cr、⁵²Cr、⁵³Cr、⁵⁴Cr，σ_{Johnson-noise}为电阻器噪声所产生的绝对误差，其与离子束强度无关，即σ_Shot-noise为散粒效应产生的绝对误差，其与离子束强度的平方根呈现正比例关系，即μ＝10V；σ_干扰主要为同质异位素和同质多原子分子的干扰，标准条件下设置为0。

由上述误差模型可获得同位素质谱分析过程中所产生的任意误差值。

步骤2，基于所述步骤1的误差模型，利用蒙特-卡罗法模拟得到样品在质谱分析过程中的模拟信号值，并计算出各个信号值间的比值，比如在Cr同位素的模拟过程中，通过该步骤可获得⁵⁰Cr/⁵²Cr，⁵³Cr/⁵²Cr，⁵²Cr/⁵²Cr的模拟信号比值，也即为测量值；

所述步骤2进一步包括以下步骤：

步骤2.1，获取双稀释剂的同位素组成；

该步骤中：若令λ表示某单稀释剂在双稀释剂中的比例，如下式所示，以⁵⁰Cr-⁵⁴Cr双稀释剂为例，λ值的变化能够产生一系列的双稀释剂离子束I：

I_s＝λ×I_50Spk+(1-λ)×I_54Spk，

其中，I_s为双稀释剂中Cr同位素在质谱分析过程中的离子束强度集，I_50Spk、I_54Spk分别为⁵⁰Cr、⁵⁴Cr单稀释剂中Cr同位素离子束强度集，这里的I值也就是质谱分析过程中的理想信号值，即为质谱分析时无电阻器噪声，散粒效应和质量歧视效应影响的模拟信号值。

其中，I_50Spk、I_54Spk的数学表达式为：

I_50Spk＝(I₅₀,I₅₂,I₅₃,I₅₄)^T,I_54Spk＝(I₅₀,I₅₂,I₅₃,I₅₄)^T，

其中，I₅₀、I₅₄、I₅₀、I₅₄分别为单稀释剂中各同位素的理想信号值。

根据双稀释剂中同位素的理想信号的比值可获得双稀释剂的同位素组成：

(ⁱCr/^kCr)_s＝I_s-i/I_s-k，

其中，k＝52,i＝50、53、54，I_S-i、I_S-k为I_S中相对应同位素的理想信号值。

步骤2.2，根据混合样品中同位素的理想信号的比值，获取自然样品与双稀释剂的混合样品的同位素组成；

其中，若令p表示双稀释剂在混合样品中的比例，则根据下式，p值的变化可获得一系列的混合样品离子束强度值：

I_M＝p×I_s+(1-p)×I_n，

其中，I_M为混合样品中的同位素在质谱分析过程中的理想信号集，I_s为双稀释剂中的同位素在质谱分析中的理想信号集，I_n为自然样品中同位素在质谱分析中的理想信号集，其可用标准溶液中的同位素组成来表示。

步骤2.3，获取所述混合样品的模拟测量值，以Cr为例，即为⁵⁰Cr/⁵²Cr，⁵³Cr/⁵²Cr，⁵²Cr/⁵²Cr的模拟测量值；

预设自然分馏因子(自然样品相对于同位素标准的分馏因子)为α，质谱仪器的质量歧视因子为β，在⁵⁰Cr-⁵⁴Cr双稀释剂、⁹⁷Mo-¹⁰⁰Mo双稀释剂的数值优化过程中α、β取值为0，则由下面指数律公式可以得到混合样品在质谱分析过程中所产生的一系列模拟信号值，此处的模拟信号值也就是没有受到电阻器噪声和散粒效应影响的信号值，所述指数律公式表示为：

I'_m＝I_MM^β，

其中，I'_m为同位素在混合样品质谱分析中的模拟信号值集合，M为某同位素的质量值。

步骤2.4，与此同时，根据所述步骤1.2得到的误差模型，可以得到误差的协方差矩阵，如⁵⁰Cr-⁵⁴Cr双稀释剂中，⁵⁰Cr、⁵²Cr、⁵³Cr、⁵⁴Cr四个同位素用于双稀释剂的分馏校正，所以建立4×4误差协方差矩阵，即为：

$V=\left(\begin{array}{cccc}{\mathrm{\σ}}_1^2& 0& 0& 0\\ 0& {\mathrm{\σ}}_2^2& 0& 0\\ 0& 0& {\mathrm{\σ}}_3^2& 0\\ 0& 0& 0& {\mathrm{\σ}}_4^2\end{array}\right)$

步骤2.5，根据所述步骤2.3得到的模拟信号值(I'_m)和所述步骤2.4得到的误差的协方差矩阵，利用蒙特-卡罗法模拟得到质谱分析过程中所产生的一系列新的混合样品的模拟信号值集合(I_m值)，此处的信号值同时受到电阻器噪声、散粒效应和质量歧视效应的影响，即为(I_m-A1、I_m-A2、I_m-A3、I_m-A4)^T，它们分别是四个同位素^A1X、^A2X、^A3X、^A4X所对应的信号值，在Cr同位素中I_m-A1、I_m-A2、I_m-A3、I_m-A4，分别为I_m-50、I_m-52、I_m-53、I_m-54。本发明中所使用的蒙特-卡罗法是将一组均匀分布的数据转变为一组正态分布的数据，如图2所示。

由公式可以获得质谱分析过程中的模拟同位素比值，其中i、k为某元素中用于双稀释剂分馏校正的同位素序号，如选k＝1,则i＝2、3、4。

步骤3，利用鸟巢迭代算法校正所述步骤2得到的模拟测量值，并得到对应的自然分馏因子值；

如某元素具有4个或4个以上的同位素时，以使用其中4个同位素为例，分别假定为^A1X,^A2X,^A3X,^A4X等，^A1X-^A2X为双稀释剂，^A4X/^A3X是从混合样品中提取出的自然样品的同位素比值。

所述步骤3进一步包括以下步骤：

步骤3.1，根据混合样品的同位素测量值双稀释剂的同位素比值自然样品的同位素比值通过下面公式可以得到同位素双稀释剂的试验值

${\left(\frac{X_{A1}}{X_{A2}}\right)}_{s,\mathrm{trial}}=\frac{{\left(\frac{X_{A1}}{X_{A3}}\right)}_s}{{\left(\frac{X_{A2}}{X_{A3}}\right)}_s}\×\frac{{\left(\frac{X_{A1}}{X_{A3}}\right)}_m-{\left(\frac{X_{A1}}{X_{A3}}\right)}_n}{{\left(\frac{X_{A1}}{X_{A3}}\right)}_s-{\left(\frac{X_{A1}}{X_{A3}}\right)}_m}\×\frac{{\left(\frac{X_{A2}}{X_{A3}}\right)}_s-{\left(\frac{X_{A2}}{X_{A3}}\right)}_m}{{\left(\frac{X_{A2}}{X_{A3}}\right)}_m-{\left(\frac{X_{A2}}{X_{A3}}\right)}_n},$

其中，混合样品的同位素测量值和双稀释剂的同位素比值均由步骤2获得，自然样品的同位素比值可由标准溶液的同位素组成替代。

步骤3.2，根据得到的试验值[(^A1X/^A2X)_s,trial，得到仪器的质量歧视因子β：

β＝ln((^A1X/^A2X)_s/(^A1X/^A2X)_s,trial)/ln(M_A1/M_A2)，

其中，M_A1、M_A2、为^A1X，^A2X同位素的质量数；

步骤3.3，根据仪器的质量歧视因子β值对于混合样品的同位素比值进行校正：

M_Ai＝m_Ai(M_Ai/M_Ak)^-β，

其中，M_Ai为校正后混合样品的同位素比值，m_Ai为校正前混合样品的同位素比值，M_Ai、M_Ak为^AiX、^AkX同位素的质量数。

步骤3.4，将校正后的混合样品的同位素比值代替校正前的混合样品的同位素比值，并基于上述同位素双稀释剂的试验值的求取公式，重新得到同位素双稀释剂的试验值(^A1X/^A2X)_s,trial，然后进行重复迭代直至β＝0或收敛，同时混合样品的校正值也收敛；

步骤3.5，根据最新的混合样品的同位素测量值、双稀释剂的同位素比值和自然样品的同位素比值，得到新的自然样品的同位素试验值(^A4X/^A3X)_n,trial：

${\left(\frac{X_{A4}}{X_{A3}}\right)}_{n,\mathrm{trial}}={\left(\frac{X_{A4}}{X_{A3}}\right)}_s+({\left(\frac{X_{A4}}{X_{A3}}\right)}_m-{\left(\frac{X_{A4}}{X_{A3}}\right)}_s)\×\frac{{\left(\frac{X_{A1}}{X_{A3}}\right)}_n-{\left(\frac{X_{A1}}{X_{A3}}\right)}_s}{{\left(\frac{X_{A1}}{X_{A3}}\right)}_m-{\left(\frac{X_{A1}}{X_{A3}}\right)}_s};$

步骤3.6，根据所述步骤3.5得到的自然样品的同位素试验值(^A4X/^A3X)_n,trial得到相对于标准样品的同位素分馏值α，也即自然分馏因子：

α＝ln((^A4X/^A3X)_n/(^A4X/^A3X)_n,trial)/ln(M_A4/M_A3)，

其中，M_A3、M_A4、为^A3X,^A4X同位素的质量数。

进而利用下式得到自然样品中其他同位素的试验值：

S_Ai＝s_Ai(M_Ai/M_Ak)^-α，

其中，S_Ai为校正后混合样品的同位素比值，s_Ai为标准溶液的同位素比值。

最后将求得的自然样品的同位素试验值代替原来的自然样品同位素试验值。

步骤3.7，重复步骤3.1至步骤3.6进行迭代，直至自然样品的同位素比值或α值收敛，完成对于模拟信号值的校正。

步骤4，计算自然分馏因子的误差；

所述步骤4进一步包括以下步骤：

步骤4.1，得到混合样品中的真实同位素比值；

该步骤中，根据步骤2获得的混合样品无电阻器噪声和散粒效应时的模拟信号值，经公式计算得到混合样品中的真实同位素比值，其中，I'_m-Ai、I'_m-Ak为混合样品中^AiX、^AkX无电阻器噪声和散粒效应干扰而产生的信号值，如选k＝1,则i＝2、3、4。

另外，将经步骤3校正得到的自然分馏因子记为：α₁。

步骤4.2，变动N次λ、p的组合值，由所述步骤2可以得到N个混合样品在电阻器噪声和散粒效应作用下的同位素测量值集合，进而由所述步骤3得到N个自然分馏因子集合；

在本发明一实施例中，每变动一次λ、p的组合值，可由步骤2得到1000组混合样品的同位素测量值，即由1000组[(^A1X/^A3X)_m、(^A2X/^A3X)_m、(^A4X/^A3X)_m]值构成的一个集合，经步骤3分别对集合中的每一组数据进行分馏校正，得到1000个自然分馏因子，记为i由1变到1000。

步骤4.3，基于所述步骤4.1和4.2得到的结果，求取自然分馏因子的误差。

该步骤中，由公式计算得到每一个α值的误差，并进一步求出这1000个误差值的平均值E_α，作为所述自然分馏因子的误差。

步骤5，以自然分馏因子误差和双稀释剂在混合样品中的比例为对象，作二维误差曲线图，所述二维误差曲线图中最小误差点所对应的相关配比参数值即为最优配比参数值。

该步骤中，以p值为横坐标，E_α值为纵坐标，作E_α-p的二维误差曲线图，如图3、图4所示，图3为⁵⁰Cr-⁵⁴Cr双稀释剂的E_α—p二维误差曲线线图，图4为⁹⁷Mo-¹⁰⁰Mo双稀释剂的E_α—p二维误差曲线图。从图中可以看出，无论⁵⁰Cr-⁵⁴Cr双稀释剂，还是⁹⁷Mo-¹⁰⁰Mo双稀释剂E_α-p曲线，均呈现U型曲线，p值靠近1或0时误差值急剧增大，这一现象在国际上认为是Overspike和Underspike现象，也就是由双稀释剂添加过多、或过少引起的，在读图时需要将这两部分排除。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种用于多种元素同位素双稀释剂技术中获取最优配比值的数据处理方法，其特征在于，该方法包括以下步骤：

步骤1，根据质谱的仪器参数建立误差模型；

步骤2，基于所述步骤1获得的误差模型，利用蒙特-卡罗法模拟得到样品在质谱分析过程中的测量值；

步骤3，利用鸟巢迭代算法校正所述步骤2得到的测量值，并得到对应的自然分馏因子值；

步骤4，计算自然分馏因子的误差；

步骤5，以自然分馏因子的误差和双稀释剂在混合样品中的比例为对象，作二维误差曲线图，所述二维误差曲线图中最小误差点所对应的相关配比参数值即为最优配比参数值。

2.根据权利要求1所述的方法，其特征在于，所述步骤1进一步包括以下步骤：

步骤1.1，获取质谱的仪器参数；

步骤1.2，根据所述步骤1.1得到的仪器参数建立误差模型。

3.根据权利要求2所述的方法，其特征在于，所述仪器参数至少包括仪器的频宽、仪器的阻抗、仪器运行时的温度，仪器的电阻器噪声和散粒效应。

4.根据权利要求2所述的方法，其特征在于，所述误差模型表示为：

其中，i表示同位素序号，σ_{Johnson-noise}为电阻器噪声所产生的绝对误差；σ_Shot-noise为散粒效应产生的绝对误差；σ_干扰为同质异位素和同质多原子分子的干扰。

5.根据权利要求1所述的方法，其特征在于，所述步骤2进一步包括以下步骤：

步骤2.1，根据双稀释剂的同位素理想信号的比值，获取双稀释剂的同位素组成；

步骤2.2，根据混合样品的同位素理想信号的比值，获取自然样品与双稀释剂的混合样品的同位素组成；

步骤2.3，获取所述混合样品的模拟信号值；

步骤2.4，根据所述步骤1的误差模型，得到误差的协方差矩阵；

步骤2.5，根据所述步骤2.3得到的模拟信号值和所述步骤2.4得到的误差的协方差矩阵，利用蒙特-卡罗法模拟得到质谱分析过程中所产生的新的混合样品的模拟信号值集合。

6.根据权利要求5所述的方法，其特征在于，所述步骤2.3中，由指数律公式得到混合样品在质谱分析过程中所产生的模拟信号值：

I′_m＝I_MM^β，

其中，I'_m为同位素在混合样品质谱分析中的模拟信号值集，I_M为混合样品中的同位素在质谱分析过程中的理想信号集，M为某同位素的质量值，β为质谱仪器的质量歧视因子。

7.根据权利要求1所述的方法，其特征在于，所述步骤3进一步包括以下步骤：

步骤3.1，根据混合样品的同位素测量值、双稀释剂的同位素比值、自然样品的同位素比值得到同位素双稀释剂的试验值；

步骤3.2，根据得到的试验值，得到仪器的质量歧视因子β；

步骤3.3，根据仪器的质量歧视因子β值对于混合样品的同位素比值进行校正；

步骤3.4，将校正后的混合样品的同位素比值代替校正前的混合样品的同位素比值，重新得到同位素双稀释剂的试验值，然后进行重复迭代直至β＝0或收敛，同时混合样品的校正值也收敛；

步骤3.5，根据最新的混合样品的同位素测量值、双稀释剂的同位素比值和自然样品的同位素比值，得到新的自然样品的同位素试验值；

步骤3.6，根据所述步骤3.5得到的自然样品的同位素试验值得到相对于标准样品的同位素分馏值，即自然分馏因子，进而得到自然样品中其他同位素的试验值，最后将求得的自然样品的同位素试验值代替原来的自然样品同位素试验值；

步骤3.7，重复步骤3.1至步骤3.6进行迭代，直至自然样品的同位素比值或α值收敛，完成对于模拟信号值的校正。

8.根据权利要求7所述的方法，其特征在于，所述步骤3.3中，利用下式对于混合样品的同位素比值进行校正：

M_Ai＝m_Ai(M_Ai/M_Ak)^-β，

其中，M_Ai为校正后混合样品的同位素比值，m_Ai为校正前混合样品的同位素比值，M_Ai、M_Ak为^AiX、^AkX同位素的质量数。

9.根据权利要求1所述的方法，其特征在于，所述步骤4进一步包括以下步骤：

步骤4.1，得到混合样品中的真实同位素比值；

步骤4.2，变动N次某单稀释剂在双稀释剂中的比例值，以及双稀释剂在混合样品中的比例值，由所述步骤2得到N个混合样品在电阻器噪声和散粒效应作用下的同位素模拟信号值集合，由所述步骤3得到N个校正自然分馏因子集合；

步骤4.3，基于所述步骤4.1和4.2得到的结果，求取自然分馏因子误差。

10.根据权利要求9所述的方法，其特征在于，所述步骤4.3中，首先利用所述步骤3校正前后的自然分馏因子值得到每一个自然分馏因子的误差，并通过所述步骤4.2得到N个自然分馏因子误差值的平均值E_α，作为所述自然分馏因子的误差。