Estimation of the Significant Wave Height from Marine Radar Images without External Reference

Abstract

In the context of the sea state monitoring by means of the X-band marine radar, the estimation of a significant wave height ($H_s)$ is, currently, one of the most challenging tasks. For its estimation, a calibration is usually required using an external reference, such as in situ sensors, and mainly buoys. In this paper, a method that allows us to avoid the need for an external reference for $H_s$ estimation is presented. This strategy is, mainly, based on the correlation between a raw radar image and the corresponding non-calibrated wave elevation image to which varying its amplitude by using a scale factor creates a mathematical model for the radar imaging. The proposed strategy has been validated by considering a simulated waves field, generated at varying sea state conditions. The results show a good estimation of the significant wave height, confirmed by a squared correlation coefficient greater than 0.70 for each considered sea state.

1. Introduction

X-band marine radar is a useful and well-assessed technology for sea state monitoring for both offshore and coastal areas. Within the last three decades, several algorithms, mainly based on spectral analysis, have been developed in order to retrieve the sea state parameters starting from the raw radar images [1,2,3,4,5,6]. Currently, the wave radar system is employed in many applicative scenarios since these algorithms provide accurate measurements of the sea states in terms of the period, length, and direction of the dominant waves, the significant wave height, and the reconstruction of the sea surface current and bathymetry fields [7,8,9,10,11,12,13].

In this context, one of the most interesting and challenging tasks is the estimation of the significant wave height ($H_s$). In literature, there are two classes of approaches for the $H_s$ estimation. The first class is based on the calibration of the wave spectra thanks to an external reference, such as the wave buoy. In this frame, the methods based on spectral analysis [14] or approaches based on time analysis [15,16,17] or other techniques, including an iterative least square approach [18] and a wavelet-based algorithm [19], are recalled.

The second class is about the possibility to calibrate $H_s$ without additional reference sensors [20,21,22,23]. In particular, in Wijaya, A.P. and van Groesen, E. 2016 [20], the estimation of $H_s$ is carried out by correlating the shadowing phenomenon to a visibility function stored in a database. In Gangeskar, R. 2014 [21], by assuming a geometric shadowing condition, shadowed areas are first extracted from the image by an edge detection algorithm. Afterward, by using the calculated illumination ratios in local areas, the Root Mean Square (RMS) surface slope is derived by curve fitting Smith’s function [24]. Lastly, $H_s$ is estimated from the RMS surface slope and the average zero-crossing wave period. In Liu, X. et al. 2016 [22], few improvements of the algorithm in Gangeskar, R. 2014 [21] have been presented. In particular, the use of the upwind information to fix the azimuthal location of the subarea is employed to calculate the RMS of the surface slope. The average zero-crossing wave period can be estimated from the sea spectra obtained by processing the X-band radar images. Recently, a shadowing mitigation method, which allows the estimation of wave elevation fields in coastal areas through the sea clutter data obtained from the X-band marine radar systems in extreme grazing incidence angles without calibration or the empirical Modulation Transfer Function (MTF) adjustments, has been proposed [23].

The purpose of this work is to describe a novel and alternative approach for estimating the $H_s$ measurements from the radar image without the use of an external reference. In particular, the basic idea is to estimate this sea wave parameter by performing a correlation procedure between raw radar images and the corresponding non-calibrated wave elevation images to which, varying its amplitude by means of a scale factor, a mathematical model for the radar imaging is applied [2]. Herein, we show a preliminary study of the proposed strategy by considering a simulated wave field provided by a Fourier domain approach [25] with the spectral proprieties of the Pierson-Moskowitz model [26]. The simulation of sea wave images is completed by applying the modulation effects (shadowing and tilt modulation) introduced from the particular acquisition geometric of the radar antenna (grazing angle) [2]. In order to evaluate the performance of the proposed approach, the wave field for different sea states in terms of significant wave height and peak period have been generated.

Therefore, the paper is organized as follows. Section 2 describes the generation of synthetic data and the proposed strategy. The results are presented in Section 3. Discussion and conclusions are at the end of the paper.

2. Materials and Methods

This section is devoted to describe both the mathematical models used for generating the synthetic data and the proposed algorithm for estimating $H_s$.

2.1. Generation of the Synthetics Data

2.1.1. Generation of the Wave Field

Synthetic sea wave images have been generated with the Fourier domain approaches, which was introduced by Mastin 1987 by employing measured spectral proprieties of real ocean waves [25].

The sea wave field $h\left(\overline{r},t\right)$ is modelled as a sum of sinusoids with time-dependent amplitudes, defined as [27]:

$$h\left(\overline{r},t\right)={\displaystyle \sum }_k\tilde{h}\left(\overline{k},t\right)e^{i\overline{k}\overline{·r}}$$

(1)

where $\overline{k}=\left(k_x,k_y\right)$ is the 2D wave vector whose amplitude $k=\sqrt{k_x^2+k_y^2}$ represents the wavenumber. $\overline{r}=\left(x,y\right)$ is the generic point and $t$ is the time. The set of complex Fourier domain amplitudes and their initial phase values for the wave elevation field at time zero is defined as:

$${\tilde{h}}_0\left(\overline{k}\right)=\frac{1}{\sqrt{2}}\left({\xi }_{Re}+i{\xi }_{Im}\right)\sqrt{S_{PM}\left(\overline{k}\right)}$$

(2)

where ${\xi }_{Re}$ and ${\xi }_{Im}$ are ordinary independent draws from a Gaussian random number generator, with mean 0 and standard deviation 1 [27].

In the present paper, the spectral properties of the wave field are modelled by the Pierson-Moskowitz (PM) model, $S_{PM}\left(\overline{k}\right)$, suitable for fully developed wind seas [26] and it is defined as follows:

$$S_{PM}\left(\overline{k}\right)=\frac{\alpha g^2}{k^4}{e}^{-\frac{5}{4}\left(k_p/k\right)}$$

(3)

where $\alpha =0.0081$ is the Phillips constant, $g$ is acceleration due to the gravity, and $k_p$ is the peak wave number. This latter is directly related to the wind speed at a height of 10 m above the sea level, ${\mathcal{W}}_{10},$ by means of the expression $k_p=g/{\mathcal{W}}_{10}{}^2$.

Assuming linear wave theory, sea surface waves are dispersive and their propagation mechanism are governed by the dispersion relation given by

$$\omega \left(\overline{k}\right)=\sqrt{gk \tanh \left(kd\right)}+\overline{k}·\overline{U}$$

(4)

where $\omega$ is the angular frequency, $d$ is the water depth, and $\overline{U}=\left(U_x,U_y\right)$ is the surface current vector. Given this dispersion relation $\omega \left(\overline{k}\right),$ the Fourier domain amplitudes of the wave elevation field at time $t$ can be expressed as:

$$\tilde{h}\left(\overline{k},t\right)={\tilde{h}}_0\left(\overline{k}\right)e^{i\omega \left(\overline{k}\right)t}+{\tilde{h}}_0{}^{\ast }\left(-\overline{k}\right)e^{-i\omega \left(\overline{k}\right)t}$$

(5)

This relation preserves the complex conjugation property ${\tilde{h}}^{\ast }\left(\overline{k},t\right)=\tilde{h}\left(-\overline{k},t\right)$ and is convenient because it will allows us to compute $\tilde{h}\left(\overline{k},t\right)$ and, hence, $h\left(\overline{r},t\right)$ via the inverse 2D Fast Fourier Transform (FFT) on demand at any time t without having to compute it at any other time [27].

2.1.2. Radar Imaging

Simplified X-band radar imaging models were developed by References [2,28,29]. These works identify two main mechanisms, which are considered responsible for the modulation of the Radar Cross Section (RCS) for X-band marine radars including the shadowing and the tilt modulation. Herein, for the generation of the synthetic radar images, the model proposed by Nieto et al. [2] is adopted and briefly summarized below. Other sources of RCS modulation have been considered to be either negligible, as in the case of hydrodynamic modulation or pre-filtered from the received radar signal and in the case of the local wind velocity modulation.

