Sentinel-3/FLEX Biophysical Product Confidence Using Sentinel-2 Land-Cover Spatial Distributions

The estimation of biophysical variables from remote sensing data raises important challenges in terms of the acquisition technology and its limitations. In this way, some vegetation parameters, such as chlorophyll fluorescence, require sensors with a high spectral resolution that constrains the spatial resolution while significantly increasing the subpixel land-cover heterogeneity. Precisely, this spatial variability often makes that rather different canopy structures are aggregated together, which eventually generates important deviations in the corresponding parameter quantification. In the context of the Copernicus program (and other related Earth Explorer missions), this article proposes a new statistical methodology to manage the subpixel spatial heterogeneity problem in Sentinel-3 (S3) and FLuorescence EXplorer (FLEX) by taking advantage of the higher spatial resolution of Sentinel-2 (S2). Specifically, the proposed approach first characterizes the subpixel spatial patterns of S3/FLEX using inter-sensor data from S2. Then, a multivariate analysis is conducted to model the influence of these spatial patterns in the errors of the estimated biophysical variables related to chlorophyll which are used as fluorescence proxies. Finally, these modeled distributions are employed to predict the confidence of S3/FLEX products on demand. Our experiments, conducted using multiple operational S2 and simulated S3 data products, reveal the advantages of the proposed methodology to effectively measure the confidence and expected deviations of different vegetation parameters with respect to standard regression algorithms. The source codes of this work will be available at https://github.com/rufernan/PixelS3.

remote sensing (RS) data currently play a fundamental role in many important application domains, e.g., [2]- [4]. Being one of the most important European initiatives, the Copernicus program [5] works for providing global monitoring data from space that are useful for environmental and security applications. Additionally, the European Space Agency (ESA) also develops the so-called Earth Explorer missions that are focused on addressing specific scientific challenges to advance the understanding of Earth systems and novel EO capabilities. As a result, different Sentinel concepts and Earth Explorer missions have been designed to guarantee the operational provision of RS data for dealing with current and future challenges, as well as societal needs [6].
Among the available Copernicus resources, Sentinel-2 (S2) and Sentinel-3 (S3) capitalize on important synergies, because both missions are focused on the continuous monitoring of the Earth surface using multispectral (MS) imagery. On the one hand, S2 [7] comprises two twin satellites (S2A and S2B) that carry the multispectral instrument (MSI), which provides 13 spectral bands (B01-B12) containing the wavelength region from 443 to 2190 nm, with a spatial resolution from 10 to 60 m. On the other hand, S3 [8] includes two identical satellites (S3A and S3B) that are equipped with the Ocean and Land Color Instrument (OLCI), which captures the Earth surface using 21 bands (Oa01-Oa21) in the spectral range from 390 to 1040 nm, with a spatial resolution of 300 m. In general, S2 and S3 missions provide operational products of vegetation, soil, and water cover. Nonetheless, the spatial-spectral differences between both instruments make their products more appropriate for specific RS applications [9]. The OLCI sensor is optimized to capture features related to oceans, inland waterways, and coastal areas [10], whereas the higher spatial resolution of MSI makes this sensor more effective for land-cover characterization applications [11].
At the same time, the Fluorescence EXplorer (FLEX) mission [12] is being developed as an ESA Earth Explorer mission for quantifying the photosynthetic activity and how this process affects the carbon and water cycles. Specifically, FLEX is envisaged to fly a satellite in 2024 that will be orbiting in tandem with one of the S3 constellation satellites to provide an integrated package of measurements for exploiting the terrestrial vegetation fluorescence from space. The solar-induced chlorophyll fluorescence [13] has shown to be a reliable indicator of the actual photosynthetic activity of plants, since it measures the emission of red and near-infrared light from green vegetation tissues in This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ response to the absorption of photosynthetic active radiation. However, measuring such emissions from space is challenging because they represent a small fraction of the total radiance acquired by the sensor, which eventually hinders the task of decoupling the fluorescence signal from the light reflected by the vegetation [14]. To overcome this limitation, FLEX incorporates the fluorescence imaging spectrometer (FLORIS) which has been specially designed to measure the fluorescence using an ultra-fine spectral sampling interval (SSI) of 0.1 nm (focused on the O 2 -B and the O 2 -A absorption bands, i.e., 500-780 nm spectral range) and a spatial resolution of 300 m. Despite the fact that this spectral configuration is adequate to extract the fluorescence signal from a technical perspective [15], it may be insufficient to characterize other important features related to the atmosphere while also generating some undesirable effects when dealing with the intrinsic land-cover heterogeneity [16]. In the case of the atmospheric information, FLORIS is supported by the OLCI instrument to perform an accurate atmospheric correction since FLEX will fly in tandem with one of the S3 satellites. Nonetheless, the pixel heterogeneity still remains as an open-ended issue because both sensors share the same coarse spatial resolution of 300 m.
In general, a spatial resolution of hundreds of meters per pixel may have a negative impact on the retrieval of certain biophysical parameters and also on their interpretation, based on a main consideration. The quantification of the fluorescence and other related variables at coarse spatial resolutions aggregates rather different canopy structures that make the percentages of soil, shadows, and background highly variable within the same pixel. Hence, some changes in the fluorescence amplitude can be masked by these highly variable scene components, which may eventually affect the link between the chlorophyll fluorescence sensed by the instrument and the actual vegetation physiology [17]. As a result, modeling such heterogeneity has become an important facet within the RS research community [18] and different strategies have been developed to account for the spatial heterogeneity when modeling biophysical variables [19], including multidimensional radiative transfer models [20], [21] or unmixing methods [22], [23]. However, this work explores a different research line based on exploiting inter-sensor synergies in order to model such spatial heterogeneity in S3/FLEX from a Level-4 data processing perspective.
The open availability of S2 data, together with the higher spatial resolution of the MSI sensor, provides widespread opportunities to characterize, at a subpixel level, the spatial heterogeneity of the pixels sensed by FLORIS and OLCI which can give us important insights about the reliability of the biophysical products generated within the S3/FLEX tandem mission context. On the one hand, the fluorescence signal (and other related parameters) has been proven to be highly affected by the variability of the spatial components contained in a ground pixel of the instrument [24], [25]. In this way, modeling the subpixel spatial heterogeneity of S3/FLEX throughout the S2 instrument can be very useful to relieve the gap between spatially coarse variable measurements and plant physiology [26]. On the other hand, the existing synergies among the Copernicus program and other missions (like FLEX) also motivate the development of high-level data products that take advantage of different instruments to refine, without additional costs, the operational products generated by the corresponding mission ground segments. Precisely, the growing development of advanced products for different Sentinel missions exemplifies this emerging trend [27]- [31]. Nevertheless, these and other relevant works available in the literature are mainly focused on the fusion of the data without considering the reliability of the distributed operational products among their objectives.
In this scenario, this article proposes a new methodology to study and quantify the confidence of different vegetation indicators, that serve as chlorophyll fluorescence proxies within the S3/FLEX tandem missions, using the high spatial resolution of S2 for this goal. Unlike other existing works, the proposed approach has been specifically designed for generating confidence maps of S3/FLEX vegetation products by exploiting unexplored inter-sensor relationships. More specifically, this work has a two-fold objective. On the one hand, we intend to highlight the advantages of developing confidence maps for operational data products from an inter-sensor perspective in the context of Sentinel and FLEX missions. On the other hand, we pretend to use the high spatial resolution of S2 to characterize the land-cover distribution of S3/FLEX pixels with the target of uncovering subpixel information and modeling the confidence of their corresponding biophysical indicators. With these objectives in mind, we design and develop a new statistical methodology to estimate the confidence levels of S3/FLEX fluorescence proxy variables while analyzing and identifying the situations where these confidence levels are applicable. In particular, we first characterize the subpixel spatial patterns of S3/FLEX using S2 data. Then, we model the influence of the uncovered spatial patterns and the OLCI/FLORIS product values in the errors of different vegetation indicators, via a multivariate analysis. Finally, we make use of these modeled distributions to infer the confidence and expected deviations for new test products under demand. The experimental part of the work, which considers multiple operational MSI products and their corresponding simulated OLCI counterparts, reveals the advantages of our contribution when focusing on chlorophyll estimations and also assesses the extension of the proposed methodology to other S3/FLEX biophysical parameters. In sum, the main contributions of this article can be summarized as follows.
1) We develop a new methodological framework to measure the confidence and expected deviations of S3/FLEX biophysical products by taking advantage of the higher spatial resolution of S2. 2) Focused on the chlorophyll biophysical parameter as a fluorescence proxy, we show the advantages of obtaining the confidence of such data products from an inter-sensor perspective over traditional regression algorithms. Besides, we also analyze the applicability of the proposed framework to other related parameters. 3) We build several confidence maps for simulated OLCI products to validate the performance of the proposed framework. The codes of this article will be available for reproducible research. 1 The rest of the article is organized as follows. Section II reviews some of the most important related works and highlights the novelty of this article. Section III presents the proposed methodology while detailing its main constituent parts and its adaptation to the FLEX mission. Section IV contains the experimental part of the work, conducted with operational MSI and simulated OLCI products for validation. Finally, Section V concludes the work and provides some hints at plausible future research directions.

