previous up next print clean


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}\end{displaymath} (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)}} .\end{displaymath} (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}.\end{displaymath} (3)
Equation (3) should hold for any values of kx 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)\end{displaymath} (4)
valid for all values of kx 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}}}\end{displaymath} (5)
where kz is constant for two given values of kx 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 .\end{displaymath} (6)
For a constant kz equation (6) has the analytic solution  
P(k_x,z_0+z,\omega)={P(k_x,z_0,\omega)}e^{ik_z z} .\end{displaymath} (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 kz. There are several restrictions on the values of kz. Equation (6) has the solution (7) only for real values of kz which imposes the condition

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

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)}} .\end{displaymath} (8)
The function

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

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

k_z z - \omega t = phase(z,t) .\end{displaymath}

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

k_z z = { \omega t + const}\end{displaymath}

for the phase of a particular plane wave. The plane wave is moving downward when kz 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 kz 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 kz 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 z0 and we want to find the wavefield at depth z0+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 z0 and we want to find the wavefield at depth z0-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 kz 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}}\end{displaymath} (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 (v1, v2, ...) in the interval [vmin,vmax] 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),...)\end{displaymath}

to obtain

(P_1^*(x,z_0-\Delta z,\omega),P_2^*(x,z_0-\Delta z,\omega),...) .\end{displaymath}

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 kx=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}\end{displaymath} (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}\end{displaymath} (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. kx = 0) we have  
k_z = \pm {\omega \over v }\end{displaymath} (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 vj.

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.

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

previous up next print clean
Stanford Exploration Project