DEVELOPING REALISTIC FDTD GPR ANTENNA SURROGATES BY MEANS OF PARTICLE SWARM OPTIMIZATION

The antenna is the most important part of a groundpenetrating radar (GPR) system and defines the employed electromagnetic pulse and how it is transferred to the ground. It is crucial to account for these coupling effects in numerical simulations and to implement realistic antenna models, e.g., for fullwaveform inversion (FWI). We present a method of developing and adapting 3D Finite-Difference Time-Domain (FDTD) models of GPR antennas, complete with electric components, dielectric material properties and feed pulse details. We exemplify this with a commercially available, shielded 400MHz GPR antenna, a model of which was set up by fitting synthetic data to an experimental signal of the antenna reflected at a metal plate in air. For this FWI we used a particle-swarm optimization (PSO) algorithm because the fitted parameters show complex individual effects on the GPR waveform. The resulting antenna model is then validated against data measured in air, water, and with a metal plate in the near-field of the antenna. Overall, the synthetic data reproduces the validation data very accurately. Signals of objects placed in the near-field of the antenna and the change of the shape and frequency content of the radiated wavelet with varying subsurface properties are emulated correctly.


DEVELOPING REALISTIC FDTD GPR ANTENNA SURROGATES BY MEANS OF PARTICLE SWARM OPTIMIZATION Sam Stadler, and Jan Igel
Abstract-The antenna is the most important part of a groundpenetrating radar (GPR) system and defines the employed electromagnetic pulse and how it is transferred to the ground. It is crucial to account for these coupling effects in numerical simulations and to implement realistic antenna models, e.g., for fullwaveform inversion (FWI). We present a method of developing and adapting 3D Finite-Difference Time-Domain (FDTD) models of GPR antennas, complete with electric components, dielectric material properties and feed pulse details. We exemplify this with a commercially available, shielded 400 MHz GPR antenna, a model of which was set up by fitting synthetic data to an experimental signal of the antenna reflected at a metal plate in air. For this FWI we used a particle-swarm optimization (PSO) algorithm because the fitted parameters show complex individual effects on the GPR waveform. The resulting antenna model is then validated against data measured in air, water, and with a metal plate in the near-field of the antenna. Overall, the synthetic data reproduces the validation data very accurately. Signals of objects placed in the near-field of the antenna and the change of the shape and frequency content of the radiated wavelet with varying subsurface properties are emulated correctly.
Index Terms-Electromagnetic propagation, Finite difference methods, Geophysics computing, Ground penetrating radar, Particle swarm optimization, Ultra wideband antennas.

