基于微观图象特征的水泥水化仿真及宏观属性预测的方法

摘要

本发明公开了基于微观图象特征的水泥水化仿真及宏观属性预测的方法，首先采用三维显微计算机断层扫描仪来获取水泥微观结构的演化同时测量其宏观属性。然后采用三维图象特征来描述水泥在每一龄期的微观结构，并利用遗传编程类算法直接建立图象特征与宏观属性间的以动力系统方式描述的动态关系。利用得到的动力系统来模拟水泥水化并预测宏观属性。本发明避免了对水泥的微观图象进行复杂的图象分割与物相分类，也不需要考虑物理及化学知识，从图象本身出发计算三维图象特征，并进一步得到图象特征与宏观属性之间的动态关系，避免了不正确的物相分析与图象分割带来的影响，从而提高对宏观属性的发展变化趋势的预测精度。

说明书

技术领域

本发明涉及水泥计算机仿真与预测领域，具体的讲，涉及基于微观图象特征的水泥水化 仿真及宏观属性预测的方法。

背景技术

水泥材料的水化反应是一种多物相、多尺寸、多孔隙的复杂的时变体系。物相是指物质 中具有特定的物理化学性质的组成部分，水泥同时含有硅酸三钙、硅酸二钙、铝酸三钙、铁 铝酸四钙、游离氧化钙等多种不同的成分。此外，水泥泥浆的组成单元从纳米级的水化硅酸 钙到毫米级的骨料，在不同尺度下其微观结构都不尽相同。水泥在与砂石、水等混合后形成 硬化体的水化、硬化过程中要经过多种物理和化学的变化，并形成具有一定强度的复杂的孔 隙结构。其内部机理的复杂性及随时间变化的特征使得高性能水泥材料的研究存在许多挑战 性的问题。由于在理论方法和设计工具方面的严重不足，使得作为世界上水泥产量最高的我 国，生产的高标号水泥及高性能混凝土所占的比例很低。为将高性能水泥材料的研究水平推 向新的高度，急需新的科学方法和高效、适用的分析设计工具。对水泥的水化反应过程进行 完整精确的仿真建模不但可以预测水泥性能的变化，解决其本身存在的问题，更可以辅助材 料学家设计出新型水泥材料，对高性能水泥材料的开发和产品质量的提高具有重大的科学价 值和广阔的应用前景。从上世纪后期开始，国际上的专家学者就开始为水泥的水化过程建立 模型来实现研究水泥微观结构与性能之间的关系进而达到改进水泥性能的目的。

Jennings与Johnson在1986年第69卷的美国陶瓷学会杂志上公开了一种水泥微观结构 仿真方法，论文名称为：水泥复合物微观结构演化的模拟(Simulation of microstructure  development during the hydration of a cement compound)。该方法用球型的硅酸三钙与水组成 的混合系统来表示微观结构的水化过程。反应的每一步中，硅酸三钙的颗粒逐步增大，并且 逐步生成水化硅酸钙与氢氧化钙。生成的水化硅酸钙环绕在硅酸三钙颗粒的四周，而氢氧化 钙则生长在孔隙空间中。该方法的问题是，在水化过程中的每一步，颗粒厚度的增长都是相 同的，并没有嵌入水化动力学方程。

荷兰代尔夫特理工大学TUD于1991年在其技术报告中公开了一种连续系统模型 HYMOSTRUC，报告名为：硬化水泥材料中的水化仿真与结构成形（Simulation of hydration and  formation of structure in hardening cement-based materials）。该系统将水泥颗粒模拟为立方体 空间参考单元内随机分布的球形颗粒，提出了一种水泥球型颗粒的体积随着水化过程的进行 逐渐增大并桥接的水化模型。经过多年的完善和发展，该模型已经可以模拟水泥材料因水化 引起的变形、强度、硬度变化等问题，并能够预测硬化水泥浆体的渗透性能，可以初步模拟 加入了石灰石等其它颗粒的水泥的水化过程。该方法的问题在于并没有给出水泥泥浆的微观 结构，水化产物也集中在单一颗粒上。

