An improvement in accuracy and spatial resolution of the pattern-matching sea ice drift from SAR imagery

ABSTRACT

Sea ice drift is a crucial parameter for sea ice flux, atmospheric and ocean circulation, and ship navigation. Pattern matching is widely used to retrieve sea ice drift from Synthetic Aperture Radar (SAR) data, but it often yields mismatched vectors and coarse spatial resolution. This study presents a framework to enhance the spatial resolution and accuracy of pattern-matching sea ice drift derived from SAR images. The framework employs the Accelerated-KAZE feature extraction method and Brute-Force feature matcher to extract feature-tracking sea ice drift vectors from SAR data, with mismatched vectors subsequently removed. The pattern-matching vectors are then refined by fusing with these feature-tracking vectors, using a Co-Kriging algorithm. Using the sea ice drift product from the Technical University of Denmark space as the pattern-matching vector field for refinement, the framework's effectiveness is evaluated by comparing the refined vectors with buoy displacements and pattern-matching vectors across five selected regions. Results show a reduction in velocity and direction root mean square error (RMSE) by 0.47 km/d (22%) and 4.97° (28%), respectively, and an enhanced spatial resolution from 10 km to 1 km. The findings demonstrate the framework's success in improving the accuracy and resolution of pattern-matching sea ice drift from SAR imagery.

KEYWORDS:

• Sea ice drift
• pattern matching
• SAR
• feature tracking
• fusion

1. Introduction

In the context of global warming, sea ice has undergone significant changes in recent years (IPCC 2021). Sea ice serves as a vital mediator in the heat exchange process between the ocean and atmosphere (Jacob et al. 2012; Li et al. 2021; Demchev et al. 2017), with rapid changes having a significant influence on regional and global climate systems (Wang and Zhao 2020; Makshtas, Shoutilin, and Andreas 2003). Sea ice drift, a key dynamic characteristic of sea ice, is vital for monitoring rapid changes in sea ice conditions (Martin and Augstein 2000; Zhang et al. 2003; Volkov, Mushta, and Demchev 2019; Petrou and Tian 2017). Estimations of sea ice drift are important at various scales. At a regional level, such as the pan-Arctic, information on ice drift is essential for measuring the transfer of mass, momentum, energy, and heat at the interface between the atmosphere, sea ice, and ocean (Kramer, Johnsen, and Brekke 2015; Lei et al. 2020; Spreen et al. 2006; Zhao and Liu 2007). Furthermore, sea ice drift data is necessary for understanding sea ice responses to wind stress and ocean currents (Fan et al. 2020; Oikkonen and Haapala 2011; Stern and Lindsay 2009; Spreen, Kwok, and Menemenlis 2011). At a local scale, sea ice drift significantly impacts ship navigation, gas and oil exploration, and fishing activities in the Arctic Ocean (Yu et al. 2021; Ducklow et al. 2007; Petrou and Tian 2017; Askne and Dierking 2008; Girard-Ardhuin and Ezraty 2012).

Sea ice drift measurements began with the utilization of ship observations, drift stations, and buoys (Girard-Ardhuin and Ezraty 2012). However, these measurements offered limited data with sparse spatial distribution (Colony and Thorndike 1984; Girard-Ardhuin and Ezraty 2012; Meier and Dai 2006). Remote sensing enables extensive spatial and temporal observations, allowing for comprehensive and regular assessments of sea ice drift on a large scale (Berg and Eriksson 2014; Hollands 2012; Komarov and Barber 2014; Lavergne et al. 2010; Shokr, Wang, and Liu 2020). Sea ice drift retrieval employs a variety of remotely sensed data, such as passive microwave data (Spreen, Kwok, and Menemenlis 2011; Meier and Dai 2006; Lavergne et al. 2010), microwave scatterometer data (Haarpaintner 2006; Girard-Ardhuin and Ezraty 2012), optical imagery(Fang et al. 2023; Ninnis, Emery, and Collins 1986; Petrou and Tian 2017) and synthetic aperture radar (SAR) data (Berg and Eriksson 2014; Demchev et al. 2017; Kwok et al. 1998; Komarov and Barber 2014; Muckenhuber and Sandven 2017; Howell, Brady, and Komarov 2022; Korosov and Rampal 2017; Li et al. 2022a; 2022b; Qiu and Li 2022). However, the coarse spatial resolution of microwave scatterometers and passive microwave radiometers, typically between 6.25 and 25 km, limits their use in high-resolution sea ice drift retrieval (Petrou and Tian 2017). Similarly, cloud cover can impair the coverage of optical sensor-derived sea ice drift vectors (Dybkjaer 2013). In contrast, SAR imagery, with its weather – and sunlight-independent monitoring capabilities and high spatial resolution (typically from 1 to 100 m), is extensively utilized for sea ice motion estimation (Petrou and Tian 2017).

Pattern matching methods are widely employed for sea ice drift retrieval using SAR data (Kwok et al. 1998; Komarov and Barber 2014; Howell, Brady, and Komarov 2022; Karvonen 2012). The well-known sea ice drift product from SAR imagery, released by the Technical University of Denmark (DTU) space, is generated based on the pattern matching method (Saldo and Hackett 2022). Generally, sea ice drift is tracked from image pairs with overlapping areas, acquired at different times. In pattern matching methods, both images are divided into subregions or templates. Template matching is then used to identify the best matches between image pairs, with the relative displacement and direction of the best-matched templates denoting sea ice drifts (Girard-Ardhuin and Ezraty 2012; Haarpaintner 2006; Spreen, Kwok, and Menemenlis 2011; Tschudi, Meier, and Stewart 2020; Ninnis, Emery, and Collins 1986; Lavergne et al. 2021; Hyun and Kim 2017). Pattern matching is simple to implement (Lavergne et al. 2010). However, as a transversal algorithm, it is susceptible to image noise, resulting in some mismatched vectors and thereby affecting the accuracy of the obtained sea ice drift. Moreover, pattern matching methods often yield a non-dense vector field with a coarser spatial resolution than the input data for efficiency purposes (Petrou and Tian 2017).

