基于CT影像的砂岩文物结构特征分析

王莹莹，徐金明, ，黄继忠

(1.上海大学土木工程系, 上海 200444；2.上海大学文化遗产保护基础科学研究院, 上海 200444)

0 引言

石窟文物大多有数百甚至上千年的历史，在水和环境污染等各种营力作用下，普遍存在风化问题、裂隙问题、水岩相互作用等问题[1]。由于文物表面变化过程取决于表面的成分和结构，使用CT(Computerized Tomography，计算机层析摄影)扫描试验和数字图像处理技术可以确定表面附近的成分、结构及其变化特点[2]。

CT技术具有无扰动、可多层面分析、机制明确等优点，近年来受到岩土工程研究者的关注。方建银等[3]将岩石CT扫描断面分为完整区、损伤区和孔洞裂纹区，分析了三轴压缩试验过程中粉砂岩孔隙率、损伤率和完整率的变化特征；尹小涛等[4]通过提取CT 图像中的裂纹，分析了加载过程中裂纹分维数的变化特征；薛华庆等[5]使用CT扫描结果来表征油砂、致密砂岩和页岩样品的微观结构，分析了常规测试方法与CT扫描表征结果的差异性；HAN等[6]基于CT扫描图，建立了合成岩石和贝雷砂岩的数值模型，利用试验结果验证了数值模型的合理性；闻磊等[7]建立了利用CT数和X射线吸收系数表示的岩石冻融衰变函数模型，分析了岩石冻融过程中的完整性损失。

一些学者结合数字图像处理技术与CT图像研究岩石的结构特征。朱健伟等[8]结合数字图像处理技术和扩展有限元方法研究了变形破坏过程中岩石试样中的应力场、位移场、不同组分损伤情况及其与裂缝萌生扩展过程的关系。刘向君等[9]以常规砂岩为研究对象，使用微CT扫描结果与图像处理技术建立了反映真实孔隙结构的三维数字岩芯模型；于庆磊等[10]利用数字图像处理技术和CT图像表征混凝土材料的结构，提出了环状分区与分割阈值自动识别相结合的CT图像分割算法；杨振琦等[11]利用数字图像处理技术和CT图像表征花岗岩试件的结构，建立了花岗岩试件的三维数值模型，研究了对单轴压缩条件下花岗岩试件的破裂过程。

现有研究中，砂岩中不同组分实际类型和真实分布的重视还不够，也没有将CT技术与砂岩表面附近风化特征很好结合起来。本文拟在使用CT技术获取三处砂岩文物(石门山、云冈石窟、乐山大佛)表面附近切片图像基础上，分析不同深度切片图像的成分、结构与数字特征参数，研究图像特征参数、砂岩所处环境、表面附近风化程度变化的关系，这对合理分析石质文物特点、有效制订石质文物保护措施都有一定的参考价值。

1 CT试验原理

CT 扫描试验所用仪器由X射线源、样品转台、探测器、重构系统组成。X射线穿透物质后，X 射线强度衰减程度与物体密度、物体厚度成正比；CT 扫描图片的灰度值直接表现了X 射线的衰减程度：物质的密度越大，X射线衰减程度就越高，扫描图片灰度值就越大。对于岩石孔隙和裂隙，由于主要成分为密度很小的空气，X 射线衰减程度低，相应CT图像为黑色、灰度值较小。

试验仪器主要性能指标为：X 射线源电压为140 kV；电流为100 μA；扫描厚度为41.5 μm；放大倍数为4.81；空间分辨率为1～2 μm；密度对比分辨率为20%；可识别最小体积为10 μm左右。主要操作参数为：扫描模式为连续扫描；扫描过程中，X 射线源和探测器位置保持不变，样品匀速从180°旋转至180°、每旋转n度拍摄1张图片。

2 砂岩文物中不同组分类型和分布的确定

砂岩文物的风化，主要取决于岩石结构。岩石结构主要包括岩石中不同组分的颗粒形态、大小和连接方式。下面将根据所得CT图像、使用数字图像处理技术研究不同组分的类型和分布，它们是岩石结构的核心要素。

2.1 砂岩文物的来源

试验样品分别来源于云冈石窟、乐山大佛和石门山石窟的砂岩，它们的外观情况见图1。

图1 不同样本的外观情况Fig.1 Appearances of specimens at various locations

2.2 砂岩文物的CT图像

