High-Resolution Elevation Model of Lop Nur Playa Derived from TanDEM-X

In this paper, a digital elevation model (DEM) was produced for Lop Nur playa produced with the data from TanDEM-X mission. The spatial resolution is 10 m. It covers an area of 38,000 km² for orthometric height from 785 m to 900 m above sea level, which is composed of 42 interferometric synthetic aperture radar (InSAR) scenes. A least-square adjustment approach was used to reduce the systematic errors in each DEM scene. The DEM produced was validated with data from other sensors including Ice, Cloud, and land Elevation Satellite (ICESat) Geoscience Laser Altimeter System (GLAS) and aerial Structure-from-Motion (SfM) DEM. The results show that global elevation root mean square error to GLAS is 0.57 m, and the relative height error to SfM DEM in complicated terrain is 3 m. The excellent height reliability of TanDEM InSAR DEM in Lop region was proved in this paper. A reliable high-resolution basic topographic dataset for researches of Lop Nur was provided.

Lop Nur used to be the terminal lake of all rivers in the Tarim Basin but had been desiccated since 1973 [1]. The remaining playa is located in the northwestern China between latitude 39°N and 41°and longitude 89°E and 92°E. It is important for studying the interaction between ancient human activities and paleoclimate change, especially after the discovery of Loulan city of Chinese Han Dynasty in 1901 by Hedin [2]. Meanwhile, his proposed wandering lake theory [2], which believed the ancient Lake Lop Nur would alternate its lake basin between depressions in this region according to elevation difference derived from erosion and accumulation [3], had been controversial among the international Earth science community for nearly a century [4]. The debate and other paleo-environment researches, such as determination of the lake boundary, discovery of paleo river channel, and building a 3D hydrodynamic model, all require high-resolution full-coverage topographic materials [5, 6].

Although some topographic surveys have been conducted, a reliable digital elevation model (DEM) over Lop Nur playa is still absent. The primitive leveling surveys [1] distributed sparely in the playa with limited height accuracy. Afterwards Differential Global Positioning System (DGPS) measurements were completed in the desiccated lake basin [7], while Ice, Cloud, and land Elevation Satellite (ICESat) Geoscience Laser Altimeter System (GLAS) elevation was estimated to be more reliable than that of DGPS in Lop Nur playa [8]. Although the height accuracy of both reached several decimeters, the spatial distribution of DGPS and GLAS was confined either along existing roads due to the reaching accessibility of vehicles or along tracks of the satellite orbits. As the DGPS and GLAS measurements show, the Lop Nur playa is quite flat with an elevation difference of 5 m from the center to the edge, and a slope of only 0.2‰. The smooth terrain requires an excellent relative height accuracy of DEM to reveal its topography. Global DEM products SRTM [9, 10] and GDEM [11, 12] both cover the Lop Nur playa, but neither the spatial resolution nor height accuracy is appropriate to interpret the fluvial landforms and the topography of lake basin.

TanDEM-X mission is designed to generate a global DEM with unprecedented accuracy [13] as the basis for various applications, such as landslide tracking [14], glacier changing [15–17], vegetation evaluating [18, 19], lava flow estimating [20, 21], and permafrost monitoring [22]. The TanDEM-X mission consists of twin satellites TerraSAR-X and TanDEM-X working on bistatic interferometry mode, which is characterized by the illumination of a target region on the ground by one transmitter and the simultaneous acquisition of the backscattered signals with two receivers [23]. This simultaneous acquisition avoids possible influences from temporal decorrelation and atmospheric disturbances [13]. The relative height accuracy is 2 m when slope is less than 20% [24]. Moreover, TanDEM-X acquisitions in StripMap mode have a ground spatial resolution of 2–3 m, which is much better than past SAR sensors previously cited, e.g., SIR-C (12.5 m), TOPSAR (10 m), ERS-1/2 (30 m), or ALOS/PALSAR (10 m). The high spatial resolution and the bistatic mode of the TanDEM-X mission therefore provide unprecedented opportunities to produce high-resolution DEM to characterize the elevation varying of Lop Nur playa.

This study aims to (1) eliminate remaining system errors with the least-square adjustment approach; (2) derive DEM of Lop Nur playa with a spatial resolution of 10 m using TanDEM-X/TerraSAR-X interferometric synthetic aperture radar (InSAR) image pairs; (3) compare the DEM results with ICESat GLAS footprints, aerial Structure-from-Motion (SfM) DEM, and SRTM DEM.

