Mapping and Assessment of Tree Roots Using Ground Penetrating Radar with Low-Cost GPS

In this paper, we have presented a methodology combining ground penetrating radar (GPR) and a low-cost GPS receiver for three-dimensional detection of tree roots. This research aims to provide an effective and affordable testing tool to assess the root system of a number of trees. For this purpose, a low-cost GPS receiver was used, which recorded the approximate position of each GPR track, collected with a 500 MHz RAMAC shielded antenna. A dedicated post-processing methodology based on the precise position of the satellite data, satellite clock offsets data, and a local reference Global Navigation Satellite System (GNSS) Earth Observation Network System (GEONET) Station close to the survey site was developed. Firstly, the positioning information of local GEONET stations was used to filter out the errors caused by satellite position error, satellite clock offset, and ionosphere. In addition, the advanced Kalman filter was designed to minimise receiver offset and the multipath error, in order to obtain a high precision position of each GPR track. Kirchhoff migration considering near-field effect was used to identify the three-dimensional distribution of the root. In a later stage, a novel processing scheme was used to detect and clearly map the coarse roots of the investigated tree. A successful case study is proposed, which supports the following premise: the current scheme is an affordable and accurate mapping method of the root system architecture.


Introduction
Environmental issues such as the conservation of natural heritage and ancient trees have become priority objectives of urgent protection [1,2]. Very limited knowledge is present on how the elements causing the rapid death of entire forests interact with each other. The quality and distribution of roots can be used as reliable diagnostic parameters for early tree decline [3][4][5]. Therefore, if the depth of the root and its extension to the surrounding area can be estimated, more effective preconditions can be provided to guide a wide range of research decisions [6].
Destructive methods can provide accurate mapping of the root systems. Nonetheless, intrusive approaches are often impossible to implement in the field and time-consuming. In addition, digging and trenching can affect the surrounding of trees and may cause irreversible damage [7]. Destructive methods are limited in space and can only provide local information on the tree conditions. Instead, non-destructive testing (NDT) methods are becoming more and more common in forestry and among other international venues and cultural facilities. The site was once the place where the U.S. military was stationed after World War II. After that, the Kawauchi Hagi Hall and the surrounding facilities were built in this area. The survey site was located in a rectangular area measuring 10 m × 4 m in the large open lawn, about 6-m beside the surrounding trees. The latter are typical Metasequoia glyptostroboides, with an average diameter at breast height of around 50 cm and a tree height that can approach 20 m. Figure 1a shows the scenario of the investigated site. Metasequoia glyptostroboides belong to a shallow-rooted species. The main roots of trees are underdeveloped, whereas the lateral roots or adventitious roots grown by radiation are much longer. Most of the roots are distributed on the surface of the soil. Lateral roots are mainly distributed between 0.2-1 m, with some tree roots being exposed to the ground. The white dashed line with an arrow in Figure 1a indicates the investigation direction crossing the root distribution that was used for the data collection stage.
For the purpose of this investigation, a 500 MHz RAMAC shielded bowtie antenna was employed. Frequency characteristics of the antenna is summarised in Table 1. Using a low-cost GlobalSat GPS receiver (Figure 1b) in combination with a GPR system, the positional coordinates of the moving antenna was able to be recorded in the collected GPR track. In order to obtain the exact location of the relevant GPR trace, the local GEONET station and the GNSS information were needed for post-processing of the position data collected by the low-cost GPR on the site.  Table 1. System parameters of the RAMAC 500 MHz Shielded Antenna system used for investigation purposes [27].

