Extraction and inversion of two-dimensional spectrum parameters based on TDIP data

Conventional TDIP can only obtain two parameters: apparent resistivity and polarizability, so it is difficult to interpret the effective information of ore bodies with high precision, especially the symbiotic and associated ore bodies with little difference in resistivity and polarizability, such as pyrite and carbonaceous limestone, while the spectral parameters (cm and τ) can quantitatively evaluate the abnormal characteristics of mineralized bodies and are rich in geoelectric information. Based on the field apparent resistivity observation mode of Cole-Cole model, the charge rate values of symmetrical quadrupole arrangement of two typical mineralized bodies are calculated by finite element method forward simulation from the 2D point source, this paper puts forward a method for extracting the corresponding spectrum parameters from the chargeability value of TDIP data and verifies the effectiveness of the method.The results show that the abnormal shape is obvious and is consistent with the actual model parameters, which can reflect the spectral parameter characteristics of each abnormal body, realize spectral imaging, and provide a favorable basis for further division of polarization abnormal bodies. This method has been applied to the TDIP data processing of a lead-zinc mine in Wuzhou City, Guangxi Province and good results have been obtained.


Introduction
Pelton researches that the complex resistivity spectrum caused by IP effect satisfies the Cole-Cole model by measuring various rocks and ores [1].Geophysicists successfully evaluate IP anomalies and distinguish the abnormal bodies that cause IP anomalies by making use of the differences spectral parameters [2,3,4,5,6].In recent years, with the deepening of research, the application field of spectrum parameters has also developed rapidly, especially in the geotechnical field, Titov [7] and Revil [8], have found that the spectrum parameters are related to the conductivity of water in rock and soil, capillary pressure and so on. The extraction of Cole-Cole model parameters has always been a hot issue in electromagnetic method, and researchers have done a lot of research on it [9,10].The spectrum parameters are extracted from the complex resistivity data, but the SIP needs to observe many (at least [12][13][14][15] complex resistivity data, and the observation time is long, the field efficiency is very low, and each frequency point is affected by electromagnetic coupling. It is not widely used in production. Relatively speaking, the TDIP observation time is short, the efficiency is high, and it is basically not affected by electromagnetic coupling after a period of power outage [11,12].Therefore, the universality of domestic application is stronger, and the abnormal shape of time-domain apparent chargeability is similar, so only selects any one of them to explain, resulting in a waste of resources.At present, in  [15]et al inverted the Cole-Cole model to extract the apparent spectrum parameters from the 1D and 2D polarization data analysis.Wei Zhong [16] and others calculate forward modeling directly based on the Cole-Cole model in the time domain. The effects of increasing time window and emission waveform are taken into account when calculating the chargeability, and the spectrum imaging is realized simply and quickly. In this paper, the chargeability of two typical mineralized bodies under symmetrical quadrupole arrangement is calculated by finite element forward modeling based on the Cole-Cole model.The corresponding spectrum parameters (c m and τ) are extracted by least square inversion from the chargeability. The abnormal morphological characteristics of the apparent spectrum parameters are analyzed and applied to the field prospecting for distinguish the polarization abnormal bodies and improve the exploration resolution.

Forward calculation method of 2D chargeability based on Cole-Cole model
According to the charge rate spectrum calculation method of Cole-Cole model, the TDIP chargeability msi is defined, and the arrangement coefficient K and current I are reduced: In order to calculate the chargeability(msi), we can first calculate the resistivity ρ0 without polarization effect;Then the equivalent resistivity ρ*(ξi) of the first observation integral period is solved;Finally, the chargeability(msi) of the I observation integral period is obtained.In order not to affect the calculation speed and efficiency, the time spectrum resistivity of each grid node can be calculated first, and then the equivalent resistivity of the i observation period on each node can be calculated. the equivalent resistivity of the I observation period and the resistivity without polarization are calculated by resistivity forward simulation, and the calculated results are brought into (2) to calculate the chargeability of the i period. The finite element forward simulation uses a given initial model and triangulation to calculate the resistivity of each node.The number of mesh points in X direction is 145 and that in y direction is 50. In horizontal direction, the spacing of division is as follows: 64, 32, 16,8,4,4,4,4,4,4,2,2,2,2,2,2,2,2,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, Then calculate the position of the model element in the grid, put the model element into the grid that has been divided, and then get the true spectrum parameter msi.

Extraction method of apparent spectrum parameters
The spectral parameters are extracted by least square method. The specific steps are as follows: ①It can be summarized as minimizing the following objective function： In formula (3), NP represents the number of chargeability；msi indicates the charge rate of the i observation period；X = {c,m,τ} represents the spectral parameters in the model;mli(X) indicates that when the model is X, chargeability of the i observation period is calculated.
②If you expand it at X by Taylor, and ignore the terms with more than quadratic terms, you can get : Where NM is the number of model parameters；ΔX is the modification of the model. In order to achieve the minimum value of the above equation, it is necessary to satisfy the following linear equation: The modification ΔX of the prediction model X is calculated by solving the above equation, and the initial model X*=X+∆X is modified and iterated until the accuracy requirement is met.The partial derivative of equation (5) is calculated, the difference method is adopted, the equation is dimensionless, and a new parameter model is obtained by singular value decomposition, so as to achieve stable convergence.  [17]. The above method is used to deduce the spectrum parameters and compare them. Figure 1 shows that spectrum parameters extracted from TDIP data are the same as from the SIP data, and the first and tail branches tend to the true spectrum parameters of the first layer and the second layer, and the values are relatively stable. The jitter in the middle transition section of the spectral parameters extracted by TDIP is slightly larger, which is due to the poor convergence of inversion due to the lack of data in inversion calculation, but it can reflect the information of spectral parameters of each layer, indicating that the program is accurate.

2D model
Aiming at two common mineralized body models, using the above principles and methods to carry out forward simulation using symmetrical quadrupole arrangement, it is calculated that the delay period of chargeability in four observation periods of each node is selected as follows:  Table 1.

Figure3. Pseudo-section map of inversion spectrum parameters about upright model
In the inversion spectral parameter pseudo-section figure 3, cs, ms and τs all present a closed circle corresponding to the model parameters below the geometric center of the plate, and the abnormal range is slightly "flat and thin", which is similar to an upright "rectangle". Its "axis" (anomaly center) is about 9 meters near the center of the polarizing body.. In the spectral parameters, the fitting accuracy of ms is higher, and the fitting accuracy of cs and τs is slightly lower, which leads to the worse shape of the anomaly. The cs surface map, there is a widening of the anomaly value at 12-16m at the anomaly center, and there are sporadic low values on the right side of the τs section. This is related to the division mode of dense in the middle and sparse on both sides when dividing the center of the grid.
Figure4. Pseudo-section map of inversion spectrum parameters about two rhombic polarized model Figure 4 shows that the relative size of cs,ms,and τs the relative position of anomaly shape can be distinguished. The closed circle below the abnormal body is clearly visible, the anomaly is obvious, and corresponding to the model parameters, it can reflect the corresponding geoelectric information, and the spectral parameter type is similar to the shape of "lentil". The information of spectrum parameters can be inverted by using the msi of different periods, which can provide more abundant

Application of measured data in TDIP
IP exploration of a lead-zinc mine in Wuzhou City, Guangxi Province is used as measured data.There are 20 sounding points, the distance between the points is 20m, and the direction of the sounding, line is east-west.Abnormal area is more obvious, the upper part basically is high , and the middle part is low. The all range is large and not obvious for abnormal, so it is difficult to distinguish the mineralized body. Figure 5 shows that it is low when the point number is 2200-2300m and depth is 100m from the cs and τs pseudo-section,however it is high from the ms pseudo-section, so it can speculated to be a lead-zinc mineralized body. However, when the point number is between 2320 and 2400m with depth below 100m also shows high polarization anomaly from the ms pseudo-section and it is very low from ts pseudo-section, it is not obvious in the cs pseudo-section. In summary, The range of abnormal II area is reduced, and the exploration accuracy is improved. However, attention should be paid to the influence of regional carbonaceous strata,through the verification of boreholes, the location of boreholes is located at 2300m. The verification results show that the surface layer is pyrite-mineralized and lead-zinc body and the lower part is carbonaceous layer, which accords with the inferred results.More abundant geoelectric information is obtained by frequency spectrum parameters.

5.Conclusions
 The values of resistivity and charging rate under symmetrical quadrupole device are calculated by forward simulation of two common mineralized body models, and the parameters of Cole-Cole model in frequency domain are extracted by least square method. The abnormal shape is obvious and corresponds well with the model parameters. The correctness of the program is verified by one-dimensional comparison, and two examples show that the abnormal shape of the low resistivity plate body is approximately "upright rectangle", and the rhombic mineralizatized body is "lentil-shaped", and the anomaly is clear and obvious, which provides a favorable basis for further division of polarization abnormal bodies.  The inversion effect of this method is good in practical application and improve exploration ability .However, inversion value may produce large errors only four apparent chargeability are used to extract apparent spectrum parameters cs, ms and τs, and there are observation errors and interference errors in the field, Fortunately, at present, TDIP can observe the secondary potential attenuation and apparent chargeability values in multiple periods, so it is only necessary to ensure the quality of field data and enough data. The more accurate apparent spectrum parameters cs, ms and ts can be extracted.