Updating Local Anisotropies with Template Matching During Geostatistical Seismic Inversion

Iterative geostatistical seismic inversion methods are widely used to predict petro-elastic rock properties from seismic reflection data. When the model perturbation technique uses two-point geostatistics, these methods struggle to reproduce complex and nonstationary geological environments such as faults, folds and highly variable depositional systems. These limitations are often due to the use of a global variogram model to express the expected spatial continuity pattern of the property of interest. In complex geological environments a global variogram model might be unable to detect local heterogeneities and rapid variations of lithology, and result in nonrealistic geological models. Local heterogeneities might be predicted from the data using seismic attribute analysis, which can be imposed during geostatistical seismic inversion as local anisotropy models. In these approaches, the information about the local spatial continuity patterns is fixed and will guide and condition the entire inversion procedure, which can lead to errors and uncertainty in areas where this approach is not appropriate due to high local uncertainty of geological features, given the poor signal-to-noise ratio of the data or the presence of important geological features below the seismic resolution. This work proposes an iterative geostatistical seismic inversion method which iteratively updates the local spatial continuity models based on the trace-by-trace misfit between observed and predicted seismic data. The update of the local spatial continuity models aims at surpassing the limitations of the seismic inversion methods that use a fixed a priori variogram model. The method is successfully illustrated in a challenging two-dimensional synthetic data set and in a real case application. The results demonstrate the benefit of updating iteratively the imposed local spatial continuity patterns based on the data misfit. The inverted models are capable of better predicting the location of faults and reproducing the continuity of sinuous channels.


Introduction
Seismic inversion methods are widely used to predict the spatial distribution of the subsurface geology. They make it possible to transform seismic amplitudes into two-or three-dimensional numerical models of elastic rock properties and facies (e.g., Bosch et al. 2010;Grana et al. 2022). The relationship between seismic data and predicted models can be summarized by where d obs is the seismic data, m is the subsurface properties of interest, and e represents the measurement errors and assumptions about the physical model under investigation during data processing (Tarantola 2005;Tompkins et al. 2011). F is the forward model that provides the mapping between the observed data and the subsurface properties (Tarantola 2005). As F −1 is unknown, predicting the petro-elastic rock properties from seismic data is a challenging nonlinear inverse problem, ill-posed and with multiple solutions (Tarantola 2005). Seismic inversion methods can be divided into two main groups with different assumptions, advantages and limitations (Francis 2006;Bosch et al. 2010;Filippova et al. 2011;Tylor-Jones and Azevedo 2022). Deterministic seismic methods (e.g., Russell and Hampson 1991) can be summarized as optimization procedures that start with an a priori elastic model, which is perturbed to minimize the mismatch (or maximize the similarity) between real and synthetic seismic data. Deterministic methods provide a single best-fit solution, which is a smooth representation of true variable and heterogeneous subsurface geology. Hence, they do not provide any uncertainty assessment of the inverted models.
The second group of seismic inversion methods are statistical-based (Grana et al. 2022). Contrary to deterministic methods, they allow assessing the uncertainty of the inverse solution. Under this framework, Bayesian linearized seismic inversion methods (e.g., Buland and Omre 2003; Grana and Della Rossa 2010;Grana 2016;Grana et al. 2017) assume a multi-Gaussian distribution of the model parameters and the errors within the data, and the linearization of the forward model. With these assumptions, the analytical description of the posterior distribution is also a multi-Gaussian distribution.
In the case of iterative geostatistical seismic inversion methods, the computational burden is due to the use of stochastic sequential simulation and co-simulation methods as model perturbation techniques. This issue has been tackled with the parallelization of the simulation algorithms to graphics processing units (GPU;Ferreirinha et al. 2015) or the use of proxies to reduce the dimensionality of the model parameter space (e.g., Nunes et al. 2019). One of the advantages of geostatistical seismic inversion methods is the reproduction of the first and second statistical moments as computed from the observed data (i.e., the well-log data) and the reproduction of a variogram model used to express the expected spatial continuity pattern of the subsurface property of interest. This capability is ensured by the stochastic sequential simulation and co-simulation methods (e.g., Soares 2001).
Nevertheless, these methods have some limitations when applied to complex sedimentary and/or structural geological environments as they may lead to inverted models that are geologically unrealistic. These limitations are mainly related to the use of a single and global variogram to describe the spatial continuity pattern of the subsurface property, which is unable to represent the spatial distribution of the true geology.
In a Bayesian inversion framework, the issue of nonstationary spatial continuity models has been addressed by including varying spatial covariance matrices (e.g., Madsen et al. 2020;Aune et al. 2013;Bongajum et al. 2013). In iterative geostatistical seismic inversion, this limitation can be overcome by imposing local variogram models, which represent the local spatial continuity pattern in the rock property prediction (e.g., Pereira et al. 2020). The local variogram anisotropy models can be computed a priori and obtained, for example, from seismic attribute analysis of the observed seismic data. However, imposing local variogram models, which are fixed during the entire iterative procedure, can lead to the propagation of errors trough the iterations if the assumed spatial continuity model is unreliable. This might happen in areas of low signal-to-noise ratio, where the local anisotropy variogram models are highly uncertain due to discontinuous seismic reflections.
This work proposes a self-update approach of the local variogram models (i.e., local anisotropy models) by integrating, at each iteration, the mismatch between predicted and observed seismic data through the iterative seismic inversion approach. In the application examples shown herein, the local variogram models are updated using the template matching technique, due to is relatively low computational cost. However, alternative methods could be applied. Iteratively updating this information avoids fixing a single local variogram model that might reflect noise areas of incoherent information. The updated variogram models constrain the stochastic sequential cosimulation of a new set of realizations in the subsequent iterations of the seismic inversion ( Fig. 1).
In the next section, the proposed iterative geostatistical seismic inversion method is described in detail, and the concept of self-updating of local anisotropies using template matching is explained.
To evaluate the performance of the proposed method, the following section shows its application to synthetic and real examples. The self-updating inversion method is benchmarked against an iterative geostatistical seismic inversion (GSI), with a global variogram model, and to geostatistical seismic inversion with local anisotropies (GSI-LA) with fixed local variogram models, estimated a priori from seismic attribute analysis (Pereira et al. 2020). The three methods were compared in terms of the reproduction of the true seismic amplitude, the spatial continuity of the synthetic seismic and acoustic impedance models and corresponding facies prediction. Facies were classified from the acoustic impedance models co-simulated during the last iteration of the geostatistical seismic inversion, using an acoustic impedance (Ip) threshold. Channels with low acoustic impedance values (i.e., sand channels) were classified as sand, and high acoustic impedance values were classified as mud (i.e., background facies). In addition, the Ip models were plotted in a multidimensional scaling space (MDS) computed with a modified Hausdorff distance (e.g., Cox and Cox 1994;Scheidt and Caers 2009;Azevedo et al. 2014).
This section is followed by the discussion of the results and ends with the presentation of the main conclusions.

