Method for measuring induced polarization and electrical resistivity of the subsoil using multiple transmitters transmitting simultaneously with CDMA coding

Abstract

The invention relates to a method for measuring the polarizability and resistivity of a subsoil comprising at least the following steps: - selection and reading of measurement data (2) from several transmitters injecting electrical signals into the subsoil and from at least one receiver capturing induced signals restored by the subsoil, the data forming a reception sequence from the receiver, the reception sequence being coded, - estimation of a spontaneous polarization (8) of the subsoil over the entire reception sequence, - spectral domain filtering of the estimated spontaneous polarization (10), - subtraction of the estimated spontaneous polarization (12), - decoding of the resistivities (14) and obtaining the measured resistivities as a function of the transmitter, and - decoding of the chargeabilities (16) and obtaining the measured chargeabilities as a function of the transmitter. Figure for abstract: figure 2

Description

Title of the invention: Method for measuring induced polarization and electrical resistivity of the subsoil using multiple transmitters transmitting simultaneously with CDMA coding

[0001] The invention relates to a method for determining the nature of a subsoil and more particularly for measuring the induced polarization and the electrical resistivity of a subsoil.

[0002] It is known, in order to determine the nature of a subsoil, to inject different types of signals into the subsoil in order to determine and quantify the reactions of the different elements composing this subsoil to the injection of the signals, which ultimately makes it possible to obtain certain physicochemical properties of the constituent materials of said subsoil. This information makes it possible to make projections regarding the exploitation of a place, for example in the context of regional planning.

[0003] For example, methods are known for determining the nature of a subsoil based on the emission of seismic waves. Measuring the propagation of waves in the subsoil makes it possible to determine the nature and structure of the subsoil.

[0004] It is also known to measure the electrical polarizability and resistivity of the subsoil by injecting into the latter, via a current transmitter connected to two electrodes (often denoted A and B), an electric current forming an emission signal. The different materials forming the subsoil will react in different ways, by attenuating the current more or less and by polarizing differently. A receiver connected to two so-called potential electrodes (often denoted M and N), measures the induced voltage resulting from the reaction of the subsoil to the emitted signal, and an analysis of the data thus obtained makes it possible to estimate the nature of the subsoil in an area between the transmitter and the receiver, in particular thanks to the induced polarization (PP) but also thanks to the electrical resistivity.

[0005] To obtain more complete data from a site, it is necessary to carry out a mapping of induced polarization and resistivity of the subsoil by means of several transmitters associated with several receivers distributed over the site. We can therefore speak of a device distributed with several transmitters / receivers. Traditionally, these measurements are carried out with successive transmitters, injecting their signals one after the other. However, to increase efficiency or to adapt to rapidly evolving processes, it may be useful to carry out simultaneous injections of the various transmitters. In this case, since the receivers receive the signals from the various transmitters at the same time, it is necessary to provide a code associated with each transmitter in order to be able to dissociate the signals coming from each of them at the receivers.

[0006] Determining the nature of a subsoil by measuring induced polarization and resistivity can pose disadvantages in terms of reliability of the measurements carried out due to the use of transmitters not regulated in current or even natural background noise (spontaneous polarization or "PS", magnetotelluric signals or "MT").

[0007] "Spontaneous polarization" means local temporal variations of the potential measured at the receiver having a frequency / significantly lower than the transmission frequency / ) ( / « / )) but with amplitudes often much higher than the signal received from the transmitter. The linear spontaneous polarization that is often discussed is an ideal case that does not exist, or only briefly: in reality, spontaneous polarization is always undulating ([Fig.3]). On the contrary, if / » / 0 we will rather speak of HF noise and this will be eliminated by averaging or "stacking".

[0008] The object of the invention is in particular to provide a method for measuring the electrical polarizability and resistivity of a subsoil, using data from several simultaneous transmitters, in order to measure in a single pass and in the most reliable manner possible the induced polarization and the resistivity associated with each of these transmitters.

[0009] To this end, the invention relates to a method for measuring the induced polarization and the resistivity of a subsoil comprising at least the following steps:

[0010] - selection and reading of measurement data from simultaneous injecting transmitters electrical signals coded in the subsoil and at least one receiver capturing an induced signal restored by the subsoil following the injections of said electrical signals by all the synchronous transmitters, the data forming a reception sequence from the receiver, the reception sequence being the superposition of coded potentials resulting from coding of the transmitters according to code division multiple access coding and forming different coding patterns configured to discriminate the signal emitted by each transmitter,

[0011] - estimation of a spontaneous polarization (8) of the subsoil at the receiver level, on the entire reception sequence,

[0012] - spectral domain filtering of the estimated spontaneous polarization (10),

[0013] - subtraction of the filtered estimated spontaneous polarization sequence from the sequence reception (12), - decoding, for each coding pattern processed in the reception sequence, of the subsoil resistance associated with this pattern for each transmitter,

[0014] - averaging of the resistances obtained for all the coding patterns processed in the reception sequence, for each transmitter used, and obtaining the average resistance then the apparent resistivity associated with each transmitter,

[0015] - decoding, for each coding pattern processed in the reception sequence, of partial chargeabilities and the overall chargeability of the subsoil associated with this pattern for each transmitter used (16), and

[0016] - averaging of the partial and global chargeabilities obtained for all the patterns of coding processed in the reception sequence, for each transmitter used, and obtaining the partial apparent chargeabilities and the global chargeability associated with each transmitter.

[0017] "Apparent resistivity" means the resistivity of a homogeneous subsoil which gives the same voltage as that measured at the level of potential electrodes for the same current injected to the current electrodes. The definition is the same for apparent chargeability. These so-called "apparent" values are in fact averages over a certain volume of soil between the transmitter and the receiver up to a depth equal to approximately half the transmitter-receiver distance.

