地壳厚度的时空变化不仅直接反映了板块边界的造山过程，也反映板内稳定克拉通减薄裂解破坏等重要的地球动力学过程(Yin and Harrison, 2000; 张旗等, 2001; 嵇少丞等, 2008; 葛粲等, 2011; 张旗, 2011; 朱日祥等, 2012)，而且地壳增厚和减薄引起的次级效应会对地球表面环境和气候改变、生物演化、矿产资源的富集来深远的影响(Chiaradia, 2014; Klein et al., 2014; Vérard et al., 2015; 刘少峰和王成善, 2016; 刘欣雨等, 2017)。因此，数字重建地壳厚度的时空变化对于研究地球系统的发展演变具有重要意义。地球物理方法能为我们提供现今地壳形态，但却无法回溯历史地壳厚度变化情况。为此，我们需要借助一些记录古代地壳厚度信息的“化石”来推测其变化过程。通过近50年的研究实践，研究者根据岛弧地区岩浆岩有限的数据集提出了几种与地壳厚度有一定相关性的地球化学指标，具有相关性的化学元素主要包括K₂O、Na₂O和CaO等主量元素以及Sr、Y、La、Ce、Sm、Dy、Yb等微量元素，提出了K60、K57.5、K55、Na6.0、Ga6.0、Ce/Y、Sm/Yb、Dy/Yb、Sr/Y、La/Yb等地球化学指标来推测古地壳的时空变化信息(Condie and Potts, 1969; Jakeš and White, 1972; Condie, 1973; Leeman, 1983; Hildreth and Moorbath, 1988; Plank and Langmuir, 1988; Haschke et al., 2002; Haschke and Günther, 2003; Lee et al., 2007; 汪洋, 2007; Mantle and Collins, 2008; Mamani et al., 2010; Chapman et al., 2015; Chiaradia, 2015; Profeta et al., 2015)。前人根据数据集统计拟合了这些指标与地壳厚度之间的对应关系，并将其外推运用于古造山带地区，根据不同年代的地化指标变化约束古地壳厚度随时间的变化情况，研究成果涉及北美加拿大前寒武地盾(Condie and Potts, 1969)、北美科迪勒拉山造山带西南部盆岭省(Chapman et al., 2015)、南美中安第斯造山带(Mamani et al., 2010)、新西兰造山带(Mantle and Collins, 2008)、燕山造山带燕辽地区(邓晋福等, 1996)、太行山造山带(罗照华等, 1997)、三江地区义敦岛弧、东南沿海岩浆弧、藏南冈底斯弧、燕山造山带辽西地区火山岩(官文慧和汪洋, 2017)。然而，前人研究所采用地化指标的含义、数据集来源、数据范围、地壳厚度约束范围、数据筛选条件、拟合公式等均不相同，更为重要的是，其使用数据集的数据量较为有限，样品地点往往局限在岛弧地区、带有区域性特色，岩浆来源相对单一，源区基本为地幔楔。本文基于全球火山岩地球化学数据集，研究全球尺度复杂构造背景下，火山岩的主量元素、微量元素比值与地壳厚度之间的统计关系；了解火山岩不同地化指标相对于现今地壳厚度的数据分布特征；分析壳源和幔源两类不同类型的火山岩在约束地壳厚度方面的差异；探讨地化指标与地壳厚度具有相关性的可能原因。

