年轻火山岩样品的K-Ar、⁴⁰Ar/³⁹Ar定年结果的准确性受到多种因素的限制。这种限制主要来自于两个方面的不确定性，一方面是样品属性，主要是样品K-Ar同位素体系的均一性，如包裹体、斑晶含量、颗粒形态、粒度大小等，也包括氩丢失或过剩带来的样品中K-Ar体系的不均一；另一方面是仪器测量过程中的不确定性，包括仪器的稳定性、准确测定的量程范围等。

权伍勋等(2013)通过对误差传递的详细数学分析，发现⁴⁰Ar/³⁹Ar年龄的相对误差与年龄本身成反比，样品年龄越年轻，测年结果受到初始⁴⁰Ar/³⁶Ar比值的偏差和测试过程的偏差的影响越大，甚至出现呈指数增长的偏差(周晶等, 2013)。但遗憾的是，样品初始⁴⁰Ar/³⁶Ar比值与Nier值之间的偏差广泛存在，即样品的K-Ar同位素体系呈现一定程度的不均一，这会导致同一件样品重复测试的结果之间存在差异(Singer et al., 2008; Matchan and Phillips, 2014)。周晶等(2013)用激光全熔⁴⁰Ar/³⁹Ar法对国际上广泛使用的10余种标准样品进行研究，发现多数标准样品的表观年龄在概率统计图上呈多峰或不对称的单峰分布，表明样品(0.1~0.3mg)同位素体系的不均一是普遍现象。

Acta Petrologica Sinica岩石学报2021, 37(8) 周鑫等: 望天鹅火山喷发时限的激光⁴⁰Ar/³⁹Ar年代学望天鹅火山是东北地区新生代晚期喷发的一座中心式火山，其规模仅次于长白山，了解望天鹅火山的喷发历史可以为研究第四纪环境变化提供重要的参考。已有的年代学数据以K-Ar年龄为主，多数集中在2~3Ma，且不同研究者获得的年龄差异较大(樊祺诚等, 1998; 陈晓雨等, 2008; Wei et al., 2007)。这种差异很可能是样品同位素体系的不均一和测试过程中的偏差所造成的，为约束火山喷发时代带来了很大的不确定性，甚至导致研究人员将测试偏差造成的年龄结果的多样，误以为是多次喷发。这一问题的根源在于，前人只对少数样品点进行分析，样品对火山喷发时限的代表性不够，无法给出准确的喷发时限。针对这一问题，本文通过火山垂向剖面的样品采集，运用激光全熔⁴⁰Ar/³⁹Ar定年方法，并引入高斯混合模型，对望天鹅火山剖面的测年数据进行更精细地处理，以期更好地限定望天鹅火山的喷发时代。

1 地质背景

望天鹅火山位于吉林省东南部的中朝边境，在长白山西南方向约32km处。望天鹅火山被认为是中心式喷发，是因为火山口向周围呈放射状分布。整座火山由火山盾和火山锥组成，主峰位于火山口的东南边缘，高2051.4m。火山锥覆压在火山盾之上，面积远小于火山盾，主体已被后期的地表过程所破坏。火山盾在我国境内东西长约100km，南北宽近50km，面积近4000km²，喷发规模仅次于长白山天池火山(图 1)。火山盾覆盖在侏罗纪地层之上，构成了望天鹅火山的主体，由溢流玄武岩所构成。台地顶面十分平坦，高度大约为1130m，厚度约300m，由中心向四周逐渐减薄(图 2)。

① 吉林省地质局. 1964.1∶200000区域地质图

火山地质调查表明，无论是火山盾，还是火山锥，其内部或二者之间均没有发现有沉积层和风化面，即明显的火山喷发间断，是短时期内喷发堆叠所形成。受东亚季风的影响，望天鹅地区水系发育，由火山锥向四周呈放射状流出，切割出V形河谷，将熔岩流垂向形态很好地暴露出来，显示火山盾分布区柱状节理构造很发育(图 3)。

