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

Abstract

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/L_GN, where L'_GN is the MODIS-retrieved normalized sun glint radiance image and L_GN 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 R_s 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 R_s, 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 satellite images.

1. 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 L_GN) are computed from the two procedures listed above, by taking their ratio R = L'_GN/L_GN, 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 R_s 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 R_s(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.

2. The Data

2.1. 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).

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.

Scene	MODIS File	Date	Oil Spill Location	Oil Spill-Water Contrast
1	A20062001040	19 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
2	T20062040835	23 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
3	A20062071045	26 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
4	A20062161040	4 August 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
5	T20101441650	24 May 2010	Gulf of Mexico	Positive
			90.00°–87.10°W, 27.50°–30.50°N
6	A20042441120	31 August 2004	Crete	Negative
			25.00°–27.50°E, 32.00°–34.00°N
7	A20101471855	27 May 2010	Gulf of Mexico	Positive/Negative
			90.00°–87.10°W, 27.50°–30.50°N
8	A20101151855	25 April 2010	Gulf of Mexico	Positive
			88.90°–87.80°W, 28.46°–29.50°N
9	T20101191655	29 April 2010	Gulf of Mexico	Positive
			88.90°–87.80°W, 28.46°–29.50°N
10	A20101451905	25 May 2010	Gulf of Mexico	Positive
			90.00°–87.10°W, 27.50°–30.50°N
11	A20132290935	17 August 2013	Arabian Gulf	Positive
			52.70°–53.50°E, 25.90°–26.30°N

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.

Scene	MODIS File	Date	Oil Spill Location	Oil Spill-Water Contrast
1	A20062001040	19 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
2	T20062040835	23 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
3	A20062071045	26 July 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
4	A20062161040	4 August 2006	Lebanon coast	Positive
			34.50°–36.00°E, 33.50°–35.00°N
5	T20101441650	24 May 2010	Gulf of Mexico	Positive
			90.00°–87.10°W, 27.50°–30.50°N
6	A20042441120	31 August 2004	Crete	Negative
			25.00°–27.50°E, 32.00°–34.00°N
7	A20101471855	27 May 2010	Gulf of Mexico	Positive/Negative
			90.00°–87.10°W, 27.50°–30.50°N
8	A20101151855	25 April 2010	Gulf of Mexico	Positive
			88.90°–87.80°W, 28.46°–29.50°N
9	T20101191655	29 April 2010	Gulf of Mexico	Positive
			88.90°–87.80°W, 28.46°–29.50°N
10	A20101451905	25 May 2010	Gulf of Mexico	Positive
			90.00°–87.10°W, 27.50°–30.50°N
11	A20132290935	17 August 2013	Arabian Gulf	Positive
			52.70°–53.50°E, 25.90°–26.30°N

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.

