本發(fā)明涉及元素測試領(lǐng)域，尤其涉及一種測定硼同位素比值的方法。

背景技術(shù)：

硼具有兩種穩(wěn)定的同位素B¹⁰和B¹¹，它們的相對豐度分別約為20％和80％。由于自然界中硼同位素組成變化很大，而且在水體中硼是易溶元素，因此硼同位素被廣泛地應(yīng)用于多種地質(zhì)環(huán)境下地球化學過程的示蹤研究?；谪撾x子熱電離飛行時間質(zhì)譜(Negative Ion Thermal Ionization Time of Flight Mass Spectrometer:NITI-TOF-MS)的分析技術(shù)是一種簡單、快捷的光譜測試和分析手段，可實現(xiàn)硼元素的實時測量，更為顯著的特點是其使用的樣品制備流程簡單，處理時間短。

但是，這種分析方法存在的不足為：由于儀器分辨率限制，B¹⁰和B¹¹兩個相鄰峰相互干擾、重疊，嚴重影響硼同位素相對豐度的準確計算。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于提供一種測定硼同位素比值的方法，從而解決現(xiàn)有技術(shù)中存在的前述問題。

為了實現(xiàn)上述目的，本發(fā)明所述測定硼同位素比值的方法，所述方法包括：

S1，獲取含硼樣本的原始質(zhì)譜圖，然后采用主成分分析法去除異常樣本，得到初始譜圖；

S2，采用S-G平滑濾波器對初始譜圖進行平滑處理，得到平滑譜圖；

S3，利用高斯線形函數(shù)對平滑譜圖進行曲線擬合，得到硼各同位素的峰特征值；

S4，根據(jù)得到的硼各同位素的峰特征值計算硼同位素豐度。

優(yōu)選地，步驟S2，具體按照下述步驟實現(xiàn)：

將初始譜圖劃分成多塊局部區(qū)域，然后采用最小二乘算法基于公式(1)依次對每個局部區(qū)域中的數(shù)據(jù)進行平滑處理：

式(1)中f表示根據(jù)已知的擬合多項式階數(shù)得到的光譜圖模型，n代表已知的擬合多項式階數(shù)，b_q表示未知的第q個擬合系數(shù)，i^q表示第q個擬合變量；q＝0,1,2,……,n。

更優(yōu)選地，為了使平滑處理后得到的譜圖與初始譜圖誤差最小，根據(jù)公式(2)對初始譜圖的平滑處理：

式(2)中，x(i)表示初始譜圖中的橫坐標為i的數(shù)據(jù)點的特征值數(shù)據(jù)，s表示在局部區(qū)域中擬合數(shù)據(jù)長度的一半。

優(yōu)選地，步驟S3，具體按照下述步驟實現(xiàn)：

S31，假設(shè)平滑譜圖含有m個實驗觀察值和p個高斯子峰；

由于平滑譜圖中的每個高斯子峰含有3個參數(shù)，所述參數(shù)包括峰強度α，中心頻率ν₀，半峰寬度σ，故p個高斯子峰共有t＝3p個待定參數(shù)，所以根據(jù)高斯線型函數(shù)，所述平滑譜圖的光譜擬合函數(shù)G(v,β₁,…,β_t)表示為：

式(3)中，v表示頻率，β₁,…,β_t表示t個待定參數(shù)，G_k(v,α_k,σ_k,v_0,k)表示第k個高斯子峰的光譜擬合函數(shù)，k＝1,2……,p；α_k表示第k個高斯子峰的峰強度，σ_k表示第k個高斯子峰的半峰寬度，v_0,k表示第k個高斯子峰的中心頻率；

S32，獲取平滑譜圖中的譜帶與公式(3)得到的擬合譜帶間的殘差平方和E的表達式(4)：

公式(4)中，y_r，v_r分別是第r個實驗觀察點的強度與頻率的觀察值，r＝1,2……,m；

S33，因，當殘差平方和E最小化時，得到的擬合譜帶最接近初始譜圖，故，將公式(4)轉(zhuǎn)換得到的公式(5)，求取當殘差平方和E最小化時，t個待定參數(shù)β₁,……β_t，即得到硼各同位素的峰特征值；

其中，β_j＝β₁,……β_t。

更優(yōu)選地，步驟S33中，求取當殘差平方和E最小化時，t個待定參數(shù)β₁,……β_t，具體為：在擬合譜帶最接近初始譜圖的條件下，采用Levenberg–Marquardt求解最優(yōu)化的t個待定參數(shù)：

A1，設(shè)定β_j的初始值，并記為設(shè)初始值與真實初始值之差為δ_j，表示為公式(6)：

