Based on the Improved Shuffled Frog Leaping Algorithm for ERT Inversion

1. Introduction

In the field of geophysical electrical exploration, high-precision resistivity data offers vital insights for the in-depth analysis of underground geological structures and the characteristics of geological bodies. A precise understanding of the resistivity distribution across various underground formations aids in comprehending the occurrence and behavior of groundwater, the distribution of aquifers, and the flow paths of groundwater. This, in turn, provides scientific support for the rational development and conservation of water resources [1,2,3,4,5]; in the field of engineering geology, it can effectively identify potential geological hazards, such as faults and karst caves [6,7,8,9,10,11], ensuring the safety and stability of construction projects. It is evident that high-precision resistivity data is the cornerstone of refined detection in electrical exploration, and its importance cannot be overstated.

Direct Current Electrical Resistivity Tomography (ERT) is widely popular in near-surface exploration due to its cost-effectiveness and strong resistance to interference. However, with the development of exploration demands toward greater precision, its shortcomings in inversion accuracy have gradually become more apparent [12,13,14,15,16,17]. ERT inversion is essentially a nonlinear problem, and traditional linear inversion methods, represented by the LS, linearize it. While this simplifies the computational process to some extent, it severely compromises the accuracy of the inversion results. Nonlinear inversion methods, represented by Genetic Algorithm (GA) and Particle Swarm Optimization (PSO) [18,19,20,21,22,23], although theoretically capable of handling nonlinear problems, often encounter the issue of getting trapped in local optima in practical applications. In addition, deep learning methods have made rapid progress in recent years, but they mainly suffer from the problem of poor generalization [24,25,26,27,28]. Given the complex and variable geological conditions, the significant differences in geological conditions across different regions often lead to the existing deep learning models being unable to adapt well to these variations, resulting in a significant reduction in the reliability and accuracy of the inversion results.

In recent years, iterative inversion methods based on random search or heuristic random search have been introduced into ERT inversion. However, extensive inversion studies using real data have also revealed their shortcomings, such as slow convergence speed and the tendency to get stuck in local extrema [29,30,31,32,33]. In light of this, some scholars have proposed a novel heuristic optimization algorithm based on the foraging behavior of frogs — the SFLA [34,35,36,37,38,39]. This algorithm combines the advantages of two population-based intelligent optimization algorithms, MA (Memetic Algorithm) and PSO (Particle Swarm Optimization). It is characterized by a simple concept, few adjustable parameters, fast computational speed, strong global search capability, good robustness, and ease of implementation. It has achieved excellent results in pattern recognition, function optimization, and signal processing, and has also seen preliminary applications in the field of geophysical inversion [40,41,42,43,44].

The overall goal of this paper is to address the shortcomings of the SFLA, particularly the decline in population diversity in the later stages of the algorithm and its tendency to get stuck in local optima. This is achieved by introducing an equilibrium grouping strategy during the local search phase, which balances the contribution weights of each subgroup to the global optimal solution and suppresses local optimum traps caused by high-quality groups dominating. Additionally, to overcome the issue of the algorithm’s single-step movement rule, which leads to local optima in later stages, an adaptive adjustment operator for the movement step is designed. This dynamically adjusts the search step size, thereby improving the overall performance of the algorithm. Furthermore, the improved SFLA is applied to real-world engineering cases to verify its reliability and guiding value.

2. Inversion of ERT Based on Shuffled Frog Leaping Algorithm

In the inverse framework of the SFLA, the fitness function directly represents the evaluation criterion of the solution's quality, which essentially is the objective function to be optimized (typically the root means square error (RMS) or residual norm). The algorithm's search process is fitness-driven: during each iteration, individual position updates are guided by both the subgroup-best and global-best solutions, whose selection is entirely determined by fitness ranking (solutions with optimal goodness-of-fit dictate the search direction). Specifically, the fitness function or inversion objective function ${\mathsf{\Phi }}_i$ is constructed as:

$${\Phi }_i=\sqrt{{\sum }_{i=1}^N ({\rho }_i^{obs}-{\rho }_i^{cal}{)}^2}, i=1, 2,\dots ,N$$

(1a)

