Solution of the Optimization Problem of Magnetotelluric Sounding in Quaternions by the Differential Evolution Method

: The article discusses the application of quaternion Fourier transforms and quaternion algebra to transform Maxwell’s equations. This makes it possible to present the problem of magne-totelluric sensing (MTS) in a more convenient form for research. Studies of the inverse MTS problem for multi-layer regions are presented using the differential evolution method, which demonstrates high convergence. For single-layer regions, a new method for solving inverse problems based on minimizing the quadratic functional using conjugate optimization methods is considered. Numerical results obtained using special Python libraries are presented, with analysis and conclusions.


Introduction
Currently, magnetotelluric sounding (MTS) is a powerful and relevant tool for the study of underground geology and a variety of geophysical phenomena.The application of this method provides valuable information about the composition and structure of the Earth, which is important for solving a variety of applied and scientific problems.The natural electromagnetic field contains fluctuations of various frequencies.As a result of a phenomenon known as the skin effect, the higher-frequency vibrations of the MT field weaken faster with increasing depth, whereas the low-frequency components of the spectrum penetrate to greater depths.Consequently, the high-frequency components of the field contain information, mainly about the surface layers of the section.With decreasing frequency, the contribution of deeper parts of the section to the observed field increases, which allows us to obtain information about the deep parts of the geoelectric section [1][2][3].
The MTS method has a number of advantages.Unlike other methods, the use of generator sets is not required, however, it allows research to be carried out to a considerable depth and is available only for studies using methods with a controlled source.
The theoretical foundations of the MTS method are based on the Maxwell equation systems, which describe the electromagnetic field.Maxwell's equations are the basis of classical electrodynamics and describe the interaction of electric and magnetic fields, as well as the propagation of electromagnetic waves in space [1].
In recent years, scientists have been increasingly interested in using quaternion Fourier transforms (QFTs) to solve partial differential equations.In [4,5], the basic properties of direct and inverse partial differential transformations using QFT were studied and proved, and the applications of QFT for solving problems of thermal conductivity and the wave Computation 2024, 12, 127 2 of 22 equation were considered.It is shown that the use of the QFT makes it possible to exclude one of the spatial or temporal variables.After the transformation, the problem is considered in a private domain, allowing us to preserve the physical nature of the equation itself.New scientific results have also been obtained to improve the modification of QFT and prove its unique properties, such as convolution and correlation.Since Hamilton's work, quaternions as a mathematical tool have been used to formulate Maxwell's equations in the theory of electromagnetism.Applying quaternions to Maxwell's equations brings them to a convenient form.
In [6,7], the use of quaternions allowed the author to present electromagnetic fields in a compact form and write down a fundamental solution.
In many works [8][9][10], the properties of integration and differentiation of quaternion functions have been proved, and tables convenient for use have been compiled.The application of QFT simplifies the system of Maxwell's equations while preserving the spectral properties of impulses in space.
The general formulation of the MTS inverse problems is to restore the distribution of the physical parameters of the environment in the Earth's crust according to the data of physical fields measured on the Earth's surface.The main purpose of this approach is to study the structure of the near-surface layer of the Earth's crust, as well as to search for mineral deposits.The inverse MTS problem makes it possible to determine a geophysical model (resistance) based on measurements of apparent resistance.Since the early 1950s, when the MTS method was proposed, methods for solving inverse problems in this area began to develop and find applications in various fields, including marine geophysics, environmental studies, analysis of the electrical structure of the earth's crust and upper mantle, geological research, and exploration of mineral deposits [11,12].
MTS inverse problems are often characterized by incorrectness, which usually leads to low quality solutions and high sensitivity to noise in the input data.The standard way to deal with these problems is to reformulate the task.Inverse problems in this field are usually nonlinear and non-unique [13][14][15][16].To solve them, geophysicists usually use linearized regularization methods [17,18].In the works [2,3,13], multi-layer models for the numerical solution of the coefficient inverse problem of MTS are considered.
When using machine learning methods, this paradigm can be changed by adding additional information.This information can either be entered directly into the algorithm, or used indirectly by taking it into account when forming a training sample.An approach based on the indirect use of additional information, based on alternative measurement methods or structural assumptions, is shown in [19].The method, which includes the direct introduction of additional information, is considered in the work [20][21][22][23].This method is associated with the integration of geophysical methods, which involves the simultaneous use of data from several geophysical methods to solve the inverse problem.
Among the widely used methods of inverse MTS problems, the following can be distinguished: the Occam inversion method, which builds a model with minimal roughness at a given degree of discrepancy, providing a stable and rapidly converging solution [24,25]; fast relaxation inversion, where horizontal derivatives are approximated based on the fields of previous iterations [26]; and nonlinear conjugate gradients that prevent the need to evaluate the complete Jacobi matrix and solve the linearized inverse problem at each iteration step [18].These methods of inverse problems, as a rule, modify the objective functional to reduce the uniqueness of the solution.
In this study, multi-layer and single-layer models are considered for determining the resistance and power of horizontal isotropic layers in the problem of MTS.For the first time, the use of quaternion Fourier transforms for solving stationary and non-stationary MTS problems described by Maxwell's equations is investigated.The application of quaternions to Maxwell's equations leads them to a convenient form for solving.Algorithms for the numerical solution of the considered models have been developed.In the proposed approaches for solving inverse problems, quadratic functionals are minimized.An adaptive algorithm of differential evolution has been developed for the multi-layer model.Numer-ous calculations of the single-layer and multi-layer MTS problem in complex variables using special Python libraries are presented.The numerical results are analyzed and the revealed patterns are presented in the form of graphs.

