THE PHASE SHIFT PLUS INTERPOLATION MODELING

The derivation of the Phase Shift Plus Interpolation (PSPI) modeling algorithm starts with the scalar wave equation  

$${\partial^2 p \over \partial z^2}+{\partial^2 p \over \partial x^2} = {1 \over v^2} {\partial^2 p \over \partial t^2}$$ (1)

where p=p(x,z,t) is the pressure field and v=v(x,z) is the earth velocity.

The pressure field p(x,z,t) is a finite function and can be therefore expressed as a double Fourier series  

$$p(x,z,t)={\sum_{k_x} \sum_{\omega} P(k_x,z,\omega) e^{i(k_x x -\omega t)}} .$$ (2)

Substituting equation (2) in equation (1) we obtain  

$${\sum_{k_x} \sum_{\omega}[ {\partial^2 P(k_x,z,\omega) \over... ...a^2 \over v^2(x,z)} P(k_x,z,\omega)] e^{i(k_x x-\omega t)}= 0}.$$ (3)

Equation (3) should hold for any values of k_x and $\omega$.This is possible only if each term inside the square parenthesis is zero. This reasoning is similar to the condition that if a polynomial is zero for any values of x, the coefficients of the polynomial are zero. Equation (3) becomes  

$${ {\partial^2 P(k_x,z,\omega)} \over {\partial z^2}}={(k^2_x- {\omega^2 \over {v^2(x,z)}})} P(k_x,z,\omega)$$ (4)

valid for all values of k_x and $\omega$.The problem is that in this form, the x-coordinate in the pressure field is Fourier transformed and there is no direct correspondence between a point (x,z) in the medium, the velocity v(x,z), and the corresponding value of p(x,z,t) at that location.

For a constant velocity we write  

$$k_z={\pm {[{\omega^2 \over v^2}-{k^2_x }]^{1 \over 2}}}$$ (5)

where k_z is constant for two given values of k_x and $\omega$.Equation (5) is the well known dispersion relation. Equation (4) becomes an ordinary differential equation  

$${\partial^2 P \over \partial z^2}=-k_z^2 P .$$ (6)

For a constant k_z equation (6) has the analytic solution  

$$P(k_x,z_0+z,\omega)={P(k_x,z_0,\omega)}e^{ik_z z} .$$ (7)

Equation (7) shows how in constant velocity media we can find the pressure field at a certain depth level $P(k_x,z_0+z,\omega)$, if we know it at any other depth level $P(k_x,z_0,\omega)$.

However there are several ambiguities in equation (7) that we need to discuss. One is the question of time propagation. If we know the pressure field (or wavefield) at a certain depth we can propagate it forward in time or backward in time. We can also propagate it up (along the z-axis) or down. To understand how we determine the propagation direction we have to analyze the values and sign of k_z. There are several restrictions on the values of k_z. Equation (6) has the solution (7) only for real values of k_z which imposes the condition

$${\omega^2 \over v^2}-{k^2_x } \geq 0 .$$

The solution represented in equation (7) is the Fourier transform of the wavefield. The general solution in time-space coordinates is obtained by summing all the Fourier coefficients obtained from equation (7)  

$$p(x,z,t)={\sum_{k_x} \sum_{\omega} P(k_x,z_0,\omega) e^{ik_z z} e^{i(k_x x -\omega t)}} .$$ (8)

The function

$$e^{i(k_z z + k_x x -\omega t)}$$

represents a plane wave. Equation (8) sums many plane waves to obtain the general solution. We examine the sign of k_z necessary to upward propagate a wave by examining each plane wave solution. If we ignore k_x x, which determines the lateral variation, we can write

$$k_z z - \omega t = phase(z,t) .$$

The phase is constant along a plane wave, and we write

$$k_z z = { \omega t + const}$$

for the phase of a particular plane wave. The plane wave is moving downward when k_z has the same sign with $\omega$ because z increases with t in order to keep the phase constant. So for the upward moving waves we need to have opposite signs of k_z and $\omega$(z is decreasing when t is increasing). We have now figured out that in order to have only upgoing waves we have to look at the sign of $\omega$ and assign to k_z the opposite sign.

Next we correlate this information with the depth level where we want to find the wavefield. If the known wavefield is at depth z₀ and we want to find the wavefield at depth z₀+z, then we have to propagate the wavefield back in time (toward t=0) because we know that the wavefield travels upward. If the known wavefield is at depth z₀ and we want to find the wavefield at depth z₀-z, then we are propagating the wavefield forward in time. This is the direction we are interested in for modeling.

However for depth varying velocity v(z) we have k_z approximately constant only for small depth intervals ($\Delta z$) where we can consider the velocity constant. Therefore equation (7) becomes  

$$P(k_x,z_0+\Delta z,\omega)={{P(k_x,z_0,\omega )}e^{ik_z \Delta z}}$$ (9)

and can be used to downward or upward extrapolate the wave field for a small depth interval.

For laterally variant media, Gazdag and Sguazzero (1984) propose to downward extrapolate the wavefield one depth interval at a time with several velocities. We can apply the same idea to upward propagate the wavefield with several velocities. We consider several velocities (v₁, v₂, ...) in the interval [v_min,v_max] and upward propagate the wavefield $P(k_x,z_0,\omega)$ to $P(k_x,z_0-\Delta z,\omega)$ with each velocity. We can afterward inverse Fourier transform in x the resulting wavefields

$$(P_1(k_x,z_0-\Delta z,\omega),P_2(k_x,z_0-\Delta z,\omega),...)$$

to obtain

$$(P_1^*(x,z_0-\Delta z,\omega),P_2^*(x,z_0-\Delta z,\omega),...) .$$

To obtain a single upward propagated wavefield, at each point $(x, z_0-\Delta z)$associated with a velocity $v(x,z_0-\Delta z)$, the value of the resulting wavefield in this point is interpolated between the two wavefields with closest velocities ($v_i(x,z_0-\Delta z)$).

In addition we need to implement a technique (Gazdag and Sguazzero, 1984) in the PSPI modeling algorithm to ensure that all the zero dips (corresponding to the case k_x=0) are upward propagated without distortion. The technique consists of multiplying the wavefield $P(x,z_0,\omega)$ with  

$$e^{-i{\omega \over v_j}\Delta z}$$ (10)

after the Fourier transformation in x and multiplying the upward propagated wavefield $P_j(k_x,z-\Delta z,\omega)$ by  

$$e^{i{\omega \over v(x,z)} \Delta z}$$ (11)

where the subscript j denotes the index of the constant velocity used in the upward propagation step. For a zero dip plane wave (i.e. k_x = 0) we have  

$$k_z = \pm {\omega \over v }$$ (12)

and the plane wave can be upward propagated with the velocity v(x,z) corresponding to the point coordinates (x,z) instead of the constant velocity v_j.

As it will be seen later when the PSPI results are compared to the Split-step Fourier results, this supplemental phase subtraction and addition improves the accuracy in the migration case but actually produces worse results in the modeling case. Figure shows the flow of the PSPI migration and modeling algorithms.

 

pspiboth
Figure 1 Phase Shift Plus Interpolation (PSPI) migration and the conjugate transpose PSPI modeling algorithm.