Mapping High-Resolution Soil Moisture over Heterogeneous Cropland Using Multi-Resource Remote Sensing and Ground Observations

High spatial resolution soil moisture (SM) data are crucial in agricultural applications, river-basin management, and understanding hydrological processes. Merging multi-resource observations is one of the ways to improve the accuracy of high spatial resolution SM data in the heterogeneous cropland. In this paper, the Bayesian Maximum Entropy (BME) methodology is implemented to merge the following four types of observed data to obtain the spatial distribution of SM at 100 m scale: soil moisture observed by wireless sensor network (WSN), Advanced Spaceborne Thermal Emission and Reflection OPEN ACCESS Remote Sens. 2015, 7 13274 Radiometer (ASTER)-derived soil evaporative efficiency (SEE), irrigation statistics, and Polarimetric L-band Multi-beam Radiometer (PLMR)-derived SM products (~700 m). From the poor BME predictions obtained by merging only WSN and SEE data, we observed that the SM heterogeneity caused by irrigation and the attenuating sensitivity of the SEE data to SM caused by the canopies result in BME prediction errors. By adding irrigation statistics to the merged datasets, the overall RMSD of the BME predictions during the low-vegetated periods can be successively reduced from 0.052 m·m to 0.033 m·m. The coefficient of determination (R) and slope between the predicted and in situ measured SM data increased from 0.32 to 0.64 and from 0.38 to 0.82, respectively, but large estimation errors occurred during the moderately vegetated periods (RMSD = 0.041 m·m, R = 0.43 and the slope = 0.41). Further adding the downscaled SM information from PLMR SM products to the merged datasets, the predictions were satisfactorily accurate with an RMSD of 0.034 m·m, R of 0.4 and a slope of 0.69 during moderately vegetated periods. Overall, the results demonstrated that merging multi-resource observations into SM estimations can yield improved accuracy in heterogeneous cropland.


Introduction
Soil moisture (SM) plays a fundamental role in land-atmosphere exchange processes [1], hydrological processes [2] and terrestrial water cycle trends [3].L-band microwave remote sensing is the most promising technique for obtaining global SM measurements with frequent revisit times, and is independent of the effects of clouds and solar illumination [4].Satellite-based microwave instruments, such as the Soil Moisture and Ocean Salinity (SMOS) satellite [5], and Aquarius [6] have demonstrated this capability.However, the spatial resolution of SM products from these instruments is lower than 40 km [7].Although it is possible to use data with such a coarse spatial resolution in the field of hydrology, most hydrological processes are better observed and modeled at resolutions of less than 1 km [8].The fine resolution requirement is also derived from agriculture-related applications and river-basin management [9].For instance, irrigation events affect the spatial distribution of SM at the field scale (100-1000 m) in croplands [10].Active microwave sensors (synthetic aperture radar (SAR)) have fine spatial resolution, but SAR (e.g., ENVISAT) backscatter data in high frequency is not sensitive to soil moisture in irrigated and vegetated areas [11][12][13][14].
With the development of the wireless sensor networks (WSNs) which make it possible to acquire simultaneous measurements of regional SM, strategies for merging WSNs and remote sensing data are currently being investigated to refine high resolution observations using a variety of optical sensors in the context of geostatistics [15,16].Visible/infrared/thermal sensors are applied to provide indirect measurements of SM, and have resolutions that are 10-500 times finer than the resolutions of microwave radiometers.WSN/optical data merging techniques for obtaining fine resolution SM information typically rely on the intrinsic relationships between vegetation indices (VIs) and land surface temperature (LST) with SM using empirical approaches [17].Based on these relationships, a range of optical SM indicators has successfully been merged with WSNs to obtain SM pattern, such as LST [18], temperature vegetation dryness index (TVDI) [19] and apparent thermal inertia (ATI) [20,21].
However, the previous studies have shown that the validity of WSN/optical data merging methods is often questionable in heterogeneous land surfaces [18,19].The non-stationary and anisotropy of SM spatial correlation structure result in problems when applying WSN measurements directly for geostatistics.In order to account for this issue, some anisotropic semi-variograms need to be applied using remote sensing variables.Although the use of some optical indicators (e.g., TVDI and LST) [19] could somewhat mitigate the anisotropy of the measured SM data, the residuals still would not be isotropic or meet the requirements of geostatistics because some optical SM indicators exhibit low sensitivity to the extreme values of SM [18,19,22].This low sensitivity increases the errors in spatial structures which are key to geostatistical methods.Moreover, optical SM indices have low sensitivity to SM under moderate to dense vegetation covers [23].Dense canopies result in saturated VIs and attenuate LST variations related to SM [10].An effective method of merging multi-sensor SM should enable the integration of sufficient information from multi-resource observations whilst accounting for issues of SM anisotropy and dense vegetation cover.As a nonlinear spatiotemporal mapping method, Bayesian Maximum Entropy (BME) is capable of integrating data with various accuracies and from various sources [24][25][26].In addition, BME is a probabilistic method that can enable the application of uncertain data to provide information of interest while accounting for the data uncertainties.BME has been successfully used in many research fields, including environment assessment [27][28][29][30], soil science [31][32][33], and public health [34,35].Furthermore, BME has shown potential as a merging method for integrating WSN measurements and remote sensing data [36].
In this paper, we attempt to present a multi-source fusion method to estimate high resolution SM data over heterogeneous cropland by merging the WSN observations, irrigation statistics, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Polarimetric L-band Multi-beam Radiometer (PLMR)-derived SM products in the BME paradigm.This paper is organized as follows: Section 2 presents the study area and data used.The methodology, including BME principles and the merging framework, is introduced in Section 3. Section 4 presents the results and evaluates the BME predictions obtained by merging multi-resource observations from the perspective of accuracy.Finally, the entire paper is discussed and conclusions are presented in Sections 5 and 6, respectively.

