CC BY
13
УДК 512.54
Numerical Evaluation of the Truncated Singular Value Decomposition Within the Seismic Traveltimes Tomography Framework
Alexandr S. Serdyukov*
Novosibirsk State University, Pirogova, 2, Novosibirsk, 630090 Russia
Andrey V. Patutin^ Tatiana V. Shilova*
Institute of Mining, Krasny, 54, Novosibirsk, 630091 Russia
Received 20.08.2013, received in revised form 25.12.2013, accepted 20.02.2014
The method of truncated singular value decomposition is a powerful tool of regularization and solution of inverse problems. Application of this method is limited by the memory requirements for the calculation of singular value decomposition of matrices. The solution is a usage of iterative procedures, that provide the possibility to evaluate only the largest singular values and corresponding vectors. For the practical application of this methodology one should answer the questions: what part of the spectrum can be determined numerically and whether the number of singular vectors is enough for the solution. These problems are considered in the present paper within the seismic traveltimes inversion framework.
Keywords: inverse problems, numerical methods, truncated singular value decomposition, geophysics, seismic traveltimes.
Introduction
The seismic traveltimes tomography is the conventional method of the elastic medium parameters estimation, using the times that are expended for the seismic waves propagation from the sources to receivers [1]. After the linearisation the problem is reduced to the solution of integral equation with the compact operator (tomography operator) in the right hand side. In the case of the least squares solution the most natural regularization approach is the truncated singular value decomposition (SVD). In practice the infinite dimensional tomography operator is replaced by the matrix approximation (tomography matrix), constructed by the discretization of the data and model spaces. Thus one comes to a question: how strong is the identity between the tomography matrix SVD, used for the numerical data inversion, and the SVD of the initial tomography operator. In the present paper we observe the behaviour of the tomography matrix SVD during the disretization step splitting via the model numerical experiment.
* AleksanderSerdyukov@yandex.ru tss3032@yandex .ru ^ ShilovaTanya@yandex.ru © Siberian Federal University. All rights reserved
1. Theory
Firstly we'll give the brief introduction to the seismic traveltimes tomography. Than, we'll introduce the r-solution concept, following the paper [2], which is the theoretical basement of the truncated SVD approach.
1.1. Seismic traveltimes tomography
Seismic tomography is the group of methods for the Earth's subsurface structure recovering, using the measurements of seismic waves. One can use the whole seismograms, recorded in the number of receivers, as the data, that leads to the inverse problem for the elastic wave equation. This approach is known as full wave form inversion or waveform tomography. As an alternative, one can pick up the traveltimes from the seismograms and use them for inversion, that leads (in the high-frequency limit approximation) to the inverse problem for eikonal equation. This kinematic approach, known as traveltimes tomography, is considered in the present paper.
It turns out that the traveltimes are more linearly related to large-scale velocity perturbations than the waveforms, so one should resort to traveltime tomography to resolve these structures. Also the traveltimes inversion is preferable when the dynamical characteristics of the seismic waves (waveforms and amplitudes) can't be recorded accurately, for instance in the case of unmeasurable conditions of wave initializing in the sources. Actually, the kinematic and dynamic inversion methods are complementing each other during the seismic data inversion [3].
The seismic traveltime tomography is conventionally considered in the high-frequency approximation framework, called geometrical seismics. The seismic waves are considered as the propagation of the weak discontinuities of the wave equation solution. Let's consider the isotropic elastic equations:
— + — ( (dUi + duk^ = d2Uk (d
0xk\ Bxi) dxi\ dxk Oxi J P dt2 '
where ui; i = 1... 3 are the components of displacement vector, A, p and p are the Lame parameters and density of the elastic medium. We use the notation of the summation over the repeated indexes in the left hand side of the equation (1).
The weak discontinuities of the solution of the elastic equation (1) can exist only on the characteristics. The characteristic relation leads to the eiqonal equation [4]:
I Vt |2 = V-2 , (2)
V p
where t(x) is the traveltime for the space point x, Vp = (A + 2p)/p. VP is the pressure wave velocity (we are going to consider only this type of elastic waves). The inverse quantity for the velocity, s = 1/VP standing in the right hand side of the eikonal equation is the slowness.
The seismic traveltimes tomography is the inverse problem for the eikonal equation (2). The limitation of the acquisition system and strong curvilinearity of the seismic rays (rays are the characteristics of eiqonal equation) in the complex velocity models make the procedures of back-projection and convolution, that are used in X-ray tomography and are based on the Fourier transform, not really suitable for seismic tomography. The most general and strong for seismic traveltimes inversion is the matrix approach [5], based on the linearization and discretization that reduce the problem to a solution of the large linear algebraic system.
Let's introduce the small disturbance Ss of the known slowness s. Consider the elastic medium with the unknown slowness that is the sum: s + Ss. The traveltimes in this medium would be the sums: t + St. The linearisation of the eikonal equation (2) leads to the following relation for every source-receiver pair:
St = I Ss dl, (3)
r
0
where the integral in the right hand side is taken over the ray trajectory (i.e the characteristic of the eikonal equation) connecting source with the receiver. The problem is to find the small disturbance of the slowness 6s, using the traveltimes residuals dr.
Let's consider the rectangular grid disretization of the target area in the model space. The disturbances of the slowness are supposed to be constant in each cell. The integral from the relation (3) for every source-receiver pair the would be the corresponding sum. Thus the problem is reduced to the determination of the unknown vector of the disturbances values 6s in the rectangular grid cells, using the traveltime data vector 6t , which is composed by the traveltimes residuals for all source-receiver pairs. That leads to the solution of the tomographic linear algebraic system:
BIS =
where B is the matrix approximation of the tomographic operator, called tomography matrix. As it is seen from the relation (3), one needs only to know the length of the ray path through each cell of the rectangular grid to construct the tomography matrix.
In the majority of the researches the effect of the discretization parameters during the matrix approximation on the solution of the tomography problem is explored just by synthetic data inversion. It is clear that such approach is not rigorous and not general enough. Actually, there are some theoretical results, connecting the disretization parameters with the resolution of the tomography method. These results are based on the usage of the Fourier transform. They don't take into account the compactness of the tomography operator and the present of the data and approximation errors [5]. That's why the matrix approximation SVD convergence during the disretization step splitting is need to be observed.
1.2. R-solution
In this section we'll introduce the r-solution concept and some theoretical facts, following the paper [2] Let's consider the operator equation:
Ax = f, (4)
where A: X to Y is linear compact operator, acting from separable Hilbert space X into the Hilbert space Y.
By virtue of compact operator singular value decomposition theorem [6], the operator A : X to Y could be decomposed. It means that there exist orthonormal basis in the model space: {xn G X, n G N} — right singular vectors, orthonormal basis in the data space: {yn G Y, n G N} — left singular vectors and set of nonnegative scalar values: {sn > 0, n G N} — singular values such that:
Axn = snyn,
A * _
A yn snxn?
span{xn} = R(A*) = N (A)±, (5)
span{yn} = R(A) = N (A*)^,
where R(A), (R(A*)) is an image and N (A), (N (A*)) is a kernel of operator A (A*), span(xn) and span(yn) are the linear subspaces, spanned by the n singular vectors. The set {sn, n G N} doesn't have any nonzero limit points.
Note that compact operators have the following important property [7]: the multiplicities of the nonzero singular values are equal to one. The singular values and corresponding vectors are ordered in the singular value descending order: si > s2 > s3 > ... > sn ... .
For any right hand side of the equation (4), satisfying the condition:
f e R(A) © N (A*), the minimal norm extremal solution is defined as follows [8] :
x** = Atf, (6)
where At is the pseudoinverse operator, that is determined by using the SVD (5):
At = ^ oyo x„.
/ sn
sn =0
Definition. r-solution is the element of the model space that is defined by the following formula:
x[r] Atf xn.
sn
n=1
If the minimal norm extremal solution (6) exist, r-solution is represented by its projection on the linear subspace, spanned by the r largest right singular vectors. The operator Atr is called the r-pseudoinverse. This operator is defined and limited on the whole space Y. Moreover, for every f e D(Ar) and VS > 0 always exist r(S), such that for all r > r(S) the following inequality holds: || Arf — Ajf ||< S. Therefore, the set of operators Ar is approximating the pseudoinverse operator Ar on its definition domain. Thus it can be considered as a regularization of the inversion [9].
Within the numerical solution of equation (4) neither the operator A, nor the right-hand-side f are known exactly. One has A and f instead. The approximation and data errors are given as follows:
|| A — A |K Si || f — f ||< S2.
The disturbed right-hand side f may not belong to the definition domain of the pseudoinverse operator even if f e D(Ar). That's why the observation of the difference || Arf — Ar/ || during S12 ^ 0 has no sense. However, since the operators Aj and AJ are defined and limited on the whole space Y, the difference || Aj f — AJ f || specifies the evaluation accuracy of the exact equation (4) solution projection on the first r right singular vectors for the given levels of operator approximation and data errors.
The solution of some inverse problems consists in finding a function, depending on several variables, in some limited target area. As usual, the number of variables is not more than three. In particular, the seismic tomography problems are consist in finding the unknown elastic velocity distribution, that is the function of spatial variables. In practice the data turns out to be a vector, so the data space is finite dimensional. Almost all numerical methods of solving such inverse problems are based on finding decomposition coefficients, corresponding to the finite set of some basis functions. The common choice is the 'step' functions, that equal to one in the grid cells and zeros outside. Thus the passage to the finite model and data dimensions X and Y is made. The equation (4) becomes the large linear algebraic system in this case. The r-solution and other related notions, can be introduced similar way in the case of algebraic system. The singular value decomposition of the matrix corresponds to the SVD (5) of initial operator. It should be noted that, the matrix, approximating the initial compact operator, derives it's properties. The linear system is ill-conditioned. In the case of exact operator the number r for which the r-solution can be stably determined for the known data error level is given by the relation:
|| Sx[r] || ^ si If, (7)
cM
r
where 5x[r] is the residual, coming out from finding the r-solution x[r], using the data f, containing the error f.
The matrix approximation of the initial operator leads to the errors. So, the estimate (7) can't be used. In the paper [2] there is the estimate of the r-solution accuracy in the case of inexact operator. From this estimate follows that the r-solutions of the linear algebraic systems, that approximate initial operator, are converging to the r-solution of the initial operator by the splitting of the disretization step. However, the estimate doesn't provide the number of singular vectors that ensure the stability of solution for the given discretization parameters. The goal of our research is to produce the numerical exploration of the SVD of the seismic traveltime tomography operator matrix approximation.
2. The numerical study of the tomography matrix SVD convergence during the discretization step splitting
As was previously pointed out, the stable recovering subspaces are spanned by the 'largest' right singular vectors. Let's observe the effect of the discretization step splitting on the tomography matrix SVD by considering the model example. The target area is the square, on the opposing sides of which N sources and N receivers are placed, see the Fig. 1. Such acquisition geometry is common for cross-well tomography. For the reason of simplicity and numerical experiment purity, we suppose that the initial undisturbed (containing) medium is homogeneous. In this medium the rays are just straight lines, so the tracing for every source-receiver pair and the ray path length computation for each grid cell are evaluated analytically and don't bring the additional bias. The initial infinite dimensional operator appears when N to the sources and receivers are placed in every point of the corresponding side of the square and velocity is the two dimensional function, defined in the considered target area.
Fig. 1. Model exaple: aqusition system and discretization
We explore the behavior of the singular spectrum via increasing the number N. Let's compare the first largest 500 singular values and right singular values for the sequence: N = 100, 200, 400, 800.
We have performed numerical SVD by using the ARPACK package. Although the tomography matrix is sparse (and stored in the special format) the memory space of usual PC was not enough in the case N = 800, so the numerical computations were made on four-processor cluster with common-pool with 256 GB memory space. The codes were automatically threaded and optimised.
On the Fig. 2 on can see the plots of the relative residuals between the singular values:
SN,2N (i) =
SN (i) - S2N (i) S2N (i) '
where sN (i) is the i-th singular value of tomography matrix for the corresponding N.
,051-1-1-1-1-1-1-1-1-1-
0 50 100 150 200 250 300 350 400 450 500
c)
Fig. 2. The convergence of singular spectrum for the cross-well tomography. The relative difference between singular values for a) N = 100 and N = 200, b) N = 200 and N = 400, c) N = 400 and N = 800
As it seen form the Fig. 2 the spectrum begins to converge from the largest singular values during the mesh splitting. As follows form the change of the singular spectrum by the duplication of the number N, the tomography matrix approximate the initial infinite dimensional operator at least with first order accuracy. Using the given comparison one can estimate the precision of singular values numerical evaluation.
The right singular vectors produce the basis in the model space. Thus the accuracy of their
computation is also important. For the convergence measurement lets consider the following relation: _
A (.) II VN (i) - V2N (i) || (9)
0N,2N(i) = -1, t/ .-.sm--(9)
11 V2N(i) 11
where V2N(i) is the i-th right singular vector for the 2N x 2N grid. VN(i) is the vector defined on the same grid 2N x 2N, produced from the vector VN (i) computed on the grid N x N by linear interpolation.
On the Fig. 3 on can see ¿Nj2N for the sequence N = 100, 200, 400 for the cross-well tomography. From these plots one can estimate what part of the singular vectors are computed reliably for the given grid. Note, that the first 'largest' (corresponding to the largest singular values) right singular vectors are computed more accurately. For instance, from the Fig. 3 it is seen, that about 70 right singular vectors are the same for the whole sequence of N
By observing the singular vector images, one can define some properties. The sequence of singular vectors Vgoo(i), i = 10, 68, 376 is shown on the Fig. 4 a), b), c). Note, that the 'pattern' structure becomes more fine during the rank number increasing. This means, that low spatial frequency velocity structures are restored better during the inversion. One can also see, that the 'pattern' is more fine in the middle part of the target area. Thus, finer velocity structures can be resolved there.
Note that resolution of the spacial grid is not enough for the imaging of singular vectors starting from some large number. As it seen from the Fig. 4 c), the grid 100 x 100 turns out to be too coarse for the approximation of 376-th right singular vector. It is different from the vector, shown Fig. 4 b). From the point of r-solution concept, it is more important to explore the convergence of the right singular values subspaces instead of single vectors. The measure of the two subspaces equivalence is the angle between them (see [10]). Let's introduce the definition of this angle. Consider the Hilbert separable space H. Let Xi, X2 be the two its N-dimensional subspaces. Introduce the matrix:
G =
N
(1) x(2)) '
(x)'', xj(
composed of the scalar products of the orthonormal bases of X1 and X2. The orthogonal substitutions of the X1 and X2 bases lead to the matrix transformation by the orthogonal matrices V, U. Using the SVD theorem for the matrix G, one can find V and U such that there would be the diagonal matrix in the result:
G = S = diag(aj )
with the elements aj > 0. This means that biorthogonal bases are chosen in the subspaces Xi and X2.
(x(1) j
(X((1) ,xj.2)) =0, i = j
(10)
where i, j = 1,..., n. By the definition oj = cos ^, where n/2 ^ ^ ^ 0 are the angles between vectors xj1^ h xj2). The increase-order sorted set of these angles ,^2,.. } is called the series of angles between the subspaces X1 and X2, the maximum of them is the angle between subspaces.
On the Fig. 5 the series of angles between subspaces, spanned by the i largest singular vectors for N = 100 and N = 200, are given. The vectors of the coarse grid are transported to the finer grid by the linear interpolation. From the Fig. 5 it is seen that the angles are so small that the accuracy of their computation is doubtful. From the similar plots for the finer grids one can see
1
a3
c)
Fig. 3. The residuals between singular vectors for different grids; a) N=100 h N=200, 6) N=200 h N=400, b) N=400 h N=800
that the order of angles is decreasing about two times by the dividing into two the mesh step. We will not present here the other plots, since they are spurious. Although the angles between subspaces seem to be small, it is hard to deduce the subspaces equivalence from that, because the dimensions of the considered spaces are huge.
Let's consider the example of the singular value decomposition usage for the solution of cross-well tomography problem. The target region is the square 400x400 meters, split by the rectangular grid with the step 1 meter. The 400 sources and 400 receivers were used. The acquisition system was similar to the previous (see Fig. 1). The synthetic velocity model and the result of the inversion are shown on the Fig. 6. The homogeneous medium with pressure wave velocity equal to 2700 m/s was used as initial model. The 400 hundred 'first' singular
value vectors was used for synthetic data inversion. The usage of the large part of the singular spectrum has no sense, cause the rest part of the spectrum is computed unreliably for the given spatial discretization step.
Fig. 4. Right singular vectors VN(i): a) N = 800, i = 10, b) N = 800, i = 68, c) N = 800, i = 376, d) N = 100, i = 376
Fig. 5. The right singular subspaces angles for N=100, N=200
2640 2660 2680 2700 2720 2740 2760 2780m/c
50 100 150 200 250 300 350 400 m
Fig. 6. The synthetic velocity model (on the top) and its tomography revovery resut
3. Conclusion
The analysis of the marix approximation SVD convergence by the spatial grid splitting for the tomography cross-well operator has been performed in the present paper.
The reliably computed part of the singular spectrum has been picked out for the given discretization and model example.
From the presented research follows that the spatial frequencies of the right singular vector are increasing with their rank. There is a number, such that the vector pattern can't be approximated using the given rectangular grid. The comparison of the singular spectrum for the sequence of splitting meshes helps to find the regularization parameter - the number of largest singular values and vector, used for the inversion.
From our point of view the analysis of the Fourier transform of the right singular vectors is
the perspective direction of the study. The left singular vectors subspaces are also need to be observed, because they form the basis in the data space. We believe, that their study would allow to develop some new algorithm of data filtration and preprocessing.
The research was supported by the Russian Federation President grant for young scientists (number ЫК-2598.20Ц.5).
References
[1] L.Hatton, M.H.Worthington, J.Makin, Seismic data processing: theory and practice, Merlin Profiles Ltd., 1986.
[2] V.A.Tcheverda, V.I.Kostin, R-pseudoinverse for a compact operator, Sib. Elektron. Mat. Izv, 7(2010), 258-282 (in Russian).
[3] P.Mora, Inversion=migration+tomography, Geophysics, 54(1989), 1575-1586.
[4] S.V.Goldin, Introduction to geometric seismic, Novosibirsky Gos. Univ., 2005 (in Russian).
[5] F.Natterer, The Mathematics of Computerized Tomography, Vieweg+Teubner Verlag, 1986.
[6] T.Kato, Perturbation theory for linear operators, Springer-Verlag, 1995.
[7] L.V.Kantorovich, G.P.Akilov, Functional analysis, Oxford etc., Pergamon Press, 1982.
[8] M.Z.Nashed, G.F.Fortuba, A Unified Theory of Generalized Inverses, Generalized Inverses and Applications, Academic Press, New York-San Francisko-London, 1976.
[9] A.N.Tikhonov, V.Y.Arsenin, Methods for solving ill-posed problems, Nauka, Moscow, 1979 (in Russian).
[10] S.K.Godunov, A.G.Antonov, O.P.Kiriluk, V.I.Kostin, Guaranteed Accuracy in Numerical Linear Algebra, Kluwer Academic Publishers, 1993.
Вычислительные аспекты применения метода усеченного сингулярного разложения при решении задачи сейсмической томографии
Александр С. Сердюков Андрей В. Патутин Татьяна В. Шилова
Метод усеченного сингулярного 'разложения является мощным методом 'регуляризации и решения обратных задач. Применение данного метода на практике сильно ограниченно объемами машинной памяти, требуемой на вычисление сингулярного разложения матриц большой размерности. В этой связи актуально использование итерационных процедур поиска старшей части сингулярного спектра. При этом возникают вопросы, какая часть спектра может быть найдена достоверно и достаточно ли найденное количество сингулярных векторов для построения решения. Данные вопросы исследуются в настоящей работе численно на примере обратной кинематической задачи геофизики.
Ключевые слова: обратные задачи, численные методы, усеченное сингулярное разложение, геофизика, времена пробега сейсмических волн.
