方思源，董少春，胡 歡

南京大學 地球科學與工程學院，南京 210023

1 前言

近年來，基礎地質(zhì)數(shù)據(jù)表現(xiàn)出了爆炸性增長的態(tài)勢，復雜且大量的地質(zhì)數(shù)據(jù)開始被系統(tǒng)性的積累整合，形成了具有相當規(guī)模的大數(shù)據(jù)集，地球科學從一個數(shù)據(jù)匱乏的時代逐步邁入大數(shù)據(jù)時代（Guo et al.，2017；吳沖龍等，2018；董少春等，2019）。地球科學大數(shù)據(jù)巨大的信息量改變了地球科學研究的方法、思路和理論，提供了多樣化的數(shù)據(jù)挖掘手段和知識發(fā)現(xiàn)途徑（翟明國等，2018；Karpatne et al.，2019）。地球科學大數(shù)據(jù)分析更加注重挖掘和分析數(shù)據(jù)之間隱藏的內(nèi)在聯(lián)系與特征，以數(shù)據(jù)驅(qū)動來發(fā)現(xiàn)規(guī)律，解釋現(xiàn)象，探求問題，減少了主觀經(jīng)驗上的人為干擾（張旗等，2015；張旗和周永章，2017；張旗等，2018；周永章等，2018；羅建民和張旗，2019）。地球科學大數(shù)據(jù)的處理也給地學研究帶來了新的挑戰(zhàn)。大數(shù)據(jù)的處理通常離不開超級計算機、云計算、并行計算等計算機技術的支撐（Wyborn and Evans，2015；廖湘科等，2016）。但是由于缺乏適合于地球科學大數(shù)據(jù)的專門算法、模型和軟件，因此難以對大規(guī)模數(shù)據(jù)集進行快速、有效的處理，一定程度上限制了地球科學大數(shù)據(jù)的應用發(fā)展。因此，根據(jù)大數(shù)據(jù)的特點和數(shù)據(jù)處理需求開發(fā)適合地球科學大數(shù)據(jù)的處理、分析和可視化算法、模型或軟件勢在必行。

以巖石地球化學研究為例，元素協(xié)變圖（Element Variation Diagram）常被用來研究元素之間的相關關系，對巖漿演化和不同構造背景下元素的地球化學行為進行分析。元素協(xié)變圖清晰、直觀地展現(xiàn)了元素在巖漿演化過程中的變化趨勢，結(jié)合地質(zhì)構造背景等信息可以幫助分析不同地質(zhì)過程對元素地球化學行為的影響以及成礦元素的時空分布規(guī)律性等問題，為揭示成礦元素的富集成礦規(guī)律，指出找礦的區(qū)域構造方向等提供依據(jù)（王登紅等，2015；蒲東，2019；羅建民等，2019；周永章等，2019）。

元素協(xié)變圖通?；跀?shù)據(jù)值在全值域采用單一的線性或非線性函數(shù)模型進行擬合，生成趨勢線來揭示元素的地球化學行為特征。但由于采用單一函數(shù)模型擬合整個巖漿演化過程中元素的變化趨勢，忽略了巖漿演化不同階段元素地球化學行為特征的差異性。當數(shù)據(jù)集較小時，這種差異性可能表現(xiàn)得并不明顯；而大數(shù)據(jù)支撐下的元素地球化學行為特征在不同階段的差異性可能非常明顯。如果仍然采用單一的函數(shù)進行擬合，勢必掩蓋掉這種差異性，從而忽略了重要的演化特征，難以真實反映演化趨勢。因此分段擬合趨勢線的思路更符合大數(shù)據(jù)支撐下元素協(xié)變圖的繪制和分析。

本文以元素協(xié)變圖為例，從地球科學大數(shù)據(jù)的角度出發(fā)，針對元素在巖漿演化中的變化特點，突破傳統(tǒng)函數(shù)理論模型的限制，采用積分學中“窮竭法”的思路，設計了一個針對元素行為趨勢研究的多區(qū)間線性回歸擬合模型，并基于此模型開發(fā)出一款基于巖石地球化學大數(shù)據(jù)分析元素在巖漿演化過程中地球化學行為的可視化軟件。利用該軟件，我們對GEOROC數(shù)據(jù)庫中的地球化學元素進行數(shù)據(jù)處理，基于大規(guī)模地球化學數(shù)據(jù)集生成不同構造背景下的元素協(xié)變圖，以進行地球化學行為的分析。

