Global altimetric mean sea surface derived from the geodetic phase of the ESA ERS-1 mission utilising a spectral least squares collocation technique

Abstract:

This paper describes a method that enables a high resolution global Mean Sea Surface (MSS) to be derived by a least squares collocation technique within the frequency domain via the use of the fast Fourier Transform (FFT). The MSS is produced from the entire two cycles of the ERS-1 Geodetic Phase consisting of 336 days of off-line Precise Ocean Product (OPR) data. This data is corrected for tidal and environmental/electronic effects. Long-wavelength radial orbit error is minimised via the fitting of cubic spines across ERS-1 and TOPEX/Poseidon dual crossover residual time series. After further editing, the remaining data is gridded at intervals of (a limit imposed by the cross-track spacing at the equator allowing a resolution of 16 kilometers) whilst non-oceanic regions are filled (a requirement of the FFT) with a reference geoid derived from the OSU91A gravity field. The spectral least squares solution draws on the method of objective mapping within the spatial domain by transforming the whole process into the frequency domain thus relinquishing the spatial requirement for matrix inversions. Finally the resultant MSS is compared against both an independently derived TOPEX/Poseidon reference MSS revealing a RMS fit of 17.9 cm and 13.6 cm against TOPEX/Poseidon altimetry itself. TOPEX/Poseidon altimetry had a 17.1 cm fit with the reference MSS.

Keywords: Marine Geoid, FFT, Optimal Interpolation

1. Introduction

2. Mean sea surface solution

• 2.1 Altimetry Preprocessing
• 2.2 Generation of Gridded MSS
• 2.3 Spectral Optimal Interpolation
  • 2.3.1 Generation of Residual MSS and Altimetric Noise
  • 2.3.2 Autocovariance of Reference geoid
• 2.4 Solution

3. Conclusions

4. References

1. Introduction

Accurate estimation of the marine geoid is of particular interest to oceanographers and geodesists but its determination from satellite altimetry poses a problem in that each measurement of the ocean surface height contains the marine geoid, quasi stationary sea surface topography, tidal phenomena and a range of errors relating to the measurement itself. Most of the measurement errors can be corrected to a high degree which enables us to obtain accurate measurements of the surface which is a composite of the marine geoid and quasi stationary sea surface topography; this surface is known as the mean sea surface (MSS). Varying methods for its estimation from satellite altimetry have been under investigation since the facility became available in the 1970's. Following regional studies conducted by  [Wunsch and Zlotnicki, 1984],  [Mazzega and Houry, 1989] and  [Blanc et al., 1990] based on the theory of objective analysis, results have shown that the technique yields accurate determinations of the MSS. Unfortunately this method is limited to small geographical regions ( ) as it requires the inversion of a covariance matrix describing the total estimated error within the MSS; which in turn requires the generation of a covariance matrix for the altimeter data errors. These limitations may be reduced by considering the problem in the spectral domain with matrix inversion in the spatial domain replaced by division at each frequency. This paper describes such a technique with applications to deriving a global high-resolution MSS from ERS-1 satellite altimetry collected during the Geodetic Mission (GM). A notable advantage of the technique is the speed of computation and, in principle the relatively ease with which different data sets can be combined. The method has been validated against a reference MSS [Basic and Rapp, 1992] and TOPEX/Poseidon sea surface heights (SSH), however a rigorous analysis of the formal errors has not been undertaken as yet.

2. Mean sea surface solution

Altimetry was taken from the ERS-1 geodetic mission fulfilling the criteria of a high resolution data set. In the first instance altimetry data were edited for erroneous points and points over land/ice, leaving data over ocean and shallow water for correction and gridding. The MSS to be spectrally optimally interpolated is transformed to the frequency domain by the 2-dimensional FFT. However since this facility requires a complete grid and void regions exist over land/ice it is necessary to augment the MSS with reference heights, in this case calculated from the OSU91A geopotential model [Rapp et al., 1991]. The spectral least squares procedure requires Power Spectral Density (PSD) knowledge of the errors resident within the MSS and noise within the altimetry. The first of these is supplied by the auto-covariance (ACV) function of a reference geoid calculated from geoid error degree variances of the OSU91A between degrees 2 and 360 and extrapolated to degree 2400 by a power decay fit. The noise PSD is derived from the altimetry itself. A combination of these power spectral densities (section 2.3) supplies us with a spectral least squares estimator which when multiplied with the complex spectral MSS gives us a spectral optimally interpolated MSS for transformation back to the spatial domain. Each of these processes is now described in detail.

