Using an Iterative Algorithm to Predict Topography From Vertical Gravity Gradients and Ship Soundings

Abstract

With the help of the analytical expression of the vertical gravity gradient (VGG) generated by a rectangular prism, an analytical algorithm for predicting seafloor topography using the VGG data has been studied. Nevertheless, ship sounding data are an essential constraint in solving the seafloor topography. This paper combines ship sounding data with VGG anomaly to predict the seafloor topography. The main research contents include the following: Using the ship soundings and VGG data in the study area, the observation equations about sea depth are established, and the stability of the equations are studied; furthermore, considering the influence of seafloor topography outside the study area on the observation equations, these effects are divided into boundary effects and far-field effects, and different processing methods are proposed. Finally, the method is tested on the East Pacific Rise, only using VGG anomaly and adding the mean value to fix the boundary region, the RMS error of the results is 108.8 m; SIO's model is added to the boundary region and the seven maximum absolute errors are replaced by ship sounding data, the RMS error of the results can reach 94.2 m and the accuracy improvement is 13.42%.

Introduction

As early as 1953, humans have summited Mount Qomolangma, Earth's highest peak. Nonetheless, 67 years later, China's Struggle submersible dived to the maximum depth of 10,909 m in 2020. This shows the difficulty of studying topography; however, the global sea depth plays an important role in the research of physical oceanography, marine ecology, marine geology, and other related geosciences, so the prediction of topography will be a very valuable topic (Becker et al., 2009; Jakobsson et al., 2008; Sandwell et al., 2006). Therefore, the topography has been studied for decades but mainly using satellite altimetry and gravity anomaly data (Baudry & Calmant, 1991, 1996; Calmant & Baudry, 1996; Hsiao et al., 2016; Hwang, 1999; Jung & Vogt, 1992; Lambeck, 1981; Ramillien & Cazenave, 1997; Wang et al., 2001).

Parker (1972) expanded the perturbing field into the form of a Fourier series, which made the calculation faster. At the same time, based on the forward calculation example of the bedrock magnetic field with a seabed depth of 1–2 km, the feasibility of the method was proven. Parker's method provides a theoretical basis for using gravity anomaly to predict topography. After that, isostatic compensation is applied to the study of topography (Detrick & Watts, 1979; McKenzie & Bowin, 1976). Dixon et al. (1983) combined flexural isostatic compensation theory and the linear response function suggested that altimeter data and topography are highly correlated. Baudry et al. (1987) used SEASAT data to precise the location of unsurveyed seamounts, and the accuracy can reach about 15 km under certain conditions. Furthermore, Smith and Sandwell (1997) found the accuracy of the gravity grid is 3–7 mGal with a resolution limit of 20–25 km, and the global topography with a horizontal resolution of 1–12 km was predicted by combining altimetry and ship sounding data.

In addition, Ibrahim and Hinze (1972) used Gravity-Geologic Method (GGM) to estimate the bedrock topography; subsequently, this method has been widely used in the prediction of topography (Hsiao et al., 2011; Kim et al., 2011; Kim & Yun, 2017; Ouyang et al., 2014; Xiang et al., 2017). Kim et al. (2010) also compared GGM with Sandwell and Smith's method, and the results show that GGM has a better effect in short-wavelength components but the latter is more suitable for medium and long wavelengths. However, based on deep sea sampling analysis, the density variation of seabed bedrock is very complex (Expedition 349 Scientists, 2014), and the two methods are also greatly affected by the fluctuation of topography (Brown, 1998; Watts, 1979).