2 方法原理

2.1 元素協(xié)變圖

元素協(xié)變圖通過數(shù)字化的方式將樣品數(shù)據(jù)中元素含量（或比值）設為橫縱坐標，繪制成二元或三元平面散點圖。其中，哈克圖解（Harker diagram）（Harker，1909）是最有代表性的元素協(xié)變圖之一。因為隨著巖漿演化的進行，通常SiO2含量升高，所以常將SiO2含量設為元素協(xié)變圖的橫坐標以指示巖漿演化的方向，并將其他主要元素的氧化物含量設為縱坐標，由此分析巖漿演化過程中主要元素含量的變化趨勢以及判斷不同礦物分離結(jié)晶對元素變化的影響。隨著研究的不斷深入，元素協(xié)變圖不再局限于分析主要元素的氧化物，而是擴展到了微量元素含量、元素比、同位素比值等地球化學參數(shù)。通過元素協(xié)變圖分析，一些元素（組合）、元素之間的相互關系等趨勢特征可以作為“指紋元素”，在判別各類巖石成因、構造背景、演化階段上具有良好的指示作用（Green，1995；李永軍等，2015）。

2.2 窮竭法的區(qū)間劃分

古希臘數(shù)學家歐多克斯最先創(chuàng)立了“窮竭法”（Method of Exhaustion）的思想，后來由歐幾里得將其成果收入到《幾何原本》中，阿基米德對其進行了進一步完善，最后被引用于積分思想中（邵明湖，1990）。近現(xiàn)代“窮竭法”在經(jīng)過不斷的完善改良后，經(jīng)常被用于求取平面坐標中一些由曲線圍成的面積。其基本思路是：通過沿X軸做切割，將X軸平均劃分成n個區(qū)間段，在區(qū)間內(nèi)近似的認為曲線是水平不變的，因此可以將其看成n個矩形。每一區(qū)間段中點對應的函數(shù)值代表了矩形的高，依次計算每個矩形的面積S1+S2+…+Sn，則這些矩形的面積之和SN可以較為精確的代表曲線圍成的面積S（圖1）。利用窮竭法求取面積S的公式如公式1所示。

在分析巖漿演化過程中元素行為特征時，通常用橫坐標的元素指示巖漿演化的方向。區(qū)別于以往在整個橫坐標采用單一函數(shù)擬合生成元素協(xié)變圖的方法，本文利用窮竭法的思想將橫坐標平均劃分為若干個區(qū)間，每個區(qū)間代表巖漿演化過程中的一段很短的過程。在每個很短的時間段內(nèi)近似認為元素含量特征基本相同或者滿足同一簡單的函數(shù)模型，因此可以取一固定值或函數(shù)值代表每個區(qū)間的元素特征。

圖1 窮竭法區(qū)間劃分示意圖Fig. 1 Schematic diagram of interval division in the method of exhaustion

2.3 區(qū)間線性擬合

區(qū)間內(nèi)擬合模型的選擇對最終趨勢線擬合效果有著直接影響。擬合模型通常需要結(jié)合區(qū)間內(nèi)的樣品數(shù)量、數(shù)據(jù)質(zhì)量、分布方式等因素綜合考慮。劃分后每個區(qū)間的跨度很小，因此對應每個較短的演化過程中的元素變化并不復雜。在保證了區(qū)間內(nèi)數(shù)據(jù)的數(shù)量與質(zhì)量的前提下，可以認為區(qū)間段內(nèi)元素的含量呈簡單的線性變化，然后在每個區(qū)間內(nèi)進行線性擬合，獲得該區(qū)間內(nèi)的回歸趨勢線。本文采用最小二乘法對區(qū)間內(nèi)數(shù)據(jù)點進行線性擬合，對每個區(qū)間段內(nèi)的一系列樣本點（y1，y1）、（x2，x2）……（xn，yn）構建一元線性回歸方程。

