Feature Extraction and Intelligent Identification of Induced Polarization Effects in 1D Time-Domain Electromagnetic Data Based on PMI-FSVM Algorithm

The transient electromagnetic method can obtain resistivity and chargeability simultaneously in polarizable medium detection. Typically, we assume that the earth may contain a chargeable medium if the electromagnetic (EM) data appear negative values or sign reversals. Unfortunately, with barely perceptible characteristics, some EM responses with the induced polarization (IP) effects are considered to be non-polarizable responses. Insufficient understanding of features and inaccurate identification of the IP responses limits the use of the IP effects for broader purposes. For these reasons, we perform 1D forward modeling to discuss the degree of EM response affected by the IP effects and to extract polarization characteristics. To identify the IP effects, we combine partial mutual information (PMI) and the fuzzy support vector machine (FSVM) methods to complete the intelligent identification algorithm. We verify the efficiency and practicality of the algorithm by building Debye loops in field experiments. From the analysis, we distinguish the strong and weak IP effects by introducing the impact ratio. The strong IP responses manifest fast decays and sign reversals, and the weak IP responses primarily show fast decays or outward concavity. The identification algorithm validation results show that the recognition accuracy reaches 90.7%. In the field experiment verification, the Debye loop successfully simulates the IP effects of different intensity, and the identification results indicate that the algorithm has potential in the measured data. With this intelligent identification algorithm, the measurements can provide access to the weak polarizable medium when the impact ratio exceeds 30%.


I. INTRODUCTION
The induced polarization (IP) effect is a crucial electrochemical phenomenon that exists primarily in the dispersive medium with high economic value, such as disseminated sulfide minerals, some aqueous media, and hydrocarbon resources. The fast development of the multiparameter interpretation method or inversion extends the application of the IP data to a broader area, such as the detection of metal mines, The associate editor coordinating the review of this manuscript and approving it for publication was Gustavo Callico . water resources, geothermal resources, and oil resources, which has made it a heavily researched issue [1]- [3].
The transient electromagnetic (TEM) method is a critical time-domain sounding method with notable performance in the search for the polarizable targets. The TEM has many advantages: the secondary electromagnetic (EM) field measured is not affected by the ground electrodes, and it contains both resistivity and chargeability information. However, the TEM system operates at much higher frequencies than the conventional IP method, making it insensitive to some economic sulfide deposits, which have time constants of many seconds or more [4]. Initially, in the TEM detection, the IP effects raised attention because the EM data revealed sign reversals or negative values in later time [5]- [7]. Therefore, the crucial step is to analyze the impact of the IP effects on the EM data under different parameters and to extract characteristic parameters of IP responses. The influences of polarization parameters on the time of sign reversals and the maximum absolute value of the negative response were studied intensively in [8]- [9]. An overall analysis reported the influences of the polarization parameters, size, and buried depth of the target on the negative responses. They also demonstrated that only moderately sized, highly-polarizable media in a shallow layer could generally be detected [10].
Furthermore, in the field experiment, the traditional way to identify the IP effects is to investigate the negative responses or sign reversals. Unfortunately, IP decay curves are easily affected by polarization parameters and earth topography, or even system parameters [11], [12], and usually display more complex polarization characteristics. For example, a survey for the IP effects in airborne TEM data (AEMIP) found that the measured EM data do not have sign reversals [13]. If we continue to identify these IP data by the appearance of sign reversal, that will lead to subjective conjecture and inaccuracy in the identification of polarizable media. After exploring the features of the IP effects on EM data, one survey found that the EM data might show only one indicator, i.e., steep decay, which can enable the AEMIP inversion [14].
Machine learning (ML) is a collection of a variety of algorithms (e.g., neural networks, support vector machines, self-organizing mapping, decision trees, random forests, case-based reasoning, and genetic programming) [15]. Among these algorithms, neural network and support vector machine (SVM) methods have been successfully applied to the classification of underground rock structures [16], geophysical parameter estimation [17], [18], and even in earthquake detection and localization [19]. EM responses can be divided into two types: with or without the IP effects, which belongs to the typical classifier category. In 1995, the SVM method was first proposed in [20], and the SVM method has many unique advantages in solving nonlinear high-dimensional pattern recognition and delivers outstanding performance in real-world classification problems and estimation problems. Traditional SVM has developed to more precise models, such as the least-squares support vector machine [21] and the fuzzy support vector machine (FSVM) [22]. We apply the FSVM method for classification problems with outliers or noises because it is robust.
Moreover, if we consider a high number of variables simultaneously, the cost in time and precision of the FSVM is substantial; therefore, we usually select the main variables from numerous variables. Partial least squares, principal component analysis, and Schmidt transform can extract critical features from all variables to reduce dimension and eliminate correlations [23]- [25]. However, these methods perform poorly in solving nonlinear problems. The partial mutual information (PMI) method is a method derived from the mutual information (MI) concept, which uses MI values to quantify the correlations between variables. Compared with the MI method, the PMI method can not only eliminate the coupling between multiple input variables but also reduce the data dimension and improve the learning algorithm's performance, especially for the nonlinear problem. The method has been applied to nonlinear variable selection in artificial neural networks and achieved excellent performance in the learning and generalization process [26]. Therefore, we currently apply an intelligent method to identify electromagnetic data with IP effects. We perform 1D forward modeling to simulate the EM response based on the homogeneous half-space, three-layered, and ten-layered models with different electrical properties. We divide EM curves into roughly three types and further summarize and supplement the polarization characteristics. By analyzing different impact ratios, we define the strong polarization response and weak polarization response. With the initial input, we select the optimal characteristic parameters and eliminate the redundant variables using the PMI method. Finally, we establish and train the FSVM model for IP response recognition and verify the algorithm's accuracy and feasibility through numerical test and field simulation experiments.

