CN102768699B - 基于ct图像精确重构异质材料微观有限元网格模型的方法 - Google Patents

Abstract

本发明提出了一种利基于CT图像精确重构异质材料微观有限元网格模型的方法，该方法通过工业CT采集序列断层图像，并在数字化和阈值分割的基础上将断层图像中的微结构信息映射到重构的有限元网格模型中，使得异质材料中的任何细节结构信息都可以在重构模型中得到重现。本发明通过对比度受限自适应直方均衡化、中值滤波和像素点插值提高重构精度，通过图像裁剪和像素点合并提高重构效率。本发明直接重构出有限元分析精度较高的矩形（二维）和长方体（三维）单元网格模型，避免了现有几何重构方法在重构和网格划分等环节中的误差积累，提高了重构精度和效率。本发明可广泛应用于异质材料的性能预测及优化设计等领域。

Description

基于CT图像精确重构异质材料微观有限元网格模型的方法

技术领域

本发明属于有限元网格模型重构技术领域，涉及异质材料有限元网格模型重构方法，尤其是一种基于CT图像精确重构异质材料微观有限元网格模型的方法。

背景技术

由多种组分构成的异质材料（如复合材料、多孔材料等）是天然材料及合成材料中最为常见的材料之一，这些材料的宏观性能（如刚度、强度和韧度等）主要由其微观结构决定。因此，透彻地研究异质材料的微观结构对其宏观性能的影响对于设计和开发新的高性能异质材料具有重要意义。有限元方法是研究材料微观结构与宏观性能之间关系的最有效的方法之一，该方法需要首先建立能够反映材料真实微观结构的模型。然而，由于异质材料微观结构的复杂性，建立其微结构模型是材料设计的一大难题。

在异质材料性能的研究中，人们提出了多种克服建模困难的方法，其中，比较传统的是单胞模型方法。该方法基于异质材料在宏观上表现出的均质特性，将其微观结构理想化为具有简单几何特征的模型。这种方法能够有效地预测异质材料的组分性能、体积分数、形状和分布对宏观性能的影响，但是该方法忽略了异质材料真实的微观几何构造，难以用于对材料进行精确的分析。

随着数字图像处理技术的不断发展，研究人员提出了利用异质材料断层扫描图像重构其微观几何模型的方法。该方法采集断面图像的基础上，通过将图像分割为不同的组分相并进而提取各相边界信息的方法建立材料的几何表面或实体模型。该几何重构方法考虑了异质材料固有的形态学、聚集态和分布特征，并且尽可能少地引入微结构假设，这些特点使得该方法被广泛地用于分析异质材料的宏观行为。

尽管近年来几何重构方法趋于完善，对异质材料进行宏观性能的有限元预测仍然需要发展新的微观结构重构方法。因为采用现有几何重构方法建立的模型需要进一步划分网格，而基于复杂表面或实体划分的网格多为三角形或四面体单元，有限元分析精度不高；另外，异质材料微观结构的复杂性限制了网格的质量，且网格剖分过程中容易丢失材料微观结构中相对细小的信息。这些缺陷使得几何重构方法难以更加广泛地应用于材料的研究设计中，因而发展能够克服这些缺陷的新型重构方法显得非常重要。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种基于CT图像精确重构异质材料微观有限元网格模型的方法，该方法通过工业CT采集材料的序列断层图像，并将断层图像中的微结构信息映射到重构的有限元网格模型中，使得异质材料中的任何细节结构信息都可以在模型中得到重现；为了提高重构的精确度，引入了对比度受限自适应直方均衡化、中值滤波和像素点插值等图像处理方法；另外，提出了图像裁剪和像素点合并等方法，用于提高后期有限元分析效率。

本发明的目的是通过以下技术方案来解决的：

这种基于CT图像精确重构异质材料微观有限元网格模型的方法，包括以下步骤：

1）利用工业CT采集异质材料中各平行截面的序列断层图像，并将图像数字化；

2）对图像进行阈值分割，并通过对比度受限自适应直方均衡化、中值滤波和像素点插值提高图像反映材料真实微观结构的精度；

