An Extended Kriging method to interpolate soil moisture data measured by wireless sensor network

In recent years, wireless sensor network (WSN) has emerged as a new technique to collect Earth observation data at a relatively low cost and minimal labour over large areas. However, WSN observations are still point data. To determine the 10 spatial distribution of a land surface parameter, interpolation of these point data is necessary. Some geostatistical interpolation methods, such as the Ordinary Kriging method, Universal Kriging method, Co-KrigingCokriging method, and Regression Kriging method, have been used in various fields. However, capturing the spatial distribution pattern of heterogeneous land surface parameters is still difficult. For example, near-surface soil moisture is a critical parameter for agriculture management, and hydrological and ecological research. However, as soil moisture is related to many factors such as topography, soil 15 type, and vegetation, even a WSN observation grid is not sufficiently dense to reflect its spatial distribution pattern. This study developed a method to interpolate WSN-measured soil moisture with the aid of remote sensing images. The underlying idea is extension of to extend the traditional Kriging algorithm by introducing spectral variables, specifically, Normalized Difference Vegetation Index (NDVI) vegetation index (VI) and albedo, from satellite imagery as supplementary information to aid interpolation. Thus, the new Extended Kriging algorithm operates on spatial and spectral combined space. The algorithm 20 has been applied to WSN-measured data in the HiWATER campaign to generate daily soil moisture maps in the 4.5 km × 5 km oasis area in the middle reaches of the Heihe River, western China, from June 10 to July 15, 2012. Visual inspections indicate that the result from the Extended Kriging algorithm shows more spatial details than that of the traditional Kriging algorithm, and the temporal variation of patch-average soil moisture is, in general, consistent with precipitation/irrigation data. Leave-one-out cross-validation was also adopted to estimate the interpolation accuracy. The Root Mean Square Error (RMSE) 25 of the Extended Kriging method was found to be smaller than that of the original Ordinary Kriging method. Analysis with minimum variance of error (σκ σk), a self-uncertainty indicator given by the Kriging algorithm, also gave the same conclusion. Further testing indicated that if high-resolution land surface temperature maps are available, they can be added to the spectral variables and further improve the interpolation accuracy.