[0018] Thus, and after selection and reading of the data from the receiver(s), the spontaneous polarization (i.e. the natural potential of the subsoil present between the receiving electrodes) is estimated and subtracted to retain only the potentials resulting from the stimulation induced by the transmitters. The cleaning of the reception sequence recorded by the receiver(s) is decisive in order to eliminate any signal not caused by the transmitters before decoding.

[0019] The decoding of the resistances and chargeabilities at the receiver level then makes it possible to obtain, for each transmitter, the apparent induced polarization induced by the injection of the electric current as well as the apparent electrical resistivity of the area investigated by the measuring device (area between the receiver and each transmitter). This determination finally makes it possible to estimate the nature of the subsoil in this same area.

[0020] According to other optional characteristics of the measuring method taken alone or in combination: - the measurement data recovered comes from current transmitters not regulated in current; - the induced polarization and resistivity are measured from several reception sequences from several receivers; - the coded sequences comprise several elementary patterns of the coded injection sequences, i.e. 24 ON / OFF emission levels in the case of three transmitters, and in which the data analysis is carried out on at least one of these coding patterns, and in which an estimation of the spontaneous polarization is carried out; - before the estimation of the spontaneous polarization, an estimation and a suppression of a continuous component and / or a slow drift of the current injected by each transmitter used; - the following steps are carried out during the spectral domain filtering step of the estimated spontaneous polarization (10): - calculation of a Fourier transform from the estimate of spontaneous polarization, - suppression of spectral components ranging from / o / 2 to 3 / 0 / 2, / 0 being the frequency of the signal from the transmitters resulting in a ripple frequency affecting the estimated spontaneous polarization, and - reconstruction of the suppressed spectral components by piecewise Hermite cubic interpolation, and - obtaining a spontaneous polarization filtered by calculating an inverse Fourier transform; - during the step of decoding the resistances (14) and the partial and global chargeabilities (16), a current matrix in which the continuous component has been subtracted, and a potential vector in which the filtered estimated spontaneous polarization has been subtracted are used; and - we determine a vector of chargeabilities from a vector of pseudoresistances observed during the time decrease.

[0021] The invention also relates to a method for determining the nature of a subsoil, comprising the following steps:

[0022] - injection of coded electrical signals into the subsoil via several simultaneous electric current transmitters, each transmitter injecting into the subsoil an electric signal comprising a code specific to the transmitter,

[0023] - capture, by at least one receiver distinct from the transmitters, of an induced signal returned by the subsoil following the injection of electrical signals, the received signal forming a reception sequence, and

[0024] - implementation of a method for measuring the induced polarization and the re subsoil sistivity according to the invention. Brief description of the figures

[0025] The invention will be better understood on reading the following description given solely by way of example and with reference to the appended drawings in which:

[0026] [Fig. 1] is a representation of the different stages of the measuring method according to the invention,

[0027] [Fig.2] is a representation of CDMA coding for three transmitters on one receiver ideal, and

[0028] [Fig.3] is an illustration of the estimation of spontaneous polarization by splines cubics followed by spectral filtering compared to classical polarization estimation spontaneous by the 3-point average, on synthetic induced resistivity-polarization data in the case of 4-second emission steps. Detailed description

[0029] We now refer to [Fig.l] schematically illustrating the different stages of a method for measuring the induced polarization and the electrical resistivity of a subsoil by the use of several simultaneous transmitters (commonly called "Tx", abbreviation used hereinafter) injecting an electric current (i.e. an electrical signal) coded into the subsoil, and at least one, but preferably several receivers (commonly called "Rx", abbreviation used hereinafter) capturing the electrical potentials or voltages restored by the subsoil following the injections of electric current from all the synchronous transmitters. The receiver, or receivers, capture(s) the reaction of the subsoil to the injection of electric current by all the simultaneous transmitters, as a function of the positioning on the ground of each Tx-Rx pair, more precisely as a function of the local electrical properties of said ground (resistivity, polarizability).

[0030] The invention focuses on the processing (a priori in deferred time) of the recordings of coded electrical currents injected simultaneously to the various Tx and of the recordings of the potentials received in reaction to the various Rx. There is therefore first of all a phase of selection and reading of the data acquired by the various simultaneous transmitters (injected currents) and by at least one receiver capturing an electrical potential signal (voltage) induced by the subsoil following the injections of electrical signals by all the transmitters, the data forming a reception sequence from the receiver, the reception sequence being the superposition of coded potentials resulting from the coding of the transmitters according to a code division multiple access (CDMA) coding, which makes it possible to discriminate the signal received from each transmitter.

[0031] First of all, and concerning the steps prior to this phase of selection and reading of data, we will address the question of the acquisition of these data. To do this, we have an acquisition device formed by several transmitters (device called "multi-Tx", abbreviation used subsequently) and by at least one receiver. These various elements are assumed to be "independent" in the sense that they are simply synchronized with each other by GPS, without requiring any wired or radio link to connect them (although such a link would not be detrimental to the invention, it would even allow our processing to be applied in real time).

[0032] Concerning the transmitters, each of them injects an electric current into the subsoil by emitting a coded square signal, each transmitter being assigned a specific code according to a code division multiple access coding (or CDMA coding for “Code Division Multiple Access” in English). This coding will be called “CDMA coding” afterwards.

[0033] CDMA coding is a coding of information from a transmitter with the aim of being able to multiplex data from different simultaneous transmitters and subsequently demultiplex them and thus find the data from each transmitter separately. The general purpose of multiplexing is therefore the transmission of several communications on the same channel (here acquisition time interval).

[0034] CDMA codes (for example the CDMA codes of [Fig.2]) are binary codes (i.e. successions of -1 and +1) whose length N (number of bits) is a power of 2 (N=2Ak). Among the large number of possible binary sequences (2AN in theory), only Nl = 2Ak-l possess the orthogonality properties required for a synchronous multiple transmitter process. These pseudo-random sequences are known as Walsh sequences.

