The ordinary nautical X-band radar is an active imaging radar system working in microwave frequency band. It transmits electromagnetic waves to the sea surface and generates images by receiving electromagnetic waves reflecting from the sea surface. Generally, the gray value of radar images represents the backscattering strength of echoes from the sea surface. The backscattering is mainly caused by centimeter-scale surface roughness elements. The long surface waves will modulate the small-scale scattering elements and affect the radar backscattering.

At moderate incidence angles, tilt and hydrodynamic modulations are well-established imaging mechanisms, while at grazing incidence, the modulation stems also from the shadowing of the radar beam due to the ocean waves (Alpers et al., 1981; Seemann et al., 1997). In addition, the non-Bragg scattering will play an important role when the incidence angle is large (≥70°) (Plant and Keller, 1990). Furthermore, wave breaking is believed to have contributions to radar backscattering at grazing incidence especially for HH polarization (Lee et al., 1995).

Since the 1960s, scientists have tried to retrieve ocean wave parameters from nautical X-band radar images. Ijima et al. (1964) and Wright (1965) are among the first to report the use of marine radar for wave observation. Subsequently, scientists began to digitize the radar images and employ two-dimensional (2-D) Fourier transform to marine radar images. They also found that the inverted spectra are similar to the spectra obtained from conventional buoy data (Hoogeboom and Rosenthal, 1982; Ziemer et al., 1983). However, owing to the point symmetry of spatial power spectra, the resulting spectra showed a 180 ambiguity in the direction of wave propagation. Later on, the problem of the 180 ambiguity is removed by using two successive images (Atanassov et al., 1985).

Young et al. (1985) developed an analysis technique for a full time series of radar images to retrieve an energy density spectrum through three-dimensional (3-D) Fourier transform. On the basis of an iterative least-squares fitting technique and an error-estimation model, Senet et al. (2001) proposed a method to improve the accuracy of the estimated current velocities. The method to obtain wave height information from SAR images was first applied to X-band radar image sequences by Ziemer and Gunther (1994). By assuming that wave fields have characteristics of spatial homogeneity and temporal stationary, Nieto-Borge et al. (2004) developed a linear modulation transfer function. Chen et al. (2015) proposed a new quadratic polynomial modulation transfer function based on the VV-polarized radar measurements taken from heterogeneous nearshore wave fields. Furthermore, the dual polarized and coherent marine X-band radars have been developed as new-generation devices (Cui et al. 2010; Trizna, 2011; Nwogu and Lyzenga, 2010). Nowadays, nautical X-band radar has been exploited in numerous oceanographic studies, including surface waves (Nieto-Borge and Soares, 2000; Hessner et al., 2001; Wang et al., 2007), near-surface currents (Senet et al., 2001; Gangeskar, 2002), bathymetry (Bell, 1999; Wang et al., 2015a), surface elevation maps (Nieto-Borge et al., 2004; Wang et al., 2015b), wind vectors (Dankert et al., 2003; Wang and Duan, 2010; Lund et al., 2012), and wave breaking (Catalan et al., 2011).

To apply Fourier transform in the image analysis, the assumption of spatial homogeneity (or stationary) within the analyzed area is necessary (Wu et al., 2011a). However, wave fields usually are non-stationary and inhomogeneous due to wave refraction, shoaling, and other coastal mechanisms in nearshore areas. The wavelet transform (WT) is now recognized as a flexible and effective technique to analyze non-stationary and inhomogeneous signals, as well as radar images of surface waves (Chandler, 2005). Doong et al. (2006) employed FT and WT on simulated images, respectively, and found that there are differences between the FT-spectra and WT-spectra. They suggested that an inhomogeneous image analysis method is necessary to analyze the nautical radar images. Chuang at al. (2008) developed an algorithm of applying a 2-D WT in simulated wave fields and showed that 2-D WT is useful for identifying non-stationary and inhomogeneous wave properties in shallow water. Feng et al. (2011) introduced 2-D WT to capture the inhomogeneity of wave fields from nearshore X-band radar images. Their results show that WT-based wave number spectra can present the shoaling of nearshore waves. Wu at al. (2011a) implemented 2-D continuous wavelet transforms (CWT) on satellite and X-band radar images to calculate local image spectra of inhomogeneous wave fields. They found that the local image spectra reveal the phenomena of wave refraction and nonlinearity in nearshore areas. Compared with Fourier spectra, the wavelet spectra present more accurate one-dimensional (1-D) spectra. Wu et al. (2011b) developed a procedure for quantifying the inhomogeneity from the ocean remote-sensing images by implementing the 2-D CWT and recommended an inhomogeneous algorithm for obtaining wave parameters from radar images. An et al. (2015) proposed a self-adaptive 2-D CWT-based algorithm for extracting wave information from nautical X-band radar images and found that the wavelet scaling parameters will affect the results of wave field analysis.