Introduction
Wireless sensor network (WSN) is a novel technique for ground data collection that is currently in high demand.It has been applied in various fields, such as hydrology, soil environment, atmospheric environment, forest meteorology, and fire disaster (Othman and Shazali, 2012).WSN is a cross-discipline technology that integrates sensor, automatic control, communication, and data analysis (Gong, 2007).In contrast to traditional ground data observation methods, the multiple sites comprising a WSN facilitate measurement of regional distribution of parameters.Further, its small size and relatively low price (Othman and Shazali, 2012) enables more nodes to be installed in the study area to meet research requirements.This helps us to obtain more extensive and much denser observation data of a specific land surface parameter.Compared with satellite remote sensing technique, WSN is a more flexible platform that can support a larger variety of sensors (Gong, 2007;Li and Gao, 2008).It can measure leaf area index (LAI), soil moisture, soil temperature, electrical conductivity, and atmospheric infrared radiance temperature at each observation node (Li et al., 2015;Dou et al., 2016).In contrast to remote sensing, in which data are obtained during overpass periods, monitoring by WSN is continuous and in real time.This complements the area data measured by remote sensing and also provides reference data for validation of satellite and aerial remote sensing (Li et al., 2016).
However, WSN only measures point data, whereas spatially continuous data over a certain region are increasingly required in environmental researches and management in order to make effective decisions or justified interpretations (Li and Heap, 2011).Consequently, to determine the spatial distribution of a land surface parameter, data interpolation is necessary.Some traditional interpolation methods, such as the inverse distance weighting method (ISW), polynomial method, and Kriging method, have been widely used to convert point data to spatial distribution.The ISW method and the polynomial method are types of non-geostatistical interpolation methods, which are relatively fast and easy to compute, but have limited estimation accuracy (Lu and Wong, 2008).By In contrast, the Kriging method uses variogram analysis to estimate the spatial variation structure, and takes spatial autocorrelation into consideration (Aalto et al., 2013).The Kriging method is a group of stochastic interpolation methods comprising the Simple Kriging method, Ordinary Kriging method, Universal Kriging method, Co-KrigingCokriging method, Regression Kriging method, and Residual Kriging method (Oliver and Webster, 1990).These methods have been used in various studies.For example, Maruyama et al. (2010) estimated peak ground velocities using the Simple Kriging method.Pokhrel et al. (2013) used the Simple Kriging method to estimate liquefaction potential and validated interpolation results.Wu et al. (2013) utilized the Residual Kriging method with the input variables of latitude, longitude, and elevation to estimate the average monthly temperature in the United States.Aznar et al. (2013) used the Co-KrigingCokriging method and the Canadian regional climate model as secondary information to spatially interpolate a mean monthly temperature time series.Liang et al. (2016) used the Co-KrigingCokriging method to estimate daily NO3-N loads in an agricultural river with the assistance of daily discharge.Through the development of geostatistical interpolation methods, more information about the rule or pattern of the distribution of the target parameter was explored, either through statistics of this parameter, or through its correlation with other parameters.
Soil moisture is a vital factor in agriculture, ecology, and hydrological studies (Dripps and Bradbury, 2007;Toth et al., 2008;Yu et al., 2011;Wang et al., 2012).Despite advances in recent interpolation algorithms, estimation of the spatial distribution of near-surface soil moisture is still unsatisfactory because soil moisture is highly heterogeneous on both spatial and temporal scales even over short distances (Gómez-Plaza et al., 2000).Many researchers have applied remote sensing to invert soil moisture.For example, Fan et al. (2015) improved the performance of Ts/VI space in retrieving soil moisture based on Compact Airborne Spectral Imager (CASI)/Thermal Airborne Spectrographic Imager (TASI) data.El Hajj et al. (2016) used neural networks (NNs) to estimate surface soil moisture from X-band SAR data over irrigated grassland areas.Ponnurangam et al. (2016) used a compact polarimetric decomposition with surface component inversion to estimate soil moisture on bare and vegetated agricultural soils.However, soil moisture imposes significant difficulty in quantitative remote sensing inversion because optical remote sensing is not directly sensitive to soil moisture, whereas thermal infrared and microwave remote sensing usually have low spatial resolution.Ground observation data, such as that from WSN measurement, are always necessary as supplementary data for remote sensing inversion of soil moisture.In other words, remote sensing information can be viewed as supplementary data to aid the interpolation of ground-measured soil moisture.In this light, Yao et al. (2013) used four spatial interpolation methods to estimate soil moisture in a complex terrain catchment and Regression Kriging proved to be optimal in producing a map with more details and better accuracy than other methods.Liao et al. (2016) analyzed analysed various sources of uncertainty, such as soil properties and terrain indices, while estimating near-surface soil moisture content with the aid of Co-KrigingCokriging at two typical hill slopes.However, as soil moisture is correlated affected by with many factors, such as elevation, vegetation, temperature, and irrigation, it is rather difficult to single out an optimally related one significantly linear correlated factor to aid interpolation in the Regression Kriging or Co-Kriging algorithm.Some researchers have also tried to use data assimilation methods in soil moisture estimation.Gao et al. (2014) estimated the spatial pattern of soil moisture using the BME method, which is based on WSN data and auxiliary information from ASTER (Terra) land surface temperature measurements.
In this paper, we propose extending the Kriging method by introducing high-resolution remote sensing imagery spectral variables into the interpolation algorithm.The spectral variables, such as vegetation index (VI), albedo, and temperature, can be either directly or indirectly relevant to soil moisture.The new Extended Kriging algorithm operates in the combined space of spectral dimension and spatial dimension.Thus, the semivariogram is also estimated in the combined space.The proposed algorithm has been tested with soil moisture data acquired by WSN in the HiWATER campaign to generate time series of a daily soil moisture maps at 30 m spatial resolution and daily temporal resolution.

