基于多线程的岩心图像超维重建快速算法

摘要： 针对基于邻域块匹配的超维算法在重建过程中每次只能对1 个待重建块进行重建，并且字典搜索过程十分耗时导致重建效率低的问题，本文提出了2 种方法来对基于邻域块匹配的超维算法进行计算速度的优化. 首先，将分区域重建应用到基于邻域块匹配的超维算法中，提出了分区域并行重建算法，实现了对多个待重建块同时进行重建；其次，在对每次字典元素搜索的过程中，使用了字典并行搜索，实现了这一过程的加速；最后将这2 种方法进行结合，并且通过高中低3 种不同孔隙度的训练图像生成的字典来对二维参考图像进行多次重建.通过将本文提出的算法、基于邻域块匹配的超维算法以及一些传统重建算法的重建结果和真实岩心三维结构的统计特征参数进行对比并且将不同重建算法的重建时间进行对比，来验证本文改进的超维算法的有效性.

关键词： 三维重建； 超维重建； 多孔介质； 并行算法； 多线程

中图分类号： TP391 文献标志码： A DOI： 10. 19907/j. 0490-6756. 2024. 043002

1 引言

在数字岩心重建领域中，利用二维图像进行三维岩心重建是热点研究方向. 重建算法主要分为传统算法和基于学习的算法. 传统算法有模拟退火［1］（ Simulated Annealing，SA）和多点地质统计［2］（Multiple-point Geostatistics，MPS）以及相对应的改进算法［3］. 在基于学习的算法的领域中，Mosser 等［6］使用了对抗网络模型［7］对砂岩进行了重建，Feng 等［8］首次提出了BicycleGan［9］这一网络模型的二维到三维的重建. He 等基于学习和模式匹配的概念，提出了超维［10］（Super-Dimension，SD）的重建算法. 该算法成功地将相似的三维结构作为先验信息来引导三维重建. SD 算法包括选择训练图像建立字典以及根据二维参考图像进行三维重建. 之后Li 等将马尔可夫过程块匹配［11，12］应用到字典建立过程中，使重建过程学习到了先验信息，从而提高了匹配的精确度. Xia 等［13］提出了基于邻域块匹配的超维重建算法，优化了字典的匹配方式，提高了重建精度.

本文基于分区域重建和多线程提出了分区域并行重建算法和字典并行搜索，进一步提高了超维重建算法的效率. 首先，将待重建区域分为多个区域，不同区域进行并行重建；其次，在字典搜索匹配的过程中，使用多线程来对字典元素进行并行搜索匹配. 实验结果表明，基于不同的参考图像，本文提出的算法可以在保证重建出较好的三维结构的基础上减少重建时间，提高超维算法的重建效率.

2 算法理论介绍

2. 1 基于邻域块匹配的超维算法介绍

2. 1. 1 字典建立过程

在训练过程中，将对样本空间的一些三维特征进行学习，比如模式信息等，并将这些三维特征信息建立二维块到三维块的映射关系，通过此映射关系建立超维字典. 在真实的岩心样本中依次截取三维块，基于块和字典的概念，将这些三维块保存为字典，用“Value”表示中心块，“Key”表示邻域块，图1 为匹配对模块的结构和其结构所对应的位置.

邻域块分为左邻域块、前邻域块和下邻域块.考虑到岩心样本的尺度和匹配过程的复杂度，若三维块尺寸越小，字典建立过程越难以学习到真实岩心样本的空间结构关系，若三维块尺寸过大，则需要很大的空间来储存字典，并且匹配时间也会很长. 基于此，我们将三维块表示为5×5×5 的正方体. 字典完整的建立过程如图2 所示. 图2中，使用邻域块匹配对模板并采用光栅路径的方式对真实岩心进行扫描，将每次扫描到的匹配对保存到字典集中作为1 个字典元素.

2. 1. 2 三维重建过程

在重建初期，先建立全空的三维结构，将二维参考图像作为待重建三维结构的第1 层，然后利用光栅路径扫描的方式对待重建的三维结构进行扫描，对于每个待重建块，基于其所有的邻域块在字典中进行搜索匹配，将匹配字典元素的中心块的数值赋予给待重建块. 在全空的三维结构中，底部第1 层和边界区域没有完整的邻域块对字典进行匹配，那么对于这些待重建块只利用那些有限的邻域条件进行匹配. 待重建块和字典元素之间的匹配标准主要参考了邻域块匹配机制和孔隙度控制机制，定义如下.