Methodology
This section details the proposed iterative geostatistical seismic inversion method with self-update local anisotropy models. The proposed method can be divided into four main steps: (1) initial data and stochastic sequential simulation of Ip models; (2) forward modeling; (3) local anisotropy update; and (4) stochastic update of the model parameters with stochastic sequential co-simulation.
As the proposed method builds upon the global stochastic seismic inversion (GSI; Soares et al. 2007), its main steps are first described. This summary is followed by a detailed description of the steps introduced under the scope of this work.

Geostatistical Seismic Inversion
The global stochastic inversion (Soares et al. 2007) inverts full-stack seismic data for Ip. It can be summarized in the following sequence of steps: 1. In the first iteration, stochastic sequential simulation (e.g., direct sequential simulation; Soares 2001) of N realizations of Ip from well-log data and using a global spatial continuity pattern expressed by a unique global variogram model. 2. Computation of N synthetic seismograms by convolving the reflection coefficients obtained from the simulated Ip models (in step 1) with a wavelet. 3. Analysis of the trace-by-trace similarity (S) between the synthetic and real seismic traces as summarized by where x and y are the real and synthetic seismic traces, respectively, with Ns seismic samples. S is bounded between −1 and 1 and expresses how two seismic traces are similar in terms of amplitude content and waveform. If synthetic and real traces are similar, then S is close to 1. 4. From the N realizations of Ip, the elastic traces that produce the highest similarity between synthetic and real seismic are stored in an auxiliary volume with the corresponding S value. 5. In the next iteration, the auxiliary volumes from the previous step are used in the stochastic sequential co-simulation (e.g., Soares 2001) of a new set of N Ip models. 6. Return to step 2 and continue until the global mismatch between the synthetic and the real seismic data is below a given threshold.
All Ip models simulated and co-simulated during the iterative procedure reproduce the well-log data at their locations, and the histogram as retrieved from the Ip-log data. This iterative geostatistical seismic inversion method assumes a single variogram model of Ip, which is normally inferred from the existing well-log data and reproduced in all the realizations simulated and co-simulated during the iterative procedure.
The GSI assumes the stationarity of the spatial continuity pattern, which is hardly valid in complex geological environments. This limitation is explored in the present work by including local variogram models during the model perturbation step. The local variogram models are updated iteratively based on the data misfit computed from the previous iteration. The proposed iterative geostatistical seismic inversion methodology is detailed next.