42 TanDEM-X/TerraSAR-X Coregistered Single look Slant range Complex (CoSSC) [25] image pairs from the TanDEM-X mission, which is supported by DLR TanDEM-X science coordination project of OTHER6906 and ATI_HYDR7333, are used to generate the target DEM. Table 1 lists the characteristics of image pairs used in this paper. They are all acquired on ascending orbit with right looking direction, and in HH polarization. The data are obtained during 2011 to 2015.

List of TanDEM-X COSSC used in this paper.

The preadjustment DEMs are processed in InSAR TanDEM-X Bistatic DEM Workflow module of SARscape© 5.2.1, where Goldstein filter [26] and Delaunay Minimum Cost Flow (MCF) [27] are assigned for noise reduction and phase unwrapping, respectively. A detailed description of the full processing of SAR data with SARscape© is provided by Sahraoui et al. [28]. The licenses of SARscape© and ENVI© are granted by Institute of Remote Sensing and Digital Earth, CAS. As the geo-location accuracy can be assumed as at least 10 m [29], the output resolution of geocoded DEM could be approximatively equivalent to it. Thus, the DEM was obtained with a ground spatial posting of 0.3 arc sec, which, in the studied area, corresponds to 7 m in longitude and 9 m in latitude. From this value and according to the range/azimuth resolution of the radar images, a multilooking factor of 4 is applied in range and in azimuth to the initial radar images to reduce the speckle. In the interferometric processing, elevation is only calculated for pixels that have a coherence value larger than 0.25 to prevent aberrant values in steep slopes where layover and shadowing effects are present.

Although instrument and system calibrations have been carried out to produce CoSSC [30], offsets of few meters and tilts in range and azimuth of some decimeters still remained [31], which means these errors should be estimated before mosaic.

The data of ICESat mission was used as the height truth value for both adjustment and validation. The GLAS on ICESat is a laser altimeter system that measures the Earth surface with a footprint of about 70 m. The footprints were successfully used as ground control points to generate DEM [23, 32, 33] for its reliable vertical accuracy of average 13 cm [34–37]. It operated from 2003 to 2009. GLA14 is the global land surface altimetry product of ICESat mission [34]. In this paper, to remove the noise in GLA14, selection rules considering cloud, attitude quality mark, energy saturation parameters, and the waveform characteristic [38] were executed to screen data. In particular, the waveform characteristic thresholds guarantee that the terrain variation in the footprint would be less than 1 meter [37]; thus the GLAS GCPs could be adaptive to InSAR DEM although their spatial resolutions are different. The data ellipsoid was transformed from TOPEX/Poseidon to WGS84 [39]. The ellipsoidal height was amended to orthometric height by geoidal undulation derived from the Earth Gravitational Model EGM2008 [40, 41]. 18,440 qualified GLAS GCPs were obtained after preprocessing, where half of them were randomly chosen as adjustment control points and the remaining were taken as the vertical height uncertainty check points.

Aerial photogrammetry SfM [42–45] is a burgeoning technology to obtain high-resolution DEM in a small scale. To evaluate the local performance of the resulting InSAR DEM, the ruined fort LK (89.667987°E, 40.094285°N) [46] near Lop Nur playa is taken as a validation site. In 2018, aerial photogrammetry SfM field work was conducted to obtain a 3D model of the ancient fort. Aerial photos were taken by DJI Phantom 3 Pro in 30 meters above land, with more than 80% overlay of an approximately 200 x 200 m area covering the fort ruin. Total of 456 photos were aligned to generate DEM within spatial resolution of 2 centimeters in PhotoScan©. A detailed PhotoScan© workflow was described by Johnson et al. [47]. The educational license of PhotoScan© is granted by Agisoft.

Since there are systematic errors in TanDEM-X DEM, majorly resulting from baseline estimation error, a polynomial function model [31] was used for least square adjustment, as follows:(1)gp,x,y=a0+a1·x+a3·ywhere g is the error function model; x is the position along the azimuth direction; y is the position along the ground range direction; p: p=a0,a1,b1, unknown error parameters; a0 denotes offset; a1 denotes tilt in azimuth; a3 denotes tilt in range.

