Roughness Change Analysis of Sea Surface From Visible Images by Fractals

This paper presents an analysis of sea surface roughness from visible images based on fractals features for the first time. We measure fractal dimensions of sea surface images by the methods of box counting, fractional Brownian motion, and area measurement. The empirical values of sea surface roughness are found from measured wind speed and an empirical relation. The correlations between wind speeds (or sea surface roughness) and fractal dimensions are evaluated based on the data from a field experiment. Further, four types of noises with different parameters are introduced to the sea surface images. Then noise suppression performances of six methods are evaluated. Our experiments have demonstrated that our fractal methods and empirical relation between wind speeds and fractal dimensions are effective for analyzing roughness changes of sea surface from visible images.


I. INTRODUCTION
Ocean surface waves significantly enhance the heat and gas exchanges between ocean and atmosphere. Sea surface roughness, a parameter of wavy surface in surface flux formulations, has a strong dependency on wind speed, and can also be modulated by large scale flow motions, such as, internal waves, tidal flows, surface currents, and rains [1]. That rain reduces the roughness of the sea surface is a belief commonly held among seafarers. Rain alters the surface roughness through the production of wavelets by raindrops, Rayleigh jets, and cavities, as well as by damping of highfrequency waves [2]. Coupled atmosphere-ocean models will potentially increase skills of the forecasts, and therefore, may simulate the large circulations more faithfully depending on the adequate surface parameterization schemes [3].
To represent the feedbacks to atmosphere which alters the coupling between atmosphere and ocean, wave dependency is introduced to the sea surface roughness parameter in The associate editor coordinating the review of this manuscript and approving it for publication was Guitao Cao . surface flux parameterization [4]. Zhang et al. [5] have evaluated different surface roughness parameterization schemes on numerical simulation of typhoon Vongfong. It is clearly shown that sea surface roughness can change the distribution of significant wave height. To direct measure sea surface roughness at sea is challenging under harsh environment conditions, the roughness is commonly computed through other wind and sea variables based on various formulas that have been proposed since Charnock (1955). Charnock's formula is based on a dimensional argument. No firm consensus has yet been reached on which of the expression is superior to others in general [6].
Over more than a half century, a number of methods for analysis and extract sea surface roughness from observations have been proposed. During the Southern Ocean Waves Experiment (SOWEX), images of sea surface were recorded from a low flying aircraft with a downwardlooking video camera. Walsh et al. [7] developed a numerical simulation to relate contrast of sea surface image to the portion of sea-surface mean square slope (MSS) due to gravity-capillary waves under light wind conditions. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Crossingham et al. [8] described a laser slope gauge that measures the deflection of a laser beam that is refracted by wavy air-sea interface. Microwave radiometer [9]- [20], scatterometer [12]- [14], [22], synthetic aperture radar [15]- [17], sonar [18] and satellite altimeter [19], [21] have also been used to measure and analyze sea surface roughness.
Here we first introduce a novel concept of using fractal dimension to quantifying roughness of wavy ocean surface.
The assumption for fractal modeling of sea surface roughness is that the fractality of sea surface can be directly expressed in the fractal spectrum, and vice versa, the fractal spectrum, in turn, relies on the same parameters that used to characterize the sea surface. Once these parameters are selected, representing measurements of the two-dimensional stochastic process, i.e. sea surface, can be constructed. Therefore, characteristics of sea surface, such as, MSS, significant sea height, roughness, may be mapped.
Through examining the properties of the sea surface fractal models at a great extend, as demonstrated from this study, it is shown that a suitable estimator can be found to extract the sea surface roughness from sea surface image. Our method is quite intuitive, and it is demonstrated empirically that the visible image of sea surface changes its fractal characteristics as sea state changes. This has also been highlighted from early studies of estimating the statistics of the sea clutter under various sea states [23].
Our study focuses on three aspects: fractal dimension measurements from visible images, estimations of wind speed (or sea surface roughness) functional depedence on fractal dimensions based on the extracted data, and effects of four types of noises and four types of image differentiates. The proposed box counting, fractional Brownian motion and area measurement approaches can accurately calculate fractal dimensions of sea surfaces from visible images, and the proposed methods can correctly estimate roughness change of sea surface from visible images by using six fractal methods under light wind conditions. One motivation of our research is ocean wave observations, measurements, 3D reconstructions and property estimations based on wave features extracted from 2D visible images. These methodologies developed can also be applied to analysis of other transparent and opaque object surfaces as well.
This paper is organized as follows: Section 2 proposes and constructs fractal dimension measurements of sea surface from visible images by six fractal descriptors. Section 3 outlines relations of sea surface roughness and wind speed. Section 4 presents the experimental results of fractal dimension measurements, empirical fittings, error analysis, effects of noises, transformation analysis, oil spill detection and so on.

