CN104599321A - 基于X-ray CT图像的真实集料颗粒离散元模型构建方法 - Google Patents

Abstract

本发明公开了一种基于X-ray CT图像的真实集料颗粒离散元模型构建方法，为应用于离散元数值分析的沥青混合料试件中的集料颗粒提供一种基于X-ray CT图像的真实集料颗粒离散元模型构建方法，实现能高度拟合真实集料几何形状兼具最小离散单元数量的集料颗粒离散元模型的自动构建。本发明基于沥青混合料试件细观结构模型的离散元分析，完成从沥青混合料试件X-ray CT断面图像到精确拟合真实集料颗粒几何形状兼具最少单元数量的集料颗粒离散元模型的自动转换，以辅助建立沥青混合料的离散元模型。

Description

基于X-ray CT图像的真实集料颗粒离散元模型构建方法

技术领域

本发明涉及数值分析领域，具体是一种基于X-ray CT图像的真实集料颗粒离散元模型构建方法。

背景技术

三维细观力学分析技术是一种通过研究沥青混合料内部细观结构来评价沥青路面宏观路用性能的重要技术。其中，建立能够反映路面试件真实细观构成的三维数值模型是获得可信分析结果的重要前提。

在数值分析方法中，离散元方法在准确模拟沥青混合料的变形、开裂等细观力学行为方面有着明显的优势。由于混合料内部的集料颗粒在混合料组分中占有极大的比例，其几何形状对混合料的细观力学性能有重要影响，用以表示集料颗粒的离散单元集合能否精确拟合真实集料颗粒的几何形状对于仿真可信度十分关键。为了使真实集料颗粒的离散元模型能够精确拟合该颗粒的几何形状，需要对真实集料颗粒进行精确三维实体建模，在此基础上完成从该颗粒的实体模型到离散元模型的合理转换。

目前比较常见的集料颗粒离散元模型构建方法主要分为两类，一类是基于混合料X-ray CT断面图像的集料颗粒离散元建模，另一类是基于既定数学模型的集料颗粒离散单元集合生成。

前者通过使用X-ray CT对混合料试件进行断面扫描，以无损方式获取其断面图像，通过设定灰度阈值以区分图像中的三相材料，将每个集料像素构建为一个离散单元。但三相材料像素的灰度范围通常存在重合部分，仅仅通过设定灰度阈值判断各相材料像素，其准确性很难保证；而离散单元与像素的一一对应关系也会造成单元数量的大幅增加，从而造成巨大的计算量和存储量需求。后者通过开发一种球形离散单元产生算法生成虚拟集料颗粒的离散单元集合，而这种虚拟生成的离散单元集合与真实集料颗粒在几何形状上存在较大差异，从而影响混合料细观力学性能分析的可靠性。

美国密歇根理工大学(Michigan Technological University)2011届博士学位 论文《Discrete Element Methods for Asphalt Concrete:Development and Application of User-defined Microstructural Models and a Viscoelastic Micromechanical Model》采用两种方式：一是随机生成具有不同球度、方位、棱角性等属性参数的虚拟集料颗粒的离散元模型；二是从混合料X-ray CT扫描图像中识别出集料像素并构建为集料球形离散单元。在此基础上构建混合料试件的离散元模型并使用PFC3D软件进行了动态模量和蠕变刚度测试的仿真。其中，通过第一种方式获得的混合料离散元模型由于虚拟集料颗粒与真实集料颗粒较大的几何形状差异，并不适合预测混合料的细观力学性能；而通过第二种方式获得的混合料离散元模型中的集料颗粒离散元模型对相应颗粒几何形状的拟合精度较低，而单元数量却较多，这将大幅增加离散元仿真的计算时间。国家发明专利《基于沥青路面试件X-ray CT图像的集料细观实体模型重构方法》(专利号：ZL201210172375.7)采用计算机图像处理和图形学三维建模技术，依据混合料X-ray CT图像重构了混合料试件中所有集料颗粒的实体模型，避免了像素与有限元单元一一对应的单元生成方法，极大地提升了沥青混合料的有限元数值仿真效率。

本发明在国家发明专利《基于沥青路面试件X-ray CT图像的集料细观实体模型重构方法》(专利号：ZL201210172375.7)的基础上，开发集料颗粒轮廓识别技术，在准确识别X-ray CT图像中集料轮廓的基础上建立集料颗粒的高精度三维实体模型，并采用模型转换算法，获得能高度拟合真实集料颗粒几何形状并拥有最小单元数量的集料颗粒离散元模型。

发明内容

