Implementation

The frequency-wavenumber log-stretch PS-AMO operator works on a regularly sampled cube. However, the data is usually recorded on an irregular mesh. To overcome this obstacle, I will first map the data that are on an irregular mesh to a regular 5-D grid with dimensions ($t, {\rm CMP}_x, {\rm CMP}_y,h_x,h_y$). For migration efficiency I decrease the dimensionality of the dataset, I reduce the data dimensions by creating a common-azimuth cube oriented along the inline direction, that is to eliminate the crossline offset (h_y) axis.

I use the nearest-neighbor interpolation operator ($\bf L'$)to map the data from an irregular mesh into a regular mesh. Then, I use the PS-AMO operator to transform data from non-zero crossline offsets ($h_y \neq 0$) to zero crossline offset (h_y=0), I refer to this operation as an operator $\bf Z'$which is a summation over h_y. I allow some mixing between h_x by expanding our summation to form h_x=a and h_y=0 as following:  

$$\sum_{a/2 - \Delta h_x}^{a/2 + \Delta h_x} \sum_{h_y},$$ (46)

where $\Delta h_x$ is small.

I combine these two operators to estimate a 4-D model ($\bf m$) from a 5-D irregular dataset ($\bf d$) through the objective function,  

$$Q(\bf m) = \Vert\bf d- \bf L\bf Z\bf m\Vert^2.$$ (47)

In the context of least-squares inversion a model regularization term can be included. This regularization could be a derivative operator ($\bf {G_h}$)applied to data cubes that have been transformed to the same offset using the PS-AMO operator. This regularization assumes that the reflection amplitudes are a smooth function of the reflection angle and azimuth. Therefore, this least-squres problem (equation ) becomes:

$$Q(\bf m) = \Vert\bf d- \bf L\bf Z\bf m\Vert^2 + \Vert\epsilon_G {\bf {G_h}} \bf m\Vert^2.$$ (48)

The adjoint solution for this inverse problem is  

$${\bf m}_{L2} \approx {\bf m}_{\rm adj} = \bf Z' \bf L' \bf d.$$ (49)

However, the adjoint solution is not ideal. The irregularity of our data can lead to artificial amplitude artifacts. A more accurate solution to this inverse problem is the weighted adjoint solution. This solution uses a diagonal operator, that is a model-space weighting function, to approximate the Hessian for the least-squares problem  Claerbout and Nichols (1994); Rickett (2001), such that  

$${\bf m}_{L2} \approx \bf m_{\rm wt} = {\bf W_{m}}^{-1} \bf Z' \bf L' \bf d$$ (50)

where the diagonal operator ${\bf W}_{\rm m}$ is an approximation to the Hessian  

$$\bf H^{-1} \approx {\bf W_{m}}^{-1} \equiv \rm diag \left[ \... ...f L\bf Zm_{{\rm ref}}}{{\bf m}_{{\rm ref}}+ \epsilon } \right],$$ (51)

where $\epsilon$ is a tiny value to avoid dividing by zero.

The solution of problem  is not feasible on a single computer. The computational requirements are onerous, but potentially manageable. However, the memory requirements are not. A full regularized 5-dimensional cube, that I create after applying $\bf L'$,can easily reach tens of gigabytes. This size of data makes it almost impossible to practically implement any algorithm for 3-D prestack seismic data-processing on a single machine.

Clapp (2004) introduces an efficient python library for handling parallel jobs. The library makes it easy for the user to take an already existing serial code and transform it into a parallel code. The library handles distribution, collection, and node monitoring, commonly onerous tasks in parallel processing.

The main prerequisite to using the python library is to build an efficient serial code, and to describe how the parallel job should be distributed on a cluster. For this problem I chose to split along the h_x axis. I created a series of tasks, each assigned to produce a single ($t,{\rm CMP}_x, {\rm CMP}_y$) volume. Each task is passed a range of h_x's defined by equation . The resulting model volumes are then recombined to form the regularized 4-D output space.