II. 1D FORWARD MODELING FOR THE IP EFFECTS
To analyze the polarization characteristics, we perform a 1D forward modeling method to numerically simulate the EM responses induced by the central circular loop. The transmitter transmits a step waveform current. When the transmitter and receiver are located on the ground surface, the frequency-domain equation of the vertical component of the EM response at the center is where H z is the vertical magnetic field intensity, I is the amplitude of the transmitter current, a is the radius of the transmitter coil, λ is the integral variable, J 1 (λ) is the first-order Bessel function, Z 0 is the wave impedance of the surface medium, and Z (1) is the total wave impedance. First, we calculate the solution in the frequency-domain, and then we obtain time-domain responses through the frequency-time transform method. Conventional frequency-time transform methods include sine, cosine transform, and the Guptasarma filtering method.
To describe the dispersion characteristics of the earth, we introduce the Cole-Cole model [4]. The conductivity expression in the frequency-domain is where σ ∞ is the conductivity at the infinite frequency and τ is the time constant. Moreover, η and c are the chargeability and frequency dependence, respectively, both of which range VOLUME 8, 2020 from zero to one. When c = 1, the Cole-Cole model reduces to the Debye model. Fig. 1 shows the equivalent circuit for a polarizable medium induced by the inductive source, where emf represents the induced electromotive force generated by the primary field during the switch-off time. In the field experiment of Shaoguo, in northeastern China, we also built a circuit loop based on the same loop to test the algorithm. The results demonstrate that when η is larger, EM data shows sign reversals; when η is smaller, EM data only shows changes within slopes in the late time. The appearance of sign reversal also depends on c and σ ∞ . We are not going to talk about how the Cole-Cole parameters affect EM decay curves. We try to summarize the characteristics of EM curves with the IP effects. Our analysis reached a similar conclusion comparing with the AEMIP survey that the EM data might show the only indicator of the IP effects, i.e., steep decay [14]. Additionally, we also find that some curves show an outward concavity (inset figure of Fig.2 (a)). Thus, we divided the curves into three types by the changing of slopes: curves with (A) steep decay and sign reversal; (B) steep decay but no sign reversal; (C) an outward concavity. Unfortunately, the latter two types of IP responses are usually regarded as non-IP responses, which leads to inaccurate interpretation. The IP effects of the rock are affected by many factors, such as mineralogy (especially the concentration of metal particles and clay minerals), porosity and tortuosity, pore saturation, and pore water salinity. The chargeability characterizes the intensity of the IP effects within the rock. However, the chargeability can only quantify the physical characteristics of the local medium specifically and cannot reflect the impact degree of the IP effects on the EM field comprehensively. Thus, we define a parameter to demonstrate an overall impact degree of the IP effects on the EM data, including the influences of earth topology, polarization parameters, and measurement system. The total field can be regarded as the sum of the fundamental EM field and the polarization field [27]. Therefore, to qualify the impact degree of the IP effects on EM data, we define the impact ratio below.
where Q is the impact ratio, V Fund represents the fundamental EM responses, and V Total represents the total EM responses that include the IP effects.
In this study, we present the impact ratio of the previous three types in Fig. 3, 4, and 5. From Fig. 3(a), we know that V Fund is approximately a straight line on the log-log axis with a negative slope. The slopes of V Total are negative in the early time, and its absolute value is increasing with time. At 0.6 ms, sign reversal appears, and at 0.9 ms, the absolute value of negative response reaches the maximum. Then, the slopes change to positive values. Q increases with time, and during  the sign reversal, it is close to 100%. After 1 ms, the impact ratio soars to 150%. For type B (Fig. 4), with the increment of time, the absolute value of the slope also increases, and Q increases from 35% to 75%. For type C (Fig. 5), the absolute value of the slope of V Total first increases and then decreases, which reveals an outward concavity. Q rises from 25% to 70%, approximately. The impact ratio of type A exceeds 100%, even twice the other two types, so we call it a strong IP response. The decay curves of type B and C are less affected by the IP effects and show polarization characteristics indistinctly; therefore, we call them weak IP responses. Actually, from the analysis of the decay curves, we know that the slopes and sign reversals are the crucial IP characteristics. Thus, apart from the characteristic parameters about sign reversals such as time of sign reversal, the maximum absolute value of negative responses [7], [10], and the duration of negative values, we also need to pay attention to the slopes of curves.

