High-precision 1′ × 1′ bathymetric model of Philippine Sea inversed from marine gravity anomalies

The Philippine Sea, located at the edge of the northwestern Pacific Ocean, possesses complex seabed topography. Developing a high-precision bathymetric model for this region is of paramount importance, as it provides fundamental geoinformation essential for Earth observation and marine scientific research, including plate motion, ocean circulation, and hydrological characteristics. The gravity–geologic method (GGM), based on marine gravity anomalies, serves as an effective bathymetric prediction technique. To further strengthen the prediction accuracy of conventional GGM, we introduce the improved GGM (IGGM). The IGGM considers the effects of regional seafloor topography by employing weighted averaging to more accurately estimate the short-wavelength gravity component, along with refining the subsequent modeling of long-wavelength gravity component. In this paper, we focus on seafloor topography modeling in the Philippine Sea based on the IGGM, combining shipborne bathymetric data with the Scripps Institution of Oceanography (SIO) V32.1 gravity anomaly. To reduce computational complexity, the optimal parameter values required for IGGM are first calculated before the overall regional calculation, and then, based on the terrain characteristics and distribution of sounding data, we selected four representative local sea areas as the research objects to construct the corresponding bathymetric models using GGM and IGGM. The analysis indicates that the precision of the IGGM models in four regions is improved to varying degrees, and the optimal calculation radius is 2^′. Based on the above finding, a high-precision ${\mathrm{1}}^{\prime }\times {\mathrm{1}}^{\prime }$ bathymetric model of the Philippine Sea (5–35° N, 120–150° E), known as the BAT_PS model, is constructed using IGGM. Results demonstrate that the BAT_PS model exhibits a higher overall precision compared to the General Bathymetric Chart of the Oceans (GEBCO), topo_25.1, and DTU18 models at single-beam shipborne bathymetric points.

1 Introduction

The Philippine Sea serves as a convergence zone for continental and oceanic plates, resulting in frequent and intense plate tectonic activity (Richter and Ali, 2015; Lallemand, 2016; Holt et al., 2018). The Philippine Sea has complex topography, including trenches, island arcs, ridges, seamounts, basins, and rifts, exhibiting typical characteristics of a trench–arc–basin system. Additionally, it stands among the highly dynamic areas for geospatial exploration globally, attracting significant attention as a focal area for international scientific research. As fundamental geoinformation for marine scientific research, the accurate acquisition and application of high-precision seafloor topography are crucial (Kunze and Smith, 2004; Jena et al., 2012; Hirt and Rexer, 2015; Hu et al., 2020). They provide necessary information to guarantee an in-depth comprehension of the marine environment and support marine development and governance (Ryabinin et al., 2019; Wolfl et al., 2019).

The Seabed 2030 project, a collaboration under the General Bathymetric Chart of the Oceans (GEBCO) and the Nippon Foundation, strives to achieve the creation of a comprehensive global mapping of the seafloor topography by 2030 (Mayer et al., 2018). In May 2023, the International Hydrographic Organization (IHO) announced that Seabed 2030 project had collected mapping data covering 24.9 % of the seabed. Additionally, multiple GEBCO models (such as GEBCO_2019, 2020, 2022, etc.) have been released, which integrate single-beam shipborne bathymetric data, high-resolution multi-beam shipborne bathymetric data, and predicted bathymetry based on satellite altimetry data. For unexplored areas that have not been covered by bathymetric data, the marine gravity field derived from satellite altimetry data (Zhu et al., 2022; Zhou et al., 2023; Wei et al., 2023) is currently the main source for constructing high-precision and high-resolution seafloor bathymetric models, compared to alternative observation methods (Bondur and Grebenyuk, 2001; Smith, 2004; Hilldale and Raff, 2008; Yeu et al., 2018; Tozer et al., 2019; Yu et al., 2022, 2023; Xu et al., 2023).

