基于图像重建方法的岩体损伤形态表征研究*

王昱博，陈 清

（合肥工业大学 土木与水利工程学院，安徽 合肥 230009）

1 引言

岩体所受外部荷载情况复杂，荷载造成的内部损伤的起裂部位、扩展路径隐藏于岩体内部，难以直接观测到损伤区域发展过程。因此，岩体损伤演化路径和裂隙空间形态的表征对于研究岩体损伤演化的力学行为具有重要意义［1］。现阶段针对岩体（石）的三维重建大体上分为两种：一类是利用数字图像相关（DIC）、扫描电子显微镜（SEM）、计算机断层扫描（CT）、声发射（AE）等无损技术获取岩石真实细观结构，经处理在软件中进行三维重建还原，这种方法精细度高，但对设施、场地要求比较高，同时往往需要处理大量数据；另一类的方法从数值模拟的角度出发，通过程序算法模拟岩体（石）细观损伤，利用二维位图重构三维形态，这种方法成本较低、具有较好的可重复性。

国内外学者对岩体（石）细观结构的三维重建开展了相关研究，取得了丰硕的成果。黄家国等采用扫描电镜、核磁共振和计算机断层扫描3种方法重建了孔隙页岩的高密度三维模型［2］；熊自明等基于ASIFT及RANSAC算法，将隧道掌子面的数字摄影进行处理，实现非接触式隧道掌子面裂隙纹理的重建［3］；Zhou X P等利用CT对含孔隙砂岩进行研究，采用分层退火算法处理层位图像重建砂岩的不规则特征［4］；张亚等利用Monte-carlo随机模拟方法将水电站坝址区裂隙几何参数导入模拟软件中形成三维裂隙网格模型，对裂隙岩体渗流的研究带来一定帮助［5］；何延龙等通过退火算法与分形算法结合的Hoshen-Kopelman算法将黏土矿物基团的主要结构特征分类，构建了含黏土矿物分布的三维重建多孔介质模型，较好地反映了主要黏土矿物的分布及产状［6］；张挺等基于多孔介质二维图像和多点统计方法，通过训练图像形成三维空间中的体素，重建出多孔介质的重建模型［7］；赵延林等利用FLAC3D实现三维裂隙的随机重建，并在重建模型中进行细观渗流计算，分析孔隙结构对渗流的影响［8］。

学者们研究成果具有重要的应用价值，但对于工程岩体而言，所研究对象基本为岩石试件，与工程尺度的岩体存在较大数量级差距，且模型重建过程相对复杂，方法较难遵循。本文在前人研究的基础上，基于简便易行的图像三维重建方法，结合数字图像处理技术、数值模拟，于FLAC3D中建立含软弱层的隧道数值模拟模型，对隧道模型掘进过程中造成的损伤进行三维重建。

2 图像重建方法

利用图像进行三维重建的方法主要有三种，直接重建法、间接重建法和混合重建法［9］。其中直接重建法和间接重建法较为常用。

2.1 直接重建法

直接重建法大致分为四种算法，光线投射法、足迹表法、错切-变形法、体元投射法［10］。这类算法是一种根据图片体素进行投影从而重构三维模型的技术，不生成中间几何曲面。现今使用的直接重建法基本由光线投射法和足迹表法发展而来。光线投射法以一条光线穿过位图序列构成的体数据，在这条线上对体素进行等距采样，设置光照模型获得、重构模型，示意图见图1。而足迹表法则以构成三维体数据的位图存储顺序采集体素信息，读取某种体素在每一幅位图中的影响范围，往一个平面投影体素，得到投影重构模型。直接重构法计算生成的三维图像质量好、稳定，对于需要观察内部结构信息的三维模型具有优势，缺点是计算量较大、计算时间长。

图1 光线投影法原理

2.2 间接重建法

间接重建法是CT技术进行快速模型重构常用的一种方法。这种方法从图像构成的数据集中根据阈值提取所需表面，并将此表面多边形化，用不同的方法拼接多边形，形成大量的中间曲面，叠加曲面以构成三维表面模型［9］。间接重建法主要有移动立方体法、轮廓拼接法、移动四面体法、剖分立方体等算法［10］。典型的间接重建法是由Lorensen［11］提出的移动立方体算法（MC），如图2所示，该法通过阈值或极值等体素在立方体中剖分三角网格来近似等值面，将三角形网格组合拼接后，计算梯度向量为每个三角形生成法线，以便渲染阴影形成光照效果，构建出物体的三维形态。这类算法简单、重构速度快，但包括MC算法在内的间接重构法存在一些通病：有图像前处理时间长、模型精细度较低和内部信息较难观察等问题。