Geostatistical Seismic Inversion with Self-Update Local Anisotropies
The main idea of the proposed geostatistical seismic inversion consists in updating the local spatial covariances (i.e., local variogram models), based on new secondary information (i.e., the mismatch between observed and synthetic seismic data). Information about local anisotropies is inferred directly from the data, and updated iteratively. The self-updating procedure aims at correcting unreliable estimates about the local spatial continuity patterns, as the iterative inversion converges. The proposed inversion approach uses the geostatistical seismic inversion framework (Sect. 2.1). The inversion starts with a global stationary covariance model C(h), assumed valid and representative of an entire area. The average impedance model of the ensemble of realizations generated at each iteration (i.e., the expected Ip value at each location x 0 ) is scanned by a pattern recognition classification technique to detect and quantify local nonstationary geological patterns, for example, local anisotropic covariances C θ (h), where θ represents the direction angles and anisotropy ratios of the local geological features (e.g., channels, fractures, shale barriers). The following section describes in detail the template matching technique, as an example for pattern recognition classification.
Once local covariance models C θ (h) are modeled, they are allocated to the corresponding spatial location x 0 . At the subsequent iteration, the geostatistical cosimulation is performed with updated local models of covariances (Horta et al. 2010;Caeiro et al. 2015).
The geostatistical seismic inversion with self-update local anisotropies (GSI-LA-Self-Update) can be summarized in the following sequence of steps ( Fig. 1): 1. In the first iteration, N Ip models are generated as described in step 1 for GSI (Sect. 2.1). During the stochastic sequential simulation, a global variogram model inferred from the available well-log data is used. Each of the geostatistical realizations is transformed into synthetic seismograms, and S (Eq. 2) is computed. 2. At the end of the first iteration, the pointwise average (i.e., E-type) model of the generated Ip realizations is scanned by the pattern recognition classification algorithm. In template matching (Sect. 2.3), an ensemble of elliptical templates, T k , with different combinations of anisotropy ratios and azimuths defined a priori, is used. Each template mimics different possible local variogram models with different parameters; the best template match (see Sect. 2.3 for details) is retained at each location and used as local anisotropy in the next iteration. 3. A new set of N Ip models is then co-simulated with local anisotropies (Horta et al. 2010) revealed by the best template obtained in step 2. 4. The iterative procedure repeats steps 2 to 4 until a given number of iterations is reached or the global S between the real and synthetic seismic data is above a given threshold. In the application examples shown herein, the number of iterations was set to six.
The inversion method is introduced to invert full-stack seismic data for acoustic impedance (Ip), but its extension to the elastic domain (e.g., geostatistical seismic amplitude-versus-angle (AVA) inversion) (e.g., Azevedo et al. 2018) is straightforward. In these methods, the templates used to estimate the local spatial continuity models would be deployed in the auxiliary models computed at the end of each iteration, for each elastic property (i.e., density, P-wave and S-wave velocities) individually.

Template Matching and Local Anisotropy Estimation
Templates are commonly used in image processing for feature selection, matching, enhancement or detection (e.g., edge or object detection) (Brunelli 2009). In general, templates correspond to kernels or vectors that can take any geometrical shape. They are used as a moving window that scans an image and tests a set of matching criteria between the image and the template. Template matching is used within the application examples shown herein due to its low computational cost and facility of implementation. Alternative methods for local anisotropy estimation could be used. This selection will impact the results obtained.
In the proposed method, the templates are used in the pointwise average model (E-type) of the N realizations of Ip models generated at a given iteration. At the first iteration, the E-type model is not informative regarding local anisotropies, as all the Ip realizations were generated without conditioning to the seismic data. However, as the iterative procedure assimilates the seismic data, it starts showing the expected spatial patterns observed in the seismic data but with higher spatial resolution, originated by the stochastic sequential simulation. These patterns learn from the data and will appear stronger or weaker depending on the local trace-by-trace similarity S.
The set of templates ( Fig. 2) scans the E-type Ip model, calculated from the ensemble of Ip models generated at a given iteration, with predefined parameters such as direction and anisotropy ratio. At each location, the variance between all the Ip values of a given template is computed. The template with the minimum variance is selected as representative of the local spatial continuity model. This approach acts as a proxy of a local experimental variogram and modeling. The template matching criterion is the minimum variance inside the template T k (Fig. 3) and is summarized by

