1. No category

Reflection tomography : how to take the spatial
correlations between traveltime data into account ?
F. Renard and P. Lailly
ABSTRACT
In the standard formulation of reflection tomography, the norm chosen in the data space for
measuring the mismatch between the observed and calculated traveltimes is the Euclidean one :
the picked traveltimes are supposed to be uncorrelated pieces of information. In addition, their
spatial sampling (discretization) is irregular and somewhat random. Thus, the solution of the inverse problem, and all the associated quantities (a posteriori confidence in particular), may be, at
least in principle, dominated by numerical artifacts.
We propose a new formulation of reflection tomography so as to obtain an intrinsic solution,
i.e. a solution that is independent of the discretization. For this purpose, we take into account the
spatial correlations between traveltimes : the data are no more only made up of the traveltimes,
but also of the derivatives of these traveltimes w.r.t. the Common Mid-Point coordinates and the
offset. We develop the calculations required to solve this new inverse problem. Then, we invert
synthetic data so as to check the software we have developped from those calculations.
INTRODUCTION
Correct imaging of the subsurface requires accurate knowledge of the propagation velocity of
seismic waves. For this purpose, the KIM research consortium has developed the reflection tomography software jerry (Jurado et al., 1996, Jurado et al., 1997) that allows to calculate such models
by inversion of reflection traveltimes picked on seismic sections. The standard formulation of this
inverse problem consists in minimizing a least-squares function that measures the distance between the data (the observed traveltimes) and the seismic response of the model (the calculated
traveltimes), in order to recover the subsurface model that “best” fits the data.
For sake of simplicity, an Euclidean norm in the data space is usually chosen : doing so, we
KIM Research Consortium - Institut Franc¸ais du P´etrole, 1 & 4 avenue de Bois-Pr´eau, 92852 Rueil-Malmaison, France.
1
consider the picked data to be uncorrelated pieces of information, which is obviously not true whenever the associated source-receiver pairs are close to one another. To overcome this difficulty, we
formulate a continuum inverse problem allowing to take the spatial correlations of information
into account. Our final goal is to develop a discretization strategy that ensures the inverse problem
with sampled kinematic data to be a genuine approximation of the continuum inverse problem.
After a short reminder concerning the standard formulation of reflection tomography, we detail
all the calculations required to solve the new inverse problem, and, in particular, we focus on
the tedious calculation of ray perturbations, necessary to solve the associated linearized inverse
problem. Finally, we give a preliminary inversion result of synthetic data.
THE STANDARD FORMULATION OF REFLECTION TOMOGRAPHY
Model description
We use a blocky model representation of the subsurface : we parametrize the interface geometries and the velocity within each layer by bi-cubic B-spline functions.
A z-explicit interface is thus defined as
"! ! #
are the $&')
where % ( $*+ % B-spline parameters in and directions, respectively, and ! ! are the B-spline base functions. The B-spline interfaces are defined for the complete model domain
(Clarke, 1996) and thus can cross each other.
Moreover, a velocity law is associated with each layer. In this paper, we deal with general ,.velocities, defined as
where
;
7
5 4 54 8 4 7
7
6 9 ! 6 ! ! 9 =3
/0
123
7
6 9:<;
6 9 are the $*' > ( $?+ > ( $?@ > B-spline parameters in the and directions.
Forward problem
The forward problem of reflection tomography is formulated as a two-point ray tracing : the
source and the receiver are fixed. This problem is solved by a bending method. First we specify the
signature of the ray (defined by the set (source, receiver, list of interfaces)), and we initialize the
trajectory, which is a path linking the source to the receiver, verifying the signature. The traveltime
along this trajectory is given by A2BDCEFBHG JI KD A K , where L is the source and L the receiver,
and A K is the traveltime between impact points L K and L K , these impact points being connected
by straight lines (Figure 1). Then we solve the two-dimensional system
M NA BDCEFBHG
MLK OFPRQSP $UT O"
moving the impact points on the interfaces. The solution trajectory satisfies Fermat’s principle, i.e.
makes the time function stationary. This trajectory is a ray.
2
P0=Source
PN=Receiver
∇
✶
Pn-1
Pn+1
Velocity (V abc )
tn
+ Mn+1
+ Mn
tn+1
Reflector (Z uv )
Pn(xn,yn,zn)
Fig. 1 Definition of a trajectory : a trajectory connects the source to the receiver passing by all impact points on
interfaces. Successive impact points on interfaces are linked by a straight line.
Inverse problem
9 6 the forward modeling operator, that is the operator that gives, for a specified
We denote by
reflector, the traveltimes associated with all source-receiver pairs. The inverse problem consists in
minimizing the objective function defined 7 by (Delprat-Jannaud and Lailly, 1993)
with
and
where 3 O 9 6 T
- 4
!#"
- /$2
&% (' ' +
- T
*)
!,"
+
- #&% (' ' M M M .- #&% 0/ M &1 / M &1 / M M
- #&1
is a subsurface model defined by $ 2 interfaces and $ > layer velocities.
is the norm in the data space and 54 is the a priori covariance operator on the data. Its
"
diagonal elements specify the uncertainties on the data, and its off-diagonal elements specify the
correlations between these uncertainties.
tomography is a non linear inverse problem since the forward modeling operator
9 6Reflection
is non linear. This non linear problem is solved with a Gauss-Newton algorithm, which relies
on successive linearizations of the forward map. The linearized objective function 76 that has to
be minimized at each iteration is
7
.9: O
9= T 9 - 9: T 7
7
>A>=@CG BDFE is the Jacobian
where is a current model, 9= is the model perturbation, ; ?
9
H
6
.
matrix evaluated at , and 9
T
Minimizing this objective function amounts to solving the linear system
7
;*I 4 ; J FK&L 9:U ;*I 4 9 (1)
86
<;
3
where J FK L is the matrix made up of the regularization terms. This linearized inverse problem is
well-posed (Delprat-Jannaud and Lailly, 1993) : if the number of data is sufficient, this problem
has a unique solution that depends continuously on the data (stability) and, from a practical point
of view, the more kinematic data we have at our disposal and the greater the regularization weights
are, the more stable is the solution.
In general, we solve the linearized inverse problem by a conjugate gradient algorithm (see for
instance Chauvier et al., 2000) allowing the velocities and the interfaces to be simultaneously inverted.
For sake of simplicity, the a priori covariance operator is chosen diagonal : the picked traveltimes are supposed to be uncorrelated pieces of information. Thus 4 is also diagonal and the
non linear objective function becomes
7
7
O
A 9 6H T&A - T
7
where $
is the number of data and is the uncertainty associated with data A .
Unfortunately, the assumption of mutually independent kinematic data is not realistic, whenever
the associated source-receiver pairs are close to one another.
Another, closely related problem is : how should the kinematic information be spatially sampled if we want the solution of the inverse problem, and the associated quantities (a posteriori
covariance operator for instance) not to change substantially if we change the spatial sampling of
the kinematic data ?
Figure 2 shows the kind of sampling we use in practice. Near offset data are quite finely sampled,
both along the and directions, whereas multi-offset data are picked along sparse profiles (pa
rallel to the direction), and are thus quite finely sampled along the and directions ( is the
offset) but very sparsely in the direction.
4
x
x
y
y
Fig. 2 A typical survey for a reflection tomography experiment : the near offset data are finely sampled along
and (the CMP coordinates) directions. The multi-offset data are picked on sparse multi-offset lines. They are
thus finely sampled sampled along and (the offset) directions. The spatial sampling of the kinematic data is thus
very anisotropic and as a consequence, the solution model and all its associated quantities are corrupted by numerical
artifacts.
This sampling of the data space is quite practical as it simplifies a lot data management and
interpretation on 3D workstations, but it is also very anisotropic, which may pose some problems.
First, there are directions along which the kinematic data are finely sampled, but although we
know that dense data involve strong correlations, those correlations of information are not taken
into account. Then, there are other directions along which the kinematic data are very sparsely
sampled : in those directions, we may wonder if we have enough information to recover all the
spatially varying features of the subsurface. Consequently, we can suspect a footprint of the discretization of the data space on the solution model. This footprint may be observed in the misfits,
i.e. in the differences between observed and calculated traveltimes. Although the RMS of these
misfits may be quite satisfactory, very often the misfits are not erratically distributed, contrary to
what the choice of an Euclidean norm seems to mean. In fact, these misfits may show organized
fluctuations of high amplitude (several times higher than the RMS value), which is not very satisfactory 1 . All our attempts to make these fluctuations decrease to reasonable amplitudes by making
use of weights in the Euclidean norm have been fruitless : this shows that this norm is not well
adapted for solving the tomographic inverse problem.
Our final goal is to obtain an intrinsic solution to the reflection tomography problem. The
strategy we propose is based on an intrinsic representation of the kinematic information : we do
not want the solution of the inverse problem to depend on the spatial sampling of the kinematic
data. For this purpose, we need first a formulation of the continuum problem taking correlations
of information into account, and then we need to develop a discretization strategy that ensures
the inverse problem with sampled kinematic data to be a genuine approximation of the continuum
problem.
1
One explanation could be the existence of systematic mis-pickings in the data. However, it seems quite unlikely that the interpreter
“misses” the seismic event with such an error !
5
TAKING CORRELATIONS BETWEEN KINEMATIC DATA INTO ACCOUNT
Formulation of a continuum inverse problem
For sake of simplicity, we consider for all the following calculations that the multi-offset lines
have a constant azimuth. As our goal is to take into account the correlated nature of the kinematic
data, a first idea consists in exploiting the Bayesian approach (see Tarantola and Valette, 1982;
Duijndam, 1988).
Then, in the case of a continuum of data ' , where 0 ( and are the CMP
7 generalized least-squares
7 function is given by
coordinates, is the offset), the
' 9 6H T
'
I 4 ' 9 6 T
'
T
I T
where (4 and are the covariance operators a priori on the data and on the model, respectively.
Nevertheless, the kinematic data considered as spatially correlated and the smoothness imposed to
the solution model are two non independent pieces of information (data could not be supposed to
be correlated if the subsurface was not supposed to be smooth). Thus, one of the basic hypothesis
of the Bayesian approach, i.e. the independency of the two states of information a priori on the
data and on the model, is not verified.
To avoid this difficulty we base the formulation of the continuum formulation problem on a determinist approach. Such an approach consists in choosing the norms in the data and the model
spaces, so as to provide “relevant” quantifications of the misfits and the deviations from the a
priori model : these norms must reflect the confidence we have in the seismic data and in the a
priori model, respectively. The norm in data space also must reflect the correlated nature of the
kinematic data.
Thus, we define a new objective function
O
0 9 6 1 T M 9 6 1
'
M 1 M 9 6
O
0
+
M 0 O 1 M 9 H6 1 M
T - O
7
0 7 ' ' ' T 7 1 ' ' T
T
'
7 1 ' ' ' 0 ' ' '
(2)
Here, we have made explicit the dependence of the data on the three spatial coordinates
and . 1 and ' + 0 are weight functions defined
0 O
0 0 Q
'
' + 1 ' + 0 reflects the uncertainty of the traveltime at point 1
typical correlation length at point 0 in the directions ' + 1
1 with
6
and ' + 0 reflect the
, respectively : the greater those
lengths, the stronger7 are the correlations.
With this formulation, the data for tomographic inversion are no longer constituted of the obser
ved traveltimes
only. The slopes of the observed traveltimes in the , and directions ;
>A@
>A@
and >A@ >
, now show up as data. The four functions , ' , + ,
allow us to quan> ' , > +
tify the confidence we can put in the a priori information concerning the traveltime data, whereas
the quantification of the confidence we can put in the a priori model still consists in an adequate
choice of the regularization weights on the interfaces and the velocities.
The inverse problem with sampled kinematic data
Of course, as the traveltime data are picked on seismic sections, we have to solve the inverse
problem with sampled kinematic data. For sake of simplicity, we consider constant sampling in
tervals , and
in directions , and respectively. Our final goal is to approximate the
above continuum objective function (2) by using the discrete objective function
O O 9 6H A T&A ' M M
7
A 9 6H TSA 7
+ M M 9 H6 T&A A
7
M M 9 H6 A
T*A 7
- T
(3)
7
' + , and ' , + and are the correlation lengths for the traveltime A where in , and directions, respectively.
in order to make (3) a genuine apThe main problem that remains is the choice of 0
proximation of the continuum objective function (2). In practice, the chosen sampling intervals
will result from a compromise : they have to be fine enough to ensure the relevance of the approximation of the continuum function, but they have to be coarse enough so as to avoid the
management of huge volumes of data.
The associated linearized inverse problem
The new norm in the data space associated with (3) induces many changes in the solution of the
inverse problem, but the main one is the reformulation of the linearized inverse problem. Indeed,
as we use three new types of data compared with the classical formulation, the Jacobian and the
Hessian of the objective function are modified. The new linearized inverse problem now consists
in solving the linear system (cf. equation (1))
7
7
7
7
+ + ; + ; I ; J FK L 9= 9 ' ; + I + 9 + ; I 9 7
(4)
; I 4
; 'I '
7
where 9 ' + > > '@ + and ; ' + >=> G > >A' @ +B DF E , and where the covariance matrices " 4 ' + 7
7
+ 7
are diagonal and defined by 4 , ' ,
and , O $ % .
9 '
9 +
,
and 9
are simply calculated from the data. The calculation of ; is classical
(see Jurado et al., 1996; Jurado et al., 1998). The calculation of the three new Jacobians ; ' ; + and
;
I 4
;
7
' ' ; '
; I 9 ; I ;
is explicited in the next section.
Computation of ray perturbation
We now detail the calculation of ; ' >> G >A@ > 'BDE ( ; + and ; are calculated similarly).
For this purpose, we consider a subsurface model made up of $ O" z-explicit interfaces defined
2 , and $ 3D velocities / 023 / 123 .
by #
As we work in CMP coordinates, we can use the important result of Adler, 1996,
M 9 6H
M T
where is the projection of (the ray parameter at the source or at the receiver) onto the
direction. As a consequence, the computation of ; ' amounts to calculating
M
M
T
" the source and at the receiver are function of the ray impact points
As the ray parameters at
2
O $ T O % , we have to calculate ray impact points perturbations due to
L
+ @ '
model perturbations : >>G , >>=G , >> G .
Once these Jacobians known, ; ' is computed using
M M+
M
M+ M
$ T O if M
where O
if and M M
+
M
+ M
M M
M + M
+
.
We now address the calculation of ray perturbations. A ray described by the trajectory the source and L the receiver satisfies Fermat’s principle, and thus veri
fies, if we denote by the subsurface model,
L L , with L -
M
F M A F K
M
M K
M A K F M K K < K K OFPRQSP $UT O"
(5)
We refer the reader to Jurado et al., 1998 for the detailed computation of the derivatives required
to make satisfy (5).
Let us now consider a model perturbation 9: : we perturb a B-spline parameter corresponding
either to a velocity or to an interface. The resulting perturbed trajectory 9 still satisfies
Fermat’s principle in the perturbed model 9: -
9
F
8
9:
and, as a consequence, using the first order Taylor expression, the calculation of the impact point
perturbations 9 amounts to solving for given 9=
M
M -
9
7
7
7
M
M
9:
(6)
Case of a velocity perturbation
"
in the case of a perturbation of a velocity B-spline parameter :
6 We
9 now 6 explicit
9 9 the
6 9 calculations
O. $ T O % . In that case, equation (6) yields
7
;
; ;
7
M
M
M A
"
O.
A
O
% 9"/ 6 9
"
9
M M M M A 9"
$UT
M M/6 9
7
7
M
M
M
M
M
9
9
A
M M / A 6 9 M 9"/ 6 9
M
M
M
M
@ M
M
M
M
A
A
9 9 M
M
M
M
M
M
M
M
M M A M 9 M M
AM M M 9 M M A M M M 9 Moreover, A I A 7 , where A is the traveltime between points L and L , and thus only
> > by - A / a function equal to one if
"
> I and > I are different from zero. Finally, if we denote
6
9
( i.e. if the ray segment L L % is within layer ), and equal to zero
A possibly depends on
otherwise, equation (6) ; becomes
9
7
7
7
M A 7
M
M M / 6 9 - A / 9"/ 6 9 M MA / 6 9 - A / 9"/ 6 9
M A 9" M A 9" M A 9 M A "9 M M
M M
M M
M M
M A 9 M A 9 M M
7 M M 7
7 7 M
M A M
M
M / 6 9 M - A / 9"/ 6 9 M M A / 6 9 M - A / 9"/ 6 9
M
M
M
M
M
A A 9 9 M
M
M
M
M
M
M A M 9"
M A M 9"
M M M
M M M
M A M 9 M A M 9 M M M
M M M M A A M M M M M 9" M M 9 (7)
" required can be found in appendix A. This expression is valid for All the derivatives
O. $ T O % , and thus we have at our disposal $ T O" equations for $ T
and , O" unknowns. Consequently, calculating ray perturbations amounts to solving a linear system of
dimension $ T O" " for each ray 'and
for each
B-spline parameter. The @ solution
of such a linear
> +
> .
O
O
% > $ T , >G and >G . Then we easily compute >G , and the Jacobian ; ' is
system gives us computed using
M
M M
M M
M M
where Finally, ;
'
O
if and ;
M + M +
$ T O if M
+
M
+ M
M
+
M
+
and .
are computed in a similar way. Indeed, still following Adler, 1996
M 9 6
M G
T
and if we consider the multi-offset lines to have been acquired in the direction
M 9 6H
M G Case of an interface perturbation
"
If we perturb an interface parameter ( 9
"
O. $ T O % ), equation (6)
yields (we have kept the same notations as for the case of a velocity perturbation)
10
M
M
M
A
A
A
O. $UT O % 9" 9 9 9 M
M
M
M
M
M
M
M
M A M 9"
A
9" M M M
M M M
M
M A
9 9 M
M
M
M M
AM M M 9" M M A M M M 9"
M M A M M M 9" Moreover, A I A , where A is the traveltime between points L and L and thus only
> different from zero. Besides, if we denote by - L a function equal to
> I and > > I are
%
one if L in on and zero otherwise, we obtain
M A 9 M A 9" M A 9 M A 9"
M M
M M
M M
M M
M A 9 M A 9 M M
M M M A - M 9. M A - M 9 M M
M M M
M L
L
%
%
M A M 9 M A M 9" M A M 9"
M A M 9"
M M M
M M M
M M M
M M M
M A M 9 M A M 9 M M M
M M M M
M
M
9. M A M - M 9. "
A - M
M
M
M
M M M
M L
L
%
%
M A A M M
M
9"
9"
M
M
M
M
M
M
M - % L 9 (8)
All the derivatives required can also be found in Appendix A. Here again, we have $ T
O" equations to recover $ T O" unknowns. Moreover, we notice that the coefficients of the
"
linear system are the same as the ones for a velocity perturbation, as soon as we deal with the
same ray. Consequently, these coefficients have to be calculated for each velocity and interface
model parameter, but only once per ray, whereas the right hand side of the linear system has to be
calculated for each ray and for each model parameter.
Solution of the linearized inverse problem
Once we have all the necessary operators at our disposal, we solve the linearized inverse problem thanks to a conjugate gradient algorithm. From a computational point of view, we decide to
11
store the Hessian matrix (see formula (4)) rather than the Jacobian, in order to avoid the storage of
four huge Jacobians (we usually invert tens of thousands of data in order to recover thousands of
parameters !) (see also Chauvier et al., 2000, this volume, for a discussion of this subject).
AN EXPERIMENT ON SYNTHETIC DATA
0.000
1.307
9.400
X
1.044
13.000
7.000
1.000
-5.000
13.000
11.000
9.000
7.000
5.000
3.000
1.000
-1.000
Z (km)
0.109
0.172
0.235
0.299
0.362
0.425
0.488
0.551
0.614
0.677
0.740
0.803
0.866
0.929
0.992
1.055
1.118
1.181
1.244
-3.000
-5.000
In order to validate the calculation, we perform an experiment on synthetic data. For this purpose, we try to invert data associated with the btert interface (see Figure 3) of the CASSIS model
(Jurado and Sinoquet, 1996).
0.109
X
0.509
2.089
0.908
3.133
1.307
4.178
6.267
9.400
6.267
0.000
5.222
3.133
Z
0.109
Y
7.311
0.509
8.356
0.908
1.307
Y
Z
Fig. 3 View of btert, from the CASSIS model : this interface is represented by a 2D cubic B-spline, and is defined by
parameters.
This interface is defined by , ( O parameters and, compared to the original CASSIS
configuration, has been raised by 500 meters to avoid triplications in the data that can hamper the
convergence of the Gauss-Newton algorithm (see Renard and Lailly, 1999). We use a constant 3D
layer velocity, defined by 023 , .
; 11200 multi-offset traveltimes organized in CMPs, these CMPs being
For our experiment, we use
distributed on a regular grid that covers the whole surface of the model. The source-receiver pairs
associated with one CMP are aligned in the direction. The CMP sampling intervals are 375 m
and 195 m for the and directions, respectively, and the offset sampling interval is 400 m, with
offsets ranging from 0 to 1600 m.
+ Moreover, we choose constant uncertainties ( ) and correlation lengths ( ' ) to perform our experiment. Again, our goal here is to validate our implementation of the
new formulation of the reflection traveltime tomography we have proposed : we want to check that
we are able to minimize the new non linear objective function and to converge to a satisfactory
12
solution.
Finally, we choose constant initial models ( K , and K ) and use a progressive
;
decrease of the regularization weights, as suggested in Renard and Lailly, 1999.
The results of the inversion show very low misfits : RMS= , O and Max= O O ; we
have managed to minimize the new objective function.
This experiment has thus allowed us to validate our implementation. Of course, we still do not
know how to make a proper choice of the different parameters that appear in our objective function.
Finding constructive answers to this basic question is the goal of our future work.
CONCLUSION
Reflection traveltime tomography is attractive for an accurate determination of the velocity
distribution within the subsurface. Nevertheless, the solution of this inverse problem, as it is classically formulated, a priori depends on the spatial sampling of the kinematic data, and may be, in
turn, at least in principle, corrupted by numerical artifacts.
In quest of an intrinsic solution to the tomographic inverse problem, we have proposed a new
objective function associated with a continuum problem taking the spatial correlations between
kinematic data into account. This new objective function takes into account not only the observed
traveltimes but also the derivatives of these observed traveltimes w.r.t. the CMP coordinates and
the offset. This new objective function makes use of four possibly spatially varying parameters : an
uncertainty and three correlation lengths for each traveltime. Adapting the solution of the associated linearized inverse problem to this formulation is the most delicate point. Indeed, it requires the
computation of three new Jacobians, which amounts to computing ray perturbations. This cumbersome calculation has been described in detail. Finally, we have performed a simple experiment
on synthetic data, which has led us to very satisfactory results (i.e. very low final misfits) : it has
allowed us to validate our implementation.
Finally, how to make a proper choice, in practice, of the different parameters appearing in the
new objective function, and how to discretize the kinematic information, remain the questions to
answer. This will be the next stage in our study. The answers will, of course, depend on the data.
ACKNOWLEDGEMENTS
This research was carried out as part of the Kinematic Inversion Methods consortium project
(KIM). The authors hereby acknowledge the support provided by the sponsors of the project.
They would also like to thank Fabrice Jurado and Delphine Sinoquet for their invaluable help and
advice, Pr M. Cuer, from Montpellier II University, for many fruitful discussions, and Andreas
Ehinger for his critical review of this paper.
REFERENCES
Adler, F., 1996, Reflection tomography from prestack migrated images : Ph.D. thesis, Universit´e
de Pau et des Pays de l’Adour.
13
Chauvier, L., Masson, R., and Sinoquet, D., 2000, Implementation of an efficient preconditioned
conjugate gradient in jerry : KIM Annual Report, Institut Franc¸ais du P´etrole, Rueil-Malmaison,
France.
Clarke, R. A., 1996, Two point raytracing in 3D blocky models : KIM Annual Report, Institut
Franc¸ais du P´etrole, Pau, France.
Delprat-Jannaud, F., and Lailly, P., 1993, Ill-posed and well-posed formulations of the reflection
travel time tomography problem : Journal of Geophysical Research, 98, 6589–6605.
Duijndam, A., 1988, Bayesian estimation in seismic inversion. Part I : principles : Geophysical
Prospecting, 36, 878–898.
Jurado, F., and Sinoquet, D., 1996, Application of reflection tomography on the 3D synthetic
dataset CASSIS : KIM Annual Report, Institut Franc¸ais du P´etrole, Pau, France.
Jurado, F., Sinoquet, D., and Lailly, P., 1996, Jerry : a 3D reflection tomography designed for
complex structures : KIM Annual Report, Institut Franc¸ais du P´etrole, Pau, France.
Jurado, F., Sinoquet, D., Ehinger, A., Lailly, P., and Monesti´e, C., 1997, Improvements in the reflection tomography software jerry : KIM Annual Report, Institut Franc¸ais du P´etrole, Pau, France.
Jurado, F., Lailly, P., and Ehinger, A., 1998, Fast 3D two-point raytracing for traveltime tomography : Proceedings of SPIE, Mathematical Methods in Geophysical Imaging V, 3453, 70–81.
Renard, F., and Lailly, P., 1999, Robust and accurate determination of complex velocity/depth
models by reflection travel time tomography : KIM Annual Report, Institut Franc¸ais du P´etrole,
Rueil-Malmaison, France.
Tarantola, A., and Valette, B., 1982, Inverse Problems
physics, 50, 159–170.
Quest for information : Journal of Geo-
APPENDIX A
The second derivatives in jerry
Our aim is here to calculate the second derivatives appearing in equations (7) and (8). Rays
are exactly straight lines if the velocity is homogeneous. For variant velocities, we approximate
"
the
constant slowness by mean of the slowness at points L K , L K and K (the mid-point of
L K L K % ) (cf. Figure 1).
O" O
K / K L and the traveltime
AK
O
/ K / L K %
between two successive impact points is given by
K K T K K T K = K T K K - K . According to equations (7) and (8), we have to
In the same way, we also have A K <
<
calculate, for , , and , ,
M AK
M- K M K
M
M M
M
O M M M M - K M M K M K M - K K M -M K K K
K K
K K
K K
K K
M
M
M
M
M
M
M- K
K
K
K
K
K
K
A
K
K
M K M K M K M K - M K M K M K M K M K M K AK K- K 14
,
M K
M K M K
M K A M K -M K M K - K
M A K M - K M K M K M K M K M K M K
M K M
M K M K
M AK M AK M K
M K M K M K M K O" M/
M K
T / K M T
M AK M M- K M
O / L K K
M K M- K
M K
M K M K K M - K M K M K M - K M - K K M K M K K M K M K
K
M
Of course, we need first to calculate the expressions of the first derivatives of the traveltimes with
respect to the impact point coordinates and the model parameters. We refer the reader to Jurado et
al, 1998, for these detailed calculations. We here just recall the main results
M AK
M AK
K
M- K K
KM M K M M / B % M - K K T K - K
K
K
We now calculate the required second derivatives
M K
M K M K O
"
O M / M / M
M
T /H K T / K
M
M
M /
/
/
O
/ L K M B M B T / L K M M
M K
M K - M K M
O
K
M K M K M M
K
M
O
K
M K M K M
O M / M M
B %
M
M K - K
9
K T K - K
T K T K K T K - K
O M / M / O M / M
M M
/ K /0 K M
T / K M- K
M - K
M K M K T M KM K 9
K
- K T K T K K T K M K - K M K
M K
M K - M K T M - K M K "
OM/ M/ O M / T / K M
T / K M M
M
M/
M/
M /
O
%
/
D L K M B M B T / L K M M B
M/
M/
M /
O
O
K O
M K M K
M
M M M
T / /
0
/
K
K
K
15
Finally, we have to calculate > >
54 4 M/
7
M 12 6 123
and > > > 123 . We easily obtain
7
7
M!
M!
M 6 ! ! 9 9 ! 6 M ! 9 7
854 6 9 "9
9 <;
and thus
4 54 58 4 7
M /
7 6 9
M M 023
6 9 ;
7
! 6 ! MM! 9 %
9
7
M! 6 M! !
M! 6! 7M! 9
9 M
M 7 9 9 M 7M %
M ! !
! 6 MM! MM! 9 %
9 ! 6 M
9
9
7
M ! 7 M ! 9 9 ! 6 ! M ! 9 %
9 ! 6 M
M
M
O
Consequently, we are now able to calculate without any difficulties expressions from to M ! 6 ! 7!
"
9 9 M
79 M! 6 M! !
"
9 9 M
" M ! M 7 M ! 9
9
9
9 M 6 !
M
M A K 7 M - K M K7
M K 7
M KM/6 9 M K M/6 9 - K M KM/6 9
M K 7
M A K 7 M - K M K 7
M K M / 6 9 M K M / 6 9 - K M K M / 6 9
The next expressions to calculate are
A K / O , we have
7
7
M K7
O
O
O
O
M / 6 9 T / L K ! 6 ! ! 9 B T / L K ! 6 ! ! 9 B ! T
If -
and thus
M
K
7
O
M KM/6 9 O
T / L
O O
T , / K
K
"
M ! !
6M
" M ! !
6M
7
! 9
!
7 B T / L K ! 9
T /0 K
7
A K / O , we have
7
7
M K7
O
O
O
O
!
!
!
!
!
M / 6 9 T / L K 6 9 B T / L K 6 ! 9 B T
If and thus
M
7
K
M K M / 6 9 O
O
T / L
O
T , / O
K
K
"
M ! !
6M
" M !
7
7
! 9
! ! !
7 B T 0/ L K 6 9 B
6! ! 9 !
M
T 0
/ K 16
! ! !
/ K 6 9 6! !
! 6 !
7
O
M/
9B M B %
7
! 9 M M / %
7
! 6 ! ! 9 / K O
M/
M B 7
%
M/
6! ! 9 M %
7
M ! 6 ! ! 9
M ! 6 ! 7!
9' M
9
M
with
We are now able to calculate and
It now remains to calculate, for
M K M
M K M K K we easily obtain
M K
M KM K and finally
9
<
K
M K M K M K
M KM M!
"
9 + 9 ' M 7M!
9
!
!
9 @ 6
M
M
and
As
7
!
M
!
9 + ! 6 M
9
"
9
M ! K !
K M
M !
"
' 9 ' M K ! K M!
K M K 9
9
M!
! K M K %
M!
M!
+ M K M K &%
9
M !
+ ! K M K &%
M K
M ! ! < K <
K M KM M
Thus, and are known : we now have at our disposal all the second derivatives necessary to
solve the linear systems (7) and (8).
17