I. INTRODUCTION
Ground-penetrating Radar (GPR) is a geophysical measurement tool that has a broad application spectrum in geoscience and engineering with variations of targets, exploration depth and instrument configurations [1]. It is commonly used for high-resolution, near-surface investigations. The GPR method uses electromagnetic (EM) waves, with frequencies ranging from 1 MHz up to 4 GHz, to image structures that are related to changes in the dielectric and electrical properties of the subsurface. This is not limited to geologic materials and can, among others, include asphalt, concrete, ice and wood [2].
Numerical simulations of GPR scenarios have been used in the past to study various aspects of GPR measurement scenarios. This includes, e.g., the optimization of measurement configurations [3] and the processing and analysis strategies for time-lapse GPR data [4]. Furthermore, complex EM wave phenomena were investigated, such as the dispersion characteristics of subsurface GPR waveguides [5], the influence of a complex subsurface on GPR signals such as fine-layered sediment sequences [6], subsurface heterogeneity [7], [8] or dispersive soil properties [9]. Lastly, numerical simulations are Leibniz Institute for Applied Geophysics, Hannover, Germany, Email: Sam.Stadler@leibniz-liag.de also the basis of inversion techniques of GPR data like fullwaveform inversion (FWI) of GPR crosshole data [10]. There are a couple of numerical methods used for the simulation of electromagnetic wave propagation, such as Method of Moments, Finite Elements and Finite-Difference Time-Domain (FDTD). The latter is widely used for GPR due to its flexibility in modelling and operation in the time-domain and thereby its ability to include a broad range of frequencies with one simulation. Also, modern FDTD software packages are able to keep the total size of the models relatively small by means of minimizing unwanted energy reflections from the finite grid borders [11], [12]. As the name implies, the FDTD method uses a finite-difference grid that is staggered in time and space to compute Maxwells equations [13].
When performing numerical studies of scenarios where antenna coupling and near-field effects appear, it is important to include 3D antenna models in the simulations for the following reasons. In GPR scenarios with shallow targets and complex subsurfaces, the majority of the transmitted energy interacts and moves within the near-field and in particular, within the reactive near-field of the antenna [14]. When the subsurface material properties inside the near-field of the antenna change, so do the antenna characteristics, ground coupling and the emitted wavelet. GPR antennas are usually placed on the ground and often targets are situated in the near-field of the antenna. This results in a complex EM wave-field in the immediate vicinity of the antenna. Another effect that can only be emulated by realistic 3D antenna models are multiple reflections of waves bouncing back into the subsurface at the bowties or metallic shield of the antenna. However, using 3D models of complex GPR antennas for numerical simulations has only recently been introduced [15], [16]. This is mainly because, when simulating in 3D and for a uniformly spaced grid, this means adopting a very fine mesh discretization (<3 mm), which is necessary for accurately implementing the geometry and the electrical components [17], [18], resulting in a great computational effort. Until now, there is a lack of numerical models of GPR antennas especially those of midrange and lower frequency (< 1 GHz), which are commonly used in geoscientific and civil engineering applications (e.g. [19]). This is due to finding a compromise between resolution and depth of investigation. The electric properties of most geological materials are affected by the presence of free water and its relaxation characteristics that cause an increased effective conductivity and reduced permittivity at higher frequencies [20]. In the frequency range between 10 MHz to 500 MHz, the dispersion and attenuation of the electromagnetic waves are the lowest. This is the frequency range called GPR plateau, which typically ensures the best GPR performance. Another factor are scattering losses that strongly increase with higher frequencies. An accurate antenna model will enable more accurate simulations for a range of different purposes. For one, any application where the signal source is in the nearfield vicinity of any other medium would profit from using a full antenna model as the near-field effects would implicitly be simulated correctly. Furthermore, the correct evaluation of signal amplitudes would also benefit from the usage of an antenna model, where the correct amplitude ratio between crosstalk and transmitted signal would be inherent. Lastly, applications like FWI where a GPR source signal has to be guessed would no longer have that need as the antenna model would correctly simulate the near-field interaction between the antenna and the surrounding material properties.
Simplified GPR sources have long been used in the past for a diverse number of applications, e.g. inverse scattering problems [21], [22]. Early works [23]- [26] have focused on implementing simple antenna types and/or parts of antenna models which do not yet include all antenna components. Recent works have progressed to using more geometrically/electrically accurate antenna models including, e.g., the shielding and absorber material [16], [27]- [30]. The difficulty with implementing antenna materials in the simulations is to find accurate parameter values that reliably reproduce the measured signals and antenna coupling effects. Generally, little to no information is available for materials used in commercially available antennas. Different approaches have been used to obtain sets of suitable material parameters. [16] used the Taguchi optimization method [31] to optimize the antenna crosstalk, while [30] used a linear/non-linear full waveform inversion to optimize the antenna crosstalk and the signal from an antenna that is directly coupled to a metal plate. In both works, antennas with frequencies above 1 GHz were used.
Computing capabilities of FDTD forward calculation strongly increased in the last years, due to advancements in hardware and software, such as code parallelization and the use of high performance GPUs. This now allows the use of global metaheuristic optimization methods, e.g., particle swarm optimization (PSO), for computationally expensive problems involving 3D GPR antennas. PSO has the advantage of simplicity, fast convergence rate, ability to tackle ill-posed and nonlinear problems, being little sensitive to local minima, not relying on a priori information and its ability to find a global optimized solution [32]- [34]. PSO has been used in the past in various geophysical fields to perform inversion studies of magnetotelluric, DC resistivity and induced polarization (IP) data [35]- [37]. In the realm of wave methods (seismics and GPR), PSO algorithms have been applied to joint GPR and P-wave traveltime inversion, crosshole traveltime tomography, Scholte-wave inversion, or the identification and characterization of buried objects from constant-offset GPR [38]- [45].
In this paper, we present a performant approach to create numerical 3D antenna models. In contrast to previous works, which optimize the antenna crosstalk and signals from anten-nas directly coupled to a metal plate, we optimize the emitted wavelet by means of analyzing the distant reflection at a metal plate. Furthermore, we implement an antenna in the mid frequency range while so far, only antennas with frequencies above 1 GHz have been modelled. Important to keep in mind is that with decreasing frequency and consequently increasing wavelength and antenna size, it becomes more difficult to have full control over the experimental setting and to collect clean data without disturbances from, e.g., side reflections. We developed a numerical model of a shielded commercially available 400 MHz GPR antenna that emits accurate wavelets in different kinds of media. First, the 3D geometry of the main antenna components (bowties, metallic shield casing, absorber foam etc.) is implemented in the FDTD model. Then, the corresponding material properties are found via fitting synthetic data to reference data of a controlled experiment by a PSO optimization algorithm. Finally, we validate this antenna model by comparing simulated and measured traces in different scenarios and subsurface media.

