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

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 km2 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.


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 1 l dv (1) 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 − 2ur cos(ψ)| u=R (2) 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 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:

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

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

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 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 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: Hence, the result of linearizing Equation (10) is 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: Equation (13) can then be transformed to the matrix form:

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]: Then, Equation (16) can be expressed as Similar to GA calculation, the linearizing result of Equation (18) is where ξ zp is the difference between the nonlinear and linear calculation results of VGG at the research point. a zi can be expressed as

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 ξ N p is the difference between the nonlinear and linear calculation results of GH at the research point. a Ni can be expressed as

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 · · · · · · · · · · · · · · · 1 a τ 1 n a τ 2 n . . . Equation (27) can be written as v = A·δb − (L − A·b 0 − ξ) (28) 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 where P L is the weight matrix of observations. The final result of the terrain adjustment iŝ 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) 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 ss accounts for the uncertainties of the measurement data, while the covariance matrix C bb accounts for uncertainties in the ST. If the measurement is regarded as independent of each other, the measurement error matrix is given by E ss = σ 2 s ·I (33) where I is the unit matrix; and σ 2 s 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 C 0 is a priori uncertainty between topographic heights; and ψ(r, r ) is the spherical distance between the locations r and r'. Ramillien [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. Equation (32) converge rapidly, the terrain covariance of each iteration can be adjusted and optimized according to the covariance propagation law:  (35) 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.

Data Preparation and Pre-Processing
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.

Data Preparation and Pre-Processing
, 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  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.  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.  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.

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

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 bb ) l 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. length (CL) of 10.5′ when the covariance value is equal to 0 / 2 C . Using the 0 C 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.

Seafloor Topography Model Inversion
Taking a priori parameters ( 0 C , 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.

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. 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.

Seafloor Topography Model Inversion
Taking a priori parameters ( 0 C , 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. 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 2 km 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 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 crosscovariance 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 leastsquare 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 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].

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. version 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 crosscovariance 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 leastsquare 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].

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  The ST models generated by the different approaches (Figures 10 and 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 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.  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 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 no- 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.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 20 ticeable 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 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 estimated 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  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.

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 2 km when the distance is zero, and the correlation length is 10.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.