WELLS NOT MATCHING THE SEISMIC MAP

Accurate knowledge comes from a well, but wells are expensive and far apart. Less accurate knowledge comes from surface seismology, but this knowledge is available densely in space and can indicate significant trends between the wells. For example, a prospective area may contain 15 wells but 600 or more seismic stations. To choose future well locations, it is helpful to match the known well data with the seismic data. Although the seismic data is delightfully dense in space, it often mismatches the wells because there are systematic differences in the nature of the measurements. These discrepancies are sometimes attributed to velocity anisotropy. To work with such measurements, we do not need to track down the physical model, we need only to merge the information somehow so we can appropriately map the trends between wells and make a proposal for the next drill site. Here we consider only a scalar value at each location. Take $\bold w$ to be a vector of 15 components, each component being the seismic travel time to some fixed depth in a well. Likewise let $\bold s$ be a 600-component vector each with the seismic travel time to that fixed depth as estimated wholly from surface seismology. Such empirical corrections are often called ``fudge factors''. An example is the Chevron oil field in Figure 6.

 

wellseis
Figure 6 Binning by data push. Left is seismic data. Right is well locations. Values in bins are divided by numbers in bins. (Toldi)

The binning of the seismic data in Figure 6 is not really satisfactory when we have available the techniques of missing data estimation to fill the empty bins. Using the ideas of subroutine mis1() , we can extend the seismic data into the empty part of the plane. We use the same principle that we minimize the energy in the filtered map where the map must match the data where it is known. I chose the filter $\bold A = \nabla'\nabla=-\nabla^2$to be the Laplacian operator (actually, its negative) to obtain the result in Figure 7.

 

misseis
Figure 7 Seismic binned (left) and extended (right) by minimizing energy in $\nabla^2 \bold s$.

Figure 7 also involves a boundary condition calculation. Many differential equations have a solution that becomes infinite at infinite distance, and in practice this means that the largest solutions may often be found on the boundaries of the plot, exactly where there is the least information. To obtain a more pleasing result, I placed artificial ``average'' data along the outer boundary. Each boundary point was given the value of an average of the interior data values. The average was weighted, each weight being an inverse power of the separation distance of the boundary point from the interior point.

Parenthetically, we notice that all the unknown interior points could be guessed by the same method we used on the outer boundary. After some experience guessing what inverse power would be best for the weighting functions, I do not recommend this method. Like gravity, the forces of interpolation from the weighted sums are not blocked by intervening objects. But the temperature in a house is not a function of temperature in its neighbor's house. To further isolate the more remote points, I chose weights to be the inverse fourth power of distance.

After the gaps are filled in the seismic information $\bold s$,we are prepared to find another map $\bold m$that matches the wells $\bold w$ exactly and tries to match the seismic information elsewhere. Now we must beware that the seismic information is likely to be inconsistent with the wells. We presume that the two should fit at high spatial frequencies but mismatch where there are wells. Think of a map of a model space $\bold m$of infinitely many hypothetical wells that must match the real wells, where we have real wells. Otherwise, we seek the minimum residual $\bold r$which is the filtered difference between the seismic data and the map $\bold m$ of hypothetical omnipresent wells. Thus we seek to fit  

$$\bold 0\quad\approx\quad \bold r \eq \bold A ( \bold m - \bold s )$$ (14)

along with the constraint

$$\bold 0 \eq (\bold I -\bold J_w) ( \bold m - \bold w )$$ (15)

which means the model matches the wells where we have wells. To solve this problem, we convert it to one that we have already solved. We define a new variable $\bold y$

$$\bold y \eq \bold m - \bold s$$ (16)

and we seek to minimize

$$\bold 0\quad\approx\quad \bold r \eq \bold A \bold y$$ (17)

Initialize $\bold y_0$ so that where(known) y0=(w-s), in other words,

$$0 \eq (\bold I-\bold J_w)( \bold y_0 \ -\ (\bold w-\bold s) )$$ (18)

Solve for $\bold y$ by the method above, and then we have our solution $\bold m = \bold y + \bold s$.