In recent years, feature tracking methods have been employed for sea ice drift retrieval from high-resolution SAR data (Muckenhuber, Korosov, and Sandven 2016; Demchev et al. 2017; Howell, Brady, and Komarov 2022). These methods involve extracting image feature points from sequential images, matching features of image pairs, and connecting them to generate sea ice displacement vectors. Some feature detection and description algorithms including Speeded-Up Robust Features (SURF), Accelerated-KAZE (A-KAZE), Scale Invariant Feature Transform (SIFT), and Oriented FAST and Rotated BRIEF (ORB) have been employed in sea ice drift retrieval (Muckenhuber, Korosov, and Sandven 2016; Demchev et al. 2017; Wang et al. 2020). Since feature tracking methods only match features, the resulting vectors may not be evenly distributed spatially, with gaps potentially occurring over sea ice areas lacking good textures (Muckenhuber and Sandven 2017). Consequently, this approach has not been adopted for generating sea ice drift products. However, the extracted features, being at pixel and sub-pixel levels, yield densely distributed vectors in areas with good textures. Furthermore, the derived drift vectors demonstrate independence regarding position, length, direction, and rotation, resulting in good performance along shear, rotation, and divergence/convergence zones (Muckenhuber, Korosov, and Sandven 2016; Muckenhuber and Sandven 2017). Generally, feature tracking methods demonstrate superior accuracy compared to pattern matching methods in extracting sea ice drift vectors (Muckenhuber, Korosov, and Sandven 2016).

Given the complementary advantages of pattern-matching and feature-tracking vector fields, recent studies have combined these two methods to estimate SAR-based sea ice motion (Muckenhuber and Sandven 2017; Korosov and Rampal 2017; Li et al. 2022a; 2022b; Qiu and Li 2022). The combined algorithm consists of three main steps: first, the feature-tracking algorithm is used to extract initial sea ice displacements; second, interpolation is performed to estimate the vectors on a regular grid; and third, pattern matching is used to refine each interpolated vector. In most studies, the ORB algorithm is used for feature tracking, while the maximum cross-correlation (MCC) is used for pattern matching (Muckenhuber and Sandven 2017; Korosov and Rampal 2017; Li et al. 2022a; 2022b; Qiu and Li 2022). Li et al. (2022a) used a quadtree ORB (Q-ORB) algorithm instead of the ORB method to obtain effective sea ice motion vectors in weak texture areas in SAR imagery. However, the initial feature-tracking vectors are unevenly distributed. Interpolation using spatially non-uniformly distributed fields generally yields high accuracy when interpolating in densely distributed areas, but may result in low accuracy in sparsely distributed areas. The subsequent MCC algorithm can only make small to medium adjustments to the interpolated vectors (Muckenhuber and Sandven 2017). Therefore, if the accuracy of the interpolated vectors is low, it may affect the accuracy of the adjusted vectors. On the other hand, pattern-matching sea ice drift vectors are evenly distributed. Thus, it is possible to consider using pattern-matching sea ice drift vectors as the background field and fusing feature-tracking vectors for adjustment.

In conclusion, pattern-matching sea ice drift vectors retrieved from SAR images are evenly distributed, yet their accuracy is relatively low, accompanied by a lower spatial resolution compared to the SAR input data. On the other hand, feature-tracking vectors from the same SAR data provide more densely distributed and accurate sea ice drift measurements, despite their uneven distribution. By fusing feature-tracking vectors with pattern-matching vectors, it is possible to improve the accuracy and spatial resolution of the latter. Thus, the objective of this study is to present a framework designed to enhance the accuracy and spatial resolution of pattern-matching sea ice drift vectors from SAR imagery, achieved by fusing them with feature-tracking vectors derived from the same source dataset.

2. Study area and data

2.1. Study area

This study uses sea ice drift product from the DTU space (referred to as DTU product hereafter) as the pattern-matching vectors, which has been generated using Sentinel-1 SAR data since December 15th, 2014. Considering the spatial heterogeneity of Arctic sea ice drifts and the availability of DTU sea ice drift data, Sentinel-1 SAR images, and buoys, five regions are selected for experimental case studies (Figure 1). Due to the Sentinel-1B anomaly starting from December 23, 2021, data from 2020 are used in this study to cover a wide study area. The Arctic annual average sea ice velocity for 2020 is analyzed and shown in Figure 1, calculated using the sea ice drift product with 25 km spatial resolution released by the National Snow and Ice Data Center (Tschudi, Meier, and Stewart 2020).

Figure 1. (a) Spatial distribution of selected DTU sea ice drift scenes; (b) Geographic locations of five experimental case regions and start positions of buoys used for validation.

Within the five regions, Region I encompasses the Chukchi Sea and Beaufort Sea (CBS), where sea ice rotates clockwise due to the influence of the Beaufort Gyre(Kwok, Spreen, and Pang 2013). The average ice velocity in CBS is 4–10 km/d in 2020, with the Chukchi sea ice moving faster than the Beaufort sea ice. Regions II, III, and IV are all situated in the central Arctic Ocean, exhibiting varying average sea ice velocities in 2020, with the highest to lowest order being Regions IV > II > III. Region II lies north of the Kara Sea and features an average sea ice velocity of 4–8 km/d in 2020. Region III, located north of the Canadian Arctic Archipelago, is influenced by both the Beaufort Gyre and Transpolar Drift, showing a low average velocity of 2–4 km/d in 2020. Region IV, positioned northeast of Greenland and including the Fram Strait, serves as the primary passage for Arctic sea ice export and the endpoint of the Transpolar Drift. The average sea ice velocity in this region (8-14 km/d in 2020) is faster than in other regions. Region V is situated in the East Siberian Sea, which is the starting point of the Transpolar Drift that transports water and ice across the Pole and along the eastern coast of Greenland. Dominated by first-year sea ice (Zhang et al. 2022), this region exhibits slow movement (averaging 2–4 km/d in 2020) near the shore and faster movement (averaging 4–10 km/d in 2020) off the shore.

2.2. Data and pre-processing

