周統(tǒng) 邱昆峰, 2 王瑀 于皓丞 侯照亮

1. 中國地質大學(北京)地球科學與資源學院，北京 100083 2. 中國地質大學地質過程與礦產資源國家重點實驗室，北京 100083 3. 維也納大學地質系，維也納 1090

磷灰石廣泛分布于火成巖、沉積巖和變質巖中(Chew and Donelick, 2012; Henrichsetal., 2018; 趙振華和嚴爽, 2019)，屬六方晶系，其分子式為M10(Z4O6)X2，式中M=Ca、Sr、Pb、Na、REE、Ba、Mn等；Z=P、As、Si、V、Cr、S、B、N、Ge等；X=F、OH、Cl、Br、I等。微量元素通常以類質同象的形式進入磷灰石晶格，在不同形成環(huán)境下其微量元素特征存在顯著差異(朱笑青等, 2004; Deng and Wang, 2016; Anderssonetal., 2019; O’Sullivanetal., 2020; Qiuetal., 2020)。因此磷灰石微量元素特征被視為追蹤物質來源和演化的有力工具，在反演巖石成因過程中具有關鍵作用。(Chewetal., 2011; Hughes and Rakovan, 2015; Webster and Piccoli, 2015; Dengetal., 2020b, 2021; Yuetal., 2021; 蘇建超等, 2021)。Belousovaetal.(2002)總結并提出了Sr-Y、Sr-Mn、Y-(Eu/Eu*)和(Ce/Yb)N-REE判別圖解(圖1)，用于區(qū)分花崗巖、花崗偉晶巖、輝綠巖、歪堿正長巖、二輝橄欖巖、碳酸巖、鈦鐵霞輝巖和鎂鐵質巖石等常見巖石中磷灰石母巖類型。張碩等(2018)進一步提出了(Gd/Yb)N-δEu、(La/Sm)N-(La/Yb)N和(La/Yb)N-REE圖解，用于判別碳酸巖、二輝橄欖巖、輝綠巖、正長巖、花崗巖、花崗偉晶巖、巖漿型鐵礦床、熱液流體相關的鐵礦床、榴輝巖、角閃巖和磷塊巖中的磷灰石母巖類型。

現(xiàn)研究表明，數據量和圖解端元的不足，極大制約了運用判別圖解對磷灰石母巖類型進行準確判別。Belousovaetal.(2002)的圖解僅針對巖漿巖中的磷灰石類型進行了區(qū)分判定，并且部分已知類型的磷灰石微量元素數據不能與其對應區(qū)域準確吻合，平均準確率僅為59%(表1)，如中-低級變質巖和高級變質巖中的磷灰石數據沒有與其相對應的判別類型區(qū)域，而與巖漿巖中的磷灰石數據嚴重重合；鎂鐵質火成巖和富堿性巖石的數據分布不僅與該圖解所劃定的判別區(qū)域差異大，而且兩種類型重合現(xiàn)象亦嚴重(圖1)。這些已有的經典判別圖解已經逐漸無法滿足學者們對于更多磷灰石母巖類型準確判別的需求。

地球大數據是一種典型的科學大數據，是具有空間屬性的地球科學領域大數據，它一方面具有海量、多源、多時相、異構、多尺度、非平穩(wěn)等大數據的一般性質，同時具有很強的時空關聯(lián)和物理關聯(lián)(郭華東, 2018)。近年來，微區(qū)測試技術日益成熟，地球化學數據大量積累，傳統(tǒng)分析方法已逐漸無法有效利用這些數據所攜帶的信息。數據的豐富和大數據技術的發(fā)展，為以大數據為依托的分析方法應用到地質學研究中提供了可能。通過開發(fā)出智能的數據處理方法可以極大提升現(xiàn)有數據判別的準確度、可信度和效率(周永章等, 2018)。

表1 已發(fā)表的磷灰石母巖類型判別圖解準確率

本文收集了全球已發(fā)表的1925個典型磷灰石樣品測試點數據，依照O’Sullivanetal.(2020)的分類方式將所有數據分為富堿性火成巖(ALK)、超鎂鐵質巖石(UM)、鎂鐵質火成巖(IM)、長英質花崗巖(S)、中-低級變質巖(LM)、高級變質巖(HM)和自生磷灰石(AUT)七種類型。在前人研究的基礎上，利用窮舉法，研究不同磷灰石微量元素對物源判別的貢獻，尋找與評估磷灰石物源判別圖解的最佳端元，建立了新的磷灰石物源判別圖解。