字典匹配时，会对待重建块的邻域块和字典元素计算匹配因数，如式（1）所示. 同时利用高斯权重分布对不同体素进行加权，高斯权重公式如式（2）所示.

其中M 为匹配因数，DicBlock 为字典元素的邻域块，RecBlock 为重建邻域块，i 为字典元素邻域块和重建邻域块的体素，| DicBlock （ i ）- RecBlock （ i ） |为匹配差的绝对值，Z 为每次匹配时体素的总数.ω （ d ） 为高斯权重分布，d 为不同体素相对于中心块的距离，通常将参数σ 设为1. 5，体素距离中心块越近时，权重越大，体素距离中心块越小时，权重越小.

孔隙度是岩心三维结构的特征参数之一，可以在一定程度上反映岩心三维结构的相似性［13］.因此超维重建要使用孔隙度控制机制来对重建过程进行控制.

在开始重建之前，首先使用式（3）计算参考图像的孔隙度，其中Φref 为参考图像的面孔率，Srefpore为参考图像孔隙相的总面积，Sref 为参考图像的面积. 之后利用式（4）对于每次的重建都计算当前已经重建好的结构的孔隙度，其中Φrec 为当前完成重建的三维结构的孔隙度，Vrecpore 为完成重建的三维结构中孔隙相的体积，Vrec 为完成重建的三维结构的体积.

之后我们使用式（5）对当前的孔隙度进行分类，并且标记孔隙度的状态，我们提出了阈值参数β，β 的定义如式（6）所示. 在重建的过程中，β 的值也随之发生改变，这个改变会使孔隙度的控制机制随重建的进行而逐渐严格. 根据重建经验，将λ的值设为0. 05，ρ 为已完成重建的体积比.

在对每个待重建块进行重建时，我们从字典中对当前的待重建块选择3 个匹配因数最高的候选块，然后从这3 个候选块中计算出孔隙度最大的候选块ElementMax 和孔隙度最小的候选块ElementMin，若当前孔隙度的状态为Low 时，则选择候选块ElementMax，若当前孔隙度的状态为High 时，则选择候选块ElementMin. 若孔隙度的状态处于Normal，则随机从3 个候选块中选择一个.

3 引入多线程的超维重建算法

3. 1 算法设计以及加速方法

基于邻域块匹配的超维算法的重建方向为固定的光栅路径方向，因此在重建过程中，必须待前1 个重建块重建完成后，才可以进行下1 个重建块的重建. 这导致每个时间点都只对1 个待重建块进行重建，从而使整个重建过程效率低. 因此可以将分区域重建引入到超维算法中，在重建过程中，将整个岩心的待重建区域分成多个区域，每个区域使用不同线程进行并行重建，从而实现对多个待重建块同时进行重建. 同时超维算法所建立的字典集含有数万个字典元素，在对三维结构进行重建时，待重建块的邻域块需要和字典的每个元素进行匹配，这导致每个待重建块都要进行上万次的匹配，这也是超维算法十分耗时的原因. 因此可以使用多线程来对字典进行并行搜索匹配，从而实现对该过程的加速. 基于上述的2 个问题，本文设计了分区域并行重建算法和字典并行搜索，并将2 种方法进行结合，具体实现方法如下.

3. 2 分区域并行重建算法

基于本文3. 1 节所述，若要同时对M 个块进行重建，则需要将整个待重建区域分成M 个区域，同时为了保证分区域后超维算法重建出的三维结构和文献［13］中的超维算法重建出的三维结构相似以及考虑到计算机本身的线程数量，本文中取M为2，即将整个待重建区域分为2 个区域，但如果区域划分过多，最终重建出的三维结构会出现一定的不连续性，为了保证分区域后超维算法重建出的三维结构和文献［13］中的超维算法重建出的三维结构相似，本文中取M 为2，即将整个待重建区域分为2 个区域，具体方法流程如图3 所示. 在图3中，将整个待重建区域分为2 个区域进行重建，对于每个扫描到的待重建块，其数值为最终通过匹配的字典元素中心块的值. 当待重建结构中的所有待重建块都完成重建时，整个重建过程就结束.

重建过程中具体的线程选择如式（7）. Thread为所选线程，sl = 5 为重建步长，YRefImage 为参考图像的长，XRefImage 为参考图像的宽. 重建刚开始时，先对0～4 层进重建，. 此时重建路径基于输入的参考图像，从左上角x = 0，y = 0 开启线程Thread1，开始对第1 行进行重建，每次重建完1 个三维块时，x + sl.

