利用三维CT 扫描技术定量计算熔岩流气泡体积的研究与实现*

张茂亮 刘真 陈德峰 郭正府＊＊ 郭文峰

ZHANG MaoLiang1，LIU Zhen2，CHEN DeFeng3，GUO ZhengFu1＊＊ and GUO WenFeng1

1. 中国科学院地质与地球物理研究所新生代地质与环境重点实验室，北京 100029

2. 中国人民大学信息学院，北京 100872

3. 首都师范大学检测成像实验室，北京 100048

1. Key Laboratory of Cenozoic Geology and Environment，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing 100029，China

2. Information School，Renmin University of China，Beijing 100872，China

3. Computerized Tomography Lab，Capital Normal University，Beijing 100048，China

2014-02-08 收稿，2014-07-10 改回.

1 引言

新生代时期，青藏高原的隆升对亚洲内陆干旱化、东亚季风演化、全球气候变冷等区域乃至全球尺度上的气候环境演变均产生了重要的影响(An et al.，2001;Guo et al.，2002;Wu et al.，2012)，利用古高度计方法回溯大陆高原隆升历史是探讨其地球动力学机制及气候环境响应的重要手段之一。火山熔岩流气泡古高度计是国外近年来逐渐兴起的一种新的探讨古高度的研究方法(Sahagian and Proussevitch，2007)，其研究思路是:在测量火山熔岩流厚度、计算熔岩流气泡体积的基础上，定量计算熔岩流喷出时的古大气压强，再通过古大气压强换算成古高度。相比于其他古高度研究方法，火山熔岩流气泡古高度计的优势包括:(1)适用范围广，可应用于任何出露气孔状玄武岩的地区;(2)由于新鲜火山岩是开展同位素测年的理想材料，因此熔岩流气泡古高程计能够给出令人信服的隆升事件的年代，这为研究高原隆升历史提供了较可靠的年代学标尺;(3)研究对象独立于生物界之外，相对易于保存，更适用于高海拔、强烈侵蚀机制下的高原地区的古高度研究;(4)计算时，不需要假设古高度与温度、纬度之间的关系，计算结果的影响因素少(Sahagian and Maus，1994)。目前，国内科学家已成功地将火山熔岩流气泡古高度计方法应用到高原隆升研究中。例如，腾冲黑空山火山全新世熔岩流的古高度计算结果与目前实际海拔高度相吻合，表明火山熔岩流气泡古高度计将成为研究青藏高原隆升历史的有效手段之一(郭正府等，2011)。

熔岩流气泡体积的计算精度直接影响着古高度计算结果的误差，这是火山熔岩流气泡古高度计实际应用中的关键和难点。熔岩流气泡体积的定量计算方法包括注胶法、岩石抛光-扫描法、体视学转换法和三维CT (Computerized Tomography)扫描法(郭正府等，2011)。其中，测量和计算精度最高的手段是基于高分辨率X 射线层析成像技术和计算机算法的三维CT 扫描法。目前国内已建立了相关的三维CT 检测成像实验室，具备开展火山熔岩流气泡古高度计研究的实验测试条件，但熔岩流气泡的三维CT 扫描图像解译与气泡体积统计的计算机算法实现仍然是亟需建立和完善的重要环节。鉴于此，本文将提出新的基于三维CT 扫描技术的熔岩流气泡体积测量与统计算法，定量计算熔岩流气泡的数量与体积，并探讨该方法在火山熔岩流气泡古高度研究中的应用意义。

2 熔岩流气泡体积统计算法

熔岩流气泡的形态、数量、体积及分布特征对于探讨岩浆气泡的形成与演化、火山喷发动力学以及利用熔岩流古高度计重建高原隆升历史等均具有重要的作用。本文提出的基于三维CT 扫描技术的熔岩流气泡体积算法旨在从熔岩流CT 体数据中分割出独立气泡，并统计独立气泡的个数和体积，包括获取CT 体数据、数据分割、标记连通区域、骨架提取、计算候选喉面面积、筛选喉面、统计气泡体积等流程(图1)。

2.1 数据获取

计算机断层成像(Computerized Tomography，CT)是基于射线与物质相互作用的无损检测技术，其原理为:X 射线穿过物体时会发生强度衰减，而衰减程度与物体的密度、原子序数以及射线能量和初始强度有关。利用这个原理，让X 射线经不同角度穿过被检测物体并被探测器接收，利用算法重建被检测物体各部分的射线衰减系数以得到其分布图像，以此反映被检测物体的密度分布。目前，三维CT 扫描技术已被广泛地应用到工业、医学、生态学、考古学等多个领域(Zubal et al.，1994;Heeraman et al.，1997;Wu et al.，2008)，并已成为获取高分辨率三维岩石数据、探讨熔岩流古高度的有效手段之一(Song et al.，2001)。

