Open AccessArticle

An Algorithm for Surface Current Retrieval from X-band Marine Radar Images

Chengxi Shen

¹,

Weimin Huang

^1,*

Eric W. Gill

¹,

Ruben Carrasco

and

Jochen Horstmann

Department of Electrical and Computer Engineering, Memorial University, St. John's, NF A1B 3X5, Canada

Centre for Materials and Coastal Research, Helmholtz-Zentrum Geesthacht, Max-Planck Strasse 1, D-21502 Geesthacht, Germany

Author to whom correspondence should be addressed.

Remote Sens. 2015, 7(6), 7753-7767; https://doi.org/10.3390/rs70607753

Submission received: 28 February 2015 / Revised: 17 May 2015 / Accepted: 8 June 2015 / Published: 11 June 2015

Download

Browse Figures

Graphical abstract
">
A block diagram of the current inversion procedure. The various quantities are defined in the text. ">
The power distribution at different frequency points for two adjacent wavenumber vectors. The upper sub-plot has only one prominent peak, and the corresponding frequency is retained as an element in the dispersion shell. The lower sub-plot shows multiple peaks that are noise contaminated and is thus discarded. ">
The curve fitting process along a particular radius. A clear sinusoidal function is estimated from the input data points, and the peak is seen to have a value close to the input current speed (2.5 m/s) and is located at 180 degrees (measured clockwise from true north), which is exactly the input current direction. ">
Comparison of input and retrieved current: (a) speed; (b) direction. ">
The FINO3 platform on which the radar is installed. ">
Comparison of acoustic Doppler current profiler (ADCP)-measured (at an 8-m depth) and radar-derived (polar current shell (PCS)) current: (a) speed; (b) direction. ">
Comparison of ADCP-measured (at an 8-m depth) and radar-derived (iterative LS (ILS)) current: (a) speed; (b) direction. ">
Significant wave height and wind speed. ">
Dependence of current speed error on SNR for the method of ILS (left) and PCS (right). ">

Versions Notes

Abstract

In this paper, a novel current inversion algorithm from X-band marine radar images is proposed. The routine, for which deep water is assumed, begins with 3-D FFT of the radar image sequence, followed by the extraction of the dispersion shell from the 3-D image spectrum. Next, the dispersion shell is converted to a polar current shell (PCS) using a polar coordinate transformation. After removing outliers along each radial direction of the PCS, a robust sinusoidal curve fitting is applied to the data points along each circumferential direction of the PCS. The angle corresponding to the maximum of the estimated sinusoid function is determined to be the current direction, and the amplitude of this sinusoidal function is the current speed. For validation, the algorithm is tested against both simulated radar images and field data collected by a vertically-polarized X-band system and ground-truthed with measurements from an acoustic Doppler current profiler (ADCP). From the field data, it is observed that when the current speed is less than 0.5 m/s, the root mean square differences between the radar-derived and the ADCP-measured current speed and direction are 7.3 cm/s and 32.7°, respectively. The results indicate that the proposed procedure, unlike most existing current inversion schemes, is not susceptible to high current speeds and circumvents the need to consider aliasing. Meanwhile, the relatively low computational cost makes it an excellent choice in practical marine applications.

Keywords:

marine radar; surface current; remote sensing

Graphical Abstract

1. Introduction

Sea surface current information is useful for off-shore activities, maritime safety and coastal protection. Unlike in situ current measurement instruments (e.g., ADCPs), X-band marine radars, which may be ship-borne or land-based, can provide current field information with very high spatial and temporal resolution within an area of about 4×4 square km surrounding the radar. Since Young et al. [1] successfully extracted ocean wave and current information from X-band marine radar image sequences in the 1980s, such radars have become useful tools for real-time measurements of surface wave and current parameters [2].