Study area
The study area, shown in Fig. 1, is located in Zhangye oasis in the middle reaches of the Heihe River Basin (HRB) in northwestern China (38.87°N,100.40°E).The HRB is the second largest inland river basin and is characterized by large areas of alpine cold and arid landscapes and a small portion of oasis agricultural land.The main oasis in the middle reaches is agricultural land with a variety of vegetation, including trees, maize, wheat, and vegetables.As the potential evaporation is very high, ranging from 1200 mm to 1800 mm per year, while the average annual precipitation in the artificial oasis is only 177 mm, irrigation is the primary source of water for crops (Li et al., 2013).
The HRB has long served as a test bed for integrated watershed studies and hydrological experiments (Cheng, 2009).
HiWATER is an ongoing watershed eco-hydrology comprehensive experiment that began in 2012 (Li et al., 2013).With the objective of improving the comprehensive observation ability, an eco-hydrology WSN was installed as a part of the HiWATER basic experiment.From May 2012 to September 2012, 50 WATERNET nodes were installed in a 5.5 km × 5.5 km forcifoci experimental area in the main oasis in the middle reaches.
In this study, we chose an area approximately 4.5 km × 5.0 km covered by WATERNET as the experimental area, which consists of 48 WATERNET nodes (Fig. 1).All the nodes were located in the cornfield.

Data resource
WSN data and remote sensing data were both used in this study.We used the data covering the period from June 10, 2012 to July 15, 2012, when the vegetation cover of the study area changed from sparse to dense.WSN data were obtained by the 48 WATERNET soil moisture sensors in the experimental area.Each of the WATERNET nodes measured soil moisture (SM) and soil temperature in two layers (4 cm, 10 cm), and the data were collected every 5 minutes (Jin et al., 2014).The soil moisture measurements are based on the frequency-domain reflectometry method using a Hydro Probe II (HP-II) sensor (Gao et al., 2014).Details of WATERNET design and other information are given in Jin et al. (2012).Because the data collected by some of the WSN nodes are affected by sensor noise and other abnormal conditions during wireless data transfer, smoothness and noise reduction treatments are necessary.Besides, when the soil moisture is close to saturation, the soil moisture sensors cannot work properly and sometimes give abnormal values (Zhang et al., 2015a).Thus, this part of data also need to be removed.
We first excluded the abnormal data by assigning the zero value, negative value and abnormally high value (soil moisture content > 50 percent) as invalid data NaN.Then, we averaged the data for the whole day, and used the average result as the final soil moisture value for each node.
Remote sensing data served as auxiliary data in estimation of the distribution of near-surface soil moisture.We used both satellite remote sensing images and air borne remote sensing images to derive spectral variables.NDVI and albedo were derived from CCD camera on the Chinese HJ satellite, which has four channels in the visible and near infrared spectral range.
The spatial resolution was 30 m and revisiting frequency approximately was 2 days.Owing to the influence of clouds, only five clear-sky images on the following dates during this period could be used: June 15, June 19, June 29, July 8, and July 13 in the year 2012.All the images used here were pre-processed with calibration, geometric rectification, and atmospheric correction.
Land surface temperature (LST) is also an important commonly used variable in monitoring soil moisture with remote sensing.The data were acquired from the airborne sensors of CASI/TASI on July 10, 2012, between 12:00 and 12:30 (local time), at an altitude of 2500 m (Xiao and Wen, 2013;Fan et al., 2015).The spatial resolution was 2.5 m.CASI acquires data covering the visible and near-infrared (VNIR) region of the spectrum.LST was derived from TASI, which is a hyperspectral thermal infrared sensor released by ITRES, Canada in 2006 (Wang et al., 2011).In this study, we introduced a soil moisture map as simulation data in calculating the semivariogram and some interpolation analysis.This soil moisture map, which is presented in Fig. 7, was also derived from CASI/TASI data (Fan et al., 2015).

Method
The new spatial interpolation method proposed in this paper is based on the traditional Kriging algorithm.The proposed method extends the traditional X and Y spatial coordinates to spatial and spectral combined coordinates, and utilizes remote sensing derived spectral variable NDVI and albedo in the interpolation algorithm as supplementary information.In this section, first, an outline of the traditional Kriging algorithm is given, then the technique employed to extend the Kriging algorithm is explained.While fitting the semivariogram of the soil moisture, as the number of WSN nodes is insufficient to gather robust statistics, we used the remote sensing-derived soil moisture map mentioned in Sect.2.2 to calculate the variance function.

Basic formula
Kriging is an interpolation method derived from regionalized variable theory, which inherited the concept from geostatistics (Oliver and Webster, 1990).It has been used to provide linear unbiased predictions at unsampled locations and depends on expression of the spatial variation of the variable in terms of the semivariogram (Burgess and Webster, 1980;Cressie, 1990).
This method quantifies and reduces the uncertainties of estimation, minimizing self-estimated prediction errors (Gao et al., 2014).The core of Kriging is an optimally linear unbiased estimator that can be expressed as follows (Journel and Huijbregts, 1978): where  * is the estimated value of the variable at location  0 ,  is the number of the closest neighboring sampled data points used for interpolation,   is the Kriging weight assigned to each observation (  ).
Optimal estimation requires the minimum variance of errors: To ensure unbiased estimation, the following constraint must satisfy the equation as follows: To solve this constrained optimization problem, the Lagrange Multiplier Method (LMM) is adopted.With Eq. ( 2) as the objective function and Eq. ( 3) as the constraint, the LMM minimizes the following cost function: where  is the Lagrange Multiplier.At the minimum point of the cost function, the differentiation of  with respect to each of its variables is zero.Thus, the optimization problem decomposes into one of solving the following set of equations: Differentiating the cost function, we have The minimum variance of error (  2 ), as is shown in Eq. ( 2), can be used as a quality indicator in estimation (Yamamoto, 2000).It can evaluate the intrinsic estimation uncertainty from the algorithm itself.

Estimating semivariance and semivariogram
Semivariance and semivariogram, containing spatial correlation information, are important concepts in geostatistics.The semivariance of variables at certain locations is estimated from the semivariogram function, which is a function of the distance between the two locations.Usually, a de-trending pre-process is applied to the observation data.After this pre-processing, the spatial distribution of the variable is assumed stationary, which means that the semivariance does not change with location.
On the basis of this assumption, the semivariance can be estimated from the data that a random variable is well correlated in space as a function of separation distance.The semivariance (γ) of Z between two data points is defined as where h is the distance between points   and  0 , and (ℎ) is the semivariogram (Webster and Oliver, 2001).
The semivariogram is usually estimated from the statistics of sample points as follows: where  is the number of pairs of sample points separated by distance ℎ (Burrought and McDonnel, 1998).
As the number of WSN nodes is insufficient to gather robust statistics, the soil moisture map retireved from airborne hyperspectral remote sensing was used here to derive the semivariogram function.This soil moisture map was derived from the airborne hyperspectral datasets of CASI/TASI (Fan et al., 2015), acquired on July 10, 2012, between 12:00 and 12:30 (local time), at an altitude of 2500 m.In the calculation, we sampled 3000 9000 random points every time from the soil moisture map to calculate the semivariance by repeating the sampling thrice.The semivariance calculated from the soil moisture map, which is shown in Fig. 2. Here, we assumed that this semivariogram could be applied to interpolate WSN measured soil moisture in the period from June 10, 2012 to July 15, 2012.
A spherical model was used in this study as the semivariogram model, which is extremely important for structural analysis and spatial interpolation (Burrought and McDonnel, 1998).
where (ℎ) is the semivariance;  0 represents a nugget, which is the minimum variability observed or the "noise" at a distance of zero;  1 is the structural variance,  0 +  1 represents the sill variance; and  is the range that signifies the correlation length in geostatistics.
The Kriging method requires the second-order stationarity for geostatistical inference and assumes it to be isotropic.As our study area was in the central part of the oasis, and no significant soil wetness spatial trend could be found, the soil moisture observations in the study area were assumed to meet the above requirements.Because the area of the oasis is limited, i.e., the closest desert is about 4 km away from the center of the study area, the semivariogram statistics beyond 4 km may be affected by the presence of desert.Therefore, the range of the spherical model was set as 2500 m. Figure 2 shows the fitting curve of the semivariogram obtained using the spherical model.To reduce error, the spatial distance was divided into ten bins.We averaged the semivariance values of each bin as the final semivariance value of each distance.As the amount of the sampling data was insufficient when the spatial distance was increased up to a certain extent, its average semivariance value was not stable.Therefore, these invalid data were removed before semivariogram fitting.Hence, we only used the data when h<=4000 m in the fitting process of Fig. 2.