$$S={\Phi }_i+\lambda \parallel {{\Delta }^{2}m}_i{\parallel }_2^2$$

(1b)

In the equation, ${\rho }_i^{obs}$​ represents the observed value of the $i$-th block, ${\rho }_i^{cal}$ represents the forward calculated value of the $i$-th block, $N$ is the total number of measurement points, $\mathcal{R}\left(i\right)$ is the regularization term, and $\lambda$ is the weight coefficient (typically taken as 0.001 based on empirical values), ${‖{{\mathsf{\Delta }}^{2}m}_i‖}^2$ is the second-order smoothness constraint applied to the model parameters, and ${‖·‖}_2$ represents the L2 norm. Additionally, to avoid negative values for apparent resistivity and resistivity due to the large range of resistivity variations in the inversion process, the natural logarithm of the model resistivity is taken, i.e., $\ln \rho$, which reduces the resistivity range and improves the inversion stability.

For the division of subgroups, to maintain the diversity of frogs within each subgroup, the frogs are sorted by their fitness values from best to worst. As per Equation (2), the sorting is done such that the first frog enters the first group, the second frog enters the second group, and so on. The $m$-th frog enters the $M$-th group, and the $m+1$-th frog goes back to the first group, continuing in this manner. Thus, the total number of frogs, M×N, is divided into M subgroups, each containing N frogs.

$$M_j=\left\{{\mathsf{\Phi }}_{j+(n-1)\times M}|1\le n\le N\right\}; 1\le j\le M$$

(2)

In addition, during the renewal of frog individuals, the position of the block represented by the frog individual remains unchanged, and its resistivity value is updated. Within each sub - group, let ${\mathit{X}}_b$and ${\mathit{X}}_w$ represent the frog individuals with the best and worst fitness values in the current sub - group, respectively, and ${\mathit{X}}_g$ represents the frog individual with the best fitness value in the entire population. During evolution, the frog individual with the worst fitness, ${\mathit{X}}_w$, jumps towards the frog individual with the best fitness in the sub - group, ${\mathit{X}}_b$.

$${∆}_w\left(t+1\right)=rand(0,1)\times \left[\left({\mathit{X}}_{\mathit{b}}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{w}}\left(\mathit{t}\right)\right)\right]$$

(3a)

$$X_w\left(t+1\right)={\mathit{X}}_{\mathit{w}}\left(\mathit{t}\right)+{∆}_w\left(t+1\right)$$

(3b)

$$R_{min}\le {∆}_w\left(t\right)\le R_{max}$$

(3c)

In the above equation, $rand\left(0, 1\right)$ represents a random number between 0 and 1, ${∆}_w\left(t\right)$ represents the jump step (random weight), and ${∆}_w\left(t\right)\in [R_{min},R_{max}]$, where $R_{max}$is the maximum allowed step size. In the above procedure, if the fitness of the new solution is better than that of the old solution, the new solution replaces the old solution. If no improvement is made to ${\mathit{X}}_{\mathit{w}}$​, then the population's global best solution ${\mathit{X}}_{\mathit{g}}$​ is used to replace ${\mathit{X}}_{\mathit{b}}$​ and recompute, i.e., Equation (4) replaces Equation (3a).

$${∆}_w\left(t+1\right)=rand(0,1)\times \left[\left({\mathit{X}}_{\mathit{g}}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{w}}\left(\mathit{t}\right)\right)\right]$$

(4)

If the new solution still fails to improve ${\mathit{X}}_{\mathit{w}}$, a new solution is randomly generated to replace ${\mathit{X}}_{\mathit{w}}$. At this point, after completing the local search within each subgroup, all individuals from the subgroups are merged back into a complete population. The algorithm then performs global information exchange (global exploration) and re-sorts all individuals based on their fitness values. Additionally, the information of the optimal individual ${\mathit{X}}_{\mathit{g}}$​ in the entire population is updated. This process continues until the predefined global maximum iteration count is reached, or the objective function value no longer changes or changes very little, at which point the program terminates.

