Oil Spill Detection in Glint-Contaminated Near-Infrared MODIS Imagery

We present a methodology to detect oil spills using MODIS near-infrared sun glittered radiance imagery. The methodology was developed by using a set of seven MODIS images (training dataset) and validated using four other images (validation dataset). The method is based on the ratio image R = L'GN/LGN, where L'GN is the MODIS-retrieved normalized sun glint radiance image and LGN the same quantity, but obtained from the Cox and Munk isotropic (independent of wind direction) sun glint model. We show that in the R image, while clean water pixel values tend to one, oil spills stand out as anomalies. Moreover, we provide a criterion to distinguish between positive and negative oil-water contrast. A pixel in an R image is classified as a potential oil spill or water via a variable threshold Rs as a function of L'GN, where the threshold values are obtained from the slicks of our training dataset. Two different fitting curves are provided for Rs, according to the contrast sign. The selection of the correct fitting curve is based on the contrast type, resulting from the criterion above. Results indicate that the thresholding is able to isolate the spills and that the spills of the validation dataset are successfully detected. Spurious look-alike features, such as clouds, and other non-spill features, e.g., large areas at the glint region border, are also detected as oil spills and must be eliminated. We believe that our methodology represents a novel and promising, though preliminary, approach towards automatic oil spill detection in optical