望天鹅火山的喷发早于长白山火山，但关于喷发的确切时限，目前仍没有得到精确的年代学结果。刘嘉麒(1987)最早对望天鹅十八道沟玄武岩用K-Ar法进行定年，得到2.24±0.23Ma的年龄。刘祥等(1989)将望天鹅火山喷发分成3期，分别是长白期玄武岩(16.40±1.49Ma)；望天鹅期玄武岩(5.56±0.22Ma)；红头山期碱性岩(3.11±0.053Ma)。樊祺诚等(1998)获得的K-Ar年龄均属上新世早期(2.12~2.87Ma)，以此为根据，将望天鹅火山期次划分为长白期(？~2.87Ma)、望天鹅期(2.69~2.41Ma)和红头山期(2.12Ma)。陈晓雨等(2008)的K-Ar测年结果介于1.87~7.04Ma。其中，十五道沟国家公园获得玄武粗安岩样品年龄为2.10~1.86Ma；上覆火山锥顶部碱性流纹岩测年结果是3.14~2.23Ma；火山锥体下部玄武粗安岩却为7.04~4.43Ma。Wei et al. (2007)获得十九道沟玄武岩的K-Ar年龄为3.66±0.05Ma，火山锥体底部的火山岩K-Ar年龄为3.03±0.02Ma和3.25±0.04Ma。文献显示部分火山锥体顶部的年龄老于国家公园内火山盾上获取样品的年龄，这一现象难以理解。

2 样品采集与前处理

野外调查显示，望天鹅火山熔岩层内部没有发现可以识别的不连续面，我们认为其喷发是在相对较短的时间范围内完成的。地质填图工作发现，火山岩自火山口向周围呈放射状分布，根据这种几何形态，火山盾和火山锥的垂向剖面能够代表望天鹅火山喷发完整的时间跨度。因此，选择在露头出露良好的十九道沟剖面自上而下系统地采集了样品(图 1)。样品详细信息见表 1。

样品的显微结构显示(图 4)，火山岩具有斑状结构。斑晶有少量橄榄石，多数为斜长石、辉石，少数岩石中斜长石斑晶可以达到粒径5mm以上。基质为隐晶质或斜长石、辉石矿物微晶。已有的研究显示，橄榄石斑晶有过剩氩因素的干扰，在样品处理时应予以剔除(周晶等, 2013)。

样品粉碎后，选择无风化边、无包体、无裂隙填充物的颗粒，用去离子水进行超声波清洗，再用钢制碎样器进一步破碎。最后在双目镜下进行手工挑选，去除橄榄石斑晶，获得180~250μm粒径的样品。

接着，先后用去离子水和丙酮浸泡样品，同时超声波振荡5min，再用去离子水冲洗干净，低温(80℃左右)烘干。每一件样品称取50~100mg，和用于校正K、Ca、Cl诱发同位素的K₂SO₄、CaF₂、KCl，特别是监测照射过程中的快中子通量的标准样品，装到高纯铝箔中，密封于石英瓶中。测试标样是ACs透长石(Turrin et al., 1994; Renne et al., 1998; Nomade et al., 2005)和ZBH-25黑云母(桑海清等, 2006)，分别用于计算样品的中子活化参数(J值)和J值在石英瓶中的变化情况。将真空密封的石英瓶置于中国原子能科学研究院49-2反应堆B4孔道接受快中子照射，照射时间为6小时7分钟，快中子瞬时通量为5.73×10¹⁷n/cm²。

3 测试流程

将中子活化后的样品在显微镜下从铝箔纸中取出，放入干净的样品座中。样品座上有215个直径小于2mm的小孔，每个孔中移入5~10个颗粒，总质量约0.1~0.3mg。通过全自动控制程序，实现对每一个样品点的全熔，并获取每个点的Ar同位素测试数据(周晶等, 2013)。

