Improved optical image matching time series inversion approach for monitoring dune migration in North Sinai Sand Sea_ Algorithm procedure, application, and validation

Sand dune migration poses a potential threat to desert infrastructure, vegetation, and atmospheric conditions. Capturing the patterns of long-term dune migration is useful for predicting probable desertification issues and wind conditions across vast desert areas. In this study, we employed optical image matching and a singular value decomposition approach to estimate the rates of dune migration in the North Sinai Sand Sea using the free Landsat 8 and Sentinel-2 archives. Our optical image matching time-series selection and inversion (OPTSI) algorithm limited the difference in the solar illumination of correlated pairs to decrease shadows and seasonal variability. We found that the maximum annual dune migration rates were 9.4 m/a and 15.9 m/a for Landsat 8 and Sentinel-2 data, respectively, and the results of time-series analysis revealed the existence of seasonal variations in dune migration controlled by wind regimes. The directions of sand movement extracted from the mean velocity solution agreed strongly with each other and with the drift directions estimated using wind data from meteorological stations. We assessed the uncertainty of each solution based on the variance of stable areas. Our results showed that the proposed inversion decreased uncertainty by up to 25% and increased the spatial coverage by up to 20%. This algorithm is also promising for the retrieval of historical time series on the ground displacements of glaciers and slow-moving landslides employing free archives that provide high-frequency


