1. Introduction
As an essential process in ocean exploration, bathymetric surveys play a fundamental role in developing marine resources, protecting the marine ecological environment, promoting marine science and technology, and safeguarding marine rights and interests. Acquiring accurate seafloor topography (ST) data is indispensable to the marine industry and is necessary for the protection, development, and management of oceans.
At present, bathymetric survey methods mainly include shipborne bathymetry, airborne lidar bathymetry (ALB), autonomous underwater vehicle (AUV), remotely operated vehicle (ROV), and bathymetry inversion using satellite altimetry gravity data. Shipborne bathymetry can yield highly accurate measurements but has several shortcomings, such as low-measurement efficiency and unevenly distributed results. Smith pointed out that deep-sea surveying and mapping would take more than 200 ship years and require tremendous workforce and material resources [
1]. ALB is mainly used to measure shallow water environments with high measurement efficiency and only has limited capabilities for deep-sea surveys and areas with serious beam propagation constraints. Altimetric satellites can provide ocean observations uninterruptedly under all weather conditions, providing real, complete, global, and precise sea surface data [
2,
3]. Global and regional gravity fields over the world’s oceans derived from satellite altimetry are now an efficient tool for modeling ST [
4]. This technology has been used for deep-sea research, exploration, and development. For example, GEBCO, supported by the Nippon Foundation-GEBCO Seabed 2030 project [
5], published the first global ST model (GEBCO_2019) with a grid size of 15″ in 2019. Scripps Institution of Oceanography (SIO), National Oceanic and Atmospheric Administration (NOAA), National Geospatial-Intelligence Agency (NGA), and other organizations released the global terrain model (SRTM15+ v2.0) with a grid size of 15″ in 2019 [
6]. Most of the ST data in GEBCO_2019 and SRTM15+ v2.0 were acquired using satellite altimetry gravity data.
Currently, ST inversion methods using satellite altimetry gravity data can be generally divided into space domain and spectral-domain approaches. ST inversion in the frequency domain is mainly based on the relationship between the mass body and gravity anomaly [
7], such as regression analysis [
8,
9,
10,
11,
12,
13], admittance function method [
14,
15,
16], simulated annealing method [
17,
18]. The admittance function method and regression analysis ignore the influence of nonlinear ST on the inversion results. Baudry et al. [
19,
20] and Fan et al. [
15] attempted to solve the problem of high-order ST inversion in the frequency domain. First proposed by Calmant [
21], the ST inversion method in the space domain can productively solve the problem of high-order ST inversion by using the nonlinear iterative least-square method. Calmant et al. used geoid height (GH) obtained from ERS-1, GEOSAT, and TOPEX/POSEIDON to build a global ST model in the space domain [
22]. Ramillien and Wright used sea surface gravity anomaly (GA) to restore ST in the New Zealand region (called RW99) in the space domain [
23].
Other studies have proposed using an approximate estimation to describe the formula in the space domain [
21,
22,
23]. The inversion parameters are mostly replaced by empirical values in the numerical analysis. For example, Ramillien and Wright set the inversion truncation wavelength to 500 km, the correlation length of the Hirvonen covariance function to 0.2°, and the terrain variance to 100 m–2000 m [
23]. For Calmant et al., the inversion truncation wavelength was set to 500 km, the correlation length of the Hirvonen covariance function to 20 km, and the terrain variance to 500 m [
22]. The inaccuracy of input parameters can seriously affect the quality of ST modeling and computational efficiency, which should be studied further. However, while some studies have carried out ST inversion in the space domain, none have used vertical gravity gradient (VGG) to predict bathymetry in the space domain. Compared to gravity anomaly, gravity gradients can be more sensitive to short bathymetric wavelengths [
24,
25], such as wavelengths ranging from 2 to 12 km [
24].
In this study, we evaluated the critical formulas in the space domain for ST inversion using the sea surface GH, GA, and VGG and described the design and implementation of ST inversion. We selected part of the Malaysia Airlines Flight MH370 search area as the experimental site. The inversion band of VGG/GA was determined through frequency domain coherence analysis [
26]. The effective initial value of the iteration of the Hirvonen covariance function was obtained by adopting the practical statistical formula of terrain covariance [
27]. We also used the GA and VGG as input data, along with sparse ship bathymetric data, to estimate ST. We then compared our models’ performance with the DTU18 (released by the Technical University of Denmark) and SIO V20.1 (released by Scripps Institution of Oceanography) bathymetry models.
2. Methods
The ST is mainly formed by (1) the cyclic process of continuous regeneration (at the mid-ocean ridge) and continuous consumption of the oceanic plates; (2) the volcanic activity within the oceanic plates; (3) the physical mechanism of cooling and contraction of the oceanic plates; and, (4) the lithosphere deformation under the action of terrain load.
In
Figure 1, the radius of the Earth is
R, the height of the ST is
b, the depth of the ocean is
h, and the reference depth of the ST is
d. When only considering single-crust structure and constant density of the crust and seawater, the disturbing potential (
T) at point (
P) caused by the mass cylinder (
M) is as follows [
28]:
where the subscript
j represents the location of the mass cylinder (
M);
u is the geocentric distance of the research point;
and
are the latitude and longitude of the research point;
G is the gravitational constant of the Earth;
v is the volume of the mass cylinder; and
is the difference between crust density and seawater density. The spatial distance (
l) between the integral element (
Q) in the mass cylinder (
M) and the research point (
P) can be expressed as
where
r is the geocentric distance of flow point
Q; and
is the angular distance between
P and
Q. Assuming the bottom area of column
M is
,
can be expressed as
where
and
are the spatial distance of the base area along the latitude and longitude; and
is the latitude at the center of the base area. Equation (1) can be expressed as
where
and
represent the lower and upper bounds of a cylinder (
M). These lower and upper bounds can be calculated using the expression:
2.1. Gravity Anomaly Calculation
The gravity disturbance produced by a cylinder (
M) at the research point (
P) can be obtained by calculating the radial direction derivative of the disturbing potential. According to Equation (4), the sea surface gravity disturbance at point
P can then be expressed as
At the sea surface, the gravity disturbance is close to the free-air gravity anomaly (FAGA), supposing that the free-air reduction is already done if one assumes the stationary sea surface to coincide with the geoid.
The GA at point
P is the sum of the influence of all the terrain cylinders around point
P, as shown in
Figure 2. So the GA is expressed as [
21,
22]:
where
i represents the ith column. The integral in Equation (7) can be expressed as
Equation (7) can then be expressed as
Equation (7) is the nonlinear function of ST (
b). Since the difference between the flow point topography (
b) and the radial (
r) is constant, it can also be considered that Equation (7) is the function of the radial (
r). Equation (7) can then be expressed as
where
and
are the ST and geocentric radius of the
i-point, respectively; and
f( ) represents a nonlinear operator. The result of linearizing (expand at
, this paper sets
) Equation (10) is
where
means that when
,
is infinitesimal. Equation (11) is further treated as follows:
Hence, the result of linearizing Equation (10) is
where
represents sea surface GA at research point;
represents
i-point ST; and
is the difference between the nonlinear and linear calculation results of GA at the research point. The linearization coefficient (
) of the nonlinear operator in Equation (10) can be calculated using the formula:
Equation (13) can then be transformed to the matrix form:
2.2. Vertical Gravity Gradient Calculation
The VGG is derived from calculating the radial direction derivative of GA. According to Equation (7), the VGG expression is as follows [
29,
30,
31]:
Equation (16) can be expressed as
Similar to GA calculation, the linearizing result of Equation (18) is
where
is the difference between the nonlinear and linear calculation results of VGG at the research point.
can be expressed as
2.3. Geoid Height Calculation
According to Brun’s formula, the GH produced by the cylinder at the research point is [
28]
where
is the mean gravity acceleration at the Earth’s surface. The GH at the research point can be expressed as
where
Equation (22) can then be expressed as
Similar to GA and VGG calculation, the linearizing result of Equation (22) is
where
is the difference between the nonlinear and linear calculation results of GH at the research point.
can be expressed as
2.4. Seafloor Topography Inversion
Assuming the initial observation value of ST at
i-point is
, the gravity observation value at point
P is
(
stands for GA or VGG or GH), and there are m gravity observations in the study area, the error equation is expressed as
Equation (27) can be written as
where
is the error matrix;
is the coefficient matrix;
is the terrain adjustment values;
is the terrain initial value matrix;
is ST observation values and gravity values; and,
is the nonlinear and linear difference vector of the terrain calculation. Adopting the least square adjustment method, the result of ST adjustment is
where
is the weight matrix of observations. The final result of the terrain adjustment is
The terrain covariance matrix and the observation covariance matrix can be expressed as follows:
Introducing the Newton iteration method, when the error of the observation value is taken into account, the iteration result of Equation (30) is
where
r and
s are the estimated location and measurement data;
and
are the vector at the
n iteration and
n−1 iteration;
is the initial value vector of ST; and
is a vector of the ST and sea surface gravity data. Matrix
accounts for the uncertainties of the measurement data, while the covariance matrix
accounts for uncertainties in the ST. If the measurement is regarded as independent of each other, the measurement error matrix is given by
where
is the unit matrix; and
is the variance of the measurement. The standard deviation (SD) of satellite altimetry GA is generally about 2–5 mGal [
24,
32]. The root mean square (RMS) of shipborne gravity is approximately 1–3 mGal after considering crossover adjustments [
33,
34]. Considering the measurement error, navigation conditions, and the high-frequency noise, the SD of acoustic bathymetry results is usually 300 m [
23]. We modeled the uncertainty associated with a priori solution through the Hirvonen function [
21,
22,
23]:
where
is a priori uncertainty between topographic heights; and
is the spherical distance between the locations
r and
r’. Ramillien and Wright set 300–2000 m for the parameter [
23], while Calmant et al. set 500 m for the parameter. For the topography correlation length (
) [
22], Calmant et al., set the value to 0.25° and 20 km [
21,
22]. To make Equation (32) converge rapidly, the terrain covariance of each iteration can be adjusted and optimized according to the covariance propagation law:
Based on the above discussion, the flow chart of ST inversion by a nonlinear iterative least square method in the space domain is drawn, as shown in
Figure 3.
4. Discussion
DTU18 bathymetry model (hereinafter referred to as DTU18) and SIO V20.1 bathymetry model (hereinafter referred to as SIO V20.1) can be derived by combining satellite altimetry and multibeam data. In this study, we evaluated and compared the DTU18 and SIO V20.1, and the results are shown in
Figure 11a,b.
The ST models generated by the different approaches (
Figure 10 and
Figure 11) were assessed and compared. In general, we found that BAT_GA_ ILS, BAT_VGG_ ILS, and BAT_GA_VGG_ ILS are in close agreement with the DTU18 and SIO V20.1 and can reflect the topographic and geomorphological characteristics of the study area. We also found that BAT_GA_ILS, BAT_VGG_ILS, BAT_GA_VGG_ILS, and SIO V20.1 provided better terrain characterization of the study area than the DTU18. For example, DTU18 is smoother than other ST models in the red box range, having a comparatively weaker ability to present terrain details. In the highlighted portion of the image, three distinct topographic uplifts can be seen in the southwest, recognizable topographic uplifts in the southeast, and noticeable topographic changes in the northern seamount area four images (i.e., BAT_GA_ILS, BAT_VGG_ILS, BAT_GA_VGG_ILS, and SIO V20.1).
In comparison, DTU18 has smooth terrain, and no topographic details, such as topographic uplifts, can be seen in the same locations. This suggests that satellite altimetry gravity data and nonlinear iterative least square method can be used in generating the ST model. Using SIO V20.1 as reference, the topographic features of BAT_GA_ILS, BAT_VGG_ILS, and BAT_ GA_VGG_ ILS are largely similar with SIO V20.1. However, careful inspection of the images shows that BAT_VGG_ILS is better than BAT_GA_ILS and BAT_GA_VGG_ILS. BAT_GA_ILS and BAT_GA_VGG_ILS show the terrain to be relatively smooth. For example, in the seamount area (101.9° E, 33.2° S), BAT_VGG_ILS shows more prominent terrain details, while BAT_GA_ILS and BAT_GA_VGG_ILS show smoother textures. In the upper part of the red box, BAT_GA_ILS and BAT_GA_VGG_ILS show smooth features at the terrain change, while BAT_VGG_ILS separates the terrain features more distinctly. BAT_VGG_ILS contains more abundant high-frequency terrain information than BAT_GA_ILS and BAT_GA_VGG_ILS. VGG contains higher frequencies than GA from coherence analysis between the VGG/GA and ST.
The statistical results of the five ST models in the effective inversion area (within the highlighted region) are shown in
Table 3.
According to the statistical results of the ST models, the maximum sea depth is about 7000 m, and the average sea depth ranges from 4500 m to 4800 m. The study area’s topography fluctuates considerably, and the SD of ST models is about 900 m. The complex topography brings particular challenges for the MH370 debris search. SM data in the effective inversion area (a total of 605 checkpoints) was used as the external check reference to evaluate the accuracy of each ST model. In the original
(
,
) area, there were 1889 points, while in the final statistical
(
,
) region, the number of checkpoints was 695. Sea depth values of the ST model were interpolated using the bilinear interpolation method, and the difference between interpolated and actual sea depth values at checkpoints was calculated. The statistical results of the validation assessment are summarized in
Table 4. The relative accuracy (RA) is defined as the ratio of RMS to the absolute value of the average sea depth in the study area.
The correlation coefficients for BAT_GA_ILS, BAT_VGG_ILS, and BAT_GA_VGG_ILS are all above 0.98, which is similar to the correlation coefficient for SIO V20.1 and higher than for DTU18. The checking accuracy (RMS) for the SIO V20.1 is 51.10 m, which is the highest among the ST models. The checking accuracy for BAT_VGG_ILS is much higher than for BAT_GA_ILS and BAT_GA_VGG_ILS, while the checking accuracy for BAT_GA_VGG_ILS is slightly lower than for BAT_GA_ILS.
The two gravity models used in the paper were derived from satellite altimetry data GA and VGG. The data sources for the two gravity models are almost the same, so we believe that the two models are not completely independent. The fusion of satellite altimetry GA and VGG did not yield excellent results since (1) they do not contain the same frequency content, and (2) high-frequency information is somehow suppressed (a larger SD value is shown). However, the results of the study are still largely meaningful. Given gravity data from different sources (e.g., shipborne, airborne, gravity satellite, etc.), nonlinear iterative least square method can be used to estimate ST by combining multiple gravity data. In this study, we explored the feasibility of combining GA and VGG to model ST in the space domain. No previous study has used the VGG and GA together to predict bathymetry in the space domain.
The checking accuracy for BAT_GA_ILS, BAT_VGG_ILS, and BAT_GA_VGG_ILS are 2.49 times, 4.54 times, and 2.45 times the value of the DTU18 model. The RA for the five ST models is smaller than 10% in the study area, and that of our inversion models is smaller than 4%. The RMS of BAT_VGG_ILS, which was constructed using VGG, is smaller than 100 m. The RA is close to 2%, while the RC of SIO V20.1 is close to 1%. Analyzing and comparing the values of the ST model, the mean and median values for SIO V20.1 are within 2 m, while the absolute mean values for BAT_GA_ILS and BAT_GA_VGG_ILS are about 25 m. The absolute mean value for BAT_VGG_ILS is about 45 m, while for DTU18 is about 270 m. These values indicate systematic and other types of error present in BAT_GA_ILS, BAT_VGG_ILS, BAT_GA_VGG_ILS, and DTU18. These errors may have been handled in the modeling process for SIO V20.1. For example, to reduce (or eliminate) noise and downward continuation factors, high-frequency information has been filtered in advance [
12]. Moreover, systematic error in the ST model has been corrected through the SM data (SIO V20.1 was constrained almost by MH370 bathymetric surveys), and the error of input modeling data has been strictly processed.
Figure 12 presents the ratio of the checkpoints within a specific range to the total number of checkpoints. Based on the figure, the performance of SIO V20.1 is better than the other four models. In the 200 m range, the ratio of checkpoints for BAT_VGG_ILS is close to 95%; about 80% for BAT_GA_ILS and BAT_GA_VGG_ILS, and less than 50% for DTU18. While the curves in BAT_GA_ILS and BAT_GA_VGG_ILS roughly coincide, BAT_GA_ILS performed slightly better than BAT_GA_VGG_ILS.
The statistical analysis of the relative error for the five ST models is summarized in
Table 5. The relative error is defined as the ratio of the ST model checking difference to the SM data at the checkpoint. The results show that the SD of relative error is smallest for SIO V20.1 and largest for DTU18, while the values for BAT_GA_ILS and BAT_GA_VGG_ILS are similar. The maximum relative errors in BAT_GA_ILS, BAT_VGG_ILS, BAT_GA_ILS, and BAT_GA_VGG_ILS are less than 20%. In contrast, the value for DTU18 is close to 40%, which suggests that our ST models are better than the DTU18 model for the given study area. To analyze the performance of the ST model in different sea depths, the proportion of relative error under different sea depth conditions is shown in
Figure 13.
To analyze the performance of the ST model at different sea depths,
Figure 13 shows the proportion of relative error under different sea depth conditions. The findings suggest that with the increase in sea depth, the relative errors of the five ST models are decreasing and converging. The relative errors for BAT_GA_ILS, BAT_VGG_ILS, BAT_GA_VGG_ILS, and SIO V20.1 fluctuate near the zero value, while DTU18 has noticeable systematic errors. The relative error variation for BAT_GA_VGG_ILS is very similar to BAT_GA_ILS. Comparing
Figure 13b with
Figure 13a,c, the relative error of BAT_VGG_ILS is mainly concentrated around the zero value, indicating an inversion advantage for VGG. This is an interesting finding which may need further investigation in future research.
To further analyze the effectiveness of the ST model in different seabed topographic configurations, the spatial distribution of the relative error for each model was generated and is presented in
Figure 14. Large relative errors appear mainly in areas near the sea mountain where the terrain fluctuates sharply, while small relative errors can be found over generally flat terrain. Possible explanations for such findings are as follows: (1) The actual terrain density may change due to the influence of the sedimentary layer, which could have affected the calculation results; (2) The seafloor geology changes considerably over highly rugged terrain, which can easily cause isostasy and anomalies in the sea surface gravity data (i.e., terrain signal); and (3) For areas with steep and undulating terrain, the inherent attenuation of the gravity field signal does not provide reliable results.
DTU18 had significantly more areas with large relative error than the other four ST models, while SIO V20.1 had the least. Furthermore, areas with high relative error were considerably greater in BAT_GA_ILS and BAT_GA_VGG_ILS than in BAT_VGG_ILS. These findings suggest that SIO V20.1 performs comparatively better compared to any other model. One possible reason is that the SIO V20.1 bathymetric model was constructed using data from the MH370 bathymetric surveys, including checkpoint data. Moreover, SIO V20.1 was constrained by the SM data after model construction [
12,
13], while the other models were not. Smith et al. published a paper suggesting [
13]
Their method was improved by adding a constrain-propagation step: grid cells constrained by data were set to the median of data values in the cell, and then a finite-difference interpolation routine was used to perturb neighboring esti-mated values toward the constrained values. Thus, in well-surveyed areas, the accuracy and resolution depended only on the grid spacing and the quality of the constraints.
Table 6 shows the statistical analysis when our models are compared with SIO V20.1. The SD for BAT_GA_ILS–SIO V20.1 and for BAT_GA_VGG_ILS–SIO V20.1 is about 200 m. When large variations are ignored, the systematic deviations for BAT_GA_ILS–SIO V20.1 and BAT_GA_VGG_ILS–SIO V20.1 go down to about 30 m. For BAT_VGG_ILS–SIO V20.1, the SD is about 120 m, and the maximum difference is about 690 m. This means that BAT_VGG_ILS is closer to SIO V20.1.
The statistical summary of differences within the different ranges is shown in
Figure 15. Assuming that the strict regulation of the normal distribution is ignored, the normal distribution of the histogram is analyzed, and the difference distribution curve is generated for each model (shown in red line). The results show that most of the disparities are kept within 200 m. For BAT_GA_ILS and BAT_GA_VGG_ILS, 86% of the deviations with SIO V20.1 are within 300 m. For BAT_VGG_ILS, 90% of the difference with SIO V20.1 is less than 200 m.
5. Conclusions
This paper thoroughly discussed the processes and rigorous analytical formulas in the space domain for ST inversion using sea surface GH, GA, and VGG. We used the Malaysia Airlines Flight MH370 Search Area as the study site to test the effectiveness of the proposed algorithm. Considering the VGG/GA in the study area as external input, ST models were constructed through nonlinear iterative least-square method in the space domain. The inversion results were then compared with DTU18 and SIO V20.1. Based on the results, the following conclusions were made:
(1) Taking the 20–200 km inversion waveband ST as input data, the terrain covariance was calculated using the empirical formula for terrain covariance. The statistical results show that the terrain variance is 0.6365 when the distance is zero, and the correlation length is 10.5′.
(2) ST inversion using nonlinear iterative least square method tends to converge after two iterations in the study area. This is true whether the sea surface GA/VGG is used alone, or whether GA and VGG are used together. Therefore, we suggest that the number of iterations in ST inversion should be set to 2 when considering calculation efficiency.
(3) The checking accuracy of the ST model using sea surface VGG was highest, which is about twice the checking accuracy using sea surface GA or using combined GA and VGG data. The checking accuracy of the ST model generated using combined GA and VGG was not much higher than the accuracy of the model generated using only GA.
(4) The relative accuracy of the ST model constructed in this study was higher than 4%, while the relative accuracy of DTU18 was 9.27%. Within the 200 m range, the percentage of checkpoints for BAT_VGG_ILS was close to 95%, about 80% for BAT_GA_ILS and BAT_GA_VGG_ILS, and less than 50% for DTU18.
Note that while the use of nonlinear iterative least square method to build the ST model can effectively alleviate the effects of high order ST, the solution process involves large-scale matrix inversion and other problems. The calculation efficiency is lower compared with the ST frequency domain inversion method, which can seriously affect the processing time for ST modeling. For future endeavors, we plan to introduce parallel calculations and explore more optimal algorithms to improve the ST modeling efficiency in the space domain.