In geophysical exploration, magnetic data inversion is generally initiated by transforming so-called ill-posed problems to optimization models constructed such that the parameters of the model explained by observation data are good descriptions of the subsurface anomaly (Mehanee et al., 1998; Liu et al., 2015; Essa & Elhussein, 2020; Mbonu et al., 2021). The process of finding a solution to an inverse problem normally begins with the supposition of a good initial model. This initial model commonly made from drilling data or reports from previous geophysical expositions is progressively smoothened in stepwise iterative runs until a subjectively suitable fit between the measured and estimated data is acquired. The smoothening is executed by making meta-heuristic forward adjustments to the parameters of the model (Gharehchopogh & Gholizadeh, 2019; Wang & Li, 2019; Hayyolalam & Pourhaji Kazem, 2020). In this research, the model parameters of interest are amplitude coefficient (K) related to the composition of the body, depth (z), location of the origin (x_o), angle of magnetization (α), and a shape factor (q) whose value is unique for each different structural shape (Figure 1)

For this study, we employ the general formula for magnetic anomaly profile -T (x_j) over simple-geometrically shaped structures meticulously extrapolated from the inductive analysis of mathematical expressions for horizontal, vertical, and total magnetic anomaly of spheres (Prakasa Rao & Subrahmanyam, 1988), horizontal cylinders (Rao et al., 1973) and thin sheets (Gay, 1963). The formula (Equation 1) is given as (Mehanee et al., 2021): T ( x j ) = K Az 2 + B ( x j − x 0 ) + C ( x j − x 0 ) 2 [ ( x j − x 0 ) 2 + z 2 ] q ,     j = 1,2,3,4 , … (1) where z is depth to the buried anomaly, K the amplitude coefficient and, α the angle of magnetization normally in the plane of the principle profile coinciding with the x-direction. This is illustrated by Prakasa Rao and Subrahmanyam, 1988) for spheres and compiled by Gay (1965) for horizontal cylinders and thin sheets (Table 1). x₀ indicates the coordinate of the center of the structure, §-the dip, I_T the effective inclination of the geomagnetic field in the vertical plane normal to the body’s strike while q is the shape factor taken as 2.5 for spheres, 2.0 for infinitely-long horizontal cylinders, and 1.0 for thin geologic sheets (Mehanee et al., 2021). FHD and SHD are the first and second derivatives of the anomaly respectively while A, B, and C constants are defined in Table 2;

The five controlling model parameters (K, z, α, xo, and q) are computationally obtained by strategically introducing procedures outlined in the succeeding section on the objective function expressed in Equation 2 (Essa & Elhussein, 2020). Puzzling of optimal parameters for the geophysical model of interest is carried out such that the misfit between the true and estimated data (calculated using the objective function - Equation 2) is minimized (Mehanee et al., 1998; Ouadfel & Taleb-Ahmed, 2016). ∑ i=1 S ( T i m - T i c ) 2 S (2) where T_i ^m and T_i ^c are respectively the magnetic anomaly from observed data and that estimated using the proposed methodology.

The MRF optimization algorithm (which is the backbone of our new methodology) is designed based on the foraging strategies of marine-living Manta rays (Figure 2). The foraging strategies include chain, cyclone, and somersault techniques (Ghosh et al., 2021). For our geophysical case, the iterating vectoral positions of the Manta-rays (search agents) represents the probable positions of physical parameters sought-for. The position of the plankton itself (also vectoral) points to the geoscientific optimization problem’s solution. For the chain foraging strategy, the agents identify and frame a head-tail formation (also referred to as foraging chain) towards a potential plankton. As a result, any plankton lost by one Manta-ray will not escape the next; thereby improving the overall rate of exploitation. In the cyclone style of foraging, the manta rays swim towards the food in design spirals such that while maintaining its foraging chain-line, each agent is concurrently advancing towards its target food. Alternatively, the solution to the magnetic problem is placed as a pivot in the summersault foraging technique; such that the agents recurrently translate about this pivot before moving into new positions. The result of this is that the new positions are not purely random but rather, circumvents of the best position found (at least up to that point). These numerical constructs for these foraging strategies as adapted from Zhao et al. (2020) and remodeled for our geophysical case are described below; Chain foraging: Equation 3 T ( x j ) i D = { T ( x j ) i D ( t ) + r . ( T ( x j ) i − 1 D ( t ) − T ( x j ) i D ( t ) ) + w . ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) )  i = 2,3 , … , N T ( x j ) i D ( t ) + r . ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) ) + w . ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) )  i = 1 (3) w=2 .rand() . | log (r) | (4) where T ( x j ) ( t ) is the ith search agent’s position at a time t. rand() –a vector randomly generated within 0 and 1; w- the coefficient of weight, N is the maximum iterations allowed; while T ( x j ) b e s t D ( t ) points to the healthiest pool of planktons representing the optimum vectoral position of the sought geomagnetic parameters. The exponentiation factor D represents the dimension of the position. For the study, D = 4 since we are dealing with four parameters Equation 4.

