Migration for a single constant-offset section

Consider a single 2-D constant-offset section p(t,y;h₀) corresponding to a half-offset h₀ in a 3-D cube where all the other constant-offset sections are zero. This is equivalent to introducing a Dirac $\delta$ function in offset coordinates. The Fourier transform over offset of such a 3-D field is

$$P(\omega,k_y,k_h)=e^{-ik_h h_0}p(\omega,k_y;h_0).$$

Substituting this wavefield in equation (1), the prestack migration equation can be written as  

$$\begin{array} {lcl} p(t=0,k_y,h=0,z) & = & \displaystyle{ {\... ...) \int dk_h \; e^{ik_z(\omega,k_y,k_h)z-ik_h h_0}}.\end{array}$$ (5)

The DSR algorithm can be rewritten, extracting the multiplication with the wavefield $P(\omega,k_y;h_0)$ outside the k_h loop.

		FFT along t,y axes $p(t,y;h_0) \rightarrow P(\omega,k_y)$ 
		forall $k_y,\omega,k_h;$ $ph(\omega,k_y,k_h)=e^{-ik_h h_0}$ 
		do z 
		 		do ky 
		 		do $\omega$ 
		 		 		do kh 
		 		 		$ph(\omega,k_y,k_h)=ph(\omega,k_y,k_h)e^{ik_z dz}$ 
		 		 		$phase(\omega,k_y)=phase(\omega,k_y)+ph(\omega,k_y,k_h)$ 
		 		$M(z,k_y)=M(z,k_y)+P(\omega,k_y)*phase(\omega,k_y)$ 

The loop in k_h becomes a mere numerical integration of the exponential term, integral which can asymptotically approximated through the use of the stationary phase (Appendix A). For depth variable velocity, we have to sum the phase terms corresponding to each depth level. In other words, for depth variable velocity the exponential term contains the sum of all the phases corresponding to previous depth levels. Just knowing the stationary phase for a depth level is not enough to carry the computation to the next depth level.

This is a terribly slow prestack phase-shift migration algorithm that migrates one constant-offset section in almost the time necessary to migrate all constant-offset sections. The algorithm can be sped up hundreds of times if a good stationary phase approximation can be found to the integral  

$$\int dk_h e^{ik_z z- ik_h h_0}.$$ (6)

Such an algorithm would become:

		FFT along t,y axes $p(t,y;h_0) \rightarrow P(\omega,k_y)$ 
		do z 
		 		do ky 
		 		do $\omega$ 
		 		compute phase fast
		 		$M(z,k_y)=M(z,k_y)+P(\omega,k_y)*phase$ 

In Appendix A, I discuss different properties of the phase function and how the stationary phase could be computed numerically. I could not find an analytical solution for the stationary phase, as it involves finding the roots of a polynomial of sixth degree. However, even if an analytical solution is found for constant velocity, for depth variable velocity the solution incorporates the sum of phases corresponding to different depth slices.