Estimation of Handheld Ground-Penetrating Radar Antenna Position with Pendulum-Model-Based Extended Kalman Filter

: Landmines and explosive remnants of war are a signiﬁcant threat in tens of countries and other territories, causing the deaths or injuries of thousands of people every year, even long after military conﬂicts. Effective technical means of remote detecting, localizing, imaging, and identifying mines and other buried explosives are still sought and have a great potential utility. This paper considers a positioning system used as a supporting tool for a handheld ground penetrating radar. Accurate knowledge of the radar antenna position during terrain scanning is necessary to properly localize and visualize the shape of buried objects, which helps in their remote classiﬁcation and makes demining safer. The positioning system proposed in this paper uses ultrawideband radios to measure the distances between stationary beacons and mobile units. The measurements are processed with an extended Kalman ﬁlter based on an innovative dynamics model, derived from the model of a pendulum motion. The results of simulations included in the paper prove that using the proposed pendulum dynamics model ensures a better accuracy than the accuracy obtainable with other typically used dynamics models. It is also demonstrated that our positioning system can estimate the radar antenna position with the accuracy of single centimeters which is required for appropriate imaging of buried objects with the ground penetrating radars.


Introduction
The presence of landmines and explosive remnants of war (ERW), such as artillery shells, grenades, rockets, bombs, and cluster munition remnants, poses a significant worldwide threat in the areas of current and past military conflicts.It results in deaths and injuries of mostly civilian victims even many years after the wars.
According to the yearly reports of the Landmine Monitor [1,2], providing a global overview of the landmine situation, tens of millions of landmines are still buried underground in at least 60 countries and other territories.Only a single year 2021 brought 7073 casualties of mines/ERW (2492 killed and 4561 injured) in 54 different countries, and 80% of the victims were civilians [1][2][3].
Considering the significance of the problem, efficient methods of mine clearance are still tough.Currently, various metal detectors (MD) are often used for this purpose, and contemporary MDs offer excellent parameters, enabling the detection of even very small and deeply buried metal objects [4][5][6][7][8][9].Paradoxically, this high sensitivity can be also their drawback leading to many false detections which lengthen the time necessary for demining.Moreover, MDs do not offer any way to initially identify or classify the detected objects and every detection must be carefully examined.What is even more problematic and dangerous, not all contemporary landmines and ERWs contain metal elements, which limits the usefulness of MDs in mine clearance operations.
In military applications, GPRs can be installed on large armored manned vehicles with enhanced immunity to nearby explosions [31][32][33][34][35].For increased safety of the crew, the radar antennas are usually attached to the end of long arms in front of the vehicle.A good alternative is mounting GPRs on remotely controlled unmanned wheeled vehicles [36] or tracked vehicles [30,35,37,38], which eliminates the risk for the crew, and reduces the costs of the purchase and the exploitation of such systems.The GPRs on vehicle platforms, however, have limited utility in difficult terrain: mountainous areas, forests, dumps, urban surfaces covered with debris, or interiors of buildings, where landmines and other explosives can be typically found.A good solution applicable in such areas is a handheld version of the ground penetrating radar (HH-GPR) [39][40][41][42].The problem of estimating the antenna position of such type of radar is addressed in this paper.
The GPR operation requires emitting electromagnetic energy in the direction of the ground.The transmitted radio waves penetrate near-surface layers of the soil and encounter on their way various objects and layers of different permittivity ε and conductivity σ, which results in reflecting and scattering back a portion of the transmitted energy.The echo signals are received, collected, and processed to detect and create images of buried objects.
Most contemporary GPRs are pulse radars [15,16,22,24,43,44], transmitting repeatable, very short, high-amplitude pulses and receiving strongly attenuated echo signals reflected or scattered back from layers' boundaries and buried objects [3,45,46] as shown in Figure 1.Time delays of subsequent peaks in the received echo signals are proportional to the depth of the detected objects or layers of different permittivity.Collecting and joint processing multiple echo signals, so-called echograms or radargrams, for a GPR moving along a predefined scanning path enables locating and imaging those objects and layers [3,15,16].
Three types of visual presentations of GPR radargrams are used in practice [3,14,16,21].A single echogram was obtained for only one GPR antenna position with coordinates (i, j) is a one-dimensional signal representation, called an A-scan (Figure 2a).Time delays of the signal peaks in the A-scan are usually converted into respective depths and the Z-axis is scaled in the distance units [21,23,39,46,47].An analysis of GPR data is typically based on a two-dimensional signal representation, called a B-scan, which is a dataset created from many A-scans acquired for various antenna locations along a usually linear scanning path, as shown in Figure 2b.It represents a radar image of a vertical surface intersecting the scanned terrain volume below the scanning path.Due to a relatively large GPR antenna beamwidth, the same buried objects are illuminated many times from different antenna locations and consequently from different distances.Therefore, the echo signals form hyperbolic structures visible in the B-scans [14,23,39,46,47].An example of such a structure for a single-point object is shown as a red hyperbole in Figure 2b.
Collecting A-scans for multiple antenna locations in the nodes of a grid span onto the OXY surface, one can create another type of GPR signal visual presentation, called a C-scan (Figure 2c).This is a three-dimensional signal representation, which is very useful in visualizing, identifying, and classifying buried objects.
The C-scans are often presented as a set of two-dimensional greyscale or color images, created as horizontal sections through the C-scan volume on various depths [3,21,48].An example of such a single image is shown in Figure 3.The presented scans were made using a pulse radar produced by IDS GeoRadar company, containing a DAD K2 control unit and an antenna with a central frequency equal to 900 MHz.This radar made 850 soundings per second, the duration of the probing pulse was about 1 ns, and the obtained spatial resolution was about 5 cm.
Knowledge of accurate positions of a GPR antenna moving along a scanning path is necessary to properly assemble all the acquired radargrams and create high-quality GPR B-scans or C-scans.Several scientific papers [44,49] and patents [50] suggest that the GPR antenna positioning accuracy should be better than one eight of a radar signal wavelength [51].As typical GPRs work at a frequency range between 400 MHz and 4 GHz (wavelengths from 7.5 to 75 cm) [49,51], the antenna positioning accuracy should be of the order of single centimeters which requires using very high-accuracy navigation systems.
As most of the listed above devices or systems are not adequate for HH-GPRs, due to their large size, weight, specific installation requirements, vulnerability to jamming or signal shadowing, and too low accuracy, the authors of this paper proposed a system based on several ultrawideband (UWB) radio modules.This concept was first described in an authors' conference paper [54], where physical models of a mobile unit and UWB beacons were presented.The mentioned paper also contained a description of an autocalibration procedure, used for self-locating the UWB beacons for quickly establishing a frame of reference before the scanning process, and presented an initial assessment of the system's accuracy which in the scanning zone reaches desired level of 2-3 cm.
In another authors' conference paper [41], it was claimed and demonstrated that the accuracy of the UWB positioning system can be further improved with a properly chosen estimation algorithm.In that paper, using an extended Kalman filter (EKF) based on a GPR antenna motion model, derived from the mathematical pendulum motion model, was proposed.The mentioned paper, however, contained a proof of concept rather than a complete and applicable positioning solution, as the proposed pendulum-based dynamics model used in the EKF was oversimplified to present the main idea only.It assumed that the attachment point of the "pendulum", which is the position of a GPR operator's arm, is initially known and that the angle of orientation of the main axis of the scanning section always equals zero degrees.These assumptions can hardly be met in practice.Moreover, the mentioned conference paper contained only a sketch of the system's model and very limited results of its simulative testing.
This paper can be considered a significantly extended version of the above-mentioned conference paper.It presents an elaborated, practically applicable version of the GPR antenna positioning system using UWB radio modules and includes a complete description of its extended mathematical model and detailed results of its simulative testing.The main novelty of this paper includes: 1.
Elaboration and detailed presentation of an advanced and practically applicable dynamics and observation model of the UWB-based GPR antenna positioning system, with relinquished simplifying assumptions of the model presented in [41]; 2.
Elaboration and detailed presentation of the estimation algorithm used in the proposed GPR antenna positioning system; 3.
Presentation of new and detailed results of simulative tests of the positioning system for various realistic system configurations.
This paper is organized as follows.A general concept of the ground penetrating radar, types of GPR data visualizations, accuracy requirements for GPR antenna positioning, technologies used for GPR positioning, previous authors' works in this field, and a discussion of the novelty of this paper are presented in Section 1.The system's description, its mathematical model, and the estimation algorithm elaborated by the authors are presented in Section 2. The methodology and the results of simulative testing of the GPR antenna positioning system are presented in Sections 3 and 4 contain a discussion.

