Computation of the gradient vector

The gradient vector defined in equation (5) is the sum of traveltime residual term ${\bf g}^{(1)}$ and condition term ${\bf g}^{(2)}$. The computation of the condition term is straightforward because matrix ${\bf B}$ is given. In principle, we can compute the traveltime residual term by first computing matrix ${\bf A}$, and then performing the matrix multiplication. However, matrix ${\bf A}$ has $N\times M$ elements that are functions of r; hence, both its calculation and storage are of high cost. If we can compute the vector ${\bf g}^{(1)}$ directly using a finite difference method, we can considerably reduce the cost in computation time and storage.

If we solve the first-order linear PDE (9) using the method of characteristics, we find that

$$a_{ij}(r)=\int^{s_i(r)}_0 \beta_j(\hat{x}_i(\xi,r),\hat{z}_i(\xi,r)) d\xi,$$ (10)

where $\hat{x}_i(s,r)$ and $\hat{z}_i(s,r)$ are the coordinates of the ray that starts from the number i source and ends at the receiver at depth r, and s is the arc-length of the ray. Clearly, the elements of matrix ${\bf A}$ are the integrations of the basis functions along rays. The component of the gradient vector can be computed as follows:  

$$\begin{array} {lll} g^{(1)}_j & = & \displaystyle{\int_{r}\s... ...int_z \beta_j(x,z) \delta t_i(\hat{r}_i(x,z))}Jdxdz,\end{array}$$ (11)

where J is the Jacobian of the coordinate transformation:

$$J= \left\vert \begin{array} {ll} \displaystyle{\partial s \o... ...al x}\right)^2+ \left({\partial r \over \partial z}\right)^2}.$$ (12)

The last step of equation (11) is an area integration that is easy to implement. Appendix A shows that function $r=\hat{r}_i(x,z)$ satisfies the first-order linear PDE, as follows:  

$${\partial \tau_i \over \partial x}{\partial \hat{r}_i \over ... ...tau_i \over \partial z}{\partial \hat{r}_i \over \partial z}=0,$$ (13)

with the initial condition that function $\hat{r}_i(x,z)$ is equal to the receiver depth at each receiver location. We can use a finite difference algorithm, similar to the one used for traveltime calculation, to solve equation (13). The algorithm should start from receiver locations, and extrapolate function $\hat{r}_i(x,z)$ in the opposite directions of the wave propagations, towards the source.