Ocean Wind and Wave Measurements Using X-Band Marine Radar : A Comprehensive Review

Ocean wind and wave parameters can be measured by in-situ sensors such as anemometers and buoys. Since the 1980s, X-band marine radar has evolved as one of the remote sensing instruments for such purposes since its sea surface images contain considerable wind and wave information. The maturity and accuracy of X-band marine radar wind and wave measurements have already enabled relevant commercial products to be used in real-world applications. The goal of this paper is to provide a comprehensive review of the state of the art algorithms for ocean wind and wave information extraction from X-band marine radar data. Wind measurements are mainly based on the dependence of radar image intensities on wind direction and speed. Wave parameters can be obtained from radar-derived wave spectra or radar image textures for non-coherent radar and from surface radial velocity for coherent radar. In this review, the principles of the methodologies are described, the performances are compared, and the pros and cons are discussed. Specifically, recent developments for wind and wave measurements are highlighted. These include the mitigation of rain effects on wind measurements and wave height estimation without external calibrations. Finally, remaining challenges and future trends are discussed.


Introduction
Ocean wind and wave information is crucial for various on-and off-shore activities such as coastal construction, ship navigation, and marine resource development (e.g., mining and oil industries).Traditionally, wind and wave information has been obtained from anemometers and buoys, respectively.However, wind flow distortion due to the superstructure or the movement of the anemometer platform may affect wind measurements [1], and the employment of buoys may be limited by the surrounding environment [2].Complementing the point measurements of such in-situ sensors, remote-sensing instruments such as synthetic aperture radar (SAR) [3][4][5], high frequency (HF) radar [6][7][8], and X-band marine radar [9][10][11][12] have been developed as alternative approaches over the last few decades.Particularly, X-band marine radar has been rapidly developed as an ocean remote sensor since it can image both the spatial and temporal variations of the sea surface with high resolutions.Furthermore, X-band marine radar is widely installed on ships for navigation purposes.If reliable wind and wave measurements can be obtained from such radars, the costs associated with traditional in-situ sensors such as anemometers and buoys can be significantly reduced.Thus, land-based or ship-borne X-band marine radar has garnered wide attention.
X-band marine radar generally operates at grazing incidence.Radar backscatter from the sea surface is primarily generated by Bragg resonance [9] interactions between the radar-transmitted electromagnetic microwaves and the high frequency gravity-capillary waves which are influenced by the effects of gravity and surface tension.For X-band radar operating with horizontal-transmit-horizontal-receive (HH) polarization, wedge scattering [13] and small-scale breaking waves [14] also play significant roles in radar backscattering.
Since the sea surface roughness that causes the radar backscatter is induced by local winds, the dependence of the radar cross section (RCS) on wind speed [15] and the angle between radar look direction and wind direction [16] can be exploited to extract wind parameters from radar images [12].On the other hand, long sea surface waves, including wind seas and swell, can be observed in radar images since they modulate radar backscatter signals.The modulations are mainly associated with hydrodynamic effects induced by interactions between capillary waves and long waves, tilt modulation due to the varying effective incidence angles of the radar signals along the slopes of long waves, and shadowing due to the obstruction of radar signals caused by higher waves [17].The shadowing and tilt modulation are depicted in Figure 1a. Figure 1b is an example of radar image.Due to the modulations, wave parameters can be extracted from radar images [18][19][20].Due to its high spatial and temporal resolutions, X-band marine radar has been used in many other oceanic applications, which involve target tracking [21], oil spill detection [22], internal wave analysis [23], wave group retrieval [24], surface current determination [25][26][27][28], upwelling observation [29], tide estimation [30], bathymetry mapping [31,32], surface elevation approximation [33,34] and prediction [35,36], coastal erosion [37], and detection of marine mammals [38].In this paper, only the work related to sea surface wind and wave measurements is reviewed.The purpose of the review is to comprehensively compare the techniques developed for such measurements and to explore remaining challenges.Section 2 contains a review of the techniques for wind parameter extraction.The techniques for wave parameter extraction are reviewed in Section 3. Conclusions and an outlook for future trends appear in Section 4.

Wind Measurements
Due to the dependence of the radar cross section (RCS) on wind speed [15] and the angle between the radar look direction and the wind direction [16], X-band marine radar has been exploited to retrieve wind parameters [12].Wind direction has generally been determined based on the fact that the backscatter of an HH-polarized X-band radar operating at grazing incidence presents only one peak in the upwind direction [39].Wind speed has generally been determined by establishing a model relating wind speed to other parameters derived from radar images.The wind estimation methods can be classified as wind-streaks-based, image-intensity-variation-over-azimuth-based, and spectral-analysis-based.

Wind-Streaks-Based Techniques
Wind streaks are linear features first observed in SAR images for ocean surface wind direction estimation [40].They are mainly caused by local turbulence, foam, and surfactants on the ocean surface.Such wind streaks are also visible in the temporally-integrated X-band marine radar images and aligned with the mean wind direction.Thus, they can also be used for wind measurements using X-band marine radar.

Local Gradients Method and Neural Network
The local gradients method (LGM) was originally proposed to extract wind-induced streaks from SAR images [41].Dankert et al. [42][43][44] used LGM to determine wind direction as the orientation of wind-induced streaks from X-band marine radar images.Wind speed was estimated using a neural network (NN) that employed the mean normalized RCS (NRCS), radar-derived wind direction, and radar look direction as the input parameters, and wind speed as the output.
For wind direction estimation, the first step is to integrate a radar image sequence over time.Ocean surface wave signatures that show high variability in the time domain can then be removed.As a result, only static and quasi-static signatures of low frequencies, such as wind-induced streaks, remain visible in the integrated image.In the next step, the integrated image is subjected to repeated smoothing and sub-sampling, in which the image resolution is thrice decreased by a factor of 2, resulting in a so-called Gaussian image pyramid.Subsequently, local gradients are obtained from all the sub-sampled images with three different resolutions using the optimized Sobel operators [41], and the orientations of local wind streaks are determined to be those that are orthogonal to the local gradients.Then, wind direction is derived as the mean or the mode of all the determined local wind streak orientations.While a 180 • ambiguity exists in the derived wind direction, it can be removed using the following techniques.The first technique is based on the cross-correlation function (CCF) of two temporally-integrated radar image sub-sequences.The moving distance and direction of image patterns can be indicated by the location of the CCF peak, thus removing the 180 • ambiguity.The second technique is based on the cross-spectrum (CS) of the CCF.The movements of image harmonic features can be derived from the phases of the CS values for the corresponding wavenumbers of the harmonic components.Since the resulting motion directions of the image harmonic features are always within ±90 • of the downwind direction, the 180 • ambiguity can be removed.Another technique, based on the extraction of wind gusts visible in the image sequence [45], is discussed in Section 2.1.2.The scheme of the LGM for wind direction estimation is shown in Figure 2.For wind speed determination, each temporally-integrated radar image over the range of 600 m or 900 m to 2100 m is divided into range-azimuth bins of 300 m in range and 5 • in azimuth.The least inputs required for the NN involve the mean NRCS of the bins within ±15 • in both cross-wind directions, the radar-retrieved wind direction and the radar look direction.The mean NRCS in the cross-wind directions is used for two reasons.Firstly, the NRCS is more sensitive to wind speeds in the cross-wind directions.Secondly, the NRCS in the upwind and downwind directions may be affected by blockage and wind shadowing due to neighbouring platforms, respectively.Improvements can be achieved if other parameters related to the sea state and atmospheric conditions are also included as the inputs of the NN.The sea state parameters include the signal-to-noise ratio (SNR), which is proportional to the square of significant wave height, and the peak wave phase speed, which can be calculated from peak wave frequency and peak wavenumber.The atmospheric conditions measured by external sensors include the air-sea temperature difference and the relative humidity.
In [42][43][44], two sets of data were used to test the techniques.The first dataset was collected from February to June 2001 on the platform Ekofisk 2/4 k located in the North Sea around 200 km off the west coast of Norway with wind speed up to 18 m/s.The second dataset was collected from August 2003 to November 2004, on the research platform FINO-I in the German Bight of the southern North Sea with wind speed up to 16 m/s.The comparison between in-situ and radar-derived wind directions showed a correlation coefficient of 0.99 and a standard deviation of 12.8∼14.2• .For wind speed estimation, by including the atmospheric conditions (air-sea temperature difference and relative air humidity) as additional inputs, the standard deviation can be reduced by 0.12∼0.24m/s.Good results for wind measurements have been obtained using this method.However, it may be difficult to apply the LGM to ship-borne radar data because ship motion may lead to less reliable estimation of the streak orientations from the temporally-integrated radar images.A preprocessing technique, such as georeferencing, may be required for its application to ship-borne radar data.For the NN technique, various input parameters are needed and a training process must be conducted.

