In-Depth Verification of Sentinel-1 and TerraSAR-X Geolocation Accuracy Using the Australian Corner Reflector Array

This article shows how the array of corner reflectors (CRs) in Queensland, Australia, together with highly accurate geodetic synthetic aperture radar (SAR) techniques—also called imaging geodesy—can be used to measure the absolute and relative geometric fidelity of SAR missions. We describe, in detail, the end-to-end methodology and apply it to TerraSAR-X Stripmap (SM) and ScanSAR (SC) data and to Sentinel-1 interferometric wide swath (IW) data. Geometric distortions within images that are caused by commonly used SAR processor approximations are explained, and we show how to correct them during postprocessing. Our results, supported by the analysis of 140 images across the different SAR modes and using the 40 reflectors of the array, confirm our methodology and achieve the limits predicted by theory for both Sentinel-1 and TerraSAR-X. After our corrections, the Sentinel-1 residual errors are 6 cm in range and 26 cm in azimuth, including all error sources. The findings are confirmed by the mutual independent processing carried out at University of Zurich (UZH) and German Aerospace Center (DLR). This represents an improvement of the geolocation accuracy by approximately a factor of four in range and a factor of two in azimuth compared with the standard Sentinel-1 products. The TerraSAR-X results are even better. The achieved geolocation accuracy now approaches that of the global navigation satellite system (GNSS)-based survey of the CRs positions, which highlights the potential of the end-to-end SAR methodology for imaging geodesy.

T HE ability to accurately geolocate an SAR image onto a predefined Earth model (e.g., a reference ellipsoid or a Digital Terrain Model) without needing ground control points is a unique feature of the SAR imaging technique known since the late 1970s and the SEASAT mission [1]. The applied computational strategy and the discussions on the determined shortcomings given in [1] are still of high importance for any of today's SAR missions: accurate orbit determination; a tight link between orbit ephemeris time and the SAR instrument; accurate handling of the timing by the SAR instrument and in the subsequent image processing; usage of the range-Doppler model for geolocation in range and azimuth; agreement between the physical SAR imaging process and the approximations made in the SAR processor; and high quality of the reference used to verify the geolocation. Other considerable absolute geolocation error sources, which were assumed to be negligible at the time of the SEASAT mission, are the meter-level path delays introduced by the Earth's atmosphere and, at a submeter scale, the dynamic effects of the solid Earth.
In this article, we present a comprehensive description of all these elements in order to achieve centimeter-level geolocation and solve the common limitations in the SAR processing of contemporary SAR missions. Assessing the SAR ALE relies on the comparison of range and azimuth coordinates of point targets extracted from the images with the expected values derived from the range-Doppler model and the known target coordinates. Current spaceborne SAR missions design the geometric accuracy to meet the resolution supported by the SAR sensor. They typically aim for an ALE of half the pixel spacing to ensure accurate geolocation of the SAR scenes for direct comparison with other geospatial data sets. For Sentinel-1, RADARSAT-2, TerraSAR-X, and COSMO-SkyMed, this translates into official requirements for the guaranteed geolocation ranging from 1 m up to a few meters [2]- [6].
However, rigorous geometric considerations in all of the computations and meticulous correction of perturbing signals enable SAR geolocation (or SAR range and azimuth measurements) far below the average resolution of the image product. This advanced approach was first studied in detail for TerraSAR-X high-resolution spotlight products using long-term stable CRs [7], [8], ultimately resulting in a one sigma geolocation accuracy of 1-2 cm for suitable targets (e.g., CRs) [9]. Meanwhile, other SAR missions have partly adopted these advanced methods. In particular, for the Sentinel-1 mission and its openly accessible image products [10], there is now a growing demand to improve the geolocation accuracy of the SAR products beyond the initially demonstrated performance [11].
During the commissioning phase of the S1A spacecraft launched in April 2014, the SM product was determined to be well within the 2.5-m geolocation specification, but the remaining azimuth offsets of about 1.8 m were larger than expected [11]. Such behavior was also detected later in the IW product and was found to increase with slant range [12]. Systematic effects were also identified in the range  [17]. results of the Sentinel-1 IW product [12], [13]. These findings were related to known bistatic SAR effects stemming from platform motion during transmission and reception of the radar signal [14]- [16]. However, up until now, a consistent explanation of their connection with the processing of the Sentinel-1 IPF SAR processor was not available. Therefore, another motivation of this article is to provide methods to correct the geometric processing of Sentinel-1 and achieve an ALE for the standard IW product that is ultimately only limited by the quality of the atmospheric corrections and the orbit, as is already the case with TerraSAR-X.
To analyze and evaluate the new methods with Sentinel-1, we use the permanent CR array located in Queensland, Australia (see Fig. 1). It consists of 40 triangular trihedral CRs distributed over a large geographic area, with each CR having an inner leg dimension of between 1.5 and 2.5 m. For each CR, accurate reference coordinates are known from GNSS surveys [18]. We first prove the integrity of this CR array with TerraSAR-X SM and SC data, where we have extensive prior experience [9]. In summary, the goals of this article are: 1) a thorough description of the SAR geolocation verification based on the practical experience of DLR/TUM and UZH; 2) documentation of impacts of common approximations applied in SAR processing with a special emphasis on Sentinel-1; and 3) demonstrating the methods with TerraSAR-X and Sentinel-1 using the Australian CR array, thereby assessing quality and possibilities of this unique SAR ground infrastructure.
This article is organized as follows. Section II describes the method for verifying the SAR image geolocation with CRs, providing recommendations for installation of new CRs and a step-by-step procedure to meticulously derive the ALE of SAR products. Section III focuses on the postprocessing of Sentinel-1 timing values to account for the platform motion effects and the topographic mismatch, currently limiting the geolocation accuracy. The section is supplemented by an appendix, which describes the details of the correction  Table I); example shows a 0.6-m CR in K-band (λ = 1.25 cm). Images adapted from [14], [Fig. 7.12] and [20], [Fig. 3].
computations based on the annotations of Sentinel-1 IW products. In Section IV, we describe the Australian CR array, and Section V gives an overview of the Sentinel-1 and TerraSAR-X data sets used in this article. Section VI presents the geolocation results of the study; first, the verification of the array with TerraSAR-X, and, second, the detailed analysis of Sentinel-1 IW product, which involved cross-validations of both UZH and DLR/TUM methodologies. In Section VII, we discuss the results in a wider context, and Section VIII contains the final conclusions and outlook to future developments.

A. Corner Reflector Installations
Triangular trihedral CRs [see Fig. 2(a)] are often considered the most practical device for calibrating radar systems: they have a well-defined geometric reflection center and a high peak RCS [see Fig. 2(b)] and are inherently stable and relatively simple to manufacture and inexpensive compared with active devices. The RCS signature does not decrease very quickly off-boresight [see Fig. 2(b)], making it useful in situations where precise alignment of the CR boresight with the SAR sensor boresight cannot be achieved or the use of multiple pass geometries is envisaged [19]. Most of the backscatter from a trihedral CR (i.e., corresponding to a 3-dB beamwidth) correspond to a range spanning nearly 40 • in both elevation and azimuth off-boresight directions [see Fig. 2(b)]. A disadvantage is that to enable a large RCS signal, it is necessary to use large CRs, especially at lower SAR frequencies, such as L-band. The peak theoretical RCS (that occurs along the boresight vector) σ peak of triangular trihedral CRs depends on the wavelength λ and the size of the reflector, given by the inner side length (inner leg dimension) a [20] σ triangle,peak = 4πa 4 3λ 2 .
In order to detect the CR in the image to determine its radar coordinates, a minimum visibility of the CR against the clutter background in a real SAR image is required. The contrast affects the ability to accurately measure the CR's peak intensity position in an SAR image and is usually expressed by the SCR, which is given by the ratio of the CR intensity I peak (i.e., radar cross section divided by pixel area) and the mean background intensity I mean clutter surrounding it [14] SCR = I peak I mean clutter .
(2) Table II lists typical SCR values (in decibels) inferred from the Sentinel-1 and TerraSAR-X product types and CR sizes used in this article. The SCR can be used to estimate the range and azimuth measurement accuracy provided by the SAR product (i.e., the limit of detecting the peak of the CR point response in the image) [21], [22] where ρ R,A denotes the product resolution in slant range and azimuth, respectively. Table II gives values of σ R,A across the different Sentinel-1 and TerraSAR-X products. We expect some accuracy limitations in the azimuth direction for the Sentinel-1 IW product and the TerraSAR-X SC product due to the coarser resolution. For slant range as well as the TerraSAR-X SM product, in general, the uncertainties of the orbit, the correction for the atmospheric path delays, and the accuracy of the reference coordinates of the Australian Array are expected to dominate the overall geolocation error. Given a CR of known dimensions and other acquisition constraints, the simplest way to guarantee maximum SCR is to place the CR in a location of low clutter (low backscatter intensity from the surroundings). Large concrete or asphalt surfaces without nearby structures are ideal, but fields with short grass generally also generate low enough backscatter to make high signal-to-clutter ratios (SCRs) possible even with reasonably sized CRs (e.g., [23]). The Australian CR array described in this article is a good example of the latter case. To ensure a "clean" impulse response function from the CR, no strong scatterers or their sidelobes should be visible within at least several samples in all directions from the planned CR location. This means that products with larger sample spacing (e.g., Sentinel-1 EW mode products) require proportionally larger low-clutter areas for ideal CR placement. In addition to CR placement considerations, optimizing the CR orientation relative to the expected imaging geometry (or geometries) is also a critical factor when designing measurement campaigns. Once installed, the CR should not be moved or reoriented as this would potentially affect the phase center position and disturb the stability of the geometric SAR time series. To orient a trihedral CR for a spaceborne mission, the azimuth (horizontal) angle should be aligned to the approximate heading angle of the orbital tracks into consideration and the elevation angle aligned to the mean satellite elevation. Usually, more than one acquisition geometry is needed to ensure geometric measurement diversity. In the case of the Australian CR array, the reflectors are all oriented toward a commonly used spaceborne SAR ascending-orbit geometry, with azimuth angles near 256 • and CR boresight elevation angles ranging between 53 • and 59 • across the array [18]. The inability to create perfect boresight alignment with the SAR viewing angles across different missions highlights the advantage of using trihedral CRs, which has a wide 3-dB RCS pattern (e.g., Fig. 2).
Given the high geometric accuracy and overall system stability of today's spaceborne SAR sensors, sensor calibration and validation experiments should include reference CRs with positions known to centimeter accuracy or better. Therefore, another potential source of error is the surveyed CR position itself. Measurement of the CR vertex position as defined by the intersection of the three orthogonal plates (corresponding to the brightest point in an SAR image, Fig. 2) using differential GNSS (DGNSS) provides positions with cm-level accuracy if the recordings are performed over a period of at least 20-30 min, as recommended in [24,  It is not possible to place a GNSS antenna at the position of the CR vertex directly, so a reference point position is surveyed instead. The reference point is usually a well-defined feature on the ground over which the CR is to be mounted. On a natural surface, such as a field, a fixed reference point may be created (e.g., by hammering a metal spike into the ground). In this situation, the CR vertex is placed directly above the surveyed mark, and the vertical offset is measured and considered in the GNSS postprocessing. An alternative approach, used in the case of the Australian CR array, involves determining the 3-D offset of the CR vertex from the surveyed reference point.
Finally, weather conditions can have a significant impact on the usefulness of the CRs as geolocation targets. Bird-transported debris can quickly clog up drainage holes, and subsequent rainfall may accumulate in the CR. During winter times, snow may accumulate in the reflector. Both of these effects are observable in SAR images as a decrease in both the SCR and RCS [18]. Monitoring one or both of these intensity indicators is highly recommended, especially for CRs that is not regularly inspected.

