InSAR‐Derived Horizontal Velocities in a Global Reference Frame

Interferometric Synthetic Aperture Radar is used to measure deformation rates over continents to constrain tectonic processes. The resulting velocity measurements are only relative, however, due to unknown integer ambiguities introduced during propagation of the signal through the atmosphere. These ambiguities mostly cancel when using spectral diversity to estimate along‐track motion, allowing measurements to be made with respect to a global reference frame. Here, we calculate along‐track velocities for a partial global data set of Sentinel‐1 acquisitions and find good agreement with ITRF2014 plate motion model and measurements from GPS stations. We include corrections for solid‐earth tides and gradients of ionospheric total electron content. Combining data from ascending and descending orbits we are able to estimate north and east velocities over 250 × 250 km areas and their accuracy of 4 and 23 mm/year, respectively. These “absolute” measurements can be particularly useful for global velocity and strain rate estimation where GNSS measurements are sparse.

The Terrain Observation with Progressive Scan (TOPS) approach that is used for the Interferometric Wideswath mode does, however, provide much greater spectral diversity in burst overlap regions, which are observed from two look angles differing by around 1°. Differencing backward-looking and forward-looking interferograms over the same burst overlap region allows estimation of along-track displacement with a precision of around 1 mm when averaged over whole overlap regions (Grandin et al., 2016).
Time series approaches have also been developed to estimate along-track velocities in the burst overlap regions (Hooper et al., 2020;Li et al., 2021). These studies treat the along-track velocities as relative measurements, but here we explore how accurately we can measure them in a global reference frame, as the along-track measurements are with respect to reference satellite positions, localized within the Earth-centered Earth-fixed framework of the ITRF2014.
We estimate the velocities in relatively large blocks containing, on average, 36 burst overlaps, in order to observe large-scale tectonic motion and characterize some of the error sources, such as the solid-earth tides and ionosphere, based on external models. In addition, we identify and correct systematic discrepancies in the along-track offset values related to a change in orbital ephemerides in mid-2020, and report differences between measurements from Sentinel-1A and Sentinel-1B satellites.

