A Simulation Method of Two-Dimensional Sea-Surface Current Field for Trajectory Crossing Spaceborne SAR

Abstract

The trajectory-crossing method is an important and effective method for spaceborne SAR (synthetic aperture radar) to detect a two-dimensional current field. However, in practical applications, there are few spaceborne SAR data that meet the requirements of close acquisition time and overlapping trajectories, which makes it difficult to comprehensively analyze the impact of different systems and environmental parameters on current measurement results from real data processing. With the proposal of a SAR constellation plan in the future, the data resources will be constantly enriched, and the trajectory-crossing method will be widely used. It is, thus, necessary to lay some theoretical foundation at present. A two-dimensional sea-surface-current-field simulation method for trajectory crossing spaceborne SAR is proposed in this paper. Based on the principles of the trajectory-crossing method, the proposed method can realize the simulation of two-dimensional sea-surface-current field under given spaceborne SAR system and environmental parameters. In this paper, the simulation process of this method is given, and the simulation experiment is performed. Compared with the current measurement results of SAR data, the simulation experiment shows that the current velocity error is less than 0.03 m/s and the direction error is less than 10 degrees, which proves the reliability of the proposed simulation method. The proposed method lays a foundation for analyzing the influence of various parameters in the application of the trajectory-crossing method.

1. Introduction

As an important ocean dynamic parameter, the sea-surface-current field has a profound impact on climate change, ocean engineering, fishery resources, energy, offshore target detection and so on [1,2]. Therefore, monitoring of ocean currents is becoming crucial. Spaceborne synthetic aperture radar (SAR) has the ability to observe the earth all day, in all weather and at high-resolution, which makes it an important remote-sensing tool that can retrieve the sea-surface-current field.

Currently, along-track interferometry (ATI) [3] and Doppler centroid analysis (DCA) [4] are the two main techniques for spaceborne SAR to detect ocean current fields; however, they can only retrieve the range velocity component of the current, i.e., they cannot directly obtain the current vector—that is, the two-dimensional current field. Generally, two sets of SAR data with crossed trajectories can be used to synthesize the velocity components in two different directions on the overlapping sea-surface area to thus obtain the current vector. This is an important two-dimensional current field detection method, which is denoted as a trajectory-crossing method in this paper.

However, due to the requirement that the overlapping area needs to cover the research scope and as the data acquisition time is close, there are fewer real spaceborne SAR data to meet this feature. The current measurement ability of the trajectory-crossing method is only verified in some airborne [5] and spaceborne [6] SAR data processing, and there is a lack of further research for the impact of different system and environmental parameters on the current measurement effect. Therefore, simulation is an effective research method. At present, to the best of knowledge, there is no specific simulation method for the trajectory crossing current measurement process.

Due to the complexity of sea-surface electromagnetic scattering, the simulation of sea-surface SAR data usually does not start from the original echo but mainly considers the modeling using wave spectrum and modulation transfer function. In 2000, Romeiser and Thompson [7] proposed a sea-surface Doppler spectrum model based on the improved combined sea-surface model, in which the wave–current interaction and the modulation transmission process of sea-surface backscattering were considered [8,9]. Then, the model was integrated into SAR ocean imaging model and widely used in sea-surface current measurement and simulation [10,11,12]. Therefore, this model provides a possibility for the simulation of two-dimensional sea-surface-current field for trajectory crossing spaceborne SAR.

Based on the principle of the sea-surface SAR imaging process and trajectory-crossing method, a two-dimensional current field simulation method is proposed in this paper. In the second part, the specific process of the proposed simulation method is given. In the third part, the simulation experiment and verification are conducted. Finally, the above is summarized.

2. Simulation Method of Two-Dimensional Sea-Surface Current Field for Trajectory Crossing Spaceborne SAR

