Pb-Sr-Nd-Hf同位素參數計算及程序設計

路遠發，李文霞

（長江大學資源與環境學院，武漢 430100）

1 前言

Pb、Sr、Nd、Hf同位素都是重穩定同位素，地質過程中不發生同位素分餾，因而其同位素組成只與來源有關而與過程無關。同時，這些穩定同位素也都是放射性成因同位素，其同位素組成與所處地質環境地球化學特征（如U/Pb比、Sm/Nd比和Lu/Hf比等）及地質演化歷史有關。由于不同區域、不同構造環境的巖石/礦物地球化學特征及演化歷史不同，因此這些同位素組成就成了物質來源與構造環境判別的重要標志，在研究殼幔分離與殼幔相互作用中起著不可替代的作用，倍受廣大地學工作者的重視，使之成為當前地球科學的熱點之一，相關成果與研究論文不勝枚舉。

隨著TIMS、LA-MC-ICP-MS、離子探針技術的發展與普及，當前同位素分析測試技術已經有了很大的提高。但是，同位素實驗室給出的結果(同位素比值)所反映的地質信息非常有限，更多的信息(參數)需要通過數據處理來提取。因此，同位素數據的精度不僅僅取決于實驗室的測試精度，同時也取決于數據處理是否正確。不正確的處理方法當然得不到正確的結果，而用這些不準確、甚至錯誤的數據來討論地質作用，就很難保證研究結果的正確性。然而，當前我國地學期刊中的同位素數據處理水平參差不齊，存在不少錯誤，如個別論文中同位素初始比值比測試值還要高、鍶同位素初始比值(校正值)甚至低于BABI值等。當然，這些顯而易見的錯誤是極少數的，更為普遍的是那些難以識別的隱性錯誤(雖與正確的結果有一定的差異，但不經過驗算不易察覺)。究其原因，是目前國內還缺少一個統一的數據處理軟件，各位專家多在按自己的理解、用自己的方法處理數據。由于同位素數據處理涉及相對復雜的公式推導，而該學科的成熟與普及時間并不長，許多科學家并未系統學習過同位素地球化學的理論，因而在理解上難免會出現偏差，導致數據處理錯誤。為此，本文對鉛-鍶-釹-鉿同位素的相關參數以及這些參數的計算公式作一個全面的介紹，在此基礎上，給出筆者設計的計算機程序，希望能給國內學者提供一個統一的數據處理平臺。

2 計算原理

2.1 鉛(Pb)同位素參數

2.1.1 鉛同位素組成與增長方程

自然界中鉛有8種同位素，其中4種是放射性同位素，4種是穩定同位素。4種放射性鉛同位素是210Pb、211Pb、212Pb和214Pb，它們分別是238U、235U、232Th三個衰變系列的中間產物，即鈾系中的210Pb和214Pb，錒系中的211Pb和釷系中的212Pb。除210Pb的半衰期較長(T1/2=22.3a)外，其他幾個放射性同位素半衰期都很短。這些放射性鉛同位素的豐度都很低，相對于穩定的鉛同位素可以忽略不計。

鉛 的 四 種 穩 定 同 位 素 是204Pb、206Pb、207Pb、208Pb，其中204Pb是非放射成因同位素，自形成至今，其同位素豐度保持不變。206Pb、207Pb、208Pb分別是238U、235U和232Th衰變系列的最終產物，豐度隨年齡和U、Th、Pb的含量變化而變化。

相應的206Pb、207Pb、208Pb的增長方程為：

式中下標S代表樣品的測定值（可以省略），i代表初始值，λ238代表238U的衰變常數，余類推，t代表礦物(巖石)的形成年齡，下同。

2.1.2 鉛同位素主要參數

2.1.2.1 μ值和υ值

定義式為:μ=238U/204Pb，υ=235U/204Pb

由(4)式可知，在一個μ值一定的體系中，從地球形成演化到現在鉛的同位素組成應為：

式中a0為地球形成時的206Pb/204Pb比值，下文中的b0和c0則分別代表地球形成時的207Pb/204Pb和208Pb/204Pb比值，T為地球的形成年齡。如果從t時刻開始演化到現在，則有：

由于鉛同位素研究的是普通鉛，即是t時刻的鉛同位素組成，即從地球形成演化到t時刻的鉛(并不是從地球形成演化到現在積累的鉛)，所以有：