As can be seen, the SFLA dynamically divides the population into multiple meme groups (subgroups), with each group independently performing local search evolution. After the evolution, the groups are recombined. This mechanism combines the local refinement of meme algorithms with the global collaboration of swarm intelligence, effectively avoiding the limitations of a single search strategy; Within the meme group, the local optimal solution ${\mathit{X}}_{\mathit{b}}$ is used to guide the improvement of the worst solution ${\mathit{X}}_{\mathit{w}}$. If this fails, the global optimal solution ${\mathit{X}}_{\mathit{g}}$ is introduced. This approach enables multi-level optimization and reduces the risk of premature convergence to some extent.

3.1. Improvement of the Grouping Strategy

To address the above issues, this paper improves the grouping strategy based on the following balanced distribution design approach: After sorting the M×N individuals based on fitness, the first M individuals are randomly assigned to each group (one individual per group). The next (M+1) to (2×M) individuals are randomly assigned to each group (one individual per group), and the subsequent allocation follows a similar method, continuing until the ((M−1) × N+1) to M×N individuals are randomly assigned to each group (one individual per group).

3.2. Adaptive Movement Operator

As shown in Equation (3), the leap distance of the frog is determined by the product of the weight coefficient rand and the distance between the best and worst individuals (${\mathit{X}}_{\mathit{b}}$-${\mathit{X}}_{\mathit{w}}$​). By controlling the size of the weight coefficient rand, the frog leap distance can be adjusted. However, since the weight coefficient rand is a random number, the frog leap distance is updated randomly, making it difficult to fully utilize the guidance of the optimal solution. This can easily cause the algorithm to get stuck in local optima or skip the global optimum. To address this issue, this paper introduces an adaptive movement operator strategy in the research process to improve the algorithm's optimization accuracy and accelerate convergence.

The basic concept behind the adaptive movement factor is that, in the early stages of the algorithm, the frog should take relatively large steps to facilitate a broad search for the global optimum and effectively avoid getting trapped in local optima. As the number of iterations increases, the algorithm narrows the search space around the global optimum, enabling the frog to make smaller, more precise movements to fine-tune the solution. This adaptive strategy enhances the algorithm's ability to balance exploration and exploitation, improving both its optimization accuracy and convergence speed. Consequently, in the later stages of the algorithm, the movement step size is reduced, allowing for a more refined search within the current search space and ensuring that the global optimum is not overlooked. This gradual reduction in step size helps the algorithm focus on fine-tuning the solution and achieving better convergence, while still retaining the ability to explore larger areas during the earlier stages.

After multiple experiments, this paper adopts an exponential function to design the adaptive step-size movement operator. The adaptive movement operator used in the design is:

$${\theta }_i=a^{\left(sqrt\left(i/ii\right)\right)}$$

(5)

Where $i$ represents the current iteration number within the subpopulation, ii represents the current iteration number of the population, and $a$ is set to 0.5. By introducing the adaptive movement operator, equations (3a) and (4) are replaced by:

$${∆}_w\left(t+1\right)={\theta }_i\times \left[\left({\mathit{X}}_{\mathit{b}}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{w}}\left(\mathit{t}\right)\right)\right]$$

(6)

$${∆}_w\left(t+1\right)={\theta }_i\times \left[\left({\mathit{X}}_{\mathit{g}}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{w}}\left(\mathit{t}\right)\right)\right]$$

(7)

In the early stages of the algorithm, the population iteration counts ii is small, resulting in a larger search step size. This allows the algorithm to perform a large-scale search for the global optimum in the solution space, helping to avoid getting trapped in a local optimum. As the iteration count ii increases, the algorithm narrows down the search space where the global optimum is located. In the later stages of the algorithm, the step size decreases, enabling the algorithm to perform a refined search within the current range, which helps to prevent missing the global optimum.

4. Numerical Simulation

To verify the accuracy and reliability of the improved SFLA for inversion in two-dimensional resistivity tomography data, three typical electrical models are established for inversion analysis. These models include a double rectangular high-resistivity anomaly body, a normal fault, and a reverse fault.

4.1. Forward Modeling Calculation

