Seismic data inversion by cross-hole tomography using geometrically uniform spatial coverage

^IInstituto Geofísico Sismológico Volponi, Facultad de Ciências Exactas Físicas y Naturales, Ignacio de Ia Roza y Meglioli. Rivadavia, CP. 5400, San Juan, República Argentina. Phone: [phone]; Fax: [phone]- E-mail: [e-mail]

^IIDepartamento de Matemática, Facultad de Ingeniería, Av. Libertador Gral. San Martin, 1109 (O) C.P 5400, San Juan, República Argentina- E-mail: [e-mail] ³Center for Applied Geomaterials Research. Atlanta, Georgia. USA. http://geosystems.gatech.edu. Postal address: Instituto Geofísico Sismológico Volponi, Ruta 12, Km 17, C.P. 5413, Rivadavia, San Juan, República Argentina - E-mail: [e-mail]

ABSTRACT

One of the inversion schemes most employed in seismic tomography processing is least squares and derived algorithms, using as input data the vector of first arrivals. A division of the whole space between sources and receivers is performed, constructing a pixel model with its elements of the same size. Spatial coverage is defined, then, as the sum of traveled length by all rays through every pixel that conform the medium considered. It is related, therefore, with the source-receiver's distribution and the form of the domain among them. In cross-hole array, rays do not evenly sample the properties of the medium, leading to non-uniform spatial coverage. It is known that this affects the inversion process. The purpose of this paper, then, was to study the problem of spatial coverage uniformity to obtain travel path matrices leading to inversion algorithms with better convergence. The medium was divided in elements of different size but with an even spatial coverage (named as 'ipixels'), and then it was explored how this improved the inversion process. A theoretical model was implemented with added noise to emulate real data; and then the vector of measured times was generated with known velocity distribution. Afterwards an inversion method using minimum length solution was performed to test the two domain divisions. The results showed that the fact of using ipixels not only improved the inversion scheme used in all cases; but in addition allowed to get convergence where it was impossible to do using pixels; particularly through the method considered. This is a direct result of the improvement of condition number of the associated matrices.

Keywords: seismic tomography, cross-hole, spatial coverage, pixel, ipixel.

RESUMEM

Uno de los esquemas de inversion más empleados en el procesamiento de datos de tomografía sísmica es ei de mínimos cuadrados y algoritmos relacionados, utilizando el vector de primeros arribos como datos de entrada. Se Neva a cabo una division del dominio completo entre emisores y receptores, con el objeto de disenar un modelo de píxeles dei mismo tamano. Se define Ia cobertura espacial como Ia suma de los tiempos de viaje de todos los rayos en cada uno de los píxeles que conforman el medio. Por Io tanto este parâmetro está relacionado con Ia distribución emisor-receptor y con la forma dei dominio entre los mismos. En ei dispositivo cross-hole los rayos no muestrean de igual forma ai medio, conduciendo a una cobertura espacial no uniforme. Se sabe que este inconveniente afecta ei proceso de inversion. El propósito de este artículo fué ei de estudiar ei problema de Ia uniformidad de Ia cobertura espacial a fin de lograr matrices de tiempo de viaje que conduzcan a algoritmos de inversion con mejor convergência. El medio se dividió en elementos de diferente tamano pero con cobertura espacial uniforme (denominados 'ipixels'). Se implemento un modelo teórico con ruido a fin de simular datos reales; y el vector de tiempos se calculo con una distribución conocida de velocidades. Luego se probo Ia convergência de Ias dos formas de division del dominio utilizando ei método de solución por mínima longitud del vector de tiempos. Los resultados demostraron que ei hecho de emplear ipixels no solo mejoró la inversion en todos los casos, sino además permitió lograr convergência en casos donde resultó imposible utilizando pixeles; particularmente con el método utilizado. Este es un resultado directo dei aumento de Ia condición de las matrices asociadas.

Palabras-clave: tomografía sísmica, cross-hole, cobertura espacial, pixeles, ipixeles.

INTRODUCTION

