Reconciling processing and inversion: multiple attenuation prior to wave-equation inversion

Claudio Guerra and Alejandro Valenciano

ABSTRACT

Seismic inversion is very sensitive to the presence of noise. In an inversion scheme, noise is any event in the data not predicted by the forward modeling because of either inadequate physics or incomplete model parameterization. Therefore, noise exists in the residual space and, if coherent, slows convergence to an acceptable result, or in the worst case, dominates the whole process, inhibiting the efficacy of inversion. The ideal solution is to include the noise in the forward modeling, but this can lead to a prohibitively expensive inversion scheme in many cases. A more practical approach is to make the data agree with the physical assumptions of the inversion scheme. In this sense, noise attenuation is a preprocessing step before inversion. Here, we illustrate the problem by applying one-way wave-equation inversion to a portion of the well known Sigsbee2b data. In the present example, noise is represented by multiples, not modeled by the one-way wave-equation. After characterizing the noise in the migrated data, we use a dip filter to estimate it and a non-stationary adaptive filter technique to subtract it from the migrated data. The results attest to the importance of agreement between the physics of the inversion and the data (or vice-versa).

INTRODUCTION

Migration with the full wave equation is, in principle, capable of creating an image of the subsurface by reversing the propagation of all wave modes (Wapenaar et al., 1987). It is a computationally demanding process, which is why the wave equation is approximated by the one-way wave equation. Migration with the one-way approximation is formulated such that, ideally, the only seismic events present in the input data are primary reflections. When this is not the case, because the pre-processing was unable to completely eliminate noise (multiple reflections, for instance), the final image presents migrated events that do not correspond to the geology. In other words, these events are not modeled by the physical theory on which the migration method is based. They are frequently called migration noise or migration artifacts.

Multiple attenuation prior to inversion  2  Guerra and Valenciano

As is widely known, any inversion problem is very sensitive to noise. In seismic, several authors have shown that coherent noise degrades the performance of the velocity inversion and yields erroneous results (Chauris and Noble, 2001; Shin and Min, 2006; Li and Symes, 2007). This problem is more severe when the physics which governs the the inverse problem does not model the noise. Therefore it is usually necessary to include some kind of a priori information about the model as an additional regularization equation in the inversion or to perform noise suppression prior to inversion.

Noise also affects least-squares inversion of migrated images into reflectivity. In offshore data, multiples are the main coherent noise. (Valenciano, 2008) addresses the problem of migrated multiple reflections in inversion by applying a pre-processing step to attenuate them. Here, using the Sigsbee2b dataset (Paffenholz et al., 2002), we characterize the multiples in the subsurface-offset domain, detail the pre-processing step and show the impact it has on inversion.

LINEAR LEAST-SQUARES INVERSION

Linear least-squares inversion provides a theoretical approach to compensate for experimental deficiencies (e.g., limited acquisition geometry) and complexities of the overburden, while maintaining consistency with the acquired data. For seismic imaging, it can be summarized as follows. Given a linear modeling operator L, compute synthetic data d using d D Lm, where m is a reflectivity model.

The quadratic cost function,

S(m) D ∥d − d_{obs}∥^{2} D ∥Lm − d_{obs}∥^{2},  (1)

is formed, where d_{obs} denotes the recorded data.

The reflectivity model m̂ that minimizes S(m) is given by

m̂ D (L′L)^{−1} L′d_{obs} D H^{−1} m_{mig},  (2)

where L′(the migration operator) is the adjoint of the linear modeling operator L, m_{mig} is the migrated image, m_{mig} D L′d_{obs}, and H D L′L is the Hessian of S(m).

Since the model space can be large, computing the inverse of the Hessian matrix is a big challenge for most geophysical imaging problems. For this reason, it is often more feasible to compute the inverse image as the solution of the linear system of equations,

Hm̂ D m_{mig},  (3)

by using an iterative inversion algorithm. In this approach to the inverse problem, only one matrix-vector multiplication of the Hessian matrix with the model vector is necessary per iteration, in contrast with other methods (Clapp, 2005) that require one migration and one modeling every iteration. Still, explicit computation of the Hessian for the entire model space is too expensive in practice. (Valenciano, 2008) discusses how exploiting the structure of the Hessian matrix and the localization in the model space makes this problem tractable.

