Kinematic GPR-TPS Model for Infrastructure Asset Identification with High 3D Georeference Accuracy Developed in a Real Urban Test Field

Abstract

This paper describes in detail the development of a ground-penetrating radar (GPR) model for the acquisition, processing and visualisation of underground utility infrastructure (UUI) in a controlled environment. The initiative was to simulate a subsurface urban environment through the construction of regional road, local road and pedestrian pavement in real urban field/testing pools (RUTPs). The RUTPs represented a controlled environment in which the most commonly used utilities were installed. The accuracy of the proposed kinematic GPR-TPS (terrestrial positioning system) model was analysed using all the available data about the materials, whilst taking into account the thickness of the pavement as well as the materials, dimensions and 3D position of the UUI as given reference values. To determine the reference 3D position of the UUI, a terrestrial geodetic surveying method based on the established positional and height geodetic network was used. In the first phase of the model, the geodetic network was used as a starting point for determining the 3D position of the GPR antenna with the efficient kinematic GPR surveying setup using a GPR and self-tracking (robotic) TPS. In the second phase, GPR-TPS system latency was quantified by matching radargram pairs with a set of fidelity measures based on a correlation coefficient and mean squared error. This was followed by the most important phase, where, by combining sets of “standard” processing routines of GPR signals with the support of advanced algorithms for signal processing, UUI were interpreted and visualised semi-automatically. As demonstrated by the results, the proposed GPR model with a kinematic GPR-TPS surveying setup for data acquisition is capable of achieving an accuracy of less than ten centimetres.

1. Introduction and Motivation

The key motivation for the study was to develop an improved kinematic ground-penetrating radar (GPR) model for the acquisition, processing, visualisation, positioning and mapping of underground utility infrastructure (UUI), compliant with present day demands and accuracy standards. Positioning and mapping of UUI in urban areas is perhaps the most complicated of GPR exercises among all types of civil engineering applications. This is because radargram patterns of the urban areas of utility orientations, diameters, depths, lateral material types and strata are often non-typically compared to other infrastructures [1]. Underground utility infrastructure includes a privately-, publicly-, or co-operatively-owned line, facility or system for producing, transmitting or distributing communications, cable television, power, electricity, light, heat, gas, oil, crude products, water, steam, waste or any other similar commodity, including any fire or police signal system or street lighting system [2].

Accurate knowledge of geometry and the horizontal and vertical positions (3D position) of UUI is a fundamental aspect of their efficient management and optimal planning, while it also reduces maintenance costs and has a minimum effect on traffic and users [3,4]. There is a large number of UUIs that have inadequate information about their 3D position and basic properties, which is unacceptable in modern times. In many databases, where information does exist, the only reference is often digitised and vectored, usually consisting of unmaintained cadastre plans of utilities that have been transferred to geographical information systems (GIS). Users usually accept the information of UUI as reliable and accurate datasets. They normally ask themselves about the accuracy and reliability of data only in circumstances whereby they realise that the information is significantly less accurate than expected, and, as a result, find themselves in an annoying situation.

Recording UUI and, consequently, acquiring much better knowledge on the geographic location in three-dimensional space, could prevent many accidents and injuries involving UUI worldwide, which are manifested in the interruptions of supply and gas explosions [5]. In the UK alone, there are millions of miles of underground utilities with often inaccurate, incomplete, or non-existent location records [6]. It is estimated that the total length of UUI in the US is in excess of 56 million kilometres [7]. More than half of the existent UUI in the US does not have a known horizontal position in the area [8]. Underground utility infrastructure is inadvertently struck every minute in the US, equating to nearly 500,000 annual utility strikes [9,10]. The societal cost associated with UUI damage in the US in 2016 is estimated at $1.5 billion [11]. That is why a systematic improvement of GIS data of UUI is necessary, based on the criteria of integrity and reliability. The investigated study showed a positive return-on-investment ranging from US $2.05 to $6.59 for every dollar spent on improving underground utility location data [12].

Thomas et al. [6] listed the degree of horizontal accuracy of the GPR method in the UK required by utility stakeholders. Stakeholders require horizontal position errors to be less than 10 cm and definitely no higher than 30 cm. The allowable position errors are guided by ASCE 38-02 from the USA [2], AS 5488-2013 from Australia [13] and ICE PAS 128:2014 from the UK standards [14]. These standards categorise the GPR method as the second-best quality level. The ICE PAS 128:2014 sub-divides the accuracy into horizontal and vertical accuracies as a function of the detected depth. The best horizontal accuracy is ±15 cm or ±15% of the detected depth, whichever is greater, whilst vertical accuracy is ±15% of the detected depth.

On the basis of the horizontal accuracy of coordinates, which is led by the consolidated cadastre of public infrastructure in the Republic of Slovenia, it is estimated that more than 73% of the UUI is determined with a low accuracy of one metre or more [15]. In Slovenia, the one-year overall direct damage has been estimated at €2.1 million or €1.04 per capita [16]. The interest among local owners and operators of UUI is shown to improve the 3D positioning accuracy and the detection of unregistered objects using a non-destructive method. The horizontal and vertical accuracy of recently registered UUI in Slovenia are introduced by the policies, recommendations and requirements of owners in tenders and contracts. The horizontal accuracy of the coordinates of points, assessed with a semi-major axis of a standard ellipse, must not be less than 12 cm. Furthermore, the vertical accuracy of the points estimated by the standard deviation of the height must not be less than 10 cm.

In general, the methods of detecting and recording the existing UUI are either partially destructive or totally non-destructive. Non-destructive GPR is one of the most effective tools for subsurface investigation, detection and recording of UUI [17]. The identification of underground structures from the raw GPR data (radargrams) is a complex, non-trivial and sometimes also intuitive task [18]. Ground-penetrating radar has the ability to detect and locate both metallic and non-metallic objects, and also depth estimation, and has a better resolution than other technologies [19,20,21,22,23]. Underground utility infrastructure in urban areas is normally within a few metres (usually up to 2.5 metres) from the surface, which falls well within the GPR survey range (200–1000 MHz) [24,25,26]. The maximum depth of investigation decreases when the electromagnetic wave (EMW) frequency increases. The central frequency and the bandwidth of the electromagnetic signal determine the system resolution, namely, the minimum physical size of the UUI can be detected.

Over the last decade, kinematic GPR surveying in relation to real-time 3D georeference methods has been studied and enhanced. In this sense, there are fieldworks, such as [27,28,29,30], which have presented kinematic acquisition strategies simultaneously acquiring GPR and positional data. The greatest difficulty in this strategy is presented by synchronisation and/or merging of GPR and positional data subsequent to the field survey.