A2，為了確定δ_j，將函數(shù)G(ν_i,β)在初值處一階泰處勒展開，得到公式(7)：

其中，公式(7)中ν_l表示第l個實驗觀察點的頻率觀察值；

A3，將G(ν_l,β_j)的一階導數(shù)逼近代入(5)式，并記G_l0＝G(ν_l,β_j⁽⁰⁾)，同時令得到公式(8)：

將公式(8)展開得到關(guān)于δ_j＝(δ₁,…,δ_t)^T的線性方程組(9)：

其中，設(shè)定：

在方程組(9)中引入了阻尼因子λ，得到：

其中，I為單位陣；

將線性方程組(11)中的單位陣I替換成海塞矩陣A^TA的對角線元素組成的對角陣，得到了設(shè)初始值與真實初始值之差為δ_j的計算公式(12)：

A4，將得到的δ_j其代入(6)式，可以得到逼近真實初始值的β_j⁽¹⁾；

A5，將β_j⁽¹⁾的值賦給β_j⁽⁰⁾，按照A1至A4的算法進行迭代計算得到下一個逼近真實初始值的參數(shù)值。

更優(yōu)選地，迭代的次數(shù)至少為10次。

更優(yōu)選地，當次迭代后計算得到的殘差平方和E_當次與前一次迭代后計算得到的殘差平方和E_前一次相比減小明顯時，則后一次迭代運算時的阻尼因子λ_后一次小于當次迭代運算時的阻尼因子λ_當次；

當次迭代后計算得到的殘差平方和E_當次與前一次迭代后計算得到的殘差平方和E_前一次相比減小不明顯時，則后一次迭代運算時的阻尼因子λ_后一次大于當次迭代運算時的阻尼因子λ_當次。

更優(yōu)選地，A1，設(shè)定β_j的初始值，具體按照下述步驟實現(xiàn)：

利用二階或更高階導數(shù)譜逼近譜峰并得到譜峰的中心頻率ν₀，譜峰符合條件公式(13)：

B⁽¹⁾(ν)＝0,B⁽²⁾(ν)＜0 (13)

其中，B⁽¹⁾(ν)代表平滑譜圖的一階導數(shù)譜，B⁽²⁾(ν)代表平滑譜圖的二階導數(shù)譜，ν代表頻率；

峰強度α的初始值估計：待p個高斯子峰的初始中心位置確定后，利用線性插值方法估計平滑譜圖B(ν)在中心頻率ν₀處的值，并將得到的值作為峰的強度α的初始估計值：

α_m＝B(ν_0m),m＝1,…,p (14)

半峰寬度σ的初始值估計：從譜峰的中心頻率ν₀開始向兩側(cè)尋找平滑譜圖B(ν)的峰強度值降至與α/2最近的兩個點ν₁,ν₂，并分別計算ν₁,ν₂與ν₀的距離，取距離最小者作為半峰寬度σ的初始估計值。

本發(fā)明的有益效果是：

本發(fā)明所述方法是基于Levenberg–Marquardt(L-M)的曲線擬合算法將B¹⁰和B¹¹重疊峰剝離開，去除相互干擾影響，實時計算出峰位點對應(yīng)的峰高、峰面積等特征參數(shù)。相比傳統(tǒng)方法中直接定位峰位點的計算方法，本發(fā)明所述方法有效提高硼同位素相對豐度的測定準確率，可望為地球化學過程的示蹤研究提供有力保障。

附圖說明

圖1是實施例中的含硼樣本的原始質(zhì)譜圖；

圖2是實施例中得到的平滑譜圖；

圖3是實施例中利用高斯線形函數(shù)對平滑譜圖進行曲線擬合后得到擬合譜圖；

圖4是實施例中含硼各同位素特征參數(shù)的譜圖；

圖5是實施例中測定硼同位素比值的方法的流程圖。

具體實施方式

為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點更加清楚明白，以下結(jié)合附圖，對本發(fā)明進行進一步詳細說明。應(yīng)當理解，此處所描述的具體實施方式僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

實施例

本實施例所述測定硼同位素比值的方法，包括：

S1，獲取含硼樣本的原始質(zhì)譜圖，然后采用主成分分析法去除異常樣本，得到初始譜圖，參照圖1；

采用核工業(yè)北京地質(zhì)研究院分析測試研究所自主研發(fā)的負離子熱電離飛行時間質(zhì)譜儀(NITI-TOF-MS)采集數(shù)據(jù)，實驗獲得3000個樣本的質(zhì)譜圖。

S2，采用S-G平滑濾波器對初始譜圖進行平滑處理，得到平滑譜圖，參照圖2；

S3，利用高斯線形函數(shù)對平滑譜圖進行曲線擬合，得到硼各同位素的峰特征值，參照圖3；

