基于真实细观尺度的沥青混合料三维重构算法

万成，张肖宁，贺玲凤，王端宜，张吉庆

(华南理工大学 土木与交通学院，广东 广州，510640)

沥青混合料具有损伤、非均匀、非连续的特性，其性能既受沥青、集料和空隙的体积含量影响，也受这些因素的空间分布影响[1]。在大多数的情况下，沥青混凝土被假设为均匀的，忽略其微细观结构的影响。造成这种现象的一个主要原因是在非均匀性或微细观结构的定量观测和评价方面所遇到的困难。例如，在很多研究中，采用统计数学描述材料非均匀性的方法(如Weibull分布等)[2]，但是，基于统计学的非均匀性和微细观结构研究其缺点在于其不能准确地反映实际材料的非均匀性和微细观结构。认识到基于沥青混合料细观尺度进行研究的重要性以来，计算机层析摄影技术或计算机层析识别技术(X-ray CT技术)在土木工程领域中的应用和发展，使得研究沥青混合料内部微细观结构成为可能。Masad等[3]在1999年采用X-ray CT获取圆柱形沥青混合料试件的内部结构组成。从此，沥青混合料细观领域的研究成为一个新的发展方向。此后，Wang等[4-6]利用X-ray CT考虑沥青混合料的实际分布形态，进行混合料内部二维结构的定量或定性分析，特别是 You等[7-8]利用高精度扫描仪采集沥青混合料试件的内部结构信息后，建立二维数值模型，并采用有限元方法建立二维黏弹力学数值模型，在细观力学数值模拟方面取得丰富的成果。陈俊[9]依据概率理论，推导沥青混合料二维截面上粗集料数量级配的计算式，运用计算机随机投放技术，并综合考虑集料不规则形状和级配特征，实现混合料二维数字试件的生成。但目前对混合料内部微细观研究多为基于二维的识别与分析。由于尚不能完全考虑集料、沥青及空隙的三维的体积特征，因而其研究结果在多单元体积元统计意义上具有代表性，而在考虑某个单元体时仍然具有一定的随机性与变异性。沥青混合料内部结构空间分布状况的三维重构将是今后研究的重点[10-11]。沥青混合料三维重构可用于分析混合料的三维体积组成，更可进行细观力学模拟，模拟各种荷载形式作用下混合料内部微细观结构的变化与发展。在此过程中，基于空间实际分布形态的真三维重构试样显得尤为重要。本文利用X-ray CT扫描实际沥青混合料试件，直接获取试样内部结构，根据沥青混合料三维重构原理，把握从三维可视化模型向三维数值模型转化的本质特点，基于X-ray CT技术对沥青混合料进行三维数值重构算法进行重点研究。

1 三维重构理论

图像三维重建就是要从二维的图像序列中提取三维对象的信息，使得用户能够直观地看到三维对象的组织结构，并且加强图像中原有的各种细节。三维重建的任务是从采样数据中恢复物体的三维结构，即建立物体的立体模型。重建的方法分为基于外表面的三维重建和直接由体素三维重构2大类。本文为基于体素的三维重构。

体数据的最基本的单元是体素，概念与像素类似。体数据集是定义在三维空间网格上的标量或向量数据，这些网络通常是正交网格。数据定义在网格节点上，相邻的8个网格节点构成1个立方体。如图1所示。假设在空间有一个长为(L-1)×VX，宽为(M-1)×VY，高为(N-1)×VZ的长方体，此长方体的面平行于坐标面。对于节点空间坐标(i, j, k)，i=0, 1, 2, …, L-1；j=0, 1, 2, …, M-1；k=0, 1, 2, …, N-1，式中 VX，VY，VZ为长度、宽度和高度方向的单位长度；i，j，k为编号。体数据与三维数组相对应，大小为L×M×N，定义在此三维空间网格的网格节点上，在节点(i, j, k)处，有体数据值f(i, j, k)[12]。

图1 体素定义示意图Fig.1 Sketch map of pixel definition

定义在此三维空间网格的各立方体上的单一体数据，称之为体素，即体素(i, j, k)有体数据值f(i, j, k)。针对断层图像序列而言，即每幅图像的1个像素被扩展为三维空间的1个匀质的体素，该体素的体数据值f(i, j, k)即为其所对应像素的灰度、颜色索引值或颜色值。可以根据1个算子，从体数据值f(i, j, k)得出体素(i, j, k)所对应的颜色分量C(i, j, k)及非透明度O(i, j,k)。根据上述原理，可知，每个体素由8个相邻节点组成，实现重建的关键技术是求出每个节点的三维坐标(i, j, k)，以及每个体素与自身相邻8个节点形成的对应关系。