This paper focuses specifically on methods that enable real-time data fusion approaches. The authors’ goal is to achieve a wireless connection between GPR and a positional instrument, for which it is intended that a terrestrial positioning system (TPS) with a high-precision self-tracking electronic tachymeter will be used. The TPS automatically tracks a 360° prism mounted on the GPR antenna’s sledge. The fundamentals were set by Lehmann and Green [26]. Kinematic GPR surveying using a precision self TPS method (GPR-TPS method) was, in terms of a wireless connection between GPR and TPS, upgraded by Böniger and Tronicke [27], who used a radio link. Such a GPR-TPS method, however, has radio-link crosstalk limitations and systematic latency.

The advantages and the application of the presented kinematic GPR-TPS methods are justified by the fact that the single-frequency global navigation satellite system (GNSS) receivers do not meet the required positioning accuracy. A significant displacement of the single-frequency kinematical GNSS fixes (up to 8 m) from the TPS trace along the same route can be observed [31]. When receivers are dual-frequency GNSS, the surveying method (real-time kinematic (RTK) or the virtual reference station network concept (VRS-RTK)) is limited due to the interference of GPR and GNSS waves, the limited signal reception and the impact of multipath on the observation in urban environmental variables, where the most UUIs are installed. Hence, accurate GNSS positioning is not always available in urban GPR surveys and an alternative method is required [4]. An option for this alternative is the rotary laser positioning system method [30], which represents a real-time approach that uses a cable for the fusion of GPR and positional data in real time.

For the purposes of this study, the so-far partially developed processes were joined into a comprehensive GPR-TPS model for the acquisition, process and visualisation of UUI. One of the important goals of this study was to demonstrate the applicability and performance using a contemporary TPS system for kinematic GPR data acquisition with the required accuracy in real urban test pools (RUTPs). The aim was to develop a model that does not require: a cable connection between GPR and TPS and/or additional hardware and set-up of a reference line (e.g., orthogonal grid) in order to navigate GPR along it. In order to avoid radio-link crosstalk, a Bluetooth wireless technology standard was applied (BT-232B-E Bluetooth/RS232 adapter and Carrier Reflector Mount). To improve kinematic GPR-TPS model 3D positioning accuracy, a systematic latency will be evaluated by using image matching of radargram pairs. Real urban test pools have been designed for this purpose and are among the rare ones in the world [32,33,34].

The rest of the paper is organised as follows: Section 2 presents the construction of an RUTP, buried UUI and the establishment of a geodetic reference network. Implementation of the proposed kinematic GPR-TPS model and estimation of systematic latency are presented in Section 3, which is followed in Section 4 by the capturing, processing and visualising of GPR data. Section 5 presents an interpretation of the results obtained with the proposed GPR-TPS model. Section 6 provides some final conclusions and directions for future work.

2. Design and Layout of Real Urban Testing Pools: The Establishment of a Geodetic Reference Network and an Estimation of the Positional and Height Accuracy of the Measured UUI

A large number of uncontrolled variables in the real environment prevent the creation of a unique sequence of processing GPR signals for efficient and reliable identification of all the different individual elements of UUI. For the purposes of this study, a test environment, consisting of all the key elements of real urban environments, was established. This was done in order to develop and assess the GPR-TPS model in representative and controlled conditions (accurate knowledge of the thickness, structure and material of the layers of subsurface and position, depth, inclination, diameter and material specifications of the types of UUI used) by design and construction of the RUTP. The constructed RUTP, with all the important contemporary UUI in Slovenia, provides a sufficient approximation of the design of pavement structures for regional and local roads as well as other pedestrian surfaces. The latter was carried out on the basis of today’s applicable technical specifications for the design of pavements in Slovenia [35,36,37,38].

The existing concrete storage blocks with a roof in an abandoned storage area were reused, (Figure 1). In the first block, with dimensions in the plan view area of 6.03 m², a pedestrian pavement was built (Figure 2a). In the second block, with area dimensions of 5.73 m², a roadway of a local road composition was built (Figure 2b), while in the third block, area dimensions of 6.21 m², a roadway of a regional road composition was built (Figure 2c). Depending on the specified (expected) traffic load bearing strength of the base and the climatic characteristics of the environment in which they are located, the thickness of the layer of the individual structural roadway and pedestrian walkway was calculated and the structure and material of each layer, as shown in Figure 2, was set. The formation level represents the soil with some rock debris, while the subbase and base course were represented by crushed rock aggregate and consisted of 100% limestone grains of a diameter of 0–125 mm and 0–32 mm. The asphalt course was represented by a layer of asphalt concrete (AC 22 base 50/70 and/or AC 8 surf 50/70).

For the purpose of this research, different types of infrastructure in the RUTP on the basis of the existing legislation, policies and technical specifications for the installation of UUI (see

Table 1. A summary of all the buried pipes, electrical cables, electronic communication cable (ELC), coax cables and cable duct in the real urban test pools (PVC—polyvinyl chloride, PE—polyethylene, Al—aluminium, Cu—copper, ND—nominal diameter, OD—outside diameter).

Object No.	Type	Manufacturer	Material	Nominal Diameter ND/OD [mm]	Depth [cm]	Pool No.
1	Gas	Drinplast	PE 100	100	120	2
2	Gas	KontiHidroplast	PE 100	63	110	3
3	Gas	TotraPlastika	PE 100	32	80	2
4	Water	TotraPlastika	PE 80	90	130	2
5	Water	TotraPlastika	PE 80	110	150	3
6	Water	Deriplast	PE 80	32	80	2
7	Water	Tubi PVC	PVC	110	150	3
8	Water	Tubi PVC	PVC	110	160	1
9	Water	SvoodnySokol	Duktil	110	180	1
10	Sewage	OsterndorfKun.	PE 200	125	170	1
11	Electrical	Elka	Al/PE HD	31	80	3
12	Electrical	SKW	Al/PE HD	31.2	80	3
13	ELC	Minerva	PE 80	125	80	1
14	ELC	Rumaplast	PE 80	2 × 50	80	1
15	ELC	Miner.- Mapitel	PE 80	40	40	1
16	ELC	Miner.- Mapitel	PE 80	50	40	2
17	ELC	Minerva	PE 80	110	80	1
18	Coax	Commscope	Cu/PE 80	24.4	80	3

) were used. When installing UUI, a sand bed with a thickness of about 8–10 cm made of grains of sand of 0–4 mm was prepared and covered with the same material with a thickness of approximately 10–12 cm above the apex point.

