## 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,^{3} 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,^{10} fractal algorithms,^{11} and multiscale wavelet representations.^{12} A different group of algorithms based on geometric and statistical features have also been used,^{13} while those based on the physical modeling of complex systems require oil viscoelastic properties and external data such as scatterometer wind fields.^{14} Finally, a recent review paper^{15} provides a comprehensive analysis of all the issues related to oil spill detection with SAR images. A relevant paper^{16} 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.^{17} 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,

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 \text{\hspace{0.17em}}\text{\hspace{0.17em}}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 \text{\hspace{0.17em}}\text{\hspace{0.17em}}i\in S$, $p({x}_{i}|{x}_{r}\forall \text{\hspace{0.17em}}\text{\hspace{0.17em}}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.

## (2)

$$\underset{\mathrm{MAP}}{\overset{\circ}{W}}=p(W|X)=\underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}\frac{p(X|{w}_{j})p({w}_{j})}{\sum _{i=1}^{m}P(X,{w}_{i})}\phantom{\rule{0ex}{0ex}}\approx \underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}[\mathrm{log}\text{\hspace{0.17em}}p(X|W)+\mathrm{log}\text{\hspace{0.17em}}p(W)].$$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

## (3)

$$\underset{\mathrm{MAP}}{\overset{\circ}{W}}=\underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}\text{\hspace{0.17em}}U(W|X)\approx \underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}[U(X|W)+U(W)],$$*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)$.

^{18}The potential function is given by

## (4)

$$V({w}_{i},{w}_{s})=\{\begin{array}{ll}\beta & \delta ({w}_{i}\ne {w}_{s})\\ -\beta & \delta ({w}_{i}={w}_{s})\end{array},$$## (5)

$${f}_{x}(x)=1/\sqrt{2\pi}\sigma \text{\hspace{0.17em}}\mathrm{exp}[-{(x-{\mu}_{i})}^{2}/{2\sigma}_{i}^{2}],$$^{19}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

^{20}given by the Metropolis sampler.

^{21}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

## (6)

$$U(W|X)=\underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}\{\sum _{s\in S}[\mathrm{ln}(\sqrt{2\pi}\sigma )+{(x-{\mu}_{i})}^{2}/{2\sigma}^{2}]+\sum _{c\in C}V(c)\},$$*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:

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.

## 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

where ${\mu}_{j}>{\mu}_{i}$. Inflection points of adjacent probability density functions are crossed in the point $({\mu}_{j}-{\mu}_{i})/2$. This is a kind of data-driven term designed to be evasive, letting the model take control of the probability decision. Simulated annealing can be applied in order to obtain binary segmentations, but the quality of the label field heavily depends on the parameters provided to the energy model. Ideally, the training windows must be assigned to prototype regions of the label field which is going to be determined. The statistical disadvantage is that the training windows use data whose joint distribution functions are the sum of several distribution functions. Pixel values are random variables and a training window gives a partial statistical representation of a given class. In the conventional parameter estimation approach, $X$ denotes the noisy observed image and $W$ is the free-noise version of $X$ (or the label field obtained by a segmentation process) with a parameter vector $\mathrm{\Theta}$. The problem is to find the labeling field $W$, and in an MAP approach, the biased solution can be described by## (8)

$$(\stackrel{\circ}{W},\stackrel{\circ}{\mathrm{\Theta}})=\underset{W,\mathrm{\Theta}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}\text{\hspace{0.17em}}p(W,X|\mathrm{\Theta}),$$^{22}Alternative solutions are the algorithms of expectation maximization and the stochastic expectation maximization.

^{23}

^{,}

^{24}

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

## (9)

$$\stackrel{\circ}{W}=\underset{W}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}\text{\hspace{0.17em}}p(W|\stackrel{\circ}{\mathrm{\Theta}}),$$## (10)

$$\stackrel{\circ}{\mathrm{\Theta}}=\underset{\mathrm{\Theta}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}\text{\hspace{0.17em}}p(\stackrel{\circ}{W}|\mathrm{\Theta}).$$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:

1. Initialization of parameters

2. Do $n$ iterations: Apply minimization by stochastic relaxation [Eqs. (6) and (9)].

3. Update vector $\stackrel{\circ}{\mathrm{\Theta}}$ from the previous $n$ label field [Eq. (10)].

4. Update neighborhood