The prediction of seafloor topography based on satellite altimetry data has undergone a progression from one-dimensional to two-dimensional approaches, and various methods have been developed, such as the gravity–geology method (GGM), the frequency domain method (Parker, 1973; Watts, 1978; Smith and Sandwell, 1994), and the simulated annealing method (Yang et al., 2018, 2020). The GGM, proposed by Ibrahim and Hinze (1972), was initially used to measure the height of bedrock beneath sediments in land areas and subsequently employed in the inversion studies of seabed topography. Based on gravity anomalies derived from satellite altimetry missions, Smith and Sandwell (1994) constructed the corresponding bathymetric model and observed a significant correlation between seafloor topography and gravity anomaly within the wavelength band of 15 to 160 km. Yang et al. (2018) predicted seafloor topography in the western Pacific Ocean using vertical gravity gradient data through simulated annealing. Compared with the shipborne measured depth, the root mean square of the prediction result was 236 m, representing a 22 % improvement in accuracy over SIO (Scripps Institution of Oceanography) bathymetric model. However, the simulated annealing in the spatial domain includes forward and inverse modeling, requiring significant computational resources for modifying the initial model through continuous iteration. While the GGM and S&S (Smith and Sandwell, 1994) method use the linear relationship between seafloor topography and gravity data to construct empirical functions, the accuracy in regions with complex topography needs improvement. The frequency domain method is based on the spectral relation between seafloor topography and marine gravity field, using the first-order term of the Parker formula to approximate bathymetry. This method omits the effect of higher-order terms and includes several complex geophysical parameters.

As one of the commonly used methods for bathymetric prediction, GGM can use shipborne bathymetry and gravity anomalies to build a bathymetric model with high accuracy (Nagarajan, 1994; Wei et al., 2021; Kim et al., 2010; Hsiao et al., 2011; Annan and Wan, 2020). The GGM divides the gravity anomaly into a short-wavelength component and long-wavelength component. The short-wavelength gravity can be calculated using the Bouguer plate formula, while the estimation of the long-wavelength gravity is particularly important. Annan and Wan (2020) adopted an adaptive mesh form to approximate the long-wavelength gravity with a prediction accuracy of 180.20 m in the Gulf of Guinea, comparable to the accuracy of the ETOPO1 model and SIO model. In the conventional GGM, the process of calculating short-wavelength gravity from sea depth according to the Bouguer plate formula is a one-to-one mapping, which ignores the gravity effect caused by surrounding seafloor topography. To overcome this limitation, An et al. (2022) proposed the improved GGM (IGGM), which considers the effect of regional seafloor topography. The IGGM finely calculates the short-wavelength gravity by weighted averaging and refines the subsequent long-wavelength gravity modeling, significantly enhancing the overall accuracy of the bathymetric model.

Due to its prominent strategic position and complex marine environment, constructing an accurate bathymetric model of the Philippine Sea is of utmost importance. Therefore, this paper focuses on constructing a higher-precision bathymetric model for this region, utilizing the IGGM in combination with shipborne bathymetric data and the SIO V32.1 gravity anomaly. To effectively enhance the calculation efficiency of the bathymetric model and be guided by the topographic characteristics and the distribution of shipborne bathymetric data, four representative areas are chosen as the research objects before the overall solution. This preliminary stage involves exploring the optimal values of required parameters and assessing the applicability of IGGM for the Philippine Sea. Subsequently, a ${\mathrm{1}}^{\prime }\times {\mathrm{1}}^{\prime }$ bathymetric model of the Philippine Sea (5–35° N, 120–150° E) is obtained using a priori values of certain parameters. Finally, the accuracy of the bathymetric model is evaluated by comparing with existing models, including the GEBCO_2022, topo_25.1, and DTU18 models.

2 Data sources

The Philippine Sea, located between the two island chains in the western Pacific Ocean, consists of the Philippine Basin, Parece Vela Basin, Shikoku Basin, Kyushu–Palau Ridge, Izu–Ogasawara Trench, Mariana Trench, and Mariana Trough. Recognizing the importance of the seafloor topography as essential geoinformation within this region, this paper focuses on the Philippine Sea (5–35° N, 120–150° E) as the study area, which presents significant challenges for shipborne bathymetry and the inversion of seafloor topography based on satellite altimetry data due to its complex terrain and geomorphology.

This paper used single-beam shipborne bathymetry data obtained from the National Centers for Environmental Information (NCEI). However, due to the large time span over which the bathymetric data were acquired, some data collected before GPS became operational were often poorly localized and contained significant measurement errors (Smith, 1993). Therefore, it is necessary to use a bathymetric model with higher accuracy as a reference (e.g., GEBCO_2022 model) to ensure the quality of the shipboard bathymetric data. This involves interpolating the GEBCO_2022 model to the shipborne bathymetry points, calculating the difference with the shipborne depth, and then excluding points where the difference exceeds 3 times the standard deviation (SD) of the whole difference. The gravity anomalies were obtained from the SIO V32.1 model with a resolution of 1 arcmin, which was released in August 2022. The V32.1 model in the study area is shown in Fig. 1.

