A 2-D (x,z) isotropic medium is assumed. The velocity model is composed of a sequence of layers separated by boundaries consisting of linked linear segments of arbitrary dip. Layer boundaries must cross the model from left to right. Layer thicknesses may be reduced to zero to model pinchouts or isolated bodies. The velocity within a layer is defined by velocity values specified at arbitrary x-coordinates along the top and bottom of the layer. The x-coordinates at which layer boundaries and upper and lower velocities are specified can be completely general and independent within and between layers. Velocity discontinuities across layer boundaries are allowed but not required. For the purposes of ray tracing, the model is automatically broken up into an irregular network of trapezoids, each with dipping upper and lower boundaries and vertical left and right sides. The velocities at the four corners of the trapezoid are used to interpolate a velocity field within the trapezoid so that the velocity varies linearly along its four sides. Therefore, horizontal as well as vertical velocity gradients may exist within a trapezoid. A simulation of smooth layer boundaries is possible in which the incident and emergent ray angles are calculated using the slope of the smoothed boundary.
The source(s) may be positioned anywhere in the model and rays may be directed any angle. The receivers are always assumed to be at the top of the model. Both P- and S-wave propagation can be considered including (multiple) conversions. A unique Poisson's ratio may be assigned to each trapezoid of the model. Refracted, reflected and head waves may be traced, each possibly containing multiple and/or surface reflections and conversions. Ray take-off angles are determined automatically by the program for those ray groups specified by the user using an iterative shooting/bisection search mode. Ray tracing is performed by numerically solving the ray tracing equations for 2-D media, a pair of first order ODE's, using a Runge Kutta method. The ray step length is automatically adjusted at each step to maximize efficiency while maintaining accuracy. Travel times are calculated by numerical integration along ray paths using the trapezoidal rule. A plot of the model and all rays traced may be produced along with a plot of reduced travel time versus distance for the observed and calculated data.
The partial derivatives of travel time with respect to those model parameters selected for adjustment are calculated analytically during ray tracing; these parameters include velocities and the vertical position of boundary nodes. The travel times correspond to any ray paths which can be traced through the model, being either first or later arrivals. The travel time residuals with respect to the observed data are also calculated. The travel times and partial derivatives are linearly interpolated to the observed seismogram locations since two-point ray tracing is not required. The partial derivatives and travel time residuals are output and used later as input to the program DMPLSTSQR which updates the model parameters by applying the method of damped least-squares to the linearized inverse problem.
The model parameterization is well suited to the inversion of refraction/reflection data since realistic earth models can be represented by a minimum number of model parameters, i.e.,. the number and position of parameters specifying each layer can be adapted to the data's subsurface ray coverage. Layer boundaries, including the surface, may be horizontal (one parameter) or consist of numerous straight line segments. A layer may have a constant velocity (one parameter) or the velocity structure may be defined by many upper and lower layer velocity points. Different velocity points may be specified above and below a layer boundary if a velocity discontinuity is required across the boundary, or a single row of velocity points may be specified if an interface with no discontinuity is needed. The vertical velocity gradient or layer thickness may be fixed in all or part of a layer during the inversion if there is insufficient ray coverage to independently determine an upper and lower layer velocity or layer thickness.