獲得區(qū)間內(nèi)線性回歸線后，在每個區(qū)間的回歸線上取中點作為特征點，代表該區(qū)間對應演化階段中元素的特征含量。將每個區(qū)間的特征點連接起來即可獲得一條代表元素含量隨巖漿演化的整體變化趨勢線。

3 方法實現(xiàn)及可視化應用

區(qū)間線性擬合方法的實現(xiàn)主要包含以下步驟（詳見圖2）：

基于上述思路，我們在Matlab GUI框架下實現(xiàn)了該模型，并開發(fā)了可視化數(shù)據(jù)操作和輸出界面，設計出一個面向大數(shù)據(jù)的元素地球化學行為分析應用 軟 件（Elemental Behavior of Magmatic Evolution，EBME）。EBME軟件界面如圖3所示。

EBME以巖石地球化學數(shù)據(jù)庫為基礎分析研究巖漿演化過程中重要的地球化學元素在不同構造背景下含量的變化趨勢，具有繪制擬合線、保存圖片、保存擬合線特征數(shù)據(jù)等多項功能。EBME支持動態(tài)加載大數(shù)據(jù)集和簡單參數(shù)設置（如選擇橫縱坐標元素、設定區(qū)間長度、選擇制定條件的數(shù)據(jù)等）。用戶使用該軟件不需要具備Matlab編程基礎，也無需安裝Matlab軟件，通過可視化操作即可非常方便地完成區(qū)間數(shù)據(jù)擬合與繪圖過程，實現(xiàn)基于大數(shù)據(jù)的元素協(xié)變圖的數(shù)據(jù)處理與分析功能。EBME軟件主要操作流程為：（1）選擇并導入巖石地球化學數(shù)據(jù)庫文件（支持Excel格式批量數(shù)據(jù)導入）；（2）選擇一種或多種區(qū)域構造或大地構造背景對數(shù)據(jù)進行篩選（默認加載全部數(shù)據(jù)集）；（3）選擇橫、縱坐標元素，設置x軸區(qū)間、劃分擬合區(qū)間、標準差數(shù)值；（4）點擊繪圖進行數(shù)據(jù)投圖與趨勢擬合，并在右側(cè)繪圖區(qū)顯示。（5）點擊“保存圖片”按鈕可導出結(jié)果圖片，點擊“導出趨勢數(shù)據(jù)”按鈕可導出趨勢擬合線特征點的坐標，點擊“清空”按鈕即可刪除軟件加載的所有數(shù)據(jù)。

圖2 分區(qū)間擬合方法的流程示意圖Fig. 2 Flow chart of dividing interval fitting method

圖3 EBME界面示意圖Fig. 3 Interface of EBME

4 軟件應用

4.1 數(shù)據(jù)來源

本文數(shù)據(jù)來源于GEOROC（Max Planck Institute for Chemistry Mainz, Germany, 2018-11-17）數(shù)據(jù)庫。該數(shù)據(jù)庫是馬克斯·普朗克化學研究所建設與維護的海洋和大陸巖石數(shù)據(jù)庫，共發(fā)布了從1883年至今公開發(fā)表的數(shù)百萬個包含100多種元素或同位素含量的巖石地球化學數(shù)據(jù)。GEOROC數(shù)據(jù)庫建設成熟，數(shù)據(jù)量大，開放性好，有堅實的理論基礎，已取得了豐富的研究成果，在巖石、火山、沉積、古海洋以及大氣等科學研究中應用廣泛（葛粲等，2018；焦守濤等，2018； Liu et al.，2018；Ueki et al.，2018；Zhao et al.，2019）。

4.2 數(shù)據(jù)的整合和清洗

GEOROC數(shù)據(jù)庫含有400多個數(shù)據(jù)表。由于每個數(shù)據(jù)表的屬性名稱、內(nèi)容格式各異，直接下載的數(shù)據(jù)無法直接利用，主要存在以下三點問題：（1）不同數(shù)據(jù)表中的表格屬性名稱（字段名）、屬性個數(shù)、排列順序不一致。（2）不同數(shù)據(jù)表中數(shù)據(jù)的單位量綱不一致，或無標注，部分數(shù)據(jù)缺少屬性標注，變量關系不清晰。（3）部分屬性名表述不清，數(shù)據(jù)格式錯誤。