Materials and Methods
This study used simulations of the variation of the Earth's electric field.Magnetic variations and electric earth currents observed on the Earth's surface follow specific relationships, as they represent different manifestations of the Earth's natural electromagnetic field [27].This article examines the theoretical relationship between magnetic variations and the Earth's currents from the point of view of Maxwell's equations.The studied theoretical issues allow numerical modeling of the electrical properties of the deep layers of the Earth's crust.
The study assumes consideration of lithospheric layers consisting of an upper conductive layer (sedimentary cover), intermediate resistive layers (consolidated crust and upper mantle), and a conductive base (asthenosphere) [1].
To solve this problem, boundary conditions of the third kind are considered, since they generalize the properties of Dirichlet and Neumann conditions.When solving the problem with the presented boundary conditions, it is assumed that the medium of a single-layer region is homogeneous and isotropic.For multi-layered regions or areas with more complex geometries, other boundary conditions may be required that take into account transitions between layers and their interaction.

Transformation of Maxwell's Electromagnetic Field Equations
Consider the system of Maxwell's equations in differential form rot It is well known that Equation (1) describes the law of total current; Equation (2) describes Faraday's law of electromagnetic induction, the third equation follows on from Gauss's law, and the fourth equation means that the flux of magnetic charges does not exist.
This introduced the notation ε a = ε 0 • ε, µ a = µ 0 • µ.Substituting ( 5)-( 7) into (1)-( 4), there will be 8 scalar equations for 6 scalar unknowns of the function rot div where j ex is the density of external currents, ρ is the volume charge density, and j con is the conductivity current density.
If the rot operator is from both sides of Equation ( 9), then it turns out rot rot According to the rules of vector analysis, the following takes place: where div → E = ρ ε a from Equation (10).From (12), taking into account (10) and ( 13), the following was obtained Similarly, a hyperbolic equation for H was obtained Equations ( 14) and (15) for → E and → H are of hyperbolic type and differ only in the right hand side.Therefore, in general, it can be written as If U is represented as a harmonic process Substituting in (16), the Helmholtz equation was obtained where k = ω c µ ε + i σ ωε 0 is complex dielectric constant.Equations ( 14) and ( 16) are reduced to an equation of the form (18) using the method of complex amplitudes or direct Fourier transform.This approach is widely used to solve applied problems.