美国国家标准技术研究院NIST的Bentz在1997年第80卷的美国陶瓷学会杂志上公开了另 一种基于细胞自动机模型与数字图像的水泥模拟系统CEMHYD3D，论文名为：硅酸盐水泥水 化及微观结构演化的三维计算机仿真（Three-Dimensional Computer Simulation of Portland  Cement Hydration and Microstructure Development）。该方法始于水泥的粒径分布和背散射电 子图像，运用体视学原理重构三维初始微结构模型，然后以代表不同的水泥矿物相的像素为 基础，运用根据水化反应定义的细胞自动机规则操纵这些像素来实现对水泥水化不同时期三 维微观结构的模拟。利用细胞自动机逐象素点的控制水泥数字图象的演化过程，从而达到仿 真的目的。

美国国家标准技术研究院NIST的Bullard在2007年第15卷的材料科学与工程中的仿真与 建模上公开了一种称作HydratiCA的随机仿真方法，论文名为：一种水成材料系统中反应与运 输的三维微观结构模型（A three-dimensional microstructural model of reactions and transport in  aqueous mineral systems）。该模型引入了动力学演化的基本原理并于细胞自动机的方法进行 了结合，可以对包括水泥水化在内的多种矿物水化过程进行模拟。HydratiCA可以直接模拟物 相的增长与溶解，溶解过程中移动物质的扩散，溶解过程中或固体表面的物质之间的络合反 应以及成核现象。把多种不同的水化动力学过程集成到了单一的模拟框架之中。

瑞士洛桑联邦理工学院的Bishnoi和Scrivener在2009年第39卷国际水泥与混凝土研究杂 志上公开了一种称为μic的微观结构仿真平台，论文名为：一种新的水泥水化仿真平台μic （μic:a new platform for modelling the hydration of cements）。μic使用大量的多层逐步增长 的球型颗粒来表达水泥微观结构的演化，该方法提供了一种基于向量方法的全新实现。μic增 加了对反应条件和每一个颗粒周围邻域的考虑，真实的表示了水泥颗粒的分布情况。考虑到 水泥水化的许多方面仍然没有被很好的理解，该模型允许用户完全定制仿真过程。

当今存在多种仪器可用于观测水泥微观结构，如扫描电子显微镜，激光共聚焦显微镜， 三维显微计算机断层扫描仪等。此外，也有研究关注微观结构特性与硬化水泥的宏观物理属 性之间的关系，并为水泥的水化过程建立模型。常规来说，为了在水泥的微观结构与宏观性 能之间建立联系，微观图象的所有象素都需要被分类为不同的物相。然后再以分类结果为基 础，建立诸如体积比等微观描述，从而为计算宏观属性奠定基础。基于图象分割的结果，我 们也可以采用BENTZ公开的方法来模拟水泥微观结构的演化。然而，尽管根据物相分类与图象 分割的结果，人们可以利用有限元分析的手段来计算物理性能。然而由于受到解析度，部分 容积效应，噪音等因素的影响，精确的分类与分割是极难实现的。尤其是仅仅采用单一类型 的图象来可靠的区分物相类别更是困难。而分类的不可靠性会影响微观结构的描述，并进一 步影响宏观物理性能和水化仿真结果的计算。

目前国内外为实现水泥水化微观结构的仿真，以及进一步实现根据微观结构预测宏观属 性，都需要首先把微观图象的所有象素分类为不同的物相，然后采用细胞自动机等模型进行 微观结构的演化，而为了计算宏观属性，则需要根据物相分类的结果，采用有限元法等进行 计算。然而由于受到解析度，部分容积效应，噪音等因素的影响，精确的分类与分割是极难 实现的。尽管有人采用X射线衍射分析来进行物相分析并建立三维数据，但是该三维数据是 根据二维分布推测出来的三维结构，并非真实三维结构。尽管三维显微计算机断层扫描仪可 获取到真实三维微观结构，但是对其得到的三维图象来可靠的区分物相类别更是困难。分类 的不可靠性会影响微观结构的描述，并进一步影响宏观物理性能和水化仿真结果的计算。

发明内容

