Smoothing 2-D model HESS for Kirchhoff migrations

This is an example how to fit the gridded velocities and the gridded interfaces by a model with interfaces, suitable for ray tracing.

We take the gridded velocities and interfaces of the data set 'hess.dat', provided by Amerada Hess, and build the model suitable for ray tracing.

Data and model parametrization

File 'hess.dat' has the form of a file with input data for the MODEL package. The dimensions of the model are 40200 times 21000 feet. Two interfaces in the model are discretized on the grid of 64 points with spacing 200 feet. Background velocities are discretized on the grid of 202*106 points with spacing 200 feet*200 feet, velocity inside a salt body is constant.
Model to be smoothed.
Fig. 1: Model 'hess.dat'.

We can see that the horizontal variation of the velocity field is less pronounced than the vertical variation, and that the vertical variation is greater in the shallow parts of the model than in its deeper parts. The interfaces are quite smooth, but the edges are sharp. According to these observations, we manually choose new parametrization of the model. Let's emphasize, that the parametrization is based on authors opinion, and is expected to be optimized during the inversion. The parametrization is as follows: for the interfaces we take irregular B-spline grid of 28 points with 26 points regularly spaced each 500 feet in the region of the salt body and 2 points at the sides of the model. For the background velocity, we take horizontally regular grid of 17 points with spacing 2500 feet and vertically irregular grid of 31 points with spacing from 200 to 800 feet. The parametrization may be found in file 'hes1-mod.dat', which is used as a starting model for the inversion.

For interfaces we have the data grid of 2*64=128 points with spacing 200 feet and the irregular B-spline grid of 2*28=56 points. For the velocities we have the data grid of 202*106=21412 points with spacing 200 feet. The irregular B-spline grid of 17*31=527 points is used for the background. The velocity of the salt body is constant (1 B-spline gridpoint). Total number of data points is thus 128+21412=21540 points, total number of B-spline points is 56+527+1=584 points. Program 'gmdmgmt.for' then requires 584*21540+21540+584*584=12941956 storage locations in array RAM and we may perform the calculations with MRAM=14000000, which is a value suitable for a computer with 64 Mbytes of memory.

Notes on history file ' hes-inv.h'

Choosing the maximum Lyapunov exponent

To estimate the maximum sum of travel times from the source and receiver, we compute average slowness 0.120611 E-3 s/feet. If we consider maximum sum of ray paths about 50000 feet, we arrive at the estimation, that the maximum sum of travel times from the source and receiver is limited by Tmax = 6 s. According to the equation (2) of the paper 'mar-inv.htm', we choose the maximum Lyapunov exponent lambda of the new model as
lambda .LE. 2/Tmax = 0.333 /s . (1)
The product of the numbers of arrivals from a source and a receiver then should not exceed exp(2), i.e., should not exceed 7.

Choosing the maximum Sobolev norm of the model

According to the equation (7) of the paper 'mar-inv.htm', making the same approximation as in equation (10) of the paper, we arrive at
||v|| .LE. sqrt(8/3) * ( 4 * (1+ln(2)) / Tmax )2 / average(v) , (2)
where ||v|| is the 2-D Sobolev norm of the velocity field given by file 'sob22.dat'. For average velocity average(v) = 8074.73 feet/s of the background, we wish the maximum Sobolev norm of the velocity field in the model
||v||max = 2.6 E-4 (s*feet)-1 . (3)

In this example, we try to apply the minimum reasonable smoothing of the interfaces. We thus do not prescribe maximum Sobolev norm of the interfaces.

Inversion and smoothing

Objective function is taken in the form of
y = (|x-x0|/(1 foot))2 + (SOBMULx*||x||)2 + (|v-v0|/(1 foot/s))2 + (SOBMULv*||v||)2 , (4)
where |x-x0| is the standard deviation of the position of the interfaces in the model from the gridded interfaces in model 'hess.dat', ||x|| is the 2-D Sobolev norm of the curvature of the interfaces given by file 'sob22.dat', |v-v0| is the standard velocity deviation of the model from gridded velocities and ||v|| is the 2-D Sobolev norm of the velocity field given by file 'sob22.dat'.

