Qinghai-Tibet Plateau Glacier-DEM product derived from the Chinese stereo mapping satellites

The Digital Elevation Model (DEM) is an important basic data for glacier research on the Qinghai-Tibet Plateau. With the rapid development of China 's stereo mapping satellite, it is possible to obtain high-precision DEM in the glacier area of the Qinghai-Tibet Plateau. In this study, the stereo images and laser altimetry data of ZY-3 and GF-7 satellites were comprehensively used, and the matching and joint adjustment methods of laser elevation control points and stereo images were proposed. A high-precision DEM product in glacier area was established, and the elevation accuracy of DEM was verified by using UAV data and ICESat-2 data. The results show that the accuracy of GF-7 DEM is better than 1.5 m when the slope is less than 6 degrees, and the overall accuracy is better than 2.8 m. The accuracy of ZY-3 DEM is better than 3.5 m when the slope is less than 6 degrees, and the overall accuracy is better than 5.0 m. Compared with the open source DEM with medium spatial resolution, the DEM elevation accuracy based on domestic stereo mapping satellite images is better, and the grid is finer and more detailed to describe the texture characteristics of the end of the glaciers. Compared with the high spatial resolution data set HMA DEM, the ZY-3 DEM is slightly worse and the GF-7 DEM is better, and both the ZY-3 DEM and GF-7 DEM are better than the HMA DEM in terms of coverage integrity in the Qinghai-Tibet Plateau. Chinese stereo mapping satellites can generate high-precision DEM in the glacier area, and can provide more accurate and reliable terrain reference data than the existing open free data for the study of glaciers in the Qinghai-Tibet Plateau.


Introduction
The Qinghai-Tibet Plateau, renowned as one of the largest glacier regions on Earth, holds immense significance for global freshwater reserves (Kang, 2010).The glaciers' stored freshwater not only maintains the stability of local ecosystems, but also replenishes numerous major Asian rivers, including the Yangtze and Yellow Rivers, through glacial meltwater.This meltwater is crucial for billions of downstream residents' lives, agricultural production, and general well-being (Immerzeel, 2010).However, with the escalating global warming, glaciers in the Qinghai-Tibet Plateau are melting rapidly, leading to not only a depletion of regional water reserves but also frequent glacier-related disasters.These disasters pose a significant threat to the safety of human lives, property, and the overall quality of life.Consequently, it's imperative to monitor glacier changes in the Qinghai-Tibet Plateau (Yao, 2012).DEM (Digital Elevation Model) plays an indispensable role in glacier research (Zhang, 2019).It can provide high-precision surface terrain information, including key parameters such as three-dimensional shape, volume, ice thickness change and dynamic characteristics of glaciers (Kääb, 2015).However, the Qinghai-Tibet Plateau is difficult to carry out large-scale field measurements due to its steep terrain and harsh climate.Therefore, the application of satellite remote sensing technology is necessary.For example, high-resolution optical stereo images, synthetic aperture radar interferometry (InSAR) technology and airborne / space-borne laser altimetry are used to generate and update the DEM of glaciers in the Qinghai-Tibet Plateau.In recent years, China 's stereo surveying and mapping satellites have developed rapidly.ZY-3 and GF-7 have been launched one after another, achieving 1:50,000 and 1:10,000 scale highprecision three-dimensional map respectively.In particular, GF-7, which was successfully launched in 2019, is equipped with a dual-line array stereo camera with sub-meter resolution and a high-precision laser altimeter, which has unique advantages in the grid size and accuracy of DEM products.Aiming at the problems of low grid value of open free DEM products, poor accuracy consistency of DEM products from different sources, and certain penetration error of microwave radar on snow cover on glacier surface (Zhou, 2018), the production of DEM derived from GF-7 and ZY-3 stereo satellites with higher precision and finer grid in the glacier area of the Qinghai-Tibet Plateau is discussed and validated preliminarily (Zhu, 2023).This paper introduces Chinese GF-7 and ZY-3 stereo mapping satellites, and introduces the DEM establishment and elevation accuracy optimization method in detail.Then, the Chinese satellite DEM and open source DEM data are comprehensively compared from elevation accuracy and performance details.Finally, the experimental comparison results are summarized and the establishment of large-scale glacier DEM by China's stereo mapping satellites is prospected.

