Comparison of Doppler-Derived Sea Ice Radial Surface Velocity Measurement Methods From Sentinel-1A IW Data

The near-instantaneous radial velocity of a target can be obtained using the Doppler effect of synthetic aperture radar (SAR), which is widely used in ocean current retrieval. However, in sea ice drift velocity measurements, only a Doppler centroid estimation algorithm in frequency domain has been studied, so whether there is a better algorithm is worth exploring. In this article, based on Sentinel-1A interferometric wide data, three Doppler centroid estimation algorithms applied to ocean current retrieval are selected. Combined with the characteristics of the Terrain Observation by Progressive Scans mode, made two applicability adjustments to each algorithm, and finally applied the three algorithms to sea ice radial surface velocity measurements. The first adjustment is to explore and determine the optimal parameters. The second adjustment is to use parallel computing to improve the efficiency, which is improved by an average of 43.55%. In addition, the deviation of Doppler centroid estimation bias correction is verified using rainforest data, and the deviation is controlled at approximately 3 Hz. Based on the three algorithms, five sets of experiments are carried out in this article. By analyzing and comparing the results of each algorithm, it is found that the results of the three algorithms are relatively consistent, among which the correlation Doppler estimation algorithm has the advantages of high efficiency and high precision and is the most suitable method for sea ice drift measurement among the three methods. However, for SAR images with abnormal speckles caused by human activities, the sign Doppler estimation algorithm can effectively remove abnormal speckles and ensure the smoothness of the image with better adaptability.


I. INTRODUCTION
I NCREASING global warming has led to accelerated Arctic sea ice melting, which continues to increase the navigation Manuscript  potential of the Arctic ice area due to the significant reduction in sea ice coverage [1]. However, complex and changeable sea ice will drift driven by wind and currents, increasing the risk to ship navigation and offshore operations [2]. At the same time, sea ice extrusion can cause the collapse of oil platforms, ship blockage, and other sea ice disasters [3]. Therefore, it is very important to study sea ice drift monitoring methods and realize sea ice drift retrieval and prediction with high spatial and temporal resolutions. Remote sensing monitoring is currently the most popular sea ice drift monitoring method. Synthetic aperture radar (SAR), which can penetrate clouds and fog, is not affected by polar bad weather and provides all-weather and all-day monitoring. It has become the most commonly used sensing method for sea ice drift monitoring. However, SAR has some limitations, including a long revisit period. If a common time-series method, such as feature tracking or pattern matching techniques, is used to measure sea ice drift, the change in sea ice will make the measurement more difficult due to a revisit period of several to 10 days [4], and it is difficult to meet the needs in emergency situations. If the requirements for temporal resolution are increased, there will be insufficient data. In the study of sea surface velocity retrieval based on SAR data, the commonly used Doppler centroid estimation algorithm can obtain the near-instantaneous radial surface velocity (RVL) of sea currents using a single SAR image [5], which makes up for the shortcomings of low SAR temporal resolution. Although its spatial resolution is not as high as that of time-series method, it is sufficient to meet the requirement of spatial resolution for large-scale sea ice drift. Therefore, this article attempts to apply this algorithm to RVL measurements of sea ice drift.
Madsen [6] proposed two time-domain algorithms, correlation Doppler estimation (CDE), and sign Doppler estimation (SDE), based on a review of frequency-domain Doppler centroid estimation algorithms, and evaluated the application of the three algorithms. Kim et al. [5] used the CDE algorithm to retrieve the sea surface velocity near Jeju Island based on raw ERS-1 SAR data and confirmed that it was consistent with the observed in situ data within the error bound limits. Since then, an increasing number of researchers have begun to retrieve the sea surface current field based on the Doppler centroid estimation algorithm [7], [8], [9], [10], [11], [12], [13], [14]. Before Sentinel-1 was deployed, most research was based on Envisat or Radarsat data acquired using the Stripmap and ScanSAR scanning methods. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ After Sentinel-1 was put into deployed, the Sentinel-1 Level 2 Ocean (OCN) RVL product has used to generate high-resolution sea surface current maps, which are conducive to ocean observation and modeling, and have become a data source of interest to researchers. For example, Johnson et al. [15] studied the sea surface current field based on Sentinel-1A (S1A) interferometric wide (IW) and extra wide (EW) mode RVL products. Due to the bias caused by platform effects and instrument anomalies contained in the Level 2 RVL product, the sea surface current retrieval based on Sentinel-1 Level 2 RVL data has shifted to the removal of wind and wave bias and the comparative evaluation with other sensors [16], [17], [18], [19], and the research and exploration of the algorithms will gradually be reduced.
In the development of Doppler centroid estimation algorithms, the measurement of sea ice drift RVL has become a major application. Hansen et al. [20] proposed that the Doppler centroid estimation algorithm can be used for sea ice drift, and compared it with traditional feature tracking methods. In 2015, based on Envisat Stripmap Mode (SM) data, Kraemer et al. [21] calculated the radial sea ice surface velocity in the Fram Strait using a Doppler centroid estimation algorithm and corrected the azimuth and range bias. Three years later, Kraemer et al. provided the theoretical background for Doppler centroid estimation based on S1A EW mode data. For the first time, they quantitatively compared the ice drift from the Doppler frequency shift with the drift derived using time-series methods over comparable time scales and pointed out that more accurate results may be obtained using IW data with higher resolution [22]. Based on Envisat ScanSAR data, Park et al. [23] calculated sea ice drift velocity using a combination of feature tracking and Doppler centroid estimation in 2019 and concluded that it is highly consistent with the root mean square error (RMSE) of the time-averaged velocity.
However, overall, there are few studies on sea ice drift measurements using the Doppler centroid estimation algorithm, and there are few studies based on Sentinel-1 data. At the same time, different from the various Doppler centroid estimation algorithms for sea surface current field retrieval research, sea ice drift measurement research is only based on a frequency-domain algorithm. Therefore, it is necessary to try additional Doppler centroid estimation algorithms for sea ice drift measurement and analyze the optimal algorithm selection. Therefore, this article attempts to apply three Doppler centroid estimation algorithms for sea surface velocity retrieval to sea ice RVL measurements. The three algorithms are the chirp phase angle (CPA) algorithm in frequency domain, CDE algorithm and SDE algorithm in time domain. Previously, the three algorithms have been based on ScanSAR mode data, while the higher precision IW data mentioned by Kraemer et al. [22] is Terrain Observation by Progressive Scans (TOPS) mode data. Therefore, research using IW data needs to adjust the algorithm according to different scanning methods to make it applicable.
In summary, based on the Level-1 single-look complex (SLC) image of the S1A IW mode, the CPA algorithm, CDE algorithm, and SDE algorithm are used to estimate the Doppler centroid frequency shift to achieve sea ice RVL measurement. In this process, parallel computing is used to improve the algorithm efficiency. Through comparative analysis and accuracy evaluation, the best parameter combination, bias correction deviation, and algorithm selection suggestions of the three algorithms for IW mode data are determined, providing detailed and complete algorithm reference and selection suggestions for subsequent research.

