Triple Collocation Analysis for Two Error-Correlated Datasets: Application to L-Band Brightness Temperatures over Land

: The error characterization of satellite observations is crucial for blending observations from multiple platforms into a unique dataset and for assimilating them into numerical weather prediction models. In the last years, the triple collocation (TC) technique has been widely used to assess the quality of many geophysical variables acquired with different instruments and at different scales. This paper presents a new formulation of the triple collocation (Correlated Triple Collocation (CTC)) for the case of three datasets that resolve similar spatial scales, with two of them being error-correlated datasets. Besides, the formulation is designed to ensure fast convergence of the error estimators. This approach is of special interest in cases such that ﬁnding more than three datasets with uncorrelated errors is not possible and the amount of data is limited. First, a synthetic experiment has been carried out to assess the performance of CTC formulation. As an example of application, the error characterization of three collocated L-band brightness temperature (TB) measurements over land has been performed. Two of the datasets come from ESA (European Space Agency) SMOS (Soil Moisture and Ocean Salinity) mission: one is the reconstructed TB from the operational L1B v620 product, and the other is the reconstructed TB from the operational L1B v620 product resulting from application of an RFI (Radio Frequency Interference) mitigation technique, the nodal sampling (NS). The third is an independent dataset, the TB acquired by a NASA (National Aeronautics and Space Administration) SMAP (Soil Moisture Active Passive) radiometer. Our analysis shows that the application of NS leads to TB error reduction with respect to the current version of SMOS TB in 80% of the points in the global map, with an average reduction of approximately 1 K over RFI-free regions and approximately 1.45 K over strongly RFI-contaminated areas.


Introduction
Triple collocation (TC) analysis has been widely used for the quality assessment of many remotely sensed geophysical variables, such as ocean winds [1,2], soil moisture [3][4][5][6], and sea surface This paper is organized as follows. Section 2 details the proposed triple collocation formulation and the classical formulation used for comparison and describes the synthetic data used to assess the performance of the method and the remotely sensed L-band brightness temperature maps used in the example of application. The main results and findings are summarized in Section 3. Conclusions and perspectives from this work are given in Section 4. The complete development of the formulation used in this work can be found in Appendix A.

Settings and Notation
We have three series of spatially and temporally collocated measurements x i of the same variable θ: where a i is the so-called scaling calibration coefficients and b i is the measurement biases. As a matter of fact, if made precise by redefining the signal θ, we can just focus on the evaluation of the intercalibration factors α 12 and α 13 , where α ij = a i /a j . Each measurement has an unknown error, and we want to estimate the standard deviation of the errors δ i that we denote by σ δ i . As we have an explicit independent term b i to account for any nonzero mean, we consider the errors to be unbiased, so δ i = 0 and the errors are uncorrelated from the physical quantity, θδ i = 0. We assume that errors 1 and 2 are correlated, with a covariance φ 12 = δ 1 δ 2 , but the errors of these two variables are completely uncorrelated from the error of measurement 3 ( δ 1 δ 3 = 0 and δ 2 δ 3 = 0). The three datasets are assumed to resolve similar spatial scales, so the representativeness error is neglected. We also assume that all the variances and covariances are finite.
Calibration coefficients a i can be estimated by independent methods previously used for triple collocation, for example, by using an iterative linear least-squares approximation [12]. Some studies have been devoted to examining various methodologies for estimating the calibration coefficients and to assessing their performances in land applications [13,14].
In this work, measurements are assumed to have intercalibration factors equal to 1. Given that the intercalibration factor α 12 = a 1 a 2 can be computed directly from the covariances of x 1 and x 2 with x 3 (α 12 = s 13 s 23 ), we can use it to recalibrate the variable x 2 , and therefore, we can assume without loss of generality that a 1 = a 2 . We also assume that a 1 = a 3 , so α 13 = 1 (see Appendix A)).
We denote by s i (i = 1, 2, 3) the variances of the measurements: We also denote by s ij , i < j the covariances between the measurements i and j:

Least Squared Error Triple Collocation (LSETC)
The classical approach to estimate error variances is to use the Least Squared Error Triple Collocation (LSETC) [5], as discussed in the Appendix A. We can write a very simple linear system relating the 6 order-2 moments of the measurements with 5 unknowns (σ 2 θ , σ 2 δ 1 , σ 2 δ 2 , σ 2 δ 3 , and φ 12 ), in the following way: Notice that, while, by construction, the LSE solution grants that the average error is minimized, it does not ensure that the convergence with respect to the number of samples is the fastest one.
We have used the classical LSETC for comparison with the performance of the new triple collocation method presented in this work (Section 2.1.3).

Correlated Triple Collocation (CTC)
In this paper, we propose an alternative method for the case of two error-correlated datasets, the Correlated Triple Collocation (CTC). The full theoretical derivation of CTC can be found in the Appendix A. Let us explain here the central idea and the underlying assumptions and summarize the final formulas.
Statistical fluctuations can make the estimates of error variances in TC negative, something that may happen when the amplitude of the statistical fluctuation is larger than the variance of any or all of the errors. A negative estimate for a variance must be interpreted as the value of the variance cannot be estimated at that level of fluctuation. The formulas have been designed precisely to ensure that a maximum number of estimated error variances are positive by construction in order to converge faster to the real values of the error variances. This design makes our method less sensitive to statistical fluctuations and therefore best suited for the analysis of limited samples of data, since we can grant that at least two of the three measurements have positive estimates of their error variances (see Appendix A.4). This allows a reliable estimation of errors with scarce sampling sizes. Our method has a faster convergence with the sampling size to the real values of the error variances than the LSETC (see the synthetic results in Section 2.2).
We have looked for a linear transformation of x 1 and x 2 into two new variables x 1 and x 2 with errors δ 1 and δ 2 that are uncorrelated. We define the auxiliary scaling parameters u and v as follows: 12 ; v = s 1 − s 12 s 1 + s 2 − 2s 12 (6) and the order-2 moments of the uncorrelated-error variables x 1 and x 2 as: Then, we can compute the values of the error variances of the original variables by applying the rule of transformation for covariance matrices. The estimates for the error variances and the error covariance using CTC are given by the following expressions:

Generation of Synthetic Data
Synthetic data has been generated in order to assess the performance of the proposed CTC formulation as a function of the following: i The sampling size of the series of triplets, N. ii The correlation ρ 12 between the errors of the measurements x 1 and x 2 . iii The differences in the standard deviations of the three error measurements σ δ 1 , σ δ 2 , and σ δ 3 : Therefore, we have generated ensembles of simulated data according to these five parameters, N, ρ 12 , σ δ 1 , σ δ 2 , and σ δ 3 . We have defined three test cases or groups of datasets depending on the values of the error standard deviations: -Case 1 ("small uncorrelated"): The measurement with uncorrelated error has an error standard deviation significantly lower than the other two measurements (σ δ 1 = 0.5, σ δ 2 = 0.25, and σ δ 3 = 0.1). -Case 2 ("equal"): All errors are considered equal (σ δ 1 = σ δ 2 = σ δ 3 = 0.5). -Case 3 ("large uncorrelated"): The measurement with uncorrelated error has an error standard deviation significantly higher than the other two measurements (σ δ 1 = 0.1, σ δ 2 = 0.25, and σ δ 3 = 0.5).
The values of the signal θ have been randomly generated following a Gaussian distribution with a standard deviation of 1 and fixed mean (taken as 0, but it is irrelevant for the calculations). The measurement errors δ 1 , δ 2 , and δ 3 are also Gaussian with standard deviations according to each one of the test cases, and δ 1 and δ 2 are correlated with a prescribed covariance obtained from φ 12 . Notice that, in all cases, the standard deviations of the errors have been set to be smaller than those of the signal. Each synthetic dataset consists of 100,000 series of triple observations. This allows us to compute the average estimators of the error standard deviations δ i and the uncertainty of each one of those estimators (calculated as the standard deviation of the estimators over the ensemble of 100,000 realizations).

Analysis on the Intercalibration Factors
In this paper, we have assumed that α 12 = α 13 = 1. The first condition is not really restrictive for the type of measurements considered here (as explained in the Appendix A, we can estimate α 12 as α 12 = s 13 s 23 and use this estimate to renormalize the measurements x 2 . However, having an intercalibration factor α 13 = 1 would lead to poor estimates of all the errors.
The statistical fluctuations in the estimation of the intercalibration factors can lead to erroneous impression of having nontrivial intercalibrations. We have performed a synthetic experiment to show the impact of statistical fluctuations on the estimations of the intercalibration factors. For that, we have generated triples of measurements with no biases and trivial calibration factors (a i = 1, b i = 0 ∀i). The errors have been generated as completely uncorrelated Gaussians and with the same error standard deviations defined in Section 2.2 (case 1: σ δ 1 = 0.5, σ δ 2 = 0.25, and σ δ 3 = 0.1; case 2: σ δ 1 = 0.5, σ δ 2 = 0.5, and σ δ 3 = 0.5; and case 3: σ δ 1 = 0.1, σ δ 2 = 0.25, and σ δ 3 = 0.5), and in all cases, we consider the signal as a Gaussian with standard deviation σ θ = 1). For each sampling size from N = 50 to N = 1000, we have generated 100,000 series of triplets, and for each series, we estimate the intercalibration factors α 12 and α 13 , as explained in the Appendix A.1:

L-Band Brightness Temperatures over Land
As an example of application, we have applied CTC to maps of L-band brightness temperatures acquired by two satellites.

Nodal Sampling: Reduction of RFI Contamination in SMOS Images
Contamination by RFIs is still an important source of error in SMOS TB [16,17]. The application of NS technique to TB images aims to mitigate the degradation produced by RFIs [10,11]. NS has been thoroughly tested and validated over oceans, both at brightness temperature and Sea Surface Salinity (SSS) levels [11,18]. In terms of TBs, NS performance has been assessed by comparing measured ocean TBs with the modeled TBs derived by evaluating the geophysical model function presented in [19] with some geophysical priors (SSS provided by World Ocean Atlas 2009 climatology and sea surface temperature and wind speed provided by ECMWF (European Centre for Medium-Range Weather Forecasts)) [18].
At the SSS level, global comparisons (60 • S-60 • N) with in situ data provided by Argo floats show a significant improvement in NS SSS wit a Root Mean Square (RMS) equal to 0.81 with respect to the SSS from the current operational TBs (RMS = 1.09). Global comparisons to data from ship tracks show also a slight improvement in terms of RMS when applying NS (RMS = 1.82 and a correlation coefficient of 0.61) with respect to the SSS from the current operational TBs (RMS = 1.87 and correlation coefficient R = 0.55).It must be pointed out that the spatial coverage of SSS retrievals from the corrected NS TBs increases by 20% on average with respect to the SSS from the current operational TBs, particularly over semi-enclosed seas and strongly RFI-contaminated regions [20].
The application of NS could be of special relevance over land, since most of the RFI sources are located there and their impact is more noticeable in soil moisture retrievals. The comparison of measured TB to a forward model over land is not as straightforward as in the case of the ocean since footprints over land are much more heterogeneous and the model depends on many parameters [21,22]. Therefore, we propose to apply TC analysis to evaluate the impact of applying NS to reduce TB errors over land. In this work, we use three collocated TB datasets: SMOS nominal TB (reconstructed from the operational L1B product v620), SMOS NS TBs (reconstructed also from L1B product v620 and applying NS), and TBs measured by a SMAP radiometer. CTC formulation is required, since we are analyzing two datasets from SMOS acquisitions, and therefore, their errors are correlated.