Data
Heihe Watershed Allied Telemetry Experimental Research (HiWATER), was conducted from 2012 to 2015, located in the Zhangye artificial oasis in the middle reaches of the Heihe River Basin in northwestern China.The implementation of HiWATER is to reveal the processes and mechanisms of an eco-hydrological system in an inland river basin [37].While a full description of the data set is given in [2], a brief overview of the most pertinent details is provided here.The data used in this study include WSN measurements, irrigation statistics, ASTER data, and PLMR-derived SM products collected throughout the experimental area during periods of low (30 May 2012; 24 June 2012) and moderate vegetation (10 July 2012).

Soil Moisture Wireless Sensor Network and Irrigation Statistics
As one part of HiWATER, an ecohydrological WSN was installed in the experimental area covering 4.2 × 4.2 km 2 as shown in Figure 1.The objective of WSNs is to capture the multi-scale spatiotemporal SM dynamics, and validate the remote sensing SM products [38].The experimental area is covered predominantly by crops and rarely artificial facility (e.g., buildings, roads).The WSN consists of 48 WATERNET nodes for collecting SM data at depths of 4 and 10 cm based on the frequency-domain reflectometry method using Hydro Probe II (HP-II) sensors.The distribution of WSN nodes was designed using an optimal spatial sampling strategy [39].Table 1 shows a summary of WSN data used in this research.WSN measurements were obtained at a high temporal resolution (in 10-minute intervals), thereby ensuring their concurrency with the remote sensing observations.

PLMR-Derived Soil Moisture Products
PLMR, an airborne microwave radiometer, is capable of dual-polarized brightness temperature measurements (TB) at a frequency of 1.413 GHz and bandwidth of 24 MHz.The instrument is equipped with six receiving beams, with fixed look angles of ±7°, ±21.5°, and ±38.5° [40,41].TB are used to derive surface SM in this study.TB -to-SM inversion was performed using a modeling approach based on a τ − ω model combined with surface roughness, single scattering albedo, canopy, and polarization mixing parameters [42].The PLMR SM product for 10 July 2012, is available with a spatial resolution of 700 m.The PLMR SM product was well correlated with the ground-observed SM data (RMSE = 0.04 m 3 • m −3 , where RMSE is the root mean square error).A full description of the retrieval and validation of PLMR SM products is presented in [43].

ASTER Data
ASTER is an advanced multispectral imager which is to fly on EOS-AM1 polar orbiting spacecraft [44].The ASTER records in 15 spectral bands covering from visible and near-infrared, shortwave-infrared, and thermal infrared ranges, with spatial resolutions of 15, 30, and 90 m, respectively [44].Three clear-sky ASTER images were acquired over the experimental region on 30 May, 24 June and 10 July 2012.The normalized difference vegetation index (NDVI) was computed using the ASTER red and near-infrared bands with a spatial resolution of 15 m, then re-sampled at a resolution of 100 m.This procedure can reduce the spatial scale effects of the NDVI [45].The LST data were derived from the ASTER data at 90 m spatial resolution using the temperature emissivity separation algorithm [46], combined with the water vapor scaling atmospheric correction method [47], then LST data were re-sampled at 100 m resolution.The RMSE of LST for 30 May, 24 June, and 10 July are 1.56 k, 1.66 k, and 1.77 k, respectively.

