基于玄武岩地球化学大数据的机器学习重建揭示地球岩石圈厚度演化及其对板块构造的启示

第一作者与机构：张振杰（中国地质大学（北京）地球科学与资源学院、地质过程与矿产资源国家重点实验室、深时数字地球前沿科学中心），陈国雄（中国地质大学（武汉）地质过程与矿产资源国家重点实验室），Timothy Kusky（中国地质大学（武汉）地质过程与矿产资源国家重点实验室），杨杰（中国地质大学（北京）地质过程与矿产资源国家重点实验室、深时数字地球前沿科学中心），成秋明（中国地质大学（北京）地质过程与矿产资源国家重点实验室、深时数字地球前沿科学中心；中山大学地球科学与工程学院）。

发表信息：本研究于2023年12月15日发表在《科学进展》（*Science Advances*）期刊上。

一、 研究背景与目标

本研究属于固体地球科学领域，聚焦于地球动力学和早期地球演化。地球岩石圈作为行星最外部的刚性层，记录了地球历史上的地质事件和过程。然而，随着时间回溯，地质记录愈发碎片化，导致我们对地球岩石圈厚度如何随时间演化，以及这种演化如何影响地球的构造系统知之甚少。传统的观点认为，由于软流圈自太古宙以来的长期缓慢冷却，全球岩石圈体积是缓慢而稳定增长的。但近年来的研究，包括一些克拉通（如华北克拉通）经历显著减薄甚至破坏的证据，以及对地壳厚度周期性变化的认识，挑战了这一观点。因此，定量重建地球历史上岩石圈厚度的演化，对于理解地球热演化、板块构造的启动与发展、超大陆旋回等根本性科学问题至关重要。

然而，直接重建岩石圈厚度演化面临巨大挑战。来自地球深部的地幔捕虏体虽然能直接揭示岩石圈的成分和热状态，但样品稀少且时间分布不连续。玄武岩作为地球上分布最广、最丰富的火成岩，其原始熔体来源于软流圈，其化学成分被认为对熔融压力（深度）高度敏感，而对地幔源区成分变化的敏感性相对较低，因此玄武岩的主量元素氧化物组成有望成为估算其形成时岩石圈厚度的有效指标。

本研究的目标是：利用全球玄武岩岩石地球化学大数据，通过机器学习方法，重建过去38亿年以来全球岩石圈厚度的演化历史；分析岩石圈厚度变化与地壳变质作用热梯度（温度/压力比，T/P）记录之间的关联；并基于这些关联，探讨地球动力学体制（特别是板块构造）的长期演化过程。

二、 详细研究流程与方法

本研究是一个典型的数据驱动型研究，整合了全球地球化学数据库、高温高压实验数据、机器学习算法和地球动力学分析。其工作流程主要包括以下几个关键步骤：

1. 数据收集与预处理 * 研究对象与样本量：研究团队从EarthChem等全球地球化学数据库中，系统收集了年龄跨度为38亿年至今的玄武岩样品数据。经过严格筛选（例如，限定MgO含量在7-17 wt%以最小化分离结晶影响，剔除氧化物总量低于95 wt%或烧失量高于3 wt%的可能蚀变样品，并利用3倍四分位距法剔除每个1亿年时间窗口内的异常值），最终构建了一个包含24,194个玄武岩全岩氧化物分析数据的全球数据库（Data S1），这些样品在全球空间分布上相对均匀。 * 实验数据训练集：为了建立玄武岩成分与熔融压力之间的关系模型，研究者从LEPR（Library of Experimental Phase Relations）数据库和相关文献中，编译了1392个与橄榄石和斜方辉石平衡的实验玄武质熔体成分数据（Data S2）。这些数据包含了已知的压力（Pf）和温度（Tf）条件，作为机器学习模型的训练和测试基础。

2. 变量选择与因果关系分析 * 在利用玄武岩成分估算岩石圈厚度（与Pf成正比）时，一个关键问题是是否需要同时考虑熔融温度（Tf）。为此，研究采用了因果森林（Causal Forest, CF） 这一机器学习方法进行因果推断分析，以评估Tf对Pf估计的因果效应（条件平均处理效应，CATE）。 * 分析方法：构建了两个CF双机器学习估计器，分别以Tf为处理变量预测Pf，和以Pf为处理变量预测Tf。使用合成少数类过采样技术（SMOTE）处理数据不平衡问题，并将数据分为训练集（80%）和测试集（20%）。 * 结果与决策：分析显示，Tf对Pf回归的CATE值仅为0.013 GPa（相当于岩石圈厚度约0.43 km），远小于Pf本身的数值范围；而Pf对Tf回归的CATE值则为44°C。这表明Pf对Tf回归影响很小，而Tf受Pf影响较大。因此，在估算岩石圈厚度时，可以忽略Tf变量，仅使用玄武岩的主量氧化物成分作为模型输入特征。这一结论与实验岩石学的观察一致。

