Full-waveform inversion of ground-penetrating radar data in frequency-dependent media involving permittivity attenuation

Full-waveform inversion (FWI) of ground-penetrating radar (GPR) data has received particular attention in the past decade because it can provide high-resolution subsurface models of dielectric permittivity and electrical conductivity. In most GPR FWIs, these two parameters are regarded as frequency independent, which may lead to false estimates if they strongly depend on frequency, such as in shallow weathered zones. In this study, we develop frequency-dependent GPR FWI to solve this problem. Using the τ -method introduced in the research of viscoelastic waves, we deﬁne the permittivity attenuation parameter to quantify the attenuation resulting from the complex permittivity and to modify time-domain Maxwell’s equations. The new equations are self-adjoint so that we can use the same forward engine to back-propagate the adjoint sources and easily derive model gradients in GPR FWI. Frequency dependence analysis shows that permittivity attenuation acts as a low-pass ﬁlter, distorting the waveform and decaying the amplitude of the electromagnetic waves. The 2-D synthetic examples illustrate that permittivity attenuation has low sensitivity to the surface multioffset GPR data but is necessary for a good reconstruction of permittivity and conductivity models in frequency-dependent GPR FWI. As a comparison, frequency-independent GPR FWI produces more model artefacts and hardly reconstructs conductivity models dominated by permittivity attenuation. The 2-D ﬁeld example shows that both FWIs reveal a triangle permittivity anomaly which proves to be a reﬁlled trench. However, frequency-dependent GPR FWI provides a better ﬁt to the observed data and a more robust conductivity reconstruction in a high permittivity attenuation environment. Our GPR FWI results are consistent with previous GPR and shallow-seismic measurements. This research greatly expands the application of GPR FWI in more complicated media. the performance of frequency-dependent GPR FWI. Transverse magnetic (TM) mode waves are used in the simulation and inversion of GPR surface recordings. In frequency-independent media, there are waveﬁeld components H x , H z and E y and reconstructed model parameters ε and σ . In frequency-dependent media, additional memory variable r yl and reconstructed model parameters τ σ and τ ε are added with L = 1. The magnetic permeability μ is constant and equal to its value of vacuum. We set τ σ = 0 s based on the analysis in the previous section.

FWI of GPR data in frequency-dependent media 505 standard linear solid (SLS) mechanisms, as shown, for example in Blanch et al. (1995). These equations use similar physical models to characterize the frequency dependence of electrical parameters and seismic moduli (Carcione & Cavallini 1995). It implies that GPR FWI can refer to viscoelastic FWI already developed in the seismic community, such as Fabien-Ouellet et al. (2017) and Jiang (2019). However, Maxwell's equations differ from the viscoelastic equations because the complex permittivity and conductivity always occur in a combined form and both cause attenuation (Turner & Siggins 1994;Bradford 2007). In viscoelastic case, the magnitude of dispersion and attenuation of a single wave (P or S waves) can be described by a single parameter (τ P or τ S ), which corresponds to the case of considering only permittivity attenuation or conductivity attenuation in EM waves (Blanch et al. 1995).
In the last decade, most GPR FWIs still used frequency-independent sensitivity kernels when updating model parameters (Klotzsche et al. 2019). There are two reasons for this. One is that simultaneous reconstruction of permittivity and conductivity is already a challenging task in frequency-independent GPR FWI, while studies of frequency-dependent GPR FWI that require estimation of more parameters take a back seat. Another reason is that frequency-independent GPR FWI, when combined with source signal estimation, can handle weak permittivity attenuation environments where quality factor (Q) ≥ 20 (Belina et al. 2012). In the study of Belina et al. (2012), since their model gradients were practically the same as those of the frequency-independent GPR FWI and the static conductivity was ignored, researchers did a permittivity-only inversion with a priori Q model. To better image the realistic materials with Q < 20 and investigate crosstalk between multiple parameters, frequency-dependent GPR FWI must be developed.
FWI is also known as an inverse scattering problem in the microwave community. Winters et al. (2006) proposed the time-domain inverse scattering technique to estimate the frequency-dependent average dielectric properties, using a short relaxation time approximation in Debye scatterers. Based on their work, Papadopoulos & Rekanos (2011) introduced an auxiliary differential equation (ADE) with the polarization current density, which extended the feasibility of microwave imaging. Deng et al. (2021) reported the high-performance computation of EM FWI in the single-pole Debye model, using the theoretical basis of the inverse scattering technique. These studies suggested that we can similarly develop frequency-dependent GPR FWI. Nevertheless, they did not give an explicit measure of permittivity attenuation in EM wave propagation, which makes the inversion algorithm less intuitive than viscoelastic FWI involving seismic attenuation.
In this paper, we use the τ -method, with reference to the viscoelastic waves (Blanch et al. 1995), to quantify the attenuation of EM waves caused by the complex permittivity. Based on this, we modify the time-domain Maxwell's equations proposed by Bergmann et al. (1998) and implement the frequency-dependent GPR FWI. With a frequency dependence analysis, we show how each model parameter affects EM wave propagation, including amplitude attenuation and velocity dispersion. Finally, we apply 2-D GPR FWI in the synthetic examples and field examples and compare the performance of the FWI considering and not considering permittivity attenuation. Since surface-based GPR plays an important role for characterizing near-surface targets, we perform a surface-based GPR FWI in this study to investigate the applicability and limitations of using reflected waves.