本文對數(shù)據(jù)庫中火成巖全巖地球化學數(shù)據(jù)進行了如下整合和清洗，以便開展與巖漿演化有關的元素地球化學行為分析：（1）對各個表格中的屬性進行分析和比較，根據(jù)應用需求對相同含義的屬性進行整合，統(tǒng)一屬性名稱與順序。例如：Fe元素的百分含量統(tǒng)一換算成FeOT；（2）對單位量綱不一致的數(shù)據(jù)進行單位換算，統(tǒng)一量綱，也對屬性格式表述方式進行統(tǒng)一調(diào)整。例如，主量元素和燒失量（LOI）的單位統(tǒng)一為wt%，微量元素的單位統(tǒng)一為（×10-6）；（3）剔除內(nèi)容錯誤、模糊、定義不清或存在異常的數(shù)據(jù)。此外，為便于分析不同構造背景下元素地球化學行為的演化特征，我們將所有數(shù)據(jù)按照Intra Continental（陸內(nèi)構造背景）、Convergent Margins（匯聚板塊邊緣構造背景）和Oceanic（海洋構造背景）進行分類。經(jīng)過上述數(shù)據(jù)整理和清洗，共計獲得巖石地球化學數(shù)據(jù)樣本261350條，其中包括111237條陸內(nèi)數(shù)據(jù)，122044條匯聚板塊邊緣數(shù)據(jù)和28069條海洋數(shù)據(jù)。每條數(shù)據(jù)均記錄有樣品采集的地理位置（緯度和經(jīng)度坐標）、巖石類別、樣品測試類型和元素含量等信息。

4.3 結(jié)果分析

4.3.1 參數(shù)設置

批量加載從GEOROC數(shù)據(jù)庫采集到的261350條樣品數(shù)據(jù)作為數(shù)據(jù)源，將MgO的百分含量作為橫坐標軸，含量限制在0～30 wt%范圍內(nèi)。MgO含量由高至低即x坐標軸從右往左代表巖漿演化的進行方向，并將MgO含量劃分成15個區(qū)間，即每段小區(qū)間中MgO含量的跨度為2 wt%，那么元素趨勢特征點橫坐標MgO含量分別為1、3、5、7……29（wt%）。由EBME根據(jù)區(qū)間內(nèi)線性趨勢線擬合方法分析不同元素在陸內(nèi)、海洋、匯聚板塊三種大地構造背景下隨巖漿演化的含量變化行為。

4.3.2 趨勢線生成及分析

根據(jù)上述參數(shù)設置，采用區(qū)間線性擬合模型生成了高場強元素、大離子親石元素等元素與MgO的元素協(xié)變圖，匯總?cè)鐖D4所示。從圖4中可以通過演化趨勢線看出不同元素在巖漿演化過程中的富集與虧損趨勢、元素之間的相關關系以及元素在不同大地構造背景下的差異性變化。

4.3.2.1 高場強元素及稀土元素

根據(jù)圖4所示，隨著巖漿演化的進行，匯聚板塊邊緣和海洋兩個環(huán)境下高場強元素總體表現(xiàn)出隨巖漿演化逐漸富集的趨勢。在不同構造背景下，部分元素表現(xiàn)出一定的差異性：匯聚板塊構造背景下的高場強元素的含量低于其他兩種構造背景，海洋構造背景下的元素富集程度更高。其中圖4（1）-（20）顯示Nb、Ta、Th以及輕稀土元素LREE在巖漿演化前期MgO含量大于20 wt%階段，其元素在陸內(nèi)構造背景下的含量明顯為高值，并且有一定的虧損趨勢，而在其他兩種構造背景下則呈現(xiàn)低值并有微弱的富集趨勢，直到演化后期MgO含量降低到8%左右時，Nb、Ta、Th和LREE元素才呈現(xiàn)出較為一致的富集趨勢。然而Zr、Hf和HREE卻沒有出現(xiàn)這樣的特征，不同構造背景下整體富集趨勢較為統(tǒng)一。