II. METHODS
The GPR measurements that provided the reference data were performed with a shielded 400 MHz nominal center frequency antenna (model 50400S, S.N.: 1536) and a SIR 4000 apparatus, both from Geophysical Survey Systems (GSSI) Inc. Raw data were acquired without applying any filter. Prior to all measurements, the system was switched on for 15 min to warm up and to reduce any temporal drift of the signal. The center frequencies of the antenna while placed in air (ε r = 1) and on water (ε r = 82) were measured to be ≈ 500 MHz and ≈ 300 MHz respectively. This corresponds to wavelengths λ air ≈ 0.6 m and λ water ≈ 0.11 m.
For the numerical simulations of the GPR scenarios, we used the open source FDTD software package gprMax [46]. It is an established open source software package that allows to set up extensive and complex models for GPR problems. It contains features that allow for fast modelling and simulation via tools such as perfectly matched layer (PML) boundaries, which efficiently minimize boundary reflections [47], and the usage of parallel computation on GPUs [48], which significantly speeds up the simulation. We performed the simulations on a NVIDIA RTX 8000 GPU with 48 GB of RAM.
The goal of our optimization problem is to find a set of parameter values that minimize the fitness function, in our case the relative root mean square error (rel. RMSE) between the amplitude values of the measured and simulated GPR traces.
Hereby, "d" represents the measured data, "G" the forward operator, i.e the FDTD calculation, "m" the model parameters and ∥ 2 the L2 norm. To solve the optimization problem and find values for the antenna model parameters we used the PSO method [49]- [51]. It uses a fixed population of particles (called a swarm) that each represent a possible solution to the given problem to iteratively find a global optimal solution. In the PSO algorithm, each particle represents one set of parameter values p, which lie within predetermined boundaries. Every particle maintains the memory of its overall personal best position P best and velocity v i . Within each generation, the particle velocities are updated jointly with their personal best positions P best and the position of the globally best particle in the whole swarm Gbest, with regard to the fitness function.
The new velocities are then used to compute the new optimal positions for each particle. Over time, the particles move towards more promising regions within the parameter space. The velocity update formula is given by (2), and the position update formula by (3): Hereby, the index i denotes the particle, the index k denotes the iteration, r 1 and r 2 are random numbers in [0, 1], w is the inertia weight, and c 1 and c 2 are the "cognitive" and "social" scaling factors, respectively. Setting w, c 1 and c 2 needs to be done in accordance to the problem at hand, largely depending on empirical selection [52]. The choice of these parameters has been made during preliminary studies on their influence on the inversion result. Herein, the default values of the used open-source software package "inspyred" [53] resulted in a good balance between localizing the minimum and a rapid convergence of the swarm, minimizing the total run-time of our problem. The "inertia" was hereby w = 0.5 to enable an equal weight between the past and the current information, and the scaling factors were c 1 = c 2 = 2.1. These values fall within a range of values explored in previous works [54], [55]. The "neighborhood size" was set to the entire population, to ensure the global best particle is considered when updating the velocity [49], [50].

A. Experimental reference data
For the reference measurements in air, a metal plate was positioned 2 m (3.3 λ air ) away from the antenna in order for the crosstalk and reflection to not overlap in time. The antenna was placed upside down on the ground, covered by a 2 m layer of styrofoam and the metal plate centrally positioned on top of the styrofoam (see Fig. 1). This setup ensured that any spurious reflections from the ground or from nearby objects were minimized. Styrofoam was used because it provides an ideal balance between mechanical stability and having no effective influence on the signal. The material is non-magnetic and has a relative permittivity negligibly different to that of air [56]. To keep the simulation times as short as possible during the optimization, the metal plate needs to be small enough to fit inside a small model space, but large enough to produce a strong enough reflection. A metal plate of 0.4 m × 0.4 m fulfilled both requirements. A DC-shift removal has been applied to center the amplitudes around zero, but no additional processing has been applied to the reference data. The full time signal, containing the crosstalk and reflected signal is shown in Fig. 2.

B. Design of the antenna model geometry
The numerical 3D model of the 400 MHz antenna (Fig. 3) was created by first implementing the geometry of the real antenna, including the shape of the bowties, printing boards, absorber foam, metallic casing, plastic casing and air gaps between particular materials. The real antenna includes a skid plate underneath the plastic casing. As both materials are made of plastic, we merged the skid plate and the plastic casing, resulting in a 8 mm thick layer of plastic below the antenna bowties. Although the absorbing foam is probably composed of several different layers, we assumed a homogeneous absorber to reduce the number of fitting parameters. The transmission lines and electronic components were replaced with electric resistances on the transmitter and receiver feed edges. The resolution of the FD grid was set to 2 mm which is a trade-off between the need for high resolution of the fine structures of the bowties and the computational capacities needed for a 3D-simulation.
The feed point of the transmitter antenna is implemented by applying a voltage source between two conducting wires, i.e. over a FD grid edge, that connect to the bowties (see Fig. 3D). The receiver point at the receiving antenna functions in the opposite way in that it records the voltage over that edge element. Furthermore, both feed points are assigned an electrical resistance, which is not the intrinsic antenna impedance at the open feed point, but a substitute of the characteristic feed line electrical resistance, balun, transmitter and receiver electronics, etc., which are not included in the model.