SMOS Brightness Temperatures
The operational SMOS L1B product (v620) has been used in this analysis [23,24]. A similar procedure as the one used in the operational SMOS Level 1C processing has been applied to obtain the brightness temperatures in the antenna reference frame from the L1B product [25]. The difference is that the operational processor uses an hexagonal grid of 128 × 128 points instead of the grid of 64 × 64 points that is used in this processing to reduce the computational cost without loss of information. In fact, the lowest resolution to transform visibilities into brightness temperatures is set by the MIRAS instrument characteristics and given by (NT × NT), with NT = 3 · NEL + 1, where NEL = 21 is the number of elements per arms. Using an hexagonal grid of 64 × 64 points ensures the maximum independence among measurements of the same snapshot [26].
In the case of applying NS, the algorithm starts also from the L1B v620 product and the output is directly the reconstructed TBs in the antenna reference frame [10]. RFI flags provided in the L1B products [24] have not been applied, since we are interested in assessing the impact of NS, particularly over RFI-contaminated regions.
TBs have been transformed from antenna frame to surface reference frame, which includes a geometrical rotation due to rotation of the polarization basis and a Faraday rotation due to the presence of the ionosphere [27]. Attenuation of the atmosphere has been corrected to obtain brightness temperatures at the Bottom of the Atmosphere (BOA), that is, at land surface. A Lambert Azimuthal equal area projection [28] is used for generation of the georeferenced TBs.
Notice that SMOS observations are multi-angular (0-65 • ), while SMAP measures at a constant incidence angle (40 • ). Therefore, the error characterization by TC can be done only for a reduced range of incidence angles taken from SMOS to fit that of SMAP. Therefore, TB measurements acquired under the incidence angle range between 37.5 • and 42.5 • and, in the alias-free region of SMOS field-of-view, have been used to generate SMOS TB maps. Since we are not exploiting the multi-angular capability of SMOS, we expect higher errors in SMOS TB than SMAP TB (due to the worse radiometric sensitivity of SMOS at that incidence angle range).
Several studies have been carried out to compare SMOS and SMAP TB measurements in RFI-free regions, showing comparable performances [29,30]. However, it is expected that SMOS TB errors increase in RFI-contaminated areas since SMOS does not have RFI mitigation on board. SMAP is much less affected by this contamination thanks to its on-board RFI mitigation back end [31,32].

SMAP Brightness Temperatures
The L1B Radiometer Half-Orbit Time-Ordered Brightness Temperatures, Version 3 product (Composite Release Identifier R13080) [33] has been used for generation of the SMAP TB maps. The horizontally and vertically polarized brightness temperatures at the surface of the Earth after on-board RFI filtering are used in this analysis, averaging the fore-and aft-look samples [34]. No additional filter has been applied to SMAP TBs, that is, all measurements have been considered (also if the on-board software was unable to correct the brightness temperature when an RFI is detected).

Spatiotemporally Collocated TB Maps
Global 3-day maps have been generated every day by interpolating the SMOS and SMAP TB data to the nearest neighbour in a 0.25 • regular longitude-latitude grid. The temporal resolution of the maps (3 days) has been selected as a trade-off between achieving global spatial coverage and having enough long time-series for triple collocation analysis. SMOS and SMAP data from ascending and descending overpasses have been used in generation of the maps. The time period considered ranges from January 2016 to September 2017, so 628 time-overlapping (i.e., 210 independent) 3-day maps are available for triple collocation analysis.
It is assumed that the spatiotemporal averages of the SMOS and SMAP data as well as the averages in a range of incidence angles in SMOS that fit the SMAP incidence angle are representative of the same geophysical quantity (the "typical" or average brightness temperature in the given 0.25 • × 0.25 • pixel and for the given 3-day period). There could be a significant difference between the resulting 3-day SMOS and SMAP TB maps because of the specific sampling time of each instrument and the geophysical variability itself. Such differences do not affect our analysis, as they are simply accounted for as part of the associated error measurement.

Effective Spatial Resolutions of SMOS and SMAP TB Maps
One of the required conditions to apply CTC is that the three measurements being compared must be defined over similar spatial scales. We have assessed the scales of the three TB datasets (SMAP, SMOS nominal, and SMOS NS) by analyzing their Power Density Spectra (PDS) [35]. We have computed their PDS for two large regions, differently orientated (zonal or meridional). The first kind of PDSs (labelled "Central Asia") is defined by the zonal transects extending from 45 • to 110 • W; the PDS of those transects are averaged over a range of latitudes going from 40 • to 60 • N and over the full time extension of the maps. The second kind of PDS (labelled "Central Africa") is defined by the meridional transects extending from 30 • S to 30 • N; the PDS of those transects are averaged over a range of longitudes going from 15 • to 30 • W and over the full time extension of the maps.
The obtained PDSs, S(k), are expected to behave as a power law of the wavenumber, k. Therefore, we expect S(k) ∼ k −β , where the value of the scaling exponent β is expected to be between 1 and 2. This means that, when representing the PDS vs. k in a log-log plot, we should observed a straight line with a slope exactly equal to −β.
The behaviour of the PDS is expected to deviate from the pure power law in the presence of resolution-limiting effects. It may happen that the measurements have a certain level of white noise, so the actual behavior of the PDS would then be S(k) ∼ k −β + η, where η is the amplitude level of the white noise, which is a constant, independent of the wavenumber k. Therefore, for wavenumbers large enough, the power-law term becomes negligible and we would then have S(k) ∼ η. We can determine a threshold wavenumber k 0 such that, for k < k 0 , we have with a good approximation S(k) ∼ k −β (therefore, the physical scales are well resolved) and, for k > k 0 , we have with a good approximation S(k) ∼ η (therefore, the measurement is dominated by noise and no physical scales are utterly resolved). This threshold k 0 is the resolution wavenumber, which defines the resolution wavelength λ 0 = k −1 0 . Since at a wavelength we can resolve two points (one assigned to the positive part of the sinus, and the other assigned to the negative part), we can define the resolution scale as Noise is the most common effect limiting the effective resolution of data, but there are others. Low-pass filtering effectively suppresses noise but lower scales as well, and therefore, it also limits the effective resolution of the data. In terms of PDS, the impact of applying a low-pass filter is evidenced by a sudden drop of value for wavelengths k beyond a given threshold value k 0 . Exactly as in the previous case, this threshold wavenumber k 0 serves to determine the effective resolution of the data.

