Electromagnetic methods are among the most important tools used in environmental and engineering studies [1, 2]. Electromagnetic inverse scattering is concerned with determining the properties of a scatterer from knowledge about scattered fields. To address the electromagnetic inverse scattering problem, the electrical properties of a scattering medium can be determined from information about the electromagnetic source and scattered data [3]. Notably, if the scatterer is a polarizable body, it can store charge when an electric field is applied to it by an electromagnetic source. This phenomenon is called induced polarization (IP) [4]. Since IP is an electromagnetic phenomenon, the presence of polarizable materials in the investigated media is an important indicator of the existence of metallic minerals [5]. These minerals may have different polarization responses that cause different IP signatures, which could be of interest to engineering and environmental applications [6–8]. However, extracting IP information from electromagnetic data is a challenging task due to the nonlinear relationship between scattered electromagnetic data and the IP properties of the scatterer [9]. Therefore, a nonlinear inversion methodology that ensures the accurate recovery of IP-related parameters from scattered electromagnetic fields is necessary. One of the requirements of inversion methodology is an efficient forward algorithm for conducting electromagnetic measurements resulted from scatterer parameters [10, 11].

Therefore, this paper is divided into two parts—the first part considers a forward formulation for modeling scattered electromagnetic fields in the presence of a polarizable medium, while the second part discusses an inversion setting that employs the neural networks approach to reconstruct IP-related parameters from scattered electromagnetic data. Several neural network design factors are studied to ensure that the best performance is achieved from an error perspective. Furthermore, the capabilities of the proposed methodology are demonstrated against different noise levels, and the results are reported.

The forward algorithm is an essential part of an inverse setting, since it (or at least a part of it) can be frequently used in inversion methodology [12]. The accuracy of inversion is significantly influenced by the accurate formulation of the forward problem. Notably, the several methods that have been developed for electromagnetic forward modeling, such as finite difference and finite element techniques, mainly rely on numerically solving Maxwell’s differential equations [13, 14]. However, these conventional numerical techniques suffer from multiple drawbacks, such as numerical dispersion and numerical polarization [15, 16]. To overcome these shortcomings, an efficient forward solver is proposed in this paper for use in an optimal setting as part of the inversion methodology. In this context, it is worth mentioning that although some approaches for the forward modeling of such an electromagnetic problem have been proposed by previous studies [17, 18], they are limited in their handling of more generalized situations characterized by alternating electrical properties per layer. In contrast, the framework proposed in this paper offers a more general solver that can handle situations in which the electrical properties of the scatterer are space-dependent on the z and/or y directions per layer. In other words, the derived formulations can be utilized for the forward modeling of both homogeneous and inhomogeneous multilayer scatterers, leading to improvements in the accuracy of extracting IP information concerning electromagnetic fields in the presence of complex structures.

The purpose of the proposed forward solver is to attain the scattered electromagnetic field arising from a vertically-layered polarizable scatterer. Here, an electric line source that is infinite and unphased in the x direction is considered the electromagnetic source that results in an electric field in the x direction. Notably, Maxwell’s equations in free space can be defined as [18]:

where J(y,z) represents the electric polarization currents and M(y,z) refers to the magnetic polarization currents. In the presence of a scattering object in free space, the equations can be reformulated as equivalent to J(y,z) and M(y,z), as follows [19]:

Furthermore, in the presence of a polarizable medium, the polarization response can be described using a constitutive formulation that connects the electric field to the electric current in the form of a frequency-dependent complex conductivity. This frequency dependence can be expressed using the Cole–Cole model [20]:

where m is a unitless number commonly referred to as charge-ability, representing the magnitude of the polarization effect. Furthermore, τ is a characteristic time constant, c denotes the frequency dependence, and σ_inf represents the conductivity in the absence of IP effects [21].

Notably, in this paper, the analysis is limited to a vertically stratified layered medium, i.e., vertical variations of the electrical properties of the medium, as shown in Fig. 1. Furthermore, since ε of most earth materials are in the order of free space permittivity, i.e., ε₀ = 8.85 × 10⁻¹² F/m [22], it is assumed that ε for all instances is equal to ε₀. It is also assumed that μ everywhere is equal to μ₀ and μ_r(z)=1, as a result of which M(y,z)=0. Furthermore, the electric field inside the scattering media is expressed as [19]:

where E_xs(y,z) stands for the total field inside the scattering object, while (y, z) and (y′, z′) are the Cartesian coordinates of an observation point and the polarization current, respectively. Furthermore, the double integral here is on the scattering object geometry, where T represents the total thickness of the scattering object along the z-direction (Fig. 1). Therefore, the incident electric field Ex0(y,z) can be defined as [19]:

where k0=wμ0ɛ0 is the free space wavenumber and H0(1)(.) is the Hankel function of the zeroth order and first kind, expressed as [19]:

where Δ and Δ_s represent the observation point and source point locations, respectively. The line source is located at Y_s = Z_s =0. By implementing the Fourier transform on y, the following integral equation was obtained:

where J_x (k_y,z)=0 at T<z<0 representing an integral equation in the z direction for each k_y where the varying layered properties are incorporated into the integral equation using:

where k22(z)=w2μ0ɛ0ɛr(z)+iwμ0σ(z).

At this point, Maxwell’s equations within the scattering object were modified into a volume integral equation based on the electric polarization currents [23]. This modification points to an intrinsic advantage of the formulations derived in this work—the modeling domain encloses only the scattering object that needs to be computed, leading to the localization of the unknown field quantities within the boundaries of the scattering object only. Therefore, the proposed method is less demanding compared to conventional numerical techniques that require the modeling domain to enclose a considerably larger region than the actual scattering object to be able to describe the inhomogeneous nature of the scatterer [24]. Effectively, the approach proposed in this paper avoids unnecessary computations of a high number of unknown field quantities spread across a wider region, which could be problematic, especially if the goal is to attain a large number of forward responses for neural network training.

The methodology used in this study to deal with the proposed volume integral equation incorporated the use of an eigensolution approach [25]. This led to the formulation of a complete, discrete orthonormal set { γ_n(k_y,z)}, which is deduced from an integral eigenvalue equation, with each element having its respective eigenvalue (β_n(k_y)). The eigenfunctions were then applied to the integral equation (9). At this point, the orthonormality of the complete set was taken into account to obtain the expansion coefficients, as follows:

This formed a linear system solution, expressed as follows:

where D(k_y) = Diag(β_n(k_y))F(k_y) and r(k_y) = ∑_nβ_n(k_y)γ_n(k_y,0)F_mn(k_y), n,m = 1,2,...N, where I is a unit matrix of size N, and N represents the number of eigenfunctions included in the representation. After computing the {J_n(k_y)} coefficients, the scattered electric field’s distribution was obtained using the following equation:

In this paper, the proposed forward solver is included in the inversion approach to retrieve the Cole–Cole parameters of the polarizable medium from scattered electromagnetic fields. Notably, inversion of electromagnetic data usually requires repeatedly evaluating a number of forward model responses through an iterative process [26, 27]. On the other hand, common numerical forward formulation approaches suffer from a computationally expensive feature—any variation in any of the medium parameters leads to unnecessary recomputation of the entire forward scheme, beginning from the first step [28, 29]. In this work, an avoidance of this drawback has been achieved to improve the forward model efficiency. To achieve this, the proposed forward formulation is implemented to be dividable into two parts. The first part deals with the computation of the complete orthonormal set for the polarization current. Notably, the calculations in this part were conducted independently of the characteristics of the scattering medium and computed only once. The second part deals with the computation of the fields for a specific choice of medium parameters. Notably, this is the only part that is recomputed based on each forward model response. This computational arrangement helped avoid the requirement of conducting many forward modeling calculations, thus offering an efficient scheme for providing numerous instances of scattered fields for various values of the medium parameters, which is a necessary step in the proposed inversion approach.

To assess the performance of the proposed forward model, its results were compared with those obtained for the frequency domain finite element analysis approach using the COMSOL software package. The scattered electric field values were compared at 10 wavenumbers, the results of which are illustrated in Fig. 2. The triangles and circles in the figure represent the scattered electric field values of the proposed forward model approach and the finite element approach, respectively.

As shown in Fig. 2, the values obtained using the proposed approach are close to those of the finite element analysis solution. In this context, the relative error term per wavenumber can be defined as follows:

where EPv is the computed scattered electric field obtained at v^th wavenumber using the proposed forward model, while EFv is the computed scattered electric field obtained at v^th wavenumber using finite element analysis. To validate the comparison results in terms of the relative error between the two approaches, Φrelv was calculated for 100 iterations of each wave-number. In addition, the average error term Φavgv was also computed, where:

The diamonds in Fig. 2 denote the results of Φavgv per wavenumber, indicating that the average error is very small at all studied wavenumbers, maintaining an error margin that does not exceed 10%.

The objective of the inversion approach employed in this paper was to invert the backscattered electromagnetic measurements to recover the physical parameters of the polarizable vertically stratified layered medium. The physical parameters of interest include conductivity at infinite frequency σ_∞, chargeability m, time constant τ, and frequency dependence c. In electromagnetic methods, the relationship between the data and the medium parameters is nonlinear [30–32]. Therefore, in the course of the inversion process in this research, medium parameter estimation was conducted by adopting a nonlinear inversion scheme based on a neural network approach.