对三种不同砂岩文物试样进行了CT扫描试验，某一切片位置的CT 扫描图像如图2 所示，为后续处理方便，将图像都转化为灰度图像。灰度图像中，灰度级是0～255中的整数。图2中，云冈石窟砂岩切片位置距表面14.822 mm，像素高度最大为477，像素宽度最大为585；乐山大佛砂岩切片位置距表面16.102 mm，像素高度最大为477，像素宽度最大为618；石门山砂岩切片位置距表面14.822 mm处，像素高度最大为411，像素宽度最大为963。

图2 不同地区砂岩文物样本的CT图像Fig.2 CT images of sandstone relics specimens in various regions

2.3 CT图像中组分类型和分布的确定

本次研究采用灰度分界阈值法、通过编制MATLAB代码来获得砂岩文物CT图像中不同组分的类型和位置。

(1)不同组分灰度分界阈值的确定

图2中砂岩文物的主要组分类型是石英和长石。这两种组分的灰度分界阈值，使用点选法确定：首先采用肉眼鉴定确定不同位置点对应的成分类型，然后提取该点的位置坐标x、y和相应的灰度值；然后获得所选点灰度值的统计结果；最后将所得统计结果应用于整幅图像[12]。

对于图2(a)，不同组分的点选位置与灰度值见表1，相应的统计结果见表2。为了将点选灰度值统计结果更好地应用于整幅图像，使用灰度值统计结果(灰度均值和方差)来确定两种组分的灰度分界阈值T。

假设两种组分灰度为相邻分布，灰度较高组分的灰度均值和方差分别是H1和σ1，灰度较低组分的灰度均值和方差分别是H2和σ2，将两种组分的灰度分界阈值T定义为：

表1 图2(a)中不同组分的位置与灰度值Table 1 Location and gray value of various compositions in Fig.2(a)

表2 表1中不同组分的统计结果Table 2 Statistics of various compositions from table 1

(1)

式中：int——取整计算。

H1——长石点选点数；

H2——石英点选点数；

σ1——长石灰度均值；

σ2——石英灰度均值，取整是因为灰度值必须为一整数。

使用这一方法，使阈值分割结果与肉眼观察所得不同组分总体百分比尽量一致，得到图2(a)(云冈石窟砂岩)长石和石英之间的灰度分界阈值为25。使用相同方法得到图2(b)(乐山大佛砂岩)长石和石英的灰度分界阈值为45，图2(c)(石门山砂岩)长石和石英的灰度分界阈值为27。

(2)不同组分的分布

使用前述所得灰度分界阈值分别对图2(a)、图2(b)、图2(c)进行分割，就可得到砂岩文物中不同组分分布的二值图像(分别见图3、图4和图5)。

由图3～图5可知，云冈石窟砂岩(图3)中长石与石英含量相差不大，乐山大佛砂岩(图4)中长石含量远多于石英含量，而石门山砂岩(图5)中石英含量远多于长石含量。

图3 由图2(a)得到的不同组分分布Fig.3 Composition distributions from CT image Fig.2(a)

图4 由图2(b)得到的不同组分分布Fig.4 Composition distributions from CT image Fig.2(b)

图5 由图2(c)得到的不同组分分布Fig.5 Composition distributions from CT image Fig.2(c)

3 不同砂岩文物结构特征的比较

3.1 不同地区砂岩文物的环境特征风化作用是砂岩文物的主要病害之一，取决于岩石中不同组分类型和环境因素[13]。不同地区砂岩文物的主要环境特点如表3所示。由表3可以看出，不同地区砂岩的环境特点差异很大。

3.2 砂岩组分的几何特征参数

砂岩组分的几何特征参数可使用圆形度R0、矩形度R、内切圆半径r(mm)、形状离散指标e及偏心率E来表征[14]，分别表示组分区域接近圆的程度、对外接矩形的充满程度、内切圆的半径、形状的复杂性及对象区域的紧凑性，计算公式分别为：

表3 不同地区砂岩的环境特点Table 3 Environmental features of sandstone relics in various regions

R0=4πS/L2

(2)

R=S/SM

(3)

r=2S/L

(4)

e=L2/S

(5)

E=A/B

(6)

式中：S——组分区域的面积/mm2；

L——组分区域的周长/mm；

SM——组分区域最小外接矩形的面积/mm2；

A、B——分别表示组分区域的长轴长度与短轴长度/mm。

使用这些计算公式，得到不同地区砂岩文物中不同组分几何特征参数的计算结果(表4)。

