Abstract
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 upwind direction 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. The modified method produces promising results by comparing radar-derived wave heights with buoy data, and the RMS difference is found be 0.59 m.
1. 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–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–7]. Another class of algorithms is based on the statistics of radar sea surface images. Through a constant threshold probability of illumination based on the theory of geometric optics [8], a model relating the significant wave height to the island-to-trough ratio extracted from the image was established in [9]. In [10], an algorithm based on [9], but with a varying , was proposed to enhance the wave height determination. In [11, 12], wave height was derived by analyzing the texture of X-band 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.
This paper proceeds as follows. In Section 2, the shadowing-analysis-based wave height algorithm is briefly reviewed and the proposed modifications are described. Section 3 contains experimental results obtained from shipborne radar data using both the original and modified algorithms, as well as the comparison with buoy data. Finally, a conclusion and future directions for this work appear in Section 4.
2. Method
2.1. The Shadowing-Analysis-Based Wave Algorithm
The algorithm was introduced in detail in [15]. It contains the following major steps.
2.1.1. 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 directions, and and are range and azimuth, respectively. In [17], this convolution results in edge-detected images given byA 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]
2.1.2. 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].
2.1.3. Estimating RMS Surface Slope
Having obtained the illumination ratios for each azimuth direction, the RMS surface slope 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 for all azimuth angles, an average RMS surface slope can be calculated.
2.1.4. Estimating Significant Wave Height
Finally, from the average RMS surface slope and the average zero-crossing wave period , the significant wave height can be determined as in [15] bywhere is the gravitational acceleration. can be derived from the radar images themselves using existing wave algorithms [1–5]. However, the average zero-crossing wave period 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 .
2.2. Modification
2.2.1. Data Quality Control
Before processing the raw radar data, low quality images, such as rain-contaminated or low-backscatter 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 zero-intensity pixels in X-band nautical radar images, the zero-pixel 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.
2.2.2. 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], , minimizeswhere 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.
2.2.3. 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 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.
3. Results
3.1. Data Overview
In order to test the modified method, data provided by Defence Research and Development 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.
The radar utilized in the sea trial was a standard HH-polarized 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).
3.2. 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 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 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 and the 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 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) and 10(b), respectively. With the modification, the correlation coefficient is increased from 0.81 to 0.91, and the RMS difference is reduced significantly from 1.82 m to 0.59 m. By excluding the data with buoy-recorded wave heights larger than 8 m, the correlation coefficients between the buoy-recorded and radar-derived wave heights are 0.47 and 0.68, respectively, for the original and modified algorithms. The corresponding RMS differences are 1.79 m and 0.50 m, respectively.
(a)
(b)
4. 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 a data quality control process to exclude rain cases and low-backscatter images; a scheme for smoothing the edge pixel intensity histogram to determine shadow threshold; and 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 zero-crossing wave periods to make it completely independent of external sensors.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The work was supported in part by Natural Sciences and Engineering Research Council of Canada (NSERC) under grants to Weimin Huang (NSERC 402313-2012) and Eric W. Gill (NSERC RGPIN-2015-05289) and by an Atlantic Innovation Fund Award (Eric W. Gill, Principal Investigator). The work was also supported by the Department of Innovation, Business and Rural Development of Newfoundland and Labrador under Grant 30-10921-008. In addition, The authors would like to thank Dr. E. Thornhill at Defence Research and Development Canada (DRDC) for the provision of the radar and buoy data.