用于多孔介质三维重构的模式密度函数模拟算法

本发明涉及利用数理统计和计算机建模方法从一张二维(2D)训练图像随机模拟多孔介质三维(3D)结构技术，属于图像处理技术领域。

3D结构有助于深入研究多孔介质的微观结构特性，但是由于多种原因，在实际生产中无法直接得到或者是难以得到多孔介质的3D结构。利用一张2D图像借助计算机技术和数学建模方法可以重建出具有相似特征的3D结构。

围绕此技术，多年来众多学者做了大量努力和探索，也诞生了多种较为经典的重建算法。按照重建方式，可以将目前产生的3D重建算法分为2D逐层重建和3D空间直接重建两种。

2D逐层重建的基本思想是将3D空间看作是由相邻2D图像序列构成。Tahmasebi和Hajizadeh所提的方法就是属于此类典型的算法。此类算法的关键点是如何控制相邻层之间的连续性和变化性，使之与符合训练图像的特征，为此Gao提出一种新颖的三步采样算法，能够很好地解决此问题。由于只是在某个面进行模式匹配，因而这类算法3D重建结果垂直重建方向上图像形态较差，一般需要后处理。

模拟退火算法SA和分级模拟退火算法HSA是3D空间直接重建方式中最为典型的一种，也属于最优化算法。该算法先将训练图像中的某个统计函数设定目标函数，然后按照训练图像孔隙度对重建的3D空间随机初始化，再不断随机选取不同相两点进行交换，使重建图像逐渐收敛于目标函数。为了避免在交换过程中出现局部最优情况，该算法引入金属冶炼过程中的退火机制。这类算法最大优点是稳定性，即当重建结果收敛于目标函数时重建结果中一定会呈现出目标函数所反映的图像特征。这种算法关键点是如何确定目标函数，目前模拟退火算法中一般是将自相关函数作为目标函数，但是实验证明基于低阶统计函数的自相关函数是无法准确描述复杂图像形态，这是由于满足同一自相关函数的图像形态可能是千差万别。

Okabe提出另一种3D空间直接重建算法，该算法首先从重建3D空间选取一条重建随机路径，然后利用单个标准方程模拟算法SNESIM分别确定3个中心过待模拟点且相互垂直模式，最后利用加权值法计算出待模拟点的值。此类算法在实际应用中仍然存在若干亟待解决的关键点：第一，应该依什么准则确定3个垂直切面的模式权值还需商榷；第二，在实际应用中SENSIM算法具有不稳定性。

基于纹理合成算法是另外一种在近年发展较为迅猛的按照3D空间直接重建算法。但是这类重建算法仅仅是在强调重建结果与原始训练图像在外观视觉上的相似性，而选择忽略两者统计特征上的一致性。

基于过程法亦属于3D空间直接重建算法，该算法模拟地质形成过程，但是复杂的地质形成过程和条件导致该类算法重建速度非常的漫长和难以控制。

当然也有其他重建算法，例如基于马尔科夫链算法、快速傅里叶变化算法、离散小波变换算法、神经网络算法等等。上述各种算法与3D精确重建要求都或多或少不足。一种精确重建应该包括三个方面要求：第一，3D重建结果与训练图像低阶统计特性相一致，如孔隙度，自相关函数，线性路径函数等；第二，3D重建结果（3个方向切面）与训练图像的形态特征相似；第三，算法具有稳定性，即在相同情况下重建结果具有较高的可重现性和对不同类型图像的普适性。本发明的目的在于为了解决上述现有技术中所存在的问题，为随机图像重建提供一种精确算法。

本发明引入概率论知识，在利用随机变量密度函数表征图像模式分布基础上原创性地提出模式密度函数概念，同时提出将模式密度函数作为迭代过程的目标函数的算法—模式密度函数模拟算法。并针对该重建算法中关键点速度问题提出邻域统计法和相邻网格反相选点法，并简化了温度控制机制。由于本发明所提出的模式密度函数属于高阶统计函数，当各级重建结果的模式密度函数趋于收敛于目标函数时不仅可以确保重建结果与训练图像的低阶统计特征一致，而且也可以保证形态特征相似。