B. SAR Geolocation Analysis Procedure
The verification of the geolocation quality of spaceborne SAR systems is based on the comparison of the range and azimuth measurements extracted from the image and the expected values derived from the geometric imaging model. The imaging model describes how the 3-D coordinates of a CR are projected to 2-D image coordinates. It includes the satellite orbit product, the Earth's dynamic effects, the SAR sensor, the wave propagation, and the methods employed by the SAR processor. We recommend using the SLC images for the analysis because it is the first level of SAR processing that results in an image with slant-range and azimuth geometry. Moreover, to ensure the best possible geolocation analysis at the centimeter level, meticulous care must be taken in all the involved procedures. We have decomposed the process into a step-by-step sequence that consists of: 1) extraction of the SAR measurements; 2) correction of the atmospheric delays; 3) computation of the CR reference position considering the solid Earth dynamics; 4) prediction of the expected range and azimuth values; and 5) evaluation of the SAR ALE. The methods build upon our earlier research, which can be found in [7], [8], [25], and [26]. In the following, we provide descriptions that address the key aspects of each of the five steps.
1) Point Target Analysis: The determination of the CR's location in the SLC image is part of the PTA, which yields not only the measured range time τ m,CR and the measured azimuth t m,CR but also provide useful quality information, such as the SCR. From a theoretical point of view, the range and azimuth coordinates of a CR may be extracted down to less than 1/100th of a pixel, only limited by the SCR of the point response, as discussed in Section II-A. Due to the finite bandwidth of the transmitter and the finite amount of time the ground is in view of the radar beam, the SAR signal is band-limited in range and azimuth [14]; thus, the point response of a CR is spread across several pixels, resulting in a cross-shaped signature as shown in Fig. 3.
To achieve such a precise extraction in practice, we use spectral zero-padding to oversample the patch of 32 × 32 pixels centered at the peak that contains the CR point response. The 2-D spectrum computed by the fast Fourier transform is extended with zeros at the spectral minimum, resulting in an oversampling equivalent to the number zeros after applying the inverse fast Fourier transform. Large numbers of zeros, such as 1024 or 2048, can provide an oversampling to directly locate the peak with sufficient accuracy, but their usage is problematic when regarding the memory consumption and computational speed. Therefore, the combination of zero-padding using factors of 32 or 64 and refinement with an analytical surface is more efficient [21]. In our PTA approach, we fit an elliptic paraboloid of the shape to the central 3 × 3 pixels of the peak in the amplitude image (see Fig. 3). i and j denote the pixel indices, a 0,1,2,3 are the coefficients of the paraboloid, and m i and m j are the refined pixel indices of the peak location. From our tests comparing zero-padding with high oversampling factors against the two-step approach, we can confirm that the latter method is capable of maintaining the required accuracy of better than 1/100th of a pixel [25].
2) Atmospheric Path Delay Corrections: The group velocity of electromagnetic waves propagating in the Earth's atmosphere is lower than the speed of light in vacuum [27]. Therefore, the SAR signals traveling through the atmosphere experience a group delay. As a result, the inferred SAR ranges become too large if the speed of light in a vacuum is assumed to convert the measured two-way travel time to the geometric range. Space geodetic techniques, such as GNSS, also make use of radio signals and, thus, have developed methods to determine and correct the atmospheric delay effects [24], [27]. Following these GNSS-based methods, the atmosphere can be decomposed into two parts for radio signals of up to 30 GHz: a nondispersive neutral part (termed troposphere) and a dispersive part consisting of free electrons and ions (termed ionosphere). The typical impact of the tropospheric delay in the SAR slant range is in the order of 2.5-4 m, whereas the ionospheric delay may reach several decimeters for X-band (9.65 GHz, TerraSAR-X), up to 1 m for C-band (5.405 GHz, Sentinel-1), and several meters for L-band (1.2 GHz, ALOS-2).
For the mitigation of tropospheric path delays in SAR geolocation, several approaches have been proposed, e.g., [8], [28], and [26]. At geodetic stations and close to permanent GNSS stations, the tropospheric zenith delay inferred from GNSS observations may be used (see, for instance, the products issued by the IGS [29]). Moreover, the delay can be computed through integration of the atmospheric state described by 4-D models, such as the operational European Centre for Medium-Range Weather Forecasts (ECMWF) analysis model [30], or by model-based methods using data recorded at close-by terrestrial meteorological stations [31]. In the following, we describe the two tropospheric corrections that are applied in this article at the Australian CR array: the ECMWF analysis model and the use of data from meteorological stations.
A straightforward way to employ the operational ECMWF analysis data is the VMF1 product. 1 It defines not only the mapping function to convert between tropospheric zenith delays and tropospheric slant-range delays in geodetic techniques, such as like GNSS or Very Long Baseline Interferometry [27], [32]; the VMF1 product also provides reliable path delay estimates in the zenith direction, which are precomputed from the operational ECMWF analysis model along with the coefficients of the mapping function. The data are distributed as global grids with a spatial resolution of 2 • × 2.5 • (latitude by longitude) and four epochs per day, i.e., 00, 06, 12, and 18 h UT. These grids enable direct computation of both the hydrostatic and the wet component of the zenith tropospheric delay at a global scale when applying a spatiotemporal interpolation for the site of interest. The details on the usage of the delay data as well as the comparison of the standalone path delays with the best estimates from GNSS stations are given in [33]. The validation shows an agreement on the order of 1-2 cm for the tested stations [33] and confirms the applicability of the VMF1 product as a tropospheric correction method. Schematically written, the tropospheric path delay in slant range of an SAR observation at the location φ, λ, h (latitude, longitude, and height) and epoch t may be summarized as where z h,w are the hydrostatic and wet zenith delays inferred from the grids and corrected for the actual CR height [33], and VMF1 h,w denote the hydrostatic and wet mapping functions depending on the zenith angle ζ of the SAR sensor to CR LOS. Alternatively, the tropospheric path delay can be inferred from terrestrial meteorological station data that are freely available online. For this article, the used data with a temporal resolution of 15 min is provided by the Australian Bureau of Meteorology [34]. The temperature, pressure, and humidity measurements from stations near the observed site and close to the acquisition time are coupled with two models to determine the hydrostatic path delay and the wet path delay in the zenith direction. The hydrostatic delay is based on the well-established model of Saastamoinen [35] in the form described in [36], whereas the wet delay follows the closed-form expression of the path delay integration that is proposed in [37]. A comprehensive summary of the method, including the scaling of the meteorological parameters for the desired station height (i.e., the height of the CR), can be found in the UNB3 model of the University of New Brunswick [38], which formed the original basis of our implementation [31]. For the mapping to SAR LOS, the total delay is converted by a multilayer extension of the mapping functions (described, for instance, in [27]), which uses multiple slant range to layer height ratios to approximate the curvature of the Earth and the tropospheric layers. The concept is illustrated in Fig. 4 for three layers. We found that 20 layers suffice to provide less than 1-mm error (compared with using more layers), even for incident angles up to 55 • . Schematically written and in analogy to (5), the 20-layer method may be stated as where P, T , and e denote the air pressure, the temperature, and the water vapor pressure at the station, and R n over H n is the ratio of the slant-range path and scale height of the respective layer. Critically, the layer boundaries are defined such that each contributes the same fraction of the total modeled delay. Layer pierce points along the vertical and slant-range rays are determined by iterating the height-dependent path delay model while keeping track of the positions in global Cartesian (XYZ) coordinates. This leads to a layer stack with increasingly thick layers at higher altitudes (as shown in Fig. 4). The application of the approach has been demonstrated to provide tropospheric path delay estimates accurate to the centimeter level for a CR near the meteorological station, especially when care is taken during the conversion from zenith to LOS delay [11], [12], [31]. Fig. 4. Multilayered "ratio stack" approach to derive the atmospheric slant delay. The concept is shown here for three layers (the final computation uses 20 layers), whereby the atmosphere is modeled as a stack of constant-density layers. Each layer contributes the same amount to the total delay through a unique conversion ratio R n /H n , as shown in (6).
Regarding the ionospheric path delay contribution, the global ionospheric maps computed under the umbrella of the IGS are a reliable source for modeling and correcting the dispersive delay [39]. The maps are inferred from the global network of permanent IGS GNSS stations and provide an estimate for the free electrons in the upper atmosphere (the so-called TEC). The maps have a resolution of 5 • × 2.5 • (longitude, latitude) and a temporal resolution of 1 h and are distributed as daily data cubes. The TEC information is vertically condensed to a single spherical layer at an altitude of 450 km, which is the typical height of the maximum electron concentration [40]. In our analysis, the vTEC computation relies on the TEC solution of the CODE (University of Berne, Switzerland), which contains not only the TEC maps but also the corresponding uncertainty maps derived from the global least-squares approach applied to compute the daily solutions [40].
The vertical TEC can be converted to a frequency-dependent slant-range path delay [40] where f is the frequency, vTEC is the vertical total electron content interpolated at the location (latitude ϕ, longitude λ) of the IPP for the epoch t of the SAR measurement, and MF(ζ ) is the mapping function, depending on the SAR sensor to CR LOS. Note that the vTEC interpolation is not performed at the CR location but at the IPP, which is the intersection of the LOS vector with the modeled spherical shell containing the vTEC; the details are available in [40]. The vTEC, as stated in (7), accounts for the signal of the overall ionosphere, roughly extending as far out as 1500 km above the Earth's surface [27], whereas TerraSAR-X and Sentinel-1 have an orbit altitude of 514 km [5] and 712 km [3], respectively. From an analysis of the electron concentration with height computed from the international reference ionosphere model IRI2007 [41], we derived a 75 % rescaling factor for the total vTEC to match the altitude of TerraSAR-X, and we assume a factor of 90 % for the altitude of Sentinel-1. Here, we used these constant scaling factors, but further improvements to better account for spatial and temporal variations in the ionospheric signal could be made.
For the SAR analysis of the Australian CR array, we benefit from several Australian GNSS receivers contributing to the global IGS network. 2 Because of these stations, the expected uncertainty of the ionospheric path delay correction derived for the SAR observations is in the order of 1-2 cm. The values apply for the Sentinel-1 C-band frequency of 5.405 GHz, which experiences a stronger delay by the dispersive ionosphere than the 9.65-GHz X-band frequency used by TerraSAR-X. For TerraSAR-X, the estimated uncertainty is at the millimeter level.
In the conclusion of this section on correcting the atmospheric path delays, we can summarize our methods as follows: DLR/TUM uses the VMF1 product to derive the tropospheric path delays, whereas UZH relied on data from nearby meteorological stations. In terms of the ionosphere, DLR/TUM has fully implemented (7), while UZH approximated the ionospheric delay by evaluating the ionospheric maps at the CR location and performing a simplified cosine-based mapping [31]. Thus, a refined approach may be tested at UZH in the future that corresponds more closely to the ionospheric model at the heart of the TEC products. We further address the impact of the different path delay computation strategies used by TUM/DLR and UZH on SAR geolocation in Section VII.
3) Computation of CR Reference Coordinates: Geodetic measurements of the Earth's surface performed by satellites are sensitive to the dynamics of the solid Earth (e.g., GNSS [24]). SAR satellites are no exception when analyzing range and azimuth with respect to the given orbit, and therefore, the state of the Earth's crust must be modeled to describe the reference position of a CR at the epoch of the SAR acquisition. The precise orbit determination of TerraSAR-X and Sentinel-1 is made possible with the onboard GNSS and is given in the ITRF (see [42] for the latest release ITRF2014) when using the GNSS products of the IGS [43], [44]. Consequently, the orbit state vectors annotated to the SAR products of TerraSAR-X and Sentinel-1 comply with the ITRF and its conventions, which allows us to use the same conventions to correctly model the CR position in the ITRF.
The models used within the ITRF are documented in the geodetic conventions issued by the IERS and comprise all the crust displacements related to tidal dynamics. Table III lists all the effects considered, along with their typical orders of magnitude. The direct tides of the solid Earth are the dominant effect that causes a daily position variation on the order of centimeters horizontally and decimeters vertically. Moreover, there are secondary crustal effects due to the tides in the ocean (ocean tidal loading) and in the atmosphere (atmospheric tidal loading), as well as the rotational deformations linked to polar motion (pole tides and ocean pole tide loading). The details and the usage of the underlying geodynamic models can be found in [45]. We have implemented all of these effects according to the conventions, which enables the definition of  [45] the instantaneous CR position as  Table III) may be inferred from the ITRF solution of the closest IGS GNSS station located on the same tectonic plate [47].