Figure 1 SIO V32.1 gravity anomaly model in the study area.

The GEBCO_2022 Grid, released in June 2022, is a continuous global terrain model with a grid spacing of 15 arcsec (GEBCO, 2022), briefly referred to as the GEBCO model. Over time, GEBCO has developed a range of bathymetric datasets and products. GEBCO provides a Type Identifier (TID) grid that specifies the type of source data used for each grid cell, including single-beam shipborne bathymetry, multi-beam shipborne bathymetry, lidar bathymetry, prediction depths derived from satellite gravity, and grid interpolation depths. By continually collecting and incorporating the latest and most relevant bathymetric data, the GEBCO series of models represents the cutting-edge level of accuracy for terrain models in the world.

Furthermore, the DTU18 and topo_25.1 bathymetric models were used in this paper, both of which have a grid resolution of 1 arcmin. The DTU18 model was released by the Technical University of Denmark (DTU) in 2019. SIO has been continuously updating its bathymetric models for an extended period, and the topo_25.1 model, released in January 2023, is the latest version of the current series.

3 Methodology

3.1 Principle of the improved gravity–geologic method

In fact, there is a nonlinear relationship between seafloor topography and marine gravity anomalies. In many geodetic calculations, the nonlinear issue can be linearized by employing a suitable reference field (Hwang, 1999). The principle is to decompose the gravity anomaly into two main components: the regional gravity field representing the long-wavelength gravity, and the residual gravity field representing the short-wavelength gravity. Among them, the long-wavelength component is generated by deeper mass variations and the short-wavelength component is derived from variations in the local bedrock under sediment (Kim et al., 2010; Hsiao et al., 2016). The gravity anomaly (Δg) is composed of the long-wavelength component (Δg_reg) and the short-wavelength component (Δg_res):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }g=\mathrm{\Delta }{g}_{\mathrm{reg}}+\mathrm{\Delta }{g}_{\mathrm{res}}.\end{array}$$

Using the known depth at control points to obtain the short-wavelength gravity component ($\mathrm{\Delta }{g}_{\mathrm{res}}^{j_n})$:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{g}_{\mathrm{res}}^{j_n}=\mathrm{2}\mathit{\pi }G\mathrm{\Delta }\mathit{\rho }(E_{j_n}-D),\end{array}$$

where G is the gravitational constant ($\mathrm{6.672}\times {\mathrm{10}}^{-\mathrm{11}}$ m³ kg s⁻²); Δρ is the density contrast (kg m⁻³) between seawater and bedrock; $E_{j_n}$ is the depth at j_n point; and D is a reference datum, which is generally the deepest depth of control points.

The GGM uses the Bouguer slab Eq. (2) to compute the short-wavelength component of the sea surface point, which represents a single-point calculation form as shown in the left part of Fig. 2 and is not considered rigorous. In response to this problem, IGGM refines the long-wavelength gravity model by considering the short-wavelength gravity effect of regional seafloor topography, as shown in the right-hand side of Fig. 2. Using the control point (j_n) as the center, R is the calculated radius for estimating the effect of seafloor topography, m is the number of sounding points within the calculation range, and $j_n^m$ is the encompassing shipborne sounding point.

Figure 2 Schematic geometries of GGM and IGGM.

Based on Eq. (2), introducing cos ^kθ_m as the weight parameter, the short-wavelength gravity effect ($\mathrm{\Delta }{g}_{\mathrm{res}}^{j_n^m})$ of the surrounding shipborne points on the control point is

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{g}_{\mathrm{res}}^{j_n^m}=\mathrm{2}\mathit{\pi }G\mathrm{\Delta }\mathit{\rho }\left({\cos }^k{\mathit{\theta }}_m\right)(E_{j_n^m}-D),\end{array}$$

where θ is calculated by the arctangent value, based on the depth of the control point and the horizontal distance between the control point and the surrounding sea surface points. k is an unknown parameter determined by the iterative algorithm and is used in conjunction with cosθ_m as a weighting parameter to quantify short-wavelength gravity effects caused by the surrounding points.