II. RELATED WORK
The spatial heterogeneity of the Earth's surface constitutes a fundamental part of the biodiversity research [32] while also having important implications on the estimation of biophysical parameters from RS data [17]. In general, the aggregation of very different canopy structures when computing vegetation variables which are intrinsically influenced by nonlinear factors generates a negative effect on their corresponding retrieval tasks since small changes on the biophysical parameters can be easily masked by highly variable scene components [25], [33].
To relieve these limitations, the development of integral retrieval techniques to effectively model the pixel heterogeneity has become an important issue and two main trends can be identified within the RS literature [19], [34]: i) physical models and ii) empirical/statistical models. On the one hand, physical approaches are based on the so-called radiative transfer models (RTM) which describe, using the physical laws, how the remotely sensed radiation is detected by the acquisition instrument after considering certain atmospheric and surface parameters, such as atmospheric absorption, aerosol/cloud scattering and surface emissions among others [35]. In this way, these methods provide the solver of RTM equations according to a given scene representation, which may include different complexities of spatial structures and radiative transfer processes to embed the corresponding pixel heterogeneity [36]. For instance, it is the case of the lookup table approach [37] which has shown to be one of the most efficient schemes for managing complex models [38].
On the other hand, empirical/statistical methods aim at obtaining the correlation between the spectral signal captured by the sensor and the ground-truth measurements in order to model the corresponding biophysical parameters [35]. In this case, the pixel heterogeneity could be handled via some sort of disaggregation techniques, such as unmixing [22], which is able to decompose pixels into a collection of constituent spectral signatures (endmembers) and their fractional proportions (abundances). In more details, one of the most popular statistical alternatives is the nonparametric regression framework [34], where different regression algorithms are initially trained (over a specific training set) and subsequently used to make predictions under demand. For instance, Coops et al. propose in [39] using a linear regression algorithm, named the partial least squares, to estimate the content of foliage nitrogen from hyperspectral data. Addink et al. [40] also make use of another linear technique, i.e. the Ridge regression, to effectively map the biomass and leaf area index (LAI). In spite of the positive results achieved by the linear scheme, other authors advocate the use of non-linear methods instead. It is the case of Karimi et al. who propose in [41] applying the support vector regression (SVR) to estimate different types of crop variables (leaf nitrogen content, leaf chlorophyll content, etc.) from hyperspectral imagery, obtaining better performance than linear algorithms. Verrelst et al. [42] make use of the Gaussian processes regression (GPR) to accurately retrieve different biophysical parameters from simulated Sentinel-2 and Sentinel-3 images. Other authors also consider regression models based on decision trees. For instance, Im et al. [43] investigate the usability of regression trees to characterize the vegetation of hazardous waste sites.
Although these and other relevant methods have shown positive results in estimating different biophysical variables within the Copernicus program context [42], the inherent errors of heterogeneous surfaces still motivates the development of new algorithmic tools to accurately quantify the confidence of the biophysical products provided by the S3/FLEX tandem missions. Note that many of the existing works are focused on the retrieval task itself [35], but they often overlook the need of modeling the expected errors associated to the operational retrieval chains. Precisely, this work aims at relieving this operational need from a Level-4 data perspective by developing a novel statistical methodology that is able to provide performance improvements with respect to state-of-the-art machine learning regression techniques when modeling S3/FLEX errors of proxy chlorophyll fluorescence parameters.

