An investigation of inversion method to measure the radial velocity of Kuroshio from Sentinel-1 SAR data

ABSTRACT

Studies on the retrieval of ocean surface current by Synthetic Aperture Radar (SAR) mainly concentrate on coastal regions with land coverage, allowing for continuous land to estimate non-geophysical Doppler shifts. As a strong western boundary in the Western Pacific, the Kuroshio flows through mostly islands or open sea, and lacks land coverage for range bias correction. We propose a non-geophysical Doppler shift correction algorithm suitable for Kuroshio observation. Three solutions matching different land coverage were adopted in the range direction, and the scalloping was removed by using average filtering in the azimuth direction. Comparison between SAR radial velocity and Global Surface Lagrangian Drifter (GLD) reveals that the correlation coefficient (R²) increase to 0.686, and root mean square error (RMSE) of 0.215 m/s after non-geophysical corrections, while the wind-wave bias corrections make R² increase to 0.806 and RMSE decrease to 0.13 m/s. There is also a good consistency between SAR inversion and Regional Ocean Circulation Model (ROMS). Sensitivity analysis illustrates that ascending pass is more suitable for monitoring Kuroshio, and high incidence angles are less affected by ocean conditions and non-geophysical factors. The algorithm can be used in regions with discontinuous or absent land.

Highlights

• Three schemes for various land coverage in the range direction, including AP, IAP, and NIAP.

• Using average filter to remove the scalloping signal in the azimuth direction.

• SAR currents compared with GLD data, $R^2$ of 0.806 and $RMSE$ of 0.13 m/s.

• The method can be used in both coastal and island/offshore regions.

KEYWORDS:

• SAR
• Doppler shift
• non-geophysical
• ocean current
• Kuroshio

1. Introduction

Ocean surface circulation is a large-scale and stable flow of seawater that affects the transport of salt and heat, and plays a crucial role in the formation and evolution of climate and weather (Gens 2008). Ocean currents observations are essential for a variety of applications, including marine debris tracing, ship navigation, fish industry management (Röhrs et al. 2021), and high-potential tidal energy development (Ferreira, Estefen, and Romeiser 2016).

Direct ocean surface current measurements include in-situ observation methods, such as drifting buoys, shipboard measurements, and moored instruments (Rascle et al. 2020). These methods have high accuracy, but it is difficult to capture the spatiotemporal variability of ocean currents because of their high cost and limited spatial coverage (Wang et al. 2021). Therefore, these are only used as validation data for large-scale ocean current observations. In recent years, microwave remote sensing techniques, such as high-frequency radar (HFR) (Paduan and Washburn 2013), radar altimetry (Yu et al. 2018), and synthetic aperture radar (SAR) (Elyouncha, Eriksson, and Johnsen 2022) have been widely applied to map ocean currents at multiple spatial scales. HFR obtains the ocean current with a wide coverage range within 400 km offshore, providing temporal and spatial resolutions of 1 hour and 3–6 km, respectively. These have been applied in oil spills (Bellomo et al. 2015), search-and-rescue operations (Shen et al. 2019) and other maritime activities. The spatial and temporal resolutions of radar altimetry are approximately 25 km and 10 days, respectively (Johannessen et al. 2014). These methods are generally unsuitable for monitoring sub-mesoscale ocean processes within 50 km of the coastline. However, understanding the variability of ocean currents in sub-mesoscale or coastal areas is important for offshore activities, including coastal port construction, trade vessel route planning, and marine disaster monitoring (Melet et al. 2020). SAR as an active imaging microwave remote sensing is capable of capturing surface ocean processes such as winds, waves, and currents in all weather conditions (Kravtsov et al. 1998; Moiseev et al. 2020).

Along-Track Interferometry (ATI) and Doppler Centroid Anomaly (DCA) are two major methods for measuring the ocean surface current from SAR data. The ATI method uses two along-track SAR data to inverse current with a high spatial resolution (~100 m) (Rashid and Gierull 2021; Yoshida, Ouchi, and Yang 2021), but it requires high imaging quality and cost (Wang et al. 2022). The DCA method is suitable for monitoring sea surface dynamics using SAR data with a single antenna. Currently, most SAR images are acquired using a single antenna, such as Sentinel-1 (Johnsen et al. 2016), ENVISAT-ASAR (Hansen et al. 2011) and GF-3 (Sun et al. 2022). Therefore, the amount of available data for DCA is much greater than ATI, making it suitable for the long-term ocean currents monitoring. Chapron et al. (2005) analyzed the feasibility of DCA method to retrieve ocean currents in terms of Doppler frequency shift theory and summarized the factors that cause Doppler shift bias (Chapron, Collard, and Ardhuin 2005). Johannessen et al. (2008) and Mouche et al. (2012) subsequently constructed the RIM model (Kudryavtsev et al. 2005), DopRIM model (Johannessen et al. 2008) and CDOP model (Mouche et al. 2012). These studies explained the effects of wind, waves, and current interactions on the Doppler shift. The DCA method has been applied in regions such as Agulhas (Krug et al. 2010), Nordic (Hansen et al. 2012), and the Yangtze River (Yang and Yijun 2023). These results had been compared with drifting buoys and HFR radar (Martin et al. 2022), and its advantages and accuracy in weak and strong ocean currents monitoring have been verified.