The data process workflow was shown in Figure 1. Due to the sparse distribution of GLA14 (Figure 2) in Lop Nur playa, how to select the primary DEMs for adjustment is a problem. The primary DEMs are picked up manually according to the conditions: (a) there are at least 2 tracks of GLA14 in the DEM scene; (b) the GLA14 GCPs distribute throughout the DEM region as X shape, V shape, or parallel lines (Figure 2). GLA14 GCPs can be the height truth value for these selected DEMs, while the rest of DEMs depend on the tie-points extracted from the adjusted DEMs. Only part of the GLA14 GCPs was valid for adjustment, and the rest of them were used for height validation. Tie-points were extracted manually from the overlap region of adjacent adjusted DEMs. In order to reduce the influence of geo-location bias between DEMs, tie-points were selected on relatively flat terrain. Since the salt crust protects the playa surface from wind erosion, the terrain was assumed to be stable during the TanDEM-X acquisition period.

PHOTO (COLOR): The workflow of DEM generation.

PHOTO (COLOR): The primary DEMs, dependent DEMs, and screened GLA14 footprints. Images with minimum and maximum coherence were marked by different color.

A block adjustment approach was applied to generate global TanDEM-X DEM [23]. Since the number of DEM scenes was only 42, a simplified least-square adjustment approach was used in this paper. For the primary DEMs, the observation equation is as follows:(2)∆H^GLA14,I=gIp^,x,ywhere I is the index of DEM; ∆H^GLA14,I is the adjusted height difference between GLA14 GCPs and primary DEM.

For the dependent DEMs, the observation equation would be as follows:(3)∆H^J,K=gKp^,x,y-gJp,x,ywhere J is the index of adjusted DEM; K is the index of dependent DEM; ∆H^J,K is the adjusted height difference between dependent DEM and adjusted DEM.

Combining with GLA14 GCPs and tie-points, the estimated offsets and tilts of every DEM can be derived successively. To decrease the residual error between adjacent adjusted DEMs, averaging on the overlap region was applied in mosaic. 42 adjusted TanDEM-X DEM scenes were mosaicked to produce the Lop Nur playa DEM.

The coherence parameter is an indicator of the stability of the ground targets. Generally, it is a function of geometric and temporal decorrelation between the two acquisitions, and it ranges from 0 (unstable) to 1 (highly stable). In the case of bistatic interferometry, as temporal decorrelation is reduced, coherence is mainly dependent on the SAR hardware performance. Average coherence of item 1025014_018 is minimum while that of item 1080470_002 is maximum, whose coherence and interferogram were illustrated in Figure 3.

PHOTO (COLOR): Coherence of (a) item 1080470_002, (b) item 1025014_018, and interferogram of (c) item 1080470_002 and (d) item 1025014_018.

The result TanDEM-X DEM is presented in Figure 4. Although differences between adjacent DEMs were reduced by averaging, steps in the mosaicked DEM occurred partly. The residual error histogram (Figure 5) shows that the errors distribute more intensively in the mean value of ~0 than the normal distribution. The global Root Mean Square Error (RMSE) of the DEM is 0.57 m, which is in the range of TanDEM-X DEM accuracy specification of 2 m when slope is less than 20% [13]. The InSAR DEM difference to GLA14 was projected on the tinting DEM to present the spatial features (Figure 4). The orthometric height of the resulting Lop Nur DEM is from 785 m above sea level (asl) to 900 m asl, which covers an area of around 35,000 km2.

PHOTO (COLOR): TanDEM-X DEM and the error distribution. The size of each error point is in direct proportion to the absolute value of its error, and its color represents its magnitude.

PHOTO (COLOR): Histogram of TanDEM-X DEM height error to screened GLA14 footprints, with normal distribution fit line in red.

On one hand, it would be appropriate to compare the global performance of TanDEM-X DEM and SRTM DEM in the study region. SRTM was the previous generation global interferometric DEM [48, 49], whose absolute height error was estimated at several meters with a spatial resolution of 3 arc sec (~90 m) in Eurasia [50]. In Figure 6, SRTM was rendered with the same color ramp with Figure 4, where the height error to GLA14 GCPs was also projected. The size of each error point was in direct proportion to the absolute value of its error with the same ratio of Figure 4. The error histogram of SRTM to GLA14 GCPs was shown in Figure 7. The mean value was 3.04 m and the RMSE was 1.74 m.

PHOTO (COLOR): SRTM DEM and the error distribution. The size of each error point is in direct proportion to the absolute value of its error with the same ratio of Figure 4, and its color represents its magnitude.

PHOTO (COLOR): Histogram of SRTM DEM height error to screened GLA14 footprints, with normal distribution fit line in red.

On the other hand, it is proper to validate the local performance of TanDEM-X DEM referencing to the aerial DEM. The aerial SfM DEM is in a spatial resolution of 2 cm with relative elevation accuracy estimated at the order of decimeter [44]. Five points in LK (location see Figure 6) were selected for quantitative comparisons, and the results were summarized in Table 2. Since they are in different elevation systems, it makes sense to compare the relative height rather than the absolute height. The flat ground to the east of the ruined fort is taken as the local base height (Figure 6).