本发明的目的是提供一种基于X-ray CT图像的真实集料颗粒离散元模型构建方法，以解决现有技术中集料颗粒离散元模型对集料颗粒真实形状的拟合精度低及模型中单元数量巨大的问题。

为了达到上述目的，本发明所采用的技术方案为：

基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：包括以下步骤：

(1)、以BMP格式、分辨率为512×512像素、扫描间隔一定的沥青混合料试件X-ray CT灰度图像文件组为数据源，选择一个集料颗粒并截取其在各扫描图像中的投影图块，并保存为BMP格式的灰度图像文件组Aggregate[n]，用 以重构该集料颗粒的实体模型，其中：

所述沥青混合料试件是指以沥青材料和集料按适当比例配制，胶结成整体的圆柱形样本，记半径为radius，高度为height；

所述X-ray CT灰度图像文件组是指使用X-ray CT设备以一定纵向间隔扫描沥青混合料试件而获得的断面扫描图像数组，每幅断面扫描图像保存为BMP格式的灰度图像文件；

所述灰度图像文件是指将彩色图像中像素的颜色分量值按照一定比例进行加权，并将加权值附于每个颜色分量而获得的图像文件，其中，灰度图像中所有像素的灰度均介于0到255之间；

所述投影图块是指一幅扫描图像中包围目标集料颗粒图像的最小矩形内的图块；

所述实体模型是指以边界表示法表示的三维几何结构，记为模型M＝(FACE,EDGE,VERTEX,R)，其中FACE表示模型中面的集合，EDGE表示边的集合，VERTEX表示顶点的集合，R表示模型中各元素之间的邻接关系；

记录着面的类别：平面、样条面，以及面的曲面方程，其中face为FACE中的元素；

记录着边的曲线方程，其中edge为EDGE中的元素；

记录着顶点的坐标p(x,y,z)，其中vertex为VERTEX中的元素；

R＝{r₁,r₂}，其中：

r₁＝{(face₁,face₂,edge)|edge∈face₁∩face₂,edge∈EDGE,face_i∈FACE,i＝1,2}，r₁表示face₁和face₂相邻于边edge；

r₂＝{(edge₁,edge₂,vertex)|vertex∈edge₁∩edge₂,vertex∈VERTEX,edge_i∈EDGE,i＝1,2}，r₂表示边edge₁与edge₂相邻于点vertex；

(2)、遍历灰度图像文件组Aggregate[n]，使用轮廓检测及建模算法，检测出每幅图像Aggregate[i](i＝1,2,…,n)中集料颗粒的轮廓像素，存入轮廓像素数组P_i(i＝1,2,…,n)，再建立相应的轮廓面实体模型OF_i(i＝1,2,…,n)并存储；

(3)、输入集料颗粒所对应的全部轮廓面实体模型{OF_i|i＝1,2,…,n}，调用ACIS7.0的底层建模功能Skinning，获得BODY类型的集料颗粒实体模型M；

(4)、输入M，利用模型转换算法，获得M以球形离散单元形式表示的实体模型DM{Sphere_k|k＝1,2,…,p}，再获得以文本文件格式TXT存储的球形离散单元信息，导入PFC3D 4.0软件获得集料颗粒的离散元模型。

所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述步骤(2)中，轮廓检测及建模算法描述如下步骤：

1)、将图像Aggregate[i]的中心像素及其八邻域像素分别作为该图像所对应灰度差分八叉树GDO_i的根结点及根结点的子结点，并置此9个像素的访问标志bV＝true及根结点所有子结点的可扩展标志bE＝true后，计算GDO_i中所有连接线段的灰度差分值d_FC并将当前层号h置为2，最大深度H置为Aggregate[i]纵、横向分辨率中较小者一半的整数部分；

2)、使用灰度差分阈值自动分析算法获得灰度差分阈值T_G，依次判断每个位于GDO_i中第h层的结点与其父结点间连接线段的d_FC是否大于T_G，若是，则该结点的对应像素是轮廓像素，将该像素存入P_i并置该结点的bE＝false，否则仅置该结点的bE＝true；

3)、遍历GDO_i中第h层的结点若的bE＝true，则将对应像素的八邻域像素中bV＝false的像素作为的子结点，并置所有子结点的bE＝true、子结点对应像素的bV＝true，使h增1；

4)、持续运行步骤2)～步骤3)直至h＝H或者GDO_i中第h层上的结点所对应像素均为轮廓像素为止，完成构建GDO_i；

5)、依据P_i中各像素在图像Aggregate[i]中的邻接关系，对P_i中的像素进行分组排序，将在图像中依次邻接的像素按邻接顺序排列并归于一组，若P_i只包含 一组像素，执行步骤7)，否则执行步骤6)；