Forward problem
The propagation of EM waves in heterogeneous media can be expressed by Maxwell's equations as follows: with three generalized constitutive relations: where B is the magnetic flux, H is the magnetic field, D is the electric displacement and E is the electric field. J c is the conduction current density and J e is the electric current sources. ∇ is the Laplace operator, × is the curl operator and * is the time convolution. The magnetic permeability μ is essentially frequency independent, while the dielectric permittivity ε and electrical conductivity σ are described as complex frequency-dependent quantities (Carcione 1996;Bergmann et al. 1998). The attenuation and dispersion processes of EM waves depend on the ε and σ of the media. Substituting eq. (2) into eq. (1), Maxwell's equations become In order to solve eq. (3) explicitly in the time domain, we must calculate the convolutions ε * ∂ t E and σ * E. We use the Debye model (eq. 4) to approximate the relaxation function of the permittivity in the time domain (Carcione 1996), i.e. the frequency-dependent function Downloaded from https://academic.oup.com/gji/article/232/1/504/6670780 by KIT Library user on 25 October 2022 in the frequency domain.
where τ El and τ Dl are the relaxation times of the electric field and displacement of the lth Debye model, respectively. τ El < τ Dl for the lth mechanism. τ Dl = 1/(2π f l ) where f l is the relaxation frequency of the lth Debye peak and L is the number of relaxation mechanisms. H(t) is the Heaviside function. ε s is the static dielectric permittivity, corresponding to the permittivity in the frequency-independent media. Following the τ -method proposed by Blanch et al. (1995), we define the permittivity attenuation τ ε = 1 − τ El /τ Dl as the attenuation level of the lth Debye model for EM waves. By doing so, ε can be rewritten as In eq. (5), one can see that τ ε is independent of the lth Debye model, which means the same τ ε can be used for different Debye models. The newly defined variable is dimensionless (0 ≤ τ ε < 1) and more straightforward for identifying attenuation and dispersion caused by the complex permittivity. Using the Debye model, the convolution ε * ∂ t E becomes Due to space limitations in the main content, we omit the computational process of convolution (readers can find details of eqs 6, 7, 10, and 11 in Appendix A). We introduce an ADE with a memory variable r l of the lth mechanism corresponding to E: Note that the definition of r l is slightly different from the memory variable in Carcione (1996) and Bergmann et al. (1998). Similar to viscoelastic wave equations described in Bohlen (2002), we incorporate ε s into r l . By doing so, we can easily develop self-adjoint wave equations for GPR FWI (shown in Section 2.2). If only one Debye model is used, the memory variable r l is related to the polarization current density J used in the microwave imaging of Papadopoulos & Rekanos (2011) by On the other hand, the relaxation function of the conductivity can be described by a Kelvin-Voigt type model (Pipkin 2012) as where τ σ is the relaxation time including the out-of-phase component of conductivity, σ s is the static conductivity, and δ(t) is the Dirac function. Then, the convolution σ * E becomes Hence, we replace two convolutions with multiplications with the effective optical permittivity ε e ∞ and effective optical conductivity σ e ∞ as below: The real effective permittivity and real effective conductivity with respect to the angular frequency ω are given by (Bergmann et al. 1998) Downloaded from https://academic.oup.com/gji/article/232/1/504/6670780 by KIT Library user on 25 October 2022 In fact, the effective optical parameters are the real effective parameters at infinite frequency. The attenuation factor α, phase velocity v and quality factor Q are described as (Turner & Siggins 1994) where the loss tangent tan δ = σ e (ω)/[ωε e (ω)]. With tan δ = 1, we obtain the transition angular frequency ω t = 2π f t = σ e (ω t )/ε e (ω t ) where f t is the transition frequency. In frequency-independent media where τ ε = 0, σ e ∞ can be simplified to σ s . Therefore the attenuation is determined by σ s alone, which is called conductivity attenuation. In frequency-dependent media, however, both σ s and τ ε contribute to attenuation. Based on eqs (7) and (11), we rewrite eq. (3) as We solve eq. (15) using the FDTD method of second order in time and fourth order in space. In the standard staggered grid of Yee (1966), the memory variable r l locates at the same position of the electric field E. The electrical parameters ε e ∞ , ε s , σ e ∞ and τ ε are averaged to the location of E, while the magnetic parameter μ is averaged to the location of the magnetic field H. The convolutional perfectly matched layer is included at the model boundary to absorb the outgoing waves (Roden & Gedney 2000).
Unlike eqs (6)-(8) in Papadopoulos & Rekanos (2011), eq. (15) uses the effective optical conductivity, resulting from different definitions of the polarization current density and memory variable (see eq. 8). Our modification in eq. (15) takes three advantages of the τ -method (Blanch et al. 1995). First, we use one τ ε for all Debye models and set τ Dl as a priori, which means less memory usage and computations in numerical modelling of EM waves than using τ Dl and τ El for L Debye models. Second, if σ s σ e ∞ , we can simply approximate a constant Q by Q ≈ 2/τ ε . Thus, the attenuation is similar to a linear function of frequency and the simulator can account for the waveform distortion of EM waves (see Section 3 for an example). Third, we use only one τ ε in eqs (13) and (14) to quantify the magnitude of attenuation and dispersion introduced by complex permittivity, which is intuitive and reduces the number of parameters reconstructed by frequency-dependent GPR FWI in Section 2.2.
For convenience, we express eq. (15) in a matrix-vector formalism with u = H x , H y , H z , E x , E y , E z , r x1 , r y1 , r z1 , ..., r x L , r yL , r zL T , where 0 3 is the 3 × 3 zero matrix and I 3 the 3 × 3 unit matrix. The superscript * is the transpose conjugate operator which is equivalent to transpose operator T for real variables. When deriving the transpose conjugate of A, we use the zero-valued boundary condition, i.e. equations 54 and 55 in Yang et al. (2016).