III. PROPOSED FRAMEWORK
This section details the proposed methodology to quantify the confidence of S3/FLEX biophysical products using S2 data. Specifically, the presented framework is composed of the three following sequential steps ( Fig. 1): (A) subpixel spatial patterns extraction, where S3/FLEX subpixel variabilities are characterized according to S2 spatial information (Section III-A), (B) multivariate error distributions, where the relations among MSI spatial patterns, OLCI/FLORIS biophysical values, and their corresponding errors are uncovered via a multivariate analysis (Section III-B), and (C) product deviation estimates, where the expected deviations for test OLCI/FLORIS products are estimated (Section III-C). In addition, Section III-D provides some extra remarks about the application of the proposed framework to the future FLEX mission. The following subsections describe all these elements in detail; however, let us start by defining the notation considered in this work and the parameters of interest.
Let X = {X 1 , . . . , X N } be a collection with N MSI bottom-of-atmosphere (BOA) reflectance products and Y = {Y 1 , . . . , Y N } their corresponding OLCI/FLORIS counterparts at the same data processing level, where X i ∈ R x 1 ×x 2 ×x 3 and Y i ∈ R y 1 ×y 2 ×y 3 ∀i ∈ [1, N]. Note that x 3 and y 3 are the spectral dimensions of the MSI and OLCI/FLORIS instruments, respectively. Let R be the spatial scaling ratio between both image collections, such that x 1 = Ry 1 and x 2 = Ry 2 . Since the proposed approach pursues to exploit S2 spatial information, let S represent the maximum entropy band of X to efficiently extract MSI structural components with this simple band selection criteria. Equation (1) shows the expression required to compute the entropy measure over the jth band of the ith product (i.e., X i j). Let F be a biophysical product estimator to calculate the corresponding product feature maps as Given an external test image pair X q and Y q , the objective is to produce a quality estimate for the intrinsic error of Y F q , which is identified by ε Y F q ∈ R y 1 ×y 2 , to characterize the confidence of the biophysical indicator F in OLCI/FLORIS using the high-resolution spatial information of MSI.
Regarding the considered biophysical estimators, this work is focused on the chlorophyll biophysical variable since it is a vegetation parameter that can estimated from actual S2/S3 data and it can also be used as proxy measurements for vegetation fluorescence in the context of the S3/FLEX tandem mission.
In this work, we make use of two different non-linear vegetation indices as chlorophyll indicators for the sake of simplicity. Specifically, the normalized difference vegetation index (NDVI) [44] and the plant senescence reflectance index with near-infrared (PSRI-NIR) [45] are considered in this article. Equations (2) and (3) show their mathematical expressions, where NIR, Red, and Blue represent the corresponding instrument bands centered at near-infrared, red and blue wavelengths, respectively. On the one hand, NDVI is an effective indicator of the biomass greenness which is closely related to the fraction of solar radiation absorbed by live plants during photosynthesis. On the other hand, PSRI-NIR is focused on the ratio between carotenoid pigments and chlorophyll which becomes particularly useful to monitor vegetation physiological states. Note that these two indicators are just used to produce chlorophyll estimations that serve as vegetation fluorescence proxies in the context of this work.