Relative elevation differences of the two DEMs.

Favorable coherence is crucial for interferometry. Item 10804700_002 is in the north of the study area (Figure 2). The landform is mainly mountain in the north of this image and solid salt crust in the south. Although the terrain varies in mountains (Figure 3(c)), the coherence is significantly high. Coherence in some parts of the mountains is relatively low at about 0.5 due to shadowing effects. The coherence of salt crust is also high. Item 1025014_018 is in the east of our region (Figure 2). The landforms are mountains and salt crust in the north of this image just like item 10804700_002, where the coherence is also near 1. In the north of item 1025014_018, the landform is desert. Due to the volume scattering of sand, the volume decorrelation effects decreased the coherence dramatically [51]. Similar volume decorrelation effects occur in other images, such as items 1130615_004 and 1043725_011, which makes the average coherence being lower than the others (see Table 1). However, in the western side of sand dunes which face the satellites, the coherence is relatively high at approximate or greater than 0.5. The local reliable features ensure interferogram available for further terrain inversion (Figure 3(d)).

The phase unwrapping method depends on the coherence features in the study region. There are three optional phase unwrapping methods in SARscape© 5.2.1, which are Region Growing, Minimum Cost Flow, and Delaunay Minimum Cost Flow, respectively.

(1) Region Growing is suggested to avoid setting a high coherence threshold in order to limit the possibility of introducing erroneous phase jump, like the unwrapping islands, in the output unwrapped phase image.

(2) Minimum Cost Flow should be adopted when the unwrapping process becomes difficult due to the presence of large areas of low coherence or other growing limiting factors; in such cases the Minimum Cost Flow algorithm enables obtaining better results than using the Region Growing method. This approach considers a square grid all over the image pixels. All pixels whose coherence is lower than the unwrapping coherence threshold are masked out.

(3) Delaunay Minimum Cost Flow is the same approach of the previous method, with the only difference that the grid does not necessarily cover all image pixels, but only those above the unwrapping coherence threshold. Moreover, it adopts the Delaunay triangular grid instead of square one. As a result, only the points with good coherence are unwrapped, without any influence from the low coherence pixels. It is especially useful when there are several areas of low coherence distributed throughout the image; in such case the other unwrapping approaches would eventually produce phase islands/jumps, while the Delaunay approach is able to minimize these jumps, which is the reason why Delaunay Minimum Cost Flow was chosen as the phase unwrapping method in this study.

However, the phase unwrapping error, including other error sources, still existed in resultant DEM. Referring to GLA14, errors in the eastern playa were small, usually less than 1 m, while the errors in the marginal desert, mountain, and Yardan region were sometimes significant. One of the possible reasons was that the geo-location of them was not rigorously aligned, or spatial scale of TanDEM-X DEM (10 m) did not match that of GLA14 (~70 m), so the elevation difference would be dramatic if the local terrain varies sharply, such as in the mountain and Yardan region. Besides, the volume scatter in desert could result in decoherence, which reduced the accuracy of phase unwrapping. Previous research showed that the high-frequency noise-like phase errors were usually below 1 m [31]. The RMSE of 0.57 m between resultant DEM and GLA14 indicated that the remnant errors are smaller than 1 m. On all accounts, the DEM was accurate enough to tint with a 1-m color ramp interval.

The visual performance of TanDEM-X DEM (Figure 4) was greater than SRTM DEM (Figure 6). There were some obvious northwest-southeast artifacts in SRTM DEM. Due to the insufficient spatial resolution and height accuracy, SRTM DEM presented the rougher terrain variations than TanDEM-X. More precisely, the spatial resolution of TanDEM-X DEM 0.3 arc sec was one order of magnitude higher than that of SRTM DEM 3 arc sec, and the standard deviation of the error reduced by a factor of 3, with a value of 0.57 m for TanDEM-X (Figure 5), against 1.74 for SRTM (Figure 7). In Lop Nur playa, TanDEM-X DEM was more reliable than previous SRTM DEM.