Parameters RAMAC 500 MHz Shielded Antenna
Limit frequency lower than -10 dB: Low In this study, to make the maximum amount of clear three-dimensional views of the tree roots and understand how the root distributed over the subsurface, high-density GPR data acquisitions are usually required. According to [27,30], the sampling rate Δ must be equal to or less than the Nyquist sampling rate ∆ along the direction of investigation and cross-investigation for the full-resolution 3D imaging: For the purpose of this investigation, a 500 MHz RAMAC shielded bowtie antenna was employed. Frequency characteristics of the antenna is summarised in Table 1. Using a low-cost GlobalSat GPS receiver (Figure 1b) in combination with a GPR system, the positional coordinates of the moving antenna was able to be recorded in the collected GPR track. In order to obtain the exact location of the relevant GPR trace, the local GEONET station and the GNSS information were needed for post-processing of the position data collected by the low-cost GPR on the site. In this study, to make the maximum amount of clear three-dimensional views of the tree roots and understand how the root distributed over the subsurface, high-density GPR data acquisitions are Remote Sens. 2020, 12, 1300 4 of 16 usually required. According to [27,30], the sampling rate ∆x must be equal to or less than the Nyquist sampling rate ∆x s along the direction of investigation and cross-investigation for the full-resolution 3D imaging: where θ is the main lobe angle of the antenna that was used in the investigation and λ is the wavelength of the subsurface material. Since the main lobe angle of the GPR antenna used in the near field is usually greater than 60 • , the maximum space interval should be equal to or less than one-quarter wavelength of the antenna bandwidth on the subsurface [27]. In order to evaluate the dielectric constant of the soil, the time-domain reflectometry (TDR) technique was used at several locations in the measurement area, as shown in Figure 2. The average value of the volumetric water content in the subsurface material was 42%, so that the subsurface could be classified as a wet soil. The relative permittivity measured near the soil surface was 27 which was calculated using the equation described in [31][32][33]. The corresponding subsurface velocity was 0.06 m/ns.
where is the main lobe angle of the antenna that was used in the investigation and is the wavelength of the subsurface material. Since the main lobe angle of the GPR antenna used in the near field is usually greater than 60°, the maximum space interval should be equal to or less than onequarter wavelength of the antenna bandwidth on the subsurface [27].
In order to evaluate the dielectric constant of the soil, the time-domain reflectometry (TDR) technique was used at several locations in the measurement area, as shown in Figure 2. The average value of the volumetric water content in the subsurface material was 42%, so that the subsurface could be classified as a wet soil. The relative permittivity measured near the soil surface was 27 which was calculated using the equation described in [31][32][33]. The corresponding subsurface velocity was 0.06 m/ns.
According to the parameters of the 500 MHz RAMAC shielded antenna in Table 1, the underground velocity is equal to 0.06 m/ns, and the spatial sampling interval is equal to or less than 5.5 cm. Therefore, the interval along the investigation direction should set less than 5.5 cm to meet the conditions of full-resolution 3D imaging. Under these circumstances, the system could provide a 3D subsurface image with 5.5 cm of horizontal resolution and 6.7 cm of vertical resolution in this survey scenario. Moreover, it is worth mentioning that an excess of 900 m of surveys with the GlobalSat GPS Receiver were performed in approximately half an hour. In particular, more than 90 surveys have been carried out across the 4-m cross survey direction. The raw radargram of the whole tree root survey is shown in Figure 3. A superposition of multiple responses could be seen up to 30 ns.   According to the parameters of the 500 MHz RAMAC shielded antenna in Table 1, the underground velocity v is equal to 0.06 m/ns, and the spatial sampling interval is equal to or less than 5.5 cm. Therefore, the interval along the investigation direction should set less than 5.5 cm to meet the conditions of full-resolution 3D imaging. Under these circumstances, the system could provide a 3D subsurface image with 5.5 cm of horizontal resolution and 6.7 cm of vertical resolution in this survey scenario. Moreover, it is worth mentioning that an excess of 900 m of surveys with the GlobalSat GPS Receiver were performed in approximately half an hour. In particular, more than 90 surveys have been carried out across the 4-m cross survey direction. The raw radargram of the whole tree root survey is shown in Figure 3. A superposition of multiple responses could be seen up to 30 ns.

Data Processing and Methodology
The processing presented in this study includes four sequential stages. In the pre-processing stage, time correction and a signal noise filter were applied to increase the clutter rate of the entire data. Then, by using the local GEONET station, Ultra-Rapid ephemeris, precise clock data, and an optimisation algorithm were able to obtain accurate location information, which was used. Next, the near-field Kirchhoff migration was used to map the distribution of the coarse roots. In the last stage, a novel processing scheme was used to identify and clearly map the coarse roots of the investigated tree. All the algorithms in each stage were implemented in MATLAB.

