Merging Microwave, Optical, and Reanalysis Data for 1 Km Daily Soil Moisture by Triple Collocation

Abstract

High-spatiotemporal resolution soil moisture (SM) plays an essential role in optimized irrigation, agricultural droughts, and hydrometeorological model simulations. However, producing high-spatiotemporal seamless soil moisture products is challenging due to the inability of optical bands to penetrate clouds and the coarse spatiotemporal resolution of microwave and reanalysis products. To address these issues, this study proposed a framework for multi-source data merging based on the triple collocation (TC) method with an explicit physical mechanism, which was dedicated to generating seamless 1 km daily soil moisture products. Current merging techniques based on the TC method often lack seamless daily optical data input. To remedy this deficiency, our study performed a spatiotemporal reconstruction on MODIS LST and NDVI, and retrieved seamless daily optical soil moisture products. Then, the optical-derived sm₁, microwave-retrieved sm₂ (ESA CCI combined), and reanalysis sm₃ (CLDAS) were matched by the cumulative distribution function (CDF) method to eliminate bias, and their weights were determined by the TC method. Finally, the least squares algorithm and the significance judgment were adopted to complete the merging. Although the CLDAS soil moisture presented anomalies over several stations, our proposed method can detect and reduce this impact by minimizing its weight, which shows the robustness of the method. This framework was implemented in the Naqu region, and the results showed that the merged products captured the temporal variability of the SM and depicted spatial information in detail; the validation with the in situ measurement obtained an average ubRMSE of 0.046 m³/m³. Additionally, this framework is transferrable to any area with measured sites for better agricultural and hydrological applications.

1. Introduction

Soil moisture (SM) is a crucial part of the water [1] and carbon cycle [2]. It has a considerable impact on the land–atmosphere interactions as one of the primary input factors for various land surface and hydrological models [3,4]. To monitor vegetation growth [5,6], weather [7], and extreme climates, such as drought and flood [8,9,10], high spatial resolution and spatiotemporally continuous SM products are essential.

The SM information is mainly acquired by in situ measurement, remote sensing retrieval, and model assimilation. Among them, the in situ measurement with the highest accuracy can be regarded as the “true value”, but is limited by high cost and weak spatial representation [11,12]. In addition, the change in SM affects vegetation growth, temperature, and evapotranspiration, which can be detected by optical remote sensing. Therefore, various indices including vegetation index, temperature index, and thermal inertia index were developed [13,14,15] to indirectly estimate SM. Optical SM products offer a high spatiotemporal resolution. For vegetated land, optical signals are only affected by vegetation, which allows for calibrating the vegetation effects on SM retrieval. However, clouds and rain contaminated the optical remote sensing signals, preventing their applications in monitoring soil and vegetation conditions.

Compared with optical remote sensing, microwave SM is based on the physical process of electromagnetic wave propagation and interaction with the soil [16,17]. Due to the stark contrast between the dielectric constants of dry soil (3) and liquid water (80), more moisture content in the soil mixture presents a higher dielectric constant. The conversion between SM and the dielectric constant is conducted by dielectric models. Microwave remote sensing can be categorized into active and passive modes. The active microwave SM retrievals include change detection algorithms [18] and iterative optimization [19]. Passive microwave SM retrieval has the following three categories: i) Reverse-order solutions based on the radiative transfer model, such as the single channel algorithm (SCA) [20] for soil moisture active passive (SMAP) mission. ii) Iterative methods based on radiative transfer model, such as multi-angle algorithm [21] for soil moisture and ocean salinity (SMOS) mission, dual-channel (DCA) [22], multi-temporal dual-channel (MT-DCA) [23], and multi-channel collaborative algorithm (MCCA) [24] for SMAP, etc. iii) Microwave remote sensing index-based methods, such as the land parameter retrieval model (LPRM) algorithm [25] in Advanced Microwave Scanning Radiometer (AMSR-E). Although passive microwaves can obtain relatively accurate SM, they are constrained by coarse spatial resolution, which is often tens of kilometers. In addition, the microwave propagates within vegetation layers with multiple paths, thus both the soil and vegetation contribute to the received signals. It is challenging to completely isolate the soil and vegetation components, resulting in poor performances over high vegetation conditions [26,27].

Alternatively, hydrological or land surface process models driven by atmospheric data can simulate SM. By assimilating remote sensing, station data into these process models, reanalysis SM products with a daily or even sub-daily temporal resolution and continuous spatial coverage, were obtained. Furthermore, reanalysis data provide the root zone SM, which is more critical than surface SM for vegetation growth, but it was limited by coarse spatial resolution and insufficient precision.

Considering the strengths and weaknesses of each product, the uncertainty can be reduced by multi-source data fusion to get an optimal estimation. Nevertheless, it is a vital issue to determine the optimal weights for each data in the fusion process. The Kalman filter is one of the most common methods for assigning weights [28], but most studies adopted prior knowledge to define the relative error of each product. Therefore, the relative weights are subjective, and the results are not necessarily optimal. In recent years, some studies were carried out on data fusion using the triple collocation (TC) method, which calculates the relative errors of three datasets that are independent and linearly connected to the true value. For example, the TC method was used to merge precipitation data [29] and estimate terrestrial water storage [30]. In terms of the SM data merging, Zeng et al. investigated the TC approach to combine ERA-Interim reanalysis data, microwave satellite data, and in situ measurement data [31]. Peng et al. merged microwave, JULES model simulation, and site data into triplets to obtain fine spatial resolution SM [32]. Moreover, the well-known ESA CCI SM products are based on the TC method, which merged active microwave (ERS1/2 SCAT, ASCAT-A/B/C) and passive microwave (SMMR, SSM/I, TMI, AMSRE, WindSat, SMOS, FY-3B/C/D MWRI, GPM, AMSR2, and SMAP) products to obtain the long-term (1978-2021) SM datasets [33,34].

However, the current TC merging of SM is mainly conducted with microwave, model simulation, and site data. Compared with these data, optical data have the advantage of high spatial resolution, but only a few studies introduced optical data under clear-sky conditions into triples [35] due to the issue of pixel contamination by clouds and rain. To solve this problem, our study proposed a temperature difference-vegetation index to retrieve SM from MODIS vegetation and land surface temperature (LST) data after applying spatiotemporal interpolation [36]. Then, using the ESA CCI SM products as a reference, the CDF matching approach was performed to match optical and reanalysis SM products. The matched three data were formed into triples, and their relative errors were estimated by the TC method. Finally, the weight of each product was calculated by the least squares algorithm, resulting in a merged SM product after a significant test. This study combined the strengths of optical, microwave, and reanalysis data to obtain seamless daily 1 km SM products with satisfactory accuracy.

2. Study Area and Materials

2.1. Ground Measurement Data

The study region is in Naqu (Figure 1), which is in the heart of the Tibetan Plateau, China. The Naqu SM observation network (CTP-Naqu) covers the longitude range of 91.5 to 92.5°E and the latitude range of 31 to 32°N [37]. At a resolution of 1 km, the entire study area includes 10,000 pixels. This study used the SM data at a depth of 0–5 cm over all stations during 2017–2020.

The land cover type in CTP-Naqu is mainly sparse alpine grassland, and the climate type is cold and semi-arid [38]. The soil freezes in November and begins to thaw in May of the following year due to the cold climate [39]. Soil freezing will affect the measurement of SM. Therefore, the data from May to October of each year were chosen for this study.

2.2. Data for Merging

2.2.1. ESA CCI SM Products

The European Space Agency (ESA) provides the SM products of the Climate Change Initiative (CCI) project, which consists of active, passive, and combined products, with a spatial resolution of 25 km and a period of more than 40 years (1978–2021). ESA CCI SM team updates the algorithm every year and has to date released 13 versions of the product [33,40]. The latest version (V 07.1) SM was selected in the study, which improves spatiotemporal coverage compared with the previous version. In addition, previous studies have shown that the combined product has the best performance in the Naqu region [41,42] among the three products, thus we chose the combined product.