Cross-hole travel-time tomography surveys consist of registering the first arrival pulses of waves that travel through subsurface between seismic sources and receivers (e.g. geophones, hydrophones) located at opposite sides between boreholes (Sheriff & Geldard, 1995) (Fig. 1).

When attempting to model anomalies between boreholes in geophysical exploration, it is necessary to manage some kind of inversion scheme after the pre-processing tasks of the survey data, to be able to assess some kind of geophysical model about the medium involved. The tomographic result consists, generally, in a slowness or velocity distribution through the domain that permits to discover variations in it, which will lead to detect the anomalies sought, if any.

In this work the ray theory will be employed, considering both the host medium and the inclusion as homogeneous (with different velocities), rendering straight-ray propagation (Kolsky, 1963). This approximation is well suited for depths greater than 10m as was demonstrated by Imhof (2007). At smaller ones, stress-dependent anisotropy and heterogeneity are present producing ray bending, therefore eliminating the independence between the model raypath matrix and the parameters investigated, and so complicating the inversion (Santamarina & Reed, 1994; Santamarina et al., 2001).

After the acquisition and through the picking and preprocessing tasks, the vector of travel times _y^<meas> is formed with Mx 1 length; being M number of measurements (i.e. equations). The problem, then, lies to apply an inversion technique to solve the system:

where S is the M x N raypath matrix that represents the model and _x^<est> the N x 1 vector of unknown slownesses (or velocities) distribution.

Considering x^<est>, it is important to quote that generally two forms of representation exist in order to face the inversion scheme:

a) Parametric based: the medium is considered with few unknowns: host velocity V_host (constant), inclusion velocity V_inc (if any, constant too), parameters that locate it: x_inc, y_inc\ and finally other that permit to assess its form and size, i.e. circular or elliptical (see Santamarina & Cesare, 1994; Imhof & Calvo, 2003; Imhof, 2007). This will lead to an over-determined system of equations that, when possibly, is solved with variational or least squares inversion algorithms, for example.

b) Pixel based: The medium is divided in N uniform (i.e. same size) pixels; each of them is an unknown with a particular estimated (after the inversion) velocity or slowness value. When the inversion process is finished, an image reconstruction will show the position of each element and its velocity value (or colour pattern related). This will show the position of inclusion through colour contrast (Tarantula, 1987; Menke, 1989; Fernandez, 2000). It is clear that the smaller the pixels, the better the image resolution obtained, but at the cost of N increment, leading so to an under-determined system of equations.

Though in (b) a uniform pixel representation was mentioned, other forms to divide the medium can be studied; since the usual form leads to uneven distribution in spatial coverage (see Theory).

The primary objective of this paper is to implement a different form of discretization of the medium with the purpose to study first the conditioning and rank of several raypath matrices made with those distinct elements assembly; and second use them to study their behavior with a typical inversion method using a theoretical model constructed to that end.

THEORY

Forward and inverse problem

Keeping in mind straight ray theory, the forward problem permits to compute the predicted (output) travel time vector _y<P^red> (M x 1) as a function of the vector of known true slownesses (in-put)x^<true> (A^x 1),andthetransformationmatrix S(MxjV), which represents the travel path that connects 'x to 'y';

In the inverse problem, the M values collected in vector y are measured (y^<meas>), and the aim is to estimate the N unknown parameters x^<est>:

In general (as in this case) the inverse S^_1 cannot be determined, and a pseudo inverse must be used instead (Penrose, 1955):

In a well posed problem the solution exists and is unique and stable (Santamarina & Fratta, 1998). The problem here is that cross-hole tomography is not a well posed problem because of the relative distribution of sources and receivers (Branham, 1990; Tarantola, 2005) that involves rays almost parallel which increase condition number of associated s matrix.

The procedure is to calculate first the appropriate matrix s of travel paths, then evaluate its pseudoinverse s^<P^seudoinv> (depending of the method of inversion considered) and last to apply Eq. (4) to find x^<est>.

Considering straight path simplifies very much the formulation at hand because the matrix s must be assembled only once (explicit form).