Up to now, VGG has been used less than gravity in predicting topography. In fact, Hubbert (1948) obtained the theoretical formula of gravity generated by a two-dimensional geological body through the line integration method. Based on Hubbert's work, Talwani et al. (1959) applied it to the inversion of the Mendocino submarine fault zone and found that the fault in the south is 3 km thinner than that in the north. Subsequently, Talwani & Ewing. (1960) used the numerical integration method to calculate the gravity anomaly generated by any three-dimensional geological body and applied it to topographic correction. Nagy et al. (2000) derived the analytical gravitational potential and its derivatives for the prism and pointed out a number of applications where they can be used to study various aspects of local gravity fields. However, this formula has been mainly used in the inversion of density interfaces for decades (Geng et al., 2014; Martinez et al., 2013; Qin & Huang, 2016; Zhdanov et al., 2004), and it is not fully used in predicting topography.

Early, Wang (2000) proved that VGG can better reflect the local characteristics of topography than gravity itself. Zhu (2007) discussed the computational method of VGG generated by topography above a geoid based on the VGG anomaly formula of the rectangular prism. In recent years, Kim and Wessel (2016) adopted a new analytic method to compute VGG anomalies generated by complex topography. Particularly, Yang et al. (2018) used an annealing algorithm to predict topography and concluded that the RMS error is 236 m compared with the measured depth data and that the prediction accuracy for topography can be improved by 22%.

In fact, topography with complex shapes can be approximately regarded as a combination of a series of rectangular prisms after appropriate partitioning, and VGG generated by a rectangular prism with constant density has a strict mathematical expression (Nagy et al., 2000; Talwani, 1973), so VGG generated by any topography can be forward computed. Calmant (1994) earlier established observation equations to simulate the prediction of topography from the gravity anomaly formula generated by prisms, and the formula is nonanalytic and needs to be discretized. Besides, the accuracy and error analysis of topography prediction based on VGG are insufficient.

In this paper, based on the analytic formula of VGG generated by rectangular prism, the analytical observation equations between VGG anomaly and topography are established. Then, the stability and resistance to random errors of the observed equations are verified by iterative calculation to further simulate the weakening of various disturbances in test area prediction. Finally, our method is applied to the calculation of the test area, and various algorithms are tested to find the best prediction result.

Method and Theory VGG Generated by Rectangular Prism

As shown in Figure 1, the length, width, and height of the rectangular prism Ω are 2a, 2b, H-h respectively, where H (H is the maximum depth, also called the reference depth) is the height from the bottom of Ω to sea level and h is the height from the top of Ω to sea level. Assuming P is any point at sea level, we need to calculate the VGG generated by Ω at point P.

Figure 1. Schematic diagram of a synthetic rectangular prism.

The local coordinate system is shown in Figure 1, thus the gravitational potential generated by Ω at point P can be expressed as follows: [Image Omitted. See PDF]where G is the gravitational constant ($$6.672\times {10}^{-11}{\text{Nm}}^{2}{\text{kg}}^{-2}$$), $$\rho $$ is the density of Ω, and $$\left({x}_{1},{y}_{1},{z}_{1}\right)$$ is the integral variable. Based on the previous work (Nagy et al., 2000; Talwani, 1973; Yu, Xu, & Wan, 2021), the VGG generated by Ω at point P can be written as follows: [Image Omitted. See PDF]where $$R=\left\{\left({x}_{1},{y}_{1}\right);-a\le {x}_{1}\le a,-b\le {y}_{1}\le b\right\}$$ is a rectangular region, [Image Omitted. See PDF]and $${I}_{R}(h,x,y)$$ is the same as $${I}_{R}(H,x,y)$$ after replacing H by h.

Thus, the expression of VGG generated by Ω at any point P is obtained.

Establishment of Observation Equations

The above discussion is based on a single prism. In fact, a study area can be approximately discretized into multiple rectangular prisms. Therefore, it is assumed that $$R=\left\{(x,y);-a\le x,y\le a\right\}$$ is a square area on the sea surface (such as Figure 2) and $${{\Omega}}_{1}=\left\{(x,y,z);(x,y)\in R,H-h(x,y)\le z\le H\right\}$$ is the corresponding topography, where $$h(x,y)$$ is the height from the top of Ω₁ to sea level. This section will study how to solve $$h(x,y)$$ from VGG anomalies generated by Ω₁.