As is well known, in direct current (DC) resistivity exploration, the forward modeling calculation of apparent resistivity simulates the potential response generated at the surface by the subsurface electrical structure using numerical methods. The Finite Difference Method (FDM), due to its simplicity in implementation and computational efficiency, is widely used to solve this problem, providing a theoretical data foundation for inversion interpretation. For comparative research, the resistivity of the anomaly body is uniformly set to 1000 Ω·m, while the resistivity of the surrounding rock is assumed to be 100 Ω·m. Figure 2 presents the apparent resistivity cross-sections obtained from the finite difference forward modeling calculations for three different models. The red areas represent high-resistivity anomaly bodies, and the distribution of these anomalies closely matches the actual model's high-resistivity anomaly distribution (indicated by the black dashed line, as used throughout). The central positions of the anomalies and the locations of the actual model anomalies are very similar, confirming the accuracy of the forward modeling results.

4.2. Inversion Results

The LS method has become the most widely used linear inversion technique in ERT due to its theoretical completeness and flexibility in regularization. Therefore, in this study, LS is selected as the comparison method for inversion. To simulate the noise encountered in field data collection, 5% random noise was added to the data obtained from the model forward calculation. These noisy data were then used as the observed data for the inversion process. To prevent negative apparent resistivity values due to large variations in resistivity, the natural logarithm of the model resistivity was applied. This transformation reduces the resistivity range and enhances the stability of the inversion process. In the inversion procedure, the parameters were set as follows: a total of 3600 frogs, divided into 60 groups, with 60 frogs in each subgroup, and a maximum of 800 iterations. The improved SFLA was then employed for data inversion. The inversion results were compared with those obtained using LS and the standard SFLA to validate the reliability of the proposed method.

In the inversion process using the improved SFLA, the resistivity value range is set from 2 to 1100 Ω·m. After applying the natural logarithm, this range is transformed to 0.693 to 7.003. The fitness value variation curve throughout the inversion process using the improved SFLA is shown in Figure 3. As depicted in the figure, following the inversion with the improved SFLA, the fitness curves of the three models initially exhibit a substantial decrease. After 20 iterations, the rate of decrease in the fitness value slows, indicating that the search step size is decreasing as iterations progress, allowing for a more refined search for the global optimum. Upon reaching 400 iterations, the fitness values of the three models show only minimal changes, suggesting that the optimization process has nearly converged to the global optimal solution.

In addition, to study the inversion performance of the SFLA, this paper uses the Root Mean Square Error (RMS) to measure the fit between the inversion model and the actual measured data. The smaller the RMS value, the smaller the difference between the predicted values of the model and the actual observed values, which provides a good measure of the inversion accuracy. The calculation formula is as follows:

$$RMS=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N {\left({\displaystyle \frac{{\rho }_i^{obs}-{\rho }_i^{cal}}{{\rho }_i^{obs}}}\right)}^2}$$

(8)

In the formula, $i$ represents the number of measurement points, $N$ is the total number of measurement points, ${\rho }_i^{obs}$ is the actual observed value at the $i$-th measurement point, and ${\rho }_i^{cal}$​ is the forward calculated value at the $i$-th measurement point.

Figure 4 compares the performance of three inversion algorithms in a typical geological model. In Model_A (uniform anomaly body model) (Figure 4a1-a3), although the LS inversion captures the general spatial distribution trend of the anomaly body, it exhibits significant shortcomings. The diffusion of the inversion interface results in blurred boundaries, and the convergence process displays considerable oscillations, ultimately yielding an RMS value of 1.92. The standard SFLA produces an edge localization error of 5 grid units for the anomaly body (RMS = 1.63). In contrast, the improved SFLA, by incorporating adaptive step-size control and a population diversity maintenance mechanism, greatly improves convergence stability. The boundary error of the anomaly body is generally confined to within 2 grid units (RMS = 0.91), and the volume effect is notably reduced.

Model_B (normal fault structure model) (Figure 4b1-b3) illustrates the differences in structural identification capabilities. The LS inversion results in excessive smoothing, making it difficult to identify the resistivity transition zone between the hanging wall and footwall of the fault (RMS = 2.18), and it completely loses the characteristic normal dip of the fault displacement. The standard SFLA inversion exhibits local distortions, with the fault plane dip angle error reaching 12° (RMS = 1.63). In contrast, the improved SFLA successfully reconstructs the fault's geometric elements, achieving a dip angle error of less than 5%, and clearly defining the boundary between the hanging wall and footwall (RMS = 0.94).