Chinese stereo mapping satellites
The GF-7 was successfully launched from the Taiyuan Satellite Launch Center on 3 November 2019.It is China's first civilian sub-meter optical transmission stereo mapping satellite, operating in sun-synchronous orbit.The two line array stereo camera can effectively acquire panchromatic stereo images with 20 km swath and better than 0.8m resolution panchromatic and 3.2m multi-spectral images as well as high-precision laser altimetry data, achieving 1:10000 scale satellite stereo mapping.The sub meter stereo images of the domestic GF-7 satellite can achieve full coverage for three years and annual coverage for local areas in the Qinghai Tibet Plateau region.The ZY-3 serial satellites, a civilian high-resolution stereo mapping satellite constellation designed for 1:50000 stereo mapping, with three linear array stereo and multi-spectral cameras.These cameras include a nadir panchromatic TDI CCD camera with a ground resolution of 2.1 meters, two panchromatic TDI CCD cameras with a ground resolution of 3.5 meters for forward and backward observation, and a multispectral camera with a ground resolution of 5.8 meters.These cameras can produce 1:50000 scale mapping products.Currently, three ZY-3 satellites, named with ZY3-01, ZY3-02, and ZY3-03 has been launched in 2012, 2016, and 2020, respectively.The forward and backward camera resolution of ZY3-02 and ZY-03 has been improved from 3.5 meters to 2.7 meters, and they are equipped with laser altimeter.With these three satellites, it is possible to achieve full coverage of 2-meter stereo images of the Qinghai Tibet Plateau at least once a year, with half-year or even quarterly coverage of local areas (Li, 2019).

Selection of laser elevation control points
The laser elevation control point is derived from the GF-7 SLA03 product.SLA03 is the laser altimetry standard product after full waveform processing, geometric positioning, atmospheric and tidal correction, etc., composed of threedimensional coordinates (longitude, latitude, and ellipsoid height with respect to the WGS84), footprint images, and waveform characteristic parameters, and is the basis for extracting elevation control points and making lake water level and other products (Li, 2021).It should be noted that not all laser altimetry data can reach the nominal accuracy due to the influence of atmospheric, cloud, surface reflection, and other factors.To ensure the accuracy requirements of elevation control, certain rules have to be applied to select and extract effective points from massive laser altimetry data.In this study, the item of elevation control point flag (ECP_Flag) in laser altimetry data metadata was selected as the main filter criterion.ECP_Flag marked 1 and 2 indicate that the elevation accuracy of a point is very high or high and can be used as a elevation control point (Zhu, 2023;Chen, 2022).Considering the glacier elevation change, if the laser point is inconsistent with the stereo image acquisition time, the laser point on the glacier area is eliminated, while the laser points on the non-glacier area is retained as the elevation control point after selecting.

Matching of ECPs and the stereo images
Each GF-7 ECP provides space 3D coordinates with accurate elevation value, as well as a footprint image and the image position where laser point on the footprint image.The accurate position of the ECPs on the stereo images can be matched between the footprint image and the stereo images.The steps are follows: 1) Image feature points extraction Feature operators include Harris, LoG, SIFT and others are commonly used.Among them, LoG (Laplace of Gaussian) has the properties of stable features and insensitive to image distortion for images of different scales, sources and depths.Therefore, LoG is used as a feature operator in footprint image and stereo image matching.Obtain the overlapping area of the footprint image and the stereo image, divide the uniform grid on the overlapping area, and use the LOG operator to obtain the uniformly distributed reference points.These feature points will be used for subsequent correlation coefficient matching and least-squares matching.