2.1 Altimetry Preprocessing

Off-line Precise Ocean Product (OPR) data for the entire ERS-1 GM between 10th April 1994 and 23rd March 1995 (MJD 49454 - 49790) were used. This period consisting of 2 cycles with a repeat period of 168 days seperated by a manoeuvre at the end of the first cycle to the extent that the second cycle ground track was offset from the first in order to create an interleaved data set. Each cycle consisted of 2411 passes giving an equatorial cross-track spacing of approximately 8 kilometers ( ) and allowing a maximum MSS resolution capability of 16 kilometers ( ) due to the Nyquist sampling theorem. The altimetric geophysical data records were edited to remove land/ice flagged data leaving ocean/shallow water information to supply as much sea surface information as possible despite tidal definition over shallow water regions. Large fluctuations in sea surface height may occur at land/ice boundaries due to errors in land/ice masks resulting in reflections from shorelines or ice boundaries. Such points were removed by comparison with a reference surface. The remaining data are corrected for ocean tide (CSR3.0), solid Earth tide, polar tide, wet/dry tropospheric delay, ionospheric delay, inverse barometric effect and bias drift. Sea state bias is calculated as 5.95% of the significant wave height [Carnochan, 1996].

Finally, long-wavelength radial orbit error is reduced by the use of a dual TOPEX/Poseidon and ERS-1 crossover generating procedure [Carnochan et al., 1994]. Here, dual satellite crossover residuals for the MJD period 49449 - 49795 were generated with a de-correlation period of 5 days resulting in 333058 crossovers for this period. After thresholding of the residuals their RMS were approximately 11-12 cm. A correction for radial orbit error was calculated by fitting cubic splines through the time series of dual satellite crossover residuals. This procedure differed from [Le Traon et al., 1995] in that single ERS-1 crossover residuals were not incorporated. It is however accepted that their inclusion will be beneficial for extreme latitudinal areas not covered by dual crossovers and regions such as near Indonesia and the Mediterranean. However in this study ERS-1 altimetry covering the high latitudinal extremities are later used as regions for windowing by the FFT such that the above mentioned geographical regions are of little use in this study. Hence, it was considered justifiable to only incorporate TOPEX/Poseidon and ERS-1 dual crossovers and use altimetry outside for FFT windowing (section 2.3). The criteria for selecting spline knot locations is similar to [Le Traon et al., 1995]; a knot is placed at the first and last dual crossovers on a pass consisting of more than 10 crossovers and a third knot is placed at the center crossover if more than 5 crossovers exist within the pass. With this procedure cubic spline RMS fits to the residuals was 6 - 7 cm.

2.2 Generation of Gridded MSS

Altimetric data was gridded using bi-linear interpolation at a longitudinal and latitudinal grid spacing of set by the cross-track spacing at the equator. This means that the grid is somewhat oversampled at all other regions covered by ERS-1. However, for this particular technique, we require a regular latitude/longitude grid and are thus left with this limit.

Since the data set contains altimetry over all water surfaces it is necessary to remove land/water - ice/water boundary points that differ from a reference surface (OSU91A plus a SST model) by more than a meter when plotted such points existed around boundaries. The procedure for calculating an optimal surface (section 2.3) via the FFT requiring a complete grid, void regions over continents, ice masses and islands therefore need filling. A maximum likelihood method was tried in order to fill these regions but failed for all but islands and peninsula's. Void regions were therefore filled with a geoid derived from OSU91A. All references to the OSU91A derived geoid incorporated a corrective offset [Rapp et al., 1994].

2.3 Spectral Optimal Interpolation

In the spatial domain the optimal estimate, , at a point with geodetic coordinates is given by [Mazzega and Houry, 1989]

where, represents our gridded final optimal estimate, and is a initial long wavelength estimate (20 by 20 OSU91A derived geoid) of the mean sea surface to be removed from the altimetric data set . This long-wavelength surface is added on again later with the premise [Wunsch and Zlotnicki, 1984] that this removal reduces the signal power leaving errors within the OSU91A and sea surface topography thus allowing an easier comparison of remaining signal and noise. There are two covariance matrices namely; , the covariance of errors describing our initial estimate and the covariance describing the errors within the altimetric signal. In practice there is difficulty implementing equation (1) if large regions are to be studied due to the inversion resulting in heavy time and computational consumption. It is therfore desirable in the circumstance of global MSS analyses that a more suitable method is obtained.

