Regional Seafloor Topography by Extended Kalman Filtering of Marine Gravity Data without Ship-Track Information

Abstract

An iterative Extended Kalman Filter (EKF) approach is proposed to recover a regional set of topographic heights composing an undersea volcanic mount by the successive combination of large numbers of gravity measurements at sea surface using altimetry satellite-derived grids and taking the error uncertainties into account. The integration of the non-linear Newtonian operators versus the radial and angular distances (and its first derivatives) enables the estimation process to accelerate and requires only few iterations, instead of summing Legendre polynomial series or using noise-degraded 2D-FFT decomposition. To show the effectiveness of the EKF approach, we apply it to the real case of the bathymetry around the Great Meteor seamount in the Atlantic Ocean by combining only geoid height/free-air anomaly datasets and using ship-track soundings as reference for validation. Topography of the Great Meteor seamounts structures are well-reconstructed, especially when regional compensation is considered. Best solution gives a RMS equal to 400 m with respect to the single beam depth observations and it is comparable to RMS obtained for ETOPO1 of about 365 m. Larger discrepancies are located in the seamount flanks due to missing high-resolution information for gradients. This approach can improve the knowledge of seafloor topography in regions where few echo-sounder measurements are available.

1. Introduction

Accurate knowledge of the Earth’s topography remains fundamental to understand surface processes. Detailed mapping of the sea floor topography (or bathymetry) is essential for studies in oceanography, biology, marine geology, natural disasters, airline crashes over ocean [1] habitat loss and marine resources. In oceanography, currents and tides physical characteristics are controlled by the shape of the seafloor as well as smaller-scale topographic features like seamounts.

So far our knowledge of seafloor topography remains imprecise. The traditional technique to map bathymetry is a tedious process made by research vessels equipped with single or multibeam echo-sounders [2,3]. Using this traditional technique only less than 18% of the worldwide seafloor has been determined with a spatial resolution of about 1 km [2]. However, reaching a complete coverage of the ocean seabed at an horizontal resolution of 100 m would take hundreds of years (as per one surveying boat) [4]. Moreover, because of the geometry of the multibeam echo-sounders, the resolution and the accuracy decreases with the increasing depth [5].

Radar altimeters aboard ERS-1 and GEOSAT satellites have measured the variations of the gravity field over all the oceans with high accuracy and medium spatial resolution. In March 1995, ERS-1 completed its dense (~8 km track spacing at the Equator) mapping of the sea surface topography between the latitudes of 81.5° North and South. The data have been combined and processed to form a global marine geoid or gravity grid [6,7]. In the wavelength band from 15 to 200 km, gravity anomaly variations are strongly correlated to sea floor topography, and thus, can be used for recovering topographic signals. The basic theory for predicting sea floor topography from satellite altimeter measurements is provided by [8], where the idea is to use sparse depth soundings to estimate the long wavelengths and short wavelengths predicted from downward continued satellite gravity measurements [9].

Bathymetric prediction was largely concerned with detecting uncharted seamounts from their gravitational signature [10,11,12]. Filters for modelling have also been designed to detect height and location of undersea structures along satellite radar altimeter tracks [13,14,15]. The first global bathymetric maps based on marine geoid (i.e., static sea surface) data were published by [16] and [17]. Inversion of these satellite data to estimate the bathymetry was initially performed along separate satellite tracks (1-D) using the so-called admittance method [8,18,19,20]. Moreover, [21] and [22] computed 1-D fitting of synthetic topographical profiles along satellite tracks. Denser coverage of the oceans by the GEOSAT (1985–1986) and ERS-1 (1994–1995) Geodetic Mission increased the horizontal resolution to better than ~10 km anywhere. These new generation of data sets (either used alone or merged with sparse echo soundings) have been inverted to compute 2-D estimates of the bathymetry [7,9,23,24,25,26,27]. Reference [28] has reviewed the computation of sea floor topography from altimetry and echo sounding data.

Most studies have been conducted in the Fourier domain where signal convolutions are replaced by simple spectral products. Reduction of the infinite power series of the density interface topography to gravity ratio, as proposed by [29], to the very first term has often been used for simplification. However, as discussed by [30] and later by [31] this truncation or linear assumption, severely underestimates the gravity to bathymetry ratio at wavelengths greater than 50 km. In particular, the gravity signature of the single-term predicted bathymetry computed with the Parker’s fully-developed series does not reproduce the gravity data from which that bathymetry was initially derived. Nevertheless, higher-order terms can be accounted in an iterative inversion process to estimate Fourier coefficients of the sea floor topography [23,31].

It remains not possible to combine gravity/geoid data to sparse in situ depth observations provide by ships for constraining the topographical solution in the Fourier domain, so that the two types of data can only be processed separately. Instead of being included in the inversion, the shipborne data are used to complete the bathymetry derived from satellite altimetry [7,9,27,31]. In the data preparation, the ship-track soundings have to be interpolated onto grids of regular sampling to be adapted for Fast Fourier Transform (FFT) techniques. The « spectral » methods consist of inverting spectra of the gridded gravity data directly by assuming no noise polluting these starting data, so the spectra of the predicted bathymetric signals include the contribution of this noise. In addition, these FFT methods require numerical stabilization [23,25] or low-pass filtering [7,9,27] to avoid or, at least, minimize the development of unrealistic short-wavelength undulations.

In the case of vast oceanic areas, the number of observations is large, the non-linear method developed by [32] and applied in the Atlantic Ocean by [33], around New Zealand [26] and finally globally [24] consists of inverting matrices of huge dimensions, and thus hardly applicable as itself. One of the strategies to decrease the dimensions of the problem, was to decompose the studied oceanic area into overlapping small rectangular cells covering the entire oceanic province, the inversion being made in each of these geographical cells of a few degree dimensions [24,26,33], but with the risk of spatial truncation errors (i.e., loss of signals energy) limited by including ship-track soundings for constraints.

Most recent studies have proposed to combine different gravity data using non-linear least square method [34]. But dealing with matrix of large dimensions is required.

Here, we propose an alternative procedure to manage the dimensions of the inversion problem better, by the use of the progressive Kalman Filtering-type assimilation of large numbers of geoid height or free-air anomaly observations. The ship-track soundings are used as a reference for validation and therefor are excluded from the inversion process.

The article is structured as follows: first, the derivation of the non-linear gravity operators for direct modelling of gravity observations from a set of topography heights is presented, as well as the Extended Kalman Filtering (EKF) procedure used for inversion of only gravity information. Secondly, the data used in this study are described (gravimetric observations, global bathymetric grids, ship-track soundings). Next, to demonstrate the feasibility of the proposed method, the seafloor topography recovery is made using simulated gravity data. Finally, our approach is applied to the well-known real case of the tabular Great Meteor (GM) guyot [W27°00′–W30°00′; N28°30′–N31°30′] in the mid-Atlantic Ocean [35,36].

2. Materials and Methods

2.1. Methodology

2.1.1. Forward Problem Using Radial Integrals

In this part, the theoretical expressions of the gravity quantities due to the presence of an interface of irregular shape are derived. A four-layer model is considered in this study (Figure 1) as proposed by [37]. We defined water density as ρ_w, seamount or load density as ρ_h, sediments density ρ_s, crust density as ρ_c and mantle density as ρ_m. This is a simplified model which allows considering different assumptions ρ_s = ρ_h or ρ_s = ρ_h = ρ_c.

For a set of discrete j = 1… M rock columns of constant density composing the topography measured at an external observation point of index i, the total potential and vertical gravity anomaly quantities are respectively [26,32]:

$$V_i={\displaystyle \sum _{j=1}^{M}V(r_i,h_j,{\phi }_{ij})=G\left({\rho }_h-{\rho }_w\right){\displaystyle \sum _{j=1}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\xi }(r_i,R+h_j,{\phi }_{ij})-\tilde{\xi }(r_i,R,{\phi }_{ij})\right]}}$$

(1)

and

$$\begin{array}{l}\Delta g_i={\displaystyle \sum _{j=1}^M\Delta g(r_i,h_j,h_i^P,{\phi }_{ij})}=2\pi G\left({\rho }_h-{\rho }_w\right)\left[\tilde{\zeta }(r_i,R+h_i^P,0)-\tilde{\zeta }(r_i,R,0)\right]\\ +G\left({\rho }_h-{\rho }_w\right){\displaystyle \sum _{\begin{array}{l}j=1\\ i\ne j\end{array}}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\zeta }(r_i,R+h_j,{\phi }_{ij})-\tilde{\zeta }(r_i,R+h_j^P,{\phi }_{ij})\right]}\end{array}$$

(2)

where r_i is the radial distance between the geocenter and the observation point P, h_i^P is the rock column height below P (i.e., the thickness of the considered spherical cap is h_i^P), R is the mean Earth’s radius (R = 6371 km), G is the gravitational constant (G~6.667 × 10⁻¹¹ m³∙kg⁻¹∙s⁻²), θ_j is the latitude, Δλ_j and Δθ_j are the dimensions of the mass element of index j in longitude and latitude, respectively, and φ_ij is the angular distance between the mass element and the observation point P. The functions $\xi$ and $\zeta$ are the estimated primitives of the radial functions of potential and gravity anomalies and were defined by [32] (Appendix A). These simple discrete forms are suitable to ease linear inversion of gravity data considering one layer without compensation effects due to flexure w of the lithosphere in the context of regional compensation (Figure 1).

Similarly as the computation of topographic effects, we can write the potential created by the regional compensation as [32]:

$$\begin{array}{l}{V}_i{}^C={\displaystyle \sum _{j=1}^{M}V(r_i,w_j,{\phi }_{ij})}\\ =G({\rho }_s-{\rho }_c){\displaystyle \sum _{j=1}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\xi }(r_i,R+w_j,{\phi }_{ij})-\tilde{\xi }(r_i,R,{\phi }_{ij})\right]}\\ +G({\rho }_c-{\rho }_m){\displaystyle \sum _{j=1}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\xi }(r_i,R_m+w_j,{\phi }_{ij})-\tilde{\xi }(r_i,R_m,{\phi }_{ij})\right]}\end{array}$$

(3)

and

$$\begin{array}{l}\Delta g_i{}^C={\displaystyle \sum _{j=1}^M\Delta g(r_i,w_j,w_i,{\phi }_{ij})}\\ =2\pi G({\rho }_s-{\rho }_c)\left[\tilde{\zeta }(r_i,R+w_i,0)-\tilde{\zeta }(r_i,R,0)\right]\\ +2\pi G({\rho }_c-{\rho }_m)\left[\tilde{\zeta }(r_i,R_m+w_i,0)-\tilde{\zeta }(r_i,R_m,0)\right]\\ +G({\rho }_c-{\rho }_s){\displaystyle \sum _{\begin{array}{l}j=1\\ i\ne j\end{array}}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\zeta }(r_i,R+w_{j,}{\phi }_{ij})-\tilde{\zeta }(r_i,R,{\phi }_{ij})\right]}\\ +G({\rho }_c-{\rho }_m){\displaystyle \sum _{\begin{array}{l}j=1\\ i\ne j\end{array}}^M\Delta {\lambda }_j\Delta {\theta }_j\cos {\theta }_j\left[\tilde{\zeta }(r_i,R_m+w_{j,}{\phi }_{ij})-\tilde{\zeta }(r_i,R_m,{\phi }_{ij})\right]}\end{array}$$

(4)

where w_j is the deflection of the lithosphere deduced from the load given by [38,39] (Appendix B) and depending upon the elastic thickness T_e.

The complete gravitational effects are computed by adding to Equations (1) and (2) the signals of deflection at the external observation point P called V_i^C and Δg_i^C. The geoid heights are derived from the disturbing potential T by considering Brun’s formula [40]:

$$N(r,\phi )=\frac{T(r,\phi )}{\gamma }$$

(5)

where γ is the gravity acceleration (~9.81 m∙s⁻²). It is worth noting that a regular discretization is not required.

2.1.2. Extended Kalman Filter for Recovery Seafloor Topography

Kalman Filtering (KF) is an optimal estimation algorithm that uses the series of sequential measurements obtained at different epochs and containing statistical noise [41,42]. It produces estimates of unknown parameters that tend to be more accurate than those based on a single set of measurements alone, by estimating joint probability distributions. At epoch number k, the observation equation is:

$$z_k=f(x_k)+{\nu }_k$$

(6)

where x_k is the vector that contains the unknown parameters, i.e., the topographic heights to be recovered, and f represents one of the non-linear operators given by Equations (1) and (2). ν_k is the observation noise which is assumed to be zero-mean Gaussian white noise [43] with covariance matrix R_k, such as ν_k ~ N(0,R_k).

In the Extended Kalman Filter (EFK), the state transition ensured by the operator f does not need to be linear, and that is the case for gravity data (see Equations (1) and (2)). Besides, F_k is defined as the Jacobian matrix formed by the partial derivatives of f versus each parameter x_k at epoch k, i.e.,:

$$F_k^i=\frac{\partial f(x_k)}{\partial x_k^i}$$

(7)

P_k is the a priori error covariance matrix of variable x_k. To describe the evolution of x_k in the estimation process, a model of propagation is built:

$$x_k=A_k{x}_{k-1}+{\omega }_{k-1}$$

(8)

where x_k−1 is the state of x_k at the previous time k − 1 and A_k is the projection matrix that is imposed as the identity matrix to simplify. ω_k−1 represents the process noise which is assumed to be drawn from a zero-mean normal distribution with covariance matrix Q_k−1, i.e., ω_k−1 ~ N(0,Q_k−1).

The two successive steps of an EKF iteration are: (1) Prediction of a priori state and associated error covariance (projection); (2) Innovation (or equivalently update or correction).

EKF works recursively, iteration-by-iteration, and the information from geoid heights or free-air anomaly data is cumulated to build and refine progressively the solution vector x_k containing the topographic heights. Schema of the iterative process is shown in Figure 2. The detailed prediction and innovation steps of an EKF iteration are:

Prediction of a priori state and associated error covariance (projection):

$$x_k=A_k{x}_{k-1}+{\omega }_{k-1}$$

(9)

$$P_k=A_k{P}_{k-1}{A}_k^T+Q_{k-1}$$

(10)

• Innovation (or equivalently update or correction):

Measurement pre-fit residuals:

$$y_k=z_k-f(x_k)$$

(11)

Pre-fit residual covariance:

$$B_k=R_k+F_k{P}_k{F}_k^T$$

(12)

Optimal Kalman gain:

$$K_k=P_k{F}_k^T{B}_k^{-1}$$

(13)

So that a posteriori state and covariance estimates are:

$$x_k^{new}=x_k+K_k{y}_k$$

(14)

$$P_k^{new}=(I-K_k{F}_k)P_k$$

(15)

where I represents the identity matrix. Once the new solution vector x_k^new and its error matrix P_k^new are obtained (Equations (14) and (15)), they are injected in Equation (9) to start a new iteration by considering a new flow of observations z_k.

2.2. Region of Study: Great Meteor Guyot

Discovered between 1925 and 1927 by the German research vessel Meteor, the Great Meteor (GM) seamount is located at 29°57′10.6′’N and 28°35′31.3′’W in the southern part of the Seewarte seamount chain on the Azores plateau (Figure 3). It is the largest guyot in the North Atlantic with a volume of 24,000 km³. The crust under the GM has an age of 85 million years deduced by the magnetic anomaly data (An34). It has a shallow and flat summit ranging between 150 and 300 m below the sea level suggesting it may have emerged in the past [36]. It is covered of 600 m-thick layer Isotopic K-Ar and Ar-Ar analysis of dredged basalts on the top and the flanks of the GM provided ages of 10.7 and 16.3 million years respectively. There are two small seamounts just southwest of GM lying at 3800 m (Figure 3): the Little Meteor (LM) and the Closs seamount (CL), as from the official GEBCO gazetteer [44].

Elastic thickness T_e in this region has been largely studied [37] and the obtained value is ~19 km.

2.3. Data

2.3.1. Geoid Height and Free-Air Anomaly Measurements

A regional grid of the geoid model EGM2008 computed by the GeoForschungsZentrum GFZ [45] is extracted over the GM seamount [N29°–N31°; W27°–W30°] by fixing the grid sampling, e.g., 1/25°. A full-degree resolution grid (maximum harmonic degree of 2190) i.e., with a maximal spatial resolution ~20 km as described by [40].