当第1 行的待重建块重建完成一半时，既当完成对x = （XRefImage /2） + （sl - （XRefImage /2）% sl ）- sl处重建时，开启线程Thread2，线程Thread2 从x =（XRefImage / 2） + （sl - （XRefImage / 2） % sl ） 处进行重建，并重建至所需宽度.

当开启线程Thread2 时，当前层的右半部分由线程Thread2 负责，不再由线程Thread1 负责. 线程Thread1 负责整个重建区域的左半部分，并按照光栅路径进行重建. 线程Thread2 负责整个重建区域的右半部分，并按照光栅路径进行重建，直至完成对整个三维结构的重建.

由于整个重建过程基于邻域块匹配机制，因此对于每行的重建，只有当线程Thread1 重建完左半部分时，线程Thread2 才可以开启. 而对于整个三维结构的重建，每层的重建要基于之前层进行，即当第l 层至l + 4 层全部重建完成时，才可以对第l + sl 层至l + 4 + sl 层进行重建.

由于整个重建过程基于孔隙度控制机制，为了保证分区域重建的三维结构整体的孔隙度和基于邻域块匹配的超维算法接近，在进行分区域重建时，仍计算整体的孔隙度来对重建进行约束，仍然使用式（4）中的Φrec 代表当前完成重建的三维结构的孔隙度，Vrecpore 此时为图3 中所有已重建块和邻域块中孔隙相的体积，Vrec 此时代表图3 中所有已重建块和邻域块的体积.

3. 3 字典并行搜索

在本文中，对字典进行多线程并行搜索匹配的策略使用了OpenMP［14，15］. 在字典搜索匹配过程中，将字典集分成K 个区域，每个区域都开启1 个独立的线程来进行搜索匹配，并且每个线程都选择1 个匹配因数最高的字典做为候选块，之后利用孔隙度控制机制在K 个候选块中选出最佳候选块.

考虑到2. 1 节中整个字典搜索匹配过程中所选择的候选块的个数以及计算机本身的线程数量，同时OpenMP 中线程的开启和关闭也需要一定的时间，因此本文中将取K 为3. 即开启3 个线程对字典进行并行搜索. 假设字典集的元素个数为N，此时我们令Q = N/K，R = N%K. 字典并行搜索的示意图如图4 所示.

如图4 所示，若N 个字典可以平均分为3 个区域，即R = 0，则每个线程对Q 个字典元素进行搜索匹配，若N 个字典不可以均分，则将最后剩余的R 个字典使用最后1 个线程进行搜索匹配.

4 实验及结果分析

利用计算机的断层扫描（Computed Tomography）可以获得岩心真实的三维微观结构，在真实的岩心三维结构序列图中抽取1 张作为二维参考图像，并通过本文提出的方法进行重建. 本章实验采用 128 GB RAM、Inte（l R） Core（TM） I912900K3. 19 GHz 的64 位处理器的计算机. 将重建结果和文献［13］提出的基于邻域块匹配的超维重建算法、模拟退火算法［1］、MPS 算法［2］重建出的结果进行对比，来验证本文提出的算法的有效性.

4. 1 算法有效性分析实验

分别对高孔隙度、中孔隙度、低孔隙度3 种不同孔隙度的岩心进行重建实验，实验中的每组字典分别是该组岩心的CT 序列图. 我们从3 种CT序列图中抽取的二维参考图像如图5 所示. 它们的面孔率分别为18. 75%、11. 27%、6. 77%，点尺寸长度为15. 0 μm，图片尺寸大小为128×128. 图片中的白色区域为孔隙相，黑色区域为岩石相.

在实验中，将3. 2 节和3. 3 节中的2 种方法进行结合来对3 张参考图像进行10 次重复实验，同时分别使用文献［13］的超维算法、SA 算法、MPS算法进行5 次重复实验来和本文方法进行对比.

4. 1. 1 三维结构对比

