Abstract
Least-squares reverse time migration (LSRTM) is powerful for imaging complex geological structures. Most researches are based on Born modeling operator with the assumption of small perturbation. However, studies have shown that LSRTM based on Kirchhoff approximation performs better; in particular, it generates a more explicit reflected subsurface and fits large offset data well. Moreover, minimizing the difference between predicted and observed data in a least-squares sense leads to an average solution with relatively low quality. This study applies L1-norm regularization to LSRTM (L1-LSRTM) based on Kirchhoff approximation to compensate for the shortcomings of conventional LSRTM, which obtains a better reflectivity image and gets the residual and resolution in balance. Several numerical examples demonstrate that our method can effectively mitigate the deficiencies of conventional LSRTM and provide a higher resolution image profile.
Introduction
Seismic migration is an inverse procedure of forward modeling, which can restore the interior of the earth medium with record data. Specifically, migration attempts to eliminate the effects caused by the process of physical propagation and obtain an image that clearly depicts the structural information of interest. Reverse time migration (RTM), a state-of-the-art seismic imaging method (; ), identifies the aforementioned acausal procedure appropriately. Based on two-way wave equation, RTM is powerful for handling complex geological settings and velocity with dramatic variation in the lateral direction. Therefore, it can deal with steep dips and salt dome better than conventional migration (; ; ). However, most migration methods, including RTM, use the adjoint operator to compute the image instead of the inverse operator (). Practical data suffers from many factors, such as irregular acquisition geometry and limited aperture of the acquisition system. These deficiencies generate artifacts and degrade the resolution. To overcome these limitations, least-squares migration (LSM) was proposed to combine with RTM (; ). Therefore, seismic imaging can be regarded as a linearized inverse problem. With a proper initial velocity model, seismic records can be inverted to a more accurate profile. LSRTM iteratively reduces the residual between predicted data and observed data in a least-squares framework; therefore, the adjoint operator can keep approaching the inverse operator. Many results have indicated that LSRTM has a better performance than conventional RTM and migration (; ; ).
The precondition of seismic inversion is forward modeling, which maps the parameter model to seismic data. There are two main approaches to build linear approximation between physical model and wavefield (). One is the most commonly used Born approximation based on small perturbation (; ). This requires that high-order scattered wavefields are much weaker than primary field. The Born operator describes a linear relationship between model perturbation and primary reflected wave. It divides the wavefield into two parts: background wavefield and perturbation wavefield. LSRTM based on Born approximation can achieve model perturbation with these two fields. In addition, an alternative scheme for modeling is Kirchhoff approximation (). Compared with Born modeling, the Kirchhoff operator delineates the connection between primary reflected wave and reflectivity. Different operators lead to distinctive results under these two physical contexts. However, neither Born nor Kirchhoff approximation can avoid the impact on seismic image in a least-squares sense. Because minimizing the L2 norm only provides an average solution (; ). It is essential to seek a balance between the residual and resolution. According to geological recognition, the earth medium usually presents a layered spatial distribution. The reflection coefficient that mirrors strata attributes should be sparsity, that is, the part of model that does not generate reflected wave ought to be zero. Therefore, the inverted model needs a sparse limitation.
This study implements a Kirchhoff modeling formula for LSRTM promoted by sparsity. The reflectivity model should be regularized with L1 norm while minimizing the residual of wavefield in the form of the L2 norm. Referring to ‘least absolute shrinkage and selection operator’ (Lasso) problem (), this reformed LSRTM can be solved by the algorithm of spectral projected gradient for L1 minimization (SPGL1), which is designed to solve sparse least squares (). Examples show that our method can effectively overcome the problems mentioned above.
Method
RTM has great advantages in imaging steep strctures such as salt dome. However, it suffers from low-frequency noise compared to conventional migration. Least-squares migration can get closer iteratively to the optimal solution and eventually obtain a relatively high signal-to-noise ratio, high resolution and amplitude equalized profile that eliminates the influence of the acquisition system. It contains three steps: constructing a linear modeling problem first, using the forward and backward propagation wavefields to image, and finally updating the physics model according to the residual.
Linear Modeling Operators
The linearization of nonlinear forward problem is essential to seismic inversion, making the physical progress more explicit; moreover, converting the medium parameter becomes easier. The choice of a linear operator will lead to different physical significance and images. It is a common way to use Born approximation to realize linearization. The real velocity model is divided into two parts: background velocity and velocity perturbation . Given a perturbation , it generates a corresponding wavefield perturbation The Born operator describes the relationship between reflected wave and model perturbation. Specifically, the incident wave interacting with model perturbation becomes a new source, namely the Huygens principle, and then the new source generates wavefield perturbations. This can be expressed as follows in time domain:where represents the background field propagating in , is the source signature located at and excited at , the model perturbation is denoted by , which describes velocity changes compared to background velocity. is a point in model. This study assumes that the density is a constant (Eq. 1) and (Eq. 2) can be rewritten in form of an integral using Green’s theorem:where is the Green’s function from to , propagates from to . Green’s function is governed by:where the is Dirac function.
Born approximation represents scattered phenomenon caused by model perturbation, which could be a means of linearizing seismic inversion. However, this approximation is accurate when scattered field is much weaker than background field (), which is a disadvantage of Born approximation. It cannot describe kinematic and dynamic information of seismic waves well with strong reflector. And studies have shown that Born approximation has limited angle validity and it cannot appropriately predict the reflections generated with large incident angle ().
Compared to the Born operator, the Kirchhoff operator relates the reflectivity to wavefield perturbation. Therefore, it depicts the interaction between the incident field and reflectivity rather than velocity perturbation. There is a relationship between reflectivity and model perturbation when the perturbation and incident angle are small ():where the is the reflection coefficient at point with incident angle between the incidence and the normal line (Figure 1). This means that we can obtain the wavefield perturbation under the Kirchhoff approximation by substituting (Eq. 6) into (Eq. 4), and we haveHere we turn Kirchhoff modeling equation into the same form as Born approximation. Then Eq. 7 can be rewritten as.
FIGURE 1
It should be noted that the term can be replaced by the generalized angle-dependent reflectivity model to get rid of the limitations of small perturbation and incident angle. Although there are some methods to solve the propagation direction of wave, such as Poynting vector and Plane Wave Decomposition (PWD), it is still tedious and time-consuming to obtain the angle term. Here we give an approximate scheme ().
Each shot can invert a reflectivity image, here we sum the images obtained by all shots. Then, we regard the summation as the final reflectivity model and use it to iterate. Approximately, we can get an averaged reflectivity model by multiple shots stacking. Therefore, we can get the predicted data by using this stacked reflectivity rather than the angle-dependent term . Note that is an averaged reflectivity over all illuminated angles.With this approximate reflectivity , we can express Eq. 9 asIn sum, with the relationship of reflectivity and model perturbation, two linear approximations have a similar form, which expresses their common ground. The difference between two approximations is also evident. From Eq. 6, approximately equals to one and can be ignored for a small incident angle. Therefore, reflectivity can be regarded as the spatial derivative of model perturbation. The inverted model after spatial derivation has a higher resolution, that is, the spectrum has been improved. More details are provided in the numerical tests.
Least-Squares With Sparse Optimization
In contrast to full waveform inversion (FWI) (), LSRTM first establishes a linear relationship between physical model and corresponding response (), then it implements the inverse problem. The least-squares method (LSM) only requires the construction of a migration operator and inverse migration operator, which is conjugated to each other. It can reduce the residual between the observed and predicted data iteratively to approach the optimal solution of the inverse problem gradually. According to the linear approximation above, we can express them in the form of a matrix:where the is predicted data, such as background or perturbation fields. represents modeling operator and is the physics model. Usually, it is assumed that the background velocity has been obtained in advance, and then the data can be predicted. Hence, the misfit function can be expressed as:The model , which makes (the Jacobian matrix) equal to 0, is the optimal solution of Eq. 13. However, the computation of the Jacobian matrix is quite time consuming, particularly for seismic exploration. We adopt the adjoint-state method to calculate the adjoint operator of modeling operator , Specifically, the gradient of can be obtained by back propagation of the wavefield residual and background field, here we give the gradient based on Kirchhoff approximation (; ):Where is the adjoint wavefield governed by:
According to Eq. 13, we can obtain a least-squares solution . Note that LSM provides a smooth solution of the model, which is determined by the properties of the L2 norm. As a result, LSM has a limited ability to improve the quality of the image. Here we give a simple model to display the impact of LSM. In this example, we use the Ricker wavelet with a center frequency of 30 Hz and a time sampling interval of 1 ms. With convolution model theory, we can get seismic records via the convolution of Ricker wavelet and reflection coefficients, which can be obtained by . Conversely, reflection coefficients can be obtained by the deconvolution of seismic records and wavelet, that is, . Figure 2C is the result of deconvolution, and it is hard to identify the reflectors. Compared to deconvolution, Figure 2D shows that LSM improves the resolution obviously. However, many oscillations caused by near the real reflection coefficients should not exist. That’s why we regard the least-squares solution as a smooth or average solution. The actual model indicates that the medium presents a layered spatial distribution, as shown in Figure 2A or Figure 2B, that is, the sub-surfaces are sparse. In Figure 2E, the inversion result with sparse constrained LSM performs quite well, and these oscillations generated by LSM are suppressed; thus, the resolution and sparsity of the reflection coefficient series are improved effectively.
FIGURE 2
Due to the feasibility and sparse property of L1-norm, we modify the objective function with L1 norm to realize sparse reconstruction of the model in this study. Generally,
Eq. 13can be reformed with two new problems.
1 Basis Pursuit (BP) problem
depicts a BP problem that comes from compressed sense theory, and it aims to seek a sparse solution that satisfies
. However, practical seismic data inevitably contain noise, and
Eq. 16can be modified as a basis pursuit denoising (BPDN) problem:
where the
describes the noise level in the data, and
Eq. 16and
Eq. 17are equivalent to each other when
.
2 Least Absolute Shrinkage and Selection Operator (LASSO) problem
where the
is an explicit limitation of sparsity on
. Problems
and
are different descriptions of the same question. They are equivalent in the sense that there exists a solution
of
for a given
, and there exists a corresponding
that makes
also be a solution of
.
Both problems mentioned above can be solved by the algorithm of spectral projected gradient for L1 minimization (SPGL1). Given a constraint , we can obtain the residual norm from Eq. 18:LetEq. 20 recasts as a problem of finding the root of a nonlinear equation and defines a continuous curve, the Pareto curve (Figure 3).
FIGURE 3
For a given , SPGL1 uses the Newton method to approach the root, and as the updates iteratively, the optimal solution of problems and can be obtained. Therefore, we balance the 2-norm of the residual against the 1-norm of the solution eventually (). From Figure 3, the question is degraded to a simple Lasso problem when the noise level factor is equal to 0. In this study, we set . Note that synthetic seismic records do not contain noise in general, so we set the noise level factor to be zero. Actually, the algorithm of SPGL1 can deal with noisy data, and we can add some random noise or set some traces to be zero in synthetic data. Besides, the determination of parameter is quite important. According to the theoretical model, we can calculate the perturbation model or reflectivity model and make a rough estimate of . In general, it is appropriate to set the value of tau to tens of times that of the calculated perturbation model or reflectivity model. Then, the parameter can be adjusted according to the inversion results.
Here we summarize the workflow of L1-regularized LSRTM as follows:
1) Obtain the predicted data with migration velocity , and get the with ;
2) Set the initial model and predict the data based on Born or Kirchhoff approximation, therefore we can get the residual and gradient operator ;
3) Input the parameters of , and set ;
4) Solve the Lasso problem with the algorithm of SPGL1, and update the , and until ;
5) Output the result .
Numerical Example
In this study, two theoretical models are used to test the validity of the proposed method, including a single diffraction point and complex fault model. Both are based on the two-way acoustic wave equation. Here we use the finite difference method on regular grid.
Single-Diffraction Point
To verify the effectiveness of this method, we first set a simple model with a diffraction point of 2000 m/s embedded in the background velocity of 1,000 m/s (Figure 4), and the entire model has been discretized into 201 × 201 grids in the horizontal and vertical directions, respectively, with the same interval of 5m. The geometry system is arranged as follows: a total of 21 shots are uniformly distributed on the surface of this model with an interval of 50m. Geophones are also placed on each grid point on the surface. We use the Ricker wavelet with a center frequency of 25 Hz for modeling, and the sampling interval is 0.5 ms. In this example, we set 1,000 m/s as the migration velocity.
FIGURE 4
As shown in Figure 5A, the image is obtained by LSRTM based on Born approximation. This is consistent with the actual situation to a certain degree. The single scatter point is blurred with a disturbing cross pattern (marked by a yellow arrow). However, it should be a dot on the image (). This is because we use the adjoint operator to migrate rather than the inverse operator in Eq. 14. Specifically, Eq. 13 defines a normal equation with . The term , Hessian matrix, is equivalent to a blur operator acting on the true image . Furthermore, includes the influence of irregular acquisition, limited acquisition aperture, band limited source, etc., which generate artifacts and degrade the resolution of the image (). Compared to LSRTM, the same method in Figure 5C with the L1 constraint performs better; it mitigates the distortion caused by the blur operator. Therefore, with the promotion of sparsity, the resolution in the least-squares method has been improved significantly, and the image looks more like a scatter point.
FIGURE 5
The LSRTM based on the Kirchhoff operator inverts the reflectivity directly from the seismic records. Figure 5B displays the image produced by Kirchhoff approximation. Compared to Figure 5A, least-squares RTM based on Kirchhoff operator suffers from the same problems. Similarly, we implement the L1 norm on LSRTM, which is shown in Figure 5D. The cross pattern is eliminated clearly, and we obtain an explicit dot rather than a blurred spot. Therefore, for a simple model, the sparsity-promoting LSRTM based on Kirchhoff approximation can effectively improve the resolution of the image. The results calculated by LSRTM in Figures 5A,B iterate 5 times. Figures 5C,D use the SPGL1 algorithm for iterating 10 times with .
Fault Model
We also test the other relative complex model. In this fault model (Figure 6A), there are some classical geological structures, including folds, fault blocks, and depressions. Therefore, it appropriately shows the complex structure of near-surface media. The maximum and minimum velocities are 4,000 m/s and 1,500 m/s, respectively. Similarly, we discretize it into 265 × 367 grids with an interval of 5 m. Thus, a total of 25 shots are uniformly located at the surface of this model. The modeling seismic wavelet is same to last experiment.
FIGURE 6
As shown in Figures 7A,B, images inverted by LSRTM based on two approximations fit the fault model well, and the contact relationship between structures can be clearly depicted. To further improve the resolution of these images, we combined L1 norm regularization with LSRTM to reconstruct the model. From Figures 7C,D, the method based on Kirchhoff approximation recovers the stratum’s sparsity more effectively. The results in Figures 7A,B are calculated by LSRTM with 10 iterations. Figures 7C,D use the SPGL1 algorithm for iterating 10 times with and , respectively. Note that the amplitude of inverted results is different because of the value of parameter .
FIGURE 7
Furthermore, we enlarge the model framed by red rectangle in Figure 6A, which has step-like strata (marked by yellow arrows in Figure 8A). After inverting, the reflectivity image in Figure 8D produced by constrained LSRTM based on Kirchhoff operator agrees with the actual situation.
FIGURE 8
Furthermore, images inverted by two different approximations have different phases. According to Eq. 6, reflectivity can be derived from model perturbation. With the assumption of a small incident angle, is roughly equivalent to 1 and can be ignored. Then, Eq. 6 can be rewritten as , where the wavenumber . Therefore, reflectivity is the spatial derivative of model perturbation . As a result, the image inverted by Kirchhoff performs sparser and sharper, and there is a phase shift of 90° between perturbation model and reflectivity model.
Figure 9 shows the amplitude spectra of the images in Figures 7B–D, respectively. The spectra curves are the sum of each trace by the spatial Fourier transform along the depth. The red one is generated from the image inverted by L1-regularized LSRTM based on Born approximation. The blue and green spectrum curves are produced by unconstrained and constrained LSRTM with Kirchhoff approximation, respectively. Because of the spatial derivative and sparse constraint, the spectrum of the image inverted by L1-LSRTM with Kirchhoff approximation has more high-wavenumber components than that of Born approximation, which explains that Kirchhoff approximation improves the resolution of the image.
FIGURE 9
Conclusion
The LSRTM recasts classical seismic inversion as a linear inverse problem. By means of linear approximation, physical model is related to the corresponding wavefield. Thereafter, we can reduce the residual between predicted and observed data iteratively to directly invert the interest parameters. This study introduces two linearization methods. Born approximation obtains the relationship between the model and physical response based on perturbation theory. With the help of the Born operator, we derive another type of linear method, namely the Kirchhoff operator, which relates the reflectivity to wavefield explicitly. Moreover, these two methods have a relationship of a spatial derivative, and there is a phase shift between perturbation model and reflectivity model. Although two operators are different physical quantities, the resolution can be improved by Kirchhoff approximation.
LSRTM can mitigate the shortcomings of other migration methods, while the solution is smooth and deviates from the true model. Specifically, there are redundant oscillatory axes in the strata that should be sparsely distributed. Therefore, we reform the question as a sparsity-promoting LSRTM. The SPGL1 algorithm can effectively solve this problem and invert a sparse image that matches the model well. Examples prove the validity of our study.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
XH-Q: Methodology, Investigation, Formal analysis, Visualization, Writing—Original Draft WX-Y: Conceptualization, Writing—Review and Editing, Project administration WC-Y: Validation, Data Curation, Writing—Review and Editing ZJ-J: Supervision, Writing—Review and Editing, Funding acquisition.
Acknowledgments
We thank the National Natural Science Fund of China (under grant 42074158) for supporting this work. And we are grateful for the comments provided by reviewers, which improve the quality of this paper a lot.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
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
BaysalE.KosloffD. D.SherwoodJ. W. C. (1983). Reverse Time Migration. Geophysics48, 1514–1524. 10.1190/1.1441434
2
BeylkinG. (1985). Imaging of Discontinuities in the Inverse Scattering Problem by Inversion of a Causal Generalized Radon Transform. J. Math. Phys.26, 99–108. 10.1063/1.526755
3
BleisteinN.CohenJ. K.HaginF. G. (1987). Two and One‐half Dimensional Born Inversion with an Arbitrary Reference. Geophysics52, 26–36. 10.1190/1.1442238
4
BleisteinN. (1987). On the Imaging of Reflectors in the Earth. Geophysics52, 931–942. 10.1190/1.1442363
5
DaiW.FowlerP.SchusterG. T. (2012). Multi-source Least-Squares Reverse Time Migration. Geophys. Prospecting60, 681–695. 10.1111/j.1365-2478.2012.01092.x
6
DuttaG.SchusterG. T. (2014). Attenuation Compensation for Least-Squares Reverse Time Migration Using the Viscoacoustic-Wave Equation. Geophysics79, S251–S262. 10.1190/geo2013-0414.1
7
JiangB.ZhangJ. (2019). Least-squares Migration with a Blockwise Hessian Matrix: A Prestack Time-Migration Approach. Geophysics84, R625–R640. 10.1190/geo2018-0533.1
8
LecomteI. (2008). Resolution and Illumination Analyses in PSDM: A ray-based Approach. The Leading Edge27, 650–663. 10.1190/1.2919584
9
LiuH. W.LiB.LiuH.TongX. ‐L.LiuQ. (2010). The Algorithm of High Order Finite Difference Pre-stack Reverse Time Migration and GPU Implementation. Chin. J. Geophys.53, 1725–1733. 10.1002/cjg2.1530(in Chinese)
10
LiuY.ChangX.JinD.HeR.SunH.ZhengY. (2011). Reverse Time Migration of Multiples for Subsalt Imaging. Geophysics76, WB209–WB216. 10.1190/geo2010-0312.1
11
LiuY.HeB.ZhengY. (2020). Controlled-order Multiple Waveform Inversion. Geophysics85, R243–R250. 10.1190/geo2019-0658.1
12
LiuY.LiuX.OsenA.ShaoY.HuH.ZhengY. (2016). Least-squares Reverse Time Migration Using Controlled-Order Multiple Reflections. Geophysics81, S347–S357. 10.1190/geo2015-0479.1
13
McMechanG. A. (1983). Migration by Extrapolation of Time-dependent Boundary Values*. Geophys. Prospect31, 413–420. 10.1111/j.1365-2478.1983.tb01060.x
14
PlessixR.-E. (2006). A Review of the Adjoint-State Method for Computing the Gradient of a Functional with Geophysical Applications. Geophys. J. Int.167, 495–503. 10.1111/j.1365-246X.2006.02978.x
15
SchusterG. T. (2017). Seismic Inversion. Tulsa: Society of Exploration Geophysicists. 10.1190/1.9781560803423
16
StoltR. H.WegleinA. B. (2012). Seismic Imaging and Inversion: Application of Linear Inverse Theory. New York: Cambridge University Press, 1–404. 10.1017/CBO9781139056250
17
TarantolaA. (1984). Inversion of Seismic Reflection Data in the Acoustic Approximation. Geophysics49, 1259–1266. 10.1190/1.1441754
18
TibshiraniR. (1996). Regression Shrinkage and Selection via the Lasso. J. R. Stat. Soc. Ser. B (Methodological)58, 267–288. 10.1111/j.2517-6161.1996.tb02080.x
19
van den BergE.FriedlanderM. P. (2009). Probing the Pareto Frontier for Basis Pursuit Solutions. SIAM J. Sci. Comput.31, 890–912. 10.1137/080714488
20
van den BergE.FriedlanderM. P. (2011). Sparse Optimization with Least-Squares Constraints. SIAM J. Optim.21, 1201–1229. 10.1137/100785028
21
WangX. Y.ZhangJ. J.XuH. Q.TianB. Q. (2021). Least-squares Reverse Time Migration with Wavefield Decomposition Based on the Poynting Vector. Chin. J. Geophys. (in Chinese)64, 645–655. 10.6038/cjg2021O0120
22
WangY. H. (2016). Seismic Inversion: Theory and Applications. Oxford: Wiley, 85–86. 10.1002/9781119258032.ch7
23
WuD.YaoG.CaoJ.WangY. (2016). Least-squares RTM with L1 Norm Regularisation. J. Geophys. Eng.13, 666–673. 10.1088/1742-2132/13/5/666
24
YangK.ZhangJ. (2019). Comparison between Born and Kirchhoff Operators for Least-Squares Reverse Time Migration and the Constraint of the Propagation of the Background Wavefield. Geophysics84, R725–R739. 10.1190/geo2018-0438.1
25
YoonK.ShinC.SuhS.LinesL. R.HongS. (2003). 3D Reverse-Time Migration Using the Acoustic Wave Equation: An Experience with the SEG/EAGE Data Set. The Leading Edge22, 38–41. 10.1190/1.1542754
26
ZhangY.DuanL.XieY. (2015). A Stable and Practical Implementation of Least-Squares Reverse Time Migration. Geophysics80, V23–V31. 10.1190/geo2013-0461.1
27
ZhuJ.LinesL. R. (1998). Comparison of Kirchhoff and Reverse‐time Migration Methods with Applications to Prestack Depth Imaging of Complex Structures. Geophysics63, 1166–1176. 10.1190/1.1444416
Summary
Keywords
least-squares reverse time migration (LSRTM), kirchhoff approximation, L1-norm regularization, sparsity constraint, born approximation
Citation
Hong-Qiao X, Xiao-Yi W, Chen-Yuan W and Jiang-Jie Z (2021) Sparse Constrained Least-Squares Reverse Time Migration Based on Kirchhoff Approximation. Front. Earth Sci. 9:731697. doi: 10.3389/feart.2021.731697
Received
28 June 2021
Accepted
17 August 2021
Published
03 September 2021
Volume
9 - 2021
Edited by
Hao Hu, University of Houston, United States
Reviewed by
Jincheng Xu, Southern University of Science and Technology, China
Chuang Li, Xi’an Jiaotong University, China
Updates
Copyright
© 2021 Hong-Qiao, Xiao-Yi, Chen-Yuan and Jiang-Jie.
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: Zhang Jiang-Jie, zhangjj@mail.iggcas.ac.cn
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.