Raypath matrix

The entries in S (Fig. 2) are calculated identifying the intersection of the individual 'm ray with the 'n' pixel boundaries and computing the Pythagorean length 'd' inside the element.

Therefore, the size of S will grow with the increment of measurements and/or discretization density.

Considering again Figure 1 and having still in mind Figure 2; the number of sources and receivers there is 7 and so the total number of rays is 49 (constant for this array). If the medium is divided in 7x7=49 pixels, S will be of 49 rows (M, number of rays) per 49 columns (N, number of pixels, whose slownesses are unknown). It is relevant to note that the increment of the pixel density (N), will lead to consider more pixels in the computation of travel paths for each receiver.

This system of equations is apparently even-determined (equal number of equations than unknowns). The word apparently means that in some cases (especially in cross-hole, see Imhof, 2007) the system is ill conditioned and the rank of S < N, M. This will give an underdetermined system leading to infinite solutions. The same will occur, for example, dividing the medium in more pixels to improve the resolution of the images, (the limitation of dividing the medium with N = M is that the resolution is coarse, because the size of the pixels is large).

Incrementing the number of transducers to improve resolution is not practical and always possible to do; first, due to the amount of survey effort needed (cost) and, second, if rays are so near, the condition of S matrix augment and not necessarily add information to the system (Santamarina & Fratta, 1998; Fernandez, 2000).

Matrix of Spatial coverage (Sc)

As quoted, the size of S is M rows (number of measurements) per N columns (number of unknown pixel slownesses). The sum of each individual column of S brings the total length traveled by the rays in one pixel (Fig. 2). Applying this summation in all columns of S, gives a 1 x N row vector that re-arranged following the geometric pattern, brings the Sc matrix of V vertical pixels by Y horizontal ones where N = s +1. It was considered here s = t.

Two examples of cross-hole Sc matrices are represented in Figure 3 for 10 source-receiver' pairs; (a) for 100 elements and (b) 400 ones. The dark zones indicate lower values of spatial coverage. This means that the Information gathered to solve the velocity or slowness of It Is lower than In other lighter zones In other words the precision for the evaluation of the pixel values will not be uniform. Due to more rays traversing one pixel get more information of it (similar to CDP concept in reflection seismics), the lightest sector at the center of Sc matrix depicts the maximum resolution and precision for the pixel/s situated there. But what happens when an inclusion being searched for is far from that position? The resolution to locate it will be poorer (Santamarina & Fratta, 1998).

An alternative form to divide the medium to improve the resolution in the dark zones is proposed: Instead of separate it in pixels of same size and different spatial coverage at each; it'll be divided in elements of equal spatial coverage and distinct individual sizes and named as ipixels. Due to this fact, the elements of Swill be different. Figure 4 shows the domain divided in two densities of ipixels. The same color tones depict same information at each element.

Any type of element discretization for the physical domain is perfectly possibly to accomplish, because any type of it is only geometrical and has for purpose to make the s matrix and arm the system of equations having the distances traveled by the rays.

Conditioning of a matrix

If there is a high contrast between numerical values of the elements inside a matrix; it is possible that the calculated rank for it has an erroneous linearly independent number of rows. Due to this, the condition number κ is a better choice to study the conditioning of a matrix (i.e. how near it is of a matrix of lower rank) than its rank (Strang, 1980).

κ is defined as the ratio between the singular values with maximum and minimum absolute values (Golub & Van Loan, 1989):

A matrix is said to be ill-conditioned if κ is very large.

Geometrically, the maximum and minimum values represent the axes of an ellipse (Branham, 1990). If the ratio tends to unity the correlation is null (independent) and if it is very large, that is almost perfect (dependent) and so ill conditioned. In other words, ill conditioning means that the number of linearly dependent rows or columns of the matrix tends to increment. This means that the number of equations/measurements decreases related to the number of unknowns.

Lastly, the singular values give an indication related to confidence and the relationship among the measurements of the phenomena. Small values suggest limited information about the parameter.

Inversion: Minimum Length Solution (MLS)