Introduction
Marine oil spill pollution, apart from rare accidents, such as ships in distress or the rupture of pipelines, is a frequent occurrence, mainly coming from illegal discharges of hydrocarbons that are intentionally released by ships into the marine environment.The capability of satellite remote sensing in the detection and monitoring of oil spills at sea has been largely demonstrated [1][2][3][4][5].Oil spill remote sensing techniques are based on visible, infrared, microwave and radar sensors.Among all satellite sensors, Synthetic Aperture Radar (SAR) is "the most applicable space-borne sensor for operational oil spill detection" [1], mainly because of its high spatial resolution and all-weather and all-day capabilities, that is, in observing the sea at night and in cloudy conditions.
SAR oil spill detection relies on the fact that an oil film has a dampening effect on the capillary waves and decreases the sea surface backscattering, resulting in a dark region that contrasts with the brightness of the surrounding spill-free sea.SAR detection has however some limitations, such as the conditions outside the optimal wind speed range for detection (i.e., 3-10 m/s) or the presence of natural films (i.e., phytoplankton, rain cells and freshwater slicks), that can dampen the short waves and create dark patches on the surface.These factors thus yield false positives (i.e., look-alikes), which cannot be discriminated from true spills by using the sole SAR data.Furthermore, SAR data are expensive, not available daily (long revisit time, several days) and cover small areas, due to their reduced swath width, i.e., 100-500 km.These limitations call for a multi-platform SAR slick detection effort, as well as for the additional use of satellites with optical sensors.The latter, with their daily global coverage, wider swaths (greater than 1000 km) and, in many cases, free-of-charge data, can complement and optimize SAR detection, also in discriminating between oil spills and look-alikes via the retrieval of their optical properties [6].Furthermore, a multi-platform observation increases oil spill detection effectiveness, since more observational data can be fed to slick displacement and transformation forecast models, which are the object of recent research efforts [7][8][9] that seek to integrate observation and forecasting in detection systems.
The first attempts to detect oil spills with optical sensors [10] are more recent with respect to SAR.This is due to the increased spatial resolution of the new generation optical sensors, i.e., MODIS (Moderate Resolution Imaging Spectroradiometer) onboard the TERRA and AQUA platforms and MERIS (Medium Resolution Imaging Spectrometer) onboard the ENVISAT satellite, the latter having ceased to function in May 2012.Indeed, the 250-m spatial resolution of some MODIS optical bands is now sufficient to resolve typical illegal discharge slicks of 1-10 km with sufficient detail, while the typical 1-km resolution of past sensors, such as the Advanced Very High Resolution Radiometer (AVHRR) sensor, was practically useless in slick detection.Several studies exploiting oil spill optical properties in MODIS data have been carried out by Hu et al. [11], Bulgarelli and Djavidnia [12] and Zhao et al. [6].However, in contrast to SAR detection, the effectiveness of which is well documented, the literature on optical oil spill retrieval is still sparse, since there is a lack of oil spill detection methodologies, if we exclude first attempts, such as those proposed by Adamo et al. [13] and [6].In particular, the work of [6] is of greater interest and provides a comprehensive analysis of oil spill spectral properties for a set of oil spill cases observed in the Arabian Gulf, both under glint and glint-free conditions.In particular, [6] devises a methodology to distinguish between oil spills and algal blooms, based on the use of the FAI (Floating Algae Index) technique.However, the oil spill detection methodology of [6] is developed and applied only for oil spill cases in the Arabian Gulf.
The lack of methodologies in optical oil spill detection is due to the difficulties emerging from the MODIS dependence on the illumination-view geometry, atmospheric conditions and cloudiness.Furthermore, a comprehensive database of "certified" oil spills, i.e., observed in situ by naval or aerial surveys, is not yet available.
In general, from the point of view of optical remote sensing, the visibility of an oil spill is due to the optical contrast between the spill and surrounding water.Oil-water contrast strictly depends on: (1) the illumination-view geometry, viz. the satellite and solar zenith angles and the relative satellite-solar azimuth angle; (2) sea state, i.e., wind speed; and (3) oil spectral properties, i.e., refractive index and absorption coefficient [14], and oil conditions, i.e., the thickness and state (e.g., fragmentation, evaporation, etc.).Summarizing Otremba and Piskozub [15][16][17][18], the remotely-sensed optical contrast of an oil-covered sea surface increases either when the sensor sees the oil slick close to the sun reflection angle, that is in glint conditions (positive contrast, oil brighter than water), or when the slick is observed at a low angle (no glint), almost horizontally (negative contrast, oil darker than water).Furthermore, in "moderate" glint regions (located at the edge of the glint "patch"), the contrast of slicks can change from positive to negative and vice versa.Finally, contrast is expected to be always positive for oil emulsions [17].Thus, parameterization of the change in contrast sign as a function of the Sun-satellite geometry is needed for a good oil spill detection algorithm, i.e., to know what contrast is to be expected.
State-of-the-art optical oil spill detection thus tends to distinguish between glint and glint-free conditions, because the mechanism of oil spill visibility is different in the two cases.Indeed, in glint conditions, which occur when the solar beam is directly reflected by the sea surface to the sensor, an oil spill can be detected as a sea surface roughness anomaly, in a similar way to the SAR detection mechanism.In glint-free conditions, instead, an oil spill can be detected only via its optical properties, i.e., differences in the spectral signatures with respect to water [6,11].
In the first case, the wave dampening produced by an oil film enhances the reflected radiance with respect to the surrounding water.Hence, the oil film itself, by acting like a mirror, results as a brighter feature in an optical image.In the latter case, a floating oil film blocks or absorbs the water-leaving radiance from below and absorbs the incoming surface radiance (the oil absorption coefficient is greater than water).This modulation in the total upward radiance can give rise to an optical contrast between oil and surrounding water, thus making an oil film visible in an optical image [12] as a dark feature.
The presence of sun glint contamination is the most favorable condition in which an oil spill can be detected, because the mirror effect of oil on water greatly aids the specular reflection (the oil refractive index is greater than the water index).We reiterate, however, that the detectability via optical properties (i.e., oil spectral signatures) is not possible in glint regions, and that spills can be detected just as sea surface roughness anomalies.This is because the sun glint radiance, without ever entering the water, is so large compared to the entire water-leaving signal from the euphotic zone that the retrieval of bio-physical products, such as chlorophyll concentration, is compromised [19].In general, the atmospheric correction procedure tends to flag the glint-contaminated regions of the images, so no optical products are retrievable (besides TOA radiances and products to which an atmospheric correction is not applied) in such conditions.However, the use of TOA radiances "as is" for oil spill detection can be inefficient, in that the oil-water contrast depends on atmospheric conditions, which may either modulate and, in the limit, hide the true spill surface signal or generate atmospheric oil spill "look-alike" features (e.g., aerosols) in the image.Therefore, a custom atmospheric correction that allows one to process glint-contaminated images must be applied.
The aim of this work is to develop a methodology for oil spill detection in MODIS near-infrared (859 nm), high resolution (250 m), sun glittered, top of atmosphere optical imagery.The idea is based on the computation of the normalized sun glint radiance, i.e., the solar radiance directly transmitted to the sensor due to specular reflection as if there were no atmosphere, following two different procedures: (1) using MODIS TOA radiances; and (2) using the Cox and Munk (1954) isotropic (wind-direction-independent) sun glint model [20].To obtain the normalized sun glint radiance ("glint radiance") from MODIS TOA radiances, it is necessary to first perform an ad hoc atmospheric correction, which does not flag glint regions, in order to remove the effects of Rayleigh scattering and aerosols.Once the two glint radiances (respectively L'GN and LGN) are computed from the two procedures listed above, by taking their ratio R = L'GN/LGN, we expect that: (1) the ratio image R will show a reduction on radiance variability, also due to the removal of atmospheric, natural variability from the TOA reflectance input images; (2) the ratio image will be independent from the solar-view geometry and atmospheric conditions; and (3) oil spills and other surface features that dampen the sea surface roughness will be detected as anomalies on the ratio image, since the Cox and Munk (1954) model glint radiance is always computed for a clean (i.e., oil-free) sea, while the experimental glint radiance obtained from the TOA radiance perceives the oil spill as "unnatural".
Since positive and negative oil-water contrast is possible in glint-contaminated images, we also develop an auxiliary procedure to characterize contrast type by means of threshold values in the illumination-view geometry parameters and in the glint radiance of the image.Finally, any given pixel in an R image will be classified as a potential oil spill or water via a variable threshold Rs function of L'GN in each pixel.
The method (Section 3) was developed using a set of seven MODIS images ("training dataset") and tested using another four images ("validation dataset").All of the images are relative to an oil spill database that is built from certified oil spill events (Section 2).Although we make use of a small number of images (a fact that makes this effort preliminary), our database is comprehensive for positive and negative oil-water contrast cases.The results (Section 4) indicate that oil spill pixels are successfully detected by comparing R for each pixel with the Rs(L'GN) threshold curve, built using R and L'GN of the known slicks of our training dataset.Spurious look-alike features, such as clouds, are also detected and must be eliminated.Conclusions are given in Section 5.

The Oil Spill Database
By collecting the available information on oil spill disasters, such as the Lebanon oil spill (Mediterranean Sea) occurring during the 2006 bombing [12,21] and the Deepwater Horizon drilling rig explosion (Gulf of Mexico, 2010) [8,[22][23][24], and lesser events (e.g., spills released by ships) occurring in the Arabian Gulf [6] and in the Mediterranean Sea [13], we built a fairly large database of in situ oil spill observations.For each of these events, we acquired the corresponding MODIS (AQUA and TERRA) Level 1A passes and retained those in which the spill was well visible and not covered by clouds.This brought us to finally consider eleven MODIS-detected oil spills (Table 1).Seven oil spills of this final database have been selected to build the algorithm (training oil spill database), while the remaining four were used to validate the methodology (validation oil spill database).Oil spill events (Table 1) include a temporal suite of 4 images of the Lebanon 2006 coastal spill event (Scenes 1-4), an oil spill south of the island of Crete (Eastern Mediterranean Sea, Scene 6), 5 images relative to the Deepwater Horizon 2010 disaster (Scenes 5, 7-10) and several slicks detected in the Arabian Gulf [6] (Scene 11).All cases include glint conditions, as defined in the SeaDAS processing software (SeaWiFS Data Analysis System, distributed by the NASA Ocean Biology Processing Group).
The slicks in the MODIS TOA radiances at 250-m resolution of the oil spill dataset were manually digitized using the ENVI Region of Interest (ROI) tool, in order to produce oil spill position text files (Figure 1).The criteria to visually identify and digitalize the slicks that compose the certified oil spills in the images were: (1) presence of similar slicks in companion SAR imagery when available; (2) an elongated/winding shape; and (3) spills present in the literature cited above.The slicks were digitalized in a "conservative" way, in that minor or uncertain features and pixels apparently belonging to a slick, but with a qualitatively too low contrast with the surrounding water, were not considered.Nevertheless, since it is impossible to certify every single slick of a spill by visiting it with a vessel and companion SAR imagery may also display look-alikes, some uncertainty is always present in the digitized results.1): oil spills are brighter than the surrounding water (positive oil-water contrast); (b) the same image with manually-digitized oil spill pixels (ROIs, green regions).The black area indicates the swath's border.
Table 1 also indicates the oil-water contrast sign, i.e., whether the contrast of the spills with respect to the surrounding water is positive or negative.We reiterate that in moderate glint areas, both positive and negative contrast can be found in the same MODIS image (Table 1, Scene 7).This distinction was necessary, because we wanted to parameterize the contrast sign as a function of the illumination-viewing geometry, in order to treat differently the positive and negative contrast cases in our procedure.Finally, this oil spill database can be considered as our initial sample for global ocean oil spill cases and must be updated as new oil spill cases become known or historical, certified oil spills are included in future research works.

Modis Data Products
MODIS Level 1A data, corresponding to the oil spill database described above (Table 1), were downloaded from the NASA website [25] and processed to Level 2 using the standard NASA ocean color data processing, SeaDAS, v. 6.4.Only the full-resolution 250-m products at λ = 859 nm wavelength were retained, because this band shows the best contrast between oil slicks and surrounding clean sea water [11][12][13] and is one of the two MODIS bands available at 250-m spatial resolution.
We reiterate that TOA radiance Lt is the radiance received by the sensor containing both the atmospheric and surface signal; Rayleigh radiance Lr and aerosol radiance La are, respectively, the radiances due to molecular and aerosol scattering [26], overlying a Fresnel-reflecting ocean.In particular, in SeaDAS, the Rayleigh look-up-tables (LUT) are built for an atmosphere bounded by a rough water surface, while the LUT for the aerosol selection are built for an atmosphere bounded by a flat sea surface.The normalized sun glint product LGN is the radiance due to the specular reflection of the direct solar beam, as estimated by the Cox and Munk (1954) isotropic model.Wind speeds are from the National Centers for Environmental Prediction (NCEP).These parameters are used as the input dataset for the detection methodology described in Section 3.
We can define the glint angle α using the Sun-satellite geometry products [11] as: The glint angle α is the angle between the satellite viewing direction and the direction of the specular reflection, with θ0, θ and Δϕ being the solar zenith (solz), sensor zenith (senz) and the Sun-satellite relative azimuth angle, respectively.α contains the full information about the solar-sensor geometry and will be discussed in Section 4.

Background
The total radiance at the top of the ocean-atmosphere system can be written as (Gordon and Wang [26]): where Lr(λ) is the radiance resulting from multiple scattering by air molecules (Rayleigh scattering) in the absence of aerosols, La(λ) is the radiance resulting from multiple scattering of aerosols in the absence of air and Lra(λ) is the interaction term between molecular and aerosol scattering [26].Lg(λ) is the glint radiance due to specular reflection of the direct solar beam, while Lwc(λ) and Lw(λ) are the whitecap and water-leaving radiances.Finally, T(λ) and t(λ) are, respectively, the direct and diffuse atmospheric transmittances in the sensor-viewing direction.Note that in the near-infrared bands (λ = 859 nm), we can neglect the water-leaving contributions, at least in the Case I (open ocean, oligotrophic) water.Indeed, in coastal (Case II) turbid waters, Lw(λ) at 859 nm are not negligible, due to inorganic suspended sediment or dissolved organic matter contributions.However, in glint-contaminated situations, this contribution in Case II waters is still negligible with respect to the glint radiance value [13].Furthermore, whitecaps, i.e., Lwc(λ), may be neglected, in that even at high wind speeds, the maximum contribution is about 0.5% of the TOA radiance and, thus, of the glint radiance, which is the main component of Lt in glint conditions [27].We hereafter drop λ in the following equations, since we will deal exclusively with 859-nm products.
If we include the Rayleigh-aerosol interaction Lra into a unique aerosol contribution LA = La + Lra [28], we can write Lt as: The sun glint radiance Lg can be therefore expressed as: where F0 and T0 are, respectively, the extraterrestrial solar irradiance and the atmospheric direct transmittance along the solar direction.
LGN is the normalized sun glint radiance, that is the glitter radiance just above the sea surface as if there were no atmosphere and with the extraterrestrial solar irradiance F0 = 1.The standard ocean color data processing computes the normalized sun glint radiance based on the Cox and Munk (1954) formulation that represents the sea surface perturbed by the wind as a collection of tilted facets, the slope statistics of which is derived from a probability distribution function.SeaDAS v6.4 adopts the isotropic version of the Cox and Munk (1954) model, which is independent of wind direction.The Cox and Munk (1954) normalized sun glint radiance can be expressed as: where r(ω) is the Fresnel reflection coefficient, ω the specular reflection angle, P(β,σ) the probability distribution of sea surface facet slopes as a function of zenith angle β of each facet and σ is an isotropic sea surface slope coefficient, dependent on wind velocity W (m•s −1 ).In practice, LGN depends on solar zenith θ0, sensor zenith θ, the relative Sun-satellite azimuth angle Δϕ and surface wind velocity W, because also β depends on the above angles [28], i.e., LGN = LGN (θ0, θ, Δϕ, W).The solar-sensor geometry (θ0, θ, Δϕ) is provided by MODIS, while the wind velocity W data by NCEP.
In SeaDAS, the normalized glint radiance is used as a mask, so pixels with LGN ≥ 0.005 are flagged, and no further processing is done on MODIS data above this LGN threshold.In areas where 0 < LGN < 0.005, the SeaDAS processing performs a glint correction in order to remove the effects of glint contamination and, thus, to obtain the desired products.Furthermore, for some authors, the SeaDAS threshold value of 0.005 can be stretched to 0.01 [29], thus increasing the correctable number of pixels in these glint-contaminated areas.However, we maintain the 0.005 threshold in order not to have the glint correction on pixels with 0.005 < LGN < 0.01 and, thus, to treat all pixels with LGN > 0.005.

Oil Spill Detection Strategy
By using Equations ( 3) and ( 4), we can express the normalized sun glint radiance as: with: where τr and τa are the Rayleigh and aerosol optical thicknesses and θ0 and θ are the solar and sensor zenith angles, respectively.We can define L'GN as the MODIS-measured normalized sun glint radiance, which represents the "true" counterpart of the Cox and Munk (1954) model LGN.Comparison studies between various existing glint models and L'GN [28] showed that the two (isotropic and anisotropic) Cox and Munk (1954) models are the most accurate for glint evaluation in MODIS images.In most cases, the isotropic approximation to the Cox-Munk model is sufficiently accurate.
Our methodology for oil spill detection is based on the computation of the ratio: From the definition of R, we expect values around 1 for pixels belonging to clean sea water and deviations from this value for oil spill pixels.Furthermore, the ratio R is much less dependent on viewing geometry than the glint radiance values of its components, in that this dependence in L'GN and LGN somewhat cancel each other in the ratio.Furthermore, R is independent from atmospheric aerosol and Rayleigh scattering, since L'GN and LGN are products derived at the sea surface.
To obtain the ratio image R, the MODIS-derived sun glint radiance L'GN is needed, since this quantity is not provided by SeaDAS, contrary to the Cox and Munk (1954) isotropic glint radiance LGN.

MODIS Sun Glint Radiance Computation (L'GN)
To obtain L'GN (Equation ( 6)), we need the aerosol radiance LA and the aerosol optical thickness τa in glint regions, since the SeaDAS processing flags glint regions and, thus, does not provide them.However, following Zhang and Wang [28], we can assume that the aerosol type and its optical thickness vary only slightly between glint and glint-free areas, since the aerosol variability scale is modulated by atmospheric systems characterized by a space scale of O (100-1000 km).Thus, we can determine a dominant aerosol model and its average <τa> over the glint-free regions and then compute LA over the glint area using such an aerosol model and this <τa>.
To obtain LA for each image of the dataset, we first found the dominant aerosol radiance models over each region [28], the values of which encompass the observed TOA radiance ("minimum" and "maximum" models), as defined by the SeaDAS 6.4 atmospheric correction procedure, in regions where glint is negligible (LGN < 10 −7 ).We reiterate that SeaDAS, since version 6.0, uses 80 aerosol model types [30] (8 humidity profiles × 10 aerosol size fractions).We then computed the average model ratio and the average aerosol optical thickness <τa> corresponding to the above aerosol models.These results are shown in Table 2 for λ = 859 nm.Thus, the minimum and maximum dominant aerosol models, the average model ratio and <τa> become the fixed inputs to the atmospheric correction model with which we can compute LA over the glint regions in each scene.Since SeaDAS does not include an atmospheric correction algorithm that uses fixed values of both τa and minimum and maximum models, we modified the SeaDAS code to accommodate such fixed parameters.Finally, T0 and T (Equation ( 7)) are computed by using the value <τa>.This allows for the computation of L'GN (Equation ( 6)).In rare cases, as in Scene 5 (Table 1), few pixels scattered in the image have L'GN < 0 (failure of atmospheric correction).A possible way to overcome this problem is to apply a "white aerosol approximation", which assumes a spectrally flat aerosol reflectance for the atmospheric correction procedure.For simplicity, in this work, we masked pixels with L'GN < 0.
To verify the consistency between L'GN and LGN, we computed the mean bias (MBE) and standard deviation of the L'GN − LGN differences in each image, including all pixels, i.e., also oil slicks, clouds and clean seawater.This comparison revealed a negative mean bias in almost all images, i.e., with L'GN slightly less than LGN.The MBEs are found to be about 10 −3 -10 −2 , i.e., 10%-15% of the typical LGN and L'GN around 0.05 (Figure 2c).Therefore, these MBEs may well be due to uncertainties in the two above glint radiance estimations.Now, since we want the ratio R = L'GN/LGN to be as close as possible to 1 for clean sea water pixels, in order to detect slicks as deviations from this value, we systematically subtract this bias from each L'GN of the dataset before computing the ratio L'GN/LGN (Equation ( 8)).Hereafter, we indicate with the same symbol L'GN the bias-corrected glint radiances, i.e., L'GN = L'GN (uncorrected) − MBE.

Evaluation of LGN, L'GN and R
We computed LGN, L'GN and the ratio images R = L'GN/LGN for all of the scenes of the dataset (Table 1).Figures 2-4 show the results of LGN, L'GN and R relative to Scenes 2 (Lebanon), 6 (Crete) and 7 (Gulf of Mexico) of the dataset.The known spills have positive (Figure 2), negative (Figure 3) and both positive and negative (Figure 4) contrast with respect to the surrounding clean sea water.A comparison between glint radiances LGN and L'GN is shown both along a transect crossing a portion of the spill and with histograms covering the full image.In the best cases (Figures 2 and 3), L'GN and LGN tend to coincide in clean sea water pixels (black line and dots in Figures 2d and 3d), giving values of R close to one (Figures 2f and 3f), while they differ (R ≠ 1) in oil spill pixels (Figures 2d and 3d, green dots) and a stronger separation occurs for cloudy pixels (high-value isolated black dots in the same figures).The overall behavior of LGN and L'GN is shown by their histograms (Figures 2e and 3e), and the match between populations is quite good, as expected by the theory.
No cloud masking was adopted in the image processing.Indeed, trials in masking clouds revealed that, in glint conditions, many oil spill pixels were flagged together with true cloud pixels, since oil and cloud radiances are both similarly high.We tried to use both the SeaDAS cloud mask flag and cloud masking procedure proposed by Wang and Shi [31] based on a reflectance thresholding using different NIR bands.Since both procedures flagged slick pixels along with clouds, we preferred to keep clouds, and we intend to develop an ad hoc cloud masking algorithm.3.
For a few cases of the dataset and only for some pixels of the images, L'GN and LGN showed strong differences also for clean sea pixels, thus yielding R ≠ 1 (see Figure 4c,d at the transect ends), with strong fluctuations of R around the reference value of one (Figure 4e,f).In particular, the Cox and Munk (1954) glint radiance LGN goes to zero in the western (left) portion of the image (Figure 4c,d), as the no-glint area east of the image scene is approached, but where LGN is still greater than 0.005.On the other hand, the experimental glint radiance L'GN is still significantly nonzero there.This behavior is a limit for the procedure in its present version, in that there might be the chance to classify as probable oil also evidently non-spill regions.This is probably due to the fact that both L'GN and LGN are derived under somewhat "strong" assumptions, i.e., a constant aerosol model and constant aerosol optical thickness for the former and a wind-direction-independent model (i.e., Cox and Munk, 1954) for the latter.Furthermore, [32] shows that, when glint radiance becomes less predominant at the edge of the glint zone, multiple scattering becomes dominant and should be taken into account to improve the estimate of model glint radiance.Another improvement can be possibly obtained by computing LGN with the Cox and Munk (1954) anisotropic (wind-direction-dependent) model.3.
The ratio image R, being almost independent of the atmospheric conditions and solar-viewing geometry, shows a reduction of the high spatial variability, usually present in TOA Lt products due to the atmospheric signals.This reduction in radiance variability results in a more uniform and flattened image, i.e., the slick-clean water contrast is enhanced in R. Indeed, if we compute the mean values and standard deviations of the digitized oil spills over the TOA radiances Lt for all of the scenes with positive contrast of the training dataset, we obtain a total Lt range of [2.50, 15.00] mw•cm −2 •µm −1 •sr −1 for the oil spill radiances.This interval includes the standard deviations at its extremes.Thus, comparing this range with the equivalent R range [0.85, 2.25], it is evident that the relative variation is of the order of 83% and 62%, respectively, indicating a successful reduction of the variability in the R images.3.

Oil Spill Separation Index
The fact that oil spills are often well separated from clean sea water in R images (e.g., Figure 2f) led us to adopt an R thresholding method for oil spill detection.For each scene of the training dataset with positive contrast, we defined an oil spill separation index, Rs, as the value above which lies approximately 70% of the R values of the oil spill ROI pixels.Then, if we highlight all of the pixels with R ≥ Rs in the same R image, we obtain a map of probable spills with positive contrast (Figure 5a, red features).For negative contrast, Rs is the value below which we find approximately 70% of the ROI R values.Thus, highlighting the pixels with R < Rs, we obtain a map of probable spills characterized by negative contrast (Figure 5b, blue features).An example of both the positive and negative oil spill contrast case is shown in Figure 5c.What is evident from Figure 5 is that also a subset of non-oil spill pixels are classified as probable oil spills using the Rs thresholds.This is mainly due to the fact that some pixels, such as cloud or clean water pixels, have R values that differ from one and thus fall in the oil spill R range (particularly evident in Figure 5a,c).Furthermore, in Figure 5c, a large area highlighted in red is particularly evident, which can easily be eliminated as a look-alike because of its extension.This occurs because of the vicinity to the no-glint area, as described in Section 4.1.
In order to limit look-alike occurrence, we have chosen the above 70% threshold to determine Rs, even though this misses some oil spill pixels of the training ROIs.This choice proves to be a good compromise, as for the results from the validation effort described below.
The Rs index values found for each scene of the training dataset are listed in Table 3.We found six Rs index values for the positive contrast cases and two for the negative ones.The separation indices are also plotted against the corresponding average L'GN values, i.e., the mean L'GN value for which R ≥ Rs (R ≤ Rs) for positive (negative) contrast cases (Figure 6, Table 3).A cubic spline was fitted to obtain an analytical relationship for positive contrast cases (Figure 6, red curve), while a linear relationship was obtained from the two negative contrast cases (Figure 6, blue line).Note that the blue line for the negative contrast cases applies only for L'GN < 0.03, so the two data points with which it was obtained at least are reasonably representative of the line.Nevertheless, novel additions of known oil spill cases to the training dataset are needed and will allow a more precise determination of these curves.  1 and 3).The red and blue curves are the cubic spline and linear interpolation, respectively, between Rs values.The red (blue) curve represents the lower (upper) R limit for oil spill detection in positive (negative) contrast situations.
It is important to note that the Rs data relative to the red threshold curve (positive contrast, Figure 6) have a smooth dependence with the sun glint radiance L'GN, i.e., they lie nicely along the fitted curve.This behavior does not hold if we look at maximum values of R in the training dataset oil spill ROIs.Indeed, for any fixed L'GN value in a given ROI, the corresponding R maximum values are highly variable with respect to any curve fit.This probably depends on oil thickness, e.g., because an increasing slick thickness enhances the damping of surface waves, thus increasing contrast in glint conditions.Therefore, we adopted only one Rs vs. the L'GN threshold curve (Figure 6), above which we classify pixels as oil, with no upper R limit.The inverse holds for negative contrast cases (blue curve, Figure 6).
In this way, when we have a new R image to inspect, we derive Rs for all L'GN values in the image using our curves (Figure 6) and highlight as probable oil spills (as in Figure 5 above and Figures 7 and 8 below) those pixels with R values above (positive contrast, in red) or below (negative contrast, in blue) our Rs vs. L'GN curves.The only problem to solve is to determine which curve to use as the detection limit, i.e., to determine the sub-regions of the image for which we expect oil spills with positive or negative contrast.For this purpose, we can use the glint angle α (Equation ( 1)) or the Cox and Munk (1954) glint radiance LGN.In fact, these parameters have already been used as a possible threshold to separate positive from negative oil-water contrast [11].Our trials (Table 3) and those of Hu et al. [11] indicate that positive contrast regions are those for which either α < 12° or LGN > 0.05, while negative contrast is to be expected for α > 17°.All other regions (those with 12° < α < 17° or LGN < 0.05) are to be regarded as possibly hosting both positive and negative contrast oil spills.The thresholding of α is also close to the results of Bulgarelli and Djavidnia [12], who found one Lebanon case in which contrast switched from negative to positive for α = 18.31°, i.e., very close to our mixed contrast interval.Table 3. Set of parameters used to characterize oil spills in glint-contaminated MODIS imagery.Rs are the lower (positive contrast spills) and upper (negative contrast) limit R values for detection (curves in Figure 6).All other values are mean values computed over the oil spill ROIs, where R ≥ Rs (R ≤ Rs) for positive (negative) contrast cases (70% of ROI pixels).Note that we split the last oil spill case (No. 7) into two rows, that is for positive and negative contrast.The glint angle (or glint radiance LGN) criterion allows us to apply one of the curves in Figure 6 for oil spill detection in R images for which either positive or negative contrast is expected, while in regions with mixed contrast, both curves are applied, thus highlighting as probable oil spills pixels with R below the blue curve, as well as those above the red curve of Figure 6.

Application
Our oil spill detection methodology has been applied to the R images of our four validation cases (Scenes 8-11 of Table 1) relative to the Deepwater Horizon oil spill accident of April, 2010, in the Gulf of Mexico, well described by Liu et al. [8,22], and several slicks detected in the Arabian Gulf [6].We understand that the number of validation scenes is too low for significance and that they belong to the same (Deepwater Horizon) spill, with the exception of Scene 11 relative to the Arabian Gulf.Nevertheless, different acquisition dates and times for the same spill event mean different solar-sensor viewing geometry, sea state, wind and atmospheric conditions, oil evaporation and dispersion conditions.All of this gives representativeness even to the validation of our detection procedure, which represents, to our knowledge, a novel step towards automatic oil spill detection.
The results discussed here (Figures 7 and 8) were obtained using our methodology, which we briefly summarize for convenience.For each pixel in a given glint-contaminated image, the procedure: (1) extracts glint radiance L'GN, R and glint angle α; (2) finds the location of L'GN on the x-axis of Figure 6; (3) finds the threshold Rs value according to the red (blue) curve of Figure 6 if α < 12° (α > 17°), since positive (negative) oil-water contrast is expected; for 12° < α < 17°, the two Rs values, relative to both curves, are considered; (4) compares R with the Rs(L'GN) value for the proper contrast case; two comparisons are made for 12° < α < 17°; (5) if R > Rs (positive contrast case, red curve of Figure 6) or R < Rs (negative contrast case, blue curve of Figure 6), then the pixel is classified as a probable oil spill pixel.In the 12° < α < 17° case, a pixel is classified as a probable spill, whether R > Rs of the red curve or R < Rs of the blue curve, since both positive and negative oil-water contrast can coexist.The pixels in Figures 7 and 8 with 12° < α < 17° and R > Rs (red curve) were colored red, while those with R < Rs (blue curve) were colored blue.Figures 7 and 8 show the R maps together with probable spills detected by our method, highlighted in red (for oil spill pixels classified with positive contrast) and blue (negative contrast), for all of the validation scenes.In the best cases (Figures 7d,e and 8c), the methodology works very well, in that the white spill patch (Figures 7a,b and 8b) is entirely highlighted, with very few and small look-alikes.
Quantitatively, we find that the percentage of the manually-digitized ROI pixels successfully detected for the validation Scenes 8-10 is 100% and 75% for Scene 11.From Figure 7 (Gulf of Mexico), it is also well evident how the contrast in the R image is highly variable within each spill (Figure 7a-c), in that R in the slick ranges from just above R = 1 to R > 2.5 in the center of the slick.This variability justifies the use of the Rs vs. L'GN thresholding curve, rather than an [Rs, maximum R] range for detection.
Figure 7f, on the other hand, seems to reveal some problems in the detection.Indeed, large, obviously non-spill areas are also highlighted, both as positive and negative contrast probable spills, respectively, in red and blue.In particular, blue areas include both small dark features (Figure 7f), probably oil spills, as well as large clean sea regions.Red areas include a wide stripe to the right (east) of the image, with very high R values (Figure 7c).The latter problem occurs due to clouds and the vicinity of the no-glint area, as described above.Finally, one notices also that regions in Figure 7f classified as negative contrast (blue) areas are too large to be spills.This is understandable, because we developed the negative contrast blue line in Figure 6 with two data points only and because the image (Figure 7) is particularly cloudy.Furthermore, for this particular case in proximity of the Mississippi Delta (Figure 7f), both positive and negative contrast non-spill features detected by our method might be due to elevated levels of suspended particulate matter discharged by the river, which alter reflectance.
Figure 8, relative to the Arabian Gulf case treated by [6], summarizes briefly the oil spill detection methodology.Indeed, from the input data, i.e., MODIS TOA radiance Lt (859 nm) (Figure 8a), the ratio image R is computed (Figure 8b).By applying the threshold curve Rs(L'GN) (red curve in Figure 8d), we obtain the map of probable spills (Figure 8c).Finally, by visually comparing Lt and R (Figure 8a,b), it is also evident how the ratio map R is more uniform and "flattened" than Lt.

Conclusions
We presented a new methodology for the detection of marine oil spills using MODIS glint-contaminated, high-resolution (250 m), near-infrared (λ = 859 nm) images.The methodology was developed using a set of seven glint-contaminated MODIS images (training dataset) and validated using four other images (validation dataset), in all of which, certified oil slicks are present.The dataset is comprehensive for positive and negative oil-water contrast cases.
Oil spill detection is based on the ratio R = L'GN/LGN, where L'GN is the MODIS-retrieved normalized sun glint radiance and LGN the same product, but obtained from the Cox and Munk isotropic (wind-direction-independent) sun glint model.Indeed, we showed that in the R ratio images, while clean water pixel values tend to one, oil spills stand out as anomalies better than in other products, such as top of atmosphere radiances or L'GN and LGN examined separately.The detection of a probable oil spill is then automatically achieved by comparing the values in an R image to threshold values Rs obtained from the known slicks of the training dataset.Another novel feature introduced, to our knowledge, for the first time, consists of the possibility of detecting slicks with positive and/or negative contrast via a criterion that assigns the expected contrast type to a pixel, given its illumination and view conditions.Two different Rs(L'GN) threshold fitting curves are then used for detection in the two contrast cases.Results indicate that the above thresholding is able to isolate the spills and that the spills of the validation dataset are successfully detected.Finally, the detected features include some obvious non-spill items, which may be pruned manually or with a clustering algorithm.
We feel that the method presented is therefore a novel and promising, though preliminary, approach towards automatic oil spill detection in optical satellite images, which advances the knowledge in oil spill detection.
The method may be improved via future research, which should include: (1) the non-spill feature elimination cited above; (2) the enrichment of the oil spill dataset for further training and validation; (3) the use of the Cox and Munk anisotropic (wind-direction-dependent) sun glint model; (4) the implementation of the multiple scattering contribution into the Cox and Munk model [32]; and (5) the use of hyperspectral sensors, such as VIISR, in order to overcome the limitations of having only two high-resolution MODIS bands.Finally, the proposed methodology is satellite-independent and is thus exportable to other sensors, such as the next generation of optical satellites, which includes SENTINEL-3.

Figure 2 .
Figure 2. (a) MODIS-derived L'GN (bias-corrected), relative to Scene 2 of Table 1 (Lebanon coast, positive oil spill-water contrast); (b) LGN from Cox and Munk; (c) R = L'GN/LGN; (d) comparison between LGN (solid line: Cox and Munk (1954) isotropic) and L'GN (dotted line) along the transect (black line in (a)).Green dots denote the oil spill pixels in the corresponding ROI; (e) LGN and L'GN histograms; (f) R values along the transect (black line in (c)) with green dots for oil spill pixels.The solid red line represents the oil spill separation index value Rs from Table3.

Figure 3 .
Figure 3. (a) MODIS-derived L'GN (bias-corrected), relative to Scene 6 of Table 1 (Crete, negative oil spill-water contrast); (b) LGN from Cox and Munk; (c) R = L'GN/LGN; (d) comparison between LGN (solid line: Cox and Munk (1954) isotropic) and L'GN (dotted line) along the transect (black line in (a)).Green dots denote the oil spill pixels in the corresponding ROI; (e) LGN and L'GN histograms; (f) R values along the transect (black line in (c)) with green dots for oil spill pixels.The solid red line represents the oil spill separation index value Rs from Table3.

Figure 4 .
Figure 4. (a) MODIS-derived L'GN (bias-corrected), relative to Scene 7 of Table 1 (Gulf of Mexico, positive and negative oil spill-water contrast); (b) R = L'GN/LGN; (c,d) comparison between LGN (solid line: Cox and Munk (1954) isotropic) and L'GN (dotted line) along the transects crossing the spills with negative (red line in (a)) and positive contrast (blue line in (a)).Green dots denote the oil spill pixels in the corresponding ROI; (e,f) R values along the transects (red and blue line in (b)) with green dots for oil spill pixels.Solid red lines represent the oil spill separation index value Rs from Table3.

Figure 5 .
Figure 5. R for the MODIS image of: (a) Lebanon (Scene 2, positive contrast); (b) Crete (Scene 6, negative contrast); (c) Gulf of Mexico (Scene 7, positive and negative contrast).In red, pixels with R ≥ Rs.In blue, pixels with R ≤ Rs.The Rs values used for detection are also indicated.

Figure 6 .
Figure 6.Oil spill index Rs as a function of L'GN for the seven scenes of the training dataset (Tables1 and 3).The red and blue curves are the cubic spline and linear interpolation, respectively, between Rs values.The red (blue) curve represents the lower (upper) R limit for oil spill detection in positive (negative) contrast situations.

Figure 7 .
Figure 7. Oil spill R maps relative to the three validation scenes of Table 1 (Deepwater Horizon oil spill disaster of 2010): (a) 25 April 2010 (Scene 8); (b) 29 April 2010 (Scene 9); (c) 25 May 2010 (Scene 10); (d-f) the regions detected as probable oil spills, same scenes; (g,h) black dots are R values along transects indicated by black lines in (a,b), while red lines are Rs values along the same transect, as inferred from the curves of Figure 6; (i) comparison between LGN (solid line: Cox and Munk (1954) isotropic) and L'GN (dotted line) along the transect in (c).Glint-free (i.e., L'GN < 0.005) areas are marked in grey in (c,f).

Figure 8 .
Figure 8. Arabian Gulf oil spill event (Scene 11).The methodology steps in brief: (a) input data, TOA radiance Lt (859 nm); (b) ratio map R; (c) oil spills detected in R; (d) black dots are R values along the black transect in (c), and the red line is Rs values along the same transect, as inferred from the curves of Figure 6.

Table 1 .
Oil spill database.Scenes 1-7 were used to develop the oil spill detection methodology (training dataset), while 8-11 were used for validation.For each scene, the MODIS file relative to the oil spill observation, acquisition date, bounding box surrounding the spill and oil-water contrast are given; e.g., A20062001040, is the standard NASA name of the MODIS pass and refers to the image from the AQUA sensor on 19 July 2006, 10:40 UTC.

Table 2 .
Retrieved dominant aerosol models, aerosol model ratio and τa (λ = 859 nm) in regions where glitter is negligible for the 11 scenes of Table 1.Aerosol model characteristics are indicated in Columns 3 and 4 (e.g., 26 r70f05v01 0.35 indicates model No. 26 in the current version of SeaDAS, with 70% air column humidity, 5% fine aerosol fraction and the Angstrom coefficient alpha = 0.35).