6)、对P_i中每个像素组的两个末端像素，分别找出与之距离最短的另一像素组的末端像素与之配对，以配对像素之一为起点，另一配对像素为终点，以起点和终点间的连线方向为检测方向，将T_G逐次减1并检测轮廓像素，直至起点和终点间的轮廓像素在两者连线方向上依次相邻；

7)、调用ACIS7.0的底层建模功能Spline_Fitting拟合P_i中的像素，获得FACE类型的轮廓面实体模型OF_i；

其中：

所述八邻域像素是指图像中一个像素的上、下、左、右、左上、右上、左下和右下八个方向上与之邻接的像素；

所述灰度差分八叉树GDO是指用以描述集料颗粒轮廓像素检测过程中所涉及像素间的邻接关系及邻接像素间的灰度差分的八叉树；

所述八叉树是指由若干结点按照结点间的父子关系以直线段连接而成、且每个结点最多有八个子结点的一种树形数据结构；

所述灰度差分值d_FC是指GDO内子结点与父结点对应像素间的灰度差；

所述当前层号h是指GDO_i在构建过程中的当前层数；

所述最大深度H是指所建GDO允许包含的最大层数；

所述轮廓面是指以P_i中像素的拟合结果为边界的平面。

所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述轮廓检测及建模算法的步骤2)中，灰度差分阈值T_G自动分析算法描述如下步骤：

1)、将集料像素灰度范围上限赋值为Aggregate[i]所对应灰度直方图中的最大灰度值；

2)、将变量t赋值为

3)、让t每次递减1，直至中像素占Aggregate[i]中全部像素的比例超过集料像素比例p_AG；

4)、让t每次递减1，直至灰度为t的像素中不再包含集料像素为止，将集料像素灰度范围下限赋值为t；

5)、让t每次递减1，直至灰度为t的像素中不再包含沥青胶浆像素为止，将沥青胶浆像素灰度范围上限赋值为t；

6)、将灰度差分阈值T_G赋值为与的差值；

其中：

所述集料像素灰度范围上限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的上限；

所述灰度直方图是指相应图像中各灰度级上像素数量的统计图；

所述集料像素比例p_AG是指集料像素占Aggregate[i]全部像素的预估比例；

所述集料像素灰度范围下限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的下限；

所述沥青胶浆像素灰度范围上限是指Aggregate[i]所对应灰度直方图中沥青胶浆像素灰度范围的上限；

所述灰度差分阈值T_G是指用以判断两任意像素是否分别为集料和沥青胶浆像素的灰度差分最小值。

所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述步骤(4)中，模型转换算法描述如下步骤：

1)、使用ACIS的底层建模功能Get_Bounding_Box获得M的包围盒B_M，以B_M的中心位置为中心，B_M最长边的长度为边长，建立BODY类型的立方体初始剖分空间IDS，并将其作为M所对应的转换八叉树TO的根结点root；

2)、使用ACIS的底层建模功能Mass_Proporties计算IDS与M的体积重合比例R_V，若R_V大于转换精度阈值T_A，那么将root的占有状态标志S赋值为占有状态fo，否则，若R_V＜1-T_A，则将S赋值为不占有状态cno；否则将S赋值为部分占有状态po；保存该信息到root，并标记root为未访问；

3)、将IDS剖分为8个等体积的立方体，作为TO中root的子结点Node₁,…,Node₈并标记为未访问，标记root为已访问；

4)、广度优先遍历TO，将每一个S＝fo或者S＝cno的Node标记为已访问，对每个未访问且S＝po的Node，若其对应立方体体积的1/8大于预设最小体积V_min，将该立方体剖分为8个新的立方体，将新增立方体作为该Node的子结点并标记为未访问后，为各立方体对应的S赋值，并将该Node标记为已访问；

5)、持续执行步骤4)直至所有Node都标记为已访问或者Node对应立方体不用再度剖分，完成建立TO；

6)、遍历TO，对所有S＝fo的Node，建立相应立方体的内切球实体模型，获得M的球形离散单元形式；

其中：

所述包围盒B_M是指将集料颗粒实体模型M完全包围起来的最小长方体；

所述初始剖分空间IDS是中心位置与B_M相同、并能包含B_M的最小立方体；

所述转换八叉树TO是指用于描述剖分产生的各立方体之间相互覆盖关系和状态的八叉树结构；

所述体积重合比例R_V是指立方体与M重合部分的体积与立方体体积的比值；

所述转换精度阈值T_A是指预设用于判定立方体被M占有的体积重合比例R_V的最小值； 

所述占有状态标志S是指用于标识立方体与M体积重合程度的标志，包括占有状态fo，不占有状态cno和部分占有状态po；

