An adaptive unstructured mesh based solution to topography least-squares reverse-time imaging

Least-squares reverse-time migration (LSRTM) attempts to invert for the broadband-wavenumber reflectivity image by minimizing the residual between observed and predicted seismograms via linearized inversion. However, rugged topography poses a challenge in front of LSRTM. To tackle this issue, we present an unstructured mesh-based solution to topography LSRTM. As to the forward/adjoint modeling operators in LSRTM, we take a so-called unstructured mesh-based “grid method”. Before solving the two-way wave equation with the grid method, we prepare for it a velocity-adaptive unstructured mesh using 5 a Delaunay Triangulation plus Centroidal Voronoi Tessellation (DT-CVT) algorithm. The rugged topography acts as constraint boundaries during mesh generation. Then, by using the adjoint method, we put the observed seismograms to the receivers on the topography for backward propagation to produce the gradient through the cross-correlation imaging condition. We seek the inverted image using the conjugate gradient method during linearized inversion to linearly reduce the data misfit function. Through the 2D SEG Foothill synthetic dataset, we see that our method can handle the LSRTM from rugged topography. 10 Copyright statement. The copyright notice will be updated ...


Introduction
Seismic imaging is a geophysical technique that investigates the subsurface using observed sesimograms.Rugged topography is an important issue in exploration geophysics, which is usually present with low-velocity weathered layers, high-velocity outcrops of steep dips.Conventional static corrections (Shtivelman and Canning, 1988) may produce imperfect imaging results when complicated geological structures with rugged topography occur.In this case, reverse-time migration (RTM) (Whitmore, 1983) comes into being a proper solution.
Based on the two-way wave equation, RTM has proven to be an important tool in imaging complex subsurface structures without dip limitation (Whitmore, 1983;Baysal et al., 1983).RTM has also got applications to ground-penetrating radar imaging (Foroozan and Asif, 2010;Liu et al., 2017a), medical imaging (Kosmas and Rappaport, 2006;Wang et al., 2016) and waveform inversion (Van Leeuwen and Mulder, 2010;van Leeuwen and Herrmann, 2013;Van Leeuwen and Herrmann, 2013;Kadu et al., 2017).However, in practice, the survey zone usually has the topography effect (Burgin et al., 2014).Several papers 1 Solid Earth Discuss., https://doi.org/10.5194/se-2019-37Manuscript under review for journal Solid Earth Discussion started: 11 March 2019 c Author(s) 2019.CC BY 4.0 License.on the topography RTM applications have been published.McMechan and Chen (1990) produce a topography RTM image by injecting the traces onto the nodes at different elevations within the finite-difference stencil.Rajasekaran and McMechan (1995) further promote this scheme to real data.However, the regular finite difference stencil cannot get rid of the staircase approximation.To overcome this limitation, several approaches considering the topography effect are proposed.Lan et al. (2014) cope with the rugged topography by flatting it to create a curvilinear mesh.Shragge (2014) deals with RTM from topography with the generalized coordinate system.However, these two techniques involve a "flattening" strategy, which cannnot handle a very rugged topography.Liu et al. (2017b) propose to deal with rugged topography RTM with adaptive, unstructured triangular meshing.RTM, however, still acts as the adjoint operator to the linearized forward operator (Born modeling) rather than the inverse operator (Claerbout, 1992), leading to imperfect images.To overcome this imperfectness, LSRTM attempts to approach the inverted operator via linearized inversion.The concept of least-squares migration (LSM) dates back to Lailly (1983).To date, LSRTM has got widespread attention and development.Wong et al. (2011) use LSRTM to image the ocean-bottom data.Dai and Schuster (2013) improve the computational efficiency of LSRTM by exploiting source-encoding.Zhang et al. (2015) propose a practical LSRTM solution by using a cross-correlation objective function.Yao and Jakubowicz (2016) formulate LSRTM in a generalized matrix form.Liu (2016) presents that eliminating the redundant effect of source wavelet in RTM can improve the convergence of LSRTM.Liu et al. (2017c) introduce a prestack cross-correlative LSRTM approach.Xu and Sacchi (2017) formulate the exact adjoint operator to the forward modeling to improve the performance of LSRTM.Liu and Peter (2018) introduce a one-step data-domain LSRTM.
In this paper, we explore to handle LSRTM from rugged topography by using a modeling solver named "grid method" (Zhang and Liu, 1999;Zhang, 2004;Gao and Zhang, 2006) based on an adaptive unstructured mesh.The computational cost and the stability of the "grid method" are similar with those of finite-difference method at O(dt 2 , dx 2 ) (Zhang and Liu, 1999).
The unstructured mesh, which is velocity-adaptive, is generated before wavefield simulation.During the mesh generation, the rugged topography serves as boundary constraints.The meshing is based on the Delaunay Triangulation plus Centroid Voronoi Tessellation (DT-CVT) algorithm (Du et al., 1999;Du and Gunzburger, 2002;Du et al., 2006).The generated adaptive mesh has coarser and finer grids for higher and lower velocities, respectively.However, for a given modeling method, the adaptive unstructured mesh can help run the seismic simulation at a lower computational cost, which alleviates the computational burden of topography LSRTM.When the mesher and solver are ready, we begin to run the LSRTM workflow with MPI (Message Passing Interface) in parallel.
We organize the paper as follows.First, we review the grid method.Then, we present the adaptive, unstructured-mesh generation, with rugged topography as constraints.Next, we introduce the topography LSRTM workflow.Finally, we validate our methods on the Foothill model. 2 Theories and Methods