Application of QFT for a Non-Stationary Problem
Further, using QFT, Cauchy problems for hyperbolic-type equations are obtained.As indicated in Section 2.1, the system of Maxwell's equations is solved by methods of complex amplitudes and direct Fourier transform.The use of QFT makes it possible to obtain a linear differential equation of the second order, the right side of which contains the spectra of pulses.
The QFT is an innovative tool in the field of signal and image processing, enriching the arsenal of data analysis and processing methods.They are a generalization of the classical Fourier transform, which is used for signal and image processing, including three-dimensional and four-dimensional samples.
Physical signals often have a three-dimensional nature, reflecting the three-dimensional space of the environment.This highlights the need to use three-dimensional samples in the analysis of such signals.It is worth noting separately that color images, in turn, consist of three components per pixel, which is due to the trichromatic nature of human vision.In this regard, QFTs are proving to be very effective tools for processing such complex data.
It is important to note that the fourth dimension of quaternions plays an essential role, allowing the representation of signals in the frequency domain using four dimensions.This not only expands the possibilities of data analysis and processing, but also allows the use of homogeneous coordinates to represent the most general set of geometric operations in three dimensions.
QFTs play a key role in translating a signal or image from the temporal or spatial domain to the frequency domain.Thanks to this representation, we can analyze signals and images in the frequency domain, which can lead to fast and efficient information extraction.
QFTs are reversible, which means that the original signals or images can be reconstructed from their representation in the frequency domain.Moreover, this representation can be modified before inverting the transformation, which opens up possibilities for modifying signals or images, for example, suppressing, attenuating, or amplifying certain frequency components.
The QFT is also a powerful tool for analyzing and processing signals and images.Their application opens up new perspectives for working with three-dimensional and four-dimensional data, enriching our knowledge and methods in the field of digital signal and image processing.
Consider the following problem: To solve the problem ( 19)-( 21), the theory of quaternions and quaternion Fourier transforms was applied.The quaternion algebra over R, denoted as H, is an associative non-commutative four-dimensional (4D) algebra, where i, j, k are imaginary units that follow the multiplication rules established by Hamilton: In [8], the following definition of QFT is given: which is defined by the following formula: , ω 2 , ω 3 ) are frequencies according to spatial variables x = (x 1 , x 2 , x 3 ).
Using the tabular formulas of the QFT, transformation of the operator ∆ from [9], and the multiplication rules (23) for Equation (19), the following was obtained: Suppose that Introducing the notation F q {U}(ω) = y, Equation ( 23) was rewritten in the following form: where g(ω, t) = F q { f }(ω).
The general solution ( 26) is represented as First, consider the following characteristic equation for a homogeneous equation: the solution of which has the form where Given the right side of Equation ( 26) and the initial conditions, the solution is represented as The application of the inverse QFT transformation allows us to write the solution of the problem (19)- (20) in the form

Application of QFT for a Stationary Problem
The application of QFT is not limited to spatial variables only, but can also be extended to a temporal variable.The use of QFT in time translates the original problem into the frequency domain, which allows for a more in-depth investigation of its nature.
The key advantage of using QFT in the time domain is the ability to analyze the frequency characteristics of a signal or image.In the frequency domain, it becomes more understandable to identify the main components of the signal and their interactions.This allows in-depth study of spectral features and variations of electromagnetic fields, which can be important in solving problems of mathematical geophysics.
Consider the following boundary value problem in the domain: The right-sided QFT transformation used for the problem ( 26)-( 30) has the following form: where ω = (ν, ω 2 , ω 3 ), q = i+j+k √ 3 ∈ H is a pure unit quaternion for which the following equality holds: Using the QFT transformation for Equation ( 26) Hence, assuming that the electromagnetic field is homogeneous in the x and y directions, the following one-dimensional problem for the Helmholtz equation with quaternionic variables was obtained: To obtain a solution to the original problem, the inverse QFT transformation is used, which is defined as The solution of the problem ( 34)- (36) in quaternion variables strongly depends on the sign of the discriminant of the characteristic equation [28,29].Therefore, the numerical solution of MTS problems is further carried out in complex variables.The authors set a goal to present a more simplified solution methodology.The proposed solution technique can be generalized to quaternionic variables, since Python has an extension to quaternionic numbers.
The QFT can process vector data more stably because quaternions naturally include information about the rotation and orientation of vectors.Combining operations into a single expression can reduce the number of intermediate calculations, which reduces the accumulation of rounding errors.

Results and Discussion
With regard to the formulation of the problem of magnetotelluric sensing, in this section, there will be a considered schematic representation of the model section (Figure 1).The model consists of several horizontally homogeneous and isotropic layers with different properties, resistivity ρ m and power h m (where m is the layer number).
Computation 2024, 12, 127 8 of 23 solution of MTS problems is further carried out in complex variables.The authors set a goal to present a more simplified solution methodology.The proposed solution technique can be generalized to quaternionic variables, since Python has an extension to quaternionic numbers.The QFT can process vector data more stably because quaternions naturally include information about the rotation and orientation of vectors.Combining operations into a single expression can reduce the number of intermediate calculations, which reduces the accumulation of rounding errors.

