Abstract
Seismic imaging and inversion become extremely challenging when dealing with salt structures. Conventional state-of-the-art full-waveform inversion (FWI) fails to recover those features in areas where salt is present. A widely used solution in industry, however, involves substantial human interpretation. In this study, a regularized isotropic full-waveform inversion that penalizes the velocity drops in depth is introduced to recover the top parts of salt bodies. Then, an automatic salt flooding is applied to reconstruct deeper parts of the salt. Finally, in order to improve the accuracy of the inverted model, which is strongly affected by anisotropy, an anisotropic FWI is used to update the velocity model. We tested the approach on both synthetic and field datasets. Our FWI results revealed satisfying salt recovery as well as detailed velocity features in areas close to salt bodies.
1 Introduction
A typical salt dome containing halite and anhydrite rocks is mostly impermeable. As it moves upwards, it bends the strata of existing rocks, which forms pockets where petroleum and natural gas are usually trapped. As a result, the delineation of salt and subsalt is of high interest in oil and gas production. Tomography is a standard practice for salt delineation where the best possible sediment velocity model is first built. Then, a top-of-salt interpretation using sediment flood migration followed by a base-of-salt interpretation using salt flood migration is carried out. For complex areas with complex salt geometry, a number of salt scenario tests are also necessary, which is not only labor-intensive and time-consuming but also prone to human misinterpretation. Alternative approaches, based on artificial intelligence (AI), are proposed to build a salt velocity model (; Rusmanugroho et al., 2022).
Seismic inversion plays a key role in delineating subsurface structures (; ; ; ). The full-waveform inversion (FWI) technique, which is considered the most promising data-driven tool to iteratively build velocity models, was first proposed by ), which aims to build high-quality and resolution velocity models with high precision for seismic imaging. However, the presence of complex salts causing high velocity contrast fails to provide high-resolution velocity model building from seismic data. This is caused by various factors such as strong multiples and insufficient subsurface illumination caused by velocity complexities in the overburden. demonstrated that it is possible to use FWI to build a salt velocity model from a smooth initial model. However, the frequency they used was as low as 0.5 Hz, which is much lower than the real case in field data. Since then, traditional processing procedures involve calculating various salt attributes (; ; ) to detect salt boundaries. These attempts often require both geophysical processing skills and geological interpretation (; ; ). and applied FWI to successfully correct some salt structure misinterpretation based on the availability of low frequency and full azimuth of ocean-bottom node (OBN) data, which greatly improved subsalt delineation in the Gulf of Mexico (GOM). Zhang et al. (2018) proposed a time-lag FWI (TLFWI) which uses a time-lag cost function (; ) to reduce amplitude errors in the standard data-mismatch cost function (). further applied TLFWI to velocity model building for areas with the presence of salt bodies.
The edge-preserving technique such as total variation (TV), first proposed by and introduced in image restoration by Vogel and Oman (1998), was used in the estimation of density contrast distribution from gravity data inversion (; ). Compared to other sediments, salt bodies are steeply dipping complex-shaped structures with higher velocity and lower density. Based on this, proposed the application of TV to acoustic velocity perturbation estimation. These edge-preserving techniques usually are very successful in salt top recovery. However, they would fail to recover the bottom of the salt body and the subsalt due to erroneous velocity estimations in the vicinity of the salt body. Therefore, some additional techniques are required to recover the rest of the salt body (especially the deeper part of the salt). proposed the automatic salt flooding into salt delineating, by which the position of bottom salt can thus be picked up precisely. further applied TV to recover the sharp boundaries of salts. Additionally, he also implemented salt flooding to recover the bottom of salt bodies.
In this paper, we further implement this technique into anisotropic media by proposing a new gradient calculation equation and applying it into the misfit function with TV regularization and flooding terms. This study solves two cascading minimization problems to determine the velocity field from the salt-affected dataset, where an edge-preserving technique is first applied to invert the top part of the salt (using the regularized FWI). Then, the output of the edge-preserving technique is used as an input to recover the rest of the salt (using the automatic salt flooding). Finally, we test this approach on both synthetic and real datasets. Our FWI results have revealed satisfying salt recovery as well as detailed velocity features in areas close to salt bodies.
2 Theory
2.1 Isotropic FWI + TV
The minimization of the misfit function with regularization can be expressed aswhere m is the discretized Earth subsurface parameters defined at all space points, Jd (m) is the objective function, d0 is the observed seismic data, dm is the modeled seismic data, and TV is the regularization parameter. In the aforementioned equation, TV is a positive regularization parameter that controls the trade-off between the data misfit term and the amount of interface-preservation in the inversion process. When we use low-frequency data for isotropic FWI, the sharp edges and discontinuities, such as the top of the salt structures, are preserved. Furthermore, ||m||TV can be denoted as
2.2 Cycle skipping
However, if low-frequency data are missing or low-wavenumber components of the initial velocity model are inaccurate, the adjoint source will suffer from cycle-skipping issue, which occurs when the observed and predicted data mismatch by more than half a cycle (; ; ; ). As a result, it leads the inversion to be trapped in a local rather than the global minimum. The adjoint source is shown in Figure 1; for example, when the arrival time is less than 40 ms, there is barely cycle skipping. With the increase in arrival time, cycle skipping can be detected and the adjoint source is severely affected.
FIGURE 1
Cycle-skipped local minima are a major challenge for robust FWI. Theoretically, low frequencies are crucial for FWI to avoid cycle skipping because at low frequencies, both observed and modeled data should match within the half wave-cycle. To overcome cycle skipping, multiscale approaches are designed to update the model using the frequency bands from low to high (; ).
Additionally, some other approaches, including modifying objective function, are also widely used to mitigate cycle skipping, such as cross-correlation, deconvolution, and Wiener filter. In this paper, the automatic cycle-skipping detection approach is applied, and then cycle-skipped energy is removed by modifying the adjoint source with a low pass filter. The adjoint source and the corresponding frequency distribution at different arrivals are shown in Figure 2. The classic Ricker wavelet is used as a reference source, and its frequency distribution is shown in Figure 2 a.1 and a.2. When the arrival time difference between the observed and modeled data is 20 ms (Figure 2 b.1–b.2) and 40 ms (Figure 2 c.1–c.2), respectively, there is no cycle-skipping frequency detected. Cycle skipping can be detected when the arrival time difference reaches 60 ms (Figure 2 d.1–d.2), where a notch appears at 16 Hz. Therefore, we design a notch detection algorithm, using which frequency components exceeding the notch point (notch frequency) can be filtered out. The highest frequency without cycle skipping can thus be estimated. The adjoint source of the L2 norm waveform difference FWI becomes
FIGURE 2
where f is a zero-phase low pass filter. It is similar to applying a low pass filter on the gradient to remove the high wavenumber local minimum. With a smoothed starting velocity model, the simulated shot gather of the first few iterations usually only contains diving waves, and the notch frequency can be easily detected from the residual between the synthetic and observed data. When reflections are involved, the automatic detection of the notch frequency becomes more complex due to conflicting events.
2.3 Automated flooding
A misfit function with Automated flooding can be expressed aswhere JF (q) is the objective function, q is the updated model, u is the model from regularized FWI (FWI + TV), and F is the regularization parameter. In the aforementioned equation, it starts from the final output of isotropic ‘FWI + TV′ and penalizes the velocity drop in the model. Since the Automated flooding is thus applied everywhere, the high velocity is smeared across any region with a drop in the velocity field rather than only the salt body region. Therefore, in areas with low velocity, Automated flooding might have a detrimental effect. In that case, a minimal human intervention is required to flood from the top of the salt.
2.4 Anisotropic TTI FWI
Based on TV regularization and automated flooding in conjunction with isotropic FWI using low-frequency data, the geological shape of the salt body is anticipated to be well recovered. To obtain more detailed layered structures, an anisotropic TTI FWI is then applied. Although the FWI procedure is, in essence, an optimization problem to find the minimum of the objective function, a precise approximation of the anisotropic wave equation with less interactive noise from different wave modes is still considered the core of the objective function, which generates the modeled data that are used to calculate the data difference between observed data and will be ultimately contributed to the objective function. The approximated anisotropic q-P wave equation () we used in this study iswhere , , and are new wavenumbers in a rotated coordinate system and the lower case p used in this equation denotes the scalar wavefield of the P wave component. vh and vp0 are P wave velocity in horizontal and vertical direction, respectively; vn is the normal move-out velocity; kx, ky, and kz are wavenumbers in the original coordinate system without rotation; and is the Thomsen parameter.
The gradient of vertical P wave velocity, , can be expressed as
A real example is an OBN dataset from a Malaysian field that suffers from anisotropy, and a gas cloud imaging problem is used to validate our TTI FWI. Shallow gas layers are stacked on top of each other, which makes the travel path of the wavefield more complicated. A wipe-out zone can be seen under the gas layers, as shown in Figure 3A. The lowest frequency in these data is 7 Hz (below which the data suffered from low signal-to-noise ratio), as shown in Figure 3B. The 3D TTI FWI starts from 7 Hz to 12 Hz with a total of 45 iterations.
FIGURE 3
Figure 4A shows the initial velocity model obtained by smoothing the tomography velocity model. The updated velocity model after running 45 iterations of TTI FWI is shown in Figure 4B. Compared with the initial velocity model, the updated velocity below the gas zone reveals more detailed geological layers, which matches the acoustic well logging data. In deeper parts, even more layered structures can be detected, where a fault-like structure is distinctively present, as shown in Figure 4B.
FIGURE 4
3 Examples
3.1 Synthetic example
The SEG Advanced Modeling (SEAM) 3D salt model was generated by SEAM Corporation as the SEAM phase I project. SEAM’s mission is to produce Earth models and simulations on those models to advance the state of the art in multiple disciplines over a multiyear period. To be useful in multiple disciplines, the SEAM Subsalt Earth Model was designed to capture as much physics and realism as possible in a 3D model that was relevant to oil and gas exploration, including stratigraphy and reservoirs. In the first example, we discrete the SEAM velocity model into 588 × 672 × 212 cells, where the cell sizes are all 60 × 60× 30 m. The orthogonal seismic acquisition has altogether 3,822 shots and 388,278 receivers, with a maximum offset of approximately 10 km. Only a part of the SEAM model is used to reduce the computation cost of 3D TTI simulation and FWI. The selected part is also re-gridded to be 20 × 20 × 10 m, with a maximum depth of 4,500 m. The size of the selected model is 10.5 × 12 × 4.5 km, where 1,176 shots are recorded and a total of 32,021 receivers are regularly placed on the surface for each shot. The starting frequency is 2 Hz. Figure 5A shows the true velocity model. The initial model is obtained by smoothing the true velocity model, as shown in Figure 5B.
FIGURE 5
As a first step, the TV-regularized isotropic FWI with salt flooding is applied to produce and preserve the top of the salt body. Figure 6A shows an inverted velocity model after 30 iterations. Compared with the initial velocity model in Figure 5B, after TV regularization and the automated salt flooding, the shape of the salt (especially the top of the salt) is mostly recovered. However, most of the detailed structures are still missing, such as the clear water bottom and sedimentary layers. Therefore, an anisotropic FWI based on the new q-P approximation equation is applied for higher frequency iterations. To reduce the computational time, we only update a specific area of the SEAM model. The inversion result after anisotropic FWI is shown in Figure 6B and compared with the inverted velocity update after TV regularization and salt flooding as shown in Figure 6A, more detailed fine-layered structures below the water bottom line and at deeper parts are delineated, as circled out in Figure 6B.
FIGURE 6
Finally, a cross-section in the middle of SEAM is illustrated in Figure 7, the inverted model is shown in Figure 7B, and more detailed fine-layered structures at deeper parts are delineated, compared with the initial model, which is shown in Figure 7A. Figure 7C shows the true velocity model as a benchmark. The inverted and true velocity models match with each other quite well. Moreover, two validation wells (well A and well B) are demonstrated in the true model for further validation of our methodology. Figure 8 shows the comparison between initial, inverted, and true velocity models for both wells A and B, where the blue curve denotes the initial velocity and the curve in red denotes the inverted velocity. The inverted velocity model matches with the true velocity model, especially in well B, and our inversion reaches a precision of 0.986 in well B.
FIGURE 7
FIGURE 8
3.2 Field data example
For the real data example, a marine seismic dataset acquired in a region of Gulf of Mexico (GOM) is applied. The area is divided into three distinct tectonic regions: compressional, extensional, and salt-dominated. The data are recorded using a short streamer of length 4.8 km. The lowest usable frequency in the dataset is approximately 6 Hz. We consider the model of interest to be 18 km × 22 km × 8 km deep. Figure 9A shows the velocity model acquired from tomography, in which the salt body is well demonstrated; however, some detailed structures are not well recovered. In that case, we will just use its salt as the reference salt body for our FWI procedure. We first smooth the velocity model to destroy the salt shape, as shown in Figure 9B. The water layer velocity is kept unchanged throughout the inversion scheme.
FIGURE 9
Furthermore, we start the inversion from 6 Hz (due to the lack of low-frequency data) up to 12 Hz, and the final inverted velocity model is shown in Figure 10. Compared with the initial model, the shapes of the salt body have been well updated. In addition, some detailed layered structures are also recovered. For this field data example for FWI, it is noted in the inverted model that the bottom of the salt might not be well captured. The ability of the inversion model to correct the velocity below the salt hinged on the availability of long model wavenumber updates capable of reversing the flooding. In this dataset, the lack of low frequencies and the limited available offsets contributed to this limitation.
FIGURE 10
4 Conclusion
In this study, we proposed a regularizing anisotropic FWI by penalizing the velocity drops with depth. For the first stage, we apply a regularized isotropic FWI (‘FWI + TV′) to match the low-frequency data, so as to preserve the top and sharp edges of the salt body. For the second stage, an automated flooding is applied to flood from the top of the salt to recover the deeper parts of the salt body. After the salt body is mostly recovered, the final stage that applies anisotropic TTI FWI using higher-frequency bands is used to invert more detailed structures. This split approach is extremely convenient because it is easy to implement in existing FWI schemes. This procedure is mathematically more robust and less cumbersome than the standard “migrate-pick-flood” method. Furthermore, the ‘TTI FWI′ ensures the delineating of subtle detailed structures near the salt body. We have shown the application of the proposed method to parts of the 3D SEAM model, in which the lowest available frequency of the observed synthetic datasets is 2 Hz. In addition to synthetic examples, this study has considered a challenging field dataset to invert for salt provinces.
Statements
Data availability statement
The data analyzed in this study is subjected to the following licenses/restrictions: The field datasets analyzed for this study solely belong to PETRONAS and cannot be made publicly available for legal or ethical reasons. Requests to access these datasets should be directed to junxiao.li@petronas.com.my.
Author contributions
JL and HR contributed to the conception and design of the study. JL contributed to anisotropic FWI and field data processing and wrote the first draft of the manuscript. MK and HR contributed to 3D TV and automated flooding. KX contributed to the cycle-skipping research. FD contributed to synthetic data application. HR, JL, and KX contributed to manuscript revision and read, and approved the submitted version.
Acknowledgments
The authors would like to thank PETRONAS for providing permission to show seismic/well data from the Malay basins and allowing the publication of this work.
Conflict of interest
Authors JL, HR, and FD were employed by the company PETRONAS Research Sdn Bhd, Exploration Technology, Group Research and Technology; author MK was employed by the company S-Cube; and author KX was employed by the company CNPC BGP.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
AlkhalifahT. (2014). Full waveform inversion in an anisotropic world: Where are the parameters hiding?EAGE Publications.
2
AnagawA. Y.SacchiM. D. (2012). Edge-preserving seismic imaging using the total variation method. J. Geophys. Eng.9 (2), 138–146. 10.1088/1742-2132/9/2/138
3
Araya-PoloM.JenningsJ.AdlerA.TaylorD. (2018). Deep-learning tomography. Lead. Edge37 (1), 58–66. 10.1190/tle37010058.1
4
BerthelotA.SolbergA. H. S.GeliusL-J. (2013). Texture attributes for detection of salt. J. Appl. Geophys.88, 52–69. 10.1016/j.jappgeo.2012.09.006
5
BorisovD.SinghS. C. (2015). Three-dimensional elastic full waveform inversion in a marine environment using multicomponent ocean-bottom cables: A synthetic study. J. Int.201, 1215–1234. 10.1093/gji/ggv048
6
BrendersA. J.PrattR. G. (2007). Full waveform tomography for lithospheric imaging: Results from a blind test in a realistic crustal model. Geophys. J. Int.168 (1), 133–151. 10.1111/j.1365-246x.2006.03156.x
7
ChapmanC. H.PrattR. G. (1992). Traveltime tomography in anisotropic media: Theory. Geophys. J. Int.109, 1–19. 10.1111/j.1365-246x.1992.tb00075.x
8
ChaventG.ClémentF.GómezS. (1994). “Automatic determination of velocities via migration-based traveltime waveform inversion: A synthetic data example,” in SEG technical program expanded abstracts 1994 (Society of Exploration Geophysicists), 1179–1182.
9
EsserE.GuaschL.HerrmannF. J.WarnerM. (2016). Constrained waveform inversion for automatic salt flooding. Lead. Edge35 (3), 235–239. 10.1190/tle35030235.1
10
FarmerP.MillerD.PieprzakA.RutledgeJ.WoodsR. (1996). Exploring the subsalt: Oilfield review, 8.
11
FarquharsonC. G. (2008). Constructing piecewise-constant models in multidimensional minimum-structure inversions. Geophysics73 (1), K1–K9.
12
JingZ.ZhangY.ChenZ.LiJ. (2007). “Detecting boundary of salt dome in seismic data with edge detection technique,” in 2007 SEG annual meeting (OnePetro).
13
JonesI. F. (2012). Tutorial: Incorporating near-surface velocity anomalies in pre-stack depth migration models. First Break30, 47–58. 10.3997/1365-2397.2011041
14
KalitaM.KazeiV.ChoiY.AlkhalifahT. (2019). Regularized full-waveform inversion with automated salt flooding. Geophysics84 (4), R569–R582. 10.1190/geo2018-0146.1
15
LaillyP. (1984). The seismic inverse problem as a sequence of before stack migration: Conference on inverse scattering: Theory and application, 206–220.
16
LeeuwenT.HerrmannF. J. (2013). Mitigating local minima in full-waveform inversion by expanding the search space. Geophys. J. Int.195, 661–667. 10.1093/gji/ggt258
17
LiJ.XinK.DzulkefliF. S. (2022). A new qP-wave approximation in tilted transversely isotropic media and its reverse time migration for areas with complex overburdens. Geophysics87 (4), S237–S248. 10.1190/geo2021-0433.1
18
LuoY.SchusterG. T. (1991). Wave-equation traveltime inversion. Geophysics56 (5), 645–653. 10.1190/1.1443081
19
MartinsC. M.WilliamsA.LimaV. C. B.Silva.J. C. (2011). Total variation regularization for depth-to-basement estimate: Part 1—mathematical details and applications. Geophysics76 (1), I1–I12. 10.1190/1.3524286
20
McCannD.ComasA.MartinG.McGrailA.LeveilleJ. (2012). Seismic data processing empowers interpretation: A new methodology serves to mesh processing and interpretation. Hart’s E&P.
21
MichellS.ShenX.BrendersA.DellingerJ.AhmedI.FuK. (2017). “Automatic velocity model building with complex salt: Can computers finally do an interpreter's job?,” in 2017 SEG international exposition and annual meeting (OnePetro).
22
PanW.InnanenK. A.MargraveG. F.FehlerM. C.FangX.LiJ. (2016). Estimation of elastic constants for HTI media using Gauss-Newton and full-Newton multiparameter full-waveform inversion. Geophysics81 (5), R275–R291. 10.1190/geo2015-0594.1
23
PanW.InnanenK. A.YuG. (2018). Elastic full-waveform inversion and parametrization analysis applied to walk-away vertical seismic profile data for unconventional (heavy oil) reservoir characterization. Geophys. J. Int.213, 1934–1968. 10.1093/gji/ggy087
24
PrattR. G.ShinC.HickG. (1998). Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion. Geophys. J. Int.133, 341–362. 10.1046/j.1365-246x.1998.00498.x
25
RudinL. I.StanleyO.Fatemi.E. (1992). Nonlinear total variation based noise removal algorithms. Phys. D. nonlinear Phenom.60 (1-4), 259–268. 10.1016/0167-2789(92)90242-f
26
RusmanugrohoH.LiJ.Daniel Davis MuhammedM.SunJ. (2022). “3D velocity model building based upon hybrid neural network,” in SEG/AAPG international meeting for applied geoscience & energy (OnePetro).
27
ShenX.AhmedI.BrendersA.DellingerJ.EtgenJ.ScottM. (2017). “Salt model building at Atlantis with full-waveform inversion,” in SEG technical program expanded abstracts 2017 (Society of Exploration Geophysicists), 1507–1511.
28
ShinC.MinD-J. (2006). Waveform inversion using a logarithmic wavefield. Geophysics71 (3), R31–R42. 10.1190/1.2194523
29
TarantolaA. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics49, 1259–1266. 10.1190/1.1441754
30
VogelC. R.OmanM. E. (1998). Fast, robust total variation-based reconstruction of noisy, blurred images. IEEE Trans. image Process.7 (6), 813–824. 10.1109/83.679423
31
WangP.ZhangZ.MeiJ.LinF.HuangR. (2019). Full-waveform inversion for salt: A coming of age. Lead. Edge38 (3), 204–213. 10.1190/tle38030204.1
32
WangZ.HegazyT.LongZ.AlRegibG. (2015). Noise-robust detection and tracking of salt domes in postmigrated volumes using texture, tensors, and subspace learning. Geophysics80 (6), WD101–WD116. 10.1190/geo2015-0116.1
33
WuLuoR. J.WuB. (2014). Seismic envelope inversion and modulation signal model. Geophysics79, WA13–WA24. 10.1190/geo2013-0294.1
34
ZhangZ.MeiJ.LinF.HuangR.WangP. (2018). “Correcting for salt misinterpretation with full-waveform inversion,” in 88th Annual International MeetingSEG, Expanded Abstracts, 1143–1147. 10.1190/segam2018-2997711.1
Summary
Keywords
full-waveform inversion, TV regularization, automatic salt flooding, anisotropy, salt dome
Citation
Li J, Rusmanugroho H, Kalita M, Xin K and Dzulkefli FS (2023) 3D anisotropic full-waveform inversion for complex salt provinces. Front. Earth Sci. 11:1164975. doi: 10.3389/feart.2023.1164975
Received
13 February 2023
Accepted
27 March 2023
Published
14 April 2023
Volume
11 - 2023
Edited by
Wenyong Pan, Institute of Geology and Geophysics (CAS), China
Reviewed by
Huaizhen Chen, Tongji University, China
Bing Wang, China University of Petroleum, China
Qi Hu, University of Calgary, Canada
Updates
Copyright
© 2023 Li, Rusmanugroho, Kalita, Xin and Dzulkefli.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Junxiao Li, junxiao.li@petronas.com.my
This article was submitted to Solid Earth Geophysics, a section of the journal Frontiers in Earth Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.