B. SLOPE FITTING BASED ON PIECEWISE NONLINEAR LEAST-SQUARES METHOD
Decay curves (−dB z /dt or B z ) in a nonpolarizable homogeneous half-space approximate a straight line under the log-log coordinate axis. These curves agree with the form of y = bx k , where k represents the slope under the log-log axis [28]. We apply a piecewise nonlinear least-squares method to fit the slopes because it is a nonlinear problem. From the theoretical analysis, we know that the slope of strong IP responses is negative before the maximum absolute value of negative responses; therefore, we extract only the negative slopes. If the curves have multiple sign reversals, we care mainly about the first one because it has the largest amplitude. We first locate the sign reversal, then record the time of the sign reversal and the maximum absolute value of negative responses, and then set the end time T end at the same point. If the curves do not have any sign reversal, we set T end at the end of the curves. The whole fitting time is where T initial is the initial time. T initial takes 100 µs for three reasons: first, most sign reversals appear after 100 µs; second, due to the limitation of the sampling frequency of the receiver, few points are recorded before 100 µs; third, the fitting error of early time is undesirable because t > t L = 100a 2 /2π10 7 ρ does not meet, where a is the radius of the transmitting coil, and ρ is the resistivity [29]. Then, we divide the T f into several segments and perform a piecewise fitting algorithm to fit the changing slopes.
The length of time window t directly affects the fitting accuracy, especially in the curves with sign reversal because they change rapidly. Figure 6 shows the fitting solutions when t takes different values. An exceptionally longer time window such as 1500 µs, 1000 µs, and 750 µs can make the slopes obtained inaccurate and may lead to failure in extracting IP features. In the fitting process of the sharply changing curve, the shorter the time window is, the better the fitting precision is. However, an exceptionally short t may create a large amount of data. To guarantee the fitting accuracy and efficiency, t = T f / n, where n is the number of parameters we desire, and t is proportional to T f . For example, if T f is 2.5 ms, when we set n as ten, the time window is 250 µs.