Synthetic Experiments on Error-Correlated Triplets
We have defined three metrics to assess the performance of the CTC and LSETC methods:

•
Fraction of valid retrievals is the ratio of the total valid retrievals (that is, nonnegative estimates of the error variances σ 2 δ i ) to the total number of realizations. The closer to 1, the better.

•
Bias is the difference between the average of all valid estimates of the error standard deviations σ δ i and the value used for the generation of the dataset. The closer to 0, the better. It provides the bias in our estimates of σ δ i . Positive bias indicates that the error is overestimated, and negative bias indicates that it is sub-estimated.

•
Uncertainty is the standard deviation of the valid estimates of error standard deviations. The closer to 0, the better. It provides the accuracy in our estimates of σ δ i .
The results of the synthetic experiments are shown in Figures 1-3. Each panel corresponds to a given value of the sampling size N (per rows). We show each one of the three metrics per columns. In all panels, the x-axis corresponds to the error correlation ρ 12 . Solid lines are for CTC, while dashed lines correspond to LSETC.  Quality assessment for test case 2 ("equal") (σ δ 1 = 0.5, σ δ 2 = 0.5, and σ δ 3 = 0.5): the results are shown as a function of the correlation parameter ρ 12 . From left to right: Fraction of valid retrievals, normalized mean, and normalized uncertainty. From top to bottom: Results for N = 50, N = 100, N = 500, and N = 1000. Grey color corresponds to the measurement x 1 , purple is for measurement x 2 , and red is for measurement x 3 . Solid lines are the results for correlated triple collocation, while dashed lines are for least squared error triple collocation. The main findings and conclusions from the analysis with synthetic data (Figures 1-3) can be summarized as follows: • The fraction of valid points is very large even for scarce samplings (N = 50) for the measurements with the largest error standard deviations. The number of valid retrievals for the "small uncorrelated case" and "large uncorrelated case" is lower for the measurement with the lowest error standard deviation: in the range of 60% for scarce sampling and increasing slowly for larger sampling sizes. CTC has in general a larger number of valid retrievals than LSETC, especially in the "small uncorrelated case" and less in the "large uncorrelated case". The fractions of valid points for LSETC and CTC are very similar in the "equal case".

•
Biases are not very large in the "small uncorrelated case" and "equal case". Even for scarce samplings (N = 50), they are at most about 10% of the largest error standard deviation for the CTC and about 20% for the LSETC. The situation is worse in the "large uncorrelated case", where it has 30% of the largest error standard deviation (both for CTC and LSETC) for scarce sampling and only attains 10% for good sampling (N = 500) or better. In most cases, the performance in terms of biases of CTC is better than that of LSETC.

•
The measurement with the smallest error standard deviation has always a positive bias in CTC, indicating that its error standard deviation is always overestimated. This bias is reduced rapidly as sampling size N increases. In the "equal case", biases are negligible for the three measurements even for scarce sampling. • Uncertainties are small in the "small uncorrelated case" and moderate in the other two cases.
In the two latest cases, we expect uncertainties to be around 10% of the largest error standard deviation even with excellent samplings (N = 1000). CTC outperforms LSETC, especially in the "small uncorrelated case". • From the experiments, we see that the dependence of all metrics on the value of the error correlation ρ 12 is weak in most cases. For the "small uncorrelated case", the bias and uncertainty decrease at high correlation values for CTC, since the two measurements with larger errors become essentially the same, but in all cases, CTC outperforms LSETC. Hence, CTC is very robust independently of the degree of correlation between those errors.
Based on all previous points, we conclude that CTC is an optimal method for assessing TC errors in the triplets formed by SMOS nominal, SMOS NS, and SMAP TBs. This corresponds to the "small uncorrelated case (since the measurement with uncorrelated error, SMAP, has an error standard deviation significantly lower than the SMOS measurements). In this case, CTC outperforms LSETC in the convergence speed and in terms of bias and standard deviation. Therefore, we will take CTC as our reference triple collocation method for the characterization of radiometric errors in satellite L-band brightness temperatures over land.

Impact of Statistical Fluctuations on the Estimation of Intercalibration Factors
For each sampling size N, we have computed the mean and standard deviation of the intercalibration factors over the ensemble of 100,000 realizations. While the means are very close to 1 (not shown), the standard deviations are not negligible (see the log-log plot in Figure 4) and decay slowly with sampling size 1/ √ N. This is expected for fluctuations according to the Central Limit Theorem. The (extrapolated) value for N = 1 determines the factor multiplying 1/ √ N, and this factor determines the vertical positioning of the curves. As shown, those factors depend on the values of the error standard deviations and significantly differ for the different cases. The standard deviation of the intercalibration factors is typically around 0.1 for N = 50 (that is, an error of ≈10% around the right value) and even around 0.02 or 2% for N = 1000. Notice that those curves have been obtained from simulated data assuming Gaussian distributions for both the errors and the signal. Dispersion in the estimation of the intercalibration factors could be even larger if the distribution of the signal significantly deviates from a Gaussian distribution and if the measurement errors are correlated.  We conclude that statistical fluctuations, even with moderately large sampling size, are important enough to give the false impression that the measurements are poorly intercalibrated.