本发明提出了一种利用二维训练图像随机重建具有相似特征三维结构方法，具体原理内容包括下述四个方面：

(1)提出模式密度函数概念和模式密度函数模拟算法

用C表示图像I中所关注的特征，例如形态特征、像素之间的位置关系等，这些特征是变化的，因此从概率论角度看C是一随机变量,ci，i=1,2…n表示特征C可能的取值。f(ci)表示图像I中C=ci时概率大小即

f(c i )=Pr{C=c i }i=1,2…n（1）

称f(ci)，i=1,2…n为特征密度函数，简记f(c)，其反映图像特征C的分布情况。图像随机重建是指利用数学建模的方法构建出与训练图像T具有相似特征分布的图像，但不是对T的拷贝。换句话说存在一个由具有相似特征分布图像构成的图像空间S={I}，T是代表S的一个样本，重建就是利用T反映的特征f(c)推测出这个图像空间S中其他可能样本I。而重建的关键点是如何挖掘与描述隐含在训练图像T中能够代表S共同特征C。

为了使重建结果与训练图像具有相似的形态特征分布，此时应选择f(c)表示图像形态特征分布，而图像的形态特征C就是模式。确定的图像具有确定的模式分布，不同图像空间的图像模式分布也不同，如附图1所示。此时特征密度函数f(c)即为模式密度函数（PatternsDensityFunction）PDF，此时PDF可以表示为

f(pat i )=Pr{Pat=pat i }i=1,2…n（2）

其中pati，i=1,2…n表示模式Pat的取值，同理模式密度函数f(pati)，i=1,2…n可简记为f(pat)。

用Tsize表示选择模板的大小，Npat表示模式变量Pat为具体模式pat时在训练图像T中出现的次数，Isize表示图像的大小，nmax表示最多的模式个数。f(pat)可以表述为：