Radargram Signal Processing
Prior to any interpretation attempt, the raw data were subject to a pre-processing step in order to reduce undesired reflections and increase the overall signal-to-clutter ratio. In particular, the pre-processing step consisted of: (a) Time-zero adjustment This step is applied to adjust the initial position of the surface reflection in GPR signals (the time when the radar signal spreads out from the antenna and enters the ground is known as "zero time"). In order to move all traces to the zero-time position before other processing methods, time-zero adjustment must be applied.

(b) Background removal
Cross-coupling between the transmitter and the receiver as well as a ringing noise and multiple reflections can mask the less dominant reflections from the root. In order to mitigate such unwanted signals, background removal methods should be applied. Subtracting the background signal proved very effective for the investigated case study.

(c) Band-pass filtering
In order to improve the quality of the data, a band-pass filter is applied to eliminate different types of noise in this stage. In addition, filters are also used to extract useful information and hidden patterns from the recorded data. In order to remove the direct current offset and suppress the high-frequency noise, a band-pass filter was employed with the cut off frequencies F Low and F High indicated in Table 1.

Post-Processing of the Tracking Position
After pre-processing the GPR raw data, the post-processing will be performed on the GPR trace position data collected by the low-cost GPS receiver. So that the high precision position of each GPR trace can be applied for the 3D full resolution imaging for further investigation. GPS provides continuous global position, time and navigation services. The satellites send navigation and observation data to the receivers [34][35][36]. As the GPS receiver is affected by a number of different parameters, estimating the receiver position is not a trivial task. These parameters include GPS receiver clock bias, troposphere, ionosphere, multipath propagation and so on [37,38].
Each of these factors can be modelled separately. The basic pseudo-range equation is given as [39]: Remote Sens. 2020, 12, 1300 where c is the speed of light, T ε is the error produced by the troposphere, I ε is the error produced by the ionosphere and ε mul is the multipath error. ∆t s is the offset generated by the satellite bias and ∆t r is the receiver clock bias. ρ g is the geometric distance from the GPS satellite to the receiver which can be given by: where (x r , y r , z r ) and (x s , y s , z s ) are the positional coordinates of the receiver and a given satellite, respectively. In Equation (2), the error term of correction is divided into two groups. The troposphere error, the ionosphere error and the satellite clock bias are set into one group, represented as R s . The receiver clock bias and the multipath error are set into another group, represented as R r . Therefore, Equation (2) can be formulated as: In Equation (4), the receiver coordinates can be accurately calculated after correctly estimating the items R s and R r . By using signals from at least four satellites at the same time, the GlobalSat GPS receiver can provide real-time single point positioning [40]. The weakness of GPS comes from the satellite position, the satellite clock offset and the ionosphere influence, amongst others [41][42][43]. Nevertheless, the accuracy can be improved by post-processing the received signal.
Three to nine hours after the observation, the observed Ultra-Rapid ephemeris and clock data can be obtained online. This information is released four times a day at 03:00, 09:00, 15:00 and 21:00 UTC. The predicted Ultra-Rapid ephemeris and clock data are also released at the same time, which can be obtained in advance. Instead, information from Rapid 1 takes longer (from 17 to 41 h) to get online. It is released every day at 17:00 UTC [44][45][46][47][48].
International GNSS Service (IGS) is a cooperation with many organisations and countries. This collects data from more than 300 continuously operating reference stations around the world. Based on these data, it develops precise satellite ephemeris and clock solutions. The processing phase involves up to eight IGS analysis centers and the results are freely distributed by IGS. To this effect, this service allows researchers to post-process the observations based on the information provided by the IGS. Compared with a broadcast ephemeris and clock correction, the aforementioned approach is more accurate since these data reflect the precise position and clock offset of the satellite during the actual measurement. These data are classified into several categories and discussed in the literature [49][50][51].
In general, differential processing techniques depend on at least two receivers standing at control stations (known as a "bases") with a known location. In Japan, there exist around 1200 GEONET stations that record data every 30 s, as shown in Figure 4a. Therefore, it is possible to know the magnitude of the basic position error of the receiver. If the difference between the deviation of the reference station and the deviation of the flow station can be found, the position error at the reference line can be calculated. With the differential process, a correction is generated. The procedure can reduce the position error of unknown points. In general, this method can provide the location with sub-meter resolution from single-frequency pseudo-range observations. However, this accuracy level still cannot meet the requirements of current applications.
It is also worth noting that this application does not require a global solution. Only the 2D relative position of each GPR acquisition point needs to be estimated accurately. For this purpose, we selected a reference GEONET station (as shown in Figure 4b) 14 km away from the station. Assuming that the troposphere error T ε , the ionosphere error I ε and the satellite clock bias of the reference station were the same as those of the experiment site, the R s terms in Equation (4) can be directly estimated. Moreover, the receiver bias ρ g and the multipath error ε mul can be estimated and suppressed by designing an advanced extended Kalman filter.
Remote Sens. 2020, 12, 1300 where is the post-processed position of each acquisition point, is the recorded data at each GPR acquisition point, is the design matrix, which can be replaced by (as described in Figure  5a) and is the receiver clock bias and other minor correction terms such as the multipath error and the residual error.   In Equation (4), R s is mainly caused by the substantial errors of the positions of the satellites, the clock offset, the ionosphere correction etc. Using the record data from the reference station, a correction matrix was designed with the accuracy of the satellite's positions, clock offset and ionosphere information, which can be downloaded afterwards from the International GNSS Service (IGS) database in order to minimise R residual . The measurement equation could be written as: where P r is the post-processed position of the reference station, x re f , y re f , z re f are the recorded coordinates of the data at the reference station, H 0 is the design matrix that contains the satellite and the reference station receiver clock biases plus the ionospheric and tropospheric effects and another minor correction term. Lastly, R residual is the residual error. Based on the aforementioned, the least square term of the residual error R s in the experiment period can be written as: Consequently, the optimally-designed matrix H 0 can be calculated as: and Equation (4) can be written as: where P m is the post-processed position of each acquisition point, x is the recorded data at each GPR acquisition point, H is the design matrix, which can be replaced by H 0 (as described in Figure 5a) and R s is the receiver clock bias and other minor correction terms such as the multipath error and the residual error.   The R s term is estimated and removed based on the conventional Kalman forward-backward filter. The estimated offset of the receiver clock located at the tracking station does not represent the actual offset due to the fact that the obtained observation data were pre-processed [52][53][54]. An apriority bias of the GPS receiver was used in the beginning to estimate the term R s . Finally, all measurements were then corrected by the terms H and R s .
In order to update the filter over time, the predictive state vector of the system model should be used. Since the GPS receiver bias is assumed to be a linear function of time, the drift and all the other parameters can be considered as constants. It is worth noting that the drift is not strictly the same, as it will change slowly with time. Figure 5b depicts the complete flowchart of the proposed methodology consisting of six distinct and sequential steps:

(a) Forward Filter Initialisation
The coarse value of the GlobalSat GPS receiver is used as a prior value for the bias and the drift, with all the other elements set to zero. In addition, the noise of the filter needs to be set in this step.

(b) GPS Bias and Multi-Pass Error Constraint
In this step, the GPS bias and multi-pass error will propagate towards the current filter state.

(c) Kalman Filter Measurement Update
After the filter run is finished, a state vector computes the mean of the forward and backward results of the filter state.

(d) Results and Covariance Matrix
In this step, the results are stored to be used in the subsequent filter at this stage. In the meantime, the covariance matrix is needed to be evaluated. Then, the smoothing factor should be calculated.

(e) Judgement
The procedure is repeated until the processing reaches the end time, or until it meets the constraint condition. Here, Allen variance σ a is used to obtain a rough estimate of the expected receiver bias and residual error. The Allan variance σ a is a measure of frequency stability in clocks, oscillators and amplifiers [37]. After all the procedures are finished, the Kalman filter parameters are stored and used to design the smooth operator. Figure 6a shows the moving trajectories of the antennas recorded by the GlobalSat GPS receiver. During the investigation, the effects of the ionosphere had little to negligible effects to the overall precision of the algorithm. In fact, the drift of the GPS sensor itself was proven to have the biggest impact on the overall accuracy. The results after applying the Kalman filter and after converting the data to a local Cartesian coordinate system are shown in Figure 6b. The shape of the trajectory was consistent with the actual GPR system movement.