4) Prediction of Reference Timings:
With the ability to compute the instantaneous CR position at the epoch of the SAR acquisition (8), the expected values for the two-way range time τ e,CR and the corresponding azimuth t e,CR can be determined by solving the range-Doppler equations that define the SAR imaging geometry [1], [48] These equations relate the orbit state vector of the satellite (position X s and velocityẊ s ) with the CR position vector X CR , and the conversion for the range time involves the speed of light in vacuum c. Note that the equations are given for the epoch of zero-Doppler. Both the Sentinel-1 IPF processor and the TerraSAR-X SAR processor TMSP generate SLC image products according to the zero-Doppler geometric convention [3], [5].
To find the CR's image coordinates (t e,CR , τ e,CR ), (10) must be evaluated in an iterative way for the time of closest approach t e,CR . The orbit is annotated in the metadata of TerraSAR-X and Sentinel-1 products as a set of discrete state vectors with time tags referenced to UTC [5], [49]. The arc of the orbit covering an image acquisition can be accurately modeled and interpolated in time by using polynomials (e.g., conventional polynomials fit by least squares [26] or the Chebyshev polynomials [50]). The brute force incremental variation of the sensor position or more refined methods, such as the Gauss-Newton iteration (described, for instance, in [51]), may be used to resolve (10) with respect to a predefined zero threshold, typically in the order of 10 −9 , which corresponds to an error of 0.7 mm when assuming a slant range of 700 km. Once t e,CR has been found, the corresponding τ e,CR is computed from (9).

5) ALE Analysis:
The ALE analysis is based on a statistical analysis of multiple data that takes over one or more CRs. The definition of the error sign is subjective. Here, we subtract the predicted timings from the observations that are corrected for the path delays. For the range, the conversion between units of time and units of length is straightforward using the speed of light in vacuum, but care must be taken to account for one-way versus two-way travel paths. The azimuth, however, is ambiguous, as it may be interpreted as an along-track error at orbit level (conversion to units of length with satellite velocity) or as location error at ground level (conversion to units of length with beam velocity sweeping the ground). For ALE definition, we use the latter.
Using the notation of Section II, the ALE in range and azimuth for a single SAR image may be given as Repeating the process for several image acquisitions results in τ CR and t CR residuals, which can be statistically quantified by the mean and the standard deviation. These key figures define the ALE, but also more detailed insights regarding different beams or different imaging modes are readily possible, provided that sufficient numbers of images have been acquired. We recommend 10-20 images to allow for a robust initial ALE analysis, whereas, for temporal stability analysis, the image series should cover at least several months. Rigorously applying the zero-Doppler model makes the ALE analysis sensitive to deviations and distortions stemming from the SAR image processing. Equations (9) and (10) are the geometric interpretation of zero-Doppler, where the azimuth time t represents the event of closest approach between the SAR antenna phase center and a target on the ground, and the corresponding distance to the target is annotated as two-way range delay time τ . This definition of the zero-Doppler coordinate system for the focused image demands orthogonal axes and timings (t, τ ) that are independent of each other. In practice, the SAR image formation algorithms implemented in operational processors of existing SAR missions may not provide exact zero-Doppler geometry for the final SAR image, as discussed in Section III. Therefore, the meticulous processing of the SAR geolocation error in combinations with a CR array provides a means to discover and ultimately correct such effects.
Finally, we want to emphasize that the ALE method is not only useful to characterize the geometric quality of an SAR system but also allows the verification of the CR ground infrastructure once the sensor has been checked and properly calibrated. In this article, we make use of both abilities: the well-calibrated TerraSAR-X mission can provide the validation of the CR coordinates of the Australian array determined with GNSS and of our methodology, whereas, for the Sentinel-1 mission, our aim is to analyze and improve the geolocation quality.

III. SAR PROCESSING EFFECTS
The TerraSAR-X and Sentinel-1 missions provide level 1 SLC products that are specified to represent radar reflectivity in the zero-Doppler coordinates [3], [5]. However, considerable offsets and distortions may still be present in the final level 1 image products if the approximations of the SAR signal model and the insufficient update rates of spatially variant focusing parameters, which are used to increase the processing efficiency, are not taken into account.
The main offset contributions arise from the widely used stop-and-go approximation, which neglects the bistatic nature of the antenna movement in azimuth as well as the Doppler shift of the signal caused by the satellite motion when focusing the range pulses [14], [48], [52]. The TerraSAR-X processor TMSP does not employ these approximations, but they are used in the Sentinel-1 IPF, which results in stop-and-go artifacts for the Sentinel-1 images. The concept underlying the stop-and-go approximation, the impact on the Sentinel-1 image products, and how to correct for in postprocessing are discussed in the Sections III-A and III-B.
Another source of potential systematic azimuth offsets is a mismatch between the modeled azimuth Doppler FM rate and the real value in TOPS or SC imaging modes [53]. The TerraSAR-X processor frequently updates the terrain height parameter used for azimuth FM-rate calculations, whereas the Sentinel-1 IPF maintains constant average values for extended areas of the observed terrain. The details of the effect regarding Sentinel-1 IW products are explained in Section III-C.
The typical orders of magnitude of these effects along with the contributions remaining in the Sentinel-1 and TerraSAR-X products are listed in Table IV.

