基于计算机断层扫描的火山岩气孔含量及大小分布特征无损快速分析

程荣， 钱生平， 孙添力， 周怀阳，*

(1.同济大学海洋与地球科学学院， 上海 200092;2.同济大学海洋地质国家重点实验室， 上海 200092;3.同济大学电子与信息工程学院， 上海 201804)

气孔构造是火山岩的主要代表性特征之一。硅酸盐岩浆在上升过程中，压力降低，挥发份过饱和后发生气体出溶，并经过气泡长大、合并、变形，或者气体逃逸的改造，形成火山岩的气孔构造[1-3]。根据岩浆中挥发份压力与环境压力平衡的原理，利用火山岩中气孔的大小分布可用来估计大陆溢流玄武岩地区的古大气压和古海拔[4-6]。对火山岩气孔特征的准确评估可以提供有关岩浆上升速度和喷发方式的信息，从理论上确定火山喷发的时间和喷发的性质，提高预测火山爆发的能力[7-9]。不仅如此，对火山岩气孔特征(气孔的大小、空间和形状分布，整体的连通性、定向性等)的研究，还可以帮助人们了解岩浆的原始挥发份含量，源岩(区)性质，岩浆的物理(黏性、流动性等)性质，以及气泡的形成和生长过程[2,10-12]。

人们很早就开始使用各种实验方法研究火山岩中气孔的物理特征。Houghton等[13]根据阿基米德原理，通过分别测量火山碎屑岩在空气和水中的质量计算岩石的孔隙率(水浸法)。Sahagian等[14]在水浸法的基础上，先使用低黏度聚合物(液体环氧树脂)充填气孔，然后用氢氟酸溶解岩石基体，最后通过对残余填充物进行计数和测量获得气孔的大小分布。Whitham等[15]使用树脂进行真空浸渍，证实浮石中所有气孔是相互连接的，并进一步通过氮吸收技术和压汞法分别得到气孔的表面积和不同孔径气孔的孔隙率分布。Shea等[16]基于二维图像数据，将气孔视为规则球体进行体视学转换来分析火山岩气孔的三维分布特征。随着光学和电子显微镜成像技术的应用和普及，人们从依据薄片的光学照片[17]和背散射电子图像[18]统计气孔大小分布和孔隙率，发展为使用立体扫描电镜(SSEM)构建样品表面的三维数字高程模型来观察气孔形态和计算气孔体积[19]。然而，这些方法通常很耗时，采集的数据量较少，并且会破坏样品的原有组成和结构，其推广应用受到一定的限制。

高分辨率X射线计算机断层扫描，也称工业CT，是一种现代化的材料检测技术。目前，该技术已经广泛应用于医学、工业、考古学、生态学等多个领域，并且正在逐步发展为地球科学领域的重要分析手段[20-21]。自Song等[22]的开创性工作以来，越来越多的研究使用工业CT讨论火山学问题，这些研究大多集中在气孔结构和大小分布的测量上，并强调相互贯通的气孔构造对岩浆渗透率和气体散失的影响，显著提高了人们对火山脱气动力学的理解[23-26]。对岩石气孔进行分析时，工业CT具有无损、快速、实时、高分辨率(微米到亚微米级)立体成像的特点，相比于传统的称重法、浸泡法、注胶法和分层磨片光学成像法具有明显的优势，可以让研究者对岩石的内部结构(即晶体、气孔和基质)进行三维空间的可视化和量化，并通过上述任何一种技术直接测量同一样品而获得可重复的结果和补充性研究[27-30]。

本文通过将工业CT扫描与计算机图像处理技术相结合，使用了三种分析火山岩中气孔含量及大小分布的软件或程序(VG Studio MAX、ImageJ、Matlab代码)，此过程能快速、便捷地提取出气孔可靠的形态学数据；同时，通过改进火山岩中气孔体积转化为挥发份质量分数的方法，定量评价挥发份在岩浆中的富集程度。

1 实验部分

1.1 样品来源与实验设备