C. Optimization of the antenna material parameters
Besides the transmitter and receiver electrical resistances on the feed points, the materials that compose the antenna were each assigned an electric conductivity σ and a relative permittivity ε r in gprMax. We assumed the magnetic permeability of every material to be µ r = 1. We also assumed every material to have σ = 0, except for the absorbing foam. Exact values for the ε r and σ of the absorbing foam, the ε r of the PCB and the ε r of the plastic casing were unknown. Furthermore, it is common for the feed pulses of GPR antennas to be asymmetric [57], which is why we chose to construct a pulse that is made up of two differently steep Gaussian-shaped flanks. These two pulse flanks, as well as the above mentioned six parameter values, were the subject of the PSO optimization.
FDTD simulations of the same scenario as during the experiment were performed within the PSO optimization sequence. The model had the dimensions 0.5 m×0.5 m×2.26 m [x, y, z]. For the simulation, the antenna was placed in air (ε r = 1), 2 m away from the metal plate, which was implemented as a perfect electric conductor (PEC). Air was used as a substitute for the styrofoam in the experimental setup (see Fig. 1), as their relative permittivities are very similar. Perfectly matched layers (PMLs) were placed at each border, with the default width of 10 grid cells, to terminate the inbounding electromagnetic waves. The time discretization is determined by gprMax, in accordance to the Courant-Friedrichs-Lewy (CFL) condition [58] and with respect to the highest ε r in the model. The simulation time was set to 24 ns to ensure that the reflected wavelet and its coda are completely recorded.
The calculation of the rel. RMSE corresponding to each particle, i.e. a particular set of material parameter values, begins with cutting out and isolating the reflected signals in the reference and simulated traces. They are then each resampled at the same rate of 1·10 -3 ns and normalized to their maximum absolute amplitude. The normalization is necessary because the GPR system records non-calibrated voltages, and therefore the measured and simulated traces cannot be compared directly. Lastly, the measured and simulated traces are stepwise shifted on to another in time until a minimum rel. RMSE is reached. The shift has to be performed to compensate for small geometry inaccuracies between the experimental setup and the numerical model. Further, varying the antenna parameters during the optimization causes changes of time-zero, when the signal is emitted by the antenna. The resulting, minimum rel. RMSE value is then assigned to each particle.
The lower bounds of the electrical resistances were set to 0.01 Ω because a value of zero would produce a hardsource in the FDTD algorithm that causes unwanted reflections from inbounding electromagnetic field variations. The upper bounds were set to a high value of 1000 Ω. The bounds of the absorber ε r were set to the natural lower limit of one and an upper limit of four, as the absorber is typically a foam which mostly consists of air and hence has a low permittivity. The electrical conductivity of the absorber was also unknown, and we chose to use a range from 0 S/m to 1 S/m. The latter value corresponds to a very high attenuation of ≈ 337 dB/m, for ε r = 2 with a center frequency of 400 MHz. For a twoway travelpath of approximately 0.4 m within the antenna case, this produces an amplitude attenuation of ≈ 135 dB, which is higher than the dynamic range of the GPR system. The bounds for the ε r of the PCB and plastic casing were also set to the natural lower limit of one, and an upper limit of seven. The latter represents a high value for typical PCB material, e.g. FSR-4, as well as for a common plastic, e.g. high-density polyethylene [59], [60].
The two pulse flanks were each assigned a rise time, which is the time the pulse rises from quasi zero to maximum amplitude (RT left ) or falls from maximum amplitude to zero (RT right ), respectively. As a Gaussian function approaches zero asymptotically, we used a lower threshold of 1.11 % of the maximum value to define the beginning of the rise times, which corresponds to three times the standard deviation of a normal Gaussian function. When t max is the time of maximum feed-pulse amplitude, we get: The two rise times had lower and upper limits of 0.12 ns and 19.5 ns. Low rise times create very steep flanks with a higher frequency content than high rise times. The value 0.12 ns was the lowest possible to implement in the simulation due to the CFL criterion. However, such a pulse rise is much steeper than the one that can typically be realized in the electronics of GPR transmitters. A feed signal with a rise time larger than 2 · 19.5 ns results in the signal energy only barely having entered the simulation at all, and produces strong variations in the crosstalk and transmitted signal. Furthermore, pulses with such long rise times contain very little to no energy around the centre frequency of a 400 MHz antenna and would therefore not be suited to drive the antenna effectively (see Fig. 6). An overview of the parameter bounds are summarized in Table I, together with the final optimized parameter values.
For the optimization of the eight antenna parameters, the PSO algorithm was set up with a population of 40 particles. We found this number to be a good compromise between the total run-time per generation and sampling the eight-dimensional parameter space. A single FDTD forward simulation runs in approx. two minutes, resulting in ≈ 1.5 h. computation time per generation. The parameter space is thereby only sparsely sampled, but, even increasing the number of particles significantly, will only result in a marginally better sampling of the individual parameter dimension, while the calculation time will dramatically increase. To counteract the relatively small number of particles per optimization, we ran 10 optimization processes with different seeds, and analyzed the spread of the results. This will provide a measure of the uncertainty of the parameters of the final antenna model.
We defined the stopping criterion for the optimization process by analyzing the long-term trend. This is because the rel. RMSE of a PSO optimization is not necessarily decreasing steadily, but is prone to temporary fluctuations, e.g., when migrating from a local to a global minimum. Comparing the best rel. RMSE of the swarm of the current generation to that of 10 generations prior proved to be a stable measure of the convergence. We stopped the optimization when the change of the best rel. RMSE in the aforementioned comparison was smaller than 0.01%.