[0035] Synchronous CDMA coding is generally referred to as "Synchronous CDMA". The message to be transmitted is directly multiplied by the binary code. On reception, the correlation of the signal with a replica of the transmitter code is calculated, which allows the message bits to be regenerated.

[0036] For the purposes of geophysics, the chosen code (succession of -1 and +1) is used to modulate an elementary emission current pattern, typically a square pattern with six steps ON+, 0, ON-, 0, ON+, 0 (i.e. succession of positive and negative square waves theoretically of the same intensity interspersed with signal interruption steps) to allow an acquisition of induced polarization (or PP, abbreviation used subsequently) followed by a compensation of spontaneous polarization (or PS, abbreviation used subsequently) "on three points" (method described later).

[0037] Within the framework of the invention, the following formalism can be adopted.

[0038] We start from a set of m simultaneous transmitters coded by CDMA coding. As explained above, the number of available codes is a power of 2. By eliminating the uniform codes, the number of transmitters is of the form m=2Ak-l=Al. In practice, k can be 2, 3 or 4, therefore m can be 3, 7 or 15. We will subsequently consider only the case with 3 transmitters (k=2, m=3, A=4), which assumes coding on 4 bits (N=m+l=2^k), i.e. 4 levels (see references 18, 20 and 22 in [Fig.2]). The case with 7 transmitters would be identical, with nevertheless 7 longer coded patterns (the coding being on 8 bits, therefore 8 levels).

[0039] In the following we will call "tip" (see reference 24 in [Fig.2]) an elementary pattern of the coded injection sequences. A tip is therefore made up of the minimal assembly of a certain number of ON+, ON- (such that ION- = -Ion+) or OFF (such that Iqff=0) steps, as discussed above, an assembly which differs according to the prospecting method. a. In the DC resistivity method without PP measurement, the traditional tip consists of two steps, ON+ and ON-, which represent a TO period of the emitted signal (this tip is used in the following document: Yamashita Y., Lebert F., Gourry JC, and Bourgeois B.,2013, A practical field experiment of multiple transmission resistivity profiling using code division multiple access. Near Surface Geoscience 2013 - Proceedings of the 19th European Meeting of Environmental and Engineering Geophysics, Expanded Abstracts, Mo S2a 10). In the context of the invention, the chosen tip will consist of three steps ON+ ON- ON+ (or its inverse ON- ON+ ON-), this pattern making it possible to compensate for the PS by a sliding average "on three points" chosen within the same tip.This is a classic method for eliminating a linear drift by averaging the signal at three times q, E, E positioned identically on three consecutive levels, and ultimately obtaining the value of the potential freed from this drift, therefore the potential 1 f. due exclusively to the transmitters in the absence of any other noise source, namely pf - p] -¼ [x - 2 x E, ) / 4, a value which is reported at the center of the triplet (i.e. point E). This extraction carried out before decoding is theoretically valid only for a linear PS. It provides the potential value VPS at the three times considered, s ~ - If, ~ F, -Fi:t and F^ - F3 ~Vn, as well as the slope a of the extracted linear PS: a ~ (K - ) / (t, - q )• .

[0040] This procedure can however be faulty in several situations, even in mono-Tx (i.e. with a single transmitter), in particular:

[0041] - if the PS is strongly non-linear (quadratic, sinusoidal, etc.);

[0042] - or even in the presence of strong variations in the current injected by one or more transmitters. In this case, the potentials VI to V3 above may be affected differently by these variations, leading to an incorrect determination of VTX. These variations should therefore be corrected before averaging over 3 points;

[0043] - or finally in the presence of high frequency noise of high amplitude: in this case also, if the previous correction is done blindly, noise can bias the determination of VTX and VpS.

[0044] The problem is particularly acute for the measurement of induced polarization since it involves measuring a much lower potential than in resistivity (because measured after the injection has been cut off), on which the errors will be reflected much more strongly than on the apparent resistivity.

[0045] The new approach proposed in the invention to estimate the PS (spline adjustment followed by spectral interpolation), makes it possible to subtract a maximum of PS before applying the average over 3 points. Similarly, averaging or "stacking" carried out in several stages allows a good part of the high frequency noise to be eliminated before the 3-point average. a. In the induced polarization method, the minimal tip consists of four steps ON+, 0, ON-, 0 (i.e. square wave with 50% of steps at 0), which is the waveform classically used for PP due to the emission dead times, during which the temporal decay of the potential can be analyzed, leading to the apparent chargeability. These 4 steps represent a period T0 of the emitted signal. This minimal tip is used in the following two documents: Yamashita Y., Lebert F., Gourry J.C., Bourgeois B. and Texier B., 2014. A method to calculate chargeability on multiple-transmission resistivity profile using code-division multiple-access. SEG Technical Program Expanded Abstracts 2014: pp. 1775-1779 et Yamashita Y. and Lebert F., 2016. Time Domain IP Profile by Multi-current Transmission Technique - Water Tank Experiment and SP Noise Effect Estimation. Near Surface Geoscience 2016 - Proceedings ofthe 22nd European Meeting of Environmental and Engineering Geophysics, Barcelona 4-8 Sept 2016, Expanded Abstracts. But, as introduced previously, the tip chosen in practice will rather be a tip with 6 steps ON+, 0, ON-, 0, ON+, 0 to be able to compensate the PS on 3 points. As for the resistivity, the choice of the pattern of one and a half periods is justified by the fact of having two positive slots and one negative slot (or vice versa) on each tip to be able to carry out, at the end of the dead times, a PS suppression by sliding average on three consecutive OFF steps within the same tip (c / [Fig.3]).

[0046] In the following, the reasoning is based on a six-step PP tip allowing PS compensation. Unless otherwise indicated, we will assume that each step has a duration of 2 seconds: in this case, the period T0 is 8 s (therefore the frequency fi is 0.125 Hz) and the duration of the elementary tip is 12 s (1.5xT0).

