2.1. Sumudu Integral Transform
The Sumudu integral transform was proposed in [
64] as an alternative to the Laplace transform. It is defined as follows
If the original function
f has an argument
t, then the Sumudu image of this function has the same notation, but its argument is replaced by
u. It should be noted that the Sumudu image of a real-valued function is also a real-valued function. Thus, in subsequent calculations, unlike using the Laplace or Fourier transform, there is no need to resort to complex numbers. In addition, the calculation of electromagnetic signals through the Fourier transform at late times after turning on the source becomes very expensive due to the need to integrate rapidly oscillating and weakly decaying integrands [
65]. So, when finding the Sumudu image of a function, computational costs and random access memory requirements are reduced.
Important properties of the Sumudu transform include preserving the dimension of a function: measurement units of the function and its image are the same [
64]. The Sumudu transform has the following properties [
64]:
where a, b and c are some constants.
In order to avoid ambiguity, further on we use the superscript to denote the number of the vector element
, and raising a number to a power is designated as
. The Sumudu image of a power function has the form:
where
is the gamma function [
66]. Knowing the expansion of an analytic function into a power series, one can easily obtain a power series representation of its Sumudu image, and vice versa:
Here
n! is the factorial of an integer
n. The Sumudu transforms for a derivative and integral are as follows [
64]:
The convolution of two functions has the Sumudu image [
64]:
The following limit relations apply:
We also note that:
where
is the delta function.
By substituting the variable in (1), one can get the relationship between the Sumudu (
S) and Laplace (
L) transforms [
67,
68]:
where
is the Laplace transform. Therefore, the described technique for numerically calculating the inverse Sumudu transform can be applied to invert a real-valued Laplace transform. Also, knowing the interrelation between the two transforms, it is feasible to acquire the Sumudu image for functions with known Laplace images.
A disadvantage of the Sumudu transform is the lack of an explicit formula for obtaining its inverse transformation. Without resorting to the Sumudu transform properties and tables with images for some functions, the transformation can be conducted through solving the corresponding first-kind Fredholm integral equation. This is an ill-posed problem: when modeling electromagnetic responses, it requires a special regularizing operator taking into consideration the specifics of the measured signal [
62]. In consequence of the ill-conditioning of the first-kind Fredholm integral equation, getting the inverse Sumudu transform by this approach requires a substantially lower error of the corresponding Sumudu image than the admissible error in the resulting function. When the Sumudu image is a solution to a boundary-value problem derived with a numerical method, the achievement of a reasonably small error in the solution may result in tangible computing efforts. Accordingly, what is needed is to devise a computation technique for the inverse Sumudu transform, less sensitive to the noise level in the Sumudu image. This would enable reducing the demands for the allowable error of the Sumudu image and, as a result, could save computational resources.
An effective tool for tackling this problem is the application of advanced machine learning technologies, namely artificial neural networks (ANN). The innate characteristics of ANNs are the capacity for generalizing complex nonlinear dependencies, and resistance to noise in input data. ANN-based algorithms are being actively employed in multiple areas, which includes resistivity prospecting [
69,
70].
2.2. Modeling TEM Signals Through Sumudu Transform
Let us consider TEM sounding of the Earth’s interior. The source will be a current pulse in a transmitter circular loop of radius
r. The sounding result is a time sweep of the electromotive force (EMF) induced in a circular receiver loop of the same radius at a distance
d from the current source:
d is the distance between the centers of the transmitter and receiver loop, with
d > 2
r. To obtain a mathematical model that describes the sounding process, we apply Maxwell’s system of equations. As the boundary of the computational domain
, we use the surface
. The latter is so far from the transmitter loop that the electromagnetic field values on it can be assumed to be zero. The source of the electromagnetic field is a current pulse described by the Heaviside step function
:
where
is the current density in the transmitter loop. Therefore, the initial electromagnetic field values are equal to zero.
By utilizing the state equations and assuming the relative dielectric permittivity and magnetic permeability equal to unity, we exclude the electric and magnetic induction vectors from Maxwell’s system of equations. As a result, we formulate the mathematical model:
where
is the electric field strength,
is the magnetic field strength,
is the external electric current density,
is the specific electrical conductivity,
is the dielectric permittivity and
is the magnetic permeability,
is the external charge density.
Via the Sumudu time transform, we reduce the mathematical model (5)–(12) to the form:
Excluding the magnetic field strength from equations (13)–(18) and assuming
, we obtain the following boundary-value problem in regard to the Sumudu image of the electric field strength vector:
The relevant boundary-value problem in the time domain has the form:
We assume the current density in the loop to be constant. Then, the right-hand side of (19) will correspond (up to sign) to the Sumudu image of the δ-function, which is equivalent to the stepwise switching on of the current at time . Thus, the solution to the boundary-value problem (19)–(20) will be the Sumudu image of the fundamental solution to problem (21)–(24) in time. By means of this solution and the convolution operation, one may calculate the electric field strength for an external current pulse with any time dependence.
Depending on how the electrical conductivity depends on spatial coordinates, there is an appropriate method for solving the boundary-value problem (19)–(20) in partial derivatives – for instance, the vector finite element method [
71]. Eventually, we find the Sumudu image of the electric field strength; integrating it along the contour of the receiver loop, we can acquire the Sumudu image of the EMF induced in this loop. To obtain the EMF as function of time, it is necessary to carry out the inverse Sumudu transformation. Since the resulting field
is the fundamental solution, we assume the EMF to differ significantly from zero only at
. In this case, to perform the inverse Sumudu transformation with respect to the EMF image, the following integral equation is to be solved:
where
is the EMF induced in the receiver loop,
is its Sumudu image.
2.3. Numerical Inverse Sumudu Transform
Let us draw attention to the first-type Fredholm integral equation:
where
.
We assume the function
to be represented as a set of its values
(possibly distorted by some noise) at points
belonging to the segment
:
In this case, we replace equation (25) with the following system of equations:
The function
is approximated with its values
at points
:
The integrals in the system of equations (26) are approximated through the quadrature formula of the trapezoidal method with nodes
:
In this manner, from solving the integral equation (25) we move on to solving a system of linear algebraic equations by the collocation method based on quadrature formulas [
72]:
which for brevity is denoted as:
System (27) is ill-conditioned, since it is a discrete analogue of the first-kind integral equation (25), ill-conditioned one [
73]. In this case, the vector elements on the right side of
g may contain some noise of varying intensity. To solve the system (28), we use Tikhonov’s regularization method [
73] and proceed to the following minimization problem:
where
is the regularization parameter,
is the regularizing non-singular matrix. The solution to problem (29) has the form:
The specific type of the matrix depends on the properties of the kernel and on the right-hand side of integral equation (25), and, as a consequence, on the properties of the proposed solution . This choice often seems to be heuristic. One common type of the matrix is the identity matrix. The choice of a specific type of , in the context of the transient electromagnetic sounding problem being solved, is discussed below.
In order to select the value of the parameter
, we use the quasi-optimal criterion [
73,
74], for which purpose we introduce the function:
As a quasi-optimal value of
, we choose the smallest of the values
that ensure the local minimum
. According to [
73], the following relationship holds:
We present this expression as the difference of two vectors:
where
.
We introduce the matrix
:
Then the difference on the right side of (31) can be represented as follows:
Let us give consideration to the simple iteration method with the preconditioning matrix
to solve the system of equations (28). The first two iterations of the method look like:
The second term in (33) is a refinement of the approximate solution obtained at the first iteration. Consequently, the difference on the right side of (31) can be interpreted as an addition to the approximate solution in the simple iteration method with the preconditioning matrix . Based on this and the quasi-optimality criterion, we deduce the choice of the parameter α to be carried out in such a way that the preconditioning matrix provides a small correction to the first approximation in the simple iteration method; this correction might be omitted without significantly deteriorating the solution accuracy.
2.4. Features of numerical inverse Sumudu transform for TEM signals
As consistent with the initial conditions (22) and (23), at
the EMF induced in the receiver loop is equal to zero. With allowance for this and the limiting property (2), we set the elements of the vectors
and
with index 0 equal to zero, since they correspond to the values of the EMF at
and the Sumudu images of the EMF at
. Based on this, we exclude from system (27) the vector element of unknowns
and the first equation. As a result, we acquire a modified system:
Further on, by the system of equations we mean the modified system (34).
As far as the integration region
contains drastically different-scale
values, the
and
points used to discretize the integral equation (25) are specified using the following relation:
where
The
and
values are selected depending on the expected properties of the time-dependent EMF induced in the receiver loop.
As described in [
75], the absolute values of the EMF induced in the receiver loop can differ by more than 6 orders of magnitude over time. Accordingly, different elements of the
vector will differ just as much. Therefore, choosing the identity matrix as the regularizing
does not provide a satisfactory result. This is because the elements of the
vector, corresponding to the values of the sought-for function at late times after turning on the current in the source, do not have a significant effect when calculating the second term in (29) compared to those associated with the function values at early times. In order for the second term in (29) to have sufficient sensitivity to all elements of the
vector, the diagonal matrix
is to be used as the matrix
, the diagonal elements being specified in the following manner:
where
is a parameter. An increase in the parameter
leads to a sensitivity growth of the second term in (29) to the values of the sought-for function at late times. At the same time, the sensitivity to the function values at early times deteriorates. Determining the value of
requires keeping the sensitivity balance at all times.
Thus, an approximate inverse transformation for the Sumudu image of the EMF can be performed through the following expression:
To specify a pair of the parameters
and
, we draw on the idea of choosing the optimal (in some sense) preconditioning matrix in the simple iteration method. In this case, the preconditioning matrix will depend on as many as two parameters –
and
:
When inverting the Sumudu image of the EMF, the direct use of the quasi-optimal criterion for choosing the parameter
also encounters difficulties. They are associated with various scales of the elements of the vector
at different times. The elements of the difference vector (31) are highly different-scale, with the sensitivity of the vector difference norm to the late-time elements being extremely low. To overcome this deficiency, we scale the elements of the difference vector from (31). We define the following function:
where
As the optimal pair of the parameters
and
we choose
and
, which provide a minimum of the function
. We search for an approximate solution to the minimization problem by completely enumerating the values
from the set where elements are given as:
Note that the above-proposed technique for the approximate solution to the integral equation (25) can be applied for calculating the inverse Laplace transform of a real function. This function can be either known or obtained through (3) from some Sumudu image. For the inverse Laplace transform to be performed, the Sumudu kernel
is to be replaced with the Laplace kernel
, and the collocation points
– with
according to the formula:
In sum, the inverse Sumudu transformation can be conducted in two ways. The first one is to utilize the above-proposed method directly to the Sumudu image. The second is to apply this method to invert the Laplace image obtained using (3). The computational properties of these two approaches depend on the different spectral properties of the linear system matrices (34) that approximate the kernels of the Sumudu and Laplace transforms.
2.5. Computational Experiment
The use of the Sumudu transform to model TEM signals is illustrated in the following example. The modeling area is divided by a horizontal plane into two half-spaces that are homogeneous in physical properties: the upper half-space is non-conductive air; the lower half-space is a conductive earth. There are circular transmitter and receiver loops located on the Earth’s surface. The current is switched off at time . It can be shown that the electric field strength obtained under similar conditions is the fundamental solution in time to problem (21)–(25). Accordingly, its Sumudu image in time satisfies (19)–(20). If the radii of the loops are sufficiently small compared to the distance between their centers, then the transmitter loop can be replaced with a vertical magnetic dipole; the EMF in the receiver coil can be taken to be proportional to the time derivative of the z-component of the magnetic field strength at the center of the receiver loop.
For the Laplace image of the function
there exists an analytic expression [
75]:
where
,
is the specific electrical conductivity of the lower half-space,
r is the distance between the centers of the loops,
M is the magnetic dipole moment. Through the connection between the Laplace and Sumudu transforms (4), we write the Sumudu image of the function
as follows:
Having performed the analytic inverse Laplace transformation for (37), we obtain the analytic expression for
[
75]:
where
,
is the error function:
Figure 1 gives the function
and its Sumudu image (38) at
r = 100 m,
σ = 0.01 S/m. By the example of
and its Sumudu image, we examine the features of performing the inverse Sumudu transformation for the TEM sounding problem via a numerical solution of the first-type integral equation (25). The function
has characteristic properties for the EMF induced in the receiver loop (
Figure 1). At small times the function does not change noticeably, but at later times it begins to vary in proportion to
.
Let us focus on the capability of the proposed computational method for the inverse Sumudu transformation of function (38). We define a set of points
and
through (35). To do this, one is to set the values of the parameters
,
and
. As computational experiments have shown, a fairly accurate result is acquired if the following heuristic rule is employed to set the values of
and
:
, where
b is the point on the time axis for which the relation holds:
Finding the value of
b requires knowing the function
in advance. An estimate (possibly not very accurate) of this function can be made via the method proposed in the presented research, by setting the value of the parameter
b a priori large. The value of the parameter
n is further chosen to be equal to 100 in all cases. Having the set of points
and
, we generate the matrix and vector of the right side of linear system (34). To solve this system, we resort to Tikhonov’s regularization method (36). We find the parameters
и
by means of the proposed modification of the quasi-optimal criterion in the following range:
The right-side vector of linear system (28) may be distorted by noise. This may happen when obtaining the elements of this vector from the solution to problem (19)–(20), which was found by grid methods (finite difference, finite element, etc.) or some other approximate methods. To simulate this effect, we modify the elements of the right-side vector
without noise:
where
is the assigned relative noise level. The noise level of a particular element of the right-side vector depends on the value of the element itself. This is because different elements of the vector have different scales. As a consequence, additive noise will greatly distort individual elements of the vector with little effect on others.
We introduce the relative error function:
where
is the exact function,
is an approximation of the function
.
Now, move on to the influence of the noise level
in the Sumudu image (38) on the relative error of inverting the Sumudu image with the proposed method.
Figure 2 illustrates the dependence of the relative error
on time and on noise level for TEM signals at a distance between the loops of 100 m and the lower half-space conductivity of 0.1 S/m. As previously noted, the proposed method for conducting the inverse Sumudu transformation, with minor modifications, can be applied to invert a real-valued Laplace image. Let us use this fact to obtain an approximation of function (39) from the Laplace image (37).
Figure 3 shows the relative errors
for the noise level
and
in the Laplace image (37) with a distance between the loops of 100 m and conductivity of the lower half-space of 0.1 S/m.
As it appears from the figures, the relative error in calculating the inverse Sumudu and Laplace transformations through the proposed method directly depends on the noise level in the original image. The lower the noise level, the smaller the error of the reconstructed function. Moreover, if the noise level is not very high, then more accurate results are obtained by inverting the Laplace transform. However, at a sufficiently high noise level, inverting the Laplace transform does not allow achieving a satisfactory result. The high level of relative error in the reconstructed functions at early times is associated with oscillations of these functions around the exact value. In order to solve the inverse problem of TEM sounding, the EMF values are required at later times – at which the level of relative error (depending on the noise level) takes acceptable values.
Thus, to calculate the inverse Sumudu transformation by the proposed method at low noise levels in the Sumudu image, one should use the connection between the images of the Laplace and Sumudu transforms, obtain the Laplace image and invert it. In contrast, at sufficiently high noise levels it is necessary to apply the proposed method directly to inverse the Sumudu transform.
2.6. Neural Network Algorithm Development
A neural network algorithm is created with the help of supervised learning. During the first stage, we acquire a data set with training examples, each representing a pair “function ” – “Sumudu image ”. This data set has been generated drawing on mathematical modeling for characteristic geoelectric situations and sounding systems. The functions and their images have the form of value vectors at predetermined points .
The acquired data pairs are then normalized by multiplying by . Such a transformation ensures the exclusion of the sounding system parameters: current strength, loop radii, number of winding turns, etc.
The training set is extended by data augmentation: extra data are created from existing ones by way of simple transformations not requiring substantial computational power. Augmentation increases both data volume and diversity. As far as this technique facilitates recognizing data patterns, it is effective for reducing overfitting of machine learning models [
76]. Owing to the linearity property of the Sumudu transform, from two data pairs
we can get a third one as their linear combination:
where
a and
b are constants. The resulting data have to be scaled in the same manner.
In order for the stability of the algorithm to increase, we add normally distributed noise to the input data . The root-mean-square deviation of the noise holds proportion to the signal level and is equal to 5%, which matches with the expected measurement accuracy of the sounding system.
In TEM exploration, the relative change in the signal appears to be more important than its absolute value. This is why we turn to training data transformation: it distorts input and output data spaces for more efficient ANN training. To recapture the required signal dimension, post-processing is applied to the ANN calculation result.
To control the neural network learning (to deter overfitting), the resulting data set (pairs “Sumudu image ” – “function ”) is separated into two subsamples: 75% directly for training (“training data”) and 25% for control (“test data”).
The inverse Sumudu transform algorithm is accomplished on the basis of ANNs – a class of universally approximating mathematical models [
77,
78]. In a generic form, an ANN is an arbitrary function of input arguments and internal parameters fitted during the learning process:
where
are the input arguments represented by the Sumudu-image vector values;
W are the internal parameters of the ANN, which are searched while training;
is the set of functions approximating the vector of the inverse Sumudu transform;
n is the number of elements in the input and output vectors.
So, the ANN consisting of a combination of linear and nonlinear operations approximates the inverse Sumudu transform by converting the Sumudu images of the functions into their original form. At each training iteration, the neural network result is checked against the known mathematical modeling result.
The architecture of the developed ANN (
Figure 4) is a multilayer perceptron (fully connected neural network). First, there is an input layer with a vector of 100 elements (Sumudu image of the function). Also, there are
N hidden layers, each with
M neurons accompanied by the nonlinear operation ReLu (“rectified linear unit”), And, finally, the output layer where the vector is finally formed – the result of the inverse Sumudu transform (also of 100 elements). The ultimate version of the ANN architecture encompasses 4 hidden layers of 64 neurons.
The weight coefficients of neurons in the ANN layers are initially set at random fashion and then fitted in the process of training. In the present case, estimating the optimal parameters of the ANN is a supervised learning task, which is addressed with the backpropagation algorithm. The ANN training is fulfilled with the Adam algorithm – a modification of stochastic gradient descent with the adaptive first-order and second-order momentum estimates [
79]. The cost function minimized during the training is the mean absolute error (MAE). This function, coupled with the technique of preprocessing data from the training set enables minimizing the relative deviation between the ANN calculation result and true sounding curves. For the learning efficiency to improve, we practice a sequential decrease in the nominal step of gradient descent (training rate) depending on its iteration number.
The architecture of the ANN is finally selected subsequent to the results of numerical experiments. The total training time is 40 minutes (the number of epochs is 2000,
Figure 5). The algorithm is trained via parallel computing on the GPU RTX 2080 graphics accelerator.
2.7. Neural Network Algorithm: Results and Discusion
After obtaining the optimal values of the weight coefficients, the trained ANN can be used for the inverse Sumudu transformation. The neural network algorithm is tested on an additional set of data that are not involved in the ANN training (
Figure 6).
Near the zero transition point of the signal to be reconstructed (characteristic minimum in the middle of the time axis), the relative error occasionally increases as opposed to the other time intervals. This is because the absolute value of the first logarithmic derivative of the given function seems to exceed those at the other time intervals.
Table 1 gives the results of evaluating the algorithmic performance on the training and test data.
Figure 7 contains overall histograms of pointwise residuals retrieved from the test subsample data not employed during the training.
As is seen from the data analyzed, our neural network algorithm allows for the inverse Sumudu transform with high accuracy, reasonably good for tackling practical issues. It is significant that when solving first-kind integral equations by traditional methods, the characteristic error in the resulting solution is found to be significantly higher owing to the ill-conditioning [
73].
Performance assessment for the algorithm concerned is performed on the central processor CPU Intel i7-8700 (scrutinizing the time needed for the inverse transformation of one Sumudu image). The computing time of the neural network algorithm averaged over ten thousand runs is 2.9·10
-2 s, whereas that of the numerical algorithm [
62] equals to 9.3 s. The neural network calculations can also be made in batch mode, making possible efficient parallel computations on multi-processor devices.
To recap, the developed high-performance neural network algorithm exhibits a higher operation speed (320 times faster than the numerical algorithm) with lower resource intensity. In view of the attained transform accuracy, the algorithm is to be implemented into software for modeling TEM signals.