A. Sentinel-1A Data
The SAR data used in this article come from the S1A satellite launched in April 2014. The satellite is in a sun synchronous orbit. The single satellite revisit period is 12 days. The antenna faces the right side of the flight direction. It operates in the Cband with horizontal and vertical polarizations. S1A operates in different modes in different regions, with SM, IW, EW, and wave mode to meet the requirements of users for different spatial resolutions and image sizes.
In this article, IW mode data are selected for the following reasons: 1) The acquisition of IW mode data adopts the TOPS mode. The TOPS mode avoids a decrease in inversion accuracy due to antenna gain in SM because the antenna beam sweeps forward along the azimuth direction. 2) Among the four S1A modes, the IW and EW data are obtained in TOPS mode. The width of the IW mode reaches 250 km, and the width of the EW mode reaches 410 km, which is much larger than the width of 80 km in the SM mode.
3) The resolution of the Level 1 SLC data in the IW mode is 2.7 × 22 to 3.5 × 22 m, while the resolution of the EW mode is 7.9 × 43 to 15 × 43 m. With the average drift velocity of 3-20 cm/s over the Arctic [29], the EW model will not be able to obtain Doppler estimates with the required accuracy in most cases [21]. Therefore, the IW SLC data have a higher spatial resolution, which is helpful for obtaining a more refined sea ice RVL.

B. NSIDC Sea Ice Drift Velocity Products
The sea ice motion dataset used for verification is the Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors, Version 4 dataset from the National Snow and Ice Data Center (NSIDC). The dataset contains daily and weekly sea ice motion vectors as well as browse images representing the weekly data. The daily and weekly sea ice motion data are georeferenced to the Northern and Southern 25 km original EASE-Grid projections. The size of the northern EASE-Grid is 361 × 361 cells, and the lower left geolocation point of the grid is 29.71270°N, 45°W. In this article, the daily sea ice motion vector products are used. The sea ice motion is stored in the dataset with u and v components corresponding to the x and y directions of the coordinates, respectively.

C. Buoy Data
The International Arctic Buoy Programme (IABP) has deployed many buoys in the Arctic Ocean to provide meteorological and oceanographic data for real-time operational requirements and research purposes including support to the World Climate Research Programme and the World Weather Watch Programme. The buoy data used in this article are from a buoy with an ID of 30023466412770 deployed by the IABP. This buoy continually logged its GPS position at 15-min intervals, and the sea ice drift average velocity every 15 min was calculated using the change in buoy position.