II. FRACTAL DIMENSION MEASUREMENTS
In analysis of surface wave dynamics, it is sometimes convenient to describe the shape of sea surface in Fourier linear expansion. However, an ocean surface wave field is a nonlinear and random process. For the purpose of extracting wave field parameters, such as surface roughness, it can be more efficient, we believe, to represent the surface image directly by nonlinear descriptors, such as, to describe the image complicity at different scales by fractal geometry, attempted in this paper. The fractal dimension (FD) of a bounded self-similarity set A in Euclidean n-dimensional space, by definition, can be expressed as [24]: where N r is the least number of distinct possible patterns in A at the scale r. The union of N r distinct patterns should be able to cover the corresponding self-similarity set completely. This section describes how to systematically construct fractal dimensions for measuring sea surface roughness from nearshore still images. Six fractal methods constructed are respectively box counting (BC), improved box counting (IBC), gray value statistic (GVS), histogram mean (HM), power spectral density (PSD) and area measurement (AM).

A. BOX COUNTING APPROACHES 1) BOX COUNTING METHOD
Our calculation of box counting dimensions is extended from [25] and described as follows. A grid size M × M graylevel image f (x, y) is partitioned into subimages of r × r size (1 < r ≤ M /2, for integer r). r is associated with the size of corresponding cubic box. Let the gray values of the four corners of the subimage (i, j) be I k (k = 1, 2, 3, 4), then we have where n r (i, j) is the number of size r boxes that are needed to fill up the covered volume of the subimage (i, j). N r is the total number of the boxes covering the whole image, If we rewrite (1) for D in an equivalent form of N r · r D = Const, the box-counting dimension D can be found by the tangent of the function curve, (log(N r ), −log(r)). It can be estimated by the least square linear fitting of log(N r ) against −log(r). There definitely exists a fractal scale in the space scale range of [r 1 , r 2 ] in which the FD estimate of the ideal fractal image can be computed through where N r 1 and N r 2 are the numbers of boxes needed to cover the ideal fractal surface at the box scales r 1 and r 2 respectively [25]. For example, a flat surface, or, a constant gradient image, D should ideally be equal to 1 from BC. D is independent on orientation of images by its definition. FD estimate error by the above scheme will increase with the increase of r. A parameter γ is introduced here to represent scale of actual volume and maximum volume. γ is defined as [24]: where I i (i = 1, 2, · · · , r × r) is the gray value of i-th pixel in the image window A 1 of size r × r (See Fig. 1), I max and I min are respectively the highest and lowest gray values in A 1 . The larger the value of γ is, the closer the gray value of each pixel is to I max . This can cause that the boxes cannot be filled into A 1 and the value of n r (i, j) calculated from (2) is too large. So with the increase of γ , the number of filled boxes should be reduced. Through the experimental test and analysis comparison, we conclude that γ = 0.8 is the critical point. The computing procedure of the improved method can be described as follows. First, I max , I min , γ of A 1 are calculated. If γ > 0.8, n r (i, j) is set to 1, if γ ≤ 0.8, n r (i, j) is calculated using (2). This method overcomes the shortcoming of the above method and improves the accuracy.

B. FRACTIONAL BROWNIAN MOTION APPROACHES 1) GRAY VALUE STATISTIC METHOD
FBM is constructed by the image gray scale increments in x and y directions [26]. For our application, the intensity variations of an image, i.e., surface slope variations of random wavy water surface, in space is measured. This method is simple and fast in computations. The algorithmic outline for FD on FBM is given below: 1. The means of increments in the x, y directions are calculated in the subwindow scale r= 1, 2, 3, 5 (r≥1) by the following equations: where I (x, y) is the pixel gray value at (x, y) in the image window of N ×N size. We randomly set N = 140.

Hurst exponent H satisfies the equation
where C is an arbitrary constant. H can be calculated by regression fitting the following system of equations [26]: Substituting three sets of parameters r 1 , r 2 , E r 1 , E r 2 into (10), we obtain three values of H and calculate average value of these three values.
4. The fractal dimension, GVS, of FBM D can be related to the Hurst exponent H by where d denotes the space dimension. We set d= 2.

