Prediction error filters to enhance differences

Methodology

PEFs attempt to capture the inverse spectrum of the data. In the 1-D case, we could calculate a filter with the inverse spectrum by transforming into the frequency domain and then doing a sample by sample division,

$$Y(\omega) = \frac{1}{D(\omega)} ,$$ (1)

where $Y(\omega)$ is the filter and $D(\omega)$ is the data in the Fourier domain. There are two problems with this approach. First $D(\omega)$ can be small or zero valued and $y(t)$ is not compact. Claerbout (1999) shows that a compact filter can be estimated by solving the least squares system

$$\bf y= ( \bf D'\bf D)^{-1} \bf D\bf0,$$ (2)

where $\bf y$ is a filter whose zero lag is fixed at one, $\bf0$ is a vector of 0s, and $\bf D$ is convolution with the data. In general the shape of filter is arbitrary but needs to be large enough to capture the spectrum of the wavelet and the dips present in the volume.

The residual $\bf r$ of the estimation procedure can be calculated by convolving the filter $\bf Y$ with the data

$$\bf r= \bf Y\bf d .$$ (3)

The residual will be large when the filter is not large enough to fully describe the stationary spectrum or the data is non-stationary. The left panel of Figure 1 shows the result of migrating the Marmousi dataset using a standard downward continuation based migration. The right panel shows the result of first estimating a series of PEFs in overlapping patches on the migrated image and then applying the PEFs to the migrated volume (applying equation 3). Note the areas of large residual generally correspond to unconformities and fault locations.

The next step is to apply this same series of PEFs to another image. The first question is what happens if we apply a filter estimated on volume `a' to volume `b' which has significantly different spatial statistics. The left panel of Figure 2 shows a simple plane wave. The right panel of Figure 2 shows the result applying the filters estimated from the Marmousi migration. The dominant feature is still the planewave. The amplitude of the residual is on average an order of magnitude higher than the residual shown in the right panel of Figure 1. If dataset `b' has a spectrum close to `a' we get a different result. We should see large values at both where the stationarity assumption of the PEF is invalid and at places where the covariance description of `a' and `b' are different. Figure 3 illustrates this point. Both the left panel of Figure 1 and 3 are calculated by a source-receiver Phase-Shift Plus Interpolation (PSPI) algorithm. The left panel of Figure 1 shows the result of using up to eight reference velocities, the left of panel of Figure 3 uses a single reference velocity at each depth step. The right panel of Figure 3 shows the result of applying equation3 using the filter calculated from the eight velocity migration. Note that in addition to the large residual locations seen in the right panel of Figure 1, we now see additional locations. Generally the large values are at and below areas of large dip, where the first order split step correction is least accurate.

base
Figure 1. The left panel the result of PSPI migration of the Marmousi dataset using 8 reference velocities. The right panel shows the result of equation 3.

different
Figure 2. The left panel is a planewave. The right panel is the result of applying the filter estimated from the 8 velocity Marmousi image.

onevel
Figure 3. The left panel is the result of PSPI migration of the Marmousi dataset using one reference velocity. The right panel shows the result of equation 3 using filters calculated from the data shown in the left panel of Figure 1.

What we really would like is just the differences caused by the change in the migration algorithm. For notational convenience we will define $\bf r_{a,b}$ as the residual of applying a filter estimated on dataset `a' to dataset `b'. Simply dividing $\bf r_{a,a}$ by $\bf r_{a,a}$ is not feasible due to the zero in $\bf r_{a,a}$. One approach to this problem is adding an epsilon term to the denominator. Another approach is smoothing. We first take the absolute value $\bf A$, and then smooth the resulting volume. As a result, we end up with an estimate of the fitting error $\bf e$,

$$\bf e=\frac{ \bf S \bf A \bf r_{a,b}}{\bf S \bf A \bf r_{a,a}}- \bf 1$$ (4)

An alternate approach is to add a scaling term that emphasizes errors where the original data is large. We can do this by applying a smoothed envelope function $\bf E$ to dataset $\bf a$,

$$\bf e_{{\rm scaled}} = \bf E \bf a \left( \frac{ \bf S \bf A \bf r_{a,b}} {\bf S \bf A \bf r_{a,a}} -\bf 1 \right) .$$ (5)

Figure 4 shows the result of applying equation 5 comparing the one and eight reference velocity images. Note how the differences are located at steep dips, where we would anticipate the single reference velocity approach failing.

error-onevel Figure 4. The result of applying equation 5 comparing the one and eight reference velocity images.

Prediction error filters to enhance differences