A. Data Preprocessing Method
Under the ScanSAR or SM, the antenna pointing is fixed, and the Doppler frequency shift is constant, but the IW mode data provided by the S1A satellite used in this article are collected using the TOPS mode in which the antenna beam sweeps forward along the azimuth direction in each burst, resulting in a linear change in the Doppler centroid with the azimuth time. This linear frequency modulation with time can be removed through deramping. A function with opposite symbols and the same frequency modulation and Doppler centroid rate is established to obtain the deramped phase, and the raw SLC complex data are multiplied by the deramped phase to complete the deramping. This is one of the TOPS mode data preprocessing steps, and the specific steps can be found in [24].
In addition, since the antenna beam is switched back and forth between the three subswaths in TOPS mode, there is overlapping and invalid information filling between bursts, so there will be azimuth discontinuity and black band problems during stitching. The overlapping and invalid information between bursts will affect the continuity and accuracy of the bursts in the Doppler centroid estimation results. To remove the invalid and overlapping information and merge the bursts in subswath based on the azimuth time, it is necessary to perform deburst processing on the data. In this article, the overlapping azimuth position is determined using the start time and last time of each burst in the metadata and the azimuth time interval. The invalid information of the overlapping part can be completed by subtracting the average difference between the two bursts from the pixel value of the azimuth time corresponding to the adjacent burst. The effective information of the overlapping part, it can be retained by the pixel average, and finally a continuous and nonoverlapping image in the azimuth direction is obtained.
In summary, preprocessing considers the difference between the TOPS mode and the ScanSAR or SM. Before Doppler centroid estimation, data deramping and deburst processing are performed for the unique IW data imaging mode.

B. Radial Surface Velocity Calculation Method
The measurement of the sea ice RVL is usually based on Doppler centroid estimation, because according to the Doppler effect of SAR, the Doppler centroid frequency shift f dm (in Hz) has the following relationship with the effective relative velocity v r between the SAR instrument and the Earth's surface : where λ is the carrier wavelength of the SAR system. For moving scatterers on the Earth's surface, the Doppler centroid shift is affected by many factors, including geophysical Doppler shift f phy , predicted geometric Doppler shift f geom , and the bias caused by scalloping and antenna electronic mispointing f bias . In this article, the overall Doppler centroid shift is called the measured Doppler centroid f dm In this article, three algorithms are selected to estimate the measured Doppler centroid f dm , including the CPA algorithm in frequency domain, CDE algorithm and SDE algorithm in time domain. The three algorithms are usually used in the retrieval of the sea surface current field. In this article, they are applied to the measurement of the sea ice drift RVL.
The Doppler centroid frequency used to measure the sea ice drift RVL is the surface scatterer geophysical Doppler shift f phy .
To obtain an accurate f phy , it is necessary to remove the predicted geometric Doppler shift f geom and the known bias f bias contained in f dm and convert f phy to the geophysical RVL v rel using where θ is the radar incidence angle. The RVL is positive when the target moves close to the radar, and negative when the target moves away from the radar. In this article, the sea ice drift RVL calculated through f phy inversion using three algorithms is compared and analyzed. The three algorithms are introduced next.

C. Measured Doppler Centroid Estimation Algorithm
Common measured Doppler centroid estimation algorithms are divided into frequency-domain and time-domain algorithms. This section describes the CPA algorithm in frequency domain and CDE algorithm and SDE algorithm in time domain. These three algorithms are based on ScanSAR model data to retrieve the sea surface current field. For the TOPS model data used in this article, some adaptive adjustments and parameter selection are needed.
The three algorithms use the direct relationship between the phase of the correlation coefficient and the Doppler centroid frequency to calculate the Doppler centroid frequency, so the different methods to obtain the correlation coefficient are the main difference among the three algorithms.
1) CPA Algorithm: For an energy-limited signal, the energy of the signal in the time domain is equal to the energy of the signal in the frequency domain; that is, the total energy of the signal remains unchanged through a Fourier transform [25]. Therefore, the algorithm first performs a Fourier transform on each image block, converts the local image to the frequency domain, and then obtains its azimuth power spectrum P (η) through (Fourier transform square/interval length). Then, uses the Fourier transform dual relationship between the correlation function and the power spectrum to find the correlation function where N is the azimuth length of the block and FFT represents the Fourier transform of P (η). Hou [25] corresponding the FFT result to [−PRF/2, PRF/2] through FFT shift function, estimated the azimuth energy of the result processed by FFT shift function.
Considering that the energy integrals on the left and right sides of the Doppler centroid frequency are equal [26], Hou performs sinusoidal curve fitting on the azimuth energy spectrum to find the center point frequency as the Doppler centroid frequency. In this article, the method is adjusted, the frequency range is set as [0, (N−1) * PRF/N], the FFT shift transform and sine curve fitting are removed, and the correlation coefficient corresponding to the first harmonic of the correlation function is directly used as the correlation coefficient corresponding to the Doppler centroid frequency.
2) CDE Algorithm: Doppler centroid estimation is also performed using the Fourier transform dual relationship between the power spectrum of the signal and the correlation function [6]. By defining a random process by x(n) = x (nT ), the correlation coefficient of the two samples with an interval of η can be expressed asR where x * (m) is the conjugate complex numbers of x(m). Since most SAR systems are only sampled slightly above the Nyquist rate, the correlation coefficient rapidly approaches zero as η increases. Therefore, η = 1 is often assumed in actual calculations.