Accurate separation of geophysical and non-geophysical Doppler shifts is essential for DCA method. Non-geophysical Doppler shifts consist of satellite geometry/attitude, scalloping, antenna miss-pointing, and all unknown errors (Elyouncha et al. 2019). Doppler shifts over the land are usually used as a reference for the estimation and elimination of antenna miss-pointing, and all unknown errors. Studies on the ocean current retrieval using Sentinel-1 C-band SAR data focus on coastal regions, which can provide continuous reference data in the range direction (Elyouncha et al. 2019; Martin et al. 2022; Moiseev et al. 2020). However, most of these take less account of the land distribution, and a systematic method for eliminating non-geophysical components in the presence of discontinuous or non-existent land cover remains lacking. When the land in SAR data is discontinuous or absent, it is easy to cause discontinuity or overestimation. The Kuroshio, an ocean current in the western Pacific, exhibits high temperature, salinity, and velocity, significantly influencing weather and climate changes in China’s coastal areas, East Asia, and even the world. It is worth noting that the Kuroshio, unlike the coastal currents, flows through mostly islands or open seas. It is crucial to effectively utilize the limited land reference data for estimation and correction of non-geophysical component. Consequently, referring to traditional non-geophysical Doppler correction methods, we propose a method suitable for Kuroshio monitoring.

In this study, Kuroshio inversion was accomplished by using the Sentinel-1 Level 2 Ocean (OCN) product, and the non-geophysical terms separately for each swath were corrected in the range and azimuth directions. We used three corresponding schemes by considering three different land distributions in the range direction, while the scalloping can be removed by mean filtering in the azimuth direction, and the Global Surface Lagrangian Drifter (GLD) data and Regional Ocean Circulation Model (ROMS) numerical simulation currents were used to evaluate the inversion results.

2. Datasets

2.1. Sentinel-1 SAR level-2 data

The SAR data were provided by the Sentinel-1 Level 2 Ocean (OCN) data with a total swath approximately 250 km wide, including ocean wind field (OWI, spatial resolution of 1 km $\times$ 1 km) and Radial Surface Velocity (RVL, spatial resolution of 0.9 km $\times$ 0.7–0.8 km) (Biron, Van Wychen, and Vachon 2018). The IW mode consisted of three sub-swaths (IW1, IW2, and IW3) with incidence angles from near to far between 30°to-36°, 36°to-42°, and 41°to-46°, respectively. Figure 1 shows the spatial locations of the adopted nine ascending scenes (2019/09/19, 2019/09/24, and 2019/09/05) and 11 descending scenes (2019/01/11, 2018/01/21, and 2019/09/06).

Figure 1. Bathymetric map of the Kuroshio region. The solid and dashed lines represent the ascending and descending SAR data passes, respectively, the red circles are GLD buoy data, and the download link for bathymetry data: https://topex.ucsd.edu/pub/srtm15_plus/.

2.2. Validation datasets