测试工作在北京大学全自动全时标高精度激光⁴⁰Ar/³⁹Ar定年系统中完成。该系统由激光熔样系统、纯化系统、VG5400质谱氩同位素分离探测系统以及全自动控制系统4部分组成。测试流程参见(周晶等, 2013)。年龄计算中的衰变常数取Steiger and Jäger (1977)的推荐值。质量歧视因子D=1.0045±0.0013(McDougall and Harrison, 1999; Renne, 2000)。全自动测样由Berkley地质年代中心Alan Deino博士编写的Mass spec v.7.84程序控制，并进行数据处理。Ca和K照射产生的干扰反应的校正参数为：[³⁶Ar/³⁷Ar]_Ca=00.0002775±0.000025，[³⁹Ar/³⁷Ar]_Ca=0.0006633±0.0003535，[⁴⁰Ar/³⁹Ar]_K=0.0039448±0.0015421。测试期间仪器本底水平为⁴⁰Ar=(5.83E-16)±(2.33E-18)mol, ³⁹Ar=(1.94E-18)±(1.67E-19)mol, ³⁸Ar=(4.49E-19)±(1.80E-19)mol, ³⁷Ar=(7.20E-19)±(1.30E-19)mol, ³⁶Ar=(2.44E-18)±(9.03E-20)mol。望天鹅火山岩激光⁴⁰Ar/³⁹Ar测试数据显示³⁹Ar数量级在~10^-14mol，远高于10^-18~10^-19摩尔本底值。

4 望天鹅火山激光⁴⁰Ar/³⁹Ar测年结果

整个测试共获得132个有效数据结果。激光⁴⁰Ar/³⁹Ar表观年龄的分布范围为0.00±0.74Ma(±1σ)~3.23±0.23Ma(±1σ)。具体定年结果见表 2，原始数据见电子版附表1。

5 高斯混合模型的数学原理

统计学上把对具有相似特征的数据分类的方法叫做聚类分析，对于含有多个正态分布的复杂数据总体，高斯混合模型(Gaussian Mixture model，GMM)是公认的最优聚类方法(Jain and Dubes, 1988; Krishnan et al., 1996)，它是多个高斯分布函数的线性组合，其表达式如下(Fraley, 1998; Fraley and Raftery, 1998)：

其中，α_k是系数，α_k≥0，∑_k=1^Kα_k=1；ϕ(y|θ_k)是高斯分布密度，，称为第k个分模型。

在GMM模型中有3组参数需要估计：α_k，μ_k，σ_k。在进行参数估计时通常选用最大期望算法(Jain and Dubes, 1988; Fraley et al., 2012)，在计算过程中，首先对每个分量k设置α_k，μ_k，σ_k的初始值，然后计算分模型k对观测数据的响应度：

再根据算得的计算新的α_k，μ_k，σ_k。其中：

重复上述过程，直到收敛。通过迭代计算过程，计算机程序可以自动识别一组数据中包含的正态分布个数，每个正态分布对应的均值和方差，以及每个数据所属的聚类。

6 讨论 6.1 样品K-Ar同位素体系的均一性

火山岩剖面自下而上的年龄应具有由老到新的变化规律，但望天鹅十九道沟剖面样品的等时线年龄并没有呈现这一特征，而且所有样品的年龄在置信区间内都有重合，未显示存在多个喷发期次(图 5a)。样品的初始氩比值在剖面上也没有趋势性的变化，同时在置信区间内也都有重合(图 5b)。这表明，在仪器测试分辨能力下，6件样品是同时同源的，年龄和初始氩比值的差异是样品K-Ar同位素体系不均一所造成的，不能认定为样品来自不同期次的火山喷发，因此也不能依据单一样品的年龄结果限定火山喷发时代。

6.2 火山喷发时限的确定

