Open-source sea ice drift algorithm for Sentinel-1 SAR imagery using a combination of feature-tracking and pattern-matching

. An open-source sea ice drift algorithm for Sentinel-1 SAR imagery is introduced based on the combination of feature-tracking and pattern-matching. Feature-tracking produces an initial drift estimate and limits the search area for the consecutive pattern-matching, that provides small to medium scale drift adjustments and normalised cross-correlation values. The algorithm is designed to combine the two approaches in order to beneﬁt from the respective advantages. The considered 5 feature-tracking method allows a computationally efﬁcient computation of the drift ﬁeld and the resulting vectors show a high degree of independence in terms of position, lengths, direction and rotation. The considered pattern-matching method on the other side allows better control over vector positioning and resolution. The pre-processing of the Sentinel-1 data has been adjusted to retrieve a feature distribution that depends less on SAR backscatter peak values. Applying the algorithm 10 with the recommended parameter setting, sea ice drift retrieval with a vector spacing of 4 km on Sentinel-1 images covering 400 km x 400 km, takes about 4 minutes on a standard 2.7 GHz processor with 8 GB memory. The corresponding recommended patch size for the pattern-matching step, that deﬁnes the ﬁnal resolution of each drift vector is 34 × 34 pixels ( 2 . 7 × 2 . 7 km). To assess the potential performance after ﬁnding suitable search restrictions, calculated drift results from 246 Sentinel-15 1 image pairs have been compared to buoy GPS data, collected in 2015 between 15 th January and 22 nd April and covering an area from 80.5 ◦ N to 83.5 ◦ N and 12 ◦ E to 27 ◦ E. We found a logarithmic normal distribution of the displacement difference with a median at 352.9 m using HV