圖1 不同巖石類型的磷灰石微量元素散點圖(底圖分區(qū)據Belousova et al., 2002)Fig.1 Scatterplots of published apatite trace element data of different rock types (base map after Belousova et al., 2002)

1 磷灰石微量元素數據預處理

1.1 數據收集

本研究收集了近三十年來，已經發(fā)表的20篇經典文獻，整理出了全球典型磷灰石微量元素數據集(圖2)。數據集中收集了不同母巖類型的610個磷灰石樣品，包括1925個測試點的微量元素數據，共計34628個元素值。其中測試點的微量元素數據包括：富堿性火成巖測試點數據263個(Zirneretal., 2015; Maoetal., 2016)，超鎂鐵質巖石測試點數據121個(O’Reilly and Griffin, 2000; Ihlenetal., 2014; Maoetal., 2016; Chakhmouradianetal., 2017)，鎂鐵質火成巖測試點數據843個(Sha and Chappell, 1999; Belousovaetal., 2002; Hsiehetal., 2008; Tollarietal., 2008; Chuetal., 2009; Tangetal., 2012; Ihlenetal., 2014; Maoetal., 2016)，長英質花崗巖測試點數據90個(Dill, 1994; Sha and Chappell, 1999; Hsiehetal., 2008; Chuetal., 2009)，中-低級變質巖測試點數據272個(Nishizawaetal., 2005; El Korhetal., 2009; Maoetal., 2016; Henrichsetal., 2018)，高級變質巖測試點數據172個(Beaetal., 1994; Bea and Montero, 1999; Belousovaetal., 2002; Nutman, 2007; Henrichsetal., 2018)和自生磷灰石測試點數據164個(Dill, 1994; Shields and Stille, 2001; Nishizawaetal., 2005; Joosuetal., 2015, 2016)。

圖3 磷灰石微量元素含量提琴圖Fig.3 Violin plot of apatite trace element data

論文數據不同的來源，不同的測試項目與檢測元素，導致每組數據的缺失值不同。綜合已整理的數據集，本文選用了各文獻中交集最多的14個元素及LREE(La+Ce+Pr+Nd)含量值進行最佳判別端元窮舉工作，服務于判別圖解的端元建立(見下文)。這14種元素包括La、Ce、Pr、Nd、Sr、Sm、Eu、Mn、Y、Gd、Dy、Yb、Th和U。該數據集La值含量變化從0.03×10-6到20488×10-6，平均值為1341×10-6；Ce值含量變化從0.22×10-6到53005×10-6，平均值為2484×10-6；Pr值含量變化從0.04×10-6到5121×10-6，平均值為259.66×10-6；Sr值含量變化從18.6×10-6到22498×10-6，平均值為1130×10-6；Nd值含量變化從0.37×10-6到18436×10-6，平均值為992.45×10-6；Sm值含量變化從0.12×10-6到2829×10-6，平均值為193.73×10-6；Eu值含量變化從0.04×10-6到278.4×10-6，平均值為23.1×10-6；Gd值含量變化從0.19×10-6到2031×10-6，平均值為187.93×10-6；Dy值含量變化從0.63×10-6到3232×10-6，平均值為158.34×10-6；Y值含量變化從10.02×10-6到19970×10-6，平均值為814.74×10-6；Yb值含量變化從0.58×10-6到1094×10-6，平均值為54.99×10-6；Mn值含量變化從3.50×10-6到26102×10-6，平均值為1097×10-6；Th值含量變化從0.0004×10-6到899.41×10-6，平均值為59.59×10-6；U值含量變化從0.001×10-6到1421×10-6，平均值為19.88×10-6；LREE值含量變化從0.83×10-6到97051×10-6，平均值為5077×10-6。微量元素含量提琴圖展示了該磷灰石微量元素數據集的分布與數據交集情況(圖3)。

1.2 數據清洗

數據清洗能夠最大程度地消除所構建數據集的統(tǒng)計學偏差。對所選取15個特征端元進行數據清洗，刪除特征元素為空值、負值和異常值的磷灰石數據。由于自生磷灰石類型缺少Sr、Th、U等多種元素的數據，因此所選取的研究對象為富堿性火成巖、超鎂鐵質巖石、鎂鐵質火成巖、長英質花崗巖、中-低級變質巖和高級變質巖共六種類型。數據集中各端元的數值大范圍分布(如，Ce：0.22×10-6～53004×10-6；U：0.001×10-6～1421×10-6)可通過對數處理使數據分布更加集中(Aitchison, 1982; Thomas and Aitchison, 2006; Wangetal., 2021)。清洗后磷灰石微量元素數據進行對數轉換，相較于原始數據，更利于數據分析和研究。