3. 机器学习模型构建与岩石圈厚度重建 * 核心算法：研究采用随机森林（Random Forest, RF） 回归算法来建立玄武岩氧化物成分（SiO₂, TiO₂, Al₂O₃, FeO_t, MnO, MgO, CaO, Na₂O, K₂O）与熔融最终压力（Pf）之间的非线性映射关系。 * 模型训练与验证：使用上述1392个实验数据训练RF模型。通过随机搜索和交叉验证进行超参数优化。为了评估模型预测的不确定性，研究者重复了“训练-预测”过程100次，每次随机将实验数据划分为70%训练集和30%测试集，并构建100个独立的RF模型。最终每个玄武岩样品的Pf值为100次预测的平均值，方差用于评估不确定性。模型在测试集上表现出很高的预测精度（R² > 0.86），平均预测方差对应的岩石圈厚度不确定性约为11.4公里。 * 厚度计算与应用：将训练好的RF模型应用于全球24,194个天然玄武岩样品，预测每个样品形成时的Pf值。岩石圈厚度通过公式“厚度 = Pf / (g * ρ)”计算，其中g为重力加速度，ρ取地幔密度3.25 g/cm³。为了获得平滑的演化曲线，对计算出的Pf值（及从文献汇编的564个变质岩T/P数据）进行时间序列分析，采用每1亿年为一个窗口、两侧重叠2亿年的滑动平均，并用局部加权散点平滑法（LOWESS）进行拟合。

4. 时间序列分析与相关性研究 * 同步性与因果关系分析：计算了重建的岩石圈厚度序列与变质岩T/P比值序列之间的全局皮尔逊相关系数（PCC）。进一步，为了探究两者在不同时间尺度上的协同变化，计算了局部PCC（采用0.6 Ga和1.2 Ga的时间窗口，分别对应超大陆旋回和超大洋旋回的周期）。 * 周期性分析：对两个时间序列进行了连续小波变换（基于Morlet小波），以检测其周期性信号及其随时间的变化。并进行了小波相干性分析，以评估两个序列在不同时间和频率尺度上的一致性及相位关系。 * 因果时序分析：通过计算时滞互相关函数（TLCC） 和窗口化TLCC，来量化岩石圈厚度与变质T/P比值这两个非平稳时间序列之间的同步性和因果关联方向，寻找是否存在一个序列领先于另一个序列的系统性时滞。

5. 模型验证与结果稳健性检验 * 与现代观测对比：将模型重建的现代岩石圈厚度平均值（~73 km）与全球岩石圈模型Litho1.0给出的现代活动岩石圈平均厚度（~71 km）进行对比，结果高度一致。 * 区域验证：将模型应用于中国东部新生代玄武岩数据，重建的岩石圈厚度空间变化趋势（从东向西增厚）与地球物理观测（如地震学得到的岩石圈-软流圈边界深度）相符。 * 模型敏感性测试：为了测试模型在高压力环境下的表现，研究者还训练了一个仅使用低压（ GPa）实验数据的RF模型。新模型预测的压力值相对较小，但与原始模型得出的相对变化趋势一致，仅在早太古代由于样本稀少出现差异，证明了主模型的稳健性。

三、 主要研究结果

1. 全球岩石圈厚度演化历史：重建结果显示，地球活动岩石圈厚度在过去40亿年间在60-140公里范围内波动，并未表现出长期、单向的增厚趋势。识别出四个显著的岩石圈减薄事件，分别发生在： * 古太古代（3.5-3.2 Ga） * 早元古代（2.5-2.0 Ga） * 新元古代早期（1.0-0.8 Ga） * 显生宙（0.5-0.3 Ga，即古生代） 在这些减薄事件之间，则存在岩石圈增厚的阶段。其中，在中元古代晚期（约13亿年前）达到平均厚度峰值约140公里。

2. 岩石圈厚度与变质作用的热耦合关系： * 全局正相关：岩石圈厚度与地壳变质热梯度（T/P）从38亿年前至今存在显著的全局正相关性（PCC = 0.37， p = 5×10⁻¹⁴）。这种相关性从约23亿年前开始增强（PCC > 0.5），并在16亿年后达到0.7的高值，表明两者耦合程度随时间增加。 * 周期性协同：小波分析揭示，岩石圈厚度和T/P序列均存在约0.3 Ga、0.6 Ga和1.2 Ga的显著周期信号，分别对应于威尔逊旋回、超大陆旋回和超大洋旋回。小波相干性分析显示，两个序列在0.6 Ga（超大陆旋回）周期上表现出强烈的同相位一致性。 * 局部相关性变化：以0.6 Ga为窗口的局部PCC分析显示，其变化与超大陆旋回耦合：在超大陆聚合期（如2.8-2.3 Ga， 1.0-0.8 Ga之后），PCC值升高（正耦合）；在超大陆裂解期，PCC值下降（解耦）。 * 因果关系与时滞：时滞互相关分析表明，T/P序列领先于岩石圈厚度序列约0.2 Ga时，两者相关性最高。这意味着当前的岩石圈厚度可以决定约2亿年后的地壳变质条件，揭示了岩石圈厚度对地壳热状态的调控作用存在一个时间延迟。

