1.

Introduction

The Deepwater Horizon (DWH) oil spill in the Gulf of Mexico, which flowed unabated for three months in 2010, was the largest accidental marine oil spill in the history of the petroleum industry; its source was a sea-floor oil gusher resulting from the April 20, 2010, DWH explosion which claimed 11 lives. The gushing wellhead was capped only after 87 days, on July 15, 2010. Beginning with the first days after the accident, all the Earth observing satellites focused their image acquisitions over the Gulf of Mexico. Among the many sensors on board the satellites, the synthetic aperture radar (SAR) is certainly the most powerful one for imaging different ocean phenomena, like waves, surface winds, oil spills, and sea-ice in all-weather conditions. Thanks to a European Space Agency (ESA) project for studying ocean phenomena, among the various SAR products covering the DWH accident, the Envisat/ASAR wide swath (WS) ones were chosen; in particular, two significative ASAR WS images were selected: the first available after the accident, on April 26, and the second one when the oil spill was already fully developed, on May 2. The interaction of the highly coherent radiation of a radar signal with the ocean elements combined with the atmospheric conditions produces a very complex backscatter.^1,2 The geophysical system collaterally creates the speckle phenomenon, which produces the characteristic grainy appearance of SAR images. While speckle can even be exploited to analyze SAR oil spills at full resolution,³ it generally causes difficulties in image interpretation and is usually removed with specialized filters^4,5 but at the risk of degrading the spatial resolution. Heavily oil-polluted ocean surfaces provide a specular reflection and a reduced Bragg scattering. Compared with the semispecular backscatter of the open sea surfaces, pixels of oil spills produce darker signatures in SAR images. Several other natural phenomena, such as low tides, low wind areas, biogenic material, and oceanic or atmospheric fronts, produce the so-called look-alikes of oil slicks.^6,7 Various papers have been presented in the literature describing semiautomatic parametric and nonparametric algorithms for oil spill detection, for instance, adaptive thresholding segmentation methods,^8,9 neural networks,¹⁰ fractal algorithms,¹¹ and multiscale wavelet representations.¹² A different group of algorithms based on geometric and statistical features have also been used,¹³ while those based on the physical modeling of complex systems require oil viscoelastic properties and external data such as scatterometer wind fields.¹⁴ Finally, a recent review paper¹⁵ provides a comprehensive analysis of all the issues related to oil spill detection with SAR images. A relevant paper¹⁶ has demonstrated that only dual- and/or full-polarimetric SAR images can give precise oil spill detection, clearly distinguishing between oil and look-alikes; in the case of single-polarimetric images, the same objective can only be attained if external information is attached to the image.¹⁷ Unfortunately, regular acquisitions by SAR observing satellites are mainly single-polarimetric, while dual- and/or full-polarimetric ASAR images are extremely rare in the ESA archive. This paper presents a processing scheme based on binary segmentation whose aim is to provide an efficient tool to measure the marine oil spill extent in SAR images; for the reasons explained above, it was decided to apply the processing scheme to single-polarimetric images, with an approach that only makes use of the radiometric information of the SAR scene. The segmentation process is modeled by taking into account the Bayes formulation. The optimization of the probability functions by means of the Markov random field (MRF) theory is defined by the maximum a posteriori (MAP) criterion. The stochastic minimization algorithm is modified in order to update both the parameters of the a priori label model and its system neighborhood.

2.

Methods

2.1.

Bayesian Framework

Many events in nature show a behavior that is difficult to predict as they are affected by different variables. During the process of acquisition of an SAR image, a given pixel may take an unpredictable value, even when it belongs to objects with a well-known nature. This random behavior justifies the application of probability theory in SAR image analysis. Making use of probability theory and statistics, in this section, different observations are integrated in order to define specific parametric models. Let $X$ be the image to be segmented and $W$ the segmentation result, where $w_i\in \{\mathrm{1,2},\dots ,m\}$ and $m$ is the number of classes (or labels). Bayes’ rule,

Eq. (1)

$$p(w_j|X)=p(X|w_j)p(w_j)/\sum _{i=1}^{m}P(X,w_i),$$