本发明要解决的技术问题是提供基于微观图象特征的水泥水化仿真及宏观属性预测的方 法，不需进行复杂的物相分析与图象分割，便可建立水泥真实三维微观结构与宏观属性之间 的动态关系并预测宏观属性的方法。该方法从图象本身出发直接描述微观与宏观之间的关系， 可充分利用三维图象的特性，避免物相分析与图象分割，提高宏观属性的预测精度。

为了实现上述目的，本发明采用如下技术方案：

基于微观图象特征的水泥水化仿真及宏观属性预测的方法，主要包括如下步骤：

步骤（1）：采用目标样品制作多份水泥泥浆样本并进行养护：其中一份水泥样本用于观 测，其他水泥样本用于测量宏观属性；

在指定的若干个龄期取出观测用水泥样本，并利用三维显微计算机断层扫描仪对观测用 水泥样本进行观测并记录水泥的微观结构三维图象数据；

每次在同一龄期，取出测量宏观属性用水泥样本进行测量，得到宏观属性数据；

将记录的水泥的微观结构三维图象数据与宏观属性数据，根据相同龄期一一对应的原则 建立对应关系，得到水泥的微观结构三维图象数据的演化数据与宏观属性数据的演化数据；

步骤（2）：针对步骤（1）得到的水泥的微观结构三维图象数据的演化数据，以其中任一 龄期的微观结构三维图象作为固定图象，其他龄期的微观结构三维图象作为浮动图象，将所 有其他龄期的微观结构三维图象根据固定图象执行三维刚体图象配准，得到经过三维刚体图 象配准的每一龄期的微观结构三维图象数据；

步骤（3）：针对步骤（2）中已得到的经过三维刚体图象配准的每一龄期的微观结构三维 图象数据，统计每幅微观结构三维图象的三维图象特征，建立起三维图象特征随时间演变的 变化数据；

步骤（4）：采用遗传编程类算法，根据步骤（1）得到的水泥宏观属性数据的演化数据和 步骤（3）得到的三维图象特征随时间演变的变化数据，逐步演化出图象特征与宏观属性动态 关系的动力系统；

步骤（5）：以建立的图象特征与宏观属性动态关系的动力系统为基础，能够预测微观结 构三维图象特征和宏观属性的变化趋势，输入初始水泥图象，统计其微观结构三维图象特征， 测量出初始宏观属性，将微观结构三维图象特征代入得到的动力系统，利用四阶龙格库塔法， 能够预测出微观结构三维图象特征和宏观属性在某一龄期的值及其变化趋势。

所述步骤（1）的宏观属性包括水泥的各种宏观物理特征：抗压强度、抗拉强度、抗折强 度。

所述步骤（3）的图象特征包括各种用于描述图象中象素值分布规律的特征，包括灰度直 方图、灰度共生矩阵、马尔可夫随机场、模式光谱。

所述步骤（4）的遗传编程类算法包括遗传编程家族的遗传编程、基因表达式编程、多表 达式编程等一类可进化程序与公式的算法。

所述步骤（1）包括如下步骤：

A1:指定水泥品种，指定水灰比，指定观测尺度，指定观测与测量龄期；

A2:制作1份观测用水泥样本，同时按照总观测天数制作多份测量用样本，对所有样本按 照相同的养护条件进行养护；

A3:在任一指定龄期，取出观测用水泥样本，利用三维显微计算机断层扫描仪获取样本的 三维微观结构图象，用后放回继续养护；在同一龄期，取出测量用样本测量宏观属性；

A4:将记录的真实三维微观结构数据与宏观属性数据根据相同龄期一一对应的原则建立 对应关系，得到三维微观结构的和宏观属性的演化数据。

所述步骤（2）包括如下步骤：

B1:从所有龄期的观测图象中，选择一幅水泥样本位置最接近图象中心的龄期的三维图象 作为固定图象，其余图象作为浮动图象；针对每一副浮动图象，执行B2至B4；

B2:设定初始变换参数，采用刚体变换，将每一副浮动图象变换到新的坐标系下，并采用 最近邻插值；

B3:计算平均均方误差来作为评价该次变换好坏的标准；

B4:若达到预定精度，则停止；否则，执行一步Powell算法，更新变换参数，并返回步骤 B2；