5. Repeat steps 2 and 3 until convergence of $\stackrel{\circ}{\mathrm{\Theta}}$.

## 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}=\mathrm{22,500}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.

## Table 1

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

Product ID | ASA WSM 1P |

Name | ASAR wide swath standard image |

Description | ASAR product generated from level 0 data collected when the instrument is in wide swath mode. The product includes slant range to ground range corrections and it covers a continuous area along the imaging swath |

Coverage | 400 km×400 km (approximately) for a scene 400 km×4000 km max |

Geometric resolution | ∼150 m×150 m |

Radiometric resolution | Product equivalent number of looks > 11.5 |

Pixel spacing | 75 m×75 m |

Size | 59 Mbytes for a scene; 584 Mbytes for a stripe |

## 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{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$ and $850\times 700\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$, covering an area of $\sim 13,500$ and $\mathrm{13,388}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{km}}^{2}$, respectively.

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.

Step 1 is just computing the parameters (mean and variance) of the Gaussian model by the training data chosen by an expert. The initial distribution of samples takes into account two training windows of $30\times 30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$, which are assigned at pattern regions of classes ${w}_{1}$ and ${w}_{2}$, where ${w}_{1}$ is the sea class and ${w}_{2}$ is the oil spill class. The variance is approached by Eq. (7). For the image of April 26, the initial parameters are ${\mu}_{1}=35$, ${\mu}_{2}=45$, and ${\sigma}^{2}=25$, while for the image of May 2, the initial parameters are ${\mu}_{1}=50$, ${\mu}_{2}=70$, and ${\sigma}^{2}=100$.

Step 2 computes the first $n$ iterations of the simulated annealing method. Each iteration contributes to minimize the energy $U(W|X)$, providing a less noisy label field.

Step 3 updates the Gaussian parameters from the $n$’th MAP label field. The updated parameters for the image of April 26 are ${\mu}_{1}=29.7$, ${\mu}_{2}=58.5$, and ${\sigma}^{2}=207.3$, and that for the image of May 2 are ${\mu}_{1}=47$, ${\mu}_{2}=79.7$, and ${\sigma}^{2}=267.3$. With this recursive method, for $n=1$, the result is just random noise with a uniform distribution. For $n<10$, the segmentation results are still strongly affected by the noise, but for $n\approx 20$, some homogeneous dominions are distinguished, which indicates that the value of the so-called critical temperature of the cooling schedule has been exceeded. By the visual inspection of the obtained field distribution, in a number of $n=30$ iterations, the existence of homogeneous regions can be verified, leaving the remaining small noisy regions to be regularized.

Step 4 aims to obtain fine detections. The neighborhood system applied in the first iterations provides good detection of objects of greater size with an eight-connexity window. The topology described above may oversegment fine structures; by empirical analysis, it was decided that a four-connexity neighborhood system is more adapted to the spatial context of the label model. The number of iterations for the simulated annealing was, thus, fixed to ${n}_{2}$, where ${n}_{2}<n$. As a consequence, vector $\stackrel{\circ}{\mathrm{\Theta}}$ is updated every ${n}_{2}=10$ iterations (during the iteration numbers 30, 40, and 50), see Fig. 1.

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.

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.

## Table 2

Oil slicks area of the two test images.

Date | April 26 | May 2 |

Size | 750×800 | 850×700 |

Total pixels | 600,000 | 595,000 |

Oil slick pixels | 100,264 | 131,084 |

Oil slick area (km2) | 2255.9 | 2949.4 |

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.

## Acknowledgments

The authors are grateful to the European Space Agency for providing ASAR imagery under the project C1F 1069.

## References

## Biography

**Miguel Moctezuma** received his BS and MS degrees in electrical engineering from the National University of Mexico (UNAM) in 1985 and 1988, respectively, and his MSc and PhD degrees in image processing from the ENST, Paris, France, in 1991 and 1995, respectively. He was a visiting fellow with the Italian National Research Council, Bologna, Italy, in 2013. He is a professor with the Telecommunications Department of the Engineering School, UNAM.

**Flavio Parmiggiani** received his MS degree in physics from the University of Milan, Milan, Italy, in 1970. From 1970 to 1982, he worked in the field of biological cybernetics at the Italian National Research Council (CNR). Since then, he has been working in the field of remote sensing and satellite image processing. From 1989 to 2012, he was responsible for the satellite receiving station of the Italian Antarctic Base for real-time operations.