[0047] [Fig.2] illustrates on the left the three CDMA 26 codes used for the coding of a multi-Tx transmission with three transmitters Txl (reference 18), Tx2 (reference 20) and Tx3 (reference 22). Each CDMA 26 code extends over a duration of 4 bits (i.e. 4 steps worth ±1). These three codes are indeed orthogonal two by two. On the right, this same figure shows the coded signal emitted by each transmitter in the case of a PP measurement. These coded signals are the result of the multiplication of each of the 3 previous codes by the elementary PP tip with 6 steps (circled on [Fig.2] for the signal from transmitter Tx2). We thus obtain for each transmitter a succession of 4 tips which form the “CDMA pattern” specific to this transmitter. The duration of 4 tips of 6 steps, or 24 steps, is the minimum time required to decode the PP signal in the case of three simultaneous transmitters. For elementary steps of 2 seconds, the minimum CDMA pattern of 4 tips therefore lasts 48 s, which corresponds to 4800 samples with a conventional sampling at 100 s / s. It will be noted in the illustrated example that the 3rd code at the bottom of [Fig.2] is a "transparent" code (different from a uniform code), that is to say that the resulting CDMA pattern is identical to the uncoded PP signal used in mono-Tx. The coded transmission sequences 28 of the different transmitters 18, 20 and 22 are represented on the right part of [Fig.2].Below these we find the coded reception sequence 30; the latter corresponding to the sum of the coded transmission sequences 28 of the different transmitters 18, 20 and 22, that is to say to the superposition of the three coded transmitters injecting the same current on the same resistance (for example if the three transmitters considered are on the same ground and if the sequence is measured in an ideal receiver giving the same geometric coefficient for the three Tx), in the absence of PS, PP and MT noise.

[0048] In practice, the injection duration is always longer than a single 4-tip CDMA pattern: in fact, the spline processing used to estimate the PS requires at least one spare tip before / after the current pattern to support the spline: thus, the minimum sequence is in practice 6 tips, or 36 steps. In addition, the CDMA pattern is repeated cyclically for as long as necessary to obtain a good signal-to-noise ratio (with an overall periodicity equal to 6xT0).

[0049] The following remarks are made: 1. Note that CDMA patterns as we define them are preserved by circular permutation, that is to say that by shifting one bit on the code (i.e. one tip on the pattern) we remain on the same code (i.e. the same pattern). Consequently, when decoding CDMA on a pattern of 4 tips, it is not necessary to shift 4 tips to process the following data. It is possible to shift only one tip: this is called "overlapping" processing (the overlap range in this case is 3 tips), a procedure which generally proves to be more robust than non-overlapping processing. 2. It should also be noted that the coding of the transmitters does not modify either the temporal position of the OFF steps (zero steps), nor their periodicity which remains 1 out of 2 (duty cycle equal to 50%). Indeed, the multiplication of the elementary tip of PP by a code consisting of +1 and -1 reverses the ON steps, but leaves the OFF steps unchanged. This is visible in [Fig.2], where we see that the OFF steps always remain in the same places and are well synchronous on the three CDMA patterns of the coded transmitters. Consequently, the combination of the three transmission patterns at the receiver will also leave the temporal position of the OFF steps unchanged, and we will always observe an OFF level framed by two ON levels.

[0050] The transmitters used are transmitters configured to inject the CDMA sequences described above. It is possible that these transmitters are only voltage regulated and therefore do not provide a constant current due to variations in the tap resistance of the AB electrodes over time (frequently up to 100%). This is also the case for certain transmitters used during field validations. For the same reason, each transmitter can give its own current value, depending mainly on the tap resistance of the AB electrodes connected to it.

[0051] On reception, the acquisition of the signals can be carried out for example by two-channel electrical receivers having a sampling frequency of 100 s / s and a useful spectrum limited to the 0-10 Hz band.

[0052] As explained above, the example described to illustrate the invention is limited to the case of three transmitters Txl, Tx2 and Tx3. In the case where only one of them emits a signal, the transmitter Txl for example, the latter injects a current Zi(t) into the subsoil over the time t. Under the usual conditions of electrical prospecting, in which the frequency is low enough to neglect any effect of induction or electromagnetic propagation, the electric potential V(t) created by the current / | (t) at a receiver Rx is given by Ohm's law, in the absence of PS and noise:

[0053] [Math.l] (1)

[0054] where Ri denotes the overall resistance of the subsoil between the transmitter Txl and the receiver Rx, itself given by:

[0055] [Math.2] (2)

[0056] where pai and Ki respectively denote the apparent resistivity of the subsoil and the geometric coefficient between the transmitter Txl and the receiver Rx. Formula (1) is purely real, which means that there is no phase shift between the injected current and the received potential. It should be noted that the electrical parameters of the subsoil (such as pal and R i) are assumed to be constant, i.e. free from temporal variation during the duration of the survey (even if very locally at the level of the transmission, the tap resistances can vary over time).

[0057] If the preceding time signals V(t) and / |(t) are discretized into n samples two by two synchronous, formula (1) can still be written in matrix form:

[0058] [Math.3] £ | i 1 X * ! • 4 =^| 4 jd Uu where we note that Ri does not depend on time. (3)

[0059] In the context of electrical prospecting, the unknowns are the resistances Rj ( / = 1, ni) associated with each Tx which ultimately allow the apparent resistivities to be calculated, knowing the geometric coefficient K, of each Tx-Rx pair.

[0060] If two transmitters Txl and Tx2 simultaneously inject a signal into the subsoil, by virtue of the principle of superposition of sources, the potential measured at the receiver Rx will be:

[0061] [Math.4] (4)

[0062] Finally, in the case of simultaneous transmission by the three transmitters Txl, Tx2 and Tx3, the potential measured at the receiver Rx will be:

[0063] [Math.5] (5)

[0064] It is possible to define the matrix of injected currents I, with n lines corresponding to the n time samples and m columns corresponding to the m simultaneous transmitters (m being equal to 3 in this example):

