High-Resolution Spaceborne SAR Geolocation Accuracy Analysis and Error Correction

Abstract

High-accuracy geolocation is crucial for high-resolution spaceborne SAR images. Most advanced SAR satellites have a theoretical geolocation accuracy better than 1 m, but this may be unrealizable with less accurate external data, such as atmospheric parameters and ground elevations. To investigate the actual SAR geolocation accuracy in common applications, we analyze the properties of different geolocation errors, propose a geolocation procedure, and conduct experiments on TerraSAR-X images and a pair of Tianhui-2 images. The results show that based on GNSS elevations, the geolocation accuracy is better than 1 m for TerraSAR-X and 2 m/4 m for the Tianhui-2 reference/secondary satellites. Based on the WorldDEM and the SRTM, additional geolocation errors of 2 m and 4 m are introduced, respectively. By comparing the effectiveness of different tropospheric correction methods, we find that the GACOS mapping method has advantages in terms of resolution and computational efficiency. We conclude that tropospheric errors and ground elevation errors are the primary factors influencing geolocation accuracy, and the key to improving accuracy is to use higher-accuracy DEMs. Additionally, we propose and validate a geolocation model for the Tianhui-2 secondary satellite.

1. Introduction

Spaceborne synthetic aperture radar (SAR) is an indispensable tool in Earth observation. SAR satellites produce images by actively transmitting and receiving microwave signals and are, therefore, able to work all day and in all weather. For a specific kind of study, when the quality of the high-resolution optical images is unsatisfactory, SAR images are almost the only available data source. When both types of images are available, SAR images are still helpful in providing distinctive radiation information, improving the results [1,2,3,4,5], which is a common benefit in landcover classification research [6,7,8].

As remote sensing technology develops, an increasing number of spaceborne satellites with excellent performance have been launched, providing a wealth of high-resolution data for research and applications. For example, the advanced TerraSAR-X and Tianhui-2 satellites offer stripmap products with a resolution better than 5 m. At this resolution level, the outlines of most ground targets are recognizable, enabling SAR images to be applied in more domains, such as target change monitoring and urban feature extraction [9]. High geolocation accuracy is crucial for utilizing the geometric information of high-resolution SAR images, enabling them to correct the geometry of optical images [10,11] and complete surveying tasks [12]. When matching the geometric features from high-geolocation-accuracy SAR images and other data, no additional processing besides geolocating is required because images are naturally registered after geolocation. Otherwise, the SAR imaging geometry must be corrected first [13,14], increasing the complexity of the data processing, especially when control point data are unavailable. In interferometric SAR (InSAR) research, geolocation inaccuracy can be compensated for using coregistration methods like the amplitude correlation method, but different types of images cannot be properly registered using these methods [15,16,17].

The factors affecting SAR geolocation accuracy are solid Earth effects, orbit error, atmospheric error, and ground elevation error [18]. The influence of solid Earth effects is relatively minor, so we only discuss this in Section 4. The orbital accuracy of advanced SAR satellites has exceeded that of early satellites by a large margin, benefitting from improvements in SAR instruments and orbit observation accuracy, thus ensuring a high upper limit for SAR geolocation accuracy. The orbit error of TerraSAR-X is at the centimeter level [18]. For Sentinel-1, Envisat, and COSMO-SkyMed, it is less than 1 m [18,19]. Therefore, we do not treat orbit error as a primary error source and instead focus on atmospheric error and ground elevation error. The findings in the Global Navigation Satellite System (GNSS) domain help correct atmospheric error, including tropospheric delay and ionospheric delay. For SAR signals, the properties of tropospheric delay are nearly the same as those of GNSS signals, so tropospheric correction algorithms and products can be applied to SAR images directly, including the ray-tracing (RT) method [20] and zenith path delay (ZPD) data [21]. Correcting ionospheric delay also benefits from GNSS products, mainly depending on the vertical total electron content (VTEC) data. To achieve the best possible correction effect, Zhang et al. proposed a relatively comprehensive model based on VTEC data [22]. Given that the ground elevation error is more closely related to the quality of the elevation data than the SAR data, researchers have conducted fewer studies on it compared with other errors. Ground elevation error is expected to introduce a geolocation error proportional to itself [23,24], but no experiments have been conducted to verify this relationship, and no typical values are available for reference.

Although the theory of SAR geolocation error correction was well developed many years ago, auxiliary data for correction have recently improved relatively quickly. For tropospheric correction, the following two new datasets are worth noting: the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 [25] and the Generic Atmospheric Correction Online Service for InSAR (GACOS) [26,27,28]. The accuracy and resolution of ERA5 have reached new heights compared with previous atmospheric reanalysis data, and the GACOS attempts to simplify tropospheric correction using the zenith delay mapping (ZDM) method. In terms of elevation data, the high-quality WorldDEM achieves much higher vertical accuracy than the widely used Shuttle Radar Topography Mission (SRTM). However, the two tropospheric datasets have not yet been introduced to SAR geolocation, and it is unknown how much improvement will be achieved by replacing the SRTM with the WorldDEM.

To fully utilize the aforementioned advanced auxiliary data and simplify the procedure while ensuring geolocation accuracy, we geolocated measurement points on TerraSAR-X and Tianhui-2 images. For the tropospheric correction, we did not use ERA5 data directly but instead input it into the Weather Research and Forecasting Model (WRF) to obtain tropospheric parameters with higher resolution. We calculated the delay using both the WRF results and the GACOS data. We studied the influence of the ground elevation error by comparing the geolocation result based on GNSS-measured height and the height data from the WorldDEM and the SRTM. By conducting experiments on Tianhui-2 images, we found that the traditional range–Doppler model is not applicable, so we proposed a bistatic model for Tianhui-2. After removing the influence of the orbit error, the Tianhui-2 geolocation results showed patterns consistent with those of TerraSAR-X, proving the effectiveness of our model and the simplified geolocation procedure. Our research demonstrates the advantages of the GACOS and the WorldDEM and can serve as a reference for other high-resolution SAR geolocation studies.