Methodology
BME is a spatiotemporal analysis and mapping method based on a combination of Bayesian law and maximum entropy theory [26].This combination permits logical incorporation of measurements (hard data and soft data) from various knowledge bases [25].In the BME framework, the measured values are defined as hard data, and specific observations are treated as soft data which are uncertain measurements expressed as interval values, probability statements, etc. [24].In this study, the BME method was used to estimate the distribution of SM at a spatial resolution of 100 m by merging WSN measurements, ASTER-derived soil evaporative efficiency (SEE) data, irrigation statistics, and PLMR-derived SM.The general scheme of this study is depicted in Figure 2. First, only WSN data was applied to predict SM at 100 m-scale, and the predictions were denoted by SM from Method I. Second, ASTER-derived SEE and WSN data were merged to predict SM at 100 m-scale, and the predictions were denoted by SM from Method II.Third, irrigation statistics were added to the merged datasets to estimate SM, and the estimation was denoted by SM from Method III.Fourth, the downscaled SM at 100 m-scale obtained from PLMR SM was further added to the merged datasets, and the predictions were denoted by SM from Method IV.

Predictions Obtained by WSN (Method I)
For comparison of BME predictions using multi-resource observations, only WSN data was used to estimate SM firstly.In the framework of BME, when only hard data are available, BME algorithm reduces to the Kriging methods [48].In this study, ordinary kriging (OK) method was applied to estimate SM using only WSN data.
OK, as an optimally-linear unbiased method, can be expressed: where  * (  ) is the optimal estimator at location s o , (  ) is the WSN measurement at the location   ,  is the number of WSN measurements.  (  ) is kriging weights at location   .
To ensure unbiased estimation, the variance between  * (  ) and (  ) is required to be minimum: where σ must follow the constraint condition Equation (3) to ensure linear unbiased estimation.
where the weight coefficients   and the optimal estimator  * (  ) are calculated using the Lagrange Multiplier Method with Equation (2) as the object function and Equation (3) as the constraint.  can be calculated by the covariances [49]: where  is the matrix of the covariances (ℎ) expressed as a function of the separated distance ℎ.   is the vector of the covariances at unknown location.(ℎ) can be calculated by the variogram(ℎ) : where   represents the nugget, which characterizes the minimum variability observed or the "noise" at ℎ = 0, and  − is the partial sill representing the structural variance.  and  − are estimated by Equation( 6)- (7).(ℎ) can be calculated using the semi-variance function: where ℎ is the separated distance between the random variables   and  +ℎ .(ℎ) is the number of pairs of points of   and  +ℎ at a given separated distance ℎ .Generally, (ℎ) generally can be modeled by variogram models such as an exponential, spherical or Gaussian models.Gao, et al. found that the variogram of WSN observations can be fitted by the nugget-exponential model better [18].Thus, the nugget-exponential model was selected to describe (ℎ), which is expressed as [50]: where  is the lagged distance representing the correlation length.

SEE Calculation
SEE represents the surface moisture availability [51], which is the ratio of the actual to potential soil evaporation.The rationale for choosing SEE as soft information in the BME procedure is based on its good correlation with near-surface SM and relative stability during the daytime on clear-sky days [52].SEE can be derived from ASTER LST and NDVI using Equation (8) [53]: where   is the soil skin temperature;   is the maximum   , representing the soil temperature at the minimum SM; and   is the minimum   , representing the soil temperature at the maximum SM.Using the triangle approach,   is estimated as: where  , is the ASTER LST,   is the vegetation skin temperature, and   is the fraction of vegetation cover.Herein,   is defined as the temperature of the bare soil.  is assumed to be uniform in the experimental area, equal to the minimum ASTER LST [53,54].The effect of root-zone SM on   was not considered in the calculation of SEE.  can be expressed as follows [55]: where   is the ASTER NDVI;   and   are NDVI values of the dense vegetation canopy and the bare soil in the study area, respectively.Table 2 presents the summary statistics for   .More details regarding the SEE calculation are given in [53,56,57].

Soft Data Estimated from the SEE
The SEE is used to provide soft information at 100m-scale in terms of interval values in the BME procedure.SEE has the moderate and high correlation with SM, but SEE represents SM of bare soil at a depth of 0-2 cm, which is inconsistent with WSN and PLMR sensed depth (0-4 cm); SEE fails to reflect the root-zone SM in the vegetated areas.Due to these certainties, SEE-derived SM was processed as soft data in this study.
Firstly, SEE was converted into SM using the Lee-Pielke model, which accounts for the transport of water from inner soil pores to the soil surface and for the soil hydraulic conductivity [22,58].The SM estimations obtained using the Lee-Pielke model represent a bare soil layer of 0-2 cm [51].Lee-Pielke model is expressed as follows [59]: where   represents the SM for bare soil.β representing the surface moisture availability is a coefficient in bulk transfer formula [51], and it is called "soil evaporative efficiency (SEE)" in [53,60].θ  is the SM at field capacity, which depends on soil type and wind speed [61].θ  as an empirical parameter was assumed to be 0.3 or 0.35 m 3 •m −3 in previous works [58,62].However, considering that θ  varies with soil type and wind speed, Merlin, et al. [53] proposed that θ  can be obtained using SEE and SM by least square regression.This is achieved by inverting Equation ( 11) and analytically expressing θ  as a function of SEE and SM as follows: where   was represented by in situ measured SM. β represented the , which was equal to ASTER-derived SEE.In this study, in situ measurements at 4cm depth were used in Equation ( 12) considering the small bias and good correlation of SM at a depth of between 0-2 cm and 0-4 cm as Figure 3 shown.Considering the uncertainties of SEE-derived SM, the prediction interval of   can be estimated as following: where   is the prediction interval corresponding to   ,  −2 is the critical value of the t-distribution with (n−2) degrees of freedom and a confidence level of 95% [18].n is the number of WSN measurements. ⋅ is the standard deviation of the regression error of Lee-Pielke model,  ̅̅̅̅̅̅ is the mean of , and   is the sum of the squares of the  deviations [18].Figure 4 displays the SM-SEE relationship and corresponding to the prediction interval for each period.The accuracy of   retrieved from  is also displayed in Figure 4.