Optical Flow-Based Motion Estimation
Dankert et al. [45] proposed an optical flow motion estimation-based technique (OFM) to retrieve wind fields from X-band marine radar images.This technique is based on the extraction of the movement of wind gusts that are visible in filtered and integrated images.Wind speed can be directly derived from the wind gusts without calibration.
Firstly, to filter out the ocean surface wave signatures, the linear dispersion relationship for surface gravity waves is applied to the 3D wavenumber-frequency spectrum of an image sequence.Next, a temporal moving average over 27 images is applied to a sequence of 32 images to remove high variability in the time domain, resulting in 5 averaged images G( γ, t), with γ = (x, y) and t being space and time coordinates, respectively.For each averaged image, the same smoothing and sub-sampling process as described in Section 2.1.1 is implemented to obtained the so-called Gaussian image pyramid.Then, wind gusts become visible in the resulting images.Based on a common assumption on optical flow, the intensity change in the averaged image is caused only by the motion of wind patterns.Thus, the total time derivative associated with G( γ, t) equals zero, which results in [46] (∇G where ∇G and G t are the spatial gradient and partial time derivative of the image intensity, respectively, T denotes the transpose of a matrix, and f o = [ f x , f y ] T is the optical flow, with f x and f y being the components in the x and y directions, respectively.f o denotes the local wind velocity vector in this case.Under the assumption that the optical flow is constant within a local area and by using a local weighted least-squares technique [46], f o can be determined.In Equation ( 1), the first-order derivative of the image intensity with respect to space and time are calculated by optimized 3D Sobel operators [46].Due to the application of 3D convolution (optimized 3D Sobel operators), the local wind velocity vector cannot be determined from the first and last images of the image sequence and the borders of images.Finally, all available estimated local wind velocity vectors are smoothed, and the mean wind velocity vector is determined as the mode of the smoothed local wind velocity vectors.A flowchart of the OFM is shown in Figure 3.To test the technique, only 1332 data samples from the same dataset collected on the platform Ekofisk 2/4 k in the North Sea were used [45].The correlation coefficient and standard deviation for wind direction measurements were 0.96 and 32.1 • , respectively, as compared to reference data from in-situ sensors.For wind speed measurements, the correlation coefficient and standard deviation were 0.64 and 2.68 m/s, respectively.It was observed that increasing the number of images involved in the moving average may improve the results.Although the technique does not require a calibration process, the accuracy of the results is relatively low compared to other methods.Also, as for the LGM, a preprocessing step, such as georeferencing, is required for ship-borne radar applications.

Intensity-Variation-Over-Azimuth-Based Techniques
As mentioned earlier, the variation of radar backscatter intensity over azimuth can also be exploited to retrieve wind direction.Other than the NN, wind speed can also be retrieved using an explicit model.

A. Single-CF
In [47], Lund et al. proposed a single-CF-based method for retrieving both wind direction and speed.Wind directions are extracted using a harmonic function that is least-squares fitted to the radar backscatter intensity as a function of azimuth.Wind speeds are determined from the average radar backscatter intensity using an empirical third-order polynomial obtained through training.
For HH-polarized X-band radar images collected at grazing incidence, only one backscatter intensity peak exists in the upwind direction [39].This is the principle for wind direction estimation.First, radar backscatter intensities are averaged along each unobstructed azimuthal direction.Then, the average radar backscatter intensity σ wDir θ as a function of the azimuthal direction θ is least-squares fitted by a cosine square function where p 0 , p 1 , and p 2 are parameters determined through least-squares fitting.Wind direction is determined as the azimuthal direction corresponding to the peak of the fitted function.
With least-squares fitting, wind directions can be retrieved even if the upwind direction appears in the obstructed portions of the image.
Wind speeds are retrieved by exploiting the dependence of radar backscatter intensity on wind speed.The overall average radar backscatter intensity σ wSpd of one radar image is calculated as the integration of the fitted function in Equation (2) over the full azimuth.This is expressed as The fitted function in Equation ( 2) is used for calculating σ wSpd because it automatically fills the gaps in obstructed image portions.Then, a third-order polynomial relationship between σ wSpd and wind speed is assumed and trained using the wind speed data that are collected by external sensors such as anemometers.Based on the obtained third order polynomial, wind speed can be retrieved from σ wSpd for each individual radar image.The basic steps of the single-CF method for wind direction and speed retrieval are illustrated in Figure 4.The datasets used to test the method were collected from R/V Roger Revelle travelling in the Philippine Sea near Taiwan island during two storms from 6 August to 9 August and from 12 September to 17 September 2010, respectively [47].Test results show that the CF-based method works well for ship-borne X-band marine radar data, even if some portions of the radar image are obstructed by ship structures.Furthermore, only one parameter derived from radar data is needed to train the model function, significantly simplifying the wind speed estimation compared to the NN-based technique.

B. Dual-CF
Under low sea states, the overall image intensities are low, and this may lead to inaccurate wind estimation when using the CF-based method.Therefore, Liu et al. [48] modified the original CF-based method [47] with a dual-CF technique to improve wind direction and speed results under low sea state conditions.
In the dual-CF technique, an initial estimation of wind direction is first obtained by using the single-CF method.Then, to update the wind direction value, the radar data within ±60 • of this estimated wind direction are employed for a second CF process.The range of ±60 • is used because it is wide enough to include the actual wind direction and it is narrow enough to exclude low image intensities due to low sea states.For wind speed retrieval, negative values of the fitted function derived from the dual-CF technique are discarded from the integration in Equation (3) for the calculation of the average radar backscatter intensity σ wSpd .
To test the technique, two sets of data [48] were collected by two radars (Decca and Furuno) installed on the Canadian Navy research ship CFAV Quest during a sea trial in 2008 in the North Atlantic Ocean off the east coast of Canada while wind speed varied from 2 m/s to 16 m/s.Both the original CF-based method and the dual CF technique were applied to the two datasets.The dual-CF technique showed trivial improvements over the single-CF method with respect to the standard deviations of both wind direction and speed.

C. Two-Model-CF
In both the single-CF method and dual-CF technique, radar data collected in the presence of rain were excluded from wind estimation.This is because rain affects the radar backscatter and leads to less reliable wind estimation.Huang and Gill [11] proposed an algorithm to improve wind speed extraction from rain-contaminated radar images.Radar data contaminated by rain are first identified and extracted based on a parameter called zero-pixel percentage (ZPP).Then, the same third-order polynomial used for rain-free data in [47] are trained separately for the rain-contaminated data.
The algorithm was tested using the Decca radar data as mentioned in Section 2.2.1-B and [11].The comparison between in-situ and radar-derived wind speeds using the two-model algorithm showed a standard deviation of 1.7 m/s, which is improved by 1.2 m/s compared to the results derived by using only one model for both rain-free and rain-contaminated radar data.However, in spite of the improvement of the wind speed estimation accuracy under rain conditions, the practicality of the algorithm is compromised by an additional procedure for distinguishing rain-free and rain-contaminated data and an extra training model.

Intensity Level Selection (ILS)-Based Methods
A. Original ILS Vicen-Bueno et al. [49] proposed an ILS-based method for retrieving wind parameters from temporally-integrated and spatially-smoothed radar images.Wind direction is determined by searching for the maximum range for a selected intensity level along the azimuth, while wind speed is derived from the selected intensity level and its corresponding maximum range by training an empirical third-order polynomial geophysical model function.
In the first step, radar images are temporally integrated and spatially smoothed to remove high variabilities in the time and space domains, respectively.Then, a sequence of intensity levels is defined based on the image intensity range.For every azimuthal direction of each integrated and smoothed image, the maximum range with an intensity level equal to or greater than each defined intensity level is identified.Next, for each defined intensity level, the identified maximum ranges from all azimuthal directions are smoothed by a moving average over 5 • , and the lowest defined intensity level with all its associated smoothed maximum ranges being greater than an inner distance boundary is designated as the selected intensity level.Finally, the wind direction is determined as the azimuthal direction corresponding to the peak of the smoothed maximum ranges associated with the selected intensity level.After applying the above procedures to the first several integrated and smoothed images, rather than all the defined intensity levels, only the selected intensity level from the previous image and its upper and lower defined intensity levels are used for wind retrieval.This makes the algorithm more efficient.
Wind speed w Spd , which is estimated from the selected intensity level and the peak of smoothed maximum ranges associated with the selected intensity level, can be expressed as with where r p max is the peak of smoothed maximum ranges associated with the selected intensity level L S , α r is the conversion rate between r p max and w Spd , and p 0 , p 1 , p 2 , and p 3 are the parameters determined by least-squares fitting.For the training process, wind speed is measured by external sensors such as anemometers.After the training is completed, wind speed can be estimated from L S and r p max for each individual integrated and smoothed radar image.Figure 5 depicts the basic steps of the method.Datasets used to test this method were collected from the research platform FLIP at approximately 30 km northwest of Bodega Bay, California, USA, during four periods in June 2010 while wind speed varied from 4 m/s to 22 m/s [49].As compared to reference values from anemometer data, the standard deviations of wind direction and speed results were 11 • and 0.6 m/s, respectively, for the training data, and 14.3 • and 0.8 m/s, respectively, for the test data.

B. Modified ILS
In the original ILS method, the obstruction of the radar field of view (FOV) and the appearance of islands may lead to less accurate results.Thus, Liu et al. [48] modified the original ILS method in [49] with blockage and island recognitions and an additional CF process, resulting in more robust wind direction and speed results when blockages or islands exist in the radar field of view.
First of all, blockages in the radar FOV due to obstruction by ship structures or other unknown factors are recognized.Blocked azimuthal directions are identified as those with ZPPs greater than 20%, and data in those directions are excluded from the ILS.Then, islands in the radar FOV are recognized as intensity bumps in the radar images.The intensity bump in each azimuthal direction is identified by applying the MATLAB built-in function, findpeaks.Data in the azimuthal directions in which islands are identified are also excluded from the ILS.Furthermore, an outer distance boundary is set for the ILS due to possible noise in the very far range.All smoothed maximum ranges associated with the selected intensity level should be smaller than the outer distance boundary.Finally, the harmonic function ( 2) is used to least-squares fit the smoothed maximum ranges to fill the gaps due to the exclusion of data.The azimuthal direction corresponding to the peak of the fitted function is determined as the wind direction and is further utilized for the wind speed determination which follows the original ILS procedure.
The same datasets as mentioned in Section 2.2.1-B were used to test the method [48].After using the modified ILS method, the standard deviations of wind direction and speed were reduced by 4.9∼7 • and 0.1∼0.3m/s, respectively.C. Texture-Analysis-Incorporated Huang et al. [50] proposed a texture-analysis-incorporated (TAI) algorithm to extract wind parameters from rain-contaminated radar images.The TAI method first requires the identification of the rain-contaminated pixels in each azimuthal direction of radar images based on texture analysis.Then, only the data from the azimuthal directions in which the rain contamination is negligible are used to extract wind parameters by applying the modified ILS method.
The texture of a radar image is generated based on the calculation of intensity RMS difference.Next, the portion contaminated by rain is identified.If the number of pixels with intensities larger than a threshold (see [50]) is higher than or equal to 20 in one azimuthal direction in the texture image, that direction is considered as the one with negligible rain contamination and is retained.By this means, wave signatures can be retained in the image texture while rain clutter can be removed.
The Decca radar data as described in Section 2.2.1-B were also used for validating this method [50].For rain-contaminated data, the standard deviations of wind direction and speed results are 19.9 • and 2.0 m/s, respectively, by using the TAI method, and 34.4 • and 3.3 m/s, respectively, without incorporating TAI.

Probability Distribution Function
Chen et al. [51] proposed a probability distribution function (PDF)-based method to retrieve nearshore sea surface wind vectors.With an adjustment for the influence of land, wind direction is also retrieved based on the azimuthal dependence of radar image intensity levels.Wind speed is determined from a model based on the PDF of radar image intensity levels and significant wave height, H S .
The method for wind direction retrieval is similar to single-CF.The first difference is that a sine function instead of (2) is used.In addition, the azimuthal direction corresponding to the peak of the fitted function, θ peak , is not regarded as wind direction directly but adjusted by an angle ∆θ.This is because the appearance of two hills in the radar image, θ peak may not be the wind direction.An average ∆θ is empirically obtained using the anemometer reference [51].
For wind speed retrieval, the expectation σ exp of the radar image intensity levels is calculated from the PDF of those levels.Then, a relationship between wind speed w Spd and σ exp as well as H S is sought through fitting where p 0 , p 1 , and p 2 are parameters determined by the least-squares fitting process.In this model, the wind speed changes approximately linearly with σ exp under low sea state conditions, while the term H S mainly takes into account swell effects under high sea states.After the model is obtained, wind speed can be determined from σ exp and H S .Datasets used to test the method were collected on Haitan Island, China, from October to December 2010 [51].The data consisted of 535 radar image sequences and simultaneous anemometer and buoy measurements, and were divided into 10 groups based on corresponding wind speed values which ranged from 0 m/s to 20 m/s in increments of 2 m/s.One third of the data randomly selected from each interval were used for training the wind speed model, and the remainder were used for validating the model.An acceptable RMS difference of 26.2 • for wind direction estimation was obtained while wind speeds ranged from 5 m/s to 20 m/s.For the data with wind speed varying from 0 m/s to 20 m/s, the RMS difference of wind speed retrieval results was 1.37 m/s.It was also found from the result that H S should be included in Equation ( 6) under moderate to high wind conditions.

Spectral-Analysis-Based Techniques
The methods described in Sections 2.1 and 2.2 are based on the analysis of the spatial and temporal variations in the images.Another category of algorithms were designed mainly based on image spectral analysis.

Background Noise
Izquierdo and Soares [52] presented a method to retrieve wind speed from the background noise (BGN) of radar image sequences.
An image frequency-wavenumber spectrum is first obtained by applying the 3D discrete Fourier transform (DFT) to the radar image sequence.Then, spectral components corresponding to ocean waves and BGN are separated using the linear dispersion relationship for gravity waves.The total spectral energy of BGN, F BGN T , is used to estimate wind speed based on the trained linear model [52].
Datasets used to test the method were collected close to the Port of Sines in Portugal, from September, 2000 to December, 2001 [52].A radar image was divided into 36 directional sectors, each of which covers 10 • .The calibration was implemented for each directional sector.The results showed that wind speeds derived from the sector along the upwind direction agreed best with the wind speed references.However, the RMS difference of the result is large compared to other methods.Due to the small dataset, the method needs to be further evaluated using more data.

1D Spectral Analysis
Most of the aforementioned methods can only provide good wind measurements in the absence of rain.However, if the radar data is contaminated by rain, less reliable wind measurements generally result.Recently, Wang and Huang [53] found that when the wind speed is high (over about 8 m/s), obtaining wind direction using the CF-based method in [47] is not significantly affected by rain.However, the accuracy may substantially decrease if the radar images are acquired during rain events under low wind speed.Thus, a method to estimate wind direction based on the 1D spectral analysis (1D-SA) of radar backscatter in the wavenumber domain is proposed in [53].Later, Huang and Wang [54] extended the 1D SA-based method for wind speed estimation from both rain-free and rain-contaminated radar data.
In the 1D-SA algorithm, a DFT is first applied to the backscatter intensity sequence in each azimuth θ of a radar image to obtain its wavenumber spectrum E θ (k) with k being wavenumber.For rain-free and high-wind-speed rain-contaminated cases, the spectral value at zero wavenumber |E θ (0)| instead of σ wDir θ is used for fitting Equation (2).Wind direction is also determined as the azimuth corresponding to the peak of the fitted function.In low-wind-speed, rain-contaminated cases, an integral of spectral value |E θ (k)| over the wavenumber range [0.01, 0.2] is calculated for each azimuth direction θ as S θ .Then, S θ is used in Equation ( 2) for least-squares fitting to determine the wind direction.
In wind speed retrieval, for both rain-free and rain-contaminated radar images, the spectral integration of |E θ (k)| with respect to θ and k is calculated for each image as S T .Then, a logarithmic relationship between S T and wind speed w Spd is sought through curve-fitting where p 0 , p 1 and p 2 are parameters determined from the curve-fitting process.After the parameters are determined based on (7), wind speed can be estimated from S T for any radar data.The Decca radar data as described in Section 2.2.1-B were used to test the method [53,54].Both the single-CF and 1D-SA methods were applied.The low-backscatter images were excluded for the wind direction retrieval but retained for the wind speed retrieval.The test results indicate that the 1D-SA wind algorithm is applicable to both rain-free and rain-contaminated radar data.Moreover, the method significantly improves wind measurements under rain conditions with a reduction of 25.1 • and 5.9 m/s, respectively, in the standard deviation of wind direction and speed over the single-CF.However, the wind direction result may not be accurate if the upwind direction is contaminated by rain, and the wind speed result is not always robust when the wind speed is high.

Ensemble Empirical Mode Decomposition-Based Algorithms
Ensemble empirical mode decomposition (EEMD) is a technique that is developed for the analysis of non-stationary and non-linear data.By using the EEMD, such data can be decomposed into several intrinsic mode functions (IMFs) which are narrow-banded in the frequency domain.Then, the analysis of the original data can be realized through the analysis of each IMF.
A. 1D-and 2D-EEMD Liu et al. [55] proposed 1D-and 2D-EEMD-based algorithms to retrieve wind direction from rain-contaminated X-band marine radar data.In the 1D-EEMD-based algorithm, assuming that one polar radar image with M rows and N columns is expressed as where i m,n is the image intensity at the mth row and the nth column.The application of the EEMD [56] to each column of I results in where g m,n,j is the element of G j , G j with j ranging from 1 to J − 1 is the jth IMF component of the 1D EEMD of I, in which J − 1 is the number of the IMF components, and G J is the residual component.
Then, the standard deviation of one IMF component, or the combination of several IMF components, as a function the azimuthal direction is least-squares fitted to the harmonic function (Equation ( 2)).The azimuthal direction corresponding to the peak of the fitted function is retrieved as the wind direction.
In the 2D-EEMD-based algorithm, the EEMD [56] is again applied to each row of G j , which results in where h m,n,j,k is the element of H j,k , H j,k with k ranging from 1 to K − 1 is the kth IMF component of 1D EEMD of G j , in which K − 1 is the number of the IMF components, and H j,K is the residual component.Then, based on a comparable minimal scale combination principle [57], the lth IMF of the 2D EEMD of I can be obtained as where l ranges from 1 to the minimum of J − 1 and K − 1, and C min(J,K) is the residual component.
As in the 1D-EEMD-based algorithm, wind direction is determined as the azimuthal direction corresponding to the peak of the harmonic function (Equation ( 2)) that is least-squares fitted to the standard deviation of one IMF component or the combination of several IMF components.
Only rain-contaminated radar data collected from the Decca radar as described in Section 2.2.1-B were used to test the 1D-and 2D-EEMD-based algorithms [55].Both of the EEMD-based algorithms were compared to the earlier 1D-SA-based algorithm.With reference to anemometer measurements, the correlation coefficients of the wind directions retrieved by both the EEMD-based algorithms were 0.98 while that retrieved by the 1D-SA-based algorithm was 0.95.The RMS differences for the 1D-EEMD-, 2D-EEMD-, and 1D-SA-based algorithms were 12.7 • , 11.4 • , and 20.1 • , respectively.Better results were obtained by the EEMD-based algorithms.However, the computational cost, especially for the 2D-EEMD-based algorithm, is high.

B. EEMD-Normalization
Later, Huang et al. [58] proposed an EEMD-normalization-based algorithm to extract wind direction as well as wind speed from both rain-free and rain-contaminated X-band marine radar data.
First, a data control strategy was proposed to distinguish rain-free, low-wind-speed rain-contaminated, and high-wind-speed rain-contaminated radar data.For each azimuthal direction of one radar image, if the intensity levels of all pixels are higher than 1, the direction is marked as a high-clutter direction (HCD).If the HCD percentage (HCDP), which is defined as the ratio of the number of HCDs to the number of total directions for one radar image, is higher than 5%, the image is classified as a rain-contaminated radar image.Otherwise, it is categorized as a rain-free radar image.For a rain-contaminated image, if the high-pixel percentage (HPP) [53], which is defined as the ratio of the number of pixels with intensity levels higher than 100 to the total number of pixels, is higher than 30%, the image is classified as a high-wind-speed rain-contaminated image.Otherwise, it is categorized as a low-wind-speed rain-contaminated image.Then, after the application of the 1D EEMD to the polar radar image I, a normalization scheme [59] is applied to each column of G 1 , which is the first IMF component of the 1D EEMD of I.With the normalization scheme [59], each column of G 1 is separated into its amplitude modulation (AM) and frequency modulation (FM) parts, which is expressed as where • denotes the element-by-element multiplication of two vectors, and a 1 (∼, n) and f 1 (∼, n) are the AM and FM parts of G 1 (∼, n), respectively.For rain-free and high-wind-speed rain-contaminated radar images, the average of the residual component G J , as a function of the azimuthal direction, is least-squares fitted to the harmonic function (Equation ( 2)).Wind direction is determined as the azimuthal direction corresponding to the peak of the fitted function.For low-wind-speed rain-contaminated radar images, the average of the AM part of the first IMF, a 1 , as a function of the azimuthal direction, is least-squares fitted to the harmonic function (Equation ( 2)).The azimuthal direction corresponding to the peak of the fitted function is determined as the wind direction.For wind speed retrieval, the total averages of a 1 and G J for each radar image are calculated as m a and m r , respectively, and a combination M com is obtained as Then, Equation ( 7) is trained to seek the relationship between wind speed and M com through curve-fitting.After the training process, wind speed can be retrieved from M com for any radar image using Equation (7). Figure 6 gives the basic steps of the EEMD-normalization-based algorithm for wind parameter estimation.Two datasets collected from a sea trial off the east coast of Canada from November 26 to December 4 using Decca and Furuno radars were used to test the EEMD-normalization-based algorithm [58].During the test, anemometer-measured wind speed varied from 2 m/s to 16 m/s.The 1D-SA-based algorithm was also applied to the datasets.For wind direction retrieval, the RMS differences compared to anemometer measurements were 11.5 • and 12.3 • for the EEMD-normalization-based and 1D-SA-based algorithms, respectively, for all data, and 12.9 • and 17.8 • , respectively, for rain-contaminated data.The RMS difference for the 1D-SA-based algorithm is different from that in Sections 2.3.2 and 2.3.3-Asince more data are included in this section.For wind speed estimation, the RMS differences were 1.31 m/s and 1.48 m/s for the EEMD-normalization-based and 1D-SA-based algorithms, respectively, for all data, and 1.41 m/s and 1.92 m/s, respectively, for rain-contaminated data.Improvements for wind parameter extraction, especially in rain cases, were obtained by the EEMD-normalization-based algorithm.It should be pointed out that, compared to the 1D-SA-based algorithm, the EEMD-based algorithms significantly improve the results for cases in which the image portions along the upwind direction are also contaminated by rain.The limitation of the EEMD-based algorithms is that the computational cost is high.

Summary
The inter-comparison of the reviewed techniques for wind measurements is listed in Table 1.It should be noted that the RMS difference values in Table 1 were obtained from different datasets.Then, it may not be appropriate to compare the performance of different methods based only on these values.It can be observed from the table that ship motion has a significant effect on the LGM and OFM; the original ILS and BGN methods may not be applicable if portions of the radar images are blocked; the two-model-CF, texture-analysis-incorporated, and 1D spectral analysis algorithms work well under rain effects; and the OFM is calibration-free for wind speed retrieval.

Wave Measurements
Wave information can be extracted from X-band marine radar sea surface images due to the visibility of ocean waves in the radar images.Algorithms for such wave measurements are usually based on a wavenumber-frequency spectral analysis of the radar image time series using a 3D DFT [9].Recently, a variety of algorithms for wave measurements using X-band marine radar have appeared.These algorithms can be categorized as spectral-analysis-based and texture-analysis-based for non-coherent radar.For coherent radar, the algorithms are mainly based on the estimation of surface radial velocity.

Spectral-Analysis-Based Techniques
For wave algorithms based on spectral analysis, wave spectra are first obtained from series of radar images.Then, wave parameters are derived from the obtained wave spectra.

A. Traditional Algorithm
The traditional 3D-DFT-based algorithm is a relatively mature technique, which has also been implemented in the commercial wave monitoring system WaMoS II [2,10].
In this technique, rectangular sub-image sequences are first extracted from the time series of consecutive radar images.Then, the intensity of each image pixel in the sub-image sequence is normalized by subtracting the temporally-averaged intensity of the corresponding pixel.This image normalization enables the mitigation of static patterns caused by radar backscatter dependence on range.After that, a 3D DFT is applied to the normalized sub-image sequence to produce an image spectrum.A high-pass filter with an angular frequency threshold ω th is then implemented on the image spectrum to eliminate non-stationary and non-homogeneous trends in the sub-image sequence, resulting in a high-pass filtered image spectrum F (3) where k = (k x , k y ) is the wavenumber vector and ω is the angular frequency.Subsequently, an encounter current velocity, defined as the combined velocity of the radar platform and sea surface current, is estimated from the high-pass filtered image spectrum F (3) I ( k, ω) using appropriate methods [25][26][27][28].The velocity is then incorporated into the linear wave dispersion relationship to extract the wave-related image spectrum F Usually, , where β 0 is the MTF exponent determined empirically using the 1D wavenumber spectra derived from radar measurements and in-situ sensors [17].With the obtained 3D wave spectrum F W ( k), 2D frequency-directional wave spectrum E (2) (ω, φ), and 1D frequency wave spectrum S(ω) can be derived [10,18].Then, wave parameters such as wave direction and periods can be deduced from the wave spectra [19,20].Figure 7 shows the basic steps of the traditional 3D-DFT-based algorithm.However, since the radar image intensities are not directly related to sea surface elevations, but to the strength of radar backscatter signals, the wave spectra obtained as above represent relative values which are not properly scaled.Therefore, wave height cannot be directly determined from the radar-derived wave spectra.For estimating wave height, a widely accepted approach is based on a technique originally developed for SAR applications [60], in which H S is assumed to be proportional to the square root of the signal-to-noise ratio (SNR) [61].The signal is assumed as the total energy of the estimated wave spectrum and the noise is computed as the energy due to the speckle caused by the sea surface roughness.The 3D BGN spectrum F where F (3) I ( k, ω) by employing the high harmonic linear wave dispersion relationship.Then, the SNR may be computed from the wave and BGN spectra [61].Finally, a linear model relating H S to the square root of SNR can be constructed as where p 0 and p 1 are parameters determined through fitting the square root of SNR and buoy-measured H S .After the linear model is obtained, H S can be derived from the SNR for each radar image sequence.
The SNR-based method has been tested by many researchers with various datasets [2,10,19,20,[61][62][63].These test results showed that the correlation coefficients between H S measured by buoy and radar were 0.71∼0.89,and the RMS differences were 0.18∼0.42m.In [63], it was found that the difference of buoy and radar positions mainly accounted for the low correlation of 0.71 due to a significant change of bathymetry.In [20], the method was applied to the data collected by HH-and vertical-transmit-vertical-receive (VV)-polarized radars.Comparison showed that the VV-polarized radar was more suitable for the H S estimation since it produced better results.Although this SNR-based method for H S estimation has been successfully applied to such stationary radar systems, it was found in [64,65] that it may be challenging to achieve a comprehensive calibration for radar operations on moving vessels on the open sea.
Later, in [66], the dependency of wave parameter retrieval on range and azimuth was analyzed and removed.This enables robust wave parameter measurements in the case where the radar FOV is partially obstructed.By retrieving the SNR, peak wave period, and peak wave direction from sub-image time series azimuthal-evenly distributed in the radar FOV, it is found that the SNR shows a higher peak in the upwave direction and a lower peak in the downwave direction, and both the peak wave periods and directions derived from different azimuths display a period of 180 • .Thus, for the H S estimation, the SNRs derived from the azimuthal-evenly distributed sub-images in the unobstructed radar FOV are least-squares fitted to a Fourier series.Then, the higher peak (upwave direction) or the lower peak (downwave direction) of the fitted curve is used to estimate H S .The peak wave period and direction are estimated as the average of those derived from the azimuthal-evenly distributed sub-images covering 180 • in the unobstructed radar FOV.

B. Multilayer Perceptrons
A limitation of the standard SNR-based method for H S estimation is that local wind speed should be higher than a threshold to produce sufficient sea surface roughness and radar backscatter.Under low wind speeds or swell-dominated sea state conditions, the BGN energy may be low, leading to large SNR and, therefore, overestimated H S .Thus, Vicen-Bueno et al. [67] proposed a non-linear model using artificial neural networks (ANN) with multilayer perceptrons (MLP) for H S estimation.
In addition to the square root of the SNR, peak wavelength and mean wave period are also used as input parameters.One hidden layer with 15 hidden neurons is selected for the MLP.
The technique was tested using two datasets [67].One dataset with 49096 measurements and manifested bimodal (wind sea and swell) sea states was collected at Ekofisk in the Norwegian sector of the North Sea from October, 2008 to May, 2009.Another dataset with 67,684 measurements was collected at FINO 1 in the German basin of the North Sea, with a swell-dominated sea state, from April, 2005 to August, 2008.Both the standard SNR-and MPL-based methods were applied to the two datasets.As compared to the reference values measured by buoy, the correlation coefficient and the standard deviation of H S derived by the MLP-based method were a little better than those of the standard SNR-based method.A more significant improvement was observed for the dataset collected at FINO 1 since this location was dominated by swell.

C. Iterative Least-Squares
Huang et al. [68] presented an iterative least-squares (LS)-based method for wave measurements, which is based on the iterative LS methods for ocean surface current determination from X-band marine radar images [25,69].This method, which does not require a band-pass filter as in the traditional 3D-DFT-based algorithm (see Figure 7), simplifies the wave measurement.
After applying the 3D DFT to the radar image series and high-pass filtering the image spectrum, the ocean surface current has to be determined.An initial guess of the surface current is first derived using a basic LS method as in [9].Then, a lower threshold than that in [9] is used to obtain spectral points that may correspond to the fundamental and harmonic modes of the linear wave dispersion relationship.The current-induced Doppler-shifted dispersion relationship at the pth harmonic mode is expressed as where g is the acceleration due to gravity, and U is the surface current vector.In each iteration, ω p ( k) is calculated using the surface current vector derived from the previous iteration or the initial guess.Due to the aliasing effect, ω p ( k) located in the qth frequency interval is folded to ω p,q ( k) through mapping according to the properties of periodicity and point symmetry about the origin.The obtained spectral points that satisfy min p,q where ω i is the angular frequency in the ith angular frequency panel of the 3D image spectrum and ∆ω is the angular frequency resolution, are used to estimate a new surface current vector by the LS fitting as in [25].The iteration is repeated until the termination criterion defined in [69] is met.When the final surface current is determined, since the spectral points used in the iterations have been automatically classified as the fundamental and different harmonic modes of the dispersion relationship, the band-pass filter in the traditional 3D-DFT-based algorithm is not required.Here, the spectral points associated with the fundamental mode (p = 0) and the first harmonic mode (p = 1) are extracted for the wave measurement.Then, after applying an MTF, wave spectra and related wave parameters can be derived.The iterative LS-based wave measurement method was tested using both vertically-and horizontally-polarized X-band radar data [68].Both the traditional 3D-DFT-based algorithm and the iterative LS-based method were applied to the data.The vertically-polarized land-based radar data were acquired at Skerries Bight near St. John's Harbour on the east coast of Canada on 15 and 20 December 2010.The results showed that compared to the buoy reference, the RMS differences of the peak wave periods and peak wave directions derived using both techniques were less than 1 s and 10 • , respectively.The horizontally-polarized ship-borne radar data were acquired about 220 km Southeast of Halifax from 14:16 on 28 November to 12:06 on 29 November 2008, with peak wave period ranging from 10 s to 12 s and peak wave direction being around 110 • .The RMS differences for the peak wave periods derived from the traditional algorithm and the iterative LS-based method were 0.87 s and 0.98 s, respectively.For the peak wave directions, the corresponding values were 6.4 • and 7.3 • , respectively.The results derived from the iterative LS-based method were similar to those from the traditional 3D-DFT-based algorithm, which indicated the capability of the method.

D. New MTFs
In the traditional 3D-DFT-based algorithm, due to radar imaging mechanisms such as tilt modulation and shadowing, a linear MTF with a constant exponent is usually required for wave spectrum correction [17].Such an MTF is developed using an HH-polarized radar under deep water conditions, where the wave field can be safely assumed as stationary and homogeneous without severe interactions of wind sea and swell, and the variation of the exponent due to range and azimuth can be neglected.However, under shallow water conditions like near-shore regions, the linear MTF with a constant exponent may not work well due to the heterogeneity of the wave field.Therefore, Chen et al. [70] proposed a new quadratic polynomial MTF for near-shore regions using a VV-polarized radar.The MTF is fitted by a quadratic polynomial as where p 0 , p 1 , and p 2 are coefficients determined by least-squares fitting.Later, Qiu et al. [71] proposed another MTF with the same linear model as in [17], but considering the dependence of the exponent on range and azimuth.The range-and-azimuth-dependent variable exponent of the MTF, β v , is given by where r a is the range, and θ v the angle between the radar look direction and peak wave direction.The methods incorporating the new MTFs were tested using datasets collected on Haitan Island, China, from 22 December 2012 to 6 January 2013 and from 12 to 18 January 2013.Peak and mean wave periods were derived using both of the new MTFs and the traditional constant-exponent linear MTF.The comparison between the quadratic polynomial MTF and the linear MTF with the exponent of -0.55 involved 520 measurements [70].Compared to the buoy measurements, the RMS difference of peak wave periods derived using the quadratic polynomial MTF was 0.88 s, which was similar to that of 0.74 s derived with constant-exponent linear MTF.The RMS difference of mean wave periods derived using the quadratic polynomial MTF was 0.53 s, which was 0.50 s better than that using constant-exponent linear MTF.For the comparison between the variable-exponent linear MTF and the linear MTF with the exponent of -1.2,only the data with H S higher than 2 m were used [71].During the test, the buoy-measured peak and mean wave periods varied from 7 s to 11 s and from 6.5 s to 9 s, respectively.The RMS differences of the peak and mean wave periods for the variable-exponent linear MTF were 0.95 s and 0.48 s, respectively, which were improved by 2.24 s and 2.21 s, respectively, over the constant-exponent linear MTF.The results indicated that the new MTFs may be more suitable for near-shore regions where the ocean wave field is heterogeneous.

E. Adaptive Recursive Positioning
In the traditional 3D-DFT-based algorithm for wave measurements, since radar backscatter is dependent on the azimuthal direction, sub-image sequences located in different azimuthal directions may lead to different wave measurement results.This problem has been mitigated by averaging wave spectra derived from several sub-image time series evenly distributed in azimuth in the radar FOV [72].
Recently, Al-Habashneh et al. [73] proposed an adaptive recursive positioning method (ARPM) to further improve wave measurements, especially under bimodal (wind wave and swell) sea states.
In the first step, an initial directional wave spectrum is derived using three sub-image time series evenly distributed in azimuth in the radar FOV as in the traditional algorithm.Then, azimuthal directions of spectral peaks can be determined from the directional wave spectrum.In the next iteration, sub-image time series located in those azimuthal directions determined from the last iteration are used to derive a new directional wave spectrum.Final wave spectra and related wave parameters are determined after three iterations.The flowchart of the ARPM is illustrated in Figure 8.The ARPM in [73] was tested using datasets collected from ship-borne Decca and Furuno X-band marine radars near Halifax, Nova Scotia, off East Coast of Canada, on 27 November, 29 November, 1 December, and 3 December 2008.The performance was validated by calculating correlation coefficients between radar-derived and buoy-measured wave spectra on a discrete frequency grid.Compared to the traditional 3D-DFT-based algorithm, the ARPM produced an improvement of 9.8% in the averaged correlation coefficient.The accuracies for wave period and peak wave direction estimations were also enhanced by 15% to 30% and 6 • , respectively.Although, due to the iterations, additional computational time is needed for the ARPM, real-time applications can still be realized.A limitation of the ARPM is that its best performance requires a full radar FOV.If the radar FOV is partially obstructed, it may not improve the wave measurement over the traditional 3D-DFT-based algorithm.

F. Geometrics-Based SNR Estimation
The standard SNR-based H S estimation requires the surface current velocity to be determined in order to separate the wave components from the noise in the image spectrum [42].However, possible inaccuracies in the determination of surface current may lead to erroneous SNR, thus affecting the H S estimation.Wang et al. [74] proposed an algorithm to derive the SNR based on the geometrics of the linear wave dispersion relationship, which does not require the surface current determination.
A 3D image spectrum is first obtained by the application of a 3D DFT to a sub-image time series.With the assumption of deep water condition, for any given point ( k 1 , ω 1 ), a set of curves can be obtained based on the geometric relationship as An example of these curves can be found in [75].Of all the curves in each ω − k plane, those satisfying are first extracted.G c (m c , n c ) is the m c th spectral value on the n c th curve, which is obtained using the nearest-neighbour interpolation.Then, the extracted curves satisfying where N c is the total number of the curves extracted using Equation (22), are determined to be associated with the dispersion relationship.Other curves are associated with noise.Therefore, the signal and noise can be separated for each ω − k plane, and the SNR can be derived for the H S estimation.
The geometrics-based SNR estimation algorithm was tested with datasets collected at Zhoushan City, Zhejiang Province, China, from January 6 to 12, 2014, with H S varying from 0.5 m to 2 m [74].Both the geometrics-based algorithm and standard SNR-based method were applied to the datasets.As compared to the buoy reference, the correlation coefficient and the RMS difference of the retrieved H S using the geometrics-based algorithm were 0.89 and 0.19 m, respectively, which were improved by 0.06 and 0.03 m, respectively, over the standard SNR-based method.

A. Original Algorithm
For the wave extraction using the 3D-DFT-based algorithms mentioned above, the homogeneity within the observed area is assumed.Such an assumption is safe under deep water conditions without severe interactions of wind sea and swell.However, in coastal areas with shallow water, heterogeneities may exist because of the varying topography.Therefore, Chuang et al. [76] and Wu et al. [77] applied a 2D continuous wavelet transform (CWT) to X-band marine radar images to obtain ocean wave spectra.
The 2D CWT of an image s( γ), where γ = (x, y) represents the coordinates in the image, may be expressed as [78] where C Ψ is a constant denoting that the mother wavelet function Ψ satisfies the admissibility condition, b = (b x , b y ) the translation parameter, a the dilation parameter, and r −ϕ the rotation matrix, which is defined as In [77], the Morlet wavelet is selected as the mother wavelet function because of its directionality.To reduce the computation time, the CWT can be calculated in discrete form and in the Fourier space [78].After dilation of wavenumber and rotation properly, the local wavenumber spectrum at position (b x , b y ) can be obtained.X-band radar images used to verify the 2D-CWT-based algorithm were collected in the southern part of Taiwan [77].The directional wave spectra derived by the 2D-CWT-based algorithm showed local ocean wave features at different locations of radar images.Furthermore, compared to the Fourier transform, the normalized 1D wavenumber spectrum derived using the 2D-CWT-based algorithm agreed better with that measured by an in-situ buoy.

B. Self-Adaptive 2D-CWT
In the original 2D-CWT-based algorithm, the effect of ocean wave conditions on the selection of wavelet parameters is not considered.It is shown in [79] that the dilation parameter may affect wave extraction, and a look-up table is constructed for selecting the dilation parameter values according to different wave conditions.However, the wave data collected by external sensors are used to make such a table.Later, An et al. [80] proposed a self-adaptive 2D-CWT-based algorithm for wave extraction by selecting the dilation parameter adaptively.
The maximum value of a is regulated by a calibration parameter β.A look-up table for determining β with different mean wave periods T 01 is provided in [79].The values of β in the table are based on an image resolution ∆x of 10.5 m and a sub-image length N x of 128.For a different image resolution ∆x n and a different sub-image length N x n , β n may be modified as In the algorithm for selecting appropriate β, a sub-image is first selected and normalized.Then, an image spectrum is obtained by applying the 2D CWT with an initial value of β (e.g., 1.4 for ∆x = 10.5 and N x = 128).Next, the 1D frequency wave spectrum E( f ) and the mean wave period T 01 are derived.If any peak above a threshold exists at the left side of the major peak of E( f ), it is regarded as a false peak.If a false peak exists, the value of β is increased by ∆β and the 2D CWT is applied to the normalized sub-image again.If no false peak exists, the calculated T 01 is compared with that in the look-up table for the corresponding β.If the difference between the two periods is larger than a threshold ξ (if T 01 < 10 s, ξ = 0.5; otherwise, ξ = 1), the value of β is increased by ∆β and the 2D CWT is applied to the normalized sub-image again.If the difference is not larger than ξ, β is selected as the appropriate value for the wave extraction using the 2D-CWT-based algorithm.The flowchart of the self-adaptive 2D-CWT-based algorithm is shown in Figure 9.
Both the traditional 3D-DFT-based and self-adaptive 2D-CWT-based algorithms were applied to the datasets that are described in Section 3.1.1-C[80].The results from both the land-based and ship-borne radar data showed that the wave spectra and wave parameters derived from both algorithms agreed well with buoy records.Furthermore, for the ship-borne radar data, taking the buoy data as ground truth, the RMS differences of H S , mean wave period, and wave direction derived from the traditional 3D-DFT-based and self-adaptive 2D-CWT-based algorithms were 0.53 m and 0.61 m, 1.91 s and 1.51 s, and 22.96 • and 24.18 • , respectively.During this experiment, the buoy-measured H S , mean wave period, and wave direction were about 2 m to 3.5 m, 7 s to 10 s, and 100 • to 120 • , respectively.Although the results indicated the applicability of the self-adaptive 2D-CWT-based algorithm, some limitations of the algorithm still exist.Firstly, the computational cost of the algorithm is expensive.Secondly, the effect of surface current is not considered.Thirdly, the directional ambiguity cannot be eliminated by the 2D CWT analysis of a single radar image.

Array Beamforming Algorithm
Ma et al. [81] proposed an array beamforming algorithm to extract the directional wave spectrum and related wave parameters from X-band radar image time series.Beamforming is a technique that extracts the signal arriving from a desired direction through constructive interference and suppresses other signals through destructive interference [82].In this algorithm, image pixels are considered as elements in an antenna array, and image pixel intensities are considered as signals received by the elements.
Assuming that X ω (t) is an A × 1 array, representing the signals with the angular frequency ω and different directions received by an antenna array of A antennas, the beamforming result for direction φ may be expressed as where H is the conjugate transpose operator, and a ω,φ is the steering vector, which consists of the time delays of the received signal between the reference antenna and other antennas.For X-band radar sea surface imaging, s ω,φ (t) may represent the sea surface waves with the amplitude u ω,φ , and X ω (t) may represent the complex image pixel intensities.A circular antenna array is selected for this algorithm, and by considering the radar image pixel as the antenna, an approximate circular array is constructed by the pixels located on a circle.X ω (t) can be obtained by applying an all-phase fast Fourier transform (ApFFT), which significantly improves the phase spectrum accuracy [83], to image intensity time series at each pixel.The time delays in a ω,φ are determined based on the wavenumber vector and the geographical positions of the reference pixel and other pixels in the selected circular array.The wavenumber vector is derived from the selected angular frequency ω and direction φ using the dispersion relationship with current effects.Thus, different current velocity values lead to different beamforming results and, in turn, result in different directional wave spectrum E (2) (ω, φ), which is where ∆ω and ∆φ are the angular frequency and direction resolutions, respectively.Current velocity is retrieved by seeking the maximum total power, which is the integration of E (2) (ω, φ).Then, wave parameters such as wave directions and periods can be derived from the directional wave spectrum E (2) (ω, φ).H S is estimated based on the SNR using Equation ( 16), where the signal energy can be calculated from E (2) (ω, φ), and the noise energy is estimated as the subtraction of the signal energy from the total energy, which is derived from the mean of the power spectrum at each pixel.Figure 10 illustrates the basic steps of the array beamforming algorithm.An experiment was conducted in Zhoushan, China, on January 8, 9, 10, and 12, 2014 [81].During the experiment, buoy-measured H S , mean wave period, and peak wave direction ranged from 0.5 m to 2 m, 4 s to 7 s, and 50 • to 250 • , respectively.For 45 buoy-radar comparison pairs, the correlation coefficient and RMS difference of the mean wave period, peak wave direction, and H S are 0.95 and 0.28 s, 0.9 and 15.14 • , and 0.95 and 0.12 m, respectively.Moreover, the wave spectrum derived by the array beamforming algorithm is similar to that derived by the traditional DFT-based method with an MTF.The novel array beamforming algorithm utilizes the beamforming technique based on individual image pixels, suppressing the effects of shadowing modulation and moving vessels in the images.The coordinate transformation and MTF can also be avoided.In the future, the validation of the algorithm for ship-borne radar applications with higher encounter current speed and for heterogeneous wave fields may be required.

Texture-Analysis-Based Techniques
Another class of algorithms for wave measurements is based on the texture analysis of radar sea surface images.Wave parameters, such as H S , are derived directly from the image texture without the estimation of wave spectra.

Illumination Probability
A model that relates H S to the crest-to-trough length ratio extracted from the radar image and a threshold probability of illumination (P 0 ) was first proposed in [84] for estimating H S .Crests and troughs in the radar image represent illuminated areas reflected from visible wave portions and shadowed areas with no reflections, respectively.In [85], Buckley and Aler validated this model by assuming a constant P 0 based on the theory of geometric optics [14].The essence of the theory is that, for a constant P 0 , the crest-to-trough ratio decreases as the wave height increases.Later, Buckley and Aler [86] enhanced the H S determination utilizing a varying P 0 .Based on the assumption that P 0 decreases as wave height increases, a linear relationship between P 0 and the mean image intensity in the up-wave direction was found.Then, a P 0 that varies with mean image intensity values was used in the model for the H S estimation.According to the test using datasets collected in the southern Labrador Sea from 30 January to 17 March 1997 [85,86] and a varying P 0 , the correlation coefficient and the mean difference between the radar-derived and wavemeter-measured H S were improved from 0.54 to 0.67, and from 28% to 20%, respectively, over that with a constant P 0 .In this method, illuminated and shadowed areas in the radar image are exploited to estimate H S .However, an external sensor is required for the calibration.Later, an algorithm also based on shadowing analysis incorporating sea surface slope and mean wave period was proposed, which does not require calibration [87].

Statistical Analysis
Gangeskar [88] estimated the wave height through analyzing the texture of radar images.The connection between H S and various parameters derived from the statistical analysis of images was investigated.Datasets used to test the algorithm were collected from the Gullfaks C (GFC) oil production platform in the North Sea, with H S ranging from 0.5 m to 9 m [88].17 statistical parameters were used to estimate H S .The average absolute correlation coefficient between the radar-derived and reference H S were 0.65 and 0.67, respectively, for the Cartesian start ranges (CSR) of 1500 m and 900 m.The best result was obtained using the parameter of the number of extrema per area in the smoothed image, with the standard deviations of 0.99 m for the CSR of 1500 m and 0.98 m for 900 m.
Later, Gangeskar [89] developed an adaptive method incorporating a scheme for detecting problematic situations such as low wave height and precipitation.In this adaptive method, different parameters were selected to estimate wave height according to the data quality identification results.Under good conditions, the parameter of extrema number is used, while under precipitation or low wave height, the parameter of gray level co-occurrence matrix (GLCM) correlation [89] is used.The result provided by the same datasets as in [88] showed that by using the adaptive method, the standard deviations were reduced to 0.68 m and 0.77 m for the CSR of 1500 m and 900 m, respectively.Although precipitation may affect radar images by inducing noise, the adaptive method still permits H S estimation.

Tilt-Based Algorithm
In most of the aforementioned cases, the algorithmatic outputs require calibration by additional reference sensors such as wave buoys.Dankert and Rosenthal [43,90] proposed an algorithm to determine ocean surface elevations based on tilt modulation without calibration.In their method, both hydrodynamic and shadowing modulations are considered to be insignificant to the radar imaging mechanism.
First, each image in a radar image sequence is corrected with a 2D antenna pattern, which is obtained from the temporal average of all the radar images.Next, a third-order polynomial f (r a , θ), where r a is the range, is least-squares fitted to the mean RCS for each azimuthal direction θ.With the given antenna height h ant , the depression angle as a function of range can be calculated as The local ocean surface tilt ϕ(r a , θ), which is assumed to be equal to the change of the local depression angle, can be obtained as where ψ(r a , θ) is the local depression angle, which can be determined as where σ r a is the local RCS at range r a , and f −1 is the inverse function of f .The local tilt angles are derived for all the images of an image sequence, resulting in a tilt image sequence.Finally, the ocean surface elevation is derived by convolving tilt images with a 2D integration filter kernel [90].Here, the convolution in the spatial domain is realized by the multiplication of the 3D DFT of the tilt image sequence and the integration transfer function in the Fourier domain.Then, the dispersion relationship with current effects and a 3D Gabor filter are used to extract the 3D wave spectrum.The ocean surface elevation is obtained by applying an inverse 3D DFT to the 3D wave spectrum.
The algorithm was tested using datasets collected in the central North Sea from February to September 2001 [43,90].During the test, H S measured by in-situ sensors varied from 0.5 m to 6 m.For the radar-derived and in-situ sensor-measured H S , a correlation coefficient of 0.94 and a standard deviation of 0.35 m were obtained.However, since shadowing modulation cannot be neglected for X-band radars operating at grazing incidence, the algorithm may not be suitable for the standard X-band marine radar mounted closer to sea level.

Shadowing-Based Algorithms
A. Original Algorithm Gangeskar [87] proposed an algorithm for the H S estimation based on shadowing analysis.This algorithm also does not require any calibrations using external sensors.Furthermore, it is suitable for the X-band marine radar operating at the grazing incidence.The basic steps of the shadowing-based algorithm are shown in Figure 11.
where i ranges from 1 to 8, representing 8 neighbouring pixel directions.Then, for each edge image I E i (m, n), a thresholding process is implemented so that the pixels with the upper B-percentile of pixel intensity levels are assigned the value of 1 and other pixels the value of 0, resulting in 8 thresholded edge images I T i (m, n).This is because I E i (m, n) represents the intensity level difference of a pixel with its ith-direction neighbouring pixel, and the difference is large for edge pixels that separate shadowed and unshadowed areas of an image.The thresholded edge images for 8 different directions are then summed as an overall edge image I T (m, n), which is given by If a pixel has edges in more than τ F directions, it is regarded as a single pixel noise and is removed from the overall edge image.Thus, the overall edge image I T (m, n) is filtered as Image pixels with the value of 1 in the filtered overall edge image I F (m, n) are identified as edge pixels.The mode of the raw image intensity levels of edge pixels is estimated as the shadow threshold τ S .Then, a shadow image is determined such that the image pixels with intensity levels lower than the shadow threshold τ S are shadowed while others are illuminated.Subsequently, the shadow image is divided into range-azimuth bins of 30 m in range and 10 • in azimuth.The illumination ratio for each bin is computed as the number of illuminated pixels divided by the total number of pixels, and the illumination ratio L(γ) as a function of grazing angle γ is obtained for each 10 • -azimuth sector.By curve fitting Smith's function [93] to the illumination ratio L(γ) using the Neilder-Mead simplex method [94], the RMS surface slope σ RMS may be derived for each sector.Then, the quadratic mean of the RMS surface slopes for all sectors of an image is computed as σ A RMS .Finally, H S is determined as where T m02 is the average zero-crossing wave period derived from the radar image sequence using existing algorithms [19,20].The algorithm was tested using datasets collected at Veslefrikk platform in the North Sea in January 2008 [87].The correlation coefficient between the radar-derived and reference H S was 0.87.The shadowing-based algorithm is totally based on theory and does not require calibration with additional sensors.However, in [87], the test was conducted using the average zero-crossing wave period T m02 measured by a reference sensor.Therefore, the algorithm needs to be further validated using the radar-derived average zero-crossing wave period T m02 in order to be fully independent of external sensors.In addition, the height of the radar used in the test [87] was 43.2 m, which was higher than the regular marine radar height of 20-30 m.Since the radar height affects the grazing angles and shadowing, the algorithm also needs to be validated using a regular marine radar.

B. Modified Algorithm
Later, Liu et al. [91,95] modified the original shadowing-based algorithm [87] by smoothing the edge pixel intensity histogram and selecting an appropriate subarea.Firstly, the edge pixel intensity histogram is smoothed using a spline function, and the shadow threshold is determined from the smoothed histogram.This results in a more robust shadow threshold estimation.Secondly, due to the dependence of radar backscatter on wind direction, shadow areas may be overestimated by using image portions far from the upwind direction, leading to overestimated wave height.Thus, after the shadow image is obtained, only the subarea over ±5 • around the upwind direction is extracted from the shadow image and used for one RMS surface slope calculation.The RMS surface slopes obtained from all images in one image sequence are averaged for one H S estimate.
The modified shadowing-based algorithm was tested using the same datasets described in Section 2.3.3-B, with H S varying from 1 m to 5 m [91,95].The standard SNR-based and original shadowing-based algorithms were also applied to the datasets.With respect to buoy-measured H S , the correlation coefficient of H S derived from the standard SNR-based, original shadowing-based, and modified shadowing-based algorithms were 0.76, 0.64, and 0.81, respectively.The RMS differences of H S derived from the three corresponding algorithms were 0.51 m, 2.17 m, and 0.47 m, respectively.Through the modifications, the H S estimation was significantly improved over that of the original shadowing-based algorithm.
C. Water-Depth-Incorporated Algorithm Since (35) is valid only under deep water conditions, Wei et al. [96] modified the shadowing-based algorithm in [87] by considering water depth for its application to a coastal area as where d w is the water depth, and λ the wavelength.x 1 = 0.098 and x 2 = 0.289 are determined by minimizing the RMS error between the function tanh(kd) and the piecewise linear function approximated to it.The water-depth-incorporated shadowing-based algorithm was tested using the datasets acquired at the coastal area of Haitan Island in Fujian Province from October to November 2010, and from November 2014 to January 2015 [96], with H S ranging from 0.5 m to 5 m.The correlation coefficient and the RMS difference between H S measured by the water-depth-incorporated algorithm and buoy are 0.68 and 0.51 m, respectively, while those for the original algorithm are 0.53 and 0.58 m, respectively.

Support Vector Regression Algorithm
Salcedo-Sanz et al. [97] estimated H S using a support vector regression (SVR) algorithm and shadowing information from simulated and field X-band marine radar images.In this algorithm, the SVR is trained from the simulated data and applied to field data.Calibrations using external sensors can thus be avoided.A classical model, ε-SVR [98], is considered in this algorithm.Shadow probability and the tangent of the local incident angle are included in the SVR training scheme using the simulated data.Shadowing probability is obtained as the ratio of the number of shadowed pixels to the number of total pixels.The tangent of the local incident angle is determined from the radar range, wave elevation, and antenna height above the sea surface.
After the training phase using the simulated data, the SVR is applied to field data collected at the German research platform of FINO 1 in the German basin of the North Sea in 2006 and 2013 [97].As compared to the ground truth, the mean absolute error (MAE) and mean square error (MSE) of H S obtained with the SVR algorithm were 0.66 m and 0.82 m, respectively, for data acquired in 2006 with H S varying from 2 m to 9 m.For data acquired in 2013 with H S varying from 2 m to 6 m, the MAE and MSE were 0.50 m and 0.38 m, respectively.

Empirical Orthogonal Function-Based Algorithms
A. Principal component Chen et al. [99] proposed an empirical orthogonal function (EOF)-based algorithm for wave measurements using X-band marine radar.Wave parameters, including H S , peak wave period, wavelength, and wave direction were determined from the principal components (PCs) of the radar image sequence derived from the EOF.The method is referred to as EOF-PC hereafter.
First, a 2D time-space image sequence can be represented by where C is the sample number in space, and D the sample number in time.After the application of EOF [100,101], I TS can be decomposed as The jth column of V, v j , is the jth eigenvector of the matrix I TS I T TS .The kth row of Z, z k , is the kth PC of the 2D time-space image sequence I TS , representing the variation of waves in the time domain.For H S estimation, a linear relationship between H S and the standard deviation of any one PC, std(z k ), is established as where p 0 and p 1 are parameters which can be determined from the least-squares fitting.With the obtained linear relationship, H S can be estimated from the standard deviation of one PC for any radar image sequence.Then, the first PC is used for the estimation of peak wave period.The maximum entropy method (MEM) [102] is used to derive the power spectral density (PSD) of the first PC.The peak wave frequency f p is determined as the frequency corresponding to the maximum PSD, and the peak wave period T p is determined as 1 f p .The process of estimating the peak wavelength is similar to that of peak wave period.The difference is that the 2D time-space image sequence to be decomposed by the EOF consists of C samples in time and D samples in space, as opposed to Equation (37).After the decomposition, PCs can be obtained in the space domain, representing the spatial variation of waves.Then, the peak wavelength can be obtained from the PSD of the first PC derived by the MEM [102].In the 2D space domain, two peak wavelengths along two orthogonal directions need to be derived, and the peak wavelength and wave direction can be obtained based on the geometric relationship.The ambiguity of 180 • of the wave direction is removed by applying the 3D DFT to the image sequence.The wave direction is selected as the one close to the peak of the image spectrum.The basic steps of the EOF-PC-based algorithm are shown in Figure 12.The EOF-PC-based algorithm was tested using the same datasets collected from two experiments as described in Section 3.1.1-D[99].During the experiments, buoy-measured H S and peak wave period ranged from 0.5 m to 3.5 m and from 5 s to 11 s, respectively.The comparison of the radar-derived and buoy-measured results showed a correlation coefficient of 0.69 and a RMS difference of 0.84 s for peak wave period and 0.93 and 0.21 m for H S .By applying the EOF, the algorithm enables wave parameter estimation from heterogeneous ocean wave field such as occur in near-shore regions.The linear wave dispersion relationship and MTF are also not required.However, a calibration process with an external sensor is still necessary for wave height estimation.

B. Rotated EOF
Since the first few PCs derived by the EOF may not represent the main modes of the wave field and wave parameters estimated from EOF analysis may be affected by the observed area and selected pixel positions, Yu et al. [103] and Zhang et al. [104] utilized a rotated EOF (REOF) to estimate H S .
After the EOF decomposition Equation ( 38), a diagonal matrix Λ is composed by the first P eigenvalues of I TS I T TS , with P being selected as the lowest value for which the ratio of the sum of the first P eigenvalues to the sum of all the eigenvalues is greater than 85%.An initial factor loading matrix is obtained as the multiplication of the first P columns of V and Λ, and normalized by dividing each element at the cth row by h 2 c , which is the sum of squares of all elements in the cth row of the initial factor loading matrix.Then, the factor loading matrix V is obtained by multiplying the normalized initial factor loading matrix with an orthogonal matrix [104] that maximizes where v cp is the element of V at the cth row and the pth column.Then, the common factor matrix Z is obtained from The qth row of Z , z q , is the qth rotated PC of the 2D time-space image sequence I TS .Finally, H S can be estimated from the standard deviation of any one rotated PC, std(z q ), by training a linear model as in Equation ( 39) or a polynomial model [104].
The REOF-based algorithm was validated using the datasets collected in the Bohai Bay of China from July 25 to 26, 2014, with H S ranging from 0.5 m to 1.8 m [103,104].The correlation coefficient and RMS difference between the buoy-measured and radar-derived H S from the REOF-based algorithm are 0.88 and 0.17 m, respectively, for a linear model, and 0.92 and 0.14 m, respectively, for a polynomial model.The EOF-PC-based algorithm was also applied to the same datasets.The corresponding correlation coefficient and RMS difference are 0.86 and 0.19 m, respectively.The results from the employed datasets indicated that the REOF-based algorithm improved the H S estimation over the EOF-PC-based algorithm.However, a calibration process using an external sensor is still required.
C. Joint PDF Later, by extending the EOF-PC-based algorithm, Chen et al. [105] developed an algorithm to retrieve wave heights from X-band marine radar image sequences based on the joint PDF of dimensionless wave periods and wave heights, without using external sensors for calibrations.
Based on Equations ( 37) and (38), a dominant wave field may be reconstructed as I 1 can be expressed in polar coordinates as I 1 (m, n).Since H S is proportional to the standard deviation of PCs, I 1 (m, n) is rescaled as • std(z 1 ).
Then, with the peak wave direction φ p determined by using the traditional 3D-DFT-based method, the radial profile I 1 (m, φ p ) is used to derive dimensionless wave periods and wave heights.Zero-crossing wavelength l 0 i is obtained from I 1 (m, φ p ) as the distance between two neighbouring crests or troughs.Similarly, zero-crossing wave heights can be obtained as where r a j and r a j+1 are ranges corresponding to two neighbouring crests and troughs, and a c is the coefficient to be determined in the following procedures.Zero-crossing wave periods t 0 i can be derived from l 0 i using the dispersion relationship.Subsequently, dimensionless wave periods are obtained as τ 0 i = t 0 i t 0 , where t 0 is the average zero-crossing wave period.Dimensionless wave heights are obtained as α 0 j = h 0 j h 0 , where h 0 is the average zero-crossing wave height.
The joint PDF of dimensionless wave periods and wave heights for random ocean waves, f (α h , τ t |υ 0 ), is defined in [106,107].The spectral width parameter υ 0 is determined by fitting the histogram of dimensionless wave period τ 0 i to the theoretical PDF of dimensionless wave periods which is derived by integrating f (α h , τ t |υ 0 ) over α h from 0 to 5 with a step of 0.001.With the derived υ 0 , the theoretical PDF of dimensionless wave heights can be obtained by integrating f (α h , τ t |υ 0 ) over τ t from 0 to 3 with a step of 0.001.The coefficient a c can be calculated as a c = p s t p s r , where p s r and p s t are the shaping parameters for the histograms of dimensionless wave height α 0 j and the theoretical PDF of dimensionless wave heights, respectively [105].Accordingly, H S can be determined as where h 0m is the wave height having the maximum probability.Overall, the basic steps of the joint-PDF-based algorithm are shown in Figure 13.The joint PDF-based algorithm was tested using datasets collected on Haitan Island of China from October, 2014 to January, 2015, with H S varying from 0.2 m to 3.5 m [105].The datasets contained 671 measurements for HH-polarized radar and 548 measurements for VV-polarized radar.The correlation coefficient and the RMS difference of H S derived from the HH-polarized radar and the buoy are 0.78 and 0.51 m, respectively, and those for the VV-polarized radar are 0.77 and 0.51 m, respectively.Wave height can be retrieved from the algorithm directly without using external sensors for calibration.However, compared to the SNR-and shadowing-based algorithms, the accuracy of the joint PDF-based algorithm needs to be improved.

EEMD-Based Algorithm
Most recently, Liu et al. [108] proposed an algorithm for the estimation of H S based on EEMD.In this algorithm, a sub-image located at range of 240 m to 1200 m and azimuth of ±30 • around the upwind direction is first extracted from each radar image for more distinct wave signatures.For a sub-image with M rows and N columns, the EEMD [56] is applied to each column of I sub .Then, the nth column can be represented as where D j (∼, n ) and R(∼, n ) are the jth IMF and the residual of I sub (∼, n ), respectively.Then, a normalization scheme [59] is applied to each column of D j to obtain where A j (∼, n ) and F j (∼, n ) are the AM and FM parts of D j (∼, n ), respectively.Subsequently, a parameter A H S is defined as Finally, a linear model that relates H S to A H S is sought through least-squares fitting.After the model is obtained, H S can be estimated from A H S for any radar image.
Only the Decca radar data from the datasets, with H S varying from 1 m to 5 m, described in Section 2.3.3-B were used to test the EEMD-based algorithm [108].With respect to buoy reference, the correlation coefficient of H S derived from the standard SNR-based, modified shadowing-based, and EEMD-based algorithms were 0.61, 0.84, and 0.91, respectively.The RMS differences for the three algorithms were 0.78 m, 0.48 m, and 0.36 m, respectively.The algorithm enables H S estimation from a ship-borne radar.However, calibration with an external sensor is still necessary.

Coherent Radar
All of the work discussed in the previous sections is based on the analysis of non-coherent radar data.Presently, X-band coherent marine radar is undergoing rapid development for wave field and sea surface current measurements [109][110][111][112][113]. The advantage of coherent radar is that it provides not only the image intensity but also the radial velocity of the sea surface scatterers.Based on this feature, algorithms have been developed to estimate wave height directly from the radar data without calibration.
Hwang et al. [114] proposed an empirical method to estimate H S from the RMS value of Doppler velocities of the scattering element and peak wave frequency.The Doppler velocity can be derived by a pulse pair approach or through the DFT.In the pulse pair approach, the Hilbert transform is applied to each radar pulse to obtain the complex radar signal.Then, the Doppler frequency, ω D , is obtained as the temporal derivative of the phase, φ D , of the complex signal.The Doppler velocity of the scattering element is calculated as where λ r is the wavelength of the transmitted radar signal and γ is the grazing angle.Alternatively, the Doppler frequency can be derived from the first moment of the Doppler frequency spectrum by using the DFT approach [115].The comparison between the pulse pair and DFT approaches in [114] shows that either works well for the determination of the Doppler velocity.After u D is obtained, the peak wave frequency ω P is derived from the surface displacement spectrum, which is generated from the Doppler velocity wavenumber spectrum [114].Finally, H S is determined as where u D RMS is the RMS value of the Doppler velocities and X is an empirical coefficient.The technique was tested using dual-polarized X-band coherent radar data collected from a tower about 60 km off the Georgia coast from 9 to 11 April 2006 [114].It was found that the empirical coefficient X was 1 for the HH polarization and 1.3 for the VV polarization.However, with the Doppler velocity information obtained from the coherent radar, an empirical relationship is still required for H S estimation.Recently, Carrasco et al. [116] proposed a semi-empirical method for retrieving H S .In that method, the Doppler frequency is also derived from the temporal derivative of the phase of the complex radar signal, and the Doppler velocity is also derived from (49).Based on the analysis of 15-min time series of the Doppler velocities obtained at 4 different ranges, it was found in [116] that the standard deviation of the Doppler velocities is almost the same as that of the sea surface heave recorded by buoy.Since H S is calculated as 4 times the standard deviation of the sea surface heave, it can be estimated from the Doppler velocity as where τ D is the standard deviation of the time series of the Doppler velocities in one range.The method was tested using the data collected from the German research platform Fino-3 in the German Bight in the southern North Sea from 6 March to 14 July 2015.During each data collection period, the radar first operated in the rotational mode to derive the peak wave direction, and then operated in the static mode (i.e., pointing in the peak wave direction) to derive the Doppler velocity and to estimate H S .The correlation coefficient and RMS difference of radar-derived and buoy-measured H S were 0.96 and 0.23 m, respectively, with H S up to 5 m.The result was shown to be comparable to the best result obtained from a non-coherent radar but does not require external calibration.It was found that H S tended to be overestimated for wind seas and underestimated for swell-dominated seas.Furthermore, the projection effects may lead to an underestimated H S , especially in multi-mode seas.Later, Carrasco et al. [117] developed a non-empirical, calibration-free algorithm for H S estimation by considering the energy loss due to the projection effects.The same pulse pair approach as in [116] is first used to achieve the Doppler velocity.Then, shadowed areas are identified and removed from the Doppler velocity time-range series by using a confidence measure.After that, a 2D DFT is applied to obtain the Doppler velocity wavenumber-frequency spectrum, and this is followed by a filtering with the linear wave dispersion relation to remove non-periodic contributions to the Doppler velocity.Then, the resulting 2D Doppler velocity spectrum is transformed to the 2D surface elevation spectrum based on linear wave theory [116].The integration of the 2D surface elevation spectrum over wavenumber and frequency domain results in m 0,P , which is the total variance of the surface elevation spectrum projected in the radar look direction.As in [116], a radar operating in both rotational and static modes is employed.The 2D wavenumber spectrum derived from the radar image sequence acquired in the rotational mode is used to obtain the projection loss ratio, r P , which is calculated as the ratio of the spectral energy projected in the radar look direction to the total spectral energy.Finally, H S is determined as For wind measurements, wind direction is generally determined from the upwind peak in images collected from HH-polarized X-band radars operating at grazing incidence; wind speed is generally determined from parameters related to the NRCS derived from the radar image.Recently, much effort has been focused on techniques for wind parameter extraction from rain-contaminated radar data, which may still be one of the future trends for wind measurements.Another remaining challenge for wind measurements is that a calibration phase using external sensors is generally required for wind speed estimation.In the future, calibration-free algorithms for wind speed estimation need to be developed.This may be realized by exploiting both the spatial and temporal variations of image patterns.Techniques such as the optical flow-based motion estimation in [45] may be further developed to measure the movements of wind gusts visible in images.
For wave measurements, various algorithms based on the spectral or texture analysis of radar images have been used.Traditionally, wave height estimation requires calibrations using external sensors.Recently, techniques have been rapidly evolving for calibration-free wave height estimation.These include tilt- [90] and shadowing-based algorithms [87], and an algorithm exploiting the joint PDF of dimensionless wave periods and wave heights [105].Calibration-free wave height estimation may still be one of the future trends for wave measurements.Coherent X-band radar has also been proved to be an effective solution for such purposes since it provides the radial velocity of ocean surface waves that can be used to derive H S [116,117].In addition, a simple model for marine radar backscatter signals from the ocean surface has been developed [118], which may be further investigated for calibration-free wave height estimation.One of the remaining issues of wave measurements is that no effective algorithms have been developed for mitigating rain effects.With regard to this, S-band radar is a promising tool since, due to its longer wavelength, it can yield wave measurements under rain conditions [119].Recently, models for SAR backscatter from the ocean surface under rainfall have been developed [120,121].X-band marine radar sea surface imaging mechanisms in rain events are also worth exploring.

Figure 2 .
Figure 2. Scheme of the local gradients method (LGM) for wind direction estimation.

Figure 3 .
Figure 3. Flowchart of the optical flow motion estimation-based technique (OFM).

Figure 4 .
Figure 4. Basic steps of the single-Curve Fitting (CF) method for wind direction and speed retrieval.

Figure 5 .
Figure 5. Basic steps of the original Intensity Level Selection (ILS)-based method for wind parameter estimation.

F
( k, ω) from the high-pass filtered image spectrum F (3) I ( k, ω).Because of the non-linear radar imaging process, a modulation transfer function (MTF) T M ( k) is developed to estimate the wave spectrum F

Figure 10 .
Figure 10.Basic steps of the array beamforming algorithm.

Figure 12 .
Figure 12.Basic steps of the empirical orthogonal function (EOF)-principal component (PC)-based algorithm for: (a) H S and peak wave period; (b) peak wavelength and wave direction.

Figure 13 .
Figure 13.Basic steps of the joint-probability distribution function (PDF)-based algorithm.

Table 1 .
Inter-Comparison of the Reviewed Wind Measurement Techniques.