The Establishment of a Geodetic Network for the Purpose of Determining the Position of the UUI and GPR

For the purpose of determining the relevant position of the UUI in the 3D space and the data captured by the GPR, a geodetic network in the vicinity of the RUTP was developed and established, which represents the geometrical basis. The ETRS89/TM national coordinate reference system of Slovenia was used [39].

In order to obtain the approximate positions of the geodetic network, the GNSS VRS-RTK method [40,41] was used. The main objective of the positioning and height adjustment of the observations was to determine the most probable positions and heights of the geodetic network points with high precision on the basis of conventional terrestrial geodetic observations (see Figure 3). A detailed description of the adjustment is found in Appendix A.

For the purpose of determining the positions and heights of UUI, a polar coordinate surveying method (classical geodetic method) for detail points was carried out. The method used directly determines the relative coordinates of the detailed points in relation to the given points of the geodetic network. Trigonometric levelling was applied to estimate the height differences between the reference point in the network and the UUI. The horizontal position and height of the object of UUI were recorded with the top point (apex) of the pipe or cable every 10–15 cm.

The 3D display of the measured UUI in the reference coordinate system with the installed cables and pipes and also layers of pavement structures as well as subbase and subgrade are shown in Figure 4. The occupancy of the 3D space is illustrated by the means of a known outer diameter of each object. The UUI connections of smaller outer diameters are often curved (Object No. 3, 6, 11) due to the obstacles in the real environment. The UUI are even more difficult to detect and record in cases where they are not in a straight line. Table 1 describes the characteristics of the installed and measured UUI in the RUTP in detail.

The accuracy of the position and height of the installed UUI was determined by the standard deviation of the measured apex points using the geodetic method. Taking into account the impact of the positioning and height precision of network reference point 1001 (Figure 3), the law of propagation of variance–covariance [42] was applied to calculate the accuracy of the position and height of the UUI in RUTP, considering the accuracy of determining the angles and lengths provided by the instrument used. Taking into account the error of centring the TPS over the point 1001 ${\sigma }_c$ and signalling the detailed point ${\sigma }_{cs}$, the following was obtained:

$${\sigma }_{E_i,N_i}^2={\sigma }_{N_i}^2+{\sigma }_{E_i}^2+{\sigma }_c^2+{\sigma }_{cs}^2$$

(1)

where ${\sigma }_{E_i,N_i}^2$ is the variance of the position of the detailed apex point $i$, ${\sigma }_{N_i}^2$ is the variance of the detailed apex point $i$ in the $N$ direction, ${\sigma }_{E_i}^2$ is the variance of the detailed apex point $i$ in the $E$ direction, ${\sigma }_c$ is the error of centering and ${\sigma }_{cs}$ is the error of signalling.

The height accuracy of the detailed points was calculated using the following equation [42]:

$${\sigma }_{H_i}^2={\sigma }_{H_{1001}}^2+{\cos }^2{z}_{1001}^i·{\sigma }_d^2+d_{1001}^i{}^{2}si{n}^2{z}_{1001}^i·{\sigma }_z^2+{\sigma }_i^2+{\sigma }_l^2$$

(2)

where ${\sigma }_{H_i}^2$ is the variance of the detailed apex point of the height, ${\sigma }_{H_{1001}}^2$ is the variance height of the reference point, $z_{1001}^i$ is the zenith distance to the detailed apex point $i$, $d_{1001}^i$ is the slope length to the detailed apex point $i$, ${\sigma }_d^2$ is the variance of length, ${\sigma }_z^2$ is the variance of zenith distance, $\mathrm{and}{\sigma }_i^2$ and ${\sigma }_l^2$ are the variances of establishing the height of the instrument and the signal.

The lowest accuracy calculated from ${\sigma }_{E_i,N_i}=2.1$ mm and ${\sigma }_{H_i}=4.3$ mm of individual detailed points has been used as an actual positional and height accuracy of the measured UUI in the RUTP.

3. Implementation of the Kinematic GPR-TPS Model

The GPR observations were performed in a time domain. An orthogonal grid of 20 × 20 cm on the surface of the RUTP only allowed the simple and precise movement of the GPR antenna. A total of 85 longitudinal and transverse GPR profiles were recorded. Longitudinal profiles took place parallel to the test pool walls. Shorter transversal profiles were parallel to each other and perpendicular on longitudinal profiles (Figure 5). All profiles were recorded in straight lines on the ground with a very small inclination (less than 0.8%). Monostatic shielded antennas with a central frequency of 270 MHz (120 ns, 12 profiles), 400 MHz (80 ns, 59 profiles) and 900 MHz (30 ns, 12 profiles) were used.

Changes in the GPR position were traced using a TPS measurement system, which uniquely determines the position. Since the position in the kinematic GPR-TPS model is a function of time, mathematical adjustments based on a set of redundant observations are not possible. The accuracy of the kinematic GPR-TPS model was determined on a priori information on the measurement system and the method chosen. The expected accuracy of the reflector position in the kinematic self-TPS method is largely dependent on the size of the motion vector (the speed and direction of movement), frequency of measurements, physical obstacles in the area, precision of the reflector, the transmission time or latency between the two systems (GPR and TPS) and the distance from the TPS. The manufacturers recommend a distance of at least 10 m from the TPS, due to the limitations in the angular speed of the TPS movement. Terrestrial positioning system is used today in many high-precision applications that require sub centimetre positioning precision [43].

Figure 6 shows the measurement of horizontal angles, zenith distances, and lengths in the kinematic method. The self-TPS Leica TCRP 1201+ with automatic target recognition, 360° Leica reflector GRZ4 and GPR instrument SIR-3000 (geophysical survey system) upgraded with a carrier reflector were combined. A positional update rate of up to 5–6 Hz was achieved at the speed of a moving reflector ≈ 0.5 m/s. The positional accuracy refers to ISO 17123-4, and the used TPS was (0.003 m; 1.5 ppm) in auto-tracking mode [44]. In order to avoid cable connections and for improved flexibility of fieldwork, a real-time wireless Bluetooth connection was established that allowed the fusion of data between the two systems. For this purpose, it was necessary to upgrade the GPR assembly with a BT-232B-E Bluetooth/RS232 (4800 baud rate, 8 bit, 1 stop bit) interface with an external dipole antenna through the use of a belt frequency between 2402 MHz and 2483.5 MHz. The fusion of data between the GPR system and TPS was enabled through a pseudo-GGA NMEA (National Marine Electronics Association) string [44].

