TRAITEMENT AUTOMATIQUE DES SONDAGES ELECTRIQUES AUTOMATIC PROCESSING OF ELECTRICAL SOUNDINGS* (1970)

A bstract 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: the field measurements are automatically processed, without loss of information, into more easily usable data; 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 h i and resistivities ρ i ). They generally started by calculating the Stefanesco (or a similar) kernel function, from the integral equation of the apparent resistivity: image where r is the distance between the current source and the observation point, S 0 the Stefanesco function, ρ( z ) the resistivity as a function of the depth z, J 1 the Bessel function of order 1 and λ the integration variable. Thicknesses and resistivities had then to be deduced from S 0 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: image where K 1 is the modified Bessel function of order I. Using dimensionless variables t = r/2h 0 and y(t) = ζ (r)/ρ 1 and subdividing the earth into layers of equal thicknesses h 0 (highest common factor of the thicknesses h i ), ø becomes an even periodic function (period 2π) and the integral takes the form: image 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 R i = h i .ρ i –as function of the longitudinal conductance C i = h i /ρ 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: image minimum with inequalities as constraints: image where t l are the values of t for the sequence of measured resistivities and p l 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: image where: image 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 proc