B5:得到配准后的在各个时间点的水泥微观结构三维图象。

所述步骤（4）包括如下步骤：

D1:设定函数集（+，-，*，÷等），设定终端结点集（包括所有选择的图象特征及需要预 测的宏观属性，以及系数）；

D2:随机生成初始化种群；

D3:为种群中的每个染色体计算适应值；首先将染色体解码为对应的表达式，然后采用四 阶龙格库塔法计算该表达式所对应的时间序列，进而采用下述公式计算该序列的适应值

$\mathrm{Fitness}=\frac{\underset{i=1}{\overset{L+1}{\Sigma }}\underset{j=1}{\overset{M}{\Sigma }}\left|\frac{x-\hat{x}}{x}\right|}{(L+1)M}+\alpha \frac{\underset{i=1}{\overset{L+1}{\Sigma }}l\left(s_i\right)}{(L+1)K}---\left(1\right)$

其中L表示图象特征的总数，M表示序列中的点的数目，即观测龄期总数。K表示每个 染色体的最大长度，x表示在某时间点的真实值，表示预测值，s_i表示染色体中第i个表达 式的子串，l(s_i)表示该子串解码出的表达式的实际长度，α表示公式中后半部分的权重。

D4:根据适应值的高低，采用遗传编程类算法的进化算子对种群中的染色体进行修改，产 生新的种群；进化算子包括交叉算子，变异算子，插串算子；

D5:若达到最大迭代次数，将适应值最佳个体解码，得到最终的动力系统并退出；否则转 向步骤D3。

所述步骤（5）包括如下步骤：

E1:针对同一水泥品种，输入指定的水灰比，指定观测尺度，指定观测与测量龄期；

E2:养护到初始龄期，利用三维显微计算机断层扫描仪对样品进行观测并记录初始三维 图象数据；

E3：根据初始三维图象数据统计三维图象特征；

E4：将三维图象特征作为初始变量，代入步骤（4）得到的动力系统，利用四阶龙格库塔 法，可预测图象特征和宏观属性到某一龄期的变化趋势。

本发明的有益效果：

1、建立了水泥微观结构的图象特征及其宏观属性之间的以动力系统方式描述的动态关 系，从而避免了不正确的物相分析与图象分割带来的影响，提高了对宏观属性发展变化趋势 的预测精度。

2、首先采用三维显微计算机断层扫描仪来无损的获取水泥水化过程中微观结构的演化同 时测量其宏观属性。然后采用三维图象特征来描述水泥在每一龄期的微观结构，并利用遗传 编程类算法直接建立图象特征与宏观属性间的动态关系。该动态关系以动力系统（微分方程 组）的方式来描述。利用得到的动力系统来模拟水泥水化并预测水泥的宏观属性。该方法避 免了对水泥的微观图象进行复杂的图象分割与物相分类，也不需要考虑物理及化学知识，只 需从图象本身出发计算三维图象特征，并进一步得到图象特征与宏观属性之间的动态关系， 避免了不正确的物相分析与图象分割带来的影响，从而提高对宏观属性的发展变化趋势的预 测精度。

附图说明

图1为本发明所采用的为样本建立图象特征动力系统的流程图；

图2为本发明所采用的根据图象特征动力系统预测微观结构演化和宏观属性的流程图；

图3为基于遗传编程类算法的进化水泥图象特征动力系统的流程图；

图4为基于遗传编程的进化水泥图象特征动力系统的流程图；

图5为基于基因表达式编程的进化水泥图象特征动力系统的流程图；

图6为基于多表达式编程的进化水泥图象特征动力系统的流程图；

图7为本发明的整体流程图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

如图1所示，本发明的建立水泥图像特征动力系统的步骤为：

步骤（101）：针对将要建模的水泥品种，预设好将要观测和测量的龄期序列，设定好水 灰比。制作一份观测用样品，然后制作测量宏观属性的样品N组，

步骤（102）：利用相同的养护条件对所有水泥样品进行养护；

步骤（103）：判断是否达到序列中某一预设龄期，如果是就进入步骤（104），如果否就 返回步骤（102）；