Results and Discussion
With regard to the formulation of the problem of magnetotelluric sensing, in this section, there will be a considered schematic representation of the model section (Figure 1).The model consists of several horizontally homogeneous and isotropic layers with different properties, resistivity  and power ℎ (where  is the layer number).Consider a right-hand coordinate system, where the  axis is pointing downwards, and the  and  axes are located in the plane separating air and Earth (Figure 1).The earth is divided into n horizontal isotropic layers with resistances  ,  , … ,  ,  and capacities ℎ , ℎ , … , ℎ .The power of the last layer (the th layer) is considered infinite.The upper half-space ( < 0) is filled with non-conducting air with a permittivity of  .The magnetic permeability  is constant and equal to one in all areas of the model [8,19].
Suppose that a primary elliptically polarized plane  homogeneous monoharmonic electromagnetic wave  ⃗ ,  ⃗ , falling from above onto the Earth's surface, propagates along the  axis.This wave can be decomposed into two independent linearly polarized plane homogeneous waves, the  ,  wave and  ,  wave.Wave  ,  causes the ground electric field  ⃗ =   ⃗ and the magnetic field  ⃗ =   ⃗.Wave  ,  causes the electric field  ⃗ =   ⃗, and the magnetic field  ⃗ =   ⃗.By superposition of these elementary waves, a field was obtained: complex vectors  ⃗ ,  ⃗ ,which are always orthogonal.Indeed, it follows from the conditions of axial symmetry of the medium that Consider a right-hand coordinate system, where the Z axis is pointing downwards, and the X and Y axes are located in the plane separating air and Earth (Figure 1).The earth is divided into n horizontal isotropic layers with resistances ρ 1 , ρ 2 , . . ., ρ n−1 , ρ n and capacities h 1 , h 2 , . . ., h n .The power of the last layer (the nth layer) is considered infinite.The upper half-space (z < 0) is filled with non-conducting air with a permittivity of ε 0 .The magnetic permeability µ is constant and equal to one in all areas of the model [8,19].
Suppose that a primary elliptically polarized plane xy homogeneous monoharmonic electromagnetic wave H i , falling from above onto the Earth's surface, propagates along the z axis.This wave can be decomposed into two independent linearly polarized plane homogeneous waves, the E l x , H l y wave and E l y , H l x wave.Wave E m x , H m y causes the ground electric field , and the magnetic field . By superposition of these elementary waves, a field was obtained: The value is called impedance.On the Earth's surface, at z = 0, the input impedance Z n was obtained, where the n index corresponds to the number of layers in the section.

Determination of the Input Impedance Using the Recurrent Formula
According to Formula (39), it is sufficient to consider the components E x and H y .In each layer, the value E x satisfies the one-dimensional Helmholtz equation: where k m is the wave number for the mth layer, defined in such a way that Re(k m ) > 0: where ρ m is the resistivity of the mth layer, respectively, ω is the frequency.The general solution of the considering equation has the form: where A m and B m are constants depending on the cut parameters and frequency.The exponent e −k m z is a set of plane homogeneous waves propagating from top to bottom, and the exponent e +k m z is a set of plane homogeneous waves propagating from bottom to top.
For the nth layer (the underlying base), the value of B n is zero, which follows from the condition of attenuation of the electromagnetic field at z → ∞ .
To determine the H y component, our study uses the formula derived from the second Maxwell equation.From Maxwell's second equation, the following expression will be obtained: According to (39), (41), and (42), the Z impedance at the z level within the m layer can be expressed as follows: By applying the identities of N.V. Lipskaya [30], we obtain a recurrent formula that expresses the dependence of the impedance on the roof of the m-th layer Z(z m−1 ) on the impedance on the m-th layer Z(z m ), the properties of the layer (k m and h m ), and the frequency ω.It allows us to express the impedance on the roof of the m-th layer through the specified parameters.
The impedance on the surface of the lower (nth) layer is equal to: This condition follows from the law of conservation of energy and electromagnetic induction.If the tangential components of the field were not continuous, then there would be sharp changes in the electric and magnetic fields at the interface, which would lead to a violation of energy conservation and electromagnetic induction.
Using the impedance at the lower boundary of the (n − 1st) layer, its value at the previous (n − 2nd) boundary will be calculated.Then at the next (n − 3rd) boundary, and so on, until the Earth's surface will be calculated.The value of the impedance on the Earth's surface will depend only on the properties of the layers (resistances and powers) and frequency.This value is a solution to the direct MTS problem in the Tikhonov-Cagniard model.
For the convenience of programming, we will bring these formulas to a certain form.The impedance on the roof of the m-th layer can be presented in the following form [3]: here R(z m−1 ) is the reduced impedance, depending on the properties of the medium: where k m = −iωµ 0 /ρ m .At relatively high frequencies, the function depends mainly on the properties of the upper layers of the section, and at low frequencies, on the properties of deep layers.Therefore, the function R(ω) can be called the frequency response of a horizontally layered section.One of the main characteristics in the study of a horizontally layered medium is the apparent resistance, which is denoted as ρ T .In the formulation of the direct problem, it is required to find the function ρ T .The formula used to sequentially recalculate the reduced impedance from the lower boundary (where R n = 1) to the Earth's surface can be written as follows in the form of recalculation of R into ρ T : (47) Expressing the hyperbolic arctangent in terms of the natural logarithm and using the following notation: the final formula for calculating the reduced impedance will be obtained: Using the resulting formula, a program can be created to solve a direct problem in the one-dimensional case.Based on this fractional-linear ratio, an algorithm has been developed to calculate the apparent resistance of a horizontally layered medium.