所建數據集中，鎂鐵質火成巖測試點數據843個，長英質花崗巖測試點數據90個，數據存在類間不平衡。本文故引入隨機抽樣法，最大限度平衡數據集中的數據。隨機抽樣法包括隨機過采樣法和隨機欠采樣法。隨機過采樣法通過復制少數樣本，并添加到原樣本中擴充原始數據集；隨機欠采樣法通過隨機選擇并刪除多數樣本的數據，以減小樣本(He and Garcia, 2009)。本文通過隨機過采樣的方法將長英質花崗巖數據完整復制2次、超鎂鐵質巖石數據完整復制1次和高級變質巖數據隨機選取一半復制1次；通過隨機欠采樣的方法處理富堿性火成巖數據、鎂鐵質火成巖數據和中-低級變質巖數據，以得到數據量相對平衡的數據集。

數據清洗后的磷灰石微量元素數據集符合大數據7V特征，即規(guī)模大(Volume)、多樣化(Variety)、動態(tài)性(Vitality)、準確性(Veracity)、價值化(Value)、高速性(Velocity)和可視化(Visualization)，更適用于大數據分析應用。該數據集包含磷灰石測試點的微量元素數據共1241個，其中富堿性火成巖測試點數據185個，超鎂鐵質巖石測試點數據186個，鎂鐵質火成巖測試點數據253個，長英質花崗巖測試點數據174個，中-低級變質巖測試點數據232個和高級變質巖測試點數據211個。此數據集為建立新判別圖解奠定了數據基礎。

2 構建磷灰石判別圖解

2.1 窮舉法篩選二元圖解端元

通過窮舉的方式，計算預處理后的磷灰石微量元素數據集中15個特征(La、Ce、Pr、Nd、Sr、Sm、Eu、Mn、Y、Gd、Dy、Yb、Th、U和LREE)任意兩特征間的比值，得到La/Ce、La/Pr、La/Nd等105種元素含量比值。聯(lián)合這15種元素含量值與所得105種元素含量比值，得到120個構建圖解的端元。再將這120個圖解端元中任意兩個端元進行組合，構建二元圖解，得到共7140個圖解。

圖4 輪廓系數排名前12的磷灰石母巖類型判別圖解Fig.4 Top 12 biplots of apatite in deifferent rock types from the dataset

圖5 磷灰石Eu/Y-Ce判別圖解的核密度估計圖Fig.5 Kernel density estimation of apatite Eu/Y vs. Ce discrimination diagram

圖6 磷灰石物源類型Eu/Y-Ce判別圖解Fig.6 Eu/Y vs. Ce discrimination diagram for apatite provenance classification

為評估上述7140個二元圖解的判別效果，本研究引入輪廓系數(Silhouette Coefficient)，以評估構建判別圖解的最適端元。輪廓系數是評價聚類效果好壞的一種量化方式(Rousseeuw, 1987)，由數據內聚度和分離度兩種因素所決定。內聚度評估同一分類簇內，數據分布的緊密程度。分離度評估不同分類簇之間，數據的分離程度。輪廓系數數學表達式為(1)：

(1)

其中a(i)為對象i與同一簇內所有其他對象的平均相異度(相異度：研究對象的差異程度)，b(i)為對象i與其最近的簇內所有對象的平均相異度。輪廓系數的范圍為[-1, 1]，當輪廓系數s(i)接近1時，“內部”差異a(i)遠小于“相鄰”差異b(i)，聚類效果良好；當s(i)接近-1時，“內部”差異a(i)遠大于“相鄰”差異b(i)，樣品i應被分類到“相鄰”簇內，聚類方式錯誤；當s(i)為0時，“內部”差異a(i)與“相鄰”差異b(i)相等，聚類無(明顯)關聯(lián)性(Rousseeuw, 1987)。

通過輪廓系數量化數據簇內部相異度與數據簇之間的相異度，對所有組合的圖解按照輪廓系數值降序排序。輪廓系數越大，表示該二元圖解判別磷灰石母巖類型的效果越好?；谳喞禂担疚暮Y選出效果最好的端元，來構建二元圖解。

2.2 核密度估計圈定數據分布