By using a priori information about the covariance of the model (model regularization), one can add regularization to solve the otherwise ill-posed inversion problem. A more customary regularization for the inversion in the poststack image domain is to add a damping factor that penalizes an increase of the values of the model. This regularization makes no use of any physical knowledge we might have about the seismic reflectors. It is implemented by adding a small value to the diagonal of H(x,x′) in equation 3:

(H C εI)m̂− m_{mig} D r ≈ 0,  (4)

where I is the identity operator, r is the residuals vector, and ε is an scalar parameter that governs the strength of the regularization.

NOISE CHARACTERIZATION

The occurrence of shadow zones below the salt body is strongly conditioned by the concavity of its base. In the Sigsbee2b model, the comparison of the diagonal of the Hessian matrix in Figure 1 with the reflectivity model in Figure 2 shows that where the base of the salt is concave down, the transmitted energy is focused and high illumination values are produced; where the base of the salt is concave up, the transmitted energy is spread along to divergent paths, producing low illumination zones. Unfortunately, for reflected energy and also for multiple energy, this geometry acts in the opposite way, concentrating energy where the base of the salt is concave up and dispersing energy where it is concave down. In the shadow zones, therefore, it is very likely that the strongest remaining amplitudes correspond to multiple reflections, while in the well-illuminated zones these events produce little or no spurious interference with the primary energy (Figure 3). Consequently, during one-way wave-equation inversion, multiple energy can dominate the residuals, surpassing primaries in the shadow zones.

Figure 1: Illumination pattern given by the diagonal of the Hessian matrix of the Sigsbee2B model. Dark gray represents low illumination and light gray represents high illumination diag-Sis [ER]

Multiple attenuation prior to inversion  4  Guerra and Valenciano

Figure 2: Reflectivity model of Sigsbee2B model. refle-new [ER]

Figure 3: Migration of Sigsbee2b data – zero-subsurface-offset section. mig-stack [ER] SEP–134

One solution to this problem is to attenuate the multiples (events that, in inversion, do not fit the model) from the migration. This requires that the multiples be separated from the primaries or modeled in such a way that during subtraction the primaries are minimally affected.

Traditional data-space demultiple schemes dealing either with periodicity or differential moveout become ineffective in sub-salt settings for two main reasons: a) complex ray paths make periodicity assumptions fail, and b) salt focusing effects concentrates sub-salt reflections in the near offset range, making Radon-type demultiple schemes ineffective. (Sava and Guitton, 2005) and (Alvarez et al., 2007) demonstrated that primaries and multiples have different behavior in image space (subsurface-offset or reflection-angle), and that multiples can be adaptively subtracted from the migration in the image domain without significantly affecting the primaries.

The Sigsbee2B model (Paffenholz et al., 2002) has a high-reflectivity water-bottom. In contrast with Sigsbee2a, this characteristic generates strong peg-leg multiples, which bounce between the salt edge and the water-bottom. In Figure 4 we use part of the zero-offset section of Sigsbee2b data (non-free-surface data) to highlight these multiples and the salt limits in the time domain. The diffractions labeled 1 correspond to the undulating top of the salt; that labeled 2 is the reflection from its base; 1^{st}-order peg-leg diffracted multiples, concurrent with the base of the salt, are labeled 3; the 1^{st}-order peg-leg multiple from the base of the salt is identified with label 4; and internal multiples generated inside the salt body are labeled 5.

Figure 4: Sigsbee2b zero-offset time section. The labels correspond to 1) top of the salt; 2) base of the salt; 3) 1^{st}-order peg-leg diffracted multiples; 4) internal multiples originated at the salt limits; and, 5) 1^{st}-order peg-leg multiple from the base of the salt. ZOtime [ER]

After prestack depth-migration, multiple reflections are distorted, because they are downward propagated beyond their actual reflection depth, using the incorrect velocities of the time-concurrent primaries. If their reflection traveltimes are coincident with sub-salt reflections, their propagation directions are severely deviated; therefore, the shape they assume shows a strong imprint of the velocity model.