步骤（104）：取出观测用样品，采用三维显微计算机断层扫描仪进行三维图象数据的获 取；

步骤（105）：在同一龄期取出一组测量宏观属性用样品，测量并记录其宏观属性；

步骤（106）：判断是否达到预设序列的最终龄期，如果是就进入步骤（107），如果否就 返回步骤（102）；

步骤（107）：选定观测得到的序列三维图象中水泥样品位置最接近中心的龄期的图象为 固定图象，其他龄期图象为浮动图象，开始执行三维刚体图象配准。

步骤（108）：判断是否所有图象配准完毕，如果是就进入步骤（114），如果否就进入步 骤（109）；

步骤（109）：设定初始三维刚体变换参数；

步骤（110）：采用最近邻插值，执行三维刚体变换；

步骤（111）：计算浮动图象与固定图象之间的平均均方误差；

步骤（112）：判断是否达到期望精度，如果是就返回步骤（108）；如果否就进入步骤（113）；

步骤（113）：执行一步Powell算法；返回步骤（110）；

步骤（114）：统计每一龄期图象的图象特征，建立起图象特征随时间变化的时间序列数 据；

步骤（115）：根据图像特征随时间变化的时间序列和宏观属性序列，采用遗传编程类算 法，逐步进化图像特征的动力系统；

步骤（116）：将种群中的适应值最好的个体解码为最终得到的水泥图象特征动力系统。

如图2所示，根据图象特征动力系统预测微观结构演化和宏观属性的流程包括如下步骤：

步骤（201）：指定同一品种同一水灰比水泥；

步骤（202）：指定观测尺度、预测的龄期；

步骤（203）：制作水泥样品并养护到初始龄期；

步骤（204）：利用三维显微计算机断层扫描仪对样品进行扫描；

步骤（205）：提取三维图象特征；

步骤（206）：测量初始宏观属性；

步骤（207）：代入已得到的动力系统，采用龙格库塔法求解；

步骤（208）：得到从初始到指定龄期的图象特征及宏观属性的变化趋势。

如图3所示，所述步骤（115）的遗传编程类算法的通用步骤为：

步骤（301）：设定函数集，终端节点集；

步骤（302）：随机生成初始化种群；

步骤（303）：将种群中所有染色体解码，生成表达式；

步骤（304）：采用龙格库塔法得到仿真时间序列，并进一步计算所有可行动力系统的适 应值；

步骤（305）：执行遗传编程类算法中的进化算子；

步骤（306）：判断是否达到迭代次数，如果达到了就对种群中的适应值最好的个体解码 为最终得到的水泥图象特征动力系统；如果没有就返回步骤（303）。

如图4所示，如果所述步骤（115）中的遗传编程类算法为遗传编程时，所述步骤（115） 的基于遗传编程的进化水泥图象特征动力系统的流程如下：

步骤（401）：设定函数集、终端结点集、种群大小、迭代次数、进化算子参数；

步骤（402）：随机生成初始化种群，种群中的每个个体含有一组表达式树，表达式树个 数为方程个数，每棵表达式树对应于图象特征动力系统中的一个方程；

步骤（403）：对种群中的所有个体的所有表达式树进行遍历解码，得到每个个体对应的 动力系统；

步骤（404）：对每一个动力系统，采用龙格库塔法计算仿真时间序列。将该时间序列与 输入的真实图象特征序列及宏观属性序列进行比对，采用D3中的方程(1)计算适应值。该适 应值即为该动力系统所对应的个体的适应值；

步骤（405）：采用选择算子，从种群中按照适应值高低，生成新一代的种群；

步骤（406）：执行变异算子，按照一定的概率从表达式树中选择节点进行变异，从函数 集和终端节点集中随机选择元素进行替换；

步骤（407）：执行交叉算子，按照一定的概率，从种群中任选两个个体进行交叉，交叉 时对被选择的两个个体下一组表达式树中的对应的每棵表达式树都进行交叉，随机的从每课 树中选择交叉点，然后交换交叉点下面的子树；

步骤（408）：判断是否达到迭代次数，如果是就输出种群中适应值最优的个体，结束； 如果否就返回步骤（403）。