The datasets utilized in this study include 366 DTU sea ice drift product scenes, 277 Sentinel-1 A/B SAR images, 359 buoys, and a physical vector dataset of land. Sentinel-1 A/B SAR data are employed for sea ice drift retrieval based on the feature-tracking method. GPS information from buoys serves to assess the accuracy of the retrieved sea ice drift vectors, while the physical vector data of land are employed for pre-processing Sentinel-1 A/B SAR data. The geographic distribution of the selected 366 DTU sea ice drift scenes can be seen in Figure 1, with data used in each region detailed in Table 1. The start dates of image pairs correspond to the dates of the first SAR image acquisition.

Table 1. Characteristics of data used in five regions.

Region	Number of DTU NRT sea ice drift data	Number of Sentinel-1 Images	Start Date of image pairs (in 2020)	Number of Buoys
I	87	62	12.05; 12.08; 12.16; 12.17; 12.27; 12.28.	41
II	34	17	01.01; 01.02; 01.03; 01.04; 01.05; 01.06.	138
III	68	40	11.12; 11.13; 11.14; 11.20; 11.21.	19
IV	152	70	05.06; 05.07; 05.08; 05.09; 05.13; 05.14; 05.15.	149
V	25	38	03.04; 03.07; 03.08; 03.09; 03.12; 03.20.	12
All	366	277	/	359