Figure 2. The 80 km × 80 km of the simulated topography Ω1, the maximum height of topography is 0.9 km and the minimum is 0.1 km. The blue center box 60 km × 60 km inner region R, the influence of 10 km boundary region D outside the R.

The step length is chosen as t and the interval $$[-a,a]$$ into 2N equal parts by the step length t. In this way, the area R can be gridded, and the typical subgrid area is $${R}_{ij}=\left[{x}_{i},{x}_{i+1}\right]\times \left[{y}_{j},{y}_{j+1}\right]$$, $${x}_{i}=it,{y}_{j}=jt$$, $$i,j=-N,{\cdots},0,{\cdots}N-1$$. If t is small enough, the $$h(x,y)$$ in the subgrid area R_ij can be treated as a constant $${h}_{ij}$$. Using Equations 2 and 3, and the principle of superposition, the VGG generated by Ω₁ on sea level can be expressed as follows: [Image Omitted. See PDF]

Removing the influence of sea water from Equation 4, then the VGG anomaly at P on the study area R generated by Ω₁ can be rewritten as follows: [Image Omitted. See PDF]Where $${\Delta}\rho ={\rho }_{s}-{\rho }_{w}$$, $${\rho }_{s}$$ and $${\rho }_{w}$$ are the density of bedrock and seawater, respectively. If the observed value of VGG anomaly on study area R are measured on a discrete grid, we have the following equation: [Image Omitted. See PDF]where $$\left({x}_{p},{y}_{q}\right)$$ is the observation point at sea level R, $$p,q=-N,{\cdots},0,{\cdots}N$$. Equation 6 is the general observation equation for solving $${h}_{ij}$$. When we consider the boundary region D (shown in Figure 2), the approach of this paper is as follows: We only use the simulated VGG anomalies data (which are generated both from inner region R and boundary region D) on inner region R to solve the topography of inner region R and boundary region D, and then the topography of boundary region D is removed and only the topography of inner region R is taken as the final prediction result. Therefore, Equation 6 can be rewritten as follows: [Image Omitted. See PDF]

Hence, Equation 7 is the equation we establish between topography and VGG anomaly. Since Equation 7 is nonlinear with respect to $${h}_{ij}$$, it can be solved in an iterative way after the initial value $${h}_{ij}^{(0)}$$ is given, and the general iteration method is as follows: [Image Omitted. See PDF]

We use the least square method to solve $${h}_{ij}^{(k+1)}$$ in Equation 8 in the case that $${h}_{ij}^{(k)}$$ are given.

Regularization Algorithm

The Equation 6 established in the previous section is derived and analyzed without consideration of uncertainties. However, when dealing with the test area, the situation below the test area is much more complicated because the observed gravity includes contribution from the masses in the boundary area. At this time, we can use Equation 7 to solve and weaken the influence of the boundary region D.

Then, Equation 7 can be linearized, then linearized equations Ah = b can be obtained, where A is the coefficient matrix, h is all the quantities to be solved, and b is the observed VGG anomaly. Since the observation points $$\left({x}_{p},{y}_{q}\right)$$ are all in R, linearized equations Ah = b may produce singularity. To ensure the stability of Equation 7, the regularization equations are introduced as follows: [Image Omitted. See PDF]Where E is the identity matrix of (2N)² × (2N)² order and $$\alpha > 0$$ is the regularization parameter (the unit is $${10}^{-18}{m}^{-2}{s}^{-4}$$). Because $$\alpha > 0$$, there is always a unique solution for Equation 9 with respect to h (Tautenhahn, 1998).

Consequently, when the VGG anomalies are given in advance, we can predict topography by solving Equation 9. Then, the convergence of the iteration method and resistance to random errors of Equation 9 are analyzed in Section 3, and the treatment method of the boundary effect is also studied.

Simulations and Error Analysis Analysis of Stability of the Observation Equations

