Enhancing Regional Quasi-Geoid Refinement Precision: An Analytical Approach Employing ADS80 Tri-Linear Array Stereoscopic Imagery and GNSS Gravity-Potential Leveling

Abstract

This research investigates precision enhancement in regional quasi-geoid refinement through ADS80 tri-linear array scanning stereoscopic imagery for aerial triangulation coupled with GNSS gravity-potential modeling. By acquiring stereoscopic imagery and analyzing triangulation accuracy using an ADS80 camera, we performed this study over the Qinghai–Tibet Plateau’s elevated, desolate terrain, collecting 593 GNSS points following high-precision stereoscopic imagery modeling. By utilizing 12 gravity satellite models, we computed geoid heights and China’s 1985 Yellow Sea elevations for 28 benchmarks and GNSS points, thereby refining the Qinghai Province Quasi-Geoid Model (QPQM) using geometric techniques. The findings reveal that POS-assisted ADS80 stereoscopic imagery yields high-precision triangulation with maximal horizontal and elevation accuracies of 0.083/0.116 cm and 0.053/0.09 cm, respectively, across five control point arrangements. The RMSE of normal heights for 1985, processed via these GNSS points, achieved decimeter precision. By applying error corrections from benchmarks to the 1985 elevation data from gravity satellites and performing weighted averaging, the precision of EGM2008, SGG-UGM-2, and SGG-UGM-1 models improved to 8.61 cm, 9.09 cm, and 9.38 cm, respectively, surpassing the QPQM by 9.22 cm to 9.99 cm. This research demonstrates that the proposed methods can significantly enhance the precision of regional quasi-geoid surfaces. Additionally, these methods offer a novel approach for rapidly establishing regional quasi-geoid models in the uninhabited areas of the Qinghai–Tibet Plateau.

1. Introduction

The geoid is a gravitational equipotential surface that closely approximates the mean sea level and accurately reflects the Earth’s internal density distribution. It serves as the reference surface for defining orthometric height systems in geodesy. The quasi-geoid is derived from the geoid and utilizes the telluroid as the boundary surface for Molodensky’s boundary value problem [1], substituting the quasi-geoid in place of the geoid in solutions that historically applied the Stokes boundary value problem [2,3]. This replacement facilitates more practical and convenient methods for determining the Earth’s shape.

While the quasi-geoid is not strictly a level surface and an equipotential surface of gravity, it closely approximates one and is employed as an auxiliary computational surface; the normal height system of 1985, which is strictly physically meaningful, has been adopted in China [4]. With the continuous advancement of modern surveying technologies, there is an increasing urgency for high-precision regional quasi-geoid models in the unpopulated regions of the Tibetan Plateau to support analyses in glaciers, permafrost, and natural ecosystems, among other fields. The traditional methods used in the construction of quasi-geoid models, such as leveling measurements [4], GNSS (Global Navigation Satellite Systems) measurements [5], and terrestrial gravity surveys, face numerous challenges in these uninhabited areas of the Tibetan Plateau. In recent years, the integration of GNSS geopotential leveling with the maturation of optical satellite remote sensing, medium- and low-altitude aerial photogrammetry, and satellite gravimetry technologies [6,7] has offered new technical pathways for enhancing the accuracy of regional quasi-geoid models in the unpopulated areas of the Tibetan Plateau.

The Leica Geosystem Airborne Digital Sensor ADS80, as cited by [8,9], serves as the data platform for acquiring three-dimensional stereo images in the remote regions of the Qinghai–Tibet Plateau in this paper. By utilizing the three-dimensional stereo imagery to produce DOMs and DEMs, as well as to collect spatial coordinates [10], the GNSS gravity-potential leveling method enables the direct acquisition of normal height coordinates based on the regional quasi-geoid model. This approach significantly reduces the workload and costs associated with traditional leveling measurements. Wei [7] proposed the GPS gravity-potential leveling method in 2007. By leveraging the complementary aspects of geometric and physical geodesy, this method achieves an integrated approach that combines the geodetic datum, which is based on the reference ellipsoid, with the vertical datum, grounded in mean sea level. This approach facilitates the integration of three-dimensional geo-positioning with altitude elevation measurements, enabling the provision of four-dimensional coordinates: latitude, longitude, geodetic height, and altitude. Employing GPS leveling data and the Eigen-cg03c Earth gravity field model enhances the precision of altitude measurements to 0.5 m and achieves better than 10 cm precision in altitude differences over distances of tens of kilometers. Previously, the elevation accuracy obtained through this method was primarily limited by the precision of GPS-measured geodetic heights within the WGS-84 system [11,12].

The geodetic heights obtained from GNSS measurements can accurately meet the high-precision positioning requirements in conventional surveying and remote sensing fields. Furthermore, several countries have sequentially released high-order global gravity models, such as gravity satellite SGG-UGM-2 [13], XGM2019e_2159 [14], GECO [15], EIGEN-6C4 [16], and EGM2008 [17]. The calculations derived from high-order gravity satellite models have achieved mgal-level precision in measuring gravity values and gravitational potentials. Moreover, Interferometric Synthetic Aperture Radar (InSAR) technology can utilize the phase of radar single-look complex image data to obtain deformation signals in the study area, providing reliable theoretical support for the stability analysis of the region [2]. These advancements offer substantial theoretical support for the sophisticated integration of aerial photogrammetry and the GNSS gravity-potential leveling method in regional quasi-geoid modeling, as presented in this paper.

Therefore, this study tackles the challenges associated with conducting GNSS, gravity, and leveling measurements in the uninhabited regions of the Qinghai–Tibet Plateau. Considering the necessity of establishing a regional quasi-geoid model in this area, we emphasize the feasibility of integrating aerial photogrammetry with a gravity-potential leveling method in data processing (Figure 1). The stability of the study area was analyzed using a Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) [2], and 12 different gravity satellites were selected for calculating the gravity-potential $W_P$ and $W_O$ with respect to the study area and the reference surface, as well as the average normal gravity ${\gamma }_m^P$ and $g_m^P$ (Equations (1) and (2)). The study employed a King Air 350 aircraft (Figure 1) equipped with the Leica ADS80 push-broom digital photogrammetry system to capture stereo imagery data of the study area. The ADS80 stereo imagery data underwent a high-precision POS-assisted aerial triangulation using a densification ray method for regional network adjustment, thereby accurately restoring the true three-dimensional surface features at the time of aerial photography. Concurrently, the three-dimensional coordinates of surface feature points were extracted in the CGCS2000 coordinate system (China Geodetic Coordinate System 2000). Based on these three-dimensional coordinates, the GNSS gravity-potential leveling method [7,11,18] was employed to calculate the normal height, $H_P^{\gamma }$, the height anomaly, $\xi$, and the difference between the geoid and the reference ellipsoid, $N$ (Figure 1, Section 2.5).