Fig. 2
Ensemble of predefined elliptical templates T k , mimicking different local variogram models with different parameters (direction, main (La) and minor ratios (Lb) ). The templates set corresponds to all possible combinations between these three parameters where x i are the grid cell values inside the template T k , μ is the mean of these values, and N T k is the total number of grid cells in template T k . The local variogram models, as revealed by the selected templates, will be used in the next iteration to generate the new ensemble of Ip models.

Application Examples
The proposed method was tested in two different cases. The first one corresponds to a two-dimensional synthetic full-stack seismic section and the second to a real full-stack seismic section. Both application examples represent complex geological environments related to deepwater sedimentation and tectonics. The seismic data show patterns of nonparallel seismic reflectors related to the presence of turbidite channels, faults and folded horizons. This complex geology results in a highly anisotropic spatial continuity pattern of the subsurface petro-elastic properties, making the reproduction of these geological features a challenge in geostatistical seismic inversion, when using a global spatial continuity pattern expressed by one global variogram model.

Synthetic Case Application
This synthetic data set mimics a deepwater sedimentary environment corresponding to a turbidite system with typical meandriform channels with low acoustic impedance values (Fig. 4). The background facies is characterized by high Ip values with a strong anisotropy, contrasting with low Ip values modeled with an isotropic variogram model within the channel facies. The spatial distribution of Ip in the background facies was modeled with a spherical variogram model oriented with a constant dip of 27°. The range used for the major direction was set to 20 grid cells and the vertical variogram to 5 grid cells. The channels were modeled with an isotropic variogram model and a large horizontal variogram range of 200 grid cells and 2 grid cells in the vertical direction.
The data set includes a synthetic two-dimensional seismic profile (i.e., true seismic data) with dimensions of 200 by 200 grid cells in the i-and k-directions, respectively (Fig. 4a). The observed seismic data do not include any noise, and the wavelet used to compute it was the same as applied in the inversion (i.e., no uncertainty about the wavelet was included in this study). Noise-free seismic data were considered because the nonstationary nature of the geological background is already a challenge for geostatistical methods based on two-point statistics (i.e., stochastic sequential simulation and co-simulation). Part of the conditioning data set are two wells at positions i = 25 and i = 185 with Ip logs from 50 to 180 ms and 60 to 170 ms, respectively.
The template matching used a set of 96 elliptical templates generated with different azimuths ranging from 0°to 360°with a step size of 15°. The set of templates also included different main and minor variogram range ratios reproducing different possible combinations of local variogram model parameters (Fig. 2). The main range varies with an interval from 12 to 24 grid cells, and the minor range varies between 4 and 8 grid cells, both with a step size of 2. The local spatial continuity models estimated with template matching and used as input in the GSI-LA-Self-Update are shown in Fig. 5. The proposed geostatistical seismic inversion was parameterized with six iterations generating 32 realizations of acoustic impedance models at each iteration. As described above, the simulation procedure was conditioned by a local variogram obtained from template matching and updated after each iteration. The results obtained were compared with those resulting from GSI and GSI-LA.
To ensure a fair comparison of the proposed method with GSI and GSI-LA, all inversions were parameterized equally in terms of data conditioning, number of iterations and simulations per iteration. The main difference between the runs lies in the variogram parameterization. GSI uses a global omnidirectional spherical variogram model with isotropic horizontal ranges of 20 grid cells and vertical range of 5 grid cells, and a predefined local variogram model in GSI-LA computed from seismic attributes (Fig. 6).
The proposed method converged after six iterations, achieving a global similarity S of 0.89 between the true and the predicted synthetic seismic (Fig. 7). The performance of GSI-LA-Self-Update is similar to GSI-LA in terms of global similarity S (Eq. 2). The global convergence of both geostatistical inversion methods with local anisotropy shows that using local information about the expected spatial continuity pattern leads to a significant improvement when comparing to the use of a global variogram model (Fig. 7).