We use unit given velocity deviation and unit given interface deviation in the selected objective function. The unit deviations are realized by proper selections of parameter ERRMUL, which should be chosen as a square root of the ratio of the whole model volume to the volume corresponding to a single data point of the considered data set.

The inversion of the interfaces and the inversion of the velocity field are mutually independent in this example. We can thus perform the inversion of the interfaces in the first step, and the inversion of the velocity in the second step.

(1) Inversion of the interfaces:

We choose SOBMULv=0 s*feet, and perform the inversion without smoothing of the velocity, which is not of our interest at the moment.

In this example, we try to apply the minimum reasonable smoothing of the interfaces. We thus start the inversion without smoothing of the interfaces, with SOBMULx=0 feet (the interfaces are smoothed just by the projection onto the B-splines). We obtain oscillations of both surfaces, which produces other salt body than that one given by data. In the next steps of inversion we thus include minimization of the curvature of the interfaces by increasing the value of SOBMULx. When we reach approximately the value of SOBMULx=2500 feet, the spurious salt body moves outside the model volume. We thus consider the value of SOBMULx=2500 feet as the value corresponding to the minimum reasonable smoothing of the interfaces. Whether the interfaces are sufficiently smooth will be tested by ray tracing.

We compute standard interfaces deviation |x-x0|=5.6 feet, stored in the file 'hes-sd1.out'. This deviation causes a travel time error dT approximately
dT = |x-x0| * |1/average(v) - 1/salt(v)| . (5)
For the average velocity average(v) = 8074.73 feet/s of the background and the velocity of the salt salt(v) = 14800 feet/s, we estimate the standard travel time deviation caused by the deviation of the interfaces dT = 3.15 E-4 s.

(2) Inversion of the velocity field:

We have selected the value of SOBMULx=2500 feet, which ensures proper position of the interfaces. We will keep the value during the velocity inversion.

(a) First iteration:
We first perform inversion with SOBMULv=0 s*feet (i.e., without smoothing), and obtain the standard velocity deviation of the model
|v-v0| = 13.94 feet/s (6)
stored in file 'hes-vd1.out'. The Sobolev norm of the velocity field smoothed just by the projection onto the B-splines
||v|| = 37 E-4 (s*feet)-1 , (7)
stored in file 'hes-vn1.out', is large, see (3).

(b) Second iteration:
We choose
SOBMULv = |v-v0| / ||v||max , (8)
as in the paper 'mar-inv.htm'. Then
SOBMULv = 13.94 / 2.6 E-4 s*feet = 53615 s*feet (9)
and we enter the rounded value of
SOBMULv = 50000 s*feet (10)
for the second iteration. The resulting Sobolev norm
||v|| = 2.9 E-4 (s*feet)-1 , (11)
of the velocity looks acceptable, see (3).

(c) Subsequent iterations:
We choose several values of SOBMULv for subsequent iterations in order to see the behavior of the model in dependance on SOBMULv. After each inversion we perform several tests, which are described below in the chapter "Notes on file 'hes-test.h'", and compute for each model the value of average Lyapunov exponent lambdaav and the maximum number of arrivals NUMmax from the source located in the left top corner of the model. We performed the following tests of SOBMULv:
model |v-v0| / (feet/s) ||v||  / ((s*feet)-1) lambdaav / (s-1) NUMmax
desired values   2.6 E-4 0.333 7
model 'hess.dat'     0.423 11
SOBMULv=0 13.94 37 E-4 0.380 7
SOBMULv=17500 15.19 4.9 E-4 0.345 6
SOBMULv=35000 17.41 3.5 E-4 0.316 5
SOBMULv=50000 19.24 2.9 E-4 0.300 5
SOBMULv=70000 21.42 2.4 E-4 0.287 5

We can see, that for SOBMULv=0 s*feet we obtain the standard velocity deviation about 14 feet/s. Desired value of lambdaav=0.333 s-1 corresponds to the velocity deviation about 16 feet/s. Most of the velocity deviation is thus caused by the projection onto the B-splines, and we should try the inversion with finer parametrization of the model.

A figure of Fourier spectrum of velocity deviation should help us to decide, in which direction and how many times we should densify the parametrization of the model. We must keep in mind, that the parametrization must not be denser than the data.
Spectrum of velocity deviation.
Fig. 2: Fourier spectrum of velocity deviation in the model obtained with SOBMULv=50000 s*feet.