Hence, the stability of Equation 9 is evaluated by a simulation example. A simulated topography Ω₁ of 80 km × 80 km is shown in Figure 2, which is composed of prisms with length and width of 2 km, whose heights are from 0.1 to 0.9 km. At the same time, the topography Ω₁ is placed under sea level with the depth H ranging from 1 to 6 km (increasing by 1 km), and the VGG anomaly at sea level in the region R generated by the simulated topography Ω₁ can be obtained from forward computations with Equation 4. At the same time, the depth of the boundary region D is fixed by the simulated true value.

The RMS error of the allowable iteration error be set to 10⁻⁴ m and the initial value of $${h}_{ij}^{(0)}$$ be selected as 10 m. Through calculation, the number of iterations and convergence conditions for different H are shown in Figure 3.

Figure 3. Diagram of iteration number and RMS error for different H.

From Figure 3, the accuracy requirements can be met through 5 to 9 iterations for different H. In other words, H only influences the number of times of iterations and does not determine the final convergence results of Equation 9. This simulation example demonstrates that Equation 9 has good stability.

Simulations Combined With Prior Data

Further, for the simulated 80 km × 80 km of Ω₁, the topography under the central 60 km × 60 km region R was taken as the unsolved study area. It is clear that topography outside the inner region R can contribute to VGG anomaly in the inner region R. As a simulation, this paper only considers the influence of 10 km region D outside the R. Yu et al. (2021) discussed the influence of VGG anomaly generated by topography undulation 10 km away, which basically presents a low-frequency feature. In this paper, the effect of low frequency noise is considered as a constant and brought into Equation 9 to be solved together.

The step size t = 2 km is also selected to grid the simulated topography Ω₁. The contribution of boundary region D to the VGG anomalies on region R is called the boundary effect of R as shown in Figure 4.

Figure 4. The schematic diagram of grid simulated topography Ω1. The R is the final unsolved study area. The 10 km region D outside the inner region R is called as the boundary region. The purple arrow indicates that the first two depth values are averaged to obtain the next depth value.

Then, Equation 9 has a total of (4N-19)² equations, where A is the matrix of (4N-19)² × (2N)², h is the (2N)² order vector of all the quantities to be solved, and b is the VGG anomalies observation matrix of (4N-19)² order vector. Because the regularization method is used, the results in the boundary region D may diverge during the iteration, which can influence the final solution. Therefore, for the boundary region D, this paper uses the mean value to fix, that is, the sea depth of the boundary region D is always obtained by taking the mean value of its first two depths (shown in Figure 4). We have [Image Omitted. See PDF]

For the selection of the optimal regularization parameter $$\alpha $$, the simulated topography is placed on the H of 1–6 km (increasing by 1 km), and the corresponding curve is shown in Figure 5 with the RMS error as the standard.

Figure 5. The RMS errors of the results corresponding to different regularization coefficients α

Figure 5 shows that the selection of regularization parameters is different for different H. Generally, the deeper the H is, the smaller the regularization parameters to be selected. The reason is that the deeper the H is, the smaller the coefficient in matrix A, so the selected regularization parameters are relatively small. From the final results, the RMS error of the simulated results can be reduced to less than 2 m by using the regularization method and selecting appropriate regularization parameters. This shows that the regularization method can effectively deal with the boundary effect caused by the boundary region D.

Analysis of Resistance to Random Errors

The simulation experiments in the previous two sections show that Equation 9 have good convergence when there is no error in VGG anomaly. In this section, we turn to discuss the resistance to random errors of Equation 9. According to estimations for the Geosat altimetric data by Khafifid, 1993 and Bouman et al. (2011), the standard deviation for the computed VGG error can be less than 0.85 Eötvös. Therefore, the random error of VGG anomaly data is set to 1 Eötvös in the following simulation arithmetic. This means that the random errors of normal distribution with mean value of 0 Eötvös and standard deviation of 1 Eötvös will be added into the observed VGG anomaly data.