BME Estimation
In BME algorithm,   as soft information was merged with WSN measurements as hard data to estimate SM.The SM predictions ̅  can be expressed as [25]: where  * (  |  ,  ℎ ) is a posterior probability density function (PDF). ℎ ,  , and   denote the values of the hard, soft data, and unknown values at the estimation locations, respectively.Here,  ℎ is WSN-observed SM, and   is equivalent to the   .In this study, the SM in each ASTER-pixel area was considered to be homogeneous, because intensive field campaigns or higher resolution (finer than 100 m) remote sensing data are unavailable to analyze the SM heterogeneity within ASTER pixels.Hence, the WSN-observed SM data were used to represent the referenced values at 100 m-scale.
In the formulation of Bayesian law,  * (  |  ,  ℎ ) at each estimation locations updates from the prior PDF   (  ) as Equation ( 15): where   consists of a vector of   ,  ℎ and   .  (  ,  ℎ ) is the joint PDF of the hard and soft data.  (  ), the prior PDF of   , is expressed by a joint PDF rather than an actual prior PDF in the BME framework.Hence,   (  ) should be obtained first, which is calculated by general knowledge  , based on the maximum entropy theory. is expressed in mathematical formula as Equation ( 16): where the expectations ℎ ∝ is the moments of general knowledge, which is mathematically expressed by the covariance function (  ) as Equation (17). ∝ (  ) is expressed by the covariances (  ,   ) of random variables at location   and   as Equation (18).The covariances can be calculated by Equations ( 5)- (7).
∝ (  ) = (  ,   ) (  ) is mathematically expressed as [24]: where  0 is a Lagrange multiplier and  is a vector of the Lagrange multiplier.In order to maximize ℎ ∝ based on the maximum entropy theory,  is calculated by substituting Equations ( 17)- (19) into Equation (16).Then  is inserted back into Equation ( 19) to obtain   (  ).

Predictions Obtained by Merging WSN, SEE and Irrigation Statistics (Method III)
Irrigation leads to the strong spatial heterogeneity of SM in the experimental area [24].However, in the framework of BME, the SM distribution is treated as a spatial random field (SRF) [65], where SM is a regionalized random variable.Hence, under the requirements of SRF and the second order stationary assumption of covariance, the anisotropic variability of SM is necessary to be removed.Irrigation statistics were used to construct SM trends for the removal of SM anisotropy.The spatial SM distribution can be decomposed into a primarily spatial trend (  ) and a random auto-correlated residual component (  ) as Equation (20) expressing: where   can be further divided into the SM trend in the irrigated portion denoted by   , and the SM trend in the non-irrigated portion denoted by  − as following: and  − can be estimated using data from the irrigation districts combined with WSN-observed SM as follows: (1) Identify the irrigated and non-irrigated areas.Although the frequencies and durations of irrigation are not uniform among different fields, the average irrigation interval of a field is five days according to the irrigation statistics.The irrigation events corresponding to ASTER overpass dates do not influence the SM pattern of those days, because each ASTER overpass occurs at noon local time rather than at night when irrigation events occur.Therefore, the irrigated areas with respect to ASTER overpass times are defined as the fields in which irrigation events occurred during the five days before the corresponding ASTER overpass.The rest of the experimental area consists of non-irrigated areas.For example, for the SM obtained from the ASTER overpass conducted on 30 May, the SM redistribution and SM trends are expected to be dominated by the irrigation events that occurred from 25 May to 29 May.The overpass dates of the ASTER data and the periods influenced by irrigation are displayed in Table 2.
(2) The values of   are represented by the average WSN observations associated with the irrigated areas according to the ASTER overpass times, and the values of  − are the averages of the synchronous WSN measurements acquired in the non-irrigated areas.SM trends constructed by irrigation statistics in the three investigated periods are displayed in Figure 5.
After performing the above procedure, the   results were applied to remove the SM trends of the WSN SM data and the   results calculated in Section 3.2.2, and the residuals of the WSN SM and   data were obtained.The residuals of the WSN SM and   data were regarded as hard and soft data, respectively.BME predictions were obtained by applying both sets of residuals in Equations ( 14)- (19).  can be obtained as the sum of the BME predictions and   .