This article is organized as follows: Section 2 analyzes the principles of the geolocation errors and introduces the geolocation procedure. Section 3 presents the variation ranges of the errors and compares the effectiveness of different correction methods. Section 4 discusses deficiencies in our experiment and other factors affecting geolocation accuracy. Section 5 contains the conclusions of this research.

2. Materials and Methods

2.1. Components of Geolocation Errors

The components of geolocation errors differ in terms of azimuth and range. Typically, the azimuth geolocation accuracy is only affected by the satellite orbit error owing to inaccuracies in satellite positions or mismatched instrumental times [18]. Coregistrating SAR images with digital elevation model (DEM)-simulated amplitude images is a widely used orbit correction method without ground control points [23]. However, it does not apply to high-resolution images and meter-level-accuracy geolocation. We discuss this method in Section 4. In this study, we assumed that there was no significant azimuth geolocation error.

In addition to the orbit error, the range geolocation accuracy is also affected by the atmospheric error and ground elevation error (Figure 1) [18,23]. Influenced by the atmosphere, the signal’s equivalent velocity, $v$, is smaller than the velocity of light in a vacuum, $c$. The propagation time will decrease by $\left({\displaystyle \frac{1}{v}}-{\displaystyle \frac{1}{c}}\right)\left|S{T}_2\right|$ when calculated using $c$, called atmospheric delay [18]. For convenience, atmospheric delay can be approximated and expressed in a distance form as $\mathsf{\Delta }{s}_{atm}={\displaystyle \frac{c-v}{c}}\left|S{T}_2\right|$. The atmospheric delay comprises the ionospheric delay and the tropospheric delay. The ionospheric delay primarily depends on the ionospheric state and the signal’s frequency [22]. In extreme conditions, ionospheric delay is at the decimeter level for X-band SAR, about 1 m for C-band SAR, and may exceed 10 m for L-band SAR [18]. Tropospheric delay is mainly influenced by the tropospheric state and the length of the propagation path. Given that the signal’s zenith angle, $\theta$, generally ranges from 20° to 70°, and the zenith tropospheric delay is typically 2 m to 3 m, the total tropospheric delay is about 2 m to 4 m.

The elevation error, $\mathsf{\Delta }h$, causes a geolocation error, $\mathsf{\Delta }{s}_{ele}=\mathsf{\Delta }h\mathit{cos}\mathsf{\theta }$ (Figure 1). Since the accuracies of the currently available DEMs are limited [29], the elevation error tends to have the largest impact on the geolocation accuracy, potentially introducing an error of approximately 10 m.

Table 1. SAR geolocation errors in the range.

Geolocation Errors	Typical Value	Factors
$\mathrm{Tropospheric} \mathrm{delay}, \mathsf{\Delta }{s}_{tro}$	2 m~4 m [30]	Atmospheric state, zenith angle
$\mathrm{Ionospheric} \mathrm{delay}, \mathsf{\Delta }{s}_{ion}$	0.1 m~20 m [22]	Electronic density, signal frequency
$\mathrm{Elevation}-\mathrm{related} \mathrm{error}, \mathsf{\Delta }{s}_{ele}$	0.01 m~10 m [29]	Elevation accuracy

lists the components of the range geolocation errors.

Based on the previous discussion, let $∆{\mathrm{s}}_{\mathrm{t}\mathrm{r}\mathrm{o}}$ represent the tropospheric delay, $∆{\mathrm{s}}_{\mathrm{i}\mathrm{o}\mathrm{n}}$ denote the ionospheric delay, and $∆{\mathrm{s}}_{\mathrm{e}\mathrm{l}\mathrm{e}}$ signify the elevation-related error. We categorized the geolocation error caused by other factors, including the satellite orbit error, into the following residual terms: $∆{\mathrm{r}}_{\mathrm{r}\mathrm{g}}$ in range and $∆{\mathrm{r}}_{\mathrm{a}\mathrm{z}}$ in azimuth. Let $\mathsf{\Delta }{\mathrm{s}}_{\mathrm{r}\mathrm{g}}$ and $\mathsf{\Delta }{\mathrm{s}}_{\mathrm{a}\mathrm{z}}$ denote the geolocation errors in range and azimuth. We obtain the following:

$$∆{\mathrm{s}}_{\mathrm{r}\mathrm{g}}=∆{\mathrm{s}}_{\mathrm{t}\mathrm{r}\mathrm{o}}+∆{\mathrm{s}}_{\mathrm{i}\mathrm{o}\mathrm{n}}+∆{\mathrm{s}}_{\mathrm{e}\mathrm{l}\mathrm{e}}+∆{\mathrm{r}}_{\mathrm{r}\mathrm{g}}$$

(1)

$$∆{\mathrm{s}}_{\mathrm{a}\mathrm{z}}=∆{\mathrm{r}}_{\mathrm{a}\mathrm{z}}$$

(2)

Next, we introduce the tropospheric and ionospheric correction methods and the geolocation procedure.

2.2. Tropospheric Correction

