The basic idea is similar to phase unwrapping Claerbout (1999), but instead of summing phase differences to get total phase, dips are summed to get total time shifts that are then used to flatten the data. To apply the shifts, a reference is held constant and all other traces are shifted vertically to match it.
The first step is to calculate dips everywhere in the 3-D seismic cube. Dip can be efficiently calculated using a plane-wave destructor as described in Claerbout (1992) or with an improved dip estimator that is described in Fomel (2002). We primarily use the latter. For each point in the data cube, two components of dip, px and py, are estimated in the x direction and y direction, respectively. These can be represented everywhere on the mesh as vectors and .
Our goal is to find a time-shift (or depth-shift) field such that its gradient approximates the dip . Our objective function is:
We apply the Euler-Lagrange equation to equation (1) to find:
We solve this using a Gauss-Newton approach by iterating over equations (6)-(8)
To dramatically improve efficiency, we solve equation (7) in the Fourier domain. We apply the divergence to the estimated dips and divide by the z-transform of the Laplacian in the Fourier domain with:
Once we have converged the resulting time-shift field, contains all of the time-shifts required to map the original unflattened data to flattened data. This is implemented by applying the time-shifts relative to a reference trace. In other words, each trace is shifted to match the reference trace.
In general, we operate on a one-time slice at a time. After iterating until convergence, we then select the next time slice and proceed down through the cube. In this way, each time slice is solved independently.
The process of flattening tends to alter the spectrum of the data by stretching and compressing in time. Even worse, it can disrupt its continuity. To insure a monotonic and continuous result, it should be sufficient to first smooth the input dips in time (or depth). In some instances, it may be necessary to enforce smoothness during the integration of the dips. This can easily be accomplished by defining a 3D gradient operator with an adjustable smoothing parameter as: