Next Article in Journal
An Adversarial Generative Network for Crop Classification from Remote Sensing Timeseries Images
Previous Article in Journal
Analysis of a Bistatic Ground-Based Synthetic Aperture Radar System and Indoor Experiments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Seafloor Topography Estimation from Gravity Anomaly and Vertical Gravity Gradient Using Nonlinear Iterative Least Square Method

1
Information Engineering University, Zhengzhou 450001, China
2
School of Geodesy and Geomatics, Wuhan University, Wuhan 430074, China
3
Institute of Geophysics, PGMF and School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
4
School of Land Science and Technology, China University of Geosciences (Beijing), Beijing 100000, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(1), 64; https://doi.org/10.3390/rs13010064
Submission received: 21 November 2020 / Revised: 22 December 2020 / Accepted: 23 December 2020 / Published: 26 December 2020

Abstract

:
Currently, seafloor topography inversion based on satellite altimetry gravity data provides the principal means to predict the global seafloor topography. Researchers often use sea surface geoid height or gravity anomaly to predict sea depth in the space domain. In this paper, a comprehensive discussion on seafloor topography inversion formulas in the space domain is presented using sea surface geoid height, gravity anomaly and introduces an approach that uses vertical gravity gradient. This would be the first study to estimate seafloor topography by vertical gravity gradient in the space domain. Further, a nonlinear iterative least-square inversion process is discussed. Using the search area for the Malaysia Airlines Flight MH370 as study site, we used the DTU17 gravity anomaly model and SIO V29.1 vertical gravity gradient to generate the seafloor topography. The results of the proposed bathymetric models were analyzed and compared with the DTU18 and SIO V20.1 bathymetric models. The experimental results show that the gravity anomaly and vertical gravity gradient in the study area are strongly correlated with the seafloor topography in the 20–200 km wavelength range. The optimal initial iteration values for seafloor topography variance and correlation length are 0.6365 km 2 and 10.5′, respectively. Shipborne measurements from SONAR data were used as external checkpoints to evaluate the bathymetric models. The results show that the RMS for BAT_VGG_ILS (inversion model constructed by vertical gravity gradient) is smaller than for BAT_GA_ILS (inversion model constructed by gravity anomaly) and BAT_GA_VGG_ILS (inversion model constructed by gravity anomaly and vertical gravity gradient). The relative accuracy of the DTU18 bathymetry model was 9.27%, while the relative accuracy of the proposed seafloor models was higher than 4%. Within the 200 m difference range, the proportion 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 the DTU18. The results show that the nonlinear iterative least square method in the space domain is feasible.