Finer B-spline parametrization of the model

We choose new parametrization of the model as follows: for the interfaces we keep the irregular B-spline grid of 28 points; for the background velocity we take horizontally regular grid of 26 points with spacing 1600 feet and vertically irregular grid of 58 points with spacing from 200 to 400 feet. The parametrization may be found in file 'hes2-mod.dat', which is used as a new starting model for the inversion.

We have now 56+1+26*58=1565 B-spline points, and the memory requirements are 1565*21540+21540+1565*1565=36180865 storage locations. This is too much and we thus cannot use all the data. We resample the data in horizontal direction and use the data grid of 67 points with spacing 600 feet horizontally and original 106 points with spacing 200 feet vertically. We have now 67*106+2*64=7230 data points and the memory requirements are 1565*7230+7230+1565*1565=13771405 storage locations.

We performed the following tests of SOBMULv with new starting model 'hes2-mod.dat':
model |v-v0| / (feet/s) ||v||  / ((s*feet)-1) lambdaav / (s-1) NUMmax
desired values   2.6 E-4 0.333 7
model 'hess.dat'     0.423 11
SOBMULv=1000 5.52 12.3 E-4 0.398 6
SOBMULv=10000 7.28 6.5 E-4 0.365 5
SOBMULv=17500 9.03 5.2 E-4 0.343 5
SOBMULv=25000 10.68 4.4 E-4 0.327 5
SOBMULv=35000 12.75 3.7 E-4 0.313 5

Note that with the new parametrization we are not able to perform the inversion with SOBMULv=0 s*feet, because we have now more background velocity B-spline gridpoints in the salt body, and the inversion becomes ill-conditioned.

Notes on history file 'hes-test.h'

History file 'hes-test.h' serves for numerical tests of model 'hes-mod.out' obtained after the inversion realized by history file 'hes-inv.h'. The figures displayed below are obtained by running 'hes-inv.h' with starting model 'hes2-mod.dat' and the values of SOBMULx=2500 feet and SOBMULv=25000 s*feet.

Statistical properties of the model

We discretize the velocity in the smoothed model 'hes-mod.out'. We calculate the gridded velocity deviation from the velocity in the model 'hess.dat', the gridded relative velocity deviation, the standard velocity deviation and the standard relative velocity deviation, the Fourier spectrum of the velocity deviation, the deviations of the interfaces, the relative deviations of the interfaces, and the standard deviation of the interfaces and standard relative deviation of the interfaces.
Velocity in the inverted model.
Fig. 3: Velocity in the smoothed model 'hes-mod.out'.

Velocity deviation.
Fig. 4: Velocity deviation.

Relative velocity deviation.
Fig. 5: Relative velocity deviation.

Standard velocity deviation stored in file 'hes-vd2.out' is 10.68 feet/s.
Standard relative velocity deviation stored in file 'hes-vr2.out' is 0.18 % (i.e. 1.8 per mille).
Spectrum of velocity deviation.
Fig. 6: Fourier spectrum of velocity deviation.

Standard deviation of the interfaces stored in file 'hes-sd2.out' is 5.6 feet.
Standard relative deviation of the interfaces stored in file 'hes-sr2.out' is 0.08 %.

Average Lyapunov exponent for the model

Directional Lyapunov exponents and the average Lyapunov exponent for the 2-D model are calculated by program 'modle2d.for'.
Lyapunov exponents.
Fig. 7: Angular dependence of directional Lyapunov exponents and the average Lyapunov exponent.

Tests of ray tracing

Three tests of ray tracing of a refracted wave followed by interpolation between rays are performed using programs 'crt.for' and 'mtt.for'. The three sources are located at the top of the model in the corners and in the middle of the model.
Numbers of arrivals from the left source.
Fig. 8: Rays and numbers of arrivals for the source in the left top corner of the model.

Numbers of arrivals from the middle source.
Fig. 9: Rays and numbers of arrivals for the source in the middle of the top of the model.

Numbers of arrivals from the right source.
Fig. 10: Rays and numbers of arrivals for the source in the right top corner of the model.

References