III. PARAMETER SCREENING AND CLASSIFICATION METHOD A. SCREENING OF CHARACTERISTIC PARAMETERS BASED ON THE PMI METHOD
Shannon first proposed the information entropy to describe the quantity of information carried by variables [30]. Furthermore, Mutual Information (MI) represents the information shared by variables and measures the interdependence between two or more variables. It has an excellent performance in parameter selection for both linear and nonlinear variables.
Usually, for the given samples of X (x 1 , x 2 , . . . , x n ) and Y (y 1 , y 2 , . . . , y n ), where x n and y n are the n th observed data, and the probability density distributions are unknown. I (X, Y) can be calculated by the probability density estimation.
where f (·) represents the estimated probability density function. The nonparametric estimation method provides an approach for unknown distribution. Following [26], we applied the kernel density estimation method to obtain an estimation of probability density. The kernel function is the Gauss kernel function.
For a multi-input system, given random input variables X 1 and X 2 , the output is Y . The variable that corresponds to the maximum value of MI will be selected. Assuming that X 2 is chosen, if there is a correlation between X 1 and X 2 , the solution of MI of the next round will result in bias. Referring to [26], we also use a Partial Mutual Information (PMI) Method to exclude the correlation between variables by calculating conditional expectation m X1 (X 2 ) and m Y (X 2 ), which can effectively improve the accuracy of variable selection.
wherem Y (X 2 ) denotes the regression estimator and n is the number of observed data. An estimatorm X 1 (X 2 ) can be similarly constructed. Therefore, the residuals u and v can subsequently be obtained using the expressions below and they work as the inputs in the next loop.
The Akaike Information Criterion (AIC) is a standard for evaluating the performance of fitting statistical models [31]. The following equation seeks the best balance between the complexity and the fitting ability of the model.
where n is the number of observed data, u j is the residual error, and p is the number of variables. When T AIC reaches the minimum value, the set of optimal independent variables is screened out.

B. CLASSIFICATION MODEL BASED ON THE FSVM METHOD
The EM responses can be divided into two types: with or without the IP effects, which belongs to the typical classifier category. The FSVM method can effectively improve the prediction precision in binary classification problems with noise and outliers [32], [33]. For the FSVM model, the training set is where x i ∈ R n , y i ∈ [−1 1], and ξ i ≤ s i ≤ 1, where ξ is a sufficiently small positive number, s i represents the fuzzy membership degree, and l is the number of observed data.
In the binary classification problem, if the IP effects exist, y i takes the value 1, otherwise −1. To begin with, we map data from sampling space to a higher dimensional space by the kernel functions ϕ (·). Thus, the nonlinear problem is converted into a linear divisible problem.
The Lagrange function for the optimization problem is: where w is the weight vector, b is the offset variable, C is the penalty factor, and ξ k is slack variables, and α = (α 1 , α 2 , . . . , α l ) T ∈ R l + is the Lagrange multiplier. Applying the Karush-Kuhn-Tucker (KKT) condition equations, we solve the dual problem. The flow-process diagram of the classification of the IP effects based on the PMI-FSVM method is shown in Fig. 7.