Graphical Abstract

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]:
T j ( u , φ p , λ p ) | u = R = G Δ ρ v 1 l d v
where the subscript j represents the location of the mass cylinder (M); u is the geocentric distance of the research point; φ p and λ p 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
l = r 2 + u 2 2 u r cos ( ψ ) | u = R
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
Δ Ω = r 2 Δ λ Δ φ cos φ
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
T j ( u , φ , λ ) | u = R = G Δ ρ Δ Ω R 1 R 2 1 l d r
where R 1 and R 2 represent the lower and upper bounds of a cylinder (M). These lower and upper bounds can be calculated using the expression:
{ R 1 = R d R 2 = R 1 + b = R h

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
Δ g j ( u , φ p , λ p ) | u = R = T u | u = R = G Δ ρ ( Δ λ Δ φ cos φ ) R 1 R 2 r 2 ( R r cos ψ ) l 3 d r
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]:
Δ g ( u , φ p , λ p ) | u = R = G Δ ρ Δ λ Δ φ i = 1 n ( cos φ i R 1 R 2 r i 2 ( R r i cos ψ i ) l i 3 d r )
where i represents the ith column. The integral in Equation (7) can be expressed as
ζ ( r i ) = r i 2 ( R r i cos ψ i ) l i 3 d r = { R 2 R r i + 2 R ln ( R r i ) ( R r i ) ; ψ i = 0 ( r i 2 + 3 R 2 ) cos ψ i + R r i ( 2 + 3 cos ( 2 ψ i ) ) l i 1 2 R ( 1 + 3 cos ( 2 ψ i ) ) ln ( r R cos ψ i + l ) ; ψ i 0
Equation (7) can then be expressed as
Δ g ( u , φ p , λ p ) | u = R = G Δ ρ Δ λ Δ φ i = 1 n ( cos φ i · ζ ( r i ) | R 1 R 2 )
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
{ Δ g ( u , φ p , λ p ) | u = R = i = 1 n f ( b i ) = i = 1 n f ( r i ) f ( r i ) = G Δ ρ Δ λ Δ φ ( cos φ i R 1 R 2 r i 2 ( R r i cos ψ i ) l i 3 d r )
where b i and r i are the ST and geocentric radius of the i-point, respectively; and f( ) represents a nonlinear operator. The result of linearizing (expand at r 0 , this paper sets r 0 = R d ) Equation (10) is
f ( r i ) = f ( r 0 ) + f ( r 0 ) · ( r i r 0 ) + o ( r i r 0 ) 2
where o ( r r 0 ) 2 means that when ( r i r 0 ) 0 , ( r i r 0 ) 2 is infinitesimal. Equation (11) is further treated as follows:
f ( r i ) = f ( r 0 ) r i + ( f ( r 0 ) f ( r 0 ) r 0 ) = a r i r i + ξ r i
Hence, the result of linearizing Equation (10) is
Δ g ( u , φ p , λ p ) | u = R = i = 1 n ( a i · b i ) + ξ p
where Δ g ( u , φ p , λ p ) | u = R represents sea surface GA at research point; b i represents i-point ST; and ξ p is the difference between the nonlinear and linear calculation results of GA at the research point. The linearization coefficient ( a i ) of the nonlinear operator in Equation (10) can be calculated using the formula:
a i = f ( r i ) = G Δ ρ Δ λ Δ φ ( cos φ i · r i 2 ( R r i cos ψ i ) l i 3 )
Equation (13) can then be transformed to the matrix form:
Δ g ( u , φ p , λ p ) | u = R = ( a 1 , a 2 , , a n ) · ( b 1 , b 2 , , b n ) T + ξ p

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]:
Δ g z ( u , φ p , λ p ) | u = R = 2 T u 2 | u = R = G Δ ρ Δ λ Δ φ i = 1 n ( cos φ i R 1 R 2 ( 3 r i 2 ( R r i cos ψ ) 2 l i 5 + r i 2 l i 3 ) d r )
Then,
ς z ( r i ) = ( 3 r i 2 ( R r i cos ψ i ) 2 l i 5 + r i 2 l i 3 ) d r = { ( 3 R 4 r i ) R ( R r i ) 2 + 2 ln ( R r ) ; ψ = 0 R ( 11 r i 2 + 3 R 2 ) cos ψ i + r i ( 2 r i 2 + 5 R 2 + ( 4 r i 2 + 6 R 2 ) cos ( 2 ψ i ) 3 r R cos ( 3 ψ i ) ) l i 3 1 2 ( 1 + 3 cos ( 2 ψ i ) ) ln ( r i + l R cos ψ i ) ; ψ 0
Equation (16) can be expressed as
Δ g z ( u , φ p , λ p ) | u = R = 2 T u 2 | u = R = G Δ ρ Δ λ Δ φ i = 1 n ( cos φ i · ς z ( r i ) | R 1 R 2 )
Similar to GA calculation, the linearizing result of Equation (18) is
Δ g z ( u , φ p , λ p ) | u = R = i = 1 n ( a z i · b i ) + ξ z p
where ξ z p is the difference between the nonlinear and linear calculation results of VGG at the research point. a z i can be expressed as
a z i = G Δ ρ Δ λ Δ φ ( cos φ i · 3 r i 2 ( R r i cos ψ ) 2 l i 5 + r i 2 l i 3 )

2.3. Geoid Height Calculation

According to Brun’s formula, the GH produced by the cylinder at the research point is [28]
N j ( u , φ p , λ p ) | u = R = T j ( u , φ p , λ p ) | u = R γ = G Δ ρ Δ Ω R 1 R 2 1 l d r γ
where γ is the mean gravity acceleration at the Earth’s surface. The GH at the research point can be expressed as
N ( u , φ p , λ p ) | u = R = G Δ ρ Δ λ Δ φ γ i = 1 n ( cos φ i R 1 R 2 r i 2 l i d r )
where
ζ N ( r i ) = r i 2 l i d r = { 1 2 ( R r i ) ( 3 R + r i ) R 2 ln ( R r i ) ; ψ i = 0 1 4 ( 2 ( r i + 3 R cos ψ i ) l i + R 2 ( 1 + 3 cos ( 2 ψ i ) ) ln ( r R cos ψ i + l i ) ) ; ψ i 0
Equation (22) can then be expressed as
N ( u , φ p , λ p ) | u = R = G Δ ρ Δ λ Δ φ γ i = 1 n ( cos φ i · ζ N ( r i ) | R 1 R 2 )
Similar to GA and VGG calculation, the linearizing result of Equation (22) is
N ( u , φ p , λ p ) | u = R = i = 1 n ( a N i · b i ) + ξ N p
where ξ N p is the difference between the nonlinear and linear calculation results of GH at the research point. a N i can be expressed as
a N i = G Δ ρ Δ λ Δ φ ( cos φ i · r i 2 l i )

2.4. Seafloor Topography Inversion

Assuming the initial observation value of ST at i-point is b i , 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
[ v 1 v 2 v n v τ 1 v τ 2 v τ m ] = [ 1 0 0 0 1 0 0 a 1 τ 1 a 1 τ 2 a 1 τ m 0 a 2 τ 1 a 2 τ 2 a 2 τ m 1 a n τ 1 a n τ 2 a n τ m ] · [ δ b ^ 1 δ b ^ 2 δ b ^ n ] + [ 1 0 0 0 1 0 0 a 1 τ 1 a 1 τ 2 a 1 τ m 0 a 2 τ 1 a 2 τ 2 a 2 τ m 1 a n τ 1 a n τ 2 a n τ m ] · [ b 1 0 b 2 0 b n 0 ] [ b 1 b 2 b n τ 1 τ 2 τ m ] + [ 0 0 0 ξ 1 ξ 2 ξ m ]
Equation (27) can be written as
v = A · δ b ^ ( L A · b 0 ξ )
where v is the error matrix; A is the coefficient matrix; δ b ^ is the terrain adjustment values; b 0 is the terrain initial value matrix; L 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
{ δ b ^ = ( A T P L A ) 1 · ( A T P L l ) l = ( L A · b 0 ξ )
where P L is the weight matrix of observations. The final result of the terrain adjustment is
b ^ = b 0 + δ b ^
The terrain covariance matrix and the observation covariance matrix can be expressed as follows:
{ C b ^ b ^ = ( A T P L A ) 1 C L L = A C b ^ b ^ A T
Introducing the Newton iteration method, when the error of the observation value is taken into account, the iteration result of Equation (30) is
b n ( r ) = b 0 ( r ) + C b b A n T ( A n C b b A n T + E s s ) 1 · { L ( s ) A n b 0 ( r ) f ( b n 1 ( r ) ) + A n b n 1 ( r ) }
where r and s are the estimated location and measurement data; b n ( r ) and b n 1 ( r ) are the vector at the n iteration and n−1 iteration; b 0 ( r ) is the initial value vector of ST; and L ( s ) is a vector of the ST and sea surface gravity data. Matrix E s s accounts for the uncertainties of the measurement data, while the covariance matrix C b b accounts for uncertainties in the ST. If the measurement is regarded as independent of each other, the measurement error matrix is given by
E s s = σ s 2 · I
where I is the unit matrix; and σ s 2 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]:
C b b = C 0 1 + ( ψ ( r , r ) ψ 0 ) 2
where C 0 is a priori uncertainty between topographic heights; and ψ ( r , r ) 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 ( ψ 0 ) [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:
C b b = C b b A n T ( A n C b b A n T + E s s ) 1 A n C b b
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.

3. Results

3.1. Data Preparation and Pre-Processing

A 2 ° × 2 ° ( 33 °   S 35 °   S , 101 °   E 103 °   E ) area of the MH370 suspected crash site was chosen as the research area for this case study. The data sources are as follows: (1) Satellite altimetry GA and VGG were obtained from the DTU17 model [3] and from SIO V29.1. The satellite altimetry-deduced GA and VGG had 1 × 1 spatial resolution, shown in Figure 4a,b. (2) Shipborne measurements (SM) from SONAR data used in the research are published by the Australian government (http://www.ga.gov.au/about/projects/marine/mh370-data-release). A total of 7556 SM data (about four-fifths of the total data) were randomly selected as control points to compensate for non-inversion waveband ST. The remaining one-fifth of the SM points were used as external checkpoints to evaluate the ST model. The distribution of control points and checkpoints in the study area is shown in Figure 4c, where the blue and red dots represent control points and checkpoints, respectively. The statistical analysis of the gravity and SM data are summarized in Table 1.
The average sea depth is 4702.02 m (SD = 783.46 m), while the minimum sea depth is 6971.99 m. These values suggest that the study area’s terrain fluctuates sharply and that improper search equipment selection may result in large losses (some search equipment have to be near the seabed during detection). The mean GA is 28 mGal (SD = 43 mGal), while the mean VGG is −1.5 Eotvos (SD = 38 Eotvos). The discrete control points are placed into a 1′× 1′ sea depth grid using the Generic Mapping Tools (GMT). Using the gridded VGG/GA results, and ST as input data, the coherence in waveband information between the input data was obtained through coherence analysis technology [25]. The results are shown in Figure 5, where the red points show the GA and ST wavelength consistency, while the blue points present the VGG and ST wavelength consistency. Fourth-degree polynomial regression was then used to obtain the fitting relation between wavelength and coherence for the two sets.
The coherence results in Figure 5 show that VGG has a stronger coherence with the ST in the high-frequency segment, as compared with the GA. Since the signal-to-noise ratio of the input data is 1:1 when the coherence is at 0.5, we think that the result of coherence is significant when the value is above 0.5. Based on coherence analysis (Figure 5), strong correlation (where coherence is greater than 0.5) between GA/ST wavebands ranges from 20 to 230 km, while for VGG/ST, the range is between 20 and 200 km. Thus, the range 20–200 km was chosen for the ST inversion band. The band-pass filter was used to filter the VGG/GA input data, and the final results of the gravity data at the inversion band are shown in Table 2 and Figure 6.
In Figure 6, a negative belt can be seen running through the seafloor along the SE-NW direction, and two positive zones are shown in the northwest portion. These observations may suggest that the study area has two topographic uplifts in the northwest and a trench running along the SE–NW direction.

3.2. Seafloor Topography Covariance Initial Value Determination

To achieve fast convergence using an iterative least square method in the ST inversion, having more accurate initial input values and the covariance matrix is crucial. In this study, we used the empirical formula of ST covariance [27] to obtain the statistical value of ST covariance based on band-pass filtering gridding ST.
For two-dimensional process, the covariance function is [27]
( C b b ) l 1 , l 2 = 1 N 1 l 1 1 N 2 l 2 n 1 = N 1 / 2 N 1 / 2 l 1 1 n 2 = N 2 / 2 N 2 / 2 l 2 1 ( b n 1 , n 2 · b n 1 + l 1 , n 2 + l 2 ) ; l 1 = 0 , , N 1 1 , l 2 = 0 , , N 2 1
where l 1 and l 2 are distances in latitude and longitude; N 1 and N 2 are total length in latitude and longitude. According to Equation (36), we can obtain covariance statistics for different distances. The Hirvonen function parameters can then be derived according to the results of covariance statistics using the least square method.
The band-pass filtering gridding ST is the dataset used for computing covariance statistics. To fully obtain the effective statistical information of the band-pass filtering gridding ST, 30 grids were left in the longitude and latitude directions, respectively, such that the maximum distance of the statistical is 90′. The final covariance statistical results are shown in Figure 7, where the blue box indicates statistical results greater than zero, while the black box indicates negative covariance statistical results. Covariance is a non-negative number, and a negative covariance is an invalid statistical value. The results show a covariance statistical value of 0.6365 km 2 ( C 0 ) when the ST distance is zero, and correlation length (CL) of 10.5′ when the covariance value is equal to C 0 / 2 . Using the C 0 and CL as input parameters, the fitting results by the Hirnoven covariance function [21,22] are shown in the red curve in Figure 7. When the statistical covariance results are compared with the fitting results, the comparative analysis shows the values are in good agreement. So it is more appropriate to use the Hirnoven covariance function to obtain the ST covariance matrix in the study area.

3.3. Seafloor Topography Model Inversion

Taking a priori parameters ( C 0 , CL) and ST distance matrix as input data (Figure 8a), the inversion band ST covariance initial matrix is constructed by the Hirvonen covariance function and is shown in Figure 8b.
A reasonable integration radius needs to be set to calculate the VGG/GA using the terrain prism integration method. Based on the recommendations of previous studies [20,21], we selected 30′ as the integration radius. The actual input range of the sea surface VGG/GA data in the inversion process was ( 33.5 °   S 34.5 °   S , 101.5 °   E 102.5 °   E ). Due to the influence of the edge effect, the effective range of the final ST model was reduced by 30′ compared to the initial input range, i.e., 1 ° × 1 ° ( 33.5 °   S 34.5 °   S , 101.5 °   E 102.5 °   E ). To implement the inversion method smoothly, the variance of the input data (i.e., SM, GA, and VGG) must be obtained. For the accuracy evaluation of the SM data, using the third-class survey accuracy in the S-44 (fourth) hydrographic survey standards, when the average depth of the study area is 4702.02 m, the SM accuracy (mean square deviation) is 108.15 m. Given the satellite altimetry GA accuracy of about 2–5 mGal [24,35], the GA accuracy was set to 3 mGal. The accuracy of the VGG derived from the altimetry datasets was not reported. Thus, based on previous research, the accuracy of VGG was set to 15 Eotvos.
We applied the least square method to invert the ST model using the satellite altimetry gravity as input data. We continuously tuned and optimized the ST covariance and initial ST values to iteratively calculate the ST. The accuracy of the ST model was evaluated against the prepared external check data. As shown in Figure 9a, the RMS value of the ST model retrieved by GA changed with the number of inversion iterations. Using the same inversion strategy, the iterative inversion results derived by VGG and integrating the VGG and GA are shown in Figure 9b,c. As shown in Figure 9, the RMS of the inversion ST model tends to converge with the increase in inversion iterations. After seven iterations to restore the ST model, the ST covariance value was 0.6069 km 2 using VGG/GA singularly or combining VGG and GA.
As shown in Figure 9, the result of ST inversion by the nonlinear iterative least square method became stable after two iterations. This suggests that when accurate initial values are used, two iterations are enough to inverse the ST model, consistent with the findings by Calmant et al. [22] and Ramillien and Wright [23]. The accuracy of the ST model using VGG was also found to be much higher than the GA inversion (comparing Figure 9a,b). However, when both VGG and GA are used as input data for ST inversion, the final checking accuracy of the ST model is less than the result of only using either VGG or GA. This result could be affected by the inversion method’s mechanism, which would require further analysis.
The ST model constructed by GA and VGG is named BAT_GA_ILS (Figure 10a) and BAT_VGG_ILS (Figure 10b), respectively, and the ST model recovered by the combined GA and VGG is named BAT_GA_VGG_ILS (Figure 10c). The red box in Figure 10 shows the effective area of ST inversion. Since there is no gravity data input in the area outside the red box, we believe that directly using this part as ST inversion results would not be suitable. The ST results in this area are closely related to the initial ST. Figure 10 shows a noticeable NW-SE trench in the study area and distinct topographic uplifts on both sides of the trench in the northwest direction, similar to the sea surface gravity data.
Compared with the ST inversion methods in the frequency domain, the nonlinear iterative least-square inversion method has the following characteristics: (1) The sea surface gravity data does not need grid processing, and the discrete data also meets the inversion requirement; (2) Multiple gravity data, such as GH, GA, and VGG. can be directly used in modeling; and (3) The influence of ST’s nonlinear term is considered. For (1) and (2), the least square collocation (LSC) can also be performed, but calculating the cross-covariance between seafloor topography and gravity data would be difficult. Unfortunately, the nonlinear iterative least-square ST inversion method involves large-scale matrix calculation in the space domain. The increase in the number of modeling inputs, modeling range, and integral terrain radius would significantly expand the corresponding matrix dimension, increase computational complexity, and reduce efficiency.
In this study, we tested the processing efficiency of the nonlinear iterative least-square inversion method using a personal computer with Intel (R) Core (TM) i7-7700hq CPU @ 2.81GHz, 8 GB RAM, 64-bit operating system. It took nearly seven hours to complete seven ST inversion iterations using GA/VGG alone and about 21.45 h to achieve seven inversion iterations using integrated GA and VGG. Based on these results, in terms of computational efficiency, regression analysis would be preferable in deriving a global ST model [12,13].

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 2 ° × 2 ° ( 33 °   S 35 °   S , 101 °   E 103 °   E ) area, there were 1889 points, while in the final statistical 1 ° × 1 ° ( 33.5 °   S 34.5 °   S , 101.5 °   E 102.5 °   E ) 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 km 2 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.

Author Contributions

D.F. and S.L. designed the research; D.F. and X.L. performed the research; D.F., J.Y., and X.W. wrote the paper. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by the Nature Science Foundation of China (41574020, 41774021, 41674026, 41774018, 41504018, 41674082) and State Key Laboratory of Geo-Information Engineering, No. SKLGIE2019-M-1-3.

Acknowledgments

The authors are thankful to Yi Lin for editing the English text of a draft of this manuscript. Special thanks to the professional English editing service from EditX. The authors are grateful to the reviewers whose valuable comments significantly improved the presentation of the material.

Conflicts of Interest

The authors declare no conflict of interest

References

  1. Smith, W.H.F. Introduction to this special issue on bathymetry from space. Oceanography 2004, 17, 6–7. [Google Scholar] [CrossRef] [Green Version]
  2. Schaeffer, P.; Faugére, Y.; Legeais, J.F.; Ollivier, A.; Guinle, T.; Picot, N. The CNES_CLS11 Global Mean Sea Surface Computed from 16 Years of Satellite Altimeter Data. Mar. Geodesy 2012, 35, 3–19. [Google Scholar] [CrossRef]
  3. Andersen, O.B.; Knudsen, P. The DTU17 global marine gravity field: First validation results. In International Association of Geodesy Symposia; Springer: Berlin/Heidelberg, Germany, 2019. [Google Scholar]
  4. Bostrom, R. Subsurface exploration via satellite: Structure visible in Seasat images of North Sea, Atlantic Continental Margin, and Australia. Int. J. Rock Mech. Min. Sci. Géoméch. Abstr. 1989, 27, A88. [Google Scholar] [CrossRef]
  5. Mayer, L.A.; Jakobsson, M.; Allen, G.; Dorschel, B.; Falconer, R.; Ferrini, V.; Lamarche, G.; Snaith, H.M.; Weatherall, P. The Nippon Foundation—GEBCO Seabed 2030 Project: The Quest to See the World’s Oceans Completely Mapped by 2030. Geoscience 2018, 8, 63. [Google Scholar] [CrossRef] [Green Version]
  6. Tozer, B.; Sandwell, D.T.; Smith, W.; Olson, C.; Beale, J.R.; Wessel, P. Global Bathymetry and Topography at 15 Arc Sec: SRTM15+. Earth Space Sci. 2019, 6, 1847–1864. [Google Scholar] [CrossRef]
  7. Parker, R.L. The Rapid Calculation of Potential Anomalies. Geophys. J. Int. 1973, 31, 447–455. [Google Scholar] [CrossRef] [Green Version]
  8. Hu, M.; Li, J.; Li, H.; Xin, L. Bathymetry predicted from vertical gravity gradient anomalies and ship soundings. Geodesy Geodyn. 2014, 5, 41–46. [Google Scholar] [CrossRef] [Green Version]
  9. Hu, M.; Li, J.; Li, H.; Shen, C.; Xing, L. A program for bathymetry prediction from vertical gravity gradient anomalies and ship soundings. Arab. J. Geosci. 2015, 8, 4509–4515. [Google Scholar] [CrossRef]
  10. Kim, K.B.; Hsiao, Y.-S.; Kim, J.W.; Lee, B.Y.; Kwon, Y.K.; Kim, C.H. Bathymetry enhancement by altimetry-derived gravity anomalies in the East Sea (Sea of Japan). Mar. Geophys. Res. 2010, 31, 285–298. [Google Scholar] [CrossRef]
  11. Hwang, C. A Bathymetric Model for the South China Sea from Satellite Altimetry and Depth Data. Mar. Geodesy 1999, 22, 37–51. [Google Scholar] [CrossRef]
  12. Smith, W.H.F.; Sandwell, D.T. Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry. J. Geophys. Res. Solid Earth 1994, 99, 21803–21824. [Google Scholar] [CrossRef]
  13. Smith, W.H.F.; Sandwell, D.T. Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings. Science 1997, 277, 1956–1962. [Google Scholar] [CrossRef] [Green Version]
  14. Dixon, T.H.; Naraghi, M.; McNutt, M.K.; Smith, S.M. Bathymetric prediction from SEASAT altimeter data. J. Geophys. Res. Space Phys. 1983, 88, 1563–1571. [Google Scholar] [CrossRef]
  15. Fan, D.; Li, S.; Meng, S.; Lin, Y.; Xing, Z.; Zhang, C.; Yang, J.; Wan, X.; Qu, Z. Applying Iterative Method to Solving High-Order Terms of Seafloor Topography. Mar. Geodesy 2020, 43, 63–85. [Google Scholar] [CrossRef]
  16. Jung, W.-Y.; Vogt, P.R. Predicting bathymetry from Geosat-ERM and shipborne profiles in the South Atlantic Ocean. Tectonophysics 1992, 210, 235–253. [Google Scholar] [CrossRef]
  17. Yang, J. Seafloor Topography Estimation from Gravity Gradients. Ph.D. Thesis, The Ohio State University, Columbus, OH, USA, 2017. [Google Scholar]
  18. Yang, J.; Jekeli, C.; Liu, L. Seafloor Topography Estimation from Gravity Gradients Using Simulated Annealing. J. Geophys. Res. Solid Earth 2018, 123, 6958–6975. [Google Scholar] [CrossRef]
  19. Baudry, N.; Calmant, S. 3-D modelling of seamount topography from satellite altimetry. Geophys. Res. Lett. 1991, 18, 1143–1146. [Google Scholar] [CrossRef]
  20. Baudry, N.; Calmant, S. Seafloor mapping from high-density satellite altimetry. Mar. Geophys. Res. 1996, 18, 135–146. [Google Scholar] [CrossRef]
  21. Calmant, S. Seamount topography by least-squares inversion of altimetric geoid heights and shipborne profiles of bathymetry and/or gravity anomalies. Geophys. J. Int. 1994, 119, 428–452. [Google Scholar] [CrossRef] [Green Version]
  22. Calmant, S.; Bergé-Nguyen, M.; Cazenave, A. Global seafloor topography from a least-squares inversion of altimetry-based high-resolution mean sea surface and shipboard soundings. Geophys. J. Int. 2002, 151, 795–808. [Google Scholar] [CrossRef] [Green Version]
  23. Ramillien, G.; Wright, I.C. Predicted seafloor topography of the New Zealand region: A nonlinear least squares inversion of satellite altimetry data. J. Geophys. Res. Space Phys. 2000, 105, 16577–16590. [Google Scholar] [CrossRef]
  24. Sandwell, D.T.; Müller, R.D.; Smith, W.H.F.; Garcia, E.; Francis, R. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure. Science 2014, 346, 65–67. [Google Scholar] [CrossRef] [PubMed]
  25. Wan, X.; Ran, J.; Jin, S. Sensitivity analysis of gravity anomalies and vertical gravity gradient data for bathymetry inversion. Mar. Geophys. Res. 2018, 40, 87–96. [Google Scholar] [CrossRef]
  26. Marks, K.M.; Smith, W.H.F. Radially symmetric coherence between satellite gravity and multibeam bathymetry grids. Mar. Geophys. Res. 2012, 33, 223–227. [Google Scholar] [CrossRef]
  27. Jekeli, C. Spectral Methods in Geodesy and Geophysics; CRC Press: Boca Raton, FL, USA, 2017. [Google Scholar]
  28. Tscherning, C.C. Advanced Physical Geodesy. Eos 1982, 63, 514. [Google Scholar] [CrossRef]
  29. Vanicek, P.; Janak, J.; Huang, J. Mean Vertical Gradient of Gravity; Springer: Berlin/Heidelberg, Germany, 2001. [Google Scholar] [CrossRef]
  30. Santos, M.C.; Vaníček, P.; Featherstone, W.E.; Kingdon, R.; Ellmann, A.; Martin, B.-A.; Kuhn, M.; Tenzer, R. The relation between rigorous and Helmert’s definitions of orthometric heights. J. Geod. 2006, 80, 691–704. [Google Scholar] [CrossRef]
  31. Fukushima, T. Taylor series expansion of prismatic gravitational field. Geophys. J. Int. 2020, 220, 610–660. [Google Scholar] [CrossRef]
  32. Sandwell, D.T.; Smith, W.H.F. Global marine gravity from retracked Geosat and ERS-1 altimetry: Ridge segmentation versus spreading rate. J. Geophys. Res. Space Phys. 2009, 114. [Google Scholar] [CrossRef] [Green Version]
  33. Huang, M.; Zhai, G.; Ouyang, Y.; Lu, X.; Liu, C.; Wang, R. Integrated Data Processing for Multi-Satellite Missions and Recovery of Marine Gravity Field. Terr. Atmos. Ocean. Sci. 2007, 32, 988–993. [Google Scholar] [CrossRef] [Green Version]
  34. Wang, J.; Guo, J.; Liu, X.; Shen, Y.; Kong, Q. Local Oceanic Vertical Deflection Determination with Gravity Data Along a Profile. Mar. Geodesy 2018, 41, 24–43. [Google Scholar] [CrossRef]
  35. Zhu, C.; Guo, J.; Hwang, C.; Gao, J.; Yuan, J.; Liu, X. How HY-2A/GM altimeter performs in marine gravity derivation: Assessment in the South China Sea. Geophys. J. Int. 2019, 219, 1056–1064. [Google Scholar] [CrossRef]
Figure 1. Relationship between the mass body and gravity information.
Figure 1. Relationship between the mass body and gravity information.
Remotesensing 13 00064 g001
Figure 2. Sketch map of regional terrain distribution.
Figure 2. Sketch map of regional terrain distribution.
Remotesensing 13 00064 g002
Figure 3. Flow chart of the seafloor topography inversion process by the nonlinear iterative least-squares method.
Figure 3. Flow chart of the seafloor topography inversion process by the nonlinear iterative least-squares method.
Remotesensing 13 00064 g003
Figure 4. Data used in the study area: (a) gravity anomaly (GA); (b) vertical gravity gradient (VGG); (c) distribution of Shipborne measurements (SM).
Figure 4. Data used in the study area: (a) gravity anomaly (GA); (b) vertical gravity gradient (VGG); (c) distribution of Shipborne measurements (SM).
Remotesensing 13 00064 g004
Figure 5. Coherence between the VGG/GA and seafloor topography (ST).
Figure 5. Coherence between the VGG/GA and seafloor topography (ST).
Remotesensing 13 00064 g005
Figure 6. Gravity, gravity gradient in the inversion band: (a) GA in the inversion band; (b) VGG in the inversion band.
Figure 6. Gravity, gravity gradient in the inversion band: (a) GA in the inversion band; (b) VGG in the inversion band.
Remotesensing 13 00064 g006
Figure 7. Covariance statistical results.
Figure 7. Covariance statistical results.
Remotesensing 13 00064 g007
Figure 8. Matrix: (a) Distance matrix; (b) Covariance matrix.
Figure 8. Matrix: (a) Distance matrix; (b) Covariance matrix.
Remotesensing 13 00064 g008
Figure 9. RMS of the bathymetric model (BAT) derived from gravity data varies with iterations: (a) GA; (b) VGG; (c) GA + VGG.
Figure 9. RMS of the bathymetric model (BAT) derived from gravity data varies with iterations: (a) GA; (b) VGG; (c) GA + VGG.
Remotesensing 13 00064 g009
Figure 10. ST model based on different gravity data: (a) The ST model constructed by GA is named BAT_GA_ILS; (b) The ST model constructed by VGG is named BAT_VGG_ILS; (c) The ST model recovered by the combined GA and VGG is named BAT_GA_VGG_ILS.
Figure 10. ST model based on different gravity data: (a) The ST model constructed by GA is named BAT_GA_ILS; (b) The ST model constructed by VGG is named BAT_VGG_ILS; (c) The ST model recovered by the combined GA and VGG is named BAT_GA_VGG_ILS.
Remotesensing 13 00064 g010
Figure 11. ST model: (a) DTU18; (b) SIO V20.1.
Figure 11. ST model: (a) DTU18; (b) SIO V20.1.
Remotesensing 13 00064 g011
Figure 12. The ratio of checkpoints in different ranges of difference.
Figure 12. The ratio of checkpoints in different ranges of difference.
Remotesensing 13 00064 g012
Figure 13. Relationship between relative error and sea depth. (a) BAT_GA_ILS; (b) BAT_VGG_ILS; (c) BAT_GA_VGG_ILS; (d) DTU18; (e) SIO V20.1.
Figure 13. Relationship between relative error and sea depth. (a) BAT_GA_ILS; (b) BAT_VGG_ILS; (c) BAT_GA_VGG_ILS; (d) DTU18; (e) SIO V20.1.
Remotesensing 13 00064 g013
Figure 14. Spatial distribution of the relative error. (a) BAT_GA_ILS; (b) BAT_VGG_ILS; (c) BAT_GA_VGG_ILS; (d) DTU18; (e) SIO V20.1.
Figure 14. Spatial distribution of the relative error. (a) BAT_GA_ILS; (b) BAT_VGG_ILS; (c) BAT_GA_VGG_ILS; (d) DTU18; (e) SIO V20.1.
Remotesensing 13 00064 g014
Figure 15. Checking difference histogram: (a) BAT_GA_ILS- SIO V20.1; (b) BAT_VGG_ILS- SIO V20.1; (c) BAT_GA_VGG_ILS- SIO V20.1.
Figure 15. Checking difference histogram: (a) BAT_GA_ILS- SIO V20.1; (b) BAT_VGG_ILS- SIO V20.1; (c) BAT_GA_VGG_ILS- SIO V20.1.
Remotesensing 13 00064 g015
Table 1. Statistical results of gravity data and shipborne data in the study area.
Table 1. Statistical results of gravity data and shipborne data in the study area.
Data TypeMaxMinMeanSD
GA (mGal)103.50−167.80−27.9843.30
VGG (Eotvos)231.41−132.22−1.4638.51
SM (m)−6971.99−1339.06−4702.02783.46
Table 2. Statistical results of gravity data at inversion band.
Table 2. Statistical results of gravity data at inversion band.
Data TypeMaxMinMeanSD
GA (mGal)128.90−119.064.64 × 10−938.86
VGG (Eotvos)170.59−114.584.33 × 10−934.42
Table 3. Statistical results of ST models (unit: m).
Table 3. Statistical results of ST models (unit: m).
ST ModelMaxMinMeanMedianSD
BAT_GA_ILS−7087.54−1181.34−4750.35−4731.711002.40
BAT_VGG_ILS−6932.85−1508.22−4853.85−4807.57905.33
BAT_GA_VGG_ILS−7092.84−1173.36−4750.37−4733.731005.73
DTU18−6729.75−1480.00−4539.33−4510.63978.95
SIO V20.1−6831.39−1402.62−4814.20−4792.52883.54
Table 4. Statistical checking results (unit: m).
Table 4. Statistical checking results (unit: m).
ST ModelMaxMinMeanMedianSDRMSCCRA (%)
BAT_GA_ILS447.43−589.43−53.92−25.20169.47177.700.98633.72
BAT_VGG_ILS330.31−246.3844.7247.7086.5797.370.99492.04
BAT_GA_VGG_ILS442.27−590.89−53.92−24.20172.41180.510.98603.78
DTU18749.07−1233.87−271.24−203.50349.86442.450.92399.27
SIO V20.1175.01−182.81−1.05−1.3551.1351.100.99821.07
Note: CC denotes the correlation coefficient.
Table 5. Relative error statistics (%).
Table 5. Relative error statistics (%).
ST Model|Max||Min|MeanSD
BAT_GA_ILS19.820.011.474.14
BAT_VGG_ILS15.610.01−1.002.24
BAT_GA_VGG_ILS19.970.011.484.21
DTU1839.860.005.848.31
SIO V20.16.840.00−0.011.23
Table 6. ST model difference statistics (Unit: m).
Table 6. ST model difference statistics (Unit: m).
Model Comparison|Max||Min|MeanMedianSDRMS
BAT_GA_ILS- SIO V20.1911.270.0163.8528.31198.48208.47
BAT_VGG_ILS- SIO V20.1689.380.03−39.65−47.93120.34126.69
BAT_GA_VGG_ILS- SIO V20.1905.970.1463.8329.53200.85210.72
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fan, D.; Li, S.; Li, X.; Yang, J.; Wan, X. Seafloor Topography Estimation from Gravity Anomaly and Vertical Gravity Gradient Using Nonlinear Iterative Least Square Method. Remote Sens. 2021, 13, 64. https://doi.org/10.3390/rs13010064

AMA Style

Fan D, Li S, Li X, Yang J, Wan X. Seafloor Topography Estimation from Gravity Anomaly and Vertical Gravity Gradient Using Nonlinear Iterative Least Square Method. Remote Sensing. 2021; 13(1):64. https://doi.org/10.3390/rs13010064

Chicago/Turabian Style

Fan, Diao, Shanshan Li, Xinxing Li, Junjun Yang, and Xiaoyun Wan. 2021. "Seafloor Topography Estimation from Gravity Anomaly and Vertical Gravity Gradient Using Nonlinear Iterative Least Square Method" Remote Sensing 13, no. 1: 64. https://doi.org/10.3390/rs13010064

APA Style

Fan, D., Li, S., Li, X., Yang, J., & Wan, X. (2021). Seafloor Topography Estimation from Gravity Anomaly and Vertical Gravity Gradient Using Nonlinear Iterative Least Square Method. Remote Sensing, 13(1), 64. https://doi.org/10.3390/rs13010064

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop