Estimating Time Series Soil Moisture by Applying Recurrent Nonlinear Autoregressive Neural Networks to Passive Microwave Data over the Heihe River Basin, China

Abstract

A method using a nonlinear auto-regressive neural network with exogenous input (NARXnn) to retrieve time series soil moisture (SM) that is spatially and temporally continuous and high quality over the Heihe River Basin (HRB) in China was investigated in this study. The input training data consisted of the X-band dual polarization brightness temperature (TB) and the Ka-band V polarization TB from the Advanced Microwave Scanning Radiometer II (AMSR2), Global Land Satellite product (GLASS) Leaf Area Index (LAI), precipitation from the Tropical Rainfall Measuring Mission (TRMM) and the Global Precipitation Measurement (GPM), and a global 30 arc-second elevation (GTOPO-30). The output training data were generated from fused SM products of the Japan Aerospace Exploration Agency (JAXA) and the Land Surface Parameter Model (LPRM). The reprocessed fused SM from two years (2013 and 2014) was inputted into the NARXnn for training; subsequently, SM during a third year (2015) was estimated. Direct and indirect validations were then performed during the period 2015 by comparing with in situ measurements, SM from JAXA, LPRM and the Global Land Data Assimilation System (GLDAS), as well as precipitation data from TRMM and GPM. The results showed that the SM predictions from NARXnn performed best, as indicated by their higher correlation coefficients (R ≥ 0.85 for the whole year of 2015), lower Bias values (absolute value of Bias ≤ 0.02) and root mean square error values (RMSE ≤ 0.06), and their improved response to precipitation. This method is being used to produce the NARXnn SM product over the HRB in China.

1. Introduction

Surface soil moisture (SM), which is one of the most important factors in regional and global water cycle processes, has played an increasingly important role in the earth sciences (e.g., meteorology, hydrology, agriculture, and biomass analysis). The estimation of SM from remote sensing data is one of the most feasible ways to generate surface SM products at regional scales and especially at global scales. Passive microwave remote sensing offers unique advantages for the retrieval and evaluation of SM [1,2,3]. Much evidence has been provided by previous studies [4,5,6,7] that passive microwave remote sensing techniques offer an efficient way to estimate surface SM at large scales.

Several spaceborne passive microwave sensors are available, including primarily the L-band radiometer onboard the Soil Moisture Active and Passive (SMAP) satellite, the Microwave Imaging Radiometer with Aperture Synthesis (MIRAS) onboard the Soil Moisture and Ocean Salinity (SMOS) satellite, the Advanced Microwave Scanning Radiometer II (AMSR2) onboard Earth from the Global Change Observation Mission-Water 1st (GCOM-W1) spacecraft, the Tropical Rainfall Measuring Mission (TRMM) Microwave Imager onboard the TRMM satellite, and the Special Sensor Microwave/Imager (SSM/I) onboard the Defense Meteorological Satellite (DMS). The ongoing operation of these major satellite sensors is intended to estimate parameters that are vital for the hydrological and meteorological sciences, such as surface SM. These SM products have played an important role and have shown great potential in climatic and meteorological studies [8,9,10]. However, data from these passive microwave sensors are not ideal for obtaining SM from regions characterized by extensive vegetation cover [11] or from regions that experience severe man-made radio frequency interference (RFI) or extraterrestrial galactic radiation [12]. In addition, surface SM products are originally discontinuous, due to the footprints and repeat intervals of the satellite orbits, restrictions in instrument design and limitations of the inversion algorithms. Thus, the defects described above present obstacles for the application of surface SM in the study of geoscience.

In terms of the performance of passive microwave remote sensing of SM in China, two products derived from AMSR2, namely the Japan Aerospace Exploration Agency (JAXA) SM product [2,13,14] and the Vrije Universiteit Amsterdam and NASA Goddard Space Flight Center (VUA-NASA) SM product [15,16], provide superior results [6,17]. The VUA-NASA SM product is hereafter referred to as the Land Surface Parameter Model (LPRM) because its retrieval methodology uses the LPRM. In recent research [6,17,18], the JAXA and LPRM SM products, along with the Global Land Data Assimilation System (GLDAS) SM products [19], were evaluated in the cold semi-arid area of the Tibetan Plateau (TP). The results showed that the three SM products have different advantages and disadvantages, which leads to the conclusion that SM products must be improved if they are to be used in the cold, semi-arid areas of China.