Predictions Obtained by Merging the WSN, SEE, Irrigation Statistics, and PLMR SM Data (Method IV)
The microwave brightness temperature data from PLMR were well correlated with SM [42,66,67].In vegetated areas, the microwave brightness temperature may be more closely related to SM than some optical dryness indices (e.g., TVDI) [19,68].Therefore, PLMR-derived SM at 700m were introduced to improve the accuracy of BME estimations.A downscaling approach proposed by Merlin et al. [53] was used to transform the PLMR SM into SM at 100 m resolution.The downscaling approach was based on high-resolution optical data derived SEE and physically-based model predictions of SEE.Considering the uncertainties caused by the downscaling method and the discrepancy between PLMRscale and 100m-scale, the downscaled SM at 100 m from PLMR data was processed as soft data.The downscaling method was formulated as [53]: where   is the PLMR SM product and   is the corresponding downscaled SM product at 100m-scale.The definition of   is the same as that ofin Equation ( 1  When considering the uncertainty resulting from the downscaling approach,   can be transformed into soft data intervals [70].The variance  2 of residuals obtained from the downscaling model were calculated.Then, the upper (  ) and lower (  ) bounds of the prediction interval were calculated using Equations ( 23) and (24): Next, the SM trends of   and   were removed using the   results constructed from the irrigation statistics shown in Section 3.3, and the residuals of   and   were obtained.The residuals of   and   were treated as soft data, and the residuals of WSN-observed SM were treated as hard data.These data were applied in Equations ( 14)-( 19) to obtain BME predictions.Next, the final Method IV results were obtained by summing these BME predictions and the   .

SM Variograms and Residuals
Figure 7 shows the variograms and residuals of the WSN-observed SM data.The corresponding parameters of the fitted variogram models are displayed in Table 3.The variograms of WSN SM for the three investigated periods exhibit irregular fluctuations as separation distance increases, which reflects the non-stationary and anisotropic of WSN SM.It can be attributed to that the local irrigation may spatially redistribute the SM, and damage the correlated structure of the SM variograms.This would result in inaccuracies of the fitted variogram models.The WSN-observed SM variograms for the three periods exhibit small nugget effects (0.01-0.02) and large sill variances (37-51.1),as shown in Table 3.The sills of WSN SM for 30 May and 24 June are larger than those of 10 July, which indicates that the spatial variations of SM for 30 May and 24 June are larger than those of 10 July.This result is also indicated by the smaller ranges of the SM data for 30 May and 24 June than that for 10 July.The smaller spatial variation of SM on 10 July indicates the stronger spatial autocorrelation of SM.However, after removing SM trends constructed from the irrigation statistics, the residual variograms for the three periods exhibit smaller ranges of spatial dependence and smaller sills compared with the variograms of the WSN SM data.Therefore, the irrigation statistics could be used to successfully remove the principal SM trends to obtain isotropic variograms of the residuals and reasonable variogram models, which can be converted into covariances for application in the BME framework.

Spatial Estimation of Soil Moisture
Figure 8 displays the four method estimations for the three periods.The white areas in the maps represent the masked areas that contain buildings and are not correlated with the surrounding SM distribution.Among the three periods, larger spatial variations in SM occurred on 30 May and 24 June, and smaller variations were observed on 10 July.Both Method I (Figure 8a-c) and Method II estimations (Figure 8d-f) show smooth changes in SM, such that the boundaries of the irrigated fields are unclear, although the Method II estimations show larger areas with high SM values.However, the irrigated regions are clearer in the Method III estimations, as shown in Figure 8g-i, because merging the irrigation statistics into the data used for the predictions provided helpful information regarding the irrigation districts.(a-c), (d-f), (g-i), j are the method I, II, III and IV estimations for three periods, respectively.The white areas in each map represent the masked areas that contain buildings.
Compared with Method I, II, III (Figure 8c,f,i) for 10 July, a larger spatial variance in SM is displayed in the Method IV estimations (Figure 8j).This result indicates that the PLMR-derived SM data provided more spatial SM information for improved BME estimation during this moderate vegetated period.Generally, Method III and Method IV predictions display more information regarding the spatial pattern of SM than Method I and Method II results, especially in large field regions with sparse WSN observations.The quantitative validation of the BME estimates is discussed in the next section.

Validation of BME Estimations
The cross validation (leave-one-out) method was used to evaluate the BME estimations [71].Forty-eight in situ measurements were used to validate the predictions.The scatterplots in Figure 9 present the cross validation results for the four methods estimations for the three periods.Quantitative results in terms of the root mean square difference (RMSD), bias, the coefficient of determination (R 2 ), and slope between the BME estimations and in situ measured SM are presented in Figure 9.
As shown in Figure 9a-c, all Method I estimations show the poor predictions with large RMSD and small R 2 and slopes.Compared with Method I estimations, Method II results obtained for the low-vegetated periods (30 May and 24 June) exhibited moderate correlations with in situ SM but low slopes, as shown in Figure 9d,e.These findings indicate the inclusion of SEE in Method II provide more information related to SM than Method I, but Method II overestimates SM at low values and underestimates SM at high values.
The Method II results for the moderately vegetated period (10 July) exhibit a poor correlation (0.45) with in situ SM, whereas the slope is only 0.26 as shown in Figure 9f.This result potentially occurred due to the failure of SEE to reflect SM in the vegetated areas.As shown in Figure 4, SEE as soft data is only weakly correlated with SM on 10 July, which increases the uncertainties in the Method II prediction.In addition, the biases of the Method II results for the three periods are −0.0083,−0.0065 and −0.002 m 3 • m −3 indicating that Method II predictions underestimate the overall SM in the study area.These findings may be explained by the overall underestimation of the SM prediction intervals calculated from SEE.The RMSD of all Method II aregreater than the acceptable RMSD (0.04 m 3 • m −3 ) for SM estimation missions [66,67].The reasons for the unacceptable RMSD may be that Method II estimations just represent SM for bare soil, rather than SM for bare and vegetated soil at 100m scale.Another reason could be the inaccurate SM variogram models, which provide biased spatial structure information as prior knowledge for auxiliary interpolation in the BME algorithm.Generally, the accuracy of Method II predictions for all three periods is not satisfactory, potentially due to the inaccuracies of the variogram models and soft data from SEE.
Compared with Method II, R 2 and slope between Method III and in situ SM measurements obtained for the low-vegetated periods (30 May and 24 June) are improved significantly, as shown in Figure 9g,f.In particular, the extreme SM values can be better estimated by Method III predictions.As shown from a comparison of the error statistics for the Method II and Method III predictions, the RMSD decreased from 0.0461 to 0.0324 m 3 • m −3 for 30 May and from 0.0571 to 0.0331 for 24 June, whilst R 2 increased from 0.39 to 0.62 and from 0.26 to 0.66, respectively.From the bias indicators, Method III for 30 May and 24 June provide excellent unbiased estimations.This result indicates that merging the irrigation trends data into the datasets used to generate the predictions improved the accuracy of BME estimations.In summary, compared with Method III in the low vegetated periods, Method III caused the overall RMSD to decrease from 0.052 to 0.033 m 3 • m −3 , and the overall R 2 and slope between the predicted and in situ measured SM increase from 0.32 to 0.64 and from 0.38 to 0.82, respectively.(a-c), (d-f), (g-i), j are the validation of predictions from method I, II, III and IV for three periods, respectively.Units for RMSD and Bias are % Nevertheless, for the moderately vegetated period (10 July), no significant improvement relative to Method II results (Figure 9f) was achieved when using the Method III (Figure 9i).However, R 2 and slope between the Method IV results and in situ SM are superior to those for the Method II and Method III results, as shown in Figure 9f,i,j.The RMSD decreased from 0.0411 m 3 • m −3 for Method III to 0.0342 m 3 • m −3 for Method IV, whereas R 2 increased from 0.18 to 0.4 and the slope increased from 0.41 to 0.69.These results indicate that merging PLMR-derived SM data into the data used for BME predictions can improve the accuracy of SM estimations.In particular, Method IV enhanced the estimation of the extreme values, as shown in Figure 9j.This suggests that PLMR-derived SM data can provide additional SM information during moderately vegetated periods when the SEE fails to sufficiently reflect the SM pattern.Finally, these results suggest that merging WSN data, ASTER data, irrigation statistics, and PLMR-derived SM data for predictions could yield better SM estimates for low-to-moderately vegetated periods.

Method I and Method II
It has been demonstrated that SM estimations can be improved by integrating multi-resource observations over heterogeneous cropland.However, the poor predictions from Method I indicate the sparse point-observations provide the limited SM information in heterogeneous cropland.An optical SM index, the SEE can be incorporated into Method II to provide high resolution SM information.Method II estimations are better related to in situ measurements than Method I, but still have poor prediction accuracies.The reasons for the poor accuracies could be that Method II is valid for bare soil and fails to estimate SM in vegetated areas.Thus, Method II could lead to the prediction errors in the vegetated areas, especially for 10 July under moderate vegetation cover.In future works, some optical indices such as evaporative fraction accounting for the root-zone SM in vegetated areas should be used in BME algorithm.
The calculation of the SEE assumes that the canopy skin temperature is relatively uniform throughout the experimental area.In fact, the canopy skin temperature varies with various factors (e.g., leaf water content), especially when SEE is applied at large scales, such as the SMOS-scale (~40 km).Additionally, the triangle approach used to estimate SEE is based on a number of assumptions, which makes it sensitive to various sources of error, as described in [72,73].

Method III
The calculated covariance when using variograms is a significant prior knowledge that can be used to express the spatial correlations and dependencies applied for interpolation in the BME algorithm.When irrigation statistics are used in BME predictions, Method III results offer the satisfactory prediction accuracy during low-vegetated periods.The predominant SM trends can be successfully removed by including irrigation statistics in the analysis, which results in reasonable variograms.Compared with the SM trends constructed based on remotely sensed variables in [19], the SM trends constructed based on irrigation statistics eliminate the problem of the failure of remotely sensed SM trends to remove extreme SM values.Moreover, the inclusion of irrigation information allows the boundaries of irrigated and non-irrigated areas to be described clearly (Figure 8), which is difficult to achieve in BME predictions that are based only on merging sparse point observations and remotely sensed variables.However, two assumptions regarding the application of irrigation statistics should be noted: (1) the period of influence of irrigation events is assumed to be five days.The length of this period should be further investigated by analyzing the temporal dynamic of SM caused by irrigation; (2) The residual SM variograms are assumed isotropic.However, the irrigation statistics can be used to remove the primary SM trends and the residuals could remain anisotropic in local regions due to the strong heterogeneity of the SM.These errors in the fitted residual variograms could be reduced using the local anisotropy model [74].

Method IV
Merging PLMR SM data into the datasets used for BME predictions successfully yields Method IV results with satisfactory accuracy for the moderately vegetated periods.In vegetated areas, the TB observed in the L-band are more sensitive to variations in SM than optical data, because vegetation attenuation is smaller at lower frequencies and at horizontal polarization [75].The ability to achieve good performance using Method IV depends on the acceptably accurate PLMR SM products and the downscaling method.In this study, the downscaled SM from PLMR SM, treated as soft data, were beneficial as auxiliary data for improving BME estimations during the moderately vegetated periods.However, the bias from the downscaled results may bring in uncertainties to Method IV estimations.As shown in Figure 6, the downscaled results showed the underestimation at low values, which could result in the underestimation of Method IV (Figure 9j).Likewise, the overestimation of the downscaled results may lead to the overestimation of Method IV.Additionally, Method IV was applied for only one day (10 July 2012) in this study, and it could be further evaluated in the days of other seasons.
In the downscaling method [57], it is assumed that the SM distribution within each PLMR pixel is represented by SEE at ASTER scale.Namely, the heterogeneity of a PLMR pixel can be described by the SEE at 100 m-scale.However, the errors incurred by this method will increase when SEE at 100 m-scale are not accurate enough to represent the heterogeneity within the PLMR pixels.Additionally, the downscaling method requires the completeness of ASTER-derived SEE within PLMR pixels.This requirement could limit the application of the method in areas where data at 100 m-scale are missing.The development of optical-based downscaling approaches of microwavederived SM is still in its infancy, and further evaluation studies are needed [57,76].
The SM values of buildings at ASTER-and PLMR-scales are considered as zero m 3 m −3 in the downscaling method.This relies on one implicit assumption that SM of buildings derived from PLMR using the L-MEB model is zero m 3 m −3 .The sensitivity analysis of the L-MEB model has indicated that SM of building (or rocks) is nearly zero m 3 m −3 [77,78].However, it is noted that the intra-pixel was assumed to be homogeneous in L-MEB model [79,80].The effects of different land types within pixels on TB values should be discussed for the appropriate application of L-MEB model over heterogeneous land surfaces.
Note also that the mismatch of sensing depth between Lee-Pielke model estimations and WSN/PLMR observations.The WSN/PLMR sensing depth is approximately 4-5 cm, whereas Lee-Pielke model only estimates the top 2 cm soil layer.Although Lee-Pielke model estimation was processed as soft data in BME method, the sensed depth mismatch between Lee-Pielke model estimation and PLMR/WSN observations could lead to the large error of Method II, III, and IV.The issues could be solved using the inversion models which can derive SM at 0-4 cm depth using optical data.It is also noted that Lee-Pielke model fails to estimate SM larger than filed capacity as Equation (11) shown.Thus, soft data from Lee-Pielke model could result in the underestimation of SM at high extreme values, particularly in the irrigated areas.
It should be noted that WSN-observed SM were considered to represent the SM at 100m-scale in this study.However, the footprint of WSN measurements is far smaller than that of ASTER.Due to the strong heterogeneity of SM, the mismatching footprint could lead to the errors in the validation of BME estimations.Thus, the representative footprint of in situ measurements should be investigated, and strategies for upscaling point observations could be further introduced into the BME algorithm [15].

Conclusions
Results indicate that BME predictions obtained by merging only SEE and WSN measurements are well correlated with in situ SM, but underestimate the overall SM in the low-vegetated periods.After irrigation statistics were added to the merged datasets, the non-stationary and anisotropic of the SM data were successfully removed, and Method III results clearly outlined the irrigation domains.During the low-vegetated periods, Method III decreased the overall RMSD from 0.052 m 3 • m −3 to 0.033 m 3 • m −3 , and increased the overall R 2 and slope between the predicted and in situ measured SM data from 0.32 to 0.64 and from 0.38 to 0.82, respectively.Furthermore, adding PLMR SM data into the merged datasets used for BME predictions improved the performance of the BME predictions during the investigated moderately vegetated period.The RMSD decreased from 0.0411 for Method III results to 0.0342 m 3 • m −3 for Method IV, whilst the R 2 and slope increased from 0.18 to 0.4 and from 0.41 to 0.69, respectively.
In summary, the issues on SM prediction in heterogeneous cropland are noticeably improved by merging WSN, ASTER-derive SEE, irrigation statistics, and PLMR data in the framework of BME.However, this study only estimates the spatial SM pattern without considering spatiotemporal auxiliary data.In future studies, the time series observations from multi-resource will be applied to obtain SM estimates with high spatiotemporal resolution in the framework of BME.

Figure 1 .
Figure 1.Overview of the study area and soil-moisture monitoring network.The right image is the ASTER L1B VNIR image covering the study area on 24 June 2012.The RGB components are channels 1 (0.56 μm), 3 (0.81 μm), and 2 (0.66 μm), respectively, with 15 m resolution.

Figure 2 .
Figure 2. Flowchart of SM estimation in irrigated areas using multi-resource observations.

Figure 3 .
Figure3.Scatters of SM measurements at a depth of 2cm versus that of 4cm.The measurements at 2 cm and 4 cm corresponding to ASTER overpass were provided by automatic meteorological stations (AWSs)[63,64].There were 8, 12, 12 AWSs measurements available in the study area for 30 May, 24 June, and 10 July, respectively.

Figure 4 .
Figure 4. SM-SEE scatter plots, linear regression lines and associated prediction intervalsfor soft data.(a-c) correspond to 30 May 2012, 24 June 2012, and 10 July 2012, respectively.The dashed lines represent the regression line between the transformed SEE and in situ measurements.The solid lines indicate the prediction interval for   at a confidence level of 95%, where   =  −1 (1 − 2 0.5 ) / .R represents the correlation coefficient between the transformed SEE and the in situ SM.RMSE is root-mean-squared error of   , which is presented in units of %•m 3 • m −3 .

Figure 5 .
Figure 5. SM spatial trends modeled from the irrigation statistics and WSN data with 90m spatial resolution.The white areas in each map represent masked areas that contain buildings.
). 〈  〉  denotes the average of   within the PLMR grid box.derivative of SM with respect to SEE.In this study, we obtained (   )  via optimization using a least square method[69].The SEE of buildings are considered zero.As shown in Figure6, the   results predicted using the downscaling method exhibit a RMSE of 0.048 m 3 • m −3 and a correlation coefficient of 0.61 between the downscaled and in situ measured SM data.

