*Article* **Unified Land–Ocean Quasi-Geoid Computation from Heterogeneous Data Sets Based on Radial Basis Functions**

**Yusheng Liu and Lizhi Lou \***

College of Surveying and Geo-Informatics, Tongji University, 1239 Siping Road, Shanghai 200092, China; lyssyl@tongji.edu.cn

**\*** Correspondence: llz@tongji.edu.cn

**Abstract:** The determination of the land geoid and the marine geoid involves different data sets and calculation strategies. It is a hot issue at present to construct the unified land–ocean quasi-geoid by fusing multi-source data in coastal areas, which is of great significance to the construction of land– ocean integration. Classical geoid integral algorithms such as the Stokes theory find it difficult to deal with heterogeneous gravity signals, so scholars have gradually begun using radial basis functions (RBFs) to fuse multi-source data. This article designs a multi-layer RBF network to construct the unified land–ocean quasi-geoid fusing measured terrestrial, shipborne, satellite altimetry and airborne gravity data based on the Remove–Compute–Restore (RCR) technique. EIGEN-6C4 of degree 2190 is used as a reference gravity field. Several core problems in the process of RBF modeling are studied in depth: (1) the behavior of RBFs in the spatial domain; (2) the locations of RBFs; (3) ill-conditioned problems of the design matrix; (4) the effect of terrain masses. The local quasi-geoid with a 1 resolution is calculated, respectively, on the flat east coast and the rugged west coast of the United States. The results show that the accuracy of the quasi-geoid computed by fusing four types of gravity data in the east coast experimental area is 1.9 cm inland and 1.3 cm on coast after internal verification (the standard deviation of the quasi-geoid w.r.t GPS/leveling data). The accuracy of the quasi-geoid calculated in the west coast experimental area is 2.2 cm inland and 2.1 cm on coast. The results indicate that using RBFs to calculate the unified land–ocean quasi-geoid from heterogeneous data sets has important application value.

**Keywords:** fusion of heterogeneous data; unified land–ocean quasi-geoid; radial basis functions; Remove–Compute–Restore; EIGEN-6C4

### **1. Introduction**

The geoid is a fundamental element in determining the shape of the Earth [1]. The calculation of a high-quality geoid requires firstly constructing the Earth's gravity field with high precision and resolution. With the enrichment of measurement means, the breadth and depth of geospatial data are being continuously improved by terrestrial, shipborne, airborne and satellite gravity surveys, etc. [2–7]. The unified land–ocean quasi-geoid can be constructed by fusing heterogeneous data sets in the boundary areas between land and ocean, which is beneficial to advance the construction of land–ocean integration. At present, calculation methods of the gravimetric geoid are mainly divided into analytical methods and statistical methods. Analytical methods include the classical Stokes, Hotine, Molodensky and Helmert integral algorithms [8–12]. The defects of analytical methods lie in the strict requirement of the boundary surface and the difficulty fusing multi-source gravity signals. Statistical methods represented by the least-squares collocation (LSC) have significant advantages in fusing heterogeneous data sets. LSC is first introduced into the study of local gravity field approximation by Krarup and Moritz. Tscherning, Rapp, Hwang and other scholars have performed a lot of later research [13–17]. The defect of LSC is that it is difficult to construct an appropriate and accurate local covariance model.

**Citation:** Liu, Y.; Lou, L. Unified Land–Ocean Quasi-Geoid Computation from Heterogeneous Data Sets Based on Radial Basis Functions. *Remote Sens.* **2022**, *14*, 3015. https://doi.org/10.3390/rs14133015

Academic Editor: Xiaogong Hu

Received: 16 May 2022 Accepted: 21 June 2022 Published: 23 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

However, as a statistical algorithm, LSC is still a good choice. This article mainly focuses on analytical methods.

Radial basis functions have been widely used in local gravity field modeling in recent years due to their simple function forms and the ability to fuse multi-source data. RBFs are first proposed as a function approximation theory in mathematics [18]. The Multiquadric (MQ) kernel function proposed by Hardy and the point mass function are first introduced into Earth's gravity field approximation [19]. Weightman is one of the first scholars to apply the point mass function in physical geodesy. This new method can replace the classical spherical harmonic function to express the Earth's gravity field. Reilly and Herbrechtsmeier apply the point mass model to the fusion of multi-source data earlier. They use simulated altimetry data to invert marine gravity anomalies and combine the measured gravity anomalies on land to construct a unified local gravity field, whose accuracy reaches 20 mGal after being verified [20]. Barthelmes designs a free-positioned point mass optimization algorithm based on the least-squares adjustment by setting four free parameters on each mass point, which effectively reduces the number of mass points and significantly improves the calculation efficiency [21]. Lehmann further studies the free-positioned algorithm, which can be used more flexibly to determine the local geoid with measured gravity data [22]. These studies lay the foundation of radial basis function modeling. Then, a series of high-order RBFs such as radial multipoles and Poisson wavelets are derived from the point mass function [23,24]. In addition, some scholars introduce the classical Blackman kernel, the Shannon kernel and the spherical spline function into the RBF model [25,26]. After years of development, scholars have performed a lot of research on the application of various RBFs in local gravity field modeling. The research contents mainly focus on the selection of RBF types, the treatment of ill-conditioned problems and the determination of the spatial position of RBFs.