In the context of ocean surface, the scattering mechanism for X-band radiation is that of Bragg scatter from capillary waves (ripples) caused by the action of wind. Since these capillary waves are usually modulated by longer waves through hydrodynamic modulation (HM), tilt modulation (TM) and shadowing (SH), the longer waves (on the order of 10 s to 100 s of meters) become “visible” in the radar images [3]. However, the recorded image is not a direct representation of the ocean surface due to the non-linearity of the modulation processes [4–6]. Thus, to reconstruct the true surface elevation profile, an inversion scheme is required [1]: the radar image sequence is first transformed into the wavenumber-frequency domain using a 3-D fast Fourier transform (3-D FFT), and the surface current speed and direction can be determined by analysing the obtained image spectrum. Next, the (linear) energy associated with the ocean waves is separated from background noise through a frequency filter based on the linear dispersion relationship incorporating the estimated current velocity. Subsequently, by applying a modulation transfer function (MTF) to the image spectrum, the corresponding wave height spectrum can be obtained. It is apparent from the above procedure that an inaccurate current estimation may cause deviation in the dispersion relation filter, and consequently, the overall inversion accuracy may be significantly reduced. To this end, several current inversion algorithms have been developed, including the basic least-squares (LS) fitting technique [1], the weighted LS method [7], the iterative LS (ILS) approach [8], the dispersive surface classificator (DiSC) method [9] and the normalized scalar product (NSP) procedure [10]. Furthermore, to refine the performance of these inversion schemes, various modifications have been proposed and applied [11]. The ILS approach is efficient, but its performance depends on the choice of a threshold for the first guess of current velocity. The DiSC algorithm overcomes the spectral resolution limitation, but the literature on the application of this method to the problem of current extraction is very limited. The NSP method works effectively under both low and high current speeds of encounter (i.e., a combination of the velocities of the radar platform and current), but requires a significant amount of computation time. The reported accuracies of these methods are typically in the range of 0.1 m/s to 0.2 m/s.

In this paper, a new strategy to extract the surface current information from X-band marine radar images is presented. While most traditional current inversion techniques introduced above are based on the concept of “correlation” and are performed using the dispersion shell in Cartesian coordinates, the algorithm proposed here assumes a deep water scenario and involves a sinusoidal curve fitting process using the so-called current shell in polar coordinates. The algorithm based on this polar current shell (PCS) is explained in detail in Section 2. In Sections 3 and 4, the strategy is tested against both simulated radar images and field data collected by a vertically-polarized radar installed on the Forschungsplattformen in Nord- und Ostsee No. 3 (FINO3) offshore research platform in the North Sea. A summary and a few constraints on the algorithm appear in Section 5.

2. The PCS Current Algorithm

As with all existing current retrieval methods associated with marine radar data, time series of M (usually a power of two) consecutive rectangular sub-areas of sea clutter images must be extracted first. Based on the raw sub-image sequence I₀(x, y, t) ((x, y) denotes the position on the sub-image, and t is the time coordinate), the current vector

\vec{U} (U, ϕ_{\vec{U}})

for the selected sub-area can then be determined. The sub-image is usually chosen as a square one with a size of N × N pixels. The typical values for N and M are 128 and 32 or 64, respectively. Figure 1 shows the flow chart for the proposed current extraction algorithm. The algorithm consists of a total of four major steps, and each of these is explained in detail in the following subsections.

2.1. Generating Image Spectrum

The first major step involves the determination of the 3-D image power spectrum as the square of the magnitude of the 3-D FFT of the raw sub-image sequence. Although this step appears in all existing current inversion algorithms, two extra procedures are applied prior to the 3-D FFT. Firstly, the sub-image sequence is multiplied by a smooth tapering function, W (n), in what may be described as a one-dimensional discrete case [12].