This study established conversion relationships among various elevation systems. These systems include the reference ellipsoid, geoid, quasi-geoid, and Earth’s surface normal height in the China 1985 Yellow Sea elevation system, $H_W$ [19,20]. For the China 1985 height data, it is the geo-potential value of the mean sea level at Qingdao tidal gauge station. The orthometric height and the normal height of the mean sea level at Qingdao tidal gauge station were both defined as zero in China, as well as the orthometric height, $H_P^g$ (WGS-84), and the geodetic height, $H_P$ (CGCS2000), as delineated in Figure 1 and following the methodologies of Li et al. [19] and Yang [21]. This approach effectively overcomes the limitations inherent in traditional leveling measurements within China. It aims to serve as a near-equivalent alternative to leveling measurements in uninhabited areas, thereby facilitating the rapid establishment or refinement of high-precision regional quasi-geoid models in these challenging regions.

2. Materials and Methods

2.1. Study Area

The study area is situated in a high-altitude uninhabited region at the border between Qinghai and Tibet Provinces, as depicted in Figure 2. The region’s average elevation is approximately 3000 m, with several significant faults, including F212, F217, F214, F216, and F221, distributed in the surrounding areas [22]. Within the study area, there are neither single nor multiple faults that span across its entirety. We compiled a dataset of seismic events from 1902 to 2023 in this region [23], which revealed only two earthquakes with magnitudes below 5.0 MW occurring within the study area. For the quantitative analysis of the vertical and horizontal displacement variations in the region, we utilized the European Space Agency’s (ESA) Sentinel-1A imagery, along with imagery from the Chinese satellites GF-2 [23] and GF-7 [24].

These analyses were conducted using InSAR [25,26] and optical image pixel offset methods [27,28]. The terrain in the study area is observed to be relatively gentle, primarily consisting of plains with an absence of large mountainous areas or canyons. The geological plate and fault activities in this region are also relatively stable. These geographical and geological conditions are conducive to enhancing the precision of the regional quasi-geoid model refinement, which is accomplished using the GNSS gravity-potential leveling model, as detailed in the subsequent sections.

2.2. Acquisition and Processing of GF2, GF7, and Sentinel-1A Imagery

In this study, high-resolution optical imagery was employed for the extraction of horizontal displacements in the study area, utilizing image matching and correlation analysis methods [27,28]. Additionally, 87 scenes of Sentinel-1A imagery were utilized to extract line-of-sight (LOS) displacements, employing the SBAS-InSAR technique [29,30]. The main image from 18 February 2019, along with all other images, was subjected to multi-view processing, registration, and interferogram coherence processing. Interferometric pairs exceeding the spatiotemporal baseline threshold were excluded, resulting in 425 interferometric pairs that met the spatiotemporal baseline threshold conditions (Figure 5c). The Minimum Cost Flow (MCF) [31] method was employed for spatiotemporal three-dimensional phase unwrapping. In order to address the issue of 2π ambiguity, the Goldstein method was used for the phase filtering of the interferograms, eliminating the estimated and residual constant phase and remaining phase ramps post-unwrapping. Ground Control Points (GCPs) were used to refine the orbit and re-leveling, and atmospheric high-pass filtering was applied to remove atmospheric phase effects. Corrections for Digital Elevation Model (DEM) errors, atmospheric delays, and orbital error interference phases were performed. The deformation rate and residual topography were estimated, and the corrected phase and unwrapped phase were combined and converted into elevation data. These data were then geocoded to the WGS84 coordinate system, and the Line-Of-Sight (LOS) deformation values were multiplied by ${\cos \theta }^{-1}$ [30] (where $\theta$ is the incidence angle) to convert them to vertical deformation values, obtaining the LOS deformation quantity (Figure 5a) and vertical deformation quantity (Figure 5b) in the overlapping area.

The high-resolution optical imagery was acquired through China’s High-Resolution Earth Observation System (abbreviated as “High-res Special Project”), using the Gaofen-2 (GF-2) and Gaofen-7 (GF-7) satellites. These satellites offer sub-meter spatial resolution, high temporal resolution, and efficient spectral resolution capabilities [32]. The spatial resolutions of GF-2 and GF-7 are 0.6 m and 0.7 m, respectively [33].

We selected satellite imagery data from the GF-2 and GF-7 satellites, covering the period from 2020 to 2022. This dataset included 32 scenes from GF-2 and 46 scenes from GF-7. The panchromatic and multispectral imagery obtained from these GF-2 and GF-7 optical satellites were processed using the Geomatica Banff satellite imagery data processing platform. The workflow entailed computing Rational Polynomial Coefficient (RPC) parameters, conducting aerial triangulation via the encrypted bundle adjustment method, and applying corrections to both panchromatic and multispectral band imagery. The subsequent stages included image fusion, enhancement, and mosaicking with cropping. This process led to the generation of multi-temporal true-color Digital Orthophoto Map (DOM) imagery data [34,35]. The DOM generated via the previously described methods was analyzed using Co-registration of Optically Sensed Images and Correlation (COSI-Corr) software (cosi-corr_pakOct19) [28]. This analysis focused on identifying changes in corresponding feature points, thereby facilitating the examination of horizontal surface deformations (Figure S1), which reflect lateral shifts in these feature positions [31,36,37].

2.3. Acquisition and Processing of ADS80 Push-Broom Tri-Linear Stereo Imagery Data

We utilized the KingAir350 aircraft (Figure 3a), equipped with the Leica ADS80 three-line scanner digital photogrammetry aerial survey system (Figure 3(a1–a3)) to conduct aerial surveys. This allowed us to acquire stereo imagery data with a resolution exceeding 0.2 m for the study area (Figure 3(a4)). The ADS80, a premium airborne digital aerial photogrammetry system developed by Leica of Switzerland, incorporates the SH92 lens and is complemented by the state-of-the-art gyro-stabilized platform PAV80. It also features an operational display unit OI40 and the CU80 camera control system. This system seamlessly integrates a high-precision Position and Orientation System (POS) derived from both an Inertial Measurement Unit (IMU) and the Global Navigation Satellite System (GNSS) [38,39].

In order to integrate the Leica ADS80 aerial survey system, modifications were made to the KingAir350 aircraft’s payload platform. The resulting configuration, as illustrated in Figure 3(a3), enabled the KingAir350 to operate in high-altitude plateau environments, with a flight ceiling ranging from 0 to 10,000 m. This ensured compliance with the altitude requirements specified in our flight plan. By using Mission-Pro flight planning software (Version number: Leica Mission Prov 12.11.0), a total of 16 east–west flight paths were designed for the study area (Figure 3b). The average ground projection spacing between images was 5 km, with an average lateral overlap of 35% between adjacent flight strips. Each flight strip was defined by three waypoints, with the initial and final waypoints serving as the starting and ending points for image capture (Figure 3b). Additionally, a pre-imaging waypoint was established 3 km before the starting point. Upon reaching this waypoint, the gyro-stabilized platform PAV80, the POS system (IMU + GNSS), and the camera control system CU80 initiated initialization and entered a pre-imaging state. The aerial surveyor used the OI40 monitor to engage in real-time communication with the pilot, discussing the current flight altitude, discrepancies from the designated flight altitude, and aircraft orientation. This facilitated adjustments to the flight path in preparation for a seamless transition into the image capture starting point.