We have already coded a 2-D gradient operator igrad2 . Let us combine it with its adjoint to get the 2-D laplacian operator. (You might notice that the laplacian operator is ``self-adjoint'' meaning that the operator does the same calculation that its adjoint does. Any operator of the form $\bold A'\bold A$ is self-adjoint because $(\bold A'\bold A)'=\bold A'\bold A''=\bold A'\bold A$. )

 

module laplac2 {                                    # Laplacian operator in 2-D
logical, parameter, private  :: T = .true., F = .false.
use igrad2
real, dimension (m1*m2*2), allocatable :: tmp
#%_init    (m1, m2)
    integer m1, m2
    call igrad2_init (m1, m2)
#%_lop (x, y)
    integer stat1, stat2
    if( adj) {
    	stat1 = igrad2_lop (  F,  F, y, tmp)   # tmp = grad y
    	stat2 = igrad2_lop (  T,  T, x, tmp)   # x = x + grad' tmp
    } else   {
	stat1 = igrad2_lop (  F,  F, x, tmp)   # tmp = grad x
    	stat2 = igrad2_lop (  T,  T, y, tmp)   # y = y + grad' tmp
    }
}

Subroutine lapfill2() is the same idea as mis1() except that the filter $\bold A$ has been specialized to the laplacian implemented by module laplac2 .

Subroutine lapfill2() solves two problems: (1) extending the seismic data to fill space, and (2) fitting the map approximately to the seismic data while exactly to the wells.

 

module lapfill {   # fill empty 2-D bins by minimum output of Laplacian operator
  use laplac2
  use cgstep_mod
  use solver_mod
contains
  subroutine lapfill2( niter, m1, m2, yy, mfixed) {
    integer,                    intent (in)     :: niter  # iterations
    integer,                    intent (in)     :: m1,m2  # data size
    logical, dimension (m1*m2), intent (in)     :: mfixed # mask for known
    real,    dimension (m1*m2), intent (in out) :: yy     # model
    real,    dimension (m1*m2)                  :: zero   # laplacian output
    call laplac2_init ( m1,m2);            zero = 0.      # initialize

call solver ( laplac2_lop, cgstep, yy,  zero, niter,
				  x0 = yy,  known = mfixed)

call laplac2_close ()                                 # garbage collection
    call cgstep_close ()                                  # garbage collection
  }
}

The final map is shown in Figure 8.

 

finalmap
Figure 8 Final map based on Laplacian roughening.

Results can be computed with various filters. I tried both $\nabla^2$ and $\nabla$.There are disadvantages of each, $\nabla$ being too cautious and $\nabla^2$ perhaps being too aggressive. Figure 9 shows the difference between the extended seismic data and the extended wells. Notice that for $\nabla$ the difference shows a localized ``tent pole'' disturbance about each well. For $\nabla^2$ there could be large overshoot between wells, especially if two nearby wells have significantly different values. I don't see that problem here.

 

diffdiff
Figure 9 Difference between wells (the final map) and the extended seismic data. Left is plotted at the wells (with gray background for zero). Center is based on gradient roughening and shows tent-pole-like residuals at wells. Right is based on Laplacian roughening.

To understand the behavior theoretically, recall that in one dimension the filter $\nabla$ interpolates with straight lines and $\nabla^2$ interpolates with cubics. This is because the fitting goal $\bold 0 \approx \nabla \bold m$,leads to ${\partial \over\partial \bold m'} \bold m'\nabla'\nabla \bold m = \bold 0$or $\nabla'\nabla \bold m = \bold 0$, whereas the fitting goal $\bold 0 \approx \nabla^2 \bold m$leads to $\nabla^4 \bold m = \bold 0$which is satisfied by cubics. In two dimensions, minimizing the output of $\nabla$gives us solutions of Laplace's equation with sources at the known data. It is as if $\nabla$ stretches a rubber sheet over poles at each well, whereas $\nabla^2$ bends a stiff plate.

Just because $\nabla^2$ gives smoother maps than $\nabla$does not mean those maps are closer to reality. This is a deeper topic, addressed in later chapters. It is the same issue we noticed when comparing figures -.