Tenzer and Klees carry out experiments in the plain area to compare the performance of the point mass kernel, the Poisson kernel, radial multipoles and Poisson wavelets. According to the results, they conclude that these RBFs can obtain almost the same accuracy of gravity field modeling when the depth of RBFs is chosen properly [27]. Bentel et al. use simulated data to model the local gravity field based on the Shannon low-pass kernel, the Shannon high-pass kernel, the Blackman low-pass kernel, the cubic polynomial kernel, the Poisson multipoles kernel, the Abel–Poisson kernel truncated and the Abel–Poisson kernel. The experimental results show that the Blackman low-pass kernel, the cubic polynomial kernel and the Abel–Poisson kernel truncated perform best [28]. In the process of RBF modeling, the design matrix may be ill conditioned due to the uneven distribution of observations and excessive number of RBFs. Tikhonov regularization is currently the mainstream method to deal with ill-conditioned problems. Wu et al. adopt zero-order and first-order Tikhonov regularization to solve ill-posed equations and prove that firstorder regularization has better performance [29]. Based on the Tikhonov regularization, Liu et al. analyze the defects of the L-curve method and variance component estimation (VCE) in determining regularization parameters and then propose two combined methods, VCE-Lc and Lc-VCE, which are proven to be superior to traditional methods [30,31]. The spatial positions of RBFs have a great influence on the modeling result. Eicker uses various spherical grids to place RBFs such as the Reuter grid, the Triangle-Center grid and the Triangle-Vertex grid, which are more evenly distributed than the geographical grid [32]. Klees and Witter propose an adaptive selection of RBFs before parameter estimation according to the distribution, signal variation and noise of observations [33]. Tenzer and Klees, respectively, use the GCV and RMS minimization methods to determine the optimal depth of RBFs and establish a linear functional relationship between the depth and the correlation length of RBFs. Experiment results prove that the optimal depth is related to the type of RBFs [27]. Tenzer et al. subsequently find in the study of mountainous areas that the RBF model solution is very sensitive to change in depth of even several hundred meters when regional shape fluctuation and gravity signals vary greatly [34]. To sum up the current research by scholars, the main problems are as follows: (1) the selection of RBF

types and the determination of spatial positions have not been standardized, which needs further research; (2) the studies carried out using real gravity data and the experiments carried out in the land–ocean boundary areas are still too few.

Considering the above problems, we propose a multi-layer RBF model to fuse heterogeneous data sets for the unified land–ocean quasi-geoid. In Section 2, two coastal areas and relevant experiment data are detailed. The modeling workflow is given, and the modeling methodologies are introduced in detail, including (1) the characteristics of RBFs in the spatial domain; (2) the spatial position of RBFs, i.e., the resolution and depth of spherical grids; (3) Tikhonov regularization to deal with ill-conditioned problems of the model design matrix; (4) residual terrain model (RTM) to represent the impact of terrain mass. The local quasi-geoid with a 1 resolution is calculated by setting up multi-layer RBF networks, respectively, on the flat east coast and the rugged west coast of the United States in Section 3. In Section 4, we discuss the results and shortcomings of the above research and propose an outlook for future work. Section 5 summarizes the main research content and conclusions of this article.

#### **2. Data and Method**

This article studies the RBF model to fuse multi-source gravity data. Numerical experiments are carried out in the land–ocean junction areas. The experiments are conducted in two topographically different areas, one on the east and the other on the west coast of the United States. The east coast of the United States is relatively flat, expending inland from the coastline into a broad plain. On the west coast, under the influence of the Rocky Mountains, the land elevation begins to rise rapidly not far from the coastline. Details of the relevant experiment data are described in this section. In addition, the general modeling process is summarized, and the core modeling methodologies are studied in depth.

### *2.1. Data Preparation*

The east coast experiment area is located in North Carolina, USA, with an average elevation of approximately 8 m. The target area range is between 34◦ and 37◦ latitude and between −78◦ and −74◦ longitude. The west coast experiment area is in the border area of Oregon and Washington State in the northwest of the United States. Its land part is rugged with an average elevation of approximately 334 m. The target area range is between 44◦ and 47◦ latitude and between −126◦ and −122◦ longitude. The regional topography is shown in Figure 1.

**Figure 1.** Regional terrain in two coastal areas of the United States: (**a**) east coast; (**b**) west coast.

The gravity data used in the experiments include measured terrestrial, shipborne, airborne gravities and gravity anomalies derived from satellite altimetry. Terrestrial and shipborne gravity data come from the NGS99 gravity data set provided by the National Geodetic Survey (NGS). NGS99 is a compilation of the measured terrestrial and shipborne gravity data across the USA (data source: https://www.ngdc.noaa.gov/mgg/gravity/19

99/data/regional/ngs99, accessed on 1 March 2021). The product provides free-air gravity anomalies (FAA), Bouguer anomalies, etc., after a series of preprocessing. In the east coast experiment area, we use 2860 terrestrial gravity points and 1678 shipborne gravity points, which are shown in Figure 2a. The distribution of shipborne gravity points is uneven, which is reflected in the dense distribution of data on the shipping route but many gaps outside the route. The shipborne gravity signals vary greatly, while the terrestrial gravity signals vary more gently and are densely distributed in most areas. In the west coast experiment area, there are 2814 terrestrial gravity points and 2288 shipborne gravity points, as shown in Figure 2b. Different from the east coast experiment area, the distribution of shipborne gravity points in the west coast area is more even.