A. Subpixel Spatial Patterns
The first step aims at characterizing the spatial patterns of the MSI sensor in order to represent Y according to the subpixel variability of X. Due to their significant spatial resolution differences, we initially propose conducting an inter-sensor registration procedure to relieve some possible spatial deviations between OLCI/FLORIS and MSI operational data. More specifically, we up-scale each Y i a factor of R with a bi-cubic kernel asỸ i . Then, we register both first principal components (consideringỸ i as master and X i as slave) using an affine transformation model together with the mutual information metric and the One Plus One Evolutionary Optimizer [28]. Finally, the obtained transformations are applied to X in order to make sure that each OLCI/FLORIS pixel better fits into a R × R MSI patch. Once this inter-sensor registration has been conducted, we extract patches from S (maximum entropy band of X) as P = {P 1 , . . . , P N }, where P i ∈ R M ×R 2 contains the M nonoverlapping image patches of a R × R size that are extracted from X. Note that, after this procedure, the mth pixel of the ith OLCI/FLORIS product, i.e. Y i(m) , matches with the mth patch of the ith MSI product, i.e., P i(m) . It is important to highlight that we make use of the maximum entropy band to isolate the spatial information of MSI in an efficient way. Besides, the spatial relationships between OLCI/FLORIS pixels and MSI patches are intrinsically given by their corresponding image positions since each pixel-patch pair represents the same area on the Earth surface.
With the purpose of modeling K representative types of patches (each one with a certain spatial distribution), we cluster P into K groups. In this work, we make use of the Gaussian mixture model (GMM) clustering algorithm [46] due to its proven capabilities to account for the spatial variability of RS data [47], [48]. Under GMM assumptions, each cluster is modeled as a Gaussian distribution and the corresponding cluster assignments are statistically distributed in a fuzzy fashion. In this case, our GMM can be formulated as where μ k and Σ k are the mean and co-variance parameters of the kth Gaussian component and α k represent the mixture coefficients that lie on the corresponding K-dimensional probability simplex. All these parameters completely characterize the GMM and they are collectively expressed by Θ. It is important to note that each P i(m) MSI spatial patch in P is associated to a Y i(m) OLCI/FLORIS pixel in Y since both elements cover the same region on the Earth surface. Hence, the mixture coefficients that are computed over P are inherently connected to Y according to their corresponding geolocations. Given a training set of OLCI/FLORIS products (Y) and their associated MSI patches (P), we formulate the estimation of the GMM parameters as a maximum likelihood problem via the expectation maximization (EM) algorithm [49]. As a result, each mth pixel of the ith OLCI/FLORIS product can be characterized by the cluster posterior probabilities of their corresponding inter-sensor MSI spatial patches as the posterior of P i(m) for the kth cluster. This process allows us to efficiently detect characteristic subpixel spatial distributions in OLCI/FLORIS from the higher resolution information of MSI. That is, Y pixels can be efficiently characterized at a subpixel level using the uncovered representative Gaussian components with the objective of identifying subtly spatial differences among potentially ambiguous product values in OLCI/FLORIS.

B. Multivariate Error Distributions
The objective of this step consists of modeling the influence of the uncovered subpixel spatial patterns (i.e., Y c ) and the corresponding OLCI/FLORIS product values (i.e., Y F ) in the intrinsic errors of the biophysical product estimates. In this way, these relationships can be employed for quantifying the confidence of biophysical OLCI/FLORIS products using their MSI counterpart image products. To meet this goal, we conduct a multivariate probabilistic analysis to learn a joint probability distribution of the form p(k, f, e), where k, f , and e represent the stochastic variables for MSI spatial patterns, biophysical OLCI/FLORIS values, and their corresponding product errors, respectively. For the sake of compactness, we overload k, f , and e notations to act as natural iterators in summations and as their corresponding random variables in probability arguments. Characterizing the joint probability distribution provides us complete information about the problem since it allows estimating any marginal or conditional probability distribution over the random variables of interest. More specifically, the p(k, f, e) probability distribution can be empirically approximated by normalizing the mass function of the data as follows: where H(·) represents the data frequency histogram which can be obtained using the following procedure. Initially, the empiri- Then, Y F and ε T F are discretized into B f and B e uniform bins as Y B and ε T B in order to frame the data into the following nonparametric density estimators: where · , max(·) and min(·) are the floor, maximum, and minimum functions. Finally, the multivariate frequency histogram can be computed as shown in (9) H , δ represents the Dirac delta function, and H(k, f, e) contains the number of times the spatial pattern k has a biophysical product value f with an error e. To build H, we make use of a 3-D array that contains in its x-axis Y B , in its y-axis ε T B , and in its z-axis the number of Gaussian components (i.e., K). In this way, for each xy plane of each Gaussian component, we have a grid defined by B f and B e cells. Then, for each discrete Gaussian component, we accumulate on each grid cell the mixture coefficients of the samples laying on this specific grid. As an illustration, Fig. 2 shows an example of OLCI's multivariate error distributions for PSRI-NIR when considering K = 4. For better visualization, we use a hard-clustering assignment based on the most likely Gaussian components. As it is possible to see, each cluster tends to concentrate specific biophysical product ranges (horizontal axis) and error values (vertical axis), which indicates that MSI supporting data can provide important insights about the confidence of biophysical OLCI/FLORIS products.
The joint distribution computed by (6) provides us the probability that each k, f , and e falls in any particular range of discrete product and error values for each spatial pattern. Nonetheless, we are specifically interested in modeling the resulting confidence through the p(e|f ) conditional probability, that is, the probability of error (e) for a given set of biophysical OLCI/FLORIS product values (f ). Therefore, marginalizing (6) over k, we obtain the following expressions: where (a) represents the probability of the MSI spatial patterns given the OLCI/FLORIS values and (b) the corresponding error probability distribution given both k and f stochastic variables.
In the case (a), we efficiently approximate this probability by taking advantage of the previously computed pixel-cluster assignments as Note that this expression indicates that we fix the cluster probabilities of a given discretized OLCI/FLORIS biophysical value to the posterior probability of its corresponding intersensor MSI spatial patch. Even though the absolute values of the same biophysical parameter can logically be different depending on the instrument, the spatial variability captured by MSI can be used to characterize the uncertainty in OLCI/FLORIS as a measure of the natural variability in biophysical phenomena. As a result, (11) allows a direct estimation of p(k|f ) by means of Y c . In the case of (b), this conditional probability distribution can be estimated by applying the Bayes' rule and marginalizing as follows: Inserting (11) and (12) into (10), we are able to characterize the p(e|f ) probability distribution from Y c , Y F , and T F training data.

