Impacts of GCP Distributions on UAV-PPK Photogrammetry at Sermeq Avannarleq Glacier, Greenland

Abstract

Real-Time/Post-Processing Kinematic (RTK/PPK) technology has been widely applied in Unmanned Aerial Vehicle (UAV) photogrammetry in glaciological research. Considering that ground control points (GCPs) cannot be set on glaciers, evaluating the impacts of one-sided distribution is essential. In this study, 8571 images were captured at Sermeq Avannarleq glacier in western Greenland from 4 August 2021 to the 6th, covering approximately 85 km², with the furthest distance being 13.22 km away from the coastline. Benefited by the meandering coastline, 11 roving stations roughly uniformly distributed on bare rock were surveyed with the RTK technique. PPK-geotagged images were processed in Agisoft Metashape Professional to derive the DSMs, utilizing eight different configurations of GCP distributions that gradually extended longitudinally (along the glacier flow direction) to the upper part of the glacier. The accuracy of DSMs was evaluated by referring to the validation points (VPs) that were not employed in the Bundle Block Adjustment (BBA). The results indicated that the RMSE values of the easting, northing, and height of the reconstruction model georeferenced by only PPK geotagging (no GCPs applied) were 0.038 m, 0.031 m, and 0.146 m, respectively. Applying four GCPs located at one side of the region but with both longitudinal and lateral distribution improved the RMSE values in easting, northing, and vertical to 0.037 m, 0.031 m, and 0.081 m, respectively, and these values were stable even when distributing four GCPs evenly or when increasing the number of GCPs to eleven. Moreover, the cross-validation with ICESat-2 and ArcticDEM performed only at an off-glacier region also suggested that vertical accuracy shows significant improvements for every configuration of GCPs compared to the reconstruction model optimized only by PPK, but such improvements were not obvious if the number of GCPs exceeded four. Moreover, no elevation ramps appeared in the UAV DSM, even for the GCP configuration with only two GCPs distributed at the terminus. Therefore, combining PPK with only a few GCPs but distributing in both directions of the surveying region can offer a viable solution for obtaining glacier DSMs at the coastline with decimeter-level accuracy.

1. Introduction

Over the past several decades, aerial photogrammetry has experienced rapid advancements. A Structure-from-Motion (SfM) algorithm is capable of generating various high-resolution photogrammetric derivatives including point clouds, 3D reconstruction models, Digital Surface Models (DSMs), Digital Terrain Models (DTMs), and orthomosaics by detecting and matching common feature points across multiple images [1]. Unmanned Aerial Vehicles (UAVs), or drones, are becoming progressively favored due to their ability to provide rapid, sensitive, cost-effective, and repeatable measurements, particularly for small-scale areas [2]. For instance, drones are being extensively utilized for monitoring fault zones [3], conducting glacial and periglacial geomorphology [4,5,6], mapping landslides [7], monitoring flood hazard [5], topographical mapping [8], agricultural investigation [9], ecological studies [10], etc. Drones have become increasingly popular because of their ability to capture detailed and precise data, thereby enhancing efficiency and effectiveness in various applications across diverse fields [11].

The Greenland Ice Sheet (GrIS) plays a crucial role as a freshwater reservoir, and it is the primary driver of the global sea level rise in recent decades. Variations in its size, thickness, and melting pattern have the potential to impact future water resources and human well-being [12]. Fast-flowing outlet glaciers terminating in the ocean account for 88% of Greenland ice discharge into the ocean, significantly influencing the evolution of the ice mass balance of the GrIS [13]. Glacier dynamics and changes in surface structures indicate that systemic imbalances within these glaciers could result in rapid changes in the future. The monitoring of glacier dynamics and the assessment of their stability rely on observations of multiple factors, including flow speed, surface melt, ice front position, and grounding line location. Therefore, high spatiotemporal glacier terrain data play a crucial role in evaluating their influence on sea-level rise and response to climate change. Common methods involve the comparative analysis of multi-temporal satellite remote sensing datasets or spatiotemporal modelling of ice dynamics based on field observations [14]. However, the traditional remote sensing observation and modelling methods applied in these studies often exhibit lower spatial and/or temporal resolution, thus hardly meeting the demands of studying rapid changes in marine-terminating glaciers, especially in the summer, and leading to an inadequate understanding of the influence of glacier evolution. High-resolution targeted observations are essential for understanding the impacts of glacier dynamics.

Unmanned Aerial Vehicle (UAV) photogrammetry offers the potential to compensate for the gap between field observations and satellite remote sensing at the GrIS. This approach yields higher spatial resolution (i.e., cm level resolution) and enhanced temporal resolution, which is critical for studying GrIS changes, especially in the summer. UAV photogrammetry has been successfully used in previous studies to survey glaciological processes (e.g., quantifying ice melting [15,16]; mapping front calving and crevasse dynamics [4,8]; measuring debris-covered glaciers’ changes [17,18]; tracking glaciers’ flowing motion [15,19]).

The precision in terms of positions and orientations of UAV-captured images is crucial to ensure the accuracy of UAV photogrammetric products. It is usually achieved through two methods: the indirect georeferencing approach using both the Bundle Block Adjustment (BBA) method and ground control points (GCPs), and the direct georeferencing method using BBA without GCPs. While the former involves performing aerial triangulation (AT) using GCPs, the latter involves adjusting the 3D reconstruction model involving the positions and orientations of arial images obtained from Real-Time/Post-Processing Kinematic (RTK/PPK) techniques onboard the drones [20]. Both the surveying of GCPs and of the camera perspective center coordinates employ the RTK method, which involves acquiring real-time position data from satellites and reference GNSS stations. However, in the case of a very large region survey, in which the drone may lose its connection to the reference station, the PPK method still provides accurate locations of the camera perspective centers [21], eliminating the reliance on communication links and avoiding latency issues.

To date, several studies have employed RTK/PPK without GCPs for direct georeferencing, yielding topographic products with high accuracy. When evaluating the accuracy of photogrammetric products within forest terrain, PPK methods demonstrated superior vertical accuracy compared to traditional indirect approaches with four GCPs, but similar accuracy to approaches with nine GCPs [22]. Žabota and Kobal [23] explored the application of UAVs on rockfalls, aiming to investigate the accuracy of photogrammetric products optimized by PPK and GCPs. Results indicated that the total RMSE values for the three sites were 0.056, 0.066, and 0.305 m with PPK, compared to 0.019, 0.036 and 0.014 m obtained with GCPs. Cho and Lee [21] employed PPK and GCPs to investigate accuracy decline in areas lacking in RTK signals. Without PPK, errors tended to escalate away from the GCPs. However, with PPK, no significant disparity was noted, with horizontal errors of up to 10 cm and vertical errors of up to 20 cm within a 1 km radius of the GCP. Consequently, with PPK, even if RTK signals were lost, achieving a vertical error of the check point (CP) within 3 ground sample distances (GSDs) (20 cm) was feasible up to 1 km from the GCP. While PPK exhibits slightly lower accuracy compared to traditional indirect methods, it remains a viable solution for mapping hazardous and exposed areas. As a result, georeferencing without GCPs emerges as a practical alternative to traditional methods.

GCPs can further improve the performance of photogrammetric products optimized with PPK. A previous study [24] examined the effects of the number and spatial distribution of GCPs on elevation accuracy. The study also evaluated how RTK/PPK positioning can enhance georeferencing while reducing the required number of GCPs. The results indicated that while a minimum of twelve GCPs is necessary for reliable georeferencing in the GCP-only scenario, the GNSS-assisted (RTK/PPK) method requires only six GCPs to achieve comparable accuracy. Another study [25] presents an experimental investigation using three different block control configurations: GCP only, RTK-obtained camera stations only, and a combination of camera stations with one GCP. The findings reveal that regardless of the RTK mode, the first and third configurations deliver the highest consistency within the DSMs, indicating the importance of GCPs. Hugenholtz et al. [26] performed field testing on three sets of orthoimages and DSMs produced for comparing spatial accuracies: two sets were created using direct georeferencing images from the RTK UAV and non-RTK UAV, and another set was created by using GCPs during the external orientation of the non-RTK UAV images. The results indicate that direct georeferencing using an RTK UAV can achieve horizontal accuracy similar to that obtained with a network of GCPs. However, for topographic measurements where the highest level of accuracy is essential, GCPs are still essential for establishing the DSM.

As discussed above, indirect georeferencing using GCPs is still a better choice due to its higher accuracy. In the case of indirect georeferencing, although the preparation and execution of a UAV flight can be conducted from a safe location, the collection of GCPs must be carried out directly within the area of interest. Most current studies have concentrated on regions with flat or easily accessible terrain, where conducting field GCP measurements is feasible. However, the interest in this research lies in understanding the applicability of UAVs in the rapidly flowing glacier terminus at the margin of the GrIS. Unlike other areas, GCPs are generally unavailable on a glacierized surface due to its moving and it being hard to access, so GCPs are usually distributed at off-glacier regions [6]. In such cases, it is crucial to evaluate the accuracy trends in regions far away from the GCPs and to evaluate the feasibility of generating a topographic model using only PPK-based direct georeferencing or by utilizing a minimal number of easily accessible GCPs.

The primary objective of this study is to analyze the impact of the number and distribution of GCPs on the accuracy of photogrammetry-derived products on glaciers. This research compares the accuracy of photogrammetric models and DSMs generated using two georeferencing methods: (1) direct georeferencing utilizing only coordinates of photo perspective centers obtained through the PPK technique, and (2) indirect georeferencing employing both PPK-based photo coordinates and field-collected GCPs (seven different GCP distribution configurations). The findings aim to balance the workload and accuracy of aerial photogrammetry in polar ice sheet studies.

2. Materials and Methods

2.1. Study Site

Sermeq Avannarleq glacier (69°24′N, 49°55′W), located in the southwestern part of Greenland (Figure 1), is a typical marine-terminating outlet glacier. It is situated immediately north of the Jakobshavn Isbrae glacier (69°09′N, 50°18′W), which is one of the largest marine-terminating glaciers in Greenland both in terms of volume discharge and flow speed [27]. However, due to its high speed and extensive area, it is challenging and costly to conduct large-scale flights at Jakobshavn Isbrae Glacier. We opted for the nearby similarly conditioned and relatively easily accessible Sermeq Avannarleq glacier as the experimental area. The total experimental area covered by aerial photos is approximately 85 km².

As shown on Landsat-8 satellite imagery (3 August 2021) of the glacier, both sides of the glacier terminus feature exposed rocks. The meandering coastline allows for the collection of GCPs (red pentagrams in Figure 1) on the stable areas with a roughly uniform distribution.

2.2. Methods

The research workflow comprises five main steps (Figure 2):

(1) Preparation for the aerial survey mission in the study area, including setting a reference GNSS station and planning a flight route.

(2) With one GNSS receiver placed at the reference station conducting static continuous observations, aerial photos were collected using an RGB camera mounted on a UAV, and the synchronized real-time GNSS information was stored in base files under the PPK mode. Another GNSS receiver on the ground collected data on 10 RTK roving stations utilizing RTK technology.

(3) Reference station GNSS information was utilized in the WingtraHub software (Version 2.0) to geotag the camera centers of the aerial photos using the PPK technique.

(4) Geotagged aerial photos were first aligned through SfM by identifying and matching feature points in the overlap between images, and then processed using two georeferencing methods (PPK direct and PPK+GCPs indirect) to reconstruct the surface model and generate Digital Surface Models (DSMs) and Digital Ortho Models (DOMs) in Agisoft Metashape photogrammetry software (Version 1.7).

(5) The precision of the 3D reconstruction models was assessed with the locations of VPs (Validation Points) in different configurations of GCP distributions (parts of the RTK roving stations were designated as GCPs and the remainder as VPs). Moreover, external cross-validation of the UAV DSMs was performed using ArcticDEM and ICESat-2 elevation products.

2.2.1. UAV Platform, Flight Planning, and Aerial Photo Acquisition

This flight mission utilized the high-accuracy surveying and mapping UAV WingtraOne (Figure 3), a hybrid between fixed-wing and multirotor designs with a 125 cm wingspan. It weighs 3.7 kg and flies for a maximum of 1 h at a 16 m/s flying speed. It vertically takes off and lands (VTOL), and it switches into a fixed-wing mode in flight. The drone carries a payload consisting of a 42 MP (megapixel) CMOS sensor Sony RX1RII full-frame camera equipped with a 32.8 mm fixed lens. Technical specifications of the flight platform are listed in

Table 1. Technical specifications of the UAV system.

UAV Body WintraOne		Payloads Sony RX1R II
Drone type	VTOL	Sensor type	Full-frame
Weight	3.7 kg	Sensor size	35 mm × 23.3 mm
Max. payload weight	800 g	Mega pixel	42.4 (7952 × 5304)
Endurance	59 min	Focal lens	32.8 mm

. The system is equipped with a Multi-Frequency Multi-Constellation PPK module which can store the position information of the camera trigger events in the output data to achieve cm level absolute accuracy after post-processing. The UAV equipped with Sony RX1R II enables absolute accuracy (RMSE) horizontally down to 1 cm and vertically down to 3 cm with PPK and without GCPs, according to the WingtraOne manufacturer (https://wingtra.com/mapping-drone-wingtraone/drone-survey-accuracy/, accessed on 26 September 2024).

The flight mission was performed over a period spanning three days, from 4 August 2021 to the 6th (Figure 1). Since the most severe height changes (melting) usually occur at the glacier terminus, which is also important for the flux gate method for evaluating glacier mass balances, the study area focuses on the glacier terminus and the neighboring off-glacier region suitable for GCP placement, also considering the limitations of UAVs and RTK in terms of communication range and field workload. The originally planned study area was rectangular, covering the whole glacier terminus; however, due to climatic conditions encountered during field operations, we had to abandon the uncompleted tasks in the upper glacier region. As a result, the final study area appeared cruciform, as shown in Figure 1. WingtraPilot software (Version 2.0) was used to construct the optimal flight path to cover the study area, with an along-path overlap of 80% and cross-path overlap of 60%, respectively. Moreover, the flight altitude remained at the pre-set 300 m above ground to achieve a GSD of 5 cm/pixel.

A reference GNSS station (50°14′09.14″N, 69°21′59.9″W, yellow pentagram in Figure 1) using a Trimble R10 GNSS receiver was set on stable bare rock to conduct a static survey session and to serve as a reference site for the UAV-PPK and RTK roving stations. Its location was calculated using Continuous GNET (Greenland GNSS Network), which includes the sites Ilulissat (ILUL) (69°22′21.38″N, 49°54′49.33″W), Kangia North (KAGA) (69°13′20.28″N, 49°48′52.56″W), Rinks Isbræ (RINK) (71°50′54.6″N, 50°59′38.4″W), and Aasiaat (AASI) (70°43′34.07″N, 55°15′26.46″W), with the D-GNSS technique.

To assess the accuracy of photogrammetric products reconstructed from PPK-corrected images at the glacier terminus, 10 RTK roving stations (Figure 1) together with a reference GNSS station were positioned on bare bedrock along the glacier edge. These stations were surveyed using a Trimble R10 RTK receiver, with locations determined in reference to the GNSS station mentioned above (Figure 3b). The RTK receiver guarantees horizontal and vertical accuracies of 8 mm + 0.5 ppm and 15 mm + 0.5 ppm, respectively, in the RTK mode, according to the manufacturer [28]. Considering that the maximum distance from the roving sites to the reference site (Figure 1) is ~8 km, the horizontal and vertical accuracy should be less than 12 mm and 19 mm, respectively.

2.2.2. Photogrammetric Processing

A total of 8571 aerial photos were captured with the raw GNSS data and base files (RINEX) data of the reference GNSS station during the flight. WingtraHub (https://knowledge.wingtra.com, accessed on 26 September 2024) software for PPK image geotagging was used to correct the GNSS data of the camera centers of the flights.

These geotagged images were then processed using Agisoft Metashape Professional (version 1.7, https://www.agisoft.com, accessed on 26 September 2024) based on the SfM algorithm. With images and their PPK-corrected locations/orientations taken as input, the software automatically extracts and matches corresponding features within the overlapping area of multiple photos, forming a network of tie points (the same point locations on the ground that appear in multiple overlapped images). The 3D coordinates of these points are then determined through a least squares network optimization technique known as BBA. In the BBA process, GCPs can be incorporated as additional external constraints to further refine the network. By minimizing the reprojection error and the difference between the projected features and their corresponding image positions, the BBA method optimizes both the internal and external orientation parameters of the images within the 3D reconstruction model. Optimal alignment can be achieved, and camera location errors were estimated in the process. In this case, 8448 of 8571 images were aligned. The magnitude of these errors served as an indicator of the reliability and accuracy of the camera locations.

To assess the accuracy of the 3D reconstruction models optimized without GCPs and with different GCP distributions, we carried out 8 GCP configurations, providing external control to optimize the model (Figure 4). The 1st configuration took none of the RTK roving stations (also the reference GNSS site; same for the following parts) as GCPs to externally control BBA, but regarded all RTK roving stations as validation points (VPs) to verify the accuracy of the 3D reconstruction model. The 2nd to 8th configurations used 2/4/7/9/11 of RTK roving stations located at the terminus of the glacier as GCPs to control the BBA process, and took the rest as VPs to explore how the error distributes relative to the distance and terrain (Figure 4). While the distribution of increasing numbers of GCPs extended from the glacier terminus inward for the 2nd to 6th configurations, two other configurations (7th and 8th) had 4 and 7 GCPs evenly distributed (referred to as 4GCPs-E and 7GCPs-E, compared to the uneven distributions, referred to as 4GCPs-T and 7GCPs-T) across the study area (Figure 4). The remaining RTK roving stations (VPs) were used to investigate the impact of the number and distribution of GCPs on the 3D reconstruction model. For the configuration of 4GCPs-T, the furthest surveyed area is about 10 km to the nearest GCP. The error between the coordinates of VPs in the 3D reconstruction model computed by Agisoft Metashape and the field-collected VPs’ coordinates is employed to assess the 3D reconstruction model’s accuracy in every configuration.

Dense cloud points were built based on the SfM-MVS (Multi-View Stereo) image matching algorithm, which operates on a single-pixel scale. In this step, we set the ‘Quality’ parameters in Agisoft Metashape to ‘High’, which allowed the software to create dense cloud points from original-resolution images. With dense cloud points taken as input, the DSM and orthomosaic could be produced. Here, we set ‘Interpolation’ to ‘Enabled’ and ‘Filtering’ to ‘Medium’, which allowed us to smooth the surface of dense cloud points and to correct the altitude of points using neighborhood points’ altitudes and interpolate altitudes of dense cloud points to produce the DSM, after which the DOM was generated through digital differential correction based on the DSM [1]. For all 8 GCP configurations, the GSD of the DSM (Figure 5a) and the DOM (Figure 6a) in GeoTIFF format was consistently set to 10.9 cm and 5.4 cm, respectively. All DSM products and subsequent analyses were based on the WGS-84/Stereographic North Polar (EPSG:3413) projection coordinate system, and the elevation was set in refererence to the WGS 1984 ellipsoid.

2.2.3. Accuracy Assessment with Different GCP Configurations

The accuracy of the 3D reconstruction models in different GCP configurations was evaluated mainly by VPs. Firstly, we marked the appearances of each RTK roving site in every photo to identify the location of each site in the 3D reconstruction model, and RTK roving sites were accordingly regarded as GCPs or VPs for 8 different GCP configurations. Then, the computed coordinates of the roving sites in the 3D reconstruction model were retrieved through Agisoft Metashape automatically, and the differences between reference coordinates measured with the GNSS RTK receiver and computed positions could be compared. $X$, $Y$, and $Z$ represent easting, northing and ellipsoidal height, respectively. The coordinate differences between reference coordinates (X_R, Y_R, Z_R) and those computed by the construction model (X_C, Y_C, Z_C) were further used to calculate the root mean square error (RMSE) following Equation (1), to evaluate the accuracy of the photogrammetric model:

$${RMSE}_D=\sqrt{{\displaystyle \frac{{\sum }_1^N{\left(D_R-D_C\right)}^2}{N}}},D\in X,Y,Z$$

(1)

where D represents directions including X, Y, Z, $D_R-D_C$ represents the difference between the reference coordinates $D_R$ from the field measurements and the photogrammetry-computed coordinates $D_C$ from the 3D reconstruction model, and N represents the number of points involved.

To assess the systematic errors among DSMs produced in different GCP configurations, we also conducted inter-comparisons of multiple photogrammetry-derived DSMs. This evaluation aimed to determine the discrepancies or the spatial ramp of height, particularly for the area far away from the GCPs.

Additionally, to explore how the number and distributions of GCPs influence the accuracy of the total photogrammetric model, we checked the spatial distribution of VP coordinates’ errors in each direction to analyze whether or not errors were correlated with the terrain or the distance between GCPs and VPs.

2.2.4. External Cross-Validation of UAV Photogrammetric DSMs

This study also employed ICESat-2 and ArcticDEM for external elevation accuracy validation of UAV photogrammetry-derived DSMs.

The ATLAS (Advanced Topographic Laser Altimeter System) onboard of the ICESat-2 satellite is a six-beam laser altimeter. Each beam employs a 532 nm wavelength (green) laser with a pulse repetition frequency of 10,000 Hz, leading to an along-track ground spacing of approximately 0.7 m and a ground segment diameter of about 11 m [29,30]. According to the NSIDC and previous studies [30], the L3A level land ice height dataset ATL06 provides geolocated land-ice surface heights (above the WGS-84 ellipsoid, ITRF2014 reference frame) with centimeter-scale precision along 40 m segments of ground track, spaced 20 m apart, and it is more accurate and precise than any prior altimetric product. We employed two ground tracks that closely coincided with the time of the UAV flight mission and overlapped with the study site to validate the accuracy of photogrammetry-derived DSMs. These two reference ground tracks, numbered 870 and 1230, were observed on 17 August 2021 and 12 September 2021, 11 days and 37 days after the UAV flight, respectively (Figure 6a). Due to the temporal fluctuation in surface elevations caused by glacier flow and dynamic imbalances, we exclusively employed off-glacier points situated in stable terrain for accuracy assessment. Points with low quality (atl06_quality_summary = 1) were also excluded from the analysis.

The ArcticDEM strip is a time-stamped DSM derived from a series of optical stereo images (Worldviews) using open-source photogrammetry software SETSM [31] (Surface Extraction with TIN-based Search-space Minimization) by the Polar Geospatial Center. We selected an ArcticDEM strip that covers the study area, which was observed 21 September 2021, 26 days after the UAV flight mission (Figure 6b). Also, only the off-glacier stable terrain was utilized to cross-compare with UAV photogrammetry-derived DSMs. Based on a previous study [31], the ArcticDEM shows an internal precision of better than 0.5 m, but some strips may exhibit systematic horizontal or vertical biases at several meters. Therefore, we applied the surface-based co-registration method proposed by Nuth et al. [32] to align the ArcticDEM strip to the photogrammetry-derived DSM. Two DSMs of the same terrain surface that are not perfectly aligned exhibit a specific relationship of elevation differences correlated to the direction of the terrain (aspect). This method employs a multiple regression model to correct this 3D shift bias between two DSMs based on their geometric relationship [32].

As for the accuracy evaluation metrics, we utilized the RMSE (Equation (1), as mentioned in Section 2.2.3), mean (Equation (2)), standard deviation (STD, Equation (3)), and mean absolute error (MAE, Equation (4)) of differences in elevation between UAV photogrammetry-derived DSMs and ICESat-2 or ArcticDEM at off-glacier regions to assess the performance of UAV photogrammetry-derived DSMs.

$${Mean}_Z={\displaystyle \frac{Z_{DSM}-Z_E}{N}}$$

(2)

$${STD}_Z=\sqrt{{\displaystyle \frac{{\sum }_1^N{\left({diff}_Z-\overline{{diff}_Z}\right)}^2}{N}}}$$

(3)

$${MAE}_Z={\displaystyle \frac{1}{N}}{\sum }_1^N\left|Z_{DSM}-Z_E\right|$$

(4)

where $Z_{DSM}$ and $Z_E$ represents elevations of photogrammetric DSMs and an external elevation data source, and in this case, $Z_E$ is the elevation of the ICESat-2 altimeter or ArcticDEM, ${diff}_Z$ represents the difference between $Z_{DSM}$ and $Z_E$, and $N$ represents the number of points or pixels involved.

3. Results

3.1. Camera Location Error

Before introducing GCPs to optimize the 3D reconstruction model, the software performs image alignment by SfM and estimates the internal and external orientation parameters of each image. The camera errors representing the RMSE of the difference between camera input coordinates (EXIF imported) and estimated coordinates are shown in

Table 2. Average camera location RMSE in BBA using only PPK image coordinates.

X RMSE (m)	Y RMSE (m)	Z RMSE (m)	XY RMSE (m)	Total RMSE (m)
0.024	0.060	0.275	0.064	0.282

. X, Y, and Z represent easting, northing, and ellipsoidal height direction, respectively. The direct georeferencing method that relies only on PPK image coordinates achieves a horizontal RMSE at the centimeter level (XY RMSE = 0.064 m), but the vertical RMSE persists at the decimeter level (Z RMSE = 0.275 m). This indicates that direct referencing based on the PPK technique is reliable in terms of horizontal accuracy, but less precise in the vertical direction. Therefore, it is essential to introduce external GCPs to optimize the 3D reconstruction model.

3.2. Accuracy Assessment with Different GCP Configurations

3.2.1. 3D Reconstruction Model Accuracy Assessed by VPs

Table 3. Accuracy of the 3D reconstruction models in different GCP configurations.

	GCPs				VPs				RTK Roving Sites
Configuration/RMSE (cm)	X	Y	Z	XYZ	X	Y	Z	XYZ	X	Y	Z	XYZ
PPK (no GCPs)	\	\	\	\	3.83	3.11	14.66	15.44	3.83	3.11	14.66	15.44
2GCPs	0.07	0.11	0.11	0.16	3.36	3.53	11.27	11.58	1.87	3.13	9.04	10.38
4GCPs-T	0.07	0.12	0.16	0.27	3.74	3.15	8.11	9.48	1.81	2.19	4.45	5.44
7GCPs-T	0.93	0.83	0.27	1.28	1.87	3.62	7.44	8.47	1.33	2.01	3.67	4.13
9GCPs	0.48	0.67	0.73	1.21	1.91	3.52	7.48	8.33	1.14	1.58	3.03	3.73
11GCPs	0.43	0.43	0.21	1.02	\	\	\	\	0.43	0.43	0.21	1.02
4GCPs-E	0.28	0.34	0.30	0.90	2.23	3.29	7.14	8.80	1.29	2.04	3.72	4.39
7GCPs-E	0.38	0.26	0.53	0.82	1.67	3.21	7.04	7.94	1.34	1.41	3.92	4.51

shows the RMSE of the 3D reconstruction model in reference to GCPs, VPs, and RTK roving stations (GCPs and VPs). As shown above, when the study area lacks GCPs and relies solely on PPK-based camera locations for georeferencing, the RMSE from 11 VPs is 0.038 m in easting, 0.031 m in northing, and 0.146 m in height. When optimized by two GCPs, the VPs’ horizontal accuracy of the 3D reconstruction model remains nearly constant, while the VPs vertical RMSE decreases to 0.112 m, indicating a limited improvement in vertical accuracy. In configuration 4GCPs-T, increasing the number of GCPs located in the terminus to four, with distributions both across and along the direction of glacier flow, results in a further reduction of the VPs’ vertical RMSE to 0.081 m. However, further increasing the longitudinal distribution (along the glacier flow direction) of GCPs does not significantly enhance the vertical accuracy of the 3D reconstruction model in the cases of 7GCPs-T and 9GCPs. Moreover, in the cases of 11GCPs, 4GCPs-T, and 7GCPs-T, the even distribution of GCPs across the study area only leads to a minor improvement of VPs’ vertical RMSE, which stabilizes at 0.07–0.08 m. As for the internal precision of the 3D reconstruction model evaluated by GCPs, the horizontal and vertical RMSE of GCPs are all at the millimeter-level.

3.2.2. Spatial Distribution of the 3D Reconstruction Model’s Accuracy

The errors of the GCPs and VPs are plotted in Figure 7, as is the spatial distribution of the accuracy of VPs’ errors in different GCP configurations. The magnitudes of the errors between different points vary significantly, so to better display the differences, the size of the vectors represents the denary logarithm of the original error values. To explore how GCP distributions influence the accuracy of the spatial pattern of the photogrammetric model, when no GCPs are used to optimize the model, the horizontal errors of VPs are almost randomly distributed across the study area and remain at the cm level, but the vertical errors are much higher at the decimeter-level. With an increase in the number of GCPs and an expanded longitudinal distribution range, when the horizontal accuracy of the same VP shows minimal and random changes, there is a noticeable improvement in the vertical accuracy. Additionally, VPs closer to the GCPs typically exhibit higher vertical accuracy. Furthermore, even though the number of GCPs located at the terminus increases to 7 and 9 with a longitudinal distribution extension, the remaining VPs show only a slight difference in accuracy compared to the 4GCP-T configuration. For GCP evenly distributed configurations 11GCPs, 4GCPs-E, and 7GCPs-E, the horizontal and vertical errors of VPs show no distinct improvement as well compared to the 4GCPs-T configuration. Consequently, with an appropriate distribution of four GCPs in both longitudinal and lateral directions within the survey area, the vertical accuracy of the 3D reconstruction model is effectively controlled. Moreover, this vertical accuracy tends to stabilize even with an increase in the number of GCPs or in a more even distribution.

3.2.3. Inter-Comparison between UAV-DSMs with Different GCP Configurations

To assess whether there was a ramp error in elevation in the glacier interior part away from the terminus where GCPs distribute, we inter-compared and analyzed DSMs produced in all GCP configurations. Figure 8 presents the spatial distribution of elevation differences between DSMs derived from various GCP configurations. The differences in these DSMs highlight the impact of the number and distribution of GCPs on the vertical accuracy of the DSM. Figure 9 shows histograms of elevation differences between various DSM configurations. The statistical information in

Table 4. Inter-comparisons of UAV photogrammetry-derived DSMs among different GCP configurations.

Comparison	MEAN (m)	MAE (m)	STD (m)	RMSE (m)
DSM11GCPs−DSMPPK	0.048	0.129	0.141	0.146
DSM11GCPs−DSM2GCPs	0.029	0.907	0.116	0.119
DSM11GCPs−DSM4GCPs-T	0.004	0.049	0.093	0.095
DSM11GCPs−DSM7GCPs-T	0.003	0.051	0.084	0.091
DSM11GCPs−DSM9GCPs	0.011	0.043	0.089	0.096
DSM11GCPs−DSM4GCPs-E	0.008	0.040	0.079	0.086
DSM11GCPs−DSM7GCPs-E	0.004	0.041	0.086	0.092
DSM4GCPs-T−DSM4GCPs-E	0.009	0.055	0.083	0.089

further supports the visual results from the figures.

Based on the elevation differences between DSM_11GCPs and DSM_PPK, as well as DSM_11GCPs and DSM_2GCPs, it is evident that the two GCPs, positioned only at the glacier terminus, improve the vertical accuracy of the DSM compared to PPK, with MAE and RMSE reducing from 0.129 m and 0.146 m to 0.097 m and 0.119 m, respectively. However, when examining the elevation differences between DSM_11GCPs and DSM_4GCPs-T, employing four GCPs distributed laterally and longitudinally is more effective in enhancing DSM elevation accuracy, lowering the MAE and RMSE further to 0.049 m and 0.095 m. Furthermore, comparisons between DSM_11GCPs and DSM_4GCPs-T as well as DSM_7GCPs-T, and DSM_9GCPs show minimal differences across the study area, with consistent MAE and RMSE values around 0.05 m and 0.095 m, even with the addition of more longitudinal GCPs. The differences are uniformly distributed, indicating no DSM ramp errors despite asymmetrical GCP placement. Additionally, when assessing the influence of relatively evenly distributed GCPs, the comparison between DSM_11GCPs and DSM_4GCPs-T as well as DSM_11GCPs and DSM_7GCPs-E reveals a similar pattern in both distribution and magnitude, with MAE and RMSE values at 0.04 m and 0.09 m, respectively. Furthermore, comparing the results between DSM_4GCPs-T and DSM_4GCPs-E, it is evident that the improvement of DSM_4GCPs-E is not substantial comparing to DSM_4GCPs-T.

From the discussion above, while increasing the number of GCPs enhances DSM accuracy, the most significant improvement is achieved by ensuring a balanced distribution of GCPs both longitudinally and laterally within the study area. Adding GCPs beyond a certain number (four GCPs in this case) does not lead to substantial gains in accuracy, as evidenced by the nearly identical elevation difference patterns and minimal changes in RMSE and MAE values. Additionally, the uniformity of the elevation differences across the entire region implies no ramp errors, even if the GCPs are asymmetrically placed. Thus, four evenly distributed GCPs ensure the accuracy of the DSM; however, when the surveying condition does not allow for an even distribution, four GCPs with certain lateral and longitudinal distributions at the terminus still effectively control the quality of the derived DSM.

3.3. Cross-Validation with ICESat-2

The statistics and distribution histograms of elevation differences in off-glacier areas between UAV photogrammetry-derived DSMs (4 August 2021 to the 6th) and ICESat-2 ATL06 land ice elevations (17 August 2021 and 12 September 2021) based on 314 samples are presented in

Table 5. Comparisons of off-glacier ICESat-2 ATL06 points’ elevations and UAV photogrammetric DSMs at off-glacier regions with different GCP configurations.

Comparison	MEAN (m)	MAE (m)	STD (m)	RMSE (m)
IS2−DSMPPK	0.201	0.217	0.689	0.916
IS2−DSM2GCPs	0.179	0.194	0.608	0.883
IS2−DSM4GCPs-T	−0.052	0.141	0.541	0.533
IS2−DSM7GCPs-T	−0.046	0.144	0.498	0.509
IS2−DSM9GCPs	0.031	0.129	0.504	0.492
IS2−DSM11GCPs	0.036	0.133	0.519	0.551
IS2−DSM4GCPs-E	0.041	0.136	0.477	0.529
IS2−DSM7GCPs-E	−0.037	0.140	0.511	0.497

and Figure 10, respectively. The difference in MAE between DSM_PPK and ICESat-2 is 0.217 m, with an RMSE of 0.916 m, showing a decimeter-level deviation. For DSM_2GCPs, an MAE of 0.194 m and an RMSE of 0.883 m demonstrate that the lateral distribution of two GCPs only slightly improves vertical accuracy compared to DSM_PPK. The differences between DSM_4GCP-T and ICESat-2 are at the decimeter level, with an MAE of 0.141 m, and the difference is concentrated more densely around 0 m, indicating a substantial accuracy improvement compared to DSM_PPK and DSM_2GCPs. As the number of longitudinally placed GCPs gradually increases in DSM_7GCPs-T and DSM_9GCPs, the differences between ICESat-2 and UAV photogrammetric DSMs show similar patterns and statistics, indicating that the addition of longitudinal GCPs brings a minimal improvement in accuracy compared to the DSM_4GCPs-T.

Considering that the GCPs in the configurations discussed above were placed asymmetrically across the study area, additional comparisons were made between DSMs derived from different GCP configurations and ICESat-2 elevations to assess whether further GCP placement improves accuracy. The MAE and RMSE differences between ICESat-2 and UAV photogrammetric DSMs DSM_11GCPs, DSM_4GCPs-E, and DSM_7GCPs-E stabilize at ~0.135 m and ~0.530 m, marginally outperforming DSM_4GCP-T. This indicates that an even distribution of four GCPs can ensure DSM vertical accuracy, but when even distribution is not feasible, placing four GCPs with both lateral and longitudinal distribution near the glacier terminus still effectively maintains the derived DSM quality, which is similar to the conclusion reached when validating the VPs.

3.4. Cross-Validation with ArcticDEM

A comparison of ArcticDEM (2 September 2021) and UAV photogrammetric DSMs (4 August 2021 to the 6th), as presented in Figure 11, reveals elevation differences, particularly between the glacier ice and off-glacier exposed bedrock areas. The glacier ice shows a noticeable elevation decrease, potentially due to summer ice thinning or tidal influences, while the off-glacier terrain remains more stable. Therefore, the comparison focuses on the off-glacier bedrock region and covers over-five-million-pixel samples. Figure 12 and

Table 6. Cross-comparison statistics of ArcticDEM and different UAV photogrammetry-derived DSMs in off-glacier region.

Comparison	MEAN (m)	MAE (m)	STD (m)	RMSE (m)
ArcticDEM−DSMPPK	−0.314	0.612	1.760	1.513
ArcticDEM−DSM2GCPs	−0.303	0.557	1.682	1.579
ArcticDEM−DSM4GCPs-T	−0.094	0.481	1.483	1.466
ArcticDEM−DSM7GCPs-T	−0.106	0.524	1.562	1.439
ArcticDEM−DSM9GCPs	−0.096	0.495	1.521	1.454
ArcticDEM−DSM11GCPs	−0.112	0.477	1.517	1.422
ArcticDEM−DSM4GCPs-E	−0.101	0.485	1.509	1.431
ArcticDEM−DSM7GCPs-E	−0.097	0.512	1.531	1.454

show the histograms and statistics of the elevation differences between ArcticDEM and UAV photogrammetry-derived DSMs in the off-glacier region.

Both DSM_PPK and DSM_2GCPs show similar systematic errors in elevation, with a difference in MAE of ~−0.6 m and in RMSE of ~1.55 m compared to ArcticDEM (Table 6). The introduction of four GCPs, distributed both laterally and longitudinally, significantly reduces the MAE to 0.481 m. This improvement is accompanied by a tighter clustering of the elevation differences around 0 m, suggesting that the use of four GCPs markedly enhances DSM elevation accuracy. Beyond this point, even with further increases in longitudinal GCPs, the MAE stabilizes at around 0.5 m across different configurations, indicating minimal additional accuracy gains. The elevation difference patterns remain consistent, demonstrating that the use of four GCPs suffices for achieving high-quality control over the elevation accuracy of DSM products, and further GCPs beyond this threshold, particularly longitudinal ones, contribute little to further improvement.

The comparison of evenly distributed GCP (DSM4_GCPs-E and DSM_7GCPs-E) configurations relative to ArcticDEM reveal that although DSM_11GCPs demonstrates the highest accuracy, the differences between DSM_11GCPs and ArcticDEM are only slightly better than DSM_4GCPs-E and DSM_7GCPs-E. The MAE in these configurations stabilizes at approximately 0.5 m, with minor improvements in RMSE across these three configurations. This indicates that increasing the number of evenly distributed GCPs beyond four does not significantly enhance DSM accuracy. Similar to the results of cross-validation with ICESat-2 and validation with VPs, 4GCPs and 7GCPs configurations with either side or even distribution suggest no obvious difference. Despite the slight advantage of even GCP distribution, the terminus-based configuration still ensures relatively high DSM quality.

4. Discussion

4.1. Influence of GCP Distributions

Previous studies [33,34,35] have consistently shown that an evenly distributed and dense network of GCPs was effective in controlling the accuracy of UAV photogrammetric products. In glacial environments, having at least one control point in the central region of the glacier can significantly enhance vertical accuracy [19]. However, due to the inherent dangers, inaccessibility, and instability of the glacier’s interior, stable GCPs are typically only obtainable on relatively flat terrain near the terminus or margins of the ice sheet.

This study demonstrates that for a survey area of ~100 km², and when only PPK-geotagged image coordinates are used for photogrammetry processing and without constraining with GCPs, the 3D reconstruction model achieves a horizontal RMSE within 4 cm and a vertical RMSE within 15 cm as verified by VPs. The inclusion of 4 GCPs located near the glacier terminus (4GCP-T) with a certain longitudinal and lateral distribution successfully improves accuracy, reducing the horizontal and vertical RMSE of VPs to within 4 cm and 8 cm, respectively. The maximum distance between the furthest area and the nearest GCP is less than 8 km. In a case in which the 4 GCPs are distributed more evenly across the survey area (4GCPs-E) and even more GCPs are incorporated (7GCsP-E), the improvement in horizontal and vertical RMSE is insignificant.

The comparison between DSM_4GCPs-T and ICESat-2 revealed that the precision is comparable to DSM_11GCPs, DSM_7GCPs-E, and DSM_4GCPs-E, exhibiting MAE and RMSE values of approximately 0.14 m and 0.53 m, respectively. This indicates that the DSM optimized using four GCPs positioned near the glacier terminus (DSM_4GCPs-T) achieves a similar level of vertical accuracy to that of DSMs utilizing more and evenly distributed GCPs. Despite the fewer GCPs and their asymmetric distribution, the quality of DSM_4GCPs-T remains on par with DSM_11GCPs, which has GCPs spread more evenly across the survey area. Moreover, no evidence of elevation ramp bias was detected across the DSM_4GCP-T-surveyed region. This suggests that even without the use of additional GCPs or an evenly distributed GCP configuration, the 4GCPs-T configuration effectively controls elevation errors, maintaining consistent accuracy throughout the entire region. These findings suggest that in cases where glacierized area does not allow for the setting of GCPs, four GCPs asymmetrically distributed across the study area (~100 km²) can still yield high-quality DSM products with minimal compromises in accuracy.

Establishing an extensive network of GCPs on marine-terminating glacial terrain can be challenging or costly. Obtaining only a few GCPs near the glacier terminus is more feasible than establishing a large, evenly distributed GCP network. This research eliminates the requirements for an evenly distributed network of GCPs or GCPs located at the center of the survey area, addressing the significant logistical challenges, especially at marine-terminating calving glacier fronts.

Usually, the glaciers on the GrIS and/or Antarctic ice sheets can be divided into land- and marine-terminating types. For the former type, the wide off-glacier region can serve as a zone for GCPs with both longitudinal and lateral distribution, while for the latter type, if the glacier is not too wide, the GCPs distributed at the off-glacier zones at two sides of the glacier can still guarantee accuracy. The DSM’s accuracy still needs to be investigated if the glacierized region for surveying is too far from the GCPs (>10 km).

In case no GCPs are available, the PPK technique still helps to significantly correct the positions of the photo perspectives, providing high accuracy, with vertical and horizontal RMSE values within 4 cm and 15 cm as verified by VPs for photogrammetric 3D reconstruction, and producing terrain products such as a DSM with a decimeter-level vertical RMSE compared to ICESat-2 (Table 5).

4.2. Limitations of the Study

The quality of photogrammetry relies on UAV image matching; most image matching algorithms are based on image features consisting of spatial distributions of local grayscale gradients, and therefore, they are sensitive to illumination intensity and angle variation resulting from changes in solar position [36]. Glacial ice surfaces often exhibit low contrast, limited texture, and high reflectance from fresh snow, and so slight variations in illumination on the ice surface can result in significant local changes in grayscale, leading to inaccurate surface construction [15]. In this study, adverse weather conditions led to insufficient illumination during UAV flights over the upper part of the glacier, resulting in uneven brightness in the orthomosaic (Figure 5a), and making color calibration challenging. Therefore, when applying UAVs to glacier investigations, careful consideration should be given to factors such as illumination conditions and photo overlapping.

Due to the instability and unreliability of the glacier ice surface, validation using GCPs directly on the ice surface was not feasible for this and similar studies. Although cross-validation of drone DSMs was performed with ArcticDEM and satellite altimetry, there was a time gap between UAV flight and ArcticDEM’s original image (WorldViews) acquisition, and ice melting is uncertain, especially in the summer season. Moreover, like UAV photogrammetry-derived DSMs, ArcticDEM originated from WorldViews stereo pairs and was also produced utilizing photogrammetric methods, leading to similar systematic errors. Therefore, here, the ArcticDEM can only be used for cross-validation or comparison with UAV photogrammetry-derived DSM and cannot serve as a completely accurate reference for assessing the vertical accuracy of UAV DSMs.

Moreover, due to the limited battery capacity and challenging field operations in glacial environments, the surveyed area in this study is approximately 85 km², with a maximum longitudinal extent of only 13.22 km, which is very small compared to the Greenland ice sheet. The RTK/PPK technique relies on differential GNSS based on signals from the same satellites under similar space conditions. In the case of a further inland extension of the study area, the distance between the UAV-mounted GNSS receiver and the reference receiver increases, leading to a degradation in PPK correction effectiveness. Moreover, the longitudinal extent covered by the distribution of 4GCPs-T is approximately 16% of the total longitudinal depth of the survey area. If the survey area extends further inland the GrIS, this proportion will decrease, and the possible accuracy loss requires further study. Therefore, the conclusions regarding accuracy control in this study are applicable only to situations with areas smaller than or equal to this surveyed area.

4.3. Suggestions for a UAV Survey Working Plan at Marine-Terminating Glaciers

The geometric accuracy of photogrammetry relies significantly on the georeferencing approach. Experimental results from this study indicate that within a ~100 km² survey area, relying on RTK/PPK-based UAV and incorporating a small number of GCPs, even if only located on one side of the survey area, is sufficient to achieve decimeter-level vertical accuracy in DSM production. Therefore, when conducting ice sheet surveys, it is highly recommended to utilize drones equipped with RTK/PPK modules to reduce the need for field collection of GCPs, ensuring economy and efficiency in fieldwork.

Based on the results of our experiment, we have summarized some economical and reliable suggestions for obtaining terrain products like DSMs in polar research:

(i) Photo coordinates correction: Although studies [21] suggest that RTK and PPK techniques yield similar improvements in GNSS measurement accuracy under favorable signal conditions, signal disconnection during the flight in the RTK technique due to long distance compromises the accuracy of photo coordinates. Unlike RTK, PPK post-processing technology, however, avoids this issue by processing the data after the flight, ensuring more reliable results [21]. For the correction of the center position coordinates of the photos, we highly recommend PPK over the RTK technique.

(ii) Flight planning: At polar regions, the low solar zenith angle results in significant illumination differences, even for short intervals during UAV flights, which can cause the direction and shape of local shade and shadow patterns to undergo considerable variation. During flight missions, it is advisable to schedule flights during periods with similar weather and lighting conditions to minimize variations in exposure and illumination across different areas. This helps prevent issues such as uneven coloration in the extraction of tie points and the resulting orthomosaic.

This study yielded the 3D reconstruction model with vertical RMSE values of ~3 GSD (0.146 m) based solely on PPK-geotagged images and ~1.5 GSM (0.081 m) with the constrain of at least 4 GCPs. Since most drone payloads select wide-angle cameras, the pixel size on the ground relies on pixel density and flight height. Then, the flight height determines the final accuracy. Here, a 300 m flight height gives a vertical accuracy of ~0.08 m.

(iii) GCPs collection: To achieve optimal accuracy in 3D reconstruction, it is recommended that GCPs be uniformly distributed across the survey area. However, in challenging conditions where collecting a sufficient number of GCPs is difficult, placing a smaller number of GCPs (four GCPs in this study) with both lateral and longitudinal distribution can still enhance accuracy. Even when GCPs are positioned on one side of the survey area, but with a longitudinal extent covering more than 16% of the surveying region, the elevation accuracy of the entire survey area can be significantly improved. In areas far from the GCPs, no significant ramp error is observed in the resulting DSM. In extremely challenging environments where it is impossible to deploy any GCPs, relying solely on the RTK/PPK technique for UAV-based glacier measurements is still feasible, with horizontal RMSE values of 0.03–0.05 m (~0.9 GSD) and vertical RMSE values of 0.146 m (~3 GSD). With the addition of GCPs at the edges of the survey area, there is minimal improvement in the already highly precise horizontal accuracy, but vertical RMSE can be enhanced from 0.146 m (~3 GSD) to 0.07–0.08 cm (~1.5 GSD). The inclusion of a few GCPs effectively controls the vertical accuracy of photogrammetric products without incurring substantial costs, which can be particularly significant for research on hard-to-reach polar ice caps.

5. Conclusions

This study conducted a photogrammetry survey at a typical marine-terminating glacier in Greenland using a VTOL UAV equipped with a PPK module. A total of 8571 images covering approximately 85 km² with a survey depth of 13.22 km at a flight altitude of 300 m were captured, achieving a GSD of 5 cm/pixel. We also obtained multiple GCPs using a GNSS receiver with RTK technology for the 3D reconstruction model quality control and validation purposes.

Evaluations of VPs (validation point) on the ground showed that the RMSEs of the 3D reconstruction model were 0.038 m (~0.9GSD) horizontally and 0.146 m (~2.6GSD) vertically without the use of any GCPs (photogrammetry only, with PPK-geotagged images for georeferencing). Considering that collecting GCPs was costly and that they should not be set on glaciers, we incrementally increased the number of GCPs from the glacier terminus to interior parts of the study sites to investigate the impacts of the numbers and distributions of GCPs on the accuracy of the 3D reconstruction models and DSMs. The results allow us to conclude that (i) the introduction of GCPs improves vertical accuracy compared to photogrammetry optimized only on images geotagged by PPK; (ii) after the number of GCPs exceeds four with a certain lateral and longitudinal distribution (longitudinal extent covering more than 16% of the surveying region), the vertical RMSE of the 3D reconstruction model stabilizes at ~0.08 m (~1.5GSD) under the flight setting of this research, and further longitudinal and/or even distribution of GCPs does not lead to a significant improvement in accuracy; and (iii) there is no DSM ramp error across the surveyed area even when four GCPs located asymmetrically at one side of the surveyed area.

Comparisons with ICESat-2 laser altimetry show that after the correction by GCPs, the MAE in elevation difference between the UAV photogrammetry-derived DSM and ICESat-2 laser points is around 0.14 m at off-glacier regions. The comparisons with ICESat-2 and ArcticDEM show that the UAV DSM vertical accuracy corrected with 4GCPs-T, 4GCPs-E and 11GCPs configurations shows imperceptible differences, and no elevation ramp error is observed across the surveyed area. Therefore, we recommend using four GCPs, located at the glacier terminus region and distributed in both directions with a longitudinal extent covering more than 16% of the surveying region, for a trade-off between DSM quality and field workload in ~100km² scale glacier terrain survey research.