Model_C (reverse fault structure model) (Figure 4c1-c3) further demonstrates the robustness of the algorithm. The LS inversion fails to capture the fault morphology, resulting in an RMS value of 2.06. The standard SFLA only captures the high-resistance characteristics of the hanging wall (RMS = 1.93). In contrast, the improved SFLA effectively characterizes the reverse fault elements, achieving a reasonable reconstruction of the dip angle. The morphology of both the hanging wall and footwall is generally clear, with no oscillatory artifacts, and the RMS is reduced to 0.87. A comprehensive quantitative analysis (Table 1) reveals that the improved SFLA achieves an RMS < 1.0 in all three models, marking a 57.3±2.1% reduction compared to the LS method and a 43.7±1.8% reduction compared to the standard SFLA.

Comprehensive quantitative analysis (Table 1) shows that the improved SFLA achieves RMS < 1.0 in all three models, which is a 57.3±2.1% reduction compared to the LS method and a 43.7±1.8% reduction compared to the standard SFLA. Its advantage lies in the balanced grouping strategy, which balances the subgroup contribution weights. Additionally, the introduction of the adaptive step-size control mechanism effectively suppresses local optimal traps. This algorithm provides a new technical approach for the fine characterization of electrical structures in complex tectonic areas.

5. The Practical Application of the Improved Hybrid Frog Leaping Algorithm

The study area is located 5 km northeast of Dongchang Street, Biancheng Town, Jurong City, Jiangsu Province (Figure 5). The survey line E1 crosses borehole ZK2 (with a bedrock depth of 6.7m) and ZK3 (with a bedrock depth of 12.0m). The bedrock is limestone (resistivity theoretical value: 10²~10⁴ Ω·m), with overlying plain fill soil (50~500 Ω·m). The stratigraphic sequence revealed by the borehole lithological logs is shown in Figure 6 and Figure 7. In the improved hybrid frog leaping algorithm inversion, the total number of frogs is set to 6400 (divided into 80 groups, with 80 frogs per group). The algorithm is run for 1000 iterations, with a resistivity constraint range of 2 to 10⁴ Ω·m (after natural logarithmic transformation, this corresponds to a range of 0.693 to 9.210). The fitness variation curve in Figure 8 shows that: at iteration 9, the fitness value decreases from 4.237 to 0.612; at iteration 400, it further decreases to 0.277; and after iteration 733, it stabilizes and converges to 0.189. This reflects that the algorithm has both fast global search capabilities and fine local optimization abilities.

The inversion results comparison in Figure 9 shows that: although the LS inversion (Figure 9a) can identify the trend of bedrock burial depth (indicated by the black dashed line in the figure), the bedrock interface is blurred, and the root mean square error (RMS) reaches 2.920. The improved SFLA inversion (Figure 9b) shows a significant enhancement in boundary clarity (indicated by the black dashed line in the figure): At a horizontal distance of 55m (corresponding to ZK3 location), the inversion depth is 14.05m (measured 12.0m, error 17%); at a horizontal distance of 96m (corresponding to ZK2 location), the inversion depth is 6.9m (measured 6.7m, error 3.0%); the inversion resistivity of the limestone layer is 8.2×10³ Ω·m, which aligns with the expected geological range; the root mean square error is reduced to 1.053 (a decrease of 64%). The algorithm accurately captures the undulating morphology of the bedrock— in the horizontal distance range of 20-58m, a continuous downward trend is observed; between 58-78m, a rising feature appears; from 78-136m, a descending trend is observed, followed by a gradual uplift. This aligns with the depth variations revealed by the boreholes, indicating the accuracy and reliability of the inversion results.

6. Discussion and Analysis

The performance of the SFLA essentially depends on the synergy between parameter configuration and the evolutionary mechanism. The population size (P) directly determines the breadth of the solution space: when P is too small, the population struggles to cover key regions of the high-dimensional parameter space, leading to the omission of the global optimum; when P is too large, redundant individuals significantly increase computational costs without linearly improving the convergence probability. This contradiction requires determining the P value based on the nonlinearity of the inversion problem, which is generally positively correlated with the degrees of freedom of the model parameters.

