next up previous print clean
Next: VIRTUAL-RESIDUAL PRECONDITIONING Up: Preconditioning Previous: INVERSE LINEAR INTERPOLATION

EMPTY BINS AND PRECONDITIONING

Chapter [*] formulates the basic principle of missing data interpolation as follows:

A method for restoring missing data is to ensure that the restored data, after specified filtering, has minimum energy.
Mathematically, this principle can be expressed by the simple equation  
 \begin{displaymath}
D m \approx 0\;,\end{displaymath} (26)
where m is the data vector, and D is the specified filter. The approximate equality sign means that equation (26) is solved by minimizing the squared norm (the power) of its left side. Additionally, the known data values must be preserved in the optimization scheme. Introducing the mask operator K, which can be considered as a diagonal matrix with zeros on the missing data locations and ones elsewhere, we can rewrite equation (26) in the more rigorous form  
 \begin{displaymath}
D (I-K) m \approx - D K m = -D m_k\;,\end{displaymath} (27)
in which I is the identity operator, and mk is the known portion of the data. It is important to note that equation (27) corresponds to the limiting case of the regularized linear system  
 \begin{displaymath}
\left\{\begin{array}
{l}
K m = m_k\;, \\ \lambda D m \approx 0\end{array}\right.\end{displaymath} (28)
for the scaling coefficient $\lambda$ approaching zero. This means that we put far more weight on the first equation in (28) and use the second equation only to constrain the null space of the solution. Applying the general theory of preconditioning, one can immediately transform system (28) to the equation  
 \begin{displaymath}
K P x \approx m_k\;,\end{displaymath} (29)
where P is a preconditioning operator, and x is a new variable, connected with m by the simple relationship

\begin{displaymath}
m = P x\;.\end{displaymath}

If D is represented by a discrete convolution, the natural choice for P is the corresponding deconvolution operator:  
 \begin{displaymath}
P = D^{-1}\;.\end{displaymath} (30)
The helix transform provides a constructive way of implementing multidimensional deconvolution by one-dimensional recursive filtering.

The following module, mis2 [*], does the actual interpolation.

 

module mis2 {
  use helicon
  use polydiv
  use mask1
  use cgstep_mod
  use solver_mod
contains
# fill in missing data by minimizing power out of a given filter.
# by helix magic works in any number of dimensions
  subroutine mis1( niter, xx, aa, lag, known,  precon, trans) {
    logical,                intent( in)     :: precon, trans
    integer,                intent( in)     :: niter
    logical, dimension( :), intent( in)     :: known
    real,    dimension( :), intent( in out) :: xx
    real,    dimension( :), pointer         :: aa
    integer, dimension( :), pointer         :: lag
    logical, dimension( :), pointer         :: mm
    real,    dimension( :), allocatable     :: yy, dd
    integer                                 :: nx, ny, pad, n1, n2
    nx = size( xx); pad = maxval( lag); ny = nx + 2*pad
    allocate( yy( ny), mm( ny)); yy = 0.; mm = .false. 
    if (trans) { n1 = 1;       n2 = nx     # emulate transient boundaries 
       mm( ny-pad+1:ny)  = .true.          # zero values beyond end.
    } else {     n1 = 2*pad+1; n2 = ny}
    yy (n1:n2) = xx; mm (n1:n2) = known
    if( precon) {                          #  KP p = K y,   m = P p
       call mask1_init( mm)
       call polydiv_init( aa, lag)
       call solver_prec( mask1_lop, cgstep, niter= niter, x= yy, dat= yy, 
                    prec= polydiv_lop, nprec= ny, eps= 0.)
    } else {                               #  KA m = 0
       allocate( dd( ny));  dd = 0.
       call helicon_init( aa, lag)
       call solver( helicon_lop, cgstep, niter= niter, x = yy, dat = dd, 
                    known = mm, x0 = yy)
       deallocate( dd)
    }
    call cgstep_close( )
    where( xx == 0.) xx = yy(n1:n2)
    deallocate( yy, mm)
  }
}

The operator mask1, apearing in mis2 [*]. is a simple masked assignment (copy under mask.)

 

module mask1 { # masking operator
  logical, dimension( :), pointer :: m
#% _init( m)
#% _lop( x, y)
  if( adj) 
	   where( m) x += y
  else  #
           where( m) y += x
}

 
fill
fill
Figure 7
Views of the bottom of the Sea of Galilee. Missing data interpolation with preconditioning (right) and without preconditioning (left) after 40 and 2 iterations (top and bottom plots). (Fomel)


view burn build edit restore


next up previous print clean
Next: VIRTUAL-RESIDUAL PRECONDITIONING Up: Preconditioning Previous: INVERSE LINEAR INTERPOLATION
Stanford Exploration Project
2/27/1998