Advanced Search

Indexed by SCI、CA、РЖ、PA、CSA、ZR、etc .

Volume 32 Issue 2
Apr.  2021
Turn off MathJax
Article Contents

Weisheng Hou, Hengguang Liu, Tiancheng Zheng, Wenjie Shen, Fan Xiao. Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s). Journal of Earth Science, 2021, 32(2): 455-467. doi: 10.1007/s12583-021-1443-x
Citation: Weisheng Hou, Hengguang Liu, Tiancheng Zheng, Wenjie Shen, Fan Xiao. Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s). Journal of Earth Science, 2021, 32(2): 455-467. doi: 10.1007/s12583-021-1443-x

Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s)

doi: 10.1007/s12583-021-1443-x
More Information
  • Multiple-point statistics (MPS) is a useful approach to reconstruct three-dimensional models in the macroscopic or microscopic field. Extracting spatial features for three-dimensional reconstruction from two-dimensional training images (TIs), and characterizing non-stationary features with directional ductility are two key issues in MPS simulation. This study presents a step-wise MPS-based three-dimensional structures reconstruction algorithm with the sequential process and hierarchical strategy based on two-dimensional images. An extension method is proposed to construct three-dimensional TIs. With a sequential simulation process, an initial guess at the coarsest scale is simulated, in which hierarchical strategy is used according to the characteristics of TIs. To obtain a more refined realization, an expectation-maximization like iterative process with global optimization is implemented. A concrete example of chondrite micro-structure simulation, in which one scanning electron microscopy (SEM) image of the Heyetang m eteorite is used as TI, shows that the presented algorithm can simulate complex non-stationary structures.
  • 加载中
  • Barnes, C., Shechtman, E., Finkelstein, A., et al., 2009. Patchmatch: A Randomized Correspondence Algorithm for Structural Image Editing. ACM Transactions on Graphics, 28: 1-11. https://doi.org/10.1145/1531326.1531330 doi:  10.1145/1531326.1531330
    Chen, G. X., Zhao, F., Wang, J. G., et al., 2015. Regionalized Multiple-Point Stochastic Geological Modeling: A Case from Braided Delta Sedimentary Reservoirs in Qaidam Basin, NW China. Petroleum Exploration and Development, 42(5): 697-704. https://doi.org/10.1016/S1876-3804(15)30065-3 doi:  10.1016/S1876-3804(15)30065-3
    Chen, Q. Y., Mariethoz, G., Liu, G., et al., 2018. Locality-Based 3-D Multiple-Point Statistics Reconstruction Using 2-D Geological Cross Sections. Hydrology and Earth System Sciences, 22(12): 6547-6566. https://doi.org/10.5194/hess-22-6547-2018 doi:  10.5194/hess-22-6547-2018
    Comunian, A., Giudici, M., Landoni, L., et al., 2018. Improving Bowen-Ratio Estimates of Evaporation Using a Rejection Criterion and Multiple-Point Statistics. Journal of Hydrology, 563: 43-50. https://doi.org/10.1016/j.jhydrol.2018.05.050 doi:  10.1016/j.jhydrol.2018.05.050
    Comunian, A., Renard, P., Straubhaar, J., 2012. 3D Multiple-Point Statistics Simulation Using 2D Training Images. Computers & Geosciences, 40: 49-65. https://doi.org/10.1016/j.cageo.2011.07.009 doi:  10.1016/j.cageo.2011.07.009
    Dimitrakopoulos, R., Mustapha, H., Gloaguen, E., 2009. High-Order Statistics of Spatial Random Fields: Exploring Spatial Cumulants for Modeling Complex Non-Gaussian and Non-Linear Phenomena. Mathematical Geosciences, 42(1): 65-99. https://doi.org/10.1007/s11004-009-9258-9 doi:  10.1007/s11004-009-9258-9
    Feng, W. J., Wu, S. H., Yin, Y. S., et al., 2017. A Training Image Evaluation and Selection Method Based on Minimum Data Event Distance for Multiple-Point Geostatistics. Computers & Geosciences, 104: 35-53. https://doi.org/10.1016/j.cageo.2017.04.004 doi:  10.1016/j.cageo.2017.04.004
    Gueting, N., Caers, J., Comunian, A., et al., 2017. Reconstruction of Three-Dimensional Aquifer Heterogeneity from Two-Dimensional Geophysical Data. Mathematical Geosciences, 50(1): 53-75. https://doi.org/10.1007/s11004-017-9694-x doi:  10.1007/s11004-017-9694-x
    Houlding, S. W., 1994. 3D Geoscience Modeling Computer Techniques for Geological Characterization. Springer-Verlag, Berlin Heidelberg
    Hu, L. Y., Liu, Y., Scheepens, C., et al., 2013. Multiple-Point Simulation with an Existing Reservoir Model as Training Image. Mathematical Geosciences, 46(2): 227-240. https://doi.org/10.1007/s11004-013-9488-8 doi:  10.1007/s11004-013-9488-8
    Ji, L. L., Lin, M., Jiang, W. B., et al., 2017. An Improved Method for Reconstructing the Digital Core Model of Heterogeneous Porous Media. Transport in Porous Media, 121(2): 389-406. https://doi.org/10.1007/s11242-017-0970-5 doi:  10.1007/s11242-017-0970-5
    Kaufmann, O., Martin, T., 2008. 3D Geological Modelling from Boreholes, Cross-Sections and Geological Maps, Application over Former Natural Gas Storages in Coal Mines. Computers & Geosciences, 34(3): 278-290. https://doi.org/10.1016/j.cageo.2007.09.005 doi:  10.1016/j.cageo.2007.09.005
    Lessenger, M., Gladczenko, T., Hardt, J., et al., 2019. Facies-Calibrated Petrophysical and Geocellular Property Modeling Using Data Analytics and Multi-Point Statistics in the Delaware Basin. SPE/AAPG/SEG Unconventional Resources Technology Conference, July 2019, Denver. https://doi.org/10.15530/urtec-2019-421
    Li, J. Q., Zhang, P. F., Lu, S. F., et al., 2019. Scale-Dependent Nature of Porosity and Pore Size Distribution in Lacustrine Shales: An Investigation by BIB-SEM and X-Ray CT Methods. Journal of Earth Science, 30(4): 823-833. https://doi.org/10.1007/s12583-018-0835-z doi:  10.1007/s12583-018-0835-z
    Li, Y., Teng, Q. Z., He, X. H., et al., 2019. Super-Dimension-Based Three-Dimensional Nonstationary Porous Medium Reconstruction from Single Two-Dimensional Image. Journal of Petroleum Science and Engineering, 174: 968-983. https://doi.org/10.1016/j.petrol.2018.12.004 doi:  10.1016/j.petrol.2018.12.004
    Ma, B. J., Wu, S. G., Mi, L. J., et al., 2018. Mixed Carbonate-Siliciclastic Deposits in a Channel Complex in the Northern South China Sea. Journal of Earth Science, 29(3): 707-720. https://doi.org/10.1007/s12583-018-0830-4 doi:  10.1007/s12583-018-0830-4
    Mariethoz, G., Caers, J., 2014. Multiple-Point Geostatistics. John Wiley & Sons, Chichester. https://doi.org/10.1002/9781118662953 doi:  10.1002/9781118662953
    Mariethoz, G., Renard, P., Straubhaar, J., 2010. The Direct Sampling Method to Perform Multiple-Point Geostatistical Simulations. Water Resources Research, 46(11): W11536. https://doi.org/10.1029/2008wr007621 doi:  10.1029/2008wr007621
    Mohammadmoradi, P., Institute of Petroleum Engineering University of Tehran Iran, Rasaeii, M., et al., 2014. Reconstruction of Non-Stationary Complex Spatial Structures by a Novel Filter-Based Multi Scale Mps Algorithm. Open Transactions on Geosciences, 2014(2): 60-71. https://doi.org/10.15764/geos.2014.02007 doi:  10.15764/geos.2014.02007
    Pyrcz, M. J., Deutsch, C. V., 2014. Geostatistical Reservoir Modeling. Oxford University Press, Oxford
    Shen, W., Hu, S., Lin, Y., et al., 2013. Chemical and Petrologic Study of the Heyetang Meteorite. Chinese Journal of Polar Research, 25(4): 386-393. https://doi.org/10.3724/sp.J.1084.2013.00386 (in Chinese wtih English Abstract) doi:  10.3724/sp.J.1084.2013.00386(inChinesewtihEnglishAbstract)
    Song, W. H., Yao, J., Ma, J. S., et al., 2018. Pore-Scale Numerical Investigation into the Impacts of the Spatial and Pore-Size Distributions of Organic Matter on Shale Gas Flow and Their Implications on Multiscale Characterisation. Fuel, 216:707-721. https://doi.org/10.1016/j.fuel.2017.11.114 doi:  10.1016/j.fuel.2017.11.114
    Strebelle, S., 2002. Conditional Simulation of Complex Geological Structures Using Multiple-Point Statistics. Mathematical Geology, 34(1): 1-21. https://doi.org/10.1023/a:1014009426274 doi:  10.1023/a:1014009426274
    Tahmasebi, P., Hezarkhani, A., Sahimi, M., 2012. Multiple-Point Geostatistical Modeling Based on the Cross-Correlation Functions. Computational Geosciences, 16(3): 779-797. https://doi.org/10.1007/s10596-012-9287-1 doi:  10.1007/s10596-012-9287-1
    Tahmasebi, P., Sahimi, M., 2012. Reconstruction of Three-Dimensional Porous Media Using a Single Thin Section. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 85(6): 066709. https://doi.org/10.1103/physreve.85.066709 doi:  10.1103/physreve.85.066709
    Tahmasebi, P., Sahimi, M., 2015. Reconstruction of Nonstationary Disordered Materials and Media: Watershed Transform and Cross-Correlation Function. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 91(3): 032401. https://doi.org/10.1103/physreve.91.032401 doi:  10.1103/physreve.91.032401
    Tang, Y. W., Atkinson, P. M., Zhang, J. X., 2015. Downscaling Remotely Sensed Imagery Using Area-to-Point Cokriging and Multiple-Point Geostatistical Simulation. ISPRS Journal of Photogrammetry and Remote Sensing, 101:174-185. https://doi.org/10.1016/j.isprsjprs.2014.12.016 doi:  10.1016/j.isprsjprs.2014.12.016
    Western, A. W., Blöschl, G., Grayson, R. B., 2001. Toward Capturing Hydrologically Significant Connectivity in Spatial Patterns. Water Resources Research, 37(1): 83-97. https://doi.org/10.1029/2000wr900241 doi:  10.1029/2000wr900241
    Yang, L., Hou, W. S., Cui, C. J., et al., 2016. GOSIM: a Multi-Scale Iterative Multiple-Point Statistics Algorithm with Global Optimization. Computers & Geosciences, 89:57-70. https://doi.org/10.1016/j.cageo.2015.12.020 doi:  10.1016/j.cageo.2015.12.020
    Zhang, T. F., Gelman, A., Laronga, R., 2017. Structure-and Texture-Based Fullbore Image Reconstruction. Mathematical Geosciences, 49(2): 195-215. https://doi.org/10.1007/s11004-016-9649-7 doi:  10.1007/s11004-016-9649-7
    Zhang, W. B., Duan, T. Z., Liu, Z. Q., et al., 2016. Application of Multi-Point Geostatistics in Deep-Water Turbidity Channel Simulation: a Case Study of Plutonio Oilfield in Angola. Petroleum Exploration and Development, 43(3): 443-450. https://doi.org/10.1016/S1876-3804(16)30051-9 doi:  10.1016/S1876-3804(16)30051-9
    Zhou, Z. C., Mei, L. F., Shi, H. S., et al., 2019. Evolution of Low-Angle Normal Faults in the Enping Sag, the Northern South China Sea: Lateral Growth and Vertical Rotation. Journal of Earth Science, 30(6): 1326-1340. https://doi.org/10.1007/s12583-019-0899-4 doi:  10.1007/s12583-019-0899-4
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Figures(15)  / Tables(3)