Then, the short-wavelength component ($\mathrm{\Delta }{g}_{\mathrm{res}}^{j_n}$) is corrected using weighted averaging:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathrm{\Delta }{g}_{\mathrm{res}}^{j_n}& =\mathrm{2}\mathit{\pi }G\mathrm{\Delta }\mathit{\rho }{\displaystyle \frac{\begin{array}{c}\left({\cos }^k{\mathit{\theta }}_{\mathrm{1}}\right)(E_{j_n^{\mathrm{1}}}-D)+\left({\cos }^k{\mathit{\theta }}_{\mathrm{2}}\right)\\ +(E_{j_n^{\mathrm{2}}}-D)\mathrm{\cdots }+\left({\cos }^k{\mathit{\theta }}_m\right)(E_{j_n^m}-D)\end{array}}{{\cos }^k{\mathit{\theta }}_{\mathrm{1}}+{\cos }^k{\mathit{\theta }}_{\mathrm{2}}+\mathrm{\cdots }+{\cos }^k{\mathit{\theta }}_m}}\\ & =\mathrm{2}\mathit{\pi }G\mathrm{\Delta }\mathit{\rho }{\displaystyle \frac{\sum _{s=\mathrm{1}}^m\left({\cos }^k{\mathit{\theta }}_s\right)(E_{j_n^s}-D)}{\sum _{s=\mathrm{1}}^m\left({\cos }^k{\mathit{\theta }}_s\right)}}.\end{array}\end{array}$$

Subtracting the refined short-wavelength component calculated by Eq. (4) from the gravity anomaly to obtain the long-wavelength component ($\mathrm{\Delta }{g}_{\mathrm{reg}}^{j_n}$) at the control point j_n:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{g}_{\mathrm{reg}}^{j_n}=\mathrm{\Delta }{g}^{j_n}-\mathrm{\Delta }{g}_{\mathrm{res}}^{j_n}.\end{array}$$

The long-wavelength component at control points is gridded using a tension spline function to obtain the long-wavelength gravity field. The long-wavelength component ($\mathrm{\Delta }{g}_{\mathrm{reg}}^i$) at any point i can be calculated through cubic spline interpolation. Subsequently, the short-wavelength component ($\mathrm{\Delta }{g}_{\mathrm{res}}^i$) at the prediction point i is

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{g}_{\mathrm{res}}^i=\mathrm{\Delta }{g}^i-\mathrm{\Delta }{g}_{\mathrm{reg}}^i.\end{array}$$

Finally, based on a variation in the Bouguer formula, the predicted depth (E_i) at point i is inversely calculated using the short-wavelength component:

$$\begin{array}{}\text{(7)}& E_i={\displaystyle \frac{\mathrm{\Delta }{g}_{\mathrm{res}}^i}{\mathrm{2}\mathit{\pi }G\mathrm{\Delta }\mathit{\rho }}}+D.\end{array}$$

3.2 Calculation process of improved gravity–geologic method

The key of IGGM is how to determine and calculate the optimal values of Δρ, k and calculation radius R, so that gravity anomalies can better characterize the basic geoinformation of seafloor topography. In this study, referring to the process of determining Δρ in GGM, the parameters Δρ, R and k are also calculated by iterative method in IGGM. The value ranges of Δρ, R and k are set as 0.1×10³ kg m⁻³–5×10³ kg m⁻³, 0^′–30^′, and 0.1–15, and their corresponding iteration steps are 0.1×10³ kg m⁻³, 0.1^′, and 0.1, respectively. The correlation coefficient and SD between the predicted depth and the shipborne bathymetry at points i are analyzed under the conditions of different parameter values. When the correlation coefficient is the largest and the SD error is the smallest between predicted depths and measured depths (i.e., ${max}_{\mathrm{1}\le i\le n}C{C}_i$ and ${min}_{\mathrm{1}\le i\le n}{\mathrm{SD}}_i$), the corresponding parameters are the optimal values. In the iterative process of IGGM, the parameters are independent of each other. The optimal values of Δρ, R, and k are determined after all the iterations of the parameters have been completed. Figure 3 presents the flowchart of IGGM, providing a detailed visual representation of the process.

Figure 3 Flowchart of the IGGM.

4 Results and discussion

4.1 Parameter determination and applicability verification of IGGM