2.2.2. CLDAS SM Product

The reanalysis SM product was derived from the China Land surface Data Assimilation System (CLDAS-V 2.0 near real-time product, 0.0625 × 0.0625°) by China Meteorological Administration. Based on the data fusion and assimilation technology, CLDAS leverages multiple forcing variables to drive the various land surface models (e.g., CLM and Noah-MP) to obtain SM [43,44]. In addition, CLDAS SM data consider five soil layers to capture the vertical profiles. To match the sensed depth of the satellite, the daily data at a depth of 0–5 cm were used in this study.

2.2.3. MODIS Product

MODIS sensors are carried on the Aqua and the Terra spacecraft. For the LST, this study adopted the MYD11A1 product from the Aqua spacecraft with local overpass times at 1:30 pm and 1:30 am, respectively. The temperature difference between the two overpasses was considered as a proxy for the maximum temperature difference on that day. In addition, an 8-day normalized differential vegetation index (NDVI) product was produced by combining two 16-day NDVI data from Aqua (MYD13A2) and Terra (MOD13A2) at 1 km resolution. Moreover, we used the MODIS land cover type product (MCD12Q1) and selected the IGBP classification layer to mask pixels with water and participate in the error calculation.

2.2.4. DEM Data

The Space Shuttle Radar Topographic Mapping Mission (SRTM) [45] provided the 3 arc-second digital elevation model (DEM) data that were used in the LST interpolation.

3. Methods

As shown in Figure 2, all data were preprocessed and resampled to 1 km at first. Soil moisture retrieval from optical data was carried out as follows: i) Daily NDVI products were generated by Savitzky-Golay (SG) filter algorithm [46] and linear interpolation. ii) The multi-temporal interpolation method was utilized to create the daily seamless LST product using daily NDVI, DEM, and MODIS LST data. iii) Temperature difference-vegetation index (TDVI) was proposed to retrieve SM from the interpolated NDVI and LST. Subsequently, we combined optical (TDVI-based), microwave (ESA CCI combined), and reanalysis (CLDAS) soil moisture products into triples. The deviation between the three products was corrected using the CDF matching algorithm, and the relative error of each product was estimated by the TC approach. Finally, the merged daily SM products at 1 km resolution were obtained based on the least square merging and relevant significance test, followed by an evaluation against station observations.

3.1. Retrieve Cloud-Free Daily Optical SM

3.1.1. Obtain Daily MODIS NDVI

This study used 16-day MODIS NDVI (MYD13A2 and MOD13A2) rather than daily NDVI. This is due to the fact that the 16-day NDVI was synthesized by selecting high-quality pixels, which considerably decreased the number of pixels that were missing or of poor quality due to cloud contamination. The MYD13A2 and MOD13A2 data were combined to obtain 8-day synthetic NDVI data.

The steps involved in interpolating NDVI data include the following. First, the pixel_reliability layer was utilized to control the quality, and then the SG filter method was conducted to reduce the noise in the NDVI image. The SG filter is a polynomial smoothing algorithm based on the least squares principle for smoothing time series data. Finally, pixel-by-pixel linear interpolation is performed on the time series to acquire daily NDVI data [47].

3.1.2. Develop Seamless MODIS LST

This study adopted a multi-temporal interpolation method to interpolate LST, which assumes that the LST variation characteristics in different pixels under similar environmental conditions (including seasons, altitudes, and vegetation) are highly correlated [36,48,49]. For two LST images acquired with a close date interval ($d_0$ and $d_1$), the following relationship exists:

$$LS{T}_{d_1}=f\left(LS{T}_{d_0},NDV{I}_{d_1},DEM\right)$$

(1)

where $NDV{I}_{d_1}$ represents the NDVI image corresponding to the $d_1$ day, and DEM is the digital elevation of the study area. With a linear assumption, the relationship $f$ becomes:

$$LS{T}_{d_1}=a·LS{T}_{d_0}+b·DEM+c·NDV{I}_{d_1}+d$$