3）建立与序列断层图像形成映射关系的有限元网格拓扑模型，利用图像中的组分相灰度信息确定有限元模型中各个单元的材料属性，并通过图像裁剪和像素点合并来减小有限元网格模型的单元规模。

上述步骤1）具体按照以下方法进行：

通过工业CT机在试件表面选取一个宽为w，高为h的能够代表异质材料微观结构特征的代表性矩形区域I₁，采集该区域图像；然后将CT机的扫描位置沿试件的轴向平移微小距离d_n，其中n=1～D-1；并采集新的图像I_n+1，得到总层数为D的一系列断层图像；断层图像用像素点矩形表示，且矩阵沿i和j方向的大小分别记为W和H，故单个像素点的尺寸为s=w/W=h/H；像素点(i,j)的灰度用离散函数g(i,j)表示，故单张断层图像可以表示为灰度矩阵[g(i,j)](i＝1～W,j=1～H)；离散函数g(i,j,k)用于表示序列断层图像中像素点(i,j,k)的灰度，故序列断层图像表示为三维数组{g(i,j,k)}(k=1～D)。

上述步骤2）具体按照以下方法进行：

（1）图像阈值分割；引入统计图像中给定灰度值的像素点数目的灰度直方图h(g)：

$h\left(g\right)=\underset{1\≤i\≤W,1\≤j\≤H}{\mathrm{count}}\{g(i,j)=g\}---\left(8\right)$

对于具有M(M≥2)相的异质材料，灰度直方图具有M-1个波谷，将这些波谷对应的灰度值记为灰度阈值T_n(n=1,2,…,M-1)，则异质材料中各组分相由这一系列灰度阈值来分割，即图像中灰度值低于阈值T₁的像素点表示密度最低的组分相，灰度值在区间(T₁,T₂]内的像素点表示密度第二低的组分相，依此类推便可将断层图像完全地分割为M个组分相；各组分相的体积分数由以下累积分布函数f(g)来确定：

$f\left(g\right)=\frac{1}{\mathrm{WH}}\underset{n\≤g}{\mathrm{\Σ}}h\left(n\right)---\left(9\right)$

从而密度最低的组分相的体积分数为 密度第二低的组分相的体积分数为 依此类推便可确定全部M个组分相的体积分数；

（2）对比度受限自适应直方均衡化；采用对比度受限自适应直方均衡化降低环境及仪器误差对断层图像的影响，从而提高图像反映异质材料真实微观结构的精度；

（3）中值滤波；对断层图像进行中值滤波处理，以降低图像噪声；

（4）像素点插值；为了能够精确地重构材料沿z轴方向的微结构特征，序列图像的层间距要尽可能地接近图像像素点的尺寸，即层间距需要满足下式：