Scanning Profiles
As has already been mentioned, creating B-scans requires moving a GPR antenna over the ground, ideally along a linear scanning path (profile) with constant velocity, to collect linearly arranged and uniformly separated radargrams.Creating C-scans requires repeating such scanning (profiling) for many equidistant lines in one direction, as shown in Figure 4a, or bi-directionally, as shown in Figure 4b, where the antenna position is marked as a letter A [7,15,16].In multichannel GPRs, with several equidistant antennas, the profiling can be realized quicker, unidirectionally (like in Figure 4a) for several scanning paths at a time.Although in favorable conditions the profiling shown in Figure 4 can be at least approximately realized with GPRs installed on vehicles (carefully driven or remotely controlled, in non-demanding terrain and with the use of an accurate supporting navigation system), this can hardly be achieved with HH-GPRs.The elements of the scanning path, in this case, are shown in Figure 5, where the letters A and S represent the positions of the antenna and the sapper.

UWB Positioning System
The structure of the HH-GPR antenna positioning system proposed in this paper is shown in Figure 6.It is composed of four stationary modules M 1 ÷ M 4 serving as radio beacons and two mobile modules M A and M S .The M A module is installed over the GPR antenna and the M S module over the sapper's shoulder.All the modules contain UWB transceivers.Distance measurements realized by these transceivers are collected and processed using estimation algorithms described in the further part of the paper.The following variables are used in Figure 6: d Aj -distance between a j-th beacon and the antenna module M A , d Sj -distance between a j-th beacon and the sapper module M S , x j , y j -coordinates of a j-th beacon position, x A , y A -coordinates of the M A module position, x S , y S -coordinates of the M S module position, l-length of the HH-GPR handle (horizontal distance between M S and M A ), θ-angle between the horizontal projection of the GPR antenna handle and the central axis of the scanning section.
We assumed that the UWB radios used in our system are PulsON P440 modules from TDSR [55].They use the two-way time-of-flight (TW-TOF) method for ranging and offer an operating range between 300 and 1100 m and a ranging accuracy of about 2 cm in line of sight (LOS) conditions.Such parameters give the potential to build a positioning system with the desired centimeter-level accuracy, required in the considered application of the HH-GPR antenna positioning.
The placement of beacons outside a potentially hazardous area, as shown in Figure 6, is only one of the possible options, suggested for quick and easy deployment of the system in terrain.Other beacons' locations are also possible, and their relative positions with respect to the mobile units M A and M S influence the accuracy of the UWB positioning system, which will be discussed in detail in the Results section of the paper.