$\[ f\left({\mathrm{pat}}_i\right)=\left\{\begin{array}{c}\mathrm{Pr}\{Pat={\mathrm{pat}}_i\}=N_{{\mathrm{pat}}_i}/I_{size}\\ \underset{i}{\Sigma }\mathrm{Pr}\{Pat={\mathrm{pat}}_i\}=1\\ n_{\max }=2^{T_{size}^2}\end{array}\right.,(i=1,2,\mathrm{...}{n}_{\max })---\left(3\right) \]$

模式的提取受模板尺寸大小控制。模板尺寸越大，获取模式越多，重建图像形态越精确，但同时重建图像所用时间会急剧增加；反之模板尺寸越小，获取模式越少，尺寸图像用时会降低，但是重建图像形态精度会降低。多级网格系统能够有效解决模板尺寸和重建图像形态之间矛盾问题。多级网格系统中各级网格之间满足下述关系

((T x -1)2 M-j +1)((T y -1)2 M-j +1)（4）

T x ×Ty是最小网格尺寸，M是多级网格的总级数，j(j=1,2,…M)是多级网格中某级网格的编号。在级数为M的多级网格形式下，f(pat)可以表述为:

$\[ f\left({\mathrm{pat}}_i^j\right)=\left\{\begin{array}{c}\mathrm{Pr}\{{\mathrm{Pat}}^j={\mathrm{pat}}_i^j\}=N_{{\mathrm{pat}}_i^j}/I_{size}^j\\ \underset{i}{\Sigma }\mathrm{Pr}\{{\mathrm{Pat}}^j={\mathrm{pat}}_i^j\}=1\\ n_{\max }=2^{T_{size}^2}\end{array}\right.,(i=1,2,\mathrm{...}{n}_{\max };j=1,2,\mathrm{...}M)---\left(5\right) \]$

相应地Isizej表示j级网格图的尺寸大小，Npatj表示第j级网格图中模式变量Patj为具体模式patij的出现次数。同一图像不同级网格图对应的模式密度函数也不尽相同，如附图2所示。

为了均衡重建的速度和精度，本发明选用模板尺寸为3×3大小。由于此模板尺寸较小，为了能够有效获取图像的形态特征，需要将训练图像分割到最小尺寸网格Gmin。Gmin是指将原始训练图像按照从大到小进行分级，当分级图中不再出现3×3大小全为孔隙相区域时对应分级图像的尺寸。此时说明利用3×3模板完全可以表达图像所有可能存在的模式。

本发明在提出模式密度函数概念和融合多级网格系统系统基础上，并提出将该函数作为迭代过程目标函数的随机结构重建算法，将其命名为模式密度函数模拟算法（PatternsDensityFunctionSimulation）PDFSIM。

(2)邻域统计法

在迭代过程中为了确定是否接受交换后的新状态，每交换一次，都需要计算重建图像与训练图像PDF的差值。以图像尺寸大小为128×128和3×3模板为例，如果采用全图统计，计算PDF需要128×128×3×3次，然后计算差值再需要循环512次，是否接受新状态都需要512次循环传递PDF。而一幅图像重建需要几百万次甚至更多次的迭代，因此使得PDFSIM实际不可操作。

幸运的是每交换两个点所引起变化的模式仅仅是这两个点对应3×3邻域18个模式（两个交换点的3×3邻域在重叠的情况下会更少），如附图3所示。因此本发明仅仅只计算交换点变化前后引起变化的函数值。此时统计PDF环节仅需要18×3×3次，同理计算差值最多也仅仅需要18次，新状态PDF传递最多也仅需要18次。理论分析和实验证明采用算法可以缩短重建时间数十倍，使得采用PDF为目标函数的迭代算法实际可操作。

（3）相邻网格邻域反相选点法

本发明提出相邻网格邻域反相选点法的目的是为了进一步加快PDFSIM重建算法速度。它的可行性是基于两点考虑。第一，随着迭代的深入不同相趋于聚合。第二，在多级网格机制中，前一网格层数据在后一层网格中作为条件数据，其极大地影响和决定周围像素值。

相邻网格反相选点法如图4所示，以上层网格已模拟点Pd为中心邻域D3内，若当前层布置的随机点Px的像素值与Pd相反，则将Px选为交换点。由于该方法排除了部分无效的待交换点，且随着迭代的深入孔隙不断聚合，有效交换点数越来越少，从而大大加快了收敛速度。

（4）简化温度控制机制

在PDFSIM算法中，需要根据每次交换前后重建图像的能量差值△E判断是否接受新状态。E1表示交换前的能量，E2表示交换后的能量，即

△E=E2-E1（6）

为了防止在迭代运算中出现局部最优，模拟退火算法引入温度控制机制，即采用Metropolis准则接受新状态，即：

（7）

T为重建结构当前的温度。但是退火机制极大地影响了重建的速度，为了进一步加快重建速度，本发明采用一种简化的温度控制机制,即仅仅对多级网格中最小网格采用温度控制,而对其余层网格放弃温度控制。即其余层对新状态接受准则为

（8）

两个理由能够支持此理论。第一，实验显示放弃温度控制而导致局部最优化现象容易发生在多级网格的最小网格。第二，在其他各级网格中，上级网格点作为条件数据能够有效阻止局部最优现象发生。实验显示放弃温度控制机制在大幅提升重建速度基础上，对随机图像重建的精度不会产生明显的影响。附图5显示为采用温度控制机制和简化温度控制机制重建对比效果。实验表面本发明采用的简化温度控制机制不会对随机图像重建的精度产生明显影响。

图1-1是本发明实施例中给定的某种电池材料二维图像；

图1-2是本发明实施例中给定的某种地质岩心二维图像；

图1-3是图1-1和图1-2对应模式密度函数分布图；

图2-1本发明实施例中给定的某种地质岩心二维原始图像；

图2-2是图2-1对应的1:2尺寸网格图像；

图2-3是图2-1对应的1:4尺寸网格图像；

图2-4是图2-1对应的1:8尺寸网格图像；

图2-5是图2-1、2-2、2-3、2-4对应的模式密度函数分布图；

图3是本发明实施例中邻域统计法示意图；

图4-1是本发明实施例中相邻网格反相选点法中前级网格分布图；

图4-2是本发明中相邻网格反相选点法中本级网格初始点分布图；

图4-3是本发明中相邻网格反相选点法中本级网格选择交换点分布图；

图5-1是本发明实施例中一个原始随机图；

图5-2是对图5-1采用原始温度控制机制重建图；

图5-3是对图5-1采用简化温度控制机制重建图；

图6-1是本发明实施例中一个原始随机图；

图6-2为采用本发明算法重建3D结构；

图6-3为6-2图的透视图；

图6-4为6-1图参考的micro-CT扫描的3D结构；

图6-5为6-4图的透视图；

图6-6为6-1、6-2、6-4的X方向两点相关函数比较图；

图6-7为6-1、6-2、6-4的Y方向两点相关函数比较图；

图6-8为6-2、6-4的Z方向两点相关函数比较图；

图6-9为6-1、6-2、6-4的X方向线性路径函数比较图；

图6-10为6-1、6-2、6-4的Y方向线性路径函数比较图；

图6-11为6-2、6-4的Z方向线性路径函数比较图；

图6-1为经过处理后二值化岩心图像，其中黑相表示岩心中的孔隙，白相表示岩心中的岩石。下面结合本案例对本发明具体实施过程作进一步详细说明，但所述实施案例仅对本发明实现方法作一个详细说明，而不应理解为是对本发明保护内容的任何限制。本发明具体实施过程包括下述几个步骤：

第一步，按照多级网格系统将原始图像从大到小分割至Gmin大小，成多级图像序列，并统计各级图像的模式密度函数f(patij)。

第二步，初始化3D重建结构。重建过程是按照网格图像从小到大的顺序进行。若重建的3D结构是最小网格，则按照对应原始网格图像的孔隙度随机初始化数据；若重建其他网格3D结构，则按照对应原始网格图像的孔隙度随机配置除上级网格点外其余网格点。同时计算此时3D结构的模式密度函数，并计算其与对应原始网格图像的模式密度函数之间的差值E1。

第三步，设置退出循环条件。该条件用来控制重建每个3D网格结构的循环退出条件。

第四步，选取交换网格点。若是重建最小网格3D结构，则将布置所有随机点作为交换点；若重建其他网格3D结构，则按照前述的相邻网格邻域反相选点法选取交换点。

第五步，从上步选取的交换点中随机选取两点进行交换，并按照邻域统计法计算交换两点后重建图像模式密度函数及其与对应原始网格图像的模式密度函数之间的差值E2。然后计算E1与E2之间的差值△E。

第六步，判断是否接受交换状态。若重建3D结构是最小网格，则按照模拟退火准则接受交换状态，即按照（9）式接受交换状态。若重建图像是其他网格图像，则按照简化退火准则接受交换状态，即按照（10）式接受交换状态。

第七步，判断是否结束某网格层结构的重建过程。若是没有到达退出循环条件，则继续交换选择的随机交换点，使重建结构的模式密度函数收敛于原始图像的模式密度函数。若是达到退出循环条件，结束本网格3D结构的重建，并将重建的结果扩展到下一层对应的网格点上，作为条件硬数据不再参与下层迭代。

第八步，在完成所有网格层重建后，完成重建内容。最终重建的3D结构如图6-2和图6-3所示。

为了展示本发明3D随机结构重建效果，本实施案例同时给出了micro-CT扫描岩心参照的3D结构，如图6-4和6-5所示。从视觉直观上看，利用本发明从一张2D训练图像重建的3D结构与参考的实际3D结构形态非常相似，从而验证了本发明的能够精确描述图像的形态结构。

为了进一步显示本发明的精确性，在该实施方案中利用随机图像重建中常用的两点相关函数和线性路径函数对重建的结果进行量化比较。比较结果如图6-6、6-7、6-8和6-9、6-10、6-11所示。由图可见，重建图像与原始图像的统计函数曲线完全吻合。从而说明采用本发明随机重建3D结构能够精确描述2D训练图像的统计特征。

上述实施例只是本发明的优选实施案例，并不是对本发明所述技术方案的限制，只要是不经过创造性劳动即可在上述实施案例的基础上实现的技术方案，均应视为落入本发明内容的保护范围内。