# Teleseismic Tomography for Imaging the Upper Mantle Beneth Northeast China
## 2 Methodology
### Subspace Method  
In many tomographic methods, there are some one-dimensional subspace methods, such as the steepest descent method and conjugate gradient method are all one-dimensional subspace methods, and these methods only perform a line minimization at each iteration. Therefore, we need to consider constructing a subspace method, which can simultaneously perform the minimization along the search direction that together span a subspace of the model space, as shown in Figure 1. According to the above reasons, the subspace iteration method is selected for inversion in this paper. The subspace method limits the minimization of the quadratic approximation of the objective function \(S(\mathrm{m})\) to an n-dimensional subspace of the model space, so that the perturbation \(\delta\mathrm{m}\) occurs in the space spanned by a set of n M-dimensional basis vectors \(\left\{a^{j}\right\}\):  
\[\delta\mathrm{m}=\sum_{j=1}^{n}\mu_{j}a^{j}=\mathrm{A}\mu \tag{7}\]  
\(\mu_{j}\) is the component which determines the length of the corresponding vector \(a^{j}\) that minimizes the quadratic form of \(S(\mathrm{m})\) in the space spanned. A is the M \(\times n\) projection matrix. If the objective function is smooth enough, the following formula can be obtained:  
\[S(\mathrm{m}+\delta\mathrm{m})\approx S(\mathrm{m})+\overset{\wedge}{\gamma} \delta\mathrm{m}+\frac{1}{2}\delta\mathrm{m}^{T}\overset{\wedge}{\mathrm{H}} \delta\mathrm{m} \tag{8}\]  
where \(\overset{\wedge}{\gamma}=\partial S/\partial\mathrm{m}\) and \(\overset{\wedge}{\mathrm{H}}=\partial^{2}S/\partial\mathrm{m}^{2}\) are gradient vector and Hessian matrix, respectively. Their expression can be expressed as:  
\[\overset{\wedge}{\gamma}=\mathrm{G}^{T}C_{d}^{-1}[g(\mathrm{m})-\mathrm{d}_{ obs}]+\varepsilon C_{m}^{-1}(\mathrm{m}-\mathrm{m}_{0}) \tag{9}\]  
\[\overset{\wedge}{\mathrm{H}}=\mathrm{G}^{T}C_{d}^{-1}\mathrm{G}+\mathrm{V}_{m} \mathrm{G}^{T}\mathrm{G}_{d}^{-1}[g(\mathrm{m})-\mathrm{d}_{obs}]+\varepsilon C _{m}^{-1} \tag{10}\]  
Figure 1: A contour plot of S(m) which is a function of two parameters of different physical dimensions [14]. where \(\mathrm{G=}\partial g/\partial\mathrm{m}\) is the Frechet derivative matrix of partial derivatives which is calculated during the solution of the forward modeling. \(\mathrm{C_{d}}\) is a data covariance matrix. \(\mathrm{C_{m}}\) is a priori model covariance matrix. \(\varepsilon\) is referred to as the damping factor. Through the above formula, the following formula can be obtained:  
\[S(\mathrm{m+}\delta\mathrm{m})=S(\mathrm{m})+\sum_{j=1}^{n}\mu_{j}\gamma^{ \mathrm{A}}\mathrm{a}^{j}+\frac{1}{2}\sum_{j=1}^{n}\sum_{k=1}^{n}\mu_{j}\mu_{k }\big{[}\mathrm{a}^{k}\big{]}^{T}\overset{\wedge}{\mathrm{H}}\big{[}\mathrm{ a}^{j}\big{]} \tag{11}\]  
and locating the minimum of \(S\) with respect to \(\mu\):  
\[\frac{\partial S(\mathrm{m})}{\partial\mu_{q}}=\gamma^{\mathrm{\overset{ \wedge}{T}}}\mathrm{a}^{q}+\sum_{k=1}^{n}\mu_{k}\big{[}\mathrm{a}^{k}\big{]}^{ T}\overset{\wedge}{\mathrm{H}}\big{[}\mathrm{a}^{j}\big{]}=0 \tag{12}\]  
for \(q=1,\ldots,n\). \(\mu\) can be expressed as:  
\[\mu=-\bigg{[}\mathrm{A}^{T}\overset{\wedge}{\mathrm{H}}\mathrm{A}\bigg{]}^{- 1}\mathrm{A}^{T}\overset{\wedge}{\gamma} \tag{13}\]  
where \(\delta\mathrm{m}=\mathrm{A}\mu\), and \(\overset{\wedge}{\mathrm{H}}=\mathrm{G}^{T}\mathrm{C}_{d}^{-1}\mathrm{G}+ \varepsilon\mathrm{C}_{m}^{-1}\):  
\[\delta\mathrm{m}=-\mathrm{A}\Big{[}\mathrm{A}^{T}\big{(}\mathrm{G}^{T} \mathrm{C}_{d}^{-1}\mathrm{G}+\varepsilon\mathrm{C}_{m}^{-1}\big{)}\mathrm{A} \Big{]}^{-1}\mathrm{A}^{T}\overset{\wedge}{\gamma} \tag{14}\]  
At present, there are two popular methods to evaluate the quality of travel time tomographic solution: detection board test and recovery test. The former adds the positive and negative phase anomalies with the same value to the initial model, and then calculates the path and travel time by forward modeling. The resultant relative travel time residual is obtained by the difference between the results and the travel time of the initial model.