Solver: the grid method
As an unstructured mesh-based approach, the grid method can model seismic wave propagation in heterogeneous media by accurately considering the surface and interface topographies (Zhang and Liu, 1999;Zhang, 2004;Gao and Zhang, 2006).Let us start with the 2D acoustic wave equation where c (x, z) is the acoustic velocity and P (x, z, t) the pressure.Fig. 1 illustrates the local mesh in the grid method.We assign the velocity parameters c(x, z) at the centroids of the cells, and the pressure P (x, z, t) at the nodes.The dash lines link the cell centroids and the edge midpoints.Given node k centered in the domain Ω k circles by contour 1-2-3-4-5-6-7-8-9-10-11-12-1, we integrate both sides of Eq. (1), yielding ¨Ωk 1 with S kl being the curve of the l cell around node k, and m the number of cells around node k.For example, in Fig. 1, node k is enclosed by contour 1-2-3-4-5-6-7-8-9-10-11-12-1 within its surrounding six cells.
After taking M k = ˜Ωk 1 c 2 dxdz inside domain Ω k , the left-hand side of Eq. ( 2) reduces to M k ∂ 2 P ∂t 2 k .For a typical cell ijk, we calculate the first-order triangular-difference operators Cook et al. (2007) as with ∆ being the area of triangle ijk, , and so on.Substituting Eq. ( 3) into the right-handside of Eq. ( 2) yields Then, we obtain a weak form of Eq. ( 1) represented around node k as With an O dt 2 approximation to ∂ 2 P ∂t 2 k , we update the pressure from t and t − dt to t + dt using In numerical implementations, we calculate and store the parameters M , a, b and ∆ in advance.We only keep the pressures on the nodes, the number of which is twice smaller than that of the cells.In terms of the computational cost and the accuracy and stability, this method is similar with the O(dt 2 , dx 2 ) finite-difference method (Zhang and Liu, 1999).