Estimating Along-Track Displacements
We extract the Level-1 single-look complex (SLC) data of the standard Sentinel-1 Interferometric Wide Swath burst units, and merge those bursts into larger spatial frames that we have defined. One typical frame consists of 13 bursts in each of the three subswaths and covers an area of approx. 250 × 250 km. We perform a standard procedure to coregister the SLC data for each acquisition toward a reference epoch . During this process, values of the mean subpixel shift in the azimuth direction, Δa, are estimated.
In detail, we use the Copernicus Sentinel-1 Precise Orbit Determination ephemerides to resample each SLC into the geometry of the reference epoch, considering topography with respect to the WGS-84 reference ellipsoid . We then run several iterations of intensity cross correlation to estimate a refined sub-pixel shift in both slant range (across-track) and azimuth (along-track), Δa ICC , directions, with a precision below 0.01 pixels. Next, we resample the SLC to an intermediate product and iteratively estimate a further refined azimuth coregistration offset Δa SD using spectral diversity of burst overlaps (De Zan & Monti Guarnieri, 2006).
Burst overlap areas are imaged with multiple lines of sight, and taking the difference between interferograms formed for each line of sight cancels the across-track displacement, leading to an estimate of the azimuth offset with expected precision below 0.0005 pixels. To ensure coherence, this offset is estimated between pairs of acquisitions close in time that have been resampled to the reference epoch, thus already shifted by their azimuth offset. The spectral diversity phase Δϕ SD (Yagüe-Martínez et al., 2016) is related to azimuth pixel offset Δa SD by where PRF = 486.486 Hz is the pulse repetition frequency (azimuth sampling rate) of the Sentinel-1 system, and Δf DC is the Doppler centroid frequency difference in the burst overlap area.
We sum the subpixel offsets calculated from cross correlation and spectral diversity to give Δa values with respect to the reference frame of the orbit ephemerides, which since 16 February 2017, i.e., since version 1.3 of the precise orbit generating system (Peter et al., 2021a(Peter et al., , 2021b, is the ITRF2014 realisation of the International Terrestrial Reference System. We extract relevant metadata per each frame, such as the average ground footprint heading angle, α, and the pixel resolution in the azimuth direction, r azi . We use r azi to convert the frame-wise azimuth pixel offsets, Δa, to the along-track displacement u az , In total, we selected 107,476 Δa values from 1,063 LiCSAR frames, that have a minimum count of 30 Δa values per frame. As the LiCSAR processing currently concentrates on the Alpine-Himalayan orogenic belt (AHB), we 3 of 11 provide outputs of our analyses for both the full global data set and the AHB subset, defined by a bounding box of 25°W−110°W and 25°N−45°N.

Contributions to the Along-Track Displacement
Several factors contribute to displacement measured in the along track direction: where u az,tecton is the motion due to tectonic displacement between acquisitions, u az,tide is displacement due to the difference in the solid-earth tides between acquisitions, u az,iono is due to the change in the along-track gradient of total electron content (TEC) between acquisitions, and ɛ is any residual due to orbit inaccuracies and noise. As our aim is to isolate u az,tecton , we estimate and subtract values for u az,tide and u az,iono .

Solid-Earth Tide Contribution
Solid-earth tides are calculated as the response to lunisolar gravitational attraction, computed from low-precision geocentric coordinates for the Moon and the Sun (Petit & Luzum, 2013), as implemented in the GMT software (Wessel et al., 2019). We expect the tide model to be sufficiently accurate for this analysis, although updated tide models exist, which also consider ocean tide loading or polar motion (Ducarme & Schüller, 2018;Martens et al., 2019).
We use the solid-earth tide model and calculate displacement in the along-track direction, u az,tide , for the centre coordinates of each frame from the tidal components d E , d N in the east and north, respectively, using Equation 4.

Ionospheric Contribution
For the ionospheric contribution to the along-track displacement, u az,iono , we apply the IRI2016 model (Bilitza et al., 2017) that we consider sufficiently accurate for investigating seasonal fluctuations of u az caused by ionospheric influence.
We use the inverse of the approach for estimating the ionospheric influence on spectral diversity values (Gomba et al., 2016), and estimate u az,iono based on modeled TEC between the satellite and the ground at the center of the frame, for two hypothetical forward-looking (A) and backward-looking (B) bursts, that is in a geometry similar to the solution of Liang et al. (2019).
Based on the look angle geometries, we calculate geographic coordinates of two ionosphere piercing points located in the line of sight toward the central point of the overlap between both hypothetical central bursts A, B at the altitude of the peak of ionospheric F2 layer, H iono , extracted for each acquisition epoch by IRI2016. We then use IRI2016 to extract vertical TEC columns (TEC v ) between the ground and the approximate satellite altitude, at those coordinates and time.
Afterward, we estimate the average incidence angle below the ionosphere piercing points Θ IPP (Ya'acob et al., 2010) and use it to convert the modeled TEC v columns to the slant (line-of-sight) direction, TEC s , as: where R is the radius of the Earth.
The SAR carrier wave of frequency f 0 is modulated by the Doppler effect from the satellite motion and the electronic beam steering during the TOPS acquisition mode, in the frequency range f A , f B = f 0 ± 0.5Δf DC at the burst edges (Grandin et al., 2016). The phase advance ϕ iono of an electromagnetic wave, oscillating in frequency f, induced during propagation through ionosphere is (Gomba et al., 2016): 4 of 11 where = 40.308193 3 2 is a constant for which we assume a commonly applied value (Hoque & Jakowski, 2012) and c is the speed of light in vacuum. Yagüe-Martínez et al. (2016) and Equations 1, 2, and 6, we can derive the ionosphere-induced along-track displacement u az,iono toward the frame reference epoch as

From Equation 5 of
where ΔTEC s (A, B) is the difference between the TEC s estimates for the reference epoch and each acquisition in the line of sight toward the ground center of the overlap between bursts A and B, respectively. The difference of Doppler centroid frequencies for those lines of sight is Δf DC = f A − f B . We use the average value of Δf DC , as it varies between swaths (Grandin et al., 2016). Figure 1 demonstrates the influence of the corrections applied to the original u az values. The time series are generated from u az values adjusted by the median value for each frame and filtered using a moving 10-sample median filter.

Impact of the Corrections
Color represents the absolute latitude of each frame center, and shows how the ionospheric gradient is strongest at low latitudes and during ascending (evening) acquisitions. The seasonal nature of the tide and ionospheric corrections is also apparent, as well as a general decrease in ionospheric gradient between 2015 and 2019 as the solar cycle waned.

Effect of Updated Orbit Ephemerides
The orbit ephemerides products of Sentinel-1A and Sentinel-1B changed on 29th and 30 July 2020, respectively, to version 1.7, incorporating a correction of the on-board GPS antenna reference point position, as identified by Peter et al. (2020). This correction implied a 3-D position change of approximately 60 mm and led to a decrease in the root mean square error (RMSE) of the 3-D position of the satellites to below 10 mm (Peter et al., 2021a(Peter et al., , 2021b. In the along-track direction of the satellite, the shift in the antenna reference point, d ARP , is equal to 39 mm (Tables 3-1  We analyze our u az values corrected for both solid-earth tides and ionosphere as described in Text S1 in Supporting Information S1, confirming the change in the orbit ephemeride products cause an offset of d ARP = −39 mm, and we apply this correction. As described in Texts S1.1 and S1.2 in Supporting Information S1, we additionally observe systematic offsets between both Sentinel-1A and Sentinel-1B satellites that are not a single value. As the apparent offsets between the 1A and 1B satellites have no consistent pattern and, as tested, using only one satellite would increase the overall standard error, we ignore these offsets. Yet, on the example of Maduo earthquake (Chen et al., 2021) in Text S1.1 in Supporting Information S1, we demonstrate that such offset must be considered to prevent misinterpretation as large-scale deformation.
Finally, we evaluate that the improved precision of satellite positioning by the updated orbits decreases RMSE of u az -based estimates by 12%, or by median of 4.3 mm, in Text S1.3 in Supporting Information S1.

Estimation of Plate Motion
For the final velocity estimation, we invert the corrected data set of u az values by linear regression with the Huber loss function (Huber, 2009), applied with parameters alpha = 1.0 and epsilon = 1.35, per frame. We calculate the mean square error, 2 , assuming two degrees of freedom. We use the mean square error as an estimate of the variance for individual measurements and then apply Equation 8 to propagate the errors, where Q m is the variance-covariance matrix of the model parameters, −1 = 1∕ 2 is the inverted N × N variance-covariance matrix of the measurements (I N being the N × N identity matrix), and = [ , .., ; 1, .., 1] is an N × 2 matrix where t is the acquisition time in years. The first element of Q m is the estimate of the variance of velocity 2 . Estimates of velocities in ascending and descending tracks are plotted in Figure 2.

of 11
Following modified procedures for line-of-sight decomposition (Wright et al., 2004), we decompose the final u az -based estimates of horizontal displacement velocity, v, in the satellite azimuth direction from descending and ascending tracks, to eastwards and northwards direction components, v E and v N , respectively. We establish a global grid with a pixel size of approx. 250 × 250 km, map overlapping frames having their centroid inside the common grid cell, and calculate = [ , ] for pixels with = (i > 0, j > 0 being numbers of overlapping frames from descending and ascending tracks, respectively), by a least squares inversion of where matrix A contains look vector transformation coefficients per each of n = i + j grouped frames as = [sin 1..sin , cos 1..cos ] . We evaluate the accuracy of the estimates using propagation of errors as before by Equation 8.
As the azimuth direction contains only a small component of eastwards motion, the eastwards velocity estimates are very sensitive to outliers. We therefore filter our data set of velocity estimates, removing 38 grid cells that have unreasonably large eastwards velocities, |v E | > 200 mm/year, or RMSE, σ v,E > 50 mm/year.

Comparison to Reference Data
We first assess the accuracy of our global estimates by comparison with the ITRF2014 plate motion model (Altamimi et al., 2017) velocities, v PMM , averaged per set of coordinates within each of the grid cells. We carry out a statistical evaluation of the quality of our estimates, noting that the extracted plate motion model velocities do not reflect actual ground motion in deforming plate boundary zones, leading to expected real differences between the two data sets.
We also perform the comparison over the Alpine-Himalayan Belt with velocities estimated from a global network of GPS measurements (Kreemer et al., 2014), v GPS . The GPS measurements are provided in ITRF2008 no-net-rotation frame. We compare plate motion model velocities of ITRF2008, and observe that they differ by less than 0.5 mm/year on average from the ITRF2014 model over the Alpine-Himalayan Belt. The averaged GPS velocities and the plate motion model velocities differ on average by 8 mm/year northwards and 4 mm/year eastwards, in this plate boundary zone.
We first convert the reference velocities to the azimuth direction of the satellite measurements as in Equation 4 and compare to our along-track velocity estimates, grouped by direction of orbital tracks (ascending or descending). We then compare east-west and north-south velocities to our derived estimates.

Evaluation of the Final Accuracy
Median values of v and corresponding RMSE, σ v , for both our full data set (253 grid cells) and the Alpine-Himalayan Belt subset (183 grid cells from 279 ascending and 255 descending frames) are given in Table 1, together with corresponding median values from reference datasets, and .
The expected accuracy of our velocity estimates is given by their RMSE estimate. The final median 2-sigma error, 2̃ , over the Alpine-Himalayan Belt is 3.6 mm/year northwards and 20.0 mm/year eastwards.
The median difference of velocities toward GPS estimates are reported in Table 1b. Our final estimates deviate from GPS by 13.4 mm/year eastwards and −8 mm/year northwards, on average. This is higher than we would expect from the supplied orbit errors and the expected precision of our measurements, but we note that these accuracy estimates also include errors in the GPS estimates, differences between ITRF2008 and 2014, errors introduced by differences in position between GPS and InSAR measurements, and residual ionospheric errors. Note. The subscript 0 refers to the original data set before the corrections, tc refers to tide-corrected and tic to tide-and ionosphere-corrected.

Visualization of the Final Results
We visualize the horizontal velocities from our final u az -based estimates after both solid-earth tide and ionosphere corrections and compare them to the ITRF2014 plate motion model in Figure 3a. We then convert our estimates to the Eurasia-fixed reference frame by subtracting the ITRF2014 no-net-rotation plate motion velocities and adding velocities in the Eurasia-fixed model, and compare to the GPS velocities provided in this frame by Kreemer et al. (2014), in Figure 3b.
The figure highlights expected spatial patterns of deviation from the plate motion model in deforming regions, for example, an increased westward motion of western Turkey ), but it also shows higher eastwards velocity of the eastern part of the Alpine-Himalayan Belt than measured by GPS, and other discrepancies.

Discussion
We extract along-track u az estimates of large-scale deformation from averaged along-track pixel offsets Δa of Sentinel-1 data. The coregistration procedure refined by spectral diversity aims to estimate Δa with a precision of 0.0005 pixels, corresponding to u az = 7 mm, although theoretically the method can reach precision below 1 mm (Grandin et al., 2016). The reliability of the estimate depends on interferometric coherence, which decreases with time. We therefore use previously resampled data close in time to estimate the spectral diversity, but estimation errors may thus propagate in time. We did not implement an approach to reduce the error propagation as in (Fattahi, Agram, & Simons, 2017).  (Kreemer et al., 2014) in Eurasia-fixed ITRF2014 reference frame. Error ellipses represent 2-σ (root mean square error) in respective directions, neglecting their possible correlation. Figure includes global active faults (Styron & Pagani, 2020) and tectonic plate boundaries (Hasterok et al., 2022).

of 11
The list of factors influencing u az is not comprehensive (Gisinger et al., 2022;Yunjun et al., 2022). We limited our experiments to general models of solid-earth tides and ionosphere, the latter of which primarily affects ascending (dusk) frames, especially at lower latitudes. The accuracy of our velocity estimates significantly improves after incorporating corrections for solid-earth tides and ionospheric propagation, although we note large contributions to u az values at the end of the Solar Cycle 24 that are not fully corrected (Figure 1).
An influence of the absolute magnitude of TEC is not expected to contribute to the along-track displacement estimate (Fattahi, Simons, & Agram, 2017;Gomba et al., 2016), yet it affects the intensity cross-correlation estimate that is then mitigated by spectral diversity, see Text S2 in Supporting Information S1. The ionospheric correction would benefit from separate estimation of its induced phase advance gradients per satellite sub-swath, which is planned future work, and also from TEC variation estimated from the interferometric phase (Gomba et al., 2016).
We also considered u az due to troposphere delay gradients at the global mean height of water vapor (Bock et al., 2013) based on GACOS data (Yu et al., 2018), and decided to neglect this term, as it should induce less than 1 mm shift (see Text S1.1 in Supporting Information S1).
Our investigation on the updated orbits in version 1.7 (see Section 4 and Text S1 in Supporting Information S1) confirms that the updated orbit ephemeride products shift u az values after the end of July 2020 by −39 mm. We estimate the offset as −35.5 ± 4.3 mm at 2-sigma error with a median value of −39 mm. Therefore, we apply this reported along-track offset of the GPS antenna reference point directly, in order to combine the u az values with mixed orbits, without having to reprocess all data formed using the old orbits.
In Text S1 in Supporting Information S1, we also observe an increased model accuracy for u az values using the updated, more accurate, orbits. The RMSE of described model decreased by 12% in median when using updated orbits. Thus, reprocessing data with the updated orbits would improve accuracy of our estimates.
In addition, we notice apparent offsets between Sentinel-1A and Sentinel-1B satellites that appear non-constant and unchanged after the update of orbits, see Figure S1 in Supporting Information S1. We observe no spatial pattern in these offsets. As our current data set contains relatively small amount of data values from Sentinel-1B, we could not reliably estimate this term, which therefore decreases accuracy of our velocity estimates.
Finally, we observe increased variance of data points before mid-2016, which often increases overall RMSE.
Although an outlier detection technique was used to estimate velocities, an extended data set would allow dropping data from this time period while further improving the estimation accuracy.

Conclusions
We demonstrate the possibility of recovering accurate estimates of large-scale horizontal motion using along-track measurements of Sentinel-1 averaged over approximately 250 × 250 km cells. Such measurements can also be made at much higher resolution, albeit with lower accuracy.
The u az estimates can be considered along-track measurements that are given in a global Earth-centered Earth-fixed no-net-rotation framework due to the GPS-based positioning of Sentinel-1 satellites. Our final u az -based velocity estimates have median 2-sigma errors of 4 mm/year northwards and 23 mm/year eastwards (see Table 1). They fit reasonably well with the ITRF2014 plate motion model, although there are discrepancies, some of which are expected at plate boundary deformation zones, and some of which are not, such as extra eastwards motion over parts of the Alpine-Himalayan Belt.
Using averaged velocity estimates from GPS measurements over the Alpine-Himalayan Belt as reference, we identify a median bias of our estimates corresponding to 13.4 mm/year eastwards and −8.0 mm/year northwards, which is larger than we would expect, and further investigation into the cause is needed.
Accuracy can be further improved by incorporating better models of error sources, for example, an improved model of ionosphere (Galkin et al., 2022), or other sources of error (Gisinger et al., 2021) such as error due to bi-static satellite configuration, ocean tide loading or polar tide (Martens et al., 2019). Finally, accuracy can be further improved by reprocessing all data with the latest orbits and incorporating large-scale across-track (range) offset measurements, which have lower precision, but are more sensitive to east-west motion.

Data Availability Statement
We provide extracted data and estimates used within this article as an indexed Supporting information Data Set S1 and in a FAIR-compliant archive at Lazecký et al. (2022) [CC-4.0]. The software code for this research is available as Supporting information Data Set S2 as well as in Lazecký et al. (2022) [GNU GPL 3.0], and will be further updated in https://github.com/comet-licsar/daz [GNU GPL 3.0]. This work contains modified Copernicus Sentinel-1 data [2014][2015][2016][2017][2018][2019][2020][2021] available from https://scihub.copernicus.eu and Copernicus Sentinel-1 Precise Orbit Determination products available from Copernicus Sentinels POD Data Hub https://scihub.copernicus. eu/gnss. The data set of GPS measurements used in this research is included in the supplementary information files of Kreemer et al. (2014). The plate motion model values for ITRF2014 and ITRF2008 were extracted from the UNAVCO plate motion calculator at https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html. We also used available open-source models to perform corrections on solid-earth tides (Wessel et al., 2019) and ionosphere (Bilitza et al., 2017). Figure 1 and Figure S1 in Supporting Information S1 were made with Matplotlib available under the Matplotlib license at https://matplotlib.org. Figure 2 was created using QGIS available under GNU GPL license at https://www.qgis.org. Figure 3 was created through PyGMT (Uieda et al., 2022) and includes data from (Styron & Pagani, 2020) and (Hasterok et al., 2022).