對于普通鉛(不含U的礦物、巖石，或雖含U但與鉛相比可以忽略不記)來說：

(206Pb/204Pb)t=(206Pb/204Pb)S，因此：

類似公式（9），我們可以得到：

因238U/235U=137.88，因此υ=μ/137.88, υ不具有獨立性，完全可以用μ來代替，現在已很少使用。

2.1.2.2 ω值

ω的定義式為：ω=232Th/204Pb

類似于公式(9)，我們可以得到

2.1.2.3 Th/U比值

首先，應該知道，同位素比值是摩爾分數比，而Th/U比值是質量分數比，因此計算過程中要特別注意將U、Th的摩爾分數轉換成質量分數w(U)和w(Th)。

由于U主要由235U、238U兩個同位素構成，且238U/235U=137.88，其原子量為238.03，由此可得：

Th雖然有6種天然放射性同位素：232Th、228Th、230Th、234Th、227Th和231Th。除232Th外，其余5種為238U、235U和232Th 三個衰變系列的中間產物。由于這些放射性同位素的豐度值非常低，半衰期又十分短，因此一般將232Th的豐度值作為100%。

2.1.2.4 單階段模式年齡(tCDT)

模式年齡作為成巖-成礦年齡來使用已經不被認可，但仍不失為一個重要的鉛同位素參數，如華南的鉛鋅礦床，成因不同其模式年齡也不同[1]。因此，模式年齡特別是單階段模式年齡仍是學者們常用的參數之一。

單階段模式年齡是從地球形成時的原始鉛開始，按一定的μ值進行演化到脫離該U-Pb體系(不再接受放射性鉛的積累)的年齡。地球原始鉛取自Canyon Diablo鐵隕石中的隕硫鐵(參數參見表1)，因此本文用tCDT表示單階段模式年齡，以便與兩階段模式年齡相區別。

聯立公式(9)和(11)，得

由上式可以看出，只要測得206Pb/204Pb和207Pb/204Pb值，就可以算出模式年齡tCDT。式(16)是一個超越方程，手工是無法獲得其解的，但用計算機計算已不再是難題。

2.1.2.5 兩階段模式年齡(tSK)

Stacey 和Kraners[2]提出了Pb同位素兩階段演化模式。他們認為第一階段從45.7億年開始至37億年，在37億年時地球發生分異事件，從37億年至現在為第二階段。

類似于單階段模式年齡，我們可以獲得兩階段模式年齡(tSK)的計算公式如下：

該公式與單階段模式年齡公式基本一致，所不同的是起始時間和初始比值不同，相關參數詳見表1。

2.1.2.6 △α、△β、△γ

△α，△β，△γ分 別 代 表206Pb/204Pb和207Pb/204Pb和208Pb/204Pb相對于同時代地幔的相對偏差，是研究鉛的來源和構造環境的重要參數[3]，其計算公式如下：

式中右下角標，S為普通鉛的測定結果，M為同時期地幔鉛。某一時間(t)地幔值可根據現代地幔值(參見表1)，按公式(4)～(6)進行計算。關于年齡t，計算時優先使用給定的年齡值，如果缺少年齡值，則以單階段模式年齡值代之，如果單階段模式年齡值＜0，則取t=0Ma，其中以給定年齡的計算結果最為可靠。如果缺少樣品的形成年齡，也可大致估計一個年齡值，如已知樣品采自燕山早期花崗巖體，可估計年齡160 Ma，可能比使用模式年齡更可靠。

2.1.2.7 V1、V2

鉛同位素組成的一大特點就是同時有206Pb/204Pb、207Pb/204Pb、208Pb/204Pb三組同位素比值，轉換成相對于同時代地幔鉛的千分差為△α、△β和△γ。在普通的直角坐標圖中，一次只能使用兩個變量。朱炳泉等[3]、常向陽等[4,5]通過拓撲分析，將△α-△β-△γ三維坐標圖中的點投影到下列回歸平面中：

在該回歸平面中△α，△β的投影點的坐標變換為：

將△值回歸平面強制通過0點(即地幔源的點)，可得a=0,b=2.0367,c=-6.143

在該回歸平面中進一步建立通過源點的新二維直角坐標系(V1，V2)，則新坐標值為[5]：

將常數a、b、c的數值代入到公式(19)和(20)中，可得如下近似計算式：

2.1.2.8 初始比值(206Pb/204Pb)i、(206Pb/204Pb)i、(206Pb/204Pb)i