The validation datasets included Global Surface Lagrangian Drifter (GLD) data, Global Ocean Multi-observation products (MOB), and Regional Ocean Circulation Model (ROMS) numerical simulation currents. GLD data were provided by the Drifter Data Assembly Center at NOAA (Atlantic Oceanographic and Meteorological Laboratory, https://www.aoml.noaa.gov/envids/). The ocean surface current velocity is calculated from the position changes captured by buoys, and then optimally interpolated to uniform 6-hour resolution (Krug et al. 2010; Wang et al. 2021). To compare with the SAR radial current velocity, GLD data recorded during 48 h were collected and projected onto the SAR range direction. A total of 67 points in the ascending scene and 27 points in the descending scene were used to validate SAR velocities.

MOB product (http://marine.copernicus.eu) belongs to global total ocean current reanalysis data with a 1 h frequency and a one-fourth deg regular grid, and the product ID is MULTIOBS_GLO_PHY_MYNRT_015_003. The total current fields consist of geostrophic and Ekman currents at the surface and 15 m depth, obtained by combining CMEMS satellite-observed geostrophic surface currents with Ekman currents simulated using ERA5 wind stress. Compared with MOB product, model currents have higher spatial resolution. The numerical simulation of the ocean surface current was based on a high-resolution 3D model established by the ROMS model. It employed terrain-following coordinates in the vertical direction with 15 layers and a grid spacing of 1.5 km. The model adopts the Bulk Formula method for atmospheric forcing, with upper boundary condition derived from the National Centers for Environmental Prediction (NCEP) data at 6-hour intervals. Vertical mixing and tidal forcing were implemented using universal length-scale turbulence mixing and the TPXO9 dataset. The initial conditions and open boundary conditions for the model have been derived from HYCOM hydrographic dataset (spatial resolution of 0.08° × 0.08°) (Jang et al. 2021). Finally, the temporal and spatial resolutions of the output data are 0.5 h and 1.5 km, respectively.

These datasets are potentially different due to differences in sources. GLD data have the highest accuracy, but the coverage and quantity are limited. MOB product can effectively show the direction of strong ocean currents, but the spatial resolution is low, suitable only for displaying the direction of the Kuroshio central current. Compared with reanalysis data, ROMS model currents enhance the spatiotemporal resolution of simulated flows by employing higher-resolution grids and time steps, serving as an effective supplement to GLD and MOB data.

3. Methodology

3.1. Doppler shift caused by currents

For Sentinel-1 OCN data, the total observed Doppler shift ($f_{dc}^{obs}$) includes geophysical and non-geophysical components (Equation (1)(1) $f_{dc}^{obs}=f_{dc}^{phys}+f_{dc}^{non-phys}$(1) ).

(1) $f_{dc}^{obs}=f_{dc}^{phys}+f_{dc}^{non-phys}$(1)

(2) $f_{dc}^{phys}=f_{dc}^{ww}+f_{dc}^{osc}$(2)

(3) $f_{dc}^{non-phys}=f_{dc}^{geo}+f_{dc}^{scal}+f_{dc}^{elec}+\mathrm{\Delta }{f}_{dc}$(3)

where $f_{dc}^{phys}$ is the geophysical term induced by wind-wave $f_{dc}^{ww}$ and current $f_{dc}^{osc}$ (Equation (2)(2) $f_{dc}^{phys}=f_{dc}^{ww}+f_{dc}^{osc}$(2) ) (Moiseev, Johannessen, and Johnsen 2022). From Equation (3)(3) $f_{dc}^{non-phys}=f_{dc}^{geo}+f_{dc}^{scal}+f_{dc}^{elec}+\mathrm{\Delta }{f}_{dc}$(3) , $f_{dc}^{non-phys}$ is the non-geophysical term, groups the geometrical term $f_{dc}^{geo}$, the antenna miss-pointing term $f_{dc}^{elec}$, the scalloping ${\mathit{f}}_{dc}^{scal},$and unknown error term $\mathrm{\Delta }{f}_{dc}$. Among these, $f_{dc}^{geo}$ is generated by the satellite velocity relative to the rotating earth, and $f_{dc}^{elec}$ is due to the deviation between the actual antenna pattern and the theoretical antenna pattern, which can be regarded as a function of the incidence angle. $f_{dc}^{scal}$ occurs due to the antenna gain variation, and $\mathrm{\Delta }{f}_{dc}$ is due to the sensor temperature compensation and inaccurate estimation of the non-geophysical Doppler shifts.

The initial calibrations for $f_{dc}^{geo}$ and $f_{dc}^{elec}$ are based on the variables (“rvlDcGeo” and “rvlDcMiss”) provided by the Sentinel-1 Level 2 OCN product. However, these two computed terms are imprecise (Johnsen et al. 2016). There are noticeable residual non-geophysical biases, as scalloping and the land signals are not equal to 0 Hz (Danilo and Bella 2018). The residual errors mentioned above, together with $f_{dc}^{scal}$ and $\mathrm{\Delta }{f}_{dc}$ are defined as $f_{dc}$ in this study. As these errors, $f_{dc}$ appears linear variation with incidence angle in the radial or range direction and jumps over sub-swaths. This is related to the antenna miss-pointing, and this “range bias” usually requires subtracting the average Doppler shift over land. In addition, $f_{dc}$ presents a scallop effect with significant periodic changes (±10 Hz) along the azimuth direction. This is an artifact that hinders the direct analysis of the spatial changes of surface current velocity (Wang et al. 2023). Therefore, $f_{dc}$ can be decomposed into two radar directions for correction: the range and the azimuth direction (see 3.2 for more details). The geophysical Doppler $f_{dc}^{phys}$ is obtained after subtracting $f_{dc}$. In order to retrieve ocean surface current $f_{dc}^{osc}$ from $f_{dc}^{phys}$, the empirical geophysical model function (GMF) called CDOP is used to estimate and remove the wind-wave bias $f_{dc}^{ww}$ (see Equation (4)(4) $f_{dc}^{osc}=f_{dc}^{obs}-f_{\mathrm{d}\mathrm{c}}-f_{dc}^{ww}=f_{dc}^{obs}-f_{\mathrm{d}\mathrm{c}}-f_{dc}^{cdop}\left(u_{10},\varphi ,\theta ,p\right)$(4) ) (Elyouncha, Eriksson, and Johnsen 2022; Moiseev et al. 2020). The input parameters of the CDOP model include wind speed ($u_{10}$) and direction ($\varphi$) at 10 m height retrieved from SAR images provided in Sentinel-1 Level 2 OWI products (“owiWindSpeed” and “owiWindDirection”), incidence angle ($\theta$) and polarization ($p$) at C-band SAR. The Doppler shift contributed by the ocean current can be expressed as:

(4) $f_{dc}^{osc}=f_{dc}^{obs}-f_{\mathrm{d}\mathrm{c}}-f_{dc}^{ww}=f_{dc}^{obs}-f_{\mathrm{d}\mathrm{c}}-f_{dc}^{cdop}\left(u_{10},\varphi ,\theta ,p\right)$(4)

The SAR radial surface velocity can be derived from the estimated $f_{dc}^{osc}$ as Equation (5)(5) $U_{osc}=-\frac{\pi \cdot f_{dc}^{osc}}{k_e\cdot sin\theta }$(5) ,

(5) $U_{osc}=-\frac{\pi \cdot f_{dc}^{osc}}{k_e\cdot sin\theta }$(5)

where $k_e$ is the electromagnetic wave number for Sentinel-1, which is equal to 113.28 m⁻¹, and $\theta$ is the incidence angle. The overall workflow (Figure 2) of radial velocity inversion includes range and azimuth calibrations.

Figure 2. The flowchart of the data processing.

3.2. Non-geophysical Doppler shift correction algorithm

3.2.1. Correction in the range direction

The biased Doppler shift $f_{\mathrm{d}\mathrm{c}}$ fails to correctly represent the non-geophysical parameters variations along the range direction, especially as observed over land. Based on the assumption that the Doppler shift over land should be 0 Hz, the residual deviation can be approximated by the averaged Doppler shift in data over land (Wang et al. 2022). However, there are three types of land coverage in the range direction according to the “rvlLandCoverage” variable from Sentinel-1 OCN data: (1) continuous land coverage, (2) non-continuous land coverage, and (3) no land coverage. We proposed three corresponding schemes.

Scheme i, the method applies to continuous land in the range direction, as seen the IW3 in Figure 3(a1) of the SAR image from 19 September 2019. An algorithm combining average filtering and polynomial fitting (AP) was proposed. Specifically, the mean range profile $f_{dcRm}$ was first obtained by averaging the Doppler shifts over land along the azimuth direction (Equation (6)(6) $\begin{array}{c}{f}_{dcRm}\left(j\right)=\frac{1}{n}\sum _{j=1}^n{f}_{dc}\left(i,j\right)\end{array}$(6) ). From Equation (7)(7) $\begin{array}{c}{f}_{dcR}=polyfit\left(j,filter\left(f_{dcRm}\right),3\right)\end{array}$(7) , the $f_{dcRm}$ is smoothed by an average filter to reduce the interference of abrupt signals, followed by polynomial fitting to obtain the range profile $f_{dcR}$. The calibration of range bias is achieved by subtracting $f_{dcR}$ from $f_{dc}$.

(6) $\begin{array}{c}{f}_{dcRm}\left(j\right)=\frac{1}{n}\sum _{j=1}^n{f}_{dc}\left(i,j\right)\end{array}$(6)

(7) $\begin{array}{c}{f}_{dcR}=polyfit\left(j,filter\left(f_{dcRm}\right),3\right)\end{array}$(7)

Figure 3. Three correction schemes in the range direction.

The “Note” in the manuscript is missing from the lower right corner.

where $i$ and $j$ denote the azimuth and range indices of Doppler shift grid, which are equal to the rows and columns in 2D matrix, respectively. $n$ is the component in the range direction in set $\left\{1,2,\cdots ,j\right\}$. $filter$ and $polyfit$ are 1D digital average filter and polynomial fit functions. The results of $f_{dcRm}$ and $f_{dcR}$ correspond to the blue and red lines in Figure 3(a2), respectively.

Scheme ii, a method for the discontinuous land in the range direction, such as IW3 in Figure 3(a2) from 21 January 2018. Using Equation (6)(6) $\begin{array}{c}{f}_{dcRm}\left(j\right)=\frac{1}{n}\sum _{j=1}^n{f}_{dc}\left(i,j\right)\end{array}$(6) fails to produce a continuous $f_{dcRm}$. To avoid the strong discontinuities in the calibration results, we combine an interpolation method with Scheme i (called IAP). Specifically, using a moving window to interpolate the missing values in $f_{dcRm}$ with the median value in the neighborhood. As for the SAR image in Figure 3(a2), the window size of 19, and the interpolated profile is denoted by a dotted line in Figure 3(b2). Then, using Equation (7)(7) $\begin{array}{c}{f}_{dcR}=polyfit\left(j,filter\left(f_{dcRm}\right),3\right)\end{array}$(7) to compute $f_{dcR}$, and subtract it from $f_{dc}$ to perform range bias calibration.

Scheme iii, a method named NIAP for the no land in SAR image or reference image (adjacent orbit and same period). As seen the IW1 in Figure 3(a3) of from 6 September 2019, neither scheme i or ii is applicable. This limits the application of SAR data to monitor ocean currents without land coverage. If the Normalized Radar Cross Section (NRCS) variations of the SAR data are not significant (i.e. non-complex sea conditions), we propose exploiting the Doppler shift over the ocean to simulate non-geophysical Doppler shift variations, following the linear variation in the range direction (Wang et al. 2022; Yang and Yijun 2023). However, it is worth noting that the Doppler shifts over the land are often smaller than in the ocean due to the sea surface motion. Therefore, the empirical conditions adopted by Danilo and Bella (2018) was introduced to calculate $f_{dcRm}$. According to the “rvlDcObsStd” variable in Sentinel-1 Level 2 OCN products, the outliers or high-frequency values with a standard deviation exceeding 4 Hz were first removed. Then, $f_{dcRm}$ is obtained by averaging along the azimuthal direction, and we also set the minimum sample for computing the profile to 10 pixels to reduce the interference from the abrupt signal. Additionally, it is necessary to apply interpolation consistent with Scheme ii to ensure $f_{dcRm}$ is continuous. Finally, $f_{dcR}$ is obtained by Equation (7)(7) $\begin{array}{c}{f}_{dcR}=polyfit\left(j,filter\left(f_{dcRm}\right),3\right)\end{array}$(7) , and subtracted from $f_{\mathrm{d}\mathrm{c}}$ to calibrate the range bias.

The above correction method in the range direction can be expressed by the following pseudocode:

Table

Range Bias Calibration Algorithm

Input: biased Doppler shit $f_{\mathrm{d}\mathrm{c}}$, landmask $Valu{e}_{Land}$, standard deviation $\delta$ Output: mean range profile $f_{dcRm}$, range profile $f_{dcR}$ $ifValu{e}_{Land}=100\mathrm{\%}continuouslandintherangedirection$ $f_{dcRm}=mean\left(f_{\mathrm{d}\mathrm{c}}\right)$ $f_{dcR}=ployfit\left(size\left(f_{\mathrm{d}\mathrm{c}}\right),2\right),filter\left(f_{dcRm}\right),3)$ $elseifValu{e}_{Land}=100\mathrm{\%}continuouslandintherangedirection$ $f_{dcRm}=fillmissing\left(mean\left(f_{\mathrm{d}\mathrm{c}}\right){,}^{\prime }media{n}^{\prime },win\right)$ $f_{dcR}=ployfit\left(size\left(f_{\mathrm{d}\mathrm{c}}\right),2\right),filter\left(f_{dcRm}\right),3)$ $elseValu{e}_{Land}=0\mathrm{\%}nolandintherangedirection$ $f_{dcRm}=fillmissing\left(mean\left(f_{\mathrm{d}\mathrm{c}}\left(\delta \le 4Hz\right)\right){,}^{\prime }media{n}^{\prime },win\right)$ $f_{dcR}=ployfit\left(size\left(f_{\mathrm{d}\mathrm{c}}\right),2\right),filter\left(f_{dcRm}\right),3)$ Note: $Valu{e}_{Land}$ is the value corresponding to the variable “rvlLandCoverage” of Sentinel-1 Level 2 OCN, 100 and 0 represent land and sea, respectively. $mean\left(f_{\mathrm{d}\mathrm{c}}\right)$ is the mean value from each column, and the pixels of each column exceed 10. $size\left(f_{\mathrm{d}\mathrm{c}}\right),2)$ and $fillmissing$ are the number of columns of the image and interpolation, respectively. As seen the histogram in Figure 3(c1,c2), Doppler shifts over land satisfies the assumption of being equal to 0 HZ.