Mathematical Model
As can be seen in Figure 5, the scanning profiles are composed of fragments that resemble arcs rather than straight sections.Moreover, the velocity of the HH-GPR antenna is more changeable than in GPRs installed on vehicle platforms, as typically an operator (sapper) performs a swinging motion, initially accelerating and finally decelerating the antenna.Therefore, the collected radargrams are not linearly arranged nor uniformly separated.Nevertheless, the acquired A-scans can be used to create two-or three-dimensional GPR visualizations of buried objects provided that the antenna positions are known for all the collected radargrams [3,15,16].
A single arc belonging to the scanning profile is shown in Figure 7.If we consider the changeable angular velocity of the antenna motion (initially accelerating and finally decelerating), such a trajectory resembles the motion of a mathematical pendulum [56], and can be described by the following formula: where: θ-angle between the horizontal projection of the GPR antenna handle and the central axis of the scanning section, a-acceleration forcing the HH-GPR antenna (M A module) motion, l-length of the HH-GPR handle (horizontal distance between M S and M A ).The acceleration a is analogous to the gravity acceleration g in the mathematical pendulum motion model.Contrary to g, which can be considered a constant, the acceleration a is more changeable and to large extent depends on the operator's strength, fatigue, style of HH-GPR operation, etc., thus we treat it as an additional variable to be estimated and we model it as a Wiener stochastic process [57][58][59].Considering the geometrical relationships shown in Figure 7, the equations describing the antenna and the sapper's arm motion can be formulated as follows: where: x A , y A -coordinates of the HH-GPR antenna (M A module) position, x S , y S -coordinates of the sapper's arm (M S module) position, u x S , u y S -Gaussian white noises representing random components of the sapper's arm (M S module) motion, l-length of the HH-GPR handle (horizontal distance between M S and M A ), θ-angle between the horizontal projection of the GPR antenna handle and the central axis of the scanning section, γ-angle between the horizontal projection of the GPR antenna handle and the OY axis of the frame of reference, ω-angular velocity of the HH-GPR antenna (M A module) motion, a-acceleration forcing the HH-GPR antenna (M A module) motion, u a -Gaussian white noise representing random changes of a.
Rewriting Equation (2) to fit it into the standard form of a nonlinear continuous dynamics model [60][61][62][63]: one obtains the following detailed version of this model, which has been further used in our estimation algorithm: The nonlinear observation model in the following standard form [60][61][62]: has been formulated assuming that at every step k the UWB positioning system realizes four distance measurements between a j-th beacon and the antenna module M A : and four distance measurements between a j-th beacon and the sapper module M S : where: d Aj -distance between a j-th beacon and the antenna module M A , d Sj -distance between a j-th beacon and the sapper module M S , x j , y j -coordinates of a j-th beacon position, x A , y A -coordinates of the M A module position, x S , y S -coordinates of the M S module position, h-sapper's arm height, v Aj , v Sj -distance measuring errors for M A and M S modules.
A detailed version of the observation model, which has been further used in our estimation algorithm, is as follows: As the antenna module M A is kept close to the soil during scanning, and the differences between slant distances d Aj and their horizontal projections are very small, we assumed that their altitude over the ground can be omitted in the observation model.On the other hand, the M S module is placed over the ground on the sapper's arm, and its altitude h is non-negligible.In our model, we assumed that it is constant, as its changes in the order of centimeters during the system's operation can be neglected for typical distances from the UWB beacons, which are in the order of tens of meters.In a real system the altitude h can be a settable constant, adjusted before using the system, based on the sapper's height.