The trajectory-crossing method needs to use two sets of spaceborne SAR data with overlapping areas, usually from ascending and descending passes, respectively, and their acquisition time should be close. Geometric relationships of different velocity components of sea-surface motion over overlapping area are shown in Figure 1, where the blue and red thick arrows indicate the flight direction of the ascending and descending satellite, respectively. The sea-surface imaging range is S and G, bounded by blue and red dashed lines, respectively. The middle is the overlapping area of S and G, where $V_{tot}$ and $D_{tot}$ represent the magnitude and direction of current velocity, respectively. Sea-surface movement produces projection components in four directions, where $V_S$ and $V_G$ represent horizontal radial velocity components along beams far from the ascending and descending satellites. $V_V$ and $V_U$ represent surface velocity components along the north–south and east–west directions, respectively. In addition, $\alpha$ represents the intersection angle between the direction of the beam center line and the east–west direction on the horizontal plane.

The simulation method proposed in this paper is mainly divided into three parts: radial velocity simulation, wave velocity simulation and current vector simulation. The flow chart is shown in Figure 2. The first step is to calculate the radial velocity using the SAR ocean imaging model, the second step is to calculate the wave motion component, and the third step is to calculate the two-dimensional current field combined with the trajectory-crossing method. The radar parameters include the signal frequency, look direction of the beam and incident angle. Platform parameters include the satellite orbit height and flight speed.

Based on the SAR ocean imaging model [7,8,9] proposed by Romeiser, University of Hamburg, Germany, the sea-surface horizontal radial velocity components along the spaceborne SAR beam look in directions with the ascending and descending passes simulated. The model supports simulation in the frequency range of 0.4 to 35 GHz. The main processes include the calculation of the wave–current interaction model, sea-surface scattering model and sea-surface SAR Doppler model.

Wave–current interaction is a physical process that must be considered in the real sea motion environment. We need to first give the sea-surface-current field and wind field parameters and then solve the action spectrum balance equation and calculate the modulated wave spectrum according to the weak hydrodynamic interaction theory. The action spectrum balance equation to be solved is [9,13,14]

$$\frac{\partial N\left(\mathit{x},\mathit{k},t\right)}{\partial t}+\left[c_g\left(\mathit{k}\right)+U\left(\mathit{x},t\right)\right]\frac{\partial N\left(\mathit{x},\mathit{k},t\right)}{\partial \mathit{x}}-\mathit{k}\frac{\partial U\left(\mathit{x},t\right)}{\partial \mathit{x}}\frac{\partial N\left(\mathit{x},\mathit{k},t\right)}{\partial \mathit{k}}=Q\left(\mathit{x},\mathit{k},t\right)$$

(1)

where t represents time, $\mathit{x}=\left(x,y\right)$ represents the sea-surface spatial position, $\mathit{k}=\left(k_x,k_y\right)$ represents the wave number, and $k_x$ and $k_y$ represent the wave numbers in the X and Y directions, respectively. N represents the action spectral density of sea-surface micro scale wave, U represents the surface current field, and $C_g$ represents the modulated wave group velocity. These variables about ocean waves are related to the conditions of sea-surface wind field. Q represents the source function:

$$Q\left(\mathit{x},\mathit{k},t\right)=\mu \left(\mathit{k}\right)N\left(\mathit{x},\mathit{k},t\right)\left(1-\frac{N\left(\mathit{x},\mathit{k},t\right)}{N_0\left(\mathit{k}\right)}\right)$$

(2)

where $N_0$ represents the action spectral density in the sea-surface equilibrium state, and $\mu$ represents the relaxation rate.

The relationship between the wave spectrum and action spectral density is [15]

$$N\left(\mathit{x},\mathit{k},t\right)=\frac{\rho {\omega }_0\left(\mathit{k}\right)}{\mathit{k}}\psi \left(\mathit{x},\mathit{k},t\right)$$

(3)

where $\psi$ represents the modulated wave spectrum. The natural angular frequency ${\omega }_0\left(\mathit{k}\right)=\sqrt{g\mathit{k}+\left({\tau }_s/\rho \right){\mathit{k}}^3}$, g is the gravitational acceleration, $\rho$ is the seawater density, and ${\tau }_s$ is the seawater surface tension.

The sea-surface scattering model is based on the improved combined sea-surface model [8,9] proposed by Romeiser and Alpers, which considers the modulation effect of large-scale wave on small-scale Bragg waves. The normalized scattering coefficient can be calculated by taking the modulated wave spectrum and the given radar parameters as the input of the scattering model [8,16]