Mesher: Delaunay Triangulation plus Centroid Voronoi Tessellation
The quality of unstructured mesh is crucial for the numerical simulation.We generate an adaptive mesh by using the Delaunay Triangulation (DT) plus Centroidal Voronoi Tessellation (CVT) (Du et al., 1999;Du and Gunzburger, 2002;Du et al., 2006) algorithm.In computational mathematics, the DT is a special triangulation for a given set of points.DT algorithm maximizes the minimum angle of all the triangle angles.This property can be used to avoid sliver triangles such that we choose DT for the generation of an initial mesh.The mesher honors the rugged topography by taking it as constraint boundaries.Then, we assign the generators onto the boundaries, with spatial intervals adaptive to the local velocities.Next, we circularly "pave" new Delaunay cells from the boundaries toward the interior with the doubly-linked data structure (Blacker and Stephenson, 1991).Note that a singly-linked list is insufficient for "circularly paving" because one cannot remove a singly-linked list entry in constant time.However, the initial mesh lacks quality control (Du and Gunzburger, 2002), and its cell sizes change dramatically, as shown in Fig. 3a.Therefore, we attempt to perform mesh optimisation by using the Centroidal Voronoi Tessellation (CVT) algorithm.
The Voronoi Tessellation of a set of points is dual to its Delaunay triangulation.CVTs are a special kind of special Voronoi Tessellations, whose tessellation generators are coincident with the "mass" centroids of the corresponding Voronoi regions regarding a given "density" function (Du and Gunzburger, 2002).Because the spatial samplings are adaptive to the local migration velocities, in the 2D case, we define the "density" function σ in CVTs as where c max is the velocity maximum, and c 0 (x, z) the migration velocity.Given a Voronoi cell V , its "mass" centroid (x * , z * ) is represented as For Voronoi tessellations {V i } k i=1 , we define the cost (energy) function concerning generators {x i , z i } k i=1 as in which F {x i , z i } k i=1 can be treated as the summed rotational inertia regarding the velocity model.For minimization of Eq. ( 9), we apply an algorithm named Lloyd's method (Du et al., 1999;Du and Gunzburger, 2002), which iteratively constructs the tessellations and centroids of the Voronois.Fig. 2 shows a toy model about how DT-CVT algorithm works.Note that while the interior generators are dynamically moving, we must keep the boundary control points fixed to maintain the topography information.With curved boundary constraints, Fig. 3 shows the comparison between the initial mesh (Fig. 3a) and the optimised mesh (Fig. 3b).With each cell approximating a normal triangle we can see that the latter has better mesh quality over the former.Compared with simpler meshing techniques (for example, only by DT), our technique performs well in controlling the mesh quality.
The unstructured mesh has advantage in depicting irregular interfaces.The unstructured mesh of optimised quality benefits the seismic modeling in terms of accuracy and computational cost.For example, given the initial mesh in Fig. 3a by a simpler mesher, we can see there are some skewed cells inside.In seismic modeling, the time sampling interval is determined by the smallest space interval.According to our experiences in inhomogeneous (also as shown in Fig. 3b), the smallest interval in an optimised mesh is twice larger than that in a simpler unstructured mesh, which helps alleviate the expensive computational cost of the seismic modeling in LSRTM.
The engine of LSRTM is the two-way wave-equation based seismic modeling.The accuracy together with flexibility for seismic modeling is up to its mesh discretization, especially for irregular strong discontinuities.The adaptive unstructured mesh has advantage in accurately depicting discontinuous interfaces.To validate this capability, we compare our unstructured Comparing with (a), the cells in (b) approach more to normal triangles.Note that both meshes have the same number of elements, and the mesher preserves the interfaces well during mesh generation.mesh based modeling method with the finite difference algorithm in the context of wavefield simulation across strong a strong velocity contrast.The model composed of two velocity blocks is shown in Fig. 4a.We first simulate the wavefield propagation using the finite difference method.With a strong displaying scale, we can observe the artificial scattering waves caused by the staircase approximation.The strength of the artifacts depends on the ratio between the central frequency of the source function and the spatial discretization.Usually the lower the ratio is, the more likely the staircase caused artifacts will occur.
As a comparison, we perform the adaptive unstructured mesh based wavefield simulation.Fig. 4c shows the generated mesh by our algorithm.While the interface is delineated well, the mesh size is adaptive to the local velocities.Fig. 4d shows the wavefield running on the generated mesh, which has no scattering artifacts caused by the staircase approximation.Although the amplitude magnitude of these artifacts are not strong, we can avoid them using our method in seismic modeling for the sake of accuracy.
We use a simplified Foothill velocity model to test the capability of our mesher in the context of rugged surface.In Fig. 5, the velocity below the rugged topography ranges from 2500 m/s to 5000 m/s, increasing gradually with depth.There are 39 constraint points in black dots to characterize the rugged topography.We first discretize the topography according to the local velocities, with the control points being fixed.Then, we generate the unstructured mesh using the DT-CVT algorithm.We can see that while conformal to the rugged topography, the grids are adaptive to local velocities.Note that Fig. 5 is only for a quick demonstration.In numerical simulations, e.g. in Fig. 9, our mesh will be more refined and of better quality.