Estimation Algorithm
An extended Kalman filter for HH-GPR antenna position estimation was designed based on the previously described dynamics and observation models and its flowchart is shown in Figure 8.
After initialization of the EKF at step k = 0 or after closing each subsequent filter's loop at steps k > 0, the filter alternately performs prediction and correction steps.The prediction step requires previous calculations of the fundamental matrix F, the transition matrix Φ, and the covariance matrix of disturbances Q at every step k.The method of calculating the F matrix (more precisely it is F k−1 but to shorten the notation the index k − 1 will be omitted in further equations) is explained in Appendix A. where: and ∆t is a period between two successive time steps k − 1 and k.
The Q c matrix from Equation (11) represents the covariance matrix of continuous disturbances which is a 3-by-3 diagonal matrix containing power spectral densities S x S , S y S and S a of the noises u x S , u y S and u a composing the disturbances vector u(t) in Equation ( 4): The predicted state vector xk|k−1 is calculated in accordance with the following general equation [64,65]: but in practical calculations we use Heun's numerical integration method [65][66][67][68] which leads to the following formulae: where xk−1|k−1 is the final state vector estimate from the previous step k − 1.
Apart from the predicted state vector xk|k−1 , the covariance matrix of prediction errors P k|k−1 is calculated based on the covariance matrix of filtration errors P k−1|k−1 from the previous step k − 1 as follows [60][61][62]: where we use the mentioned matrices Φ and Q.The correction step requires previous calculations of the observation matrix H at every step k, and the method of its calculation is explained in Appendix B. This step involves a calculation of the Kalman gains matrix K k , a correction of the predicted state vector xk|k−1 based on the current measurement vector z k , which produces the final estimate xk|k at the step k, as well as calculation the covariance matrix of filtration errors P k|k , and these operations are realized as follows [60][61][62]: The R k matrix in Equation ( 16) is the covariance matrix of measurement errors [60,61] which is formed as an 8-by-8 diagonal matrix, containing the variances of all eight distance measurements performed between pairs of UWB modules in the positioning system: where σ 2 Aj and σ 2 Sj represent variances of distance measurement between a j-th beacon and the antenna module M A or the sapper module M S .