C. Product Deviation Estimates
The last step of the proposed framework is focused on estimating the expected deviations for a query biophysical OLCI/FLORIS product. That is, given an external image pair X q and Y q , the final objective consists in estimating ε Y F q , which quantifies the deviations of the Y F q biophysical test product. In order to reach this goal, we initially conduct the inter-sensor registration procedure described in Section III-A but, in this case, between X q and Y q . Then, we extract R × R nonoverlapping patches from the S band in X q as P q , corresponding to every pixel in Y q . Subsequently, we model test cluster assignments as Y c q by computing the corresponding posterior probabilities to each one of the previously estimated Gaussian components. Note that we use the q subindex to refer to the test data and their distributions. As a result, we define the following expression: . . .
where μ k , Σ k , and α k represent the mean, co-variance matrix, and mixture coefficients of the kth Gaussian component calculated in Section III-A. Accordingly, we can also redefine p(k|f ) and p(e|f ) conditional distributions for the test data as In contrast to (10), p q (e|f ) represents that the knowledge is transferred from training to exploitation scenarios since it connects the p(e|k, f ) distribution computed over the training set with p q (k|f ) that contains the cluster probabilities for the test data. Once Y c q and p q (e|f ) have been characterized, it is necessary to discretize the test biophysical product values using the data range considered in Section III-B. Hence, we define where Y F q = F(Y q ). Finally, we propose using the expectation of the error as general confidence metric for the corresponding test products. Specifically, we estimate the final expected error for a given test biophysical value as follows: where, according to the aforementioned notation, e, k, and f are overloaded to act as natural iterators in summations (and their corresponding random variable values otherwise). Note that variable e represents the discretized error value associated to the eth error bin, and it serves to translate the conditional error probabilities into error value estimates that can be compared to the ground-truth data. In this work, we use the expected value as confidence figure; however, alternative statistics computed from p q (e|f ) (such as variance or mode) could be also used instead. Eventually, the proposed approach training and test
computations are summarized in Algorithms 1 and 2 to clarify both steps and their implementation details.

D. FLEX Peculiarities
Although the proposed methodology has been designed for managing OLCI/FLORIS products, it is important to highlight some important peculiarities of the fluorescence sensed by FLEX when compared to other chlorophyll-related variables. More specifically, the chlorophyll concentration and other vegetation products that can be computed from S3 tend to be stable in short time periods of hours or even days, which eventually makes the local observation time negligible from a practical perspective. In this way, the operational biophysical products provided by OLCI can be directly used within the proposed framework since they are able to globally reflect the seasonal vegetation state. However, the dynamic nature of the fluorescence makes it necessary to take additional considerations into account.
In the case of FLEX, the local solar time along its orbit highly affects the fluorescence captured by FLORIS due to the intraday variations in the vegetation biochemical constituents, canopy structures, and even sun-observer geometry that finally cause important physical interferences in the fluorescence retrieval process [50]. Consequently, it becomes necessary to consider some sort of auxiliary model to alleviate these fluorescence variations which are not related to the vegetation physiology  I  SID DATASET: TRAINING AND TEST MSI PRODUCT DESCRIPTIONS itself. In this work, we emphasize the need of generating, when FLEX data will be available, an empirical function based on real measurements to model how the orbit point and local time affect the sensed fluorescence. Then, each FLORIS product could be preprocessed according to its acquisition conditions to make the fluorescence values globally consistent with respect to the inter-sensor spatial information.

IV. EXPERIMENTS
In this section, the experimental part of the work is detailed, including the considered datasets (Section IV-A), the experimental configuration (Section IV-B), as well as the corresponding results and discussion (Section IV-C).

A. Datasets
In this work, we consider two different data collections of S3 OLCI and S2 MSI products, which comprise different European regions of interest, covering natural parks, urban areas, countryside, and coastal regions, among others. The first archive (Section IV-A1) is used for training and testing purposes with operational MSI imagery and synthetic OLCI products in order to relieve the lack of actual ground-truth information. The second collection (Section IV-A2) serves as external testing archive with real coupled OLCI and MSI data.

1) Synthetic Inter-Sensor Dataset (SID):
This collection includes a total of 22 real MSI images and their corresponding simulated OLCI counterparts in order to sort out the lack of actual ground-truth biophysical measurements. Table I summarizes the selected scenes where 11 images have been employed for training and the other ones for testing. Additionally, Fig. 3 displays the geographic location of the considered training (in blue) and test (in red) products. All the images are cloud-free operational MSI Level-1 C data products which have been downloaded from the Copernicus Open Access Hub platform 2 and processed using the Sentinel Application Platform (SNAP) software. In detail, they have been resampled to 20-m spatial resolution and atmospherically corrected using the Sen2Cor processor. For quantitative validation purposes, we simulate S3/FLEX and ground-truth data from these S2 images. Specifically, we employ the simulation procedure utilized in [51] to synthesize OLCI data from MSI. First, a Gaussian-like point spread function (PSF) with a full-width at half-maximum fixed  to OLCI/FLORIS spatial resolution (300 m) is defined. Then, this function is applied over the MSI data products. Finally, the resulting images are resampled to a 300-m spatial resolution using an average interpolation kernel. Additionally, we calculate the S3/FLEX ground-truth values by computing the vegetation indicators at S2 spatial resolution and then averaging those values at 300 m. Note that, the simulation of S3/FLEX data and their corresponding ground-truth information allow us to avoid inter-sensor spectral differences in the index computations due to the lack of actual field biophysical measurements for S3/FLEX over the considered regions. After all these steps, the corresponding spatial image sizes are 5490 × 5490 pixels in S2 and 366 × 366 in S3/FLEX. Note that the spatial scaling factor between the two considered types of images is R = 15×. Besides, according to the aforementioned notation, S2 and S3 training products are identified by X and Y, whereas the coupled test images are X q and Y q , respectively.