Due to the large size of the Philippine Sea and the substantial amount of bathymetric data, improving the calculation efficiency of the bathymetric model becomes essential. If certain unknown parameters can be determined before the overall model calculation, it can effectively achieve the purpose. In view of this, according to the different topographic features and bathymetric data distribution, four representative local areas within the Philippine Sea were first selected as the research subregions for analyzing and determining the unknown parameters in the IGGM. Then the optimal parameter values chosen afterwards were used to construct the final bathymetric model of the Philippine Sea. The four selected areas were as follows: A (21–27° N, 130–137° E); B (24–28° N, 138–145° E); C (8–15° N, 128–138° E); and D (10–16° N, 140–148° E), as shown in Fig. 4a. Area A is situated in the Daito Basin and the northern part of the Philippine Basin. Area B includes parts of the Shikoku Basin, the Izu–Ogasawara Islands arc and Trench, along with other adjacent sea areas. Area C is mainly situated in the Philippine Basin and the Parece Vela Basin, with less topographic relief, except for the Palau Ridge at the junction of the two basins. Area D is similar to area B and has complex terrain, including the Parece Vela Basin, and the Mariana Trench with depths exceeding 10 000 m. Table 1 provides the removal statistics of the four areas, and Fig. 4b illustrates the distribution of shipborne bathymetric data after gross error removing.

Table 1 Removal statistics of shipborne bathymetric data, according to the 3σ criterion, with the GEBCO model serving as a reference.

Figure 4 Distribution of four local areas and shipborne bathymetric data (red points are control points, black points are checkpoints, and the basemap is the GEBCO model).

To verify the applicability of IGGM in the Philippine Sea, corresponding bathymetric models were constructed using GGM and IGGM in four areas. Multiple independent cruises were selected within each area, which were not involved in the model calculations, for independent verification purposes (represented by the black points in Fig. 4b) to ensure the reliability of the accuracy evaluation between GGM and IGGM. In both GGM and IGGM, the control points were divided into two parts. The first part was utilized for model calculation under different parameter values, while the second part was used for iteratively selecting the local optimal solutions of required parameters. Once the optimal values of each parameter were determined through iteration, the final bathymetric model was calculated using all the control points. This two-step process allowed for a comprehensive evaluation of the parameter values and ensured that the final bathymetric model was calculated by the best possible combination of shipborne bathymetric data, gravity anomalies, and optimized parameter values. For the division of control points in the iterative method, the conventional GGM still adopted the original proportional selection method, where one control point for iteration was selected at an interval of three points on each cruise, as shown in Fig. 5. On the other hand, IGGM divided control points using the cruise selection method, based on the distribution of shipborne cruises, as shown in Fig. 6. The proportional selection method tended to consider the entire area for obtaining the optimal values of the parameter. Contrastingly, the cruise-distribution-based selection weakened the effects caused by the interpolation of neighboring points.

Figure 5 Proportional selection in GGM. The green control points were used to calculate bathymetric models corresponding to different unknown parameter values, and the white control points were employed for iteratively selecting the optimal value of Δρ. The ratio of green dots to white dots was 3:1.

Figure 6 Selection based on the distribution of shipborne cruises in IGGM. The green control points are used to calculate bathymetric models corresponding to different unknown parameter values, and the white control points were employed for iteratively selecting the optimal parameter values of Δρ, R, and k.

This section mainly involves the following steps. First, the bathymetric models corresponding to the initial values of each parameter were calculated using green control points and were interpolated to obtain the predicted depth at white control points. Second, the correlation coefficient and SD of predicted depths and shipborne measured depths at white control points were calculated. The parameter values that yield the largest correlation coefficient and smallest SD (i.e., ${max}_{\mathrm{1}\le i\le n}C{C}_i$ and ${min}_{\mathrm{1}\le i\le n}{\mathrm{SD}}_i$) were considered optimal. Subsequently, the final bathymetric model was constructed by combining all control points and using the optimal parameter values. The predicted depths at checkpoints were then obtained through cubic spline interpolation, and the prediction accuracies for both GGM and IGGM were evaluated based on the shipborne bathymetric data. The optimal parameter values for GGM and IGGM in the four regions were presented in Table 2. It reveals that the optimal computational radius R of IGGM is 2^′ in different areas, but there are differences in the optimal values of other Δρ and k. Additionally, the density contrast value in area D shows a significant difference. Figure 7 displays the total number of surrounding points within a 2^′ radius centered on each shipborne point, with black points indicating no other single-beam bathymetry point distribution. Statistical analysis shows that the numbers of the black points in the four areas are 0, 454, 131, and 130, respectively. The percentages of surrounding shipborne points exceeding 20 are 72.1 %, 59.1 %, 72.9 %, and 92.8 % for the four areas, respectively. The optimal selection of density contrast in both GGM and IGGM enables the gravity anomalies to better characterize the distribution of seafloor topography (Nagarajan, 1994; Hu et al., 2012; Kim and Yun, 2018). However, it should be noted that the density contrast is used as an empirical parameter only for obtaining the optimal accuracy of the bathymetry model, thus diminishing its original physical significance.