3) SDE Algorithm:
This is like the idea of the CDE algorithm. The difference is that the algorithm can obtain the correlation function of a signal using only the symbol information of the signal [6]. When filtering circularly symmetric Gaussian processes in linear systems, the output is also a circularly symmetric Gaussian. However, due to the Gaussian property, the correlation properties of the process can be derived by analyzing the sign alone. The theoretical basis of the algorithm is the "arcsine law" of a Gaussian process [27], which describes the correlation of a Gaussian process. Because the SAR echo signal is a quasi-Gaussian process, the algorithm can be applied to SAR data.
Suppose x(t) and y(t) are two real Gaussian processes with zero mean. We can define a sign function such that the value of x and y greater than or equal to 0 is 1, and the value of x and y less than 0 is −1. By representing the SLC data of S1A as h(k) = I(k) + jQ(k), the sign correlation of the SLC data is estimated The "arcsine law" is used to derive the normalized correlation coefficients The complex correlation coefficient of the SLC data is further derivedρ

D. Bias Correction Methods
Bias is caused by many nongeophysical anomalies. The S1A data mainly include the predicted geometric Doppler shift f geom caused by the relative motion between the satellite platform and a stationary target on the surface of the rotating earth, the azimuth bias f scal caused by scalloping, the range bias f em caused by antenna electronic mispointing, and the residuals f caused by inaccurate estimation of the nongeophysical terms and unknown biases, where the residuals f cannot be removed.

1) Predicted Geometric Doppler Shift Estimation:
The predicted geometric Doppler shift is caused by the relative motion between the satellite platform and a stationary target on the surface of the rotating earth, which can be obtained from the stable orbit and attitude information of the satellite. The Doppler centroid prediction model needs to be established by (5) is the Doppler centroid polynomial coefficient, t 0 is the standard slant range time, and t is the slant range time image. The measured Doppler centroid frequency that removes the predicted geometric Doppler shift is called the Doppler centroid anomaly.
2) Scalloping Correction: Scalloping is a periodic signal in the azimuth direction in the Doppler centroid anomaly, which is caused by a scalloping signal. This hinders the analysis of the spatial variation in the Doppler centroid and affects the surface velocity, so it should be corrected. The data can be converted from the image domain to the frequency domain using a 2-D Fourier transform and corrected in the spectrum [17], [28]. Since scalloping is a bias in the azimuth direction, the dominant frequency is found along the central axis in the azimuth direction. Because most geophysical information is contained in the dominant frequency, a low-pass filter can be designed to filter out all frequencies greater than the first harmonic.

3) Antenna Electronic Mispointing Correction:
The degradation of the transmit/receive modules of the satellite antenna array over time will lead to a change in the Doppler centroid, resulting in an antenna electronic mispointing bias that varies with the range unit. The bias also varies with the polarization mode and the subswath. By calculating and analyzing the Doppler centroid anomaly in a scene with homogeneous backscatter, the antenna electronic mispointing bias can be estimated. Rainforest scenes are usually selected, but if the study area has no rainforest scenes that are close to it in time, the adjacent land area in the same-orbit or adjacent-orbit with relatively homogeneous backscattering can be selected as the reference data. If the image contains enough land pixels in the range direction, the land pixels can be directly used as the reference data.
First, processing is performed on the land scenes to obtain the Doppler centroid frequency after scalloping correction. Second, the azimuth mean of each range position is calculated to obtain the range distribution curve. Curve fitting is used to find the relationship between the range Doppler centroid and the incident angle. Then the relationship is substituted into the research area to calculate the antenna electronic mispointing bias. Finally, the  bias is subtracted from each pixel along the range direction to obtain an accurate geophysical Doppler frequency shift.

A. Experimental Data Selection
This article selects six sets of data, including the first five sets of experimental data and corresponding reference data acquired in the Bering Strait region, the Baffin Bay region, the Fram Strait region, and the Spitsbergen Islands region, respectively, and the sixth set of Amazon rainforest data for bias analysis. The important data parameters are shown in Table I, including data stop time, product unique ID, absolute orbit number, subswath number used in the experiment, pass, polarization, incident angle, data use, and research area. The reference data are the land data selected according to the principle of time closeness, space adjacent, and incidence angle similarity to the experimental data, which are used to correct the antenna electronic mispointing bias. Fig. 1 shows the intensity images of the five sets of sea ice experimental data.