3.2.2. Correction in the azimuth direction

There still exists $\pm$10 Hz periodic variation along the azimuth direction in the SAR image after range bias calibration. To remove the scalloping, we define the Doppler shift after range calibration as $f_{dc-dcR}$. We aim to identify the burst period of the scalloping by calculating $f_{dcAm}$, $f_{dcAs}$ and $f_{dcA}$. Specifically, $f_{dcAm}$ is the mean azimuth profile obtained by averaging along the range direction in each sub-swath data, and an average filter smoothing is applied to $f_{dcAm}$ to get $f_{dcAs}$. The period variation signal $f_{dcA}$ then is obtained by subtracting $f_{dcAs}$ from $f_{dcAm}$. In Figure 4,we select one of the SAR images in Figure 3(a1) to show the azimuth calibration process. Figure 4(a) shows IW1 of the selected data, and the blue, black, and red lines correspond to $f_{dcAm}$, $f_{dcAs}$ and $f_{dcA}$, respectively. The $f_{dcA}$ in IW2 and IW3 are also obtained through this method (Figure 4(b)). Finally, $f_{\mathrm{d}\mathrm{c}\mathrm{A}}$ is subtracted from $f_{\mathrm{d}\mathrm{c}-\mathrm{d}\mathrm{c}\mathrm{R}}$ to complete the azimuth calibration.