Cyclone foraging: T ( x j ) i D ( t + 1 ) = { T ( x j ) best D ( t ) + rand ( ) . ( T ( x j ) i − 1 D ( t ) − T ( x j ) i D ( t ) ) + µ  .   ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) )  i = 2,3 , … , N T ( x j ) best D ( t ) + rand ( ) . ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) ) + µ  .   ( T ( x j ) best D ( t ) − T ( x j ) i D ( t ) )  i = 1 (5) where µ, known as the cyclone foraging’s coefficient of weight, is generated stochastically using Equation 6; µ=2 e rand() T-t+1 T()   .sin 2πrand() (6) It is noteworthy to add that after executing Equation 5 for our structural problem, we uncomfortably observe that the optimizing parameters for our geophysical structure tend to vacillate around a certain position in a loose loop-like manner. To resolve this particularity, we force each of our agents to find new positions far from its current best position in the search space by randomly injecting a wild function q ^D _rand with the whole search taken as the reference position (Ben et al., 2021b). q ^D _rand was constructed numerically using Equation 7. q rand D =L B D +rand() . (U B D -L B D ) (7) as such, Equation 7 became; T ( x j ) i D ( t+1 ) = { q rand D ( t ) +rand() .(T ( x j ) i-1 D ( t ) -T ( x j ) i D ( t ) )+µ  .( q rand D ( t ) -T ( x j ) i D ( t ) ) i=2,3,…,N q rand D ( t ) +rand() . q rand D ( t ) -T ( x j ) i D ( t ) )+µ  .( q rand D ( t ) -T ( x j ) i D ( t ) ) i=1 (8) This strategy that the algorithm extensively performs global search; substantially improving general exploration. Equation 8.