B. Parameter Selection
According to the Doppler centroid estimation algorithm principle, each block calculates a Doppler centroid frequency. The block size determines the number of pixels used to estimate the Doppler centroid frequency, so it affects the computational efficiency of a single block and the smoothness of the Doppler centroid frequency to a certain extent. The sea ice movement trend in a certain range is relatively consistent. If a block is too small, it will cause RVL redundancy and Doppler centroid frequency roughness. If a block is too large, it cannot guarantee that the RVL is consistent within the block range, and it is easy to lose motion details. Therefore, it is very important to select an appropriate block size.
The step size controls the resolution of the Doppler centroid frequency. The shorter the step size is, the larger the dimension of the Doppler centroid image is. In the RVL calculation, the number of pixels that need to be interpolated will be smaller, and the resolution will be higher. However, step size that is too short will increase the number of iterations, which is not conducive to improve computational efficiency. Therefore, it is particularly important to select an appropriate step size.
An ocean area of 1 km × 1 km can be regarded as a single homogeneous target [30]. The number of pixels in the range direction after preprocessing is approximately twice that in the azimuth direction. Considering that the block size needs to be proportional to the input image size, the block size should meet the requirement of the length in the range direction being twice that in the azimuth direction. However, the resolution of the S1A SLC IW data used in the experiment is approximately 13.86 m × 2.33 m (azimuth × range), so it is difficult to achieve a resolution For the selection of the block and step size parameters, this article is based on geophysical Doppler shift estimation results and computational efficiency analysis. Taking the results of the first set of data under the CDE algorithm as an example, the geophysical Doppler frequency shift estimation results (local) obtained with different block sizes and step sizes are shown in Table II, the pixel mean and standard deviation of this result are also given in Hz.
Taking the results of the first set of data under the three Doppler centroid estimation algorithms as an example, the parallel computing time for different block sizes and step sizes is shown in Fig. 2.
According to the image analysis (Table II) and computational efficiency (Fig. 2), the block size has a great influence on the Doppler centroid frequency. At 64 × 128 pixels, although the overall computational time is shorter, the Doppler centroid frequency is rough and the standard deviation is the largest among the three, reaching approximately 27.7 Hz. At 256 × 512 pixels, the Doppler centroid frequency is smoother, and the standard deviation is only approximately 3.7 Hz, but it also loses many motion details. At the same time, it enlarges the influence of abnormal speckles caused by human activities, such as ships on the Doppler centroid frequency, and the computational time is much larger than that of the other two sizes. At 128 × 256 pixels, the image smoothness and detail retention are moderate, and the computational time is also maintained in the middle level, so it is the optimal block size.
This article uses 128 × 256 pixels to analyze the influence of the step size on the image and computational efficiency. The mean and standard deviation of the Doppler frequency of the three step sizes increase with increasing step size (Table II), and the image dimension decreases with increasing step size, which confirms that the smaller the step size is, the higher the accuracy. In addition, the difference between the standard deviations of step sizes of 2 × 4 and 4 × 8 pixels is 0.203, while the difference between step sizes of 2 × 4 and 8 × 16 pixels is 1.006. Combined with the computational time (Fig. 2), for a step size of 4 × 8 pixels relative to a step size of 2 × 4 pixels, the efficiency increased by 73.81%, and for a step size of 8 × 16 pixels relative to a step size of 2 × 4 pixels, the efficiency increased by 92.08%. In principle, a step size of 2 × 4 pixels is the best choice as the minimum step size, but its computational efficiency is too low, and it may be difficult to meet the computational time requirements in batch processing or emergency applications. The accuracies of step sizes of 2 × 4 and 4 × 8 pixels are not much different, but the computational efficiency is greatly improved. Therefore, this article selects 4 × 8 pixels as the optimal step size.

C. Parallel Computing
The Doppler centroid frequency measured in this article is calculated using a sliding block, and the calculation is carried out inside a single block in a sliding order. The 128 × 256 pixels block slides upward in the range with a step of 8 pixels and slides upward in the azimuth with a step of 4 pixels. The whole image is roughly calculated 6000 × 6000 times. If only sliding gradually, the number of iterations is too large and the computational efficiency is low.
Parallel computing is a computing mode in which many instructions are performed at the same time. The idea is to decompose the calculation process into smaller pieces under the premise of simultaneous execution and solve it in a concurrent manner. That is, according to the number of CPU cores, the processor cores participating in the calculation are established and simultaneously perform several processes to reduce the processing delay. While greatly improving the computational efficiency, it will not affect the result of Doppler centroid estimation.
However, to achieve calculation results that are completely consistent with the single machine calculation results, the entire preprocessed SLC image cannot be simple blocked considered in each block. According to the block sliding law, before computation based on the principle of parallel computing. In this article, the block and step sizes are the image was divided into 5 × 5 image blocks in the way shown in Fig. 3 to ensure that Fig. 3. Block diagram. The 8 × 9 color blocks represent the sliding window, the area circled by the red box and the blue box represents the two adjacent blocks when the block is divided, the red arrow points to the last sample of the first block, and the blue arrow points to the first sample of the second block. the subsequent accurate splicing can be achieved, and the result is completely consistent with the single machine calculation. In the nonboundary part, the first and last lines/samples of each block are calculated as follows (for example): In the formula, Samples represents the overall sample number of SLC images (line corresponding to azimuth, sample corresponding to range), Samples block represents the block range length, x step represents the range step, Samples num represents the number of blocks in the range direction, and m represents the block range cell.
The purpose of parallel computing is to improve the computational efficiency. Based on the first set of data, the running time of single machine and parallel computation are counted in nine sets of comparative experiments with a 128 × 256 pixels block size, four CPU cores, combined three different step sizes and three different measured Doppler centroid frequency estimation algorithms, respectively (Fig. 4).
The parallel computing improves the computational efficiency in different algorithms with different step sizes (Fig. 4). According to the detailed percentage of efficiency improvement (parallel time/single machine time × 100%) in Table III, the best effect is the SDE algorithm with a 2 × 4 pixels step size, which is increased by 68.34%, and the worst effect is the CPA algorithm with a 8 × 16 pixels step size. From the comparison of the means of the different step sizes, it can be seen that the smaller the step size or the longer the single line time is, the better the parallel effect. From the comparison of the means of  the different algorithms, it can be seen that the computational efficiency is improved in the following order SDE algorithm > CDE algorithm > CPA algorithm. The computational efficiency of parallel computing in the nine sets of experiments increased by an average of 43.55%.