During the aerial survey flights, the aircraft executed a “Figure-eight” flight pattern before reaching the initial waypoint and again after completing the final waypoint. This maneuver was employed to mitigate potential linear drift errors in the IMU resulting from extended flight durations. Additionally, each individual flight path was limited to approximately 30 min. The acquisition of stereo imagery data in the study area was contingent on weather conditions and airspace availability, necessitating two flight missions, as shown in Figure 3c. Prior to the flights, ground control points were established and marked for measurement. The study area adopted a grid-based layout for the distribution of 46 ground control points, as shown in Figure 3d. Given the challenging and high-altitude terrain of the measurement area, we utilized the Precise Point Positioning (PPP) method for ground control point measurements. The Gamit/Globk software (Version number: Gamit/Globk 10.71), in conjunction with International GNSS Service (IGS) reference stations, was employed for precise coordinate calculations, with all measurements referenced to the CGCS2000 coordinate system [21]. In order to achieve high-precision POS data, a post-processing differential mode was employed for POS data computation. In the study area, three GNSS reference stations were established, ensuring that the baselines of the GNSS reference stations were confined within a 50 km radius. This setup guaranteed that synchronized observations covered the entirety of the region, as shown in Figure 3e. During the flight of the Leica ADS80 aerial survey system, these stations were activated synchronously, working collaboratively with the onboard differential GNSS receiver. This synchronous observation aimed to counteract errors from factors such as ionospheric and tropospheric disturbances during the computation of onboard differential POS data.

Following the completion of the flight, data processing steps such as downloading, flight quality checks, and pre-processing, were carried out using the Xpro software (Version number: XPro 6.4.7) to obtain Level 1 (L1) data. Next, the initial exterior orientation elements of the aerial camera were calculated using Inertial Explorer software (Version number: Inertial Explorer 8.7). Finally, combining internally processed ground control points, L1 images, and exterior orientation elements, aerial triangulation was performed using Xpro software (Version number: XPro 6.4.7). This process involved applying the collimation equation for geometric correction, producing L1 image data optimized for stereo mapping [10]. After completing these steps, the Mapmatrix Digital Stereoplotting Photogrammetric Workstation [40] was utilized for the stereo editing and collection of three-dimensional coordinate points [41]. Subsequent processing involved generating Digital Elevation Models (DEMs) from orthophoto images, as well as creating DOM.

2.4. Data Extraction of GNSS Gravity-Potential Control Points in Dimensional Environments

In Section 3.3.1, the precision of aerial triangulation densification for control point scheme five, as shown in Figure 7g,h, achieved centimeter-level accuracy in both horizontal and vertical directions. By leveraging this precision, the Chinese digital photogrammetric mapping system MapMatrix3D was utilized to construct and restore a digital 3D stereoscopic model of the study area. With the support of a stereoscopic graphics card and wearing 3D glasses, we authentically recreated the stereoscopic scene of the study area, as shown in Figure 4a. Through a combination of different stereoscopic image display modes—namely Pattern 1 (Dual Stereoscopic), Pattern 2 (Anaglyphic), and Pattern 3 (Polarized Light)—and manual interpretation, a total of 594 points coordinates based on the CGCS2000 coordinate system were extracted (pertaining to GNSS gravity-potential control points) within the stereoscopic environment. In order to control the precision effectively during point extraction, the point positions in the study area were referenced to the decimeter-precision regional quasi-geoid model (Li 2012) of Qinghai Province (Qinghai Province Quasi-Geoid Model, QPQM), as depicted in Figure 2. Point coordinates were collected near the QPQM model grid points, with an average grid spacing of 2′30″ (Figure 4b).

The QPQM model was constructed using 213,862 terrestrial gravity data points within the Qinghai Province territory and 101 high-precision GPS leveling data points. The EGM96 was employed as the reference gravity field, with terrain corrections applied using Shuttle Radar Topography Mission (SRTM-3) data. The height of the gravitational quasi-geoid was interpolated and extrapolated via the remove-restore technique incorporating Molodensky’s method. The resultant quasi-geoid of Qinghai Province, with a resolution of 2′30″ and an accuracy of better than 0.186 m [42], was determined by removing the vertical deflections between GPS leveling and the gravity quasi-geoid [43]. This model also serves as the foundational grid digital model for refining the regional quasi-geoid using the GNSS gravimetric leveling method in subsequent sections, and it provides a benchmark for precision assessment.

2.5. Acquisition and Processing of GNSS/Leveling Data

In order to verify the accuracy of the normal heights calculated using the GNSS gravity-potential leveling model [42,43], this study selected two leveling lines traversing the study area in an east–west direction. These lines were connected to the national second-order leveling lines. GNSS and leveling measurements were conducted at six and 22 specific points along these two routes, respectively, as illustrated in Figure 9. Orthometric height measurements (referenced to the 1985 orthometric height datum) were conducted at the same leveling points [43,44], adhering to the standards of the national second-class leveling route [44]. Concurrently, static GNSS measurements were performed to obtain the CGCS2000 geodetic height at these points [45].

For the computation of orthometric heights at the leveling points, all points were initially constrained to the national benchmarks, which were in line with the specifications of the traverse leveling route. Following this, error distribution was executed utilizing the least squares method, based on the principles of indirect adjustment under conditional constraints. Static GNSS measurements at each leveling point were processed using the Gamit/Globk software (Version number: Gamit/Globk 10.71) [46,47]. This processing incorporated precise ephemerides, clock corrections, broadcast ephemerides, and the accurate coordinates and velocity fields of IGS stations, as published by the official National Aeronautics and Space Administration (NASA), for baseline solution and adjustment [48,49]. The processed coordinates were then transformed to correspond with the 2000th epoch within the ITRF97 framework, resulting in ellipsoidal heights based on the CGCS2000 geodetic height data.

2.6. Acquisition and Processing of Gravity Satellite Data

We integrated mean gravity and gravity anomalies, utilizing the GNSS-derived geoid model to compute the orthometric heights for all geodetic height points acquired as described in Section 2.3 and Section 2.4 We selected 12 global gravitational field models spanning degrees/orders from 360 to 2196 [50]. These models include EGM2008, SGG-UGM-1, SGG-UGM-2, Eigen-6C4, XGM2019e-2159, EGM96 [51,52], GECO [15,53], GGM05C [17], GO-CONS-GCF-2-DIR-R6 [13], Tongji-GMMG2021S [54], Eigen-CG03C [16], and Tongji-GRACE02 [54]. All the satellite gravity data were sourced from the International Centre for Global Earth Models (ICGEMs). The geodetic coordinates of the test points, which were obtained by using the method described in Section 2.3 (CGCS2000 coordinate system), were utilized to calculate their gravitational potential based on the Earth’s gravitational model [55].