The most common atmospheric correction methods are the RT method and the ZDM method [20]. The RT method models the actual signal propagation process and, thus, possesses superior accuracy [31]. The ZDM method approximates the propagation process and has advantages in computational complexity. Many InSAR studies have demonstrated the effectiveness of the ZDM method for phase correction [32]. We validate its efficacy for geolocation in this study.

The RT method divides the propagation path into multiple segments, calculates the tropospheric delay corresponding to each segment, and sums them (Figure 2). Since the result is virtually unchanged when segments are more than 20 [18], dividing the path into 20 segments is sufficient to meet the accuracy requirement for geolocation. The bending in the propagation path due to refraction has nearly no impact on the results, so approximating the path as a straight line introduces only negligible errors [22].

Implementing the RT method requires the pressure, $P_{i }(hPa)$; temperature, $T_{i }(K)$; vapor partial pressure, $e_i(hPa)$; and path length, $h_i\left(m\right),$ for each segment, as shown in Equations (3) and (4) [33,34], where $\theta$ is the zenith angle, and $N$ is the radio refractivity, equaling ${10}^{-6}\left(n-1\right)$, and $n$ is the refractivity.

$$N_i=77.6{\displaystyle \frac{P_i}{T_i}}+71.6{\displaystyle \frac{e_i}{T_i}}+3.75\times {10}^5{\displaystyle \frac{e_i}{T_i^2}}$$

(3)

$$\mathsf{\Delta }{s}_{tro}={\displaystyle \frac{{10}^{-6}}{\cos \theta }}\sum _i{N}_i·h_i+{\displaystyle \frac{{10}^{-6}}{\cos \theta }}77.6\times 287.05\times {\displaystyle \frac{P_{top}}{9.81}}$$

(4)

We added an extra term to the conventional formula in Equation (4) to correct the delay above the tropopause. The vapor partial pressure is virtually zero; hence, the additional delay can be approximated as hydrostatic delay, which can be estimated using the pressure at the tropopause.

The ZDM method estimates the tropospheric delay in the propagation path using the zenith delay. The distance between the zenith path and the propagation path of the signal increases as the altitude increases. However, the thickness of the troposphere is limited, generally not exceeding 30 km [22]. Therefore, the maximum distance is around 30 km at the tropopause, at which the scales of the tropospheric parameters do not significantly differ. The rationality of replacing the parameters along the propagation path with those along the zenith path is the theoretical foundation of the ZDM method. The mapping function is needed to transform the length of the zenith path into that along the propagation path, which can be $\sec \theta$ or from related products, such as VMF1. The most significant advantage of the ZDM method is simplifying the computation, allowing the zenith tropospheric delay to be potentially released as a standard product, such as the GACOS and other station observation products.

2.3. Ionospheric Correction

A common practice for ionospheric correction is to approximate the entire ionosphere as a thin layer located at a height of 450 km, referring to the intersection point of the signal, with the thin layer as the piercing point. The ionospheric delay is calculated using the vertical total electron content (VTEC) at the piercing point (Equation (5)) [22].

$$\left|∆s_{ion}\right|={\displaystyle \frac{40.31\times {10}^{16}·VTEC}{f_0^2\sqrt{1-{\left(\frac{R}{R+H}\mathit{sin}\mathsf{\theta }\right)}^2}}}$$

(5)

The unit of the VTEC is the total electron content unit (TECU), equaling ${\displaystyle \frac{{10}^{16}}{m^2}};$ $f_0$ (Hz) represents the signal frequency; R is the radius of the Earth; H stands for the height of the ionosphere 450 km; and $\mathsf{\theta }$ denotes the zenith angle on the ground. As most SAR satellites orbit within the ionosphere, the result is larger than the actual ionospheric delay, so theoretically, an adjustment based on the satellite’s orbital altitude is needed [18]. Gisinger et al. proposed that the coefficient is 75% for TerraSAR-X and 90% for Sentinel-1 [18], while Zhang et al. considered the coefficient to be around 70% for Sentinel-1 [22]. The adjustment contributes to centimeter-level improvements for Sentinel-1 but may worsen the ionospheric correction accuracy for ALOS-2 [22]. Most of the time, the adjustment is insignificant—less than 5 cm for TerraSAR-X based on our calculations. Therefore, we do not discuss the adjustment for the ionospheric delay further.

2.4. Geolocation Procedure

A geolocation procedure consists of the following two parts: coordinate system transformation and error correction [18]. The transformation from the geographic coordinate system into the SAR image coordinate system is called indirect geolocation, while the reverse process is known as direct geolocation [35,36]. Given that our research can be applied to target monitoring and data registration, we focus solely on indirect geolocation. We used the range–Doppler model [37] to perform the transformation. This study does not discuss other models, such as the rational polynomial model, which is based on photogrammetric principles with relatively lower accuracy [38]. The only error we correct in our procedure is atmospheric. We omit orbit correction and correction for solid Earth effects, which might be applied in other studies to reduce the complexity of the SAR geolocation.

Assuming that one of the geographic coordinates to be geolocated has been converted into the geocentric coordinate ${\overrightarrow{R}}_T$, the range vector, $\overrightarrow{R}(t)={\overrightarrow{R}}_T-{\overrightarrow{R}}_S(t)$, can be obtained with the satellite’s position, ${\overrightarrow{R}}_S(t)$, at the imaging time, $t$. The fundamental idea behind geolocation is to determine $t$ and $|\overrightarrow{R}(t)|$, correct the tropospheric delay and ionospheric delay for $|\overrightarrow{R}(t)|$, calculate the image coordinates (row and col), and resample the image (Figure 3).

