SciELO - Scientific Electronic Library Online

Home Pagelista alfabética de periódicos  

Serviços Personalizados



Links relacionados

  • Em processo de indexaçãoCitado por Google
  • Não possue artigos similaresSimilares em SciELO
  • Em processo de indexaçãoSimilares em Google


CT&F - Ciencia, Tecnología y Futuro

versão impressa ISSN 0122-5383

C.T.F Cienc. Tecnol. Futuro vol.8 no.2 Bucaramanga jul./dez. 2018 

Artículos originales



César Ariasa  *  , David Abreob  , Sergio-Alberto Abreob  , Luis Duquea  , Ana-B Ramírezb 

a Centro de Ciencias de la Computación, Instituto Tecnológico Metropolitano, calle 73 No. 76A - 354 - Campus Robledo, Medellín, Colombia.

b Centro de Procesamiento Sísmico, Universidad Industrial de Santander, carrera 27 calle 9 Bucaramanga, Colombia.


Full waveform inversion (FWI) has been recently used to estimate subsurface parameters, such as velocity models. This method, however, has a number of drawbacks when applied to zones with rugged topography due to the forced application of a Cartesian mesh on a curved surface. In this work, we present a simple coordinate transformation that enables the construction of a curved mesh. The proposed transformation is more suitable for rugged surfaces and it allows mapping a physical curved domain into a uniform rectangular grid, where acoustic FWI can be applied in the traditional way by introducing a modified Laplacian. We prove that the proposed approximation can have a wide range of applications, producing precise near-surface velocity models without increasing the computing time of the FWI.

Key words: Full Wave Form Inversion; Reverse Time Migration; Rugged topography; Velocity estimation; Acoustic wave equation


La inversión de onda completa (FWI) ha sido usada recientemente para estimar parámetros del subsuelo, tal como modelos de velocidad. Este método sin embargo, tiene diferentes inconvenientes cuando se aplica a zonas con topografía abrupta debido a la aplicación forzada de una malla cartesiana sobre una superficie curvada. En este trabajo, presentamos una transformación de coordenadas simple que posibilita la construcción de una malla curvada. La transformación propuesta es más apropiada para superficies abruptas y permite mapear un dominio físico curvado a una malla rectangular uniforme, donde la FWI acústica puede ser aplicada de la manera tradicional introduciendo un Laplaciano modificado. Nosotros sugerimos que la aproximación propuesta puede tener un amplio rango de aplicaciones, produciendo modelos de velocidad precisos cerca a la superficie sin incrementar el tiempo de cómputo.

Palabras-clave: Inversión de Onda Completa; Migración en Tiempo Inverso; Topografía abrupta; Estimación de velocidades; Ecuación de onda acústica


Artificial seismic waves are used to study the Earth's inner structure. These waves travel through the subfurface and their echoes are collected by sensors placed at the surface, known as geophones. This information is used to infer the properties of the subsurface such as velocity of propagation of seismic waves or density. Subfurface velocity models were traditionally obtained with first arrival travel time tomography, a technique that uses part of the recorded information, like early arrivals [1]. The advances of computational power in hardware and software enable current methods to consider the full information of the waveform to produce high-resolution properties maps; this is the case of FWI [2]- [5]. After applying this method, migration algorithms can find an accurate location of subsurface reflectors.

FWI is a data-fitting method that uses a rough velocity model of the subsurface and iteratively finds better estimates. The purpose is to iteratively reduce the squared error between the modeled seismogram (computed with the estimated velocity model) and the observed seismogram. FWI can be applied in time domain [3,4,6] or frequency domain [7]; each option implies advantages and disadvantages. The first scheme is more straightforward, faster and requires less memory. It is also sensitive to time sampling, requires more computational time, and may present cycle-skipping problems. On the other hand, in the frequency domain, only one or a few frequencies can be used for FWI. If all the frequencies are simultaneously inverted [7], this approach is equivalent to its time domain counterpart. There are several ways to update the estimate of the velocity model in each iteration, but one of the simplest is the gradient descent method. Furthermore, the strategy proposed by Plessix [8] can be used to compute the gradient via the ad-joint state method. We refer the reader to a comprehensive review for the FWI method written by J. Virieux and S. Operto [5].

Currently, FWI is an active research topic from a mathematical, computational and geophysical standpoint. Several drawbacks of this method still remain as open problems. Other issues are t the inconvenience of avoiding local minima, the sensitivity of the method to amplitude errors [5], and its computational efficiency-related complications. In addition, the problem of incorporating a curved topography into the method is also of great interest as good estimates of model parameters at the near surface are required. This is important in the oil industry because the extractive potential may be significant in onshore zones with complex geology that hide complicated structures. The traditional approach to handle topography is the datuming method [9], which adjusts the amplitudes of the computed wave-field so that the wave-fronts seem to be emitted from points in the curved surface of a complex topography. However, datuming introduces amplitude errors to which the FWI algorithm is very sensitive, leading to potentially inaccurate images of the Earth's subsurface.