D. Bias Correction Deviation Analysis
The analysis of bias correction deviation needs to be verified using a large-area homogeneous backscatter image relative to the geostationary because the Doppler centroid anomaly of such images should be zero, but for sea ice, such images are not easy to obtain, so a set of rainforest data (set 6) are selected to analyze the bias correction deviation by referring to the experiment conducted by Kraemer et al. [21] and [22].
In Fig. 5(a), two biases that need to be corrected are shown: the bias caused by scalloping in the azimuth direction and bias caused by the electronic mispointing of the antenna in the range direction. Assuming that the two biases are independent of each other, they can be corrected separately. In this section, the CDE algorithm applied to the set 6 data is taken as an example. The block size is 128 × 256 pixels and the step size is 4 × 8 pixels. Fig. 5(b) is a low-pass filter constructed according to Section III-D, where the blue line is the normalized magnitude of the center point of the 2-D spectrum distributed along the azimuth direction, and the red line is the constructed low-pass filter that retains effective information. The low-pass filter filters most of the high-frequency signals, effectively removing the scalloping [ Fig. 5(c)]. Fig. 5(d) displays the range mean before and after scalloping. The jagged structure of the range mean along the azimuth direction is effectively removed after correction, and the mean distribution is smoother.
The antenna electronic mispointing bias can be fitted in the range direction as described in 3.4. Fig. 6(a) shows the bias estimation of the rainforest data, and its change in the range direction is 14-16 Hz. Fig. 6(b) shows the bias correction result, that is, the final geophysical Doppler shift. By comparing the azimuth mean [ Fig. 6(c)], it can be seen that the bias caused by the antenna electronic mispointing is effectively removed, and the Doppler centroid frequency is close to zero.
For the three different measured Doppler centroid frequency estimation algorithms, the effect of bias corrections is also different. Table IV shows the geophysical Doppler frequency shift of the rainforest data obtained using the three algorithms. In removing the null boundary, the number of samples (N), mean (Mean), median (Median), standard deviation (Std), mean absolute error (MAE), and RMSE of 50-6100 samples are taken, the unit is Hz, the velocity corresponding to the RMSE at an incident angle of 39°is v, and the unit is meters/second. As shown in Table IV, the means and medians of the geophysical Doppler shift corrected using the three algorithms are close to 0. The CDE algorithm has the smallest absolute MEA and RMSE, but the overall RMSE is approximately 3 Hz, corresponding to a v of 0.07-0.09 m/s.