与传统的二维扇束CT 相比，锥束CT 具有射线利用率高、扫描时间短、各向分辨率一致等优点。实际应用中，通常采用立式平板探测器系统作为锥束CT 扫描的工作平台(图2)。检测时，X 射线源和二维探测器都是固定的，待测物体在样品台上围绕Z 轴旋转，从不同角度接受X 射线的照射。本文所采用的三维岩石数据采集装置为首都师范大学检测成像工程中心研制的微焦点锥束CT 实验设备，射线源能量范围为10 ～225kV，系统最高空间分辨率为3μm。扫描的气孔状玄武岩岩芯样本直径为25mm，扫描电压能量为90kV。待测样品投影数据的分辨率为1024 ×1024 ×720，即转台每旋转0.5°采集一幅投影数据，其中每幅投影数据尺寸为1024 ×1024。

图1 熔岩流气泡体积统计的算法流程图Fig.1 Algorithm flowchart for vesicle volume calculation of lava flows

图2 锥束CT 扫描示意图Fig.2 Diagram of the cone-beam CT scanning

图3 气孔状玄武岩岩芯的三维CT 体数据(a)-熔岩流体数据三维视图;(b)-熔岩流体数据二维断层切片Fig.3 CT data of vesicular basaltic lava flows(a)-three-dimensional view of lava flow data;(b)-slice of lava flow data

FDK(Feldkamp-Davis-Kress)型算法是锥束CT 设备使用的主流重建算法之一(Feldkamp et al.，1984)。目前已在该套微焦点CT 设备上实现了基于GPU(Graphics Processing Unit)硬件加速的FDK 重建算法，包括滤波和反投影全过程(赵星等，2008)。重建的气孔状玄武岩三维CT 体数据的分辨率为512 ×512 ×512，其中每个体素的尺寸为49μm。图3a为重建后的气孔状玄武岩三维CT 体数据，图3b 为岩芯样品的二维断层切片体数据，三维CT 扫描图像清楚地展示了熔岩流在冷凝、固结成岩过程中形成的气泡及其形态、大小、分布等特点。

2.2 数据分割

图4 气孔状玄武岩三维CT 体数据的阈值分割结果Fig.4 Results of threshold segmentation on 3D CT data of vesicular basalt

在采集并获得气孔状玄武岩岩芯的三维CT 体数据后，利用不同物质对射线吸收衰减系数的不一致性，将体数据分割为具有明显灰度差异的岩石、气泡和无效(背景)部分，以进一步分析玄武岩中气泡的形态、大小并计算气泡体积。基于这样的三维CT 体数据特点，采用阈值分割方法处理气孔状玄武岩图像，其原理为:参照不同物质的灰度分布，选取恰当的阈值H0、H1。选定阈值后，遍历体数据中的所有体素，将所有体素的灰度值与阈值H0、H1相比较。所有灰度值小于灰度阈值H0的体素，将被判定为无效(背景)部分;相反，所有灰度值大于或等于灰度阈值H1的体素，将被判定为岩石部分;其余体素则为气泡部分(图4)。因此，选取合适的灰度阈值是准确判定气泡形态、大小与体积的关键(Song et al.，2001)。

2.3 标记连通区域

经过阈值分割的体数据，各气泡区域被分割成彼此分离的连通区域。所谓连通区域是指体数据在空间中的一个最大连通子集。连通性的定义如下:在一个连通集中，任意两个元素之间都存在一条完全由这个集合中的元素构成的路径，可以表示为:两个元素P、Q 是连通的，当且仅当存在一条路径P1-P2-P3-…-Pn，使得P1=P，Pn=Q，∀1≤i≤n-1 有Pi与Pi+1相邻。

标记连通区域就是寻找体数据中的所有气泡对象，并且将属于同一连通区域的所有体素用唯一的标记值进行标记。本文采用基于深度优先搜索原理的区域生长法对连通区域进行标记。首先对体数据进行顺序扫描，每遇到一个气泡体素就分配一个标记值，然后在该气泡体素的邻域进行检测。如果有尚未标记的气泡体素，则赋予相同的标记值，反复进行此操作，直到对体数据中所有体素扫描结束。标记后的连通区域可能是独立的单个气泡，也可能是融合的多个气泡(图5)。为了准确地提取气泡连通区域的骨架，还需要将融合的多个气泡分割开来，以提高气泡数量和体积的计算精度。