Numerical Results of Solving the Direct and Inverse MTS Problem
The model problem under consideration includes three layers, each of which is characterized by the values of resistivity ρ m and power h m (where m = 1, 2, 3 is the layer number) [31].In each of these layers, the magnetic permeability µ m is equal to 1.In this model, the Earth's crust is divided into n = 3 horizontal isotropic layers, each of which has its own values of resistivity ρ 1 , ρ 2 , ρ 3 and power h 1 , h 2 .It is assumed that the power of the last layer (the third layer) is infinite.
During exploratory studies, fluctuations in the magnetic telluric field are recorded, having a period from fractions of a second to several minutes, known as short-period oscillations (SPO).Short-period oscillations occur under the influence of charged particles emitted by the sun entering near-Earth space.To visually represent the graph of apparent resistance on the time axis, √ T is often used, which allows you to obtain information about the characteristics of the geoelectric section.
The relationship between frequency and period is determined in accordance with [1]: For numerical calculation, certain period values were selected, which were determined based on the physical considerations described in [1].
For calculations, the frequency values in the repartitions of 10 −2 < ω i < 10 2 were considered.As can be seen from Figure 2a, the values of ω i are defined for 200 points and are inversely proportional to the values of the period T i (Figure 2b).emitted by the sun entering near-Earth space.To visually represent the graph of apparent resistance on the time axis, √ is often used, which allows you to obtain information about the characteristics of the geoelectric section.The relationship between frequency and period is determined in accordance with [1]: For numerical calculation, certain period values were selected, which were determined based on the physical considerations described in [1].
For calculations, the frequency values in the repartitions of 10 <  < 10 were considered.As can be seen from Figure 2a, the values of  are defined for 200 points and are inversely proportional to the values of the period  (Figure 2b).In the direct problem, it is necessary to calculate the function  , which represents the apparent resistance, using the Formula (38).The reduced impedance is determined as follows: Figure 3 shows a graph of the apparent resistance for this medium model.In the direct problem, it is necessary to calculate the function ρ T , which represents the apparent resistance, using the Formula (38).The reduced impedance is determined as follows: Figure 3 shows a graph of the apparent resistance for this medium model.The apparent resistance graph demonstrates the correspondence of the reduced curve with theoretical concepts of apparent resistance.It is noticeable that at high frequencies (small periods), the apparent resistance function depends on the upper layer, and as the period increases, we get information about the second and third layers.The apparent resistance graph demonstrates the correspondence of the reduced curve with theoretical concepts of apparent resistance.It is noticeable that at high frequencies (small periods), the apparent resistance function depends on the upper layer, and as the period increases, we get information about the second and third layers.

Mathematical Formulation of the Inverse MTS Problem
The inverse MTS problem involves a function of many variables, where some variables are fixed, such as the number of layers and the frequency of the incident wave, while other variables vary, such as the thickness of the layers and their resistivity.The purpose of the inverse problem is to find the values of layer power and resistivity that minimize the sum of the squares of the deviations between the theoretical values of the apparent resistance modulus and the measured values for a given set of frequencies.

The Method of Differential Evolution
In the field of scientific and engineering research, global minimization of target functionality is one of the key tasks.The main goal is to find the optimal values of system parameters to optimize certain system properties.Usually, the system parameters are represented as a vector.The objective function is a criterion that we strive to minimize in order to achieve an optimal solution.
There are many methods of global minimization, each of which has its advantages and limitations.Some of the most common methods include evolutionary algorithms, genetic algorithms, particle swarm methods, annealing simulation optimization methods, and more [32][33][34][35].And when the objective function is nonlinear and undifferentiable, direct search methods are preferred.The most well known of these methods are the Nelder-Mead algorithms, the Hooke-Jeeves method, genetic algorithms (GAs), and evolutionary