Vertical component of the gravity anomaly (free-air anomaly) was derived globally from radar altimetry and made available by the DTU 2015 [46] with a one minute (1′) spatial resolution. This spatial resolution is also reduced to 1/25° in our study. In fact, the correlation between gravity anomalies and topography generally decreases under 20 km of wavelengths due to satellite altimetry limitations [9]. Regional deep Earth signal, e.g., density contrast at Moho depth, is extracted from the geoid heights and free-air anomalies using Generic Mapping Tools (GMT) function “grdtrend” [47]. The geoid height and free-air anomaly grids to be inverted are plotted on Figure 4 and reveal the relevant gravitational signatures of the GM, LM and CL seamounts with seafloor depth of a few meters (or hundreds of mGal).

2.3.2. Bathymetric Grids and Ocean Depth Surveys

For estimating the performance of the methodology, our solution is compared to the global bathymetric models ETOPO1 [48,49] and GEBCO [50,51].

ETOPO 1 provides land and seafloor topography with a spatial resolution of 1 arc-minute. This model is built using global and digital datasets including in situ depth surveys.

The GEBCO grids are computed with a spatial resolution of 15 arc-second (~450 m) using bathymetric grids and surveys from many data contributors around the world. In this study the GEBCO2020 version is considered.

We have also used single beam depth observations provided by National Center for Environmental Information [52] available in the GM region.

3. Results

3.1. Application to Simulated Data and Recovery Test

In this section, we will present the recovery tests concerning our approximations, the methodology performance and the choice of inversion parameters. The accuracy of the radial integral approximation in the forward problem is evaluated. Secondly the independent recovery tests were made to choose the best inversion parameters (a priori variance of observations errors, a priori variance of the topography heights) for each gravity data type. Finally, we checked the contribution of a combined-data inversion of different observations, geoid heights and free-air anomalies, in order to improve the seafloor topography assessment.

3.1.1. Forward Problem: Validation of Non-Linear Operators