图2 移动立方体算法的基本组合［11］

对于岩体的三维损伤重构模型来说，具有良好的可视性、精细度和易于观察内部损伤状态等特性是十分必要的。matlab具有强大的可视化能力，间接重构法和直接重构法都可以在其中实现，但直接重构法计算生成的三维模型质量更好、更加精细，模型内部信息均可呈现，用来展现岩体损伤的三维形态是一种较好的方法。因此，本文选择在matlab中运用直接重构算法进行基于图像的岩体损伤重建。

3 重建图像的获取与损伤辨识

3.1 图像切制

在计算机断层扫描技术中，位图的切制是按照一定顺序或一定方向（或者理解为与位置有关）通过扫描或剖切研究对象得到的一系列二维平面图。本方法中，首先在数值模拟软件中编制损伤特征参量函数、图片自动切制函数，实现位图的自动切制。切制的方法类似于计算机断层扫描层位图的获取方式，在特征参量指定量值范围的云图上沿模型高程方向，以一定的间距dd为距离切制多个平面，形成位图图序。在切制位图时，减小间距dd以获得更多位图图像，利于详细描述损伤的重构特征。图3为损伤部位的提取与重建示意。

图3 损伤部位的提取与重建

3.2 损伤部位辨识

位图中损伤特征参量的辨识，其本质就是通过技术手段对损伤特征参量位图中关键像素点的位置和颜色信息进行辨识。为此，我们需要寻找合适的量化评价岩体裂隙宏细观力学效应的特征参量，获取损伤参量的位图。

YAI（yield approache index）［12-13］是一种基于Mohr-Coulomb强度破坏准则的评价岩体危险性程度的定量指标。YAI估计了应力诱发岩体屈服，出现破裂的风险程度，YAI值在0～1内时岩体处于发生损伤的塑性状态，裂纹逐渐蔓延，当大于1时即宣告岩体的破坏。本文以此指标量化岩体损伤程度并于FLAC3D中建立数值模拟模型，在数值模拟模型的YAI云图上显示岩体损伤的范围及程度，计算完成后切制多张岩体损伤状态的位图，供辨识和提取。

由于位图是由一系列带有颜色特征像素点进行不同的排列和染色所组成的，位图中不同的颜色表示着相应像素点带有不同的颜色阈值，故进行图像处理时，位图只需要显示发生损伤部分。包含损伤信息的YAI位图和大多数图片一样是RGB颜色格式，但实际上RGB不能很好地反应图像的形态特征，计算时也存在数据量较大的问题，所以我们所获得的位图需要进行灰度处理［14-16］。灰度处理的方法分为线性处理和非线性处理两大类，常用的线性处理算法是经典心理学算法，这种算法通过确定三种颜色的权重将位图灰度化，图像处理过程中采用式（1）的RGB分量加权平均算法将彩色图处理成灰度图，此种RGB分权方式得出的灰度图效果最为理想，可以将岩体损伤区域与其他区域清晰地区分出来。流程见图4。

其中，f为处理后的灰度值，R为红色灰度值，G为绿色灰度值，B为蓝色灰度值。

图4 图像处理流程

对于一幅处理完成的灰度位图，图中像素点的灰度值为0～255之间。移动鼠标至matlab中灰度图的某一区域时，界面上会显示此区域的灰度值，此时，以此区域的灰度值将当前区域与其他区域分割开，这种成为图像分割界限的灰度值被称为阈值。因此，分割出位图中损伤区域与非损伤区域的关键是找到损伤区域的阈值。如图4中损伤区域的灰度值为76，以76为阈值灰度图中的损伤区域分割出来，使用式（2）进行逻辑计算，将位图分割成由黑、白色的像素点所组成的二值图。后期在重构损伤模型的步骤中辨识和提取二值图中的损伤部分即可。处理好的二值图的颜色阈值只有单个变量，在后期识别过程中控制单个变量即可实现辨识过程，简化了辨识和提取的难度。

其中g（i,j）为逻辑运算符，i、j为像素点坐标，f（i,j）为像素点灰度值，T为阈值。

由此，可以利用Matlab精确地识别出二值图中岩体的损伤区域。此外，在数值模拟软件的自动出图阶段应对每一幅图进行有规律地命名，在Matlab中根据同样命名规律运用循环语句将选用的位图序读入到矩阵中，以实现图像的逐幅提取，达到自动处理图像的功能，如图5所示。

图5 二值图位图图序

4 基于位图图序的三维重构

4.1 数值模型描述

内含软弱带的非均质隧道进行损伤重建。马蹄形隧道直径10m，隧道左上方含一条软弱带。隧道轴线沿地应力的大主应力方向，模型的尺寸为长30m，宽15m，高30m，模型单元数为338333，节点数为58282。数值模拟模型如图6所示。