$$\begin{array}{c}〈\sigma 〉={\sigma }^{\left(0\right)}+〈{\sigma }^{\left(2\right)}〉\\ ={\left.\sigma \right|}_{\mathit{s}=0}+\int \int \left({\left.\frac{{\partial }^2\stackrel{\wedge \vee }{\sigma }}{\partial \stackrel{\wedge }{s_p}\partial \stackrel{\vee }{s_p}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}+{\left.\frac{{\partial }^2\stackrel{\vee \wedge }{\sigma }}{\partial \stackrel{\vee }{s_p}\partial \stackrel{\wedge }{s_p}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}\right)·k_p^2\psi \left(\mathit{k}\right){\mathrm{d}}^2\mathit{k}\\ +\int \int \left({\left.\frac{{\partial }^2\stackrel{\wedge \vee }{\sigma }}{\partial \stackrel{\wedge }{s_n}\partial \stackrel{\vee }{s_n}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}+{\left.\frac{{\partial }^2\stackrel{\vee \wedge }{\sigma }}{\partial \stackrel{\vee }{s_n}\partial \stackrel{\wedge }{s_n}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}\right)·k_n^2\psi \left(\mathit{k}\right){\mathrm{d}}^2\mathit{k}\\ +\int \int \left({\left.\frac{{\partial }^2\stackrel{\wedge \vee }{\sigma }}{\partial \stackrel{\wedge }{s_p}\partial \stackrel{\vee }{s_n}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}+{\left.\frac{{\partial }^2\stackrel{\wedge \vee }{\sigma }}{\partial \stackrel{\vee }{s_n}\partial \stackrel{\wedge }{s_p}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}{\left.\frac{{\partial }^2\stackrel{\vee \wedge }{\sigma }}{\partial \stackrel{\vee }{s_p}\partial \stackrel{\wedge }{s_n}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}+{\left.\frac{{\partial }^2\stackrel{\vee \wedge }{\sigma }}{\partial \stackrel{\vee }{s_n}\partial \stackrel{\wedge }{s_p}}\right|}_{\stackrel{\wedge }{\mathit{s}}=0}\right)·k_p{k}_n\psi \left(\mathit{k}\right){\mathrm{d}}^2\mathit{k}\end{array}$$

(4)

where the symbol $〈·〉$ represents the statistical average; ${\sigma }^{\left(0\right)}$ represents the normalized backscattering coefficient when the sea-surface slope is zero; $〈{\sigma }^{\left(2\right)}〉$ is the sum of second-order Bragg scattering caused by sea-surface slope; $\mathit{s}=\left(s_p,s_n\right)$ represents the sea-surface slope vector; $k_p$ and $k_n$, respectively, represent the Bragg wave number components parallel and perpendicular to the radar look direction; $\psi \left(\mathit{k}\right)$ represents the wave number spectrum; ∧ and ∨ denote the Fourier transform and its conjugate of $\sigma$ to wave number $\mathbf{k}$; and $\wedge \vee$ and $\vee \wedge$ denote the Fourier transform and its conjugate of $\sigma$ to combined wave number $k_1-k_2$.

The Doppler spectrum and variance can be calculated using the double Gaussian Doppler spectrum model [7] proposed by Romeiser and Thompson, which divides the sea-surface echo into the superposition of the Bragg wave Doppler spectrum far from and towards the radar. Each Doppler spectrum component is defined as a Gaussian shape and expressed as

$$\begin{array}{c}\hfill S\left(f_D\right)=\frac{〈{\sigma }_{+}〉}{\sqrt{2\pi {\gamma }_{D+}^2}}{e}^{-{\left(f_D-{〈f_{D+}〉}_{\sigma }\right)}^2/\phantom{-{\left(f_D-{〈f_{D+}〉}_{\sigma }\right)}^2{\gamma }_{D+}^2}{\gamma }_{D+}^2}\\ \hfill +\frac{〈{\sigma }_{-}〉}{\sqrt{2\pi {\gamma }_{D-}^2}}{e}^{-{\left(f_D-{〈f_{D-}〉}_{\sigma }\right)}^2/\phantom{-{\left(f_D-{〈f_{D-}〉}_{\sigma }\right)}^2{\gamma }_{D-}^2}{\gamma }_{D-}^2}\end{array}$$

(5)

where ± represents the Bragg wave components away from and towards the radar, respectively; ${〈f_{D\pm }〉}_{\sigma }$ is the Doppler center weighted by the normalized backscattering coefficient; and ${\gamma }_{D\pm }$ is the variance of Doppler spectrum. The specific calculation formulas of ${〈f_{D\pm }〉}_{\sigma }$ and ${\gamma }_{D\pm }$ can be found in [7]. Then, the horizontal Doppler velocity $V_{dop}$ along the beam radial direction can be expressed as

$$V_{dop}=\frac{\pi ·{〈f_{D\pm }〉}_{\sigma }}{k_{e}sin\theta }$$

(6)

where $k_e$ is the wave number of the electromagnetic waves, and $\theta$ is the incident angle. Therefore, the radial horizontal sea-surface Doppler velocity along the ascending and descending satellite beam represented by $V_S$ and $V_G$ in Figure 2 is calculated successively by the above process.

The radial velocity is actually the Doppler velocity that reflects the average motion state of scatterers in the sea-surface resolution unit. In the range of medium incident angles, the radial Doppler velocity $u_r$ includes not only the current velocity but also the component of the sea-surface wave motion [17].

$$u_r=u_c+u_o+u_b$$

(7)

where $u_c$ is the current velocity, $u_o$ is the large-scale wave orbital velocity, and $u_b$ is the net Bragg wave phase velocity.

At present, the related wave motion model is not perfect, and thus it is difficult to extract the current velocity component by directly calculating the accurate wave velocity [10]. The authors in [18] stated that when there is no current field, the sea-surface velocity obtained by using the SAR ocean imaging model comes from the wave motion. Therefore, in the second step, the input current field velocity is set to zero, combined with the same wind field, and through the calculation process similar to the radial velocity, the wave motion velocities in the ascending and descending beam direction are obtained, which are $V_S^w$ and $V_G^w$, respectively. Since the influence of large-scale waves is considered in the previously used wave spectrum and scattering model, the simulated wave velocities include the sum of the last two velocity components of Equation (7).

By making a difference between the radial Doppler velocity and the wave velocity, the current velocity components $V_S^c$ and $V_G^c$ in the ascending and descending beam direction can be obtained:

$$V_S^c=V_S-V_S^w$$

(8)

$$V_G^c=V_G-V_G^w$$

(9)

Then, combined with the trajectory crossing current measurement geometric model in Figure 1, the east–west velocity component $V_U$ and north–south velocity component $V_V$ of the sea-surface current can be calculated:

$$V_U=\frac{V_S^c-V_G^c}{2cos\alpha }$$

(10)

$$V_V=\frac{V_S^c+V_G^c}{2sin\alpha }$$

(11)

Finally, the magnitude and direction of the current vector are, respectively, expressed as:

$$V_{tot}=\sqrt{V_U^2+V_V^2}=\frac{\sqrt{V_S^2+V_G^2-2V_S{V}_{G}cos\left(2\alpha \right)}}{sin\left(2\alpha \right)}$$

(12)

$$D_{tot}=arctan\left(\frac{V_V}{V_U}\right)=arctan\left(\frac{V_S+V_G}{\left(V_S-V_G\right)tan\alpha }\right)$$

(13)

3. Simulation Experiment and Verification

In this section, a simulation experiment was conducted based on the spaceborne SAR system and sea environment parameters. The effectiveness of the simulation method was verified by comparing the simulation results with the measurement results using SAR data. First, the spaceborne SAR data, system parameters and input sea environment data involved in the experiment are introduced. The geographical location of the spaceborne SAR data used in this paper is shown in Figure 3. This is the Mozambique strait between southeast Africa and Madagascar, where the sea surface is affected by the warm current of Mozambique, and the motion state is relatively stable. The red box indicates the stripmap SAR data range of ESA sentinel-1 satellite’s ascending pass, the acquisition time is 8 October 2019, 15:35 UTC (orbit number: 018385), the black box indicates the stripmap SAR data range of China Gaofen-3 (GF-3) satellite’s descending pass, and the acquisition time is 9 October 2019, 02:43 UTC (orbit number: 016651). Although these two SAR data were not acquired at the same time, the time interval is short, the current motion status is basically unchanged here, and these two SAR images overlap partially in this sea area, which meets the requirements of the trajectory-crossing method.

Figure 4 shows the SAR amplitude image of Sentinel-1 satellite and the horizontal radial Doppler velocity image of its secondary product. The dimmer brightness in the lower left corner of Figure 4a is due to the lower wind speed and the weaker echo energy due to the oil spill on the sea, which also results in errors in the Doppler velocity, as shown in the corresponding position in Figure 4b. The absolute values are small in the lower right part of the Doppler velocity map, and the overall velocity includes the influence of sea-surface wind and wave motion, which needs to be further removed.

Figure 5 shows the SAR amplitude image of the GF-3 satellite and the Doppler velocity estimated from the single look complex (SLC) data. The dark part of the right corner of Figure 5a is obviously affected by the sea oil film, which reduces the signal-to-noise ratio of the image. The Doppler velocity of Figure 5b is estimated from the Doppler processing in reference [19]. It can be seen that the absolute value in the lower left corner of the image is small, in addition to the anomaly in the oil film position, and the overall velocity overlaps with the influence of the wind and wave motion.

Although the overlap range of the above SAR data is small, it has clear velocity characteristics, which is conducive to the validation of simulation method. Therefore, a simulation experiment was performed in this sea area. The input system parameters are consistent with the real SAR data parameters, representing the ascending Sentinel-1 satellite and the descending GF-3 satellite, respectively, as shown in

Table 1. Simulation parameters.

Parameters	Sentinel-1 (Ascending)	GF-3 (Descending)	Units
Radar Frequency	5.4	5.4	GHz
Incident Angle	26.3	35	deg
Beam Direction	15	165	deg
Platform Speed	7590	7568	m/s
Orbit Altitude	707	713	km

. It should be noted that, at present, there are few ascending and descending data from a single satellite that can meet the requirements of irradiating the same sea area at a short time, and it is difficult to obtain them directly. In future applications, two or more satellites must cooperate to obtain the ascending and descending passes data of the same space-time sea area for current measurement. Therefore, the following dual satellite SAR data processing and system simulation are of practical significance.

The hybrid coordinate ocean model (HYCOM) reanalysis current field data and the European centre for medium-range weather forecasts (ECMWF) wind field data are obtained for the simulation input, respectively, as shown in Figure 6. The HYCOM model reanalysis data, which combines satellite altimeter data with temperature and salt data obtained by a buoy, is widely used in ocean current research [20,21,22]. The date selected for the HYCOM data is the same as Sentinel-1 SAR data acquisition date with spatial resolution of 1/12 degree. Wind field data at the same time and place were obtained from the ECMWF website with spatial resolution of 1/4 degree, and interpolated to ensure the same number of data points as the input current field. From the input data, it is known that the velocity in the overlapping area is small, and the sea-surface wind speed is mainly in the range of 2–7 m/s, which belongs to the moderate wind speed.

The given current and wind data are input into the SAR ocean imaging model to calculate the radial Doppler velocity $V_S$ under the parameters of ascending Sentinel-1 satellite and the radial Doppler velocity $V_G$ under the parameters of descending GF-3 satellite, as shown in Figure 7a,b, respectively. It can be seen that, due to the same orbital inclination of these two satellites, the distribution of the relative value of the velocity in the look direction of the two beams is approximately symmetrical; however, the absolute value is different due to the influence of wind and wave motion. Figure 7c,d shows the wave velocity when the current field is set to zero under the same wind field conditions, in which the absolute value of the velocity changes due to the change of the relationship between wind direction and beam look direction. When the wind direction is perpendicular to the beam look direction, the wave velocity is zero, which is consistent with the research conclusion in the literature [11].

According to the flow chart in Figure 2, the corresponding sea-surface radial velocity and wave velocity of ascending and descending passes in Figure 8 are subtracted, and the current vectors are calculated using the trajectory-crossing method. The simulated two-dimensional current field is shown in Figure 8. Comparing with the input current field given in Figure 6a, we found that they were in good agreement. Note that, as for the computational demand, the calculation process of this method takes only a few minutes, the calculation amount of wave spectrum is the largest. The calculation times will be different according to the amount of input current and wind data.

To quantitatively validate the experimental results, we selected eight geographic locations in the overlap area (Figure 9) and compare the simulation results with the SAR data measurement results, as shown in

Table 2. Comparison between the simulated and measured current data.

Location Point	Coordinate	Simulated Data		Measured Data
		Velocity	Direction	Velocity	Direction
1	(42.32,−18.04)	0.260 m/s	−9.66$°$	0.250 m/s	−16.74$°$
2	(42.32,−18.12)	0.293 m/s	−12.68$°$	0.261 m/s	−17.59$°$
3	(42.32,−18.20)	0.322 m/s	−19.78$°$	0.299 m/s	−24.05$°$
4	(42.32,−18.28)	0.312 m/s	−33.06$°$	0.275 m/s	−37.45$°$
5	(42.32,−18.36)	0.317 m/s	−42.84$°$	0.294 m/s	−47.56$°$
6	(42.24,−18.36)	0.295 m/s	−55.31$°$	0.263 m/s	−58.59$°$
7	(42.40,−18.36)	0.285 m/s	−47.69$°$	0.256 m/s	−59.21$°$
8	(42.40,−18.44)	0.297 m/s	−72.36$°$	0.268 m/s	−80.77$°$

. As the size of the simulated current field resolution unit is about 8 km and that estimated from Sentinel-1 and GF-3 SAR data is about 1 km, the 8 × 8 Doppler velocities around the corresponding eight latitude and longitude positions in Figure 4b and Figure 5b are averaged. After subtracting the simulated wave velocity in Figure 7 from the average data, the two-dimensional current field results of the SAR data are calculated in combination with Equations (8)–(13). It can be seen that the corresponding current velocity and direction are close. However, because the current and wind data input by the simulation will be different from the instantaneous environment in the real SAR data, there will also be discrepancies between the data in Table 2.

The statistical analysis of the experimental data is shown in Figure 10. Figure 10a shows a scatterplot of the velocity of the simulated current field compared with that of the measured current field from SAR data, with a correlation coefficient greater than 0.9 and a root mean square error (RMSE) of 0.028 m/s. Figure 10b shows the scatterplot comparing the direction of the simulated current compared with that of the measurement data. The correlation coefficient is greater than 0.98 and the root mean square error (RMSE) is 6.05 degrees. It can be seen that the simulated current field is in good agreement with the measured current of SAR data.

We also notice that there is a visible deviation between these two data, which will not only come from the difference of input data but also be limited by the model. For example, the wave spectrum is different from the real sea-surface state. Nevertheless, we believe that the simulation results are consistent with the measurement results. Therefore, the proposed method in this paper can effectively simulate the two-dimensional sea-surface-current field in the given spaceborne SAR system and marine environment parameters. The results also show that the proposed simulation method is conducive to the parameter analysis process of the trajectory-crossing method in the application of current measurement in the future.

4. Conclusions

With the continuous enrichment of spaceborne SAR data in the future, the trajectory-crossing method to measure the two-dimensional current field on the sea surface will be widely used. However, at present, it is difficult to comprehensively analyze the influence of system and environmental parameters on the current measurement effect from the real SAR data. In this paper, a simulation method of two-dimensional sea-surface-current field for trajectory crossing spaceborne SAR is proposed, which can realize two-dimensional current field simulation under the given spaceborne SAR system and environmental parameters. In the simulation process, the wave–current interaction and the influence of sea-surface backscattering on the SAR Doppler centroid frequency are considered, and the influence of the wave motion on the current measurement results is removed.

The simulation experiment selected the real spaceborne SAR system and sea-surface environment parameters, and the results are in good agreement with the SAR data processing results, which shows the effectiveness of the simulation method in this paper. Therefore, the proposed method is conducive to the parameter analysis process of the trajectory-crossing method in the application of current measurement. However, because the influence of the wave-breaking scattering mechanism under high wind speed and sea conditions was not considered, the simulation method in this paper is not suitable for two-dimensional current field simulation in extreme marine environments. In future work, more scattering mechanisms will be considered to improve the scope of application and promote the development of application research for two-dimensional sea-surface-current field measurements.