2) Real Inter-Sensor Dataset (RID):
This archive is made of four coupled S2 MSI and S3 OLCI operational data products in order to validate the proposed approach performance with actual inter-senor data. That is, the objective of using this collection is to assess the transfer of knowledge from simulated to real scenarios. Table II describes the selected OLCI/MSI test scenes that have been downloaded from the Copernicus Open Access Hub platform. Fig. 4 shows in yellow their corresponding geographic locations. Note that MSI products are fully considered whereas OLCI data have to be processed to crop the overlapping area between both sensors. In the case of MSI, we follow the same data processing pipeline described in Section IV-A1. In the case of OLCI, we make use of the SNAP software to correct the products via the radiance to reflectance processor. Besides, we reproject each OLCI image to its associated S2 tile in order to obtain the intersection between both images. Since OLCI ground-truth information is logically unavailable when considering real intersensor data, we approximate S3 ground-truth biophysical values by calculating the vegetation indicators at S2 spatial resolution and then averaging those values at 300 m. After these steps, the corresponding spatial image sizes also are 5490 × 5490 pixels in S2 and 366 × 366 in S3. Note that, unlike SID, this dataset may contain some inter-sensor spatial misalignments which will be taken into account by the inter-sensor registration procedure integrated within the proposed framework.

B. Experimental Settings
The experiments conducted in this work pursue to validate the proposed approach performance to estimate the confidence of S3/FLEX products using S2 spatial information. Since the presented method deals with the problem from an inter-sensor perspective, regression algorithms can also be applied to generate such confidence estimations in a straightforward way. Consequently, the experimental part of the work is aimed at providing evidence that the proposed approach is able to improve general regression techniques when modeling intrinsic S3/FLEX biophysical product errors. To achieve this goal, the presented framework is compared to five different regression algorithms, i.e., linear regression [52], ridge regression [52], SVR with radial basis function [53], Gaussian process regression with squared exponential [54], and binary decision tree for regression [55]. Additionally, we consider the NDVI and PSRI-NIR vegetation indicators to evaluate the considered methods. Specifically, we make use of these two indexes because both can be computed from real S2/S3 imaging data and they also allow us to approximate the chlorophyll content as proxy measurements for the vegetation fluorescence framed in the S3/FLEX tandem mission context. Under this scheme, we conduct the two following experiments.
1) Experiment 1: In this experiment, we employ the SID collection for fully training and testing the proposed approach and all the considered competitors. More specifically, all the regression techniques have been trained using the simulated biophysical products of S3/FLEX (Y F ) and their corresponding spatial MSI characterizations (Y c ) to predict the associated ground-truth errors, that is, ε T F = |Y F − T F |. Note that these ground-truth errors are also required for training any supervised model (likewise the proposed approach and regression models) or even validating unsupervised ones. Additionally, empirical ground-truth measurements acquired over specific regions of interest could be used for computing these ground-truth errors from real data instead of simulated. In all cases, we make use of the MATLAB R2019a implementations with automatic scale, data standardization, and the default settings for the other available parameters. For the sake of efficiency, we train the regressors with a random subset of 5 × 10 4 samples using a 1% percentile threshold for outlier removal. In addition, we consider an extended subset of 5 × 10 5 samples for computing the spatial characterizations used by the proposed approach and the considered competitors. Once the training phase is completed, each test product is raised as an external image pair (Y F q and Y c q ) to estimate its expected deviations (ε Y F q ). Then, we quantify the effectiveness of the methods by computing the mean squared error (MSE) between the ground-truth deviations ε T F q and ε Y F q . Specifically, we employ this figure of merit due to its quadratic error computation that penalizes highly differing predictions. To generate the ground-truth biophysical maps in S3/FLEX, i.e., T F and T F q , we downscale X F and X F q using an average R × R kernel. In the case of the proposed approach, it is important to highlight that the two first stages are considered for training and the last one for testing. Besides, we also analyze the following parameter settings B f , B e = [64 128 256, 512] and K = [4,8,12], for the number of discrete product values and the number of spatial patterns in the clustering, respectively.
2) Experiment 2: This experiment makes use of the RID collection as an external test collection to assess the proposed approach performance with real inter-senor data. In more details, we evaluate how the models trained in Section IV-B1 perform over RID for investigating the transfer of knowledge from simulated to real scenarios. For the sake of simplicity, we use the same testing experimental settings described in Section IV-B1, but only considering the best number of clusters. Note that, in this experiment, there may appear some inter-sensor spatial misalignments that the proposed approach relieves with the incorporated inter-sensor registration procedure. In order to conduct an experimental comparison as fair as possible, we also use this registration mechanism (Section III-C) for all the tested competitors.

