Modified FILTERSIM Algorithm for Unconditional Simulation of Complex Spatial Geological Structures

Facies and fracture network modeling need robust, realistic and multi scale methods that can extract and reproduce complex relations in geological structures. Multi Point Statistic (MPS) algorithms can be used to model these high order relations from a visually and statistically explicit model, a training image. FILTERSIM as a pattern based MPS method attracts much attention. It decreases the complexity of computation, accelerates search process and increases CPU performance compare to other MPS methods by transferring training image patterns to a lower dimensional space. The results quality is not however as satisfactory. This work presents an improved version of FILTERSIM in which pattern extraction, persisting and pasting steps are modified to enhance visual quality and structures continuity in the realizations. Examples shown in this paper give visual appealing results for the reconstruction of stationary complex structures.


Introduction
Static modeling of complex reservoirs calls for more than two point statistics.Certain features of these objects such as curvilinearity cannot be expressed via two-point relations [1].Multi Point Statistics (MPS) is used to produce realistic realizations of complex structures after the pioneer work of Strebelle who applied image processing concepts for this purpose [2].MPS based on image processing techniques and proposes image construction approach.In this approach, the patterns of the final image to be constructed (facies, fracture network…) are obtained from a training image that depicts relevant geological patterns expected to be found in the subsurface [1].Furthermore, MPS modeling has the potential to improve estimation in addition to simulation and the so called E-type average value at each node is taken as a local estimated value if training image is appropriate [3].
Training Images (TIs) provide spatial relations, existing patterns and degree of connectivity which must be reproduced in the reconstructed images.The goal is not to reproduce the training image, but to simulate a model that shares some of the multivariate characteristics of the true distribution [4].Achieve the most appropriate TI rise several challenges like: statistical richness of the training image, the scale of the training image, the grid definition of the simulated model, and the univariate distribution of the training image.In this area, Boisvert selects the most data representative TI by comparing the distribution of runs and the multiple point density function from TIs [5].Boogaart applied variography as an auxiliary element to choose TIs which represent the true joint distribution of the field under consideration [6].
Pattern to point SNESIM method is the first algorithm for MPS modeling [7,8].SNESIM is a combination of traditional sequential simulation with a search tree concept.Multiple point statistics (observed patterns) of a certain size of training images are stored in a tree structure; Node properties of realizations are then assigned one by one in a search process loop.Search tree is extremely memory demanding.Straubhaar proposed a list approach modification of SNESIM [9].Several case studies of successful application of the method are reported [10][11][12][13][14]. Okabe also modified SNESIM for pore network reconstruction [15].
Pattern to pattern methods such as SIMPAT [1] FIL-TERSIM [16] and DisPAT [17] are completely isolated from two-point statistics and eliminate the probabilistic paradigm in MPS algorithms.These methods inherently take the probabilities of the whole multiple point patterns conditioned to the same multiple point data event from the training image [18].SIMPAT produces visual appealing realizations; however complete search in pattern database is needed.This makes the process of finding the most similar pattern too complex and CPU time consuming.FILTERSIM reduces pattern dimensions by applying spatial filters and transferring patterns to a lower dimensional space.Coded patterns are then clustered and a prototype is chosen for each bin.This accelerates the search process and decreases run time compared to SIMPAT by means of reducing number of analogy loops.Structures continuity in FILTERSIM method however is not preserved as it is in SIMPAT.
Dimensional reduction in patterns is done by applying directional filters [16] or wavelet decomposition on training image [19].Honarkhah proposed a distance based algorithm DisPAT and automatic information theory based algorithms for selection of optimum template size and dimension of patterns target space [17].
In this paper, we proposed an improved version of FILTERSIM method.Modifications are made to the template size and pattern selection algorithm so that features connectivity enhance in realizations.Although computation cost increases marginally compare to the original method, results quality is far better in our Modified FILTERSIM approach.Frequency component of T has the advantage of accelerating selection of appropriate pattern in simulation process.This achieves by elimination of repeated patterns and consequently repeated calculations.Pattern stores larger scale spatial information as the template size increases with the cost of simulation uncertainty and complexity.To find the optimum template size in a single grid simulation, Shannon's entropy versus template size plot is used [20].In order to calculate entropy in different template sizes, the following equation is applied:

Training Image Processing
where K = number of possible outcomes of the random variable, and p i represents the probability mass function.
Two different behaviors are expected with increasing the template size.In the first stage, entropy will sharply increase since the average number of nodes needed to guess the remaining nodes in the patterns is increased.At a later stage, where the template size has passed the optimal range in an elbow that represents the stationary features of the training image, entropy would increase slowly.This is because the information needed for encoding a large pattern with some stationary features is approximately equal the information of the stationary feature of the training image itself [17].Mathematical details are discussed extensively by Honarkhah [17].As shown in results of FILTERSIM automation in , , , , , , Frequency Copyright © 2012 SciRes.GM to six directional filters are measured.These filters are listed in Table 1.Each filter score is divided into five categories in order to cluster the patterns into bins.Many of these created bins are empty and will be eliminated in the search process which is only done on selected prototypes for the filled bins.Each prototype is binary average pattern of patterns in a specific bin.Table 2 shows the TI processing steps of our Modified FILTERSIM approach.
Having processed the TI and filled the pattern database as just described, reconstruction of realizations can performed in unconditional and conditional situations.Unconditional simulation procedure is discussed in detail in the next section.

Unconditional Simulation Algorithm
Simulation begins by defining a random path on realization nodes.Data event in the template is extracted in each node and its similarity to the existing patterns is then measured.

Data event vector
T is a set of previously informed nodes in unconditional simulation while hard data is also added to T in conditional simulation.We used the following single point manhattan distance function: Having the random path defined, a search among the prototypes is done for each visited node to select the most similar prototype which shows minimum distance function , d x y with data event as below:

dev u Prototype h d dev u h Prototype h
The bin of selected prototype contains the most similar patterns to the visited data event.FILTERSIM routine is to select a random pattern among these patterns and set it as realization node.This obviously is not the best procedure and may lose the feature connectivity as is a known drawback of this method.Here we modified the procedure and performed a second searching step to find the most similar pattern to data event in the selected bin (Figure 3).Pattern selection in candidate cluster is based on the following distance function: 0 ( ), ( ), ( ) Only the pattern nodes with known values are considered in the above relation.d x y x y The pattern with minimum distance function is then pasted entirely on realization nodes.If some patterns have the same distance to the data event vector, the one with the most frequency is selected.Random selection of patterns is only allowed when they all have the same frequency.This pattern selection strategy has improved the FILTERSIM performance and connectivity features of reconstructed realizations compare to the original method.This will be fully explored in the results section.Table 3 shows the stepwise procedure of Modified FIL-TERSIM approach described above.
The last but not the least point is the updatable region of pasted pattern on realization nodes which is used for the first time by Arpat in SIMAPT [1].The selected pattern after imposing on the visited node of realization is not entirely fixed and some outer parts of it can be replaced by future pasted patterns as shown in part 4 of Table 3.This provides the opportunity of superior realizations with closer features of connectivity to those of TI. (Figure 4).

2.
Define a random path on re G grid of realization re.

3.
At every node u, extract the data event dev ( ) T u from realization re and find the most similar prototype that minimizes dev ( ), Once the most similar bin is found, search most similar pattern in that bin to maximize similarity s d , if some patterns have the same distance to the data event vector, most frequent of them is selected and if all of them have the same frequency, selection will be done randomly.dev T u +

5.
Move to the next updatable node of the random path and repeat the above steps until all grid nodes are fixed.

Results and Discussion
Modifications are made to the original FILTERSIM algorithm in several ways of optimizing template size, considering pattern frequency component in database and pattern distance ranking in candidate bin instead of random selection.Performance of this modified algorithm is investigated in this section by applying the method on several binary 2D training images with grid    6, the best template size of 13 × 13 is obtained from entropy plot (Figure 8).This template size preserves the TI features more precisely in the generated realization compare to other template sizes of 9 × 9 and 17 × 17 as shown in Figure 6.
Performance of modified FILTERSIM is also investigated by applying the method for a continuous 2D training image with grid size of 128 × 128 (Figure 9).In this case better visual continuity of our algorithm is clear which is obtained the optimum template size of 15.