Figure 1. (a) MODIS TERRA TOA radiance image, Lt(859 nm). 23 July 2006, 08:35 UTC, Lebanon coast (Scene 2 of Table 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.

Figure 1. (a) MODIS TERRA TOA radiance image, Lt(859 nm). 23 July 2006, 08:35 UTC, Lebanon coast (Scene 2 of Table 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.

2.2. 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 therefore obtained the following MODIS L2 parameters (using SeaDAS nomenclature), all at λ = 859 nm (250-m spatial resolution):

Observed radiance products:

• Top of atmosphere (TOA) radiance: L_t;

Atmospheric and other intermediate products:

• Rayleigh radiance: L_r;
• Aerosol radiance: L_a;
• Aerosol optical thickness: taua (τ_a used hereafter);
• Aerosol models: aer_model_min, aer_model_max, aer_model_ratio;
• Normalized sun glint radiance: glint_coef (L_GN hereafter).

Solar-viewing geometry products:

• Sun and solar zenith angles: solz, senz (θ₀, θ hereafter);
• Sun and solar azimuth angles: sola, sena.

Ancillary products:

• Wind speed: wind_speed.

We reiterate that TOA radiance L_t is the radiance received by the sensor containing both the atmospheric and surface signal; Rayleigh radiance L_r and aerosol radiance L_a 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 L_GN 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:

cosα = cosθ₀cosθ − sinθ₀sinθcosΔϕ

(1)

The glint angle α is the angle between the satellite viewing direction and the direction of the specular reflection, with θ₀, θ 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.

3. Oil Spill Detection Methodology

3.1. Background

The total radiance at the top of the ocean-atmosphere system can be written as (Gordon and Wang [26]):

L_t(λ) = L_r(λ) + L_a(λ) + L_ra(λ) + T(λ)L_g(λ) + t(λ)L_wc(λ) + t(λ)L_w(λ)

(2)

where L_r(λ) is the radiance resulting from multiple scattering by air molecules (Rayleigh scattering) in the absence of aerosols, L_a(λ) is the radiance resulting from multiple scattering of aerosols in the absence of air and L_ra(λ) is the interaction term between molecular and aerosol scattering [26]. L_g(λ) is the glint radiance due to specular reflection of the direct solar beam, while L_wc(λ) and L_w(λ) 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, L_w(λ) 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., L_wc(λ), 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 L_t 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 L_ra into a unique aerosol contribution L_A = L_a + L_ra [28], we can write L_t as:

L_t = L_r + L_A + TL_g

(3)

The sun glint radiance L_g can be therefore expressed as:

L_g = F₀T₀L_GN

(4)

where F₀ and T₀ are, respectively, the extraterrestrial solar irradiance and the atmospheric direct transmittance along the solar direction. L_GN 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 F₀ = 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:

L_GN = r(ω)P(β, σ)/(4cosθ₀cosθcos⁴β)

(5)

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⁻¹). In practice, L_GN depends on solar zenith θ₀, sensor zenith θ, the relative Sun-satellite azimuth angle Δϕ and surface wind velocity W, because also β depends on the above angles [28], i.e., L_GN = L_GN (θ₀, θ, Δϕ, W). The solar-sensor geometry (θ₀, θ, Δϕ) 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 L_GN ≥ 0.005 are flagged, and no further processing is done on MODIS data above this L_GN threshold. In areas where 0 < L_GN < 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 < L_GN < 0.01 and, thus, to treat all pixels with L_GN > 0.005.

3.2. Oil Spill Detection Strategy

By using Equations (3) and (4), we can express the normalized sun glint radiance as:

L'_GN = (L_t − L_r − L_A)/F₀T₀T

(6)

with:

T₀T = exp{−[ τ_r + τ_a] (1/cosθ₀ + 1/cosθ)}

(7)

where τ_r and τ_a are the Rayleigh and aerosol optical thicknesses and θ₀ 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 L_GN. 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:

R = L'_GN/L_GN

(8)

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 L_GN somewhat cancel each other in the ratio. Furthermore, R is independent from atmospheric aerosol and Rayleigh scattering, since L'_GN and L_GN 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 L_GN.

3.3. MODIS Sun Glint Radiance Computation (L'_GN)

To obtain L'_GN (Equation (6)), we need the aerosol radiance L_A 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 L_A over the glint area using such an aerosol model and this <τ_a>.

To obtain L_A 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 (L_GN < 10⁻⁷). 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.

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).

Scene	Location	Aerosol Model Min	Aerosol Model Max	Aerosol Model Ratio	τa(859)
1	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.49	0.18
		r70f05v01	r70f10v01
		0.35	0.68
2	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.55	0.12
		r70f05v01	r70f10v01
		0.35	0.68
3	Lebanon coast 29.00°–34.00°E, 33.50°–35.00°N	26	25	0.38	0.21
		r70f05v01	r70f10v01
		0.35	0.68
4	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.70	0.18
		r70f05v01	r70f10v01
		0.35	0.68
5	Gulf of Mexico 90.00°–84.50°W, 27.50°–30.50°N	57	56	0.78	0.30
		r85f02v01	r85f05v01
		0.15	0.45
6	Crete 24.00°–29.00°E, 32.00°N–34.00°N	24	23	0.51	0.08
		r70f20v01	r70f20v01
		1.12	1.40
7	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	26	25	0.50	0.18
		r70f05v01	r70f10v01
		0.35	0.68
8	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	26	25	0.60	0.10
		r70f05v01	r70f10v01
		0.35	0.68
9	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	27	26	0.50	0.15
		r70f02v01	r70f05v01
		0.07	0.35
10	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	25	24	0.46	0.15
		r70f10v01	r70f20v01
		0.68	1.12
11	Arabian Gulf 49.00°–56.00°E, 24.00°–27.00°N	6	5	0.80	0.35
		r30f05v01	r30f10v01
		0.29	0.61

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).

Scene	Location	Aerosol Model Min	Aerosol Model Max	Aerosol Model Ratio	τa(859)
1	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.49	0.18
		r70f05v01	r70f10v01
		0.35	0.68
2	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.55	0.12
		r70f05v01	r70f10v01
		0.35	0.68
3	Lebanon coast 29.00°–34.00°E, 33.50°–35.00°N	26	25	0.38	0.21
		r70f05v01	r70f10v01
		0.35	0.68
4	Lebanon coast 30.00°–34.00°E, 33.50°–35.00°N	26	25	0.70	0.18
		r70f05v01	r70f10v01
		0.35	0.68
5	Gulf of Mexico 90.00°–84.50°W, 27.50°–30.50°N	57	56	0.78	0.30
		r85f02v01	r85f05v01
		0.15	0.45
6	Crete 24.00°–29.00°E, 32.00°N–34.00°N	24	23	0.51	0.08
		r70f20v01	r70f20v01
		1.12	1.40
