Evaluation of induced polarization measurements using a new inversion method

In this paper, a new inversion method is proposed to process laboratory-measured induced polarization (IP) data. In the new procedure, the concept of the series expansion-based inversion is combined with a more general definition of the objective function. The time constant spectrum of the IP effect is assumed a line spectrum approximated by a series of Dirac’s delta function resulting in a square-integrable forward problem formula. This gives the applicability of the generalized objective function. The expansion coefficients as unknowns represent the model parameters of the inversion procedure. We use the new inversion procedure on an apparent polarizability dataset measured on a rock sample originated from the Recsk ore complex, northeast Hungary. The inversion results was compared to those of three additional laboratory datasets, which were measured on samples rich in ore minerals collected from the same area. The results are compared to those given by the traditional series expansion-based least squares method. It is shown that the newly proposed method gives more accurate and stable parameter estimation.


Introduction
The demand for mineral raw materials is constantly increasing due to the expanding population of the world and in this context the needs of the population and industry. Replenishment of mineral resources (energy carriers, ores, building materials), drinking water and the reduction of the risk posed by geological hazards require more and more accurate and detailed 1 3 knowledge and mapping of the underground by earth scientists. One of the most suitable tools for this is geophysical measurements, including geoelectric methods. Related to direct current (DC) test methodology, the measurement of induced polarization (IP) is a very efficient method. Previously, it was used primarily in connection with ore exploration (Seigel 1959;Wait 1959;Keller and Frischknecht 1966;Tavakoli et al. 2016) due to the significant effect of metallic polarization because of metallic minerals. In addition, it plays a role in the study and characterization of environmental pollution (Viezzoli et al. 2006;Turai 1985Turai , 2011, exploration of geological structures, coal and graphite, separation of sandy-clay rocks and archaeological geophysics (Abu Zeid et al. 2016).
The IP phenomenon appears as an exponentially decaying effect in time of polarization currents induced in subsurface polarizable rock formers or human-derived structures (buried wastes). The form of these mathematical relationships in the field of geophysics and physics is similar to the mathematical description of many other phenomena, so the inversion methods of IP data evaluation include more general applicability in other fields.
The application of the IP method to direct hydrocarbon research is an untapped opportunity (Turai and Dobróka 2002). The basis for this possibility is that oil-eating bacteria release hydrogen sulphide if the geological trap contains oil or natural gas. At the edges of the hydrocarbon trap, the hydrogen sulphide is released and pyrites the iron oxide content of the cover rocks from the depth of the reservoir to the surface, thus in the edge zone of the hydrocarbon reservoir so-called pyrite chimneys are created. As pyrite is a highly polarizable ore mineral, ring-like IP anomalies can be measured over productive reservoirs (containing oil and/ or natural gas) by surface measurements. If no IP anomaly can be measured at the edges of a reservoir, then that structure does not contain hydrocarbons, only thermal water. The inexpensive way of direct hydrocarbon exploration described above could be of great importance in locating new hydrocarbon reservoirs in the future, and thus in satisfying the growing energy demand.
In this paper, a new spectral approach of inversion processing of IP measurements is presented. The main point of this is that the metric measuring the difference of the direct problem describing the phenomenon and measurement results are defined as a distance function valid in the space of square-integrable functions. Since the model function during the inversion can be integrated into a closed form, with this choice the numerical uncertainties of the inversion can be reduced.

Induced polarization measurements
The phenomenon of induced polarization observed in DC measurements was first described by Schlumberger (1930), then Dakhnov (1941) and Bleil (1953), and the first mathematical model was given by Seigel (1959). The essence of measuring IP in the time domain is that when the excitation current is switched on, the polarizable structures are first electrically saturated. After the voltage has stabilized during the excitation, the excitation current is switched off and the monotonically decreasing voltage between the measuring electrodes is recorded at given sampling intervals. The quotient of the maximum voltage during secondary (measured during decay-U(t) ) and primary excitation (U 0 ) gives the value of the apparent polarizability

3
A frequently used quantity is also the apparent chargeability (m a ) which is proportional to the area under the decay curve. Its value can be specified in the interval t 2 − t 1 as The following effects can be identified in the background of the polarizability of rocks. Membrane polarization occurs in ionically conductive porous rocks due to the negative surface charge of clay particles. During electrode polarization, an electrochemical interaction occurs between the electronically conductive particles in the rock and the ionic solution in the pore space. Filtration polarization occurs in sedimentary assemblages bound to the conductive fluids filling the pores due to the difference in ion concentration because of the different mobility of negative and positive ions. Redox polarization occurs in soils and rocks containing oxidative or reductive chemical components, for example, due to chemical contamination.
Thus, excitation can result in several polarization mechanisms in rocks, even simultaneously which can be characterized by different characteristic time constants. The combined effect of these appears in measurement data, and if they can be separated, they can provide basic information about the structure, lithology and inhomogeneities of the investigated area. In a real physical environment, the separate polarization processes can be countable, but the phenomenon can also be understood as a polarization system characterized by an infinite number of time constants (Turai 1981). One will see that in the inversion evaluation of IP series expansion based inversion will be used, however, as a result of the above, this approximation can also be considered as a proper description of the IP problem. In this paper, the IP decay curves recorded in the time domain will be evaluated by a new spectral inversion method called the G_LSQ algorithm and compared with the previously used traditional method (T_LSQ method). Our goal is to present an inversion methodology that increases the accuracy and reliability of series expansion based inversion solutions if the basis functions can be integrated into a closed form.

Data processing by series expansion based inversion method
The Department of Geophysics of the University of Miskolc has been successfully developing geophysical inversion methods for several decades, during which a new research direction the series expansion based inversion has played an increasing role. The essence of this is that the processing and interpretation of data measured on complex (laterally and vertically inhomogeneous) geological structures is performed by series expansion based discretization and the inversion procedure formulated for expansion coefficients. The greatest advantage of the method is that even with the introduction of relatively few (several times ten) expansion coefficients, a suitable resolution can be achieved so that the problem to be solved leads to an overdetermined inverse problem. The method has been used in several fields: gravitational Völgyesi 2008, 2010), DC geoelectric (Gyulai et al. 2010(Gyulai et al. , 2017, magnetotelluric (Dobróka et al. 2013), borehole geophysics (Dobróka and Szabó 2010;Dobróka et al. 2016), data processing (Vass and Dobróka 2010;Dobróka et al. 2015), induced polarization (Turai et al. 2010;Turai and Dobróka 2011). In this paper the method of series expansion based inversion is used to evaluate induced polarization data measured on ore-bearing rock cores in laboratory, to determine the time constant spectrum and to isolate the individual polarization effects.
The direct problem was formulated by Turai (1981) with the following integral equation where t is the time elapsed since the excitation current was switched off, is the time constant and w( ) is the time constant spectrum. Our goal is to determine the time constant spectrum for which Turai (1981) introduced the TAU transformation operation which can be performed by several methods, e.g., Turai (1985) developed polynomial interpolation and Fourier series solutions. In the following the TAU transformation as an inverse problem will be solved.
Since the time constant spectrum w( ) is an unknown function, discretization is required in the first step, i.e., the spectrum must be defined by a finite number of parameters. The discretization can be performed in the form of a series expansion according to a properly chosen basis function system Φ q where expansion coefficients B q represent the unknowns of the problem to be determined and Φ q ( ) denotes the known basis functions, Q is the number of series expansion parts. Substituting the series expansion notation of the time constant spectrum formulated in Eq. (5) into the answer equation Eq. (3) representing the direct problem, one can get the following relation which gives the (theoretical) value of the apparent polarizability at the moment t k after the excitation current is switched off. Introduce the notation G kq for the integral (Jacobian matrix) defined in Eq. (6), i.e.
It can be seen that by introducing Eq. (7), Eq. (3) takes a very simple form, the calculated values of the apparent polarizability are obtained as a linear expression of the expansion coefficients The appropriate choice of basis functions depends on the task, considering also any prior knowledge (assumptions). Series expansion according to Dirac delta functions fits the previous practice of the IP method, which is suitable for describing the so-called "line" time constant spectrum. Then Eq. (5) can be written as where q is the time constant belongs to the q-th discrete polarization mechanism. Equation (6) formulates the direct problemwhich due to the integration property of Dirac delta can be written as from which one can see that the meaning of expansion coefficients is the spectral amplitude belongs to the given time constant q . Then the elements of the Jacobian matrix can be produced as where q = −1 q . The vector containing the difference between the measured and calculated data can be written as The objective function of the inverse problem can be formulated as a suitably chosen norm of this vector.

Traditional T_LSQ algorithm
The following choice of the objective function is common in the inversion of discrete data (least-squares method) where N is the number of measured data. The solution leads to the known system of normal equations The solution of Eq. (10) to the expansion coefficients is As a solution, the elements of the model parameter vector estimated in the last iteration step are accepted, with which the calculated values of data are closest to the measured data (according to the Gaussian least squares method).
With coordinate writing the matrix of the system of equations and the transformed data vector are with which the normal equation can be written in this way After determining the vector of the model parameters, substituting the expansion coefficients into Eq. (5), the time constant spectrum can be calculated. By plotting the time constant spectrum as a function of the time constants, the individual polarization effects become distinguishable.

Generalized G_ LSQ algorithm
The simplicity of Eq. (8) where Δt is the time elapsed between two consecutive measurements (assuming equidistant sampling). Let us introduce the notations and Then Eq. (15) leads to the next result using Eqs. (16) and (17) or in matrix form This inhomogeneous system of linear equations can be solved for series expansion coefficients and thus the spectral coefficients can be determined The TAU transformation as an inverse problem can also be given in general form using the G_LSQ method. Let us introduce a column vector defined by the following elements from which, using Eqs. (19)-(21) one can obtain the expansion coefficients representing the approximation of w( )-which can be interpreted as an approximation of TAU transformation-as As can be seen, the G_LSQ method produces the matrix of the system of equations analytically, therefore it is expected to give a more accurate inversion result compared to the T_LSQ method. To investigate the relationship between the two methods, write the matrix in Eq. (12) using Eq. (9) as where assuming equidistant sampling t k = kΔt . Introducing the notation q = exp(− q + l Δt) the matrix can be produced as a sum of geometric series If the quantity q + l Δt is small, by replacing the denominator with the first member of its Taylor series the following result is obtained which is equal to the matrix given in Eq. (19) if time is measured in Δt units ( Δt = 1 ). If the condition of uniform sampling and Taylor series approximation is met, the two methods give the same result.
However, the G_LSQ procedure can be further improved by refining Eq. (17) if data are approximated linearly between two samples: With partial integration, Eq. (23) takes the following form Then the improved system of equations is where A ql is given by Eq. (19), r l by Eq. (24) and the solution in matrix form is With the produced expansion coefficients, a linear time constant spectrum = −1 can be given.

Application to laboratory-measured dataset
In the late 1970s and early 1980s, intensive method and instrument development took place at the Department of Geophysics (University of Miskolc) for induced polarization laboratory measurements. This resulted in a high-precision measurement system and datasets sampling fully the polarizability curve. At this time, the measured data were processed by the methods of calculating the logarithmic derivative apparent polarizability curve and decomposition into exponential components, as well as by spectral analysis. In this paper, using one of these datasets, we will show that the proposed series expansion based inversion G_LSQ method can be used to process and interpret laboratory measurement data, as well as several IP effects can be isolated on the generated time constant spectra.
The development of an instrument suitable for induced polarization measurements in the laboratory required persistent experimentation and careful preparation. First, it is important to cut the core samples along parallel planes for good electrical coupling. Before measurement, rock samples were placed in distilled water for 24 h. In addition, it is important to reduce the phenomenon of electrode polarization at metal electrodes connected to rock, which would add noise to the measurements. For this purpose, non-polarizing electrodes were used. Namely, the copper electrodes were soaked in a saturated solution of their salt (100% copper sulphate solution rendered appropriately gelatinized with agar) and the solution provided electrical coupling to the sample by dissociated ion exchange through porous plates. The constructed sample holder is shown in Fig. 1, which was placed in a metal box acting as a Faraday cage to isolate it from electrical impacts in the measurement environment. A detailed experimental study has shown that this electrode configuration gives the best signal to noise ratio with excellent sensitivity. The power electrodes A and B indicated in Fig. 1 were connected to a bipolar current source, thus creating three excitation states: I g = 0, I g < 0, I g > 0 (Turai 1979), here I g is the excitation current.
Preliminary studies have shown that it is expedient to use a current of 10 μA for IP measurements, which in this case corresponded to a current density of 1.04 μA/cm 2 . At the beginning of the measurement, the natural potential baseline was determined, i.e., there was an awaiting time until voltage measured in the unexcited state was stabilized in the measuring circuit. To determine the apparent polarizability, after an excitation time of 15 min, the power source was turned off and the potential difference between the measuring electrodes M and N was recorded for nearly 1000 s. Sampling was based on a combined lin-log time series according to the algorithm where t 0,0 = 0.125 s , a = 2 , = 0.1 , i = 0, 1, … , 9 , k = 0, 1, … , 12 . Reference times are given in Table 1. A research report (Turai 1984) summarized IP data measured on 21 rock samples. In our investigations, the sample marked 1 (Sample 1. in the following) was first selected for discussion. The sample is a highly polarizable scarnic iron ore from the Recsk basin.
A graphical representation of the measured polarizability values is shown in Fig. 2. It can be seen that Sample 1 was well polarized due to its iron content. Its measured polarizability value decreased from 92.7 to 1.6% in 972.8 s.
During data evaluation, the time constant domain was divided into Q = 100 equal parts, and the spectral amplitude B q was considered to belong to the time constant value defined by the middle of the q-th interval (q = 1, …, 100). The number of measurement data was 130, i.e. the inverse problem was overdetermined. In order to ensure the positivity of the spectral amplitudes, the parameters b q = log(B q ) were introduced as new inversion variables. In this nonlinear inversion, the number of iterations was uniformly chosen to be 40 in each case, as further iterations no longer resulted in a significant change in the value of the estimated model parameters. The quality of the inversion (26) t k,i = t 0,0 a k (1 + i), Fig. 1 Schematic picture of the instrument used for laboratory IP measurements result was characterized by the relative distance in data space, the estimation errors of parameters and the correlation norm. The relative data distance is given by the following formula Table 1 Reference times for the IP measurements in seconds (i = 0, 1, … , 9, k = 0, 1, … , 12) The covariance matrix in parameter space for the variables b q = lg(B q ) can be calculated as (Menke 1984), where − is the generalized inverse matrix (built from the derivatives of polarizability data with respect to the b q elements), while ( ) is the covariance matrix in data space. It was supposed that The estimation error of the expansion coefficients is calculated as the square root of the relevant element in the main diagonal of . The correlation matrix of expansion coefficients is given byfrom which the correlation norm (mean spread) can be calculated as where ij is the Kronecker delta.

Inversion results
The inversion of the measured data was performed first with the traditional T_LSQ algorithm. As a starting model to all (logarithmic) inversion variables were given a value of 0.15. Then according to Eq. (27), the distance of measured and calculated (on the starting model) data is D = 21.14, which is a very large value. After 40 iterations the distance is D = 0.0377. This good fit is demonstrated in Fig. 3.
The line spectrum obtained as a result of the inversion method is shown in Fig. 4. In our experience, there is no excitation mechanism in the time constant domain above 500 s. Since Q = 100 unknowns were chosen, the possible time constant intervals were equally defined as 5 s long. The spectral lines (expansion coefficients) shown in Fig. 4 appear in the middle of these intervals. The associated values (time constant, spectral density, estimation error and relative error) for each spectral line (greater than 0.001) can be seen in Table 2. The sample shows remarkable excitation at time constants of 5, 45 and 350 s, smaller amplitude occurs at 110 and 115 s. The average correlation of the inversion variables is indicated by the value of the mean spread S = 0.87, the average relative estimation error is 0.11. The inversion of the measured data was also performed with the G_LSQ algorithm. As a start model, the logarithmic inversion variables were set to 0.15 again. After 40 iterations, the distance between the resulting model and the measured data was D = 0.000137, which is a significant improvement over the T_LSQ procedure. The linear spectrum obtained as a result of the inversion is shown in Fig. 5, where each tau interval was chosen as before. In Table 3. one can see the associated values (time constant, spectral density and estimation error) for each spectral line (greater than 0.001). The sample shows significant excitation at time constants of (5,10), (60,65) and (340,345) s, while the smaller amplitude obtained at T_LSQ does not occur at (110, 115) s.
The duplicate spectral lines are the consequence of the discretization we used. There is no reason to expect that a certain time constant in a real inversion procedure coincides with any of the pre-defined equidistantly sampled i values. Instead, can range between two neighbouring positions i ≤ ≤ i+1 and the inversion procedure gives spectrum lines at i and i+1 . The real amplitude w is shared between the two positions depending on the distance of from the fixed neighbours ( i , i+1 ) . We assume for the two spectral amplitudes giving the equivalent time constant as amplitude weighted mean  1 3 and the equivalent spectral amplitude w = w i + w i+1 . The time constant values estimated by the inversion method are significantly above 1 s, which means the physical polarization type is definitely metallic. The correlation of the G_LSQ inversion variables is indicated by the value of the correlation norm S = 0.76, while the average relative estimation error is 0.00246. It can be stated that  1 3 compared to T-LSQ this result shows a significantly more accurate inversion estimate with less correlation. The next sample selected for analysis is chalcopyrite, a highly polarizable copper ore from the Recsk basin (Sample 2). (The measured apparent polarizability changes in the interval [96.3, 7.5]%.) In what follows only the results of the G_LSQ algorithm will be shown. As a starting model, the same value was used as in the case of Sample 1 to all (logarithmic) inversion variables resulting in. According to Eq. (27), the distance of measured and calculated (on the start model) data is D = 15.13, which is again a very large value. After 40 iterations the data distance is D = 0.000901, a very good fit demonstrating also the initial model independence of the inversion algorithm. The fit in the data space is shown in Fig. 6.
The line spectrum obtained as a result of the D_LSQ inversion method is shown in Fig. 7. In our experience, there is no excitation mechanism in the time constant domain above 600 s. Since Q = 100 unknowns were chosen, the possible time constant intervals were equally defined as 6 s long. The spectral lines shown in Fig. 7 appear in the middle of these intervals. The associated values (time constant, spectral density, estimation error and relative error) for each spectral line (with amplitude greater than 0.001) can be seen in Table 4. The sample shows significant excitation at time constants of (18,24), (96,102) and (498,504) s resulting in three equivalent time constants at 20.3 (s), 100.0 (s) and 500.5 (s) with 0.111, 0.251 and 0.526 spectral amplitudes, respectively. The average correlation of the inversion variables is indicated by the value of the mean spread S = 0.87, the average Fig. 6 Decay curve of the measured (dots) and calculated (solid line) apparent polarizabilities found by using D_LSQ method (Sample 2) relative estimation error is 0.00725. It can be seen that the sample has a very large spectral amplitude at an extremely high time constant (500.5 s).
The third sample in our study is the moderately polarizable Sample 3 with galenite, sphalerite and enargite constituents (Recsk basin). (The measured polarizability changes in the interval [40.3, 0.16]%.) The starting model of the inversion is the same as in the previous cases. After 40 iterations the data distance is D = 0.000556, representing again a good fit as demonstrated in Fig. 8. The line spectrum obtained using the D_LSQ method is shown in Fig. 9. Preliminary investigations showed that there is no excitation mechanism  above 400 s. Since Q = 100 unknowns were chosen, the possible time constant intervals were equally defined as 4 s long. The spectral data values can be seen in Table 5. The sample shows only two significant excitation mechanisms at time constants of (20,24) and (312,316) s resulting in two equivalent time constants at 28.3 (s) and 313.7 (s) with 0.113 and 0.095 spectral amplitudes, respectively. The average correlation of the inversion variables is indicated by the value of the mean spread S = 0.89, the average relative estimation error is 0.00525. In the fourth example, we show a less polarizable Sample 4 with galenite and sphalerite constituents (Recsk basin). (The measured polarizability changes in the interval [16.6, 0.06]%.) With the previously applied starting model, the G_LSQ inversion results D = 0.00114 data distance after 40 iterations, representing a moderate fit in the data space. The line spectrum given by the D_LSQ method is shown in Fig. 10. As in the case of Sample 3, there is no excitation mechanism above 400 s, thus the constant intervals were again equally defined as 4 s long. The spectral data values can be seen in Table 6. The sample shows only two excitation mechanisms at equivalent time constants of 7.88 (s) and 189.3(s) with very low 0.0611 and 0.0368 equivalent spectral amplitudes, respectively.
From the results, it can be seen that by processing the apparent polarizability data by a series expansion based inversion method, the time constant spectrum can be adjusted with the appropriate resolution and accuracy by both the T_LSQ and G_ LSQ methods for practical use. Thus, the investigated methods are equally suitable for isolating polarization effects. The G_LSQ method shows significantly more accurate inversion results (with an order of magnitude for the average relative estimation error) with less correlation. The results found in the four samples prove that the G_ LSQ inversion method is applicable in a very broad range of polarizability. The more accurate and stable inversion was made possible by the calculation of the matrix of the system of normal equations by analytical integration. Thus, it can be concluded that the application of the objective function defined in Eq. (14) leads to a significant inversion advantage.

Summary
Time-domain induced polarization measurements are often used to solve ore exploration tasks. Numerous publications have been published for the field use of this method. In this paper laboratory, IP apparent polarizability data were processed using the series expansion based inversion method instead of decomposing them into previous exponential components. After inversion estimation of expansion coefficients, the time constant spectrum was calculated, which was plotted graphically to isolate several polarization effects. The time constant values found by the new inversion methods are significantly above 1 s, which proves that the physical polarization type is definitely metallic. The generalization of the  inversion objective function was proposed and by minimizing it the G_LSQ algorithm was defined. The traditional T_LSQ and the new G_LSQ algorithm were tested on laboratory data measured on four samples with different polarizabilities and mineral constituents. It was shown that the new G_LSQ method resulted in significantly more accurate inversion (one order of magnitude concerning the average relative estimation error) and less correlated model parameters.