Figure 4. Doppler shifts in the azimuth direction, with is a $\pm$10 Hz periodic variation. Figure (a) shows the identification process of the scalloping (only IW1 swath), while figure (b) shows the method applied to the full image.

4. Results

4.1. Kuroshio retrieval

As shown in Figure 5, the $U_{osc}$ is positive when the ocean surface current moves away from the radar sensor, while negative velocity is toward the radar sensor (Van Wychen et al. 2018). The Appendix contains the results of current inversions of SAR data for other time phases. The velocity of the central Kuroshio is between 0.5 and 1 m/s and meanders along the Okinawa Trough from southwest to northeast.

Figure 5. Result of current retrieval. Figures (a) and (b) show the radial velocities on September 24, 2019, and January 21, 2018, respectively. The Kuroshio central current is the blue line, and the Kuroshio counter Current is the red line. The arrows are derived from the MOB products.

$U_{osc}>0$ corresponds to red in Figure 5(a), indicating the Kuroshio central current at (28°N-30°N) meanders northeastward along the Ryukyu Trench with a width of 40–60 km. At 28.5°–30.5° N, 127.5°–130° E, the velocity in the Kuroshio central current region is significantly higher than in the neighborhoods, forming a straight and fast-flowing current zone that bends and flows toward the adjacent sea regions, rather than a single flow direction. $U_{osc}<0$ corresponds to red in Figure 5(b) and clearly shows the characteristics of the Kuroshio central current, but the width is narrower, between 20 and 30 Km. Additionally, the Kuroshio Counter Current is more significant in Figure 5(b), corresponding to the blue region with $U_{osc}>0$, which is located at 26.5°–27.8°N, 127.5°E-128.5°E, and the velocities are between 0.4 and 0.6 m/s. Both SAR radial velocity and MOB products observe the Kuroshio central current (blue line) and the Kuroshio Counter Current (red line).

