Wave Height Estimation from Shipborne X-Band Nautical Radar Images

A shadowing-analysis-based algorithm is modified to estimate significant wave height from shipborne X-band nautical radar images. Shadowed areas are first extracted from the image through edge detection. Smith’s function fit is then applied to illumination ratios to derive the root mean square (RMS) surface slope. From the RMS surface slope and the mean wave period, the significant wave height is estimated. A data quality control process is implemented to exclude rain-contaminated and low-backscatter images. A smoothing scheme is applied to the gray scale intensity histogram of edge pixels to improve the accuracy of the shadow threshold determination. Rather than a single full shadow image, a time sequence of shadow image subareas surrounding the upwinddirection is used to calculate the average RMS surface slope. It has been found that the wave height retrieved from the modified algorithm is underestimated under rain and storm conditions and overestimated for cases with low wind speed.Themodifiedmethod produces promising results by comparing radar-derived wave heights with buoy data, and the RMS difference is found be 0.59m.


Introduction
Marine radars can image both temporal and spatial variations of the sea surface while buoys provide only temporal point measurements.The radar signature of the sea surface, also known as sea clutter, is undesirable and generally suppressed in navigation purposes, but it is useful in monitoring the sea state [1].The clutter is, in general, generated by the Bragg scattering of the near-grazing incidence radar signal with short wind-induced sea surface ripples.Longer waves become visible in radar images due to their modulations of the short ripples, mainly via hydrodynamic modulation, tilt modulation, and shadowing [2].Thus, analysis of time series of X-band nautical radar sea surface images allows the estimation of directional wave spectra and integrated sea state parameters [3][4][5].Algorithms for such purposes have been largely developed during the last several decades.
A widely accepted method of wave height estimation for X-band radar is based on the signal-to-noise ratio (SNR) derived from the image spectrum [5][6][7].Another class of algorithms is based on the statistics of radar sea surface images.Through a constant threshold probability of illumination  0 based on the theory of geometric optics [8], a model relating the significant wave height to the island-totrough ratio extracted from the image was established in [9].In [10], an algorithm based on [9], but with a varying  0 , was proposed to enhance the wave height determination.In [11,12], wave height was derived by analyzing the texture of Xband radar sea surface images.Other techniques, including an iterative least square approach [13] and a wavelet-based algorithm [14], have been similarly developed.In all cases, the algorithmic outputs require calibration by additional reference sensors such as wave buoys.
Recently, a shadowing-analysis-based wave height algorithm has been proposed [15].Assuming a geometric shadowing condition, shadowed areas are first extracted from the image by edge detection.Then, using the calculated illumination ratios in local areas, the RMS surface slope is derived by curve-fitting Smith's function [16].Finally, the significant wave height is estimated from the RMS surface slope and the average zero-crossing wave period.Unlike earlier approaches, this algorithm does not require calibrations using additional reference sensors.Therefore, it shows promise for improving the ease of implementation and reducing operational cost.Here, such a method is modified and applied to data acquired from a radar on a moving ship.

Method
2.1.The Shadowing-Analysis-Based Wave Algorithm.The algorithm was introduced in detail in [15].It contains the following major steps.