Thus, we propose hereunder a method based on a previous work that used generalized meshes for seismic modeling [10]. This work used a general physical domain with curved upper surface and mapped it into a rectangular domain, where a uniform grid was used to propagate and retro-propagate wavefields. The same transformation domain was used to obtain imaging conditions. The general transformation consists of modifying the Laplacian operator by introducing space-dependent coefficients known as metric coefficients. Wave propagation in this domain is similar to that of a Riemannian space, where the metric is not Euclidean. Once the wave-fields and velocity models are computed in the transformed domain, these are translated back into the (curved) physical domain.

Reverse time migration in curved meshes [11] and FWI methods for rugged topography [12] have been proposed before. Unfortunately, their computational cost is greater than that of their Cartesian counterparts. The approaches proposed in the literature, although very clever, are rather unpractical as transformations drastically reduce the spatial step size, which results in limiting the time step size due to the Courant-Friedrichs-Lewy condition.

In this work, we introduce a type of transformation that preserves the spatial size and, therefore, the time step. Hence, the execution time of the algorithm is shorter compared to the same method using Cartesian coordinates. Experimental results show that our approach enables to perform FWI in domains with a general rugged surface. It also estimates satisfying velocity models even for the near-surface region, where obtaining good results is usually difficult.


The transformation that maps a rectangular domain with coordinates (ξ 1, ξ 2) (named computational domain) into the physical domain of coordinates (x 1, x 2) is

where φ(ξ 1)=φ(x1) is a smoothed function that represents the curved upper boundary of the physical domain. This transformation is depicted in Figure 1.

Figure 1 The straight lines of the rectangular domain are transformed into curved lines in the physical domain. For example, the horizontal line on top of the computational domain is mapped into the curved dark line in the physical domain, which represents the topographic profile. 

Note from (1) that Δx1= Δξ 1 and from (2) that if we vary ξ 2 along a vertical line with ξ 1 =const, we have φ=const, which implies Δx2= Δξ 2, i.e., the step size in space is not affected by the transformation.

When the variables in (1) and (2) are changed, a new expression is found for the Laplacian operator by using the chain rule. The acoustic wave equation for the transformed domain is given by

Equation 4 gives the generalized Laplacian, where |g| is the absolute value of the determinant of the metric tensor g ij, given by

In order to re-write the wave equation, we need the contravariant representation of the metric tensor g ij =(g ij) -1 [13], and a sum over repeated indexes is implied. Note that c 2 is the square of the velocity vector, which is a scalar value and, therefore, it is not transformed. However, its arguments are transformed.

Expanding the Laplacian we can re-write Equation 4 in a more convenient way, as


The elements ζ i are also geometric coefficients that are computed only once, as well as g ij . For our specific transformation, given in Equations 1 and 2, we have

(11) introduced a heuristic stability condition for this equation, given by

Now, Equation 3 can be numerically solved in the computational domain, i.e., in a rectangular domain, and the classic FWI algorithm can be implemented.

The purpose of the traditional FWI algorithm is to minimize the quadratic error given by

where ξs=(ξ 1s,ξ 2s=0) are the source positions, that according to Equation 2 corresponds to

x2= φ(ξ 1)=φ(x1), i.e., the mountain border, and the ξ r are the locations of the receivers. In Equation 11, P obs is the wave-field measured by the sensors (geophones in a real on-shore seismogram), and P mod is the modeled seismogram, obtained with the velocity model, c k, estimated at iteration k. The algorithm requires an initial velocity model c0, and it updates the estimate at each iteration by using the gradient descent direction. This is

where G(ck) is the gradient given by Plessix [8]:

In Equation 13, ΔP(T-t,ξ 1,ξ 2) is the retro-propagation of the residual field P mod - P obs . The parameter λ k is a step size that can be optimized by Liu, et. al. [14]


Once the FWI method retrieves the estimated velocity models, it can be mapped into the physical domain by using the transformation given in Equation 2.