7	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	26	25	0.50	0.18
		r70f05v01	r70f10v01
		0.35	0.68
8	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	26	25	0.60	0.10
		r70f05v01	r70f10v01
		0.35	0.68
9	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	27	26	0.50	0.15
		r70f02v01	r70f05v01
		0.07	0.35
10	Gulf of Mexico 90.00°–85.50°W, 27.50°–30.50°N	25	24	0.46	0.15
		r70f10v01	r70f20v01
		0.68	1.12
11	Arabian Gulf 49.00°–56.00°E, 24.00°–27.00°N	6	5	0.80	0.35
		r30f05v01	r30f10v01
		0.29	0.61

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 L_A 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, T₀ 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 L_GN, we computed the mean bias (MBE) and standard deviation of the L'_GN − L_GN 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 L_GN. The MBEs are found to be about 10⁻³–10⁻², i.e., 10%–15% of the typical L_GN 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/L_GN 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/L_GN (Equation (8)). Hereafter, we indicate with the same symbol L'_GN the bias-corrected glint radiances, i.e., L'_GN = L'_GN (uncorrected) − MBE.

4. Results

4.1. Evaluation of L_GN, L'_GN and R

We computed L_GN, L'_GN and the ratio images R = L'_GN/L_GN for all of the scenes of the dataset (Table 1). Figure 2, Figure 3 and Figure 4 show the results of L_GN, 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 L_GN 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 (Figure 2 and Figure 3), L'_GN and L_GN tend to coincide in clean sea water pixels (black line and dots in Figure 2d and Figure 3d), giving values of R close to one (Figure 2f and Figure 3f), while they differ (R ≠ 1) in oil spill pixels (Figure 2d and Figure 3d, green dots) and a stronger separation occurs for cloudy pixels (high-value isolated black dots in the same figures). The overall behavior of L_GN and L'_GN is shown by their histograms (Figure 2e and Figure 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.

Figure 2. (a) MODIS-derived L'_GN (bias-corrected), relative to Scene 2 of Table 1 (Lebanon coast, positive oil spill-water contrast); (b) L_GN from Cox and Munk; (c) R = L'_GN/L_GN; (d) comparison between L_GN (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) L_GN 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 R_s from Table 3.

Figure 2. (a) MODIS-derived L'_GN (bias-corrected), relative to Scene 2 of Table 1 (Lebanon coast, positive oil spill-water contrast); (b) L_GN from Cox and Munk; (c) R = L'_GN/L_GN; (d) comparison between L_GN (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) L_GN 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 R_s from Table 3.

For a few cases of the dataset and only for some pixels of the images, L'_GN and L_GN 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 L_GN 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 L_GN 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 L_GN 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 L_GN with the Cox and Munk (1954) anisotropic (wind-direction-dependent) model.

Figure 3. (a) MODIS-derived L'_GN (bias-corrected), relative to Scene 6 of Table 1 (Crete, negative oil spill-water contrast); (b) L_GN from Cox and Munk; (c) R = L'_GN/L_GN; (d) comparison between L_GN (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) L_GN 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 R_s from Table 3.

Figure 3. (a) MODIS-derived L'_GN (bias-corrected), relative to Scene 6 of Table 1 (Crete, negative oil spill-water contrast); (b) L_GN from Cox and Munk; (c) R = L'_GN/L_GN; (d) comparison between L_GN (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) L_GN 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 R_s from Table 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 L_t 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 L_t for all of the scenes with positive contrast of the training dataset, we obtain a total L_t range of [2.50, 15.00] mw∙cm⁻²∙µm⁻¹∙sr⁻¹ 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.

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/L_GN; (c,d) comparison between L_GN (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 R_s from Table 3.

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/L_GN; (c,d) comparison between L_GN (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 R_s from Table 3.

4.2. 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, R_s, 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 ≥ R_s in the same R image, we obtain a map of probable spills with positive contrast (Figure 5a, red features). For negative contrast, R_s is the value below which we find approximately 70% of the ROI R values. Thus, highlighting the pixels with R < R_s, 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.

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 ≥ R_s. In blue, pixels with R ≤ R_s. The R_s values used for detection are also indicated.

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 ≥ R_s. In blue, pixels with R ≤ R_s. The R_s values used for detection are also indicated.

What is evident from Figure 5 is that also a subset of non-oil spill pixels are classified as probable oil spills using the R_s 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 R_s, 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 R_s index values found for each scene of the training dataset are listed in Table 3. We found six R_s 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 ≥ R_s (R ≤ R_s) 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.

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

Figure 6. Oil spill index R_s as a function of L'_GN for the seven scenes of the training dataset (Table 1 and Table 3). The red and blue curves are the cubic spline and linear interpolation, respectively, between R_s 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 R_s 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 R_s 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 R_s 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 Figure 7 and Figure 8 below) those pixels with R values above (positive contrast, in red) or below (negative contrast, in blue) our R_s 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 L_GN. 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 L_GN > 0.05, while negative contrast is to be expected for α > 17°. All other regions (those with 12° < α < 17° or L_GN < 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. R_s 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 ≥ R_s (R ≤ R_s) 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.

Scene	Rs (L'GN/LGN) (Positive Contrast)	Rs (L'GN/LGN) (Negative Contrast)	LGN	L'GN	Glint Angle α	Solar Zenith Angle θ 0	Sensor Zenith Angle θ	Sun-Sensor Azimuth Difference Δϕ	Wind Speed
					(°)	(°)	(°)	(°)	(m/s)
1	1.15	-----	0.074	0.100	10.00	18.17	18.32	147.71	2.18
2	1.10	-----	0.060	0.070	10.55	20.14	20.70	149.51	3.23
3	1.02	-----	0.035	0.035	15.32	20.17	6.20	147.30	4.62
4	1.05	-----	0.044	0.045	13.18	21.36	16.95	141.50	3.56
5	1.20	-----	0.10	0.150	5.32	15.91	15.54	164.68	2.79
6	-----	0.80	0.012	0.010	24.13	29.41	8.61	134.52	5.72
7	1.12	-----	0.067	0.075	9.07	17.42	25.32	167.75	3.97
7 bis	-----	0.75	0.025	0.018	17.07	16.25	32.97	171.35	2.32

Table 3. Set of parameters used to characterize oil spills in glint-contaminated MODIS imagery. R_s 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 ≥ R_s (R ≤ R_s) 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.

Scene	Rs (L'GN/LGN) (Positive Contrast)	Rs (L'GN/LGN) (Negative Contrast)	LGN	L'GN	Glint Angle α	Solar Zenith Angle θ 0	Sensor Zenith Angle θ	Sun-Sensor Azimuth Difference Δϕ	Wind Speed
					(°)	(°)	(°)	(°)	(m/s)
1	1.15	-----	0.074	0.100	10.00	18.17	18.32	147.71	2.18
2	1.10	-----	0.060	0.070	10.55	20.14	20.70	149.51	3.23
3	1.02	-----	0.035	0.035	15.32	20.17	6.20	147.30	4.62
4	1.05	-----	0.044	0.045	13.18	21.36	16.95	141.50	3.56
5	1.20	-----	0.10	0.150	5.32	15.91	15.54	164.68	2.79
6	-----	0.80	0.012	0.010	24.13	29.41	8.61	134.52	5.72
7	1.12	-----	0.067	0.075	9.07	17.42	25.32	167.75	3.97
7 bis	-----	0.75	0.025	0.018	17.07	16.25	32.97	171.35	2.32

The glint angle (or glint radiance L_GN) 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.

4.3. 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 (Figure 7 and Figure 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 R_s 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 R_s values, relative to both curves, are considered;

(4)

compares R with the R_s(L'_GN) value for the proper contrast case; two comparisons are made for 12° < α < 17°;

(5)

if R > R_s (positive contrast case, red curve of Figure 6) or R < R_s (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 > R_s of the red curve or R < R_s of the blue curve, since both positive and negative oil-water contrast can coexist. The pixels in Figure 7 and Figure 8 with 12° < α < 17° and R > R_s (red curve) were colored red, while those with R < R_s (blue curve) were colored blue.

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 R_s values along the same transect, as inferred from the curves of Figure 6; (i) comparison between L_GN (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 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 R_s values along the same transect, as inferred from the curves of Figure 6; (i) comparison between L_GN (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. Arabian Gulf oil spill event (Scene 11). The methodology steps in brief: (a) input data, TOA radiance L_t (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 R_s values along the same transect, as inferred from the curves of Figure 6.

Figure 8. Arabian Gulf oil spill event (Scene 11). The methodology steps in brief: (a) input data, TOA radiance L_t (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 R_s values along the same transect, as inferred from the curves of Figure 6.

Figure 7 and Figure 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 (Figure 7d,e and Figure 8c), the methodology works very well, in that the white spill patch (Figure 7a,b and Figure 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 R_s vs. L'_GN thresholding curve, rather than an [R_s, 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 L_t (859 nm) (Figure 8a), the ratio image R is computed (Figure 8b). By applying the threshold curve R_s(L'_GN) (red curve in Figure 8d), we obtain the map of probable spills (Figure 8c). Finally, by visually comparing L_t and R (Figure 8a,b), it is also evident how the ratio map R is more uniform and “flattened” than L_t.

5. 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/L_GN, where L'_GN is the MODIS-retrieved normalized sun glint radiance and L_GN 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 L_GN examined separately. The detection of a probable oil spill is then automatically achieved by comparing the values in an R image to threshold values R_s 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 R_s(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.