The mean value is used to fix the boundary region D to solve Equation 9, where the mean value is calculated by Equation 10. In addition, in order to analyze or compare the effect by the fixation of mean value for the boundary region D, the true value of the boundary region D is also fixed in iterative computations. Therefore, two fixing methods for the boundary region D are discussed: One is to use mean value and the other is to use the simulated true value. In Section 4, we can approximate the topography model as the true value to verify the effectiveness.

For different H, the RMS error curves of the simulated results from different fixing methods for boundary region D are shown in Figure 6. It can be known from Figure 6 that the RMS error of the simulated results increases with increasing H and its antirandom error decreases obviously because the signal-to-noise ratio is lower. With the boundary region D fixed by the mean value, when the H is 1 and 6 km, the RMS errors of the simulated results are 17.9 and 116.0 m, respectively, and the difference between is nearly 6.5 times, which is similar to the 6 times change in H.

Figure 6. The RMS error curves of the simulated results from different fixation methods.

When the boundary region D is fixed by the simulated true value, the RMS error of the simulation results is relatively reduced; that is, the accuracy of the boundary region D can effectively improve the final result. Compared with the use of mean value fixation, the accuracy of true value fixation is improved by 6.1–33.4 m with increasing of the H from 1 to 6 km, and the relative improvement rate ranges from 24.5% to 34.1%.

Moreover, we take the eight absolute maximum errors in the simulation results as the fixed points and compute the final solution again. The statistical results of the H of 3 and 4 km are shown in Table 1. Because the results of other sea depths are basically similar, the maximum sea depth of the actual computation of this paper is 3.5 km in Section 4, which is beneficial for comparison and explanation. In Table 1, Max_T, Min_T, and RMS_T represent the maximum error, minimum error, and RMS error of the fixation results of the simulated true value only for the boundary region D, while Max_C, Min_C, and RMS_C represent the maximum error, minimum error, and RMS error after taking the eight absolute maximum errors as the fixed points, and μ represents the accuracy improvement rate of fixing eight points relative to only true value fix: [Image Omitted. See PDF]

Table 1 The Statistical Table of Results With the Eight Absolute Maximum Errors as Fixed Points

H(km)	True value fixed (m)			Max-error fixed (m)			Proportion
	MaxT	MinT	RMST	MaxC	MinC	RMSC	μ
3	181.6	−148.9	33.9	103.2	−102.2	31.2	8.0%
4	140.9	−200.3	48.9	124.4	−119.8	46.8	4.3%

Table 1 shows that after fixing points, the final results have been improved to various degrees. In terms of RMS error, the H values of 3 and 4 km are improved by 8.0% and 4.3%, respectively, corresponding to the addition of only eight absolute maximum error points. This calculation example shows that fixing the boundary region D with the true value and the eight absolute maximum error points can improve the accuracy of simulated results. This algorithm will also be regarded as the basis for the actual computation in Section 4.

Actual Computation Ship Soundings and VGG Anomalies Data

The previous sections have theoretically verified the effectiveness of the proposed algorithm in predicting topography from the VGG data on sea level. In this section, the East Pacific Rise, whose longitude from 252.8° to 253.6° East and latitude from − 6.4° to − 5.8° South, is selected to predict its topography, and the ship sounding data from Figure 7a shows that the measured topography is relatively complex, including mountains, valleys, and ridges with large topographic fluctuations, as well as relatively flat areas.

Figure 7. (a) The ship sounding data from NOAA's National Geophysical Data Center and (b) distribution of the ship sounding data after a 2 km × 2 km grid. The purple box is the inner region and the outer region is the boundary region.

The ship sounding data are chosen from NOAA's National Geophysical Data Center (NGDC, https://www.ncei.noaa.gov/maps/geophysics/) and shown in Figure 7a. There are 6580 ship survey points in the selected test area, and the maximum sea depth, the shallowest sea depth, and the average sea depth are −3425.0, −1912.0, and −3013.3 m, respectively. Similarly, a grid step of 2 km is selected, and the average value of all ship sounding points in each grid is taken as the sea depth of the corresponding grid, and the grid results are shown in Figure 7b. There are 1249 ship sounding grid points after the grid, including 777 ship sounding grid points in the inner region, accounting for 60% (777/1296) of all 1,296 grid points in the test area. In other words, 40% of the test area still has no bathymetric data. One of the significant purposes of topographic prediction is to better fill the 40% of the test area without ship sounding data.