2.4 骨架提取

图5 连通区域标记结果示意图Fig.5 Tagging results of the connected regions

骨架是原始图形的一种压缩表示，与原始图形保持了相同的拓扑结构，并且存在于图形的对称轴上，能够反映图形的拓扑与形状信息。准确地对标记后的气泡连通区域进行三维骨架提取是生成和处理喉面的关键步骤。骨架提取算法通常可分为:拓扑细化法、距离变换法、几何分析法和广义势场法等(Cornea et al.，2007)。其中，拓扑细化法的优点在于可以较好保留原始图形的拓扑特征，即保证骨架连通性，同时也能保持骨架的单像素宽度。

本文采用Palágyi and Kuba(1999)提出的拓扑细化法对气孔状玄武岩三维CT 图像中的标记连通区域进行骨架提取。该算法的基本思想是从边界开始向内演化，利用反复迭代方法对三维模型进行层层剥离。在每次迭代过程中，将边界上的每个体素点及其26 邻域的体素与3 ×3 ×3 模板进行匹配，符合模板匹配的体素将被保留至下一次迭代过程，而不符合模板匹配的体素将被删除(图6)。以此类推，算法将逐步搜索到中轴骨架的位置，直至没有任何点可以被删除，最终获得单像素宽的骨架(图7)。

该算法是基于边界的算法，只对边界上的体素进行模板匹配，极大地减少了计算量。由于算法的并行性，可在一次迭代过程中从边界上删除若干点，因而适合使用GPU 等并行处理器硬件提高运行速率。

2.5 喉面处理

气泡连通区域的骨架提取完毕后，计算各个骨架点在骨架线上的切线方向，也就是候选喉面的法平面方向。连通区域的骨架可以视为一条空间中的曲线，利用插值型的求导公式来计算骨架线在该点的切线方向。令x0，x1，x2，x3，x4为骨架上连续的5 个节点，则对应每一点的求导公式为:

图6 骨架提取算法匹配模板Fig.6 Matching template for skeleton extraction algorithm

图7 连通区域骨架提取示意图Fig.7 Skeleton extraction of connected regions

对于连通区域骨架线长度小于5 的情况，则利用两点求导公式或三点求导公式(此处不再赘述)。切线方向确定后，也就能确定过该点且垂直于切线向量的法平面，称为候选喉面。计算候选喉面面积是统计玄武岩气泡体积的重要环节，本文提出一种新的方法——逐点法，用于候选喉面的面积计算。与其他计算方法相比，该方法的计算结果更为精确，且不受连通区域形状的约束。对于计算得到的所有候选喉面面积，需从中筛选出真正的喉面用于气泡体积计算，即确定融合气泡的分割面。本文采用基于区域划分的阈值计算方法，减少了传统单一阈值分割方法导致的小连通区域过度分割的误差，提高了气泡数量和体积的统计精度。

2.6 统计结果

当喉面确定后，即可对三维CT 体数据进行重新标记。在连通区域标记值的基础上，将独立气泡从连通区域中分离出来，并赋予新的标记值。当所有气泡均被标记后，则可统计独立气泡的数量与体积。不同的标记值数量代表独立气泡的数量;被标记为同一标记值的体素则归属于同一个气泡，因此可以统计出该气泡的体素个数。最后，将数据获取环节计算得到的气泡分辨率物理尺寸换算为每个独立气泡对应的实际物理体积，对气泡体积进行排序并分析其众数分布特征，计算得到气孔状玄武岩岩芯的气泡体积。

3 算法的主要优势

在本文提出的算法中，影响气泡体积最终计算结果的关键步骤包括:(1)分割三维CT 体数据，以判定三维CT 扫描图像中的岩石、气泡和无效(背景)部分;(2)提取气泡连通区域的骨架，以此构建气泡连通区域中单像素宽的连通线段;(3)计算候选喉面面积，即以骨架切线为法向量的截面面积;(4)筛选有效喉面，即从候选喉面中选取能够代表气孔状玄武岩气泡形态的分割面。其中，准确选取气泡连通区域分割面的操作难度最大，因此本文主要针对喉面面积的计算过程与有效喉面筛选方法进行了优化。

3.1 候选喉面面积计算