allows the computation of the posterior probability $p(w_j|X)$ of class $w_j$. The conditional probability $p(X|w_j)$ must be derived and the a priori probability of event $w_j$ must be known before solving the problem. MRFs have the property that the conditional distribution of a particular pixel given the values of all the other pixels in the whole image, is equal to the conditional distribution obtained by only considering the pixel values in its neighborhood.

In order to incorporate the pixel information to the terms of probability of Bayes’ rule, the use of some topological definitions is required. The sites in a grid $S$ are related to one another via a neighborhood system. A neighborhood system for $S$ is defined as $N=\{N_i|\forall i\in S\}$, where $N$ is the set of sites of neighboring pixel $i$. A random region $X$ on a lattice $S$ with neighborhood system $N$ is said to be an MRF if $\forall i\in S$, $p(x_i|x_r\forall r\ne i)=p(x_i|x_{Nr})$. A clique $C$ for $(S,N)$ is defined as a subset of sites in $S$. It consists either of a single site $C_1=\{i\}$, or a pair of neighboring sites $C_2$, or a triplet of neighboring sites $C_3$, and so on. The collections of single-site, pair-site, and triple-site cliques will be denoted by $C_1$, $C_2$ and $C_3$, respectively. The collection of all cliques for $(S,N)$ is $C=C_1\cup C_2\cup C_3\dots$. The Bayes’ rule can be applied in order to define the equations of the random field model. The joint probability of radiometric information is $p(X,W)=p(X|W)p(W)$ and is derived from the observed data. The image acquisition process is described by means of the equation $p(X|W)=\prod _{s\in S}p(x_s|w_s)$. The MAP approach uses the posterior probability in order to obtain an optimization of the segmentation process.

Eq. (2)

$$\begin{gathered} \underset{\mathrm{MAP}}{\overset{\circ }{W}}=p(W|X)=\underset{W}{\arg \max }\frac{p(X|w_j)p(w_j)}{\sum _{i=1}^{m}P(X,w_i)} \\ \approx \underset{W}{\arg \max }[\log p(X|W)+\log p(W)]. \end{gathered}$$

The law of the total probability $\sum _{i=1}^{m}P(X,w_i)$ can be ignored and the MAP equation can be approached using the numerator term. In a pixel-based segmentation approach, every $x_s$ pixel is assigned to the $w_j$ class that maximizes Eq. (2). Given the local heterogeneity of SAR images, a good strategy is to take into account the information of the set of sites neighboring pixel $x_s$. In a contextual-based approach, a discrete MRF is applied for modeling the segmentation problem. The a posteriori energy function is derived from Eq. (2) and results in

Eq. (3)

$$\underset{\mathrm{MAP}}{\overset{\circ }{W}}=\underset{W}{\arg \min }U(W|X)\approx \underset{W}{\arg \min }[U(X|W)+U(W)],$$

where $U(W)$ is the a priori energy function and $U(X|W)$ is the energy function of the observed data model. The Bayesian formulation requires that assumptions are made about prior probabilities $P(W)$. A simple Potts model is used to define the a priori energy function $U(W)$.¹⁸ The potential function is given by

Eq. (4)

