Abstract
3D SIP inversion algorithm involves multiple parameters, and the key is the calculating speed and memory. A whole set of quasi-linear (QL) theories has taken shape in recent years, including the QL approximation method proposed by Zhdanov, quasi-analytic approximation, and localized quasi-linear (LQL) approximation. They are characterized by high speed and accuracy in electromagnetic field numerical modeling. The above-based 3D QL inversion algorithm, boasting quicker calculating speed plus more stable and favorable inversion effect, has been adopted profoundly in electromagnetic prospecting, whereas its frequent source conversion requires recalculating the dyadic Green’s function and primary field each time, thus delaying the 3D SIP modeling speed. This study makes use of the spatial symmetry in the primary field and Green function to propose an effective and quicker QL forward modeling method, which has the hallmark of higher calculating speed owing to less calculating times, and makes feasible the 3D SIP conjugate gradient inversion algorithm with Cole–Cole parameter range constraints.
1. Introduction
SIP, a branch of the induced polarization (IP) method, is capable of working out the resistivity and polarizability parameters and frequency-related coefficient, as well as a time constant. The latter two can distinguish 3D abnormal bodies based on the structure. This means a lot in identifying deep minerals, defining local enriched mineralized areas, and assessing the prospect [1]. Both electromagnetic effect and IP effect are observed in SIP. The previous SIP method overlooked the electromagnetic effect’s impact or separated it in the approximation formula. For the past few years, the SIP forward and inversion method, taking the frequency domain Maxwell equation and Cole–Cole model as basis, or to say, taking both electromagnetic effect and IP effect into account, has deepened the research and understanding on the electromagnetic effect by leaps and bounds [2–4]. Nevertheless, difficulties still exist in both theoretical study and practical application, of which the electromagnetic effect and 3D inversion are of the greatest concern [5–7].
A 3D electromagnetic method is impractical for actual application due to large memory consumption and slow calculating speed. The Born approximation method was first introduced to solve the integral equation, which assumed the scattered field as zero; then, the incident field of the abnormal body is the total field. By this means, the forward modeling time is shortened because of no need to solve the equation. However, the poor calculating accuracy and limited scope of application restrict the Born method only applicable for the condition of weak reflection [8]. The QL approximation method, proposed by Zhdanov and Fang compared its calculation result with that of the integral equation and Born method, is featured with higher calculation speed and accuracy, as well as less memory consumption. The 3D QL inversions of magnetotelluric data made by Zhdanov et al. observed 161 points at 7 frequencies in Minamikayabe, Japan, gaining a result similar to that of 2D inversion made by Takasugi et al., and it provided more geological structure details for geological interpretation of the aquifer area [9–11]. Zhdanov et al. [12], in 2006, made SIP nonlinear conjugate gradient inversion based on the Cole–Cole model and trial inversion calculation based on a theoretical model, but yielded little other than zero-frequency resistivity.
Compared with the above geophysical prospecting methods, the SIP method, by making high-density measurements in spatial domain and frequency domain with a stronger capacity of resisting disturbance, can obtain more electrical parameters and reveal richer abnormal data through multiparameter comparative interpretation. Nevertheless, with peculiar transmitting and observing patterns, it needs to keep altering the transmitting and receiving positions and changing the transmitting-receiving distance in the measurement of depth. The QL approximation method, if being used in 3D SIP forward modeling, needs repetitive computations of reflection coefficient and of the dyadic Green’s function plus recalculation of the primary field in each time of source conversion. In consequence, the forward modeling speed is horribly slow. Now that the speed and accuracy of calculation are decisive for the method’s application prospect, this study introduces the QL approximation method into 3D SIP forward and inversion modeling research by taking advantage of the spatial symmetry of the primary field and Green’s coefficient to explore a speedy QL approximation method. This method optimizes and quickens the algorithms and thus greatly promotes the accuracy of high-frequency spectrum electromagnetic modeling. It can effectively measure the complicated geo-electric structures to apply in large-scale IP data inversion. Section 2 discusses the forward calculation process of three-dimensional SIP. Section 3 introduces the calculation process and parameter determination of the SIP fast QL method. Section 4 gives a theoretical model inversion test.
2. QL Forward Modeling Theory
Figure 1 shows a 3D abnormal body on the homogenous Earth. The volume integral equation can be drawn from Maxwell equations [9, 10]:where refers to the total field, is the primary field, D is the anomaly region, j is the current element number, is the current element, r is any random element, is Green’s coefficient, and is the abnormal conductivity. Total field minus primary field is the secondary field:wjere refers to the secondary field.
QL approximation theory assumes that [11], in the anomaly region, the abnormal filed and setting field bear a linear relation. Plugged into equation (2), the approximate quasi-linear solution is worked out:where is the electric reflection coefficient, to be figured out by solving the minimum of the following equation:
In the QL approximation method, the reflection coefficient is calculated in correlation with abnormal conductivity, primary field, and Green’s coefficient. In particular, Green’s coefficient is most time-consuming in forward modeling due to the laborious calculation amount involved, whereas Green’s function spatial symmetry can be fully utilized to reduce workload and speed up the calculation for the reflection coefficient.
Dissect the abnormal body at the directions of x, y, and z in the number of elements of , , and , respectively. Green’s coefficient matrix at the direction of x can be expressed as
herein represents the field generated by element j as the source at the center of element i or to say Green’s coefficient of element j at element i. Assuming that the relative position of element i and j stay invariant, the value of should be the same like ; thus, this matrix conforms to the structure of the Toeplitz Matrix.
It is adequate to store the elements of the first row and line, by which the whole matrix can be expressed. Hence, it is sufficient in the quick QL approximation method to only calculate and store the elements of the first line in the matrix of equation (1), which, though claiming more memory, involves much less calculation amount than the QL approximation method, , approximately. Though each time of source conversion requires recalculation of reflection coefficient, Green’s coefficient needed is all the same. In consequence, it is feasible to calculate and store Green’s coefficient matrix only once to figure out all reflection coefficients in all source conversions, which is particularly applicable for multisource electromagnetic forward modeling [13].
Referring to the conclusion drawn by Pelton et al. [14], the complex resistivity of ores can be shown by the Cole–Cole model:where is the complex resistivity, , , , and are, respectively, the zero-frequency resistivity, polarizability, frequency-based coefficient, and time constant, , and the transmitting frequency.
The 3D SIP forward modeling with IP effect can be achieved by substituting the complex resistivity in equation (6) into equation (3).
3. QL Inversion Theory
Based on QL approximation theory [15], the secondary field at the ground receiving point is indicated aswhere refers to the medium’s property parameter:
The relation of and (reflection coefficient) is built after working out from equation (8), which is as follows in the abnormal region:
3.1. The Inversion of Matter Property Coefficient m
For a given setting conductivity , the secondary field can be calculated according to the observed electric field amplitude and its phase position:where represents the secondary field data, the theoretical data, the Green linear operator, the element primary field, the abnormal conductivity, and the reflection coefficient.
Due to the severe multiplicity of solutions in 3D inversion modeling, the smooth constraint inversion method is adopted to set an objective function as follows:where represents complex-valued conjugate, is the data variance vector, is the Lagrange factor, R is the smoothness matrix, and . Make derivation of from the above and take the minimum value to draw the inversion equation:where is the residual vector of observed data and forward theoretical data and is the model correction.
The fitting precision is defined aswhere is the fitting precision and is the quantity of data.
Given the large calculated amount involved in 3D inversion modeling, the conjugate gradient method is selected, which occupies relatively less memory, and the inversion process is as follows:(1)Preset the initial model and Lagrange factor (2)Calculate the right-hand item , making (3)Calculate , holding for later use(4)(5)Calculate the step size (6)Renew the model , figure up the fitting difference, and then jump out if it is up to the given accuracy or is less than the threshold value, or otherwise continue the cycle(7)(8)(9)(10)Make i = i + 1, and then skip to (4)
3.2. Determination of Cole–Cole Model Parameter
The inversion of the Cole–Cole model parameter is carried out element by element. The Cole–Cole model parameter of any random element is defined by the following equation [16, 17]:where , represents the setting complex conductivity, while is the abnormal complex conductivity.
The Cole–Cole model parameter can be worked out by conjugate gradient iteration in nonlinear equations at four frequencies. Since Cole–Cole model parameters are all positive real number changing within a certain range, the constrain function provided by Commer and Newman is introduced here [16, 18–20]:where represents the Cole–Cole model parameter, is the correction of the inversion parameter, and is the initial model. In this way, is limited within the range, which will effectively reduce the solution multiplicity in inversion modeling.
In brief, use the data at four frequencies to make conjugate gradient inversion of , then work out (complex conductivity) of each abnormal element through equations (8) and (9), and subsequently, make conjugate gradient inversion of the Cole–Cole parameter based on at four frequencies.
4. Theoretical Data Inversion Measurement
In this paper, we make the inversion measurement with an intermediate gradient for observing at the equator direction of the source. The relatively long transmitting source stimulates a stronger electromagnetic field, which facilitates signal detection and enhances the signal-to-noise ratio (SNR). Yet, since the transmitting source penetrates the measured area, the apparent resistivity will be heavily raised due to the severe variation of field intensity near the source plus the interference from the high-frequency electromagnetic effect. A theoretical model inversion test is to be carried out below to detect whether it can yield precise results from intermediate gradient data.
The model is set, as shown in Figure 2(a) (left), with a polarized abnormal body of low resistivity in reversed T-shape in the homogenous half-space nonpolarized setting. Under the resistivity of , the zero-frequency resistivity, polarizability, frequency-related coefficient, and time constant of the abnormal body are, respectively, , 0.5, 0.3, and 20 s. The 2000 m long given transmitting source is located at the line of the ground at the center coordinates of (0, 0, 0).
(a)
(b)
13 lines are set right above the abnormal body, with 41 measuring points at each line, 533 points in total. The line spacing and point spacing are both 40 m. The forward modeling is carried out at four frequencies, i.e., 1 Hz, 4 Hz, 16 Hz, and 64 Hz, with a given current of 10 A to stimulate the transmitting source. The apparent resistivity obtained from the equation is shown in Figure 2(b) (right), where IP abnormality can be identified at 1 Hz, and along with frequency increasing, the apparent resistivity of the whole region becomes higher, while the abnormalities vanish.
Convert the data of apparent resistivity and phase position from forward modeling above into electric field for inversion modeling. By subdividing the inversion region evenly in squares with the side length of 40 m in elements, the inversion modeling focuses on three aspects: (1) the feasibility of SIP inversion in intermediate gradient, (2) the impact of the selected initial model upon inversion result and the reason, and (3) the intermediate gradient QL inversion’s resistance against the noise. Refer to Table 1 for the correlation of inversion result with the initial model and noise and Figure 3 for the inversion result under different initial models with each random noise level.
(a)
(b)
(c)
As for intermediate gradient SIP inversion under the noiseless condition, if the selected initial model is close to the background model, it is very favorable to 3D inversion, as it only needs to fit the abnormity produced by the T-shape abnormal body. Otherwise, when it comes to the large difference between the initial model and the background model, large abnormal conductivity will arise in all subdivided elements. As shown in Figure 3(a), except for the area deeper than 100 m keeping initial status as a result of low sensitivity, most other areas can ultimately recover the real setting status. However, regarding T-shape abnormal body, the resistivity property can only be roughly recovered. As for Figures 3(b) and 3(c), for the result in noise condition, under 3% noise condition, the inversion result can vaguely present the abnormal property of resistivity, polarizability, and frequency-related coefficient, while under 5% noise condition, the inversion can only roughly regain the abnormal property of resistivity, but is unreliable in other parameters, for they are impacted by noise with some false abnormities arising. Known from the forward modeling data, the abnormal field arising from “T”-shape abnormal body accounts for a maximum of 8% in the total field. This means, in case the noise exceeds half of the maximum value of an abnormal field, the QL inversion will fail to recognize the abnormity. In other words, the existence of noise not merely causes false abnormities but also decreases the resolution ratio in inversion.
In light of this, the result of the conjugate gradient-based QL inversion is greatly influenced by the selection of the initial model, but this is inevitable in any kind of linear inversion modeling. Still, if with a proper initial model selected, it can achieve the desirable results, with quicker calculation speed besides.
As for the calculation speed, the quick QL approximation method, similar to the integral equation method, the Born approximation method, and NL approximation method, is positively related to the quantity of divisional elements, i.e., the more the elements, the longer the time needed for calculation. The time taken in the QL and in the quick QL approximation method are, respectively, listed and compared in Table 2:
As it shows, the quick QL approximation method, by taking advantage of the spatial symmetry of Green’s function, can conduct the forward calculation far more quickly than the QL approximation method, at an efficiency equivalent to that of hundreds of computers, thus making it feasible to carry out myriad 3D inversion calculations at an ordinary computer.
5. Conclusion
SIP inversion, involving enormous memory, tremendous calculations with a severe multiplicity of solutions, and especially complicated physical parameters, is a challenging task highly demanding for the performance of forward and inversion methods. This paper proposes a QL approximation method with higher calculating speed to lay a foundation for thorough research on multisource electromagnetic forward and inversion modeling and makes the 3D SIP inversion algorithm doable with the Cole–Cole model parameter range constraints.(1)In the QL approximation method, the spatial properties of the primary field and dyadic Green’s function can serve effectively to reduce the calculation amount and the time needed for forward modeling. Exploring more the physical essences of Green’s function and utilizing them will further speed up the QL approximation algorithm.(2)The inversion is conducted in two phases based on the characteristics of the QL approximation method and SIP method. First, adopt the material property parameter to convert the nonlinear inversion equation into a linear one, by which the inversion iteration speed gets much faster, and by adding the smooth constraint, the multiplicity of solutions is improved. Second, calculate based on the above and its relation with reflection coefficient and complex conductivity , and add the Cole–Cole parameter range constraint to fulfill the inversion.(3)As a theoretical trial for the data with intermediate gradient observation, the inversion algorithm can evaluate the all-around performance of QL inversion modeling through several representative model inversion results. It can quickly and steadily present the zero-frequency resistivity, polarizability, and frequency-related coefficients with the preferable inversion effect.
Since some calculation procedures can go ahead independently of each other in the process of QL integral equation forward modeling, in the next phase, we plan to utilize parallel algorithm to further improve the corresponding modeling speed.
Data Availability
The data used to support the findings of the study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the science and technology research project of Jiangxi Provincial Department of Education (GJJ191061).