Inverse problem
The objective function that we use in GPR FWI is where the synthetic data are extracted from the forward wavefield u at the receiver position by the restriction operator R; d = Ru(m) − d obs is the residual between the synthetic data and observed data d obs . The FWI seeks to iteratively reconstruct the model parameters m by minimizing as follows: where the step length λ is calculated by line search (Pica et al. 1990). Considering the different sensitivities of model parameters to the objective function, we use an individual step length for each parameter (Ernst et al. 2007b). The model update direction for the (k + 1)th iteration m k+1 is computed by the conjugate-gradient method, where γ is the scale factor (Polak & Ribiere 1969) and P is the pre-conditioner (Plessix & Mulder 2004). We use a multiscale strategy to avoid cycle skipping in FWI (Bunks et al. 1995). In the following, we show how to calculate the model gradient ∂ /∂m. Due to that M 1 , M 2 and A are self-adjoint, we obtain the adjoint-state equations as follows (Plessix 2006): where the waveform residual R * d is used as the adjoint sources for backpropagation andû is the adjoint wavefield. In contrast to the forward problem, the computation of eq. (21) is done backward from the recording time T to 0. Hence we reverse the time by substituting t = T − t, u (t ) =û(T − t) and d (t ) = d(T − t) in the above equation and get the self-adjoint-state equations as below: Then we cross-correlate the forward wavefield u and the back-propagating wavefieldû to compute the model gradient Eqs (22) and (23) indicate two advantages of our modification in eq. (15). One advantage is that we can use the same forward solver for the forward and back-propagating wavefield simulations without any additional programming works. Another advantage is that all model parameters are distributed at the diagonal positions of M 1 and M 2 , which means that the form of model gradients is very simple (see Appendix B for details).
Theoretically, in eq. (23), m can be any model parameter of eq. (15). However, we are only interested in the electrical parameters m = (ε e ∞ , σ e ∞ , ε s , τ ε ) in GPR FWI. According to eq. (12) and the chain rule, we convert the gradient of the objective function from the In this study, we use the real effective permittivity ε e and real effective conductivity σ e at the reference angular frequency ω 0 as the input and output parameters of GPR FWI, where ω 0 is set as the peak frequency of the source wavelet. Thus we ensure the physical consistency of the results obtained by frequency-dependent GPR FWI and frequency-independent GPR FWI, similar to the correction for the phase velocity in viscoacoustic FWI (Kurzmann et al. 2013). In the following sections, ε e (ω 0 ) and σ e (ω 0 ) will be referred to as the permittivity ε (or relative permittivity ε r ) and the conductivity σ for convenience if not explicitly stated. In frequency-dependent GPR FWI, one can update the static parameters (ε s and σ s ) using the gradients calculated by eq. (24), or update the effective parameters (ε and σ ) by applying the chain rule again based on eq. (13) (see eq. B8 for details). Corresponding amplitude spectra. The spectra are divided by the maximum amplitude of the thick black line. The translucent grey area represents the amplitude spectrum of the source wavelet, and the thin dashed line marks the peak frequency. One relaxation mechanism is used when τ ε = 0.2 (the blue line).

F R E Q U E N C Y D E P E N D E N C E A N A LY S I S
In order to analyse the effect of different model parameters on EM wave propagation, we make a comparison experiment using the 1-D analytical solution of EM wave in homogeneous media (Blanch et al. 1995). As shown in Fig. 1, we discuss three kinds of typical media in this section. All media have the same relative permittivity (ε r = 6) so that the signals at the reference frequency propagate with the same phase velocity. The latter two attenuating media have the same conductivity (σ = 2 mS m −1 ) in order to generate the same attenuation level at the reference frequency. We take the wave propagation in a non-attenuating medium as a reference (the thick black line). One relaxation mechanism is employed in frequency-dependent medium where τ ε = 0.2 and σ s = 0.1 mS m −1 . Thus, according to eq. (13), the permittivity attenuation τ ε dominates the conductivity (95 per cent) and contributes to the permittivity slightly (10 per cent) at the reference frequency. We set the relaxation time of the conductivity τ σ = 0 s in this paper so that we can focus more on the effect of τ ε on the data. The value of τ ε is chosen to approximate a constant Q (= 10) which has been observed in some geological materials (Turner & Siggins 1994).
Compared to the reference medium, the attenuating media lead to a significant decrease in amplitude in Fig. 1(b). We observe different waveform and amplitude decreases in the frequency-dependent medium (the blue line) where τ ε works as a low-pass filter, with more attenuation at frequencies above the peak frequency and less attenuation in the other frequency ranges. This is presented as a shift of the amplitude spectrum towards lower frequencies in Fig. 1(b), corresponding to the waveform deformation in Fig. 1(a). In contrast, the conductivity attenuation medium scales the amplitude of different frequencies equally, without changing the waveform shape.
To find the reasons of these phenomena, we show the phase velocity v, attenuation factor α and quality factor Q with respect to frequency in Figs 2(a)-(c). Those three characteristics are computed from the real effective parameters shown in Figs 2(d) and (e) through eq. (14). The reference medium has infinite Q value, no attenuation, and a constant velocity. Along with the dominant frequency range of the source wavelet (25-80 MHz), the conductivity attenuation medium has a linear increase of the quality factor, almost the constant attenuation level, and similar velocity with the reference medium due to that the transition frequency (≈6 MHz) is much smaller than the peak frequency of the source wavelet (50 MHz). That results in almost identical travel time and amplitude changes of different frequency components in Fig. 1. Unlike the frequency-independent media, the frequency-dependent medium causes both phase velocity and attenuation factor to increase with frequency in the dominant frequency region, which is the reason for distorting the waveform and amplitude spectrum in Fig. 1. Although we use only one relaxation mechanism, the Q approximated is close to the desired value. With more relaxation mechanisms, the Q approximation can be further improved. As shown in Figs 2(d) and (e), dielectric permittivity and electrical conductivity are frequency-independent in the reference medium and conductivity attenuation medium and become frequency-dependent in the permittivity attenuation medium. When the frequency is greater than 25 MHz, the conductivity is proportional to the attenuation factor, and the permittivity is inversely proportional to the square of phase velocity.