A. VARIABLES SCREENING
One-dimensional forward modeling is performed in 880 models that contain 440 homogeneous half-space and 440 threelayered models. Here, the second layer of the three-layered model is polarizable, and its parameters vary like the halfspace model. The first and third layers are nonpolarizable media, with the conductivities fixed to 0.01 S/m. The transmitter is a square loop of 40 m × 40 m, and the current is a 50 A step waveform. We select 200 values of conductivity, i.e., 0.001: 0.001: 0.2 (S/m) for the nonpolarizable half-space models and the middle layer. The number of chargeable halfspace and chargeable layered models is 240, respectively. We choose three values for conductivity, i.e., 0.002, 0.02 and 0.2 S/m, four values for chargeability, i.e., 0.2, 0.4, 0.6, and 0.8, five values for time constant, i.e., 0.1, 1, 10, 100, and 1000 ms, and four values for the frequency component, i.e., 0.2, 0.5, 0.7, and 1. The depths of the first and second layers are 100 m, while the depth of the last layer is infinite. In total, 480 samples are affected by the IP effects, and 400 samples are free from the IP effects. The responses for -dBz/dt are calculated.
Through the analysis of polarization characteristics, we extract the maximum absolute value of negative responses as x 1 , time of sign reversal as x 2 , duration of negative responses as x 3 , and the ten slopes along with the time as x 4 , x 5 , . . . , x 13 , respectively. The values for each parameter are normalized to the [0 1] range. To screen the main characteristic variables, we use the dataset of 880 samples as the input data. After three rounds, we finally get the solution until T AIC no longer decreases. A part of the results is shown in Table 1. From the first circle results, we can see that PMI between x 13 and y achieves the maximum, and T AIC is −829.29. Thus, the slope in the latest time is an essential factor in identifying the IP effects. The parameter x 13 is added to set S. In the second round, the PMI between x 2 (time of sign reversal) and y achieves a maximum with T AIC decreasing to −2404.32. Thus, x 2 is added to set S. Then, in round three, x 10 cannot be added to S though PMI between x 10 and y is the highest because T AIC no longer decreases. Hence, the time of sign reversal (x 2 ) and the last fitting slope (x 13 ) are screened out by the PMI algorithm.

B. RECOGNITION RESULTS OF THE TEST SET
Next, we test the accuracy of the algorithm by adding another 400 samples of the ten-layered model. The depths of the previous nine layers are 50 m, while the depth of the last layer is infinite. We randomly select the Cole-Cole parameters of each layer by a program within the range of common values, VOLUME 8, 2020 where σ ∞ ranges from 0.0001 S/m to 0.1 S/m; η ranges from 0 to 0.8, τ takes 0.1 s, and c takes 1 or 0.5. Thus, we randomly divided the 1280 samples into a training set and a test set. The radial basis function (RBF) maps the nonlinearinseparable data to a high dimensional space. In choosing σ in RBF, a large value will lead to terrible fitting results; however, a minimal value will result in overfitting. In the training process, we select σ =0.01 to achieve a relatively good performance. The optimal set of variables S obtained by the PMI method is the input set, and whether the IP effect occurs or not is the output. To describe the performance of the algorithm, we apply the accuracy index P A : where z i is the number of samples that have been recognized correctly, and z m is the total number of samples. Fig. 8 shows the results of the recognition of EM response with or without the IP effects. The black curve in the inset figure shows the classification hyperplane that demonstrates the boundary between non-IP responses (−1) and IP responses (+1). Inside the hyperplane are the non-IP responses, whereas outside the range are IP responses. The slope of non-IP responses ranges from −2.7 to −2. From Table 2, we can determine that the accuracy of the algorithm is 90.7%, which can make a good performance in classification. The algorithm performed notably in the non-IP medium samples recognition with the accuracy reaching 94.7%, whereas a small imperfection in the identification of IP responses with an accuracy rate of 86.5%. Therefore, we further investigate the reasons for errors in the IP response category. We extract the impact ratio of the synthetic model to analyze the relationship between identification accuracy and the impact degree. Among 319 IP samples, the strong IP effects (impact ratio > 100%) account  for 38%, and they are recognized correctly; the weak IP effects account for 62%. From Fig. 9, we can see these misrecognized samples (approximately 13%) mainly concentrate on the weak IP effects, especially in cases when the impact ratio is below 30%. The results manifest that the algorithm can successfully identify the EM response affected by the IP effects greater than 30%.

