The interests of the scientific community working on the Soil Moisture and Ocean Salinity (SMOS) ocean salinity level 2 processor definition are currently focused on improving the performance of the retrieval algorithm, which is based on an iterative procedure where a cost function relating models, measurements, and auxiliary data is minimized. For this reason, most of the effort is currently focused on the analysis and the optimization of the cost function.
Within this framework, this study represents a contribution to the assessment of one of the pending issues in the definition of the cost function: the optimal weight to be given to the radiometric measurements with respect to the weight given to the background geophysical terms.
A whole month of brightness temperature acquisitions have been simulated by means of the SMOS-End-to-End Performance Simulator. The level 2 retrieval has been performed using the Universitat Politècnica de Catalunya (UPC) level 2 processor simulator using four different configurations, namely, the direct covariance matrices, the two cost functions currently described in the SMOS literature, and, finally, a new weight (the so-called effective number of measurement).
Results show that not even the proposed weight properly drives the minimization, and that the current cost function has to be modified in order to avoid the introduction of artifacts in the retrieval procedure. The calculation of the brightness temperature misfit covariance matrices reveals the presence of very complex patterns, and the inclusion of those in the cost function strongly modifies the retrieval performance. Worse but more Gaussian results are obtained, pointing out the need for a more accurate modeling of the correlation between brightness temperature misfits, in order to ensure a proper balancing with the relative weights to be given to the geophysical terms.
a. The SMOS mission
In May 1999 the European Space Agency (ESA) approved the Soil Moisture and Ocean Salinity (SMOS) Mission as the second of its Living Planet Programme Earth Explorer Opportunity Missions to provide global and frequent soil moisture and sea surface salinity (SSS) maps.
SMOS was launched on 2 November 2009, and after the first calibration and checkout period (the so-called “commissioning phase”), SSS level 3 products will be distributed; the expected accuracy is 0.1–0.4 psu over 100 × 100–200 × 200 km2 in 30–10 days, respectively (Font et al. 2004).
The single payload embarked on SMOS is the Microwave Imaging Radiometer by Aperture Synthesis (MIRAS; McMullan et al. 2008); it is a 2D interferometric radiometer operating at the protected L band with a nominal frequency of 1413.5 MHz and a bandwidth of 27 MHz. It consists of three deployable arms connected to a central hub (8-m-diameter radiometer when completely deployed). The arms are equally spaced with an angular separation of 120°. Each arm encompasses three segments, each one containing six L-band radiometers [Lightweight Cost-Effective Frontend (LICEF)]; four more radiometers are placed in the central hub, for a total of 66 radiometers. In addition, there are three noise injection radiometers located in the central hub, each of which consists of two LICEF receivers coupled to a single antenna. The total number of elements is, therefore, 69 antennas and 72 receivers.
b. The measurements acquisition
SMOS is an interferometric radiometer. The basic concept of interferometric radiometry is to synthesize a large aperture using a number of small antennas. The output voltage of each pair of antennas [e.g., antenna 1 and 2, located at (x1, y1) and (x2, y2)] are cross correlated to obtain the “visibility samples” as expressed by the following equation:
where u and υ are the spatial frequencies of the visibility sample , kB is the Boltzmann constant , B1 and B2 the receivers’ noise bandwidth, G1 and G2 the available power gains, and b1(t) and b2(t) are the signals measured by the elements 1 and 2, respectively.
The complete set of visibility samples is called the visibility map, and it is approximately the Fourier transform of the brightness temperature distribution of the scene. To invert this process either the inverse Fourier transform can be applied as a first approximation (Camps et al. 1997) or a more sophisticated G-matrix inversion (Camps et al. 2008; Anterrieu and Camps 2008) can be used. The major advantage of interferometric radiometry is the multiangular measurement capability: the output of an interferometric radiometer is, in fact, an image; this permits several views, under different incidence angles, of the same point on the earth before it exits the field of view (FOV).
According to the MIRAS instrument design, the distance (d) between antennas does not satisfy the Nyquist criterion (Camps et al. 1997), and part of the FOV is affected by aliasing. The SMOS FOV always contains both the earth and the sky; because the sky is a very stable and well-known target, both its direct contribution and the alias that it induces can be estimated and removed by the visibility map. The resulting so-called Extended Alias-Free (EAF)-FOV has the shape of a distorted hexagon. Figure 1 shows the EAF-FOV and the variation of the incidence angle (dashed line) and the spatial resolution inside it (dash–dot). A point in the boresight (the center of the swath for the case of SMOS) of the satellite is observed approximately 150 times (75 for each polarization) under an incidence angle ranging from 0° to 65°, and with a spatial resolution from 30 to 100 km.
c. SSS retrieval in SMOS
The SMOS level 2 (from brightness temperatures to SSS) retrieval algorithm has been defined according to a Bayesian approach: it embodies prior information to ease the retrieval. Assuming normal statistics on both the a priori information and the observations, the general maximum likelihood estimation (MLE) reduces to a least squares problem, the solution of which can be found through the minimization of a “cost function” (χ2) expressed by
where is the misfit in the observations (measurement minus model) and is the misfit covariance matrix. Until a proper estimation of is obtained in the official ocean salinity level 2 processor, it is defined as being diagonal and the misfits are considered completely uncorrelated; this consideration is equivalent to writing Eq. (2) as
where is the nth element of , which is a function of the sea surface temperature, salinity, and roughness; and is the radiometric noise for the nth observation.
Previous studies (e.g. Gabarró et al. 2009) showed that defining the optimal cost function is not straightforward and that auxiliary external information, in particular, wind speed (U10), sea surface temperature (SST), and possibly modeled or climatological SSS must be added to Eq. (3).
The following two different cost functions are currently used within the SMOS community (Zine et al. 2008):
In both formulations the cost function is composed of the following two main contributions (or information providers):
the first term is representative of MIRAS measurements, which is a function of the original (true) geophysical parameters , and modeled observables, which is a function of the parameters that are going to be retrieved , weighted by the radiometric noise of the nth observation as in Eq. (3); and
- the constraints for the auxiliary SSS, SST, and 10-m-height wind speed U10 (as sea surface roughness descriptor) are the second, third, and fourth terms, respectively, ; these are weighted by the inverse of the variance of the misfit existing considering the corresponding auxiliary field with respect to the original one, as defined in Eq. (6):
In the real SMOS case (not the simulation), Porig is not known and allows weighting of the a priori information with the value of the geophysical parameter. The value of is representative of the reliability of this information: large indicates that the estimation is not reliable, leading to a very small weight within the total χ2 minimization, and vice versa.
The difference between the two formulations lies in the factor 1/Nobs weighting of the observables term in Eq. (5). Actually, Eqs. (4) and (5) represent two extreme cases; in the first each misfit is assumed to provide the maximum information content, and their contributions are thus summed up (once squared and normalized by the radiometric noise) to construct the brightness temperatures term in the cost function [Eq. (4)]. The second option [Eq. (5)] is appropriate for the case of completely redundant misfit samples: the average contribution is used to define the cost function.
d. Expected correlation in the brightness temperatures
Because of MIRAS’ characteristics, some correlation is expected between the brightness temperature errors of different grid points within the same snapshot and among consecutive snapshots. Any imaging radiometer, in fact, is affected by the following three types of noise (Font et al. 2008):
the radiometric resolution (Randa et al. 2008) (ΔT) (the temporal standard deviation of the zero mean random error resulting from the finite integration time);
the radiometric bias (the spatial average of all of the systematic errors);
and the radiometric accuracy [the spatial standard deviation of all of the systematic errors; see Torres et al. (2005)].
The first type of noise is random within the same snapshot as well as from snapshot to snapshot. The second and third types are random from pixel to pixel within the same snapshot, but they are systematic from one snapshot to another, and are responsible for the above-mentioned correlation. In addition, concerning SMOS, the following two other sources of spatial correlation can be identified:
the finite spatial resolution of the instrument (Fig. 1), which is larger than the icosahedral Snyder equal area hexagonal grid of aperture 4 and resolution 9 (ISEA 4H9; see Snyder 1992; Suess et al. 2004) in which level 2 data will be projected (SMOS 50-km-average resolution pixels are oversampled to the 15-km ISEA 4H9 grid size).
Thus, to summarize, correlation is present in the SMOS processing chain at level 0 (for radiometric bias and accuracy) and at level 1 (for image reconstruction and projection of the brightness temperatures), as shown in Fig. 2.
Considering the earth’s reference frame instead of the satellite’s reference frame, each grid point “observes” the satellite and samples its antenna pattern with a frequency related to the time between snapshots (whose upper limit is fixed by the ISEA grid spacing). As it is explained, the SMOS synthetic antenna pattern presents a correlation length that is larger than this sampling frequency, inducing correlation among the errors on the various measurements of the same grid point.
As remarked in section 1c, according to the SMOS level 2 retrieval procedure, the misfit in the SMOS measured brightness temperatures is assumed to be uncorrelated, which is equivalent to considering the matrix in Eq. (2) as being diagonal, and thus leading to Eqs. (3) and (4). Taking into account Eq. (4), the presence of correlation between different misfits results in a loss of the information provided by the observables with respect to the background terms.
Two different approaches have been followed so far concerning salinity retrieval from MIRAS brightness temperatures: misfits can be, once squared and normalized, summed up [Eq. (4); see Zine et al. (2008)] or averaged [Eq. (5); see Camps et al. (2005) and Talone et al. (2007)]. However, the correlation induced by the instrument generates an intermediate and more complex situation between Eqs. (4) and (5).
Aiming at evaluating the impact of the correlation among the misfits on the SMOS-measured brightness temperatures, a whole month of overpasses has been simulated.
The simulation scenario is presented in section 2. The level of correlation of the measurement errors, and thus the weight to be given to the observables term in the cost function, is assessed in section 3. In this section, the covariance matrices are estimated, and a new weight regarded as the “effective number of observations” is introduced. The comparison of the retrieval results using the four formulations is considered in section 4, and, finally, the main conclusions of this work are summarized in section 5.
2. Simulation scenario
Because at the time of this study SMOS output was not yet fully calibrated, SMOS-like brightness temperatures were simulated using the SMOS End-to-End Performance Simulator (SEPS; see Camps et al. 2003; SEPS 2006a) in its full mode (including copolarized and cross-polarized measured antenna patterns for each antenna, all instrument errors, and G-matrix image reconstruction). To model sea surface emission, the Klein and Swift model for the seawater dielectric constant (Klein and Swift 1977) and the linear fit to Hollinger’s measurements (Hollinger 1971) for the wind speed contribution to brightness temperature have been used. Because the objective of this study is the estimation of the correlation induced by the instrument, in order to avoid further contributions from other sources the radiometric sensitivity has been set to zero in the simulations, increasing the integration time to very large values (Randa et al. 2008), that is, . The radiometric sensitivity is, in fact, according to its definition, already taken into account by the term in Eqs. (3)–(5). Sixty-four ascending and descending overpasses have been simulated during the month of March 2007 (SEPS time) consisting, on average, of more than 200 snapshots each. As mentioned, the measurement acquisition has been simulated using the measured MIRAS antenna patterns, instrument drifts, and current G-matrix inversion algorithm. Simulations output has been projected onto the ISEA grid as real SMOS data.
Concerning the geophysical parameters, the following two databases are defined:
original data (used to feed SEPS and generate the brightness temperatures): daily outputs of a 0.5° configuration of the Nucleus for European Modelling of the Ocean (NEMO)–Océan Parallélisé (OPA) ocean model (Madec 2008; Mourre et al. 2008) used as original SSSorig and SSTorig, while U10orig fields come from 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40; Uppala et al. 2005); and
auxiliary data [used in the level 2 cost function; see Eqs. (4) and (5)]: SSSaux and SSTaux come from the Levitus climatology (Levitus 1998), and U10aux are extracted from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (Kalnay et al. 1996);
Simulations have been carried out following these steps:
SEPS-generated brightness temperatures (Fig. 3a) have been masked to eliminate the transition areas at the beginning and at the end of the sequence, and the remaining grid points are shown in Fig. 3b. As an example, the number of observations for one of the ascending overpasses is shown in Fig. 3c as a function of the distance to the ground track.
Selected brightness temperatures have been compared to the ones resulting from running only the forward model (using the same geophysical and orbital parameters); the difference between them is the instrument-induced radiometric error (radiometric bias plus radiometric accuracy, because radiometric sensitivity has been set to zero).
The calculated differences have been sorted and grouped by a number of observations. For each one of the bins, the covariance matrix has been computed as follows: for each one of the nth observation samples, the Galton–Pearson’s correlation coefficient (Rodgers and Nicewander 1988) has been calculated between all the possible pairs to construct the covariance matrices [estimate of , in Eq. (2)]. Bins with a number of samples smaller than the number of observations have not been taken into account because they do not provide representative results.
Sea surface salinity retrievals for all the 64 overpasses have been performed using the estimated covariance matrices as expressed by Eq. (7):
Results have been considered as a master case to be compared with the various approximations of the cost function.
At this stage, aiming at adapting the current cost functions to the characteristics of the misfits’ covariance matrix, a new weight is defined. To do so, an analysis of the estimated covariance matrices has been carried out. The eigenvector decomposition has been applied to the inverse covariance matrices and the number of eigenvectors describing 99% of the variance has been defined as the effective number of measurements Neff.
To test the impact of introducing the Neff in the cost function, the SSS for the simulated scenario has been retrieved by means of the level 2 processor (Talone et al. 2007), the cost function in Eq. (8):
First, the covariance matrices have been calculated and analyzed. Figure 4 presents an example of a covariance matrix (78 pairs of TH–TV observations), the matrix is plotted for (a) horizontal polarization (H pol), (b) vertical polarization (V pol), and (c) the first Stokes parameter in brightness temperature (TI) (Randa et al. 2008); the color scale is in decibels to enlarge the observed dynamic range. According to Butora and Camps (2003), the error in the radiometric measurements is correlated; correlation clusters are evident in Fig. 4 and characterized by a very complex pattern.
As explained in section 3, the eigenvector decomposition has been applied to the inverse covariance matrices; in Fig. 5 an example of the eigenvalue spectrum is shown for the case using the first Stokes parameter in brightness temperatures, with (a) 22 and (b) 78 pairs of observations, respectively. As can be observed, the trends are very different, in the first case (22 observation pairs) the spectrum is almost constant and sharply decreases in the last three eigenvectors; in the second case, instead, the decrease is more uniform. This change of regime can be also noticed in the trend of the effective number of observations (Neff), derived by the analysis of the covariance matrices as described in section 3. Results are shown in Fig. 6 for (a) H pol, (b) V pol, and (c) the first Stokes parameter in brightness temperature (TI). The number of (TH, TV) pairs (approximately half the number of observations) is represented in the abscissas, and the ordinates shows the ratio between the effective number of observations Neff and the total Nobs. The ratio Neff/Nobs is shown as a density plot, with color showing the occurrence of any particular pair Nobs–Neff/Nobs simulated along the whole month; the solid line is the linear fit of Neff. Figure 6d is the normalized histogram of the number of observations (the sum of all of the bins is equal to 1).
Even though some differences can be noticed between TH, TV, and TI (H gives higher results), the trend is very similar, and two different regimes can be observed. For Nobs ∈ (1, 30], Neff changes from being equal to Nobs to the asymptotic value of 0.82Nobs for H pol, 0.76Nobs for V pol, and 0.79Nobs for TI; when is Nobs ≥ 31 (with a distance to ground track <350 km) the ratio Neff/Nobs remains almost constant. The steep change in the regime around 30 observations is basically due to the lack of samples between 30 and 70 pairs of observations. This behavior is apparent from Figs. 6d and 3c: in the former figure the normalized histogram of the number of observations is shown, whereas in the last figure the number of observations is presented as a function of the cross-track distance. The very steep increase of Nobs shown in Fig. 3c is related to the along-track dimension of the SMOS FOV, shown in Fig. 1, and is the cause of the lack of estimates of Neff in Figs. 6–c, as explained in section 2. In fact, when only a few measurements are available for a certain Nobs, the covariance matrix is not calculated because it is not representative. This change of regime suggests an objective way of defining the useful swath width of an SMOS overpass at approximately 700 km, which is a bit larger than the official Q swath or the “narrow swath” (631 and 640 km; see SEPS 2006b; Barré et al. 2008). Another feature to be highlighted in Figs. 6a–c is the banded structure of the ratio Neff/Nobs for low values of Nobs, which was not expected, and is probably due to the specific shape of the SMOS FOV and the consequent distribution of Nobs in the cross-track dimension.
Figure 7 shows W/Nobs as a function of the total number of observations; the Neff calculated for the first Stokes parameter in the brightness temperatures is marked with the solid line, while the dashed and dash–dot lines stand, respectively, for the cases in Eqs. (4) and (5),
where the symbol ( )′ indicates the fitting result and not the calculated Neff/Nobs.
As can be observed, Neff takes intermediate values between 1 and Nobs, as expected.
The four configurations [using the cost functions defined in Eq. (7): (master), Eq. (4): (W = Nobs), Eq. (5): (W = 1), and Eq. (8): (W = Neff)] have been tested and results have been compared. SSS has been retrieved using the brightness temperature resulting from simulating the same scenarios, but now the SMOS nominal value for the integration time (τ = 158 ms) is applied.
Error statistics considering only the grid points fully observed for all 64 overpasses are summarized in Table 1 through its mean value (μ), standard deviation (σ), rms (defined as , with N being the total number of grid points taken into account), and the X2 factor.
The latter is defined as the quadratic sum of the difference between the retrieved SSS error–normalized histogram (observable) and the pdf of a normal distribution with the same mean and variance (model), weighted by the uncertainty associated with each observation, as expressed in Eq. (11) (Barlow 1989),
The normalized histograms of the SSS retrieval error (the sum of all of the bins is equal to 1) using Eqs. (7), (4), (5), and (8) are shown in Figs. 8a–d, respectively. In order to not alter the results, σi in Eq. (11) has been considered constant and equal to 1. To calculate X2, both the SSS error and the Gaussian pdf have been quantisized in 0.1-psu bins (which is the expected resolution of SMOS at level 2 in one overpass); moreover the sum, which should be calculated in the interval (−∞, ∞), has been computed only for the interval (−10 psu, 10 psu).
According to Table 1, retrieval results are only slightly affected by the change of W, because for all of the configurations the rms error is constant at 2.39 psu, with very high X2 (~5.2 − 5.4). The difference is, instead, noticeable if it is compared with the case of directly using the covariance matrices in the retrieval. In this case, the rms error is equal to 3.78 psu [mostly resulting from the increase of the error standard deviation σ (3.75 psu)]; on the other hand, error statistics are much more Gaussian presenting a X2 = 0.89. Comparing Figs. 8a–d it is evident that results are better when using whichever of the Eqs. (4), (5), or (8); nevertheless, the low Gaussianity of the error is also manifest, indicating that some artifacts are introduced in the retrieval procedure resulting from the high weight given to the constraints of the auxiliary parameters. Previous studies (Sabia et al. 2010; Gabarró et al. 2009) remarked on the necessity of correctly balancing the different terms in the cost function. The result of the minimization of the cost function, in fact, strongly depends on the relative weight given to each of the factors of χ2. In particular, because all of the elements of the covariance matrices, calculated as explained in section 3, are positive when Eqs. (4), (5), or (8) are used, the contribution of the constraints for SST and U10 to the total cost function is much larger than the contribution of the measured brightness temperature , if compared to the case using Eq. (7). The consequence is a good, but fake, retrieval, when the retrieved parameters drift toward the reference parameters.
According to the results of this study, the inclusion of the brightness temperature misfit covariance matrices [Eq. (7)] gives a very different result with respect to the case of using the approximated cost functions [Eqs. (4), (5), and (8)] and should be taken into account in the choice of the relative weights, which should be updated.
To improve the characterization of the cost function used in the SMOS ocean salinity level 2 processor, the correlation between measurement misfits has been analyzed using simulated data. Correlation is expected, due to the intrinsic nature of any imaging radiometer, to the possible structures induced by the image reconstruction algorithm, and, finally, to the projection of the brightness temperatures onto the ISEA grid.
To assess this point, one complete month of overpasses (64 in total) in the North Atlantic Ocean have been simulated using the SMOS End-to-End Performance Simulator (SEPS) in its full mode. The SMOS level 2 processor simulator (SMOS-L2PS) has been used to retrieve SSS from the brightness temperatures calculated by SEPS. As geophysical input parameters (original and auxiliary data), a North Atlantic configuration of the NEMO-OPA ocean model and the Levitus climatology have been used for SSS and SST, while ERA-40 and NCEP–NCAR products have been chosen for U10.
The SEPS-simulated brightness temperatures have been compared to the ones obtained by directly forwarding a brightness temperature model to estimate the correlation of the radiometric errors induced by the instrument. To do so, the covariance matrices of the misfit between the SEPS-retrieved and the forward model brightness temperatures, sorted and grouped by the number of observations (Nobs), have been computed. In addition, as a test, eigenvalue decomposition has been applied and the number of eigenvectors required to describe the 99% of the variance has been defined as the effective number of measurements (Neff). Its trend as a function of the number of observations has been analyzed, and the results suggest the presence of two regimes: the first one is noise dominated, where Neff is almost equal to Neff; and the second one is where Neff increases with Nobs according to the constant slope of 0.8.
Introducing Neff in the cost function resulted in applying a weight to the average residual term of the observational part of the SMOS ocean salinity level 2 cost function equal to the factor . The consequent impact has been assessed by either comparing the retrieval performance with that obtained using the estimated covariance matrices directly or with both of the cost functions present in the SMOS literature.
Conclusions can be summarized by the following three points:
Based on the two regimes of Neff, a threshold can be established to define objectively the useful swath of SMOS as 700 km centered on the satellite ground track, where the relation Neff/Nobs is constant.
The three approximated cost functions [Eqs. (4), (5), and (8)] give very similar performances. The analysis of the cost functions suggests that both the current configurations [Eqs. (4) and (5)] and the proposed weight [Neff, see Eq. (8)], although ensuring better performance, may be introducing nonlinearities in the retrieval procedure if compared to the results obtained using the misfit covariance matrices directly. According to previous studies (Sabia et al. 2010), nonlinearities may be due to a nonoptimum balancing of the cost function that should be modified.
Furthermore, the inclusion of the brightness temperature misfit covariance matrices strongly modifies the error statistics, revealing the need for a more accurate modeling of the correlation in the brightness measurement misfits. This must be introduced in the cost function and its impact on the relative weights must be applied to the assessed auxiliary parameters.
The authors would like to acknowledge the anonymous reviewers whose valuable suggestions and remarks contributed to the strengthening of this paper, as well as Dr. Baptiste Mourre for providing NEMO-OPA geophysical auxiliary data. This study has been founded by the Spanish National Program on Space through Project ESP2005-06823-C05 and the Spanish Ministry of Science and Innovation through the Formación de Personal Investigador (FPI) Fellowship ESP2005-06823-C05-02.
Current affiliation: Serco SpA, Frascati, Italy.
Current affiliation: European Space Agency-ESRIN, Frascati, Italy.