由于区分不同地质环境下的岩浆活动对于岩浆形成过程研究的重要性，岩浆构造环境判别已经成为地球化学领域一个重要的研究方向。在野外，岩石遭受蚀变是常见现象，而蚀变强烈的岩石不利于判别构造环境。但是，组成岩石的某些矿物抗风化和蚀变的能力较强，常可以见到在强烈蚀变的岩石中还有单斜辉石的残晶保留，这些残晶可能保留了原始未受蚀变的岩石的某些信息，因此可利用单斜辉石的残留恢复原岩的成分、岩石成因及其形成的构造环境^[1-2]。Nisbet等收集了329个单斜辉石数据^[3]，开创了利用单斜辉石矿物成分恢复岩浆岩构造环境的先河，并取得了较好的结果，此后，许多学者也对单斜辉石产出的构造环境进行了研究，结果表明，在组成岩石的造岩矿物和副矿物中，以单斜辉石和尖晶石的效果较好，橄榄石、斜方辉石、斜长石、角闪石及副矿物锆石、榍石、磷灰石、帘石类（绿帘石、褐帘石、斜黝帘石等）等的效果较差^[3-5]。其中又以单斜辉石最受学术界的青睐，因为：①单斜辉石是镁铁-超镁铁岩中最常见的造岩矿物；②单斜辉石结晶顺序较晚，保留了岩浆演化时所经历的成分、温度、压力等的变化；③单斜辉石的化学成分往往与构造环境密切相关^[3]。因此，在矿物学研究中对单斜辉石的分析较深入。但是，由于研究手段和思路的限制，利用单斜辉石作为构造环境判别指标的研究虽然取得了一些进展，总体效果并不明显。早期的判别图只考虑少量矿物元素^[6-7]，地球化学研究由于条件限制只能采用抽样或典型的方法进行^[8-10], 如图 1可见常用的判别图^[11-14]。而基于少量特征的判别图方法在面对全体数据时就显得捉襟见肘。笔者先前的研究表明：①目前的玄武岩构造环境判别图都存在问题，尤其是主量元素判别图存在多解性、矛盾性和不准确性。②由于组成成分扩展和大多数早期的识别标志被突破，大洋中脊玄武岩（MORB）和洋岛玄武岩（OIB）不易区分。③岛弧玄武岩（IAB）和其他2种岩石的判别图存在重叠，只有包含Th，Ta(Nb)，Sr的判别图依然起效^[16-18]。与此同时，全部数据的判别图表明各类构造背景中的镁铁-超镁铁岩组成成分变化很大，先前的判别图表现出了它们的局限性^[19-20]。

图  1  文献中常用经典若干玄武岩构造环境判别图^[15]

a— Th/Yb-Ta/Yb图解^[11]；b—Ti/Y-Nb/Y图解^[11]；c—V-Ti图解^[12]；d—Hf/3-Th-Nb/16图解^[13]；e—TiO₂-10MnO-10P₂O₅图解^[14]

Figure  1.  Some classical magmatic environment discrimination diagrams

下载: 全尺寸图片 幻灯片

一些用于判别岩浆构造环境和从全部数据中找出最具有区分度元素的方法应时而出，基于全球聚合数据的机器学习方法进一步推动了这一领域的研究^[21-22]。在这些机器学习方法中，基于决策树的方法由于较高的可解释性和良好的分类性能被广大研究人员广为采用。Vermeesch^[2]和韩帅等^[21]分别使用决策树算法和随机森林算法对玄武岩构造环境判别问题进行了分析，朱林奇等^[24]和韩启迪等^[25]分别将AdaBoost算法和GBDT算法应用于岩性识别问题。基于决策树的方法在地球化学研究中具有良好的表现，但是对决策树类方法（尤其是boosting和bagging集成方法）在地球化学领域的比较分析仍是亟需研究的课题。