I N V E R S I O N O F S Y N T H E T I C DATA
We use several synthetic examples of 2-D models (x-z plane) to analyse the performance of frequency-dependent GPR FWI. Transverse magnetic (TM) mode waves are used in the simulation and inversion of GPR surface recordings. In frequency-independent media, there are wavefield components H x , H z and E y and reconstructed model parameters ε and σ . In frequency-dependent media, additional memory variable r yl and reconstructed model parameters τ σ and τ ε are added with L = 1. The magnetic permeability μ is constant and equal to its value of vacuum. We set τ σ = 0 s based on the analysis in the previous section.  The EM model space shown in Fig. 3 is 10 m × 36 m, and the grid spacing is 0.1 m. The three-layer background parameters are set according to Table 1. As the 2-m-thick air layer above the ground is not updated during the inversion, we do not show it in the figures of this paper. The true model for the synthetic examples consists of the background model and triangular perturbations, where the permittivity perturbations are relatively weak, so that we can better study the reconstruction of the conductivity and permittivity attenuation models in the inversion. Receivers are placed on the ground at 0.2 m intervals to record the radargrams of the electric field component with a time sampling of 0.15 ns and a time window of 300 ns. Eighteen transmitters generate the electric field at 2 m intervals on the ground (see red stars in Fig. 3). The source is a shifted Ricker wavelet with a centre frequency of 50 MHz. For each source, the receivers record data with an offset of 1-20 m.  The observed data are simulated in the true model and are the same for the frequency-dependent GPR FWI and frequency-independent GPR FWI (τ ε = 0). The minimum wavelength of the EM waves observed in the soil layer is about 1 m.
We use the background model (Table 1) as the initial model to perform FWI. The model parameters are updated simultaneously in multiparameter inversion examples. Apart from the pre-conditioner shown in eq. (20), we also multiply the gradient by a user-defined taper to mitigate source and receiver artefacts. We use a multiscale strategy with five inversion stages and up to 15 iterations per stage. From the first stage to the fifth stage, we sequentially implement frequency band variations from 5 to 15, 25, 35, 50, and 80 MHz in a Butterworth bandpass filter. At the beginning of each stage, we apply the source-time function (STF) inversion ). If the relative misfit change between the current iteration and the previous second iteration is less than 1 per cent, the inversion will switch to the next stage or stop if it is the last stage.

Noise-free data
We build a true model consisting of the background model and several trench anomalies (triangles) spatially uncorrelated on each model. In Fig. 3, variations in the static parameters lead to anomalies in permittivity and conductivity models. To insulate the permittivity and conductivity models from permittivity attenuation anomalies, we adjust these static parameters at locations where τ ε anomalies exist. The strong crosstalks between ε r , σ , and τ ε occur in the frequency-dependent GPR FWI results. In the conductivity model, we observe original conductivity anomaly and severe footprints generated by the anomalies of permittivity and τ ε . Frequency-independent GPR FWI suffers even more from the crosstalks, which makes the original conductivity anomaly indistinguishable. In frequency-independent GPR FWI, more artefacts appear in the non-anomalous area of the permittivity and conductivity models because of the absence of permittivity attenuation. Although permittivity attenuation is weakly sensitive to the data, it is helpful for the correct reconstruction of the conductivity anomaly (the last trench) in the frequency-dependent GPR FWI.
To further investigate the performance of two FWIs, we compare the observed data and inverted data in Figs 4(a) and (b). In the observed data, the ground wave travels faster in the first source than in the last source due to the permittivity anomaly on the left side. The reflected wave of the first source propagates along the low permittivity attenuation anomaly (τ ε = 0) and shows a lower amplitude than that of the last source at an offset of 10-20 m, similar to the observation in Fig. 1(b). In GPR FWI, the STF inversion can provide an "effective" source wavelet that is low-pass filtered (Belina et al. 2012). Thus it can account for the frequency-dependent effects of EM waves in weak permittivity attenuation environments (Q ≥ 20 and τ ε ≤ 0.1). Nevertheless, after travelling through the high permittivity attenuation regions, such an 'effective' source may generate waveforms different from the observed data, e.g., the reflected wave with offsets greater than 10 m in Fig. 4(b). As a consequence, frequency-independent GPR FWI attributes the waveform differences to perturbations throughout the model space, resulting in some unwanted artefacts in Fig. 3. On the contrary, frequency-dependent GPR FWI can reconstruct trench anomalies in the permittivity and conductivity models and reduce artefacts due to its ability to describe permittivity attenuation. Frequency-dependent GPR FWI, therefore, agrees highly with the observed data. This example illustrates that the combination of frequency-independent GPR FWI and STF estimates can only partially explain the waveform distortions caused by permittivity decay. It is necessary to consider frequency-dependent GPR FWI when reconstructing the conductivity model.

