Kriging-Based Robotic Exploration for Soil Moisture Mapping Using a Cosmic-Ray Sensor

Soil moisture monitoring is a fundamental process to enhance agricultural outcomes and to protect the environment. The traditional methods for measuring moisture content in soil are laborious and expensive, and therefore there is a growing interest in developing sensors and technologies which can reduce the effort and costs. In this work, we propose to use an autonomous mobile robot equipped with a state-of-the-art non-contact soil moisture sensor that builds moisture maps on the fly and automatically selects the most optimal sampling locations. The robot is guided by an autonomous exploration strategy driven by the quality of the soil moisture model which indicates areas of the field where the information is less precise. The sensor model follows the Poisson distribution and we demonstrate how to integrate such measurements into the kriging framework. We also investigate a range of different exploration strategies and assess their usefulness through a set of evaluation experiments based on real soil moisture data collected from two different fields. We demonstrate the benefits of using the adaptive measurement interval and adaptive sampling strategies for building better quality soil moisture models. The presented method is general and can be applied to other scenarios where the measured phenomena directly affects the acquisition time and needs to be spatially mapped.

into the kriging framework. We also investigate a range of different exploration strategies and assess their usefulness through a set of evaluation experiments based on real soil moisture data collected from two different fields. We demonstrate the benefits of using the adaptive measurement interval and adaptive sampling strategies for building better quality soil moisture models. The presented method is general and can be applied to other scenarios where the measured phenomena directly affects the acquisition time and needs to be spatially mapped.

