Electromagnetic Scattering and Statistic Analysis of Clutter from Oil Contaminated Sea Surface

In order to investigate the electromagnetic (EM) scattering characteristics of the three dimensional sea surface contaminated by oil, a rigorous numerical method multilevel fast multipole algorithm (MLFMA) is developed to preciously calculate the electromagnetic backscatter from the two-layered oil contaminated sea surface. Illumination window and resistive window are combined together to depress the edge current induced by artificial truncation of the sea surface. By using this combination, the numerical method can get a high efficiency at a less computation cost. The differences between backscatters from clean sea and oil contaminated sea are investigated with respect to various incident angles and sea states. Also, the distribution of the sea clutter is examined for the oil-spilled cases in this paper.


Introduction
Detection of oil spills is a long-term objective in sea remote sensing [1], [2].The oil contaminated sea water is a threat to the marine environment.The vertical polarization of high frequency surface wave radar (HFSWR) has the feature of small energy attenuation, all weather operating and over-horizon.Most importantly, targets can be detected beyond the range of visibility.Thus HFSWR can be used to detect oil spills in far distance.When the sea water is contaminated by oil, the sea dynamics as well as the surface tension are changed [3].Accordingly, the EM scattering characteristics from the oil contaminated sea surface are changed, such as radar echo and sea clutter distributions.The difference between backscatters from clean sea and oil contaminated sea is a useful feature and desirable in Synthetic Aperture Radar (SAR) imagery simulation of ocean scene which is often used to detect oil pollution on the ocean surface.Therefore, how to quantify the scattering differences is an open problem for the radar engineering.
Many researches have been made to examine the EM backscatter from the oil contaminated sea surface, but most of them are limited to two dimensional (2D) sea surface or focus on approximate method [3], [4].However, the hypothetical 2D model lames their applications to three dimensional (3D) cases in practical engineering.Meanwhile, the oil contaminated sea surface is a two-layered dielectric problem.Therefore, the simulation in computational electromagnetics (CEM) becomes much complicate than the homogenous dielectric medium case, and needs further researches.It is our motivation.
In this paper, the electric field integral equation (EFIE) on the interfaces is derived using the wave equation combining the boundary condition.In order to exactly simulate the EM backscatter from the sea surface, the rigorous numerical multilevel fast multipole algorithm (MLFMA) [5], [6] is applied to obtain reliable results.In order to accelerate the EM simulation process when using MLFMA at low-grazing-angle (LGA) cases, a combination model of Thorsos illumination window [7] and resistive loading window [8] is applied to avoid the edge effect induced by artificial truncation of the sea surface in range and azimuth directions.These treatments can greatly improve the computation efficiency thus the simulation of the sea clutter is possible using relatively less computation resources.It is noted that this proposed method is more universal than these existing analytical methods, and unlimited by the 3D sea states or larger incident angles.This rigorous numerical method and its treatments make the Monte Carlo simulation at LGA incidence possible.Then, the electromagnetic backscatter from the oil contaminated sea surfaces is computed with the variation of incident angles and sea states.As a comparison, the backscatter reduction and the statistic characteristics of the two kinds of sea clutter are analyzed compared with the clean sea surface.Some interesting conclusions are summarized in our experiments.

Sea Model
The sea surface can be generated using a spectral method which considers it as a superposition of harmonics.
The amplitudes of the harmonics are proportion to a certain wind-dependent surface-roughness spectrum W(K, ).
Then the sea surface can be obtained by inverse fast Fourier transform (IFFT).In this paper, the 2D Pierson-Moskowitz (PM) spectrum is adopted to generate the 3D sea surface [9].The PM spectrum is defined by where w K is the spatial wavenumber, U is the wind speed at the height of 19.5 m above the sea surface, two constants α = 8.1 ×10 -3 and β = 0.74.The angle φ is measured in the horizontal x-y plane with respect to the x-axis and φ w is the wind direction.In our simulation, the wind direction is towards positive x-axis.
When the sea surface is contaminated by oil, the sea surface tension as well as the sea movement is changed.In the research of Lombardini et al. [10], it is demonstrated that the height of oil contaminated sea surface is damped comparing with the clean sea surface.This damping effect can be expressed by an attenuation coefficient, called Marngoni viscous damping coefficient [11].In this paper, the oil layer thickness of the contaminated sea is 1 mm. Figure 1 shows the comparison of the geometrical profile between the clean sea surface and the contaminated sea surface.The simulation is performed at the wind speed of U = 7 m/s which is considered as the average wind speed over the ocean.It is observed that the surface heights of the contaminated sea are indeed damped comparatively to the clean sea surface.Moreover, the small capillarityscale variations are strongly damped, which implies a strong damping of the surface slopes.

