|
1.IntroductionThe 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 filters4,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 paper15 provides a comprehensive analysis of all the issues related to oil spill detection with SAR images. A relevant paper16 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.Methods2.1.Bayesian FrameworkMany 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 be the image to be segmented and the segmentation result, where and is the number of classes (or labels). Bayes’ rule, allows the computation of the posterior probability of class . The conditional probability must be derived and the a priori probability of event 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 are related to one another via a neighborhood system. A neighborhood system for is defined as , where is the set of sites of neighboring pixel . A random region on a lattice with neighborhood system is said to be an MRF if , . A clique for is defined as a subset of sites in . It consists either of a single site , or a pair of neighboring sites , or a triplet of neighboring sites , and so on. The collections of single-site, pair-site, and triple-site cliques will be denoted by , and , respectively. The collection of all cliques for is . The Bayes’ rule can be applied in order to define the equations of the random field model. The joint probability of radiometric information is and is derived from the observed data. The image acquisition process is described by means of the equation . The MAP approach uses the posterior probability in order to obtain an optimization of the segmentation process. The law of the total probability can be ignored and the MAP equation can be approached using the numerator term. In a pixel-based segmentation approach, every pixel is assigned to the 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 . 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 where is the a priori energy function and is the energy function of the observed data model. The Bayesian formulation requires that assumptions are made about prior probabilities . A simple Potts model is used to define the a priori energy function .18 The potential function is given by where is a model parameter which indicates the convergence of the solution, and is the Dirac’s delta function which has a unitary value when the classes of pixels and 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 . with mean values and variance . The Hammersley-Clifford theorem defines the probability density for an MRF under the form of a Gibbs distribution,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 annealing20 given by the Metropolis sampler.21 The probability to belong to a particular class is calculated by means of the posterior energy function , and considering Eqs. (4) and (5), the operation is where 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 , 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 and the label field is obtained. The minimization method is summarized as follows:where is the neighboring pixel system, is the number of iterations, is a proposed initial random solution, and 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 is just a symbolic representation of temperature. In a recursive minimization, has a high value and the probability distribution of the label field is uniform. Labels are randomly assigned with uniform probability, where and is the number of classes (or labels). After a number of iterations, when ;, , is the segmentation result and Eq. (6) attains a minimum value. The configuration of the 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, , and the result of the previous segmentation field, , are required. The update stage is described in the following section and is part of the proposal of this paper. 3.Parameter EstimationThe Gaussian model parameters are unknown but the ergodic property of the stochastic process can be assumed so that the mean values can be estimated from training windows. The variance for every class is approached by where . Inflection points of adjacent probability density functions are crossed in the point . 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, denotes the noisy observed image and is the free-noise version of (or the label field obtained by a segmentation process) with a parameter vector . The problem is to find the labeling field , and in an MAP approach, the biased solution can be described by where are the estimated parameters. To work out a joint solution using incomplete data from Eq. (8) is a difficult task.22 Alternative solutions are the algorithms of expectation maximization and the stochastic expectation maximization.23,24A 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 parameters of the Gaussian model and to perform an iterative adaptive segmentation in order to obtain the field. The pseudo code for solving Eqs. (9) and (10) is as follows: 4.MaterialsAs 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 . 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 1Wide swath medium-resolution image (ASA WSM 1P) product.
5.ResultsFrom 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 and , covering an area of and , respectively. The processing scheme can now start with the adaptive algorithm parameters fixed to and . With regard to the algorithm and to its parameters, some elucidations are needed.
Finally, steps 2 and 3 are repeated until vector 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 , an initial temperature parameter , and second-order cliques. The interaction coefficient 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 is applied: , where 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 2Oil 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.ConclusionsAn 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. AcknowledgmentsThe authors are grateful to the European Space Agency for providing ASAR imagery under the project C1F 1069. ReferencesF. W. Leberl, Radargrammetric Image Processing, Artech House, Norwood, Massachusetts
(1990). Google Scholar
H. ChenD. LiX. Li,
“Mathematical modeling of oil spill on the sea and application of the modeling in Daya Bay,”
J. Hydrodynam. B, 19
(3), 282
–291
(2007). http://dx.doi.org/10.1016/S1001-6058(07)60060-2 JOUHEI 1001-6058 Google Scholar
M. Migliaccioet al.,
“A physically consistent speckle model for marine SLC SAR images,”
IEEE J. Oceanic Eng., 32
(4), 839
–847
(2007). http://dx.doi.org/10.1109/JOE.2007.903985 IJOEDY 0364-9059 Google Scholar
A. BaraldiF. Parmiggiani,
“A refined gamma MAP SAR speckle filter with improved geometrical adaptivity,”
IEEE Trans. Geosci. Remote Sens., 33
(5), 1245
–1257
(1995). http://dx.doi.org/10.1109/36.469489 IGRSD2 0196-2892 Google Scholar
A. BaraldiF. Parmiggiani,
“An alternative form of the Lee filter for speckle suppression in SAR images,”
Graph. Models Image Process., 57
(1), 75
–78
(1995). http://dx.doi.org/10.1006/gmip.1995.1008 CGMPE5 1049-9652 Google Scholar
K. N. Topouzelis,
“Oil spill detection by SAR images: dark formation detection, feature extraction and classification algorithms,”
Sensors, 8
(10), 6642
–6659
(2008). http://dx.doi.org/10.3390/s8106642 SNSRES 0746-9462 Google Scholar
Y. LiJ. Li,
“Oil spill detection from SAR intensity imagery using a marked point process,”
Remote Sens. Environ., 114
(7), 1590
–1601
(2010). http://dx.doi.org/10.1016/j.rse.2010.02.013 RSEEA7 0034-4257 Google Scholar
A. Solberget al.,
“Automatic detection of oil spills in ERS SAR images,”
IEEE Trans. Geosci. Remote Sens., 37
(4), 1916
–1924
(1999). http://dx.doi.org/10.1109/36.774704 IGRSD2 0196-2892 Google Scholar
D. Meraet al.,
“Adaptive thresholding algorithm based on SAR images and wind data to segment oil spills along the northwest coast of the Iberian Peninsula,”
Mar. Pollut. Bull., 64
(10), 2090
–2096
(2012). http://dx.doi.org/10.1016/j.marpolbul.2012.07.018 MPNBAZ 0025-326X Google Scholar
S. SinghaT. J. BellerbyO. Trieschmann,
“Detection and classification of oil spill and look-alike spots from SAR imagery using an artificial neural network,”
in Proc. of IEEE Int. Geoscience and Remote Sensing Symp.,
5630
–5633
(2012). Google Scholar
M. MarghanyA. CracknellM. Hashim,
“Modification of fractal algorithm for oil spill detection from RADARSAT-1 SAR data,”
Int. J. Appl. Earth Obs. Geoinf., 11
(2), 96
–102
(2009). http://dx.doi.org/10.1016/j.jag.2008.09.002 0303-2434 Google Scholar
W. PeiY. Zhu,
“Wavelet transform based edge detection of marine oil spill image,”
in IEEE Int. Conf. on Automation and Logistics,
470
–475
(2008). Google Scholar
C. BrekkeA. Solberg,
“Oil spill detection by satellite remote sensing,”
Remote Sens. Environ., 95
(1), 1
–13
(2005). http://dx.doi.org/10.1016/j.rse.2004.11.015 RSEEA7 0034-4257 Google Scholar
M. MigliaccioM. TranfagliaS. A. Ermakov,
“A physical approach for the observation of oil spills in SAR images,”
IEEE J. Oceanic Eng., 30
(3), 496
–507
(2005). http://dx.doi.org/10.1109/JOE.2005.857518 IJOEDY 0364-9059 Google Scholar
A. Solberg,
“Remote sensing of ocean oil-spill pollution,”
Proc. IEEE, 100
(10), 2931
–2945
(2012). http://dx.doi.org/10.1109/JPROC.2012.2196250 IEEPAD 0018-9219 Google Scholar
M. MigliaccioA. GambardellaM. Tranfaglia,
“SAR polarimetry to observe oil spills,”
IEEE Trans. Geosci. Remote Sens., 45
(2), 506
–511
(2007). http://dx.doi.org/10.1109/TGRS.2006.888097 IGRSD2 0196-2892 Google Scholar
A. Gambardellaet al.,
“One-class classification for oil spill detection,”
Pattern Anal. Appl., 13
(3), 349
–366
(2010). http://dx.doi.org/10.1007/s10044-009-0164-z 1433-7541 Google Scholar
R. KindermannJ. L. Snell, Markov Random Fields and Their Applications, 7th ed.American Mathematical Society, Providence, Rhode Island
(1980). Google Scholar
D. GemanD. Geman,
“Stochastic relaxation Gibbs distributions and the Bayesian restoration of images,”
IEEE Trans. Pattern Anal. Mach. Intell., PAMI 6
(6), 721
–741
(1984). http://dx.doi.org/10.1109/TPAMI.1984.4767596 ITPIDJ 0162-8828 Google Scholar
S. Z. Li, Markov Random Field Modeling in Image Analysis (Advances in Pattern Recognition), 3rd ed.Springer, London
(2009). Google Scholar
N. Metropoliset al.,
“Equations of state calculations by fast computing- machines,”
J. Chem. Phys., 21
(6), 1087
–1092
(1953). http://dx.doi.org/10.1063/1.1699114 JCPSA6 0021-9606 Google Scholar
S. LakshmananH. Derin,
“Simultaneous parameter estimation and segmentation of Gibbs random fields using simulated annealing,”
IEEE Trans. Pattern Anal. Mach. Intell., 11
(8), 799
–813
(1989). http://dx.doi.org/10.1109/34.31443 ITPIDJ 0162-8828 Google Scholar
N. XinB. Yifang,
“An adaptive contextual SEM algorithm for urban land cover mapping using multitemporal high-resolution polarimetric SAR data,”
IEEE J. Sel. Topics Appl. Earth Obs. Remote Sens., 5
(4), 1129
–1139
(2012). http://dx.doi.org/10.1109/JSTARS.2012.2201448 1939-1404 Google Scholar
P. MassonW. Pieczynski,
“SEM algorithm and unsupervised statistical segmentation of satellite images,”
IEEE Trans. Geosci. Remote Sens., 31
(3), 618
–633
(1993). http://dx.doi.org/10.1109/36.225529 IGRSD2 0196-2892 Google Scholar
BiographyMiguel 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. |