C. Results and Discussion
1) Experiment 1: Tables III-V provide the quantitative assessment of the PSRI-NIR error estimates when considering different number of clusters: 4, 8, and 12, respectively. Specifically, each table shows the MSE metric between the estimated errors and their corresponding ground-truth values of the PSRI-NIR product to evaluate the tested methods from a quantitative perspective. In rows, the results obtained for test products are shown, being the last row the corresponding average. In columns, all the considered regression algorithms together with the proposed approach are displayed. Note that, in the case of   the proposed approach, different values for the B f,e parameters are tested in the last four columns. It should also be mentioned that the best metric result for each test product is highlighted in italic font to facilitate analysis. Analogously, Tables VI-VIII present the MSE results of the NDVI vegetation index using the same format.
In addition to the quantitative evaluation provided by the MSE metric, some visual results of the corresponding product error estimates are also presented to evaluate the methods from a qualitative perspective. In particular, Fig. 5 shows the estimated product errors for the 32ULB test tile when considering the PSRI-NIR biophysical parameter and K = 8. As it can be observed, Fig. 5(a) shows the simulated OLCI/FLORIS PSRI-NIR product, Fig. 5(b) contains the corresponding ground-truth errors, and from Fig. 5(c)-5(h), the error maps generated by   all the considered methods are displayed, i.e., LIN, RID, SVR, GPR, BRT, and the proposed approach (with B f,e = 128). Note that the color palette used for the biophysical product visualization has been fitted into the PSRI-NIR parameter range, i.e., [−1,1]. Additionally, all the error maps make use of the same color codes, ranged between the minimum and maximum global values in order to easy the comparison. In all the visual results, a small image detail has been also magnified to better show even small performance deviations among the methods with homogeneous and heterogeneous pixel areas. Following the same configuration, Fig. 6 presents the visual results obtained over the 32TQR test product with K = 8 when considering the NDVI biophysical parameter. In this case, we show the proposed approach estimated errors when B f,e = 512 since it provides the best quantitative performance.
According to the results presented in Tables III-V, there are several aspects which deserve to be mentioned regarding to the PSRI-NIR error estimates. The first remarkable point is related to the general performance of the considered methods. Note that the proposed approach obtains the best MSE-based quantitative result followed by GPR, RID, SVR, LIN, and BRT in the corresponding performance order. In the case of the proposed approach, the best performance is achieved when considering B f,e = 128; nonetheless, the existing differences are certainly small, which indicates its robustness to the discretization parameter. When focusing on the other methods, the proposed approach is able to significantly reduce the average MSE index by a 13% with respect to GPR, 19% for RID/SVR, 23% for LIN, and a 47% for BRT. Another important aspect can be observed when analyzing the experiments in more details based on the number of clusters, i.e., K = [4,8,12]. Specifically, the obtained results show that considering a higher number of clusters seems to provide some performance advantages when estimating the errors; however, such improvements tend to be unsteady while being reduced. As an example, RID, SVR, and BRT obtain the best results with K = 8, whereas LIN, GPR, and the proposed approach generate the most accurate error estimations with K = 12. In many of the cases, the performance variations are generally small, which indicates that from a certain number of clusters, they rather affect the final error estimations with the exception of GPR. From an intuitive perspective, this behavior could be justified by the type of clustering we carry out in this work. Note that we use GMM with soft-clustering assignments and, in this situation, using many mixture components may lead to a higher variance in the resulting estimates [56].
When considering the NDVI error estimations (Tables VI-VIII), a similar trend can be observed but with some differences. Although the global average results are quite similar to the PSRI-NIR case, we can observe that the proposed approach improvements with respect to RID, SVR, and GPR are reduced for the NDVI parameter. Specifically, the presented method is able to decrease the average MSE metric by a 49% for BRT, a 21% for LIN, a 12% for RID, a 12% for RID, 4% for SVR, and 3% for GPR. These variations can be generated by the different nature of the NDVI index since the corresponding ground-truth errors are more than twice smaller than those of PSRI-NIR. Precisely, an evidence of this fact can be found in the discretization parameters. According to the PSRI-NIR results, the proposed approach consistently achieves the best performance when considering B f,e = 128. However, the optimal value for NDVI is B f,e = 512, which reveals that the corresponding multivariate error distribution [i.e., H(k, f, e)] requires a higher precision for the discretization because the objective errors are significantly smaller. Regarding the number of clusters, it is also possible to see that the NDVI index has a slightly different trend. In this case, the best results are always achieved with K = 12 but the performance differences with respect to K = 8 are very small (except for LIN).
With all these considerations in mind, the proposed approach achieves the best global quantitative error estimation followed by GPR, SVR, RID, LIN, and BRT. However, in order to complete the experimental comparison, an additional qualitative analysis has been conducted with a particular focus on those methods with similar performances. From the visual results displayed in Nonetheless, the high amount of output noise often makes this method overestimate product errors. When comparing the details of Fig. 5(b) and (g), we can see that BRT is able to correctly predict many errors along the river shape but it also introduces important inaccuracies that eventually penalize the final result.   that also make the proposed approach generally better from a qualitative perspective. It is worthwhile noting that the proposed approach is able to estimate the errors along the river in a more accurate way while generating the same limited amount of noise. In the case of the NDVI index (Fig. 6) , we can see that these methods obtain the two best visual results since they produce the most similar outputs to the corresponding ground-truth error map [ Fig. 6(b)]. However, the proposed approach results are certainly the best ones since the errors predicted in the coastal area definitely have a more accurate shape than the GPR ones.