Mathematical Formulation of the Inverse MTS Problem
The inverse MTS problem involves a function of many variables, where some variables are fixed, such as the number of layers and the frequency of the incident wave, while other variables vary, such as the thickness of the layers and their resistivity.The purpose of the inverse problem is to find the values of layer power and resistivity that minimize the sum of the squares of the deviations between the theoretical values of the apparent resistance modulus and the measured values for a given set of frequencies.

The Method of Differential Evolution
In the field of scientific and engineering research, global minimization of target functionality is one of the key tasks.The main goal is to find the optimal values of system parameters to optimize certain system properties.Usually, the system parameters are represented as a vector.The objective function is a criterion that we strive to minimize in order to achieve an optimal solution.
There are many methods of global minimization, each of which has its advantages and limitations.Some of the most common methods include evolutionary algorithms, genetic algorithms, particle swarm methods, annealing simulation optimization methods, and more [32][33][34][35].And when the objective function is nonlinear and undifferentiable, direct search methods are preferred.The most well known of these methods are the Nelder Mead algorithms, the Hooke-Jeeves method, genetic algorithms (GAs), and evolutionary strategies (ESs).These methods are effective in solving optimization problems with a nonlinear and undifferentiable objective function.All these methods use different search strategies and adaptive mechanisms to achieve global convergence.
An important aspect of global minimization is the selection and adjustment of algorithm parameters.Different algorithms have their own parameters that can affect the effectiveness of the search.Incorrect parameter selection can lead to convergence to a local minimum or convergence that is too slow.Therefore, it is important to conduct experiments with different parameter values and analyze their impact on the minimization process.
The method of differential evolution (DE) is one of the simple but effective algorithms for global minimization of the target functional.It was proposed by Rainer Storn in 1995 and has since been widely used in various fields of science and engineering.
The main idea of the differential evolution method is to efficiently explore the parameter space using evolutionary operators.The algorithm is based on the evolution of a population of parameter vectors, where each vector represents a potential solution to the optimization problem.
One of the key advantages of the differential evolution method is its simplicity and relative independence from parameters.The algorithm does not require gradient information or analytical expressions of the objective function, which allows it to be used for a wide class of tasks.It also has good convergence to the global minimum in the case of multimodal functions.
The differential evolution (DE) method is a parallel direct search method that uses parameter vectors as a set for each generation.The set consists of new population (NP) vectors of parameters, denoted as q iG , where i = 1, 2, . . ., NP, and G is the generation number.NP remains unchanged during the minimization process.The initial population may be randomly generated, especially if nothing is known about the system.By default, a uniform probability distribution is assumed for all random solutions, unless otherwise specified.In the case of a preliminary solution, the initial population is often generated by adding normally distributed random deviations to the nominal solution q nom,0 .
The main idea of DE is a scheme for generating trial vectors of parameters.DE creates new parameter vectors by adding a weighted vector of the difference between two elements of the population to the third element.If the resulting vector gives a lower value of the objective function than the corresponding element of the population, then it replaces this element in the next generation.The comparison vector may be, but is not necessarily, part of the trial vector generation process described above.In addition, the best parameter vector q best,G is estimated for each generation to track progress in the minimization process.
Using distance and direction information from the aggregate to generate random deviations leads to the creation of an adaptive circuit with excellent convergence properties.

The mutation is created by donor vectors
i f rand i,j ≤ Cr or j = j rand , . Stop condition: J(q G+1 ) < ε [26].

