Characterization of Induced Polarization Parameters from Electromagnetic Data using Evolutionary Approach

Article information

J. Electromagn. Eng. Sci. 2023;23(4):303-309
Publication date (electronic) : 2023 June 19
doi :
1Department of Electronics, Nuclear Materials Authority, Cairo, Egypt
2Advanced Industrial, Technical and Engineering Center, Cairo, Egypt
*Corresponding Author: Mohamed Elkattan (
Received 2022 August 23; Accepted 2023 February 13.


Electromagnetic methods are one of the most important tools for exploring the physical properties of scatterers. When the scattering object is polarizable, more information can be gained from electromagnetic measurements to determine the scatterer’s electrical properties. In this paper, an inversion scheme is introduced to extract induced polarization information from electromagnetic measurements. The inverse scattering problem is reformulated as an optimization problem that utilizes an efficient forward model solver to compute the scattered field. Furthermore, a simulated annealing approach is considered in the presented scheme, wherein different cooling schedules are evaluated to solve the inverse problem. The results of applying the proposed scheme to various case studies are reported.

I. Introduction

Interest among researchers in gaining information from interactions of electromagnetic waves with a scatterer has recently increased due to its importance in numerous science and engineering applications [13]. When an electrically polarizable medium is present within the scatterer, the electromagnetic (EM) response is affected by what is called “induced polarization” (IP) [4]. Induced polarization is a dynamic electromagnetic phenomenon that is directly sensitive to the electrical properties of the scatterer [5]. Extracting the induced polarization parameters from electromagnetic data provides additional information about the scatterer’s conductivity [6]. Consequently, this facilitates solutions to substantially more problems that have many applications in science and technology-related fields [79]. In this work, an inversion methodology is proposed to retrieve induced polarization parameters from a scattered electromagnetic field. A major challenge in inverting EM data for IP parameters is the requirement of an accurate forward model [10]. Therefore, one of the objectives of this paper is to propose a forward modeling algorithm to compute the scattered EM field in a manner that accounts for the induced polarization effect.

The other objective of this work is to propose a procedure to solve the inverse problem without affecting its nonlinear nature. This is accomplished through handling the inverse problem via an optimization procedure in which the proposed forward model can be used iteratively for different values of the unknown induced polarization parameters. Formulated as an optimization problem, the inverse problem can be tackled with one of the evolutionary approaches. Evolutionary approaches are global optimization algorithms that possess direct exploration capabilities of the entire parameters’ space. Consequently, they offer the advantages of enhancing results for highly nonlinear problems [11]. Different types of evolutionary approaches have been successively developed. Among them, the most popular are genetic algorithms (GA) [12] and simulated annealing (SA) techniques [13].

Simulated annealing is a global search approach that is able to approximate a global optimum of a given function using a limited number of control parameters. It is potentially used for approaches wherein the objective is to find an acceptable global optimum solution in a relatively large and multidimensional search space rather than to reach the exact global optimum solution [14]. Furthermore, due to its stochastic searching mechanism, it is less likely to be trapped in local minima [15]. In this work, the simulated annealing technique has been applied to solve the investigated inverse scattering problem. Different simulated annealing algorithms were designed within the proposed inversion methodology. To validate the efficiency of the inversion methodology, various different settings were considered, and the results were evaluated in terms of parameters’ estimation errors.

II. Forward Problem Formulation

In electromagnetic modeling (i.e. forward problem), one considers the problem of computing the electromagnetic data, given the physical properties of the scatterer [16]. In this paper, the forward problem is an integral part of the inversion methodology. Therefore, the efficient formulation of the forward problem is crucial for the accuracy of the inversion’s results. Here, we introduce semi-analytical forward calculations to compute the scattered electromagnetic field from the scatterer in the presence of the induced polarization effect. First, we assume a transmitter of an infinite unphased electric line source in the y direction to produce an electric field in the same direction. In this paper, we will consider a 2D scatterer model, as presented in Fig. 1, wherein the polarizable scatterer consists of an object (Ob) embedded inside a layered medium (LM). For this 2D problem, the presence of the scatterer can be represented as electric polarization and magnetic polarization currents [17]:

Fig. 1

The general model of a polarizable scatterer with an overall thickness of L.

(1) J(x,z)=(σ(x,z)-iωɛ0(ɛr(x,z)-1))E(x,z),
(2) M(x,z)=-iωμ0(μr(x-z)-1)H(x,z),

Here, we will limit the current analysis to a case wherein the permeability μ of the scatterer is the free space permeability μ0. Hence, for such a y-independent problem, the electric field inside the scatterer is calculated as follows [17]:

(3) Eys(x,z)=Ey0(x,z)-wμ040Ldz-dxJy(x,z)H0(1)(x,z,x,z),
(4) Ey0(x,z)=-ωμ04H0(1)1(k0r-rs),
(5) H0(1)(k0r-r)=1π-dkxe(ikx(x-x)+ikzz-z)kz,
(6) kz=k02-kx2,         given that imgkz>0

where Eys(x, z) is the total field inside the scatterer, and Ey0(x,z) is the incident electric field. H0(1)(.) is the Hankel function of the first kind and zeroth order, k0=ωμ0ɛ0 is the free space wavenumber, and r and rs are the observation and source points’ locations respectively. (x, z), and (x′, z′) are observation points and electric polarization current Cartesian coordinates, respectively. Here, the double integral is on the geometry of the scatterer with L referring to the layered medium thickness along the z-direction. We can now solve the integral equation for the electric polarization current as follows:

(7) J(x,z)=Q(x,z)2π-dkxi2kzeikxxeikzz+Q(x,z)2π--0Ldkxdxdzi2kzeikx(x-x)eikzz-zJ(x.z),

where Q(x, z) can be written as follows:

(8) Q(x,z)=iωμ0σ(x,z)+w2μ0ɛ0(ɛr(x,z)-1),

Using the eigen solution in [18] yields Eq. (7) in the following form:

(9) n-dkxeikxxan(kn)ψn(kx,z)=Q(x,z)2π-dkxi2kzeikxxeikzz+Q(x,z)2π--0Ldkxdxdzi2kzeikx(x-x)eikzz-z[Σn-dxxeikxxan(kx)ψn(kx,z)],

where an stands for the current expansion coefficient and ψn stands for the eigenfunction. The approach used here to solve this forward problem involves calculating the inverse Fourier transform for both sides and then multiplying by the complete set { ψm(kx,z)} [19]. Subsequently, by making use of the orthogonality of the complete discrete orthonormal set {ψn (kx, z)}, and its properties [20], we can calculate the current expansion coefficients in the following form:

(10) am(kx)=n-dkxF(m,n,kx,kx)λn(kx)ψn(kx,0)+n-dkxλn(kx)an(kx)F(m,n,kx,kx),


(11) F(m,n,kx,kx)=0L-dxdzψm(kx,z)e-ikxxQ(x,z)2πeikxxψn(kx,z),

Furthermore, for the investigated scatterer, we can define Q(x, z) as Q(x, z) = Q1 + Q2, where

(12) Q1=klayer2-k02,

everywhere in the layered medium except for the object and

(13) Q2=kobj2-k02,

everywhere inside the object:

(14) klayer2(x,z)=w2μ0ɛ0ɛr(layer)+iwμ0σ(layer),
(15) kobj2(x,z)=w2μ0ɛ0ɛr(obj)+iwμ0σ(obj),

In the presence of a polarizable scatterer, induced polarization effects can be described in terms of frequency-dependent complex conductivity, where the frequency dependence can be calculated using the most commonly used Cole–Cole model as follows [21]:

(16) σ(w)=σ(1-η1+(1-η)(iwτ)c),

where η is a physical parameter called “chargeability,” which quantifies the intensity of the induced polarization effects within the medium. τ is the characteristic time constant, and c is the Cole–Cole exponent term, which refers to the degree of frequency dependence. σis the conductivity at infinite frequency or in the absence of IP effects [22].

By approximating -dkxΣsδkx with S terms in equation (10), and by truncating Σn with N terms, we can reformulate the scatterer polarization current expansion coefficients in the form of a linear system:

(17) a(kx)=(I-D(kx))-1u(kx),

where a(kx) is the NS length vector, and u(kx) is the NS length vector. D(kx) is a square matrix of dimensions NS x NS. I is a unit matrix of size NS. Here, N stands for the total number of eigenfunctions to be involved in the representation, while S represents the horizontal wavenumbers. Once the coefficients {an (kx)} are computed, the scattered electric field’s distribution attained at the receiver can be written as follows:

(18) Eyr(kx,0)=iwμ0Σnan(kx)   λn(kx)ψn(kx,0).

III. Inversion Approach

The purpose of the inversion problem is to retrieve accurate values for the Cole–Cole parameters that describe the IP phenomenon, from the scattered EM field. As the EM response is a nonlinear function of the Cole–Cole model parameters [23], an estimation of these parameters should be performed using a nonlinear inverse approach. This is accomplished by rearranging the nonlinear inversion problem into a global optimization process, wherein optimization is conducted by minimizing the objective function, ζ, defined as follows:

(19) ζ=Σw=1WEmw-Ecw2,

where Emw is the electric field obtained at the wth wavenumber, and W represents the total wavenumbers to be included. Ecw is the electric field calculated for the wth wavenumber based on an estimation of the 4O induced polarization properties of the scatterer. Therefore, the objective function is the 4O dimensional function in the following form:

(20) ζ=({ηo},{σo},{τo},{co}),   o=1:O;

where ηo, σo, τo, and co are the chargeability, the conductivity at an infinite frequency, the time constant, and the Cole–Cole exponent of the oth object respectively. Here, the developed forward approach is used as an integral part of the inversion approach to iteratively minimize the multidimensional objective function through a simulated annealing approach. The simulated annealing algorithm starts with a random initial solution derived from the 4O Cole–Cole parameter’s estimates, which are the current estimates. The algorithm then starts to explore the search space to propose new estimates from the neighborhood to generate a proposed solution, wherein the distance between the proposed estimates and the current estimates in the search space are selected based on a probability distribution that depends upon a control variable called temperature [24]. The algorithm then begins to evaluate both the current and the proposed estimates using the objective function, ζ, to determine whether the proposed solution is better or worse than the current solution.

If the proposed estimates yield a solution with an objective function value lower than the objective function value of the solution obtained by the current estimates, the proposed estimates then replace the current estimates in the new iteration. On the other hand, the algorithm may opt to accept the proposed estimates as the new current estimates, even if the proposed estimates lead to an objective function value higher than the objective function value obtained from the current estimates. By allowing these types of choices that reflect an increase in the objective function, the opportunity to explore more areas in the search space can be increased, while the probability of solution being trapped in a false local minimum can be decreased. However, these choices are made according to a specific probability, which is based on an acceptance function:

(21) X=11+e(ζproposed-ζcurrentT)

The temperature T is essential to regulate the acceptance probability of a new solution. During the simulated annealing process, to determine the optimal parameters’ estimates, T is decreased systematically according to a designed cooling schedule Tc. The cooling schedule is applied to reduce the probability of accepting inferior estimates gradually throughout the algorithm iterations to determine the ideal estimate [25]. This results in a less favorable value of the objective function and hence to an optimal or a suboptimal solution to the inverse problem. In this paper, three cooling schedules are designed and applied: a fast, Boltzmann, and exponential cooling schedule. The fast cooling schedule is defined as follows:

(22) Tcf=T0r+1

The Boltzmann cooling schedule is defined as follows:

(23) Tcb=T01+ln(r+1)

Lastly, the exponential cooling schedule is defined as follows:

(24) Tce=γrT0,

where T0 is the initial temperature, γ is the cooling rate, and r = 1, 2, 3 . . . is the number of iterations.

IV. Simulation And Results

Here, the problem formulation is limited to estimate the Cole–Cole parameters for each of the objects and the layered medium, given the scattered electromagnetic fields. The scattered electromagnetic fields are collected at 11 wavenumbers. The forward model is depicted in Fig. 2, where the thickness of the layered medium is 150 m, while the thickness of the embedded object is 50 m. The electromagnetic source is located at xs=zs=0.

Fig. 2

The scatterer model in this paper.

Using the simulated annealing technique in the inversion scheme permits the direct inclusion of a priori information about the possible range for each of the Cole–Cole parameters. Generally, the chargeability η can take values between 0 and 1, while the time constant value τ is within the following range: 10−6 s ≤ τ ≤ 1 s [26]. The Cole–Cole exponent c can have values from 0.1 to c = 1, while the conductivity at infinite frequency σ is within the range from 10−4 S/m to 1 S/m [27]. In this paper, we selected 10 polarizable scatterers’ configurations that are constructed from the possible ranges of the Cole–Cole parameters. The investigated configurations are summarized in Table 1.

Induced polarization properties of the 10 investigated configurations in this paper

To study the results of the simulated annealing approach, an error term is defined as follows:

(25) Λ=|Pa-PePa|×100,

where Pa, and Pe are the actual and the estimated values of the parameter, respectively. Furthermore, to compare the three cooling schedules’ performances, an error margin is computed 100 times for each estimated Cole–Cole parameter per cooling schedule and defined as follows:

(26) em=ΣuΛ100,u=1,2,3,..,100.

The three designed cooling schedules use 100°C as the initial temperature and 0.95 as a cooling rate. Furthermore, in this paper, the simulated annealing algorithm is designed to seek an optimum error of 0% between the estimated, and the actual Cole–Cole parameter value. If the algorithm couldn’t achieve this optimum error value, still a stopping mechanism is imposed, where after 10, 000 iterations, the algorithm stops, and outputs the estimated Cole–Cole parameters that yields the nearest error value to the optimum.

Figs. 3 and 4 summarize the comparison between the three cooling schedules’ results and the investigated configurations in terms of the error margin of each of the Cole–Cole estimated parameters for the layered medium and the object, respectively. The comparisons in Figs. 3 and 4 reveal that, among the three proposed cooling schedules, the exponential cooling schedule yields em values of about 10% as the worst error value for an estimated parameter within the investigated configurations.

Fig. 3

Error margins of the 10 configurations’ layered mediums using exponential, Boltzmann, and fast cooling schedules: (a) chargeability; (b) conductivity at infinite frequency; (c) time constant and (d) exponent term.

Fig. 4

Error margins of the 10 configurations’ objects using exponential, Boltzmann, and fast cooling schedules: (a) chargeability; (b) conductivity at infinite frequency; (c) time constant and (d) the exponent term.

Conversely, fast and Boltzmann cooling schedules yield lower em values in some cases when compared to the exponential cooling schedule. However, both of fast and Boltzmann cooling schedules yield worst em values that are higher than that for the exponential cooling schedule for an estimated parameter, which could reach values > 20%. According to the results’ analysis, it can be determined that implementing the exponential cooling schedule for the proposed inversion approach yields a reasonably low error percentage for the Cole–Cole estimated parameters when compared with the other two cooling schedules.

V. Conclusion

In this paper, an inversion approach has been presented to estimate the induced polarization information from the electromagnetic fields using information about the electromagnetic source and the scattered electromagnetic data. The proposed approach involves a semi-analytical formulation to perform the forward calculations of the electromagnetic fields that are affected by a polarizable medium. Furthermore, the nonlinear inversion problem is handled using a simulated annealing technique with different cooling schedule designs.

The simulated annealing algorithms were applied to 10 different scatterer configurations to evaluate the effects of the different cooling schedules on the estimated induced polarization parameters. It has been determined that implementing an exponential cooling schedule within the simulated annealing algorithm is the proper inversion scheme to utilize for this problem, given a worst estimation error of about 10%. The obtained results in this research validate the accuracy of the proposed inversion scheme in determining reasonably proper knowledge about the scatterer’s induced polarization properties in light of the non-unique nature of such problems.


1. Khan A. S., Mukerji S. K.. Electromagnetic Fields: Theory and Applications 1st edth ed. Florida: CRC Press; 2020.
2. LoVetri J., Asefi M., Gilmore C., Jeffrey I.. Innovations in electromagnetic imaging technology. IEEE Antennas and Propagation Magazine 62:33–42. 2020;
3. Elkattan M., Kamel A.. Multiparameter optimization for electromagnetic inversion problem. Advanced Electromagnetics 6(3):94–97. 2017;
4. Viezzoli A., Manca G., Wjins C.. Causes and effects of the AIP trap in AEM data. Journal of Applied Geophysics 175:103970. 2020;
5. Luo Y., Zhang G.. Theory and Application of Spectral Induced Polarization Tulsa, OK: Society of Exploration Geophysicists; 1998.
6. Kratzer T., Macnae J.. Induced polarization in airborne EM. Geophysics 77(5):E317–E327. 2012;
7. Slater L.D., Lesmes D.. IP interpretation in environmental investigations. Geophysics 67(1):77–88. 2002;
8. Kemna A., Binleyz A., Slater L.. Crosshole IP imaging for engineering and environmental applications. Geophysics 69(1):97–107. 2004;
9. Kemna A., Binley A., Cassiani G., Niederleithinger E., Revil A., Slater L., Williams K. H., Flores Orozco A., Haegel F.-H., Hördt A., Kruschwitz S., Leroux V., Titov K., Zimmermann E.. An overview of the spectral induced polarization method for near-surface applications. Near Surface Geophysics 10:453–468. 2012;
10. Manca G., Viezzoli A.. A thorough synthetic study on IP effects in AEM data from different systems. ASEG Extended Abstracts 2018(1):1–7. 2018;
11. Coello C.A.C., Lamont G.B., Veldhuizen D.A.V.. Evolutionary Algorithms for Solving Multi-objective Problems 2nd edth ed. New York: Springer; 2007.
12. Haupt R. L., Werner D. H.. Genetic Algorithms in Electromagnetics New Jersey: John Wiley & Sons; 2007.
13. Ledesma S., Ruiz J., Garcia G.. Simulated annealing evolution. Simulated Annealing-Advances, Applications and Hybridizations Intech Publishing; 2012.
14. Fang Y., Ming H., Li M., Liu Q., Pham D. T.. Multiobjective evolutionary simulated annealing optimisation for mixed-model multi-robotic disassembly line balancing with interval processing time. International Journal of Production Research 58(3):846–862. 2019;
15. Aydin M. E., Fogarty T. C.. A Distributed Evolutionary Simulated Annealing Algorithm for Combinatorial Optimisation Problems. Journal of Heuristics 10:269–292. 2004;
16. Sheng X.Q., Song W.. Essentials of Computational Electromagnetics Singapore: John Wiley & Sons; 2012.
17. Volakis J.L., Sertel K.. Integral Equation Methods for Electromagnetics Raleigh, NC: SciTech Publishing; 2012.
18. Elkattan M., Kamel A. H.. Estimation of Electromagnetic Properties for 2D Inhomogeneous Media Using Neural Networks. Journal of Electromagnetic Engineering and Science 22(2):152–161. 2022;
19. Elkattan M., Kamel A.. Characterization of electromagnetic parameters through inversion using metaheuristic technique. Inverse Problems in Science and Engineering 29(4):567–585. 2021;
20. Elkattan M.. An efficient technique for solving inhomogeneous electromagnetic inverse scattering problems. Journal of Electromagnetic Engineering and Science 20(1):64–72. 2020;
21. Tarasov A., Titov K.. On the use of the Cole–Cole equations in spectral induced polarization. Geophysical Journal International 195:352–356. 2013;
22. Seigel H.O., Nabighian M., Parasnis D.S., Vozoff K.. The early history of the induced polarization method. Leading Edge 26:312–321. 2007;
23. Dias C. A.. Developments in a model to describe low frequency electrical polarization of rocks. Geophysics 65(2):437–451. 2000;
24. Siddique N., Adeli H.. Simulated annealing, its variants and engineering applications. International Journal on Artificial Intelligence Tools 25(6):1630001. 2016;
25. Li H., Landa-Silva D.. An adaptive evolutionary multi-objective approach based on simulated annealing. Evolutionary computation 19(4):561–595. 2011;
26. Zorin N.. Spectral induced polarization of low and moderately polarizable buried objects. Geophysics 80(5):E267–E276. 2015;
27. Viezzoli A., Manca G.. On airborne IP effects in standard AEM systems: tightening model space with data space. Exploration Geophysics 51(1):155–169. 2020;


Mohamed Elkattan

Dr. Mohamed Elkattan has been a lecturer at the Egyptian Nuclear Materials Authority since August 2013. He received his Ph.D. degree from Ain Shams University in Egypt in 2013. He received his M.Sc. degree in electronics and communications engineering from Cairo University in Egypt in 2007. His research interests include antennas and wave propagation.

Aladin Kamel

Prof. Aladin H. Kamel received a B.Sc. in electronics and communications engineering in 1975 (Distinction with Honors) and a B.Sc. in pure mathematics and theoretical physics in 1978 (Distinction with Honors), both from Ain Shams University in Cairo, Egypt. He conducted his M.Sc. studies in electromagnetic boundary value problems at Ain Shams University and in solid state physics at the American University in Cairo. He received a Ph.D. in Electrical Engineering from New York University in New York, USA in 1981. He won the best paper award from the Antennas and Propagation Society of the IEEE in 1981. He was an Assistant Professor of Electrical Engineering at New York University, a Manager of Research with the IBM Europe Science and Technology Division, and an executive manager with the Regional Information Technology and Software Engineering Center in Cairo, Egypt. His current research interests encompass analytical and numerical techniques related to the radiation, scattering, and diffraction of waves, as well as inverse scattering problems.

Article information Continued

Fig. 1

The general model of a polarizable scatterer with an overall thickness of L.

Fig. 2

The scatterer model in this paper.

Fig. 3

Error margins of the 10 configurations’ layered mediums using exponential, Boltzmann, and fast cooling schedules: (a) chargeability; (b) conductivity at infinite frequency; (c) time constant and (d) exponent term.

Fig. 4

Error margins of the 10 configurations’ objects using exponential, Boltzmann, and fast cooling schedules: (a) chargeability; (b) conductivity at infinite frequency; (c) time constant and (d) the exponent term.

Table 1

Induced polarization properties of the 10 investigated configurations in this paper

Scatterer Configuration no. η σ (S/m) τ (s) c
1 LM1 0.1 0.005 10−2 0.3
Ob1 0.3 0.01 10−3 0.5
2 LM2 0.2 0.01 0.1 0.7
Ob2 0.6 0.1 10−2 0.3
3 LM3 0.4 0.001 0.5 0.1
Ob3 0.75 0.018 0.15 0.25
4 LM4 0.1 1 0.08 0.1
Ob4 0.5 2 × 10−2 10−4 0.6
5 LM5 0.3 0.1 1 0.2
Ob5 0.95 10−3 0.005 1
6 LM6 0.5 10−4 0.01 0.35
Ob6 0.8 1 10−4 0.7
7 LM7 0.25 0.05 10−4 0.25
Ob7 1 10−2 10−3 0.9
8 LM8 0.6 10−3 10−3 0.5
Ob8 0.5 0.1 0.1 0.8
9 LM9 0.7 2 × 10−2 0.05 0.6
Ob9 0.9 10−3 0.08 0.5
10 LM10 0.8 0.5 10−6 1
Ob10 0.1 0.001 0.5 0.4