Connecting kinematic TPS surveying and GPR observations were enabled by a set of our own software solutions (DST_converter, Parse and Run) developed for the purpose of this research and written in scripting languages (Python, PHP and Bash). Only the position, height and time of the initial and final track on the GPR profile were recorded in the time mark file (TMF) of the GPR control unit. All other measured positions, heights and times are recorded in the on-board TPS data storage. The DST_converter reads the binary files distance tracenumber (DST) ASCII format and TMF of the individual profile and records one trace in each row. It registers the initial and the final trace and their corresponding time, which is corrected for the value of the calculated latency (see sub-Section Latency). Parse software solution in the first step uses a linear interpolation to attribute time to individual traces without any data. Then the programme, with an accuracy of 10⁻³ s, finds and attributes the position and height of every individual trace of data on the objective time in the TPS data file. The traces are attributed a position and height in a new DST file. The format and the content of the DST files recording are enabled by the software for processing Reflexw 7.5 [45], which is attributed to the radargram. Figure 7 shows 61 recorded points with the GPR-TPS method, along with the profile at a measurement frequency of 5–6 MHz and the speed of movement of the reflector ≈ 0.5 m/s.

If the position and height are not known at a given time, then they are not attributed to the trace. Empty traces (the spaces between the blue points in Figure 7) were attributed to coordinate values between the two closest traces with a known position and height using linear interpolation.

Latency

The positional accuracy of the kinematic GPR-TPTS model mainly depends on two factors: the accuracy of the TPS measurement system used and systematic latency. In this case, latency refers to the time delay between the actual measured position and its record on the GPR trace [27].

In the kinematic GPR-TPS model shown in Figure 6, the time delay can be attributed to the processing of data in TPS and GPR as well as to the transfer of data via a wireless interface. On the basis of these time delay causes, a total latency can be derived, which is the sum of all the individual latencies of the hardware presented with the equation:

$$L_{\mathrm{total}}=L_{\mathrm{TPS}}+L_{\mathrm{Bluetooth}}+L_{\mathrm{GPR}}$$

(3)

where $L_{\mathrm{total}}$ is the total latency of the GPR-TPS model, $L_{\mathrm{TPS}}$ is the latency system of the electronic tachymeter, $L_{\mathrm{Bluetooth}}$ is the latency of the wireless interface and $L_{\mathrm{GRP}}$ is the latency of the GPR system.

Linking the position measured with the TPS and GPR data using time is incorrect for the time delay that is represented by total latency. Individual traces in the GPR data attributed the position at the time of the measurement $t_m$, which is incorrect because, due to the time delay, this is the position of one of the previously measured traces [27]. The position of the trace should be written at the correct time $t_p$. The correct time $t_p$ and consequently the correct position of the trace can be calculated when total latency, represented by $L_{\mathrm{total}}=t_m-t_p$, is known.

In order to estimate total latency, the same profile on the S–N (forward $r_f$) and on the N–S (reverse $r_r$) directions, assuming the same rate of GPR walking speeds, were recorded. The result of latency is laterally offset between radargrams $r_f,r_r$ for the latency value as shown in Figure 8.

Assuming the identical depth of the target object, the image matching when detecting lateral displacement between the pair of radargrams $r_f,r_r$ was used. Objective methods of image matching were divided into three main categories: feature-based matching, area-based matching, and structural matching. The relative simplicity of the method places area-based matching among the most commonly used [46]. Using an area-based matching method and selecting the targeted object in the form of a hyperbole on the forward radargram $r_f$ and automatically locating the most probable homologous reflection on the reverse radargram $r_r$, an alignment of the radargrams or calculation of their lateral offset was performed. When matching images, there may be some pixel differences in vertical directions. In such cases, it was necessary to move one of the images up or down in a vertical direction before a lateral offset was performed. Before using the method of area-based matching, the radargrams were treated with a series of procedures of processing phase III.

One of the methods to describe the alignment or match of a corrected radargram of a profile pair is called the correlation technique [47], which is based on matching the selected area of two images. The areas are typically square-shaped. This is basically a movement of one image towards the other and a calculation of their mutual matching in individual movements. The place with the highest correlation value corresponds (1 represents a perfect match and 0 indicates a complete mismatch) to the solution in which the images are most aligned and defines the parameters searched (offset in the number of pixels). The accuracy of the procedure is the size of one pixel; however, it can be improved by means of interpolation methods. Similarity can also be looked for among the images in which there are minor discrepancies because of the rotation or size, which is common in radargrams in an urban environment [47,48].

Different fidelity measures were used as a unit of the image alignment and matching of the radargram pairs $r_f,r_r$. The first measure was the correlation coefficient Corr while the second was the mean squared error (MSE). Matching is based on a similarity examination of individual pixels on the radargram pairs $r_f,r_r$. A calculation of the correlation coefficient is derived from the standard deviation ${\sigma }_{r_f}$ and ${\sigma }_{r_r}$ of radiometric values in the section of the radargram pairs $r_f,r_r$ and their covariance ${\sigma }_{r_f{r}_r}$:

$$r_{r_f{r}_r}=\frac{{\sigma }_{r_f{r}_r}}{{\sigma }_{r_f}{\sigma }_{r_r}}$$

(4)

$${\mu }_{r_f}=\frac{1}{NM}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M{r}_f\left(i,j\right)$$

(5)

$${\mu }_{r_r}=\frac{1}{NM}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M{r}_r\left(i,j\right)$$

(6)

$${\sigma }_{r_f}=\sqrt{\frac{1}{\left(N-1\right)\left(M-1\right)}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M\left(r_f\left(i,j\right)-{\mu }_{r_f}\right)}$$

(7)

$${\sigma }_{r_r}=\sqrt{\frac{1}{\left(N-1\right)\left(M-1\right)}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M\left(r_r\left(i,j\right)-{\mu }_{r_r}\right)}$$

(8)

$${\sigma }_{r_f{r}_r}=\frac{1}{\left(N-1\right)\left(M-1\right)}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M\left(r_f\left(i,j\right)-{\mu }_{r_f}\right)\left(r_r\left(i,j\right)-{\mu }_{r_r}\right)$$

(9)

where N and M represent the radargram dimensions, $\left(i,j\right)$ are the values of the (i,j)th samples in $r_f$ and $r_r$, $r_{r_f{r}_r}$ is the correlation coefficient, ${\sigma }_{r_f{r}_r}$ is the covariance of radiometric values of a pair of the compared radargrams, ${\sigma }_{r_f}$ and ${\sigma }_{r_r}$ are standard deviations, and ${\mu }_{r_f}$ and ${\mu }_{r_r}$ are the median values of individual images [49].