Numerical Results of the Inverse Problem by the Method of Differential Evolution
In this section, we present solutions to the multi-layer nonlinear inverse problem of magnetotelluric sensing by reducing it to minimizing the functional (40) according to the approach proposed in Section 4.2.This model problem was also considered in [2,3].The numerical calculations below demonstrate that our method is significantly superior to others in terms of accuracy and speed of execution.
Figure 4 shows an image demonstrating that the proposed algorithm exhibits a high rate of convergence.The value of the functional decreases to the level of 10 −6 , in just 100 iterations.We chose this accuracy for comparison with the results of the work [3].From the data in Table 1, it can be seen that the value of the functional is rapidly decreasing.
the data in Table 1, it can be seen that the value of the functional is rapidly decreasing.The differential evolution method provides high convergence in solving the inverse MTS problem due to its ability to globally optimize complex and multimodal functions.The convergence criteria used, such as changing the values of the objective function and changing the population, ensure that the solutions found are close to the global optimum, which is especially important when working with multi-layer models in MTS [36].
The differential evolution (DE) method is highly resistant to interference in the input data when solving the inverse problem of magnetotelluric sensing (MTS).This stability is  The differential evolution method provides high convergence in solving the inverse MTS problem due to its ability to globally optimize complex and multimodal functions.The convergence criteria used, such as changing the values of the objective function and changing the population, ensure that the solutions found are close to the global optimum, which is especially important when working with multi-layer models in MTS [36].
The differential evolution (DE) method is highly resistant to interference in the input data when solving the inverse problem of magnetotelluric sensing (MTS).This stability is related to its stochastic nature and the ability to optimize globally.The DE uses a population of potential solutions, which allows the algorithm to be less sensitive to noise in the data [36].The mutation and crossover stage helps to create a variety of solutions that can smooth out the effects of interference.Random components in the process of mutation and crossover provide exploration of various regions of the solution space, which helps to avoid the traps of local minima caused by noise.The selection process at each step ensures that only the best solutions that minimize the target function are preserved, which reduces the impact of abnormal data.
By considering the expression (58) and multiplying by an arbitrary function ψ(z, ω), by integrating resulting expression by z and ω after that, the following expression will appear: Given the boundary conditions ( 59) and (60), the following expressions will be obtained: By definition, the main part of the functional increment is the gradient and then the formulation of the conjugate problem and the gradient of the functional will be as follows: 5.2.Algorithm for Solving the Inverse Problem with the Landweber Method 1.The initial approximation is σ (0) ; 2. Solution of the direct problem (50)-( 52) according to σ (n) ; 3. Calculation of J σ (n) according to the Formula (54); 4. Solution of the conjugate problem (61)-( 63) according to u 0, ω; σ (n) ; 5. Calculation of the gradient of the functional J ′ σ (n) by the Formula (64); 6. Finding the next approximation σ (n+1) = σ (n) − α • J ′ σ (n) and switching to 2.

Numerical Experiments
Computational where ω 0 = 1.13 • 10 8 [13].The inverse MTS problem is numerically solved by the Landweber's method.It is necessary to restore the value of the conductivity coefficient of the The initial approximation  = 1 is given.Figures 11 and 12 show graphs of the of the inverse problem () based on the restored values of the medium cond and Figures 13 and 14 show graphs of the exact solution, the () function.As can be seen from Table 2, 276 iterations were performed, while the value functional decreased by 30.562283 times, and the norm of the difference between th and approximate values of the medium conductivity coefficient takes the 5.5 × 10 .In the Landweber method, the use of minimization of the quadratic function conjugate optimization methods makes it possible to achieve high accuracy of so due to the effective use of gradient information and faster convergence to a minimu Landweber method is known for its ability to efficiently search for solutions in the eter space, which makes it possible to more accurately reconstruct the distribution rameters in single-layer models.The use of the quadratic functional simplifies c tions, as it allows linear approximation of changes in the objective function, which r computational costs.

Conclusions
The magnetotelluric sounding method is one of the classical methods of el research aimed at studying the alternating electromagnetic field caused by m spheric and ionospheric activity in order to obtain information about the structur As can be seen from Table 2, 276 iterations were performed, while the value of the functional decreased by 30.562283 times, and the norm of the difference between the exact and approximate values of the medium conductivity coefficient takes the value 5.5 × 10 −3 .In the Landweber method, the use of minimization of the quadratic functional with conjugate optimization methods makes it possible to achieve high accuracy of solutions due to the effective use of gradient information and faster convergence to a minimum.The Landweber method is known for its ability to efficiently search for solutions in the parameter space, which makes it possible to more accurately reconstruct the distribution of parameters in single-layer models.The use of the quadratic functional simplifies calculations, as it allows linear approximation of changes in the objective function, which reduces computational costs.