The imaging time, $t,$ is solved by the Doppler equation (Equation (6)), where ${\overrightarrow{V}}_S(t)$ represents the velocity of the satellite, $f_0$ denotes the signal frequency, and $c$ is the speed of light in a vacuum.

$$f_D=-2{\displaystyle \frac{f_0}{c}}{\displaystyle \frac{{\overrightarrow{V}}_S(t)·\overrightarrow{R}(t)}{\left|\overrightarrow{R}(t)\right|}}$$

(6)

Most SAR images are produced in zero-Doppler geometry, $f_D\equiv 0$, so ${\overrightarrow{V}}_S(t)$ is orthogonal to $\overrightarrow{R}(t)$. We solve for $t$ using an iterative method (Equation (7)).

$$t_{k+1}=t_k+{\displaystyle \frac{{\overrightarrow{V}}_S\left(t_k\right)·\overrightarrow{R}\left(t_k\right)}{{\left|{\overrightarrow{V}}_S\left(t_k\right)\right|}^2}}$$

(7)

The principle is as follows (Figure 4a): ${\displaystyle \frac{{\overrightarrow{V}}_S\cdot \overrightarrow{R}}{\left|{\overrightarrow{V}}_S\right|}}$ is the projection of $|\overrightarrow{R}|$ in the direction of ${\overrightarrow{V}}_S$. If ${\overrightarrow{V}}_S$ is orthogonal to $\overrightarrow{R}$, the projection should be zero. We need to increase the imaging time, $t$, to meet the condition, with the increment value being ${\displaystyle \frac{{\overrightarrow{V}}_S\cdot \overrightarrow{R}}{{\left|{\overrightarrow{V}}_S\right|}^2}}$.

Since the satellite keeps moving between the signal-transmitting and -receiving moments, Equation (6) is only an approximation, but studies have shown that it is sufficiently accurate [39]. With $t_{start}$ representing the azimuth start time of the image, $dt$ denoting the azimuth line time, $R_0$ representing the nearest range of the image, $dR$ representing the range pixel spacing, and $∆r_{tro}$ and $∆r_{ion}$ representing the tropospheric delay and ionospheric delay, the image coordinates $(row \mathrm{a}\mathrm{n}\mathrm{d} col)$ can be calculated by Equations (8) and (9). The signs of $∆r_{tro}$ and $\mathsf{\Delta }{r}_{ion}$ are positive, so the range coordinate, $col$, is relatively smaller without atmospheric correction.

$$row={\displaystyle \frac{t-t_{start}}{dt}}$$

(8)

$$col={\displaystyle \frac{\left|\overrightarrow{R}(t)\right|-R_0+∆r_{tro}+\Delta r_{ion}}{dR}}$$

(9)

Owing to the difference in the image standards, the azimuth coordinates of the Tianhui-2 images were calculated using Equation (10) [40].

$$row={\displaystyle \frac{t-t_{start}+\frac{\left|\overrightarrow{R}(t)\right|}{c}}{dt}}$$

(10)

For the Tianhui-2 secondary images, the single-satellite model in Equation (6) is not applicable. The range vector, ${\overrightarrow{R}}^{\prime }$, and the velocity, ${\overrightarrow{V}}_S^{\prime }$, of the reference satellite must be considered, as shown in Equations (11)–(13) and Figure 4b, where $\mathsf{\Delta }t={\displaystyle \frac{\left|\overrightarrow{R}\left(t\right)\right|}{c}}$. The azimuth coordinate, $row$, follows the same formula as the reference images, as shown in Equation (10) [40].

$${\displaystyle \frac{{\overrightarrow{V}}_S\left(t+\mathsf{\Delta }t\right)·\overrightarrow{R}\left(t+\mathsf{\Delta }t\right)}{\left|\overrightarrow{R}\left(t+\mathsf{\Delta }t\right)\right|}}+{\displaystyle \frac{{\overrightarrow{V}}_S^{\prime }\left(t-\mathsf{\Delta }t\right)·{\overrightarrow{R}}^{\prime }\left(t-\mathsf{\Delta }t\right)}{\left|{\overrightarrow{R}}^{\prime }\left(t-\mathsf{\Delta }t\right)\right|}}=0$$

(11)

$$t_{k+1}=t_k+{\displaystyle \frac{{\overrightarrow{V}}_S\left(t_k+\Delta t\right)·\overrightarrow{R}\left(t_k+\mathsf{\Delta }t\right)}{2{\left|{\overrightarrow{V}}_S\left(t_k+\mathsf{\Delta }t\right)\right|}^2}}+{\displaystyle \frac{{\overrightarrow{V}}_S^{\prime }\left(t_k-\mathsf{\Delta }t\right)·{\overrightarrow{R}}^{\prime }\left(t_k-\mathsf{\Delta }t\right)}{2{\left|{\overrightarrow{V}}_S^{\prime }\left(t_k-\Delta t\right)\right|}^2}}$$

(12)

$$col={\displaystyle \frac{\left|\overrightarrow{R}\left(t\right)\right|+\left|{\overrightarrow{R}}^{\prime }\left(t\right)\right|}{2dR}}+{\displaystyle \frac{-R_0+∆r_{tro}+\Delta r_{ion}}{dR}}$$

(13)

The final step is the sampling of the pixel (row and col) and placing it in the original geographic coordinate.

The error, $\mathsf{\Delta }h$, in the elevation will lead to a geolocation error, $\mathsf{\Delta }{s}_{ele}=\mathsf{\Delta }h\mathit{cos}\theta$, in the range (Figure 5). If the elevation with error is more significant than the actual elevation, the range coordinate $col$ with error will be relatively minor. We discuss the range geolocation accuracy in the image coordinate system for the rest of the article; so here, we use $\mathsf{\Delta }\mathrm{h}\cos \mathsf{\theta }$ to describe the range error. The distance $|OB|$ on the ground is $\mathsf{\Delta }\mathrm{h}\cot \mathsf{\theta }$, which is more commonly seen in direct geolocation research.