Estimating Shadow Threshold.
The edges in the image that separate shadowed areas from illuminated areas are identified using an edge detection technique.Here, this involves the convolution of a raw radar image (, ) with a simple pixel difference operator   (, ) for each of  = 1, 2, . . ., 8 directions, and  and  are range and azimuth, respectively.In [17], this convolution results in  edge-detected images (  (, )) given by   (, ) =  (, ) ⊗   (, ) . ( A thresholding process is then applied to the eight edge images using a threshold value equal to the highest percentile of the pixels.Edge image pixels whose intensity levels are higher than the threshold are assigned the value of 1, and the remaining pixels take the value of 0. The process results in eight thresholded edge images   (, ).
Subsequently, an overall edge image   (, ) is obtained by combining the eight thresholded edge images   (, ) and a filtering process given by [15]   (, ) = The filtering is implemented to remove the single pixel noise that has edges in more than   directions.
Using the raw radar image pixels corresponding to the pixels of intensity value of 1 in   (, ), a statistical distribution   () of gray level values  is created.From the distribution, the gray level threshold   used to identify shadow can be determined as [15]   = mode (  ()) . (4)

Calculating Illumination Ratio.
Using the shadow threshold   determined in Section 2.1.1,the shadow image can be derived.Pixels with the values less than   are regarded as shadowed, and the remaining pixels are understood to be illuminated.Next, the shadow image is divided into segments along the range  and the azimuth , and the illumination ratio () as a function of grazing angle  for each segment is calculated [15].

Estimating RMS Surface Slope.
Having obtained the illumination ratios for each azimuth direction, the RMS surface slope  RMS of a random rough surface described by a one-dimensional Gaussian surface height probability density function (PDF) can be derived by curve-fitting Smith's function [16] for that direction.The curve fitting is implemented by the Neilder-Mead simplex method in one dimension [18].After deriving  RMS for all azimuth angles, an average RMS surface slope   RMS can be calculated.
2.1.4.Estimating Significant Wave Height.Finally, from the average RMS surface slope   RMS and the average zerocrossing wave period  02 , the significant wave height  0 can be determined as in [15] by where  is the gravitational acceleration. 02 can be derived from the radar images themselves using existing wave algorithms [1][2][3][4][5].However, the average zero-crossing wave period  02 measured by well-accepted buoy data instead of radar was used here since the radar employed in this study did not produce a good value for  02 .

Data Quality Control.
Before processing the raw radar data, low quality images, such as rain-contaminated or lowbackscatter cases, should be discarded.The data quality control procedure described in [19] is employed here.Rain leads to the change of the normalized radar cross section (NRCS) and significantly affects wave height retrieval.Due to the strong impact of rain on the number of zerointensity pixels in X-band nautical radar images, the zeropixel percentage (ZPP), which is defined as the ratio of the number of zero-intensity pixels to the total number of pixels in an image, is identified as a parameter for rain recognition [20].For the data used here, pixels with gray scale intensity lower than 5 are regarded as zero-intensity pixels [19].Images with ZPP less than 10% are considered as rain-contaminated and are not used for wave height retrieval.
Low-backscatter images that appear almost completely black due to low wind speed or unknown system errors contain little or no wave information.A parameter called low-clutter direction percentage (LCDP), which is defined as the ratio of the number of low-clutter directions to total number of directions in an image, is used for identifying low-backscatter images [19].If the ZPP of a single azimuth direction is higher than 40% (empirical and also varies with systems), the direction is regarded as a low-clutter direction.Then, the images with LCDP higher than 90% are excluded from subsequent processing.

Edge Pixel Intensity Histogram Smoothing.
In [15], the shadow threshold is directly determined as the intensity value corresponding to the highest occurrence of the intensity histogram of edge pixels.This is viable if the distribution is smooth.However, the data used in this work has a small gray scale depth (8-bit, i.e., 0-255) and a relatively small number of pixels.The shadow threshold may not be correctly determined by seeking the highest occurrence of the histogram.In order to improve the accuracy of shadow threshold determination, a smoothing process using a spline function is applied to the edge pixel intensity histogram.The smoothing spline function [21], , minimizes where  is the smoothing parameter determining the tradeoff between fidelity to the data and smoothness of the function.
When  approaches 0, the function converges to a simple linear least squares regression.When  approaches infinity, the function converges to the interpolating spline.Here,  is selected as 9 since the tendency can be maintained and the outlier of the distribution can be removed with this choice.This smoothing process is performed using the MATLAB built-in curve-fitting function "fit".It should be noted that other smoothing methods (e.g., moving average, median filtering) can remove outliers but may change the details of the tendency of a distribution.Thus, an inaccurate mode of the edge pixel intensity histogram may be obtained, leading to an inaccurate shadow threshold.However, the spline function does not cause such a problem.The shadow threshold is estimated from the smoothed intensity histogram.

Subareas Selection.
In [15], all the azimuth directions in the shadow image were used for the derivation of the average RMS surface slope.However, the image portion in the azimuths far from the upwind direction usually has low sea clutter intensities and overestimated shadowed areas.Including such a portion may result in wave height overestimation.
Here, for each shadow image, only a subarea selected from the portion ±5 ∘ around the upwind direction is used for RMS surface slope estimations since the clutter signal is stronger in those directions and more robust results may be obtained.The technique for determining the upwind direction is based on a dual-curve-fitting algorithm found in [19].Then, a RMS surface slope is derived from the subarea for each image.Moreover, an average RMS surface slope is calculated using a time sequence of images instead of all the subareas in a single image.Since the RMS surface slope is calculated from the portion around the upwind direction in each single image, the variation of RMS surface slope obtained from the upwind direction is negligible between consecutive images.Thus, it does not require two consecutive radar images to be perfectly overlapped, and the ship motion will not affect the result significantly.As with most existing works on the topic, a typical value of 32 images is used here to obtain an average RMS surface slope for wave height estimation.November 27, 0:00 November 28, 0:00 November 29, 0:00 0 Canada (DRDC) is used.The data was collected from 23:43 November 26 to 12:06 November 29 (2008), in a sea trial approximately 300 km south-southeast of Halifax, Nova Scotia, Canada.Three free-floating Triaxys directional wave buoys were deployed about 4 to 15 km apart to measure the wave field.The distances between the ship and the buoys were generally less than 10 km, but occasionally up to 15 km throughout the trial.The water depth around the ship and the buoys is about 200 m. Figure 1 depicts the ship's course and the track of the three wave buoys.Figure 2 displays the time sequence of ship speed throughout the sea trial.

Data
The radar utilized in the sea trial was a standard HHpolarized Decca marine radar which operated at 9.41 GHz with a sampling frequency of 20 MHz.The radar covered 360 ∘ in azimuth with a beam width of 2 ∘ and an antenna rotation speed of 28 rpm.The antenna was installed at a height of 21.9 m above the sea level, covering a range from 240 m to 2160 m with a range resolution of 7.5 m.The radar was connected to a Wave Monitoring System II (WaMoS II) [22].As indicated in Section 2.2.2, the system scaled and stored the radar backscatter power in gray levels from 0 to 255 (8-bit unsigned integers), with 0 corresponding to lowest radar return (black colour in the radar image) and 255 to the maximum radar return (white colour in the radar image).

Experimental Results.
One example of the B-scan (i.e., polar coordinate) raw radar image is shown in Figure 3.
The corresponding edge image obtained by edge detection and filtering is depicted in Figure 4, in which edges are shown in black.The threshold values used for the edge detection and filtering in this study are 20% and 5 (  in (3)), respectively.The requirement of a higher threshold for the edge detection than that (10%) in [15] in order to produce robust results was likely due to the different operational parameters including lower antenna height and lower range and azimuth resolutions used here.In [23], an additional filtering process with a constant threshold was introduced to remove those edges located too far away from shadow in order to obtain the reasonable shadow threshold from the edge pixel intensity histogram.However, from the result in [23], it can be seen that filtering with a constant threshold was not robust to variation in sea state.Here, the smoothing process is used instead.The intensity distributions of the edge pixels (red circles), along with the spline-fitted curve (red line), and the entire set of image pixels (blue points) for the image in Figure 3 are shown in Figure 5, in which pixels with zero intensity level are excluded, and the gray level shadow threshold is illustrated by the dash-dot.Note that if no smoothing is used, the shadow threshold will be determined by the outlier for the histogram curve, corresponding to the level intensity of 63 rather than a correct value (40 in Figure 5).It should also be noted that the shadow threshold is estimated for each image.Thus, the ship motion does not have a significant effect on the shadow threshold estimation.After thresholding the raw image in Figure 3, the corresponding shadow image is obtained and shown in Figure 6.In that figure, shadowed areas are shown as black and the subarea used for RMS surface slope calculation is the portion between the dashes-dots.Figure 7 depicts the illumination ratio as a function of grazing angle and the corresponding Smith's function fit for one single subarea (the portion between the dashes-dots in Figure 6).The RMS surface slope  RMS is estimated by the curve fitting.The threshold that distinguishes the usable data and eliminated data in Figure 7 is sought by gradually eliminating the illumination ratio data used for the  RMS calculation from the longest range towards the direction of decreasing ranges (i.e., increasing grazing angles).This threshold is determined as the range beyond which the corresponding illumination ratio data is excluded from calculating  RMS and the  RMS obtained from the remaining data is the smallest [15].This is done because the radar backscatter from long ranges may be weak simply due to the decay law for the electromagnetic energy.This causes an overestimation of shadowed areas at long ranges, leading to a corresponding overestimated  RMS and wave height.
The original and modified shadowing-based algorithms described above are both applied to the quality-controlled Decca radar data, and the results are compared with the reference data measured by the buoys.The comparison of the time sequences of significant wave heights is displayed in Figure 8.It should be noted that a storm appeared between 2:30 and 12:00 on November 28, and radar data was not recorded for most of this period.Moreover, low quality images were discarded by the data quality control process.It can be observed that the wave heights obtained using the original algorithm are consistently overestimated.The overestimation is mainly due to the overestimated shadow threshold caused by the outliers in the intensity distribution of the edge pixels (see Figure 5).However, the radar results derived from the modified algorithm agree well with the buoy data for most of the period.Differences are observed over some periods.The wave heights were underestimated from 4:20 to 6:30 on November 27 and overestimated from 2:00 to 2:30 on November 29.During these periods, light rain and relatively low wind speed occurred.These conditions were not detected by the data quality control process.Since rain enhanced the image intensity, shadowed areas were reduced.Thus, wave heights were underestimated.However, low wind speed resulted in overestimation of shadowed areas and wave heights.From 2:50 to 3:50 and from 11:20 to 11:50 on November 28, wave heights were high but also underestimated.During this period, a strong swell signature was observed in addition to the wind wave component.An example of the wave frequency spectrum derived from buoy data during this period is given in Figure 9, in which the dual-mode wave field is clearly seen.The peak frequency in Figure 9 is 0.03 Hz which corresponds to a swell wave period of 33.3 s.Therefore, wave height estimation for such a complex sea state needs to be further analyzed.The corresponding scatter plots of the retrieved significant wave heights using the original and modified algorithms with the reference data are shown in Figures 10(a

Conclusions and Outlook
In this paper, a modified shadowing-analysis-based wave height estimation method has been applied to X-band nautical radar data.The modifications include (1) a data quality control process to exclude rain cases and low-backscatter images; (2) a scheme for smoothing the edge pixel intensity histogram to determine shadow threshold; and (3) employing of a time sequence of subareas around the upwind direction to calculate the average RMS surface slope.By comparing the radar-derived results and the buoy-measured data, it has been found that the wave height retrieved from the algorithm is underestimated under rain and storm conditions and overestimated under low wind speed.Still, the proposed method produces promising results, with a RMS difference of 0.59 m and a correlation coefficient of 0.91.However, in order to improve the robustness of this wave height algorithm, the effects of rain, low wind speed, and storm conditions with dual-mode wave fields need to be further analyzed.This will be more intensively investigated in the next phase of the work.Moreover, the algorithm needs to be further validated using radar data that can produce good estimation of average zerocrossing wave periods to make it completely independent of external sensors.

Figure 2 :
Figure 2: Time sequence of ship speed throughout the sea trial.

Figure 5 :
Figure 5: distribution of edge pixels and image pixels.Smoothed histogram is indicated by solid line and shadow threshold is illustrated by dash-dots.

Figure 6 :Figure 7 :
Figure 6: Shadow image.Shadowed areas are black.Subarea is the portion between the dash-dots.

Figure 8 :
Figure 8: Comparison of the time sequences of significant wave heights derived by the original and the modified shadowing-based algorithm and buoys.

Figure 9 :Figure 10 :
Figure 9: Wave frequency spectrum containing wind wave and strong swell components.