D. Experiment 2
Tables IX and X contain the quantitative assessment of the PSRI-NIR and NDVI error estimates when considering K = 12. Likewise, in the previous experiment, test data products are shown in rows whereas columns display the considered methods. According to the reported results, it is possible to make some important observations. From a general perspective, we can see a global quantitative performance decrease produced by the use of actual OLCI products for testing the biophysical indices. Note that, we use the computation of PSRI-NIR and NDVI over MSI as OLCI's ground-truth (since we do not have actual ground-truth information for OLCI). However, the existing spectral differences between both sensors make these vegetation indicators do not obtain the same result across the instruments, which logically affects the quantitative evaluation. When analyzing the results in more details, we can observe that the proposed approach obtains the best average performance followed by BRT, GPR, SVR, LIN, and RID. Despite the fact that some of the considered regression algorithms show some deviations with respect to the use of synthetic test data, the proposed approach is able to obtain the best quantitative results, which reveals its ability to manage real inter-sensor data.
In general, the conducted experiments reveal that the proposed approach provides competitive advantages when estimating OLCI biophysical product errors from a Level-4 data processing perspective. Precisely, these improvements are based on two main aspects: the inter-sensor scheme and its probabilistic nature. On the one hand, the presented inter-sensor scheme allows taking advantage of the higher spatial resolution of S2 to uncover S3/FLEX subpixel information that becomes very useful to characterize the biophysical product errors generated by the lack of spatial resolution. Under the context of the Copernicus program and other missions like FLEX, inter-satellite synergies can play a fundamental role to relieve sensing limitations by exploiting the strengths of each individual instrument. In this sense, the proposed inter-sensor scheme takes and implements these ideas to fulfill the inherent requirements of estimating OLCI/FLORIS biophysical product errors in a novel way using inter-sensor MSI data. On the other hand, the proposed methodology is able to produce more accurate uncertainty estimates than traditional regression algorithms because it deals with the problem from a multivariate analysis perspective which has been specially designed to exploit S3/FLEX and S2 properties. That is, jointly modeling biophysical product values, errors, and subpixel spatial patterns as empirical distributions allows us to better relate all these components based on the richer spatial information provided by S2 while also accounting for the statistical significance of the data that eventually becomes essential to reduce noisy error predictions as displayed in the qualitative results.
Traditional regression algorithms pursue to project sample points from an initial feature space made of the concatenation of MSI spatial patterns (k) and OLCI/FLORIS biophysical values (f ) onto a target domain with their corresponding errors (e). Although each method logically follows its own intuitions for conducting such projection, all of them share a common aspect: the inter-sensor information provided by MSI and the intrasensor data of OLCI/FLORIS are integrated into a single characterization by means of concatenation. That is, both information sources are jointly considered as a single observable entity. However, the proposed framework provides a more powerful scheme since two different random variables are considered in this regard. More specifically, the proposed approach defines a probabilistic model where the corresponding errors (e) depend on k and f random variables in a way that the conditional probability distribution p(e|k, f ) is able to control the balance between MSI spatial patterns (k) and OLCI/FLORIS biophysical values (f ). In other words, both k and f are considered observable variables that affect the error but this influence is controlled for each observed pair, which eventually becomes more suitable for managing the inter-sensor spatial variability than using the standard feature concatenation of regular regression algorithms.

V. CONCLUSION
This work has proposed a novel statistical approach to study and quantify the variability of the biophysical products provided by the S3/FLEX tandem missions by taking advantage of the higher spatial resolution of S2. Specifically, OLCI/FLORIS subpixel spatial patterns are initially characterized using inter-sensor MSI data. Then, the influence of these patterns and biophysical values in the operational product errors is modeled using a multivariate analysis perspective. Finally, the uncovered distributions are employed to infer the confidence level and expected deviations for new test products. A comprehensive experimental comparison, including multiple operational MSI products and their simulated OLCI counterparts, has been conducted to validate the proposed approach with respect to several regression methods using different biophysical parameters.
One of the first conclusions that arises from this work is the importance of the inter-sensor data within the S3/FLEX context and how complementary S2 information can support managing the spatial heterogeneity problem of biophysical parameters. That is, this work shows that considering a Level-4 data processing perspective allows for the exploitation of inter-sensor synergies that may contribute to relieve some important sensing limitations without additional costs.
Another relevant aspect is related to the operational reliability of S3/FLEX biophysical parameters. The presented methodology has been designed to improve the certainty of operational S3/FLEX products by using a new statistical approach that is able to model the relationships between biophysical errors and inter-sensor patterns. In this regard, the proposed approach is able to outperform traditional regression algorithms, with both synthetic and real imagery, by considering two independent inter-sensor observable variables unlike standard feature concatenation. Precisely, such improvements may become very useful in operational environments to produce more relevant S3/FLEX biophysical products for different downstream RS applications. In our future work, we will investigate the extension of the proposed framework to manage multimodal data from other missions while developing additional models to better characterize subpixel information by means of spatial-spectral patterns. The study of RS data processing techniques to deal with inter-sensor cloud contamination will be also examined as future work. He has been a Visiting Researcher with the University of Bristol, Bristol, U.K., the University of Cáceres, Cáceres, Spain, and the Technische Universität Berlin, Berlin, Germany. He is a Postdoctoral Researcher with the Computer Vision Group, Universitat Jaume I, where he is also a member of the Institute of New Imaging Technologies. His research interests include multimedia retrieval, spatio-spectral image analysis, pattern recognition techniques applied to image processing, and remote sensing.
Dr. Fernandez-Beltran is also a member of the Spanish Association for Pattern Recognition and Image Analysis (AERFAI), which is part of the International Association for Pattern Recognition (IAPR). He has received the Outstanding Ph.D. Dissertation Award from Universitat Jaume I, in 2017.
Filiberto Pla received the B.Sc. and Ph.D. degrees in physics from the Universitat de València, Spain, in 1989 and 1993, respectively.
He is now a Full Professor with the Departament de Llenguatges i Sistemes Informàtics, University Jaume I, Castellón, Spain. He has been Visiting Scientist with Silsoe Research Institute, the University of Surrey, University of Bristol, U.K.; CEMAGREF, France; the University of Genoa, Italy; Instituto Superior Técnico of Lisbon, Portugal; the Swiss Federal Institute of Technology ETH-Zurich, the IDIAP Research Institute, Switzerland; the Technical University of Delft, The Netherlands; and the Mid Sweden University, Sweden. He has been the Director of the Institute of New Imaging Technologies, University Jaume I. His research interests include color and spectral image analysis, visual motion analysis, 3D image visualization, and pattern recognition techniques applied to image processing.
Dr. Pla is member of the Spanish Association for Pattern Recognition and Image Analysis (AERFAI), which is a member of the International Association for Pattern Recognition (IAPR).