在测量学上，同一对象的多次测量结果呈正态分布，同位素年代学方法也符合这一规律。如果火山只有一次喷发，表观年龄的概率分布一定呈正态分布(胡燕等, 2018)，如果有多次喷发，表观年龄的概率分布是多个不同均值的正态分布的叠加，呈混合高斯分布(Zhou et al., 2021)。因此，考虑到采样只是人为获取代表火山喷发时间的样品的方式，可以打破传统意义上样品独立意义的界限，在统一的时间尺度标准下，对有代表性意义的所有样品做统计分析。通过高斯混合模型对数据进行聚类分析，可以识别和恢复代表不同火山喷发期次的多个正态分布，用以确认不同喷发期次的年龄。基于这一思路，可以对望天鹅火山十九道沟剖面样品的132个表观年龄数据进行统计，探究其喷发期次。

借助R 3.4.0程序“mclust”软件包的“densityMclust”函数，对望天鹅火山的表观年龄用GMM模型进行了聚类分析，其结果如图 6所示。所有表观年龄包含两个正态分布，其均值分别为1.68±0.70Ma和2.58±0.32Ma。后者占数据点多，正态分布形态好，是望天鹅火山喷发事件年龄真值的反映。均值较小的一个正态分布可能蕴含了长白山火山喷发热事件导致的氩丢失的信息，这也与部分样品等时线年龄的初始⁴⁰Ar/³⁶Ar比值小于Nier值的特征相一致(表 2)。

2.58 ±0.32Ma是表观年龄的统计结果，在⁴⁰Ar/³⁹Ar方法中，表观年龄的计算直接使用Nier值代替初始⁴⁰Ar/³⁶Ar比值，会带来年龄结果的偏差，因此需要求取等时线年龄来代表火山喷发的时限(Deino, 2012; 权伍勋等, 2015)。能够构成正态分布是样品符合等时线计算所要求的同位素体系“同时”、“同源”、“封闭”的数学表达，同一正态分布的数据可以用于计算等时线年龄。望天鹅火山岩样品⁴⁰Ar含量为10^-13~10^-14mol，反等时线在表达样品真实年龄上更有计算优势。对构成正态分布的101个数据做反等时线，年龄为2.73±0.18Ma，初始氩比值为293.84±0.97(图 7)。由于数据包含了整个火山岩剖面的年龄信息，这一结果比单一样品的年龄更有代表意义，能够代表望天鹅火山喷发的时限。

令人惊奇的是，望天鹅火山喷发时间恰好与北半球大冰期事件(Northern Hemisphere Glaciation, NHG)相吻合。2.73Ma前后北半球高纬地区冰川急剧扩张，深海氧同位素快速变重，在亚欧大陆北部、北美洲北部、太平洋和大西洋北部等地都留下了许多生物地球化学和古海洋沉积记录(Maslin et al., 1998; Sosdian and Rosenthal, 2009; Rohling et al., 2014)。关于NHG的成因至今没有定论，Kennett and Thunell (1975)和Rea et al. (1995)曾提出由于新生代晚期北半球火山喷发频率增加，向大气中释放含硫气体，引起全球降温。张文防等(2019)对北太平洋风尘沉积的研究发现，在2.73Ma前后，北太平洋1208钻孔出现了大量火山灰沉积，可能蕴藏着火山喷发带来全球变冷的信息。望天鹅火山覆盖面积广，喷发规模大，很有可能是钻孔沉积物中火山灰的主要来源。

7 结论

野外地质调查和同位素年代学测试结果都显示望天鹅火山在现有仪器测试精度下只有一次喷发，其对应的时限是2.73±0.18Ma，与NHG事件的启动时间相一致。望天鹅火山分布面积大，其喷发很有可能是同时期北太平洋深海沉积物中火山灰的主要来源，也可能是NHG事件的气候动力因素之一，但关于二者之间的具体联系，还需要从火山喷发规模、释气量、岩浆成分，火山喷发气候效应的数值模拟，年龄结果的优化等方面开展进一步研究工作。

致谢      感谢两位审稿人提出的详细而有建设性的意见和建议。

通讯作者感谢26年前李继亮老师的帮助和指导，师恩难忘。