A. Optimization results & analysis
Of the 10 optimization runs launched with different seeds, five converged at rel. RMSE of ca. 10.5 %, while the other five converged at ca. 20 %. The latter would result in signal reflections that show significant differences to the measured curves. For the best five optimization runs, the optimization progression of the rel. RMSE values of the best particle per generation is shown in Fig. 4. All parameters converge within the given search interval, with the exception of the ε r casing, which converges very gradually to the lower limit of one. For the aforementioned five best optimization runs, the simulated reflections of the best antenna models, as well as the reference data, are shown in Fig. 5. Optimization run # 5 contains the particle, i.e. set of antenna model properties with the lowest final rel. RMSE of 10.3%. The other four optimization runs have similarly low, final rel. RMSE values, and their reflected signals are barely visually discernable. The final antenna model produces a reflection from a metal plate in air that very accurately fits to the measured trace in terms of the amplitudes and zero-crossings (Fig. 5). Also the frequency spectrum of the reflected wavelet is very well reproduced by the antenna model. It is noteworthy that for comparison, all data curves in the time domain are shifted to the point of minimum misfit, a process that is repeated for all data curves in the time-domain from hereon out. Table I summarizes the optimized antenna parameters including the value of the best PSO inversion run (optimized value), the bounds of the search range, as well as the spread of the results of the best five inversion runs. These resulted in similarly good results with a misfit between 10.3 % and 10.5 %. The pulse-form that is produced by the final RT left and RT right values is plotted in Fig. 6. Therein also depicted for comparison are the sharpest and smoothest pulses that are produced by the parameter bounds.
In Fig. 7 we show all particles of all 10 optimization runs, thus regardless of the generation, per pair of parameters. This provides a visualisation of the fitness function within the search range of the PSO. One has to note, that this figure shows the projection of an eight-dimensional parameter space on the plane of the depicted pair of two parameters. The data were plotted in order of decreasing rel. RMSE, i.e., the good antenna models are plotted on top and hide less good parameter combinations that lie behind them. We can thereby see the overall distribution of particles with low rel. RMSE in terms of the space that each parameter pair encompasses. Relative RMSE values above 16 % (the maximum is 189 %)   Reflections from metal plates in air, produced by the best antenna models from each of the best five optimization runs (see Fig. 4). The orange curve is the measured, reference data, and the blue curve results from the overall best antenna model. Only the peak frequency of the reference data and the overall best, # 5 (blue) are specified on the right. are plotted in yellow and not differentiated in the colorbar, because such antenna models generated signals that show noticeable differences to the measured signal in terms of amplitude and zero-crossings. Particles with hues of blue exhibit the best fits with little to no visual differences between each other and the reference trace. It is obvious that for some antenna parameters as, e.g., the absorber conductivity and RT left , the minimum of the rel. RMSE is sharp, which means that the fitness function is very sensitive to this parameter. On the other hand, the fitness function is less sensitive to, e.g., the source and receiver electrical resistances or RT right , which both exhibit a rather broad minimum plateau. We can also see some dependencies between individual antenna parameters. For instance, the rel. RMSE is small, if either one or both source or receiver electrical resistance are at around 250 Ω, but the rel. RMSE is high, if both electrical resistances are high. Fig. 7 gives a rather general overview of the sensitivity of the fitness function to the individual parameters. A more detailed insight in the vicinity of the minima is given in Fig. 8, which shows the sensitivity of the rel. RMSE to a single parameter when the other parameters are kept constant at their optimum value. To clarify, the overall 8-dimensional parameter-space may look more complex. The changes to  single parameters are therefore dependent on the remaining parameter values that are taken from the final antenna model. The rel. RMSE range is limited to 100 % and does not show the maximum values, to gain more visible information in the range of "good" fits. We can see that the source and receiver electrical resistances have the lowest impact on the emitted wavelet shape, where the impact of σ absorber and RT right is very strong, with a very sharp and pronounced minimum of the fitness function.

B. Validation of the optimized antenna model
In this section, the final antenna model obtained from the parameter optimization (see Table I) is validated for different scenarios. We show the comparison of simulations and measurements for: the crosstalk and the amplitude ratio between crosstalk and reflection with the antenna operated in air, the antenna operated on water and the signal caused by a metal plate in the near-field of the antenna. A good fit of simulated and experimental data for this broad band of scenarios is a good indication that the antenna model produces accurate results, regardless of the medium in use.