本文笔者基于在线GEOROC数据库(http://georoc.mpch-mainz.gwdg.de/georoc/)和PetDB数据库(https://www.earthchem.org/petdb) ^[26]中单斜辉石的主量元素，使用决策树、随机森林、Adaboost和GBDT方法对来自3种不同构造背景（OAB，IAB和MORB）的镁铁-超镁铁岩数据进行了对比研究，为进一步判别岩浆环境提供技术支撑。本文比较了4种树类算法使用单斜辉石地球化学元素数据在镁铁-超镁铁岩环境判别问题上的性能，为地球化学分类研究的树类机器学习算法选择提供参考依据，同时得出在本数据集中区分不同构造环境的5种主要元素。

1.   树类算法

大量机器学习方法在分类问题上有良好表现，如支持向量机（SVM）、极限学习机（ELM）、人工神经网络（ANN）、决策树及集成算法或深度学习算法^[27-29]。但是，大多数方法，如SVM，ANN，ELM和深度学习算法将原始输入映射到更高维的空间，以黑箱模型的方式工作，具有较低的可解释性^[30]。而决策树方法的树形结构生成的规则符合人类的决策方式，易于理解，因此决策树方法具有较高的可解释性^[31-32]，适用于岩浆构造环境判别问题的研究。

随着集成技术的出现，为进一步提高性能，结合决策树和集成技术的判别方法也得以发展。基于决策树的集成算法主要可分为2类：bagging算法和boosting算法。bagging算法，通过结合多棵协作决策分类树进行投票决策，通常能取得更好的泛化性能。随机森林方法是常见的bagging扩展变体。boosting算法是一种串行的集成算法，其在每轮迭代产生弱分类器后，加权求和所有的弱分类器来获得最终分类结果。该算法关注于上一轮错分样本，因此常表现出更好的分类性能。广泛使用的boosting算法是AdaBoost算法和GBDT算法。以下简单描述四种经典的树类算法，即决策树、随机森林、AdaBoost和GBDT。

1.1   决策树

决策树是一种具有强鲁棒性、低计算消耗、高可解释性等优点的机器学习方法，其被广泛应用于分类任务^[33]。近年来，许多决策树算法得以快速发展，如ID3，C4.5，CART等。在决策树算法中，决策规则以树的形式组织构造，其中非叶结点表示决策条件，将结点数据分割为更小的结点；叶结点为终止结点，表示决策规则的终点。图 2展示了CART算法的一个示例。图 2-a为决策树的一种分割模式，其中a₁和a₂表示输入特征，m₁, m₂, m₃, n₁, n₂为分割点。图 2-b表征图 2-a相应的判别图，当决策条件得到满足时，左子树工作，否则右子树工作。本次采用的决策树方法为CART算法。

图  2  CART示例

Figure  2.  An example of CART

下载: 全尺寸图片 幻灯片

1.2   随机森林

随机森林是一种实用的分类算法，该分类算法的类别确定通过多棵子树投票实现。作为一类基于树的bagging集成算法，随机森林构建多棵不同的决策树，其中每一棵子树拥有相同规模的训练数据。为避免过拟合，一方面随机森林采用依赖替换实例的bootstrap方法^[34]从原始数据采样获取每棵子树的训练数据，另一方面，随机森林通过随机采样划分每棵子树的训练特征。

1.3   AdaBoost

起初用于特征分类和回归的AdaBoost算法已经成为机器学习领域中最成功的分类算法之一。同bagging算法相比，AdaBoost算法更关注于错分类的样本，其通过串行的方式来构造学习器^[35]。即一个弱学习器先被构造起来，然后下一个分类器的样本权重基于在先前构造中错分类样本来调整。通过重复“构建与调整”过程，一系列侧重于不同样本的学习器被构建起来，而最终的分类结果通过先前的弱分类器加权得到。不同于随机森林，集成AdaBoost算法中的多棵决策树分别拥有不同的权重。

原始的AdaBoost算法只关注于二分类问题，而SAMME和SAMME.R算法是2种常见的可用于多分类任务的AdaBoost变种算法，这2种算法的性能都已在实践中得以验证^[35]。本文使用性能更好的SAMME.R方法作为AdaBoost的实现进行三分类构造环境判别问题分析。

1.4   GBDT

梯度提升决策树（GBDT）在许多领域都属于一种高泛化性能技术。不同于AdaBoost，GBDT方法主要关注损失函数而不是错分类的样本，最终的强分类器只与前一个学习器有关。需要注意的是，GBDT中的弱学习器必须是决策回归树。

梯度提升决策树由多棵决策树组成，所有树结果的累计输出作为最终决策结论。梯度提升决策树在每次构造新的弱分类器时，关注于之前所有树的结论和残差，通过拟合得到当前残差回归树。与单决策树相比，梯度提升决策树通过弱化训练数据的微小变动对弱分类器的影响来支持最终强分类器的决策结论。

2.   实验

本文对4种树类算法，即决策树、随机森林、AdaBoost和GBDT在岩浆判别环境问题上的性能进行评估。

2.1   岩石数据

本文中使用的实验数据主要来自于2个在线数据库：GEOROC和PetDB，它们富含全球新生代镁铁-超镁铁质岩石单斜辉石数据^[32]。数据主要包含来自3种环境的玄武岩，即岛弧玄武岩、洋岛玄武岩和大洋中脊玄武岩，详细数据的地理分布如图 3所示；数据的降维可视化可由基于t-SNE的散点图得出，如图 4所示。实验数据集包含13种主量元素：SiO₂, TiO₂, Al₂O₃, Cr₂O₃, Fe₂O₃, TFeO, FeO, CaO, MgO, MnO, NiO, K₂O, Na₂O。实验使用python3.6环境在一台Intel(R) Core(TM) i5- 4200M, 2.50GHZ, 8.00GB的机器上进行。

图  3  镁铁-超镁铁岩中单斜辉石分布

IAB—岛弧玄武岩；MORB—大洋中脊玄武岩；OIB—洋岛玄武岩

Figure  3.  Distribution for clinopyroxene of mafic-ultramafic rocks

下载: 全尺寸图片 幻灯片

图  4  基于t-SNE的散点图

Figure  4.  The t-SNE-based scatter diagram

下载: 全尺寸图片 幻灯片

在进行分类任务前，需要对原始数据进行有效的清洗^[17]。本次实验采用的数据集在初步清洗后仍然存在部分数据缺失，因此需要进行数据填充。由于空缺位的实际元素含量远低于测量仪器范围，通过结合专家知识和常用策略，本文使用0.001填充缺失值^[36]。经过初步数据清洗后，共获得大洋中脊玄武岩有效数据795条，洋岛玄武岩有效数据198条，岛弧玄武岩有效数据329条，各项数据统计信息见表 1。

表  1  镁铁-超镁铁岩中单斜辉石主量元素统计信息

Table  1.  Major element content of clinopyroxene in mafic-ultramafic rocks in the dataset

主量兀素	IAB(岛弧玄武岩）				OIB(洋岛玄武岩）				MORB (大洋屮脊玄武岩）
	数据量	平均数/%	屮位数/%		数据量	平均数/%	屮位数/%		数据量	平均数/%	屮位数/%
SiO2	329	52.09	52.30		198	48.95	49.92		795	51.87	51.73
TiO2	324	0.36	0.20		198	1.72	1.24		784	0.16	0.10
Al2O3	329	3.63	3.62		198	4.45	3.67		795	4.34	4.55
Cr2O3	296	0.68	0.69		135	0.46	0.42		790	1.20	1.23
Fe2O3	52	1.27	0.99		1	3.34	3.34		10	1.46	1.65
TFeO	254	3.60	2.81		184	6.13	6.78		225	2.55	2.49
FeO	75	3.28	3.11		14	5.99	5.85		570	2.80	2.69
CaO	329	22.49	22.56		198	22.22	22.33		795	22.16	22.32
MgO	329	16.37	16.52		198	14.38	15.34		795	17.07	17.18
MnO	307	0.10	0.10		192	0.12	0.13		789	0.09	0.09
NiO	171	0.05	0.03		12	0.03	0.01		601	0.05	0.05
K2O	213	0.01	0.00		70	0.01	0.00		180	0.01	0.01
Na2O	320	0.43	0.39		198	0.55	0.34		778	0.31	0.18

下载: 导出CSV 

| 显示表格

由表 1可见，数据集中各类玄武岩的样本比例不平衡，为缓解样本不平衡对模型分类性能的影响，本文主要从以下角度考虑问题: ①决策树方法通常在不平衡数据上表现良好，而集成算法自身的加权特性使集成方法能更好地训练不平衡数据，本文使用的方法对于处理不平衡数据具有天然的优势；②对训练集的分类进行类加权处理，即每个类的输出概率均嵌入成本误差信息；③在性能度量时，使用多种度量指标，混淆矩阵、F1值等能更好地衡量模型在不平衡数据集上的表现。

2.2   性能测量

准确性（Accuracy）是分类问题性能评估的直观选择。此外，宏平均精度（MaP）、宏平均召回率（MaR）和宏平均F1值（MaF）也是多分类问题常用的指标。MaP表示每一类样本精确度的算术平均值，MaR表示每一类样本召回率的算术平均值，MaF是基于MaP和MaR的调和平均值。本次实验使用以上指标评估4种树类方法的性能，此外，本文使用混淆矩阵和散点图来更直观地展示实验分类效果。

2.3   参数选择

模型参数选择是机器学习算法建模过程中不可避免的过程。考虑到本文中模型参数数量较多，使用贝叶斯优化方法进行参数调整^[37]。实验最终确定参数如表 2所示。

表  2  镁铁-超镁铁岩中单斜辉石主量元素参数设置

Table  2.  Parameter settings of clinopyroxene in maficultramafic Rocks

参数	决策树	随机森林	AdaBoost	GBDT
max features	0.55	0.36	0.08	0.03
max depth	21	6	39	27
min samples split	2	5	3	4
min samples leaf	4	3	2	2
n estimators	-	90	210	710
learning rate	-	-	0.123	0.008
subsample	-	-	-	0.52

下载: 导出CSV 

| 显示表格

其中：“max_features”为决策树在进行最优分割时使用的特征占总特征的比例，参数搜索范围为[0.1, 1]；“max_depth”为树的最大深度，搜索范围为[5,50]；“min_samples_split”为分割决策树结点所需最小样本量，参数搜索范围为[2,6]；“min_samples_leaf”为叶结点即一条决策规则终止时此结点所需的最小样本量，参数搜索范围为[1,5]；“n_estimators”为集成算法中子树的数量，参数搜索范围为[20,1000]；“learning_rate”控制boosting算法中每棵树的贡献度，搜索范围为[0.1, 1]；“subsample”表征GBDT算法中用于拟合基学习器所使用的样本比例，搜索范围为[0.5, 0.8]；“-”代表该算法没有相应的参数需要调整^[38]。

3.   实验结果

本部分探讨树类模型在单斜辉石主量元素数据集进行10次独立重复实验的对比结果。首先展示基于决策树算法生成的决策图，之后基于随机选择的一次实验结果绘制混淆矩阵，然后列出了4种模型的特征重要性排序，并基于最重要的2个特征绘制散点图以更直观地展示模型的性能，之后采取10次重复实验的平均指标对比了4种模型的分类性能。

3.1   基于决策树的决策图

同其他基于神经网络机器学习的算法相比，决策树算法的一个巨大优势是其优秀的可解释性。决策树图可以通过树形结构清晰地展示其决策过程。考虑到生成决策树的规模，图 5只展示了部分决策树的决策图。

图  5  基于主量元素的决策树（部分）

Figure  5.  Part of decision tree graph

下载: 全尺寸图片 幻灯片

图 5中的每个结点数据表征当前的统计信息。以根节点为例，TFeO≤4.575表示分类决策条件，这个判别条件将样本分割为“True”和“False”两部分。samples = 100.00%表示当决策到当前结点时所剩余的样本占全体样本的比例。gini = 0.667表示基于当前样本集的Gini指数。values = [0.333, 0.333, 0.333]是当前样本集中每类样本的比例。需要注意，Gini指数为0的结点代表当前结点为叶结点，它们表征决策终点。

3.2   混淆矩阵

混淆矩阵是分类结果的数值化表示，每个类别（IAB / OIB / MORB）的实际样本数量和预测得出的样本数量都易从混淆矩阵中得出，混淆矩阵的判别结果是进行分类性能比对和散点图绘制的基础。图 6为4种方法在当前数据集建立模型所得的混淆矩阵，其中每一类的样本数量清晰可见，如图 5中对于决策树算法，有7个标签为IAB类的样本被决策树模型错误地预测为OIB类，而有55个标签为IAB类的标签被决策树模型预测正确。

图  6  基于主量元素的混淆矩阵（代号同表 1）

Figure  6.  Confusion matrix on major element data

下载: 全尺寸图片 幻灯片

3.3   特征重要性排序和散点图

为进一步探讨影响构造环境判别的主要特征元素，本文基于排列重要性方法^[39]对4种方法中的特征重要性进行了排序，排序结果如图 7所示。根据特征重要性排序的结果，在图 8中用2种最主要成分分别作为横、纵坐标轴构建散点图，展示4种算法在岩浆环境判别问题上的性能。

图  7  特征重要性排序

Figure  7.  Feature importance on major element data

下载: 全尺寸图片 幻灯片

图  8  基于主要成分散点图（代号同表 1）

Figure  8.  Feature-based scatter diagram on major element data

下载: 全尺寸图片 幻灯片

由图 7可知，决策树方法得出的最重要主量元素，为TFeO，Al₂O₃，TiO₂，Cr₂O₃，FeO；随机森林方法认为，影响算法进行构造环境判别最重要的主量元素为Cr₂O₃，TFeO，TiO₂，FeO，Al₂O₃；AdaBoost方法认为，影响算法判别最重要的主量元素依次为Cr₂O₃，TiO₂，TFeO，Al₂O₃，FeO；GBDT方法认为，最重要的主量元素为TiO₂，TFeO，Al₂O₃，FeO，Cr₂O₃。

由图 8可见，基于最主要的2种主量元素可以区分出大部分测试数据的构造环境，但是仍有小部分的测试数据在二维坐标系中交叠混杂。这种现象一方面说明了4种模型在区分岩石构造环境问题上的有效性，另一方面也说明了只使用2种特征来区分不同的构造环境会存在部分数据难以区分的现象，使用全数据模式研究地质学构造势在必行。

3.4   性能评估

为了评估树类算法的性能，本文设计了10次独立重复实验。表 3中的宏平均指数表征了决策树类算法的性能。从表 3的MaF和Accuracy的结果对比可见，boosting算法相较于bagging算法具有更好的分类性能，单决策树算法的效果劣于集成算法。

表  3  基于主量元素的性能指标

Table  3.  Performance indexes on major element data

测量指标	决策树	随机森林	AdaBoost	GBDT
MaP	0.8393(+/-0.0357)	0.9120(+/-0.0268)	0.9219(+/-0.0179)	0.9224(+/-0.0292)
MaR	0.8416(+/-0.0347)	0.8904(+/-0.0398)	0.9023(+/-0.0345)	0.9057(+/-0.0302)
MaF	0.8389(+/-0.0294)	0.8997(+/-0.0305)	0.9108(+/-0.0241)	0.9130(+/-0.0243)
Accurcy	0.8715(+/-0.0199)	0.9212(+/-0.0280)	0.9302(+/-0.0193)	0.9315(+/-0.0227)

下载: 导出CSV 

| 显示表格

4.   讨论

上述实验结果从决策图、混淆矩阵、性能评估等方面展示了树类机器学习方法在单斜辉石构造背景判别问题上的良好表现，验证了树类机器学习方法有效处理地质环境构造等有监督问题的可行性。实验结果表明，相较于传统的判别图方法，依托于全数据模式的机器学习分类模型能更充分地利用地质特征，挖掘数据深层次隐含信息，得出可靠结论。但同时，机器学习为了更大限度利用提供的数据信息，需要综合考虑更多的输入特性，使其模型结构通常较复杂。相对于判别图在地质构造环境研究中的直观性，机器学习模型更不易解释，因此机器学习模型在地质领域的深层次应用，还需要进一步探讨。毫无疑问的是，随着全数据模式在地质领域的兴起，机器学习方法在地质应用上的研究将更加广泛且深入。

本文中决策图、混淆矩阵和性能评估部分的结果也展示了不同树类方法在单斜辉石数据上的分类效果。集成方法基于多棵决策树综合考虑，随机森林通过多棵子树投票确定最终分类结果，boosting方法基于上一轮分类效果，改进模型以构建基分类器，相对于单决策树模型，随机森林和boosting方法均能明显提高分类准确度。随机森林基于随机划分的数据集和随机自助采样特征对每棵子树进行训练，最终投票选择分类结果，其通过降低模型方差提高模型性能。boosting方法关注上一轮分类器模型的分类结果，以迭代优化基分类器，最终加权多个基分类器获得分类结果，boosting方法通过降低模型偏差提高模型性能。2种方式均能有效提升模型性能，然而如图 7所示，本文使用的地质数据中不同的主量元素对判断样本所属类别的影响差异较大，因此相较于boosting分类方法，无权重投票的随机森林方法效果略差。毋庸置疑，多棵子树的集成方法的时间消耗较决策树方法更长，而集成方法中随机森林的多棵子树间并无联系，可以并行执行计算过程，因此随机森林的建模耗时比boosting方法短。

此外，图 7展示了4种方法得到的特征重要性排序，因为使用同一数据集，整体上4种方法表现出类似的排序结果，均认为TFeO, Al₂O₃, TiO₂, Cr₂O₃, FeO是区分3种构造环境的最主要元素，基于每种方法判定的2种最重要特征的数据分布图展示了特征重要性排序的结果可靠性和4种分类模型的性能。由图 8可见，得到的特征可以在一定程度上区分3种构造环境，但是因为只使用了2种主要的特征绘制散点图，不同类别的数据存在重叠部分，可见简单使用2种特征并不能完全划分出不同样本所处的构造环境，全数据模式应当在地质领域发挥更重要的作用。本文的特征重要性排序也展示了机器学习在地质方面的另一个应用，即特征重要性判别，虽然不同的机器学习方法得出的结果不同，但是结合专家经验，机器学习方法能够给地质领域提供更多的辅助功能。

5.   结论

本文基于单斜辉石的矿物元素数据，采用树类算法从可解释性、分类性能、特征重要性排序3个方面对镁铁-超镁铁岩的构造环境进行了分析。考虑到本文所用树类算法均基于决策规则，理论上可以完整绘制4种算法的决策树图，然而受限于集成树的规模和复杂性，本文只展示决策树算法的决策树图。就可解释性和构造过程的易理解性而言，决策树算法更胜一筹。就分类准确性而言，GBDT在当前数据集上表现更好。当对于模型复杂度和准确率均有所要求，随机森林方法对于岩浆环境判别问题更适宜。此外，4种方法对主量元素数据集的特征重要性排序在一定程度上表现出相似性，Cr₂O₃, TFeO, TiO₂, FeO和Al₂O₃可以被认为是识别构造背景的最重要主量成分。

本次研究取得了较好的结果，说明采用单斜辉石来鉴别镁铁-超镁铁岩形成的构造环境是可行的，这就提出了一个判别图研究的新方向。但是，如何进一步将研究结果应用于解决实际地球化学应用中，还需要进一步探讨。