Noise-contaminated data
In the synthetic example above, the observed data are free of noise. It is unrealistic and may lead us to overestimate the performance of the GPR FWI approaches. Therefore, starting from this section, we add Gaussian noise with a signal-to-noise ratio (SNR) of 20 dB, with respect to the strongest amplitude, to the observed data. As a result, the reflected waves in the observed data shown in Figs 4(c) and (d) are heavily disturbed and below the noise level after a 15 m offset. We repeat the synthetic example using noisy data and present the inversion results in Fig. 5.
Compared to Fig. 3, we see fewer model updates of the two FWIs in Fig. 5. It is reasonable as noise slows down the convergence of data misfits. On the one hand, the reconstructions of the four trenches more or less deteriorate. On the other hand, the crosstalks and artefacts are also suppressed. Consequentially, frequency-independent GPR FWI shows clearer results of the conductivity model where the original conductivity anomaly becomes distinguishable. Fig. 4 also indicates that the existence of noise increases the difficulty of data fitting. In the last radargram, the two FWIs using noisy data have fewer differences in the waveform than those using noise-free data. To some extent, noise stabilizes frequency-independent GPR FWI and allows it to converge to better results than the noise-free case. Overall, frequencydependent GPR FWI still possesses fewer data misfits and model artefacts than frequency-independent GPR FWI when the observed data are contaminated by noise.

Three-parameter model
In the true model of this example, the perturbations of the permittivity and conductivity models are located at the same position and have the same values in all trenches. However, the static permittivity and static conductivity in each trench are different due to the variations of the permittivity attenuation model. In Fig. 6, we see that frequency-dependent GPR FWI and frequency-independent GPR FWI have similar Downloaded from https://academic.oup.com/gji/article/232/1/504/6670780 by KIT Library user on 25 October 2022 Figure 5. Models of the three-parameter synthetic example where all anomalies are spatially uncorrelated. Gaussian noise (SNR = 20 dB) is added to the observed data. The three columns are the true models, the frequency-dependent GPR FWI results and the frequency-independent GPR FWI results (ε r and σ ), respectively. Figure 6. Models of the three-parameter synthetic example where all anomalies are spatially correlated. Gaussian noise (SNR = 20 dB) is added to the observed data. The three columns are the true models, the frequency-dependent GPR FWI results and the frequency-independent GPR FWI results (ε r and σ ), respectively. performance in terms of reconstructing the permittivity model. The estimation of the conductivity perturbation deteriorates from left to right as permittivity attenuation increases. Besides, frequency-dependent GPR FWI shows more in-trench reconstructions than frequency-independent GPR FWI. In the frequency-independent GPR FWI results, the upper boundary of the last conductivity trench becomes discontinuous and only the two sides can be distinguished. On the contrary, frequency-dependent GPR FWI still reconstructs the four conductivity trenches with good resolution, although the first two trenches in the permittivity attenuation model are incorrectly estimated.

Two-parameter model
We perform another two tests to investigate the crosstalk of different parameters in the GPR FWI. Fig. 7 shows the same true models of ε and τ ε as Fig. 6. The static conductivity model is adjusted so that the conductivity model is unchanged. The increasing values of τ ε from the left to the right trench means that the percentage of τ ε in permittivity increases from 0 to 10 per cent. We observe heavy crosstalks between τ ε and ε in the results of frequency-dependent GPR FWI. The four permittivity attenuation trenches are described as high-value anomalies, and the four permittivity trenches are recovered to different degrees. For example, the last trench is reconstructed better than the first one, even Figure 7. Models of the two-parameter synthetic example where all anomalies are caused by the static permittivity ε s and permittivity attenuation τ ε . Gaussian noise (SNR = 20 dB) is added to the observed data. The three columns represent the true models, the frequency-dependent GPR FWI results and the frequency-independent GPR FWI result (ε r ), respectively. Gaussian noise (SNR = 20 dB) is added to the observed data. The three columns represent the true models, the frequency-dependent GPR FWI results and the frequency-independent GPR FWI result (σ ), respectively. though they have the same value in the true model. This should attribute to the more realistic dispersion and attenuation in the last trench. These differences also occur in the permittivity result of frequency-independent GPR FWI. Since the small effect of τ ε on the phase velocity (see Fig. 2c), the permittivity results obtained from frequency-dependent GPR FWI has only a slight improvement compared to that from frequency-independent GPR FWI (Fig. 7).
When it comes to the synthetic example of σ and τ ε shown in Fig. 8, we adapt the static permittivity to make permittivity the same as the background values. Although four conductivity anomalies look the same in the true model, they are combinations of different static conductivity and permittivity attenuation. The percentage of τ ε in these conductivity anomalies increase from 0 in the first one to 75 per cent in the last one. We observe severe interference of conductivity on permittivity attenuation in the results of frequency-dependent GPR FWI, where all τ ε anomalies are reconstructed as high values. Similar to that shown in Fig. 6, it is difficult for frequency-independent GPR FWI to recover the conductivity anomalies dominated by permittivity attenuation. In the permittivity attenuation model reconstructed by frequency-dependent GPR FWI, the conductivity model introduces crosstalks that become weaker with depth (Fig. 8), while the permittivity model causes crosstalks around the trench boundaries and soil-rock interface (Fig. 7). It is due to different sensitivity kernels of permittivity and conductivity in FWI (Meles et al. 2011). Therefore, in Fig. 6, we observe a superposition of these effects when reconstructing three parameters simultaneously.