In Figure 5 the migration result of Sigsbee2b is shown in a magnified view. The frame on the left is the zero-subsurface-offset section, and the one on the right is a subsurface offset gather – ODCIG – taken at CMP position 38825 ft. The 1^{st}-order peg-leg diffracted multiples, labeled 1, show the effect of over-migration, because they are migrated with salt velocity; labels 2 and 4 stand for migrated peg-leg multiples from the base of the salt, and immediately above them, labeled 3 and 5, are the migrated internal multiples. Some reflectors show focusing close to zero-offset. Notice the strong correlation between the shape of migrated multiples and the base of the salt caused by the deviation of the propagation direction mentioned above. For comparison, Figure 6 shows the migration of Sigsbee2a. In this case, multiples related to the water-bottom (labels 1, 2 and 4) are much weaker or absent, because the water-bottom in this model is characterized by a “soft” interface.

Figure 5: Sigsbee2b shot-profile migration using the cross-correlation imaging condition. The labels correspond to migrated events from 1) 1^{st}-order peg-leg diffracted multiples; 2) 1^{st}- order peg-leg multiple from the base of the salt; 3) internal multiples originated at the salt limits; 4) 1^{st}-order peg-leg multiple from the base of the salt; and, 5) internal multiples originated at the salt limits. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at x D 38600 ft. sig2b01 [ER]

Multiples are more evident in regions of low illumination. In Sigsbee data, these regions are associated with a concave-up base of the salt, in which downgoing energy is defocused and upgoing energy is focused. Much of the remaining energy in these regions corresponds to multiples. Figure 7 exemplifies this problem. The ODCIG on the right frame was selected from a low illumination area. The horizontal to upward-curved events are migrated multiples, Figure 6: Sigsbee2a shot-profile migration using the cross-correlation imaging condition. The labels correspond to migrated events from 1) 1^{st}-order peg-leg diffracted multiples; 2) 1^{st}- order peg-leg multiple from the base of the salt; 3) internal multiples originated at the salt limits; 4) 1^{st}-order peg-leg multiple from the base of the salt; and, 5) internal multiples originated at the salt limits. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at x D 38600 ft. Compare with Figure 5. sig2a01 [ER] which contaminate the entire sub-salt section. Because of the low illumination, some reflectors appear as dipping events below 14000 ft. Notice that, in the zero-subsurface-offset section, because of the imprint of the salt velocity, all the multiple modes show different dips than the reflectors.

The different dip behavior of primaries and multiples can be used as a criterion to separate them using a dip filter in the k_{x} − k_{h} wavenumber domain in a pre-processing step. After discrimination, and prior to subtraction from the original data, amplitude and phase of the estimated multiples must be adjusted. We perform this correction by using non-stationary filters, according to the strategy of (Alvarez and Guitton, 2007). They advocate that the best adjustment between the estimated multiples and the multiples in the data is achieved by simultaneously matching the estimation of primaries and multiples to the data containing both. The strategy uses small overlapping patches of the input data to compute local filters in a least-squares inverse problem.

Figure 7: Sigsbee2b shot-profile migration (subsurface-offset) using the cross-correlation imaging condition. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at x D 36800 ft. Poorly illuminated areas are dominated by the presence of migrated multiples. migoff1 [ER]

The next section presents the results of pre-processing and its impact on the inversion output.

RESULTS

The Sigsbee2B model Hessian matrix (diagonal shown in Figure 1) and the migration (Figure 3) were computed using one-way wave-equation propagators. As result, the multiples were neither included in the modeling nor propagated with the correct velocity/direction during migration. For this reason, the multiples will be treated as noise for the inversion. In the following, we will show how they prevent the inversion fail from recovering the primary signal in the shadow zones.

Figure 7 shows the migration result in the subsurface-offset domain. The front face corresponds to the image at zero subsurface offset, and the side face corresponds to the subsurfaceoffset gather at a position inside the shadow zone (x D 36800 ft). Since the primaries are illuminated at few reflection angles and from a dominant slanted wave-propagation direction, their signature in the subsurface-offset domain is a slanted line. Because of the multiple bounces on dipping interfaces (15^{∘} watter-bottom and 30^{∘} salt flank), multiple energy is spread horizontally in the ODCIG, resembling a zero-reflection-angle event. Figure 8 shows the migration result transformed to the reflection-angle domain. The front face corresponds to the image at the 0^{∘} reflection angle, and the side face corresponds to the reflection-angle gather at a position inside the shadow zone (x D 36800 ft). Notice how multiples are concentrated around the 0^{∘} reflection angle, and primaries are mapped at higher reflection angles. This behavior of the multiples and primaries can also be found in field datasets (Valenciano, 2008).