S4，根據(jù)得到的硼各同位素的峰特征值計算硼同位素豐度，參照圖4。

更詳細的解釋說明：

(一)通常情況下，通過NITI-TOF-MS采集的初始譜圖會受到噪聲干擾，導致出現(xiàn)一些“假峰”，這些特征會影響后續(xù)硼同位素豐度計算。采用S-G平滑濾波器(Savitzky-Golay filter)先對初始譜圖平滑，去除不相關(guān)噪聲。S-G平滑濾波器基于局部區(qū)域最小二乘算法，通過多項式平滑局部窗中數(shù)據(jù)。

步驟S2，具體按照下述步驟實現(xiàn)：將初始譜圖劃分成多塊局部區(qū)域，然后采用最小二乘算法基于公式(1)依次對每個局部區(qū)域中的數(shù)據(jù)進行平滑處理：

式(1)中f表示根據(jù)已知的擬合多項式階數(shù)得到的光譜圖模型，n代表已知的擬合多項式階數(shù)，b_q表示未知的第q個擬合系數(shù)，i^q表示第q個擬合變量；q＝0,1,2,……,n。

其中，為了使平滑處理后得到的譜圖與初始譜圖誤差最小，根據(jù)由公式(1)轉(zhuǎn)換的公式(2)對初始譜圖的平滑處理：

式(2)中，x(i)表示初始譜圖中的橫坐標為i的數(shù)據(jù)點的特征值數(shù)據(jù)，s表示在局部區(qū)域中擬合數(shù)據(jù)長度的一半。

(二)對采集的樣品譜圖，為了準確計算B¹⁰和B¹¹的峰位，利用曲線擬合將原始譜圖解析為一個或一個以上單峰之和，并且使單個子峰加起來擬合的譜圖與實測譜圖之間誤差平方和達到最小。通常利用高斯線形函數(shù)擬合譜圖數(shù)據(jù)。高斯線型函數(shù)G(ν)可寫成：

其中，α是高斯子峰的強度，ν是頻率點，ν₀是中心頻率，σ是半峰寬度。

步驟S3，具體按照下述步驟實現(xiàn)：

S31，假設(shè)平滑譜圖含有m個實驗觀察值和p個高斯子峰；

由于平滑譜圖中的每個高斯子峰含有3個參數(shù)，所述參數(shù)包括峰強度α，中心頻率ν₀，半峰寬度σ，故p個高斯子峰共有t＝3p個待定參數(shù)，所以根據(jù)高斯線型函數(shù)，所述平滑譜圖的光譜擬合函數(shù)G(v,β₁,…,β_t)表示為：

式(3)中，v表示頻率，β₁,…,β_t表示t個待定參數(shù)，G_k(v,α_k,σ_k,v_0,k)表示第k個高斯子峰的光譜擬合函數(shù)，k＝1,2……,p；α_k表示第k個高斯子峰的峰強度，σ_k表示第k個高斯子峰的半峰寬度，v_0,k表示第k個高斯子峰的中心頻率；

S32，獲取平滑譜圖中的譜帶與公式(3)得到的擬合譜帶間的殘差平方和E的表達式(4)：

公式(4)中，y_r，v_r分別是第r個實驗觀察點的強度與頻率的觀察值，r＝1,2……,m；

S33，因，當殘差平方和E最小化時，得到的擬合譜帶最接近初始譜圖，故，將公式(4)轉(zhuǎn)換得到的公式(5)，求取當殘差平方和E最小化時，t個待定參數(shù)β₁,……β_t，即得到硼各同位素的峰特征值；

其中，β_j＝β₁,……β_t。

另，步驟S33中，求取當殘差平方和E最小化時，t個待定參數(shù)β₁,……β_t，具體為：在擬合譜帶最接近初始譜圖的條件下，采用Levenberg–Marquardt求解最優(yōu)化的t個待定參數(shù)：

A1，設(shè)定β_j的初始值，并記為設(shè)初始值與真實初始值之差為δ_j，表示為公式(6)：

A2，為了確定δ_j，將函數(shù)G(ν_i,β)在初值處一階泰處勒展開，得到公式(7)：

其中，公式(7)中ν_l表示第l個實驗觀察點的頻率觀察值；

A3，將G(ν_l,β_j)的一階導數(shù)逼近代入(5)式，并記G_l0＝G(ν_l,β_j⁽⁰⁾)，同時令得到公式(8)：

將公式(8)展開得到關(guān)于δ_j＝(δ₁,…,δ_t)^T的線性方程組(9)：

其中，設(shè)定：

在方程組(9)中引入了阻尼因子λ，得到：

其中，I為單位陣；