2 三维重构算法

图2 CT系列连续图像Fig.2 Series of continuous CT images

利用CT扫描混合料试件后可以得到N幅一系列连续CT图像，如图2所示。对图像进行叠加，利用图像处理软件如VGStudio，IDL和Image Pro Plus等软件实现试样的三维可视化重构，如图3所示。必须指出的是三维可视化模型不能用来进行力学数值模拟，开发基于 CT扫描的三维图像技术来获取沥青混合料试件的内部结构，除了用于分析试样的内部体积组成，更重要的是为进行力学数值模拟，即虚拟力学试验，这就需要将三维可视化试样转换成三维数值试样。由于没有专用于沥青混合料的图像处理软件，并且处理过程相对复杂，实际上现有工业CT所附带的图像处理软件主要用于工业生产领域，如工业产品缺陷检查等。本文重点研究考虑空间真实细观结构的沥青混合料三维数值模型的重构算法，根据提出的算法开发沥青混合料试样三维重构程序，建立反映空间实际分布形态的沥青混合料三维数值模型，为虚拟力学试验奠定基础。

CT图像数量可以通过参数设定控制，理论上，只要获得的CT图像数量越多，相邻间距越小，越接近真实混合料试样。例如，对于一个标准马歇尔试件，其标准高度为63.5 mm，若获得1 001幅CT图像，则相邻2幅CT图像之间的间距为0.063 5 mm，这样小的图像间距对于建立的三维结构误差可以忽略不计。CT图像可以转换成256级的灰度图像，即8位黑白图像。对于每幅CT图像的灰度像素矩阵，假设其矩阵大小为mn，即像素个数为mn，因此1个像素可以代表1个长方体单元，每个长方体单元具有8个节点，多幅连续CT图像转换成的灰度像素矩阵，叠加形成的灰度单元体网格，如图4所示。取其中任意一幅CT图像的灰度像素矩阵，如图 5所示，其 1个像素由1个字节描述。0表示黑色，255为白色。从图5的像素矩阵示意图可见：每1个像素点是1个微小的正方形区域，可以用来表征材料的二维平面网格。根据上述原理和重构思想，可以得出结论：实现重建的关键技术是求出每个节点的三维空间坐标(i, j, k)，以及每个单元体与自身相邻8个节点的匹配逻辑关系。

图3 连续CT图像实现三维可视化试样重构Fig.3 3D reconstruction of visualization model with continuous CT images

图4 多幅CT图像的灰度单元体网格Fig.4 Gray unit net of multiple CT images

图5 单幅CT图像像素矩阵示意图Fig.5 Sketch map of pixel matrix of one image

2.1 像素单元体的节点编号

取图4中第1幅CT灰度图像的像素矩阵，其单元和节点编号示意图如图6所示。单元体编号I=1, 2,3, …, mn，节点编号J=1, 2, 3, …, (m+1)(n+1)。单元和节点编号规则为：从矩阵左下角第1个单元、节点起，从左至右、从下往上编号，每一行的末尾编号与下一行的起始单元、节点编号相连，从图 6可见：1幅CT灰度图像的像素矩阵的单元数目为mn，节点数目为(m+1)(n+1)。

取图4中第2幅CT灰度图像的像素矩阵，其单元和节点编号示意图如图7所示，单元、节点编号续接第1幅CT图像像素矩阵最末1个单元、节点编号，所以该图的第 1个单元、节点编号分别为 mn+1，(m+1)(n+1)+1，其余编号规则与第1幅相同。

比较第1幅与第2幅示意图的单元、节点编号可知：相邻2幅图片对应位置的单元编号相差mn，节点编号相差(m+1)(n+1)，所以对于图4中第k幅图片，其单元和节点编号如图8所示，为方便示意图中标记，记 P=(k-1)mn，Q=(k-1)(m+1)(n+1)。

2.2 单元体与节点对应

提取相邻两层图片中相邻的 8 个节点，构成 1个单元体。根据2.1节中单元和节点编号规则，可方便地用数学描述出每个单元体与自身8个节点之间的对应逻辑关系。

图6 第1幅CT图像像素矩阵单元和节点编号示意图Fig.6 Sketch map of elements and nodes label for the first CT image pixel matrix

图7 第2幅CT图像像素矩阵单元和节点编号示意图Fig.7 Sketch map of elements and nodes label for kth CT image pixel matrix

