### Abstract

We present a new approach to the joint simultaneous inversion of *PP* and *PS* angle gathers for the estimation of Pimpedance, S-impedance and density. Our algorithm is based on three assumptions. The first is that the linearized approximation for reflectivity holds. The second is that PP and PS reflectivity as a function of angle can be given by a set of linearized equations. The third is that the background trend can be described by a linear relationship between the logarithm of P-impedance and both S-impedance and density. Given these three assumptions, we show how a final estimate of P-impedance, S-impedance and density can be found by perturbing an initial low-frequency model. After a description of the algorithm, we then apply our method to both a synthetic and real data set.

### Introduction

The goal of pre-stack seismic inversion is to obtain reliable estimates of P-wave velocity (*V _{P}*), S-wave velocity (

*V*), and density (

_{S}*ρ*) from which to predict the fluid and lithology properties of the subsurface of the earth. One solution to this problem was discussed by Russell et al. (2006) in a special edition of the Recorder, which was based on the SEG expanded abstract by Hampson et al. (2005). Figure 1, taken from Russell et al. (2006) shows the fundamental seismic concept behind pre-stack inversion, in which we show the mode-converted reflected and transmitted waves from an incident seismic P-wave at the boundary between two elastic layers in the subsurface of the earth, with an incident angle equal to

*θ*

_{1}.

The computation of the amplitudes of the reflected P and Swaves was first solved by Zoeppritz (1919). Using the Zoeppritz equations, and a knowledge of the angle of incidence of the input P-wave as well as the velocities and densities of the two layers (noting that the other angles can be computed using the angle of incidence and the velocities using Snell’s law), we can easily perform forward modeling at the layer boundary to compute the amplitudes of the reflected and transmitted waves. This can also be generalized to multiple layers for both the primaries-only case and also the more complete case in which transmitted, converted and multiple reflections are taken into account (Fuchs and Müller, 1971, Kennett, 1983). However, the inverse problem of finding the velocities given the amplitudes is much more difficult than the forward problem, given the nonlinear nature of the Zoeppritz equations. For that reason, Hampson et al. (2005) perform pre-stack inversion by starting with the linearized Aki-Richards equation (Aki and Richards, 2002, Richards and Frasier, 1976) which was re-formulated algebraically by Fatti et al. (1994) as

where *c*_{1}=1+tan^{2}*θ* , *c*_{2}= –8*γ*^{2}sin^{2}*θ* , *γ*=*V _{S}*/

*V*,

_{P}*c*

_{3}= –0.5tan

^{2}

*θ*+

^{2}

*γ*

^{2}sin

^{2}

*θ*, and the linearized

*P*-reflectivity (

*R*), S-reflectivity (

_{P}*R*) and density reflectivity (

_{S}*R*) are given by

_{D}It should be pointed out that the inversion problem was also discussed by Simmons and Backus (1996), who inverted for the linearized P-reflectivity (*R _{P}*), S-reflectivity (

*R*) and density reflectivity (

_{S}*R*) terms of equations 2 through 4, using the Aki- Richards linearized approximation of equation 1. Simmons and Backus (1996) also assumed that

_{D}*ρ*and

*V*are related by Gardner’s relationship (Gardner et al. 1974), given by

_{P}and that *V _{S}* and

*V*are related by Castagna’s equation (Castagna et al., 1985), given by

_{P}The authors then use a linearized inversion approach to solve for the reflectivity terms given in equations 2 through 4.

Buland and Omre (2003) use a similar approach which they call Bayesian linearized AVO inversion. Unlike Simmons and Backus (1996), their method is parameterized by the three terms, Δ*V _{P }*/

*V*, Δ

_{P}*V*/

_{S }*V*, and Δ

_{S}*ρ*/

*ρ*, again using the Aki-Richards approximation. The authors also use the small reflectivity approximation to relate these parameter changes to the original parameter itself, or

where ln represents the natural logarithm. Similar terms can be given for changes in both *S*-wave velocity and density. This logarithmic approximation allows Buland and Omre (2003) to invert for velocity and density, rather than reflectivity, as in the case of Simmons and Backus (1996). Hampson et al. (2005) inverted directly for the *P* and *S*-impedance terms and density, using a similar small reflectivity approximation to the one used by Buland and Omre (2003). They also used modified relationships similar to those of equations 5 and 6.

In the present study, we extend the simultaneous inversion method presented by Hampson et al. (2005) to the joint inversion of *PP* and *PS* gathers, so we refer to the method as joint simultaneous inversion. In the next section we will start by reviewing post stack inversion, then we will review the full theory of prestack *PP* simultaneous inversion, and then we will extend the method to *PP* and *PS* data.

### Theory