Simulated geoid height and gravity anomaly grids are used to check the validity and the robustness of the proposed approach of recovery bathymetry. The corresponding geoid height or free-air anomaly observed at sea surface are easily computed using the operator defined in Equations (1) and (2) on regular grids. These operators are calculated by the analytical radial integration and used in our inversion processing (see Equation (11)). Then, it is important to analyze the accuracy of our approximation with respect to other previous methods: (i) 2D FFT Parker method that corresponds to a planar approximation (GMT function “gravfft”, [53]; (ii) tesseroids approach using Taylor expansion proposed by [54] for estimated geoid heights or the software developed by [55] for estimating free-air anomalies.

The simple case of a 2500 m-amplitude conical non-compensated seamount of constant density 2700 kg∙m⁻³ lying under a seafloor depth of 5000 m and a water density of 1000 kg∙m⁻³ is considered to compute gravitational effects (geoid heights and free-air anomalies).

Figure 5 shows the results from different methodologies to estimate the gravitational signal of such a seafloor structure. We notice that the discrepancies between planar (Parker) and spherical approximations are significant, especially for the geoid height anomalies ~40 cm at the seamount peak and ~−40 cm at the seafloor level. Simulated anomalies are in very good agreement with gravitational effects estimated by other methodologies (differences are lesser than 1.5 cm and 1 mGal in absolute value). Large differences are located around the seamount summit.

3.1.2. Inverse Problem Analysis by EKF. Simulated Case

As the gravity observation errors are assumed to be spatially decorrelated, as previous studies [32]:

$$R_{k-1}={\sigma }_d^{2}I$$

(16)

where σ²_d represents the variance of the observation errors. Obviously, the real altimeter data are probably spatially correlated with no white noise. Further studies would need to determine a realistic structure of matrix R_k−1. σ_d should take a realistic value representing the true observation’s accuracy. The actual observations techniques ensure that geoid heights and free-air anomalies achieve an accuracy of ~5–10 cm [56] and 1–5 mGal [57], respectively, then the values of 10 cm and 5 mGal are adopted for σ_d.

The topographic solution to be recovered is defined as a set of juxtaposed rock columns of unknown heights. The starting a priori covariance of the topographic heights is of Hirvonen-type [58,59]:

$$P_{k-1=0}=\frac{{\sigma }_m^2}{1+{(\phi /{\phi }_0)}^2}$$

(17)

where σ²_m is the a priori variance of the topographic heights (i.e., the elements of the unknown vector x_k−1 to fit), ϕ is the angular distance between these topographic heights, and φ₀ is the angular length of decorrelation. The analysis of ϕ₀ by [32] has demonstrated that the best inversion results are obtained for ϕ₀ = 0.2°. Preliminary numerical tests have confirmed that topography recovery is optimal for this value of 0.2°. The elements of this matrix P_k will evolve during the EKF estimation process.

The accuracy of the method to the a priori σ_m needs to be checked. The a priori standard deviation of topographic heights, σ_m, varies between 50 and 3000 m. The absolute error of the difference between the input model (a 2500 m-amplitude conical non-compensated seamount) and the solution (estimated topography derived from simulated geoid height and gravity anomaly observations) are considered as an indicator of the accuracy to this parameter. Numerical tests show that, as expected, this growth parameter σ_m has to be large enough to accelerate the convergence to a stable solution in the case of “cold start” conditions (the first iteration topography is initiated as 0 m everywhere). A rule-of –thumb is that σ_m takes values greater than 2500 m (Figure 6), i.e., greater than the expected maximum height. Geoid height data inversion shows a quicker convergence than free-air anomaly inversion. In fact, at iteration 20, the errors obtained for the solution derived from the geoid heights inversion with respect to the reference bathymetry, range between 12 to 30 m as σ_m decrease from 3000 m to 500 m and the worst result is ~190 m of error for a σ_m = 50 m. For the case of the solution estimated using free-air anomalies, the same conclusion is found but errors are larger, for σ_m > 500 m the errors are ~170 m. The convergence stability is obtained around 10 iterations. In consequence, for the next recovery tests, the first 10 iterations are retained. A white noise is also added for simulating the influence of noisy observations with amplitudes of 10 cm for the geoid heights, and 5 mGal for free-air anomalies. Evidently, topography recovery quality decrease with the noise level. The topography of the seamount is found with an RMS of around 50 m after 10 iterations (Figure 7). Geoid height inversion is more affected by noise fluctuations and up to 2 iterations noise propagation remains important causing an unrealistic variations on the seabed of ±500 m (Figure 7b) when the seamount peak is well reconstructed with a maximum value of 2500 m. In order to prove the performance of EKF, the solution derived from the combined inversion geoid heights/free-air anomalies is computed.

Results show lower absolute error and RMS, which is interpreted as a better reconstruction of the seamount and a minimization of noise propagation after three iterations (RMS < 100 m and absolute error < 300 m; Figure 7a). The observed oscillations of the seafloor for combined inversion geoid heights/free-air anomalies have weaker amplitude ±100 m (Figure 7b) and the seamount peak is also well reconstructed with a minimum water depth of ~2500 m like for geoid height inversion. In the seafloor plain no oscillation is observed when no white noise is added.

Finally, we perform numerical tests including compensation effects due to flexure of the lithosphere (see Figure 1). The synthetic gravity observations are computed considering increasing T_e values (5 km, 10 km, and 20 km). We set constant Young modulus and Poisson’s ratio values of 7∙10¹⁰ N∙m⁻² and 0.25 [37], respectively. Moho depth is set to 6 km, as expected for a mature oceanic crust [60]. In the case of the no noisy data, we also obtain the typical decrease for recovery error versus the iteration number, where a couple of iterations are efficient enough to reach a stable solution. However, convergence to the reference bathymetry is slowed down in comparison with the non-compensated case, and the seamount amplitude is underestimated after 10 iterations. The RMS reaches 150–200 m for both types of input gravity data. When realistic amplitudes of noise (10 cm and 5 mGal respectively) are considered (Figure 8), free-air anomaly inversion gives less error of recovery (RMS range between 150–200 m) than geoid height inversion (RMS of 300 m). Numerical tests show that the combined inversion accelerates the convergence (Figure 8 a) with a reduction of absolute error of 25% in the case of T_e = 5 km, and of 15% in the case of T_e = 10 km with respect to only geoid height (or free-air anomaly) inversions. In the optimal case of T_e = 20 km, the RMS of the combined inversion and the only free-air anomaly inversion are similar, RMS of 150 m (Figure 8b). However, when the two types of data are combined, the profile of the estimated bathymetry show an improvement in recovering of the seamount summit (reference peak 2500 m; only free-air anomaly inversion peak at 2271 m; combined inversion gives a more realistic peak with a value of 2453 m close to the reference).

3.2. Inversion of a Real Case

For further demonstration, the method of EKF prediction is now challenged to recover the more complicated shape of a real undersea relief. For this purpose, the example of a huge tabular topography that has been widely discussed by many authors [35,36,61], so that available historical datasets, e.g., ship-track soundings, is used for validating the EKF solution obtained from geoid height and free-air anomaly inversion. Accurate depth data over GM seamount are not included in our seafloor topography computation, as they are used as a reference and compared to the solutions obtained by EKF. As estimated in previous works [36,60] load density is considered 2700 kg∙m⁻³ in this region, crust density is 2900 kg∙m⁻³ and mantle density is 3350 kg∙m⁻³. The sediment thickness is considered negligible [62]. To take the regional compensation effects into account the elastic thickness T_e is fixed at 19 km [37]. Mean ocean depth is assumed to be 4150 m in the Great Meteor region. A priori inversion parameters are set to σ_m = 4000 m, σ_d = 10 cm for geoid heights and σ_d = 5 mGal for free-air anomalies.

Figure 9 presents the bathymetric solutions provided by the EKF inversion (INV) using geoid heights (GEOID), free-air anomaly data (FAA) and as well as their combination (GEOID + FAA), respectively. The structures of the GM seamount are almost perfectly reconstructed if we consider that we used only satellite data and no ship-track information. This is very important for a worldwide reconstruction of the seafloor where few ship-track data exist.

Structures of other small undersea features are also well distinguished in the EKF solution and the LM and CL are identified (Figure 9). Besides the associated a posteriori maps of uncertainties, the solutions in Figure 9 reveal the development of edge effects especially in the geoid height inversion. Extension of these edge effects represents only twice or three times the sampling rate. For this reason, a larger region around the GM is considered for computing the topography. The rapid convergence of the geoid inversion (as see in Section 3.1.2) allows us to stop the estimation process at the second iteration in order to reduce edge effects. In fact, these effects quickly increase after two iterations giving a noisier solution. In addition, when the regional compensation is taken into account, the edge effects are propagated to the seabed, as this is well represented by the large formal errors obtained at seabed larger than 450 m.

In the case of the combined inversion, the EKF solution may benefit from long-wavelength component from geoid height anomalies (proportional to r⁻¹) and shorter wavelength component brought by the free-air anomalies that depend upon r⁻². Geoid heights are used to estimate the initial bathymetry after two iterations and free-air anomalies complete the solution up to 10 iterations.

The free-air anomaly inversion yields errors of 50–150 m at the seabed, and lesser than 30 m around the seamounts. The combined inversion of different gravity data type produces similar formal errors or a posteriori uncertainties. Figure 10 shows the difference between the EKF estimates bathymetry and global bathymetric models ETOPO-1 and GEBCO. Amplitudes of seafloor topography solutions are very realistic compared to the ones extracted from ETOPO-1 and GEBCO grids that include the sparse depth measurements and digitalized level curves. Weak differences are of 50 m or less located at the top of the GM. The combination of two types of data (geoid heights/free-air anomalies) leads to low errors at the intermediate spatial wavelengths. However, some important differences of ±400 m still exist at long-wavelengths indicating an inaccurate reference seabed depth variability. Larger discrepancies (>1400 m) are located at seamounts flanks that are gradient areas and come from missing information at short-wavelengths. It is well-known that these signals are better reconstructed when the echo-sounder ship-track observations are available and included in the estimate [32].

In addition, Figure 11 shows the difference with respect to ship-track bathymetric observations. The function “grdtrack” of GMT [47] is used to derive these differences. We notice the good reconstruction of the GM shape especially at the top of this seamount with differences of 40 m or less, representing 28% of the surface. As expected, larger differences are in the mountain flanks due to missing high resolution information and some over large differences are observed in the seafloor area and are due to the oscillations associated to the edge effects.

Table 1. Values of Global Bathymetric Models minus ship-track soundings and EKF solutions minus ship-track soundings. Units are meters.

Model Type Used to Compute Diff.		MAX	MIN	MEAN	RMS
Global grids	ETOPO1	4681.8	−1858.6	−12.5	364.5
	GEBCO	4634.1	−1866.83	13.5	323.5
EKF solution uncompensated case	GEOID	4069.8	−2328.8	−295.2	540.3
	FAA	3974.7	−2503.1	−298.8	543.3
	COMB	3809.4	−2472.0	−374.8	565.4
EKF solution Compensated case	GEOID	4649.9	−2168.9	22.7	402.9
	FAA	4924.5	−1892.5	200.9	445.8
	COMB	4586.7	−2137.2	−28.4	401.8

shows the statistical analysis of the differences between global bathymetric grids, the EKF seafloor topography solutions (uncompensated and compensated cases) and the measured ship-track soundings considered as a reference. The RMS strongly decreases when compensation effects are taken into account, i.e., 400–450 m, compared to the uncompensated case, i.e., 540–565 m. Moreover, the seafloor topography derived from a combined inversion of geoid heights and free-air anomalies gives slightly better results and remains comparable to the RMS of global bathymetric grids, especially for ETOPO-1. For GEBCO, the RMS is lower, about 323 m, because this grid already contains seafloor depth information at a higher spatial resolution.

Finally, EKF method enables to invert data onto irregularly-sampled grids and thus can reduce to half the number of the parameters (rocks columns) and consequently the computation time. Numerical tests have displayed a gain factor in computation time of about 3. We have used the irregularly-sampled grid shown in Figure 12 to invert the seafloor topography over the GM. Inversion results are shown in Figure 13. If we compare the difference between the combined inversion (GEOID+ FAA) regularly sampled versus bathymetry from ship-track soundings (Figure 11) and combined inversion irregularly sampled versus bathymetry from ship-track soundings (Figure 13a), the resulting maps have a similar shape of the differences which are concentrated in areas of strong topographic gradient. The RMS obtained for the regularly and irregularly sampled cases is similar with a value of ~401.9 m (see Table 1 and Figure 13b).

4. Discussion and Conclusions

We have proposed an improved approach based on iterative integration of continuous incoming gravity and geoid data by using an iterative EKF prediction to recover a regional set of topographic heights to characterize seafloor topography and apply it successfully for deriving oceanic areas without ship-track information. The matrix expression of the method offers the possibility of combining different types of data with various accuracy in input such as geoid heights and free-air anomalies. Compared to previous studies, the non-linear operators related to gravity quantities have been refined and numerically tested. The improvements of the inversion method that makes the inversion process faster are: (i) the use of a progressive integration of data by the EKF to deal with large numbers of observations, and (ii) the possibility to use irregular geographical grids (that could be Delaunay decomposition strategies) to reduce the number of parameters (or topographic heights) to determine and have the benefit of higher resolution in regions of real interest. Complementary geophysical information such as elastic thickness, sediment thickness and densities have been included in our EKF inversion scheme. The method has been tested by inverting simulated gravity data corresponding to an idealistic conical seamount, and the RMS is of 50 m if these starting data have a ~10 cm of uncertainty on the geoid heights or ~5 mGal in terms of free-air anomalies. Undersea reliefs in the GM region have been fully retrieved from gravity data with an RMS of 400 m. This error remains acceptable since it is comparable to the difference between global bathymetric models combined with echo-sounders observations. Large differences are located in regions of important gradient due to missing high resolution information. Long-wavelengths still need to be constrained by independent observations of seafloor depth which generate in our case an offset of ~50 m. The GM top structure and depth are well estimated with discrepancies of less than 40 m with respect to the depth measurements observed along ship-tracks. This discrepancies could come from the choice of the densities values, the determination of T_e or a combination of these factors.

In conclusion, we have shown that it is possible to reconstruct the topography of the seabed with a very comparable accuracy using only gravity data. This should allow us to reconstruct the topography of large ocean trenches and in areas that are not covered by acoustically-measured water depth.