Summersault foraging: T ( x j ) i D ( t+1 ) =T ( x j ) ( t ) +P .(rand() .T ( x j ) i D ( t ) -rand() .T ( x j ) i D ( t ) ,  i=1,2,…,N (9) where the term P - a constant known as the somersault factor determines the range at which the Manta rays somersault. Equation 9.

Generally, the chain strategy encourages individual positioning based on short-term historical positioning and current global best. The cyclone foraging strategy makes each individual update its position with respect to both its preceding neighbor and the reference position while the summersault foraging technique enforces thorough adaptive search.

For exhaustive details about the adoption of the MRF algorithm for optimization scenarios, the readers are referred to Zhao et al. (2020). In our new method, we start by initiating a random population of Manta-rays in the search domain (designed around the upper bound (UB) and lower bound (LB) of our magnetic model parameters) with each agent position in the space representing characteristic parameters of the geologic structure. In our geophysical case, the selection of bounds are done based on priori geological knowledge such as information from maps, geophysical data and drilling reports. At each iterative step, the individuals are coerced to update these positions based on the referential position and the individual out front. This reference position depends on t/T which decreases from 1/T to 1 (Xu et al., 2020; Mbonu & Ben, 2021). With t/T less than rand (), exploitation is implemented else, exploration is performed. Based on this, our MRF optimization algorithm allows the individuals to decide and when necessary, switch between chain and cyclone foraging behaviors. Then through somersault foraging, the agents adaptively update their positions within the space and with respect to the most optimal position so far found.

It is noteworthy to add that these calculations and updates explained are interactively performed step-wise until the stopping criterion is met. After each step, the quality is assessed. This is done by gauging the misfit between the measured and estimated anomalies. The misfit is calculated using the root mean technique (Abdelrahman et al., 2012; Mehanee, 2014).

Eventually, convergence is achieved and the four required model parameters are returned. The pseudo-code for the whole optimization process is shown in Figure 3.

The proposed method is subjected to controlled experiments using data generated synthetically for idealized spheres, horizontal cylinders, and thin sheet models. First, analysis is done on the synthetic anomaly constructed such that it is noise-free; then the data is consciously corrupted with different levels of noise and reanalyzed.

We now seek to assess the new procedure’s performance with anomalies from complicated and interfering subsurface structures. To simulate this scenario, we once again generate model data synthetically but, in this scenario, from multiple source bodies (with multifarious parameters) placed at proximity. For the sake of this example, we adopt a horizontal cylinder model with K = 55,000 nT m², z = 15 m, α = 25, q = 2, x_o = 30 m, and a thin sheet model with K = 75 nT m, z = 9 m, α = -15, q = 2, x_o = 120 m. The profile length is 400 m. The magnetic field anomaly generated is as shown in Figure 7.

We then apply steps consistent with the proposed methodology (as with those in the two preceding sections). However, in this case, the structures are modeled together with parameter bounds as displayed in Table 4 selected. From results (Figure 7; Table 4), the parameters can be observed as being well recuperated from and showing admirable consistency with the actual field anomaly confirming the congruency of the methodology.

In this section, we compare the optimization performance of MRF with those from other methods commonly used for geophysical inversion. In our analysis, we are concerned with the convergence rate, quality/depth of exploration, and accuracy.

For the analysis, three common inversion techniques were selected namely; PSO, SA, and GA.

The performance of the new technique is further investigated using three field cases.

In our final case study, we will be analyzing a multi-model anomaly from Hamrawein field–an Egyptian mining field situated around the western throat of the Red Sea in North Africa. This geologic region is a constituted part of the East-African rift system believed to have been structurally architected from the anti-clockwise rotational divergence of the Arabian plate from the African tectonic Plate with the central Mediterranean Sea as the pole of rotation (Noweir & Fheel, 2015). The field is petrologically dominated with volcanic rocks such as pillow tholeiitic basalts overly genetically on layers of ultramafic and gabbroic strata. These layers of strata are unevenly covered with volcanic and sedimentary rocks of calco-alkaline geochemistry (Salem et al., 2005; Smith & Salem, 2005). Figure 12 is a geological map of Qusier area showing the Hamrawein field.

Figures 13, 14 show a 15 km magnetic anomaly profile (MN) extracted from a magnetic map designed from an extensive high definition aeromagnetic survey data originally carried out by Salem et al. (1999). For this study, the profile is digitized at a 200 m interval. From priori geologic studies, we can ascertain that the profile is characterized by two prominent anomalies. However, the structure or characteristics of these anomalies cannot be distinctly defined initially and as such could be interpreted broadways as a combination of any of the simple shapes. This situation normally renders inherent non-uniqueness around solutions for the anomaly.

Fudging in the peculiarity of this problem and to resolve this technical complication, we run the algorithm for all possible dual-combinable shapes. The least RMSE error of 3.341 is achieved and accepted for model shapes with q values of 1.032 and 0.981 respectively which proximate synonymously with models of thin sheet structure. Then too, our algorithm eliminants K = 115,018.987 nT m, z = 494.74 m, α = 68.276° and x_o = 4571.033 m as optimal model parameters for the first buried structure while the parameters of the second structure are computed to be K = 61,836.497 nT m, z = 458.023 m, α = 51.383° and x_o = 14897 m (Table 10). From Figure 14, it can be observed that the fit between the observed and estimated anomalies is consistently excellent.

Table 11 compares results obtained in this present study with those obtained from similar studies reported in literature. Salem et al., 2005) modified the enhanced local wavenumber method and used it to interpret the two subsurface anomalies as thin sheet structures buried at depths of 555.7 and 441.2 m respectively. Salem et al., 2005) did heterological computation of the anomaly using analytic signal derivative but reported depths of 540 and 447 m respectively. Salem (2011) on the other hand conducted his analysis of the Hamrawein anomaly using both total gradient (TG) and local wave number (LW) methods. The LG methodology outputted depths of 486.5 and 440.4 m while the LW method out-turned 432.6 and 422.8 m as respective depths to the two anomalies. Essa & Elhussein, 2018 and Essa & Elhussein (2020) used particle swarm optimization (PSO) to carry out magnetic inversion for the area. While the former employed robust PSO for this task and reported the subsurface structures as thin sheets with depths of 623.05 and 494.14 m, the latter who emphasized the use of the optimization technique to infer second moving average residual magnetic anomalies also interpreted the structures as thin sheets but emplaced them at depths of 604 and 500 m respectively. Comparing results from this study with these previous findings, it can be generally deduced that resolutions using this new methodology stand in good agreement with those from previous reports. Figure 15 shows the convergence behavior of the investigated case studies (Pima, Parnaiba, and Hamrawein anomalies).

Evaluation of results obtained using this new methodology has revealed that the technique exhibits admirable stability in the presence of noise, remarkable flexibility especially when confronted with interfering anomalies, and great pervasiveness in the quantitative resolution of magnetic inversion problems. The consistency of the results obtained from the analysis of the field examples when compared with background results from other similar studies conducted with other methods further affirms the reliability of the new methodology. While this consistency is continuous for shape and depth parameters, Tables 6, 8, 10 show that there are significant variations in the amplitude parameter for all field cases. This large variation could be attributed to two main reasons. The first, which is the mathematical role of K in the forward model (Equation 1) has already been explained in Synthetic Examples Section. Ambiguity makes up the second reason. The K parameter—directly related to the magnetic susceptibility is a physical property that points to the type of materials making up the buried anomaly. As has been explained severally in literature (Clark, 1996; Crowther, 2003; Marques et al., 2014; Teixeira et al., 2018; César de Mello et al., 2020), magnetic susceptibility is a very non-unique rock property and could sometimes range up to factors of 10⁴. It is therefore not surprising to find authors reporting different values for K. Howasmuch, it will be rather skewed to conduct appraisal of this parameter based on these literature reports alone. The efficiency of the new method in the estimation of K would only be assessed without bias if the result were gauged against laboratory findings for cores harvested from these field sites. Unfortunately, this information is not available for any of the three field cases. Nonetheless, the new technique still exhibited comparative edge over other well-known and conventional techniques especially on the grounds of convergence rate, and quality of the anomaly parameters (depth and shape) resolved.