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 $*' > ( $?+ > ( [email protected] > 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>[email protected] 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 ; >[email protected] >[email protected] and >[email protected] > , 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 >[email protected] > '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

© Copyright 2020