The future Surface Water Ocean Topography (SWOT) mission aims to observe water bodies and short-scale ocean surface topography with unprecedented spatial resolution and accuracy. However, the topography estimates will be contaminated by errors of various signals (geophysical and instrumental) featuring, in large part, strong dependencies on the radar range direction (cross track). This study shows that a cross-spectral analysis performed along track for all cross-track combinations can detect most of these errors and can provide estimates of their power spectral densities. From a series of outputs of the SWOT science team simulator, a cross-spectral method was developed to simulate the estimation of the error budget compared to the actual error budget in the simulator. The study determined that the error spectra of the dominant terms can be estimated at very high accuracy. Beyond the obvious applications for the future SWOT data calibration and validation, the spectral estimates of the error budget will have applications for state estimate problems using SWOT data.
The Surface Water Ocean Topography (SWOT) mission (Fu et al. 2012) is expected to provide two-dimensional measurements of ocean and inland water surface topography at unprecedented spatial resolution of a few kilometers (Esteban-Fernandez 2013), over a 120-km-wide swath. The main instrument of the payload’s mission is a Ka-band Radar Interferometer (KaRIN), which has a 10-m baseline (Rodríguez 2016). The satellite will follow a 21-day exact repeat cycle, allowing full coverage of Earth. The first 60 days of the mission will be dedicated to calibration and validation (CalVal) of this new altimetry data. For this purpose, the satellite will first follow a 1-day exact repeat orbit, providing samples of daily repeated observations (Rodríguez 2016). Despite the high frequency of revisits during this phase, the ocean surface topography is expected to vary significantly day to day (e.g., Dibarboure and Morrow 2016; Qiu et al. 2017). This will be challenging for the data validation. since the observed differences will not be necessarily attributed to measurement errors (Fu and Ubelmann, 2014). Moreover, as described in Esteban-Fernandez (2013), the measurement error is the sum of various terms (e.g., KaRIN noise, roll errors, phase drift errors) with potentially comparable amplitudes. Measuring the separate error terms during the CalVal phase is therefore a real challenge.
In this study we propose an approach based on the internal consistency of the 2D SWOT sea surface height (SSH) to measure and to sort out all range-dependent errors. A cross-spectral analysis is applied, followed by a decomposition into isotropic- and cross-track-dependent signals. If standard along-track spectral analysis provides only a spectrum of the sum of the components (SSH plus errors), we will show that cross spectra bring additional information that allows for distinct separations. Indeed, cross spectra are sensitive to the cross-track structure and to the two-dimensional properties of the signal. The use of cross-spectral analysis for synthetic aperture radar at different ranges has already been successfully applied (Engen et al. 1994; Bao and Alpers 1998; Dowd et al. 2001) to retrieve ocean wave spectra. In altimetry, Yale et al. (1995) and Smith (2015) applied cross-spectral analysis to assess the coherence between repeat tracks.
To implement and test the cross-spectral approach proposed in this study, we used a synthetic SWOT dataset from the SWOT Ocean Science Team simulator (Gaultier et al. 2016), presented in section 2. The cross-spectral signatures of the different SWOT signals and error components from the simulator will be highlighted in section3. Then, section 4 will describe the decomposition method; last, the results will be presented and discussed in sections 5 and 6.
2. Experimental setup
The SWOT Ocean Science Team simulator (Gaultier et al. 2016) is used to generate synthetic observations of SSH as sampled by the future satellite, featuring random and systematic instrumental errors. The ocean input scene and error budget scenario are described below. For technical details about the software, the reader can refer to the simulator website (https://github.com/SWOTsimulator). In this present study, we consider 60 days of the mission’s fast sampling phase (1-day exact repeat period), as illustrated in Fig. 1.
a. Sea surface height inputs
We use an SSH scene from a numerical run of the MITgcm (Marshall et al. 1997), a primitive equation model. This numerical run, called llc4320, resolves mesoscale eddies, internal tides, and other hydrostatic processes at scales as small as about 5 km. The llc4320 is the highest resolution (1/48°) of a hierarchy of simulations. Successful validations of the run have been recently published (e.g., Rocha et al. 2016), suggesting that the variability, in particular for the surface topography, is realistic for both the low- and high-frequency dynamics. This is the first global run covering such a large spectrum of temporal and spatial scales; therefore, it is very interesting for providing realistic scenes of SSH for SWOT, to the extent of modeling capabilities. We use the 84 days of llc4320 outputs spanning September 2011 through December 2011.
b. Error budget scenario
The systematic and random errors added to the sampled SSH follow the spectral allocations given in Esteban-Fernandez (2013). These allocations, expressed in along-track power spectral density, are shown in Fig. 2, which was converted in equivalent SSH and averaged in the KaRIN swath (10–60-km cross track). For each error component, the simulator generates a random parameter in the physical space, composed of a series of harmonics whose amplitudes follow the prescribed error spectrum. Homogeneous and independent phase distributions are assumed for all harmonics. The random parameters are then projected into the topography error in the swath, with specific geometries such as a linear cross-track slope for the roll and phase errors or a parabolic shape for the baseline dilation error.
c. 1-day differences
The 1-day differences of SWOT synthetic observations during the fast sampling phase are considered in this study. Indeed, working on the differences will cancel out a significant part of the medium- to large-scale SSH, making the cross-spectral analyses more efficient for error budget estimation (Dibarboure and Morrow 2016). Note that the fields of differences contain twice the error budget as the error fields are supposed to vary at short time scales (less than an hour). The estimations performed in this study will therefore apply to the 1-day error differences that are twice the spectra shown in Fig. 2.
A sample of the difference field is shown in Fig. 3 (upper panel) along a 700-km-long segment of swath in the North Atlantic Ocean. One can clearly identify some patterns owing to the 1-day difference field of SSH (δSSH) shown in the last panel, as well as a clear contribution of the error sample differences shown in the other panels. This illustration confirms the ocean validation challenge mentioned in the introduction: the separation of the observed field into the different components does not appear obvious in physical space, even during the fast sampling phase. In spectral space, the direct power spectral density computation would reveal the sum of the component spectra (black line in Fig. 2), but the separation would not be trivial. However, we will see that a cross-spectral analysis can exploit the 2D specificity of the error terms, allowing a separation.
3. The along-track cross-spectral density of the SWOT observations
a. The SWOT “XSD-cube”
Given the geometry of the SWOT grid [two ribbon-shaped images (50 km wide in the range direction) with a 20-km gap near the subsatellite point] and the geometrical dependencies of errors in the range direction (e.g., linear for roll), it is relevant to perform the 2D signal analysis in Fourier space in the along-track direction while being in physical space in the cross-track direction. This can be achieved by considering the cross-spectral density (XSD) between along-track signals at different range positions.
If and are two monodimensional discrete series of along-track SSH ( being the along-track coordinate) at different range positions and , respectively, then the XSD between and is defined as
where FFT is the Fourier transform and the asterisk (*) denotes the conjugate. The term is also monodimensional, real, and expressed in the along-track wavenumber kx dimension. Note that for the particular case , the XSD is the power spectral density (PSD) of .
From the SWOT swath, the “XSD cube” can be defined by the concatenation of XSD between all cross-track values, as represented in Fig. 4. The result is a function of ri, rj, and kx where (ri, rj) is the matrix of range combinations. If is the range dimension, then there are combinations [the combination matrix (ri, rj) being symmetric]. Practically, considering a 1-km square grid (the expected posting resolution of ocean KaRIN products), the range dimension is 102, which makes individual XSD computations for a given SWOT sample. For any given along-track wavenumber kx, a single 2D slice (ri, rj) of the 3D XSD cube can be interpreted as a “covariance matrix” (symmetric positive) for a particular wavelength. The diagonal terms (ri = rj) are the PSDs (or variances at kx). In the following, we will show the specificity of these slices for the different components of the SWOT signal.
b. The distinct XSD signatures of the SWOT signals
The different components of the SWOT signal have distinct signatures in cross-spectral density. In this section we highlight the existence of a specific analytical model for each signal, allowing for the distinction methodology in section 4.
The table in Fig. 5 gives an overview of all the analytical models. The first column is the topography signature of each signal in physical space. These expressions are derived into the XSD as a function of wavenumber and a couple of cross-track values (ri = rj), as shown in the second column. In each case we show (in the following) that the XSD can be decoupled as a product of a function of wavenumber only by a function (noted as g) of (ri = rj). This is what we call the analytical model, later, which is plotted in the third column.
1) Roll, phase, baseline dilation, and timing errors
The roll signature can be defined in the swath by a single along-track parameter multiplied by a geometrical function of range , as shown by Eq. (1a) in Fig. 5. This decomposition, decoupled in and , makes it possible to derive the XSD analytical expression between two series at ranges ri and rj, by applying Eq. (1). Assuming that has an along-track PSD of and taking out of the FFT operators in Eq. (1) (independent of ), the XSD analytical expression [Eq. (1b) in Fig. 5] is directly obtained. For a given , this XSD is the product of the roll parameter PSD at by a signature function plotted in the third column. The function is also the covariance, in the range space, of the roll signal normalized by the along-track roll parameter. For and of opposite sign, note the negative values expressing the negative correlation of the roll signal between opposite swaths (as a result of the antisymmetric nature of roll topographic signature in physical space).
The XSD of phase, baseline dilation, and timing signatures can be similarly expressed as shown in the second, third and fourth lines of Fig. 6, respectively. The parameters are noted as , and , respectively.
The signature functions are all specific and independent (i.e., none of them is a linear combination of others). The phase function is the same as roll when and have identical signals (i.e., same swath), but zero elsewhere because the phase errors are uncorrelated between the left and right swaths (no correlations for the opposite range). This difference of signature functions between the roll and phase will be the basis for their separation. The baseline dilation function is quadratic, and the timing function is unitary because of its constant value with range.
2) KaRIN noise
The XSD of the KaRIN noise can also be expressed under a similar decomposition. In physical space (in the swath), the KaRIN noise is defined by independent Gaussian realizations of zero mean and variance . The variance is directly proportional to the 2D flat spectrum of level and the grid size area [as shown by Eq. (5a) in Fig. 5]. Note that the are modulated by the range and the sea state essentially, as described in Esteban-Fernandez (2013). For , the independence implies that the XSD is zero [see Eq. (5b)]. For , the XSD (which is the PSD) is nonzero, of value . To follow a similar decomposition of XSD in a product of a PSD (to be estimated later) by a signature function , knowing that is range dependent, we split the XSD expression into the sum for each range, as shown in the table. For a given range , the XSD of the range of the KaRIN noise between and is zero everywhere except for , as shown in the signature function in the third column. The XSD of the full KaRIN noise (sum for all ) results in a diagonal signature function.
3) Isotropic signals
The rest of the SWOT signal (the 1-day topography difference and potentially other geophysical error differences) can be reasonably well approximated as isotropic with respect to range and azimutal directions. Indeed, we will stack SWOT images from ascending and descending orbit segments and a wide range of headings (0°–77°). So, on average geophysical signals should not have a strong azimuth/range signature. Under this isotropy approximation, it is possible to express the XSD between two along-track transects as a function of the PSD (noted as ) and write its expression under a similar form as the abovementioned signals, that is, the product of the PSD with a signature function. As shown in the second column of Fig. 5 and demonstrated in appendix A, this function, noted as , is defined as
c. Qualitative detection of specific signature functions in the observed XSD cube
Figure 6 shows some slices of the XSD cube at three given wavelengths. For the long wavelength, at 5000 km, the XSD has a clear dependence on , stronger for and of identical sign (see the red patterns in the diagonal blocks, where intensity increases with range). This originates in the dominance of phase drift and roll errors at long wavelengths. For the intermediate wavelength, at 200 km, the XSD slice has less range dependencies and seems to be closer to a function of the absolute range distance . This is the characteristic of the range-independent signal, dominated by δSSH at this wavelength. For the short wavelength, at 10 km, the diagonal terms are clearly dominant, showing the presence of an uncorrelated signal whose intensity is minimum between 20 and 45 km, and increases at near and far range. This signal results from the KaRIN instrument speckle noise.
If the dominating signatures show up visually in the XSD slices, then we will demonstrate in the following section that a least squares inversion can provide a quantitative estimation of all individual signatures from the XSD slice.
4. Methodology for estimating the PSDs from the SWOT XSD cube
a. Step 1: Computation of the SWOT XSD cube
The first step of the method is the computation of the XSD cube from the input dataset of the SWOT 1-day topography difference. The data are split into a 3000-km-long segment over the ocean excluding land. This length is a trade-off between the maximum wavelength until we want to estimate the error budget (about half the segment, here 1500 km) and the number of possible land-free segments. From each individual data segment, an XSD cube is computed as defined in section 2. A detrending followed by a 10% tapering (sufficient from sensitivity tests) is applied to the series, mitigating the effect of nonperiodicity as generally performed in 1D spectral analysis of altimetry data. The cubes are then averaged for all segments.
b. Step 2: XSD least squares decomposition and PSD estimations
The observed SWOT signal is the sum of the component described in section 1. If we assume that the error terms are uncorrelated (see the discussion in section 5), then the observed XSD is the sum of each component XSD that writes under the signature function form
where is the residual not projecting onto the signature functions.
The goal is to estimate , and the from the observed cube. The problem is solved separately for each slice (given ) as illustrated in Fig. 7, where the different terms of Eq. (8) are represented. If we note , then the slice of the cube [ is flattened in a one-dimensional vector for all couples], η is the vector of concatenated parameters to estimate (at ) of dimension n = 108 if the range dimension is 104, and is the matrix of concatenated functions, of dimension , Eq. (8) writes
We verified that the determinant of is nonzero, meaning there are no linear dependencies between the functions. In these conditions, it is possible to compute a least squares inversion to estimate the η minimizing :
where ηa is the η estimate. This inversion, computed separately for all , provides an estimation of the 1D parameter along-track spectra describing the SWOT error budget.
As mentioned above, the signature function , for a given , has a dependence on the spectral slope at all wavenumbers above (shorter absolute wavelengths projecting onto in the direction). To handle this dependency on , an ensemble of inversions is performed at different between −5 and −1 (covering a wide ensemble of possible spectral slopes for the isotropic signals). The different residual are compared and the slope solution of the minimum RMS for is chosen. We verified that the slope solution obtained for each is in good agreement with the slope of the reconstructed spectrum of the isotropic signal.
The uncertainty can be provided for each estimation. Indeed, the covariance matrix of the estimation error can be related to the least squares mismatch as follows:
where is the sum of the squared mismatches on the XSD slice divided by the degree of freedom,
where , the vector size minus the number of parameters to estimate, is the degree of freedom. The error bars on the parameters are defined from the diagonal of the matrix. Note that these error bars stand for the estimation of the spectral parameter owing to the ensemble of data samples. However, the spectral parameter has its own uncertainty with respect to the actual spectral parameter (here, the one defined in the simulator), given by the standard χ law based on the number of independent data samples. It relies on the assumption that the SWOT time series are independent (not true if the samples are too close in time with respect to SSH variability); therefore, they potentially lead to underestimated uncertainty. The combination of these two uncertainties will be considered in the error bounds presented in section 4. We will verify that these error bounds appear realistic with respect to the true spectra.
The results of the cross-spectral method applied to the 60 days’ worth of 1-day fast sampling phase SWOT synthetic data are presented in this section.
a. Effective reconstruction of spectra
The roll and phase are the dominant error terms at medium and long wavelengths, as shown in the error budget in Fig. 2. The estimation of their PSDs resulting from the least squares inversion is very accurate at all wavenumbers, as shown in the left panel of Fig. 8 with respect to the simulator inputs. Although both the roll and phase have a linear cross-track signature on topography, they have been spectrally separated thanks to their distinct XSD signature on the off-diagonal block of the cube slice (phase error correlation being zero for the opposite range). At 50-km wavelength, where the phase error is 10 times less energetic than the roll error, the estimation of phase is still successful, despite the relative increase in uncertainty indicated by the error bars.
The estimation of the baseline dilation, which is about 10 times less energetic than the roll and phase, is still accurate above 100 km and below 10 km. Between 100 and 10 km, the spectral values are more uncertain, as indicated by the error bars. Note this is the spectral band where the baseline dilation happens to be the smallest with respect to the sum of the SWOT signal (between two and three decades, a minor contributor of the error budget), as shown in Fig. 2.
The estimation of the topography difference is also successful and very accurate down to about 15 km. In this experiment the barotropic tides and the other short-frequency signals of the long wavelength (such as inverse barometer) have not been filtered, explaining why the spectrum of the 1-day topography difference is still red at the long wavelength. This was not a problem for its estimation. Below 10–15 km, the estimation becomes impossible because the amplitude of noise largely dominates the XSD signal.
The timing errors, as weak as they are in the simulator, cannot be accurately estimated as shown in the first panel of Fig. 9. Its uniform XSD signature is certainly less distinctive with respect to other signatures. In particular, at the long wavelength, the signature is getting similar to δSSH. However, the upper limit of the confidence interval is in this case useful information (upper gray area in Fig. 9a): if the spectral values cannot be detected, then we know that they are below a certain limit that is very small with respect to the total error budget. If the timing error happens to be higher, up to the level of threat with respect to the global error budget (a multiplication of 10 or more, in amplitude), then it would be accurately detected as indicated by the right panel if Fig. 9. Sensitivity experiments (not shown) revealed that this detection was possible down to a factor of 10 below δSSH.
b. KaRIN noise budget
The KaRIN noise budget estimation is very accurate. As shown in Fig. 10, the integration of the spectral estimation for all wavelengths gives the expected noise variance corresponding to the gridcell area. The near-range and far-range noise increases are well captured. The KaRIN noise being dominant by more than a factor of 10 above all other components at the shortest wavelength (4 km in the simulator) and uncorrelated by assumption, such an accurate estimation, was expected. Note that if the local KaRIN noise PSD were below the PSD of δSSH, it could not be detected. There are no range-dependency-allowing estimations at lower energy as for the roll, phase, and baseline dilation terms discussed above.
c. Can we capture rapid changes in SWOT errors?
The results presented above were obtained for a series of 60 days, which is the expected duration of the mission’s fast sampling phase, and are averaged over the whole ocean. However, although unlikely, some error components might be unstable or might not meet the requirements. So, it is important to be able to detect and quantify it as fast as possible. Figure 11 shows the roll, phase, and baseline estimation after 2 days, to be compared with Fig. 8. Indeed, the 2-day estimation is degraded, but the dominant terms (here, roll and phase) are still fairly captured. This result indicates the capability to quickly detect variations in dominant error terms. Therefore, in the nominal case (if the error budget is actually close to present expectations) the budget would be confirmed after a few days for roll and phase, and slightly more for baseline, which has very low expected energy.
Similarly, if the weak components (e.g., baseline and timing) were increased by a factor of 10, we found that they would be accurately measured after a few days. The cross-spectral method has the advantage of quickly providing, after a few days, some upper-interval estimates of each component even significantly below the spectrum of the total SWOT signal. Then, if these components are actually very weak as per the specifications, then more time is needed to reduce the interval of confidence around the actual spectral values.
a. Validity of assumptions
The spectral parameter estimation in this study relied on a few assumptions, which are discussed below.
First, we assumed no correlation between error components. In the Fourier decomposition, the spectral phases of the signals were independent in Eq. (8). Note that this assumption is not violated in the presence of a correlation between the amplitudes of the errors, which may likely happen during events such as eclipse transitions (thermal effects) or maneuvers that would increase all error terms at the same time. However, the phases of the error terms must remain uncorrelated. Regarding roll, phase drift, baseline dilation, and timing errors, the assumption is certainly reasonable, since the signals are generated by different physics (platform, mechanical vibration, electronics). Regarding residual errors as a result of sea-state bias, possible correlations are anticipated with SSH (as existing in conventional altimetry; Zlotnicki 1994). This point will deserve particular attention when SWOT sea-state bias simulations are available. Note that residual correlations between two terms would have specific signatures on the SWOT cross-spectral cube, which could be studied and detected as well in the future.
Another important assumption is the isotropy for geophysical signals at 1-day differences. From our simulations, which contain only modeled high-resolution ocean surface topography among the geophysical signals, the isotropy stood very well, even in the presence of internal tides, and had no impact on the estimated spectra. A major geophysical component not simulated in this study is the wet tropospheric signal. However, at 1-day differences, only the small spatial scales of the wet troposperic signal should be highlighted. These latter are reasonably well isotropic as discussed in Ubelmann et al. (2014). Note that if the cross spectra are computed from the corrected data using the two-beam radiometer, then some range dependencies may occur. The XSD signature would therefore be specific, but it could be anticipated.
b. Method extension
Beyond the results of our simulations considering a supposedly known ensemble of error geometries, the cross-spectral analysis can be constituted as a tool to study errors not anticipated in simulations. Indeed, the SWOT XSD cube itself highlights a wealth of information on the SWOT topography signal, in particular on its range dependencies. It could prove insightful to analyze the cube slices for all wavenumbers to determine the potential signatures on which the SWOT signal may project. This cube analysis should be very instructive and help with understanding the SWOT error budget.
Finally, other variables may extend the cube analysis. The nadir altimeter data cross spectrum with KaRIN data would surely bring useful information for CalVal purposes. Beyond the height variable, we might also consider the radar brightness or spectral coherence from the SWOT interferograms in order to detect range dependencies.
This study has shown the potential of a cross-spectral analysis on the future SWOT data to measure the spectra of different error and signal components. In particular, structured error terms such as roll, phase, and baseline dilation errors can be measured at very high accuracy even if the energy is not dominant with respect to the total signal and if they are not detectable from SWOT images in physical space. Since the SWOT performance requirements are defined in spectral space, the method presented in this study should be useful for validating the error budget during the CalVal phase. A few days’ worth of data will be sufficient to characterize the major potential error sources, and the full 60 days of the CalVal phase should allow for accurate error quantifications and separations.
Once the mission is on its nominal orbit, the error spectral estimates will be a useful quantity for signal denoising and for parameterizing optimal filters for higher-level products (Ruggiero et al. 2016). Indeed, the error spectral estimates can be easily converted into error covariances used for any state estimate problem. The mapping of SWOT data in a continuous gridded product will benefit from such accurate error covariance matrices to perform optimal inversions. Also, more sophisticated systems, such those combining SWOT data and ocean circulation models through data assimilation (M. Benkiran et al. 2018, unpublished manuscript), will require these matrices.
This study was funded by the Centre National d’Etudes Spatiales (CNES) as part as the Algorithm Design Team and CalVal support (Contract CLS.DOS.BL-16.072).
Relation between the XSD and 1D PSD for 2D Isotropic Signals
Let us consider a 2D isotropic signal , where and are two orthogonal directions, respectively, of 2D power spectral density . The XSD between two transects, noted as , is the 1D inverse Fourier transform of the 2D PSD along the axis,
where is the distance separating the transects. Isotropy assumption Is shown below:
The integration of is performed for absolute wavenumbers equal to or higher than the given . In the case of a red spectrum (most of the geophysical signals), say, with a slope steeper than −2, only absolute wavenumbers slightly higher than will have significant weight in the integration of Eq. (A2). Over such a limited range of wavenumbers above , geophysical signals can be reasonably well approximated by the tangent linear of the PSD in log scale (assuming the spectral slope variations are weak). We can therefore write as
Using appendix B where is the 1D isotropic power spectral density,
This Eq. (A5) links the cross-spectral density in the direction between two transects separated by in the direction, noted as , with the 1D along-track spectral density. We have not found any obvious literal expression for this integral. For our purposes, it will be computed numerically.
Relation between 1D PSD and 2D PSD for 2D Isotropic Signals
Let us consider a 2D isotropic signal , where and are two orthogonal directions, respectively, of 2D power spectral density . Its 1D isotropic power spectral density along axis writes
which is defined only for ; this is why we have a factor of 2.
Since is isotropic, can be defined by a function of the absolute wavenumber
Case of spectra following a power law
Let us assume follows a power law, writing under the form
Then Eq. (B3) writes
It can be shown that
Using Euler transformation, we have
Using the Gaussian theorem, we have
where is the gamma function, the extension of the factorial function for any real number.
And finally using Eq. (4), we have the following expression for as a function of :
The spectral slope of the 1D spectrum is therefore higher than one unit compared to the 2D spectral slope. If a 2D isotropic signal obeys a −4 power law in 2D spectral space, then the 1D spectrum computed along 1D transects will follow a −3 power law. Note that a function of the slope modulates the magnitude of the 1D spectrum.