初始比值，也稱校正值，是扣除礦物(巖石)形成后放射性成因鉛后的鉛同位素比值。由公式(4)～(6)可得：

由上式可以看出，計算初始值關鍵是要求得238U，235U，232Th和204Pb的值。計算步驟如下：

①n(204Pb)值(簡作204Pb，下同)的計算

鉛的四個同位素之間遵循下列關系

式中p(204Pb)，p(206Pb)，p(207Pb)，p(208Pb)分別代表204Pb，206Pb，207Pb，208Pb在總鉛中所占的比例(摩爾分數)，由公式(19)可得：

②235U、238U的計算

③232Th的計算

由上述三步，我們就可以求得樣品的235U/204Pb、238U/204Pb和232Th/204Pb，然后根據公式(25)～(27)求得各同位素比值的初始值。至于公式(25)～(27)中的t，一定要用樣品形成的年齡而不是模式年齡。

2.2 鍶(Sr)同位素參數

2.2.1 Sr同位素組成及87Sr增長方程

自然界Rb有85Rb和87Rb兩種同位素，分別占72.15%和27.85%。Sr有88Sr、87Sr、86Sr和84Sr 四個同位素，其所占的比例分別為：82.56%，7.02%，9.86%和0.56%；其中87Sr是放射性成因的，其由母體同位素87Rb通過一個β-衰變而來：

其對應的增長方程為：

2.2.2 Sr同位素組成的主要參數

2.2.2.1 初始比值(87Sr/86Sr)i

由公式(20)可得：

2.2.2.2 εSr

下標UR代表均一地幔庫。

2.2.2.3fRb/Sr

該參數反映Rb/Sr比值相對于均一地庫的富集與虧損程度，計算公式如下：

2.3 釹(Nd)同位素參數

2.3.1 Nd的同位素組成及143Nd增長方程

自 然 界Nd有142Nd、143Nd、144Nd、145Nd、146Nd、148Nd和150Nd七個同位素。其中144Nd具有天然α衰變，150Nd有β-衰變，但是它們的半衰期都很長(＞1015年)，自地球形成以來衰減的量是當前測試水平無法準確測出，且也是可以忽略不計。144Nd的相對豐度為23.7938%。

自然界Sm有144Sm、147Sm、148Sm、149Sm、150Sm、152Sm和154Sm七個同位素，其中147Sm的相對豐度為14.9957%。147Sm因α衰變引起子體同位素143Nd增長。

其對應的增長方程為

2.3.2 Nd同位素主要參數[6]

2.3.2.1 初始比值(143Nd/144Nd)i

由公式(32)得：

2.3.2.2 模式年齡TCHUR、TDM和T2DM

實測一組Sm-Nd同位素(143Nd/144Nd)S和(147Sm/144Nd)S, 當給定初始值后就可計算公式(33)中的年齡值，即模式年齡。如果以同時代(t)球粒隕石地幔庫(CHUR)的值作為樣品的初始值，我們會得到相對于球粒隕石地幔庫的模式年齡，即TCHUR；如果采用同時代虧損地幔庫(DM)的值作為樣品的初始值，我們就得到了虧損地幔模式年齡TDM。以下以TCHUR為例：

其中：

將(27)代入(26)得：

同理，我們可以得到相對于虧損地幔(DM)的模式年齡:

如果先計算出TDM，則T2DM可按下式計算[7]：

式 中：fSa=-1，fCC=,t為部分熔融發生的時間，即第二階段開始的時間[8-9]，通常用巖體的侵位年齡來表示。

一般來說，對于太古代巖石，由于當時的殼幔還沒有充分分離，因此我們通常以球粒隕石地幔庫為參照；元古代以后，殼幔已經完全分離，此時的地幔多為虧損地幔，因此常以虧損地幔為參照。需要注意的是球粒隕石地幔庫沒有單階段與兩階段之分，是一個理想的均一、穩定的同位素體系。

2.3.2.3 εNd

2.4 Hf同位素參數

2.4.1 Hf的同位素組成及176Hf增長方程

在自然界中Hf有174Hf、176Hf、177Hf、178Hf、179Hf與180Hf 六個同位素，Lu有175Lu和176Lu兩個同位素，其中176Lu通過一個β-衰變轉變為176Hf：

依據放射性衰變方程，可獲得Lu-Hf年代學方程：

2.4.2 Hf同位素主要參數[10]