The signatures of ocean waves in nautical X-band radar images mainly exhibit as strip shapes and show obvious scale and directional characteristics. Owing to the effect of refraction, shoaling, and other coastal mechanisms, the wave fields usually are irregular and inhomogeneous in nearshore areas. The wave crest lines usually are curved and broken. In recent years, Candes and Donoho (2004) and Candes et al. (2006) developed a new geometric multiscale transform, the so-called curvelet transform (CT), which allows an optimal sparse representation of objects with C²-singularities. The needle-shaped elements of the transform own very high directional sensitivity and anisotropy. The transform represents edges and singularities along curves much more efficiently than traditional WT. Therefore, CT has good capability to analyze ocean wave fields. This paper proposes a new method to determine wave directions from nautical X-band radar images based on the CT. The procedure can be briefly described as follows: (1) a subimage is cut from an original nautical X-band image; (2) the subimage is decomposed at different scales, in different directions and at different positions by using CT, and curvelet coefficients are obtained; (3) at each scale, the mean values of the coefficient matrices in each direction are calculated; (4) the direction that has the largest mean value correspond to the direction of wave orientation; and (5) the 180 ambiguity is removed by the calculating cross-correlation coefficients (CCCs) between successive collected images. The effectiveness of the method is demonstrated by experiments.

This section gives a brief introduction to the second-generation CT (Candes and Donoho, 2004; Candes et al., 2006). It is implemented through a band-pass filter bank in the Fourier domain.

Let x=(x_1,x_2)\in R^2 be a variable in spatial domain, \omega =({\omega }_1, {\omega }_2)\inR² be a variable in frequency domain, and (r,\theta ) be the polar coordinates in frequency domain. Define W(r) and V(t) representing a pair of smooth, real-valued, and non-negative window functions in the frequency domain. They are supported on r\in (1/2,2) and t\in [-1,1], respectively. Usually, W(r) and V(t) are called radial window and angular window, respectively. They need satisfy the following admissibility condition:

These conditions can be satisfied by taking the scaled Meyer windows (Ma and Plonka, 2007; Candes and Donoho, 2004),

where \nu is a smooth function,

For j\ge 0 scale, a frequency window U_j can be defined in Fourier domain as

where \left[j/2\right] is the integer part of j/2. According to the support of W(r) and V(t), W(2^{-j}·)=[2^{j-1},2^{j+1}], V(2^{\left[j/2\right]}{}^{-j}·)=[-2^{-\left[j/2\right]}, 2^{-\left[j/2\right]}], the support of U_j(r,\theta ) is a polar “wedge”, which is highly anisotropic.

Define a function {\phi }_j(x) by means of its Fourier transform {\widehat{\phi }}_j(\omega )=U_j(\omega ). Let {\phi }_j be a mother curvelet. Then the system of curvelets can be indexed by three parameters: a scale 2^{-j}; an equispaced sequence of rotation angles {\theta }_{j,l}=2\pi l\cdot 2^{-\left[j/2\right]}, 0\le l\le 2^{-\left[j/2\right]}-1, ; and a position x_k^{(j,l)}=R_{{\theta }_{j,l}}^{-1}(k_1{2}^{-j},k_2{2}^{-\left[j/2\right]}{)}^{\text{T}}, (k_1,k_2)\in Z^2, where R_{{\theta }_{j,l}} denotes a rotation matrix,

The curvelets are defined as

{\phi }_{j,l,k} is of rapid decay away from a 2^{-j} by 2^{-j/2} rectangle with center x_k^{(j,l)}, and orientation {\theta }_{j,l} with respect to the vertical axis in x.