[0065] [Math.6] Txl Tx2 Tx3

[0066]

[0067] (6) It is possible to note V the vector of the electric potential with n elements observed at the receiver Rx: [Math.7]

[0068]

[0069] (7), and R the vector of the resistances of the subsoil R / j = l,m) between the receiver Rx and each transmitter Tx: [Math. 8]

[0070]

[0071]

[0072] (8) Formula (5) above can therefore be written in matrix form as follows: [Math.9] V - IR ( !Twi ) (9) Formula (9) corresponds to Ohm's law extended to several sources and allows to express that the potential V, measured at a time sample i is a linear combination of the currents injected by the three transmitters at the same time sample (without taking into account noise and spontaneous polarization), that is to say:

[0073] [Math. 10] m Vze[U], r = ^44=44+44+44 7=1 (10)

[0074] At this stage of the processing, the number of samples n represents all of the data selected for CDMA decoding (i.e. a subset of the time sequence common to the 4 recorders) and can therefore be at least 4800 for 2 s steps sampled at 100 s / s (i.e. length of a CDMA pattern of 4 tips, the 2 support tips of the splines not being counted in the number ri).

[0075] The invention aims firstly to solve the system of linear equations corresponding to formulas (9) or (10) to determine the vector of resistances R of the m Tx-Rx couples, and this even in the presence of significant disturbances in the received signals V (affected by PS and noise) and emitted I (affected by low frequency drifts and injection instabilities if the Tx are not current regulated).

[0076] In the presence of spontaneous polarization, it is imperative to eliminate the latter as much as possible before determining the resistivity and the induced polarization of the subsoil. As mentioned above, the known method known as the "three-point average" is not effective in the case of a non-linear PS or in the presence of HF noise.

[0077] For the first of these disturbances, a preliminary treatment has been developed aimed at estimating the PS by fitting piecewise polynomial functions (or splines), followed by spectral domain filtering aimed at rejecting the frequency / 0, and this PS estimate is then subtracted from the measured signal.

[0078] As for the HF noise, we eliminate it as much as possible by averaging the time samples selected on each useful range of each step, thus reducing the number of samples per step to 1 (for each tip, we thus obtain 3 samples on the ON steps, used to decode the resistances, and 3 samples on the OFF steps, used to decode the chargeabilities).

[0079] After these averages, the residual PS is finally eliminated by filtering on three points, thus reducing the number of samples to one per tip for the resistance and one per tip for the PP, a number necessary and sufficient for the CDMA decoding of the current pattern.

[0080] To have synchronous data of the same size between the receiver Rx and the three transmitters, a similar treatment is applied to the time series of the currents injected into the subsoil: a) adjustment of splines to eliminate the continuous components and the possible slow drifts of the emission sequences (but without filtering spectral because the cutoff of the Tx is independent of the polarity of the current), b) averaging of the ION signals over each useful time range at the end of each ON level.

[0081] Let us now go into the details of the processing allowing the cleaning of the measured signals before decoding.

[0082] The first step advantageously consists of analyzing the signals from the transmitters Tx and the receiver(s) Rx (step 4 of [Fig.l]) to find the position of the rise ramp at the start of the chosen tip (time t0). This involves, for each time series, searching for the order numbers (indices) of the points corresponding to the beginnings and ends of each emission level: this is done by derivation of the signal. In a non-automated implementation, the main difficulty for the operator consists of locating the start of each tip, invariably characterized by a doublet of consecutive slots of the same polarity on at least one of the coded transmitters: it is then necessary to point to the start of the second slot of this doublet (see [Fig.2]). From this t0, the operator chooses the number m0 (integer) of patterns that he wishes to process, with m0>l or at least one pattern of 4 tips.At this stage, the number n of samples in formulas 1 to 10 will therefore be worth 48OOxmo for 2 s steps sampled at 100 s / s. Note that the time t0 and the number of patterns selected must allow for a complete tip to be available to the left and right of the time sequence to support the splines.

[0083] Following this first step, it is advantageous to carry out an estimation and a suppression of the low frequency component of the transmission signals (step 6): the latter sometimes have a continuous component, or even a slow drift which is somewhat analogous to a PS. This “transmission drift” is likely to distort the CDMA decoding. It is therefore advantageous to estimate this component, for example using cubic splines based on all the zero steps of the transmission signals. This component is sought for each transmitter j over the entire selected transmission duration, then it is subtracted from the signal of each transmitter to obtain a cleaned transmission signal, stored in the matrix S_clean: the latter has the same structure as the matrix I seen above (n rows and 3 columns, where n can be 4800xmo).

[0084] Then, it is advantageous to carry out an estimation of the PS affecting the reception signal (step 8). In a first phase, the PS is estimated over the entire reception duration selected by fitting a cubic spline based on the end of the steps to zero. To reduce the influence of ambient HF noise, the potentials of the interpolation support points are obtained by performing the arithmetic mean of the potentials of the last ten points of the OFF steps. The result of this processing is a time series denoted PS_est, of duration equal to that of the selected time window, is at least a pattern of 4 tips. The vector PS_est has the same structure as the initial vector V (in particular the same number of samples n) and is synchronous with it.

[0085] Following this estimation of the PS, a spectral domain filtering of the time series described by the vector PS_est is carried out (step 10). Indeed, the PS estimated in the previous step is affected by a frequency / 0 ripple, i.e. the frequency of the emission signal, which is explained by the fact that the support points of the splines are shifted alternately upwards or downwards relative to the real value of the spontaneous polarization (shift due to the induced polarization residue).