Such a solution may be found by the transformation of equation (1) into the frequency domain by the use of the FFT. Here matrix inversions translate to division at each matrix element, hence it should be possible to analyse regions of varying size and resolution. Following [Bracewell, 1978] and [Schwarz et al., 1990], given two signals a and b, their power spectral density may be calculated either directly from their initial form or alternatively from their associated covariance structure . The former method is represented as

where the Fourier operator is given by and denotes the complex conjugate. The second property is attained by

where, and are the signal means. If either or are zero then . Using the facilities given by equations (2) and (3), equation (1) may be converted to the frequency domain giving,

or,

where, denotes the inverse FFT operator for the transformation from frequency to spatial domains. Also , and are transformations of , and respectively.

The associated covariance of the error estimation may be calculated from

and its spectral equivalent is

It is therefore possible, in theory, to calculate the optimal estimate by constructing the spectral equivalent of equation (1) through four FFT transformations giving , , and .

2.3.1 Generation of Residual MSS and Altimetric Noise

A requirement of equation (4) is the generation of a grid representing the errors resident within the altimetric measurements. Deriving such a grid poses a problem in the sense that it is necessary the have knowledge of noise contributing sources on both a temporal and global scale; such a huge task is undesirable. In previous studies, (for example  [Wunsch and Zlotnicki, 1984], [Mazzega and Houry, 1989] and [Blanc et al., 1990]) over small geographical regions, the altimetric noise budget is calculated as the sum of a number of modelled covariance functions describing each of the predominant noise sources. More recently [Tsaoussi and Koblinksy, 1994] developed a model for the calculation of the error covariance for sea surface topography by incorporating scaled biases for each of the altimetric corrections.

The method used in this study requires the calulation (equation 7) of a residual surface obtained by removing a reference geoid and sea surface topography from our initial MSS thus leaving a surface that contains errors within the OSU91A, the SST and altimetric measurement errors. This surface is then scaled to reduce the long-wavelength component (Fig.1(top)) and used as a model describing altimetric errors. Here the use of altimery in this model is clearly an advantage however the existence of long wavelength component detracts from it. The spatial noise surface is therfore given by

where, A is a scaling constant. The PSD ( ) is then calculated using equation 2, a 1-dimensional version is displayed in Figure 1(bottom).

2.3.2 Autocovariance of Reference geoid

The autocovariance function of the error within the initial estimate of the MSS is given by

where, is the mean radius of the Earth, n the degree, the maximum degree, are Legendre polynomials and is the spherical distance between two points. are the geoid error degree variances at degree n calculated between degrees 2 and 360 from the OSU91A potential coefficient standard deviations and by

where m is the order. For n=361 and above the geoid error degree variances are given by the power decay law (see figure 2 for the extended geoid error degree variances).

Given this arrangement is calculated over the globe at the same grid specifications as previously described. Since the gridded ACV is perfectly symmetrical it need not be windowed. is calculated using equation (2) with zero mean ( ).

2.4 Mean Sea Surface Solution

All grids are periodic in both latitude and longitude with the exception of the autocovariance which is purely symmetric. For consistency all grids are Hanning windowed over latitudinal extremities ( ) to avoid spectral leakage. Following the production of spectral grids , , and the computaion of can take place (eq. 4) the result of which is complex though the imaginary component is of the order of due to edge effects.

In order to validate the solution it is necessary to make a comparison with other MSS models. At the time of writing limited validations were made against TOPEX/Poseidon sea surface heights and a MSS model derived from GEOS-3, SEASAT and GEOSAT data [Basic and Rapp, 1992] available on AVISO distributed CDROM s. Data for both TOPEX/Poseidon and the reference MSS were taken for cycle 66 of the TOPEX/Poseidon mission. A comparison of our MSS revealed a global fit of 17.9 cm with the reference MSS and a fit of 13.6 cm against TOPEX/Poseidon data; The TOPEX/Poseidon sea surface height fit with the reference MSS was calculated to be 17.1 cm. Differences are plotted in Figures 3 and 4, in Fig. 3 there are some large fluctuations which may arise as a result of model differences or the existence of variability infomation due to difference in time between observations. Also as comparison is made over shallow waters differences may occur purely from localised tidal characteristics. Comparison with TOPEX/Poseidon in Figure 4 is more promising. With the collection of stacked repeat passes it would be hoped that these fluctuations reduce to a more desirable level. In order to fully calibrate the model it is necessary to calculate a formal error covariance for the model which may be obtained from the use of equation 6. Finally two examples of the optimally interpolated MSS or presented for two regions with wavelengths less than 2000 km. The first of these is for the North Atlantic (Figure 5) with the Mid-Atlantic ridge system shown in great detail and the second (Figure 6) the Southern ocean.

