TRAITEMENT AUTOMATIQUE DES SONDAGES ELECTRIQUES AUTOMATIC PROCESSING OF ELECTRICAL SOUNDINGS*

I. Introduction

In this section the problem is stated, its physical and mathematical difficulties are indicated, and the way the authors try to overcome them are briefly outlined.

Made up of a few measurements of limited accuracy, an electrical sounding does not define a unique solution for the variation of the earth resistivities, even in the case of an isotropic horizontal layering.

Interpretation (i.e. the determination of the true resistivities and thicknesses of the ground‐layers) requires, therefore, additional information drawn from various more or less reliable geological or other geophysical sources. The introduction of such information into an automatic processing is rather difficult; hence the authors developped a two‐stage procedure:

• a)   the field measurements are automatically processed, without loss of information, into more easily usable data;
• b)   some additional information is then introduced, permitting the determination of several geologically conceivable solutions.

The final interpretation remains with the geophysicist who has to adjust the results of the processing to all the specific conditions of his actual problem.

II. Principles of the procedure

In this section the fundamental idea of the procedure is given as well as an outline of its successive stages.

Since the early thirties, geophysicists have been working on direct methods of interpreting E.S. related to a tabular ground (sequence of parallel, homogeneous, isotropic layers of thicknesses hi and resistivities ρi). They generally started by calculating the Stefanesco (or a similar) kernel function, from the integral equation of the apparent resistivity:

where r is the distance between the current source and the observation point, S0 the Stefanesco function, ρ(z) the resistivity as a function of the depth z, J1 the Bessel function of order 1 and λ the integration variable. Thicknesses and resistivities had then to be deduced from S0 step by step. Unfortunately, it is difficult to perform automatically this type of procedure due to the rapid accumulation of the errors which originate in the experimental data that may lead to physically impossible results (e.g. negative thicknesses or resistivities) (II. 1).

The authors start from a different integral representation of the apparent resistivity:

where K1 is the modified Bessel function of order I. Using dimensionless variables t = r/2h0 and y(t)=ζ (r)/ρ1 and subdividing the earth into layers of equal thicknesses h0 (highest common factor of the thicknesses hi), ø becomes an even periodic function (period 2π) and the integral takes the form:

The advantage of this representation is due to the fact that its kernel ø (function of the resistivities of the layers), if positive or null, always yields a sequence of positive resistivities for all values of θ and thus a solution which is surely convenient physically, if not geologically (II.3). Besides, it can be proved that ø(θ) is the Fourier transform of the sequence of the electric images of the current source in the successive interfaces (II.4).

Thus, the main steps of the procedure are: a) determination of a non‐negative periodic, even function ø(θ) which satisfies in the best way the integral equation of apparent resistivity for the points where measurements were made; b) a Fourier transform gives the electric images from which, c) the resistivities are obtained. This sequence of resistivities is called the “comprehensive solution”; it includes all the information contained in the original E.S. diagram, even if its too great detail has no practical significance.

Simplification of the comprehensive solution leads to geologically conceivable distributions (h, ρ) called “particular solutions”. The smoothing is carried out through the Dar‐Zarrouk curve (Maillet 1947) which shows the variations of parameters (transverse resistance Ri= hi.ρi–as function of the longitudinal conductance Ci=hi/ρi) well suited to reflect the laws of electrical prospecting (principles of equivalence and suppression). Comprehensive and particular solutions help the geophysicist in making the final interpretation (II.5).

III. Computing methods

In this section the mathematical operations involved in processing the data are outlined.

The function ø(θ) is given by an integral equation; but taking into account the small number and the limited accuracy of the measurements, the determination of ø(θ) is performed by minimising the mean square of the weighted relative differences between the measured and the calculated apparent resistivities:

minimum with inequalities as constraints:

where tl are the values of t for the sequence of measured resistivities and pl are the weights chosen according to their estimated accuracy.