We start by reviewing the principles of model-based post-stack inversion (Russell and Hampson, 1991). First, by combining equations 1 and 7, we can show that the small reflectivity approximation for the P-wave reflectivity is given by

where i represents the interface between layers i and i+1. If we consider an N sample reflectivity, equation 8 can be written in matrix form as

where L_{Pi} = ln(Z_{Pi}).

Next, if we represent the seismic trace as the convolution of the seismic wavelet with the earth’s reflectivity, we can write the result in matrix form as

where S_{i} represents the i^{th} sample of the seismic trace and w_{j} represents the j^{th} term of an extracted seismic wavelet. Combining equations 9 and 10 gives us the forward model which relates the seismic trace to the logarithm of P-impedance:

where W is the wavelet matrix given in equation 10 and D is the derivative matrix given in equation 9. If equation 11 is inverted using a standard matrix inversion technique to give an estimate of L_{P} from a knowledge of S and W, there are two problems. First, the matrix inversion is both costly and potentially unstable. More importantly, a matrix inversion will not recover the low frequency component of the impedance. An alternate strategy, and the one adopted in our implementation of equation 11, is to build an initial guess impedance model and then iterate towards a solution using the conjugate gradient method.

We can now extend the theory to the pre-stack inversion case. For a given angle trace S(θ) we can therefore extend the zero offset (or angle) trace given in equation 11 by combining it with equation 1 to get

where L_{Pi} = ln(Z_{S}) and L_{D} = ln(ρ). Note that the wavelet is now dependent on angle. Equation 12 could be used for inversion, except that it ignores the fact that there is a background relationship between L_{P} and L_{S} and between L_{P} and L_{D}. Because we are dealing with impedance rather than velocity, and have taken logarithms, our relationships are different than those given by Simmons and Backus (1996) and are given by

and

That is, we are looking for deviations away from a linear fit in logarithmic space. This is illustrated in Figure 2.

Combining equations 12 through 14, we get

where

(Note that when we differentiate equations 13 and 14, the constant terms k_{c} and m_{c} go to zero and are therefore not present in equation 15). Equation 15 can be expressed in matrix form as

If equation 16 is solved by matrix inversion methods, we again run into the problem that the low frequency content cannot be resolved. A practical approach is to initialize the solution to [L_{P} ΔL_{S} ΔL_{D}]^{T} = [ln(Z_{P0}) 0 0]^{T}, where Z_{P0} is the initial impedance model, and then to iterate towards a solution using the conjugate gradient method.

We will now extend the previous derivation to include pre-stack converted-wave measurements (*PS* gathers) that have been converted to *PP* time. We will use the linearized form of the equation that was developed by Aki, Richards, and Frasier (Aki and Richards, 2002, Richards and Frasier, 1976). It has been shown by Margrave et al. (2001) that this equation can be written as

where

and

The reflectivity terms R_{S} and R_{D} given in equation 17 are identical to the terms given in equations 3 and 4. Using the small reflectivity approximation, we can therefore re-write equation 17 as:

Next, using the relationships between S-impedance, density and P-impedance given in equations 13 and 14, equation 18 can be further re-written as

where

Note that equation 19 allows us to express a single *PS* angle stack as a function of the same three parameters given in equation 15. Also, equation 19 is given at a single incident angle *θ*. When we generalize this equation to *M* angle stacks, we can combine this relationship with equation 16 and write the general matrix equation as

Equation 20 gives us a general expression for the simultaneous inversion of *N PP* angle stacks and *M PS* angle stacks. Note that we extract a different wavelet for each of the *PS* angle stacks, as was done for each of the *PP* angle stacks.

### Model Example

We will now apply this method to a model data example, shown in Figure 3. The model for these inversion tests is the Ostrander model (Ostrander, 1984), which consists of a gas sand between two identical shales, in which both the *V _{P}* /

*V*ratio and acoustic impedance of the gas sand are less than that of the shales (this is called a class 3 AVO anomaly). The synthetic data is calculated using the Aki-Richards equations given in equations 1 (for the

_{S}*PP*synthetic) and 17 (for the

*PS*synthetic). Traces are modeled from 0 to 45 degrees. For the display shown in Figure 3, the PS data has been converted to PP time. The wavelet is a band-pass filter with cutoff frequencies of 5/10-50/60.

The reason for using such a simple synthetic example is that we know exactly what the correct answer is, so that as we vary the key parameters in the inversion we can quantify the effect of these parameters on the inversion result. In particular, we have found two issues which can affect the inversion result:

- Constant gamma (
*V*/_{S}*V*ratio) versus variable gamma, and_{P} - Whether the PS data has been transformed to PP time or left in the original PS time.