Crosstalk and amplitude ratio reflection/crosstalk
As we only used the signal reflected at a metal plate in air for optimizing the antenna parameters, the antenna crosstalk is an independent property that can be used to validate the model. The peak-to-peak amplitudes and zero-crossings conform very well to the measured data and there are only small signal deviations in the coda, between 3.75 ns and 8 ns (Fig. 9). There is also a good fit in the amplitude spectra of both curves and the observation that the crosstalk of the antenna has a lower frequency than the emitted wavelet is well reproduced by the model. The amplitude ratio of a reflection to that of the antenna crosstalk is a characteristic of the distance of the reflecting object, its radar cross section and the antenna characteristics. The amplitude ratio in the measured trace is ≈ 0.0408, and that of the simulated trace is ≈ 0.0389, which corresponds to a deviation of < 5% and confirms the validity of the antenna model for applications in air.

Antenna operated on water
We tested the performance of the antenna model on water, which is an extreme case as water has the highest permittivity of all materials that have to be expected in the field. We compare the simulated and measured reflection from a metal plate (with the dimensions 0.3 m×0.4 m) placed at 0.35 m depth in water (3.2 λ water ), with the antenna resting on the water surface (see Fig. 10). The metal plate was connected to the antenna via nylon strings and positioned centrally below the antenna. In the simulation, the frequency-dependent complex dielectric properties of water are implemented via a Debye-pole, for which we obtained the temperature-dependent relaxation parameters provided by [61]. For a temperature of 15 • C, this resulted in a static relative permittivity ε 0 = 82.1, the high-frequency relative permittivity ε ∞ = 6.0 and the relaxation time τ = 10.8 · 10 −3 ns. The DC conductivity of the water was measured with a conductometer to be σ = 0.0259 S/m and was directly implemented in the model. In Fig. 11 we can see that the reflection is very well simulated with only small amplitude variations in the coda and almost identical frequency spectrum and peak frequencies. However, the simulated antenna crosstalk does not fit as well and there are significant differences in the peak-to-peak amplitude ratio and the zero crossings. Also the amplitude ratio of crosstalk and reflected wave of the measured trace (≈ 5.85) is not reproduced accurately in the simulation (≈ 2.03). The peak frequency is emulated quite well, even though the simulated trace contains higher frequencies than the measured one.

Metal plate within the near-field of the antenna
The calibration data was obtained with a metal plate at 2 m distance to the antenna and characterizes the emitted wavelet. We performed three additional measurements with the metal plate in the reactive near-field of the antenna, and compared them to simulations with the optimized antenna model (Fig. 12). The measurements were done at distances of 0 m, which corresponds to a capacitive short-circuiting of the antenna, 0.06 m (0.1 λ air ) and 0.13 m (0.22 λ air ). The resulting signals change strongly with varying distance of the metal plate to the antenna. This is due to a direct influence of the plate on the antennas characteristics and interference patterns resulting from the signals bouncing between the plate and the antenna. The simulated and measured signals are very similar, with some small deviations in the peak-topeak amplitude ratio, especially when the antenna is placed directly on the metal plate. At 0 m, there are also some small deviations of the zero-crossings, and the frequency spectrum of the simulated trace has a stronger low-frequency content. At 0.06 m and 0.13 m, the overall shape, the zero-crossings and the frequency spectra, with two notable peaks, fit well. The simulated traces have slightly higher peak frequencies than the measured traces, but generally fit well.

V. DISCUSSION
This work builds upon the previous advances in antenna optimization, in particular by [16] and [30]. There, GPR an-  tennas with high center frequencies were optimized (1.2 GHz and 1.5 GHz) by means of fitting the crosstalk of the antennas in air and the signal of an antenna terminated with a metal plate with a Taguchi and a hybrid linear/non-linear full-waveform inversion. In contrast, we optimized a lower frequency antenna and chose to fit the emitted wavelet rather than the nearfield behaviour, which we used as model validation. In further contrast, we chose to use the global optimization scheme PSO because, much more than the aforementioned methods in [16] and [30], it is better able to tackle and explore the unknown, potentially multi-modal multi-dimensional solution space in our problem. Being insensitive to a lack of a priori information as well as to local minima, the swarm behavior of the method is better suited for exploring the entirety of the solution space and thus has a higher chance of finding a global optimal solution. Despite these advantages, the PSO brings along a higher computational cost, however, current computational capabilities allow the use of such expensive optimization strategies. In addition, the effort for this kind of optimization procedure is only necessary once, as a successfully optimized antenna model can subsequently be used indefinitely for simulating/inverting various GPR scenarios. It shall be mentioned that the usage of other methods, such as hybrid variants of PSO and gradient methods, could be explored in this scenario in the future, and may be able to find good solutions faster. We opted to use a PSO algorithm with multiple seeds due to the unpredictable and non-linear influences of the parameters on the emitted wavelet. The results show that the similarly good, final optimized antenna models (see Fig. 5) do not strongly depend on the seed of the starting model as they almost all converge on similar values (see Fig. 4).
Creating an antenna model by fitting the parameters to calibration data will not necessarily result in an exact 1:1 replica of the real antenna, but yields a model that shows a similar behavior and produces signals very close to that of the real antenna. Therefore, we did not expected that every parameter of the antenna model will converge to realistic values as a model simplification at one part of the antenna can cause a counterbalance in the properties of another part. In the end, our optimization routine yielded an antenna model that very accurately replicates measured signals in various scenarios. The model was validated in air and water, which cover the span of most extreme cases that can be expected in the field. It is able to replicate main aspects such as frequency shifts of the emitted wavelet when the antenna is operated on different media, the amplitudes and zero-crossings of the signals, as well as near-field effects such as interference patterns between crosstalk and reflected signals, multiple reflections and antenna ringing.