Over the past few decades, intelligent inversion algorithms have been used in several studies for the retrieval of SM from microwave radiometers [20]. Neural network is a popular and effective tool that is based on pre-computed and prior information, which has been widely used in many fields, such as pattern recognition [21,22], pattern classification [23,24,25] and parameter estimation [26,27,28]. One of the benefits of neural networks is their capacity for classification and function approximation [29]. The use of neural networks in parameter estimation is generally accepted because neural networks can simultaneously handle nonlinear dependencies and complex statistical problems [30]. Liou et al. [31] integrated a dynamic neural network, the dynamic learning neural network (DLNN), with a biophysically based land surface process/radio-brightness (LSP/R) model to retrieve SM and other land surface parameters from simulated data. New error propagation learning back propagation (EPLBP) neural networks were then built for SM retrieval from 1.4 GHz, 6.9 GHz and 10.65 GHz multi-band dual polarization brightness temperature (TB) after training with passive microwave data [32,33]. Del Frate et al. [34] used BP neural networks to retrieve SM and other agricultural variables from multi-frequency radiometer TB. Chai et al. [35] use land cover, soil texture, soil temperature, and canopy temperature as auxiliary data to provide a priori information. Additionally, the integration of independent Normalized Difference Vegetation Index (NDVI) and land surface temperature (LST) data sets into the input of the neural network matrix generated substantially improved results. Santi et al. [36] utilized BP neural networks to obtain a SM map of the entire world from AMSR-E TB. The results were tested in Mongolia and Australia with a root mean square error (RMSE) no lower than 0.065 m³/m³ and a robust determination coefficient (R²) equal to 0.82. In recent years, joint inversions using active and passive microwave data have been attempted [37,38]. These attempts used both TB and the backscattering coefficient as the main variables of the input data set. Although all previous research has shown the promise of neural networks, there is still a lack of applications of this method in practice. In addition, all of the existing studies have been conducted using BP neural networks, which have shortcomings including slow learning speeds, a tendency to find locally optimal solutions, and other challenges. Thus, the application of a new feedback neural network to more input variables deserves consideration.

In this work, the possibility of combining AMSR2 X-band dual polarization TB with other auxiliary data to reconstruct spatially and temporally continuous, high-quality time series SM was investigated within the Heihe River Basin (HRB) in China. Furthermore, neural network based on nonlinear auto-regressive models with exogenous inputs (NARXnn) was implemented and validated. The predicted SM from the trained NARXnn was validated with in situ measurements of SM obtained using data collectors. The predictions were also compared with GLDAS SM during the entire validation period. This paper is organized as follows. In Section 2, the study area and data set are described. In Section 3, the details of the architecture of the NARXnn and expositions of the calculation of the fused SM are provided. In Section 4, the development of the NARXnn is summarized, and the results are validated in Section 5. Section 6 draws conclusions and presents several analyses.

2. Study Area and Data

2.1. Study Area

The study area was bounded by the latitudes 37.5°–43°N and the longitudes 97°–102°E and was located in the HRB in China (Figure 1). It covered an area of approximately 450 km × 555 km, with elevations ranging from 790 m to 5713 m. The area is substantially influenced by the southwestern monsoon and has an annual average temperature of approximately −1 °C and annual precipitation of 300–500 mm. Thus, the local climate is semiarid and cold alpine [39]. The land cover types (i.e., water, low vegetation, and barren and sparse vegetation) were obtained, rescaled and reclassified from MODIS land cover data (MCD12C1). The proportion of low vegetation was 65.5%, and the proportion of barren areas and sparse vegetation was 34.1%. In total, a rectangular area that contained four hundred and forty 0.25° × 0.25° grids was selected as the study area.

2.2. In Situ Measurements

A major research plan entitled “Integrated research on the eco-hydrological process of the HRB” has been launched by the National Natural Science Foundation of China (NSFC) and the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) was initiated in 2010. Three key experimental areas within the HRB (a cold experimental area in the mountain cryosphere of the upper reaches; an artificial oasis experimental area in the middle reaches; and a natural oasis experimental area downstream) were selected in which to conduct intensive and long-term measurements [40]. In situ measurements within the three key experimental areas were selected for validation, and these data were kindly provided by the Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/).

The first network, named the “Watershed Allied Telemetry Experimental Research observation network” (WATER-NET), was located within the upstream regions of a sub-basin (the Babao River basin) of the Heihe River [41]. The network extended approximately from 37.75°N to 38.5°N and from 100°E to 101.25°E. This area contained predominantly bare soil and low vegetation. There were 40 sites within this network where SM measurements were automatically collected at 5-min intervals and at a depth of 4 cm. The in situ SM measurements collected by this network covered the period from 1 July 2013, to 31 December 2014 [41,42,43]. The probe on each site is a Hydro Probe II (HP-II) sensor. Four typical 0.25° × 0.25° grids (hereafter referred to as WATER-NET 1–4) that contained low vegetation were chosen for validation. The grids are outlined in dark colors and numbered in Figure 1. These four grids were selected because each one contained at least five in situ measurement sites.

The second network, named the “Hydrometeorological observation network” (HM-NET), was established within the middle and downstream regions of the HRB. The automatic data collectors collected SM data using in situ probes at a depth of 4 cm and 10-min intervals during the period from 1 December 2013 to 31 December 2014 [44,45,46]. The soil moisture sensors for major sites are CS616. For minor sites, the sensors are ECH2O-5 (Decagon Device). Five of these sites were located in the middle stream region of the HRB, all within a 0.25° × 0.25° grid (hereafter referred to as HM-NET 1) that extended from 38.75°N to 39°N and from 100.25°E to 100.5°E. Within the lower stream area, there were five sites located around the border of two 0.25° × 0.25° grids that extended from 41.75°N to 42.25°N and from 101°E to 101.25°E. Grid average values were used in this study to better cope with the problems associated with poor site distributions. Hereafter, HM-NET 2 is used to represent the average value of the two grids. The grids are outlined in dark colors and numbered in Figure 1.

2.3. Data Sets for NARXnn

2.3.1. The Products of AMSR2