Figure 6 .
Figure 6.The scatter plots between in situ measurements and the downscaled SM from PLMR SM products for 10 July.R represents the correlation coefficient between in situ observations and the downscaled SM.

Figure 7 .
Figure 7. SM semi-variograms (red circles), the correspondingly fitted SM variograms (red solid line), the semi-variograms of SM residual (blue squares) and the correspondingly fitted residual variograms (blue solid line).

Figure 8 .
Figure 8. SM spatial distribution estimated by four methods with 90 m spatial resolution.(a-c),(d-f), (g-i), j are the method I, II, III and IV estimations for three periods, respectively.The white areas in each map represent the masked areas that contain buildings.

Figure 9 .
Figure 9. Cross-validation and accuracy statistics of the results from the four methods.(a-c),(d-f), (g-i), j are the validation of predictions from method I, II, III and IV for three periods, respectively.Units for RMSD and Bias are %

Table 1 .
Summary of WSN data used in the study.
Note: The units for Max, Min, Mean and SD are m 3 • m −3 .

Table 2 .
Summary of the ASTER, PLMR and irrigation statistics.
Note: Units for   and  − are m 3 • m −3 .Max_Fr, Min_Fr, and Mean_Fr represent the maximum, minimum, and mean of fraction vegetation cover, respectively, in the study area.

Table 3 .
Parameters of the variogram models.