所述广度优先遍历是指按照从上到下、从左到右的顺序依次遍历八叉树中结点的方式；

所述预设最小体积V_min是指预设用于判断剖分是否终止的体积最小值；

所述球形离散单元形式是指M经模型转换算法转换而成的由尺寸不一的球 组成的实体模型；

所述球形离散单元信息是指DM{Sphere_k|k＝1,2,…,p}中每个球Sphere_k的球心和半径；

所述离散元模型是指可用于离散元数值仿真的计算模型。

本发明中，ACIS7.0是由美国Spatial公司生产的基于面向对象软件技术的三维几何造型引擎，PFC3D 4.0是由ITASCA咨询集团开发的一种主要用于模拟任意形状、大小的二维圆盘或三维球体集合体的运行及其相互作用的颗粒流分析程序。

本发明的有益效果如下：

本发明基于沥青混合料试件细观结构模型的离散元分析，完成从沥青混合料试件X-ray CT断面图像到精确拟合真实集料颗粒几何形状兼具最少单元数量的集料颗粒离散元模型的自动转换，以辅助建立沥青混合料的离散元模型。

附图说明

图1a为某沥青混合料试件X-ray CT扫描图像序列。

图1b为选取的单个集料颗粒在试件X-ray CT扫描图像序列中的投影图块。

图2a为集料颗粒在各投影图块中轮廓像素的检测结果图。

图2b为集料颗粒在各投影图块中对应的轮廓面模型。

图3为集料颗粒对应的实体模型。

图4a为集料颗粒实体模型转换而成的以球形离散单元形式表示的实体模型。

图4b为集料颗粒对应的离散元模型。

具体实施方式

基于X-ray CT图像的真实集料颗粒离散元模型构建方法，包括以下步骤：

(1)、以BMP格式、分辨率为512×512像素、扫描间隔一定的沥青混合料试件X-ray CT灰度图像文件组为数据源，选择一个集料颗粒并截取其在各扫描图像中的投影图块，并保存为BMP格式的灰度图像文件组Aggregate[n]，用以重构该集料颗粒的实体模型，其中：

沥青混合料试件是指以沥青材料和集料按适当比例配制，胶结成整体的圆柱形样本，记半径为radius，高度为height；

X-ray CT灰度图像文件组是指使用X-ray CT设备以一定纵向间隔扫描沥青混合料试件而获得的断面扫描图像数组，每幅断面扫描图像保存为BMP格式的灰度图像文件；

灰度图像文件是指将彩色图像中像素的颜色分量值按照一定比例进行加权，并将加权值附于每个颜色分量而获得的图像文件，其中，灰度图像中所有像素的灰度均介于0到255之间；

投影图块是指一幅扫描图像中包围目标集料颗粒图像的最小矩形内的图块；

实体模型是指以边界表示法表示的三维几何结构，记为模型M＝(FACE,EDGE,VERTEX,R)，其中FACE表示模型中面的集合，EDGE表示边的集合，VERTEX表示顶点的集合，R表示模型中各元素之间的邻接关系；

记录着面的类别：平面、样条面，以及面的曲面方程，其中face为FACE中的元素；

记录着边的曲线方程，其中edge为EDGE中的元素；

记录着顶点的坐标p(x,y,z)，其中vertex为VERTEX中的元素；

R＝{r₁,r₂}，其中：

r₁＝{(face₁,face₂,edge)|edge∈face₁∩face₂,edge∈EDGE,face_i∈FACE,i＝1,2}，r₁表示face₁和face₂相邻于边edge；

r₂＝{(edge₁,edge₂,vertex)|vertex∈edge₁∩edge₂,vertex∈VERTEX,edge_i∈EDGE,i＝1,2}，r₂表示边edge₁与edge₂相邻于点vertex；

(2)、遍历灰度图像文件组Aggregate[n]，使用轮廓检测及建模算法，检测出每幅图像Aggregate[i](i＝1,2,…,n)中集料颗粒的轮廓像素，存入轮廓像素数组P_i(i＝1,2,…,n)，再建立相应的轮廓面实体模型OF_i(i＝1,2,…,n)并存储；

(3)、输入集料颗粒所对应的全部轮廓面实体模型{OF_i|i＝1,2,…,n}，调用ACIS7.0的底层建模功能Skinning，获得BODY类型的集料颗粒实体模型M；

(4)、输入M，利用模型转换算法，获得M以球形离散单元形式表示的实体模型DM{Sphere_k|k＝1,2,…,p}，再获得以文本文件格式TXT存储的球形离 散单元信息，导入PFC3D 4.0软件获得集料颗粒的离散元模型。