The VGG anomaly model is from the German Research Center for Geosciences (GFZ, https://icgem.gfz-potsdam.de/tom_longtime) and is shown in Figure 8a, and the maximum, the minimum, and the standard deviation are 51.15, −21.0, and 10.43 Eötvös, respectively. Based on the previous simulation calculations, 10 km outside the solved test area is also selected as the boundary region, as shown in Figure 8a, that is, the solved test area is inside the white dotted line (the final prediction result is shown in Figure 8b) and the boundary region is outside the white dotted line.

Figure 8. (a) Distribution of vertical gravity gradient (VGG) anomalies model from the German Research Center for Geosciences, the white box is the inner region and the outer region is the boundary region; (b) the predicted topography with the boundary region fixed by the mean value (The boundary region has been intercepted.); (c) the predicted topography with the boundary region is fixed by the Scripps Institution of Oceanography (SIO's) model; and (d) the predicted topography with the boundary region is fixed by the SIO's model and the seven absolute maximum errors are used as the fixed points. The red dotted line with the latitude of −6 is the profile line we selected from east to west, and its corresponding topography and VGG anomaly will be presented in Figure 10.

Estimation of Algorithm Results

Based on the VGG anomaly model data in Figure 8a, fix the boundary region (shown in Figure 8a) with the mean value of each iterative results, then use Equation 9 for iterative computations, and finally select the inner area as the prediction results shown in Figure 8b. The prediction results are compared with all 777 ship sounding points in the inner area after gridding. The statistical results are shown in the second row of Table 2, and the maximum, minimum, and RMS errors are 498.3, −570.9, and 108.8 m, respectively. The average depth of the predicted results is −2970.2 m, and the prediction validity can reach 3.66% (108.8 m/2970.2 m) of the average depth.

Table 2 The Statistics of the Results of the Three Algorithms (The Units Are Meters)

Methods	Average depth	Max error	Min error	RMS	τ	β
Mean value fixed (A1)	−2970.2	489.3	−570.9	108.8	3.66%	0%
Model fixed (A2)	−2964.4	512.6	−563.6	101.4	3.42%	6.80%
Max error fixed (A3)	−2960.2	414.7	−411.3	94.2	3.18%	13.42%

Furthermore, the boundary region is fixed with model data from 1 ′ × 1′ global topography model of the Scripps Institution of Oceanography (SIO), this global topography model is obtained by combining altimetry data and ship sounding data (Sandwell et al., 2014; Smith & Sandwell, 1997). The same iterative method is used for the solution. The predicted topography is shown in Figure 8c, and the accuracy of the results is shown in the third row in Table 2. We draw the error distribution of the prediction results without fixed sea depths in this paper as shown in Figure 9a, at the same time, make statistics on the absolute error and the corresponding number as shown in Figure 9b. It can be known from Figure 9b that the absolute error is mainly concentrated within 200 m, accounting for 88.54% (688/777); while the absolute error greater than 400 m has 7 points, accounting for 0.90% (7/777), and the 7 points are all concentrated within the purple circle in Figure 9a.

Figure 9. (a) Statistical error distribution of prediction results compared with ship surveys and (b) statistical absolute error and corresponding number.

We choose the 7 absolute maximum errors for the following reasons: The absolute maximum errors generally appear at the places with the largest topographic fluctuations. If this area with the largest fluctuations can be fixed by the ship sounding data, the prediction accuracy can be effectively improved. Then, in the actual prediction, we only need to measure the depth of the sea area similar to seamounts (corresponding to the place with the largest fluctuation of VGG anomaly), which can effectively save the cost of ship surveying.

If the sea depths with the seven absolute maximum errors are used as the fixed points, the predicted topography is shown in Figure 10d. The results of the comparison with ship sounding points are shown in the fourth row of Table 2.

Figure 10. (a) The pink line is the vertical gravity gradient (VGG) anomalies data along the selected section of GFZ model. The cyan line is the VGG anomalies data along the selected profile of the prediction results. The shaded part is the relative error of the two. (b) The cyan star points are the ship sounding data along the selected profile. The pink line is predicted bathymetric depths.

Table 2 shows the result statistics of the three algorithms, where τ represents the ratio of the RMS error to the average depth and β represents the RMS improvement ratio of the A2 and A3 algorithms relative to algorithm A1. According to Table 2, the average depth of the results of the three algorithms is not much different. The boundary region is fixed by the SIO's model in algorithm A2, which compared with the mean value fixed of algorithm A1, the accuracy of algorithm A2 is higher, and the RMS error of the final result is improved by nearly 6.80%. This actual calculation result is consistent with the conclusion of the simulation calculation in Figure 6. This shows that the accuracy of the boundary region has a great impact on the improvement of the accuracy of the algorithm in this paper. The reason is that the VGG anomalies generated by the boundary region have an impact on the inner region and present the characteristics of high frequency. Because the prediction of topography mainly also utilizes the high-frequency characteristics of VGG anomalies, the effects of boundary region will lead to the reduction of the signal-to-noise ratio and affect the prediction results.

The absolute maximum error of the result of algorithm A2 is mainly concentrated in the red circle area of Figure 8c, which is also the test area with large VGG anomalies changes in Figure 8a, so the topographic fluctuation is more intense. Furthermore, if the seven absolute maximum errors in algorithm A2 are fixed, that is, algorithm A3, the maximum and minimum errors of the prediction results are reduced accordingly, and the RMS is reduced by 13.42% compared with algorithm A1. Compared with the algorithm A1, the τ values of the other two algorithms are also relatively reduced, which indicates that the relative error is decreasing, that is, the prediction accuracy is improving.

To analyze the prediction results more clearly, we take the profile line selected in Figure 8d as the section (the longitude and latitude are converted into length units), and the corresponding graph is shown in Figure 10. Figure 10a shows the VGG anomalies curve of the model and prediction results (along the profile line selected in Figure 8d), and the maximum and minimum differences are 4.11 Eötvös and −5.33 Eötvös, respectively. At the same time, the statistical difference shows that the mean value of the difference between the two is −3.23 Eötvös and the standard deviation is only 0.16 Eötvös, indicating that the fluctuation of the difference between the two is very small.

There is a large difference between model predictions and measured VGG at the peak and trough (black box in Figure 10a), which corresponds to the place with large topographic fluctuation in Figure 10b. There are two main reasons for the large error: The first is that the algorithm uses the 2-km grid, that is, the horizontal resolution of the final result is 2 km, and the mean depth will be less representative as each grid increases with the topographic fluctuation; on the other hand, the accuracy of the VGG anomaly model is also limited, and the signal-to-noise ratio is even worse at greater topographic fluctuations (Yu, Xu, & Wan, 2021).

From the comparison between the ship sounding data distribution in Figure 10b and the prediction result curve in this paper, it can be known that the similarity between the prediction result and the ship sounding depth is very high between 0 and 50 km. However, in the black box of Figure 10b, the difference between the two becomes larger. The large difference might be caused by comparison of data sets with different resolutions. Limited by resolution, the VGG is unable to recover short-wavelength topography. In addition, the ship sounding data used in this paper are single beam from decades ago, and its measurement accuracy may have errors (Smith, 1993). However, with the development of multibeam measurement, this problem will be improved to a certain extent, but we are unable to obtain the multibeam data of the test area in this paper.

Discussion and Conclusions

Compared with previous studies that use bilinear interpolation function to weaken low-frequency effects (Yu, Xu, & Wan, 2021), this paper refines the sources of interference errors as “boundary effects” and “far-field effects,” and the far-field effects are regarded as a constant and introduced into the observation equations to solve. Comparison of the VGG anomaly of prediction results and the model data is shown in Figure 10a, the mean value of the difference between the two is −3.23 Eötvös and the standard deviation is only 0.16 Eötvös. In other words, the difference value is closer to a systematic error characteristic, which may be caused by far region topographic fluctuation or flexural isostatic, and we are studying and trying to weaken this impact.

The method of treating it as a constant in this paper may not be optimal but this constant can be introduced directly into the observation equation to solve. This algorithm is simple and improves the prediction efficiency. Of course, the use of filtering to better eliminate the influence caused by far field is also worthy of further study. In addition, this paper sets the boundary area as 10 km, but the selection of this size needs to consider the topography in the far region. If the topographic fluctuations in the far region are severe, the boundary region can be selected to be larger to effectively weaken the influence of high-frequency errors.

In this paper, the step length of the test area grid is 2 km; that is, the topography fluctuation within 2 km × 2 km is only represented by an average sea depth, which means that the fluctuation within 2 km is difficult to distinguish. Sandwell and Smith (1997) recovered marine gravity anomaly from Geosat and ERS-1 with an accuracy of about 4–7 mGal, with the development of altimetry satellites, the accuracy of a new marine gravity anomaly can reach about 2 mGal (Sandwell et al., 2014). Nevertheless, in order to improve the accuracy of seafloor topography prediction, we should not only improve the accuracy of gravity data but also improve the resolution. The Next-Generation SWOT (Surface Water Ocean Topography) Satellites may be a revolutionary improvement in marine gravity recovery (Morrow et al., 2019; Yu, Hwang, et al., 2021), and the SWOT ocean products will also provide the 1 km spatial resolution over oceans (Neeck et al., 2012).

In summary, the main results of this paper are as follows: (a) This paper develops a new analytical algorithm for predicting topography through an iterative solution from VGG anomalies and verifies the effectiveness of this method through simulation; (b) in order to address the boundary effect, the boundary region is introduced, and the regularization is adopted when solving the linearized equations to ensure the stability of the equations; and (c) the theoretical method in this paper is applied to the test area, by adding the fixed algorithm to the boundary region and the seven maximum absolute errors, and the RMS error of the prediction result reaches 94.2 m, which verifies the effectiveness of this method. Even if we do not use any ship sounding data and the topography model, the RMS error of the final result can reach 108.8 m.

This paper provides an analytic iterative algorithm for predicting topography. At the same time, for the boundary effect problem of traditional geophysical inversion, more regions are solved and only the middle part is taken as the final result. It may be a new idea to use the regularization method to deal with the boundary effect in this paper. However, there are still great challenges in the prediction topography, such as gravity signal separation and flexural isostatic problems.

Acknowledgments

The authors thank GFZ, NGDC, and SIO for supplying actual data. I would also like to thank Dr. Tang He for his help in the calculation and Dr. Wang Qiuyu for his advice on the drawing. At the same time, we also have to thank the reviewer Junjun Yang and another anonymous reviewer for their professional comments on this paper. This work is funded by the National Nature Science Funds of China (Grant 41774089).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The ship sounding data are from NOAA's National Geophysical Data Center (NGDC, https://www.ncei.noaa.gov/maps/geophysics/), we can get the required data by entering the latitude and longitude; VGG anomalies model can be obtained from the German Research Center for Geosciences (GFZ, http://icgem.gfz-potsdam.de/tom_longtime), and the model selected is “dV_ell_RET2012_plusGRS80” in this paper; and the topography model can be downloaded from the Scripps Institution of Oceanography (SIO, https://topex.ucsd.edu/cgi-bin/get_data.cgi), we can get the required data by entering the latitude and longitude. The data used in this paper are given in the attachment. Software Instructions: In this paper, Fortran is used for the calculation and Matlab and Origin are used for the drawing, for which we would like to express our gratitude.