(2)

where a, b, c, and d represent regression coefficients.

Based on the above principles, quality control was performed on all LST images by the QC layer at first, then, the images with more than 90% clear-sky LST pixels were selected as the reference images ($LS{T}_{d_0}$), and the other images were interpolated images ($LS{T}_{d_1}$). In addition, the date interval between $d_0$ and $d_1$ should be within 30 days to avoid the seasonal variation of LST. If the date interval is greater than 30 days, we searched for reference images of adjacent years. The specific reference image information is shown in Tables S1 and S2 in the Supplementary Materials.

Next, clear-sky pixels were selected in the interpolated image, and the LSTs in the reference image corresponding to these clear-sky pixels were found. Subsequently, the regression coefficients in Equation (2) are fitted using the ancillary remote sensing dataset (NDVI, DEM) of the corresponding location. Finally, based on Equation (2) with known regression coefficients, the LST estimation of missing pixels in the interpolated images was completed.

Through the above process, the percentage of valid LST pixels for each image was above 90%. To produce a seamless LST image, the remaining pixels were interpolated using the inverse distance weight (IDW) technique.

3.1.3. Generate Cloud-Free Daily SM Product

In optical remote sensing, many studies have used surface temperature, vegetation index, and albedo as indicators of SM, such as the temperature vegetation dryness index (TVDI) [14], apparent thermal inertia (ATI) [50], and vegetation supply water index (VSWI) [51]. Inspired by the above index, we proposed a temperature difference-vegetation index (TDVI), defined as the ratio between diurnal temperature difference ($LS{T}_{day}-LS{T}_{night}$) and NDVI:

$$TDVI=\frac{LS{T}_{day}-LS{T}_{night}}{NDVI}$$

(3)

This index and the measured SM have a strong correlation. Therefore, cloud-free daily TDVI data can be calculated using cloud-free daily LST and NDVI data, and then the cloud-free daily SM product over the study area was retrieved by the fitted relationship between TDVI and the measured SM.

3.2. Multi-Source SM Products Merging

To merge multi-source SM products, all data were adjusted to a unified numerical range through bias correction, and the relative error of each data was estimated by the TC method. Then, the weight of each product was calculated by the least squares technique according to the relative error. Finally, the merging of multi-source data was completed based on the significance test.

3.2.1. Bias Correction using CDF Matching

Due to discrepancies in sensor frequencies, observation principles, and retrieval algorithms, different SM products have certain systematic biases that need to be removed before data merging. This study adopted the CDF matching technique, which is considered as an enhanced non-linear method for eliminating the statistical moment difference between two datasets [52,53], and does not change the trend of the data. Additionally, the reference data quality influences the accuracy in merging the product [54]. We used ESA CCI combined SM data as the reference data since previous studies [55,56] demonstrated its excellent accuracy. The CLDAS and TDVI-based SM products were normalized to have the same range and distribution as the ESA CCI combined SM by CDF matching.

3.2.2. Error Estimation by Triple Collocation

Without relying on the true value, the TC technique [57] possesses the ability to calculate the error of each product, which requires that the errors of the three data are uncorrelated. We selected optical, microwave, and CLDAS data to build triples. The data sources and algorithms are different, thus the precondition of error irrelevance is satisfied.

The TC method assumes a linear correlation between the SM product and the true value:

$$\left\{\begin{array}{c}S{M}_x={\alpha }_x+{\mathsf{\beta }}_x\theta +{\epsilon }_x\\ S{M}_y={\alpha }_y+{\mathsf{\beta }}_y\theta +{\epsilon }_y\\ S{M}_z={\alpha }_z+{\mathsf{\beta }}_z\theta +{\epsilon }_z\end{array}\right.$$

(4)

where $\theta$ is the true SM; $S{M}_i$ ($i\in \left\{x,y,z\right\}$) represent three SM products matched by the CDF method; the additive and multiplicative deviations from the true value are represented by $\alpha$ and $\mathsf{\beta }$, respectively; and $\epsilon$ is the random error.

Since the errors between products and true values are independent, according to the covariance calculation formula, the following equation can be obtained:

$${\mathrm{C}}_{ij}=Cov\left(S{M}_i,S{M}_j\right)=\left\{\begin{array}{c}{\beta }_i{\beta }_j{\sigma }_{\theta }^2, for i\ne j\\ {\beta }_i^2{\sigma }_{\theta }^2+{\sigma }_{{\epsilon }_i}^2, for i=j\end{array}\right.$$

(5)

where ${\sigma }^2$ is the random error variance. Let ${\mathsf{\lambda }}_i={\beta }_i{\sigma }_{\theta }$, the above equation transforms into:

$${\mathrm{C}}_{ij}=\left\{\begin{array}{c}{\mathsf{\lambda }}_i{\mathsf{\lambda }}_j, for i\ne j\\ {\mathsf{\lambda }}_i^2+{\sigma }_{{\epsilon }_i}^2, for i=j\end{array}\right.$$

(6)

By multiplying and subtracting Equation (6) in pairs, the error variance value of each data can be calculated as:

$$\left\{\begin{array}{c}{\sigma }_{{\epsilon }_x}^2={\sigma }_x^2-\frac{{\sigma }_{xz}{\sigma }_{xy}}{{\sigma }_{yz}}\\ {\sigma }_{{\epsilon }_y}^2={\sigma }_y^2-\frac{{\sigma }_{yz}{\sigma }_{xy}}{{\sigma }_{xz}}\\ {\sigma }_{{\epsilon }_z}^2={\sigma }_z^2-\frac{{\sigma }_{xz}{\sigma }_{yz}}{{\sigma }_{xy}}\end{array}\right.$$

(7)

Note that only pixels with more than 100 samples are retained in TC error calculation, as the results are unreliable when the number is less than 100 [58]. Furthermore, when the computed error variance ${\sigma }_{{\epsilon }_i}^2$ of any product is negative, the error variances of other products in the triplet are invalid. To fill the invalid pixels generated by the above two cases, we averaged all error variances of valid pixels under the same land cover type, and then assigned this average value to the invalid pixels for the same cover type.

3.2.3. Weight Estimation via Least Square Merging

A least square merging technique introduced by Yilmaz et al. was successfully adopted to merge thermal infrared, passive microwave, and model SM products [59]. It is advantageous to combine various SM products since it minimizes random errors. The least squares merging form is as follows:

$$S{M}_{merge}=w_{x}S{M}_x+w_{y}S{M}_y+w_{z}S{M}_z$$

(8)

where $S{M}_{merge}$ is merged SM and $w_i$ is the weight of each SM product. To obtain an unbiased merged estimate, the sum of $w_i$ is equal to 1.

The purpose of the least square method is to minimize the error variance of the merged SM product (${\sigma }_{{\epsilon }_m}^2$), which has a cost function (J) form:

$$J={\sigma }_{{\epsilon }_m}^2=w_x^2{\sigma }_{{\epsilon }_x}^2+w_y^2{\sigma }_{{\epsilon }_y}^2+w_z^2{\sigma }_{{\epsilon }_z}^2$$

(9)

$$J={\sigma }_{{\epsilon }_m}^2=w_x^2{\sigma }_{{\epsilon }_x}^2+w_y^2{\sigma }_{{\epsilon }_y}^2+{\left(1-w_x-w_y\right)}^2{\sigma }_{{\epsilon }_z}^2$$

(10)

Setting $\frac{\partial J}{\partial x}=0$ and $\frac{\partial J}{\partial y}=0$, $w_i$ can be obtained:

$$\left\{\begin{array}{c}{w}_x=\frac{{\sigma }_{{\epsilon }_y}^2{\sigma }_{{\epsilon }_z}^2}{{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_y}^2+{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_z}^2+{\sigma }_{{\epsilon }_y}^2{\sigma }_{{\epsilon }_z}^2}\\ w_y=\frac{{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_z}^2}{{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_y}^2+{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_z}^2+{\sigma }_{{\epsilon }_y}^2{\sigma }_{{\epsilon }_z}^2}\\ w_z=\frac{{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_y}^2}{{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_y}^2+{\sigma }_{{\epsilon }_x}^2{\sigma }_{{\epsilon }_z}^2+{\sigma }_{{\epsilon }_y}^2{\sigma }_{{\epsilon }_z}^2}\end{array}\right.$$

(11)

Equation (11) indicates that the weight of one product is proportional to the uncertainty of the other two products. When only two products are available, the formula for $w_i$ is:

$$\left\{\begin{array}{c}{w}_x=\frac{{\sigma }_{{\epsilon }_y}^2}{{\sigma }_{{\epsilon }_x}^2+{\sigma }_{{\epsilon }_y}^2}\\ w_y=\frac{{\sigma }_{{\epsilon }_x}^2}{{\sigma }_{{\epsilon }_x}^2+{\sigma }_{{\epsilon }_y}^2}\end{array}\right.$$

(12)

Gruber et al. [60] proposed the defining of weights according to the signal-to-noise (SNR) ratio attributes of dataset rather than error variances. However, all the data were unified into the same data space by CDF matching (Section 3.2.1), thus Equations (11) and (12) are feasible.

3.2.4. Merging Based on Correlation Significance Level

The reliability of TC is weak when one or more of the products in the triplet have low temporal coverage or poor quality. Therefore, if a product in a triple is not strongly linked with the other products, TC results are frequently regarded as inaccurate. A scheme based on a correlation significance level (p-value) of 0.05 [32,60] was adopted to increase the spatial coverage of the merged SM products.

As shown in

Table 1. Merging scheme based on Pearson correlation with a significance level of 0.05.

If (p-Value < 0.05)?			Decision
C–E *	C–T	E–T
✓	✓	✓	TC weighted average
✓	✓	✗	C
✗	✓	✓	T
✓	✗	✓	E
✓	✗	✗	Arithmetic mean (C,E)
✗	✓	✗	Arithmetic mean (C,T)
✗	✗	✓	Arithmetic mean (E,T)
✗	✗	✗	Disregard pixel (NaN)

* C represents CLDAS SM data, T refers to TDVI-based SM data, and E is ESA CCI combined SM data.

, the merging scheme on each pixel contains four cases. In the first case, the TC weighted average approach (Section 3.2.3) was applied when all three products are significantly correlated (p-value < 0.05). In the second case, when one product in the triplet is closely related to the other two products, we chose this product. In the third case, the correlation significance level is present in only two products. We calculated the arithmetic mean of these two products. Finally, the pixel was ignored if there is no strong association between the three products.

3.3. Spatiotemporal Analysis Method

To investigate how the four SM products (CLDAS, ESA CCI combined, TDVI-based, and merged SM) change in temporal and spatial dimensions over the study period (every May to October from 2017 to 2020), we employed the Hovmöller diagram [61], which averages all the values in a longitude or latitude row, places the average values on one axis, and the other axis represents time. The longitude or latitude axis reveals the spatial distribution of SM, while the time axis shows the seasonal and inter-annual variation of SM.

4. Results and Discussions

4.1. TDVI-Based SM Retrieval

Since the land cover is quite homogeneous (sparse alpine grassland) in this study, building a station-by-station relationship is not necessary. We averaged the SM of all stations and the corresponding TDVI values across the study region, and then fitted the relationship between them. As seen in Figure 3, TDVI and the measured SM exhibit a negative correlation. The mechanism between TDVI and SM is explained as follows:

First, the amount of soil moisture influences the soil thermal inertia since water has a large specific heat capacity. Additionally, soil thermal inertia reflects the ability of soil to resist temperature changes; the more soil moisture, the higher the thermal inertia, and the lower the soil temperature. Therefore, the soil temperature difference and SM within a day are negatively correlated.

Moreover, SM is closely related to vegetation growth. When the vegetation coverage is high, it reduces the evaporation of surface SM, and the vegetation index reflects the growth of vegetation; the more luxuriant the vegetation, the higher the index, thus SM and vegetation index have a positive association. Eventually, the TDVI was obtained by dividing the soil diurnal temperature difference by the vegetation index, which can better describe the variation of soil moisture.

A logarithmic relationship exists between site-averaged SM and pixel-averaged TDVI. Utilizing the logarithmic equation obtained by fitting, the SM value for each pixel could be acquired from the TDVI value.

4.2. Triple Collocation Analysis

4.2.1. Effect of LST Interpolation on Triple Collocation

Figure 4 illustrates that despite the 4-year time period of this study (2017–2020), only 39.68% of the 10,000 pixels in the study area had triples with days greater than 100, and the maximum number of days is only 138. If the LST is not interpolated, more than half of the weights derived by the TC approach are not available, and the resulting merged products will have a considerable number of missing pixels. The triplet ratio is increased to 99.93% after interpolation, and most pixels are available for more than 300 days.

The ratio of available triples before and after interpolation has significantly increased. Given this, LST interpolation efficiently addresses the issue of a severe lack of optical data pixels and enables the optical SM data under all weather conditions to participate in the TC analysis.

4.2.2. Error and Weight Analysis

Figure 5 exhibits the errors and weights of CLDAS, ESA CCI combined, and TDVI-based SM products. Due to the CDF matching operation, the error here is the relative error between each SM data. Among the three products, ESA CCI combined SM has the smallest overall error, which supports our decision to use it as the reference data for CDF matching. The error of TDVI-based SM product ranks second, which proves that the SM retrieved through our developed TDVI index is reliable. The CLDAS error is the largest, and we also discovered that the CLDAS product contains unreasonable values on some stations, which may be caused by inappropriate parameter settings (detailed analysis in Section 4.3.1).

The weight of each pixel for each product was determined based on the results of the error. The ESA product had an average weight of 0.418, which had a greater weight in the central and marginal regions. The TVDI product had a superior weight in the rest of the region, with an average weight of 0.378. Whereas, the CLDAS product was assigned an average weight of 0.204, which had a low weight throughout the study area.

4.3. Merged Soil Moisture Results

4.3.1. Spatiotemporal Variation of SM

The Hovmöller diagrams in Figure 6 and Figure 7 demonstrate that ESA CCI combined SM product is heavily missing in May and October (dark blue), while the merged products have the highest spatiotemporal coverage. In addition, we observed that the spatiotemporal distribution of CLDAS data is entirely different from the other three data throughout 2017 and from May to June in 2018. To further explore this issue, time series plots were produced between four products and in situ SM data, selecting two of the sites as examples (Figure 8), and the time series plots for all sites are shown in the Supplementary Materials. In Figure 8, the CLDAS data deviate from other data in the above period, which further confirms our findings through the Hovmöller diagram. In addition, the CLDAS soil moistures were lower than the other products at site BC03 in 2019–2020, which may be due to the high soil organic matter content (SOC) of the top soil in the Naqu region; therefore, the porosity and water retention capacity are large, resulting in high soil water content [44,62]. The surface process models driven by CLDAS products may not consider the effect of SOC and therefore underestimated the soil moisture. Nevertheless, the CLDAS data also exhibited a specific cumulative trend, thus we suggested that the probability density matching step needs more verification in the Naqu region.

It should be emphasized that the anomaly of CLDAS data had no impact on the study results since the TC technique was used to calculate the relative error of each product. Section 4.2 shows that the CLDAS error obtained by the TC method was the largest and it was assigned a low weight. This phenomenon also proves the robustness of our proposed framework.

The spatiotemporal distributions of four SM data followed similar patterns except for the anomalous part. In terms of time variation, SM in each year was relatively low in May, reached a high level in July and August, then began to decline, and reached a low level in October. This temporal characteristic is consistent with the time pattern of soil freezing and thawing on the Tibetan Plateau [63]. Regarding the spatial distribution, SM varied slightly in longitude or latitude, partly since the averaging process may reduce the spatial variability. Additionally, there was insufficient spatial heterogeneity since the study area was small (1° × 1°) and had a uniform land cover.

In summary, the merged product captures the temporal variability of SM well, greatly improves the spatial resolution, and expands spatial and temporal coverage.

4.3.2. Compared with the In Situ Data

This study evaluated four products using all of the available measured sites during the study period and selected ubRMSE, RMSE, Correlation, and Bias as the performance metrics [64]. After removing the stations with no observation data and the total number of observations less than 50, 48 stations were finally involved in the calculation. The boxplots in Figure 9 display the metric statistics of the four products. Overall, the ESA CCI combined SM performs the best, followed by the merged product, TDVI-based SM is the third, and the CLDAS product performs the poorest. Although the merged results rank second, the difference with ESA CCI is minimal. The average ubRMSE of the merged SM is close to 0.04 m³/m³ (0.046 m³/m³), achieving satisfactory accuracy. More importantly, the merged product has a higher spatiotemporal resolution compared with the other products.

Furthermore, Figure 10 shows the metrics distribution of the merged SM product at each site. The product performed better in the northwest part with dense sites and worse in the northeast part with sparse sites. The reason may be that dense sites will have a greater impact on the TDVI fitting results. The effect of the number of in situ stations on the merged results will be further investigated in a subsequent study.

4.4. Uncertainty Analysis

The uncertainties in this study are primarily attributed to the following three aspects:

First, when retrieving soil moisture from optical data, there may be an oversaturation in NDVI, and the effect of vegetation indices, such as enhanced vegetation index (EVI) and leaf area index (LAI) should be explored in future studies. In addition, the impact of the number of measured sites on the retrieval results should be further considered.

Second, in estimating the product error by the TC method, this study estimated the time-invariant error, but the error of SM product varies with vegetation phenology, surface roughness, and environmental conditions. Therefore, the time-variant error should be calculated to combine soil moisture products [65,66]. In addition, the CDF matching method was used in this study to eliminate the bias between products. Although this method is commonly used, different rescaling techniques will affect the accuracy of TC error estimation, thus the effectiveness of methods, such as variance and mean matching and normalization can be further investigated [67].

Finally, in product merging, the nearest neighbor method was adopted to resample all products to 1 km; however, the original resolutions of ESA CCI and CLDAS are 25 and 6.25 km, respectively. There is a scale difference between them and the 1 km optical products, which will lead to uncertainty in the results.

5. Conclusions

This paper proposed a framework to generate seamless 1 km daily soil moisture products by the spatiotemporal interpolation technique and the triple collocation method. This framework merged optical, microwave, and reanalysis data to leverage their strengths.

First, the interpolation method was used to improve the spatiotemporal coverage of optical data, which tackled the pain point of a serious lack of optical data and realized the operation of introducing all-sky optical data into TC triples, which is rare in previous TC-based merging studies.

Second, a temperature difference-vegetation index was established, which has a close relationship with the measured SM. This index is a reference for soil moisture retrieval in the Naqu area. It should be noted that this index serves as a proxy for soil moisture in this framework, and the soil moisture index can be adjusted in other regions; therefore, our method is flexible and transferable.

In addition, the merged framework based on the TC method can not only obtain the merged soil moisture products, but also evaluate the performance of input products. Through the study, the CLDAS products in the Naqu area were found to exhibit anomalies in 2017 and the first half of 2018, and the causes of the anomalies were further explored through site validation. The CLDAS product needs further validation to check its accuracy in the Tibetan Plateau.

Finally, spatiotemporal analysis and metrics evaluation revealed that the merged products captured the dynamic changes in soil moisture well, depicted spatial details better, and achieved satisfactory accuracy.

In conclusion, this research produced a set of high spatiotemporal resolution soil moisture data that can be utilized as an input variable in atmospheric or land surface models to facilitate climate and vegetation analysis in the Tibetan Plateau as well as assist in regional water resource management. Furthermore, the ESA CCI SM and MODIS (or its successor VIIRS) products have a long time series, allowing us to extend the merged products to longer periods in future studies.