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

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. Remote Sens. 2015, 7 7754


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][5][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.

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 0 (x, y, t) ((x, y) denotes the position on the sub-image, and t is the time coordinate), the current vector U (U , φ 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.

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].
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 0 (k x , k y , ω), where k x and k y are the wavenumber components in the xand y-directions, respectively, and ω is the angular frequency of the wave components.

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: where g is the gravitational acceleration, d is the water depth, U is the current velocity and k is the wave vector whose length is: In this analysis, the water is assumed to be sufficiently deep for all wavenumbers, so that Equation (2) can be simplified as: 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 x p = (p − 128)∆k and k x q = (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 0 (k x , k y , ω).The output of this HP filter can be described as: 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 1 (k xp , k yq , ω).If the maximum energy of this vector is lower than a threshold P , defined by P = I 1 (k x , k y , ω)| max /2000, all data points within the vector are excluded from further processing, i.e., This process eliminates possible interference from low-energy background noise spikes and reduces the number of iterations required in subsequent steps.If the threshold 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.
is an appropriate choice.
• Again, for a certain column vector along the angular frequency axis in the spectrum I 2 (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: Finally, by grouping all extracted (k xp , k yq , ω j ) together, a two-dimensional matrix, ω 0 (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 0 (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.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.

Converting to the Polar Current Shell
The current-included dispersion relation in Equation ( 4) may be rewritten as: where θ = |θ k − φ U | denotes the intersection angle between the current direction φ U and the wave vector θ k and ω U (k x , k y ) = k • U represents the frequency shift induced by a surface current U of magnitude U and is, thus, referred to as the "current shell" in Cartesian coordinates (i.e., Cartesian current shell).The term √ gk (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 ω 0 (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, θ k ).This allows inspection of the data points along a specific radius or direction, where either θ k (i.e., θ) 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.

Determining Current Parameters
Interference due to noise may still fall in the polar current shell ω U (k, θ 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 θ is fixed so that ω 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 worth mentioning that all of the points corresponding to ω 0 (k x , k y ) = 0 should be removed from the vector k 128 ] before the test.After the interference in all radial directions is treated as discussed above, a least-squares curve fitting along the circumferential direction for each definite radius k is implemented using the model function: The current speed U and direction φ U can be obtained through least-squares curve fitting by minimizing the error cost function, 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, θ) 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.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.

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 31 = 2.5 m, mean period T 01 = 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.

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 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.

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.

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

Figure 2 .
Figure2.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 .
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 5 .
Figure 5.The FINO3 platform on which the radar is installed.

Figure 8 .
Figure 8. Significant wave height and wind speed.

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

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