The MSE is based on the calculation of the average square difference of pixels between the compared radargram pairs $r_f,r_r$. The method is attractive because of the simplicity of the calculation and the clear physical meaning [49,50,51]:

$$MS{E}_{r_f{r}_r}=\frac{1}{NM}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^M{\left(r_f\left(i,j\right)-r_r\left(i,j\right)\right)}^2$$

(10)

The procedures may only be used if the radargrams are of the same resolution. In assessing the latency, 12 longitudinal profiles measured with a 400 MHz antenna in both directions were used. Assuming the same speed in both directions, the position of the reflection of the installed pool target object should be the same or within the limits of precision of the used kinematic GPR-TPS model.

Figure 9 shows the estimated latency of the 12 profiles at walking speeds of between 0.49 and 0.34 m/s. Due to the nature of walking speed, uneven movement can occur. The selection of reflection of the target object for the calculation of the latency estimation was conditioned by the smallest deviation of speed movement back and forth over a distance of 0.5 m within each pair of profiles. When calculating the latency, both speeds were averaged. Both matching methods resulted in statistically almost identical results of the estimated latency $L_{\mathrm{total}}$, namely, 0.515 s with the correlation coefficient and 0.511 s with the MSE with a standard deviation of ±0.031 s to ±0.029 s.

When calculating the corrections of the kinematic GPR-TPS model, the mean latency was calculated with the MSE fidelity measure due to the smaller dispersion. An example can be found in profile 10, where the impact of the estimated latency was 0.511 s, at an average walking speed of 0.475 m/s and a size of 0.24 m with a standard deviation of 0.014 m.

Our software solution was developed in order to calculate the matching and alignment of radargrams and to evaluate the results and effectiveness of image matching. The “Settlement” programme was written in the C++ programming language. Sections of radargrams were subtracted in an attempt to evaluate the results and the impact of image matching. Individual sections or parts of the radargram pairs were subtracted by subtracting the value of the grey levels of the individual pixels of the radargram pairs (forward $f_{rf}\left(i,j\right)$ and reverse $g_{rr}\left(i,j\right)$), where the difference was formed as $s\left(i,j\right)=f_{rf}\left(i,j\right)-g_{rr}\left(i,j\right)$. The impact of latency was presented, where the images of the subtracted sections were with the latency under-corrected, with the latency corrected and with the latency over-corrected by two (Figure 10, red box). The best match, where the reflection of the target object was subtracted, is given by the middle Figure 10b obtained by subtracting the sections corrected for the estimated latency by the MSE fidelity measure. Full subtraction cannot be expected for the measurements under real circumstances.

4. The GPR-TPS Model: Capturing, Processing and Visualising Data

The proposed model is divided into several stages, as shown in the schematic flow diagram in Figure 11. The model was based on standard and alternative methods and algorithms for processing the GPR data. The proposed model upgrades the classic GPR methodologies in terms of introducing an upgrade of the software and hardware of the existing GPR system and the one adapted for obtaining metric data in real time. In phase I, the model software was upgraded with an automatic integration of geophysical and geodetic datasets by time. The software was also upgraded in terms of the latency estimation of the proposed model used in phase II and upgrades in phase I. In this model, the software solutions developed for this research in phase I and phase II were introduced (see Section 3 and sub-Section Latency). Phase III consisted of a sequence of selected processing procedures of GPR data, which were already applied in phase II, as shown in Figure 11 (blue box). Profile processing (Phase III) was based on the authors’ recommendations [1,52,53]. An identification basis of the analysis in the model was introduced by selecting the target objects. A 3D interpretation and visualisation of data results were carried out, where the spatial relationships between the selected targeted objects were shown. This was represented in phase IV.

The final quality of the model, in addition to the chosen processes and specific targeted objects, is influenced by the assessment of adequate data processing flow, types of processes and appropriate set of parameters for each process [34]. Yilmaz [52] determined the success of the processing with the appropriate selection and sequence of procedures, selecting the correct set of parameters and also evaluating the results and the disadvantages of each procedure. Each site was determined by its features and effects that have been recognised by an experienced operator and considered in the model in such a way that it expands and improves it. The GPR data were processed using the following procedures.

For further 3D processing, the software used required the same number of traces within the set of the profiles of the same length. Unification of the traces was also recommended for the unified determination of the velocity field of the whole set of profiles. The conversion of two-way travel time to depth and Kirchhoff’s migration with 2D velocity field was dependent on the number of traces of the processed profile. Due to the changes in the thickness of the pavement layers along the profile, a 2D velocity field was used for each profile separately. The estimation of the EMW propagation velocity in the medium was based on generalisations and assessment of the medium properties, which was, consequently, one of the most complex tasks (see sub-Section Determining the (Velocity) Depth). The conversion of time to depth with the 2D velocity field was anticipated in the model as a basis for determining the depths of the target objects. The range of the depth axis was conditioned by the depth of the UUI and depth capability of the used antenna, which was dictated by the frequency range of the antenna. The model included Kirchhoff’s time migration after which a manual signal gain correction was applied. The model anticipates a preliminary procedure of a trace time cut, which is conditioned by the assessment of the EMW velocity in the medium. Having no knowledge of the estimated velocity, it was unreasonable to apply a trace time cut without the risk of losing important information. Removal of the initial DC bias was applied by a DC shift filter. The DC-shift filter converts GPR signals to a form with zero mean, subtracting their constant component. The model takes into account a time zero correction, where the manual method of determining and the delay of the start time are used. It is necessary to place the zero time at a clearly definable and stable location in the early wavelet. This is commonly at either the negative or positive maximum peaks of the first wavelet, or the zero amplitude point between these two peaks.

The tests show that the use of manual determination is reasonable, due to the presence of direct and air echoes, which could mislead the automatic process of the time zero correction and indirectly considerably influence determination of the depth of the target objects. The size of the time zero intervals is in correlation with the central frequency of the antenna used. Removing the background is a widely used method for removing noise and is, as such, indispensable in any model of data processing [46]. The model anticipates the removal of the background along the entire profiles. Following a thorough examination of the three gain functions, a manual gain (y) function [45] was applied, the gain values (db) of which were manually entered. Manual gain allows for any gains along the profile with a focus on the target object: however, for proper implementation it requires a greater number of iterations and the knowledge of the signal strength decay curve, which is calculated for each central frequency of the antenna. If necessary, the automatic gain functions (AGC) or energy decay can also be used. A band-pass frequency with a tapered cosine window was applied to sharpen and smooth the noise in the model after using two band-pass filters.