Kirchhoff Migration
In order to detect the three-dimensional distribution of the tree roots, a high-resolution and accurate imaging algorithm was implemented in this section. The underground 3D geometry can be reconstructed from the results of dense measurements. Three-dimensional migration is required to focus reflection and diffraction and collapse the recorded hyperbolas to their origin [55,56].
In this research, we applied the Kirchhoff migration, which is based on the generalised waveequation [57,58]. Mathematically, the Kirchhoff migration method can be formulated in the time domain as follows: where represents a source-free surface; is the unit vector normal to ; is the positional vector of the estimated wave field; is the positional vector of the integration point; is the propagation speed in media, represents the wave field within the surface and = .
Normally, the far-field approximation / ≫ 1 ⁄ is applied in (9). Therefore, the second term of Equation (9) can be omitted. Since we were oriented in the near-field target, such an assumption was no longer suitable for the purpose of this research. To mitigate that, the entirety of Equation (9) was used in the current study.
As described in Section 3, TDR measurements were conducted at different points within the surveying area. The average of the results turned out to provide a wet soil with a volumetric water content of 42% found using the Topp equation. The measured relative permittivity in the near-surface of the soil was 27 with a corresponding subsurface velocity of 0.06 m/ns. Since most of the roots were located in the shallow region of the subsurface, a constant velocity of the medium was assumed in this paper.

Kirchhoff Migration
In order to detect the three-dimensional distribution of the tree roots, a high-resolution and accurate imaging algorithm was implemented in this section. The underground 3D geometry can be reconstructed from the results of dense measurements. Three-dimensional migration is required to focus reflection and diffraction and collapse the recorded hyperbolas to their origin [55,56].
In this research, we applied the Kirchhoff migration, which is based on the generalised wave-equation [57,58]. Mathematically, the Kirchhoff migration method can be formulated in the time domain as follows: where S 0 represents a source-free surface; n is the unit vector normal to S 0 ; r is the positional vector of the estimated wave field; r 0 is the positional vector of the integration point; v is the propagation speed in media, U represents the wave field within the surface S 0 and Normally, the far-field approximation ω/Rv 1/R 2 is applied in (9). Therefore, the second term of Equation (9) can be omitted. Since we were oriented in the near-field target, such an assumption was no longer suitable for the purpose of this research. To mitigate that, the entirety of Equation (9) was used in the current study.
As described in Section 3, TDR measurements were conducted at different points within the surveying area. The average of the results turned out to provide a wet soil with a volumetric water content of 42% found using the Topp equation. The measured relative permittivity in the near-surface of the soil was 27 with a corresponding subsurface velocity of 0.06 m/ns. Since most of the roots were located in the shallow region of the subsurface, a constant velocity of the medium was assumed in this paper.

Results and Discussion
The previous section presented the processing methodology of our proposed strategy for tree root mapping and assessment. In this section, the results concerning the 3D distribution of the roots system from the field measurements will be discussed.

3D Migration Result
The 3D migration results were obtained by applying the Kirchhoff migration discussed above. One-twentieth of the resolution was applied as the grid size during the migration progressing. Figure 7 shows the migrated vertical profiles along the survey direction at different cross-survey directions