2) HISTOGRAM MEAN METHOD
A random function I (x) is said to be a fractional Brownian function if for all x and x [27]: where F(y) is a cumulative distribution function and Pr is the probability. Note that x can be interpreted as vector quantities, thus can be extended to two or more topological dimensions as we will do late in the section.
If I (x), in our case the image intensity, is a fractal Brownian function, we can rewrite (12) by taking ensemble average and found that second-order statistics of the image change with scale r = x [27]: where Histogram mean based FD, HM, is constructed as following: Four new image blocks g m (i, j) (we set m = 1, 2, 3, 5) (see Fig. 2 (a)-(d)) are obtained from the original image block f (i, j) by the following equation:

3) POWER SPECTRAL DENSITY METHOD
A fractal in spectral domain can be defined by the property that power spectral density [28] has a function form of where ρ is power spectral density in a region within power spectral radius r (see Fig. 3 and Fig. 4).  Taking logarithm on both sides of (16), then we have where C is a constant.
In a log-log plot, it is assumed to be a linear function of spatial scale, with slope 2H − 1. The corresponding PSD fractal dimension D can then be found from the fitting slope.

C. AREA MEASUREMENT APPROACH
The area fractal here is associated with a property of bar graph (see Fig. 5). The AM FD estimate of the ideal fractal surface can be calculated by selecting two scales r 1 and r 2 as where A(r) is visible surface area of all adjunct rectangular boxes with scale r, and D is fractal dimension [24]. The surface of each rectangular box includes a top and four sides. In order to avoid repeated computation, we only need to calculate vertical surface areas of two sides (see Fig. 5(b)).

A(r) can be calculated as
where A H , A v1 and A v2 are the areas of top, right and front sides. f (i, j) is the height (gray level) of a rectangular box with a scale r. M and N are the box numbers in the x, y directions.

III. ROUGHNESS PARAMETERIZATION
In fluid dynamics, the law of the wall states that the average velocity of a turbulent flow at a certain point is proportional to the logarithm of the distance from wall, or the boundary of the fluid region. The concept of an ''equivalent surface roughness'' or dynamical roughness over the ocean is to extend the ''law of the wall'' to air-sea interface for the empirical relation between wind speed (at some reference height) and the net air-sea momentum flux. There are different approaches of modeling roughness length by scaling beginning with Charnock [29] in that roughness length is scaled by the friction velocity. To account for sea state dependency, wave-age scaling (Kitaigorodskii and Volkov 1965, Donelan 1990) and the steepness of the dominant wave scaling (Hsu 1974, and Taylor and Yelland 2001) have been proposed. The steepness of the dominant wave here is commonly defined as the ratio of significant wave height to peak wavelength and wave-age is defined as the ratio of wave phase speed to the friction velocity. For modeling the roughness length changes due to wind fetch or duration, the variations modeled by the wave steepness scaling are smaller, an order of 10%, than by the method of wave-age scaling. Under mixed sea conditions, the wave steepness scaling method is preferred over waveage scaling, while under ''young'' sea conditions, the waveage method is the best performer. Sea-state dependent models (steepness or wave age) are seen to not perform as well as the steepness scaling model in swell-dominated conditions [30].
In the bulk parameterization, the surface fluxes are proportional to the measured mean variables with transfer coefficients to be determined [31], i.e., where ρ and c are the density and isobaric specific heat of air; τ , H and E are the stress, heat and moisture fluxes; U , T and Q are the wind speed, potential temperature and specific humidity at a reference height above sea surface and U s , T s and Q s are the wind speed, temperature and specific humidity at the surface. The momentum transfer coefficient, C D , is also known as the surface drag coefficient. The approximate relationship between C D and u 10 (wind speed at 10m height) has the equation form of: where the values of p, q, r for different ranges of wind speed are different.
In COARE algorithm's parameterization, the roughness length is separated into two parts [6]: where z smooth 0 accounts for roughness of the ocean surface when the surface is aerodynamically smooth, when the surface stress is supported only by viscous stress at surface. The second part z rough 0 accounts for the portion of roughness elements that caused by wave stress due to wind action on the surface gravity waves.
Visible image method only measure small-scale waves due to the limitation on image size and resolution of cameras. It has been found that, at the first order, sea-state parameterizations that based on swell components do not improve much up in estimating surface fluxes [30]. This is because that the large-scale swells may dominate in surface elevation changes, however, to the first-order, it is still the shorter sea waves carry the major part of wind stress. Indeed, the deficiency in roughness length estimations under limited duration or fetch are found to be on an order 10% or less. Thus for open ocean deep water, or pure wind and seas, the momentum transferred from the air to water is at a rate close to that predicted for a fully developed sea regardless of fetch or duration. This has been confirmed by an early study [32].
In general, sea surface roughness and wind speed scaling have mutually dependent relations from various parameterization approaches that have to be solved interactively. Practically, wind speed and sea surface roughness is often fitted into an empirical relationship: For example, from [33], a = 0.15, b = −1.6 and c = 0.366. Here u 10 is the wind speed at 10m height above water level. The particular values of a, b and c may all depend on sea states and is beyond subject of this work. The choice of a particular set of values is not critical in demonstration of our new methodology of measuring sea surface roughness from imagery.
The actual vertical height between our anemometer and sea surface is 5.28m. Conversion the 10m wind from 5.28m wind is by [34]: u 10 = 1.08u 5.28 (24) For this experiment, the actual expression of roughness in term of wind speed becomes: z 0 = 0.15(1.08u 5.28 − 1.6) 2 + 0.366 (25) where the unit of z 0 is 10 −4 m. In reality, there is not a single wind speed-dependent parameterization scheme that can work perfectly well under all conditions for the complicity of the interface coupling systems. For example, the majority of the surface stress is supported by wind-driven waves, while these wind-driven waves are modulated by long waves and water currents. The slopes and amplitudes of wind-driven waves, and therefore the wind stress is modulated by longer surface waves. In practice, the inclusion of additional parameters for coupling mechanisms with their own measurement uncertainties in the bulk flux algorithm tends to increase the uncertainties in estimating fluxes due to difficulties in oceanic observations. Therefore, the potential improvements from the wave age-and wave slope-dependent parameterizations may be better utilized in applications where higher quality wave measurements are available.