A. Bistatic Offset Effects in Azimuth
With only a few exceptions (e.g., time-domain backprojection, [54]), most SAR image formation algorithms have to rearrange the 1-D stream of echo packets (see Fig. 5) as a 2-D raw signal matrix because this allows direct access of the synthetic apertures along the azimuth time coordinate t during SAR focusing [48]. This is achieved by introducing the two-way round trip time τ as a second coordinate.  As shown in Fig. 5, the timing of SAR pulse transmission and echo data reception is composed of a sequence of PRI. A single PRI starts with the transmission of a pulse followed by an echo window where the instrument is switched to the receive mode. The target echoes of a pulse transmitted in the i th PRI are received within the echo window of the (i + rank)th PRI. The parameter rank specifies the number of traveling pulses. Taking the leading edge of the pulse as reference, the echo of a pulse transmitted at t T x and scattered back from a target on ground is registered at receive time t Rx . The two-way pulse round trip time τ as computed from these events reads which is the equivalence of the target's range to the sensor (9). At the start of each PRI, the range time τ is reset to the initial value The projection of the received echoes onto diagonal lines in the t, τ plane reveals the direct coupling of both time coordinates (see Fig. 6). Processing the raw signal data such that this coupling is taken into account for an image focused on zero-Doppler requires careful consideration of both the quasi-bistatic imaging situation and the implications of the 2-D data arrangement.
After sampling and rearranging of the data, as shown in Fig. 7, the echo window is labeled with an azimuth index i , and each sample within the echo window is labeled with a range index k. In accordance with the timelines of Fig. 6, the corresponding receive time annotation (t Rx , τ ) of a single where t PRI#0 is the receive time label of the first recorded PRI, t SWST is the sampling window start time, and f r is the range sampling frequency. After pulse compression along each column, the echoes of targets at different ranges do not overlap any longer, and the intensity of the kth sample corresponds to the range distance For equally precise azimuth geolocation in the final SAR image, we have to keep in mind that each range measurement is performed in a quasi-bistatic SAR imaging configuration where the antenna moves along the orbit during transmit and receive. In the case of Sentinel-1, the satellite travels up to 40 m between a pulse transmit and the closing of the echo window. Thus, neither the transmit time t T x nor the receive time t Rx are the correct events in the timeline to consistently annotate the range measurements. In fact, the instant in time exactly in the middle between transmit and receive, referred to, in the following, as t bistatic , has to be used as the correct raw data timing in order to finally obtain zero-Doppler time for the focused SAR image In contrast to this exact bistatic annotation of the SAR raw data using the time computation (14), the stop-and-go approximation assumes that the satellite stops at the pulse transmit and waits until the echo signal is received and annotated. Thus, the stop-and-go approximation neglects the bistatic nature of the acquisition and labels all the range measurements within one echo window with a single tag t stop−go , which is inferred from the receive PRI holding the echo window (e.g., the event of pulse receive) The deviation in the azimuth raw data timing introduced by the stop-and-go approximation (18) leads to an equivalent deviation in the zero-Doppler annotation of the focused SAR image. Accordingly, either the raw data timing or the final azimuth annotation of focused SAR image has to be corrected by the difference of (17) and (18) when inserting (14) and (15) for t Rx (i, k) and τ (k), respectively. The range-dependent azimuth timing correction t corr is applied in full to all TerraSAR-X image products generated by the TMSP [16]. In contrast, the Sentinel-1 IPF only applies a constant shift (the so-called "bulk-correction") to the stopand-go timing (18), which equals to half of the range delay at mid swath of the TOPS center beams (IW2 or EW3) [3] Comparing (20) with (19), it is obvious that even for the range sample equivalent to τ mid , the azimuth time t IPF is not the zero-Doppler azimuth timing of this range sample. Therefore, the IPF bulk correction has to be reversed, and the correct zero-Doppler time annotation is obtained by addition of the bistatic correction (19) Both the rank and the PRI defining the τ 0 of a particular beam (IW1-3 or EW1-5) are provided in the Sentinel-1 product annotation [49]. The applicable τ mid is readily computed from the IW2 or the EW3 annotation. Fig. 8 summarizes the azimuth time correction steps required to obtain the zero-Doppler time annotation for a focused SAR image. Starting from the stop-and-go annotation of the samples (i, k) at t stop−go (red circles), first, the receive times t Rx have to be reestablished (blue arrows), followed by the bistatic delay corrections (green arrows), resulting in offset-free zero Doppler annotations (green squares), as described in (19). The bulk shift (red dashed arrow) used by the Sentinel-1 IPF yields an azimuth time annotation as indicated by the red squares, which becomes corrected for zero-Doppler times applying (21). From the sketch, it can be concluded that the IPF bulk shift would have to account for the stop-and-go annotation in order to provide the correct azimuth time annotation at τ mid , which is also the reason why the residual bistatic correction proposed in [11] could not resolve the Sentinel-1 azimuth timings.

B. Effect of Doppler Shift in Range
The principle of SAR is based on the variation of the distance between the sensor and a target [14], [48]. Consequently, there is a velocity component in the LOS direction, and the transmitted and received radar pulses are affected by a Doppler frequency shift stemming from the motion of the satellite. The target's range history can be modeled by the trigonometric relation where r 0 represents the closest approach range and v e the effective velocity [14]. The sensor-to-target velocity is found by differentiationṙ The effective velocity v e is calculated from the satellite state vector (position X s , velocityẊ s , and accelerationẌ s ) and the target position vector X CR as follows [14]: Finally, the observed Doppler shift for the radar carrier frequency f radar amounts to whereas the shift of the range bandwidth is negligible, and the effect on the location of the correlation peak is noticeable when performing the range compression of the linear frequency modulated range pulses received from the target. The range shift effect is illustrated in Fig. 9. It amounts to where K r denotes the FM-rate of the range chirp. Due to the fact that a single transmitted radar pulse is reflected by different targets with different sensor-to-target geometries, several Doppler shifts are superimposed within one received range line. Therefore, a compensation of the Doppler effect immediately after receive is not possible. However, the individual azimuth viewing angles become accessible in the course of processing when the SAR data are represented in the azimuth-Doppler domain. Therein, the Dopplerdependent range shift corrections may be readily applied as linear range frequency phase ramps in the 2-D frequency representation of the SAR signal The TMSP compensates this intrapulse motion effect in the 2-D Fourier domain, whereas, in the case of Sentinel-1, the range shifts are not accounted for by the IPF SAR processor. Here, it must be emphasized that the authors of [52] and [53] suggest the opposite sign for the spectral correction phase in (27), which is not correct. In fact, implementation used inside the TMSP follows the sign as derived from (25) and (26) for the correction given in (27). If this was not the case, we would observe distortions in the range measurements of up to 15 cm in SC for the TSX results presented in this article. For Sentinel-1 SM data that are characterized by azimuth spectra almost centered at zero-Doppler, the Doppler shift (26) changes from positive to negative during sensor approach and sensor recede along the aperture. Therefore, the effect cancels almost entirely during range compression and only a negligible range defocussing remains. In contrast to that, the TOPS raw data and its focused SAR images are characterized by the Doppler centroids strongly varying in along-track direction [55]. Accordingly, azimuth-dependent range shifts become visible in the Sentinel-1 TOPS burst images. In order to compute the corresponding range shift corrections τ corr in postprocessing, the range chirp rate K r and the Doppler centroid frequency f DC within the focused burst images at location (τ, t) have to be derived from the Sentinel-1 product annotation The details of the Sentinel-1 computations taking into account the TOPS beam steering are given in Appendix A. As depicted in Fig. 9, the correction needs to be added to the annotated range times to compensate for the effect across the range times provided by the Sentinel-1 annotation

C. Topographic Induced Shift in Azimuth Due to Doppler FM-Rate Mismatch
The azimuth location accuracy may be also affected by inaccurate modeling of the azimuth Doppler FM-rate parameter of a target's range history. For spaceborne curved orbits, the sensor-to-target distance as a function of time does not only depend on the closest approach range but also on the orientation of the LOS vector with respect to the orbital plane. This introduces a dependence of the FM-rate on the height of a target. In fact, the same effect can be used to determine the height of a target from a single SAR image [56]. For TOPS and SC mode data, a local FM-rate mismatch may introduce significant azimuth offsets.
Spectral-domain SAR processing algorithms and the formulation of their transfer function are based on the straight flight pass approximation (22) and the concept of the effective velocity, as outlined in (24). Since the value of the scalar product of the satellite's acceleration vector and the LOS vector (24) depends on the target 3-D coordinates X CR , the effective velocity v e and, furthermore, the derived Doppler azimuth FM-rate k a depend on the local terrain height as well For SAR data with the azimuth spectra centered around zero (e.g., in yaw-steered SM mode), a mismatch of k a only leads to defocussing and a phase offset, an effect known as quadratic phase error [14]. As illustrated in Fig. 10 and described in [53], the small local target bandwidth B target and a local Doppler centroid f DC of the target lead to an offset t FM-mismatch if the true geometrically conditioned Doppler FM-rate (dashed light blue) and the modeled FM-rate (red) do not coincide. The corresponding azimuth offset caused by the Doppler FM-rate mismatch is given as and the correct geometric zero-Doppler time is obtained by In the case of TerraSAR-X, the TMSP closely follows the terrain height variations when modeling the Doppler FM-rate on a burst by burst basis, which corresponds to an update rate of about 0.4 s (= SC burst duration in azimuth). This keeps the FM-rate mismatch-induced azimuth error at the millimeter level for SC. In contrast, the Sentinel IPF maintains a constant average terrain height across all TOPS subswaths for the Doppler FM-rate calculation, and the update rate is lower because of the longer azimuth duration (about 3 s) of the TOPS bursts [57]. Therefore, target height deviations of 1000 m, possibly occurring in mountainous regions, can lead to azimuth offsets on the order of 1 m [53]. The details on how to compute the correction using (31) and the Sentinel-1 product annotation can be found in Appendix B.

IV. AUSTRALIAN CORNER REFLECTOR ARRAY
In this article, we take advantage of a unique array of CR situated in Australia, which is located in the State of Queensland, approximately 180 km west of Brisbane, on the eastern side of the Australian continent (see Fig. 1). The Australian array consists of 40 permanent CR installations that are distributed across an area spanning about 130 × 130 km, and the construction was completed on November 21, 2014 [18]. The array was constructed for two main purposes: 1) to support geometrical and radiometric calibration/ validation of international satellite SAR missions; 2) to enable the comparison of interferometric measurements with colocated GNSS campaign measurements for the ongoing monitoring of ground surface displacements in a producing gas field. Each CR uses a triangular trihedral design for long-term structural rigidity. The majority of the deployed CR (34 of the 40) has an inner leg dimension of 1.5 m, while three reflectors have an inner leg dimension of 2 and 2.5 m, respectively. These six larger CRs are concentrated in the eastern portion of the array (see Fig. 1). Note that the CR numbering is given according to the geodetic survey marks covering an even wider area [18], of which the marks 7 and 33 are not equipped with a CR; thus, the highest CR label reads 42.
Deployment sites are generally within the open agricultural pasture, which provides excellent LOS visibility to the orbiting SAR sensors and low radar backscatter (clutter). SCR values derived from the deployed CR are given in Table II. The CRs are installed on concrete foundations that consist of a slab suspended above the ground level on four concrete pillars sunk to a depth of up to 3 m to promote local stability in the vertisol soils found in this part of eastern Australia (see Fig. 11). Those CRs deployed on grazing properties are enclosed by plastic perimeter fencing to protect them from animal interference. All 40 of the deployed CRs have been oriented for ascending orbital passes (assuming right-looking SAR satellites, e.g., TerraSAR-X or Sentinel-1) since the time of installation, i.e., the CR boresight is facing a westerly direction.
The ITRF coordinates of the CR were obtained via a GNSS survey that was carried out during the installation. The average uncertainties (2σ ) of the GNSS solutions are reported as 2.2 cm horizontally and 4.4 cm vertically [18]. A recent reprocessing of the GNSS data has corrected erroneous coordinates at some of the sites (see the Supplementary Material published online for [18]). The provided coordinates refer to a reference point beneath the CR structure, located on the mounting stand (see Fig. 11).
For the ALE analysis, it is required to know the position of the CR apex in three dimensions (i.e., the intersection of the three reflector plates within the open CR aperture), see Section II-A. Since that position has not currently been measured for the CR of the Australian array, we estimate the position using a model. For the known current orientation of the CR (azimuth and elevation), the Cartesian XY Z offset of the CR apex with respect to the known reference point is calculated using a computer-generated model of the CR design. These XY Z offsets are then used to obtain a coordinate for the CR apex. The impact of these phase center offsets on the ALE analysis is shown in the results of TerraSAR-X (see Section VI-A1).
The GNSS survey of the CR coordinates was processed for the ITRF epoch of January 1, 2015. In order to transform the coordinates from this epoch to the epochs of the TerraSAR-X and Sentinel-1 SAR acquisitions, the linear velocity was inferred from the ITRF solution of the GNSS stations at Sydney (IGS station SYDN) and Townsville (IGS station TOW2), both located on Australia's east coast. Considering the spatial distances to the Australian array (700 km to SYDN and 975 km to TOW2) and expressing the velocities in the local frames of the stations yields an average velocity at the array location of −0.0327 m/y, −0.0086 m/y, and 0.0496 m/y when expressed the global ITRF X, Y , and Z , respectively. This is equivalent to the typical horizontal movement of the Australian plate: 0.0553 m/y northward and 0.0233 m/y eastward [42].

A. Sentinel-1
Sentinel-1 uses a C-Band SAR payload (center frequency 5.405 GHz, equivalent to a wavelength of 5.6 cm), and its primary acquisition mode is the IW mode, which is the Sentinel-1 implementation of the TOPS SAR mode [55]. It provides three parallel swaths, termed IW1, IW2, and IW3, that are scanned in an interleaved burst pattern by switching the beam from swath to swath [see Fig. 12(a)]. The total swath width on the ground is approximately 250 km, and the average pixel resolution is 3 m × 22 m in slant range and azimuth, respectively [3]. The bursts have cross-track and along-track overlaps to ensure seamless coverage. This may result in a double coverage of targets in the overlapping areas, and the different bursts would have different Doppler frequencies. Some of the CRs of the Australian array happen to be located in such overlap areas, which makes them particularly suited to investigate the bistatic effects discussed in Section III. It is important to note that the bursts are accurately synchronized such that the pattern shown in Fig. 12(a) is closely repeated in every acquisition to allow the interferometric alignment of the scenes [3]. Therefore, any systematic behavior related to data processing can be assumed identical in repeat pass acquisitions.
After having completed the commissioning phase of the first satellite S1A in autumn 2014, the Sentinel-1 mission regularly captures data of the CR array in the standard IW mode on ascending-pass orbits once every 12 days, which is the repeat cycle of the Sentinel-1 orbit [10]. At the beginning of the mission, a few images were acquired from the pass of track 9 (the number is the relative orbit number within the 12-day repeat cycle), but since then, only the track 111 has been used. Therefore, almost all the Sentinel-1 data that we use in this article correspond to the track 111 acquisition geometry (see Table V). With the addition of the second spacecraft, namely, S1B, declared operational in October 2016, the mission now usually performs acquisition of the CR array every six days. The images are freely accessible after registering at the Sentinel-1 open access hub. 3 For our study, we make use of all the data available for the period of November 2014 to March 2018, which amount to 104 acquisitions (see the details given in Table V). 3 https://sentinel.esa.int/web/sentinel/sentinel-data-access The systematic processing scenario and the product timeliness impose that all the Sentinel-1 data-takes have to be processed within a timeframe of 24 h [10], while, in practice, most of the data are available within 4 h of acquisition. SAR processing depends on several auxiliary products, for instance, the orbit product. The detailed list of the underlying auxiliary product configuration can be found in the annotation of each Sentinel-1 product. Typically, the so-called restituted orbit that becomes available a few hours after downlink is used in the IPF to generate the SLC product [11]. For higher accuracy requirements, the user may also download POE product released 21 days after the acquisition from the Sentinel-1 quality control website. 4 Cross-comparisons of orbit solutions computed with different software packages and methods indicate an absolute accuracy of 5-cm RMS for the official POE product [58]. Moreover, the timeliness of the data flow results in the generation of the products with the latest version of the Sentinel-1 IPF operated at the time of acquisition. As reprocessing of older data is not necessarily foreseen and only performed in dedicated campaigns [10], any long-term data series will eventually contain products formed by different IPF versions, e.g., IPF v2.35 to IPF v2.84 and different orbit types, which may have an impact on the ALE analysis if there have been changes in the IPF processing baseline. Due to the tight consistency requirements regarding the interferometry, we do not expect any major impacts, but to ensure a consistent analysis, the Sentinel-1 data of the Australian CR array were reprocessed with the IPF v2.84 and the latest processing baseline (as of November 2017), including the precise orbit product. Hence, our ALE analysis is entirely based on SLC products using IPF v2.84 and using the final precise orbit ephemeris product to replace the orbit annotated in the products.

B. TerraSAR-X
The TerraSAR-X mission consists of two satellites TSX and TDX that operate in X-band with a frequency of 9.65 GHz (wavelength of 3.1 cm). Contrary to Sentinel-1, data acquisition is mainly order driven, and the user can select any of the available SAR modes. For our study, we use SM images ordered from the TerraSAR-X archive that have been captured between September 2014 and September 2015, as well as SC images acquired in the period April 2017 to February 2018 (see Table VI). Note that the overall Australian CR array was officially completed in November 2014, but individual CR installations already began as early as September 2014. As for the coverage, the SM scenes cover two smaller regions of the array, whereas the SC data contain all the CR but at the expense that the acquisitions had to be divided into two parts [see Fig. 12(b)]. Like for the Sentinel-1 IW mode, overlaps in the SC burst maintain the seamless ground coverage. The data involve two passes (relative orbits 64 and 140), of which the relative orbit 64 is limited to the eastern part of the array because this steeper geometry is already at the 20 • lower boundary of the data access range [5].
The TerraSAR-X data were generated with the latest TMSP v4.10 and the PSO product, which is the default orbit if the order does not require near-real-time delivery [5]. Long-term analysis of the PSO product with independent SLR observations confirms the very high quality of the orbits, which was assessed to be better than 2 cm [59]. Both satellites, TSX and TDX, are involved in the acquisitions used (see Table VI), but there is no particular reason to distinguish the individual sensors because we have recalibrated and verified their ALE at the 1-3-cm level [26], [60]. The long-term stable CR located at the geodetic observatory of Metsaehovi, Finland, is used to determine the range and azimuth geometric calibration constants of TSX and TDX. By using these constants in the analysis of the Australian CR array, we expect the remaining ALE offsets to be close to zero, i.e., within the 2-4-cm limit set by the accuracy of the ITRF coordinates (see Section IV) and the limit defined by the SCR of the CR in the SM and SC image products (see Section II-A).

C. Summary on Processing Setups
The important elements concerning the different steps of the ALE analysis (see Section II) and the summary on the data products as used in the processing are listed in Table VII. For both missions, the image products have been prepared with SAR processors employing consistent configurations and baselines of auxiliary products (final orbits, instrument parameters, and so on) but keep in mind that, for the TerraSAR-X mission, range and azimuth are accurately geometrically calibrated. The range and azimuth residuals presented in Section VI are the observed timings corrected for the atmospheric path delays (and the additional timing corrections required for Sentinel-1) minus the expected timings inferred from the CR reference coordinates (see Section II-B5). To convert the timing differences to units of meters, the speed of light in a vacuum is used for the range residuals (divided by two to obtain one-way range residuals), and the zero-Doppler ground track velocity is used to convert the azimuth residuals. For the latter, we assume constant values for TerraSAR-X and Sentinel-1 (see Table VII) because the variation across the swaths of the used products is only some 2-3 m/s, which means that the conversion error for the azimuth residuals found in this article does not exceed 2 mm.

A. TerraSAR-X Array Verification
The geometrical quality of the TerraSAR-X products is used to verify two aspects of the Australian CR array: first, the modeling of the CR apex required to describe the reference coordinate of a reflector (see Section IV), and, second, the geometrical fidelity of the entire array as measured with SAR. For the first task, we use the SM scenes covering two subsets of the array, while, for the second objective, all of the CRs are analyzed with the SC data.