Article Metrics

Article views(4) PDF downloads(2) Cited by()

Related
Proportional views

Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s)

doi: 10.1007/s12583-021-1443-x

Abstract: Multiple-point statistics (MPS) is a useful approach to reconstruct three-dimensional models in the macroscopic or microscopic field. Extracting spatial features for three-dimensional reconstruction from two-dimensional training images (TIs), and characterizing non-stationary features with directional ductility are two key issues in MPS simulation. This study presents a step-wise MPS-based three-dimensional structures reconstruction algorithm with the sequential process and hierarchical strategy based on two-dimensional images. An extension method is proposed to construct three-dimensional TIs. With a sequential simulation process, an initial guess at the coarsest scale is simulated, in which hierarchical strategy is used according to the characteristics of TIs. To obtain a more refined realization, an expectation-maximization like iterative process with global optimization is implemented. A concrete example of chondrite micro-structure simulation, in which one scanning electron microscopy (SEM) image of the Heyetang m eteorite is used as TI, shows that the presented algorithm can simulate complex non-stationary structures.

Weisheng Hou, Hengguang Liu, Tiancheng Zheng, Wenjie Shen, Fan Xiao. Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s). Journal of Earth Science, 2021, 32(2): 455-467. doi: 10.1007/s12583-021-1443-x
Citation: Weisheng Hou, Hengguang Liu, Tiancheng Zheng, Wenjie Shen, Fan Xiao. Hierarchical MPS-Based Three-Dimensional Geological Structure Reconstruction with Two-Dimensional Image(s). Journal of Earth Science, 2021, 32(2): 455-467. doi: 10.1007/s12583-021-1443-x
  • Three-dimensional geological model is a kind of mathematical model with geological significance, which can provide an intuitive perspective for geological phenomena. At present, three-dimensional geological models can provide multi-perspective understanding, data foundation, and decision support for the geological objects studied in the macroscopic (e.g., Zhou et al., 2019; Ma et al., 2018; Kaufmann and Martin, 2008; Houlding, 1994) and microscopic fields (e.g., Li J Q et al., 2019; Li Y et al., 2019; Song et al., 2018). The key issue of reconstructing three-dimensional structures is to find and clarify the spatial relationship between known data, and infer the attribute values on unsampled locations.

    Combined with the theory of regionalization variable and spatial correlation analysis, geostatistical methods can be applied to estimate the local spatial correlation of the data, and are widely used to construct geological models. These kinds of methods can be roughly divided into two-point statistics (TPS) methods, object-based methods and multiple-point statistics (MPS) methods (Mariethoz and Caers, 2014). Object-based methods use geological knowledge to reconstruct geological bodies with complex geometry. However, the acquisition and application of prior knowledge are limited and subjective. With the assumption of independence of the targets' bodies, simulating the spatial relationship between the geological bodies and matching the hard data become a difficult task in the modeling process (Pyrcz and Deutsch, 2014). TPS methods perform the local optimal unbiased estimation of the unsimulated area, in which the spatial structure relationship of the known data is obtained with the variogram (e.g., Mariethoz et al., 2010). Compared to object-based methods, TPS methods are more efficient to generate models that are faithful to hard data. However, the spatial relationship function extracted by the variogram only represents the pattern between two known points in space. The lack of the measure of the overall spatial correlation of the known data makes it difficult for TPS methods to reconstruct complex geological structures. MPS methods focus on the expression and random recombination of spatially related features among multiple points in training images (TIs), and are more suitable for three-dimensional stochastic reconstruction of complex geological structures (Zhang et al., 2016; Mariethoz et al., 2010). The key issue of MPS methods is to capture spatial structure features from TI(s) and reconstruct these structures using stochastic processes.

    MPS simulation has been widely used in both macroscopic and microscopic fields, such as oil and gas reservoir exploration (Hu et al., 2013), hydrogeology (Comunian et al., 2018), remote sensing mapping (Zhang et al., 2017; Tang et al., 2015), petrophysics (Lessenger et al., 2019), porous media simulation (Chen et al., 2018; Tahmasebi and Sahimi, 2015). Obtaining and selecting TI(s) is one of the key factors for constructing reasonable realization in MPS. However, obtaining three-dimensional TI(s) is time consuming and costly. An ideal initial three-dimensional training model is not easy to build, but sketching two-dimensional geological data can be obtained much easier. Therefore, using two-dimensional image(s) as TI(s) in MPS simulation is an alternative path. Directly extracting the conditional probability distribution from single or multiple two-dimensional TIs is one solution to reconstruct three-dimensional models. This has been implemented in algorithms such as DS (Mariethoz et al., 2010), SNESIM (Strebelle, 2002). Another solution is to use MPS algorithm to generate a series of two-dimensional slices, which are stacked to form a three-dimensional model, such as s2Dcd (Comunian et al., 2012), CCSIM algorithm (Tahmasebi and Sahimi, 2015, 2012; Tahmasebi et al., 2012), CCSIM-TSS method (Ji et al., 2017). As slice number increases in the simulation grid (SG), the accumulation of errors will cause artifacts and structural discontinuities in final results. The hybrid method (Gueting et al., 2018) generates multiple two-dimensional slices as conditional data and TIs, and the remaining areas are simulated by DS algorithm (Mariethoz et al., 2010). But choosing a suitable switch point between the two methods has become a new problem.

    The other challenge is extracting and characterizing non-stationary features with directional ductility in the TI(s). Most MPS methods treat non-stationary features as one of the attributes of TI(s). Some methods use conditional probability (Strebelle, 2002), high-order statistics (Mohammadmoradi et al., 2014; Dimitrakopoulos et al., 2009), data events (Mariethoz et al., 2010) to extract features from TI(s), and ignore the conflicts between the structural properties of geological bodies and the stochastic properties of the simulation process. And some undefined hidden features are difficult to obtain, such as the morphology of strata, which are difficult to characterize with probability. Some methods divide the non-stationary TI(s) into several sub-regions with stationary features, and extract statistical information from the sub-regions to obtain the non-stationary model (Chen et al., 2015). In this scenario, simulation that is implemented in each sub-region separately will lead to discontinuous artifacts along the contacts between the sub-regions. Therefore, the method of partition modeling should be considered more carefully.

    Aiming to solve the two challenges mentioned above, based on GOSIM algorithm (Yang et al., 2016), we propose an MPS-based three-dimensional simulation method using two-dimensional section(s) as TI(s). The presented method combines a hierarchical sequential simulation process and GOSIM iterative algorithm. Auxiliary data based on the features of the simulation objects is used for the sequential simulation process, which helps to build an initial guess for the realization to overcome the convergence problem in the iterative process. Iterative optimization reduces the accumulation of errors caused by sequential simulations to obtain more realistic simulation results (Yang et al., 2016). The following sections of this paper are organized as follows: Section 1 describes the basic idea and the details of the proposed algorithm, and how we combine the hierarchical strategy based on the object's feature with the GOSIM method. In Section 2, a concrete example of the chondrite granules reconstruction is given in which SEM image of Heyetang meteorite slice is used as TI. Discussion of the testing results is implemented in Section 3. In Section 4, a conclusion is drawn.

  • The presented method, in this study, is partially based on the GOSIM (Yang et al., 2016), of which a multiscale iterative strategy is used with global optimization. One of the advantages of the GOSIM algorithm is that error accumulation in the simulation can be corrected by an iterative procedure. Note that GOSIM is suitable for both two-dimensional and three-dimensional simulations with two-dimensional TI(s) and three-dimensional TI(s) respectively. Here, a multiple-scale iterative method combined with a sequential process is proposed, in which two-dimensional images are used as TIs.

    Five main steps are implemented in the presented algorithm (in Fig. 1). The first step is presetting simulation parameters including the SG sizes at each scale, the iterative number of each scale, the pattern size of sequential simulation, and the pattern size of the iterative process etc. Using two-dimensional image(s), the three-dimensional TI(s) at the coarsest scale is constructed for the next step. In the third step, an auxiliary three-dimensional model and sub-pattern databases are built based on three-dimensional TI(s). In the fourth step, an initial guess of the realization is built at the coarsest scale. The final step is to refine the model with the GOSIM algorithm. Several iterations are performed to simulate the structures. This step will be done with several iterative times to obtain the reconstruction at the finest scale.

    Figure 1.  The schematic chart of the presented algorithm.

  • High-quality simulation results are seriously impacted by reasonable TI(s) in three-dimensional MPS-based modeling (Feng et al., 2017; Mariethoz and Caers, 2014). In this study, with two-dimensional image(s), a three-dimensional TI reconstruction method is proposed. The process of constructing three-dimensional TI(s) is described as follows:

    At first, two-dimensional images used for simulation are imported into the SG as shown in Fig. 2a. Then, a random path is constructed to access the unassigned area around the training image, and the thickness of this area is greater than or equal to the preset pattern size. Here, a moving window with a size of 3×3×3 is used to visit each voxel in a random path. The attribute value of the voxel to be assigned is randomly selected from the currently visited voxel by the moving window. Figures 2c and 2d show the process of assigning values to the two voxels q and p in planar view, where the blue frame is the moving window, the gray blocks are the assigned voxels, and the white blocks are the voxels to be assigned. Finally, the above process is implemented repeatedly until values at all voxels in the random path have been assigned, and then three-dimensional TIs with the thickness of about one pattern size are obtained (as shown in Fig. 2b).

    Figure 2.  The process of constructing three-dimensional TIs from two-dimensional images. The three-dimensional prospective views (a) and (b) respectively show that two-dimensional TIs are imported into a SG and expanded to three-dimensional TIs. The planar view (c) and (d) show the process of assigning values to voxels p and q in a random simulation path, where the blue wire frame represents the 3×3×3 moving window, the gray blocks represent the assigned voxels, and the white blocks represent the unassigned voxels.

    Definitely, even if the images have been expanded to a three-dimensional TI, a single two-dimensional TI contains only the two-dimensional patterns. Especially, when anisotropic features exist in the simulation object, whole structures cannot be reasonably characterized in three-dimensional space. To solve the problem that a single TI cannot effectively provide unevenness and variability of three-dimensional spatial features, MPS methods with two-dimensional TI(s) reconstruct three-dimensional models generally using two or more perpendicular profiles (Gueting et al., 2017; Comunian et al., 2012). Therefore, realizations with a pair of mutually perpendicular TIs are obtained for analysis in this paper.

  • The GOSIM algorithm randomly generates an initial model before multiscale iteration. For binary images without complex geometric features, such as river channels, GOSIM algorithm can obtain reasonable simulation results. However, the initial model generated by the sequential simulation, when reconstructing non-stationary ductility features, will lead to more iterations to obtain a reasonable realization. Therefore, a reasonable initial model is particularly important in the simulation.

    A hierarchical strategy, in this study, is introduced in the process of building the initial model, and auxiliary data is used to guide the pattern selection. According to the characteristics of the simulation object and prior knowledge, including information obtained from TI(s) and expert knowledge, a three-dimensional model is built firstly as auxiliary data that divides the SG into different regions. Then in the sequential simulation process of each sub-region, only the sub-pattern database corresponding to the sub-region is used.

    To build the sub-pattern databases by pattern extraction from TI(s), it is necessary to classify these patterns. For example, classifying different layer properties in a sedimentary cycle into the same category can effectively reduce complexity of modeling. By scanning the three-dimensional TI(s), according to the pattern classification rules that are given artificially, three-dimensional sub-pattern databases of multiple categories are constructed for subsequent sequential simulation process.

    The goal, here, is to reconstruct the three-dimensional geometry of the central chondrite in Fig. 3, in which the TI is a SEM image of the chondrite slice (Fig. 3a). Note that the central olivine granules are similar to the spheroid in three-dimensional space, the maximum radius of the central olivine granules r can be obtained by scanning the two-dimensional TI(s). Therefore, according to the idea above and prior knowledge, to build an auxiliary three-dimensional model, we assumed that the central olivine granules values are distributed in a sphere with a radius r where the center is located on the center of the SG. Thus the region outside the sphere is classified as sub-region a of the auxiliary three-dimensional model. In order to reduce the splitting artifacts due to hierarchical modeling, the sphere area at the center of the auxiliary three-dimensional model is divided into two sub-regions. The region where the distance from the center of the sphere is less than t×r is used as region c, and the rest region of the center sphere is marked as region b. The t is a preset parameter. The red circles in Fig. 3b show the divided range of three sub-regions in a two-dimensional perspective. Then, three sub-pattern databases of A, B, and C of the three sub-regions a, b, and c respectively are obtained from the three-dimensional TI(s). The patterns in the A database do not contain the value of the central olivine granules. The difference between B and C lies in the proportion of central olivine granules in patterns. The patterns of which central olivine granules proportions are less than the threshold y will be assigned to database B, and other patterns will be assigned to database C.

    Figure 3.  (a) The SEM image of chondrite slice; (b) the manual interpretation image of chondrite slice. The different gray value represents a different mineral component. The red circles represent the divisions of the TI by each sub-region in a two-dimensional perspective.

    Note that, the patterns in the sub-pattern database A do not contain the value of the central olivine granules, which limits the distribution of the central olivine granules in the reconstructions. It ensures that other attribute values do not be pasted too much in the hypothesized chondrite simulation area. After testing, when the threshold y is set to 0.8 and t is set to 0.5 in this study, more reasonable results can be obtained.

  • After constructing three-dimensional TIs and the sub-pattern databases for simulation, in this study, the SG is initialized by a sequential simulation process with a hierarchical strategy. In order to speed up the sequential simulation, each sub-pattern database for storing three-dimensional pattern is constructed by scanning the constructed three-dimensional TI(s), in which extracted patterns are clustered. The geometry of the possible overlapping area according to the simulation path (Fig. 4) and the similarity value are mainly factors that are considered in patterns clustering. In this study, Hamming distance is used to calculate the similarity value, which counts the number of voxels that are different in the same overlap zone between two patterns.

    Figure 4.  Basic types of the overlapping zone in the simulation. ps is pattern size, os is the thickness of the overlap region, and gray regions represent the overlap zone. The more complex overlapping zones are combined by the above overlapping zones.

    Based on the hierarchical strategy, the initialization process that is implemented with a sequential simulation includes mainly four steps.

    Firstly, a random path is defined to scan the unassigned voxels in the SG. In this step, it should be ensured that there is at least one overlap zone between the area already simulated in the SG and the area to be simulated around each voxel in the random path.

    Secondly, according to the type of overlap zone (Fig. 4) of the unassigned voxel on the path and the sub-region of the auxiliary three-dimensional model, the similarity calculation is performed with the patterns in the sub-pattern database corresponding to the sub-region. The top f similarity patterns will be selected as candidates, where f is a preset parameter.

    Thirdly, one of the candidate patterns that satisfy the similarity requirement between the current data event and the overlap region, is randomly selected and pasted on the position of the current data event. In each search step, candidate patterns are obtained based on the probability determined with the Inverse Distance Weighted Method.

    Finally, the algorithm repeats the above steps until the whole simulation path is scanned.

  • Although the initial results above are a kind of realization, there are still many discontinuities and artifacts caused by pattern pasting and hierarchical modeling. Hence, a GOSIM-based iterative step is used to refine the realization.

    The initialization is implemented at the coarsest scale. Realizations with finer scale are upsampled from the results at the coarse-scale. This is an iterative process as shown in Fig. 1. In the simulation at each scale, several iteration processes are executed with the GOSIM algorithm. Searching (E-Step) and updating (M-Step) are included in each iteration of each scale.

    In the E-step, each of the voxels ui,j, k∈R(i=1, 2, ..., n1; j= 1, 2, ..., n2; k=1, 2, ..., n3) except the conditional data in the SG (with the size n1×n2×n3) is randomly assigned a candidate pattern randomly extracted from the TI(s), and the distance between the patterns and the data events in the initial model is calculated. Let $ {p_{i', j', k'}} $∈TI$ \left({i' = 1, 2, \ldots, {m_1};j' = 1, 2, \ldots, {m_2};k' = 1, 2, \ldots, {m_3}} \right) $ denote the voxels in the TI(s) (with the size m1×m2×m3), $ C\left({{p_{i', j', k'}}} \right) $ is the pattern in TI(s) that is centered on voxel $ {p_{i', j', k'}} $. Referring to the Patch Match Method (Barnes et al., 2009), the coordinates of $ {p_{i', j', k'}} $ in the TI(s) corresponding to each voxel except the conditional data in SG are temporarily stored in a SG F with the same size of the SG. That is $ F\left({{u_{i, j, k}}} \right) = \left[ {i', j', k'} \right] $. Then, the search is implemented with two steps. The first step is the propagation process, which scans the voxels in the SG F. When accessing voxel F(ui, j, k), first found neighboring voxels in F according to the preset offset, such as F(ui+1, j, k) and F(ui, j-1, k), are corresponding to the candidate pattern C(F(ui+1, j, k)) and C(F(ui+1, j, k)). After that, offset patterns of the patterns stored in the adjacent voxels, and the original candidate pattern C(F(ui, j, k)), are compared with the current data event. The pattern with the largest similarity is selected as the candidate pattern to replace the original candidate pattern. The coordinate of this is assigned to the current voxel F(ui, j, k). The second step of the search is a random process, in which a search window is set in the three-dimensional TI(s) on the current corresponding $ {{p_{i', j', k'}}} $ coordinate, and the random extracted pattern in the window is compared with the current candidate pattern $ C\left({{p_{i', j', k'}}} \right) $. If there is a pattern more similar to the current data event, then the F(ui, j, k) is updated with the coordinate of this pattern. After each window is searched, the search is continued by reducing the window size according to a preset parameter.

    After the search process in the E-step, the values of the SG are updated via the M-step, to find an appropriate pattern based on the patterns found in the E-step. When updating the attribute values in the M-step, the preponderant area method is used instead of the weighted average method. Because each attribute value in the SG is only used as a code for the formation attribute or granule attribute, the weighted average method may account for incorrect properties.

    Candidate patterns from multi-TIs also need to be integrated. Assuming that the number of TIs is W, Let Cw(Fw(ui, j, k)) denotes the candidate pattern from the TIw(w=1, 2, ..., W) corresponding to the voxel ui, j, k in the SG. Cv(Cw(Fw(ui, j, k))) represents the attribute value of the candidate pattern central voxel. Then, the updated value of the current grid R(ui) is the property of the maximum frequency of all candidate patterns from different TIs, and can be denoted as

    where freq(·) is the frequency of the properties of related patterns. In each EM-step, the E-step usually executes several times when searching for the most similar pattern and the M-step is implemented once.

  • In this study, the presented algorithm is tested with reconstructing three-dimensional chondrite granules, in which a SEM image (Fig. 3) of the Heyetang meteorite is used as TI (Shen et al., 2013). The Heyetang meteorite, which fell in October 1988, is the most primitive meteorite among non-Antarctic ordinary chondrites in China. It retains the information that is of great significance to study the origin and evolution of the solar system. The selected image is one of backscatter electron images of Heyetang Meteorite. In this image, the chondrites have clear structures that are spherical or sub-spherical. It is composed of various chondrites, opaque minerals, and matrix filling between these granules.

    The SEM image of the chondrite section has been manually interpreted as a two-dimensional TI. The interpreted image (Fig. 3b) is a grayscale BMP image with a size of 270×400. The actual size of one pixel is 2.5 μm×2.5 μm. In the manually interpreted image, different attribute is represented with different gray values. The primary purpose of this modeling is to simulate the olivine granules in the center of the image. So the attributes of olivine granules in the SEM image are divided into central olivine granules and peripheral olivine granules. The gray value 33 is the olivine granules in the center, and 56, 64, 74 are the olivine granules in the periphery. The gray value 135 is FeO, 150 is FeS, 194 is Fe-Ni, and 115 is the matrix. The algorithm is implemented based on a python environment, using a desktop workstation with a 20-cores CPU and 128 G of memory for calculations. A multi-process parallel acceleration is used in the GOSIM-based iterative step. Two scales are used in the simulation. The finest scale of the SG is 270×400×400. We tested the presented algorithm under different simulation parameters as shown in Table 1. The time for calculating one realization is no more than 10 h. The visualization of the results is based on Paraview Software.

    Pattern size for initialization Pattern size for iteration Using hierarchical strategy Using vertical TIs
    Para 1 7×7×7 5×5×5
    Para 2 7×7×7 5×5×5 -
    Para 3 7×7×7 5×5×5 - -
    Para 4 7×7×7 5×5×5 -
    Para 5 5×5×5 5×5×5
    Para 6 9×9×9 5×5×5
    Para 7 7×7×7 7×7×7
    Para 8 7×7×7 9×9×9
    - means that the parameter is not used in the simulation.

    Table 1.  The simulation parameters in testing the presented algorithm

  • The size of the initial SG is 135×200×200. Assuming that the samples have certain transverse isotropy in space, two TIs will be perpendicular to each other and introduced into the SG, as shown in Fig. 5a. The copy of the original TI is imported as soft data, and provides constraints perpendicular to the original TI. In the sequential simulation, the pattern size is 7×7×7. The initial model generated by sequential simulation based on the hierarchical strategy is shown in Fig. 6a. The pattern size used in the iterative step is 5×5×5, and the number of executions of the propagation and random search processes is three times in each iteration. Compared with the initial model, artifacts are reduced in the realization after the iterative process as shown in Fig. 6b.

    Figure 5.  (a) TIs for simulation; (b) the manual interpretation image of chondrite slice.

    Figure 6.  Simulation results using parameters of Para 1. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

    From the visual perception, the central olivine spheroids in the results are ellipsoidal, which meets the expectations under the assumptions. The olivine and other granules in the periphery are also reasonably distributed (Figs. 6b-6e), but the shape of some granules is different from the granules in the TI.

  • As shown in Fig. 7, the realization reveals reasonable geometry and topology of chondrite granules. When the copy TI as the soft data is removed in the simulation with Para 2, only one TI that extended from the image of Fig. 3b is imported in the SG (Fig. 5b). The size of granules in the initial results in Fig. 7a is larger than that in Fig. 6a. Also, the geometry is stretched along the y-axis. Compared with the simulation results in Figs. 7b-7e, the distribution and morphology of the outer olivine granules and other granules in the results are similar. But, the central olivine spheroids as shown in Fig. 7c are flatter in the direction perpendicular to the hard data, though the ductility of the direction is preserved along the x-axis.

    Figure 7.  Simulation results using parameters of Para 2. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

  • When the simulation is implemented without any auxiliary three-dimensional TIs and hierarchical strategy, as the Para 3 in Table 1, the realization will be fully decided by the pattern size and stochastic reconstructing process. The construction of the initial model did not use the hierarchical strategy and directly used all the patterns extracted from the three-dimensional TI. Granules in all results including initial (Fig. 8a) and final results (Figs. 8b-8e) are distorted along the y-axis direction because of lacking hierarchical strategy. The geometry of granules along the x-axis direction parallel to the TI direction is more reasonable, as shown in Figs. 8a-8e. When two TIs (Fig. 5a) are used as the Para 4 in Table 1, also, the geometry of granules in realizations (Figs. 9a-9e) is stretched. Note that the central olivine granules distribute outside of the central area in the TIs in Fig. 8c and Fig. 9c. It means that the hierarchical strategy plays a major role in simulating the reasonable granule geometry, although more TIs provide more constraints in the simulation.

    Figure 8.  Simulation results using parameters of Para 3. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

    Figure 9.  Simulation results using parameters of Para 4. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

  • The pattern size in the algorithm is related to the size of the SG. When a multiple-scale strategy is used in the implementation, the changing of the SG size will ensure that the multiple-scale structure can be captured (Yang et al, 2016). In the presented algorithm, the pattern size can be different in the initialization and the iterative process. When changing the pattern size for initialization and other parameters are the same as Para 1 (Para 5 in Table 1), the size of granules in the initial realization is much larger than that in Fig. 6a. In the final realization, the central olivine granules (Fig. 10c) become more concentrated. Other granules in Figs. 10d-10e) are larger than those in Figs. 6d-6e and show a stronger continuity. When the pattern size increases to 9×9×9 as the Para 5 and other parameters are the same as the Para 1, granules size of the initial and final realization (in Figs. 11a-11e) become lager, and numbers of granules decrease. When the pattern size for the iterative process increases (e.g., Para 7 and Para 8) and other parameters is the same as Para 1, the volume of the central granules realization becomes more concentrated (Figs. 12c and 13c) and granule size increases as the pattern size increases as shown in Figs. 12d, 12e and Figs. 13d, 13e.

    Figure 10.  Simulation results using parameters of Para 5. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

    Figure 11.  Simulation results using parameters of Para 6. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

    Figure 12.  Simulation results using parameters of Para 7. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

    Figure 13.  Simulation results using parameters of Para 8. (a) The initial results; (b) the final realization; (c) the olivine granules in the center; (d) the olivine granules distributed on the periphery; (e) FeO, FeS, Fe-Ni granules.

  • As a statistical feature of micro-structure, connectivity analysis is of great significance for analyzing the physical properties, geometric characteristics, and distribution of geological structures (Western et al., 2001). In this study, the two-point connectivity probability function is calculated for the quantitative analysis of the geometric characteristics of granules in TI and realizations. The shapes of the two-point connectivity probability functions of realizations and TI appear similar (in Fig. 14). In the realizations without hierarchical strategy (Fig. 8 and Fig. 9), the granules, especially the central olivine, exhibited abnormal ductility (Para 3 and Para 4 in Fig. 14). The volume proportion of the granules in the realizations with different parameters is different as shown in Table 2.

    Figure 14.  The curve of the two-point connectivity probability function for the TI and the results

    Proportion of component (%) Para 1 Para 2 Para 3 Para 4 Para 5 Para 6 Para 7 Para 8
    Central olivine granules 8.62 7.84 13.62 10.57 8.77 8.67 8.55 8.60
    Peripheral olivine granules with high iron 15.07 21.53 16.79 11.39 10.70 15.63 13.46 13.41
    Peripheral olivine granules with medium iron 5.65 4.38 6.16 5.97 6.52 5.45 3.93 3.05
    Peripheral olivine granules with low iron 0.70 0.50 2.02 1.03 0.53 0.55 0.53 0.48
    FeO granules 0.17 0.08 0.81 0.63 0.11 0.32 0.12 0.09
    FeS granules 0.87 0.75 1.42 1.26 0.55 1.05 0.75 0.79
    Fe-Ni granules 0.40 0.24 1.16 0.35 0.50 0.28 0.24 0.21
    Matrix 68.52 64.68 58.02 68.80 72.32 68.05 72.42 73.37

    Table 2.  Granules proportion in realizations with different parameters

    Unlike the macro-structures, the micro-structures have more singularities. As a quantitative index, the fractal dimension can describe the complexity and irregularity of the micro-structure. Here, Minkowski-Bouligand dimension is used to describe the geometric characteristics of the central olivine granules. Table 3 shows the fractal dimension of the central olivine granules in the TI, the results and two-dimensional slices of which the azimuth angle are 45° and 135° in each result. The fractal dimension of central olivine granules in the three-dimensional model with hierarchical strategy is between 2.434-2.507. And the fractal dimension of the central olivine granules in the two-dimensional slices of each result is slightly smaller than the fractal dimension in the TI. Note that the fractal dimension of central olivine granules in the results without hierarchical strategy (Para 3 and Para 4) is slightly smaller than in the results with hierarchical strategy.

    TI Para 1 Para 2 Para 3 Para 4 Para 5 Para 6 Para 7 Para 8
    Realizations 2.495 2.459 2.323 2.383 2.498 2.434 2.485 2.507
    Slice in 45° 1.709 1.650 1.660 1.557 1.620 1.677 1.662 1.675 1.669
    Slice in 135° 1.709 1.679 1.659 1.610 1.651 1.689 1.694 1.682 1.695

    Table 3.  Central olivine granules' fractal dimension in realizations with different parameters

    In essence, the pattern used for the GOSIM-based iterative process only extracts local information. For the ductile structures that are larger than the pattern size, the optimization process does not perform well. It illustrates that the hierarchical strategy in this paper is effective in simulating non-stationary targets.

    Simulation results show that the shape of the object is affected by the pattern size used in sequential simulation and iterative process. When the pattern size for sequential simulation is too large, it is not conducive to reconstruct the spatial structures with a small scope. If the pattern size is too small, the extension characteristics of the spatial structure will be exaggerated. Therefore, the pattern size should be carefully set in a proper value range.

    The proportion of the central olivine granules in TI is 25.1%. When the central granule of the TI is an ideal circle, the volume ratio of the sphere with the same diameter in the two-dimensional image is about 7.75%. Therefore, in the simulation results of the example using hierarchical strategy, the volume proportion of the central olivine granules is about 7.55% to 8.77%. This proportion is basically in line with expectations. According to the proportion of different components in realizations as shown in Table 2, pattern size has a limited impact on the proportions in the final results.

    The volume proportion of the simulation granules except the central olivine granules, such as FeO, FeS, and Fe-Ni granules, shows relatively large fluctuations. In the simulation results, the maximum volume proportion of FeO granules is four times the minimum volume proportion. Note that the distance calculation for pattern matching using Hamming distance has the same weight for all attributes. As a result, the attribute values with small volume proportion in the patterns are likely to be ignored during the calculation, and uncertainty of these attribute values in the reconstruction increases.

    Some of the olivine granules in the periphery of the simulation results with the hierarchical strategy (Fig. 15) are still unnatural in shape and show excessive continuity. Essentially, the three-dimensional TI used in this study is a simple extension of the two-dimensional image, and it does not guarantee that sufficient morphological information is provided in the direction of expansion. To solve this problem, it is necessary to add constraints on the morphology of peripheral granules in the sequential simulation process, or to adjust the two-dimensional TI extension method based on the morphological characteristics of the peripheral olivine granules.

    Figure 15.  The peripheral olivine granules with abnormal geometry in the simulation results with parameters of Para 1.

    Labeling attribute of each grid is one of the problems in the presented method. The presented method, in essence, is implemented on the gray value of the interpreted image. Without interpretation, structure boundaries cannot be simulated and revealed effectively with the proposed algorithm. Because the original image is not strict in discriminating attributes based on gray values, the same category will show different attribute values in different positions. In the original SEM image, the olivine granules in the center of the image have no special gray value. Therefore, it is still costly to classify granules and matrix by gray value.

  • A novel MPS simulation method that combines the sequential simulation and iterative process for constructing three-dimensional structures is presented, in which two-dimensional images are used as TIs. Three-dimensional TI(s) for simulation is constructed by extending two-dimensional images. A hierarchical strategy is introduced in the initialization process according to the non-stationary ductility characteristics of the simulation object. The auxiliary three-dimensional model and sub-pattern databases guide the selection and paste of candidate patterns, which includes expert knowledge and morphological information. An optimization GOSIM-based iterative process is implemented for obtaining a fine realization with the initialization.

    The example of chondrite granules simulation, in which the TI is a two-dimensional SEM image, illustrates that the micro-structures can be effectively simulated with the presented algorithm. The hierarchical strategy divides the SG into several sub-regions to constrain the impact of non-stationary ductility to the realization. Pattern sizes for initialization and iterative process are key parameters to constrain the geometry of granules' structures. The granules' proportions, the fractal dimension and the curve of the two-point connectivity probability function of different realizations ensure the randomness of geometry of the results and illustrate the spatial uncertainties of realizations.

  • The authors greatly appreciate the valuable comments of the anonymous reviewers in refining this manuscript. This research was substantially supported by the National Natural Science Foundation of China (NSFC) Program (Nos. 41972302, 41772345). The final publication is available at Springer via https://doi.org/10.1007/s12583-021-1443-x.

Reference (32)

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return