Figure. 6: Optimally interpolated mean sea surface for the Southern Ocean with wavelengths greater than 2000 km. removed .

3. Conclusions

This paper has described a technique to derive a global high resolution MSS using a spectral optimal interpolation technique. The obvious advantages of this method over its spatial analogue are that of the computational speed of the FFT and the ability to solve for a global MSS rather than smaller regions. However, with the use of the FFT we inherit its shortfalls which in this case reduce the latitudinal extremities of ERS-1 data (though there may be a simple solution to this in the form of mirroring data) and the generation of land/ocean (and ice/ocean) boundary effects which have not been analysed in this paper. Our method, as with other studies, for deriving MSS noise characteristics clearly need further investigation as does the generation of the noise surface derived purely from altimetry. Despite these problems the method will open the path for the integration of mean sea surfaces from a range of altimeter missions from the past and the future (e.g. ENVISAT, JASON and GEOSAT follow-on) as spectral techniques cater for this with the intention of providing a global high resolution, high accuracy mean sea surface.

4. References

Basic and Rapp, 1992

T Basic and R H Rapp. Ocean wide prediction of gravity anomalies and sea surface heights using GEOS-3, Seasat and Geosat altimeter data and ETOPO5U bathymetric data. Report 416, Ohio State Univ., Dept. of Geod. Sci. and Surv., Columbus, 1992.

Bracewell, 1978

R Bracewell. The Fourier Transform and its applications. McGraw-Hill, New York, 2nd ed. edition, 1978.

Blanc et al., 1990

F Blanc S Houry P Mazzega and J F Minster. High-resolution, high-accuracy altimeter derived mean sea surface in the norwegian sea. Marine Geodesy, 14:57-76, 1990.

Carnochan et al., 1994

S Carnochan, P Moore, S Ehlers, C Lam, and P Woodworth. Improvement of radial positioning of ERS-1 through dual crossover analysis with TOPEX-Poseidon. SP 361, ESA, 1994.

Carnochan, 1996

S Carnochan. Orbit and altimetric corrections for the ERS satellites through the analysis of dual and single satellite crossovers. PhD thesis, University of Aston in Birmingham, 1996.

Le Traon et al., 1995

P Le Traon, P Gaspar, F Bouyssel, and H Makhmaraa. Use of Topex/Poseidon data to enhance ERS-1 data. J. Atm. Ocean. Tech., 12:161-170, 1995.

Mazzega and Houry, 1989

P Mazzega and S Houry. An experiment to invert seasat altimetry for Mediterranean and Black sea mean surfaces. Geophys. J., 96:259-272, 1989.

Rapp et al., 1991

R Rapp, Y M Wang, and N K Pavlis. The Ohio State 1991 geopotential and sea surface topography harmonic coefficient models. Report 410, Ohio State Univ., Dept. of Geod. Sci. and Surv., Columbus, 1991.

Rapp et al., 1994

R Rapp, Y Yuchan, and Y M Wang. Mean sea surface and geoid gradient comparisons with TOPEX altimeter data. Journal of Geophysical Research, 89:246757-24667, 1994.

Schwarz et al., 1990

K P Schwarz, M G Sideris, and R Forsberg. The use of FFT techniques in physical geodesy. Geophys. J. Int., 100:485-514, 1990.

Tsaoussi and Koblinksy, 1994

L S Tsaoussi and C J Koblinsky. An error covariance model for sea surface topography and velocity derived from topex/poseidon altimetry. Journal of Geophysical Research, 1994.

Wunsch and Zlotnicki, 1984

C Wunsch and V Zlotnicki. The accuracy of altimetric surfaces. Geophys. J. Royal Astron. Soc., 78:795-808, 1984.

Acknowledgements

The authors would like to thank the Natural Environment Research Council for their financial support for this study.

Keywords: ESA European Space Agency - Agence spatiale europeenne, observation de la terre, earth observation, satellite remote sensing, teledetection, geophysique, altimetrie, radar, chimique atmospherique, geophysics, altimetry, radar, atmospheric chemistry