步骤(2)中，轮廓检测及建模算法描述如下步骤：

1)、将图像Aggregate[i]的中心像素及其八邻域像素分别作为该图像所对应灰度差分八叉树GDO_i的根结点及根结点的子结点，并置此9个像素的访问标志bV＝true及根结点所有子结点的可扩展标志bE＝true后，计算GDO_i中所有连接线段的灰度差分值d_FC并将当前层号h置为2，最大深度H置为Aggregate[i]纵、横向分辨率中较小者一半的整数部分；

2)、使用灰度差分阈值自动分析算法获得灰度差分阈值T_G，依次判断每个位于GDO_i中第h层的结点与其父结点间连接线段的d_FC是否大于T_G，若是，则该结点的对应像素是轮廓像素，将该像素存入P_i并置该结点的bE＝false，否则仅置该结点的bE＝true；

3)、遍历GDO_i中第h层的结点若的bE＝true，则将对应像素的八邻域像素中bV＝false的像素作为的子结点，并置所有子结点的bE＝true、子结点对应像素的bV＝true，使h增1；

4)、持续运行步骤2)～步骤3)直至h＝H或者GDO_i中第h层上的结点所对应像素均为轮廓像素为止，完成构建GDO_i；

5)、依据P_i中各像素在图像Aggregate[i]中的邻接关系，对P_i中的像素进行分组排序，将在图像中依次邻接的像素按邻接顺序排列并归于一组，若P_i只包含一组像素，执行步骤7)，否则执行步骤6)；

6)、对P_i中每个像素组的两个末端像素，分别找出与之距离最短的另一像素组的末端像素与之配对，以配对像素之一为起点，另一配对像素为终点，以起点和终点间的连线方向为检测方向，将T_G逐次减1并检测轮廓像素，直至起点和 终点间的轮廓像素在两者连线方向上依次相邻；

7)、调用ACIS7.0的底层建模功能Spline_Fitting拟合P_i中的像素，获得FACE类型的轮廓面实体模型OF_i；

其中：

八邻域像素是指图像中一个像素的上、下、左、右、左上、右上、左下和右下八个方向上与之邻接的像素；

灰度差分八叉树GDO是指用以描述集料颗粒轮廓像素检测过程中所涉及像素间的邻接关系及邻接像素间的灰度差分的八叉树；

八叉树是指由若干结点按照结点间的父子关系以直线段连接而成、且每个结点最多有八个子结点的一种树形数据结构；

灰度差分值d_FC是指GDO内子结点与父结点对应像素间的灰度差；

当前层号h是指GDO_i在构建过程中的当前层数；

最大深度H是指所建GDO允许包含的最大层数；

轮廓面是指以P_i中像素的拟合结果为边界的平面。

轮廓检测及建模算法的步骤2)中，灰度差分阈值T_G自动分析算法描述如下步骤：

1)、将集料像素灰度范围上限赋值为Aggregate[i]所对应灰度直方图中的最大灰度值；

2)、将变量t赋值为

3)、让t每次递减1，直至中像素占Aggregate[i]中全部像素的比例超过集料像素比例p_AG；

4)、让t每次递减1，直至灰度为t的像素中不再包含集料像素为止，将集料像素灰度范围下限赋值为t；

5)、让t每次递减1，直至灰度为t的像素中不再包含沥青胶浆像素为止，将沥青胶浆像素灰度范围上限赋值为t；

6)、将灰度差分阈值T_G赋值为与的差值；

其中：

集料像素灰度范围上限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的上限；

灰度直方图是指相应图像中各灰度级上像素数量的统计图；

集料像素比例p_AG是指集料像素占Aggregate[i]全部像素的预估比例；

集料像素灰度范围下限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的下限；

沥青胶浆像素灰度范围上限是指Aggregate[i]所对应灰度直方图中沥青胶浆像素灰度范围的上限；

灰度差分阈值T_G是指用以判断两任意像素是否分别为集料和沥青胶浆像素的灰度差分最小值。

步骤(4)中，模型转换算法描述如下步骤：

1)、使用ACIS的底层建模功能Get_Bounding_Box获得M的包围盒B_M，以B_M的中心位置为中心，B_M最长边的长度为边长，建立BODY类型的立方体初始剖分空间IDS，并将其作为M所对应的转换八叉树TO的根结点root；

2)、使用ACIS的底层建模功能Mass_Proporties计算IDS与M的体积重合比例R_V，若R_V大于转换精度阈值T_A，那么将root的占有状态标志S赋值为占有状态fo，否则，若R_V＜1-T_A，则将S赋值为不占有状态cno；否则将S赋值为部分占有状态po；保存该信息到root，并标记root为未访问；