Fig. 5
Evolution of the local spatial continuity models estimated with template matching at the end of each iteration and used in GSI-LA-Self-Update for the next iteration: a-e first iteration to fifth iteration, respectively  Furthermore, the synthetic seismic predicted with the proposed method reproduce the location of the main seismic reflections and the amplitude content of the true seismic data (Fig. 8). However, there are some areas where the proposed method struggled to converge, and hence are areas with higher uncertainty, corresponding to lower local S associated with abrupt changes in the spatial continuity of the true seismic reflections. These regions are also associated with poorer estimations of the local continuity patterns. Figure 5 shows the evolution of the local variogram models through the six iterations. While the first iteration is non-informative, the model predicted at the end of the inversion does reproduce the main spatial features observed in the true one (Fig. 6), with larger dip values.
In seismic inversion, a good match between observed and synthetic seismic is not sufficient to ensure a good model of the subsurface properties. Because of the nonunique nature of the seismic inverse solution, considerably different elastic models might originate similar data. Figure 9 compares the average Ip model computed for the three geostatistical seismic inversion methods here studied. For illustration purposes, the pointwise average Ip model (i.e., E-type model) computed from the stochastic Fig. 9 Vertical sections of Ip mean models from last iteration for the three tested methods: a GSI; b GSI-LA; c GSI-LA-Self-Update. The lower Ip values correspond to meandriform sand-rich channels, and the high Ip values to mud-rich sedimentary sequences. Black arrows point to low Ip values, related to sand channel facies, highlighting the improvement introduced with the proposed method realizations, co-simulated during the last iteration of the inversion, is shown (Fig. 9). The mean Ip model shows a good match in terms of geological spatial continuity pattern, particularly of the meandriform channels, represented by low Ip values. At the large-scale, all the tested methods retrieve Ip models similar to the true one. However, the spatial continuity and connectivity within the channels is more evident with the proposed method of self-update. Also, sharp discontinuities are better predicted when using local information updated iteratively based on the data mismatch (Fig. 9).
To evaluate the predicted models in terms of their spatial continuity, all realizations of Ip from the last iteration resulting from the three inversion methods were classified into two facies, using a single threshold in the Ip domain. Ip values below or equal to 8,500 kPa s/m were classified as sand, corresponding to channel facies, and Ip values above 8,500 kPa s/m were classified as mud, corresponding to muddy background facies. Predicting reliable facies is of outmost importance for volumes calculation and reserve estimation.
The comparison of the most likely facies models with the true facies model for the three methods is shown in Fig. 10. This simple exercise enhances the interpretation of the pointwise Ip models (Fig. 9). The proposed method shows a better local spatial continuity and connectivity of the geological sedimentary patterns, particularly for the low acoustic impedance values, related to sand. In this case, the proposed geostatistical seismic inversion model reproduces channels with increased continuity, a smaller number of artifacts and enhanced connectivity of the predicted channels when compared with the global and fixed local spatial continuity models. While the differences in the classified facies models are not large, they are highly relevant as channel connectivity is of outmost importance for fluid flow modeling and reserve calculation.
To quantify the differences between the predicted ensembles, the Ip volumes cosimulated in the last iteration were plotted in a two-dimensional MDS space using a dissimilarity matrix built with the modified Hausdorff distance. The modified Hausdorff distance selected as it is sensitive to the shape of the features within the model (e.g., Scheidt and Caers 2009). To facilitate the interpretation of the results, the Euclidean distance was computed between each realization generated during the last iteration for the geostatistical seismic inversion methods considered and the true Ip (Cox and Cox 1994;Scheidt and Caers 2009;Azevedo et al. 2014) (Fig. 11). The Ip models predicted by the proposed method show two distinct features: they are closer to the true Ip; and they have a similar range between the minimum and maximum distances. These results illustrate that iteratively updating the local variogram models with information from the data mismatch has a positive contribution (i.e., the models become more reliable), without compromising the exploratory capabilities of the original method (i.e., the variability within the ensemble of predicted models is preserved).
To evaluate the classified models and to explore their misclassification, Table 1 shows a confusion matrix. In this matrix, the diagonal values are the percentage of pixel values correctly classified (i.e., the predicted class is equal to the true class), while off-diagonal values are the percentage of misclassified pixel values. The confusion  The bottom row displays the overall accuracy for each class The values in bold are the percentage of pixel values correctly classified matrix shows that the GSI-LA-Self-Update performs better, estimating the highest total percentage of correctly classified pixel values (89.15%) and the smallest percentage (10.85%) of misclassified pixel values. Among these values, 76.16% of pixel values were correctly classified as mud and 12.99% of pixel values were correctly classified as sand. However, and as expected, the misclassification occurs mostly for sand, where 8.24% of pixel values were classified incorrectly as mud, since this class exhibits a more complex spatial distribution pattern. For the mud class, 2.61% of pixel values were misclassified as sand.