2) Rough matching based on correlation coefficient
The object coordinate is calculated from the image point coordinate (x0,y0)of the feature point, and the initial position of the feature point on the stereo image (x,y) is obtained by substituting it into the stereo image imaging model.Calculate the correlation coefficient between the stereo image area centred on (x,y) and the footprint image area centred on (x0,y0) with the same search window size point by point according to Equation (1), the search window size is consistent with the laser spot,and take the maximum correlation coefficient point (x1 ' ,y1 ' ) as the pixel-level registration point: where g and g are the grey values of the footprint image and stereo image, respectively, g and g are the average grey value of the matching window, and m and n represent the width and height of the search window, respectively.Through the analysis of laser footprint image and stereo image, it can be seen that the difference between the two is about 4 ~5 times in spatial resolution.In terms of radiation, the gray value range of footprint image is quite different from that of stereo image.In order to reduce the influence of the large difference in the range on the matching, the footprint image and the stereo image should be stretched uniformly first in the matching.
3) Accurate matching based on least-squares In the actual calculation, the size of the image matching window is small, so the affine transformation relationship can be used to represent the geometric distortion between images, as shown in Equation (2).
(2) (x0,y0) and (x1,y1) is the coordinate of the corresponding points on the laser footprint image and the stereo image, (a0,a1,a2,b0,b1,b2) are the affine transform parameters.At the same time, considering the radiation deformation between images, the linear gray distortion of the laser footprint image relative to the stereo image can be expressed by the gray distortion coefficients h0 and h1 as Equation (3): is the gray value of the laser footprint image at the is the gray value of the stereo image at the (x1,y1) position.Considering the influence of noise on the images, the transformation relationship between the laser footprint image and the stereo image can be expressed as Equation ( 4) are the noise value of the corresponding position.According to Equation (4), the affine transformation parameters are iteratively calculated by the least square method.The least squares image matching method makes full use of the information in the image window for adjustment calculation, which can further improve the accuracy of image matching to the sub-pixel level and improve the accuracy of the initial matching of the correlation coefficient.Based on the correlation coefficient matching results as the initial value, the convergence of the least squares method system is solved.

Block adjustment model based on rational function model
The rational function model (RFM) uses a ratio polynomial to express the relationship between the normalized object-space coordinates (X,Y,Z) and corresponding normalized imagespace coordinates (x,y).Given its simple form and high substitution accuracy, it has become a common geometric imaging model for satellite images (Liu, 2022).Its general expression can be written as follows where P1, P2, P3, and P4 are general polynomials, the power of which generally takes three.Therefore, there are 80 rational polynomial coefficients.After obtaining the position of the ECPs on the stereo image by matching the ECPs with the stereo image, the mathematical mapping relationship between the image side and its corresponding object side is established by using the rational polynomial.According to previous scholars' research, the image compensation method has a stronger ability to eliminate systematic errors.Therefore, this paper uses ECPs and stereo images to construct an affine transformation model to realize the image compensation of stereo images and update the RPC information of stereo images.After obtaining the stereo image with improved accuracy, the DEM is extracted by image dense matching.Among them, the ZY-3 satellite stereo image generates 5m grid ZY-3 DEM, and the GF-7 satellite stereo image produces 2m grid GF-7 DEM.

3.1Glacier DEM product
The GF-7 sub-meter stereo images and laser elevation control points on the Qinghai-Tibet plateau glacier are collected to generate the high accuracy DEM product with 2 meter grid size, based on the method of this article.Almost 200,307 square kilo-meter DEM product covering the glacier and relative region has been generated by the GF-7 satellites.Which is the first DEM with 2 meter grid size in the Qinghai-Tibet plateau glacier derived from the Chinese stereo satellites.The coverage of the completed Glacier DEM products is shown on the Table 1.

Elevation accuracy validation based on UAV data
A total of 271 images of the glacier tongue of Weigele Dangxiong Glacier with an area of about 0.14 square kilometers were collected by unmanned air vehicle (UAV), and the DEM were established.The UAV DEM resolution is 1.5 cm.Compared with the GNSS data collected at the same time, the UAV DEM elevation accuracy RMSE is about 0.27 m.The UAV data is shown in the Figure 1 .
Figure 1.The UAV data of Weigele Dangxiong Glacier We collected a set of GF-7 image data and a set of ZY-3 image data in this area.The acquisition time of the GF-7 data is in September 2021, and the ZY-3 data is in December 2021.Using laser elevation control points and stereo images, the regional DEM products are established by using the method in this paper.Among them, the GF-7 DEM product uses 16 ECPs, and the GF-7 DEM resolution is 2m.While the ZY-3 DEM uses 35 ECPs, and the ZY-3 DEM resolution is 5 m.The DEM results of the glacier area are shown in the Figure 2 and Figure 3.