E. Analysis of the Radial Surface Velocity Results
Based on the data in Sets 1-5, the three Doppler centroid estimation algorithms are used to measure the sea ice drift RVL and the results are analyzed and compared.
As described in Section III, the three algorithms combined with parallel computing are used to estimate the Doppler centroid frequency. In this process, the block size is 128 × 256 pixels and the step size is 4 × 8 pixels, so that the resolution and computational efficiency can be guaranteed. Then, the bias corrections are performed to obtain the geophysical Doppler shift. After that, the image calculated using the sliding block is restored to the size before the Doppler centroid estimation (edge removal), through bilinear interpolation to achieve pixel-by-pixel incident angle and latitude and longitude matching. Finally, the geophysical Doppler shift is converted to the sea ice RVL using (3).
The results are compared below, where positive and negative velocities represent the movement of sea ice toward and away from the radar, respectively. The S1A as a right looking radar, the S1A antenna points to the right side of the flight direction. If we do not consider the heading relative to the south direction of the offset, for the descending data, toward the radar movement is moving from west to east, away from is moving from east to west, and ascending data is the opposite. In Table V, the comparison between different Doppler centroid estimation algorithms is presented by the computational time (Time), number of samples (N), mean, median, and standard deviation (Std). A comparison among the RVL obtained using the Doppler centroid estimation algorithms and the NSIDC product and the IABP buoy data is   presented using the MAE and RMSE, where the time unit is seconds and the speed unit is meters/second. Fig. 7 shows the first set of data acquired on February 2, 2021, in the Bering Strait. In the intensity image [ Fig. 7(a)], some land can be seen in the lower right corner. In the corresponding RVL in Fig. 7(b) and (c), some bright red and dark blue speckles can be seen along the land coast, but no such speckles exist in Fig. 7(d). This is because when there are ships and man-made structures in the study area, abnormal speckles will be generated in an SAR image. The CPA algorithm and CDE algorithm retain these abnormal speckles, showing bright red and dark blue spots in the image. The SDE algorithm only uses the sign of the signal to calculate the correlation function, and does not assign a large weight to bright targets, so it effectively removes these abnormal speckles. The overall image looks smoother, and the Std also decreases (Table V). However, the other hand, CPA algorithm and CDE algorithm can retain more detailed information, and the measurement effect is clearer and more  delicate. As seen in Fig. 7(e), the RVL means of the sea ice drift retrieved using the three algorithms are all positive; that is, the drift of the sea ice in this scene is from west to east in the range direction. Fig. 8 shows the second set of data obtained in the Bering Strait on December 2,2021. Two longer ice cracks are visible in the intensity image [ Fig. 8(a)], and their traces are faintly visible in Fig. 8(b)-(d). This is because seawater has a higher surface velocity than sea ice and therefore appears redder or bluer than the sea ice area in the RVL images. In Fig. 8(b)-(d) at lines 2000-4000 and 8000-10 000, the sea ice velocity is obviously different from that in other regions. In addition to the influence of the scalloping correction deviation, the RVL definition should also be considered. The RVL is only the horizontal component of the real sea ice drift velocity in the range direction. If the drift direction of the ice is vertical to the antenna direction, it is difficult to observe using the Doppler centroid estimation algorithm. Therefore, the velocity values very close to zero in these two regions do not mean that the sea ice does not actually drift. Fig. 8(e) shows that the mean RVL of the sea ice drift  retrieved using the three algorithms is mostly negative; that is, the drift of sea ice in this scene is from east to west in the range direction. Fig. 9 shows the results of the third set of data acquired on 16 December 2021 along the Greenland coast in Baffin Bay. In the upper left part of Fig. 9(b) and (c), there are some bright red or dark blue speckles. In Fig. 9(d), the speckles are effectively removed, which further illustrates the effectiveness of the SDE algorithm for abnormal speckle removal. In the intensity image [ Fig. 9(a)], there is a very small part of land area in the upper left corner, which is more than the sea ice cover on the surface of the coastal water. In the velocity image [ Fig. 9(b)-(d)], this part of the velocity value is significantly higher than other sea ice cover areas in the image. In the ice crack on the right side at lines 1000-2500 and the sea water region in the middle at lines 4000-5000, the velocity value in the RVL image is greater than that in the sea ice cover region, although the direction is not the same. In addition, it can be seen in Fig. 9(e) that the mean value of the RVL range direction of the sea ice drift measured using the three algorithms is mostly negative, that is, the sea ice drift in the scene can be regarded as moving from east to west in the range. Fig. 10 shows the relevant content of the fourth set of data obtained on December 27, 2021. This area is located on the northwest side of the Fram Strait, close to Greenland. Fig. 10(a) shows that there is a large amount of floating ice in this area. Except for complete floating ice on the right side at lines 6000-8000, there are no other obvious texture features. The floating ice can be clearly observed in Fig. 10(b)-(d). Combined with Fig. 10(e), on the whole, the direction of the sea ice drift velocity in this range is from west to east at lines 1-10 000, and from west to the east at lines 10 000-12 000.
In comparison with the NSIDC product, the temporal resolution of the Polar Pathfinder Daily 25 km EASE-Grid Sea  Ice Motion Vectors, Version 4 dataset is one day, the spatial resolution is 25 × 25 km, and the spatial resolution is approximately 360 × 360 m according to x, y coordinate interpolation. The NSIDE product is not a raster dataset, but a point-vector dataset, with each point containing its x and y coordinates, its velocity components u and v along x and y, and its latitude and longitude coordinates. The RVL results also contain latitude and longitude coordinates. The two coordinates are unified into the same coordinate system, and the velocity values of the corresponding coordinate points between them are compared. Table V shows that the sea ice velocity provided by the NSIDC product is less than the sea ice drift RVL measured by the Doppler centroid estimation algorithms. The velocity difference in Sets 1, 3, and 4 is nearly 10 times, and in Sets 1 and 3, the velocity direction is opposite. This is mainly because the NSIDC product provides the average daily velocity, while the RVL is a near-instantaneous velocity. Sea ice is a feature with strong velocity variability. The longer the time interval is, the more different the calculated sea ice average velocity is from the real value. In addition, the original resolution of the NSIDC product is 25 km × 25 km, and the velocity value used for verification is obtained through interpolation. Therefore, the average velocity of the sea ice drift provided by the NSIDC product is not enough to verify the near-instantaneous sea ice RVL. The RMSE of 0.3-0.5 m/s in Sets 1, 3, and 5 also illustrates this problem. Set 2 is the closest velocity, and the RMSE still reaches 0.2 m/s. However, on the whole, from the comparison of the Std and RMSE of the different algorithms on each set of data, it can be seen that the CDE algorithm has the highest accuracy and the SDE algorithm has better image smoothness.
To evaluate the accuracy of the sea ice RVL more comprehensively, this article also adds a fifth set of experiments (Fig. 11). The experimental scene is located in the eastern sea area of the Spitsbergen Islands. The RVL estimated using the three algorithms is compared with the measured data collected by buoys. In 10 000 lines of land coast, there are obvious dark blue areas on the left side, and obvious bright red areas on the middle and right side [ Fig. 11(b)-(d)]. These two areas are located along the land coast and are mainly composed of seawater and crushed ice, so the velocity values are relatively large.
The buoy is marked with yellow dots in the figure. The average velocity every 15 min is approximately 0.1237 m/s, and the direction is from east to west. The radial velocity value of the corresponding position is also positive (yellow dots), and the direction is from east to west. The overall mean value of the RVL calculated using the three algorithms is distributed at approximately 0.19-0.21 m/s, and the average velocity of the corresponding buoy position is approximately 0.18 m/s, which is greater than the buoy velocity. However, because the velocity calculated using the buoy is the mean within 15 min, while the RVL is a near-instantaneous velocity and only represents the upward component of the real velocity in the range direction, it is reasonable to have a certain gap between the two. Limited by a lack of measured data, there is only one buoy point for validation data. Although the validation data are not sufficient, it is still consistent with the previous conclusion. The MAE of the sea ice RVL obtained using the CDE algorithm is the smallest, which is 0.0602 m/s, and the RMSE is also the smallest, which is 0.0630 m/s, so the accuracy is the highest ( Table V).
The results of the sea ice drift RVL for the five sets of data are comprehensively analyzed to compare the three Doppler centroid estimation algorithms. The mean and median values of the RVL obtained using the three algorithms are very close, and the velocity image also shows that the sea ice movement trend is consistent. The main difference among the three algorithms is that the CDE algorithm is more efficient than the other two algorithms. At the same time, compared with the NSIDC product and the IABP buoy data, the MAE and RMSE are the smallest. It is an algorithm with high efficiency and high precision. However, in data with SAR image speckle anomalies caused by human activities, the SDE algorithm can effectively remove abnormal speckles and obtain smoother velocity images. Therefore, the Std is the smallest and has better adaptability. The CPA algorithm has no obvious advantage over the other two algorithms.