Sensitivity Analysis of the Estimated Error Variances to Changes in the Intercalibration Factors
We have carried out an experiment to analyze the impact of errors in the intercalibration factors on the retrievals obtained by the application of CTC. As discussed earlier, the factors a 1 and a 2 can be adjusted to be coincident, so the only factor missing is a 3 (or the intercalibration factor α 13 ).
We have generated a series of synthetic data to assess the impact of having an intercalibration factor α 13 significantly differing from one. For that, we have computed a dataset of series for each value of α 13 between 0.8 (20% lower) and 1.2 (20% larger) in intervals of 0.01. We have generated 100,000 series of synthetic signals, each one with a sampling size of N = 500, and the correlation coefficient between measurements 1 and 2 (ρ 12 ) has been kept fixed to 0.5.
To assess the results, we have taken as metrics the fraction of valid retrievals and the biases of the estimates.
As shown in Figure 5, the bias is a decreasing function of α 13 for the error-correlated measurements x 1 and x 2 , while it increases for the uncorrelated measurement x 3 . For x 3 , the relation is almost perfectly linear when σ δ 3 = 0.5 (that is, the "equal" and "large uncorrelated" cases in Section 2.2). Probably, the excess or lack of signal associated to the change in the intercalibration factor is attributed to the error of x 3 measurement, with a similar change in magnitude but of opposite sign for measurements x 1 and x 2 . Therefore, the impact on the estimates of the error standard deviations will be larger when the variance of the signal (here, conventionally fixed at 1) is larger. Bottom: Biases on the retrieved error standard deviations as a function of α 13 . The columns from left to right are for the "small uncorrelated", "equal", and "large uncorrelated" cases defined in Section 2.2.
As evidenced by the curves of valid retrievals (bottom plots in Figure 5), when the bias induced by the change in α 13 is close to a critical value when σ δ i would go to zero, there is a sudden drop in the fraction of valid retrievals. Therefore, the number of valid retrievals is a good indicator of having problems with the intercalibration factors. The transition points depend on the ratio between the error standard deviation and one of the signal. In our figures, the drops happen at percentages of variation of the intercalibration factor that are equal to the standard deviation of the error of the affected measurement; this is not a coincidence but a consequence of the variance 1 for the signal. A signal with larger variance will have transition points at lower values of the percentage of miscalibration. When working with real signals, as the variance of the errors is typically much smaller than that of the signal, the range of admissible errors in the intercalibration factor is quite narrow in practice.

Error Characterization of Satellite L-Band Brightness Temperatures over Land
The triple collocation analysis has been performed for TB maps at H-polarization, V-polarization, and for the half of the first Stokes parameter. Only results for the latter are shown, but similar results have been obtained for the H and V polarizations.
In order to assess if the datasets resolve similar spatial scales, we have computed the PDS for the Central Asia and Central Africa regions for the three datasets ( Figure 6). The three curves are very close to each other for large scales and are in relatively good agreement with the expected power law. The behaviour of the zonal transect (Central Asia) and of the meridional transect (Central Africa) are slightly different probably due to the RFI contamination over Asia, which is filtered on-board in the case of SMAP and not in SMOS. For Central Asia, the PDSs start to slightly diverge at about k = 0.3 (degrees) −1 . Then, close to 1 (degrees) −1 , the three curves reach their limit: both SMOS PDSs become essentially horizontal while SMAP experiences a serious drop of value, indicating that it is a smoother product. In the case of Central Africa, the three curves are packed closer together, but again, at about 1 (degrees) −1 , both SMOS curves become rather horizontal while SMAP undergoes a drop, although less sharp than in the other case. We conclude that, for the three datasets, the threshold wavenumber is about 1 (degrees) −1 , which means that their effective resolution is approximately the same, corresponding to 0.5 • . This is consistent with the theoretical spatial resolutions for SMAP (37 km × 49 km) and SMOS (35 to 50 km). The intercalibration factors have been estimated by the ratio of the average TBs for the time-series between the two corresponding datasets. The intercalibration factor between SMOS nominal and SMOS NS (α 12 ) is very close to the one for most of the points in the map (see the left column in Figure 7, mean value: 0.997). Differences are mainly concentrated around coastlines and ice edges, where the NS TBs are systematically higher than the nominal ones (due to known residual systematic biases present in NSv2 [36]) and in the location of RFI sources, as reported in [11]. Analyzing the ratio between SMOS NS and SMAP (α 13 ), differences are more noticiable over Europe and China due to the higher RFI contamination in SMOS TB (see the right column in Figure 7). In any case, the global mean of α 13 is 1.014, indicating a good overall agreement between SMOS NS and SMAP measurements. Notice that the intercalibration factors are close to 1, with deviations in the expected range due to statistical fluctuations (see Section 3.1.1). The error correlation parameter ρ 12 between the two SMOS datasets (see Figure 8) has been estimated following the equations summarized in Section 2.1.3. It can be appreciated that errors are highly correlated for most of the points: the mean correlation is 0.71, and over wide areas of the globe, ρ 12 is very close to 1 and, in all cases, positive. Those places where the correlation significantly drops seem to be related to the presence of intense RFIs, which can be corroborated looking at the map of SMOS RFI probability in Figure 9. This RFI probability has been computed following the procedure detailed in [37] for the period January 2016-September 2017.
Global maps of the error standard deviation for each one of the three TB datasets are shown in Figure 10. Notice the low estimated errors over Antarctica, which can be explained by the very low temporal and spatial TB variability over that region. The application of NS has led to an overall significant reduction of TB errors with respect to SMOS nominal TB, particularly over strongly RFI-contaminated regions. NS errors are more concentrated around the locations of the RFI sources (see Figure 9). This reduction over Europe, Arabian, India, and China can be more clearly appreciated in the zoom of Figure 11. It is important to point out that the error estimation of SMAP was not valid (negative values; see Appendix A) over 15% of the gridpoints (Figure 10c). In the case of the SMOS nominal and SMOS NS, the number of data gaps (nonvalid estimates) is very low (1.02% and 1.34%, respectively). This is consistent with the findings in our experiment with synthetic data case 1 in Section 2.2. As we show below, SMAP error standard deviation is smaller than that of the other two and the number of valid retrievals is approximately 85% for 210 samples, which is in between 60% of retrievals for N = 100 samples and approximately 80-90% for N = 500 obtained in the synthetic experiments.  The map of the differences between the error standard deviation of SMOS nominal and SMOS NS (Figure 12a) confirms that the application of NS effectively reduces the RFI contamination over land. The exception is for gridpoints close to the coastlines and those in the locations of RFI sources (see the blue spots and their correspondence to the RFI probability map in Figure 9), which present larger errors in NS TBs. These are known limitations of NSv2, as published in [11,36]. The residual contamination in coastal (and in ice edges) gridpoints is expected to be largely reduced by applying a refined version of the NS algorithm, very recently developed to improve the performances in land/ocean/ice transitions. Looking at the difference between the error standard deviation in SMOS nominal and SMAP (Figure 12b), the highest differences are found over RFI-contaminated regions but are also noticiable over RFI-free areas. These differences are significantly reduced when analyzing the difference between the error standard deviation in SMOS NS and SMAP (Figure 12c). Errors in NS TBs are larger than those of SMAP in RFI-contaminated regions, which is expected since RFI contamination in SMAP is much more moderate thanks to its on-board digital detector back end [31,32]. For some RFI-free regions, NS leads to lower errors than SMAP (see bluish regions in Figure 12c). The smaller error standard deviation of SMAP with respect to SMOS was expected because the radiometric accuracy in the averaged SMOS incidence angle range (37.5-42.5 • ) is worse than that of SMAP. However, it is important to notice that SMOS final retrieved error is comparable to SMAP when considering the whole incidence angle range. What is actually significant is that, with the application of NS (which is known to lead to an improvement of the effective radiometric accuracy [11]) over certain areas, the SMOS data acquired close to SMAP incidence angle can overcome SMAP accuracy.  Figure 10.
Histograms of the differences between error standard deviations of SMOS nominal and SMOS NS are shown in Figure 13a for the global map and in Figure 13b for the region that is more affected by RFI contamination, covering Europe, Asia, and north Africa (longitude range (20 W, 180 E), latitude range (10 N, 80 N)) (see Figure 9). Statistics show that NS reduces SMOS TB errors over strongly RFI-contaminated regions (mean error reduction approximately 1.5K) and generally improves TB quality (approximately 75% of the gridpoints over land present lower errors for NS TBs than for the nominal ones). Similar statistics are obtained when considering the global map, although the mean error reduction is decreased. The impact of NS is lower over RFI-free regions, but it is still beneficial since it reduces overall ripples in TB images.
In the case of SMAP, error estimates for 15% of the gridpoints have not converged to a meaningful value (white shading in Figure 10c). A refinement of the TC analysis is proposed in Section 3.2.1 to obtain error estimates for the global map.

Inferring SMAP Errors Overestimate Gaps
Looking at the maps of the standard deviation of the SMAP TBs ( Figure 14a) and their error standard deviation (Figure 10c), it is evident that both variables are related by an increasing function. Over 85% of points have a valid estimate of σ δ 3 ; we computed the histogram of σ δ 3  As a strategy to improve the estimation of σ δ 3 over the gaps, we computed the parameters of a linear regression between σ δ 3 and σ 3 . We performed a weighted least square fit, where each value of σ 3 was weighted by the inverse of the conditioned variance. The resulting linear relation is σ δ 3 = aσ 3 , where the empirical value of a is a = 0.42.
Assuming the linear relation between the error variance and the measurement variance, the variance of the geophysical signal can be estimated as follows: Then, the covariance of the combined SMOS and SMAP variables, s 23 , which is in fact an estimate of σ 2 θ , can be substituted in Equation (8) by the value yielded by Equation (10). Error standard deviations for the three measurements have been computed (hereafter, adjusted error standard deviations) with the proposed refinement ( Figure 15). Adjusted error estimates are very similar to the error standard deviations estimated in the initial analysis (Section 3.2) but with almost complete spatial coverage. Note that mean adjusted errors have slightly increased with respect to the previous mean error estimates (see the mean values in the captions of Figures 10 and 15, respectively) and with an improved coverage (less gaps). NS TBs show an overall decrease of errors both in RFI-contaminated regions and over RFI-free regions (see Figure 16). Statistics in the histograms of Figure 17 confirm that the mean error reduction of SMOS NS with respect to SMOS nominal TBs is very similar to the one obtained in the first iteration of CTC (Section 3.2, Figure 13). NS leads to an improvement of the TB quality in 79% of the gridpoints over land.

Conclusions
A new formulation of the triple collocation analysis has been developed for the specific case of three datasets resolving similar spatial scales and presenting correlated errors for two of them. The CTC method has been designed to be as less data demanding as possible, and it has been first assessed using synthetic data. The performance of this method has been compared with the LSETC: CTC attains a higher fraction of valid retrievals and outperforms the classical LSETC in terms of biases and uncertainties of the error estimates, specially for scarce to moderate sampling sizes. In the "equal case" (similar errors in the three measurements) and in the "large uncorrelated case" (the measurement with uncorrelated error has an error standard deviation significantly higher than the other two), both methods perform quite similar, although CTC attains a higher fraction of valid retrievals, specially in scarce samplings.
This methodology can be particularly beneficial for the error characterization of variables for which getting measurement systems with uncorrelated errors is challenging or not feasible. The method is thought to be especially suited for its application to triplets of remote sensing maps that resolve similar spatial scales, where having three collocated, completely independent remote sensing maps is unlikely.
As an example of the application of this method, we have used it for the characterization of radiometric errors in L-band brightness temperatures over land. The triplets formed by SMAP, SMOS nominal, and SMOS NS TBs have been specifically analyzed for evaluating NS performances over land as an RFI mitigation technique. CTC method has been shown to provide spatial maps of errors within a reasonable margin of uncertainty even with a moderate sampling size.
The triple collocation analysis has revealed that TB errors are significantly reduced when NS is applied with respect to SMOS nominal data, attaining a mean error reduction of approximately 1.5 K in strongly RFI-contaminated regions and approximately 1 K globally. The application of NS has been shown to improve the quality of TBs in 79% of the gridpoints of the map. Most of the gridpoints where nominal TBs present lower errors than NS are in the precise locations of the RFI sources (measurements that in any case will be discarded before geophysical retrieval) and along the coastline. The increase of the systematic biases very close to the coasts and ice edges is a known performance issue with NSv2 [36]. Precisely, a refined version of NS focused on reducing these systematic biases in the coasts has been recently developed, and it is currently under extensive validation.
The extension of the present analysis to global TB maps (including ocean and ice) is ongoing, with the aim to assess the NS impact on TB independently on the target. The same CTC methodology can be applied to the retrieved geophysical parameters (both soil moisture and ocean salinity). We are currently working in the error characterization of different L-band sea surface salinity products (including SMAP, SMOS, and reanalysis data) with promising results [38].
The CTC technique is currently being used by the SMOS Level 1 Expert Support Laboratories as a metric for assessing the performance at TB level of the next versions of the SMOS L1 processor with respect to the current one. The TC presented in this work is a very valuable tool since it can provide insights about the improvements/limitations of the changes proposed at L1 without the need for a forward model or for further retrieval of the geophysical parameters and validation of them. In this context, the performances of a new extrapolating technique of the TB frequencies outside the star coverage of SMOS to mitigate RFI contamination effects are also being evaluated [39].  A very coarse estimate of the errors could be obtained if we consider that all calibration factors a i are equal to 1 (after a proper normalization) and just estimating the signal variance σ 2 θ from any measurement covariance, for instance, s 12 ; we would therefore have σ 2 δ i = s i − s 12 . However, statistical fluctuations due to finite sampling would make this kind of estimator noisy and less reliable. A better estimation of the signal variance could be obtained by averaging all the covariances,σ θ 2 = 1 3 (s 12 + s 13 + s 23 ) and then by applying σ 2 δ i = s i −σ θ 2 . Although this is a better estimate, it is also very affected by the fluctuations associated to finite sampling and particularly by errors in the calibration factors a i . If we assume that the calibration factors are not known, it is still possible to estimate the error variances, and thus, any possible error in the calibration factors can be neglected. In such a case, we have more unknowns than equations. It is easy to verify that the final expression for the most general case is as follows: An additional advantage of this expression is that it makes optimal use of the statistics, so the fluctuations due to finite sampling are smaller and the estimates are more accurate.
We can estimate the intercalibration factors α ij = a i /a j by dividing the covariances of both variables with the third one:

. Two Measurements with Correlated Errors but Uncorrelated from a Known Third Measurement: Least Squared Error Triple Collocation
One of the approaches that allows us to apply the TC in the case of two measurements with correlated errors but uncorrelated from the error of the third one is the least squares solution, which minimizes the squared error of the estimates [9]. Assuming that all calibration factors a i are equal to 1, we can write a very simple linear system relating the 6 order-2 moments of the measurements with 5 unknowns (σ 2 θ , σ 2 δ 1 , σ 2 δ 2 ,σ 2 δ 3 , and φ 12 ) in the following way: Defining the vectors s t = (s 1 s 2 s 3 s 12 s 13 s 23 ) and Ω t = (σ 2 θ σ 2 δ 1 σ 2 δ 2 σ 2 δ 3 φ 12 ), the equation above can be written as follows: where A is a 6 × 5 matrix given by the following: Notice that this is an overdetermined system, so only if the order-2 moments of the measurements are perfectly well determined will there be a solution. In general, fluctuations associated with having finite samples of data will lead to an incompatible linear system with no algebraic solution. The least squared error estimate is based on minimizing the error of the estimates [5], and it equates to applying a simple pseudo-inverse to Equation (A6): Applying the least squared error estimate, Equation (A8), to our matrix A, Equation (A7), we have and from where we obtain the least squared error solution, given by As a matter of fact, the least squared error solution, Equation (A11), is just one the possible solutions to the system of equations posed in Equation (A5) because it is an overdetermined system in which one of the equations must be eliminated (directly or by a proper combination of all equations); in fact, the LSE equation is just the solution of the determined system obtained by averaging the last two equations of Equation (A5). While, by construction, the LSE solution grants that the average error is minimized, it does not ensure that the convergence with respect to the number of samples is the fastest one. This is because some of the estimated error variances can be negative and must then be discarded, slowing the convergence. A possible way to circumvent this problem is to look for a different solution in which we aim to have as many positive estimates of the error variances as possible. This is precisely the basis of our CTC method.

