A Simple Strategy to Mitigate the Aliasing Effect in X-band Marine Radar Data: Numerical Results for a 2D Case

For moderate and high speed values of the sea surface current, an aliasing phenomenon, due to an under-sampling in the time-domain, can strongly affect the reconstruction of the sea surface elevation derived from X-band radar images. Here, we propose a de-aliasing strategy that exploits the physical information provided by the dispersion law for gravity waves. In particular, we utilize simplifying hypotheses and numerical tests with synthetic data are presented to demonstrate the effectiveness of the presented method.


Introduction
It is well known that images collected by nautical X-band radar enclose sea state information [1][2][3]. In particular, the space-time evolution of the sea surface elevation can be reconstructed from the X-band radar images arising by the backscattering phenomena of the electromagnetic fields with the sea surface ripples (e.g., capillary waves) in the microwave regime [4,5]. In fact, due to the modulation effects induced on the signal backscattered by the capillary waves, the longer surface gravity waves OPEN ACCESS may be visible by the radar. The echo received by the radar undergoes to a modulation due to the change of the angle under which the sea wave is viewed (tilt modulation) and to the electromagnetic shadowing of the sea surface induced by the higher waves preceding the wave under investigation. These modulation effects introduce a non linear spectral distortion, mainly at high wave number [1,2].
Therefore, data processing is necessary to estimate the sea-wave spectrum starting from the spectrum of the radar image (image spectrum). In particular, data processing is formulated as a linear inverse problem where, starting from the radar images sequence, the sea surface elevation is obtained as a function of two spatial variables (covering to the zone investigated by the radar) and of the time. The linear inversion scheme [1][2][3]6] can be summarized as follows: after a Fourier Transform of the time-space data, a spectral filter is used to reduce the effect of all the undesired phenomenon via a dispersion relation [1,6]. Furthermore, the use of the Modulation Transfer Function (MTF) allows the passage from the image (radar) spectrum to the wave spectrum [7,8]; finally, the resulting spectrum is Fourier anti-transformed to return to the space-time domain.
Being the inversion procedure strongly based on Fourier transform techniques, an aliasing phenomenon may arise in presence of an inappropriate time steps adopted to sample the radar image. In fact, due to the slow repetition time Δt (Δt is of order 1-2 s) of a nautical radar, it arises that the waves can be temporally under-sampled (less than twice each period, i.e., the Nyquist frequency is exceeded) and the aliasing effect occurs [9][10][11].
In this paper, by considering the simpler case of a 2D problem, i.e., the sea surface elevation is a function of only one spatial variable (denoted by x) and of the time, first we analyze how the under-sampling in time domain affects the reconstruction procedure in both the cases of presence and absence of a surface current moving in the same direction of the waves (head waves). Here, we mean the surface current as given by two contributions: the first one is the "physical" sea surface current; the second one is an "equivalent" sea surface current that arises due to the relative movement of the radar system with respect to the sea surface.
Therefore, the problem at hand is relevant to the naval context, where the characterization of the sea state with X-band radar system becomes a challenge, when a semi-displacement or a planning vessel is considered, even with head sea conditions. The high speed of the ship emphasizes the aliasing phenomenon of the data in the time domain, making it necessary to tackle the problem.
The aliasing problem has been already addressed in [9][10][11][12], where a strategy was also proposed to solve this problem. In particular, the reconstruction procedure is based on the accurate knowledge and exploitation of the gravity dispersion relation and aims at reconstructing the signal also when the frequency of the spectral energy exceeds the Nyquist frequency. This is performed thanks to the application of spectral symmetries [10] and the sampling theorem is not violated because the dispersion relation as additional physical information is used for the unambiguous reconstruction of the spectrum [11].
As general comment, it is worth noting that, here, de-aliasing is possible since we exploit information on the physical model of the phenomenon and therefore the case at hand is completely different from the classical one reported in the text books [13] where, in absence of any a priori and model information, it is not possible to estimate the aliasing and neither apply any mitigation strategy. As a matter of fact, de-aliasing techniques have been already applied in other applications as video images processing [14,15], optical imaging [16].
In this paper, we move in the same framework of [10][11][12] and we present a strategy to mitigate the effect of the aliasing phenomenon, which accounts for the gravity waves dispersion relation but in a different way compared to the strategy presented in [10][11][12]. In particular, differently from the procedure (exploiting spectral symmetries) presented in [10][11][12], here the folding problem is mitigated thanks to a procedure where the exploitation of a "virtual surface current" permits one to go beyond Nyquist"s limit and reconstruct the aliased dispersion shell. More specifically, the reconstruction of the aliased dispersion shells is obtained by simply exploiting a variable change based on the gravity waves dispersion relation. The effectiveness of the proposed strategy is shown with results against synthetic clutter and clutter-removed data and also a comparison is presented with the strategy outlined in [12].
Therefore, the paper is organized as follows. Section 2 reports briefly the data processing approach for sea state monitoring from X-band radar data. In Section 3, the aliasing problem is sketched by referring to head waves whereas a strategy to address this problem is proposed in Section 4. Reconstruction results against synthetic data are shown in Section 5 and finally, the Conclusions follow.

The Data Processing Approach
This Section briefly presents the scheme commonly exploited to reconstruct the time-spatial evolution of the sea surface elevation starting from the X-band radar images. For sake of simplicity, the inversion approach is presented in the 2D case with a sea surface elevation that is considered as a function of the time t and of only the spatial variable x. The proposed method can be generalized to the 3D full radar scene considering the spatial variables x and y. The reconstruction procedure is schematized in the block diagram of Figure 1 where each block is detailed below. The initial step consists in Fourier transforming (2D Fast Fourier Transform, 2D-FFT) the raw data images sequence so to obtain the 2D image spectrum F(k,ω). A High-Pass (HP) filter is then applied to the image spectrum in order to remove the effects due to the received signal power decay along the spatial range (x-variable) [2].
In the second step the expected linear gravity wave components are extracted from the HP filtered image spectrum F I (k,ω). More in detail, the dispersion relation of linear gravity waves, which relates the wave number k to the pulsation ω(k) and the sea surface current U, is considered: where g is the gravity constant and h the water depth [6]. In the following we assume infinite depth conditions, i.e., tanh(kh)  1.
In order to properly apply the filtering process based on the dispersion relation, an accurate estimation of the current value U is necessary. The current determination represents one of the key steps in the full inversion procedure; in fact, a bad estimate of U arises an incorrect filtering of the wave components. Several strategies can be adopted to the purpose, like the ones proposed in [17][18][19].
Once the current U has been estimated, the band-pass (BP) filter   U k G , , is built according to Equation (1) and applied to the image spectrum F I (k,ω) (see block diagram of Figure 1). As a result of the filtering procedure, the spectrum [18,19]. The minimization of the radar modulation effect introduced in the sea wave backscattered signal, allows to get the expected sea-wave spectrum F W (k,ω) starting from the filtered image spectrum as a good estimation [2,7,8]. Finally, the knowledge of the sea wave spectrum F W (k,ω) allows the calculation of the typical sea state parameters and the time-space evolution of the sea surface elevation   t x,  can be recovered by performing an inverse 2D-FFT on the function F W (k,ω).

The Effect of the Aliasing
This Section is devoted to investigating the effect of the aliasing phenomenon on the performance of the inversion procedure described in the Section 2. This analysis provides the basis to propose an effective de-aliasing strategy that will be presented in the Section 4. Let us first consider the case of absence of sea surface current, i.e., U = 0; in this case, the gravity waves dispersion relation becomes: By defining B k as the maximum wave-number for which the spectrum of the radar image    , k F can be assumed different from zero (in the test cases considered in this paper, we assume the spectrum is assumed to be different from zero when its modulus is larger than 0.1 times the peak value), . Apart from the modulation effects introduced by the specific radar sensing mechanisms [2,8], the radar image is a sampled version of the sea surface elevation ) , ( t x  , collected at uniform sampling steps t  and x  .
By defining t s     2 as sampling frequency and x k s    2 as the spatial sampling frequency, and according to the well known Nyquist"s law [13], the reconstructed wave elevation ) , ( t x  will be not affected by aliasing if: and: For typical X-band radar images, the spatial step x  is suitable to follow the marine waves; therefore in our analysis we assume condition (5) to be satisfied and in the following only the under-sampling in the time domain will be considered.
When relation (4) does not hold, the replicas of the dispersion relation centered at . This phenomenon is clearly indicated in Figure 2, where the aliasing of the spectrum is caused by the non proper choice of the measurement parameters as given in Table 1.

rad/s
In particular, the undesired folding phenomenon occurs for the range of the wave-number Figure 2), is determined through: that provides: According to Equation (7), all the waves with a wavelength will be affected by the aliasing phenomenon. For sake of simplicity, only progressive waves, i.e., waves whose spectrum has support in the first and third quadrant of the spectral plane (ω/k), are considered hereinafter. This entails that the phase velocity ω/k of these waves is positive, i.e., they propagate along the positive x-axis. In this specific case, the folding phenomenon arises two undesired shells in the second and fourth quadrant of the spectral domain (ω/k) (see Figure 2). These folded (aliased) shells are characterized by a negative phase velocity and accordingly they are representative of "wrong" (with not physical meaning) waves propagating in the direction opposite to the true ones (see Figure 2). Let us turn now to consider the case of a non null sea surface current. In particular, we consider a surface current moving in the same direction of the waves (head waves). In this case, we have to consider the dispersion law (see Figure 3):  (4), stating the absence of aliasing, is more and more difficult to be satisfied as long as the surface current value U increases. Therefore, we can state that the presence of such a current emphasizes the effect of the aliasing by amplifying the folding phenomenon. In fact, the bandwidth BP  , corresponding to the maximum value of the wave-number B k , becomes: and relation (9) permits to determine the minimum value Û of the surface current U for which the spectrum folding arises. In fact, the folding problem arises when Equation (4) does not hold and this gives: that after simple manipulations provides: When relation (11) holds, the collected data will be affected by aliasing and folded (aliased) shells arise in the second and fourth quadrant; this phenomenon is shown in the upper panel of Figure 4 (the measurement parameters are the same of Table 1). Accordingly, for a value of the surface current U larger than Û , the wave spectrum components affected by aliasing are located at Figure 3, is given by: leading to:  (13) is equal to the one specifically evaluated in the case without surface current [compare Equation (13) with Equations (6) and (7)].

The Strategy to Mitigate the Folding Problem
This Section is devoted at presenting a strategy with the aim of mitigating the effect of the spectrum folding on the sea-state reconstruction results. The proposed strategy moves within the same framework of the techniques presented in [10][11][12], where the reconstruction of the aliased dispersion shells is achieved by exploiting spectral symmetries. Here, the difference of the proposed approach compared to the ones in [10][11][12] is concerned with the procedure used to reconstruct the folded (aliased) dispersion shells. In particular, the here proposed approach can be summarized as made up of three main steps.
The first step begins with the accurate determination of the sea surface current by means of the approach proposed in [18]; such a sea surface current estimation approach has shown good rejection performances, even against data affected by the aliasing. The knowledge of the sea surface current is exploited for the evaluation of the folding point defined through Equations (12) and (13); after, a "virtual surface current" is applied with the aim of turning from the folded spectrum to a non folded spectrum.
The second step, is the same adopted in the approach described in [10][11][12] and is concerned with the exploitation of a zero-padding technique. This second step is preparatory to the third step, that consists of the passage from the zero padded/assembled spectrum to the original one by "erasing the virtual surface current".
As said above, the first step aims at assembling the sea wave spectrum in the allowable Nyquist band . Such an assembling procedure is achieved by exploiting a "virtual negative surface current v U " in order to consider the "virtual dispersion relation ) ( k  " expressed by: where v U U U  . In particular, the choice of the surface virtual current v U is performed with the aim of minimizing the quantity | ) , that defines the bandwidth of the assembled spectrum.
From a mathematical point of view, the virtual current v U has to be chosen as the one that globally minimizes the cost function | ) , With the aim of addressing the problem of the global minimization of ) (U G , we investigate the behavior of ) , ( U k  [see Equation (14)] for different values of the virtual current; Figure 6 Figure 7, where a virtual surface current v U equal to 11 m/s is used. , we perform the choice of the v U so to achieve minimum bandwidth condition given by: (this choice corresponds to the curve at v U =12.2 m/s in Figure 6).
The solution of the Equation (16) in the unknown Ũ provides the solution: The choice of the virtual current according to (17) is numerically justified from Figure 8 that depicts the function | ) , as function of the looked for current value Ũ . The value where the function ) (U G attains the minimum is equal to the one provided by Equation (17) and states the goodness of the choice in (16) and (17). Let us turn now to exploit the above reasoning to set-up the de-aliasing strategy. First of all, we note that, for the values of the current Ũ given by Equation (17), the bandwidth is equal to ; such a relation permits to state that the bandwidth is dependent only on the maximum value B k of the wave-number that one is interested to investigate. The lower panel of Figure 4 shows the overall result of the change of variable applied to the folded spectrum (reported in the upper panel of the same Figure) as result of the "application" of the virtual current chosen on the basis of Equation (17).
In this way, we can state that the unfolding spectrum technique is successfully applied when the holds; this inequality entails that the spectrum is reliably unfolded for the range of the wave-number ] . Accordingly, as long as the value of the Nyquist frequency s  decreases, the reliable wave-number range, where the unfolding of the spectrum is reliably performed, decreases. The value  permits to define the maximum bandwidth within which it is possible to obtain the unfolding of the spectrum; in fact, by substituting in the relation (8) and after simple manipulations, we obtain the value of the maximum bandwidth within which it is possible to obtain the unfolding of the spectrum as: Once the change of variables in (14) has been performed, a zero padding on the image spectrum ), (  can be applied in order to get the new sampling frequency: The final step consists in the addition, by means of a change of variables, of the (inverted) virtual surface current v U to the X band radar images in order to shift the zero padded unfolded spectrum to the original one (see Figure 9): where the the pedix N indicates the zero padded spectrum and As highlighted in Figure 9, the image spectrum   Figure 4 by applying the proposed strategy to solve the aliasing phenomenon.

Numerical Results
This Section aims at presenting results conforming the effectiveness of the strategy described in Section 4. To this end, synthetic data, generated by a fully 2D numerical wave-maker [20] reproducing the physical conditions existing in a real wave tank (see Figure 10) will be used. The spatial and time evolution of a wave train propagating on the free-surface SL  and generated through a hinged paddle c  moving with angular velocity a(t), has been investigated by a non viscid model. In this framework the mathematical statement is described through the Laplace equation for the velocity potential function.
Once the velocity potential is computed on the boundary domain, the nonlinear free-surface equations are stepped forward by a fourth-order Runge-Kutta scheme and the motion of the wavemaker is updated.
A domain decomposition technique has been used [21], with the aim to save computational time and memory effort for the simulation of long time wave evolution. Further details of the numerical model used as well as of the treatment of the free-surface can be found in [20]. A JONSWAP sea spectrum with H1/3 = 0.094 m and T0 = 1.97 s has been simulated in the numerical wave tank; a scale factor of 20 has been used to reconstruct the wave elevations corresponding to the full scale sea state (H1/3 = 1.88 m and T0 = 8.8 s). Here, H1/3 represents the significant sea surface elevation, and T0 the modal period associated with the prescribed spectrum. The JONSWAP sea spectrum is used in this paper, but the proposed approach is also suitable for other kinds of sea spectra.
The sea-state data have been generated by performing an average of the spectra of three sea states with the duration of 6 m. The data are sampled with a time-step of 0.34 s and a spatial step of 0.6 m. In particular, a total number of Nx = 6,306 spatial samples has been used leading to an extent of 3,750 m. For the time discretization Nt = 1,066 time-samples have been considered by achieving an overall acquisition time of 6 min.
The averaged data has been decimated; in fact, the samples actually used to perform the reconstruction are Nx = 630 and Nt = 32 with a step of Δx = 5.9 m and Δt = 2.4 s, according to the spectral parameters of Table 1. This leads to a spatial and time extent of 3,717 m and 76.8 s (in full scale), respectively. A surface current U of 7 m/s and directed (head waves) as the wave motion has been added to the data.
The corresponding radar data have been generated by exploiting the procedure proposed in [2], where the model of the electromagnetic scattering exploits the geometrical optics approximation and the shadowing and tilt modulation are accounted for. In particular, we simulated the electromagnetic scattering in the case of a X-band radar radiating at 9.3 GHz and located at the height of 20 m with respect to the sea zero-level.
The upper panel of Figure 11 shows the image spectrum of the data affected by the folding phenomenon: the quantity ω B (see Table 1) does not satisfy Equation (5). Clearer shells in Figure 11 represent the higher modes of the image spectrum due to modulation phenomena generated during the radar image formation [1].
The change of variable in (14) is performed by exploiting the value of the virtual current m/s 12.2 U v  m/s that allows us to achieve the value m/s 5.2 Ũ  according to Equation (17); this allows us to achieve the unfolded spectrum as shown in the lower panel of Figure 11.  (14). Figure 12 depicts the sea-wave spectrum after use of the filtering procedure described in the block diagram of Figure 1; the result is in good agreement with the reference (unfolded) sea-wave spectrum (see Figure 13) obtained from the same synthetic data at the finest spatial and time resolution (Nx = 6,306 spatial samples and Nt = 1,066 time-samples).   Figure 14 for the spatial range [2,700,3,200] m is shown in Figure 15.   Figure 13 in the spatial range of [2,700,3,200] m.
In order to give a quantitative evaluation of the effectiveness of the proposed strategy the normalized quadratic norm error: which correspond to an improvement of 30% in the quality of the reconstruction. Finally, we present a comparison between the proposed procedure and the strategy implemented according to the scheme in [12] exploiting a "shifting" procedure based on spectral symmetries. In particular, we assume to make an error of 3.5 m/s in the knowledge of the sea surface current so that we assume erroneously that 5 . 3  U m/s. Figure 16. The assembled spectrum according to the procedure in [12]. Figure 16 depicts the re-constructed spectrum according to the procedure outlined in [12]; it can be noted that in this case, due to the fact that the performances of the strategy is strongly dependent on the accurate estimation of the sea surface current, the spectrum is not well reconstructed and its support exhibits undesired "jumps". Figure 17 depicts the spectrum assembled by means of the proposed strategy. It shows how, even in presence of inaccurate knowledge of the sea surface current, it is possible to assemble the spectrum so as to avoid jumps so that the support of the spectrum behaves as a continuous function ) (k  At this point, by using this unfolded spectrum, it is simpler to re-evaluate the sea surface current and after to reconstruction of the sea surface elevation is possible. Figure 17. The assembled spectrum according to the procedure proposed herein.

Conclusions
This paper has dealt with the question of the aliasing problem that may arise in the reconstruction of the sea surface elevation starting from the images collected by a X-band radar system. In particular, the effect of a non appropriate time-step, adopted in radar data acquisition, has been thoroughly analyzed for the progressive waves. The investigation was concerned with both the cases of the absence and presence of the surface current and we have shown, by simple theoretical considerations, how the effect of an increasing value of the surface current badly affects the aliasing problem.
Based on these theoretical considerations, a strategy was proposed to mitigate the aliasing problem and its effectiveness has been shown by an analysis of synthetic data. The proposed method has been presented for unimodal sea states. As future development of this research activity, more challenging cases where the aliasing problem occurs, such as regressive and mixed waves, should be addressed and where the further question of the ambiguities arises [12]. This investigation would permit applying the method in an operational way for different sea state conditions. Finally, as future activity, we will address the operational capability of the proposed technique to the more realistic 3D case and in the real world by processing experimental data.