Introduction
Sand covers approximately 25% of Earth's desert regions, where sand dunes are the most common landforms (Ghadiry et al., 2012), and migratory dunes pose hazards to infrastructure projects in vast deserts (Al-Mutiry et al., 2016;Bruno et al., 2018;Dabboor et al., 2013;Gad, 2016). Over the long-term, sandstorms can affect the health of vegetation and lead to desertification (Ahmady-Birgani et al., 2017;Wang et al., 2003). Dust storms also affect atmospheric conditions, leading to climatic changes that can impact human health (Middleton, 2019;Middleton and Kang, 2017). Therefore, monitoring the migration of sand dunes is important during the planning stages of new cities and the construction of desert infrastructure (Hermas et al., 2012). Moreover, dynamic dune information can be used as a proxy for wind direction, especially in vast deserts where no meteorological records are available (Vermeesch and Leprince, 2012).
Previous studies involving dune migration monitoring have been performed using field measurements and remote sensing techniques (Hermas et al., 2012;Hugenholtz et al., 2012). Field measurements provide high accuracy data; however, they are time-consuming, expensive, and typically focused on a single dune field. Since the 1970s, and especially since Landsat Thematic Mapper images have been captured, remote sensing has greatly contributed to the transition from such field-based studies to studying complex dunes across several aspects, such as boundary conditions, dune patterns, hierarchies, and dune activity (Hugenholtz et al., 2012). Early attempts at remotely sensing dune migration, either through satellite images or aerial photographs, were conducted by measuring the distance between two lines representing the crests of dunes at two different epochs (Hugenholtz et al., 2012). However, the large spatial areas covered by these measurements in order to extract crest lines from satellite images are exposed to different sources of error arising from sensor geometries, topographic shadowing, and different illumination conditions (Hugenholtz et al., 2012). All of these factors can affect the accuracy and reproducibility of classical remote sensing solutions. Additionally, this method has a tendency to result in the overestimation of dune migration (Hugenholtz et al., 2012). Change detection techniques were introduced by Mohamed and Verstraeten (2012), who applied write function memory insertion and red, green, and blue (RGB) clustering methods to monitor dune dynamics at the scale of dune fields and sand seas. These techniques allow for a high spatial coverage but they do not provide a quantitative estimate of dune migration. Some researchers have highlighted the possibility of using high-resolution digital elevation models (DEMs), acquired from light imaging detection and ranging (LIDAR) data, to monitor volumetric changes and dune migration (Dong, 2015;Reitz et al., 2010). However, the high accuracy that is obtainable using LIDAR data, is considered costly, especially when covering large areas.
The coregistration of optically sensed images and correlation (COSI-Corr) (Leprince et al., 2007) provides an alternative approach to monitoring dune migration with less human interaction. Vermeesch and Drake. (2008) first employed such a method to estimate the migration of barchan dunes in the Bodélé Depression using freely available Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images. They reported on the efficacy of the technique to estimate the sand flux during dune migration by integrating horizontal displacements with the dune heights extracted from the stereopairs of the ASTER images. Subsequently, COSI-Corr has been used in many studies to monitor the migration of terrestrial dunes (Al-Ghamdi and Hermas, 2015;Al-Mutiry et al., 2016;Baird et al., 2019;Hermas, 2015;Hermas et al., 2012Hermas et al., , 2019Necsoiu et al., 2009;Sam et al., 2015;Vermeesch and Leprince, 2012) and extraterrestrial dunes (Ayoub et al., 2014(Ayoub et al., , 2014Bridges et al., 2012). Necsoiu et al. (2009) employed a radiometric combination by merging ASTER images with Satellite for the Observation of Earth-5 (SPOT-5) data, which spanned five years, to monitor dune migration in the subarctic Great Kobuk Sand Dunes, Alaska. They proposed the possibility of a radiometric combination of images from different archives and introduced a robust migration direction scheme to better filter noise, especially for slowly migrating dunes. Scheidt and Lancaster (2013) validated the results of COSI-Corr matching using Geographic Information System (GIS) measurements; their study was conducted to monitor dunes in southern Namibia from 2001 to 2009, employing ASTER images. They tested seven pairs with different time separation ranges, from 1 to 9.6 years, and found that pairs with maximum time separations were not effective in retrieving migration rates, especially when the area under consideration had dunes with migration rates greater than the spacing between the dunes. More recently, Baird et al. (2019) investigated the potential of employing the Landsat archives as inputs to the correlation algorithm COSI-Corr, where they used Landsat Level-1 Precision Terrain (L1TP) products, especially Landsat-5 TM to monitor the dune celerity of the fast moving dune filed in the Bodélé Depression, Chad. They reported that two main factors are critical in the selection of Landsat scenes as an input to the correlation algorithm. The first is that the migration rates should exceed the misregistration errors and the second is the distribution of the GCPs should be over stable areas.
Monitoring dune migration over time can be considered as a proxy to the windiness condition over the vast deserts where the metrological observations are rare. Information on the windiness condition can help in addressing the activity of the dust storms. Vermeesch and Leprince (2012) used the historical records of the barchans dune migration as a proxy to monitor the status of the windiness in Sahara Desert. They applied the optical image matching employing seven optical images from different archives (ASTER, and Landsat 4) and reported on the stability of the wind regime in Sahara Desert. Ayoub et al. (2014) introduced a time series to monitor dune dynamics in the Nili Patera Dune Field employing High Resolution Imaging Science Experiment (HiRISE) images and applied weighted principal component analysis (PCA) to decrease noise in the correlated results. They reported on the seasonal variability of the sand flux and suggested a threshold of sand motion to be applied for large scale model wind filed.
Several studies have exploited the redundancy of correlated pairs to decrease uncertainty and enhance spatial coverage using free archives (e.g., Landsat 7/8 or Sentinel-2) that provide images with high revisit times (Altena et al., 2019;Bontemps et al., 2018;Dehecq et al., 2015;Fahnestock et al., 2016;Lacroix et al., 2019). Bontemps et al. (2018) first introduced the inversion of time series in the optical image domain, employing the full redundancy of all possible pairs generated from 16 SPOT-5 images. A simple least-squares inversion of the design matrix was applied owing to the full rank associated with the full redundancy. Their proposed algorithm effectively decreased the uncertainties and enhanced the spatial coverage of the solution. Additionally, they recommended the selection of pairs with similar illuminations to further limit uncertainties. Lacroix et al. (2019) applied the same rationale as Bontemps et al. (2018) by employing Landsat 8 images to study the spatiotemporal dynamics of slow-moving landslides. Their study revealed the existence of seasonal signals associated with time series results; however, the study area was mountainous and vegetation-free and they attributed seasonal variations to differences in the solar illumination of correlated pairs. Altena et al. (2019) used the product of GoLIVE dataset (https://nsidc.org/data/golive) derived from Landsat 8 to monitor the spatiotemporal patterns of the glacier velocity over southern Alaska. They used the least squares inversion to generate time series and reported on the effectiveness of the method in providing an overview about the location and the time of glacier dynamic events.
The inversion algorithm introduced by Bontemps et al. (2018) is considered suitable for commercial satellites with a limited number of images, but it would not be effective when employing free archives that provide images with a high revisit time (e.g., Landsat 8 and Sentinel-2), and the full network combines all pairs without any consideration of the quality of the images. This can affect the results of such time series, especially for pairs with higher variability in solar illumination, leading to seasonal signals, such as those reported by Lacroix et al. (2019). According to our knowledge, there are two other conspicuous gaps in the existing literature: (1) few studies have included the matching of results from either Landsat 8 or Sentinel-2 to monitor dune migration while considering the acquisition geometry, radiometric resolution, and revisit time of these satellites; (2) higher degrees of redundancy have not been employed when establishing dune migration time series, which have instead depended upon adjacent pairing schemes (Ayoub et al., 2014;Vermeesch and Leprince, 2012). Ding et al. (2020) used the free archives of Landsat 8 and Sentinel-2 to study the spatiotemporal patterns of dune migration near Minqin Oasis in northwestern China. The method is considered an extension to the studies related to the inversion of the optical imagery (Altena et al., 2019;Bontemps et al., 2018;Lacroix et al., 2019).
Here, we propose an optical image matching time-series selection and inversion (OPTSI) algorithm. This algorithm can be used to automatically select pairs that are less affected by differences in the illumination and inversion of time-series data. We conduct empirical tests to investigate the relationship between different baselines and uncertainties. A total of 576 pairs, belonging to different tiles, are shared in the baseline test. We conclude that differences in the elevation and azimuth of the Sun are highly correlated with uncertainty, when compared to temporal and spatial baselines. The selection criteria is considered effective for decreasing the time needed to process all pairs, as well as the probable seasonal variations associated with differences in solar illumination. The time series is then inverted using singular value decomposition (SVD). We apply the OPTSI algorithm to monitor dune migration in the North Sinai Sand Sea (NSS) using the free Landsat 8 and Sentinel-2 archives. It is worth noting that we use these two satellites together in order to validate the magnitude and direction of dune migration by comparing both solutions. Additionally, we compare the direction of dune migration to the resultant drift directions extracted from meteorological stations.

Study area
The Sinai Peninsula is an arid environment in which the rate of evaporation exceeds the rate of precipitation. The annual temperature ranges from~14.0°C to~40.0°C (Hermas et al., 2012). Such aridity and high temperatures cause the paucity of vegetation cover observed and provide a suitable environment for dune activity. The NSS is located in the northwestern corner of the Sinai Peninsula (Fig. 1A). It extends in the north-south (N-S) direction from the Mediterranean coast and Lake Bardawil to the side slopes of Um Khushaib Mountain. It extends in the east-west (E-W) direction from the Suez Canal and the Bitter Lakes to the mountains of El-Maghara and Yelleq. The NSS is considered the only sand sea outside the Western Desert and the most complicated dune field in Egypt (Embabi, 2018). All dune patterns and forms are found in the NSS and they are exposed to highly varying wind regimes. The wind speed and direction vary both seasonally and by location. Effective wind speeds > 6 m/s prevail during winter and spring (Embabi, 2018;Hermas et al., 2012). There are two prevailing directions in the effective wind regimes: northwest during summer and spring, and southwest during autumn and winter (ElBanna, 2004;Embabi, 2018;Hermas et al., 2012).

Optical images
We used the optical images from the free archives of Landsat 8 and Sentinel-2 as inputs to the correlation algorithm. All images were downloaded using the bulk downloader application from the United States Geological Survey (USGS) website (https://earthexplorer.usgs. gov/). We used level 1 T images from three orbits of Landsat 8 and level 1C images from one orbit of Sentinel-2 to perform the baseline test (Fig. 1B). After the baseline test, we generated time series of dune movement over part of the NSS using 31 images from orbit (P: 176, R: 39) and 31 images from orbit (TR36VU) (Fig. 1C). The Landsat 8 images cover the period from April of 2013 to March of 2018, while the Sentinel-2 images cover the period from December of 2015 to February of 2019 (Table S1). Both Landsat 8 and Sentinel-2 images were orthorectified using Planet DEM90 (http://www.planetobserver.com) and DTED Level 1 from the National Imagery and Mapping Agency, respectively (Kääb et al., 2016), which were sufficient for our application. Both archives have been enhanced in their radiometric resolutions, where Landsat 8 had 12 bits scaled to 16 bits, compared with the 8 bits of Landsat 7 (Fahnestock et al., 2016). This enhancement allows for the improved monitoring of variations in contrast over bright targets, such as glaciers and dunes (Fahnestock et al., 2016).  (Embabi, 2018;Hereher, 2000). (B) Landsat 8 and Sentinel-2 footprints used for baseline tests. Black boxes show the stable areas of each tile used for estimating uncertainty. (C) Landsat 8 and Sentinel-2 footprints used to study dune movement. The magenta polygon represents the area of interest for Landsat 8, which was used to cut the entire footprint after matching to maintain the same number of rows and columns in all images. The black rectangle represents the stable area for assessing the uncertainty of both solutions. The green polygon is the intersection between the two footprints and was used to cross-validate the two archives. The inset shows the geographic location of the study area with respect to Africa. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Surface wind data
Surface wind data from two meteorological stations in our study area (Bir El-Abd and El-Malease) were acquired from the Egyptian Meteorological Authority. The recorded wind data were provided on a mean annual basis from 2013 to 2015. Wind data were measured at 10 m above the ground and arranged into 12 wind speed classes in 12 directions. For each station, the drift potential (DP), resultant sand drift potential (RDP), and resultant drift direction (RDD) were estimated (Table S2) through vector analysis of the drift estimated for each of the 12 directions (Hereher, 2018(Hereher, , 2010. Data with velocities < 6 m/s were ignored, as such velocities are ineffective for transporting sand (Fryberger, 1979). We used the wind records to generate the sand roses (Fig. S1), and we compared the RDD extracted from wind records to the prevailing directions extracted from optical images.

Methodology
The proposed inversion algorithm included five steps: (1) the automatic selection of pairs, (2) image correlation, (3) the filtering of deformation fields, (4) the inversion of time series, and (5) the estimation of uncertainty and cross-validation of the results (Fig. 2).

Baseline test
Optical image matching results were exposed to different sources of errors, such as topographic shadowing and stripe artifacts, orthorectification and co-registration residuals, and temporal decorrelation stems from surface variations over time. Some of these errors can be suppressed before the correlation process to limit the degree of uncertainty (e.g., orthorectification residuals and the cast shadows). The orthorectification residuals can be mitigated by matching images from the same orbit (Kääb et al., 2016). Also, Lacroix et al.(2018) compared the uncertainties in case of matching images belong to same orbit and neighboring orbits; they found that the uncertainty in the later increased by around 27%.
Noteworthy artifacts stem from large differences in illumination conditions, especially on steep slopes (Bontemps et al., 2018;Ding et al., 2016;Lacroix et al., 2018;Stumpf et al., 2014). These artifacts originate from the interaction between light and three-dimensional (3D) objects; the resulting shadows significantly interact with the topography, thereby affecting the frequencies of satellite images and confusing the matching algorithm (Lacroix et al., 2019). Bontemps et al. (2018) recommended matching images from the same season to control for uncertainties. Meanwhile, Lacroix et al. (2019) reported that seasonal variations were found in their time-series results; however, the area under consideration was a vegetation-free, mountainous region. Such seasonal variations often result from differences in illumination conditions, with rugged topography magnifying these effects. Ding et al. (2016) defined the term "topographic shadowing artifact" (TSA) as a function of differences in the elevation and azimuth of the Sun, and topographic height. They recommended choosing pairs with lower TSAs to decrease the effect of shadows in correlations.
Temporal baselines have also been discussed as one of the dominant factors affecting matching results (Dehecq et al. 2015;Fahnestock et al. 2016). Choosing a suitable time separation is considered a dilemma. On the one hand, it is recommended to select pairs with a short time span to avoid mismatch-induced noise, which can contaminate correlation results (Bontemps et al., 2018;Lacroix et al., 2019). However, Fahnestock et al. (2016) reported that selecting short-span pairs may lead to an increase in the geolocation error of the results, especially in cases where there is insufficient displacement. Consequently, balancing the temporal baseline is critical and depends upon the target under consideration.
Spatial baselines are not commonly used in optical image matching, for which the notion was first introduced by Ding et al. (2016), who defined it as the distance between the geometric centers of the master and slave images. They reported that selecting pairs with spatial baselines of < 200 m is promising for decreasing stripe artifacts, and could be useful as an alternative to the mean subtraction method followed in previous studies (Bontemps et al., 2018;Leprince et al., 2008Leprince et al., , 2007. Previously, Bontemps et al. (2018) reported the relationship between Sun elevation difference (SED) and the uncertainty before and after the inversion. Also, Lacroix et al. (2018) investigated the relationship between the temporal baseline and the uncertainty in different cases such as neighboring orbits, same orbit, with cloud cover and free of clouds. Similarly, in this study, we performed an empirical test to investigate the relationship between different baselines and uncertainties. We exploited that the study area included a large extension of mountains in the south; the geological formation is displayed in Fig. S2. The baselines of the correlated pairs were defined as radiometric, temporal, and spatial baselines and the uncertainty of the deformation field was estimated based on the standard deviation of the stable areas.
Three tiles of Landsat 8 (P: 174, R: 39, P: 175, R: 39, P: 176, R: 39) and one of Sentinel-2 (TR36VU) were included in the baseline test (Fig. 1B). We selected pairs from metadata files by setting the values of TSA and spatial baseline (SB), where the TSA could be estimated as a function of the elevation and azimuth of the Sun in the master and slave images. We selected pairs with TSA < 0.5*h, where h is the topographic height, and SB < 200 m. An exception was tile (P: 174, R: 39), where these values increased to 1*h and 250 m to better investigate the relationship between the radiometric baselines across the rugged topography covered by this tile. Additionally, pairs generated from the adjacent pairing criteria were included in this test, where adjacent pairing was performed by pairing each image with the consecutive one. Table 1 shows the number of pairs, the criteria for their selection, and the area of the stable region.
The batch processing of COSI-Corr was used to correlate all pairs with the parameters and bands discussed in Section 2.4.2. We masked the stable area according to the displayed polygons ( Fig. 1B) and then the records of these stable areas were filtered out, first, by removing the pixels with lower signal-to-noise ratios (SNR < 0.95) in both their E-W and N-S displacements and second, by removing the values outside the range of μ ± 3σ. The expected σ from the correlation was 1/10 of the pixel size (Leprince et al., , 2007Scherler et al., 2008) by assuming μ = 0; values outside the range of ± 5 m were discarded. We estimated the uncertainty of each pair as the standard deviation of the remaining points. Scatter plots of the relationships between the different baselines and the standard deviation were generated. Additionally, the percentage of mismatches, as discussed by Necsoiu et al. (2009), was plotted against the different baselines. Pearson's correlation coefficients (R 2 ) were used to represent the strengths of the relationships between each baseline and either the uncertainty or the level of mismatch.

Pairing criteria
Our selection algorithm was mainly oriented to select high-quality pairs to reduce error and time costs. We constructed our selection algorithm on the following five decisive factors: (a) Limiting the error derived from orthorectification residuals, where these residuals stem from the vertical errors in the DEM that translate to horizontal displacements, especially in pairings from neighboring tiles (Kääb et al., 2016). Previously, Lacroix et al. (2018) compared the levels of uncertanty when mtaching images in cases of the same and neighboring orbits and found that the uncertanity in the later increased by~27%. Consequently, we matched pairs from the same tile in our study. (b) Limiting the cloud cover associated with correlated pairs. Lacroix et al. (2018) reported on the effects of cloud cover in correlated images on the resultant uncertainties, and found that they reached a maximum of~4 m. Therefore, we limited the cloud cover to < 1% to decrease the level of mismatch, especially over moving targets and control the uncertainty. (c) Limiting the differences in solar illumination between correlated pairs. The difference in the illumination between correlated pairs causes shadows and seasonal variations. In previous studies, researchers have recommended the selection of pairs with little variation in illumination conditions (Bontemps et al., 2018;Ding et al., 2016;Lacroix et al., 2019;Stumpf et al., 2014) to decrease the associated artifacts. Our baseline test results (see Section 3.1) also supported limiting the differences in the angle of the Sun to further limit uncertainty and levels of mismatch. (d) Limiting the time separation. The selection of a suitable time separation for the correlated images is a problem. The processing of short-and long-time separations was introduced by Bontemps et al.
(2018), who used a full network to establish time series. Here, we mixed the processing of short-and long-time separations and limited the maximum separation based on the maximum displacement over the study area and the window size used for correlations. It has Table 1 Information on the baseline test for the different case studies, including the number of pairs shared in the test of different tiles, area of the stable region, method of pairing selection, and Pearson's correlation coefficient (R 2 ) between the different baselines and uncertainties in both directions (E-W and N-S). For more details about the calculation of TSA, please refer to (Ding et al., 2016). been reported that the migration rates over our study area range from 1 to 27 m/a (Hermas et al., 2012). (e) Ensuring connectivity of the pairing network. According to the single baseline subset (SBAS) method applied in interferometric synthetic aperture radar (InSAR) studies, the selection approach always makes the images occur in different subsets (Berardino et al., 2002). Consequently, the conditions of the design matrix are affected, which leads to the rank deficiency of the design matrix; the inversion is no longer valid using the least-square scene (Berardino et al., 2002). Therefore, in our selection, we focused on ensuring a better condition for the design matrix and less oscillation of the SVD solution by limiting the number of subsets (Reinisch et al., 2017).
The selection procedure was applied before downloading the images. We first downloaded the metadata file of the tile under consideration. The cloud cover was then set to 1% and the SED, SAD, TB, and SB were estimated for all possible pairs. Finally, the thresholds were set iteratively to limit the Sun elevation difference (SED), Sun azimuth difference (SAD), TB, and the number of subsets. The thresholds of the selection algorithm were determined after several trials. The SED and SAD were set to less than 10°and 8°for Landsat 8 and Sentinel-2, respectively, while the maximum time separation was set as less than 4.5 and 3.05 years for Landsat 8 and Sentinel-2. The SB threshold was determined flexibly to adjust the network when we found that there was no considerable impact of the SB on uncertainty. The SB was set as less than 2500 m for Landsat 8, while no limitation was set for Sentinel-2. Consequently, 90 and 73 pairs were selected for establishing the time series for Landsat 8 and Sentinel-2, respectively. Information on pairs can be found in (Tables S3 and S4).

Image correlation
The selected pairs were correlated using the frequency correlator of COSI-Corr (Leprince et al., 2007), which retrieves the horizontal displacement from the phase shift in the low-frequency content, wherein the software estimates the displacement iteratively in two steps by applying a variable window size. A large-scale window size maximizes the correlation between two images, while a small window size provides a finer estimation of ground displacement (Lacroix et al., 2019). We employed the panchromatic band 8 (15 m) of Landsat 8 as an input to the correlation algorithm to better capture the details of the sand dunes. Sentinel-2 provides a 10-m resolution for the blue (B2), green (B3), red (B4), and infrared (B8a) bands. We used the near-infrared band for correlations after testing and comparing the performance of the four bands (see Supporting Information).
A batch correlator was used to automate the correlation process. To perform the correlation, several parameters must be selected, including the initial and final window sizes, the step size, the number of iterations, and the mask threshold. The initial and final window sizes were selected after performing a sensitivity analysis to investigate the relationship between the window size and uncertainty (see Supporting Information). We selected a 128 × 64 window size for Landsat 8 and a 64 × 32 window size for Sentinel-2.
Step sizes of 4 × 4 pixels for Landsat 8 and 6 × 6 for Sentinel-2 were selected to ensure similar ground resolutions of the generated maps. A frequency mask was applied to ensure the selection of the locations of the cross-spectrum, where the phase information was valid (Sun et al., 2017). The mask threshold was set to 0.9, and four robustness iterations were used.
Correlations always produced three maps: E-W displacement, N-S displacement, and SNR maps. The first two maps represented the displacement in the east-west and north-south directions, respectively, where positive signs represented displacement toward the east and north. The SNR map represented a function to evaluate the quality of the correlation, with values ranging from 0 to 1; higher values were assumed to represent better matches. The definition of the SNR produced by COSI-Corr was different from the classic definition, where it can be estimated as a function of the normalized cross-spectrum and a weighting matrix for different frequencies (Bontemps et al., 2018).

Filtering of displacement fields
We performed the following four steps to mitigate the effect of the associated errors originating from orthorectification residuals, co-registration, and the misalignment of charge couple device (CCD) arrays (Bontemps et al., 2018;Ding et al., 2016;Scherler et al., 2008): (a) We discarded the potential outliers by removing pixels with lower SNRs from both displacement fields (E-W and N-S directions). The determination of a suitable SNR threshold is critical for balancing the real deformation while discarding correlation noise. To balance between removing the mismatches and the spatial coverage, we set the SNR threshold to 0.9. It is worth noting that removing the pixels with lower SNR values may inadvertently eliminate the correct displacement, whereas retaining those with higher SNR values will retain mismatches (Paul et al., 2017). Then, we removed the outliers from the displacement fields outside the range of ± 20 m. (b) The second step involved removing the linear ramps by applying a polynomial fit to the stable areas and then subtracting the generated plane from the raw results. (c) We then removed the associated stripes and jitter artifacts from the deformation fields. These stripes were mainly caused by the misalignment of the CCD arrays of the sensor (Ding et al., 2016;Leprince et al., 2008Leprince et al., , 2007Scherler et al., 2008). This step was performed through the following procedures (Stumpf et al., 2018): (1) determination of the azimuth angle of the stripes measured from the north, (2) rotation of all of the deformation fields to be in vertical alignment, (3) masking of the stable areas, (4) estimation of the median value of each column in the stable area (Bontemps et al., 2018), (5) subtraction of the value of each pixel in the column from the median value of each column, and (6) rotation of the deformation field back to the original location. (d) Finally, the remaining values were denoised using a non-local mean filter (Ayoub et al., 2009).

Time-series inversion
If we have (N + 1) images belonging to the same tile spanning a period from t t [ , , ] N 0 , where t 0 is the reference starting date, the selection of pairs will always lead to the location of the images in different subsets, where the number of pairs (M) can be estimated according to the following relationship: (Berardino et al., 2002) The incremental displacement between successive dates can be described as: (2) and the relationship between the observed matching measurements ( d tj ) and the displacements (d ti ) can be described according to the following equation: where the two vectors, IE and IS, are the time indices, as each pair includes the master image IE ( ) and slave image IS ( ). It is worth noting that master and slave are ordered chronologically: > IE IS and = j M 1, , (Berardino et al., 2002) Consequently, the relationship between the set of M observations and the displacements can be described as: Therefore, Berardino et al. (2002) suggested the calculation of the E. Ali, et al. ISPRS Journal of Photogrammetry and Remote Sensing 164 (2020) 106-124 mean velocity of displacement between each adjacent time epoch. The mean velocity can be described as follows: From Eq.
(2), the relationship between the mean velocity v i , and the time epochs can be described as follows: Rearranging Eq. (7) to be in matrix form: where Matrix B is an incidence matrix with the dimensions of [M × N], such that the elements in the matrix are i IE 1 j j and otherwise, 0. Matrix B is similar to an incidence matrix, wherein the conditions of the matrix vary according to the configuration of the selected pairs (Berardino et al., 2002). Design matrix B has two conditions based on the number of subsets of image acquisitions. First, in the case of all acquisitions belonging to one subset, B will be an N-rank matrix with a well (M = N) or over-determined condition (M > N) (Berardino et al., 2002). In this case, the inversion of matrix B can be obtained in the least squares sense as follows: Generally, the selection procedure affects the temporal acquisition, which usually do not belong to the same subset. Accordingly, matrix B is exposed to a rank deficiency, and the term (B B) T 1 is a singular matrix. The relationship between the rank and number of subsets can be described as follows (Berardino et al., 2002): where L is the number of subsets.
The pseudo-inverse of the rank-deficient matrix can be obtained by SVD, and this inverse provides a minimum normal least-squares solution. The decomposition of matrix B using SVD is as follows (Berardino et al., 2002): Accordingly, the solution of Eq. (8) in the case of a rank deficiency can be represented as follows: To better consider the effects of decorrelation with time, we applied the weight, W of each pair based on the time separation between the master image IE ( ) and the slave image IS ( ). This weight can be estimated as follows (Lacroix et al., 2019): Incremental displacements were acquired by multiplying the mean velocity between each epoch by the time separation. The inversion was performed pixel-by-pixel in both the E-W and N-S directions. Because of the filtering process, some pixels were discarded, leading to partially stacked pixels. To balance the redundancy level and spatial coverage, these pixels were included in the inversion if the percentage of the remaining pixels in the time series was larger than 70% and 75% of the full stack for Landsat 8 and Sentinel-2 data, respectively.

Uncertainty estimation
Previously, researchers have evaluated the uncertainty of matching results depending on their variance over stable areas (e.g., Bontemps et al., 2018;Dehecq et al., 2015;Kääb et al., 2016;Lacroix et al., 2018). Lacroix et al. (2019) reported that such evaluations of the uncertainty over the stable area may be not representative of the uncertainty in all regions of the image, especially in cases of the heterogeneous mass with small moving targets. However, the stable area would be useful in case of the paucity of in situ measurements to assess the uncertainty of the moving targets. Dehecq et al. (2015) evaluated the uncertainty of moving glaciers and found that, over the moving targets, the uncertainty was almost twice that over the stable area, as stated by Lacroix et al. (2019). We estimated the uncertainty in the deformation fields based on the standard deviation of stable areas. A stable area of 1086 km 2 was selected to be fully covered by the footprints of both Landsat 8 and Sentinel-2 (black rectangle, Fig. 1C). We compared the uncertainties of individual pairs before and after inversion following the same rationale of Bontemps et al. (2018) and estimated the uncertainty of each pair in three stages: (1) raw pairs, (2) after the filtering process, and (3) after inversion. We used the points that came from the inversion of the full stack and then estimated the percentage of improvement in the uncertainty thanks to the filtering and the inversion. To explore the effects of the inversion algorithm and the redundancy level on uncertainties, we independently compared the cumulative displacements from both the inversion and non-inversion solutions. Both solutions provided displacements at the same time epochs, wherein the latter was performed by matching consecutive images. The non-inversion results were filtered by applying the procedure described in Section 2.4.3. To ensure equality in the comparison, we masked out a similar stable region (Fig. 1C) and then filtered out the records over the stable area by discarding values outside the range ± 5.0 m. Due to the filtering process, the number of remaining points varied between dates, which could affect the comparison depending on the standard deviation. Additionally, the standard deviation would overestimate the uncertainty, as it does not account for the reduction of error due to averaging over larger regions (Bolch et al., 2011). The standard error (SE) can be used as an alternative to estimate uncertainties. According to the rules of error propagation, we estimated the uncertainty as a function of the SE, and the mean of the displacement over stable areas (MED). The error of displacement was estimated as (Sun et al., 2017): where N eff is the number of independent measurements, PS is the pixel size, and D is the distance of spatial autocorrelation, where D is approximately 20 times D (Bolch et al., 2011)

Post-processing
After the filtering process, further smoothing of the results can be useful for decreasing the effects of seasonal signals encountered in the results of time series. The incremental time-series maps of both solutions were filtered by applying the direction filter introduced by Necsoiu et al. (2009). The robust migration direction of each date included negative values, wherein the noise level exceeded the signal of dune movement; we discarded the values of migration at these locations. We introduced two approaches to analyze the time series of dune movement. The first was the cumulative mobility of dunes (Section 3.3.1), wherein we estimated the total cumulative movement without considering the direction of movement. This approach can be considered to reflect the total energy of the wind triggered the dune mobility from the start to the end date. The second approach involved the cumulative migration of dunes along the prevailing direction (see Section 3.3.2), wherein we assumed that the dune only had one direction of movement. The mean direction of movement was extracted for each pixel by applying the azimuth equations after estimating the E-W and N-S mean velocity maps. After estimating the prevailing direction, we projected the direction of each date into the prevailing direction by estimating the differences between the migration direction of each date and the prevailing direction. The magnitude of movement was multiplied by the cosine of the difference between these two directions. The second approach can be considered to reflect the work done by the wind to force dunes to move along the prevailing direction. We analyzed time-series trends by capturing dune migration over a group of points belonging to different morphologies (Fig. 5C). We then applied a statistical filter to remove outliers, depending on the direction of dune movement. The time series of both solutions (i.e., Landsat 8 and Sentinel-2) were compared in terms of the degree of correlation and standard deviation of the difference between the migration values at the common epochs. Table 1 shows the Pearson's correlation coefficients between the uncertainties and the different baselines. For Landsat 8, the radiometric baselines, especially the SED, displayed a positive linear trend for all cases in the E-W (Fig. 3A1-B4) and N-S (Fig. 3C1-D4) directions. The correlation coefficients ranged between 0.23 and 0.79 in the E-W direction and between 0.18 and 0.94 in the N-S direction. A similar trend for SAD was observed, but with a lower degree of correlation, where the correlation coefficients ranged between 0.24 and 0.62 in the E-W and 0.29 and 0.91 in the N-S directions. For the Sentinel-2 test, we observed that the radiometric baseline components, as well as the TB, did not record any positive linear relationships (Fig. S3). It is worth noting that observations in the N-S direction displayed a higher degree of correlation with the radiometric baselines than in the E-W direction. These significant differences can be attributed to the orbit of Landsat 8 being ordinated along the N-S direction in the study area as reported by Lacroix et al. (2019).

Baseline test results
For Landsat 8, the TB of the selected pairs ranged from 16 d to 1600 d. A slight positive relationship was observed in the E-W direction (Fig. 4A1-A3), while a negative correlation was observed in the N-S direction (Fig. 4B1-B3). The correlation coefficients between the TB and uncertainty ranged between 0.26 and 0.80 in the E-W direction and between 0.00 and 0.80 in the N-S direction, where the relationship can be better described by a cyclical waveform with a frequency of 1 year. It is interesting to note that higher correlations between the TB and uncertainty were recorded for the adjacent pairings of the Landsat 8 tile (P: 174, R: 39), where the time separation ranged from 16 to 128 d. The TB and uncertainty showed a positive linear relationship with correlation coefficients of 0.80 and 0.64 in the E-W and N-S directions, respectively (Fig. S4). These strong linear relationships may support restricting the time separation to limit the uncertainty, as suggested by Lacroix et al. (2018). However, small time separations lead to increases in geolocation errors in the displacement field, especially for slowmoving targets (Fahnestock et al., 2016).
We found that pairs with the same time separation recorded different values of uncertainties, which agrees with the findings of Lacroix et al. (2018). A detailed investigation of these pairs revealed that for smaller TBs (16 d), both SED and SAD were lower than 10°and the uncertainties reached a maximum of 0.75 m in both directions at the maximum values of SED and SAD. There was no clear linear trend between either SED or SAD and TB, where the relationship can be better described by a cyclical variation (Fig. S5). We further explored the relationship between the radiometric baselines (SED and SAD) and the uncertainties after imposing the selection criteria. The uncertainties of both SED and SAD lost a positive linear relationship. However, a slightly increasing linear trend was observed in the relationship between TB and uncertainty (Figs. S6 and S7). It can be concluded that when the differences in solar illumination were limited, the TB began to control the uncertainties, which would support applying pair-based weight according to the time separation of each pair to better include the pairs in the time series as recommended by Lacroix et al. (2019). The relationship between the mismatches and baselines revealed that the SED and SAD were < 10°, leading to mismatch levels of < 10% (Fig. S8), which would support limiting the values of SED and SAD to control the mismatches level and enhancing the spatial coverage. The relationship between SB and uncertainty is shown in Fig. 4C1-D3, wherein SB shows the weakest relationship among the four baselines. However limiting the spatial baseline to smaller values (< 200 m) as recommended by Ding et al. (2016) would be useful in decreasing the stripes artifact comes from the misalignment of CCD arrays, especially for tiles without sufficient extension of stable targets to perform the stripes artifact corrections.

Spatial analysis of dune patterns
The NSS comprises different dune morphologies, including barchan, linear, and transverse dunes. Linear dunes are the most common landform in the eastern part of the sand sea near the Bitter Lakes. To the east of the mountainous region, the wind regime is bidirectional and both barchan and transverse dunes are present (Hermas et al., 2012). Previous studies (Embabi, 2018;Hermas et al., 2012) have shown the variability of linear and barchan dune migration over the NSS. Linear dune migrations range from 0.7 m/a to 27 m/a, while those of barchan dunes range from 3.6 m/a to 17.3 m/a (Hermas et al., 2012). These large variabilities in dune migration are attributed to the different techniques used for capturing them, as well as the dune type, dune size, duration of measurement, and study period (Hermas et al., 2012).
Long-term annual rates of dune migration were extracted from Landsat 8 over the period from April of 2013 to March of 2018 and from Sentinel-2 from February of 2015 to December of 2019. Fig. 5 shows the spatial variability of the mean annual rates extracted from both solutions, where the migration rates ranged from 0.5 to 9.4 m/a for Landsat 8 and from 0.5 to 15.6 m/a for Sentinel-2. The lower limit was set to be 0.5 m/a, as we discarded pixels showing changes of less than 0.5 m/a and considered them as interdune areas, which were defined as areas that do not experience significant migration rates. We extracted the migration rates of each dune morphology based on the polygons displayed in Fig. 5C. The mean net migrations of the linear dunes were 1.57 m/a and 1.63 m/a for Landsat 8 and Sentinel-2, respectively, while for barchan dunes, they were 1.22 m/a and 1.36 m/a. Maximum migration rates were captured in the southern part of the study area (white box in Fig. 5), where linear dunes had mean net migrations of 4.14 m/a and 4.44 m/a for Landsat 8 and Sentinel-2, respectively. We extracted the mean migration direction for each pixel by applying the azimuth equations to the mean annual rates for both the E-W and N-S directions. It can be seen that the migration directions were mainly toward the east and northeast. An exception was seen in the area enclosed by the white rectangle in Fig. 5, where the dune migration direction was toward the south; this area included linear dunes arranged parallel to the N-S direction. The prevailing directions are consistent with previous reports on prevailing dune migration directions in this region (Embabi, 2018;Hermas et al., 2012) Fig. 6 shows a comparison between the cumulative displacement captured by Landsat 8 and Sentinel-2 at different locations (see also Fig. 5C). The time series of both solutions were interpolated in the overlapping periods to obtain the values of dune migration at the overlapping time epochs. The Sentinel-2 epochs were then shifted upward with a value equal to the perpendicular distance between the two linear fits of both solutions (magenta and green lines). We investigated the standard deviations of the differences between both solutions, and the values ranged between 1.94 m and 2.88 m. Moreover, the correlation coefficients between the overlapping epochs for all cases exceeded 0.97, which indicates a strong agreement between the trends captured by both solutions. This robust correlation reveals that both solutions captured the same trends of dune migration. Through a detailed investigation of the time series, we found that both solutions experienced variations in the trend of dune movement over time, where the trend varied between being stable, sharply increasing, moderately increasing, and decreasing. This variability in the temporal trend is also supported by previous reports of large variability in wind speed and direction (Hermas et al., 2012).

Cumulative dune mobility
The slope of the linear fit extracted by the least squares method can be used to describe the yearly movement rate at each locations. However, it cannot be considered as the migration rate of dunes, as this approach does not consider the variability of the direction of dune movement through the time series. The rates of movement ranged from 4.79 to 12.95 m/a for Landsat 8 and between 8.36 and 13.25 m/a for Sentinel-2. It is worth noting that the comparison between both solutions was exposed to several sources of noise, including the difference in the performance of the correlation algorithm when correlating Landsat 8 and Sentinel-2 because of their spatial and radiometric resolutions, the performance of the inversion algorithm owing to differences in the redundancy levels of the pairing networks between Landsat 8 (90/30) and Sentinel-2 (73/30), and the complex wind regime of dune migration in the study area. Fig. 7 shows the cumulative time series of the same dune morphologies as in Fig. 6, but considers one direction to represent dune migration. The linear fit of each solution was estimated using a least squares approach, wherein the migration rates along the prevailing direction ranged from 1.97 to 2.92 m/a for Landsat 8 and from 1.63 to 3.30 m/a for Sentinel-2. It is worth noting that negative values indicated the epochs in which the migration direction was opposite to the prevailing direction of movement. This occurred when the difference between the migration angle and the prevailing migration angle ranged between 90°and 270°. We also observed significant differences between the cumulative migration with and without consideration of the direction due to variability in the migration direction between different epochs. The projected cumulative displacements along the prevailing movement direction can be assumed as the work done by wind to trigger barchans dunes to migrate and linear dunes to elongate, where the dune dynamic patterns are mainly controlled by the prevailing wind direction. Furthermore, seasonal signals were still encountered in the results of the time series, regardless of the selection criteria, inversion method, and post-processing technique. These seasonal signals were assumed to express the nature of the dune patterns in the study area, and this assumption is supported by previous reports on the variability of the wind regime in northern Sinai (Hermas et al., 2012). The prevailing wind directions in the NSS are to the northwest, north, and northeast, and the wind direction also varies seasonally.  E. Ali, et al. ISPRS Journal of Photogrammetry and Remote Sensing 164 (2020) 106-124 The study area is exposed to pronounced variations in the percentage of the effective wind regime with different seasons, which is considered effective for transferring sand particles (Hermas et al., 2012). The temporal pattern of the dune migration recorded over different morphologies is useful for interpreting the conditions of sand dunes over time and hence in predicting climatic changes and atmospheric conditions; however, this should be validated using the monthly values of the RDP/DP extracted from wind data. Unfortunately, such monthly values are not available for interpretation in the study area. Additionally, medium-resolution matching results of dune migration may be insufficient, especially when focusing on a single dune field.

Pair-based uncertainty
Figs. 8 and 9 show comparisons between the levels of uncertainty before filtering, after filtering, and after inversion for Landsat 8 and Sentinel-2, respectively. It can be seen that the uncertainties for both Landsat 8 and Sentinel-2 raw pairs after removing the outliers were within the range of the COSI-Corr accuracy at 1σ, which ranges from 1/ 10 to 1/5 of the pixel size (Leprince et al., , 2007. However, some pairs had much higher uncertainties that exceeded 1.5 m in both directions for Landsat 8. A similar finding was observed for some Sentinel-2 pairs, especially in the N-S direction. The filtering process reduced the uncertainties on average by 25% and 24% in both directions for Landsat 8 and Sentinel-2, respectively. The uncertainties after inversion ranged from 0.32 to 0.99 m in the E-W direction and from 0.33 to 1.02 m in the N-S direction for Landsat 8. For Sentinel-2, they ranged from 0.20 to 0.60 m in the E-W direction and from 0.40 to 1.50 m in the N-S direction. Inversion decreased the uncertainties by 20% on average in both directions for Landsat 8 and by 33% and 27% in the E-W and N-S directions, respectively, for Sentinel-2. 3.4.2.1. Cumulative solution uncertainty. Fig. 10 shows the comparison between the inverted and non-inverted solutions in terms of the root mean squared error (RMSE). For Landsat 8, the RMSEs were 3.2 m and 4.1 m in the E-W and N-S directions, respectively, for the inverted solution, while for the non-inverted solution, the RMSEs were 5.0 m and 5.7 m. For Sentinel-2, the RMSEs were 2.24 m and 2.75 m in the E-W and N-S directions, respectively, for the inverted solution, while for the non-inverted solution, they were 2.94 m and 3.20 m. Further investigations were conducted depending upon the error of displacement in which, for Sentinel-2, the mean errors of displacement over the entire time series were 0.19 m and 0.82 m for the inverted and non-inverted solutions, respectively. It is worth noting that displacement error is considered efficient, as it takes into consideration the number of points, the mean of the records, and the standard deviation. The dependence on the standard deviation is considered ineffective, especially for data with a skewed distribution, as the records are no longer distributed normally around a mean.

Mean velocity uncertainty.
We assessed the uncertainty of the mean velocity solution based on the standard deviation of the stable areas (Fig. 1C). The histograms of the noise over the stable area are shown in Fig. 11D and E. It is notable that the noise histogram displays an approximately Gaussian distribution and does not have a Rayleighlike shape; this provides evidence that the mean velocity solution is an unbiased measurement, as reported by Necsoiu et al. (2009). The mean of the records was almost zero, with a standard deviation of 0.15 m/a and 0.13 m/a in the E-W and 0.13 m/a and 0.12 m/a in the N-S directions for Landsat 8 and Sentinel-2, respectively. It is clear that the recorded uncertainties with the inversion of one rate were less than the uncertainties in the case of an incremental solution for both Landsat 8 and Sentinel-2, where the average standard deviation over all dates records values of 0.63 m and 0.36 m for Landsat 8 and Sentinel-2, respectively. This supports the previously reported effects of redundancy on uncertainty. These findings also agree with those of Dehecq et al. (2015), who reported that increasing the number of pairs leads to decreasing the uncertainties.

Validation
We compared the long-term mean migration from Landsat 8 (15 m) and Sentinel-2 (10 m) to validate our results. Significant spatial variability among the migration rates existed in the study area (Hermas et al., 2012); the range of variation captured in previous studies was 1-27 m/a. Such a wide range of migration rates can be attributed to the variability of the wind speed and direction, duration of the measurements (i.e., the capture interval), and time interval over which the measurements were performed (i.e., the study period) (Hermas et al., 2012). Nevertheless, it is interesting that similar trends of dune migration rates and patterns were observed, with an increase in the rates recorded in the Sentinel-2 solution (Fig. 11). The differences between the two solutions can be attributed to the different periods used to generate their mean velocities. A cross-comparison between them at common points (see intersection in Fig. 1C) reveals a positive linear relationship with correlations of 0.80, 0.89, and 0.85 for the E-W direction, N-S direction, and magnitude, respectively (Fig. 11). Additionally, the RMSEs of the differences between the results of both solutions were 0.44, 0.59, and 0.57 m/a for the E-W direction, N-S direction, and magnitude, respectively.
In terms of the direction of movement, the histogram shown in Fig. 11F reveals the frequency of the movement of both solutions. The direction of movement displayed a bimodal distribution with two peaks: the first at~55°and the second at~155°. These agree with the spatial distribution of the direction of movement, wherein the major direction was toward the east and northeast in all locations. An exception was observed in the southern NSS (white rectangle in Fig. 5C), where the direction of movement was toward the southeast.
To further investigate the spatial variability and the agreement between the two solutions, the migration rates across six transects at different locations (Fig. 5C) were extracted and compared. Fig. 12 illustrates the migration rates along different transects from P1-P6. These transects differed in length and passed through different dune morphologies. It can be seen that the location of the interdune areas (low migration rates) and areas of dune crests (high migration rates) can be identified clearly for profiles passing through the barchan and linear dunes. It is also interesting to note that almost the same trend of dune migration was captured; however, the rates of migration recorded by the Sentinel-2 solution were higher than those of Landsat 8. The correlation coefficients between both solutions ranged from 0.68 to 0.88, wherein the lowest correlation recorded for the profiles passed through barchan dunes (P1, P3, and P4). In contrast, transects passing through linear dunes indicated a good correlation (P2 and P6). This comparison revealed that both archives can be used to retrieve the magnitudes and directions of dune migration with high reliability and reproducibility. Furthermore, the coarse resolution provided by these archives should be sufficient for studies at the dune-field-scale. However, narrowing the scope of study to a single dune field (e.g., to study the physics of the dune) may require higher resolution images.

Performance of the proposed algorithm
Our method is an improvement of the inversion of the optical image matching techniques introduced in previous studies (Bontemps et al., 2018;Lacroix et al., 2019). In our algorithm, we first selected pairs automatically to take advantage of the processing of redundant pairs to reduce uncertainty and enhance spatial coverage (Dehecq et al., 2015). The selection criteria were ordinated mainly to decrease the difference in the solar illumination between pairs. Limiting the difference in the angles of the Sun should be useful for decreasing shadows and seasonal signals (Bontemps et al., 2018;Lacroix et al., 2019); additionally, our baseline test revealed a strong relationship between solar illumination and uncertainty. Because of the selection, the inversion was no longer valid using the least squares method, similar to previous studies in which full networks were established (Altena et al., 2019;Bontemps et al., 2018;Lacroix et al., 2019). Therefore, an SVD inversion was used in our study.
The contribution of our algorithm was assessed by investigating the effect of the algorithm on spatial coverage and uncertainty (Bontemps et al., 2018). First, we explored the effect of our algorithm on enhancing the spatial coverage by first comparing the mean velocities from all pairs before and after inversion (Figs. S11 and S12) and then comparing the cumulative displacements extracted with and without inversion (Figs. S13 and S14). Enhancement in the spatial coverage was estimated by calculating the percentage of the non-masked pixels before and after inversion (Table 2). This comparison revealed that inversion decreased the uncovered area when averaging all pairs by 16% and 25% for Landsat 8 and Sentinel-2, respectively. The uncovered area decreased by an average of 12% for Landsat 8 and Sentinel-2 when comparing the cumulative displacements of the inverted and non-inverted solutions. It is worth noting that the algorithm introduced in the study of Bontemps et al. (2018) significantly reduced the gaps by 14% when averaging 240  Fig. 1C) for both Landsat 8 (blue) and Sentinel-2 (orange) in the E-W and N-S direction, respectively. Panel (F) displays the histograms of the direction of dune movement for Landsat 8 (blue) and Sentinel-2 (orange) in the overlap area (see Fig. 1C). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) pairs of SPOT-5 data and compared before and after the inversion. Our proposed algorithm significantly enhances the spatial coverage when compared to this previous algorithm. It is worth noting that ensuring good spatial coverage is critical to enable capturing the whole picture of the spatiotemporal patterns of the surface displacements.
We examined the effects of the proposed algorithm on uncertainties. The algorithm helped decrease the uncertainty of the individual pairs tõ 7% of the pixel size for both Landsat 8 and Sentinel-2 data. Based on the inversion solution, the maximum uncertainties recorded were approximately 10% and 7% of the pixel resolutions for Landsat 8 and Sentinel-2, respectively. Bontemps et al. (2018) reported that uncertainty levels of approximately 10-20% of the pixel size were recorded in previous studies (Heid and Kääb, 2012;Michel and Avouac, 2002;Scherler et al., 2008;Taylor et al., 2008) employing a mediumresolution of 10-15 m when applied over different time periods ranging from 1 month to 3 years. However, the full-rank inversion performed by Bontemps et al. (2018) achieved uncertainty levels of approximately 16-20% of the pixel size when employing a 10-m resolution to monitor deformation over 35 years. In summary, our algorithm effectively decreases the levels of uncertainty and generally outperforms those proposed in other studies in which medium resolutions were employed. Additionally, our method was applied to monitor deformation over 5.5 years.

Prevailing migration direction from mean velocity vs. resultant drift direction
The wind regime in North Sinai is complex because of the variability in both the spatial and temporal domains (Embabi, 2018;Hermas et al., 2012). In previous studies, researchers have analyzed sand movement by estimating the sand drift based on wind data records from meteorological stations (Hereher, 2014(Hereher, , 2010Philip et al., 2004). To determine the reliability of our solution in extracting the direction of dune migration, we compared the extracted prevailing migration direction to the RDD computed from the wind data of two meteorological stations, at Bir El-Abd and El-Malease. Additionally, we compared our results with the RDD reported in previous studies (ElBanna, 2004;Hereher, 2000Hereher, , 2014Philip et al., 2004) (Fig. 13).
The mean annual records of the two meteorological stations were available from 2013 to 2015. We estimated the potential sand drift according to Fryberger's (1979) equation, where wind speeds > 6 m/s were considered effective for transporting sand particles. As shown in Fig. 13, dune migration direction ranged from 54 to 160°for Bir El-Abd and 46-166°for El-Malease. This wide range may be attributed to the spatial variability of the wind speed and the length of the measurement period. Despite this high variability, previous studies (Embabi, 2018;Hermas et al., 2012) have reported that sand in the NSS tends to move toward the east and southeast. Our mean velocity results agreed strongly with such findings, as shown in Fig. 5. Here, the migration was almost eastward in the majority of locations, and toward the south in Fig. 12. Comparison between the spatial variability of the migration rates at different locations and through different transects (see Fig. 5C). The migration rates for Sentinel-2 are in red and Landsat 8 are in blue. P2 and P6 pass through a group of linear dunes. P1, P3, and P4 pass through a group of barchan dunes, and P5 passes along one of the linear dunes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2
Comparison between the spatial coverage of two cases: (1) between the mean migration rates extracted from all pairs before and after inversion, and (2) between the cumulative displacment extracted from the conventional and inverted solutions. the southern part of the study area, where the prevailing dune morphology was longitudinal and structured from north to south. For a better comparison between the prevailing direction and the RDD estimated from the wind data, we extracted small areas at the locations of the two meteorological stations. Then, we statistically analyzed the data to identify and discard the outliers. The prevailing directions of the Landsat 8 mean velocity solution recorded 96.4°± 12°and 102°± 15.6°at Bir El-Abd and El-Malease, respectively. The RDD extracted from recent wind data ranged from 59 to 98°, with a mean value of 74°± 18.2°, at Bir El-Abd and from 83 to 102°, with a mean value of 95.4°± 8.6°, at El-Malease. The standard deviations of the directions at Bir El-Abd were more significant than those at El-Malease, which may be considered as an indicator of greater local variability in the wind regime. It can be observed that good agreement between the values extracted from the prevailing migration direction and the average RDD, especially for Bir El-Abd, where the differnces between the prevailing direction and the mean RDD recorded 22°and 7°, at Bir El-Abd and El-Malease, respectively. Theses difference are considered acceptable, where the two values were extracted from different methods, and over different time interval (i.e., the study period). To further investigate the potential of the mean velocity for estimating a reliable prevailing directions of movement, we compared the sand roses of each dune morphology extracted from both solutions (i.e. Landsat 8 and Sentinel-2), as shown in Fig. 14. Generally, there was a strong agreement between both solutions, especially for linear dunes and sand sheets. A quantitative analysis of the prevailing movement directions revealed that the maximum difference between the two solutions as 9°, which is considered promising, especially in this complex wind regime. Fryberger (1979) used the ratio between DP and RDP to express the directional variability of the prevailing wind direction in a field. The ratio of RDP/DP can be separated into three classes: the first is when RDP/DP is high (> 0.7), which means that a unidirectional wind regime exists and the transverse dune is dominant, the second is when (0.7 ≤ RDP/DP ≤ 0.3), which occurs when a bidirectional wind regime prevails and linear dunes occur, and the third is when the ratio (RDP/ DP < 0.3) occurs for regions with multidirectional wind regimes, in which star dues are common (Hereher, 2018). The three-year average RDP/DP values were 0.53 and 0.58 at Bir El-Abd and El-Malease, respectively, which considered second class with linear dune prevailing dune morphology. Together, these findings demonstrate that the prevailing direction derived from the mean velocity solution can effectively represent prevailing wind direction, especially when there is a paucity of metrological stations, such as in vast deserts (e.g., the Sahara and Arabian deserts).

A comparison with previous studies
Previous studies have been conducted on dune migration over different parts of the NSS. Hermas et al. (2012) reported on the migration rates extracted from previous studies. It is worth noting that in most previous studies, field measurements or classical remote sensing techniques were employed to focus on either the linear or barchan dunes. The migration rates have exhibited a wide range of variation, which supports the complexity of the dune filed there. Hermas et al. (2012) applied the matching technique of COSI-Corr, employing two SPOT-4 images spanning 292 days from 2007 to 2008, to monitor dune migration over part of the NSS. We compared our velocity results extracted from Sentinel-2 to those of Hermas et al. (2012) at the same locations, where both satellites have the same resolution (10 m). The mean velocity solution covered 1150 days from 2015-2019. Hermas et al. (2012) reported that the migration rates ranged from 4.0 to 20 m/ a, with mean of 6.8 m/a and 7.7 m/a for linear and barchan dunes, respectively. Meanwhile, the mean velocities from Sentinel-2 data recorded migration rates of 1.68 m/a and 1.22 m/a for linear and barchans dunes, respectively. These values were extracted after removing the pixels with values less ≤0.5 m/a, which were considered to represent interdune areas. We also compared the displacement values over the mountainous areas, where the study by Hermas et al. (2012) showed mean displacement of 5.5 m, compared to 0.027 ± 0.19 m/a for the mean velocity results of Sentinel-2. It is obvious that there were large differences in the extracted rates between both solutions. This variability can be attributed to the differences in the time of application and the complexity of the dune field. Additionally, there were large differences between the records of displacement over the mountainous area between our study and that of Hermas et al. (2012), wherein the latter introduced large values without any interpretation. Hermas et al.
(2012) compared their extracted results to previous records, which may be ineffective, especially in this study area where a wide range of migration is present over a complex dune field. However, the mean velocity results of Sentinel-2 showed lower uncertainties, and agreed well when compared to those values extracted from the mean velocity results of Landsat 8. This supports the use of our methods for representing the annual migration rates of dunes with high reliability.

Conclusions
We adapted a time-series inversion algorithm employing the matching results of both Landsat 8 and Sentinel-2. Our algorithm comprises two procedures: (1) the automatic selection of pairs with little differences in solar illumination to decrease the effects of shadows and seasonal variability, (2) the inversion of the time series. Owing to redundancy, the algorithm effectively decreases the uncertainties and enhances the spatial coverage. It is worth noting that our algorithm recommended limiting the sun angles as a first priority to control the cast shadows and the seasonal signals where other parameters used in this study may not be suitable for monitoring different targets, especially the cloud cover percentage. Therefore, the careful tuning of parameters would aid in balancing the number of images shared during inversion and the connectivity of the established network.
We applied our algorithm to monitor dune movement over a part of the NSS. The time-series displacements revealed the existence of seasonal signals in dune migration, which agreed with the existing wind regimes. The maximum annual migration rates recorded were 9.4 m/a and 15.9 m/a for Landsat 8 and Sentinel-2, respectively. We extracted the prevailing migration directions for both solutions, which also agreed with each other and with the RDD extracted from the wind data. These prevailing directions can be used as a proxy for wind direction, especially in vast deserts (e.g. Sahara desert and Arabian deserts), where no meteorological records are available. The proposed algorithm, which employed medium-resolution imagery, showed excellent performance in capturing the patterns of dune dynamics at the dune-field-scale in both the spatial and temporal domains. However, focusing on a single dune would require the use of high-resolution archives. Understanding the temporal patterns of dune movements and correlating them with continuous wind records can aid in predicting climatic changes and atmospheric conditions associated with dust transport. The method proposed here can be used to monitor the spatiotemporal patterns of glaciers and slow-moving landslides with good spatial coverage and smaller uncertainties, especially for archives that provide images with high revisit times.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  Fig. 5C; A-D are for the Landsat 8 solution, whereas E-H are for Sentinel-2. Each column represents different dune morphologies arranged, from left to right, as follows: barchan, linear, sand sheets, and the linear dunes in the white box in Fig. 5C. The magenta arrows represent the prevailing movement directions.