**Figure 2.** Distribution of terrestrial and shipborne gravity points: (**a**) east coast; (**b**) west coast.

Due to the data gaps between ship lines, satellite altimetry and airborne gravities with more uniform distribution are needed as supplements. Satellite altimetry used in this article are DTU15 gravity anomalies with a 2 resolution (data source: https://ftp.space.dtu. dk/pub/DTU15, accessed on 1 March 2021). Under the influence of the complex terrain on the ground, the accuracy of gravity anomalies inversed from satellite altimetry in inshore areas is not so good, but their high resolution and uniform distribution can make up for the deficiency of shipborne gravity to a certain extent. The distribution of DTU15 on sea is shown in Figure 3.

**Figure 3.** Distribution of DTU15 gravity points: (**a**) east coast; (**b**) west coast.

Airborne gravity data come from the GRAV-D project, which is developed by the NGS to redefine the vertical datum of the USA. GRAV-D currently covers most areas of the USA, providing full-field gravity data, which can be turned into free-air gravity disturbances (FAD) or free-air gravity anomalies (data source: https://www.ngs.noaa.gov/GRAV-D, accessed on 1 March 2021). In the east coast experiment area, 5367 airborne gravity points are selected, and their average flight altitude is approximately 5460 m. In the west coast

experiment area, the number of airborne gravity points is 7458, and the average flight altitude is approximately 6938 m. The distribution of airborne gravity data is shown in Figure 4. Airborne gravity data have high resolution and uniform distribution, so they are not limited by ground terrain conditions. The shortcoming of airborne data is that the high flight altitude leads to low data accuracy, thus it is difficult to simulate the full-band gravity signal on ground points.

**Figure 4.** Distribution of airborne gravity points: (**a**) east coast; (**b**) west coast.

GPS/leveling data are used to verify the accuracy of the unified land–ocean quasigeoid (data source: https://geodesy.noaa.gov/GPSonBM, accessed on 1 March 2021). In the east coast experiment area, due to the flat land terrain, there are many measured GPS/leveling points with a total of 807, which is 177 in the west coast experiment area, as shown in Figure 5. GPS and leveling surveys are impossible to be carried out offshore. So, this article uses GPS/leveling points along the coast to check the accuracy level of the unified land–ocean quasi-geoid on sea.

**Figure 5.** Distribution of GPS/leveling points: (**a**) east coast; (**b**) west coast.

#### *2.2. RBF Modeling Strategies*

In this article, the local quasi-geoid is calculated based on the RBF model following the basic framework of the RCR technique [35]. After removing GGM and terrain effects from gravity observations, the residual gravity is included in the RBF model as input data. Helmert VCE is used to evaluate various observations [36]. Tikhonov regularization is used to solve ill-conditioned problems that may occur in the design matrix. A small number of points are selected from GPS/leveling data as control points, and the remaining data are as internal checkpoints. Under the constraint of control points, the optimal RBF network is determined to construct the model. The quasi-geoid calculated is compared with internal checkpoints, and the standard deviation (STD) of the difference between them is taken as the accuracy indicator. The flowchart is shown in Figure 6.

**Figure 6.** Flowchart of RBF modeling.

When the spherical grid is used to place the RBFs, not all grid points should be included in the model, which is because the RBFs have strong localization characteristics, and their energies are mainly concentrated around the center. So, the observations near RBFs are the main contribution to the simulated gravity signal. If the RBFs without enough observations around are included in the model, the local gravity field will be over-parameterized, and the reliability of the model solution will be reduced. In order to locate and eliminate the redundant RBFs, it is usually necessary to introduce a filter radius *RI*:

$$R\_I = \chi \cdot \upsilon\_{0.5} \tag{1}$$

where *I* is the number of RBFs; *ν*0.5 denotes the correlation length of RBFs, which specifically refers to the straight-line distance between the RBFs and observations when the absolute value of RBFs decays to half of the maximum value. *χ* represents correlation length parameter. With the center of RBFs as the spherical center, if the number of all observations within the radius *RI* is greater than *q*, the RBFs will remain; otherwise, they will be excluded. *χ* and *q* need to be determined based on the actual situation. In this article, we specify *q* = 0. *RI* ranges between 50 and 90 km. We determine the suitable value of *χ* through some trial calculations, ensuring *RI* meets the modeling requirements. The general adaptive screening process is shown in Figure 7.

**Figure 7.** Processing flow of adaptive screening technique.

Considering the altitude of airborne gravity data is much higher than the terrestrial, shipborne and satellite altimetry data, the one-layer RBF network cannot simulate the full gravity signals well. Therefore, we design a multi-layer RBF network to fuse these four types of data, which are shown in Figure 8. We are currently only using two layers, and more layers are of course encouraged to be used.

**Figure 8.** Multi-layer RBF networks.

In Figure 9, surface gravity denotes terrestrial, shipborne and satellite altimetry data; *dA* and *dB* denote the depth of grid A and B. The specific modeling process is as follows:

1. GGM and RTM are removed from the terrestrial, shipborne and satellite altimetry observations to obtain residual gravity anomalies, *ΔgA* = *Δg* − *ΔgGGM* − *ΔgRTM*. The RBF network A is determined by the STD minimization. After the three kinds of gravity data are fused to calculate the RBF model parameters, the airborne gravity points are taken as prediction points and the corresponding model gravity disturbance *δgA* will be calculated.


**Figure 9.** Behavior of RBFs in the spatial domain: (**a**) the IMQ kernel; (**b**) the Poisson kernel; (**c**) radial multipoles of order 1; (**d**) Poisson wavelets of order 1.

#### *2.3. RBF Modeling Methodology*

RBF is a nonlinear function with the local characteristic and radial symmetry. In essence, the RBF model only depends on the relative position relationship between the computed points and the center of RBFs. Represent the spatial position of RBFs as *y*(*y*1, *y*2, *y*3), which is usually placed on a sphere inside the Earth such as the Bjerhammar sphere. The point *x*(*x*1, *x*2, *x*3) outside the sphere is used to represent the spatial position of the observations. RBFs can be expressed as:

$$\Psi(x,y) = \sum\_{m=0}^{\infty} \Psi\_m(2m+1) \left(\frac{R\_B}{|x|}\right)^{m+1} P\_m\left(\pounds^T \mathcal{Y}\right)\_{\prime} |y| \langle R\_{B\prime}|x| \rangle R\_B \tag{2}$$

where *RB* denotes the radius of the Bjerhammar sphere; *Pm* denotes the Legendre polynomial of degree m; *ψ<sup>m</sup>* is the Legendre coefficient, which is the shape factor of RBFs, determining their properties in the spatial and frequency domain; *x*ˆ and *y*ˆ are the unit vector of *x* and *y*. Define the depth of RBFs as: *ds* = *RB* − |*y*|.

According to the Runge–Krarup principle, the external disturbing potential *T* can be approximated by a harmonic function completely embedded inside the Earth. So, the linear combination of RBFs can be used to approximate *T*:

$$T(\mathbf{x}) = \frac{GM}{R\_B} \sum\_{i=1}^{K} \beta\_i \Psi(\mathbf{x}, y\_i) \tag{3}$$

where *K* denotes the number of RBFs; *β<sup>i</sup>* denotes the unknown coefficient of the RBF model. Based on the RCR technique, *T* can be decomposed into three components:

$$T = T\_{GGM} + T\_{Terrian} + T\_{res} \tag{4}$$

where *TGGM* denotes long-wave signals represented by GGM; *TTerrian* denotes high-frequency gravity information implied by terrain masses, used to represent short-wave signals; *Tres* denotes residual disturbing potential. According to the functional relationship between the observations and *T*, the linear combination of RBFs can be further used to simulate the observations:

$$\begin{aligned} \Delta \mathbf{g}\_{res}(\mathbf{x}) &= \sum\_{i=1}^{K} \mu\_{i} \left( -\frac{\partial}{\partial |\mathbf{x}|} \Psi(\mathbf{x}, y\_{i}) - \frac{2}{|\mathbf{x}|} \Psi(\mathbf{x}, y\_{i}) \right) \\ \delta \mathbf{g}\_{res}(\mathbf{x}) &= -\sum\_{i=1}^{K} \mu\_{i} \frac{\partial}{\partial |\mathbf{x}|} \Psi(\mathbf{x}, y\_{i}) \\ \zeta\_{res}(\mathbf{x}) &= \sum\_{i=1}^{K} \mu\_{i} \frac{\Psi(\mathbf{x}, y\_{i})}{\gamma} \end{aligned} \tag{5}$$

where *μ<sup>i</sup>* = *GM RB <sup>β</sup>i*; *<sup>Δ</sup>gres*(*x*) denotes residual gravity anomaly; *<sup>δ</sup>gres*(*x*) denotes residual gravity disturbance; *ζres*(*x*) denotes residual height anomaly. In this article, all geoid undulations have been converted to height anomalies (GPS/leveling observations), ensuring the unity of elevation datum:

$$
\mathbb{Z} \approx \mathcal{N} - \frac{\Delta \mathbb{g}\_B}{\overline{\mathcal{Y}}} H \tag{6}
$$

where *ΔgB* denotes bouguer gravity anomaly; *γ* denotes mean normal gravity; *H* denotes the topographic height. It is worth noting that the above formula is an approximate formula.

Equation (5) can be abstracted as a linear observation model:

$$y\_p = A\_p \mathfrak{x} + \varepsilon\_p \tag{7}$$

where *yp* denotes the observation vector of class *p*; *x* is the vector of model coefficients; *ε <sup>p</sup>* is the error vector. The equation can be solved using the least-squares estimation. Due to the lack of the prior information related to the error of observations, the weight of all kinds of observations can be determined by the VCE technique.

#### 2.3.1. Characteristics of RBFs in the Spatial Domain

Several types of RBFs widely used in local gravity field modeling are as follows: the IMQ kernel, the Poisson kernel, the radial multipoles kernel and the Poisson wavelets kernel, represented by Ψ*IMQ*, Ψ*pk*, Ψ*rm* and Ψ*pw* respectively. We use these RBFs' analytical expressions to construct the model. See the related research for their specific formulas [37]. Convert the analytical formula of RBFs to the unit sphere, and then, respectively, draw their pictures in the spatial domain when *ds* = 300 km, as shown in Figure 9. It can be seen from the figure that RBFs have obvious localization characteristics. In the position far from the center of RBFs, the signals' energies rapidly decay, which is conductive to the concentration of local gravity signals.

#### 2.3.2. RBF Networks

The design of the RBF network is divided into the determination of the sphere position and the buried depth. There are mainly two strategies to determine the sphere position of RBFs. One is to directly place RBFs under observations, which is extremely dependent on the spatial distribution of observations. The other is to place RBFs on the regular grid, which can make RBFs evenly distributed in the computation areas. At present, the second strategy is more commonly used, which is also adopted in this article. The simplest regular spherical grid is a geographical grid, which is one kind of isogonal grids. Geographical

grid nodes space along the longitude circle and latitude circle with equal-angle intervals. Figure 10 shows its distribution globally and in the Arctic.

**Figure 10.** Distribution of the geographical grid.

It can be clearly seen from Figure 11 that the geographical grid is very unevenly distributed in high-latitude regions. The higher the latitude is, the denser the grid dots are. In a slightly large area, such a defect will become more obvious. In view of this phenomenon, scholars propose a kind of spherical isometric grid, the Reuter grid, which ensures that the spherical distances between grid dots in meridian and latitude direction are equal, so that RBFs can be distributed evenly in high-latitude regions, as shown in Figure 11. The specific formulas can be found in the relevant articles [32].

**Figure 11.** Distribution of the Reuter grid.

The determination of the optimal RBF network is always the focus and difficulty in RBF modeling. When RBFs are buried deep, the gravity signals simulated by the model are mainly concentrated at low frequencies. This kind of model has strong stability but leads to the simulated gravity signals being coarse due to omission errors. The shallower the buried depth of RBFs is, the more sensitive they are to the high-frequency signals such as terrain masses. The stability of the model will be reduced, and ill-conditioned problems may occur. When the number of RBFs is too large caused by high grid resolution, the overlap between RBFs will increase, which leads to the RBFs model being over-parameterized. On the contrary, too few RBFs are insufficient to simulate full gravity signals.

When the RBF model is constructed using discrete gravity points, it is difficult to uniquely determine the appropriate position of RBFs by theories before modeling. Scholars generally use posterior statistical methods such as GCV, RMS or STD minimization technique to screen out the optimal RBF network [27]. In this article, the STD minimization technique is adopted to design the RBF network. This method is simple and efficient, but it lacks strict a theoretical basis. The STD minimization technique divides a certain range of grid resolutions and depths into different combinations according to a certain step. These combinations are, respectively, used for modeling and compared with the control points. The combination achieving the smallest standard deviation is determined as the optimal RBF network.

#### 2.3.3. Tikhonov Regularization Technique

When directly using discrete data for modeling, the ill condition of the design matrix may occur due to the uneven distribution of observations and the excessive number of RBFs. Considering the lack of the sufficient prior information of observations, Tikhonov regularization is introduced to deal with the ill-conditioned problems [38]. Tikhonov regularization meets the following estimation criteria:

$$\min\_{\mathbf{x}} : \Phi(\mathbf{x}) = \|\mathbf{y} - A\mathbf{x}\|\_P^2 + \mathfrak{a} \|\mathbf{x}\|\_Q^2 \tag{8}$$

The general estimation formula is written as:

$$\pounds\_{\mathbb{K}} = \left(A^T P A + \kappa Q\right)^{-1} A^T P y \tag{9}$$

where *Q* is the regularization matrix; *α* is the regularization parameter. In this article, *Q* is assumed to be the identity matrix, i.e., *Q* = *I*. After singular value decomposition (SVD), the regularization solution can be obtained as:

$$
\overline{\chi}\_a = \sum\_{i=1}^n \frac{\sigma\_i^2}{\sigma\_i^2 + a} \frac{\mu\_i^T \overline{\mathcal{Y}}}{\sigma\_i} v\_i \tag{10}
$$

The regularization parameter is determined by the following two methods: the MSE and L-curve techniques. The basic principle of the MSE technique is as follows:

$$\min\_{\alpha} : \text{Tr}(MSE) = \sum\_{i=1}^{n} \frac{\sigma\_0^2 \sigma\_i^2 + a^2 z\_i^2}{\left(\sigma\_i^2 + a\right)^2} \tag{11}$$

where *zi* = *v<sup>T</sup> <sup>i</sup> x*. *x* denotes the truth value of unknown parameter *x*. We usually use the estimated result *x*ˆ to replace *x*, which will affect the optimality of regularization parameter *α* to some extent.

The L-curve technique is used to find the best one from different regularization parameters to achieve the optimal balance between the residual sum of squares *VTPV* and the smoothness of regularization function *x*ˆ*TQx*ˆ. Construct a curve with *μ*(*α*) as abscissa and *λ*(*α*) as ordinate. The point with maximum curvature (inflection point) corresponds to the optimal regularization parameter.

$$\mu(\mathfrak{a}) = \log(\|\mathfrak{y} - A\mathfrak{x}\|\_P), \ \lambda(\mathfrak{a}) = \log\left(\|\hat{\mathfrak{x}}\|\_{\mathbb{Q}}\right) \tag{12}$$

#### 2.3.4. Residual Terrain Model

The high-frequency gravity signals implied by the terrain mass have an important influence on local gravity field modeling, especially in mountainous areas [39]. RTM can be used to represent terrain effects [40–42]. The RTM technique is essentially the differences between the real terrain and the reference terrain, as shown in Figure 12. The real terrain surface usually denotes digital terrain model (DEM) with a high resolution. The reference terrain is relatively smooth with a lower resolution, which can be computed from the spherical harmonic terrain model such as RET2014 [43]:

$$z^{RET2014}(\phi,\lambda) = \sum\_{n=0}^{n\_{\text{max}}} \sum\_{m=0}^{n} \left(\overline{HC}\_{nm} \cos m\lambda + \overline{HS}\_{nm} \sin m\lambda\right) \overline{P}\_{nm}(\sin \phi) \tag{13}$$

where *nmax* = 2160; *HCnm*, *HSnm* denotes full normalized height coefficients. Based on Forsberg's classic TC program, RTM results can be well calculated.

**Figure 12.** Residual terrain model.

#### **3. Results**

This section describes in detail the experiment process and the results of calculating the unified land–ocean quasi-geoid from heterogeneous data sets based on RBFs. Compared with the results of the Stokes integral, the advantages of RBFs in multi-source data fusion are proved.

#### *3.1. The East Coast Experiment Area of the USA*

Based on the RCR technique, EIGEN-6C4 (degree 2190) gravity anomalies are removed from the measured terrestrial, shipborne and DTU15 data, as shown in Figure 13. Detailed statistics of the residual gravity anomalies are shown in Table 1, which shows that the residual terrestrial gravity signals are the smoothest. The amplitude of shipborne residual gravity anomalies is the largest, whose STD reaches 5.261 mGal. The terrain of the experiment area is quite flat, so the impact of terrain masses can be ignored, seen from the value of *Δgterrestrial* − *ΔgEIGEN*−6*C*<sup>4</sup> and *Δgterrestrial* − *ΔgEIGEN*−6*C*<sup>4</sup> − *ΔgRTM* in Table 1.

**Figure 13.** Distribution of residual gravity anomalies in the east coast experiment area: (**a**) terrestrial and shipborne; (**b**) DTU15.



Firstly, we use terrestrial and DTU15 gravity anomalies for modeling, considering the shipborne data in this zone are few and unevenly distributed. Set the grid resolution range to 0.1◦~0.9◦ and the step size be 0.1◦; set the depth range to 10~50 km and the step size be 1 km. The (part of) condition number of the design matrix is counted as shown in Figure 14. It can be seen from the figure that the lower the grid resolutions and the deeper the depths are, the smaller the condition number is. Otherwise, the condition number will be larger, making the ill condition more serious.

**Figure 14.** Condition number of the design matrix.

Taking the combination with resolution of 0.1◦ and depth of 11 km as an example, the number of RBFs is 2324 and the condition number of the design matrix is 4.35 × 106. The accuracy of the ill-conditioned model is 2.454 × 104 m checked by control points. Use the MSE method and the L-curve method, respectively, to determine the regularization parameter, as shown in Figure 15. The results show that the optimal regularization parameters determined by the MSE and L-curve methods are all 1.995 × <sup>10</sup><sup>−</sup>4, and the accuracy of the geoid calculated using this parameter is 6.9 cm, which is a great improvement compared with the original modeling result.

**Figure 15.** Determination of the optimal regularization parameter: (**a**) MSE; (**b**) L-curve.

The calculation results of all combinations of resolutions and depths are shown in Figure 16. Figure 16a shows the accuracy of original modeling without regularization; the accuracy of RBF modeling using Tikhonov regularization is shown in Figure 16b. As can be seen in Figure 16, when the condition number is greater than 1000, the serious ill-conditioned problems lead to the gradual deviation of the modeling results from the normal values. In particular when the grid resolution is over 0.3◦, the accuracy of the model geoid is lower than the meter level. Using the Tikhonov regularization technique, the accuracy of the ill-conditioned model is greatly improved to the centimeter level, which


improves that Tikhonov regularization can effectively solve ill-conditioned problems in the RBF model.

**Figure 16.** Accuracy of RBF modeling: (**a**) without regularization; (**b**) with regularization.

Then shipborne gravities are added to terrestrial and DTU15 data sets. The weight ratio of terrestrial, DTU15 and shipborne data is 1:0.9410:0.0989 determined by the Helmert VCE technique. Fuse the three data sets to construct the RBF model with a grid resolution of 0.4◦ and a depth of 45 km, calculating regional gravity field at airborne points.

The distribution of airborne gravity disturbances after removing EIGEN-6C4 values is shown in Figure 17a. The fluctuation of airborne gravity signals is relatively stable. Gravity disturbances *δgA* at airborne points are calculated based on the network A determined by terrestrial, DTU15 and shipborne data. Then *δgA* are removed from the residual gravity in Figure 17a and as input data for the RBF modeling, as shown in Figure 17b. The statistical information of the two residual gravity disturbances is shown in Table 2.

**Figure 17.** Distribution of residual airborne gravity disturbances in the east coast experiment area: (**a**) *δgairborne* − *δgEIGEN*−6*C*4; (**b**) *δgairborne* − *δgEIGEN*−6*C*<sup>4</sup> − *δgA*.

**Table 2.** Residual airborne gravity disturbances in the east coast experiment area (mGal).


Based on the STD minimization technique, the network B is determined with a resolution of 0.8◦ and a depth of 16 km using *δgairborne* − *δgEIGEN*−6*C*<sup>4</sup> − *δgA* as input data. The spatial distribution of networks A and B is shown in Figure 18, where the bottom blue dots

refer to the network A; the red dots in the middle refer to the network B; the top yellow dots refer to the geoid grid points computed by the RBF model.

**Figure 18.** Distribution of networks A and B.

The RBF model is then constructed based on the network B, which is added together with the modeling result calculated by terrestrial, DTU15 and shipborne data based on the network A to obtain the final quasi-geoid, as shown in Figure 19a.

**Figure 19.** The unified land–ocean quasi-geoid from heterogeneous data sets in the east coast experiment area: (**a**) terrestrial + shipborne + DTU15 + airborne gravity; (**b**) terrestrial + shipborne + airborne gravity.

The quasi-geoid calculated by the RBF model is essentially a gravimetric geoid. It has a different datum than the GPS/leveling geoid, resulting in significant systematic errors. We use simple polynomial fitting to convert the datum of the gravimetric geoid to GPS/leveling geoid. The formula is as follows:

$$
\Delta N = a\_0 + a\_1(\varphi - \varphi\_m) + a\_2(\lambda - \lambda\_m) + a\_3(\varphi - \varphi\_m)^2 + a\_4(\varphi - \varphi\_m)(\lambda - \lambda\_m) + a\_5(\lambda - \lambda\_m)^2 + \cdots \tag{14}
$$

where *α*0, *α*1, *α*<sup>2</sup> ··· denote fitting coefficients, i.e., bias parameter and tilt parameter; *ΔN* denotes the differences between the gravimetric geoid and the GPS/levelling geoid. (*ϕm*, *λm*) denotes the center longitude and latitude. The order of the formula can be taken very high, but generally two orders is enough.

Considering the DTU15 data have already been introduced by the EIGEN-6C4 model, we also calculate the results without DTU15 gravity data, shown in Figure 19b. The calculation methods of the two geoids are the same.

The accuracy of the two quasi-geoids shown in Figure 19 is checked by the GPS/leveling data, respectively, on inland and coast, as shown in Table 3. On inland, the accuracy of the quasi-geoid fusing terrestrial, shipborne, DTU15 and airborne data is 1.9 cm, which is 1.3 cm on sea. After the DTU15 gravities are removed, the modeling accuracy is 1.9 cm inland and 1.2 cm on coast, not much different from the previous result.


**Table 3.** Accuracy of the unified land–ocean quasi-geoid in the east coast experiment area (cm).

#### *3.2. The West Coast Experiment Area of the USA*

The west coast of the USA has quite different topographic features from the east coast. Affected by the overall terrain high in the west and low in the east, the mountains in the west extend to the coastal zone, making the land–ocean boundary areas rugged, which brings many difficulties to the geoid calculation. Based on the RCR technique, EIGEN-6C4 gravity anomalies are removed from the measured terrestrial and shipborne gravity, as shown in Figure 20a. It can be seen from the figure that after the EIGEN-6C4 value is removed, shipborne gravity signals become relatively flat. There are serious omission errors in EIGEN-6C4 on terrestrial gravity points affected by the topographic relief. This article calculates RTM based on the prism integral to simulate high-frequency gravity signals and compensate for the omission errors existing in GGM.

**Figure 20.** Distribution of residual terrestrial and shipborne gravity anomalies in the west coast experiment area: (**a**) without RTM; (**b**) with RTM.

In the calculation of RTM, the DEMs are usually expanded by approximately 2◦ from the edge of the gravity data area to avoid edge effects, as shown in Figure 21. The calculation efficiency can be improved by setting inner and outer computing regions. For example, the inner region within a 100 km radius uses high-resolution DEMs and the outer region within a 200 km radius uses low-resolution DEMs. The impact of terrain mass outside the outer region can be ignored due to the oscillating nature of RTM elevations [42].

**Figure 21.** DEMs for the RTM technique: (**a**) SRTM15+; (**b**) RET2014.

The residual gravity anomalies obtained by removing RTM from the terrestrial observations are shown in Figure 20b, and their STD is reduced from 15.983 mGal to 4.134 mGal, indicating that the compensation for EIGEN-6C4 omission errors by RTM reaches approximately 74%. The residual DTU15 gravity anomalies are shown in Figure 22. The statistical information of various residual gravity data is shown in Table 4, which indicates that the variation range of residual DTU15 gravity anomalies is the smallest. The variation amplitude of residual terrestrial gravity anomalies is close to shipborne data.

**Table 4.** Residual gravity anomalies in the west coast experiment area (mGal).


Firstly, terrestrial and shipborne gravity data are used for modeling. The weight ratio of the two observations is 1:1.004 computed by the Helmert VCE technique, which shows that their accuracy level is close. Based on the STD minimization technique, the optimal RBF grid resolution is determined as 0.6◦ and the optimal depth is 22 km. Then we add residual DTU15 gravity anomalies to the terrestrial and shipborne gravity to construct the RBF model with an optimal grid resolution of 0.8◦ and an optimal depth of 26 km.

The distribution of airborne gravity disturbances after removing the EIGEN-6C4 value is shown in Figure 23a. The variation amplitude of residual airborne gravities is relatively small. Based on the RBF network A determined by terrestrial, shipborne and DTU15 data, the model gravity disturbances are calculated on airborne points. Then *δgA* are removed from the residual gravities in Figure 23a and as input data for the RBF modeling, as shown in Figure 23b. The statistical information of the two kinds of residual gravity disturbances is shown in Table 5.

**Table 5.** Residual airborne gravity disturbances in the west coast experiment area (mGal).


**Figure 23.** Distribution of residual airborne gravity disturbances in the west coast experiment area: (**a**) *δgairborne* − *δgEIGEN*−6*C*4; (**b**) *δgairborne* − *δgEIGEN*−6*C*<sup>4</sup> − *δgA*.

Based on the STD minimization technique, network B is determined with a resolution of 0.8◦ and a depth of 16 km. The RBF model is then constructed based on the network B, which is added with the modeling result based on the network A to obtain the final quasi-geoid, as shown in Figure 24a. We also calculate the results without DTU15 gravity data, shown in Figure 24b.

**Figure 24.** The unified land–ocean quasi-geoid from heterogeneous data sets in the west coast experiment area: (**a**) terrestrial + shipborne + DTU15 + airborne gravity; (**b**) terrestrial + shipborne + airborne gravity.

Check the accuracy of the three quasi-geoids shown in Figure 24 by 64 GPS/leveling points inland and 36 GPS/leveling points on coast, as shown in Table 6. The accuracy of the RBF model quasi-geoid is 2.2 cm inland, which is 2.1 cm on coast. After the DTU15 data are removed, the modeling accuracy is 2.0 cm inland and 2.1 cm on coast.

**Table 6.** Accuracy of the unified land–ocean quasi-geoid in the west coast experiment area (cm).


According to the results of the two experiments on the east and west coast, we consider that the accuracy of the unified land–ocean quasi-geoid can be improved by fusing heterogeneous data sets. Fusing terrestrial, shipborne, DTU15 and airborne gravities based on a multi-layer RBF network achieves great modeling accuracy both inland and on coast.

#### **4. Discussion**

Due to the complex topography of land–ocean junction areas, terrestrial and shipborne gravity measurements cannot be fully carried out. As it is also affected by the terrain, the accuracy of satellite altimetry data inshore is reduced. The airborne gravity is not limited by the terrain, but its accuracy will be lost with downward continuation. Therefore, the important premise of constructing a high-precision unified land–ocean quasi-geoid is to effectively fuse heterogeneous data sets. Considering the shortcomings of the commonly used integral algorithm, we use the RBFs to calculate the unified land–ocean quasi-geoid, the accuracy of which is improved by fusing more types of data sets. However, there are still some problems to be discussed and further studied, as follows:


good modeling results by screening lots of RBF networks, but it increases too much redundant calculations, which hinders the solution efficiency of the RBF model to a great extent. So, it is very important to develop a rigorous and logical method to determine the RBF networks.

(3) The accuracy of the quasi-geoid fusing terrestrial, shipborne and DTU15 data is quite high, but the improvement is not obvious after adding airborne gravity. The main reason is that the gravity signals on airborne points simulated by network A is insufficient, resulting in little change in residual gravity disturbances after removing *δgA*, as shown in Tables 2 and 5. Theoretically, if the result of *δgB* = *δg* − *δgGGM* − *δgA* is significantly reduced, the function of the network B will be more obvious. In the future, we can further improve the multi-layer RBF network and try to set more layers of RBF grids, simulating the gravity signals more accurately.

**Figure 25.** Behavior of RBFs in the spectral domain: (**a**) the IMQ kernel; (**b**) the Poisson kernel; (**c**) radial multipoles of order 1; (**d**) Poisson wavelets of order 1.

#### **5. Conclusions**

This article designs a multi-layer RBF model to construct the unified land–ocean quasigeoid fusing the measured terrestrial, shipborne, satellite altimetry and airborne gravity data in coastal areas. Several core problems in the process of RBF modeling are studied in depth. The local quasi-geoid with a 1 resolution is calculated, respectively, on the flat east coast and the rugged west coast of the United States. Several conclusions are summarized as follows:

(1) The behavior of four types of RBFs—the IMQ kernel, the Poisson kernel, radial multipoles and Poisson wavelets—is analyzed in the spatial domain. The figures show that RBFs have significant localization characteristics in the spatial domain, which is helpful to concentrate more gravity signals in local gravity field approximation. Placing RBFs on the geographic or the Reuter grid, the optimal RBF network, i.e., the optimal grid resolution and depth can be effectively determined based on the STD minimization technique.


Through the above research results, we can roughly summarize the advantages of RBFs compared with other geoid calculation methods: (a) RBFs are simple and easy for multi-source data fusion; (b) RBFs can directly deal with discrete observations without gridding or downward continuation; (c) RBFs have a strong adaptive ability. In general, the RBF model has an important application value to calculate the unified land–ocean quasi-geoid from heterogeneous data sets.

**Author Contributions:** Conceptualization, Y.L. and L.L.; methodology, Y.L.; funding acquisition and project administration, L.L.; investigation and software, Y.L.; supervision and validation, L.L.; writing—original draft, Y.L.; writing—review and editing, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 41874004.

**Data Availability Statement:** All data for this paper are properly cited and referred to in the reference list.

**Acknowledgments:** The authors are very grateful to the NGS, the DTU, the NOAA and the ICGEM for providing gravity, GPS/leveling, geoid and GGM data. The authors express thanks to Hirt and Tozer for providing DEM data. The software Gravsoft shared by Forsberg also helped a lot.

**Conflicts of Interest:** The authors declare no conflict of interest.