A shadowing effect can be interpreted, to a first approximation, as a geometrical phenomenon. At low grazing angles, the wave free surface, in a few areas, prevents the illumination of some other parts of the wave field causing the shadowing of the nearby wave facets. As a result, the shadowing phenomenon can be summarized and modelled as if the radar antenna is not actually receiving any signal from the shadowed parts of the sea surface [2,20,30,31]. Moreover, this phenomenon will be more severe with increase of the distance from the radar antenna.

Since, in the real data, the shadowed area is not equal to zero, areas are modelled as Gaussian white noise. The simulated shadowed image $M_{sh}\left(\overline{r},t\right)$ is given by the equation below.

$$M_{sh}\left(\overline{r},t\right)=\{ \begin{array}{c}1\\ noise\end{array} \begin{array}{c}if no shadowing occurs\\ otherwise\end{array}$$

(6)

Regarding the tilt modulation, the power received by the radar antenna depends on the local slope of the observed surface. This causes a modulation of the received radar signal, which depends on the angle between the radar illumination ray and the normal vector to the wave surface. The received radar signal is simulated by carrying out the scalar product between $\overline{u}\left(\overline{r},t\right)$ and $\overline{n}\left(\overline{r},t\right)$ defined as 3D vectors in a Cartesian axes system, which refers to the viewpoint of the radar and the normal vector to the surface, respectively. Therefore, the tilt simulated image $M_{tilt}\left(\overline{r},t\right)$ is given by factor $\mathcal{T}\left(\overline{r},t\right)=\overline{n}\left(\overline{r},t\right)· \overline{u}\left(\overline{r},t\right)$.

$$M_{tilt}\left(\overline{r},t\right)=\{\begin{array}{c}T\left(\overline{r},t\right)\\ noise \end{array} \begin{array}{c}if \mathcal{T}\left(\overline{r},t\right)>0 and M_{sh}\left(\overline{r},t\right)\ne noise\\ otherwise\end{array}$$

(7)

A pictorial view of the tilt model and the shadowing effect is reported in Figure 1.

2.2. Significant Wave Height Estimation Methodology

As previously mentioned, the shadowing phenomenon depends on both the $H_s$ and the radar antenna height. Specifically, assuming given and known antenna height, the effect of shadowing on the radar image increases when the significant wave height grows. Consequently, the shadowing area contains intrinsic information of the $H_s$. Based on these assumptions, an algorithm that allows us to estimate $H_s$ without an external reference was developed. In particular, the proposed algorithm is organized in three main steps whose flow diagram is shown in Figure 2.

(1)

Reconstruction of non-calibrated wave elevation by means of the inversion procedure [2].

(2)

Estimation of the gray–level threshold in order to identify the areas affected by shadowing [21,22].

(3)

Estimation of the significant wave height by means of the correlation between the raw radar image after the thresholding procedure (step 2) and the corresponding “reconstructed” radar image, which is achieved by applying the modulation effects to the non-calibrated wave elevation image retrieved during step 1.

2.2.1. Step 1: Inversion Procedure

The retrieval of the sea state parameters is carried out by exploiting the inversion procedure described by Nieto et al. 2004 [2].

Specifically, starting from the 3D Fourier Transform of the raw radar sequence, $I\left(\overline{r},t\right)$, the image spectrum $F\left(\overline{k},\omega \right)$ is achieved. In order to remove the non-stationary and non-homogenous trends in the image radar sequence, a High Pass filter is applied to the image spectrum and the spectrum $F_I\left(\overline{k},\omega \right)$ is obtained.

The second step of the procedure aims at extracting the linear components of the gravity waves from the image spectrum $F_I\left(\overline{k},\omega \right)$. To this aim, the dispersion relationship in Equation (4), which includes the parameters associated with the sea surface current and bathymetry, is used to build a Band Pass filter [32]. Therefore, the knowledge of the surface current $\overline{U}=\left(U_x,U_y\right)$ and bathymetry $d$ is a key point in the inversion procedure. In fact, an incorrect estimation of these parameters translates to an incorrect spectral filtering with a detrimental effect on estimating the sea state parameters. In this work, the normalized scalar product (NSP) [3] is used to estimate the bathymetry and the surface current. Once these parameters are estimated, it is possible to properly construct the band-pass filter (BP) $G\left(\overline{k},\omega ,d,\overline{U}\right)$ and to filter the image spectrum $F_I\left(\overline{k},\omega \right)$. The filtered image spectrum ${\tilde{F}}_I\left(\overline{k},\omega \right)$ is not yet the desired sea wave spectrum because it is affected by the radar modulation effects. Thus, a modulation transfers function (MTF) is used to reduce these kinds of distortions. This function was empirically determined by Nieto et al. [2] and performed herein. From the wave spectrum $F_w\left(\overline{k},\omega \right)$, it is possible to determine the characteristic sea state parameters in terms of direction, length, period of the dominant waves, and the significant wave height. As far as $H_s$ is concerned, this parameter is tied to the non-calibrated spectral moment of zero order $m_0$ by the following equation.

$$H_s=4\sqrt{m_0}$$

(8)

The last step of the inversion algorithm is to establish the non-calibrated wave elevation sequence $\eta \left(\overline{r},t\right)$ as a function of the time and spatial variables by applying an inverse 3D FFT to the sea wave spectrum.

2.2.2. Step 2: Estimating Shadow Threshold

As mentioned in Section 2.1.2, the radar images are affected by the shadowing, for which no information can be recovered for the sea spots that are not in the line of sight (LOS). Moreover, at given and known radar antenna height $Z_a,$ this phenomenon becomes more prominent when $H_s$ increases and the distance from the radar antenna grows. Hence, in order to discriminate the shadowed areas from the wave pattern, an edge detection procedure is applied to the radar image $I\left(x,y\right)$. Such a technique involves the convolution of $I\left(x,y\right)$ with a pixel difference operator $D_n\left(x,y\right)$ for each of $n=1, 2,\dots , 8$ directions [21,22]. The kernel of the difference operator is given by

$$D_1\left(x,y\right)=\begin{array}{ccc}0& 0& 0\\ 0& 1& -1\\ 0& 0& 0\end{array}$$

(9)

The operator takes a single kernel mask and rotates it in 45-degree increments through all eight compass directions: N, NW, W, SW, S, SE, E, and NE. Furthermore, the convolution outcome is given by:

$$I_{E_n}\left(x,y\right)=I\left(x,y\right)⨂D_n\left(x,y\right)$$

(10)

To each edge image $I_{E_n},$ a threshold value equal to the upper N-percentile of the pixels is used and the thresholding of the edge image number n is given by

$$I_{E_n}\left(x,y\right)=\{\begin{array}{c}1,\\ 0,\end{array} \begin{array}{c}{I}_{E_n}\left(x,y\right)>upper N\text{-}percentile\\ otherwise\end{array}$$

(11)

Subsequently, an overall edge image $I_T\left(x,y\right)$ is obtained by summing the eight edge images by the following two-step procedure.

$$I_T\left(x,y\right)={{\displaystyle \sum }}_{n=1}^8{I}_{E_n}\left(x,y\right)$$

(12)

$$I_F\left(x,y\right)=\{\begin{array}{c}1,\\ 0,\end{array} \begin{array}{c}{I}_T\left(x,y\right)\in \left[1,{\tau }_F\right]\\ otherwise\end{array}$$

(13)

The filtering is implemented to remove the single pixel noise that has edges in more than ${\tau }_F$ directions. The pixels of intensity value of 1 in $I_F\left(x,y\right)$ can now be used to identify the pixels that are not in the shadow in the raw radar image.

The corresponding pixel values $\rho$ from the radar image is used to create a histogram $F_D\left(\rho \right)$ of the statistics of image pixels located between the shadow and no shadow, from which the shadow threshold ${\tau }_s$ can be determined by means of the histogram mode as follows [14].

$${\tau }_s=mode\left(F_D\left(\rho \right)\right)$$

(14)

The mode in Equation (14) is the data value that occurs the most in the considered radar image. The statistical distribution as well as the shadow threshold is determined for each radar image [22].

2.2.3. Step 3: Estimating Wave Height

The idea behind the estimation of the true (not scaled) significant wave height is to correlate the radar image after the thresholding procedure (step 2), $I_{th}\left(x,y\right)$, with a “reconstructed” radar image, $\tilde{I}\left(x,y\right)$, achieved by applying the modulation effects, described in Section 2.1.2, to the corresponding reconstructed wave elevation image $\tilde{\eta }\left(\overline{r},t\right)$. Additionally, $\tilde{\eta }\left(\overline{r},t\right)$ is normalized through a dimensionless factor equal to the maximum of the considered sequence.

By assuming the radar antenna height $Z_a \mathrm{is}\mathrm{known}$, the correlation can be performed at a single time $t$ of the considered $\tilde{\eta }\left(\overline{r},t\right)$ by varying its amplitude through a scaling factor $\epsilon$.

The best scaling factor $\widehat{\epsilon }$, which will allow us to calibrate $\tilde{\eta }\left(\overline{r},t\right)$, can be estimated by looking for the maximum of the correlation function, i.e.,

$$\widehat{\epsilon }=\arg \underset{\epsilon }{\max }\frac{COV{\left(I_{th}\left(x,y\right), \tilde{I}\left(x,y\right)\right)}^2}{Var\left(I_{th}\left(x,y\right)\right)·Var\left(\tilde{I}\left(x,y\right)\right)}$$

(15)

where $\tilde{I}\left(x,y\right)=M_{sh,tilt}\left(\epsilon ·\tilde{\eta }\left(x,y\right)\right)$ for each $\epsilon =0.2, 0.4,\dots .{\epsilon }^{max}$ and $I_{th}\left(x,y\right)$ is the thresholded radar image. Therefore, the scaling factor is used to calibrate the wave elevation sequence retrieved in step 1 and ${\widehat{H}}_s$ is calculated by using Equation (8) from the corresponding calibrated spectral moment of zero order.

3. Results

In order to validate the proposed strategy, seven datasets were generated by considering a different peak period T_p, peak wavelength λ_p, significant wave height H_s, and wind speed ${\mathcal{W}}_{10}$ (see

Table 1. Input parameters for the generation of simulated sea states.

Sea State	Sea State Parameters
	Hs (m)	Tp (s)	λp (m)	${\mathcal{W}}_{10} (\mathbf{m}/\mathbf{s})$
SS1	1.05	4.38	30.02	6.84
SS2	2.20	6.20	60.05	9.68
SS3	3.22	7.59	90.08	11.85
SS4	4.33	8.77	120.11	13.69
SS5	5.59	9.80	150.13	15.31
SS6	6.51	10.74	180.16	16.77
SS7	7.25	11.60	210.19	18.11

for the details). A deep-water condition and the absence of surface current are assumed (i.e., $\tanh \left(kd\right)=1$ and $\overline{U}=0 m/s$ in Equation 4).

The time histories were generated considering a space-time grid corresponding to a typical configuration of a realistic wave radar device. Specifically, each dataset was constructed in a full circular area with an outer radius of 1004 m, and with the spatial spacing equal to $\mathsf{\Delta }\overline{r}=4m$ and $N_t=128$ time instants with the interval between two successive images being equal to $\mathsf{\Delta }t=2s$. The main propagation direction for all considered sea states is from the north ($0$ degrees) to the south. The synthetic radar images have been generated by applying the shadowing and tilt modulation models reported in Equations (6) and (7), respectively. For all considered datasets, the radar antenna was located in the middle of the radar image with given height equal to $Z_a=20 m$, which corresponds to a typical real installation [11,12,20], above the mean sea level (see Figure 1). In the following, the images related to the sea state SS3 (see Table 1) at time instant ${\mathrm{t}}_{64}=128\mathrm{s}$ are shown.

Figure 3a,b show the synthetic sea wave image, obtained by the Fourier domain approach described in Section 2.1.1, and the corresponding radar image, respectively.

For each considered dataset, the three steps of the proposed procedure, described in Section 2.2, are applied. In particular, step 1 provides the reconstruction of the non-calibrated wave elevation sequence (see Figure 4).

Afterward, in step 2, the edge detection technique is applied to the radar image to separate the shadow area from the remaining areas, according to the procedure described in Section 2.2.2.

Figure 5a shows the statistical distributions of the edge pixels (red circles) compared to the distributions of all image pixels (black circles). The shadow threshold is indicated by the dashed line. The shadowed radar image, represented in binary scale, is obtained thresholding the radar image and is depicted in Figure 5b.

The modulation phenomena are then applied to the non-calibrated wave elevation image by varying its amplitude, which is correlated to the thresholded radar image (see Equation (14)). The behaviors of the R² of each considered dataset when varying the scaling factor $\epsilon$ are plotted in Figure 6a and $\widehat{\epsilon }$ is retrieved from their maximum value, which is used to calibrate the corresponding wave elevation sequence. Figure 6b shows the calibrated wave elevation image.

In order to evaluate the robustness of the proposed strategy, it is applied to 20-time instants of each considered sea state. Specifically, the radar and non-calibrated wave elevation images from time instant ${\mathrm{t}}_{54}=108\mathrm{s}$ to ${\mathrm{t}}_{73}=146\mathrm{s}$ are considered. The behavior of the estimated ${\widehat{H}}_s$ for each time instant is depicted in Figure 7, from which it is possible to note that, with increasing significant wave height, there is more variability for the estimated values.

This trend can be explained in terms of the standard deviation and mean squared correlation coefficient that are depicted in Figure 8a,b, respectively. It is evident that the standard deviation increases and the squared correlation decreases with increasing ${\widehat{H}}_s$.

At the end of the procedure, the significant wave height is calculated using Equation (8) and considering the calibrated moment of the zero order of the sea wave spectrum and each calibrated wave elevation sequence. In addition, in Figure 9, the generated and calibrated wave elevation profiles from the south to the north for each considered sea state are shown.

4. Discussion and Conclusions

The present work proposes a strategy for estimating the significant wave height from X-band radar images without the use of external sensors in order to calibrate the sea wave spectrum. Specifically, the strategy is based on the research of the maximum correlation, in time domain and at the same time instant, between a raw radar image and the “reconstructed” radar image.

To reach such a purpose, the wave field is generated through baseline simulations aimed at reproducing different wave conditions. In particular, herein, as a preliminary study of the proposed strategy, a fully developed sea state in the deep-water condition (PM spectrum) was adopted.

The radar imaging was performed by applying the modulation effects (i.e., shadowing and tilt modulation) to the wave images in order to recreate the radar data sequence. The inversion procedure described in Section 2.2.1 was applied and the non-calibrated wave elevation sequence was provided. In order to perform the proposed strategy, two main steps were considered: i) a procedure based on the estimation of the gray-level threshold to discriminate the region of the image radar affected by the shadowing phenomenon, ii) application of the modulation effects to the corresponding non-calibrated wave elevation image when varying its amplitude by means of a scale factor. Such images are then exploited to estimate the $\widehat{\epsilon }$ value when the maximum of the correlation is reached. For each considered sea state, the proposed strategy provides a squared correlation coefficient value (see Figure 6a) greater than 0.70 and the estimated $\widehat{\epsilon }$ values were used to calibrate the corresponding wave elevation sequences $\tilde{\eta }\left(\overline{r},t\right)$.

Furthermore, it is noted that the behavior of R² decreases when increasing the ${\widehat{H}}_s$ (see Figure 8b). This behavior is likely related to two phenomena that affect the reconstruction of the calibrated wave elevation sequence: i) the sea state monitoring by means of the wave radar system deals with a finite-length time series. By performing the Fast Fourier Transform, some errors, named edge effects, occur at the beginning and at the end of the calibrated wave elevation (Figure 9); ii) the shadowing phenomenon becomes more prominent when increasing the distance from the radar antenna [33].

Lastly, the present work can be considered as a proof of concept for estimating the significant wave height without the use of an external sensor, which allowed us, then, to avoid the operation of calibration during the installation of a wave radar system.

The application of the proposed approach on real-world data is expected in the near future. Lastly, the refinement of the method and testing it in various sites is required.