Selection of spectral variables
In this study, we selected NDVI and albedo as the auxiliary information to aid the interpolation.It is based on the following two main considerations: (1) NDVI and albedo are the spectral indexes which are fairly easy to be obtained from almost all high resolution remote sensing data sources; (2) NDVI and albedo represent most of the remote sensing information in the visible and near infrared spectral range.Although selecting NDVI and albedo as spectral variables is out of practical considerations, it is still necessary to analyze the correlation between soil moisture and NDVI/albedo, which are shown in Fig.

3.
In Fig. 3 (a) and (b), the NDVI and albedo data were collected from five available HJ satellites images to compare with the soil moisture data observed by WSN of the corresponding dates; in Fig. 3 (c) and (d), the NDVI and albedo data were on July 13, comparing with the whole soil moisture inversion map after being masked.As can be seen from the comparison results, the NDVI and albedo are correlated to soil moisture to a certain extent, but the correlation coefficient is not significant (absolute value less than 0.5).An explanation of this result may be that NDVI can reflect the vegetation status, and vegetation usually grows good when soil moisture is abundant.However, in irrigated land, all the fields get enough irrigation; then the correlation between soil moisture and NDVI are weakened.In sparsely vegetated land, the albedo of dry soil is usually higher than that of wet soil.So, there is a correlation between soil moisture and albedo.But this correlation can be disturbed by soil type and the presence of vegetation.Actually, it is usually considered impossible to estimate soil moisture from visible and near infrared remote sensing data.
As land surface temperature map was acquired on July 10, 2012, with airborne remote sensing technique, we also analysed the correlation between soil moisture and LST.In order to get enough sample points, we use the soil moisture map from remote sensing inversion to correlate with LST, and the result is presented in Fig. 4.This analysis shows that LST has a better correlation with soil moisture (R=-0.53772)than that of NDVI and albedo, which explains why thermal infrared remote sensing is emphasised in drought monitoring.However, the correlation is not stable because of the complexity in land surface processes.It also should be mentioned that high resolution thermal infrared remote sensing image are currently rarely available data source, so, its potential can only benefit practical applications in the future.

Extending the Kriging method to incorporate remote sensing information
To reflect more details of the spatial distribution pattern of soil moisture, we propose a new algorithm that incorporates remote sensing variables, i.e., NDVI and albedo, into the basic Kriging method.The traditional interpolation space is the spatial space depicted by x and y coordinates.The new algorithm extends the interpolation space to the combined spatial and spectral space, in which NDVI and albedo are treated as coordinates, just like  and y.The distance in the combined space is characterized by the spatial distance and the spectral distance, as follows: where ℎ is the spatial distance, Δ and Δ are the coordinate differences between two sampled points,  represents the spectral distance, ΔNDVI and Δalbedo are the differences of NDVI and albedo values between two sampled points, and σ DNVI and σ albedo are two normalization factors (in this study, we simply set their values as 0.1, 0.1).
If LST is also used together with NDVI and albedo, the Eq. ( 11) will be expended as Eq. ( 12).
where, ΔLST represents the difference value of temperature between two sampled points from the airborne LST map, which is aggregated to 30 m resolution to match with NDVI and albedo, and σ LST is the normalization factors, the value of which is one degree.
Correspondingly, the semivariogram model was extended to the spatial and spectral combined space, as signified in the equations below: where  1 and  2 are the semivariogram values with respect to ℎ and ,  is the overall semivariogram, and  1 and   are the lag distances of the spatial and spectral variables.
We also used the soil moisture map to derive the semivariance statistics.To reduce the error, the spatial/spectral distance was divided into ten bins as in that of traditional Kriging.The bin values with insufficient sample numbers were removed before semivariogram fitting.We only used the data when h<=4000 m and s<=4 in the fitting process.The semivariance, as a function of ℎ and , is shown in Fig. 5.
Here, the X-axis is spatial distance, the Y-axis is spectral distance, and the color in each grid represents the average semivariance value of the soil moisture.
When h > 4000 and s > 4, the sampled data quantity is not sufficient to satisfy statistical significance.Therefore, we divided ℎ into six intervals ranging from 0 to 4000 m, and s into six intervals ranging from zero to four; the semivariance in Fig. 4 is the average value in these intervals.
Using the semivariagram model as in the above Eq.( 13), Eq. ( 14), and Eq. ( 15), the fitting semivariance diagram can be obtained, as shown in Fig. 56.The  1 of spatial distance was pre-set as 2500 m, and the  2 of spectral distance was preset as 3.5.Then, the fitted parameters  0 ,  1 , and  2 were 6.4468, 2.8762, and 2.9972, respectively.