Pseudoinverse matrix

Data and model resolution matrices

When the system is over-determined the least squares inversion method find a set of values _x^<est> that minimizes L₂-norm of the vector of residuals:

But, with under-determined systems, this is not possibly to do, since the problem has infinite solutions.

One form to face this difficulty is to find the minimum length of the vector _x^<est> that satisfies the u constraints _y^<meas> = S • x (Santamarina & Fratta, 1998; Imhof, 2007).

To solve Eq. (4) it is necessary to calculate the pseudoinverse s^<pseudoinv>. This is a matrix that fulfills four conditions (see Penrose, 1955) and whose products S • s^<pseudoinv> and/or s^<pseudoinv> . s are not necessarily equal to the identity matrix I.

For underdetermined and consistent MLS method the following expression gives the pseudoinverse (Santamarina & Fratta, 1998):

The data (D) resolution matrix renders the difference between y<pred> _ancj y<meas> :

where:

Finally, the model (M) resolution matrix depicts the difference between _x^<est> and _x^<tme>:

where:

The ideal case develops when D = I and/or M = I which means that y^<pred> = _y^<meas> (nu|| data error) and/or x^<est>_=x^<true> (nu|| model error).

METHODOLOGY

Two Matlab programs were developed to construct Sc and Sp matrices for any size of regularly distributed rectangular pixels and represent them in graphical form (examples in Fig. 3).

Then, with other two Matlab programs, several Sc of irregular pixels derived from the previous high density uniform pixel ones were constructed grouping them with the characteristic of bringing the same spatial coverage at each one (Fig. 4). Finally, new raypath matrices Si were conformed to those Sc of ipixels.

Afterwards the conditioning of Sp and Si was tested. Table 1 bring the values of κ and rank for some of them and Figure 5 shows a graphical representation of singular values.

Figure 5

Theoretical model and Minimum Length Solution (MLS) method

The designed model consisted in a high velocity inclusion (V_inc = 2000 m/s) located at known coordinates and size, insert in a medium of background velocity V_back = 400 m/s (Fig. 6). At the left border of the medium ten sources were considered. At the opposite side were positioned equal number of receivers, rendering a total number of 100 rays.

Afterwards, the medium was discretized; S was constructed with a ray tracing algorithm and the vector x^<true> formed, following the geometry. Then y^<Pred> was calculated applying Eq. (2): y^<Pred> = s x^<true>. Finally random noise was added to form y^<meas>.

To solve the under-determined and consistent problem in matrix form the minimum length solution was applied, using the following equations.

Replacing Eq. (6) in Eq. (4) renders:

For this solution the corresponding generalized inverse s^<pseudoinv>. Data (D) and Model (M) resolution matrices were defined in Eqs. (6); (8) and (10), considering S as Sp or Si depending on the case.

Afterwards, x^<est> (slownesses) was inverted, calculating the velocity vector _v^<est> and represented graphically when possible, assigning the velocities to the pixel/ipixel values at the corresponding x-y positions.

RESULTS AND ANALYSES

Conditioning and rank of the raypath matrices Sp and Si

Dividing the medium in 10x10=100 pixels (or ipixels) imply that the raypath matrix Sp (or Si) will be of M rays (rows representing the measurements) per 100 pixels (columns). If the medium discretization is denser, the number of columns will be greater. For example a 20x20=400 pixels will give a matrix Sp of M x 400 size. As mentioned, in this work M = loo (10 sources and 10 receivers).

Table 1 depicts the condition number κ and rank of several matrices Sp and Si.

The κ parameter is lower with T than with 'p'. Besides the rank is greater in the first case. All these suggest that the inversion scheme would function better using Si than Sp, as predicted earlier using SVD decomposition.

Inversion results. Model parameters

Table 3 summarize the inversion results for several matrix configurations. Six parameters were determined, including height and width of the rectangular inclusion. Figure 7 shows the inversion solution of four cases with ipixel used.

In all cases, it was impossible to get model parameters with MLS method using Sp matrices (without additional information, which is not the case studied here). This was not the fact using Si ipixel ones.