Introduction
Sea ice drift has a strong impact on sea ice distribution on different temporal and spatial scales.Motion of sea ice due to wind and ocean currents causes convergence and divergence zones, resulting in formation of ridges and opening/closing of leads.On large scales, ice export from the Arctic and the capability of optical sensors for reliable year-round sea ice monitoring.Unlike optical sensors, Space-borne Synthetic Aperture Radar (SAR) are active sensors, operate in the microwave spectrum and can produce high resolution images regardless of solar illumination and cloud cover.Since the early 1990's SAR sensors are delivering systematic acquisitions of sea ice covered oceans and Kwok et al. (1990) showed that sea ice displacement can be calculated from consecutive SAR scenes.The geophysical processor system from Kwok et al. (1990) has been used to calculate sea ice drift fields in particular over the Western Arctic (depending on SAR coverage) once per week with a spatial resolution of 10-25 km for the time period 1996-2012.This extensive dataset makes use of SAR data from RADARSAT-1 operated by the Canadian Space Agency, and from ENVISAT (Environmental Satellite) ASAR (Advanced Synthetic Aperture Radar) operated by ESA (European Space Agency).
To resolve drift details on a finer scale, a high-resolution sea ice drift algorithm for SAR images from ERS-1 (European Remote-sensing Satellite from ESA) based on pattern-matching was introduced by Thomas et al. (2008), that allowed drift calculation with up to 400 m resolution.Hollands and Dierking (2011) implemented their own modified version of this algorithm to derive sea ice drift from ENVISAT ASAR data.
To provide drift estimates also in areas where areal matching procedures (like cross and phase correlation) fail, Berg and Eriksson (2014) introduced a hybrid algorithm for sea ice drift retrieval from ENVISAT ASAR data using phase correlation and a feature based matching procedure that is activated if the phase correlation value is below a certain threshold.
The current generation of SAR satellites including RADARSAT-2 and Sentinel-1 are able to provide images with more than one polarisation.Komarov and Barber (2014) and Muckenhuber et al. (2016) have evaluated the sea ice drift retrieval performance with respect to the polarisation using a combination of phase/cross-correlation and feature-tracking based on corner detection respectively.Muckenhuber et al. (2016) has shown that feature-tracking provides on average around four times as many vectors using HV polarisation compared to HH polarisation.
After the successful start of the Sentinel-1 mission in early 2014, high-resolution SAR images are delivered for the first time in history within a few hours after acquisition as open-source data to all users.This introduced a new era in SAR Earth observation with great benefits for both scientists and other stack holders.Easy, free and fast access to satellite imagery facilitate the possibility to provide products on an operational basis.The Danish Technical University (Pedersen et al. (2015), http://www.seaice.dk/)provides an operational sea ice drift product based on Sentinel-1 data with 10 km resolution as part of the Copernicus Marine Environment Monitoring Service (CMEMS, http://marine.copernicus.eu).
The sea ice covered oceans in the European Arctic Sector represent an important area of interest for the Sentinel-1 mission and due to the short revisit time in the Arctic, our area of interest is monitored by Sentinel-1 on a daily basis (ESA, 2012).This paper follows up the work from Muckenhuber et al. (2016), who published an open-source feature-tracking algorithm to derive computationally efficient sea ice drift from Sentinel-1 data based on the open-source ORB algorithm from Rublee et al. (2011), that is included in the OpenCV Python package.We aim to improve the feature-tracking approach by combining it with pattern-matching.
Unlike Berg and Eriksson (2014), the feature-tracking step is performed initially and serves as a first guess to limit the search area of the pattern-matching step.
From a methodological point of view, algorithms for deriving displacement vectors between two consecutive SAR images are based either on feature-tracking or pattern-matching.
Feature-tracking detects distinct patterns (features) in both images and tries to connect similar features in a second step without the need for knowing the locations.This can be done computationally efficient and the resulting vectors are often independent of their neighbours in terms of position, lengths, direction and rotation, which can potentially be an important advantage for resolving shear zones, rotation and divergence/convergence zones.The considered feature-tracking approach identifies features without taking the position of other features into account and matches features from one image to the other without taking the drift and rotation information from surrounding vectors into account (Muckenhuber et al., 2016).However, due to the independent positioning of the features, very close features may share some pixels and since all vectors from the resolution pyramid are combined, the feature size varies among the matches, which implies a varying resolution.In addition, the resulting vector field is not evenly distributed in space and large gaps may occur between densely covered areas, which can eventually lead to missing a shear or divergence/convergence zone.
Pattern-matching, on the other hand, takes a small template from the first image at the starting location of the vector and tries to find a match on a larger template from the second image.Simple pattern-matching methods based on normalised cross-correlation often demand a considerable computational effort.Nevertheless, this approach is widely used, since it allows to define the vector positions.For practical reasons, a pyramid approach is generally used to derive high-resolution ice drift.This speeds up the processing, but potentially limits the independence of neighbouring vectors, since they depend on a lower resolution estimate (Thomas et al., 2008).
The objective of this paper is to combine the two approaches in order to benefit from the respective advantages.The main advantages of the considered feature-tracking approach are the computational efficiency and the independence of the vectors in terms of position, lengths, direction and rotation.
The considered pattern-matching method on the other side allows better control over vector positioning and resolution, which is a necessity for computing divergence, shear and total deformation.The presented algorithm, all necessary software requirements (python incl.Nansat, openCV and SciPy) and the satellite data from Sentinel-1 are open-source.A free and user friendly implementation shall support an easy distribution of the algorithm among scientists and other stakeholders.
The paper is organised as follows: The used satellite products and buoy data are introduced in Section 2. The algorithm description including data pre-processing is given in Section 3, together with tuning and performance assessment methods.Section 4 presents the pre-processing, parameter tuning and performance assessment results and provides a recommended parameter setting for the area and time period of interest.The discussion including outlook can be found in Section 5.

Data
The Sentinel-1 mission is a joint initiative of the European Commission and the European Space Agency (ESA) and represents the Radar Observatory for the Copernicus Programme, a European system for monitoring the Earth with respect to environmental and security issues.The mission includes two identical satellites, Sentinel-1A (launched in April 2014) and Sentinel-1B (launched in April 2016), each carrying a single C-band SAR with a centre frequency of 5.405 GHz and dual-polarisation support (HH+HV, VV+VH) also for wide swath mode.Both satellites fly in the same near-polar, sun-synchronous orbit and the revisit time is less than 1 day in the Arctic (ESA, 2012).The main acquisition mode of Sentinel-1 over sea ice covered areas is Extra Wide Swath Mode Ground Range Detected with Medium Resolution (EW GRDM) and the presented algorithm is built for processing this data type.The covered area per image is 400 km × 400 km and the data are provided with a pixel spacing of 40 m × 40 m in both HV and HH polarisation.The introduced algorithm can utilise both HV and HH channel.However, the focus of this paper is put on using HV polarisation (mainly acquired over the European Arctic and the Baltic sea), since this channel provides in our area of interest on average four times more feature tracking vectors than HH (Muckenhuber et al., 2016), representing a better initial drift estimate for the combined algorithm.
To illustrate the algorithm performance and explain the individual steps, we use an image pair acquired over Fram Strait.The acquisition times of the two consecutive images are 2015-03-28 07:44:33 (UTC) and 2015-03-29 16:34:52 (UTC), and the covered area is shown in Figure 1.This image pair covers a wide range of different ice conditions (multiyear ice, first-year ice, marginal ice zone etc.) and the ice situation is representative for our area and time period of interest.
To evaluate suitable search limitations and assess the potential algorithm performance, we use GPS data from drift buoys that have been set out in the ice covered waters north of Svalbard as part of the Norwegian Young Sea Ice Cruise (N-ICE2015) project of the Norwegian Polar Institute (Spreen and Itkin, 2015).The ice conditions during the N-ICE2015 expedition are describe on the project website (http://www.npolar.no/en/projects/n-ice2015.html) as challenging.The observed ice pack, mainly consisting of 1.3-1.5 m thick multiyear and first-year ice, drifted faster than expected and was very dynamic.Closer to the ice edge, break up of ice floes has been observed due to rapid ice drift and the research camp had to be evacuated and re-established four times.This represents a good study field, since these challenging conditions are expected in our area and time period of  interest.The considered GPS data have been collected in 2015 between 15 th January and 22 nd April, and cover an area ranging from 80.5 • N to 83.5 • N and 12 • E to 27 • E. The buoys recorded their positions either hourly or every three hours.In the later case, the positions have been interpolated for each hour.

Data pre-processing
To process Sentinel-1 images within Python (extraction of backscatter values and corresponding geolocations, reprojection, resolution reduction etc.), we use the Python toolbox Nansat (Korosov et al., 2016), that builds on the Geospatial Data Abstraction Library (http://www.gdal.org).As done in Muckenhuber et al. (2016), we change the projection of the the provided ground control points (latitude/longitude values given for certain pixel/line coordinates) to stereographic and use spline interpolation to calculate geographic coordinates.This provides a good geolocation accuracy also at high latitudes.The pixel spacing of the image is changed by averaging from 40 m to 80 m, which is closer to the sensor resolution of 93 m range × 87 m azimuth, and decreases the computational effort.
For each pixel p, the Sentinel-1 data file provides a digital number DN p and a normalisation coefficient A p , from which the normalised radar cross section σ 0 raw is derived by the following equation: The normalised radar cross section σ 0 raw reveals a logarithmic distribution and the structures in the sea ice are mainly represented in the low and medium backscatter values rather than in the highlights.
Therefore, we change the linear scaling of the raw backscatter values σ 0 raw to a logarithmic scaling and get the backscatter values σ 0 = 10 * lg(σ 0 raw ) [dB].A representative backscatter distribution over sea ice is shown in Figure 2. Using a logarithmic scaling provides a keypoint distribution for the feature tracking algorithm that depends less on high peak values, while the total number of vectors increases.To apply the feature-tracking algorithm from Muckenhuber et al. (2016), the SAR backscatter values σ 0 have to be converted into intensity values i with 0 ≤ i ≤ 255 for i ∈ R.This conversion is done by using Eq. ( 2) and setting all values outside the domain to 0 and 255.
The upper brightness boundary σ 0 max is set according to the recommended value from Muckenhuber et al. (2016), i.e. -18.86 dB and -10.97 dB for HV and HH respectively.The lower boundary σ 0 min was chosen to be -32.5 dB (HV) and -25.0 dB (HH), since this was found to be a reasonable range of expected backscatter values.Figure 3 shows the image pair Fram Strait after the conversion into intensity values.For the sake of computational efficiency, the same intensity value scaling is used for the pattern-matching step.

Sea ice drift algorithm
The presented sea ice drift algorithm is based on a combination of feature-tracking and patternmatching, and is designed to utilise the respective advantages of the two considered approaches.
Computationally efficient feature-tracking is used to derive a first estimate of the drift field.The concept of the FAST keypoint detector (Rosten and Drummond, 2006) to find corners on several resolution levels.The patch around each corner is then described using an modified version of the binary BRIEF descriptor from Calonder et al. (2010).To ensure rotation invariance, the orientation of the patch is calculated using the intensity-weighted centroid.Muckenhuber et al. (2016) applies a Brute Force matcher that compares each feature from the first image to all features in the second 225 image.The comparison of two features is done using the Hamming distance, that represents the number of positions in which the two compared binary feature vectors differ from each other.The best match is accepted if the ratio of the shortest and second shortest Hamming distances is below a certain threshold.Given a suitable threshold (and unique features), the ratio test will discard a high number of false matches, while eliminating only a few correct matches.Muckenhuber et al. (2016) found a suitable parameter setting for our area and time period of interest, including a Hamming distance threshold of 0.75, a maximum drift filter of 0.5 m/s, a patch size of 34 × 34 pixels and a resolution pyramid with 7 steps combined with a scaling factor of 1.2.
Due to the resolution pyramid, the considered feature area varies from 2.7 × 2.7 km to 9.8 × 9.8 km and the resulting drift field represents a resolution mixture between these boundaries.
We adjust the algorithm from Muckenhuber et al. ( 2016) by applying a logarithmic scaling for the SAR backscatter values σ 0 instead of the previous used linear scaling (Section 3.1).In addition, we extract for each vector the rotation information α, i.e. how much the feature rotates from the first to the second image.
Applying the adjusted feature-tracking algorithm provides a number of un-evenly distributed vectors (e.g.blue vectors in Figure 4) with start positions x 1f , y 1f on the first image (SAR 1 ), end positions x 2f , y 2f on the subsequent image (SAR 2 ) and corresponding rotation values α raw f .The index f represents a feature-tracking vector and ranges from 1 to F , with F being the total number of derived feature-tracking vectors.For the sake of computational efficiency, the vectors from all resolution pyramid levels are treated equally.
To avoid zero-crossing issues during the following filter and inter-/extrapolation process (in case the image rotation δ between SAR 1 and SAR 2 is close to 0 • ), a factor |180 − δ| is added to the raw rotation values α raw f using the following Equation: This centres the reasonable rotation values in the proximity of 180 • .After applying the filter and inter-/extrapolation process, the estimated rotation α is corrected by subtracting |180 − δ|.

II Filter
To reduce the impact of potentially erroneous feature-tracking vectors on the following steps, outliers are filtered according to drift and rotation estimates derived from least squares solutions using a third degree polynomial function.Considering a matrix A, that contains all end positions x 2f , y 2f in the following form If the simulated start position, derived from f (x 2f , y 2f , b), deviates from the feature-tracking start position x 1f , y 1f by more than 100 pixels, the vector is deleted.The same accounts for rotation outliers.If the simulated rotation deviates from the feature-tracking rotation α f by more than 60 • , the vector is deleted.We found a third degree polynomial function to be a good compromise between allowing for small to medium scale displacement and rotation discontinuities, and excluding very unlikely vectors, that eventually would disturb the following steps.The parameters for the filter process, i.e. 100 pixels (displacement) and 60 • (rotation), have been chosen according to visual interpretation using several representative image pairs.Figure 5 illustrates the filter process by depicting the results from image pair Fram Strait.

III First guess
The remaining feature-tracking vectors are used to estimate the drift incl.rotation on the entire first image, i.e. estimated x 2 , y 2 and α values are provided for each pixel on SAR 1 (Figure 6).The  3).Red points were identified as outliers and deleted.
constructed by p and the two opposite corners, e.g.A 1 is the area between p, and the corners with value v 2 and v 3 .
To provide a first guess for the surrounding area, values are estimated based on the least squares solutions using a linear combination of x 1 and y 1 .Considering a matrix C, that contains all start positions x 1f , y 1f in the following form As mentioned above, the rotation estimates α are now corrected for the adjustment applied in Equation 3, by subtracting |180 − δ|.

IV Pattern-matching
The estimated drift field derived from feature-tracking provides values for x 2 , y 2 and α at any location on SAR 1 .The representativeness of this estimate however, depends on the distance d to the closest feature-tracking vector.Therefore, small to medium scale adjustments of the estimates are necessary, depending on the distance d (NB: the representativeness also depends on the variability of the surrounding vectors, but for the sake of computational efficiency, we only consider the distance d as representativeness measure).We apply pattern-matching at chosen points of interest to adjust the drift and rotation estimate at these specific locations.
The used pattern-matching approach is based on the maximisation of the normalised crosscorrelation coefficient.Considering a small template t 1 around the point of interest from SAR 1 with size t 1s × t 1s and a larger template t 2 around the location x 2 , y 2 (defined by the corresponding first guess) from SAR 2 with size t 2s × t 2s , the normalised cross-correlation matrix NCC is defined as (Hollands , 2012): x ,y t 1 (x , y ) 2 x ,y t 2 (x + x , y + y )) 2 (9) given by dx and dy. ( To restrict the search area t 2s to a circle, we set all values of NCC that are further than t 2s /2 away from the centre position to zero.This limits the distance from the first guess to a constant value, rather than to an arbitrary value depending on the looking angle of the satellite.c To account for rotation adjustment, the matrix NCC is calculated several times: template t 1 is rotated around the initially estimated rotation α from α − β to α + β with step size ∆β.The angle β is the maximum additional rotation and represents therefore the rotation restriction.The NCC matrix with the highest cross-correlation value M CC is returned.
To illustrate the pattern-matching process, an example, taken from image pair Fram Strait, is shown in Figure 7.
The described process demands the specification of four parameters: t 1s , t 2s , β and ∆β.
The size of the small template t 1s × t 1s defines the considered area that is tracked from one image to the next and hence, affects the resolution of the resulting drift product.Sea ice drift might be different on different resolution scales.This is particularly an issue in the case of rotation.The feature-tracking vectors provide the first guess and this vector field should represent the same drift resolution as considered by the pattern-matching step.In order to be consistent with the resolution of the feature-tracking step and achieve our goal of a sea ice drift product with a spatial scaling of less than 5 km, we use the size of the feature-tracking patch of the pyramid level with the highest resolution to define the size of t 1 .That means, we use t s1 = 34 pixels (2.7 km).
The size of the larger template t 2s × t 2s restricts the search area on SAR 2 , i.e. how much the first guess can be adjusted geographically, and the angle β restricts the rotation adjustment of the  tion with the least time difference to the second satellite acquisition has been calculated using the following equation: where u and v represent eastward and northward drift components of the displacement vector derived by the algorithm, and U and V the corresponding drift components of the buoy.

Search restrictions evaluation
To find suitable values for restricting the size of the search window t 2s and the rotation range defined by β, we calculated drift vectors, that can be compared to the considered GPS buoy dataset, using restrictions that are computationally more demanding than we anticipate for the recommended setting, i.e.Based on an automatic search, we found 244 matching Sentinel-1 image pairs (consisting of 111 images), that allowed for comparison with 711 buoy vectors (buoy locations are shown in Figure 9).The distance D (Equation 15) between the buoy location at the time of the second image SAR 2 and the corresponding algorithm result, represents the error estimate for one vector pair.To identify algorithm results that are more likely erroneous, vector pairs with a value D above 1000 m are marked with red dots in Figure 10   The blue lines indicate the recommended setting for β (Equation 14) with dmin = 10 and dmax = 100.

Performance assessment
Using the recommended search restrictions from above, the algorithm has been compared to the N-ICE2015 GPS buoy data set (Figure 9) to assess the potential performance after finding suitable search restrictions for the area and time period of interest.The automatic search provided 246 image pairs (consisting of 111 images) and 746 vectors for comparison for the considered time period (15 th January to 22 nd April) and area (80.5 • N to 83.5 • N and 12 • E to 27 • E).NB: this is a higher number of vectors than found for the evaluation of the search restrictions, since the used search windows t 2 are smaller and vectors closer to the SAR edge may be included.
The results of the conducted performance assessment are shown in Figure 12.We found that the probability for a large D value (representative for the error) decreases with increasing maximum cross-correlation value M CC.Therefore we suggest to exclude matches with a M CC value below a certain threshold M CC min .This option is embedded into the algorithm, but can easily be adjusted or turned off by setting M CC min = 0. Based on the findings shown in Figure 12, we recommend a cross-correlation coefficient threshold M CC min = 0.4 for our time period and area of interest.
Using the suggested threshold reduces the number of vector pairs from 746 to 588 for the HV channel and to 478 for the HH channel.
The conducted performance assessment also reveals a logarithmic normal distribution of the distance D (Equation 15) that can be expressed by the following probability density function (solid red line in Figure 12): with µ and σ being the mean and standard deviation of the variable's natural logarithm.We found the mean and variance of the distribution lnN to be µ = 5.866 and σ 2 = 1.602 for HV polarisation and µ = 6.284 and σ 2 = 2.731 for HH polarisation (solid red lines in Figure 12).The medians of the logarithmic normal distribution are e µ = 352.9m for HV polarisation and e µ = 535.7 m for HH polarisation (dashed red lines in Figure 12).takes the same amount of computational effort for all image pairs consisting of Sentinel-1 images with 400x400 km coverage.

Recommended parameter setting
The process 'I Feature-tracking' depends on the setting of the feature-tracking algorithm and varies strongly with the chosen number of features.Using the recommended setting from Muckenhuber et al. (2016), that includes the number of features to be 100000, the presented computational effort can be considered representative for all image pairs, independent of chosen points of interest and overlap of the SAR scenes.
The last process 'II Pattern-matching and III Combination' however, depends on the considered image pair and the chosen drift resolution.The computational effort is proportional to the number of chosen points of interest.Given a evenly distributed grid of points of interest, the computational effort increases with overlapping area of the SAR scenes, since pattern-matching adjustments are only calculated in the overlapping area.The effort potentially decreases with a higher number of well distributed feature-tracking vectors, since the size of the search windows t 2 (and slightly the range of the angle β) increases with distance d to the closest feature-tracking vector.5 Discussion and outlook To estimate the potential performance of the introduced algorithm for given image pairs, given ice conditions, given region and given time, we compared drift results from 246 Sentinel-1 image pairs with corresponding GPS positions from the N-ICE2015 buoy data set.We found a logarithmic error distribution with a median at 352.9 m for HV and 535.7 m for HH (Figure 12).The derived error values represent a combination of the following error sources: -Timing: Buoy GPS data were collected every 1-3 hours and the timing does not necessarily match with the satellite acquisition time.
-Resolution: The algorithm returns the drift of a pattern (recommended size = 34 pixels, see Table 1), whereas the buoy measures the drift at a single location.
-Conditions: The ice conditions around the buoy is not known well enough to exclude the possibility that the buoy is floating in a lead.In this case, the buoy trajectory could represent a drift along the lead rather then the drift of the surrounding sea ice.
actual error of the algorithm.
A main advantage of the combined algorithm compared to simple feature-tracking, is the user defined positioning of the drift vectors.The current algorithm setup allows the user to choose whether the drift vectors should be positioned at certain points of interest or on a regular grid with adjustable spacing.Constricting the pattern-matching process to the area of interest minimises the computational effort according to the individual needs.
The recommended parameters shown in Table 1 are not meant as a fixed setting, but should rather give a suggestion and guideline to estimate the expected results and the corresponding computational effort.The parameters can easily be varied in the algorithm setup and should be chosen according to availability of computational power, needed resolution, area of interest and expected ice conditions (e.g.strong rotation).
The presented combination of feature-tracking and pattern-matching can be applied to any other application that aims to derive displacement vectors computationally efficient from two consecutive images.The only restriction is that images need to depict edges, that can be recognised as keypoints for the feature-tracking algorithm, and the conversion into intensity values i (Equation 2) needs to be adjusted according to the image type.
The remote sensing group at NERSC is currently developing a new pre-processing step to remove thermal noise on HV images over ocean and sea ice.First tests have shown a significant improvement of the sea ice drift results using this pre-processing step before applying the presented algorithm.This is ongoing work and will be included into a future version of the algorithm.
The European Space Agency is also in the process of improving their thermal noise removal for Sentinel-1 imagery.Noise removal in range direction is driven by a function that takes measured noise power into account.Until now, noise measurements are done at the start of each data acquisition, i.e. every 10-20 minutes, and a linear interpolation is performed to provide noise values every 3 seconds.The distribution of noise measurements showed a bimodal shape and it was recently discovered that lower values are related to noise over ocean while higher values are related to noise over land.This means, that Sentinel-1 is able to sense the difference of the earth surface brightness temperature similar to a passive radiometer.When the data acquisition includes a transition from ocean to land or vice versa, the linear interpolation fails to track the noise variation.The successors of Sentinel-1A/B are planned to include more frequent noise measurements.Until then, ESA wants to use the 8-10 echoes after the burst that are recorded while the transmitted pulse is still travelling and the instrument is measuring the noise.This will provide noise measurements every 0.9 seconds and allows to track the noise variations in more detail.In addition, ESA is planning to introduce a change in the data format during 2017 that shall remove the noise shaping in azimuth.These efforts are expected to improve the performance of the presented algorithm significantly (Personal Communication with Nuno Miranda, January 2017).
Having a computationally efficient algorithm with adjustable vector positioning allows not only to provide near-real time operational drift data, but also the investigation of sea ice drift over large Strait.The trackers send their position every 5-30 min to deliver drift information with high temporal resolution.This efforts shall help to gain a better understanding of short-term drift variability and by comparison with calculated sea ice drift, we will investigate how displacement vectors from subsequent satellite images relate to sea ice displacements with higher temporal resolution.
The focus of this paper in terms of polarisation was put on the HV channel, since this polarisation provides on average four times more feature-tracking vectors (using our feature-tracking approach) than HH and therefore delivers a finer initial drift for the first guess.We found our area of interest well covered with HV images, but other areas in the Arctic and Antarctic are currently only monitored in HH polarisation.Considering the four representative feature-tracking image pairs from Muckenhuber et al. (2016), the relatively best HH polarisation performance (i.e.most vectors from HH, while at the same time fewest vectors from HV) was provided by the image pair that had the least time difference, i.e. 8 h, compared to 31 h, 33 h and 48 h.Therefore, we assume that the HV polarisation provides more corner features that are better preserved over time.And more consistent features could potentially also favour the performance of the pattern-matching step, but this is only an assumption and has not been tested yet.Another argument is that the presented feature-tracking approach identifies and matches corners, which represent linear features.The linear features on HH images are more sensitive to changes in incidence angle, orbit and ice conditions than the linear features on HV images.This could explain the better feature-tracking performance of the HV channel.
However, pattern-matching is less affected by changing linear features and more sensitive to areal pattern changes.This could potentially mean that the HH channel performs better than HV when it comes to pattern-matching.However, at this point, these are just assumptions and will be addressed in more detail in our future work.
Utilising the advantage of dual polarisation (HH+HV) is certainly possible with the presented algorithm, but increases the computational effort.A simple approach is to combine the feature tracking vectors derived from HH and HV and produce a combined first-guess.Pattern-matching can be performed based on this combined first-guess for both HH and HV individually and the results can be compared and eventually merged into a single drift product.Having two drift estimates for the same position, from HH and HV pattern-matching respectively, would also allow to disregard vectors that disagree significantly.However, this option would increase the computational effort by two, meaning that the presented Fram Strait example would need about 8 min processing time.
After implementing the presented algorithm into a super-computing facility, we aim to test and compare the respective performance of HV, HH and HH+HV on large datasets to identify the respective advantages.
The current setting of the feature-tracking algorithm applies a maximum drift filter of 0.5 m/s.We found this to be a reasonable value for our time period and area of interest.However, when consider-ing extreme drift situations in Fram Strait and a short time interval between image acquisitions, this threshold should be adjusted.
As mentioned above, we deployed three GPS tracker in Fram Strait and they recorded their positions with a temporal resolution of 5-30 min between 8 th July until 9 th September 2016 in an area covering 75 • N to 80 • N and 4 • W to 14 • W. Considering the displacements with 30 min interval, we found velocities above 0.5 m/s on a few occasions, when the tidal motion adds to an exceptionally fast ice drift.
The GPS data from the hovercraft expedition FRAM2014-2015 In case the estimated drift from feature-tracking reaches velocities close to 0.5 m/s, the patternmatching step might add an additional degree of freedom of up to 8 km, which could eventually lead to a higher drift result than 0.5 m/s, depending on the time interval between the acquisitions.
The smaller the time difference, the larger is the potentially added velocity.In order to be consistent when combining the drift information from several image pairs with different timings, one should apply a maximum drift filter on the final drift product of the presented algorithm that has the same maximum velocity as the feature-tracking filter.The corresponding function is implemented in the distributed open-source algorithm.As an alternative, one could adjust the search window according to the time span.However, this would add additional complexity to both the algorithm and the parameter evaluation and needs more research on how the search window should be adjusted depending on the time span.For the sake of computational efficiency, we suggest the simple approach to remove final drift vectors above the maximum speed. 30km

Figure 1 .
Figure 1.Coverage of image pair Fram Strait that is used as representative image pair to explain the algorithm approach.The dashed rectangle depicts the area shown in Figure 4 to illustrate the vector distribution of the algorithm steps.

210Figure 4 .
Figure 4.The flowchart on the left depicts the five main steps of the algorithm.The right column illustrates the evolution of the drift results using image pair Fram Strait in HV polarisation and a grid with 4 km spacing.(NB: the part of the image pair that is depicted here is marked with a dashed rectangle in Figure 1.) Blue vectors are derived applying an adjusted version of the feature tracking algorithm from Muckenhuber et al. (2016).Black vectors indicate the initial drift estimate (first guess) based on filtered feature-tracking vectors.The final drift product (yellow to red vectors) are derived from combining the first guess with pattern-matching adjustment and applying a minimum cross-correlation value.A total of 4725 vectors have been found on image pair Fram Strait with a M CC value above 0.4 in 4 min.

Figure 5 .
Figure 5. Filter process applied on image pair Fram Strait in HV polarisation.The x-axis represent the simulated start position and rotation, derived from f (x 2f , y 2f , b) and the y-axis represent the feature-tracking start position x 1f , y 1f and rotation α f .NB: the image rotation is δ = 129.08 • , which means the rotation was adjusted by 50.92 • (Equation3).Red points were identified as outliers and deleted.

Figure 6 .
Figure 6.Example of estimated drift and rotation (first guess) based on filtered feature-tracking vectors using image pair Fram Strait in HV polarisation.The three panels show the components x2, y2 of the estimated end positions and the estimated rotation α for each pixel on on the coordinate system x1, y1 of the first image (SAR1).

Figure 7 .
Figure 7. Pattern-matching using initial drift estimate from feature-tracking: The small template t1 (left) around the point of interest on SAR1 is rotated from α−β to α+β and matched with the large template t2 (middle) from SAR2, that has its centre at the estimated end position x2, y2.The right contour plot shows the normalised crosscorrelation matrix NCC of the rotation β * that provided the highest maximum cross-correlation coefficient M CC.The estimated end position x2, y2 of this example has to be adjusted by dx = −21 pixels, dy = 32 pixels to fit with the location of M CC = 0.71.Rotation adjustment β * was found got be 3 • .NB: X and Y -axis represent pixel coordinates.

Figure 8 .
Figure 8. Example to illustrate the distribution of distance d to the closest feature-tracking vector using image pair Fram Strait in HV polarisation.Values outside the range dmin ≤ d ≤ dmax are set to dmin = 10 and dmax = 100.The points with value dmin represent the start positions x 1f , y 1f of the feature-tracking vectors on the coordinate system x1, y1 of SAR1.The figure depicts the matrix that the algorithm considers for the distribution of d.

Figure 9 .
Figure 9. Considered buoy locations from the N-ICE2015 expedition that were used for comparison with algorithm results.Green and blue colour indicates start locations (on SAR1) to which the algorithm provided vectors with a M CC value above and below 0.4 using (a) HV and (b) HH polarisation.
and Figure 11.Vector pairs with D < 1000 m are plotted with 410 black dots.

Figure 10 andFigure 10 .Figure 11 .
Figure 10 and Figure 11 show the resulting pattern-matching adjustment of location (dx, dy) and rotation (dβ) using the computationally demanding restrictions.The values are plotted against distance d to the next feature tracking vector in order to identify the dependence of the parameters on d.The blue lines in Figure 10 and Figure 11 indicate the recommended restrictions.This represents a 415

Figure 12 .
Figure 12.Calculated ice drift using recommended search restrictions compared to buoy GPS data using (a,b,c) HV and (x,y,z) HH polarisation.Light grey represents vectors with maximum cross-correlation values M CC < 0.4 and results after using the suggested threshold M CCmin = 0.4 are shown in black.(a,x) M CC values against distance D (Equation 15) between algorithm and buoy end position.The blue line indicates the recommended setting for M CCmin = 0.4.(b,y) Logarithmic histogram of distance D with 100 bins between 10 m and 10 5 m including two logarithmic normal distributions that were fitted to all results (grey) and to the filtered results with M CC > 0.4 (solid red line).(c,z) Comparison of drift distance derived from algorithm against buoy displacement for the filtered results with M CC > 0.4.
areas and long time periods.Our next step is to embed the algorithm into a super-computing facility to further test the performance in different regions, time periods and ice conditions and evaluate and combine the results of different polarisation modes.The goal is to deliver large ice drift datasets and open-source operational sea ice drift products with a spatial resolution of less than 5 km.This work is linked to the question how to combine the different timings of the individual image pairs in a most useful way.Having more frequent satellite acquisitions, as we get with the Sentinel-1 satellite constellation, enables to derive displacements for shorter time gaps and the calculated vectors will reveal more details e.g.rotational motion due to tides.As part of a scientific cruise with KV-Svalbard in July 2016, we deployed three GPS trackers on loose ice floes and pack-ice in Fram (https://sabvabaa.nersc.no),that was collected with a temporal resolution of 10 s between 31 st August 2014 until 6 th July 2015, did not reveal a single 30 min interval during which the hovercraft was moved by ice drift more than 0.45 m/s.The hovercraft expedition started at 280 km south from the North Pole towards the Siberian coast, crossed the Arctic Ocean towards Greenland and was picked up in the north-western part of Fram Strait.
(x , y ) and t 2 (x , y ) representing the value of t 1 and t 2 at location x ,y .The summations are done over the size of the smaller template, i.e. x , y , x and y go from 1 to t 1s .Template t 1 is moved with step size 1 pixel over template t 2 both in horizontal (x) and vertical (y) direction and the cross-correlation values for each step are stored in the matrix NCC with size (1 + t s2 − t s1 ) × (1 + t s2 − t s1 ).The highest value in the matrix NCC, i.e. the the maximum normalised cross-correlation value M CC, represents the location of the best match and the corresponding location adjustment is

Table 1 .
Recommended parameter setting for sea ice drift retrieval from Sentinel-1 using the presented algo-The processing time depends on the parameter setting and the chosen vector distribution.Using the recommended parameter setting from Table1, allows high-resolution sea ice drift retrieval from a Sentinel-1 image pair within a few minutes.Figure4depicts calculated ice drift vectors for the image pair Fram Strait on a grid with 4 km (50 pixels) spacing.The corresponding processing times are shown in Table2.The calculations have been done using a MacBook Pro from early 2013 with a 2.7 GHz Intel Core i7 processor and 8 GB 1600 MHz DDR3 memory.The total processing time for 4725 vectors with a normalised cross-correlation value above 0.4, is about 4 minutes.This can be considered a representative value for an image pair with large overlap, good coverage with featuretracking vectors and 4 km grid spacing.The initial process in Table2'Create Nansat objects from Sentinel-1 image pair and read matrixes'

Table 2 .
Processing time for sea ice drift retrieval from image pair Fram Strait on a grid with 4 km (50 pixels) spacing using HV polarisation (Figure4).Representative for an image with large overlap and good coverage with feature-tracking vectors.