Anton G. Kolonin, Institute of Mathematics, Siberian Branch of Russian Academy of Science

The hardest problem of seismic inversion within diffraction approach is that most of techniques are based on assumption of weak velocity anomalies in geological medium. It is tricky to deal with contrast velocity anomalies and distinct areas of velocity variations from zones of absorption and dissipation of seismic energy. In this paper attempt is made to provide simply way to perform seismic wave inversion targeted on localization of solitary scatterers or emitters of seismic energy. Inversion approach assumes that absorption is present and plays its role together with velocity in forming complex dissipation factor with active (absorption) and reactive (delay) components. Reconstruction of distribution of this factor in geological medium considered being the main task of seismic diffraction tomography.

Seismic, inversion, tomography, diffraction, emission, dissipation

Inversion of seismic field is studied in details since pioneer works of Petrashen and Nahamkin [5] and Claerbout [2]. Idea to use seismic inversion within diffraction tomography approach for localization of geological heterogeneities was strongly promoted by Devaney [3], Wu and Toksos [8] and developed by many other authors. Most of approaches developed for absolutely elastic medium without of absorption or decay accounting and within Born's approximation of weak scattering. At the same time there are many geological objects that are not weak but strong scatterers with absorption features. For example, hydrothermal ore deposits, kimberlite intrusions and karst objects are contrast scatterers with anomalous absorption properties. Those objects may be small enough to be investigated using traditional ray-based techniques and diffraction tomography approach based on seismic wave field inversion is necessary to deal with them. Interesting results of absorption and decay of seismic wave studying are presented by Phaizullin and Shapiro [6] and Ampilov [1]. It is shown in [6] that within most of actual conditions absorption may be described by exponential function. Those results are used in non-diffractional ray tomography technique for amplitudes of seismic waves processing and absorption sections reconstruction developed by author [7]. At the same time, results of physical modeling carried out by author earlier [4] showed that velocity and absorption anomalies participates together in forming complex diffractional wave field. Actually, it is not easy to distinct wave filed anomalies caused by variations of velocity and by changes of absorption factor within scatterer. This paper introduces approach that incorporates absorption or decay factor of seismic wave propagation and diffraction applied to localization of absorbent or complex (absorbent and decelerating) scatterers.

Side result is technique of inversion of seismic emission field with accounting of absorption. Problem of location of seismic emitter is related to mine-holes or boreholes discovering or their orientation control during drilling process, centers of seismic danger control and earthquakes or mine-shocks prognosis.

Let us consider period of seismic wave (seconds) with frequency (Hertz). Velocity (meters/seconds) will me used as decay (seconds/meters) bounded with wavelength (meters).

, (1)

Actually, delay could be named by slowness or relative delay or specific delay. But term compact "delay" will be used in the rest part of this paper.

Amplitude of seismic oscillations excited at source with initial amplitude will depend on distance from source and time since moment of excitation .

Thus, absorption factor (1/meters) will be expressed as

(2)

and decay factor (1/seconds) will be expressed as

. (3)

Relation between (2) and (3) is clear by (1) as .

Decrement of decay

(4)

could be expressed as value of decay (seconds/meters)

(5)

extensively used in the rest of this paper.

When considering plane waves and fixed frequency, amplitude of seismic oscillations in the source may be written as below

(6)

in assumption that initial phase is known. Let us introduce time-dependent function of oscillations in source

, (7)

function of oscillations at any distance of given source in any time

(8)

and function of oscillations synchronized with source at any distance of it.

  (9)

In (9) we introduced value of complex dissipation that incorporates time active (decay or attenuation) and reactive (delay or slowness) components.

(10)

Let us express change of complex function of plane wave passing from point 1 to point 2 (Figure 1a)

(11)

on the path of length

.

Relative change of wave function could be expressed as

(12)

and complex dissipation factor (CDF) could be introduced as

. (13)

Delay and decrement of decay could be extracted from CDF as below.

(14)

When and , domain of definition of CDF is .

When (), CDF value turns into .

For absolutely elastic medium () with no delay (, ) complex dissipation and CDF value turns into as well. For absolutely elastic medium () with real values of delay CDF modulus turns into .

Let us consider any wave scattered on given area of medium with size (Figure 1b). We may introduce original wave field before scattering as

and immediately after scattering as

. (15)

Than, change of wave function that happens during scattering could be expressed as .

Let us consider that average values of and in medium surrounding given area are known. Thus, we introduce average value of and anomaly variation of in area of size .

(16)

Having that, we rewrite (13) as

(17)

and extract anomaly variations of components similar to (14) using known value of .

From (5) we express variation of decay and of its decrement as below.

(18)

Let us consider harmonic signal with given frequency where is number of harmonic. Signal is excited in source as initial wave function.

(19)

It is considered that average values of complex dissipation , delay and decay are frequency-independent. Thus, we express value of wave field registered at distance of source as below.

(20)

Now, we introduce receiver of oscillations identified by and source identified by . Length of path from source to receiver is and wave field within ray approach is

 . (21)

Let us assume that scattering happens once for each area of medium (Figure 1b). We perform discretization of entire medium within Cartesian coordinate system indexing by and per each axis (Figure 1c). Two-dimensional discretization is proposed for simplicity. Accordingly, cylindrical waves should be assumed in further considerations. But, it is not big complication to consider three-dimensional discretization and spherical waves. This way, we introduce wave single-scattering field as sum of multiple single scattering acts performed once by each area with size surrounding node of discretization matrix.

(22)

Full path of single-scattered wave will be equal to and angle of scattering will be expressed as . Factor of indicatrix of scattering could be expressed as

and assumed to be