1) Testing of CR Apex Transformation:
In total, the two stacks of TerraSAR-X SM scenes cover 11 reflectors of the array, of which five reflectors are larger CR having an inner leg dimension of 2.0 or 2.5 m. These CRs are expected to show larger residuals in the ALE analysis when not applying the transformation because the offset between the survey point and the location of the CR apex is larger than for the smaller 1.5-m reflectors. Moreover, we expect to find the residuals to be in the negative range, as the reference range predicted using the survey point (instead of the apex) is too long. This is because the difference between apex and survey point is mainly in height (see Fig. 11), which means that the apex location is closer to the sensor when observing in SAR slant-range direction. The results without the apex transformation confirm that this is indeed the case [see Fig. 13(a)]. When computing the difference of measured timings and the expected reference timings, we observe that the range residuals are all negative, and the larger CR have larger residuals [see, for instance, CR5 (dark blue signatures)]. In addition, there is a clear separation between the two passes, track 64 and track 140, because the images of the steeper geometry of track 64 are more sensitive to the difference in height. Applying the apex transformation resolves both the separation of the two geometries as well as the differences related to the CR dimension, and the overall measurements become almost perfectly centered [see Fig. 13(b)]. The remaining ALE (1σ ) across all reflectors amounts to −0.9 ± 3.5 cm in range and −2.4 ± 3.4 cm in azimuth. These results are another independent validation of our ALE methods, the high geometrical quality of the TerraSAR-X mission across the different products, and the applicability of the precisely determined TSX and TDX calibration constants. They are also a confirmation of the reprocessed CR reference coordinates and the method for modeling the CR apexes (see the discussion in Section IV) because the overall standard deviations remain well within the 2-4-cm uncertainty given for results of the GNSS survey. For a certain CR, we notice slightly larger residuals: the CR5 azimuth measurements are clustered at −7 cm, or the range measurements of CR8 are concentrated at −5 cm, which are useful indicators that the coordinates and the apex transformation of these CR could be further reexamined.
2) ALE Results With ScanSAR: The ALE results for all 40 CR using the TerraSAR-X SC data are visualized in Fig. 14(a). The appearance is very similar to the results obtained with SM data. In terms of the slant-range resolution, the SC product is almost identical to SM [5], which explains the close agreement in range ALE of 0.9 ± 4.6 cm. As for the azimuth, the result is perfectly centered at 0.4 cm, but the standard deviation (1σ ) increases to 10 cm (versus the 3.4 cm determined for SM). However, the latter is explicable when looking at the values derived from the SCR of the CRs in an SC product with a 19-m azimuth resolution (see Section II-A). Interestingly, our experimental result is even better than the approximately 20 cm estimated for a 1.5-m CR, assuming spatiotemporal random clutter. Another notable feature is the increased spread of a few measurements in azimuth [scattered signatures at the top of Fig. 14(a)]. Closer inspection of the data revealed that these measurements are almost entirely related to two products acquired in December 2017; hence, we assume some problems with these specific images even though the estimated SCRs are in line with the other data and the atmospheric corrections are comparable to other acquisition dates in the series.
By visualizing the mean values in the range and azimuth for each CR [see Fig. 14(b)], we may conclude the following: first, 75 % of the CRs are within 5 cm of the reference coordinates in both range and azimuth, which agrees with the conclusion already formed with regard to the analysis of the SM data; second, there are a few CRs with relatively large mean offsets, for instance, the reflectors 4, 9, and 42. The corresponding range and azimuth standard deviations are visualized in Fig. 15. The number of observations per CR is somewhat inhomogeneous because TerraSAR-X SC is unable to cover the whole array with one product. Despite this, we note that the CRs displaying larger offsets in Fig. 14(b) remain fairly inconspicuous in terms of their standard deviations.
Regarding the azimuth errors, we are able to give further explanations. The array offers the unique possibility to identify spatial-dependent effects related to the SAR processing, in particular, for burst related modes, such as TerraSAR-X SC or Sentinel-1 IW, because the reflectors are distributed across large parts of the product and at different positions within bursts. When arranging the range and azimuth ALE residuals according to the CR locations (via timings t and τ ) within a burst, systematic effects induced by the processing will show up as distinctive patterns. Moreover, the repeat pass configuration allows for a CR-wise temporal averaging of the residuals, which reduces the random error by the square root of available scenes when assuming normal distribution [61].
For the SC range residuals, the analysis yields only random outcomes, but the combination of azimuth residuals and azimuth time in burst reveals a clear correlation (see Fig. 16). Because of the averaging, the results of track 140 have smaller errors as this track has up to four times more samples. At first glance, the plot appears somewhat busy, but the underlying pattern becomes immediately clear if only a single reflector is traced. CR3 appears at close to t = 0 s and at t = 0.36 s, which means that the reflector is located in an along-track overlap region of two consecutive SC bursts. Moreover, there is another data point at t = 0.27 s, stemming from the adjacent burst of the next subswath; therefore, the CR3 is effectively processed three times in a single SC product acquired from track 140, and it is also processed once more in the data of track 64. The overall average of CR3 shown in Fig. 14(b) is not affected because the offsets almost cancel when combining the entire data of CR3, but, for other CRs, such as the reflectors 6, 9, or 42, this effect is the main reason for the azimuth offsets observed in the ALE results.
Such azimuth effects are often related to the bistatic correction discussed in Section III-A. For TerraSAR-X, however, the correction is fully implemented in the TMSP and consistently applied during the processing of each SAR mode (spotlight SAR, Stripmap SAR, and ScanSAR). In our previous studies with TerraSAR-X [7]- [9], we found no  TerraSAR-X ScanSAR azimuth offsets averaged for each CR according to its azimuth time t in the ScanSAR burst. The bars mark the 1σ error of the derived mean values. Labeled data points are discussed in the text.
indications for any azimuth-related distortions; thus, we suspect the effect to be specifically related to the SC mode. Also, the topographic mismatch (see Section III-C) can be ruled out, as, in the TMSP, the effective velocity parameter is defined for mid-azimuth of each burst and quadratically scaled across the range, which maintains a good approximation of the real situation because of the short duration of the SC bursts. Tests carried out for CR3 confirm that the remaining error caused by the FM-rate is only a few millimeters. As of today, we have not yet identified the cause, and further investigations will be performed in the future.

B. Sentinel-1 Analysis
The main motivation for analyzing the Sentinel-1 IW data of the CR array is the removal of the effects not considered by the IPF. First, we illustrate the behavior of the three effects discussed in Section III on behalf of the measurements of S1B. Subsequently, ALE results of S1A and S1B are compared, independently computed by DLR/TUM and UZH, and taking into account all of the IPF-related corrections.
1) Processor Effects: By applying the concept of sorting the ALE residual according to the timings within a burst (as already demonstrated with TerraSAR-X SC) and refining the results by multitemporal analysis, one can study, in detail, the postprocessing corrections required for the IPF-generated range and azimuth timings. In the following discussion, we only use the data of S1B because the data of S1A are additionally affected by an antenna event that occurred in June 2016 [62] as well as the images alternating between HH and VV polarizations during the early acquisition of the array. We address these two points, in more detail, in Section VI-B2.
The Sentinel-1 bistatic azimuth correction (21), as defined in Section III-A, is determined by the round trip time τ . Therefore, the limitation of the IPF data is clearly visible in the plot showing averaged azimuth residuals against the fast time within a burst (see Fig. 17). The IW data cover the entire array, and every data point is the outcome of averaging 30 measurements, which leads to typical error estimates of about 5 cm (1σ ) for the azimuth mean values (hence, the error bars are too small to be perceived in the figure). Due to the simplified bulk correction employed by the IPF, there is a separation between the two subswaths IW1 and IW2. This was already observed in earlier studies [12] and is caused by the fact that the receive time event and the start of the echo window in the timeline change for the subswaths due to the adaptation of the PRI. The linear trend results from the bistatic nature of the data acquisition not accounted for by the stop-andgo approximation. As it can be seen from the dashed lines, the method of removing the IPF bulk correction combined with the proper bistatic correction perfectly reproduces the situation observed with the array of CRs. Note that the correction (21) is given such that the values need to be added to the measured azimuth, but, for Fig. 17, we reversed the sign to superimpose the correction on the data. The same applies to the other corrections visualized in the subsequent figures. The small offsets between the corrections and the data are due to the lack of a refined cm-level geometrical calibration for the Sentinel-1 sensor, which is required to fine-tune and center the SAR measurements in accordance with the precise correction methods [9].
For the Doppler-induced distortion of the range measurement, we employ the same analysis method, but, here, it is the combination of slow time t in burst and the range residuals that make the effect visible. The corresponding result obtained from averaging is shown in Fig. 18. The error estimates (1σ ) for the range mean values are in the order of 1 cm. The strongest impact is found in the CRs located in the along-track overlap of two burst; thus, the three signatures at t = 0.1 s and at t = 2.9 s that belong to the CRs 2, 28, and 40, respectively, indicate the total signal in the Sentinel-1 IW data amounting to 0.8 m (see Fig. 18). The correction for the Doppler shift, as described in (28), is again closely matching the pattern found in the uncorrected IPF data. The difference in the slopes for the subswaths IW1 and IW2 is caused by the antenna steering rate and the FM-rate of the range chirp K r , which are adapted by the instrument for the different subswaths.
The third effect, i.e., the topographically induced mismatch of the azimuth FM-rate, depends on the azimuth location within a burst and the height of the CR with respect to the height assumed when computing the equivalent velocity parameter. For the CRs of the array, which is distributed across ellipsoidal elevations of 316-432 m (see Fig. 1), the error is on the order of a few centimeters (see Fig. 19), but it may become as large as 15 cm for CR2, which has the highest elevation of 432 m and is located at the very border in a burst overlapping area. Because of the overlap situation, the total mean of CR2 is hardly affected, which is identical to the impact of an FM-rate mismatch in SAR modes, such as SM or spotlight, where a mismatch does not cause a shift and leads only to defocussing [56]. The standard deviation  of the CR2 azimuth residuals, however, is reduced from 24.8 to 19.9 cm when applying the correction. Therefore, the correction is mainly noticed for individual reflectors, but it is less prominent when looking at the overall ensemble of the array reflectors.
In summary, we conclude from these results that the timing corrections, as described in the theory of Section III, are closely matched by the data of the 40 CRs, and in our next step, we can apply them to analyze the ALE of both sensors S1A and S1B.
2) ALE Results With Sentinel-1 IW: Fig. 20 illustrates the quality of the IW SLC geolocation estimation at the Australian CR array when the radar timings are processed as provided in the product. The ALE scatter in Fig. 20(a) shows the estimated geolocation accuracy of the delivered products, before correction of the perturbations and the now-understood timing errors. A mean range bias on the order of one range sample (roughly 3 m) exists, almost entirely due to signal path delay in the troposphere and the additional delay in the ionosphere. Correcting for this path delay and the signals originating from the solid Earth dynamics centers the residuals on zero range, as shown in Fig. 20(b). The overall result is well within the official 7-m specification for the geolocation of IW products [3], but we clearly perceive the systematic separation for the subswaths IW1 to IW3 in azimuth due to the bistatic effects and the skewing of the range caused by the Doppler shift.
The ALE results in Fig. 21 show the impact of the additional IPF timing refinements on range and azimuth, as described in Section III and illustrated in Section VI-B1. The results have been independently generated by UZH and DLR/TUM, using separate implementations for the SAR ALE analysis and the Sentinel-1 specific corrections (see Section V-C). S1A and S1B are presented in separate figures to analyze the difference between both sensors. For comparison, Fig. 21(a), (c), (e), and (g) shows the ALE for the level 1 SLC products as provided by the Sentinel-1 IPF without our IPF specific postprocessing, which are basically a decomposition of Fig. 20(b) for the individual Sentinel-1 sensors. Note that there was a small difference in temporal coverage for S1A because the openly accessible meteorological data used by UZH to compute the tropospheric delays (see Section II-B2) were only available from March 2015 onward. Because of this, UZH extended the period to approximately match the amount of processed S1A data.
After applying the IPF-specific corrections for range and azimuth, the scatter is significantly reduced, and the point clouds associated with the different subswaths become superimposed, with the azimuth now being much closer to zero [see Fig. 21(b), (d), (f), and (h)]. Across all the figures, we achieve very good visual correspondence between our independently computed UZH and DLR/TUM results and also observe a consistent impact of the IPF specific corrections. Looking at the numbers, the standard deviations confirm this agreement. In the case of S1A and the timing corrections applied, the values for range and azimuth (1σ ) read approximately 7 cm and 36 cm, while, for S1B, the numbers are 7 and 26 cm, respectively. Compared with the results without the additional timing corrections, this corresponds to an improvement in standard deviation by roughly a factor of 2 in azimuth and a factor of 4 in range.
There remain small mean offsets in range and azimuth, which may be interpreted as refined geometric constants for our updated processing methods, provided that they can be confirmed by ALE offsets of other validation sites as done in the tests performed with TerraSAR-X [9]. Both our analyses agree in the azimuth offsets, for which we determined about 2 cm (S1A) and −19 cm (S1B), whereas, in range, there is a systematic difference of approximately 13 cm between all the UZH and DLR/TUM results. This difference is mainly due to the methods used by our groups to correct atmospheric delay (see Section V-C), and we address this aspect, in more detail, in the subsequent discussion of the results.
An important point of discussion is the significantly larger spread in the azimuth standard deviation of S1A, which remains even after applying all of the corrections. The cause for this is most probably the outage of the SAR antenna tile #11 of S1A that occurred on June 16, 2016 [62] because the temporal history of the azimuth residuals displays a discontinuity that coincides with that event. Fig. 22 shows a temporal view of the ALE results of Fig. 21(f) and (h), for which the individual ALE results of the 40 CRs were summarized as mean range and azimuth of each acquisition. When looking at the S1A azimuth time series, the jump in June 2016 is clearly visible, matching the time when the antenna tile outage took place. The results before this event are very similar to the averaged outcomes of S1B, suggesting that the outage has affected the delay pattern of the antenna, as discussed in [62]. Interestingly, the observed change in the data of about 0.35 m is somewhat larger than the 0.1 m predicted in the report [62]. Additional contributions could be caused by a deformation of the point response, which is not considered by our PTA modeling a symmetric paraboloid surface (see Section II-B1). Moreover, there is an increased spread in the early S1A azimuth data before July 2015, which we also assume to be related to the S1A sensor, because, in the TerraSAR-X SM data of 2015, there are no signs of such disturbances.
Another interesting detail is the difference observed in the range time series regarding the HH and VV polarizations. From our verification with TerraSAR-X, we know that differences in the order of 4 cm are caused by the staggered structure of the VV and HH transmit/receive elements of the antenna and the different signal routing of the polarization channels [9]. We assume that the same holds true for Sentinel-1, and ultimately, the mission will require two sets of geometrical calibration constants to ensure consistent range measurements.
For the subsequent ALE analysis, we decided to reduce the S1A data and only use the consistent measurements after October 2016. This has the additional benefit of perfectly comparable S1A and S1B results because each sensor is now regularly capturing the array once per 12-day repeat cycle, resulting in the same amount of 30 acquisitions for each sensor for the given period. The S1A ALE for the reduced data period is shown in Fig. 21(e) and (j), and the numerical results in range and azimuth amount to 6 ± 7 cm (19 ± 6 cm for DLR/TUM) and 21 ± 29 cm, respectively. Again, we observe the already mentioned range difference between UZH and DLR/TUM, but, except for the individual offsets in range and azimuth, the results of S1A and S1B are now in much closer agreement, with S1A having only a slightly larger azimuth standard deviation.
In order to assess the now fully corrected results for any potentially unresolved timing effects, we performed, once more, an analysis of the residuals according to the slow time and fast time within the bursts (see Fig. 23). For range, there is an impressive correspondence between both sensors, and we can easily relate the individual data points of S1A and S1B for each subswath by visual inspection. Moreover, there   are no immediate indications of any systematic errors in the range data, and when looking at the scale of the spread across the CRs, we see that the sensors have now reached a level that is already sensitive to the centimeter accuracy of the individual CR reference coordinates. For azimuth, we notice a similar correspondence between identical beams of S1A and S1B, but, among the beams IW1 and IW2 of S1A, there exists a systematic offset. Again, we assume this to be related to the aforementioned antenna event because if the same analysis is performed using the consistent S1A HH data before the event (October 2015 to June 2016, see Fig. 22), we obtain a homogeneous result for both beams, which is very similar to what we observe for S1B in Fig. 23(b).