V. CONCLUSION
At present, in sea ice drift measurement, only the Doppler centroid estimation algorithm based on the frequency domain has been studied. Therefore, this article selected three Doppler centroid estimation algorithms applied to sea surface current retrieval and combined the unique TOPS method of the S1A IW data to adjust the applicability of the algorithms. The optimal parameters of the three algorithms were explored, parallel computing was added to improve the computational efficiency of each algorithm, and then they were applied to the RVL measurement of sea ice drift. The three algorithms are the CPA algorithm in frequency domain, CDE algorithm and SDE algorithm in time domain. Based on these three algorithms, this article carried out five sets of experiments, analyzed and compared the results of each algorithm, and provided optimal algorithm adaptation suggestions.
The selection of algorithm parameters considered the balance between computational efficiency and resolution. Through experimental analysis, this article determined that the optimal algorithm parameters for sea ice applications are a block size of 128 × 256 pixels and a step size of 4 × 8 pixels (azimuth × range). Using parallel computing for Doppler centroid estimation improves the computational efficiency by an average of 43.55%. Among the three algorithms, parallel computing was the best for the SDE algorithm, with an average increase of 59.39%. The deviation of scalloping correction and antenna electronic mispointing bias correction were verified using rainforest data. The deviation was approximately 3 Hz, and the corresponding velocity at a 39°incident angle was 0.07-0.09 m/s. By analyzing and comparing the results of each algorithm in the five experiments, this article found that the results of the three algorithms are relatively close. Among them, the CDE algorithm has the advantages of high computational efficiency and high precision, and is the best Doppler centroid estimation algorithm among the three. However, for data with abnormal speckles in SAR images caused by human activities, the SDE algorithm can effectively remove abnormal speckles and directly obtain smooth Doppler centroid estimation images, which has better adaptability. The parameter selection scheme, parallel computing method, and Doppler centroid estimation algorithm applicability suggestions for different data types discussed in this article can provide a reference for subsequent research on radial surface velocity measurements of sea ice drift.