Figure 1 demonstrates two curvelets in spatial domain. The big signal demonstrates a curvelet of the 1st direction at the 3rd scale. The small signal demonstrates a curvelet of the 7th direction at the 4th scale. Furthermore, in this case, there are 32 and 32 directions at the 3rd and 4th scales, respectively.

The curvelet coefficient c\left(j,l,k\right) of a function f\in L^2(R^2) can be given by an inner product,

where j, l and k are the scale, direction, and spatial location parameters, respectively. The inverse CT can be given as

Two fast digital CT algorithms are presented by Candes et al. (2006). The wrapping-based algorithm is employed in this paper. In the wrapping approach, both the forward and inverse digital CT algorithms have the same computational cost of O(N^2\log N) for an N\times N image.

Let f(x, y) be a small image with size m×n, g(x, y) be a large image with size M\times N, here m\le M, n\le N. Let S_{x,y} be a subimage in g(x, y) with top-left corner coordinate (x, y) and size m×n,

Let \rho (x,y) be the CCC between f(x,y) and S_{x,y},

where \mathrm{cov}(f,S_{x,y}) is the covariance between f(x,y) and S_{x,y},

\overline{f} and \overline{S_{x,y}} are the mean value of f(x,y) and S_{x,y}, respectively. D_f and D_S are the variances of f(x,y) and S_{x,y}, respectively,

Usually, a larger CCC represents higher similarity between two images.

A calibration experiment was implemented on the northwestern coast of the Zhangzi Island in the northern area of the Yellow Sea of China from December 2008 to April 2009 (Cui et al., 2010). A nautical X-band radar was mounted on the edge of a cliff at a height approximate 35 m above the sea level. The rotation period of the antenna was 2.5 s, and 32 images were recorded in a single radar image sequence. Radar pulse length and pulse repetition frequency were 70 ns and 3 000 Hz, respectively. The resolution of the images was 3.75 m. The water depth of the images area was 35 m. At the same time, a pitch-roll buoy (SZF), developed by Ocean University of China, was deployed about 1 km off the coast. Wave parameters were collected once every 3 h. During each collection, the buoy recorded about 17 min of continuous data consisting of a time series of 2 048 points with a sampling interval of 0.5 s. The buoy could measure 16 directions with a directional resolution of 22.5°.

Additionally, radar images collected on the eastern coast of Haitan Island, Fujian Province, China were also used in this paper (Chen et al., 2014a, b). During the experiment, the same radar was mounted on the top of a hill by the sea. The height of the hill is about 40 m.

Figure 2 is a nautical X-band radar image collected at 13:58 on 19 September, 2010 on the eastern coast of Haitan Island. A subimage with size of 825 m×825 m is chosen and labeled as I1. Taking I1 as a representative example, the method to determine wave directions will be introduced in the following part.

First, Subimage I1 is decomposed at different scales, in different directions, and at different positions by CT. The curvelet coefficients are obtained. The decomposition is conducted at 5 scales, and in 64, 128, and 128 directions at the 2nd, 3rd, and 4th scales, respectively. The 1st and 5th scales do not have directionality. To be more specific, there are 64 curvelet coefficient matrices at the 2nd scale, each coefficient matrix corresponding to a specific direction. For a specific coefficient matrix, each element in the matrix is a curvelet coefficient that corresponds to a specific position. The 1st scale is the coarsest scale. From the 1st to the 5th scale, the coefficients of higher scales reflect the information of higher frequencies.

To exhibit the structure of the curvelet coefficients more clearly, the coefficients are displayed as several concentric-quadrate rings as shown in Fig. 3. The most centralized square corresponds to the curvelet coefficient matrix of the 1st scale, and each element in the matrix is a curvelet coefficient that corresponds to a specific position. From the center to the periphery, the 1st concentric-quadrate ring corresponds to the coefficient matrices of the 2nd scale. The Subimage I1 is decomposed in 64 directions at the 2nd scale. The concentric-quadrate ring is formed by 4 big rectangles, and each big rectangle is formed by 16 small rectangles, respectively. Because Fig. 3 has been downsized, the gaps between these small rectangles become blurred. The direction referred here means the orientation of the curvelet basis function. The 2nd concentric-quadrate ring corresponds to the coefficients of the 3rd scale. It is formed by 128 small rectangles. Judging from the brightness level of these small rectangles, one can find that the coefficients in several directions are larger than the coefficients of other directions.