Table 2 Optimal values of parameters for GGM and IGGM.

Figure 7 Total number of surrounding points within 2^′ of each shipborne point as the center (black points represent no other single-beam shipborne points within a 2^′ radius).

The corresponding bathymetric models constructed by GGM and IGGM were interpolated to obtain the predicted depths at the checkpoints. Table 3 displays the accuracy comparison between the predicted depth and the shipborne bathymetry in each area. Results indicate that the IGGM models in the four areas show varying degrees of accuracy improvement. The most significant enhancement is observed in areas A and D, where the SD is reduced by approximately 16.36 % (25.88 m) and 11.05 % (33.39 m). In areas B and C, the improvements of IGGM are limited to within 10 m, which is only slightly better than GGM models. The shipborne bathymetric data are evenly distributed in areas A and B, but area A is located in the Daito Basin with relatively slow topographic changes and exhibits the most substantial accuracy improvement. In contrast, there are Shikoku Basin and Izu–Ogasawara Island arc in area B, whose eastern side is located at the junction of the Philippine Sea plate and the Pacific Plate, and the topographic drop can reach up to thousands of meters. Terrain fluctuation is the main factor that affects the inversion accuracy, so the complex terrain in area B leads to limited accuracy improvement of IGGM. Area C, positioned in the western Philippine Basin and the Parece Vela Basin, has relatively gentle terrain. However, the scarcity of sounding bathymetric data in this region, results in a large error when predicting checkpoint depths using gridding and interpolation processing. Although area D includes the Mariana Trench with drastic topographic variations, it has a large number of shipborne bathymetric data.

The mean absolute percentage error (MAPE) is regarded as an indicator of relative accuracy and defined as the average value of the ratio (positive) of the prediction error to the measured depth:

$$\begin{array}{}\text{(8)}& \mathrm{MAPE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left|{\displaystyle \frac{E_i-H_i}{H_i}}\right|\times \mathrm{100}\mathit{\%},\end{array}$$

where H_i is the shipborne measured depth, and E_i represents the predicted depth. A smaller MAPE indicates a higher accuracy of the bathymetric model. The result shows that IGGM models in the four regions have higher relative accuracy and can better characterize the geoinformation of seafloor topography. In conclusion, compared to GGM models, IGGM models effectively improve accuracy, with the degree of improvement affected by the distribution of shipborne bathymetry and terrain fluctuation. The significant improvement is observed in areas with flat terrain or sufficient and uniform distribution of bathymetric data. Conversely, less data and larger terrain fluctuation result in reduced accuracy improvement. Additionally, gridding and interpolation also introduce varying degrees of errors.

Table 3 Statistics of GGM and IGGM models at checkpoints (units in meters).

Table 4 Optimal values of Δρ and k in each subregion during the construction of the BAT_PS model.

In areas A, B, and C, the GGM model and IGGM models show relatively consistent bathymetry performance overall, as shown in Fig. 8. However, notable differences are observed in area D, particularly concentrated in the vicinity of the Mariana Trench and Mariana Trough, such as 10–12° N, 143–145° E; 12–13° N, 140–141° E; 10–11° N, 140–141° E; and 15–16° N, 143–144° E. The GEBCO model, which incorporates the latest shipborne bathymetric data with high accuracy, was used as the standard to evaluate the accuracy levels of GGM and IGGM models. Figure 9a shows the distribution of seafloor topography in the GEBCO model. Figure 9b and c illustrate the differences between this model, the GGM model, and IGGM model, respectively. Since the GEBCO model has a grid interval of 15 arcsec, it needs to be interpolated to the corresponding grid points using the cubic spline. It is evident that the IGGM model exhibits smaller errors than the GGM model and is closer to the GEBCO model within four ranges marked by black boxes in Fig. 9a. Especially, in the vicinity of the Mariana Trench marked in the black box of Fig. 9a, the GGM model shows significant errors, further indicating the higher-precision performance of the IGGM model in area D as depicted in Fig. 8. Within the red box range of Fig. 9c, the IGGM model exhibits larger errors than the GGM model. This problem may be attributed to the lack of shipborne bathymetric data in this range, resulting in the optimal parameters selected by IGGM being less applicable to a certain local area. According to statistics, differences between the GEBCO model and the GGM and IGGM models within the range of 0–300 m account for 86.2 % and 88.8 %, while the proportions of the differences above 500 m are 6.1 % and 3.5 %, respectively. The GEBCO model, despite having a grid interval of 15 arcsec, exhibits lower true spatial resolution in marine regions that lack of measured bathymetric data, as illustrated in the basemap of area C in Fig. 6.

Figure 8 Comparison of the GGM model and IGGM model in four test seas.

Figure 9 The GEBCO model and the absolute value of its differences in comparison with the GGM and IGGM models.

Based on the above statistics, it is evident that the IGGM, which considers the effect of regional seafloor topography, demonstrates higher accuracy at checkpoints and achieves a significant improvement in the Philippine Sea. These results also confirm that the IGGM has better applicability. For all parameters required in IGGM, the optimal calculation radius in the four areas remains consistent at 2^′, which is not affected by the distribution of seafloor topography and shipborne bathymetric data, but the optimal values of Δρ and the index k in the weight factor differ in each area. Given these findings, a predetermined calculation radius of 2^′ is adopted for the final bathymetry modeling of the Philippine Sea (5–35° N, 120–150° E).

4.2 Bathymetric model based on IGGM in the Philippine Sea

Considering the huge amount of shipborne bathymetric data in the entire Philippine Sea and the limitations in computer capacity, the region (5–35° N, 120–150° E) is divided into nine subregions of 10°×10° each, as shown in Fig. 10a. To ensure data quality, the shipborne bathymetric data within each subregion are preprocessed by removing gross errors that exceed 3σ, with the GEBCO model used as a reference for comparison. Following the procedure of IGGM, the bathymetric models of each subregion are inverted. Finally, the ${\mathrm{1}}^{\prime }\times {\mathrm{1}}^{\prime }$ bathymetric model of the Philippine Sea (BAT_PS) is finally obtained, as depicted in Fig. 10b. During the construction of the BAT_PS model, the optimal parameters values required in IGGM are given in Table 4.

Figure 10 (a) Subregions of the Philippine Sea and distribution of shipborne bathymetric data. (b) BAT_PS bathymetric model constructed by IGGM in the Philippine Sea.

The accuracy of the BAT_PS model was evaluated using three reference models: GEBCO, topo_25.1, and DTU18. The BAT_PS, GEBCO, topo_25.1, and DTU18 models were interpolated to all shipborne bathymetry points within the study area. As all the single-beam bathymetric data were involved in the model calculation, they could serve as an evaluation of internal consistency accuracy. The statistical results were presented in Table 5, with two categories of assessment: primary check (PC) and second check (SC). The primary check is the analysis of the difference between the predicted depth and the measured depth at all shipborne measured points. In the second check, the points with poor prediction quality based on the 3σ criterion are removed before conducting the statistical analysis. The results demonstrate that the BAT_PS model exhibits better agreement with the shipborne bathymetric data in the primary check, with SD and MAPE values of 89.56 m and 9.79 %, outperforming the other three models. MAE represents the absolute magnitude of prediction errors, while MAPE reflects the percentage magnitude of prediction errors relative to the shipborne bathymetry. In the second check after removing gross errors, the accuracy levels of the BAT_PS model and the GEBCO_2022 model are approximately equal, surpassing topo_25.1 model and DTU18 model. Overall, the accuracy ranking from high to low is BAT_PS model, GEBCO model, topo_25.1 model, and DTU18 model.

Table 5 Statistics of BAT_PS, GEBCO, topo_25.1, and DTU18 models at all shipborne bathymetric points (units in meters).

In the statistics of two check accuracies, both before and after gross error elimination, the DTU18 model outperforms the topo_25.1 model, with MAPE values of 20.88 % and 20.26 %, compared to the topo_25.1 model's values of 29.97 % and 27.19 %. In the second checks, the topo_25.1 model has excluded poor-quality points, with both the SD and RMS showing better performance compared to the DTU18 model. However, it is noteworthy that the MAPE of topo_25.1 is still worse. The topo_25.1 model may exhibit a large relative error in shallow seas. Figure 11 shows the relationship between the relative errors between four models and different depths. The relative error is described as the percentage of the prediction error and measured depth. MAPE can be regarded as a comprehensive statistical measure of relative error. As the depth increases, the relative errors within each model exhibit a decreasing trend.

Figure 11 Relationship between relative errors and sea depths.

Table 6 presents the statistics of MAPE at different depths after the second check. Specifically, the topo_25.1 model has a large relative error within a depth of 2000 m, which explains the phenomenon discussed above. After reaching a depth greater than 2000 m, the relative error in the topo_25.1 model gradually becomes smaller than that of the DTU18 model. The construction of the BAT_PS model only relies on single-beam data without the inclusion of multi-beam bathymetric data. Consequently, the accuracy in some deep-sea areas is slightly lower. Additionally, the relative errors in each model increases when the depth exceeds −8000 m. This can be attributed to most of these areas being near trenches, where both the accuracy of shipborne bathymetry and the precision of inversion methods face challenges due to the abrupt terrain.

Table 6 MAPEs for four models at different depths after the second check.

Figure 12 (a) Depth correction grid constrained by shipborne bathymetry. (b) CBAT_PS model constrained by shipborne bathymetry.

To further improve the accuracy of the BAT_PS model and scale down its prediction error, we use shipboard bathymetry data for the final constraints in the following steps (SRTM, 2019): (1) interpolate the BAT_PS model to obtain predicted depths at ship measurement points and calculate the difference from the actual measured depths; (2) supplement grid points located 5^′ away from the ship measurement points with zero-depth differences; (3) use the GMTs (Generic Mapping Tools) module surface to generate a corrected grid (Fig. 12a) by combining the depth differences at ship measurement points and the zero values at grid points; and (4) restore the corrected grid to the PS model to obtain the CBAT_PS (constrained BAT_PS) model, as shown in Fig. 12b. About 96.29 % of the corrected values were counted to be within the range of −100–100 m, indicating that the overall accuracy of the BAT_PS model is still reliable. Finally, after being constrained by shipborne bathymetry, the CBAT_PS model has also been uploaded.

5 Conclusions

Topography, especially bathymetric topography, has been a significant research focus, with unique applications in the field of geoscience. Although multiple techniques for bathymetric inversion and their corresponding improvement methods are temporarily constructed, a fully rigorous bathymetric inversion theory is yet to be established. Consequently, it is necessary to refine and perfect the original method based on the current theory, serving as the foundation and focus of future bathymetric inversion research. In this paper, the IGGM, which considers the short-wavelength component effect of regional seafloor topography, is used to predict bathymetry in the Philippine Sea. The primary conclusions are summarized as follows:

1. Four local areas in the Philippine Sea were selected as research subjects, and corresponding bathymetric models were calculated using both GGM and IGGM. Overall, the accuracy of IGGM at checkpoints outperformed that of GGM, with improvements of 16.36 %, 3.31 %, 3.80 %, and 11.05 %. This verifies the applicability of IGGM in the Philippine Sea. As a result, the optimal calculation radius of IGGM was set at 2^′ for the subsequent construction of the overall bathymetric model of the Philippine Sea.

2. A high-precision BAT_PS model, with a grid of 1 arcmin, was established for the Philippine Sea (5–35° N, 120–150° E) using IGGM. The SD error between the predicted depths of the BAT_PS model and shipborne-measured depths reached 89.56 m at the single-beam shipborne points, and the accuracy of the BAT_PS model was 59.91 m without poor quality points. This accuracy is essentially equivalent to that of the GEBCO_2022 model and is significantly better than the topo_25.1 and DTU18 models. In addition, we use the shipborne bathymetric data to constrain the BAT_PS model and make the final CBAT_PS model.

3. In this study, only single-beam bathymetric data were used for the construction of the BAT_PS model, while multi-beam bathymetric data are abundant in certain areas, such as the western Philippine Sea basin, the Mariana Trench, and Mariana Trough. Future research will explore integrating multi-beam data into the model construction to further enhance the accuracy of the BAT_PS model. Additionally, the accuracy of seafloor terrain inversion heavily relies on the ocean gravity field. The observation data from the Surface Water and Ocean Topography (SWOT) swath altimetry satellite are hoped to further enhance the resolution and accuracy of the marine gravity field, thereby greatly improving the prediction precision of bathymetry.