Comparison with open free DEM
In the study area of Gangqin Glacier and Puruogangri Glacier, ICESat-2/ATL08 products were used to evaluate the accuracy of five DEM products, as showed in Figure 6, including SRTM (Rodriguez, 2005), AW3D30 (Tadono, 2016), TanDEM, ZY-3 DEM and GF-7 DEM.   8.It can be seen that the GF-7 DEM in the study area shows the best accuracy level in all kinds of terrain; the accuracy of GF-7 DEM and ZY-3 DEM is very stable, and the error floating between different terrains is the smallest, while the TanDEM stability of 90 m grid is not high.In summary, the DEM based on GF-7 and ZY-3 satellite data has higher accuracy and stability than the open source medium spatial resolution terrain reference data.The data with lower resolution has a poor ability to express the real terrain.GF-7 DEM, ZY-3 DEM, and HMA have higher resolutions.The ability to express the real terrain is stronger, but it is obvious from the map that the HMA data in the two regions has a large part of the missing coverage is limited, and there are also many holes in the data area.

Conclusion
In this paper, the combination of laser altimetry data and stereo images is applied to the generation of DEM in glacier area for the first time, which provides a new idea for the study of glaciers in the Qinghai-Tibet Plateau.The main conclusions are as follows: (1)The results show that the accuracy of GF-7 DEM in glacier area is better than 1.5m when the slope is less than 6 degrees, and the overall accuracy is better than 2.8m.The accuracy of ZY-3 DEM is better than 3.5m when the slope is less than 6 degrees, and the overall accuracy is better than 5.0m.
(2) Compared with the existing open free DEM data in the glacier area, the Chinese satellite DEM is superior to the open free DEM (SRTM, AW3D, TanDEM), and the glacier detail texture features are more obvious.Relying on its own laser altimetry data, the DEM elevation accuracy of GF-7 glacier area is better than that of HMA DEM.Although the accuracy of ZY-3 DEM is slightly lower than that of HMA DEM, it has better coverage than HMA DEM with sparse coverage and holes.
The existing Chinese stereo mapping satellite can provide stable data for glacier scientific research in the Qinghai-Tibet Plateau.
With the accumulation of long periods, the stereo mapping satellite can also form long-term sequence data to better serve the glacier change monitoring in the Qinghai-Tibet Plateau.

Figure 6 .
Figure 6.Distribution of ICESat-2 / ATL08 laser points as reference for validation in Gangqin area A total of 1178 ATL08 laser points are used to validate the elevation accuracy of ZY-3 DEM, GF-7 DEM and SRTM, AW3D30 and TanDEM .The results are shown in the Table 2, and the slope information is derived from DEM products.The elevation accuracy of GF-7 DEM, ZY-3 DEM, SRTM, AW3D30 and TanDEM in the non-ice region is 0.77 m, 1.24 m, 3.06 m, 3.57 m and 1.72 m, respectively.GF-7 DEM shows the best accuracy level in different slope and different terrain.The 0-2°flat area and the 2-6°hilly area belong to the area with gentle terrain slope.The accuracy of GF-7 DEM is better than 1m, and the accuracy of ZY-3 DEM is better than 2.5m.With the increase of slope, the accuracy of the five DEMs in the study area has decreased, but the accuracy of GF-7 DEM is still better than that of other DEMs, followed by ZY-3 DEM, AW3D30 and SRTM.The accuracy of TanDEM is the worst, and the maximum error is 23.12 m.Table2.Elevation accuracy validation result of five kinds of DEM in Gangqin Glacier

Figure 8 .
Figure 8.The relationship between RMSE and terrain slope derived from the five kinds of DEM Figure 9 is the mountain shadow maps of GF-7 DEM, ZY-3 DEM, SRTM, TanDEM, AW3D30 and HMA(Shean,2017) in the study area, which represents the terrain surface information

Figure 9 .
Figure 9. Shadow maps among the six kinds of DEM in Gangqin Glacier

Table 1 .
The covering information of the Glacier DEM product

Table 2 .
Elevation accuracy validation result of five kinds of DEM in Gangqin Glacier