Prestack Split-Step Migration

Figure B.1: Flowcharts for prestack split-step modeling and migration. $\mathcal{F}$ , $\mathcal{F}^{-1}$ , $\mathcal{P}$ , and $\mathcal{C}$ denote the Fourier transform, the inverse Fourier transform, the phase shift operator, and the phase correction operator, respectively (see text for details). (a) The source wavefield $P(x,z,\omega)$ is propagated from the surface of the earth to depth $z$ in steps $\Delta z$ . (b) At each depth $z$ of the earth, a reflected wave is generated by $m(x,z)P(x,z,\omega)$ , where $m(x,z)$ is a reflectivity model. The wave is then propagated upward to the surface $z=0$ . The total reflected wavefield $R(x,z,\omega)$ consists of the superposition of the reflected and propagated waves originating from below. The total reflected wavefield collected at the surface is the data, i.e., $R(x,z=0,\omega) \equiv D(x,z=0,\omega)$ . (c) The data are then downward continued from the surface back to depth $z$ in steps $\Delta z$ . Finally, the migration image $I(x,z)$ (not shown) is constructed by applying the imaging condition: $I(x,z) = \sum_{\omega}P^*(x,z,\omega) D(x,z,\omega)$ , or $I(x,z) = 2\sum_{\omega>0} \textrm{Re} \{ P^*(x,z,\omega) D(x,z,\omega) \}$ , assuming the DC component is 0.

The steps of prestack split-step migration are described in the flowcharts shown in Figure B.1. This presentation closely follows Kuehl and Sacchi (1999), and is included here for convenience because analysis of computational complexity refers to it.

Consider first the forward propagation of a wavefield. The split-step operator $\mathcal{L}$ per layer can be decomposed into a succession of four linear operators $\mathcal{F}$ , $\mathcal{P}$ , $\mathcal{F}^{-1}$ , and $\mathcal{C}$ . First, the seismic wavefield $P(x,z,\omega)$ at $z$ is transformed to the wavenumber $k_x$ domain by the Fourier operator $\mathcal{F}$ . Second, the phase-shift operator $\mathcal{P}$ is applied to the wavefield in the $(k_x,\omega)$ domain:

$$P_1(k_x,z,\omega) = P(k_x,z,\omega)e^{-i\Delta z\sqrt{(\omega u_0)^2 - k_x^2}},$$ (7.1)

where $u_0$ is the mean slowness for the current layer. Third, $P_1(k_x,z,\omega)$ is transformed to the space $x$ domain by the inverse Fourier operator $\mathcal{F}^{-1}$ . Fourth, the phase correction operator $\mathcal{C}$ is applied in the $(x,\omega)$ domain. This accounts for the lateral slowness variation $\Delta u(x) = u(x) - u_0$ :

$$P(x,z+\Delta z,\omega) = P_1(x,z,\omega) e^{-i\omega\Delta z\Delta u(x)}.$$ (7.2)

Altogether, it is given that

$$\mathcal{L} = \mathcal{C}\mathcal{F}^{-1}\mathcal{P}\mathcal{F},$$ (7.3)

of which the adjoint is

$$\mathcal{L}^\dagger = \mathcal{F}^\dagger \mathcal{P}^\dagger {\mathcal{F}^{-1}}^\dagger \mathcal{C}^\dagger$$

$$= \mathcal{F}^{-1} \mathcal{P}^* \mathcal{F} \mathcal{C}^*,$$ (7.4)

The adjoint operator $\mathcal{L}^\dagger$ applies to the case of `backward propagation', or downward continuation of the data, as illustrated in Figure B.1(c). This ensures that the migration operator is the adjoint of the forward modeling counterpart.