When the integral in the above expression is conveniently replaced by a finite sum, the problem of minimization becomes one known as quadratic programming. Moreover, the geophysicist may, if it is considered to be necessary, impose that the automatic solution keep close to a given distribution (h, ρ) (resulting for instance from a preliminary interpretation). If φ(θ) is the ø‐function corresponding to the fixed distribution, the quantity to minimize takes the form:

where:

The images are then calculated by Fourier transformation (III.2) and the resistivities are derived from the images through an algorithm almost identical to a procedure used in seismic prospecting (determination of the transmission coefficients) (III.3).

As for the presentation of the results, resorting to the Dar‐Zarrouk curve permits: a) to get a diagram somewhat similar to the E.S. curve (bilogarithmic scales coordinates: cumulative R and C) that is an already “smoothed” diagram where deeper layers show up less than superficial ones and b) to simplify the comprehensive solution.

In fact, in arithmetic scales (R versus C) the Dar‐Zarrouk curve consists of a many‐sided polygonal contour which múst be replaced by an “equivalent” contour having a smaller number of sides. Though manually possible, this operation is automatically performed and additional constraints (e.g. geological information concerning thicknesses and resistivities) can be introduced at this stage. At present, the constraint used is the number of layers (III.4).

Each solution (comprehensive and particular) is checked against the original data by calculating the E.S. diagrams corresponding to the distributions (thickness, resistivity) proposed. If the discrepancies are too large, the process is resumed (III.5).

IV. Examples

Several examples illustrate the procedure (IV). The first ones concern calculated E.S. diagrams, i.e. curves devoid of experimental errors and corresponding to a known distribution of resistivities and thicknesses (IV. 1).

Example I shows how an E.S. curve is sampled. Several distributions (thickness, resistivity) were found: one is similar to, others differ from, the original one, although all E.S. diagrams are alike and characteristic parameters (transverse resistance of resistive layers and longitudinal conductance of conductive layers) are well determined. Additional informations must be introduced by the interpreter to remove the indeterminacy (IV.1.1).

Examples 2 and 3 illustrate the principles of equivalence and suppression and give an idea of the sensitivity of the process, which seems accurate enough to make a correct distinction between calculated E.S. whose difference is less than what might be considered as significant in field curves (IV. 1.2 and IV. 1.3). The following example (number 4) concerns a multy‐layer case which cannot be correctly approximated by a much smaller number of layers. It indicates that the result of the processing reflects correctly the trend of the changes in resistivity with depth but that, without additional information, several equally satisfactory solutions can be obtained (IV. 1.4).

A second series of examples illustrates how the process behaves in presence of different kinds of errors on the original data (IV.2).

A few anomalous points inserted into a series of accurate values of resistivities cause no problem, since the automatic processing practically replaces the wrong values (example 5) by what they should be had the E.S. diagram not been wilfully disturbed (IV.2.1).

However, the procedure becomes less able to make a correct distinction, as the number of erroneous points increases. Weights must then be introduced, in order to determine the tolerance acceptable at each point as a function of its supposed accuracy. Example 6 shows how the weighting system used works (IV.2.2).

The foregoing examples concern E.S. which include anomalous points that might have been caused by erroneous measurements. Geological effects (dipping layers for instance) while continuing to give smooth curves might introduce anomalous curvatures in an E.S. Example 7 indicates that in such a case the automatic processing gives distributions (thicknesses, resistivities) whose E.S. diagrams differ from the original curve only where curvatures exceed the limit corresponding to a horizontal stratification (IV.2.3).

Numerous field diagrams have been processed (IV. 3). A first case (example 8) illustrates the various stages of the operation, chiefly the sampling of the E.S. (choice of the left cross, the weights and the resistivity of the substratum) and the selection of a solution, adapted from the automatic results (IV.3.1). The following examples (Nrs 9 and 10) show that electrical prospecting for deep seated layers can be usefully guided by the automatic processing of the E.S., even when difficult field conditions give original curves of low accuracy. A bore‐hole proved the automatic solution proposed for E.S. no 10, slightly modified by the interpreter, to be correct.