Figure 8: Sigsbee2b shot-profile migration (reflection angle) using the cross-correlation imaging condition. The front face corresponds to the image at the 0^{∘} reflection angle, and the side face corresponds to the reflection-angle gather at x D 36800 ft. migang1 [ER]

We used a three-stage pre-processing strategy to attenuate the multiple energy before inversion. It includes discriminating the multiple energy in the k_{x}−k_{h} plane (to generate a model of the multiples); performing amplitude and phase correction of the discriminated multiples in a least-square sense; and subtracting them from the original data.

The discrimination of the multiple energy relies on the differences in dip patterns as described previously: multiples show up in the k_{x}−k_{h} plane at low k_{h} and high k_{x}, and primaries at high k_{h} and low k_{x}. Therefore, we build a model of the primaries by submitting the migrated data to a k_{x} − k_{h} filter at every depth step. Figure 9 shows a depth slice at z D 14150 ft of the migration in the k_{x} − k_{h} plane, and Figure 10 shows the result of filtering. The model for the multiples (Figure 11) is obtained after subtraction of the filtered result (Figure 10) from the migration (Figure 9).

Figure 9: Depth slice at z D 14150 ft of the migration in the k_{x} − k_{h} plane. not-filtered [ER]

Since the separation of primaries and multiples is not perfect, some of the multiples are present in the model of primaries, and some primary energy leaks to the model of multiples. Because of this cross-talk between the estimates of primaries and multiples, it is desirable to control the amount of attenuation. Therefore, the impact on the primaries of the multipleattenuation process can be ameliorated if the amplitude and phase of the multiples model are adjusted, minimizing cross-talk and the differences from the multiples present in the original data in a least-squares sense. To perform this adjustment, we use simultaneous adaptive matching of primaries and multiples, formulated by (Alvarez and Guitton, 2007).

After pre-processing, the multiples have been largely attenuated, making the migration image more suitable for inversion. Figure 12 shows the filtered migration in the subsurfaceoffset domain. As in Figure 7, the front face corresponds to the image at zero subsurface offset, and the side face corresponds to the subsurface-offset gather at a position inside the shadow zone (x D 36800 ft). This result is corroborated in the angle domain (Figure 13). The front face corresponds to the image at the 0^{∘} reflection angle, and the side face corresponds to the reflection-angle gather at a position inside the shadow zone (x D 36800 ft).

Inversion was computed using equation 3 for two different right-hand-side vectors. One is the migrated image at zero subsurface offset without filtering (front panel, Figure 7), and the other is the migrated image at zero subsurface offset after filtering (front panel, Figure 12).

The results of the inversion are shown in Figures 14 and 15 (unfiltered and filtered input), and should be compared with the migrations shown in Figures 7 and 12, respectively. The inversion with the unfiltered migration is more unstable. This should require the use of a high value of the regularization parameter, which reduces the effectiveness of the inversion outside the shadow zones.

Figure 10: Primaries after bandpassing the migration (Figure 9) in the k_{x} − k_{h} plane. filtered [ER]

Figure 11: The model for the multiples in the k_{x} − k_{h} plane, obtained after subtraction of the filtered result (Figure 10) from the migration (Figure 9). diff-filtered [ER] Figure 12: Sigsbee2b shot-profile migration (subsurface-offset) using the cross-correlation imaging condition, after multiple attenuation. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at x D 36800 ft. Compare with Figure 7. migoff1-filt [ER]

Figure 13: Sigsbee2b shot-profile migration (reflection angle) using the cross-correlation imaging condition, after multiple attenuation. The front face corresponds to the image at the 0^{∘} reflection angle, and the side face corresponds to the reflection-angle gather at x D 36800 ft. Compare with Figure 8. migang1-filt [ER] Figure 14: Inversion in the poststack image domain with the unfiltered migration as the input (equation 3) after seven iterations of a conjugate-gradient iterative solver and no regularization. invz [ER]