We will address both of these issues in the following tests. For point number 1, note that the scaling terms c_{2} and c_{3} in equation 1 and c_{4} and c_{5} in equation 17 are dependent on γ, the* V _{S}* /

*V*ratio, and that this ratio should vary as a function of depth or time. Thus, approximating this ratio as a constant for the sake of stability in the algorithm will actually produce a slightly incorrect result. In our first implementation of this method, the variable γ was assumed to be constant and supplied by the user. A second approximation in this implementation was that the angle θ was taken as the incident angle, when actually it should be the average between the incident and transmitted angle, as discussed above. As just mentioned, the motivation for both of these approximations was that they allowed the c terms to be constants with time, which sped up and stabilized the calculation. To address these problems, we therefore added an option to calculate γ automatically instead of letting the user input the value. When this option is used, the following things happen:

_{P}- A time-variant gamma is calculated from the initial guess model and used to start the inversion.
- Every five iterations during the inversion process, the timevariant gamma is updated with the current values from the inversion. This is an approximate solution to the non-linear problem created with one of the “unknowns” being in the coefficients.
- The correct average angle is used for all calculations, using the current time-variant velocity.
- The process runs much slower than the constant gamma assumption.

To illustrate the difference between constant and variable values of gamma, Figure 4 shows the results on the logs of inverting the synthetic using a constant value of γ = 0.5. The figure shows the results (in red) after 50 iterations, starting from a smooth initial guess (shown as the smooth black line). The correct answer is shown as the blocky blue line on the same plot as both the initial guess and the final result. Before the inversion, the *PS* data had been transformed to PP time.

In Figure 4, it is obvious that the final result is reasonably good but not exact, and is different on each well log. The fits to the *P*-impedance, density and *V _{P}* /

*V*ratio are actually much better than the fit to the S-impedance, where there is an obvious decrease, rather than increase, in the gas sand (that is the event between 300 and 350 ms). To quantify the fit, the RMS error has been written at the top of each well log. Note that the P-impedance has an RMS error of 432 (m/s)*(g/cc), the S-impedance has an RMS error of 197 (m/s)*(g/cc), the density has an RMS error of 0.015 (g/cc) and the

_{S}*V*/

_{P}*V*ratio has an RMS error of 0.049. Figure 5 shows the synthetic seismogram results from inverting the synthetic using a constant value of γ = 0.5. The figure shows the

_{S}*PP*and

*PS*synthetics on the left in red, the

*PP*and

*PS*original seismic gathers in the centre in black, and the

*PP*and

*PS*error in red on the right, where the error was computed by subtracting the seismic gathers from the synthetic gathers. Notice that the

*PP*error is largest on the event at 350 ms at the far offsets and that the

*PS*error is largest at the deepest event and is fairly uniform across all offsets.

Next, Figure 6 shows the results on the logs of inverting the synthetic using a variable value of γ. The figure shows the results (in red) after 50 iterations, starting from a smooth initial guess (shown as the smooth black line). The correct answer is shown as the blocky blue line on the same plot as both the initial guess and the final result. Before the inversion, the *PS* data had been transformed to *PP* time. This time, the fit to all four logs is quite reasonable visually, and the fit to the S-impedance after inversion shows a much improved response through the gas sand zone. Again, we have quantified the fit by posting the RMS errors at the top of each well log. As expected, the error on the S-impedance has dropped dramatically. However, the error on the other three logs after inversion has gone up.

Figure 7 shows the synthetic seismogram results from inverting the synthetic using a variable value of γ. The figure shows the *PP* and *PS* synthetics on the left in red, the *PP* and *PS* original seismic gathers in the centre in black, and the *PP* and *PS* error in red on the right, where the error was computed by subtracting the seismic gathers from the synthetic gathers. Notice that when compared with Figure 5, the data error using a constant gamma value, the error has been dramatically reduced, especially at the far offset.

Next, we will turn our attention to the second problem identified above, which is whether the *PS* data has been transformed to *PP* time or left in the original *PS* time. Our inversion code so far has assumed that the *PS* data has been converted to *PP* time. This was convenient for a variety of reasons. However, the transformation of the *PS* data to *PP* time is generally a stretch/squeeze process, unless we have the unrealistic case of constant *V _{P}* /

*V*ratio background. This produces an effective time-variant wavelet, as pointed out by Gaiser (2013). The current implementation approximates the time-variant wavelet by a constant wavelet extracted from

_{S}*PS*data in

*PP*time. An alternate formulation leaves the PS data in PS time, and performs the inversion on both data sets simultaneously. This is illustrated in Figure 8, in which the two synthetics are shown on the left side of the figure in

*PP*time and are shown on the right side of the figure in

*PP*and

*PS*time, respectively.