由表4可以看出，就颗粒面积和周长来说，都是石门山砂岩最大，乐山大佛砂岩次之，云冈石窟砂岩最小；对于颗粒圆形度，所有三处文物长石圆形度都远小于石英圆形度；对于矩形度，乐山大佛砂岩较大，石门山砂岩较小，云冈石窟砂岩比较离散；对于内切圆半径来说，石英都远大于长石；就离散指标来说，石英都远小于长石；对于偏心率来说，乐山大佛砂岩中长石远大于石英，且云冈石窟砂岩的偏心率小于石门山砂岩的偏心率。

这些几何特征参数变化与砂岩文物环境密切相关。云冈石窟与乐山大佛都临水而建，乐山大佛更是濒临大渡河、青衣江和岷江三江汇流处，因此，乐山大佛砂岩受水的影响最大、颗粒形状变化最小，云冈石窟砂岩次之，石门山颗粒形状最不规则。

表4 不同地区砂岩文物中不同组分的几何特征参数Table 4 Geometrical parameters of compositions in sandstone relics at various regions

3.3 基于灰度直方图的数字图像特征参数

基于灰度直方图的纹理特征参数，也可用于描述颗粒组分区域的形状特征。这些参数包括均值μ、方差σ2、偏度S、峰值K、能量E和熵H，计算公式分别为：

(7)

式中：p(ri)——灰度值ri出现的概率；

n(ri)——图像中灰度值为ri的像素点总数；

N——图像的总像素点数。

均值表示一幅图像或区域的平均亮度；方差反映了纹理的粗糙程度；偏度是对灰度分布偏离对称情况的度量；峰值描述了灰度分布的倾向；能量是图像所具有的信息量的度量；熵反映了图像中纹理的非均匀程度或复杂程度。

不同样本的灰度直方图(横坐标是灰度级，纵坐标是灰度级出现的概率)见图6。

图6 不同地区样本的灰度直方图Fig.6 Gray histogram of specimens in different regions

三处砂岩文物基于图6(灰度直方图)的图像纹理特征参数计算结果如表5所示。

由表5可以看出，均值是云冈石窟砂岩最大、石门山砂岩次之、乐山大佛砂岩最小，表明云冈石窟砂岩在灰度直方图上平均亮度最大；方差是云冈石窟砂岩最大、石门山砂岩次之、乐山大佛砂岩最小，表明云冈石窟砂岩的纹理最为粗糙；偏度是乐山大佛砂岩最大、石门山砂岩次之、云冈石窟砂岩最小，表明云冈石窟对灰度分布偏离对称情况最小；峰值是乐山大佛砂岩大于石门山砂岩、云冈石窟砂岩，表明乐山大佛砂岩灰度分布散布于尾端；能量是乐山大佛砂岩最大、云冈石窟砂岩次之、石门山砂岩最小；熵是石门山砂岩最大、云冈石窟砂岩次之、乐山大佛砂岩最小，表明石门山砂岩纹理最复杂、乐山大佛砂岩最为简单。

表5 不同地区砂岩CT图像基于灰度直方图的图像纹理特征参数Table 5 Textural parameters based on gray histogram of sandstone CT images in different regions

基于灰度直方图的这些CT图像纹理特征参数也与砂岩文物环境密切相关。云冈石窟的砂岩为杂砂岩，分选性差、磨圆较差、结构成熟度低，一般富含石英、含量可达50%；乐山大佛两侧岩石是红砂岩，质地疏松、易于风化；石门山石窟多雕刻于含泥质砂岩，砂岩孔隙大、吸水性强。因此，云冈石窟砂岩在灰度直方图上平均亮度最大、纹理最粗糙、灰度分布偏离对称情况最小；乐山大佛砂岩灰度分布散布于尾端；石门山砂岩纹理最为复杂，乐山大佛砂岩最简单均匀。

3.4 基于灰度共生矩阵的数字图像特征参数

灰度直方图不能完整反映图像像素的相对位置，可以根据灰度共生矩阵、通过计算对比度f1、相关度f2、角二阶矩f3、逆差矩f5和熵f4等纹理参数来描述图像灰度值的空间位置关系，图像对比度反映了图像的清晰程度和纹理的深浅程度；相关度反映了图像纹理的精致程度；角二阶矩反映了灰度分布均匀程度和纹理粗细程度；逆差炬反映了图像纹理局部变化的大小。

计算公式分别为

(8)

式中：i，j——分别为两个像素的灰度；

p(i，j)——共生矩阵归一化后的元素。

不同地区砂岩文物CT图像中不同组分基于灰度共生矩阵图像纹理特征参数的计算结果见表6。

表6 不同地区砂岩CT图像基于灰度共生矩阵的图像纹理特征参数Table 6 Textural parameters based on grayscale co-occurrence matrix of sandstone CT images in various regions