[0086] To remove these ripples, we move into the spectral domain by Fourier transform and we perform a suppression of the spectral components ranging from / 0 / 2 to 3 / 0 / 2 (limits included). We then reconstruct these components by “pchip” interpolation (for “Piecewise Cubic Hermite Interpolation” or piecewise cubic Hermite interpolation) based on the frequencies located on either side of the suppressed spectral band. The interpolation is carried out, on the abscissa on the logarithm of the frequencies and, on the ordinate on the logarithm of the spectrum expressed in polar coordinates, that is to say on the decimal logarithm of the amplitude and on the spectral phase of the estimated PS.We then return to the time domain by inverse Fourier transform to obtain the filtered spontaneous polarization estimate, which is a new time series noted PS_est* having the same structure as, and synchronous with, the sequence of potentials currently being processed V (therefore always n time samples where n is for example 4800xmo).

[0087] [Fig.3] illustrates an estimation of the PS by cubic splines (curve 32) followed by spectral filtering (curve 34), compared to a classic estimation of the PS by the average over 3 points (curve 36) on synthetic resistivity-PP data in the case of emission steps of 4 seconds. This figure illustrates the great importance of the spectral filtering according to the invention, which makes it possible to pass from curve 32 to curve 34, (which is almost perfectly superimposed on the true PS, barely visible in dashed line under curve 34), but also the significantly lower result of the classic average over 3 points (curve 36) thus clearly showing the benefit of the invention on a textbook case.

[0088] Finally, after estimation and filtering of the PS in the spectral domain, the PS vector obtained PS_est* is subtracted from the potential vector V currently being processed (step 12). This gives a new potential vector cleaned of the PS, denoted P_clean, which will then be used with S_clean for decoding the resistance and the chargeability.

[0089] At this stage P_clean and S_clean comprise a large number of lines (eg 4800xmo) corresponding to the n samples of the time window selected for the processing (subset of the sequence common to the 4 recorders, including at least one CDMA pattern of 4 tips).

[0090] At the end of the previous cleaning steps, and having made the experimental observation that the CDMA decoding propagated noise and errors in an unstable manner, we then proceed, before decoding, to smoothing the data (step 13 of [Fig.l]) in two successive reduction steps:

[0091] a) arithmetic averages (“stacking”), aimed at eliminating HF noise as much as possible, averages carried out on the useful ranges at the end of each level, ranges which differ according to the type of data to be decoded, resistances or chargeabilities. For the decoding of the resistances, the ranges in question are formed of the ten final samples of each ON level. For the decoding of the twenty chargeability time channels (see below), the stacking ranges are smaller and comprise only four samples per decay window.

[0092] b) three-point filtering carried out on three analogous average values from the previous step, positioned identically on three levels of the same type within the same tip, aiming to eliminate PS residues at reception and BF drift residues at transmission. The value obtained is reported at the center of the tip.

[0093] These two steps drastically reduce the number of samples: in the end, only one value of each species is retained per tip, namely, for each tip: i) a value of Ion per Tx and a value of V0N at the Rx, measuring the average height of the tip emission slots (ON steps), and ii) a value of VOff(L) at the Rx for each time decay window defined on the OFF steps (the time channel index k can vary for example from 1 to 20).

[0094] For each physical data to be decoded (resistances or chargeabilities) we thus obtain four values for the CDMA pattern being decoded: we thus move from a very large number of time samples (n=4800 as standard) to a very reduced number (A=4) of reliable data. These four values, which represent a summary of the information of the 4-tip pattern considered, are in the necessary and sufficient number to carry out the CDMA decoding of a multi-Tx acquisition with 3 transmitters (for comparison, in mono-Tx it would be sufficient for each Tx to have an analysis duration 4 times shorter to obtain an estimate of the same parameters, but to the detriment of the S / N ratio since the data would be less averaged: we can therefore retain that for the same total acquisition duration with 3 transmitters, the multi-Tx implementation brings added value in terms of S / N ratio due to the necessary extension of the analysis duration).

[0095] The above processing operations can be formalized using the following formulas (11) and (12), which relate to a CDMA pattern being decoded.

[0096] First, from P_clean, it is possible to define a new vector of po potentials, noted P, containing the four average potential values at the end of the ON stage (one per tip), defined as follows for i=l to N (N being equal to 4 in the present case for an emission with three transmitters, but remember that it would be equal to 8 in the case of 7 simultaneous transmitters).

[0097] [Math. 11] ViefLA], = - [^,,(1)+ 1+O)-2yJ^,^ (11)

[0098] where the horizontal bar designates the arithmetic mean over the 10 useful points at the end of the ON level and where the symbol "..." placed above it designates the average over three points carried out on these averages. The figures in parentheses designate the numbers of the unaffected levels on each tip.

[0099] Similarly, from S_clean it is possible to define a new current matrix noted S, always concerning the pattern being decoded, for the m transmitters:

[0100] [Math. 12] (12)

[0101] For yV=4 and m=3, the vector of useful potentials P and the matrix of filtered currents S can be written, for the pattern being processed:

[0102] [Math. 13] Rx Txl % S2, Tx2 Tx3 Ç\s,. P. P- s- and S ~ R s32 U 7 S y ï -SC, (13)

[0103] If the window selected for decoding contains several patterns (m0>l), the processing will be repeated by shifting by one tip (for so-called "overlapping" processing) until reaching the last pattern of the window. The data obtained at each iteration will then be averaged after CDMA decoding, as will be seen later.

[0104] Following the suppression of spontaneous polarization and ambient noise, the vector of useful potentials P and the matrix of filtered currents S are linked by the following matrix Ohm's law:

[0105]

[0106]

[0107]

[0108]

[0109]