INTRODUCTION
Management of water resources is of considerable concern in different parts of the world, with many areas facing prolonged droughts, while others experience devastating floods. The availability of water in the soil is essential for vegetation. In an agricultural setting, crop health depends greatly on soil moisture. It is precisely for this reason that soil moisture monitoring is key to improving agricultural processes. Perhaps the most obvious advantage of technologies for obtaining high-resolution soil moisture maps is that they would enable highly efficient irrigation planning, for example, providing an accurate estimate of the quantity of water that should be put into a field and its required spatial distribution across the field.
Soil moisture is typically assessed either by a direct but lengthy procedure involving collecting physical soil samples followed by lab measurements, or by hand-held instruments used to measure moisture indirectly through proxies such as surface tension (manometers), or changes in soil conductivity (e.g. time-domain reflectometry (TDR) [1]. All of these methods are very laborious, time consuming and expensive. Recent advances in sensing technology introduced a new, non-contact method for measuring soil moisture using fast neutron detectors ( [2]). The neutrons are generated by cosmic rays and are reflected from the soil. The reflected neutron count is directly proportional to soil moisture content. Such sensors were successfully deployed at static locations covering large areas of land [3] but also as high-resolution variants with reduced field of view and increased sensitivity [4].
The most common method for creating soil moisture maps is to use data that are manually collected at pre-determined locations in the field and extrapolate the expected measurements for unvisited regions using kriging or Gaussian Process Regression [5,6]. This is a costly and laborious process, especially in the case of soil moisture monitoring, where the methods and instruments used to take measurements across the field require a high amount of labour and post-processing. For this reason there is a growing interest in developing instruments and methodologies to help reduce the effort and costs, while improving the quality of the resulting soil moisture models.
In this work, we propose to use an autonomous mobile robot equipped with a non-contact soil moisture sensor that builds soil moisture maps on the fly and automatically selects the most optimal sampling locations. The robot is guided by an autonomous exploration strategy driven by the quality of the soil moisture model (i.e. Kriging Variance) which indicates areas of the field where the information is less precise, improving overall model quality. The employed fast neutron counting sensors provide a special category of measurements in which the acquisition time directly depends on the intensity of the phenomenon: in our case, the sensor registers more neutrons in drier soils. We model the sensor using the Poisson distribution and use a special kriging variant for this type of measurements. As a result, the exploration strategy plans not only the optimal sampling location but also the required acquisition time at each sampling location.
The contributions of this work are as follows: • application of a novel fast neutron counting sensor for robotic-assisted spatial mapping of soil moisture; • integration of the Poisson measurement model into the kriging estimation and exploration framework, which devises optimal spatial locations and measurement intervals, improving the resulting moisture models; • evaluation and validation of the proposed framework on data collected from two different field environments.
The remainder of the paper is structured as follows: Section 2 presents related work in soil moisture surveying and robotic exploration, followed by Section 3, which details our approach to Poisson kriging and exploration for soil moisture mapping using a mobile robot. The experimental framework is presented in Section 4, followed by results and their analysis in Section 5, and final conclusions in Section 6.

RELATED WORK
Robotic environmental monitoring applications have attracted a lot of attention in the last few years [7]. One of the advantages of using robots for environmental modelling and monitoring is that they can build models on the fly. At the same time, many authors have discussed how to use the model itself to plan new observations for data acquisition that improve the overall model. For example, Kerry et al. [8] demonstrated that kriging semivariograms are highly useful for sampling planning in precision agriculture. They proposed to use ancillary information to estimate a semivariogram and thus determine the spatial frequency of sampling based on the semivariogram parameters.
Other researchers [9] propose the generation of an initial set of samples to obtain a semivariogram that can be extrapolated to find new sample positions. Marchant & Lark [10] proposed an adaptive approach for optimizing reconnaissance surveys. They sampled at pre-planned positions, and calculated the probability density function of the sampling density required for the main survey in a Bayesian framework. If the requirements were not met, the number and location of observations within further phases were selected to reduce the uncertainty of the required sampling density. However, the effort required to survey a soil variable and simultaneously build and analyse the variance of the kriging model of the soil meant that these authors stopped short of planning the entire sampling procedure based on kriging models.
Robots, on the other hand, are able to create and update models of their operational environments through robotic exploration. A common approach is to plan trajectories that completely cover the area assuming some prior knowledge of the environment [11]. Other well-known exploration techniques drive the robot towards unmapped areas of the environment. For example, greedy approaches such as [12] drive the robot towards the nearest location where new information can be gained. In frontier-based exploration [13], the robot is driven towards the boundary between the known and unknown parts of the environment, while information driven 'next-bestview' methods use reward functions to predict the utility of an unexplored location [14].
Many authors have proposed informative path planning (IPP) techniques for modelling physical phenomena with an unknown spatial distribution. These techniques address how to plan a path that maximizes sensor information [15] and can be classified into two approaches: those that depend only on a priori information about the environment [16], and adaptive sampling techniques that can be modified depending on the observations made [17]. More recently, Popovic et al. [18] proposed an adaptive informative path planning methodology to map green biomass in an agricultural setting.
Other authors have opted to use different model properties to plan robot actions. For example, Gao et al. [19] propose the use of an informative sampling technique to minimize the total distance travelled by a fleet of phenotyping robots. To do this, they model the environment using Gaussian Processes and use the model variance to plan the most informative paths for the fleet. Marchant and Ramos [20] use Gaussian Processes to plan the paths that guarantee both to observe the phenomenon of interest and improve the modelling of the same phenomenon for environmental monitoring applications such as ozone concentration across the USA.
Other authors have chosen to use Ordinary Kriging to model in-field phenomena. Glaser et al. [21] use it to model soil properties perceived with a multi-spectral camera, and then use the resulting model to improve the robot localisation. Kim and Shell [22] proposed an augmentation of Ordinary Kriging to enable modelling of ocean current dynamics which they use for adaptive path planning in the field in ocean multi-robot scenarios. Pulido Fentanes et al. [23] proposed a robotic exploration methodology aimed at building soil condition maps using ordinary kriging variance as a reward function for exploration. The current work builds upon this approach to model soil moisture measured with a novel sensor that does not follow a normal distribution. To achieve this we combine Poisson kriging with a kriging-based exploration methodology.

METHODOLOGY
In this work, we propose a kriging-based exploration pipeline for agricultural mobile robots to facilitate efficient mapping of soil moisture. The framework combines a unique sensor model, an on-line spatial mapping component and an exploration strategy to guide the robot to the next best sampling location.
We consider a special category of measurements which are based on counting, and hence follow a Poisson distribution. An inherent property of such measurements is that their uncertainty directly depends on the length of the measurement interval. In our scenario, we use a robot-mounted soil moisture sensor (see Sec. 3.1) which counts low energy neutrons as a proxy for soil moisture. Therefore the soil moisture level will affect the amount of time the robot spends at each sampling location. For the spatial mapping we use a version of ordinary kriging which incorporates measurements following a Poisson distribution (see Sec. 3.3). We use the Kriging Variance (KV) as a reward function for the exploration strategy to plan the optimal location for each subsequent measurement. Section 3.4 discusses the different exploration strategies that have been applied in this work.
The original kriging framework was presented in our previous work for mapping soil compaction [23]. In this paper, we generalise and extend the approach to take into account measurements following a Poisson distribution. This results in exploration strategies which not only consider the optimal sampling location but also adjust the measurement duration for each reading to ensure a high-quality model.

Soil Moisture Measurement Using a Cosmic-Ray Sensor
The main sensor used in this work is based on measuring fast neutrons, which are generated by cosmic rays and reflected from the soil [2]. The intensity of the reflected neutrons is affected by the hydrogen in the soil, and hence provides an indication of the soil moisture content. A neutron detector is a tube containing a gas that can convert thermal neutrons into detectable electrons by ionisation. Since the detectors are sensitive to fast neutrons only, the low energy neutrons (after colliding with the hydrogen atoms) are not counted. As a result, a higher neutron count means more fast neutrons and corresponds to dryer soil. To improve the sensitivity of the detector to fast neutrons, a polyethylene shield is used as a moderator.
Several correction procedures need to be applied on the acquired neutron counts (which we refer to as the raw neutron count N raw ) in order to account for variations in background cosmic ray intensity, atmospheric pressure and humidity [3]. The reference values for the corrections are established during a calibration procedure which requires reference soil moisture values to be established by direct soil moisture measurements using traditional equipment. The correction factors include: • Cosmic ray intensity: where C is the measured neutron count rate (from the nearest monitoring station) and C 0 is the value measured during calibration.
• Pressure: where P is the measured barometric pressure (using a barometer), P 0 is an arbitrary reference value (e.g. 1010 hPA) and β is the barometric pressure coefficient established during calibration.
• Humidity where Q is the measured humidity (derived from temperature measurements) and Q 0 is the average humidity during calibration.
The corrected neutron count N crr is obtained by multiplying the raw neutron counts by the correction factors: N crr can then be used to calculate Volumetric Water Content (VWC), which provides the final measure of the soil moisture. Since in this paper we mainly work with the corrected neutron counts N crr , we refer the interested readers to [3] for further detail of the exact conversion procedure.
The summarised methodology for measuring soil moisture has been used successfully by [3], who have established a network of soil moisture monitoring stations in the UK covering an area of 12 hectares. Although this coverage is useful for large scale soil moisture assessment, its application to individual fields in agriculture is limited. To achieve higher spatial resolutions, we have employed a high-sensitivity version of the sensor consisting of 12 neutron detectors with a bespoke polyethylene shield to limit the detection footprint of the sensor to ∼10 m. The sensor mounted on our agricultural mobile robot Thorvald can be seen in Fig. 3.

Poisson Distribution Measurements and Sampling Regime
Our soil moisture sensor provides the corrected neutron counts N crr . The appropriate probabilistic model for modelling count data and events is the Poisson distribution, with parameter λ representing the average count rate over a period of ten seconds. However the uncertainty σ in the measurement depends directly on total neutron count over the measurement time, and is calculated as follows: Figure 1 shows the histogram reading for a sensor measurement and the evolution of the λ and σ parameters for the same measurement over time. Figure 1b shows how the standard error and variance decrease over time, meaning that readings with longer duration achieve higher quality.
The sampling regime is the criterion used to decide how long each measurement should last. In this scenario, the quality of the measurement is directly correlated to the number total number of observed events (N corr ). For this reason, we propose to use two different methodologies:  The Adaptive Measurement Intervals (AMI) regime uses a threshold typically defined in terms of σ m (see Eq. 5) to determine the duration of a measurement. In practice, this means that in this case the robot will stay at each location until the normalised standard error falls below a pre-determined percentage of the total amount of counts, so that the robot will stay longer in places were the count rates are lower (or the soil is wetter in this scenario) and spend less time in locations with higher count rates. Figure 2 illustrates the evolution of the normalised standard error (σ m ) over time for different rates (λ), where the dashed lines indicate thresholds that can be used for this sampling regime, the time at which the threshold lines intersect the standard error lines, represents the point at which the measurement is considered complete. This guarantees a maximum incertitude limit for each measurement which adapts to the actual neutron rate forcing the robot to stay longer at places where the rate of events is lower than usual or to leave as soon as possible in places with higher rates.

Poisson-Kriging
Ordinary kriging (OK) has proven to be an effective method for interpolating spatial data when the data's main source of error is intrinsic to the measurement technique, for example, when it depends on the precision of an instrument. However, when the variance of the measurement depends on the phenomenon itself, as in the case of events that can be modelled using a Poisson distribution, Ordinary Kriging does not have a way to incorporate the different variances from each data point.
For this reason, different authors have proposed specific implementations of kriging methods that deal with data that is not normally distributed. Monestiez et al. [24] presented a kriging methodology to model whale populations using data from observers on ferries and cargo ships, which can be modelled using a poisson distribution. This approach is known as Poisson-Kriging (PK) and has since been used to model phenomena as diverse as Cancer mortality [25] and gamma-ray spectral mapping [26]. For this reason, we have chosen this methodology for the current work.
PK provides an estimateẐ(x 0 ) for a variable Z at unknown location x 0 while assuming a constant unknown mean over its neighbourhood, although in this case the observations Z(x i ) are dependent on some underlying mean count rate and the amount of time spent at each location. The estimate is a weighted linear combination of the available observation z i = Z(x i ) and the amount of time spend at each location t i from a set of locations x i . The estimator is thus described as follows:Ẑ where n i=1 w i = 1 to ensure unbiased estimates. To correctly estimate the values at x 0 the weights w = [w 1 , . . . , w n ] T must be calculated. This can be achieved by solving the Poisson-Kriging system, which is a linear system of n + 1 equations. n j=0 w ij C ij + w im t i + µ = C ix 0 f or i = 1, ..., n where C ij is the covariance of the observed values, C ix 0 is the covariance at the prediction location x 0 , and µ is a Lagrange factor which ensures the optimal solution. Finally,m is estimated from the data as a weighted average of the count rates, where the weights correspond to the observation times.
Once this system is solved, the estimated values at location x 0 can be found using Eq. 6, and the associated variance of the prediction σ 2 can be calculated using the same equation as in ordinary kriging:

Semivariogram
In empirical scenarios, it is possible to use a semivariogram created from the real-world data to express the relation between locations and estimate the weights for each observation. However, unlike Ordinary Kriging in this case it is necessary to account also for the observation times for each data point. For this reason PK uses a weighted variogram estimator, which takes into account the different observation times: where h is the chosen distance,m is the same mean as in Eq. 7 and I d ij ∼h is a gating function that takes a value of 1 when i and j are roughly distance h apart, and 0 otherwise. N (h) is a normalising factor calculated as follows: The semivariograms γ(h) can take multiple forms, but are generally characterised by an equation that can be parametrised. We use the following Gaussian semivariogram model in our work: with the following three parameters: nugget p 0 , range p 1 and sill p 2 [27]. The parameters for this equation are automatically fitted from the semivariogram of the sampled data using the soft L 1 norm minimization scheme [28].

Exploration Strategies
Our proposal is to use the variance of the kriging (KV) process (see Eq. 8) as a measurement of information gain. The use of Kriging Variance as a reward function for robotic exploration has been previously studied in [23,27]. In this work, we compare some well-known exploration strategies and how they interact with the sampling regime. The methods to be tested can be classified into: Next-Best-View methods and Adaptive Sampling methods.

Next-Best-View (NBV)
These methods update the environment model every time a new sample is acquired, and then choose a new location depending on the distribution of the KV across the field. Location selection is done using one of the following strategies: • Greedy: The next sampling point is the point with the highest KV in the set of candidate locations.
• Monte Carlo: a set of candidate sampling locations is generated each time, and each candidate location is allocated a weight depending on its KV. The next sampling location is selected randomly, but in a way that guarantees that the probabilities are distributed according to the weight of each candidate.

Adaptive Sampling
In this category, strategies generate an initial plan that is modified depending on the KV after each model update. In this case, the robot will plan a sampling regime based on a random trajectory and a mission time horizon, which depends on the minimum expectations of measurements to be made in each case. Every new sample taken triggers a model update, which removes sampling points based on their KV. The targets whose KV is below the overall mean of the model are removed, and then as many new points as necessary to meet the minimum expectation of measurements in the remaining mission time are added by drawing a new waypoint from a set of candidates weighted by their KV. Finally, a new route is re-planned through the new set of points using a Traveling Salesperson (tsp) algorithm.

Hardware Setup
Our experimental set-up consists of an autonomous mobile robot Thorvald [29] equipped with a custom-made, high sensitivity soil moisture sensor based on fast neutron counting principle manufactured by Hydroinnova (see Fig. 3). The 12 neutron detectors are accompanied by temperature and humidity sensors which are used for providing the corrected neutron counts every 10 s. The sides and top of the sensor are shielded by using a 50 mm polyethylene shield to limit the detection footprint of the sensor to 10 m. The total weight of the sensor is around 300 kg. The sensor is interfaced with the robot through an Ethernet link. The robot is controlled through an in-built PC running Linux OS and Robot Operating System (ROS). The platform is equipped with a GNSS sensor, which enables robot localisation and geo-tagging of the collected data samples. The navigation component uses a graph-based representation, allowing the robot to move between a pre-determined set of waypoints.

Datasets
Evaluating the performance of robotic exploration strategies is inherently difficult and previous work in that domain often relies on simulated experiments (e.g. [30]). In our case, we propose to use the 'surrogate' models of soil moisture, based on data collected from two real fields with the described equipment. We used the collected data in off-line 'simulations' to compare different exploration strategies and understand their overall performance. Simulations using a surrogate model are a useful tool to compare exploration methods [14,23], providing the 'ground truth' for the exploration results.
The two data collection sites include an airfield at the Lincolnshire Aviation Heritage Centre in East Kirkby, UK and a wheat stubble field near Volos, Greece. Both fields were prepared in such a way so that they had equal parts of dry and wet land. Such an arrangement enabled us to systematically test the effectiveness of kriging-driven exploration strategies under a significant gradient between dry and wet areas akin to a step response.
The airfield site (see Fig. 4) features a hard border between the grass field and concrete airstrip. Since concrete contains low levels of hydrogen, the airstrip provides a perfect replacement for dry conditions (5% Volumetric Water Content (VWC)). The data collection took part in March 2018 and therefore the grass field was in a relatively wet condition (20% VWC). 13 The wheat stubble field in Greece (see Fig. 4) covered a rectangular area of approx. 7 ha. The data collection took part in June 2018 under dry weather conditions. To create a wet area, the field was irrigated prior to data collection resulting in a wet/dry border with VWC of 18 % for the dry part and 24 % for the wet area, representing a fairly low gradient between the two parts. The whole field was meshed into a grid of 72 sampling locations with a spatial resolution of 30 × 30 m. The measurement interval for all the points was set to 10 min.
Both datasets were used to create a set of testing models which were used to verify multiple hypotheses presented in Sec. 5. Each one of this testing models has neutron rates as inputs for the measurement model which were then extrapolated across the testing area using Ordinary Kriging (OK). This way an estimated rate can be produced for every location on the field. Once this is done the extrapolated rate is used as λ to produce simulated readings every 10 s (real sensor's update rate) at specific locations using a Poisson distribution, resulting in a high density models used as a reference. The models include: • Synthetic model which is based on real sensor rates recorded from the airfield (see Fig. 5a).
We generated two models representing a high and low gradient between sensor rates for wet and dry soil respectively. The high gradient corresponds to rates of 2.5 and 5.0 counts/s for the wet and dry part respectively. For the low gradient the values are 3.0 and 4.0 counts/s respectively.
• Simulated model is based on the real data recorded in the airfield and extrapolated into multiple lines covering a rectangular area. The rates recorded along the single line crossing the wet/dry border are used to generate additional 5 parallel lines 10 m apart.
• Validation model in which the real data from the wheat stubble field is used (see Fig. 5b). This model represents the most realistic soil moisture conditions and is used to validate the proposed algorithms.
(a) (b) Figure 5: The high gradient synthetic model generated from the airfield (a) and the validation model generated from the wheat stubble field (b).

EXPERIMENTS
To evaluate our framework, we have devised a set of experiments to test multiple hypotheses. First, that the robot will focus on sampling the area with the highest uncertainty, i.e. the border between the soil and concrete parts of the field and borders of the field. Second, we want to verify how much does the rate difference between the wet/dry parts of the field influence the exploration process (we call this a step response). Finally, we want to analyse the different impact of having a Fixed Measurement Interval (FMI) and an Adaptive Measurement Interval (AMI) which warrants a minimum measurement uncertainty before moving on to the next sampling point. Because our sensor follows the Poisson distribution model, we believe that the robot will require less time to sample the dry area of the field as it would have observed a higher number of events in the same time reducing the measurement σ.
The results presented in this section were obtained using simulated runs over the testing models presented in Section 4.2. The performance of the exploration methods presented in this section is evaluated in terms of travelled distance and model error. For assessing the quality of the resulting model, we compare the model produced against the surrogate model used for the exploration. To compare any two resulting models A and B we use Mean Square Error (MSE):

Fixed vs Adaptive Measurement Interval
To compare the influence that the sampling regime has over the exploration process, a greedy strategy was tested in the synthetic experimental set-up following four different sampling regimes, two Fixed Measurement interval (FMI) and two Adaptive Measurement Interval (AMI) experiments. For the FMI case one experiment was set to 10 minute intervals (FMI-long) and the other one to 5 minute intervals (FMI-short). For the AMI case one experiment was set to a 2.5% measurement σ threshold (AMI-long) and the other one to a 3% threshold (AMI-short). Short and long cases should have comparable measurement times between them.  Figure 6: High-gradient synthetic scenario: comparison of performance for methods using Fixed vs Adaptive Measurement Interval in terms of (a) travel distance, and (b) Mean Square Error. Average results over ten runs, coloured areas represent standard deviation for each case. Figure 6a shows that distance is driven mainly by measurement time. This was predictable given that the amount of time that the robot spends collecting data is inversely proportional to the amount of time the robot spends navigating from one location to another. In figure 6b it can be seen that AMI regimes lead to faster convergence than their FMI counterparts.
Adaptive Measurement Interval strategies achieve better quality in shorter times because they can optimise the sampling time and drive exploration considering the conditions of the field (for example, the robot will spend less time in drier places as it will observe a higher number of events and achieve higher levels of confidence for the readings). These gains are highly dependant on the variability of the soil moisture in the field, for example, in a predominantly wet field the gains from adaptive sampling interval strategies will be less noticeable. To verify this hypothesis, this analysis was also performed on a simulation with lower gradient between the wet and dry parts. Figure 7 shows a comparison of the performance of both regimes in scenarios with different gradients. From this figure it can be seen that the difference in performance between both regimes is not as noticeable in the scenario with the lower gradient. However, the travelled distance for the Adaptive regime is slightly higher in both cases, indicating that sampling regimes are not important for controlling the travelled distance, and that this is a factor that is probably driven by the exploration strategy.  Figure 8, presents a comparison between the performance of both sampling regimes in the simulated model. Comparing these results to the ones obtained with the synthetic model (Sec-tion 5.1), it is possible to see that the results are almost identical for both cases. This indicates that despite the fact that variability is just slightly higher than the lower gradient synthetic scenario, the sampling regime does have an influence, this could indicate that it is preferable to use an AMI regime in every case as it becomes influential with medium gradients but its cost in travel distance is not much higher.

Comparison of the Exploration Strategies
To verify the influence of different exploration strategies over the exploration process, we ran a series of simulations with 3 different strategies namely: Greedy, Monte Carlo and Adaptive Sampling. In all cases we used AMI as the measurement interval regime to isolate the effects of the exploration strategy only.  Figure 9: Synthetic airfield scenario: performance for different strategies using Adaptive Measurement Intervals in terms of (a) distance, and (b) Mean Square Error. Coloured areas represent standard deviation over ten runs. Figure 9 shows the performance of the different exploration strategies. From these results, it can be seen that exploration strategies have a high influence on the distance travelled by the robot. In particular, it can be noticed that adaptive sampling strategies can achieve models that are just slightly worse than the other two strategies but travel much lower distances. This is a very important consideration, because in larger fields long travel distances can translate into significant amount of time not spent on gathering data, which at the end might end up with the lower quality models. Figure 10 presents the outputs of the exploration result for the high-gradient synthetic model. The figure shows the resulting models for a field after two hours of autonomous exploration with the trajectories followed by the robot. One interesting thing is that the greedy strategy drives the robot mostly to the edges of the field. This is mainly because the kriging methods are better at interpolation than extrapolation, so the highest variances are always around the limit areas. This has the advantage that it can drive the model's variance down very quickly. It might also mean, however, that it can miss relevant infield information. In comparison, the adaptive sampling followed a much smoother and shorter trajectory and took samples that were better distributed across the field.
To verify this findings we performed the same test on the simulate airfield scenario. Figure  11 shows that the performance of the different strategies is similar to that exhibited in Figure 9. This indicates that the behaviour of each strategy is consistent and does not tend to vary much across testing scenarios.
The outputs (see Figure 12) show again that greedy strategies follow very long paths and outer sampling points contrasting to the adaptive sampling method which follows a more balanced approach, that seems to linger around areas that are either drier or wetter than usual. The Monte Carlo approach shows an interesting behaviour, it appears to be going back and forwards around the border between the grass and concrete, this seems to be because there higher vari-  ances around the border area however the paths are very random and the travelling distance is heavily penalised. In that sense the adaptive sampling strategy has a big advantage over Monte Carlo because it follows the same principle for choosing targets but at the same time it reduces travel distance.

Validation on the Surrogate Model
To validate the methodology, several experiments were executed simulating an exploration task of four hours. Figure 13 presents the resulting models for three experiments using different exploration strategies and AMI as sampling regime. It is possible to see by simple visual inspection that the resulting models do not reflect perfectly the validation model. We believe that this is mainly due to two factors: first ,the gradient between wet and dry parts in this environment was very low which made it hard for the methods to identify areas of high variance. And second, the size of the environment limits how many samples per hectare the robot can achieve. This means in practice that the maps had much lower resolution than the validation model, hence each sample represents a much broader area.
This being said, it is worth noting that all the strategies generated models whose wetter areas and dryer areas correspond to those of the validation model. Also the soil moisture maps produced give a very good estimation of the areas were water deficit and concentration are in the field. Most likely, the miss-alignment between the validation model and the model outcome could have been overcome by having a longer mission. The fact that resulting model can discriminate wet and dry areas in such a short time (the validation model required more than 60 hours of work) is very encouraging.

Conclusion
In this paper, we proposed an exploration framework for autonomous mobile robots equipped with a soil moisture sensor to create high quality soil moisture maps. The sensor is a novel device based on fast neutron counting which enables non-contact measurements of soil moisture. Such a class of sensors can be modelled by the Poisson distribution and we demonstrated how to integrate such measurements into the kriging framework. We also investigated a range of different exploration strategies and assessed their usefulness in different scenarios. The proposed framework was evaluated on a range of datasets based on real soil moisture data collected from two different fields.
One of the important findings of the paper is the fact that the sampling regime's contribution to the overall exploration process is highly dependant on the characteristics of the field. In fields with high variability and less uniform distribution of soil moisture, the use of Adaptive Measurement Interval shows significant improvements in model quality compared to a Fixed Measurement Time regime. We also demonstrated that adaptive sampling strategies guarantee lower navigation times and spend more time obtaining samples leading to models of comparative quality to the non-adaptive strategies but with a much lower travel distance. This is especially important in large fields where travelling takes a significant proportion of the exploration time. Greedy methods tend to sample the outer border of the environments, which is where the kriging variance is usually higher. They tend to miss localised patches, although their overall model quality is comparable. For small fields with uniform soil moisture distributions these might be preferable exploration strategies.
Although the presented framework was demonstrated for the soil moisture mapping, it is a general approach which can be used to map other soil properties such as compaction, chemical composure, etc. It is a framework that would be particularly suitable in scenarios where the measured phenomena directly affects the acquisition time and needs to be spatially mapped. This includes applications such as rainfall measurements, people and animal counting, gas detection etc. One of the follow up questions arising from this research is if changing the time measurement regime on the fly could improve the resulting models even further. Future work could also address the additional path planning constraints caused by the layout of typical agricultural fields which feature soil beds and rows. Finally, the framework will be extended to map multiple soil properties at the same time.