C. RECOGNITION RESULTS OF TEM MEASURED DATA
Shaoguo town is located in northeastern China, where the human and external electromagnetic noises are at a low level, and we can record high SNR data. Thus, we lay a central transmitter-receiver system and conduct TEM experiments here. The transmitter is a square loop of 50 m × 50 m with an amplitude of 10 A. The transmitter waveform is a bipolar trapezoidal waveform, in which the transmitter frequency is 12.5 Hz, the duty cycle is 50%, and the switch-off time is 150 µs with the linearity of 94%. The peak-to-peak values of the maximum input voltage and the receiver's noise floor are ±5 V and ±50 µV, respectively. The current SNR is 100 dB, and it is proportional to the square root of the number of stacking. Thus, we apply 128 times stacking to guarantee the SNR can reach 120 dB. The location and TEM device settings are shown in Fig. 10. When carrying out field experiments, we sometimes face a lack of accurate geologic information, and we also need to consider efficiency limitations. Therefore, we first do preliminary method verification, which can narrow the gap between theoretical experiments and field experiments. In this experiment, we propose a method to simulate the EM responses with strong or weak IP effects by building the Debye loop in hardware. We wind a ten-turn coil with a radius of 0.25 m, and the self-inductance is about 0.1 mH. We also use adjustable resistors (R 0 and R 1 ) and a capacitor (C) to build the Debye model. Then, we connect it in series at two ends of the multiturn coil. To reduce the distributed resistance of the coil, the material of the wire is silver, and the resistance is 0.25 approximately. In this way, its internal resistance can be neglected compared with R 0 and R 1 . In addition, the distributed capacitance between the current loop and the ground, turn and turn is roughly the pF level and even smaller, and they are almost negligible. The device allows us to set different parameters according to different media. Thus, we build loops of strong IP effects and weak IP effects. Parameter settings for the strong IP effects: R 0 = 20 , R 1 = 100 , C = 100 µF, τ = (R 0 + R 1 ) C = 0.006 s, and η = R 0 / (R 0 + R 1 ) = 0.167; and for the weak IP effects: R 0 = 40 , R 1 = 400 , C = 50 µF, τ = 0.022 s, and η = 0.091.
Based on the superposition principle [27], we regard the responses of the ground and the Debye loop as the background response and the IP response. The total field measured is the superposition of the two responses [34]. We must note that in real detection of the target, the total response is inseparable, and the total response is no longer the responses of a homogeneous half-space or a layered structure specifically. We first measure the background response B 0 without the Debye loop using the HT dc SQUID sensor, then we measure the total response B T with the Debye loop coupled in the center. The HT dc SQUID sensor's noise level is 100fT/ √ Hz@10KHz, the slew rate is 30 mT/s, and the bandwidth is 100 kHz. The magnetic field response generated by the Debye loop can be decoupled by subtracting the background field from the total field, that is, B D = B T − B 0 . Fig. 11 and 12 show the responses measured for the Debye loop 1 and Debye loop 2, where black, blue, and red curves represent the background responses, IP responses, and the total responses, respectively. Figures 11 and 12 show the strong IP responses and weak IP responses when the impact ratios are approximately 200% and 75%, respectively. The total response in Fig. 11 decays rapidly in late time and changes to negative values. The background response gradually decays to the noise level. After subtracting the background field from the total response, the pure IP response reveals polarization characteristics of fast decay and sign reversal. In Fig.12, the total response decays to the noise level, but no sign reversal appears. However, after subtracting the background field from the total response, the pure IP response reveals an ideal polarization curve with a sign reversal. We extract the polarization characteristics parameters from the six curves as the input and show FIGURE 11. TEM measurement responses with (red curve) and without (black curve) of the test Debye loop 1. The difference is pure IP response (blue curve).  identification results in Table 3. The FSVM algorithm classifies these curves correctly. The total responses are identified correctly because their impact ratios are larger than 30%. The background responses, which are free from the IP effects, are classified as non-IP responses. The results verify the feasibility of the FSVM recognition algorithm. The experimental results are highly consistent with the simulation results. The method has the following advantages: 1) it can simulate the IP effects with different polarization intensities; 2) it can be coupled to the total field through the superposition principle; 3) it provides a verification approach and narrows the gap between theoretical experiments and field experiments.

