TIME-VARIABLE DECONVOLUTION

Filters sometimes change with time. We sometimes observe signals whose spectrum changes with time. Simple modifications to earlier convolution subroutines such as icaf1() lead to subroutine nicaf1(), which allows for a different filter on each output data point.

 

module nicaf1 { # Nonstationary Internal Convolution, Adjoint is the Filter
  real,    dimension (:), pointer :: yy
  integer                         :: n, na
#% _init (yy, n, na)
#% _lop (aa (n,na), rr (n))
     integer a, r, y
     do a = 1, na {
     do r = na, n {                     y = r - a + 1 
        if( adj)
           aa (r, a) += rr (r)    * yy (y)
        else
           rr (r)    += aa (r, a) * yy (y)
     }}
}

Compared to those in icaf1(), these filter coefficients take much computer memory because there is a separate filter for every data point. Normally, the number of filter coefficients is many fewer than the number of data points, but here we have very many more. Indeed, there are na times more. Allowing time variation takes us from overdetermined to very undetermined. We can estimate all these filter coefficients by the usual deconvolution formulation (18)

$$\bold 0 \quad\approx\quad \bold r \eq \bold Y \bold K \bold a +\bold r_0$$ (55)

but we need to supplement it with some damping equations, say  

$$\begin{array} {lll} \bold 0 &\approx& \bold Y \bold K \bold ... ...d r_0 \\ \bold 0 &\approx& \epsilon\ \bold R \bold a\end{array}$$ (56)

where $\bold R$ is a roughening operator.

Experience with missing data in Chapter shows that when the roughening operator $\bold R$ is a differential operator, the number of iterations can be large. We can speed the calculation immensely by ``preconditioning''. When we define a new variable $\bold m$ by $\bold a=\bold S\bold m$and insert it into (56) we get

$$\begin{eqnarray} \bold 0 &\approx & \bold Y \bold K \bold S\bold m + \bold r_0 \\ \bold 0 &\approx & \epsilon\ \bold R \bold S\bold m\end{eqnarray}$$ (57) (58)

The idea is to find $\bold m$.If we wanted to see $\bold a$ (normally we are only interested in $\bold r$)we could easily find it with $\bold a=\bold S\bold m$.

Before viewing the full program, we can reduce clutter in two ways. First, because the smoothing and roughening operators are somewhat arbitrary, we may as well replace (58) by the simpler $\bold 0\approx \epsilon\bold I\bold m$.The second way we can reduce clutter is to simply drop the damping (58) and keep only (57); then to control the null space, we need only to start from a zero solution and limit the number of iterations. As a practical matter, without (58) we must find a good number of iterations and with it we must find a good value for $\epsilon$.

The fitting (57) uses the operator $\bold Y\bold K\bold S$.For $\bold Y$ we can use subroutine nicaf1() , for the smoothing operator $\bold S$ we can use leaky integration with subroutine leakint() , and for the constraints that a₀=1 for all t, we can throw in a code fragment to zero those $\Delta a$and a few more for gapped filters. The result is subroutine tvdecop().  

module tvdec { # Linear operator for time-variable deconvolution
  use leakint2
  use nicaf1
  real, dimension( ny*na), allocatable :: aa
#% _init( yy, ny, na, rho)
    real, dimension( :), pointer :: yy
    integer, intent( in)         :: ny, na
    real,    intent( in)         :: rho 
    call nicaf1_init( yy, ny, na)
    call leakint2_init( ny, na, rho)
#% _lop( mod, dat)
    integer stat1, stat2
    if( adj) {
       stat2 = nicaf1_lop  ( adj, .false., aa, dat)
       stat1 = leakint2_lop( adj, .true.,  mod, aa)
    } else {
       stat1 = leakint2_lop( adj, .false., mod, aa)
       stat2 = nicaf1_lop  ( adj, .true.,  aa, dat)
    }
}

 

module leakint2 {
  use leakint
  integer :: n1, n2
#% _init (n1, n2, rho)
  real, intent (in) :: rho
  call leakint_init (rho)
#% _lop (x (n1,n2), y (n1,n2))
  integer :: i2, stat1
  do i2 = 1, n2
    stat1 = leakint_lop (adj, .true., x (:,i2), y (:, i2))
}

Finally, we turn to the initialization. Something new and slightly tricky arises here. It is an example of a fitting goal with preconditioners. If we were planning to start iterating from $\bold a=\bold 0$,there would be no problem because $\bold a=\bold S\bold m$tells us that we would be starting from $\bold m=\bold 0$.In fact, we must start from filters $\bold a_0$ that have a leading coefficient equal unity. The correct starting $\bold m_0$ must be the solution to the equation $\bold a_0=\bold S\bold m_0$.Luckily, in this case, we can easily confirm that the leaky integral of $\bold m_0(t) = (1, 1-\rho, 1-\rho, 1-\rho, 1-\rho, \cdots )$is $\bold a_0(t) = (1, 1, 1, 1, 1, \cdots )$.Otherwise, the filter solver subroutine tvdcon() is as uneventful as other simple solvers.  

module tv_decon { # Time Variable (TV) Deconvolution
  use tvdec
  use cgstep_mod
  use solver_mod
contains
  subroutine tvdecon( rho, gap, na, yy, rr, niter) {
    integer,                   intent( in)  :: gap, na, niter
    real,                      intent( in)  :: rho
    real,    dimension( :),    pointer      :: yy
    real,    dimension( :),    intent( out) :: rr
    real,    dimension( size( yy)*na)       :: mm                
    logical, dimension( size( yy)*na)       :: mask
    real,    dimension( size( rr))          :: dat
    integer                                 :: nt
    nt = size( yy) ; dat = 0.                         
    mm = 0. ; mm( 1) = 1. ; mm( 2:nt) = 1.- rho
    mask = .false.; mask( :gap*nt) = .true.
    call tvdec_init( yy, nt, na, rho)
    call solver( tvdec_lop, cgstep, x0 = mm, x = mm, dat = dat,
                 known = mask, niter = niter, res  = rr) 
    call cgstep_close()
    call tvdec_close()
  }
}

Figure 16 shows a synthetic data example using these programs. As we hope for deconvolution, events are compressed. The compression is fairly good, even though each event has a different spectrum. What is especially pleasing is that satisfactory results are obtained in truly small numbers of iterations (about three). The example is for two free filter coefficients and a gap of 6. Observe that the first 6 points of waveforms are preserved in the early iterations only.

 

tvdecon90
Figure 16 Time variable deconvolution with two free filter coefficients and a gap of 6. Four iterations is better than many.

I hope also to find a test case with field data, but experience in seismology is that spectral changes are slow, which implies unexciting results. Many interesting examples should exist in two- and three-dimensional filtering, however, because reflector dip is always changing and that changes the spatial spectrum.

EXERCISES:

1. Sketch the matrix corresponding to operator nicaf1 . HINTS: Do not try to write all the matrix elements. Instead draw short lines to indicate rows or columns. As a ``warm up'' consider a simpler case where one filter is used on the first half of the data and another filter for the other half. Then upgrade that solution from two to about ten filters.