2.4.2.1 初始比值 (176Hf/176Lu)i

由式（45）得：

2.4.2.2 模式年齡(TCHUR、TDM、TDMC)

Hf同位素模式年齡的推導過程與Nd同位素模式年齡是一致的，此處從略。

(注：上述各參數中，上標Hf是為了與Nd同位素的相關參數相區別，單獨使用時可不加上標。在文獻中很少見到，本文參照Nd的TCHUR計算公式建立了的計算公式，可用于殼幔還沒有充分分離的太古代巖石的研究)

2.4.2.3 εHf

2.5 同位素常數

上述許多公式中都使用了各種同位素的常數，表1中匯總了Pb-Sr-Nd-Hf同位素的主要常數。

表1 Pb-Sr-Nd-Hf同位素的主要常數Table 1 Main constants of Pb-Sr-Nd-Hf isotopes

3 程序設計

3.1 Pb-Sr-Nd-Hf同位素參數計算程序簡介

Pb-Sr-Nd-Hf同位素參數計算是作者開發的Geokit軟件中的一個組成部分[25]。十多年來，在廣大用戶的支持下，作者先后發布了多個Geokit升級版本，當前版本命名為GeokitPro。該版本的小號以更新的日期命名，如一個完整的版本號Geokit-Pro20210328代表2021年3月28日更新。有關同位素參數計算，在GeokitPro中均集成在“地球化學數據庫”模塊中。用戶可以點擊Geokit菜單中的相關菜單項，就可以打開數據處理工作薄，然后根據需要切換需要的工作表(圖1)。詳細的使用方法請參閱隨軟件發布的的幫助文件。

表2 兩階段模式年齡(tSK)計算實例Table 2 An Example for two-stage model ages of lead isotopes

圖1 地球化學數據庫的工作薄的工作表標簽Fig. 1 The menu of Geochemical database

進入到工作表后，用戶可以在工作表中直接進行相關數據庫管理與數據處理工作。由于工作表有很多快速輸入功能，也可以使用復制-粘貼的功能，因此錄入數據較工作窗口方便得多，數據準備工作可以在表格中進行，用戶可以直接在相應的表格中按要求輸入相關數據(少量的數據修改也可以在窗口中進行)。如果數據庫中已有相關數據，可以直接調入到工作表中，方法是點擊相關工作表中第一行左側的圖形(按鈕)，彈出相應的工作窗口，在窗口中先選擇項目、再選擇數據集，選擇數據集后，該數據集下的數據就會調入到工作表中。圖2是Sr-Nd同位素工作窗口，鉛同位素、鉿同位素工作窗口的結構與之類似。當在工作表中準備好數據后，點擊窗口或工作的按鈕“參數計算”，即可完成當前工作表中的相關參數的計算，并將計算結果顯示在工作表中。

圖2 Sr-Nd同位素數據管理與參數計算窗口Fig. 2 The window of data manage and parameter calculating for Sr-Nd isotopes

3.2 計算實例與結果檢驗

3.2.1 鉛同位素參數

⑴ 輸入參數：必需參數 (206Pb/204Pb)S、(207Pb/204Pb)S、(208Pb/204Pb)S（下標S代表樣品測試值，以下省略）。

可選參數(計算初始值時為必須)：參考年齡 (t/Ma)，U，Th，Pb含量。

⑵ 輸出參數：206Pb/207Pb、tCDT、tSK、μ、ω、Th/U、V1、V2、△α、△β、△γ、(206Pb/204Pb)i、(207Pb/204Pb)i、 (208Pb/204Pb)i。

此外，為保證數據的完整性，在數據庫中需要錄入樣品信息等相關數據（下同）。

不同的作者因研究目的不同，使用的參數也有所不同。因此，這些參數通常不會出現在同一篇論文中。以下我們選擇了幾篇論文中的數據，分別對這些參數進行驗證。

⑶ 計算實例

①模式年齡

模式年齡包括單階段模式年齡(tCDT)和兩階段模式年齡(tSK)。由式（16）和式（17）可知，兩個模式年齡的公式是一樣的，只是參數不同。我們利用Stacey and Kramers[2]原文中的前6組數據對本軟件進行考核，可以發現，所有計算結果均只在末位存在微小差異，這種差異可能與計算過程中所采用不同的數據精度有關，而不是由計算錯誤所引起。因此，可以認為本文的計算結果與原作者的計算結果是一致的。本軟件除了Stacey and Kramers[2]的兩階段模型外，還設計了用戶自定義兩階段模型，用戶只需要給出自己樣品的第二階段起始時間和第二階段的初始值即可。