Figure 9 shows the results on the logs of inverting the synthetic using *PP* time for the *PP* data and *PS* time for the P*S* data. The figure shows the results (in red) after 50 iterations, starting from a smooth initial guess (shown as the smooth black line). The correct answer is shown as the blocky blue line on the same plot as both the initial guess and the final result. The fit to all four logs is quite reasonable visually, and the fit to the S-impedance after inversion shows the best response through the gas sand zone of all of our tests.

Again, we have quantified the fit by posting the RMS error at the top of each well log. This shows a drop in error on all four logs from the result in Figure 6. Our conclusion is that leaving the *PS* data in *PS* time gives an improved result, as predicted by the theory. However, the run-time can be significantly longer.

### PP-PS Data Example

In this real data example, we will show a joint PP-PS inversion using a case study from northeastern Alberta. The workflow for this example involved the following steps:

- Correlate the
*PP*and*PS*data to wells - Pick the corresponding horizons on both datasets.
- Use horizon based event matching to establish the mapping between
*PP*and*PS*time. - Invert
*PP*and*PS*data using simultaneous inversion.

Figure 10(a) shows the correlated *PP* data and Figure 10(b) shows the correlated *PS* data from our example, where the *V _{P}*,

*V*, and density logs are shown on the left of each figure, and the stacked data (

_{S}*PP*or

*PS*) is shown on the right. The synthetic trace (blue) to seismic trace (red) correlation is shown between the logs and the stacked data. Note that the ties are quite reasonable in both cases shown. The wavelets used to create the synthetic were extracted from the seismic data.

Finally, we applied joint inversion to this real dataset from Alberta, using the theory discussed in earlier in this paper. In this case we have only two inputs: the *PP* stack and the *PS *stack. We assume that the *PP* stack is equivalent to a *PP* angle gather at zero degrees and that the *PS* stack is equivalent to a *PS* angle stack at 15 degrees. The value of 15 degrees was arrived at by cross-correlating the traces in the *PS* stack with model traces at various angles and selecting the maximum correlation angle.

The quality of the inversion was evaluated by displaying inversion results at the well location overlaid on the corresponding real log curves. Figure 11 shows an example of such a display. Here, we are displaying P-impedance, S-impedance, Density, and *V _{P}* /

*V*ratio. In each case, the real log is in blue, while the inversion result is in red. At the top of each track, the RMS error, the average difference between the real and estimated logs, is listed. For example, we can see that the average error in Pimpedance is 419 (m/s)*(g/cc).

_{S}In total, four different inversions were run to test the parameters described in the previous section:

- Variable background gamma (Vs/Vp) with the
*PS*data transformed to*PP*time. - Variable background gamma (Vs/Vp) with the
*PS*data in its original*PS*time. - A constant background gamma (Vs/Vp = 0.6), with the
*PS*data in*PP*time. - A constant background gamma (Vs/Vp = 0.6), with the
*PS*data in*PS*time.

The constant background gamma of 0.6 was selected as the best from a series of test values for this parameter.

For each inversion test, the RMS errors were calculated at the well location using the method described in Figure 11. These results are tabulated in Figure 12, showing the errors for the estimated Zp, Zs and Density. For each case, the error is the RMS difference between the real curve and the inverted curve.

The test cases in Figure 12 are ordered with decreasing error, or improved inversion result, from top to bottom. We can see that performing the inversion with the *PS* data in *PS* time is consistently better than transforming the *PS* data to *PP* time. This conforms with the theoretical expectation. On the other hand, surprisingly, the variable gamma tests are not as good as the constant gamma of 0.6. This probably means that the actual *V _{P}* /

*V*background is slowly varying and different from that calculated from the initial model.

_{S}The result of applying joint inversion to the seismic datasets is shown in Figure 13, where we are displaying the inverted *V _{P}* /

*V*ratio. Note that the

_{S}*V*/

_{P}*V*ratio in colour on the inserted well log matches the colour on the inversion results quite well. The colour scale in shown on the right hand side of the plot.

_{S}### Conclusions

We have presented a new approach to the joint simultaneous inversion of pre-stack seismic data which produces estimates of P-impedance, S-Impedance and density. This method allows us to incorporate both *PP* and *PS* data into the solution, if we have first calibrated both datasets. The method is based on three assumptions: that the linearized approximation for reflectivity holds, that reflectivity as a function of angle can be given by the Aki-Richards equations, and that the background trend can be described by a linear relationship between the logarithm of Pimpedance and both S-impedance and density. Our approach was shown to work well for a modelled gas sand and for a joint *PP-PS* dataset from northeast Alberta. In the first example, prestack *PP* and *PS* synthetic data was used. In the joint *PP-PS* dataset the prestack data was not available and we applied the method to full stacks.

## Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

## Share This Article