Data acquisition and pre-processing
We conducted field measurements at the northeast corner of the gliding airfield in Rheinstetten, Germany. This test site is well known from previous GPR and seismic studies (Schaneng 2017;Wegscheider 2017;Pan et al. 2018Pan et al. , 2021Wittkamp et al. 2019;Gao et al. 2020;Irnaka Downloaded from https://academic.oup.com/gji/article/232/1/504/6670780 by KIT Library user on 25 October 2022 Figure 9. Models of the field data example in the Rheinstetten test site. The three columns are the initial models, the frequency-dependent GPR FWI results and the frequency-independent GPR FWI results, respectively. The white dashed triangle in the initial model outlines the target trench, known as the Ettlinger Line. The red stars are the transmitters, and the dense black triangles are the receivers of the 16th transmitter.  al. 2022). It is covered by sediments consisting of gravel and sand from the Rhine river. The ground layer is composed of partially saturated soil. At the test site, a defensive "V" shaped trench named Ettlinger Line (dashed triangle in Fig. 9) was built in the early 17th century and has been refilled and is now completely flattened to the surface. The existence and shape of the Ettlinger Line have been delineated in detail via 3-D GPR migration imaging (Wegscheider 2017), 3-D Rayleigh-wave dispersion inversion (Schaneng 2017;Pan et al. 2018), 2-D joint elastic FWI of Rayleigh and Love waves (Wittkamp et al. 2019), 2-D viscoelastic FWI of Rayleigh waves (Gao et al. 2020), and 3-D viscoelastic FWI of the surface waves Irnaka et al. 2022). The surface GPR profile is positioned perpendicular to the Ettlinger Line. The two ends of the profile (from southwest to northeast) are in the same location as 'C' and 'B' in fig. 1(a) in Pan et al. (2018). The acquisition settings for the surface GPR data are listed in Table 2. Our GPR data were recorded using a single-channel pulseEKKO Pro GPR system equipped with a pulseEKKO Ultra receiver. The Ultra receiver stacked the records 256 times to obtain a higher SNR by reducing the random noise. The nominal centre frequency of the transmitter is 200 MHz. We deployed the transmitter-receiver orientation in HH mode to acquire TM wave data. The receiver was mounted on a sledge for smooth movement and tracked at the centimetre level by employing a real-time kinematic (RTK) positioning using a self-tracking total station as presented by Boniger & Tronicke (2010). To obtain multioffset GPR data for one source, we fixed the transmitter and moved the receiver towards the transmitter. Then we changed the transmitter location and moved the receiver away from the transmitter to produce the next gather. It took us two minutes to record one radargram and six hours for all 165 radargrams. Our measurements were slower than those of Lavoué (2014) who extracted hundreds of multioffset gathers from 15 common-offset gathers with an offset interval of 0.5 m. However, each trace in the multioffset GPR data used in Lavoué (2014) might be generated by a different transmitter-ground coupling. In contrast, our measurement had the same transmitter-ground coupling in one gather and had denser receiver intervals (∼0.1 m).
In order to apply the FWI, we have to pre-process the GPR data first. The steps for data pre-processing are listed in Table 3, similar to those used in Domenzain et al. (2021). Due to the uneven walking speed of the worker when moving the sledge, the data we acquired has irregular trace spacing. To ensure a balanced illumination in the measurement area, we apply the data gridding, that is 2-D spline interpolation in the time-offset domain at regular trace spacing. Finally, we transform the data acquired in the 3-D world into 2-D line-source data because Table 3. Surface GPR data (200 MHz) pre-processing steps.
(1) Data resampling in the frequency domain (2) Interpolation of clipped direct-arrival amplitudes (3) DC-shift removal and dewow (4) Bandpass filtering (5-400 MHz) (5) Bad traces removal and offset limitation (6) Data gridding in the time-offset domain (7) 3-D / 2-D transformation we use a 2-D forward solver in GPR FWI. We use the transformation of the reflected waves to correct for phase and amplitude differences between the 3-D and 2-D data .

Data inversion
We build an initial model (7 m × 45.2 m with a grid spacing of 0.04 m) in Fig. 9, which does not show the air layer of 1 m thickness. The topographical variations along the survey line are minor and therefore ignored. The initial relative permittivity is 9 at the ground (the velocity of ground wave is 0.1 m ns −1 ) and decreases slightly to 8 at a depth of 6 m. On the other hand, the initial conductivity is 6 mS m −1 at the ground and decreases gradually to 2 mS m −1 at a depth of 6 m. The initial permittivity attenuation is 0.1. We use similar inversion settings as the synthetic examples presented in the previous section, except that the frequency bands vary from 5 to 30, 40, 50, 70, and 100 MHz. These bands are chosen to avoid cycle skipping in GPR FWI since our initial model is relatively simple (Bunks et al. 1995). In frequency-dependent GPR FWI, the relaxation frequency of the Debye model f l is 50 MHz and the reference angular frequency ω 0 = 2π f l , ensuring that the conclusions we drew in Sections 3 and 4 are still applicable to the field data example. To save computational time, we select 18 gathers for FWI with a source interval of 2 m (see Table 2). The reconstructed permittivity models of the two FWIs in Fig. 9 illustrate the presence of a triangular anomaly, the Ettlinger Line, which is in high agreement with the 3-D GPR migration results of Wegscheider (2017). On the right of the trench, we observe a strong permittivity contrast (a slightly right-tilted reflector) at 1.1-1.5 m depth, described as a low S-wave anomaly in the 3-D shallow-seismic FWI results of Pan et al. (2021) and Irnaka et al. (2022). Unlike the shallow-seismic FWI with a resolution of 1 m, GPR FWI has a much higher resolution (∼0.25 m). Consequently, the reflector appears to be connected to the trench in the GPR FWI results, whereas their connection is more ambiguous in the shallow-seismic FWI results. One cannot observe the same reflector on the left side of the trench because the ground surface on this side was higher than on the right side and has been levelled [see schematics in fig. 3 in Irnaka et al. (2022)]. This implies that the reflector might be the original ground surface. Besides, both FWIs reveal a conductive layer near the ground and ranging from a thickness of 1.0 m on the left to 1.5 m on the right. More importantly, frequency-dependent GPR FWI reconstructs a permittivity attenuation model consistent with the permittivity model, where the Ettlinger Line is also visible.
The main differences between the two FWIs come from the right part of the highly conductive layer. This part is discontinuous and more similar to the permittivity model in frequency-independent GPR FWI. To understand these differences, we show the radargrams of the 16th source and the estimated source signals in Fig. 10. For both FWIs, the air and ground waves are difficult to fit because we use a 3-D / 2-D transformation of the reflected wave, and the line source cannot describe the radiation patterns and antenna-ground coupling in the real world. In the 16th radargram, the reflected waves become dominant at offsets greater than 2 m due to the strong permittivity contrast on the right side of the trench. Frequency-independent GPR FWI matches well to the reflection events with offsets shorter than 4 m, beyond 4 m its performance degrades. Frequency-dependent GPR FWI better fits the amplitude of reflection events at offsets greater than 4 m. It results from the high permittivity attenuation layer on the right of the trench, similar to the case in Fig. 4. The two FWIs have similar final data misfits (difference less than 1 per cent) due to strong decay characteristics of the surface recordings along the offset direction. Nevertheless, from the perspective of reflection fitting, frequency-dependent GPR FWI performs better and thus probably produces a more reliable conductivity model. Note that we acquired the first 11 radargrams on the first day and the last 7 radargrams on the second day. The estimated source signals of the first 11 sources are different from the last 7 sources in Fig. 10(b), probably due to coupling differences caused by slight near-surface moisture variations resulting from drying or wetting at night. However, the source signals used on the same day show similar waveform shapes and travel times, which indicates the stability of the two FWIs and source wavelet estimation.

D I S C U S S I O N
In the examples of frequency dependence analysis and inversion, we ignore the relaxation time of conductivity (τ σ = 0 s) and use only one relaxation mechanism (L = 1) for costing purposes. The importance of considering τ σ > 0 s and L > 1 deserves further study. In the field example, the reconstructed permittivity attenuation model delineates the subsurface targets and provides a meaningful geological interpretation. Its high heterogeneous distribution also demonstrates the necessity of including it in the inversion. Based on our observation, we regard it still to be challenging to accurately estimate the permittivity attenuation model due to its weak sensitivity to the data and the crosstalk between multiple parameters. To better image the permittivity attenuation, two strategies are available. One strategy is to decouple the radiation patterns of the reconstructed parameters. It requires choosing a suitable parametrization so that different gradients have opposite behaviour with respect to azimuth (Yao et al. 2018). Another strategy is to use the Hessian operator computed by Newton's method in optimization (Métivier et al. 2013;Operto et al. 2013;Gao et al. 2021). Truncated Newton method has been employed in frequency domain GPR FWI by Pinard et al. (2015). When it is introduced to the time domain GPR FWI, we have to find a balance between the elimination of crosstalk and the increase in computational cost, which needs further investigation.
Previous studies have shown the potential of joint inversion because other data can provide some complementary information to the GPR inversion (Linde et al. 2008;Domenzain et al. 2020;Qin et al. 2022). For example, the electrical resistivity (ER) data are sensitive to electrical conductivity. Therefore, the joint inversion of GPR and ER data may help to improve the reconstruction of conductivity (Domenzain et al. 2021). Note that the estimated conductivity is frequency-dependent in our GPR FWI, but frequency-independent in the ER inversion. Transformation of the two can be achieved by eq. (13), but depends on the reliability of the permittivity attenuation reconstruction, which requires further improvement in the future.

C O N C L U S I O N
In this paper, we quantified permittivity attenuation for EM wave propagation simulator using the τ -method introduced from the seismic community, which not only saves memory and reduces computations in forward modelling but also simplifies the constant Q approximation and frequency-dependent GPR FWI. By defining the parameter τ ε , we proposed a new form of time-domain Maxwell's equations to describe the propagation of EM waves in frequency-dependent media. These equations have two advantages in GPR FWI. The first one is their self-adjoint property, which allows performing frequency-dependent GPR FWI without changing the backpropagation engine. The second advantage is that the model gradients derived from these equations are fairly simple, ensuring as few modifications as possible to the existing frequency-independent GPR FWI code. Frequency dependence analysis revealed that, in the GPR spectrum range, the attenuation caused by the static conductivity is frequency-independent, while that caused by τ ε is frequency-dependent. The permittivity attenuation acts as a low-pass filter and leads to waveform deformation and less amplitude decay of EM waves.
The 2-D synthetic examples confirmed the effectiveness and limitations of frequency-dependent GPR FWI. The spatially uncorrelated examples indicated that surface multioffset GPR data are weakly sensitive to permittivity attenuation compared to permittivity and conductivity. They also demonstrated that both frequency-dependent GPR FWI and frequency-independent GPR FWI suffer heavy crosstalk of different parameters. The spatially correlated examples showed that frequency-dependent GPR FWI works well in estimating the permittivity and conductivity models coupled with different permittivity attenuation levels, and frequency-independent GPR FWI produces comparable permittivity results as the dispersion effect is weak. Although combined with STF estimation, frequency-independent GPR FWI cannot handle the reconstruction of conductivity anomalies dominated by permittivity attenuation. In this case, frequency-dependent GPR FWI can provide more reliable conductivity results even if the permittivity attenuation models fail to be estimated and the observed data contain a degree of noise. These synthetic examples suggested the need for frequency-dependent GPR FWI in the presence of dispersion effects and variable attenuation.
When applying to field data acquired at the Rheinstetten test site, both frequency-dependent GPR FWI and frequency-independent GPR FWI successfully outlined the Ettlinger Line as a triangle permittivity anomaly. On the right of the trench, we found a slightly right-tilted layer which exhibits strong permittivity attenuation. Frequency-dependent GPR FWI showed better data fitting and more continuous conductivity structure in this region, and was likely to be more robust than frequency-independent GPR FWI. Previous GPR migration imaging and shallow-seismic FWIs have verified our GPR FWI results. Crosstalk mitigation and joint inversion with other geophysical approaches need to be studied in the future. The use of other models describing relaxation phenomena of permittivity and conductivity also deserves further investigation.

A C K N O W L E D G M E N T S
This work is financially supported by the China Scholarship Council (No. 201806260258). Tan Qin would like to thank Tilman Steinweg and Mark Wienöbst for their help in developing the WAVE-Toolbox. The authors sincerely thank Lars Houpt, Felix Bögelspacher, Leon Merkel, Michael Mayer, Hagen Steger, Roland Helfer, Philipp Koyan and the master students at Geophysical Institute, Karlsruhe Institute of Technology for their help in field data acquisition, and the editor Rene-Edouard Plessix, the reviewer Hansruedi Maurer and another anonymous reviewer for their constructive comments.

DATA AVA I L A B I L I T Y
An open-source software (GPL) package containing the source code, models and data used in this paper is provided in the WAVE-Toolbox (https://github.com/WAVE-Toolbox). 520 T. Qin, T. Bohlen and N. Allroggen  The convolution of the dielectric permittivity ε and the electric field E is Therefore, we define the memory variable r l of the lth mechanism corresponding to E as The convolution of the conductivity σ and E is We sum the two convolutions and get where ε e ∞ = ε s (1 − τ ε ) + σ s τ σ , σ e ∞ = σ s + ε s τ ε  (2011), we present the gradients as the zero-lag cross-correlation of forward wavefield u and adjoint wavefieldû in the following. Note that we actually replaceû(t) in eq.(23) byû (T − t) which is given by the same forward solver. By doing so, the cross-correlation ofû(t) and u(t) becomes the convolution ofû (t) and u(t). Thus, we obtain the gradients of the electrical parameters as follows: where r εσ = T 0 L l=1 τ Dl (τ Dl r εl + r σ l ) dt, r εl = r xl ∂ t r xl +r yl ∂ t r yl +r zl ∂ t r zl , r σ l = r xl r xl +r yl r yl +r zl r zl .
Eq. (B2) can be simplified by the following equation, based on eq.(A2) where τ Dl (τ Dl ∂ t r l + r l ) = − εs τε L E: r xl E x +r yl E y +r zl E z dt.
Then the gradient of ε s and τ ε are rewritten as: r xl E x +r yl E y +r zl E z dt, r xl E x +r yl E y +r zl E z dt, According to eq. (A5) and the chain rule, we convert the gradients of the objective function from the parameter class m = (ε e ∞ , σ e ∞ , ε s , τ ε ) to m = (ε s , σ s , τ ε , τ σ ) where ε s = ε s and τ ε = τ ε .