实验样品为2018年“嘉庚号”无人遥控深潜航次在南海珍贝海山通过水下机器人(ROV)抓取到的气孔状玄武岩(图1a)。样品表面呈深灰色，块状构造，斑状结构，具轻微蚀变，基质主要由微晶组成，常见的斑晶矿物有橄榄石、辉石、长石等。由X射线荧光光谱仪(XRF熔片法)测试的全岩主量成分含量为：SiO244.63%，Al2O315.87%，Fe2O311.11%，MgO 6.89%，CaO 8.69%，Na2O 3.96%，K2O 2.76%，P2O50.83%，MnO 0.28%，TiO22.94%。样品中的气孔构造十分发育，大多数气孔为扁圆形或不规则，且相互独立，与一般玄武岩中多为偏圆形气孔有所不同。少数气孔有合并现象，在空间中随机分布，既没有聚类也没有规则的间距。气孔尺寸相差较大，大的可达2～3mm，小至微米级。考虑到样品大小与CT扫描的分辨率密切相关，以及获取一系列规则切片的需要，本实验将检测样品切割成19mm×22mm×43mm的四方柱体。

a—“嘉庚号”无人遥控深潜航次ROV水下采样现场情况，采样水深1488m；b—微焦点X射线CT扫描成像系统示意图，XYZ线代表三维空间坐标轴，探测器用于接收投影数据。图1 气孔状玄武岩样品和检测设备Fig.1 Vesicular basalt samples and testing equipment

a—样品三维CT重建效果； b—岩石中的气孔提取过程，蓝色为渲染的气孔； c—去除岩石基体后的气孔三维立体分布； d—CT切片的灰度直方图，气孔与岩石基体灰度差异大，容易区分； e—CT切片图，样品区域内黑色为气孔，灰色为基质和低密度矿物斑晶，白色为高密度矿物斑晶； f—切片中的气孔分布图； g—批量处理CT切片的Matlab代码运行框架和操作流程。图2 基于工业CT扫描的三种气孔含量及大小分布统计方法Fig.2 Types of methods used to estimate content and size distribution of vesicles based on X-ray CT

玄武岩样品在中国石油化工股份公司石油物探技术研究院使用微焦点锥束CT实验设备(美国通用电气公司)进行测试，设备型号为Phenix v|tome|x s，测试电压150kV，束流220mA，检测样品的系统分辨率为19.9μm。典型的CT布局系统如图1b所示，X射线源和探测器分别置于转台两侧，转台带动样本360°旋转，每转动一个微小的角度后，探测器会接收由X射线照射样本获得的投影图数据。其原理是利用X射线在穿过物体的断面时，会因物质密度的不同而发生不同程度的衰减，计算机系统通过将衰减系数转为灰度值，再经过一系列校正后就可以重构样品的三维像素体图像[29]。因此，CT图像可以反映被检测物体内部物质密度的分布关系。对于具有较大气孔和简单气孔结构的气孔状玄武岩，由于气孔与基体密度差异大，CT图像上的灰度会明显不同(图2d)。一般气孔呈现灰度值较低的黑色，图像软件很容易借助灰度阈值分割，边缘检测，形态学运算(开闭、腐蚀膨胀等)等方法提取出岩石内部的气孔形貌，这一转化过程所带来的误差也可忽略[31-32]。

1.2 气孔含量及大小分布统计方法