IV. EXPERIMENT RESULTS
In this section we present the experimental results on the relationship between wind speeds (or surface roughness) and fractal dimensions, error analysis of the six fractal methods, VOLUME 8, 2020 effects of four noises (Gaussian, salt and pepper, random and multiplicative noises), and transformation analysis.

A. ROUGHNESS ANALYSIS
The data are acquired at a site in Lian island coast of Lianyungang city (China). The site is chosen for the advantages of high fluctuations of both winds and waves, less impurity in seawater and easy access for observation. The observation setup is illustrated by Fig. 6. Cameras on shore and on an unmanned aerial vehicle are used to collect images. A total of 3635 sets of data were collected in this experiment. Each set of data includes one image of sea surface and its corresponding wind speed value at 5.28 meter height. Image data are collected using single camera with a resolution of 3264 × 2468 pixels and wind speed data are obtained by a SMART SENSOR anemometer (AR-816) with an accuracy of 0.1m/s. The measured values of wind speed are between 0∼3.9m/s. The measured vertical distance between the anemometer and the sea level is 5.28 meters. The horizontal length of the field of view is approximately 6 meters and the vertical width is about 5 meters. The aperture of the camera's lens is F1.8 at 2.0× optical magnification.
In the experiments of roughness analysis of sea surface by the six fractal methods (See Fig. 7), the values of fractal dimension are extracted from sea surface images first. Then relations between wind speeds (or sea surface roughness) and fractal dimensions are constructed by using Whisker charts based on the extracted data and equation (25). The maximum value of the upper whisker line is the point at which there are 25% data points whose values lie above the maximum value. The minimum value of the lower whisker line is the point at which there are 75% data points whose values lie above the minimum value. The red solid point is the average value of each set of data. Most FD values of IBC, GVS, HM, PSD and AM methods are between 2 and 3, while only a few FD values of IBC method are between 3 and 6 (See Fig. 7(b)).
As shown in Fig. 7, the mean value of wind speed (or sea surface roughness) increases monotonically with all six extracted fractal dimensions. Mean square error (MSE) and relative error (RE) are used as algorithm evaluation criterion. As shown in TABLE 1, the fitting accuracy of the GVS method is the highest. The above analysis showed that the fractal dimension of visible sea image is proportional to the sea surface roughness in shallow water and the fractal dimension are capable for identifying sea surface roughness changes caused by the wind changes.  Experimental results show that fractal dimensions change with the changes of noise parameters. The changes of fractal dimensions in TABLE 2 and Fig. 9 show that the histogram mean, the improved box counting and the area measurement methods are less affected by noises and the box counting method is more affected by noises. The fractal dimension of the box counting method is calculated only according to the gray values of the four corners of all sub-images so that this method is easily affected by noises.