Analysis based on simulation data
In this study, we used two or more auxiliary variables to aid the interpolation.Multivariate geostatistical technique algorithms, such as Cokriging, also provide various methods for combining auxiliary information.Here, we compared the proposed Extended Kriging with Cokring, using NDVI, albedo and land surface temperature as covariates.The Cokriging interpolation was implemented with the GSTAT package in R programming language.
In order to conduct this analysis and also ensure the amount of sampling points, we used the simulation data, the soil moisture map derived from airborne hyperspectral remote sensing, to extract sampling points and validation points.Firstly, as in Sect.3.1.2,we sampled 9000 random points to calculate the semivariogram of Cokriging method.Then we used different numbers of points, ranging from 300 points to 30 points, to interpolate soil moisture.The numbers of sampling points are shown in Table 1.All the points in the soil moisture inversion map were used as validation points to calculate RMSE (Root Mean Square Error) and average   , which is introduced in Sect.3.1.1as self assessment of interpolation uncertainty in Kriging algorithm, for each interpolation.As the locations of sampling points may influence the interpolation result and accuracy, we sampled randomly and repeated enough times to decrease the disturbing of point locations, and calculated the average RMSE value.For example, when the number of interpolation points is 30, we randomly sample 30 points for 100 times, interpolate respectively, and calculate the mean value of RMSE.The times of repetition and the interpolation accuracy comparison results are shown in Table 1.
From Table 1, we can see that the interpolation uncertainty (indicated both by σ κ and RMSE) decreases while the number of sampling points increases, and the estimator σ κ can reflect the variation trend of the actual RMSE.We also find out that the RMSE of Extended Kriging and Cokriging are very close, except that Cokriging performs a little better when the number of sample points is less than 50.The performance of Extended Kriging is similar with that of Cokriging in the quantitative perspective.However, the advantage of the Extended Kriging can be perceived when the interpolation results of the two methods are visually compared with the reference soil moisture map (Fig. 7). Figure 8

Results and discussion
The proposed Extended Kriging method was applied to the WATERNET measurements in the HiWATER campaign to generate daily soil moisture maps in the 4.5 km × 5 km oasis area from June 10 to July 15, 2012.To evaluate the proposed performance of Extended Kriging method, we first visually inspected the interpolation results and compared them with the results obtained using the tradition Kriging method.Then, we examined the precipitation and irrigation data to ascertain whether the temporal variation of the interpolated soil moisture was consistent with the water input.Subsequently, a quantitative uncertainty analysis was conducted with both Kriging's internal uncertainty estimator   , and the RMSE from leave-one-out cross-validation method.Finally, we demonstrated that if a high-resolution remote sensing image of surface temperature is available, the interpolation accuracy could be significantly improved with this additional spectral information.

