Water-bottom and pegleg multiples can contaminate primary reflections and distort AVO amplitude information. We develop an inversion method for suppressing multiple reflections while simultaneously preserving amplitudes along primary reflections. An iterative time-domain conjugate gradient scheme is used to find a velocity scan which ``fits'' its associated CMP gather to within a few percent misfit energy, when a hyperbolic forward modeling operator is applied to that inverted velocity scan. This ensures that most amplitude and AVO information is preserved in the transform pair. A primary velocity trend is then automatically picked in the scan by a Monte Carlo method, and a mask is designed on this basis to isolate multiple energy in velocity space. An estimate of the multiple reflection energy is generated by forward modeling the masked velocity scan. An amplitude-preserved multiple-suppressed CMP gather is then obtained by subtracting the estimated multiples from the input CMP gather. Current test results with the Mobil AVO dataset are very encouraging.
Multiple reflections contaminate primary reflection data and distort AVO amplitude information. Many previous multiple-suppression techniques have focused on kinematic multiple-suppression, e.g., predictive deconvolution, without necessarily preserving primary amplitude information. Recently, the general increase in use of amplitudes in such applications as AVO analysis has promoted the need for amplitude-preserved multiple-suppression.
Candidates for amplitude-preserved multiple-suppression include surface-consistent techniques derived from wave-equation theory (, ). Although theoretically attractive, these techniques are expensive, and require estimation of unrecorded near and far offset traces, in addition to a good estimate of the source wavefield.
In contrast, velocity-stack operators are not strictly derived from the wave equation, but are valuable processing tools (). They rely on an assumption of weak lateral velocity variation, which results in rapid midpoint-independent processing. An ideal velocity-stack operator for the multiple-suppression application is one that (1) has an efficient and stable inverse , and (2) can best separate the multiple and primary reflection energy in velocity-stack space.
Time-invariant velocity-stack operators
Time-invariant velocity-stack operators, such as the -p slant-stack () and the parabolic Radon transform (), have Toeplitz structures in the frequency domain, and are therefore easily invertible (). This feature of the inverse transform satisfies the first design criterion for multiple suppression. Unfortunately, time-invariant slant-stack and parabolic Radon transforms do not best fit the time-varying hyperbolic nature of reflection arrivals, and therefore tend to smear together the respective energy of multiple and primary reflection events in velocity-stack space, violating the second criterion. Additionally, frequency-domain Toeplitz inversions tend to suffer from artifacts caused by space-time data truncations caused by the mute zone and dead traces. Foster and Mosher presented an efficient frequency-domain transform pair for a family of time-shifted hyperbolic operators which improves discrimination of primary and multiple reflection velocity energy, and frequency-dependent damping to help suppress mute artifacts.
Time-variant velocity-stack operators
Ideally, we would like to use a time-variant hyperbolic operator that stacks and scatters along Dix rms-velocity hyperbolas. This time-variant hyperbolic operator offers excellent separation of primary and multiple reflection energy in the corresponding velocity-stack space. Unfortunately, time-variant hyperbolic trajectories do not give rise to frequency-domain Toeplitz structure, and these operators are not analytically invertible (). However, given an operator and its adjoint, a least-squares problem can be formulated to invert the transform using an iterative time-domain conjugate gradient (CG) method ().
We present a method for amplitude-preserved multiple suppression using a time-variant hyperbolic operator and an iterative time-domain conjugate gradient velocity-stack inversion. The time-domain nature of this inversion naturally accommodates data truncations arising from the mute zone and dead traces, and allows nonstationary data- and model-space operator weighting (, ). Given this forward and inverse transform pair, we isolate the multiple reflection energy in the inverted velocity space based on moveout velocity difference. This mask is automatically designed from a Monte Carlo pick of the primary velocity trend (), and is based on an approximation of the corresponding rms velocity of the first water-column multiple. The masked velocity scan is forward-modeled with the Dix hyperbola operator to estimate the multiple reflections, which are then subtracted from the input gather to produce our best attempt at an amplitude-preserved multiple-suppressed CMP gather.
This paper proceeds as follows. First, we review the theory of the velocity-stack operator. Next, we demonstrate the least-squares velocity-stack inversion for both synthetic and marine field data. Finally, we mask the velocity scans and model the coherent multiples, and subtract to obtain amplitude-preserved multiple-suppressed CMP gathers. In both synthetic and field data tests, multiple reflections are favorably suppressed with apparently good preservation of primary reflection amplitude information.
The theory for the velocity stack inversion is addressed from two different perspectives by Lumley and Nichols . More details regarding the conjugate gradient velocity stack inversion methodology can be found therein. In this paper we will review the theory briefly and then concentrate on the multiple-suppression application.
The basis of the inversion is a time-domain hyperbolic summation/scatter operator . In this paper we will follow the convention of Nichols that is a scattering operator which maps points in the velocity space model to a CMP gather in the seismic data space :
where h is source-receiver offset, t is two-way traveltime, v is Dix rms velocity, and is two-way vertical traveltime (pseudodepth). The companion paper by Lumley assumes a hyperbolic summation operator which is the adjoint to such that .The hyperbolic operators sum/scatter along Dix rms velocity hyperbolas. The operator maps only within the unmuted data zone and recognizes dead traces. Summation and scattering use linear interpolation of trace values; offset, time and cosine operator weighting; and optional rho-filter spectral shaping.
The least-squares estimate of from (1) requires minimizing the L2 norm of the residual error,
with respect to . The formal solution to (2) is
Unfortunately, for large-scale problems typical of exploration seismic applications, neither the Hessian operator nor its inverse can be practically computed. Furthermore, this inverse may be unstable, and so an alternate solution is desirable. Nichols gives a good discussion of weighted least-squares and various Lp norm solutions to the velocity-stack problem.
The method of conjugate gradients () provides an iterative solution for which is suitable for large-scale problems and may be truncated prior to absolute convergence when a desirable inverse has been estimated. The first step of a conjugate gradient solution is an application of the adjoint operator to the data:
and each subsequent iteration proceeds to estimate the inverse Hessian . For the velocity-stack problem, is a linear version of the industry-standard semblance velocity scan (). The semblance coherency measure is not linear since the weights are a function of the squared data.
The linear gradient solution is unlikely to be a good fit to the data when forward modeled:
unless is pseudounitary (). Therefore, unlike the industry-standard velocity scan, the CG method iterates until a velocity scan is found which does fit the data after forward modeling, to within some error tolerance. For linear operators, the conjugate gradient method is guaranteed to converge in N iterations, where N is the number of unknowns . Fortunately, a satisfactory CG solution can often be estimated in a small fraction of N steps. At this point of ``practical convergence'', the L2 velocity scan can be forward-modeled with the scattering operator to accurately match the input data. Hence, an amplitude-preserving velocity transform has been found, and if multiple energy can be isolated from the velocity scan without distorting primary energy, an amplitude-preserved multiple-suppression technique will result.
SYNTHETIC AND MARINE DATA
We will demonstrate our multiple suppression methodology on marine data and a full-waveform elastic synthetic. The marine data are part of the Mobil AVO workshop, and the synthetic was modeled to approximate the marine data using a sonic and density log at a well located along the marine survey ().
Mobil AVO data
Figure shows a marine CMP from the Mobil dataset located at Well A. The data are rich in multiple reflection energy which contaminate the kinematics and AVO amplitude information along the primary reflection arrivals. The accompanying velocity semblance scan shows the primary velocity trend and several multiple velocity trends. The Mobil CMP gathers are 60-fold with a minimum and maximum offset of 262 m and 3.2 km respectively. The data are sampled at 4 msec intervals with a 6 second record length. With the multiple contamination, there is not much coherent signal energy above 60 Hz. At this location, the target horizons are at about 2.1 and 2.6 seconds traveltime respectively.
Haskell synthetic seismograms
In order to test our algorithm, we created a synthetic CMP gather using the sonic and density log measurements from Well A. We generated seismograms using a 1-D full-wave Haskell-Thompson reflectivity matrix propagator method, which includes all shear conversions, free-surface and interbed multiple reflections, and transmission effects (). Gaussian-distributed random noise was added to each synthetic dataset with a S/N amplitude ratio of 2:1, and bandpassed to the same spectrum as the Mobil data.
We created one dataset without free-surface multiples by extending the water column upwards to infinity. This represents our ideal dataset after multiple suppression, as shown in Figure . Note the late arriving shear conversions PSSP which are recorded as compressional waves in the water column. The shear conversions have zero amplitude at near offset in the CMP gather, and show up as slow velocity peaks at around 3 seconds in the semblance scan.
We also created a second dataset including free surface multiples, as shown in Figure , with which to test our multiple-suppression algorithm. A comparison of the semblance spectra show that most of the multiple energy correlates with water-column multiples which reverberate in the water layer and have reflected once at a deeper target. This phenomenon gives rise to the characteristic family of multiple velocity trends which are delayed from each other by the water depth traveltime (0.45 seconds) and shifted to slightly lower rms velocities. Three water-column reverberations are clearly visible beneath the primary trend in the velocity scan of Figure .
VELOCITY SCAN INVERSION
The first step in the multiple-suppression process is to invert the CMP data for a velocity scan that fits the input when forward modeled back to data space. This is the velocity inversion step. We show that we can obtain a good estimate of the velocity scan which fits the input CMP data on both the Haskell synthetics and Mobil marine data.
Velocity inversion of the Haskell data
We used 12 iterations of conjugate gradients to estimate the velocity scan for the Haskell data. We used the pseudounitary hyperbolic operator described by Lumley with an offset and time-variant operator weighting:
where h and t are in units of [km] and [seconds] respectively. A straight L2 norm was applied with no model weighting. Stacking velocities were calculated in the range 1.2-3.0 km/s in 30 m/s increments. The velocity stack inversion takes about 70 CPU seconds for 12 CG iterations on an HP 750 workstation per CMP gather, including all network I/O.
Figure shows the input gather containing free-surface multiples, the inverted velocity scan, and the predicted data obtained by forward modeling from the inverted velocity scan. The input and predicted data are compared in Figure alongside of the data residual, obtained by subtracting the predicted data from the input data. Figure shows that at the 12th iteration, the CG inversion has explained about 88% of the energy of the input data. The data residuals in Figure show that most of the unexplained information is random noise. However, there are some faint deep events in the residual gather which correspond to late shear arrivals which are not completely mapped in the chosen P-wave velocity range, and therefore cannot be completely modeled in the predicted data.
Velocity inversion of Mobil data
The CG algorithm was applied to the Mobil data in the same manner as the Haskell data. The operator weights were the same as in (6), and 12 iterations were deemed sufficient to reduce the residual error to an acceptable level. The velocity stack inversion takes about 70 CPU seconds for 12 CG iterations on an HP 750 workstation per CMP gather, including all network I/O.
Figure shows the raw Mobil input gather, the velocity scan inversion, and the forward-modeled predicted data. Note the strong ``sideswipe'' artifacts in the velocity scan. These correspond to the bright reflection energy near the mute line between 2 and 3 seconds in the CMP gather. We believe these events to be multiple reflections of a high-amplitude post-critical reflection at or near the seabottom. Although these artifacts are strong, the primary velocity trend is clearly present, and relatively uncontaminated.
Figure compares the Mobil input data to the predicted data, along with the data residuals. There is a faint trace of reflection energy which is spatially coherent in the residuals, but most of the residual energy is uncorrelated noise. Figure shows that at the 12th iteration, the CG inversion has explained about 93% of the energy of the input data. It is interesting to note that the inversion of the Mobil data converges to a lower relative error than the inversion of the Haskell synthetic data. This is partly because some shear energy in the Haskell data falls outside the P-wave range chosen for the scans, and therefore they cannot be modeled accurately. The Haskell data residual energy effectively plateaus at the level of the irreducible shear-wave energy and uncorrelated noise.
The next step in the procedure is to design a mask to separate primary and multiple energy in the inverted velocity scans. We tested two methods and obtained superior results when we isolated multiple energy in the velocity scan, forward-modeled multiples, and subtracted them from the input data, as opposed to direct forward modeling of primary reflections from the masked velocity scan. Modeling primaries suppresses non-hyperbolic noise along with multiples, but a large range of velocities is needed to accurately capture all the primary energy. In contrast, subtracting modeled multiples does not suppress non-hyperbolic noise, but a good range of low velocities is practical to compute in order to accurately capture all the multiple reflection energy.
The mask to isolate the multiple reflection energy is designed based on an estimate of the primary velocity trend and the first water-column multiple velocity trend. First, we need a velocity curve that best fits the primary velocity trend. This can be supplied as an initial velocity estimate, or as in our case, generated automatically by a Monte Carlo velocity picking routine applied to the smoothed envelope of the velocity scan ().
Next, we use the primary velocity estimate, and knowledge of the water depth and water velocity, to approximate the rms velocity function describing the first water-column multiples. The rms velocity of a primary reflection, , is given by the simple 1-D relation:
where vi is the interval velocity of the ith layer, is the two-way traveltime within the ith layer, and t is the total two-way traveltime from the recording surface to the reflector.
The rms velocity of the first water-column multiple reflection associated with that primary, , can be written as:
where tw is the two-way traveltime within the water layer, and vw is the water velocity. The latter equation can be Taylor-expanded as:
yielding a second order approximation:
The velocity mask is designed to isolate the multiple reflection energy such that the mask is zero for velocities above the curve, unity for velocities below the curve, and ramps in between the two curves. We use a power-law ramp that removes more energy close to the curve than the multiple trend.
Haskell data mask
Figure shows the inverted velocity scan for the Haskell multiples data, and the masked gather in which the multiple reflection energy has been isolated. The primary and multiple velocity functions are overlain on the unmasked scan to display the mask design criteria. Note that most of the multiple energy has been isolated without including much of the primary energy. Because some of the primary velocity peaks have been smeared laterally to model the near offset data truncation, a small amount of primary energy has been included in the multiple mask. This shows the need for a velocity inversion which tries to compress the energy of velocity peaks as much as possible, either by operator weighting or a good initial model estimate (), or by model weighting or a norm change (, ).
Mobil data mask
Figure shows the inverted velocity scan for the Mobil data, and the masked gather in which multiple reflection energy has been isolated. The primary and multiple velocity curves are overlain on the unmuted scan. Note that the mask velocity very nicely isolates the multiples and the strong artifacts of the multiply-reflected postcritical waterbottom reflection, leaving the primary reflections relatively well-separated.
After velocity scan masking, the masked scans can be forward-modeled to predict an estimate of the multiple reflection data. This involves one application of the hyperbolic scattering operator .The multiple-suppression is completed by subtracting the modeled multiples from the input CMP gather.
Haskell data multiple suppression
Figure shows the input Haskell multiples gather, the multiple-suppressed gather, and the removed multiples. Note that the multiples have been well-suppressed, without much apparent impact on the offset amplitude-variation along the primary reflections. Figure shows an enlarged view of the Haskell data multiple suppression results. Figure shows semblance velocity scans made on the input, multiple-suppressed, and removed multiples CMP gathers. There is a strong improvement in separations of multiples from primaries in these plots.
The toughest test, unavailable with field data, is shown in the enlarged Figure . Here, we compare the ideal Haskell synthetic modeled without free-surface multiples to the multiple-suppressed data. Aside from the shear arrivals, the multiple-suppressed data matches the Haskell primary data very well. This is evidenced in that the primaries residual data shows that most P-wave events are at or below the noise level, indicating a good kinematic match, and that there are no significant discrepancies in AVO behavior. There is a slight tendency for the residual error to be somewhat larger at near offset events, due to the inclusion of some lateral primary smear in velocity space. This could be improved by a model-space weighting in the inversion which further compacts primary and multiple velocity peaks. The strongest non-shear event in the primary residual is the partially removed multiple reflection at 1.0 seconds. This is the first water-bottom multiple and it is not well-separated in velocity space from shallow primary reflection moveout, and hence, cannot be completely suppressed.
Mobil data multiple suppression
Figure shows the results of our multiple suppression to the Mobil data, and Figure shows an enlarged view of this comparison. The input data are compared with the multiple-suppressed data and the removed multiples. There is especially good resolution of uncovered primary events at 1.8, 2.1, and 2.6 seconds, the latter two of which are reportedly hydrocarbon-bearing target zones. The velocity semblance scans of Figure show that the multiples have been well-separated from the primaries in velocity space. In fact, since semblance is a nonlinear coherency measure, it appears that a stronger primary velocity trend has been imaged in the multiple-suppressed scan as compared to the input scan. Ideally, new primary velocity curves could be estimated and the entire multiple-suppression process could be iterated. In practice, we suggest subsequent prestack processing (migration, AVO) with the multiple-suppressed gathers after an single extra updated velocity analysis estimate.
We have presented a conjugate gradient method for velocity-stack inversion and amplitude-preserved multiple suppression. A time-domain hyperbolic sum/scatter operator is used to find a least-squares estimate of a velocity scan which, when forward modeled, fits the input data to within some prescribed error tolerance. A Monte Carlo algorithm is used to automatically pick the primary velocity trend in the inverted scan, and a mask is designed on this basis to isolate multiple reflection energy from the velocity space. The masked velocity scan is forward modeled to predict the multiple reflections, which are then subtracted from the input data to produce an amplitude-preserved multiple suppression. Current tests with full-waveform synthetic data and marine data from the Mobil AVO workshop appear very encouraging.
We intend to further explore the practical utility of other norms and model-space weights for further separating primary and multiple velocity energy. We plan to apply a tuned version of this algorithm to the entire 2000 CMP gathers in the Mobil AVO dataset and present our results at the SEG workshop in Los Angeles.
We thank Bob Keys and Stew Levin at Mobil Oil for supplying the AVO workshop data.