A two-layer structure was constructed as the forward model for the scattering object. The forward model is displayed in Fig. 1, with the source and the receiver located at the same position above the scattering object. Furthermore, without compromising generality, it was assumed that the structure thickness is fixed, with both Layer 1 and Layer 2 being 100 m thick. Notably, the values of the electrical properties of each layer were mainly derived from the ranges of the electrical properties of geological configurations [33].

The input of the designed neural network was the scattered electromagnetic fields collected at nine spatial harmonics, represented by real and imaginary parts at those spatial harmonics. This resulted in an input vector of 18 values. Meanwhile, the output of the designed neural network was the corresponding IP properties of the scattering object, represented as Cole–Cole parameters. Notably, since the analysis conducted in this work was limited to a two-layer structure with known layer thicknesses, an output vector of eight values was attained (Fig. 3).

After the input and output patterns were identified, the neural network was developed through a learning process that involved making iterative changes to the scatterer’s properties. This was conducted by introducing a proper dataset to the designed neural network, comprising various input vectors as well as with their accompanying output vectors [34]. Notably, the required dataset for the learning process was obtained using the forward solver presented in this paper.

For this purpose, 25^2S forward responses were considered to represent alternating combinations of IP parameters. Furthermore, to guarantee fair representation, the box-car random number generator Rand (1, 1) was employed to explore the model-to-data space more robustly. The required dataset was then developed as follows:

with s=1:S, where S represents the number of layers of the structure. Hence,

Furthermore, while computing the forward responses for the developed dataset, the Cole–Cole parameters were chosen to vary within certain ranges coinciding with their possible values in nature. Typically, m can take values between 0 and 1, while c can change from 0.1, representing poorly sorted chargeable sources, to 1, representing ice [35]. Furthermore, σ_∞ ranges from 10⁻⁴ S/m to 1 S/m, while 10⁻⁶ s ≤ τ ≤ 5×10⁻² s [36]. Therefore, for the neural network to achieve the required performance, several design factors need to be tuned to achieve the proper network structure for addressing the inverse problem. The performance of the network with regard to each design factor was evaluated by measuring the upper and lower error bounds for each of the predicted Cole–Cole parameters, as follows:

where R stands for one of the estimated Cole–Cole parameters, f represents the parameter estimated by the neural network, v refers to the parameter values on which the model was trained, e_ub denotes the maximum overestimation of the parameter, and e_lb stands for the maximum underestimation of the parameter.

To achieve the best possible performance, the structure of the neural network had to be optimized by tuning the different design parameters. The first design parameter considered for tuning was the number of hidden layers and neurons in each layer of the neural network structure. Notably, various techniques have previously been proposed to find the optimal size of the hidden layer as a function of the input and output vector sizes or in relation to the number of training examples [37, 38]. However, all of these techniques are problem-dependent. Moreover, defining an optimal number of neurons per hidden layer continues to be an open subject [39, 40]. In this work, a hidden layer with 20 neurons was selected to model the nonlinear inverse problem based on a compromise between achieving good generalization results and avoiding probability of a linear estimate (Fig. 3).

The second design parameter that needed to be considered was the size of the different datasets used in the various processes during neural network implementation. This study involved three datasets whose sizes had to be defined—the training, testing, and validation datasets. Notably, training data are utilized during the training process to update the network’s weights, testing data are utilized in the learning process to examine the network’s reaction to untrained data, and validation data are used to inspect the network’s accuracy after the completion of the learning process [41]. Different approaches have been presented to determine the size of each dataset or to relate it to a number of weights in the network. However, since these relationships are simply rules of thumb derived from experience, they cannot be commonly used for all problems [42]. In this work, the examples generated by the forward model were partitioned into the three required datasets. The examples for each dataset were chosen randomly from the pool of generated data while accounting for the compromise between satisfying training needs and ensuring better generalization of the network. In this study, the training dataset comprised 40% of the generated examples, while the testing and validation datasets each consisted of 30% of the generated examples, which differed from those chosen for the training dataset.

The third design parameter that had to be chosen was the neural network’s training algorithm, which is an integral part of neural network implementation because it defines the methodology for updating network weights to ensure correct mapping between inputs and outputs [43]. In this study, five training algorithms representing different methodologies were implemented and compared in terms of error: the gradient descent with momentum algorithm “Tgdm”; the conjugate gradient with Powell-Beale restarts “Tcgb”; the Broyden, Fletcher, Goldfarb, and Shanno quasi-Newton algorithm “Tbfg”; the one-step secant algorithm “Toss”; and the Bayesian regularization algorithm “Trbr” [44–46].

Figs. 4 and 5 illustrate the neural network performance observed on implementing the different training algorithms, with the other network design factors kept fixed. Neural network performance is represented as the upper and lower error bounds for each of the predicted Cole–Cole parameters of each layer. For each estimated parameter, it is observed that utilizing the Bayesian regularization training algorithm “Trbr” to predict the Cole–Cole parameters result in smaller values of the upper and lower error bounds compared to the other training algorithms.