if .

Actually, discretization, limited aperture of summation and replacing spherical waves with cylindrical distort actual field. So distortion recovery factor could be used to recover distortions

(23)

and complete (with recovered distortion) field could be obtained,

(24)

Finally, let us assume distribution of variations of within investigated medium and write formula for field scattered on those variations.

(25)

To perform inversion, let us replace , with , (Figure 1c) and consider full path of single-scattered wave as . Inversion performed within approach of [1] gives inverted wave field as

(26)

Having assumed assuming anomaly (16) present in , within area of size and having (23) and (25) accounted, we obtain scattered wave field.

(27)

Let us rewrite (17) to be used in given case.

(28)

Substitution of (26) into (27) gives final formula connecting CDF distribution in medium with registered wave function.

(29)

Let us introduce two complex coefficients independent on CDF variations with help of (10).

(30)

(31)

Having that, we can simplify (29) using (30), (31).

(32)

Unknown amplitude of signal excited in source per each frequency harmonic may be evaluated having inversed field in given number of discrete values within discretization byand numbers of discretization nodes.

(33)

Value of CDF desired as characteristic of complex dissipation variations could be obtained as

(34)

or in another form quite similar to (13) and (17)

. (35)

with definition of intermediate complex functions for "difference" field

(36)

and "basic" field

. (37)

Assuming linearity of all transformations above, we could derive alternative interpretation of described inversion schema from (26), (34).

(38)

This could be interpreted as prior determination of difference field and its posterior inversion and normalization.

For absolutely elastic medium () we obtain that

(39)

Attempt to simplify inversion process could be done as follows. We may consider from (25) and (17) that

(40)

(41)

Let us assume that anomaly variation of CDF is present in one area surrounding only. In other areas surrounding we have if and if . Thus, we discover that

(42)

and, substituting (37),

. (43)

It is interesting to compare (43) with (13), (17) and (35).

If assume that CDF is frequency independent () and variations of delay are relatively small (,) we may use results of inversions within number of different frequency harmonics to obtain average CDF value.

(44)

Interesting partial case when simplified forms of derived formulas could be used is reconstruction of emission sources. That means, only distances from hypothetical emitter of seismic waves to receiver are considered to reconstruct initial emitter oscillations per given frequency harmonic. Wave function of oscillations registered by receiver and simplified inversion field are given below.

(45)

(46)

More simplified case may be obtained in assumption that there is only one unknown emitter in whole medium, i.e. if and if . This case,

(47)

and

. (49)

Latter formula may be simplified in case of absolutely elastic medium .

. (50)

Inversion was carried out for through-passing waves scattered on strongly absorbent heterogeneity () imbedded into weakly absorbent medium () as shown at Figure 2a. Model data was calculated using (40). Results of inversion made by (26) and (44) shown as average modulus of inversed field (Figure 2b) and average phase delay of inversed field (Figure 2c). It is clear that heterogeneity well resolved using modulus (i.e. amplitude) because of active component of complex dissipation presented by absorption. At the same time, small false anomalies of phase delay are present as it is caused by reactive (delay) component of complex dissipation. That shows, that two components of complex dissipation are not resolved well one from another.

Another experiment was carried out for reverse scattered field formed on boundary with complex relief and several local extensive scatterers above the boundary (Figure 3b). Boundary and scatterers were modeled as areas of high value of CDF modulus; i.e. only active component was present. Original wave field shown at Figure 1a was calculated using (40) for impulse signal. Impulse signal was converted to frequency domain in range of 50-250 Hz. Spectral presentation of wave field was inversed by (26) and (44) using different steps of frequency discretization. Distribution of CDF modulus (Figure 3c) obtained with step equals to 2.5 Hz shows that reconstruction of boundary and solitary refractors is well. It works well despite that used formulas are derived in assumption of single-time scattering when only one scatterer is actually present. At the same time, distribution of CDF modulus obtained with frequency step equals to 10 Hz doesn't allow even to resolve boundary properly.

Another sample is related to problem of borehole orientation control during drilling process. Model of emission data that might be produced by drill bit is obtained on two receiver profiles on surface (Figure 4a). Data gathered on each profile was inversed using (46). This inversion allows determining position of emission source within plane of each profile. Because the problem is three-dimensional, it is needed to take together results obtained by several profiles (e.g. two profiles in given case) to perform spatial determination of emitter location (Figure 4b).

Single scattering approach could be used for localization of even non-weak but contrast heterogeneities with anomaly values of absorption and complex dissipation. Two basic forms of use are present - inversion of straightforward-scattered field applicable to cross-hole measurements or vertical seismic profiling and inversion of backward-scattered field applicable to surface measurements of reflection data and vertical seismic profiling as well.

Complex dissipation factor introduced gives a good opportunity to abstract inversion process from nature of scatterer and use obtained complex value of inversed field to resolve active and reactive components of dissipation related to velocity or decay variations. At the same time, this resolution is not easy and may be incomplete especially in case of limited amount of data and insufficient spatial and spectral discreteness.

Natural seismic sources, producers of technological noise peculiar to mining or drilling processes may be considered as seismic emitters. Location of those emitters is possible within three-step process. First, linear measurements of seismic emission data by surface profiles, wells or mine holes should be done. Next, inversion should be done to obtain integral sections of relative seismic emission power. Finally, results of several two-dimensional inversions should be taken together to perform three-dimensional interpretation.

Author is very grateful to Dr. Sergey V. Goldin for helpful discussions and useful critique. This research was carried out in Institute of Mathematics, Siberian Branch of Russian Academy of Science under supervision of Dr. Mikhail M. Lavrentiev.