4.3.2.2 大離子親石元素

圖4（21）～（24）顯示大離子親石元素（Ba、Na2O、K2O、Rb）隨著巖漿演化總體也表現(xiàn)出逐漸富集的趨勢。其中Ba元素同樣也表現(xiàn)出在陸內(nèi)構造背景下，演化初期富集程度很高，在其他兩種背景下則明顯偏低的特征。此外，Na2O、K2O、Rb與富集趨勢存在差異：K2O、Rb元素在演化初期富集趨勢微弱較為平緩，但是到MgO含量為8%左右時，富集趨勢徒增，而Na2O演化趨勢基本呈線性富集，演化后期海洋構造背景下其富集程度相對高于另兩種構造背景。

4.3.2.3 其他

圖4（25）顯示Zn元素在巖漿演化過程中具有微弱的富集特征，但在不同的構造背景下，巖漿演化后期表現(xiàn)出顯著差異性趨勢。海洋構造背景下Zn含量明顯高于其他兩類構造背景，而在匯聚板塊構造背景下Zn含量最低。在巖漿演化后期，在海洋構造背景下Zn含量呈現(xiàn)上升趨勢，而在匯聚板塊構造背景下，則呈現(xiàn)出下降的趨勢。

圖4 基于區(qū)間線性擬合方法繪制的元素協(xié)變圖Fig. 4 Element variation diagram based on the dividing interval linear fitting method

由此可見，通過不同構造背景下分區(qū)間擬合的巖石地球化學大數(shù)據(jù)元素協(xié)變圖分析，揭示出不同元素在巖漿演化過程中表現(xiàn)出差異性的地球化學行為，較傳統(tǒng)認識更能體現(xiàn)出演化趨勢和特征，為更好的認識巖漿演化過程中的元素地球化學行為提供了更多證據(jù)。

5 總結(jié)

本文針對以往元素協(xié)變圖生成方式無法適應地球化學大數(shù)據(jù)集的弊端，引入窮歇法的思路，采用分區(qū)間線性擬合方法繪制巖漿演化過程中元素變化的趨勢圖，以直觀、清晰的方式揭示巖漿演化不同階段元素富集與虧損特征，并在Matlab GUI平臺上開發(fā)了面向巖石地球化學大數(shù)據(jù)的元素協(xié)變圖自動繪制軟件——EBME。該軟件具有操作界面友好，參數(shù)設置簡單，運行速度快等特點，大大降低了大數(shù)據(jù)處理分析的門檻。通過從GEOROC數(shù)據(jù)庫中收集的261350條火成巖樣品數(shù)據(jù)進行整合分析，利用該軟件自動繪制生成了新型元素協(xié)變圖。結(jié)果表明其擬合的趨勢線很好的揭示了元素在不同構造背景和不同演化階段下的演化趨勢，為深度挖掘和分析巖漿演化過程中元素地球化學行為提供了新的工具和思路。但現(xiàn)階段EBME軟件的數(shù)據(jù)處理模塊還比較簡單，僅限于基于窮歇法的區(qū)間線性擬合模型，尚不具備數(shù)據(jù)預處理功能，數(shù)據(jù)處理模型和算法還有待進一步豐富和完善。

地球科學大數(shù)據(jù)的快速發(fā)展，為面向大數(shù)據(jù)處理需求的專業(yè)軟件的研制與開發(fā)帶來了新的機遇和挑戰(zhàn)。根據(jù)大數(shù)據(jù)特點和處理需求精準定制的分析處理軟件必將在大數(shù)據(jù)驅(qū)動下的地球科學研究中發(fā)揮更大的作用。

致謝：感謝馬克斯·普朗克化學研究所提供了GEOROC數(shù)據(jù)集，感謝中國科學院地球化學研究所高劍峰研究員對程序設計和本文寫作提供的寶貴意見，感謝南京大學地球科學與工程學院的周會群教授對論文算法和寫作給予的修改意見。感謝兩位審稿人和編輯對本文提出的寶貴意見和建議。文章系“深時數(shù)學地球”（Deep-time Digital Earth）大科學計劃系列成果之一。