VII. DISCUSSION
The analysis performed for the two SAR missions TerraSAR-X and Sentinel-1 at the Australian CR array illustrates the usefulness in calibrating and cross-comparing two SAR systems with different payloads (X-band versus C-band) and different resolutions-provided both are carefully processed using the same standards. Table VIII summarizes the results of Section VI derived with TerraSAR-X SM/SC and Sentinel-1 IW. For the Sentinel-1 estimates, the results from TUM/DLR and UZH computed according to the setups specified in Table VII are provided side by side.
The higher resolution SM data of TerraSAR-X have an advantage due to the corresponding gain in SCR. If we recall the 2-4 cm accuracy of the CR reference coordinates as well as the additional uncertainties introduced by the atmospheric path delays (2 cm for the troposphere and sub-cm for the ionosphere, see Section II-B), the accuracy of the TerraSAR-X precise science orbit (on the order of 1.5 cm [59]), and the SCR contribution (estimated as 0.7 and 2 cm for range and azimuth, see Section II-A), these TerraSAR-X SM results can be interpreted as the geolocation limit presently achievable at the Australian CR array. By combining the assumed uncertainties as uncorrelated errors [61], we can theoretically gauge the These estimates provide a plausible description of what we found in practice for the TerraSAR-X Stripmap data. Regarding the range measurement in SC, we can expect some minor degradation due to the increase in the SCR contribution, and we may assume a slightly larger uncertainty for the CR reference coordinates (e.g., 3 cm) because, in the SC analysis, all the CRs are contributing to the overall result, and some of the CRs are additionally measured with the steeper geometry that is more sensitive to height errors. Thus, we consider the experimentally determined range result of 4.6 cm to be in line with our expectations. For azimuth, the error is dominated by the SCR and the coarser azimuth resolution, but, interestingly, we would have to assume an error of about a factor of 2 smaller for σ SCR to end up with the 10 cm observed for the array-provided that our assumption of uncorrelated clutter is true. On the observation side, we still have to account for the systematic effect found in the azimuth residuals (see Section VI-A2).
Similar considerations can be made for Sentinel-1 IW. Compared with TerraSAR-X, it is not possible to externally validate the Sentinel-1 orbits with SLR because the spacecrafts are not equipped with laser retroreflectors [58], but an orbit accuracy of 3 cm appears a reasonable assumption when looking at the cross-validation of the independently generated orbit solutions [58]. For the larger impact of the ionosphere in C-band, one may conclude an uncertainty of 2 cm using the RMS of the TEC maps (see Section II-B2). Consequently, we obtain 8.3 and 49.2 cm for σ R and σ A , respectively, when inserting the SCR-based estimates for IW and the 1.5 m CR (see Section II-A) into (33). Again, these overall estimates are driven by the SCR contribution, and as with TerraSAR-X SC, they appear too conservative compared with our experimental findings (see Table VIII). On the other hand, we can safely conclude that we are attaining the limit set by the IW product and SCR typically provided by 1.5-m CRs (or equivalent point scatterers).
Comparing the S1A and S1B ALE estimates from TUM/DLR and UZH in Table VIII, both teams are in remarkable agreement on the azimuth ALE-both the mean biases and standard deviations are in very good agreement. This is especially encouraging considering that each team has developed completely independent implementations of the various postprocessing ALE corrections. The techniques used are different regarding certain aspects (e.g., point target analysis of the image products and path delay correction, see Table VII), and yet the final results are now comparable. The agreement reflects the common theoretical framework for these corrections that now exists in the published literature.
However, one notable exception to the agreement between results from the two teams is the mean range ALE for Sentinel-1. While UZH estimates an offset of 6 cm for the most recent S1A data and −9 cm for S1B, the corresponding values from TUM/DLR value are 19 and 3 cm. In other words, a roughly 13-cm difference in the mean offsets exists between the two teams. As the range ALE is mainly affected by the atmospheric path delay, it is likely that the different methods used by the teams are the cause of the discrepancy. A comparison of the averaged total atmospheric path delay (sum of troposphere and ionosphere) underlying the results of S1A (data after October 2016, see Table VIII) yielded values of 3.257 m (UZH) and 3,089 m (TUM/DLR), resulting in a 17-cm difference induced by the path delay correction. This is a clear indicator that the atmospheric path delay accounts for most of the observed differences between the two groups, and future collaboration between the teams may help determine which aspects of the atmospheric models are most significant. For instance, there are approximations in the mapping applied by UZH when computing the ionospheric path delay correction (see Section II-B2), which will be refined in a future update.
The difference also points out a problem when performing a geometrical calibration of SAR sensors. The geometrical calibration includes not only the unknown contribution from internal electronic delays but also any bias introduced to the SAR system by the applied correction methods, a situation that was also encountered with the initial geometrical calibration of the TerraSAR-X mission [9]. Therefore, it is important to Fig. 24. Average range offset for the individual CRs computed from range residuals of the experimentally calibrated Sentinel-1 IW data (see Table VIII) and the TerraSAR-X SC data (range results of Fig. 14).
compare the ALE results of a sensor across several test sites if the aim is a centimeter-level geometrical calibration that ensures consistent geolocation accuracy at a global scale. This is the next important step for Sentinel-1 in order to confirm the mean offset determined for S1A and S1B at the Australian CR array. The comparison of the analysis with additional sites will verify the applicability of the found ALE offsets as geometrical calibration constants. Especially, S1A needs to be carefully examined because of the consequences of the antenna event on the azimuth ALE.
As a first outlook on the upcoming possibilities with such a type of data, we may use the overall range offsets of S1A (reduced data set) and S1B, as given in Table VIII, to correct the measurements, which should make the results comparable with TerraSAR-X, because, as an ensemble, the CRs were found to be centered (see the SC results in Table VIII). Fig. 24 shows the comparison of the TerraSAR-X SC (see Fig. 14) and the experimentally calibrated Sentinel-1 sensors using the mean DLR/TUM results across the different reflectors. Due to the possibility of averaging 30 measurements, the Sentinel-1 range results can be considered as reliable as the TerraSAR-X findings, and the side-by-side comparison confirms the very good agreement for almost all of the reflectors, especially for the CRs with the largest offsets. We notice that the same CRs already pointed out in the TerraSAR-X analysis, namely, the reflectors 4, 5, 9, 14, and 22, are also identified with S1A and S1B, which supports our earlier assessment with TerraSAR-X that these reflectors are the primary candidates for which the reference coordinates (and possibly the apex transformation) should be reexamined.
With postprocessing of the IPF timings, the Sentinel-1 IW products have reached a very high level of geolocation accuracy. Given that the geometrical calibration constants can be confirmed stable at a global scale, the products are comparable with the TerraSAR-X SC data offering a similar resolution, and in the long run, the Sentinel-1 mission will generate significant benefit from the very large amount of data that are continuously acquired for the entire Earth. It also paves the way for usage of different sensors for common applications, e.g., monitoring of commonly accessible point targets in range and azimuth. The example of CR 15 in Fig. 25 illustrates this by showing the joint time series of the range residuals of Sentinel-1A/B and TerraSAR-X. Since all the data were  Table VIII have been removed  to center the Sentinel-1 results. consistently processed and the averaged mean values of S1A and S1B have been used as preliminary calibration constants to generate these Sentinel-1 results, the ALE residuals of the different sensors become directly aligned. Moreover, we see the advantage of having continuous coverage guaranteed by the Sentinel-1 acquisition plan. For TerraSAR-X, we had to order the images and not all the acquisitions could be immediately fulfilled due to a shorter duty cycle and narrower swath and conflict with other orders in the same region, while S1A and S1B allow a regular repeat pass sampling of up to six days, as it is the case for the Australia CR array. This regular coverage makes Sentinel-1 well suited for long-term monitoring applications, even with the IW product and its seemingly coarse resolution, because it is readily possible to apply monthly or bimonthly averaging to generate reliable ranging time series. The IW azimuth will probably play a less important role in such type of applications because of the SCR versus resolution limitation, but, for the Sentinel-1 SM product with a typical slant range by azimuth resolution of 3 × 4.5 m, we can expect the ALE to be on the order of 4-5 cm in both directions. Therefore, the verification of the IPF specific corrections and homogenization of the geometrical calibration across Sentinel-1 IW and Sentinel-1 SM products should be envisaged.
VIII. CONCLUSION Highly accurate methods for computing and analyzing the geolocation of SAR satellites were described and applied to TerraSAR-X and Sentinel-1 products acquired for the Australian CR array. Presented as a comprehensive step-bystep computational approach, precise point target analysis, atmospheric path delay correction, and the modeling of geodynamic effects in the CR coordinates were coupled with reliable ground installations, accurate GNSS survey, and orbit products of the highest available quality.
For TerraSAR-X SM products, the combination of coordinates and apex transformation was found to be consistent with the 2-4-cm accuracy reported for the GNSS survey results compared with the TerraSAR-X geolocation of the 11 CRs covered by SM images. An overall evaluation of the array with TerraSAR-X SC confirmed these findings; the corresponding SC ALE (1σ ) across all 40 CRs was determined with −0.009 ± 0.046 and −0.004 ± 0.100 m in range and azimuth, respectively. In addition, we used the array to assess the geometric consistency of the TerraSAR-X SC bursts. A systematic distortion of approximately 10 cm in azimuth was identified, which was previously unknown and that will be subject to further investigations.
For Sentinel-1, IPF-specific corrections to remove the bistatic effect in azimuth and the Doppler-based distortion in the range have been derived from general radar theory. Both have been verified by the geolocation processing carried out at the Australian CR array. These corrections strongly reduce the swath-dependent azimuth offsets found in earlier studies of the IW product and a skewing of about 1 m in the range residuals across the different CRs. Moreover, the azimuth shifts caused by the mismatch of the FM-rate in azimuth focusing have been taken into account and were found to be as large as 15 cm for specific CRs of the array.
With all of the corrections applied, we determined the ALE of the Sentinel-1 IW product with 0.033 ± 0.060 m in range and −0.181 ± 0.264 m in azimuth when using the data of S1B. The remaining offsets can be compensated by the geometrical recalibration, once these findings have been verified with other test sites. These Sentinel-1 ALE results were generated by UZH and TUM/DLR using independent implementations, as summarized in Table VII. The very good agreement of the results confirms our common understanding of SAR geolocation processing theory, especially in terms of the postprocessing required for the IPF-generated timings. The only significant difference discovered in our results is a shift of approximately 13 cm between the range results, which could be traced to the different strategies used to compute the atmospheric path delays. Possible refinements have been identified in the ionospheric path delay computation at UZH, and further comparisons between the groups may help to isolate the other contributions to the difference, once the update has been completed.
The remaining range and azimuth standard deviations determined for Sentinel-1 are explained by the different uncertainties stemming from the CR coordinates, the orbit, the atmospheric corrections, and the SCR of the reflectors in the IW products. The latter is the main reason for the 26-cm azimuth standard deviation, whereas, for range, we assume similar contributions by all of the aforementioned aspects. Processing the S1A data for the same period as S1B yields almost identical results, but, in the overall S1A data dating back to October 2014, we discovered systematic effects in azimuth. These are most likely related to the outage of tile #11 of the S1A SAR antenna and will require the special emphasis when performing a refined geometrical calibration of the S1A satellite to ensure consistent accuracy in azimuth throughout the mission. In addition, we found different offsets in the S1A range residuals that are linked to the VV and HH polarization channels. Similar differences were observed with TerraSAR-X in earlier studies and may be compensated by dedicated geometric range calibration constants for each polarization.
When averaging multitemporal observations, both Sentinel-1 and TerraSAR-X are now at a level where we expect to see the quality of the GNSS-based reference coordinates of the individual CRs. Five reflectors with range offset larger than 5 cm were identified by SAR. While, for the overall analysis of SAR products, this was of minor consequence, given that the total ensemble of CRs was still close to zero means, it can be noticed in the burst-specific analysis of the range measurements. Therefore, the reference coordinates of these reflectors could be resurveyed to see if they still match the initially determined coordinates or if there are differences that may account for the SAR results. New terrestrial surveys that directly observe the apex have already been conducted at all 40 CRs in the Australian array and updated reference coordinates will be available from Geoscience Australia in due course.
These latest results with Sentinel-1 IW products confirm that the mission's geolocation capability has become very close to the comparable SC products of TerraSAR-X when applying the refinements to the IPF timings, as shown in the comparison of all ALE results given in Table VIII. In principle, these refinements could now be implemented in the operational IPF to simplify the geolocation of the standard product, but this would cause inconsistencies with already processed data in the archives. As the reprocessing of the entire Sentinel-1 archive and a renewed distribution of updated products to the users is therefore not considered feasible, it was decided to keep these IPF refinements as a postprocessing step to be applied by the user and to provide supportive documentation to make the implementation as straightforward as possible. Nevertheless, the IPF is currently being extended to support the refined timing computation for testing and validation purposes across the different Sentinel-1 products. Finally, we would like to emphasize that the level of the Sentinel-1 geolocation results presented in this article is already far beyond the initially specified 7-m geolocation requirement of the IW product [3].
Our results underline the possibility of accurate range and azimuth observations if modern spaceborne SAR sensors are processed with meticulous care in all the steps involved in geolocation. Maintaining the accuracy in each step is the key to the proposed methodology. Furthermore, ensuring such quality across differences missions is the next step toward global SAR imaging geodesy which requires dedicated reference sites. The unique wide-area Australian CR array is now confirmed to be such an accurate ground infrastructure for testing SAR satellites. These permanent CRs not only enable the monitoring and calibration of different SAR sensors but ultimately a comprehensive usage of SAR range and azimuth observations across different wavelengths, as demonstrated for Sentinel-1 and TerraSAR-X. We actively encourage other SAR satellite operators to join these activities.