在确保工业CT的X射线束空间分辨率足以解析火山岩气孔特征的前提下，使用配套商用软件VG Studio MAX(https://www.volumegraphics.com/en/products/vgstudio-max.html)重构样品的三维立体图像，并调节灰度阈值范围分离出气孔，测量其结构形态等参数(图2a～c)。然后，由该3D数据处理软件输出样品主视、侧视、俯视三个方向的二维CT切片，并通过开源软件ImageJ(http://rsb.info.nih.gov/ij/)对CT切片逐个进行图像处理和气孔形态学运算，从而在二维平面上更为精细地展现气孔细节，研究气孔在岩石内部的空间分布差异(图2d～f)。

尽管使用ImageJ处理切片可以获取所需的气孔含量、大小及形态数据，但这对于处理高精度CT扫描获得的大批量CT切片来说：①操作较为繁杂，耗时费力；②批量处理时占用电脑内存过大，不适用一般配置的个人电脑；③软件本身对单次处理文件也有一定的限制。为此，作者编写了一套专门处理CT切片数据的Matlab程序。本程序在Win10 X64位操作系统、Intel(R) Core(TM) i5-7200U CPU@2.5GHz2.71GHz、RAM 4GB的运行环境下，使用MATLAB R2018b批量处理900张1266×2204分辨率TIF格式的主视方向切片图仅耗时6.8min(其中CT Process为4.5min)，而且操作简便，占用内存小，可以随时中断然后继续运行(源代码见附件)。

该程序的运行框架和操作流程如图2g所示。首先，创建CT_image、binary_image、mask_image三个文件夹，并把CT_process.m、area_calculation.m、hole_stats.m三个Matlab代码置于同级目录下，然后将CT切片放在CT_image文件夹中，依次运行三个代码，最终输出每张切片中气孔的统计结果。其中，CT_process算法是该计算程序主要的图像处理部分。具体过程是：先将CT切片图转换为8位的二值化灰度图，在区分边界后，创建蒙版保护背景部分不被处理，再使用图像闭操作算法获取样本区域的大小，最后输出二值化气孔图和样本区域图，并依次循环处理下一张切片。在此基础上，area_calculation和hole_stats算法通过对气孔图和样本区域图中像素点的运算、统计，可以计算出感兴趣区的气孔含量和气孔大小分布情况。因为CT切片中灰度值的绝对大小还与原子序数和X射线能量相关，像元灰度值只能反映各物质密度的差异，并没有严格的物理意义[33]，所以在此之前需要基于VG Studio MAX和ImageJ的气孔含量测定结果(与之相等或相当)，采用逆分析的办法手动设定统一的灰度阈值分割点，以此减少因为阈值选取不恰当导致气孔提取偏少或部分基体也被提取出来的影响。

1.3 气孔体积-挥发份质量分数转换方法

大量研究表明，岩浆中的挥发份主要是二氧化碳和水蒸气，其含量和存在形式能够显著影响岩浆的流变学性质，在控制火山喷发行为中起着重要作用[34-36]。为定量分析挥发份在岩浆中的富集程度，将一定体积的气孔含量转化为火山岩中挥发份的质量分数是必要的。在假定溶解在熔体中的气体遵循“亨利定律”的前提下，Head等[37]将火山岩气孔的体积分数视为岩浆中逸出CO2或H2O的理想气体的含量，计算海底火山在临界爆破条件下不同水深的气体质量分数。但是，该方法采用的气体密度是通过理想气体状态方程得到的，随着水深加深压力增大，气体性质会逐渐偏离理想气体，最终得到的气体质量分数也会偏离实际情况。

已知平衡状态下一种气体的密度只与饱和温压有关，我们可以通过调用美国国家标准与技术局(NIST)热物性数据库中给定温压下的气体密度，对火山岩气孔体积-挥发份质量分数转换方法加以改进。该NIST热物性数据库(https://www.nist.gov/srd/refprop)使用的是当前可用的、最精确的气体状态方程和模型，收集的数据经过了严格的评估和审查，理论上会更为准确。如表1所示，本文对方法作了改进之后，转换结果与Head等[37]的方法相比，在同样的条件下，两者计算值的相对标准偏差(RSD)均在合理的范围内，在水深不超过3500m时，CO2含量小于4.30%，H2O含量小于0.95%，CO2偏差更大，并且随着水深加深压力增加，气体状态逐渐偏离理想气体的缘故，RSD有进一步扩大的趋势。因此，改进后的方法并非基于简单理想气体状态方程，而是采用数据库气体密度计算的方法得到的气体质量分数更加准确，特别是对于水深较大的海底火山岩样品。

2 结果与讨论

2.1 气孔的三维空间分布特征

在使用VG Studio MAX重构样品的立体空间后，可以直接获得气孔的个数、体积、表面积、直径等参数。通过该方法得到玄武岩样品的气孔体积分数为12.32%。如图3a所示，气孔尺寸分布呈现出对数正态分布的特点，相关系数R2皆大于0.96，这种分布可能反映岩浆中的气泡是由简单的成核和生长事件引起的，偏右分布的大气泡是许多小气泡发生聚集合并的结果[39]。其中，与气孔体积相等的等效球直径分布在100～800μm区间内，在180～200μm之间的数量最多；气孔最大外接圆直径(即最大直径)分布在140～1500μm区间内，分布区域相对较广，在340～360μm之间的数量最多。最大直径比等效球直径的分布明显偏右，这应该是气泡生长受到岩浆性质(黏度、流动性)影响导致气泡变形的结果[11,40]。如图3b所示，气泡从成核到不断长大的过程中，比表面积会不断减小，在达到一定尺度后趋于稳定，稳定值约0.02μm-1。这种变化特征除了反映气泡比表面积与直径近似呈反比例函数的关系，也在一定程度上反映了气泡生长向规则球形转变的趋势，因为气泡在应力松弛时本质上有利于岩浆-气泡界面上应力的均匀分布，从而减少其表面能[16,41]。

表1 本文与其他学者获得的气孔体积-质量分数转换结果对照

Table 1 Comparison of conversion results of vesicle volume to mass fraction between this paper and another scholar

水深(m)气体全部为CO2(%)气体全部为H2O(%)本文Head等(2003)[37]RSD(%)本文Head等(2003)[37]RSD(%)5001.9131.9350.570.8010.8020.0610003.6783.7611.121.5751.575015005.3435.5201.632.3352.3370.0420006.9147.2162.143.0843.0870.0525008.4018.8522.613.8213.8250.0530009.80910.4323.084.5454.5530.09350011.14512.1474.305.2585.3590.95

注：气体体积分数ηCO2/H2O(vol/%)=ηvesicle(vol/%)=75%，深水压强P(MPa)=1026(kg/m3)×9.8(m/s2)×Depth(采样水深,m)/106+0.101325(MPa,标准大气压)，岩浆温度T=1255℃，熔体密度ρlava=2700(kg/m3)；气体密度ρCO2/H2O(kg/m3)采用的数值不同，Head等[37]由理想气体状态方程计算得到，本文是先在系统中安装好Nist Refprop9.1，并将其加载到Excel宏中，然后在单元格中输入气体的密度函数[38]调用NIST数据库特定温压下的气体密度得到；岩石密度ρrock(kg/m3)=ρlava(kg/m3)×[1-ηCO2/H2O(vol/%)]+ρCO2/H2O(kg/m3)×ηCO2/H2O(vol/%)，气体质量分数ηCO2/H2O(wB/%)=ρCO2/H2O(kg/m3)×ηCO2/H2O(vol/%)/ρrock(kg/m3)×100%。RSD为两种方法CO2质量分数计算值的相对标准偏差。

a—气孔尺寸分布直方图，红线为体积等效球直径拟合曲线，绿线为最大外接圆直径拟合曲线；b—气孔尺寸与比表面积关系，比表面积为气孔表面积与体积之比。图3 气孔的三维空间分布情况Fig.3 Three-dimensional spatial distribution of vesicles

2.2 气孔的二维切片面分布特征

通过ImageJ和Matlab处理主视、左视、俯视三个方向上的二维切片，获得了气孔含量和气孔数密度在剖面上的变化特征(图4)。二维切片统计的气孔含量在各个剖面上围绕三维体积分数12.32%上下波动。两种方法在三个方向上所有切片气孔的平均含量分别为主视12.69%/12.48%、左视12.66%/12.50%、俯视12.66%/12.51%，与三维统计结果的相对标准偏差(RSD)均小于1.48%，并且两种方法统计的气孔含量的变化趋势和数值大小基本重合，同步性近乎一致。这说明基于切片图作数据再处理的方法在统计精度上是稳定可靠的。

从剖面方向来看，随着深度加深，气孔含量在主视、左视方向上的切片波动较大，但无明显分段，在俯视方向上的切片大体可以分为三段：10mm以下13%～16%，10～35mm之间11%～13%，35mm以上11%～14%，气孔含量整体呈现先减小后增加的趋势。而且，三个剖面上的气孔含量变化基本与气孔数密度(单位面积的气孔个数)的变化正向同步，表明气孔数密度是决定气孔含量大小的决定性因素。少数不同步的区间段(图4a，8～16mm；图4b，4～10mm)是因为局部切片中大气孔与小气孔占比差异较大的缘故，气孔数密度比较一致的情况下，大气孔多的气孔含量会有所升高。

a—主视图方向剖面；b—左视图方向剖面；c—俯视图方向剖面。

图4 二维切片中气孔含量和气孔数密度的面分布特征

Fig.4 Profile distribution features of vesicle content and vesicle number density in two-dimensional slices

火山岩中气孔的大小，形状和空间分布记录了岩浆脱气过程中气泡的形成和生长历史，最终的气泡尺寸分布取决于成核的早晚和次数，并受气泡生长、合并和变形的影响[3,14,39]，这些过程与岩浆性质、气体扩散速率、温压条件变化等因素密切相关[1,11,42]，所以该样品中气孔面分布特征(含量、数量、大小)的较大差异，表明气泡的成核和生长过程受到了多种因素的影响。

2.3 挥发份的质量分数

根据第1.3节改进后的气孔体积-挥发份质量分数转换方法，计算南海气孔状玄武岩样品的挥发份质量分数。将采样水深近似等于古水深Depth=1488m，岩浆温度和熔体密度根据全岩主量成分数据，使用Kurt Hollocher设计的CIPW norm Program计算(Fe3+/Fe=0.12)，其中T=1286℃，ρlava=2971kg/m3，气孔体积分数采用三维统计值η(vol/%)=12.32%，最终计算出气体的质量分数为ηCO2(wB/%)=0.233%，ηH2O(wB/%)=0.099%。若岩浆逸出的气体全部为CO2，挥发性气体过饱和点立刻成核，根据玄武质岩浆中CO2的溶解度公式，可以得到岩浆开始脱气释放出CO2，大致发生的最大深度为海底以下3829m(图5)。

岩浆的脱气过程与挥发份在岩浆中溶解度的变化密切相关，在恒定的温度下，挥发份的溶解度受压力和熔体成分的控制[34]。当岩浆沿着裂隙通道上升侵位或喷出地表时，由于玄武质岩浆中CO2含量高，溶解度远比H2O小(图5)，因此在海底高压(>10MPa)环境下，岩浆脱气进入气泡中的挥发份绝大部分是CO2[43-45]。尽管CO2的扩散系数相对较小，但仍可快速实现熔体溶解的CO2与气泡之间的平衡，所以对于一般情况下的熔体上升速率(未爆破)，岩浆能够保持近似平衡脱气的过程[44,46]。

溶解度公式：n=0.0023P(CO2)，n=0.1078P0.7(H2O)。式中：n(wB/%)代表溶解气体的量，压力P(MPa)=Depth(m)×9.8(m/s2)×2700(kg/m3)/106。

图5 气体在玄武质岩浆中的溶解度与压力的关系(据Parfitt等[47]修改)Fig.5 Relationship between gas solubility in basaltic magma and pressure (Modified according to Parfitt,et al[47])

2.4 方法的适用性和有待改进之处

基于高分辨率工业CT扫描，本文使用的三种统计气孔的软件程序各有其利弊。其中：①商业化软件VG Studio MAX可视化功能强大，可以进行样品的三维重建，但无法对内部细节逐个分析，功能也相对有限，且数据处理成本较高；②ImageJ是开源软件，免费使用，功能齐全，对图像格式限制少，但操作繁杂，效率较低，占用内存大，不适用于大量切片的批量处理；③Matlab代码的图像处理功能和气孔提取算法，可以快速获取岩石内部不同方向上各个剖面的气孔分布情况，批量处理，操作简便，占用内存小，针对个性化的需求，还可再度编译。

上述气孔统计方法仍存在一些局限性，如无法区分气孔中充填的杏仁体，在通过灰度阈值分割提取切片中的气孔时会存在视觉误差。这对于气孔尺寸大且结构简单的玄武岩样品来说影响可以忽略不计，但如果是微气孔含量高、气孔结构复杂的样品，最直接的解决办法是使用更高精度的工业CT进行测试，或从切片重构算法、三维数据逆处理、多尺度单元体表征三个方面对CT检测和数据处理过程进行改进[28]。对于岩石中气孔相互连接组成的多孔网络，则需要基于“骨骼”的概念设计算法将单个气孔分离、标记和计数，从而了解多孔空间的拓扑结构[48]。这些信息可以与岩浆的物理、流变和化学特性相结合，为模拟和理解岩浆过程提供基础。

火山岩气孔体积转化为挥发份质量分数的方法假定的是封闭岩浆体系，岩浆从洋壳深部到海底表层发生平衡脱气，气体全部为单一气体的理想情况。对于实际的岩浆脱气过程则需要考虑更多的因素；如果能选取具有代表性的多个火山岩(确保低中高气孔含量的样品都能采集到)，并将气体之间的比例、气泡定型时上覆岩层的压力、古水深的变化代入运算，那么由火山岩气孔得到的气体质量分数才可能代表真正意义上的岩浆脱气量。

3 结论

工业CT能够对岩石进行三维成像，可视化样品内部的空间结构，具有分辨率高和无损样品的特点。三种气孔统计方法能够快速提取岩石中的气孔，获取气孔的含量和大小分布情况。其中，商用软件VG Studio MAX适用于气孔的三维重构和体积测量，开源软件ImageJ适用于少量CT切片的图像处理和气孔的二维形态学运算，自主开发的程序代码使用方便快捷，适用于大量CT切片的批量处理。通过调用美国国家标准与技术局热物性数据库的气体密度，改进的海底表层环境下火山岩中气孔体积-挥发份质量分数转换方法，在深水条件下气体质量分数计算结果更加准确。

本文方法还可通过计算机算法分割岩石中相连的气孔和去除充填的杏仁体，并根据熔体包裹体成分和各气体的溶解度关系推测每种气体所占的比例来进行改进，从而反映更加真实的气孔构造特征和挥发份含量的变化，帮助人们了解火山喷发和岩浆脱气的动力学过程。