E(d)+σ(d)≤ξs其中， $\left\{\begin{array}{c}E\left(d\right)=\frac{1}{D}\underset{n=1}{\overset{D}{\mathrm{\Σ}}}{d}_n\\ \mathrm{\σ}\left(d\right)=\sqrt{\frac{1}{D}\underset{n=1}{\overset{D}{\mathrm{\Σ}}}{(d_n-E\left(d\right))}^2}\end{array}\right.---\left(10\right)$

参数ξ依照精度要求而定，推荐该值的取值范围为1.0≤ξ≤2.0；如果E(d)+σ(d)超过了上限ξs，序列图像便需要通过像素点插值的方法来提高重构精度；原始图像I_n和I_n+1之间插值的图像层数D_n以及新的序列图像之间的层间距 由下式确定：

将插值形成的新图像和原始图像按照他们的z坐标进行排序，则插值图像各像素点的灰度值由下式确定：

g(i,j,k)=P(Z_k,g(i,j,n)),k∈Ni,n∈Na

(11)

式中，Ni和Na是分别指向插值新图像和原始图像序号的集合，Z_k为新序列图像中第k张图像的z坐标，P(Z,g)是以(Z_n,g(i,j,n)),n∈Na为插值点的插值函数。

上述步骤3）具体按照以下方法进行：

（1）二维重构；建立与单层图像各像素点形成映射的二维有限元网格拓扑模型，该有限元模型由四节点的矩形单元构成，并且单元尺寸与图像像素点的尺寸s相等，模型中第n个节点的坐标由下式确定：

式中，“%”为整数除法取余， 取小于运算对象的最大整数；模型中第n个单元的四个节点为：

第n个单元的材料属性编号为：

$p_n=\left\{\begin{array}{cc}1,& g(i,j)\&le;T_1\\ r,& g(i,j)\&Element;(T_{r-1},T_r],1<r<M\\ M,& g(i,j)>T_{M-1}\end{array}\right.$ 其中，

（2）三维重构；建立与序列断层图像像素点形成映射关系的三维有限元网格模型，该模型由八节点的长方体单元构成，单元沿x和y方向的尺寸与图像像素点的尺寸s相等，沿z方向的尺寸与序列图像的层间距相等；模型中第n个节点的坐标为：

模型中第n个单元沿逆时针方向排列的八个节点由下式确定：

模型中第n个单元的材料属性编号由下式确定：

$p_n=\left\{\begin{array}{cc}1,& g(i,j,k)\&le;T_1\\ r,& g(i,j,k)\&Element;(T_{r-1},T_r],1<r<M\\ M,& g(i,j,k)>T_{M-1}\end{array}\right.$ 其中，

（3）图像裁剪及像素点合并；重构得到的模型中包含的单元数目非常大，使得后续有限元分析非常耗时；由于断层图像中靠近图像边缘的像素点受到的噪声干扰最多，故采用裁剪图像周界像素点的方式来减小重构模型的单元规模；然而，裁剪的像素点数目越多，图像中包含的材料聚集团簇的数目也就越少，使得该方法的应用具有一定的局限性；这个缺陷可以通过像素点合并来予以克服，即通过下式来合并图像中相邻的像素点：

上式中， 为合并之后新图像的灰度值，α和β分别为i和j方向合并的像素点数目；如果W和H不能分别被α和β整除，则原图像右侧和上侧缺失的像素点通过复制附近边界像素点的方式来予以填补。

本发明具有以下有益效果：

（1）本发明直接精确地重构出反映异质材料真实微观结构的有限元网格模型，避免了现有几何方法在重构和划分网格等环节中的误差积累，节约了重构时间，提高了重构效率；

（2）本发明采用了能够使断层图像更加精确地反映材料真实微观结构的图像处理方法，并且重构出的网格分别为有限元分析精度较高的矩形（二维）和长方体（三维）单元，相对于现有几何重构方法，大大提高了有限元分析的精度；

（3）本发明通过采用能够减小有限元网格模型单元规模的图像处理方法，有效地平衡了有限元分析精度与计算时间之间的矛盾，在保证分析精度的前提下提高了计算效率；

（4）本发明重构的有限元网格模型可以方便地用于分析异质材料中各组分相的形态、聚集特征以及分布等微观结构与宏观性能之间的关系，从而为材料的开发设计提供依据。

附图说明

图1为本发明的流程图；

图2为序列断层图像的采集；

图3为图像的数字化；

图4为与单张断层图像形成映射关系的二维有限元网格模型的建立；

图5为二维有限元网格模型的节点和单元编号；

图6为与序列断层图像形成映射关系的三维有限元网格模型的建立；

图7为三维有限元网格模型的节点和单元编号；

图8为某两相异质材料的单张断层图像；

图9为图8中图像的灰度直方图和累积分布函数；

图10为经过CLAHE和MF处理后的单张断层图像；

图11为图10中图像的灰度直方图和累积分布函数；

图12为对某两相异质材料重构出的二维有限元网格模型；

图13为对图10中图像进行的裁剪处理及其阈值分割结果；

图14为对图10中图像进行的像素点合并处理及其阈值分割结果；

图15为经图像裁剪和像素点合并处理之后重构得到的二维有限元网格模型；

图16为某两相异质材料的序列断层图像；

图17为图16中序列断层图像的灰度直方图和累积分布函数；

图18为对某两相异质材料重构出的三维有限元网格模型；

图19为经像素点合并处理之后重构得到的三维有限元网格模型。

具体实施方式

本发明首先通过工业CT采集异质材料的序列断层图像，并将采集的图像数字化；然后通过灰度直方图和累积分布函数对图像进行阈值分割，并通过对比度受限自适应直方均衡化、中值滤波和像素点插值提高图像反映材料真实微观结构的精度；进而在建立有限元网格拓扑模型的基础上，通过拓扑模型与断层图像之间的映射关系确定模型的单元材料属性，并采用图像裁剪和像素点合并减小有限元网格模型中的单元规模。整个有限元网格模型重构的具体实施流程如图1所示，下面根据该流程详细地描述具体技术问题。

1.序列图像采集及数字化

本发明通过工业CT采集异质材料的序列断层图像，采集过程如图2所示。首先将打磨抛光后的异质材料试件固定在工业CT机的转台上，调节转台及平移系统，在试件表面选取一个宽为w，高为h的能够代表异质材料微观结构特征的代表性矩形区域I₁。该区域的选择不仅要使采集的图像中包含足够多的相聚集簇，同时又要限制图像中相聚集簇的数目以保证较高的图像分辨率和较快的重构速度。确定代表性矩形区域I₁后，用CT机采集该区域的图像，然后将CT机的扫描位置沿试件的轴向平移微小距离d₁（注意，在扫描位置的平移过程中必须保证在x和y方向不产生位移），采集新的矩形区域I₂的图像，重复以上过程便可采集总层数为D的一系列断层图像。上述过程中，相邻两层图像之间的距离d_n(n=1～D-1)要尽可能接近图像像素点的尺寸，以使序列断层图像能够精确地反映材料在z方向的微结构特征。另外，层数D的选择需要保证由序列图像重构出的模型在z方向包含足够数目的组分相聚集簇。

每一张断层图像都是以矩形排列的一个个像素点构成的，因此这些图像可以在图3所示的像素点坐标中以矩阵的形式予以表示。如图，图像矩阵沿i和j方向的大小分别记为W和H，每个像素点在宽度和高度方向的尺寸相等，并记为s=w/W=h/H。在每个像素点中，存储着反映图像视觉效果的信息，这些信息由名为灰度的整数表示，因此像素点(i,j)的灰度可以用离散函数g(i,j)予以表示，单张断层图像可以表示为灰度矩阵[g(i,j)](i=1～W,j=1～H)。对于序列断层图像，我们建立像素点左手坐标系下的k坐标，用以记录这一系列图像的序号，则离散函数g(i,j,k)可用于表示像素点(i,j,k)的灰度，三维数组{g(i,j,k)}(k=1～D)可用于表示序列断层图像。

2.图像阈值分割及图像处理

（1）图像阈值分割

工业CT采集的断层图像中各像素点的灰度值主要由异质材料各组分相的密度决定，即组分相密度越高，其吸收的X射线越多，该相占据的图像区域中的灰度值也相应地越高。根据这一规律，图像中的不同组分相可以由图像的灰度值予以区分。为了分割图像中不同的组分相，引入灰度直方图h(g)，该直方图通过统计图像中给定灰度值的像素点数目来确定，即：

$h\left(g\right)=\underset{1\&le;i\&le;W,1\&le;j\&le;H}{\mathrm{count}}\{g(i,j)=g\}---\left(18\right)$

对所有的灰度值进行遍历统计便得到图像的灰度直方图。

对于具有M（M≥2）相的异质材料，其不同的组分相具有不同的密度，从而在CT图像中由不同的灰度值表征。灰度直方图具有M-1个波谷，这些波谷便是各相材料的分界点。将M-1个波谷处对应的灰度值记为灰度阈值T_n(n=1,2,…,M-1)，则异质材料中各组分相由这一系列灰度阈值来分割，即图像中灰度值低于阈值T₁的像素点表示密度最低的组分相，灰度值在区间(T₁,T₂]内的像素点表示密度第二低的组分相，依此类推便可将断层图像完全地分割为M个组分相。

将断层图像分割之后，各组分相的体积分数便可以由累积分布函数f(g)来确定。累积分布函数由以下经过像素点总数目归一化处理的公式计算，即：

$f\left(g\right)=\frac{1}{\mathrm{WH}}\underset{n\&le;g}{\mathrm{\&Sigma;}}h\left(n\right)---\left(19\right)$

利用该累积分布函数，可以确定密度最低的组分相的体积分数为 密度第二低的组分相的体积分数为 依此类推便可确定全部M个组分相的体积分数。

（2）对比度受限自适应直方均衡化

断层图像的采集容易受到环境及仪器误差等一系列因素的影响，因此采集到的图像不能精确地反映异质材料真实的微观结构，而这一缺陷可以通过对比度增强和降噪等图像处理方法得以改善。

对比度增强有助于重现异质材料的真实微观结构，因为该方法能够增强图像中组分材料之间的对比度。本发明采用对比度增强系列方法中应用最为广泛的对比度受限自适应直方均衡化（CLAHE），该方法的数学描述详见Reza的综述文章（A.M.Reza,J VLSI Sig Proc Syst,38(2004)35-44）。

（3）中值滤波

除了对比度增强，降噪也能提高图像反映材料真实微观结构的精度。图像降噪方法很多，包括均值滤波、中值滤波和自适应滤波等，本发明选用其中降噪效果较好的中值滤波（MF）。对于特定形式的随机噪声，中值滤波能够在尽可能少地模糊边界和其他细节的情况下提供优异的降噪性能。中值滤波的理论和实现方法详见Gonzalez关于图像处理的书籍（R.C.Gonzalez,R.E.Woods,Digital Image Processing,Second Ed.,Prentice Hall,New Jersey,2001）。

（4）像素点插值

对于序列断层图像，为了能够精确地重构材料沿z轴方向的微结构特征，相邻图像之间的距离要尽可能地接近图像像素点的尺寸，也就是说，层间距需要满足下式：

E(d)+σ(d)≤ξs (20)

式中，E(d)和σ(d)分别为层间距的数学期望和标准差，即

$\left\{\begin{array}{c}E\left(d\right)=\frac{1}{D}\underset{n=1}{\overset{D}{\mathrm{\&Sigma;}}}{d}_n\\ \mathrm{\&sigma;}\left(d\right)=\sqrt{\frac{1}{D}\underset{n=1}{\overset{D}{\mathrm{\&Sigma;}}}{(d_n-E\left(d\right))}^2}\end{array}\right.---\left(21\right)$

参数ξ依照具体的精度要求而定，但是本发明推荐该值的取值范围为1.0≤ξ≤2.0。

公式(22)表明只要大多数的层间距不超过规定的上限ξs，序列图像便可直接用于重构异质材料的微观模型。相反，如果E(d)+σ(d)超过了上限ξs，序列图像便需要通过像素点插值的方法来提高重构精度。在像素点插值的过程中，原始图像之间将产生一系列新的图像，而这些图像用于重现原始图像之间的微结构细节。原始图像I_n和I_n+1之间插值的图像层数D_n以及新的序列图像之间的层间距 由下式确定：

插值形成的新图像和原始图像需要按照他们的z坐标进行重新排序，并且新序列图像的记号与原始图像保持一致。插值图像各像素点的灰度值由以下插值公式确定：

g(i,j,k)=P(Z_k,g(i,j,n)),k∈Ni,n∈Na (24)

式中，Ni和Na是分别指向插值新图像和原始图像序号的集合，Z_k为新序列图像中第k张图像的z坐标，P(Z,g)是以(Z_n,g(i,j,n)),n∈Na为插值点的插值函数。

3.有限元网格模型重构

（1）二维重构

本发明可以根据序列断层图像中的任何一张图像I_n来重构异质材料的二维真实微观结构。建立与单层图像各像素点形成映射的二维有限元网格拓扑模型，如图4所示。该有限元模型由四节点的矩形单元构成，并且单元沿x和y方向的尺寸与图像像素点的尺寸s相等。为了便于建立有限元模型与图像像素点之间简单的数学映射关系，对有限元模型节点和单元编号如图5所示，则模型第n个节点的坐标由下式确定：

式中，x_n和y_n分别为节点n沿x和y方向的坐标，“%”为整数除法取余， 取小于运算对象的最大整数。模型第n个单元逆时针方向排列的四个节点由下式确定：

在有限元数据文件中，不同的单元材料属性由其属性编号确定，且属性编号用从1开始的自然数予以表示，编号越小，代表单元的材料密度越小。从而根据阈值分割后的断层图像和映射关系可确定第n个单元的材料属性编号为：

$p_n=\left\{\begin{array}{cc}1,& g(i,j)\&le;T_1\\ r,& g(i,j)\&Element;(T_{r-1},T_r],1<r<M\\ M,& g(i,j)>T_{M-1}\end{array}\right.---\left(27\right)$

其中，i和j与n之间的关系为：

（2）三维重构

建立与序列断层图像像素点形成映射关系的三维有限元网格模型，如图6所示。该有限元网格模型由八节点的长方体单元构成，单元沿x和y方向的尺寸与图像像素点的尺寸s相等，沿z方向的尺寸与序列图像的层间距相等。对有限元网格模型的节点和单元编号如图7所示，则模型第n个节点的坐标为：

模型中第n个单元沿逆时针方向排列的八个节点由下式确定：

与二维重构类似，模型中第n个单元的材料属性编号可由下式确定：

$p_n=\left\{\begin{array}{cc}1,& g(i,j,k)\&le;T_1\\ r,& g(i,j,k)\&Element;(T_{r-1},T_r],1<r<M\\ M,& g(i,j,k)>T_{M-1}\end{array}\right.---\left(31\right)$

其中，i、j和k与n之间的关系为

（3）图像裁剪及像素点合并

根据以上二维和三维重构方法建立的有限元网格模型能够保留异质材料详细而真实的微结构特征，但是，重构的模型中包含的单元数目非常大，而且随着图像分辨率的提高，模型中将包含更多的单元，这将使得后续有限元分析非常耗时。由于断层图像中靠近图像边缘的像素点受到的噪声干扰最多，故可以通过裁剪图像周界像素点的方式来减小重构模型的单元规模。

然而，裁剪的像素点数目越多，图像中包含的材料聚集团簇的数目也就越少，使得该方法的应用具有一定的局限性。这个缺陷可以通过像素点合并来予以克服，即通过下式来合并图像中相邻的像素点：

上式中， 为合并之后新图像的灰度值，α和β分别为i和j方向合并的像素点数目。如果W和H不能分别被α和β整除，则原图像右侧和上侧缺失的像素点通过复制附近边界像素点的方式来予以填补。

4.实施例

（1）二维重构

首先利用工业CT机采集如图8所示的某种两相异质材料的单张断层图像，其基体相材料的密度低于增强相材料的密度，故图像中趋近黑色的部分为基体相材料，趋近白色的为增强相材料。该图像i和j方向的像素点数目分别为W=606和H=597，x和y方向的实际尺寸分别为w=16.12μm和h=15.88μm，从而每个像素点的尺寸为s=0.0266μm。图9给出了该图像的灰度直方图和累积分布函数，由图可知，灰度阈值T₁=154，基体相的体积分数为 实验测得材料基体相的体积分数为55.00%，即通过图像分割得到的结果具有2.64%的相对误差。为了减小建模误差，对断层图像实施以下四种图像处理：

（a）CLAHE处理，处理后T₁=140，v₁=54.92%，相对误差为0.08%；

（b）MF处理，处理后T₁=154，v₁=53.63%，相对误差为2.49%；

（c）先CLAHE后MF处理，处理后T₁=140，v₁=54.98%，相对误差为0.04%；

（d）先MF后CLAHE处理，处理后T₁=140，v₁=54.54%，相对误差为0.84%。

由以上这一系列图像处理结果来看，先对断层图像进行CLAHE处理，然后进行MF处理时，图像能够最为精确地反映异质材料的真实微观结构。经过这种处理之后的断层图像如图10所示，其灰度直方图和累积分布函数如图11所示。利用本发明提出的二维重构方法对处理后的图像进行重构，得到如图12所示的二维有限元网格模型。

由于重构的有限元模型中单元数目过多，对图像进行进一步的图像裁剪和图像合并处理。首先，对图像实施如图13所示的一系列图像裁剪处理，由结果可知，处理后的图像的相对误差均不超过2.79%。图像裁剪的程度可以在平衡重构精度和计算效率的基础上予以确定，如果要求相对误差不超过1.00%，那么较好平衡计算效率的处理方法是裁减掉图像28个像素点的周界，在这种情况下，图像的相对误差为0.68%。

图像裁剪之后单元数目仍然过多，即单一的图像裁剪无法适当地减小模型单元规模。尝试对图像开展如图14所示的一系列像素点合并处理，由结果可知，处理后的图像相对误差均不超过1.94%，仍然可以通过衡量重构精度和计算效率来确定像素点合并的数目。

经过一系列试凑发现，首先裁剪掉图像28个像素点的周界，然后对图像进行5个像素点的合并处理能够较好地平衡重构精度和计算效率。经过处理后W=110，H=108，T₁=140，v₁=54.81%，相对误差为0.352%。由该图像重构得到的二维有限元网格模型如图15所示。

（2）三维重构

利用工业CT机采集如图16所示的上述两相异质材料的序列断层图像，其W=H=256，w=h=6.8096μm，s=0.0266μm，D=11，d_n=0.158μm。为了提高各层图像重现材料微观结构的精度，首先对序列图像进行CLAHE处理，然后对其进行MF处理。处理之后的图像的灰度直方图和累积分布函数如图17所示，由图可知，T₁=139，v₁=54.66%，相对误差为0.62%。

取参数ξ为1.0，则E(d)+σ(d)超过了上限ξs=0.0266μm，故先对序列图像进行自然样条插值，得到新的序列图像中包含61层图像，且新序列图像相邻两层之间的距离为0.0263μm，这足以重构材料沿z方向的微结构细节信息。由新序列图像重构得到的三维有限元网格模型如图18所示，重构时T₁=139，v₁=55.21%%，相对误差为0.38%。

重构得到的三维模型中包含过多的单元，需要对其进行处理以减小单元规模。由于该序列图像沿i和j方向的像素点数目不多，这里只对其进行像素点合并处理，且处理过程中沿两个方向只合并两个像素点。经过图像处理之后重构出的三维有限元网格模型如图19所示，其相对误差为0.49%。

上述建模重构实例是针对两相异质材料而进行的，但是本发明所提出的技术并不仅限于两相异质材料，由前文述及的技术原理可知，本发明可以容易地对任意多相的异质材料进行有限元网格重构。

Claims (1)

1.一种基于CT图像精确重构异质材料微观有限元网格模型的方法，其特征在于，包括以下步骤：

1)利用工业CT采集异质材料中各平行截面的序列断层图像，并将图像数字化；步骤1)具体按照以下方法进行：

通过工业CT机在试件表面选取一个宽为w，高为h的能够代表异质材料微观结构特征的代表性矩形区域I₁，采集该区域图像；然后将CT机的扫描位置沿试件的轴向平移微小距离d_n，其中n＝1～D-1；并采集新的图像I_n+1，得到总层数为D的一系列断层图像；断层图像用像素点矩阵表示，且像素点矩阵沿i和j方向的大小分别记为W和H，故单个像素点的尺寸为s＝w/W＝h/H；像素点(i,j)的灰度用离散函数g(i,j)表示，故单张断层图像表示为像素点灰度矩阵[g(i,j)]，其中i＝1:W,j＝1:H；离散函数g(i,j,k)用于表示序列断层图像中像素点(i,j,k)的灰度，故序列断层图像表示为三维数组{g(i,j,k)}，其中k＝1～D；

2)对图像进行阈值分割，并通过对比度受限自适应直方均衡化、中值滤波和像素点插值提高图像反映材料真实微观结构的精度；步骤2)按照以下方法进行：

(1)图像阈值分割；引入统计图像中给定灰度值的像素点数目的灰度直方图h(g)：h(g)是灰度g的函数：

$h\left(g\right)=\underset{1\&le;i\&le;W,1\&le;j\&le;H}{count}\{g(i,j)=g\}---\left(1\right)$

对于具有M相的异质材料，其中M≥2，灰度直方图具有M-1个波谷，将这些波谷对应的灰度值记为灰度阈值T_n，n＝1,2,…,M-1，则异质材料中各组分相由这一系列灰度阈值来分割，即图像中灰度值低于阈值T₁的像素点表示密度最低的组分相，灰度值在区间(T₁,T₂]内的像素点表示密度第二低的组分相，依此类推便可将断层图像完全地分割为M个组分相；各组分相的体积分数由以下累积分布函数f(g)来确定：

$f\left(g\right)=\frac{1}{WH}\underset{n\&le;g}{\&Sigma;}h\left(n\right)---\left(2\right)$

其中，n为原始图像序号；g为灰度；

从而密度最低的组分相的体积分数为密度第二低的组分相的体积分数为依此类推便可确定全部M个组分相的体积分数；

(2)对比度受限自适应直方均衡化；采用对比度受限自适应直方均衡化降低环境及仪器误差对断层图像的影响；

(3)中值滤波；对断层图像进行中值滤波处理，以降低图像噪声；

(4)像素点插值；序列图像的层间距需要满足下式：

E(d)+σ(d)≤ξs其中，

E(d)和σ(d)分别为层间距的数学期望和标准差；s为图像像素点的尺寸；参数ξ依照精度要求而定；如果E(d)+σ(d)超过了上限ξs，序列图像便需要通过像素点插值的方法来提高重构精度；原始图像I_n和I_n+1之间插值的图像层数D_n以及新的序列图像之间的层间距由下式确定：

将插值形成的新图像和原始图像按照他们的z坐标进行排序，则插值图像各像素点的灰度值由下式确定：

g(i,j,k)＝P(Z_k,g(i,j,n)),k∈Ni,n∈Na (5)

式中，Ni和Na是分别指向插值新图像和原始图像序号的集合，Z_k为新序列图像中第k张图像的z坐标，P(Z)是以(Z_k,g(i,j,n)),k∈Ni，n∈Na为插值点的插值函数，即：

P(Z)＝Interpolation(Z_k,g(i,j,n)),k∈Ni,n∈Na (6)

3)建立与序列断层图像形成映射关系的有限元网格拓扑模型，利用图像中的组分相灰度信息确定有限元模型中各个单元的材料属性，并通过图像裁剪和像素点合并来减小有限元网格模型的单元规模，步骤3)按照以下方法进行：

(1)二维重构；建立与单层图像各像素点形成映射的二维有限元网格拓扑模型，该有限元模型由四节点的矩形单元构成，并且单元尺寸与图像像素点的尺寸s相等，模型中第n个节点的坐标由下式确定：

式中，“％”为整数除法取余，取小于运算对象的最大整数；模型中第n个单元的四个节点为：

第n个单元的材料属性编号为：

其中，

(2)三维重构；建立与序列断层图像像素点形成映射关系的三维有限元网格模型，该模型由八节点的长方体单元构成，单元沿x和y方向的尺寸与图像像素点的尺寸s相等，沿z方向的尺寸与序列图像的层间距相等；模型中第n个节点的坐标为：

模型中第n个单元沿逆时针方向排列的八个节点由下式确定：

模型中第n个单元的材料属性编号由下式确定：

其中，

(3)图像裁剪及像素点合并；采用裁剪图像周界像素点的方式来减小重构模型的单元规模；即通过下式来合并图像中相邻的像素点：

上式中，为合并之后新图像的灰度值，α和β分别为i和j方向合并的像素点数目；如果W和H不能分别被α和β整除，则原图像右侧和上侧缺失的像素点通过复制附近边界像素点的方式来予以填补。