Horizontally Uniform Structures

The 2N_x by 2N_x matrix A can be partioned according to

$$[A] =\pmatrix{0 & I_{N_x}\cr A_{21} & 0}$$ (5)

where I_Nx denotes the N_x by N_x identity matrix and the N_x by N_x submatrix A₂₁ is given by

$$[A_{21}] =\pmatrix{-{\omega^2 \over c^2}-{\partial^2 \over \p... ... \ldots &-{\omega^2\over c^2}-{\partial^2 \over \partial x^2} }$$ (6)

The horizontal discretization requires an approximation for the derivative term in equation (6). As with the ordinary phase-shift method (Gazdag, 1978), this derivative can be calculated by Fourier method

$${\partial^2 P\over \partial x^2} = -k_x^2 P$$ (7)

Since multiplication in the space domain is equivalent to convolution in the wavenumber domain, the multiplication of the sloth (squared slowness) becomes convolution. For horizontally uniform structures, however, the Fourier transform of the velocity function is a delta function. Therefore, the submatrix A₂₁ has the form

$$[A_{21}] =\pmatrix{-{\omega^2 \over c^2}+k_x^2 & 0 & \ldots &... ...dots & \vdots \cr 0 & 0 & \ldots &-{\omega^2\over c^2}+k_x^2 }$$ (8)

The diagonal elements of this matrix are the eigenvalues of A₂₁ and a matrix in which the columns consist of the eigenvectors form an identity matrix. By letting the eigenvalues of submatrix A₂₁ as

$$\lambda_i=-{\omega^2 \over c^2}+k^2_{x_i},$$

the eigenvalues of matrix A become $\pm \sqrt{\lambda _i}$, for i=1,...,N_x-1. The corresponding eigenvectors are given by

$$\pmatrix{I_0 \cr \sqrt{\lambda_0}I_0},\pmatrix{I_1 \cr \sqrt... ...ldots,\pmatrix{I_{N_x-1} \cr -\sqrt{\lambda_{N_x-1}}I_{N_x-1}}.$$

where I_i denotes i-th column of identity matrix. With these eigenvectors we define a 2N_x by 2N_x matrix Q in which the columns consist of the eigenvectors of A. We then have the relation

$$\pmatrix{P\cr {\partial P\over \partial z}}_{z+dz} = \exp[Ad... ...mbda dz}Q^{-1} \pmatrix{P\cr {\partial P\over \partial z}}_{z}$$ (9)

where $\Lambda$ is the diagonal matrix

$$\Lambda = \pmatrix { \sqrt{\lambda_0} & 0 &\ldots & 0 & 0 & ... ...\cr 0 & 0 & \ldots & 0 & 0 & \ldots & -\sqrt{\lambda_{N_x-1}}}$$ (10)

and

$$Q^{-1}={1\over 2}\pmatrix { I_0^T &\ldots &{I_0^T/ \sqrt{\la... ...I_{N_x-1}^T &\ldots &-{I_{N_x-1}^T/ \sqrt{\lambda_{N_x-1}}}\cr}$$ (11)

When both sides of equation (9) are multiplied by Q^-1 from the left we obtain

$$Q^{-1}\pmatrix{P\cr {\partial P\over \partial z}}_{z+dz} = \... ...bda dz}]Q^{-1} \pmatrix{P\cr {\partial P\over \partial z}}_{z}$$ (12)

The equation (12) is a matrix representation for a set of equations

$$(P_i + {1 \over {\sqrt {\lambda_i}}}{\partial P_i \over \par... ...{1 \over {\sqrt{ \lambda_i}}}{\partial P_i \over \partial z})_z$$ (13)

and

$$(P_i - {1 \over {\sqrt {\lambda_i}}}{\partial P_i \over \par... ...{1 \over {\sqrt {\lambda_i}}}{\partial P_i \over \partial z})_z$$ (14)

The solutions of equations (13) and (14) represent upgoing and downgoing waves respectively. When an eigenvalue $\pm \sqrt{\lambda _i}$ is real, the component is propagated by a multiplication by a real exponential, and results in an evanescent component which should be removed for numerical stability. When an eigenvalue $\pm \sqrt{\lambda _i}$ is imaginary, the corresponding coefficient is propagated by phase shift.