[0110] [YES] [Math. 14] P = SR (14) To solve equation (14) above (step 14), it is known to multiply both sides of the equality by the transposed matrix of S. We then obtain the following equation: [Math. 15] S T P - S T SR (15) In the case of perfect transmission data where the columns of S are strictly orthogonal codes, the decoding stops there, the STS product being diagonal and reducing to the square norms of each CDMA code. The determination of the three elements of the vector of resistances R is then obvious knowing that the left term STP of equation (15) is already a 3-element vector (when m=3). This is what is practiced in the articles of Yamashita et al. (2013, 2014, 2016). For the case of imperfect transmission data (especially in the case of non-current-regulated transmitters), the STS matrix is not diagonal. However, it is a square matrix (and symmetric and diagonal dominant), so it is possible to invert it and multiply the two terms of the equality by this inverse. By rearranging the previous equation, the CDMA decoding of the resistors for the common CDMA pattern then amounts to solving: [Math. 16] R = (S T S) 1 S'P (16) $i is a matrix of size mxN, known as a pseudo inverse of the matrix S. Its physical dimension is 1 / A, so equation (16) consists of dividing potentials by currents (i.e. inverted Ohm's law). Since the R vector has been obtained for a pattern of 4 tips, we can shift the sequence by one tip and repeat the same operation on the next pattern of 4 tips (so-called "overlapping" processing). We proceed in this way until we reach the end of the chosen injection sequence, then we average the values obtained on the different sets of 4 tips processed. The final vector of resistances R thus calculated will finally allow, by the formula term

[0112]

[0113] (2), to obtain the apparent resistivities of the subsoil.

[0114] Once the decoding of the resistors has been carried out and the vector R obtained, it is possible to find the vector of the theoretical coded potentials V / received at the receiver from each transmitter j (i.e. as if the transmitter j emitted alone) using the scalar Ohm's law. This is possible since we know, for each sample i of the time series, the current ly injected by each transmitter j:

[0115] [Math. 17] (17)

[0116] Similarly, on the reduced data (i.e. after averaging and 3-point filtering), for the N tips of a given CDMA pattern, we can define the vector of average theoretical coded potentials W, coming from transmitter j at the end of the ON steps:

[0117] [Math. 18] Vz e V / e [1, m], (18) (N being 4 in the present case where m=3).

[0118] Knowledge of these theoretical potentials is important for decoding chargeability. Indeed, the partial chargeabilities M(t) and global M are generally defined by the following formulas, given for a single emitter:

[0119] [Math. 19] (19)

[0120] [Math.20] a / --......-.....vy.w)A / = ( L -- / , ).. jN (20)

[0121] where Von is the (average) end-of-charge potential and VOff(0) the discharge potential observed at a time t after the current has been cut off (formulas to be multiplied by 1000 to obtain the results in mV / V). The overall chargeability M is none other than the average of the partial chargeabilities over the overall time decay window [tl,t2] defined by the operator.

[0122] The different time chargeability windows defined for the range of electrical devices used are used. The positions h of the window centers and their widths depend on the window definition mode selected (arithmetic, logarithmic) and the chosen step duration, itself adapted to the relaxation times of the polarization phenomenon studied.

[0123] We can now approach the CDMA decoding of chargeabilities (step 16). The procedure for decoding chargeabilities essentially follows that of resistances. The only conceptual difference is that decoding must now be done in the absence of transmission current: it will therefore not be able to involve synchronous signals as previously where the currents and potentials were pointed at the same time on the same transmission level. To allow CDMA decoding of potentials measured after the transmitters have been switched off, it is known to extrapolate the current injected by each Tx before the switch-off onto the dead time (Yamashita et al., 2014, 2016). Decoding chargeabilities Mj^) is then possible with the use of these fictitious currents.

[0124] Formula (19) suggests that the determination of the chargeability M(L) is analogous to the determination of the resistance R from formula (1), i.e. R ~ 1 in which the potential VOn is replaced in the numerator by the potential Voff(L), and in the denominator the current Z0N injected at the end of charging by the potential RxJqn produced by this current.

[0125] In CDMA mode, a possible approach to obtain the chargeability M / L) induced by each transmitter j is therefore to replace in formula (16) the matrix of currents S by the matrix W of the theoretical potentials V0N of the various Tx before cutoff, matrix constructed as follows from formula (18) for the particular case where m =3:

[0126] [Math.21] (21)

[0127] where Wy represents the average theoretical potential at the end of the ON steps created by the transmitter j during the tip i of the current pattern (average carried out on the same useful range of 10 samples used to determine the Rj in the previous step, then expurgated of the PS residue by the 3-point method).

[0128] In this approach, the decoding of the chargeability then amounts to solving the following equation, for each CDMA pattern processed and for each time window k:

[0129] [Math.22] M(It) = (wTwy' WT Vo„«,) (A4 (A xi) (22) where the N values of the vector VOff(C) are defined by the formula (23) below, for the N tips of the current pattern (A=4 in our example):

[0130] [Math.23] Vi e [I, V], V^) = .1 ) + .3)- 2 x 447(^.2)] / 4 (23)

[0131] This formula is analogous to formula (11) seen above for the ON steps: for each time decay window k, on each of the three OFF steps existing inside tip i of the current pattern, we calculate the arithmetic mean of the elements of P_clean belonging to this window, then we carry out the three-point filtering of these three means. We thus obtain, for each time channel k, an OFF potential value per tip.

[0132] It will be noted that the term s | ' WT 'formulated (22) is the pseudo-inverse of W, of size mxN and physical dimension 1 / Volt: equation (22) therefore consists of dividing potentials by potentials.