Experimental limitations
As the calibration data are the basis of the deduced antenna model, the reproducibility of the experimental data was striven for, but is variable for a number of reasons. Firstly, measurements with different GPR antennas but of the same type, probably produce slightly different signals caused by, e.g., manufacturing tolerances of the electronic components. Secondly, some experimental uncertainties remain, even though much care was taken to accurately set up the geometries, to avoid spurious reflections from the side as well as from improper lain cables and to minimise the time-drift of the antenna and GPR apparatus. In summary, the experimental uncertainty is probably larger than the variance between the five best results of our PSO optimization runs, that show a misfit between 10.3-10.5 % to the calibration data. More accurate reference data from experiments under controlled conditions, e.g., in a box with varying homogeneous materials would be desirable. However, realizing these measurements, with decreasing antenna frequencies and consequently rising wavelength and larger antennas, is costly and poses a big challenge.

Model limitations
All the geometric details we could obtain, were implemented in a numerical model (Fig. 3) with a 2 mm spatial discretization, which could not be lowered due to the FDTD calculation capacities. Consequently, delicate structures as, e.g., the feed points of the bowties, PCB boards or air gaps may not be resolved adequately. An improvement to these limiting circumstances would be to incorporate subgridding routines in the FDTD code in the future to be able to refine the model grid where the implementation of fine geometric details is necessary while keeping the overall model size and simulation run-time at an acceptable level. Overall, not all antenna parts could be replicated completely and had to be either simplified, e.g., the absorber foam, or neglected, e.g., the transmitter and receiver electronics. Notably the antenna electronics were implemented as transmitter and receiver resistances. Integrating the electronics of the antenna into the model and describing the interaction between electronics, transmission lines and the antenna, seems to be the most challenging point and currently, there is no practicable solution in sight.

Optimization performance
Of the 10 optimization runs we started, and with our stopping criterion, five converged at rel. RMSE values of ≈ 10, 5 % (Fig. 4). The "failure-rate" of five out of ten optimizations converging at higher rel. RMSE values of ≈ 20 % is most likely due to the entire parameter space being sparsely sampled, and the complexity of the fitness function. These circumstances increase the probability of the PSO swarms converging in local minima. The simulated reflections of the best antenna models proved to have a very good fit to the measured trace, with very little deviations in the amplitudes, zero-crossings and frequency spectra (Fig. 5).
The aforementioned five optimization runs show slightly different progressions. This was expected, as each optimization started with a different seed and the total eight-dimensional parameter space is large compared to the population size. Of the eight fitted antenna parameters, the absorber conductivity and RT left proved to be well defined in all optimization runs and tended to have a large influence on the wavelet (Fig.s 4  and 7). The remaining six parameter are also found by all optimization runs (Table I), albeit at slightly different stages (Fig.s 4), which means that these parameters have a smaller influence on the wavelet.

Fitted antenna parameters
Even though the absorber was implemented with only one layer, as supposed to several layers in the real antenna, its material components converged on nearly identical values in the five best optimization runs. Figures 4 and 8 both show, that for ε r and especially σ absorber, similar values are found fast, and that they lie within a relatively well defined minimum with regard to the rel. RMSE The two rise-time parameters, of the left and right side of feed pulse, are assigned two very different values, with a sharp rise from the start and a relatively slow decline after the maximum amplitude (Fig. 6). Not only was a difference in the two parameters expected, but the sharp rise and slow fall in amplitude also coincides with the typical form of antenna feed pulses [57]. In all five optimizations, RT left converges on almost the same value while RT right shows some variation ( Table I).
The rel. RMSE as a function of RT right shows a rather broad minimum plateau (Fig. 8) which indicates that it only has a small impact on the wavelet. In the future, the pulse form could be created in a more flexible manner, by increasing the parameterization of the pulse, albeit with a higher parameter number and therewith the expense of a longer optimization duration.
The source electrical resistance, and especially the receiver electrical resistance, are the parameters with the widest range in their final values. This indicates, that these two parameters have a relatively small influence on the shape of the transmitted antenna signal. The absolute amplitude of the recorded signal is also influenced by the electrical resistances. This however cannot be considered in the optimization, as the amplitude is measured in an arbitrary scale by the GPR system and not given in calibrated voltages. Hence, we had to normalize the experimental and simulated data during the optimization.
The permittivity of the plastic casing is the only parameter that converged on one of its bounds, the lower bound of ε r = 1, which is a non-realistic value for plastic. We attribute this to be possibly an effect of geometric discrepancies between the real and numerical antenna due to the limited spatial discretization. These discrepancies include the small air gap between bowties and plastic casing, and the grooves of the skid plate which were not implemented in the model. Fig. 8 also shows, that the rel. RMSE, as a function of the ε r of the plastic casing, does not have a steep descent in the vicinity of the minimum. This indicates, that a slight change in ε r value would only slightly change the performance of the antenna, which would in turn however represent a more realistic ε r . A possible solution for finding a realistic value for the ε r of the casing, could be to weigh the match to a measured signal in a media other than air, into the fitness function. This however will remove the option of an independent validation of the antenna model or require further measurements under controlled conditions, e.g., by using oilwater emulsions as in [16]. The latter seem unfeasible or at least extremely expensive when targeting the emitted wavelet for lower-frequency antennas.