Topography LSRTM workflow
When the mesher and solver are ready, we start to run the LSRTM workflow.Generally, LSRTM consists of three circular sections: (i) demigration (Born modeling), (ii) migration (RTM) and (iii) linearized inversion.Both the demigration and migration sections are based on two-way wave-equation.In AppendixA, we derive in detail the demigration and migration sections.We note that both the first equation in Eq. (A-3) and the equation in Eq. (A-6) can run with Eq. ( 5) directly.The computation of 5 Eq. (A-3) involves the point-source term s (t; x S ).We prepare it with s (t; x S ) M k (x S ) for Eq. ( 5).To run the second equation in Eq. (A-3), we need some small modifications to Eq. ( 5) as follows that this mesh is just for demonstration.In numerical applications, e.g.Fig. 9, we have a refined mesh, which is of better quality.
where m k denotes the local reflectivity model at node k.In Appendix A, we express the Born modeling in a compact matrix form as with L being the demigration operator denoted in Eq. (A-3), m the logarithmic reflectivity model (Operto et al., 2000;Tromp et al., 2005) and d the demigrated data.When m is the true reflectivity, d can be the observed data.In practice, the RTM profile is inferior to the true reflectivity model due to the imperfectness of the adjoint operator.To overcome its imperfectness, we attempt to approach an inverted reflectivity model via a least-squares inversion with the following L2 misfit: The regularization term in Eq. ( 12) is to prevent overfitting when solving an ill-posed problem.The scaling factor, which is used to balance the importance of the regularization versus the data misfit, can be obtained by trial and error.The term W is to penalize roughness in the solution.LSRTM chooses to seek the inverted image via linearized inversion by using the gradient optimisation methods, such as the conjugate gradient algorithm.As stated in Appendix A, we can use the adjoint method to obtain the gradient as Note that we are using the adjoint operator derived on a continuous level.We refer to Xu and Sacchi (2017) to provide a detailed derivation about the exact adjoint operator of the discretised system in the prestack case.Because unstructured mesh based seismic modeling involves an irregular stencil, which is not easy to seek the exact adjoint operator as that in the finite Solver (1) denotes that when the input control paramter is "1", the executable runs as the forward operator.Solver (2) denotes that when the input control paramter is "2", the executable runs as the adjoint operator.Note that in the CG optimisation part, we run "1" in the Solver when we need to determine the line-search step length using the conjugate gradient algorithm.
difference stencil, we leave it for our future research.It is strongly recommended to use the exact adjoint operator.By doing so, the pair of migration/demigration operators will be symmetric, which is a prerequisite for conjugate gradient method to have a better convergence behavior in LSRTM Xu and Sacchi (2017).
We design the topography LSRTM workflow following the "SPECFEM" Komatitsch and Tromp (1999) style.Fig. 6 shows the flow chart.Because the background velocity stays the same throughout LSRTM, we run the M esher once for all.For the Solver part, we give it two control options: "1" for the forward operator and "2" for the adjoint operator.To invert for the reflectivity model, we employ the conjugate gradient algorithm as the optimisation method, in which we need one forward operator and one adjoint operator to seek the line search step length (Dai and Schuster, 2013).The computational cost of iterative LSRTM is expensive, so it is necessary to run with MPI in parallel.The work is distributed as embarrassingly parallel tasks, with one core per task.The CG optimisation script acts as a coordinator to communicate with the Solver to distribute the forward/adjoint operators over different nodes.