為了進一步精確和量化磷灰石物源判別圖解，在窮舉法得到的最優(yōu)端元圖解基礎上，本文運用核密度估計，圈定磷灰石微量元素數據分布。核密度估計(Kernel Density Estimation, KDE)，屬于非參數檢驗方法，由Rosenblatt (1956)和Parzen (1962)提出，又名Parzen-Rosenbaltt窗(Parzen-Rosenblatt window)，是概率論中估計位置密度的函數。核密度由核和帶寬決定，數學表達式為(2)：

(2)

其中K為非負函數，代表核；h為平滑參數，代表帶寬。核密度估計方法不對數據分布增加任何假定，而是根據數據本身的特點和性質進行分布擬合(Rosenblatt, 1956)?；诤嗣芏裙烙嫞筛鼫蚀_的圈定最優(yōu)判別圖解中不同磷灰石母巖類型的位置分布。

3 結果

通過評估輪廓系數，得到磷灰石物源判別效果最好的前12個圖解分別為Eu/Y-Ce圖解(SC=0.26929)、Eu/Y-LREE圖解(SC=0.26847)、Eu/Y-Pr圖解(SC=0.26826)、Eu/Y-La圖解(SC=0.26585)、Eu/Y-La/Sr圖解(SC=0.25959)、Eu/Dy-La/Sr圖解(SC=0.25036)、Sr/Y-La圖解(SC=0.25001)、Ce/Gd-La/Sr圖解(SC=0.24855)、Eu/Y-Pr/Eu圖解(SC=0.24772)、Eu/Y-Nd圖解(SC=0.24713)、Pr/Gd-La/Sr圖解(SC=0.24684)和Eu/Dy-Pr圖解(SC=0.24654)(圖4)。其中磷灰石物源判別效果最好的圖解為Eu/Y-Ce圖解，其輪廓系數為0.2629。為了更詳細的展示Eu/Y-Ce圖解中數據的分布情況，繪制了Eu/Y-Ce圖解的核密度估計圖(圖5)。

4 討論

4.1 數據處理

隨機抽樣方法的使用，能有效的提高判別圖解的準確率。當數據集存在類間不平衡時，分類結果更偏向于樣本容量大的類型，因此樣本容量小的類型分類結果不準確(Zhongetal., 2021)。He and Garcia (2009)提出，對于大多數不平衡的數據集，采樣技術的應用有助于提高分類器的準確性。但在數據預處理階段，過采樣與欠采樣的方法都可能會存在相應的問題。欠采樣方法從多數樣本隨機刪除對象時可能會導致多數類相關信息缺失，過采樣方法將少數樣本復制的數據添加到原始數據集中可能會導致過度擬合(He and Garcia, 2009)。本研究針對清洗后的鎂鐵質火成巖類型數據使用欠采樣方法時，按照磷灰石區(qū)域與類型，隨機刪除三分之二的數據，確保每種區(qū)域與類型均有信息保留。數據平衡前后的Eu/Y-Ce輪廓系數從0.1750增加到0.2693，磷灰石物源判別顯示出更好的結果，證明過采樣與欠采樣方法在本研究中的有效性。

4.2 篩選判別圖解端元

磷灰石物源判別圖解的端元窮舉結果顯示(圖4)，相較于單元素，元素比值具有更好的表現(xiàn)性。本研究中，單元素作為端元區(qū)分效果最優(yōu)的組合為Sr-La，但其輪廓系數0.2202，在所有端元組合中排名僅為68位，且各類型磷灰石均有嚴重重合現(xiàn)象，尤其是高級變質巖與長英質巖石和鎂鐵質火成巖幾乎完全重合(圖1)。加入元素比值后得到磷灰石物源判別圖解端元最優(yōu)組合為Eu/Y-Ce，輪廓系數為0.2693，Eu/Y-Ce圖解顯示出了更好的磷灰石物源判別效果(圖5)。

在窮舉法評估指標問題上，耿廳等(2019)在窮舉不同礦床類型鋯石圖解工作中引入了Calinski-Harabasz(CH)指標，CH指標亦為聚類效果評價常用的一種指標。本研究在進行此工作時發(fā)現(xiàn)CH指標在不同物源類型簇之間分離程度評判上表現(xiàn)較差，內部緊密但分離程度較差的結果會被給予過高分數，此結果不能充分利用有限的二維空間劃分數據分布。相比之下，輪廓系數在不同類型之間則顯示出更好的分離度。

4.3 核密度估計

