Chebyshev Approach

For an alternative spectral approach, I adopted the Chebyshev-$\tau$method Gottlieb and Orszag (1977); Lanczos (1956). The Chebyshev-$\tau$ velocity continuation algorithm consists of the following steps:

1.

Transform the regular grid in t to Gauss-Lobato collocation points, required for the fast Chebyshev transform. First, a new variable $\xi$ is introduced by the shift transform:  

$$\xi = 1 - \frac{2\,t^2}{T^2}$$ (10)

so that the domain $0 \leq t \leq T$ is mapped into the domain $1 \geq \xi \geq -1$. Second, the $\xi$ grid points are distributed regularly in the cosine projection: $\xi_j = \cos(\frac{\pi j}{N}), j = 0,1,2,\ldots,N$.

2.

Transform the initial image P₀ (x,t) into the Chebyshev space in $\xi$ and Fourier transform in x, using the FFT algorithm. The Chebyshev-Fourier representation of P₀ (x,t) is  

$$P_0 (x,t) = \sum_{k=-N_x/2}^{N_x/2-1}\sum_{j=0}^{N_t} \hat{P}_{kj} T_j (\xi) e^{i k x}\;,$$ (11)

where T_j denotes the Chebyshev polynomial of degree j.

3.

Apply equation (1) to advance the image in velocity v. It is convenient to rewrite this equation in the form  

$${{\partial P} \over {\partial v}} = \frac{v T^2}{4}\,\int d\xi\, {{\partial^2 P} \over {\partial x^2}}\;.$$ (12)

In the Chebyshev-$\tau$ domain, the double differentiation in x is performed by multiplying the Fourier transform of P by -k², and integration in $\xi$ is performed as a direct operations on the Chebyshev coefficients. In particular, if $\sum_{j=0}^{N} a_j T_j (\xi)$ is the Chebyshev representation of the function $f (\xi)$, then the coefficients b_j of $\int f (\xi) d\xi$ are defined by the relation  

$$2 j\,b_j = c_{j-1} a_{j-1} - a_{j+1}$$ (13)

where c₀ = 2, c_j = 0 for j < 0, and c_j = 1 for j > 0. The constant of integration (and, correspondingly, the coefficient b₀) can be found at each velocity step from the boundary conditions (2), which are transformed to the form  

$$\left.P\right\vert _{\xi=-1} = \sum_{j=0}^{N_t} \hat{P}_{kj} (-1)^j = 0\;.$$ (14)

For the velocity advancement I used an implicit Crank-Nicolson scheme, which is unconditionally stable independent of the velocity step size. By writing equation (12) in the matrix form  

$$\frac{\partial \bold{P}}{\partial v} = \bold{A} \bold{P}\;,$$ (15)

the Crank-Nicolson advancement is represented by the equation  

$$\bold{P}_{v+dv} = \left(\bold{I} - \bold{A}\,\frac{dv}{2}\r... ... \left(\bold{I} + \bold{A}\,\frac{dv}{2}\right) \bold{P}_v\;,$$ (16)

where $\bold{I}$ is the identity matrix. The inverted matrix $\left(\bold{I} - \bold{A}\,\frac{dv}{2}\right)$ has a tridiagonal structure, except for the first row, implied by the boundary condition (14). A careful treatment of the boundary condition by the matrix-bordering method Boyd (1989); Faddeev and Faddeeva (1963) allows for an efficient inversion at a tridiagonal solver speed.

4.

Transform the result of the velocity advancement back to the physical domain.

5.

Transform the grid back to being regularly space in t.

 

cheb1
Figure 6 Synthetic seismic data before (left) and after (right) transformation to the Chebyshev grid in squared time.

 

cheb1-inv
Figure 7 The left plots show the reconstruction of the original data after transforming back from the Chebyshev grid to the original t grid. The right plots show the difference with the original model. Top: using the original grid size (N_t = 200). Bottom: increasing the grid size by a factor of three.

The first advantage of the Chebyshev approach comes from the better conditioning of the grid transform. Figure 6 shows the synthetic data before and after the grid transform. Figure 7 shows a reconstruction of the original data after transforming back from the Chebyshev grid (Gauss-Lobato collocation points). The difference with the original image is negligibly small.

 

cheb1-fft
Figure 8 Left: Synthetic data after Chebyshev transform. Right: the real part of the Fourier transform in the space coordinate.

The second advantage is the compactness of the Chebyshev representation. Figure 8 shows the data after the decomposition into Chebyshev polynomials in $\xi$ and Fourier transform in x. We observe a very rapid convergence of the Chebyshev representation: a relatively small number of polynomials suffices to represent the data.

 

cheb-impl
Figure 9 Impulse responses (Green's functions) of velocity continuation, computed by the Chebyshev-$\tau$ method. Top: without zero padding, bottom: with zero padding on the x axis. The left plots correspond to continuation to a larger velocity (+1 km/sec); the right plots, smaller velocity, (-1 km/sec).

The third advantage is the proper handling of the non-periodic boundary conditions. Figure 9 shows the velocity continuation impulse responses, computed by the Chebyshev method. As expected, no wraparound artifacts occur on the time axis, and the accuracy of the result is noticeably higher than in the case of finite differences (Figure 1).