Alternative Positioning Algorithms
Apart from the proposed pendulum-model-based EKF, simpler algorithms can also be used to estimate the HH-GPR antenna position.One possible solution is a non-linear least squares (NLS) algorithm [57,69,70] which processes a vector z(k) of distance measurements collected at each step k without using the previously estimated state vector and without filtration.Such an algorithm does not use any dynamics model either.The NLS requires an initialization by assigning at least coarse values to the antenna coordinates x A and y A and subsequently, it improves their estimates iteratively.This algorithm is simple but due to lack of filtration, its accuracy is not high.
Better estimation results can be obtained by using EKF filters based on nearly-constantvelocity (CV) or nearly-constant-acceleration (CA) dynamics models [57,[71][72][73][74], which are typically applied in navigation and radiolocation.The CV model (20) assumes a rectilinear uniform motion, whereas the CA model ( 21) assumes a uniformly accelerated motion, and, in both cases, small disturbances of these ideal movements are modeled by the vector u(t): 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 where: x A , y A -coordinates of antenna position, v x , v y -components of antenna velocity, a x , a y -components of antenna acceleration, u v x , u v y -Gaussian white noises representing random disturbances of CV motion, u a x , u a y -Gaussian white noises representing random disturbances of CA motion.
Clearly, the CV and CA models do not fit ideally the actual HH-GPR antenna motion but nevertheless, they can be used for prediction in EKFs.Such filters are not optimal, but they are simpler than the EKF presented in Figure 8 because both dynamics models are linear, and the prediction of the state vector is realized like in a linear Kalman filter [60][61][62]: Moreover, the transition matrix Φ and the covariance matrix of disturbances Q can be calculated in advance before the filter implementation using the simple formula [62] and they remain constant during the filters' operation.Thus, such EKFs do not require in-run calculation of the Jacobian matrix F and the matrices Φ and Q.
All the mentioned algorithms, NLS and EKFs based on the CV and CA models, have been implemented by the authors and tested to compare them with the previously described pendulum-model-based EKF.Further in the paper, the following acronyms will be used for these algorithms: NLS, EKF-CV, EKF-CA, and EKF-PND.
Although the EKF-CV and the EKF-CA are not optimal, their accuracy can be maximized by choosing appropriate power spectral densities of disturbances S v x , S v y of u v x , u v y noises in the CV model or S a x , S a y of u a x , u a y noises in the CA model.Their choice affects the values of the Q matrices and consequently influences the information quality [75], notably estimation errors of the filters.The process of choosing filters' parameters and minimizing their estimation errors is called "tuning Kalman filter" [76] and it was realized in the case of the EKF-CV and EKF-CA designed in our research.

Results
The described HH-GPR antenna positioning system and the proposed pendulummodel-based EKF were implemented and simulatively tested in MATLAB ® version R2022b.The simulations included an assessment of the dependence of the system's accuracy on the positions of UWB stationary modules M 1 ÷ M 4 deployed around the area of interest, where the demining process is going to be performed.The results of these analyses are presented in Section 3.1.In further experiments, the accuracy of the EKF-PND filter was analyzed for chosen scanning sections.This accuracy was also compared with the accuracy of the NLS, EKF-CV, and EKF-CA algorithms.The results of these tests and accuracy comparisons are presented in Section 3.2.

Influence of UWB Beacons' Locations on System's Accuracy
Possible locations of the UWB stationary modules M 1 ÷ M 4 are to large extent dependent on the terrain characteristics and the obstacles present around the scanning area.The sapper usually cannot place it freely when he deploys the system's elements in a previously unsearched and potentially hazardous terrain.When approaching a minefield, he usually knows which part of the terrain is free of explosives and where the search should start.Thus, the most typical and safe locations for placing the UWB beacons lay in front and on the sides of the minefield, as shown in Figure 6.Such a system's geometry is certainly not optimal from the accuracy point of view, but even under the mentioned limitations, the actual placement of M 1 ÷ M 4 may significantly influence the positioning accuracy in various areas of the minefield.
To verify the mentioned dependence of the positioning accuracy on the locations of the UWB stationary modules, three system configurations with different locations of the M 1 ÷ M 4 modules were considered.The assumed positions of the modules are given in Table 1 and are graphically presented in Figure 9.
Table 1.Locations of the UWB stationary modules for different system configurations.

Configuration Number
Coordinates of the UWB Modules  We assumed that an area of 500 m × 500 m lying in front of the UWB beacons is divided by a grid with cells of 10 m × 10 m each.For every node of this grid, a set of ten thousand UWB measurements was generated in MATLAB ® , and its position was estimated using a simple iterative NLS algorithm, without any Kalman filtration.Based on the parameters of P440 modules, declared by their producer [55], we assumed that UWB measurement errors have zero-mean Gaussian distribution with a standard deviation of 2 cm.Next, RMS errors (RMSE) for each node were calculated and the obtained results for the three system's configurations are shown as colormaps in Figure 10.As the RMSE is very large in the vicinity of the UWB modules, the colormap is presented for the area where the y coordinate is larger or equal to 50 m.In practice, it means that the actual placement of the UWB beacons should be in the foreground of the minefield, far enough ahead of its border, to ensure that the positioning accuracy in the planned search zone is high.As can be seen, the smallest positioning errors are achievable in front of the place, where the UWB stationary modules M 1 ÷ M 4 are located.The high-accuracy zone is wider and deeper for a more extended baseline of the positioning system.
Mean and maximal RMSE values for the whole area of 450 m × 500 m and for smaller areas, limited to the nearest 100 m × 100 m and 50 m x 50 m respectively, in front of the UWB stationary modules M 1 ÷ M 4 are given in Table 2.As can be seen, the proposed positioning system can provide a centimeter level of accuracy in areas large enough for practical demining tasks and for reasonable and practically realizable systems' configurations.The high-accuracy zone could certainly be extended if the UWB stationary modules were more distributed around the area to be scanned, however for safety reasons it cannot always be achieved.

Positioning Accuracy
The accuracy of the EKF-PND filter was assessed and compared with the accuracy of EKF-CV and EKF-CA filters as well as with the accuracy of an NLS algorithm in MATLAB ® for the C1 configuration of the system.The dynamics and observation models given by Equations ( 4) and (8) were used to generate the antenna trajectories and the parameters chosen during the simulations are given in Table 3.The choice of these parameters was done in such a way that the shape and duration of the resulting antenna trajectory resemble typical HH-GPR antenna trajectories.The same standard deviations of all the distance measuring errors σ Aj = σ Sj = σ and power spectral densities S x S , S y S , S a given in Table 3 were used in the EKF-PND for setting the values of the R and Q matrices.The EKF-CV and EKF-CA also use σ Aj = σ Sj = σ as given in Table 3, but as their dynamics models are different, their Q matrices required finding power spectral densities of different noises S v x , S v y or S a x , S a y .This was done during the mentioned tuning process and the obtained values are as follows: and S a x = S a y = 6.1•10 −3 m 2 s 5 .Firstly, the results of the antenna position estimation with the EKF-PND and the NLS algorithm were compared for various orientations of scanning sections and chosen results of these tests are presented in Figure 11.These experiments confirmed that the EKF-PND filter works properly and achieves a similar accuracy for various orientations of the central axis of the scanning section.
Next, a closer inspection of the estimation results was done for all the implemented algorithms for a chosen orientation of the central axis of the scanning section equal to 45 • .
A comparison of HH-GPR antenna positions estimated with NLS, EKF-CV, EKF-CA, and EKF-PND is presented in Figure 12.As can be seen, all these algorithms are capable of properly estimating the antenna position, however, their accuracy is noticeably different and required further analysis, which will be presented further.
At this step of the simulations, the results of the estimation of other elements of the state vector x from the dynamics model given by Equation ( 4) were analyzed and they are presented in Figures 13-15.The angle θ between the horizontal projection of the antenna handle and the central axis of the scanning section, estimated with the EKF-PND filter, is shown in Figure 13.The results of angular velocity estimation are presented in Figure 14. Figure 15 contains an estimate of the acceleration a forcing the HH-GPR antenna movement.All these figures contain only results for the EKF-PND, as other algorithms do not estimate variables such as θ, ω, and a.Based on the above results we created a histogram of RMS antenna position errors for NLS, EKF-CV, EKF-CA, and EKF-PND and it is presented in Figure 17.From this and the previous figure, one can conclude that the EKF-PND is more accurate than all other tested algorithms and the EKF-CV and EKF-CA perform similarly, but still better than the NLS algorithm.The EKF-CA is slightly more accurate than the EKF-CV.A comparison of numerical values of average RMS antenna position errors for all the realizations and for the NLS, EKF-CV, EKF-CA, and EKF-PND algorithms are given in Table 4.This table also presents percentage improvements of accuracy for EKF-CV and EKF-CA vs. NLS and EKF-PND versus all other algorithms.As can be seen, in the chosen simulation scenario, the EKF-PND provides positioning results about 40% more accurate than other tested EKFs and about 60% better than NLS.

Discussion
In this paper, an accurate positioning system dedicated as a supporting tool for a handheld ground penetrating radar was presented.The system uses ultrawideband radio technology for accurate distance measurements and processes them to estimate the GPR antenna position.Various estimation algorithms were used for this purpose, from NLS, through simple EKFs (EKF-CV and EKF-CA), based on those typically used in radiolocation and navigation CV and CA dynamics models, to the EKF-PND, based on the proposed by the authors' dynamics model derived from the model of a pendulum motion.
The results of simulations included in the paper have demonstrated that the proposed positioning system can provide a desired centimeter level of accuracy in areas large enough for practical demining tasks.They also have shown how the actual placement of UWB beacons influences the system's accuracy.It occurs that the smallest positioning errors are achievable in some distance in front of the area where the beacons are located and that the high-accuracy positioning zone is wider and deeper for a more extended baseline of the system.
Further experiments have confirmed that the EKF-PND filter works properly for various orientations of the central axis of the scanning section and have proved that using the proposed pendulum dynamics model ensures a better accuracy than the accuracy obtainable with other typically used dynamics models CV and CV.The simulations have shown that the EKF-PND provides positioning results about 40% more accurate than other tested EKFs (EKF-CV and EKF-CA) and about 60% better than NLS.The final shape of the fundamental matrix F k−1 can be obtained by placing all its individual elements given by the Equations (A2)-(A10) at appropriate positions in (A1) and it is given below as Equation (A11).The final shape of the observation matrix H k , obtained by placing all its individual elements given by the Equations (A13)-(A16) at appropriate positions in (A12) is given below as Equation (A17).

Figure 1 .
Figure 1.General idea of GPR operation.

Figure 3 .
Figure 3. Examples of a horizontal section through a C-scan: (a) color image; (b) grayscale image.

Figure 7 .
Figure 7. Part of HH-GPR antenna trajectory (a single arc of the scanning profile).

Figure 8 .
Figure 8. Flowchart of the Extended Kalman Filter used for HH-GPR antenna position estimation.Using the calculated F matrix and the G matrix from the equation (4), Φ k,k−1 and Q k−1 matrices are obtained as follows [60-62]:

Figure 9 .
Figure 9. Locations of the UWB stationary modules for different system configurations.

Figure 13 .
Figure 13.Angle between the horizontal projection of the antenna handle and the central axis of the scanning section estimated with EKF-PND.

Figure 14 .
Figure 14.Angular velocity of the antenna motion estimated with EKF-PND.

Figure 15 .
Figure 15.Acceleration estimated with EKF-PND.To better compare the accuracy of estimation with various algorithms, we conducted a series of ten thousand simulations and calculated average RMS antenna position errors for the whole scanning sections for each realization of the simulations.The obtained RMSE values are shown in Figure16.Single points in various colors are RMS antenna position errors obtained with NLS, EKF-CV, EKF-CA, and EKF-PND.Although they are changeable in various simulation runs, they form bands on noticeably different levels.

Table 2 .
Mean and maximal RMSE values for areas of various sizes.

Table 3 .
Parameters used in dynamics and observation models during simulations.