The content of this article comes from the Journal of Surveying and Mapping, Issue 2, 2024 (drawing review number: GS Jing (2024) No. 0297)
A Gaussian surface function estimation method for inverting seabed topography using perturbation gravity data
Zhai Zhenhe
, Sun Zhongmiao, Guan Bin, Ma Jian, Li Duan
State Key Laboratory of Geographic Information Engineering, Beijing 100029
Funds: National Natural Science Foundation of China (42174001;41774018) and State Key Laboratory of Geographic Information Engineering (SKLGIE2023-22-5)
Abstract:The 1'×1' grid disturbance gravity data in the South China Sea is obtained by using the SARAL satellite altimetry data from January 2017 to December 2020, and the accuracy reaches 5.5 mGal by comparing with the shipborne gravity data. In this paper, a method of using gravity data and Gaussian surface function to solve regional characteristic parameters and then estimate the seabed topography model is proposed, and the gravity data obtained by the SORAL satellite inversion is used to carry out computational experiments in the South China Sea. The results show that the accuracy of the 1'×1' seabed terrain estimated under the division condition of a group of 10×10 grids is improved by about 10 m compared with the prior model by using the ETOPO-1 prior model, and the accuracy of the 1'×1' seabed topography estimated under the division condition of a group of 9×9 grids is improved by about 9 m compared with the prior model by using the DTU18 prior model. The calculation results show that the gravity data obtained by satellite altimetry inversion and the five characteristic parameters obtained by Gaussian surface function can represent the surface characteristics of the seabed topography in the corresponding region to a certain extent, and then the prior seabed model can be refined without relying on external measured data. For surface estimation, the smaller the grid division, the more the surface function reflects the regional variation characteristics. Therefore, for the future development of satellite altimetry technology, higher resolution gravity field detection technology is expected to continue to improve the ability to invert seabed terrain details.
Keywords: perturbation gravity, seabed topography, SARAL satellite, Gaussian surface estimation, ETOPO-1 DTU18
翟振和, 孙中苗, 管斌, 等. 利用扰动重力数据反演海底地形的高斯曲面函数估计方法[J]. 测绘学报，2024，53(2)：231-238. DOI: 10.11947/j.AGCS.2024.20230204ZHAI Zhenhe, SUN Zhongmiao, GUAN Bin, et al. An estimation method of seabed topography based on Gauss surface function using ocean gravity data[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(2): 231-238. DOI: 10.11947/j.AGCS.2024.20230204 阅读全文：http://xb.chinasmp.com/article/2024/1001-1595/20240203.htm
INTRODUCTION
At present, seabed topography information is mainly obtained by shipborne multi-beam surveying, airborne laser bathymetry, multispectral inversion, and satellite altimetry inversion [1-6]. Shipborne multi-beam measurement has the highest accuracy and has become the main means of high-precision seabed topography mapping in the sea area, but it is difficult for shipborne measurement technology to achieve global uniform coverage due to the limitation of the survey platform and the objective environmental conditions of the ocean. Ocean satellite altimetry can obtain global gravity data with a resolution of up to 1′ in about 2.5 years [7-9], and then can invert and obtain a certain resolution and evenly distributed seabed topography, which has become one of the important means for the construction of 3D models of seabed topography. The mainland's self-developed new-generation ocean altimetry satellite was launched in March 2023, which is expected to further improve the level of fineness of seabed topography information in the wide sea area. At present, research institutions in various countries widely use data including shipborne bathymetry and satellite almetry to build seabed topography models, such as the National Centers for Environmental Information (NCEI), the Technical University of Denmark (DTU) and the Scripps Institution of Oceanography of Oceanography (SIO) and other research institutions have successively published global seabed topography grid models such as ETOPO-1, DTU18, and SIO V19.1 [10-11]. Based on satellite altimetry data and shipborne bathymetric data, domestic scholars have also obtained local or global seabed topography, such as BAT_WHU2020 [12]. In terms of seabed topography inversion methods, the gravity anomaly data obtained from satellite altimetry are widely used to invert the seabed topography through gravity geology, analysis, frequency domain and least squares configuration [13-25], and the above methods generally use linearization to take the first-order term for calculation, and Ref. [26] takes into account the influence of the higher-order term to construct a nonlinear second-order correlation model, which can be used to invert the seabed topography with a certain accuracy under the condition of knowing a small amount of shipborne water/gravity data. Although great progress has been made in the inversion of seabed topography using satellite altimetry data, it is still necessary to increase research and expand inversion ideas considering the complex functional relationship between gravity field and seabed topography. Based on the stochastic correlation between gravity field and seabed topography, this paper proposes to estimate the new seabed topography based on the prior model by constructing surface features, which provides a reference for further expanding the application value of satellite altimetry data.
1 Gravity inversion model of ocean disturbance considering singular influences
Considering the potential advantages of perturbation gravity in the geodetic boundary value problem, the perturbation gravity delta δgp of the grid point is calculated by using the accurate and fast inversion method of perturbation gravity [27]. (1)
where γ denotes normal gravity, ξ and η denote the meridian and unitary components of the perpendicular deviation, respectively, α denotes the azimuth from the calculated point to the observation point, and the derivative of the kernel function K(ψpq) is denoted as
(2)
When the angular distance of the spherical surface approaches 0 ψ, the singularity problem occurs in Eq. (2), in order to avoid singularity, the calculation formula of perturbation gravity at the singularity point is deduced under the condition of ignoring the influence of the second order and above terms in Ref. [2].
(3)
where ξy and ηx represent the gradient of the meridian component of the perpendicular deviation in the latitude direction and the gradient of the unitary component of the perpendicular deviation in the longitude direction, respectively, and s0 represents the size of the radius of the grid.
2 Estimation method of Gaussian surface function for inversion of seabed topography from sea gravity data
Therefore, this paper proposes the hypothesis that the surface features of the gravity data in a small area are in good agreement with the surface features of the seabed topography in the corresponding region, and if the accurate gravity field feature information is known, it is possible to recover the corresponding seabed topography features, and the seabed topography in the small area can be estimated. In order to quantitatively analyze the surface variation characteristics, the Gaussian surface function is used to describe the variation characteristics. The two-dimensional Gaussian surface function can be used to describe the data distribution of two dimensions, and plays an important role in image filtering and edge detection, and can generate Gaussian convolution kernels, which can be used for Gaussian filtering in image processing, which can effectively remove Gaussian noise, and can also play the role of classification, clustering and regression in machine learning kernel data mining, so as to predict data more accurately. The Gaussian surface function is represented by five characteristic parameters, namely, the regional peak, the longitude and latitude corresponding to the regional peak, and the variation amplitude of the regional data in the longitude and latitude directions.
Suppose that the gravity data of the regional grid obtained by satellite altimetry inversion or other means is δg(i) (the grid size is set to be 1′×1′), i=1,2, ..., n. The data are divided into small regions (m>3) consisting of m×m 1′ ×1′ grids from north to south and west to east. In each small region, the five characteristic parameters GT, GxT, GyT, CGx, and CGy are obtained by using the known gravity data according to the Gaussian surface function model (Eq. (4)) and the principle of least squares
(4)
In the formula, GT represents the peak value of the gravity field (such as gravity anomaly and perturbation gravity) in the local area, GxT and GyT represent the longitude and latitude information corresponding to the peaks, CGx and CGy represent the amplitude of the change of the gravity field in the direction of longitude and latitude, and xi and yi represent the longitude and latitude information corresponding to each gravity field variable.
The prior seabed topography model Hs(i) is introduced, i=1,2, ..., n, and the grid size is 1′×1′. Using the prior seabed topography data, five characteristic parameters HT, HxT, HyT, CHx and CHy are obtained according to the Gaussian surface function model (Eq. (5)) and the principle of least squares
(5)
In the formula, HT represents the peak value of the seabed topography in the local area, HxT and HyT represent the longitude and latitude information corresponding to the peak, CHx and CHy represent the amplitude of the seabed topography change in the longitude and latitude directions, and xi and yi represent the longitude and latitude information corresponding to each seabed topography variable.
According to the analysis of the existing land gravity, elevation, ocean gravity and bathymetry data (including grid data and survey line data), the two types of data are basically the same in the changing peak position, so GxT and GyT can be used to represent the longitude and latitude corresponding to the peak of the seabed topography to be estimated, and CGx and CGy can be assigned a certain scale factor (Equation (6) and Equation (7)), and then used as the variation range of the seabed topography to be estimated in a small area
(7)
The peak HTnew of the seabed topography to be estimated is cycled in the range of (HT-1000, HT+1000) and the scale is cycled in the range of (-1/10, 1/10). Then the solution is carried out according to Eq. (8).
(8)
When the mean square deviation between Hest-i and the prior seabed topography data Hs(i) is less than a certain threshold, the cycle ends and the small-area estimation of seabed topography is output.
Repeating the above steps, all the small areas can be calculated, and the final output is the seabed topography of the complete area, as shown in Figure 1.
| |
| 图 1 海底地形反演计算流程Fig. 1 Flowchart of seabed topography inversion |
| Diagram options |
3 Calculation test
3.1 Gravity inversion of sea disturbances
The test data is based on the GDR-F version of the SARAL satellite released by the French Space Agency, and the time span is from January 2017 to December 2020, during which the SORAL satellite is in a drifting orbit, so the data is relatively dense. The regional range of the inversion data used is 0°N-20°N, 105°E-125°E, and the regional range of the output gravity product is 111.4°N-112.9°N, 10°E-12.4°E×. In order to analyze the sea surface distribution density of the SARAL satellite data, the official 1 Hz sea surface height data was first plotted and displayed in the local sea area calculated by the gravity field, as shown in Figure 2 and Figure 3, respectively.
| |
| 图 2 SARAL卫星1 Hz海面高数据在1′×1′格网内的分布(蓝色圆点代表测量值)Fig. 2 Distribution in 1′×1′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements) |
| Diagram options |
| |
| 图 3 SARAL卫星1 Hz海面高数据在3′×3′格网内的分布(蓝色圆点代表测量值)Fig. 3 Distribution in 3′×3′ grid of 1 Hz sea surface height data from the SARAL satellite (blue dots represent measurements) |
| Diagram options |
As can be seen from Fig. 2 and Fig. 3, the distribution of 1 Hz sea surface height data of SARAL satellite in the 1′×1′ grid is very uneven, and it cannot fully meet the resolution requirements of 1′×1′, and the resolution of the measurement point grid is better than that of 3′×3′ in general. The actual gravity field inversion is calculated according to the size of 1′×1′ and 3′×3′ grids, respectively, and the following three reasons are mainly considered: first, it is to compare with the existing 1-point grid ship-borne gravity data, so as to give the accuracy of the SARAL satellite inversion gravity field; DTU18 model;Thirdly, considering that the proposed method needs to be solved in a small area by using known grid observations, and theoretically the finer the grid division, the better the effect, and considering that the known shipborne water depth grid model in the experimental area is 1′×1′ grid, therefore, the calculation of 1′×1′ grid is conducive to the verification of the method. In order to avoid the influence of the far-zone term, the removal and recovery method is adopted, and the 2160-order EGM2008 model is used in the prior model, and the gravity of the disturbance in the South China Sea is calculated according to Eq. (1)-Eq. (3), as shown in Fig. 4.
| |
| 图 4 SARAL卫星测高数据反演获得的1′格网扰动重力数据Fig. 4 The 1′ grid disturbing gravity data obtained from SARAL satellite altimetry data |
| Diagram options |
Compared with the shipborne gravity data in this area, the accuracy of the 1′×1′ grid disturbance gravity data obtained by using the SARAL satellite altimetry data and the EGM2008 model reaches 5.5 mGal, and the accuracy of the 3′×3′ grid disturbance gravity data reaches 4.8 mGal (Table 1), which is in line with the normal level of the existing satellite altimetry inversion gravity field.
表 1 SARAL卫星测高数据反演扰动重力精度Tab. 1 Estimation accuracy of disturbing gravity inversion from SARAL satellite
3.2 Regional seabed topography inversion
The test area is consistent with the inversion range of gravity data, specifically 10°N-12.4°N, 111.4°E-112.9°E, the average sea depth of the test area is 3828 m, the maximum depth is 4316 m, the minimum depth is 1549 m, and the variation range is 562 m. The gravity data used are the perturbation gravity data of the 1′ ×1′ grid obtained in the previous section. The prior seabed topography model was based on the 1′×1′ grid ETOPO-1 bathymetric model published by the National Environmental Information Center of the United States (see Table 2 for comparison results with shipborne water depth) and the 1′×1′ grid DTU18 bathymetric model published by the Technical University of Denmark (see Table 2 for comparison results with shipborne water depth), and the shipborne bathymetric check grid data was obtained on the basis of the measured data of the South China Sea from 2018 to 2019 (the survey line spacing is about 2 km) of the Ministry of Natural Resources, with a resolution of 1′.
表 2 先验水深模型的精度统计Tab. 2 Estimation accuracy of seabed topography based on prior bathymetric model
Since the surface function needs to estimate 5 parameters, the gravity observation grid data has the minimum amount of calculation, and the accuracy of the parameters estimated by too few observations is low, and too many observations will reduce the adaptability range of the parameters, thereby reducing the recovery of high-frequency information. In order to investigate the influence of different grid size division on the calculation results, the regional grid division was set as 3×3 grids (i.e., each small area has 9 1′×1′ grids), 4×4 grids, 5×5 grids, 6×6 grids, 9×9 grids, 10×10 grids and 12×12 grids, and the calculation results under different grid division conditions are shown in Table 3, in which the three-dimensional graphics of seabed topography inverted based on 10×10 grids are shown in Figure 5.
表 3 基于ETOPO-1先验水深模型的海底地形估计精度Tab. 3 Estimation accuracy of seabed topography based on ETOPO-1 prior bathymetric model
| |
| 图 5 基于ETOPO-1先验水深模型估计得到的海底地形Fig. 5 Seabed topography estimated based on ETOPO-1 prior bathymetric model |
| Diagram options |
The estimated seabed topography obtained based on the above method and the DTU18 prior bathymetric model is shown in Table 4 when compared with the shipborne data, and the 3D graphics of the seabed topography based on the 9×9 mesh inversion are shown in Figure 6.
表 4 基于DTU18先验水深模型的海底地形估计精度Tab. 4 Estimation accuracy of seabed topography based on DTU18 prior bathymetric model
| |
| 图 6 基于DTU18先验水深模型估计得到的海底地形Fig. 6 Seabed topography estimated based on DTU18 prior bathymetric model |
| Diagram options |
From the above results, it can be seen that although the gravity data with more than 5 observations can be estimated by the Gaussian surface in theory, in fact, due to the resolution of the gravity data itself and the nonlinear characteristics of the Gaussian surface function, the solution of the 5 parameters is still greatly affected by the amount of data and the distribution of data. In addition, if there is a complex trend of seabed topography, such as a non-smooth and discontinuous interface, the Gaussian surface function is limited by the resolution of gravity observation data, and it is difficult to reconstruct complex seabed topography features. According to the comprehensive analysis, for the ETOPO-1 prior model, the 10×10 grid division has the best effect, and the estimated seabed terrain accuracy is improved by about 10 m compared with ETOPO-1. For the DTU18 model, the 9×9 grid division achieves the best effect, and the estimated seabed terrain accuracy is improved by about 9 m compared with that of DTU18. On the whole, for the inversion of the seabed topography of the 1′×1′ grid, the partition solution is carried out according to the 9×9 grid division, and the better results can be obtained without relying on the support of external sea depth data.
4 Conclusion
Based on the stochastic correlation characteristics of gravity data and bathymetric data, this paper explores the method of estimating seabed topography by using satellite altimetry inversion gravity data and surface characteristic function, and obtains the following useful conclusions.
(1) The overall grid density of the 1 Hz data measured by the SARAL satellite in the South China Sea is better than 3′×3′, but it cannot meet the requirements of the 1′×1′ grid resolution, and the accuracy of the inverted 3′×3′ grid gravity field is 4.8 mGal, and the inverted 1′×1′ grid gravity field accuracy is better than 5.5 mGal.
(2) ETOPO-1 and DTU18 are both internationally published authoritative bathymetry models, and their model construction itself uses data from various sources such as satellite inversion bathymetry and shipborne bathymetry, and the above model is used as a priori model, and the gravity data obtained by satellite altimetry is used to obtain the characteristic parameters through the Gaussian surface function in the form of 9×9 grid division to estimate the seabed topography, and its accuracy is about 10 compared with the prior bathymetry model The increase of m level shows that the existing models such as ETOPO-1 and DTU18 can still be refined without relying on external measured bathymetric data by using the latest gravity observation data, and also shows that the idea of estimating by the characteristic parameters of Gaussian surface function is feasible.
(3) In general, the proposed method is more suitable for those sea areas where the shipborne survey data cannot be effectively covered or the coverage is uneven, and the high-quality gravity data can be used in these areas and combined with prior information for optimization and improvement. In practice, the estimation accuracy of the Gaussian surface function is actually related to the density and accuracy of the observation data, if the observation data resolution is high, the adaptability of the Gaussian function will be strong, and the surface function can reflect the detailed characteristics of the seabed topography area change (high-frequency information). Therefore, the finer the gravity field, the more obvious the detailed features that can theoretically reflect the topography of the seabed. For the future development of satellite altimetry technology with ocean surveying and mapping as the mission goal, it is necessary to improve the high real physical resolution of the sea surface while improving the accuracy, so that it is possible to use the high-resolution gravity field to improve the inversion ability of the fine seabed topography in the unsurveyed sea area.
About the Author
About the first author:ZHAI Zhenhe (1980—), male, Ph.D., associate researcher, research direction is space geodesy and vertical datum. E-mail: [email protected]
First trial: Zhang Lin review: Song Qifan
Final Judge: Jin Jun
information