C. TRANSFORMATION ANALYSIS
Similar images have similar FD values, but images with same FD value may not be exactly the same. In order to distinguish different images of similar FD's, we can check FD of transformed images. We tested on different images constructed by the following four image transformations: where I (x, y) is original image, I 1 (x, y), I 2 (x, y), I 3 (x, y) and I 4 (x, y) are the images after the transformations, L 1 = I min + I m /2, L 2 = I max − I m /2, I max , I min and I m are the highest, lowest and average gray values in I (x, y), (2w + 1) is a scalar and is set to 3. VOLUME 8, 2020 An original image and the corresponding four transformed images are shown in Fig. 10. FD values of seven sets of data (See TABLE 3) are obtained by using the PSD method. As seen in TABLE 3, each set of data has two images with the same FD value, while the FD values are different among data sets. In addition to I 1 , the FD values of I 2 , I 3 and I 4 change distinctly with each other because their pixel gray values have changed and differ from one another after the transformations. So image transformations can be used to distinguish different images with same FD values.

D. OIL SPILL DETECTION
Oil is one of the major pollutants of the marine environment. It may be introduced in diverse ways, such as natural sources, offshore production, sea traffic, tanker accidents, atmospheric deposition, river run off and dumping. Oil spills can seriously affect the marine ecosystem and cause social and scientific concerns since they seriously effect fragile marine and coastal ecosystem. The amount of pollutant discharges and associated effects on the marine environment are important parameters in evaluating sea water quality. Oil spills can easily be  seen but it is difficult, even for an expert, to specify the scope of oil spill at sea [35].
To investigate kinetic characteristics of an oil film at sea, we conducted an experiment west of Point Conception, California on November 3, 2015. We dropped frozen canola oil from a plane flying at an altitude of 500m. The aircraft then flew upwind and downwind over the oil patches at 4 to 5 minute intervals, varying the altitude. We recorded sea surface images with cameras of various angles of view (See Fig. 11). The oil patches formed by the oil are elongated in the direction of the wind, indicating different diffusivities along and across the wind direction.
The wind direction was NWN and the speed was forecast by the National Oceanic and Atmospheric Administration (NOAA) to be 9m/s during our experiment. The closest NOAA windspeed data are from buoys (NDBC 46053 and 46054) located at 34.252N 119.853W and 34.265N 120.477W, which measured winds of average We then calculated that the drifting speed of the oil patches was 4% of wind. The six fractal methods are examined in order to evaluate their performance in surface oil film detection. Fig. 11 contains two aerial images with several oil patches. Six regions in Fig. 11(a) and six in Fig. 11(b) are selected for analysis. In Fig. 11(a), Part 4, 5 and 6 are areas with oil and Parts 1, 2 and 3 are areas without oil. In Fig. 11(b), Part 5 and 6 are areas with oil and Parts 1, 2, 3 and 4 are areas without oil. TABLE 4 shows that the GVS, IBC and PSD methods have the most discriminating ability because their roughness values of oil contaminated parts are quite different from that without oil. The discriminating abilities of the HM method are the worst because their roughness values are not very different. From the above it can be concluded that the GVS, IBC and PSD methods have better performance in oil spill detection than others.

V. CONCLUSION
We have established a new robotic technique for extracting sea surface roughness measurements from visible images. It is based on a novel concept that sea surface random field can be represented by fractal geometry. Specifically, the complicity of random wave image can be represented by the fractal dimension. Our fundamental assumption is that the complicity of the sea surface geometry is linked to the surface VOLUME 8, 2020 drag, and therefore surface roughness. By this connection, a suits of well developed methods, such methods from fractal geometry, are successfully attempted here. Through this study, we have shown: how to construct fractal dimension measurements from sea surface images, what are the empirical relations between wind speeds and fractal dimensions, and how effective of these methods under noises. The algorithms presented in this paper include box counting, fractional Brownian motion and area measurement approaches. Our experiments have demonstrated that our fractal methods can analyze roughness change of sea surface from visible images properly, and the HM and the PSD methods are less affected by noises and the AM and the BC methods are more sensitive to the noises. The GVS, IBC and PSD methods have a better performance in oil spill detection. In future, we are working towards to reconstruct 3D sea surfaces and transparent objects from 2D images [36].