Theory

Automatic timeshift estimation algorithms maximize the similarity between baseline and monitor seismic traces (Hale, 2006; Fomel and Jin, 2009). Alternatively, one can minimize the residual between baseline and shifted monitor seismic traces (Phillips and Fomel, 2016; Hale, 2009). This condition can be represented in the $Z$-transform notation as

$$D(a,p)=1-aZ_4W(p) \ ,$$ (1)

where $W$ is a warping filter, $Z_4$ is a phase shift between time-lapse seismic images, $a$ is an amplitude weight, and $p$ is the 4D timeshift. This assumption fails when layers are below seismic resolution and significant interference between reflection sidelobes exists. We partially alleviate the problem of interference by decomposing seismic images into discrete frequency components using the local time-frequency transform (Liu and Fomel, 2013).

The local-time-frequency transform is based on the idea of non-stationary regression. A digital signal $f(t)$ can be represented as Fourier series

$$f(t)=\sum_{n=-\infty}^{\infty}b_n(t)\psi_n(t)$$ (2)

where the Fourier coefficients ${b}_n$ are allowed to vary temporally and $\psi_n$ is a vector consisting of complex exponentials at each corresponding Fourier frequency. The Fourier coefficients $b_n$ are estimated by regularized least-squares inversion (Fomel, 2008).

We then use amplitude-adjusted plane-wave destruction filters (Phillips and Fomel, 2016) to estimate timeshifts at each frequency. In the linear operator notation, equation (1) is modified to

$$\mathbf{D}(\mathbf{a},\mathbf{p})\mathbf{d}=\mathbf{B}_l(\mathbf{p})\mathbf{d}- (\mathbf{a})\mathbf{B}_r(\mathbf{p})\mathbf{d}$$ (3)

where $\mathbf{B}_l$ and $\mathbf{B}_r$ are the left- and right-hand side of the plane-wave destruction filter (Fomel, 2002), as described by Phillips and Fomel (2016), and $\mathbf{d}$ is the local time-frequency transform of the time-lapse seismic data. Our objective is to minimize the plane-wave residual between the time-lapse seismic images at each frequency ( $\mathbf{D}(\mathbf{a},\mathbf{p})\mathbf{d}\approx\mathbf{0}$).

The dependence of $\mathbf{D}$ on $\mathbf{a}$ is linear; however, $\mathbf{p}$ enters in a non-linear way (Fomel, 2002). We separate this problem into linear and non-linear parts using the variable projection technique (Golub and Pereyra, 1973; Kaufman, 1975). The algorithm is described below (Phillips and Fomel, 2016):

1. Set $\mathbf{p}_0 = \mathbf{0}$ and $\mathbf{a}_0 = \mathbf{1}$.
2. Hold the scale $\mathbf{a}$ constant and compute the shift $\mathbf{p}_n$ using accelerated plane-wave destruction (Chen et al., 2013).
3. Shift the monitor image using $\mathbf{p}_n$
4. Hold the shift $\mathbf{p}_n$ constant and compute the amplitude ratio $\mathbf{a}_n$ by the smooth division of the left and right side of the plane-wave destruction filter $\mathbf{D}$ in equation(3):

   $$\mathbf{a}_n = \left<\frac{\mathbf{B}_l(\mathbf{p}_n)\mathbf{d}}{\mathbf{B}_r(\mathbf{p}_n)\mathbf{d}}\right>$$ (4)

5. Scale the monitor image using $\mathbf{a}_n$
6. Iterate until convergence (return to step 2).

If timeshifts are large ( $\max \vert\mathbf{p}_n\vert\geq$ 10 samples), rather than setting $\mathbf{p}_0 = \mathbf{0}$, it may be necessary to instead provide a low frequency estimate of the timeshift calculated using another algorithm, such as local similarity (Fomel and Jin, 2009). So long as this “small timeshift" condition is satisfied, this algorithm provides an improved estimate of true 4D timeshifts from spectrally decomposed time-lapse seismic images.