由表6可以看出，三个地区砂岩文物中，长石的对比度都小于石英，说明石英图像的纹理和沟纹特征都大于长石；云冈石窟砂岩的对比度>乐山大佛砂岩的对比度>石门山砂岩的对比度，说明云冈石窟砂岩的纹理和沟纹特征最大，乐山大佛砂岩居中，石门山砂岩最小；砂岩中长石和石英的相关度都相差不大，说明砂岩图像纹理的精致程度相似；砂岩的角二阶矩都为0.514 0～0.641 2，说明砂岩中不同组分的纹理一般、不够细致均匀；砂岩中长石的熵都等于石英的熵，云冈石窟砂岩的熵>石门山砂岩的熵>乐山大佛砂岩的熵；砂岩中长石和石英的逆差矩，都是石门山砂岩>乐山大佛砂岩>云冈石窟砂岩，说明石门山砂岩不同组分变化规律性最强，乐山大佛砂岩居中，云冈石窟砂岩最杂乱。

基于灰度共生矩阵的这些CT图像纹理特征参数与砂岩文物环境密切相关。云冈石窟建造年代最为久远，经受了较多的人为与生物破坏；石门山建造年代距今最近；乐山大佛和石门山石窟气候湿润，风力较小，常年处于髙湿状态，水岩作用显著。因此，云冈石窟砂岩的纹理和沟纹特征参数最大，乐山大佛砂岩居中，石门山砂岩最小；三个地区砂岩图像纹理的精致程度相似，不同组分的纹理不够均匀；石门山砂岩图像区域的变化规律最强，乐山大佛砂岩居中，云冈石窟砂岩最差。

3.5 砂岩文物纹理参数随深度的变化情况

随着CT图像扫描深度的变化，不同地区岩样纹理特征也会发生变化[15]。对所有CT图像，每隔一定距离(乐山大佛、云冈石窟、石门山砂岩分别为415、382、382 μm)取一幅CT图像计算基于灰度直方图的CT图像纹理特征参数，得到纹理参数随深度的变化情况(图7)。图中，UZ1和UZ2分别为均值和一致性指标。

图7 不同地区砂岩文物基于灰度直方图纹理参数随深度的变化Fig.7 Textural parameters of specimen CT images from gray histogram at various depths in three regions

由图7可知，对于云冈石窟砂岩，距表面25 mm内参数变化比较平稳，但25 mm后纹理参数大小明显增加；对于乐山大佛砂岩来说，距表面7.5 mm内变化比较明显，但7.5 mm后趋于稳定；对于石门山砂岩来说，距表面7.5 mm内纹理参数急剧增加，7.5 mm到30 mm纹理参数急剧下降，30 mm后纹理参数又急剧上升，变化比较明显。

因此，随着距表面距离大小的增加，乐山大佛砂岩纹理变化最小，云冈石窟砂岩次之，石门山砂岩形状最不规则，三处砂岩纹理复杂大小的顺序是石门山石窟>云冈石窟>乐山大佛。这说明，乐山大佛砂岩表面风化严重，随深度增加风化程度变化不明显；云冈石窟砂岩是随着深度的增加，风化起初不明显，但深度达到25 mm后，风化程度骤增；石门山砂岩由于各种营力作用与岩性的复杂性，随深度增加、风化程度变化比较复杂。

4 结论与展望

本文在获取石门山石窟、云冈石窟、乐山大佛砂岩表面附近CT图像基础上，分析了图像中不同组分的分布特征、数字特征参数、所处环境及其与表面附近风化程度变化的关系，得到如下结论：

(1)云冈石窟砂岩中，长石与石英含量相差不大；乐山大佛砂岩中，长石含量远多于石英含量；石门山砂岩中，石英含量远多于长石含量。

(2)乐山大佛砂岩颗粒形状变化最小，云冈石窟砂岩次之，石门山砂岩形状最不规则；石门山砂岩纹理最复杂，乐山大佛砂岩最为简单均匀；石门山砂岩图像区域的变化规律最强，乐山大佛砂岩居中，云冈石窟砂岩最为杂乱。

(3)乐山大佛砂岩表面风化严重，但随深度增加、风化变化不明显；云冈石窟砂岩风化起初不明显，但深度达到25 mm后，风化程度骤增；石门山砂岩随深度增加、风化程度变化比较复杂。

确定两种或多种组分(包括出现基质和胶结物)时的分界灰度值，是本项研究的基础问题，目前还没有很好的解决方法；使用深度学习解决这一问题，联合使用灰度分布特征、不同组分的形态特征、深度学习的泛化能力，可能是一种比较好的解决方法。