Summary
We proposed a Modified FILTERSIM algorithm to simualte 2D binary stationary structures.We modified the original algorithm with introducing optimum template size for each training image, a pattern frequency parameter in database, and additional search among stored patterns in the selected bin.Shannon entropy concept was applied in FILTERSIM and used to infer optimum template size that can capture local structures in each template.Frequency    component increase the chance of more repeated patterns to be selected; and prototypes and patterns are chosen based on their minimum manhatan distance function with data events within the realization under construction.These modifications considerable improved the algorithm performance and produced more visual appealing images compare to FILTERSIM.These improvements add marginal computation cost to FILTERSIM and yet still much faster than SIMPAT.
Our algorithm, SIMPAT and original FILTERSIM have been implemented in C#.We have made our code stand alone and provided a user interface for uploading training and hard data to display simulation results.In the future, we plan to extend Modified FILTERSIM algo-rithm to reconstruction of non-stationary images and 3D pore space structures.The later will be based on 2D real thin section images.

Figure 1 .
Figure 1.Schematic of training image preprocessing steps consist of gridding, scanning and pattern extraction (As template size increases, pattern stores larger scale spatial information).

Figure 2 ,
spatial statistics of the results are fixed with the passing of the Elbow in entropy plot.By applying frequency component and optimum template size concepts, training image is scanned with the optimum template after which the responses of patterns

Figure 2 .
Figure 2. Effect of template size on simulation results (where the template size has passed the optimal window, stationary features in realizations are approximately equal the information of the stationary feature of the training image).

Figure 3 .
Figure 3. Schematic of searching between prototypes in Step#1 and between selected cluster patterns in Step#2 (In Modified FILTERSIM , Step#1 choose 2 nd column and Step#2 chooses the red pattern while in the FILTERIM Step#2 is performed in random).

4 .
Paste the entire selected pattern to updatable data event:   * = pat T α h .Set inner nodes as fixed and outer nodes as updateable.

Figure 4 .
Figure 4. Effect of updatable nodes on structures continuity.(a) Training image; (b) Realization with updatable region; (c) Realization without updatable region.

Figure 5 . 8 .
Figure 5. CPU Time Qualitative Comparison of SIMPAT, FILTERSIM and Modified FILTERSIM methods.size of 128 × 128.Figures 6 and 7 compare reconstructed realizations of 2 different TIs at 3 template sizes with unconditional algorithms of FILTERSIM, SIMPAT and Modified FILTERSIM.Better visual continuity of the Modified FILTERSIM results is clear in all these cases.Entropy changes of the two TIs at various template sizes are shown in Figure 8. Elbows are around the optimum template size which is consistent with the superior results of these sizes in Figures 6 to 7. For example, for the first training image in Figure6, the best template size of 13 × 13 is obtained from entropy plot (Figure8).This template size preserves the TI features more precisely in the generated realization compare to other template sizes of 9 × 9 and 17 × 17 as shown in Figure6.Performance of modified FILTERSIM is also investigated by applying the method for a continuous 2D training image with grid size of 128 × 128 (Figure9).In this case better visual continuity of our algorithm is clear which is obtained the optimum template size of 15.

Figure 6 .
Figure 6.Simulation results with three different template sizes by applying SIMPAT, FILTERSIM and Modified FILTESIM methods.

Figure 7 .
Figure 7. Simulation results with three different template sizes by applying SIMPAT, FILTERSIM and Modified FILTESIM methods (2D thin section image is used as a training image).

Figure 8 .
Figure 8. E-Type realization obtained from 14 realizations by optimum sized template.

Figure 9 .
Figure 9. 2D continuous training image and simulation result (Optimum Template Size is 15).

Table 2 . Training image preprocessing steps of the Modified FILTERSIM approach.
2. Find optimum template size (Entropy plot).3. Scan TI and store all patterns in patdb along with the frequency of each pattern as the last component of