Substituting this wavefield in equation (1), the prestack migration equation can be written as
FFT along t,y axes forall do z do ky do do kh
The loop in kh 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
FFT along t,y axes do z do ky do compute phase fast
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.