EM Modeling
Clean sea water is homogeneous dielectric medium with a relative permittivity  1 and permeability  1 .However, when the sea water is contaminated by oil whose relative permittivity is  2 and permeability is  2 (generally,  1 =  2 = 1), the contaminated sea surface becomes a twolayered dielectric material with its thickness d, as shown in Fig. 2.
Generally, the HFSW can work at LGA for some extreme cases.Therefore, the multipath effect is complicated for EM simulation.In this section, the rigorous numerical method MLFMA and its treatments are proposed to get the backscatters from the two-layered medium problem.The computation efficiency can be greatly improved, thus makes the Monte Carlo simulation of the sea surface possible.

MLFMA
Using the wave equation and the boundary conditions, EFIE can be derived as functions of the electric equivalent electric current J and equivalent magnetic current M which exist on the interface S 0 between the oil and air as well as the interface S 1 between oil and water.The two sets of surface currents on the boundaries S 0 and S 1 (J 0 , M 0 , J 1 , M 1 ) can be determined using the following numerical simulation procedure.Here the equivalent impedance boundary condition is used to deal with the boundary condition scattering problem [12].It means the electric current J and equivalent magnetic current M satisfy , where n is the outer normal vector at each field point.Therefore, the unknown number in (2) will be halved.The unknown equivalent electric current can be obtained by the Galerkin's method in method of moment (MoM), as following.
To investigate the statistic characteristics of the backscatter of the sea surface, Monte Carlo simulation is necessary.However, it is rather computationally expensive and time-consuming.In order to accelerate the computing process, MLFMA is used to solve the above equation.In the process of MLFMA, the interactions between the elements are classified as near-region and the far-region.The nearregion matrix elements are calculated directly using MoM, while the far-region elements are acquired by using MLFMA.
where NR represents the near-region and FR represents the far-region, l means the basis function, n I are the expansion coefficients needed to solve using MoM, and nm Z are the elements of the impedance matrix.More details can be referenced in [5].

Treatments
Numerical methods can guarantee sufficient precision when solving EM computing problems even at LGA.However, it will bring edge diffraction which is induced by the artificial truncation to the sea surface.The electric current J will change suddenly at the edge of the truncated surface.This phenomenon is called edge effect.In the research of EM scattering in those numerical methods, how to avoid the edge effect is an important issue.The most popular way to suppress the edge effect is using the illumination tapered window (Thorsos window) [7] which is got by modulating a tapered function to incident electric field, as following.
where g is a constant that controls the width of the illumination beam,  i is the incident angle.
However, according to the studies by Jin [13], the length L of rough surface has to satisfy the condition: , where  is the EM wavelength.
Accordingly, at LGA the scale of the computing sea surface is very large.For example, if  i = 85°, then L  932.
And it leads to a large amount of unknown number on the sea surface.Therefore, instead of the illumination window, the resistive loading window is applied to suppress the edge effect under the plane wave illumination.It can be accomplished by simply adding RJ to the EFIE formula, taking (2a) as an example, as following 0 ( ) where R is the resistance, and ( ) L  is the matrix-vector multiply operator and denotes the procedure in left side of (2a).The resistance is loaded within the resistive loading area, which was discussed in detail in reference [14].Once the resistive loading method is used, the length L of sea surface only needs to be as large as 10 times of the correlation length of the rough surface.It is much smaller than the condition of Thorsos window case.However, it is found that the resistive loading will bring the convergence problem because that the impedance matrix nm Z is a non-diagonally dominant matrix for this two-layered dielectric medium scattering problem in our simulations.Since the resistance R will attenuate the value of self-impedance element, the condition number of the final discretized impedance matrix will become worse.Therefore, more computational CPU time will be taken in the impedance matrix iteration processes.
A method combining the two different windows is proposed to accelerate its convergence.Resistive loading window is applied in range direction while a Gaussian illumination window is used in azimuth direction.In azimuth direction, the incident plane wave is modulated by a Gaussian function p y e   instead of the Thorsos window.Meanwhile, an important work in range direction is to extend the original sea surface with smoothly curved sections (for example, we used radius of 10λ) that join to planar sections which are angled 30° down from the horizon.So that all the points on the extension surfaces are shadowed by the original rough surface.The resistive loading of the edges has the advantage of not requiring the surface to be modified at LGA.As a comparison, a 600 m (L)600 m (W) PM sea surface realization at U = 3 m/s was generated and simulated under those system parameters: radar frequency 15 MHz, incident angle 60°, vertical polarization.As shown in Fig. 4, the convergence iteration is 228 steps for the traditional resistive window [8] in range and azimuth directions whereas only 50 steps are needed for the proposed combination of the two different windows.Therefore, the proposed method to suppress the edge effect is more effective and can greatly reduce the CPU time, which is a big improvement in Monte Carlo simulation.