2.5. Materials and Experiments

Our study area is located within the Changping District of Beijing, China, spanning latitudes 40.15°N to 40.27°N and longitudes 116.20°E to 116.35°E, covering an area of 10 km × 10 km. The geomorphological types from northeast to southwest are a mountainous region, a piedmont alluvial fan, and a plain. The measurement points were obtained by GNSS differential measurement with centimeter-level 3D positioning accuracy. A total of 164 points were utilized (Figure 6), of which 135 are situated in flat locations within the plain or piedmont alluvial fan areas, and 29 are in the mountainous region.

This study used 11 TerraSAR-X images and a pair of Tianhui-2 images. The TerraSAR-X satellite constellation operates in the X-band at a frequency of 9.65 GHz. It comprises two identical satellites, TerraSAR-X and TanDEM-X, which exhibit nearly the same performance. The two satellites were launched in June 2007 and June 2010. The TerraSAR-X images we used were acquired by the TanDEM-X satellite, from 3 July 2018 to 29 May 2019, and are Level 1B single-look complex (SLC), slant-range, and HH-polarized images in the stripmap mode. The ground resolution is azimuth 3.3 m × range 2.2 m, the pixel size is azimuth 1.8 m × range 0.9 m, and the swath is azimuth 50 km × range 30 km.

The Tianhui-2 satellite constellation, similar to TerraSAR-X, comprises two identical satellites that operate in the X-band at a frequency of 9.60 GHz. The reference and secondary satellites acquired the selected images on October 28th, 2019, and are Level 1B SLC, slant-range, and HH-polarized images. They share the same image parameters. The ground resolution is 3.0 m in both azimuth and range, the pixel size is azimuth 2.0 m × range 0.9 m, and the swath is azimuth 30 km × range 30 km. The orbit accuracy of the Tianhui-2 satellites is lower than that of TerraSAR-X, with potential errors reaching up to 2 m [41,42].

We first verified the completeness of the geolocation error decomposition in Equations (1) and (2). In Section 3.1, we geolocate the measurement points based on the GNSS elevations, correct the tropospheric delay through the RT method, implement the ionospheric correction, and evaluate the results. When performing tropospheric correction, we utilize ERA5 data and the WRF program. ERA5 is a reanalysis dataset for the troposphere derived from global meteorological observations [25], dividing the troposphere into 37 layers according to the pressure, providing the pressure, temperature, and vapor partial pressure for each layer, with a horizontal resolution of 0.25° × 0.25° and a temporal resolution of 1 h. The WRF is a high-accuracy numerical weather model widely applied in InSAR research [43]. We input the ERA5 data as the initial condition into the WRF to obtain atmospheric parameters at the imaging time with a resolution of 0.01° × 0.01°. For the ionospheric correction, we utilize the TEC products from the International GNSS Service (IGS). The TEC product, recorded in units of 0.1 TECU (${10}^{15}/{\mathrm{m}}^2$), comes with a spatial resolution of 2.5° × 5° and a temporal resolution of 2 h [44].

Section 3.2 checks the validity of ERA5, analyzes the characteristics of the tropospheric delay, and compares the performances of the RT method and the ZDM method. Owing to a lack of in situ data, we verify ERA5 using the ZPD data from the BJFS station near the study area. The GACOS is an emerging source of ZPD data; it is quasi-real-time and globally acquirable with a high resolution of 90 m × 90 m. In this section, it proves to perform as well as ERA5. In Section 3.3, we calculate the upper limit and variation ranges of the ionospheric delay for the whole year.

Section 3.4 analyzes the impact of the DEM elevation error on the geolocation accuracy. The SRTM and the WorldDEM are used in this section. The SRTM has a resolution of 30 m × 30 m; its vertical datum reference is the EGM96 geoid model, and raw data were collected between 11 and 22 February 2000. The WorldDEM has a higher resolution of 12 m [45]; its vertical datum reference is the EGM2008 geoid model, and raw data were acquired in January 2011 and ceased around mid-2015.

3. Results

3.1. Geolocation Accuracy

We chose the image acquired on 17 January 2019, as the reference image for the TerraSAR-X. First, we visually selected the measurement points in the TerraSAR-X reference image and then coregistered it to other TerraSAR-X images to control the quality of the point selection. For each measurement point (Figure 7a), we selected a 65 m × 65 m area surrounding it, subdividing it into 1 m × 1 m grids, and transferring the image onto this area based on the GNSS elevation. Given that the elevation-related geolocation error, $\mathsf{\Delta }{s}_{ele}$, is at the centimeter level, the range geolocation error, $\mathsf{\Delta }{s}_{rg}$, is approximately equivalent to the atmospheric delay, $\mathsf{\Delta }{s}_{tro}+\mathsf{\Delta }{s}_{ion}$. The geolocation results without atmospheric correction have smaller range coordinates (Figure 7b), consistent with the model in Equation (9). We checked the geolocation results with the atmospheric correction for all 164 measurement points and found that every point was accurately geolocated to the expected coordinate (Figure 7c), indicating that the geolocation residuals, $∆s_{rg}$ and $∆s_{az}$, were less than $1$ m. The results substantiate the completeness of the geolocation error decomposition in Equations (1) and (2) and the effectiveness of the atmospheric correction. Through the geolocation procedure proposed in Section 2.4, the TerraSAR-X images can achieve an accuracy better than $1$ m.