4.2. Comparison with GLD data

To compare with different ocean current data, the GLD data and ROMS ocean model were projected onto the SAR radial direction (Rio, Mulet, and Picot 2014). The statistical analysis is found by comparing the calibrated radial velocities from 6 Sentinel-1 scenes with 94 collocated GLD points. In Figure 6(a), $U_{osc}$ is only corrected the scalloping without range bias calibration, and the root mean square error (RMSE) and correlation coefficient (R²) are 0.736 m/s and 0.519, respectively. After correcting in the range direction, R² increased from 0.519 to 0.686, while RMSE decreased from 0.736 m/s to 0.215 m/s (Figure 6(b)). After removing the wind-wave $f_{dc}^{ww}$, R² ≈0.806 and RMSE ≈0.13 m/s (Figure 6(c)). These all passed the significance test with a P-value less than 0.001. The correlation coefficients are significant at a confidence level of 99%. Therefore, the ocean current velocities retrieved from SAR image have a reliable accuracy that can be applied to the comprehensive observation of ocean hydrodynamic parameters.

Figure 6. Comparing SAR radial surface velocity with the GLD data. Figure (a) is only the azimuth correction completed, while figure (b) shows the results of range and azimuth calibration. Figure (c) shows the results of wind-wave bias correction.

4.3. Comparison with ocean model

Figure 7(a) shows the ocean current provided by the ROMS model on 24 September 2019, with 1.5 km × 1.5 km. This is generally consistent with the SAR inversion, but there are also local differences, especially in the coastal regions (black circles). There is a coastal current from northeast to southwest at approximately −0.3 to −0.5 m/s, but this was not detected by SAR in Figure 5(a). We infer that this is related to the interference of land signals in SAR images. There is a significant change in the gradient of the land echo signal along the coastal and land boundaries. Current research usually masks a distance range (approximately 10 km) along the coast (Moiseev, Johannessen, and Johnsen 2022; Wang et al. 2021). For statistical analysis, we randomly selected 1949 points within the Kuroshio Current of each swath. Figure 7(b) shows an R² of 0.664 and RMSE of 0.12 m/s, while Figure 7(c) demonstrates that approximately 70% of the values fall within the range of −0.3 to 0.3 m/s. The visualizations of the two types of ocean current data are generally consistent, but the retrieved velocity is typically higher than the velocity provided by the ROMS model. This is related to disparity in data sources. The spatial resolution of the reanalysis data selected for the initial field of the ROMS model is low, and the grid is rough. This may result in an underestimation of the velocity in some regions.

Figure 7. Comparison with ROMS model. (a) shows the ROMS model projected onto the SAR radial direction. (b) depicts a scatter, with different colors representing matching points from different swaths. (c) shows the histogram of the velocities difference (D-Value) between SAR and ROMS.

5. Discussion

5.1. Influence of proposed non-geophysical corrections

In Figure 8, we take two SAR scenes in Figure 5 on 24 September 2019, and 21 January 2018, as examples to show the changes before and after non-Doppler shift correction. From Figure 8(a1,a2), it is difficult to directly observe the features of the ocean currents due to the non-geophysical Doppler effect, such as Doppler shifts over land are not equal to 0 Hz and scalloping. The empirical strategies were proposed to reduce the impact of non-geophysical parameters on ocean current observations. The mean value of Doppler shifts over land after bias correction close to 0 Hz after range bias calibration (see Figure 8b1,b2), and this satisfies the “zero motion” assumption. However, the significant scallop signals can still be observed in the azimuth direction. From Figure 8(c1,c2,c3), it is found that the correction method in Section 5.2.1 is suitable for scalloping calibration and improving the quality of Sentinel-1 OCN data. We subsequently proceed to assess the effectiveness of the proposed calibration scheme.

Figure 8. Non-geophysical Doppler correction. (a1) and (a2) show the biased Doppler shift $f_{\mathrm{d}\mathrm{c}}$. (b1) and (b2) show the results after correction in the range direction. (c1) and (c2) show the effective geophysical Doppler shift $f_{dc}^{phys}$ after removing the scalloping.

The first significant artifact is that the Doppler shift over land is not equal to 0 Hz, which is related to inaccurate estimation of $f_{dc}^{elec}$ and $\mathrm{\Delta }{f}_{dc}$. A common solution is averaging over land to estimate the range bias. For scenes with land, Schemes i and ii in Section 3.2.2 can be selected to correct the absolute bias. In Figure 9, both ascending (Figure 9(a)) and descending (Figure 9(b)) passes exhibit similar variations in the range direction. From IW1 to IW3, the range bias shows a linear downward as the incidence angles increase and jump over swaths. The values of IW1 are higher than the other two swaths, so the impact of non-geophysical bias is more significant in the near range direction. The results indicate that the range bias can be approximated as a function of the incidence angle.

Figure 9. Non-geophysical Doppler shift variation in range direction. Figures (a) and (b) show the estimated range bias for the ascending (19 September 2019) and descending (11 January 2019), respectively. Blue line is the mean Doppler shift over land along azimuth, and red line is the error fitting value.

For scenes without land, it is possible to estimate by interpolating from reference images with the adjacent orbits and the same period. However, there are some extreme phenomena without appropriate reference data in the open sea. Therefore, Scheme iii in Section 3.2.2 attempts to use data over the ocean under non-complex sea conditions. From Figure 10(a), Area 1 and Area 2 are selected to compare different range bias correction schemes. For Area 1, we compared Scheme i and iii, and the red lines in Figure 10(b) show the range bias fitting of the two methods, with a mean deviation of about 8 Hz. The blue lines in Figure 10(b) display the range bias fitting of Scheme i and 3 in Area 2, with a mean deviation of about 6 Hz. In summary, the error between Scheme iii and Scheme i or ii is within one scalloping period. The statistical analysis in Figure 10(c) reveals the correlation of range calibration results between Scheme iii and Scheme i, as well as between Scheme iii and Scheme ii, and R² equal to 0.84 and 0.98, respectively. Therefore, in the absence of a land reference image, if the variations in NRCS are not significant, that is non-complex sea conditions such as typhoons, strong tidal currents, and eddies. Scheme iii is suitable for range error correction without land coverage.

Figure 10. Comparison of three range bias estimation schemes.

5.2. Sensitivity analysis for SAR current results

5.2.1. Influence of SAR geometry parameters

The satellite track and incidence angles are the two important geometric parameters that impact the quality of radial velocity retrieval. The Sentinel-1 track angle of ascending and descending passes are ∼350° and ∼190° (clockwise rotation from the north), respectively. In the study area, these passes are, respectively, aligned with the offshore (from land to sea) and onshore (from sea to land) wind direction, as well as correspond to the direction of Kuroshio (southwest-northeast) and Kuroshio Counter Current (northeast-southwest). When the wind-wave is propagating near the satellite track or azimuth direction, the impact on gathering Doppler signals collected in the range direction will be lessened, and this is more significant in high incidence angle. Therefore, in the ascending pass on 24 September 2019, the observed trajectory of the Kuroshio central current reveals a higher degree of overlap with GLD buoy data in the same region compared to the descending pass. We conclude that ascending passes are more suitable for monitoring and mapping ocean currents for the Kuroshio.

Both radar backscatter and Doppler shift estimations from wind/waves are affected by the incidence angle. The lower incidence angle receives stronger echo signals. This implies that low incidence angles are more sensitive to high flow velocity and rough ocean surfaces and thus have a greater impact on ocean current inversion accuracy. Based on the collected Sentinel-1 SAR images, we found that the observed Doppler shifts in IW1 are usually higher than those in IW2 and IW3, and the error of $f_{\mathrm{d}\mathrm{c}}$ was correspondingly more significant. When the land coverage is discontinuous, SAR images with land coverage at neighboring orbits and times can be used as a reference and combined with moving average and polynomial fitting for error estimation. When there is no land coverage, the Doppler shift in the range direction can be estimated by interpolating the neighboring values in the neighborhood; however, this scheme is not optimal and is only applicable to calm ocean conditions. Therefore, schemes ii and iii in this study may be prone to overestimating the error in the range direction, as observed in IW1. The incidence angle increases from 30° to 46° when compared to GLD data, and the RMSE is 0.16 m/s, 0.12 m/s, and 0.11 m/s, respectively. Furthermore, the agreement of velocity with GLD and ROMS was better in IW2 and IW3 than in IW1, with IW3 showing the best effect. In general, applying our suggested methods to remove non-geophysical components can reveal the underlying geophysical signatures and improve the quality of the data.

5.2.2. Influence of wind-wave bias on RVL

The SAR radial surface velocity (RVL, also called Doppler velocity) generated by wind/wave can be estimated and eliminated by CDOP model (Elyouncha, Eriksson, and Johnsen 2022). To distinguish between the RVL generated by wind/wave and current, we define the former as $U_{ww}$, which is related to radar parameters (incidence angle and polarization) and ocean conditions (wind speed and direction), as shown in Figure 11.

