
Citation: | Wei Ouyang, Weijian Mao, Xuelei Li, Wuqun Li (2014). Seismic inversion with generalized Radon transform based on local second-order approximation of scattered field in acoustic media. Earthq Sci 27(4): 433-439. DOI: 10.1007/s11589-014-0092-x |
Sound velocity inversion problem based on scattering theory is formulated in terms of a nonlinear integral equation associated with scattered field. Because of its nonlinearity, in practice, linearization algorisms (Born/single scattering approximation) are widely used to obtain an approximate inversion solution. However, the linearized strategy is not congruent with seismic wave propagation mechanics in strong perturbation (heterogeneous) medium. In order to partially dispense with the weak perturbation assumption of the Born approximation, we present a new approach from the following two steps: firstly, to handle the forward scattering by taking into account the second-order Born approximation, which is related to generalized Radon transform (GRT) about quadratic scattering potential; then to derive a nonlinear quadratic inversion formula by resorting to inverse GRT. In our formulation, there is a significant quadratic term regarding scattering potential, and it can provide an amplitude correction for inversion results beyond standard linear inversion. The numerical experiments demonstrate that the linear single scattering inversion is only good in amplitude for relative velocity perturbation (δc/c0) of background media up to 10 %, and its inversion errors are unacceptable for the perturbation beyond 10 %. In contrast, the quadratic inversion can give more accurate amplitude-preserved recovery for the perturbation up to 40 %. Our inversion scheme is able to manage double scattering effects by estimating a transmission factor from an integral over a small area, and therefore, only a small portion of computational time is added to the original linear migration/inversion process.
In exploration geophysics, the core task is to identify oil and gas resources in the subsurface. Mechanical properties inversion is a principal objective of this task, which provides not only geometrical imaging of a medium but also recoverable estimation of its physical parameters. Rigorously, direct inversion method based on wave equation which controls seismic wave propagation in actual heterogeneous medium, until now, is less understood, because the inversion is exactly nonlinear. In practice, simplified inversion algorithms are widely probed in order to obtain an approximate solution of properties of the actual medium. The most popular of them is the linearization of inverse problem (Born/single scattering forward approximation) and solving the linearized inverse problem (e.g., Cohen and Bleistein 1979; Raz 1981; Clayton and Solt 1981; Bleistein and Gray 1985; Blesitein et al. Bleistein et al. 1985; Symes 2008; Weglein et al. 2009).
In the real world of the seismic inversion, the main obstacle is that recorded data consist of imperfect frequency and amplitude information which will result in parameters inversion distortion. From this point, one does not concern the complete restoration of the physical parameters, but how much they jump (their gradient) at their discontinuous locations. Beylkin (1985) paper firstly embodied this inversion strategy. Beylkin made use of powerful mathematical tools, such as pseudo differential operators and generalized Radon transforms (GRT), and reduced the inverse problem to the GRT inversion, which relates to Beylkin determinant (BD). This inversion procedure assumed that the scattered field is approximated by Born approximation, and BD is nonsingular (i.e., finite, nonzero, and single ray path). Bleistein (1987) proposed a modification of GRT method of Beylkin (1985). The modification formulates the forward problem as Kirchhoff single scattering from interfaces instead of Born single scattering from volumes and results in reconstruction of reflection coefficient. Beylkin and Burridge (1990) further extended the GRT inversion to an elastic isotropic medium and proposed a linear multiparameter inversion scheme based on how to exploit the angle dependence of the single scattering. The extended GRT method does not need to recover the angle information at the specular point which is required by the method of Bleistein (1987). The GRT approach to the linear inversion of seismic data for anisotropic elastic parameters has been investigated by de Hoop et al. (1994), (1999) and Burridge et al. (1998).
Our preference in this paper is to relax the constraint of Born approximation which is only valid for small scattering potential and single scattering propagation pattern. We present an approach to handle nonlinear inverse problems in acoustic medium by taking into account the second-order Born approximation. Direct inversion solution for second-order Born approximation is not easily obtained. We here propose a scheme to approximately deal with second scattering problem. Assuming the main contribution of double scattering comes from a local area around the first scattering point, the scattered field of second-order Born approximation can be also expressed by GRT about quadratic scattering potential. With the GRT formula of the quadratic scattering potential, we invert scattering potential according to Beylkin (1985) method. The numerical experiments demonstrate that our quadratic inversion can give more accurate recovery for amplitude, and is compatible with relative velocity perturbation (δc/c0) up to 40 %, compared with Born linear inversion which is only good in the perturbation up to 10 %. We also find that the choice of radius of double scattering domain is closely related to perturbation intensity of velocity, and the radius decreases as the velocity perturbation increases, which is preliminarily verified by numerical tests based on quadratic inversion for second-order Born seismic data at the end of the section.
In a constant density medium the acoustic wave equation can be expressed as
|
(1) |
where U is the pressure field, ω the angular frequency, x a point in the imaging domain, s a fixed point source, c the variable velocity, and Δx the Laplacian operator at the point x. If the scattering potential is defined by
|
(2) |
where c0(x) is the background velocity (c=c0+δc), U can be written as
|
(3) |
where y is located in the imaging domain or on its boundary, U0 the Green function for background velocity c0, and Usc the scattering field. With these definitions, analogous to the Lippman-Schwinger equation, we arrive at
|
(4) |
Equation (4) is a nonlinear integral equation for scattering potential. When the scattering potential is small, we can obtain the linearized integral equation by replacing the pressure field U inside the integral by background field U0, i.e.,
|
(5) |
Equation (5) is conventionally called a Born approximation or single scattering approximation, and its approximate solution for f can be obtained by the inversion of a causal GRT (Beylkin 1984, 1985; Miller et al. 1987). We should note that Born approximation is only valid for small scattering perturbation. Until now, how to directly solve scattering potential from integral Eq. (4) is unknown. The most constructive solution may belong to inverse scattering series (Weglein et al. 2003), but this series is extremely difficult to achieve in practical depth migration inversion implementation.
In order to reduce the limitation of the small perturbation and make inversion algorithms as realistic as possible, we consider second-order Born approximation for nonlinear integral Eq. (4) which can be exactly written
|
(6) |
where r represents the receiver position (Fu and Bouchon 2004). In general, there is no direct inversion solution for f from integral Eq. (6) so far. Here, we attempt to use the GRT method to asymptotically invert the scattering potential from Eq. (6) under appropriate assumptions.
For any fixed point x, we only consider the interaction with another point y located within a small area Bx which is surrounding x in the imaging domain. Note that the background velocity c0 is generally obtained by smoothing real velocity model. Within the local volume area Bx, we can appropriately assume f(y)≈f(x), which is crucial for relationship between scattered field and GRT on quadratic scattering potential.
Therefore, Eq. (6) can be recast as
|
(7) |
Let
|
(8) |
where E(x, s, ω) is an integral associated with the background Green function and can be called 'transmission factor'.
Now we compute the term E(x, s, ω) in 2D, the 3D case can be similarly achieved. Given an integral area Bx with center at x and radius εequal to a percentage of the wave length, i.e.,
|
(9) |
|
(10) |
where A(x, y) and τ(x, y) satisfy transport equation and eikonal equation, respectively, (Miller et al. 1987). Substituting for U0 in Eq. (8) with Eqs. (9) and (10) gives
|
(11) |
We also assume
|
(12) |
within the local integral area Bx. Then the transmission factor becomes
|
(13) |
which is obtained by polar coordinates transform, and the fact that
Using Eqs. (9) and (13), Eq. (7) can be written as
|
(14) |
where AT(s, x, r) and τT(s, x, r) are the total amplitude and traveltime function, respectively, i.e.,
|
(15) |
Integral Eq. (14) can be related to a casual GRT about 'quadratic scattering potential' defined by
|
(16) |
To see this, consider the transform
|
(17) |
Therefore, we have
|
(18) |
where the Fourier transform
So far we have established the relationship between scattered field and causal GRT based on second-order Born approximation and aforementioned assumptions, in the next section we will consider the problem of finding the high-frequency asymptotic solution f(x) in Eq. (18).
Following Beylkin 1985's derivation in 2D, we define the backprojection operator correspondent with the GRT appeared in Eq. (17)
|
(19) |
where u(t, r, s) is a smooth function of time t and receiver r, s is a fixed point source, z is imaging point, and h(z, r, s) is Beylkin determinant (see its exact definition in Beylkin 1985).
If we choose u(t, r, s) as follows
|
(20) |
where Usc is given in Eq. (14), then apply the back projection operator to Eq. (20); we obtain the quadratic equation of scattering potential
|
(21) |
where T is a smooth pseudo differential operator (Hörmander 1970) and 〈·〉 is an operator of partial reconstruction (Beylkin 1985).
Neglecting the smooth term in Eq. (21), we obtain a high-frequency asymptotic reconstruction algorithm of scattering potential in terms of second-order Born approximation
|
(22) |
Since the left side of Eq. (22) is a migrated image, it contains all information about the discontinuities of the perturbations. By solving the quadratic Eq. (22), we can not only ascertain the location of discontinuities of scattering potential f but also obtain the jump of scattering potential at the discontinuous position of the perturbation. Note that if there is no quadratic term in the Eq. (22), we just obtain the approximate solution for the linearized inverse problem which is corresponding to first-order Born approximation Eq. (5).
In order to simply testify the merits of the quadratic inversion Eq. (22) and the limitations of linear inversion, we conduct a set of typical numerical experiments and will compare the normalized inversion results in the aspect of amplitude preservation.
Seven simple 2D models with different relative velocity perturbations (δc/c0, where δc=c-c0) are set up, each with two flat layers and the interface at the depth of 1, 250 m. The velocity in the top layer is 2, 000 m/s, and the velocity in the bottom layer is, respectively, 2040, 2200, 2400, 2600, 2800, 3200, and 3600 m/s. The synthetic data are generated from finite-difference modeling for a fixed shot located at the center of the surface, and 401 receivers (Fig. 1). The source function is chosen to be Ricker wavelet with a dominant frequency of 20 Hz. The receiver spacing is 10 m. We used a constant 2, 000 m/s as background velocity, so the relative perturbations below the interface are exactly 2, 10, 20, 30, 40, 60, and 80 % for corresponding models.
The inverted scattering potential images from first-order and second-order Born approximation (we can choose integral radius ε to be within a range between 20 and 50 % of main wavelength; for simplicity, here let l=0.2, ε=0.4πc0/ω) are shown in Figs. 2 and 3. Since the first-order solution from small perturbation (2 %) model is accurate, the source wavelet amplitude can be obtained from the ratio of the true scattering potential to the inverted one. By using it, we remove the source wavelet effects. Furthermore, we normalize the inverted results by corresponding true one. When the normalized result is close to unitary, it means the inverted solution is accurate. Figure 4 displays the normalized results picked along the interface from 7 models for comparison. The results show that the inverted parameters from GRT linear inversion are different from the true models significantly when the perturbation is greater than 10 %, which is accordant with the validity of first-order Born approximation to the volume-scattering wave studied by Fu and Bouchon (2004). For the inversion from our GRT quadratic inversion, the results are close to the true ones, even if the perturbation is up to 40 %. If the velocity perturbation exceeds 40 %, the wide angle error is unacceptable in both GRT linear and quadratic inversion. However, the results from quadratic inversion are still accurate at near-offsets.
The forward modeling by Born approximation is simply a superposition of the scattered fields from all individual scattering points, and omits the innate nonlinear relationship between the alternation of medium properties and the concomitant change in the scattered field. As pointed by Hudson and Heritage (1981), the applications of Born approximation beyond its validity may be accepted in phase (traveltime), but do not conserve energy (in amplitude). This has been verified for heterogeneous media by Wu (1996) and Fu (2006).
To solve nonlinear inverse scattering problem in acoustic medium, one also generally linearizes integral Eq. (4) and obtains the asymptotical linear solution. The drawbacks of the linearized inversion procedure are obvious; especially it neglects multiple scattering effects and is only valid for small scattering potential. Many relative research works are now concentrated on how to abandon the assumption of Born (single scattering) approximation. The most creative work may originate from inverse scattering series algorithms associated to Born series (Weglein et al. 2003). However, this direct inversion has its innate complexity in dealing with realistic model.
To obtain a direct asymptotical solution for complicate inverse scattering imaging is difficult, as a first attempt, we propose to solve the acoustic inverse problem based on second-order Born approximation. We have focused on double scattering solution and asymptotic calculation of the corresponding transmission factor within a local area. The local second-order approximation of scattered field is naturally related with a causal GRT defined on second-order scattering potential. After standard GRT inversion, we obtain a quadratic inversion formula. The synthetic examples we presented further confirm that our quadratic approximation can deal with double scattering problem from the local area and reconstruct large perturbation model accurately. Compared with the linear inversion, our quadratic inversion formula involves only an additional term associated with transmission factor, which results that only a small portion of computation time is added to original migration inversion processing.
Although Eq. (22) is an approximate inversion solution of integral Eq. (6), it has more actual significance and gives better results than linear inversion method. The parameter l is an important factor in Eq. (22), which is proportionally related to the radius of double scattering domain. In the previous section, we obtain a reasonable solution for numerical tests by roughly choosing l=0.2. Here, we further discuss the behaviors of l in the inversion in order to give a guide of how to set the parameter. In terms of scattering propagation, the reasonable radius (parameter l) should decrease as velocity perturbation increases. We preliminarily conduct a set of numerical experiments to verify this claim. The experiment configurations are the same as one appeared in the last section (i.e., the velocity in the top layer is 2, 000 m/s, and the velocity in the bottom layer is, respectively, 2040, 2200, 2400, and 2800 m/s). We generate the first-order and second-order Born synthetic datasets, and apply the linearized inversion to the two datasets. The results show that the linear inversion solutions for the first-order Born data are perfectly match the true models, but the linear inversion solutions for the second-order Born data are different from the true models significantly (Fig. 5). Then we perform the quadratic inversion for second-order Born data by parameter l=0.2, 0.19, 0.16, 0.135, respectively, the quadratic inversion solutions match the true models very well (Fig. 5). The results from this experiment are consistent with what we expected and further validate the quadratic inversion formula.
The work presented here is a first and important step for true-amplitude nonlinear migration/inversion with GRT. The encouraging numerical results from simple model with two homogenous flat layers guide us to further comprehensively investigate the effectiveness of quadratic nonlinear inversion in complex medium.
This research was supported by Innovation Project of Chinese Academy of Sciences and State Key Laboratory of Marine Geology, Tongji University (No. MGK1408). Authors would like to thank two reviewers for their constructive suggestions.
Beylkin G (1984) The inversion problem and applications of the generalized Radon transform. Commun Pure Appl Math 37:579-599 doi: 10.1002/(ISSN)1097-0312
|
Beylkin G (1985) Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform. J Math Phys 26:99-108 doi: 10.1063/1.526755
|
Beylkin G, Burridge R (1990) Linearized inverse scattering problem of acoustic and elasticity. Wave Motion 12:15-52 doi: 10.1016/0165-2125(90)90017-X
|
Bleistein N (1987) On the imaging of reflectors within the earth. Geophysics 52:931-942 doi: 10.1190/1.1442363
|
Bleistein N, Gray SH (1985) An extension of the Born inversion procedure to depth dependent velocity profiles. Geophys Prospect 33:999-1022 doi: 10.1111/gpr.1985.33.issue-7
|
Bleistein N, Cohen JK, Hagin FG (1985) Computational and asymptotic aspect of velocity inversion. Geophysics 50:1253-1265 doi: 10.1190/1.1441996
|
Burridge R, de Hoop MV, Spencer C (1998) Multiparameter inversion in anisotropic elastic media. Geophys J Int 134:757-777 http://cn.bing.com/academic/profile?id=30b49b71a912d92401317de06c60405f&encoded=0&v=paper_preview&mkt=zh-cn
|
Clayton RW, Solt RH (1981) A Born-WKBJ inversion method for acoustic reflection data. Geophysics 46:1559-1567 doi: 10.1190/1.1441162
|
Cohen JK, Bleistein N (1979) Velocity inversion procedure for acoustic waves. Geophysics 44:1077-1087 doi: 10.1190/1.1440996
|
de Hoop MV, Burridge R, Spencer C, Miller DE (1994) GRT/AVA migration/inversion in anisotropic media. Proc SPIE 2301:15-17 doi: 10.1117/12.187482
|
de Hoop MV, Spencer C, Burridge R (1999) The resolving power of seismic amplitude data: An anisotropic inversion/migration approach. Geophysics 64:852-873 doi: 10.1190/1.1444595
|
Fu LY (2006) Comparison of different one-way propagators for wave forward propagation in heterogeneous crustal wave guides. Bull Seismol Soc Am 96:1091-1113 doi: 10.1785/0120050159
|
Fu LY, Bouchon M (2004) Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media-I. Theory of two-dimensional SH case. Geophys J Int 157:481-498 http://cn.bing.com/academic/profile?id=37b113705d50511460cb708467ceadeb&encoded=0&v=paper_preview&mkt=zh-cn
|
Hörmander L (1970) Fourier integral operators. I. Acta Mathematica 127:79-183 http://cn.bing.com/academic/profile?id=aa512809fe97d41379a3b42f3c004a27&encoded=0&v=paper_preview&mkt=zh-cn
|
Hudson JA, Heritage JR (1981) The use of the Born approximation in seismic scattering problems. Geophys J R Astr Soc 66:221-240 doi: 10.1111/j.1365-246X.1981.tb05954.x
|
Miller DE, Oristaglio M, Beylkin G (1987) A new slant on seismic imaging: migration and integral geometry. Geophysics 52:943-964 doi: 10.1190/1.1442364
|
Raz S (1981) Three-dimensional velocity profile inversion form finite offset scattering data. Geophysics 46:837-842 doi: 10.1190/1.1441221
|
Symes WW (2008) Approximate linearized inversion by optimal scaling of prestack depth migration. Geophysics 73:R23-R35 doi: 10.1190/1.2836323
|
Weglein AB, Araújo FV, Carvalho PM, Stolt RH, Matson KH, Coates RT, Corrigan D, Foster DJ, Shaw SA, Zhang H (2003) Inverse scattering series and seismic exploration. Inverse Probl 19:R27-R83 doi: 10.1088/0266-5611/19/6/R01
|
Weglein AB, Zhang H, Ramírez AC, Liu F, Lira JEM (2009) Clarifying the underlying and fundamental meaning of the approximate linear inversion of seismic data. Geophysics 74:WCD1-WCD13 doi: 10.1190/1.3256286
|
Wu RS (1996) Synthetic seismogram in heterogeneous media by one-return approximation. Pure Appl Geophys 148:155-173 doi: 10.1007/BF00882059
|