表1 岩石力学参数

图6 开挖隧道计算模型

该隧道轴线沿地应力中间主应力方向，大主应力为0.6MPa；中间主应力为侧向压力0.45MPa；最小主应力为垂直方向0.26MPa。隧道的开挖方式采用全断面一次性开挖，不考虑支护结构，隧道进尺为3m一个工况，共分5个工况开挖。

4.2 损伤模型重建

依据开挖工况及已介绍的方法截取重建损伤模型所需位图，截取间距为0.05m，自下而上截取，获取模型的YAI云图特征参量位图600幅，为便于辨识损伤区域，将位图进行二值化处理，损伤区域为位图中灰度值为0（黑色）的部分，部分位图见图7及图8。

图7 部分灰度位图

图8 部分二值化位图

利用二值化的位图进行岩体损伤区域的三维重建，直观展现损伤区域的空间形态，具体步骤如下：

（1）封装成三维体数据。

将由FLAC3D导出的n幅损伤位图以命名序列读入matlab中，利用cat函数对选取的N幅图像数据进行三维体数据集的构造，得到X×Y×N的三维矩阵。

（2）优化三维体数据字段量。

所形成的三维体数据由多幅图像构成，为了防止计算时内存溢出，采用reducevolume和smooth函数减少数据字段和平滑处理，减少的数据量比例越大体数据集中数据越少，因此在保证运算顺利的基础上应尽量抽取更多的数据以保证三维重构的完整性。

（3）转换位图坐标。

三维体数据矩阵D的数据量减少时，模型像素点坐标可能会发生改变，需要检查模型高度和宽度方向所占的像素点坐标与其在模型对应点的实际坐标的一致性，参考模型在实际坐标中的位置对选取云图的像素点矩阵x、y、z进行坐标转换。

（4）计算岩体损伤百分率。

通过识别每张位图损伤区域的像素面积乘以高度值，计算损伤总体积，显示于重建模型旁。

（5）提取等值面进行三维构造。

在完成模型像素点坐标转换后，调用函数isosurface从体数据提取等值数据绘制等值面。二值图中损伤区域为黑色，黑色的灰度值为0，提取出三维体数据中灰度值为0的网格点，连接仰角相同的点以形成等值面，在显示平面累计投影，将等值面的顶点和面的信息传递给patch函数，定义和构造所选区域图像的颜色和光线。

（6）优化显示效果。

最后，基于三维体数据调用isonormal、lighting和view等函数及命令分别定义所选区域的光照效果、光线阴影以及观察视角等属性，通过绘制模型边界轮廓图显示出破损三维体在模型中的具体位置，重构模型如图9所示。

图9 YAI云图及隧道损伤特征三维空间形态

图9中的分图分别为不同开挖步对应的岩体损伤特征三维空间形态，伴随着开挖的进行，隧道各开挖步造成的岩体损伤在三维空间中逐步累积，由于受到隧道左上方的软弱层的影响，岩体损伤区域靠近软弱层一侧的发展和分布较多，损伤重建模型很好地反应了这一点。结合图9和图10可以看出，初步开挖3m时岩体的损伤体积(DR)损伤率(DV)分别为11.4183m3和0.090621%，至隧道开挖完成，DR和DV分别为223.4181m3和1.7732%，岩体损伤程度得到了量化，损伤的累积变化近似看作线性关系。

图10 岩体损伤体积和损伤率

因此，根据上述数值模拟和三维重建的结果，可以认为该方法较好地量化了工程岩体的损伤程度、预测了工程作业过程可能造成的岩体损伤，形象地重现了模型中岩体损伤的形态和分布情况。

5 结论

本文基于图像三维重建技术，结合数字图像处理技术、数值模拟，进行了岩体损伤的表征，得出以下结论：

（1）通过设置不同的特征参数，图像重建方法能够重建诸如位移、应力、YAI、单元安全度等三维形态图，有助于综合分析岩体损伤机制，从而实现岩体内部损伤三维形态的有效表征与直观展示。

（2）隧道损伤重建案例表明：通过辨识和重构的岩体内部损伤空间结构模型，能够清晰地展现工程岩体损伤区的空间分布形态，有效分析、推测后续工况中工程岩体损伤的发生位置、扩展路径、演化特征等岩体损伤规律。

（3）图像重建方法不受工程尺度效应限制，可作为传统岩体无损检测方式的补充方案，辅助分析工程岩体稳定性相关问题，为工程灾害控制的优化设计提供一定的参考，具有较好的适用性和工程价值。