如图5所示，如果所述步骤（115）中的遗传编程类算法为基因表达式编程时，所述步骤 （115）的基于基因表达式编程的进化水泥图象特征动力系统的流程如下：

步骤（501）：设定函数集、终端结点集、种群大小、迭代次数、进化算子参数；

步骤（502）：随机生成初始化种群，种群中的每个个体为一个线性的超级染色体。超级 染色体包含多个按照从左至右顺序排列的线性的子染色体，个数为方程个数。每个子染色体 的编码规则为K表达式。每个子染色体对应于图象特征动力系统中的一个方程；

步骤（503）：对种群中的所有个体的所有超级染色体进行解码，得到每个个体对应的动 力系统；

步骤（504）：对每一个动力系统，采用龙格库塔法计算仿真时间序列。将该时间序列与 输入的真实图象特征序列及宏观属性序列进行比对，采用D3中的方程(1)计算适应值。该适 应值即为该动力系统所对应的个体的适应值；

步骤（505）：采用选择算子，从种群中按照适应值高低，生成新一代的种群；

步骤（506）：执行变异算子，按照一定的概率从超级染色体中选择位点进行变异，变异 值从函数集和终端节点集中随机选择元素进行替换；

步骤（507）：执行单点交叉算子，按照一定的概率，从种群中任选两个超级染色体进行 交叉。交叉时从两个超级染色体中随机选择一个相同的位点，然后互相交换位点之后的线性 字串；

步骤（508）：执行单点交叉算子，按照一定的概率，从种群中任选两个超级染色体进行 交叉。交叉时从两个超级染色体中随机选择两个相同的位点，然后互相两点之间的线性字串；

步骤（509）：执行插串算子，按照一定的概率，对种群中的每一个超级染色体的每一个 子染色体，随机选择两个位点，将一个位点后面的随机长度的字串插入到另一个位点中；

步骤（510）：执行根插串算子，按照一定的概率，对种群中的每一个超级染色体的每一 个子染色体，随机选择一个位点，将一个位点后面的随机长度的字串插入到其最开头的字串；

步骤（511）：判断是否达到迭代次数，如果是就输出种群中适应值最优的个体，结束， 如果否就返回步骤（503）。

如图6所示，如果所述步骤（115）中的遗传编程类算法为多表达式编程时，所述步骤（115） 的基于多表达式编程的进化水泥图象特征动力系统的流程如下：

步骤（601）：设定函数集、终端结点集、种群大小、迭代次数、进化算子参数；

步骤（602）：随机生成初始化种群，种群中的每个个体为包含有多个多基因表达式个体 的超级染色体。每个超级染色体包含有多个多基因表达式染色体(MEP染色体)，MEP染色体 数目为动力系统的方程个数。每个MEP染色体对应于图象特征动力系统中的一个方程；

步骤（603）：随机从种群选择两个个体；

步骤（604）：对两个个体中的所有MEP染色体，随机选择一个交叉位点，交换两个MEP 染色体中交叉位点之后的基因；

步骤（605）：对交叉后的新个体中的所有基因按照一定的概率执行变异操作；

步骤（606）：对变异之后两个新个体进行解码，对每一个个体的每一个基因都解码为动 力系统中的一个方程。采用龙格库塔法为所有方程计算仿真时间序列。将该时间序列与输入 的真实图象特征序列及宏观属性序列进行比对，采用D3中的方程(1)计算适应值。具有最高 适应值的基因所对应的方程所组成的动力系统为该个体的动力系统，该最高适应值为该个体 适应值；

步骤（607）：将两个新个体中适应值较高的个体的适应值与当前种群最坏的个体的适应 值进行比较，所小于最坏个体的适应值，则用新个体替换最坏的个体；

步骤（608）：判断是否达到迭代次数，如果是就输出种群中适应值最优的个体，结束， 如果否就返回步骤（603）。

如图7所示，本发明的整体流程图，步骤（116）的下一步就是步骤（201），图1得到了 图像特征动力系统，图2利用图1得到的图象特征动力系统预测微观结构演化和宏观属性。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限 制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付 出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。