② V1、V2

由于V1、V2兩個參數是在△α、△β、△γ三個參數的基礎上完成的，因此驗證了V1、V2的正確性也就自然驗證了△α、△β、△γ三個參數的正確性。Geokit對V1、V2兩參數的計算采用的是公式(21)和(22)。表3中的計算實例數據取該參數建立者的原文[4-5]，由對比可以看出，本文的計算結果與原文是一致的。

表3 V1、V2參數計算實例Table 3 An Example for V1 and V2 of lead isotopes

③ 初始值計算

由于大多數學者在研究鉛同位素時多采用礦石鉛，不需要進行鉛同位素校正，因此在鉛同位素研究中，初始值計算并不多。表4中的計算實例數據為秦嶺古MORB型巖石的鉛同位素組成[26]。作者原文中有兩組年齡值，一組是本文引用的MORB形成年齡（1000 Ma），另一組是變質年齡(350 Ma)。從表4可以看出：用本軟件的計算結果與原文給出的結果都是一致的。

表4 初始鉛計算實例Table 4 An Example for Initial ratios of lead isotopes

我們對高庚等[27]報道的安徽銅陵白芒山輝石閃長巖體的鉛同位素數據也進行了驗算，結果也是除個別數據在末位有±1的差異外，絕大部分計算結果是一致的。

3.2.2 Sr-Nd同位素參數計算

⑴ 輸入參數

必 要 參 數：87Rb/86Sr，87Sr/86Sr，147Sm/144Nd，143Nd/144Nd，t。

選擇參數：Rb、Sr、Sm、Nd，當缺失87Rb/86Sr時，Rb、Sr為必須參數，當缺失147Sm/144Nd時，Sm、Nd為必須參數。

⑵ 輸出參數

(87Sr/86Sr)i，εSr(0)，εSr(t)，f Rb/Sr，(143Nd/144Nd)i，TDM，T2DM，Tchur，εNd(0)，εNd(t)，fSm/Nd。其中Tchur一般用于殼幔還沒有明顯分離的太古代巖石。

⑶ 計算實例

表5和表6中的數據為拉薩地塊西段中新世查加寺鉀質火山巖的Sr-Nd同位素組成[28]，計算時年齡采用24 Ma。限于篇幅，本文只列出其中前6組數據（下列各例相同）。

表5 Sr同位素參數計算實例Table 5 An Example for Sr isotope parameters

表6 Nd同位素參數計算實例Table 6 An Example for Nd isotope parameters

另外，對田世洪等[29]所給青海玉樹東莫扎抓鉛鋅礦床Sr-Nd同位素數據進行驗證，本文給出的結果與原文也是一致的。

3.2.3 Hf同位素參數計算

⑴ 輸入參數

(176Lu/177Hf)S，(76Hf/177Hf)S，t。此外，樣品信息及元素含量不參與計算，不作為計算的輸入參數，但為保證數據的完整性，在數據庫中需要錄入相關信息及數據。

⑵ 輸出參數

(176Hf/177Hf)i，TDM，TDMC，TCHUR，εHf(0)，εHf(t)，fLu/Hf。

⑶ 計算實例

表7中的數據為西藏岡底斯東部察隅花崗巖體(CY1-01)鋯石部分Hf同位素數據[30]。與原文數據對比，除個別初始值、TDM、TDMC在末位有±1的差異外，其它數值完全一致。

表7 Hf同位素參數計算實例Table 7 An Example for Hf isotope parameters

以上計算實例及其與前人計算結果的對比表明，本軟件設計是合理的，計算結果是正確的。

4 結語

許多地質-地球化學研究的最終結論來源于同位素地球化學數據的支持，因此同位素地球化學數據處理的正確與否關系到這些結論的準確與正確性，必須引起有關科學家及相關期刊編輯部的高度重視。

本文系統地給出了Pb-Sr-Nd-Hf各同位素參數的計算公式及其導出過程，在此基礎上作者編寫了一個計算機軟件，并將本軟件的數據處理結果與前人的數據處理結果進行了對比，結果證明程序設計是合理的，計算結果是正確的。希望該軟件能成為國內學者處理相關數據的工具，避免再出現數據處理錯誤。