3D Migration Result
The 3D migration results were obtained by applying the Kirchhoff migration discussed above. One-twentieth of the resolution was applied as the grid size during the migration progressing. Figure  7 shows the migrated vertical profiles along the survey direction at different cross-survey directions  Based on the migration results, it was also obvious that most of the reflections from the roots were located in the upper region, as shown in Figure 8. In the figures, some areas with a strong amplitude showed the bright blocks feature. These bright blocks appear to be randomly distributed, and they can be mainly attributed to the heterogeneity of the soils, rocks and similar features in the subsurface. Hence, it is difficult to discriminate coarse roots by using GPR from these output types.
There are many reasons that limits the applicability and the accuracy of GPR for detecting tree roots. For example, if a cluster of fine roots surrounds a deeper coarse root, it is difficult to detect the coarse root. Moreover, if the distance between two coarse roots is less than the horizontal and vertical resolution of the observation system at the same time, then two coarse roots will be considered as a unique root. In addition, roots that are parallel to the polarisation of the antenna will result in weak reflections that will be masked by unwanted clutter and noise. Based on the aforementioned factors, the root can be detected according to the reflection intensity, direction and surrounding environment from the migrated slices. The distribution of tree roots must be mapped on a three-dimensional environment, as the reflection amplitude in the same depth slice varies along the root. Therefore, a coherent 3D approach is necessary in order to effectively locate and track the distribution of the roots.

3D Roots Detection Result
The root distribution can be tracked from the image index. The coarse tree roots normally can be displayed continuously in the 3D migrated data. Therefore, these image indexes can be associated with the root. In this paper, high-energy regions inside a 3D cube above a certain threshold were used to extract root locations. In Figure 10a, a B-scan profile which contains a tree root is shown. The profile is perpendicular to the root extension direction. In this migrated cube, a 1.5 m × 1.5 m × 0.3 m (survey direction × cross-survey direction × depth) block was selected. The peak position of the cube was found and a vertical curve (depth) and two horizontal curves (survey and cross-survey directions) of the peak were used for next step analysis. Additionally, the −3 dB width of the peak is defined as the root in each direction, as it is shown in Figure 10b,c.
Based on the migration results, it was also obvious that most of the reflections from the roots were located in the upper region, as shown in Figure 8. In the figures, some areas with a strong amplitude showed the bright blocks feature. These bright blocks appear to be randomly distributed, and they can be mainly attributed to the heterogeneity of the soils, rocks and similar features in the subsurface. Hence, it is difficult to discriminate coarse roots by using GPR from these output types.
There are many reasons that limits the applicability and the accuracy of GPR for detecting tree roots. For example, if a cluster of fine roots surrounds a deeper coarse root, it is difficult to detect the coarse root. Moreover, if the distance between two coarse roots is less than the horizontal and vertical resolution of the observation system at the same time, then two coarse roots will be considered as a unique root. In addition, roots that are parallel to the polarisation of the antenna will result in weak reflections that will be masked by unwanted clutter and noise. Based on the aforementioned factors, the root can be detected according to the reflection intensity, direction and surrounding environment from the migrated slices. The distribution of tree roots must be mapped on a three-dimensional environment, as the reflection amplitude in the same depth slice varies along the root. Therefore, a coherent 3D approach is necessary in order to effectively locate and track the distribution of the roots.