Real Case Application
The proposed method was applied to a real case study aiming at the prediction of the spatial distribution of acoustic impedance from full-stack seismic data. The real seismic data were acquired over a complex geological sedimentary and tectonic setting.
The study area is characterized by sedimentary layering sequences with the presence of turbidite systems, with typical meandering channels. Additionally, tectonic structures are also present in the seismic data and are related to faulting and folding, making it possible to observe several normal faults, tilt and deformation of the sedimentary sequences in the seismic data. The data set corresponds to a two-dimensional vertical seismic section, the inversion grid is composed of 100 by 1 by 206 cells in the i-, j-and k-directions, respectively. Additionally, a well with Ip-log data (Fig. 12b) is located in grid cell i = 74, with a length from grid cell k = 55 to 187, and a 101 ms wavelet (Fig. 12a). For template matching, an ensemble of 108 elliptical templates with azimuths varying between 0°a nd 55°dip with a step angle of 10°was created. The main and minor ellipse ratios were scanned within an interval of 12 to 32 grid cells and 4 to 12 grid cells with a twocell step. The seismic inversion was performed with six iterations, with 32 realizations models of Ip in each iteration.
To evaluate the performance of the proposed method, the results were compared with the results from GSI and GSI-LA, using the same initial data. At the end of the seismic inversion procedure, the proposed method converged, reaching a global similarity S between the real seismic data and the synthetic data of 0.82. GSI and GSI-LA also converged as expected, reaching a smaller global similarity of 0.76 and 0.80, respectively (Fig. 13). As expected, the synthetic seismic data of the proposed method reproduce the real seismic data particularly well (Fig. 14), not only in terms of seismic reflection positions and continuity, seismic textures and corresponding anisotropy variations, but also in terms of the amplitude content. Regarding the sedimentary and tectonic geological structures, from the synthetic seismic data it is possible to interpret the sedimentary sequences including the seismic pattern related to stacked channels linked to a turbidite system. The dipping of the seismic reflectors with some tilting and fault patterns is also well reproduced in the synthetic seismic predicted by the proposed method.
Comparing the proposed method with the synthetic seismic data resulting from the other geostatistical seismic inversion methods (Fig. 14), there is an improvement on the reproduction of seismic reflectors and their discontinuities, providing a more consistent and realistic geological meaning, either sedimentary or tectonic, to the seismic patterns.
Nevertheless, there is some random residual noise, in areas of lower similarity (Fig. 15). The residuals are mainly related to high seismic amplitude values at well log location, which are associated with a poor seismic-to-well tie (Fig. 15). The residuals away from the well location are most likely linked to the representativeness of the wavelet in those locations and to the fact that the data is broadband seismic, introducing some additional uncertainty during wavelet extraction.
The predicted Ip models obtained with the proposed method have some aspects worthwhile to discuss in detail (Fig. 16). The use of local and more informative spatial models allows to highlight the presence of tectonic features, such as normal faults, when compared with the results obtained with GSI considering a global variogram model. Also, the detritical sedimentary patterns, such as turbidite channels, are better identified, characterizing and revealing the location of potential targets. Figure 16c shows a vertical section of the predicted Ip with a fixed local anisotropy model. This Fig. 14 a Vertical section of the true seismic data with well position (black dashed line); b synthetic seismic data calculated from Ip mean models of the sixth iteration GSI-LA-Self-Update; c GSI; d GSI-LA. Black arrows point to areas where it is possible to verify the similarity between the real seismic data and the synthetic data from all tested methods model is greatly discontinuous, indicating the poor a priori estimation of the local spatial continuity patterns, and illustrating the drawbacks of such approach. If the a prior estimation is not reliable, the quality of the inverted models will be poor.
An important issue is related with the ability of the proposed method to update the local spatial patterns (Fig. 17) while keeping the reproduction of the statistics of the true Ip-log data. Figure 12b shows a comparison between both histograms. In Fig. 17, showing the evolution of the local spatial patterns, it is possible to see the update and improvement of the dip azimuth model between iterations with the proposed method, using template matching.