Spatio-Spatial and temporal trend analysis
A spatial map of the interpolated soil moisture can directly reflect the capability of the interpolation methods to obtain the parameter's variation pattern.Compared with the original Kriging method, the advantages of the new Extended Kriging method can be visually observed in the Spatio-spatial and temporal pattern of the interpolation results.We used the semivariogram obtained from Sect.3.3, soil moisture data (38 nodes) of the WATERNET at a depth of 10 cm on July 10, 2012, and the remote sensing data (NDVI and albedo) on July 13, 2012, to estimate the soil moisture distribution of the experimental area.The interpolation results obtained for these two methods are shown in Fig. 9, where (a) is the result of the original Kriging method and (b) is the result of the Extended Kriging method.Both result maps are masked to exclude residential areas.
It is clear that the result for the Extended Kriging method shows more details of the field soil moisture distribution than that for the original Kriging method.The interpolation result for the original Kriging method can only show a smooth and continuously changing trend surface (Fig. 9a).In contrast, from the Extended Kriging interpolation result, it is clear that the soil moisture is not smooth and continuously changing.Figure 9b reveals a few rough veins in the ground surface and subtle changes in the field soil moisture.
In addition to measuring data from multiple sites simultaneously, WSN is also able to obtain real-time and continuous observations.This enabled us to estimate the continuous soil moisture distribution of the study area on both spatial and temporal domains.We interpolated the soil moisture distribution for every day from June 10, 2012 to July 15, 2012, and the average number of WSN nodes could be used in this period is about 39, ranging from 36 to 43.Consequently, we obtained 36 interpolation maps.Six of the maps, with seven-day intervals, are shown in Fig. 10.As can be seen, the spatial distribution of the soil moisture content changed significantly over time, indicating the necessity for dynamic monitoring of soil moisture.
The drought areas and their borders can be clearly identified in the interpolated maps, which is very useful information for agricultural management.
variations to soil moisture.However, there is a second explanation: because the remote sensing image was not available around the precipitation dates, we had to use an image that was acquired on an earlier date.For example, the interpolation result for June 28 is based on the image acquired on June 19; thus, this large uncertainty resulted from temporal mismatch between remote sensing data and ground data.

Effect of adding temperature information to the interpolation method
Land surface temperature (LST) is an important index connected with soil moisture, usually appearing as negatively correlated with soil moisture.Therefore, it is a potential spectral index that can support interpolation of soil moisture.However, high-resolution temperature remote sensing data are not as widely available as those of NDVI and albedo are.Fortunately, an LST map at 2.5 m resolution is available for July 10, 2012, as hyperspectral thermal infrared airborne images were acquired by the TASI instrument on that date.Therefore, we added the LST as the third spectral index, together with NDVI and albedo, in the interpolation of this date, and compared the result with that obtained prior.Equation ( 12) is revised as Eq. ( 16).
Estimation maps are shown in Fig. 1215, where (a) is the result with the indexes of NDVI and albedo and (b) is the result with the indexes of NDVI, albedo, and LST.The comparison of uncertainty is shown in Table 2.
Here, ΔLST represents the difference value of temperature between two sampled points from the LST map, which is aggregated to 30 m resolution, and σ LST is the normalization factors, the value of which is one.
The result indicates that introducing a new spectral temperature index can further improve the accuracy of soil moisture content: the value of   is reduced from 9.20703.0343to 8.00432.8292and the value of RMSE from 1.6406 to 1.3958.Hence, if high-resolution LST data are available in a long time series, the future of soil moisture interpolation task can incorporate the LST information to boost accuracy.

