Research on 3 DSIPConjugate Gradient InversionAlgorithmwith Parameter Range Constraints

3D SIP inversion algorithm involves multiple parameters, and the key is the calculating speed and memory. A whole set of quasilinear (QL) theories has taken shape in recent years, including the QL approximation method proposed by Zhdanov, quasianalytic approximation, and localized quasi-linear (LQL) approximation. (ey are characterized by high speed and accuracy in electromagnetic field numerical modeling. (e 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.(is 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.


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. e latter two can distinguish 3D abnormal bodies based on the structure. is 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. e 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][3][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][6][7].
A 3D electromagnetic method is impractical for actual application due to large memory consumption and slow calculating speed. e 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]. e 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.
e 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][10][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.
e 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. is method optimizes and quickens the algorithms and thus greatly promotes the accuracy of high-frequency spectrum electromagnetic modeling. It can effectively measure the complicated geoelectric structures to apply in large-scale IP data inversion. Section 2 discusses the forward calculation process of threedimensional SIP. Section 3 introduces the calculation process and parameter determination of the SIP fast QL method. Section 4 gives a theoretical model inversion test. Figure 1 shows a 3D abnormal body on the homogenous Earth. e volume integral equation can be drawn from Maxwell equations [9,10]:

QL Forward Modeling Theory
where E(r j ) refers to the total field, E b (r j ) is the primary field, D is the anomaly region, j is the current element number, r j is the current element, r is any random element, G E is Green's coefficient, and Δσ is the abnormal conductivity. Total field minus primary field is the secondary field: wjere E a 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 n x , n y , and n z , respectively. Green's coefficient matrix at the direction of x can be expressed as (1,2) . . .
G (i,j) 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 G (i,j) should be the same like G (1,2) � G (2,3) ; 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, 1/(n x n y ), approximately. ough 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 ρ(iω) is the complex resistivity, ρ 0 , η, c, and τ are, respectively, the zero-frequency resistivity, polarizability, frequency-based coefficient, and time constant, ω � 2πf, and f the transmitting frequency. e 3D SIP forward modeling with IP effect can be achieved by substituting the complex resistivity in equation (6)

QL Inversion Theory
Based on QL approximation theory [15], the secondary field at the ground receiving point is indicated as where m refers to the medium's property parameter: e relation of m and λ (reflection coefficient) is built after working out m from equation (8), which is as follows in the abnormal region: 3.1. e Inversion of Matter Property Coefficient m. For a given setting conductivity σ b , the secondary field can be calculated according to the observed electric field amplitude and its phase position: where d represents the secondary field data, E a ql the theoretical data, G E the Green linear operator, E b 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, W d is the data variance vector, λ L is the Lagrange factor, R is the smoothness matrix, and G � G E (E b (r)). Make derivation of m from the above and take the minimum value to draw the inversion equation: where Δd is the residual vector of observed data and forward theoretical data and Δm is the model correction. e fitting precision is defined as where RMS is the fitting precision and k 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 m 0 and Lagrange factor λ L (2) Calculate the right-hand item figure up the fitting difference, and then jump out if it is up to the given accuracy or λ L is less than the threshold value, or otherwise continue the cycle 10) Make i � i + 1, and then skip to (4)

Determination of Cole-Cole Model Parameter.
e inversion of the Cole-Cole model parameter is carried out element by element. e Cole-Cole model parameter of any random element is defined by the following equation [16,17]: where x � x where x represents the Cole-Cole model parameter, Δx is the correction of the inversion parameter, and x 0 is the initial model. In this way, x is limited within the (x min , x max ) range, which will effectively reduce the solution multiplicity in inversion modeling.

Mathematical Problems in Engineering
In brief, use the data at four frequencies to make conjugate gradient inversion of m, 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.

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. e 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.
e model is set, as shown in Figure 2 e 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. e apparent resistivity obtained from the equation R a � K MN (U/I) 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 E X for inversion modeling. By subdividing the inversion region evenly in squares with the side length of 40 m in x × y × z � 20 × 12 × 10 � 2400 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.
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. is 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. e 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.

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. is 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) e inversion is conducted in two phases based on the characteristics of the QL approximation method and SIP method. First, adopt the material property parameter m 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 σ(iω) based on the above m and its relation with reflection coefficient λ and complex conductivity σ(iω), and add the Cole-Cole parameter range constraint to fulfill the inversion.  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
e data used to support the findings of the study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.