For the C-band SAR data, $U_{ww}$ in the HH polarization is higher than that in the VV polarization, making it more sensitive to signals produced by the wind wave. In Figure 11(a), $U_{ww}$ reduces as the incidence angle increases when the wind speed $u_{10}=7\hspace{0.17em}m/s$, and the relative wind direction $\varphi ={30}^{\circ }$, indicating the influence of wind and waves on the ocean current is lessened and the accuracy of radial velocity retrieval is improved. In Figure 11(b), when the incidence angle $\theta ={30}^{\circ }$ and relative wind direction $\varphi ={30}^{\circ }$, $U_{ww}$ increases with the wind speed, showing that there is a more significant effect of wind/wave on current at high wind speed. Figure 11(c) shows $U_{ww}$ variation with relative wind direction for incidence angle $\theta ={30}^{\circ }$ and wind speed $u_{10}=7\hspace{0.17em}m/s$, and $U_{ww}$ is not symmetrical to the wind direction. $U_{ww}$ is greater in upwind ($\varphi =0^{\circ }$) than in downwind ($\varphi =\pm {180}^{\circ }$), while $U_{ww}=0\hspace{0.17em}m/s$ in the crosswind. Therefore, the influence of wind/wave bias is more significant in the upwind than in the downwind and crosswind ($\varphi =\pm {90}^{\circ }$). However, assessment of the contribution of the sea state to the geophysical Doppler shift remains a challenging task (Kudryavtsev et al. 2023). In the future, we will consider the influence of wave parameters such as ocean significant wave height, peak frequency, and direction under complex ocean conditions (Yurovsky et al. 2019).

Figure 11. Variation of the Doppler velocity contributed by the ocean condition for different incidence angles, wind speeds, wind directions, and polarization parameters.

6. Conclusion

In this study, we developed a non-geophysical Doppler shifts correction algorithm. The advantage of the algorithm lies in resolving the absence of land coverage for range bias calibration, when using Sentinel-1 IW mode SAR to monitor Kuroshio. In addition, we compared and analyzed the SAR current inversion with the ocean current derived from Global Surface Lagrangian Drifter (GLD) and Regional Ocean Circulation Model (ROMS).

The proposed algorithm includes two steps: range and azimuth direction correction. In the range direction, we considered various land coverages and proposed three schemes (AP, IAP, and NIAP) to minimize residual antenna miss-pointing and unknown error. The core of the schemes includes average filter, polynomial fit, and interpolation. The 1D digital average filtering is also used to eliminate the scalloping in the azimuth direction. The non-geophysical Doppler shifts in these two directions are theoretically independent; therefore, the two steps can be exchanged when processing the data. The non-geophysical Doppler shifts correction algorithm can significantly improve the Sentinel-1 OCN data, resulting in a higher quality of Kuroshio observation and mapping.

The analysis of the SAR current inversion shows that the velocity of the central Kuroshio is between 0.5 and 1 m/s. The ascending and descending passes are suitable for the observation of Kuroshio central current and Kuroshio counter current, respectively. After non-geophysical and wind-wave Doppler calibration, the R² between SAR radial velocity and GLD range velocity are 0.686 and 0.806, with RMSE of 0.215 and 0.13, respectively. Furthermore, there was good agreement between the SAR current inversion and the ROMS ocean model. The regions where the Kuroshio flows are mainly situated in the open sea, lacking HF radar, gliders, and other observation. This limitation hinders further studies to conduct more statistical and comparative analyses. Nevertheless, the non-geophysical Doppler shift calibration performed in the case of Kuroshio demonstrates the feasibility and robustness of the method. In subsequent research, we will apply the method to the Gulf Stream and conduct a more systematic analysis with the collected public HF Radar data.

At present, wind-wave bias estimation and 2D ocean surface current vector inversion are still two major challenges in providing the reliable current fields by SAR Doppler shift. If we can obtain more satellite or in-situ observation data, these can be used to analyze the contribution of nonlinear hydrodynamic factors to ocean currents, and enhance the practicality of the method under complex ocean conditions. In addition, the current vector might be attempted by combining the Doppler Centroid Anomaly (DCA) method with sub-aperture decomposition. In conclusion, Sentinel-1 IW mode SAR data hold great potential in ocean current mapping. This study serves as an effective complement to non-geophysical Doppler calibration and provides support for the construction of ocean current time series data sets, making it crucial for remote sensing observations in ocean dynamics.

Acknowledgments

The authors sincerely thank the ESA for providing Sentinel-1 SAR products and all anonymous reviewers who provided detailed and valuable comments or suggestions to improve this manuscript.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement

The Sentinel-1 SAR OCN products that support the findings of this study are openly available at https://scihub.copernicus.eu/dhus/#/home.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2024.2304956.

Additional information

Funding

This research was supported in part by the National Natural Science Foundation of China [No. 42176174], Ningbo Yongjiang Talent Introduction Programme [No. 2021A-136-G] and in part by the Zhejiang Provincial Natural Science Foundation of China [Nos. LY22D010002, LR19D010001], in part by the Zhejiang Province “Pioneering Soldier” and “Leading Goose” Leading Goose’ R&D Project [No. 2023C01027].