Hybrid Geoid Modeling with AI Enhancements: A Case Study for Almaty, Kazakhstan

1. Introduction

The creation of a precise geoid model is one of the most important tasks of modern geodesy, as it forms the foundation for establishing a unified national height system. This is particularly relevant for countries with vast territories and complex topography. The Republic of Kazakhstan, the ninth largest country in the world by area, is characterized by pronounced elevation contrasts, ranging from the Caspian Lowland to the mountain ranges of the Tien Shan and Altai. Under such conditions, the use of global geopotential models (GGMs), such as EGM2008 or XGM2019e, provides only an averaged description of the gravity field and does not allow achieving high accuracy level required for engineering and geodetic applications. Therefore, a tailored and customized model is required for such terrains. Conventionally, geodetic data obtained using leveling techniques and the global navigation satellite system (GNSS) is used for geoid modelling. However, they do have some disadvantages: Field surveying is expensive especially in remote and difficult-to access areas, and there are very few GNSS stations and it is time consuming to level up after each leveling operation. Also, data gaps might result from processes [1,2].

The development of high precision model is one of the key tasks within the framework of ongoing national program on establishment of state coordinate and height system. This will enable the replacement of labor-intensive and expensive geometric leveling with GNSS technologies, where the orthometric height is determined as the difference between the ellipsoidal height and the geoid height. Such a transition will ensure a significant reduction in costs while maintaining high accuracy and efficiency of geodetic works.

Over the last two decades, significant progress has been made in developing refined gravimetric geoid determination methods. The Least-Squares Modification of Stokes’ Formula with Additive corrections (LSMSA, also known as the KTH method) has become a widely accepted framework for regional geoid modelling [3,4,5,6]. This method rigorously integrates heterogeneous data-gravity anomalies, DEM-derived terrain corrections, and GGMs-while accounting for truncation and spectral errors. However, residual discrepancies between gravimetric solutions and GNSS/leveling observations often remain due to datum inconsistencies, measurement noise, or local heterogeneities [7]. Traditionally, such residuals have been minimized using low-parameter Helmert transformations; yet, these approaches may overfit and cannot adequately represent complex local biases.

Recent advances in machine learning (ML) have opened new opportunities for geodetic modelling, offering flexible, data-driven alternatives for residual correction. Methods such as Gaussian Process Regression (GPR), Support Vector Regression (SVR), and ensemble learners (e.g., LSBoost) have demonstrated the ability to capture nonlinear relationships and provide robust predictive performance across geospatial domains [6,8,9,10,11]. Their application to geoid modelling remains relatively new, but early studies suggest that ML-based corrector surfaces can significantly improve hybrid geoid accuracy compared to conventional techniques [12,13,14].

In this study we aim to develop a high-resolution hybrid geoid model for the Almaty region of southeastern Kazakhstan using LSMSA method and machine learning. We integrate digitized terrestrial gravity data, the global model XGM2019e_2159, the FABDEM 30 m DEM, and a network of GNSS/leveling benchmarks. The geoid is first computed using the LSMSA method to obtain a physically consistent gravimetric baseline. Subsequently, residual corrections are modeled using both classical Helmert fitting and advanced machine-learning regressors (GPR, SVR, LSBoost), enabling a comparative assessment of accuracy, robustness, and generalization. The final outcome is a refined hybrid geoid surface (N_ALM2025), designed to support precise geodetic, engineering, and scientific applications in one of the most challenging regions of Kazakhstan.

2. Study Area

The study area lies between latitudes 42.5°- 44° N and longitudes 76°-77.5° E, covering the southeastern part of Kazakhstan, specifically the Almaty region (Figure 1). This area includes the Zailiyskiy Alatau Range, a prominent subrange of the northern Tien Shan Mountains. The region is characterized by highly variable topography, with elevations ranging from approximately 600 meters in the lowlands to over 5000 meters in the high mountain peaks. Such dramatic elevation differences within a relatively small spatial extent create a complex geophysical environment.

Geologically, the region is seismically active and dynamic due to ongoing tectonic convergence between the Eurasian and Indian plates. These tectonic activities result in pronounced crustal deformations and localized gravity anomalies. The combined effect of rugged terrain and active tectonics introduces significant challenges for geoid determination, especially when relying solely on global geopotential models (GGMs). The spatial variation of gravity and topographic masses in this region often exceeds the resolution and accuracy capabilities of global models, necessitating the application of regional hybrid geoid modelling techniques.

Furthermore, the Almaty region holds strategic importance for geodetic and surveying activities due to its urban expansion, infrastructure development, and the presence of numerous scientific institutions. Accurate geoid models are critical for precise height systems and geospatial applications in such a topographically and geophysically complex region. Therefore, this area serves as an ideal testbed for evaluating the performance of hybrid geoid modelling methods that integrate satellite-based models, terrestrial gravity data, and digital elevation models (DEMs).

3. Information Base, Required for Local Geoid Modelling

The hybrid geoid model in this study is constructed using four complementary datasets: terrestrial gravity measurements, a global geopotential model (GGM), a high-resolution digital elevation model (DEM), and GNSS/levelling points. In present study, the long-wavelength component of the gravity field is represented by the global gravity model XGM2019e_2159, since it is one of the most up-to-date combined models integrating satellite (GRACE, GOCE) and terrestrial gravity data[15]. Topographic information and terrain corrections are derived from the FABDEM 30 m dataset, offering improved representation of bare-earth elevations [16]. The primary sources of information, however, are the gravity data—compiled from digitized historical surveys, validated by modern ground measurements, and harmonized with global models—and the GNSS/levelling data, which serve both for AI-supported geoid correction and independent validation of the geoid.

3.1. Gravity Data

For this study, gravimetric data are sourced from Soviet-era gravity anomaly maps at a scale of 1:200,000. These maps provide extensive coverage and detailed gravity information across Kazakhstan, with over 90% of the national territory surveyed using gravimeters such as SN-3, GAK, and GR-2K. In priority regions, higher-resolution surveys (1:50,000 or finer) were conducted [15]. The maps are based on Bouguer anomaly reductions using densities of 2.30 g/cm³ and 2.67 g/cm³ and include associated metadata such as free-air anomalies, elevation values, and terrain corrections. Many of the 1:200,000-scale maps incorporated data from earlier surveys conducted at scales ranging from 1:25,000 to 1:50,000.

Historical gravity database was created by digitizing 14,061 gravity observation points from the following nomenclature sheets: L-43-XXVIII, L-43-XXIX, L-43-XXX, L-44-XXV, L-43-XXXIV, L-43-XXXV, L-44-XXXI, K-43-IV, K-43-V, K-43-VI, K-44-I, K-43-X, K-43-XI, K-43-XII, and K-44-VII. (Figure 2).

The digitization process includes scanning of original paper maps at a resolution of 600 dpi to preserve fine details, such as gravity isolines and control points. Subsequently, each scanned sheet was georeferenced to the WGS84 coordinate system using at least four control points, typically at coordinate grid intersections. The georeferencing was initially performed in the Pulkovo 1942 datum within the relevant Gauss–Krüger projection zone. Polynomial transformations (second or third order) were applied where necessary, and positional accuracy was maintained within 0.1-0.3 mm at the original map scale.

After georeferencing, vectorization of gravity anomaly contours and observation points was performed, and relevant attributes—such as Bouguer anomalies, terrain corrections, and elevations—were entered into structured attribute tables. Priority was given to maps with Bouguer reductions using 2.67 g/cm³ with terrain correction. Maps with 2.30 g/cm³ density reductions were used only in areas lacking higher-quality datasets.

A quality control procedure was implemented by generating a digital gravity anomaly surface using Kriging interpolation (250 × 250 m grid) and comparing it to the original raster images. Any discrepancies, often due to digitization errors, were corrected, and affected points were flagged in the database. Final data transformation from Gauss–Krüger coordinates to WGS84 was conducted using the transformation parameters Pulkovo_1942_To_WGS_1984_16. All digitized sheets were then merged into a unified spatial dataset.

To address data gaps in remote or mountainous areas lacking terrestrial surveys, global gravity data from the WGM2012 model (5′ × 5′ resolution) were integrated (Figure 3). These supplementary data enhanced the spatial continuity of the gravity field and were essential for accurate interpolation in regions where local measurements were unavailable.

To assess the reliability of the gravity data, a comparison was conducted with modern absolute gravity measurements from QazGRF (Kazakhstan Gravimetric Reference Frame) stations [17,18,19] (Figure 4). Bouguer anomaly values were interpolated from the gravity data and corrected to estimate absolute gravity, which was then compared against observed values.

The statistical evaluation shows a mean difference of 3.23 mGal and a standard deviation of 2.95 mGal, with over 90% of discrepancies within ±5 mGal. This indicates the absence of systematic bias and confirms the generally high quality of the gravity data. The dataset was subsequently converted into a free-air gravity anomaly grid, which served as the primary input for the LSMSA method in local geoid determination.

3.2. GNSS/Leveling Data

A hybrid geoid model was derived and validated from a combined dataset of GNSS observations and second- to third-order leveling benchmarks, integrated with first-order triangulation points (Figure 5).

The network covers both lowland and mountainous areas in southeastern Kazakhstan, which allows for accounting the diversity of topographic and geophysical conditions essential for precise geoid modelling. Normal heights, defined in the BHS-77 system, range from approximately 626 m to 1466 m and were compared with ellipsoidal heights derived from GNSS measurements, providing a basis for quantifying and minimizing systematic differences between the physical and geometric height systems. The GNSS surveys were conducted in static mode using dual-frequency GPS/GLONASS receivers, with observation sessions lasting no less than 24 hours to suppress atmospheric effects, including iono-spheric and tropospheric delays. The raw data were collected in RINEX format and processed with the GAMIT/GLOBK software package. The GAMIT module applied precise IGS orbits, atmospheric delay models, and multipath corrections based on double-difference carrier phase solutions, while GLOBK was used for network adjustment and integration into the global IGS reference frame through nearby permanent stations (BADG, LHAZ, URUM, NOVM). The full GNSS post-processing cycle is summarized in four main stages: (i) synchronization of data to a single time scale, (ii) elimination of atmospheric distortions, (iii) construction of baseline vectors between observation points, and (iv) least squares adjustment of coordinates to minimize random errors (Figure 6).

The quality of the adjusted solution was assessed using weighted root mean square (WRMS) values for the east, north, and vertical components, which in most cases were below 0.050 m, with the most stable stations, such as KOTU and KURS, achieving WRMS values as low as 0.004 m. Although some sites initially exhibited higher RMS values due to multipath effects, the final adjusted coordinates remained well within accepted geodetic precision.

4. Methods of the Local Geoid Modelling

The proposed methodology integrates the classical Least-Squares Modification of Stokes’ Formula with Additive corrections (LSMSA, also known as the KTH method) and advanced machine learning techniques to improve the accuracy of local geoid determination. The workflow is divided into three main stages (Figure 7).

The first stage of our workflow follows the classical LSMSA method [20] This approach provides a rigorous theoretical framework to integrate heterogeneous datasets — terrestrial gravity anomalies, global geopotential models (GGMs), digital elevation models, and GNSS/leveling observations — into a single regional geoid model.

The LSMSA method optimizes modification parameters of the Stokes kernel in a least-squares sense, thereby minimizing the mean square error of the derived geoid heights[17]. It accounts for:

-

truncation and spectral errors of gravity anomalies;

-

regularization of the poorly conditioned system for high-degree harmonics;

-

additive corrections for topography, atmosphere, ellipsoidal effects, and downward continuation.

The output is a gravimetric geoid model $N_{LSMSA}$, which serves as the baseline solution for further refinement.

In the second stage, GNSS/leveling observations are preprocessed to geoid undulations compatible with the gravimetric baseline. Kazakhstan’s national height system is realized in normal heights—we first form height anomalies ζ and convert them to geoid undulations via the Moritz [21] geoid–quasigeoid separation formula, ensuring strict commensurability with the gravimetric geoid $N_{LSMSA}$ and an unbiased assessment of GNSS/Leveling results.

$$N_{\mathrm{G}\mathrm{N}\mathrm{S}\mathrm{S}/\mathrm{L}\mathrm{e}\mathrm{v}}=\zeta +{\delta N}_{Moritz}$$

(1)

The resulting $N_{\mathrm{G}\mathrm{N}\mathrm{S}\mathrm{S}/\mathrm{L}\mathrm{e}\mathrm{v}}$ values are then split into training (70%) and test (30%) subsets for subsequent Helmert and AI comparisons.

The third stage consists of building a data-driven corrector surface from the pointwise residuals between GNSS/leveling undulations and the gravimetric baseline, and then superimposing this field on the LSMSA surface to obtain the final hybrid geoid. Traditionally, this is achieved with a 3-, 5- ,7-parameters Helmert transformation, which aligns the two datasets by modelling shifts, rotations, and scale distortions. However, Helmert fitting can introduce systematic errors when datasets use different vertical datums or contain topographic inconsistencies. To improve the flexibility and accuracy of residual modelling, we propose the use of AI-based regression methods to approximate the difference between GNSS/leveling values and the gravimetric LSMSA model.

Specifically, let$\chi =\left[\phi ,\hspace{0.17em}\lambda ,{\hspace{0.17em}N}_{LSMSA} \right]$denote the feature vector composed of standard geodetic latitude and longitude (in decimal degrees) and the gravimetric baseline undulation interpolated at the point location. For each control point, we form the pointwise residual

$$\mathit{y}\left(\mathit{\chi }\right)=\mathit{\Delta }\mathit{N}\left(\mathit{\chi }\right)={\mathit{N}}_{\mathit{G}\mathit{N}\mathit{S}\mathit{S}/\mathit{L}\mathit{e}\mathit{v}}\left(\mathit{\chi }\right)-{\mathit{N}}_{\mathit{L}\mathit{S}\mathit{M}\mathit{S}\mathit{A}}\left(\mathit{\chi }\right)$$

(2)

and learn a data-driven regressor$\widehat{f}\left(\chi \right)$ such that $\widehat{f}\left(\chi \right)\approx \mathsf{\Delta }N\left(\chi \right)$. The final hybrid geoid is then obtained by superimposition,

$$N_{hyb}\left(\chi \right)=N_{LSMSA}\left(\chi \right)+f\left(\chi \right).$$

(3)

To approximate $\mathsf{\Delta }N$ over the study area while preserving smooth, geophysically plausible behavior, we evaluate three complementary regression families: Gaussian Process Regression (GPR), Random Forest (RF), and Support Vector Regression (SVR). These methods are trained on 70% of the GNSS/leveling dataset and validated on the remaining 30%, ensuring reliable generalization.

GPR was implemented using an Automatic Relevance Determination Squared Exponential (ARD) kernel (GPR, ARD kernel). We place a Gaussian process prior on the latent field $f(\cdot )$ with i.i.d. Gaussian observation noise,

$$f\left(\chi \right)\sim \mathcal{G}\mathcal{P}\left(0, k_{\theta }\left(\chi ,{\chi }^{\text{'}}\right)\right), y_i=f\left({\chi }_i\right)+{\epsilon }_i, {\epsilon }_i\sim \mathcal{N}\left(0, {\sigma }_n^2\right). \left(4\right)$$

(4)

An automatic relevance determination (ARD) kernel is used so that the model learns separate length-scales for φ, λ, and N_LSMSA. In practice we adopt a Matérn 3/2 covariance (with fallbacks to Matérn 5/2 or squared-exponential when required). Writing

$$r\left(\chi ,{\chi }^{\text{'}}\right)=\sqrt{\sum _j{\displaystyle \frac{{(x_j-{x\text{'}}_j)}^2}{l_j^2}}}$$

(5)

The kernel read:

$$k_{Matèrn-3/2}={\sigma }_{f }^2\left(1+\sqrt{3}r\right)e^{-\sqrt{3}r}, k_{SE}= {\sigma }_f^2\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{1}{2}}\sum _j{\displaystyle \frac{{(x_j-x_j^{\text{'}})}^2}{l_j^2}},)$$

(6)

Given training inputs X and targets y, the posterior predictor at x* has mean and variance

$${\mu }_{*}\left(x_{*}\right)=k_{*}^T{\left(K+{\sigma }_n^{2}I\right)}^{-1}y,{\sigma }_{*}^2\left(x_{*}\right)=k_{**}{-k_{*}^T\left(K+{\sigma }_n^{2}I\right)}^{-1}{k}_{*},$$

(7)

with hyperparameters $\theta =\left\{{\sigma }_f\right.,\left\{l_j\right\},\left.{\sigma }_n\right\}$ fitted by maximizing the log marginal likelihood. Using φ and λ as standard degrees is handled by ARD and z-standardization, which absorb scale differences between the angular and metric features.

Support Vector Regression (SVR, RBF kernel)

SVR seeks a function with small RKHS norm while keeping residuals within an ε -insensitive tube:

$$\begin{array}{c}min\\ w,b,\mathsf{\xi },{\mathsf{\xi }}^{*}\end{array}{\displaystyle \frac{1}{2}}{‖w‖}^2+C\sum _{i=1}^n({\mathsf{\xi }}_i+{\mathsf{\xi }}_i^{*}) s.t. \left\{\begin{array}{c}{y}_i-f\left({\chi }_i\right)\le \epsilon +{\mathsf{\xi }}_i,\\ f\left({\chi }_i\right)-y_i\le {\epsilon }_i^{*},\\ {\mathsf{\xi }}_i,{\mathsf{\xi }}_i^{*}\ge 0.\end{array}\right.$$

(8)

By the representer theorem,

$$f\left(x\right)=\sum _{i=1}^n\left(a_i-a_i^{*}\right)k\left({\chi }_i,\chi \right)+b,$$

(9)

with an RBF kernel $k\left(\chi ,{\chi }^{\mathrm{\text{'}}}\right)=\mathrm{e}\mathrm{x}\mathrm{p}(-{‖\chi -{\chi }^{\mathrm{\text{'}}}‖}^2/(2^2\left)\right)$. The hyperparameters C, ε, and σ (KernelScale) govern smoothness and robustness. As with GPR, inputs (φ,λ,NLSMSA) are kept in their standard units and standardized prior to training.

Least-Squares Boosting (LSBoost)

We fit an additive ensemble of regression trees by steepest descent on the squared loss. Starting from

$$f_0\left(\chi \right)=\begin{array}{c}arg \mathrm{m}\mathrm{i}\mathrm{n}\\ c\end{array}\sum _i(y_i-c{)}^2=\overline{y}$$

(10)

we iterate for m=1,…,M:

$$r_{im}=y_i-f_{m-1}\left({\chi }_i\right), h_m\leftarrow tree fit to\left\{\left({\chi }_i\right)\}, r_{im})\},\right. f_m\left(\chi \right)=f_{m-1}\left(\chi \right)+h_m\left(\chi \right),$$

(11)

with learning rate $\in \left(\mathrm{0,1}\right]$. Regularization is controlled by (M,ν) and tree complexity (e.g., minimum leaf size). This yields a flexible, piecewise-smooth approximation of ΔN(φ,λ,N_LSMSA) that complements the smooth, kernel-based GPR/SVR models.

To assess the effectiveness of the proposed approach, we perform a comparative analysis of Helmert 7-parameter fitting versus AI-based regression fitting (GPR, RF, SVR).

The comparison is carried out on the 30% evaluation dataset using RMSE, MAE, correlation coefficients, and residual distribution statistics. The model achieving the best generalization performance and stability was selected as the final regression method for constructing the hybrid geoid model $N_{ALM2025}$ for the Almaty region.

5. Results and Discussion

5.1. Classical Geoid Modelling by LSMSA Method

Geoid model calculations included the use of XGM2019 spherical harmonic coefficients to powers of 120, 200, 300, 400, 500, 630, 760, with a combination of ground gravity data error variance C (0) - 0.5, 1, 3, 6, 9, 16 mGal². The results are presented in Table 1, Figure 8.

The initial parameters for calculating the modification parameters were:

• Degree of modification L=M= 300;
• Variance of errors in ground gravity data C (0) =1 m Gal ²;
• The size of the integration coverage is Ѱ= 1^o.

In the computational scheme of the LSMSA method for determining the geoid, the surface gravity anomalies and the GGM are used to determine Napp (Figure 9)., and then necessary corrections are added (Figure 10, Table 2).

Figure 11 shows the result of calculating the geoid heights for the study area using the method “Least Squares Modification of Stokes formula with additional corrections (LSMSA)”.

The calculated geoid heights for Almaty region using the LSMSA method vary in the range from -46.485 m to -37.644 m, the average value is -42.577 m, and the standard deviation is 2.873 m.

5.2. Corrector Surface by AI

We evaluated residuals between observed geoid undulations $N_{GNSS/Lev}$ and the ${N }_{LSMSA}$ grid interpolated to 40 control points (train = 28, test = 12). Table 3 summarizes distributional statistics for the baseline (“Before”) and for four fitted models (“After”).

Table 4 presents the residual error statistics of geoid correction for the entire dataset (ALL), the training subset (TRAIN), and the validation subset (TEST).

In the baseline solution, the residuals exhibit a systematic negative bias (Mean ≈ –0.08 m) with an RMS of 0.097 m, indicating the presence of unmodeled systematic errors. Both training and test subsets show similar bias levels, confirming that the discrepancy is not random but systematic.

Application of the Helmert transformation reduces the training errors significantly (RMS = 0.031 m); however, the test subset deteriorates markedly (RMS = 0.248 m), revealing strong overfitting and poor generalization. Although effective in aligning datasets under controlled conditions, the method fails to provide stable corrections across independent data.

In contrast, the Gaussian Process Regression (GPR, ARD-SE kernel) achieves nearly unbiased residuals (Mean = –0.003 m) with the lowest RMS over the full dataset (0.025 m). The test subset maintains relatively high accuracy (RMS = 0.045 m), demonstrating superior generalization. The Support Vector Regression (SVR, RBF kernel) also provides stable results with low RMS values (0.039 m overall; 0.032 m on test), though slightly less accurate than GPR.

The LSBoost ensemble method yields the smallest overall RMS (0.022 m) and near-zero bias. While the training errors approach zero (RMS = 0.005 m), the test subset retains acceptable accuracy (RMS = 0.040 m), similar to GPR.

5.3. Hybrid Geoid (LSMSA + AI)

Given its superior held-out performance and balanced residual distribution, SVR (RBF) was selected for the final hybrid grid

$$N_{ALM2025} ={ N}_{LSMSA} + f_{SVR}\left(∆X\right),$$

(8)

The resulting grid, computed at 2′ × 2′ resolution, preserves the long-wavelength physical consistency of LSMSA while effectively removing local biases.

Residual statistics confirm a marked improvement that, the initial RMS of ~0.10 m was reduced to ~0.03–0.04 m, and the systematic offset (≈−0.08 m) was eliminated. Figure 12 illustrates the shift from a negatively biased distribution (before correction) to a centered and compact one after SVR-based adjustment.

6. Conclusions

This study successfully developed and evaluated a high-resolution hybrid geoid model for the Almaty region of southeastern Kazakhstan, one of the most topographically and geophysically complex areas in Central Asia. The work integrated digitized Soviet-era gravity data, the global geopotential model XGM2019e_2159, high-resolution DEMs, and GNSS/leveling benchmarks within the framework of the Least-Squares Modification of Stokes’ Formula with Additive corrections (LSMSA).

The baseline gravimetric solution demonstrated the ability of the LSMSA method to capture the long- and medium-wavelength features of the gravity field; however, systematic residuals of approximately –0.08 m remained when compared against GNSS/leveling observations. Classical Helmert transformation reduced these discrepancies in the training dataset but introduced instability and overfitting in the independent validation subset, confirming its limited suitability for robust hybrid modelling.

By contrast, machine-learning approaches provided substantial improvements. Gaussian Process Regression (GPR) and LSBoost yielded nearly unbiased residuals with RMS values in the range of 0.022–0.045 m, while Support Vector Regression (SVR, RBF kernel) offered a balance of accuracy and stability across training and test subsets. The final hybrid model, constructed using LSMSA combined with AI-based regression correction, achieved centimetre-level consistency with GNSS/leveling data and eliminated the systematic bias present in the baseline solution.

These results highlight two key conclusions. First, high-quality hybrid geoid models can be developed for regions with challenging topography by combining legacy gravimetric surveys with modern global and satellite-based datasets. Second, machine-learning regressors outperform classical parametric adjustments in modelling local residuals, offering a flexible and generalizable framework for future geoid determination efforts in Kazakhstan and comparable mountainous environments.

Future work will focus on expanding the methodology to a national scale, integrating additional GNSS/leveling networks.