从不同算法的重复实验中随机抽取1 组，生成三维结构进行对比，如图6所示. 在图6 中，（a）三维结构A、（f）三维结构B、（k）三维结构C，分别代表了a、b、c 3 张参考图像所对应的真实岩心三维结构，（b）～（e）分别代表了参考图像a 基于不同算法所重建出的三维结构，（g）～（j）分别代表了参考图像b 基于不同算法所重建出的三维结构，（l）～（o）分别代表了参考图像c基于不同算法所重建出的三维结构. 从重建出的三维结构看，本文所改进的超维算法和文献［13］所提出的超维算法在视觉上基本一致，并且和真实岩心的三维结构相似，但SA 算法和MPS 算法所重建的结果和真实的岩心三维结构有一定的差异，从整体结果来看，本文改进的超维算法和文献［13］所提出的超维算法在重建出的三维结构上基本保持一致，并且和真实岩心的三维结构基本保持一致.

4. 1. 2 统计特征参数对比

对于数字岩心，统计特征参数可以反映重建结构的微观信息，常用的统计特征参数有两点相关函数［16］、线性路径函数［17］、两点簇函数等［18］. 而对于三维岩心孔隙相的空间分布，局部孔隙度［19］这一参数可以很好地对其进行反映. 在本节中使用两点相关函数、线性路径函数、两点簇函数、局部孔隙度来对不同算法重建出的三维结构和真实岩心进行分析对比. 对比结果如图7～图9 所示. 图7～图9 分别表示了3 张参考图像基于不同算法重建出的三维结构和真实岩心统计特征参数的对比.

从统计特征参数来看，本文所改进的算法的统计特征参数和文献［13］中的超维算法的统计特征参数十分接近，并且和真实岩心的统计特征参数较吻合，SA 算法和MPS 算法的统计特征参数和真实岩心的统计特征参数的差异要大于本文所改进的算法的统计特征参数. 从曲线走势和达到平稳时的参数来看，本文所改进的算法的重建结果也十分接近于真实岩心.

从局部孔隙度来看，本文所改进的算法和文献［13］中的超维算法相对于真实岩心，局部孔隙度峰值具有一定的差异，但是局部孔隙度分布的展宽基本保持一致. 这表明2 种算法重建出的三维结构和真实岩心孔隙分布情况基本保持一致，但均质性方面有一定的差异. 从整体结果来看，本文所改进的算法重建出的三维结构参数和文献［13］中的超维算法基本保持一致，并且接近于真实的岩心.

4. 1. 3 孔隙度对比

对重建的三维结构进行孔隙度计算，其中对本算法重建10 次的结果计算平均值，对文献［13］中的超维算法、SA 算法、MPS 算法重建5 次的结果计算平均值，计算结果如表1 所示. 从表中可以看出，本文改进的算法和文献［13］中的超维算法的计算结果基本一致.

4. 1. 4 时间对比

为说明引入多线程后超维算法的效率，我们对3 张参考图像的重建时间进行平均对比，其中本文算法进行10 次重建求平均，文献［13］中的超维算法、SA 算法、MPS 算法进行5 次重建求平均，时间对比如表2 所示. 从实验结果可以看出，本文改进的算法相比于文献［13］中的超维算法在速度上有明显提升，提升效率约为75％.同时，本文改进的算法相比于MPS 算法、SA 算法也有时间上的优势. 在保证重建结果的同时，本文改进的算法在重建效率上有明显提升.

表3 为分别单独使用分区域并行重建算法和字典并行搜索对参考图像a 进行5 次重建的平均时间对比. 从表 3 可以看出，分区域并行重建算法相比于文献［13］中的超维算法在时间上的提升约为49％，而字典并行搜索的提升约为60％.

4. 2 泛化性实验

为了证明本文提出的方法具有泛化性，选用不同于3 组参考图像孔隙度的岩心CT 序列图生成字典，再对3 张参考图像分别进行5 次重建. 本组实验采用孔隙度为20. 47% 的岩心建立字典. 图5中的3 张参考图像基于本文算法、文献［13］中的超维算法、SA 算法、MPS 算法所重建出的三维结构如图10 所示，（c）（h）（m）分别代表了3 张参考图像按照新字典所建出的三维结构. 从结果来看，使用不同于参考图像孔隙度的岩心生成的字典来进行重建出的三维结构和文献［13］中的超维算法在视觉上基本一致，并且和真实岩心的三维结构相似.

图11～图13 为重建结构的统计特征参数和局部孔隙度对比结果，从对比结果可以看出本文所改进的算法在使用新的字典所生成的三维结构和真实岩心三维结构相似，参数曲线都较为吻合. 从视觉以及参数分析来看，本文所改进的超维算法具有泛化性.

5 结论