The cases 2-6 show that the resolution in Vhost increases when more pixels are considered to form Si. V_inc is more difficult to define precisely, especially considering that here was not added any a-priori information to the inversion scheme. Nevertheless the variations rounded approximately 8% from the Vinc^<true>.

The pair Xinc, Yinc presented little variations and was well determined in all cases, but important is to note that increasing cases in Table 3 turned more difficult to define precisely the visual location of the anomaly (analyze four cases in Fig. 7). These applied to calculate the height and width of the inclusion.

Besides, the increment in number of unknows (more ipixels) brought non-definitions in the position and size of the inclusion. The minimum length solution method consider the minimum length of the solution vector x^<est> and caution must be taken because it has an infinite number of solutions and the adequate must be chosen.

The best results were obtained employing up to N ~ 4 . M unknowns.

Inversion results. Data and model errors

Table 2 shows L₂ norms of Data and Model errors for several cases of matrices Sp and Si. Minor the error, better the fitting. It is clearly seen that the data error is minimum in all cases, due to MLS method resolved the data perfectly. Therefore, the model error column defined the fitness quality of the inversion.

Considering, for example 20x20 pixel and ipixel sc matrices to form respectively Sp and Si and resolve the inversion, it was appreciated (Table 2, Cases 1-6):

a) The fact to divide the medium in pixels of the same size increase the model error of the inversion, as it is seen in case 1 respect to 2-6.

b) Considering now only cases 2-6; the model error decreases when more number of former pixel were considered to construct the subsequent ipixel Si matrix. For example consider the cases when those were 100x100 pixel with Me=2.5979; and 357x357, that brought a Me=1.0417. Probably this is related to an increase in precision to evaluate the ipixels.

c) Denser the original Sc matrix and coarser the ipixel Sc derived from if will signify that minor will be the spatial coverage difference between the ipixels. Conversely, coarser the former pixel matrix, more difficult will be to obtain precision in the evaluation of the irregular pixels.

d) The analysis made above is demonstrated again, in cases 7 and below.

CONCLUSIONS

The MLS method functioned properly when the spatial coverage was uniform (as it was in all ipixel cases), bringing to acceptable results inclusive when the density of ipixels was very high. In all the other cases with pixels conforming the different raypath matrices, it was impossible to obtain a convergence to a solution.

The best inversion results were obtained with a maximum number of 20x20=400 ipixels, that is four times the data number (i.e. number of measurements). Although it is possible to obtain a solution with greater ipixel densities, nevertheless diminishes the spatial resolution of the image. In these cases and if the ipixel size is important to image a smaller inclusion, will be necessary to add additional information (such as a vector of background velocities Vhost) to get a better image.

Finally, it was demonstrated that (at least in MLS) using ip improves convergence in the inversion process.

ACKNOWLEDGEMENTS

This work has been partially supported by Consejo de Investigates Científicas y Tecnológicas (CICITCÂ), Universidad Nacional de San Juan, Argentina.

Recebido em 25 fevereiro, 2008 /Aceito em 18 abril, 2009

Received on February 25, 2008 /Accepted on April 18, 2009

NOTES ABOUT THE AUTHORS

Armando Luis Imhof. Geophysicist graduated at Universidad Nacional de San Juan (UNSJ, Argentina) in 1989. Specialist on Environmental Studies from Universidad Nacional de San Luis (Argentina), 2002. PhD degree from Engineering Faculty at Universidad Nacional de Cuyo, Mendoza, Argentina, 2007. First Prize winner of XIV Convocation, Garcia Sinériz Foundation (Spain) to the Best PhD Thesis in Geophysics from Spain and Latin America. Researcher from the staff of Instituto Geoffsico Sismológico Volponi (Facultad de Cs. Exactas, Fcas y Naturales - UNSJ) and Adjoint Professor of the asignature Electrical Prospecting, Geophysics and Astronomy Department (UNSJ) since 1990. Member of Asociación Argentina de Geofísicos y Geodestas (AAGG). Topics of main Research: Geophysics applied to mining, geotechnics and hidrogeological studies.