APPENDIX A SENTINEL-1 RANGE CORRECTIONS FOR DOPPLER SHIFT IN TOPS IW DATA
The correction of the Doppler-related shifts for any location τ, t in a Sentinel-1 IW burst requires the modulation rate K r of the range chirp and the Doppler centroid frequency f DC , at the location (τ, t) of the CR (see details in Section III-B) While the chirp rate is directly given in the annotation of a Sentinel-1 SLC product (xml field name txPulseRampRate [49]), the Doppler centroid frequency comprises the geometric Doppler and the contributions from the TOPS antenna steering and has to be reconstructed by combining several elements of the product annotation. Parts of these computations were already described in detail by the ESA technical on the TOPS single-look complex data deramping [63]; hence, we can reuse the applicable equations for our description. Note that some of the symbols were adapted to match the notation of our article, but each of the equations is exactly referred to corresponding equation in [63] to avoid confusion. The geometric Doppler centroid frequency is annotated as a second-order polynomial for each burst contained in the product. The polynomials refer to the mid-azimuth time t mid of the bursts and are evaluated with respect to a reference slant-range time τ 0 (not to be confused with the τ 0 used in Section III). Both t mid and τ 0 are stored in XML annotation along with the polynomial coefficients. The Doppler centroid frequency for a CR observed at τ = τ m,CR (as given by the PTA) is, thus, computed by [63, equation 13] The contribution stemming from the antenna steering is computed from the Doppler centroid rate in the focused TOPS data k t , which combines the azimuth Doppler FM-rate k a and the Doppler centroid rate k s introduced by the antenna steering [63, equation 2] k t (τ ) = k a (τ ) · k s k a (τ ) − k s .
Like the Doppler centroid frequency, the azimuth Doppler FM-rate k a is given by polynomials, which follows the same computation scheme [ where v s is space craft velocity computed from the annotated orbital state vectors, the f c is the radar frequency of Sentinel-1, the k ψ is the antenna steering rate (required in radians per second but annotated in degrees per second), and c is the speed of light in a vacuum. Because the annotated polynomials refer to t mid of a burst, we can use the linear Doppler centroid rate of (36) to calculate the Doppler centroid required for the correction with t = t m,CR as given by the PTA.

APPENDIX B SENTINEL-1 AZIMUTH CORRECTIONS FOR TOPOGRAPHY INDUCED FM-RATE MISMATCH
The correction of the topographic mismatch, as documented in Section III-C requires the azimuth Doppler FM-rate k a used by the IPF to perform the azimuth focusing, which was already described for the range correction in Appendix A. The same holds true for the Doppler centroid frequency f DC (τ, t). Both values have to be calculated for the timings of the CR, which are known from the PTA, i.e., τ = τ m,CR and t = t m,CR . Finally, the true azimuth Doppler FM-rate is given by the sensor-totarget geometry, which can be computed from the satellite state vector and the position of the CR k geo (t) = − 2 λ · |X s (t) − X CR | . . . ((X s (t) − X CR ) ·Ẍ s (t) +Ẋ s (t) ·Ẋ s (t)). (40) λ is the wavelength of the Sentinel-1 carrier (the exact value can be inferred from the radar frequency of the product annotation); X CR is the position vector of the CR known from the survey; and X s ,Ẋ s , andẌ s denote the satellite position vector, velocity vector, and acceleration vector, respectively, which need to be computed for the azimuth time of the CR. The values may be calculated with interpolation methods from the position and velocity vectors annotated to the Sentinel-1 IW product, as it was already proposed during the computation of the true zero-Doppler time (see Section II-B4).