Second, the wave information’s centralized scales are determined. From the 1st to 5th scales, curvelet coefficients of only one scale were kept, respectively, while the coefficients of other scales were set to 0. Then, the inverse CT was employed to the modified curvelet coefficients, respectively. Therefore, one can obtain five images, which were constructed by the curvelet coefficients of only one scale, say by the curvelet coefficients of only the 1st, 2nd, 3rd, 4th, and 5th scales, respectively. Next, the standard deviations of the five images are calculated, respectively. Ocean waves have obvious scale characteristic. Therefore, the wave energy will concentrate in the curvelet coefficients of several specific scales. According to the wave theory, a large standard deviation indicates large energy in a signal (Gan et al., 2007). As shown in Fig. 4, the image, which was constructed by the coefficients of the 2nd scale, has the biggest standard deviation. The image, which was constructed by the coefficients of the 3rd scale, has the second largest standard deviation. Therefore, one can determine that the wave information is mainly centralized in the 2nd and 3rd scales for Subimage I1.

Third, the wave orientation (with the 180 ambiguity) is determined. From the 2nd to 4th scales, the mean values of the coefficient matrices in each direction are calculated, respectively. For example, Subimage I1 is decomposed in 64 directions at the 2nd scale, and there are 64 coefficient matrices at the scale. The mean values of the 64 coefficient matrices are calculated, respectively. The rest scales are manipulated in the same way. The ocean wave exhibits obvious directional characteristic, and the CT adopts anisotropic basis function. Therefore, it can distinguish signals in different directions. As shown in Fig. 5, the mean values of some directions are obviously greater than the mean values of other directions. These directions, which have the maximum mean value, correspond to the wave orientation, and these directions, which have the minimum mean values, correspond to the direction perpendicular to the wave orientation. Logically, we should determine the wave orientation by using the curvelet coefficients of the 2nd scale. However, given the fact that the 3rd scale presents the second largest standard deviation and the directional resolution of the 3rd scale is better than the 2nd scale, the wave orientation will be determined by using the coefficients of the 3rd scale.

As shown in Fig. 5, the 43th and 107th directions present the biggest mean values at the 3rd scale. There are 128 directions at the scale. The curvelet basis function rotates 127 times at this scale. The angle step is 2.81° (360°/128). For the 1st direction, the curvelet corresponds to 46.4° (45°+2.812 5°/2), the normal corresponds to 136.4° (46.4°+90°), the corresponding wave direction is 316.4° (136.4°+180°). The 43th direction corresponds to 74.4° (316.4°+(43-1)×2.81°), the 107th direction 254.3° (316.4°+(107-1)×2.81°). Therefore, one can determine that the waves propagate toward 74.4° or 254.3°. There is a 180° ambiguity.

Fourth, the 180° ambiguity is removed, and the wave direction is determined by calculating CCCs between continuous collected images. Figure 6a is the same image as shown in Fig. 2. To exhibit it more clearly, only the lower-left part of the figure is given in Fig. 6a. Furthermore, Figs 6a to c are continuous collected radar images. A subarea with size 450 m×450 m is chosen in Fig. 6a and denoted with a red rectangle. For convenience, the subarea will be called SI1 in the following part. Next, CCCs are calculated between SI1 and subareas in Figs 6b to c, respectively. Given that the wave orientation has been determined, it is not necessary to compute the CCCs in whole coverage in Figs 6b to c. The computation will be constrained in small areas, which are denoted with blue dotted rectangles. For convenience, they are called as search areas in the following part. Next, CCCs are calculated between SI1 and subareas within the search areas in Figs 6b to c, respectively.

As shown in Fig. 7, the blue and red lines depict the CCCs between SI1 and subareas in Figs 6b to c, respectively. The black line depicts the corresponding central position of SI1 in the search areas. The maximum CCCs on the blue and red lines are marked as MA1 and MB1, respectively. As discussed above, the maximum CCCs indicate that the subareas have the most similar spatial structure with SI1. As shown in Fig. 7, MA1 appears on the left side of the black line, MB1 appears on the left side of MA1. Therefore, one can determine that the wave direction is 74.4° in this area.