將線性方程組(11)中的單位陣I替換成海塞矩陣A^TA的對角線元素組成的對角陣，得到了設(shè)初始值與真實初始值之差為δ_j的計算公式(12)：

A4，將得到的δ_j其代入(6)式，可以得到逼近真實初始值的β_j⁽¹⁾；

A5，將β_j⁽¹⁾的值賦給β_j⁽⁰⁾，按照A1至A4的算法進行迭代計算得到下一個逼近真實初始值的參數(shù)值。迭代的次數(shù)至少為10次。

另：當次迭代后計算得到的殘差平方和E_當次與前一次迭代后計算得到的殘差平方和E_前一次相比減小明顯時，則后一次迭代運算時的阻尼因子λ_后一次小于當次迭代運算時的阻尼因子λ_當次；

當次迭代后計算得到的殘差平方和E_當次與前一次迭代后計算得到的殘差平方和E_前一次相比減小不明顯時，則后一次迭代運算時的阻尼因子λ_后一次大于當次迭代運算時的阻尼因子λ_當次。

其中，A1，設(shè)定β_j的初始值，具體按照下述步驟實現(xiàn)：

利用二階或更高階導數(shù)譜逼近譜峰并得到譜峰的中心頻率ν₀，譜峰符合條件公式(13)：

B⁽¹⁾(ν)＝0,B⁽²⁾(ν)＜0 (13)

其中，B⁽¹⁾(ν)代表平滑譜圖的一階導數(shù)譜，B⁽²⁾(ν)代表平滑譜圖的二階導數(shù)譜，ν代表頻率；

峰強度α的初始值估計：待p個高斯子峰的初始中心位置確定后，利用線性插值方法估計平滑譜圖B(ν)在中心頻率ν₀處的值，并將得到的值作為峰的強度α的初始估計值：

α_m＝B(ν_0m),m＝1,…,p (14)

半峰寬度σ的初始值估計：從譜峰的中心頻率ν₀開始向兩側(cè)尋找平滑譜圖B(ν)的峰強度值降至與α/2最近的兩個點ν₁,ν₂，并分別計算ν₁,ν₂與ν₀的距離，取距離最小者作為半峰寬度σ的初始估計值。

對圖3中的圖譜進行運算，B¹⁰、B¹¹完全分離開。同時獲得相應(yīng)特征參數(shù)如圖4所示，其中，B¹⁰對應(yīng)峰位9.976、峰高9139、峰寬0.2637、峰面積2565；B¹¹對應(yīng)峰位11.04、峰高40060、峰寬0.265、峰面積11300。

直接尋找極大值點作為峰位的常規(guī)方法與本實施例高斯擬合算法計算硼同位素比值結(jié)果如

表1表1所示，實驗結(jié)果表明：相比常規(guī)計算方法，利用高斯擬合質(zhì)譜圖獲得的硼同位素比值精密度更高(小于1％)。

表1常規(guī)方法與高斯擬合算法計算硼同位素比值結(jié)果

綜上，本實施例為了解決硼同位素(B¹⁰和B¹¹)豐度測定問題，采用核工業(yè)北京地質(zhì)研究院分析測試研究所自主研發(fā)的負離子熱電離飛行時間質(zhì)譜儀(NITI-TOF-MS)采集數(shù)據(jù)。實驗獲得3000個樣本的質(zhì)譜圖，首先應(yīng)用主成分分析剔除異常樣本點，接著采用S-G平滑濾波器消除噪聲影響。為了準確獲得B¹⁰和B¹¹的峰位、峰高、半峰寬、峰面積等參數(shù)，利用高斯線型對該質(zhì)譜圖擬合，采用Levenberg–Marquardt(L-M)算法求解優(yōu)化問題獲取特征。最后根據(jù)特征信息計算硼同位素豐度，硼同位素(B¹⁰和B¹¹)豐度測定流程如圖5所示。

通過采用本發(fā)明公開的上述技術(shù)方案，得到了如下有益的效果：

本發(fā)明所述方法是基于Levenberg–Marquardt(L-M)的曲線擬合算法將B¹⁰和B¹¹重疊峰剝離開，去除相互干擾影響，實時計算出峰位點對應(yīng)的峰高、峰面積等特征參數(shù)。相比傳統(tǒng)方法中直接定位峰位點的計算方法，本發(fā)明所述方法有效提高硼同位素相對豐度的測定，可望為地球化學過程的示蹤研究提供有力保障。

以上所述僅是本發(fā)明的優(yōu)選實施方式，應(yīng)當指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明原理的前提下，還可以做出若干改進和潤飾，這些改進和潤飾也應(yīng)視本發(fā)明的保護范圍。