作为选取气泡连通区域分割面的基础，候选喉面面积计算是整个算法中最为重要的环节。多三角形近似法是目前应用较多的计算喉面面积的方法之一(Shin et al.，2005)，其基本思想是以气泡连通区骨架上的每个体素为顶点，在以该点切线方向为法向量的平面上，作多个固定角度的射线。以相邻两条射线与岩芯边界的交汇点作为两个顶点，与该骨架点构成了一个三角形(图8)。计算每个连通区域的单像素宽度的连通骨架线e，对于每条骨架线e 上的点Ve，都用多三角形近似法计算该点的法平面面积。根据三个顶点的空间坐标计算出单个三角形的面积，所有三角形的面积之和即为该截面的面积。

在单个骨架点法平面的截面中(图9)，以骨架点为中心、45°为固定间隔角度向岩芯作射线。当射线与岩芯边界相交时，就得到一个截面边界点Gj。以此类推，最终得到8个截面边界点。利用骨架点与8 个截面边界点组成的8 个三角形，计算所有三角形的面积之和，并将其作为候选喉面面积的近似值。多三角形近似法的精度可以通过不断增加三角形的数量来逐步逼近截面的面积，所选则的角度越小，则得到的三角形越多，那么三角形面积之和也就越近似于候选喉面面积的实际值。

图8 熔岩流候选喉面面积计算过程示意图Fig.8 Calculation process of candidate throat surface in lava flow

图9 多三角形近似法示意图Fig.9 Diagram of the multi-triangular approximation

多三角形近似法适用于候选喉面轮廓比较规则的情况，但实际上熔岩流气泡连通区域的截面往往是很不规则的(图10)，若使用多三角形近似法计算候选喉面面积，则会产生较大的误差。多三角形近似法的误差主要由两部分组成:(1)部分三角形计算结果偏大，对于气泡连通区域内凹的情形，统计时会将该部分岩石计入候选喉面面积;(2)部分三角形计算结果偏小，对于连通区域外凸的情形，统计时会忽略凸出的连通区域部分，不计入候选喉面面积。例如，图10 中左侧闭合部分为气孔状玄武岩的气泡连通区域，闭合区域以外的部分为玄武岩。当采用多三角形近似法计算气泡连通区域面积时，得到的计算结果(即斜线阴影部分)将会包括多算的部分和少算的部分，造成难以估计的误差(图10)。

图10 采用多三角形近似法计算截面面积导致的误差示意图Fig.10 Errors of cross section area resulting from multitriangular approximation

图11 逐点法示意图Fig.11 Diagram of traversal-points algorithm

鉴于此，我们提出一种新的计算候选喉面面积计算的方法——逐点法，以解决应用多三角形近似法计算不规则气泡连通区域截面面积时产生的不确定性误差问题。首先，我们指定一个阈值Hpointonplane作为筛选有效体素点的判断标准。然后，计算体素点到平面的距离，并与指定的阈值比较。

点到平面的距离公式:

其中，(a，b，c)为截平面的法向量，(x，y，z)为三维空间中一个点的坐标。

当点到平面的距离小于或者是等于指定的阈值时，认为该点在平面上，即被判定为有效体素点。如图11 所示，实心点代表落在平面上的有效体素点，而空心点代表非平面上的体素点。候选喉面面积越大，落在该平面上点的数目也会越多;反之亦然。

利用逐点法计算候选喉面面积的基本步骤为:

记Px为编号为x 的所有标记为气泡的体素集合;Ex为编号为x 的连通区域的骨架线;Vex，j为Ex上的第j 个骨架点;Vpx，j为落在第x 区域第j 个点的截面的体素点集合;编号Ax，j为x 区骨架上第j 个点的截面面积。气泡体素v∈Px，且Dv，Vex，j≤Hpointonplane，则v∈Vpx，j。将落在平面上的有效体素个数来作为候选喉面面积，则Ax，j= |Vpx，j|(表1)。

表1 应用逐点法计算候选喉面面积的流程Table 1 Procedure of traversal-points algorithm for calculation of candidate throat surface area

虽然在应用多三角形近似法计算候选喉面面积过程中，可以通过减小试探角度来提高精度，但如何选取试探角度以及确定射线与岩芯边界的交点都是较难解决的问题，需要通过多次实验来调整合适参数以得到理想的计算结果。此外，由于不同成因熔岩流的气泡形态、截面、分布特点等差别悬殊，若要通过多三角形近似法获得较高精度的熔岩流气泡体积，只能针对不同的熔岩流分别调整参数，因此多三角形近似法的适用性较差。