The AMSR2 was launched aboard the GCOM-W1 satellite on 18 May 2012 and has provided observational data since 3 July 2012 [17]. AMSR2 offers highly accurate measurements of the intensity and scattering of microwave emissions (http://suzaku.eorc.jaxa.jp/GCOM_W/index.html). The equatorial crossing times of GCOM-W1 are approximately 01:30 (descending pass) and 13:30 (ascending pass).

In this work, the JAXA and X-band LPRM SM products were selected as training data sets of NARXnn, whereas the AMSR2 X-band (10.7 GHz) dual polarization TB (hereafter referred to as TB-X-H and TB-X-V, respectively) and the Ka-band (37 GHz) V polarization TB (hereafter referred to as TB-Ka-V) were used as part of the training data sets. The latest versions of all of these products were used. The TB products and the JAXA SM product were downloaded from http://gcom-w1.jaxa.jp/index.html, and the LPRM SM product was downloaded from http://disc.gsfc.nasa.gov/hydrology/data-holdings. The ascending Level 3 gridded 0.25° JAXA SM product and the X-band LPRM SM product were chosen for use in this study.

2.3.2. Other Auxiliary Data

Global Land Satellite product (GLASS), TRMM, GPM and GTOPO-30 products were used in this work. The GLASS LAI data were downloaded from http://glcf.umd.edu. The TRMM and GPM products were downloaded from the EARTHDATA website of NASA’s Earth Science Data Systems Program (https://search.earthdata.nasa.gov/), and the GTOPO-30 data were obtained from https://lta.cr.usgs.gov/GTOPO30.

The GLASS product was used to estimate daily Leaf Area Index (LAI) values. The GLASS products contain long time series of global land surface characteristics based on various remotely sensed data. These products include LAI, land surface albedo, broadband emissivity, downward shortwave radiation and photosynthetically active radiation [47]. The GLASS LAI product has an eight-day temporal resolution with a spatial resolution of 1 km in a sinusoidal projection system. Thus, the GLASS LAI data were resampled and clipped according to the gridded 0.25° resolution of the latitude-longitude system.

In this study, precipitation (PRC) values were extracted from the TRMM and GPM data. TRMM is a joint mission between NASA and JAXA that studies rainfall for weather and climate research. The TRMM 3B42 product (http://trmm.gsfc.nasa.gov/3b42.html), which has a 3-h temporal resolution and a spatial resolution of 0.25° by 0.25°, was selected. Daily averaged processing was performed before the training process was undertaken. The GPM IMERG L3 precipitation data (https://www.nasa.gov/mission_pages/GPM/overview/index.html) were delivered at a resolution of 0.1° × 0.1°. Here, the L3 final run daily precipitation data were extracted and then resampled to a 0.25° × 0.25° grid to enable coordination with the AMSR2 observations.

Elevation values were extracted directly from the GTOPO-30 product. GTOPO-30 is a global digital elevation model (DEM) that is the result of a collaborative effort led by the U.S. Geological Survey’s EROS Data Center. The elevations in GTOPO-30 are regularly spaced at 30 arc seconds (approximately 1 km). The GTOPO-30 product has been used in many regional and continental studies, such as studies of climate modeling, continental-scale land cover mapping, extraction of drainage features for hydrologic modeling [48], and geometric and atmospheric correction of medium- and coarse-resolution satellite image data [49]. Here, the GTOPO-30 data were stitched, resampled and clipped in line with the 0.25° × 0.25° grid.

2.4. GLDAS

The GLDAS SM data were also available from http://disc.gsfc.nasa.gov/hydrology/data-holdings. The GLDAS_NOAH025SUBP_3H (Version 001) Noah Land Surface Model L4 3 hourly 0.25° × 0.25°-degree subset (Noah V2.7) data product was used as the reference SM.

GLDAS aims to generate optimal fields of land surface states and fluxes with methods that integrate satellite-based and ground-based measured data products using advanced land surface modeling and data assimilation techniques [50]. GLDAS drives multiple, offline land surface models (LSMs), integrates many measurement-based data, and executes globally at different resolutions (0.25°, 0.1° and 40 km) [51]. Currently, GLDAS drives four LSMs: Catchment, Noah, the Community Land Model, and the Variable Infiltration Capacity model. Only the data from the Noah model were selected, partly because Noah provides a 0.25° × 0.25° product, whereas the other three products use a 1° × 1° system. Lastly, the daily GLDAS SM values were calculated by averaging the 12:00 and 15:00 values to approximate the 13:30 value.

3. Methodology

In passive microwave remote sensing, the electromagnetic energy received by a radiometer is mainly influenced by the observation frequency $F$, the observation angle $\theta$, the polarization $p$, the surface soil moisture $SM$, the surface temperature $LST$, the vegetation parameter $Veg$, and the terrain parameter $T$. This relationship can be expressed as:

$$TB=f\left(F, \theta ,p,SM,LST,Veg,T \right)$$

(1)

AMSR2 TB-X-H and TB-X-V were selected because the JAXA and LPRM SM products adopted in the current study are both obtained from X-band TB. The LST was replaced by the AMSR2 TB-Ka-V, which is explained in detail in Section 3.2. The GLASS LAI was chosen to represent $Veg$. The elevation ($H$) was selected to represent $T$ because elevation is the most significant factor determining SM variations over the HRB. In addition, as for the AMSR2 X-band TB, $\left(F,\theta ,p\right)$ were constants. Thus, for the purpose of SM retrieval, Equation (1) can be rewritten as follows:

$$SM=f\left(TB,LST,LAI, H\right)$$

(2)

NARXnn was used to explore the possibility of inferring time series SM from AMSR2 TB and other auxiliary data. Only daytime ascending orbit data were used because of their robust performance over areas experiencing droughts and semi-droughts in previous studies [6]. Unlike previous SM retrieval research that has used neural networks, in this study, the NARXnn was trained with the fused time series SM data derived from the JAXA SM and LPRM SM data products, instead of simulated or measured data, using a Bayesian regularization backpropagation training function.

3.1. NARXnn and its Configuration for SM Retrieval

NARXnns have been shown to perform well on problems involving long-term dependencies, and they are capable of simulating universal dynamical systems [52,53]. The feasibility of NARXnns as a nonlinear tool for time series modeling and prediction is key for carrying out long-term time series predictions [54].

The NARX model can be expressed as follows:

$$y\left(t\right)=f\left(y\left(t-1\right), y\left(t-2\right),\dots ,y\left(t-n_y\right),u\left(t\right),u\left(t-1\right),\dots u\left(t-n_u\right) \right)$$

(3)

where $t$ represents a certain time point. The next value of the dependent output signal $y\left(t\right)$ is regressed on previous values of the output signals $y\left(t-1\right), y\left(t-2\right),\dots ,y\left(t-n_y\right)$ and previous values of the independent (exogenous) input signals $u\left(t\right),u\left(t-1\right),\dots u\left(t-n_u\right)$. The NARX model was implemented using a feedforward neural network to approximate the function $f$. The output $y\left(t\right)$ feeds back to and impacts the function $f$, with the input $u\left(t\right)$ moving to the next time point $t+1$ with respect to a given time point $t$ [55]. The architecture of the NARXnn is shown in Figure 2. It incorporates a hidden layer of weights $W$ and $i$ neurons connected to an input layer and an output layer.

The trapped delay line (TDL) is a discrete element that allows a signal to be delayed by a number of samples. By configuring the TDL in advance, such as fixing the values of $i$ and $j$ in Equation (3) with the input data, the NARXnn can be established for time series forecasting [56]. Here, the TDL was set as a two-dimensional vector (0, 1), rendering the input signals $t-1$ and the target signal $t$. This means that the $SM$ at a time point $t$ is determined by its current state at $t$ as well as by the original $SM$ at $t-1$. As a result, Equation (3) can be expressed as follows:

$$y\left(t\right)=f\left(y\left(t-1\right),u\left(t\right)\right)$$

(4)

As for the inputs ($U$) and the output ($Y$), several steps were performed beforehand to transform all of the variables involved into a successive time series form:

$$\{\begin{array}{c} U= \left\{\left[u_1\right] \left[u_2\right]\dots \left[u_t\right]\right\}\\ Y= \left\{\left[y_1\right] \left[y_2\right]\dots \left[y_t\right]\right\}\end{array}$$

(5)

where $\left[u_t\right]$ and $\left[y_t\right]$ represent the input ($U$) and output ($Y$) values of the NARXnn at time point $t$, respectively. Since $U$ has an eight-day temporal resolution because of the eight-day temporal resolution of the GLASS LAI data product, day of year (DOY) was also considered as a variable in $U$. Thus, the NARXnn has a seven-dimensional input vector (i.e., TB-X-H, TB-X-V, LAI, TB-Ka-V, PRC, DEM and DOY) and a one-dimensional output vector (SM):

$$U_i= \left[\begin{array}{c}\begin{array}{c}TB-X-H\\ TB-X-V\end{array}\\ LAI\\ TB-Ka-V\\ PRC\\ DEM\\ DOY\end{array}\right]$$

(6)

$$Y_i= \left[SM\right]$$

(7)

Therefore, Equation (4) can be rewritten as:

$$SM\left(t\right)=f\left(SM\left(t-1\right), \left[\begin{array}{c}\begin{array}{c}TB-X-H\\ TB-X-V\end{array}\\ LAI\\ TB-Ka-V\\ PRC\\ DEM\\ DOY\end{array}\right] \left(t\right)\right)$$

(8)

The samples in the U and Y sequences were randomly and temporally divided to avoid overfitting or underfitting the NARXnn as follows: 70% were assigned to the training data set, and 30% were assigned to the test data set. The proportions (70% and 30%) were set because, under this circumstance, NARXnn performed more robustly than for other proportions. The training data set was presented to the network during training. The network was adjusted according to its error and was halted when the adaptive weight was minimized. The testing data set had no effect on training and, thus, provided an independent measure of network performance during and after training. The training function used here was Bayesian regularization backpropagation, which updates the weight and bias values according to the Levenberg–Marquardt optimization algorithm, minimizes a composite of squared errors and weights, and subsequently determines the correct combination to produce a network that generalizes well. Unlike the Levenberg–Marquardt function, which is widely used for network training, the validation process of the Bayesian regularization function stops and is disabled by configuring maximum validation failures to be zero so that training can continue until an optimal combination of errors and weights is recovered [57]. The weight/bias learning function here is a gradient descent algorithm with momentum weighting and a bias learning function, which calculates the weight change for a given neuron from the neuron’s input and error, the weight (or bias), learning rate, and momentum constant. The transfer functions of hidden layers and the output layer are set as a tan-sig transfer function [58] and a linear transfer function, respectively. The tan-sig transfer function and the linear transfer function are adapted here for the simple reason that they are useful for describing the dynamic varying range of SM by varying the input variables.

3.2. Replacement of LST by AMSR2 TB-Ka-V

Several studies [59,60,61] have shown that passive microwave TB-Ka-V have high transmission through the atmosphere and have a good linear relationship with LST. Thus, LST can be calculated from TB-Ka-V. Statistical empirical algorithms are a common method for finding linear relationships using a large number of in situ data. This method performs well under the assumption that the atmosphere and the underlying surface are relatively uniform.

Presently, there are several main statistical empirical algorithms for estimating LST using TB-Ka-V, such as that of Owe et al. [59]:

$$Ts=1.271\times T{B}_v^{37\mathrm{GHz}}-50.1$$

(9)

where $Ts$ represents LST, and $T{B}_v^{37\mathrm{GHz}}$ represents TB-Ka-V. These cases are the same as those described by De Jue et al. [60] and Holmes et al. [61], except that the coefficients are different. All of the three statistical empirical algorithms show that TB-Ka-V has a linear relationship with LST. Therefore, over the HRB, the quantitative SM estimates from the TB-Ka-V were as follows.

For the ascending orbit:

$$Ts=0.6552\times T{B}_v^{37\mathrm{GHz}}+100.3626$$

(10)

For the descending orbit:

$$Ts=0.6140\times T{B}_v^{37\mathrm{GHz}}+116.8851$$

(11)

These coefficients were obtained by comparing the AMSR2 TB-Ka-V with the in situ data collected during the period from 1 July 2013 to 31 December 2014 in the upper HRB (i.e., at WATER-NET 1–4), from 1 January 2013 to 31 December 2014 in the middle HRB (i.e., at HM-NET 1), and from 1 January 2014 to 31 December 2014 in the lower HRB (i.e., at HM-NET 2). However, in the current study, LST was replaced directly by the TB-Ka-V because the linear relationship does not change the NARXnn’s working pattern. The fitting coefficients were obtained only for the analysis in Section 3.3.

3.3. Fused SM from JAXA and LPRM

In this work, efforts were made to estimate the fused SM from a linear combination of the values stored in the JAXA and LPRM SM products as follows:

$$S{M}_{fused}= {\omega }_{JAXA}S{M}_{JAXA}+ {\omega }_{LPRM}S{M}_{LPRM}+a$$

(12)

where $S{M}_{fused}$ is the fused SM, $S{M}_{JAXA}$ represents the JAXA SM, and $S{M}_{LPRM}$ represents the LPRM SM. ${\omega }_{JAXA}$ and ${\omega }_{LPRM}$ are the weights of the JAXA SM and the LPRM SM, respectively, and $a$ is a constant coefficient.

Previous studies [6,17,18] have found that SM exhibits distinct seasonal variations in the cold, semi-arid areas of China. Figure 3 depicts a comparison of the JAXA and LPRM SM products, along with the in situ data and PRC, from four grids within the upper portion of the HRB (i.e., at WATER-NET 1–4) from 1 July 2013 to 31 December 2015. Overall, JAXA was fixed to a restricted range and displayed limited variability, whereas LPRM had a large range and always acted as a “high-jumper”. It should be noted that both products included periods of missing or incomplete data. Specifically, during the unfrozen season, LPRM showed robust correspondence with PRC and captured the entire seasonal cycle at different levels, whereas JAXA diverged from the in situ data but still bore a dynamic connection with it. During the freezing/thawing season, JAXA was nearly invariant, whereas LPRM was negatively correlated with the in situ data but was still positively correlated with PRC. During the frozen season, LPRM performed poorly and JAXA performed adequately at estimating the unfrozen water content.

Thus, the unfrozen seasons and frozen seasons, as well as the freezing/thawing seasons, were considered separately in this work. The boundaries of the unfrozen seasons, the freezing/thawing seasons and the frozen seasons were determined using the LST profiles of the HRB obtained from Equations (10) and (11). The yearly mean LST profiles (yellow lines for the ascending orbit and cyan lines for the descending orbit) from 2013 to 2015 and the mean yearly LST profile (the bold red line for the ascending orbit and the bold blue line for the descending orbit) are shown in Figure 3. According to these profiles, the time nodes of the unfrozen seasons and frozen seasons were set to the DOYs described in

Table 1. Coefficients of the unfrozen seasons and frozen seasons.

	${\mathit{\omega }}_{\mathit{J}\mathit{A}\mathit{X}\mathit{A}}$	${\mathit{\omega }}_{\mathit{L}\mathit{P}\mathit{R}\mathit{M}}$	$\mathit{a}$	RMSE (m3/m3)	R	Beginning (DOY)	End (DOY)
Frozen seasons	0.55795	0 1	0.04743	0.019	0.49	319	84
Unfrozen seasons	0.38058	0.35237	0.14116	0.067	0.76	104	298

¹ Fused SM was calculated only from JAXA SM during frozen seasons.

. It should be noted that the beginning of each unfrozen season was set to DOY = 84 because, based on the situation of the HRB, the timing of the melting point is not the intersection of the mean yearly LST profile and LST = 0 °C. Therefore, the point was set using comprehensive information from the LST profiles, the in situ SM profiles and the data recorded at the sites.

As for the unfrozen seasons and frozen seasons, the least squares method was used to calculate the coefficients. The coefficients were obtained by comparing JAXA and LPRM with in situ data during the period from 1 July 2013 to 31 December 2014 from the upper HRB (i.e., WATER-NET 1–4), from 1 January 2013 to 31 December 2014 from the middle HRB (i.e., HM-NET 1) and from 1 January 2014 to 31 December 2014 from the lower HRB (i.e., HM-NET 2). Coefficients and time nodes of the unfrozen seasons and frozen seasons are shown in Table 1. For the freezing/thawing seasons, a quadratic curve smoothing method was used to join the unfrozen seasons and frozen seasons.

The fused SM is shown in Figure 4 in orange and GLDAS SM is displayed in blue. Compared with JAXA and LPRM SM, the fused SM profile fit the in situ measurement profile fairly well and exhibited distinct seasonal features. GLDAS SM is consistent with the in situ data in profile during the unfrozen season despite some obvious underestimations, which is in agreement with previous research [6]. Therefore, in Section 5.2, NARXnn SM was compared with GLDAS SM for the unfrozen season across the whole HRB.

4. Development of the NARXnn

First, the data organization of NARXnn and the preparations made before training are stated. Second, the details of training and testing of NARXnn are explained. Lastly, the steps of the implementation of NARXnn SM are sketched briefly.

4.1. Preparations before Training

Some preparations were made before training. First, the fused SM was obtained from JAXA and LPRM. Then, the fused SM was transformed to a time series form, i.e., the output (Y). At the same time, all of the relevant variables were transformed into another time series form, i.e., the inputs (U). Given that all of the data were daily and continuous, except for the GLASS LAI product, which had an eight-day time step, the successive time series Y and U were both set as having an eight-day temporal sampling. Hence, the lengths of the original Y and U data series were both 92. A filtering process was then executed to determine the final successive time series Y, U and to ensure that they contained no missing data. In detail, the DOYs in the original time series with grids containing amounts greater than or equal to 300 were drawn out to generate a new final time series. In addition, the lengths of the final time series were 46. For the final time series, the grids in the same position over the HRB and with no missing data during the 46 DOYs were selected. Correspondingly, the number of grids filtered was 172. After these processes, the NARXnn was trained using the data from 172 grids of the 46-day final time series.

4.2. Training and Testing

The NARXnn was built for the two-year training period (2013–2014). After training of the NARXnn, the corresponding NARXnn SM during the training period was used to test the network performance. The training data (70% proportion) and the test data (30% proportion), which were described in Section 3.1, were used to evaluate the performance of the network, and the results are presented in

Table 2. Training performance of the NARXnn.

	Training Data	Test Data
Fused SM Max (m3/m3)	0.5152	0.4618
Corresponding NARXnn SM Max (m3/m3)	0.5095	0.4653
Fused SM Min (m3/m3)	0.0340	0.0360
Corresponding NARXnn SM Min (m3/m3)	0.0383	0.0383
Fused SM Mean (m3/m3)	0.1731	0.1733
Corresponding NARXnn SM Mean (m3/m3)	0.1730	0.1729
RMSE (m3/m3)	0.0098	0.0099
R	0.9934	0.9936

. All of the statistics showed a good relationship between the corresponding NARXnn SM and the fused SM.

Moreover, the scatter plots of the gridded values and the line profiles of all of the mean values of the grids between the corresponding NARXnn SM and the fused SM during the training period were constructed, as shown in Figure 5. The scatter diagram (Figure 5a) shows a “transverse line” near y = x (i.e., 1:1). The temporal plot (Figure 5b) of the corresponding NARXnn SM shows a typical “bell” shape with obvious seasonal features. In addition, the distribution of points (the 46 successive time series) revealed a uniform and homogeneous pattern, indicating that the filtering was appropriate.

4.3. Implement

After training, the NARXnn was implemented to retrieve SM for 2015. To generate successive time series of SM within the HRB, the input series (U) needed to be daily and continuous. To meet this requirement, a temporal interpolation process was executed to obtain a daily temporal time series of U. In other words, TB-X-V, TB-X-H, TB-Ka-V and LAI were adjusted to be daily and successive because H, PRC and DOY are already successive. After this, the NARXnn was implemented to generate SM.

5. Results and Validation

In this section, the NARXnn SM is compared with the in situ SM measurements, the original JAXA and LPRM SM products, and the GLDAS SM product. A comparison of the spatial distribution of the NARXnn SM and the GLDAS SM over the whole HRB is then presented.

5.1. Direct Validation

Generally, direct validation is performed using time series plots, scatter plots and performance summary tables. Time series plots can provide insight into the seasonal performance of the algorithm and the magnitude of its variance [4]. The latter two methods aim to analyze their R and RMSE values using in situ measurements. An SM product with a low RMSE value and a high R value is desirable. PRC data were also used in the analyses to detect the response of NARXnn SM to precipitation flags.

Time series plots (Figure 6) depict the dynamic range and seasonal features of the JAXA, LPRM, GLDAS and NARXnn SM.

In general, the NARXnn SM agreed well with the in situ data during the unfrozen season. A gradual rise occurred near 1 April and a decline occurred before 1 November. The NARXnn SM remained stable during the frozen season. JAXA was stable at all times. LPRM had a large dynamic range of 0.8 m³/m³ and always exhibited pronounced fluctuations. The GLDAS SM profiles still had a robust tendency to resemble the in situ data, if the exceptionally high values during the frozen season were ignored. In addition, during the unfrozen season, the GLDAS SM profiles seldom reached the height of the in situ SM. The NARXnn, LPRM and GLDAS SM products captured the PRC information with different amplitudes, whereas JAXA failed to capture the variation in PRC. As shown in

Table 3. Summary of soil moisture algorithm performance during the validation period (2015).

		Product	RMSE (m3/m3)	R	Bias	No. 1
The whole year	WATER-NET 1–4 (mean)	NARXnn	0.060	0.90	−0.02	365
		JAXA	0.153	0.62	−0.19	281.25
		LPRM	0.221	0.44	0.12	243
		GLDAS	0.100	0.36	−0.01	365
	HM-NET 1	NARXnn	0.043	0.85	0.01	365
		JAXA	0.076	0.35	−0.06	284
		LPRM	0.110	−0.50	0.06	276
		GLDAS	0.063	0.34	0.04	365
	HM-NET 2	NARXnn	0.034	0.87	0	365
		JAXA	0.091	0.40	−0.03	307
		LPRM	0.104	−0.73	−0.07	308
		GLDAS	0.074	0.13	−0.02	365
From May to September	WATER-NET 1–4 (mean)	NARXnn	0.063	0.50	0.02	153
		JAXA	0.206	0.24	−0.20	115.5
		LPRM	0.187	0.09	0.14	116.25
		GLDAS	0.086	0.37	−0.07	153
	HM-NET 1	NARXnn	0.048	0.60	0.04	153
		JAXA	0.097	0.32	−0.09	114
		LPRM	0.056	0.23	−0.01	117
		GLDAS	0.045	0.35	0.02	153
	HM-NET 2	NARXnn	0.032	0.58	−0.01	153
		JAXA	0.121	0.33	−0.12	127
		LPRM	0.116	0.20	−0.11	126
		GLDAS	0.092	0.42	−0.08	153

¹ No. denotes the number of days with recorded values.

, for WATER-NET 1–4 over the whole validation period (2015), the NARXnn SM was associated with a mean RMSE of 0.060 m³/m³, a mean Bias as low as −0.02, and a mean R value of 0.90. The comprehensive statistical data revealed that NARXnn SM showed substantial agreement with the in situ data. At HM-NET 1 and HM-NET 2, the values of R (R = 0.85 and 0.87, respectively) were slightly lower than those at WATER-NET 1–4, although a strong correlation with the in situ data is still evident. Compared with the other SM products, the NARXnn SM values were closer to the in situ data, as reflected by their lower RMSE and Bias values (RMSE = 0.043 and 0.034; Bias = 0.01 and 0, respectively). As a further sensitive test, the performance from May to September 2015 was analyzed. NARXnn SM had the highest R value (R = 0.50, 0.60 and 0.58 for WATER-NET 1–4, HM-NET 1 and HM-NET 2, respectively), revealing a relatively strong relationship. The values of RMSE and Bias were similar to the statistical figures obtained for the whole validation period (2015). Consequently, the NARXnn SM was in fairly good agreement with the in situ data compared with the other SM products.

Specifically, for WATER-NET 1–4, the NARXnn SM profiles showed a similar tendency as the in situ data during the unfrozen season, and they displayed a better sensitivity and response to PRC. Especially at WATER-NET 2 (Figure 6b), the NARXnn SM profile matched the in situ data remarkably well. JAXA generally maintained a constant value of approximately 0.1 m³/m³. LPRM showed the best response to PRC, but it extremely overestimated SM. GLDAS underestimated the SM values, but it had a similar amplitude of variation as the in situ data. For HM-NET 1 (Figure 6e), an analogous conclusion can be drawn, aside from the fact that the values obtained from all of the SM products were smaller. NARXnn SM performed the best of all of the SM products, and its profile was strongly sensitive to PRC. It should be noted that there was a relatively large gap between the unfrozen season and frozen season. Thus, the NARXnn SM profile did not show good agreement with the in situ data during the freezing/thawing season. The situation was different in the lower HRB. For HM-NET 2 (Figure 6f), the in situ SM profile remained almost constant most of the time, until external factors, such as rainfall and irrigation, came into play. We believe that this is, in part, due to the different topographic conditions in the flat, low-altitude, semiarid area of the lower HRB. Against this background, during the unfrozen season, the NARXnn SM profile had a shape that generally resembled that of the in situ data. The GLDAS SM profile paralleled the data but underestimated the in situ SM values substantially.

As a final reminder, the JAXA and LPRM SM products contained differing degrees of missing data in their time series (Table 3); the NARXnn was able to fix this problem. GLDAS failed to estimate the unfrozen water content during the frozen season.

Scatter plots (Figure 7) and performance summaries (Table 3) provide insight into the accuracy and correlation of each SM product during the entire validation period. As shown in Figure 7, the range of in situ SM data was 0–0.5 m³/m³. The JAXA values occurred in a line parallel to the X-axis and clustered in the lower left corner, indicating that JAXA had a restricted range and underestimated SM with a large bias. In contrast to JAXA, the point cloud of LPRM resembled a line parallel to the Y-axis. Given that LPRM SM fluctuated “globally”, it was unable to detect SM variations as a result of its huge range. To conclude, both JAXA and LPRM failed to capture the variations in SM conditions or account for the dynamic changes in SM. Moreover, the JAXA/LPRM SM underestimated/overestimated in situ measurements and displayed high negative/positive bias, respectively. The points of GLDAS SM were seldom within the ±20% error lines. The NARXnn SM values represented the closest approximation to the 1:1 line when plotted against the in situ data.

In terms of the actual SM difference between WATER-NET 1–4 and HM-NET 1–2, we believe that this difference is, in part, due to the differences in the hydrological, meteorological and topographic characteristics of the upper, middle and lower portions of the HRB. Specifically, WATERNET 1–4 were located in mountainous and high-altitude areas composed of plateaus and basins (Figure 1), which increase the water-holding capacities of these sites. By contrast, HM-NET 2 was located on the plains (Figure 1), which are covered by barren areas and sparse vegetation. These large distinctions probably cause differences in the surface SM values over the long run. In addition, during the validation period, the PRC falling within HM-NET 2 was less than one-fourth of the WATER-NET 1–4 values, which only exacerbated the error.

5.2. Indirect Validation

Indirect validation was performed using a method that involved comparing the NARXnn SM and the GLDAS SM with respect to their spatial distributions on selected typical days (the 4th day of each month in 2015) during the validation period.

The spatial distribution of the HRB SM maps (Figure 8 and Figure 9) provided spatial information on the differences between the NARXnn SM and GLDAS SM on the 4th day of each month in 2015. During the unfrozen season, the NARXnn SM showed an obvious spatial distribution according to GLDAS SM, despite the differences in their codomains. In the northeastern portion of the NARXnn SM maps, the lower course of the Heihe River was easily detected. In the northwestern portion of the NARXnn SM maps, the positions of the relatively humid and arid areas were roughly consistent with those shown on the GLDAS SM maps. The complexity of the spatial distribution of SM in the southern mountainous area was also apparent. During the frozen season, the SM differences exhibited relatively large values because the GLDAS SM were inaccurate.

As has been shown above, the NARXnn SM showed a distribution and seasonal change similar to those of the in situ SM data within the WATER-NET and HM-NET areas. Therefore, we have reason to believe that the NARXnn SM maps provide a relatively accurate representation of the actual distribution of SM.

6. Conclusions

As stated in Section 1, the goal of this study was to obtain a spatially and temporally continuous, high-quality, time series SM product. Therefore, the estimation of time series SM using NARXnn over the HRB in China was explored. The output variable was the fused SM, and the input variables were X-band dual polarization TB, as well as LAI, Ka-band V polarization TB, PRC, DEM and DOY. Only daytime data at 13:30 were used because of their robust performance in previous studies. The NARXnn was then trained using data covering two years (1 January 2013–31 December 2014) with the trainbr training function and the learngdm weight/bias learning function. The NARXnn was well trained, and the R values for both the training and test data sets were greater than 0.99.

To obtain more accurate input sequences, additional steps were taken. Fused SM was calculated from the AMSR2 JAXA and LPRM SM products using mathematical statistical methods. Before this process, a long-term time series analysis was performed to characterize the seasonal characteristics of the JAXA and LPRM SM products. The unfrozen, frozen and freezing/thawing seasons were considered separately, and different fit coefficients were obtained.

Traditional direct validation (using time series plots and scatter plots) and additional indirect validation (employing spatial distribution maps) during the validation period (2015) were performed. Among all of the SM products tested (i.e., NARXnn, JAXA, LPRM and GLDAS SM), the NARXnn SM product performed the best; the NARXnn SM had the highest R value and the lowest values of RMSE and Bias compared with the in situ data, as well as the best response to PRC. The NARXnn SM profile agreed better with the in situ SM profile than the other SM data products did. Moreover, the NARXnn SM had a spatial distribution that was similar to that of the GLDAS SM, and the Heihe River was identifiable from them.

Although time series SM data covering only one year were retrieved and validated by NARXnn due to the limitations of remote sensing observations and in situ data, this algorithm can be used to obtain longer, spatially and temporally continuous, high-quality SM products with sufficient resources.

The following are major shortcomings of this study, which we hope will be addressed in the future:

• The NARXnn has a seven-dimensional input vector (i.e., TB-X-H, TB-X-V, LAI, TB-Ka-V, PRC, DEM and DOY) and a one-dimensional output vector (SM). Efforts must focus on connecting more correlated variables (e.g., soil texture and land surface roughness) with SM.
• The time series SM predicted by NARXnn was based on a priori information, which likely leads to inaccurate estimation during the freezing/thawing seasons.
• This retrieval method is now being implemented to generate time series SM over the HRB in China from AMSR2 data. It could be applied to other areas and even the globe. In addition, it could also be applied to satellite data from other passive microwave sensors, such as SMAP and SMOS, in the future.

Moreover, although the SM values estimated by NARXnn cannot be considered to depict a real SM distribution due to the absence of inadequate in situ measurements within the HRB, it is of vital importance to note that the NARXnn SM product is capable of monitoring the variations and dynamic changes in SM at different latitudes and under various types of land cover.

To conclude, spatially and temporally continuous SM products with high precision covering the HRB were obtained at the scale of the HRB. With long-term remote sensing observations and other auxiliary data, the NARXnn algorithm can be used to produce long time series of spatially and temporally continuous SM products that have high precision, thereby benefiting future research and applications in hydrology, meteorology, agriculture and forestry. In addition, more effort should be put into utilizing additional in situ measurements and evaluating the performance of the NARXnn algorithm.