德国美因茨马克斯普朗克化学研究所建立并维护了GEOROC (Geochemistry of Rocks of the Oceans and Continents)数据库(http://GEOROC.mpch-mainz.gwdg.de/GEOROC/Start.asp)。该数据库综合收集汇总了15380篇发表刊物、来自11个不同的地质背景、46.3万个火山岩石和地幔捕虏体岩石样品、107.5万个分析结果的分析资料。数据中包括地理位置、纬度和经度、岩石类别和岩石类型、分析方法、实验室、参考资料和引用来源、主量和微量元素的浓度、放射性和非放射性同位素比值，以及整个岩石、玻璃、矿物和包含物的分析年龄。本文从GEOROC数据库中提取了11个不同地质背景火山岩的主量元素、微量元素、样品年龄、样品经纬度信息和5种不同矿物的微量元素数据，作为分析数据。

本文将数据集分为年轻壳源、古老壳源、年轻幔源和古老壳源4个子集，分别统计不同子集下数据随地壳厚度变化的分布规律。壳源岩浆与幔源岩浆在主量元素含量、同位素组成、微量元素和稀土元素组成上都有差别，但是尚未有绝对可信的判别准则。一般认为基性地壳发生脱水熔融产生的熔体一般具有高硅含量(Sen and Dunn, 1994; Qian and Hermann, 2013)，作为加厚下地壳熔体的典型代表，大别低镁埃达克岩的SiO₂平均值也接近70 %(He et al., 2011)。研究者普遍认为壳源岩浆以一般中、酸性为主，幔源岩浆形成的火山岩以超基性、基性为主，但并非所有中性火山岩都来自壳源，幔源岩浆也可以形成安山岩。本文粗略地以SiO₂含量57%为界，将数据集分为壳源部分和幔源部分，分别进行统计。由于本文使用的地壳厚度是现今地壳厚度，在远古时期地壳厚度相对于现今可能经过了剧烈的变化，当时生成的火山岩可能与当时的地壳厚度相关，但与现今地壳厚度无关。为了尽可能的排除这些数据的影响，同时保留足够多的数据进行数据统计，本文规定了一个时间分界，以岩浆岩年龄23Ma(中新世)为界，将数据分成年轻火山岩和古老火山岩。只分析年轻火山岩地化指标与现今地壳厚度的关系，古老火山岩的地化指标作为参考。

crust 1.0模型(Laske et al., 2013)是一个全球的三维地壳模型，是由crust 2.0模型(Bassin et al., 2000)和crust 5.1模型(Mooney et al., 1998)升级而来。模型将全球划分为1度×1度的网格，共计64800个单元，每个单元都提供了包括冰层、水层、上、中、下沉积层、上、中、下地壳、上地幔的波速和密度分层模型。地表高程、海底高程和冰层来自ETOPO1模型(Amante and Eakins, 2009)，莫霍埋深是利用人工地震研究、远震接收函数研究和重力约束共同确定。本文使用的地壳厚度模型是根据crust 1.0模型排除水层和冰层后的固体地壳厚度模型计算获得。岩石样品所在地的现今的地壳厚度，是根据样品采样地点的经纬度，通过全球地壳厚度模型插值获得，通过这种方法将地球化学数据与地球物理模型联系起来。

由于地球化学元素含量同时受到多种因素的影响，单个样本的元素含量值看似随机，但是大量样本情况下有明显的统计规律，直方图是一种直观的展现数据样本分布规律的方法，但是有一定缺陷，直方图中统计分组界限的选择对结果有很大影响，并且展现的分布曲线并不平滑，不利于进行统计参数估计。由于真实数据分布与正态分布有很大差距，利用全体数据的均值和标准差分析数据分布的集中趋势和离散程度有很大偏差(图 1)。为了分析不同地壳厚度区域中火山岩地球化学指标的分布规律，本文引入了核密度估计的方法。核密度估计方法是一种从数据样本本身出发研究数据分布特征的方法，在概率论中用来估计未知的密度函数，在统计学理论和应用领域均受到高度的重视(Danese et al., 2008)。本文利用正态分布函数作为核函数，利用核密度估计方法统计每个地壳厚度区间的地球化学概率密度函数，并通过概率密度函数主峰位置和主峰宽度，来描述数据分布的集中趋势和离散性(图 1a)。

图 1

Fig. 1

图 1 真实数据分布情况与统计方法对比 (a)直方图所统计的数据集为地壳厚度25~30km之间、SiO2含量小于57%、年龄小于23Ma的火山岩K2O含量统计样本集；红色曲线是利用正态分布函数描述该样本集，红色误差棒表示数据集的均值和标准差；蓝色曲线是利用核函数估计的概率密度分布函数，蓝色误差棒表示该密度分布函数的最高峰位置和峰宽; (b)归一化概率分布色谱图，利用颜色代表概率密度的相对高低，误差棒标志与(a)相同 Fig. 1 Comparison of real data distribution with statistical methods (a) histogram data sets are statistical K2O contents of volcanic rock samples with a crustal thickness of 25~30km, SiO2 content less than 57% and age less than 23Ma; red curve is described by normal distribution function, red error bar represents the mean and standard deviation of the data set; blue curve is the probability density function estimated using kernel density estimation method, the blue error bar represents the peak position and the peak width of the density function; (b) normalized probability distribution map, using color to represent the relative probability density, the error bar marker is the same as the left side

为了研究全球火山岩地球化学指标与地壳厚度的关系，本文将地壳厚度以5km为分组宽度，从将地壳厚度从15km到75km，每隔2.5km为分组中心来确定分组区间，分别统计全体数据集、年轻-壳源子集、年轻-幔源子集、古老-壳源子集和古老幔源子集中落入分组组区间数据概率密度函数，由于需要观察概率密度函数随地壳厚度的变化情况，因此将核函数估计的概率密度函数转化成以色彩显示的概率密度分布图(图 1b)。将不同地壳分组区间的概率密度函数归一化后合并在一起，形成地化指标随地壳厚度变化归一化联合概率密度分布图，简称联合分布图(图 2和图 3)。通过联合分布图可以直观的透视数据的分布情况，有利于掌握数据宏观分布状态，总结古、今、壳源、幔源火山岩地化元素含量的分布规律。全体数据集的地壳厚度-SiO₂含量的联合分布图显示(图 2a)，不同厚度区域的火山岩SiO₂数据分布有不同的集中趋势。在地壳厚度15km到52.5km之间，火山岩样品主要聚集在基性玄武岩区域，地壳厚度55km以上地区的火山岩样品SiO₂含量明显高于普通地壳区域的含量，但数据较为离散。对比SiO₂和MgO的联合分布图可以看到，壳源联合分布图以高SiO₂低MgO为特点(图 2b, c, g, h)，幔源的联合分布图以低SiO₂高MgO为特点(图 2d, e, i, j)。以SiO₂含量57%为界，划分壳源子集的MgO含量普遍较低，数据分布较为集中，一般分布在4%以下(图 2g, h)，而幔源子集的MgO含量普遍偏高，数据分布分散，一般在4%以上(图 2i, j)。

图 2

Fig. 2

图 2 不同数据集主量元素含量随地壳厚度变化的归一化联合概率密度分布图 (a-y)显示了地壳厚度-主量元素指标归一化联合概率密度分布情况，色标与图 1(b)一致.每行代表不同主量元素的统计结果，从上到下分别是SiO2、MgO、CaO、Na2O和K2O；每列代表不同数据集统计结果，从左到右分别是全体数据集、年轻-壳源数据集、古老-壳源数据集、年轻-幔源数据集和古老-幔源数据集 Fig. 2 Normalized joint probability density distributions of major element contents with the variation of crustal thickness in different data sets (a-y) shows the normalized joint probability density distribution of crustal thickness-principal elements. The color scale is consistent with Fig. 1b. Each row represents the statistical results of different principal elements, from top to bottom are SiO2, MgO, CaO, Na2O and K2O; each column represents the statistical results of different data sets, from left to right are the whole data sets, young-crust data sets, old-crust data sets, young-mantle data sets and old-mantle data sets, respectively

图 3

Fig. 3

图 3 不同数据集的微量元素比随地壳厚度变化归一化联合概率密度分布图 (a-y)显示了地壳厚度-地化指标归一化联合概率密度分布情况，色标与图 1b一致.每行代表不同微量元素比的统计结果，从上到下分别是Sr/Y、Ce/Y、La/Yb、Dy/Yb和Sm/Yb；每列代表不同数据集统计结果，从左到右分别是全体数据集、年轻-壳源数据集、古老-壳源数据集、年轻-幔源数据集和古老-幔源数据集 Fig. 3 Normalized joint probability density distributions of trace element ratios with the variation of crustal thickness in different data sets (a-y) shows the normalized joint probability density distribution of crustal thickness and geochemical index. The color scale is consistent with Fig. 1b. Each row represents the statistical results of different trace element ratios, from top to bottom are Sr/Y, Ce/Y, La/Yb, Dy/Yb and Sm/Yb; each column represents the statistical results of different data sets, from left to right are the whole data sets, young-crust data sets, old-crust data sets, young-mantle data sets and old-mantle data sets, respectively

地球化学和地球物理学研究发现，垂直于相同年龄的弧-沟体系走向方向上的岩浆岩的岩性从拉斑玄武岩系列逐渐转变为钙碱性亚系列和钾玄岩亚系列，在相同SiO₂含量条件下，其K₂O含量等随着距离海沟的远近而显示出系统的变化特征，这种岩浆岩化学极性被认为是俯冲带深度增加的函数(Moore, 1962; Dickinson and Hatherton, 1967; Hatherton and Dickinson, 1968, 1969)。由此将K₂O含量与地壳厚度定量化地联系起来(Condie and Potts, 1969; Condie, 1973, 1976; Leeman, 1983；汪洋, 2007)，Condie (1973)利用SiO₂归算到60%的K60指标表现出与地壳厚度呈现线性递增关系，汪洋(2007)研究表明K60指标与地壳厚度呈现的对数函数关系。Plank and Langmuir (1988)统计了25个岩浆弧的100座活火山玄武岩分析数据，发现MgO含量为6%时的Na₂O含量(Na6.0)与地壳厚度呈正相关关系，CaO含量(Ca6.0)与地壳厚度呈负相关关系(Plank and Langmuir, 1988)，虽然已有众多研究，但地壳厚度如何控制火山岩中的K₂O、Na₂O和CaO含量的背后机制尚不明确。

通过观察全球火山岩CaO、Na₂O、K₂O的归一化联合概率密度分布图(图 2k-y)，发现年轻-壳源子集没有任何统计规律，而年轻-幔源子集部分存在一定的统计规律。年轻-幔源部分数据的CaO的联合分布图显示，分布随着地壳厚度增厚存在明显的线性减小趋势(图 2n)；Na₂O分布不存在单调规律，以50km为界，前半段呈现线性递增，后半段呈现线性递减趋势，误差相对较大，规律不明显(图 2s)；K₂O含量存在明显的随地壳厚度指数增加的趋势(图 2x)。

前人研究认为，Ce/Y和La/Yb代表轻、重稀土分异程度(Mantle and Collins, 2008; Chiaradia, 2015)；Sm/Yb、Dy/Yb的比值，代表中、重稀土分异程度(Mamani et al., 2010); Sr/Y的比值，代表斜长石、石榴子石和角闪石的部分熔融和分异过程(Chiaradia, 2015; Profeta et al., 2015)。前人做区域性研究发现，Sr/Y随地壳厚度增加呈现线性增加规律(Chapman et al., 2015)，La/Yb和Ce/Y随地壳厚度增加呈现指数增加规律(Mantle and Collins, 2008；Profeta et al., 2015)。但是高Sr/Y和La/Yb也可能由多种过程导致(He et al., 2011)，其指示作用有很多争议。本文通过全球火山岩样品的联合概率密度分布图分析(图 3)，年轻-壳源子集的Ce/Y，La/Y和Sm/Yb可能存在随地壳厚度增加而指数增长的规律(图 3g, l, v)；年轻-幔源子集的Sr/Y有随地壳厚度增加而线性增加的趋势，但是由于数据离散度较大，规律不明显(图 3d)。

本文根据统计图发现的规律，进一步建立地化指标与地壳厚度的函数模型。由于数据集中异常离群数据较多，各个地壳厚度区间数据量差异较大，无法直接用全体数据来拟合函数关系，因此，本文将各个地壳厚度区间统计出来的概率密度函数的峰值提取出来建立与地壳厚度的函数关系(图 4)，通过回归分析建立了对应的拟合方程。各个地化指标的回归方程、回归系数(a、b)、确定系数(R²)、拟合标准差(RMSE)见统计表 1。确定系数R²是判别拟合方程参考价值的重要指标，衡量了回归方程整体的拟合度，R²越接近1，说明回归方程对观测值的拟合程度越好。通过对比6个回归方程的R²(表 1)，除幔源Sr/Y回归方程的R²为0.6较低外，其他5个回归方程的R²均达到0.7以上，在统计学中可以认为是具有显著相关关系。本文将进一步讨论这5种地化指标与地壳厚度存在相关性的原因。

图 4

Fig. 4

图 4 地化指标与地壳厚度关系拟合曲线 实线为拟合曲线，虚线为95%置信边界，误差棒为拟合数据和数据估计误差 Fig. 4 Fitting curve between geochemical index and crustal thickness The solid line is the fitting curve, the dotted line is 95% confidence boundary, the error bar is the fitting data and data estimation error

表 1

Table 1

表 1(Table 1)

表 1 地化指标与地壳厚度的回归方程和拟合优度参数 Table 1 Regression equation of geochemical index and crust thickness and goodness of fit 地化指标 拟合方程 a b R2 RMSE 来源 CaO y=a×x+b -0.07952 11.7 0.8458 0.6382 幔源 K2O y=a×exp(b×x) 0.2023 0.0497 0.9553 0.5215 幔源 Sr/Y y=a×x+b 0.6352 2.244 0.6016 9.7156 幔源 Ce/Y y=a×exp(b×x) 0.6482 0.03248 0.7435 1.0556 壳源 La/Yb y=a×exp(b×x) 3.582 0.0325 0.7548 5.7633 壳源 Sm/Yb y=a×exp(b×x) 0.6523 0.03092 0.8319 0.7408 壳源

表 1 地化指标与地壳厚度的回归方程和拟合优度参数 Table 1 Regression equation of geochemical index and crust thickness and goodness of fit

前人通过岛弧幔源火山岩发现CaO、K₂O和Na₂O可以很好的指示地壳厚度，且给出了相关关系公式。本文通过对全球幔源火山岩地球化学数据分析，同样发现其主量元素CaO，K₂O与地壳厚度有显著的相关关系，而Na₂O与地壳厚度不存在单调的函数关系。对于基性火山岩全岩CaO，K₂O含量与地壳厚度的相关性内在原因，仍然存在较大争论。Condie and Potts (1969)、Condie(1973, 1976)等认为与地壳混染作用有关；Dickinson and Hatherton (1967)、Hatherton and Dickinson(1968, 1969)、Plank and Langmuir (1988)等则认为主要与不同压力下幔源源区部分熔融的比例有关。然而考虑到幔源岩浆产生的部位可能是岩石圈地幔或者软流圈地幔，其源区熔融的深度与地壳厚度并无直接关系，因此源区熔融导致的CaO、K₂O与地壳厚度的相关性可能不受源区深度的影响，而与幔源岩浆上升过程中和地壳的相互作用有关。本文提出一种可能的解释，即地壳与幔源岩浆的作用过程是导致CaO、K₂O与地壳厚度相关显著的主要原因。地壳越厚，幔源岩浆到达地表前与途经地壳的穿透路径和作用时间将越长。Ca与地壳厚度呈现负相关；而K则呈现指数正相关。其可能的解释是：幔源岩浆具有黏度低、温度高的特点，容易与地壳物质发生化学混合(张旗等，2007)。全球平均地壳(含量1.8%)相对于全球平均地幔(含量200×10^-6)富集K元素，因此在幔源岩浆上升过程中地壳物质的混合作用导致熔体K含量升高，越厚的地壳导致增长的作用时间和路径进一步促进熔体的K含量提升。相对于K，地幔与地壳的Ca含量在同一数量级范围(Rudnick and Gao, 2003；Salters and Stracke, 2004)，幔源岩浆上升过程中虽然有部分的地壳Ca加入，但是幔源岩浆的低黏性同时有利于结晶分异的发生，富钙矿物，如辉石，斜长石的结晶分异过程会极大的降低熔体中Ca的含量，因而在幔源岩浆在穿透厚地壳时将增大和延长结晶分异矿物的比例和时间，导致Ca与地壳厚度的负相关性。

由于大多数壳源岩浆都起源于下地壳，源区深度与地壳厚度存在对应关系，因此，壳源岩浆的地球化学特征可能直接反映地壳厚度问题。本文研究发现全球的壳源岩浆Ce/Y，La/Yb和Sm/Yb与地壳厚度的相关性较为显著，并呈现为指数相关。

为了进一步分析微量元素与地壳厚度相关性的内在原因，本文统计了地壳中常见造岩矿物(石榴子石、角闪石、斜长石、单斜辉石、斜方辉石)中Sr、Y、La、Ce、Sm、Dy和Yb这些微量元素含量(图 5a)，统计信息见表 2。其中，斜方辉石各微量元素含量均较低，不会对熔体中微量元素变化带来较大影响；斜方辉石和角闪石，曲线形态相近，对岩浆中熔体相矿物有类似的影响；石榴子石、斜长石和角闪石，曲线形态差异较大，岩浆固体相中这些矿物的比例会直接对熔体中这些微量元素含量和比值造成影响。本文进一步统计了石榴子石、斜长石、角闪石、单斜辉石和斜方辉石中微量元素比值对数指标的分布(图 5b-f)。对比各种矿物的各项微量元素比值对数指标可以发现，石榴子石中各项指标显著低于斜长石，角闪石和单斜辉石指标分布较为接近，且居于石榴子石和斜长石中间；其中，Sr/Y指标(图 5b)，各种矿物含量差异最大，石榴子石中指标分布峰值比角闪石中的峰值低了4~5个数量级，角闪石和单斜辉石体现出差异，峰值相差约1个数量级；La/Yb指标和Ce/Y指标分布相近(图 5c, d)，石榴子石和斜长石峰值相差3~4个数量级，角闪石和单斜辉石分布相近；Sm/Yb指标(图 5e)，石榴子石和斜长石峰值相差约2个数量级，角闪石和单斜辉石分布相近。Dy/Yb指标(图 5f)，石榴子石、斜长石、角闪石和单斜辉石的分布峰值相差不到1个数量级，分布近乎一致。

图 5

Fig. 5

图 5 部分矿物的微量元素统计关系 (a)常见矿物中各微量元素含量统计；Sr/Y比值(b)、La/Yb比值(c)、Ce/Y比值(d)、Sm/Yb比值(e)和Dy/Yb比值(f)分布规律统计.图(b-e)图例与图(f)同 Fig. 5 Statistical relationships of trace elements in some minerals (a) statistics of trace elements in common minerals; regular statistics of the distribution of Sr/Y ratio (b), La/Yb ratio (c), Ce/Y ratio (d), Sm/Yb ratio (e) and Dy/Yb ratio (f) in common minerals; legends in Fig. 5b-e are same as in Fig. 5f

表 2

Table 2

表 2(Table 2)

表 2 常见造岩矿物部分微量元素含量(×10-6)数量级统计 Table 2 Statistics of magnitudes of trace element contents(×10-6)in common minerals 矿物名称 Sr Y La Ce Sm Dy Yb 角闪石 均值 2.34 1.39 0.92 1.40 0.70 0.65 0.25 标准差 0.66 0.46 0.65 0.63 0.50 0.48 0.46 单斜辉石 均值 1.78 0.97 0.42 0.90 0.30 0.29 -0.06 标准差 0.70 0.55 0.72 0.71 0.56 0.52 0.55 斜长石 均值 2.51 -0.53 0.62 0.73 -0.51 -0.89 -1.17 标准差 0.81 0.61 0.66 0.66 0.53 0.51 0.63 石榴子石 均值 -0.34 1.05 -1.22 -0.40 -0.01 0.30 0.18 标准差 0.58 0.51 0.62 0.62 0.41 0.50 0.46 斜方辉石 均值 -0.37 -0.19 -1.23 -0.89 -1.16 -0.75 -0.56 标准差 0.78 0.86 0.93 0.90 0.81 0.88 0.98

表 2 常见造岩矿物部分微量元素含量(×10-6)数量级统计 Table 2 Statistics of magnitudes of trace element contents(×10-6)in common minerals

一般在温度800℃环境下，压力超过1GPa，会有石榴子石出现。在地壳厚度较大地区，下地壳出现榴辉岩化(Wittlinger et al., 2009)，形成大量石榴子石，会极大改变岩浆熔体中Y和Yb的含量。在压力较低的地壳区域，角闪石、斜长石和单斜辉石的比例多少，也有能力影响岩浆熔体的微量元素比值指标。但是考虑到结晶相/残留相矿物组合，不仅受到温度和压力的控制，还极大的受到含水度的影响(Guerri et al., 2015)，可能会导致浅部地化指标的波动性变化，最终导致壳源火山岩微量元素比值指标构建回归方程的确定系数偏低。

本文研究所建立的地化指标-地壳厚度回归模型是基于大数据统计得出的，在数据处理时，我们观察到数据分布是具有多峰特性，不符合正态分布特点。前人研究中也发现了该问题，所以有的研究提出利用中位数代替平均值进行研究(Chiaradia，2014)，或者利用各种技术手段，如修改后的Thompson-tau统计方法剔除离群数据。本文将统计学中的核函数估计方法引入到地化大数据研究中，无需剔除异常数据，直接建立数据的概率密度分布函数，从而准确把握数据的分布规律，并测量概率密度分布函数的主峰位置和峰宽，来估计数据的集中趋势和离散特征。因此，如果其他研究者需要利用本文研究获得的地化指标-地壳厚度回归模型重建一个地区的地壳厚度变化时，同样需要大量样本数据来获得数据的概率密度分布，将数据集中分布的主峰位置对应的地化指标带入回归函数，利用单一或少量岩石样本的地化指标的均值来追寻地壳厚度变化可能带来较大的风险。此外，由于本文的统计结果是基于全球数据的统计结果，区域上也可能存在偏离的情况，例如，本文按照全球平均模型认为地幔是贫钾的，区域性富钾地幔产生的岩浆形成的火山岩则可能偏离本文的回归曲线。我们观察到在图 2x中，同样在30km地壳厚度附近出现了6%~8%高钾亮点区域，这有可能是与特殊区域的富钾地幔有关。

本文结合GEOROC数据库中岩石与矿物大数据和crust 1.0地球物理模型中的全球地壳厚度模型，分析了年轻壳源和幔源火山岩中主量元素指标K₂O、CaO、Na₂O和微量元素指标Sr/Yd、La/Yb、Ce/Y、Sm/Yb、Dy/Yb相对于现今地壳厚度的关系，通过建立各个地化指标随地壳厚度变化的归一化联合概率密度分布图，直观地发现了部分指标与地壳厚度之间的相关关系，建立了回归函数模型，并讨论了确定系数R²大于0.7的可靠地化指标可以指示地壳厚度的可能原因。通过统计、分析和讨论获得了以下结论：

(1) 幔源火山岩用以指示地壳厚度的地化指标包括K₂O和CaO，其中K₂O随地壳厚度增加呈指数增加，确定系数R²为0.96，CaO随地壳厚度增加呈线性减少，确定系数R²为0.85。其相关关系背后的机制与幔源岩浆到达地表过程中穿透路径的长度相关，而与岩浆的源区深度关系无关。幔源岩浆穿透地壳上升到地表过程中，地壳的薄厚控制了贫K幔源熔体从富K地壳中提取K元素的富集过程，同样控制了熔体中富钙矿物通过分异结晶不断减少的过程。

(2) 壳源火山岩中Sm/Yb, La/Yb和Ce/Y均表现出随地壳厚度指数增加的相关关系，确定系数R²分别为0.83，0.75和0.74。其相关关系背后的机制可能与岩浆源区深度有关，不同源区深度下的斜长石、角闪石、石榴子石的稳定域控制了部分熔融的岩浆中微量元素比值的高低。

(3) 相同地壳厚度情况下产生的火山岩地化指标统计分布与正态分布差异较大，在运用地化指标确定地壳厚度的研究中，需要较多的样本来确定数据的概率密度分布函数，数据集中分布主峰的地化指标而非平均值进行地壳厚度限定。依据单一或少量地化数据回溯历史地壳厚度变化情况，需谨慎对待。