图8 第k幅CT图像像素矩阵单元和节点编号示意图Fig.8 Sketch map of elements and nodes label for kth CT image pixel matrix

2.3 节点的空间三维坐标

单元体与节点的匹配关系确定以后，如果能确定节点的空间三维坐标，相应的单元体位置也就得到确定。由2.1中节点编号可知，对于任意第k幅图片，其节点编号在(k-1)(m+1)(n+1)+1到k(m+1)(n+1)之间，设

则对于任意第k幅图片的节点三维坐标(i, j, k)可以表示为

至此，单元和节点的编号、对应关系已经确立，节点的三维空间坐标也已经确定。混合料的三维数值模型初步建立起来，如图9所示。必须要指出的是，此时建立的模型包含所有背景像素(CT图像中的黑色区域)和混合料像素所代表的单元和节点，所以呈现的是立方体形状，而实际混合料试件是圆柱体形状。

图9 三维数值初始模型Fig.9 3D initial digital model

2.4 混合料区域内单元的识别和提取

初始三维数值模型建立以后，需要对属于混合料区域内的单元进行识别和提取。背景区域内像素灰度值都为0，混合料区域内像素灰度在0～255之间，如果单纯以阈值分割的方法进行单元提取，如(0, T1)属于空隙，[T1, T2)属于砂胶，[T2, 255]属于集料，则会把混合料区域内属于空隙的单元(灰度为 0)遗漏或者把背景区域的像素包括进去。为此，采用指针扫描和阈值分割相结合的办法，对属于混合料区域内的单元进行提取。

所谓指针扫描即采用一定半径的指针，以试样中心为圆心点，不断旋转扫描，根据指定的判别标准，对区域内单元像素进行识别，区分出属于混合料区域的单元和背景区域的单元，并对区分出来的背景区域的单元像素赋予大于255的灰度，例如1 000，这样2个区域内的灰度分布完全区分开来，如图10所示。

2.5 阈值分割

混合料区域内的单元提取出来以后，需要对单元进行识别，区分出集料、砂胶和空隙。能否准确区分3者，阈值T1和T2的选取显得尤为重要。在有限元法构造实体模型中，最重要的是确定不同材料的几何边界，采用的图象处理技术包括图像的三值化和不同组分的边界提取技术。

CT图像的三值化首先要确定的是区分空隙与沥青砂胶的阀值T1及区分沥青砂胶和集料的阀值T2。阀值 T1和 T2的确定是数字转换过程中的重点问题，较为常用的方法有双峰法[13]导数法[14]。双峰法认为：灰度分布图中，2个峰值之间的应为同一介质；而导数法认为：灰度一阶偏导数分布图中，2个峰值之间的应为同一介质。李晓军等[15]分别利用双峰法和导数法进行了试算，结果发现对于CT图像，这2种方法差别不大，本文采用双峰法进行了计算。

还可以采用物质密度对比的方法进行阀值分割。分别对集料颗粒与配制的沥青砂胶进行CT扫描，获得的图像中确定2种物质的灰阶范围。

图10 指针扫描Fig.10 Pointer scanning

初步确定 T1和 T2后，采用反复试验方法微调阀值T1和T2，图11为三值化图CT图像。图像中红色部分分别为集料、砂胶和空隙。

2.6 三维数值试样的建立

利用开发的 MATLAB程序实现上述过程，最后利用MATLAB中的DLWRITE命令将数据文件写入构造有限元模型的INP文件，该INP文件包含了连续断层图像的细观结构信息，最终建立三维细观模型，如图12所示，区别于传统的均质体有限元数值模型，该模型真实地反映实际试样中集料、砂胶和空隙的三维细观分布，即考虑3者的空间真实分布形态，只需对3者分别赋予不同的材料属性，为以后的沥青混合料虚拟力学试验奠定基础。

图11 3组分的分割Fig.11 Image segmentation of three constituents

图12 三维数值模型Fig.12 3D digital specimen

3 结论

(1) 对基于X-ray CT技术的沥青混合料三维重构方法进行研究。分析沥青混合料三维重构原理，提出了对沥青混合料进行三维数值重构的基本算法，开发沥青混合料试样三维重构程序，该数值试样区别于传统的均质体有限元模型，而是基于集料、沥青砂胶和空隙的真实三维空间分布，因而更加贴近实际情况。

(2) 基于沥青混合料真实细观尺度建立的三维数值模型，是整个沥青混合料虚拟力学试验的前提与基础，有助于以后从三维细观尺度研究沥青混合料的力学行为。