2.6.1. Height System

Assuming that the gravity-potential at a known height reference point (the Qingdao tide gauge zero-level datum, Li et al. [56]) is $W_0$ (62,636,854.2 ± 0.5 m²/s²) and that the gravity-potential of the geoid surface is $W_0$ (62,636,856.00 ± 0.5 m²/s²), which is obtained through leveling and gravity measurements for the point of interest $P$ (Figure 1), the relative difference in gravity-potential relative to the height reference point, denoted as $C=W_0-W_P$, is determined. For the normal height system [7,57], the elevation of point $P$ is given by

$$H_P^{\gamma }={\displaystyle \frac{C}{{\gamma }_m^P}}={\displaystyle \frac{W_0-W_P}{{\gamma }_m^P}}$$

(1)

In the equation, $H_P^{\gamma }$ represents the normal height of point $P$; $H_P^g$ represents the orthometric height of point $P$; $W_P$ stands for the gravity-potential at point $P$; ${\gamma }_m^P$ denotes the average normal gravity between point $P$ and the quasi-geoid surface along the plumb line; $g_m^P$ denotes the average gravity along the plumb line from point $P$ to the geoid surface.

$$H_P^g={\displaystyle \frac{C}{g_m^P}}={\displaystyle \frac{W_0-W_P}{g_m^P}}$$

(2)

2.6.2. Gravity-Potential

The gravitational potential, $W_P$, at a point $P\left(r,\phi ,\lambda \right)$ is defined as the sum of the Earth’s gravitational potential, $V_P$, and the centrifugal potential, ${\Phi }_P$, induced by the Earth’s rotation [7,58,59]. This can be expressed as follows:

$$W_P=V_P+{\mathsf{\Phi }}_P=\mathrm{V}\left(\mathrm{r},\phi ,\lambda \right)+\mathsf{\Phi }\left(\mathrm{r},\phi ,\lambda \right)$$

(3)

$$\mathrm{V}\left(\mathrm{r},\phi ,\lambda \right)={\displaystyle \frac{GM}{r}}\left[1+{\sum }_{n=1}^{N_{max}}{\left({\displaystyle \genfrac{}{}{0pt}{}{\underset{\_}{a}}{r}}\right)}^n\right]\sum _{m=0}^n\left[{\overline{C}}_{nm}\cos m\lambda +{\overline{S}}_{nm}\sin m\lambda \right]{\overline{P}}_{nm}\left(cos\phi \right)$$

(4)

$$\mathsf{\Phi }\left(\mathrm{r},\phi ,\lambda \right)={\displaystyle \frac{1}{2}}{\omega }^2\left(x^2+y^2\right)={\displaystyle \frac{1}{2}}{\omega }^{2}r{sin}^2\phi$$

(5)

${\overline{C}}_{nm}$ and ${\overline{S}}_{nm}$ are the fully normalized geopotential coefficients of degree n and order m; ${\overline{P}}_{nm}\left(sin\phi \right)$ represents the fully normalized associated Legendre functions of degree n and order m; GM is the gravitational constant of the Earth; a is the major semi-axis of the reference ellipsoid; $N_{max}$ denotes the maximum degree of the gravitational potential model. $\omega$ is the angular velocity of the Earth’s rotation; r, $\phi$, and $\lambda$ are the geocentric radius, geocentric latitude, and longitude, respectively. The calculation formula is