The model also includes the process of subtracting the average, which needs to be handled with extreme caution. It is useful in smoothing the impact of horizontal reflectors. Introducing the length of the window depends on the desired degree of smoothing and the central frequency of the antenna. For the removal of linear reflections with an inclination, which were caused by wooden barriers between the RUTP, an f-k filter was used in the model. The definition of the velocity interval parameters for the positive and negative parts was important when trying to remove unwanted reflections of this type. The parameters varied among profiles, so it was impossible to set unique parameters, which could be useful in a series of profiles. Since all the profiles were set with a height component, topographic correction was also included in the model. The importance and the results did not have a significant impact on the GPR data.

Assuming the data has small inclines (up to 3%), the topographic correction gives satisfactory results. Changing the design and integrating topographic migration seems reasonable with inclines steeper than 3%. An analysis for the recognition of objects was applied by a manual picking of UUI and sub-horizontal horizons of pavements. The envelope of distinct GPR echoes was calculated using a quadrature Hilbert transformation, which often assures better results of visualisation at spots of relatively weaker signals. It makes sense to display the results with time slices or with 3D modules in three perpendicular plains $x$, $y$, or $z$, as individual rasters and an interactive display in video format. Three-dimensional data visualisation makes it easier to show the spatial relationships between the selected UUI (see Supplementary Materials, Video S1).

The whole procedure of the proposed model and the schematic representation of the sequence of capturing processes, processing and visualisation, is shown in Figure 12.

Dependent and independent processes related to the central frequency of the antenna used were included in the model. The datasets and parameters independent of the central frequency of the antenna were: the integration of geodetic and geophysical datasets, linear interpolation of the position and the height of the blank traces, latency evaluation, reintegration geodetic and geophysical datasets taking into account the latency, trace extract, time cut, topographic correction, background removal, picking objects and 3D modelling.

Table 2. The procedures and the parameters used in the proposed model of GPR dependent on the central frequency of the antenna.

Process	Parameters 900 (MHz)	Parameters 400 (MHz)
DC shift, interval (ns)	20–25	50–70
Time zero correction (ns)	3.90	5.73
Manual signal gain, gain factor (dB)	0–36	0–33
Band-pass frequency with tapered cosine window (MHz)	420/630/1370/1700	230/320/580/750
f-k filtering limited by reflections for the positive and negative directions (m/ns)	+0.105 to +0.065 −0.055 to −0.085	+0.098 to +0.057 −0.043 to −0.072
Subtracting average (traces)	65	50
Determination of 2D velocity field, interval (m/ns)	0.080–0.132	0.080–0.131
Kirchhoff 2D time migration, ∑ width (Number of traces)	20	30
Manual signal gain, gain factor (dB)	0–25	0–25
Time to depth conversion, max depth axis (m)	1.6	2.1

shows the parameters used in the proposed GPR-TPS model, which were dependent on the central frequency of the antenna.

Determining the (Velocity) Depth

The basis for determining the depths of the UUI is the depth migration of radargrams from the time $x$, $t$ to the depth $x$, $z$ domain, for which it is advisable to know the 2D velocity field or the velocity of EMW in each separate layer of the pavement. There are several methods of EMW velocity estimation in the medium or its parts [54]. The velocity model consists of: a suitable method for EMW velocity estimation in the whole subsurface or in individual layers, a suitable time depth conversion, and a geometric correction in the form of the migration process. The use of constant mean velocity in the pavement structure, where dealing with multiple layers and UUI at different depths, is rather rough and often proves to be quite an inaccurate estimation. The problem occurs with layered media, where the velocity of EMW is required for every single layer with a different composition. It is reasonable to use the 2D velocity model, as depth migration, time migration and depth conversion can be used. Velocity function is normally provided by a 2D velocity file.

The approximation of the velocity function can roughly be estimated by hyperbolas fitting, which occurs in cylindrical or point objects’ reflections, by using the hyperbolic velocity analysis and changing the speed which determines its shape or slope (Figure 13). Geometry, or the inclination of hyperboles, is a function of the velocity of EMW in the medium and the diameter of the target object [55]. A condition for the use is clean and well-expressed hyperboles in the GPR radargram. The method is limited only to those radargrams where the hyperboles are very visible and velocity changes across the medium can, therefore, be estimated on the basis of different hyperbola spatial distribution. This may be a problem in radargrams with few or no hyperboles. The estimated velocity of the 2D velocity file with hyperbola fitting are burdened with an error of unknown diameter between 0.135 and 0.09 m/ns.

In this study, both lateral and smaller vertical structural changes are visible on the GPR radargrams. The differences in the vertical direction are sub-horizontal GPR echoes between separate layers of pavements. The layer velocity method was used for estimation of the pavement depth. The Kirchoff’s time migration with an estimated average rate of EMW velocity (0.13 m/ns) was applied, the result of which is the source horizontal reflection between layers. Knowing the approximate thickness of the two-way travel time, the mean EMW velocity can be estimated for individual layers of pavements. In assessing the velocity of the lowest layer, where lateral reflections do not appear, the known depth of the referential target objects No. 5 and 7 were used (see Table 1). Figure 14 shows the sub-horizontal GPR echoes indicated by lines which represent the boundaries between the pavement layers (a) and the velocity field of the estimated 2D velocity files of the pool profiles, which range between 0.132 and 0.080 m/ns (b).

Figure 15 shows a time-to-depth conversion with the velocity field obtained with the layer velocity method. In the last part of the profile (pool 3), the positions and depth of three distinct reflections were labelled. These are the echoes from buried power lines (objects No. 11 and 12, red circle) and water pipes (objects No. 5 and 7, blue circle).

Figure 16 shows a time to depth conversion of the velocity field obtained with the hyperbola fitting method. Certain discontinuity is clearly visible (red square) on part of the GPR profile. This indicates the incorrect selection of the EMW velocity (Figure 13). The positions and the depths of three distinct reflections (the same as in Figure 15) of the buried infrastructure are marked.

Both methods of time-to-depth conversions have their advantages and disadvantages. The first one is conditioned with the number, spatial distribution, and the diameter of GPR reflections in the form of a hyperbole, while the other is conditioned by knowing the actual thickness of the pavements and the accuracy of determining linear reflections. The thickness of the pavements is often different from those in the project documentation. Both methods are significantly dependent on subjective assessments by the operator. However, the method of layer velocity has an advantage in terms of an easier identification of lateral reflections along the pavement. In Figure 15 and Figure 16 the calculated difference of the apparent depth using the two methods is presented.