Validation and performance of the antenna model
The validation studies support the accuracy and versatility of the antenna model. In the analyzed scenarios, the simulated traces showed characteristics that are very similar to those of the measured traces ( Fig. 9-12). The shape of the crosstalk in air is almost perfectly reproduced and the amplitude ratio of the reflection to the crosstalk ends up being very close, with a deviation of < 5 %. A possible future improvement could be to include and weigh the amplitude ratio in the fitness function of the optimization. When the antenna is operated on water, the emitted wavelet is also accurately reproduced by the model, including the significant shift in center frequency when operated in water (≈ 300 MHz) and in air (≈ 500 MHz).
Only the shape of the crosstalk of the antenna on water could not be simulated precisely. This could be caused by the unrealistic permittivity of the casing (ε r = 1) of our model or further geometric inaccuracies of the model, as discussed above. Additionally, water flowing between the plastic casing and the skid plate mounted directly below it, could possibly change the crosstalk signal in the experiment. This effect is not reproduced in the simulation, as the skid plate is not separately included in the antenna model. Lastly, the differences could be caused by the feed cables and antenna electronics, which are both not included in the model, but might cause effects when the antenna is operated on material, for which it has not been designed leading to a mismatch of impedances. In the future, to reduce the misfits caused by simulating antennas on subsurface media different to that used in the optimization, a joint optimization on multiple media could be employed. This of course, will require further accurate experimental data and increase the computation time of the optimization.
Despite the aforementioned differences of the crosstalk shape, when the antenna is operated on water, its peak frequency and total length are well reproduced. As the latter defines the time window of the radar trace that is masked, and as the wave emitted in water is emulated correctly, we deem the limitation of the simulated crosstalk acceptable. The near-field effects of the antenna are very well reproduced by the model and show only very small differences. They can be explained by minimal deviations of the geometry, i.e. the distance of the metal plate to the bowties, as the measured signal is the result of a complex interference pattern of crosstalk and reflection in the near-field of the antenna.

Transferability of the approach
Generally speaking, the presented optimization methodology could be applied to antennas of different center frequencies. For antennas with higher center frequencies, the main limitation would be the fine numerical spatial discretization that might be necessary to achieve an adequately accurate numerical antenna replica. This point might become less impairing in the future with the usage of subgridding. On the plus side, higher frequency antennas are smaller in size, resulting in an overall smaller numerical model. For antennas with lower center frequencies, the spatial discretization need not be that fine, however, gaining high quality reference data without any side reflections becomes more difficult.

VI. CONCLUSION
Generating realistic FDTD antenna models can be done by implementing the geometry of the main antenna components and fitting the unknown material properties and feed pulse to experimental data. The approach of using a set of particle swarm optimizations coupled with a FDTD forward calculation proved to be effective. Even though the model does not include all details due to the limited spatial discretization or simplifications as replacing the electronics of the antenna by an electrical resistance, all but one of the fitted parameters were in the range of realistic values. Analyzing the fitness function, i.e., the misfit of the model response, provides an insight into the sensitivity of the signal of the antenna to the single antenna components. The absorber properties and the feed pulse shape were the most influencing factors, whereas the other material properties and the electrical resistance of the electronics proved to be less important. The final optimized model is able to reproduce the main features of the antenna very precisely. This includes frequency shifts of the emitted wavelet when subsurface properties change as well as complex interference patterns, when an object is placed in the near-field of the antenna. Including full antenna models and accounting for the described effects is therefore superior to using simple dipole point sources and is crucial for realistic simulations used for, e.g., full-waveform inversion (FWI). The presented 400 MHz model is a replica of an antenna that is widely used in geoscientific and engineering applications and can directly be used by the GPR-community. Further antenna models can be designed by following the presented optimization approach and incorporating the required unknown parameters of that particular antenna. Besides information on the geometry of the major components of the antenna, it would be helpful if calibration data of antennas on different material were available, e.g. by the manufacturers. These data could form the basis to build up an open-source antenna model database, from which the whole GPR-community will benefit.