Discussion
The application examples of this study illustrate the improvements obtained in the inverted models with the self-updating of the local spatial continuity pattern of Ip. In such examples, the local continuity models was estimated with template matching, but alternative methods could be used. For example, Pereira et al. (2020) used automatic Fig. 15 Vertical section extracted from the a residuals between synthetic and real seismic; b best-fit Ip; c best-fit similarity models from sixth iteration, with the proposed method. Red arrows point to low Ip values related to turbidite channels or faults variogram fitting at each cell of the inversion grid to estimate a priori local variogram models. However, this is a very computational expensive estimation method when compared to with the template matching proposed here. Template matching proved to be an efficient and easy way to estimate the local variogram models.
The proposed seismic inversion method presents a clear advantage over the alternative methods that use a fixed spatial continuity model. The application examples show that iteratively updating the local spatial continuity models, particularly when the available seismic does not allow to extract reliable local models of anisotropy, results in better predictions in terms of elastic and rock properties. The proposed method improved the geological consistency and plausibility of the inverted models. The local information highlights the presence of geological discontinuities such as faults, folding or tilting of sedimentary sequences, and improves the spatial continuity in the presence of anisotropic depositional features, such as the meandriform channels typical of fluvial-deltaic sedimentary environments or deepwater turbidite systems.
The proposed method was tested in a challenging two-dimensional synthetic data set characterized by sinuous channels built with object modeling. An alternative approach to the proposed seismic inversion method would be the direct inversion of seismic data for facies and the integration of geostatistical simulation conditioned to multi-point geostatistics (e.g., González et al. 2008). These techniques have proven their value in these kinds of geological environments. However, they present additional difficulties related to the need of a reliable training image. The synthetic case application is noisefree. However, the proposed method was also applied to a real data set that has noise Both examples used the inversion method in full-stack seismic for acoustic impedance (Ip), but its extension to the elastic domain (e.g., geostatistical seismic AVA inversion; Azevedo et al. 2018) is straightforward.
Despite all the advantages of GSI-LA-Self-Update, it is important to discuss some limitations. To avoid unnecessary computational costs and to better estimate the local variogram model parameters, it is of utmost importance to choose the most relevant templates combinations at the beginning. To optimize this choice, it is important to start by analyzing the shape, dimension and direction of the main geological features in seismic data. This will indicate the best azimuthal ranges and anisotropy ratios to begin with. It is also important to start with bigger step sizes for azimuthal ranges and with very different anisotropy ratios, and then adjust parameters. It is recommended to use a relatively small, but representative area, to fine-tune the parameters before applying to the full seismic volume. Comparing time consumption for the synthetic case as an Fig. 17 a Vertical section of local spatial patterns-dip azimuth (°) for the second iteration; b vertical section of local spatial patterns-dip azimuth (°) for the fifth iteration, with the proposed method example, the cost of one iteration to compute 32 realizations was 115 min using GSI-LA and 93 min using GSI-LA-Self-Update. Both inversion methods ran on an Intel Core i7 workstation with 32 GB of RAM. Two-dimensional examples were shown here, to facilitate the interpretation of the results, as well as the comparison with the other methods. However, the method can be applied directly to any three-dimensional volume.

Conclusion
This paper proposes an update of the GSI algorithm through a sequential updating of models based on new experimental data. Local spatial covariances are updated by conditioning the seismic inversion process to experimental data-seismic reflection data and local anisotropy models inferred directly from the data.
The new iterative geostatistical seismic inversion method is based on global stochastic seismic inversion, and uses an updatable local spatial continuity model, expressed by a local variogram model, self-updated at the end of each iteration to improve the petro-elastic properties estimation, and thus retrieve more consistent and realistic geological models with less artifacts, particularly for complex tectonic and detrital sedimentary environments.
The method was successfully applied to a synthetic and a real case, with acceptable computational costs. In the examples presented, there is a clear improvement when using not only a local variogram model, but also a self-updatable one, based on the data and during the seismic inversion procedure; this benefits the estimation of local anisotropies of complex geological environments, leading to more robust predictions of the subsurface geology.