Regional Gravity Field Modeling Using Band-Limited SRBFs: A Case Study in Colorado

1. Introduction

Gravity data is crucial for studying the Earth's gravity field and for establishing the International Height Reference System (IHRS). To evaluate the repeatability of gravity potential values and compare different gravimetric geoid models, the Joint Working Group (JWG) 2.2.2 of the International Association of Geodesy (IAG) initiated the Colorado geoid experiment in 2017.The experiment involved the use of various modeling methods and the comparison of corresponding geoid models. In this context, the National Geodetic Survey (NGS) provided terrestrial gravity data , airborne gravity data, GPS/leveling data, and a digital elevation model (DEM) for a mountainous region in Colorado, which allowed the comparison of different methods for geoid computation using the same input dataset. Currently, global gravity field models (GGMs) are routinely used to provide long wavelength information of gravity fields, while other data types provide medium and short wavelengths. In addition to GGMs, terrestrial gravity data are the most widely used data type in gravity field modeling. However, in areas that are challenging to access (e.g., polar regions, mountains, forests, coastal areas, oceans) or in very large areas exceeding 10,000 km², airborne gravity almost becomes the only viable option [1].

When using airborne gravimetry for the local geoid improvement, there are two main challenges that need to be addressed. Firstly, airborne gravity observations are taken at scattered points along flight trajectories, and a downward continuation technique must be employed [2-4]. However, the downward continuation process is an unstable procedure that amplifies high-frequency noise, leading to a reduction in the quality of airborne data. Secondly, the spectral range of airborne gravity data is uncertain. The flight altitude determines the resolution of the data, but the topographic heights beneath the flight lines often vary significantly. For example, in the mentioned study, the flight lines ranged from 5200 m to 7200 m above the ellipsoid, while the corresponding topographic heights varied from 932 m to 4376 m. This variation makes it challenging to obtain a reasonable resolution from the airborne gravity data. Secondly, airborne gravity data containing high-frequency noise and a series of filtering and estimation techniques are usually needed, which further exacerbate the uncertainty in the spectral information of airborne gravity data. To date, it is widely believed that the spectral range of airborne gravity data is approximately 2-3000, more detailed investigation is still necessary. It has been observed that using airborne gravity signals beyond their corresponding spherical harmonic (SH) degree can introduce modeling errors [5]. Therefore, finding and modeling this useful information provided by airborne gravimetry has become a critical issue in local geoid enhancement.

Several methods can be employed to model airborne gravity data for local geoid refinement, including the well-known least square collocation (LSC) method [6], least squares spectral combination method [7,8], and Slepian functions [9]. In this study, we have chosen to use band-limited spherical radial basis functions (SRBFs) for local geoid refinement. The justifications for this choice are as follows. First, the SRBF expression includes a continuation factor, allowing the establishment of observation equations directly at scattered points of observation without the need for gridding or downward continuation. This feature enables the assignment of observation errors to specific observables in space. Li utilized a set of band-limited SRBFs to model airborne gravity data and demonstrated that the SRBF method could effectively recover the gravity field while exhibiting filtering properties due to its design in the frequency domain [10]. Second, the band-limited SRBFs can be adjusted flexibly based on the spectrum information of the gravity data, which justifies their use in this study. This flexibility allows for better adaptation to the specific characteristics of the data. Third, SRBFs can express almost all types of gravity data, which form the basis for combining multisource gravity data to build a high-resolution regional gravity field. Additionally, experiments have shown that SRBFs are computationally efficient and do not require a covariance model like the LSC method, which is sometimes time-consuming[11]. Other examples and explanations of SRBFs can also be found in [12-21] and references therein.

In this study, our main focus is on refining the regional gravity field model in Colorado using band-limited SRBFs. To determine the optimal spectral information of the gravity data, we propose a residual and a-prior accuracy comparative analysis method. The organization of this work is as follows: This work is organized as follows: In Section 2, we introduce the study area in Colorado and describe the available data. We also outline the procedure for data preprocessing. In Section 3, we present the fundamental concepts of SRBFs and the parameter estimation procedure. This includes the formulation of observation equations, estimation of unknown coefficients, and calculation of gravitational functionals. In Section 4, we explain the computation procedure, which involves the RCR procedure, the residual and a-prior accuracy comparative analysis method, and how we combine different datasets. Section 5 presents our models as well as the validation of the results. Finally, Section 6 provides the conclusions and present potential future research directions.

2. Theory and Methods

A residual harmonic function ${\mathit{T}}_{\mathbf{r}\mathbf{e}\mathbf{s}}\left(\mathit{x}\right)$ outside a sphere of radius R can be written as a finite sum of SRBFs:

$$T_{res}\left(\mathit{x}\right)=\sum _{i=n}^N{\alpha }_i{\Psi }_i\left(\mathit{x},{\mathit{x}}_i\right)$$

(1)

where $\mathit{x}=r\bullet \mathit{r}=r{\left[cos\phi cos\lambda ,cos\phi sin\lambda ,sin\phi \right]}^{\mathrm{T}}$ is the position vector of the observation point P, $\lambda$ is the spherical longitude, $\phi$ is the spherical latitude, and $r=R+h$, with $h$ being the spherical height of P above sphere ${\mathsf{\Omega }}_R$ with radius R.${\mathit{x}}_i=R\bullet {\mathit{r}}_i$ are the position vectors of the SRBFs. $N$ is the number of SRBFs, $T_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)$ is the residual disturbing potential, and ${\alpha }_i$ is the unknown coefficient. $\Psi \left(\mathit{x},{\mathit{x}}_{\mathrm{i}}\right)$ is an SRBF, the specific expression for which is [22,23]

$$\Psi \left(\mathit{x},{\mathit{x}}_i\right)=\sum _{n_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{n_{\mathrm{m}\mathrm{a}\mathrm{x}}}{k}_n\frac{2n+1}{4\pi }{\left(\frac{R}{r}\right)}^{n+1}{P}_n\left({\mathit{r}}^{\mathrm{T}}{\mathit{r}}_i\right)$$

(2)

where $k_n$ is the kernel function that defines the spectral properties of the SRBFs,$R$ is the mean radius of Earth,$P_n$ is the Legendre polynomial of degree $n$, and${\mathit{r}}^{\mathrm{T}}{\mathit{r}}_i$ is the spherical distance between $\mathit{x}$ and ${\mathit{x}}_i$. $n_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and $n_{\mathrm{m}\mathrm{a}\mathrm{x}}$ are the minimum and maximum degrees of expansion, respectively. The SRBF in Eq. (2) is band-limited since the kernel functions are zero for each degree beyond $n_{\mathrm{m}\mathrm{a}\mathrm{x}}$ or below $n_{\mathrm{m}\mathrm{i}\mathrm{n}}$.

This general expression (Eq. (2)) needs to be adapted for describing different gravitational functionals. We use observations that are given in terms of gravity disturbances $\delta g$, which can be expressed as the gradient of the disturbing potential $T$. In spherical approximation, the magnitude of the gravity disturbance can be written as [24]

$$\delta g=-\frac{\partial T}{\partial r}$$

(3)

From the SRBF expression for the disturbing potential (Eq. (1)), the residual terrestrial gravity disturbances at the Earth’s surface and the residual airborne gravity disturbances at flight altitude can be written as

$$\left\{\begin{array}{c}{\delta g}_{\mathrm{r}\mathrm{e}\mathrm{s}}^{\mathrm{t}\mathrm{e}\mathrm{r}}\left({\mathit{x}}_{\mathrm{t}\mathrm{e}\mathrm{r}}\right)=\sum _{\mathrm{i}=1}^N{\alpha }_i\sum _{n_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{n_{\mathrm{t}\mathrm{e}\mathrm{r}}}{k}_n\frac{\left(n+1\right)}{r_{\mathrm{t}\mathrm{e}\mathrm{r}}}\frac{\left(2n+1\right)}{4\pi }{\left(\frac{R}{r_{\mathrm{t}\mathrm{e}\mathrm{r}}}\right)}^{\mathrm{n}+1}{\mathrm{P}}_n\left({\mathit{r}}_{\mathrm{t}\mathrm{e}\mathrm{r}}^{\mathrm{T}}{\mathit{r}}_i\right)=\sum _{i=1}^N{\alpha }_i{\Psi }_{\delta g}\left({\mathit{x}}_{\mathrm{t}\mathrm{e}\mathrm{r}},{\mathit{x}}_i\right)\\ {\delta g}_{\mathrm{r}\mathrm{e}\mathrm{s}}^{\mathrm{a}\mathrm{i}\mathrm{r}}\left({\mathit{x}}_{\mathrm{a}\mathrm{i}\mathrm{r}}\right)=\sum _{\mathrm{i}=1}^N{\alpha }_i\sum _{n_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{n_{\mathrm{a}\mathrm{i}\mathrm{r}}}{k}_n\frac{\left(n+1\right)}{r_{\mathrm{a}\mathrm{i}\mathrm{r}}}\frac{\left(2n+1\right)}{4\pi }{\left(\frac{R}{r_{\mathrm{a}\mathrm{i}\mathrm{r}}}\right)}^{\mathrm{n}+1}{P}_n\left({\mathit{r}}_{\mathrm{a}\mathrm{i}\mathrm{r}}^{\mathrm{T}}{\mathit{r}}_i\right)=\sum _{i=1}^N{\alpha }_i{\Psi }_{\delta g}\left({\mathit{x}}_{\mathrm{a}\mathrm{i}\mathrm{r}},{\mathit{x}}_i\right)\end{array}.\right.$$

(4)

where ${\mathit{x}}_{\mathrm{t}\mathrm{e}\mathrm{r}}$ and ${\mathit{x}}_{\mathrm{a}\mathrm{i}\mathrm{r}}$ represent the position vectors of terrestrial gravity disturbances and airborne gravity disturbances, respectively. $n_{\mathrm{t}\mathrm{e}\mathrm{r}}$ and $n_{\mathrm{a}\mathrm{i}\mathrm{r}}$ are the maximum degrees of expansion of the terrestrial gravity disturbances and airborne gravity disturbances, respectively.${\Psi }_{\delta g}\left(\mathit{x},{\mathit{x}}_i\right)$ is the adapted SRBF of $\delta g\left(\mathit{x}\right)$.

For terrestrial gravity observations, it is difficult to obtain a more appropriate expansion degree due to their uneven distribution. For airborne gravity data, however, due to their high altitudes and the low-pass filtering process, it is also not a simple task to determine optimal spectral information either. In this study, we provide a way to determine the optimal degree of series expansion for gravity data ($n_{\mathrm{t}\mathrm{e}\mathrm{r}}$ and $n_{\mathrm{a}\mathrm{i}\mathrm{r}}$), which will be introduced in Section 3.1.

In Eq. (4), terrestrial gravity disturbances and airborne gravity disturbances have the same SRBF coefficients ${\alpha }_i(i=1,\dots ,N)$, and we can thus combine these expressions to form our observation equations (Gauss–Markov model)

$$\left[\genfrac{}{}{0pt}{}{{\mathit{y}}_1}{{\mathit{y}}_2}\right]=\left[\genfrac{}{}{0pt}{}{{\mathit{A}}_1}{{\mathit{A}}_2}\right]\mathsf{\alpha }+\left[\genfrac{}{}{0pt}{}{{\mathit{e}}_1}{{\mathit{e}}_2}\right]$$

(5)

where ${\mathit{y}}_1$ and ${\mathit{y}}_2$ are the vectors of the terrestrial gravity disturbances and airborne gravity disturbances, and ${\mathit{A}}_1$ and ${\mathit{A}}_2$ are the design matrices.$\mathsf{\alpha }$ are the unknown SRBF coefficients, and ${\mathit{e}}_1$ and ${\mathit{e}}_2$ are the vectors of stochastic observation errors with expectation $\mathrm{E}\left\{\mathit{e}\right\}=0$ and dispersion $\mathrm{D}\left\{\mathit{e}\right\}=\mathit{C}$. $\mathit{C}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{C}}_1,{\mathit{C}}_2)$ is the variance‒covariance matrix of the observed error and meets ${\mathit{C}}_l={\sigma }_l^2{\mathit{P}}_l^{-1}(l=\mathrm{1,2})$, ${\mathit{P}}_1$ and ${\mathit{P}}_2$ is the observation weight matrix for group $l(l=\mathrm{1,2})$.

Usually, the normal equation matrix of Eq. (5) is ill conditioned and a stabilization is needed to make the computation of a solution possible. Moreover, since the true accuracy of the two datasets is not exactly known, it is especially necessary to determine reasonable weights when building the gravity field model by combining terrestrial and airborne gravity data. All of these requirements can be met using the variance component estimation method (VCE) [25]

$$\left[{\frac{1}{{\sigma }_1^2}{\mathit{A}}_1^T\mathit{P}}_1{\mathit{A}}_1+{\frac{1}{{\sigma }_2^2}{\mathit{A}}_2^T\mathit{P}}_2{\mathit{A}}_2+\frac{1}{{\sigma }_{\mu }^2}\mathit{I}\right]\mathit{\alpha }={\frac{1}{{\sigma }_1^2}{\mathit{A}}_1^T\mathit{P}}_1{\mathit{y}}_1+{\frac{1}{{\sigma }_2^2}{\mathit{A}}_2^T\mathit{P}}_2{\mathit{y}}_2$$

(6)

where ${\sigma }_1^2$ and ${\sigma }_2^2$ are the variance factors of terrestrial gravity disturbances and airborne gravity disturbances, respectively. ${\sigma }_{\mu }^2$ is the unknown prior variance factor, which is also a regularization parameter, and $\mathit{I}$ is the unit matrix. The variance factor ${\mathsf{\sigma }}_l^2$ is usually determined based on an iterative method. That is, the initial variance factors ${\sigma }_{l,0}^2$ and ${\sigma }_{\mu ,0}^2$ are first given, and then the initial SRBF coefficients are estimated with the least squares estimation method. Then, the residuals of the least square values of the observation groups are further obtained, and the new variance factor can be calculated by

$${\sigma }_{l,i}^2=\frac{{\widehat{\mathit{e}}}_{l,i}^T{\mathit{P}}_l{\widehat{\mathit{e}}}_{l,i}}{r_{l,i}}$$

(7)

$${\widehat{\sigma }}_{\mu ,i}^2=\frac{{\widehat{\mathit{e}}}_{\mu ,i}^T{\mathit{P}}_{\mu }{\widehat{\mathit{e}}}_{\mu ,i}}{r_{\mu ,i}}$$

(8)

where ${\widehat{\mathit{e}}}_{l,i}^T$and ${\widehat{\mathit{e}}}_{\mu ,i}$ are the residuals of the observation groups after the $\mathrm{i}$th iteration, and $r_{l,i}$ and $r_{\mu ,i}$ are redundant numbers, which can be expressed as

$$r_{l,i}=n_l-\frac{1}{{\sigma }_{l,i}^2}tr\left({\mathit{A}}_{l,i}^T{\mathit{P}}_l{\mathit{A}}_{l,i}{\mathit{N}}_i^{-1}\right)=n_l-\frac{1}{{\sigma }_{l,i}^2}tr\left({\mathit{N}}_{l,i}{\mathit{N}}_i^{-1}\right)$$

(9)

$$r_{\mu ,i}=u-\frac{1}{{\sigma }_{\mu }^2}tr\left({\mathit{P}}_{\mu }{\mathit{N}}_i^{-1}\right)$$

(10)

where $n_l$ represents the number of gravity observations in group $l$, $u$ is the number of coefficients of SRBFs, ${\mathit{A}}_{l,i}$ is the design matrix at the $\mathrm{i}$th iteration of the observations in group $l$, ${\mathit{N}}_{l,i}$ is the matrix of the normal equation at the $\mathrm{i}$th iteration of the observations in group $l$, and ${\mathit{N}}_i$ is the total normal equation matrix of all observations at the $i$th iteration. $\mathrm{t}\mathrm{r}$ stands for trace operation.

After the new variance factor is calculated, the next loop is entered until the convergence condition of the following equation is satisfied

$$\mathrm{m}\mathrm{a}\mathrm{x}\frac{{|\widehat{\sigma }}_{l,i}^2-{\widehat{\sigma }}_{l,i-1}^2|}{{\widehat{\sigma }}_{l,i}^2}\le \tau \left(l=\mathrm{1,2}\right)$$

(11)

where $\tau$ is the convergence threshold, and we set $\tau =0.01$. ${\widehat{\sigma }}_{l,i}^2$ is the variance factor of the $l$th group of observations after $i$th iteration.

Applying the least-squares method to Eq. (6), the unknown coefficients are estimated as

$$\widehat{\mathit{\alpha }}={\left[{\frac{1}{{\sigma }_1^2}{\mathit{A}}_1^T\mathit{P}}_1{\mathit{A}}_1+{\frac{1}{{\sigma }_2^2}{\mathit{A}}_2^T\mathit{P}}_2{\mathit{A}}_2+\frac{1}{{\sigma }_{\mu }^2}\mathit{I}\right]}^{-1}\left[{\frac{1}{{\sigma }_1^2}{\mathit{A}}_1^T\mathit{P}}_1{\mathit{y}}_1+{\frac{1}{{\sigma }_2^2}{\mathit{A}}_2^T\mathit{P}}_2{\mathit{y}}_2\right]$$

(12)

Inserting the coefficients of Eq. (12) back into Eq. (1), we can obtain the residual disturbing potential $T_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)$.

Residual height anomaly is related to the residual disturbing potential $T_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)$ by

$${\zeta }_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)=\frac{T_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)}{\gamma }$$

(13)

where $\gamma$ is the normal gravity and ${\zeta }_{\mathrm{r}\mathrm{e}\mathrm{s}}\left(\mathit{x}\right)$ is the residual height anomaly.

Then, during the recovery procedure, we add back the height anomaly component of the GGM and the topographic effects, and a full-wavelength height anomaly $\zeta \left(\mathit{x}\right)$ can be obtained.

Finally, the geoid undulation can be obtained by the following conversion:

$$N\left(\mathit{x}\right)=\zeta \left(\mathit{x}\right)+\frac{{\Delta g}_B\left(\mathit{x}\right)}{\gamma }H,$$

(14)

where $H$ is the orthometric height, $N\left(\mathit{x}\right)$ is the geoid undulation, and ${\Delta g}_B$ is the simple Bouguer anomaly. Higher-order terms of the separation term [26] are not considered here to be consistent previously published procedures.

3. Gravity Data Processing and Analysis

3.1. Area of Interest

The Colorado geoid computation experiment is a collaborative effort involving the international geodetic community [27]. The study region is located between -110° and -102° longitude and between 35° and 40° latitude, primarily within the state of Colorado, USA. The area is characterized by rugged terrain in most of its part, which is evident in the topographic map displayed in Figure 1(a). The average elevation in this region is approximately 2,017 m. The highest peak reaches an elevation of 4,376 m, while the lowest point is at 932 m above sea level. In the central part of the study area lies the St. Louis Valley, which is an extensive depositional basin situated at high altitude. It covers an approximate area of 21,000 km². To the west of the St. Louis Valley are the San Juan Mountains, which is a rugged mountain range in the Rocky Mountains. The eastern part of the study area consists of a plateau where the terrain is relatively flatter compared to other regions. This geographical information is important in understanding the characteristics and challenges of the study area for refining the regional gravity field model in Colorado.

3.2. Datasets Used for Modeling

The available data for the Colorado geoid computation experiment include terrestrial and airborne gravity data, GSVS17 GPS/leveling data, and Global Geopotential Models (GGMs). The terrestrial gravity data were provided by the National Geodetic Survey (NGS) and consist of 59,303 values obtained from over 137 survey campaigns, covering the entire area of the Colorado geoid project (see Figure 1 (b)). Grigoriadis et al. studied the spatial density of the point values for the area, and the results found that only 21.25% of the total cells of the $1^{\mathrm{\text{'}}}$ grid contain at least 1 point value, and approximately 50% of the cells of the $2^{\mathrm{\text{'}}}$ grid contain at least one value [28]. Therefore, it is more reasonable to build a geoid model with a resolution of less than $2^{\mathrm{\text{'}}}\times 2^{\mathrm{\text{'}}}$.

The raw airborne data use the Gravity for the Redefinition of the American Vertical Datum (GRAV-D) Block MS05. The north-south track spacing is 10 km, while the east-west cross track spacing is 80 km. The mean flight altitude is approximately 6,200 m. Figure 1 (b) illustrates the distribution of the survey lines in the area. As shown in Figure 1 (b), the lines are relatively evenly distributed and cover most of the study area. However, the coverage of airborne and terrestrial gravity data does not overlap completely. The airborne gravity data cover the latitudes between $35°\le \mathsf{\phi }\le 38.3$ and longitudes between $-109°\le \mathsf{\lambda }\le -102°$, There are no airborne data available for latitudes between $38.3°<\mathsf{\phi }\le 40°$ and longitudes between $-110°\le \mathsf{\lambda }<-109°$ (Figure 1 (b)). Moreover, the raw gravity data from MS05 are subject to biases between different tracks. Therefore, data from MS05 were debiased using the maximum median technique [29] and then resampled to 1 Hz. Finally, a dataset consisting of 283,716 data points was obtained. In addition, Childers et al. define the Fourier wavelength of the smallest recoverable feature as [30]

$${\lambda }_f=1.54\times Z_c$$

(15)

where $Z_c$ is the height of the measurement point above the center of the spherical body that causes the anomaly. In this study, the average altitude of the aircraft is 6.2 km, and the average topographic height is approximately 2 km; then, the distance of $Z_c$ is approximately 4.2 km. By substituting $Z_c=4.2$ into Equation (15), the minimum resolvable wavelength is calculated to be 6468 m, corresponding to a spherical harmonic degree of approximately 3000. However, this is just a rough estimate of the resolution since the study area is mountainous, and the topographic heights vary significantly.

Additionally, a total of 222 benchmarks of GSVS17 (Geoid Slope Validation Survey 2017) profile data [31] were provided by the NGS. These benchmarks are spaced approximately 1.6 km along Highway US 160 from Durango to Walsenburg. The accuracy of the geoid undulations inferred from the GPS-based ellipsoidal heights and the leveling-based orthometric heights is estimated to be approximately 1.5 cm. In addition to the comparisons made for the GSVS17 locations, entire geoid grids (at a $1^{\mathrm{\text{'}}}$ resolution) were compared based on the mean of all submitted models to assess the model’s performance over the entire area. Furthermore, the reference frame used for 3-D positioning of all three types of data (terrestrial gravity, airborne gravity, and GSVS17) is the IGS Reference Frame 2008 (IGS08) [32] on the Geodetic Reference System 1980 (GRS80) ellipsoid [33].

3.3. Data Preprocessing

In the original terrestrial data, there are cases in which several gravity observations are located at the same position, for which we only use the first of these observations (1090 points deleted). However, if the points sharing the same location but had a height difference of more than 3 m or a gravity difference of 2 mGal, then these points were both removed. Hence, 1215 points in total were excluded from the dataset (see Figure 1 (b)). Finally, a total of 58088 terrestrial gravity observations in the modeling area were obtained. Moreover, since the measurements are provided in orthometric heights H, we transform them to ellipsoidal heights using the geoid model 'GEOID 18B' provided by the National Geodesy Survey.

The airborne data had a dense distribution, resulting in a large design matrix. To simplify the computations, the sampling interval was reduced from 1 Hz to 1/8 Hz, and an average spatial interval of approximately 1 km along-track was obtained. This resulted in a total of 31,527 airborne gravity observations in the modeling area.

Then, for both types of observations, the following data preprocessing steps are performed:

① The observations were converted from absolute gravity $g$ to gravity disturbance $\mathsf{\delta }\mathrm{g}$ by subtracting the normal gravity $\gamma$ at the ellipsoidal height $h$ of the observations

$$\delta g=g-\gamma$$

(16)

② The atmospheric correction for the gravity disturbance ${\delta g}_{\mathrm{A}\mathrm{T}\mathrm{M}}$ was calculated using the method proposed by Torge [34]

$${\delta g}_{\mathrm{A}\mathrm{T}\mathrm{M}}=0.874-9.9\times {10}^{-5}H\times 3.56\times {10}^{-9}{H}^2$$

(17)

③ The atmospheric correction was added to the observations.

$$\delta g={\delta g}_0+{\delta g}_{\mathrm{A}\mathrm{T}\mathrm{M}}$$

(18)

3.4. Remove Procedure

Gravity data obtained in a limited area cannot resolve long wavelengths in the gravity field. Therefore, the remove–compute–restore (RCR) procedure is usually utilized to remove the low-frequency and high-frequency signals of gravity data before modeling. However, when applying the RCR procedure, it is important to assess the performance of the Global Gravity Models (GGMs) used, as the error of the GGMs accumulates as expansion degrees increase. In this study, five combined GGMs are investigated: EGM2008 [35], XGM2019e_2159 [36], SGG-UGM-2 [37], XGM2016 [38], and GECO [39]. These models were chosen because they have relatively high accuracies and better spatial resolutions compared to other GGMs. The error degree variance offers insight into the error of a GGM, and it supplies a rudimentary quality assessment [40]. Figure 2 shows the degree-wise accumulated geoid errors of the five GGMs. EGM2008 exhibits a relatively large error, which rapidly increases from degrees 30-220, reaching approximately 7.3 cm at degree 220, and then slowly increases to 8.2 cm at the highest degree. On the other hand, models like XGM2019e_2159, SGG-UGM-2, and GECO, which have similar expansion degrees, show reduced errors of approximately 3.1 cm, 1.9 cm, and 4.3 cm, respectively. The significant error in EGM2008 in the frequency bands between degrees 30-220 is mainly due to the lack of data from the Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite. In contrast, the other three models are developed using GOCE data, resulting in better qualities in these frequency bands. Among the investigated models, XGM2016 exhibits the minimum cumulative geoid height error. Up to degree 719, the cumulative geoid height error of XGM2016 is approximately 1.1 cm, which is significantly lower than the corresponding values of the other four GGMs. Therefore, the XGM2016 model up to degree 719 is chosen for the RCR procedure. It is important to note that the error degree variances only provide a global mean of the internal error and cannot be considered as a realistic error estimate [41].

Another important factor that influences the performance of geoid modeling is topographic effects, which play a key role, especially in mountainous areas, since they can smooth the input observations and are of utmost importance for obtaining a good least-square fit. The topographic model dV_ELL_Earth2014 [42] from degree 720 to degree 2159 and a residual terrain model ERTM2160 [43] from degree 2160 to degree ~80,000 are removed from the terrestrial data. Similarly, for the airborne data, the effect of dV_ELL_Earth2104 from degree 720 to degree $n_{\mathrm{a}\mathrm{i}\mathrm{r}}$ is removed. It is worth noting that using different topographic models above degree 2160 for the terrestrial and airborne data is reasonable because both the dV_ELL_Earth2014 and ERTM2160 models are derived from the same original data and contain the same signal [44]. This approach ensures consistency and avoids introducing additional errors due to the use of different models. Notably, for airborne data, the effect of dV_ELL_Earth2104 from degree 720 to degree $n_{\mathrm{a}\mathrm{i}\mathrm{r}}$ is removed. In contrast, previous studies have utilized different degrees for removing topographic contributions in airborne data. However, we propose that removing topographic effects with degrees exceeding the optimal expansion degree of the gravity data may introduce modeling errors. Therefore, in this study, we remove the topographic contributions up to the expansion degree of the airborne data.

Figure 3 illustrates the distributions of the terrestrial and airborne gravity disturbances before and after the removal of the GGM and topographic effects. Figure 3(c) and 3(d) clearly show that the gravity disturbances become significantly smoother in the middle part of the study area after subtracting the contributions of the GGM and topographic effects.

Table 1 presents the statistical results of the gravity observations. The standard deviations of the terrestrial and airborne gravity disturbances are reduced from 38.72 to 6.90 mGal and from 29.56 to 3.35 mGal, respectively, after removing the contributions of the GGM and topographic effects. The standard deviation of the terrestrial data is reduced by 42.15% after subtracting the GGM and further reduced by 82.18% after including the topographic effects. Similarly, the airborne observations are reduced by 72.40% after subtracting the GGM and further reduced by 88.67% after including the topographic effects. The significant reduction in standard deviation indicates that the removal of the GGM and topographic effects leads to a smoother gravity field, enabling a better least-square fit in the gravity field modeling using SRBFs. It is worth noting that the smoothing effects are more pronounced in the terrestrial observations compared to the airborne observations after including the topographic model (Table 1). This is due to the fact that the short wavelength signal of the gravity field decreases with height, and the airborne gravity data are less sensitive to the short wavelength part compared to the terrestrial gravity data. Consequently, subtracting the topographic effects has a greater influence on the terrestrial gravity data compared to the airborne gravity data. Illustrating the smoothing effects and quantifying the reduction in standard deviation demonstrates the importance of removing the GGM and topographic contributions in order to improve the quality and accuracy of gravity field modeling.

4. Gravity Field Modeling

In the gravity field modeling based on band-limited SRBFs, several parameters and procedures were considered. The following configurations were used:

① Geocentric gravitational constant ($\mathrm{G}\mathrm{M}$): 3.986 004 415×${10}^{14}{\mathrm{m}}^3/{\mathrm{s}}^2$.

② Average density of topographic masses (${\mathsf{\rho }}_0$): 2670 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$.

③ Geoid computations were performed in the tidefree system. GRS80 was used as the reference ellipsoid, and the conventional constants are referred to Moritz (2000).

Apart from these fundamental parameters, four important factors that influence the performance of SRBFs were also taken into account, namely, (1) the type of the SRBFs, (2) the location of the SRBFs, (3) the extensions of the data area and (4) the bandwidth parameters. Considering these parameters and procedures helps to ensure an effective and accurate gravity field modeling process based on band-limited SRBFs. Next, we introduce the model configuration for gravity field modeling based on band-limited SRBFs.

4.1.1. The Type of The SRBFs

There are many types of SRBFs that can be used for gravity field modeling, such as the Shannon functions, the Blackman functions, the Poisson functions, and the cubic polynomial functions. The Shannon function is a reproducing kernel of the space spanned by all solid spherical harmonics of degrees ${n=n}_{\mathrm{m}\mathrm{i}\mathrm{n}}\cdots n_{\mathrm{m}\mathrm{a}\mathrm{x}}$ [45]. As a result, it does not smooth the gravity signal in this spectral band. However, other types of SRBFs (e.g., Blackman functions, Poisson functions and cubic polynomial functions), with nonconstant kernel functions might smooth some gravity field harmonics. Though some harmonic properties are still remaining in the recovered signal, they might be significantly suppressed. Therefore, one has to be careful when interpreting results based on the SRBFs with smoothing kernel functions. Moreover, Bentel et al. studied the performance of different types of SRBFs, and the results show that the Shannon SRBFs lead to the same results as with spherical harmonics both in theory and in practice. Therefore, in this study, the Shannon kernel functions are used, the specific form of which is

$$k_n=\left\{\begin{array}{c}1nϵ[n_{\mathrm{m}\mathrm{i}\mathrm{n}},n_{\mathrm{m}\mathrm{a}\mathrm{x}}]\\ 0\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e}\end{array}\right.$$

(19)

4.1.2. Type of Grids

The location of the SRBFs depends on the type and number of grid points. There are many types of grids that can be used for regional gravity field modeling, such as the Driscoll-Healy grid, the triangle vertex grid and the Fibonacci grid. Eicker analyzed the advantages and disadvantages of different types of grids, and the results show that both the Reuter grids [46] and the triangle vertex grids are suitable for gravity field modeling with SRBFs. Bentel et al. noted that different types of grids do not differ significantly, especially compared to the other three factors that influence the modeling results. In addition, the control parameter $\beta$ of the Reuter grid is related to the maximum expansion degree of SRBFs, which also facilitates gravity field modeling. We set ${\beta =n}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, where $n_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum degree of the series expansion of SRBFs.

4.1.3. The Extensions of the Target, Observation, and Computation Area

In regional gravity modeling, the extension of the target area $\partial {\Omega }_{\mathrm{T}}$, the observation area$\partial {\Omega }_{\mathrm{O}}$, and the computation area $\partial {\Omega }_{\mathrm{C}}$ needs to be defined carefully. At the edge of the observation area$\partial {\Omega }_{\mathrm{O}}$, the unknown coefficients cannot be appropriately estimated, which provokes edge effects. In general, the target area $\partial {\Omega }_{\mathrm{T}}$ should be smaller than the observation area $\partial {\Omega }_{\mathrm{O}}$ because the modeling accuracy gradually increases as the target area decreases. However, this operation reduces the target area, which sometimes resulting in a significant waste of gravity data. Moreover, the computation area $\partial {\Omega }_{\mathrm{C}}$ should be larger than the observation area$\partial {\Omega }_{\mathrm{O}}$ [47]. Otherwise, related computations may fail in gravity field modeling. Therefore, to minimize the edge effects in the computation, $\partial {\Omega }_{\mathrm{T}}\subset \partial {\Omega }_{\mathrm{O}}\subset \partial {\Omega }_{\mathrm{C}}$ should be satisfied. In addition, the margins $\eta$ need to be defined between the three areas. In our case, the target area $\partial {\Omega }_{\mathrm{T}}$ is given to be between −109$°$ and −103$°$ and 36$°$ and 39$°$. Therefore, the margin size ${\eta }_{\mathrm{T},\mathrm{O}}$ between $\partial {\Omega }_{\mathrm{T}}$ and $\partial {\mathsf{\Omega }}_{\mathrm{O}}$ is fixed to 1$°$

The determination of the margin size ${\eta }_{\mathrm{C},\mathrm{O}}$ between $\partial {\Omega }_{\mathrm{C}}$ and $\partial {\Omega }_{\mathrm{O}}$ follows the empirical formula described by Liu et al. [20]

$${\eta }_{\mathrm{C},\mathrm{O}}\approx \frac{360}{n_{\mathrm{m}\mathrm{a}\mathrm{x}}cos\left(\left|{\phi }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|\right)}$$

(20)

where ${\phi }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum latitude value of the target area. In this study, the value ${\eta }_{\mathrm{C},\mathrm{O}}\approx 0.1°$ follows, and the reason will be explained in section 4.1.4. Figure 4 visualizes the target area, the observation area, the computation area, and the margins.

4.1.4. Maximum and Minimum Degree of Expansion

The maximum degree of expansion of SRBFs is usually determined by high-resolution gravity data, which can be determined by empirical rules [48]

$$n_{\mathrm{m}\mathrm{a}\mathrm{x}}\le \frac{\pi R}{\mathrm{R}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{l}}$$

(21)

where Resol represents the resolution of the gravity data. In this study, the spatial resolution of terrestrial data is higher than that of airborne data. Therefore, the maximum expansion degree of SRBFs is mainly decided by the terrestrial data. Moreover, Liu et al. claimed that the average spatial interval of terrestrial gravity data is approximately 3.5 km, so they expanded the SRBFs to degree 5600, and signal frequencies higher than the maximum degree were not modeled. However, the maximum expansion degree determined by Liu et al. may be problematic, since Grigoriadis et al. stated that only about 50% of the cells in the $2^{\mathrm{\text{'}}}\times 2^{\mathrm{\text{'}}}$ grid contained at least one value, which suggests that the data distribution may not be suitable to support such a high maximum degree of expansion. Therefore, this current study suggests that better results may be achieved by expanding the terrestrial gravity data to degrees less than 5400. However, determining the optimal spectral range for both terrestrial and airborne gravity data is difficult due to factors such as high altitude and uneven distribution of observations. To address this issue, we propose the residual and a-priori accuracy comparative analysis method, which will be discussed in detail in section 4.2.

For the minimum degree of expansion, Bucha et al. noted that there is still useful gravity information below the cutoff degree of the Global Geopotential Model (GGM). They concluded that starting the series expansion of SRBFs at a lower degree, such as 0, can help capture this useful part of the signal. This method has been commonly used in local gravity field modeling. However, the current study notes that this method is not always effective. In this study, since the contributions of XGM2016 up to degree 719 has been subtracted from the original signals, we set $n_{\mathrm{m}\mathrm{i}\mathrm{n}}=720$. By using this minimum degree starting from 720 instead of 0, we obtained better results, at least in the Colorado region of this study.

4.2. Residual and A-prior Accuracy Comparative Analysis Method

To efficiently use the gravity data for gravity field modeling, the residual and a-prior accuracy comparative analysis method is proposed to determine the optimal expansion degree of the gravity data. The procedure is as follows.

Step 1: In a large search interval, certain degrees with a large step size are selected as the maximum degree of expansion for the data. Using the gravity data available, the terrestrial-only and airborne-only SRBF models are built and the coefficients of SRBFs are estimated using the Variance Component Estimation (VCE) method. The modeling residuals of the gravity data are then statistically analyzed. The degree that yields the minimum root mean square (RMS) values of the residual with respect to the a-priori accuracy is identified as the initial degree.

Step 2: A smaller range around the initial degree is defined with a smaller degree step size. New SRBF models are built within this range. The differences between these models and GSVS17 data are calculated to obtain the STD value for each model. Finally, the model that has the smallest STD value is determined, and the corresponding expansion degree is considered as the optimal expansion degree of the gravity data.

It is important to note that Step 1 helps save calculation time and improve efficiency. Due to the large amount of gravity data, calculating the STD value for the geoid height difference for each model would require significant computational time. By initially selecting a larger degree range, the computational effort is reduced. In Step 2, external data (GSVS17) is used to assess the quality of the SRBF model within a smaller range. This is because there is some uncertainty in both the a-priori accuracy of terrestrial and airborne data. By incorporating external data within a smaller degree range, the computational effort is further reduced.

4.3. Gravity Field Modeling with SRBFs

We combine the terrestrial data and the airborne data using the Gauss–Markov model (Eq. (5)). Saleh et al. studied the qualities of the terrestrial gravity database of NGS, and the results show that the accuracy of terrestrial data is approximately 2.2 mGal [49]. Moreover, Varga et al. estimated the accuracy of terrestrial gravity data assuming that the errors are uncorrelated and obtained an accuracy of 2.6 mGal [50]. In addition, other scholars have estimated terrestrial gravity data at 1 mGal [51] or 3 mGal [52,53]. However, it seems that the accuracy of 1 mGal may be overstated, while the accuracy of 3 mGal may be understated. Therefore, we set the a-prior accuracy of the terrestrial gravity data to 2.2 mGal, which is essentially consistent with the results of Saleh et al. Moreover, the a-priori accuracy of the airborne gravity data is approximately 1.6-2 mGal [53], and we take the median value of 1.8 mGal as its a-priori accuracy.

The determination of the optimal expansion degree of terrestrial and airborne gravity data is conducted using the residual and a prior accuracy comparative analysis method. Figure 5 illustrates the procedure for finding the optimal degree of expansion. To fully demonstrate the difference caused by the degree of expansion of SRBFs, the STD of the differences between the geoid heights of the SRBF model and the GSVS17-derived geoid heights were also calculated. Figure 5a (Step 1) shows that a search range of 3600 to 5400, with a step size of 300, is selected for the terrestrial gravity data. For the airborne gravity data, a search range of 1600 to 3000, with a step size of 100, is chosen. It is observed that when the expanding degrees of the terrestrial and airborne data are 5100 and 1800, respectively, the RMS values of the residuals closely match their a-priori accuracy. Hence, these degrees are set as the initial values of expansion for the terrestrial and airborne gravity data (depicted by the black line in Figure 5a and b).Then, in Step 2, a smaller search interval of 4800 to 5400, with a step size of 100, is selected for the terrestrial data. For the airborne data, a search interval of 1800 to 1850, with a step size of 10, is chosen. It is observed that when the expanding degrees of the terrestrial and airborne data are 5200 and 1840, respectively, the geoid height differences reach their minimum value (depicted by the red line in Figure 5a and b).

In addition, as shown in Figure 5, although the residuals of the gravity disturbances decrease with an increase in the degree of expansion, the corresponding geoid accuracy does not exhibit the same trend. In fact, the increase of the expansion degree only improves the agreement between the recovered signal and the original signal, but it can also lead to overfitting of the observed values when the degree exceeds the optimal expansion degree of the SRBFs. Thus, the quality of the geoid does not always improve with increasing expansion degree. Additionally, the modeling error of airborne gravity data rapidly increases after the expansion degree surpasses 2200, whereas the corresponding change in terrestrial gravity data is not as significant. This indicates that the expansion degree of airborne gravity data may be more sensitive than that of terrestrial gravity data. In addition, according to the estimated resolution formula given by Eq. (15), the spherical harmonic degree of airborne gravity data is approximately 3000, whereas the optimal degree detected in our study is 1840. The STD values of geoid height differences calculated using these two degrees are approximately 0.158 m and 0.030 m, respectively (Figure 5 (b)), resulting in a difference of approximately 0.13 m. This highlights the influence of the expansion degree of the airborne gravity data on the modeling results. Moreover, the resolution of the airborne data explored in this study is approximately 10.6 km, which is larger than that of the NGS at 6.5 km. However, in terms of external accuracy, the resolution determined for the airborne data in this study seems to be more reliable than the aforementioned approach.

5. Evaluation of the Combined Solution

In this section, we compute two sets of output gravity functionals. The first is the geoid heights at the GSVS17 benchmarks, at which fourteen groups have provided geoid height results. The second is a grid model from 36$°$ to 39$°$ and −109$°$ to −103$°$ with a resolution of $1^{\text{'}}\times 1^{\text{'}}$. Fourteen groups have provided the geoid height grid model, and the comparison between our models and the mean of all the models are given.

5.1. Comparing to Models with Different Expanding Degrees of SRBFs

According to Eq. (6), there are two variables that need to be determined: the relative weight $\omega$ between the two types of observations and the regularization parameter $\lambda$. Both of these variables are determined through the VCE method, which is an iterative process to estimate the variance factors of different observation types and the prior information. Figure 6 illustrates how the variance factors of the gravity data change during the iteration process.

As seen from Figure 6, the iteration starts with the initial values of ${{\sigma }_1^2=6.25\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$, ${\sigma }_2^2=3.24{\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$ and ${\sigma }_{\mu }^2=1\times {10}^{-13}$, and after approximately 4 iterations, the variance factors tend to stabilize, and convergence is reached. The result of the finally obtained variance factors is ${\sigma }_1^2=5.56{\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$, ${\sigma }_2^2=4.15{\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$, and ${\sigma }_{\mu }^2=1.17\times {10}^{-13}$, respectively. The relative weight $\omega$ between the airborne data and the terrestrial data is approximately 1.34:1, and the regularization parameter $\lambda$ is approximately 4752 upon converting mGal to $\mathrm{m}/s^2$($1\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}=1\times {10}^{-5}\mathrm{m}/s^2$). It is worth noting that this result differs from that of Liu et al., who obtained variance factors of ${\sigma }_1^2=6.13{\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$, ${\sigma }_2^2=1.61{\mathrm{m}\mathrm{G}\mathrm{a}\mathrm{l}}^2$, and ${\sigma }_{\mu }^2=8.17\times {10}^{-14}$. Such a result is possible if both the terrestrial gravity data and the airborne gravity data are expanded to the same higher degree of 5600 [20]. Since the standard deviation of the residual airborne gravity data is small, the fitting effect is definitely better than that of terrestrial gravity data. However, it should be noted that the estimated variance factors actually represent the quality of the least-squares fit rather than the accuracy of the data itself.

To assess how our geoid model benefits from the residual and a prior accuracy comparative analysis method, our optimal computed model is compared to the model obtained with different maximum and minimum degrees of expansion of gravity data. Each model was then compared to the GSVS17 validation data, and the differences between the geoid height results and the validation data were recorded in Table 2.

The results of the investigation revealed that the geoid models obtained by setting the expansion degrees to 5600 (Model I) or 5200 (Model III) for both terrestrial and airborne gravity data, with a minimum expansion degree of 0, performed the worst with a standard error (STD) of 3.1 cm. However, when the minimum degree of expansion was changed to 720 (Models II and IV), the STD reduced to 2.6 cm. This suggests that the accuracy of the XGM2016 model up to degree 719 is higher than the residual components of the gravity data. Consequently, including the residual signals from degrees 0 to 719 does not contribute to improving the accuracy of the geoid model. This finding highlights the importance of considering the minimum expanding degree in regional gravity field modeling based on SRBFs. We note that this factor is often overlooked in previous studies. The quality of the XGM2016 model and its contribution to the overall accuracy of the geoid model should be taken into account when determining the appropriate expanding degrees for gravity data.

The optimal solution for our geoid model is Model VI, which is obtained by setting the expansion degrees to $n_{\mathrm{t}\mathrm{e}\mathrm{r}}=5200$, $n_{\mathrm{a}\mathrm{i}\mathrm{r}}=1840$ and $n_{\mathrm{m}\mathrm{i}\mathrm{n}}=720$. When compared to the GSVS17 validation data, the geoid height differences have a STD value of approximately 2.3 cm, which is the smallest among all six selected models (Table 2). The improved performance of Model VI can be attributed to the introduction of a residual and a-priori accuracy comparative analysis method, as well as the special treatment of the minimum expansion degree of the SRBFs. Comparing Model VI to Model II and Model V, which have the same expansion degree of terrestrial data but different expansion degrees of airborne data, we find that Model II has an STD of 2.6 cm, while Model V has an STD of 2.5 cm, which are increased by 3 mm and 2 mm, respectively, compared to the STD value of Model VI. Therefore, when aiming to construct a geoid model with an accuracy of 1 cm, it is necessary to consider the potential impact of different expansion degrees of the gravity data.

5.2. GSVS17 Comparisons

To further evaluate the quality of the optimal geoid model (hereinafter referred to as ColSRBF2023), five geoid models are downloaded from the International Service for the Geoid (ISG) for accuracy comparison. The five models are ColFFTW2020 from Aristotle University of Thessaloniki (AUTh), ColSRBF2019 from the Deutsches Geodätisches Forschungsinstitut (DGFI), ColUNBSH2019 from the Geospatial Information Authority of Japan (GSI), ColRLSC2019 from the Institute for Astronomical and Physical Geodesy (IAPG) and ColWLSC2020 from the Politecnico di Milano (Polimi). The reason for choosing the five models is that all five models use the same RCR procedure. For example, they all use the same reference gravity field model (XGM2016 up to degree 719) and the same methods for removing the topographic effects [53]. Therefore, compared with the above models, they can better reflect the differences caused by different modeling methods. It is important to note that there is a bias of approximately 87 cm between the geopotential value of the reference level used in our experiment and that defined by the NAVD88 datum, which is tied to the tide gauge at Rimouski [54]. This bias has been subtracted in the comparison. Figure 7 shows the differences between the geoid heights obtained by different models and GSVS17-derived geoid heights, and Table 3 lists the geoid height differences of each model along the GSVS17 profile.

Figure 7 shows that the differences between each model and GSVS17 generally follow a similar trend. The STD values of the geoid height differences range from 2.3 cm to 3.9 cm, as listed in Table 3. The model with the smallest STD value is ColSRBF2023, which has a difference of approximately 2.3 cm. Table 3 also provides the range of differences for each model relative to GSVS17. The differences range from 11.3 cm for the ColSRBF2023 model to 22.8 cm for ColWLSC2020, with an average value of 15.7 cm. These values are quite significant considering the target accuracy of 1 cm for geoid modeling. However, considering the uneven distribution of observations, the result is still quite satisfactory.

In addition, as shown in Figure 7, the differences are clear around the point numbers (Nos.) of 100 and 200, with maximum and minimum deviations of 7.9 cm and -14.9 cm, respectively. This may be attributed to the scarcity of terrestrial gravity data in this section (see Figure 1). Notably, around Nos. 100 and 200, the terrain changes are more significant, and the ColSRBF2023 model is significantly improved compared with the other five models, which indicates the effectiveness of the residual and a-prior accuracy comparative analysis method proposed in this study. These values are quite significant considering the target accuracy of 1 cm for geoid modeling.

5.3. Area Comparison of Geoid Grids

Since GSVS17 is located in a relatively flat area, comparison with GSVS17 may misrepresent the error of the entire study area. Therefore, we employ the average value of the geoid models submitted to ISG as a calibration model to further evaluate the accuracy of ColSRBF2023. Specifically, the geoid heights of the target area (36$°$N-39$°$N, 109$°$W-103$°$W) with a resolution of $1^{\text{'}}\times 1^{\text{'}}$ are first calculated by ColSRBF2023, and then the geoid heights are compared with the calibration model. Figure 8a visualizes the geoid model for the whole target area ${\partial \Omega }_{\mathrm{T}}$, with a grid resolution of $1^{\text{'}}\times 1^{\text{'}}$, and Figure 8b shows the geoid height differences w.r.t. the calibration geoid model. Table 4 provides the statistics of these differences.

From Figure 8b, we can see that large differences are mainly concentrated on the east and west sides of the target area (dark blue and dark red blocks in Figure 8b). These areas are characterized by drastic elevation variations and sparse distribution of gravity data, highlighting the importance of gravity data density and quality. In other areas (yellow and light blue regions in Figure 8b), the differences are relatively small, generally not exceeding ±5 cm. Furthermore, in terms of RMS of the differences, ColSRBF2023 has a value of 2.4cm, while the other three models are 3.0cm, 3.0 cm and 2.9 cm, respectively. ColSRBF2023 has the smallest RMS value among all the models compared, and the accuracy improved by 5-6 mm. Additionally, the model with the largest range of differences is ColUNBSH2019 (GSI), reaching 48.1 cm. ColSRBF2023 performs slightly better, with an average range of 31.8 cm (Table 4). The overall average range of differences for all models is 41.8 cm, which reflects the significant challenge of achieving 1 cm accuracy in geoid modeling in mountainous areas.

In general, the optimal model in this study is the combined geoid model ColSRBF2023 obtained by setting the expansion degree of terrestrial and airborne gravity data $n_{\mathrm{t}\mathrm{e}\mathrm{r}}=5200$, $n_{\mathrm{a}\mathrm{i}\mathrm{r}}=1840$ and $n_{\mathrm{m}\mathrm{i}\mathrm{n}}=720$, respectively (Table 2), the STD and RMS errors of which are 2.4 cm and 2.4 cm, respectively. The distribution of the difference is relatively uniform, which indicates the effectiveness of the residual and a-prior accuracy comparative analysis method and the treatment of the minimum degree of expansion of SRBFs.

6. Summary and Outlook

After the earlier work undertaken by Schmidt et al., Eicker, and Klees et al., SRBFs have been widely used for regional gravity field modeling in the last two decades. In this study, we combine terrestrial gravity data, airborne gravity data, GGMs and topographic models to calculate a high-resolution geoid model in Colorado based on band-limited SRBFs. Detailed explanations are given regarding the determination of the optimal expansion degree of gravity data: the optimal maximum expansion degree of terrestrial and airborne gravity data and the minimum expansion degree of gravity data. The main conclusions are as follows.

(1) The residual and a-priori accuracy comparative analysis method has proven to be effective in determining the expansion degree of gravity data. Using this method, we were able to determine the optimal expansion degree of terrestrial and airborne data to be 5200 and 1840, respectively. Additionally, the accuracy of our optimal geoid model, ColSRBF2023, was found to be 2.3 cm, which is the smallest compared to the other geoid models in this study. This implies that the expansion degrees of gravity data have influence on the modeling results, and cannot be overlooked when aiming for geoid modeling at the centimeter-level accuracy. Furthermore, the residual and a-priori accuracy comparative analysis method can also be utilized to detect the spectral information of gravity signals. This means that the method can provides valuable insights into the nature and characteristics of the gravity signals, enhancing our understanding of the geoid modeling process.

(2) The minimum degree of expansion of the gravity data also plays a role in the modeling result. If the expansion degree is not properly determined, it can lead to a geoid height error of ~5 mm. This highlights the importance of accurately determining the minimum expansion degree to achieve the desired level of accuracy in geoid modeling.

(3) By combining gravity data with different spectra using band-limited SRBFs, favorable modeling results can be achieved. Compared with ColSRBF2019, ColUNBSH2019, ColRLSC2019 and ColWLSC2020 submitted to ISG by different institutions, the ColSRBF2023 model showed a reduction in STD value of approximately 0.2-1.6 cm.

(4) Comparisons are also made with the mean results of all the contributions. Our geoid grid model for the whole study area delivers an STD value of 2.4 cm, which is also a small value for all participants. Moreover, there is favorable agreement between the STD error with regards to the area mean solutions and the STD error w.r.t. the GSVS17 GPS/leveling data, which shows the reliability of our geoid model ColSRBF2023.

All in all, in this study, we aimed to enhance the accuracy of the geoid model in Colorado using band-limited SRBFs. To achieve this, an investigation and determined the optimal expansion degrees for both terrestrial and airborne gravity data were conducted. The highest accuracy geoid model was obtained when the maximum expansion degrees for terrestrial and airborne gravity data were set to 5600 and 1840, respectively, while the minimum expansion degrees for both data types were set to 720.The study demonstrated the significance of accurately determining the spectral information of gravity data, particularly in regions with high topography, to improve the accuracy of geoid models. In future studies, we plan to perform further analyses of the remaining error sources in the determination of geoid models based on band-limited SRBFs, including the effect of the cutoff degree of the GGM, topographic reductions, and different approaches for handling airborne data, thereby furthering our contribution to the 1-cm geoid experiment of Colorado.