The fourth design parameter that had to be selected is the neural network’s learning function. Generally, the performance of a learning function is difficult to predict before training starts because it is affected by the problem being solved [47].

In this work, five learning functions representing different approaches were implemented and compared in terms of error: the gradient descent weight and bias learning function “Lgd,” the conscience bias learning function “Lcon,” the outstar weight learning function “Los,” the Widrow-Hoff weight/bias learning function “Lwh,” and the Hebb with decay weight learning function “Lhd” [48–50].

Figs. 6 and 7 illustrate the neural network performance for the different learning functions, with the other network design factors kept fixed. The performance was evaluated in terms of the upper and lower error bounds for each of the predicted Cole-Cole parameters of each layer. For each estimated parameter, it is observed that utilizing the Widrow-Hoff weight/bias learning function “Lwh” to predict the Cole–Cole parameters results in relatively small values for the upper and lower error bounds compared to the other learning functions.

The fifth design parameter that was optimized is the learning rate. Generally, the choice of a learning rate value involves a tradeoff between search instability and slow convergence [51]. In other words, adjustment of the learning rate depends on the complexity of the problem, and a trial-and-error procedure is normally preferred to determine a suitable value [52].

The effect of different learning rate values on neural network performance is presented in Figs. 8 and 9 for the upper and lower error bounds of each of the predicted Cole–Cole parameters for each layer. It is noted that for each predicted Cole–Cole parameter, implementing a learning rate of 0.01 results in relatively small values for both error bounds.

Based on the design presented in Fig. 3, the proposed neural network was trained and tested using the training, testing, and validation datasets, maintaining the previously selected sizes for each dataset. To avoid network bias in favor of specific output data values, the data subsets’ percentages, which were randomly chosen from the generated dataset and exhibited a box-car distribution, were kept similar to the original dataset. Furthermore, drawing on the results of the design factor analysis described in the previous section, the Bayesian regularization training algorithm “Trbr” was employed to train the proposed neural network, while the Widrow-Hoff weight/bias learning function “Lwh” was implemented during the learning process. The learning rate was set to 0.01. Fig. 10 illustrates the performance of the proposed neural network for the selected design parameters in terms of the overall mean square error (MSE). It is observed that the mean square errors of the training and testing subsets decrease gradually until they drop to relatively small values. This indicates that the chosen design parameters for the neural network resulted in a proper learning process that improved both the trained data and the generalization of the testing data.

Furthermore, considering the chosen design parameters, the neural network’s robustness in the presence of noise was evaluated by introducing random noise levels to the network input according to a box-car distribution, such that different percentages of noise were added to the training and test subsets. For performance evaluation, upper and lower error margin terms were introduced for each of the estimated Cole–Cole parameters to obtain the average upper and lower error bounds at each noise level, as follows:

Figs. 11 and 12 illustrate the neural network performance for the average upper and lower error bounds at different levels of noise, showing that both increase steadily but slowly as the percentage of noise increases. These results highlight the efficient interpolation ability of the designed neural network, resulting in average error bounds of around 10% only when the percentage of noise reaches 25%.

In this study, an inversion approach capable of recovering the Cole–Cole parameters related to the IP phenomena of a vertically stratified, polarizable, and layered medium from scattered electromagnetic data was developed. The paper maintains a two-fold structure—the first part deals with forward model computations of electromagnetic fields, while the second part demonstrates the inversion technique adopted to extract IP parameters from the fields.

In the first part, a forward semianalytical formulation for computing scattered electromagnetic fields is presented. This process involved the formation of a complete orthonormal set, whose orthogonality was used to recast a volume integral equation into a linear system for calculating scattered electromagnetic fields in the presence of a polarizable medium. To facilitate the inversion process, forward calculations were carried out in such a way that the computation of the complete orthonormal set was independent of the characteristics of the polarizable medium, which were later used separately to produce the fields for specific choices of medium parameters. This arrangement prevented unnecessary recomputing of the entire forward scheme when used as part of the inversion process, in which numerous instances of scattered fields needed to be calculated for various values of IP parameters to generate the required datasets for the proposed neural network.

In the second part, an artificial neural network is designed to perform the inversion process. The neural network received the scattered electromagnetic fields as input, and produced Cole–Cole parameters as the IP properties of the polarizable scatterer. Several neural network design factors were tuned, and comparisons were made between different training algorithms, learning functions, and learning rates in terms of error bounds. Furthermore, the performance of the proposed neural network was tested using selected design parameters. It was observed that the resulting mean square error declined significantly to reach very small values within an acceptable number of iterations. Moreover, the behavior of the proposed neural network in response to noise was studied in terms of average error bounds. It was found that the neural network trained on the selected design parameters possessed good interpolation ability, achieving accurate results with relatively small error margins in the presence of considerable noise.