Numerical example
We test our methods on the Foothill model (Gray and Marfurt, 1995)

Conclusions
We present an unstructured mesh-based workflow which can handle the LSRTM from rugged topography.Our approach is free of static correction.It mainly consists of three principal parts: mesher, solver, and LSRTM.We generate the unstructured mesh with the DT-CVT algorithm, with the rugged topography serving as boundary constraints.Our mesher honors details on the topography without approximations such as flattening and so on.Also, provided a modeling method, a mesh of adaptive cells allows a lower computational cost compared with a mesh of fixed cells.The solver to simulate the wavefield propagation is an unstructured mesh-based modeling scheme named the grid method.The LSRTM is in its conventional style.In this workflow, we do not need to do static correction before imaging and inversion, but directly inject wavefield downwards through the receivers on the rugged topography.A test on the Foothill model validate our scheme.In comparison with the conventional RTM from rugged topography, the LSRTM from rugged topography promises an imaging profile with more balanced amplitude

Figure 1 .
Figure 1.Local mesh of the grid method.The basic cells are in solid lines.The dashed lines link the centroids of the cells and the midpoint of the edges.The pressure is defined on the nodes in black dots, and the velocity and density at the center of cells in empty circles.

Figure 2 .
Figure 2. A toy model about how DT-CVT works.The Delaunay Tessellations are in solid lines and the Voronoi Tessellations in dashed lines.After the mesh optimisation by DT-CVT, we can see that the blue node moves to the position of the brown node.The optimised mesh on the right-hand-side has better mesh quality.That is, each cell approaches more to a regular triangle that those on the left-hand-side.

Figure 3 .
Figure 3.The comparison between (a) the original mesh by Delaunay Triangulation and (b) the optimised mesh by DT-CVT algorithm.

6Figure 4 .
Figure 4. Comparison between the unstructured mesh and finite difference based wavefield simulations.(a) Two block velocity model with a tilt interface.(b) Finite-difference based wavefield simulation.We have a strong displaying scale for it to intentionally show the staircase caused scattering artifacts.(c) Adaptive unstructured mesh generated in our method.The interface is kept well without staircase approximation.(d) Wavefield simulation on the generated mesh.No scattering artifacts occur around the tilt interface.

Figure 5 .
Figure 5. Unstructured mesh generated by DT-CVT algorithm.The 31 dark spots on the surface control the topography character.As the depth increases, the velocity model becomes faster.While conformal to the topography, the cells are adaptive to the local velocities.Note

Figure 6 .
Figure 6.Topography LSRTM workflow.Because the background velocity in LSRTM stays the same, we run the M esher once for all.
in Fig. 7, which has the rugged topography with a rapid variation of elevations.The model contains some steep faults and dips, typical of mountainous onverthrust regions.The total relief of the rugged topography is approximately 1.4 km.To avoid inverse crime, we use the true velocity model rather than the reflectivity model to generate the observed data.Because the free-surface multiples take no part in the conventional LSRTM, we impose no free surface onto the rugged topography.We have 100 shots spacing at 200 m.The source function is a 25 Hz Solid Earth Discuss., https://doi.org/10.5194/se-2019-37Manuscript under review for journal Solid Earth Discussion started: 11 March 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 7 .
Figure 7.The Foothill velocity model with rugged topography.

Figure 8 .
Figure 8. Rugged topography of the Foothill model.

Figure 9 .Figure 10 .
Figure 9. Local meshes of (a) near-surface and (b) deep parts generated by DT-CVT for the Foothill velocity model.We can see that the mesh in (a) honors the near-surface details (as shown by the black topography control points).Also, from a panoramic view, we can see that the cells in (b) are sparser than (a), meaning that the grids in (b) are larger than (a) thanks to adaptive meshing.

10Figure 13 .
Figure 13.Zoom-in comparison views for detailed investigations.(a) Part of true reflectivity.(b) Part of the topography RTM image.(c) Part of the topography LSRTM image.We see that in reference to (a), both (b) and (c) have the reflectors in the right places, but the image in (c) looks sharper with more balanced amplitude distributions.