3D Roots Detection Result
The root distribution can be tracked from the image index. The coarse tree roots normally can be displayed continuously in the 3D migrated data. Therefore, these image indexes can be associated with the root. In this paper, high-energy regions inside a 3D cube above a certain threshold were used to extract root locations. In Figure 10a, a B-scan profile which contains a tree root is shown. The profile is perpendicular to the root extension direction. In this migrated cube, a 1.5 m × 1.5 m × 0.3 m (survey direction × cross-survey direction × depth) block was selected. The peak position of the cube was found and a vertical curve (depth) and two horizontal curves (survey and cross-survey directions) of the peak were used for next step analysis. Additionally, the −3 dB width of the peak is defined as the root in each direction, as it is shown in Figure 10b,c.
Overall, the accuracy of the root detection using pixel intensities in the image was slightly worse than extracting the feature directly from the waveform. The -3-dB threshold did not remove all the weak reflections, such as those from the heterogeneous soils, rocks, and fine roots. To avoid false detection results, future processing should also include discarding those targets that do not have continuous distribution.  Overall, the accuracy of the root detection using pixel intensities in the image was slightly worse than extracting the feature directly from the waveform. The -3-dB threshold did not remove all the weak reflections, such as those from the heterogeneous soils, rocks, and fine roots. To avoid false detection results, future processing should also include discarding those targets that do not have continuous distribution.
The reconstructed 3D root system is illustrated from two different viewpoints in Figure 11. Six coarse roots across the survey direction can be clearly detected. The false detection results were discarded. Root 1, 2, 3 and 4 were located at the 30-cm region below the subsurface and their distributions were confirmed with a thin measuring flower pole, which can be inserted into the soil and record the depth of these targets. Root 5 and 6 were located at 0.5 m and 0.7 m, respectively. In the field data set, it was impossible to detect all the coarse roots. The characteristics of tree roots radiating along the trunk provided an assumption for identifying the extension of the roots. On the contrary, fine roots always gather on heterogeneous soil and stone. This is because it is not conducive to the growth of trachyte in these areas. Since there was no marked difference between the fine roots and the surrounding environment, it was difficult to find these targets.
Remote Sens. 2020, 01, x FOR PEER REVIEW 13 of 16 The reconstructed 3D root system is illustrated from two different viewpoints in Figure 11. Six coarse roots across the survey direction can be clearly detected. The false detection results were discarded. Root 1, 2, 3 and 4 were located at the 30-cm region below the subsurface and their distributions were confirmed with a thin measuring flower pole, which can be inserted into the soil and record the depth of these targets. Root 5 and 6 were located at 0.5 m and 0.7 m, respectively. In the field data set, it was impossible to detect all the coarse roots. The characteristics of tree roots radiating along the trunk provided an assumption for identifying the extension of the roots. On the contrary, fine roots always gather on heterogeneous soil and stone. This is because it is not conducive to the growth of trachyte in these areas. Since there was no marked difference between the fine roots and the surrounding environment, it was difficult to find these targets.

Conclusions
In this paper, a low-cost GPS receiver combined with a 500 MHz RAMAC shielded antenna was used to obtain a full-resolution 3D image of the root system of several larch trees. By using the lowcost GPS receiver and advanced processing algorithms, an effective and affordable way for mapping tree root systems was proposed. This equipment combination was used to carry out a full-resolution survey over a 10 m × 4 m area in approximately half an hour, which proved the time efficiency of the approach proposed. In addition, a dedicated post-processing methodology was implemented to enhance the accuracy of each GPR trace recorded by the low-cost GPS receiver. After that, a Kirchhoff migration approach considering the near field effects was applied to accurately map the investigated root system. Lastly, the 3D distribution of coarse roots was mapped by considering the local magnitude in a 3D migrated cube. As a result, six large coarse roots were successfully identified. Although it is not possible to have full detection of the entire root system (i.e., coarse and fine roots), Figure 11. The reconstructed 3D root system: (a) View from the starting point; (b) view from the ending point.

Conclusions
In this paper, a low-cost GPS receiver combined with a 500 MHz RAMAC shielded antenna was used to obtain a full-resolution 3D image of the root system of several larch trees. By using the low-cost GPS receiver and advanced processing algorithms, an effective and affordable way for mapping tree root systems was proposed. This equipment combination was used to carry out a full-resolution survey over a 10 m × 4 m area in approximately half an hour, which proved the time efficiency of the approach proposed. In addition, a dedicated post-processing methodology was implemented to enhance the accuracy of each GPR trace recorded by the low-cost GPS receiver. After that, a Kirchhoff migration approach considering the near field effects was applied to accurately map the investigated root system. Lastly, the 3D distribution of coarse roots was mapped by considering the local magnitude in a 3D migrated cube. As a result, six large coarse roots were successfully identified. Although it is not possible to have full detection of the entire root system (i.e., coarse and fine roots), this research demonstrates the potential of combining low-cost GPS and GPR antenna systems to achieve a full-resolution 3D GPR survey in a time-effective and affordable manner.
Author Contributions: L.Z. analysed the data, proposed the algorithm and wrote the paper. Y.W. contributed to data analysis and provided useful suggestions. I.G., F.T. and A.M.A. contributed in structuring the focus of the paper and the presentation of the results, as well as in editing the language. M.S. contributed much to the experimental design, data collection, and language correction. All authors have read and agreed to the published version of the manuscript.