The calculated depths of the UUI marked with the consecutive numbers 11, 12 and 5 and 7 add up to a total of 0.75 m, 0.82 m and 1.48 m according to the hyperbole fitting and 0.73 m, 0.79 m and 1.39 m, according to the layer velocity method. The difference is attributed to the small number of well-visible hyperboles, disregarding the diameter of UUI and the subjective assessment of the velocity obtained with the hyperbole fitting method. The heights of UUI defined with geodetic surveying methods were taken as a reference. Based on the depth deviations $\Delta h$, it can be argued that the method of layer velocity gives more accurate results based on the measured reference depth by geodetic surveying methods (see sub-Section The Establishment of a Geodetic Network for the Purpose of Determining the Position of the UUI and GPR).

5. Results

When interpreting the results obtained with the proposed kinematic GPR-TPS model, it is important to be familiar with the level of accuracy. This is assessed more accurately with reference measurements. The results of measurements using the geodetic polar method, which were carried out during installation of the UUI, were taken as a reference position and heights. The reliability of reference values is somewhat questionable because of the possible movement and subsidence of UUI and/or the medium during, as well as after, the installation. The pipes and cables are presented by line segments. They are bound by points determined by the polar method of detailed measurements. However, they at least offer some insight into the accuracy of the results that were obtained on the basis of the proposed GPR-TPS model. In the RUTP, 269 apex points of the installed UUI were recorded. The sample size is satisfactory as well as that for spatial distribution of the acquired GPR apex points of UUI. For the purpose of assessing the accuracy, an assumption has been made that the nearest apex point on the reference cylindrical object represents the point obtained by GPR-TPS model. In evaluating the results, the graphic displays of the position differences of apex points of UUI, histograms of distribution and Q–Q plots of coordinate and height deviations, were of immense help.

Great attention was paid to taking the measurements and the already tested and calibrated measurement equipment. Robust methods for the derivation of accuracy measures should, therefore, be applied if the differences are not normally distributed and/or contain skewness, kurtosis or an excessive number of outliers [56]. For this purpose, the tests of the normal distribution of both coordinate and height deviations within the accuracy assessment were consulted. The normal distribution of the positional and height differences was checked and analysed graphically/visually and by using statistical tests, the details of which can be found in Appendix B. If a normal distribution can be assumed and no outliers are present in the data set, differences of positional coordinates and heights obtained with the proposed GPR-TPS model are written as:

$$\Delta E_i=E_i-{E^{\prime }}_i$$

(11)

$$\Delta N_i=N_i-{N^{\prime }}_i$$

(12)

$$\Delta h_i=h_i-{h^{\prime }}_i$$

(13)

where ${E^{\prime }}_i$ and ${N^{\prime }}_i$ are the coordinates and ${h^{\prime }}_i$ is the height of the point $i$ obtained with the proposed kinematic GPR-TPS model, $E_i$ and $N_i$ are the coordinates and $h_i$ is the height of the reference point $i$, $\Delta E_i$ and $\Delta N_i$ are the coordinates and $\Delta h_i$ is the height difference from the reference data for point $i$.

5.1. Accuracy Assessment of the Recorded Position and Height of UUI with the Proposed GPR-TPS Model

According to the results of the testing distribution, the accuracy measure root mean square error (RMSE) and their positioning (radial) RMSE [21] was used in order to assess the accuracy of the position and height of the proposed GPR-TPS model:

$$RMS{E}_{\Delta E_i}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\Delta E_i^2}{n}}=0.037\mathrm{m}$$

(14)

$$RMS{E}_{\Delta N_i}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\Delta N_i^2}{n}}=0.066\mathrm{m}$$

(15)

$$RMS{E}_{\Delta h_i}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\Delta h_i^2}{n}}=0.115\mathrm{m}$$

(16)

$$RMS{E}_{pos}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\left(\Delta E_i^2+\Delta N_i^2\right)}{n}}=0.076\mathrm{m}$$

(17)

The analysis included a calculation of the central tendency and the dispersion of the obtained pairs of the coordinate and height differences. As a standard measure of central tendency, the arithmetic mean ${\mu }_{\Delta i}$, also called mean error (ME), of individual coordinate components was calculated. In assessing the dispersion of the coordinate differences, the standard deviation ${\sigma }_i$ of coordinate differences and the heights, where $n$ is the number of all tested points in the sample (sample size), was used:

$${\mu }_{\Delta E}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\Delta E_i=0.002\mathrm{m};{\sigma }_{\Delta E}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(\Delta E_i-{\mu }_{\Delta E}\right)}^2}{n-1}}=0.037\mathrm{m}$$

(18)

$${\mu }_{\Delta N}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\Delta N_i=0.012\mathrm{m};{\sigma }_{\Delta E}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(\Delta N_i-{\mu }_{\Delta N}\right)}^2}{n-1}}=0.065\mathrm{m}$$

(19)

$${\mu }_{\Delta h}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\Delta h_i=0.006\mathrm{m};{\sigma }_{\Delta E}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(\Delta h_i-{\mu }_{\Delta h}\right)}^2}{n-1}}=0.115\mathrm{m}$$

(20)

The positional (radial) standard deviation was calculated as:

$${\sigma }_{pos}=\sqrt{{\sigma }_{\Delta E}^2+{\sigma }_{\Delta N}^2}=0.075\mathrm{m}$$

(21)

Graphics dispersion and the distribution of points depending on the size of the positional and height differences are shown in Figure 17 and Figure 18. The estimated positional accuracy of the proposed GPR-TPS model in the RUTP was 0.076 m (positional RMSE), while the height accuracy was 0.115 m. The estimated central tendency of differences, the barycentre of difference vectors, lies 0.012 m from the origin in a northerly direction. The three times RMSE threshold, which provides an observation of gross errors, i.e., an error that was classified as an outlier if the coordinates or height differences > 3 RMSE [56], were not exceeded by any differences.

5.2. Assessment of the Incline Accuracy

An assessment of the incline accuracy of the UUI was carried out by comparing the reference incline measured with the reference geodetic method at the time of installation and the incline of the proposed GPR-TPS model. By using linear regression, the best-fitting straight line through height values and the horizontal length of the UUI (Figure 19) was found. Despite some low levels of the explained variance or the coefficient of determination $r^2$, which determine the low degree of linear-correlation of the length of the pipe or cable and the height obtained by the proposed GPR-TPS model, the inclines for assessing the accuracy were used.