Conclusions
The magnetotelluric sounding method is one of the classical methods of electrical research aimed at studying the alternating electromagnetic field caused by magnetospheric and ionospheric activity in order to obtain information about the structure of the upper layers of the Earth.The theoretical foundations of the MTS method are based on systems of Maxwell's equations that describe the electromagnetic field.Maxwell's equations in classical form are four differential equations describing the laws of conservation of charge and electromagnetic induction.This article studied the relationship between magnetic variations and Earth currents based on the system of Maxwell's equations.In the framework of the study, the Earth was considered to consist of horizontal layers with different electrical resistance, while the resistance inside each layer is considered constant, and the resistance of the upper half-space (air) is infinite.
It was proposed to solve stationary and non-stationary MTS problems using the quaternionic Fourier transform.The use of QFT makes it possible to exclude one of the spatial and temporal variables, leading to a more convenient form of solution.The application of quaternion Fourier transforms in time translates the original problem into the frequency domain, which allows for a more in-depth investigation of its nature.In the frequency domain, it becomes more clear to identify the main components of the signal and their interactions, which allows a deeper analysis of the physical aspects of the task.
The numerical solution of the direct problem of magnetotelluric sounding has shown that at relatively high frequencies the function significantly depends on the properties of the upper layers of the section and at low frequencies, on deeper horizons.Conclusions from the results of the numerical solution of the inverse problem associated with a multi-layer structure have confirmed that the differential evolution method is a very effective tool for global minimization of the objective function.
The Landweber method, based on minimizing the quadratic function using conjugate optimization methods, has been applied to solve inverse problems in single-layer MTS domains.For the case of a single-layer medium, the boundary conditions used provide all possible options for modeling the process.This makes it possible to model the field in single-layer areas with a realistic interpretation of surface and deep boundaries and their applicability to real geophysical scenarios.The Landweber method demonstrates improved accuracy and high computational efficiency.These advantages make it the preferred choice for solving complex MTS tasks, providing more accurate and efficient restoration of the parameters of geophysical models.

E
is the electric field vector, unit [V/m], → H is the vector of the magnetic field, unit [A/m], → D is the vector of electric induction, unit [C/m 2 ], → B is the vector of magnetic induction, unit [Wb/m 2 = T], → j is the vector of density of electric current, unit [A/m 2 ], -conductivity current density, and → j ex is the density of external currents.ρ is the volume density of the external electric charge, [C/m 2 ].

Figure 1 .
Figure 1.Geoelectric model of a horizontally layered medium used in MTS.

Figure 1 .
Figure 1.Geoelectric model of a horizontally layered medium used in MTS.
(a)  frequency graph (b)  period graph

Figure 2 .
Figure 2. Frequency (a) and period (b) graphs for the points taken for the computational experiment.

Figure 2 .
Figure 2. Frequency (a) and period (b) graphs for the points taken for the computational experiment.

Figure 4 .
Figure 4. Comparison of graphs of exact and approximate solutions of apparent resistance.

Figure 4 .
Figure 4. Comparison of graphs of exact and approximate solutions of apparent resistance.

Figures 5 -
7 show graphs of the real parts of the complex function u(z) for solving a direct problem, and Figures 8-10 show imaginary parts.As can be seen from the graphs, with increasing frequency, the difference in the values of the function decreases.The imaginary part of the function u(z) decreases monotonously in all cases. 2 ⋅  ; 5 ⋅  ; 7 ⋅  ; 10 ⋅  where  = 1.13 ⋅ 10 [13].Figures5-7show graphs of the real parts of the complex function () for s direct problem, and Figures 8-10 show imaginary parts.As can be seen from the with increasing frequency, the difference in the values of the function decreases.aginary part of the function () decreases monotonously in all cases.

Figure 10 .
Figure 10.Graph of imaginary part of the solution Im(U(z)) for ω = {5 • ω 0 ; 7 • ω 0 ; 10 • ω 0 }.The inverse MTS problem is numerically solved by the Landweber's iterative method.It is necessary to restore the value of the conductivity coefficient of the medium.The initial approximation σ 0 = 1 is given.Figures11 and 12show graphs of the solution of the inverse problem u(z) based on the restored values of the medium conductivity, and Figures13 and 14show graphs of the exact solution, the u(z) function.

Figure 11 .
Figure 11.Graph of the real part Re(U n (z)).
method.It is necessary to restore the value of the conductivity coefficient of the The initial approximation  = 1 is given.Figures11 and 12show graphs of the of the inverse problem () based on the restored values of the medium cond and Figures13 and 14show graphs of the exact solution, the () function.

Figure 13 .
Figure 13.Graph of the exact solution Re(U ex (z)).

Figure 14 .
Figure 14.Graph of the exact solution Im(U ex (z)).
which are always orthogonal.Indeed, it follows from the conditions of axial symmetry of the medium that → E, → H,

Table 1 .
Numerical results of solving the inverse problem.

Table 1 .
Numerical results of solving the inverse problem.

Table 2 .
Results of numerical calculations of the inverse problem.

Table 2 .
Results of numerical calculations of the inverse problem.