We coregistered the TerraSAR-X images using version 2.6.3 of the Interferometric Synthetic Aperture Radar Scientific Computing Environment (ISCE) software. The ISCE resampled the secondary image to the coordinate system of the reference image before coregistration to mitigate terrain effects, which is equivalent to geolocation without atmospheric correction (Figure 7b) but in the radar center coordinate system. The ISCE uses a constant offset in each direction to describe the coregistration result (Figure 8). In azimuth, the offsets are all less than 0.1 pixels, indicating virtually no geolocation error. In range, the offsets are less than 0.2 pixels, except for the image acquired on 5 August 2018, for which the tropospheric delay is about 0.3 pixels larger than those of the reference image. Therefore, we infer that the main cause of the offset is the difference in the tropospheric delay. The offsets affected by other factors, including the difference in ionospheric delays, are usually less than 0.2 pixels.

By contrast, the Tianhui-2 image geolocation results exhibit noticeable residuals. The geolocation error of the reference image is approximately $-1$ m in azimuth and $-1.5$ m in range (Figure 9a). For the secondary image, the azimuth geolocation error is $-1$ m, too, while the range error is estimated to be −2.5 m (Figure 9c). Because the geolocation error aligns with those of other studies [40], we infer that the primary source of the geolocation residual is the orbit error of the Tianhui-2 satellites. The secondary image’s lower geolocation accuracy may be attributable to deficiencies in time synchronization between the two satellites. We shifted the reference image by −1 m in azimuth and −1.5 m in range to correct the results and performed a similar correction for the secondary image (Figure 9b,d). We found that both images could geolocate all of the measured points accurately, indicating that the residual is independent of the imaging time, range, and Doppler frequency. The results also substantiate the effectiveness of the geolocation model for the secondary satellite described in Section 2.4.

3.2. Characteristics of the Tropospheric Delay

We calculated the zenith delay using ERA5 and compared the result with the ZPD in the BJFS station (Figure 10). The difference is larger in summer than in winter, possibly because of higher vapor pressure. The maximum difference is about 0.06 m, corresponding to less than 0.1 pixels, so we believe that ERA5’s accuracy is sufficient for geolocation. The difference is not correlated with the delay value (Figure 10b).

The tropospheric delay In the study area is calculated using the RT method and ranges from $2.6$ m to $2.8$ m, demonstrating a strong correlation with the ground elevation (Figure 11a,b). The delay in the plain varies on a smaller scale, from $2.75$ m to $2.8$ m, while in the mountainous region, it is typically 2.6 m to $2.75$ m. This is because the signal must traverse a longer distance before it reaches lower altitude areas, and there tends to be more tropospheric perturbations in mountainous regions.

We calculated the zenith tropospheric delay in the study area using the same tropospheric parameters as the RT method and then mapped it toward the line-of-sight (LOS) direction. The difference between it and the RT result is less than 0.3 cm (Figure 11c). We examined the relationship between the radio refractivity, $N$, and elevations less than 1200 m (Figure 11d). Assuming N is 300, only 10 m is needed to form a delay of 0.3 cm, as per Equation (4). Considering that the accuracy of the elevation inside ERA5 is usually less than 10 m, the difference between the two methods is of little reference value. This distribution of N supports this deduction (Figure 11d). At the same elevation, the average fluctuation of N is at most five when the elevation is less than 600 and becomes less than one after that, corresponding to a difference of only 0.3 cm. Compared with the other method, the path length accuracy is more important. Therefore, we believe that the two methods perform equally well. The results suggest that the ZDM method generally meets the accuracy requirements for geolocation.

We compared the results of the GACOS-based ZDM method and the ERA5-based ZDM method (Figure 12). The GACOS-based result exhibits the same characteristics as the ERA5-based result; that is, the tropospheric delay was about 5 cm larger in the east plain than in the west plain (Figure 12a). The difference exhibits a noticeable texture associated with the terrain (Figure 12b). It tends to be positive in valleys and negative at ridges, which is related to the WRF-processed ERA5’s 1 km lower resolution than the 90 m for GACOS. The discrepancy is about −1 cm on the plain, whereas it barely reaches extremes of around −12 cm in mountainous regions. The difference between the two data sources is much larger than between the two methods. However, owing to its advantages in resolution, although we have verified the high precision of ERA5, we still believe that the GACOS results are at least as good. Compared with the total geolocation error in range, we believe that all three methods demonstrate enough tropospheric accuracy.

The relationship between the ground elevation and tropospheric delay is linear, although the degree is notably weaker in the mountainous region than for the plain (Figure 13a). A linear regression analysis shows that for every 100 m increase in elevation within the study area, the tropospheric delay decreases by approximately 2.9 cm. In most cases, the residual of the fitting is less than 5 cm (Figure 13b). Given that the resolution of the raw ERA5 data is only 0.25° × 0.25°, comparable in scale to the extent of the study area, we only obtain a single value if tropospheric correction is directly implemented based on the ERA5 data, leading to a residual of 10 cm ~ 20 cm in the tropospheric correction (Figure 13a). Therefore, the accuracy of the tropospheric correction utilizing the raw ERA5 data is inferior to that achieved using the GACOS.

Owing to the continuity of the tropospheric parameters and the strong correlation between the tropospheric delay and ground elevation, the accuracies of the different tropospheric methods are close. Considering that the GACOS has the highest resolution, ensuring it performs well even under complex terrain conditions and is easy to access and use, we propose that the GACOS-based ZDM method is relatively optimal.