W (n) = {\begin{array}{l} \frac{1}{2} (1 + \cos (π (\frac{2 n}{a L} - 1))), & 0 \leq n \leq \frac{a L}{2} \\ 1, & \frac{a L}{2} \leq n \leq L (1 - \frac{a}{2}) \\ \frac{1}{2} (1 + \cos (π (\frac{2 n}{a L} - \frac{2}{a} + 1))), & L (1 - \frac{a}{2}) \leq n \leq L \end{array}

(1)

where L is the number of grids in one dimension, and the parameter a is empirically chosen. In our case, L is extended to three dimensions (L being N, N and M for the X, Y and t direction, respectively), and a has a value of 0.1. This process can mitigate the Gibbs’ phenomenon.

After applying the taper, zero padding is used beyond the original data size of N × N × M, so that the 3-D FFT involves 256 × 256 × 256 points. This process may reduce the noise level and produce a smoother image spectrum, and it also increases the number of data points for the later curve fitting process, which, in turn, significantly improves the performance of the algorithm. The resulting image power spectrum is denoted by I₀(k_x, k_y, ω), where k_x and k_y are the wavenumber components in the x- and y-directions, respectively, and ω is the angular frequency of the wave components.

2.2. Extracting the Dispersion Shell

The next major step is to locate the dispersion shell (the surface in the ω − k_x − k_y space represented by the dispersion relationship) from the 3-D image spectra based on the linear (fundamental) mode wave dispersion relationship, which is given by:

ω = \sqrt{g k \tanh (k d)} + \vec{k} \cdot \vec{U}

(2)

where g is the gravitational acceleration, d is the water depth,

\vec{U}

is the current velocity and

\vec{k}

is the wave vector whose length is:

k = \sqrt{k_{x}^{2} + k_{y}^{2}}

(3)

In this analysis, the water is assumed to be sufficiently deep for all wavenumbers, so that Equation (2) can be simplified as:

ω = \sqrt{g k} + \vec{k} \cdot \vec{U}

(4)

Equation (4) indicates that for each pair of (k_xp, k_yq) (here, p, q = 1, 2, …, 256) in the wavenumber plane, there exists a corresponding ω_j associated with the ocean wave components. It should be noted that k_xp = (p − 128)Δk and k_xq = (q − 128)Δk, where Δk is the wavenumber resolution of the image spectrum. Due to the influence of current, ω_j contains the Doppler shift that causes the wave energy in the 3-D image spectra to deviate from the ideal (current-free) dispersion shell. Most existing current algorithms determine the current velocity as that minimizing the square sum of the difference between the calculated and observed deviations from the spectral signal. Here, current parameters are estimated based on curve fitting the spectral points on the dispersion shell. In order to locate a precise shell of ω_j(k_xp, k_yq) from the image power spectrum, the following procedure is implemented:

In order to eliminate the non-stationary and non-homogeneous trends in the sub-image sequence, a high-pass (HP) filter is first applied to I₀(k_x, k_y, ω). The output of this HP filter can be described as:

$I_{1} (k_{x}, k_{y}, ω) = {\begin{array}{l} 0, & ω < ω_{c u t} \\ I_{0} (k_{x}, k_{y}, ω) & ω \geq ω_{c u t} \end{array}$

(5)

where the cut-off frequency is usually taken as ω_cut = 0.03 × 2π radian/s [13].
For every pair of (k_xp, k_yq), p, q = 1, 2, …, 256, in the wavenumber plane of the 3-D image spectrum, there exists a column vector along the angular frequency (ω) axis that can be denoted by I₁(k_xp, k_yq, ω). If the maximum energy of this vector is lower than a threshold $\bar{P}$ , defined by $\bar{P} = I_{1} (k_{x}, k_{y}, ω) |_{\max} / 2000$ , all data points within the vector are excluded from further processing, i.e.,

$I_{2} (k_{x p}, k_{y q}, ω) = {\begin{array}{l} 0, & I_{1} (k_{x p}, k_{y q}, ω) |_{\max} < \bar{P} \\ I_{0} (k_{x p}, k_{y q}, ω) & I_{1} (k_{x p}, k_{y q}, ω) |_{\max} \geq \bar{P} \end{array}$

(6)

This process eliminates possible interference from low-energy background noise spikes and reduces the number of iterations required in subsequent steps. If the threshold $\bar{P}$ is too high, many spectral points corresponding to the wave contribution may be lost. If it is too low, points associated with noise may be included in the subsequent processing. $\bar{P} = I_{1} (k_{x}, k_{y}, ω) |_{\max} / 2000$ is an appropriate choice.
Again, for a certain column vector along the angular frequency axis in the spectrum I₂(k_xp, k_yq, ω), only one frequency point, ω_j, is associated with the ocean wave component, i.e., locates on the dispersion shell. To determine this ω_j, all energy peaks within the vector are identified by examining the derivative of each point. If only one prominent peak is detected (the criterion here being that no other peak reaches one third the energy of the main peak within the column vector), the corresponding ω is extracted as ω_j; otherwise, ω_j is set to zero. For instance, Figure 2 demonstrates how energy is distributed along angular frequency for two adjacent wave vectors. In the upper sub-figure, the only dominant peak appears at 71 out of all 128 frequency points, indicating that ω_j = (71 − 0.5) · Δω, where Δω is the Fourier transform angular frequency resolution. The lower sub-figure, on the other hand, shows multiple irregular peaks, suggesting severe distortion caused by background noise or aliases. For this case, ω_j is set to zero. Generally, we have:

$ω_{0} (k_{x}, k_{y}) = {\begin{array}{l} ω_{j}, & if I_{2} (k_{x}, k_{y}, ω_{j}) is the only peak, \\ 0 & otherwise \end{array}$

(7)

Finally, by grouping all extracted (k_xp, k_yq, ω_j) together, a two-dimensional matrix, ω₀(k_x, k_y), is constructed. The non-zero elements in this matrix (of size 256 × 256) can be viewed as a coarse estimation of the dispersion shell.

Having implemented the aforementioned three steps, a rough dispersion shell is extracted from the image power spectrum, I₀(k_xp, k_yq, ω). It should be noted that both the location and the number of extracted shell points affect the retrieval accuracy. Whether the extracted dispersion is good enough to produce a current result will be discussed in Section 2.4.

2.3. Converting to the Polar Current Shell

The current-included dispersion relation in Equation (4) may be rewritten as:

ω_{U} (k_{x}, k_{y}) = \vec{k} \cdot \vec{U} = k U \cos \bar{θ} = ω_{0} (k_{x}, k_{y}) - \sqrt{g k}

(8)

where

\bar{θ} = | θ_{\vec{k}} - ϕ_{\vec{U}} |

denotes the intersection angle between the current direction

ϕ_{\vec{U}}

and the wave vector

θ_{\vec{k}}

and

ω_{U} (k_{x}, k_{y}) = \vec{k} \cdot \vec{U}

represents the frequency shift induced by a surface current

\vec{U}

of magnitude U and is, thus, referred to as the “current shell” in Cartesian coordinates (i.e., Cartesian current shell). The term

\sqrt{g k}

(for deep water waves) arises from the ocean wave components. It should be noted that ω(k_x, k_y) in Equation (4) has been replaced with the ω₀(k_x, k_y) obtained by the procedure mentioned in Section 2.2, since only the extracted dispersion shell is of interest here.

Next, the Cartesian current shell ω_U(k_x, k_y) is converted into polar coordinates, the result being referred to as the polar current shell,

ω_{U} (k, θ_{\vec{k}})

. This allows inspection of the data points along a specific radius or direction, where either

θ_{\vec{k}}

(i.e.,

\bar{θ}

) or U is fixed, so that only one unknown may be dealt with at a time. In order to obtain reasonable accuracy, the polar shell should have at least a size of 128 × 360; i.e., 128 points in the radial direction and 360 angular elements.

2.4. Determining Current Parameters

Interference due to noise may still fall in the polar current shell

ω_{U} (k, θ_{\vec{k}})

. The interference can be regarded as outliers that may severely affect the later least-squares fit and, thus, must be removed. According to Equation (8), for a particular radial direction in the wavenumber plane, the intersection angle

\bar{θ}

is fixed so that

\frac{ω_{U}}{k}

should remain constant. This immediately suggests the application of Grubbs’ test [14] to detect one outlier at a time, and this outlier is removed and the test iterated until no outliers are detected [15]. In addition, data points due to aliasing may be eliminated by the process. It is mentioning that all of the points corresponding to ω₀(k_x, k_y) = 0 should be removed from the

[\frac{ω_{U 1}}{k_{1}} \frac{ω_{U 2}}{k_{2}} \dots \frac{ω_{U 128}}{k_{127}}]

before the test.

After the interference in all radial directions is treated as discussed above, a least-squares curve along the circumferential direction for each definite radius k is implemented using the model function:

f (U, \bar{θ}) = \frac{ω_{U}}{k} = U \cos \bar{θ} = U \cos (θ_{\vec{k}} - ϕ_{\vec{U}})

(9)

The current speed U and direction

ϕ_{\vec{U}}

can be obtained through least-squares curve fitting by minimizing the error cost function,

E = \sum_{i = 1}^{N_{c}} {[\frac{ω_{U i}}{k} - U \cos ({(θ_{\vec{k}})}_{i} - ϕ_{\vec{U}})]}^{2}

(10)

in which N_c is the number of extracted shell points along the circumferential direction for a specific k. Using this process, robust results can still be obtained even when some dispersion shell points are not extracted. It should be noted that the data points on a few short radii will not be used for the fitting procedure, since those points are less stable due to the previous coordinate conversion. Moreover, if N_c for a specific k on the extracted current shell is smaller than 10, the curve fitting for that k is aborted, and no result is retrieved for that k. If the data points are fewer than 10 for every k in a current shell, no result is retrieved for the associated image sequence. The algorithm does not require a minimum number of valid k’s, but does require a minimum number of data points for any k that is to be used in the retrieval process. The angle corresponding to the maximum of the estimated function

f (U, \bar{θ})

is determined to be the current direction, and the amplitude of this sinusoidal function is the current speed.

After the results are averaged over a band of varying radii, a good estimate of current parameters may be obtained. An example demonstrating this curve fitting process is shown in Figure 3. It can be inferred from Figure 3 that current parameters can be retrieved successfully without the necessity of acquiring all of the current shell data points.

3. Simulation Results

The proposed inversion algorithm is first tested using 30 sets of simulated X-band radar data. Each set has 32 continuous sea clutter images, representing a constant local sea state, but different current parameters. To simulate the time-varying ocean surface, the modified Pierson–Moskowitz (P-M) spectrum [16] and a cardioid directional distribution are employed [17], and the current contribution is added in accordance with the dispersion relation in Equation (4). Specifically, for the P-M wave spectrum, a wave height of H₃₁ = 2.5 m, mean period T₀₁ = 8 s and a dominant wave propagation direction of 330° are chosen; for the surface current, the speed is varied from 0.5–15 m/s in steps of 0.5 m/s, and the direction is fixed at 180°. In addition, by setting the antenna rotation speed and height, respectively, as 48 rpm and 20 meters above the sea level, the effects of tilt modulation and shadowing are calculated and applied to the surface elevation profile (e.g., see [18]). Due to the aliasing caused by low sampling rate or non-linear effects imposed by shadowing and tilt modulation involved in the simulation, outliers may also appear in advance of the curve fitting. Moreover, the dispersion shell cannot be fully populated even for simulated data, since the spectral points associated with low power wave components are eliminated by thresholding.

The comparison of the retrieved and input current parameters is illustrated in Figure 4. It may be observed from this figure that the algorithm may be used to successfully extract current information from the radar images. It is particularly worth noticing that the inversion accuracy is not compromised at large current speeds (e.g., U > 6 m/s). Moreover, the proposed routine requires a relatively short computation time (5–10 seconds for each image set) and is, thus, able to provide real-time measurements. These promising features are further examined in the next section, where results obtained from field tests are considered.

4. Experimental Results

A vertically-polarized (V-pol) X-band marine radar was mounted by the Helmholtz-Zentrum Geesthacht Centre on the offshore FINO3 research platform (55°11.7′N, 7°9.5′E; see Figure 5) in the North Sea in late July and August of 2014. The water depth within the radar coverage is 28 m and can be roughly treated as deep water for the dispersion relationship, since the observed peak wavelength is about 56 m. The antenna that was mounted 45 m above sea level was operated with a rotation period of approximately 2.14 seconds and a range resolution of 7.5 m. One radar sequence containing 128 images (about 4.5 minutes) was collected every half an hour. A dataset of about 10 days in length was used here. The size of the sub-images that were extracted from the up-wave direction for current information retrieval is 128 × 128. The sub-image is chosen from the up-wave direction, since the wave signature in the up-wave direction is usually stronger than that in other look directions on the radar image.

The comparison between the bottom-mounted ADCP data and the radar results obtained using the routine presented here and the ILS approach is shown in Figures 6 and 7, respectively. In this work, the ADCP data collected for a depth of 8 m were used since the measurement close to the surface is unreliable. It can be seen from the ADCP data that the current speed in the area was lower than 0.5 m/s and a strong tidal signal with a 12-hour period exists in the area. It should be noted that the ADCP produced one measurement every 30 minutes, but to minimize power consumption, each ADCP measurement represents an averaged value over only one minute. Each radar-derived current was an averaged result over 4.5 minutes. This may be the reason why the ADCP data look less smooth than the radar results. Good agreement of the estimated current information with that recorded by the ADCP was obtained for most of the observation period. Although the retrieved results appear to be less stable at the end of experiment, this is actually due to the low sea state. According to the wave buoy and anemometer data collected at the time (see Figure 8), the local wind speed was lower than 2 m/s, and the significant wave height was less than 0.6 m. Under such conditions, the signal-to-noise ratio (SNR) of the radar images tends to be low. Thus, the number of data points qualifying for curve fitting is less than the predefined threshold, and the performance of the inversion process is likely to suffer. In [13], it was suggested that the sea surface wind be higher than 3 m/s for an horizontally-polarized (H-pol) X-band radar to provide reliable current and wave measurements. Here, the radar-derived current speed and direction using the PCS algorithm from the data without enough useful points for curve fitting were set to zero (outliers), and only six out 494 retrieved results were found as outliers. For the ILS approach, the current results derived from radar data with an SNR lower than 4 dB were regarded as unreliable, resulting in eight outliers.

From Figure 6, it can be observed that the radar-derived results using the PCS algorithm also agree well with the ADCP data for most of the low sea state period (the beginning and end of the dataset). This indicates that the PCS algorithm is less influenced by low sea states than some of the typical inversion schemes mentioned in the Introduction. By excluding the outliers, the bias and root mean square (RMS) difference of radar-derived current speed are −0.8 cm/s and 7.3 cm/s, respectively. For the current direction, they are 0.3° and 32.7°, respectively. The correlation coefficients of current speed and direction are 0.75 and 0.76, respectively. The bias and RMS differences for the current speed and direction associated with the ILS approach are −0.3 cm/s and 7.5 cm/s and 2.5° and 33.2°, respectively. The results indicate the performances of the two methods to be very similar. The results show that the accuracies of these two current methods are very similar. The dependence of error on SNR is shown in Figures 9 and 10 for current direction and speed, respectively. It may be generally observed from Figure 9 that, not surprisingly, the speed estimation error decreases when SNR increases. However, no obvious dependence of the direction estimation error on SNR can be seen from Figure 10. This may be because only the data with an SNR greater than 4 dB are depicted, and the SNR is high enough to obtain a good result. It should also be noted that the direction estimation error depends on the current magnitude. As noted above, for the data used here, the current speed is less than 0.5 m/s.

5. Conclusion

A new algorithm to retrieve current information from X-band marine radar images is presented. In particular, a robust sinusoidal curve fitting process is involved, and the major portion of the current algorithm is performed based on the current shell in polar, rather than in the traditional Cartesian coordinates. The technique has been first tested using simulated radar images and then validated through a comparison of the derived current results from the field data with those recorded by an ADCP. Current velocities and directions have been successfully retrieved from both simulation tests and field experiments. The RMS difference for current speed and direction obtained based on field radar data are 7.3 cm/s and 32.7°, respectively.

By applying curve fitting to the polar current shell data, the algorithm is robust, even when only a portion of the dispersion shell is extracted. Furthermore, aliasing does not need to be considered. Moreover, unlike most existing inversion algorithms, which do not work properly in cases associated with high current speeds (usually when the antenna is mounted on a fast-moving vessel) [8,10], the one presented here has no clear restrictions on the current speed or direction. It is also less affected by low sea states due to the involvement of a robust curve fitting process. In light of this, the algorithm in this paper can serve as a good alternative in X-band marine radar remote sensing. More importantly, due to its low computational cost and almost real-time performance, this scheme shows good promise for practical marine applications.

A couple of limitations should be noted. Firstly, the water depth effect is ignored for this algorithm. For the case of shallow water, the dispersion shell described in Equation (4) needs to be modified. Application of the PCS algorithm to shallow water data should be undertaken in the future. Secondly, the validation of the algorithm is based only on stationary platform radar data and simulated data here. Further testing should be undertaken by using radar data collected from a fast moving ship. In further validating the robustness of the PCS algorithm, it may be useful to pursue a comparison of the tidal ellipse parameters from the ADCP data and those from the radar-derived current results.

Acknowledgments

This work was supported by the Research and Development Corporation (RDC) Industrial Research and Innovation Fund (IRIF) Ignite grant (207765), a Natural Sciences and Engineering Research Council of Canada Discovery Grant (NSERC 402313-2012) to Weiming Huang, as well as NSERC Discovery Grant (238263-2010) and an Atlantic Innovation Fund (AIF) to Eric W. Gill. The authors would like to thank Joerg Seemann at Helmholtz-Zentrum Geesthacht in Germany for his time and input to help improve the manuscript.

Author Contributions

Chengxi Shen, under the supervision of Weimin Huang and Eric Gill, proposed the PCS algorithm, provided the initial program code and prepared part of the manuscript. Weimin Huang modified the PCS code, generated the simulation and experiment results, interpreted the results and wrote the manuscript. Eric Gill coordinated revisions of the manuscript. Ruben Carrasco and Jocehn Horstmann collected field data, generated plots from field data, interpreted the results and coordinated revisions of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Young, I.R.; Rosenthal, W.; Ziemer, F. Three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents. J. Geophys. Res. 1985, 90, 681–682. [Google Scholar]
Ziemer, F. An instrument for the survey of the directionality of the ocean wave field, Proceedings of The WMO/IOC Workshop Operational Ocean Monitoring Using Surface-Based Radars, Geneva, Switzerland, 6–9 March 1995; pp. 81–87.
Plant, W.J. The modulation transfer function: Concept and applications. In Radar Scattering from Modulated Wind Waves; Komen, G.J., Oost, W.A., Eds.; Springer Netherlands: Berlin, Germany, 1989; pp. 155–172. [Google Scholar]
Lee, P.H.Y.; Barter, J.D.; Beach, K.L.; Hindman, C.L.; Lade, B.M.; Rungaldier, H.; Shelton, J.C.; Williams, A.B.; Yee, R.; Yuen, H.C. X band microwave backscattering from ocean waves. J. Geophys. Res. 1995, 95, 2591–2611. [Google Scholar]
Plant, W.J.; Keller, W.C. Evidence of Bragg scattering in microwave Doppler spectra of sea return. J. Geophys. Res. 1990, 95, 16299–16310. [Google Scholar]
Wenzel, L.B.; Electromagnetic, Scattering. from the Sea at Low Grazing Angles. In Surface Waves and Fluxes; Geernaert, G.L., Plant, W.J., Eds.; Kluwer: Norwell, MA, USA, 1990; pp. 41–108. [Google Scholar]
Gangeskar, R. Ocean current estimated from X-band radar sea surface images. IEEE Trans. Geosci. Remote Sens. 2002, 40, 783–792. [Google Scholar]
Senet, C.M.; Seemann, J.; Ziemer, F. The near-surface current velocity determined from image sequences of the sea surface. IEEE Trans. Geosci. Remote Sens. 2001, 39, 492–505. [Google Scholar]
Senet, C.M.; Seemann, J.; Flampouris, S.; Ziemer, F. Determination of bathymetric and current maps by the method DiSC based on the analysis of nautical X-band radar image sequences of the sea surface. IEEE Trans. Geosci. Remote Sens. 2007, 46, 2267–2279. [Google Scholar]
Serafino, F.; Lugni, C.; Soldovieri, F. A novel strategy for the surface current determination from marine X-Band radar data. IEEE Geosci. Remote Sens. Lett. 2010, 7, 231–235. [Google Scholar]
Huang, W.; Gill, E. Surface current measurement under low sea state using dual polarized X-band nautical radar. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2012, 5, 1868–1873. [Google Scholar]
Tukey, J.W. An introduction to the calculations of numerical spectrum analysis. In Spectral Analysis of Time Series; Wiley: New York, NY, USA, 1967; pp. 25–46. [Google Scholar]
Nieto-Borge, J.C.; Guedes-Soares, C. Analysis of directional wave fields using X-band navigation radar. Coastal Eng. 2000, 40, 375–391. [Google Scholar]
Grubbs, F.E. Procedures for detecting outlying observations in samples. Technometrics 1969, 11, 1–10. [Google Scholar]
Shen, C.; Huang, W.; Gill, E. An alternative method for surface current extraction from X-band marine radar images, Proceeding of 2014 IEEE International Geosience and Remote Sensing Symposium (IGARSS), Quebec, QC, Canada, 13–18 July 2014; pp. 4370–4373.
Nomiyama, D.H.; Hirayama, T. Evaluation of marine radar as an ocean-wave-field detector through full numerical simulation. J. Marine Sci. Technol. 2003, 8, 88–98. [Google Scholar]
Huang, W.; Gill, E. Simulation analysis of sea surface current extraction from microwave nautical radar images, Proceedings of 2012 19th IEEE International Conference on Image Processing (ICIP), Orlando, FL, USA, 30 September–3 October 2012; pp. 2673–2676.
Nieto-Borge, J.C.; Rodriguez, R.G.; Hessner, K.; Gonzales, I.P. Inversion of marine radar images for surface wave analysis. J. Atmos. Ocean. Technol. 2004, 21, 1291–1300. [Google Scholar]

Figure 1. A block diagram of the current inversion procedure. The various quantities are defined in the text.

Figure 2. The power distribution at different frequency points for two adjacent wavenumber vectors. The upper sub-plot has only one prominent peak, and the corresponding frequency is retained as an element in the dispersion shell. The lower sub-plot shows multiple peaks that are noise contaminated and is thus discarded.

Figure 3. The curve fitting process along a particular radius. A clear sinusoidal function is estimated from the input data points, and the peak is seen to have a value close to the input current speed (2.5 m/s) and is located at 180 degrees (measured clockwise from true north), which is exactly the input current direction.

Figure 4. Comparison of input and retrieved current: (a) speed; (b) direction.

Figure 5. The FINO3 platform on which the radar is installed.

Figure 6. Comparison of acoustic Doppler current profiler (ADCP)-measured (at an 8-m depth) and radar-derived (polar current shell (PCS)) current: (a) speed; (b) direction.

Figure 7. Comparison of ADCP-measured (at an 8-m depth) and radar-derived (iterative LS (ILS)) current: (a) speed; (b) direction.

Figure 8. Significant wave height and wind speed.

Figure 9. Dependence of current speed error on SNR for the method of ILS (left) and PCS (right).

Figure 10. Dependence of current direction error on SNR for the method of ILS (left) and PCS (right).

© 2015 by the authors; licensee MDPI, Basel, Switzerland This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Shen, C.; Huang, W.; Gill, E.W.; Carrasco, R.; Horstmann, J. An Algorithm for Surface Current Retrieval from X-band Marine Radar Images. Remote Sens. 2015, 7, 7753-7767. https://doi.org/10.3390/rs70607753

AMA Style

Shen C, Huang W, Gill EW, Carrasco R, Horstmann J. An Algorithm for Surface Current Retrieval from X-band Marine Radar Images. Remote Sensing. 2015; 7(6):7753-7767. https://doi.org/10.3390/rs70607753

Chicago/Turabian Style

Shen, Chengxi, Weimin Huang, Eric W. Gill, Ruben Carrasco, and Jochen Horstmann. 2015. "An Algorithm for Surface Current Retrieval from X-band Marine Radar Images" Remote Sensing 7, no. 6: 7753-7767. https://doi.org/10.3390/rs70607753

APA Style

Shen, C., Huang, W., Gill, E. W., Carrasco, R., & Horstmann, J. (2015). An Algorithm for Surface Current Retrieval from X-band Marine Radar Images. Remote Sensing, 7(6), 7753-7767. https://doi.org/10.3390/rs70607753

Article Menu