本文在基于邻域块匹配的超维算法的基础上提出了分区域并行重建算法和字典并行搜索. 在重建阶段，将待重建区域分成为个区域，每个区域使用不同线程进行并行重建，实现了多个待重建块同时重建. 在字典匹配阶段，使用OpenMP 开启多个线程来对字典元素进行并行搜索匹配. 在后续的实验中，使用了不同字典来对参考图像进行重建. 从重建结果来看，本文所改进的超维算法可以在保证重建效果的前提下，有效地减少重建所需时间，并且本文所改进的超维算法具有泛化性.

参考文献：

［1］ Yeong C L Y，Torquato S． Reconstructing randommedia［ J］. Phys Rev E， 1998， 57： 495

［2］ Guardiano F B， Srivastava R M. Multivariate geostatistics：Beyond bivariate moments［M］//GeostatisticsTróia’92.［S. l.］： Springer Dordrecht， 1993，1： 133.

［3］ Song S. An improved simulated annealing algorithmfor reconstructing 3D large-scale porous media［J］. JPetrol Sci Eng， 2019， 182： 106343.

［4］ Xu S H， Teng Q Z， Feng J X， et al. Adaptive directsampling core 3D reconstruction algorithm ［J］. J SichuanUniv（Nat Sci Ed）， 2019， 56： 260.［许诗涵，滕奇志， 冯俊羲， 等. 自适应直接取样岩心三维重建算法［J］. 四川大学学报（ 自然科学版）， 2019，56： 260.］

［5］ Bai H， Mariethoz G. A fast edge-based two-stage directsampling method［J］. Computers amp; Geosciences，2021， 150： 104742.

［6］ Mosser L， Dubrule O， Blunt M J. Reconstruction ofthree-dimensional porous media using generative adversarialneural networks ［J］. Phys Rev E， 2017，96： 043309.

［7］ Goodfellow I， Pouget-Abadie J， Mirza M， et al.Generative adversarial networks ［J］. CommunACM， 2020， 63： 139.

［8］ Feng J， Teng Q， Li B， et al. An end-to-end threedimensionalreconstruction framework of porous mediafrom a single two-dimensional image based ondeep learning［J］. Comput Method Appl M， 2020，368： 113043.

［9］ Jia P， Huang Y， Cai B， et al. Solar image restorationwith the CycleGAN based on multi-fractal propertiesof texture features ［J］. Astrophys J Lett，2019， 881： L30.

［10］ He X H， Li Y， Teng Q Z， et al. Learning-basedsuper-dimension （SD） reconstruction of porous mediafrom a single two-dimensional image［C］//2016IEEE International Conference on Signal Processing，Communications and Computing （ICSPCC）.IEEE， 2016： 1.

［11］ Li Y， He X， Teng Q， et al. Markov prior-basedblock-matching algorithm for superdimension reconstructionof porous media［J］. Phys Rev E， 2018，97： 043306.

［12］ Li Y， Teng Q， He X， et al. Super-dimension-basedthree-dimensional nonstationary porous medium reconstructionfrom single two-dimensional image ［J］.J Petrol Sci Eng， 2019， 174： 968.

［13］ Xia Z， Teng Q， Wu X， et al. Three-dimensional reconstructionof porous media using super-dimensionbasedadjacent block-matching algorithm［J］. PhysRev E， 2021， 104： 045308.

［14］ Ayguade E， Chapman B， Mattson G T. How goodis OpenMP ［J］. Sci Programming-NETH， 2003，11： 81.

［15］ Hoffmann R B， Löff J， Griebler D， et al. OpenMPas runtime for providing high-level stream parallelismon multi-cores［ J］. J Super Comput， 2022， 78： 1.

［16］ Kerscher M， Szapudi I， Szalay A S. A comparisonof estimators for the two-point correlation function［J］.Astrophysi J， 2000， 535： L13.

［17］ Singh H， Gokhale A M， Lieberman S I， et al. Imagebased computations of lineal path probability distributionsfor microstructure representation［J］. MatSci Eng A-Struct， 2008， 474： 104.

［18］ Torquato S， Beasley J D， Chiew Y C. Two-pointcluster function for continuum percolation ［J］. JChem Phys， 1988， 88： 6540.

［19］ Hilfer R. Local-porosity theory for flow in porous media［ J］. Phys Rev B， 1992， 45： 7115.

（责任编辑： 白林含）

基金项目： 国家自然科学基金（62071315）