3)、将IDS剖分为8个等体积的立方体，作为TO中root的子结点Node₁,…,Node₈并标记为未访问，标记root为已访问；

4)、广度优先遍历TO，将每一个S＝fo或者S＝cno的Node标记为已访问，对每个未访问且S＝po的Node，若其对应立方体体积的1/8大于预设最小体积 V_min，将该立方体剖分为8个新的立方体，将新增立方体作为该Node的子结点并标记为未访问后，为各立方体对应的S赋值，并将该Node标记为已访问；

5)、持续执行步骤4)直至所有Node都标记为已访问或者Node对应立方体不用再度剖分，完成建立TO；

6)、遍历TO，对所有S＝fo的Node，建立相应立方体的内切球实体模型，获得M的球形离散单元形式；

其中：

包围盒B_M是指将集料颗粒实体模型M完全包围起来的最小长方体；

初始剖分空间IDS是中心位置与B_M相同、并能包含B_M的最小立方体；

转换八叉树TO是指用于描述剖分产生的各立方体之间相互覆盖关系和状态的八叉树结构；

体积重合比例R_V是指立方体与M重合部分的体积与立方体体积的比值；

转换精度阈值T_A是指预设用于判定立方体被M占有的体积重合比例R_V的最小值； 

占有状态标志S是指用于标识立方体与M体积重合程度的标志，包括占有状态fo，不占有状态cno和部分占有状态po；

广度优先遍历是指按照从上到下、从左到右的顺序依次遍历八叉树中结点的方式；

预设最小体积V_min是指预设用于判断剖分是否终止的体积最小值；

球形离散单元形式是指M经模型转换算法转换而成的由尺寸不一的球组成的实体模型；

球形离散单元信息是指DM{Sphere_k|k＝1,2,…,p}中每个球Sphere_k的球心和半径；

离散元模型是指可用于离散元数值仿真的计算模型。

本发明中，用C++语言，基于ACIS内核，实现了本发明所描述的算法，并 且以某沥青混合料试件的X-ray CT扫描图像为数据源，进行了真实集料颗粒离散元模型的构建。

(1)以直径为150mm，高度为164mm的圆柱形沥青混合料试件扫描步长为1mm的X-ray CT断面扫描图像文件组为数据源，选择如图1a所示的一个集料颗粒并截取如图1b所示的该集料颗粒在各扫描图像中的投影图块，保存为Aggregate[12]，用以重构该集料颗粒的实体模型；

(2)遍历灰度图像文件组Aggregate[12]，对于每幅图像，设定集料像素比例p_AG，最大深度H，并使用灰度差分阈值自动分析算法获得灰度差分阈值T_G，其中p_AG、T_G和H的值如表1所示，使用轮廓检测及建模算法，检测出如图2a所示的Aggregate[i](i＝1,2,…,12)中集料颗粒的轮廓像素，存入轮廓像素数组P_i(i＝1,2,…,12)，再建立如图2b所示的P_i中像素对应的轮廓面实体模型OF_i(i＝1,2,…,12)并存储。

表1各灰度图像文件对应的集料像素比例p_AG、最大深度H和灰度差分阈值T_G灰度图像文件Aggregate[i] 集料像素比例p_AG 最大深度H 灰度差分阈值T_G

(3)输入集料颗粒对应的全部轮廓面实体模型{OF_i|i＝1,2,…,12}，调用 ACIS7.0的底层建模功能Skinning，获得如图3所示的BODY类型的集料颗粒实体模型M；

(4)输入M，设定转换精度阈值T_A＝0.9，预设最小体积V_min＝8mm³，利用模型转换算法，获得如图4a所示的M以球形离散单元形式表示的实体模型DM{Sphere_k|k＝1,2,…,164}，其中Sphere_k(k＝1,2,…,164)对应的球心三维坐标及半径如表2所示，再获得以文本文件格式TXT存储的球形离散单元信息，导入PFC3D 4.0软件获得如图4b所示的集料颗粒离散元模型。

表2集料颗粒的球形离散单元形式实体模型中各球的球心坐标及半径

Claims (4)

1.基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：包括以下步骤：

(1)、以BMP格式、分辨率为512×512像素、扫描间隔一定的沥青混合料试件X-ray CT灰度图像文件组为数据源，选择一个集料颗粒并截取其在各扫描图像中的投影图块，并保存为BMP格式的灰度图像文件组Aggregate[n]，用以重构该集料颗粒的实体模型，其中：

所述沥青混合料试件是指以沥青材料和集料按适当比例配制，胶结成整体的圆柱形样本，记半径为radius，高度为height；