$$V(w_i,w_s)=\{\begin{array}{ll}\beta & \delta (w_i\ne w_s)\\ -\beta & \delta (w_i=w_s)\end{array},$$

where $\beta$ is a model parameter which indicates the convergence of the solution, and $\delta$ is the Dirac’s delta function which has a unitary value when the classes of pixels $x_i$ and $x_s$ are equal. The joint image model integrates both the label field to be obtained and the acquisition process of the image. Data can be modeled by a Gaussian function $f_x(x)=N({\mu }_i,{\sigma }_i^2)$.

Eq. (5)

$$f_x(x)=1/\sqrt{2\pi }\sigma \exp [-{(x-{\mu }_i)}^2/{2\sigma }_i^2],$$

with mean values ${\mu }_i$ and variance ${\sigma }_i^2$. The Hammersley-Clifford theorem defines the probability density for an MRF under the form of a Gibbs distribution,¹⁹ but the ill-posed problem is to calculate Eq. (3). The solution by a stochastic approach needs a minimization method, for example, the simulated annealing²⁰ given by the Metropolis sampler.²¹ The probability to belong to a particular class is calculated by means of the posterior energy function $U(W|X)$, and considering Eqs. (4) and (5), the operation is

Eq. (6)

$$U(W|X)=\underset{W}{\arg \min }\{\sum _{s\in S}[\ln (\sqrt{2\pi }\sigma )+{(x-{\mu }_i)}^2/{2\sigma }^2]+\sum _{c\in C}V(c)\},$$

where $C$ is the clique system. Equation (6) exposes an ill-posed problem: the class field is degraded by noise such that the observed SAR images are considered as incomplete data, so Eq. (6) can generate a family of nonoptimal solutions. In a classical solution, a set of a priori information is provided. As a first step, the parameters belonging to the set $\mathrm{\Theta }=\{{\mu }_i,{\sigma }_i^2\}$, $i=\mathrm{1,2},\dots ,m$ are derived from the existing data and are set as constants for the rest of the process. Then, the MAP approach is applied to the image $X$ and the label field $W$ is obtained. The minimization method is summarized as follows:

• 1. Initialization of parameters ($N$, $n$, $W^0$, $T^0$)

• 2. Do $n$ iterations:

  • (a) Sample $X$ in order to obtain $U(X|W)$

  • (b) Sample $W^i$ in order to obtain $U(W)$

  • (c) Compute $U(W|X)$

  • (d) Apply the Metropolis criterion

  • (e) Decrease $T^i$ after each iteration

• 3. Return $W^n$,

where $N$ is the neighboring pixel system, $n$ is the number of iterations, $W^0$ is a proposed initial random solution, and $T^0$ is the initial temperature parameter. Related to thermodynamic systems, the applied segmentation process takes into account a cooling schedule where the image pixels go slowly from high temperature to zero temperature distributions. For this reason, the method is known as simulated annealing and $T$ is just a symbolic representation of temperature. In a recursive minimization, $T^0$ has a high value and the probability distribution of the label field $W^0$ is uniform. Labels $w_i$ are randomly assigned with uniform probability, where $w_i\in \{\mathrm{1,2},\dots ,m\}$ and $m$ is the number of classes (or labels). After a number of iterations, when $n\to \infty$;, $T^n\to 0$, $W^n$ is the segmentation result and Eq. (6) attains a minimum value. The configuration of the $W^i$ label field is generated with the Metropolis criterion and assumes a Gibbs distribution. Figure 1 shows a conceptual block scheme. The segmentation is obtained by iterations: at a given iteration, both the original SAR data for deriving the energy term, $U(X|W)$, and the result of the previous segmentation field, $U(W)$, are required. The update stage is described in the following section and is part of the proposal of this paper.

Fig. 1

Conceptual block scheme of the recursive segmentation process.

3.

Parameter Estimation

The Gaussian model parameters are unknown but the ergodic property of the stochastic process can be assumed so that the mean values ${\mu }_i$ can be estimated from training windows. The variance for every class is approached by

A simpler solution of Eq. (8), obtained by a suboptimal criterion, is preferred.

In SAR images, oil slicks produce a low backscatter, which is expected to be different from the open sea signature. Our proposal is to simultaneously update the $\stackrel{\circ }{\mathrm{\Theta }}$ parameters of the Gaussian model and to perform an iterative adaptive segmentation in order to obtain the $\stackrel{\circ }{W}$ field.

The pseudo code for solving Eqs. (9) and (10) is as follows:

4.

Materials

As stated in the Introduction, among the many ASAR WS images of the DWH oil spill, two of them were selected for testing our processing scheme: the first image available after the accident (April 26) and the second one when the oil spill was already fully developed (May 2). Figures 2(a) and 2(b) show the two images after being ingested into the commercial software TeraScan of Seaspace Corp® and projected onto an equidistant cylindrical geographical reference. The geo-referenced images are shown with the lat/long grid labels to provide their geo-location. The size of the resampled pixel is $150\times 150\mathrm{m}=\mathrm{22,500}{\mathrm{m}}^2$. Table 1 gives the technical data of the ASAR WS images, which are the standard products systematically generated by the ESA Payload Data Handling Station using the ScanSAR technique.

Fig. 2

Geo-referenced ASAR wide swath images of the Deepwater Horizon oil spill after projection onto an equidistant cylindrical projection: (a) April 26 and (b) May 2.

Table 1

Wide swath medium-resolution image (ASA WSM 1P) product.

5.

Results

From the two speckled images of Fig. 2, a window covering a large part of the oil slick was extracted. The two windows, after the correction of the slant range effect on the backscatter, became our test images [Figs. 3(a) and 3(b)]; their sizes are $750\times 800\text{pixels}$ and $850\times 700\text{pixels}$, covering an area of $\sim 13,500$ and $\mathrm{13,388}{\mathrm{km}}^2$, respectively.

Fig. 3

Test images as subsets of the two ASAR raw images after correction of the slant range effect on the backscatter: (a) April 26 and (b) May 2.

The processing scheme can now start with the adaptive algorithm parameters fixed to $n=30$ and $n_2=10$. With regard to the algorithm and to its parameters, some elucidations are needed.

Finally, steps 2 and 3 are repeated until vector $\stackrel{\circ }{\mathrm{\Theta }}$ convergence, which happens in a maximum of 50 iterations. Partial results, shown in Figs. 4 and 5, allow the comparison of the performances of the two algorithms: the conventional MRF (CMRF) segmentation and the new proposed MRF (NMRF) scheme. CMRF and NMRF algorithms were run with the same set of initial parameters on subsets of the original test images, Figs. 4(a) and 5(a). Using a constant second-order neighborhood system, results of CMRF segmentation are shown in Figs. 4(b) and 5(b), where a satisfactory binary discrimination is obtained but detection of fine structures is absent. Results of the NMRF are shown in Figs. 4(c) and 5(c), where oil slick particle pathways are preserved and an improved detection of fine structures and details is clearly observed. The recurring sampling simulation was carried out with $\beta =0.35$, an initial temperature parameter $T^0=2$, and second-order cliques. The interaction coefficient $\beta$ determines the contribution of each element of the neighborhood to maintain the convergence of the solution. The Boltzmann annealing regards the relaxation process, to which, in order to allow a faster convergence, the following schedule for decreasing $T$ is applied: $T^n=q{T}^{n-1}$, where $q$ takes the value of 0.95. The numerical simulations were conducted using a PC with Pentium Core(TM) i7-2670QM CPU of 2.2 GHz and memory of 4 GB. The language used to implement the algorithms was MATLAB®. In each simulation, 51 iterations were applied. The computation time of the proposed algorithm was 860.482 and 849.994 s, for the images of April 26 and May 2, respectively.

Fig. 4

Performance of the Markov random fields (MRF) algorithms using an image subset of $150\times 180\text{pixels}$: (a) subset of the original synthetic aperture radar (SAR) window of April 26, (b) result of a conventional MRF segmentation, and (c) result of the proposed new MRF scheme.

Fig. 5

Same as Fig. 4 but with the original SAR window of May 2.

Figures 6 and 7 show the segmentation result obtained with NMRF. In Figs. 6(a) and 7(a), the detected contours are overlapped to the original test images: contours are closed forms and the derived label field is integrated into homogeneous regions, while both regularities of the oil spill distribution and the contextual constraints of the open sea regions are adequately modeled. Figures 6(b) and 7(b) show the binary map of the oil spill detection, while Table 2 gives the figures of the oil slick pixels and area in the two test images.

Fig. 6

Segmentation results with the image of Fig. 2(a): (a) with the contour between the two regions of interest, oil spill and sea, overlapped to the original image and (b) binary detection with data of Fig. 2(a).

Fig. 7

Same as Fig. 6 but with the image of Fig. 2(b).

Table 2

Oil slicks area of the two test images.

It is confirmed that the contextual functionality of the NMRF model provides a great flexibility in facilitating the reduction of the effects of speckle random behavior.

6.

Conclusions

An MRF scheme to model the contextual labeling of oil-polluted ocean surfaces is presented. In a conventional MRF approach, the set of parameters of the related model is fixed before the computation of the minimization solution and they remain unchanged along the recurring cycles. A model with constant parameters is only partially suited to face the statistical characteristics of SAR data. In the proposed scheme, the conditional distributions are estimated from training sets and their parameters are updated during the stochastic relaxation of the posterior energy function; only the pixel values are considered for deriving the probabilistic dependencies of the Bayesian approach. By incorporating an updated neighborhood system, the NMRF model offers better performances than the conventional one in detecting fine structures and provides very satisfactory results in tackling the stochastic features of the sea phenomena.