(3) 限于计算机运算速度限制，对 CT 图像只能采取降低分辨率的方法，从而影响到了集料、砂胶和空隙的分割精度，随着计算机技术的快速发展以及算法的优化，可以预见将获得更为高精度的三维数值模型，但本研究从方法上解决基于真实三维细观分布的沥青混合料重构问题。

[1] TAN Yi-qiu, XU Hui-ning, LI Xiao-min. Is normal distribution the most appropriate statistical distribution for volumetric properties in asphalt mixtures?[J]. Journal of Testing and Evaluation, 2009, 27(5): 458-478.

[2] 邱延峻, 闫常赫, 艾长发. 非均质沥青混合料劈裂试验全过程数值模拟[J]. 交通运输工程学报, 2009, 9(2): 12-16.QIU Yan-jun, YAN Chang-he, AI Chang-fa. Numerical simulation of split test process for asphalt mixture under heterogeneous state[J]. Journal of Traffic and Transportation Engineering, 2009, 9(2): 12-16.

[3] Masad E, Muhunthan B, Shashidhar N, et al. Quantifying laboratory compaction effects on the internal structure of asphalt concrete[C]//Proceeding of Transportation Research Record(TRB) 78th Annual Meeting. Washington D C: Transportation Research Board, 1999: 179-185.

[4] Wang L B, Wang Y, Mohammad L, et al. Voids distribution and performance of asphalt concrete[J]. International Journal of Pavements, 2002, 1(3): 22-33.

[5] Tashman L, Masad E, Angelo J, et al. X-ray tomography to characterize air void distribution in superpave gyratory compacted specimens[J]. International Journal of Pavement Engineering, 2002, 3(1): 19-28.

[6] Kim H, Wagoner M P, Buttlar W G. Numerical fracture analysis on the specimen size dependency of asphalt concrete using a cohesive softening model[J]. Construction and Building Materials, 2009, 23(5): 2112-2120.

[7] You Z P, Dai Q L. Dynamic complex modulus predictions of hot-mix asphalt using a micromechanical-based finite element model[J]. Canadian Journal of Civil Engineering, 2007, 34(12):1519-1528.

[8] Dai Q L. Prediction of dynamic modulus and phase angle of stone-based composites using a micromechanical finite-element approach[J]. Journal of Materials in Civil Engineering, 2010,22(6): 618-627.

[9] 陈俊. 基于离散元方法的沥青混合料虚拟疲劳试验研究[D].南京: 东南大学交通学院, 2010: 101-120.CHENG Jun. Virtual fatigue tests of asphalt mixture based on discrete element method[D]. Nanjing: Southeast University.School of Transportation, 2010: 101-120.

[10] 汪海年, 郝培文. 沥青混合料微细观结构的研究进展[J]. 长安大学学报: 自然科学版, 2008, 28(3): 11-14.WANG Hai-nian, HAO Pei-wen. Advances in microstructure study on asphalt mixture[J]. Journal of Changan University:Natural Science Edition, 2008, 28(3): 11-14.

[11] You Z P, Dai Q L. Review of advances in micromechanical modeling of aggregate–aggregate interactions in asphalt mixtures[J]. Canadian Journal of Civil Engineering, 2007, 34(2):239-252.

[12] 王旭, 赵歆波, 桂业英, 等. 医学图像序列的直接体绘制算法研究[J]. 生物医学工程学杂志, 1999, 16(3): 324-328.WANG Xu, ZhAO Xin-bo, GUI Ye-ying, et al. A new algorithm for direct volume rendering of medic image series[J]. Journal of Biomedical Engineering, 1999, 16(3): 324-328.

[13] Kranz R L. Crack-crack and crack-pore interactions in stressed granite[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1979, 16(1): 37-47.

[14] Wu F J, Thomsen L. Microfracturing and deformation of westerly granite under creep condition[J]. International Journal of Rock Mechanics and Mining Sciences, 1975, 12(5): 167-173.

[15] 李晓军, 张肖宁, 武建民. 沥青混合料单轴重复加卸载破损CT识别[J]. 哈尔滨工业大学学报, 2005, 37(9): 1228-1230.LI Xiao-jun, ZHANG Xiao-ning, WU Jian-min. Internal structure discrimination of asphalt mix with X-ray computerized tomography under repetitive uniaxial loadings[J]. Journal of Harbin Institute of Technology, 2005, 37(9): 1228-1230.