核密度估計圖是觀測研究數據分布的有效方法，在數據可視化中，以一個或多個維度的連續(xù)概率密度曲線來表示數據分布。本研究將所得的最優(yōu)磷灰石物源判別圖解，即Eu/Y-Ce圖解，通過二維核密度估計，確定不同物源類型的數據分布，得到更加準確的磷灰石微量元素數據分布范圍。為了最大程度排除異常值干擾，本文選取了82%的密度進行磷灰石母巖類型劃分(圖6)。所得磷灰石物源判別圖解，相較于其他圖解更加全面，更加準確。

4.4 地球化學分析

隨著巖漿分異程度的增加，磷灰石微量成分Sr含量降低，Y含量增加，但是結果顯示，Sr和Y作為端元時對不同磷灰石母巖類型的區(qū)分效果較差，在Belousovaetal.(2002)的判別圖解(圖1)中Sr與Y所在的維度也有很大的重疊部分，說明Sr和Y作為單元素端元所攜帶的信息不能有效區(qū)分磷灰石母巖類型。斜長石是強烈富集Eu的礦物，隨著巖漿分異程度的增加，體系中斜長石的大量生成，Eu大量進入斜長石晶格中，導致磷灰石中Eu含量的降低(Dengetal., 2020a; Qiuetal., 2021)。巖漿分異過程中，Eu和Y這種負相關的變化表明Eu/Y能有效地區(qū)分不同磷灰石母巖類型。Belousovaetal.(2002)認為Y和Eu的比值是有用的磷灰石判別指標。Ce與Eu均為變價稀土元素，具有不同的化合價態(tài)：Ce3+/Ce4+和Eu2+/Eu3+，可以反映形成體系的氧逸度變化(Ballardetal., 2002; Watsonetal., 2006; 張紅等, 2018; 邢凱等, 2018; Dengetal., 2018)。Ce3+和Eu3+作為親磷灰石元素價態(tài)，更容易替代磷灰石中的Ca2+進入晶格中，磷灰石顯示出更高的Ce和Eu含量。Eu/Y-Ce判別圖(圖6)中可以看出，相比于高級變質巖，中-低級變質巖的形成處于偏氧化的環(huán)境，Ce3+轉變?yōu)镃e4+，從磷灰石中釋放出來，顯示出偏低的Ce含量；相比于長英質花崗巖，鎂鐵質火成巖與超鎂鐵質巖石的形成處于偏還原的環(huán)境，Ce4+轉變?yōu)镃e3+進入磷灰石晶格中，磷灰石顯示出更高的Ce含量。邢凱等(2018)認為磷灰石中Ce和Eu這兩個具有相反分配行為的元素，對說明母巖漿的氧化還原環(huán)境具有重要意義。所以本文構建的Eu/Y-Ce圖解能夠顯示最優(yōu)的磷灰石物源判別效果。

5 結論

(1)磷灰石Eu/Y-Ce判別圖解是通過窮舉所有特征元素組合，引入輪廓系數篩選最優(yōu)端元組合，輔之以核密度估計方法，選取82%的密度劃分不同物源區(qū)域所得到的最優(yōu)圖解。相比于其他判別圖解，本文提出的磷灰石Eu/Y-Ce判別圖解更加準確，可涵蓋更加全面的巖石類型，是有效的磷灰石物源判別圖解。

(2)本次研究結果表明，影響磷灰石物源判別效果的主要因素可能包括磷灰石形成時的氧化還原狀態(tài)和微量元素在共生礦物組合變化過程中的不同配分行為。氧化環(huán)境中，Eu2+轉變?yōu)镋u3+，進入磷灰石晶格中，顯示出更高的Eu含量；Ce3+轉變?yōu)镃e4+，從磷灰石中釋放出來，顯示出更低的Ce含量。同時，斜長石的生成導致大量Eu進入斜長石，磷灰石也會顯示出更低的Eu含量。

(3)本文通過大數據技術，分析指示礦物的性質來識別礦物物質來源。本研究對磷灰石微量元素的數據挖掘工作，是將大數據技術運用在地球科學研究中初步探索。隨著未來磷灰石地球化學數據的更加豐富，結合更多算法，高緯度元素判別圖解的建立值得進一步探索。

致謝論文的完成得益于鄧軍院士的指導。感謝俞良軍老師對本文細心的審閱；感謝兩位匿名審稿人對本文提出了詳細的建設性意見；感謝中國地質大學(北京)的李珊珊博士后、龍政宇博士、朱紫怡和周飛為本文提供寶貴的修改意見。