Carlos Adolfo Calvo. Electromechanical Engineer graduated at Universidad Nacional de San Juan (UNSJ, Argentina) since 1978. MSc in Applied Mathematics and Informatics from Moscow Physics Institute (MEPHI). Main Professor of the asignatures Mathematics Analysis, Applied Mathematics and Numerical Methods; Engineering Faculty (UNSJ). Researcher on topics covering applied mathematics to solve engineering problems. Topics of main research: Applied mathematics.

Juan Carlos Santamarina. Civil Engineer from Universidad Nacional de Cordoba, 1982; MSc in Geotechnical Engineering at University of Maryland, 1984. PhD in Geotechnical Engineering at Purdue University, 1987. Co-director at the Center for Applied Geomaterials Research (Georgia Institute of Technology); Goizueta Foundation Professor at Georgia Institute of Technology, Atlanta, USA.

• BRANHAM RL. 1990. Scientific Data Analysis. An Introduction to Over-determined Systems. Springer Verlag. NY. 237 pp.
• FERNANDEZ AL. 2000. Tomographic Imaging the State of Stress. PhD Dissertation. Georgia Institute of Technology. Atlanta. USA. 298 pp.
• GOLUB GH & VAN LOAN CF 1989. Matrix Computations. Baltimore: John Hopkins University Press, 642 pp.
• IMHOF AL. 2007. Caracterización de Arenas y Gravas con Ondas Elásticas: Tomografía Sísmica en Cross-Hole. PhD Thesis. Universidad Nacional de Cuyo. Mendoza. Argentina. First Prize Award Fundación Garcia Siñériz for the Best Doctoral Thesis in Pure or Applied Geophysics. XIV Convocation. Universidad Politécnica de Madrid. 299 pp.
• IMHOF A & CALVO C. 2003. A Variational Formulation to Image Inclusions in 2D Travel Time Cross-Hole Tomography. Brazilian Journal of Geophysics, 21(3): 269-274.
• KOLSKY H. 1963. Stress Waves in Solids. Dover Publications Inc. NY. 213 pp.
• MENKE W. 1989. Geophysical Data Analysis: Discrete Inverse Theory. Academic Press, Inc. NY. USA. 289 pp.
• PENROSE RA. 1955. A Generalized Inverse for Matrices. Proceedings of Cambridge Phil. Society, 51:406-413.
• SANTAMARINA JC & CESARE MA. 1994. Velocity Inversion in the Near Surface: Vertical Heterogeneity and Anisotropy Internal Report. University of Waterloo. Canada.
• SANTAMARINA JC & REED AC. 1994. Ray Tomography: Errors and Error Functions. Journal of Applied Geophysics, 32: 347-355.
• SANTAMARINA JC & FRATTA D. 1998. Introduction to Discrete Signals and Inverse Problems in Civil Engineering. ASCE Press. USA. 327 pp.
• SANTAMARINA JC, KLEIN KA & FAM MA. 2001. Soils and Waves. Wiley & Sons. Ltd. England. 488 pp.
• SHERIFF RE & GELDARD LP. 1995. Exploration Seismology. 2^nd Edition. Cambridge University Press, NY. 592 pp.
• STRANG G. 1980. Linear Algebra and Its Applications. Academic Press, NY. 414 pp.
• TARANTOLA A. 1987. Inverse Problem Theory. Elsevier. Amsterdam. 613 pp.
• TARANTOLA A. 2005. Inverse Problem Theory and Model Parameter Estimation. Society for Industrial and Applied Mathematics. SIAM. Philadelphia. USA. 333 pp.

Seismic data inversion by cross-hole tomography using geometrically uniform spatial coverage

Armando Luis Imhof^I; Carlos Adolfo Calvo^II; Juan Carlos Santamarina^III

Publication Dates

• Publication in this collection
  01 Dec 2010
• Date of issue
  Mar 2010

History

• Received
  25 Feb 2008
• Accepted
  18 Apr 2009