所述X-ray CT灰度图像文件组是指使用X-ray CT设备以一定纵向间隔扫描沥青混合料试件而获得的断面扫描图像数组，每幅断面扫描图像保存为BMP格式的灰度图像文件；

所述灰度图像文件是指将彩色图像中像素的颜色分量值按照一定比例进行加权，并将加权值附于每个颜色分量而获得的图像文件，其中，灰度图像中所有像素的灰度均介于0到255之间；

所述投影图块是指一幅扫描图像中包围目标集料颗粒图像的最小矩形内的图块；

所述实体模型是指以边界表示法表示的三维几何结构，记为模型M＝(FACE,EDGE,VERTEX,R)，其中FACE表示模型中面的集合，EDGE表示边的集合，VERTEX表示顶点的集合，R表示模型中各元素之间的邻接关系；

记录着面的类别：平面、样条面，以及面的曲面方程，其中face为FACE中的元素；

记录着边的曲线方程，其中edge为EDGE中的元素；

记录着顶点的坐标p(x,y,z)，其中vertex为VERTEX中的元素；

R＝{r₁,r₂}，其中：

r₁＝{(face₁,face₂,edge)|edge∈face₁∩face₂,edge∈EDGE,face_i∈FACE,i＝1,2}，r₁表示face₁和face₂相邻于边edge；

r₂＝{(edge₁,edge₂,vertex)|vertex∈edge₁∩edge₂,vertex∈VERTEX,edge_i∈EDGE,i＝1,2}，r₂表示边edge₁与edge₂相邻于点vertex；

(2)、遍历灰度图像文件组Aggregate[n]，使用轮廓检测及建模算法，检测出每幅图像Aggregate[i](i＝1,2,…,n)中集料颗粒的轮廓像素，存入轮廓像素数组P_i(i＝1,2,…,n)，再建立相应的轮廓面实体模型OF_i(i＝1,2,…,n)并存储；

(3)、输入集料颗粒所对应的全部轮廓面实体模型{OF_i|i＝1,2,…,n}，调用ACIS7.0的底层建模功能Skinning，获得BODY类型的集料颗粒实体模型M；

(4)、输入M，利用模型转换算法，获得M以球形离散单元形式表示的实体模型DM{Sphere_k|k＝1,2,…,p}，再获得以文本文件格式TXT存储的球形离散单元信息，导入PFC3D 4.0软件获得集料颗粒的离散元模型。

2.根据权利要求1所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述步骤(2)中，轮廓检测及建模算法描述如下步骤：

1)、将图像Aggregate[i]的中心像素及其八邻域像素分别作为该图像所对应灰度差分八叉树GDO_i的根结点及根结点的子结点，并置此9个像素的访问标志bV＝true及根结点所有子结点的可扩展标志bE＝true后，计算GDO_i中所有连接线段的灰度差分值d_FC并将当前层号h置为2，最大深度H置为Aggregate[i]纵、横向分辨率中较小者一半的整数部分；

2)、使用灰度差分阈值自动分析算法获得灰度差分阈值T_G，依次判断每个位于GDO_i中第h层的结点与其父结点间连接线段的d_FC是否大于T_G，若是，则该结点的对应像素是轮廓像素，将该像素存入P_i并置该结点的bE＝false，否则仅置该结点的bE＝true；

3)、遍历GDO_i中第h层的结点若的bE＝true，则将对应像素的八邻域像素中bV＝false的像素作为的子结点，并置所有子结点的bE＝true、子结点对应像素的bV＝true，使h增1；

4)、持续运行步骤2)～步骤3)直至h＝H或者GDO_i中第h层上的结点所对应像素均为轮廓像素为止，完成构建GDO_i；

5)、依据P_i中各像素在图像Aggregate[i]中的邻接关系，对P_i中的像素进行分组排序，将在图像中依次邻接的像素按邻接顺序排列并归于一组，若P_i只包含一组像素，执行步骤7)，否则执行步骤6)；

6)、对P_i中每个像素组的两个末端像素，分别找出与之距离最短的另一像素组的末端像素与之配对，以配对像素之一为起点，另一配对像素为终点，以起点和终点间的连线方向为检测方向，将T_G逐次减1并检测轮廓像素，直至起点和终点间的轮廓像素在两者连线方向上依次相邻；

7)、调用ACIS7.0的底层建模功能Spline_Fitting拟合P_i中的像素，获得FACE类型的轮廓面实体模型OF_i；

其中：

所述八邻域像素是指图像中一个像素的上、下、左、右、左上、右上、左下和右下八个方向上与之邻接的像素；

所述灰度差分八叉树GDO是指用以描述集料颗粒轮廓像素检测过程中所涉及像素间的邻接关系及邻接像素间的灰度差分的八叉树；