Figure 15: Inversion in the poststack image domain with a filtered migration as the input to equation 3. invz-filt [ER]

The inversion of the filtered migration shows a better convergence behavior. Again, no regularization was applied. Notice the well-collapsed diffractor at (35000,17000), in contrast to the unfiltered version (Figure 14). For the latter, conjugate-gradient iterations focused on reducing the multiple energy in the space of the residuals, decreasing the overall efficacy of inversion.

Figures 16 and 17 show the residuals of the inversion with the unfiltered migration and with the filtered migration, respectively, after seven iterations of the conjugate-gradient solver. Except for the residuals in the salt body, the highest residual amplitudes of the unfiltered inversion correspond to the multiples in the shadow zone. Notice that some amount of the diffractor energy (35000,17000) is still present in the residuals. On the other hand, the residuals of the inversion with the filtered data do not show the diffractor energy. Additionally, the residuals of the fault close to the diffractor are smaller than that of the unfiltered version. Unfortunately the multiples were not completely removed by filtering, so again the inversion procedure increases their amplitudes.

Overall, the inversion results after pre-processing have more balanced amplitudes, allowing the continuation of the reflector inside the shadow zones with much-improved kinematics. The reflectors also gain more vertical and horizontal resolution, particularly seen at the two faults present in the reflectivity model (Figure 2) and at the diffractor mentioned above.

CONCLUSIONS

We showed how different modes of multiple reflections present in Sigsbee2b data behave after migration. The multiple reflections we faced are not surface-related; therefore, they can still persist after multiple attenuation during processing. As the model used to compute this dataset is based on a realistic geological setting, the same behavior is likely to happen in real data. We explored the separability of primaries from multiples after migration by applying simple k_{x}−k_{h} filters. Once discriminated, primaries and multiples were submitted to a simultaneous, adaptive, non-stationary filtering to adjust amplitude and phase and decrease the cross-talk between them. The adapted multiples were the subtracted from the migrated data.

The inversion results show that, if the physics of the wave-propagation is not adequate to account for all the propagation modes in the migrated data, coherent noise – in our case multiples – will dominate the residual space, decreasing the efficiency of inversion. It is imperative, therefore, to make data satisfy the physical approximations by means of pre-processing that, in some sense reconciles processing and inversion.

REFERENCES

Alvarez, G. and A. Guitton, 2007, Simultaneous adaptive matching of primaries and multiples with non-stationary filters: SEP-125, pages 61–77.

Alvarez, G., B. Biondi, and A. Guitton, 2007, Attenuation of specular and diffracted 2D multiples in image space: Geophysics, 72, no. 5, V97–V109.

Figure 16: Residuals of the inversion in the poststack image domain with the unfiltered migration as the input after seven iterations. resinvz [ER]

Figure 17: Residuals of the inversion in the poststack image domain with a filtered migration as the input after seven iterations. resinvz-filt [ER] Chauris, H. and M. Noble, 2001, Two-dimensional velocity macro model estimation from seis-

mic reflection data by local differential semblance optimization: applications to synthetic and real data sets: Geophysical Journal International, 144, 14–26.

Clapp, M. L., 2005, Imaging Under Salt: Illumination Compensation by Regularized Inversion: Ph.D. thesis, Stanford University.

Li, J. and W. Symes, 2007, Interval velocity estimation via NMO-based differential semblance: Geophysics, 72, U75–U88.

Paffenholz, J., B. McLain, J. Zaske, and P.Keliher, 2002, Subsalt multiple attenuation and imaging: Observations from the Sigsbee2b synthetic data set: 72nd Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts, pages 2122–2125.

Sava, P. and A. Guitton, 2005, Multiple attenuation in the image space: Geophysics, 70, no. 1, V10–V20.

Shin, C. and D. Min, 2006, Waveform inversion using a logarithmic wavefield: Geophysics, 71, R31–R42.

Valenciano, A., 2008, Imaging by Wave-Equation Inversion: Ph.D. thesis, Stanford University.

Wapenaar, C. P. A., N. A. Kineging, and A. J. Berkhout, 1987, Principle of prestack migration based on the full elastic two-way wave equation: Geophysics, 52, 151–173.