$$\left\{\begin{array}{c}r=\sqrt{x^2+y^2+z^2}\\ \phi ={\tan }^{-1}\left({\displaystyle \frac{1}{\sqrt{x^2+y^2}}}\right)\\ \lambda ={\tan }^{-1}\left({\displaystyle \frac{y}{x}}\right)={\cos }^{-1}\left({\displaystyle \frac{x}{\sqrt{x^2+y^2}}}\right)={\sin }^{-1}\left({\displaystyle \frac{y}{\sqrt{x^2+y^2}}}\right)\end{array}\right.$$

(6)

where x, y, and z represent the geocentric Cartesian coordinates.

2.6.3. Normal Gravity

For the normal height system, the average normal gravity from point $P$ to the quasi-geoid surface is given by [7,11,60]:

$${\gamma }_m^P={\gamma }_{p·}\left[1-\left(1+f+{\displaystyle \frac{{\omega }^{2}a}{{\gamma }_p}}-2f{\sin }^{2}B\right){\displaystyle \frac{H_P^{\gamma }}{a}}+{\displaystyle \frac{{H_P^{\gamma }}^2}{a^2}}\right]$$

(7)

$f$ denotes the ellipsoidal flattening; B represents the geodetic latitude of point $P$; ${\gamma }_p$ signifies the normal gravity at the projection point, $P$, on the ellipsoidal surface.

$${\gamma }_p={\displaystyle \frac{{\gamma }_e\left(1+k{\sin }^{2}B\right)}{\sqrt{1-e^2{\sin }^{2}B}}},k={\displaystyle \frac{b{\gamma }_a}{a{\gamma }_e}}-1$$

(8)

${\gamma }_e$ represents the normal gravity at the equator; ${\gamma }_a$ denotes the normal gravity at the pole; $e$ stands for the first eccentricity of the meridian ellipse. For CGCS2000, it simplifies to (units: $\mathrm{m}{\mathrm{s}}^{-2}$)

$${\gamma }_p=9.7803253361\left(1+0.00530244{\sin }^{2}B-0.00000582{\sin }^{2}2B\right)$$

(9)

The error of this formula is less than $0.1\times {10}^{-5}\mathrm{m}{\mathrm{s}}^{-2}$.

3. Results

3.1. Analysis of Displacement Characteristics in the Study Area

Between 2018 and 2022, we selected a total of 87 descending Sentinel-1A images. The SBAS-InSAR technique was employed to investigate the LOS displacement and vertical displacement characteristics within the study area. Interferometric pairs exceeding the spatiotemporal baseline threshold were excluded, resulting in 425 interferometric pairs that met the spatiotemporal baseline conditions (Figure 5c). After estimating deformation rates and residual topography, the corrected and unwrapped phases were merged and transformed into elevation data. These data were then geocoded under the WGS84 coordinate system, with the LOS deformation values subsequently being converted to vertical deformation values [61]. The outcomes reveal the time-series deformation quantities for the Line-of-Sight (LOS) (Figure 5b) and vertical directions (Figure 5a), as well as the deformation quantities curve in Figure 5d. An analysis of Figure 5a, b indicates that the vertical and LOS deformation magnitudes and intervals exhibit a high degree of consistency. The LOS deformation interval ranges from −51 mm to 30 mm, while the vertical deformation interval ranges from −64 mm to 37 mm. The average deformation in the study area fluctuates between −5 mm and 4 mm. This deformation exhibits a predominantly stable trend, with significant variations in localized regions [62,63,64].

Figure 5. InSAR data processing baseline distribution and deformation information. (a,b) Illustrate the vertical and LOS deformation rates. (b) The ellipses colored in red, gray, magenta, pink, purple, and blue each signify six distinct micro-regions representing varying terrain deformation patterns. (c) The red lines delineate interferometric image pairs, while the yellow pentagram marks the master image and the green pentagrams represent respective temporal Sentinel-1A images. (d) Showcases six unique curves that correspond to the deformation changes in the six regions depicted in (b).

Figure 5. InSAR data processing baseline distribution and deformation information. (a,b) Illustrate the vertical and LOS deformation rates. (b) The ellipses colored in red, gray, magenta, pink, purple, and blue each signify six distinct micro-regions representing varying terrain deformation patterns. (c) The red lines delineate interferometric image pairs, while the yellow pentagram marks the master image and the green pentagrams represent respective temporal Sentinel-1A images. (d) Showcases six unique curves that correspond to the deformation changes in the six regions depicted in (b).

Upon further scrutiny of six subregions, deformation quantities for Region1 through Region6 were discerned to range from −6 to 2 mm, 0 to 12 mm, −5 to 22 mm, −1 to 4 mm, −1 to 1 mm, and 0 to 45 mm, respectively. Analyzing these in conjunction with on-site topography and high-resolution optical remote sensing images (Figure S1) revealed that Regions 2 and 6, characterized as lakes and marshlands, exhibited pronounced deformations. Conversely, Regions 3 and 4, identified as dune clusters, exhibited noticeable deformation changes. Regions 1 and 5 encapsulate the overarching environmental conditions of the study area, demonstrating a relatively stable deformation trend, as illustrated in Figure 9(a1–a4). Surface horizontal deformation displacement for the periods 2020–2021, 2021–2022, and 2020–2022 was obtained using the COSI-Corr software (cosi-corr_pakOct19), as shown in Figure S1.

By comprehensively analyzing the displacements in the LOS, vertical, and horizontal directions (Figure S1) obtained from the integration of InSAR and high-resolution optical remote sensing imagery, it can be inferred that, although multiple faults are distributed within the study area and there have been earthquakes of magnitude less than 5.0 MW, the displacements in the LOS, vertical, and horizontal directions have been relatively small in recent years. This suggests a relatively stable geological plate and fault activity in the region [31]. Moreover, the predominant landforms in the area consist of plains and small hills, with no extensive high mountains or canyons. These conditions collectively contribute to the favorable enhancement of precision for the subsequent refinement of the quasi-geoid surface using the GNSS gravity-potential leveling model, as detailed in the following sections.

3.2. High-Precision POS Data Post-Differential Processing Accuracy Analysis

After the completion of two flight missions, a Post-Processing Kinematic (PPK) approach was employed [65,66,67] to differentially process the synchronized observations from the GNSS reference stations [67,68] with the acquired POS data. The IE software (Version number: Inertial Explorer 8.7) was used to perform both forward and reverse recursive Kalman filtering fusion for all exterior orientation elements, including the IMU parameters (roll φ, pitch Ω, and yaw κ) and the GNSS coordinates (X, Y, and Z). The GNSS reference station data were sampled at 1-s intervals, while the IMU was sampled at 200 Hz. The processing epochs ranged from 196,000 to 212,000, as shown in Figure 6a. Observations from four types of satellites, namely GPS, GLONASS, BeiDou, and Galileo, were collected. The Double Difference Dilution of Precision (DD_DOP) values for all observation epochs were between 0.8 and 1.8 (Figure 6b). The partial satellite observation residuals, as listed in Figure 6c, were within the range of −0.04 m to 0.04 m, indicating good and continuous dynamic observations of the satellites during the flight without significant cycle slips or data gaps. This ensured high data quality. Furthermore, an analysis of the PAV80 three-axis stabilized platform’s performance during flight revealed timely corrections made to the roll and pitch angles. The aircraft’s attitude angles were between −30° and 25° during the turn when switching flight lines. After the correction, the roll and pitch angles were within the range of −3° to 3°. Upon transition to the East (E), North (N), and Up (U) directions, the residuals in the E and N directions fluctuated between 0.04 and 0.06 m, while in the U direction, the residuals varied between 0.1 and 0.26 m, as shown in Figure 6h. The evident corrective effect demonstrated that the PAV80 three-axis stabilized platform effectively compensated for the instability in the aircraft’s attitude caused by external wind direction and airflow during flight, providing solid support for obtaining high-precision IMU data (Figure 6d).

Figure 6e decomposes the aircraft’s velocity during flight into horizontal (X, Y) and vertical (Z) components. In the X and Z directions, velocity fluctuations are observed around 0 m/s, while in the Y direction, representing the aircraft’s forward motion, the velocity component remains between 110 m/s and 125 m/s. After transformation to E, N, and U coordinates, the residuals fluctuate between 0.12 cm/s and 0.16 cm/s (Figure 6i). Additionally, the aircraft’s flight altitude consistently remains between 7420 m and 7445 m, closely matching the originally designed altitude of 7430 m (Figure 6f), with a deviation of approximately 25 m. These observations affirm that the aircraft maintains high precision in both speed and altitude throughout the flight, providing assurance for the accurate acquisition of stereo imagery. Upon transforming the horizontal (X, Y) and vertical (Z) components into E, N, and U directions, we conducted an analysis of the integrated POS accuracy from IMU (φ, Ω, κ) and GNSS (X, Y, Z) data fusion [38,68]. We found that the fusion error in the E direction fluctuates between −0.1 m and 0.04 m; in the N direction, it ranges from −0.06 m to 0.06 m, and in the U direction, it fluctuates between −0.26 m and 0.45 m. All of these errors were found to be below the empirical threshold of 0.5 m, as shown in Figure 6g. For the first and second flights, a total of 312,353 (Figure 3c) and 72,116 (Figure 3c) POS coordinate data points were calculated, respectively. The aforementioned high-precision POS data provide robust support for subsequent aerial triangulation densification, stereo image model reconstruction, and the collection of three-dimensional coordinates, which are detailed in the following sections.

3.3. Airborne Triangulation Accuracy Enhancement Analysis

3.3.1. Analysis of Ground Control Point Layout Schemes

In order to integrate the high-precision POS data calculated in Section 3.2 for aerotriangulation refinement and the establishment of stereoscopic imagery models of the study area, a uniform GCP layout scheme, as described in Section 2.3, was employed across the survey region. A total of 46 GCPs were deployed, with the baseline distance between adjacent control points maintained within 50 km, as shown in Figure 7a,g–i, depicting the challenging environments within the study area where GCP deployment and measurements were conducted. Upon completion of the measurements, all GCPs were processed to the CGCS2000 coordinate system using the Gamit/Globk software (Version number: Gamit/Globk 10.71), in conjunction with IGS station data, employing the un-differenced ionosphere-free model for PPP [68,69], achieving centimeter-level precision for all resolved GCPs. By leveraging the coordinates of these GCPs, aerotriangulation refinement was performed on the XPro photogrammetric workstation to establish a three-dimensional stereoscopic imagery model of the study area.

In order to evaluate the impact of different GCP layout schemes on the precision of aerotriangulation refinement, and to select the most optimal and cost-effective layout, the 46 GCPs were categorized into five distinct schemes. Scheme one (Figure 7b) involved using the POS data alone for network adjustment, treating all 46 GCPs as checkpoints to assess the adjustment accuracy. Scheme two (Figure 7c) utilized eight GCPs, with the remaining 38 serving as checkpoints. Scheme three (Figure 7d) increased the number of GCPs to 16, with 30 as checkpoints, while Scheme four (Figure 7e) and Scheme five (Figure 7f) incrementally increased the number of GCPs to 31 and 40, respectively, leaving 15 and 6 as checkpoints.

The intent of Scheme two was to distribute points around the outer perimeter and central region of the study area to control the precision of the aerotriangulation, and schemes three to five progressively added more control points to ascertain whether a greater number of GCPs indeed contributes to enhanced precision in aerotriangulation. The rationale behind selecting these five schemes was to analyze the minimum number of GCPs required for the study area. This approach aimed to minimize the fieldwork in uninhabited regions under adverse environmental conditions, ensuring the safety of the survey personnel, while still satisfying the precision requirements for aerotriangulation refinement.

Figure 7. Spatial distribution of GCPs and checkpoints for aerial triangulation. (a–f) Illustrate the various GCP layout schemes, with red triangles representing the GCPs and green triangles denoting the checkpoints. Depict the deployment of GCPs in various environments within the study area. (g–i) Measurement conditions in three different geographical environments.

Figure 7. Spatial distribution of GCPs and checkpoints for aerial triangulation. (a–f) Illustrate the various GCP layout schemes, with red triangles representing the GCPs and green triangles denoting the checkpoints. Depict the deployment of GCPs in various environments within the study area. (g–i) Measurement conditions in three different geographical environments.

3.3.2. Analysis of Airborne Triangulation Accuracy under Various GCP Layout Schemes

By incorporating the GCP deployment strategy discussed in Section 3.3.1, we conducted a beam method regional network adjustment for aerial triangulation densification using the XPro photogrammetric workstation [38]. Following the adjustment, the errors for the control points and junction checkpoints (JCPs) for the aforementioned deployment schemes were obtained. For the control point errors, we performed coordinate stabbing, which involved transferring the actual measured locations to their corresponding positions on the images, and carried out a combined regional network adjustment for all control points. This process resulted in obtaining the positional Mean Square Error (MSE) for each control point. Regarding checkpoint errors, upon the combined adjustment of the control point regional network, we transferred the coordinates of the points not involved in the adjustment to their corresponding positions on the images. Then, we conducted a combined adjustment for all checkpoints to determine their positional MSE [38,41,65].

As illustrated in Figure 8, in Scheme I, control points were not included in the adjustment, resulting in a zero MSE, while checkpoints exhibited a horizontal MSE of 0.189 m and a vertical MSE of 0.123 m (Figure 8a,b; Table S1). In Scheme II, the inclusion of ground control points increased the horizontal MSE for control points and checkpoints to 0.320 m and 0.235 m, respectively, and increased the vertical MSE to 0.235 m and 0.114 m (Figure 8c,d; Table S2). This initial estimate suggests that the introduction of GCPs in the regional network adjustment may introduce a systematic bias during the rectification of images from the image space coordinate system to the CGCS2000 coordinate system of the control points. In Scheme III, with an increased number of ground control points, the horizontal MSE for the control points and checkpoints was further reduced to 0.119 m and 0.184 m, respectively, with vertical MSE values of 0.115 m and 0.106 m. (Figure 8e,f; Table S3). This indicates an enhancement in the precision of aerial triangulation densification with the increasing number of control points. Continuing the increment in Scheme IV, the horizontal MSE for control points and checkpoints was 0.170 m and 0.150 m, respectively, with a vertical MSE of 0.100 m for control points and 0.133 m for checkpoints (Figure 8g,h; Table S4). The MSE for control points in the horizontal direction slightly increased, whereas it decreased in the vertical direction.

3.4. Leveling Points Accuracy Enhancement Analysis

In order to evaluate and verify the precision of the GNSS Gravity-Potential Leveling proposed in this study, level measurements were carried out in the study area. Two leveling routes were established and measured, as shown in Figure 9a. The first leveling route was located in the southwest of the study area, where six leveling points were measured using a second-order leveling standard. The second route, located in the northeast, traversed the entire study area, and it consisted of 22 leveling points measured at the same second-order leveling standard. Upon completion of the leveling measurement, PPP was employed at the same leveling points to conduct geodetic coordinate measurements, using the data processing model proposed in Section 2.6. This resulted in two sets of elevation coordinates for the same leveling points [70,71,72]. The first set is based on the 1985 elevation datum normal height coordinates system, and the second set is based on the CGCS2000 national elevation datum geoidal height coordinate system.

In the application of the GNSS gravity-potential leveling delineated in Section 2.5, gravity-potential values of 593 GNSS points were extracted in a stereo environment compared with the average normal gravity value (Table S6). These calculated normal heights were then subjected to an RMS accuracy analysis against the measured normal heights, as per the methodologies discussed in [72]. Prior to this, the gravity-potential values at the leveling points were computed using the gravity-potential model described in Section 2.6.2 (Figures S2 and S3). The geopotential value at the elevation reference point, specifically the sea-level datum at the Qingdao tide gauge, as shown in Figure 1, was set at 62,636,854.205 m²/s², in accordance with [56,70]. Further, the mean gravity values at each leveling point were calculated using the average normal gravity model outlined in Section 2.6.3. In order to analyze the impact of various global gravity models on the precision of normal height calculations for leveling points in the study area, we selected 12 different global gravity models. The spherical harmonic coefficients of the five models (EGM2008, SGG-UGM-1, SGG-UGM-2, Eigen-6C4, and GECO) were extended up to degree and order 2159, as supported by the works of Liang et al. [13], Zingerle et al. [14], Gilardoni et al. [15], and Chisenga et al. [16].

Figure 9. Comparison and analysis of normal heights from leveling route measurements and GNSS gravity-potential calculations. (a) The black lines represent the two leveling routes surveyed in the study area, with the leveling points denoted by yellow circles. The red lines indicate the designed flight paths, while the pale yellow ellipses mark the representative leveling point locations. Subfigures (a1–a4) illustrate GNSS measurements taken simultaneously at leveling points across different geographic environments. (b–m) represent the Root Mean Square Error (RMSE) values of leveling points computed using twelve geopotential models, namely EGM2008, SGG-UGM-2, SGG-UGM-1, Eigen6C4, GECO, XGM2019e-2159, EGM96, GGM05C, GO-CONS-GCF-2-DIR-R6, Tongji-GMMG2021S, EigenCG03C, and Tongji-GRACE02, respectively.

Figure 9. Comparison and analysis of normal heights from leveling route measurements and GNSS gravity-potential calculations. (a) The black lines represent the two leveling routes surveyed in the study area, with the leveling points denoted by yellow circles. The red lines indicate the designed flight paths, while the pale yellow ellipses mark the representative leveling point locations. Subfigures (a1–a4) illustrate GNSS measurements taken simultaneously at leveling points across different geographic environments. (b–m) represent the Root Mean Square Error (RMSE) values of leveling points computed using twelve geopotential models, namely EGM2008, SGG-UGM-2, SGG-UGM-1, Eigen6C4, GECO, XGM2019e-2159, EGM96, GGM05C, GO-CONS-GCF-2-DIR-R6, Tongji-GMMG2021S, EigenCG03C, and Tongji-GRACE02, respectively.

Figure 9 illustrates the residual values for the leveling points computed by each global gravity model. Notably, the EGM2008 global gravity model yielded the lowest RMSE value of 14.83 cm, as shown in Figure 9b, followed by the SGG-UGM-2 and SGG-UGM-1 models with RMSE values of 17.79 cm and 18.15 cm, respectively, as shown in Figure 9c,d. Models with an RMSE value below 20 cm included Eigen6C4 and GECO, registering values of 18.75 cm and 18.97 cm, as shown in Figure 9e,f. The models XGM2019e-2159, EGM96, and GGM05C had RMSE values of 21.88 cm, 22.88 cm, and 24.96 cm, respectively, as shown in Figure 9g–i. The GO-CONS-GCF-2-DIR-R6, Tongji-GMMG2021S, and EigenCG03C models exhibited larger RMSEs of 36.11 cm, 37.05 cm, and 39.00 cm, respectively (Figure 9j–l). A preliminary estimate indicates that the lower accuracy of these models may be attributed to their relatively lower degree of spherical harmonic expansion. The Tongji-GRACE02 model showed the largest RMSE of 43.02 cm (Figure 9m) due to its harmonic coefficients being expanded only to the 180th degree, the lowest among the models analyzed.

3.5. GNSS Gravity-Potential Accuracy Enhancement Analysis

As concluded in Section 3.4, the normal heights computed using the GNSS gravitational geopotential height leveling in this study exhibit superior RMSE precision over the QPQM regional geoid model (18.6 cm). This is particularly evident with the EGM2008 (14.83 cm), SGG-UGM-2 (17.79 cm), and SGG-UGM-1 (18.15 cm) global gravity models, as shown in Figure 9b–d. This finding validates the feasibility of the GNSS leveling method proposed in this paper, which is grounded on the actual precision of the leveling points. By utilizing the theoretical model of GNSS gravitational geopotential height leveling as described in Section 2.6 this study calculates the normal heights for all 594 points (pertaining to GNSS gravity-potential control points) collected under the three-dimensional models outlined in Section 2.4, based on the CGCS2000. In order to enhance the accuracy of the GNSS leveling method, a grid error correction approach is employed [7,34,56], wherein the normal heights obtained from the leveling measurements, H_1985, are used to correct the normal heights derived from the GNSS leveling models, H_GNSS.

After corrections, we calculated the residuals of the normal height (H_GNSS) and the RMSE values for the three-dimensional solid model collected GNSS points using 12 gravity satellite models. Post-correction, both the residuals (Figure 10a) and the precision of the RMSE values exhibited significant improvements. The EGM2008 gravity satellite model showed the minimum residuals and RMSE values (Figure 10b), which were 6.90 cm and 8.61 cm, respectively. Subsequently, the SGG-UGM-2 and SGG-UGM-1 gravity satellite models presented residuals and RMSE values of −8.10 cm, −7.70 cm, 9.09 cm, and 9.38 cm, respectively, as shown in Figure 10a,d. Gravity satellite models with residuals in the −20 cm to 20 cm range included Eigen6C4 and GECO, while those with RMSE values in the 10 cm to 20 cm range comprised Eigen6C4, GECO, and XGM2019e-2159, as shown in Figure 10e,f. The EGM96 and GGM05C models yielded residuals and RMSE values of 16.70 cm, −22.30 cm, 23.74 cm, and 25.59 cm, respectively, as shown in Figure 10h,i. The models G0-CONS-GCF-2-DIR-R6, Tongji-GMMG2021S, and EigenCG03C recorded residuals and RMSE values of −14.60 cm, −15.70 cm, −38.20 cm, 29.98 cm, 36.18 cm, and 39.91 cm, respectively, as shown in Figure 10. Lastly, the Tongji-GARCE02 gravity satellite model had residuals and RMSE values of 42.90 cm and 43.02 cm, respectively (Figure 10m).

Accordingly, it can be deduced that the EGM2008, SGG-UGM-2, and SGG-UGM-1 gravity geoid models, after refinement with level data corrections, achieve centimeter-level accuracy, which is superior to the original regional quasi-geoid model, QPQM. This demonstrates the feasibility and significant improvement in the precision of the methodology proposed in this paper. The methodology employs the ADS80 tri-linear array push-broom stereo imaging process for aerial triangulation enhancement in conjunction with GNSS gravimetric geoid model data fusion for the further refinement of the quasi-geoid in the study area. Furthermore, based on precision requirements, other global gravity models can be selectively utilized to establish a regional quasi-geoid surface. For example, when the accuracy requirement is between 10 cm and 20 cm, the Eigen6C4, GECO, and XGM2019e-2159 models are suitable candidates. For precision between 20 cm and 30 cm, the EGM96, GGM05C, and G0-CONS-GCF-2-DIR-R6 models may be adopted. In the range of 30 cm to 40 cm accuracy, the Tongji-GMMG2021S and EigenCG03C models could be selected. Lastly, for accuracy demands above 40 cm, the Tongji-GRACE02 gravity satellite model would be appropriate.

4. Discussion

4.1. Influence of Different Gravity Field Models on the Accuracy of Regional Quasi-Geoid

The method proposed in this article, which utilizes the ADS80 three-line array push-broom stereo imaging aerial triangulation technique integrated with GNSS geopotential-level model data processing for orthometric height computation, to a certain extent, can replace traditional leveling measurements, and it facilitates the rapid restoration and establishment of regional quasi-geoid surfaces. It can be inferred from Table S7 that the higher the degree of spherical harmonic expansion in the adopted gravity satellite model (ranging from degrees 180 to 2190), the higher the accuracy of the solution provided by the GNSS geopotential leveling method [17,19,21]. An accuracy analysis using leveling data measured on-site within our study area revealed that the precision of the orthometric heights computed by this method surpass that of the existing QPQM quasi-geoid surface. Moreover, after the correction of the leveling data, the accuracy of GNSS points collected via the stereo model based on the EGM2008, SGG-UGM-2, and SGG-UGM-1 gravity satellite models can achieve centimeter-level precision. This enables the rapid computation and restoration of the 1985 orthometric height datum in the study region. This is achieved without relying on ground leveling or terrestrial gravity data [19,22,71]. This is particularly crucial for providing high-accuracy regional quasi-geoid models in the uninhabited, high-altitude areas of the Tibetan Plateau, thereby offering robust theoretical and precision support.

4.2. Impact of Various Errors on GNSS Gravity-Potential Leveling

The GNSS gravity-potential leveling discussed in this article is subject to various errors during the computation process [18,70]. In particular, the King Air 350 aircraft, equipped with an ADS80 photogrammetric camera, incurs flight errors due to multiple flight factors, including wind speed, wind direction, air currents, and flight attitude, especially when navigating a “figure-eight” trajectory, as described in Section 2.3 Errors also arise from the measurement of GCPs and the processing inaccuracies within the Gamit/Globk software (Version number: Gamit/Globk 10.71), as well as from the discrepancies in post-processing differential PPK POS data calculations (as set out in Section 3.2, with a 50 km baseline among three neighboring GNSS reference stations). Errors associated with POS-assisted aerial triangulation densification and pinpointing (various control point deployment strategies are detailed in Section 3.3.1), the reconstruction of three-dimensional stereo model point collection (three modes of stereo image display are analyzed in Section 2.4), and GNSS geopotential leveling resolution errors (including zero-order and tidal terms) are further considered. A detailed accuracy analysis concerning these various errors was conducted to ensure that the precision of the original data can substantiate and bolster the computational accuracy of the GNSS geopotential leveling method proposed herein [71,72]. As the accuracy of flight platforms, photogrammetric cameras, GNSS satellites, and gravity satellite hardware platforms improves, along with the stabilization and refinement of modeling and computational approaches, a reduction in the aforementioned errors is anticipated. This reduction is expected to provide even more favorable assurance for enhancing the precision of the GNSS geopotential leveling method elaborated in this paper.

4.3. Characteristics of Refined Regional Quasi-Geoid Models

Upon refining the normal heights calculated by the GNSS gravimetric geoid determination method using leveling data, we derived the regional elevation anomaly values $\xi =H_P-H_P^{\gamma }$ (Figure 1). Subsequently, these ξ values facilitated the reconstruction of the quasi-geoid model within our study area, as shown in Figure 11. It is evident that the surface of elevation anomalies modeled by QPQM is consistent with the surface derived from the GNSS gravimetric geoid method using gravity satellite data, as seen in Figure 11; this consistency is observed both in terms of trend distribution and magnitude. They exhibit a pattern where the eastern region is lower and the western region higher, with elevation anomalies ranging between −60 m and 54 m, indicating minimal variation. This aligns with the conclusion in Section 3.1 that our study area is dominated by desolate, arid ‘yardang’ landforms. It is sparsely populated and relatively flat, with stable LOS and vertical ground deformation trends over time and minimal fault-related horizontal displacements.

The quasi-geoid model within the research area, recovered using $\xi$ values, underwent the extraction and analysis of residual values (Figure 12). It was observed that the interval of residual values for elevation anomalies falls within the range of −0.556 m to −0.595 m. For the three gravity satellite models, EGM2008, SGG-UGM-2, and SGG-UGM-1, the residuals of the computed elevation anomalies lie between −0.102 m and 0.104 m. These are lower than those derived from other gravity satellite models, according to their respective accuracies, as discussed in Section 3.5. This further underscores that, for the study area in this paper, utilizing the GNSS gravity-potential leveling to analyze the regional refinement precision enhancement effect of the quasi-geoid, EGM2008, SGG-UGM-2, and SGG-UGM-1 are the optimal choices among the three gravity satellite models. Other gravity satellite models can be selectively chosen based on varying precision requirements.

As our focus is on the normal height system under the 1985 elevation datum of China, an analysis of the ξ values was conducted. However, the analysis revealed that the GNSS gravimetric geoid model proposed in this paper is equally applicable to the calculation of orthometric heights within the WGS84 global elevation datum framework [7,57]. The quasi-geoid undulation is expressed as the height from the quasi-geoid to the reference ellipsoid, denoted as $N=H_P-H_P^g$, and the gravitational potential of the WGS84 geoid, $W_0$ (62,636,856.00 ± 0.5 m²/s²), along with the mean gravity, $g_m^P$ [7,58,59,60], can be calculated using the aforementioned 12 global gravity models. This facilitates the establishment of a transformation relationship with the global elevation datum.

5. Conclusions

In our study area situated in the unmanned, high-altitude region of the Tibetan Plateau, we conducted an accuracy analysis of stereoscopic imagery acquisition and aerial triangulation densification using the ADS80 aerial camera. This analysis was coupled with empirical data from 26 leveling points and 12 GNSS gravity-potential models to synthesize and analyze the feasibility and accuracy improvement effects of this method in refining the regional quasi-geoid. It was found that among these gravity-potential models, including EGM2008, SGG-UGM-1, EIGEN-6C4, and GECO, the models were able to achieve decimeter-level RMSE in the orthometric heights. These heights were derived from the 593 GNSS points collected by the stereoscopic image model after aerial triangulation densification with respect to the China 1985 Yellow Sea Datum and the leveling point heights. Specifically, the EGM2008, SGG-UGM-2, and SGG-UGM-1 gravity-potential models displayed the highest precision, with RMSEs of 14.83 cm, 17.79 cm, and 18.15 cm, respectively, which are comparable to the precision of the QPQM model.

Upon comprehensive analysis, it has been determined that for regional quasi-geoid refinement without correction through leveling data, the following gravimetric models can be selected based on the desired accuracy levels. For precision requirements between 10 cm and 20 cm, the Eigen-6C4, GECO, and XGM2019e-2159 models are suitable; for a precision range of 20 cm to 30 cm, the EGM96, GGM05C, and GO-CONS-GCF-2-DIR-R6 models can be employed; for precision between 30 cm and 40 cm, the Tongji-GMMG2021S and Eigen-CG03C models are recommended; and for accuracies exceeding 40 cm, the Tongji-GRACE02 model is advisable. After the leveling data are corrected using the weighted average method, the computed 1985 elevation accuracy of each gravimetric model significantly improves, with the EGM2008, SGG-UGM-2, and SGG-UGM-1 models achieving precisions of 8.61 cm, 9.09 cm, and 9.38 cm, respectively, all surpassing the precision of the QPQM model. The enhancements range from 9.22 cm to 9.99 cm in magnitude, indicating that the methods proposed in this paper can significantly improve the accuracy of the regional quasi-geoid. This also provides a new approach for the rapid establishment of regional quasi-geoid models in the unmanned areas of the Tibetan Plateau.