所述八叉树是指由若干结点按照结点间的父子关系以直线段连接而成、且每个结点最多有八个子结点的一种树形数据结构；

所述灰度差分值d_FC是指GDO内子结点与父结点对应像素间的灰度差；

所述当前层号h是指GDO_i在构建过程中的当前层数；

所述最大深度H是指所建GDO允许包含的最大层数；

所述轮廓面是指以P_i中像素的拟合结果为边界的平面。

3.根据权利要求2所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述轮廓检测及建模算法的步骤2)中，灰度差分阈值T_G自动分析算法描述如下步骤：

1)、将集料像素灰度范围上限赋值为Aggregate[i]所对应灰度直方图中的最大灰度值；

2)、将变量t赋值为

3)、让t每次递减1，直至中像素占Aggregate[i]中全部像素的比例超过集料像素比例p_AG；

4)、让t每次递减1，直至灰度为t的像素中不再包含集料像素为止，将集料像素灰度范围下限赋值为t；

5)、让t每次递减1，直至灰度为t的像素中不再包含沥青胶浆像素为止，将沥青胶浆像素灰度范围上限赋值为t；

6)、将灰度差分阈值T_G赋值为与的差值；

其中：

所述集料像素灰度范围上限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的上限；

所述灰度直方图是指相应图像中各灰度级上像素数量的统计图；

所述集料像素比例p_AG是指集料像素占Aggregate[i]全部像素的预估比例；

所述集料像素灰度范围下限是指Aggregate[i]所对应灰度直方图中集料像素灰度范围的下限；

所述沥青胶浆像素灰度范围上限是指Aggregate[i]所对应灰度直方图中沥青胶浆像素灰度范围的上限；

所述灰度差分阈值T_G是指用以判断两任意像素是否分别为集料和沥青胶浆像素的灰度差分最小值。

4.根据权利要求1所述的基于X-ray CT图像的真实集料颗粒离散元模型构建方法，其特征在于：所述步骤(4)中，模型转换算法描述如下步骤：

1)、使用ACIS的底层建模功能Get_Bounding_Box获得M的包围盒B_M，以B_M的中心位置为中心，B_M最长边的长度为边长，建立BODY类型的立方体初始剖分空间IDS，并将其作为M所对应的转换八叉树TO的根结点root；

2)、使用ACIS的底层建模功能Mass_Proporties计算IDS与M的体积重合比例R_V，若R_V大于转换精度阈值T_A，那么将root的占有状态标志S赋值为占有状态fo，否则，若R_V＜1-T_A，则将S赋值为不占有状态cno；否则将S赋值为部分占有状态po；保存该信息到root，并标记root为未访问；

3)、将IDS剖分为8个等体积的立方体，作为TO中root的子结点Node₁,…,Node₈并标记为未访问，标记root为已访问；

4)、广度优先遍历TO，将每一个S＝fo或者S＝cno的Node标记为已访问，对每个未访问且S＝po的Node，若其对应立方体体积的1/8大于预设最小体积V_min，将该立方体剖分为8个新的立方体，将新增立方体作为该Node的子结点并标记为未访问后，为各立方体对应的S赋值，并将该Node标记为已访问；

5)、持续执行步骤4)直至所有Node都标记为已访问或者Node对应立方体不用再度剖分，完成建立TO；

6)、遍历TO，对所有S＝fo的Node，建立相应立方体的内切球实体模型，获得M的球形离散单元形式；

其中：

所述包围盒B_M是指将集料颗粒实体模型M完全包围起来的最小长方体；

所述初始剖分空间IDS是中心位置与B_M相同、并能包含B_M的最小立方体；

所述转换八叉树TO是指用于描述剖分产生的各立方体之间相互覆盖关系和状态的八叉树结构；

所述体积重合比例R_V是指立方体与M重合部分的体积与立方体体积的比值；

所述转换精度阈值T_A是指预设用于判定立方体被M占有的体积重合比例R_V的最小值；

所述占有状态标志S是指用于标识立方体与M体积重合程度的标志，包括占有状态fo，不占有状态cno和部分占有状态po；

所述广度优先遍历是指按照从上到下、从左到右的顺序依次遍历八叉树中结点的方式；

所述预设最小体积V_min是指预设用于判断剖分是否终止的体积最小值；

所述球形离散单元形式是指M经模型转换算法转换而成的由尺寸不一的球组成的实体模型；

所述球形离散单元信息是指DM{Sphere_k|k＝1,2,…,p}中每个球Sphere_k的球心和半径；

所述离散元模型是指可用于离散元数值仿真的计算模型。