The DTU sea ice drift product, widely utilized in numerous studies (Qiu and Li 2022; Hwang 2013; Muckenhuber, Korosov, and Sandven 2016), has a time series commencing on March 9th, 2010. Initially, Envisat ASAR data were used until 2012, when they were replaced by RADARSAT-2 SAR data. Since December 15th, 2014, the product has been generated using Sentinel-1 SAR data (Saldo and Hackett 2022). The DTU product comprises two datasets: one set consists of near-real-time (NRT) SAR swath data-derived sea ice drift, with files named according to the start and end times of the SAR image pairs used for drift retrieval; the other set is a 24-hour mean composite, calculated from the NRT sea ice drift dataset. The latter dataset features a daily temporal resolution and a 12-hour update frequency, providing two types of sea ice drift data with start times of 00:00 and 12:00 every day. The spacing of the DTU sea ice drift product grid is approximately 10 km, and it is set in a Polar Stereographic projection (Saldo and Hackett 2022). In this study, the DTU NRT dataset is used as the pattern-matching sea ice drift. The quality status of data file in DTU product is classified into four categories: ‘nominal_quality’, ‘low_quality’, ‘rejected_by_filter’, and ‘no_input_data’, corresponding to quality parameters 0, 1, 2, and 3, respectively (Saldo and Hackett 2022). Only data files showing a quality parameter of 0 (i.e. nominal_quality) are used. Access to the DTU sea ice drift product is available through the Copernicus Marine Environment Monitoring Service (CMEMS) website (https://marine.copernicus.eu/).

Sentinel-1 SAR data are selected based on the timestamps of the chosen NRT DTU sea ice drift data in five regions. The chosen Sentinel-1A and 1B SAR data comprise Extra Wide (EW) mode (swath width 410 km) images at HH and HV polarization from Level-1 Ground Range Detected medium resolution (GRDM, 40 m spatial resolution) products. Both HH and HV channels are used in this study. Sentinel-1 data can be easily and freely accessed from the Alaska Satellite Facility Distributed Active Archive Centers (ASF DAAC) website (https://asf.alaska.edu/). For high-accuracy data geocoding in sea ice drift retrieval, the Precise Orbit Ephemerides (POE) data of Sentinel-1 products are collected for preprocessing from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). The preprocessing of Sentinel-1 SAR consists of seven successive steps: precise orbit calibration, thermal noise removal, radiometric correction, multilooking, dB conversion, ellipsoid correction, and land masking. Precise orbit calibration and ellipsoid correction utilize POE orbit files and ellipsoid coefficients for accurate SAR data geocoding. The removal of thermal noise from channel intensities is achieved by utilizing the module provided by the European Space Agency (ESA) within its Sentinel Application Platform (SNAP) software. Radiometric correction involves converting digital pixel values to radiometrically calibrated SAR backscatter (Filipponi 2019). Subsequently, 3 by 3 looks are conducted along range and azimuth to reduce speckle effects, which upscales SAR imagery to 120 m. DB conversion is to convert unitless backscatter to dB using logarithmic transformation (Filipponi 2019). Land masking is to mask land from feature tracking for computational efficiency. Land information is obtained from the physical vector dataset of land with a 1:50 million scale, accessed from the Nature Earth website (https://www.naturalearthdata.com/).

Data from buoys on the sea ice are collected from the Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC) project (https://www.meereisportal.de/en/mosaic/) and the International Arctic Buoy Program (IABP, https://iabp.apl.uw.edu/). Various buoy types, such as Surface Velocity Profiler (SVP) and Ice-Tethered Profiler (ITP), are included. Only the geographic location (latitude and longitude) and timestamp from each buoy are used for validation. Specifically, IABP buoys provide locations at 0:00 and 12:00 (UTC) daily, while MOSAiC buoys offer locations on an hourly basis. The start locations of buoys used for validation are shown in Figure 1.

3. Method

In this study, a framework is proposed to enhance the accuracy and spatial resolution of pattern-matching sea ice drift derived from SAR imagery. This framework involves two key steps. Initially, feature-tracking sea ice drift vectors are extracted from pre-processed SAR images, utilizing a suite of algorithms: A-KAZE-based feature extraction (Alcantarilla, Nuevo, and Bartoli 2013), Brute-Force (BF) matcher-based feature matching (Lowe 2004), and bad matches filtering. Following this, the pattern-matching sea ice drift is refined by fusion with the feature-tracking vectors using the Co-Kriging algorithm (Goovaerts 1997). The flowchart of the proposed framework is presented in Figure 2.

Figure 2. Proposed framework flowchart for refining the pattern-matching sea ice drift from SAR imagery.

3.1. Retrieval of the feature-tracking sea ice drift

The retrieval of feature-tracking sea ice drift from preprocessed Sentinel-1 SAR image pairs consists of four successive steps. First, feature points on HH or HV image pairs are detected and described using the A-KAZE method (Alcantarilla, Nuevo, and Bartoli 2013). Secondly, feature points in SAR image pairs are aligned by utilizing the BF matcher based on their corresponding feature descriptors (Lowe 2004). Thirdly, bad matches are filtered out, and the remaining matched feature points of image pairs are connected as sea ice drift vectors. Finally, the filtered results from the distinct polarization image pairs are combined to produce the final feature-tracking result.

The A-KAZE algorithm (Alcantarilla, Nuevo, and Bartoli 2013) is selected over other feature detection methods, such as SIFT (Lowe 2004), SURF (Bay, Tuytelaars, and Van Gool 2006), and ORB (Rublee et al. 2011), for feature detection and description as it outperforms them in several aspects, including the speed of feature detection and description, as well as the number of extracted feature points (Demchev et al. 2017). The A-KAZE method constructs a nonlinear scale space to identify feature points and subsequently generates feature descriptors using the Modified-Local Difference Binary (M-LDB) descriptor (Alcantarilla, Nuevo, and Bartoli 2013). The M-LDB descriptor offers high efficiency, scale and rotation invariance, and requires low storage space (Alcantarilla, Nuevo, and Bartoli 2013).

The BF matcher (Lowe 2004) is utilized for feature matching because it surpasses other matchers, such as Fast Library for Approximate Nearest Neighbors (FLANN) (Muja and Lowe 2009), in terms of matching a larger number of features with higher accuracy (Noble 2017). BF matcher (Lowe 2004) compares feature descriptors between two images and matches the closest features by utilizing Hamming distance (Hamming 1986). Subsequently, the Nearest Neighbor Distance Ratio (NNDR) test (Lowe 2004) is applied to the initial matches from the BF matcher. The NNDR test considers a match incorrect if the ratio between the distance of the second-best match and the best match is greater than $r$ (Lowe 2004). In this study, $r$ is set to 0.8, as this threshold can remove 90% of the incorrect matches and retain less than 10% of the correct ones (Lowe 2004). The remaining matches from the NNDR test are connected as initial sea ice drift distance vectors. These distance vectors are divided by the time interval of image pairs to obtain sea ice velocity vectors (unit: km/d).

The assumption for bad match filtering in this study is that there should be coherence in sea ice drift within a small area (Lavergne 2016). By considering the correlation with neighboring vectors, each initial vector is examined. The neighborhood scope in this study is set to 25 km. First, a vector is discarded if no vector is present within a 25 km radius (e.g. the blue vector in Figure 3). Next, the velocity and direction differences are calculated between an individual vector and the simulated vector, using neighboring vectors within a 25 km circle based on the Inverse Distance Weighted (IDW) model. The vector being tested is validated (e.g. the yellow vector in Figure 3) if both velocity and direction differences are below three times the standard deviation of its neighboring vectors (Li et al. 2022). Otherwise, the testing vector is deleted (e.g. the red vector in Figure 3).

Figure 3. Example cases of bad match filtering. Red, yellow, blue, and black vectors represent retrieved vectors after the NNDR test. Green vectors are simulated by neighboring vectors using IDW.

Previous studies suggest that utilizing ice motion vectors from both HH and HV channels together leads to a more comprehensive result as opposed to using them separately to derive ice motion information (Komarov and Barber 2014; Muckenhuber and Sandven 2017). Therefore, in this study, the feature-tracking sea ice drift extraction is implemented on both HH and HV polarization image pairs, and the resulting sea ice drifts from these distinct polarization image pairs are combined to generate the final feature-tracking sea ice drift result. To combine the results from different polarization image pairs, a union step is employed. In situations where two vectors are present on the same feature point, the vector that exhibits a smaller difference from its neighboring vectors within a 25 km radius is identified as the final result. For feature points that only have a single vector, the original values are retained.

3.2. Refinement of the pattern-matching sea ice drift

In this study, the challenge of improving the accuracy and spatial resolution of pattern-matching sea ice drift is approached as an issue of estimating fine-resolution sea ice drift, where pattern-matching vectors serve as background values and the feature-tracking sea ice drift vectors extracted are viewed as observations. The Kriging method (Goovaerts 1997) is employed, as it is an effective estimation method utilized in various studies, such as estimating rainfall data (Adhikary, Muttil, and Yilmaz 2017) and population density (Wu and Murray 2005). Sea ice drift can be estimated using the pattern-matching sea ice drift vectors with the ordinary Kriging method, based on the spatial autocorrelation information among sea ice drifts. However, the coarse spatial resolution of the pattern-matching sea ice drift results in the loss of finer details at higher resolutions. The Co-Kriging method is employed to improve the estimation of fine-resolution sea ice drift by utilizing the spatial autocorrelation and correlation information between pattern matching and feature-tracking vectors. The estimated sea ice drift from Co-Kriging is derived using the following equation: (1) ${\hat{D}}_{CK}(L_0)=\sum _{i_1=1}^n{w}_{i_1}F(L_{i_1})+\sum _{i_2=1}^m{w}_{i_2}P(L_{i_2})$(1) with $\sum _{i=1}^n{w}_{i_1}=1$and $\sum _{j=1}^m{w}_{i_2}=0$

where ${\hat{D}}_{CK}(L_0)$is the estimated value of sea ice drift at the $L_0$ location. $w_{i_1}$ and $w_{i_2}$ are the Co-Kriging weights assigned to feature-tracking drift $F$ at the location $L_{i_1}$ and the pattern-matching drift $P$ at the location $L_{i_2}$, respectively. The total number of feature-tracking and pattern-matching vectors within the surrounding region (i.e. 25 km-radius) of the location $L_0$ are represented by $n$ and $m$.

The Co-Kriging weights $w_{i_1}$ and $w_{i_2}$ can be solved using the spatial autocorrelation information and cross-dependency between feature-tracking and pattern-matching sea ice drifts. A mathematical function, i.e. a variogram model, is utilized to indicate the level of spatial autocorrelation present in the input data. Cross-variogram models can quantify the degree of cross-dependency between two datasets (Goovaerts 1997). Experimental variograms ${\hat{\gamma }}_{FF}$, ${\hat{\gamma }}_{PP}$, and experimental cross-variogram ${\hat{\gamma }}_{FP}(d)$ (equal to ${\hat{\gamma }}_{PF}(d)$) are derived from equations (2), (3), and (4), respectively. (2) ${\hat{\gamma }}_{FF}(d)={\displaystyle \frac{1}{2N(d)}\sum _{k=1}^{N(d)}[F(L_k+d)-F(L_k){]}^2}$(2) (3) ${\hat{\gamma }}_{PP}(d)={\displaystyle \frac{1}{2N(d)}\sum _{k=1}^{N(d)}[P(L_k+d)-P(L_k){]}^2}$(3) (4) ${\hat{\gamma }}_{FP}(d)={\hat{\gamma }}_{PF}(d)={\displaystyle \frac{1}{2N(d)}\sum _{k=1}^{N(d)}[F(L_k+d)-F(L_k)]\times [P(L_k+d)-P(L_k)]}$(4) where ${\hat{\gamma }}_{FF}(d)$ and ${\hat{\gamma }}_{PP}(d)$ are the experimental variogram values for feature-tracking drifts and pattern-matching drifts, respectively, at a distance of $d$. ${\hat{\gamma }}_{FP}(d)$ and ${\hat{\gamma }}_{PF}(d)$ are the values of the experimental cross-variogram between sea ice drifts from feature tracking and pattern matching with a distance of $d$. The count of data pairs located at a distance apart, $d$, is indicated by $N(d)$.

Once the experimental variograms and cross-variograms are obtained, theoretical models are employed to fit them, generating variogram and cross-variogram models of ${\gamma }_{FF}$, ${\gamma }_{PP}$, ${\gamma }_{FP}$, and ${\gamma }_{PF}$. Exponential, Gaussian, and spherical models are widely used (Wu and Murray 2005; Adhikary, Muttil, and Yilmaz 2017). In this study, the experimental variograms or cross-variograms are fitted with these three models, and the fitted models with the minimum residual sum of squares (RSS) are selected.

Based on the fitted variogram and cross-variogram models, the Co-Kriging weights $w_{i_1}$ and $w_{i_2}$ can be derived from a set of $(n+2)$ linear equations (Goovaerts 1997), which are expressed by: (5) $\begin{array}{rl}& \sum _{i_1=1}^n{\gamma }_{FF}(L_{i_1}-L_{j_1})w_{i_1}+\sum _{i_2=1}^m{\gamma }_{FP}(L_{i_2}-L_{j_1})w_{i_2}+{\mu }_1={\gamma }_{FF}(L_{j_1}-L_0)\\ & \sum _{i_1=1}^n{\gamma }_{PF}(L_{i_1}-L_{j_2})w_{i_1}+\sum _{i_2=1}^m{\gamma }_{PP}(L_{i_2}-L_{j_2})w_{i_2}+{\mu }_2={\gamma }_{PF}(L_{j_2}-L_0)\\ & \mathrm{for}{j}_1=1,\dots \dots ,n{j}_2=1,\dots \dots ,m\end{array}$(5) where ${\mu }_1$ and ${\mu }_2$ are Lagrange multiplier parameters corresponding to the two unbiased conditions. Further details about Co-Kriging can be found in Goovaerts (1997) and Journel and Huijbregts (1978).

Sea ice drift is a 2D vector with velocity $v$ and direction $d$. In this study, $v$ and $d$ are estimated with Co-Kriging separately. The spatial resolution of the fused sea ice drifts resulting from Co-Kriging is set to 1 km. The choice of a 1 km resolution in this study is motivated by two main factors: scale and application. In terms of scale, the output resolution of the fusion (1 km) is roughly in the middle of two input resolutions, with the pattern-matching results downscaled by a factor of 10 (from 10 km to 1 km) and the feature tracking results upscaled by a factor of around 10 (from 120 m to 1 km). It is worth noting that the output feature tracking results are assumed to maintain a consistent resolution with the input SAR image resolution, which is 120 m. While feature tracking algorithm can theoretically achieve higher-resolution (i.e. subpixel-level) sea ice drift results, inconsistencies in the subpixel levels of different feature points in SAR images make it challenging to accurately estimate the resolution of the output feature tracking results. From an application perspective, although the resolution of sea ice models is increasing and can reach 1 km (Hutter, Losch, and Menemenlis 2018), the resolution of most satellite-based sea ice drift products still remains at the 10 km level (Wang et al. 2022). Thus, 1 km satellite-based sea ice drift information is needed to evaluate and assimilate the models.

3.3. Validation

Acquisition times of the SAR image pairs are utilized to select temporally closest buoys within 12 h for validation. For instance, IABP buoy observations from 2020-11-12 12:00–2020-11-13 12:00 are chosen to validate sea ice drifts extracted from SAR image pairs acquired on 2020-11-12 13:47:37 and 2020-11-13 14:28:47. Some GPS records may be inaccurate due to instrument deployment issues or damage (Hwang 2013). In this study, velocity records exceeding 60 km/d are excluded from the validation set, as sea ice typically does not move at such high speeds.

For validation purposes, the interpolation method is used (Wang et al. 2022). Sea ice drift vectors within a 25 km radius surrounding a buoy are first interpolated to the buoy's starting position. Subsequently, the root mean square error (RMSE) is calculated by comparing the velocities and azimuth angles (i.e. direction) of the buoys and interpolated vectors in the geocentric coordinate system.

4. Results and discussion

4.1. Results of the feature-tracking sea ice drift

Figure 4 illustrates the initial feature-tracking drift vectors (Figure 4 (a)) and vectors after bad match filtering (Figure 4 (b)), utilizing an image pair from Region III on November 20 15:07:16 (UTC) and 21 17:26:24 (UTC), 2020 as a representative case. To reduce visual clutter caused by dense vector distribution, the vector density in Figure 4 is decreased. Figure 4 reveals the presence of many mismatches among initial feature-tracking vectors (Figure 4 (a)). Following filtering, isolated vectors and those inconsistent with their neighbors are effectively discarded, while most properly matched vectors are retained (Figure 4 (b)), demonstrating the efficacy of the filtering method employed in this study.

Figure 4. Sea ice drift vectors (density reduction ratio: 1:20): (a) initial feature-tracking vectors from SAR data and (b) vectors from SAR data after bad match filtering.

Table 2 provides accuracy assessment results for the DTU product and the retrieved feature-tracking sea ice drift, compared to buoy data, across regions I-V. As indicated by Table 2, the retrieved feature-tracking sea ice drift displays lower Velocity RMSE (0.09-1.53 km/d less) and Direction RMSE (1.56°−24.26° less) than the DTU product for all regions. The largest differences are observed in regions I (Velocity RMSE) and V (Direction RMSE). Overall, the retrieved feature-tracking sea ice drift exhibits a Velocity RMSE of 1.39 km/d and a Direction RMSE of 10.21°, which are 0.70 km/d and 7.59° lower than the DTU product, respectively. This implies that the accuracy of the feature-tracking vectors surpasses that of the pattern-matching vectors, when utilizing the same SAR data. Consequently, fusing feature-tracking retrieved sea ice drift with pattern-matching sea ice drift could enhance the accuracy of the latter. Furthermore, Table 2 indicates that the velocity RMSE of Region IV exceeds that of other regions, regardless of whether the accuracy assessment is based on DTU products or feature-tracking results. This disparity may be attributed to two factors. Firstly, the data in Region IV corresponds to the melting season, during which the SAR backscattering is attenuated by the wet ice cover, thereby may cause a decline in the accuracy of sea ice drift extraction (Karvonen 2012). Secondly, the sea ice drift velocity in Region IV is comparatively higher, which may lead to a larger velocity RMSE (Wang et al. 2022).

Table 2. Accuracy assessment results for the DTU product and the retrieved feature-tracking sea ice drift in comparison with buoys.

Region	DTU Product		Feature-tracking Sea Ice Drift
	Velocity RMSE (km/d)	Direction RMSE (°)	Velocity RMSE (km/d)	Direction RMSE (°)
I	2.50	20.92	0.97	13.64
II	1.29	8.62	1.20	7.06
III	1.65	13.56	1.37	11.06
IV	3.04	12.40	2.25	10.24
V	1.94	33.35	1.18	9.09
All	2.09	17.80	1.39	10.21

4.2. Results of the refined sea ice drift

Figure 5 presents examples of sea ice drift fields from the DTU product, retrieved feature-tracking results, and refined results in the five experimental regions. All vectors in Figure 5 are shown at the same resolution (i.e. 30 km) for comparison. In Figure 5, the SAR image pairs utilized in the DTU product and for retrieving feature-tracking sea ice drift are as follows: December 16 16:33:03 (UTC) and 17 16:26:00 (UTC), 2020 for Region I; January 03 06:35:26 (UTC) and 04 02:22:07 (UTC), 2020 for Region II; November 14 17:34:48 (UTC) and 15 18:16:22 (UTC), 2020 for Region III; May 15 05:04:17 (UTC) and 16 06:21:08 (UTC), 2020 for Region IV; and March 20 18:20:16 (UTC) and 21 18:12:51 (UTC), 2020 for Region V. As seen in Figure 5, the DTU product (Figure 5 (b)) contains some mismatched vectors (i.e. vectors showing relatively large speed and angle differences from its neighboring vectors), possibly due to the noise effect of SAR images. The feature-tracking vectors show good neighborhood consistency (Figure 5(a)), but uneven distribution (Figure 5(c)). After fusing with feature-tracking vectors, the refined results exhibit higher accuracy than the DTU product (i.e. mismatched vectors are corrected) (Figure 5 (d)). Furthermore, some DTU vectors (e.g. the orange vectors in the top left of Region II in Figure 5 (b)), exhibiting velocity differences compared to the feature-tracking vectors (Region II in Figure 5 (a)), are rectified following refinement (Region II in Figure 5 (d)).

Figure 5. Sea ice drift vectors in the five experimental regions: (a) retrieved feature tracking results, (b) DTU product, and (d) refined results; (c) spatial distribution of DTU product and feature-tracking vectors.

Table 3 displays the accuracy assessment results for the refined sea ice drift, compared with buoy data, across regions I-V. The buoy validation times in different regions are also listed in Table 3, totaling 7,908 validation times for all regions. A comparison between Tables 2 and 3 reveals that the fused results outperform the DTU product in accuracy. The overall velocity and direction RMSEs of the refined sea ice drift results decrease by 0.47 km/d (22%) and 4.97° (28%) compared to the DTU product, with velocity RMSE reductions in different regions ranging from 0.01-1.45 km/d (1%−58%) and direction RMSE reductions between 0.68-17.42° (7%−52%). Furthermore, in terms of RMSE reduction, regions I and V exhibit better refinement performance than the other three regions. Figure 6 displays boxplots showing the absolute errors of drift velocity and direction for both the DTU product and the refined results. The boxplots reveal that, in terms of median, lower quartile, upper quartile, and maximum values of absolute errors for both velocity and direction, the refined results exhibit a higher degree of accuracy than the DTU product. The success of the proposed framework in enhancing the accuracy of the pattern-matching sea ice drift from SAR imagery is demonstrated by both Table 3 and Figure 6.

Figure 6. Boxplots depicting absolute errors of drift velocity (a) and direction (b) for the DTU product and the refined results. The upper and lower quartiles are represented by the right and left boundaries of the ‘box’, respectively. The maximum and minimum values are represented by the right and left lines, respectively. The medians are shown as red lines.

Table 3. Accuracy assessment results for the refined sea ice drift in comparison with buoys.

Region	RMSE		RMSE Reduction (Compared to DTU product)		Buoy Validation Times
	Velocity (km/d)	Direction (°)	Velocity (km/d) / Percentage	Direction (°) / Percentage
I	1.05	16.43	1.45 / 58%	4.49 / 21%	267
II	1.28	7.94	0.01 / 1%	0.68 / 8%	3,209
III	1.58	12.57	0.07 / 4%	0.99 / 7%	268
IV	2.74	11.29	0.30 / 10%	1.11 / 9%	4,107
V	1.42	15.93	0.52 / 27%	17.42 / 52%	57
All	1.62	12.83	0.47 / 22%	4.97 / 28%	7,908

Figure 7 presents scatterplots of buoy displacements versus sea ice drift vectors retrieved from the DTU product and the refined results, with r representing the correlation coefficient, ranging from −1 to 1. A value of r closer to 1 indicates greater similarity between sea ice drift vectors and buoy displacements. As observed, all r values increase after refinement. Specifically, the r values between the DTU sea ice drift product and buoy displacements are 0.818 for velocity and 0.989 for direction. After fusing with the feature-tracking vectors, the r values become 0.919 for velocity and 0.997 for direction, increasing by 0.101 and 0.008, respectively. A comparison of Figure 7(a) and (b), as well as 7(c) and 7(d), reveals more scatter points clustering along the one-to-one line with reduced dispersion after fusing the DTU product with the feature-tracking vectors. In particular, most overestimated sea ice drift velocities in the DTU product (i.e. scatter points below the one-to-one line in Figure 7(a)) are corrected after refinement (Figure 7(b)). The results displayed in Figure 7 show that the proposed framework improves the accuracy of the DTU product, indicating its effectiveness.

Figure 7. Comparison of sea ice drift vectors with buoy displacements: (a) density scatterplot of buoy-based ice velocity vs. DTU ice velocity, (b) density scatterplot of buoy-based ice velocity vs. refined ice velocity, (c) density scatterplot of buoy-based ice drift direction vs. DTU ice drift direction, and (d) density scatterplot of buoy-based ice drift direction vs. refined ice drift direction.

The previous studies indicate that in the pattern-matching sea ice drift vectors, velocity error increases with ice velocity, and large direction errors are predominantly found in the very low ice velocity regime (Hwang 2013; Wang et al. 2022). Consequently, this study further examines the relationship between errors and ice velocity, with results shown in Figure 8 and Figure 9. Figure 8 reveals that the DTU product exhibits numerous vectors with large velocity and angle absolute errors in the low ice velocity regime; however, these large errors are eliminated in the refined results (as seen in the scatter points within the red circles in Figure 8 (a)-(d)). Figure 9 demonstrates that velocity and direction RMSEs decrease in all ice velocity sections after refinement. Moreover, when ice velocity is below 8 km/d, the differences in velocity and direction RMSEs between the refined results and the DTU product are notably large (Figure 9 (a) and (c)). Figure 8 and Figure 9 indicate that the DTU product contains some large errors in the low ice velocity regime, and the proposed framework effectively reduces these errors. Additionally, the proposed framework exhibits superior refinement performance in the low ice velocity regime (below 8 km/d) compared to the high ice velocity regime (8-22 km/d). This may be attributed to the DTU product employing the Maximum Cross-Correlation pattern matching algorithm, which is unable to effectively track sea ice drift at a sub-pixel scale (Hwang 2013). In contrast, feature tracking methods can retrieve sea ice drift at a sub-pixel level, thus effectively tracking sea ice drift with short displacements (Muckenhuber and Sandven 2017).

Figure 8. Scatterplots between absolute errors and buoy-based ice velocity: (a) scatterplot of the absolute error of the DTU ice velocity against the buoy-based ice velocity, (b) scatterplot of the absolute error of the refined ice velocity against the buoy-based ice velocity, (c) scatterplot of the absolute error of the DTU ice drift direction against the buoy-based ice velocity, and (d) scatterplot of the absolute error of the refined ice drift direction against the buoy-based ice velocity.

Figure 9. RMSEs in different buoy velocity sections: velocity RMSEs for the DTU product and the refined results in different buoy velocity sections ((a) below 22 km/d and (b) above 22 km/d), and direction RMSEs for the DTU product and the refined results in different buoy velocity sections ((c) below 22 km/d and (d) above 22 km/d).

4.3. Discussion

The experimental findings in the previous subsections demonstrate the capability of the proposed framework to enhance the accuracy and resolution of the pattern-matching sea ice drift from SAR imagery. Notably, some mismatched vectors are identified in the pattern-matching vectors, and several large errors exist in the low ice velocity regime. The proposed framework effectively corrects these mismatched vectors and reduces the large errors.

The novelties of the proposed framework can be summarized as follows: Firstly, the proposed framework fuses feature-tracking vectors with pattern-matching vectors to improve the accuracy and resolution of the latter. It combines pattern matching and feature tracking methods in a different way compared to existing methods. Previous combined methods utilize unevenly distributed feature-tracking vectors for interpolation, which may result in low accuracy in areas where feature-tracking vectors are sparsely distributed, despite the use of pattern matching algorithms to adjust the interpolated vectors. In contrast, the proposed framework uses evenly distributed pattern-matching vectors as the background field and adjusts them through fusing with feature-tracking vectors. This approach generally ensures that the accuracy in any area is not lower than the accuracy of the pattern-matching result. Secondly, it introduces a general methodological framework, applicable not only to refining pattern-matching sea ice drift from Sentinel-1 C-band SAR imagery, but also from other C-band data (e.g. RADARSAT-2), as well as X-band (e.g. TerraSAR-X) and L-band (e.g. ALOS PALSAR) SAR datasets. Furthermore, combining pattern-matching or feature-tracking sea ice drift information from multi-frequency SAR imagery may further enhance the accuracy of the final fused result, as the different frequency bands complement each other in sea ice drift retrieval (Lehtiranta, Siiriä, and Karvonen 2015). Moreover, the proposed framework is not limited to improving pattern-matching sea ice drift solely based on SAR imagery, but it can also be applied to other types of remotely sensed data, such as passive microwave and microwave scatterometer data. However, it should be noted that passive microwave and microwave scatterometer data typically exhibit coarse spatial resolution. Consequently, fine-resolution feature-tracking vectors cannot be retrieved from these coarse datasets. Some fine-resolution SAR (e.g. Sentinel-1 SAR and RADARSAT-2) and optical (e.g. AVHRR and MODIS) datasets (Petrou and Tian 2017; Li et al. 2022a; Fang et al. 2023) may be employed to retrieve feature-tracking vectors for refining the pattern-matching fields from passive microwave and microwave scatterometer data. Nevertheless, temporal differences, as well as differences in spatial resolutions and geographical errors among multi-sensor data, should be considered, and the cloud effect on optical sensors must be addressed.

Although the proposed framework demonstrates good performance in refining the pattern-matching sea ice drift from SAR imagery, some issues are observed in the experimental results. Figure 7 indicates that most absolute errors in drift velocity and direction for both the DTU product and the refined results are comparable. Two possible reasons can be considered. First, as Figure 9 reveals that the proposed framework exhibits better refinement performance in low ice velocity regimes than in high ice velocity regimes, it is possible that the majority of the buoys used for validation are in high ice velocity regimes. As shown in Figure 1 and Table 1, most buoys (287 out of 359) are situated in Region II and Region IV, where the sea ice exhibits relatively faster movement compared to other regions. Consequently, the buoy velocity distributions are further analyzed, with the results displayed in Figure 10. From Figure 10, it becomes evident that the velocity of most buoys used for validation is relatively high (i.e. buoys with velocities greater than 8 km/d account for 74.5% of the total number of buoys), and in the high ice velocity regime, there is not a significant difference in error between the DTU product and the refined results. Second, Figure 5 shows that mismatched vectors are corrected in the refined results, but there are no buoys in some of those mismatched areas (e.g. Region III-V in Figure 5).

Figure 10. Frequency histogram of buoy-based ice velocity: (a) below 22 km/d and (b) above 22 km/d.

Moreover, Figure 9 illustrates that the buoy-based ice velocity greater than 22 km/d yields significantly higher velocity and direction RMSEs for both the DTU product and the refined results compared to those below 22 km/d. Therefore, validation is further conducted for the aforementioned velocity regime. The results show that the velocity and direction RMSEs for the DTU product are 12.57 km/d and 40.83°, whereas the velocity and direction RMSEs for the refined vectors are 12.40 km/d and 23.30°. These values are much larger than the overall velocity and direction RMSEs presented in Tables 2 and 3 (i.e. 2.09 km/d and 17.80° for the DTU product; 1.62 km/d and 12.83° for the refined vectors). As observed in Figure 10, the validation times using buoys with ice velocities exceeding 22 km/d are quite small (58 out of 7841). These validations may be conducted on the mismatched vectors in the DTU product. As shown in Figure 5 (b), the mismatched vectors all exhibit high velocity (i.e. red vectors).

Furthermore, Table 3 reveals that the refinement effect varies across different regions. Compared to the other three regions, Regions I and V exhibit superior refinement performance. This disparity may be attributed to the differing refinement performances in low and high ice velocity regimes. To investigate this, the buoy velocity distributions in various regions are analyzed, with the results presented in Figure 11. The mean values of buoy ice velocity are 7.31, 9.61, 9.78, 13.40, and 6.51 km/d for Regions I to V, respectively, indicating that the mean values of buoy ice velocity for Regions I and V are lower than those in other regions. Moreover, Figure 11 reveals that the majority of buoy ice velocities in Regions I and V are below 8 km/d, falling within the ice velocity regime where refinement performance is relatively good.

Figure 11. Boxplots illustrating buoy ice velocities in different regions. The upper and lower quartiles are represented by the right and left boundaries of the ‘box’, respectively. The maximum and minimum values are represented by the right and left lines, respectively. The medians are shown as red lines. Green triangles signify the means of buoy ice velocities.

Additionally, several other issues are identified during the experimental process. First, the AKAZE-based feature tracking relies on handcrafted local feature descriptors and traditional similarity measurement methods. Recent research has demonstrated that employing deep learning techniques for deep feature extraction and similarity measurement can yield improved matching performance (Ma et al. 2021). Consequently, future studies could consider applying deep neural networks in the generation of feature-tracking vectors. Second, In this study, the resolution of the fused result is set to 1 km, and the neighboring scope for bad matches filtering and fusion is set to 25 km. Further studies are needed to investigate the impact of the output fusion resolution and neighboring scope on accuracy. Third, the Co-Kriging-based fusion accounts solely for spatial correlation, neglecting temporal correlation, which could also provide valuable information for vector fusing to enhance the accuracy of the refined drifts. As such, future research could explore the application of spatiotemporal information for fusion.

5. Conclusions

This study proposes a framework aimed at enhancing the resolution and accuracy of pattern-matching sea ice drift derived from SAR imagery. Specifically, feature-tracking vectors are first extracted from SAR data through the use of an A-KAZE-based feature detection method and a Brute-Force matcher-based feature matching algorithm, with an added process for filtering out mismatched vectors. Then, a Co-Kriging-based technique is applied to refine the pattern-matching sea ice drift by fusing it with the feature-tracking vectors. The DTU product is used as the pattern-matching sea ice drift. Five regions are selected as case studies, with the refined vectors evaluated using buoy GPS displacements and compared to the DTU product. The results demonstrate that the refined displacement field outperforms the DTU product in terms of drift accuracy, with overall velocity and direction RMSE reductions of 0.47 km/d (22%) and 4.97° (28%), respectively, and enhanced spatial resolution (from 10 km to 1 km). This highlights the effectiveness of the proposed framework. Notably, the pattern-matching vectors from SAR imagery contain some large errors in the low ice velocity regime and several mismatched vectors; the proposed framework effectively reduces these large errors and corrects the mismatched vectors. The proposed framework presents two novel aspects, combining pattern-matching and feature-tracking methods in a distinctive way from existing combined approaches, and offering a general methodological approach for refining pattern-matching sea ice drift derived from not only SAR imagery but also other remotely sensed data types (e.g. passive microwave and microwave scatterometer data). Future work will involve using optical and SAR datasets, as well as published sea ice drift products, to generate daily sea ice drift fields across the entire polar regions. Additionally, efforts will be made to develop deep learning techniques for sea ice drift tracking, to examine the impact of the output fusion resolution and neighboring scope on accuracy, and to include spatiotemporal correlation into the fusion.

Acknowledgements

The authors would like to acknowledge the open-source A-KAZE code (https://github.com/pablofdezalc/akaze) and BF matcher code (https://github.com/opencv/opencv/blob/4.x/modules/features2d/src/opencl/brute_force_match.cl). The authors would like to thank two anonymous reviewers for their constructive comments and suggestions, which greatly improved the quality of the manuscript.

Data availability statement

Data will be made available on request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research is supported by the National Key Research and Development Program of China [grant no 2021YFC2800705], the National Natural Science Foundation of China [grant no 42206247], Guang Dong Basic and Applied Basic Research Foundation [grant no 2021A1515110779] and Fengyun Application Pioneering Project [grant no FY-APP-2022.0201].