Appendix A.3. Two Measurements with Correlated Errors but Uncorrelated from a Known Third Measurement: Correlated Triple Collocation
Let us now suppose that φ 12 = 0. We will apply the expression for estimating the error variances when they are completely uncorrelated, Equation (A3), to obtain a first estimate of the error variances. We will denote those estimates by ω i . Applying Equation (A3), we have the following: Given that the intercalibration factor α 12 = a 1 a 2 can be computed directly from the covariances of x 1 and x 2 with x 3 (α 12 = s 13 /s 23 ), we can use it to recalibrate the variable x 2 , so that we can assume without loss of generality that a 1 = a 2 and the error estimates are then simpler: We look for a linear transformation of x 1 and x 2 into two new variables x 1 and x 2 with errors δ 1 and δ 2 that are uncorrelated. The linear transformation will be defined by means of a 2 × 2 matrix U, The variables x 1 and x 2 are defined as follows: Substituting the expressions of the variables x 1 in functions of θ and δ i , we obtain where the square brackets in the center equations stand for the corresponding variables in the rightmost equations. The variables x 1 and x 2 will be error-decorrelated if δ 1 and δ 2 are uncorrelated. Therefore, let us analyze the covariance of δ 1 and δ 2 ; it can be expressed in terms of σ δ 1 , σ δ 2 , and φ 12 as follows: Using the expression for the error variance estimates in Equation (A13), we obtain the following: δ 1 δ 2 = u 11 u 21 ω 1 + u 12 u 22 ω 2 + (u 11 u 21 + u 12 u 22 + u 11 u 22 + u 12 u 21 )φ 12 (A18) If we impose that then the covariance of δ 1 and δ 2 depends only on known quantities (ω 1 and ω 2 ). Therefore, the elements of the matrix U can be tuned to ensure that this error covariance is cancelled. One possible choice of U is the following one: which is quite convenient as its determinant is exactly equal to 1, and then, the resulting variables after applying the linear transformation have the same order of magnitude. Besides, with this selection of variables, x1 = x1 − x2 (see Equation (A15)) only contains the error terms.
Therefore, the inverse of the U matrix yields the following: If we now define the variables x 1 , x 2 , and x 3 such that the first two are given by Equation (A16) with the particular choice of U given by Equation (A20) and x 3 = x 3 , we will be able to estimate the error variances of those three variables, obtaining the error variances of δ 1 , δ 2 , and δ 3 = δ 3 . If we denote these three variances by σ 2 δ 1 , σ 2 δ 2 , and σ 2 δ 3 , we can compute the values of the error variances of the original variables by applying the rule of transformation for covariance matrices: and Notice that, as far as σ 2 δ 1 and σ 2 δ 2 are positive, the shape of Equation (A22) ensures that σ 2 δ 1 and σ 2 δ 2 are also positive. As a bonus, we obtain an estimate of φ 12 , the error covariance of variables x 1 and x 2 .
In order to compute the error variances of x i , we cannot directly apply Equation (A3) because our choice of matrix U verified that u 11 + u 12 = 0 and, thus, according to Equation (A16), we have that a 1 = 0, that is, x 1 does not depend on θ. This makes s 13 = s 23 = 0, which means some fractions in Equation (A3) are ill-defined. Looking at the second-order moments of the variables x i , we have that they can be related to the prime error variances: Notice that, for our choice of U, we have that u 21 + u 22 = 1 and this implies, reading Equation (A16), that a 2 = a 2 = a 1 (the last equality assumed at the beginning of this section, although as we explained, a 2 can be adjusted to verify this condition). Therefore, we can estimate the prime error variances from the prime moments: where we have used a 3 = a 3 and α 23 = a 2 /a 3 . We can only solve for σ δ 2 and σ δ 3 if we assume that a 1 = a 2 = a 2 = a 3 = a 3 . Therefore, the problem of two correlated error variances can be solved if we assume that the three variables have the same scale calibration. Notice also that, by construction, the three estimates of error variances given in Equation (A24) must be positive.
Hence, once the prime moments are known, we can compute the prime error variances using Equation (A24) and we can finally compute the original error variances and the covariance between errors 1 and 2 using Equations (A22) and (A23). Therefore, to proceed, we need to express the moments of the prime variables in terms of those of the original variables by looking at the definition of Equation (A16): 11 s 1 + u 2 12 s 2 + 2u 11 u 12 s 12 s 2 = u 2 21 s 1 + u 2 22 s 2 + 2u 21 u 22 s 12 s 3 = s 3 s 12 = u 11 u 21 s 1 + u 12 u 22 s 2 + (u 11 u 22 + u 12 u 21 )s 12 s 13 = u 11 s 13 + u 12 s 23 s 23 = u 21 s 13 + u 22 s 23 (A24) For our particular choice of the matrix U, Equation (A20), we have (ω 1 +ω 2 ) 2 s 2 + 2 ω 1 ω 2 (ω 1 +ω 2 ) 2 s 12 s 3 = s 3 s 12 = ω 2 ω 1 +ω 2 s 1 − ω 1 ω 1 +ω 2 s 2 + ω 1 −ω 2 ω 1 +ω 2 s 12 s 13 = s 13 − s 23 s 23 = ω 2 ω 1 +ω 2 s 13 + ω 1 ω 1 +ω 2 s 23 (A25) The notation can be simplified to ease the use of collocated triple collocation; we can work out Equations (A22)-(A25) to a more usable form, expressing them in terms of the variances of the measurements s i and their covariances s ij . Let us define the variables u and v (that correspond to the previous u 21 and u 22 ): u = s 2 − s 12 s 1 + s 2 − 2s 12 ; v = s 1 − s 12 s 1 + s 2 − 2s 12 (A26) which is the result of substituting the actual values of ω 1 and ω 2 . Then, we have When data is scarce, we can have significant statistical fluctuations depending on the distribution of the error and of the signal itself. A statistical fluctuation is a large deviation from the mean value of a variable. From a statistical point of view, a statistical fluctuation can be viewed as a random variable caused by the particular sample of the data we are taking. Fluctuations tend to be smaller as we average larger and larger amounts of data, with the typical amplitude of a fluctuation as the inverse of the square root of the amount of averaged data.
Statistical fluctuations can make the estimates of error variances in TC negative, something that may happen when the amplitude of the statistical fluctuation is larger than the variance of any or all of the errors. A negative estimate for a variance must be interpreted as that value of the variance not being estimated at that level of fluctuation. This makes the estimation of the error variance harder when this error variance is much smaller than the variance of the signal itself or the variances of the other errors. In extreme cases, due to the coupling among all moments appearing in TC equations, the estimate of all error variances can be impossible.
The equations introduced in Appendix A.3 have been designed precisely to ensure that a maximum number of estimated error variances are positive by construction. This design makes our method less sensitive to statistical fluctuations and therefore best suited for the analysis of limited samples of data. Otherwise said, our method has fast convergence with the sampling size to the real values of the error variances (see the synthetic results in Section 2.2). As we will show in the following, we can grant that at least two of the three prime variables have positive estimates of their error variances.
First of all, it is worth noting that the estimators of the original error variances will always be positive as far as the estimated prime error variances are positive, disregarding the actual values of ω 1 and ω 2 . This is a consequence of Equations (A22) and (A23).
By construction, the estimate of σ 2 . This may happen because the determination of the second-order moments of the variables are computed with finite samples of data, and as the variance of the error is typically much smaller than that of the signal, a small finite-size fluctuation can lead to a slightly negative value.
In the case of σ 2 δ 2 and σ 2 δ 3 , it happens that, sometimes, s 23 is slightly greater than s 2 or s 3 , making one of the estimators of prime error variances negative. While this problem decreases as the sampling size grows, it would be convenient to filter out invalid points. In the case of α 23 = 1 (that we will assume hereafter), σ 2 > 0, so at least one of the two last estimators of the prime error variances will be positive. This means that, out of the three prime error variance estimators, at most, one can be negative. In practice, this reduces the number of cases in which no estimate is found and makes our estimates more robust.
As a filtering strategy, when presenting a map of error variances, we will discard those points for which the estimate of the error variance (not the prime error variance) is negative. Therefore, the points in the error variance maps may be different. Typically the variable with the lower error variance is more affected by this problem than the others.