相较而言，逐点法减少了不规则气泡连通区域形状对候选喉面面积计算的干扰。更重要的是，影响逐点法计算结果的因素只有一个，即像素点到平面的距离阈值，而距离阈值的设定与气泡截面的变化无关，使得逐点法可被方便地应用于不同成因的熔岩流气泡体积计算，其适用范围比多三角形近似法更广泛。此外，逐点法易于利用计算机算法实现，可以方便地移植到GPU 等并行计算平台上。因此，在实际应用中逐点法比多三角形近似法更具有优势。

3.2 喉面选取

喉面选取的原则是从骨架线上的截面中选择最有可能成为气泡之间分割面的平面。通用的作法是为所有的平面选择一个全局阈值，小于阈值截面面积的平面将被视为分割面。然而，这种阈值选取方式忽略了一个重要的问题:不同连通区域的截面面积存在较大差异，而全局阈值方法可能导致气泡被过度分割。例如，图12 中标记值为9 的连通区域的宽度明显大于其它区域。如果对所有连通区域都采用同一个全局阈值，将会使一些体积较小的连通区域(如标记值为1、2、6、8 的区域)分割过碎，导致无法准确地选择合适的喉面。解决过度分割问题的思路就是对不同宽度的气泡连通区域采用不同的喉面阈值。因此，在本文中我们提出将阈值分段应用于不同连通区域的方法——基于区域划分的阈值计算方法。相对于传统的仅采用单一阈值的判断方法，基于区域划分的阈值计算方法自适应地调整了不同宽度连通区域的阈值，提高了分割精度。

图12 气泡连通区域的过度分割示意图Fig. 12 Diagram of excessive separation on connected regions

首先，设置算法阈值Dneck(0≤Dneck≤1)。将不同连通区域骨架线上点的截平面划分为不同的组，并对不同分组的候选喉面应用阈值Dneck。通过这种设定，体积较大的连通区域所选用的喉面阈值就会相应大一些，而体积较小的连通区域将会得到较小的阈值，从而缓解了仅采用全局单一阈值所产生的较小连通区域被过分分割的问题。然后，将Dax定义为标记值x 对应的气泡连通区域的喉面阈值，其实质为不同连通区域的骨架截面积之和乘以一个人工设定的阈值Dneck，计算公式如下:

如果∀a∈Ax且a≤Dax(a ＞0)，则此时的截面a 即为有效的气泡连通区域分割面。

4 程序运算及结果

统计程序采用C+ +语言编写，包括8000 余行代码。通过程序运算得到气孔状玄武岩岩芯的气泡数量及体积，然后将气泡体积进行排序并分析其众数分布特征(图13)。

在生成众数统计结果时，参照Sahagian et al. (2002)对横坐标取对数的方法。为了更好地分析统计结果，对于取对数后的横坐标乘以3，即拉伸3 倍以观察气泡体积众数分布。从图13 中可以得知，独立气泡体积Ln (volume)×3 的众数分布区间为［6.5，7.5］，则独立气泡体积Volume 众数分布区间为［9，12］(单位:体素)。在熔岩流气泡众数统计软件计算结果的基础上，利用火山熔岩流气泡古高度计原理，郭正府等(2011)恢复了腾冲黑空山火山全新世喷发形成的气孔状玄武质熔岩流的古高度，结果表明基于三维CT 扫描技术的熔岩流气泡体积定量计算方法是有效和可行的。

图13 气泡众数统计结果Fig.13 Results of vesicle volume mode

在算法运算效率方面，对于512 ×512 ×512 分辨率的体数据，目前算法的计算耗时是约为6 小时(测试平台为CPU:Intel Core2 Duo P8700@2.53GHz，内存:4GB，硬盘320GB@5400RPM)。

5 结论

本文提出了基于计算机断层成像技术的熔岩流气泡众数统计新方法，并针对影响熔岩流气泡数量和体积统计精度的喉面面积计算与喉面选择两个关键技术环节提出了改进方法，结果显示:

(1)利用逐点法计算熔岩流体素点落在平面上的数量，以统计截面面积大小，能够有效地克服多三角形近似法对于截面边界形状敏感的缺点，并提高喉面面积计算的精度;

(2)本文提出的基于区域划分的阈值计算方法，解决了传统喉面选择方法中单一阈值所产生的小连通区域过度分割的问题，减小了气泡数量统计的误差。