The essence of the grouping strategy is the key link in balancing local development and global exploration. If the number of groups (m) is too small, it will weaken the independence of subpopulations, leading to homogenization of information between groups, causing the algorithm to degrade into a single population optimization; whereas if m is too large, the subpopulation size (n) will become too small, and there will be insufficient diversity within the group to support deep search. The core contradiction lies in the fact that standard grouping assigns individuals based on fitness in descending order, causing high-quality subgroups (with lower group numbers) to monopolize the update rights of the global optimum ${\mathit{X}}_{\mathit{g}}$. When these subgroups fall into local optima, low-quality subgroups (with higher group numbers) lose the ability to escape because their evolutionary direction is dominated by ${\mathit{X}}_{\mathit{g}}$—this is the root cause of premature convergence. This paper breaks this limitation through a balanced grouping strategy, forcing equal contribution weights to ${\mathit{X}}_{\mathit{g}}$​ across all subgroups, thereby ensuring population diversity from a mechanistic perspective.

The number of iterations within a group (Ni) regulates the 'exploitation-exploration' transition frequency. A small Ni will interrupt local deep search, forcing the algorithm to frequently perform inefficient global jumps; whereas a large Ni will cause the subpopulation to overly focus on local areas, losing the global perspective. Ideally, Ni should decay adaptively with the iteration process: allowing sufficient local search in the early stages, and gradually decreasing in later stages to accelerate global information integration.

The physical meaning of the maximum step size (Rmax) is the search radius of the solution space. The standard SFLA uses a uniform random step size factor in the range [0,1], whose constant variance characteristic contradicts the entropy reduction law of the search process. The piecewise adaptive step size designed in this paper maintains a large step size in the early iterations to promote exploration, and then exponentially decays to a fine-tuning mode in the later stages. This mechanism significantly improves the convergence speed of the algorithm.

The global number of iterations (K) needs to match the problem complexity. For high-dimensional non-convex problems like ERT inversion, K must at least satisfy: $K\ge {\displaystyle \frac{ln\left(1/\epsilon \right)}{ln\left(1-\delta \right)}}$ where ε is the convergence threshold, and δ is the improvement probability per iteration. Empirical evidence shows that when the model parameters exceed 100, K<500K<500 increases the risk of under-convergence by more than 60%.

Additionally, it is important to note that the computational efficiency of this improved algorithm still has a bottleneck: when the population size is 6400, on a conventional Intel Core i7-10700/32GB RAM setup, a single inversion takes about 10 hours. This is primarily due to the high complexity of the adaptive step size calculation and subpopulation collaborative optimization. The next phase of research will focus on three main directions: (1) designing a GPU-based parallel population evolution architecture to accelerate fitness evaluation and step size updates using CUDA cores; (2) developing a premature convergence prediction mechanism to dynamically prune ineffective iterations; and (3) integrating transfer learning to create a parameter initialization strategy, thereby reducing the number of iterations required for convergence. Through the above optimizations, the goal is to improve computational efficiency while maintaining accuracy.

7. Conclusion

The improved SFLA proposed in this paper employs a balanced grouping strategy to optimize subgroup contributions and integrates an adaptive step size control mechanism, effectively preventing local optimum traps. Experimental results under 5% noise conditions demonstrate that the enhanced algorithm outperforms both the LS method and standard SFLA in detecting anomalous body boundaries, achieving significantly lower root mean square error (RMS) and substantially reduced volumetric artifacts.

In field validation at Baota Mountain mining area, the improved SFLA's inversion results strongly align with borehole data: along the 55m line, inverted bedrock depth is 14.05m (measured: 12.0m, 17% error); at 96m, depth is 6.9m (measured: 6.7m, 3.0% error). The algorithm accurately resolves the southward-deepening bedrock trend. Compared to the LS method (RMS=2.920), it achieves 64% higher accuracy (RMS=1.053).

The improved algorithm demonstrates marked superiority in boundary resolution and noise resistance, confirming its strong suitability for fine-scale ERT inversion in complex geological settings. This further establishes its engineering value for detailed detection of complex structures, providing a practical tool for precision exploration.