[0133] Alternatively, thanks to the prior determination of the S y (which correspond to the I0N of formula 19) in step 13 and the subsequent calculation of their pseudo-inverse (S*S) ! S1 in step 14, it is possible to calculate directly from VOff(O a new vector having the physical dimension of a resistance, noted R*(h) and defined by:

[0134] [Math.24] (wri) (w A4 < V. 51 (24)

[0135] The elements Rj*(h) of this vector (j=l,m) are called pseudo-resistances associated with each emitter j, for each window k of the time decay.

[0136] If we reduce the problem to a single transmitter, formula (24) can still be written: then, by manipulating formula (19), we see that: - M(i)xR

[0137] therefore, in the case of a single transmitter, we obtain:

[0138] [Math.25] = (25)

[0139] Extended to a distributed transmission from several simultaneous transmitters according to the invention, this formula becomes, combined with (24):

[0140] [Math.26] W*)OR = R'(y - (s'sf S T V OFF ( / A ) (26)

[0141] where ® denotes the Hadamard product (element-by-element matrix product).

[0142] The method for decoding chargeabilities according to the invention is therefore as follows: once the vector of pseudo-resistances R*(fû is obtained by formula (24), we divide the result by the vector of resistances R (element by element division) to deduce the vector of partial chargeabilities M(t,0 for each window k, i.e. in the case where m=3:

[0143] [Math.27] (27)

[0144] The overall chargeabilities Mj associated with each transmitter are then calculated from the M / L) using formula (20).

[0145] As for the resistors, the previous processing is repeated for each processed CDMA pattern from the beginning to the end of the chosen time sequence, then the obtained values are averaged to provide the final estimates of the partial and global apparent chargeabilities associated with each of the m Tx-Rx pairs. These estimates are stored in the vectors M(L) and M provided to the user.

[0146] Thanks to the CDMA decoding according to the invention, we ultimately obtain the apparent resistivities and the partial and global apparent chargeabilities of the subsoil at the level of

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157]

[0158]

[0159]

[0160]

[0161]

[0162]

[0163]

[0164]

[0165] each transmitter-receiver pair. This decoding is made reliable by the prior estimation and suppression of spontaneous polarization and ambient noise. The geoelectric parameters obtained subsequently make it possible to make hypotheses on the nature of the subsoil in the prospected area. List of references 2: Measurement data recovery step 4: Analysis step 4 of Tx and Rx signals 6: step of estimation and removal of the continuous component 8: Spontaneous polarization estimation step 10: spectral domain filtering step of the estimated spontaneous polarization 12: Subtraction step of the spontaneous polarization estimate 13: Data reduction / smoothing step 14: CDMA decoding step of the resistance vector 16: CDMA decoding step of the chargeability vector 18: Txl transmitter 20: Tx2 transmitter 22: Tx3 transmitter 24: tips 26: CDMA codes 28: coded emission sequences 30: coded reception sequence with three transmitters 32: estimation of spontaneous polarization by cubic splines 34: spectral filtering 36: classic estimate of spontaneous polarization by the average over 3 points

Claims

1. Claims Method for measuring the induced polarization and resistivity of a subsoil comprising at least the following steps: - selection and reading of measurement data from simultaneous transmitters injecting coded electrical signals into the subsoil and from at least one receiver capturing an induced signal restored by the subsoil following the injections of said electrical signals by all the synchronous transmitters, the data forming a reception sequence from the receiver, the reception sequence being the superposition of coded potentials resulting from coding of the transmitters according to code division multiple access coding and forming different coding patterns configured to discriminate the signal emitted by each transmitter, - estimation of a spontaneous polarization (8) of the subsoil at the receiver level, over the entire reception sequence, - spectral domain filtering of the estimated spontaneous polarization (10), - subtraction of the filtered estimated spontaneous polarization sequence from the reception sequence (12), - decoding, for each coding pattern processed in the reception sequence, of the resistance of the subsoil associated with this pattern for each transmitter (14), - averaging the resistances obtained for all the coding patterns processed in the reception sequence, for each transmitter used, and obtaining the average resistance then the apparent resistivity associated with each transmitter, - decoding, for each coding pattern processed in the reception sequence, the partial chargeabilities and the overall chargeability of the subsoil associated with this pattern for each transmitter used (16), and - averaging the partial and global chargeabilities obtained for all the coding patterns processed in the reception sequence, for each transmitter used, and obtaining the partial apparent chargeabilities and the global chargeability associated with each transmitter.

2. A measuring method according to claim 1, wherein the retrieved measurement data is from non-current regulated current transmitters.

3. A measuring method according to any one of the preceding claims, wherein the induced polarization and the resistivity are measured from several reception sequences from several receivers.

4. Measuring method according to any one of the preceding claims, in which the coded sequences comprise several elementary patterns of the coded injection sequences, i.e. 24 ON / OFF emission levels in the case of three transmitters, and in which the data analysis is carried out on at least one of these coding patterns, and in which an estimation of the spontaneous polarization (8) is carried out.

5. Measuring method according to any one of the preceding claims, in which, before the estimation of the spontaneous polarization (8), an estimation and a suppression of a continuous component and / or of a slow drift of the current injected by each emitter used are carried out.

6. A measuring method according to any one of the preceding claims, wherein, during the step of filtering in the spectral domain of the estimated spontaneous polarization (10), the following steps are carried out: - calculation of a Fourier transform from the estimate of the spontaneous polarization, - suppression of spectral components ranging from / o / 2 to 3 / o / 2, / o being the frequency of the signal from the transmitters resulting in a ripple frequency affecting the estimated spontaneous polarization, and - reconstruction of the suppressed spectral components by piecewise Hermite cubic interpolation, and - obtaining a filtered spontaneous polarization by calculation of an inverse Fourier transform.

7. Measuring method according to any one of the preceding claims, in which, during the step of decoding the resistances (14) and the partial and global chargeabilities (16), a current matrix in which the DC component has been subtracted, and a potential vector in which the filtered estimated spontaneous polarization has been subtracted.

8. Measuring method according to any one of the preceding claims, in which a vector of chargeabilities is determined from a vector of pseudo-resistances observed during the time decrease.

9. Method for determining the nature of a subsoil, comprising the following steps: - injection of coded electrical signals into the subsoil via several simultaneous electric current transmitters, each transmitter injecting into the subsoil an electrical signal comprising a code specific to the transmitter, - capture, by at least one receiver distinct from the transmitters, of an induced signal restored by the subsoil following the injection of electrical signals, the received signal forming a reception sequence, and - implementation of a method for measuring the induced polarization and the resistivity of the subsoil according to any one of the preceding claims.

Images