Conclusion
With the rapid development of ground-based Earth observing techniques such as wireless sensor network, we are now able to monitor environmental parameters in real time, continuously, and with multiple sample points.However, interpolation is still needed to extend the point measurement to spatial distribution of the corresponding parameter in an area.As satellite remote sensing is an efficient way of acquiring area Earth observing data, it is desirable to combine information from remote sensing and from ground-based observation networks.
The Extended Kriging method proposed in this study introduces the remote sensing image spectral information into the traditional interpolation method.NDVI and albedo are the spectral variables used in the algorithm.These spectral variables are treated in the same manner as the spatial variables, i.e.,  and .Therefore, the interpolation is fundamentally the same Kriging algorithm, but operating on the combined space of spatial dimension and spectral dimension.The semivariogram model is also extended to the combined space.A remote sensing derived soil moisture map is used in this paper to fit the semivariogram model.However, this soil moisture map can be replaced by other sources of samples as long as the dataset is sufficiently large to derive robust statistics about the semivariance.
The proposed algorithm was applied to the soil moisture dataset acquire by the soil moisture sensors network (WATERNET) in the oasis agricultural areas, which is the forcifoci experimental area of the HiWATER campaign.As the WATERNET provides continuous near-surface soil moisture measurement over 48 scattered points, the interpolation results are daily soil moisture maps from June 10, 2012 to July 15, 2012, covering an area approximately 4.5 km × 5.0 km in size.
Visual inspections indicate that the interpolation result from the proposed Extended Kriging algorithm presents much more spatial details than that of the traditional Kriging algorithm.The field-average soil moisture of several irrigation fields for long time series are associated with the precipitation data and irrigation data, and the temporal variation of soil moisture can be well explained by these water inputs.The quantitative uncertainty analysis with both the leave-one-out method and   indicate that the Extended Kriging algorithm, which operates in the spectral and spatial combined space, produces more accurate interpolation results than that of the traditional Kriging algorithm, which operates only in the spatial domain.Currently, NDVI and albedo are recommended as the spectral variables to aid interpolation because they can be easily derived from most highresolution satellite images.However, we demonstrated in the discussion that more relevant spectral variables, such as land surface temperature, could be incorporated into this algorithm to improve its performance.However, how to choose the informative spectral variables remains an open topic for this algorithm.
There are other methods that can combine information from ground measurement and information from remote sensing, e.g., the Co-KrigingCokriging, the Kriging with External Drift, and the data assimilation method (Gao et al., 2014).Although some results from the Cokriging are presented in Sect.4, we prefer not to compare the proposed algorithm with other sophisticated algorithms in terms of accuracy for the following two considerations: in the first place, the Extended Kriging algorithm is much simpler than the Co-KrigingCokriging, the Kriging with External Drift, and the data assimilation method.
Thus, it is possibly applicable in situations where the pre-conditions of other sophisticated algorithms are not satisfied.Second, as a new algorithm, the Extended Kriging algorithm utilizes easily available variables such as NDVI and albedo, which are difficult to use in the Co-Kriging method and the data assimilation method because NDVI and albedo do not have a clear physical connection with soil moisture.The Extended Kriging algorithm has potential, as well as space for still needs improvements.For example, the normalization factors for spectral variables could be refined; and the semivariogram model in the combined space is too simple.Nevertheless, we prefer to present the simplest form of the algorithm to the reader, and leave these improvements to future researchers.
There are other aspects that we could not discuss deeper in this short paper.One of them is the quality, or accuracy, of WSN-measured soil moisture.As we know, it is technically difficult to install Hydro Probe II (HP-II) sensors exactly at 4 cm and 10 cm below the surface; and the WSN sensors are more or less infected by its internal voltage and temperature; and the sensors can go to saturation in extremely low or high soil moisture values (Gao et al., 2014;Zhang et al., 2015b).It is also true that the remote sensing data may suffer from inaccurate calibration and atmospheric correction.Besides the flaws of data, the scale mismatch between the footprint of WSN nodes and 30 m resolution remote sensing pixels should be considered.
Nevertheless, currently we do not have solid data to support a study on these aspects.
Another related interesting topic is the desired density of the WSN nodes, and how to optimize the location of the nodes.
Fortunately, a similar topic has been addressed by other researchers (Wu et al., 2016), although their research area and target parameter are different.In this paper, we simply used all the valid WATERNET data.
The proposed algorithm is a development of the classical Kriging method.Although it is proposed in this paper to interpolate the soil moisture data, it is potentially applicable to other environmental parameters.As new observation technologies are being applied wider, more and more high-quality measurement data, at multiple sites in a small area, will become available to the public, and the potential of the Extended Kriging can be further explored.
and E[(  )(  )], then the equations can be solved.These values are estimated by the semivariogram function in Sect.3.1.2.
shows the interpolated soil moisture map by Extended Kriging and Cokring, both with the aid of NDVI and albedo: (a) is the result of Extended Kriging interpolated with 300 points; (b) is the result of Extended Kriging interpolated with 30 points; (c) is the result of Cokriging interpolated with 300 points; (d) is the result of Cokriging interpolated with 30 points.We can find that the interpolation results of Extended Kriging presented more detailed information for the spatial distribution of soil moisture than the results of Cokriging.

Figure 2 Figure 4 Figure 5 Figure 6 Figure 9 5Figure 12 Figure 14
Figure 2 Semivariance of sampled soil moisture data and fitting curve of the soil moisture semivariogram: the semivariance was calculated by the 9000 random sampling points from the soil moisture map; the semivariogram was fitted by spherical model with the data of h<=4000 m, and the range was set as 2500 m

Table : 5
Table1Comparison of the interpolation results for the Extended Kriging method and Cokriging method, using simulated data and the aid

Table 2
Comparison of the interpolation results for the traditional Kriging method and the Extended Kriging method on July 10, 2012