A. INFLUENCE OF CONDUCTIVITY DISTRIBUTION OF THE 1D SYNTHETIC MODEL
We found that the algorithm can successfully identify more than 80% of weak IP effects and 100% of strong IP effects. However, we also found that some errors are not from the polarizable samples but the nonpolarizable samples. We have to explore the factors that lead to identification errors. In the synthetic models, the parameters of each layer are selected randomly. We extract each layer's conductivity of the nonpolarizable ten-layered models that are wrongly and correctly identified, respectively, and show them in Fig. 13. The conductivities of the near-surface layers (when the depth<200 m) and deep layers (> 400 m) in (a) and (b) are very close. The only reason for the errors is the conductivity distribution between the depth 200 m and 400 m. Taking (b) as a reference, the mean of conductivities in (a) is smaller than that in (b), and the cyan dots skew toward the high resistivity at the same depth. Fig. 13(a) can roughly be divided into 3-layers, i.e., low-resistivity layer, middle-high-resistivity layer, and high-resistivity layer. The comparison demonstrates that the conductivity distribution is a critical error factor. If a resistive nonpolarizable medium appears in the medium depth (200 m -400 m), the EM data also show a fast decay in the late time, which may be erroneously recognized as IP responses.

B. INFLUENCE OF SNR
The SNR of the instrument has a significant influence on the precision of the signal we measured. In extreme cases, the actual measured SNR is less than 40 dB. In an excellent experimental environment, the SNR can reach more than 120 dB. We add Gaussian noise signals with different SNR  according to this range to simulate the measurement signals. Then, we extract the characteristic parameters to test the accuracy of the identification algorithm. Figure 14 shows the identification accuracy with different SNR values. The result demonstrates that SNR affects the recognition accuracy of the IP responses. When the SNR is greater than 80 dB, the recognition accuracies of the three sample sets are about 90%, while when the SNR decreases to 60 dB, the identification accuracy of non-IP samples reduces to 51.35%. At the same time, the accuracy of IP samples also nearly reduced by 10%, and total accuracy is 67.35%. Therefore, SNR achieving 80 dB is one of the necessary conditions to identify the IP effects.

C. INFLUENCES OF OTHER FACTORS
When applying the TEM method in the detection of the polarizable medium, researchers found that the transmitter frequency, anti-aliasing filter, and logarithmic sampling can also affect the sign reversals and tried to determine a way to eliminate these interferences [35], [36]. Therefore, we need to consider these factors before carrying out experiments.
In our experiment, the Debye loop lies in the center of the transmitter and receiver coils. To measure a 2D or 3D target body, we need to consider the spatial relationship between the target body and the transmitter and receiver coils. Researchers found that when the target body deviates from the receiver coil, the measured data may reveal fast decays or sign reversals that may be identified as the IP responses. In that case, we can exclude them using local geological information or the same loop mode [36]. Our current attempts to identify the IP effects have achieved good results on the macroscale layered structure, and the influence of local 2D or 3D objects on the recognition accuracy is also one of the problems that we need to solve.

VI. CONCLUSIONS
We provide an IP response identification algorithm based on the PMI-FSVM method. This method can effectively identify the weak polarization effects and establish a foundation for the high-precision interpretation of EM data. The total identification accuracy is 90.7%, whereas the identification accuracy of the IP responses is 86.5%. It is easy to identify the IP response as the non-IP response when the impact ratio is less than 30%. Conductivity distribution, SNR of the system, transmitter frequency, anti-aliasing filter, logarithmic sampling, and spatial location of the target also affect the accuracy. This method can also be applied to the identification of EM responses in other dispersive media. In addition, the simulation of IP responses in field experiments provides ideas for verifying the strategy proposed in the IP research and significantly narrows the gap between theoretical analysis and practical application.
However, we must clarify that the method is restricted to the 1D model in application because the work focused on the structure that can be regarded as layered models macroscopically. It is essential to include 2D or 3D structures to expand the application area because the IP responses may exhibit some different features.