根据上述新方法，我们设计并实现了熔岩流气泡众数统计软件，其有效性和可行性得到了气孔状玄武岩岩芯样品测试结果的验证。该方法有望提高熔岩流气泡数量、体积以及火山熔岩流气泡古高度计的计算精度。在后续工作中，需要通过进一步的实验调整系统参数，减小算法误差，以达到最佳的计算精度。在算法性能优化方面，可利用并行计算硬件进行算法加速(如多核CPU、GPU 等)，以提高统计运算的效率。通过软件与硬件两方面的进一步优化，使熔岩流气泡众数统计方法成为高原隆升研究的快速分析工具。

致谢 熔岩流气泡体积统计与计算方法由第二作者刘真完成;气孔状玄武岩三维CT 扫描测试与图像处理由第三作者陈德峰完成;熔岩流气泡测量与统计结果的地质解译由通讯作者和第一作者完成;样品测试过程中得到首都师范大学张朋教授的热心帮助;在此一并表示感谢。

An ZS，Kutzbach JE，Prell WL and Porter SC. 2001. Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan Plateau since Late Miocene times. Nature，411(6833):62 -66

Cornea ND，Silver D and Min P. 2007. Curve-skeleton properties，applications，and algorithms. IEEE Transactions on Visualization and Computer Graphics，13(3):530 -548

Feldkamp LA，Davis LC and Kress JW. 1984. Practical cone-beam algorithm. Journal of the Optical Society of America A，1(6):612-619

Guo ZF，Zhang ML，Cheng ZH，Liu JQ，Zhang LH and Li XH. 2011. A link of measurements of lava flows to Palaeoelevation estimations and its application in Tengchong volcanic eruptive field in Yunnan Province (SW China). Acta Petrologica Sinica，27(10):2863 -2872 (in Chinese with English abstract)

Guo ZT，Ruddiman WF，Hao QZ，Wu HB，Qiao YS，Zhu RX，Peng SZ，Wei JJ，Yuan BY and Liu TS. 2002. Onset of Asian desertification by 22Myr ago inferred from loess deposits in China.Nature，416(6877):159 -163

Heeraman DA， Hopmans JW and Clausnitzer V. 1997. Three dimensional imaging of plant roots in situ with X-ray computed tomography. Plant and Soil，189(2):167 -179

Palágyi K and Kuba A. 1999. A parallel 3D 12-Subiteration thinning algorithm. Graphical Models and Image Processing，61(4):199-221

Sahagian DL and Maus JE. 1994. Basalt vesicularity as a measure of atmospheric pressure and palaeoelevation. Nature，372(6505):449-451

Sahagian DL，Proussevitch AA and Carlson WD. 2002. Analysis of vesicular basalts and lava emplacement processes for application as a paleobarometer/paleoaltimeter. The Journal of Geology，110(6):671 -685

Sahagian DL and Proussevitch AA. 2007. Paleoelevation measurement on the basis of vesicular basalts. Reviews in Mineralogy ＆Geochemistry，66(1):195 -213

Shin H，Lindquist WB，Sahagian DL and Song SR. 2005. Analysis of the vesicular structure of basalts. Computers ＆ Geosciences，31(4):473 -487

Song SR，Jones KW，Lindquist WB，Dowd BA and Sahagian DL. 2001.Synchrotron X-ray computed microtomography: Studies on vesiculated basaltic rocks. Bulletin of Volcanology，63(4):252-263

Wu GX，Liu YM，He B，Bao Q，Duan AM and Jin FF. 2012. Thermal controls on the Asian summer monsoon. Scientific Reports，2:Article Number:404，doi:10.1038/srep00404

Wu XJ，Liu W，Dong W，Que JM and Wang YF. 2008. The brain morphology of Homo Liujiang cranium fossil by three-dimensional computed tomography. Chinese Science Bulletin，53(16):2513-2519

Zhao X，Zhang HT，Chen M et al. 2008. Progress in CT research at the computer tomography laboratory of Capital Normal University.Chinese Journal of Stereology and Image Analysis，13(3):158 -165 (in Chinese with English abstract)

Zubal IG，Harrell CR，Smith EO，Rattner Z，Gindi G and Hoffer PB.1994. Computerized three-dimensional segmented human anatomy.Medical Physics，21:299 -302

附中文参考文献

郭正府，张茂亮，成智慧，刘嘉麒，张丽红，李晓惠. 2011. 火山“熔岩流气泡古高度计”及其在云南腾冲火山区的应用. 岩石学报，27(10):2863 -2872

赵星，张慧滔，陈明等. 2008. 首都师范大学检测成像实验室CT 研究进展. 中国体视学与图像分析，13(3):158 -165