Furthermore, according to the position of MA1, the corresponding subarea, which is used to calculate the CCC, can be determined in Fig. 6b. The subarea is denoted with a yellow rectangle in Fig. 6b. Similar to the procedure, the corresponding subareas, which are used to calculate the maximum CCCs, can be found in Fig. 6c and is also denoted with a yellow rectangle. Furthermore, the corresponding original position of Subimage SI1 is denoted with red dotted rectangles in Figs 6b to c, respectively. One can determine that the waves propagate from northeast to southwest.

Owing to the limitation of the actual situations, the synchronous buoy data have not been obtained. The estimated wave direction is depicted by a red arrow in Subimage I1 in Fig. 2. By reading the successive collected radar images, one can find that the retrieved wave direction is consistent with the direction of wave propagation.

Taking as another representative example, Fig. 8 is a nautical X-band radar image collected at 15:35 on 25 March, 2009 on the northwestern coast of the the Zhangzi Island in the Yellow Sea of China. A subimage with size of 675 m×675 m is chosen and labeled as I2. The processing procedure of Subimage I2 is similar to the above procedure.

It is decomposed at 5 scales, and in 64, 128, and 128 directions at the 2nd, 3rd and 4th scales, respectively. First, the standard deviations of the five images, which are constructed by the curvelet coefficients of the 1st, 2nd, 3rd, 4th and 5th scales, are calculated, respectively. As shown in Fig. 9, the 2nd scale presents the biggest standard deviation, the 3rd scale presents the second biggest standard deviation. As discussed above, the wave direction will be determined by the curvelet coefficients of the 3rd scale.

Second, the mean values of the curvelet coefficient matrices in each direction at the 2nd, 3rd, and 4th scales are calculated, respectively. As shown in Fig. 10, the 38th and 102nd directions at the 3rd scale present the biggest mean values. Therefore, the wave orientation can be determined as 60.4° or 240.4°.

Third, three continuous collected radar images are given in Fig. 11, a subimage with size 450 m×450 m is chosen and labeled as SI2 in Fig. 11a, and two search areas are chosen along the retrieved orientation of wave propagation and are denoted with blue dotted rectangles in Figs 11b and c, respectively. On the basis of the above work, the CCCs are calculated between SI2 and subimages in the search areas in Figs 11b and c, respectively. As shown in Fig. 12, the blue and red lines depict the CCCs between SI2 and subimages in Figs 11b and c, respectively. According to the positions of the maximum values MA2 and MB2, the wave direction can be determined as 240.4° in this area. As shown in Fig. 8, the estimated wave direction is depicted by a red arrow in Subimage I2. According to the synchronous buoy observations, the wave direction was 247°. The estimated wave direction is consists with the buoy collected wave direction.

The Nautical X-band radar images and the synchronous buoy data observed on the coast of the Zhangzi Island in the northern area of the Yellow Sea are used to validate the effectiveness of the proposed method. Northeast or southwest winds prevail in the experimental area. Therefore, ocean surface waves are mainly propagate from the two directions. Forty-nine radar image sequences collected from 5 March, 2009 to 25 April, 2009 were used in this part. The significant wave height distributes from 0.7 to 2.3 m for the 49 samples.

As shown in Fig. 13, the wave directions observed by the buoy are plotted as black circles, the retrieved wave directions are plotted as black stars. The corresponding circles and stars are connected by the black lines. Two quantitative measurements proposed by the Ocean Wave GmbH Company for evaluating the performance of the WaMos II are employed here, Bias is 14.4°, RMSE is 18.6°. The accuracy of this method can meet the engineering needs.

The method proposed in this paper can be automatically accomplished by machine itself when the subarea is chosen, and the computation is high efficient. The directional resolution is related to the number of decomposed direction of the wave information’s centralized scales. For example, the Subimage I1 is decomposed in 128 directions at the 3rd scale. Therefore, the directional resolution for this method is 2.81° (360°/128).

This paper proposes a method to determine wave directions from nautical X-band radar images. The CT is a multiscale and multi-directional analytical tool. It is very efficient in analyzing wave signal in radar images. A subimage is chosen in an original image and decomposed by the CT. By judging the wave information’s centralized scales and directions, the wave orientation can be determined. The 180 ambiguity is removed by calculating CCCs between successive collected images. Compared with spectra-based methods, the proposed method requires only a few radar images, and the wave parameters are determined by identifying the geometric stripes.

The authors acknowledge the help received from the Curvelet.org team.