Harris Hawks Optimization for Soil Water Content Estimation in Ground-Penetrating Radar Waveform Inversion

Abstract

Ground-penetrating radar (GPR) has emerged as a promising technology for estimating the soil water content (SWC) in the vadose zone. However, most current studies focus on partial GPR data, such as travel-time or amplitude, to achieve SWC estimation. Full waveform inversion (FWI) can produce more accurate results than inversion based solely on travel-time. However, it is subject to local minima when using a local optimization algorithm. In this paper, we propose a novel and powerful GPR waveform inversion scheme based on Harris hawks optimization (HHO) algorithm. The proposed strategy is tested on synthetic data, as well as on field experimental data. To further validate our approach, the results of the HHO algorithm are also compared with those of partial swarm optimization (PSO) and grey wolf optimizer (GWO). The inversion results from both synthetic and real experimental data demonstrate that the proposed inversion scheme can efficiently invert both SWC and layer thicknesses, thus achieving very fast convergence. These findings further confirm that the HHO algorithm can be effectively applied for the quantitative interpretation of GPR data.

1. Introduction

Accurately characterizing the soil water content (SWC) in the vadose zone is essential crucial for various disciplines, including hydrology [1,2], agriculture [3,4], soil science, and environmental science. Traditional methods, like time-domain reflectometry (TDR) [5], exhibit significant limitations, such as being labor-intensive, time-consuming, and causing soil damage. Consequently, ground-penetrating radar (GPR) has emerged as a promising geophysical technique for SWC estimation due to its high sensitivity to changes in SWC [6]. It offers a non-destructive and efficient alternative to traditional methods.

Several review articles have been published to scrutinize the evolution of GPR applications in SWC estimation [6,7,8,9]. When estimating SWC, three types of GPR systems are utilized: surface GPR, borehole GPR, and off-ground GPR. Surface GPR is the most frequently employed technique, encompassing a single-offset mode [10,11,12,13,14,15,16,17,18,19] and multi-offset mode [17,18,20,21,22,23,24]. The multi-offset is an alternative when the air wave and the ground wave cannot be separated. However, the multi-offset measurement modes, including common mid-point and wide-angle reflection and refraction, are typically time-consuming and require significant measurement effort. The advancement of multi-channel GPR array antennas has helped to mitigate these drawbacks [25,26]. Surface GPR has a limited detection depth. To gain deeper information, borehole GPR is another measurement system that includes three types of measurement modes: zero-offset profile [27,28,29,30,31], multi-offset profile [32,33], and vertical radar profile [34]. Among these three modes, the vertical radar profile is less invasive and more cost-effective because it requires only one borehole. The third system is off-ground GPR, which has been developed because it is less invasive than surface and borehole GPR. Off-ground GPR is usually mounted on a vehicle or a drone, making it more suitable for rapid surface SWC estimation [35,36,37,38,39,40,41,42,43]. However, most current studies rely on GPR travel-time data to estimate the SWC, which offers limited resolution to characterize the subsurface because only some parts of the information contained in the GPR data are considered.

Full waveform inversion (FWI), which utilizes the entire waveform information, was initially proposed for seismic reflection exploration [44]. Compared with travel-time data-based inversion, FWI can provide higher resolution to characterize the subsurface. FWI has been rapidly developed for estimating subsurface properties using GPR data, as both seismic and GPR methods are founded on solving the wave equation [45,46,47]. In the context of SWC estimation, several studies have demonstrated the efficacy of GPR FWI [37,40,48]. These studies constructed an objective function by comparing Green’s function between the calculated data and observed data, subsequently estimating the SWC based on off-ground reflection data [49].

In the field of FWI, the inversion techniques are mainly classified into two main groups: those based on local optimization and those based on global optimization algorithms. Many researchers have conducted studies on GPR FWI using local optimization algorithms [50,51,52]. However, the FWI inherently poses as a highly nonlinear and multi-parameter inversion challenge. As such, the problem space is characterized by the presence of numerous extreme points. A notable issue with local optimization algorithms is their tendency to become trapped in local minima. To overcome this obstacle, one of the strategies is to obtain a reasonably accurate model derived from ray-traced inversion, which will serve as the initial model for FWI. However, consistently obtaining an initial model that is proximate to the global minima within the problem space is a challenging endeavor. To mitigate the risk of falling into local minima, the inversion based on global optimization is a superior option. For the global optimization algorithm, a variety of models are initially defined in the problem space. These models are then updated together, and they gradually move towards the global minima. Global optimization methods are often inspired by the social or biological behavior of species. The literature presents a multitude of successful implementations of global optimization. By considering a wide range of initial models and employing sophisticated search strategies, the global optimization methods can significantly reduce the chances of being trapped in local minima and enhance the robustness of GPR FWI. In global optimization algorithms, the grey wolf optimizer (GWO) was used to estimate SWC from GPR waveform data [53]. However, the GWO has limitations when dealing with high-dimensional and multi-modal functions.

The Harris hawks optimization (HHO) algorithm simulates the hunting behavior of Harris’s hawks, and it is an emerging global optimization algorithm proposed to solve optimization problems [54]. As a newcomer of population-based swarm intelligence optimization algorithms, the HHO has strong global search capability, which can efficiently explore complex search space and reduce the risk of falling into local minima. Additionally, the HHO shows a faster convergence rate on many test problems and allows for faster identification of the best solutions. It also can better maintain the diversity of the population, thus improving the comprehensiveness of the search by simulating Harris’s hawks’ hunting and strategy changes. Therefore, the HHO algorithm has begun to be applied in many areas [55,56,57].

In this paper, we demonstrate the performance of a HHO application on the inversion of GPR waveform data for the estimation of layer thicknesses and SWC in the unsaturated zone. We test the proposed scheme on data with and without noise. We also test it on field experimental data. Furthermore, in order to verify the performance of HHO, we compare its results with those of the PSO and GWO algorithms. The inversion results derived from both synthetic and experimental field data demonstrate that HHO exhibits high capacity and very fast convergence to invert SWC and layer thicknesses.

2. Materials and Methods

2.1. Harris Hawks Optimization (HHO)

The HHO algorithm was developed by Heidari et al. in 2019, which is inspired by the cooperative behavior and hunting process of Harris’s hawks in nature [54]. The evolution of the population within this algorithm is achieved by mathematical modeling of various predation strategies. There are mainly three stages in the HHO algorithm: exploration, the transition from exploration to exploitation, and exploitation. In the subsequent sections, we will mathematically model the exploratory and exploitative phases of the HHO algorithm according to the Harris’s hawks’ techniques of exploring a prey, executing a surprise pounce, and employing diverse attacking strategies. Figure 1 displays all the phases for HHO.

2.1.1. Exploration Phase

In this part, the exploration mechanism is first modeled mathematically. Harris’s hawks track and detect the prey with their sharp vision, yet at times prey detection proves challenging. Therefore, there are two strategies for Harris’s hawks to await prey detection. In HHO, the Harris’s hawks are regarded as the candidate solutions and the optimal candidate solution at each stage is considered as the targeted prey. Considering both strategies, and assuming each perching strategy has an equal likelihood, the perch based on the locations of other family members and the prey can be modeled as follows:

(1)$\overrightarrow{X}(t+1)=\left\{\begin{array}{c}{\overrightarrow{X}}_{rand}(t)-r_1\left|{\overrightarrow{X}}_{rand}(t)-2r_2\overrightarrow{X}(t)\right|,q\ge 0.5\\ ({\overrightarrow{X}}_{prey}(t)-{\overrightarrow{X}}_m(t))-r_3(LB+r_4(UB-LB)),q<0.5\end{array}\right.$

where $\overrightarrow{X}(t+1)$ represents the position vector of hawks in the next iteration $t$, ${\overrightarrow{X}}_{rand}(t)$ is a randomly selected hawk from the current population, $\overrightarrow{X}(t)$ represents the current position vector of hawks, ${\overrightarrow{X}}_{prey}(t)$ stands for the position of rabbit, $r_1$, $r_2$, $r_3$, and $r_4$ are random numbers between 0 and 1, $q$ is also random number in range of 0 and 1, which is used to decide which strategy to adopt and is updated in each iteration, and ${\overrightarrow{X}}_m(t)$ is the average position of the current population of hawks, which is calculated as follows:

(2)${\overrightarrow{X}}_m(t)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{X}_i(t),$

where $X_i(t)$ represents the position of each Harris’s hawk in iteration $t$ and $t$ indicates the total number of hawks.

The first formula in Equation (1) updates the positions of individuals based on random position ${\overrightarrow{X}}_{rand}(t)$ in the population, considering the impact of random individuals in the population on Harris’s hawks’ search for prey. The second formula in Equation (1) uses the difference between the optimal position ${\overrightarrow{X}}_{prey}(t)$ of the current population and the average position ${\overrightarrow{X}}_m(t)$ of the population, augmented by a random scaling component of a variable range $[LB,UB]$, to update the position of the individual. The $r_3$ is a scaling coefficient. In this rule, $r_3$ and $r_4$ provide additional site alternatives for population individuals, aiding Harris’s hawks in exploring different areas of the search area.

2.1.2. Transition Phase

The HHO algorithm can shift between the exploration and the exploitation according to the escaping energy of the prey. The energy of a prey decreases dynamically during the escaping process. The following equation is used to model the escaping energy of the prey:

(3)$E=2E_0(1-{\displaystyle \frac{t}{T}}),$

where $E$ indicates the escaping energy for the prey, $E_0$ represents the initial state of the escaping energy of the prey, and $T$ is the maximum number of iterations.

In each iteration, when the $E_0$ decreases from 0 to −1, it means that the prey is physically weakening, whereas when the $E_0$ increases from 0 to 1, it indicates that the prey is strengthening. During the iterations, the dynamic escaping energy $E$ presents a declining trend. When the escaping energy $\left|E\right|$ is greater than or equal to 1, the HHO is in the exploration phase, where Harris’s hawks search varied areas to explore the location of prey. Conversely, when $\left|E\right|$ is less than 1, the algorithm tries to exploit the neighboring area of the current solution to locate the best solution, implying that the HHO is in the exploitation phase.

2.1.3. Exploitation Phase

When the $\left|E\right|$ is less than 1, the HHO algorithm enters the exploitation phase. There are four potential strategies for modeling the attacking stage according to the escaping behaviors of the prey and chasing strategies of Harris’s hawks in real-world conditions. A random number that represents the chance of a prey successfully escaping is defined. When $r<0.5$, it means that the prey escapes successfully. Conversely, when $r\ge 0.5$, it indicates that the prey does not escape successfully.

1.. Soft besiege

When $0.5\le \left|E\right|<1$ and $r\ge 0.5$, Harris’s hawks will quietly encircle their prey softly to make the prey more exhausted before launching a sudden assault. This strategy is called soft besiege, which can be modeled as follows:

(4)$\overrightarrow{X}(t+1)=({\overrightarrow{X}}_{prey}(t)-\overrightarrow{X}(t))-E\left|J{\overrightarrow{X}}_{prey}(t)-\overrightarrow{X}(t)\right|,$

where $J=2(1-r_5)$ represents the random jump strength of the prey over the procedure of escaping, and $r_5$ is a random number in the range of (0, 1). In each iteration, the value of $J$ will change randomly to simulate the nature of prey motions.

2.. Hard besiege

When $\left|E\right|<0.5$ and $r\ge 0.5$, Harris’s hawks hardly encircle their intended prey to finally launch a sudden assault. This strategy is called hard besiege, which can be modeled as follows:

(5)$\overrightarrow{X}(t+1)={\overrightarrow{X}}_{prey}(t)-E\left|{\overrightarrow{X}}_{prey}(t)-\overrightarrow{X}(t)\right|,$

3.. Soft besiege with progressive rapid dives

When $0.5\le \left|E\right|<1$ and $r<0.5$, the prey possesses sufficient vitality to successfully escape and still a soft besiege is formed before launching a sudden assault. This procedure is more intelligent compared to the previous case. Here, the levy flight (LF) is introduced to mathematically model the current escape behaviors of the prey. The LF is used to emulate the real zigzag deceptive movements of the prey during the escape phase and the irregular, sudden, and rapid dives of the hawks around the prey. This strategy is called soft besiege with progressive rapid dives, which can be modeled as follows:

(6)$\overrightarrow{Y}={\overrightarrow{X}}_{prey}(t)-E\left|J{\overrightarrow{X}}_{prey}(t)-\overrightarrow{X}(t)\right|,$

Equation (6) bears resemblance to Equation (4), yet differences do exist. Subsequently, Harris’s hawks will compare the result from Equation (6) with the previous dive to investigate its feasibility as a dive (a candidate solution). If it is unreasonable (indicating Harris’s hawks consider the movements of the prey as deceptive), Harris’s hawks will not use Equation (6) to finish diving. Instead, they will perform irregular, sudden, and rapid dives as they near the prey. According to the LF, the dive will be defined as follows:

(7)$\overrightarrow{Z}=\overrightarrow{Y}+\overrightarrow{S}\times LF(D),$

where $\overrightarrow{Y}$ is calculated based on Equation (6), $D$ is the dimension of problem space, $\overrightarrow{S}$ is a random vector with dimension, and $LF$ represents the levy flight function, which can be calculated as follows:

(8)$LF(x)=0.01\times {\displaystyle \frac{\mu \times \sigma }{{\left|\upsilon \right|}^{\frac{1}{\beta }}}},$

(9)$\sigma ={\left({\displaystyle \frac{\Gamma \left(1+\beta \right)\times \sin \left(\frac{\pi \beta }{2}\right)}{\Gamma \left(\frac{1+\beta }{2}\right)\times \beta \times 2^{\left(\frac{\beta -1}{2}\right)}}}\right)}^{{\displaystyle \frac{1}{\beta }}},$

where $\mu$ and $\upsilon$ are random numbers between 0 and 1, and $\beta$ is a constant value which is usually set to be 1.5.

In other words, when Harris’s hawks discover that their prey is engaging in deceptive behavior, the action they take is adjusted according to their expected action, which is represented by Equation (6). The adjustment is achieved through the levy flight, as outlined in Equation (7). The selection of this particular diving mode is designed based on the fitness of the Harris Hawk’s current position as follows:

(10)$\overrightarrow{X}\left(t+1\right)=\left\{\begin{array}{c}\overrightarrow{Y},\hspace{1em}\hspace{1em}if\hspace{0.33em}F(\overrightarrow{Y})<F(\overrightarrow{X}(t))\\ \overrightarrow{Z},\hspace{1em}\hspace{1em}if\hspace{0.33em}F(\overrightarrow{Z})<F(\overrightarrow{X}(t))\end{array}\right.$

where $\overrightarrow{Y}$ and $\overrightarrow{Z}$ are derived from Equations (6) and (7), respectively. If we consider the optimization problem, then the $F\left(\cdot \right)$ is the objective function for that optimization problem.

4.. Hard besiege with progressive rapid dives

When $\left|E\right|<0.5$ and $r<0.5$, the prey has not got enough energy to escape and Harris’s hawks will build a hard besiege before launching a sudden assault to catch and kill the prey. Harris’s hawks endeavor to shorten the distance between their average position and the escaping prey. Equation (10) is still used as the rule for the hard besiege condition. However, the $\overrightarrow{Y}$ is calculated using the new rules as follows:

(11)$\overrightarrow{Y}={\overrightarrow{X}}_{prey}(t)-E\left|J{\overrightarrow{X}}_{prey}(t)-{\overrightarrow{X}}_m(t)\right|,$

where ${\overrightarrow{X}}_m(t)$ is calculated by Equation (2).

2.2. HHO for GPR FWI

We develop a software package for GPR FWI to estimate SWC through HHO based on MATLAB 2020a. In this study, we only focus on inversion for SWC and layer thicknesses in the vadose zone by fixing attenuation factor. There are six steps for inverting SWC for GPR waveform data, which are as follows:

Define the parameters for the algorithm (population size, maximum number of iteration, upper and lower boundaries for the search space of model);

Initialize the parameters for the models (SWC and layer thicknesses);

Calculated the soil relative dielectric permittivity according to the Topp equation [58];

Conduct the forward modeling for all the initial models of the entire population;

Calculate the fitness for every model and update the initial models by HHO algorithm;

If reaching the maximum number of iterations or meeting the stopping criterion, the iterations will stop and output the best solution. If not, a new iteration will begin from stage 3. In this study, the stopping criterion can be expressed as follows:

(12)${‖\overrightarrow{d}-\overrightarrow{G}(\overrightarrow{m})‖}_2<tol,$

where $\overrightarrow{d}$ represents the observed data, $\overrightarrow{m}$ indicates model parameters including the SWC ($\theta$) and the layer thicknesses ($z$), $\overrightarrow{G}$ is a nonlinear forward operator, and $tol$ stands for the predetermined threshold value. It means that the iteration will stop when the $L_2$ norm between the observed data and calculated data is lower than the predetermined value.

2.2.1. Objective Function

The purpose of inverse problem is to invert the model parameters based on observed data within a predetermined solution space. In general, the inverse problems are solved by an iteration algorithm. During the iteration, an objective function needs to be defined.

In this study, we first calculate the $L_2$ norm between the observed and calculated data. Then, we calculate the $L_2$ norm for the observed data. Finally, we build the objective function according to the equation as follows:

(13)$fitness=100\times {\displaystyle \frac{{‖{\overrightarrow{E}}_{obs}-{\overrightarrow{E}}_{cal}‖}_2}{{‖{\overrightarrow{E}}_{obs}‖}_2}},$

where ${\overrightarrow{E}}_{obs}$ and ${\overrightarrow{E}}_{cal}$ represent the observed data vector and the calculated data vector, respectively. Figure 2 shows the flow for inverting the SWC from GPR waveform inversion with HHO algorithm.

2.2.2. Petrophysical Relationship

In order to obtain models of electrical properties for GPR forward modeling, we need a petrophysical relationship to convert the SWC to electrical properties. In this paper, we only focus on the relative dielectric permittivity by ignoring the conductivity. Therefore, we select the widely used Topp equation to calculate the relative dielectric permittivity from the volumetric SWC [58]. With the classical parameters, the Topp equation can be expressed as follows:

(14)$\epsilon =3.03+9.30\times \theta +146.00\times {\theta }^2-76.70\times {\theta }^3,$

where $\epsilon$ stands for the relative dielectric permittivity of the soil, and $\theta$ represents the volumetric SWC. In this paper, we ignore the errors in the petrophysical relationship.

2.2.3. GPR Forward Modeling

In GPR inversion, the forward modeling is a great important link. At each iteration, the response of a subsurface model is required to be calculated. Subsequently, the parameters for the model are inverted by gradually narrowing the difference between the observed data and the calculated forward response. Maxwell’s equations are the basic theory for the modeling of electromagnetic waves. Combining Maxwell’s equations and constitutive relations, we can gain the propagation and diffusion equations of the electromagnetic waves, which are expressed as follows:

(15)$\Delta \overrightarrow{E}=\mu \sigma {\displaystyle \frac{\partial \overrightarrow{E}}{\partial t}}+\mu \epsilon {\displaystyle \frac{{\partial }^2\overrightarrow{E}}{\partial t}},$

(16)$\Delta \overrightarrow{H}=\mu \sigma {\displaystyle \frac{\partial \overrightarrow{H}}{\partial t}}+\mu \epsilon {\displaystyle \frac{{\partial }^2\overrightarrow{H}}{\partial t}},$

where $\overrightarrow{E}$ stands for the electrical field, $\overrightarrow{H}$ stands for the magnetic field, $\mu$, $\sigma$, and $\epsilon$ represent the magnetic permeability, the electrical conductivity, and the dielectric permittivity of the medium, respectively, and $t$ is the time. In this study, we only focus on the electric field. Therefore, we will not mention Equation (16) in the following section.

In this study, we assumed the electrical field to be a time-harmonic wave ($\overrightarrow{E}={\overrightarrow{E}}_0{e}^{-i\omega t}$), and so Equation (15) can be written as follows:

(17)$\Delta \overrightarrow{E}+k^2\overrightarrow{E}=0,$

where $k$ is the complex wavenumber, $\omega$ is the angular frequency. Equation (17) is also called the Helmholtz equation. The $k$ can be expressed as follows:

(18)$k^2=\mu \omega (\epsilon \omega +i\sigma ),$

When considering the attenuation ($\alpha$), Equation (18) can be written as follows [59]:

(19)$k=\omega \sqrt{\mu (\epsilon +i{\displaystyle \frac{\sigma }{\omega }})}=\beta +i\alpha ,$

where $\alpha$ and $\beta =\omega /V$ represent the absorption coefficient and phase velocity, respectively, and $V$ stands for the velocity of electromagnetic wave in the medium.

In this study, we realize the forward modeling by using the solution of Equation (17) given by Banoas follows [59]:

(20)$E(z,\omega )\exp (ikz)=E(0,\omega )\exp (i\beta z)\exp (-\alpha z),$

Let $E_0(\omega )$ represents the complex spectrum of the electrical source $e_0(t)$ at $z=0$. Assuming a reflecting layer exists at depth of $z=z_1$, the complex spectrum resulting from a vertically incident wavelet, after it has traversed a homogeneous absorbing layer from $z=0$, is given as follows:

(21)$E(\omega ,z_1)=G(z_1)R(\omega )E_0(\omega )\exp \left(i{\displaystyle \frac{\omega }{V}}2z_1\right)\exp \left(-\alpha 2z_1)\right),$

where $G\left(z_1\right)$ stands for the geometrical spreading and $R\left(\omega \right)$ is the reflection coefficient at $z=z_1$. $\alpha$ is equal to $\omega /2VQ$, where $Q$ represents the quality factor. Finally, we can apply an inverse Fourier transform to Equation (21) to yield a GPR response [59].

3. Results and Discussions

3.1. Numerical Examples

In order to test the performance of the proposed inversion scheme with HHO algorithm, we designed three synthetic examples. In Model A, the soil was assumed to be homogeneous, there was only one kind of soil, and there was a water table at 0.85 m. For Model B and Model C, we assumed that the subsurface was divided into two layers, in other words, there were two types of soil. In Model B, the SWC increased from the surface to the boundary of the two types of soil, and then a sudden decrease occurred at about 0.5 m. Below the boundary, the SWC increased with the increasing of depth until it reached the water table at 0.85 m. In Model C, we assumed that the surface reached saturated first and the SWC decreased from the surface to the boundary of the two types of soil. Below the boundary, the trend was similar to the change in SWC in Model B. The SWC increased with the increasing of depth until the water table at 0.85 m. The parameters of the models are listed in Table 1. Figure 3, Figure 4 and Figure 5 illustrate these three 1D surface models (blue solid lines in Figure 3, Figure 4 and Figure 5b) and their forward modeling responses (blue solid lines in Figure 3, Figure 4 and Figure 5a). For all the 1D forward modeling in this paper, the center frequency for the source was set to be 500 MHz.

3.1.1. Synthetic Data Without Noise

The performance of the proposed inversion scheme is first investigated on synthetic data without noise. For all these three models, the population size is set to be $N=10\times D$, where $D$ is the dimension of parameter space. Based on a few previous tests, the maximum number of iterations for data without noise is set to be 1000. The upper and lower boundaries for the search space depart 50% or more from the values of the true models for all the three models.

Figure 3, Figure 4 and Figure 5 provide a comprehensive illustration of the detailed inversion results for Model A, B, and C, respectively, when there is no noise in the data. Figure 3, Figure 4 and Figure 5a show the fitness between the calculated traces (depicted by red dashed lines) that were derived from the final models inverted using the proposed scheme and the observed traces (represented by blue solid lines) originating from the synthetic models. It can be seen that the calculated traces fit with the observed traces very closely. The differences between the calculated traces and the observed traces for each model are so minimal that they are almost imperceptible, which indicates that the proposed inversion scheme is highly accurate and reliable for applying to data without noise. Figure 3, Figure 4 and Figure 5b present a detailed comparison of the inverted SWC (represented by red dashed lines) against the synthetic SWC model (depicted by blue solid lines). It is evident that all the three true models are almost accurately recovered by the proposed inversion scheme at each layer. These comprehensive results not only demonstrate the effectiveness of the method but also validate that the proposed inversion scheme possesses the ability to accurately invert the SWC and layer thicknesses directly from the GPR waveform data. From Table 2, Table 3 and Table 4, it can also be seen that the inverted SWC and layer thicknesses match with the synthetic SWC and layer thicknesses almost perfectly at every layer, which further demonstrate the accuracy and the capability of the proposed inversion scheme for estimating the SWC and layer thicknesses based on the GPR waveform data with HHO.

3.1.2. Synthetic Data with Noise

When dealing with real data that has been collected under real-world field conditions, it is an unavoidable fact that interference (including noise and clutter) is inevitable. In this study, we only consider the influence of the noise. Consequently, we deliberately introduced 10% white Gaussian noise to the GPR data from Models A, B, and C to further verify the performance of the proposed inversion scheme. After this addition of noise, we tested the proposed inversion scheme using the data with noise to ensure that our method remains reliable under a less-than-ideal condition.

Figure 6, Figure 7 and Figure 8 illustrate the inversion results achieved through the proposed inversion scheme within a realistic setting. It becomes evident that the proposed inversion scheme still exhibits a robust performance. Although the data are contaminated with noise, the scheme is still capable of obtaining the acceptable solutions. In Figure 6, Figure 7 and Figure 8a, it can be clearly observed that the calculated traces (indicated by red dashed lines) derived from the inverted models still maintain a satisfactory alignment with the observed GPR traces (represented by blue solid lines). However, we can still see the impact of noise within the data. Specifically, when the amplitudes of the signal are diminished towards the end of the traces, it becomes challenging to identify the signal and the noise. Consequently, the fitness between the calculated traces and observed traces is observed to be less optimal towards the terminal parts of the traces compared to the earlier sections. This can also be demonstrated in Figure 6, Figure 7 and Figure 8b.

Figure 6, Figure 7 and Figure 8b present the results of SWC and layer thicknesses that have been inverted by the proposed inversion scheme based on the data with noise. It can be observed that the SWC and layer thicknesses still can be effectively and accurately recovered by the proposed inversion scheme even with the presence of noise. The inverted models (represented by red dashed lines) match well with the true models (indicated by blue solid lines) at each layer. Nevertheless, with the increases in the depth, the fitness becomes worse, particularly at the depths where the SWC almost reaches saturated. This trend corresponds to the end part of the GPR traces, which also suggests that the calculated traces do not align perfectly with the observed traces due to the disruptive effects of noise. Despite this, the overall results still serve to confirm the efficacy of the proposed inversion scheme in determining the SWC and the layer thicknesses by GPR waveform data, even if the data are impacted by noise. However, it is important to note that when compared to the inversion results derived from data without noise, the inversion results derived from noisy data show a discernible decrease in accuracy and reliability.

3.1.3. Analysis for Convergence Behavior

In order to examine the convergence behavior of HHO, Figure 9 displays the fitness evolution during HHO iterations using synthetic data both without and with noise from Model C. This offers significant understanding of the efficiency of the proposed inversion scheme. For data without noise, the fitness curves gradually converge to zero. The fitness curves gradually converge to a constant value for data with noise, demonstrating the resilience of the inversion scheme in the presence of disturbances. These findings indicate that the HHO algorithm simultaneously exhibits a high capacity for avoiding local optima and rapid convergence, which are essential characteristics for optimization problems. The ability to maintain a steady convergence rate even when faced with noisy data demonstrate the algorithm’s versatility and potential for application.

3.2. Comparisons with Other Algorithms

In order to verify the effectiveness evaluation of the proposed inversion scheme with HHO, we also performed a comparative analysis for the results of the HHO algorithm with other algorithms including the particle swarm optimization (PSO) and grey wolf optimizer (GWO) [53].

The synthetic data without and with noise which originated from Model C used in HHO inversion were, respectively, utilized to estimate the SWC and layer thicknesses from the GPR waveform data using PSO and GWO. For the sake of comparison and to ensure consistency, we applied the same search space and parameters, including the population size and the maximum number of iterations, as those that were used for the HHO inversion. Specifically, when conducting the inversion with the PSO algorithm, the inertia weighting factor was set to be $w=0.8$, and the scaling factors were selected to be $c_1=1.8$ and $c_2=2.0$, respectively.

Figure 10 and Figure 11 display the comparison of the inversion results obtained from three distinct optimization algorithms: HHO, PSO, and GWO. As depicted in Figure 10a, it can be observed that the calculated trace derived from the inversion results inverted by the HHO algorithm (represented by red solid line) exhibits a superior fit to the observed trace (indicated by green star line) in comparison to the calculated trace derived from the inversion results inverted by the PSO algorithm (represented by yellow solid line) and the GWO algorithm (represented by blue solid line). The red solid line almost perfectly aligns with the green star line at all significant peaks and valleys, demonstrating a high degree of precision in capturing the dynamic variations in the trace. In contrast, the blue solid line tracks the green star line with less accuracy than the red solid line, yet it still outperforms the yellow solid line in terms of following the observed trace. The phenomena suggests that the HHO algorithm is adept at accurately capturing the dynamic changes in the trace. When examining the entire trace, the red solid line consistently shows a higher degree of conformity with the shape of the trace compared to both the blue solid line and the yellow solid line, which indicates that the HHO algorithm possesses a higher approximation accuracy than both the GWO and PSO algorithm, highlighting its potential superiority in handling such optimization tasks. In Figure 11a, results similar to those in Figure 10a can be seen. For trace fitting, the HHO algorithm still exhibits a superior fit to the observed trace, even though the noise exists. Compared to both the blue and yellow solid lines, the red solid line matches better than the other two solid line and aligns with the green star line more accurately than the blue and yellow solid lines. This phenomena further demonstrates the potential superiority of the HHO algorithm although the existence of the noise.

Figure 10b shows a detailed comparison between the true SWC model and the SWC models inverted by using three different optimization algorithms: HHO, PSO, and GWO. Upon inspecting the degree of overlap among the four curves, it becomes evident that the red dashed line, representing the SWC inverted by the HHO, aligns more closely with the true SWC model. Conversely, the yellow dashed curve representing the SWC inverted by PSO and the blue dashed curve representing the SWC inverted by GWO exhibit more significant deviations from the true models. This finding further demonstrates the superior accuracy of the HHO algorithm in determining the SWC and layer thicknesses from GPR waveform data. Furthermore, this also implies that in terms of the correlation between SWC and depth, the SWC and layer thicknesses model inverted by the HHO algorithm is more capable of simulating real-world conditions with a higher level of precision. In Figure 11b, the inversion results are worse compared to the inversion results based on data without noise for the PSO and GWO algorithms. However, the inversion results of HHO still match closely with the true SWC model. This further demonstrates the superiority of the HHO algorithm for inverting SWC and layer thicknesses from data with noise.

At the same time, we also calculated the standard deviation (STD) and the overall average error (OAE) between the true model and the inverted model for these three optimization algorithms in Table 5 and Table 6. After thorough analysis, it can be seen that the HHO algorithm exhibited the lowest values for both the STD and the OAE when compared to the other two algorithms. This result provides compelling evidence supporting the performance of the HHO algorithm in estimating the SWC and the layer thicknesses through the GPR waveform inversion. However, we can also see that the STD and OAE are all larger than that for data without noise.

Figure 10c and Figure 11c display the comparison of the convergence behavior among these three optimization algorithms. It can be observed that the red solid line descends to a stable value more rapidly than the blue solid line and the yellow line, which means that the HHO algorithm can achieve convergence faster than the GWO and PSO algorithms. In the case of PSO in Figure 10c, the yellow solid line first reaches a stable value approximately after the 60th iteration. For GWO, the blue solid line nears stability around the 1900th iteration; however, the blue solid line can reach a lower value than the yellow solid line after approximately the 1500th iteration. In the context of optimization problems, a lower fitness value usually implies a solution that is closer to the optimum or a superior solution. The PSO algorithm may become stuck in a local minimum. This underscores that the GWO algorithm can gain more precise inversion results compared to the PSO algorithm, but the efficiency of the GWO is lower than that of the PSO. The red solid line demonstrates that the HHO not only reaches a stable value faster than the other two algorithms, after approximately the 400th iteration, but that it also attains a lower value compared to the GWO and the PSO. However, in Figure 11c, the yellow solid line first reaches a stable value after approximately the 1000th iteration. The existence of noise has a great influence on the PSO algorithm. The HHO and GWO algorithm both show good noise resistance. These results indicate that the HHO algorithm is capable of finding solutions of higher quality.

Overall, it is evident that the HHO algorithm performs better than both the GWO and the PSO algorithms regarding precision, model inversion, and convergence speed. We think that the mechanism of HHO is superior to that of GWO and PSO. First, the GWO is superior to PSO because the search of the GWO is based on a strict internal hierarchy, unlike the PSO which randomly initializes the particles of the population. And the grey wolf group is led by the three best leaders to approach the prey, which ensures that the GWO can better search for the true value [53]. For the HHO algorithm, it has a strong dynamic searching ability. By switching between the exploration phase and the development phase, the algorithm can effectively find the optimal solution in the solution space. In addition, the HHO introduced a variety of different hunting strategies, including multiple encirclement strategies, tracking, and raiding. This diversified search strategy enables HHO to better adapt to complex optimization problems and enhance the global search capability. Therefore, in the early stage, the HHO can achieve convergence faster and maintain a lower fitness value than the other two algorithms. This analysis offers significant insights for selecting the appropriate optimization algorithm.

3.3. Experimental Field Data Application

In order to assess the potential and practical use of the proposed inversion scheme, we implemented it on GPR data acquired at SCERES, an experimental site situated on the CNRS campus in Strasbourg-Cronenbourg, France. Figure 12 displays the underground structure of SCERES. It features a controlled artificial aquifer within a lined concrete basin. The edges of the basin are sealed to avoid any leakages to the outside. The inner dimensions of the basin are 24 m $\times$ 12 m $\times$ 3 m. The basin is equipped with the necessary tools for aquifer system and water table management. The basin is filled with three types of sand, and the specific parameters of the basin are detailed in Table 7. The water table lies roughly 0.85 m beneath the surface of the basin.

We used a MALA RAMAC system (which is produced by MALA company located at Stockholm, Sweden) equipped with 500 MHz antennas to collect GPR data across the basin, from one side to the other. Table 8 lists the technical index for the 500 MHz antennas. The entire length of the survey line is 24 m and we configured the system to collect one trace every 0.02 m. A DC filter and a shift of 20 samples are applied to the GPR profile. Figure 13 displays the processed GPR radargram [53]. We selected the 700th, 750th, 800th, and 850th traces as the observed traces for subsequent inversion based on experimental field data. To simulate synthetic GPR traces, it is essential to give a source wavelet. In this process, we first calculate the source wavelet by averaging the first arrivals from all traces, applying a Blackman window for weighting [61]. Subsequently, we employed a Butterworth band-pass filter to refine the source wavelet further. We took the 700th trace as an example. Figure 14 displays both the source wavelet and its corresponding spectrum.

For the experimental field data, we utilized a comparable inversion scheme as that employed for the synthetic data. The parameters that required determination were the SWC and the layer thickness. The SWC measured by Sentek sensors and the depth of each Sentek sensor were used as prior information [62]. Based on this prior information, the upper and lower boundaries for the search areas were established. Table 9 outlines the search space and inversion results for the experimental field data. Throughout the inversion, we adopted an 11-layer subsurface configuration, based on the placement of the Sentek sensors, to execute the inversion using the proposed inversion scheme.

Figure 15, Figure 16, Figure 17 and Figure 18 show the inversion results for the field data. From Figure 15, Figure 16, Figure 17 and Figure 18a, it can be seen that all the calculated GPR traces (represented by red dashed line) derived from the results inverted by the proposed inversion scheme align quite well the observed GPR traces (represented by blue solid line). Figure 15, Figure 16, Figure 17 and Figure 18b contrasts the inverted SWC and layer thicknesses obtained through the proposed inversion scheme, with the measured SWC by Sentek sensors. As illustrated in Figure 15, Figure 16, Figure 17 and Figure 18b, the inverted models (indicated by red dashed line) based on the proposed inversion scheme correspond with the measured SWC captured by Sentek sensors (indicated by blue solid line). They exhibit a comparable trend. Until a depth of approximately 0.55 m (marking the boundary between the two sandy layers), the SWC increases with the increase in depth. Approximately at 0.55 m, the SWC diminishes abruptly, potentially due to the influence of the boundary and the alterations in the medium. In the subsequent layers, the SWC continues to ascend with the increase in depth until it reaches the water table. There is a slight variation in the inverted model around 0.85 m. This may stem from the noise being regarded as signals that were fitted during the inversion. Nevertheless, the overall trend remains similar. Additionally, Figure 15, Figure 16, Figure 17 and Figure 18c demonstrate that all the fitness curves diminish at the beginning of the iterations, and then subsequently converge to a relatively constant value following the 1400th iteration. This convergence behavior is indicative of the robustness and stability of the inversion algorithm as it iteratively refines the calculated traces to better fit the observed traces. All inversion results further demonstrate the effectiveness of the proposed inversion scheme. However, we observed variations in the results obtained from different traces. We believe that the SWC exhibits not only vertical heterogeneity but also horizontal heterogeneity. Future studies should take this heterogeneity into account.

4. Conclusions

This study proposed a novel and powerful GPR FWI scheme based on the HHO algorithm to estimate the SWC and the layer thicknesses with a fixing quality factor. The HHO algorithm is inspired by the social behavior of Harris’s Hawks, mimicking their hunting behavior. Three main steps of exploration, the transition from exploration to exploitation, and exploitation, are used to realize the SWC and layer thicknesses estimation based on GPR waveform inversion. The proposed inversion scheme was first tested on three synthetic cases. Each case includes data without and with noise. In order to simulate more real-world scenarios where there is limited a priori information, we set wider boundaries for inversion search space. For further verification, we also compared the inversion results inverted from the HHO algorithm with the results from the PSO and GWO algorithms. The inversion results from the synthetic data demonstrate that the HHO algorithm is highly capable and converges rapidly to estimate the SWC and layer thicknesses from GPR waveform data. Subsequently, we further applied the proposed inversion scheme to field data collected from an experimental site at the CNRS campus in Strasbourg-Cronenbourg, France. For the experimental field data, we adopted the SWC measured by Sentek sensors and the depth of each Sentek sensor as a priori information. During the test with the experimental field data, the calculated trace closely matched the observed trace, and the inverted SWC showed a similar trend to the SWC measured by the Sentek sensors. The fitness curves from all tests on both synthetic and experimental field data demonstrates that the proposed scheme effectively avoids high local optima and achieves very fast convergence. Hence, we strongly suggest using the HHO algorithm for the SWC estimation based on GPR FWI.

This research is the initial endeavor to apply the HHO algorithm for GPR waveform inversion to invert the SWC and layer thicknesses. It also offers significant insights into the effectiveness of HHO in various geophysical inversions. In the future, we plan to improve our method to enhance the robustness and efficiency of HHO. At the same time, we also plan to consider the influence of more interference like clutter and further test the proposed inversion scheme on data with the interference of clutter. In addition, we will also continue to evaluate other cutting-edge algorithms to address the same issue. Furthermore, this technique can also be extended to challenges in other areas, including engineering, hydrology, and archeology, leveraging GPR technology.

Author Contributions

Conceptualization, H.Q. and M.Z.; methodology, H.Q., M.Z. and M.B.; software, H.Q., M.Z. and M.B.; validation, H.Q.; formal analysis, H.Q. and M.Z.; investigation, H.Q. and M.Z.; resources, H.Q.; data curation, H.Q. and M.Z.; writing—original draft preparation, H.Q.; writing—review and editing, M.Z. and M.B.; visualization, H.Q.; supervision, M.Z.; project administration, H.Q.; funding acquisition, H.Q. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

The authors thank Nolwenn Lesparre, Benjamin Belfort, Pascal Friedmann, Jean-François Girard, François Lehmann, and Olivier Razakarisoa for helping us to acquire GPR data. This experiment was funded by the ITES laboratory via an internal call. The authors are very grateful to the anonymous reviewers for their valuable suggestions and comments that helped to improve the quality of the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 The description for different phases of HHO [54].

Figure 2 The flow for inverting the SWC from GPR waveform inversion with HHO algorithm.

Figure 3 Inversion results of data without noise for Model A: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 4 Inversion results of data without noise for Model B: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 5 Inversion results of data without noise for Model C: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 6 Inversion results of data with noise for Model A: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 7 Inversion results of data with noise for Model B: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 8 Inversion results of data with noise for Model C: (a) observed GPR trace (blue solid line) and calculated GPR response (red dashed line); (b) true SWC (blue solid line) and inverted SWC (red dashed line).

Figure 9 Fitness behavior of Model C using the HHO algorithm: (a) fitness behavior for data without noise; (b) fitness behavior for data with noise.

Figure 10 Comparison of inverted results of data without noise for Model C with HHO, PSO, and GWO: (a) comparison of calculated traces from inverted models; (b) comparison of fitness behavior; (c) comparison of inverted results.

Figure 11 Comparison of inverted results of data with noise for Model C with HHO, PSO, and GWO: (a) comparison of calculated traces from inverted models; (b) comparison of fitness behavior; (c) comparison of inverted results.

Figure 12 (a) Illustration for experimental site; (b) illustration for the subsurface structure of the experimental site (Modified by [60]).

Figure 13 Processed observed GPR profile from experimental site. The red solid lines represent the position of the selected traces, which are regarded as the observed traces (the 700th, 750th, 800th, and 850th traces).

Figure 14 The source wavelet for real data modeling and the corresponding frequency spectrum for the 700th trace signal.

Figure 15 The inversion results of field data for 700th trace: (a) observed trace (blue solid line) and calculated trace (red dashed line); (b) measured SWC with Sentek (blue solid line) and inverted SWC (red dashed line); (c) fitness behavior of iteration.

Figure 16 The inversion results of field data for 750th trace: (a) observed trace (blue solid line) and calculated trace (red dashed line); (b) measured SWC with Sentek (blue solid line) and inverted SWC (red dashed line); (c) fitness behavior of iteration.

Figure 17 The inversion results of field data for 800th trace: (a) observed trace (blue solid line) and calculated trace (red dashed line); (b) measured SWC with Sentek (blue solid line) and inverted SWC (red dashed line); (c) fitness behavior of iteration.

Figure 18 The inversion results of field data for 850th trace: (a) observed trace (blue solid line) and calculated trace (red dashed line); (b) measured SWC with Sentek (blue solid line) and inverted SWC (red dashed line); (c) fitness behavior of iteration.

Table 1

The soil property parameters and search space for Model A, Model B, and Model C.

Layer	Parameters of Model A		Search Space of Model A		Parameters of Model B		Search Space of Model B		Parameters of Model C		Search Space of Model C
	SWC 1	z 2	SWC	z	SWC	z	SWC	z	SWC	z	SWC	z
1	0.05	0.15	0.03–0.08	0.08–0.23	0.37	0.15	0.19–0.56	0.08–0.23	0.07	0.15	0.04–0.11	0.08–0.23
2	0.05	0.20	0.03–0.08	0.10–0.30	0.21	0.20	0.10–0.31	0.10–0.30	0.18	0.20	0.09–0.27	0.10–0.30
3	0.07	0.10	0.03–0.10	0.05–0.15	0.25	0.10	0.13–0.38	0.05–0.15	0.31	0.10	0.15–0.46	0.05–0.15
4	0.10	0.20	0.05–0.16	0.10–0.30	0.07	0.20	0.03–0.10	0.10–0.30	0.06	0.20	0.03–0.08	0.10–0.30
5	0.20	0.10	0.10–0.30	0.05–0.15	0.20	0.10	0.10–0.30	0.05–0.15	0.16	0.10	0.08–0.23	0.05–0.15
6	0.43	0.10	0.22–0.65	0.05–0.15	0.43	0.10	0.22–0.65	0.05–0.15	0.43	0.10	0.22–0.65	0.05–0.15
7	0.43	0.10	0.22–0.65	0.05–0.15	0.43	0.10	0.22–0.65	0.05–0.15	0.43	0.10	0.22–0.65	0.05–0.15
8	0.43	$\infty$ 3	0.22–0.65	$\infty$	0.43	$\infty$	0.22–0.65	$\infty$	0.43	$\infty$	0.22–0.65	$\infty$

¹ SWC is the abbreviation of soil water content. ² z represents the layer thickness. ³ $\infty$ represents the infinite half space.

Table 2

The inversion results for Model A.

Layer	Model Parameters		Inversion Results for Data Without Noise		Inversion Results for Data with Noise
	SWC 1	z 2	SWC	z	SWC	z
1	0.05	0.15	0.05	0.15	0.05	0.15
2	0.05	0.20	0.06	0.20	0.06	0.20
3	0.07	0.10	0.07	0.10	0.06	0.10
4	0.10	0.20	0.10	0.20	0.10	0.21
5	0.20	0.10	0.20	0.10	0.19	0.10
6	0.43	0.10	0.43	0.10	0.40	0.09
7	0.43	0.10	0.43	0.10	0.43	0.10
8	0.43	$\infty$ 3	0.43	$\infty$	0.43	$\infty$

¹ SWC is the abbreviation of soil water content. ² z represents the layer thickness. ³ $\infty$ represents the infinite half space.

Table 3

The inversion results for Model B.

Layer	Model Parameters		Inversion Results for Data Without Noise		Inversion Results for Data with Noise
	SWC 1	z 2	SWC	z	SWC	z
1	0.37	0.15	0.37	0.15	0.37	0.15
2	0.21	0.20	0.20	0.20	0.20	0.19
3	0.25	0.10	0.25	0.10	0.25	0.09
4	0.07	0.20	0.07	0.20	0.07	0.20
5	0.20	0.10	0.20	0.10	0.20	0.10
6	0.43	0.10	0.43	0.10	0.43	0.09
7	0.43	0.10	0.43	0.09	0.43	0.10
8	0.43	$\infty$ 3	0.43	$\infty$	0.43	$\infty$

¹ SWC is the abbreviation of soil water content. ² z represents the layer thickness. ³ $\infty$ represents the infinite half space.

Table 4

The inversion results for Model C.

Layer	Model Parameters		Inversion Results for Data Without Noise		Inversion Results for Data with Noise
	SWC 1	z 2	SWC	z	SWC	z
1	0.07	0.15	0.07	0.15	0.07	0.14
2	0.18	0.20	0.18	0.20	0.18	0.20
3	0.31	0.10	0.31	0.10	0.29	0.10
4	0.06	0.20	0.05	0.19	0.05	0.20
5	0.16	0.10	0.16	0.10	0.16	0.10
6	0.43	0.10	0.43	0.10	0.42	0.09
7	0.43	0.10	0.43	0.10	0.43	0.09
8	0.43	$\infty$ 3	0.43	$\infty$	0.43	$\infty$

¹ SWC is the abbreviation of soil water content. ² z represents the layer thickness. ³ $\infty$ represents the infinite half space.

Table 5

The comparison of STD and OAE of three algorithms for data without noise from Model C.

	HHO		PSO		GWO
	STD 1	OAE 2	STD	OAE	STD	OAE
Model C	12.99%	0.03%	15.24%	6.86%	14.95%	4.29%

¹ STD is the standard deviation. ² OAE is the overall average error.

Table 6

The comparison of STD and OAE of three algorithms for data with noise from Model C.

	HHO		PSO		GWO
	STD 1	OAE 2	STD	OAE	STD	OAE
Model C	13.07%	0.17%	20.52%	12.97%	12.67%	1.56%

¹ STD is the standard deviation. ² OAE is the overall average error.

Table 7

Parameters for the basin at the experimental site.

Layer	Thickness (m)	Porosity (%)
1	0.5	43
2	2.5	40
3	0.5	38

Table 8

Technical indexes for the 500 MHz antennas.

Items	Technical Index
Frequency of antennas	500 MHz
Maximum detection depth	6 m
Radial resolution	0.05 m
Offset	0.18 m
Size	0.5 × 0.3 × 0.16 m
Weight	5.0 kg

Table 9

The upper and lower limitations for search space and inversion results for field data.

Layer	Search Space Limitations		Inversion Results
	SWC 1	z 2	SWC	z
1	0.04–0.11	0.08–0.23	0.07	0.18
2	0.07–0.21	0.05–0.15	0.08	0.15
3	0.10–0.30	0.05–0.15	0.22	0.05
4	0.13–0.38	0.05–0.15	0.32	0.07
5	0.04–0.05	0.05–0.15	0.04	0.08
6	0.04–0.11	0.05–0.15	0.10	0.14
7	0.11–0.32	0.05–0.15	0.14	0.08
8	0.20–0.43	0.05–0.15	0.37	0.08
9	0.20–0.43	0.05–0.15	0.21	0.12
10	0.20–0.43	0.05–0.15	0.36	0.05
11	0.20–0.43	$\infty$ 3	0.41	$\infty$

¹ SWC is the abbreviation of soil water content. ² z represents the layer thickness. ³ $\infty$ represents the infinite half space.

References

1. Vereecken, H.; Huisman, J.A.; Bogena, H.; Vanderborght, J.; Vrugt, J.A.; Hopmans, J.W. On the value of soil moisture measurements in vadose zone hydrology: A review. Water Resour. Res.; 2008; 44, W00D06. [DOI: https://dx.doi.org/10.1029/2008WR006829]

2. Dobriyal, P.; Qureshi, A.; Badola, R.; Hussain, S.A. A review of the methods available for estimating soil moisture and its implications for water resource management. J. Hydrol.; 2012; 458–459, pp. 110-117. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2012.06.021]

3. Blonquist, J.M.; Jones, S.B.; Robinson, D.A. Precise irrigation scheduling for turfgrass using a subsurface electromagnetic soil moisture sensor. Agric. Water Manag.; 2006; 84, pp. 153-165. [DOI: https://dx.doi.org/10.1016/j.agwat.2006.01.014]

4. Cai, G.; Vanderborght, J.; Langensiepen, M.; Schnepf, A.; Hüging, H.; Vereecken, H. Root growth, water uptake, and sap flow of winter wheat in response to different soil water conditions. Hydrol. Earth Syst. Sci.; 2018; 22, pp. 2449-2470. [DOI: https://dx.doi.org/10.5194/hess-22-2449-2018]

5. Robinson, D.A.; Campbell, C.S.; Hopmans, J.W.; Hornbuckle, B.K.; Jones, S.B.; Knight, R.; Ogden, F.; Selker, J.; Wendroth, O. Soil Moisture Measurement for Ecological and Hydrological Watershed-Scale Observatories: A Review. Vadose Zone J.; 2008; 7, pp. 358-389. [DOI: https://dx.doi.org/10.2136/vzj2007.0143]

6. Klotzsche, A.; Jonard, F.; Looms, M.C.; van der Kruk, J.; Huisman, J.A. Measuring Soil Water Content with Ground Penetrating Radar: A Decade of Progress. Vadose Zone J.; 2018; 17, 180052. [DOI: https://dx.doi.org/10.2136/vzj2018.03.0052]

7. Huisman, J.A.; Hubbard, S.S.; Redman, J.D.; Annan, A.P. Measuring Soil Water Content with Ground Penetrating Radar: A Review. Vadose Zone J.; 2003; 2, pp. 476-491. [DOI: https://dx.doi.org/10.2136/vzj2003.4760]

8. Liu, X.; Chen, J.; Cui, X.; Liu, Q.; Cao, X.; Chen, X. Measurement of soil water content using ground-penetrating radar: A review of current methods. Int. J. Digit. Earth; 2019; 12, pp. 95-118. [DOI: https://dx.doi.org/10.1080/17538947.2017.1412520]

9. Zhang, M.; Feng, X.; Bano, M.; Xing, H.; Wang, T.; Liang, W.; Zhou, H.; Dong, Z.; An, Y.; Zhang, Y. Review of Ground Penetrating Radar Applications for Water Dynamics Studies in Unsaturated Zone. Remote Sens.; 2022; 14, 5993. [DOI: https://dx.doi.org/10.3390/rs14235993]

10. Ercoli, M.; Di Matteo, L.; Pauselli, C.; Mancinelli, P.; Frapiccini, S.; Talegalli, L.; Cannata, A. Integrated GPR and laboratory water content measures of sandy soils: From laboratory to field scale. Constr. Build. Mater.; 2018; 159, pp. 734-744. [DOI: https://dx.doi.org/10.1016/j.conbuildmat.2017.11.082]

11. Grote, K.; Hubbard, S.; Rubin, Y. GPR monitoring of volumetric water content in soils applied to highway construction and maintenance. Lead. Edge; 2002; 21, pp. 482-504. [DOI: https://dx.doi.org/10.1190/1.1481259]

12. Irving, J.; Knight, R.; Holliger, K. Estimation of the lateral correlation structure of subsurface water content from surface-based ground-penetrating radar reflection images. Water Resour. Res.; 2009; 45, W12404. [DOI: https://dx.doi.org/10.1029/2008WR007471]

13. Klenk, P.; Jaumann, S.; Roth, K. Quantitative high-resolution observations of soil water dynamics in a complicated architecture using time-lapse ground-penetrating radar. Hydrol. Earth Syst. Sci.; 2015; 19, pp. 1125-1139. [DOI: https://dx.doi.org/10.5194/hess-19-1125-2015]

14. Lunt, I.A.; Hubbard, S.S.; Rubin, Y. Soil moisture content estimation using ground-penetrating radar reflection data. J. Hydrol.; 2005; 307, pp. 254-269. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2004.10.014]

15. Moua, R.; Lesparre, N.; Girard, J.F.; Belfort, B.; Lehmann, F.; Younes, A. On the interest of ground penetrating radar data for the estimation of unsaturated soil parameters. EGUsphere; 2022; 2022, pp. 1-24. [DOI: https://dx.doi.org/10.5194/egusphere-2022-936]

16. Pettinelli, E.; Di Matteo, A.; Beaubien, S.E.; Mattei, E.; Lauro, S.E.; Galli, A.; Vannaroni, G. A controlled experiment to investigate the correlation between early-time signal attributes of ground-coupled radar and soil dielectric properties. J. Appl. Geophys.; 2014; 101, pp. 68-76. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2013.11.012]

17. Steelman, C.M.; Endres, A.L. An examination of direct ground wave soil moisture monitoring over an annual cycle of soil conditions. Water Resour. Res.; 2010; 46, W11533. [DOI: https://dx.doi.org/10.1029/2009WR008815]

18. Yu, Y.; Klotzsche, A.; Weihermüller, L.; Huisman, J.A.; Vanderborght, J.; Vereecken, H.; van der Kruk, J. Measuring vertical soil water content profiles by combining horizontal borehole and dispersive surface ground penetrating radar data. Near Surf. Geophys.; 2020; 18, pp. 275-294. [DOI: https://dx.doi.org/10.1002/nsg.12099]

19. Huisman, J.A.; Snepvangers, J.J.J.C.; Bouten, W.; Heuvelink, G.B.M. Monitoring Temporal Development of Spatial Soil Water Content Variation: Comparison of Ground Penetrating Radar and Time Domain Reflectometry. Vadose Zone J.; 2003; 2, pp. 519-529. [DOI: https://dx.doi.org/10.2136/vzj2003.5190]

20. Allroggen, N.; van Schaik, N.L.M.B.; Tronicke, J. 4D ground-penetrating radar during a plot scale dye tracer experiment. J. Appl. Geophys.; 2015; 118, pp. 139-144. [DOI: https://dx.doi.org/10.1016/j.jappgeo.2015.04.016]

21. Bradford, J.H. Measuring Water Content Heterogeneity Using Multifold GPR with Reflection Tomography. Vadose Zone J.; 2008; 7, pp. 184-193. [DOI: https://dx.doi.org/10.2136/vzj2006.0160]

22. Strobbia, C.; Cassiani, G. Multilayer ground-penetrating radar guided waves in shallow soil layers for estimating soil water content. Geophysics; 2007; 72, pp. J17-J29. [DOI: https://dx.doi.org/10.1190/1.2716374]

23. Hubbard, S.; Grote, K.; Rubin, Y. Mapping the volumetric soil water content of a California vineyard using high-frequency GPR ground wave data. Lead. Edge; 2002; 21, pp. 552-559. [DOI: https://dx.doi.org/10.1190/1.1490641]

24. Huisman, J.A.; Sperl, C.; Bouten, W.; Verstraten, J.M. Soil water content measurements at different scales: Accuracy of time domain reflectometry and ground-penetrating radar. J. Hydrol.; 2001; 245, pp. 48-58. [DOI: https://dx.doi.org/10.1016/S0022-1694(01)00336-5]

25. Iwasaki, T.; Kuroda, S.; Saito, H.; Tobe, Y.; Suzuki, K.; Fujimaki, H.; Inoue, M. Monitoring Infiltration Process Seamlessly Using Array Ground Penetrating Radar. Agric. Environ. Lett.; 2016; 1, 160002. [DOI: https://dx.doi.org/10.2134/ael2016.01.0002]

26. Saito, H.; Kuroda, S.; Iwasaki, T.; Sala, J.; Fujimaki, H. Estimating infiltration front depth using time-lapse multioffset gathers obtained from ground-penetrating-radar antenna array. Geophysics; 2021; 86, pp. WB51-WB59. [DOI: https://dx.doi.org/10.1190/geo2020-0590.1]

27. Binley, A.; Cassiani, G.; Middleton, R.; Winship, P. Vadose zone flow model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol.; 2002; 267, pp. 147-159. [DOI: https://dx.doi.org/10.1016/S0022-1694(02)00146-4]

28. Galagedara, L.W.; Parkin, G.W.; Redman, J.D.; Endres, A.L. Assessment of Soil Moisture Content Measured by Borehole GPR and TDR Under Transient Irrigation and Drainage. J. Environ. Eng. Geophys.; 2003; 8, pp. 77-86. [DOI: https://dx.doi.org/10.4133/JEEG8.2.77]

29. Klotzsche, A.; Lärm, L.; Vanderborght, J.; Cai, G.; Morandage, S.; Zörner, M.; Vereecken, H.; van der Kruk, J. Monitoring Soil Water Content Using Time-Lapse Horizontal Borehole GPR Data at the Field-Plot Scale. Vadose Zone J.; 2019; 18, 190044. [DOI: https://dx.doi.org/10.2136/vzj2019.05.0044]

30. Kowalsky, M.B.; Finsterle, S.; Rubin, Y. Estimating flow parameter distributions using ground-penetrating radar and hydrological measurements during transient flow in the vadose zone. Adv. Water Resour.; 2004; 27, pp. 583-599. [DOI: https://dx.doi.org/10.1016/j.advwatres.2004.03.003]

31. Rucker, D.F.; Ferré, T.P.A. Near-Surface Water Content Estimation with Borehole Ground Penetrating Radar Using Critically Refracted Waves. Vadose Zone J.; 2003; 2, pp. 247-252. [DOI: https://dx.doi.org/10.2136/vzj2003.2470]

32. Looms, M.C.; Jensen, K.H.; Binley, A.; Nielsen, L. Monitoring Unsaturated Flow and Transport Using Cross-Borehole Geophysical Methods. Vadose Zone J.; 2008; 7, pp. 227-237. [DOI: https://dx.doi.org/10.2136/vzj2006.0129]

33. Wijewardana, Y.G.N.S.; Galagedara, L.W. Estimation of spatio-temporal variability of soil water content in agricultural fields with ground penetrating radar. J. Hydrol.; 2010; 391, pp. 24-33. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2010.06.036]

34. Strobach, E.; Harris, B.D.; Dupuis, J.C.; Kepic, A.W. Time-lapse borehole radar for monitoring rainfall infiltration through podosol horizons in a sandy vadose zone. Water Resour. Res.; 2014; 50, pp. 2140-2163. [DOI: https://dx.doi.org/10.1002/2013WR014331]

35. Cheng, Q.; Su, Q.; Binley, A.; Liu, J.; Zhang, Z.; Chen, X. Estimation of Surface Soil Moisture by a Multi-Elevation UAV-Based Ground Penetrating Radar. Water Resour. Res.; 2023; 59, e2022WR032621. [DOI: https://dx.doi.org/10.1029/2022WR032621]

36. Jadoon, K.Z.; Lambot, S.; Scharnagl, B.; van der Kruk, J.; Slob, E.; Vereecken, H. Quantifying field-scale surface soil water content from proximal GPR signal inversion in the time domain. Near Surf. Geophys.; 2010; 8, pp. 483-491. [DOI: https://dx.doi.org/10.3997/1873-0604.2010036]

37. Jonard, F.; Weihermüller, L.; Vereecken, H.; Lambot, S. Accounting for soil surface roughness in the inversion of ultrawideband off-ground GPR signal for soil moisture retrieval. Geophysics; 2012; 77, pp. H1-H7. [DOI: https://dx.doi.org/10.1190/geo2011-0054.1]

38. Lambot, S.; Rhebergen, J.; van den Bosch, I.; Slob, E.C.; Vanclooster, M. Measuring the Soil Water Content Profile of a Sandy Soil with an Off-Ground Monostatic Ground Penetrating Radar. Vadose Zone J.; 2004; 3, pp. 1063-1071. [DOI: https://dx.doi.org/10.2136/vzj2004.1063]

39. Lambot, S.; Weihermüller, L.; Huisman, J.A.; Vereecken, H.; Vanclooster, M.; Slob, E.C. Analysis of air-launched ground-penetrating radar techniques to measure the soil surface water content. Water Resour. Res.; 2006; 42, W11403. [DOI: https://dx.doi.org/10.1029/2006WR005097]

40. Minet, J.; Bogaert, P.; Vanclooster, M.; Lambot, S. Validation of ground penetrating radar full-waveform inversion for field scale soil moisture mapping. J. Hydrol.; 2012; 424–425, pp. 112-123. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2011.12.034]

41. Serbin, G.; Or, D. Ground-penetrating radar measurement of soil water content dynamics using a suspended horn antenna. IEEE Trans. Geosci. Remote Sens.; 2004; 42, pp. 1695-1705. [DOI: https://dx.doi.org/10.1109/TGRS.2004.831693]

42. Tran, A.P.; Ardekani, M.R.M.; Lambot, S. Coupling of dielectric mixing models with full-wave ground-penetrating radar signal inversion for sandy-soil-moisture estimation. Geophysics; 2012; 77, pp. H33-H44. [DOI: https://dx.doi.org/10.1190/geo2011-0100.1]

43. Wu, K.; Rodriguez, G.A.; Zajc, M.; Jacquemin, E.; Clément, M.; De Coster, A.; Lambot, S. A new drone-borne GPR for soil moisture mapping. Remote Sens. Environ.; 2019; 235, 111456. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111456]

44. Tarantola, A. Inversion of seismic reflection data in the acoustic approximation. Geophysics; 1984; 49, pp. 1259-1266. [DOI: https://dx.doi.org/10.1190/1.1441754]

45. Klotzsche, A.; Vereecken, H.; van der Kruk, J. Review of crosshole ground-penetrating radar full-waveform inversion of experimental data: Recent developments, challenges, and pitfalls. Geophysics; 2019; 84, pp. H13-H28. [DOI: https://dx.doi.org/10.1190/geo2018-0597.1]

46. Kruk, J.v.d.; Liu, T.; Mozaffari, A.; Gueting, N.; Klotzsche, A.; Vereecken, H.; Warren, C.; Giannopoulos, A. GPR full-waveform inversion, recent developments, and future opportunities. Proceedings of the 2018 17th International Conference on Ground Penetrating Radar (GPR); Rapperswil, Switzerland, 18–21 June 2018; pp. 1-6.

47. Slob, E.; Sato, M.; Olhoeft, G. Surface and borehole ground-penetrating-radar developments. Geophysics; 2010; 75, pp. 75A103-75A120. [DOI: https://dx.doi.org/10.1190/1.3480619]

48. Lambot, S.; Slob, E.; Chavarro, D.; Lubczynski, M.; Vereecken, H. Measuring soil surface water content in irrigated areas of southern Tunisia using full-waveform inversion of proximal GPR data. Near Surf. Geophys.; 2008; 6, pp. 403-410. [DOI: https://dx.doi.org/10.3997/1873-0604.2008028]

49. Lambot, S.; Slob, E.; Rhebergen, J.; Lopera, O.; Jadoon, K.Z.; Vereecken, H. Remote Estimation of the Hydraulic Properties of a Sand Using Full-Waveform Integrated Hydrogeophysical Inversion of Time-Lapse, Off-Ground GPR Data. Vadose Zone J.; 2009; 8, pp. 743-754. [DOI: https://dx.doi.org/10.2136/vzj2008.0058]

50. Busch, S.; van der Kruk, J.; Bikowski, J.; Vereecken, H. Quantitative conductivity and permittivity estimation using full-waveform inversion of on-ground GPR data. Geophysics; 2012; 77, pp. H79-H91. [DOI: https://dx.doi.org/10.1190/geo2012-0045.1]

51. Lavoué, F.; Brossier, R.; Métivier, L.; Garambois, S.; Virieux, J. Two-dimensional permittivity and conductivity imaging by full waveform inversion of multioffset GPR data: A frequency-domain quasi-Newton approach. Geophys. J. Int.; 2014; 197, pp. 248-268. [DOI: https://dx.doi.org/10.1093/gji/ggt528]

52. Wang, X.; Feng, D. Multiparameter Full-Waveform Inversion of 3-D On-Ground GPR With a Modified Total Variation Regularization Scheme. IEEE Geosci. Remote Sens. Lett.; 2021; 18, pp. 466-470. [DOI: https://dx.doi.org/10.1109/LGRS.2020.2976146]

53. Zhang, M.H.; Feng, X.; Bano, M.; Liu, C.; Liu, Q.; Wang, X. Soil water content estimation by using ground penetrating radar data full waveform inversion with grey wolf optimizer algorithm. Vadose Zone J.; 2025; 24, e20379. [DOI: https://dx.doi.org/10.1002/vzj2.20379]

54. Heidari, A.A.; Mirjalili, S.; Faris, H.; Aljarah, I.; Mafarja, M.; Chen, H. Harris hawks optimization: Algorithm and applications. Future Gener. Comput. Syst.; 2019; 97, pp. 849-872. [DOI: https://dx.doi.org/10.1016/j.future.2019.02.028]

55. Ouyang, C.; Liao, C.; Zhu, D.; Zheng, Y.; Zhou, C.; Li, T. Integrated improved Harris hawks optimization for global and engineering optimization. Sci. Rep.; 2024; 14, 7445. [DOI: https://dx.doi.org/10.1038/s41598-024-58029-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38548845]

56. Tripathy, B.K.; Reddy Maddikunta, P.K.; Pham, Q.-V.; Gadekallu, T.R.; Dev, K.; Pandya, S.; ElHalawany, B.M.; Sarfraz, S. Harris Hawk Optimization: A Survey onVariants and Applications. Comput. Intell. Neurosci.; 2022; 2022, 2218594. [DOI: https://dx.doi.org/10.1155/2022/2218594]

57. Hussien, A.G.; Abualigah, L.; Abu Zitar, R.; Hashim, F.A.; Amin, M.; Saber, A.; Almotairi, K.H.; Gandomi, A.H. Recent Advances in Harris Hawks Optimization: A Comparative Study and Applications. Electronics; 2022; 11, 1919. [DOI: https://dx.doi.org/10.3390/electronics11121919]

58. Topp, G.C.; Davis, J.L.; Annan, A.P. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resour. Res.; 1980; 16, pp. 574-582. [DOI: https://dx.doi.org/10.1029/WR016i003p00574]

59. Bano, M. Modelling of GPR waves for lossy media obeying a complex power law of frequency for dielectric permittivity. Geophys. Prospect.; 2004; 52, pp. 11-26. [DOI: https://dx.doi.org/10.1046/j.1365-2478.2004.00397.x]

60. Zhang, M.H.; Bano, M.; Feng, X.; Lesparre, N.; Girard, J.F.; Razakarisoa, O.; Belfort, B.; Lehmann, F.; Friedmann, P. Estimating water content of unsaturated sandy soils by GPR during a drainage experiment. IOP Conf. Ser. Earth Environ. Sci.; 2021; 660, 012018. [DOI: https://dx.doi.org/10.1088/1755-1315/660/1/012018]

61. Guillemoteau, J.; Bano, M.; Dujardin, J.-R. Influence of grain size, shape and compaction on georadar waves: Examples of aeolian dunes. Geophys. J. Int.; 2012; 190, pp. 1455-1463. [DOI: https://dx.doi.org/10.1111/j.1365-246X.2012.05577.x]

62. Jellali, S.; Muntzer, P.; Razakarisoa, O.; Schäfer, G. Large Scale Experiment on Transport of Trichloroethylene in a Controlled Aquifer. Transp. Porous Media; 2001; 44, pp. 145-163. [DOI: https://dx.doi.org/10.1023/A:1010652230922]