3. 对地球动力学过程的启示： * 正耦合机制：研究认为，在板块俯冲主导的超大陆聚合期，岩石圈的水平缩短和增厚会产生“热毯效应”，阻止地幔热量散失，导致地温梯度升高，从而产生高T/P的变质作用。因此，岩石圈增厚与高变质热梯度正相关。 * 解耦机制：在超大陆裂解期，热的地幔柱上涌会使岩石圈减薄，但同时也会提高地温梯度，导致两者解耦。更重要的是，研究指出在中太古代（>2.8 Ga）之前，岩石圈厚度与T/P呈负相关。这可能是由于早期地球以大洋内俯冲-增生体制为主：虽然俯冲会导致岩石圈局部增厚，但当时的大陆尚未大规模形成，岩石圈是破碎的，缺乏连续的热毯，俯冲带入的冷地壳物质反而促进了软流圈冷却，因此岩石圈增厚并未伴随地温梯度的同步升高。

四、 研究结论与意义

本研究通过创新的“大数据+机器学习”方法，首次定量重建了地球过去38亿年的岩石圈厚度连续演化历史，并揭示了其与变质作用记录的深层联系，对理解板块构造的启动与演化提出了新见解：

1. 板块构造启动时间：岩石圈厚度在早于28亿年前就表现出显著波动（120-80公里），这与停滞盖层（stagnant lid）模式下地幔柱引起的单向、波动较小的厚度变化模式不符。结合早期存在的负相关到正相关的转变（~2.8 Ga），研究支持某种形式的动态盖层（即板块构造）在冥古宙或始太古代就已开始运作，而非从停滞盖层逐渐演化而来。
2. 俯冲体制的转变：~2.8 Ga是一个关键的转折点。在此之前，岩石圈厚度与变质T/P负相关，指示以大洋内俯冲-增生为主的构造体制。在此之后，两者转为强正耦合，标志着现代威尔逊旋回（伴随超大陆聚合与裂解）开始占据主导地位，这与新太古代-古元古代碰撞造山带广泛分布的地质记录一致。2.8-2.3 Ga期间的第一个高PCC区间，可能代表了促进最早超大陆（Kenorland）聚合的现代威尔逊循环俯冲。
3. 岩石圈演化与超大陆旋回的关联：四个显著的岩石圈减薄事件分别对应于已知或推测的超大陆/超克拉通裂解期（如Kenorland, Columbia, Rodinia, Pangea的裂解），而增厚期则对应聚合期。特别是中元古代晚期岩石圈增厚至峰值（~140 km）产生的强烈热毯效应，可能触发了导致Rodinia超大陆裂解的超地幔柱事件。随后的减薄打破了热毯，促进了地球内部热量的释放，冷却的软流圈使得刚性、冷的俯冲带成为可能，形成了以中、低T/P变质岩对为特征的现代板块构造体制。
4. 方法论价值：本研究展示了将全球地球化学大数据、实验岩石学约束与先进的机器学习算法（RF, CF）相结合，解决经典地质学难题的强大能力。所建立的基于玄武岩成分的岩石圈厚度估算模型，为研究地球和其他类地行星的长期演化提供了新的定量工具。

五、 研究亮点

1. 开创性的数据驱动重建：首次利用全球尺度的玄武岩地球化学大数据，系统重建了地球岩石圈厚度的长期演化曲线，提供了超越局部地质记录的全新全球视角。
2. 创新的方法融合：成功将随机森林回归模型应用于岩石学反演问题，并引入因果森林进行变量筛选和机理分析，体现了数据科学与传统地球科学交叉的前沿性。
3. 揭示关键转折点：明确了~2.8 Ga是地球动力学体制从早期大洋内俯冲为主向现代威尔逊旋回为主转变的关键时间节点，为早期板块构造样式的争论提供了重要的定量约束。
4. 建立跨圈层耦合关系：通过严谨的时间序列分析（PCC、小波分析、TLCC），不仅证实了岩石圈厚度与地壳变质作用之间的统计关联，还揭示了其因果关系（岩石圈厚度决定约2亿年后的变质条件）和与超大陆旋回的协同周期，深化了对固体地球系统内部耦合过程的理解。
5. 挑战传统观念：研究结果不支持“岩石圈自太古宙以来持续缓慢增厚”的传统观点，而是揭示了其动态的、周期性的演化历史，并与超大陆旋回紧密相关，强调了全球构造事件的节律性。