LEAST SQUARES INVERSE SLANT STACKS

The amplitude spectra displaying in the previous spectrum were created using a ``direct'' slant stack. The transformation to the slant stack domain was performed using the traditional ``sum along lines'' operator plus a rho filter.

$$U(p,\tau) = \rho(\tau) * \sum_{i=1}^{nx} U(x_i,t=\tau-px)$$

Where the rho filter is defined as $\rho(\omega)=\vert\omega\vert$.

The problems with using this approach have been well documented Kostov (1990). The limited aperture and sampling in the x-t domain cause artifacts. Irregular sampling in x is not well handled by the direct stacking procedure. The alternative is to perform a least squares inverse slant stack. This method assumes that we can exactly perform the mapping from $p-\tau$ to x-t and casts the mapping from x-t to $p-\tau$ as an inverse problem. This is most conveniently calculated in the frequency domain.

The operator that maps from $p-\omega$ to $x-\omega$ is L,

$$L(p_i,x_j,\omega) = e^{-i \omega p_i x_j }$$

$$U(x,\omega) = L\, U(p,\omega)$$

Given the data in the x-t domain we must construct a mapping to the $p-\tau$ domain that ``best'' reconstructs the x-t data. The misfit function for the reconstruction is defined as a least-squares error.

$$\vert\vert L\, U(p,\omega)\ -\ U(x,\omega) \vert\vert^2$$

Minimization of this misfit is equivalent to solving the inverse problem, $L\,U(p,\omega)=U(x,\omega)$. The solution for the overdetermined problem is,

$$U(p,\omega) = (L^HL)^{-1}L^HU(x,\omega),$$

and for the underdetermined problem it is,

$$U(p,\omega) = L^H(LL^H)^{-1}U(x,\omega).$$

Where L^H is the conjugate transpose of L.

The operator (L^HL) has a Toeplitz form when the data is regularly sampled in p. This makes the inversion cheap to solve using Levinson recursion. A problem with this method is that the operator L^HL becomes singular when the data is aliased. Kostov suggests two ways to overcome this problem, ``Aliasing can be overcome by using prior information, either in the form of an interval of wavenumbers smaller than twice the Nyquist wavenumber, or by introducing ``a priori'' information via a model-covariance matrix.'' (I would suggest that the first statement can be modified to allow multiple wavenumber intervals as long as the total width does not exceed twice the Nyquist wavenumber and no two intervals occur at positions that are the alias of one another.)

My first solution to the problem is to constrain the solution to be non-zero only in those places that I have defined as unaliased. I define a masking operator, M, that is one where the data is not aliased and zero where it is aliased. I then solve the least squares problem,

$$min \vert\vert L\,M\,U(p,\omega)\ -\ U(x,\omega) \vert\vert^2.$$

The equivalent inverse problem with the solutions,

$$U(p,\omega) = (M^HL^HLM)^{-1}M^HL^HU(x,\omega),$$

and

$$U(p,\omega) = M^HL^H(LMM^HL^H)^{-1} U(x,\omega),$$

does not have the same Toeplitz structure as the unmasked problem so a different solution method must be used. I use a simple conjugate gradient algorithm Hestenes and Steifel (1952) to solve the minimization problem. This problem is solved separately for each frequency so the size of each problem is only $np\times nx$ and the convergence is rapid. The masking accelerates convergence at high frequencies because the problem is now better determined, there is no ambiguity between aliased and unaliased dips at each frequency. Figure  shows the $p-\omega$ spectrum of the least squares inverse slant stack using the masking function. All the energy is concentrated at the correct dip and the amplitude spectrum is flat. Figure  has three frames. The top frame is the original x-t domain data. The middle frame is the reconstruction of the x-t data from the masked least squares inverse slant stack. The bottom frame is the difference between the two panels plotted at the same scale. The reconstruction has been completely successful. I am not limited to reconstructing the data at the original trace spacing, the data can be modeled in the x-t domain at an arbitrary regular or irregular trace spacing.

 

single-lsq Figure 5 $p-\omega$ spectrum of the least squares transform using the masked operator.

 

single-lsqrec
Figure 6 Top: original data. Middle: reconstructed data using the masked least squares operator. Bottom: difference between the two plots.