Locally, the TanDEM-X DEM presented roughly the topography features of the ruined ramparts (Figure 8(a)). The relative heights of three points on the tops of the ramparts were all underestimated, which were higher than the assigned base height, and the difference was up to -2.8 m at the SW Wall (Table 2). A 20-m (2-pixel) displacement at the SE Wall was observed, which could result from the layover effect of the SAR system. The relative elevation of the two depressions inside the ruin was also underrated, even the depression characteristics were hardly recognized. Regardless, the worst relative error in such complicated terrain could be evaluated as 3 m, as shown in Table 2. The difference might be ascribed to the limitations of the radar system. On one hand, the 10-m resolution of InSAR DEM limits its spatial availability compared to SfM DEM. On the other hand, both walls and depressions are isolated islets with striking slope on their margin, where the interferometric fringes would discontinue easily on an interferogram [52]. Hence, it requires us to be cautious when interpreting local isolated landforms. In any case, the relative height accuracy of 3 m was in the range of TanDEM-X DEM specification (4 m when slope>20%) [13]. Although TanDEM-X DEM was not as accurate as aerial SfM DEM in both horizontal and vertical dimensions, its pivotal superiority was the full coverage of Lop Nur playa with a relatively reliable precision. It could be hardly to achieve an aerial SfM DEM throughout the study region in consideration of the amount of work. Besides, it was also difficult to find homologous points for SfM in homogeneous landforms (salt crust, desert, etc.). In current stage, spaceborne InSAR was the most proper way to derive a full-coverage DEM for Lop Nur playa.

PHOTO (COLOR): Comparison of DEM from different resources. (a) TanDEM-X DEM; (b) aerial SfM DEM. Color ramp is the same as that in Figure 3. Both are projected on WGS84 UTM Zone 45N.

The latest high-resolution TanDEM-X DEM could contribute to the paleo-environment and anthropology researches around Lop Nor playa. Previously, it was believed that there were only two main inflowing riverways for lake Lop Nor [1]. One was along the north of the Tarim Basin (region A in Figure 9), while the other was in the southeast of lake (region B in Figure 9). A wide riverway (region C in Figure 9) in the middle of the known riverways could be interpreted from the resultant DEM. Especially in the upstream of this riverway, a depression (region D in Figure 9) to the southeast of Loulan ruin could be identified, which was the Loulan depression relating with the evolution of ancient Lop Nur lake [2]. There were several hydrological channels in Loulan depression. Region C and Loulan depression might be used to be a near-lake wetland and affect the daily lives of residents in Loulan city, according to its geographical location. Geometry society paid little attention to these areas due to insufficient topographic data before. Besides, the altitude of the ancient ruins, channels and basins, could also be an important clue to investigate the relationship between human beings and nature. A key strategy was that the lake surface level could not be higher than the altitude of the ancient ruins when these ruin cities were prosperous. It might be a novel evidence for understanding sediment profiles in Lop Nur.

PHOTO (COLOR): Interpretation of the TanDEM-X DEM. The known ruin ancient forts, LE, Loulan, and LK, are labeled in the figure. River channels are marked with blue lines.

This paper describes the derivation of a highly accurate DEM from TanDEM-X CoSSC images of Lop Nur playa via a least-square adjustment approach.

The resultant DEM covers an area of around 35,000 km2 and has a horizontal resolution of 0.3 arc sec (~10 m) and a height accuracy of 0.57 m, representing the most accurate DEM in Lop Nur playa until now. The horizontal spatial resolution of the derived DEM is one magnitude higher than that of SRTM DEM, and the height accuracy is improved by a factor of 3. Validation of the derived DEM demonstrates its reliability for complicated terrains though it is not as good as for the flat region.

This highest quality DEM derived for Lop Nur can provide a basic material for general paleo-environment researches. A newfound river mouth of Lop Nur near Loulan ruin is interpreted from the resultant DEM. More filed work in this area will be expected to reveal the relations between anthropology and nature.

GLA14 data used to support the findings of this study have been deposited in the National Snow & Ice Data Center website (https://nsidc.org/data/icesat/). TerraSAR-X/TanDEM-X CoSSC and derived DEM used to support the findings of this study were supplied by DLR under license and so cannot be made freely available. Requests for access to these data should be made to TanDEM-X Science Coordination, tandemscience@dlr.de. SRTM DEMs used to support the validations of this study have been deposited in CGIAR-CSI SRTM 90m DEM Digital Elevation Database (http://srtm.csi.cgiar.org/srtmdata/). The SfM DEM data used to support the findings of this study are available from the corresponding author upon request.

The authors declare that there are no conflicts of interest regarding the publication of this paper.

This work was supported by the National Natural Science Foundation of China [41671353, U1303285, 41571363, 41431174, 61471358, 41201346, 41301394, 41301464] and the CAS Knowledge Innovation Program [KZCX2-EW-320].