Comparison of Backscatters
By using this proposed method, some experiments are presented in this section.The HFSWR operates at the frequency of 15 MHz, vertical polarization.Monte Carlo simulation is applied to get average radar backscatter.A total number of 1000 sea surfaces are generated for each sea state.The simulation sea area is 600 m long with a width of 600 m.The simulation is performed when wind speed is 3 m/s and 9 m/s, corresponding to Douglas sea state 2 and 4.  As shown in Fig. 5, the radar backscatters of the contaminated sea surface is much lower than that of the clean sea surface.The reduction of the backscatter coefficient corresponding to different sea states is shown in Fig. 6.A reduction of more than 10 dB at moderate incident angles is found through our simulation, and even 15 dB at LGA angle.In SAR systems, the oil polluted area is seen as "dark" sea.Capillary waves on the sea surface reflect electromagnetic energy.However, when the oil layer on the sea surface, the capillary waves are dampened.As a consequence, the backscatter is much weaker than clean sea surface.

Sea Clutter Distribution
By using the radar echo data acquired in the above experiments, the distribution of the sea clutter is also examined.For radar system, statistical models (Rayleigh, Lognormal, Weibull and K-distributions) are often used to describe the sea clutter.In this paper, the probability density functions (PDF) of the backscatter of the sea clutter for clean sea and oil-contaminated sea are studied.The parameters of Rayleigh, Lognormal and Weibull distribution are estimated by the maximum likelihood (ML) method.It is a standard approach to parameter estimation and provides optimum estimates in the sense that these estimates are the most probable parameter values.However, ML is computationally expensive to estimate the parameters of K distribution.A new method based on higher order and fractional moments is proposed in these references [15], [16].It can significantly reduce the computational requirement.Details about the formulas of these four distributions and the estimation methods are summarized in reference [16].Take this sea state of U = 9 m/s as an example, the sea clutter distribution is shown for different incident angles in Fig. 7. max ( ) ( ) max / ( )  As can be seen from Tab. 2 and Tab. 3, the distribution of clutter from oil contaminated sea surface is different from that of clean sea surface.For most cases, the sea clutter satisfies K distribution and Weibull distribution.When the wind speed above the sea surface is moderate (3 m/s), K distribution describes the sea clutter quite well at small incident angle.As the incident angle getting large, it is Weibull distribution for the contaminated sea whereas it is Lognormal distribution for the clean sea.When the wind speed gets higher (9 m/s), Weibull distribution will be dominant at LGA.The distribution of clutter from oil contaminated sea is first mentioned in this paper.Rigorous numerical EM simulation method can ensure sufficient accuracy.Moreover, with the treatments this method can obtain a high efficiency.Consequently, Monte Carlo simulation can be performed on the electric large problem.The real radar data from HFWSR system is usually very expensive.Instead of the real data from the radar, the simulation results by using rigorous numerical computational electromagnetic method as shown in this paper is workable as well.Moreover, it is more cost-effective and available.

Conclusions
In this paper, the sea clutter distribution is examined when there is oil pollution on the sea surface.In HFSWR system, the observation area is very large; accordingly the accurate EM simulation of the large sea surface is rather difficult.The integral equations both on the clean sea surface and the oil contaminated sea surface are presented.Numerical method with treatments on it can ensure suffi-cient accuracy to get reliable results of the backscatter from sea surface.Meanwhile, the computational resources are saved a lot, and then Monte Carlo simulation can be performed.It is found that the backscatter from oil contaminated sea surface is much weaker than that from clean sea surface.That is because the oil damped the capillary waves.Accordingly, the reflected electromagnetic energy is weaker than before.This effect can be seen in SAR image that the oil contaminated area is dark.The distribution of the sea clutter which is an important statistic characteristic is investigated by using the simulated results.Differences of the sea clutter distribution are found between the oil contaminated sea surface and the clean sea surface.

Fig. 1 .
Fig. 1.Comparison of clean sea surface and contaminated sea surface.
d) where 0 S ds   denotes principle integration, n Z is the characteristic impedance, n k is the wave number and Green function.The subscript n = 1, 2, 3 denotes the three cases in free space, sea water and oil medium, respectively.

Fig. 7 .
Fig. 7.The backscatter distribution of sea clutter.empirical distribution functions of two samples.The K-S test statistic is defined as De which is shown

F
t (x i ) is the cumulative distribution of the distribution being tested and F (x i ) is the cumulative distribution of the simulated sea clutter.De means the maximum deviation between the cumulative distribution of the tested distribution and the statistic distribution of the simulation result.The best fit distribution is with the smallest De .Tab. 1 is the K-S test of the sea clutter distribution corresponding to Fig. 7.And the sea clutter distributions are also examined under different sea states in Tab. 2 and Tab. 3.
The distribution of sea clutter for contaminated sea.
Tab. 3. The distribution of sea clutter for clean sea.