3.3. Characteristics of the Ionospheric Delay

Given the significant impact of the ionospheric delay on the X-band satellites, including TerraSAR-X and Tianhui-2, we primarily focus on the numerical characteristics of the ionospheric delay. Given that SAR satellites are typically in dawn–dusk or noon–midnight orbits, we collected the VTEC of the study area at 12 PM and 6 PM local time in 2019 (Figure 14a), finding that the variation ranges of the two moments are approximately the same, with a maximum value of around 20 TECU. Owing to direct sunlight in the tropics, the average VTEC in that area is more significant than in the study area (Figure 14b), reaching up to 25 TECU.

We used the SAR satellites TerraSAR-X, Sentinel-1, and ALOS-2 across three signal bands as examples and calculated the ionospheric delay when the VTEC was 25 TECU and θ was 45° ((Equation (5);

Table 2. Parameters associated with ionospheric delay in the three bands.

Satellite	Band	Frequency	Ionospheric Delay Corresponding to 25 TECU	VTEC Leads to 1 m of the Delay
TerraSAR-X	X	9.65 GHz	0.19 m	134.6 TECU
Sentinel-1	C	5.41 GHz	0.59 m	42.3 TECU
ALOS-2	L	1.26 GHz	10.89 m	2.3 TECU

). Generally, the ionospheric delay for X-band satellites is comparable in magnitude to the residual of the tropospheric correction and, thus, minimally impacts the geolocation accuracy, rendering ionospheric correction unnecessary. The variation in the ionospheric delay is even smaller (discussed in Section 3.1). For the C-band, the ionospheric delay is usually less than 1 m. The accuracy of the ionospheric correction is acceptable even if a residual remains. However, in the L-band, among all factors, the ionospheric delay is the most significant influence on geolocation, potentially causing an error as high as 10 m. Employing higher-resolution data such as Jet Propulsion Laboratory High-Resolution Global Ionospheric Maps (JHR) enhances the positioning precision significantly [22].

3.4. Characteristics of the Elevation-Related Error

We calculated the elevation-related error by subtracting the GNSS-based result from the DEM-based result. When geolocation was implemented on the WorldDEM and the SRTM, the elevation-related error, $\mathsf{\Delta }{s}_{ele},$ corresponding to the peak probability density, is negative in both cases (Figure 15a), indicating a tendency for the DEM elevation to be larger than the GNSS elevation at the measurement points. Specifically, the peak-density error is −1.5 m for the SRTM, and the absolute error is less than 4.0 m for 83% of the measured points. For the WorldDEM, the peak-density error is −0.6 m, and the absolute error is less than 2.0 m for 86% of the points (Figure 15b). Considering that the atmospheric correction residual is typically less than 1 m and often around 0.1 m, we conclude that the range geolocation error, $∆s_{rg}$, can generally be controlled to within 5 m when using the SRTM and to within 3 m when using the WorldDEM. There is no relatively effective method for correcting the elevation-related error, so the influence of this error can be considered the most severe among all contributing factors.

We analyzed the distribution of $\mathsf{\Delta }{s}_{ele}$ and found that it exhibits a significant correlation with the slope of the measurement points (Figure 16). As the slope increases, the distribution of $\Delta s_{ele}$ becomes less focused, and the average gradually drifts away from the origin point. For the points with the same slope intervals, $\mathsf{\Delta }{s}_{ele}$ is smaller overall for the WorldDEM than the SRTM. By comparing the optical images at the time the SAR images were acquired and when the DEMs were obtained, we discover that the terrain around the points with the largest $\mathsf{\Delta }{s}_{ele}$ has undergone significant changes, indicating that utilizing the DEM close to the time the image is obtained can help prevent extreme conditions.

4. Discussion

4.1. Generalizability of the Results

Limited by the availability of data, we only conducted our experiment in one area. The properties of the tropospheric error and the ground elevation error may be different in other places. However, we can predict the performance of our procedure based on relevant results from the GACOS and SRTM.

Yu et al. calculated the ZPD in eight areas with different topographic variations using the GACOS [26]. They considered the ZPD retrieved using the InSAR method as the benchmark, obtaining a root mean square (RMS) of around 1 cm and a maximum residual of around 6 cm, consistent with our results after transforming the ZPD into tropospheric delay (Figure 12b). Therefore, we believe that the quality of the GACOS is stable in most parts of the world. The remaining issue is the stability of the ZDM method. The maximum elevation in our study area is about 700 m (Figure 11a), while the distribution of N steadies when the elevation is larger than 600 m (Figure 11d). The difference between the ZDM method and the RT method lies in the fluctuation of the radio refractivity, N. Assuming that N converges when the elevation is greater than the maximum ground elevation, the fluctuation of N remains five when the elevation is less than the maximum, the same as in our research. When the maximum ground elevation is 3000 m, a tropospheric correction residual of only 1.5 cm is introduced by the ZDM method. Therefore, we believe that the GACOS-based ZDM method applies in most areas. The residual may be larger when the tropospheric condition is more complex. Further experiments are needed.

The degree of terrain flatness significantly affects the quality of DEMs and eventually affects the geolocation accuracy. Since the related geolocation error is proportional to the DEM error, we can predict the performances of the WorldDEM and SRTM in other areas by referring to studies on DEM quality. Han et al. studied two DEMs in four areas with different slopes [29]. The terrain of their study area in Inner Mongolia, China, is similar to our study area, where 90% of the WorldDEM errors are less than 2 m, and 90% of the SRTM errors are less than 4 m. They inspected the DEM accuracy in a flatter study area in Xinjiang, China, and a more complex area in Sichuan, China. For the WorldDEM and SRTM, the errors were 0.8 m and 5 m in Xinjiang and 12 m and 16 m in Sichuan. The performance of the WorldDEM is always better than the SRTM, but the advantage may not be significant in areas with complex terrain. There are no effective algorithms to achieve high-accuracy geolocation without high-quality elevation data.

Theoretically, our geolocation procedure applies to other satellite images, especially high-resolution images, such as those from Sentinel-1 and ALOS-2. However, it is hard to recognize our measurement points in Sentinel-1 images, so we did not conduct experiments on them. When processing C-band and L-band images, a more highly accurate ionospheric model may be needed [22]. For medium- and low-resolution satellite images, the focus of geolocation processing differs from that of our research. In those cases, geolocation errors related to atmospheric delays and ground elevation errors are usually less than 10 m, while orbit errors can be greater than 10 m [23]. Therefore, the effectiveness of our procedure will be insignificant without orbit correction.

4.2. Limitations of Orbit Correction Based on DEM-Simulated Image Coregistration

The SRTM used to be a reference for estimating geolocation errors in SAR images owing to its high horizontal accuracy, and a common method was to register SAR images and DEM-simulated amplitude images. We tested the effectiveness of this method using the ISCE (Figure 17).

We projected the SRTM onto the SAR image coordinate system using the Topo module and refer to the result as z data, which the Simamplitude module can further utilize to generate a simulated amplitude image with the same horizontal accuracy as the SRTM. Finally, the Ampcor module calculated the range shift between the two images.

Given the substantial differences between the SAR image and the simulated image, the Ampcor result generally has a much lower signal-to-noise ratio (SNR) than the results from the two real SAR images. The SNR corresponding to the peak of the probability density is 1.3 for the TerraSAR-X image (Figure 18a) and 2.9 for the Tianhui-2 reference image (Figure 18b), with most of them weaker than 5.0, which is a common threshold for SAR image coregistration. Therefore, we do not believe that this method applies to high-resolution SAR geolocation.

4.3. Other Geolocation Error Sources

Solid Earth effects and approximations in image production may also introduce geolocation errors in addition to satellite orbit error, atmospheric error, and elevation-related error. Among solid Earth effects, solid Earth tides (SETs) have the most significant impact on geolocation accuracy. These are induced by the tidal forces of the Moon and Sun [46], potentially causing elevation changes of up to half a meter in equatorial regions and up to 40 cm in mid-latitude areas. Apart from SETs, ocean tidal loading and atmospheric tidal loading can also introduce a centimeter-level error [47,48]. Commonly, studies employ programs written by Milbert to eliminate these effects [49]. Owing to the relatively minor influence of the solid Earth effects and the fact that they are hardly considered when processing other data, we have not discussed them.

Approximations applied to SAR image production can lead to geolocation errors. For Sentinel-1 images, the bistatic effect, Doppler shift in range, and Doppler FM-rate mismatch are three typical error sources [18,50]. The bistatic effect arises when the “stop-and-go” model assumes that the sensor receives all signals corresponding to the same LOS direction simultaneously. However, they are not simultaneous in reality, introducing an error of about 2 m to 3 m in azimuth. Doppler shift results from neglecting the variation in the Doppler frequency matching the same LOS, leading to a range error as large as 0.4 m. Doppler FM-rate mismatch affects the azimuth geolocation accuracy, and the error is typically less than 1.5 m. We have not discussed these factors because these approximations are not applied to TerraSAR-X images [18].

5. Conclusions

SAR satellites can acquire data all day and in all weather conditions to periodically provide images for research areas. The radiometric features of SAR images make them highly complementary to optical images and other data. With improvements in SAR image resolution, progressively finer geometric features make SAR applicable to an increasing number of domains. High-accuracy geolocation is crucial for high-resolution SAR images to ensure correct feature matching. This study analyzed the properties of SAR geolocation errors and examined correction methods, proposing an effective procedure for SAR geolocation. We employed a TerraSAR-X image and a pair of Tianhui-2 images to validate the geolocation procedure, with the results demonstrating the following:

• The theoretical SAR geolocation accuracy meets the requirement for common SAR applications. The TerraSAR-X images in this research achieved a geolocation accuracy better than 1 m based on the GNSS elevations. Affected by the orbital error and other factors, the geolocation accuracy values of the Tianhui-2 reference and secondary images are approximately 2 m and 4 m, respectively. Considering that most advanced SAR satellites have orbit data that are more accurate than 1 m, the TerraSAR-X result is of greater reference value.
• The accuracy of the elevation data primarily constrains the SAR geolocation accuracy, and using high-accuracy elevation data is critical to enhancing geolocation accuracy. When geolocating the TerraSAR-X image using the SRTM, the accuracy was better than 5 m for most measurement points. Upon switching to the higher-accuracy WorldDEM, the geolocation error for most measurement points fell below 3 m.
• The GACOS mapping method is the most optimal tropospheric correction method, because the ZDM method provides accuracy comparable to that of the RT method. At the same time, the GACOS has advantages in resolution and ease of use.

The Tianhui-2 secondary satellite receives signals transmitted by the reference satellite. The secondary-image geolocation results validate the effectiveness of the secondary-image geolocation model proposed in this research. Possibly because of the inaccuracy of the time synchronization, the geolocation accuracy of the secondary image was inferior to that of the reference image.

The difficulty in obtaining high-accuracy DEMs is the principal bottleneck in high-accuracy SAR geolocation. In the future, it might be possible to eliminate dependence on DEMs through 3D geolocation methods and to enhance the applicability of SAR geolocation.