To measure the inclination accuracy of the proposed GPR-TPS model, the RMSE, which was evaluated on the basis of the incline differences $\Delta n$, was used. The incline differences were only calculated for pipes, while the cables and cable ducts, where constant inclination cannot be assumed due to the installation procedure and their physical properties, were left out. Assessment of the inclination accuracy included determining the central tendency and the dispersion of the incline differences of the obtained pairs. As a standard measure of the central tendency of incline differences pairs, the arithmetic mean ${\mu }_{\Delta n}$ was used and the standard deviation ${\sigma }_{\Delta n}$ for the differences.

$$\Delta n_i=n_i-{n^{\prime }}_i$$

(22)

$$RMS{E}_{\Delta n}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n\Delta n_i^2}{n}}=1°{32}^{\prime }{07}^{″}$$

(23)

$${\mu }_{\Delta n}=\left|\frac{1}{n}{\displaystyle \sum }_{i=1}^n\Delta n_i\right|=0°{59}^{\prime }{10}^{″}$$

(24)

$${\sigma }_{\Delta n}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(\Delta n_i-{\mu }_{\Delta n}\right)}^2}{n-1}}=1°{12}^{\prime }{50}^{″}$$

(25)

where ${n^{\prime }}_i$ is the inclination angle of the UUIi captured by GPR-TPS model, $n_i$ is the inclination angle reference of UUIi, $\Delta n_i$ is the inclination angle differences of UUIi.

5.3. Graphic Presentation of the Final Results in the Real Urban Testing Pools (3D)

Figure 20 shows a 3D visualisation of the measured UUI in the RUTP with the proposed GPR-TPS model, based on 269 apex points of the installed infrastructure obtained by GPR data. The pipes based on the inherent apexes of the individual infrastructure with the use of linear regression are visualised. In phase I the lines and line segments through the points represented by the coordinates or numerical variables $N$, $E$ in a horizontal direction were determined. In phase II the regression line was determined that best fits the numerical variables ($h$—height, $d$—length), through which the inclinations and positions were defined in 3D. For visualisation of non-linear objects, a second-degree polynomial regression through the points presented by the coordinates $N$, $E$ in the horizontal direction was used. This was followed by linear regression ($h$—height, $d$—length), through which the endpoints height and, consequently, the inclination of the cables and ducts can be defined. Reference diameters were used for the purpose of visualisation. See Table 1 for a detailed description and the properties of the visualised UUI in the RUTP.

Figure 21 and Figure 22 show the position and height RMSE of the pipes in 3D. The reference UUI (green colour) are displayed with the consecutive numbers 9, 17 and 13 and obtained on the basis of the proposed GPR-TPS model (red colour) visualised by a regression analysis in the first testing pool.

6. Conclusions

The main achievement of the presented research is usage of the proposed kinematic GPR and self-tracking (robotic) terrestrial positioning system (GPR-TPS) model of data capture and processing. This proved that it is possible to recognise and record underground utility infrastructure (UUI) objects with the recommended/required high horizontal (positional) and vertical (height) accuracy demands of 10 cm. This model also introduced methods of 3D visualisation, which enable the spatial relationship of the UUI to be defined and provide a more illustrative and improved interpretation. Throughout the model, the advanced processing of radargrams is key in the recognition and identification of UUI, while latency plays a key role in determining the position with the kinematic GPR-TPS model with the use of pseudo-NMEA strings method.

Additionally, the usage of wireless communication standards (e.g., Bluetooth) in the model decreased crosstalk effects. The presentation and testing latency assessment using a simple GPR measuring technique of the same profiles in forward and reverse directions was performed. System latency evaluation depends on the components of the total composition or the hardware and software used. Evaluated latency was valid for the tested GPR-TPS configurations; in the event of replacement of any component from the configuration, latency can change dramatically, which has a significant impact on the final results. In this research, disregarding latency with a walking speed of 0.5 m/s leads to positioning errors in the range of 0.25 m, which can be crucial in determining the position of UUI. Based on the height (depth) differences $\Delta h$, it can be argued that the method of the layer velocity gives more accurate and realistic results in the testing pool environment than the method of hyperbole fitting. The difference between the method of hyperbole fitting and layer velocity was attributed to the small number of well-visible reflection hyperboles and disregards the diameter of UUI. All processes in the model were defined as an experimental and iterative way, which is inevitable in order to achieve good results. The processing flow of the GPR data was not fixed and depends on the GPR measurements that were determined by the physical properties of the medium and the objectives to be achieved.

A successive series of the proposed processing flow and parameters cannot be completely objective and therefore, according to some authors’ experience, is more like art than science [52]. The proposed kinematic GPR-TPS model mostly meets the recommended positional and height accuracy or at least gets close ($RMS{E}_{pos}=7.6$ cm; $RMS{E}_{\Delta h}=11.5$ cm). The model proved to be a useful tool in determining the inclination of UUI, which is restricted by the height accuracy of the method ($RMS{E}_{\Delta n}=1.5°$). When determining the inclination, the linear regression analysis into the model was introduced.

Questions and gaps arising with the interdisciplinary kinematic GPR-TPS model will have to be examined in the future in a real urban environment and improved and refined in terms of a more accurate estimation of the latency of the selected configuration. To achieve higher accuracy of the latency, the method of experimental laboratory latency of the measured composition using horizontal reference tracks, which would ensure the identical position of profile or pairs of radargrams (forward and reverse), should be used. The experimental laboratory latency evaluation would also improve the final result of the determination of position as well as filtering and smoothing of the data obtained in kinematic terrestrial processes. It would be reasonable to research the impact of the Kalman filter on the final results and include TPS, which would allow measurements to be captured with a greater frequency than 10 Hz. Such TPS will increase the density of the covered position and the height of GPR and, consequently, improve the positional and height accuracy of the UUI.

The proposed model was created for high-accuracy underground utilities mapping and only minor adaptations would be necessary for transfer to other fields of GPR investigation in urban environments. The most evident are applications for GPR surveys in arbitrarily oriented GPR transects for geotechnical engineering, protection of archaeological remains in urban contexts prior to the construction works and environmental studies with determination of the pollution due to the uncontrolled leakage and hidden waste landfills. Furthermore, this model enables high-accuracy GPR surveys on unevenly shaped surfaces. The exact heights of the GPR antenna position, obtained by the kinematic GPS-TPS model, enables accurate topographic correction of radargrams required for reliable estimation of depth and geometrical properties of diverse subsurface dielectrically contrast bodies.