We evaluate the performance of the FWI method in generalized coordinates using the Canadian Foothills model, a synthetic velocity model for a zone in British Columbia that shows several complex structures common in that region of Canada. The computational grid of the proposed experiment is 334x200, with ∆x1=∆x2=0.03 km. In every grid point, a receiver is located. The minimum offset is 0.03 km, and the maximum is 8 A km. We used seven point sources in shape of a Ricker wavelet, placed every 40 grid points starting in the grid point number 60. This FWI follows a multi-scale strategy where the sources have central frequencies of f1=5 Hz and f1=15 Hz. In the multi-scale strategy, the estimated velocity model obtained for the first frequency is then used as starting-point for the second frequency. For each wave-propagation, we use a sampling time of 1ms, a total acquisition time of 2.5 s, and the FWI included 200 iterations. The initial velocity model is shown In Figure 2a. The final estimated models for one and two frequencies are shown in Figures 2b and 2c, respectively. As the velocity model originally presents a flat bottom, the curved mesh must go beyond. Therefore, the missing area of the curved mesh Is filled with a constant velocity (6 Km/s) which does not negatively affect the results because those points are mainly immersed in the attenuation layer and do not imply a significant variation of velocities. Additionally, the topographic profile of the original model was smoothed before making any computation. The observed data was generated synthetically using the original model shown in Figure 2d.

Figure 2 (a) The starting velocity model. (b) final velocity model for a source of 5 Hz. (c) final velocity model using a source of 15 Hz and taking the model in (b) as the starting model. (d) original model. 

Figure 3 shows the convergence curve for the final result corresponding to Figure 2c. The residual decreases around 90%, when 100 iterations are performed. A quick decrease is observed for the first 10 iterations and then it flattens up. Figure 4 shows the difference between the inverted model and the original model. Small differences between the two models are observed, except at a few points near the borders of the model. It shows the percentage residual for the last stage in FWI. A difference of less than 10% is observed at the center of the model. Note that although the error bar is between -20% and +20%, these extremes are reached only in small spurious locations.

Figure 3 Convergence curve associated to the 15 Hz stage. 

Figure 4 Percentage difference between final inverted model and the original model (Figure 2d). The color bar shows a scale between -20% and +20%. 


The proposed transformation, in combination with the classical FWI algorithm, produces results with small errors as compared to the original synthetic model.

The curved grid used allows the FWI algorithm to map, in a clear way, the near-surface structures of the velocity model.

Other curved grids applied to simple velocity models were reported previously by other authors. Our method was applied to a complex model that exhibits imbricate geological structures and curved topography. The resulting velocity model shows good agreement with the original, within a small percentage difference.


This work is supported by ECOPETROL and COLCIENCIAS, as a part of research project grant No. 0266-2013.


[1] Kubota, R., Nishiyama, E., Murase K. & Kasahara J. Travel time estimation of first arrivals and later phases using the modified graph method for a crustal structure analysis. Exploration Geophysics. 2009, 40(1), 105-113. [ Links ]

[2] Lailly, Patrick, and J. B. Bednar. (1983)."The seismic inverse problem as a sequence of before stack migrations." Conference on inverse scattering: theory and application. Philadelphia, PA: Siam. [ Links ]

[3] Tarantola, A. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 1984, 49(8), 1259-1266. [ Links ]

[4] Mora, P. Nonlinear two-dimensional elastic inversion of multioffset seismic data. Geophysics , 1987, 52(9), 1211-1228. [ Links ]

[5] Virieux, J., & Operto, S. (2009). An overview of full waveform inversion in exploration geophysics. Geophysics , 2009, 74(6), WCC1-WCC26, [ Links ]

[6] Bunks, C Saleck, F.M., Zaleski, S & Chavent, G. Multiscale seismic waveform inversion, Geophysics , 1995, 60(5), 7-1473. [ Links ]

[7] Pratt, R. G. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics , 1999, 64(3), 888-901. [ Links ]

[8] Plessix, R.-E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysics Journal International, 2006, 167(2), 495-503. [ Links ]

[9] Berryhill, J. Wave-equation datuming. Geophysics , 1979, 44(8), 1329-1344 [ Links ]

[10] Shragge, J. Solving the 3D acoustic wave equation on generalized structured meshes: A finite-difference time-domain approach. Geophysics , 2014, 79(6), 1-16. [ Links ]

[11] Shragge, J. Reverse time migration from topography. Geophysics , 2014, 79(4), 1-12. [ Links ]

[12] Han, Yutong, et al. (2015). Acoustic full waveform inversion method research in areas with a rugged surface. En SEG Technical Program Expanded Abstracts 2015. Society of Exploration Geophysicists. [ Links ]

[13] Synge, John Lighton, and Alfred Schild. Tensor calculus. 1978, Vol. 5: Courier Corporation. [ Links ]

[14] Liu, Faqi, et al. (2012). 3-D time-domain full waveform inversion of a Valhall OBC dataset. En SEG Technical Program Expanded Abstracts 2012. Society of Exploration Geophysicists. [ Links ]

Received: April 13, 2018; Revised: July 06, 2018; Accepted: September 13, 2018

* email:

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License